diff --git a/src/mptrac.c b/src/mptrac.c index 261b7d8b0..eb43ece58 100644 --- a/src/mptrac.c +++ b/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]); } /*****************************************************************************/ diff --git a/src/mptrac.h b/src/mptrac.h index 6b90db1ad..5d55521c4 100644 --- a/src/mptrac.h +++ b/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. diff --git a/src/trac.c b/src/trac.c index 1e5b61eec..579d333b8 100644 --- a/src/trac.c +++ b/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) diff --git a/tests/trac_test/run.sh b/tests/trac_test/run.sh index de2148dd5..cbcc4494e 100755 --- a/tests/trac_test/run.sh +++ b/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