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>

gckpp.kpp
- Added Gam_NO3, NO3hypsisClonSALA, NO3hypsisClonSALC functions
- Added coth and reactodiff_corr functions
- Removed useless computations from Gam_NO3
- Trimmed trailing whitespace

aerochem.eqn
- Replace HET(ind_NO3,2) with function NO3hypsisClonSALA
- Replace Het(ind_NO3,3) with function NO3hypsisClonSALC

gckpp_HetRates.F90
- Removed GAMMA_NO3, HetNO3_Cl_SALA, HetNO3_Cl_SALC functions
- Removed COTH, REACTODIFF_CORR functions

gckpp_Rates.F90
- Rebuilt with KPP 2.3.2_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

gckpp.kpp
- Added CloudHet and SafeDiv functions
- Remove srMw argument from HO2uptk1stOrd function
- Add NO2uptk1stOrd and NO2uptk1stOrdAndCloud functions
- Add NO3uptk1stOrd and NO3uptk1stOrdAndCloud functions

aerochem.eqn
- Remove SR_MW(ind_HO2) from HO2uptk1stOrd function
- Replace HET(ind_NO2,1) with NO2uptk1stOrdAndCloud function
- Replace HET(ind_NO3,1) with NO3uptk1stOrdAndCloud function

gckpp_HetRates.F90
- Remove CloudHet, Het1stOrderUptakeNO3, NO3_Uptake functions
- Remove Het1stOrderUptakeNO2, NO2_Uptake functions

gckpp_Rates.F90
- Rebuilt with KPP 2.3.2_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

…eqn)

gckpp.kpp
- Minor updates for numerical stability in function CloudHet
- Add function Is_SafeDiv
- Consolidaed NO2 and NO3 cloud uptake functions
- Rename HetMinLife constant to HET_MIN_LIFE
- Rename HetMinRate constant to HET_MIN_RATE

aerochem.eqn
- Replace HET(ind_N2O5,3) with call to N2O5uptkByCloud

gckpp_Global.F90
gckpp_Rates.F90
- Rebuilt with KPP 2.3.2_gc

gckpp_HetRates.F90
- Removed code that is now in gckpp.kpp

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

gckpp.kpp
- Modify routine kIIR1Ltd for computational efficiency
- Add routine N2O5uptkByStratHCl

aerochem.eqn
- Replace HET(ind_N2O5,2) by function call to N2O5uptkByStratHCl

gckpp_Rates.F90
- Rebuilt with KPP 2.3.2_gc

gckpp_HetRates.F90
- Remove function HetN2O5_HCl

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

We have now migrated the rate-law functions that compute the update
of N2O5 by H2O from gckpp_HetRates.F90 into KPP.

gckpp.kpp
- Add functions N2O5uptkByH2O, N2O5_InorgOrg, ClNO2_BT
- Clean up N2O5uptkByH2O from original HETN2O5 function

aerochem.eqn
- Replace HET(ind_N2O5,1) with N2O5uptkByH2O function call

gckpp_HetRates.F90:
- Remove HETN2O5, N2O5, N2O5_InorgOrg, ClNO2_BT

gckpp_Global.F90
gckpp_Rates.F90
- Rebuilt with KPP 2.3.2_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

gckpp.kpp
- Added functions N2O5uptkBySALACL and N2O5uptkBySALCCL
- Updated comments

aerochem.eqn
- Replace HET(ind_N2O5,4) with call to N2O5uptkBySALACL
- Replace Het(ind_N2O5,5) with call to N2O5uptkbySALCCL

gckpp_HetRates.F90
- Remove obsolete code

gckpp_Rates.F90
- Rebuilt with KPP 2.3.2_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

gckpp.kpp
- Add rate-law functions OHlossBySALACl and OHlossBySALCCl
- Also add functions for O3 loss (but these are still being validated)
- Add BrConc_Cld and ClConc_Cld to State_Het

aerochem.eqn
- Replace HET(ind_OH,1) with call to OHuptkBySALACl
- Replace HET(ind_OH,2) with call to OHuptkBySALCCl
- Zero out HET(ind_O3,{1,2,3}) in preparation for next update

gckpp_HetRates.F90
- Removed code that has now been peeled into gckpp.kpp

gckpp_Global.F90
gckpp_Rates.F90
- Rebuilt with KPP 2.3.2_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

