include reactions from the MCM

RolfSander commented 2 years ago

Create a tool that facilitates the inclusion of a set of reactions from the MCM.

yantosca commented 2 years ago

Was looking at the MCM. There are some reactions that need hand editing:

{2.}     O + O3 = :     8.0D-12*EXP(-2060/TEMP)     ;    # should be O + O3 = 2O2
{19.}    OH + HO2 = :   4.8D-11*EXP(250/TEMP)   ;   # should be OH + HO2 = H2O + O2

yantosca commented 2 years ago

Also there is this section where the MCM assumes that you are in-lining these species:

#INLINE F90_GLOBAL 
 REAL(dp)::M, N2, O2, RO2, H2O 
#ENDINLINE {above lines go into MODULE KPP_ROOT_Global}

but wondering if we should add these to the #DEFVAR section (esp. O2 and H2O), because KPP won't parse it otherwise.

RolfSander commented 2 years ago

Yes, these are two of the problems with the MCM-generated files. There are more:

I don't think that "USE constants" and "CALL mcm_constants" are needed. Actually, I currently can't even find where I can download constants.f90 from the MCM web page.
negative exponents should be in parentheses, e.g.

KC0 = 3.28D-28M(TEMP/300)**-6.87

should be:

KC0 = 3.28D-28M(TEMP/300)**(-6.87)
In the definition of KMT06, the term "H2O" must be replaced by "C(ind_H2O)".
Photolysis reactions don't include "hv" as a reagent.
When upgrading to KPP-3.0.0, the line #INCLUDE atoms must be changed to #INCLUDE atoms.kpp.

I currently have (MECCA-specific) awk and tcsh scripts that take care of all these problems. I would like to replace these files with a new, MECCA-independent, python script. Depending on our timing, it may or may not be ready when we want to submit our manuscript.

yantosca commented 2 years ago

Also adding Killian Murphy (@kilicomu) to this thread. He is the maintainer of the MCM at U. York, UK. He and I had a long talk about KPP at IGC10. The MCM group is updating the mechanism and we can feed them any updates about the .kpp file.

RolfSander commented 2 years ago

Hello @kilicomu!

When I downloaded reactions in KPP format from the MCM a few days ago, I noticed that one of the problems I used to have is now fixed: The MCM-generated files used to be a mixture of DOS and UNIX format. Thanks for fixing this!

There are a couple of other problems (listed above in my previous post), and I currently have to run a shell script to fix them in the downloaded MCM files.

It would be great if we can have a look at this together and check how to fix it directly on the MCM web page. The MCM is a wonderful collection of reactions, and you can make a lot of atmospheric chemistry modelers happy if the MCM-generated KPP files can be used out-of-the-box :-)

RolfSander commented 2 years ago

I'm not sure if @kilicomu is reading this discussion. Are the any other ways to contact him?

yantosca commented 2 years ago

I've pinged him on my Slack.

kilicomu commented 2 years ago

@RolfSander I am indeed reading this discussion! I'm away from the office until 2022-07-18 minimum but will be getting to this soon after I'm back.

RolfSander commented 2 years ago

@kilicomu: Thanks for letting us know. I look forward to our collaboration when you're back!

RolfSander commented 2 years ago

@kilicomu: I'm back from vacation and ready to start working on the compatibility of KPP and the MCM!

kilicomu commented 1 year ago

@RolfSander I haven't forgotten about this! I've got a lot of catch-up work on plate at the moment due to being away from the office for an extended period of time. I will get back round to this once I've finished up a few other things.

RolfSander commented 1 year ago

@kilicomu: Sounds good, thanks! Let me know when you have time to work on this project!

BTW: I'm currently working on the analysis and visualization of chemical mechanisms. Once we have the MCM in a fully compatible KPP format, it will be easy to plot all reactions of a selected MCM species. Or all reactions going from A to B. Or all MCM species containing sulfur. Or...

RolfSander commented 1 year ago

Hello @kilicomu,

Last year, we discussed a couple of changes that are necessary in the MCM in order to generate valid KPP files. If you can find the time, it would be great if we can work on this together.

