Annual sowing and harvest date outputs

[samsrabin]

Description of changes

Added annual outputs of sowing and harvest dates (SDATES and HDATES, respectively). This should simplify the determination of sowing and harvest date for postprocessing.

• Sowing dates are on new dimension mxgrowseas (maximum number of growing seasons allowed to begin in a year; currently hard-coded to 1).
• Harvest dates are on new dimension mxharvests (maximum number of harvests allowed in a year), which is mxgrowseas+1. This is needed because in year Y you might harvest a field that was planted in year Y-1, then plant and harvest again.
• The lengths of these dimensions are public constants of clm_varpar.

Additionally, I removed cropplant as discussed here.

The mxharvests concept enables the addition of more such outputs to further simplify crop postprocessing—for example, yield per growing season as a direct output rather than needing to cross-reference daily grain mass against the day of harvest.

Specific notes

CTSM Issues Fixed (include github issue #): #1537

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: None yet. Should I do the same tests as I did here?

[samsrabin]

Hmm, I only meant to include those last 2 commits… Not sure what the procedure is.

[ekluzek]

@samsrabin your whole branch will show up in the PR. If you only want a part of the branch, you'd need to make a separate branch with just those commits. Is there a problem with your whole branch coming to main?

[billsacks]

Thanks a lot for your work on this @samsrabin ! It seems like these new outputs will be really useful for analyzing crop planting & harvest. This looks to be implemented well. I have a few inline comments / questions.

[billsacks]

• 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?

[samsrabin]

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.

[billsacks]

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.

[samsrabin]

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?

[slevis-lmwg]

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.

[billsacks]

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.

[danicalombardozzi]

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.

[billsacks]

• [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.

[samsrabin]

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).

[billsacks]

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.

[samsrabin]

@ekluzek (#1616 (comment)):

your whole branch will show up in the PR. If you only want a part of the branch, you'd need to make a separate branch with just those commits. Is there a problem with your whole branch coming to main?

I looked over the changes and it's all what I had expected, except for the addition of .vscode/ to .gitignore. If you'd like, I can revert that.

[samsrabin]

In a test I just did, I'm noticing a bug where sdates_thisyr doesn't get assigned when the planting date is the first day of the simulation. Let's hold off on merging until I figure out what the deal is.

[samsrabin]

Okay, that's resolved now.

[billsacks]

Thank you very much for your careful work on this, and especially for thinking about what might be needed for backwards compatibility and other edge cases!

I have some comments / questions about this that I think would be easiest to talk through, so I'll send you a note about scheduling a meeting. Mainly I'd like to better understand what's needed solely for backwards compatibility and what's needed long-term, and also understand what pieces of this can be removed if we bring in #1623 first.

Other than that, my main comment relates to adding the new variables to the restart file to allow mid-year restarts. I think that's needed, unless you see a reason it's not that I may be missing. I'm happy to talk more about that as well.

[billsacks]

• 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.)

[samsrabin]

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.

[samsrabin]

Wrapping it like so seems to do the trick:

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

[billsacks]

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)

[billsacks]

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.

[billsacks]

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.

[samsrabin]

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

[billsacks]

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.

[samsrabin]

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.

[billsacks]

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

[billsacks]

• 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?

[billsacks]

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.

• Can you please add some details to this comment like, "This block of code is needed until we can rely on all restart files having been generated with Change CropPhenology to look at the day of year as of the start of the time step, not the end of the time step #1623 resolved"

[samsrabin]

Okay, done (ca90b78).

[billsacks]

Thanks!

[billsacks]

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.

[billsacks]

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.

• I am putting a place-holder checkbox here to remind us that this will need to be resolved in some way. I'm hoping that we can fix this in the course of addressing Change CropPhenology to look at the day of year as of the start of the time step, not the end of the time step #1623, probably by changing the conditional around the CropPhenology call to no longer check doalb. But if not, I think we'll need to come up with a different fix for this issue.

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.

[samsrabin]

Yes, I think that logic works.

[danicalombardozzi]

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

[billsacks]

@samsrabin I have opened #1628 with the fixes we discussed. When you get a chance, it would be great if you could merge that branch into this one (please don't commit/push that merge to this PR branch, though) and do any checks you feel are appropriate to confirm that crop phenology is working right with those changes. (I have tested that the changes in #1628 compile and run successfully, but haven't done any checks of crop phenology with those changes in place.) (Note that you'll probably get some merge conflicts because I changed one or two lines in #1628 that you have removed in this PR. Let me know if you want any help resolving the conflicts.)

Assuming it looks good to you, my plan would be to bring #1628 to master first, then bring in this PR shortly afterwards. (I'm hopeful that that plan would make this branch bit-for-bit, though I'm not sure about that.)

[billsacks]

Great sleuthing! Your explanation would explain why root_depth changes but this doesn't impact anything else: it sounds like root_depth is changing only at a time when it doesn't matter.

I would suggest making a branch off of master with just your proposed change to the root depth calculation. Then we can run the three tests that exercise the dynroots code:

ERS_D_Ld3.f19_g17.I1850Clm50BgcCrop.cheyenne_intel.clm-clm50dynroots
SMS_Lm1.f19_g17.I1850Clm50Bgc.cheyenne_intel.clm-clm50dynroots
SMS_Ly3_Mmpi-serial.1x1_numaIA.I2000Clm50BgcCropQianRs.izumi_intel.clm-clm50dynroots

and verify that the only thing that changes in these tests (relative to baselines) is ROOT_DEPTH. If so, then I think it's safe to change this as you suggest.

@danicalombardozzi I'd also welcome any thoughts you have on whether it's reasonable to change ROOT_DEPTH to be 0 in the non-growing season, as @samsrabin suggests.

Edit: as noted below (#1616 (comment)), let's also run a multi-year global test with dynroots to be more confident that this doesn't impact anything other than the ROOT_DEPTH diagnostic.

[danicalombardozzi]

Thanks for tracking this down @samsrabin! It does seem strange to have rooting depth other than 0 when crops are not alive. Your suggested solution (set root depth to 0 when crops are not live) seems reasonable to me. The only potential issue that comes to mind is whether the root depth affects the root growth when crops are planted. For example, if root depth is set to 0 (Sam's suggestion) compared to non-zero (current code), does that change the root depth when C is allocated to roots?

[billsacks]

Maybe we should run a longer global test to verify that @samsrabin 's proposed change only impacts root depth. I would assume that, if only the root depth diagnostic field changes but nothing else does, then it's safe to assume that this is an okay change - because presumably any changes to the evolution of root depth through the growing season would impact other fields as well. @danicalombardozzi and @samsrabin does this seem like a reasonable assumption? If so, in addition to the above three tests, we could add something like a 3-year coarse-resolution global test, like ERP_P72x2_Ly3.f10_f10_mg37.I2000Clm50BgcCrop.cheyenne_intel.clm-clm50dynroots.

[samsrabin]

Good idea. I just did ERP_P72x2_Ly3.f10_f10_mg37.I2000Clm50BgcCrop.cheyenne_intel.clm-clm50dynroots with the version that had been throwing diffs, and unexpectedly it passed. Does that test actually check root depth?

Path: /glade/scratch/samrabin/tests_0308-141505ch

[billsacks]

Thanks a lot for going ahead and doing that, @samsrabin ! I'm glad to see that you seem to be getting proficient with the test system!

I agree that that is a surprising result given the earlier comparison failure that we saw. I double-checked, and ROOT_DEPTH is indeed compared (you can see that by the inclusion of ROOT_DEPTH in the clm2.h0 cprnc.out files in /glade/scratch/samrabin/tests_0308-141505ch/ERP_P72x2_Ly3.f10_f10_mg37.I2000Clm50BgcCrop.cheyenne_intel.clm-clm50dynroots.C.0308-141505ch/run/). My guess is that this has to do with the fact that the test system only compares output from the final history file in the run, rather than comparing all history files. Could it be that ROOT_DEPTH differs early in the runs, but then comes back into agreement at some later point?

Anyway, I think it probably is not worth spending time trying to figure out why ROOT_DEPTH was the same when we expected it to differ: as long as your change (1) fixes the original problem with the ERS_D_Ld3.f19_g17.I1850Clm50BgcCrop.cheyenne_intel.clm-clm50dynroots test, and (2) we're not seeing any answer changes (beyond possibly ROOT_DEPTH) in the other tests, then I think it isn't super-important to fully understand what's going on here – unless your intuition is that it's worth digging into this more.

[samsrabin]

It looks like it is indeed passing because it's only looking at the last file. Other timesteps do differ in ROOT_DEPTH, but thankfully nothing else. (See ~samrabin/rootdepth_test_20220310.py and .log.) So I think we're good to go!

It's just the conflict that needs to be resolved, although I don't see that on my fork and I'm not sure how to fix it.

[billsacks]

Okay, awesome, thank you for your continued careful work to investigate this! I am satisfied with the ROOT_DEPTH fix based on what you have found. If I understand correctly, the fix you made for ROOT_DEPTH (setting it to 0 outside the growing season) should be moved onto this branch. Can you do that, or let me know if I'm misunderstanding?

Besides that, I wanted to make sure that the temporary fixes you put in place to reconcile the baseline differences indeed handle all of the differences we were seeing, or if there still might be some differences that we don't expect / understand (e.g., in the SSP test, which I don't understand well). Here was my earlier comment on this:

At this point, my feeling is: As long as this has given you confidence that the differences between your branch and the baseline are all due to oddities in the baseline rather than issues in the new code – i.e., if you're confident that your temporary code is all compensating for errors in the baseline rather than issues in the new code – then I'm happy to trust you on this, because I can tell that you have given this a lot of thought. In this case, I feel like we should rerun the full test suite with your temporary changes in place, to verify that this solves the baseline comparison issues for all of the tests. I can do that if you point me to a branch with these temporary changes, or I can tell you what to do – whichever you'd prefer.

Please let me know if you'd like me to run the full test suite with your temporary changes or if you'd like to do that yourself (or if you already have).

It's just the conflict that needs to be resolved, although I don't see that on my fork and I'm not sure how to fix it.

This indicates that there is a conflict that will appear when merging up to the latest version of master. I see this when I do a test merge of master into this branch; it is just some use statements that should be easy to resolve.

[samsrabin]

I am confident that the changes we're seeing are due to weirdness in the original code that this branch fixes, yeah.

I had been planning to make the "root depth is zero outside the growing season" change in a separate PR, but I can do it here if you prefer.

How about I do this:

1. On my feature branch (1537-crop-date-outputs3), set root depth to zero outside the growing season.
2. Merge master in to that branch and resolve the conflict.
3. Merge that branch in to my "temporary changes" branch (1537-crop-date-outputs3-ts).
4. Run the tests.

If all looks good (you might need to help with interpretation—or if you think it's easier, I'm happy for you to run the tests), then my feature branch will be ready for merge.

If you want me to run the full test suite, just let me know which tests it includes and I'll get it going.

[billsacks]

I am confident that the changes we're seeing are due to weirdness in the original code that this branch fixes, yeah.

Awesome.

I had been planning to make the "root depth is zero outside the growing season" change in a separate PR, but I can do it here if you prefer.

Normally I like to keep things separate, but that change is small enough (in code changes and impact) and related enough to the rest of this PR that I'd say go ahead and combine it.

Your plan sounds great. In step (4) (running tests on the temporary branch that puts in place changes to get bit-for-bit results), this will involve running the full aux_clm test suite on izumi and cheyenne according to the instructions in https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#notes-for-integrators . You won't be generating final baselines at this point, but I'd probably go ahead and generate baselines with some temporary name in case they're useful (you can always delete them shortly afterwards if they aren't useful).

Assuming that all looks good now, this will be essentially good to go (pending one more round of final testing). There are a couple of small changes I want to make: fixing the Fortran unit tests and adding your new output variables to one or more tests. I'll wait to make these changes until after you do step (4). Please let me know if it's okay for me to push to this branch – and also if you'd prefer to take a crack at those final steps yourself (in which case I'm happy to provide guidance).

[billsacks]

Okay, I have run the full aux_clm test suite on branch 1537-crop-date-outputs3-ts at commit 453a64bc6, with comparisons against ctsm5.1.dev083. Almost all baseline comparisons now pass. Two dynroots tests show differences just in ROOT_DEPTH, which I believe is expected.

There are now just these tests with unexplained differences:

ERS_Ly5_Mmpi-serial.1x1_smallvilleIA.I2000Clm50BgcCropQianRs.izumi_gnu.clm-ciso_monthly
SMS_Ly5_Mmpi-serial.1x1_smallvilleIA.IHistClm50BgcCropQianRs.izumi_gnu.clm-gregorian_cropMonthOutput
ERS_Ly6_Mmpi-serial.1x1_smallvilleIA.IHistClm50BgcCropQianRs.izumi_intel.clm-cropMonthOutput
SMS_D_Ly6_Mmpi-serial.1x1_smallvilleIA.IHistClm45BgcCropQianRs.izumi_intel.clm-cropMonthOutput

All of these are single-point cases using the artificial smallville dataset. The aspect of this dataset that I think is relevant here is that the surface dataset has a little bit of every possible crop. Some (many) crops still don't end up with any area in the run because we don't have parameters for those crops (so they get "collapsed" into other crops for which we do have parameters), but I think we still end up with some crop types in these tests that aren't present in our standard global tests. So I think these tests are flagging real differences that, while they don't show up in our typical runs, could show up if people included some additional, non-standard crop types.

I looked more carefully at the test ERS_Ly5_Mmpi-serial.1x1_smallvilleIA.I2000Clm50BgcCropQianRs.izumi_gnu.clm-ciso_monthly, and found that differences start in month 22 of the run (2001-10). (I haven't looked closely at the Hist – i.e., transient – tests, but based on the pattern of which tests pass and which fail, I'm guessing differences start somewhat later in the Hist tests – probably in year 4.)

I hate to say this, because I know it would be nice to be able to just wrap this up, but I think it's worth making sure that these remaining differences are also understood and expected. Please let me know if you'd like to talk through this to strategize about how to do this investigation. (Briefly, I would probably start with ERS_Ly5_Mmpi-serial.1x1_smallvilleIA.I2000Clm50BgcCropQianRs.izumi_gnu.clm-ciso_monthly, or a simpler SMS version of this test. I would probably look at the patch-level output (which should be present in the h1 file in that test) to see which crop(s) have differences.)

[samsrabin]

@billsacks Are those only failing on Izumi, or do they also fail on Cheyenne? I have Izumi access now, but I don't know anything about my permissions to / how to run there. Could you point me to where the output files are?

[billsacks]

We run different tests on izumi vs. cheyenne. I would strongly suspect that these tests would fail baseline comparisons on cheyenne, too, but you would need to create your own baselines for them. If you have izumi access, I think it should Just Work to run there; if not, let me know. But in the meantime, here are the output files from one of the tests (others are in a similar place):

/scratch/cluster/sacks/tests_0315-124122iz/ERS_Ly5_Mmpi-serial.1x1_smallvilleIA.I2000Clm50BgcCropQianRs.izumi_gnu.clm-ciso_monthly.GC.0315-124122iz_gnu/run

Here are the output files from my testing for the just-created ctsm5.1.dev085 – which isn't exactly the right baseline tag, but that tag is bit-for-bit with dev083, which is what your branch is up-to-date with:

/scratch/cluster/sacks/tests_0315-114032iz/ERS_Ly5_Mmpi-serial.1x1_smallvilleIA.I2000Clm50BgcCropQianRs.izumi_gnu.clm-ciso_monthly.GC.0315-114032iz_gnu/run

[samsrabin]

Looks like this is a problem for winter wheat only. I amended part of the reproduction-test code in 1537-crop-date-outputs3-ts to allow winter wheat planting in a calendar year even if it was alive at the beginning of the year, and it should now pass: commit 1b7ec6111.

Note: Nothing needs to be changed in the actual branch (1537-crop-date-outputs3).

[billsacks]

Awesome, thanks a lot for tracking that down. I looked at the changes you have made recently to the 1537-crop-date-outputs3-ts branch, and just wanted to check on one other change you made there: samsrabin/CTSM@e8299fdde. Is that change needed on the main branch? If not, is it needed in order to get bit-for-bit results with the baselines? Or was that just a temporary change that isn't actually needed?

[samsrabin]

@billsacks It's just needed for the bit-for-bit reproduction test. In the actual 1537-crop-date-outputs3 branch, idop is never ≤0.