gckpp.kpp
- Add function O3uptkByBrInTropCloud
- Add function O3uptkByHBrInCloud
- Add function O3uptkByCloudAndBrSALA
- Add function O3uptkByCloudAndBrSALC
- Add function Gamma_O3_Br
- CloudHet now skips computing rates if the branch ratio is zero
- Rename several fields in State_Het (e.g. Brconc_* --> Br_conc_*, etc)
- Move parameter HSTAR_EPOX into EpoxUptkGamma, it's only needed there

aerochem.eqn
- Replace HET(ind_O3,1) with call to O3uptkByBrInTropCloud
- Replace HET(ind_O3,2) with call to O3uptkByCloudAndBrSALA
- Replace HET(ind_O3,3) with call to O3uptkByCloudAndBrSALC

gckpp_HetRates.F90
- Use H instead of State_Het for clarity in several routines
- Remove H%HetTemp
- Compute Br and Cl branching ratios in Halide_Conc (so that we don't
  have to keep recomputing them in the functions)

flexchem_mod.F90
- Use new field names of State_Het where necessary

species_database.yml
- Add Henry_CR and Henry_K0 for O3

gckpp_Global.F90 and gckpp_Rates.F90
- Rebuilt with KPP 2.3.2_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

gckpp.kpp:
- Add functions HBrUptkBySALA and HBrUptkBySALC
- Add parameters for aerosol types (instead of hardwired numbers)
- Add parameters for fine/coarse (for SSALK, aWATER etc)
- Add parameters for KHETI_SLA indices (instead of hardwired numbers)
- Updated array references with these parameters accordingly

aerochem.eqn:
- Replaced HET(ind_HBr,1) with call to HBrUptkBySALA
- Replaced HET(ind_HBr,2) with call to HBrUtpkBySALC

gckpp_Global.F90
gckpp_Rates.F90
- Rebuilt with KPP 2.3.2_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

gckpp.kpp:
- Added functions BrNO3uptkByH2O and BrNO3uptkByHCl
- Renamed some index parameters to reduce chance of namespace collisions
- Updated comments

aerochem.eqn:
- Replaced HET(ind_BrNO3,1) with call to BrNO3uptkByH2O
- Replaced HET(ind_BrNO3,2) with call to BrNO3uptkByHCl

gckpp_HetRates.F90
- Remove obsolete code

gckpp_Global.F90
gckpp_Rates.F90
- Rebuilt with KPP 2.3.2_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

In the hetchem rate-law functions there are many expressions that begin
with 4.0 * R * TEMP, etc.  These can be pulled out and stored into
constants so that they only need to be computed once.  We have added
as KPP global variables in gckpp.kpp:

INV_TEMP        = 1 / TEMP
FOUR_R_T        = 4 * CON_R * TEMP
FOUR_RGASLATM_T = 4 * RGASLATM * TEMP
EIGHT_RSTARG_T  = 8 * RSTARG * T

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

gckpp.kpp:
- Add functions Br2_Yield, Gam_HOBr_Cld, Gam_HOBr_Ice,
  and HOBrUptkByHBrAndTropCloud
- Add State_Het%Br_over_Cl_Cld field to the HetState derived type
- Renamed HBr_th to HBr_theta
- Renamed HCl_th to HCl_theta
- Renamed HNO3_th to HNO3_theta

aerochem.eqn
- Replace HET(ind_HOBr,1) with call to HOBrUptkByHBrAndTropCloud

gckpp_HetRates.F90:
- Streamline subroutine Get_Theta_Ice
- Assign State_Het%Br_over_Cl_Cld field

gckpp_Global.F90
gckpp_Rates.F90
- Rebuilt with KPP 2.3.2_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

gckpp.kpp
- Added functions Br2_Yield, Gam_HOBr_Cld, Gam_HOBr_Ice,
  HOBrUptkByHBrAndTropCloud, HOBrUptkByHClAndTropCloud.
- Remove branching (i.e. IF (X==1) THEN etc.) from these functions

aerochem.eqn:
- Replace HET(ind_HOBr,1) with call to HOBrUptkByHBrAndTropCloud
- Replace HET(ind_HOBr,2) with call to HOBrUptkByHClAndTropCloud

gckpp_HetRates.F90
- Remove more obsolete routines that are now in gckpp.kpp

gckpp_Rates.F90
- Rebuilt with KPP 2.3.2_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

…Chem into feature/msl/kpp