Something else: I have now submitted a manuscript about MEXPLORER, my analysis and visualization software for chemical mechanisms. It is available here:

https://doi.org/10.5194/egusphere-2023-1577

The manuscript contains some examples taken from the MCM mechanism. Figure 3 shows all pathways from isoprene to methacrolein in the MCM, automatically color-coded based on number of carbon atoms. Figure 5 shows all sources and sinks of VINOH in the MCM.

kilicomu commented 1 year ago

@RolfSander Hi - it's been a while. Firstly, MEXPLORER looks great, congratulations on completing the project.

We are currently under quite a lot of time pressure to get a minimally viable replacement for the existing MCM application completed before the old applications get shut down. At the moment, FACSIMILE export is our primary export target, but KPP export is on our list of things to reimplement when we can make time (we are doing this work unfunded).

So, I propose that when this new application goes live, we record a GitHub issue in the new application repository that describes the changes you would like to see, and we go from there. You would also be more than welcome to contribute to the application to get the KPP export back online and updated!

We are looking at early - mid September for the new application to go online, just to give you an idea of timescale.

RolfSander commented 1 year ago

Hello @kilicomu,

I look forward to see your new MCM application, and I'll be happy to help you with the KPP export. What computer language and what database are you using for the new system?

At the moment, FACSIMILE export is our primary export target, but KPP export is on our list of things to reimplement when we can make time

Indeed, FACSIMILE has always been the main target language of the MCM. However, I think that the usage of FACSIMILE by the scientific community has decreased a lot for a couple of reasons, e.g:

It is commercial, and a license must be bought.
It is not easy to integrate it into the code of large atmospheric fortran models.

KPP is free, open source, and generates fortran code that can directly be included into the code base of larger models. KPP is used by many models now:

https://en.wikipedia.org/wiki/Kinetic_PreProcessor

kilicomu commented 1 year ago

Glad to hear it. The new application has been written in Ruby, with a SQLite database behind the scenes.

I'm well aware of KPP's usefulness and increasing uptake! We just focused on FACSIMILE as it has been the primary export format from the MCM application for some time.

We intend to provide first class support for KPP mechanisms as soon as is feasible - we don't expect an initial implementation of a KPP export to take too long, we're just thin on time and have to focus on replicating the core features of the current application.

RolfSander commented 1 year ago

Thanks for the info! I have no experience with Ruby yet but this can of course change :-)

Regarding the database, I am familiar with postgreSQL. Using SQLite should be no problem.

yantosca commented 1 year ago

Thanks @RolfSander. I'll take a look at your MS when I get a chance. Very exciting stuff!

RolfSander commented 11 months ago

Hello @kilicomu, I noticed that the new MCM web page with a fresh design is now online. Nice! Unfortunately, the export in KPP format has not been activated yet :-( If there is anything I can contribute, e.g., as an alpha- or beta-tester, please let me know!

kilicomu commented 11 months ago

Hi @RolfSander - we're working on getting the KPP export going, sorry for the delay. We'll let you know when there is something ready to test.

RolfSander commented 11 months ago

OK, thanks. I look forward to seeing it...

stulacy commented 10 months ago

Hi @RolfSander and @yantosca :wave:, I'm the developer who has been working on the MCM redesign. Thank you for your patience while we launched the new website. I'm now looking at getting the KPP export functionality up and running and would welcome your input in the PR thread in the MCM repo. I've got a basic version working and am now going through the modifications mentioned in this thread. Cheers! Stuart

RolfSander commented 10 months ago

Hello @stulacy,

Great, thanks! I've posted a comment on your MCM repo page.

RolfSander commented 9 months ago

Hello @yantosca, @jimmielin, @msl3v and @laestrada,

As you may have seen already, an improved export in KPP syntax has recently been activated on the new MCM web pages:

https://mcm.york.ac.uk/MCM/export

This solves most of the problems that we have discussed here. However, there are still a few remaining questions (see also https://github.com/wacl-york/mcm-web/pull/210). As the rate constants in the MCM depend on several external parameters (M, H2O, j-values, etc.), I see two possibilities to define them:

1) We define them via #INLINE F90_GLOBAL and #INLINE F90_RCONST. Then everything will be in the MCM-generated *.eqn file.

2) We place all declarations and definitions into SUBROUTINE mcm_constants inside a separate file called mcm_constants.f90. Then we only have to put USE mcm_constants and CALL mcm_constants into #INLINE F90_RCONST.

