Are we using Update_RCONST() correctly?

[RolfSander]

Hello,

The switch ICNTRL(15) allows us to call Update_RCONST() several times during the KPP sub-timesteps. Some of the rate constants depend on other concentrations. In GEOS-Chem for example, k(IONO+BrSALA) depends on C(ind_IONO).

Why is Update_RCONST() using concentrations from the C array, which still has the old values before the integration started? Don't we have to use the intermediate values from the YNEW arrray inside Update_RCONST()?

[RolfSander]

@adrian-sandu: Can you help us with this question?

[adrian-sandu]

Hello,

This is from the legacy code design (one set of coefficients is fixed during each time step, except for the time dependent photolysis rates).

Perhaps we need a version Update_RCONST(Y) that takes current concentrations as argument, with double precision,optional :: Y.

[RolfSander]

Hello Adrian,

Thanks for confirming that we should use Y instead of C!

A colleague of mine (Simon Rosanka) has already written some code to tackle this problem. It is currently implemented in the MECCA chemistry module. After testing, I hope that we can transfer the changes to the KPP repository so that they will be available in the next KPP release.

[RolfSander]

Hello @adrian-sandu,

I have a related question: Is it correct that calling Update_RCONST() several times during the KPP sub-timesteps only makes sense in the non-autonomous mode?

If that is true, KPP should print a warning if a user sets the questionable combination of ICNTRL(1)=1 and ICNTRL(15)>0.

[RolfSander]

I've added @srosanka as an "Outside Collaborator" to our KPP repo. He has written a patch to deal with this Update_RCONST problem, and he offered to commit his code into a new branch here.

[yantosca]

Thanks @RolfSander @adrian-sandu @srosanka for looking into this. FYI, in GEOS-Chem, we only update rate constants before calling the integrator (otherwise it would be too time consuming) by setting ICNTRL(15) = -1. So this update wouldn't affect GEOS-Chem.

We can also get a sense of what impact this change would have by looking at the output of the C-I tests.

[srosanka]

Hi, I just pushed the required changed to the issue93 branch. In this commit I also updated all Rosenbrock integrator to take this modifications into account. Let me know if you agree to these changes and I will open a related pull request.

[srosanka]

@yantosca That is good to know. As far as I recall, the default option in MESSy also only updates all rate constants before calling the integrator. I think that this is a valid approach for most gas phase mechanisms, especially when considering the performance improvement in global model applications. However, when we implemented aqueous phase chemistry in deliquescent aerosols in MESSy (AERCHEM), we noticed that not updating rate constants that depend on variable species during the integration can significantly influence the predicted results.

[adrian-sandu]

Hi Bob, Rolf, Simon,

It makes sense to call Update_RCONST every time we call the function/Jacobian for non-autonomous systems. Perhaps very important for wet chemistry. The shortcut with calling it every operator split time step works well for gas phase, and it saves a ton of computations in 3D. So keeping both options is good.

[yantosca]

Thanks @adrian-sandu. I think so too!

[yantosca]

We can now close this issue now that PR #106 has been merged.