Update for the FSD changes in icepack. (#975)

Update fsd implementation and interfaces.  

This changes answers for all results due to changes in Icepack. Initial ice and snow volume values at the beginning of the lateral melt routine are now used for updating the snow and ice enthalpy as well other tracers for lateral melting.

Change calls to icepack_step_therm1 and icepack_step_therm2 for the new FSD rearrangement.  Update calls to step_dyn_wave, step_therm2, step_therm1 to be consistent with new Icepack interfaces

Update fsd arguments floe_rad_l, floe_rad_c, floe_binwidth, c_fsd_range to Icepack.  Add floe_rad_l, floe_rad_c, floe_binwidth, and c_fsd_range out arguments to icepack_init_fsd_bounds calls.

Update rside/fside arguments to Icepack

Update unittest and coupler layers as needed

Update driver code including remove zsal from mct/cesm.

Update Icepack to #286630f

---------

Co-authored-by: apcraig <[email protected]>

cicecore/cicedyn/general/ice_flux.F90

@@ -368,15 +368,14 @@ module ice_flux
       real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
          uatmT   , & ! uatm on T grid (m/s)
          vatmT   , & ! vatm on T grid (m/s)
-         rside   , & ! fraction of ice that melts laterally
-         fside   , & ! lateral heat flux (W/m^2)
          wlat    , & ! lateral heat rate (m/s)
          fsw     , & ! incoming shortwave radiation (W/m^2)
          coszen  , & ! cosine solar zenith angle, < 0 for sun below horizon
          rdg_conv, & ! convergence term for ridging (1/s)
          rdg_shear   ! shear term for ridging (1/s)
 
       real (kind=dbl_kind), dimension(:,:,:,:), allocatable, public :: &
+         rsiden    ,&   ! fraction of ice that melts laterally
          salinz    ,&   ! initial salinity  profile (ppt)
          Tmltz          ! initial melting temperature (^oC)
 
@@ -546,13 +545,12 @@ subroutine alloc_flux
          fsalt_da   (nx_block,ny_block,max_blocks), & ! salt flux to ocean due to data assimilation(kg/m^2/s)
          uatmT      (nx_block,ny_block,max_blocks), & ! uatm on T grid
          vatmT      (nx_block,ny_block,max_blocks), & ! vatm on T grid
-         rside      (nx_block,ny_block,max_blocks), & ! fraction of ice that melts laterally
-         fside      (nx_block,ny_block,max_blocks), & ! lateral melt flux (W/m^2)
          wlat       (nx_block,ny_block,max_blocks), & ! lateral melt rate (m/s)
          fsw        (nx_block,ny_block,max_blocks), & ! incoming shortwave radiation (W/m^2)
          coszen     (nx_block,ny_block,max_blocks), & ! cosine solar zenith angle, < 0 for sun below horizon
          rdg_conv   (nx_block,ny_block,max_blocks), & ! convergence term for ridging (1/s)
          rdg_shear  (nx_block,ny_block,max_blocks), & ! shear term for ridging (1/s)
+         rsiden     (nx_block,ny_block,ncat,max_blocks), & ! fraction of ice that melts laterally
          dardg1ndt  (nx_block,ny_block,ncat,max_blocks), & ! rate of area loss by ridging ice (1/s)
          dardg2ndt  (nx_block,ny_block,ncat,max_blocks), & ! rate of area gain by new ridges (1/s)
          dvirdgndt  (nx_block,ny_block,ncat,max_blocks), & ! rate of ice volume ridged (m/s)

cicecore/cicedyn/general/ice_forcing.F90

@@ -5473,8 +5473,6 @@ subroutine wave_spec_data
       logical (kind=log_kind) :: wave_spec
       character(len=*), parameter :: subname = '(wave_spec_data)'
 
-
-
       debug_n_d = .false.  !usually false
 
       call icepack_query_parameters(secday_out=secday)

cicecore/cicedyn/general/ice_step_mod.F90

@@ -224,9 +224,9 @@ subroutine step_therm1 (dt, iblk)
           fswsfcn, fswintn, Sswabsn, Iswabsn, meltsliqn, meltsliq, &
           fswthrun, fswthrun_vdr, fswthrun_vdf, fswthrun_idr, fswthrun_idf
       use ice_calendar, only: yday
-      use ice_domain_size, only: ncat, nilyr, nslyr, n_iso, n_aero
-      use ice_flux, only: frzmlt, sst, Tf, strocnxT_iavg, strocnyT_iavg, rside, fbot, Tbot, Tsnice, &
-          meltsn, melttn, meltbn, congeln, snoicen, uatmT, vatmT, fside, wlat, &
+      use ice_domain_size, only: ncat, nilyr, nslyr, n_iso, n_aero, nfsd
+      use ice_flux, only: frzmlt, sst, Tf, strocnxT_iavg, strocnyT_iavg, rsiden, fbot, Tbot, Tsnice, &
+          meltsn, melttn, meltbn, congeln, snoicen, uatmT, vatmT, wlat, &
           wind, rhoa, potT, Qa, zlvl, zlvs, strax, stray, flatn, fsensn, fsurfn, fcondtopn, &
           flw, fsnow, fpond, sss, mlt_onset, frz_onset, fcondbotn, fcondbot, fsloss, &
           frain, Tair, strairxT, strairyT, fsurf, fcondtop, fsens, &
@@ -269,7 +269,7 @@ subroutine step_therm1 (dt, iblk)
 
       integer (kind=int_kind) :: &
          ntrcr, nt_apnd, nt_hpnd, nt_ipnd, nt_alvl, nt_vlvl, nt_Tsfc, &
-         nt_iage, nt_FY, nt_qice, nt_sice, nt_aero, nt_qsno, &
+         nt_iage, nt_FY, nt_qice, nt_sice, nt_aero, nt_qsno, nt_fsd, &
          nt_isosno, nt_isoice, nt_rsnw, nt_smice, nt_smliq
 
       logical (kind=log_kind) :: &
@@ -304,7 +304,7 @@ subroutine step_therm1 (dt, iblk)
       call icepack_query_tracer_indices( &
          nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, &
          nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, nt_Tsfc_out=nt_Tsfc, &
-         nt_iage_out=nt_iage, nt_FY_out=nt_FY, &
+         nt_iage_out=nt_iage, nt_FY_out=nt_FY, nt_fsd_out=nt_fsd, &
          nt_qice_out=nt_qice, nt_sice_out=nt_sice, &
          nt_aero_out=nt_aero, nt_qsno_out=nt_qsno, &
          nt_rsnw_out=nt_rsnw, nt_smice_out=nt_smice, nt_smliq_out=nt_smliq, &
@@ -412,6 +412,7 @@ subroutine step_therm1 (dt, iblk)
                       ipnd         = trcrn       (i,j,nt_ipnd,:,iblk),                   &
                       iage         = trcrn       (i,j,nt_iage,:,iblk),                   &
                       FY           = trcrn       (i,j,nt_FY  ,:,iblk),                   &
+                      afsdn        = trcrn       (i,j,nt_fsd:nt_fsd+nfsd-1,:,iblk),      &
                       rsnwn        = rsnwn       (:,:),        &
                       smicen       = smicen      (:,:),        &
                       smliqn       = smliqn      (:,:),        &
@@ -464,8 +465,7 @@ subroutine step_therm1 (dt, iblk)
                       Tbot         = Tbot        (i,j,  iblk), &
                       Tsnice       = Tsnice      (i,j,  iblk), &
                       frzmlt       = frzmlt      (i,j,  iblk), &
-                      rside        = rside       (i,j,  iblk), &
-                      fside        = fside       (i,j,  iblk), &
+                      rsiden       = rsiden      (i,j,:,iblk), &
                       wlat         = wlat        (i,j,  iblk), &
                       fsnow        = fsnow       (i,j,  iblk), &
                       frain        = frain       (i,j,  iblk), &
@@ -611,12 +611,12 @@ subroutine step_therm2 (dt, iblk)
 
       use ice_arrays_column, only: hin_max, ocean_bio, wave_sig_ht, &
           wave_spectrum, wavefreq, dwavefreq, &
-          first_ice, bgrid, cgrid, igrid, floe_rad_c, floe_binwidth, &
+          first_ice, bgrid, cgrid, igrid, &
           d_afsd_latg, d_afsd_newi, d_afsd_latm, d_afsd_weld
       use ice_calendar, only: yday
       use ice_domain_size, only: ncat, nilyr, nslyr, nblyr, nfsd
       use ice_flux, only: fresh, frain, fpond, frzmlt, frazil, frz_onset, &
-          fsalt, Tf, sss, salinz, fhocn, rside, fside, wlat, &
+          fsalt, Tf, sss, salinz, fhocn, rsiden, wlat, &
           meltl, frazil_diag
       use ice_flux_bgc, only: flux_bio, faero_ocn, &
           fiso_ocn, HDO_ocn, H2_16O_ocn, H2_18O_ocn
@@ -695,9 +695,8 @@ subroutine step_therm2 (dt, iblk)
                       Tf         = Tf        (i,j,  iblk), &
                       sss        = sss       (i,j,  iblk), &
                       salinz     = salinz    (i,j,:,iblk), &
-                      rside      = rside     (i,j,  iblk), &
+                      rsiden     = rsiden    (i,j,:,iblk), &
                       meltl      = meltl     (i,j,  iblk), &
-                      fside      = fside     (i,j,  iblk), &
                       wlat       = wlat      (i,j,  iblk), &
                       frzmlt     = frzmlt    (i,j,  iblk), &
                       frazil     = frazil    (i,j,  iblk), &
@@ -724,9 +723,7 @@ subroutine step_therm2 (dt, iblk)
                       d_afsd_latg= d_afsd_latg(i,j,:,iblk),&
                       d_afsd_newi= d_afsd_newi(i,j,:,iblk),&
                       d_afsd_latm= d_afsd_latm(i,j,:,iblk),&
-                      d_afsd_weld= d_afsd_weld(i,j,:,iblk),&
-                      floe_rad_c = floe_rad_c(:),          &
-                      floe_binwidth = floe_binwidth(:))
+                      d_afsd_weld= d_afsd_weld(i,j,:,iblk))
          endif ! tmask
 
       enddo                     ! i
