Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

User defined top snow layer #792

Merged
merged 23 commits into from
Sep 9, 2019
Merged

User defined top snow layer #792

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
63b2c81
Changes from @olyson to fix #780
ekluzek Aug 16, 2019
20cbf74
Merge tag 'ctsm1.0.dev059' into snow_layers_user_defined
slevis-lmwg Aug 25, 2019
3e5856a
Allow users to specify snow layer structure in namelist
slevis-lmwg Aug 26, 2019
986cf49
Merge commit 'f2869fe202f71c4fd8c854fce47089630b4396c5' into clm50cos…
ekluzek Aug 28, 2019
7f67efc
Remove the commit to update the params file, which shouldn't be done …
ekluzek Aug 28, 2019
ee939a7
Merge tag 'ctsm1.0.dev061' into clm50coslatfirefix
ekluzek Sep 2, 2019
b2cf3a6
Changes for #787 for dwt_slash_flux
ekluzek Sep 2, 2019
f293334
Make sure urbantv is set for clm4_5 #175
ekluzek Sep 2, 2019
f49dac1
Revisions in response to @billsacks code review
slevis-lmwg Sep 2, 2019
47cd874
Merge tag 'ctsm1.0.dev061' into user_defined_top_snow_layer
slevis-lmwg Sep 2, 2019
ae415da
Merge tag 'ctsm1.0.dev062' into clm50coslatfirefix
ekluzek Sep 4, 2019
d31b2ee
Fix urbantv to point to different file for clm4_5 than clm5_0 as it s…
ekluzek Sep 4, 2019
ffe027c
Remove DWT_SLASH_* from output usermod at scales finer than gridcell,…
ekluzek Sep 6, 2019
ec4d91a
Update changelog
ekluzek Sep 6, 2019
5bcf4d6
Merge tag 'ctsm1.0.dev063' into user_defined_top_snow_layer
slevis-lmwg Sep 7, 2019
e00da7c
Revisions in response to code review
slevis-lmwg Sep 7, 2019
d7292f3
First draft of ChangeLog entry
slevis-lmwg Sep 7, 2019
e478244
Updated temporary test code checking bit-for-bit
slevis-lmwg Sep 7, 2019
db4bb80
Updated temporary test code checking bit-for-bit (part 2)
slevis-lmwg Sep 7, 2019
97cef78
Removing temporary test code
slevis-lmwg Sep 7, 2019
bd7f647
Updates needed for snowhydrology unit tests to pass
slevis-lmwg Sep 9, 2019
4c19afc
Updated ChangeLog/Sum
slevis-lmwg Sep 9, 2019
1d5ea12
Minor edit in ChangeLog
billsacks Sep 9, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 42 additions & 0 deletions bld/CLMBuildNamelist.pm
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3521,6 +3521,48 @@ sub setup_logic_snowpack {
add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'reset_snow');
add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'reset_snow_glc');
add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'reset_snow_glc_ela');
add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'snow_dzmin_1');
add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'snow_dzmin_2');
add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'snow_dzmax_l_1');
add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'snow_dzmax_l_2');
add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'snow_dzmax_u_1');
add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'snow_dzmax_u_2');

my $dzmin1 = $nl->get_value('snow_dzmin_1');
my $dzmin2 = $nl->get_value('snow_dzmin_2');
my $dzmax_l1 = $nl->get_value('snow_dzmax_l_1');
my $dzmax_l2 = $nl->get_value('snow_dzmax_l_2');
my $dzmax_u1 = $nl->get_value('snow_dzmax_u_1');
my $dzmax_u2 = $nl->get_value('snow_dzmax_u_2');