This commit combines the work done by Mike Long and Bob Yantosca
in bringing the code from sulfate_mod.F90 and gckpp_HetRates.F90
into KPP.

Resolved conflicts in:
1. GeosCore/flexchem_mod.F90
2. Headers/species_database_mod.F90
3. KPP/aerochem/aerochem.eqn
4. KPP/aerochem/gckpp.kpp

Rebuilt these files with KPP 2.3.2_gc
1. KPP/aerochem/gckpp_Function.F90
2. KPP/aerochem/gckpp_Global.F90
3. KPP/aerochem/gckpp_Jacobian.F90
4. KPP/aerochem/gckpp_Monitor.F90
5. KPP/aerochem/gckpp_Rates.F90

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

KPP/build_mechanism.sh
- Make sure we do not delete gckpp_AqRates.F90 and gckpp_Sulfate.F90

KPP/gckpp_AqRates.F90
KPP/gckpp_Sulfate.F90
- Restored from prior commit

KPP/gckpp_HetRates.F90
- Make ARSL1K, KIIR1LTD, and KIIR1R2L public
- Comment out the remaining code not already pulled into gckpp.kpp

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

Headers/species_mod.F90:
- Add nOmitted and nRealSpc to the SpcIndCt type
- Add OmittedId in the Species type
- Skip defining a new species in Spc_Create if Is_Omitted=T

Headers/species_database_mod.F90:
- Zero SpcCount%nOmitted and SpcCount%nRealSpc
- Test if a species is omitted.  If it is, increment SpcCount%nOmitted
  and then cycle to the next species
- Subtract SpcCount%nOmitted from ThisSpc%modelId
- Don't throw an error for omitted species

Headers/state_chm_mod.F90:
- Add State_Chm%nOmitted field
- Now define State_Chm%nSpecies with SpcCount%nRealSpc
- Allocate enough slots in the KPP mapping vectors so that the indexing
  (which includes the omitted species) comes out correct

GeosCore/flexchem_mod.F90:
- Update the loop where KPP variables MW, SR_MW, HENRY_K0, HENRY_CR
  are defined to skip "dummy" placeholder species

KPP/aerochem/gckpp_Sulfate_mod.F90
- Make sure is_HMS > 0 before using it
- Uncomment call to HET_DROP_CHEM
- Make sure FC > 0 before entering an IF block where it is
  found in a denominator

run/shared/species_database.yml:
- Add dummy species with Is_Omitted: true tags

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

This can be undone later.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

… into feature/msl/kpp

We will restore this later

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

aerochem.eqn:
  - Remove leftover ALKCLD and ALKDST{1-7} species.  These are for the
    acid dust uptake but are not needed here

gckpp_Util.F90: R
  - Rebuilt with KPP 2.3.2_gc

run/GCClassic/input.geos.templates/input.geos.fullchem
   - Added HMS as an advected species

run/shared/species_database.yml:
   - Remove DUMMY, HpA, HpC species, these were removed from aerochem.eqn
   - Add entries for HMS, HSO3m, SO3mm

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

(1) Only remove *.F90 files that are generated by KPP before calling
    the kpp executable.  This will allow us to leave other files
    (such as include files) untouched.

(2) Throw an error and exit if the gckpp.kpp file is not found.

(3) Update comments

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

gckpp.kpp
   - Add State_Het%wetArea(14) field to gckpp.kpp global definitions.
   - Add Buf_Cap_SALAAL and Buf_Cap_SALCAL

flexchem_mod.F90:
   - State_Het%wetArea is defined from State_Chm%WetAeroArea
     in Set_Kpp_Gridbox_Values (in flexchem_mod.F90)
   - Use KPP indices ind_HSO3m and ind_SO3mm in flexchem_mod.F90 instead
     of using the Ind_() function, which is much slower.
   - Define Buf_Cap_SALAAL = MW(ind_SALAAL) * 7.0e-5
   - Define Buf_Cap_SALCAL = MW(ind_SALCAL) * 7.0e-5

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

This should have been done at the last commit but evidently the
changes never got staged.

gckpp.kpp
   - Add Buf_Cap_SALAAL and Buf_Cap_SALACL

flexchem_mod.F90
   - Set Buf_Cap_SALAAL = MW(ind_SALAAL) * 7.0e-5_dp
   - Set Buf_Cap_SALCAL = MW(ind_SALCAL) * 7.0e-5_dp

gckpp_Global.F90:
   - Rebuild with KPP 2.3.2_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

