c13/c14 isotope time-series is awkward for 1850 control (convert to streams file)

[ekluzek]

Erik Kluzek < erik > - 2017-11-10 16:15:08 -0700
Bugzilla Id: 2533
Bugzilla CC: cdkoven, dlawren, klindsay, rfisher,

Currently if you do a 1850 control, the assumed c13/c14 is for present day (delta of 0). You can get around this by asking for the c13/c14 time-series files, and if your start year is 0001, it'll use the 1850 value from the file. It's a little awkward since if you run a control for more than one thousand eight hundred and fifty years, it'll start running through the time-series file. For future scenarios it'll use the last value in the file, which is the most reasonable thing to do.

I think the best way around it, is to change the handling of this file from custom code to using the streams mechanism. This will actually shorten the amount of code needed in CLM as well as give more control, options, and flexibility in it's handling. The control for how the data was handled would be in the namelist for it and more obvious.

I think the sectors can be handled in streams by just representing them as three grid cells. This would also move the description of the data from baked inside the code to an outside grid file that could be changed with new datasets.

Definition of done:

• Create script to convert CMIP7 to streams format
• Include a mesh file for the C14 dataset since it has a latitude dependence
• Move datasets to inputdata
• rimport datasets @slevis up to this point, when we confirm that the file work correctly as streams
• Change code to use the new files @ekluzek

[ekluzek]

Erik Kluzek < erik > - 2017-11-10 16:16:48 -0700

This also means the default for 1850 and transient cases when use_c13/use_c14 are on should be to also turn on the time-series files for these cases. That'll do the right thing most of the time.

[ekluzek]

Keith Lindsay < klindsay > - 2017-11-13 09:38:34 -0700

Erik,

When I looked at using streams to handle the CMIP6 dataset, it looked to me like the shr code didn't handle the time units in the dataset.

time:units in the dataset is "years since 1850-01-01 0:0:0.0", and shr_string_parseCFtunit in shr_string_mod.F90 has the comment
! Input string is like "days since 0001-06-15 15:20:45.5 -6:00"
! - recognizes "days", "hours", "minutes", "seconds"
! - must have at least yyyy-mm-dd, hh:mm:ss.s is optional
! - expects a "since" in the string
! - ignores time zone part

Assuming I'm interpreting the above correctly, it seems like one of the following needs to happen:

1. extend shr code to handle this dataset's time:units
2. change the dataset to have time:units that shr code handles

Option 1) seemed unappealing to me because of the shr code's support for both no_leap and gregorian calendars. It seemed like it might be complicated to implement new functionality in shr code while not being consistent with the existing shr code.

If you go the route of option 2), I'll need to modify POP's code that processes this dataset, assuming that we want to continue having CLM and POP read the same file for this forcing. Please let me know if decide to pursue this path.

Keith

[ekluzek]

Keith Lindsay < klindsay > - 2017-11-13 09:59:10 -0700

