Shortwave redistribution option added to namelist. (#326)

* Fix nt_zbgc_frac

* Add sw_redist namelist option

* Fix typo

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
 

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.

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
 

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)

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

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.

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

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'

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", ""

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

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