… into feature/msl/kpp

The length of the rate law function library (included in the
#INLINE F90_RATES section of gckpp.kpp) was becoming too long
for KPP to inline into gckpp_Rates.F90.

As a solution, we have moved all of the rate law functions to a new
F90 module (aerochem_RateLawFuncs.F90).  In gckpp.kpp, we place a
"USE aerochem_RateLawFuncs" statement within the #INLINE F90_RCONST
block.  In this way, the rate law functions are directly included
in the subroutine Update_Rconst() in gckpp_Rates.F90.

We have rebuilt the KPP/aerochem/gckpp_*.F90 files with a preliminary
version of KPP 2.3.3_gc (currently in the dev branch of the KPP repo).

Also, we save the prior value of C(ind_SALAAL) and C(ind_SALCAL)
in State_Het%SALAAL_save and State_Het%SALCAL save.  This avoids
us having to use State_Chem%Species within the aerochem_RateLawFuncs.F90
module.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

KPP evidently has a limit as to how long a rate-law library that is
inlined with F90_RATES can be.  To get around this issue, we have
moved the rate-law functions to a module named aerochem_RateLawFuncs.F90,
which is inlined with F90_RCONST.  This will get referenced directly
from subroutine Update_RConst in gckpp_Rates.F90.

We can also move all of the the defined parameters in gckpp_Global.F90
(except for CON_R) can be moved to aerochem_RateLawFuncs.F90, where they
are now declared PRIVATE.  This will avoid namespace collisions.

NOTE: There seems to be a KPP parsing error that sometimes adds I2O4
to the end of lines.  Continuing to investigate.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

We have migrated the code that computes the hetchem reaction rate of
HOBr + BrSALA = Br2 from gckpp_HetRates.F90 into the new rate law
library module aerochem_RateLawFuncs.F90.  This involved the following
updates:

aerochem_RateLawFuncs.F90:
- add IMPLICIT NONE
- add functions Gam_HOBr_Aer and HOBrUptkByBrSALAandTropCloud
- cosmetic & formatting updates

aerochem.eqn:
- Remove obsolete rxns
- Replace HET(ind_HOBr,5) with call to HOBrUptkByBrSALAandTropCloud

gckpp.kpp
- add Br_over_Cl_SALA and Br_over_Cl_SALC scalar fields to the
  HetState derived type definition

gckpp_HetRates.F90:
- In routine Halide_Conc, compute State_Het%Br_over_Cl_SALA field
- In routine Halide_Conc, compute State_Het%Br_over_Cl_SALC field

flexchem_mod.F90:
- Remove obsolete indices id_SALA and id_SALC

gckpp_Global.F90
gckpp_Rates.F90
- Rebuilt with KPP 2.3.3_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

The rate-law function for HOBr + BrSALC = Br2 has now been peeled out
of gckpp_HetRates and into aerochem_RateLawFuncs.F90:

aerochem_RateLawFuncs.F90:
- Add function HOBrUptkByBrSALCandTropCloud (which is a copy of
  function HOBrUptkByBrSALAandTropCloud) with various arguments updated

aerochem.eqn
- Replace HET(ind_HOBr,6) with call to HOBrUptkByBrSALCandTropCloud

gckpp_Rates.F90
- Rebuilt with KPP 2.3.3_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

Some of the variable names in State_Het used the _SALA and _SALC suffixes
(e.g. Br_conc_SALA, Cl_conc_SALA).  But this can be misleading because
in the case of Br_Conc_SALA, this is Br- concentration in sea salt
aerosol which comes from BrSALA.  Therefore, to make things more
generic and to also line up with the SSA and SSC parameters for
the aerosol types, we now use _SSA and _SSC suffices for fields
of the State_Het object.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

aerochem_RateLawFuncs.F90:
- Added function HOBrUptkBySALACLandTropCloud
- Added function HOBrUptkBySALCCLAndTropCloud
- Removed unused variables

aerochem.eqn
- Replace HET(ind_HOBr, 3) with call to HOBrUptkBySALACLAndTropCloud
- Replace HET(ind_HOBr, 4) with call to HOBrUptkBySALCCLAndTropCloud

gckpp_Rates.F90
- Rebuilt with KPP 2.3.3_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

….F90

We have created a new module to contain common functions that many
of the rate-law computations use.  These have been split off from
aerochem_RateLawFuncs.F90 into rateLawUtilFuncs.F90.  These include
the functions ARS_L1K, GCARR_*, CloudHet, kIIR1Ltd, ReactoDiff_Corr,
Safe_Div, and Is_Safe_Div.

Right now we are storing this in the aerochem folder but this could
be moved to another folder later on.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

k_HOBr_HSO3 is now k_HOBr_HSO3m    (HSO3m = HSO3-)

k_HOBr_HSO3_2 is now k_HOBr_SO3mm  (SO3mm = SO3--)

Also removed brIce, dummy, area variables where not needed.

Also, when computing the branching ratios in the HOBrUptk* functions,
reorder the logic so that the function call to Br2_Yield will only
occur if the Br/Cl- ratio is <= 5.0e-4.  This avoids unnecessary
function calls.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

gckpp.kpp
- In the State_Het derived type definition:
  (1) mHNO4 is now HNO4_molal
  (2) mNO3 iis now NO3_molal
  (3) mSO4 is now SO4_molal
  (4) Removed fupdateHOBr
  (5) Removed fupdateHOCl

flexchem_mod.F90
- Use _dp instead of _f8
- Assign HSO4_molal, NO3_molal, NO4_molal fields of State_Het
- Remove assignment of fupdateHOBr, fupdateHOCl fields of State_Het
(4) Updated comments

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

aerochem_RateLawFuncs.F90:
- Add function HOBrUptkByHSO3mAndTropCloud
- Add function HOBrUptkBySO3mmAndTropCloud

aerochem.eqn:
- Replace HET(ind_HOBr,7) with call to HOBrUptkByHSO3mAndTropCloud
- Replace Het(ind_HOBr,8) with call to HOBrUptkBySO3mmAndTropCloud

gckpp_Global.F90
gckpp_Rates.F90
- Rebuilt with KPP 2.3.2_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

aerochem_RateLawFuncs.F90:
- Add function Gam_ClNO3_Aer
- Add function Gam_ClNO3_Ice
- Add Function ClNO3uptkByH2OandTropCloud

aerochem.eqn
- Replace HET(ind_ClNO3,1) with call to ClNO3uptkByH2OandTropCloud

gckpp_HetRates.F90
- Define State_Het%frac_SALACL in routine Halide_Conc

gckpp_Global.F90
gckpp_Rates.F90
- Rebuilt with KPP 2.3.3_gc

NOTE: Still need to pull out the function calls from the
old HetClNO3_TCld routine...

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

aerochem_RateLawFuncs.F90
- Added function ClNO3uptkByHClAndTropCloud, but disabled the
  algorithm from the hold HetClNO3_TCld, which produces unexplained
  numerical diffs.

aerochem.eqn
- Replace Het(ind_ClNO3,2) with ClNO3uptkByHClAndTropCloud

gckpp_Rates.F90
- Rebuilt with KPP 2.3.3_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

aerochem_RateLawFuncs.F90:
- Add function ClNO3uptkByHBrandTropCloud, but skip the
  stratopsheric computation for now.
- Add function ClNO3uptkByBrSALAandTropCloud, which is not yet
  activated

aerochem.eqn
- Replace Het(ind_ClNO3,3) with call to ClNO3uptkByHBrandTropCloud
- Replace het(ind_ClNO3,4) with 0.0_dp for now

gckpp_Rates.F90
- Rebuilt with KPP 2.3.2_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

aerochem_RateLawFuncs.F90:
- Add function HOClUptkByHBr

aerochem.eqn
- Replace HET(ind_HOCl,2) with call to HOClUptkByHBr

gckpp_Rates.F90
- Rebuilt with KPP 2.3.3_gc

flexchem_mod.F90:
- Explicitly copy State_Chm%KHETI_SLA(1:11) to State_Het%KHETI_SLA
  in order to avoid any issues

NOTE: Multiplying State_Het%KHETI_SLA(9) * State_Het%xArea(13)
introduces numerical differences.  Perhaps from a problem upstream.
Will look at this later.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

In GeosCore/ucx_mod.F90:

(1) In routine CALC_STRAT_AER, make sure that all PRIVATE variables
    are initialized or zeroed in the !$OMP PARALLEL loop.  This will
    help to prevent numerical instability issues.

(2) In routine CALC_SLA_GAMMA, set RXNGAMMA(2), RXNGAMMA(5),
    RXNGAMMA(9), and RXNGAMMA(11) to 0.0_f8 instead of TINY(1e0_f8).
    The TINY function returns the smallest positive number representable
    by REAL*8 (which is ~1e-308), which is very close to the range
    of denormal numbers.  Using 0.0_f8 will set these elements to
    exactly zero, which should eliminate any numerical noise or
    numerical instability.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

We had previously commmented out the term

  k = k + H%xArea(SLA) * H%KHETI_SLA(HOCl_plus_HBr)

because this was introducing numerical noise that could not be
matched to the reference code.  The problem was due to the field
State_Chm%KHETI_SLA(9) being set to TINY(1.0e_f8) in ucx_mod.F90.
The TINY function returns the smallest possible positive value for
REAL*8 in this case, which is ~1e-308.  This value is very close to
being denormal, and can cause numerical instability.  We fixed that
issue by setting State_Chm%KHETI_SLA(9) = 0.0_f8, which avoids the
instability.

With the instability fixed, we are able to restore this term to
function HOClUptkByHBr, and we have obtained identical results
to the reference code (with the same fix in ucx_mod.F90).

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

….F90

aerochem_RateLawFuncs.F90:
- Add tropopsheric cloud uptake (following the algorithm of the old
  HetClNO3_TCld routine) into ClNO3UptkByHBrAndTropCloud

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

We have added the algorithm for ClNO3 + HCl uptake in tropospheric ice
cloud to rate law function ClNO3uptkByHClAndTropCloud.

Previously, we had only implemented the stratospheric uptake in
this function, as we had seen numerical noise.  This numerical noise
has since been fixed.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

This will allow us to peel out rate-law functions one by one, in order
to better avoid errors.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

We have restarted peeling out the rate law functions for ClNO3 + HBr,
ClNO3 + HCl, and ClNO3 + BrSALA.  The ClNO3 + BrSALA function yields
some very small numerical differences in the clear-sky part of the
computation, this still needs to be investigated.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

aerochem.eqn:
- Replace HET(ind_ClNO3,4) by call to ClNO3uptkByBrSALA (prev commit)
- Replace HET(ind_ClNO3,5) by call to ClNO3uptkByBrSALC (prev commit)
- Replace HET(ind_ClNO3,6) by call to ClNO3uptkBySALACL
- Replace HET(ind_ClNO3,7) by call to ClNO3uptkBySALCCL

aerochem_RateLawFuncs.F90
- Add function ClNO3uptkByBrSALA (prev commit), partially disabled
- Add function ClNO3uptkByBrSALC (prev commit)
- Add function ClNO3uptkBySALACL, partially disabled
- Add function ClNO3uptkBySALCCL
- NOTE: "partially disabled" means that we have commented out the
  fine aerosol uptake parts of these functions, pending investigation
  of numerical differences

gckpp_HetRates.F90:
- In routine Get_Halide_SsaConc, don't multiply by 0.3333333, but instead
  divide by 3.0_dp in expression for v_tot.  This will avoid numerical
  precision issues.

gckpp_Rates.F90
- Rebuilt with KPP 2.3.3_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

aerochem.eqn:
- Replace HET(ind_HOCl,1) with call to HOClUptkByHClAndTropCloud
- Replace HET(ind_HOCl,2) with call to HOClUptkByHBr

aerochem_RateLawFuncs.F90:
- Add function Gam_HOCl_Cld
- Add function HOClUptkByHClAndTropCloud (partially implemented)
- Add function HOClUptkByHBr

gckpp_Rates.F90
- Rebuilt with KPP 2.3.3_gc

NOTE: The rate computations in tropospheric cloud are causing some
numerical differences.  These have been disabled, pending investigation.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

aerochem_RateLawFuncs.F90:
- Fix error in routine Gam_HOCl_Cld.  The H_X (Henry's law term) should
  have been included in the term for gb_tot, but was omitted.

This results in zero diffs w/r/t the reference code.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

We now compute branchCl = k_Cl / k_tot in the hetchem routine
Gam_HOCl_Cld instead of computing it externally.

This update gives zero diffs w/r/t the reference code.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

aerochem.eqn
- Add function call to HOClUptkBySALACLandTropCloud for HOCl + SALACL rxn

gckpp_Rates.F90
- Rebuilt with KPP 2.3.3_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

aerochem.eqn
- Add call to HOClUptkBySALCCLandTropCloud for HOCl + SALCCL rxn

gckpp_Rates.F90
- Rebuilt with KPP 2.3.3_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

gckpp.kpp
- Rename HSO3_conc_Cld to HSO3_aq_Cld in TYPE(HetState)
- Rename SO3_conc_Cld to SO3_aq_Cld in TYPE(HetState)
- Add TSO3_aq_Cld to TYPE(HetState) (will hold SO3 + HSO3 sum)

flexchem_mod.F90:
- Assign values to the HSO3_aq_Cld, SO3_aq_Cld, and TSO3_aq_Cld fields
  of the State_Het object

aerochem.eqn:
- Add call to HOClUptkByHSO3mAndTropCloud for HOCl + HSO3m rxn
- Add call to HOClUptkBySO3mmAndTropCloud for HOCl + SO3mm rxn

aerochem_RateLawFuncs.F90:
- Renamed H%HSO3_conc_Cld to H%HSO3_aq_Cld
- Renamed H%SO3_conc_Cld to H%SO3_aq_Cld
- Added function HOClUptkByHSO3mAndTropCloud
- Added function HOClUptkBySO3mmAndTropCloud
- Now pass branchSO3 out of subroutine Gam_HOCl_Cld

gckpp_Global.F90
gckpp_Rates.F90
- Rebuilt with KPP 2.3.3_gc

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

gckpp.kpp:
- Renamed Br_branch_{CldA,CldC,CldG} to frac_Br_{CldA,CldC,CldG}
- Renamed Cl_branch_{CldA,CldC,CldG} to frac_Cl_{CldA,CldC,CldG}
- Renamed HSO3_conc_Cld to HSO3_aq
- Renamed SO3_conc_Cld to SO3_aq
- Renamed TSO3_conc_Cld to TSO3_aq
- Added frac_HSO3_aq
- Added frac_SO3_aq
- Updated comments

aerochem_RateLawFuncs.F90:
gckpp_HetRates.F90
- Updates for consistency with State_Het fields

gckpp_Global.F90
gckpp_Rates.F90
- Rebuilt with KPP 2.3.3_gc

These updates result in zero diffs w/r/t the prior commit.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

We have simplified the names of most hetchem rate law functions by
removing "andTropCloud" or "AndTropCloud".  Updated both the
aerochem.eqn and aerochem_RateLawFuncs.F90 modules, and rebuilt
gckpp_Rates.F90 with KPP 2.3.3_gc.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

aerochem.eqn:
- Add function call to ClNO2uptkBySALACL for ClNO2+SALACL rxn
- Add function call to ClNO2uptkBySALCCL for ClNO2+SALCCL rxn
- Add function call to ClNO2uptkByHCl    for ClNO2+HCl rxn

aerochem_RateLawFuncs.F90:
- Add function ClNO2uptkBySALACL
- Add function ClNO2uptkBySALCCL
- Add function ClNO2uptkByHcl
- Add function ClNO2uptkByBrSALA (for next step)
- Remove H%is_UCX, it is redundant

flexchem_mod.F90:
- Updated comments

gckpp.kpp
- Updated comments

gckpp_Global.F90
gckpp_Rates.F90
- Rebuilt with KPP 2.3.3_gc

This update gives zero-diffs w/r/t the reference code.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

aerochem.eqn
- Add call to ClNO2uptkByBrSALA for ClNO2 + BrSALA hetchem rxn
- Add call to ClNO2uptkByBrSALC for ClNO2 + BrSALC hetchem rxn
- Add call to ClNO2uptkByHBr for ClNO2 + HBr hetchem rxn

aerochem_RateLawFuncs.F90:
- Add function ClNO2uptkByBrSALC
- Add function ClNO2uptkByHBr
- Now pass H%pHSSA(1) to the pH argument of ClNO2uptkByBrSALA
- Now pass H%pHSSA(1) to the pH argument of ClNO2uptkBySALACL

gckpp_Rates.F90
- Rebuilt with KPP 2.3.3_gc

This commit gives zero-diffs w/r/t the reference code

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

aerochem.eqn:
- Restore call to ClNO3uptkByH2O for ClNO3 + H2O hetchem rxn

aerochem_RateLawFuncs.F90:
- Temporarily comment out the ClNO3 + H2O uptake in trop cloud
  while investigating the source of numerical differences

gckpp_Rates.F90:
- Rebuilt with KPP 2.3.3_gc

This results in zero-diffs w/r/t the equivalent modifications in the
reference code.

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>