From 1274e0add1c3b8782e31d7fdacce04fe607db224 Mon Sep 17 00:00:00 2001 From: Chloe Date: Fri, 8 Mar 2024 11:27:10 -0800 Subject: [PATCH 1/4] adjusted use_extrasnowlayers parameters based on Schneider et al., 2021 --- .../elm/src/biogeophys/SnowHydrologyMod.F90 | 12 +++---- .../elm/src/biogeophys/SnowSnicarMod.F90 | 2 +- .../elm/src/data_types/ColumnDataType.F90 | 31 +++++++++++++------ components/elm/src/main/elm_instMod.F90 | 20 ++++++++++-- components/elm/src/main/elm_varcon.F90 | 3 +- 5 files changed, 48 insertions(+), 20 deletions(-) diff --git a/components/elm/src/biogeophys/SnowHydrologyMod.F90 b/components/elm/src/biogeophys/SnowHydrologyMod.F90 index 62f2ba24242a..6354dd598b0c 100644 --- a/components/elm/src/biogeophys/SnowHydrologyMod.F90 +++ b/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 diff --git a/components/elm/src/biogeophys/SnowSnicarMod.F90 b/components/elm/src/biogeophys/SnowSnicarMod.F90 index fae56d8927d7..abaa72d0a5f8 100644 --- a/components/elm/src/biogeophys/SnowSnicarMod.F90 +++ b/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), parameter :: snw_rds_refrz = 1500._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 ) diff --git a/components/elm/src/data_types/ColumnDataType.F90 b/components/elm/src/data_types/ColumnDataType.F90 index bc6c935efa5c..1793628598ef 100644 --- a/components/elm/src/data_types/ColumnDataType.F90 +++ b/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,26 @@ 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 + !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 + !end if - if (create_glacier_mec_landunit) then + !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 + !end if - if (create_glacier_mec_landunit) then + !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 + !end if ! 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 diff --git a/components/elm/src/main/elm_instMod.F90 b/components/elm/src/main/elm_instMod.F90 index ce5fb9bda4e4..d14d6ba3f101 100644 --- a/components/elm/src/main/elm_instMod.F90 +++ b/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 + 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. + ! 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 diff --git a/components/elm/src/main/elm_varcon.F90 b/components/elm/src/main/elm_varcon.F90 index 3a95fc05c52e..2fcf09dbf8a6 100644 --- a/components/elm/src/main/elm_varcon.F90 +++ b/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) @@ -130,7 +131,7 @@ module elm_varcon real(r8) :: ac_wasteheat_factor = 0.0_r8 !wasteheat factor for urban air conditioning (-) real(r8) :: wasteheat_limit = 100._r8 !limit on wasteheat (W/m2) - real(r8) :: h2osno_max = 1000._r8 ! max allowed snow thickness (mm H2O) + real(r8) :: h2osno_max = 30000._r8 ! max allowed snow thickness (mm H2O) real(r8), parameter :: lapse_glcmec = 0.006_r8 ! surface temperature lapse rate (deg m-1) ! Pritchard et al. (GRL, 35, 2008) use 0.006 real(r8), parameter :: glcmec_rain_snow_threshold = SHR_CONST_TKFRZ ! temperature dividing rain & snow in downscaling (K) From 67625371c49535701ece5252d318044f443601ec Mon Sep 17 00:00:00 2001 From: Chloe Date: Thu, 21 Mar 2024 15:17:36 -0700 Subject: [PATCH 2/4] stylistic code mods --- .../elm/src/data_types/ColumnDataType.F90 | 58 +++++++++---------- components/elm/src/main/elm_instMod.F90 | 22 +++---- 2 files changed, 37 insertions(+), 43 deletions(-) diff --git a/components/elm/src/data_types/ColumnDataType.F90 b/components/elm/src/data_types/ColumnDataType.F90 index 1793628598ef..2d79412e568c 100644 --- a/components/elm/src/data_types/ColumnDataType.F90 +++ b/components/elm/src/data_types/ColumnDataType.F90 @@ -1744,19 +1744,19 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_ do j = -nlevsno+1, 0 if (j > col_pp%snl(c)) then 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) + ! 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 + 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 @@ -5843,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 - - !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(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') + + 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 diff --git a/components/elm/src/main/elm_instMod.F90 b/components/elm/src/main/elm_instMod.F90 index d14d6ba3f101..b8afca87e241 100644 --- a/components/elm/src/main/elm_instMod.F90 +++ b/components/elm/src/main/elm_instMod.F90 @@ -319,17 +319,17 @@ subroutine elm_inst_biogeophys(bounds_proc) h2osno_col(c) = 0._r8 endif 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. - ! + 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 From db67eb039f0da97fa141c9650c262f6d2274fdf2 Mon Sep 17 00:00:00 2001 From: Chloe Date: Tue, 26 Mar 2024 12:47:41 -0700 Subject: [PATCH 3/4] mod to use_extrasnowlayers tuning to be BFB --- components/elm/src/biogeophys/SnowSnicarMod.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/components/elm/src/biogeophys/SnowSnicarMod.F90 b/components/elm/src/biogeophys/SnowSnicarMod.F90 index abaa72d0a5f8..b0965762207b 100644 --- a/components/elm/src/biogeophys/SnowSnicarMod.F90 +++ b/components/elm/src/biogeophys/SnowSnicarMod.F90 @@ -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)) From 0bbb3b9153886a90aca685fc53b773a95e2d56a4 Mon Sep 17 00:00:00 2001 From: Chloe Date: Mon, 1 Apr 2024 13:08:57 -0700 Subject: [PATCH 4/4] bug fixes to get B4B --- components/elm/src/biogeophys/SnowSnicarMod.F90 | 2 +- components/elm/src/main/elm_varcon.F90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/components/elm/src/biogeophys/SnowSnicarMod.F90 b/components/elm/src/biogeophys/SnowSnicarMod.F90 index b0965762207b..d6d6de75ffaa 100644 --- a/components/elm/src/biogeophys/SnowSnicarMod.F90 +++ b/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 = 1500._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 ) diff --git a/components/elm/src/main/elm_varcon.F90 b/components/elm/src/main/elm_varcon.F90 index 2fcf09dbf8a6..01df6ce0afd2 100644 --- a/components/elm/src/main/elm_varcon.F90 +++ b/components/elm/src/main/elm_varcon.F90 @@ -131,7 +131,7 @@ module elm_varcon real(r8) :: ac_wasteheat_factor = 0.0_r8 !wasteheat factor for urban air conditioning (-) real(r8) :: wasteheat_limit = 100._r8 !limit on wasteheat (W/m2) - real(r8) :: h2osno_max = 30000._r8 ! max allowed snow thickness (mm H2O) + real(r8) :: h2osno_max = 1000._r8 ! max allowed snow thickness (mm H2O) real(r8), parameter :: lapse_glcmec = 0.006_r8 ! surface temperature lapse rate (deg m-1) ! Pritchard et al. (GRL, 35, 2008) use 0.006 real(r8), parameter :: glcmec_rain_snow_threshold = SHR_CONST_TKFRZ ! temperature dividing rain & snow in downscaling (K) @@ -249,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