icepack_parameters: optionally compute 'dragio' from under-ice roughness

In some applications, we want the ice-ocean drag coefficient 'dragio'
to be computed using an under-ice roughness length and the thickness of
the first ocean level instead of imposing its value directly.

Add a new boolean flag 'calc_dragio' and real parameters 'iceruf_ocn',
'thickness_ocn', where 'iceruf_ocn' is the under-ice roughness length,
and 'thickness_ocn' is the thickness of the first ocean level, to module
'icepack_parameters'.

With 'calc_dragio = .true.', compute the ice-ocean drag coefficient
following [1]. Update the documentation accordingly.

[1] https://doi.org/10.1002/2014JC010677

columnphysics/icepack_parameters.F90

@@ -224,10 +224,15 @@ module icepack_parameters
          Cstar    = 20._dbl_kind     ,&! constant in Hibler strength formula
                                        ! (kstrength = 0)
          dragio   = 0.00536_dbl_kind ,&! ice-ocn drag coefficient
+         thickness_ocn = 2.0_dbl_kind,&! thickness of first ocean level (m)
+         iceruf_ocn = 0.03_dbl_kind  ,&! under-ice roughness (m)
          gravit   = 9.80616_dbl_kind ,&! gravitational acceleration (m/s^2)
          mu_rdg = 3.0_dbl_kind ! e-folding scale of ridged ice, krdg_partic=1 (m^0.5)
                                        ! (krdg_redist = 1)
 
+      logical (kind=log_kind), public :: &
+         calc_dragio     = .false.     ! if true, calculate dragio from iceruf_ocn and thickness_ocn
+
 !-----------------------------------------------------------------------
 ! Parameters for atmosphere
 !-----------------------------------------------------------------------
