import os
import re
from os.path import exists, getsize
import numpy as np
import pandas as pd
[docs]
def create_mandchem(chemistry, mandatory_files):
"""Write the four mandatory TM5 chemistry text files from the YAML config.
Writes (in order):
1. ``REACTIONS.{schemeid}`` — one reaction string per line.
2. ``PRESCRIBED_SPECIES.{schemeid}`` — prescribed species names.
3. ``PRODLOSS_SPECIES.{schemeid}`` — production/loss species names.
4. ``DEPO_SPECIES.{schemeid}`` — depositing species names.
Files are created empty when the corresponding species list is absent.
Args:
chemistry: TM5 chemistry plugin instance.
mandatory_files (list[str]): four file paths in the order above.
"""
for mfile in mandatory_files:
os.system(f"touch {mfile}")
# Chemical reactions
with open(mandatory_files[0], "w") as f:
if chemistry.nreacs > 0:
for attr in chemistry.reactions.attributes:
reac = getattr(chemistry.reactions, attr)
f.write(reac + "\n")
if hasattr(chemistry, "prescrconcs"):
with open(mandatory_files[1], "w") as f:
for attr in chemistry.prescrconcs.attributes:
f.write(attr + "\n")
if hasattr(chemistry, "prodloss3d"):
with open(mandatory_files[2], "w") as f:
for attr in chemistry.prodloss3d.attributes:
f.write(attr + "\n")
if hasattr(chemistry, "deposition"):
with open(mandatory_files[3], "w") as f:
for attr in chemistry.deposition.attributes:
f.write(attr + "\n")
[docs]
def create_optchem(chemistry, filer, fileps):
"""Build the derived TM5 scheme files from reactions and species lists.
Reads ``REACTIONS`` and ``PRESCRIBED_SPECIES`` files, parses stoichiometry
and rate coefficients (same type codes as CHIMERE), then writes:
* ``STOICHIOMETRY`` / ``CHEMISTRY`` / ``REACTION_RATES``
* ``PHOTO_RATES`` / ``FAMILIES`` (empty)
* ``ALL_SPECIES`` / ``ACTIVE_SPECIES`` / ``ANTHROPIC`` / ``BIOGENIC``
Unlike the CHIMERE variant, the active-species file stores TM5-specific
``restart_id`` and ``mass`` columns instead of transport-scheme flags.
Args:
chemistry: TM5 chemistry plugin instance.
filer (str): path to the ``REACTIONS`` text file.
fileps (str): path to the ``PRESCRIBED_SPECIES`` text file.
Returns:
tuple[int, int]: ``(nallqmax, nphoto_rates)`` — total species count
and number of photolysis reactions.
"""
workdir = chemistry.workdir
dirchem_ref = chemistry.dirchem_ref
mecachim = chemistry.schemeid
files = f"{dirchem_ref}/STOICHIOMETRY.{mecachim}"
filec = f"{dirchem_ref}/CHEMISTRY.{mecachim}"
filerr = f"{dirchem_ref}/REACTION_RATES.{mecachim}"
filej = f"{dirchem_ref}/PHOTO_RATES.{mecachim}"
filef = f"{dirchem_ref}/FAMILIES.{mecachim}"
fileals = f"{dirchem_ref}/ALL_SPECIES.{mecachim}"
fileas = f"{dirchem_ref}/ACTIVE_SPECIES.{mecachim}"
fileanth = f"{dirchem_ref}/ANTHROPIC.{mecachim}"
filebio = f"{dirchem_ref}/BIOGENIC.{mecachim}"
os.system(f"rm -f {files} {filec} {filerr} {filef}")
# Read reactions
df_reac = pd.read_csv(
filer,
index_col=False,
header=None,
comment="#",
engine="python",
sep="\s+",
)
# Read prescribed species
prescribed_species = np.array([])
if exists(fileps) and getsize(fileps) > 0:
df_prescr = pd.read_csv(
fileps,
index_col=False,
header=None,
comment="#",
engine="python",
sep="\s+",
)
prescribed_species = df_prescr[0]
nonactive = np.append(prescribed_species, ["O2", "X", None])
# Create stoechiometry and chemistry
parts = df_reac[0].str.split("->", n=2, expand=True)
left_hand_sides = parts[0]
right_hand_sides = parts[1]
losses = left_hand_sides.str.split("+", expand=True)
nlosses = left_hand_sides.str.split("+").apply(len)
losses.insert(0, -1, nlosses)
prods = right_hand_sides.str.split("+", expand=True)
nprods = pd.Series(np.sum(~pd.isnull(prods).values, axis=1))
prods.insert(0, -1, nprods)
prods_spec = np.copy(prods.values)
stoichiometry = np.zeros((1, 6))
for (i, j), value in np.ndenumerate(prods.values[:, 1:]):
try:
prods_list = prods.values[i, j + 1].split("*")
if len(prods_list) > 1:
arr = np.array(
[
prods_list[1],
prods_list[0],
prods_list[0],
prods_list[0],
prods_list[0],
i + 1,
]
)
stoichiometry = np.append(stoichiometry, arr[np.newaxis, ...], axis=0)
if prods_list[-1] not in nonactive:
prods_spec[i, j + 1] = prods_list[-1]
else:
prods_spec[i, j + 1] = None
prods_spec[i, 0] -= 1
except AttributeError:
pass
file_chemistry = np.append(losses.values, prods_spec, axis=1)
file_chemistry = pd.DataFrame(file_chemistry)
stoichiometry = pd.DataFrame(stoichiometry[1:])
# Create reactions
reactions = df_reac[1].values
reactions_rates, nphoto_rates = read_react(reactions, filej)
reactions_rates = pd.DataFrame(reactions_rates)
reactions_rates = reactions_rates.astype({0: 'int32'})
# Active species
chemistry.nspecies = len(chemistry.acspecies.attributes)
# Create active_species (output_species) and all_species
output_species = np.array(chemistry.acspecies.attributes)
all_species = np.append(output_species, prescribed_species, axis=0)
type_spec = np.array([["type"]])
acinfos = np.array([0, 0])
type_spec = np.broadcast_to(type_spec, (all_species.shape[0], 1))
acinfos = np.broadcast_to(acinfos, (output_species.shape[0], 2))
all_species = np.append(all_species[:, np.newaxis], type_spec, axis=1)
output_species = np.append(output_species[:, np.newaxis], acinfos, axis=1)
for i in range(output_species.shape[0]):
outspec = getattr(chemistry.acspecies, output_species[i, 0])
output_species[i, 1] = str(getattr(outspec, "restart_id"))
output_species[i, 2] = str(getattr(outspec, "mass"))
for i in range(all_species.shape[0]):
if i < chemistry.nspecies:
all_species[i, 1] = "ac"
else:
all_species[i, 1] = "pr"
all_species = pd.DataFrame(all_species)
output_species = pd.DataFrame(output_species)
# Create files to create
stoichiometry.to_csv(files, sep=" ", header=False, index=False)
file_chemistry.to_csv(filec, sep=" ", header=False, index=False)
reactions_rates.to_csv(filerr, sep=" ", header=False, index=False)
all_species.to_csv(fileals, sep=" ", header=False, index=False)
output_species.to_csv(fileas, sep=" ", header=False, index=False)
output_species.to_csv(fileanth, sep=" ", header=False, index=False)
output_species.to_csv(filebio, sep=" ", header=False, index=False)
return all_species.shape[0], nphoto_rates
[docs]
def read_react(reactions, filej):
"""Parse reaction-rate strings and write the ``PHOTO_RATES`` file (TM5 variant).
Identical to the CHIMERE :func:`read_react` but returns a 2-tuple
``(reactions_rates, nphoto_rates)`` without the photo_rates DataFrame
(TM5 uses a different photolysis interface).
Classifies each reaction into one of five rate types (same codes as
CHIMERE: 1=constant, 2=simplified Arrhenius, 3=full Arrhenius,
4=pressure, 5=photolysis) and writes photolysis index ↔ name pairs
to *filej*.
Args:
reactions (np.ndarray): 1-D array of rate-formula strings.
filej (str): path to the ``PHOTO_RATES`` output file.
Returns:
tuple[np.ndarray, int]: ``(reactions_rates, nphoto_rates)`` where
*reactions_rates* is a ``(nreac, 10)`` float array with the type
code in column 0 and kinetic parameters in columns 1–9.
"""
nreacts = reactions.shape[0]
reactions_rates = np.zeros((nreacts, 10))
reactions_rates.fill(None)
photo_rates = np.zeros((1, 2))
nphoto_rates = 0
for i, value in np.ndenumerate(reactions):
# Type 1 : Constant rate
if re.search(r"k=", value):
val_list = re.split(r"=", value)
reactions_rates[i, 0] = 1
reactions_rates[i, 1] = val_list[1]
# Type 2 : Simplified Arrhenius
if re.search(r"k\(T\)=Aexp\(-B\/T\)", value):
val_list = re.split(r"=|,", value)
reactions_rates[i, 0] = 2
reactions_rates[i, 1] = val_list[3]
reactions_rates[i, 2] = val_list[5]
# Type 3 : Complete Arrhenius
if re.search(r"k\(T\)=Aexp\(-B\/T\)\(300\/T\)\*\*N", value):
val_list = re.split(r"=|,", value)
reactions_rates[i, 0] = 3
reactions_rates[i, 1] = val_list[3]
reactions_rates[i, 2] = val_list[5]
reactions_rates[i, 3] = val_list[7]
# Type 4 : Relative pressure
if re.search(r"k\(P\)=A\(B\+C\*P/Pref\)", value):
val_list = re.split(r"=|,", value)
reactions_rates[i, 0] = 4
reactions_rates[i, 1] = val_list[3]
reactions_rates[i, 2] = val_list[5]
reactions_rates[i, 3] = val_list[7]
# Type 5 : Simple Photolysis
if re.search(r"J=", value):
val_list = re.split(r"=|,", value)
reactions_rates[i, 0] = 5
reactions_rates[i, 1] = val_list[1]
photo_rates = np.append(
photo_rates, [[str(i[0] + 1), "j" + val_list[1]]], axis=0
)
nphoto_rates += 1
photo_rates = pd.DataFrame(photo_rates[1:])
photo_rates.to_csv(filej, sep=" ", header=False, index=False)
return reactions_rates, nphoto_rates
