Added option to set diffusivities in the PBL.

src/mptrac.c

@@ -2779,8 +2779,7 @@ void module_decay(
   PARTICLE_LOOP(0, atm->np, 1, "acc data present(ctl,clim,atm,dt)") {
 
     /* Get weighting factor... */
-    const double w =
-      tropo_weight(clim, atm->time[ip], atm->lat[ip], atm->p[ip]);
+    const double w = tropo_weight(clim, atm, ip);
 
     /* Set lifetime... */
     const double tdec = w * ctl->tdec_trop + (1 - w) * ctl->tdec_strat;
@@ -2898,6 +2897,8 @@ void module_diffusion_meso(
 void module_diffusion_turb(
   const ctl_t *ctl,
   const clim_t *clim,
+  met_t *met0,
+  met_t *met1,
   atm_t *atm,
   const double *dt) {
 
@@ -2918,13 +2919,19 @@ void module_diffusion_turb(
   /* Loop over particles... */
   PARTICLE_LOOP(0, atm->np, 1, "acc data present(ctl,clim,atm,dt,rs)") {
 
-    /* Get weighting factor... */
-    const double w =
-      tropo_weight(clim, atm->time[ip], atm->lat[ip], atm->p[ip]);
+    /* Get weighting factors... */
+    const double wpbl = pbl_weight(met0, met1, atm, ip);
+    const double wt = tropo_weight(clim, atm, ip) * (1 - wpbl);
 
     /* Set diffusivity... */
-    const double dx = w * ctl->turb_dx_trop + (1 - w) * ctl->turb_dx_strat;
-    const double dz = w * ctl->turb_dz_trop + (1 - w) * ctl->turb_dz_strat;
+    const double dx = wpbl * ctl->turb_dx_pbl + wt * ctl->turb_dx_trop
+      + (1 - wpbl - wt) * ctl->turb_dx_strat;
+    const double dz = wpbl * ctl->turb_dz_pbl + wt * ctl->turb_dz_trop
+      + (1 - wpbl - wt) * ctl->turb_dz_strat;
+
+
+    //printf("diffws: %g %g %g %g %g\n", Z(atm->p[ip]), wt, wpbl, dx, dz);
+
 
     /* Horizontal turbulent diffusion... */
     if (dx > 0) {
@@ -3526,8 +3533,7 @@ void module_mixing_help(
       /* Set mixing parameter... */
       double mixparam = 1.0;
       if (ctl->mixing_trop < 1 || ctl->mixing_strat < 1) {
-	double w =
-	  tropo_weight(clim, atm->time[ip], atm->lat[ip], atm->p[ip]);
+	double w = tropo_weight(clim, atm, ip);
 	mixparam = w * ctl->mixing_trop + (1 - w) * ctl->mixing_strat;
       }
 
@@ -4238,6 +4244,32 @@ double nat_temperature(
 
 /*****************************************************************************/
 
+double pbl_weight(
+  met_t *met0,
+  met_t *met1,
+  const atm_t *atm,
+  const int ip) {
+
+  /* Get PBL pressure... */
+  double pbl;
+  INTPOL_INIT;
+  INTPOL_2D(pbl, 1);
+
+  /* Get pressure range... */
+  const double p1 = pbl * 0.971832875032981;
+  const double p0 = pbl / 0.971832875032981;
+
+  /* Get weighting factor... */
+  if (atm->p[ip] > p0)
+    return 1;
+  else if (atm->p[ip] < p1)
+    return 0;
+  else
+    return LIN(p0, 1.0, p1, 0.0, atm->p[ip]);
+}
+
+/*****************************************************************************/
+
 int read_atm(
   const char *filename,
   const ctl_t *ctl,
@@ -5205,10 +5237,14 @@ void read_ctl(
     ERRMSG("Set ADVECT to 0, 1, 2, or 4!");
 
   /* Diffusion parameters... */
+  ctl->turb_dx_pbl =
+    scan_ctl(filename, argc, argv, "TURB_DX_PBL", -1, "0", NULL);
   ctl->turb_dx_trop =
     scan_ctl(filename, argc, argv, "TURB_DX_TROP", -1, "0", NULL);
   ctl->turb_dx_strat =
     scan_ctl(filename, argc, argv, "TURB_DX_STRAT", -1, "0", NULL);
+  ctl->turb_dz_pbl =
+    scan_ctl(filename, argc, argv, "TURB_DZ_PBL", -1, "0", NULL);
   ctl->turb_dz_trop =
     scan_ctl(filename, argc, argv, "TURB_DZ_TROP", -1, "0", NULL);
   ctl->turb_dz_strat =
@@ -7321,15 +7357,15 @@ void read_met_pbl(
 #pragma omp parallel for default(shared) collapse(2)
   for (int ix = 0; ix < met->nx; ix++)
     for (int iy = 0; iy < met->ny; iy++) {
-
+      
       /* Check minimum value... */
       double pbl_min =
 	met->ps[ix][iy] + DZ2DP(ctl->met_pbl_min, met->ps[ix][iy]);
       met->pbl[ix][iy] = MIN(met->pbl[ix][iy], (float) pbl_min);
-
+      
       /* Check maximum value... */
       double pbl_max =
-	met->ps[ix][iy] + DZ2DP(ctl->met_pbl_max, met->ps[ix][iy]);
+       met->ps[ix][iy] + DZ2DP(ctl->met_pbl_max, met->ps[ix][iy]);
       met->pbl[ix][iy] = MAX(met->pbl[ix][iy], (float) pbl_max);
     }
 }
@@ -8470,24 +8506,23 @@ double time_from_filename(
 
 double tropo_weight(
   const clim_t *clim,
-  const double t,
-  const double lat,
-  const double p) {
+  const atm_t *atm,
+  const int ip) {
 
   /* Get tropopause pressure... */
-  const double pt = clim_tropo(clim, t, lat);
+  const double pt = clim_tropo(clim, atm->time[ip], atm->lat[ip]);
 
   /* Get pressure range... */
   const double p1 = pt * 0.866877899;
   const double p0 = pt / 0.866877899;
 
   /* Get weighting factor... */
-  if (p > p0)
+  if (atm->p[ip] > p0)
     return 1;
-  else if (p < p1)
+  else if (atm->p[ip] < p1)
     return 0;
   else
-    return LIN(p0, 1.0, p1, 0.0, p);
+    return LIN(p0, 1.0, p1, 0.0, atm->p[ip]);
 }
 
 /*****************************************************************************/

src/mptrac.h

@@ -2634,12 +2634,18 @@ typedef struct {
   /*! Random number generator (0=GSL, 1=Squares, 2=cuRAND). */
   int rng_type;
 
+  /*! Horizontal turbulent diffusion coefficient (PBL) [m^2/s]. */
+  double turb_dx_pbl;
+
   /*! Horizontal turbulent diffusion coefficient (troposphere) [m^2/s]. */
   double turb_dx_trop;
 
   /*! Horizontal turbulent diffusion coefficient (stratosphere) [m^2/s]. */
   double turb_dx_strat;
 
+  /*! Vertical turbulent diffusion coefficient (PBL) [m^2/s]. */
+  double turb_dz_pbl;
+
   /*! Vertical turbulent diffusion coefficient (troposphere) [m^2/s]. */
   double turb_dz_trop;
 
@@ -5072,6 +5078,8 @@ void module_diffusion_meso(
 void module_diffusion_turb(
   const ctl_t * ctl,
   const clim_t * clim,
+  met_t * met0,
+  met_t * met1,
   atm_t * atm,
   const double *dt);
 
@@ -5685,6 +5693,30 @@ double nat_temperature(
   const double h2o,
   const double hno3);
 
+/**
+ * @brief Computes the weighting factor based on the planetary boundary layer (PBL) pressure.
+ *
+ * This function determines a weighting factor for a specific pressure value in relation to 
+ * the planetary boundary layer (PBL) pressure. The weighting factor is calculated as follows:
+ * - Returns 1 if the pressure is greater than a calculated upper limit.
+ * - Returns 0 if the pressure is less than a calculated lower limit.
+ * - Linearly interpolates between 1 and 0 within the range defined by the upper and lower limits.
+ *
+ * @param[in] met0 First meteorological input data structure.
+ * @param[in] met1 Second meteorological input data structure.
+ * @param[in] atm Pointer to the atmospheric data structure.
+ * @param[in] ip Index of the pressure value to evaluate within the atmospheric data.
+ * 
+ * @return Weighting factor (double) in the range [0, 1].
+ *
+ * @author Lars Hoffmann
+ */
+double pbl_weight(
+  met_t * met0,
+  met_t * met1,
+  const atm_t * atm,
+  const int ip);
+
 /**
  * @brief Reads air parcel data from a specified file into the given atmospheric structure.
  *
@@ -7271,37 +7303,26 @@ double time_from_filename(
   const int offset);
 
 /**
- * @brief Computes the weighting factor for a given pressure with respect to the tropopause.
+ * @brief Computes the weighting factor based on the tropopause pressure.
  *
- * The `tropo_weight` function calculates a weighting factor that
- * indicates how much a given pressure `p` is influenced by the
- * tropopause based on climatological data.
+ * This function calculates a weighting factor for a given pressure value in relation to 
+ * the tropopause pressure. The weighting factor is determined as follows:
+ * - Returns 1 if the pressure is greater than a calculated upper limit.
+ * - Returns 0 if the pressure is less than a calculated lower limit.
+ * - Linearly interpolates between 1 and 0 within the range defined by the upper and lower limits.
  *
- * @param clim A pointer to a `clim_t` structure containing climatological tropopause data.
- * @param t The time parameter, typically representing the time of year or specific temporal context.
- * @param lat The latitude for which the weighting factor is being calculated.
- * @param p The pressure for which the weighting factor is to be computed.
- * 
- * @return A double representing the weighting factor.
+ * @param[in] clim Pointer to the climatology data structure.
+ * @param[in] atm Pointer to the atmospheric data structure.
+ * @param[in] ip Index of the pressure value to evaluate within the atmospheric data.
  * 
- * The function performs the following steps:
- * - Calculates the tropopause pressure `pt` using the `clim_tropo` function.
- * - Determines the pressure range around the tropopause: `p1` (lower bound) and `p0` (upper bound).
- * - Computes the weighting factor based on the given pressure `p`:
- *   - Returns `1` if `p` is greater than `p0` (troposphere).
- *   - Returns `0` if `p` is less than `p1` (stratosphere).
- *   - Linearly interpolates between `1` and `0` for pressures between `p0` and `p1`.
- *
- * @note The `clim_tropo` function is assumed to provide the tropopause pressure based on climatological data.
- * @note The constants used in the pressure range calculation (`0.866877899` and its reciprocal) are specific to the function's logic.
+ * @return Weighting factor (double) in the range [0, 1].
  *
  * @author Lars Hoffmann
  */
 double tropo_weight(
   const clim_t * clim,
-  const double t,
-  const double lat,
-  const double p);
+  const atm_t * atm,
+  const int ip);
 
 /**
  * @brief Writes air parcel data to a file in various formats.

src/trac.c

@@ -232,7 +232,7 @@ int main(
 	  /* Turbulent diffusion... */
 	  if (ctl.turb_dx_trop > 0 || ctl.turb_dz_trop > 0
 	      || ctl.turb_dx_strat > 0 || ctl.turb_dz_strat > 0)
-	    module_diffusion_turb(&ctl, clim, atm, dt);
+	    module_diffusion_turb(&ctl, clim, met0, met1, atm, dt);
 
 	  /* Mesoscale diffusion... */
 	  if (ctl.turb_mesox > 0 || ctl.turb_mesoz > 0)

tests/trac_test/run.sh

@@ -44,6 +44,7 @@ BOUND_MASS = 0.0
 CONV_CAPE = 0.0
 H2O2_CHEM_REACTION = 1
 TRACER_CHEM = 1
+TURB_DX_PBL = 50.0
 TURB_DX_TROP = 50.0
 TURB_DZ_STRAT = 0.1
 TURB_MESOX = 0.16