if ($dzmin1 != 0.01 || $dzmin2 != 0.015 || $dzmax_u1 != 0.02 || $dzmax_u2 != 0.05 || $dzmax_l1 != 0.03 || $dzmax_l2 != 0.07) {
$log->warning("Setting any of the following namelist variables to NON DEFAULT values remains untested as of Sep 1, 2019: snow_dzmin_1 & 2, snow_dzmax_u_1 & 2, snow_dzmax_l_1 & 2." );
$log->warning("Leave these variables unspecified in user_nl_clm in order to use the default values." );
}
if ($dzmin1 <= 0.0 || $dzmin2 <= 0.0 || $dzmax_u1 <= 0.0 || $dzmax_u2 <= 0.0 || $dzmax_l1 <= 0.0 || $dzmax_l2 <= 0.0) {
$log->fatal_error('One or more of the snow_dzmin_* and/or snow_dzmax_* were set incorrectly to be <= 0');
}
if ($dzmin2 <= $dzmin1) {
$log->fatal_error('snow_dzmin_2 was set incorrectly to be <= snow_dzmin_1');
}
if ($dzmax_l2 <= $dzmax_l1) {
$log->fatal_error('snow_dzmax_l_2 was set incorrectly to be <= snow_dzmax_l_1');
}
if ($dzmax_u2 <= $dzmax_u1) {
$log->fatal_error('snow_dzmax_u_2 was set incorrectly to be <= snow_dzmax_u_1');
}
if ($dzmin1 >= $dzmax_u1) {
$log->fatal_error('snow_dzmin_1 was set incorrectly to be >= snow_dzmax_u_1');
}
if ($dzmin2 >= $dzmax_u2) {
$log->fatal_error('snow_dzmin_2 was set incorrectly to be >= snow_dzmax_u_2');
}
if ($dzmax_u1 >= $dzmax_l1) {
$log->fatal_error('snow_dzmax_u_1 was set incorrectly to be >= snow_dzmax_l_1');
}
if ($dzmax_u2 >= $dzmax_l2) {
$log->fatal_error('snow_dzmax_u_2 was set incorrectly to be >= snow_dzmax_l_2');
}

if (remove_leading_and_trailing_quotes($nl->get_value('snow_overburden_compaction_method')) eq 'Vionnet2012') {
# overburden_compress_tfactor isn't used if we're using the Vionnet2012
Expand Down
7 changes: 7 additions & 0 deletions bld/namelist_files/namelist_defaults_ctsm.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -298,6 +298,13 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
<h2osno_max phys="clm5_0" structure="fast" >5000.0</h2osno_max>
<h2osno_max phys="clm4_5" >1000.0</h2osno_max>

<snow_dzmin_1>0.010d00</snow_dzmin_1>
<snow_dzmin_2>0.015d00</snow_dzmin_2>
<snow_dzmax_l_1>0.03d00</snow_dzmax_l_1>
<snow_dzmax_l_2>0.07d00</snow_dzmax_l_2>
<snow_dzmax_u_1>0.02d00</snow_dzmax_u_1>
<snow_dzmax_u_2>0.05d00</snow_dzmax_u_2>

<int_snow_max phys="clm5_0">2000.</int_snow_max>
<!-- For clm4_5, make this effectively unlimited -->
<int_snow_max phys="clm4_5">1.e30</int_snow_max>
Expand Down
27 changes: 27 additions & 0 deletions bld/namelist_files/namelist_definition_ctsm.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2440,6 +2440,33 @@ Changes in this value should possibly be accompanied by changes in:
glc_snow_persistence_max_days > 0.
</entry>

<entry id="snow_dzmin_1" type="real" category="clm_physics"
group="clm_snowhydrology_inparm" valid_values="" >
Min snow thickness of layer 1
</entry>
<entry id="snow_dzmin_2" type="real" category="clm_physics"
group="clm_snowhydrology_inparm" valid_values="" >
Min snow thickness of layer 2
</entry>

<entry id="snow_dzmax_l_1" type="real" category="clm_physics"
group="clm_snowhydrology_inparm" valid_values="" >
Max snow thickness of layer 1 when no layers beneath
</entry>
<entry id="snow_dzmax_l_2" type="real" category="clm_physics"
group="clm_snowhydrology_inparm" valid_values="" >
Max snow thickness of layer 2 when no layers beneath
</entry>

<entry id="snow_dzmax_u_1" type="real" category="clm_physics"
group="clm_snowhydrology_inparm" valid_values="" >
Max snow thickness of layer 1 when layers beneath
</entry>
<entry id="snow_dzmax_u_2" type="real" category="clm_physics"
group="clm_snowhydrology_inparm" valid_values="" >
Max snow thickness of layer 2 when layers beneath
</entry>

<entry id="int_snow_max" type="real" category="clm_physics"
group="scf_swenson_lawrence_2012_inparm" valid_values="" >
Limit applied to integrated snowfall when determining changes in snow-covered fraction during melt
Expand Down
144 changes: 124 additions & 20 deletions src/biogeophys/SnowHydrologyMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -101,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) /)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In reference to this comment, the original code assigned huge(1._r8) to dzmax_l(12) and dzmax_u(12). So to match the original behavior, I think I need to do the following:

