From 251ca48aa8797ac75329ec569cec84e39a094b0a Mon Sep 17 00:00:00 2001 From: Denise Worthen Date: Mon, 7 Nov 2022 19:02:57 -0500 Subject: [PATCH] Add wave-ice coupling to nuopc/cmeps driver (#782) * merge latest master (#4) * Isotopes for CICE (#423) Co-authored-by: apcraig Co-authored-by: David Bailey Co-authored-by: Elizabeth Hunke * updated orbital calculations needed for cesm * fixed problems in updated orbital calculations needed for cesm * update CICE6 to support coupling with UFS * put in changes so that both ufsatm and cesm requirements for potential temperature and density are satisfied * Convergence on ustar for CICE. (#452) (#5) * Add atmiter_conv to CICE * Add documentation * trigger build the docs Co-authored-by: David A. Bailey * update icepack submodule * Revert "update icepack submodule" This reverts commit e70d1abcbeb4351195a2b81c6ce3f623c936426c. * update comp_ice.backend with temporary ice_timers fix * Fix threading problem in init_bgc * Fix additional OMP problems * changes for coldstart running * Move the forapps directory * remove cesmcoupled ifdefs * Fix logging issues for NUOPC * removal of many cpp-ifdefs * fix compile errors * fixes to get cesm working * fixed white space issue * Add restart_coszen namelist option * update icepack submodule * change Orion to orion in backend remove duplicate print lines from ice_transport_driver * add -link_mpi=dbg to debug flags (#8) * cice6 compile (#6) * enable debug build. fix to remove errors * fix an error in comp_ice.backend.libcice * change Orion to orion for machine identification * changes for consistency w/ current emc-cice5 (#13) Update to emc/develop fork to current CICE consortium Co-authored-by: David A. Bailey Co-authored-by: Tony Craig Co-authored-by: Elizabeth Hunke Co-authored-by: Mariana Vertenstein Co-authored-by: apcraig Co-authored-by: Philippe Blain * Fixcommit (#14) Align commit history between emc/develop and cice-consortium/master * Update CICE6 for integration to S2S * add wcoss_dell_p3 compiler macro * update to icepack w/ debug fix * replace SITE with MACHINE_ID * update compile scripts * Support TACC stampede (#19) * update icepack * add ice_dyn_vp module to CICE_InitMod * update gitmodules, update icepack * Update CICE to consortium master (#23) updates include: * deprecate upwind advection (CICE-Consortium#508) * add implicit VP solver (CICE-Consortium#491) * update icepack * switch icepack branches * update to icepack master but set abort flag in ITD routine to false * update icepack * Update CICE to latest Consortium master (#26) update CICE and Icepack * changes the criteria for aborting ice for thermo-conservation errors * updates the time manager * fixes two bugs in ice_therm_mushy * updates Icepack to Consortium master w/ flip of abort flag for troublesome IC cases * add cice changes for zlvs (#29) * update icepack and pointer * update icepack and revert gitmodules * Fix history features - Fix bug in history time axis when sec_init is not zero. - Fix issue with time_beg and time_end uninitialized values. - Add support for averaging with histfreq='1' by allowing histfreq_n to be any value in that case. Extend and clean up construct_filename for history files. More could be done, but wanted to preserve backwards compatibility. - Add new calendar_sec2hms to converts daily seconds to hh:mm:ss. Update the calchk calendar unit tester to check this method - Remove abort test in bcstchk, this was just causing problems in regression testing - Remove known problems documentation about problems writing when istep=1. This issue does not exist anymore with the updated time manager. - Add new tests with hist_avg = false. Add set_nml.histinst. * revert set_nml.histall * fix implementation error * update model log output in ice_init * Fix QC issues - Add netcdf ststus checks and aborts in ice_read_write.F90 - Check for end of file when reading records in ice_read_write.F90 for ice_read_nc methods - Update set_nml.qc to better specify the test, turn off leap years since we're cycling 2005 data - Add check in c ice.t-test.py to make sure there is at least 1825 files, 5 years of data - Add QC run to base_suite.ts to verify qc runs to completion and possibility to use those results directly for QC validation - Clean up error messages and some indentation in ice_read_write.F90 * Update testing - Add prod suite including 10 year gx1prod and qc test - Update unit test compare scripts * update documentation * reset calchk to 100000 years * update evp1d test * update icepack * update icepack * add memory profiling (#36) * add profile_memory calls to CICE cap * update icepack * fix rhoa when lowest_temp is 0.0 * provide default value for rhoa when imported temp_height_lowest (Tair) is 0.0 * resolves seg fault when frac_grid=false and do_ca=true * update icepack submodule * Update CICE for latest Consortium master (#38) * Implement advanced snow physics in icepack and CICE * Fix time-stamping of CICE history files * Fix CICE history file precision * Use CICE-Consortium/Icepack master (#40) * switch to icepack master at consortium * recreate cap update branch (#42) * add debug_model feature * add required variables and calls for tr_snow * remove 2 extraneous lines * remove two log print lines that were removed prior to merge of driver updates to consortium * duplicate gitmodule style for icepack * Update CICE to latest Consortium/main (#45) * Update CICE to Consortium/main (#48) Update OpenMP directives as needed including validation via new omp_suite. Fixed OpenMP in dynamics. Refactored eap puny/pi lookups to improve scalar performance Update Tsfc implementation to make sure land blocks don't set Tsfc to freezing temp Update for sea bed stress calculations * fix comment, fix env for orion and hera * replace save_init with step_prep in CICE_RunMod * fixes for cgrid repro * remove added haloupdates * baselines pass with these extra halo updates removed * change F->S for ocean velocities and tilts * fix debug failure when grid_ice=C * compiling in debug mode using -init=snan,arrays requires initialization of variables * respond to review comments * remove inserted whitespace for uvelE,N and vvelE,N * Add wave-cice coupling; update to Consortium main (#51) * add wave-ice fields * initialize aicen_init, which turns up as NaN in calc of floediam export * add call to icepack_init_wave to initialize wavefreq and dwavefreq * update to latest consortium main (PR 752) * add initializationsin ice_state * initialize vsnon/vsnon_init and vicen/vicen_init Co-authored-by: apcraig Co-authored-by: David Bailey Co-authored-by: Elizabeth Hunke Co-authored-by: Mariana Vertenstein Co-authored-by: Minsuk Ji <57227195+MinsukJi-NOAA@users.noreply.github.com> Co-authored-by: Tony Craig Co-authored-by: Philippe Blain --- cicecore/cicedynB/general/ice_state.F90 | 6 + cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 | 9 +- .../drivers/nuopc/cmeps/ice_comp_nuopc.F90 | 22 ++-- .../drivers/nuopc/cmeps/ice_import_export.F90 | 112 ++++++++++++++++-- 4 files changed, 129 insertions(+), 20 deletions(-) diff --git a/cicecore/cicedynB/general/ice_state.F90 b/cicecore/cicedynB/general/ice_state.F90 index 7b20879bc..862f0a8bc 100644 --- a/cicecore/cicedynB/general/ice_state.F90 +++ b/cicecore/cicedynB/general/ice_state.F90 @@ -188,6 +188,12 @@ subroutine alloc_state n_trcr_strata = 0 nt_strata = 0 trcr_base = c0 + aicen = c0 + aicen_init = c0 + vicen = c0 + vicen_init = c0 + vsnon = c0 + vsnon_init = c0 end subroutine alloc_state diff --git a/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 b/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 index f9b5116d0..3d5e5cc2a 100644 --- a/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 +++ b/cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90 @@ -78,7 +78,7 @@ subroutine cice_init2() use ice_calendar , only: dt, dt_dyn, istep, istep1, write_ic, init_calendar, calendar use ice_communicate , only: my_task, master_task use ice_diagnostics , only: init_diags - use ice_domain_size , only: ncat, nfsd + use ice_domain_size , only: ncat, nfsd, nfreq use ice_dyn_eap , only: init_eap, alloc_dyn_eap use ice_dyn_shared , only: kdyn, init_dyn use ice_dyn_vp , only: init_vp @@ -94,10 +94,12 @@ subroutine cice_init2() use ice_restoring , only: ice_HaloRestore_init use ice_timers , only: timer_total, init_ice_timers, ice_timer_start use ice_transport_driver , only: init_transport + use ice_arrays_column , only: wavefreq, dwavefreq logical(kind=log_kind) :: tr_aero, tr_zaero, skl_bgc, z_tracers logical(kind=log_kind) :: tr_iso, tr_fsd, wave_spec, tr_snow character(len=char_len) :: snw_aging_table + real(kind=dbl_kind), dimension(25) :: wave_spectrum_profile ! hardwire for now character(len=*), parameter :: subname = '(cice_init2)' !---------------------------------------------------- @@ -177,6 +179,11 @@ subroutine cice_init2() endif endif + if (wave_spec) then + call icepack_init_wave(nfreq=nfreq, & + wave_spectrum_profile=wave_spectrum_profile, wavefreq=wavefreq, dwavefreq=dwavefreq) + end if + ! Initialize shortwave components using swdn from previous timestep ! if restarting. These components will be scaled to current forcing ! in prep_radiation. diff --git a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 index 8920ea386..182308973 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90 @@ -25,6 +25,7 @@ module ice_comp_nuopc use ice_calendar , only : force_restart_now, write_ic use ice_calendar , only : idate, mday, mmonth, myear, year_init use ice_calendar , only : msec, dt, calendar, calendar_type, nextsw_cday, istep + use ice_calendar , only : ice_calendar_noleap, ice_calendar_gregorian use ice_kinds_mod , only : dbl_kind, int_kind, char_len, char_len_long use ice_fileunits , only : nu_diag, nu_diag_set, inst_index, inst_name use ice_fileunits , only : inst_suffix, release_all_fileunits, flush_fileunit @@ -80,9 +81,6 @@ module ice_comp_nuopc character(len=*) , parameter :: orb_variable_year = 'variable_year' character(len=*) , parameter :: orb_fixed_parameters = 'fixed_parameters' - character(len=*),parameter :: shr_cal_noleap = 'NO_LEAP' - character(len=*),parameter :: shr_cal_gregorian = 'GREGORIAN' - type(ESMF_Mesh) :: ice_mesh integer :: nthrds ! Number of threads to use in this component @@ -216,7 +214,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) type(ESMF_Time) :: stopTime ! Stop time type(ESMF_Time) :: refTime ! Ref time type(ESMF_TimeInterval) :: timeStep ! Model timestep - type(ESMF_Calendar) :: esmf_calendar ! esmf calendar type(ESMF_CalKind_Flag) :: esmf_caltype ! esmf calendar type integer :: start_ymd ! Start date (YYYYMMDD) integer :: start_tod ! start time of day (s) @@ -339,7 +336,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call get_component_instance(gcomp, inst_suffix, inst_index, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - inst_name = "ICE"//trim(inst_suffix) +! inst_name = "ICE"//trim(inst_suffix) + inst_name = "ICE" !---------------------------------------------------------------------------- ! start cice timers @@ -470,9 +468,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call ESMF_TimeGet( currTime, calkindflag=esmf_caltype, rc=rc ) if (ChkErr(rc,__LINE__,u_FILE_u)) return if (esmf_caltype == ESMF_CALKIND_NOLEAP) then - calendar_type = shr_cal_noleap + calendar_type = ice_calendar_noleap else if (esmf_caltype == ESMF_CALKIND_GREGORIAN) then - calendar_type = shr_cal_gregorian + calendar_type = ice_calendar_gregorian else call abort_ice( subname//'ERROR:: bad calendar for ESMF' ) end if @@ -581,9 +579,11 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) end if call icepack_query_parameters( tfrz_option_out=tfrz_option) if (tfrz_option_driver /= tfrz_option) then - write(errmsg,'(a)') trim(subname)//'error: tfrz_option from driver '//trim(tfrz_option_driver)//& - ' must be the same as tfrz_option from cice namelist '//trim(tfrz_option) - call abort_ice(trim(errmsg)) + write(errmsg,'(a)') trim(subname)//'WARNING: tfrz_option from driver '//trim(tfrz_option_driver)//& + ' is overwriting tfrz_option from cice namelist '//trim(tfrz_option) + write(nu_diag,*) trim(errmsg) + call icepack_warnings_flush(nu_diag) + call icepack_init_parameters(tfrz_option_in=tfrz_option_driver) endif ! Flux convergence tolerance - always use the driver attribute value @@ -594,7 +594,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) read(cvalue,*) atmiter_conv_driver call icepack_query_parameters( atmiter_conv_out=atmiter_conv) if (atmiter_conv_driver /= atmiter_conv) then - write(errmsg,'(a,d13.5,a,d13.5)') trim(subname)//'warning: atmiter_ from driver ',& + write(errmsg,'(a,d13.5,a,d13.5)') trim(subname)//'WARNING: atmiter_ from driver ',& atmiter_conv_driver,' is overwritting atmiter_conv from cice namelist ',atmiter_conv write(nu_diag,*) trim(errmsg) call icepack_warnings_flush(nu_diag) diff --git a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 index 5a6ce7572..d95a4d9b2 100644 --- a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 +++ b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 @@ -3,12 +3,13 @@ module ice_import_export use ESMF use NUOPC use NUOPC_Model - use ice_kinds_mod , only : int_kind, dbl_kind, char_len, log_kind + use ice_kinds_mod , only : int_kind, dbl_kind, char_len, char_len_long, log_kind use ice_constants , only : c0, c1, spval_dbl, radius use ice_constants , only : field_loc_center, field_type_scalar, field_type_vector use ice_blocks , only : block, get_block, nx_block, ny_block use ice_domain , only : nblocks, blocks_ice, halo_info, distrb_info use ice_domain_size , only : nx_global, ny_global, block_size_x, block_size_y, max_blocks, ncat + use ice_domain_size , only : nfreq, nfsd use ice_exit , only : abort_ice use ice_flux , only : strairxT, strairyT, strocnxT_iavg, strocnyT_iavg use ice_flux , only : alvdr, alidr, alvdf, alidf, Tref, Qref, Uref @@ -23,9 +24,10 @@ module ice_import_export use ice_flux , only : fsnow, uocn, vocn, sst, ss_tltx, ss_tlty, frzmlt use ice_flux , only : send_i2x_per_cat use ice_flux , only : sss, Tf, wind, fsw - use ice_state , only : vice, vsno, aice, aicen_init, trcr + use ice_arrays_column , only : floe_rad_c, wave_spectrum + use ice_state , only : vice, vsno, aice, aicen_init, trcr, trcrn use ice_grid , only : tlon, tlat, tarea, tmask, anglet, hm - use ice_grid , only : grid_type, grid_average_X2Y + use ice_grid , only : grid_type use ice_mesh_mod , only : ocn_gridcell_frac use ice_boundary , only : ice_HaloUpdate use ice_fileunits , only : nu_diag, flush_fileunit @@ -34,8 +36,10 @@ module ice_import_export use ice_shr_methods , only : chkerr, state_reset use icepack_intfc , only : icepack_warnings_flush, icepack_warnings_aborted use icepack_intfc , only : icepack_query_parameters, icepack_query_tracer_flags + use icepack_intfc , only : icepack_query_tracer_indices use icepack_intfc , only : icepack_liquidus_temperature use icepack_intfc , only : icepack_sea_freezing_temperature + use icepack_parameters , only : puny, c2 use cice_wrapper_mod , only : t_startf, t_stopf, t_barrierf #ifdef CESMCOUPLED use shr_frz_mod , only : shr_frz_freezetemp @@ -112,6 +116,7 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam character(char_len) :: stdname character(char_len) :: cvalue logical :: flds_wiso ! use case + logical :: flds_wave ! use case logical :: isPresent, isSet character(len=*), parameter :: subname='(ice_import_export:ice_advertise_fields)' !------------------------------------------------------------------------------- @@ -148,6 +153,17 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam write(nu_diag,*)'flds_wiso = ',flds_wiso end if + flds_wave = .false. + call NUOPC_CompAttributeGet(gcomp, name='wav_coupling_to_cice', value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + read(cvalue,*) flds_wave + end if + if (my_task == master_task) then + write(nu_diag,*)'flds_wave = ',flds_wave + end if + !----------------- ! advertise import fields !----------------- @@ -192,6 +208,14 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam ! from atm - dry dust deposition fluxes (4 sizes) call fldlist_add(fldsToIce_num, fldsToIce, 'Faxa_dstdry', ungridded_lbound=1, ungridded_ubound=4) + ! the following are advertised but might not be connected if they are not advertised in the + ! in the cmeps esmFldsExchange_xxx_mod.F90 that is model specific + ! from wave + if (flds_wave) then + call fldlist_add(fldsToIce_num, fldsToIce, 'Sw_elevation_spectrum', ungridded_lbound=1, & + ungridded_ubound=25) + end if + do n = 1,fldsToIce_num call NUOPC_Advertise(importState, standardName=fldsToIce(n)%stdname, & TransferOfferGeomObject='will provide', rc=rc) @@ -225,6 +249,10 @@ subroutine ice_advertise_fields(gcomp, importState, exportState, flds_scalar_nam call fldlist_add(fldsFrIce_num, fldsFrIce, 'ice_fraction_n', & ungridded_lbound=1, ungridded_ubound=ncat) end if + if (flds_wave) then + call fldlist_add(fldsFrIce_num, fldsFrIce, 'Si_thick' ) + call fldlist_add(fldsFrIce_num, fldsFrIce, 'Si_floediam' ) + end if ! ice/atm fluxes computed by ice call fldlist_add(fldsFrIce_num, fldsFrIce, 'stress_on_air_ice_zonal' ) @@ -292,7 +320,7 @@ subroutine ice_realize_fields(gcomp, mesh, flds_scalar_name, flds_scalar_num, rc type(ESMF_State) :: exportState type(ESMF_Field) :: lfield integer :: numOwnedElements - integer :: i, j, iblk, n + integer :: i, j, iblk, n, k integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain type(block) :: this_block ! block information for current block real(dbl_kind), allocatable :: mesh_areas(:) @@ -403,11 +431,10 @@ subroutine ice_import( importState, rc ) ! local variables integer,parameter :: nflds=16 integer,parameter :: nfldv=6 - integer :: i, j, iblk, n + integer :: i, j, iblk, n, k integer :: ilo, ihi, jlo, jhi !beginning and end of physical domain type(block) :: this_block ! block information for current block real (kind=dbl_kind),allocatable :: aflds(:,:,:,:) - real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: work real (kind=dbl_kind) :: workx, worky real (kind=dbl_kind) :: MIN_RAIN_TEMP, MAX_SNOW_TEMP real (kind=dbl_kind) :: Tffresh @@ -559,6 +586,29 @@ subroutine ice_import( importState, rc ) end do !$OMP END PARALLEL DO + ! import wave elevation spectrum from wave (frequencies 1-25, assume that nfreq is 25) + if (State_FldChk(importState, 'Sw_elevation_spectrum')) then + if (nfreq /= 25) then + call abort_ice(trim(subname)//": ERROR nfreq not equal to 25 ") + end if + call state_getfldptr(importState, 'Sw_elevation_spectrum', fldptr=dataPtr2d, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + do k = 1,nfreq + n = 0 + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo; ihi = this_block%ihi + jlo = this_block%jlo; jhi = this_block%jhi + do j = jlo, jhi + do i = ilo, ihi + n = n+1 + wave_spectrum(i,j,k,iblk) = dataPtr2d(k,n) + end do + end do + end do + end do + end if + if ( State_fldChk(importState, 'Sa_ptem') .and. State_fldchk(importState,'air_density_height_lowest')) then !$OMP PARALLEL DO PRIVATE(iblk,i,j) do iblk = 1, nblocks @@ -845,7 +895,7 @@ subroutine ice_export( exportState, rc ) ! local variables type(block) :: this_block ! block information for current block - integer :: i, j, iblk, n ! incides + integer :: i, j, iblk, n, k ! indices integer :: n2 ! thickness category index integer :: ilo, ihi, jlo, jhi ! beginning and end of physical domain real (kind=dbl_kind) :: workx, worky ! tmps for converting grid @@ -859,7 +909,11 @@ subroutine ice_export( exportState, rc ) real (kind=dbl_kind) :: tauxo (nx_block,ny_block,max_blocks) ! ice/ocean stress real (kind=dbl_kind) :: tauyo (nx_block,ny_block,max_blocks) ! ice/ocean stress real (kind=dbl_kind) :: ailohi(nx_block,ny_block,max_blocks) ! fractional ice area + real (kind=dbl_kind) :: floediam(nx_block,ny_block,max_blocks) + real (kind=dbl_kind) :: floethick(nx_block,ny_block,max_blocks) ! ice thickness real (kind=dbl_kind) :: Tffresh + logical (kind=log_kind) :: tr_fsd + integer (kind=int_kind) :: nt_fsd real (kind=dbl_kind), allocatable :: tempfld(:,:,:) real (kind=dbl_kind), pointer :: dataptr_ifrac_n(:,:) real (kind=dbl_kind), pointer :: dataptr_swpen_n(:,:) @@ -877,6 +931,9 @@ subroutine ice_export( exportState, rc ) ! tr_FY_out=tr_FY, tr_pond_out=tr_pond, tr_lvl_out=tr_lvl, & ! tr_zaero_out=tr_zaero, tr_bgc_Nit_out=tr_bgc_Nit) + call icepack_query_tracer_indices(nt_fsd_out=nt_fsd) + call icepack_query_tracer_flags(tr_fsd_out=tr_fsd) + call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=u_FILE_u, line=__LINE__) @@ -890,8 +947,10 @@ subroutine ice_export( exportState, rc ) tauya(:,:,:) = c0 tauxo(:,:,:) = c0 tauyo(:,:,:) = c0 + floediam(:,:,:) = c0 + floethick(:,:,:) = c0 - !$OMP PARALLEL DO PRIVATE(iblk,i,j,workx,worky, this_block, ilo, ihi, jlo, jhi) + !$OMP PARALLEL DO PRIVATE(iblk,i,j,k,workx,worky, this_block, ilo, ihi, jlo, jhi) do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) ilo = this_block%ilo @@ -904,6 +963,27 @@ subroutine ice_export( exportState, rc ) ! ice fraction ailohi(i,j,iblk) = min(aice(i,j,iblk), c1) + if (tr_fsd) then + ! floe thickness (m) + if (aice(i,j,iblk) > puny) then + floethick(i,j,iblk) = vice(i,j,iblk) / aice(i,j,iblk) + else + floethick(i,j,iblk) = c0 + end if + + ! floe diameter (m) + workx = c0 + worky = c0 + do n = 1, ncat + do k = 1, nfsd + workx = workx + floe_rad_c(k) * aicen_init(i,j,n,iblk) * trcrn(i,j,nt_fsd+k-1,n,iblk) + worky = worky + aicen_init(i,j,n,iblk) * trcrn(i,j,nt_fsd+k-1,n,iblk) + end do + end do + if (worky > c0) workx = c2*workx / worky + floediam(i,j,iblk) = MAX(c2*floe_rad_c(1),workx) + endif + ! surface temperature Tsrf(i,j,iblk) = Tffresh + trcr(i,j,1,iblk) !Kelvin (original ???) @@ -1054,6 +1134,22 @@ subroutine ice_export( exportState, rc ) call state_setexport(exportState, 'Si_snowh' , input=tempfld , lmask=tmask, ifrac=ailohi, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! ------ + ! optional floe diameter and ice thickness to wave + ! ------ + + ! Sea ice thickness (m) + if (State_FldChk(exportState, 'Si_thick')) then + call state_setexport(exportState, 'Si_thick' , input=floethick , lmask=tmask, ifrac=ailohi, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + + ! Sea ice floe diameter (m) + if (State_FldChk(exportState, 'Si_floediam')) then + call state_setexport(exportState, 'Si_floediam' , input=floediam , lmask=tmask, ifrac=ailohi, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + ! ------ ! ice/atm fluxes computed by ice ! ------