KPP export working

[stulacy]

This is a high priority feature and will be worked on once the new release is live. This in theory should just require a few modifications to the code that runs on the old site.

[stulacy]

A basic KPP export was added in #210 , although there are several outstanding issues:

• [ ] When upgrading to KPP-3.0.0, the line #INCLUDE atoms must be changed to #INCLUDE atoms.kpp, but this could break functionality for KPP 2 users
• [ ] M, N2, O2, RO2, H2O are inlined, but should they be included in #DEFVAR?
• [ ] "USE constants" and "CALL mcm_constants" have been removed as there is no evidence they are being used, but it would be helpful to double check this in practice
• [ ] One note, using E exponents will force REAL4 precision. Better to use e.g. 250.0_dp, where dp is the KPP double precision kind parameter.
• [ ] H2O is now declared as a species in the #DEFVAR section, but KPP ignores it because it does not occur in any equation. This is resolved by adding it as a product of reaction <19> HO2 + OH = H2O + O2, however, this requires the user to export inorganic reactions and means we can't use the generic PRODS. A more general solution would be nice

[RolfSander]

Hello Stuart,

Thanks for all your work of implementing the KPP export! I have made a quick test on the new MCM web page, and it works fine.

My next plans are:

• I will create a minimum working example which potential users can download for testing the MCM-generated code with KPP.

• I will discuss with the other KPP developers how we tackle the remaining problems.

• I will tell my colleagues about the new export function and solicit feedback.

I think we can finish these tasks until early next year, maybe January.

[RolfSander]

I just noticed that the MCM home page needs to be updated as well. It currently says:

(currently only FACSIMILE is supported, although KPP will be added in an upcoming release).

[stulacy]

Good spot Rolf, I've updated the home page accordingly.

Your plan sounds sensible, I look forward to working with you to polish off this functionality. Please update this issue when you have an idea for the next development phase.

[RolfSander]

Hello Stuart,

I have now created a minimum working example (MWE) which shows how to use the MCM-generated files for KPP. You can find it in our KPP repository in the mcm branch:

https://github.com/KineticPreProcessor/KPP/tree/mcm/examples/mcm

We have also discussed the most flexible file structure for external input. We now suggest to move all external input out of mcm_isoprene.eqn and into a separate file called constants_mcm.f90. The new file constants_mcm.f90 is a static file that is the same for all extracted mechanisms. It only needs to be adjusted when there are any changes in the MCM Generic Rate Parameters (https://mcm.york.ac.uk/MCM/rates/generic, https://mcm.york.ac.uk/MCM/rates/complex, or https://mcm.york.ac.uk/MCM/rates/photolysis). Ideally, mcm_isoprene.eqn and constants_mcm.f90 should look like this:

https://github.com/KineticPreProcessor/KPP/blob/mcm/examples/mcm/mcm_isoprene.eqn

https://github.com/KineticPreProcessor/KPP/blob/mcm/examples/mcm/constants_mcm.f90

Would you be able to generate mcm_isoprene.eqn automatically in this format with your export function?

With respect to our open questions, the following approaches have been implemented in the MWE:

1) Although previously USE constants and CALL mcm_constants were dysfunctional, the idea behind them is probably still be the best approach. Minor name changes were necessary, and we now only have these two lines inside the F90_RCONST block in mcm_isoprene.eqn:

  #INLINE F90_RCONST
    USE constants_mcm
    CALL define_constants_mcm
  #ENDINLINE {above lines go into the SUBROUTINES UPDATE_RCONST and UPDATE_PHOTO}

The F90_GLOBAL block can be deleted completely from mcm_isoprene.eqn.

The variables M, N2, O2, RO2, H2O and zenith are now all declared in constants_mcm.f90. It is up to the user to assign meaningful values to them (except for RO2 which is already defined).

2) H2O and PROD:

Indeed, we can never be sure if a user has H2O as a chemical species in their mechanism or not. Even though reaction <19> produces H2O, the user may not export these inorganic reactions. Therefore, it is probably best to leave H2O in the rate constants unchanged, i.e., do not convert it to C(ind_H2O). Instead, it is up to the user to define H2O properly. Either they set H2O = C(ind_H2O), or they define it in another way, e.g., H2O = myfunc(humidity).

Reactions <2> and <19> should simply produce PROD.

3) Mapping of photolysis rate constants (j-values)

All solutions discussed so far are problematic. Even if the current MCM numbering remains unchanged in the future, the mapping will be different for different mechanisms. Let's say that a user downloads the mechanism for CH4 first, but later they decide that they want to have isoprene as well. When adding the isoprene reactions, they will have the same index in the J array being used for different photolysis reactions.

