From 154065bced11ee47346f8173a4418a2d474c4751 Mon Sep 17 00:00:00 2001
From: Chloe <cwhicker@umich.edu>
Date: Thu, 13 Jun 2024 09:31:57 -0700
Subject: [PATCH] Added new snow Thermal Conductivity in SoilTemperatureMod.F90
 based on  Fourteau, et al. (2021) - is now density and temp dependent

---
 .../namelist_files/namelist_definition.xml    |  4 ++
 .../elm/src/biogeophys/SoilTemperatureMod.F90 | 61 ++++++++++++++++---
 components/elm/src/main/controlMod.F90        | 10 ++-
 components/elm/src/main/elm_varctl.F90        |  1 +
 4 files changed, 67 insertions(+), 9 deletions(-)

diff --git a/components/elm/bld/namelist_files/namelist_definition.xml b/components/elm/bld/namelist_files/namelist_definition.xml
index 16bca299ddd7..b111c1ff4282 100644
--- a/components/elm/bld/namelist_files/namelist_definition.xml
+++ b/components/elm/bld/namelist_files/namelist_definition.xml
@@ -728,6 +728,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_snow_thk" 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"
diff --git a/components/elm/src/biogeophys/SoilTemperatureMod.F90 b/components/elm/src/biogeophys/SoilTemperatureMod.F90
index fa3e822dd2fc..7d262e295a10 100644
--- a/components/elm/src/biogeophys/SoilTemperatureMod.F90
+++ b/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_snow_thk
     !
     ! !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,11 @@ 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)
+    data k_snw_tmps(:) /223.0, 248.0, 263.0, 268.0, 273.0/
     !-----------------------------------------------------------------------
     event = 'SoilThermProp'
     call t_start_lnd( event )
@@ -943,14 +948,54 @@ subroutine SoilThermProp (bounds,  num_nolakec, filter_nolakec, &
                   endif
                endif
             endif
+            
+            if (use_snow_thk) then ! chose 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))
+                    !write(iulog,*)"CAW bw(c,j)",bw(c,j)
+
+                       do i = 1, 5
+                        if (i == 1) then
+                            k_snw_vals(i) = 2.564 * (bw(c,j)/ rho_ice)**2 - 0.059 * (bw(c,j)/ rho_ice) + 0.0205
+                        else if (i == 2) then
+                            k_snw_vals(i) = 2.172 * (bw(c,j)/ rho_ice)**2 + 0.015 * (bw(c,j)/ rho_ice) + 0.0252
+                        else if (i == 3) then
+                            k_snw_vals(i) = 1.985 * (bw(c,j)/ rho_ice)**2 + 0.073 * (bw(c,j)/ rho_ice) + 0.0336
+                        else if (i == 4) then
+                            k_snw_vals(i) = 1.883 * (bw(c,j)/ rho_ice)**2 + 0.107 * (bw(c,j)/ rho_ice) + 0.0386
+                        else if (i == 5) then
+                            k_snw_vals(i) = 1.776 * (bw(c,j)/ rho_ice)**2 + 0.147 * (bw(c,j)/ rho_ice) + 0.0455
+                        end if
+                       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
+
+                     ! 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
+
+                     !  write(iulog,*)"CAW snow layer:",j
+                     !  write(iulog,*)"CAW snow temp:",t_soisno(c,j)
+                     !  write(iulog,*)"CAW snow density:",bw(c,j)
+                     !  write(iulog,*)"CAW snow thk:",thk(c,j) 
+               end if
 
-            ! 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
 
+            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
 
diff --git a/components/elm/src/main/controlMod.F90 b/components/elm/src/main/controlMod.F90
index 02c3ccf0867f..90f33c603ef4 100755
--- a/components/elm/src/main/controlMod.F90
+++ b/components/elm/src/main/controlMod.F90
@@ -282,7 +282,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_snow_thk, use_vancouver, use_mexicocity, use_noio
 
     ! cpl_bypass variables
     namelist /elm_inparm/ metdata_type, metdata_bypass, metdata_biases, &
@@ -723,7 +723,11 @@ subroutine control_spmd()
     call mpi_bcast (use_vertsoilc, 1, MPI_LOGICAL, 0, mpicom, ier)
     call mpi_bcast (use_extralakelayers, 1, MPI_LOGICAL, 0, mpicom, ier)
     call mpi_bcast (use_extrasnowlayers, 1, MPI_LOGICAL, 0, mpicom, ier)
+<<<<<<< HEAD
     call mpi_bcast (use_firn_percolation_and_compaction, 1, MPI_LOGICAL, 0, mpicom, ier)
+=======
+    call mpi_bcast (use_snow_thk, 1, MPI_LOGICAL, 0, mpicom, ier)
+>>>>>>> aa0f583905 (Added new snow Thermal Conductivity in SoilTemperatureMod.F90)
     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)
@@ -1024,7 +1028,11 @@ subroutine control_print ()
     write(iulog,*) '    use_lake_wat_storage = ', use_lake_wat_storage
     write(iulog,*) '    use_extralakelayers = ', use_extralakelayers
     write(iulog,*) '    use_extrasnowlayers = ', use_extrasnowlayers
+<<<<<<< HEAD
     write(iulog,*) '    use_firn_percolation_and_compaction = ', use_firn_percolation_and_compaction
+=======
+    write(iulog,*) '    use_snow_thk = ', use_snow_thk
+>>>>>>> aa0f583905 (Added new snow Thermal Conductivity in SoilTemperatureMod.F90)
     write(iulog,*) '    use_vichydro = ', use_vichydro
     write(iulog,*) '    use_century_decomp = ', use_century_decomp
     write(iulog,*) '    use_cn = ', use_cn
diff --git a/components/elm/src/main/elm_varctl.F90 b/components/elm/src/main/elm_varctl.F90
index b9ddd4b24397..29e05cb090e7 100644
--- a/components/elm/src/main/elm_varctl.F90
+++ b/components/elm/src/main/elm_varctl.F90
@@ -373,6 +373,7 @@ module elm_varctl
   logical, public :: use_snicar_ad       = .false.
   logical, public :: use_extrasnowlayers = .false.
   logical, public :: use_firn_percolation_and_compaction  = .false.
+  logical, public :: use_snow_thk        = .false.
   logical, public :: use_vancouver       = .false.
   logical, public :: use_mexicocity      = .false.
   logical, public :: use_noio            = .false.