I think I like the 2nd solution better. Especially when the external parameters come from a larger model (e.g., GEOS-Chem or MECCA), this approach looks more flexible to me. Prviding your own mcm_constants.f90 as an interface between the large model and KPP should be easier than modifying each MCM-generated *.eqn file.

I'd like to hear your opinion as well. Also, let me know in case you need more information!

yantosca commented 9 months ago

@RolfSander: I also like the 2nd solution better. For GEOS-Chem this is how we inline all of the rate-law functions into KPP:

#MINVERSION   3.0.0                  { Need this version of KPP or later     }
#INTEGRATOR   rosenbrock_autoreduce  { Use Rosenbrock integration method     }
#AUTOREDUCE   on                     { ...w/ autoreduce enabled but optional }
#LANGUAGE     Fortran90              { Generate solver code in Fortran90 ... }
#UPPERCASEF90 on                     { ...w/ .F90 suffix (instead of .f90)   }
#DRIVER       none                   { Do not create gckpp_Main.F90          }
#HESSIAN      off                    { Do not create the Hessian matrix      }
#MEX          off                    { MEX is for Matlab, so skip it         }
#STOICMAT     off                    { Do not create stoichiometric matrix   }

#INCLUDE fullchem.eqn     { Chemical reactions for fullchem mechanism  }

#FAMILIES                 { Chemical families for prod/loss diagnostic }
POx : O3 + NO2 + 2NO3 + PAN + PPN + MPAN + HNO4 + 3N2O5 + HNO3 + BrO + HOBr + BrNO2 + 2BrNO3 + MPN + ETHLN + MVKN + MCRHN + MCRHNB + PROPNN + R4N2 + PRN1 + PRPN + R4N1 + HONIT + MONITS + MONITU + OLND + OLNN + IHN1 + IHN2 + IHN3 + IHN4 + INPB + INPD + ICN + 2IDN + ITCN + ITHN + ISOPNOO1 + ISOPNOO2 + INO2B + INO2D + INA + IDHNBOO + IDHNDOO1 + IDHNDOO2 + IHPNBOO + IHPNDOO + ICNOO + 2IDNOO + MACRNO2 + ClO + HOCl + ClNO2 + 2ClNO3 + 2Cl2O2 + 2OClO + O + O1D + IO + HOI + IONO + 2IONO2 + 2OIO + 2I2O2 + 3I2O3 + 4I2O4;
LOx : O3 + NO2 + 2NO3 + PAN + PPN + MPAN + HNO4 + 3N2O5 + HNO3 + BrO + HOBr + BrNO2 + 2BrNO3 + MPN + ETHLN + MVKN + MCRHN + MCRHNB + PROPNN + R4N2 + PRN1 + PRPN + R4N1 + HONIT + MONITS + MONITU + OLND + OLNN + IHN1 + IHN2 + IHN3 + IHN4 + INPB + INPD + ICN + 2IDN + ITCN + ITHN + ISOPNOO1 + ISOPNOO2 + INO2B + INO2D + INA + IDHNBOO + IDHNDOO1 + IDHNDOO2 + IHPNBOO + IHPNDOO + ICNOO + 2IDNOO + MACRNO2 + ClO + HOCl + ClNO2 + 2ClNO3 + 2Cl2O2 + 2OClO + O + O1D + IO + HOI + IONO + 2IONO2 + 2OIO + 2I2O2 + 3I2O3 + 4I2O4;
PCO : CO;
LCO : CO;
PSO4 : SO4;
LCH4 : CH4;
PH2O2 : H2O2;

#INLINE F90_RATES
  ! All rates are included in fullchem_RateLawFuncs.F90, which
  ! gets referenced directly from subroutine Update_Rconst.
#ENDINLINE

#INLINE F90_RCONST
  ! Inline an include file containing rate law definitions, which
  ! will be inserted directly into subroutine Update_Rconst().
  ! This is necessary as a workaround for KPP not being able to
  ! include very large files ( > 200000 chars) directly.
  !  -- Bob Yantosca (11 Jun 2021)
  USE fullchem_RateLawFuncs
