Skip to content

Commit

Permalink
Allow users to specify snow layer structure in namelist
Browse files Browse the repository at this point in the history 
Instead of hardcoding, we now make
dzmin(1) = dzmin_1
dzmax_l(1) = dzmax_l
dzmax_u(1) = dzmax_u
where the right-hand-side variables get namelist-defined values and
the model calculates dzmin and dzmax_* for the remaining snow layers.
The calculation is based on formulas written out in a file titled
"A firn model for CLM" that has been linked to ESCOMP#729.
  • Loading branch information
@slevis-lmwg
slevis-lmwg committed Aug 26, 2019
1 parent 20cbf74 commit 3e5856a
Show file tree
Hide file tree
Showing 2 changed files with 55 additions and 21 deletions.
68 changes: 48 additions & 20 deletions src/biogeophys/SnowHydrologyMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ module SnowHydrologyMod
use clm_varpar , only : nlevsno, nlevsoi, nlevgrnd
use clm_varctl , only : iulog, use_subgrid_fluxes
use clm_varcon , only : namec, h2osno_max, hfus, denice, rpi, spval, tfrz
use clm_varcon, only: dzmin, dzmax_l, dzmax_u
use atm2lndType , only : atm2lnd_type
use AerosolMod , only : aerosol_type
use TemperatureType , only : temperature_type
Expand Down Expand Up @@ -100,20 +101,6 @@ module SnowHydrologyMod
integer, parameter, public :: LoTmpDnsSlater2017 = 2 ! For temperature below -15C use equation from Slater 2017
integer, parameter, public :: LoTmpDnsTruncatedAnderson1976 = 1 ! Truncate low temp. snow density from the Anderson-1976 version at -15C

! Definition of snow pack vertical structure
! Hardcoded maximum of 12 snowlayers, this is checked elsewhere (controlMod.F90)
! The bottom layer has no limit on thickness, hence the last element of the dzmax_*
! arrays is 'huge'.
real(r8), parameter :: dzmin(12) = & ! minimum of top snow layer
(/ 0.010_r8, 0.015_r8, 0.025_r8, 0.055_r8, 0.115_r8, 0.235_r8, &
0.475_r8, 0.955_r8, 1.915_r8, 3.835_r8, 7.675_r8, 15.355_r8 /)
real(r8), parameter :: dzmax_l(12) = & ! maximum thickness of layer when no layers beneath
(/ 0.03_r8, 0.07_r8, 0.18_r8, 0.41_r8, 0.88_r8, 1.83_r8, &
3.74_r8, 7.57_r8, 15.24_r8, 30.59_r8, 61.3_r8, huge(1._r8) /)
real(r8), parameter :: dzmax_u(12) = & ! maximum thickness of layer when layers beneath
(/ 0.02_r8, 0.05_r8, 0.11_r8, 0.23_r8, 0.47_r8, 0.95_r8, &
1.91_r8, 3.83_r8, 7.67_r8, 15.35_r8, 30.71_r8, huge(1._r8) /)

!
! !PRIVATE DATA MEMBERS:

Expand All @@ -123,6 +110,7 @@ module SnowHydrologyMod
! If true, the density of new snow depends on wind speed, and there is also
! wind-dependent snow compaction
logical :: wind_dependent_snow_density ! If snow density depends on wind or not
real(r8) :: dzmin_1, dzmax_l_1, dzmax_u_1 ! namelist-defined top snow layer information
integer :: overburden_compaction_method = -1
integer :: new_snow_density = LoTmpDnsSlater2017 ! Snow density type
real(r8) :: upplim_destruct_metamorph = 100.0_r8 ! Upper Limit on Destructive Metamorphism Compaction [kg/m3]
Expand Down Expand Up @@ -187,11 +175,15 @@ subroutine SnowHydrology_readnl( NLFilename)
wind_dependent_snow_density, snow_overburden_compaction_method, &
lotmp_snowdensity_method, upplim_destruct_metamorph, &
overburden_compress_Tfactor, min_wind_snowcompact, &
reset_snow, reset_snow_glc, reset_snow_glc_ela
reset_snow, reset_snow_glc, reset_snow_glc_ela, dzmin_1, dzmax_l_1, &
dzmax_u_1

! Initialize options to default values, in case they are not specified in the namelist
wind_dependent_snow_density = .false.
snow_overburden_compaction_method = ' '
dzmin_1 = 0.01_r8
dzmax_l_1 = 0.03_r8
dzmax_u_1 = 0.02_r8

if (masterproc) then
unitn = getavu()
Expand All @@ -218,6 +210,9 @@ subroutine SnowHydrology_readnl( NLFilename)
call shr_mpi_bcast (reset_snow , mpicom)
call shr_mpi_bcast (reset_snow_glc , mpicom)
call shr_mpi_bcast (reset_snow_glc_ela , mpicom)
call shr_mpi_bcast (dzmin_1, mpicom)
call shr_mpi_bcast (dzmax_l_1, mpicom)
call shr_mpi_bcast (dzmax_u_1, mpicom)

if (masterproc) then
write(iulog,*) ' '
Expand All @@ -243,6 +238,19 @@ subroutine SnowHydrology_readnl( NLFilename)
errMsg(sourcefile, __LINE__))
end if

if (dzmin_1 < 0.01_r8) then
call endrun(msg="ERROR dzmin_1 cannot be less than 0.01"// &
errMsg(sourcefile, __LINE__))
end if
if (dzmax_l_1 < 0.03_r8) then
call endrun(msg="ERROR dzmax_l_1 cannot be less than 0.03"// &
errMsg(sourcefile, __LINE__))
end if
if (dzmax_u_1 < 0.02_r8) then
call endrun(msg="ERROR dzmax_2_1 cannot be less than 0.02"// &
errMsg(sourcefile, __LINE__))
end if

