Update update_ocn_f implementation, Add cpl_frazil namelist (CICE-Con…

…sortium#889)

* Update update_ocn_f implementation

Add cpl_frazil namelist

Add update_ocn_f and cpl_frazil to icepack_init_parameters call, set these
values inside Icepack at initialization.

Remove update_ocn_f argument from icepack_step_therm2 call

Update runtime_diags and accum_hist to account for new Icepack and
cpl_frazil implementation.  These may need an addition update later.

* Update documentation

cicecore/cicedyn/analysis/ice_diagnostics.F90

@@ -123,7 +123,7 @@ subroutine runtime_diags (dt)
       use ice_flux, only: alvdr, alidr, alvdf, alidf, evap, fsnow, frazil, &
           fswabs, fswthru, flw, flwout, fsens, fsurf, flat, frzmlt_init, frain, fpond, &
           fhocn_ai, fsalt_ai, fresh_ai, frazil_diag, &
-          update_ocn_f, Tair, Qa, fsw, fcondtop, meltt, meltb, meltl, snoice, &
+          update_ocn_f, cpl_frazil, Tair, Qa, fsw, fcondtop, meltt, meltb, meltl, snoice, &
           dsnow, congel, sst, sss, Tf, fhocn, &
           swvdr, swvdf, swidr, swidf, &
           alvdr_init, alvdf_init, alidr_init, alidf_init
@@ -722,8 +722,9 @@ subroutine runtime_diags (dt)
          ! frazil ice growth !! should not be multiplied by aice
          ! m/step->kg/m^2/s
          work1(:,:,:) = frazil(:,:,:)*rhoi/dt
-         if (ktherm == 2 .and. .not.update_ocn_f) &
+         if (.not. update_ocn_f .and. ktherm == 2 .and. cpl_frazil == 'fresh_ice_correction') then
             work1(:,:,:) = (frazil(:,:,:)-frazil_diag(:,:,:))*rhoi/dt
+         endif
          frzn = c0
          frzs = c0
          frzn = global_sum(work1, distrb_info, &

cicecore/cicedyn/analysis/ice_history.F90

@@ -2134,11 +2134,10 @@ subroutine accum_hist (dt)
           fcondtop, fcondbot, fsurfn, fcondtopn, flatn, fsensn, albcnt, snwcnt, &
           stressp_1, stressm_1, stress12_1, &
           stresspT, stressmT, stress12T, &
-          stressp_2, &
-          stressp_3, &
-          stressp_4, sig1, sig2, sigP, &
+          stressp_2, stressp_3, stressp_4, sig1, sig2, sigP, &
           mlt_onset, frz_onset, dagedtt, dagedtd, fswint_ai, keffn_top, &
-          snowfrac, alvdr_ai, alvdf_ai, alidr_ai, alidf_ai, update_ocn_f
+          snowfrac, alvdr_ai, alvdf_ai, alidr_ai, alidf_ai, update_ocn_f, &
+          cpl_frazil
       use ice_arrays_column, only: snowfracn, Cdn_atm
       use ice_history_shared ! almost everything
       use ice_history_write, only: ice_write_hist
@@ -3238,11 +3237,11 @@ subroutine accum_hist (dt)
               if (aice(i,j,iblk) > puny) then
 !                Add in frazil flux
                  if (.not. update_ocn_f) then
-                 if ( ktherm == 2) then
-                    dfresh = -rhoi*(frazil(i,j,iblk)-frazil_diag(i,j,iblk))/dt
-                 else
-                    dfresh = -rhoi*frazil(i,j,iblk)/dt
-                 endif
+                    if ( ktherm == 2 .and. cpl_frazil == 'fresh_ice_correction') then
+                       dfresh = -rhoi*(frazil(i,j,iblk)-frazil_diag(i,j,iblk))/dt
+                    else
+                       dfresh = -rhoi*frazil(i,j,iblk)/dt
+                    endif
                  endif
                  if (saltflux_option == 'prognostic') then
                     sicen = c0
@@ -3266,14 +3265,13 @@ subroutine accum_hist (dt)
            do j = jlo, jhi
            do i = ilo, ihi
               if (aice(i,j,iblk) > puny) then
-!                Add in frazil flux
 !                Add in frazil flux
                  if (.not. update_ocn_f) then
-                 if ( ktherm == 2) then
-                    dfresh = -rhoi*(frazil(i,j,iblk)-frazil_diag(i,j,iblk))/dt
-                 else
-                    dfresh = -rhoi*frazil(i,j,iblk)/dt
-                 endif
+                    if ( ktherm == 2 .and. cpl_frazil == 'fresh_ice_correction') then
+                       dfresh = -rhoi*(frazil(i,j,iblk)-frazil_diag(i,j,iblk))/dt
+                    else
+                       dfresh = -rhoi*frazil(i,j,iblk)/dt
+                    endif
                  endif
                  worka(i,j) = aice(i,j,iblk)*(fresh(i,j,iblk)+dfresh)
               endif

cicecore/cicedyn/general/ice_flux.F90

@@ -241,8 +241,7 @@ module ice_flux
          alvdf_init, & ! visible, diffuse  (fraction)
          alidf_init    ! near-ir, diffuse  (fraction)
 
-      real (kind=dbl_kind), &
-         dimension(:,:,:,:), allocatable, public :: &
+      real (kind=dbl_kind), dimension(:,:,:,:), allocatable, public :: &
          albcnt       ! counter for zenith angle
 
        ! out to ocean
@@ -270,6 +269,9 @@ module ice_flux
          l_mpond_fresh   ! if true, include freshwater feedback from meltponds
                          ! when running in ice-ocean or coupled configuration
 
+      character (char_len), public :: &
+         cpl_frazil      ! type of coupling for frazil ice, 'fresh_ice_correction','internal','external'
+
       real (kind=dbl_kind), dimension (:,:,:,:), allocatable, public :: &
          meltsn      , & ! snow melt in category n (m)
          melttn      , & ! top melt in category n (m)

cicecore/cicedyn/general/ice_init.F90

@@ -82,7 +82,7 @@ subroutine input_data
       use ice_history_shared, only: hist_avg, history_dir, history_file, &
                              incond_dir, incond_file, version_name, &
                              history_precision, history_format, hist_time_axis
-      use ice_flux, only: update_ocn_f, l_mpond_fresh
+      use ice_flux, only: update_ocn_f, cpl_frazil, l_mpond_fresh
       use ice_flux, only: default_season
       use ice_flux_bgc, only: cpl_bgc
       use ice_forcing, only: &
@@ -265,7 +265,7 @@ subroutine input_data
         highfreq,       natmiter,        atmiter_conv,  calc_dragio,    &
         ustar_min,      emissivity,      iceruf,        iceruf_ocn,     &
         fbot_xfer_type, update_ocn_f,    l_mpond_fresh, tfrz_option,    &
-        saltflux_option,ice_ref_salinity,                               &
+        saltflux_option,ice_ref_salinity,cpl_frazil,                    &
         oceanmixed_ice, restore_ice,     restore_ocn,   trestore,       &
         precip_units,   default_season,  wave_spec_type,nfreq,          &
         atm_data_type,  ocn_data_type,   bgc_data_type, fe_data_type,   &
@@ -444,6 +444,7 @@ subroutine input_data
       ktransport = 1           ! -1 = off, 1 = on
       calc_Tsfc = .true.       ! calculate surface temperature
       update_ocn_f = .false.   ! include fresh water and salt fluxes for frazil
+      cpl_frazil = 'fresh_ice_correction' ! type of coupling for frazil ice
       ustar_min = 0.005        ! minimum friction velocity for ocean heat flux (m/s)
       hi_min = p01             ! minimum ice thickness allowed (m)
       iceruf = 0.0005_dbl_kind ! ice surface roughness at atmosphere interface (m)
@@ -1071,6 +1072,7 @@ subroutine input_data
       call broadcast_scalar(natmiter,             master_task)
       call broadcast_scalar(atmiter_conv,         master_task)
       call broadcast_scalar(update_ocn_f,         master_task)
+      call broadcast_scalar(cpl_frazil,           master_task)
       call broadcast_scalar(l_mpond_fresh,        master_task)
       call broadcast_scalar(ustar_min,            master_task)
       call broadcast_scalar(hi_min,               master_task)
@@ -2123,6 +2125,7 @@ subroutine input_data
             tmpstr2 = ' : frazil water/salt fluxes not included in ocean fluxes'
          endif
          write(nu_diag,1010) ' update_ocn_f     = ', update_ocn_f,trim(tmpstr2)
+         write(nu_diag,1030) ' cpl_frazil       = ', trim(cpl_frazil)
          if (l_mpond_fresh .and. tr_pond_topo) then
             tmpstr2 = ' : retain (topo) pond water until ponds drain'
          else
@@ -2510,8 +2513,8 @@ subroutine input_data
          floediam_in=floediam, hfrazilmin_in=hfrazilmin, Tliquidus_max_in=Tliquidus_max, &
          aspect_rapid_mode_in=aspect_rapid_mode, dSdt_slow_mode_in=dSdt_slow_mode, &
          phi_c_slow_mode_in=phi_c_slow_mode, phi_i_mushy_in=phi_i_mushy, conserv_check_in=conserv_check, &
-         wave_spec_type_in = wave_spec_type, &
-         wave_spec_in=wave_spec, nfreq_in=nfreq, &
+         wave_spec_type_in = wave_spec_type, wave_spec_in=wave_spec, nfreq_in=nfreq, &
+         update_ocn_f_in=update_ocn_f, cpl_frazil_in=cpl_frazil, &
          tfrz_option_in=tfrz_option, kalg_in=kalg, fbot_xfer_type_in=fbot_xfer_type, &
          saltflux_option_in=saltflux_option, ice_ref_salinity_in=ice_ref_salinity, &
          Pstar_in=Pstar, Cstar_in=Cstar, iceruf_in=iceruf, iceruf_ocn_in=iceruf_ocn, calc_dragio_in=calc_dragio, &

cicecore/cicedyn/general/ice_step_mod.F90

@@ -617,7 +617,7 @@ subroutine step_therm2 (dt, iblk)
       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, &
-          update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside, fside, wlat, &
+          fsalt, Tf, sss, salinz, fhocn, rside, fside, wlat, &
           meltl, frazil_diag
       use ice_flux_bgc, only: flux_bio, faero_ocn, &
           fiso_ocn, HDO_ocn, H2_16O_ocn, H2_18O_ocn
@@ -708,7 +708,6 @@ subroutine step_therm2 (dt, iblk)
                       fresh      = fresh     (i,j,  iblk), &
                       fsalt      = fsalt     (i,j,  iblk), &
                       fhocn      = fhocn     (i,j,  iblk), &
-                      update_ocn_f = update_ocn_f,         &
                       bgrid      = bgrid,                  &
                       cgrid      = cgrid,                  &
                       igrid      = igrid,                  &

doc/source/cice_index.rst

@@ -116,6 +116,7 @@ section :ref:`tabnamelist`.
    "cosw", "cosine of the turning angle in water", "1."
    "coszen", "cosine of the zenith angle", ""
    "Cp", "proportionality constant for potential energy", "kg/m\ :math:`^2`/s\ :math:`^2`"
+   "cpl_frazil", ":math:`\bullet` type of frazil ice coupling", ""
    "cp_air", "specific heat of air", "1005.0 J/kg/K"
    "cp_ice", "specific heat of fresh ice", "2106. J/kg/K"
    "cp_ocn", "specific heat of sea water", "4218. J/kg/K"

doc/source/user_guide/ug_case_settings.rst

@@ -629,6 +629,9 @@ forcing_nml
    "``calc_strair``", "``.false.``", "read wind stress and speed from files", "``.true.``"
    "", "``.true.``", "calculate wind stress and speed", ""
    "``calc_Tsfc``", "logical", "calculate surface temperature", "``.true.``"
+   "``cpl_frazil``", "``external``", "frazil water/salt fluxes are handled outside of Icepack", "``fresh_ice_correction``"
+   "", "``fresh_ice_correction``", "correct fresh-ice frazil water/salt fluxes for mushy physics", ""
+   "", "``internal``", "send full frazil water/salt fluxes for mushy physics", ""
    "``default_season``", "``summer``", "forcing initial summer values", "``winter``"
    "", "``winter``", "forcing initial winter values", ""
    "``emissivity``", "real", "emissivity of snow and ice", "0.985"