We need a solution that works out-of-the-box for the MWE but can be adjusted for improved efficiency, e.g., for global atmospheric chemistry models. I think we can achieve this with named constants in a static file that contains all j values. I have now assigned a name to each j values, you can see the names in constants_mcm.f90. What do you think about this approach?

4) Don't worry about #INCLUDE atoms.kpp. For backward compatibility, I have created a link from atoms to atoms.kpp in the KPP repository.

5) My last point is pretty unimportant, it's just my personal preference: The header section of mcm_isoprene.eqn contains 48 lines which are written as a KPP comment, using the curly braces syntax with {...}. When looking at individual lines, it may not be obvious that we're inside a comment block. Note that KPP can also define comments in preprocessor style: All lines starting with // are treated as comments. Maybe that's a better way to define a large number of comment lines?

[stulacy]

Hi @RolfSander

Thanks for getting back to me with these suggestions, they largely look straight forward to implement and I'm happy to change the block comment notation to the preprocessor style for the header.

Having the complete set of photolysis rates makes sense, but would you mind if I changed the constant names to be the MCM J numbers? I.e. J_56 = 34 rather than J_NOA = 34. This would be much easier to automate on our side, and we can assure that those J numbers don't change in the future.

[RolfSander]

Hello @stulacy,

Thanks for getting back to me with these suggestions, they largely look straight forward to implement and I'm happy to change the block comment notation to the preprocessor style for the header.

Thanks!

Having the complete set of photolysis rates makes sense, but would you mind if I changed the constant names to be the MCM J numbers? I.e. J_56 = 34 rather than J_NOA = 34. This would be much easier to automate on our side, and we can assure that those J numbers don't change in the future.

If you prefer J_56 over J_NOA, that's perfectly fine. Technically (memory usage for the J array etc.) this should make no difference. My idea behind the named constants was that chemists could understand the meaning of J_NOA directly, while looking up J_56 in a table is more of an effort.

Anyway, the only really important point is that the J numbers don't change in the future.

[stulacy]

We can go with the human readable J rates, that's fine.

I've just noticed something else that might be a potential issue. In the constants_mcm.f90 file, the RO2 sum is using the peroxy radicals that are in the isoprene mechanism. However, these peroxies are specific to the submechanism being exported, so should they be included within the global constants file, or within the exported .eqn?

If it was kept within the constants then it should comprise all peroxies within the full MCM, but (from my limited understanding of KPP), that will cause problems when some of the C(ind_x) values aren't in the mechanism. Or am I misunderstanding?

[RolfSander]

We can go with the human readable J rates, that's fine.

Thanks!

In the constants_mcm.f90 file, the RO2 sum is using the peroxy radicals that are in the isoprene mechanism. However, these peroxies are specific to the submechanism being exported, so should they be included within the global constants file, or within the exported .eqn? If it was kept within the constants then it should comprise all peroxies within the full MCM, but (from my limited understanding of KPP), that will cause problems when some of the C(ind_x) values aren't in the mechanism. Or am I misunderstanding?

Good point, thanks for spotting this! Yes, you're right, having unused C(ind_x) would cause problems. The definition of RO2 has to be moved back into the exported .eqn file.

The best place is just before CALL define_constants_mcm inside the F90_RCONST block, i.e.:

#INLINE F90_RCONST
  USE constants_mcm
  RO2 = C(ind_CH3O2) + C(ind_C51O2) + C(ind_CH3COCH2O2) + C(ind_HOCH2CO3) + &
    C(ind_CO2C3CO3) + C(ind_IPRHOCO3) + C(ind_BIACETO2) + C(ind_CH3CO3) + &
    [...]
  CALL define_constants_mcm
#ENDINLINE {above lines go into the SUBROUTINES UPDATE_RCONST and UPDATE_PHOTO}

[stulacy]

I've updated the web-app to generate both mcm_export.eqn and constants_mcm.f90 automatically. Do they look suitable to you @RolfSander ?

[RolfSander]

Hello @stulacy,

Thanks for the updated implementation!

I had to make a few changes to the files, and then I was able to reproduce exactly the same results with my MWE as before. I think all of the issues should be very easy to fix:

1) The letter R is missing in SUBOUTINE:

   END SUBOUTINE define_constants_mcm

2) You have moved the comment about peroxy radicals into the F90_RCONST block. This is good because that's where it should be. However, please note that this block is transfered verbatim into the KPP-generated Fortran files. Therefore, comments must follow Fortran syntax, not KPP syntax. Could you change the lines to:

  ! Peroxy radicals
  ! WARNING: The following species do not have SMILES strings in the database. 
  !          If any of these are peroxy radicals the RO2 sum will be wrong! 
  ! O O3 NO NO2 NO3 O1D N2O5 OH HO2 H2 CO H2O2 HONO HNO3 HO2NO2 SO2 SO3 HSO3 NA
  ! SA CL