@@ -863,7 +860,7 @@ end subroutine update_state
       subroutine step_dyn_wave (dt)
 
       use ice_arrays_column, only: wave_spectrum, &
-          d_afsd_wave, floe_rad_l, floe_rad_c, wavefreq, dwavefreq
+          d_afsd_wave, wavefreq, dwavefreq
       use ice_domain_size, only: ncat, nfsd, nfreq
       use ice_state, only: trcrn, aicen, aice, vice
       use ice_timers, only: ice_timer_start, ice_timer_stop, timer_column, &
@@ -911,8 +908,6 @@ subroutine step_dyn_wave (dt)
                                            aice        = aice           (i,j,    iblk), &
                                            vice        = vice           (i,j,    iblk), &
                                            aicen       = aicen          (i,j,:,  iblk), &
-                                           floe_rad_l  = floe_rad_l     (:),            &
-                                           floe_rad_c  = floe_rad_c     (:),            &
                                            wave_spectrum = wave_spectrum(i,j,:,  iblk), &
                                            wavefreq    = wavefreq       (:),            &
                                            dwavefreq   = dwavefreq      (:),            &

cicecore/drivers/mct/cesm1/CICE_InitMod.F90

@@ -158,10 +158,10 @@ subroutine cice_init(mpicom_ice)
           file=__FILE__,line= __LINE__)
 
       if (tr_fsd) call icepack_init_fsd_bounds ( &
-         floe_rad_l = floe_rad_l,    &  ! fsd size lower bound in m (radius)
-         floe_rad_c = floe_rad_c,    &  ! fsd size bin centre in m (radius)
-         floe_binwidth = floe_binwidth, &  ! fsd size bin width in m (radius)
-         c_fsd_range = c_fsd_range,   &  ! string for history output
+         floe_rad_l_out = floe_rad_l,    &  ! fsd size lower bound in m (radius)
+         floe_rad_c_out = floe_rad_c,    &  ! fsd size bin centre in m (radius)
+         floe_binwidth_out = floe_binwidth, &  ! fsd size bin width in m (radius)
+         c_fsd_range_out = c_fsd_range,   &  ! string for history output
          write_diags=(my_task == master_task))  ! write diag on master only
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
@@ -199,6 +199,8 @@ subroutine cice_init(mpicom_ice)
       if (trim(runtype) == 'continue' .or. restart) &
          call init_shortwave    ! initialize radiative transfer
 
