Adds a new snow compaction & routing scheme in ELM

components/elm/bld/namelist_files/namelist_definition.xml

@@ -709,6 +709,14 @@ snowpack conditions on glaciers and ice sheets, including a semi-empirical
 firn densification model.
 </entry>
 
+<entry id="use_firn_percolation_and_compaction" type="logical" category="physics"
+       group="elm_inparm" valid_values="" value=".false.">
+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>
+
 <entry id="use_noio" type="logical" category="default_settings"
        group="elm_inparm" valid_values="" value=".false." >
 Toggle to turn all history output completely OFF (possibly used for testing)

components/elm/cime_config/testdefs/testmods_dirs/elm/firn_percolation_compaction/user_nl_elm

@@ -0,0 +1 @@
+use_firn_percolation_and_compaction = .true.

components/elm/src/biogeophys/AerosolMod.F90

@@ -5,7 +5,7 @@ module AerosolMod
   use shr_log_mod      , only : errMsg => shr_log_errMsg
   use decompMod        , only : bounds_type
   use elm_varpar       , only : nlevsno 
-  use elm_varctl       , only : use_extrasnowlayers
+  use elm_varctl       , only : use_firn_percolation_and_compaction
   use elm_time_manager , only : get_step_size
   use atm2lndType      , only : atm2lnd_type
   use AerosolType      , only : aerosol_type