if (nlevsno >= 12) then
dzmax_l = huge(1._r8)
dzmax_u = huge(1._r8)
end if

Is this what you were suggesting @billsacks ?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

More precisely...
if (j >= 12) then
dzmax_l(j) = ...
dzmax_u(j) = ...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Changing mind:
I prefer the first if statement with these corrections:
dzmax_l(nlevsno) = ...
dzmax_u(nlevsno) = ...

This replicates the current behavior for nlevsno = 12. If some day we allow nlevsno > 12, this just sets the bottom layer to huge.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, the last comment is what I was thinking (where the ... should be huge(1._r8)).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i.e., the bottom layer should be set to huge, regardless of the number of snow layers. I now realize that that isn't exactly what the original code did, but I think that was the intent.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess it would be good to verify that that does the right thing. Either (a) make that change and verify that it doesn't change answers for, say, a 5-layer case (a CLM45 case), or (b) don't make that change, and ensure that having dzmax_l and dzmax_u NOT set to huge for nlevsno gives the same answers as setting it to huge for a 12-layer case. I prefer (a), but don't feel strongly about it.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I also prefer (a).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess it would be good to verify that that does the right thing. Either (a) make that change and verify that it doesn't change answers for, say, a 5-layer case (a CLM45 case), or [...]

The hobart test-suite includes a few clm45 cases and they all PASS.


!
! !PRIVATE DATA MEMBERS:

Expand All @@ -124,6 +110,17 @@ 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) :: snow_dzmin_1, snow_dzmax_l_1, snow_dzmax_u_1 ! namelist-defined top snow layer information
real(r8) :: snow_dzmin_2, snow_dzmax_l_2, snow_dzmax_u_2 ! namelist-defined 2nd snow layer information
real(r8), private, allocatable :: dzmin(:) !min snow thickness of layer
real(r8), private, allocatable :: dzmax_l(:) !max snow thickness of layer when no layers beneath
real(r8), private, allocatable :: dzmax_u(:) !max snow thickness of layer when layers beneath

! FOR TEMPORARY TESTING. DELETE WHEN DONE. slevis diag
real(r8), private, allocatable :: dzmin_orig(:)
real(r8), private, allocatable :: dzmax_l_orig(:)
real(r8), private, allocatable :: dzmax_u_orig(:)

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 @@ -171,6 +168,7 @@ subroutine SnowHydrology_readnl( NLFilename)
use shr_nl_mod , only : shr_nl_find_group_name
use spmdMod , only : masterproc, mpicom
use shr_mpi_mod , only : shr_mpi_bcast
use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=)
!
! !ARGUMENTS:
character(len=*), intent(in) :: NLFilename ! Namelist filename
Expand All @@ -188,11 +186,19 @@ 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, &
snow_dzmin_1, snow_dzmax_l_1, snow_dzmax_u_1, &
snow_dzmin_2, snow_dzmax_l_2, snow_dzmax_u_2

! Initialize options to default values, in case they are not specified in the namelist
wind_dependent_snow_density = .false.
snow_overburden_compaction_method = ' '
snow_dzmin_1 = nan
snow_dzmin_2 = nan
snow_dzmax_l_1 = nan
snow_dzmax_l_2 = nan
snow_dzmax_u_1 = nan
snow_dzmax_u_2 = nan

if (masterproc) then
unitn = getavu()
Expand All @@ -219,6 +225,12 @@ 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 (snow_dzmin_1, mpicom)
call shr_mpi_bcast (snow_dzmin_2, mpicom)
call shr_mpi_bcast (snow_dzmax_l_1, mpicom)
call shr_mpi_bcast (snow_dzmax_l_2, mpicom)
call shr_mpi_bcast (snow_dzmax_u_1, mpicom)
call shr_mpi_bcast (snow_dzmax_u_2, mpicom)

if (masterproc) then
write(iulog,*) ' '
Expand Down Expand Up @@ -801,12 +813,12 @@ subroutine BuildFilter_SnowpackInitialized(bounds, num_c, filter_c, &
! maintain answers as before.
snowpack_initialized(c) = ( &
snl(c) == 0 .and. &
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@billsacks I think the next four 0.01 values are what you wanted me to look for, right?

More generally speaking, I think this code is ready for review. Also I would like to know your preference in this case between new Unit Test and new standard test-suite test. I will be working on that next.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, these are the four places where I know this value is hard-coded. If you haven't already, thenI also thought it would be worth grepping the code for 0.01 to see if there are any other references that I'm not aware of.

As for testing: With some upcoming suggestions that I'm going to make (in which you don't have hard-coded defaults for these values, but rather have them always explicitly set in the namelist), your new code will be sufficiently covered by existing system tests, so no need to add additional tests.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I looked for all the instances of 0.01 at the time, and there were some more that looked unrelated, so I left them alone.

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 @@ -1541,7 +1553,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 thickness of snow layer (local)
real(r8):: dtime !land model time step (sec)

!-----------------------------------------------------------------------
Expand Down Expand Up @@ -1725,8 +1737,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 @@ -2259,6 +2271,7 @@ subroutine InitSnowLayers (bounds, snow_depth)
!
! !USES:
use clm_varcon , only : spval
use spmdMod, only : masterproc
!
! !ARGUMENTS:
type(bounds_type) , intent(in) :: bounds
Expand All @@ -2279,6 +2292,97 @@ subroutine InitSnowLayers (bounds, snow_depth)
zi => col%zi & ! Output: [real(r8) (:,:) ] interface level below a "z" level (m) (-nlevsno+0:nlevgrnd)
)

allocate(dzmin(1:nlevsno))
allocate(dzmax_l(1:nlevsno))
allocate(dzmax_u(1:nlevsno))

! FOR TEMPORARY TESTING. DELETE WHEN DONE. slevis diag
allocate(dzmin_orig(1:nlevsno))
allocate(dzmax_l_orig(1:nlevsno))
allocate(dzmax_u_orig(1:nlevsno))