(In reply to Erik Kluzek from comment #1)

This also means the default for 1850 and transient cases when
use_c13/use_c14 are on should be to also turn on the time-series files for
these cases. That'll do the right thing most of the time.

It sounds like this will lead to using the 1850 value for these forcings, in spinup runs, and in the 1850 control runs.

The 1850 values (permil) of Delta14C in the CMIP6 dataset in the different latitude bands are
90S-30S: -5.8
30S-30N: -4
30N-90N: -2.3

Is that what we want to be using in our BGC spinups and in the 1850 control runs? I'm used to using 0 permil for PI spinups, though this does create a jump at 1850.

I can see arguments for using 0 or for using the 1850 values.

I definitely think we should use the same value in the spinups as in the 1850 control runs. I also think we should use the same values in CLM and in POP.

Unfortunately, I don't see a specification for this in the C4MIP GMD paper, Jones et al. 2016.

[ekluzek]

charlie < cdkoven > - 2017-11-13 10:11:23 -0700

Agreed with Keith that there are arguments for using either 0 or the PI values during spinup and PI controls. I suppose I'd vote for using the PI values, even though that's not what my original code does, in the interest of avoiding artificial discontinuities.

[ekluzek]

Keith Lindsay < klindsay > - 2017-11-13 10:18:54 -0700

(In reply to Keith Lindsay from comment #2)

Option 1) seemed unappealing to me because of the shr code's support for
both no_leap and gregorian calendars. It seemed like it might be complicated
to implement new functionality in shr code while not being consistent with
the existing shr code.

I meant "... while being consistent ...".
The not slipped in from a different chain of thought.

[ekluzek]

David Lawrence < dlawren > - 2017-11-13 10:29:45 -0700

It makes more sense to me to use 1850 values for spinup and 1850 control. Question: if we were using zero per mil, then can you run with isotopes off during spinup and in the 1850 control? ... and only turn them on for the transient, thereby saving some computer time?

In the end, I think the decision on what to do should come from Charlie and Keith since they have more experience with isotopes.

[ekluzek]

charlie < cdkoven > - 2017-11-13 10:34:51 -0700

re Dave's Question: using zero or the PR value shouldn't really help either way for either isotope. For 14C, since the spinup problem is equilibration with respect to radioactive decay, you still have to do that (albeit in an accelerated way under AD). For 13C, the spinup problem there is with respect to isotopic fractionation, so, again, you still have to actually do it. The atmospheric value just becomes the boundary condition for the system in both cases, and in both cases you need to do the equilibration in AD mode in order to let it propagate rapidly through all the pools.

[ekluzek]

Erik Kluzek < erik > - 2017-11-13 11:11:58 -0700

(In reply to Keith Lindsay from comment #2)

Assuming I'm interpreting the above correctly, it seems like one of the
following needs to happen:

1. extend shr code to handle this dataset's time:units
2. change the dataset to have time:units that shr code handles

Option 1) seemed unappealing to me because of the shr code's support for
both no_leap and gregorian calendars. It seemed like it might be complicated
to implement new functionality in shr code while not being consistent with
the existing shr code.

I looked into this as well. The meaning is a bit awkward for the Gregorian calendar. And because of this the CF contention warns against using a unit of years. But, it says you define the interval as 365 for no_leap or 365.2425 for Gregorian. So how exact the end result would be depends on the reference date.

It wouldn't be hard to extend the code in cime/src/share/util/shr_string_mod to handle this. But, because of this issue we probably shouldn't. I'm bringing it up as an issue to CSEG, we'll see what feedback they have.

Here's the caution in the CF convention...

http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/build/ch04s04.html

We recommend that the unit year be used with caution. The Udunits
package defines a year to be exactly 365.242198781 days (the interval
between 2 successive > passages of the sun through vernal equinox). It
is not a calendar year. Udunits> > includes the following definitions for
years: a common_year is 365 days, a leap_year is 366 days, a Julian_year
is 365.25 days, and a Gregorian_year is 365.2425 days.

For similar reasons the unit month, which is defined in udunits.dat to be
exactly year/12, should also be used with caution.

[ekluzek]

Keith Lindsay < klindsay > - 2017-11-17 15:53:01 -0700

Regarding using 1850 values or 0 permil for 1850 control runs and spinup, I've sought input from Alex Jahn, who developed the C isotopes in POP, and Bette Otto-Bliesner and Esther Brady, as representatives of the paleo and isotope users. Everyone has either no preference or a slight preference for using 1850 values. So counting those three, Charlie and myself, no one has expressed a strong preference for using 0 permil.

So I think we've reached a consensus to use the 1850 values.

Note that this is independent of switching to shr code.

I do not plan on changing POP to use shr code, unless the dataset's time:units is changed. I will change POP to use the 1850 values in 1850 controls and spinup runs. (This is a little quirky for 14C because of the latitude bands in the dataset, but that's the way it is.)

[ekluzek]

Note, that this mechanism was inadvertently turned on in for 1850 control simulations. Since, the default is for 1850 conditions, it really didn't need to be. It only causes a problem if you run for more than 1850 years though. But, there were some simulations that ran for that long.

[ekluzek]

I've worked on this more and I want to list the issues with the current implementation. The starting year (1850) is hardcoded in, as well as that the time units will be in years. The number of sectors is hardcoded, as are the width of the bands. It will die if the number of sectors is wrong, but won't detect a change in the latitude band width. The time unit of years is considered a bad unit, because you can't detect a difference between Gregorian or no-leap calendars (like you can for other units). It does check the units, so will die if the starting year is wrong, or time isn't years, or the wrong units are used for C13 or C14. So that's good. But, the user also has no control over what years are used in the file, nor spatial nor temporal interpolation. All of these things come out of the box with a streams implementation. And it's trusted shared code, rather than custom code, so better from a maintenance standpoint. The streams implementation will automatically adjust to changes in both the temporal and spatial description of the file. The current mechanism also can't choose the time-period, which means the same file can't be used for Modern as for 20th Century. And as already shown this implementation has #592.

[ekluzek]

One other problem (with the current dataset not necessarily in the code itself) is that you can't do ncview on the file: $DIN_LOC_ROOT/lnd/clm2/isotopes/atm_delta_C14_CMIP6_3x1_global_1850-2015_yearly_v2.0_c171012.nc. This is because of the metadata on the "sector" variable is incorrect for CF. Using flag_* attributes allows you to view the file in ncview.