Add sw_redist option to CICE (CICE-Consortium#497)

This integrates the sw_redist changes from Icepack into CICE.

Add three new namelist, sw_redist, sw_frac, and sw_dtemp.

The alt03 and alt04 test cases have been updated so that ktherm=1 is bfb. The default settings are such that ktherm=2 is bfb. This is the other piece to addressing issue CICE-Consortium#485 and CICE-Consortium/Icepack#280.

cicecore/cicedynB/general/ice_init.F90

@@ -118,15 +118,17 @@ subroutine input_data
         ahmax, R_ice, R_pnd, R_snw, dT_mlt, rsnw_mlt, emissivity, &
         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
+        phi_c_slow_mode, phi_i_mushy, kalg, atmiter_conv, Pstar, Cstar, &
+        sw_frac, sw_dtemp
 
       integer (kind=int_kind) :: ktherm, kstrength, krdg_partic, krdg_redist, natmiter, &
         kitd, kcatbound
 
       character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, &
         tfrz_option, frzpnd, atmbndy, wave_spec_type
 
-      logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair, wave_spec
+      logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair, wave_spec, &
+                                 sw_redist
 
       logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_pond
       logical (kind=log_kind) :: tr_iso, tr_aero, tr_fsd
@@ -182,7 +184,8 @@ subroutine input_data
       namelist /thermo_nml/ &
         kitd,           ktherm,          conduct,     ksno,             &
         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/ &
         kdyn,           ndte,           revised_evp,    yield_curve,    &
@@ -439,6 +442,11 @@ subroutine input_data
       phi_c_slow_mode   =    0.05_dbl_kind ! critical liquid fraction porosity cutoff
       phi_i_mushy       =    0.85_dbl_kind ! liquid fraction of congelation ice
 
+      ! shortwave redistribution in the thermodynamics
+      sw_redist = .false.
+      sw_frac   = 0.9_dbl_kind
+      sw_dtemp  = 0.02_dbl_kind
+
       !-----------------------------------------------------------------
       ! read from input file
       !-----------------------------------------------------------------
@@ -723,6 +731,9 @@ subroutine input_data
       call broadcast_scalar(dSdt_slow_mode,     master_task)
       call broadcast_scalar(phi_c_slow_mode,    master_task)
       call broadcast_scalar(phi_i_mushy,        master_task)
+      call broadcast_scalar(sw_redist,          master_task)
+      call broadcast_scalar(sw_frac,            master_task)
+      call broadcast_scalar(sw_dtemp,           master_task)
 
 #ifdef CESMCOUPLED
       pointer_file = trim(pointer_file) // trim(inst_suffix)
@@ -966,6 +977,12 @@ subroutine input_data
          endif
       endif
 !tcraig
+      if (ktherm == 1 .and. .not.sw_redist) then
+         if (my_task == master_task) then
+            write(nu_diag,*) subname//' WARNING: ktherm = 1 and sw_redist = ',sw_redist
+            write(nu_diag,*) subname//' WARNING:   For consistency, set sw_redist = .true.'
+         endif
+      endif
 
       if (formdrag) then
          if (trim(atmbndy) == 'constant') then
@@ -1244,6 +1261,9 @@ subroutine input_data
             write(nu_diag,1007) ' ksno             = ', ksno,' snow thermal conductivity'
             if (ktherm == 1) &
             write(nu_diag,*) 'conduct          = ', trim(conduct),' ice thermal conductivity'
+            write(nu_diag,1012) ' sw_redist        = ', sw_redist,' redistribute internal shortwave to surface'
+            write(nu_diag,1002) ' sw_frac          = ', sw_frac,' fraction redistributed'
+            write(nu_diag,1002) ' sw_dtemp         = ', sw_dtemp,' temperature difference from freezing to redistribute'
             if (ktherm == 2) then
                write(nu_diag,1002) ' a_rapid_mode     = ', a_rapid_mode,' brine channel diameter'
                write(nu_diag,1007) ' Rac_rapid_mode   = ', Rac_rapid_mode,' critical Rayleigh number'
@@ -1630,7 +1650,8 @@ 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)
+         Pstar_in=Pstar, Cstar_in=Cstar, &
+         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, &
          tr_fsd_in=tr_fsd, tr_pond_in=tr_pond, &

configuration/scripts/ice_in

@@ -108,6 +108,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
 /
 
 &dynamics_nml

configuration/scripts/options/set_nml.alt03

@@ -14,6 +14,9 @@ tr_aero      = .true.
 calc_Tsfc = .false.
 kdyn  =  2
 ktherm = 1
+sw_redist = .true.
+sw_frac = 0.9d0
+sw_dtemp = 0.02d0
 tfrz_option = 'linear_salt'
 revised_evp = .false.
 Ktens          = 0.

configuration/scripts/options/set_nml.alt04

@@ -13,6 +13,9 @@ tr_pond_lvl  = .true.
 tr_aero      = .true.
 kitd = 0
 ktherm = 1
+sw_redist = .true.
+sw_frac = 0.9d0
+sw_dtemp = 0.02d0
 conduct = 'MU71'
 kdyn  =  1
 kevp_kernel = 102

doc/source/user_guide/ug_case_settings.rst

@@ -344,6 +344,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", "redistribute internal shortwave to surface", "``.false.``"
+   "``sw_frac``", "real", "fraction redistributed", "0.9"
+   "``sw_dtemp``", "real", "temperature difference from melt to start redistributing", "0.02"
    "", "", "", ""
 
 dynamics_nml