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

Fix interp. of surfdat soil layers so we can use interp. for 10SL case #771

Merged
merged 7 commits into from
Aug 2, 2019
Merged

Fix interp. of surfdat soil layers so we can use interp. for 10SL case #771

Changes from 2 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
d885079
Fix interp. of surfdat soil layers so we can use interp. for 10SL case
billsacks Jul 31, 2019
40a5881
Add some comments
billsacks Aug 1, 2019
9986578
Change organic_frac_squared to private
billsacks Aug 1, 2019
aadb3c5
Record merge of commit 'f40533cd' into soil_layer_def_cleanup_fixinterp
billsacks Aug 1, 2019
3cf9ae7
Merge branch 'master' into soil_layer_def_cleanup_fixinterp
billsacks Aug 1, 2019
d8f31b3
Request close to max wallclock time on cheyenne
billsacks Aug 1, 2019
789554e
Update ChangeLog
billsacks Aug 2, 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
80 changes: 37 additions & 43 deletions src/biogeophys/SoilStateInitTimeConstMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -148,7 +148,6 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
type(file_desc_t) :: ncid ! netcdf id
real(r8) ,pointer :: zsoifl (:) ! Output: [real(r8) (:)] original soil midpoint
real(r8) ,pointer :: zisoifl (:) ! Output: [real(r8) (:)] original soil interface depth
real(r8) ,pointer :: dzsoifl (:) ! Output: [real(r8) (:)] original soil thickness
real(r8) ,pointer :: gti (:) ! read in - fmax
real(r8) ,pointer :: sand3d (:,:) ! read in - soil texture: percent sand (needs to be a pointer for use in ncdio)
real(r8) ,pointer :: clay3d (:,:) ! read in - soil texture: percent clay (needs to be a pointer for use in ncdio)
Expand Down Expand Up @@ -303,22 +302,22 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
! get original soil depths to be used in interpolation of sand and clay
! --------------------------------------------------------------------

allocate(zsoifl(1:nlevsoifl), zisoifl(0:nlevsoifl), dzsoifl(1:nlevsoifl))
do j = 1, nlevsoifl
zsoifl(j) = 0.025*(exp(0.5_r8*(j-0.5_r8))-1._r8) !node depths
enddo

dzsoifl(1) = 0.5_r8*(zsoifl(1)+zsoifl(2)) !thickness b/n two interfaces
do j = 2,nlevsoifl-1
dzsoifl(j)= 0.5_r8*(zsoifl(j+1)-zsoifl(j-1))
! Note that the depths on the file are assumed to be the same as the depths in the
! model when running with 10SL_3.5m. Ideally zsoifl and zisoifl would be read from
! the surface dataset rather than assumed here.
!
! We need to specify zsoifl down to nlevsoifl+1 (rather than just nlevsoifl) so that
! we can get the appropriate zisoifl at level nlevsoifl (i.e., the bottom interface
! depth).
allocate(zsoifl(1:nlevsoifl+1), zisoifl(0:nlevsoifl))
do j = 1, nlevsoifl+1
zsoifl(j) = 0.025_r8*(exp(0.5_r8*(j-0.5_r8))-1._r8) !node depths
Copy link
Member Author

Choose a reason for hiding this comment

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

Diff (1) from my main PR comment: addition of _r8

Copy link
Contributor

Choose a reason for hiding this comment

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

Why do we calculate these depths if they are on the file and could be read? Seems risky to be calculating depths that are assumed to be on a file, when the file may change and we could be calculating the wrong depths at that point.

Copy link
Member Author

Choose a reason for hiding this comment

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

Yup, I completely agree that these should be on the file. But they aren't on there right now. Long-term, the best solution would absolutely be to put these levels on the surface datasets, but I felt that was too much to fold in with #759 or this PR.

enddo
dzsoifl(nlevsoifl) = zsoifl(nlevsoifl)-zsoifl(nlevsoifl-1)

