diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index ca21eeec1..ad571ab87 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/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 !======================================================================= diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90 index 3146c0a9a..781f49b22 100644 --- a/configuration/driver/icedrv_init.F90 +++ b/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, & diff --git a/configuration/scripts/icepack_in b/configuration/scripts/icepack_in index 0b8ba3d47..04f8aa408 100644 --- a/configuration/scripts/icepack_in +++ b/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. diff --git a/configuration/scripts/options/set_nml.calcdragio b/configuration/scripts/options/set_nml.calcdragio new file mode 100644 index 000000000..cf86664bf --- /dev/null +++ b/configuration/scripts/options/set_nml.calcdragio @@ -0,0 +1 @@ +calc_dragio = .true. diff --git a/configuration/scripts/tests/base_suite.ts b/configuration/scripts/tests/base_suite.ts index 194ffbf33..3f1c4572a 100644 --- a/configuration/scripts/tests/base_suite.ts +++ b/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 diff --git a/doc/source/icepack_index.rst b/doc/source/icepack_index.rst index e1c1816fd..601c54310 100755 --- a/doc/source/icepack_index.rst +++ b/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" diff --git a/doc/source/master_list.bib b/doc/source/master_list.bib index 1a1e1b73e..6aba6a416 100644 --- a/doc/source/master_list.bib +++ b/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}", diff --git a/doc/source/science_guide/sg_boundary_forcing.rst b/doc/source/science_guide/sg_boundary_forcing.rst index 763a757d0..e3006ed5f 100755 --- a/doc/source/science_guide/sg_boundary_forcing.rst +++ b/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 diff --git a/doc/source/user_guide/interfaces.include b/doc/source/user_guide/interfaces.include index bb9f65556..907959bf4 100644 --- a/doc/source/user_guide/interfaces.include +++ b/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)