C isotope fields can have garbage values starting in ctsm5.1.dev008

[billsacks]

Brief summary of bug

Ever since ctsm5.1.dev008, the C isotope fields in this test ERP_D_Ld10_P36x2.f10_f10_musgs.IHistClm51BgcCropGs.cheyenne_intel.clm-ciso_decStart or its equivalent (it has changed names a couple of times since then) have been garbage, at least in some grid cells.

General bug information

CTSM version you are using: ctsm5.1.dev008 and later

Does this bug cause significantly incorrect results in the model's science? Yes, for C isotopes

Configurations affected: I haven't investigated how widespread this issue is.

Details of bug

@slevisconsulting ran across this issue (#1301 (comment)). I then ran cprnc on various pairs of tags for this test (ERP_D_Ld10_P36x2.f10_f10_musgs.IHistClm51BgcCropGs.cheyenne_intel.clm-ciso_decStart or its equivalent) to figure out where the problem is coming from. It looks like C13 and C14 values were reasonable up until ctsm5.1.dev007, but comparing ctsm5.1.dev007 and ctsm5.1.dev008 I see:

 RMS C13_AR                             Infinity            NORMALIZED    Infinity
 RMS C13_GPP                          1.3693E-07            NORMALIZED  5.9904E-01
 RMS C13_HR                             Infinity            NORMALIZED    Infinity
 RMS C13_NBP                            Infinity            NORMALIZED    Infinity
 RMS C13_TOTECOSYSC                     Infinity            NORMALIZED    Infinity
 RMS C13_TOTLITC                        Infinity            NORMALIZED    Infinity
 RMS C13_TOTSOMC                        Infinity            NORMALIZED    Infinity
 RMS C13_TOTVEGC                        Infinity            NORMALIZED    Infinity
 RMS C14_AR                             Infinity            NORMALIZED    Infinity
 RMS C14_GPP                          1.2535E-17            NORMALIZED  5.9652E-01
 RMS C14_HR                             Infinity            NORMALIZED    Infinity
 RMS C14_NBP                            Infinity            NORMALIZED    Infinity
 RMS C14_SOILC_vr                       Infinity            NORMALIZED    Infinity
 RMS C14_TOTECOSYSC                     Infinity            NORMALIZED    Infinity
 RMS C14_TOTLITC                        Infinity            NORMALIZED    Infinity
 RMS C14_TOTSOMC                        Infinity            NORMALIZED    Infinity
 RMS C14_TOTVEGC                        Infinity            NORMALIZED    Infinity

Ever since then, the average values of C isotope fields are enormous according to cprnc (e.g., 1e+229), and answer changes periodically lead to RMS diffs of Infinity.

Definition of done. Isotope cases should die when initial conditions don't have isotope values.

[billsacks]

We noticed that the initial conditions file (/glade/p/cesmdata/cseg/inputdata/lnd/clm2/initdata_map/clmi.I2000Clm50BgcCrop.2011-01-01.1.9x2.5_gx1v7_gl4_simyr2000_c190312.nc) doesn't appear to have C isotope variables on it. There is a bug (I think) with going from a non-ciso to ciso run, but I think that should just result in C isotope variables being at cold start values, not being infinity. Not sure if the problem is coming from that or from the dribble_crophrv_xsmrpool_2atm change made in dev008.

[wwieder]

@klindsay28 this is the land isotope bug.

[wwieder]

Naive question on this, Bill, but can we just start a new carbon isotope run from cold start to generate non-garbage values?

[billsacks]

Not a naive question at all. I'm not sure of the answer – I don't think I ever investigated it carefully. There is #67 that is specific to starting a ciso run from a non-ciso finidat when using init_interp. This issue (1358) seems to be different, though, in that we're getting total garbage values. But I'm not positive.

[billsacks]

Oh, and I see that I noted before (#1358 (comment)) that this issue might be from problems with initialization or might be from dribble_crophrv_xsmrpool_2atm.

[wwieder]

sorry, I wanted to bring this up again before too long. Is this something we can / should investigate? If so, the next question is who (@olyson?) and how should we proceed?

[billsacks]

Yes, I feel like this should be investigated. I tracked it down to the differences between ctsm5.1.dev007 and ctsm5.1.dev008, but didn't get as far as determining the actual cause for the difference in the ctsm5.1.dev008 tag. It would be great if someone could dig further into that to see what in that tag caused this bug. Besides being critical for the sake of having correct C isotopes for simulations, this bug also means that any unintended answer changes to C isotopes are currently going undetected.

[olyson]

I can take a look at it the next couple of weeks.

[olyson]

I ran the above referenced test with ctsm5.1.dev007 and found the isotopes values to be reasonable in the history files (as expected I believe).
That test used the initial file:

/glade/p/cesmdata/cseg/inputdata/lnd/clm2/initdata_map/clmi.I1850Clm50BgcCrop-ciso.1366-01-01.0.9x1.25_gx1v7_simyr1850_c200428.nc

I then ran the same test with the initial file that is used in ctsm5.1.dev008:

/glade/p/cesmdata/cseg/inputdata/lnd/clm2/initdata_map/clmi.I2000Clm50BgcCrop.2011-01-01.1.9x2.5_gx1v7_gl4_simyr2000_c190312.nc

and got garbage results in the isotope fields.
It would be interesting (but maybe not useful?) to find out why that new initial file produces garbage values, but maybe we could just do a new spinup to generate a new initial file as suggested by @wwieder ?

[billsacks]

Thank you @olyson !

[wwieder]

+1, thanks Keith

[olyson]

Just to double-check, I ran the test with ctsm5.1.dev008 using the initial file from ctsm5.1.dev007 and got reasonable isotope values.

[wwieder]

Just circling back to this issue. Does this just mean that we just need new initial conditions files, or is there a larger bug at play?

[billsacks]

From @olyson 's testing, it sounds like switching to initial conditions files with C isotopes will solve this particular issue. It may still be that there's a larger additional bug at play, though: I know that we don't expect things to work right when going from a non-ciso initial conditions file to a ciso case (see #67 ), but I thought that would just lead to C isotope values being at cold start values and so being out of sync with the bulk C... I didn't think that would actually lead to garbage values... but it may be that starting ciso values at cold start when the bulk C is already spun up will end up resulting in garbage, so maybe #67 is the entire explanation for what's going on here.

I think that for now we should go ahead with just using initial conditions with ciso and that will resolve this issue, then we should do more investigation in the context of #67 .

[billsacks]

It could also be helpful to know: do we avoid garbage C isotope values if we start with an initial conditions file without C isotopes in a run that does not use init_interp? (#67 should only impact runs that use init_interp.) This could be helpful to understand if there's another issue at play.

[billsacks]

@olyson will do the above test (#1358 (comment))

[olyson]

Using ctsm5.1.dev008, run the ERP test without isotopes, by setting use_c13, use_c14 to .false. in cime_config/testdefs/testmods_dirs/clm/ciso/user_nl_clm:
/glade/scratch/oleson/ERP_D_Ld10_P36x2.f10_f10_musgs.IHistClm51BgcCropGs.cheyenne_intel.clm-ciso_decStart.20230808_162834_9vt7tc/

Use the resulting restart file (/glade/scratch/oleson/ANALYSIS/ERP_D_Ld10_P36x2.f10_f10_musgs.IHistClm51BgcCropGs.cheyenne_intel.clm-ciso_decStart.20230808_162834_9vt7tc.clm2.r.2002-01-05-00000.nc) in a new ERP case with isotopes on and use_init_interp=.false., by adding the following to user_nl_clm:

finidat = '/glade/scratch/oleson/ANALYSIS/ERP_D_Ld10_P36x2.f10_f10_musgs.IHistClm51BgcCropGs.cheyenne_intel.clm-ciso_decStart.20230808_162834_9vt7tc.clm2.r.2002-01-05-00000.nc'
use_init_interp = .false.
check_finidat_year_consistency = .false.

The new case is:
/glade/scratch/oleson/ERP_D_Ld10_P36x2.f10_f10_musgs.IHistClm51BgcCropGs.cheyenne_intel.clm-ciso_decStart.20230808_165537_0f9d7x

The result is garbage isotope values, e.g., C13_NBP ranges from -1.12e+10 to 1.01e+17.

On the other hand, I wonder if this should be done a bit differently, as this new run is starting Dec 30 with an initial file of Jan 5. Suggestions welcome.

[billsacks]

Thanks a lot for doing this test, @olyson . I suppose it's possible that the problem in your latest case is due to a mismatch in the start time relative to the time of the initial conditions file... I would trust your intuition more than mine on whether that is likely to be a problem here, but my guess is that this indicates problems with the logic to initialize the carbon isotope values from the bulk carbon when using a restart file that doesn't contain C isotopes. There is logic like this that is supposed to apply to each field:

CTSM/src/biogeochem/CNVegCarbonStateType.F90

Lines 1616 to 1629 in f47ce5c

	if ( carbon_type == 'c13') then
	call restartvar(ncid=ncid, flag=flag, varname='leafc_13', xtype=ncd_double, &
	dim1name='pft', long_name='', units='', &
	interpinic_flag='interp', readvar=readvar, data=this%leafc_patch)
	if (flag=='read' .and. .not. readvar) then
	if ( masterproc ) write(iulog,*) 'initializing this%leafc with atmospheric c13 value'
	do i = bounds%begp,bounds%endp
	if (pftcon%c3psn(patch%itype(i)) == 1._r8) then
	this%leafc_patch(i) = c12_cnveg_carbonstate_inst%leafc_patch(i) * c3_r2
	else
	this%leafc_patch(i) = c12_cnveg_carbonstate_inst%leafc_patch(i) * c4_r2
	endif
	end do
	end if

My initial guess is that this logic is missing or buggy for one or more fields. Off hand, I can think of two ways to determine where the problem is:

(1) Examine which C isotope fields show issues in a setup like the one you did after just one or a few time steps.

(2) If (1) doesn't provide clear answers, we could try to bisect tags to find when this problem started. This might be a pain because we'd need to have compatible non-c-iso initial conditions files for each tag we run from, which could require making multiple non-c-iso initial conditions files like you did above.

I would personally be fine deferring this further investigation to a separate issue if you and others don't have time to dig into it more now: We at least know we can get reasonable C iso values if we have isotopes on the initial conditions file to begin with, and we can change our C isotope tests to use initial conditions with C isotopes on them. But we will need to solve this problem if we want to continue supporting going from non-C-iso runs to C-iso runs.

[olyson]

I dug into this a bit more but was unsuccessful at finding the problem. I took a look at the logic for all of the isotope fields. The ones that are in place look fine to me. There were three variables that did not seem to be initialized; leafc_storage_xfer_acc, storage_cdemand, and leafcmax. I added code to deal with those three but it didn't make a difference in the results.
It seems to take 6 or 7 timesteps for the values of the fields, e.g., C13_NBP, to start looking suspicious, but I've been unable to figure out what's causing that.
I hate to drop this but I don't really have time right now to pursue this further. It sounds like you are maybe suggesting opening a new issue to flag the specific issue of going from non-C-iso to C-iso runs...
To close the current issue, it looks like we need to provide an initial file from a historical run with isotopes on it.

[billsacks]

@olyson - thanks a lot for your continued investigation here. Yes, let's resolve this current issue just by pointing to initial conditions with C isotopes on them. I have opened #2119 to track the broader issue with going from an initial condition file without C isotopes to a run with C isotopes.

[olyson]

Fyi, the spunup isotope file below works fine with cesm2_3_alpha16b (ctsm5.1.dev130) for this test (no garbage values):
ERP_D_Ld10_P36x2.f10_f10_mg37.IHistClm51BgcCrop.cheyenne_intel.clm-ciso_decStart

/glade/p/cgd/tss/people/oleson/CLM5_restarts/ctsm51_cesm23a16bctsm51d130_1deg_GSWP3V1_1850pAD.clm2.r.0521-01-01-00000.nc

However, note that there are new variables on the restart file as of ctsm5.1.dev145, so NCAR/LMWG_dev#31 is generating a new spunup isotope file for later tags.

[wwieder]

It sounds like this is resolved. @olyson tried setting up a run that didn't have isotopes on IC files and the case failed.