Fix calls to random number generator in atm_init.

slcs-jsc · Mar 13, 2024 · d5ba28e · d5ba28e
1 parent a24a815
commit d5ba28e
Show file tree

Hide file tree

Showing 2 changed files with 27 additions and 19 deletions.
diff --git a/src/atm_init.c b/src/atm_init.c
@@ -87,21 +87,25 @@ int main(
 	  for (irep = 0; irep < rep; irep++) {
 
 	    /* Set position... */
-	    atm->time[atm->np]
-	      = (t + gsl_ran_gaussian_ziggurat(rng, st / 2.3548)
-		 + ut * (gsl_rng_uniform(rng) - 0.5));
-	    atm->p[atm->np]
-	      = P(z + gsl_ran_gaussian_ziggurat(rng, sz / 2.3548)
-		  + uz * (gsl_rng_uniform(rng) - 0.5));
-	    atm->lon[atm->np]
-	      = (lon + gsl_ran_gaussian_ziggurat(rng, slon / 2.3548)
-		 + gsl_ran_gaussian_ziggurat(rng, DX2DEG(sx, lat) / 2.3548)
-		 + ulon * (gsl_rng_uniform(rng) - 0.5));
+	    double rg = gsl_ran_gaussian_ziggurat(rng, st / 2.3548);
+	    double ru = ut * (gsl_rng_uniform(rng) - 0.5);
+	    atm->time[atm->np] = (t + rg + ru);
+
+	    rg = gsl_ran_gaussian_ziggurat(rng, sz / 2.3548);
+	    ru = uz * (gsl_rng_uniform(rng) - 0.5);
+	    atm->p[atm->np] = P(z + rg + ru);
+
+	    rg = gsl_ran_gaussian_ziggurat(rng, slon / 2.3548);
+	    double rx =
+	      gsl_ran_gaussian_ziggurat(rng, DX2DEG(sx, lat) / 2.3548);
+	    ru = ulon * (gsl_rng_uniform(rng) - 0.5);
+	    atm->lon[atm->np] = (lon + rg + rx + ru);
+
 	    do {
-	      atm->lat[atm->np]
-		= (lat + gsl_ran_gaussian_ziggurat(rng, slat / 2.3548)
-		   + gsl_ran_gaussian_ziggurat(rng, DY2DEG(sx) / 2.3548)
-		   + ulat * (gsl_rng_uniform(rng) - 0.5));
+	      rg = gsl_ran_gaussian_ziggurat(rng, slat / 2.3548);
+	      rx = gsl_ran_gaussian_ziggurat(rng, DY2DEG(sx) / 2.3548);
+	      ru = ulat * (gsl_rng_uniform(rng) - 0.5);
+	      atm->lat[atm->np] = (lat + rg + rx + ru);
 	    } while (even && gsl_rng_uniform(rng) >
 		     fabs(cos(atm->lat[atm->np] * M_PI / 180.)));
 

diff --git a/src/trac.c b/src/trac.c
@@ -1681,10 +1681,12 @@ void module_oh_chem(
       k = k0 * M / (1. + k0 * M / ki) * pow(0.6, 1. / (1. + c * c));
     }
 
-    /* Correction factor for high SO2 concentration...*/
+    /* Correction factor for high SO2 concentration... */
     double cor;
     if (ctl->qnt_Cx >= 0)
-      cor = atm->q[ctl->qnt_Cx][ip]>1.9e-9 ? 1.67971e-8 * pow(atm->q[ctl->qnt_Cx][ip], -0.891564) : 1;
+      cor =
+	atm->q[ctl->qnt_Cx][ip] >
+	1.9e-9 ? 1.67971e-8 * pow(atm->q[ctl->qnt_Cx][ip], -0.891564) : 1;
     else
       cor = 1;
     /* Calculate exponential decay... */
@@ -1749,11 +1751,13 @@ void module_h2o2_chem(
 
     /* Henry constant of H2O2... */
     double H_h2o2 = 8.3e2 * exp(7600 * (1 / t - 1 / 298.15)) * RI * t;
-    
-    /* Correction factor for high SO2 concentration...*/
+
+    /* Correction factor for high SO2 concentration... */
     double cor;
     if (ctl->qnt_Cx >= 0)
-      cor = atm->q[ctl->qnt_Cx][ip]>1.45e-9 ? 1.20717e-7 * pow(atm->q[ctl->qnt_Cx][ip], -0.782692) : 1;
+      cor =
+	atm->q[ctl->qnt_Cx][ip] >
+	1.45e-9 ? 1.20717e-7 * pow(atm->q[ctl->qnt_Cx][ip], -0.782692) : 1;
     else
       cor = 1;
     double h2o2 = H_h2o2