! These three variables determine the vertical structure of the snow pack:
! dzmin: minimum snow thickness of layer
! dzmax_l: maximum snow thickness of layer when no layers beneath
! dzmax_u: maximum snow thickness of layer when layers beneath
dzmin(1) = snow_dzmin_1 ! default or user-defined value from namelist
dzmax_l(1) = snow_dzmax_l_1 ! same comment
dzmax_u(1) = snow_dzmax_u_1 ! same comment
dzmin(2) = snow_dzmin_2 ! default or user-defined value from namelist
dzmax_l(2) = snow_dzmax_l_2 ! same comment
dzmax_u(2) = snow_dzmax_u_2 ! same comment
do j = 3, nlevsno
dzmin(j) = dzmax_u(j-1) * 0.5_r8
dzmax_u(j) = 2._r8 * dzmax_u(j-1) + 0.01_r8
dzmax_l(j) = dzmax_u(j) + dzmax_l(j-1)
! error checks
if (dzmin(j) <= dzmin(j-1)) then
write(iulog,*) 'ERROR at snow layer j =', j, ' because dzmin(j) =', dzmin(j), ' and dzmin(j-1) =', dzmin(j-1)
call endrun(msg="ERROR dzmin(j) cannot be <= dzmin(j-1)"// &
errMsg(sourcefile, __LINE__))
end if
if (dzmax_u(j) <= dzmax_u(j-1)) then
write(iulog,*) 'ERROR at snow layer j =', j, ' because dzmax_u(j) =', dzmax_u(j), ' and dzmax_u(j-1) =', dzmax_u(j-1)
call endrun(msg="ERROR dzmax_u(j) cannot be <= dzmax_u(j-1)"// &
errMsg(sourcefile, __LINE__))
end if
if (dzmax_l(j) <= dzmax_l(j-1)) then
write(iulog,*) 'ERROR at snow layer j =', j, ' because dzmax_l(j) =', dzmax_l(j), ' and dzmax_l(j-1) =', dzmax_l(j-1)
call endrun(msg="ERROR dzmax_l(j) cannot be <= dzmax_l(j-1)"// &
errMsg(sourcefile, __LINE__))
end if
if (dzmin(j) >= dzmax_u(j)) then
write(iulog,*) 'ERROR at snow layer j =', j, ' because dzmin(j) =', dzmin(j), ' and dzmax_u(j) =', dzmax_u(j)
call endrun(msg="ERROR dzmin(j) cannot be >= dzmax_u(j)"// &
errMsg(sourcefile, __LINE__))
end if
if (dzmax_u(j) >= dzmax_l(j)) then
write(iulog,*) 'ERROR at snow layer j =', j, ' because dzmax_u(j) =', dzmax_u(j), ' and dzmax_l(j) =', dzmax_l(j)
call endrun(msg="ERROR dzmax_u(j) cannot be >= dzmax_l(j)"// &
errMsg(sourcefile, __LINE__))
end if
end do
if (nlevsno >= 12) then ! max nlevsno is 12
dzmax_u(nlevsno) = huge(1._r8) ! the original code set layer 12 to huge
dzmax_l(nlevsno) = huge(1._r8) ! same comment
end if

if (masterproc) then
write(iulog,*) 'dzmin =', dzmin
write(iulog,*) 'dzmax_l =', dzmax_l
write(iulog,*) 'dzmax_u =', dzmax_u
end if

! BEGIN TEMPORARY TESTING. DELETE WHEN DONE. slevis diag
dzmin_orig = (/ 0.01_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 /)
if (maxval(abs(dzmin(1:nlevsno) - dzmin_orig(1:nlevsno))) > 1.e-13_r8) then
write(iulog,*) 'dzmin_orig =', dzmin_orig
call endrun(msg="ERROR dzmin_orig /= dzmin"// &
errMsg(sourcefile, __LINE__))
else
dzmin(1:nlevsno) = dzmin_orig(1:nlevsno)
end if
dzmax_u_orig = (/ 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) /)
if (maxval(abs(dzmax_u(1:nlevsno) - dzmax_u_orig(1:nlevsno))) > 1.e-13_r8) then
write(iulog,*) 'dzmax_u_orig =', dzmax_u_orig
call endrun(msg="ERROR dzmax_u_orig /= dzmax_u"// &
errMsg(sourcefile, __LINE__))
else
dzmax_u(1:nlevsno) = dzmax_u_orig(1:nlevsno)
end if
dzmax_l_orig = (/ 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.30_r8, huge(1._r8) /)
if (maxval(abs(dzmax_l(1:nlevsno) - dzmax_l_orig(1:nlevsno))) > 1.e-13_r8) then
write(iulog,*) 'dzmax_l_orig =', dzmax_l_orig
call endrun(msg="ERROR dzmax_l_orig /= dzmax_l"// &
errMsg(sourcefile, __LINE__))
else
dzmax_l(1:nlevsno) = dzmax_l_orig(1:nlevsno)
end if
! END TEMPORARY TESTING slevis diag

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

Expand Down