From 77dfcaeca14b091da4efd501366dc7399225e474 Mon Sep 17 00:00:00 2001 From: Lars Hoffmann Date: Mon, 30 Dec 2024 09:54:40 +0100 Subject: [PATCH] Added control parameter for depth of transition layer above the PBL. --- src/mptrac.c | 50 ++++++++++++++++++++++++-------------------------- src/mptrac.h | 38 ++++++++++++++++++++++---------------- 2 files changed, 46 insertions(+), 42 deletions(-) diff --git a/src/mptrac.c b/src/mptrac.c index 1618b5cb6..d5f104c4f 100644 --- a/src/mptrac.c +++ b/src/mptrac.c @@ -3064,8 +3064,14 @@ void module_diffusion_turb( /* Loop over particles... */ PARTICLE_LOOP(0, atm->np, 1, "acc data present(ctl,clim,atm,dt,rs,rs2)") { + /* Get PBL and surface pressure... */ + double pbl, ps; + INTPOL_INIT; + INTPOL_2D(pbl, 1); + INTPOL_2D(ps, 0); + /* Get weighting factors... */ - const double wpbl = pbl_weight(met0, met1, atm, ip); + const double wpbl = pbl_weight(ctl, atm, ip, pbl, ps); const double wtrop = tropo_weight(clim, atm, ip) * (1.0 - wpbl); const double wstrat = 1.0 - wpbl - wtrop; @@ -3091,32 +3097,27 @@ void module_diffusion_turb( /* Vertical mixing in the PBL... */ if (ctl->diff_mix_pbl) { - /* Get PBL pressure... */ - double pbl, ps, tpbl, ts; - INTPOL_INIT; - INTPOL_2D(pbl, 1); + /* Get top of transition layer... */ + const double pbl_ext = pbl - ctl->diff_pbl_trans * (ps - pbl); /* Check pressure... */ - if (atm->p[ip] >= pbl) { // TODO: add delta p to mix over inversion - - /* Get surface pressure... */ - INTPOL_2D(ps, 0); + if (atm->p[ip] >= pbl_ext) { /* Get temperature... */ + double tpbl_ext, ts; intpol_met_time_3d(met0, met0->t, met1, met1->t, atm->time[ip], ps, atm->lon[ip], atm->lat[ip], &ts, ci, cw, 1); intpol_met_time_3d(met0, met0->t, met1, met1->t, atm->time[ip], - pbl, atm->lon[ip], atm->lat[ip], &tpbl, ci, cw, 1); - - /* Get density... */ - const double rhos = ps / ts; - const double rhopbl = pbl / tpbl; + pbl_ext, atm->lon[ip], atm->lat[ip], &tpbl_ext, ci, + cw, 1); /* Get new air parcel density... */ - const double rho = rhos + (rhopbl - rhos) * rs2[ip]; + const double rhos = ps / ts; + const double rhopbl_ext = pbl_ext / tpbl_ext; + const double rho = rhos + (rhopbl_ext - rhos) * rs2[ip]; /* Get new air parcel pressure... */ - atm->p[ip] = LIN(rhos, ps, rhopbl, pbl, rho); + atm->p[ip] = LIN(rhos, ps, rhopbl_ext, pbl_ext, rho); } } } @@ -4426,19 +4427,14 @@ double nat_temperature( /*****************************************************************************/ double pbl_weight( - met_t *met0, - met_t *met1, + const ctl_t *ctl, const atm_t *atm, - const int ip) { - - /* Get PBL pressure... */ - double pbl, ps; - INTPOL_INIT; - INTPOL_2D(pbl, 1); - INTPOL_2D(ps, 1); + const int ip, + const double pbl, + const double ps) { /* Get pressure range... */ - const double p1 = pbl - 0.2 * (ps - pbl); + const double p1 = pbl - ctl->diff_pbl_trans * (ps - pbl); const double p0 = pbl; /* Get weighting factor... */ @@ -5431,6 +5427,8 @@ void read_ctl( ERRMSG("Set DIFFUSION to 0, 1 or 2!"); ctl->diff_mix_pbl = (int) scan_ctl(filename, argc, argv, "DIFF_MIX_PBL", -1, "0", NULL); + ctl->diff_pbl_trans + = scan_ctl(filename, argc, argv, "DIFF_PBL_TRANS", -1, "0", NULL); ctl->turb_dx_pbl = scan_ctl(filename, argc, argv, "TURB_DX_PBL", -1, "50", NULL); ctl->turb_dx_trop = diff --git a/src/mptrac.h b/src/mptrac.h index fd93bd7d9..c4f136766 100644 --- a/src/mptrac.h +++ b/src/mptrac.h @@ -2667,6 +2667,9 @@ typedef struct { /*! Vertical mixing in the PBL (0=off, 1=on). */ int diff_mix_pbl; + /*! Depth of PBL transition layer (fraction of PBL depth). */ + double diff_pbl_trans; + /*! Horizontal turbulent diffusion coefficient (PBL) [m^2/s]. */ double turb_dx_pbl; @@ -5808,28 +5811,31 @@ double nat_temperature( const double hno3); /** - * @brief Computes the weighting factor based on the planetary boundary layer (PBL) pressure. + * @brief Computes a weighting factor based on planetary boundary layer 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. + * This function calculates a weighting factor that determines the + * contribution of a pressure level to processes within the planetary + * boundary layer. The factor is based on the relative position of the + * pressure within a linear transition range defined by `pbl` and `ps`. * - * @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]. + * @param ctl Pointer to the control structure containing configuration parameters. + * @param atm Pointer to the atmospheric data structure containing pressure levels. + * @param ip Index of the pressure level in the atmospheric data array. + * @param pbl Pressure at the planetary boundary layer. + * @param ps Surface pressure. + * @return Weighting factor for the specified pressure level: + * - Returns 1.0 if the pressure is above the upper boundary (`p0`). + * - Returns 0.0 if the pressure is below the lower boundary (`p1`). + * - Returns a linearly interpolated value between 1.0 and 0.0 for pressures within the transition range. * * @author Lars Hoffmann */ double pbl_weight( - met_t * met0, - met_t * met1, + const ctl_t * ctl, const atm_t * atm, - const int ip); + const int ip, + const double pbl, + const double ps); /** * @brief Reads air parcel data from a specified file into the given atmospheric structure. @@ -7417,7 +7423,7 @@ double time_from_filename( const int offset); /** - * @brief Computes the weighting factor based on the tropopause pressure. + * @brief Computes a weighting factor based on tropopause pressure. * * This function calculates a weighting factor for a given pressure value in relation to * the tropopause pressure. The weighting factor is determined as follows: