Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Some uses of get_days_per_year should use the average number of days in a year, not the number of days in the current year #1612

Closed
billsacks opened this issue Jan 24, 2022 · 10 comments · Fixed by #1625
Assignees
Labels
bug something is working incorrectly

Comments

@billsacks
Copy link
Member

Brief summary of bug

There are a number of places in the code where a parameter that is expressed in units of something per-year is converted to units of per-second using the result of the get_days_per_year function. However, when running with a calendar that uses leap years, it looks like these should really be using the average number of days per year, not the number of days in the current year.

General bug information

CTSM version you are using: master

Does this bug cause significantly incorrect results in the model's science? No – just small errors

Configurations affected: BGC/CN runs with a Gregorian calendar

Details of bug

This emerged from my analysis documented in #1593 (comment) . The impact is minor, but this looks easy enough to fix that I'm going to go ahead and fix it.

Currently, for things like decay rates and mortality rates that are initially expressed in per-year units, the decay or mortality rate is slightly slower in a leap year than in a non-leap year – or, if the conversion happens just once at model initialization time, then the resulting rate would depend on whether the model starts in a leap year or a non-leap year. Even though this doesn't cause large errors, it's weird....

This impacts the following subroutines:

  • C14Decay: This should probably use a fixed, hard-coded constant, converting the decay rate from a per-year value to a per-second value: as currently written, C14 decays more slowly on leap years.
  • CNGapMortality: This should probably use a fixed, hard-coded constant, converting the mortality rate parameter from a per-year value to a per-second value: as currently written, mortality is a bit slower on leap years.
  • CNPhenologyMod.F90: CNEvergreenPhenology and CNStressDecidPhenology and the multiplier of leaf_long in CropPhenology: This should probably use a fixed, hard-coded constant, converting a litterfall rate from a per-year value to a per-second value: as currently written, litterfall happens more slowly on leap years.
  • decomp_rate_constants_bgc: As for the others, this should probably use a fixed, hard-coded constant
  • ndep_interp: This one is questionable, and may depend in part on whether the ndep input file and the model are consistent in their use of a GREGORIAN vs. NOLEAP calendar. But since I think this will typically come from model run with a NOLEAP calendar, I'm thinking that it should use a fixed, hard-coded constant (so that a constant annual ndep rate on the file would result in a constant annual ndep rate in CTSM, rather than a rate that varies slightly in leap vs. non-leap years). (But from spot-checking some recent files, it looks like rates are typically already specified as per-second, so this conversion from a per-year rate to a per-second rate is never done. I think that this implies that, in a GREGORIAN run, there is slightly more ndep than what occurred in the run that created the forcing data, since a year is slightly longer.)
    • Actually, upon further reflection, I think it's best not to change this one.

My plan is to introduce a new get_average_days_per_year function that returns the average number of days per year for the given calendar (which can be obtained via a query to an appropriate ESMF calendar routine), and then use that function in the above subroutines rather than using the number of days in the current year.

@billsacks billsacks added the bug something is working incorrectly label Jan 24, 2022
@billsacks billsacks self-assigned this Jan 24, 2022
@ekluzek
Copy link
Collaborator

ekluzek commented Jan 24, 2022

This makes sense to me. I like your solution of average_days_per_year, I see how you get that from an ESMF_CalendarGet call.

@billsacks
Copy link
Member Author

billsacks commented Jan 26, 2022

Although I suspect this is generally a minor issue, it could potentially be important for exact restart in the following scenario: If the conversion from per-year to per-second happens in model initialization, currently the resulting parameter will differ depending on whether the model is starting on a leap year or non-leap year. So then answers would differ if (say) you started the model in year 2000 and ran for 2 years in a row, vs. starting in 2000, running for a year, restarting in 2001 and running for another year (in the former, you would use the same parameter value for both years; in the latter case, the parameter value for the second year would differ from that used in the first year). I'm not sure if this is currently the situation for any parameters, but I think it's worth fixing this problem to avoid this situation arising in the future.

@billsacks
Copy link
Member Author

Unfortunately, ESMF_CalendarGet isn't available in the esmf_wrf_timemgr. I realize that is only needed in mct cases, but it seems like we still need to support that usage for now. I feel like it would be useful to get this fix in place while it's still fresh in my mind, so for now I propose to have the implementation in the clm time manager use hard-coded values dependent on the calendar (since we currently only support two calendars, this shouldn't be too bad). Then I'll open an issue to clean this up once we can rely on having the real esmf library available.

@ekluzek let me know if you object to that plan.

@ekluzek
Copy link
Collaborator

ekluzek commented Jan 26, 2022

Your solution makes sense to me. This does clean it up significantly from what what is was before.

An alternative would be to require USE_ESMF_LIB=TRUE even for MCT cases right now. And maybe that would be a good way to go. Then you wouldn't have to change this later. And we could add a check for that in buildlib to tell you that you are required to set USER_ESMF_LIB==TRUE. I think I'd prefer this option and actually think this might be easier for you to do as well. But, I think either should be fine.

@billsacks
Copy link
Member Author

I don't think CTSM has the freedom to require USE_ESMF_LIB=TRUE for MCT cases: that seems like a CESM-wide decision, and it seems like that could be challenging for support on other machines for people who don't yet have ESMF available to them and so still want to run with the mct driver. Maybe I'm missing something....

@ekluzek
Copy link
Collaborator

ekluzek commented Jan 26, 2022

We can bring it up to the CESM level, it would be good to discuss at that level. But, also note, that WACCM currently requires ESMF already, so I think some components do take this choice separately. Since, it's already required by at least one component we wouldn't be breaking existing standards. And also it is a hard requirement for NUOPC, which is the future, so I don't see that breaking it now for CTSM is that big of a deal for a future where it's required. We would just be moving the requirement a little bit ahead of the required schedule -- but not that much really.

Any machines that don't have ESMF available -- are going to need ESMF available with the switch to NUOPC. So it actually would be good to find out if this is a problem, by requiring it. We would find out if there are machines that this is a problem for and be able to start addressing it. In the meantime they can use older versions, or we could give them some small code changes to do it. It's possible there really aren't any such machines that are going to be a problem, but they will be a problem with the eventually requirement to use NUOPC anyway. So it's just upping the time-line rather than really adding any new requirements.

@ekluzek
Copy link
Collaborator

ekluzek commented Jan 26, 2022

Although -- back peddling a bit -- this will require us to have CAM and CESM change MCT testing for their test lists using CTSM to require this. So you are right we need to talk to CAM and CESM about this....

@billsacks
Copy link
Member Author

You raise good points, but I'd suggest that this decision not be driven by this particular issue / change. I do see the value long-term in having more general code, but right now our time manager only supports noleap and gregorian calendars anyway, so more extensive changes would be needed if we wanted to support additional calendars – and as long as we are still just supporting those two calendars, hard-coding the result feels like a reasonably safe option. What I can do to make it even more robust is to compare the hard-coded value with the result of get_curr_days_per_year and abort if the difference is more than 1 (just in case someone hacked the esmf calendar code to model some other planet, for example, before we update the code to be more robust).

@ekluzek
Copy link
Collaborator

ekluzek commented Jan 26, 2022

Again, I'll let you decide. I agree it's relatively safe to keep it hardcoded. We do know that people are starting from CESM to model other planets though.

I did look for MCT tests in CESM and CAM. CESM only has three tests, so it wouldn't be a big deal for them to add USE_ESMF_LIB=TRUE to shell_commands for the default tests. CAM is mostly testing with MCT though, so it would be a bigger deal there. But, still they only test on cheyenne, izumi, and hobart which we know has ESMF support. And CESM MCT tests are only cheyenne, bluewaters and izumi (and they run NUOPC on bluewaters) so we know those machines are fine as well.

For our own testlist lawrencium-lr3 would be a good one to check on. I'll add that to tomorrows discussion.

@billsacks
Copy link
Member Author

Well it turns out we need the uglier solution for now anyway: When I run with the real esmf library (via a nuopc case), with a GREGORIAN calendar, ESMF_CalendarGet tells me that the calendar has exactly 365 days per year, even though it confirms that it's a Gregorian calendar. I'll contact ESMF support about that, but in the meantime, it seems like using our own, hard-coded values is the way to go.

billsacks added a commit that referenced this issue Feb 24, 2022
Use avg days per year when converting param units

When converting parameter units from per-year to per-second, use average
days per year instead of current number of days per year. This is
relevant when running with a Gregorian calendar.

See #1612 for details.

Resolves #1612 (Some uses of get_days_per_year should use the
average number of days in a year, not the number of days in the current
year)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug something is working incorrectly
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants