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

Annual sowing and harvest date outputs #1616

Merged
merged 29 commits into from
Mar 28, 2022
Merged

Annual sowing and harvest date outputs #1616

Show file tree
Hide file tree
Changes from 9 commits
Commits
Show all changes
29 commits
Select commit Hold shift + click to select a range
3f2acc2
Merge remote-tracking branch 'upstream/master'
samsrabin Oct 11, 2021
7d0acfe
Git now ignores .vscode/.
samsrabin Oct 28, 2021
65b0295
Merge branch 'master' of https://github.com/ESCOMP/CTSM
samsrabin Nov 1, 2021
ac85009
Merge tag 'ctsm5.1.dev073'
samsrabin Jan 26, 2022
9360778
Added annual outputs of sowing and harvest dates.
samsrabin Jan 26, 2022
3472154
Compile fix: Forgot to remove one line with an unused variable during…
samsrabin Jan 26, 2022
b857d8c
Use real() instead of DBLE().
samsrabin Jan 27, 2022
cc42589
Ensure that new variables are correctly set when resuming from a run …
samsrabin Jan 27, 2022
bc88cf9
REALLY ensure that new variables are correctly set when resuming from…
samsrabin Jan 28, 2022
6491189
Moved initialization of [sh]dates_thisyr to match convention.
samsrabin Feb 1, 2022
86a09e0
Improved a comment.
samsrabin Feb 2, 2022
ac01b14
CropType: Added [sh]dates_thisyr to restart files, with [sowing,harve…
samsrabin Feb 2, 2022
d1ca586
ncid now inout in CallRestartvarDimOK().
samsrabin Feb 4, 2022
7601ec8
More descriptive variable names.
samsrabin Feb 4, 2022
948bcc9
CallRestartvarDimOK now uses check_dim() instead of check_var_or_dim().
samsrabin Feb 4, 2022
b878b33
Git now ignores .vscode/.
samsrabin Oct 28, 2021
708d8e7
Merge branch 'master' into 1537-crop-date-outputs3
samsrabin Feb 4, 2022
afc1924
Merge remote-tracking branch 'billsacks/cropphenology_jday_startoftim…
samsrabin Feb 4, 2022
b89ea27
Minor .gitignore rearrange.
samsrabin Feb 4, 2022
7a2c816
Renamed mxgrowseas to mxsowings; added explanatory comment on mxharve…
samsrabin Feb 14, 2022
071b545
Merge remote-tracking branch 'billsacks/cropphenology_jday_startoftim…
samsrabin Feb 18, 2022
ca90b78
Elaboration of a comment related to backwards compatibility.
samsrabin Feb 18, 2022
8517646
Generalized backwards compatibility setting of sowing_count and sdate…
samsrabin Mar 7, 2022
3aa7d65
Now sets rootdepth for crops to 0 outside growing season.
samsrabin Mar 11, 2022
d810b0f
Merge branch 'master2' into 1537-crop-date-outputs3
samsrabin Mar 14, 2022
0c67c80
Add new crop history variables to some tests
billsacks Mar 24, 2022
a8eca51
Merge remote-tracking branch 'escomp/master' into 1537-crop-date-outp…
billsacks Mar 24, 2022
5de54d1
Fix unit test build
billsacks Mar 25, 2022
c50409a
In cropMonthOutput testmod, change the new hist file to be vector
billsacks Mar 26, 2022
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
3 changes: 3 additions & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -112,3 +112,6 @@ core.*
*.log !run.log
*.pyc
Depends

# ssr
.vscode/
78 changes: 28 additions & 50 deletions src/biogeochem/CNPhenologyMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ module CNPhenologyMod
use shr_log_mod , only : errMsg => shr_log_errMsg
use shr_sys_mod , only : shr_sys_flush
use decompMod , only : bounds_type
use clm_varpar , only : maxveg, nlevdecomp_full
use clm_varpar , only : maxveg, nlevdecomp_full, mxgrowseas, mxharvests
use clm_varpar , only : i_litr_min, i_litr_max
use clm_varctl , only : iulog, use_cndv
use clm_varcon , only : tfrz
Expand Down Expand Up @@ -1644,7 +1644,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , &
! handle CN fluxes during the phenological onset & offset periods.

! !USES:
use clm_time_manager , only : get_curr_date, get_curr_calday, get_curr_days_per_year, get_rad_step_size
use clm_time_manager , only : get_curr_calday, get_curr_days_per_year, get_rad_step_size, is_beg_curr_year
use pftconMod , only : ntmp_corn, nswheat, nwwheat, ntmp_soybean
use pftconMod , only : nirrig_tmp_corn, nirrig_swheat, nirrig_wwheat, nirrig_tmp_soybean
use pftconMod , only : ntrp_corn, nsugarcane, ntrp_soybean, ncotton, nrice
Expand Down Expand Up @@ -1673,15 +1673,12 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , &
type(cnveg_carbonstate_type) , intent(inout) :: c14_cnveg_carbonstate_inst
!
! LOCAL VARAIBLES:
integer kyr ! current year
integer kmo ! month of year (1, ..., 12)
integer kda ! day of month (1, ..., 31)
integer mcsec ! seconds of day (0, ..., seconds/day)
integer jday ! julian day of the year
integer fp,p ! patch indices
integer c ! column indices
integer g ! gridcell indices
integer h ! hemisphere indices
integer s ! growing season indices
integer idpp ! number of days past planting
real(r8) :: dtrad ! radiation time step delta t (seconds)
real(r8) dayspyr ! days per year
Expand Down Expand Up @@ -1718,8 +1715,9 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , &
leafout => crop_inst%gddtsoi_patch , & ! Input: [real(r8) (:) ] gdd from top soil layer temperature
harvdate => crop_inst%harvdate_patch , & ! Output: [integer (:) ] harvest date
croplive => crop_inst%croplive_patch , & ! Output: [logical (:) ] Flag, true if planted, not harvested
cropplant => crop_inst%cropplant_patch , & ! Output: [logical (:) ] Flag, true if crop may be planted
vf => crop_inst%vf_patch , & ! Output: [real(r8) (:) ] vernalization factor
sowing_count => crop_inst%sowing_count , & ! Inout: [integer (:) ] number of sowing events this year for this patch
harvest_count => crop_inst%harvest_count , & ! Inout: [integer (:) ] number of harvest events this year for this patch
peaklai => cnveg_state_inst%peaklai_patch , & ! Output: [integer (:) ] 1: max allowed lai; 0: not at max
tlai => canopystate_inst%tlai_patch , & ! Input: [real(r8) (:) ] one-sided leaf area index, no burying by snow

Expand Down Expand Up @@ -1751,7 +1749,6 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , &
! get time info
dayspyr = get_curr_days_per_year()
jday = get_curr_calday()
call get_curr_date(kyr, kmo, kda, mcsec)
dtrad = real( get_rad_step_size(), r8 )

if (use_fertilizer) then
Expand All @@ -1776,36 +1773,30 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , &
! from AgroIBIS subroutine planting
! ---------------------------------

! in order to allow a crop to be planted only once each year
! initialize cropplant = .false., but hold it = .true. through the end of the year

! initialize other variables that are calculated for crops
! on an annual basis in cropresidue subroutine

if ( jday == jdayyrstart(h) .and. mcsec == 0 )then

! make sure variables aren't changed at beginning of the year
! for a crop that is currently planted, such as
! WINTER TEMPERATE CEREAL = winter (wheat + barley + rye)
! represented here by the winter wheat pft

if (.not. croplive(p)) then
cropplant(p) = .false.
idop(p) = NOT_Planted

! keep next for continuous, annual winter temperate cereal crop;
! if we removed elseif,
! winter cereal grown continuously would amount to a cereal/fallow
! rotation because cereal would only be planted every other year

else if (croplive(p) .and. (ivt(p) == nwwheat .or. ivt(p) == nirrig_wwheat)) then
cropplant(p) = .false.
! else ! not possible to have croplive and ivt==cornORsoy? (slevis)
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.

What happens to the winter cereals by deleting this code? It seems like this was removed due to the deletion of cropplant above, but it's not clear to me the impact it might have.

Copy link
Contributor

Choose a reason for hiding this comment

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

Does the sowing_count logic replace the cropplant logic? If so, do we need to update the winter cereal code to use the sowing_count logic to allow for continuous annual winter temperate cereals (as the comment implies)?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

The sowing_count logic is indeed intended to replace the cropplant logic; setting sowing_count equal to 0 at the beginning of the calendar year should have the same effect as this bit of code for winter wheat.

However, you and that comment have got me wondering: It is possible for corn, soy, etc. to have growing seasons spanning the (hemisphere's) New Year. Since they would miss their only chance to get cropplant set to false, doesn't that mean they would not be planted again? E.g., for a northern-hemisphere gridcell:

  • Sown Oct. 1, 2022 (cropplant set to true)
  • Jan. 1, 2023: croplive true, so cropplant remains true
  • Harvested Jan. 3, 2023; croplive set to false
  • Sowing window 2023: Sowing block not entered because cropplant is still true
  • Jan. 1, 2024: croplive false, so cropplant set to false

Copy link
Member

Choose a reason for hiding this comment

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

Good questions. @samsrabin I think your latest comment implies that there might have been situations where the old code did the wrong thing but your new code would do the right thing; is that correct? I'm wondering if it would be important to do a long-ish run (e.g., 20 years or so?) comparing the new and old code to see if this new code changes answers in some situations. I'm not sure whether this is necessary or not; @samsrabin @danicalombardozzi what do you think? The purposes I would see of this are (1) understanding whether this code change fixes any bugs with the old code (intentionally or inadvertently) and if so, how often those bugs appeared; and (2) making sure that the new code isn't changing any answers in an undesirable way.

If we do this, we'll want to compare against #1628 . An easy starting point would be to run the test suite (or maybe first just a couple of tests from the test suite). But that only includes coarse-resolution crop tests (10x15 degrees) and only runs for a few years, so might not capture rare answer changes. Maybe that's enough – maybe we don't care if there are some rare situations where this branch changes answers – but I want to see what both of you think. Again, I'm happy going either way here.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I agree, a test like that seems like a good idea.

Copy link
Contributor

Choose a reason for hiding this comment

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

@samsrabin Thanks for pointing this out. I agree that it could have caused the original code to do the wrong thing. A test does sound like a good idea. Does running two 20-year 2-degree simulations, one with this new code and another with the original code sound reasonable? Checking the changes in crop yields would be interesting. I'm hoping that the differences aren't huge, but the test will tell!

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

@billsacks I think this test is the only thing left to do before this PR is ready to go. Is there a nice automated way to perform these runs and/or compare their results for important outputs?

! Second condition ensures everything is correctly set when resuming from a run with old code
if ( is_beg_curr_year() .or. crop_inst%sdates_thisyr(p,1) == spval ) then
Copy link
Member

@billsacks billsacks Feb 3, 2022 •
edited
Loading

Choose a reason for hiding this comment

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

I actually think there is another possible subtle issue here: CropPhenology is currently only called when doalb is .true. Based on a test I just ran, in an I compset, doalb is .false. for both time steps 0 and 1, so CropPhenology is first called in the second time step of Jan 1 – meaning this conditional wouldn't be triggered in the first year of the run. If we're starting from cold start or an old restart file, then the spval check here should handle that situation – doing this start-of-year resetting the first time CropPhenology is called, even though is_beg_curr_year will be false. However, I think this will currently cause problems if you start a new run from a restart file created from this branch (not via a real restart, i.e., CONTINUE_RUN, but a new startup case with finidat set to point to a restart file created from this branch). If I'm thinking about things right, in this situation, sdates_thisyr will typically not be spval (at least, not after you add it to the restart file, which I think is needed for other scenarios) but is_beg_curr_year won't be triggered in the first year either – and so sowing_count, harvest_count, etc. won't get reset in the first year in this scenario.

I think the easiest way to fix this is to get rid of the doalb check around the call to CropPhenology, but I want to get some confirmation from others that that is an okay thing to change. @samsrabin please comment if you have any thoughts on this.

Update (2022-02-16): The plan is to fix this via #1628 (which has now been merged in to this branch), so I'm marking this as done.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Yes, I think that logic works.

Copy link
Contributor

Choose a reason for hiding this comment

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

I agree that this logic should work. Looks as though some of it has already been included.

sowing_count(p) = 0
harvest_count(p) = 0
do s = 1, mxgrowseas
crop_inst%sdates_thisyr(p,s) = -1._r8
Copy link
Member

@billsacks billsacks Jan 27, 2022 •
edited
Loading

Choose a reason for hiding this comment

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

  • [OPTIONAL] For this and hdates_thisyr: would it make sense to use spval here instead of -1, so that years that don't use all of the array of possible dates show up as missing value rather than -1? I don't have a clear sense of what's best here, but just wanted to raise that idea.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

My thinking was that it could be useful to distinguish "cell-years that had no area of this crop" (which would save as missing) from "cell-years where this crop had area prescribed but sowing criteria were never met" (which would save as -1).

Copy link
Member

Choose a reason for hiding this comment

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

My thinking was that it could be useful to distinguish "cell-years that had no area of this crop" (which would save as missing) from "cell-years where this crop had area prescribed but sowing criteria were never met" (which would save as -1).

Ah, that reasoning makes sense to me - thanks for the explanation.

end do
do s = 1, mxharvests
crop_inst%hdates_thisyr(p,s) = -1._r8
end do
end if

! When resuming from a run with old code, may need to manually set these
if (jday == 1 .and. croplive(p) .and. idop(p) == 1 .and. sowing_count(p) == 0) then
sowing_count(p) = 1
crop_inst%sdates_thisyr(p,1) = 1._r8
end if
Copy link
Member

@billsacks billsacks Jan 31, 2022 •
edited
Loading

Choose a reason for hiding this comment

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

  • Thank you very much for thinking about what's needed for backwards compatibility and other edge cases! But I'm having trouble understanding the details of when this is needed, and would like to talk more about it. Is this really about resuming a run with old code, or might it be needed for a different reason? (I'm pushing on this because, for the sake of long-term maintenance, it helps us to know what was just needed for backwards compatibility and can eventually be dropped vs. what needs to be kept in place long-term.) Specifically, I'm wondering if this might be needed due to subtleties with how a non-restart run starts and how jday is determined: A non-restart run runs an initial time step 0 whose time values are the time just before the intended start of the run (see also Stop running 0th time step if possible, or at least clean up its handling in history output #925 ). If you start a non-restart run on the year boundary, this time step 0 will have time values corresponding to the last time step of Dec 31. jday in this subroutine is determined by a call to get_curr_calday, which looks at the time as of the end of the time step – which in this case I think will yield jday = 1. So I can imagine a situation where a crop whose planting window starts on Jan. 1 gets planted in this 0th time step. Then on the next time step (time step 1; i.e., the first time step of the year), your new reset code above this will be executed, resetting sowing_count, sdates_thisyr, etc. I can see how you might need a correction for this situation. And actually, it seems like this same situation could occur in later years as well. Does what I'm saying ring true? If so, I think you should either:
    1. Keep the current logic, but change the comment to give some of this explanation.
    2. Put in an explicit check that prevents planting from happening on the last time step of the year to avoid this weirdness. (This would probably be answer changing so I'd like to do that in its own branch that would come to master before this branch.)
    3. Do a more significant rework to avoid this issue. The first solution that comes to mind is changing jday in this subroutine to be based on get_prev_calday instead of get_curr_calday: that way, we'll be looking at the date as of the start of the time step rather than the end of the time step, which would be more consistent with the use of is_beg_curr_year to determine when to reset these variables. This is likely to be answer-changing, though, and so should be done in a separate PR; probably the best way to do that would be to make this change to jday first (and we could bring that to master), then update this branch with that change and remove this block of code.

I am inclined to go with option (3), and that was the motivation for #1623 . But I'm okay with the other options as well if you think they'd be better or just as good. It seems like option (1) has some elements of option (2) for years after the first model year, in that I think it will generally prevent planting from happening on the last time step of the year (i.e., the first time step that is currently labeled Jan 1), since sowing_count will still be 1 at that point. (And actually, as I'm writing that, I'm realizing that that aspect of the current PR could be answer changing to a small extent; we should be aware of that for testing if we stick with the current implementation.) Other than that, one possible issue I see with option (1) (i.e., the current implementation) is: if there is a year in which planting doesn't happen at all (so sowing_count has remained 0), but then planting happens on the last time step of the year (which is viewed as Jan 1), then the diagnostics would be wrong: both this year's and next year's diagnostics will say that planting happened on Jan 1.

So I guess as I'm writing this all out I'm starting to lean more towards option (3) if you're happy with it, but I'm happy to discuss further.

Thoughts?

Copy link
Member

@billsacks billsacks Feb 2, 2022 •
edited
Loading

Choose a reason for hiding this comment

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

In our discussion yesterday, @samsrabin helped me see the scenarios where this is needed. Our tentative plan is to do #1623 (my option (3) in the above comment), which will help with this and other issues that @samsrabin has run into. But even with that change, this block of code is still needed to handle our typical end-of-year restart files that were generated with code prior to #1623 being fixed: In the runs that generated these restart files, the last time step of the run was considered to be Jan 1, and so some crop patches have crops that have already been planted, with a planting date of Jan 1. As long as we are still using these restart files, we need this block of code to correctly handle these already-planted-on-Jan-1 crops.

The one remaining thing I'd like here, then, is some more details in this comment so that we know how long we need to maintain this block of code.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Okay, done (ca90b78).

Copy link
Member

Choose a reason for hiding this comment

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

Thanks!


if ( (.not. croplive(p)) .and. (.not. cropplant(p)) ) then

! Once outputs can handle >1 planting per year, remove 2nd condition.
if ( (.not. croplive(p)) .and. sowing_count(p) == 0 ) then
Copy link
Member

@billsacks billsacks Jan 27, 2022 •
edited
Loading

Choose a reason for hiding this comment

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

  • Is there potentially a behavior change here because now you're using the condition sowing_count(p) == 0 instead of (.not. cropplant(p))? At a glance it looks like the logic for sowing_count differs from the logic for the old cropplant, but I haven't looked closely. If there is a potential behavior change, can you explain what it is? One particular behavior of the new code appears to be that you only allow a single planting in a given calendar year. Did the old code have a similar restriction, including for the Southern Hemisphere?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

The intended behavior of cropplant was to restrict plantings to once per year. In CropPhenology() (note the comment at the top here):

CTSM/src/biogeochem/CNPhenologyMod.F90

Lines 1779 to 1806 in bf365e0

! in order to allow a crop to be planted only once each year
! initialize cropplant = .false., but hold it = .true. through the end of the year
! initialize other variables that are calculated for crops
! on an annual basis in cropresidue subroutine
if ( jday == jdayyrstart(h) .and. mcsec == 0 )then
! make sure variables aren't changed at beginning of the year
! for a crop that is currently planted, such as
! WINTER TEMPERATE CEREAL = winter (wheat + barley + rye)
! represented here by the winter wheat pft
if (.not. croplive(p)) then
cropplant(p) = .false.
idop(p) = NOT_Planted
! keep next for continuous, annual winter temperate cereal crop;
! if we removed elseif,
! winter cereal grown continuously would amount to a cereal/fallow
! rotation because cereal would only be planted every other year
else if (croplive(p) .and. (ivt(p) == nwwheat .or. ivt(p) == nirrig_wwheat)) then
cropplant(p) = .false.
! else ! not possible to have croplive and ivt==cornORsoy? (slevis)
end if
end if

This is the only place in the code (other than the initialization of the CropType instance) where cropplant or cropplant_patch is ever set to false.

Although the Southern Hemisphere's "year" spans across calendar years, all sowing date windows are restricted to a single calendar year (as discussed in #1484). I think it follows (although check me on this) that the Southern Hemisphere is also restricted to one planting per calendar year.

Copy link
Member

@billsacks billsacks Jan 30, 2022 •
edited
Loading

Choose a reason for hiding this comment

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

Okay, I think I understand. Does this sound right?: Previously, the code itself did not explicitly restrict planting dates to being once per calendar year (for the southern hemisphere, anyway). However, this restriction existed in practice due to a combination of (1) the current crop parameters, and (2) issue #1484 (which prevents the code from working right if the crop parameters were ever changed to try to have a planting window crossing the new year boundary).

  • [no] If that's correct, I'm happy with the change, but I would suggest changing the comment ("Once outputs can handle >1 planting per year, remove 2nd condition.") to give more details about the reasoning here and what we would like in practice. Something like the following:
! Currently we only allow one planting per calendar year.
!
! Ideally, we would not have this restriction, particularly in the Southern Hemisphere.
! However, we have this restriction right now for the following reasons:
! (1) For currently-parameterized crops, planting windows never cross the new year boundary
! (2) https://github.com/ESCOMP/CTSM/issues/1484 would need to be resolved to allow planting windows to cross the new year boundary
! (3) Outputs currently cannot handle more than one planting per calendar year

Update (2022-02-16) I am marking this as not needed, because upon further reflection, I don't think this comment change is necessary.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

You know, I'm wondering if it even makes sense to leave this extra condition in. If we changed

if ( (.not. croplive(p)) .and. sowing_count(p) == 0 ) then

to just

if ( (.not. croplive(p)) ) then

the restriction would still effectively be in place both with the sowing window method (due to them not spanning multiple calendar years) and the prescribed sowing date method I'm working on (due to hard-coded mxgrowseas = 1). If we do want to leave in such a restriction, it might make more sense to throw an actual warning or error message in PlantCrop() when it's called in a patch that already had a growing season started that year (i.e., sowing_count(p) > 0). This is wrong. Theoretically it could be possible for a crop to be planted at the beginning of its sowing window, then harvested before the end of the sowing window, thus opening up the possibility of it being planted again.

But I'm not sure I'm right that current outputs can't handle multiple growing seasons in a year. Any relevant daily outputs shouldn't break in-run, although I'm sure it would mess up a lot of people's postprocessing scripts. @slevisconsulting, 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.

Thank you for including me in this @samsrabin . I haven't looked at the crop code in long while, so I would need a bit of a refresher before expressing an opinion. I'm open to having a short meeting to discuss if you like.

Copy link
Member

Choose a reason for hiding this comment

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

I'm happy to let you, @slevisconsulting and @danicalombardozzi hash this out. However, reading through your explanation, and especially your point that, without this condition, there could theoretically end up being multiple plantings in a year, I would be inclined to leave it as you currently have it. It sounds like removing this restriction would involve (maybe among other things):

  1. Figuring out a way to prevent your theoretical possibility
  2. Resolving Bug in planting window logic (no effect given current parameters?) #1484
  3. Figuring out how to handle outputs sensibly if there can be two plantings in a given calendar year

If all of that is needed, then it feels to me like it warrants a separate PR (if it's worth doing at all).

But I don't feel strongly on this, so am happy to let you all work out what makes the most sense here.

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 not entirely sure I follow the old logic, but I do think it makes sense to constrain the code to one planting per year at this point. Multiple plantings is something that happens and the research community would like to see added, but it might need additional constraints to work realistically (e.g., second sowing only happens a certain amount of time after harvest).

I would think that current outputs should be able to handle this, but probably would need some testing to verify.


! gdd needed for * chosen crop and a likely hybrid (for that region) *
! to reach full physiological maturity
Expand All @@ -1826,22 +1817,6 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , &

if (ivt(p) == nwwheat .or. ivt(p) == nirrig_wwheat) then

! add check to only plant winter cereal after other crops (soybean, maize)
! have been harvested

! *** remember order of planting is crucial - in terms of which crops you want
! to be grown in what order ***

! in this case, corn or soybeans are assumed to be planted before
! cereal would be in any particular year that both patches are allowed
! to grow in the same grid cell (e.g., double-cropping)

! slevis: harvdate below needs cropplant(p) above to be cropplant(p,ivt(p))
! where ivt(p) has rotated to winter cereal because
! cropplant through the end of the year for a harvested crop.
! Also harvdate(p) should be harvdate(p,ivt(p)) and should be
! updated on Jan 1st instead of at harvest (slevis)

! Are all the normal requirements for planting met?
do_plant_normal = a5tmin(p) /= spval .and. &
a5tmin(p) <= minplanttemp(ivt(p)) .and. &
Expand Down Expand Up @@ -2073,6 +2048,8 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , &

else if (hui(p) >= gddmaturity(p) .or. idpp >= mxmat(ivt(p))) then
if (harvdate(p) >= NOT_Harvested) harvdate(p) = jday
harvest_count(p) = harvest_count(p) + 1
crop_inst%hdates_thisyr(p, harvest_count(p)) = real(jday, r8)
croplive(p) = .false. ! no re-entry in greater if-block
cphase(p) = 4._r8
if (tlai(p) > 0._r8) then ! plant had emerged before harvest
Expand Down Expand Up @@ -2241,8 +2218,8 @@ subroutine PlantCrop(p, leafcn_in, jday, &

associate( &
croplive => crop_inst%croplive_patch , & ! Output: [logical (:) ] Flag, true if planted, not harvested
cropplant => crop_inst%cropplant_patch , & ! Output: [logical (:) ] Flag, true if crop may be planted
harvdate => crop_inst%harvdate_patch , & ! Output: [integer (:) ] harvest date
sowing_count => crop_inst%sowing_count , & ! Inout: [integer (:) ] number of sowing events this year for this patch
idop => cnveg_state_inst%idop_patch , & ! Output: [integer (:) ] date of planting
leafc_xfer => cnveg_carbonstate_inst%leafc_xfer_patch , & ! Output: [real(r8) (:) ] (gC/m2) leaf C transfer
leafn_xfer => cnveg_nitrogenstate_inst%leafn_xfer_patch , & ! Output: [real(r8) (:) ] (gN/m2) leaf N transfer
Expand All @@ -2253,9 +2230,10 @@ subroutine PlantCrop(p, leafcn_in, jday, &
! impose limit on growing season length needed
! for crop maturity - for cold weather constraints
croplive(p) = .true.
cropplant(p) = .true.
idop(p) = jday
harvdate(p) = NOT_Harvested
sowing_count(p) = sowing_count(p) + 1
crop_inst%sdates_thisyr(p,sowing_count(p)) = jday

leafc_xfer(p) = initial_seed_at_planting
leafn_xfer(p) = leafc_xfer(p) / leafcn_in ! with onset
Expand Down
55 changes: 25 additions & 30 deletions src/biogeochem/CropType.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,11 +27,9 @@ module CropType
! Crop state variables structure
type, public :: crop_type

! Note that cropplant and harvdate could be 2D to facilitate rotation
integer , pointer :: nyrs_crop_active_patch (:) ! number of years this crop patch has been active (0 for non-crop patches)
logical , pointer :: croplive_patch (:) ! patch Flag, true if planted, not harvested
logical , pointer :: cropplant_patch (:) ! patch Flag, true if planted
integer , pointer :: harvdate_patch (:) ! patch harvest date
integer , pointer :: harvdate_patch (:) ! most recent patch harvest date; 999 if currently (or never) planted
real(r8), pointer :: fertnitro_patch (:) ! patch fertilizer nitrogen
real(r8), pointer :: gddplant_patch (:) ! patch accum gdd past planting date for crop (ddays)
real(r8), pointer :: gddtsoi_patch (:) ! patch growing degree-days from planting (top two soil layers) (ddays)
Expand All @@ -41,6 +39,10 @@ module CropType
character(len=20) :: baset_mapping
real(r8) :: baset_latvary_intercept
real(r8) :: baset_latvary_slope
real(r8), pointer :: sdates_thisyr (:,:) ! all actual sowing dates for this patch this year
real(r8), pointer :: hdates_thisyr (:,:) ! all actual harvest dates for this patch this year
integer , pointer :: sowing_count (:) ! number of sowing events this year for this patch
integer , pointer :: harvest_count (:) ! number of sowing events this year for this patch
Comment on lines +42 to +45
Copy link
Member

@billsacks billsacks Jan 31, 2022 •
edited
Loading

Choose a reason for hiding this comment

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

  • I think all four of these new variables need to be on the restart file in order to support a mid-year restart. Do you agree? If so, can you please add them to the Restart subroutine in this module? (Let me know if you need any guidance with this. In particular, adding the multi-level variables may be a bit tricky.) Actually, maybe you can get away with leaving sowing_count and harvest_count off of the restart file and recomputing those after reading sdates_thisyr and hdates_thisyr based on a loop over all points in those arrays? (If that isn't too hard, that would be preferable, because we struggle to keep the restart file from growing unnecessarily large.)

Copy link
Collaborator Author

@samsrabin samsrabin Feb 2, 2022 •
edited
Loading

Choose a reason for hiding this comment

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

I gave this a shot and ran into a problem. Here's what I added for reading sdates_thisyr (similar for hdates_thisyear):

call restartvar(ncid=ncid, flag=flag, varname='sdates_thisyr', xtype=ncd_double,  & 
     dim1name='pft', dim2name='mxgrowseas', switchdim=.true., &
     long_name='crop sowing dates for this patch this year', units='day of year', &
     scale_by_thickness=.false., &
     interpinic_flag='interp', readvar=readvar, data=this%sdates_thisyr)

But then I get the following error:

Abort with message NetCDF: Invalid dimension ID or name in file /glade/u/home/samrabin/ctsm/libraries/parallelio/src/c
lib/pio_nc.c at line 812

I'm assuming this is because mxgrowseas isn't in the restart file, right? Is there any way around this?

Cross-posted to the CESM Forum, although it's awaiting approval.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Wrapping it like so seems to do the trick:

call check_var_or_dim(ncid, 'mxgrowseas', is_dim=.true., exists=found)
if (found) then

Copy link
Member

Choose a reason for hiding this comment

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

I don't think you want that check_var_or_dim call. Instead, I think you need to add an ncd_defdim call in restFileMod.F90, similar to this:

CTSM/src/main/restFileMod.F90

Line 546 in 8d9f988

call ncd_defdim(ncid , 'numrad' , numrad , dimid)

Copy link
Member

Choose a reason for hiding this comment

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

Oh wait, I get it: maybe you already added that defdim call but the problem came in reading an earlier restart file? If so, then something like your proposal makes sense – to do a check before making the restartvar call.

Copy link
Member

@billsacks billsacks Feb 3, 2022 •
edited
Loading

Choose a reason for hiding this comment

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

Awesome, this looks great. Thanks for introducing a shared subroutine to do this check.

I have a few minor comments on your latest commit, which I'm marking as optional because if you don't have time to do them I won't hold up the PR for them:

  • [OPTIONAL] If you declare ncid as intent(inout) in your new subroutine then you should be able to get rid of the local ncid_local copy. (I'm not sure if it's reliably safe to make a copy of the ncid variable like you do here, so this would probably be good to change.)
  • [OPTIONAL] You could use check_dim rather than check_var_or_dim; the latter is mainly intended for use by callers that could have either a dimension or variable name (such as the general-purpose find_var_on_file subroutine).
  • [OPTIONAL] If you can come up with more descriptive variable names than c and d, that would be preferable.

Copy link
Collaborator Author

@samsrabin samsrabin Feb 4, 2022 •
edited
Loading

Choose a reason for hiding this comment

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

Re: check_dim() (which I only see in river routing files): Wouldn't it fail when the dimension doesn't exist? I'm not sure whether the call of pio_inq_dimid() itself would cause a failure, but if not it'd result in dimlen /= value, right?

For reference, from components/mosart/src/riverroute/RtmIO.F90:

  subroutine check_dim(ncid, dimname, value)

    ! !DESCRIPTION:
    ! Validity check on dimension
    !
    ! !ARGUMENTS:
    implicit none
    type(file_desc_t),intent(in) :: ncid      ! PIO file handle
    character(len=*), intent(in) :: dimname   ! Dimension name
    integer, intent(in)          :: value     ! Expected dimension size
    ! !LOCAL VARIABLES:
    integer :: dimid, dimlen    ! temporaries
    integer :: status           ! error code
    character(len=*),parameter :: subname='check_dim' ! subroutine name
    !-----------------------------------------------------------------------

    status = pio_inq_dimid (ncid, trim(dimname), dimid)
    status = pio_inq_dimlen (ncid, dimid, dimlen)
    if (dimlen /= value) then
       write(iulog,*) subname//' ERROR: mismatch of input dimension ',dimlen, &
            ' with expected value ',value,' for variable ',trim(dimname)
       call shr_sys_abort()
    end if

  end subroutine check_dim

Copy link
Member

Choose a reason for hiding this comment

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

You're right that check_dim isn't called directly from anywhere else in CTSM, but it is what you want here; in CTSM, check_var_or_dim just delegates to either check_dim or check_var, and the implementation of check_dim differs from the one in RtmIO.F90:

CTSM/src/main/ncdio_pio.F90.in

Lines 304 to 348 in 8d9f988

!-----------------------------------------------------------------------
subroutine check_dim(ncid, dimname, dimexist)
!
! !DESCRIPTION:
! Check if dimension is on netcdf file
!
! !ARGUMENTS:
class(file_desc_t) , intent(inout) :: ncid ! PIO file descriptor
character(len=*) , intent(in) :: dimname ! dimension name
logical , intent(out) :: dimexist ! if this dimension exists or not
!
! !LOCAL VARIABLES:
integer :: dimid
character(len=*), parameter :: subname = 'check_dim'
!-----------------------------------------------------------------------
call ncd_inqdid(ncid, dimname, dimid, dimexist)
end subroutine check_dim
!-----------------------------------------------------------------------
subroutine check_var_or_dim(ncid, name, is_dim, exists)
!
! !DESCRIPTION:
! Check if variable or dimension is on netcdf file
!
! !ARGUMENTS:
class(file_desc_t) , intent(inout) :: ncid ! PIO file descriptor
character(len=*) , intent(in) :: name ! variable or dimension name to check
logical , intent(in) :: is_dim ! if true, check for dimension; if false, check for variable
logical , intent(out) :: exists ! whether the given variable or dimension exists on file
!
! !LOCAL VARIABLES:
character(len=*), parameter :: subname = 'check_var_or_dim'
!-----------------------------------------------------------------------
if (is_dim) then
call check_dim(ncid, name, exists)
else
call check_var(ncid, name, exists, print_err=.false.)
end if
end subroutine check_var_or_dim

While it's okay to be using check_var_or_dim here, it's just an unnecessary extra level of indirection.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Good point; I missed that in my search because I was only looking in files ending with .F90.

Okay, those three tasks are done:

  • d1ca586: ncid now inout in CallRestartvarDimOK().
  • 7601ec8: More descriptive variable names.
  • 948bcc9: CallRestartvarDimOK now uses check_dim() instead of check_var_or_dim().

And I've redone the merges, hopefully correctly this time.

Copy link
Member

Choose a reason for hiding this comment

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

Awesome, looks good; thanks for fixing those things; and the merge looks good now.


contains
! Public routines
Expand Down Expand Up @@ -177,6 +179,8 @@ end subroutine ReadNML
subroutine InitAllocate(this, bounds)
! !USES:
!
use clm_varpar, only : mxgrowseas, mxharvests
!
! !ARGUMENTS:
class(crop_type) , intent(inout) :: this
type(bounds_type), intent(in) :: bounds
Expand All @@ -191,22 +195,28 @@ subroutine InitAllocate(this, bounds)

allocate(this%nyrs_crop_active_patch(begp:endp)) ; this%nyrs_crop_active_patch(:) = 0
allocate(this%croplive_patch (begp:endp)) ; this%croplive_patch (:) = .false.
allocate(this%cropplant_patch(begp:endp)) ; this%cropplant_patch(:) = .false.
allocate(this%harvdate_patch (begp:endp)) ; this%harvdate_patch (:) = huge(1)
allocate(this%fertnitro_patch (begp:endp)) ; this%fertnitro_patch (:) = spval
allocate(this%gddplant_patch (begp:endp)) ; this%gddplant_patch (:) = spval
allocate(this%gddtsoi_patch (begp:endp)) ; this%gddtsoi_patch (:) = spval
allocate(this%vf_patch (begp:endp)) ; this%vf_patch (:) = 0.0_r8
allocate(this%cphase_patch (begp:endp)) ; this%cphase_patch (:) = 0.0_r8
allocate(this%latbaset_patch (begp:endp)) ; this%latbaset_patch (:) = spval
allocate(this%sdates_thisyr(begp:endp,1:mxgrowseas))
allocate(this%hdates_thisyr(begp:endp,1:mxharvests))
allocate(this%sowing_count(begp:endp)) ; this%sowing_count(:) = 0
allocate(this%harvest_count(begp:endp)) ; this%harvest_count(:) = 0

this%sdates_thisyr(:,:) = spval
this%hdates_thisyr(:,:) = spval

end subroutine InitAllocate

!-----------------------------------------------------------------------
subroutine InitHistory(this, bounds)
!
! !USES:
use histFileMod , only : hist_addfld1d
use histFileMod , only : hist_addfld1d, hist_addfld2d
!
! !ARGUMENTS:
class(crop_type), intent(inout) :: this
Expand Down Expand Up @@ -247,6 +257,16 @@ subroutine InitHistory(this, bounds)
ptr_patch=this%latbaset_patch, default='inactive')
end if

this%sdates_thisyr(begp:endp,:) = spval
call hist_addfld2d (fname='SDATES', units='day of year', type2d='mxgrowseas', &
avgflag='I', long_name='actual crop sowing dates; should only be output annually', &
ptr_patch=this%sdates_thisyr, default='inactive')

this%hdates_thisyr(begp:endp,:) = spval
call hist_addfld2d (fname='HDATES', units='day of year', type2d='mxharvests', &
avgflag='I', long_name='actual crop harvest dates; should only be output annually', &
ptr_patch=this%hdates_thisyr, default='inactive')

end subroutine InitHistory

subroutine InitCold(this, bounds)
Expand Down Expand Up @@ -468,31 +488,6 @@ subroutine Restart(this, bounds, ncid, flag)
end if
deallocate(temp1d)

allocate(temp1d(bounds%begp:bounds%endp))
if (flag == 'write') then
do p= bounds%begp,bounds%endp
if (this%cropplant_patch(p)) then
temp1d(p) = 1
else
temp1d(p) = 0
end if
end do
end if
call restartvar(ncid=ncid, flag=flag, varname='cropplant', xtype=ncd_log, &
dim1name='pft', &
long_name='Flag that crop is planted, but not harvested' , &
interpinic_flag='interp', readvar=readvar, data=temp1d)
if (flag == 'read') then
do p= bounds%begp,bounds%endp
if (temp1d(p) == 1) then
this%cropplant_patch(p) = .true.
else
this%cropplant_patch(p) = .false.
end if
end do
end if
deallocate(temp1d)

call restartvar(ncid=ncid, flag=flag, varname='harvdate', xtype=ncd_int, &
dim1name='pft', long_name='harvest date', units='jday', nvalid_range=(/1,366/), &
interpinic_flag='interp', readvar=readvar, data=this%harvdate_patch)
Expand Down
4 changes: 4 additions & 0 deletions src/main/clm_varpar.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,8 @@ module clm_varpar
integer, public, parameter :: dst_src_nbr = 3 ! number of size distns in src soil (BGC only)
integer, public, parameter :: sz_nbr = 200 ! number of sub-grid bins in large bin of dust size distribution (BGC only)
integer, public, parameter :: mxpft = 78 ! maximum number of PFT's for any mode;
integer, public, parameter :: mxgrowseas = 1 ! maximum number of crop growing seasons to begin in any year;
integer, public :: mxharvests ! maximum number of crop harvests in any year;
! FIX(RF,032414) might we set some of these automatically from reading pft-physiology?
integer, public, parameter :: nlayer = 3 ! number of VIC soil layer --Added by AWang
integer, public :: nlayert ! number of VIC soil layer + 3 lower thermal layers
Expand Down Expand Up @@ -142,6 +144,8 @@ subroutine clm_varpar_init(actual_maxsoil_patches, actual_numcft)
cft_lb = natpft_ub + 1
cft_ub = cft_lb + cft_size - 1

mxharvests = mxgrowseas + 1

! TODO(wjs, 2015-10-04, bugz 2227) Using actual_numcft in this 'max' gives a significant
! overestimate of max_patch_per_col when use_crop is true. This should be reworked -
! or, better, removed from the code entirely (because it is a maintenance problem, and
Expand Down
Loading