+      if (write_ic) call accum_hist(dt) ! write initial conditions
+
 !      call advance_timestep()
 
    !--------------------------------------------------------------------
@@ -238,8 +240,6 @@ subroutine cice_init(mpicom_ice)
       call init_flux_atm        ! initialize atmosphere fluxes sent to coupler
       call init_flux_ocn        ! initialize ocean fluxes sent to coupler
 
-      if (write_ic) call accum_hist(dt) ! write initial conditions
-
       call dealloc_grid         ! deallocate temporary grid arrays
       if (my_task == master_task) then
          call ice_memusage_print(nu_diag,subname//':end')
@@ -274,7 +274,7 @@ subroutine init_restart
           restart_iso, read_restart_iso, &
           restart_aero, read_restart_aero, &
           restart_hbrine, read_restart_hbrine, &
-          restart_zsal, restart_bgc
+          restart_bgc
       use ice_restart_driver, only: restartfile
       use ice_restart_shared, only: runtype, restart
       use ice_state ! almost everything
@@ -285,7 +285,7 @@ subroutine init_restart
       logical(kind=log_kind) :: &
           tr_iage, tr_FY, tr_lvl, tr_pond_lvl, &
           tr_pond_topo, tr_snow, tr_fsd, tr_iso, tr_aero, tr_brine, &
-          skl_bgc, z_tracers, solve_zsal
+          skl_bgc, z_tracers
       integer(kind=int_kind) :: &
           ntrcr
       integer(kind=int_kind) :: &
@@ -301,7 +301,7 @@ subroutine init_restart
           file=__FILE__, line=__LINE__)
 
       call icepack_query_parameters(skl_bgc_out=skl_bgc, &
-           z_tracers_out=z_tracers, solve_zsal_out=solve_zsal)
+           z_tracers_out=z_tracers)
       call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, &
            tr_lvl_out=tr_lvl, tr_pond_lvl_out=tr_pond_lvl, &
            tr_pond_topo_out=tr_pond_topo, tr_aero_out=tr_aero, tr_brine_out=tr_brine, &
