About stashCValues and year fraction handling

[pjbarjhoux]

Hi,

If I understand correctly, when rapid growth in the ocean-atmosphere carbon flux triggers the ocean component to enforce a sub-annual maximum step size, the ODE solver splits the year into several sub-intervals and calls SimpleNbox::stashCValues() once after each sub-interval converges.

The question concerns the stashCValues definition : each call adds the full annual LUC flux
(luc_e − luc_u) × veg_frac to cum_luc_va, without scaling by the sub-interval fraction yf :

hector/src/simpleNbox-runtime.cpp

Lines 392 to 393 in 662e6ca

	cum_luc_va =
	cum_luc_va + unitval((luc_e - luc_u) * c[SNBOX_VEG] / total, U_PGC);

I am wondering if this full annual flux should be multiplied by yf ?

  cum_luc_va =
      cum_luc_va + unitval((luc_e - luc_u) * c[SNBOX_VEG] / total * yf, U_PGC);

exactly as the biome-loop fluxes immediately below in the same function already do, eg:

fluxpool luc_fva_biome_flux =
    yf * veg_c[biome].flux_from_fluxpool(luc_e_untracked * veg_frac);

[pjbarjhoux]

In fact, since the pool fluxes are multiplied by the sub-step year fraction yf, the biome-level LUC updates are internally consistent.

If there is an inconsistency, it is local to the cum_luc_va variable, where an annual LUC term is added without yf when the solver has split the year into multiple substeps : current_luc_e & current_luc_u (used to compute cum_luc_va) are annual fluxes fixed in slowparameval() at the begining of each year.

So biome values are correct, but the forcing term (npp_luc_adjust) may become biased because it is driven by an (over)accumulated cum_luc_va.

From what I understand, it means that the model may compute lower NPP and higher atmospheric CO2 over time due to this behavior.

[pjbarjhoux]

Attached are the comparisons, by scenario, between version 3.5 and 3.5 after yf scaling of cumulative values of cum_luc_va