@@ -107,7 +107,7 @@ subroutine AerosolMasses(bounds, num_on, filter_on, num_off, filter_off, aerosol
             ! layer mass of snow:
             snowmass = h2osoi_ice(c,j) + h2osoi_liq(c,j)
 
-            if (.not. use_extrasnowlayers) then
+            if (.not. use_firn_percolation_and_compaction) then
                ! Correct the top layer aerosol mass to account for snow capping. 
                ! This approach conserves the aerosol mass concentration
                ! (but not the aerosol amss) when snow-capping is invoked

components/elm/src/biogeophys/BalanceCheckMod.F90

@@ -9,7 +9,7 @@ module BalanceCheckMod
   use shr_log_mod        , only : errMsg => shr_log_errMsg
   use decompMod          , only : bounds_type
   use abortutils         , only : endrun
-  use elm_varctl         , only : iulog, use_var_soil_thick, use_extrasnowlayers
+  use elm_varctl         , only : iulog, use_var_soil_thick, use_firn_percolation_and_compaction
   use elm_varcon         , only : namep, namec
   use GetGlobalValuesMod , only : GetGlobalIndex
   use atm2lndType        , only : atm2lnd_type
@@ -445,7 +445,7 @@ subroutine ColWaterBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
                      + qflx_snwcp_ice(c) + qflx_snwcp_liq(c) + qflx_sl_top_soil(c)
 
                 if (lun_pp%itype(l) == istdlak) then 
-                   if (.not. use_extrasnowlayers) then
+                   if (.not. use_firn_percolation_and_compaction) then
                       if ( do_capsnow(c)) then
                          snow_sources(c) = qflx_snow_grnd_col(c) &
                               + frac_sno_eff(c) * (qflx_dew_snow(c) + qflx_dew_grnd(c) ) 
@@ -473,7 +473,7 @@ subroutine ColWaterBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
                 endif
 
                 if (lun_pp%itype(l) == istsoil .or. lun_pp%itype(l) == istcrop .or. lun_pp%itype(l) == istwet ) then
-                   if (.not. use_extrasnowlayers) then
+                   if (.not. use_firn_percolation_and_compaction) then
                       if (do_capsnow(c)) then
                          snow_sources(c) = frac_sno_eff(c) * (qflx_dew_snow(c) + qflx_dew_grnd(c) ) &
                               + qflx_h2osfc_to_ice(c) + qflx_prec_grnd(c)

components/elm/src/biogeophys/CanopyHydrologyMod.F90

@@ -16,7 +16,7 @@ module CanopyHydrologyMod
   use shr_sys_mod       , only : shr_sys_flush
   use decompMod         , only : bounds_type
   use abortutils        , only : endrun
-  use elm_varctl        , only : iulog, tw_irr, extra_gw_irr, irrigate, use_extrasnowlayers
+  use elm_varctl        , only : iulog, tw_irr, extra_gw_irr, irrigate, use_firn_percolation_and_compaction
   use LandunitType      , only : lun_pp
   use atm2lndType       , only : atm2lnd_type
   use AerosolType       , only : aerosol_type
@@ -433,7 +433,7 @@ subroutine CanopyHydrology(bounds, &
 
           qflx_prec_grnd(p) = qflx_prec_grnd_snow(p) + qflx_prec_grnd_rain(p)
 
-          if (.not. use_extrasnowlayers) then
+          if (.not. use_firn_percolation_and_compaction) then
              if (do_capsnow(c)) then
                 qflx_snwcp_liq(p) = qflx_prec_grnd_rain(p)
                 qflx_snwcp_ice(p) = qflx_prec_grnd_snow(p)
@@ -491,7 +491,7 @@ subroutine CanopyHydrology(bounds, &
 
        ! Determine snow height and snow water
 
-       if (use_extrasnowlayers) then
+       if (use_firn_percolation_and_compaction) then
           call NewSnowBulkDensity(bounds, num_nolakec, filter_nolakec, &
                                   top_as, bifall(bounds%begc:bounds%endc))
        end if
@@ -517,13 +517,13 @@ subroutine CanopyHydrology(bounds, &
              swe_old(c,j)=h2osoi_liq(c,j)+h2osoi_ice(c,j)
           enddo
 
-          if (do_capsnow(c) .and. .not. use_extrasnowlayers) then
+          if (do_capsnow(c) .and. .not. use_firn_percolation_and_compaction) then
              dz_snowf = 0._r8
              newsnow(c) = qflx_snow_grnd_col(c) * dtime
              frac_sno(c)=1._r8
              int_snow(c) = 5.e2_r8
           else
-             if (.not. use_extrasnowlayers) then
+             if (.not. use_firn_percolation_and_compaction) then
                 if (forc_t(t) > tfrz + 2._r8) then
                    bifall(c)=50._r8 + 1.7_r8*(17.0_r8)**1.5_r8
                 else if (forc_t(t) > tfrz - 15._r8) then
@@ -662,7 +662,7 @@ subroutine CanopyHydrology(bounds, &
           ! as the surface air temperature
 
           newnode = 0    ! flag for when snow node will be initialized
-          if (.not. use_extrasnowlayers) then
+          if (.not. use_firn_percolation_and_compaction) then
              if (snl(c) == 0 .and. qflx_snow_grnd_col(c) > 0.0_r8 .and. frac_sno(c)*snow_depth(c) >= 0.01_r8) then
                 newnode = 1
                 snl(c) = -1

components/elm/src/biogeophys/HydrologyNoDrainageMod.F90

@@ -7,7 +7,7 @@ Module HydrologyNoDrainageMod
   use shr_kind_mod      , only : r8 => shr_kind_r8
   use shr_log_mod       , only : errMsg => shr_log_errMsg
   use decompMod         , only : bounds_type
-  use elm_varctl        , only : iulog, use_vichydro, use_extrasnowlayers
+  use elm_varctl        , only : iulog, use_vichydro, use_extrasnowlayers, use_firn_percolation_and_compaction
   use elm_varcon        , only : e_ice, denh2o, denice, rpi, spval
   use atm2lndType       , only : atm2lnd_type
   use lnd2atmType       , only : lnd2atm_type
@@ -300,7 +300,7 @@ subroutine HydrologyNoDrainage(bounds, &
       endif           
 #endif
 
-      if (use_extrasnowlayers) then
+      if (use_firn_percolation_and_compaction) then
          call SnowCapping(bounds, num_nolakec, filter_nolakec, num_snowc, filter_snowc, &
                           aerosol_vars)
       end if

components/elm/src/biogeophys/LakeFluxesMod.F90

@@ -52,7 +52,7 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep,
     use elm_varpar          , only : nlevlak
     use elm_varcon          , only : hvap, hsub, hfus, cpair, cpliq, tkwat, tkice, tkair
     use elm_varcon          , only : sb, vkc, grav, denh2o, tfrz, spval, zsno
-    use elm_varctl          , only : iulog, use_lch4, use_extrasnowlayers
+    use elm_varctl          , only : iulog, use_lch4, use_firn_percolation_and_compaction
     use LakeCon             , only : betavis, z0frzlake, tdmax, emg_lake
     use LakeCon             , only : lake_use_old_fcrit_minz0
     use LakeCon             , only : minz0lake, cur0, cus, curm, fcrit
@@ -695,7 +695,7 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep,
          qflx_dirct_rain(p) = 0._r8
          qflx_leafdrip(p) = 0._r8
 
-         if (.not. use_extrasnowlayers) then
+         if (.not. use_firn_percolation_and_compaction) then
             ! Because they will be used in pft2col initialize here.
             ! This will be overwritten in LakeHydrology
             qflx_snwcp_ice(p) = 0._r8

components/elm/src/biogeophys/LakeHydrologyMod.F90

@@ -74,7 +74,7 @@ subroutine LakeHydrology(bounds, &
       !$acc routine seq
     use elm_varcon      , only : denh2o, denice, spval, hfus, tfrz, cpliq, cpice
     use elm_varpar      , only : nlevsno, nlevgrnd, nlevsoi
-    use elm_varctl      , only : iulog, use_extrasnowlayers, use_lake_wat_storage
+    use elm_varctl      , only : iulog, use_extrasnowlayers, use_lake_wat_storage, use_firn_percolation_and_compaction
     use elm_time_manager, only : get_step_size
     use SnowHydrologyMod, only : SnowCompaction, CombineSnowLayers, SnowWater, BuildSnowFilter
     use SnowHydrologyMod, only : DivideSnowLayers, DivideExtraSnowLayers, SnowCapping
@@ -256,7 +256,7 @@ subroutine LakeHydrology(bounds, &
 
          qflx_dirct_rain(p) = 0._r8
          qflx_leafdrip(p) = 0._r8
-         if (.not. use_extrasnowlayers) then
+         if (.not. use_firn_percolation_and_compaction) then
             if (do_capsnow(c)) then
                qflx_snwcp_ice(p) = qflx_prec_grnd_snow(p)
                qflx_snwcp_liq(p) = qflx_prec_grnd_rain(p)
@@ -288,7 +288,7 @@ subroutine LakeHydrology(bounds, &
          ! U.S.Department of Agriculture Forest Service, Project F,
          ! Progress Rep. 1, Alta Avalanche Study Center:Snow Layer Densification.
 
-         if (do_capsnow(c) .and. .not. use_extrasnowlayers) then
+         if (do_capsnow(c) .and. .not. use_firn_percolation_and_compaction) then
             dz_snowf = 0._r8
          else
             if (forc_t(t) > tfrz + 2._r8) then
@@ -400,7 +400,7 @@ subroutine LakeHydrology(bounds, &
             ! Update the pft-level qflx_snowcap
             ! This was moved in from Hydrology2 to keep all pft-level
             ! calculations out of Hydrology2
-            if (do_capsnow(c) .and. .not. use_extrasnowlayers) then
+            if (do_capsnow(c) .and. .not. use_firn_percolation_and_compaction) then
                qflx_snwcp_ice(p) = qflx_snwcp_ice(p) + qflx_dew_snow(p) 
                qflx_snwcp_liq(p) = qflx_snwcp_liq(p) + qflx_dew_grnd(p)
             end if
@@ -422,7 +422,7 @@ subroutine LakeHydrology(bounds, &
             ! Update snow pack for dew & sub.
 
             h2osno_temp = h2osno(c)
-            if (do_capsnow(c) .and. .not. use_extrasnowlayers) then
+            if (do_capsnow(c) .and. .not. use_firn_percolation_and_compaction) then
                h2osno(c) = h2osno(c) - qflx_sub_snow(p)*dtime
                qflx_snwcp_ice(p) = qflx_snwcp_ice(p) + qflx_dew_snow(p)
                qflx_snwcp_liq(p) = qflx_snwcp_liq(p) + qflx_dew_grnd(p)
@@ -437,7 +437,7 @@ subroutine LakeHydrology(bounds, &
             end if
          end if
 
-         if (.not. use_extrasnowlayers) then
+         if (.not. use_firn_percolation_and_compaction) then
             qflx_snwcp_ice_col(c) = qflx_snwcp_ice(p)
             qflx_snwcp_liq_col(c) = qflx_snwcp_liq(p)
          end if
@@ -485,7 +485,7 @@ subroutine LakeHydrology(bounds, &
            num_shlakesnowc, filter_shlakesnowc, num_shlakenosnowc, filter_shlakenosnowc, &
            atm2lnd_vars,  aerosol_vars)
 
-      if (use_extrasnowlayers) then
+      if (use_firn_percolation_and_compaction) then
          call SnowCapping(bounds, num_lakec, filter_lakec, num_shlakesnowc, filter_shlakesnowc, &
                           aerosol_vars)
       end if
@@ -719,7 +719,7 @@ subroutine LakeHydrology(bounds, &
          if (use_lake_wat_storage) then
             qflx_qrgwl(c)     = 0._r8
             if (wslake(c) >= 5000._r8) then
-               if (.not. use_extrasnowlayers) then
+               if (.not. use_firn_percolation_and_compaction) then
                   qflx_snwcp = qflx_snwcp_ice(p)
                else
                   qflx_snwcp = qflx_snwcp_ice_col(c)
@@ -732,7 +732,7 @@ subroutine LakeHydrology(bounds, &
                 (endwb(c) - begwb(c))
             endwb(c) = endwb(c) + wslake(c)
          else
-            if (.not. use_extrasnowlayers) then
+            if (.not. use_firn_percolation_and_compaction) then
                qflx_snwcp = qflx_snwcp_ice(p)
             else
                qflx_snwcp = qflx_snwcp_ice_col(c)

components/elm/src/biogeophys/SnowHydrologyMod.F90

@@ -18,7 +18,7 @@ module SnowHydrologyMod
   use decompMod       , only : bounds_type
   use abortutils      , only : endrun
   use elm_varpar      , only : nlevsno
-  use elm_varctl      , only : iulog, use_extrasnowlayers
+  use elm_varctl      , only : iulog, use_extrasnowlayers, use_firn_percolation_and_compaction
   use elm_varcon      , only : namec, h2osno_max
   use atm2lndType     , only : atm2lnd_type
   use AerosolType     , only : aerosol_type
@@ -233,7 +233,7 @@ subroutine SnowWater(bounds, &
          c = filter_snowc(fc)
          l=col_pp%landunit(c)
 
-         if (do_capsnow(c) .and. .not. use_extrasnowlayers) then
+         if (do_capsnow(c) .and. .not. use_firn_percolation_and_compaction) then
             wgdif = h2osoi_ice(c,snl(c)+1) - frac_sno_eff(c)*qflx_sub_snow(c)*dtime
             h2osoi_ice(c,snl(c)+1) = wgdif
             if (wgdif < 0._r8) then
@@ -623,7 +623,7 @@ subroutine SnowCompaction(bounds, num_snowc, filter_snowc, &
 
        ! Begin calculation - note that the following column loops are only invoked if snl(c) < 0
 
-       if (use_extrasnowlayers) then
+       if (use_firn_percolation_and_compaction) then
           do fc = 1, num_snowc
              c = filter_snowc(fc)
              burden(c) = 0._r8
@@ -637,7 +637,7 @@ subroutine SnowCompaction(bounds, num_snowc, filter_snowc, &
        do j = -nlevsno+1, 0
           do fc = 1, num_snowc
              c = filter_snowc(fc)
-             if (use_extrasnowlayers) then
+             if (use_firn_percolation_and_compaction) then
                 t = col_pp%topounit(c)
              end if
              if (j >= snl(c)+1) then
@@ -647,7 +647,7 @@ subroutine SnowCompaction(bounds, num_snowc, filter_snowc, &
                 ! If void is negative, then increase dz such that void = 0.
                 ! This should be done for any landunit, but for now is done only for glacier_mec 1andunits.
 
-                if (.not. use_extrasnowlayers) then
+                if (.not. use_firn_percolation_and_compaction) then
                    ! I don't think the next 5 lines are necessary (removed in CLMv5)
                    l = col_pp%landunit(c)
                    if (ltype(l)==istice_mec .and. void < 0._r8) then
@@ -665,7 +665,7 @@ subroutine SnowCompaction(bounds, num_snowc, filter_snowc, &
                    dexpf = exp(-c4*td)
 
                    ! Settling as a result of destructive metamorphism
-                   if (.not. use_extrasnowlayers) then
+                   if (.not. use_firn_percolation_and_compaction) then
                       ddz1 = -c3*dexpf 
                       if (bi > dm) ddz1 = ddz1*exp(-46.0e-3_r8*(bi-dm))
                    else
@@ -685,7 +685,7 @@ subroutine SnowCompaction(bounds, num_snowc, filter_snowc, &
                    if (h2osoi_liq(c,j) > 0.01_r8*dz(c,j)*frac_sno(c)) ddz1=ddz1*c5
 
                    ! Compaction due to overburden
-                   if (.not. use_extrasnowlayers) then
+                   if (.not. use_firn_percolation_and_compaction) then
                       ddz2 = -(burden(c)+wx/2._r8)*exp(-0.08_r8*td - c2*bi)/eta0 
                    else
                       p_gls = max(denice / bi, 1._r8) * grav * (burden(c) + wx/2._r8)
@@ -724,7 +724,7 @@ subroutine SnowCompaction(bounds, num_snowc, filter_snowc, &
                       ddz3 = 0._r8
                    end if
 
-                   if (use_extrasnowlayers) then
+                   if (use_firn_percolation_and_compaction) then
                       ! Compaction occurring due to wind drift
                       call WindDriftCompaction( &
                            bi = bi, &

components/elm/src/biogeophys/SnowSnicarMod.F90

@@ -11,7 +11,7 @@ module SnowSnicarMod
   use shr_kind_mod    , only : r8 => shr_kind_r8
   use shr_sys_mod     , only : shr_sys_flush
   use shr_log_mod     , only : errMsg => shr_log_errMsg
-  use elm_varctl      , only : iulog, use_extrasnowlayers
+  use elm_varctl      , only : iulog, use_firn_percolation_and_compaction
   use elm_varcon      , only : namec 
   use shr_const_mod   , only : SHR_CONST_RHOICE
   use abortutils      , only : endrun
@@ -1435,13 +1435,13 @@ subroutine SnowAge_grain(bounds, &
             !               RE-FREEZING
             !
             ! new snowfall [kg/m2]
-            if (do_capsnow(c_idx) .and. .not. use_extrasnowlayers) then
+            if (do_capsnow(c_idx) .and. .not. use_firn_percolation_and_compaction) then
                newsnow = max(0._r8, (qflx_snwcp_ice(c_idx)*dtime))
             else
                newsnow = max(0._r8, (qflx_snow_grnd_col(c_idx)*dtime))
             endif
 
-            if (use_extrasnowlayers) then
+            if (use_firn_percolation_and_compaction) then
                 snw_rds_refrz   = 1500._r8
             else
                 snw_rds_refrz   = 1000._r8 

components/elm/src/biogeophys/SoilFluxesMod.F90

@@ -8,7 +8,7 @@ module SoilFluxesMod
   use shr_log_mod	, only : errMsg => shr_log_errMsg
   use decompMod		, only : bounds_type
   use abortutils	, only : endrun
-  use elm_varctl	, only : iulog, use_extrasnowlayers
+  use elm_varctl	, only : iulog, use_firn_percolation_and_compaction
   use perfMod_GPU
   use elm_varpar	, only : nlevsno, nlevgrnd, nlevurb, max_patch_per_col
   use atm2lndType	, only : atm2lnd_type
@@ -349,7 +349,7 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, &
          ! This was moved in from Hydrology2 to keep all pft-level
          ! calculations out of Hydrology2
 
-         if (col_pp%snl(c) < 0 .and. do_capsnow(c) .and. .not. use_extrasnowlayers) then
+         if (col_pp%snl(c) < 0 .and. do_capsnow(c) .and. .not. use_firn_percolation_and_compaction) then
             qflx_snwcp_liq(p) = qflx_snwcp_liq(p)+frac_sno_eff(c)*qflx_dew_grnd(p)
             qflx_snwcp_ice(p) = qflx_snwcp_ice(p)+frac_sno_eff(c)*qflx_dew_snow(p)
          end if

components/elm/src/data_types/VegetationDataType.F90

@@ -5494,20 +5494,20 @@ subroutine veg_wf_init(this, begp, endp)
          avgflag='A', long_name='excess rainfall due to snow capping', &
          ptr_patch=this%qflx_snwcp_liq, c2l_scale_type='urbanf', default='inactive')
 
-    if (use_cn) then
+    !if (use_cn) then
        this%qflx_rain_grnd(begp:endp) = spval
        call hist_addfld1d (fname='QFLX_RAIN_GRND', units='mm H2O/s', &
             avgflag='A', long_name='rain on ground after interception', &
             ptr_patch=this%qflx_rain_grnd, default='inactive', c2l_scale_type='urbanf')
 
-    end if
+    !end if
 
-    if (use_cn) then
+    !if (use_cn) then
        this%qflx_snow_grnd(begp:endp) = spval
        call hist_addfld1d (fname='QFLX_SNOW_GRND', units='mm H2O/s', &
             avgflag='A', long_name='snow on ground after interception', &
             ptr_patch=this%qflx_snow_grnd, default='inactive', c2l_scale_type='urbanf')
-    end if
+    !end if
 
     if (use_cn) then
        this%qflx_evap_grnd(begp:endp) = spval
@@ -5537,12 +5537,12 @@ subroutine veg_wf_init(this, begp, endp)
             ptr_patch=this%qflx_dew_grnd, default='inactive', c2l_scale_type='urbanf')
     end if
 
-    if (use_cn) then
+    !if (use_cn) then
        this%qflx_sub_snow(begp:endp) = spval
        call hist_addfld1d (fname='QFLX_SUB_SNOW', units='mm H2O/s', &
             avgflag='A', long_name='sublimation rate from snow pack', &
             ptr_patch=this%qflx_sub_snow, default='inactive', c2l_scale_type='urbanf')
-    end if
+    !end if
 
     if (use_cn) then
        this%qflx_dew_snow(begp:endp) = spval