generalize logic for coupling frazil fluxes (#458)

* cherry-pick update_ocn_f commit: generalize update_ocn_f logic

* update interfaces.include

* declare update_ocn_f as logical

* move optional arg check and alphabetize doc

columnphysics/icepack_parameters.F90

@@ -147,8 +147,9 @@ module icepack_parameters
                          ! 2 = mushy layer theory
 
       character (char_len), public :: &
-         conduct = 'bubbly', &      ! 'MU71' or 'bubbly'
-         fbot_xfer_type = 'constant' ! transfer coefficient type for ice-ocean heat flux
+         conduct = 'bubbly', &          ! 'MU71' or 'bubbly'
+         fbot_xfer_type = 'constant', & ! transfer coefficient type for ice-ocean heat flux
+         cpl_frazil = 'fresh_ice_correction' ! type of coupling for frazil ice
 
       logical (kind=log_kind), public :: &
          calc_Tsfc     = .true. ,&! if true, calculate surface temperature
@@ -438,7 +439,8 @@ subroutine icepack_init_parameters(   &
          argcheck_in, puny_in, bignum_in, pi_in, secday_in, &
          rhos_in, rhoi_in, rhow_in, cp_air_in, emissivity_in, &
          cp_ice_in, cp_ocn_in, hfrazilmin_in, floediam_in, &
-         depressT_in, dragio_in, thickness_ocn_layer1_in, iceruf_ocn_in, albocn_in, gravit_in, viscosity_dyn_in, &
+         depressT_in, dragio_in, thickness_ocn_layer1_in, iceruf_ocn_in, &
+         albocn_in, gravit_in, viscosity_dyn_in, &
          Tocnfrz_in, rhofresh_in, zvir_in, vonkar_in, cp_wv_in, &
          stefan_boltzmann_in, ice_ref_salinity_in, &
          Tffresh_in, Lsub_in, Lvap_in, Timelt_in, Tsmelt_in, &
@@ -451,6 +453,7 @@ subroutine icepack_init_parameters(   &
          qqqice_in, TTTice_in, qqqocn_in, TTTocn_in, &
          ktherm_in, conduct_in, fbot_xfer_type_in, calc_Tsfc_in, dts_b_in, &
          update_ocn_f_in, ustar_min_in, a_rapid_mode_in, &
+         cpl_frazil_in, &
          Rac_rapid_mode_in, aspect_rapid_mode_in, &
          dSdt_slow_mode_in, phi_c_slow_mode_in, &
          phi_i_mushy_in, shortwave_in, albedo_type_in, albsnowi_in, &
@@ -543,8 +546,9 @@ subroutine icepack_init_parameters(   &
                             ! 2 = mushy layer theory
 
       character (len=*), intent(in), optional :: &
-         conduct_in, &      ! 'MU71' or 'bubbly'
-         fbot_xfer_type_in  ! transfer coefficient type for ice-ocean heat flux
+         conduct_in, &        ! 'MU71' or 'bubbly'
+         fbot_xfer_type_in, & ! transfer coefficient type for ice-ocean heat flux
+         cpl_frazil_in        ! type of coupling for frazil ice
 
       logical (kind=log_kind), intent(in), optional :: &
          calc_Tsfc_in    , &! if true, calculate surface temperature
@@ -896,6 +900,7 @@ subroutine icepack_init_parameters(   &
       if (present(conduct_in)           ) conduct          = conduct_in
       if (present(fbot_xfer_type_in)    ) fbot_xfer_type   = fbot_xfer_type_in
       if (present(calc_Tsfc_in)         ) calc_Tsfc        = calc_Tsfc_in
+      if (present(cpl_frazil_in)        ) cpl_frazil       = cpl_frazil_in
       if (present(update_ocn_f_in)      ) update_ocn_f     = update_ocn_f_in
       if (present(dts_b_in)             ) dts_b            = dts_b_in
       if (present(ustar_min_in)         ) ustar_min        = ustar_min_in
@@ -1156,7 +1161,7 @@ subroutine icepack_query_parameters(   &
          zref_out, hs_min_out, snowpatch_out, rhosi_out, sk_l_out, &
          saltmax_out, phi_init_out, min_salin_out, salt_loss_out, &
          min_bgc_out, dSin0_frazil_out, hi_ssl_out, hs_ssl_out, &
-         awtvdr_out, awtidr_out, awtvdf_out, awtidf_out, &
+         awtvdr_out, awtidr_out, awtvdf_out, awtidf_out, cpl_frazil_out, &
          qqqice_out, TTTice_out, qqqocn_out, TTTocn_out, update_ocn_f_out, &
          Lfresh_out, cprho_out, Cp_out, ustar_min_out, a_rapid_mode_out, &
          ktherm_out, conduct_out, fbot_xfer_type_out, calc_Tsfc_out, dts_b_out, &
@@ -1261,8 +1266,9 @@ subroutine icepack_query_parameters(   &
                             ! 2 = mushy layer theory
 
       character (len=*), intent(out), optional :: &
-         conduct_out, &     ! 'MU71' or 'bubbly'
-         fbot_xfer_type_out ! transfer coefficient type for ice-ocean heat flux
+         conduct_out, &        ! 'MU71' or 'bubbly'
+         fbot_xfer_type_out, & ! transfer coefficient type for ice-ocean heat flux
+         cpl_frazil_out        ! type of coupling for frazil ice
 
       logical (kind=log_kind), intent(out), optional :: &
          calc_Tsfc_out    ,&! if true, calculate surface temperature
@@ -1651,6 +1657,7 @@ subroutine icepack_query_parameters(   &
       if (present(conduct_out)           ) conduct_out      = conduct
       if (present(fbot_xfer_type_out)    ) fbot_xfer_type_out = fbot_xfer_type
       if (present(calc_Tsfc_out)         ) calc_Tsfc_out    = calc_Tsfc
+      if (present(cpl_frazil_out)        ) cpl_frazil_out   = cpl_frazil
       if (present(update_ocn_f_out)      ) update_ocn_f_out = update_ocn_f
       if (present(dts_b_out)             ) dts_b_out        = dts_b
       if (present(ustar_min_out)         ) ustar_min_out    = ustar_min
@@ -1856,6 +1863,7 @@ subroutine icepack_write_parameters(iounit)
         write(iounit,*) "  conduct    = ", trim(conduct)
         write(iounit,*) "  fbot_xfer_type = ", trim(fbot_xfer_type)
         write(iounit,*) "  calc_Tsfc  = ", calc_Tsfc
+        write(iounit,*) "  cpl_frazil = ", cpl_frazil
         write(iounit,*) "  update_ocn_f = ", update_ocn_f
         write(iounit,*) "  dts_b      = ", dts_b
         write(iounit,*) "  ustar_min  = ", ustar_min

columnphysics/icepack_therm_itd.F90

@@ -27,7 +27,7 @@ module icepack_therm_itd
       use icepack_parameters, only: rhosi, conserv_check, rhosmin, snwredist
       use icepack_parameters, only: kitd, ktherm
       use icepack_parameters, only: z_tracers, hfrazilmin
-      use icepack_parameters, only: saltflux_option
+      use icepack_parameters, only: cpl_frazil, update_ocn_f, saltflux_option
       use icepack_parameters, only: icepack_chkoptargflag
 
       use icepack_tracers, only: ntrcr, nbtrcr
@@ -1306,7 +1306,6 @@ subroutine add_new_ice (ncat,      nilyr,      &
                               aice0,     aice,       &
                               frzmlt,    frazil,     &
                               frz_onset, yday,       &
-                              update_ocn_f,          &
                               fresh,     fsalt,      &
                               Tf,        sss,        &
                               salinz,    phi_init,   &
@@ -1379,9 +1378,6 @@ subroutine add_new_ice (ncat,      nilyr,      &
          phi_init     , & ! initial frazil liquid fraction
          dSin0_frazil     ! initial frazil bulk salinity reduction from sss
 
-      logical (kind=log_kind), intent(in) :: &
-         update_ocn_f ! if true, update fresh water and salt fluxes
-
       ! BGC
       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
          bgrid              ! biology nondimensional vertical grid points
@@ -1617,29 +1613,29 @@ subroutine add_new_ice (ncat,      nilyr,      &
       !       is NOT included in fluxes fresh and fsalt.
       !-----------------------------------------------------------------
 
-      if (update_ocn_f) then
-         dfresh = -rhoi*vi0new/dt
+      dfresh = c0
+      dfsalt = c0
+      if (cpl_frazil == 'external') then
+         ! do nothing here, calculations are in the coupler or elsewhere
+      else
+         if (update_ocn_f) then
+            dfresh = -rhoi*vi0new/dt
+         elseif (cpl_frazil == 'fresh_ice_correction' .and. ktherm == 2) then
+            ! correct frazil fluxes for mushy
+            vi0tmp = fnew*dt / (rhoi*Lfresh) ! ocn/cpl assumes frazil volume is pure, fresh ice
+            dfresh = -rhoi*(vi0new - vi0tmp)/dt
+            frazil_diag = frazil - vi0tmp
+!        else 
+!           do nothing - other correction options could be implemented in the future             
+         endif
+
          if (saltflux_option == 'prognostic') then
             dfsalt = Si0new*p001*dfresh
          else
             dfsalt = ice_ref_salinity*p001*dfresh
          endif
          fresh  = fresh + dfresh
-         fsalt  = fsalt + dfsalt
-      else ! update_ocn_f = false
-         if (ktherm == 2) then ! return mushy-layer frazil to ocean (POP)
-            vi0tmp = fnew*dt / (rhoi*Lfresh)
-            dfresh = -rhoi*(vi0new - vi0tmp)/dt
-            if (saltflux_option == 'prognostic') then
-               dfsalt = Si0new*p001*dfresh
-            else
-               dfsalt = ice_ref_salinity*p001*dfresh
-            endif
-            fresh  = fresh + dfresh
-            fsalt  = fsalt + dfsalt
-            frazil_diag = frazil - vi0tmp
-         ! elseif ktherm==1 do nothing
-         endif
+         fsalt  = fsalt + dfsalt      
       endif
 
       !-----------------------------------------------------------------
@@ -2009,6 +2005,8 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr,           &
                                      d_afsd_latm,  d_afsd_weld,   &
                                      floe_rad_c,   floe_binwidth)
 
+      use icepack_parameters, only: icepack_init_parameters
+
       integer (kind=int_kind), intent(in) :: &
          ncat     , & ! number of thickness categories
          nltrcr   , & ! number of zbgc tracers
@@ -2019,7 +2017,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr,           &
       integer (kind=int_kind), intent(in), optional :: &
          nfsd         ! number of floe size categories
 
-      logical (kind=log_kind), intent(in) :: &
+      logical (kind=log_kind), intent(in), optional :: &
          update_ocn_f ! if true, update fresh water and salt fluxes
 
       real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
@@ -2139,6 +2137,9 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr,           &
       ! Check optional arguments and set local values
       !-----------------------------------------------------------------
 
+       if (present(update_ocn_f)) then
+          call icepack_init_parameters(update_ocn_f_in=update_ocn_f)
+       endif
        if (icepack_chkoptargflag(first_call)) then
           if (tr_iso) then
              if (.not.(present(fiso_ocn)   .and. &
@@ -2237,7 +2238,6 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr,           &
                            aice0,         aice,         &
                            frzmlt,        frazil,       &
                            frz_onset,     yday,         &
-                           update_ocn_f,                &
                            fresh,         fsalt,        &
                            Tf,            sss,          &
                            salinz,        phi_init,     &

configuration/driver/icedrv_flux.F90

@@ -19,7 +19,6 @@ module icedrv_flux
       use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
       use icepack_intfc, only: icepack_query_parameters
       use icepack_intfc, only: icepack_query_tracer_flags, icepack_query_tracer_indices
-      use icepack_intfc, only: icepack_query_parameters
       use icedrv_system, only: icedrv_system_abort
 
       implicit none
@@ -193,7 +192,6 @@ module icedrv_flux
          scale_factor! scaling factor for shortwave components
 
       logical (kind=log_kind), public :: &
-         update_ocn_f, & ! if true, update fresh water and salt fluxes
          l_mpond_fresh   ! if true, include freshwater feedback from meltponds
                          ! when running in ice-ocean or coupled configuration
 

configuration/driver/icedrv_init.F90

@@ -64,7 +64,7 @@ subroutine input_data
       use icedrv_calendar, only: npt, dt, ndtd, days_per_year, use_leap_years
       use icedrv_history, only: history_format
       use icedrv_restart_shared, only: restart, restart_dir, restart_file, restart_format
-      use icedrv_flux, only: update_ocn_f, l_mpond_fresh, cpl_bgc
+      use icedrv_flux, only: l_mpond_fresh, cpl_bgc
       use icedrv_flux, only: default_season
       use icedrv_forcing, only: precip_units,    fyear_init,      ycycle
       use icedrv_forcing, only: atm_data_type,   ocn_data_type,   bgc_data_type
@@ -98,9 +98,10 @@ subroutine input_data
          natmiter, kitd, kcatbound
 
       character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, &
-         tfrz_option, saltflux_option, frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table
+         cpl_frazil, tfrz_option, saltflux_option, &
+         frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table
 
-      logical (kind=log_kind) :: sw_redist, use_smliq_pnd, snwgrain
+      logical (kind=log_kind) :: sw_redist, use_smliq_pnd, snwgrain, update_ocn_f
       real (kind=dbl_kind)    :: sw_frac, sw_dtemp
 
       ! Flux convergence tolerance
@@ -174,7 +175,7 @@ subroutine input_data
         formdrag,        highfreq,        natmiter,        &
         atmiter_conv,    calc_dragio,                      &
         tfrz_option,     saltflux_option, ice_ref_salinity, &
-        default_season,  wave_spec_type,  &
+        default_season,  wave_spec_type,  cpl_frazil,      &
         precip_units,    fyear_init,      ycycle,          &
         atm_data_type,   ocn_data_type,   bgc_data_type,   &
         lateral_flux_type,                                &
@@ -216,7 +217,7 @@ subroutine input_data
            pndaspect_out=pndaspect, hs1_out=hs1, hp1_out=hp1, &
            ktherm_out=ktherm, calc_Tsfc_out=calc_Tsfc, &
            floediam_out=floediam, hfrazilmin_out=hfrazilmin, &
-           update_ocn_f_out = update_ocn_f, &
+           update_ocn_f_out = update_ocn_f, cpl_frazil_out = cpl_frazil, &
            conduct_out=conduct, a_rapid_mode_out=a_rapid_mode, &
            Rac_rapid_mode_out=Rac_rapid_mode, &
            aspect_rapid_mode_out=aspect_rapid_mode, &
@@ -788,6 +789,7 @@ subroutine input_data
          if (trim(atm_data_type)=='default') &
          write(nu_diag,*)    ' default_season            = ', trim(default_season)
 
+         write(nu_diag,1030) ' cpl_frazil                = ', trim(cpl_frazil)
          write(nu_diag,1010) ' update_ocn_f              = ', update_ocn_f
          write(nu_diag,1010) ' wave_spec                 = ', wave_spec
          if (wave_spec) &
@@ -980,6 +982,7 @@ subroutine input_data
            floediam_in=floediam, hfrazilmin_in=hfrazilmin, &
            ktherm_in=ktherm, calc_Tsfc_in=calc_Tsfc, &
            conduct_in=conduct, a_rapid_mode_in=a_rapid_mode, &
+           update_ocn_f_in=update_ocn_f, cpl_frazil_in=cpl_frazil, &
            Rac_rapid_mode_in=Rac_rapid_mode, &
            aspect_rapid_mode_in=aspect_rapid_mode, &
            dSdt_slow_mode_in=dSdt_slow_mode, &

configuration/driver/icedrv_step.F90

@@ -437,7 +437,7 @@ subroutine step_therm2 (dt)
       use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, nblyr, &
                                     nltrcr, nx, nfsd
       use icedrv_flux, only: fresh, frain, fpond, frzmlt, frazil, frz_onset
-      use icedrv_flux, only: update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside, fside, wlat
+      use icedrv_flux, only: fsalt, Tf, sss, salinz, fhocn, rside, fside, wlat
       use icedrv_flux, only: meltl, frazil_diag, flux_bio, faero_ocn, fiso_ocn
       use icedrv_flux, only: HDO_ocn, H2_16O_ocn, H2_18O_ocn
       use icedrv_init, only: tmask
@@ -505,7 +505,7 @@ subroutine step_therm2 (dt)
                          frzmlt=frzmlt(i), frazil=frazil(i),          &
                          frain=frain(i),   fpond=fpond(i),            &
                          fresh=fresh(i),   fsalt=fsalt(i),            &
-                         fhocn=fhocn(i),   update_ocn_f=update_ocn_f, &
+                         fhocn=fhocn(i),                              &
                          bgrid=bgrid,      cgrid=cgrid,               &
                          igrid=igrid,      faero_ocn=faero_ocn(i,:),  &
                          first_ice=first_ice(i,:),                    &

configuration/scripts/icepack_in

@@ -107,6 +107,7 @@
     calc_dragio     = .false.
     emissivity      = 0.985
     fbot_xfer_type  = 'constant'
+    cpl_frazil      = 'fresh_ice_correction'
     update_ocn_f    = .false.
     l_mpond_fresh   = .false.
     tfrz_option     = 'mushy'

configuration/scripts/options/set_nml.alt04

@@ -14,6 +14,7 @@
     krdg_redist     = 0
     calc_Tsfc       = .true.
     highfreq        = .true.
+    cpl_frazil      = 'internal'
     update_ocn_f    = .true.
     l_mpond_fresh   = .true.
     restore_ocn     = .true.

doc/source/icepack_index.rst

@@ -92,6 +92,7 @@ either Celsius or Kelvin units).  Deprecated parameters are listed at the end.
    "cp_wv", "specific heat of water vapor", "1.81x\ :math:`10^3` J/kg/K"
    "cp063", "diffuse fresnel reflectivity (above)", "0.063"
    "cp455", "diffuse fresnel reflectivity (below)", "0.455"
+   "cpl_frazil", ":math:`\bullet` type of frazil ice coupling", ""
    "Cs", "fraction of shear energy contributing to ridging", "0.25"
    "Cstar", "constant in Hibler ice strength formula", "20."
    "**D**", "", ""

doc/source/science_guide/sg_boundary_forcing.rst

@@ -278,9 +278,20 @@ In addition to runoff from rain and melted snow, the fresh water flux
 and water frozen (a negative flux) or melted at the bottom surface of
 the ice. This flux is computed as the net change of fresh water in the
 ice and snow volume over the coupling time step, excluding frazil ice
-formation and newly accumulated snow. Setting the namelist option
+formation and newly accumulated snow. 
+
+Setting the namelist option
 ``update_ocn_f`` to true causes frazil ice to be included in the fresh
-water and salt fluxes.
+water and salt fluxes. Some ocean models compute the frazil ice fluxes, 
+which then might need to be corrected for consistency with mushy physics.
+This behavior is controlled using a combination of ``update_ocn_f``, 
+``cpl_frazil`` and ``ktherm``.  In particular,
+``cpl_frazil = 'external'`` assumes that the frazil ice fluxes are 
+handled entirely outside of Icepack. When ``ktherm=2``, 
+``cpl_frazil = 'fresh_ice_correction'``
+sends coupling fluxes representing the difference between the mushy frazil 
+fluxes and fluxes computed assuming the frazil is purely fresh ice.
+Otherwise the internally computed frazil fluxes are sent to the coupler.
 
 There is a flux of salt into the ocean under melting conditions, and a
 (negative) flux when sea water is freezing. However, melting sea ice