@@ -447,8 +447,6 @@ subroutine init_restart
       if (trim(runtype) == 'continue') then
          if (tr_brine) &
              restart_hbrine = .true.
-         if (solve_zsal) &
-             restart_zsal = .true.
          if (skl_bgc .or. z_tracers) &
              restart_bgc = .true.
       endif
@@ -458,7 +456,7 @@ subroutine init_restart
          if (tr_brine .and. restart_hbrine) call read_restart_hbrine
       endif
 
-      if (solve_zsal .or. skl_bgc .or. z_tracers) then ! biogeochemistry
+      if (skl_bgc .or. z_tracers) then ! biogeochemistry
          if (tr_fsd) then
             write (nu_diag,*) 'FSD implementation incomplete for use with BGC'
             call icepack_warnings_flush(nu_diag)

cicecore/drivers/nuopc/cmeps/CICE_InitMod.F90

@@ -129,10 +129,10 @@ subroutine cice_init2()
          file=__FILE__,line= __LINE__)
 
     if (tr_fsd) call icepack_init_fsd_bounds ( &
-         floe_rad_l = floe_rad_l,    &  ! fsd size lower bound in m (radius)
-         floe_rad_c = floe_rad_c,    &  ! fsd size bin centre in m (radius)
-         floe_binwidth = floe_binwidth, &  ! fsd size bin width in m (radius)
-         c_fsd_range = c_fsd_range,   &  ! string for history output
+         floe_rad_l_out = floe_rad_l,    &  ! fsd size lower bound in m (radius)
+         floe_rad_c_out = floe_rad_c,    &  ! fsd size bin centre in m (radius)
+         floe_binwidth_out = floe_binwidth, &  ! fsd size bin width in m (radius)
+         c_fsd_range_out = c_fsd_range,   &  ! string for history output
          write_diags=(my_task == master_task))  ! write diag on master only
     call icepack_warnings_flush(nu_diag)
     if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &

cicecore/drivers/nuopc/dmi/CICE_InitMod.F90

@@ -167,10 +167,10 @@ subroutine cice_init(mpi_comm)
           file=__FILE__,line= __LINE__)
 
       if (tr_fsd) call icepack_init_fsd_bounds ( &
-         floe_rad_l = floe_rad_l,    &  ! fsd size lower bound in m (radius)
-         floe_rad_c = floe_rad_c,    &  ! fsd size bin centre in m (radius)
-         floe_binwidth = floe_binwidth, &  ! fsd size bin width in m (radius)
-         c_fsd_range = c_fsd_range,   &  ! string for history output
+         floe_rad_l_out = floe_rad_l,    &  ! fsd size lower bound in m (radius)
+         floe_rad_c_out = floe_rad_c,    &  ! fsd size bin centre in m (radius)
+         floe_binwidth_out = floe_binwidth, &  ! fsd size bin width in m (radius)
+         c_fsd_range_out = c_fsd_range,   &  ! string for history output
          write_diags=(my_task == master_task))  ! write diag on master only
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &

cicecore/drivers/standalone/cice/CICE_InitMod.F90

@@ -155,10 +155,10 @@ subroutine cice_init
           file=__FILE__,line= __LINE__)
 
       if (tr_fsd) call icepack_init_fsd_bounds ( &
-         floe_rad_l = floe_rad_l,    &  ! fsd size lower bound in m (radius)
-         floe_rad_c = floe_rad_c,    &  ! fsd size bin centre in m (radius)
-         floe_binwidth = floe_binwidth, &  ! fsd size bin width in m (radius)
-         c_fsd_range = c_fsd_range,   &  ! string for history output
+         floe_rad_l_out = floe_rad_l,    &  ! fsd size lower bound in m (radius)
+         floe_rad_c_out = floe_rad_c,    &  ! fsd size bin centre in m (radius)
+         floe_binwidth_out = floe_binwidth, &  ! fsd size bin width in m (radius)
+         c_fsd_range_out = c_fsd_range,   &  ! string for history output
          write_diags=(my_task == master_task))  ! write diag on master only
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &

cicecore/drivers/unittest/gridavgchk/CICE_InitMod.F90

