Revised OH chemistry code.

slcs-jsc · Jul 28, 2021 · 088ada9 · 088ada9
1 parent e6bf661
commit 088ada9
Show file tree

Hide file tree

Showing 3 changed files with 30 additions and 24 deletions.
diff --git a/src/libtrac.c b/src/libtrac.c
@@ -2455,7 +2455,7 @@ void read_ctl(
     ctl->tdec_strat =
       scan_ctl(filename, argc, argv, "TDEC_STRAT", -1, "0", NULL);
     ctl->oh_chem_reaction =
-      (int)scan_ctl(filename, argc, argv, "OH_CHEM_REACTION", -1, "0", NULL);
+      (int) scan_ctl(filename, argc, argv, "OH_CHEM_REACTION", -1, "0", NULL);
     for (int ip = 0; ip < 4; ip++)
       ctl->oh_chem[ip] =
 	scan_ctl(filename, argc, argv, "OH_CHEM", ip, "0", NULL);

diff --git a/src/libtrac.h b/src/libtrac.h
@@ -859,7 +859,10 @@ typedef struct {
   /*! Life time of particles (stratosphere)  [s]. */
   double tdec_strat;
 
-  /*! Coefficients for OH chemistry (k0, n, kinf, m). */
+  /*! Reaction type for OH chemistry (0=none, 2=bimolecular, 3=termolecular). */
+  int oh_chem_reaction;
+
+  /*! Coefficients for OH reaction rate (A, E/R or k0, n, kinf, m). */
   double oh_chem[4];
 
   /*! Coefficients for dry deposition (v). */
@@ -1039,9 +1042,6 @@ typedef struct {
   /*! Search radius around station [km]. */
   double stat_r;
 
-  /*! Reaction type of OH chemical module*/
-  int oh_chem_reaction;
-
 } ctl_t;
 
 /*! Atmospheric data. */

diff --git a/src/trac.c b/src/trac.c
@@ -378,7 +378,7 @@ int main(
 	module_decay(&ctl, atm, dt);
 
       /* OH chemistry... */
-      if (ctl.oh_chem_reaction!=0&&ctl.oh_chem[0] > 0 && ctl.oh_chem[2] > 0)
+      if (ctl.oh_chem_reaction != 0)
 	module_oh_chem(&ctl, met0, met1, atm, dt);
 
       /* Dry deposition... */
@@ -1313,27 +1313,33 @@ void module_oh_chem(
 
       /* Get temperature... */
       double t;
-      double k;
       INTPOL_INIT;
       INTPOL_3D(t, 1);
 
-      if (ctl->oh_chem_reaction == 3){
-
-      /* Calculate molecular density (IUPAC Data Sheet I.A4.86 SOx15)... */
-      double M = 7.243e21 * (atm->p[ip] / 1000.) / t;
-
-      /* Calculate rate coefficient for X + OH + M -> XOH + M
-         (JPL Publication 15-10) ... */
-      double k0 = ctl->oh_chem[0] *
-	(ctl->oh_chem[1] > 0 ? pow(t / 298., -ctl->oh_chem[1]) : 1.);
-      double ki = ctl->oh_chem[2] *
-	(ctl->oh_chem[3] > 0 ? pow(t / 298., -ctl->oh_chem[3]) : 1.);
-      double c = log10(k0 * M / ki);
-       k = k0 * M / (1. + k0 * M / ki) * pow(0.6, 1. / (1. + c * c));
-    } else if (ctl->oh_chem_reaction == 2){
-       k =   ctl->oh_chem[0] * exp(- ctl->oh_chem[1] *  t);    
-      } else
-      ERRMSG("No such reaction type!");
+      /* Calculate bimolecular reaction rate... */
+      double k;
+      if (ctl->oh_chem_reaction == 2)
+	k = ctl->oh_chem[0] * exp(-ctl->oh_chem[1] / t);
+
+      /* Calculate termolecular reaction rate... */
+      if (ctl->oh_chem_reaction == 3) {
+
+	/* Calculate molecular density (IUPAC Data Sheet I.A4.86 SOx15)... */
+	double M = 7.243e21 * (atm->p[ip] / 1000.) / t;
+
+	/* Calculate rate coefficient for X + OH + M -> XOH + M
+	   (JPL Publication 19-05) ... */
+	double k0 = ctl->oh_chem[0] *
+	  (ctl->oh_chem[1] > 0 ? pow(298. / t, ctl->oh_chem[1]) : 1.);
+	double ki = ctl->oh_chem[2] *
+	  (ctl->oh_chem[3] > 0 ? pow(298. / t, ctl->oh_chem[3]) : 1.);
+	double c = log10(k0 * M / ki);
+	k = k0 * M / (1. + k0 * M / ki) * pow(0.6, 1. / (1. + c * c));
+      }
+
+      /* Undefined... */
+      else
+	ERRMSG("No such reaction type!");
 
       /* Calculate exponential decay... */
       double aux