icepack: optionally compute 'dragio' using under-ice roughness length (…

…CICE-Consortium#612)

* icepack: optionally compute 'dragio' using under-ice roughness length

In CICE-Consortium/Icepack@a80472b (icepack_parameters: optionally
compute 'dragio' from under-ice roughness (CICE-Consortium/Icepack#366),
2021-06-22), Icepack was updated to optionally compute the ice-ocean
drag coefficicent 'dragio' using an under-ice roughness length and the
thickness of the first ocean level.

Leverage this new feature in CICE by adding 'calc_dragio' and
'iceruf_ocn' to the CICE namelist. Add the new variables to the index in
the documentation and add a test with the new feature (using default
values for 'iceruf_ocn' and 'thickness_ocn_layer1').

As this new feature will mostly be useful in a coupled context, we do
not add 'thickness_ocn_layer1' to the namelist as it is expected that
the ocean model will pass this information to CICE.

* ice_grid: set 'thickness_ocn_layer1' if using 'calc_dragio'

In the previous commit we updated Icepack to allow computing the
ice-ocean drag coefficient 'dragio' using an under-ice roughness length
and the thickness of the first ocean layer, 'thickness_ocean_layer1', a
new Icepack parameter.

In some situations, we have access in CICE to the thicknesses of the
ocean levels, either hard-coded (use_bathymetry = false,
bathymetry_format = default), read from a file (use_bathymetry = true,
bathymetry_format = pop), or generated from the kmt_file
(use_bathymetry = false, bathymetry_format = pop). In these situations,
for consistency set 'thickness_ocean_layer1' in Icepack to the thickness
of the first ocean level, 'thick(1)' if 'calc_dragio' is active.

cicecore/cicedynB/general/ice_init.F90

@@ -128,7 +128,7 @@ subroutine input_data
         mu_rdg, hs0, dpscale, rfracmin, rfracmax, pndaspect, hs1, hp1, &
         a_rapid_mode, Rac_rapid_mode, aspect_rapid_mode, dSdt_slow_mode, &
         phi_c_slow_mode, phi_i_mushy, kalg, atmiter_conv, Pstar, Cstar, &
-        sw_frac, sw_dtemp, floediam, hfrazilmin, iceruf
+        sw_frac, sw_dtemp, floediam, hfrazilmin, iceruf, iceruf_ocn
 
       integer (kind=int_kind) :: ktherm, kstrength, krdg_partic, krdg_redist, natmiter, &
         kitd, kcatbound, ktransport
@@ -137,7 +137,7 @@ subroutine input_data
         tfrz_option, frzpnd, atmbndy, wave_spec_type
 
       logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair, wave_spec, &
-                                 sw_redist
+                                 sw_redist, calc_dragio
 
       logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_pond
       logical (kind=log_kind) :: tr_iso, tr_aero, tr_fsd
@@ -228,8 +228,8 @@ subroutine input_data
 
       namelist /forcing_nml/ &
         formdrag,       atmbndy,         calc_strair,   calc_Tsfc,      &
-        highfreq,       natmiter,        atmiter_conv,                  &
-        ustar_min,      emissivity,      iceruf,                        &
+        highfreq,       natmiter,        atmiter_conv,  calc_dragio,    &
+        ustar_min,      emissivity,      iceruf,        iceruf_ocn,     &
         fbot_xfer_type, update_ocn_f,    l_mpond_fresh, tfrz_option,    &
         oceanmixed_ice, restore_ice,     restore_ocn,   trestore,       &
         precip_units,   default_season,  wave_spec_type,nfreq,          &
@@ -384,6 +384,8 @@ subroutine input_data
       update_ocn_f = .false. ! include fresh water and salt fluxes for frazil
       ustar_min = 0.005      ! minimum friction velocity for ocean heat flux (m/s)
       iceruf = 0.0005_dbl_kind ! ice surface roughness at atmosphere interface (m)
+      iceruf_ocn = 0.03_dbl_kind ! under-ice roughness (m)
+      calc_dragio = .false.  ! compute dragio from iceruf_ocn and thickness of first ocean level
       emissivity = 0.985     ! emissivity of snow and ice
       l_mpond_fresh = .false.     ! logical switch for including meltpond freshwater
                                   ! flux feedback to ocean model
@@ -749,6 +751,8 @@ subroutine input_data
       call broadcast_scalar(l_mpond_fresh,        master_task)
       call broadcast_scalar(ustar_min,            master_task)
       call broadcast_scalar(iceruf,               master_task)
+      call broadcast_scalar(iceruf_ocn,           master_task)
+      call broadcast_scalar(calc_dragio,          master_task)
       call broadcast_scalar(emissivity,           master_task)
       call broadcast_scalar(fbot_xfer_type,       master_task)
       call broadcast_scalar(precip_units,         master_task)
@@ -1553,6 +1557,15 @@ subroutine input_data
          endif
          write(nu_diag,1030) ' fbot_xfer_type   = ', trim(fbot_xfer_type),trim(tmpstr2)
          write(nu_diag,1000) ' ustar_min        = ', ustar_min,' : minimum value of ocean friction velocity'
+         if (calc_dragio) then
+            tmpstr2 = ' : dragio computed from iceruf_ocn'
+         else
+            tmpstr2 = ' : dragio hard-coded'
+         endif
+         write(nu_diag,1010) ' calc_dragio   = ', calc_dragio,trim(tmpstr2)
+         if(calc_dragio) then
+            write(nu_diag,1002) ' iceruf_ocn       = ', iceruf_ocn,' : under-ice roughness length'
+         endif
 
          if (tr_fsd) then
             write(nu_diag,1002) ' floediam         = ', floediam, ' constant floe diameter'
@@ -1823,7 +1836,7 @@ subroutine input_data
          wave_spec_type_in = wave_spec_type, &
          wave_spec_in=wave_spec, nfreq_in=nfreq, &
          tfrz_option_in=tfrz_option, kalg_in=kalg, fbot_xfer_type_in=fbot_xfer_type, &
-         Pstar_in=Pstar, Cstar_in=Cstar, iceruf_in=iceruf, &
+         Pstar_in=Pstar, Cstar_in=Cstar, iceruf_in=iceruf, iceruf_ocn_in=iceruf_ocn, calc_dragio_in=calc_dragio, &
          sw_redist_in=sw_redist, sw_frac_in=sw_frac, sw_dtemp_in=sw_dtemp)
       call icepack_init_tracer_flags(tr_iage_in=tr_iage, tr_FY_in=tr_FY, &
          tr_lvl_in=tr_lvl, tr_iso_in=tr_iso, tr_aero_in=tr_aero, &

cicecore/cicedynB/infrastructure/ice_grid.F90

@@ -34,7 +34,7 @@ module ice_grid
       use ice_exit, only: abort_ice
       use ice_global_reductions, only: global_minval, global_maxval
       use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
-      use icepack_intfc, only: icepack_query_parameters
+      use icepack_intfc, only: icepack_query_parameters, icepack_init_parameters
 
       implicit none
       private
@@ -2371,6 +2371,9 @@ subroutine get_bathymetry
       real (kind=dbl_kind) :: &
          puny
 
+      logical (kind=log_kind) :: &
+         calc_dragio
+
       real (kind=dbl_kind), dimension(nlevel), parameter :: &
          thick  = (/ &                        ! ocean layer thickness, m
             10.01244_dbl_kind,  10.11258_dbl_kind,  10.31682_dbl_kind, &
@@ -2390,7 +2393,7 @@ subroutine get_bathymetry
 
       character(len=*), parameter :: subname = '(get_bathymetry)'
 
-      call icepack_query_parameters(puny_out=puny)
+      call icepack_query_parameters(puny_out=puny, calc_dragio_out=calc_dragio)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
          file=__FILE__, line=__LINE__)
@@ -2417,6 +2420,14 @@ subroutine get_bathymetry
             enddo
          enddo
 
+         ! For consistency, set thickness_ocn_layer1 in Icepack if 'calc_dragio' is active
+         if (calc_dragio) then
+            call icepack_init_parameters(thickness_ocn_layer1_in=thick(1))
+            call icepack_warnings_flush(nu_diag)
+            if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
+               file=__FILE__, line=__LINE__)
+         endif
+
       endif ! bathymetry_file
 
       end subroutine get_bathymetry
@@ -2440,6 +2451,9 @@ subroutine get_bathymetry_popfile
          depth          , & ! total depth, m
          thick              ! layer thickness, cm -> m
 
+      logical (kind=log_kind) :: &
+         calc_dragio
+
       character(len=*), parameter :: subname = '(get_bathymetry_popfile)'
 
       ntmp = maxval(nint(KMT))
@@ -2509,6 +2523,15 @@ subroutine get_bathymetry_popfile
          enddo
       enddo
 
+      ! For consistency, set thickness_ocn_layer1 in Icepack if 'calc_dragio' is active
+      call icepack_query_parameters(calc_dragio_out=calc_dragio)
+      if (calc_dragio) then
+         call icepack_init_parameters(thickness_ocn_layer1_in=thick(1))
+      endif
+      call icepack_warnings_flush(nu_diag)
+      if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
+         file=__FILE__, line=__LINE__)
+
       deallocate(depth,thick)
 
       end subroutine get_bathymetry_popfile

