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.

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

Dynamic lakes - without tools changes #1109

Merged
merged 47 commits into from
Sep 29, 2020

Conversation

billsacks
Copy link
Member

@billsacks billsacks commented Aug 17, 2020

Description of changes

Dynamic lakes changes from @Ivanderkelen . I have split her initial branch into two separate PRs: one containing just the tools changes (#1073) and one with everything EXCEPT the tools changes (this one).

Specific notes

Contributors other than yourself, if any: @Ivanderkelen

CTSM Issues Fixed (include github issue #):
Resolves #200
Resolves #1140

Are answers expected to change (and if so in what way)? No

Any User Interface Changes (namelist or namelist defaults changes)?

Testing performed, if any:
(List what testing you did to show your changes worked as expected)
(This can be manual testing or running of the different test suites)
(Documentation on system testing is here: https://github.com/ESCOMP/ctsm/wiki/System-Testing-Guide)
(aux_clm on cheyenne for intel/gnu and izumi for intel/gnu/nag/pgi is the standard for tags on master)

NOTE: Be sure to check your coding style against the standard
(https://github.com/ESCOMP/ctsm/wiki/CTSM-coding-guidelines) and review
the list of common problems to watch out for
(https://github.com/ESCOMP/CTSM/wiki/List-of-common-problems).

Ivanderkelen and others added 21 commits August 17, 2020 14:52
…dynConsBiogeophysMod.f90

cv is kept constant at water cv for simplicity (altough there is a cv for each lake layer)
TEMPORARY change
… as an history field.

The history field only contains the heat content of the lake water itself, and is a column output.
…namic lake land unit. Use similar to do_transient_crops
(This commit was split off from 60c0aaf using interactive rebase,
discarding the tools changes and just keeping the relatively small
source code changes in that commit.)
This is a double precision variable on the file, so I'm making it double
precision in the code. I think ncdio_pio used to handle data type
conversions like this, but now it seems to be corrupting memory: the
original code led to errors in decompInitMod, due to what appears to be
corruption in some unrelated variable.
I'm treating do_transient_lakes similar to do_transient_pfts and
do_transient_crops in this respect. I'm not sure if all of these checks
are truly important for transient lakes (in particular, my guess is that
collapse_urban could probably be done with transient lakes - as well as
transient pfts and transient crops for that matter), but some of the
checks probably are needed, and it seems best to keep transient lakes
consistent with other transient areas in this respect.
The new function setup_logic_do_transient_lakes is based on
setup_logic_do_transient_crops.

As with my previous commit: I'm not sure if all of these checks are
truly important for transient lakes (in particular, my guess is that
collapse_urban could probably be done with transient lakes - as well as
transient pfts and transient crops for that matter), but some of the
checks probably are needed, and it seems best to keep transient lakes
consistent with other transient areas in this respect.
I didn't see any reason why the default value should be set in the code
rather than coming from namelist_defaults_ctsm.xml. Getting it from the
latter will make it easier for the default value to differ with
different physics versions.
@billsacks billsacks added the PR status: ready PR: this is ready to merge in, with all tests satisfactory and reviews complete label Aug 17, 2020
@billsacks billsacks self-assigned this Aug 17, 2020
@billsacks
Copy link
Member Author

I created this new branch by doing an interactive rebase and picking only the commits that changed files outside of the tools directory. Nearly all commits touched just tools or just files outside of tools (thank you @Ivanderkelen for keeping those changes separate!); for the one commit that touched both, I edited the commit so that only the non-tools changes ended up here. I then backed out all non-tools changes in the original branch (using git checkout to get the version of the src and tools directory from ctsm1.0.dev101) (see Ivanderkelen#2).

I have confirmed that a merge of these two separate branches results in an identical set of code to the original changes (i.e., merging 7fb5ce9 into 5a64ed8 is identical to 29af7cf).

Copy link
Member Author

@billsacks billsacks left a comment

Choose a reason for hiding this comment

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

@Ivanderkelen this is looking close to ready! I have a few inline comments below. Note that this branch (dynlakes_master_notools) is on my fork. What I'd suggest for a workflow is that you pull this down from my fork, make any changes, push it back to your fork, then open a PR into the branch on my fork. Something like:

git remote add billsacks https://github.com/billsacks/ctsm.git
git fetch billsacks
git checkout -b dynlakes_master_notools billsacks/dynlakes_master_notools
# Make any changes and commit them
git push -u Ivanderkelen dynlakes_master_notools

Then open a PR, selecting my fork and dynlakes_master_notools as the base.

Let me know if you have any questions about any of this or would like my help with any of it.

I'm going to work on getting a single-point test in place that exercises the new dynamic lake code. Once that is in place and the below comments are addressed, I think this is ready for final testing.

Comment on lines 1007 to 1010
! write(iulog,*)'Energy content of lake after calculating lake temperature (J/m²)', ncvts

! IV: currently commented out: caused crash. To do: look at this part of the code!!!
! lake_heat(c) = ncvts(c)
Copy link
Member Author

@billsacks billsacks Aug 23, 2020

Choose a reason for hiding this comment

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

  • @Ivanderkelen is there anything you still wanted to look at here, or can we remove these commented out lines? Actually, it looks like currently temperature_inst%lake_heat is never set, so if we leave this out, it looks like we should also remove temperature_inst%lake_heat.

Copy link
Contributor

@Ivanderkelen Ivanderkelen Aug 24, 2020

Choose a reason for hiding this comment

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

You're right: we can safely remove those lines, and the temperature_inst%lake_heat.

I wrote those lines to perform a double check on lake heat, as calculated in the LakeTemperatureMod.F90. As it is now included in the AccumulateHeatLake routine in TotalWaterHeatMod.F90, I removed those comments (e44c19d).

@@ -117,6 +117,9 @@ module TemperatureType
real(r8), pointer :: fact_col (:,:) ! used in computing tridiagonal matrix
real(r8), pointer :: c_h2osfc_col (:) ! heat capacity of surface water

! lake heat
real(r8), pointer :: lake_heat (:) ! total heat of lake water (J/m²)
Copy link
Member Author

@billsacks billsacks Aug 23, 2020

Choose a reason for hiding this comment

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

  • As I note in LakeTemperatureMod, it looks like this currently isn't set... unless that's changed, this lake_heat variable should be removed.

  • If this variable is kept, it should be renamed to lake_heat_col.

Copy link
Contributor

@Ivanderkelen Ivanderkelen Aug 24, 2020

Choose a reason for hiding this comment

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

I agree and removed the variable. For checking I originally also wrote lake_heat to a history field, these parts are now removed too (e44c19d and e37d269, also in TotalWaterAndHeatMod).

end do
end do

! add ice heat and liquid heat together (look at this part)
Copy link
Member Author

@billsacks billsacks Aug 23, 2020

Choose a reason for hiding this comment

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

  • @Ivanderkelen does "look at this part" mean that you still wanted to give this more thought (or that you wanted me to particularly look at it)? If you're satisfied with this, we should remove the comment in parentheses.

Copy link
Contributor

Choose a reason for hiding this comment

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

I am satisfied on this, and removed the "(look at this part)". (2bfb8d7)

It was a remnant to notify myself when doing the actual coding, but I forgot to remove it along the way.

Comment on lines 1138 to 1139
! ADD HERE NOTES ABOUT LAKES!!! lake water not accounted. because baselines
!
Copy link
Member Author

@billsacks billsacks Aug 23, 2020

Choose a reason for hiding this comment

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

  • @Ivanderkelen do you remember what you wanted to add here? If not, I could see if I can dig up any notes.

Copy link
Contributor

Choose a reason for hiding this comment

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

Yes, I think I figured out what I wanted to add:

With dynamical lakes, this routine does not correct for the heat carried within the lake water itself.
When a lake column grows, extra water content is added to the grid cell, but as water content is expressed in kg/m² and lake depth is constant, the size of the lake does not matter. Because of this, after substracting baselines, the change in lake water content will always be zero. I guess could be solved once the temperature of runoff is tracked, and lake depths can change too.

Does this makes senses to you?

I adjusted the description as follows (but feel free to change or update to make it more clear):

Suggested change
! ADD HERE NOTES ABOUT LAKES!!! lake water not accounted. because baselines
!
! With dynamical lakes, the adjusted delta_heat does not account for the added lake
! water content due to growing lakes. This is because lake depth is constant, the
! total lake water content (kg/m^2) does not change. The change in water content of
! the snow and soil in the lake column are accounted for.`

Copy link
Member Author

Choose a reason for hiding this comment

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

Thanks! To be honest, my head isn't enough in this to give this a critical read, but I think this seems good.

Comment on lines 452 to 478


! IV: add here call to subroutine to read in lake mask (necessary for initialisation of dynamical lakes)
! First open landuse.timeseries file for this.

if (masterproc) then
write(iulog,*) 'Attempting to read landuse.timeseries data .....'
if (lfsurdat == ' ') then
write(iulog,*)'fdynuse must be specified'
call endrun(msg=errMsg(sourcefile, __LINE__))
endif
endif


! open landuse_timeseries file
fdynuse = get_flanduse_timeseries()

call getfil(fdynuse, locfn, 0 )

call ncd_pio_openfile (ncid_dynuse, trim(locfn), 0)

! read the lakemask
call surfrd_lakemask(begg, endg, ncid_dynuse)

! close landuse_timeseries file again
call ncd_pio_closefile(ncid_dynuse)

Copy link
Member Author

@billsacks billsacks Aug 23, 2020

Choose a reason for hiding this comment

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

  • surfrd_lakemask (and all of the code you added in this subroutine) should only be done if we're running with transient lakes

  • Also, I think it would be cleanest if you moved all of this new code into surfrd_lakemask

So the changes to surfrd_get_data should just look like:

if (get_do_transient_lakes()) then
   call surfrd_lakemask(begg, endg)
end if

Copy link
Contributor

Choose a reason for hiding this comment

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

I applied the proposed changes and moved the added code to the surfrd_lakemask routine, see 30394f8

Comment on lines 454 to 455
! IV: add here call to subroutine to read in lake mask (necessary for initialisation of dynamical lakes)
! First open landuse.timeseries file for this.
Copy link
Member Author

@billsacks billsacks Aug 23, 2020

Choose a reason for hiding this comment

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

  • This comment should be removed

Copy link
Contributor

Choose a reason for hiding this comment

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

The comment is removed (30394f8)


if (masterproc) then
write(iulog,*) 'Attempting to read landuse.timeseries data .....'
if (lfsurdat == ' ') then
Copy link
Member Author

@billsacks billsacks Aug 23, 2020

Choose a reason for hiding this comment

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

  • This check should be on fdynuse, not lfsurdat. For this to work, you'll need to move this error check code a few lines down, to after fdynuse = get_flanduse_timeseries()

Copy link
Contributor

Choose a reason for hiding this comment

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

This is now updated (30394f8)

@billsacks
Copy link
Member Author

billsacks commented Aug 27, 2020

I've looked a bit more into the dz_lake issue mentioned in my last comment. I noticed this:

cice_eff = cpice*denh2o !use water density because layer depth is not adjusted
!for freezing

  • So I think, for consistency, the new accounting of lake ice in TotalWaterAndHeatMod should also use denh2o rather than denice.

@Ivanderkelen , do you agree? I'm not sure whether this will require some other changes in this branch as well for consistency, or if you can just make that one change. If you agree, can you make that change?

@Ivanderkelen
Copy link
Contributor

@billsacks I agree with this change - makes most sense. I don't think this requires additional changes (my basic test case compiles). The change is made in a672692 and I created a new PR (billsacks#2) to merge into your dynlakes_master_notools branch.

@billsacks
Copy link
Member Author

billsacks commented Aug 27, 2020

I have confirmed that the differences in SMS_D_Ld1.f09_g17.I1850Clm50Sp.cheyenne_intel.clm-default disappear (except for HEAT_CONTENT1, with these diffs:

diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90
index b774fe976..f52228e5d 100644
--- a/src/biogeophys/TotalWaterAndHeatMod.F90
+++ b/src/biogeophys/TotalWaterAndHeatMod.F90
@@ -449,11 +449,13 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, &
           ice_mass(c) = ice_mass(c) + h2osoi_ice(c,j)
        end do
     end do
-    
-    ! Lake water content
-    call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, &
-       lakestate_inst, liquid_mass, ice_mass)
-    
+
+    if (.false.) then
+       ! Lake water content
+       call AccumulateLiqIceMassLake(bounds, num_lakec, filter_lakec, &
+            lakestate_inst, liquid_mass, ice_mass)
+    end if
+
     ! Soil water content of the soil under the lake
     do j = 1, nlevgrnd
        do fc = 1, num_lakec

If instead of the above diffs I use these diffs:

diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90
index b774fe976..87d3256bb 100644
--- a/src/biogeophys/TotalWaterAndHeatMod.F90
+++ b/src/biogeophys/TotalWaterAndHeatMod.F90
@@ -518,7 +518,7 @@ contains
           c = filter_c(fc)
           ! calculate lake liq and ice content per lake layer first
           h2olak_liq = dz_lake(c,j) * denh2o * (1 - lake_icefrac(c,j))
-          h2olak_ice = dz_lake(c,j) * denice * lake_icefrac(c,j)
+          h2olak_ice = dz_lake(c,j) * denh2o * lake_icefrac(c,j) ! use water density of liquid water as layer depth is not adjusted
           
           liquid_mass(c) = liquid_mass(c) + h2olak_liq
           ice_mass(c) = ice_mass(c) + h2olak_ice

I am back to getting some differences, but the runoff diffs are much smaller than before:

 RMS ERRH2O                           1.0973E-20            NORMALIZED  3.5358E-06
 RMS HEAT_CONTENT1                    4.7075E+09            NORMALIZED  2.7918E+00
 RMS ICE_CONTENT1                     7.6008E+01            NORMALIZED  2.2591E-02
 RMS LIQUID_CONTENT1                  1.2883E+04            NORMALIZED  7.5430E+00
 RMS QRGWL                            5.7617E-18            NORMALIZED  1.5036E-11
 RMS QRUNOFF                          5.7617E-18            NORMALIZED  1.2437E-12
 RMS QRUNOFF_TO_COUPLER               5.7617E-18            NORMALIZED  1.2437E-12
 RMS TWS                              1.2888E+04            NORMALIZED  6.8277E-01
 RMS VOLR                             3.5449E-06            NORMALIZED  2.4071E-12
 RMS VOLRMCH                          3.5449E-06            NORMALIZED  2.9088E-11

and

 RMS l2x_Flrl_rofgwl                  3.8857E-16            NORMALIZED  8.2604E-10
 RMS x2l_Flrr_volr                    8.0947E-16            NORMALIZED  3.2027E-12
 RMS x2l_Flrr_volrmch                 8.0947E-16            NORMALIZED  2.5766E-11
 RMS r2x_Forr_rofl                    2.4679E-17            NORMALIZED  2.1198E-10
 RMS r2x_Flrr_volr                    5.7932E-16            NORMALIZED  5.9428E-12
 RMS r2x_Flrr_volrmch                 5.7932E-16            NORMALIZED  4.7794E-11
 RMS x2r_Flrl_rofgwl                  2.5315E-17            NORMALIZED  2.1745E-10

as opposed to:

 RMS ERRH2O                           2.2571E-20            NORMALIZED  7.2729E-06
 RMS HEAT_CONTENT1                    4.7075E+09            NORMALIZED  2.7918E+00
 RMS ICE_CONTENT1                     6.9700E+01            NORMALIZED  2.0718E-02
 RMS LIQUID_CONTENT1                  1.2883E+04            NORMALIZED  7.5430E+00
 RMS QRGWL                            6.0332E-07            NORMALIZED  1.4765E+00
 RMS QRUNOFF                          6.0332E-07            NORMALIZED  1.2949E-01
 RMS QRUNOFF_TO_COUPLER               6.0332E-07            NORMALIZED  1.2949E-01
 RMS TWS                              1.2887E+04            NORMALIZED  6.8276E-01
 RMS VOLR                             1.1477E+05            NORMALIZED  7.7577E-02
 RMS VOLRMCH                          1.1477E+05            NORMALIZED  8.9267E-01

and

 RMS l2x_Flrl_rofgwl                  6.2745E-07            NORMALIZED  1.2633E+00
 RMS x2l_Flrr_volr                    2.4199E-05            NORMALIZED  9.5066E-02
 RMS x2l_Flrr_volrmch                 2.4199E-05            NORMALIZED  7.2855E-01
 RMS r2x_Forr_rofl                    1.4068E-07            NORMALIZED  1.2134E+00
 RMS r2x_Flrr_volr                    1.7075E-05            NORMALIZED  1.7393E-01
 RMS r2x_Flrr_volrmch                 1.7075E-05            NORMALIZED  1.3330E+00
 RMS x2r_Flrl_rofgwl                  2.2763E-07            NORMALIZED  1.8302E+00

I think these new runoff diffs are expected: we are adding the total lake water to begwb and endwb, but not subtracting the baselines, so we end up with much bigger magnitude of the states, and I think this swamps the flux terms, leading to less precision in the flux terms and thus these greater-than-roundoff-level diffs, in a relative sense. I'm hopeful that, if we combined this with subtracting the baselines, we'd be back to only roundoff-level diffs (but I don't think I'm going to pursue that for now).

These findings support my earlier hypotheses; namely:

  • The big runoff diffs seem to be arising from the use of denice where we should be using denh2o for consistency
  • Even after correcting that, we have diffs, but they can be removed by reverting to the previous accounting of lake column total water. My next step will be make this accounting depend on a flag that differs in the different places from which these routines are called.

Ivanderkelen and others added 3 commits August 28, 2020 00:50
Adjust ice to liquid density in ice mass calculation

Adjust ice to liquid density in ice mass calculation, as lake layer is not adjusted for change in density.
We have changed the definition of total column water for lake columns,
so the baseline values for lakes are incorrect on old initial conditions
files. This commit adds some code to check if we're using an old initial
conditions file, and if so, resets the dynbal baseline values for lakes
to use the new definition.
I went back and forth about whether we should do this, but I actually
feel that it's best if we do reset the lake baselines in a branch or
continue run, if using an older restart file. If we didn't do this, we'd
want to add some logic for writing out the issue-fixed metadata for any
further restart files written from these runs, to note that this issue
isn't actually fixed yet on these restart files.
These seem to have been missing for a while (forever?).
@billsacks
Copy link
Member Author

billsacks commented Sep 9, 2020

@Ivanderkelen please see my proposed changes in billsacks#3, including my detailed comments in that PR explaining my changes and describing what I'd like from you. Please let me know if you'd like any help with how to do this, or if you'd like to talk more about any of this.

  • Once those changes are approved and tested, I'll merge that branch into this one

I was getting a gridcell-level methane balance check error when running
ERS_Lm25.1x1_smallvilleIA.IHistClm50BgcCropQianRsGs.cheyenne_gnu.clm-smallville_dynlakes_monthly.
I'm thinking it was due to ESCOMP#43. So for now, I'm disabling
methane in this testmod. (I thought about changing it to an Sp case, but
(1) currently you're not allowed to run smallville in an Sp case, and
(2) it would be nice to exercise as much of the BGC code as possible in
this test, to exercise the dynamic landunit conservation code for BGC
variables.)
Dynlakes: fix some subtle issues

### Description of changes

Fix a variety of subtle issues with dynamic lakes - particularly the accounting of total water and energy.

### Specific notes

This branch contains the following commits:

- a9fa875: This is needed to avoid counting lake water in the begwb and endwb terms, which is needed because these are used to calculate gridcell total water store (TWS), which in turn influences the methane code. Because the methane code was tuned around old values of TWS, changing TWS would lead to unintentional – and potentially large – changes in methane terms. Eventually we'd like to remove methane's dependence on TWS, but for now this workaround is needed to avoid changing behavior too much. See ESCOMP#659 (comment) for more details.

- 52105c4: This minor fix is needed for the sake of water tracers / water isotopes. It shouldn't have any impact outside of that (because the tracer_ratio of bulk water is 1)

- de3e12c and acf0984: This one is especially subtle; it is needed for backwards compatibility with old restart files. The main changes are in de3e12c; acf0984 is just a minor tweak on top of that. The problem is that, on existing initial conditions files, there can be already-existing DYNBAL_BASELINE variables (for LIQUID, ICE & HEAT). But these pre-existing variables will have baseline values of 0 for lake. Before this commit, when you started up from an old initial conditions file, the code would use these 0 values for lake baselines (because baselines are only reset if the user explicitly asks them to be reset with a namelist flag). This commit adds some code to detect if the initial conditions file is old, and if so, recomputes dynbal baselines for lake using the new definition. Note that some even older initial conditions files didn't have the DYNBAL_BASELINE variables at all; those would have been okay before this change: the problem is with initial conditions files that are somewhat but not very old - so have DYNBAL_BASELINE variables on them that use the old definition (where lake baseline values were 0).

- 8088c3c: Minor fix for a pre-existing issue

- a31875d: I'm not sure if this is actually needed, but I thought it would be good to group together the lake water content and the roughly equal-but-opposite baselines, so that these can cancel to near zero before adding the smaller terms. In principle, this should help maintain precision in these smaller terms. I thought this might help resolve some of the larger-than-expected answer changes I was seeing in testing, but I don't think it actually does... but I still thought this would be good to keep in place. **I have double and triple checked these changes, but it would be good to have an extra set of eyes on them to make sure I did this reordering correctly. In particular, I think there were some subtleties about when a term should accumulate on top of an existing value vs. setting the initial value of a variable.**
The important part of this is that it has PCT_LAKE + PCT_CROP <= 100%
for all times.
@billsacks billsacks merged commit d9b4972 into ESCOMP:master Sep 29, 2020
@billsacks billsacks deleted the dynlakes_master_notools branch September 29, 2020 16:48
@samsrabin samsrabin added enhancement new capability or improved behavior of existing capability science Enhancement to or bug impacting science labels Aug 8, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement new capability or improved behavior of existing capability PR status: ready PR: this is ready to merge in, with all tests satisfactory and reviews complete science Enhancement to or bug impacting science
Projects
Status: Done (non release/external)
Development

Successfully merging this pull request may close these issues.

Add lake water to dynbal baselines Count energy of water in lakes in total gridcell heat content
3 participants