[FEATURE REQUEST] Continued cleanup of user rate functions and code in the util/ folder

[yantosca]

This issue is a place where we can discuss further updates to cleaning up the rate-law functions that get inlined into the Rates module. A couple of action items have been identified by @RolfSander and @jimmielin

1. We still use -B0/TEMP in ARR and +B0/TEMP in ARR2. This is a terrible inconsistency. It won't be easy to fix this without backward incompatibilities, though. a. NOTE: It seems that none of the sample mechanisms use ARR2. If this is the case we can maybe remove it.

2. What is going to happen with ARR_ac, ARR_ab, etc., as mentioned in @jimmielin's paper? Let's try to find a consistent structure for all these ARR* functions.

[yantosca]

I have pushed a cleanup_util branch where we can add modifications.

[RolfSander]

I agree: Let's delete ARR2.

[RolfSander]

More questions:

• What about EP2, EP3, and FALL? Is anyone actually using them?

• Even the _sp variants of the functions return dp values. This looks strange. Would it be possible to return KPP_REAL, depending on the setting of the KPP command #DOUBLE ON/OFF?

• Some functions use temp as a dummy parameter, some use temp from _kpp_Global. Can we get some consistency here?

[yantosca]

1. FALL, EP2, and EP3 are used in models/saprc99.eqn and models/saprcnov.eqn files.

2. Also do we want to replace all of the e exponents with d exponents in the sample model equation files?

3. I can modify the functions to return KPP_REAL for the output type. That's not a problem.

4. I agree with you about temp. Not sure why that takes it as an argument when it's a global variable. I left it as is because I didn't know if that would break anything.

5. It also looks like none of the k_3rd* or k_arr* rate law functions are used either. We could probably remove those too. For GEOS-Chem we don't rely on the included user rate law functions, as we provide our own. Would that break anything in your MECCA code?

[RolfSander]

For GEOS-Chem we don't rely on the included user rate law functions, as we provide our own. Would that break anything in your MECCA code?

In addition to some functions from UserRateLaws.f90, MECCA also uses its own rate law functions from a separate file. We use the following KPP command for this:

#INLINE F90_RCONST
  USE messy_main_tools_kinetics
#ENDINLINE

What is our general strategy? Do we want to offer all kinds of possible functions in UserRateLaws.f90? Or do we prefer that the KPP users write their own functions?

It also looks like none of the k_3rd or k_arr rate law functions are used

They are used by MECCA. However, it would be no problem for me to move them out of UserRateLaws.f90 and into our MECCA-specific messy_main_tools_kinetics.f90.

FALL, EP2, and EP3 are used in models/saprc99.eqn and models/saprcnov.eqn files.

Okay, let's keep FALL, EP2, and EP3.

I agree with you about temp. Not sure why that takes it as an argument when it's a global variable.

When MECCA is run inside our global model, the code is vectorized and temp becomes a 1D-array. This works fine when we supply temp(i) as a function parameter but I'm not sure if we could take temp from _kpp_Global.

Also do we want to replace all of the e exponents with d exponents in the sample model equation files?

I'm already happy if rate constants have an uncertaincy of less than about 10 %. Defining them as dp instead of sp doesn't make much sense from a numerical point of view. The EXP function for the temperature dependence is computationally quite expensive. I assume that it uses much more CPU time in dp than in sp. This is a good reason for using sp inside SUBROUTINE Update_RCONST.

However, this mix of dp and sp creates so much trouble that it would be much simpler if everything is dp. This could be done by replacing the e exponents with d exponents, as you suggest.

So, the question probably is: Is the EXP function in dp sufficiently expensive that we want to go through the trouble of mixed precision variables inside SUBROUTINE Update_RCONST?

I can modify the functions to return KPP_REAL for the output type. That's not a problem.

Thanks for the offer. However, I'm not so sure if this is really so easy. Depending on whether KPP_REAL=sp or KPP_REAL=dp, we will need the DBLE() commands in ARR_sp or not.

[yantosca]

Thanks @RolfSander. I think we would maybe only want to keep the rate law functions that are needed for the out-of-the-box builds of the mechanisms that ship with KPP (FALL, ARR, EP2, EP3), and to have users specify their other rate law functions. That certainly would be simplest.

