User control over snow thermal conductivity scheme over glaciers

bld/CLMBuildNamelist.pm

@@ -2196,6 +2196,16 @@ sub setup_logic_snow {
     $log->fatal_error("$var is incorrect entry for the namelist variable snow_thermal_cond_method; expected Jordan1991 or Sturm1997");
   }
 
+  my $var = $nl->get_value('snow_thermal_cond_glc_method');
+  if ( $var ne "'Jordan1991'" && $var ne "'Sturm1997'" ) {
+    $log->fatal_error("$var is incorrect entry for the namelist variable snow_thermal_cond_glc_method; expected Jordan1991 or Sturm1997");
+  }
+
+  my $var = $nl->get_value('snow_thermal_cond_lake_method');
+  if ( $var ne "'Jordan1991'" && $var ne "'Sturm1997'" ) {
+    $log->fatal_error("$var is incorrect entry for the namelist variable snow_thermal_cond_lake_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

@@ -484,6 +484,14 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
 <snow_thermal_cond_method>Jordan1991</snow_thermal_cond_method>
 <snow_thermal_cond_method phys="clm6_0">Sturm1997</snow_thermal_cond_method>
 
+<!-- TODO change clm6 to Jordan with answer changing tag for glacier mass ballance -->
+<snow_thermal_cond_glc_method>Jordan1991</snow_thermal_cond_glc_method>
+<snow_thermal_cond_glc_method phys="clm6_0">Sturm1997</snow_thermal_cond_glc_method>
+
+<!-- TODO change clm6 to Sturm with answer changing tag -->
+<snow_thermal_cond_lake_method>Jordan1991</snow_thermal_cond_lake_method>
+<snow_thermal_cond_lake_method phys="clm6_0">Jordan1991</snow_thermal_cond_lake_method>
+
 <!-- Default glacier behavior is:
      Mountain glaciers: single_at_atm_topo
      Currently-unused region 1: UNSET

bld/namelist_files/namelist_definition_ctsm.xml

@@ -2960,6 +2960,16 @@ SwensonLawrence2012: Swenson and Lawrence 2012
 Parameterization to use for snow thermal conductivity
 </entry>
 
+<entry id="snow_thermal_cond_glc_method" type="char*25" category="clm_physics"
+       group="clm_inparm" valid_values="Jordan1991,Sturm1997" >
+Parameterization to use for snow thermal conductivity over glacier land units
+</entry>
+
+<entry id="snow_thermal_cond_lake_method" type="char*25" category="clm_physics"
+       group="clm_inparm" valid_values="Jordan1991,Sturm1997" >
+Parameterization to use for snow thermal conductivity over lake land units
+</entry>
+
 <!-- ========================================================================================  -->
 <!-- Namelist options related to water tracers                                                 -->
 <!-- ========================================================================================  -->

src/biogeophys/LakeFluxesMod.F90

@@ -96,6 +96,7 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep,
     use clm_varcon          , only : beta_param, nu_param, b1_param, b4_param
     use clm_varcon          , only : meier_param1, meier_param2, meier_param3
     use clm_varctl          , only : use_lch4, z0param_method, use_z0m_snowmelt
+    use clm_varctl          , only : snow_thermal_cond_lake_method
     use LakeCon             , only : betavis, z0frzlake, tdmax, emg_lake
     use LakeCon             , only : lake_use_old_fcrit_minz0
     use LakeCon             , only : minz0lake, cur0, cus, curm, fcrit
@@ -480,7 +481,19 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep,
             else
                !Need to calculate thermal conductivity of the top snow layer
                bw = (h2osoi_ice(c,jtop(c))+h2osoi_liq(c,jtop(c)))/dz(c,jtop(c))
-               tksur(c) = tkair + (7.75e-5_r8 *bw + 1.105e-6_r8*bw*bw)*(tkice-tkair)
+               select case (snow_thermal_cond_lake_method)
+               case ('Jordan1991')
+                  tksur(c) = tkair + (7.75e-5_r8 *bw + 1.105e-6_r8*bw*bw)*(tkice-tkair)
+               case ('Sturm1997')
+                  if (bw <= 156) then
+                     tksur(c) = 0.023 + 0.234*(bw/1000)
+                  else
+                     tksur(c) = 0.138 - 1.01*(bw/1000) +(3.233*((bw/1000)*(bw/1000)))
+                  end if
+               case default
+
+               end select
+
                tsur(c) = t_soisno(c,jtop(c))
             end if
 

src/biogeophys/LakeTemperatureMod.F90

@@ -1085,6 +1085,7 @@ subroutine SoilThermProp_Lake (bounds,  num_lakec, filter_lakec, tk, cv, tktopso
      use clm_varcon  , only : denh2o, denice, tfrz, tkwat, tkice, tkair
      use clm_varcon  , only : cpice,  cpliq, thk_bedrock, csol_bedrock
      use clm_varpar  , only : nlevsno, nlevsoi, nlevgrnd
+     use clm_varctl  , only : snow_thermal_cond_lake_method
      !
      ! !ARGUMENTS:
      type(bounds_type)      , intent(in)  :: bounds  
@@ -1181,7 +1182,18 @@ subroutine SoilThermProp_Lake (bounds,  num_lakec, filter_lakec, tk, cv, tktopso
              ! 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 = (h2osoi_ice(c,j)+h2osoi_liq(c,j))/dz(c,j)
-                thk(c,j) = tkair + (7.75e-5_r8 *bw + 1.105e-6_r8*bw*bw)*(tkice-tkair)
+                select case (snow_thermal_cond_lake_method)
+                case ('Jordan1991')
+                   thk(c,j) = tkair + (7.75e-5_r8 *bw + 1.105e-6_r8*bw*bw)*(tkice-tkair)
+                case ('Sturm1997')
+                   if (bw <= 156) then 
+                      thk(c,j) = 0.023 + 0.234*(bw/1000) 
+                   else
+                      thk(c,j) = 0.138 - 1.01*(bw/1000) +(3.233*((bw/1000)*(bw/1000))) 
+                   end if
+                case default
+                end select
+
              end if
 
           end do

src/biogeophys/SoilTemperatureMod.F90

@@ -624,7 +624,7 @@ subroutine SoilThermProp (bounds, num_urbanc, filter_urbanc, num_nolakec, filter
     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, snow_thermal_cond_method
+    use clm_varctl      , only : iulog, snow_thermal_cond_method, snow_thermal_cond_glc_method
     !
     ! !ARGUMENTS:
     type(bounds_type)      , intent(in)    :: bounds 
@@ -742,24 +742,46 @@ subroutine SoilThermProp (bounds, num_urbanc, filter_urbanc, num_nolakec, filter
             ! 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))
-               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
+               l = col%landunit(c)
+
+               ! Select method over glacier land unit 
+               if (lun%itype(l) == istice) then
+                  select case (snow_thermal_cond_glc_method)
+                  ! TODO, this code duplication isn't ideal and should likely be in it's own subroutine
+                  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')
+                     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 
+                        thk(c,j) = 0.138 - 1.01*(bw(c,j)/1000) +(3.233*((bw(c,j)/1000)*(bw(c,j)/1000))) 
+                     end if
+                  case default
+                     write(iulog,*) subname//' ERROR: unknown snow_thermal_cond_glc_method value: ', snow_thermal_cond_glc_method
+                     call endrun(msg=errMsg(sourcefile, __LINE__))
+                  end select
+
+               else
+                  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/
+                     ! See also https://tc.copernicus.org/articles/19/1539/2025/
+                     ! 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 
+                        thk(c,j) = 0.138 - 1.01*(bw(c,j)/1000) +(3.233*((bw(c,j)/1000)*(bw(c,j)/1000))) 
+                     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 ! close land unit if statement 
+          end if
 
          end do
       end do

src/main/clm_varctl.F90

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

src/main/controlMod.F90

@@ -207,7 +207,8 @@ subroutine control_init(dtime)
          create_crop_landunit, nsegspc, co2_ppmv, &
          albice, soil_layerstruct_predefined, soil_layerstruct_userdefined, &
          soil_layerstruct_userdefined_nlevsoi, use_subgrid_fluxes, &
-         snow_thermal_cond_method, snow_cover_fraction_method, &
+         snow_thermal_cond_method, snow_thermal_cond_glc_method, &
+         snow_thermal_cond_lake_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, &
@@ -926,6 +927,8 @@ subroutine control_spmd()
     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_thermal_cond_glc_method, len(snow_thermal_cond_glc_method), MPI_CHARACTER, 0, mpicom, ier)
+    call mpi_bcast (snow_thermal_cond_lake_method, len(snow_thermal_cond_lake_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)