zisoifl(0) = 0._r8
do j = 1, nlevsoifl-1
do j = 1, nlevsoifl
Copy link
Member Author

Choose a reason for hiding this comment

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

Diff (3) from my main PR comment (and related changes around here): changing definition of zisoifl for lowest level

Copy link
Contributor

Choose a reason for hiding this comment

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

Considering zsoifl represents the depths in the soil texture file, we introduce an abstraction here by defining a zsoifl for a soil layer that doesn't really exist, i.e. nlevsoifl + 1. Do you feel comfortable with that?

Copy link
Member Author

Choose a reason for hiding this comment

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

I agree that this feels weird, and certainly I'd feel less comfortable with this if zsoifl were used for anything more than defining zisoifl. But since zsoifl is just used to define zisoifl, I figured that the best thing to do here is to define zsoifl so that the zisoifl values end up as desired - which involves setting an assumed value for zsoifl(nlevsoifl+1). Unless you see a cleaner way to accomplish this goal?

Copy link
Contributor

@slevis-lmwg slevis-lmwg Aug 1, 2019 •
edited
Loading

Choose a reason for hiding this comment

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

Sounds good. To help with the long term development and institutional memory of what is going on here, how about a comment that explains how this code might change if we were reading the depths from the file (a file); and which specific file we are referring to; and the danger of switching to a new file given the current implementation. If such a comment exists somewhere else, I suggest copying it here, too.

Copy link
Member Author

Choose a reason for hiding this comment

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

Comment added in 40a5881

zisoifl(j) = 0.5_r8*(zsoifl(j)+zsoifl(j+1)) !interface depths
enddo
zisoifl(nlevsoifl) = zsoifl(nlevsoifl) + 0.5_r8*dzsoifl(nlevsoifl)

! --------------------------------------------------------------------
! Set soil hydraulic and thermal properties: non-lake
Expand All @@ -330,7 +329,6 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
! roof, sunwall and shadewall are prescribed in SoilThermProp.F90
! in SoilPhysicsMod.F90


do c = begc, endc
g = col%gridcell(c)
l = col%landunit(c)
Expand Down Expand Up @@ -367,36 +365,32 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
else

do lev = 1,nlevgrnd
! DML - this if statement could probably be removed and just the
! top part used for all soil layer structures
if (.not. organic_frac_squared) then
if (lev .eq. 1) then
clay = clay3d(g,1)
sand = sand3d(g,1)
om_frac = organic3d(g,1)/organic_max
else if (lev <= nlevsoi) then
do j = 1,nlevsoifl-1
if (zisoi(lev) >= zisoifl(j) .AND. zisoi(lev) < zisoifl(j+1)) then
clay = clay3d(g,j+1)
sand = sand3d(g,j+1)
om_frac = organic3d(g,j+1)/organic_max
endif
end do
else
clay = clay3d(g,nlevsoifl)
sand = sand3d(g,nlevsoifl)
om_frac = 0._r8
endif
if (lev .eq. 1) then
clay = clay3d(g,1)
sand = sand3d(g,1)
om_frac = organic3d(g,1)/organic_max
else if (lev <= nlevsoi) then
do j = 1,nlevsoifl-1
! NOTE(wjs, 2019-08-01) It appears that the code currently doesn't set
! clay, sand and om_frac explicitly under some conditions. It probably
! should. My understanding is that currently things work okay, though,
! because clay, sand and om_frac will remain set at their previous
! values, which is probably reasonable enough. See also
! <https://github.com/ESCOMP/ctsm/pull/771#discussion_r309509596>.
if (zisoi(lev) > zisoifl(j) .AND. zisoi(lev) <= zisoifl(j+1)) then
Copy link
Member Author

Choose a reason for hiding this comment

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

Diff (2) from my main PR comment: changing the logic for when zisoi == zisoifl

Copy link
Contributor

