Skip to content

Commit

Permalink
Add wave-ice coupling to nuopc/cmeps driver (CICE-Consortium#782)
Browse files Browse the repository at this point in the history 
* merge latest master (NOAA-EMC#4)

* Isotopes for CICE (CICE-Consortium#423)

Co-authored-by: apcraig <[email protected]>
Co-authored-by: David Bailey <[email protected]>
Co-authored-by: Elizabeth Hunke <[email protected]>

* 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. (CICE-Consortium#452) (NOAA-EMC#5)

* Add atmiter_conv to CICE

* Add documentation

* trigger build the docs

Co-authored-by: David A. Bailey <[email protected]>

* update icepack submodule

* Revert "update icepack submodule"

This reverts commit e70d1ab.

* 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 (NOAA-EMC#8)

* cice6 compile (NOAA-EMC#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 (NOAA-EMC#13)

Update to emc/develop fork to current CICE consortium 

Co-authored-by: David A. Bailey <[email protected]>
Co-authored-by: Tony Craig <[email protected]>
Co-authored-by: Elizabeth Hunke <[email protected]>
Co-authored-by: Mariana Vertenstein <[email protected]>
Co-authored-by: apcraig <[email protected]>
Co-authored-by: Philippe Blain <[email protected]>

* Fixcommit (NOAA-EMC#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 (NOAA-EMC#19)

* update icepack

* add ice_dyn_vp module to CICE_InitMod

* update gitmodules, update icepack

* Update CICE to consortium master (NOAA-EMC#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 (NOAA-EMC#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 (NOAA-EMC#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 (NOAA-EMC#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 (NOAA-EMC#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 (NOAA-EMC#40)

* switch to icepack master at consortium

* recreate cap update branch (NOAA-EMC#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 (NOAA-EMC#45)

* Update CICE to Consortium/main (NOAA-EMC#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 (NOAA-EMC#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 <[email protected]>
Co-authored-by: David Bailey <[email protected]>
Co-authored-by: Elizabeth Hunke <[email protected]>
Co-authored-by: Mariana Vertenstein <[email protected]>
Co-authored-by: Minsuk Ji <[email protected]>
Co-authored-by: Tony Craig <[email protected]>
Co-authored-by: Philippe Blain <[email protected]>
  • Loading branch information
@DeniseWorthen @apcraig
@dabail10 @eclare108213 @MinsukJi-NOAA @phil-blain
8 people authored Nov 8, 2022
1 parent aa1e066 commit 251ca48
Show file tree
Hide file tree
Showing 4 changed files with 129 additions and 20 deletions.
6 changes: 6 additions & 0 deletions cicecore/cicedynB/general/ice_state.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
9 changes: 8 additions & 1 deletion cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)'
!----------------------------------------------------

Expand Down Expand Up @@ -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.
Expand Down
22 changes: 11 additions & 11 deletions cicecore/drivers/nuopc/cmeps/ice_comp_nuopc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down
112 changes: 104 additions & 8 deletions cicecore/drivers/nuopc/cmeps/ice_import_export.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)'
!-------------------------------------------------------------------------------
Expand Down Expand Up @@ -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
!-----------------
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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' )
Expand Down Expand Up @@ -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(:)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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(:,:)
Expand All @@ -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__)
Expand All @@ -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
Expand All @@ -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 ???)

Expand Down Expand Up @@ -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
! ------
Expand Down

0 comments on commit 251ca48

Please sign in to comment.