In my experience with rewriting the GEOS-Chem rate-law expressions, I find that dp precision is needed internally to avoid underflow. This is certainly true in heterogeneous chemistry, where you have a lot of EXP and LOG statements. Also, in the saprc99 and saprcnov equations, there is a rate constant of 1e-54, which would produce an underflow since that number can't be represented in sp.

Ah, I didn't know that about how k_arr* functions worked in MECCA. I think it'd be good if you can port those to an inlcude file outside of KPP.

In another (local) branch I've started splitting up the ARR into ARR_abc, ARR_ab, and ARR_ac functions, to avoid calling EXP for terms that would go to 1. Will keep you posted as to how that goes.

[jimmielin]

In addition to some functions from UserRateLaws.f90, MECCA also uses its own rate law functions from a separate file. We use the following KPP command for this:

#INLINE F90_RCONST
  USE messy_main_tools_kinetics
#ENDINLINE

What is our general strategy? Do we want to offer all kinds of possible functions in UserRateLaws.f90? Or do we prefer that the KPP users write their own functions?

Hi @RolfSander, this question came up when I was working on the manuscript as well. I was going to add that we are now offering k_3rd_iupac and k_3rd_jpl, but emphasizing that new rate laws were added may give users the confusion that rate laws need to be bundled into KPP to be supported (which is certainly not true!)

So if the strategy is that users should provide their own (as in MECCA using F90_INLINE and as in GEOS-Chem where a lot of rates come from the GEOS-Chem-to-KPP interface and updated externally), I'll mention this in the manuscript in addition to the new rate laws, to clear confusion. Does this sound good?

Thanks!

[RolfSander]

I think we would maybe only want to keep the rate law functions that are needed for the out-of-the-box builds of the mechanisms that ship with KPP (FALL, ARR, EP2, EP3), and to have users specify their other rate law functions.

Yes, fine with me.

I noticed that GEOS-Chem also uses F90_RCONST to make the additional functions available to KPP:

#INLINE F90_RCONST
  USE Hg_RateLawFuncs
#ENDINLINE

I think we should explain in the user manual how F90_RCONST can be used for this purpose. It is important that the USE command in F90_RCONST comes first, before any executable f90 code.

In my experience with rewriting the GEOS-Chem rate-law expressions, I find that dp precision is needed internally to avoid underflow. This is certainly true in heterogeneous chemistry, where you have a lot of EXP and LOG statements. Also, in the saprc99 and saprcnov equations, there is a rate constant of 1e-54, which would produce an underflow since that number can't be represented in sp.

Good point. The underflow is another reason why sp is problematic.

If we say that we want dp everywhere even though it needs more CPU time, then it makes sense to change the e exponents with d exponents in the sample model equation files.

Ah, I didn't know that about how k_arr* functions worked in MECCA. I think it'd be good if you can port those to an inlcude file outside of KPP.

Yes, if we are going to keep only FALL, ARR, EP2, EP3, I will remove these functions from UserRateLaws.f90.

[RolfSander]

@jimmielin: To adjust ARR and k_3rd in your manuscript, it's probably best to wait until we have merged cleanup_util into dev. At the moment we are still discussing some minor details...

@yantosca: Thanks for adding the ARR_[abc] functions! I think we need to mention in the user manual that the functions ARR and ARR2 don't exist anymore in KPP-3.0.0.

[yantosca]

@RolfSander: Yes, I'll push a commit to update the docs in the cleanup_util branch and make sure that all instances of ARR are updated with the new functions. Then we can resolve any conflicts with @jimmielin's updates in the autoreduce branch.

[RolfSander]

All functions in UserRateLaws.f90 consistently use temp from kpp_Global now. This is good! However, I don't think we should remove the functions for 3rd-order reactions completely because the formulas from JPL and IUPAC are widely used. We could add new versions of them under slightly different names, e.g., k3rd_iupac, k3rd_jpl and k3rd_jpl_activation. For consistency with ARR_*, these new functions should also take temp from kpp_Global. If you agree, I will take care of this.

[yantosca]

That's fine with me!

[yantosca]

We can close out this thread as PR #52 is merged.