Choose a reason for hiding this comment

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

What happens here for, let's say, lev = 20 if zisoi(20) happens to be greater than zisoifl(10). I think that clay, sand, and om_frac do not get valid values for lev = 20 in that case. Even if somehow they do get valid values, maybe we need an explicit statement, so that there is no doubt about their values. I'm thinking something like this:
else if (zisoi(lev) > zisoifl(j+1)) then
clay = clay3d(g,j+1)
sand = sand3d(g,j+1)
om_frac = organic3d(g,j+1)/organic_max

Copy link
Contributor

Choose a reason for hiding this comment

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

Oh, and what if zisoi(lev) happens to be less than or equal to zisoifl(j). That case is not accounted for, either, I think.

Copy link
Member Author

Choose a reason for hiding this comment

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

That's a good point. I believe these problems existed before this PR; do you agree? (I'm just trying to distinguish between possible new problems in this PR and pre-existing problems.)

For this

Oh, and what if zisoi(lev) happens to be less than or equal to zisoifl(j). That case is not accounted for, either, I think.

I'm thinking you mean: zisoi(lev) is <= zisoifl(1), for lev = 2 or more.

Now, that said: I believe the code is actually working in a reasonable way right now, it just isn't as explicit as it could be: if these variables (clay, sand, om_frac) aren't explicitly set, they'll remain at their values from the last iteration of the level loop. So, for the case where zisoi(lev) is <= zisoifl(1), I believe we'll end up using values from level 1 of the file (because these will have been set in the lev .eq. 1 block); this seems reasonable. And, for the first case you raise, I believe we'll end up using values equal to the those in the model layer above this one; this may not be the right thing to do (e.g., if the last model layer took values from file layer 8: I think we'd really want to use file layer 10 in this case), but it doesn't seem terrible.

We could fix this to be more explicit; if so, I think this code would become cleaner - and unit testable! - if we extracted this logic for finding the right layer into its own routine, tied in with the parenthetical under "Caveat c" in this comment #759 (comment) .

However, it might make the most sense to wait and do this as part of resolving #772 , once we agree on how we want to solve that one - rather than fixing this implementation just to throw it out and replace it with something more correct.

What do you think?

Copy link
Contributor

Choose a reason for hiding this comment

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

I'm fine with that. Again though, how about adding a comment in the code with your best understanding of how this works now versus how it should work ideally.

Copy link
Contributor

Choose a reason for hiding this comment

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

moving the logic matching soil layers in the model with soil layers on the file to outside the column loop, since this matching is independent of column

Regarding the parenthetical under caveat c in #759, if there's a chance that some day we will end up with varying soil layer depths, it may not be worth changing. I'm just raising as a potential concern; @swensosc and others may have better insight.

Copy link
Member Author

Choose a reason for hiding this comment

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

Comment added in 40a5881

clay = clay3d(g,j+1)
sand = sand3d(g,j+1)
om_frac = organic3d(g,j+1)/organic_max
endif
end do
else
if (lev <= nlevsoi) then ! duplicate clay and sand values from 10th soil layer
clay = clay3d(g,lev)
sand = sand3d(g,lev)
om_frac = (organic3d(g,lev)/organic_max)**2._r8
else
clay = clay3d(g,nlevsoi)
sand = sand3d(g,nlevsoi)
om_frac = 0._r8
endif
clay = clay3d(g,nlevsoifl)
sand = sand3d(g,nlevsoifl)
om_frac = 0._r8
endif

if (organic_frac_squared) then
om_frac = om_frac**2._r8
end if

if (lun%urbpoi(l)) then
Expand Down Expand Up @@ -591,7 +585,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
! --------------------------------------------------------------------

deallocate(sand3d, clay3d, organic3d)
deallocate(zisoifl, zsoifl, dzsoifl)
deallocate(zisoifl, zsoifl)

end subroutine SoilStateInitTimeConst

Expand Down