3) In constants_mcm.f90, there are now duplicate declarations of KRO2HO2, KRO2NO and KRO2NO3. These need to be removed.

4) You have added the new complex rate constant K16ISOM which depends on HO2, NO and NO3. I think we can treat these species in the same way as H2O. We declare the variables, and if a user needs it, they can assign suitable values themselves. Thus, the line

  REAL(dp) :: M, N2, O2, RO2, H2O, zenith

needs to be changed to:

  REAL(dp) :: M, N2, O2, RO2, H2O, HO2, NO, NO3, zenith

5) Finally, I encountered a problem with the preprocessor style of KPP comments. After some time, I figured out that KPP needs a space after the //, otherwise KPP produces an error. This seems to be a bug in KPP which we need to fix. For now, it would help us if you could add a space after the two slashes in the comments.

[stulacy]

Cheers Rolf. The syntax/typos are straight forward, but issues 3&4 worry me a bit as this is auto-generated from the database and should be correct, which suggests that there's a flaw in my logic somewhere that I'll need to fix. I'm busy today and am only at work for 2 days next week so it might not be until after Christmas that I fix this.

[RolfSander]

No problem, take your time. Let's continue next year...

Have a nice Christmas time!

[stulacy]

Merry Christmas to you too!

[stulacy]

Hi Rolf, I found some time to look at this this morning and the complex rate issue was far simpler than I'd expected (I'd forgotten to take into account the website now hosts both the MCM and the CRI). Have a look at the updated exported files and let me know if they look OK now. constants_mcm.f90 mcm_export.eqn

[RolfSander]

Great, thanks!

The files are (almost) perfect now.

The only remaining issue is O2 in reactions <2> and <19>. We've discussed this earlier, but I forgot to mention it again in my recent posts. The easiest (and my preferred) solution would be to change the products in these reactions to PROD.

However, if you would like to keep the products, the line

O2 = IGNORE ;

would be necessary in the #DEFVAR section. We would then also have to check carefully what we do with the product H2O in reaction <19>.

[stulacy]

Oops I missed the PROD part of your original comment. I've made the change here so eqns 2 & 19 use PROD: mcm_export.eqn

I'm a bit confused with how the species are handled and the difference between defining a species in the mechanism file (with = IGNORE) and defining it as a variable in the Fortran constants_mcm.f90 file. Currently H2O is defined in the mechanism file with H2O = IGNORE ; and used in some reaction rates (e.g. 13), but as H2O rather than C(ind_H2O). It's also initialised in the Fortran code as a REAL and used in a complex rate definition (e.g. KMT06). O2 isn't defined in the mechanism file but is used in reaction rates (e.g. 1), and is also initialised in the constants file and used in complex rates (e.g. KMT18).

Is that correct / as expected?

[RolfSander]

Oops I missed the PROD part of your original comment. I've made the change here so eqns 2 & 19 use PROD: mcm_export.eqn

Thanks! Now it is (completely) perfect :-)

I'm a bit confused with how the species are handled and the difference between defining a species in the mechanism file (with = IGNORE) and defining it as a variable in the Fortran constants_mcm.f90 file. Currently H2O is defined in the mechanism file with H2O = IGNORE ; and used in some reaction rates (e.g. 13), but as H2O rather than C(ind_H2O). It's also initialised in the Fortran code as a REAL and used in a complex rate definition (e.g. KMT06). O2 isn't defined in the mechanism file but is used in reaction rates (e.g. 1), and is also initialised in the constants file and used in complex rates (e.g. KMT18). Is that correct / as expected?

Yes, this can be confusing. It's important to note that KPP species and Fortran variables are completely different things. For most KPP species, there is no Fortran variable with the same name. To access the chemical XYZ in Fortran, you have to use the corresponding element C(ind_XYZ) from the C array.

Only for N2, O2 and H2O, the MCM uses Fortran variables with the same name. They of course refer to the same chemical. However, from a programming point of view, they are not related at all. If a user wants to use the water from KPP as input for humidity-dependent rate constants, they have to define in their Fortran code:

H2O = C(ind_H2O)

Alternatively, they could adopt the humidity from the meteorological part of their model, e.g.:

H2O = myfunc(humidity)

Part of the confusion may result from the fact that a reaction in the *.eqn file contains two different syntaxes: It starts with KPP syntax until the colon, and then it switches to Fortran syntax until the semicolon. For example, O3 in reaction <1> is a KPP species, whereas O2 in the same line is a Fortran variable.

HTH :-)

[stulacy]

Brilliant, I've made those changes live. Feel free to reopen this ticket if you spot anything else!

Also thank you for that explainer, that's helped to solidify my mental model of how the 2 files interact.