icepack_parameters: optionally compute 'dragio' from under-ice roughn…

…ess (#366)

* 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_layer1', where 'iceruf_ocn' is the under-ice roughness
length, and 'thickness_ocn_layer1' 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

* Add 'calc_dragio' to namelist and add a test

Add the 'calc_dragio' flag introduced in the previous commit to the
Icepack namelist.

Add a test for the new feature in the base_suite.

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_layer1 = 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_layer1
+
 !-----------------------------------------------------------------------
 ! 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_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, &
@@ -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_layer1_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_layer1
+
 !-----------------------------------------------------------------------
 ! 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_layer1_in) ) thickness_ocn_layer1 = thickness_ocn_layer1_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_layer1_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_layer1_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_layer1
+
 !-----------------------------------------------------------------------
 ! 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_layer1_out) ) thickness_ocn_layer1_out = thickness_ocn_layer1
+      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_layer1 = ",thickness_ocn_layer1
         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_layer1/iceruf_ocn))**2 ! dragio at half first layer
+           lambda = (thickness_ocn_layer1 - iceruf_ocn) / &
+                    (thickness_ocn_layer1*(sqrt(dragio)/vonkar*(log(c2) - c1 + iceruf_ocn/thickness_ocn_layer1) + c1))
+           dragio = dragio*lambda**2
+        endif
+
       end subroutine icepack_recompute_constants
 
 !=======================================================================

configuration/driver/icedrv_init.F90

@@ -102,7 +102,7 @@ subroutine input_data
       ! Flux convergence tolerance
       real (kind=dbl_kind) :: atmiter_conv
 
-      logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair
+      logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair, calc_dragio
       logical (kind=log_kind) :: conserv_check
 
       integer (kind=int_kind) :: ntrcr
@@ -158,7 +158,7 @@ subroutine input_data
         update_ocn_f,    l_mpond_fresh,   ustar_min,       &
         fbot_xfer_type,  oceanmixed_ice,  emissivity,      &
         formdrag,        highfreq,        natmiter,        &
-        atmiter_conv,                        &
+        atmiter_conv,    calc_dragio,                      &
         tfrz_option,     default_season,  wave_spec_type,  &
         precip_units,    fyear_init,      ycycle,          &
         atm_data_type,   ocn_data_type,   bgc_data_type,   &
@@ -186,7 +186,7 @@ subroutine input_data
            albicev_out=albicev, albicei_out=albicei, ksno_out = ksno,   &
            albsnowv_out=albsnowv, albsnowi_out=albsnowi, &
            natmiter_out=natmiter, ahmax_out=ahmax, shortwave_out=shortwave, &
-           atmiter_conv_out = atmiter_conv, &
+           atmiter_conv_out = atmiter_conv, calc_dragio_out=calc_dragio, &
            albedo_type_out=albedo_type, R_ice_out=R_ice, R_pnd_out=R_pnd, &
            R_snw_out=R_snw, dT_mlt_out=dT_mlt, rsnw_mlt_out=rsnw_mlt, &
            kstrength_out=kstrength, krdg_partic_out=krdg_partic, &
@@ -578,6 +578,7 @@ subroutine input_data
          write(nu_diag,1005) ' atmiter_conv              = ', atmiter_conv
          write(nu_diag,1010) ' calc_strair               = ', calc_strair
          write(nu_diag,1010) ' calc_Tsfc                 = ', calc_Tsfc
+         write(nu_diag,1010) ' calc_dragio               = ', calc_dragio
          write(nu_diag,1005) ' floediam                  = ', floediam
          write(nu_diag,1005) ' hfrazilmin                = ', hfrazilmin
 
@@ -757,7 +758,7 @@ subroutine input_data
            albicev_in=albicev, albicei_in=albicei, ksno_in=ksno, &
            albsnowv_in=albsnowv, albsnowi_in=albsnowi, &
            natmiter_in=natmiter, ahmax_in=ahmax, shortwave_in=shortwave, &
-           atmiter_conv_in = atmiter_conv, &
+           atmiter_conv_in = atmiter_conv, calc_dragio_in=calc_dragio, &
            albedo_type_in=albedo_type, R_ice_in=R_ice, R_pnd_in=R_pnd, &
            R_snw_in=R_snw, dT_mlt_in=dT_mlt, rsnw_mlt_in=rsnw_mlt, &
            kstrength_in=kstrength, krdg_partic_in=krdg_partic, &

configuration/scripts/icepack_in

@@ -87,6 +87,7 @@
     natmiter        = 5
     atmiter_conv    = 0.0d0
     ustar_min       = 0.0005
+    calc_dragio     = .false.
     emissivity      = 0.985
     fbot_xfer_type  = 'constant'
     update_ocn_f    = .false.

configuration/scripts/options/set_nml.calcdragio

@@ -0,0 +1 @@
+calc_dragio = .true.

configuration/scripts/tests/base_suite.ts

@@ -16,6 +16,7 @@ smoke          col     1x1        debug,run1year,leap,dt30min
 smoke          col     1x1        debug,run1year,dyn
 smoke          col     1x1        debug,run1year,fsd12
 smoke          col     1x1        debug,run1year,fsd1
+smoke          col     1x1        debug,run1year,calcdragio
 restart        col     1x1        debug
 restart        col     1x1        diag1
 restart        col     1x1        pondcesm

doc/source/icepack_index.rst

@@ -66,6 +66,7 @@ either Celsius or Kelvin units).
    "bphi", "porosity of ice layers on bio grid", ""
    "bTiz", "temperature of ice layers on bio grid", ""
    "**C**", "", ""
+   "calc_dragio", ":math:`\bullet` if true, calculate ``dragio`` from ``iceruf_ocn`` and ``thickness_ocn_layer1``", "F"
    "calc_strair", ":math:`\bullet` if true, calculate wind stress", "T"
    "calc_Tsfc", ":math:`\bullet` if true, calculate surface temperature", "T"
    "Cdn_atm", "atmospheric drag coefficient", ""
@@ -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_layer1", "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

@@ -26,6 +26,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*"
@@ -301,6 +302,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

@@ -731,7 +731,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_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, &
@@ -751,7 +751,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, &
@@ -909,6 +909,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_layer1_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)
 
@@ -924,6 +926,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_layer1
+
   !-----------------------------------------------------------------------
   ! Parameters for atmosphere
   !-----------------------------------------------------------------------
@@ -1079,7 +1084,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_layer1_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, &
@@ -1099,7 +1104,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, &
@@ -1266,6 +1271,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_layer1_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)
 
@@ -1281,6 +1288,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_layer1
+
   !-----------------------------------------------------------------------
   ! Parameters for atmosphere
   !-----------------------------------------------------------------------
@@ -2855,7 +2865,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)