From 76764c93cac9a5737ddd02d2fe8b1596ec05dec5 Mon Sep 17 00:00:00 2001 From: "David A. Bailey" Date: Wed, 15 Jul 2020 12:03:08 -0600 Subject: [PATCH] Shortwave redistribution option added to namelist. (#326) * Fix nt_zbgc_frac * Add sw_redist namelist option * Fix typo --- columnphysics/icepack_parameters.F90 | 41 ++++++++++++++++++- columnphysics/icepack_therm_bl99.F90 | 31 ++++++-------- columnphysics/icepack_therm_mushy.F90 | 41 +++++++++++++++++++ configuration/driver/icedrv_init.F90 | 15 +++++-- configuration/scripts/icepack_in | 3 ++ configuration/scripts/options/set_nml.alt03 | 3 ++ configuration/scripts/options/set_nml.alt04 | 3 ++ configuration/scripts/options/set_nml.thermo1 | 3 ++ doc/source/icepack_index.rst | 11 +++-- doc/source/science_guide/sg_thermo.rst | 6 +++ doc/source/user_guide/ug_case_settings.rst | 3 ++ 11 files changed, 133 insertions(+), 27 deletions(-) diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index 1e346223c..b3c0f1347 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -355,6 +355,18 @@ module icepack_parameters t_sk_conv = 3.0_dbl_kind , & ! Stefels conversion time (d) t_sk_ox = 10.0_dbl_kind ! DMS oxidation time (d) + +!----------------------------------------------------------------------- +! Parameters for shortwave redistribution +!----------------------------------------------------------------------- + + logical (kind=log_kind), public :: & + sw_redist = .false. + + real (kind=dbl_kind), public :: & + sw_frac = 0.9_dbl_kind , & ! Fraction of internal shortwave moved to surface + sw_dtemp = 0.02_dbl_kind ! temperature difference from melting + !======================================================================= contains @@ -400,7 +412,8 @@ subroutine icepack_init_parameters( & op_dep_min_in, fr_graze_s_in, fr_graze_e_in, fr_mort2min_in, & fr_dFe_in, k_nitrif_in, t_iron_conv_in, max_loss_in, & max_dfe_doc1_in, fr_resp_s_in, conserv_check_in, & - y_sk_DMS_in, t_sk_conv_in, t_sk_ox_in, frazil_scav_in) + y_sk_DMS_in, t_sk_conv_in, t_sk_ox_in, frazil_scav_in, & + sw_redist_in, sw_frac_in, sw_dtemp_in) !----------------------------------------------------------------- ! parameter constants @@ -527,6 +540,13 @@ subroutine icepack_init_parameters( & rsnw_mlt_in , & ! maximum melting snow grain radius (10^-6 m) kalg_in ! algae absorption coefficient for 0.5 m thick layer + logical (kind=log_kind), intent(in), optional :: & + sw_redist_in ! redistribute shortwave + + real (kind=dbl_kind), intent(in), optional :: & + sw_frac_in , & ! Fraction of internal shortwave moved to surface + sw_dtemp_in ! temperature difference from melting + !----------------------------------------------------------------------- ! Parameters for dynamics !----------------------------------------------------------------------- @@ -835,6 +855,9 @@ subroutine icepack_init_parameters( & if (present(t_sk_conv_in) ) t_sk_conv = t_sk_conv_in if (present(t_sk_ox_in) ) t_sk_ox = t_sk_ox_in if (present(frazil_scav_in) ) frazil_scav = frazil_scav_in + if (present(sw_redist_in) ) sw_redist = sw_redist_in + if (present(sw_frac_in) ) sw_frac = sw_frac_in + if (present(sw_dtemp_in) ) sw_dtemp = sw_dtemp_in call icepack_recompute_constants() if (icepack_warnings_aborted(subname)) return @@ -887,7 +910,8 @@ subroutine icepack_query_parameters( & T_max_out, fsal_out, op_dep_min_out, fr_graze_s_out, fr_graze_e_out, & fr_mort2min_out, fr_resp_s_out, fr_dFe_out, & k_nitrif_out, t_iron_conv_out, max_loss_out, max_dfe_doc1_out, & - y_sk_DMS_out, t_sk_conv_out, t_sk_ox_out, frazil_scav_out) + y_sk_DMS_out, t_sk_conv_out, t_sk_ox_out, frazil_scav_out, & + sw_redist_out, sw_frac_out, sw_dtemp_out) !----------------------------------------------------------------- ! parameter constants @@ -1023,6 +1047,13 @@ subroutine icepack_query_parameters( & rsnw_mlt_out , & ! maximum melting snow grain radius (10^-6 m) kalg_out ! algae absorption coefficient for 0.5 m thick layer + logical (kind=log_kind), intent(out), optional :: & + sw_redist_out ! redistribute shortwave + + real (kind=dbl_kind), intent(out), optional :: & + sw_frac_out , & ! Fraction of internal shortwave moved to surface + sw_dtemp_out ! temperature difference from melting + !----------------------------------------------------------------------- ! Parameters for dynamics !----------------------------------------------------------------------- @@ -1375,6 +1406,9 @@ subroutine icepack_query_parameters( & if (present(Lfresh_out) ) Lfresh_out = Lfresh if (present(cprho_out) ) cprho_out = cprho if (present(Cp_out) ) Cp_out = Cp + if (present(sw_redist_out) ) sw_redist_out = sw_redist + if (present(sw_frac_out) ) sw_frac_out = sw_frac + if (present(sw_dtemp_out) ) sw_dtemp_out = sw_dtemp call icepack_recompute_constants() if (icepack_warnings_aborted(subname)) return @@ -1547,6 +1581,9 @@ subroutine icepack_write_parameters(iounit) write(iounit,*) " t_sk_conv = ", t_sk_conv write(iounit,*) " t_sk_ox = ", t_sk_ox write(iounit,*) " frazil_scav = ", frazil_scav + write(iounit,*) " sw_redist = ", sw_redist + write(iounit,*) " sw_frac = ", sw_frac + write(iounit,*) " sw_dtemp = ", sw_dtemp end subroutine icepack_write_parameters diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index e5483cb18..1c3a3aabb 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -18,6 +18,7 @@ module icepack_therm_bl99 #endif use icepack_parameters, only: rhoi, rhos, hs_min, cp_ice, cp_ocn, depressT, Lfresh, ksno, kice use icepack_parameters, only: conduct, calc_Tsfc, solve_zsal + use icepack_parameters, only: sw_redist, sw_frac, sw_dtemp use icepack_warnings, only: warnstr, icepack_warnings_add use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted @@ -201,8 +202,7 @@ subroutine temperature_changes (dt, & Iswabs_tmp , & ! energy to melt through fraction frac of layer Sswabs_tmp , & ! same for snow dswabs , & ! difference in swabs and swabs_tmp - frac , & ! fraction of layer that can be melted through - dTemp ! minimum temperature difference for absorption + frac logical (kind=log_kind) :: & converged ! = true when local solution has converged @@ -272,26 +272,22 @@ subroutine temperature_changes (dt, & !----------------------------------------------------------------- !mclaren: Should there be an if calc_Tsfc statement here then?? -#ifdef CESMCOUPLED - frac = c1 - dTemp = p01 -#else - frac = 0.9_dbl_kind - dTemp = 0.02_dbl_kind -#endif - if (solve_zsal) dTemp = p1 ! lower tolerance with dynamic salinity + if (sw_redist) then + + if (solve_zsal) sw_dtemp = p1 ! lower tolerance with dynamic salinity + do k = 1, nilyr Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc - if (Tin_init(k) <= Tmlts(k) - dTemp) then + if (Tin_init(k) <= Tmlts(k) - sw_dtemp) then if (l_brine) then ci = cp_ice - Lfresh * Tmlts(k) / (Tin_init(k)**2) Iswabs_tmp = min(Iswabs(k), & - frac*(Tmlts(k)-Tin_init(k))*ci/dt_rhoi_hlyr) + sw_frac*(Tmlts(k)-Tin_init(k))*ci/dt_rhoi_hlyr) else ci = cp_ice Iswabs_tmp = min(Iswabs(k), & - frac*(-Tin_init(k))*ci/dt_rhoi_hlyr) + sw_frac*(-Tin_init(k))*ci/dt_rhoi_hlyr) endif endif if (Iswabs_tmp < puny) Iswabs_tmp = c0 @@ -304,16 +300,13 @@ subroutine temperature_changes (dt, & enddo -#ifdef CESMCOUPLED - frac = 0.9_dbl_kind -#endif do k = 1, nslyr if (l_snow) then Sswabs_tmp = c0 - if (Tsn_init(k) <= -dTemp) then + if (Tsn_init(k) <= -sw_dtemp) then Sswabs_tmp = min(Sswabs(k), & - -frac*Tsn_init(k)/etas(k)) + -sw_frac*Tsn_init(k)/etas(k)) endif if (Sswabs_tmp < puny) Sswabs_tmp = c0 @@ -326,6 +319,8 @@ subroutine temperature_changes (dt, & endif enddo + endif + !----------------------------------------------------------------- ! Solve for new temperatures. ! Iterate until temperatures converge with minimal energy error. diff --git a/columnphysics/icepack_therm_mushy.F90 b/columnphysics/icepack_therm_mushy.F90 index 45418a727..b379e20ac 100644 --- a/columnphysics/icepack_therm_mushy.F90 +++ b/columnphysics/icepack_therm_mushy.F90 @@ -9,6 +9,7 @@ module icepack_therm_mushy use icepack_parameters, only: hs_min use icepack_parameters, only: a_rapid_mode, Rac_rapid_mode use icepack_parameters, only: aspect_rapid_mode, dSdt_slow_mode, phi_c_slow_mode + use icepack_parameters, only: sw_redist, sw_frac, sw_dtemp use icepack_mushy_physics, only: icepack_mushy_density_brine, enthalpy_brine, enthalpy_snow use icepack_mushy_physics, only: enthalpy_mush_liquid_fraction use icepack_mushy_physics, only: icepack_mushy_temperature_mush, icepack_mushy_liquid_fraction @@ -141,6 +142,8 @@ subroutine temperature_changes_salinity(dt, & qpond , & ! melt pond brine enthalpy (J m-3) Spond ! melt pond salinity (ppt) + real(kind=dbl_kind) :: Tmlt, Iswabs_tmp, dt_rhoi_hlyr, ci, dswabs + integer(kind=int_kind) :: & k ! ice/snow layer index @@ -201,6 +204,44 @@ subroutine temperature_changes_salinity(dt, & call conductivity_mush_array(nilyr, zqin0, zSin0, km) if (icepack_warnings_aborted(subname)) return + !----------------------------------------------------------------- + ! Check for excessive absorbed solar radiation that may result in + ! temperature overshoots. Convergence is particularly difficult + ! if the starting temperature is already very close to the melting + ! temperature and extra energy is added. In that case, or if the + ! amount of energy absorbed is greater than the amount needed to + ! melt through a given fraction of a layer, we put the extra + ! energy into the surface. + ! NOTE: This option is not available if the atmosphere model + ! has already computed fsurf. (Unless we adjust fsurf here) + !----------------------------------------------------------------- +!mclaren: Should there be an if calc_Tsfc statement here then?? + + if (sw_redist) then + + do k = 1, nilyr + + Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc + + Tmlt = liquidus_temperature_mush(zSin(k)) + + if (zTin(k) <= Tmlt - sw_dtemp) then + ci = cp_ice - Lfresh * Tmlt / (zTin(k)**2) + Iswabs_tmp = min(Iswabs(k), & + sw_frac*(Tmlt-zTin(k))*ci/dt_rhoi_hlyr) + endif + if (Iswabs_tmp < puny) Iswabs_tmp = c0 + + dswabs = min(Iswabs(k) - Iswabs_tmp, fswint) + + fswsfc = fswsfc + dswabs + fswint = fswint - dswabs + Iswabs(k) = Iswabs_tmp + + enddo + + endif + if (lsnow) then ! case with snow diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90 index 221dfeee4..3fecc57cc 100644 --- a/configuration/driver/icedrv_init.F90 +++ b/configuration/driver/icedrv_init.F90 @@ -96,6 +96,9 @@ subroutine input_data character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, & tfrz_option, frzpnd, atmbndy, wave_spec_type + logical (kind=log_kind) :: sw_redist + real (kind=dbl_kind) :: sw_frac, sw_dtemp + ! Flux convergence tolerance real (kind=dbl_kind) :: atmiter_conv @@ -131,7 +134,8 @@ subroutine input_data namelist /thermo_nml/ & kitd, ktherm, ksno, conduct, & a_rapid_mode, Rac_rapid_mode, aspect_rapid_mode, & - dSdt_slow_mode, phi_c_slow_mode, phi_i_mushy + dSdt_slow_mode, phi_c_slow_mode, phi_i_mushy, & + sw_redist, sw_frac, sw_dtemp namelist /dynamics_nml/ & kstrength, krdg_partic, krdg_redist, mu_rdg, & @@ -203,7 +207,8 @@ subroutine input_data phi_i_mushy_out=phi_i_mushy, conserv_check_out=conserv_check, & tfrz_option_out=tfrz_option, kalg_out=kalg, & fbot_xfer_type_out=fbot_xfer_type, puny_out=puny, & - wave_spec_type_out=wave_spec_type) + wave_spec_type_out=wave_spec_type, & + sw_redist_out=sw_redist, sw_frac_out=sw_frac, sw_dtemp_out=sw_dtemp) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__, line=__LINE__) @@ -559,6 +564,9 @@ subroutine input_data write(nu_diag,1005) ' phi_c_slow_mode = ', phi_c_slow_mode write(nu_diag,1005) ' phi_i_mushy = ', phi_i_mushy endif + write(nu_diag,1010) ' sw_redist = ', sw_redist + write(nu_diag,1005) ' sw_frac = ', sw_frac + write(nu_diag,1005) ' sw_dtemp = ', sw_dtemp write(nu_diag,1030) ' atmbndy = ', & trim(atmbndy) @@ -766,7 +774,8 @@ subroutine input_data phi_i_mushy_in=phi_i_mushy, conserv_check_in=conserv_check, & tfrz_option_in=tfrz_option, kalg_in=kalg, & fbot_xfer_type_in=fbot_xfer_type, & - wave_spec_type_in=wave_spec_type, wave_spec_in=wave_spec) + wave_spec_type_in=wave_spec_type, wave_spec_in=wave_spec, & + sw_redist_in=sw_redist, sw_frac_in=sw_frac, sw_dtemp_in=sw_dtemp) call icepack_init_tracer_sizes(ntrcr_in=ntrcr, & ncat_in=ncat, nilyr_in=nilyr, nslyr_in=nslyr, nblyr_in=nblyr, & nfsd_in=nfsd, n_iso_in=n_iso, n_aero_in=n_aero) diff --git a/configuration/scripts/icepack_in b/configuration/scripts/icepack_in index c4ec7fabf..553ee26e4 100644 --- a/configuration/scripts/icepack_in +++ b/configuration/scripts/icepack_in @@ -44,6 +44,9 @@ dSdt_slow_mode = -5.0e-8 phi_c_slow_mode = 0.05 phi_i_mushy = 0.85 + sw_redist = .false. + sw_frac = 0.9d0 + sw_dtemp = 0.02d0 / &shortwave_nml diff --git a/configuration/scripts/options/set_nml.alt03 b/configuration/scripts/options/set_nml.alt03 index b8cbaa370..a03f0d90a 100644 --- a/configuration/scripts/options/set_nml.alt03 +++ b/configuration/scripts/options/set_nml.alt03 @@ -3,6 +3,9 @@ tr_pond_topo = .true. tr_pond_cesm = .false. ktherm = 1 + sw_redist = .true. + sw_frac = 0.9d0 + sw_dtemp = 0.02d0 tfrz_option = 'linear_salt' conduct = 'bubbly' conserv_check = .true. diff --git a/configuration/scripts/options/set_nml.alt04 b/configuration/scripts/options/set_nml.alt04 index 0d83553ce..a112b235c 100644 --- a/configuration/scripts/options/set_nml.alt04 +++ b/configuration/scripts/options/set_nml.alt04 @@ -6,6 +6,9 @@ tr_pond_cesm = .false. frzpnd = 'cesm' ktherm = 1 + sw_redist = .true. + sw_frac = 0.9d0 + sw_dtemp = 0.02d0 tfrz_option = 'linear_salt' conduct = 'bubbly' krdg_partic = 0 diff --git a/configuration/scripts/options/set_nml.thermo1 b/configuration/scripts/options/set_nml.thermo1 index 69002ee19..39d204492 100644 --- a/configuration/scripts/options/set_nml.thermo1 +++ b/configuration/scripts/options/set_nml.thermo1 @@ -1,6 +1,9 @@ kcatbound = 0 kitd = 0 ktherm = 1 +sw_redist = .true. +sw_frac = 0.9d0 +sw_dtemp = 0.02d0 conduct = 'MU71' tfrz_option = 'linear_salt' atm_data_type = 'clim' diff --git a/doc/source/icepack_index.rst b/doc/source/icepack_index.rst index da4dff188..a88c0b363 100755 --- a/doc/source/icepack_index.rst +++ b/doc/source/icepack_index.rst @@ -185,10 +185,10 @@ either Celsius or Kelvin units). "fswint", "shortwave absorbed in ice interior", "W/m\ :math:`^2`" "fswpenl", "shortwave penetrating through ice layers", "W/m\ :math:`^2`" "fswthru", "shortwave penetrating to ocean", "W/m\ :math:`^2`" - "fswthru_vdr", "vis dir shortwave penetrating to ocean", "W/m\ :math:`^2`" - "fswthru_vdf", "vis dif shortwave penetrating to ocean", "W/m\ :math:`^2`" - "fswthru_idr", "nir dir shortwave penetrating to ocean", "W/m\ :math:`^2`" - "fswthru_idf", "nir dif shortwave penetrating to ocean", "W/m\ :math:`^2`" + "fswthru_vdr", "visible direct shortwave penetrating to ocean", "W/m\ :math:`^2`" + "fswthru_vdf", "visible diffuse shortwave penetrating to ocean", "W/m\ :math:`^2`" + "fswthru_idr", "near IR direct shortwave penetrating to ocean", "W/m\ :math:`^2`" + "fswthru_idf", "near IR diffuse shortwave penetrating to ocean", "W/m\ :math:`^2`" "fswthru_ai", "grid-box-mean shortwave penetrating to ocean (fswthru)", "W/m\ :math:`^2`" "fyear", "current data year", "" "fyear_final", "last data year", "" @@ -437,6 +437,9 @@ either Celsius or Kelvin units). "swndr", "incoming shortwave radiation, near IR, direct", "W/m\ :math:`^2`" "swvdf", "incoming shortwave radiation, visible, diffuse", "W/m\ :math:`^2`" "swvdr", "incoming shortwave radiation, visible, direct", "W/m\ :math:`^2`" + "sw_redist", "option to redistribute shortwave", ".false." + "sw_frac", "fraction of redistributed shortwave", "0.9" + "sw_dtemp", "temperature threshold from melting for redistributed shortwave", "0.02" "**T**", "", "" "Tair", "air temperature at 10 m", "K" "tday", "absolute day number", "" diff --git a/doc/source/science_guide/sg_thermo.rst b/doc/source/science_guide/sg_thermo.rst index 8bd14320c..fdffa7ac6 100755 --- a/doc/source/science_guide/sg_thermo.rst +++ b/doc/source/science_guide/sg_thermo.rst @@ -50,6 +50,12 @@ the :ref:`thermo-growth` section. We begin by describing the melt ponds and sur forcing parameterizations, which are closely related to the ice and snow surface temperatures. +Sometimes instabilities can arise when the temperature is close to the melt point in +the snow and sea ice and there is abundant internal shorwave absorbed. One can choose +to "move" the excess internal shortwave in this case up to the top surface to be reabsorbed. +The namelist parameters for this option are ``sw_redist``, ``sw_frac``, and ``sw_dtemp``. +By default, ``sw_redist`` is set to ``.false.`` + .. _ponds: Melt ponds diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 6e67c67f0..95fab3f47 100755 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -182,6 +182,9 @@ thermo_nml "``phi_c_slow_mode``", ":math:`0<\phi_c < 1`", "critical liquid fraction", "0.05" "``phi_i_mushy``", ":math:`0<\phi_i < 1`", "solid fraction at lower boundary", "0.85" "``Rac_rapid_mode``", "real", "critical Rayleigh number", "10.0" + "``sw_redist``", "logical", "shortwave redistribution", ".false." + "``sw_frac``", "real", "fraction of shortwave redistribution", "0.9" + "``sw_dtemp``", "real", "temperature from melt for sw_redist", "0.02" "", "", "", "" dynamics_nml