Introduce namelist option for thermal conductivity of snow: Jordan vs. Sturm

bld/CLMBuildNamelist.pm

@@ -2031,6 +2031,13 @@ sub setup_logic_snicar_methods {
 sub setup_logic_snow {
   my ($opts, $nl_flags, $definition, $defaults, $nl) = @_;
 
+  add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'snow_thermal_cond_method' );
+
+  my $var = $nl->get_value('snow_thermal_cond_method');
+  if ( $var ne "'Jordan1991'" && $var ne "'Sturm1997'" ) {
+    $log->fatal_error("$var is incorrect entry for the namelist variable snow_thermal_cond_method; expected Jordan1991 or Sturm1997");
+  }
+
   my $numrad_snw = $nl->get_value('snicar_numrad_snw');
   add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'fsnowoptics',
                 'snicar_numrad_snw' => $numrad_snw);

bld/namelist_files/namelist_defaults_ctsm.xml

@@ -445,6 +445,7 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
 <reset_snow_glc_ela>1.e9</reset_snow_glc_ela>
 
 <snow_cover_fraction_method>SwensonLawrence2012</snow_cover_fraction_method>
+<snow_thermal_cond_method>Jordan1991</snow_thermal_cond_method>
 
 <!-- Default glacier behavior is:
      Mountain glaciers: single_at_atm_topo

bld/namelist_files/namelist_definition_ctsm.xml

@@ -2845,6 +2845,11 @@ NiuYang2007: Niu and Yang 2007
 SwensonLawrence2012: Swenson and Lawrence 2012
 </entry>
 
+<entry id="snow_thermal_cond_method" type="char*25" category="clm_physics"
+       group="clm_inparm" valid_values="Jordan1991,Sturm1997" >
+Parameterization to use for snow thermal conductivity
+</entry>
+
 <!-- ========================================================================================  -->
 <!-- Namelist options related to water tracers                                                 -->
 <!-- ========================================================================================  -->

src/biogeophys/SoilTemperatureMod.F90

@@ -47,7 +47,7 @@ module SoilTemperatureMod
   !     o The thermal conductivity of soil is computed from
   !       the algorithm of Johansen (as reported by Farouki 1981), and the
   !       conductivity of snow is from the formulation used in
-  !       SNTHERM (Jordan 1991).
+  !       Sturm (1997) or Jordan (1991) p. 18 depending on namelist option.
   !     o Boundary conditions:
   !       F = Rnet - Hg - LEg (top),  F= 0 (base of the soil column).
   !     o Soil / snow temperature is predicted from heat conduction
