include reactions from the MCM

[RolfSander]

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

[yantosca]

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]

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]

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]

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]

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]

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

[yantosca]

I've pinged him on my Slack.

[kilicomu]

@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]

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

[RolfSander]

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

[kilicomu]

@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]

@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]

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]

@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]

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]

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]

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]

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

[RolfSander]

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]

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]

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

[stulacy]

Hi @RolfSander and @yantosca 👋, 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]

Hello @stulacy,

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

[RolfSander]

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
wacl-york/mcm-web#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]

@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]

@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]

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]

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 mcm_constants.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]

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

[RolfSander]

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]

@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]

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]

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

[RolfSander]

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:

wacl-york/mcm-web#177

[RolfSander]

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 #92.