[PULL REQUEST] Clean up heterogenous chemistry rate-law functions -- excluding sulfur chemistry updates to KPP #897
Conversation
Resolved conflicts in: - GeosCore/isorropiaII_mod.F90 Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
Configuring with cmake ... -DCUSTOMMECH=y now points to the chemical mechanism files in KPP/aerochem (instead of KPP/custom). This is for development of the KPP mechanism, and will be restored to the original functionality later. Also rebuilt the KPP/aerochem source code with KPP 2.3.1_gc. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
In the heterogeneous chemistry, the square root of molecular weight is repeatedly computed on every (I,J,L). However, this can be done outside of the loop over grid boxes. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
This is the start of the process to clean up gckpp_HetRates.F90. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
(1) Now use KIND parameter dp instead of fp. dp is in gckpp_Precision.F90.
(2) The heterogeneous chemistry rates for LVOC, IHN{1-4}, INP{B,D},
MCRHN, NCRHNB, MVKN, R4N2, IDN, ITHN, ITCN, MONITS, MONITU, and
HONIT can now be computed by a single function, HetVOC.
(3) The hetchem rates for IONITA and MONITA are constant, and can
just be set equal to a numeric constant w/o a function call.
(4) Cosmetic + whitespace changes
Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
In KPP/fullchem/gckpp_HetRates.F90, we have split function HetIxCycleSSA into 4 functions HetIxCycleBrSALAinstead of using a single function HetIxCycleSSA, we have split that function into 4 individual functions: 1. HetIxCycleBrSALA (formerly HetIxCycleSSA with N=1) 2. HetIxCycleBrSALC (formerly HetIxCycleSSA with N=2) 3. HetIxCycleSALACL (formerly HetIxCycleSSA with N=3) 4. HetIxCycleSALCCL (formerly HetIxCycleSSA with N=4) This allows us to remove branching and looping within the functions, which should allow for more efficient execution. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
We are able to remove the rate-law functions for computing uptake
of iodine species on SO4, SALA, SALC, BrSALA, BrSALC, SALACL, SALCCL
into 3 functions:
(1) HetIUptakebySulf : Uptake by tropospheric sulfate
(2) HetIUptakebySALA : Uptake by accumulation mode sea-salt
(3) HetIUptakebySALC : Uptake by coarse-mode sea-salt
IBr produced on acidic accum-mode seasalt (when SSAlk(1) <= 0.5):
Uptake by BrSALA = 0.15 * HetIUptakeBySALA
Uptake by SALACL = 0.85 * HetIUptakeBySALA
ICl produced on acidic accum-mode seasalt (when SSAlk(2) <= 0.5):
Uptake by BrSALC = 0.15 * HetIUptakeBySALC
Uptake by SALCCL = 0.85 * HetIUptakeBySALC
Thus we can retire the HetIxCycle{BrSALA,BrSALC,SALACL,SALCCL}.
Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
In KPP/fullchem/gckpp_HetRates.F90:
(1) Updates to routine CLD_PARAMS:
- General cleanup
- Return argument ClearFr = 1 - CldFr
- Now reference most met variables from State_Met
- Updated comments & whitespace updates
(2) Updates to routine GAMMA_NO3:
- Whitespace alignment and readability updates
- Now use H%NO3%SrMw instead of taking the square root
- Return value was renamed to "gamma"
(3) Updates to routine HetNO3_Cl
- Whitespace alignment and readability updates
- Now pass H%NO3%SrMw to ARSL1K instead of taking the
square root of the molecular weight of NO3
Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
We now have documented the reactions that are stored in each element of the HET array in gckpp.kpp. We have also rebuilt the KPP files with KPP-for-GEOS-Chem 2.3.1_gc, so that these comments are now also included in gckpp_Global.F90. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
The HETOH routine calls routine GAMMA_NO3 and ARSL1K twice. These have now been inlined into the main body of routine SET_HET and the HETOH routine has been removed. NOTE: This update will cause differences w/r/t the prior code, because a hardwired molecular weight for OH (17 g/mole) was used there... but the value from the species database is now 17.01 g/mole. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
…hanges gckpp.kpp and gckpp_Global.F90: - Add comments to HSTAR_EPOX to denote its value is from Ngyuen et al 2014 gckpp_HetRates.F90: - Add comments describing each aerosol type - Remove looping & branching in 1st-order uptake functions - Now call subroutine ArsL1k for each aerosol type in 1st-order uptake fns - ArsL1k now returns 0.0 instead of 1e-30 if rate cannot be computed print_het.H - Also print out temperature and relative humidity at the debug box Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
(1) For now, ARSL1K returns 1.0e-30_dp, to avoid div-by-zero issues
until we can better trap these errors.
(2) Moved Het1stOrderUptakeHO2 so that it is in alphabetical order
(3) Reorderd calls to VOC 1st order subroutine to be in alphabetical order.
Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
Move the nullifcation of the H (HetInfo) object after the print_het.H print statements. Also cosmetic changes, remove REVISION HISTORY sections from comment headers. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
This is a cosmetic change. All 1st order uptake reactions now match the reference code 13.0.1 (excluding numerical noise diffs smaller than +/- 2.0e-5%). Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
(1) Use KPP version 2.3.2_gc (in beta) to rebuild the aerochem
mechanism files.
(2) Remove XDENA, as this is a duplicate of NUMDEN.
(3) Replace 300.0_dp/TEMP with K300_OVER_TEMP variable.
(4) Comment out the UPDATE_SUN function in gckpp_Rates.F90.
(5) Keep gckpp_Rates.F90 in aerochem separate from fullchem
(i.e. these are no longer symlinked).
Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
…2_gc GeosCore/flexchem_mod.F90: - Initialize TEMP_OVER_K300, K300_OVER_TEMP and SR_TEMP in routine Set_Kpp_Gridbox_Values - Initialize MW and SR_MW in Init_FlexChem KPP/aerochem/gckpp_HetRates.F90: - Pass SR_MW to the 1st-order uptake functions - Remove TEMP and NUMDEN from CloudHet, these are global variables gckpp_Rates.F90 - Rebuilt with KPP 2.3.2_gc Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
The computation of branch assumes that kI is always greater than zero. But we have changed ARSL1K to return 0 if the reaction probability is undefined. This would cause a div-by-zero in Cloudhet. We have now prevented this by only allowing the division if kI > 0, otherwise we set branch = 0. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
gckpp.kpp - Added T_298 = 298.15_dp parameter, as this is reused - Added HENRY_K0 and HENRY_CR parameters gckpp_HetRates.F90: - Define HENRY_K0 and HENRY_CR with values from the species database - Remove TEMP, NUMDEN from HetNO3_CL (these are global) - Remove T, H from GAMMA_NO3 (these are not needed) gckpp_Global.F90 - Rebuilt with KPP 2.3.2_gc
aerochem.eqn
- Het rxn HO2 = O2 is now changed to HO2 = H2O
gckpp.kpp:
- Now use INV_T298 as the inverse of T_298, to prevent repeated divisions
gckpp_HetRates.F90:
- Streamlined HetNO3_Cl, HetBrNO3, HetClNO3_TCld routines
gckpp_{Function, Global, Jacobian, Monitor}.F90
- Rebuilt with KPP 2.3.2_gc
Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
In function CloudHet (in gckpp_HetRates.F90): - Remove DO loop; now compute liquid and ice branches w/ inlined code - Arguments branchLiq and branchIce are no longer OPTIONAL, this should result in more efficient execution Also make sure to pass branchLiq and branchIce with placeholder arguments (1.0_dp) in calls to CloudHet where necessary. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
(1) Pass branchliq, branchice to CloudHet as non-optional arguments
(2) Remove PRESENT statements from CloudHet
(3) Add parameters to Gamma_Cld_{Aer, Ice} to avoid repeated computations
Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
We have been able to peel off many of the gas-phase 1st-order rate-law fucnctions from gckpp_HetRates.F90 and add them into gckpp.kpp. From there, they will be included the KPP-generated source code whenever the mechanism is rebuilt. flexchem_mod.f90: - Use keyword arguments in calls to Set_Het and Set_Kpp_Gridbox_Values state_chm_mod.F90: - Removed the HetInfo object and HetState type gckpp.kpp - Add new rate law functions from gckpp_HetRates.F90 - Define derived type object State_Het to pass fields to rate-law functions - Define new HetState type for the State_Het object aerochem.eqn - Now reference rate-law functions defined in gckpp.kpp gckpp_HetRates.F90 - Moved several 1st-order rate-law functions to gckpp.kpp - Moved Cld_Params and the computation of halide concs to flexchem_mod.F90 gckpp_Global.F90 gckpp_Rates.F90 - Rebuilt with KPP 2.3.2_gc Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
We have moved the rate-law functions that compute 1st-order uptake of HI, I2O2, I2O3, and I2O4 from gckpp_HetRates.F90 to gckpp.kpp. We have also made the corresponding updates in aerochem.eqn, and have regenerated the solver source code files with KPP 2.3.2_gc. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
gckpp.kpp: - Add logicals ssFineisAlk, ssFineisAcid, ssCoarseisAlk, ssCoarseIsAcid to the State_Het object. These are used to determine if the fine/coarse sea salt is acidic or alkaline in rate-law functions. This allows us to compute these terms once per grid box. - Update rate-law functions to use common variables from gckpp_Global.F90 flexchem_mod.F90: - Define the ssFineIsAlk, ssFineIsAcid, ssCoarseIsAlk, ssCoarseIsAcid fields of State_Het for each grid box. aerochem.eqn - Now use rate-law functions for the first-order uptake of iodine from SALA and SALC. gckpp_Global.F90 gckpp_Rates.F90 - Rebuilt with KPP 2.3.2_gc gckpp_HetRates.F90 - Remove calls to iodine uptake rate-law functions Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
gckpp.kpp: - Add function kIIR1Ltd (and optimize it to minimize IF/ELSE branches) - Remove comments about elemetns of HET that are deleted - Trim ending whitespace aerochem.eqn - Call the various iodine uptake rate-law functions gckpp_Global.F90 gckpp_Rates.F90 - Rebuilt with KPP 2.3.2_gc Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
KPP/fullchem/rateLawUtilFuncs: - The line "IF ( concEduct < 100.0_dp ) RETURN" was evidently intended to prevent division by zero. However, this will flush rate law computations to zero for very small concentrations. - Given that we now use the Is_SafeDiv and SafeDiv functions to prevent div-by-zero errors, this error trap is no longer necessary. Comment it out for now. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
We will add this back once the HMS reactions are activated in KPP. For now, remove this in order to avoid a missing species error. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
Moved the computation of M_X and c_avg into the IF ( k_tot > 0.0_dp ) block. This will skip computation of these quantities if the total rate k_tot is zero -- then we just exit. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
Headers/state_chm_mod.F90 - Set Do_SulfateMod_Cld = T always - Set Do_SulfateMod_Seasalt = T if it is a fullchem simulation GeosCore/fullchem_mod.F90: - Add subroutine PrintFirstTimeInfo to print an informational message about which options are being used KPP/fullchem/fullchem.eqn: KPP/aciduptake/aciduptake.eqn - Disable reactions using K_CLD(1:6) - Disable reactions for HMS KPP/fullchem/gckpp_*.F90 - Rebuilt with KPP 2.3_3_gc KPP/aciduptake/gckpp_*.F90 - Removed for present. This will be activated at a later date Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
We have modified the createRunDir.sh script so that restart file GEOSCHEM_RESTARTS/v2021-09/GEOSChem.restart.fullchem.20190701_0000z.nc4 will be copied to full-chemistry run directories. This was necessary because we needed to include new species HSO3m and SO3mm. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
…e_mod It is not possible to use the hetchem reactions with HSO3mm and SO3mm while computing sulfur cloud and seasalt chemistry in sulfate_mod.F90. Therefore, we have restored the SO4H1, SO4H2, SO4H3, and SO4H4 species in fullchem.eqn and updated the rate law functions for HOCl and HOBr accordingly. GeosCore/fullchem_mod.F90 - Comment out references to ind_HSO3m and ind_SO3mm, as HSO3m and SO3mm species are not included in KPP. Headers/state_chm_mod.F90 - Set Do_SulfateMod_Seasalt = TRUE and Do_SulfateMod_Cld = TRUE. This will ensure that sulfur seasalt & cloud chemistry will still be done by sulfate_mod.F90, while we evaluate the KPP sulfur chemistry code. KPP/fullchem/fullchem.kpp KPP/aciduptake/aciduptake.kpp - Add fupdateHOBr and fupdateHOCl to TYPE(HetState) KPP/fullchem/fullchem.eqn KPP/aciduptake/aciduptake.eqn - Restore SO4H1, SO4H2, SO4H3, SO4H4 species and reactions KPP/fullchem/fullchem_HetStateFuncs.F90 - Copy fupdateHOBr and fupdateHOCL from State_Chm to State_het KPP/fullchem/fullchem_RateLawFuncs.F90 - Comment out functions HOBrUptkByHSO3m and HOBrUptkBySO3mm - Comment out functions HOClUptkByHSO3m and HOClUptkBySO3mm - Add functions HOBrHetSO4H1 and HOBrHetSO4H2 - Add functions HOClHetSO4H4 and HOBrHetSO4H3 - Remove AqRate_* functions, these are now done in fullchem_SulfurChemFuncs.F90, stored in K_MT(1:6). KPP/fullchem/gckpp.map KPP/fullchem/gckpp_*.F90 - Rebuilt with KPP version 2.3.3_gc Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
We will restore this once the sulfur chemistry is able to be computed within KPP. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
In KPP/fullchem/fullchem_SulfurChemFuncs.F90, there is an unnecessary reference to apm_mod.F90. We have removed this. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
yantosca
added
category: Feature Request
yantosca
changed the title
|
GEOS-Chem Classic integration tests all passed: |
There was a problem hiding this comment.
I am reviewing these changes now. I see that flexchem_mod.F90 is renamed to fullchem_mod.F90. I like this change, but I'm wondering if this means that for any new mechanism that we introduce we will also need to add a new module or routine to interface between GEOS-Chem and KPP?
In this commit, we are not adding new species HSO3- and SO3--, so we do not need to copy the restart file from the v2021-09 folder of the restart file path. Therefore we have restored the original restart file copy command in run/GCClassic/createRunDir.sh. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
The metals simulations were running for an entire month when executing a full integration test. A fix has been added to commonFunctionsForTests.sh to automatically modify the end date of metals simulations so they only run for 20 minutes. Signed-off-by: Melissa Sulprizio <mpayer@seas.harvard.edu>
|
GCHP integration test results: |
That was our idea. For full-chemistry mechanisms we can link to fullchem mod (as eventually the aciduptake will be a new mechanism with dust aciduptake rxns handled by KPP). For the Hg code, the interface between GEOS-Chem and the solver is done in mercury_mod.F90. If we make a CO-CO2-CH4-OCS mechanism, then we'll need to write a new interface module because there wouldn't be the need to burden |
There was a problem hiding this comment.
Hi @yantosca. I finished reviewing the files. All changes look OK to me, but I had some comments and questions.
-
Is there a reason State_Chm%Do_SulfateMod_Seasalt and State_Chm%Do_SulfateMod_Cld are stored in State_Chm and not Input_Opt with other switches? It seems strange to have logical switches in State_Chm and set them there, but maybe it's fine since these are temporary?
-
What is the
Is_Forgefield added in species_mod.F90? It doesn't seem to be used. -
I notice a new KPP mechanism for
aciduptake. Theaciduptake.eqnfile appears to be identical tofullchem.eqnwith the exception of additional reactions needed for acid uptake. We have been moving away from multiple full-chemistry mechanisms. Are there plans to eventually only include the relevant additional acid uptake reactions/species in this mechanism and then combine with the fullchem mechanism? The idea would be to avoid having to update reactions in multiple locations if a change is made anywhere infullchem.eqn. Maybe it's fine to leave this unaddressed for now if it's not activated here, but I just wanted to jot down my thoughts. -
There are some hardcoded paths (
/home/ilcentro/Work/Harvard/GC/v13/KPP/fullchem/...) in many files inKPP/fullchem/CMakeFiles/KPPFirstPass_fullchem.dir/. These should be removed or made more flexible. If this is output produced by CMake, we should consider adding them to .gitignore. -
There are also hardcoded paths to
/home/ilcentro/Work/Harvard/GC/v13inKPP/fullchem/Makefile. -
KPP/fullchem/aciduptake_DustChemFuncs.F90appears to be a stub. Is this a work in progress or still needed here? -
In fullchem.eqn reaction HO2 = O2 is changed to HO2 = H2O. Is this a bug fix or science update? Are impacts from this change quantified?
|
@yantosca Could you also please link to the 1-month benchmark results here for record keeping and to confirm that changes from this update are mainly due to numerical noise? |
These are temporary switches. They could go into Input_Opt if you wish. They won't be around forever, but right now we need them to tell the code to use sulfate_mod.F90 instead of KPP for certain reactions.
I have no idea. Maybe Mike put that in. I can take it out.
So eventually aciduptake will be a new chemical mechanism with dust uptake reactions handled by KPP. This was done in order to avoid burdening the user with having all of the 12 extra dust uptake species in the standard "fullchem" mechanism. The KPP/aciduptake folder is left for further development. I can put a README.md file there to denote that. Right now the aciduptake chemistry is handled by sulfate_mod.F90.
Yeah, that should probably not be there. It might be a CMake-generated file. I'll take a look.
I'll delete that.
That is a stub module that is needed to get the code to compile, since both the fullchem and aciduptake KPP code will share fullchem_mod.F90.
Daniel told me to do that. That had been a long-standing typo. I believe the impacts are not large. |
I'm fine with leaving as-is since they are temporary. I suppose they should not be easily toggled since users aren't meant to use those options right now.
Thanks!
OK. Sounds good. Perhaps we can discuss and clean it up / consolidate later, but getting it working for now is fine.
Thanks!
Thanks!
Sounds good.
OK. Thanks for confirming. @yantosca After you make the minor changes indicated above, please feel free to merge into |
These are not needed. Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
Headers/species_mod.F90 - Removed IS_FORGE tag in Species derived type KPP/build_mechanism.sh - Remove Makefile_gckpp after each time KPP rebuilds the gckpp_* files. Makefile_gckpp is for GNU Make but we now use CMake. KPP/aciduptake/Makefile_gckpp - Removed Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>
This PR will supersede #827. This contains modifications for heterogenous chemistry reactions and rate law functions WITHOUT the sulfur chemistry modifications of Mike Long (as these need further testing & validation).
Other chemistry updates will be added after this PR has been merged into dev.
Integration tests are pending. I will report the results here.