end subroutine SnowHydrology_readnl

!-----------------------------------------------------------------------
Expand Down Expand Up @@ -800,12 +808,12 @@ subroutine BuildFilter_SnowpackInitialized(bounds, num_c, filter_c, &
! maintain answers as before.
snowpack_initialized(c) = ( &
snl(c) == 0 .and. &
frac_sno_eff(c)*snow_depth(c) >= (0.01_r8 + lsadz) .and. &
frac_sno_eff(c)*snow_depth(c) >= (dzmin(1) + lsadz) .and. &
qflx_snow_grnd(c) > 0.0_r8)
else
snowpack_initialized(c) = ( &
snl(c) == 0 .and. &
frac_sno_eff(c)*snow_depth(c) >= 0.01_r8)
frac_sno_eff(c)*snow_depth(c) >= dzmin(1))
end if

end do
Expand Down Expand Up @@ -1496,7 +1504,7 @@ subroutine CombineSnowLayers(bounds, num_snowc, filter_snowc, &
real(r8):: h2osno_total(bounds%begc:bounds%endc) ! total snow water (mm H2O)
real(r8):: zwice(bounds%begc:bounds%endc) ! total ice mass in snow
real(r8):: zwliq (bounds%begc:bounds%endc) ! total liquid water in snow
real(r8):: dzminloc(size(dzmin)) ! minimum of top snow layer (local)
real(r8):: dzminloc(nlevsno) ! minimum of top snow layer (local)
real(r8):: dtime !land model time step (sec)

!-----------------------------------------------------------------------
Expand Down Expand Up @@ -1680,8 +1688,8 @@ subroutine CombineSnowLayers(bounds, num_snowc, filter_snowc, &
c = filter_snowc(fc)
l = col%landunit(c)
if (snow_depth(c) > 0._r8) then
if ((ltype(l) == istdlak .and. snow_depth(c) < 0.01_r8 + lsadz ) .or. &
((ltype(l) /= istdlak) .and. ((frac_sno_eff(c)*snow_depth(c) < 0.01_r8) &
if ((ltype(l) == istdlak .and. snow_depth(c) < dzmin(1) + lsadz ) .or. &
((ltype(l) /= istdlak) .and. ((frac_sno_eff(c)*snow_depth(c) < dzmin(1)) &
.or. (h2osno_total(c)/(frac_sno_eff(c)*snow_depth(c)) < 50._r8)))) then

snl(c) = 0
Expand Down Expand Up @@ -2234,6 +2242,26 @@ subroutine InitSnowLayers (bounds, snow_depth)
zi => col%zi & ! Output: [real(r8) (:,:) ] interface level below a "z" level (m) (-nlevsno+0:nlevgrnd)
)

! These three variables determine the vertical structure of the snow pack:
! dzmin: minimum of top snow layer
! dzmax_l: maximum thickness of layer when no layers beneath
! dzmax_u: maximum thickness of layer when layers beneath
dzmin(1) = dzmin_1 ! set in the namelist with default or user-defined value
dzmax_l(1) = dzmax_l_1 ! same comment
dzmax_u(1) = dzmax_u_1 ! same comment
do j = 2, nlevsno
dzmin(j) = dzmax_u(j-1) * 0.5_r8
if (dzmin(j) <= dzmin(j-1)) then
dzmin(j) = (dzmin(j-1) + dzmax_u(j-1)) * 0.5_r8
end if
dzmax_u(j) = 2._r8 * dzmax_u(j-1) + 0.01_r8
if (j == 2) then
dzmax_l(j) = dzmax_u(j) + dzmax_l(j-1) - 0.01_r8
else
dzmax_l(j) = dzmax_u(j) + dzmax_l(j-1)
end if
end do

loop_columns: do c = bounds%begc,bounds%endc
l = col%landunit(c)

Expand Down
8 changes: 7 additions & 1 deletion src/main/clm_varcon.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,6 +213,9 @@ module clm_varcon
! The values for the following arrays are set in routine iniTimeConst
!------------------------------------------------------------------

real(r8), public, allocatable :: dzmin(:) !min of top snow layer
real(r8), public, allocatable :: dzmax_l(:) !max snow thickness of layer when no layers beneath
real(r8), public, allocatable :: dzmax_u(:) !max snow thickness of layer when layers beneath
real(r8), public, allocatable :: zlak(:) !lake z (layers)
real(r8), public, allocatable :: dzlak(:) !lake dz (thickness)
real(r8), public, allocatable :: zsoi(:) !soil z (layers)
Expand Down Expand Up @@ -266,7 +269,7 @@ subroutine clm_varcon_init( is_simple_buildtemp )
! MUST be called after clm_varpar_init.
!
! !USES:
use clm_varpar, only: nlevgrnd, nlevlak, nlevdecomp_full, nlevsoifl, nlayer
use clm_varpar, only: nlevsno, nlevgrnd, nlevlak, nlevdecomp_full, nlevsoifl, nlayer
!
! !ARGUMENTS:
implicit none
Expand All @@ -276,6 +279,9 @@ subroutine clm_varcon_init( is_simple_buildtemp )
! Created by E. Kluzek
!------------------------------------------------------------------------------

allocate( dzmin(1:nlevsno ))
allocate( dzmax_l(1:nlevsno ))
allocate( dzmax_u(1:nlevsno ))
allocate( zlak(1:nlevlak ))
allocate( dzlak(1:nlevlak ))
allocate( zsoi(1:nlevgrnd ))
Expand Down

0 comments on commit 3e5856a

Please sign in to comment.