#ENDINLINE

#INLINE F90_GLOBAL
#include "commonIncludeVars.H"
#ENDINLINE

For GEOS-Chem, we also use an include file via #INLINE F90_GLOBAL... the include file contains common variables for all 3 of our mechanisms (fullchem, Hg, carbon). Part of the reason we use both methods is that at compile time you have to have all variables present to satisfy dependencies., and this seemed like it was an easy way to do it. But for the MCM I think you can put everything into the mcm_constants.f90.

yantosca commented 9 months ago

@RolfSander @jimmielin @mslv: Also note, I had to add comments in the F90_RATES section to avoid a KPP parsing error.

#INLINE F90_RATES
  ! All rates are included in fullchem_RateLawFuncs.F90, which
  ! gets referenced directly from subroutine Update_Rconst.
#ENDINLINE

RolfSander commented 9 months ago

Thanks for the explanations! If there are no objections, I will create a small MWE (minimal working example) based on the MCM isoprene mechanism, and applying the 2nd method.

RolfSander commented 9 months ago

I have now created a small minimal working example based on the MCM isoprene mechanism. The code is in the examples/mcm/ directory which you can find in the new "mcm" branch:

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

It should be easy to run the code as explained in the README file.

A few notes on the code:

All external input for KPP has been moved out of the *.eqn file and into the new file constants_mcm.f90.
The new file is called constants_mcm.f90 (not mcmconstants.f90 anymore) because with the old name it didn't survive gmake distclean which deletes all mcm*.f90 files.
The KPP-generated Makefile_mcm doesn't work because it doesn't consider the new file. There is a separate file (simply called Makefile) which must be used instead.
Instead of using index numbers to select j-values from the J() array, I have defined several INTEGER parameters. For example, to get the photolysis rate constant for H2O2 + hv, you can now write J(J_H2O2). The original numbers (https://mcm.york.ac.uk/MCM/rates/photolysis) were problematic because there are gaps, and also because they could change when selecting a different mechanism.
A mechanism may or may not contain a chemical species called H2O. For humidity-dependent rate constants, the variable H2O is now used, not C(ind_H2O) anymore. It is up to the user to define this variable in a suitable way. In the MWE, I'm using a constant value of 0.01.
If all works well, we may add a MCM-generated mechanism to our CI tests in the future.

Please let me know what you think! If you agree, I will discuss with Stuart (the MCM developer), if the MCM can generate the format suggested by us...

yantosca commented 9 months ago

Thanks @RolfSander. I'll look at this when I get a chance. Happy to add the C-I test for the MCM example.

RolfSander commented 9 months ago

Hello @yantosca,

Thanks, I look forward to reading your comments!

Also, thanks for your offer to add an MCM CI test! As we have an extra file here, we probably cannot use function run_ci_test(). Instead, I think we'll need a new function, e.g., run_mcm_ci_test.

yantosca commented 9 months ago

@RolfSander: I just pushed a PR #83 that develops a C-I test for the MCM plus also adds some Makefile updates. Let me know what you think.

RolfSander commented 9 months ago

Thanks, @yantosca, for adding the CI test!

If there are no objections, I will contact the MCM developer soon to ask if mcm_isoprene.eqn and constants_mcm.f90 can be generated automatically on the MCM web page so that they can be used by KPP directly (without the need of any postprocessing).

yantosca commented 9 months ago

Thanks @RolfSander. Let me know when the MCM stuff is complete, we can then release it as KPP 3.1.0.

RolfSander commented 9 months ago

Yes, I will mention it here when everything is ready.

If anyone is interested in the work in progress, you can find the details in this issue:

https://github.com/wacl-york/mcm-web/issues/177

RolfSander commented 8 months ago

Thanks to @stulacy, the MCM can now export a mechanism in KPP format. It works out-of-the-box with the model that we have in examples/mcm. I've started a pull request, see https://github.com/KineticPreProcessor/KPP/pull/92.

KineticPreProcessor / KPP

include reactions from the MCM #4