-
Notifications
You must be signed in to change notification settings - Fork 321
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
Gridcell-level balance checks for carbon and nitrogen #984
Conversation
Initial mods to get things set up
src/biogeochem/CNBalanceCheckMod.F90
Outdated
grc_begcb => this%begcb_grc , & ! Output: [real(r8) (:)] (gC/m2) gridcell carbon mass, beginning of time step | ||
grc_begnb => this%begnb_grc , & ! Output: [real(r8) (:)] (gN/m2) gridcell nitrogen mass, beginning of time step | ||
! totgrcc => cnveg_carbonstate_inst%totc_grc , & ! Input: [real(r8) (:)] (gC/m2) total gridcell carbon, incl veg and cpool | ||
! totgrcn => cnveg_nitrogenstate_inst%totn_grc, & ! Input: [real(r8) (:)] (gN/m2) total gridcell nitrogen, incl veg |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In this first attempt to set up the framework for gridcell-level balance checks for carbon, I have commented out my own references to totgrcc
and totgrcn
throughout, because I may not need them. Instead, I think I just need:
- all gridcell-level pools and fluxes not accounted for at the column level and
- the column-level balance error converted to the gridcell level using
call c2g
src/biogeochem/CNBalanceCheckMod.F90
Outdated
! calculate the total gridcell-level carbon balance error | ||
! for this time step | ||
grc_errcb(g) = (grc_cinputs - grc_coutputs) * dt - & | ||
(grc_endcb(g) - grc_begcb(g)) + grc_errcb(g) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In reference to my previous comment,
this where I picture using all the gridcell-level terms not accounted for at the column level PLUS the column-level error converted to the gridcell level. Pls let me know if you think that this is the wrong way to think about it @billsacks
Either way, I have encountered a small problem at this stage:
the gridcell-level error calculated just from
call c2g( bounds = bounds, carr = col_errcb(bounds%begc:bounds%endc), garr = grc_errcb(bounds%begg:bounds%endg), c2l_scale_type = 'unity', l2g_scale_type = 'unity')
(all other terms set to 0) exceeds the 1e-7 threshold (largest value 8.6E-5) and aborts. Does this just reflect the level of accuracy of the c2g calculation? I assume so.
@slevisconsulting thanks a lot for starting to work on this. I'm not sure whether it would work to do your suggestion of just using the column-level error rather than the actual column-level terms. However, I'd say that the number one point of adding these gridcell-level checks is to give us an easy way to be confident that the system is conserving mass. To me, using the column-level error is too subtle and clever to give me this easy confidence, even if it turns out to work. Actually, as I think about it more, I don't think your suggestion will work in the dynamic landcover case, because the starting gridcell amount is counted at a different point from the starting column amount. (Also, maybe a related issue: I think your suggestion becomes problematic if there are fluxes between columns that we want to ignore in the gridcell-level checks, since the gridcell-level checks should only consider fluxes into or out of the gridcell as a whole.) Maybe others see things differently, but the way that I feel confident about conservation is to treat the system as a black box: Count up all of the fluxes in and out of the box, and the total change in mass of the box, and compare the two. So I would personally start like that. Does that make sense? (I want to draw a picture, but it's hard to do in a github comment.) Now, as I mentioned in #314 , if you do this in the simplest possible way, you will end up with a lot of duplicated logic between the column-level and gridcell-level checks. So an important don't-repeat-yourself cleanup step would be to only have one place in the code where we add up all of the inputs & outputs that are used in both the column-level and gridcell-level checks. However, you could start by duplicating some code and then do this cleanup later if that's an easier way to approach the problem. As for the error being exceeded after the c2g call: No, I would not expect this, and it could indicate a problem. I'd suggest printing the individual column-level values and subgrid weights for a given gridcell where you're seeing this, and then we can try to understand what's going on. The first thing that comes to mind is wondering whether special landunits might be contaminating the average. |
Thank you @billsacks , your response helps me. And regarding the error being exceeded after the c2g call, I suspected the same, i.e. special landunits possibly contaminating the average. I'm leaning against investigating this unless you feel strongly in favor. |
I'm fine with you not investigating, given that you're going to rework your approach. |
Increase tolerance on near-zero truncation for a snow state update In UpdateState_TopLayerFluxes, the tolerance of 1.e-13 was occasionally exceeded. Although I haven't done a careful analysis, it seems okay to me to increase this tolerance slightly. Here I increase it to 1.e-12. Previously, if top-layer h2osoi_ice or h2osoi_liq were reduced to less than 1e-13 times the original value (in an absolute value sense), these masses were set to 0; now we set these masses to 0 if they are reduced to less than 1e-12 times the original. So we can now occasionally set a value to exactly 0 when before it was left at slightly different from zero. If the previous code led to a small positive value, between 1e-13 and 1e-12 times the original, this tag will change answers slightly. If the previous code led to a small negative number, it would cause the model to abort, leading to the issue reported in ESCOMP#988; this change should fix those occasional aborts. This tag also introduces the general ability to set tolerances to a custom value in calls to truncate_small_values. See ESCOMP#988 for details. Resolves ESCOMP#988
Based on 4/21/2020 conversation with @billsacks.
Now just getting C balance error.
@billsacks I requested a quick review of my work to ensure that we are on the same page. |
This test PASSes: ERI_D_Ld9.f10_f10_musgs.I1850Clm50Bgc.cheyenne_gnu.clm-default Seems too good to be true, so now I am interested in a complete review of my changes.
@slevisconsulting thanks for pointing me to your next iteration of this. I have just looked quickly at this point, but from what I see, this looks basically like what I imagined now. However, I have some questions about some of the terms you include in
this%landuseflux_grc(g) = &
this%dwt_conv_cflux_dribbled_grc(g) + &
product_closs_grc(g)
I could easily be missing something with some / all of these, but it would be great to have an explanation of why these are needed. |
@billsacks thank you for catching these. I do not have good justification for including them other than my incomplete grasp of what each term is for. I will remove them in the next commit. |
I should be clear that I haven't done a careful analysis myself. The addition of these balance checks warrants some thinking about which terms should and shouldn't be included; I haven't done that thinking, so it would be great if you could. |
Will do. |
Followed my template that worked for the gridcell-level C balance check. This test PASSes ERI_D_Ld9.f10_f10_musgs.I1850Clm50Bgc.cheyenne_gnu.clm-default Next must correct C balance check for transient simulations.
...though still getting a C balance error. Removed a line that was repeating in CNVegCarbonFluxType.F90
This test now PASSes ERP_D_Ld5.f10_f10_musgs.IHistClm50BgcCrop.cheyenne_intel.clm-allActive
Test continues to PASS as in last commit.
This test now PASSes: ERS_Lm40_Mmpi-serial.1x1_numaIA.I2000Clm50BgcCropQianGs.cheyenne_gnu.clm-monthly
@billsacks |
Thinking through this again, it's doubtful that I'm right, and I'm continuing to investigate. Sorry to ping you about it @billsacks. |
@slevisconsulting I'll reply briefly even though you said I didn't need to: When a patch disappears, its states should still exist in memory, but it is removed from the list of "active" points, and therefore is removed from the standard filters. If you're summing up weighted totals over the grid cell, this shouldn't matter (since a patch that has disappeared will have 0 area, and thus contribute 0 to the grid cell total), but I wanted to mention this in case it helps explain anything you're seeing. And note that, if you really need to access patch-level states including patches that have disappeared you can use the |
This commit includes an answer-changing modification in CNVegCarbonFluxType's subroutine Summary_carbonflux: For carbon and nitrogen to balance at the gridcell level, we must account for the terms hrv_xsmrpool_to_atm_grc and dwt_conv_[c,n]flux_grc in a single timestep rather then by dribbling over the course of a year, unless allows_non_annual_delta = .true. The modification affects the corresponding dribbled terms, which in turn affect nee_grc, landuseflux_grc, nbp_grc, and fco2, all of which are diagnostic in the land model. As a result, all other variables in history output appear unchanged from their baseline values.
@billsacks this is a good stopping point for a code review. Next in line is the gridcell-level balance of methane #315 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@slevisconsulting - Thanks a lot for your work on this! At a high level this looks good, and I think it's just about ready to go if you feel so, but I have a number of requests - please see my specific comments below.
It might be good to make a tag just with these C&N balance changes, though if you want to also put the gridcell-level methane balance checks into the same tag, that also seems okay, if those changes are relatively easy. But I'd at least suggest doing the gridcell-level water & energy balance checks in a separate tag, since they'll likely be pretty involved, and this PR is already getting big.
src/biogeochem/CNProductsMod.F90
Outdated
@@ -33,10 +33,10 @@ module CNProductsMod | |||
class(species_base_type), allocatable :: species ! C, N, C13, C14, etc. | |||
|
|||
! States | |||
real(r8), pointer :: cropprod1_grc(:) ! (g[C or N]/m2) grain product pool, 1-year lifespan | |||
real(r8), pointer, public :: cropprod1_grc(:) ! (g[C or N]/m2) grain product pool, 1-year lifespan |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- Please move this a few lines above, to the section labeled "Public instance variables"
src/biogeochem/CNProductsMod.F90
Outdated
real(r8), pointer :: prod10_grc(:) ! (g[C or N]/m2) wood product pool, 10-year lifespan | ||
real(r8), pointer :: prod100_grc(:) ! (g[C or N]/m2) wood product pool, 100-year lifespan | ||
real(r8), pointer :: tot_woodprod_grc(:) ! (g[C or N]/m2) total wood product pool | ||
real(r8), pointer, public :: tot_woodprod_grc(:) ! (g[C or N]/m2) total wood product pool |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- Please move this a few lines above, to the section labeled "Public instance variables"
src/biogeochem/CNProductsMod.F90
Outdated
call restartvar(ncid=ncid, flag=flag, & | ||
varname=this%species%rest_fname('tot_woodprod', suffix='_g'), & | ||
xtype=ncd_double, dim1name='gridcell', & | ||
long_name='', units='', & | ||
interpinic_flag='interp', readvar=readvar, data=this%tot_woodprod_grc) | ||
! Backwards compatibility mentioned below is not applicable for this var. | ||
! If field not found in restart, then set from template if provided | ||
if (flag == 'read' .and. .not. readvar .and. template_provided) then | ||
call set_missing_from_template(this%tot_woodprod_grc, & | ||
template_for_missing_fields%tot_woodprod_grc, & | ||
multiplier = template_multiplier) | ||
end if | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- I don't think
tot_woodprod
really needs to be a restart variable. We're trying to avoid adding variables to the restart file unless they're really needed, to keep the file size down and improve initialization performance. In this case, I think what you could do is to put a block of code at the end of the restart routine, like this:
if (flag == 'read') then
call this%ComputeSummaryVars(bounds)
end if
@@ -3956,7 +3955,6 @@ subroutine SetValues ( this, & | |||
this%litfire_col(i) = value_column | |||
this%somfire_col(i) = value_column | |||
this%totfire_col(i) = value_column | |||
this%fire_closs_col(i) = value_column |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- Can you please explain why you needed to remove this setting of
fire_closs_col
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
An exact copy of the removed line appears 7 lines down.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Got it - thanks for the clean up!
hrv_xsmrpool_to_atm_dribbled_grc = hrv_xsmrpool_to_atm_dribbled_grc | ||
this%dwt_conv_cflux_dribbled_grc = this%dwt_conv_cflux_dribbled_grc |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- It looks like these lines of code aren't doing anything: they're just setting variables to themselves. Please remove them. And actually, together with the comment below, I'm suggesting removing this whole larger block that you added.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Removed these two lines.
! In this case overwrite the dribbled array with the non-dribbled array | ||
! so as to conserve carbon and nitrogen at every timestep | ||
hrv_xsmrpool_to_atm_dribbled_grc = hrv_xsmrpool_to_atm_grc | ||
this%dwt_conv_cflux_dribbled_grc = this%dwt_conv_cflux_grc |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- This doesn't seem right to me. It seems like this is going to change answers in an undesirable way. I think I see why you did this, but the solution differs from this. The key point is that we have these "flux dribblers" that effectively have some new state that needs to be accounted for in the gridcell-level balance checks. You can think of them like product pools: there is some flux into the dribblers at the start of each year, then there is some gradual flux out of them. Even though they don't have an explicit state variable associated with them, there is a pseudo-state representing the amount left to dribble over the year, and this pseudo-state needs to be included in the balance checks.
See this comment at the top of src/utils/AnnualFluxDribbler.F90
:
! And, for the sake of checking conservation:
!
! - To get gridcell water (or whatever) content at the start of the time step:
!
! call mydribbler%get_amount_left_to_dribble_beg
!
! - To get gridcell water (or whatever) content at the end of the time step:
!
! call mydribbler%get_amount_left_to_dribble_end
!
! These both return the pseudo-state representing how much of the original delta
! still needs to be dribbled. The 'beg' version includes the amount left to dribble
! in the current time step; the 'end' version does not.
These get_amount_left_to_dribble
subroutines aren't currently used in the code, but I put them in because I knew they would be needed for the sake of balance checks. Now is the time when they are becoming needed. So you'll need to introduce calls to these get_amount_left_to_dribble
routines for each dribbler in the code, adding that amount to the beginning and ending carbon balance.
As a related side-note, my hope is that this PR can be done in a way that doesn't change answers for any tests in the test suite.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Working on this next.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@billsacks fyi: I am pulled away from this PR again as I investigate WRF-CTSM simulations.
src/biogeochem/CNBalanceCheckMod.F90
Outdated
grc_begcb => this%begcb_grc , & ! Output: [real(r8) (:)] (gC/m2) gridcell carbon mass, beginning of time step | ||
grc_begnb => this%begnb_grc , & ! Output: [real(r8) (:)] (gN/m2) gridcell nitrogen mass, beginning of time step | ||
totgrcc => cnveg_carbonstate_inst%totc_grc , & ! Input: [real(r8) (:)] (gC/m2) total gridcell carbon, incl veg and cpool | ||
totgrcn => cnveg_nitrogenstate_inst%totn_grc, & ! Input: [real(r8) (:)] (gN/m2) total gridcell nitrogen, incl veg | ||
cropprod1c_grc => c_products_inst%cropprod1_grc , & ! Input: [real(r8) (:)] (gC/m2) carbon in crop products | ||
cropprod1n_grc => n_products_inst%cropprod1_grc , & ! Input: [real(r8) (:)] (gC/m2) nitrogen in crop products | ||
tot_woodprodc_grc => c_products_inst%tot_woodprod_grc, & ! Input: [real(r8) (:)] (gC/m2) total carbon in wood products | ||
tot_woodprodn_grc => n_products_inst%tot_woodprod_grc & ! Input: [real(r8) (:)] (gC/m2) total nitrogen in wood products |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- [OPTIONAL] Can you rename the variables on the left-hand side of the associates so that the variable naming is more consistent here?
First, it seems like you don't need the grc
on the variable names here, since this whole routine operates on gridcell-level quantities. But if you do keep the grc
, then please be consistent about its placement in variable names; I think it would be best if you kept the _grc
suffix at the end of variable names.
Second, related to the products variable... and this point is mainly important if you keep the _grc
: For the sake of searching the code, I find it's best to avoid breaking up a variable name, if possible. For example, if someone wanted to search for cropprod1_grc
, this wouldn't show up if the variable is named cropprod1c_grc
, but if you rename it like c_cropprod1_grc
it would show up – so I would recommend something like the latter.
Update: I see that you probably just followed what was done in BeginCNColumnBalance
and in other existing routines. Probably names should be changed there, too, to keep things consistent. However, I don't feel strongly that this absolutely needs to be done, here or elsewhere.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Going with both suggestions. However, I left the c
and n
unchanged in totc, totn, begcb, begnb
because these variables are more widespread. For the sake of easier searching and finding of these variables, I think leaving them unchanged is better for now.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Good point, thanks.
src/biogeochem/CNBalanceCheckMod.F90
Outdated
! slevis: Not including seedc_grc in grc_begcb and grc_endcb because | ||
! seedc_grc equals | ||
! -1 * (dwt_seedc_to_leaf_grc(g) + dwt_seedc_to_deadstem_grc(g)) | ||
! and we account for the latter fluxes as inputs below; the same | ||
! fluxes have entered the pools earlier in the timestep. For true | ||
! conservation we would need to add a flux out of npp into seed. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thank you for adding these notes. I'm trying to get my mind back around this point, so please correct me if I'm wrong here: My understanding is that seedc is basically an accounting term that balances the carbon we have pulled out of thin air to seed new patches. So we could either include seedc in the gridcell balance OR include dwt_seedc_to_leaf_grc(g) + dwt_seedc_to_deadstem_grc(g)
, but we should NOT include both.
- [OPTIONAL] If you find my explanation helpful, you could include it in this comment. If not, though, you can leave the comment as is.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Added some of your wording to the comment.
src/biogeochem/CNBalanceCheckMod.F90
Outdated
! calculate total gridcell-level inputs | ||
grc_ninputs(g) = forc_ndep(g) + nfix_to_sminn_grc(g) + & | ||
supplement_to_sminn_grc(g) + & | ||
dwt_seedn_to_leaf_grc(g) + & | ||
dwt_seedn_to_deadstem_grc(g) | ||
|
||
if (use_fun) then | ||
grc_ninputs(g) = grc_ninputs(g) + ffix_to_sminn_grc(g) | ||
endif | ||
|
||
if (use_crop) then | ||
grc_ninputs(g) = grc_ninputs(g) + fert_to_sminn_grc(g) + soyfixn_to_sminn_grc(g) | ||
end if | ||
|
||
! calculate total gridcell-level outputs | ||
grc_noutputs(g) = denit_grc(g) + grc_fire_nloss(g) + & | ||
dwt_conv_nflux_grc(g) + product_loss_grc(g) - & | ||
som_n_leached_grc(g) | ||
|
||
if (.not. use_nitrif_denitrif) then | ||
grc_noutputs(g) = grc_noutputs(g) + sminn_leached_grc(g) | ||
else | ||
grc_noutputs(g) = grc_noutputs(g) + f_n2o_nit_grc(g) + & | ||
smin_no3_leached_grc(g) + smin_no3_runoff_grc(g) | ||
end if |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- There is a lot of duplication in logic between the setting of these gridcell-level variables and the settings of the equivalent column-level variables above. Can you please remove this duplication by just averaging
col_ninputs
andcol_noutputs
to gridcell level? From a quick glance, it looks like all of the terms incol_ninputs
andcol_noutputs
are also included in the gridcell level terms, but if there are differences, you could save an intermediate version of the column-level terms after adding up all terms that apply to both the col and gridcell level, then add any terms that just apply to the column level. Similarly, you could add any terms that just apply to gridcell level later.
By the way: I thought about suggesting something similar for carbon, but I actually like that, for carbon, you used the existing nbp term.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Done.
I agree about making this a separate tag, before even going into the methane balance. |
src/biogeochem/CNBalanceCheckMod.F90
Outdated
real(r8):: sminn_leached_grc(bounds%begg:bounds%endg) | ||
real(r8):: f_n2o_nit_grc(bounds%begg:bounds%endg) | ||
real(r8):: smin_no3_leached_grc(bounds%begg:bounds%endg) | ||
real(r8):: smin_no3_runoff_grc(bounds%begg:bounds%endg) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
...just noticed I need to delete a bunch of these here and elsewhere!
These tests PASS: @billsacks I wanted you to check my work around the use of |
Minor bug fixes needed for latest cime (1) Refactor a conditional to avoid divide by 0 (ESCOMP#1013) (2) Remove some FATES coordinate variables that were writing garbage (temporary workaround for ESCOMP#1014) (3) Minor update to manage_externals (4) Document authorship in README.md and .zenodo.json
@slevisconsulting this looks very good to me now. Thank you! I'll send an email to discuss tag planning. |
Update cime and cmeps externals; rework initialization of CNFire object Main change is to update cime and cmeps externals. cime is now essentially at the version used in cesm2_2_beta05 (but with one change backed out); cmeps is at latest master. Also, reworks initialization of CNFire object to avoid occasional segmentation faults.
@billsacks pls let me know if my ChangeLog entry seems sufficient to you. Other than that the tests have now passed and this is ready. |
Looks good to me; just to confirm: should I go ahead and merge this now? |
Description of changes
Bracket the entire loop to calculate balance checks at the gridcell level, as explained in detail in #314
Specific notes
Avoid duplication by using calculations from the column-level checks.
Contributors other than yourself, if any:
@billsacks
CTSM Issues Fixed (include github issue #):
#314
Are answers expected to change (and if so in what way)?
No
Any User Interface Changes (namelist or namelist defaults changes)?
No
Testing performed, if any:
Initially using a single test to assess where things stand:
./create_test ERI_D_Ld9.f10_f10_musgs.I1850Clm50Bgc.cheyenne_gnu.clm-default