AUSTAL (modifiziert)
pluris.c
/**
* Plume Rise Model PLURIS
*
* For theory see:
* Janicke, U., Janicke, L.: A three-dimensional plume rise model for
* dry and wet plumes, Atmospheric Environment 35 (2000), 877-890.
*
* Copyright (C) Janicke Consulting, Ueberlingen, Germany, 1998-2024
* email: info@janicke.de
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation; either version 2 of
* the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
**/
/*==================================================================== PLURIS.C
*
* 2016-12-07 uj 3.0.3 imported from lj
* 2017-05-16 uj 3.0.4 additional options adjust_step, min_us
* 2018-09-10 uj 3.1.0 additional parameter zq, take care for T<0 C and T>100 C
* 2019-01-16 uj 3.1.1 var-opt
* 2019-02-09 uj 3.1.2 const-opt = reversed var-opt (default is ust(t), ra(z))
* 2019-03-02 uj 3.1.3 security checks added
* 2021-03-01 uj 3.1.4 const-opt replaced by const-ust to reproduce BESMIN/BESMAX
* 2022-12-21 uj security check (no version change)
* 2024-01-17 uj possibility to modify break factor (no version change)
* 2024-01-17 uj 3.2.0 reduction of stack-tip downwash correction (rf)
* limit lq
*
*
*===========================================================================*/
#include
#include
#include
#include
#include
#include "pluris.h"
static int DEBUG = 0;
#define false 0
#define true 1
#define LOG fprintf(L, "PLR: "), fprintf
#define NSTORE 500
#ifndef PI
#define PI 3.14159265358979323846
#endif
float PlrMaxHeight = PLRMAXHEIGHT;
float PlrBreakFactor = PLRBREAKFACTOR; //-2024-01-17
int PlrStDownwash = 1;
int PlrVerbose = 0;
FILE *PlrMessage = NULL;
double PlrLqMax;
static char PRGNAME[] = "IBJpluris";
static char VERSION[] = "3.2.0";
static FILE *L;
static double TZERO = 273.15; // 0 Deg Celsius
static double TZEROT = 273.16; // 0 Deg Celsius
static double GRAV = 9.8066; // gravitational acceleration (m/s2)
static double RGAS = 8.314472; // gas constant (J/(mol K))
static double MDRY = 28.96546e-3; // molar mass of dry air (kg/mol)
static double MVAPOUR = 18.01528e-3; // molar mass of vapour (kg/mol)
static double RDRY; // gas constant of dry air
static double RVAPOUR; // gas constant of water vapour
static double LHEAT = 2454300.0; // latent heat (J/kg)
static double CPRESS = 1004.1; // specific heat capacity at constant pressure (J/(kgK))
static double CFAC = 3.; // factor to convert top-heat to Gaussian profile
static double ETR_ALPHA = 0.15; // entrainment constant alpha
static double ETR_BETA = 0.6; // entrainment constant beta
static double ETR_GAMMA = 0.38; // entrainment constant gamma
static double USTTIME = 120.;
static double GRAVOCPRESS;
static double LHEATOCPRESS;
static double RVORD;
static double LN10;
//bp static int NVAR = 13; no VLA for vs-compiler
static double hFinal;
static double xFinal;
static double tstore[NSTORE], hstore[NSTORE]; //-2015-12-30
static int useSimProfiles = false;
static char ID[256]; // Kurzinformation
static double HQ; // Hoehe der Quelle (m)
static double DQ; // Durchmesser der Quelle (m)
static double UQ; // Austrittsgeschwindigkeit (m/s)
static double TQ; // Temperatur der Abluft (C)
static double QQ; // Waermestrom der Quelle (MW)
static double ZQ; // Spezifische Feuchte (kg/kg)
static double RQ; // Relative Feuchte (?berschreibt sq) (1)
static double SQ; // Spezifischer Wasserdampfgehalt (kg/kg)
static double LQ; // Spezifischer Fluessigwassergehalt (kg/kg)
static double RF; // Reduktion der stack-tip downwash Korrektur (1) -2024-01-17
static double A1; // Austrittswinkel gegen die Horizontale (Grad)
static double A2; // Austrittswinkel gegen x-Achse (West-Ost) gegen Uhrzeigersinn (Grad)
static double CQ; // Anfangskonzentration (ME/m3)
static double FQ; // Quadrat der Froudezahl beim Quellaustritt (ueberschreibt tq) (1)
static double IQ; // Turbulenzintensität (1)
static double XQ; // X-Koordinate der Quelle/Mittelpunkt (m)
static double YQ; // Y-Koordinate der Quelle/Mittelpunkt (m)
static double CA; // Konzentration in der Umgebung (ME/m3)
static double DA; // Windrichtung gegen Nord im Uhrzeigersinn (Grad)
static double EF; // Reduktionsfaktor fuer die Entrainment-Funktion (1)
static double IA; // Turbulenzintensität (1)
static int NZ; // Anzahl der vertikalen Stuetzpunkte
static double PA; // Standarddruck in der untersten Messhoehe (Pa)
static double RA; // Relative Feuchte (1)
static double VS; // Sedimentationsgeschwindigkeit (m/s)
static double US; // Schubspannungsgeschwindigkeit (m/s)
static double SC; // Skalierung (m)
static double SD; // Skalierte Rechen-Schrittweite (sc)
static double SE; // Skalierter Endpunkt (sc)
static int SN; // Anzahl der Ausgabeschritte
//
static double ds; // calculation step
static double dp; // output step
static double* zv; // ambient profiles, mesh (m)
static double* uv; // ambient profiles, wind speed (m/s)
static double* dv; // ambient profiles, direction (against north, rad)
static double* xv; // ambient profiles, vx (m/s)
static double* yv; // ambient profiles, vx (m/s)
static double* wv; // ambient profiles, vz (m/s)
static double* iv; // ambient profiles, turbulence intensity (1)
static double* jv; // ambient profiles, velocity fluctuations (m/s)
static double* tv; // ambient profiles, temperature (K)
static double* pv; // ambient profiles, pressure (Pa)
static double* rv; // ambient profiles, relative humidity (1)
static double* qv; // ambient profiles, specific humidity (kg/kg)
static double* cv; // ambient profiles, concentration (MU/m3)
static double* sv; // ambient profiles, sigma_w (m/s)
static double* av; // ambient profiles, density (kg/m3)
static double xx; // plume, x-coordinate (m)
static double yy; // plume, y-coordinate (m)
static double zz; // plume, z-coordinate (m)
static double ss; // plume, s-coordinate (path along the plume axis) (m)
static double ll; // plume, l-coordinate (path along the horizontal projection) (m)
static double zt; // plume, transport time (s)
static double kk[3]; // plume, cartesian velocity components (m/s)
static double uu; // plume, average velocity (m/s)
static double jj; // plume, turbulent velocity fluctuations (m/s)
static double tt; // plume, temperature (K)
static double dd; // plume, density (kg/m3)
static double cc; // plume, concentration (MU/m3)
static double qi; // plume, specific humidity (kg/kg)
static double ei; // plume, specific liquid water content (kg/kg)
static double zi; // plume, specific total water content (kg/kg)
static double rr; // plume, radius (m)
static double aa; // plume, visible radius (m)
static double w1; // plume, angle to the horizontal (rad)
static double w2; // plume, angle against north (rad)
static double kl[3]; // ambient at plume height, cartesian velocity components (m/s)
static double ul; // ambient at plume height, wind speed (m/s)
static double jl; // ambient at plume height, turbulent velocity fluctuations (m/s)
static double tl; // ambient at plume height, temperature (K)
static double dl; // ambient at plume height, density (kg/m3)
static double cl; // ambient at plume height, concentration (MU/m3)
static double ql; // ambient at plume height, specific humidity (kg/kg)
static double el; // ambient at plume height, specific liquid water content (kg/kg)
static double zl; // ambient at plume height, specific total water content (kg/kg)
static double pl; // ambient at plume height, pressure (Pa)
static double sl; // ambient at plume height, vertical velocity fluctuations (m/s)
static double uH; // wind speed at source height (m/s) (internal)
//
// result for a particle model
//
static double ptlVq, ptlTs;
static int adjust_step = 0;
static int const_ust = 0;
static double min_us = 0;
static double ta_C = 10.;
//=======================================================================
static void init() {
RDRY = RGAS/MDRY;
RVAPOUR = RGAS/MVAPOUR;
GRAVOCPRESS = GRAV/CPRESS;
LHEATOCPRESS = LHEAT/CPRESS;
RVORD = RVAPOUR/RDRY;
LN10 = log(10.);
}
//=========================================================================
PLRPVL* _pvl(double ha, double ua, double da, double ta, double ia, double ra,
double ca) {
PLRPVL* p = (PLRPVL*) calloc(1, sizeof(PLRPVL));
p->ha = ha;
p->ua = ua;
p->da = da;
p->ta = ta;
p->ra = ra;
p->ca = ca;
return p;
}
PLRPRM* _prm(double sc, double sd, double se, int sn) {
PLRPRM* p = (PLRPRM*) calloc(1, sizeof(PLRPRM));
p->sc = sc;
p->sd = sd;
p->se = se;
p->sn = sn;
return p;
}
//=========================================================================
char* str_version(void) {
static char s[256];
sprintf(s, "%s", VERSION);
return s;
}
char* str_prm(PLRPRM prm) { //-2024-01-17
static char s[256];
sprintf(s, "Prm[sc=%5.1f, sd=%7.3f, sm=%6.1f, sn=%4d, adjust_step=%4d, min_us=%4.2f, const_ust=%d, bf=%1.2f]",
prm.sc, prm.sd, prm.se, prm.sn, prm.adjust_step, prm.min_us, prm.const_ust, prm.bf);
return s;
}
char* str_pvl(PLRPVL pvl) {
static char s[256];
sprintf(s, "Pvl[ha=%6.1f, ua=%8.3f, da=%8.3f, ta=%8.3f, ia=%6.4f, ra=%6.4f, ca=%9.3e]",
pvl.ha, pvl.ua, pvl.da, pvl.ta, pvl.ia, pvl.ra, pvl.ca);
return s;
}
char* str_src(PLRSRC src) { //-2024-01-17
static char s[256];
sprintf(s, "Src[hb=%5.1f, dm=%7.1f, ve=%4.1f, te=%5.1f, zq=%7.4f, sq=%7.4f, rh=%6.3f, lw=%7.4f, rf=%4.2f, us=%6.4f, lm=%6.0f, ta=%5.1f]",
src.hb, src.dm, src.ve, src.te, src.zq, src.sq, src.rq, src.lq, src.rf, src.us, src.lm, src.ta);
return s;
}
char* str_rsl(PLRRSL rsl) {
static char s[256];
sprintf(s, "Rsl[vr=%7.3f, tr=%5.2f, dr=%6.2f]", rsl.vr, rsl.tr, rsl.dr);
return s;
}
static char* str_plr(void) {
static char s[256];
sprintf(s, "Plr[tt=%5.1f, sq=%7.4f, lq=%7.4f]", tt-TZERO, qi, ei);
return s;
}
//==========================================================================
static void PrmSetDefaults() {
strcpy(ID, "Test");
HQ = 10.0;
DQ = 1.0;
UQ = 10.0;
TQ = 10.0;
QQ = -1.0;
ZQ = 0.0;
SQ = 0.0;
RQ = 0.0;
LQ = 0.0;
RF = 1.0; //-2024-01-17
A1 = 90.00;
A2 = 0.0;
CQ = 0.0;
FQ = -1.0;
IQ = 0.1;
XQ = 0.0;
YQ = 0.0;
CA = 1.0;
DA = 270.0;
EF = 1.0;
IA = 0.1;
NZ = 300;
PA = 101300.0;
RA = 0.70;
VS = 0.0000;
US = 0.0;
SC = 1.0;
SD = 0.01;
SE = 1000.0;
SN = 100;
}
/**
* returns the pressure for the specified height index
* @param ih the heigth index
* @return the pressure (Pa)
*/
static double getPressure(int ih) {
double tint = 0.;
double p;
int i;
//
if (ih < 1)
return PA;
for (i=1; i<=ih; i++)
tint += 0.5*(zv[i]-zv[i-1])*(1./tv[i] + 1./tv[i-1]);
p = PA*exp(-(GRAV/RDRY)*tint);
return p;
}
/**
* Saturation pressure for the specified temperature and pressure
* (implementation of the Goff-Gratch equations)
* (for simplicity TZERO=273.15 used instead of triple temperature 273.16)
*
* @param t the temperature (K)
* @return the saturation pressure
*/
static double getPs(double t) {
double tot0, t0ot, ps;
if (t > TZEROT + 100.0)
t = TZEROT + 100.;
tot0 = t/TZEROT;
t0ot = TZEROT/t;
if (t > TZEROT)
ps = LN10*(10.79574*(1-t0ot) + 0.000150474*(1-pow(10, -8.2969*(tot0-1)))
+ 0.00042873*(pow(10, 4.76955*(1-t0ot))-1))
- 5.02800*log(tot0)+ LN10*0.78614;
else
ps = LN10*(-9.09718*(t0ot-1) + 0.876793*(1-tot0)) - 3.56654*log(t0ot)
+ log(6.1071); //-2018-09-10
ps = exp(ps);
return ps*100;
}
/**
* Saturation humidity for the specified temperature and pressure
* using a more accurate relation.
*
* @param t the temperature (K)
* @param p the pressure (Pa)
* @param eta the specific liquid water content
* @return the specific saturation humidity (kg per kg wet air)
*/
static double getQs(double t, double p, double eta) {
double psp = getPs(t)/p;
double k = (1. - eta)/(1. - psp*(1. - 1./RVORD));
double qs = psp*k/RVORD;
// LOG("&&& getQs(%e, %e, %e) = %e, psp=%e, k=%e\n", t, p, eta, qs, psp, k);
return qs;
}
/**
* Returns the density of wet air
* @param t the temperature (K)
* @param q the specific humidity (kg per kg wet air)
* @param eta the liquid water content (kg per kg wet air)
* @param p the pressure (Pa)
* @return the density (kg/m3)
*/
static double getDensity(char *tag, double t, double q, double eta, double p) {
double zeta;
double d;
//
if (q > 1.) q = 1.;
if (q < 0.) q = 0.;
zeta = eta + q;
if (zeta > 1.) zeta = 1.;
if (zeta < 0.) zeta = 0.;
d = p/(RDRY*t*(1 - zeta + RVORD*q));
// if (DEBUG) LOG("getDensity(%s, %e, %e, %e, %e) = %e\n", tag, t, q, eta, p, d);
return d;
}
static double toRadians(double deg) {
return deg*PI/180;
}
static int defAmbPrf(int NZ, PLRPVL* prf) {
int i;
double qs0, ra0;
//
zv = calloc(NZ, sizeof(double));
uv = calloc(NZ, sizeof(double));
dv = calloc(NZ, sizeof(double));
xv = calloc(NZ, sizeof(double));
yv = calloc(NZ, sizeof(double));
wv = calloc(NZ, sizeof(double));
iv = calloc(NZ, sizeof(double));
jv = calloc(NZ, sizeof(double));
tv = calloc(NZ, sizeof(double));
pv = calloc(NZ, sizeof(double));
rv = calloc(NZ, sizeof(double));
qv = calloc(NZ, sizeof(double));
cv = calloc(NZ, sizeof(double));
sv = calloc(NZ, sizeof(double));
av = calloc(NZ, sizeof(double));
//
qs0 = getQs(ta_C+TZERO, PA, 0);
ra0 = prf[0].ra;
for (i=0; i zv[i] = prf[i].ha;
uv[i] = prf[i].ua;
dv[i] = prf[i].da;
dv[i] = toRadians(dv[i]);
xv[i] = -uv[i]*sin(dv[i]);
yv[i] = -uv[i]*cos(dv[i]);
wv[i] = 0.;
iv[i] = prf[i].ia;
jv[i] = sqrt(3)*iv[i]*uv[i]; //-2015-12-27
tv[i] = prf[i].ta + TZERO;
pv[i] = getPressure(i);
if (pv[i] < 0)
return -1;
double qs = getQs(tv[i], pv[i], 0);
qv[i] = ra0*qs0;
if (qv[i] > qs)
qv[i] = qs;
rv[i] = qv[i]/qs;
prf[i].ra = rv[i];
/* } */
cv[i] = prf[i].ca;
av[i] = getDensity("defAmbPrf", tv[i], qv[i], 0.0, pv[i]);
}
return 0;
}
static void setAmb(double x, double y, double z) {
int iz;
int outside = false;
if (z < zv[0]) {
iz = 0;
outside = true;
}
else if (z >= zv[NZ-1]) {
iz = NZ-1;
outside = true;
}
else {
for (iz=0; iz if ((z >= zv[iz]) && (z < zv[iz+1]))
break;
}
// LOG("&&& setAmb(%1.1f, %1.1f, %1.1f): qv[%d]=%e\n", x, y, z, iz, qv[iz]);
tl = tv[iz];
kl[0] = xv[iz];
kl[1] = yv[iz];
kl[2] = wv[iz];
jl = jv[iz];
pl = pv[iz];
cl = cv[iz];
ql = qv[iz];
sl = sv[iz];
dl = av[iz];
if (!outside) {
double d = (z - zv[iz])/(zv[iz+1] - zv[iz]);
tl += d*(tv[iz+1]-tv[iz]);
kl[0] += d*(xv[iz+1]-xv[iz]);
kl[1] += d*(yv[iz+1]-yv[iz]);
kl[2] += d*(wv[iz+1]-wv[iz]);
ql += d*(qv[iz+1]-qv[iz]);
jl += d*(jv[iz+1]-jv[iz]);
pl += d*(pv[iz+1]-pv[iz]);
cl += d*(cv[iz+1]-cv[iz]);
sl += d*(sv[iz+1]-sv[iz]);
dl += d*(av[iz+1]-av[iz]);
}
ul = sqrt(kl[0]*kl[0] + kl[1]*kl[1] + kl[2]*kl[2]);
zl = ql;
el = 0;
dl = getDensity("setAmb", tl, ql, el, pl); // not using av, because in subtile runs linear interpolation is too rough
}
/**
*
* @param hh enthalpy/cp
* @param zz specific total water content
* @param q0 specific humidity centre line
* @param t0 temperature centre line
* @return temperature
*/
static double getTmpImplicitly(double hh, double zz, double q0, double t0) {
// note: hh = enthalphy/cp;
// zz = specific total water content;
// what happens:
// Equation h = cp*T - hv*(zeta - qs(T));
// Clausius-C. qs = q0 exp(hv*(T-T0)/(R*T0*T0));
// Expansion of qs to second order -> quadratic equation for T;
double A = LHEAT/(RVAPOUR*t0*t0);
double AA = A*A;
double a2 = 0.5*LHEATOCPRESS*q0*AA;
double a1 = 1 + LHEATOCPRESS*q0*(A - AA*t0);
double a0 = LHEATOCPRESS*q0*(1 - A*t0 + 0.5*AA*t0*t0) - hh - LHEATOCPRESS*zz;
double arg = a1*a1 - 4*a2*a0;
//System.out.printf("### hh=%f, zz=%f, q0=%f, t0=%f, arg=%f\n", hh, zz, q0, t0, arg);
//if (true) System.exit(0);
if (arg < 0)
return -999;
return (-a1 + sqrt(arg))/(2.*a2);
}
/** calculates the plume temperature, liquid water content, and visible radius
*
* @param hs enthalpy/cp
* @param zeta specific total water content
* @return { eta, a, tmp )
*/
static int getEtaATmp(double hs, double zeta, double* res) {
double t0, q0, ti;
double cfac, hst, zst, tst, b, eint, a, dr, t;
double f, hr, zr;
double r;
//
// if the plume is dry, piece of cake
//
if (zeta <= 0.) {
res[0] = 0.; // eta
res[1] = 0.; // a
res[2] = hs; // tmp
return 0;
}
//
// top hat profiles
//
if (!useSimProfiles) {
t0 = tt;
q0 = getQs(t0, pl, 0.);
ti = getTmpImplicitly(hs, zeta, q0, t0);
if (ti >= 100+TZERO) { //-2018-09-10
res[0] = 0.;
res[1] = 0.;
res[2] = hs;
}
else {
res[0] = (ti > hs) ? (ti - hs)/LHEATOCPRESS : 0;
res[1] = (ti > hs) ? DBL_MAX : 0;
res[2] = (ti > hs) ? ti : hs;
}
return 0;
}
//
// similarity profiles
//
cfac = CFAC + exp(-ss/(0.5*DQ))*(1 - CFAC);
hst = cfac*(hs - tl); if (hst < 0) hst = 0;
zst = cfac*(zeta - zl); if (zst < 0) zst = 0;
tst = cfac*(tt - tl); if (tst < 0) tst = 0;
t0 = tst + tl;
q0 = getQs(t0, pl, 0.);
b = rr/sqrt(cfac);
//
// start from outside the plume, decrease the distance to the plume axis,
// find the visible plume radius, and integrate the liquid water content
//
eint = 0.;
a = 0.;
dr = b/20.;
t = t0;
// double ti = t0;
// double hr = tl + hst;
// double zr;
for (r=2*b; r>-dr; r-=dr) {
f = (b > 0.) ? exp(-r*r/(b*b)) : 1.;
hr = tl + hst*f;
zr = zl + zst*f;
ti = getTmpImplicitly(hr, zr, q0, t0);
if ((ti-hr) > 0. && a == 0.)
a = r;
if (a > 0.) {
if (ti > hr)
t = ti;
if (t > hr)
eint += (t-hr)*r;
}
}
res[0] = eint * 2*dr/(LHEATOCPRESS*rr*rr);
res[1] = a;
if (res[0] < 1.e-10) {
res[0] = 0.;
res[1] = 0.;
}
res[2] = hs + LHEATOCPRESS*res[0];
return 0;
}
static double getEntrainment() {
double sw = kk[0]*kl[0] + kk[1]*kl[1] + kk[2]*kl[2];
double wp = sw/uu - uu;
double wp2 = wp*wp;
double ws2 = fabs(ul*ul - sw*sw/(uu*uu));
double if2 = (GRAV*fabs(dl-dd)*rr)/(uu*dl);
double etr = ETR_GAMMA*if2;
if ((wp2 + ws2) > 0) {
etr += (0.5*ETR_ALPHA*wp2 + ETR_BETA*ws2)/sqrt(wp2+ws2);
}
//
// LOG("&&& getEntrainment: etr=%e, rr=%e\n", etr, rr);
//
etr *= EF*rr;
return etr;
}
static void derivs(double snew, double* v, double* dv) {
double res[3];
double v0i, ui, dt, r2, sw, eps;
//
// calculate current parameters
//
v0i = 1./v[0];
kk[0] = v[1]*v0i;
kk[1] = v[2]*v0i;
kk[2] = v[3]*v0i;
uu = sqrt(kk[0]*kk[0] + kk[1]*kk[1] + kk[2]*kk[2]);
ui = 1./uu;
dt = (snew - ss)/uu;
xx = v[8] + kk[0]*dt;
yy = v[9] + kk[1]*dt;
zz = v[10] + kk[2]*dt;
setAmb(xx, yy, zz);
ll = v[11];
zi = v[5]*v0i;
jj = sqrt(fabs(v[7]));
zt = v[12];
getEtaATmp(v[4]*v0i, v[5]*v0i, res);
ei = res[0];
aa = res[1];
tt = res[2];
qi = zi - ei;
if (qi < 0) qi = 0; //-2018-09-10
dd = getDensity("derivs", tt, qi, ei, pl);
rr = sqrt(fabs(v[0]*ui/dd));
if (aa == DBL_MAX)
aa = rr;
r2 = rr*rr;
cc = v[6]*ui/(r2);
//
// calculate current derivatives
//
if (dv != NULL) {
sw = (kk[0]*kl[0] + kk[1]*kl[1] + kk[2]*kl[2]);
eps = 2.*dl*getEntrainment();
dv[0] = eps; // AE 2001, Eq. 10
dv[1] = eps*kl[0]; // AE 2001, Eq. 11
dv[2] = eps*kl[1]; // AE 2001, Eq. 11
dv[3] = eps*kl[2] - r2*(dd-dl)*GRAV; // AE 2001, Eq. 11
dv[4] = -(r2*dd*uu)*GRAVOCPRESS*(kk[2]*ui)*(dl/dd) // AE 2001, Eq. 35
+ eps*(tl - LHEATOCPRESS*el);
dv[5] = eps*zl; // AE 2001, Eq. 34
dv[6] = eps*(cl/dl); // AE 2001, Eq. 19
dv[7] = (eps*v0i)*(uu*uu + ul*ul - 2*sw + jl*jl - jj*jj); // AE 2001, Eq. 18
dv[8] = kk[0]*ui;
dv[9] = kk[1]*ui;
dv[10] = (kk[2] - VS)*ui;
dv[11] = sqrt(kk[0]*kk[0] + kk[1]*kk[1])*ui;
dv[12] = ui;
//
// {
// int i;
// LOG("&&& derivs(ss=%e, snew=%e):\n", ss, snew);
// for (i=0; i<13; i++)
// LOG("&&& %02d v=%e, dv=%e\n", i, v[i], dv[i]);
// }
}
}
static void rk(double* y, double* dydx, int n, double x, double h, double* yout) {
int i;
double xh, h1, h2;
double* dym = calloc(n, sizeof(double));
double* dyt = calloc(n, sizeof(double));
double* yt = calloc(n, sizeof(double));
h1 = h*0.5;
h2 = h/6.;
xh = x+h1;
for (i=0; i yt[i] = y[i] + h1*dydx[i];
derivs(xh, yt, dyt);
for (i=0; i yt[i] = y[i] + h1*dyt[i];
derivs(xh, yt, dym);
for (i=0; i yt[i] = y[i] + h*dym[i];
dym[i] += dyt[i];
}
derivs(x+h, yt, dyt);
for (i=0; i yout[i] = y[i] + h2*(dydx[i] + dyt[i] + 2.*dym[i]);
free(dym);
free(dyt);
free(yt);
}
static int initialize(PLRSRC* src, int nz, PLRPVL* prf) {
int rc;
int i;
//
// set input parameters
//
NZ = nz;
TQ = src->te;
DQ = src->dm;
HQ = src->hb;
UQ = src->ve;
RQ = src->rq;
ZQ = src->zq;
SQ = src->sq;
LQ = src->lq;
RF = src->rf; //-2024-01-17
US = src->us;
ta_C = src->ta;
if (DQ <= 0 || UQ <= 0) //-2019-03-02
return -9;
//
// convert to internal units
//
TQ += TZERO;
A1 = toRadians(A1);
A2 = toRadians(A2);
DA = toRadians(DA);
//
QQ = PI*0.25*DQ*DQ*UQ*(TZERO/TQ)*1.36e-3*(TQ-10.-TZERO);
//
// set ambient profiles
//
rc = defAmbPrf(NZ, prf);
if (rc < 0)
return -1;
//
// set parameters
//
ds = SD*SC;
dp = SE*SC/SN;
xx = XQ;
yy = YQ;
zz = HQ;
ss = 0;
ll = 0;
zt = 0;
setAmb(xx, yy, zz);
hFinal = -1;
xFinal = -1;
ptlTs = -1;
ptlVq = -1;
uH = ul;
rr = DQ*0.5;
aa = 0.;
qi = 0;
ei = 0;
zi = 0;
w1 = A1;
w2 = A2;
cc = CQ;
if (FQ > 0.) {
TQ = tl/(1. - UQ*UQ/(FQ*GRAV*DQ*0.5));
if ((TQ < TZERO-20) || (TQ > TZERO+600))
return -2;
QQ = (TQ < 10+TZERO) ? 0 : PI*0.25*DQ*DQ*UQ*(TZERO/TQ)*1.36e-3*(TQ-10.-TZERO);
}
tt = TQ;
if (tt < TZERO)
return -3;
//
// set sq and lq
if (LQ < 0) LQ = 0;
if (LQ > 1) LQ = 1;
if (SQ < 0) SQ = 0;
if (SQ > 1) SQ = 1;
if (TQ >= 100+TZERO) {
LQ = 0;
RQ = 0;
if (ZQ > 0) {
SQ = ZQ/(1+ZQ);
}
}
else {
if (ZQ > 0) {
double zeta = ZQ/(1+ZQ);
double qs0 = getQs(tt, pl, 0);
LQ = (zeta > qs0) ? (zeta - qs0)/(1-qs0) : 0.;
SQ = (zeta > qs0) ? qs0*(1-LQ) : zeta;
}
else {
if (SQ == 0 && RQ > 0)
SQ = RQ * getQs(tt, pl, LQ);
}
}
if (SQ+LQ > 1) {
if (PlrVerbose > 1) LOG(L, "adjusting water content (lq=%1.6e->%1.6e, sq=%1.6e)\n", LQ, 1-SQ, SQ);
LQ = 1 - SQ;
}
//
// limit LQ -2024-01-17
if (LQ > PlrLqMax)
PlrLqMax = LQ;
if (LQ > PLRLQMAX)
LQ = PLRLQMAX;
//
qi = SQ;
ei = LQ;
zi = qi + ei;
//
//
aa = (ei > 0.) ? rr : 0.;
dd = getDensity("initialize", tt, qi, ei, pl);
FQ = (dl != dd && DQ != 0) ? UQ*UQ*dl/(fabs(dl-dd)*DQ*0.5*GRAV) : 1.e20;
uu = UQ;
jj = IQ * UQ;
kk[0] = fabs(cos(w1)*uu) * cos(w2);
kk[1] = fabs(cos(w1)*uu) * sin(w2);
kk[2] = sin(w1)*uu;
for (i=0; i<3; i++)
if (fabs(kk[i]) < 1.e-10)
kk[i] = 0.;
return 0;
}
static int checkFinalRise(int evaluate) {
double h, t, d, w, _us;
int i0, i;
double red = 1;
//
if (evaluate) {
if (xFinal < 0) {
if (PlrVerbose > 0) LOG(L, "can't determine final rise!\n");
return false;
}
if (hFinal <= 0.01) { //-2019-03-05
ptlTs = 0.;
ptlVq = 0.;
}
else if (xFinal > 0) {
//
// find half time //-2015-12-30
//
h = hFinal/2.;
i0 = -1;
for (i=0; i if (hstore[i] <= h && h < hstore[i+1]) {
i0 = i;
break;
}
}
if (i0 < 0) { //-2019-03-02
LOG(L, "*** Can't determine half rise, using t=10*log(2)\n");
t = 10.*log(2.);
}
else {
t = tstore[i0];
d = (h - hstore[i0])/(hstore[i0+1] - hstore[i0]);
t += d*(tstore[i0+1] - tstore[i0]);
}
if (PlrVerbose>3) {
for (i=0; i LOG(L, "%3d: hstore=%6.2f, tstore=%6.2f\n", i, hstore[i], tstore[i]);
if (hstore[i] < 0)
break;
}
LOG(L, "h/2=%6.2f at t=%6.2f\n", h, t);
}
ptlTs = t/log(2.);
ptlVq = hFinal/ptlTs;
}
else {
ptlTs = (UQ > 0) ? hFinal/UQ : hFinal/1.;
ptlVq = (UQ > 0) ? UQ : 1.;
}
//
// stacktip downwash
//
if (PlrStDownwash) {
red = fmin(1.0, (UQ/uH)/(1.5/(1. + 2*pow(FQ, -0.3333))));
if (red < 1.0 && RF >= 0. && RF < 1.) //-2024-01-17
red *= RF;
ptlVq *= red;
}
if (PlrVerbose>1 && red<1)
LOG(L, "downwash reduction = %1.2f\n", red);
if (PlrVerbose > 1)
LOG(L, "Vq=%1.2f; Ts=%1.2f; Dh=%1.2f;\n", ptlVq, ptlTs, ptlVq*ptlTs);
return true;
}
//
// done already?
if (xFinal >= 0)
return true;
//
// set final rise conditions
w = sqrt(uu*uu + ul*ul - 2*(kk[0]*kl[0] + kk[1]*kl[1] + kk[2]*kl[2]));
_us = (US < min_us) ? min_us : US;
if (!const_ust) {
_us *= pow(fmax(zt,USTTIME)/USTTIME, 0.18);
}
if (xFinal < 0. && w < PlrBreakFactor*_us) {
hFinal = zz - HQ;
xFinal = ll;
}
//
// ultimate zmax criterion
if (xFinal < 0. && (zz <= 0. || zz > PlrMaxHeight)) {
hFinal = (zz <= 0) ? 0. : zz - HQ; //-2019-03-02
xFinal = ll;
}
return (xFinal >= 0);
}
/*
* Run pluris to calculate plume rise.
*
* Arguments:
* src pointer to source definition
* NZ number levels in meteorological profile
* prf meteorological profile (array with NZ elements)
* prm calculation parameters
* rsl pointer to result: Vptl and Tptl
*
* Returns:
* 0 on success
* -1 on initialization error
*/
int run_pluris(PLRSRC* src, int NZ, PLRPVL* prf, PLRPRM* prm, PLRRSL* rsl) {
int rc, nstore, np;
//bp double var[NVAR];
//bp double dvds[NVAR];
double var[13]; //bp
double dvds[13]; //bp
double M, ps;
int i;
rsl->vr = 0; //-2024-01-17
rsl->tr = 0; //-2024-01-17
rsl->dr = 0; //-2024-01-17
//
L = (PlrMessage == NULL) ? stdout : PlrMessage;
int count = 0;
int step = 100;
init();
if (PlrVerbose > 1) {
LOG(L, "%s\n", str_src(*src));
LOG(L, "%s\n", str_prm(*prm));
}
PrmSetDefaults();
if (prm != NULL) {
if (prm->sc > 0)
SC = prm->sc;
if (prm->sd > 0)
SD = prm->sd;
if (prm->se > 0)
SE = prm->se;
if (prm->sn > 0)
SN = prm->sn;
if (prm->bf > 0) //-2024-01-17
PlrBreakFactor = prm->bf;
min_us = prm->min_us;
adjust_step = prm->adjust_step;
const_ust = prm->const_ust; //-2021-03-01
}
rc = initialize(src, NZ, prf);
if (PlrVerbose > 3) {
LOG(L, "profile:\n");
for (i=0; i LOG(L, "%3d %s\n", i, str_pvl(prf[i]));
}
if (rc < 0) {
if (PlrVerbose > 3)
LOG(L, "initialize() returns %d\n", rc);
return -1;
}
if (PlrVerbose > 1)
LOG(L, "%s\n", str_plr());
//
// set integration parameters to their initial values
//
M = var[0] = rr*rr*dd*uu; // M, mass flow density / PI
var[1] = M*kk[0]; // M*ux, momentum flow density / PI
var[2] = M*kk[1]; // M*uy, momentum flow density / PI
var[3] = M*kk[2]; // M*uz, momentum flow density / PI
var[4] = M*(tt-LHEATOCPRESS*ei); // M*(T-(hv/cp)*eta), enthalpy/cp flow density / PI
var[5] = M*zi; // M*zeta, water flow density / PI
var[6] = M*cc/dd; // M*c/rho, scalar quantity flow density / PI
var[7] = jj*jj; // sigma*sigma, turbulence
var[8] = xx; // x-coordinate of the plume
var[9] = yy; // y-coordinate of the plume
var[10] = zz; // z-coordinate of the plume
var[11] = 0.; // path along the horizontal projection
var[12] = 0.; // travel time
//
// do the calculation
//
nstore = 0;
tstore[nstore] = 0;
hstore[nstore] = 0;
np = 0;
ps = 0;
while (np < SN) {
if (ss < 0 || ds <= 0.) { //-2022-12-21
LOG(L, "unexpected error in plume rise calculation!\n");
return -7;
}
if (adjust_step > 0) {
ds = rr/adjust_step;
}
if (PlrVerbose>3 && count%step == 0) {
LOG(L, "ss=%10.2f, zz=%10.2f, ll=%10.2f, uu=%10.4f, tt=%5.1f, d/d0=%5.1f\n", ss, var[10], var[11], uu, var[12], 2*rr/DQ);
count = 0;
}
count++;
//
// advance one step ds
//
//bp for (i=0; i for (i=0; i<13; i++) dvds[i] = 0; //bp // initialization
if (zz > zv[NZ-1]) {
LOG(L, "plume axis exceeds available profile height (%1.1f m), calculation aborted!\n", zv[NZ-1]);
return -6;
}
setAmb(xx, yy, zz);
if (DEBUG)
//bp for (i=0; i for (i=0; i<13; i++) //bp
LOG(L, "A %02d %14.6e %14.6e\n", i, var[i], dvds[i]);
derivs(ss, var, dvds);
if (DEBUG)
//bp for (i=0; i for (i=0; i<13; i++) //bp
LOG(L, "B %02d %14.6e %14.6e\n", i, var[i], dvds[i]);
//bp rk(var, dvds, NVAR, ss, ds, var);
rk(var, dvds, 13, ss, ds, var); //bp
if (DEBUG)
//bp for (i=0; i for (i=0; i<13; i++) //bp
LOG(L, "C %02d %14.6e %14.6e\n", i, var[i], dvds[i]);
if (DEBUG && count>=1) break;
//
// transfer calculated values at ss+ds from var[] to individual parameters
//
derivs(ss, var, NULL);
ss += ds;
ps += ds;
//
if (checkFinalRise(false)) {
if (PlrVerbose > 3) LOG(L, "final rise reached, break.\n");
break;
}
if (zz < 0.) {
if (PlrVerbose > 3) LOG(L, "\nplume axis at ground, break.\n");
break;
}
if (nstore < NSTORE-3 && fabs(hstore[nstore] - (zz-HQ)) > DQ) {
nstore++;
tstore[nstore] = zt;
hstore[nstore] = zz - HQ;
}
if ((int)(ps/ds + 0.5)*ds >= dp) {
np++;
ps = 0;
}
}
nstore++;
tstore[nstore] = zt;
hstore[nstore] = zz - HQ;
nstore++;
tstore[nstore] = -1;
hstore[nstore] = -1;
//
if (PlrVerbose > 2)
LOG(L, "final rise: dh=%5.1f m at x=%4.0f m\n", hFinal, xFinal);
if (false == checkFinalRise(true))
return -2;
rsl->vr = ptlVq;
rsl->tr = ptlTs;
rsl->dr = ptlVq*ptlTs;
free(zv);
free(uv);
free(dv);
free(xv);
free(yv);
free(wv);
free(iv);
free(jv);
free(tv);
free(pv);
free(rv);
free(qv);
free(cv);
free(sv);
free(av);
return 0;
}
//**************************************************************************
#ifdef MAIN
static int parse_line(char *s, double** ppv) {
char *ptok, *ptail;
double *pv = NULL;
int i, m = 0;
//
if (s == NULL)
return 0;
ptok = s;
while (1) {
strtod(ptok, &ptail);
if (ptail == ptok)
break;
m++;
ptok = ptail;
}
if (m == 0)
return 0;
pv = (double*) calloc(m, sizeof(double));
*ppv = pv;
ptok = s;
for (i=0; i pv[i] = strtod(ptok, &ptail);
ptok = ptail;
}
return m;
}
static int read_dmna(char* fn, PLRSRC *src, PLRPVL **pprf, double *ve, double *ts) {
char name[512];
char line[1024];
PLRPVL* prf = NULL;
FILE *f;
int nv=0;
char *pl;
{
int i, n;
double* zv = NULL; // ha
double* uv = NULL; // ua
double* dv = NULL; // da
double* tv = NULL; // ta
double* rv = NULL; // ra
//
strcpy(name, fn);
strcat(name, "/IBJpluris.dmna");
f = fopen(name, "r");
while (NULL != (pl = fgets(line, 1024, f))) {
if (!strncmp(line, "--- profiles", 12))
break;
if (strlen(line) < 3)
continue;
if (!strncmp(line, "hq", 2))
src->hb = strtod(line+3, NULL);
else if (!strncmp(line, "dq", 2))
src->dm = strtod(line+3, NULL);
else if (!strncmp(line, "uq", 2))
src->ve = strtod(line+3, NULL);
else if (!strncmp(line, "tq", 2))
src->te = strtod(line+3, NULL);
else if (!strncmp(line, "rq", 2))
src->rh = strtod(line+3, NULL);
else if (!strncmp(line, "lq", 2))
src->lw = strtod(line+3, NULL);
else if (src->us <= 0 && !strncmp(line, "us ", 3))
src->us = strtod(line+3, NULL);
else if (!strncmp(line, "ve ", 3))
*ve = strtod(line+3, NULL);
else if (!strncmp(line, "ts ", 3))
*ts = strtod(line+3, NULL);
}
if (pl == NULL)
return -1;
//
while (NULL != (pl = fgets(line, 1024, f))) {
char *pn = NULL;
n = -1;
if (!strncmp(line, pn="zv ", 3))
n = parse_line(line+3, &zv);
else if (!strncmp(line, pn="uv ", 3))
n = parse_line(line+3, &uv);
else if (!strncmp(line, pn="dv ", 3))
n = parse_line(line+3, &dv);
else if (!strncmp(line, pn="tv ", 3))
n = parse_line(line+3, &tv);
else if (!strncmp(line, pn="rv ", 3))
n = parse_line(line+3, &rv);
if (n == 0) {
LOG(L, "no data for '%s'\n", pn);
return -2;
}
else if (n > 0 && n != nv) {
if (nv == 0)
nv = n;
else {
LOG(L, "'%s' has %d elements, expected: %d\n", pn, n, nv);
return -3;
}
}
} // while
fclose(f);
NZ = nv;
//
if (zv == NULL)
LOG(L, "missing 'zv'\n");
if (uv == NULL)
LOG(L, "missing 'uv'\n");
if (dv == NULL)
LOG(L, "missing 'dv'\n");
if (tv == NULL)
LOG(L, "missing 'tv'\n");
if (rv == NULL)
LOG(L, "missing 'rv'\n");
if (src->us <= 0)
LOG(L, "missing 'us'\n");
if (!zv || !uv || !dv || !tv || !rv || src->us <= 0)
return -4;
prf = calloc(nv, sizeof(PLRPVL));
for (i=0; i prf[i].ha = zv[i];
prf[i].ua = uv[i];
prf[i].da = dv[i];
prf[i].ta = tv[i];
prf[i].ia = 0.0;
prf[i].ra = rv[i];
prf[i].ca = 0;
}
}
*pprf = prf;
return nv;
}
int main( int argc, char *argv[] ) {
char dir[256];
PLRSRC src;
PLRPRM prm;
PLRPVL *prf = NULL;
PLRRSL rsl = { 0, 0 };
int i;
int nv;
int rc;
double ve=0, ts=0;
src.us = 0;
PlrVerbose = 3;
//
L = stdout;
LOG(L, "%s %s started\n", PRGNAME, VERSION);
if (argc > 1)
strcpy(dir, argv[1]);
else
*dir = 0;
nv = read_dmna(dir, &src, &prf, &ve, &ts);
LOG(L, "read_dmna(%s, ...) = %d\n", dir, nv);
if (nv <= 0) {
return 1;
}
// LOG(L, "%s\n", str_src(src));
// for (i=0; i// LOG(L, "%3d: %s\n", i, str_pvl(prf[i]));
prm.sc = src.dm;
prm.sd = 0.1; //0.01; //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
prm.se = 1000;
prm.sn = 100;
rc = run_pluris(&src, nv, prf, &prm, &rsl);
LOG(L, "run_pluris=%d, Vptl=%6.3f, Tptl=%6.3f, Hptl=%6.3f (DMNA: ve=%6.3f, ts=%6.3f, dh=%6.3f)\n",
rc, rsl.vr, rsl.tr, rsl.dr, ve, ts, ve*ts);
//
LOG(L, "%s finished\n", PRGNAME);
return 0;
}
#endif
//*****************************************************************************
* Plume Rise Model PLURIS
*
* For theory see:
* Janicke, U., Janicke, L.: A three-dimensional plume rise model for
* dry and wet plumes, Atmospheric Environment 35 (2000), 877-890.
*
* Copyright (C) Janicke Consulting, Ueberlingen, Germany, 1998-2024
* email: info@janicke.de
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation; either version 2 of
* the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
**/
/*==================================================================== PLURIS.C
*
* 2016-12-07 uj 3.0.3 imported from lj
* 2017-05-16 uj 3.0.4 additional options adjust_step, min_us
* 2018-09-10 uj 3.1.0 additional parameter zq, take care for T<0 C and T>100 C
* 2019-01-16 uj 3.1.1 var-opt
* 2019-02-09 uj 3.1.2 const-opt = reversed var-opt (default is ust(t), ra(z))
* 2019-03-02 uj 3.1.3 security checks added
* 2021-03-01 uj 3.1.4 const-opt replaced by const-ust to reproduce BESMIN/BESMAX
* 2022-12-21 uj security check (no version change)
* 2024-01-17 uj possibility to modify break factor (no version change)
* 2024-01-17 uj 3.2.0 reduction of stack-tip downwash correction (rf)
* limit lq
*
*
*===========================================================================*/
#include
#include
#include
#include
#include
#include "pluris.h"
static int DEBUG = 0;
#define false 0
#define true 1
#define LOG fprintf(L, "PLR: "), fprintf
#define NSTORE 500
#ifndef PI
#define PI 3.14159265358979323846
#endif
float PlrMaxHeight = PLRMAXHEIGHT;
float PlrBreakFactor = PLRBREAKFACTOR; //-2024-01-17
int PlrStDownwash = 1;
int PlrVerbose = 0;
FILE *PlrMessage = NULL;
double PlrLqMax;
static char PRGNAME[] = "IBJpluris";
static char VERSION[] = "3.2.0";
static FILE *L;
static double TZERO = 273.15; // 0 Deg Celsius
static double TZEROT = 273.16; // 0 Deg Celsius
static double GRAV = 9.8066; // gravitational acceleration (m/s2)
static double RGAS = 8.314472; // gas constant (J/(mol K))
static double MDRY = 28.96546e-3; // molar mass of dry air (kg/mol)
static double MVAPOUR = 18.01528e-3; // molar mass of vapour (kg/mol)
static double RDRY; // gas constant of dry air
static double RVAPOUR; // gas constant of water vapour
static double LHEAT = 2454300.0; // latent heat (J/kg)
static double CPRESS = 1004.1; // specific heat capacity at constant pressure (J/(kgK))
static double CFAC = 3.; // factor to convert top-heat to Gaussian profile
static double ETR_ALPHA = 0.15; // entrainment constant alpha
static double ETR_BETA = 0.6; // entrainment constant beta
static double ETR_GAMMA = 0.38; // entrainment constant gamma
static double USTTIME = 120.;
static double GRAVOCPRESS;
static double LHEATOCPRESS;
static double RVORD;
static double LN10;
//bp static int NVAR = 13; no VLA for vs-compiler
static double hFinal;
static double xFinal;
static double tstore[NSTORE], hstore[NSTORE]; //-2015-12-30
static int useSimProfiles = false;
static char ID[256]; // Kurzinformation
static double HQ; // Hoehe der Quelle (m)
static double DQ; // Durchmesser der Quelle (m)
static double UQ; // Austrittsgeschwindigkeit (m/s)
static double TQ; // Temperatur der Abluft (C)
static double QQ; // Waermestrom der Quelle (MW)
static double ZQ; // Spezifische Feuchte (kg/kg)
static double RQ; // Relative Feuchte (?berschreibt sq) (1)
static double SQ; // Spezifischer Wasserdampfgehalt (kg/kg)
static double LQ; // Spezifischer Fluessigwassergehalt (kg/kg)
static double RF; // Reduktion der stack-tip downwash Korrektur (1) -2024-01-17
static double A1; // Austrittswinkel gegen die Horizontale (Grad)
static double A2; // Austrittswinkel gegen x-Achse (West-Ost) gegen Uhrzeigersinn (Grad)
static double CQ; // Anfangskonzentration (ME/m3)
static double FQ; // Quadrat der Froudezahl beim Quellaustritt (ueberschreibt tq) (1)
static double IQ; // Turbulenzintensität (1)
static double XQ; // X-Koordinate der Quelle/Mittelpunkt (m)
static double YQ; // Y-Koordinate der Quelle/Mittelpunkt (m)
static double CA; // Konzentration in der Umgebung (ME/m3)
static double DA; // Windrichtung gegen Nord im Uhrzeigersinn (Grad)
static double EF; // Reduktionsfaktor fuer die Entrainment-Funktion (1)
static double IA; // Turbulenzintensität (1)
static int NZ; // Anzahl der vertikalen Stuetzpunkte
static double PA; // Standarddruck in der untersten Messhoehe (Pa)
static double RA; // Relative Feuchte (1)
static double VS; // Sedimentationsgeschwindigkeit (m/s)
static double US; // Schubspannungsgeschwindigkeit (m/s)
static double SC; // Skalierung (m)
static double SD; // Skalierte Rechen-Schrittweite (sc)
static double SE; // Skalierter Endpunkt (sc)
static int SN; // Anzahl der Ausgabeschritte
//
static double ds; // calculation step
static double dp; // output step
static double* zv; // ambient profiles, mesh (m)
static double* uv; // ambient profiles, wind speed (m/s)
static double* dv; // ambient profiles, direction (against north, rad)
static double* xv; // ambient profiles, vx (m/s)
static double* yv; // ambient profiles, vx (m/s)
static double* wv; // ambient profiles, vz (m/s)
static double* iv; // ambient profiles, turbulence intensity (1)
static double* jv; // ambient profiles, velocity fluctuations (m/s)
static double* tv; // ambient profiles, temperature (K)
static double* pv; // ambient profiles, pressure (Pa)
static double* rv; // ambient profiles, relative humidity (1)
static double* qv; // ambient profiles, specific humidity (kg/kg)
static double* cv; // ambient profiles, concentration (MU/m3)
static double* sv; // ambient profiles, sigma_w (m/s)
static double* av; // ambient profiles, density (kg/m3)
static double xx; // plume, x-coordinate (m)
static double yy; // plume, y-coordinate (m)
static double zz; // plume, z-coordinate (m)
static double ss; // plume, s-coordinate (path along the plume axis) (m)
static double ll; // plume, l-coordinate (path along the horizontal projection) (m)
static double zt; // plume, transport time (s)
static double kk[3]; // plume, cartesian velocity components (m/s)
static double uu; // plume, average velocity (m/s)
static double jj; // plume, turbulent velocity fluctuations (m/s)
static double tt; // plume, temperature (K)
static double dd; // plume, density (kg/m3)
static double cc; // plume, concentration (MU/m3)
static double qi; // plume, specific humidity (kg/kg)
static double ei; // plume, specific liquid water content (kg/kg)
static double zi; // plume, specific total water content (kg/kg)
static double rr; // plume, radius (m)
static double aa; // plume, visible radius (m)
static double w1; // plume, angle to the horizontal (rad)
static double w2; // plume, angle against north (rad)
static double kl[3]; // ambient at plume height, cartesian velocity components (m/s)
static double ul; // ambient at plume height, wind speed (m/s)
static double jl; // ambient at plume height, turbulent velocity fluctuations (m/s)
static double tl; // ambient at plume height, temperature (K)
static double dl; // ambient at plume height, density (kg/m3)
static double cl; // ambient at plume height, concentration (MU/m3)
static double ql; // ambient at plume height, specific humidity (kg/kg)
static double el; // ambient at plume height, specific liquid water content (kg/kg)
static double zl; // ambient at plume height, specific total water content (kg/kg)
static double pl; // ambient at plume height, pressure (Pa)
static double sl; // ambient at plume height, vertical velocity fluctuations (m/s)
static double uH; // wind speed at source height (m/s) (internal)
//
// result for a particle model
//
static double ptlVq, ptlTs;
static int adjust_step = 0;
static int const_ust = 0;
static double min_us = 0;
static double ta_C = 10.;
//=======================================================================
static void init() {
RDRY = RGAS/MDRY;
RVAPOUR = RGAS/MVAPOUR;
GRAVOCPRESS = GRAV/CPRESS;
LHEATOCPRESS = LHEAT/CPRESS;
RVORD = RVAPOUR/RDRY;
LN10 = log(10.);
}
//=========================================================================
PLRPVL* _pvl(double ha, double ua, double da, double ta, double ia, double ra,
double ca) {
PLRPVL* p = (PLRPVL*) calloc(1, sizeof(PLRPVL));
p->ha = ha;
p->ua = ua;
p->da = da;
p->ta = ta;
p->ra = ra;
p->ca = ca;
return p;
}
PLRPRM* _prm(double sc, double sd, double se, int sn) {
PLRPRM* p = (PLRPRM*) calloc(1, sizeof(PLRPRM));
p->sc = sc;
p->sd = sd;
p->se = se;
p->sn = sn;
return p;
}
//=========================================================================
char* str_version(void) {
static char s[256];
sprintf(s, "%s", VERSION);
return s;
}
char* str_prm(PLRPRM prm) { //-2024-01-17
static char s[256];
sprintf(s, "Prm[sc=%5.1f, sd=%7.3f, sm=%6.1f, sn=%4d, adjust_step=%4d, min_us=%4.2f, const_ust=%d, bf=%1.2f]",
prm.sc, prm.sd, prm.se, prm.sn, prm.adjust_step, prm.min_us, prm.const_ust, prm.bf);
return s;
}
char* str_pvl(PLRPVL pvl) {
static char s[256];
sprintf(s, "Pvl[ha=%6.1f, ua=%8.3f, da=%8.3f, ta=%8.3f, ia=%6.4f, ra=%6.4f, ca=%9.3e]",
pvl.ha, pvl.ua, pvl.da, pvl.ta, pvl.ia, pvl.ra, pvl.ca);
return s;
}
char* str_src(PLRSRC src) { //-2024-01-17
static char s[256];
sprintf(s, "Src[hb=%5.1f, dm=%7.1f, ve=%4.1f, te=%5.1f, zq=%7.4f, sq=%7.4f, rh=%6.3f, lw=%7.4f, rf=%4.2f, us=%6.4f, lm=%6.0f, ta=%5.1f]",
src.hb, src.dm, src.ve, src.te, src.zq, src.sq, src.rq, src.lq, src.rf, src.us, src.lm, src.ta);
return s;
}
char* str_rsl(PLRRSL rsl) {
static char s[256];
sprintf(s, "Rsl[vr=%7.3f, tr=%5.2f, dr=%6.2f]", rsl.vr, rsl.tr, rsl.dr);
return s;
}
static char* str_plr(void) {
static char s[256];
sprintf(s, "Plr[tt=%5.1f, sq=%7.4f, lq=%7.4f]", tt-TZERO, qi, ei);
return s;
}
//==========================================================================
static void PrmSetDefaults() {
strcpy(ID, "Test");
HQ = 10.0;
DQ = 1.0;
UQ = 10.0;
TQ = 10.0;
QQ = -1.0;
ZQ = 0.0;
SQ = 0.0;
RQ = 0.0;
LQ = 0.0;
RF = 1.0; //-2024-01-17
A1 = 90.00;
A2 = 0.0;
CQ = 0.0;
FQ = -1.0;
IQ = 0.1;
XQ = 0.0;
YQ = 0.0;
CA = 1.0;
DA = 270.0;
EF = 1.0;
IA = 0.1;
NZ = 300;
PA = 101300.0;
RA = 0.70;
VS = 0.0000;
US = 0.0;
SC = 1.0;
SD = 0.01;
SE = 1000.0;
SN = 100;
}
/**
* returns the pressure for the specified height index
* @param ih the heigth index
* @return the pressure (Pa)
*/
static double getPressure(int ih) {
double tint = 0.;
double p;
int i;
//
if (ih < 1)
return PA;
for (i=1; i<=ih; i++)
tint += 0.5*(zv[i]-zv[i-1])*(1./tv[i] + 1./tv[i-1]);
p = PA*exp(-(GRAV/RDRY)*tint);
return p;
}
/**
* Saturation pressure for the specified temperature and pressure
* (implementation of the Goff-Gratch equations)
* (for simplicity TZERO=273.15 used instead of triple temperature 273.16)
*
* @param t the temperature (K)
* @return the saturation pressure
*/
static double getPs(double t) {
double tot0, t0ot, ps;
if (t > TZEROT + 100.0)
t = TZEROT + 100.;
tot0 = t/TZEROT;
t0ot = TZEROT/t;
if (t > TZEROT)
ps = LN10*(10.79574*(1-t0ot) + 0.000150474*(1-pow(10, -8.2969*(tot0-1)))
+ 0.00042873*(pow(10, 4.76955*(1-t0ot))-1))
- 5.02800*log(tot0)+ LN10*0.78614;
else
ps = LN10*(-9.09718*(t0ot-1) + 0.876793*(1-tot0)) - 3.56654*log(t0ot)
+ log(6.1071); //-2018-09-10
ps = exp(ps);
return ps*100;
}
/**
* Saturation humidity for the specified temperature and pressure
* using a more accurate relation.
*
* @param t the temperature (K)
* @param p the pressure (Pa)
* @param eta the specific liquid water content
* @return the specific saturation humidity (kg per kg wet air)
*/
static double getQs(double t, double p, double eta) {
double psp = getPs(t)/p;
double k = (1. - eta)/(1. - psp*(1. - 1./RVORD));
double qs = psp*k/RVORD;
// LOG("&&& getQs(%e, %e, %e) = %e, psp=%e, k=%e\n", t, p, eta, qs, psp, k);
return qs;
}
/**
* Returns the density of wet air
* @param t the temperature (K)
* @param q the specific humidity (kg per kg wet air)
* @param eta the liquid water content (kg per kg wet air)
* @param p the pressure (Pa)
* @return the density (kg/m3)
*/
static double getDensity(char *tag, double t, double q, double eta, double p) {
double zeta;
double d;
//
if (q > 1.) q = 1.;
if (q < 0.) q = 0.;
zeta = eta + q;
if (zeta > 1.) zeta = 1.;
if (zeta < 0.) zeta = 0.;
d = p/(RDRY*t*(1 - zeta + RVORD*q));
// if (DEBUG) LOG("getDensity(%s, %e, %e, %e, %e) = %e\n", tag, t, q, eta, p, d);
return d;
}
static double toRadians(double deg) {
return deg*PI/180;
}
static int defAmbPrf(int NZ, PLRPVL* prf) {
int i;
double qs0, ra0;
//
zv = calloc(NZ, sizeof(double));
uv = calloc(NZ, sizeof(double));
dv = calloc(NZ, sizeof(double));
xv = calloc(NZ, sizeof(double));
yv = calloc(NZ, sizeof(double));
wv = calloc(NZ, sizeof(double));
iv = calloc(NZ, sizeof(double));
jv = calloc(NZ, sizeof(double));
tv = calloc(NZ, sizeof(double));
pv = calloc(NZ, sizeof(double));
rv = calloc(NZ, sizeof(double));
qv = calloc(NZ, sizeof(double));
cv = calloc(NZ, sizeof(double));
sv = calloc(NZ, sizeof(double));
av = calloc(NZ, sizeof(double));
//
qs0 = getQs(ta_C+TZERO, PA, 0);
ra0 = prf[0].ra;
for (i=0; i
uv[i] = prf[i].ua;
dv[i] = prf[i].da;
dv[i] = toRadians(dv[i]);
xv[i] = -uv[i]*sin(dv[i]);
yv[i] = -uv[i]*cos(dv[i]);
wv[i] = 0.;
iv[i] = prf[i].ia;
jv[i] = sqrt(3)*iv[i]*uv[i]; //-2015-12-27
tv[i] = prf[i].ta + TZERO;
pv[i] = getPressure(i);
if (pv[i] < 0)
return -1;
double qs = getQs(tv[i], pv[i], 0);
qv[i] = ra0*qs0;
if (qv[i] > qs)
qv[i] = qs;
rv[i] = qv[i]/qs;
prf[i].ra = rv[i];
/* } */
cv[i] = prf[i].ca;
av[i] = getDensity("defAmbPrf", tv[i], qv[i], 0.0, pv[i]);
}
return 0;
}
static void setAmb(double x, double y, double z) {
int iz;
int outside = false;
if (z < zv[0]) {
iz = 0;
outside = true;
}
else if (z >= zv[NZ-1]) {
iz = NZ-1;
outside = true;
}
else {
for (iz=0; iz
break;
}
// LOG("&&& setAmb(%1.1f, %1.1f, %1.1f): qv[%d]=%e\n", x, y, z, iz, qv[iz]);
tl = tv[iz];
kl[0] = xv[iz];
kl[1] = yv[iz];
kl[2] = wv[iz];
jl = jv[iz];
pl = pv[iz];
cl = cv[iz];
ql = qv[iz];
sl = sv[iz];
dl = av[iz];
if (!outside) {
double d = (z - zv[iz])/(zv[iz+1] - zv[iz]);
tl += d*(tv[iz+1]-tv[iz]);
kl[0] += d*(xv[iz+1]-xv[iz]);
kl[1] += d*(yv[iz+1]-yv[iz]);
kl[2] += d*(wv[iz+1]-wv[iz]);
ql += d*(qv[iz+1]-qv[iz]);
jl += d*(jv[iz+1]-jv[iz]);
pl += d*(pv[iz+1]-pv[iz]);
cl += d*(cv[iz+1]-cv[iz]);
sl += d*(sv[iz+1]-sv[iz]);
dl += d*(av[iz+1]-av[iz]);
}
ul = sqrt(kl[0]*kl[0] + kl[1]*kl[1] + kl[2]*kl[2]);
zl = ql;
el = 0;
dl = getDensity("setAmb", tl, ql, el, pl); // not using av, because in subtile runs linear interpolation is too rough
}
/**
*
* @param hh enthalpy/cp
* @param zz specific total water content
* @param q0 specific humidity centre line
* @param t0 temperature centre line
* @return temperature
*/
static double getTmpImplicitly(double hh, double zz, double q0, double t0) {
// note: hh = enthalphy/cp;
// zz = specific total water content;
// what happens:
// Equation h = cp*T - hv*(zeta - qs(T));
// Clausius-C. qs = q0 exp(hv*(T-T0)/(R*T0*T0));
// Expansion of qs to second order -> quadratic equation for T;
double A = LHEAT/(RVAPOUR*t0*t0);
double AA = A*A;
double a2 = 0.5*LHEATOCPRESS*q0*AA;
double a1 = 1 + LHEATOCPRESS*q0*(A - AA*t0);
double a0 = LHEATOCPRESS*q0*(1 - A*t0 + 0.5*AA*t0*t0) - hh - LHEATOCPRESS*zz;
double arg = a1*a1 - 4*a2*a0;
//System.out.printf("### hh=%f, zz=%f, q0=%f, t0=%f, arg=%f\n", hh, zz, q0, t0, arg);
//if (true) System.exit(0);
if (arg < 0)
return -999;
return (-a1 + sqrt(arg))/(2.*a2);
}
/** calculates the plume temperature, liquid water content, and visible radius
*
* @param hs enthalpy/cp
* @param zeta specific total water content
* @return { eta, a, tmp )
*/
static int getEtaATmp(double hs, double zeta, double* res) {
double t0, q0, ti;
double cfac, hst, zst, tst, b, eint, a, dr, t;
double f, hr, zr;
double r;
//
// if the plume is dry, piece of cake
//
if (zeta <= 0.) {
res[0] = 0.; // eta
res[1] = 0.; // a
res[2] = hs; // tmp
return 0;
}
//
// top hat profiles
//
if (!useSimProfiles) {
t0 = tt;
q0 = getQs(t0, pl, 0.);
ti = getTmpImplicitly(hs, zeta, q0, t0);
if (ti >= 100+TZERO) { //-2018-09-10
res[0] = 0.;
res[1] = 0.;
res[2] = hs;
}
else {
res[0] = (ti > hs) ? (ti - hs)/LHEATOCPRESS : 0;
res[1] = (ti > hs) ? DBL_MAX : 0;
res[2] = (ti > hs) ? ti : hs;
}
return 0;
}
//
// similarity profiles
//
cfac = CFAC + exp(-ss/(0.5*DQ))*(1 - CFAC);
hst = cfac*(hs - tl); if (hst < 0) hst = 0;
zst = cfac*(zeta - zl); if (zst < 0) zst = 0;
tst = cfac*(tt - tl); if (tst < 0) tst = 0;
t0 = tst + tl;
q0 = getQs(t0, pl, 0.);
b = rr/sqrt(cfac);
//
// start from outside the plume, decrease the distance to the plume axis,
// find the visible plume radius, and integrate the liquid water content
//
eint = 0.;
a = 0.;
dr = b/20.;
t = t0;
// double ti = t0;
// double hr = tl + hst;
// double zr;
for (r=2*b; r>-dr; r-=dr) {
f = (b > 0.) ? exp(-r*r/(b*b)) : 1.;
hr = tl + hst*f;
zr = zl + zst*f;
ti = getTmpImplicitly(hr, zr, q0, t0);
if ((ti-hr) > 0. && a == 0.)
a = r;
if (a > 0.) {
if (ti > hr)
t = ti;
if (t > hr)
eint += (t-hr)*r;
}
}
res[0] = eint * 2*dr/(LHEATOCPRESS*rr*rr);
res[1] = a;
if (res[0] < 1.e-10) {
res[0] = 0.;
res[1] = 0.;
}
res[2] = hs + LHEATOCPRESS*res[0];
return 0;
}
static double getEntrainment() {
double sw = kk[0]*kl[0] + kk[1]*kl[1] + kk[2]*kl[2];
double wp = sw/uu - uu;
double wp2 = wp*wp;
double ws2 = fabs(ul*ul - sw*sw/(uu*uu));
double if2 = (GRAV*fabs(dl-dd)*rr)/(uu*dl);
double etr = ETR_GAMMA*if2;
if ((wp2 + ws2) > 0) {
etr += (0.5*ETR_ALPHA*wp2 + ETR_BETA*ws2)/sqrt(wp2+ws2);
}
//
// LOG("&&& getEntrainment: etr=%e, rr=%e\n", etr, rr);
//
etr *= EF*rr;
return etr;
}
static void derivs(double snew, double* v, double* dv) {
double res[3];
double v0i, ui, dt, r2, sw, eps;
//
// calculate current parameters
//
v0i = 1./v[0];
kk[0] = v[1]*v0i;
kk[1] = v[2]*v0i;
kk[2] = v[3]*v0i;
uu = sqrt(kk[0]*kk[0] + kk[1]*kk[1] + kk[2]*kk[2]);
ui = 1./uu;
dt = (snew - ss)/uu;
xx = v[8] + kk[0]*dt;
yy = v[9] + kk[1]*dt;
zz = v[10] + kk[2]*dt;
setAmb(xx, yy, zz);
ll = v[11];
zi = v[5]*v0i;
jj = sqrt(fabs(v[7]));
zt = v[12];
getEtaATmp(v[4]*v0i, v[5]*v0i, res);
ei = res[0];
aa = res[1];
tt = res[2];
qi = zi - ei;
if (qi < 0) qi = 0; //-2018-09-10
dd = getDensity("derivs", tt, qi, ei, pl);
rr = sqrt(fabs(v[0]*ui/dd));
if (aa == DBL_MAX)
aa = rr;
r2 = rr*rr;
cc = v[6]*ui/(r2);
//
// calculate current derivatives
//
if (dv != NULL) {
sw = (kk[0]*kl[0] + kk[1]*kl[1] + kk[2]*kl[2]);
eps = 2.*dl*getEntrainment();
dv[0] = eps; // AE 2001, Eq. 10
dv[1] = eps*kl[0]; // AE 2001, Eq. 11
dv[2] = eps*kl[1]; // AE 2001, Eq. 11
dv[3] = eps*kl[2] - r2*(dd-dl)*GRAV; // AE 2001, Eq. 11
dv[4] = -(r2*dd*uu)*GRAVOCPRESS*(kk[2]*ui)*(dl/dd) // AE 2001, Eq. 35
+ eps*(tl - LHEATOCPRESS*el);
dv[5] = eps*zl; // AE 2001, Eq. 34
dv[6] = eps*(cl/dl); // AE 2001, Eq. 19
dv[7] = (eps*v0i)*(uu*uu + ul*ul - 2*sw + jl*jl - jj*jj); // AE 2001, Eq. 18
dv[8] = kk[0]*ui;
dv[9] = kk[1]*ui;
dv[10] = (kk[2] - VS)*ui;
dv[11] = sqrt(kk[0]*kk[0] + kk[1]*kk[1])*ui;
dv[12] = ui;
//
// {
// int i;
// LOG("&&& derivs(ss=%e, snew=%e):\n", ss, snew);
// for (i=0; i<13; i++)
// LOG("&&& %02d v=%e, dv=%e\n", i, v[i], dv[i]);
// }
}
}
static void rk(double* y, double* dydx, int n, double x, double h, double* yout) {
int i;
double xh, h1, h2;
double* dym = calloc(n, sizeof(double));
double* dyt = calloc(n, sizeof(double));
double* yt = calloc(n, sizeof(double));
h1 = h*0.5;
h2 = h/6.;
xh = x+h1;
for (i=0; i
derivs(xh, yt, dyt);
for (i=0; i
derivs(xh, yt, dym);
for (i=0; i
dym[i] += dyt[i];
}
derivs(x+h, yt, dyt);
for (i=0; i
free(dym);
free(dyt);
free(yt);
}
static int initialize(PLRSRC* src, int nz, PLRPVL* prf) {
int rc;
int i;
//
// set input parameters
//
NZ = nz;
TQ = src->te;
DQ = src->dm;
HQ = src->hb;
UQ = src->ve;
RQ = src->rq;
ZQ = src->zq;
SQ = src->sq;
LQ = src->lq;
RF = src->rf; //-2024-01-17
US = src->us;
ta_C = src->ta;
if (DQ <= 0 || UQ <= 0) //-2019-03-02
return -9;
//
// convert to internal units
//
TQ += TZERO;
A1 = toRadians(A1);
A2 = toRadians(A2);
DA = toRadians(DA);
//
QQ = PI*0.25*DQ*DQ*UQ*(TZERO/TQ)*1.36e-3*(TQ-10.-TZERO);
//
// set ambient profiles
//
rc = defAmbPrf(NZ, prf);
if (rc < 0)
return -1;
//
// set parameters
//
ds = SD*SC;
dp = SE*SC/SN;
xx = XQ;
yy = YQ;
zz = HQ;
ss = 0;
ll = 0;
zt = 0;
setAmb(xx, yy, zz);
hFinal = -1;
xFinal = -1;
ptlTs = -1;
ptlVq = -1;
uH = ul;
rr = DQ*0.5;
aa = 0.;
qi = 0;
ei = 0;
zi = 0;
w1 = A1;
w2 = A2;
cc = CQ;
if (FQ > 0.) {
TQ = tl/(1. - UQ*UQ/(FQ*GRAV*DQ*0.5));
if ((TQ < TZERO-20) || (TQ > TZERO+600))
return -2;
QQ = (TQ < 10+TZERO) ? 0 : PI*0.25*DQ*DQ*UQ*(TZERO/TQ)*1.36e-3*(TQ-10.-TZERO);
}
tt = TQ;
if (tt < TZERO)
return -3;
//
// set sq and lq
if (LQ < 0) LQ = 0;
if (LQ > 1) LQ = 1;
if (SQ < 0) SQ = 0;
if (SQ > 1) SQ = 1;
if (TQ >= 100+TZERO) {
LQ = 0;
RQ = 0;
if (ZQ > 0) {
SQ = ZQ/(1+ZQ);
}
}
else {
if (ZQ > 0) {
double zeta = ZQ/(1+ZQ);
double qs0 = getQs(tt, pl, 0);
LQ = (zeta > qs0) ? (zeta - qs0)/(1-qs0) : 0.;
SQ = (zeta > qs0) ? qs0*(1-LQ) : zeta;
}
else {
if (SQ == 0 && RQ > 0)
SQ = RQ * getQs(tt, pl, LQ);
}
}
if (SQ+LQ > 1) {
if (PlrVerbose > 1) LOG(L, "adjusting water content (lq=%1.6e->%1.6e, sq=%1.6e)\n", LQ, 1-SQ, SQ);
LQ = 1 - SQ;
}
//
// limit LQ -2024-01-17
if (LQ > PlrLqMax)
PlrLqMax = LQ;
if (LQ > PLRLQMAX)
LQ = PLRLQMAX;
//
qi = SQ;
ei = LQ;
zi = qi + ei;
//
//
aa = (ei > 0.) ? rr : 0.;
dd = getDensity("initialize", tt, qi, ei, pl);
FQ = (dl != dd && DQ != 0) ? UQ*UQ*dl/(fabs(dl-dd)*DQ*0.5*GRAV) : 1.e20;
uu = UQ;
jj = IQ * UQ;
kk[0] = fabs(cos(w1)*uu) * cos(w2);
kk[1] = fabs(cos(w1)*uu) * sin(w2);
kk[2] = sin(w1)*uu;
for (i=0; i<3; i++)
if (fabs(kk[i]) < 1.e-10)
kk[i] = 0.;
return 0;
}
static int checkFinalRise(int evaluate) {
double h, t, d, w, _us;
int i0, i;
double red = 1;
//
if (evaluate) {
if (xFinal < 0) {
if (PlrVerbose > 0) LOG(L, "can't determine final rise!\n");
return false;
}
if (hFinal <= 0.01) { //-2019-03-05
ptlTs = 0.;
ptlVq = 0.;
}
else if (xFinal > 0) {
//
// find half time //-2015-12-30
//
h = hFinal/2.;
i0 = -1;
for (i=0; i
i0 = i;
break;
}
}
if (i0 < 0) { //-2019-03-02
LOG(L, "*** Can't determine half rise, using t=10*log(2)\n");
t = 10.*log(2.);
}
else {
t = tstore[i0];
d = (h - hstore[i0])/(hstore[i0+1] - hstore[i0]);
t += d*(tstore[i0+1] - tstore[i0]);
}
if (PlrVerbose>3) {
for (i=0; i
if (hstore[i] < 0)
break;
}
LOG(L, "h/2=%6.2f at t=%6.2f\n", h, t);
}
ptlTs = t/log(2.);
ptlVq = hFinal/ptlTs;
}
else {
ptlTs = (UQ > 0) ? hFinal/UQ : hFinal/1.;
ptlVq = (UQ > 0) ? UQ : 1.;
}
//
// stacktip downwash
//
if (PlrStDownwash) {
red = fmin(1.0, (UQ/uH)/(1.5/(1. + 2*pow(FQ, -0.3333))));
if (red < 1.0 && RF >= 0. && RF < 1.) //-2024-01-17
red *= RF;
ptlVq *= red;
}
if (PlrVerbose>1 && red<1)
LOG(L, "downwash reduction = %1.2f\n", red);
if (PlrVerbose > 1)
LOG(L, "Vq=%1.2f; Ts=%1.2f; Dh=%1.2f;\n", ptlVq, ptlTs, ptlVq*ptlTs);
return true;
}
//
// done already?
if (xFinal >= 0)
return true;
//
// set final rise conditions
w = sqrt(uu*uu + ul*ul - 2*(kk[0]*kl[0] + kk[1]*kl[1] + kk[2]*kl[2]));
_us = (US < min_us) ? min_us : US;
if (!const_ust) {
_us *= pow(fmax(zt,USTTIME)/USTTIME, 0.18);
}
if (xFinal < 0. && w < PlrBreakFactor*_us) {
hFinal = zz - HQ;
xFinal = ll;
}
//
// ultimate zmax criterion
if (xFinal < 0. && (zz <= 0. || zz > PlrMaxHeight)) {
hFinal = (zz <= 0) ? 0. : zz - HQ; //-2019-03-02
xFinal = ll;
}
return (xFinal >= 0);
}
/*
* Run pluris to calculate plume rise.
*
* Arguments:
* src pointer to source definition
* NZ number levels in meteorological profile
* prf meteorological profile (array with NZ elements)
* prm calculation parameters
* rsl pointer to result: Vptl and Tptl
*
* Returns:
* 0 on success
* -1 on initialization error
*/
int run_pluris(PLRSRC* src, int NZ, PLRPVL* prf, PLRPRM* prm, PLRRSL* rsl) {
int rc, nstore, np;
//bp double var[NVAR];
//bp double dvds[NVAR];
double var[13]; //bp
double dvds[13]; //bp
double M, ps;
int i;
rsl->vr = 0; //-2024-01-17
rsl->tr = 0; //-2024-01-17
rsl->dr = 0; //-2024-01-17
//
L = (PlrMessage == NULL) ? stdout : PlrMessage;
int count = 0;
int step = 100;
init();
if (PlrVerbose > 1) {
LOG(L, "%s\n", str_src(*src));
LOG(L, "%s\n", str_prm(*prm));
}
PrmSetDefaults();
if (prm != NULL) {
if (prm->sc > 0)
SC = prm->sc;
if (prm->sd > 0)
SD = prm->sd;
if (prm->se > 0)
SE = prm->se;
if (prm->sn > 0)
SN = prm->sn;
if (prm->bf > 0) //-2024-01-17
PlrBreakFactor = prm->bf;
min_us = prm->min_us;
adjust_step = prm->adjust_step;
const_ust = prm->const_ust; //-2021-03-01
}
rc = initialize(src, NZ, prf);
if (PlrVerbose > 3) {
LOG(L, "profile:\n");
for (i=0; i
}
if (rc < 0) {
if (PlrVerbose > 3)
LOG(L, "initialize() returns %d\n", rc);
return -1;
}
if (PlrVerbose > 1)
LOG(L, "%s\n", str_plr());
//
// set integration parameters to their initial values
//
M = var[0] = rr*rr*dd*uu; // M, mass flow density / PI
var[1] = M*kk[0]; // M*ux, momentum flow density / PI
var[2] = M*kk[1]; // M*uy, momentum flow density / PI
var[3] = M*kk[2]; // M*uz, momentum flow density / PI
var[4] = M*(tt-LHEATOCPRESS*ei); // M*(T-(hv/cp)*eta), enthalpy/cp flow density / PI
var[5] = M*zi; // M*zeta, water flow density / PI
var[6] = M*cc/dd; // M*c/rho, scalar quantity flow density / PI
var[7] = jj*jj; // sigma*sigma, turbulence
var[8] = xx; // x-coordinate of the plume
var[9] = yy; // y-coordinate of the plume
var[10] = zz; // z-coordinate of the plume
var[11] = 0.; // path along the horizontal projection
var[12] = 0.; // travel time
//
// do the calculation
//
nstore = 0;
tstore[nstore] = 0;
hstore[nstore] = 0;
np = 0;
ps = 0;
while (np < SN) {
if (ss < 0 || ds <= 0.) { //-2022-12-21
LOG(L, "unexpected error in plume rise calculation!\n");
return -7;
}
if (adjust_step > 0) {
ds = rr/adjust_step;
}
if (PlrVerbose>3 && count%step == 0) {
LOG(L, "ss=%10.2f, zz=%10.2f, ll=%10.2f, uu=%10.4f, tt=%5.1f, d/d0=%5.1f\n", ss, var[10], var[11], uu, var[12], 2*rr/DQ);
count = 0;
}
count++;
//
// advance one step ds
//
//bp for (i=0; i
if (zz > zv[NZ-1]) {
LOG(L, "plume axis exceeds available profile height (%1.1f m), calculation aborted!\n", zv[NZ-1]);
return -6;
}
setAmb(xx, yy, zz);
if (DEBUG)
//bp for (i=0; i
LOG(L, "A %02d %14.6e %14.6e\n", i, var[i], dvds[i]);
derivs(ss, var, dvds);
if (DEBUG)
//bp for (i=0; i
LOG(L, "B %02d %14.6e %14.6e\n", i, var[i], dvds[i]);
//bp rk(var, dvds, NVAR, ss, ds, var);
rk(var, dvds, 13, ss, ds, var); //bp
if (DEBUG)
//bp for (i=0; i
LOG(L, "C %02d %14.6e %14.6e\n", i, var[i], dvds[i]);
if (DEBUG && count>=1) break;
//
// transfer calculated values at ss+ds from var[] to individual parameters
//
derivs(ss, var, NULL);
ss += ds;
ps += ds;
//
if (checkFinalRise(false)) {
if (PlrVerbose > 3) LOG(L, "final rise reached, break.\n");
break;
}
if (zz < 0.) {
if (PlrVerbose > 3) LOG(L, "\nplume axis at ground, break.\n");
break;
}
if (nstore < NSTORE-3 && fabs(hstore[nstore] - (zz-HQ)) > DQ) {
nstore++;
tstore[nstore] = zt;
hstore[nstore] = zz - HQ;
}
if ((int)(ps/ds + 0.5)*ds >= dp) {
np++;
ps = 0;
}
}
nstore++;
tstore[nstore] = zt;
hstore[nstore] = zz - HQ;
nstore++;
tstore[nstore] = -1;
hstore[nstore] = -1;
//
if (PlrVerbose > 2)
LOG(L, "final rise: dh=%5.1f m at x=%4.0f m\n", hFinal, xFinal);
if (false == checkFinalRise(true))
return -2;
rsl->vr = ptlVq;
rsl->tr = ptlTs;
rsl->dr = ptlVq*ptlTs;
free(zv);
free(uv);
free(dv);
free(xv);
free(yv);
free(wv);
free(iv);
free(jv);
free(tv);
free(pv);
free(rv);
free(qv);
free(cv);
free(sv);
free(av);
return 0;
}
//**************************************************************************
#ifdef MAIN
static int parse_line(char *s, double** ppv) {
char *ptok, *ptail;
double *pv = NULL;
int i, m = 0;
//
if (s == NULL)
return 0;
ptok = s;
while (1) {
strtod(ptok, &ptail);
if (ptail == ptok)
break;
m++;
ptok = ptail;
}
if (m == 0)
return 0;
pv = (double*) calloc(m, sizeof(double));
*ppv = pv;
ptok = s;
for (i=0; i
ptok = ptail;
}
return m;
}
static int read_dmna(char* fn, PLRSRC *src, PLRPVL **pprf, double *ve, double *ts) {
char name[512];
char line[1024];
PLRPVL* prf = NULL;
FILE *f;
int nv=0;
char *pl;
{
int i, n;
double* zv = NULL; // ha
double* uv = NULL; // ua
double* dv = NULL; // da
double* tv = NULL; // ta
double* rv = NULL; // ra
//
strcpy(name, fn);
strcat(name, "/IBJpluris.dmna");
f = fopen(name, "r");
while (NULL != (pl = fgets(line, 1024, f))) {
if (!strncmp(line, "--- profiles", 12))
break;
if (strlen(line) < 3)
continue;
if (!strncmp(line, "hq", 2))
src->hb = strtod(line+3, NULL);
else if (!strncmp(line, "dq", 2))
src->dm = strtod(line+3, NULL);
else if (!strncmp(line, "uq", 2))
src->ve = strtod(line+3, NULL);
else if (!strncmp(line, "tq", 2))
src->te = strtod(line+3, NULL);
else if (!strncmp(line, "rq", 2))
src->rh = strtod(line+3, NULL);
else if (!strncmp(line, "lq", 2))
src->lw = strtod(line+3, NULL);
else if (src->us <= 0 && !strncmp(line, "us ", 3))
src->us = strtod(line+3, NULL);
else if (!strncmp(line, "ve ", 3))
*ve = strtod(line+3, NULL);
else if (!strncmp(line, "ts ", 3))
*ts = strtod(line+3, NULL);
}
if (pl == NULL)
return -1;
//
while (NULL != (pl = fgets(line, 1024, f))) {
char *pn = NULL;
n = -1;
if (!strncmp(line, pn="zv ", 3))
n = parse_line(line+3, &zv);
else if (!strncmp(line, pn="uv ", 3))
n = parse_line(line+3, &uv);
else if (!strncmp(line, pn="dv ", 3))
n = parse_line(line+3, &dv);
else if (!strncmp(line, pn="tv ", 3))
n = parse_line(line+3, &tv);
else if (!strncmp(line, pn="rv ", 3))
n = parse_line(line+3, &rv);
if (n == 0) {
LOG(L, "no data for '%s'\n", pn);
return -2;
}
else if (n > 0 && n != nv) {
if (nv == 0)
nv = n;
else {
LOG(L, "'%s' has %d elements, expected: %d\n", pn, n, nv);
return -3;
}
}
} // while
fclose(f);
NZ = nv;
//
if (zv == NULL)
LOG(L, "missing 'zv'\n");
if (uv == NULL)
LOG(L, "missing 'uv'\n");
if (dv == NULL)
LOG(L, "missing 'dv'\n");
if (tv == NULL)
LOG(L, "missing 'tv'\n");
if (rv == NULL)
LOG(L, "missing 'rv'\n");
if (src->us <= 0)
LOG(L, "missing 'us'\n");
if (!zv || !uv || !dv || !tv || !rv || src->us <= 0)
return -4;
prf = calloc(nv, sizeof(PLRPVL));
for (i=0; i
prf[i].ua = uv[i];
prf[i].da = dv[i];
prf[i].ta = tv[i];
prf[i].ia = 0.0;
prf[i].ra = rv[i];
prf[i].ca = 0;
}
}
*pprf = prf;
return nv;
}
int main( int argc, char *argv[] ) {
char dir[256];
PLRSRC src;
PLRPRM prm;
PLRPVL *prf = NULL;
PLRRSL rsl = { 0, 0 };
int i;
int nv;
int rc;
double ve=0, ts=0;
src.us = 0;
PlrVerbose = 3;
//
L = stdout;
LOG(L, "%s %s started\n", PRGNAME, VERSION);
if (argc > 1)
strcpy(dir, argv[1]);
else
*dir = 0;
nv = read_dmna(dir, &src, &prf, &ve, &ts);
LOG(L, "read_dmna(%s, ...) = %d\n", dir, nv);
if (nv <= 0) {
return 1;
}
// LOG(L, "%s\n", str_src(src));
// for (i=0; i
prm.sc = src.dm;
prm.sd = 0.1; //0.01; //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
prm.se = 1000;
prm.sn = 100;
rc = run_pluris(&src, nv, prf, &prm, &rsl);
LOG(L, "run_pluris=%d, Vptl=%6.3f, Tptl=%6.3f, Hptl=%6.3f (DMNA: ve=%6.3f, ts=%6.3f, dh=%6.3f)\n",
rc, rsl.vr, rsl.tr, rsl.dr, ve, ts, ve*ts);
//
LOG(L, "%s finished\n", PRGNAME);
return 0;
}
#endif
//*****************************************************************************