configuration/scripts/ice_in

@@ -207,6 +207,8 @@
     atmiter_conv    = 0.0d0
     ustar_min       = 0.0005
     iceruf          = 0.0005
+    calc_dragio     = .false.
+    iceruf_ocn      = 0.03
     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

@@ -68,3 +68,4 @@ restart        gx3     8x2        zsal
 restart        gx3     8x2        gx3ncarbulk,debug
 restart        gx3     4x4        gx3ncarbulk,diag1
 restart        gx1    24x1        gx1coreii,short
+smoke          gx3     4x1        calcdragio

doc/source/cice_index.rst

@@ -93,6 +93,7 @@ either Celsius or Kelvin units).
    "**C**", "", ""
    "c<n>", "real(\ :math:`n`)", ""
    "rotate_wind", ":math:`\bullet` if true, rotate wind/stress components to computational grid", "T"
+   "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", ""
@@ -322,6 +323,7 @@ either Celsius or Kelvin units).
    "ice_ref_salinity", "reference salinity for ice–ocean exchanges", "4. ppt"
    "icells", "number of grid cells with specified property (for vectorization)", ""
    "iceruf", ":math:`\bullet` ice surface roughness at atmosphere interface", "5.\ :math:`\times`\ 10\ :math:`^{-4}` m"
+   "iceruf_ocn", ":math:`\bullet` under-ice roughness (at ocean interface)", "0.03 m"
    "icetmask", "ice extent mask (T-cell)", "" 
    "iceumask", "ice extent mask (U-cell)", ""
    "idate", "the date at the end of the current time step (yyyymmdd)", ""