@@ -383,7 +388,7 @@ subroutine icepack_init_parameters(   &
          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, albocn_in, gravit_in, viscosity_dyn_in, &
+         depressT_in, dragio_in, thickness_ocn_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, &
@@ -403,7 +408,7 @@ subroutine icepack_init_parameters(   &
          ahmax_in, R_ice_in, R_pnd_in, R_snw_in, dT_mlt_in, rsnw_mlt_in, &
          kalg_in, kstrength_in, krdg_partic_in, krdg_redist_in, mu_rdg_in, &
          atmbndy_in, calc_strair_in, formdrag_in, highfreq_in, natmiter_in, &
-         atmiter_conv_in, &
+         atmiter_conv_in, calc_dragio_in, &
          tfrz_option_in, kitd_in, kcatbound_in, hs0_in, frzpnd_in, &
          floeshape_in, wave_spec_in, wave_spec_type_in, nfreq_in, &
          dpscale_in, rfracmin_in, rfracmax_in, pndaspect_in, hs1_in, hp1_in, &
@@ -561,6 +566,8 @@ subroutine icepack_init_parameters(   &
          Pstar_in,      & ! constant in Hibler strength formula 
          Cstar_in,      & ! constant in Hibler strength formula 
          dragio_in,     & ! ice-ocn drag coefficient
+         thickness_ocn_in, & ! thickness of first ocean level (m)
+         iceruf_ocn_in, & ! under-ice roughness (m)
          gravit_in,     & ! gravitational acceleration (m/s^2)
          iceruf_in        ! ice surface roughness (m)
 
@@ -576,6 +583,9 @@ subroutine icepack_init_parameters(   &
          mu_rdg_in         ! gives e-folding scale of ridged ice (m^.5) 
                            ! (krdg_redist = 1) 
 
+      logical (kind=log_kind), intent(in), optional :: &
+         calc_dragio_in    ! if true, calculate dragio from iceruf_ocn and thickness_ocn
+
 !-----------------------------------------------------------------------
 ! Parameters for atmosphere
 !-----------------------------------------------------------------------
@@ -727,6 +737,9 @@ subroutine icepack_init_parameters(   &
       if (present(cp_ocn_in)            ) cp_ocn           = cp_ocn_in
       if (present(depressT_in)          ) depressT         = depressT_in
       if (present(dragio_in)            ) dragio           = dragio_in
+      if (present(iceruf_ocn_in)        ) iceruf_ocn       = iceruf_ocn_in
+      if (present(thickness_ocn_in)     ) thickness_ocn    = thickness_ocn_in
+      if (present(calc_dragio_in)       ) calc_dragio      = calc_dragio_in
       if (present(albocn_in)            ) albocn           = albocn_in
       if (present(gravit_in)            ) gravit           = gravit_in
       if (present(viscosity_dyn_in)     ) viscosity_dyn    = viscosity_dyn_in
@@ -885,7 +898,7 @@ subroutine icepack_query_parameters(   &
          p333_out, p666_out, spval_const_out, pih_out, piq_out, pi2_out, &
          rhos_out, rhoi_out, rhow_out, cp_air_out, emissivity_out, &
          cp_ice_out, cp_ocn_out, hfrazilmin_out, floediam_out, &
-         depressT_out, dragio_out, albocn_out, gravit_out, viscosity_dyn_out, &
+         depressT_out, dragio_out, thickness_ocn_out, iceruf_ocn_out, albocn_out, gravit_out, viscosity_dyn_out, &
          Tocnfrz_out, rhofresh_out, zvir_out, vonkar_out, cp_wv_out, &
          stefan_boltzmann_out, ice_ref_salinity_out, &
          Tffresh_out, Lsub_out, Lvap_out, Timelt_out, Tsmelt_out, &
@@ -905,7 +918,7 @@ subroutine icepack_query_parameters(   &
          rsnw_mlt_out, dEdd_algae_out, &
          kalg_out, kstrength_out, krdg_partic_out, krdg_redist_out, mu_rdg_out, &
          atmbndy_out, calc_strair_out, formdrag_out, highfreq_out, natmiter_out, &
-         atmiter_conv_out, &
+         atmiter_conv_out, calc_dragio_out, &
          tfrz_option_out, kitd_out, kcatbound_out, hs0_out, frzpnd_out, &
          floeshape_out, wave_spec_out, wave_spec_type_out, nfreq_out, &
          dpscale_out, rfracmin_out, rfracmax_out, pndaspect_out, hs1_out, hp1_out, &
@@ -1072,6 +1085,8 @@ subroutine icepack_query_parameters(   &
          Pstar_out,      & ! constant in Hibler strength formula 
          Cstar_out,      & ! constant in Hibler strength formula 
          dragio_out,     & ! ice-ocn drag coefficient
+         thickness_ocn_out, & ! thickness of first ocean level (m)
+         iceruf_ocn_out, & ! under-ice roughness (m)
          gravit_out,     & ! gravitational acceleration (m/s^2)
          iceruf_out        ! ice surface roughness (m)
 
@@ -1087,6 +1102,9 @@ subroutine icepack_query_parameters(   &
          mu_rdg_out         ! gives e-folding scale of ridged ice (m^.5) 
                             ! (krdg_redist = 1) 
 
+      logical (kind=log_kind), intent(out), optional :: &
+         calc_dragio_out    ! if true, compute dragio from iceruf_ocn and thickness_ocn
+
 !-----------------------------------------------------------------------
 ! Parameters for atmosphere
 !-----------------------------------------------------------------------
@@ -1279,6 +1297,9 @@ subroutine icepack_query_parameters(   &
       if (present(cp_ocn_out)            ) cp_ocn_out       = cp_ocn
       if (present(depressT_out)          ) depressT_out     = depressT
       if (present(dragio_out)            ) dragio_out       = dragio
+      if (present(iceruf_ocn_out)        ) iceruf_ocn_out   = iceruf_ocn
+      if (present(thickness_ocn_out)     ) thickness_ocn_out= thickness_ocn
+      if (present(calc_dragio_out)       ) calc_dragio_out  = calc_dragio
       if (present(albocn_out)            ) albocn_out       = albocn
       if (present(gravit_out)            ) gravit_out       = gravit
       if (present(viscosity_dyn_out)     ) viscosity_dyn_out= viscosity_dyn
@@ -1452,6 +1473,9 @@ subroutine icepack_write_parameters(iounit)
         write(iounit,*) "  cp_ocn = ",cp_ocn
         write(iounit,*) "  depressT = ",depressT
         write(iounit,*) "  dragio = ",dragio
+        write(iounit,*) "  calc_dragio = ",calc_dragio
+        write(iounit,*) "  iceruf_ocn = ",iceruf_ocn
+        write(iounit,*) "  thickness_ocn = ",thickness_ocn
         write(iounit,*) "  albocn = ",albocn
         write(iounit,*) "  gravit = ",gravit
         write(iounit,*) "  viscosity_dyn = ",viscosity_dyn
@@ -1609,6 +1633,8 @@ subroutine icepack_recompute_constants()
 
 !autodocument_end
 
+      real (kind=dbl_kind) :: lambda
+
       character(len=*),parameter :: subname='(icepack_recompute_constants)'
 
         cprho  = cp_ocn*rhow
@@ -1619,6 +1645,13 @@ subroutine icepack_recompute_constants()
         pi2    = c2*pi
         rad_to_deg = c180/pi
 
+        if (calc_dragio) then
+           dragio = (vonkar/log(p5 * thickness_ocn/iceruf_ocn))**2 ! dragio at half first layer
+           lambda = (thickness_ocn - iceruf_ocn) / &
+                    (thickness_ocn*(sqrt(dragio)/vonkar*(log(c2) - c1 + iceruf_ocn/thickness_ocn) + c1))
+           dragio = dragio*lambda**2
+        endif
+
       end subroutine icepack_recompute_constants
 
 !=======================================================================

doc/source/icepack_index.rst

@@ -68,6 +68,7 @@ either Celsius or Kelvin units).
    "**C**", "", ""
    "calc_strair", ":math:`\bullet` if true, calculate wind stress", "T"
    "calc_Tsfc", ":math:`\bullet` if true, calculate surface temperature", "T"
+   "calc_dragio", "if true, calculate ``dragio`` from ``iceruf_ocn`` and ``thickness_ocn``", "F"
    "Cdn_atm", "atmospheric drag coefficient", ""
    "Cdn_ocn", "ocean drag coefficient", "" 
    "Cf", ":math:`\bullet` ratio of ridging work to PE change in ridging", "17."
@@ -227,6 +228,7 @@ either Celsius or Kelvin units).
    "ice_stderr", "unit number for standard error output", ""
    "ice_ref_salinity", "reference salinity for ice–ocean exchanges", "4. ppt"
    "iceruf", "ice surface roughness", "5.\ :math:`\times`\ 10\ :math:`^{-4}` m"
+   "iceruf_ocn", "under-ice roughness", "0.03 m"
    "idate", "the date at the end of the current time step (yyyymmdd)", ""
    "idate0", "initial date", ""
    "igrid", "interface points for vertical bio grid", ""
@@ -446,6 +448,7 @@ either Celsius or Kelvin units).
    "Tf", "freezing temperature", "C"
    "Tffresh", "freezing temp of fresh ice", "273.15 K"
    "tfrz_option", ":math:`\bullet` form of ocean freezing temperature", ""
+   "thickness_ocn", "thickness of first ocean level", "2.0 m"
    "thinS", "minimum ice thickness for brine tracer", ""
    "time", "total elapsed time", "s"
    "time_forc", "time of last forcing update", "s"

doc/source/master_list.bib

@@ -606,6 +606,17 @@ @Article{Roberts15
   pages = {211-228},
   url = {http://dx.doi.org/10.3189/2015AoG69A760}
 }
+@article{Roy15,
+author = "F. Roy and M. Chevallier and G.C. Smith and F. Dupont and G. Garric and J.-F. Lemieux and Y. Lu and F. Davidson",
+title = {Arctic sea ice and freshwater sensitivity to the treatment of the atmosphere-ice-ocean surface layer},
+journal = JGRO,
+volume = {120},
+number = {6},
+pages = {4392-4417},
+doi = {https://doi.org/10.1002/2014JC010677},
+url = {https://doi.org/10.1002/2014JC010677},
+year = {2015}
+}
 @Article{Duarte17
   author = "P. Duarte and Coauthors",
   title = "{Sea ice thermohaline dynamics and biogeochemistry in the Arctic Ocean: Empirical and model results}",

doc/source/science_guide/sg_boundary_forcing.rst

@@ -25,6 +25,7 @@ Atmosphere and ocean boundary forcing
    ":math:`T_w`", "Sea surface temperature", "From *ocean model* to *sea ice model*"
    ":math:`S`", "Sea surface salinity", "From *ocean model* to *sea ice model*"
    ":math:`\nabla H_o`", "Sea surface slope", "From *ocean model* via flux coupler to *sea ice model*"
+   ":math:`h_1`", "Thickness of first ocean level (optional)", "From *ocean model* to *sea ice model*"
    ":math:`\vec{U}_w`", "Surface ocean currents", "From *ocean model* to *sea ice model* (available in Icepack driver, not used directly in column physics)"
    ":math:`\vec{\tau}_a`", "Wind stress", "From *sea ice model* to *atmosphere model*"
    ":math:`F_s`", "Sensible heat flux", "From *sea ice model* to *atmosphere model*"
@@ -296,6 +297,22 @@ the ocean. If the resulting sea surface temperature falls below the
 salinity-dependent freezing point, then new ice (frazil) forms.
 Otherwise, heat is made available for melting the ice.
 
+The ice-ocean drag coefficient, :math:`c_w`, can optionally be computed from the thickness of the first ocean level, :math:`h_1`, and an under-ice roughness length, :math:`z_{io}`.
+The computation follows :cite:`Roy15` :
+
+.. math::
+   c_w = c_w^* \lambda^2
+   :label: dragio
+
+where
+
+.. math::
+   \begin{aligned}
+   c_w^* &= \frac{\kappa^2} {\ln^2\left( h_1 /  z_{io} \right)}, \\
+   \lambda &= \frac{h_1 - z_{io}} {h_1 \left[ \sqrt{c_w^*} \kappa^{-1} \left( \ln(2) - 1 +  z_{io} / h_1 \right) + 1 \right] }
+   \end{aligned}
+   :label: dragio-defs
+
 .. _formdrag:
 
 Variable exchange coefficients

doc/source/user_guide/interfaces.include

@@ -730,7 +730,7 @@ icepack_init_parameters
            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, albocn_in, gravit_in, viscosity_dyn_in, &
+           depressT_in, dragio_in, thickness_ocn_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, &
@@ -750,7 +750,7 @@ icepack_init_parameters
            ahmax_in, R_ice_in, R_pnd_in, R_snw_in, dT_mlt_in, rsnw_mlt_in, &
            kalg_in, kstrength_in, krdg_partic_in, krdg_redist_in, mu_rdg_in, &
            atmbndy_in, calc_strair_in, formdrag_in, highfreq_in, natmiter_in, &
-           atmiter_conv_in, &
+           atmiter_conv_in, calc_dragio_in, &
            tfrz_option_in, kitd_in, kcatbound_in, hs0_in, frzpnd_in, &
            floeshape_in, wave_spec_in, wave_spec_type_in, nfreq_in, &
            dpscale_in, rfracmin_in, rfracmax_in, pndaspect_in, hs1_in, hp1_in, &
@@ -908,6 +908,8 @@ icepack_init_parameters
            Pstar_in,      & ! constant in Hibler strength formula 
            Cstar_in,      & ! constant in Hibler strength formula 
            dragio_in,     & ! ice-ocn drag coefficient
+           thickness_ocn_in, & ! thickness of first ocean level (m)
+           iceruf_ocn_in, & ! under-ice roughness (m)
            gravit_in,     & ! gravitational acceleration (m/s^2)
            iceruf_in        ! ice surface roughness (m)
 
@@ -923,6 +925,9 @@ icepack_init_parameters
            mu_rdg_in         ! gives e-folding scale of ridged ice (m^.5) 
                              ! (krdg_redist = 1) 
 
+        logical (kind=log_kind), intent(in), optional :: &
+           calc_dragio_in    ! if true, calculate dragio from iceruf_ocn and thickness_ocn
+
   !-----------------------------------------------------------------------
   ! Parameters for atmosphere
   !-----------------------------------------------------------------------
@@ -1078,7 +1083,7 @@ icepack_query_parameters
            p333_out, p666_out, spval_const_out, pih_out, piq_out, pi2_out, &
            rhos_out, rhoi_out, rhow_out, cp_air_out, emissivity_out, &
            cp_ice_out, cp_ocn_out, hfrazilmin_out, floediam_out, &
-           depressT_out, dragio_out, albocn_out, gravit_out, viscosity_dyn_out, &
+           depressT_out, dragio_out, thickness_ocn_out, iceruf_ocn_out, albocn_out, gravit_out, viscosity_dyn_out, &
            Tocnfrz_out, rhofresh_out, zvir_out, vonkar_out, cp_wv_out, &
            stefan_boltzmann_out, ice_ref_salinity_out, &
            Tffresh_out, Lsub_out, Lvap_out, Timelt_out, Tsmelt_out, &
@@ -1098,7 +1103,7 @@ icepack_query_parameters
            rsnw_mlt_out, dEdd_algae_out, &
            kalg_out, kstrength_out, krdg_partic_out, krdg_redist_out, mu_rdg_out, &
            atmbndy_out, calc_strair_out, formdrag_out, highfreq_out, natmiter_out, &
-           atmiter_conv_out, &
+           atmiter_conv_out, calc_dragio_out, &
            tfrz_option_out, kitd_out, kcatbound_out, hs0_out, frzpnd_out, &
            floeshape_out, wave_spec_out, wave_spec_type_out, nfreq_out, &
            dpscale_out, rfracmin_out, rfracmax_out, pndaspect_out, hs1_out, hp1_out, &
@@ -1265,6 +1270,8 @@ icepack_query_parameters
            Pstar_out,      & ! constant in Hibler strength formula 
            Cstar_out,      & ! constant in Hibler strength formula 
            dragio_out,     & ! ice-ocn drag coefficient
+           thickness_ocn_out, & ! thickness of first ocean level (m)
+           iceruf_ocn_out, & ! under-ice roughness (m)
            gravit_out,     & ! gravitational acceleration (m/s^2)
            iceruf_out        ! ice surface roughness (m)
 
@@ -1280,6 +1287,9 @@ icepack_query_parameters
            mu_rdg_out         ! gives e-folding scale of ridged ice (m^.5) 
                               ! (krdg_redist = 1) 
 
+        logical (kind=log_kind), intent(out), optional :: &
+           calc_dragio_out    ! if true, compute dragio from iceruf_ocn and thickness_ocn
+
   !-----------------------------------------------------------------------
   ! Parameters for atmosphere
   !-----------------------------------------------------------------------
@@ -2852,7 +2862,6 @@ icepack_init_bgc
 .. code-block:: fortran
 
   !
-
         subroutine icepack_init_bgc(ncat, nblyr, nilyr, ntrcr_o, &
            cgrid, igrid, ntrcr, nbtrcr, &
            sicen, trcrn, sss, ocean_bio_all)