Merge branch 'cwhicker/elm/use_extrasnowlayers-tuning' (PR #6298)

Makes the following changes for `use_extrasnowlayers = .true.`:
- Edits to tuning parameters to follow recommendations from Schneider et al., 2021,
- Changes the snow refrozen grain radius (from 1000um to 1500um), and
- Increase the max depth of snow (from 1m to 30m)
- Adds three additional fields in ELM output by default.

[BFB] except for tests longer than one month

components/elm/src/biogeophys/SnowHydrologyMod.F90

@@ -571,14 +571,14 @@ subroutine SnowCompaction(bounds, num_snowc, filter_snowc, &
      ! parameters
      real(r8), parameter :: c2 = 23.e-3_r8       ! [m3/kg]
      real(r8), parameter :: c3 = 2.777e-6_r8     ! [1/s]
-     real(r8), parameter :: c3_ams = 5.8e-7_r8   ! Schneider et al., 2020 [1/s]
+     real(r8), parameter :: c3_ams = 0.83e-6_r8   ! Schneider et al.,(2021),Table 2 [1/s]
      real(r8), parameter :: c4 = 0.04_r8         ! [1/K]
      real(r8), parameter :: c5 = 2.0_r8          !
      real(r8), parameter :: dm = 100.0_r8        ! Upper Limit on Destructive Metamorphism Compaction [kg/m3]
-     real(r8), parameter :: rho_dm = 150.0_r8    ! Upper limit on destructive metamorphism compaction [kg/m3] (Anderson, 1976; Schneider et al., 2020)
+     real(r8), parameter :: rho_dm = 150.0_r8    ! Upper limit on destructive metamorphism compaction [kg/m3] (Anderson, 1976;Schneider et al., 2021)
      real(r8), parameter :: eta0 = 9.e+5_r8      ! The Viscosity Coefficient Eta0 [kg-s/m2]
-     real(r8), parameter :: k_creep_snow = 1.4e-9_r8 ! Creep coefficient for snow (bi < 550 kg / m3) [m3-s/kg]
-     real(r8), parameter :: k_creep_firn = 1.2e-9_r8 ! Creep coefficient for firn (bi > 550 kg / m3)
+     real(r8), parameter :: k_creep_snow = 9.2e-9_r8 ! Creep coefficient for snow (bi < 550 kg / m3) [m3-s/kg]
+     real(r8), parameter :: k_creep_firn = 3.7e-9_r8 ! Creep coefficient for firn (bi > 550 kg / m3) [m3-s/kg]
      !
      real(r8) :: p_gls                           ! grain load stress [kg / m-s2]
      real(r8) :: burden(bounds%begc:bounds%endc) ! pressure of overlying snow [kg/m2]
@@ -693,12 +693,12 @@ subroutine SnowCompaction(bounds, num_snowc, filter_snowc, &
                          ddz2 = (-k_creep_snow * (max(denice / bi, 1._r8) - 1._r8) * &
                                  exp(-60.e6_r8 / (rgas * t_soisno(c,j))) * p_gls) / &
                                  (snw_rds(c,j) * 1.e-6_r8 * snw_rds(c,j) * 1.e-6_r8) - &
-                                 2.02e-10_r8
+                                 1.0e-10_r8
                       else ! High density, i.e. firn
                          ddz2 = (-k_creep_firn * (max(denice / bi, 1._r8) - 1._r8) * &
                                  exp(-60.e6_r8 / (rgas * t_soisno(c,j))) * p_gls) / &
                                  (snw_rds(c,j) * 1.e-6_r8 * snw_rds(c,j) * 1.e-6_r8) - &
-                                 2.7e-11_r8
+                                 1.0e-10_r8
                       endif
                    endif
 

components/elm/src/biogeophys/SnowSnicarMod.F90

@@ -80,7 +80,7 @@ module SnowSnicarMod
   integer,  parameter :: snw_rds_max_tbl = 1500          ! maximum effective radius defined in Mie lookup table [microns]
   integer,  parameter :: snw_rds_min_tbl = 30            ! minimium effective radius defined in Mie lookup table [microns]
   real(r8), parameter :: snw_rds_max     = 1500._r8      ! maximum allowed snow effective radius [microns]
-  real(r8), parameter :: snw_rds_refrz   = 1000._r8      ! effective radius of re-frozen snow [microns]
+  real(r8)            :: snw_rds_refrz   = 1000._r8      ! effective radius of re-frozen snow [microns]
   !$acc declare copyin(snw_rds_max_tbl)
   !$acc declare copyin(snw_rds_min_tbl)
   !$acc declare copyin(snw_rds_max    )
@@ -1441,6 +1441,12 @@ subroutine SnowAge_grain(bounds, &
                newsnow = max(0._r8, (qflx_snow_grnd_col(c_idx)*dtime))
             endif
 
+            if (use_extrasnowlayers) then
+                snw_rds_refrz   = 1500._r8
+            else
+                snw_rds_refrz   = 1000._r8 
+            endif
+
             ! snow that has re-frozen [kg/m2]
             refrzsnow = max(0._r8, (qflx_snofrz_lyr(c_idx,i)*dtime))
 

components/elm/src/data_types/ColumnDataType.F90

@@ -18,7 +18,7 @@ module ColumnDataType
   use elm_varpar      , only : nlevdecomp_full, crop_prog, nlevdecomp
   use elm_varpar      , only : i_met_lit, i_cel_lit, i_lig_lit, i_cwd
   use elm_varcon      , only : spval, ispval, zlnd, snw_rds_min, denice, denh2o, tfrz, pondmx
-  use elm_varcon      , only : watmin, bdsno, zsoi, zisoi, dzsoi_decomp
+  use elm_varcon      , only : watmin, bdsno, bdfirn, zsoi, zisoi, dzsoi_decomp
   use elm_varcon      , only : c13ratio, c14ratio, secspday
   use elm_varctl      , only : use_fates, use_fates_planthydro, create_glacier_mec_landunit
   use elm_varctl      , only : use_hydrstress, use_crop
@@ -28,6 +28,7 @@ module ColumnDataType
   use elm_varctl      , only : hist_wrtch4diag, use_century_decomp
   use elm_varctl      , only : get_carbontag, override_bgc_restart_mismatch_dump
   use elm_varctl      , only : pf_hmode, nu_com
+  use elm_varctl      , only : use_extrasnowlayers
   use elm_varctl      , only : use_fan
   use ch4varcon       , only : allowlakeprod
   use pftvarcon       , only : VMAX_MINSURF_P_vr, KM_MINSURF_P_vr, pinit_beta1, pinit_beta2
@@ -1742,8 +1743,20 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_
           end do
           do j = -nlevsno+1, 0
              if (j > col_pp%snl(c)) then
-                this%h2osoi_ice(c,j) = col_pp%dz(c,j)*250._r8
-                this%h2osoi_liq(c,j) = 0._r8
+                    if (use_extrasnowlayers) then
+                       ! amschnei@uci.edu: Initialize "deep firn" on glacier columns
+                       if (lun_pp%itype(l) == istice .or. lun_pp%itype(l) == istice_mec) then
+                         this%h2osoi_ice(c,j) = col_pp%dz(c,j)*bdfirn
+                         this%h2osoi_liq(c,j) = 0._r8
+                       else
+                         this%h2osoi_ice(c,j) = col_pp%dz(c,j)*bdsno
+                         this%h2osoi_liq(c,j) = 0._r8
+                       end if
+                    else ! no firn model (default in v2)
+                       ! Below, "250._r8" should instead be "bdsno", which is 250 kg m^3 by default
+                       this%h2osoi_ice(c,j) = col_pp%dz(c,j)*250._r8
+                       this%h2osoi_liq(c,j) = 0._r8
+                    end if
              end if
           end do
        end if
@@ -5830,26 +5843,20 @@ subroutine col_wf_init(this, begc, endc)
           avgflag='A', long_name='column-integrated snow freezing rate', &
            ptr_col=this%qflx_snofrz, set_lake=spval, c2l_scale_type='urbanf', default='inactive')
 
-    if (create_glacier_mec_landunit) then
-       this%qflx_glcice(begc:endc) = spval
-        call hist_addfld1d (fname='QICE',  units='mm/s',  &
-             avgflag='A', long_name='ice growth/melt', &
-              ptr_col=this%qflx_glcice, l2g_scale_type='ice')
-    end if
-
-    if (create_glacier_mec_landunit) then
-       this%qflx_glcice_frz(begc:endc) = spval
-        call hist_addfld1d (fname='QICE_FRZ',  units='mm/s',  &
-             avgflag='A', long_name='ice growth', &
-              ptr_col=this%qflx_glcice_frz, l2g_scale_type='ice')
-    end if
+    this%qflx_glcice(begc:endc) = spval
+     call hist_addfld1d (fname='QICE',  units='mm/s',  &
+          avgflag='A', long_name='ice growth/melt', &
+           ptr_col=this%qflx_glcice, l2g_scale_type='ice')
 
-    if (create_glacier_mec_landunit) then
-       this%qflx_glcice_melt(begc:endc) = spval
-        call hist_addfld1d (fname='QICE_MELT',  units='mm/s',  &
-             avgflag='A', long_name='ice melt', &
-              ptr_col=this%qflx_glcice_melt, l2g_scale_type='ice')
-    end if
+    this%qflx_glcice_frz(begc:endc) = spval
+     call hist_addfld1d (fname='QICE_FRZ',  units='mm/s',  &
+          avgflag='A', long_name='ice growth', &
+           ptr_col=this%qflx_glcice_frz, l2g_scale_type='ice')
+
+    this%qflx_glcice_melt(begc:endc) = spval
+     call hist_addfld1d (fname='QICE_MELT',  units='mm/s',  &
+          avgflag='A', long_name='ice melt', &
+           ptr_col=this%qflx_glcice_melt, l2g_scale_type='ice')
 
     ! As defined here, snow_sources - snow_sinks will equal the change in h2osno at any
     ! given time step but only if there is at least one snow layer (for all landunits

components/elm/src/main/elm_instMod.F90

@@ -262,7 +262,7 @@ subroutine elm_inst_biogeophys(bounds_proc)
     !
     use shr_scam_mod                      , only : shr_scam_getCloseLatLon
     use landunit_varcon                   , only : istice, istice_mec, istsoil
-    use elm_varcon                        , only : h2osno_max, bdsno
+    use elm_varcon                        , only : h2osno_max, bdsno, bdfirn
     use domainMod                         , only : ldomain
     use elm_varpar                        , only : nlevsno, numpft
     use elm_varctl                        , only : single_column, fsurdat, scmlat, scmlon, use_extrasnowlayers
@@ -318,8 +318,18 @@ subroutine elm_inst_biogeophys(bounds_proc)
            else
               h2osno_col(c) = 0._r8
            endif
-       else ! With a deeper firn model, we can no longer depend on "h2osno_max," because this will
-           !  potentially be large, resulting in a lot of artificial firn at initialization.
+           snow_depth_col(c) = h2osno_col(c) / bdsno
+       else! With a deeper firn model, we can no longer depend on "h2osno_max," because this will
+           ! potentially be large, resulting in a lot of artificial firn at initialization.          ! 
+           ! amschnei@uci.edu: UPDATE - By including a new parameter, bdfirn, we can
+           ! initialize a dense (e.g, 730 kd m^-3) "deep firn" layer for glacier
+           ! land units. This deep firn layer can be set to half of the total
+           ! maximum allotted snowpack mass (i.e., "h2osno_max"), which effectively
+           ! sets the Greenland and Antarctic ice sheets' (and mountain glaciers')
+           ! climatic (aka surface) mass balance (SMB) initial condition to 0. With this
+           ! cold start condition, a multi-decadal (100 to 300 years or more)
+           ! spin up is first required to prognose nonzero SMB.
+           ! 
            ! However... (below docstring from CLMv5)
            ! In areas that should be snow-covered, it can be problematic to start with 0 snow
            ! cover, because this can affect the long-term state through soil heating, albedo
@@ -328,14 +338,18 @@ subroutine elm_inst_biogeophys(bounds_proc)
            ! feedback may not activate on time (or at all). So, as a compromise, we start with
            ! a small amount of snow in places that are likely to be snow-covered for much or
            ! all of the year.
+           ! amschnei@uci.edu: Initializing "deep firn" for glacier columns
            if (lun_pp%itype(l)==istice .or. lun_pp%itype(l)==istice_mec) then
               ! land ice (including multiple elevation classes, i.e. glacier_mec)
-              h2osno_col(c) = 50._r8
+              h2osno_col(c) = 0.5_r8*h2osno_max   ! start with half full snow column, representing deep firn
+              snow_depth_col(c)  = h2osno_col(c) / bdfirn
            else if (lun_pp%itype(l)==istsoil .and. grc_pp%latdeg(g) >= 44._r8) then
               ! Northern hemisphere seasonal snow
               h2osno_col(c) = 50._r8
+              snow_depth_col(c) = h2osno_col(c) / bdsno
            else
               h2osno_col(c) = 0._r8
+              snow_depth_col(c) = h2osno_col(c) / bdsno
            endif
        endif
        snow_depth_col(c)  = h2osno_col(c) / bdsno

components/elm/src/main/elm_varcon.F90

@@ -67,6 +67,7 @@ module elm_varcon
   real(r8) :: oneatm = 1.01325e5_r8                         ! one standard atmospheric pressure [Pa]
 
   real(r8) :: bdsno = 250._r8                               ! bulk density snow (kg/m**3)
+  real(r8) :: bdfirn = 730._r8                              ! bulk density of deep firn (kg/m**3)
   real(r8) :: alpha_aero = 1.0_r8                           ! constant for aerodynamic parameter weighting
   real(r8) :: tlsai_crit = 2.0_r8                           ! critical value of elai+esai for which aerodynamic parameters are maximum
   real(r8) :: watmin = 0.01_r8                              ! minimum soil moisture (mm)
@@ -248,7 +249,7 @@ subroutine elm_varcon_init()
     allocate( dzsoifl(1:nlevsoifl            ))
 
     if (use_extrasnowlayers) then
-       h2osno_max = 10000._r8
+       h2osno_max = 30000._r8
     end if
 
   end subroutine elm_varcon_init