Merge branch 'cwhicker/elm/use_use_T_rho_dependent_snowthk_20240719' …

…(PR #6524)

Implemented the new snow thermal conductivity calculation method introduced in
[Fourteau, et al. (2021)](https://tc.copernicus.org/articles/15/2739/2021/) in ELM.
This improves the current method, introduced in [Jordan (1991)](https://erdc-library.erdc.dren.mil/jspui/handle/11681/11677)
which only depends on the snow density.

relevant nl changes
`use_T_rho_dependent_snowthk = .true. `
[ELM]
[BFB]
[NML]

components/elm/bld/namelist_files/namelist_definition.xml

@@ -766,6 +766,10 @@ This option enables more realistic
 snowpack conditions on glaciers and ice sheets, including a semi-empirical
 firn densification model. uses the same physics as use_extrasnowlayers but
 uses the original 5 snow layer scheme
+
+<entry id="use_T_rho_dependent_snowthk" type="logical" category="physics"
+	group="elm_inparm" valid_values="" value=".false.">
+Toggle to use new snow thermal conductivity that relies on snow temperature and density. 
 </entry>
 
 <entry id="use_noio" type="logical" category="default_settings"

components/elm/src/biogeophys/SoilTemperatureMod.F90

@@ -834,7 +834,7 @@ subroutine SoilThermProp (bounds,  num_nolakec, filter_nolakec, &
     use elm_varcon      , only : denh2o, denice, tfrz, tkwat, tkice, tkair, cpice,  cpliq, thk_bedrock
     use landunit_varcon , only : istice, istice_mec, istwet
     use column_varcon   , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv, icol_road_imperv
-    use elm_varctl      , only : iulog
+    use elm_varctl      , only : iulog, use_T_rho_dependent_snowthk
     !
     ! !ARGUMENTS:
     type(bounds_type)      , intent(in)    :: bounds
@@ -847,7 +847,7 @@ subroutine SoilThermProp (bounds,  num_nolakec, filter_nolakec, &
     type(soilstate_type)   , intent(inout) :: soilstate_vars
     !
     ! !LOCAL VARIABLES:
-    integer  :: l,c,j                     ! indices
+    integer  :: l,c,j,i                     ! indices
     integer  :: nlevbed                   ! # levels to bedrock
     integer  :: fc                        ! lake filtered column indices
     real(r8) :: dksat                     ! thermal conductivity for saturated soil (j/(k s m))
@@ -856,6 +856,17 @@ subroutine SoilThermProp (bounds,  num_nolakec, filter_nolakec, &
     real(r8) :: satw                      ! relative total water content of soil.
     real(r8) :: zh2osfc
     character(len=64) :: event
+
+    real(r8), parameter :: rho_ice     = 917._r8
+    real(r8) :: k_snw_vals(5)
+    real(r8) :: k_snw_tmps(5)
+    real(r8) :: k_snw_coe1(5)
+    real(r8) :: k_snw_coe2(5)
+    real(r8) :: k_snw_coe3(5)
+    data k_snw_tmps(:) /223.0_r8, 248.0_r8, 263.0_r8, 268.0_r8, 273.0_r8/
+    data k_snw_coe1(:) /2.564_r8, 2.172_r8, 1.985_r8, 1.883_r8, 1.776_r8/
+    data k_snw_coe2(:) /-0.059_r8, 0.015_r8, 0.073_r8, 0.107_r8, 0.147_r8/
+    data k_snw_coe3(:) /0.0205_r8, 0.0252_r8, 0.0336_r8, 0.0386_r8, 0.0455_r8/
     !-----------------------------------------------------------------------
     event = 'SoilThermProp'
     call t_start_lnd( event )
@@ -943,14 +954,39 @@ subroutine SoilThermProp (bounds,  num_nolakec, filter_nolakec, &
                   endif
                endif
             endif
+
+            if (use_T_rho_dependent_snowthk) then ! choose which snow thermal conductivity to use 
+               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))
+
+                       do i = 1, 5
+                            k_snw_vals(i) = k_snw_coe1(i) * (bw(c,j) / rho_ice)**2 - k_snw_coe2(i) * (bw(c,j) / rho_ice) + k_snw_coe3(i)
+                       end do
+
+                       do i = 1, size(k_snw_tmps) - 1
+                        if (k_snw_tmps(i) <= t_soisno(c,j) .and. t_soisno(c,j) <= k_snw_tmps(i + 1)) then
+                            thk(c,j) = k_snw_vals(i) + (t_soisno(c,j) - k_snw_tmps(i)) * (k_snw_vals(i + 1)-k_snw_vals(i)) / (k_snw_tmps(i + 1) - k_snw_tmps(i))
+                        end if
+                       end do
 
-            ! Thermal conductivity of snow, which from Jordan (1991) pp. 18
-            ! 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)
-            end if
+                     ! Handle edge cases if t_soisno(c,j) is outside the given range
+                       if (t_soisno(c,j) < k_snw_tmps(1)) then
+                           thk(c,j) = k_snw_vals(1)
+                       else if (t_soisno(c,j) > k_snw_tmps(size(k_snw_tmps))) then
+                           thk(c,j) = k_snw_vals(size(k_snw_tmps))
+                       end if
 
+               end if
+
+
+            else 
+                    ! Thermal conductivity of snow, which from Jordan (1991) pp. 18
+                    ! 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)
+                    end if
+            endif
          end do
       end do
 

components/elm/src/main/controlMod.F90

@@ -285,7 +285,7 @@ subroutine control_init( )
          use_nofire, use_lch4, use_vertsoilc, use_extralakelayers, &
          use_vichydro, use_century_decomp, use_cn, use_crop, use_snicar_frc, &
          use_snicar_ad, use_firn_percolation_and_compaction, use_extrasnowlayers,&
-         use_vancouver, use_mexicocity, use_noio
+         use_T_rho_dependent_snowthk, use_vancouver, use_mexicocity, use_noio
 
     ! cpl_bypass variables
     namelist /elm_inparm/ metdata_type, metdata_bypass, metdata_biases, &
@@ -727,6 +727,7 @@ subroutine control_spmd()
     call mpi_bcast (use_extralakelayers, 1, MPI_LOGICAL, 0, mpicom, ier)
     call mpi_bcast (use_extrasnowlayers, 1, MPI_LOGICAL, 0, mpicom, ier)
     call mpi_bcast (use_firn_percolation_and_compaction, 1, MPI_LOGICAL, 0, mpicom, ier)
+    call mpi_bcast (use_T_rho_dependent_snowthk, 1, MPI_LOGICAL, 0, mpicom, ier)
     call mpi_bcast (use_vichydro, 1, MPI_LOGICAL, 0, mpicom, ier)
     call mpi_bcast (use_century_decomp, 1, MPI_LOGICAL, 0, mpicom, ier)
     call mpi_bcast (use_cn, 1, MPI_LOGICAL, 0, mpicom, ier)
@@ -1031,6 +1032,7 @@ subroutine control_print ()
     write(iulog,*) '    use_extralakelayers = ', use_extralakelayers
     write(iulog,*) '    use_extrasnowlayers = ', use_extrasnowlayers
     write(iulog,*) '    use_firn_percolation_and_compaction = ', use_firn_percolation_and_compaction
+    write(iulog,*) '    use_T_rho_dependent_snowthk = ', use_T_rho_dependent_snowthk
     write(iulog,*) '    use_vichydro = ', use_vichydro
     write(iulog,*) '    use_century_decomp = ', use_century_decomp
     write(iulog,*) '    use_cn = ', use_cn

components/elm/src/main/elm_varctl.F90

@@ -380,6 +380,7 @@ module elm_varctl
   logical, public :: use_mexicocity      = .false.
   logical, public :: use_noio            = .false.
   logical, public :: use_var_soil_thick  = .false.
+  logical, public :: use_T_rho_dependent_snowthk     = .false.
   logical, public :: use_atm_downscaling_to_topunit  = .false.
   character(len = SHR_KIND_CS), public :: precip_downscaling_method  = 'ERMM' ! Precip downscaling method values can be ERMM or FNM
   logical, public :: use_lake_wat_storage = .false.