@@ -155,10 +155,10 @@ subroutine cice_init
           file=__FILE__,line= __LINE__)
 
       if (tr_fsd) call icepack_init_fsd_bounds ( &
-         floe_rad_l = floe_rad_l,    &  ! fsd size lower bound in m (radius)
-         floe_rad_c = floe_rad_c,    &  ! fsd size bin centre in m (radius)
-         floe_binwidth = floe_binwidth, &  ! fsd size bin width in m (radius)
-         c_fsd_range = c_fsd_range,   &  ! string for history output
+         floe_rad_l_out = floe_rad_l,    &  ! fsd size lower bound in m (radius)
+         floe_rad_c_out = floe_rad_c,    &  ! fsd size bin centre in m (radius)
+         floe_binwidth_out = floe_binwidth, &  ! fsd size bin width in m (radius)
+         c_fsd_range_out = c_fsd_range,   &  ! string for history output
          write_diags=(my_task == master_task))  ! write diag on master only
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &

cicecore/drivers/unittest/halochk/CICE_InitMod.F90

@@ -155,10 +155,10 @@ subroutine cice_init
           file=__FILE__,line= __LINE__)
 
       if (tr_fsd) call icepack_init_fsd_bounds ( &
-         floe_rad_l = floe_rad_l,    &  ! fsd size lower bound in m (radius)
-         floe_rad_c = floe_rad_c,    &  ! fsd size bin centre in m (radius)
-         floe_binwidth = floe_binwidth, &  ! fsd size bin width in m (radius)
-         c_fsd_range = c_fsd_range,   &  ! string for history output
+         floe_rad_l_out = floe_rad_l,    &  ! fsd size lower bound in m (radius)
+         floe_rad_c_out = floe_rad_c,    &  ! fsd size bin centre in m (radius)
+         floe_binwidth_out = floe_binwidth, &  ! fsd size bin width in m (radius)
+         c_fsd_range_out = c_fsd_range,   &  ! string for history output
          write_diags=(my_task == master_task))  ! write diag on master only
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &

cicecore/drivers/unittest/opticep/CICE_InitMod.F90

@@ -155,10 +155,10 @@ subroutine cice_init
           file=__FILE__,line= __LINE__)
 
       if (tr_fsd) call icepack_init_fsd_bounds ( &
-         floe_rad_l = floe_rad_l,    &  ! fsd size lower bound in m (radius)
-         floe_rad_c = floe_rad_c,    &  ! fsd size bin centre in m (radius)
-         floe_binwidth = floe_binwidth, &  ! fsd size bin width in m (radius)
-         c_fsd_range = c_fsd_range,   &  ! string for history output
+         floe_rad_l_out = floe_rad_l,    &  ! fsd size lower bound in m (radius)
+         floe_rad_c_out = floe_rad_c,    &  ! fsd size bin centre in m (radius)
+         floe_binwidth_out = floe_binwidth, &  ! fsd size bin width in m (radius)
+         c_fsd_range_out = c_fsd_range,   &  ! string for history output
          write_diags=(my_task == master_task))  ! write diag on master only
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &

cicecore/drivers/unittest/opticep/ice_init_column.F90

@@ -581,8 +581,7 @@ end subroutine init_snowtracers
 
       subroutine init_fsd(floesize)
 
-      use ice_arrays_column, only: floe_rad_c, floe_binwidth, &
-         wavefreq, dwavefreq, wave_sig_ht, wave_spectrum, &
+      use ice_arrays_column, only: wavefreq, dwavefreq, wave_sig_ht, wave_spectrum, &
          d_afsd_newi, d_afsd_latg, d_afsd_latm, d_afsd_wave, d_afsd_weld
       use ice_domain_size, only: ncat, max_blocks, nfsd
       use ice_init, only: ice_ic
@@ -636,8 +635,6 @@ subroutine init_fsd(floesize)
 
          ! initialize floe size distribution the same in every column and category
          call icepack_init_fsd(ice_ic = ice_ic, &
-            floe_rad_c    = floe_rad_c,    &  ! fsd size bin centre in m (radius)
-            floe_binwidth = floe_binwidth, &  ! fsd size bin width in m (radius)
             afsd          = afsd)             ! floe size distribution
 
          do iblk = 1, max_blocks