@@ -100,7 +100,7 @@ subroutine SoilTemperature(bounds, num_urbanl, filter_urbanl, num_urbanc, filter
     ! o The thermal conductivity of soil is computed from
     !   the algorithm of Johansen (as reported by Farouki 1981), and the
     !   conductivity of snow is from the formulation used in
-    !   SNTHERM (Jordan 1991).
+    !   Sturm (1997) or Jordan (1991) p. 18 depending on namelist option.
     ! o Boundary conditions:
     !   F = Rnet - Hg - LEg (top),  F= 0 (base of the soil column).
     ! o Soil / snow temperature is predicted from heat conduction
@@ -611,18 +611,20 @@ subroutine SoilThermProp (bounds, num_urbanc, filter_urbanc, num_nolakec, filter
     !
     ! (2) The thermal conductivity of soil is computed from the algorithm of
     !     Johansen (as reported by Farouki 1981), and of snow is from the
-    !     formulation used in SNTHERM (Jordan 1991).
+    !     formulation used in Sturm (1997) or Jordan (1991) p. 18 depending on
+    !     namelist option.
     ! The thermal conductivities at the interfaces between two neighboring
     ! layers (j, j+1) are derived from an assumption that the flux across
     ! the interface is equal to that from the node j to the interface and the
     ! flux from the interface to the node j+1.
     !
     ! !USES:
+    use shr_log_mod     , only : errMsg => shr_log_errMsg
     use clm_varpar      , only : nlevsno, nlevgrnd, nlevurb, nlevsoi, nlevmaxurbgrnd
     use clm_varcon      , only : denh2o, denice, tfrz, tkwat, tkice, tkair, cpice,  cpliq, thk_bedrock, csol_bedrock
     use landunit_varcon , only : istice, istwet
     use column_varcon   , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv, icol_road_imperv
-    use clm_varctl      , only : iulog
+    use clm_varctl      , only : iulog, snow_thermal_cond_method
     !
     ! !ARGUMENTS:
     type(bounds_type)      , intent(in)    :: bounds 
@@ -647,6 +649,8 @@ subroutine SoilThermProp (bounds, num_urbanc, filter_urbanc, num_nolakec, filter
     real(r8) :: fl                        ! volume fraction of liquid or unfrozen water to total water
     real(r8) :: satw                      ! relative total water content of soil.
     real(r8) :: zh2osfc
+
+    character(len=*),parameter :: subname = 'SoilThermProp'
     !-----------------------------------------------------------------------
 
     call t_startf( 'SoilThermProp' )
@@ -734,11 +738,27 @@ subroutine SoilThermProp (bounds, num_urbanc, filter_urbanc, num_nolakec, filter
                endif
             endif
 
-            ! Thermal conductivity of snow, which from Jordan (1991) pp. 18
+            ! Thermal conductivity of snow
             ! Only examine levels from snl(c)+1 -> 0 where snl(c) < 1
             if (snl(c)+1 < 1 .AND. (j >= snl(c)+1) .AND. (j <= 0)) then  
                bw(c,j) = (h2osoi_ice(c,j)+h2osoi_liq(c,j))/(frac_sno(c)*dz(c,j))
-               thk(c,j) = tkair + (7.75e-5_r8 *bw(c,j) + 1.105e-6_r8*bw(c,j)*bw(c,j))*(tkice-tkair)
+               select case (snow_thermal_cond_method)
+               case ('Jordan1991')
+                  thk(c,j) = tkair + (7.75e-5_r8 *bw(c,j) + 1.105e-6_r8*bw(c,j)*bw(c,j))*(tkice-tkair)
+               case ('Sturm1997')
+                  ! Implemented by Vicky Dutch (VRD), Nick Rutter, and
+                  ! Leanne Wake (LMW)
+                  ! https://tc.copernicus.org/articles/16/4201/2022/
+                  ! Code provided by Adrien Dams to Will Wieder
+                  if (bw(c,j) <= 156) then !LMW or 0.156 ?
+                     thk(c,j) = 0.023 + 0.234*(bw(c,j)/1000) !LMW - units changed by VRD
+                  else !LMW
+                     thk(c,j) = 0.138 - 1.01*(bw(c,j)/1000) +(3.233*((bw(c,j)/1000)*(bw(c,j)/1000))) ! LMW Sturm I think
+                  end if
+               case default
+                  write(iulog,*) subname//' ERROR: unknown snow_thermal_cond_method value: ', snow_thermal_cond_method
+                  call endrun(msg=errMsg(sourcefile, __LINE__))
+               end select
             end if
 
          end do

src/main/clm_varctl.F90

@@ -221,6 +221,8 @@ module clm_varctl
 
   ! which snow cover fraction parameterization to use
   character(len=64), public :: snow_cover_fraction_method
+  ! which snow thermal conductivity parameterization to use
+  character(len=25), public :: snow_thermal_cond_method
 
   ! atmospheric CO2 molar ratio (by volume) (umol/mol)
   real(r8), public :: co2_ppmv     = 355._r8            !

src/main/controlMod.F90

@@ -199,7 +199,8 @@ subroutine control_init(dtime)
          clump_pproc, &
          create_crop_landunit, nsegspc, co2_ppmv, &
          albice, soil_layerstruct_predefined, soil_layerstruct_userdefined, &
-         soil_layerstruct_userdefined_nlevsoi, use_subgrid_fluxes, snow_cover_fraction_method, &
+         soil_layerstruct_userdefined_nlevsoi, use_subgrid_fluxes, &
+         snow_thermal_cond_method, snow_cover_fraction_method, &
          irrigate, run_zero_weight_urban, all_active, &
          crop_fsat_equals_zero, for_testing_run_ncdiopio_tests, &
          for_testing_use_second_grain_pool, for_testing_use_repr_structure_pool, &
@@ -850,6 +851,7 @@ subroutine control_spmd()
     ! physics variables
     call mpi_bcast (nsegspc, 1, MPI_INTEGER, 0, mpicom, ier)
     call mpi_bcast (use_subgrid_fluxes , 1, MPI_LOGICAL, 0, mpicom, ier)
+    call mpi_bcast (snow_thermal_cond_method, len(snow_thermal_cond_method), MPI_CHARACTER, 0, mpicom, ier)
     call mpi_bcast (snow_cover_fraction_method , len(snow_cover_fraction_method), MPI_CHARACTER, 0, mpicom, ier)
     call mpi_bcast (z0param_method , len(z0param_method), MPI_CHARACTER, 0, mpicom, ier)
     call mpi_bcast (use_z0m_snowmelt, 1, MPI_LOGICAL, 0, mpicom, ier)