AUSTAL (modified)
TalWrk.c
/* ==================================================================== TalWrk.c
*
* Particle driver for AUSTAL
* ==========================
*
* Copyright (C) Umweltbundesamt, Dessau-Roßlau, Germany, 2002-2024
* Copyright (C) Janicke Consulting, 88662 Überlingen, Germany, 2002-2024
* Email: info@austal.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.
*
* last change: 2024-01-17 uj
*
* ===========================================================================*/
#include
#ifdef __linux__ //-2021-07-03
#include
#endif
#define STDMYMAIN WrkMain
#include "IBJmsg.h"
#include "IBJary.h"
#include "IBJtxt.h"
#include "IBJdmn.h"
#include "IBJstd.h"
static char *eMODn = "TalWrk";
//=========================================================================
STDPGM(tstwrk, WrkServer, 3, 3, 0);
//=========================================================================
#include "genutl.h"
//#include "genio.h"
#include "TalWnd.h"
#include "TalPrf.h"
#include "TalGrd.h"
#include "TalMat.h"
#include "TalMod.h"
#include "TalNms.h"
//#include "TalPlm.h"
#include "TalPpm.h"
#include "TalPrm.h"
#include "TalPtl.h"
#include "TalSrc.h"
#include "TalTmn.h"
#include "TalDef.h" //-2011-12-05
#include "TalWrk.h"
#include "TalWrk.nls"
#include "pluris.h" //-2018-10-04
#include "TalBlm.h" //-2018-10-04
#define LOC 0x02000000L
#ifndef PI
#define PI 3.1415926
#endif
#define WRK_REFB 0x0001
#define WRK_REFT 0x0002
#define WRK_USE_H 0x0100
#define WRK_USE_AZ 0x0200
#define WRK_DRIFT_VDROP 3.6 //-2023-07-17
#define WRK_DRIFT_VDROPEXP 0.16 //-2023-07-17
#define WRK_DRIFT_VRED 0.8 //-2023-07-17
#define INSOLID(a,b,c) ((Pbn) && *(long*)AryPtrX(Pbn,(a),(b),(c))&GRD_VOLOUT)
void** WetDepoList; // needed in TALdos() //-2023-07-17
typedef struct {
int i, j, k, is, js, ks;
float ax, ay, az;
} LOCREC;
typedef struct xfrrec { int ir, ic; float val; //-2011-10-13
} XFRREC;
/*======================== Funktions-Prototypen ============================*/
static int ClcAfu( // Berechnung der Abgasfahnen-Ueberhoehung
SRCREC *ps, // Pointer auf die Quelle
float v2, // Geschwindigkeitsquadrat
float stab, // Stabilitaetsklasse oder Temperatur-Gradient
float *pu, // Ergebnis: Ueberhoehung.
float *pt ); // Ergebnis: Anstiegszeit.
static int ClcPtlLoc( // Teilchen im Netz einordnen, evtl. reflektieren
PTLREC *pa, // Pointer auf alte Teilchen-Koordinaten
PTLREC *pp, // Pointer auf neue Teilchen-Koordinaten
LOCREC *pl, // Pointer auf Index-Werte
int h2z ); // h-Wert nach z-Wert umwandeln
static int ClcLocMod( // Feld am Ort des Teilchens berechnen
LOCREC *pl, // Pointer auf Index-Werte
MD3REC *p3, // Pointer auf Feld-Werte
int clceta ); // <=0, wenn Eta nicht interpoliert werden soll
static int RunPtl( // Teilchen-Trajektorien berechnen.
void );
static int Try2Escape( // In einem Gebaeude erzeugte Partikel herausdraengen
PTLREC *pa,
LOCREC *pl );
/*=================== Alle Groessen nur intern fuer Worker ====================*/
static long start_sec = 0; //bp
static long Flags = 0;
static int Valid = 1;
static int Nptl = 0;
static int Ndsc = 0;
static int NumRun = 0;
static int NumRemove = 0;
static int MaxRetry = 50;
static int MinRetry = 6;
static int Nx, Ny, Nz, Nzm, PrfMode, DefMode, TrcMode, Nc;
static int OnlyAdvec, NoPredict, ExactTau, WetDrift; //-2023-07-17
static long T1, T2, DosT1, DosT2;
static long PrmT2, PpmT2, WdpT2, SrcT2, XfrT2, PtlT2, ModT2;
static float Dd, Xmin, Xmax, Gx, Ymin, Ymax, Gy, Xbmin, Xbmax, Ybmin, Ybmax;
static float Zref, Zscl, Stab, *Sk, Hm, Lm, Z0, D0, Us; //-2018-10-04
static float Ta, Rh; //-2018-10-04
static int LogPluris; //-2018-10-04
static int Gl, Gi, Gp;
static char AltName[256];
static FILE *TrcFile;
static GRDPARM *Pgp; /* Grid-Parameter */
static PRMINF *Ppi; /* General Parameter */
static ARYDSC *Pma; /* Model-Array */
static ARYDSC *Ppp; /* Particle-Parameter */
static ARYDSC *Ppa; /* Particle-Array */
static ARYDSC *Pda; /* Dose-Array */
static ARYDSC *Pwa; /* Wdep-Array */
static ARYDSC *Pxa; /* Xfr-Array (Chem) */
static ARYDSC *Psa; /* Source-Array */
static ARYDSC *Pbn; /* Boundaries-Array */
static int bPerx, bPery, bDdmet;
static int TotalTime, TotalDrop, AverTime, AverDrop, InterTime, InterDrop, DropT1;
static char PhaseSymbols[] = "|/-\\";
static int CountStep, CountPhase, CountMax=50000;
static int ShowProgress, Progress=-1; //-2001-07-10
static int Nxfr; //-2011-10-13
static XFRREC *Pxfr;
static int *Ixfr;
static float precep; //-2023-07-17
static float vdrop; //-2023-07-17
static int numnet; //-2023-07-17
static double WetXmin, WetXmax, WetYmin, WetYmax; //-2023-07-17
//-------------------------------------------------------------
static unsigned long WrkRndSeed = 123456;
static int WrkRndIset = 0;
static unsigned long WrkRndULong() {
WrkRndSeed = 1664525*WrkRndSeed + 1013904223;
return WrkRndSeed;
}
static unsigned long WrkRndGetSeed() {
return WrkRndSeed;
}
static void WrkRndSetSeed(unsigned long seed) {
if (seed != 0) WrkRndSeed = seed;
else WrkRndSeed = 123456;
WrkRndIset = 0;
}
static float WrkRndFloat() {
unsigned long r = WrkRndULong();
float f = ((r>>8) & 0x00ffffff)/((float)0x01000000);
return f;
}
#define WrkRnd01() WrkRndFloat() // random numbers
#define WrkRnd() WrkRndNrm() // use normal distribution
//================================================================== WrkRndNrm
//
static float WrkRndNrm() // normal distribution, average=0, variance=1 -2008-12-11
{
static float gset;
float fac, r, v1, v2;
if (WrkRndIset == 0) {
do {
r = WrkRndFloat();
v1 = 2.0*r - 1.0;
r = WrkRndFloat();
v2 = 2.0*r - 1.0;
r = v1*v1 + v2*v2;
} while (r >= 1.0);
fac = sqrt( -2.0*log(r)/r );
gset = v1*fac;
WrkRndIset = 1;
return v2*fac;
}
else {
WrkRndIset = 0;
return gset;
}
}
static double WrkClock( // Get the clock in seconds
double base) // reference base
{ // RETURN clock - base
double clk;
#ifdef __linux__
static long tv_sec = 0;
static long tv_usec = 0;
struct timeval tval;
gettimeofday(&tval, NULL);
clk = tval.tv_sec + tval.tv_usec/1000000.0;
#else
clk = clock()/(double)CLOCKS_PER_SEC;
#endif
return clk - base;
}
/*--------------------------------------------------------------------------*/
#define a (*pa)
static void prnPtl( FILE *prn, char *txt, int np, PTLREC *pa, double tp )
{
int i;
fprintf(prn, "%s %4d: t=%1.0lf flg=%02x rfl=%02x src=%ld off=%d rnd=%ld tp=%1.1lf\n", //-2021-08-05
txt, np, a.t, a.flag, a.refl, a.srcind, a.offcmp, a.rnd, tp );
fprintf(prn, " x=(%11.3lf,%11.3lf,%11.3lf/%11.3lf) v=(%7.3f,%7.3f,%7.3f) \n",
a.x, a.y, a.z, a.h, a.u, a.v, a.w );
fprintf(prn, " afu=(%10.2e,%10.2e) \n", a.afuhgt, a.afutsc );
fprintf(prn, " g=(");
for (i=0; i fprintf(prn, " )\n");
}
static void trcPtl( int np, PTLREC *pa, double tp, MD3REC m )
{
int i;
fprintf(TrcFile,
"# %3d %ld %d %02x %02x: %8.2lf | %8.2lf %8.2lf %7.3lf | "
"%6.3f %6.3f %6.3f | %6.3f %6.3f %6.3f |",
np, a.srcind, a.offcmp, a.flag, a.refl, tp, a.x, a.y, a.z,
m.vx, m.vy, m.vz, a.u, a.v, a.w );
for (i=0; i fprintf(TrcFile, "\n");
}
#undef a
static void trcInit( int gp, int nl, int ni, long t1, long t2 )
{
fprintf( TrcFile,
"#### %d %d %d [%ld,%ld]\n", gp, nl, ni, t1, t2 );
}
/*===================================================================== ClcAfu
*/
static int ClcAfu( /* Berechnung der Abgasfahnen-Ueberhoeung */
SRCREC *ps, /* Pointer auf die Quelle */
float v2, /* Geschwindigkeitsquadrat */
float stab, /* Stabilitaetsklasse oder Temperatur-Gradient */
float *pu, /* Ergebnis: Ueberhoeung. */
float *pt ) /* Ergebnis: Anstiegszeit. */
{
dP(ClcAfu);
//float uq, dh, xe;
//int m, rc;
if (!ps) eX(1);
if (ps->ts >= 0) {
*pt = ps->ts;
*pu = ps->ts*ps->vq;
}
else if (ps->hrise >= 0) { //-2018-10-04
*pt = ps->trise;
*pu = ps->hrise;
}
/*
else {
m = PLM_VDI3782_3;
uq = sqrt(v2);
rc = TalPlm(m, &dh, &xe, ps->qq, ps->h, uq, stab,
ps->vq, ps->dq, ps->rhy, ps->lwc, ps->tmp); //-2002-12-10
if (rc < 0) {
if (rc == -1 && (PlmMsg)) eX(3); //-2002-09-23
eX(2);
}
*pu = dh;
if (uq < 1.0) uq = 1.0;
if (ps->qq>0 || ps->tmp>10) *pt = 0.4*xe/uq; //-2002-12-10
else if (ps->vq > 0) *pt = dh/ps->vq;
else *pt = 0;
}
*/
return 0;
eX_1:
eMSG(_no_source_data_);
/*
eX_3:
vMsg("%s", PlmMsg);
PlmMsg = NULL;
eX_2:
eMSG(_cant_calculate_plume_rise_);
* */
}
/*================================================================= ClcAfuRxyz
*/
static int ClcAfuRxyz( SRCREC *ps, float *pvs, float *prx, float *pry, float *prz )
{
dP(ClcAfuRxyz);
float sif, cof, sig, cog, dq;
float rx, ry, rz, ax, ay, az, r, vs, fq, gq, sh, sv, sl;
*prx = 0;
*pry = 0;
*prz = 1;
if (!ps) eX(1);
if (ps->tmp>0 && !(Flags & FLG_PLURIS)) //-2018-10-04
return 0;
fq = ps->fq;
sh = ps->sh;
sv = ps->sv;
sl = ps->sl;
if (fq==0 && sh<=0 && sv<=0 && sl<=0)
return 0;
vs = *pvs;
dq = ps->wq;
cof = cos(fq*PI/180);
sif = sin(fq*PI/180);
if (ps->gq > -999) {
gq = ps->gq; // gq: meterological orientation
sig = cos(gq*PI/180); // sig: mathematical orientation
cog = sin(gq*PI/180);
}
else {
gq = (ps->aq >= 0) ? 90-dq : dq;
cog = ps->co;
sig = ps->si;
}
rz = vs*cof; //-99-08-15
ry = vs*sif*sig;
rx = vs*sif*cog;
if (sh > 0) { //-99-08-05
ax = sig;
ay = -cog;
r = sh*WrkRndNrm();
rx += r*ax;
ry += r*ay;
}
if (sv > 0) {
ax = -cof*cog;
ay = -cof*sig;
az = sif;
r = sv*WrkRndNrm();
rx += r*ax;
ry += r*ay;
rz += r*az;
}
if (sl > 0) {
ax = sif*cog;
ay = sif*sig;
az = cof;
r = sl*WrkRndNrm();
rx += r*ax;
ry += r*ay;
rz += r*az;
}
r = sqrt(rx*rx + ry*ry + rz*rz);
*pvs = r;
if (r > 0) {
*prx = rx/r;
*pry = ry/r;
*prz = rz/r;
}
return 1;
eX_1:
eMSG(_no_source_data_);
}
/*======================================================================= ClcPlr
*/ //-2018-10-04
static int ClcPlr( // Call PLURIS to get the plume rise for this source
SRCREC *ps) // the source
{
dP(ClcPlr);
PLRSRC src;
PLRPRM prm;
PLRRSL rsl;
PLRPVL *ppvl=NULL;
PR3REC *p0, *p1;
ARYDSC *pprf=NULL;
char name[512];
//
int se = 100000;
int sn = 1000;
float sd = 0.01;
int adjust_step = 200;
float min_us = 0.05;
//
int k, k0, k1, nh, rc;
float h, vx, vy, v, d, tmp;
double *hhprf = NULL; // layer height at source position
long iprf=0, last_valid;
double c_0, c_1, c_2;
//
// explicit plume rise, nothing to do
if (ps->ts >= 0)
return 0;
//
// last calculated plume rise still valid, nothing to do
if (ps->valid >= fmin(SrcT2, ModT2)) //-2018-10-04
return 1;
//
// initialize to zero plume rise
ps->trise = 0;
ps->hrise = 0;
ps->sh = 0;
ps->sl = 0;
ps->sv = 0;
rsl.vr = 0; //-2024-01-17
rsl.tr = 0; //-2024-01-17
rsl.dr = 0; //-2024-01-17
//
// no valid meteorology
if (!Valid)
return 0;
//
// no volume flow, zero plume rise
last_valid = ps->valid;
if (last_valid < 0) last_valid = 0;
if (ps->vq <= 0) {
ps->valid = fmin(SrcT2, ModT2);
if (LogPluris) {
char t1s[40], t2s[40];
strcpy(t1s, TmString(&last_valid));
strcpy(t2s, TmString(&(ps->valid)));
vLOG(1)("WRK: [%s,%s] pluris('%s') = %s", t1s, t2s, ps->name, str_rsl(rsl));
}
return 0;
}
//
// prepare plume rise calculation with PLURIS
//
// check source parameters
if (ps->dq <= 0) eX(10);
if (ps->vq <= 0) eX(11);
if (ps->fq != 0) eX(12);
if (ps->rt < 0 || ps->rt > 1) eX(14);
if (ps->lwc < 0 || ps->lwc > 1) eX(21);
if (ps->vwc < 0 || ps->vwc > 1) eX(24);
if (ps->rhy < 0 || ps->rhy > 100) eX(22);
if (ps->rhy > 0 && ps->rhy <= 1) eX(15);
if (Rh < 0 || Rh > 100) eX(23);
if (Rh > 0 && Rh <= 1.0) eX(16);
if (Ta < 0 || Ta > 50) eX(17); //-2024-01-17
if (ps->rf < 0 || ps->rf > 1) eX(25); //-2024-01-17
//
c_0 = WrkClock(0);
vLOG(5)("WRK: ClcPlr('%s')", ps->name);
//
// profile heights Sk[] above ground (possibly reduced)
if (Gl > 1) { // nested grids
NETREC *pn = (NETREC*)GrdPtab->start; // outermost grid
iprf = IDENT(PRFarr, 0, (int)(pn->level), (int)(pn->index));
}
else { // single grid
iprf = IDENT(PRFarr, 0, Gl, Gi);
}
pprf = TmnAttach(iprf, NULL, NULL, 0, NULL); eG(106);
if (!pprf) eX(107);
k = (bDdmet) ? 2 : 0;
nh = pprf->bound[k].hgh - pprf->bound[k].low + 1;
TmnDetach(iprf, NULL, NULL, 0, NULL); eG(108);
pprf = NULL;
ppvl = (PLRPVL*)ALLOC(nh*sizeof(PLRPVL)); if (!ppvl) eX(3);
//
// source parameters
tmp = ps->tmp;
if (tmp < Ta)
tmp = Ta;
src.dm = ps->dq;
src.hb = ps->h;
src.zq = ps->wld;
src.lq = ps->lwc;
src.sq = ps->vwc;
src.rq = ps->rhy * 0.01;
src.te = tmp;
src.ve = ps->vq;
src.us = Us;
src.lm = Lm;
src.ta = Ta;
src.rf = ps->rf; //-2024-01-17
//
// calculation parameters
prm.sc = src.dm; // scale = source diameter
prm.sd = sd; // scaled step width
prm.se = se; // scaled maximum distance (ss)
prm.sn = sn; // number of log output lines
prm.adjust_step = adjust_step; // >0: dynamic step rr/adjust_step
prm.min_us = min_us; // minimum u* for break condition
prm.const_ust = 0;
prm.bf = PLRBREAKFACTOR;
if (MI.flags & FLG_BESMAX) { //-2024-01-17
prm.const_ust = 1;
prm.bf = PLRBREAKFACTOR_BESMAX;
}
vLOG(4)("WRK: ClcPlr, sc=%1.3f, sd=%1.4f, se=%1.0f, sn=%d, adjust_step=%d, min_us=%1.1f, bf=%1.2f, const_ust=%d\n",
prm.sc, prm.sd, prm.se, prm.sn, prm.adjust_step, prm.min_us, prm.bf, prm.const_ust);
//
// insert profile values
iprf = 0;
//
// 1D meteorology
if (!bDdmet) {
float a=0;
int nk;
//
if (Gl > 1) { // nested grids
NETREC *pn = (NETREC*)GrdPtab->start; // outermost grid
iprf = IDENT(PRFarr, 0, (int)(pn->level), (int)(pn->index));
}
else { // current (single) grid
iprf = IDENT(PRFarr, 0, Gl, Gi);
}
pprf = TmnAttach(iprf, NULL, NULL, 0, NULL); eG(6);
if (!pprf) eX(7);
nk = pprf->bound[0].hgh - pprf->bound[0].low + 1;
for (k=0; k a = 0;
//
// for 1D meteorology, use vertical grid points Sk
// for defining the vertical profiles for PLURIS
h = Sk[k];
//
// find the vertical interval [k0, k1] containing h
for (k1=1; k1 p1 = (PR3REC*)AryPtr(pprf, k1); if (!p1) eX(100);
if (p1->h >= h)
break;
}
if (k1 >= nk) eX(5);
k0 = k1 - 1;
p0 = (PR3REC*)AryPtr(pprf, k0); if (!p0) eX(101);
a = (h - p0->h)/(p1->h - p0->h);
vx = p0->vx + a*(p1->vx - p0->vx);
vy = p0->vy + a*(p1->vy - p0->vy);
v = sqrt(vx*vx + vy*vy);
d = (v == 0) ? 0 : atan2(-vx, -vy); // clockwise from north
d *= 180/PI;
if (d < 0) d += 360.;
ppvl[k].ha = h;
ppvl[k].ua = v;
ppvl[k].da = d;
}
TmnDetach(iprf, NULL, NULL, 0, NULL); eG(102);
pprf = NULL;
}
//
// 3D meteorology
else {
int in, nesting;
int ni, nj, nk; // number of points
double xmin, xmax, ymin, ymax, dd;
double xq, yq, ax, ay;
double zbot, ztop;
int iq, jq; // point index
int knext; // next valid k-value of ppvl
double hnext, hk; // next required height of PVL-array
int grid_found = 0;
NETREC *pn;
PR3REC *p000, *p100, *p010, *p110;
PR3REC *p001, *p101, *p011, *p111;
//
xq = ps->x;
yq = ps->y;
nesting = (GrdPtab != NULL);
hnext = 0;
//
in = GrdPprm->numnet;
knext = 0;
while(1) { // step through all grids, starting with the finest grid
in--;
vLOG(5)("WRK: ClcPlr, grid index = %d", in);
if (nesting) {
if (in < 0) {
nh = knext;
vLOG(5)("WRK: ClcPlr, grid restricts profiles to height %1.1f", Sk[nh-1]);
break;
}
pn = ((NETREC*)GrdPtab->start) + in;
xmin = pn->xmin;
xmax = pn->xmax;
ymin = pn->ymin;
ymax = pn->ymax;
if (xq < xmin || xq > xmax || yq < ymin || yq > ymax)
continue; // try the next coarser grid
grid_found = 1;
dd = pn->delta;
iprf = IDENT(PRFarr, 0, (int)(pn->level), (int)(pn->index));
}
else { // no nesting, use the current (single) grid
if (in < -1) {
nh = knext;
vLOG(5)("WRK: ClcPlr, grid restricts profiles to height %1.1f", Sk[nh-1]);
break;
}
xmin = GrdPprm->xmin;
xmax = GrdPprm->xmax;
ymin = GrdPprm->ymin;
ymax = GrdPprm->ymax;
if (xq < xmin || xq > xmax || yq < ymin || yq > ymax)
break;
grid_found = 1;
dd = GrdPprm->dd;
iprf = IDENT(PRFarr, 0, 0, 0);
}
strcpy(name, NmsName(iprf));
vLOG(5)("WRK: ClcPlr, using profile %s", name);
pprf = TmnAttach(iprf, NULL, NULL, 0, NULL); eG(103);
if (pprf == NULL) eX(104);
ni = pprf->bound[0].hgh - pprf->bound[0].low + 1;
nj = pprf->bound[1].hgh - pprf->bound[1].low + 1;
nk = pprf->bound[2].hgh - pprf->bound[2].low + 1;
// find the grid cell containing the source position
iq = (int)floor((xq - xmin)/dd);
if (iq == ni-1)
iq = ni-2;
jq = (int)floor((yq - ymin)/dd);
if (jq == nj-1)
jq = nj-2;
// just in case, should not happen
if (iq < 0 || iq > ni-1 || jq < 0 || jq > nj-1) eX(201); //-2022-10-04
ax = (xq - xmin - iq*dd)/dd;
ay = (yq - ymin - jq*dd)/dd;
p000 = AryPtr(pprf, iq, jq, 0);
p100 = AryPtr(pprf, iq+1, jq, 0);
p010 = AryPtr(pprf, iq, jq+1, 0);
p110 = AryPtr(pprf, iq+1, jq+1, 0);
zbot = (1-ax)*(1-ay)*p000->h + ax*(1-ay)*p100->h + (1-ax)*ay*p010->h
+ ax*ay*p110->h;
p001 = AryPtr(pprf, iq, jq, nk-1);
p101 = AryPtr(pprf, iq+1, jq, nk-1);
p011 = AryPtr(pprf, iq, jq+1, nk-1);
p111 = AryPtr(pprf, iq+1, jq+1, nk-1);
ztop = (1-ax)*(1-ay)*p001->h + ax*(1-ay)*p101->h + (1-ax)*ay*p011->h
+ ax*ay*p111->h;
vLOG(5)("WRK: ClcPlr, zbot=%e, ztop=%e, hnext=%e", zbot, ztop, hnext);
if (hnext > ztop-zbot) { // profile not suited
if (nesting) {
TmnDetach(iprf, NULL, NULL, 0, NULL); eG(105);
pprf = NULL;
continue;
}
else
break;
}
//
// Calculate the height above ground of the profile layers at the
// source position.
//
hhprf = ALLOC(nk*sizeof(double));
hhprf[0] = 0;
hhprf[nk-1] = ztop - zbot;
for (k=1; k p000 = AryPtr(pprf, iq, jq, k);
p100 = AryPtr(pprf, iq+1, jq, k);
p010 = AryPtr(pprf, iq, jq+1, k);
p110 = AryPtr(pprf, iq+1, jq+1, k);
hhprf[k] = (1-ax)*(1-ay)*p000->h + ax*(1-ay)*p100->h + (1-ax)*ay*p010->h
+ ax*ay*p110->h - zbot;
}
//
// Interpolate profile values at the heights of the PVL array.
//
for (k=knext; k //
// for 3D meteorology, use the mid points of the vertical
// grid points Sk for defining the vertical profiles for PLURIS;
// this is consistent with the definition of the Arakawa-C wind
// components as area averages.
hk = (k == 0) ? 0 : 0.5*(Sk[k]+Sk[k-1]); //-2018-10-04
if (hk == 0) {
k1 = 1;
}
else {
for(k1=1; k1 if (hhprf[k1] >= hk)
break;
}
if (k1 >= nk) {
knext = k;
hnext = hk;
break;
}
k0 = k1 - 1;
p110 = AryPtr(pprf, iq+1, jq+1, k0);
if (p110->vz <= WND_VOLOUT) {
v = -1;
d = 0;
}
else {
p101 = AryPtr(pprf, iq+1, jq, k1);
p011 = AryPtr(pprf, iq, jq+1, k1);
p111 = AryPtr(pprf, iq+1, jq+1, k1);
vx = (1-ax)*p011->vx + ax*p111->vx; //-2018-10-04
vy = (1-ay)*p101->vy + ay*p111->vy; //-2018-10-04
v = sqrt(vx*vx + vy*vy);
d = (v == 0) ? 0 : atan2(-vx, -vy); // clockwise from north
d *= 180/PI;
if (d < 0) d += 360.;
}
ppvl[k].ha = hk;
ppvl[k].ua = v;
ppvl[k].da = d;
} // for k
FREE(hhprf);
TmnDetach(iprf, NULL, NULL, 0, NULL); eG(102);
pprf = NULL;
if (k >= nh)
break;
} // while
if (!grid_found) eX(200); //-2022-10-04
for (k=nh-1; k>=0; k--) { // insert wind speed in region with body cells
if (ppvl[k].ua >= 0)
continue;
ppvl[k].ua = ppvl[k+1].ua;
ppvl[k].da = ppvl[k+1].da;
}
} // 3d meteorology
//
// set meteorological parameters beside wind
for (k=0; k ppvl[k].ta = BlmTemp(ppvl[k].ha, Z0, D0, Us, Lm, Ta, 0);
ppvl[k].ca = 1;
ppvl[k].ra = Rh*0.01;
}
//
// set wind direction at the ground if required
if (ppvl[0].ua == 0 && ppvl[0].ra == 0) {
ppvl[0].da = ppvl[1].da;
}
//
// run pluris
PlrVerbose = StdLogLevel - 3 + 3*LogPluris;
PlrMessage = MsgFile;
c_1 = WrkClock(c_0);
rc = run_pluris(&src, nh, ppvl, &prm, &rsl); if (rc < 0) eX(4); //-2024-01-17
c_2 = WrkClock(c_0);
//
// rise parameters
if (MI.flags & FLG_BESMAX) { //-2024-01-17
rsl.tr = 0.0;
}
ps->hrise = rsl.dr;
ps->trise = rsl.tr;
//
// set source induced turbulence
ps->sh = rsl.vr * ps->rt;
ps->sl = ps->sh;
ps->sv = ps->sh;
//
if (MI.flags & FLG_BESMAX) { //-2024-01-17
ps->sh = 0.;
ps->sl = 0.;
ps->sv = 0.;
}
ps->valid = fmin(SrcT2, ModT2); //-2018-10-04
if (LogPluris) {
char t1s[40], t2s[40];
strcpy(t1s, TmString(&last_valid));
strcpy(t2s, TmString(&(ps->valid)));
vLOG(1)("WRK: [%s,%s] pluris('%s') = %s Set[rt=%5.3f, sh|sl|sv=%7.3f] Time[%1.4f]",
t1s, t2s, ps->name, str_rsl(rsl), ps->rt, ps->sh, c_2-c_1);
}
FREE(ppvl);
ppvl = NULL;
return 0;
eX_3:
eMSG("can't create metorological profile for source '%s'!", ps->name);
eX_4:
eMSG("PLURIS failed for source '%s' (%d)!", ps->name, rc);
eX_5:
eMSG("vertical extension of meteorological profile '%s' not sufficient!", name);
eX_106:
eMSG("can't get outer meteorological profile '%s'", name);
eX_6:
eMSG("can't get outer meteorological profile '%s'", name);
eX_107:
eMSG("missing outer meteorological profile '%s'", name);
eX_7:
eMSG("missing outer meteorological profile '%s'", name);
eX_10:
eMSG("plume rise calculation requires diameter > 0 for source '%s'!", ps->name);
eX_11:
eMSG("plume rise calculation requires exit velocity > 0 for source '%s'!", ps->name);
eX_12:
eMSG("Fq must be 0 (only plume rise for vertical exit implemented) for source '%s'!", ps->name);
eX_14:
eMSG("Rt (%1.1f) not in range [0,1] for source '%s'!", ps->rt, ps->name);
eX_15:
eMSG("Rh must be specified in percent (0 or in range [1,100]) for source '%s'!", ps->name);
eX_16:
eMSG("Rh of ambient air must be specified in percent (0 or in range [1,100])!");
eX_17:
eMSG("Ta of ambient air must be in range [0,50])!");
eX_21:
eMSG("Lw (%1.3f) not in range [0,1] for source '%s'!", ps->lwc, ps->name);
eX_22:
eMSG("Rh (%1.1f) not in range [0,100] for source '%s'!", ps->rhy, ps->name);
eX_23:
eMSG("Ambient Rh (%1.1f) not in range [0,100]!", Rh);
eX_24:
eMSG("Vw (%1.3f) not in range [0,1] for source '%s'!", ps->vwc, ps->name);
eX_25:
eMSG("Rf (%1.3f) not in range [0,1] for source '%s'!", ps->rf, ps->name);
eX_200:
eMSG("Source '%s' outside grid(s)!", ps->name);
eX_100: eX_101: eX_102: eX_103: eX_104: eX_105: eX_108: eX_201:
eMSG("internal error");
}
//===================================================================ClcPtlLoc
//
static int ClcPtlLoc( /* Teilchen im Netz einordnen, evtl. reflektieren */
PTLREC *pa, /* Pointer auf alte Teilchen-Koordinaten */
PTLREC *pp, /* Pointer auf neue Teilchen-Koordinaten */
LOCREC *pl, /* Pointer auf Index-Werte */
int h2z ) /* h-Wert nach z-Wert umwandeln */
{
// dP(ClcPtlLoc);
int i, j, k, flg, m;
float r, x, y, z, dz, xi, eta, zeta;
float zg;
flg = GRD_REFBOT;
z = (pa) ? pa->z : 0; //-2002-04-10
if ((pa) && (Zref) && (z<=Zref)) flg |= GRD_REFTOP; //-2002-04-10
zeta = 0;
if (h2z == WRK_USE_H) { z = pp->h; flg |= GRD_PARM_IS_H; }
else {
z = pp->z;
if (h2z == WRK_USE_AZ) flg |= GRD_USE_ZETA;
zeta = pl->az;
}
x = pp->x;
i = -1;
if (x < Xmin) {
if (bPerx) {
x = Xmax - fmod(Xmin-x, Gx); //-2001-03-23 lj
pp->x = x;
}
else i = 0;
}
else if (x > Xmax) {
if (bPerx) {
x = Xmin + fmod(x-Xmax, Gx); //-2001-03-23 lj
pp->x = x;
}
else i = 0;
}
y = pp->y;
j = -1;
if (y < Ymin) {
if (bPery) {
y = Ymax - fmod(Ymin-y, Gy); //-2001-03-23 lj
pp->y = y;
}
else j = 0;
}
else if (y > Ymax) {
if (bPery) {
y = Ymin + fmod(y-Ymax, Gy); //-2001-03-23 lj
pp->y = y;
}
else j = 0;
}
xi = (x - Xmin)/Dd;
eta = (y - Ymin)/Dd;
pl->i = 0;
pl->j = 0;
k = pl->k;
pl->k = 0;
if (bDdmet) {
if (i==0 || j==0) { pp->flag |= SRC_GOUT; return 0; }
m = GrdLocate(Pma, xi, eta, &z, flg, Zref, &i, &j, &k, &zg, &zeta);
if (m & GRD_REFTOP) pp->refl |= WRK_REFT;
if (m & GRD_REFBOT) pp->refl |= WRK_REFB;
if (!(m & GRD_INSIDE)) pp->flag |= SRC_GOUT;
}
else {
zg = 0;
if (i==0 || j==0) {
if (x<=Xbmin || x>=Xbmax) pp->flag |= SRC_XOUT;
if (y<=Ybmin || y>=Ybmax) pp->flag |= SRC_YOUT;
}
else {
i = 1 + xi; if (i > Nx) i = Nx;
j = 1 + eta; if (j > Ny) j = Ny;
}
if (z < 0) { z = -z; pp->refl |= WRK_REFB; }
if (flg&GRD_REFTOP && z>Zref) { z = 2*Zref-z; pp->refl |= WRK_REFT; }
for (k=1; k<=Nz; k++)
if (z <= Sk[k]) { //-2002-04-10
dz = Sk[k] - Sk[k-1];
zeta = (z - Sk[k-1])/dz;
break;
}
if (k > Nz) {
k = 0; pp->flag |= SRC_ZOUT; }
}
if (i < 0) i = 0;
if (j < 0) j = 0;
pl->i = i;
pl->j = j;
pl->k = k;
pp->z = z;
if ((i) && (j)) pp->h = z - zg;
if (pp->flag & SRC_GOUT) return 0;
pl->ax = xi - (i-1);
pl->ay = eta - (j-1);
pl->az = zeta;
if (bDdmet) {
r = WrkRnd01();
if (r > pl->ax) pl->is = i-1;
else pl->is = i;
r = WrkRnd01();
if (r > pl->ay) pl->js = j-1;
else pl->js = j;
pl->ks = k;
}
else {
pl->is = i;
pl->js = j;
pl->ks = k;
}
return 1;
}
//================================================================== ClcLocMod
//
static int ClcLocMod( /* Feld am Ort des Teilchens berechnen */
LOCREC *pl, /* Pointer auf Index-Werte */
MD3REC *p3, /* Pointer auf Feld-Werte */
int clceta ) /* <=0, wenn Eta nicht interpoliert werden soll */
{
dP(ClcLocMod);
MD2REC *p2a, *p2b;
MD3REC *p3a, *p3b;
MD3REC *p000, *p001, *p010, *p100, *p011, *p101, *p110, *p111;
float a000, a001, a010, a100, a011, a101, a110, a111;
float ax, bx, ay, by, az, bz;
int i, j, k, is, js;
float a, b;
i = pl->i;
j = pl->j;
k = pl->k;
is = pl->is;
js = pl->js;
if (k == 0) eX(1);
if (bDdmet) {
if (i==0 || j==0) eX(2);
p3a = (MD3REC*) AryPtr(Pma, i, j, k-1); if (!p3a) eX(4);
if (p3a->vz <= WND_VOLOUT) return 0; /* Innerhalb eines Koerpers */
p3b = (MD3REC*) AryPtr(Pma, i, j, k ); if (!p3b) eX(5);
if (p3b->ths <= WND_VOLOUT) return 0; /* Innerhalb eines Koerpers */
if (p3b->wz <= WND_VOLOUT) return 0; /* Innerhalb eines Koerpers */
memset(p3, 0, sizeof(MD3REC)); /* Zuerst alles auf 0 setzen */
if (!OnlyAdvec) {
p3->ths = p3b->ths; /* Gradient der potentiellen Temperatur */
p3->wx = p3b->wx; /* Kompensation, x-Komponente */
p3->wy = p3b->wy; /* Kompensation, y-Komponente */
p3->wz = p3b->wz; /* Kompensation, z-Komponente */
}
if (DefMode & GRD_3LIN) { /* Geschwindigkeiten 3*linear */
p000 = (MD3REC*) AryPtr(Pma, i-1, j-1, k-1); if (!p000) eX(10);
p001 = (MD3REC*) AryPtr(Pma, i-1, j-1, k ); if (!p001) eX(11);
p010 = (MD3REC*) AryPtr(Pma, i-1, j , k-1); if (!p010) eX(12);
p100 = (MD3REC*) AryPtr(Pma, i , j-1, k-1); if (!p100) eX(13);
p011 = (MD3REC*) AryPtr(Pma, i-1, j , k ); if (!p011) eX(14);
p101 = (MD3REC*) AryPtr(Pma, i , j-1, k ); if (!p101) eX(15);
p110 = (MD3REC*) AryPtr(Pma, i , j , k-1); if (!p110) eX(16);
p111 = (MD3REC*) AryPtr(Pma, i , j , k ); if (!p111) eX(17);
bx = pl->ax;
ax = 1.0 - bx;
by = pl->ay;
ay = 1.0 - by;
bz = pl->az;
az = 1.0 - bz;
a000 = ax*ay*az;
a001 = ax*ay*bz;
a010 = ax*by*az;
a100 = bx*ay*az;
a011 = ax*by*bz;
a101 = bx*ay*bz;
a110 = bx*by*az;
a111 = bx*by*bz;
p3->vx = a000*p000->vx + a001*p001->vx + a010*p010->vx + a100*p100->vx +
a011*p011->vx + a101*p101->vx + a110*p110->vx + a111*p111->vx;
p3->vy = a000*p000->vy + a001*p001->vy + a010*p010->vy + a100*p100->vy +
a011*p011->vy + a101*p101->vy + a110*p110->vy + a111*p111->vy;
p3->vz = a000*p000->vz + a001*p001->vz + a010*p010->vz + a100*p100->vz +
a011*p011->vz + a101*p101->vz + a110*p110->vz + a111*p111->vz;
}
else {
p3a = (MD3REC*) AryPtr(Pma, i-1, j, k); if (!p3a) eX(20);
p3->vx = (1-pl->ax)*p3a->vx + pl->ax*p3b->vx; /* Vx interpolieren */
p3a = (MD3REC*) AryPtr(Pma, i, j-1, k); if (!p3a) eX(21);
p3->vy = (1-pl->ay)*p3a->vy + pl->ay*p3b->vy; /* Vy interpolieren */
p3a = (MD3REC*) AryPtr(Pma, i, j, k-1); if (!p3a) eX(22);
p3->vz = (1-pl->az)*p3a->vz + pl->az*p3b->vz; /* Vz interpolieren */
}
if (PrfMode == 3) { /* Anisotrope Modellierung */
a = 1 - pl->az;
b = pl->az;
p3a = (MD3REC*) AryPtr(Pma, is, js, k-1); if (!p3a) eX(23);
p3b = (MD3REC*) AryPtr(Pma, is, js, k ); if (!p3b) eX(24);
p3->tau = a*p3a->tau + b*p3b->tau;
if (OnlyAdvec) return 1;
if (clceta < 0) return 1;
p3->pxx = a*p3a->pxx + b*p3b->pxx;
p3->pxy = a*p3a->pxy + b*p3b->pxy;
p3->pxz = a*p3a->pxz + b*p3b->pxz;
p3->pyx = a*p3a->pyx + b*p3b->pyx;
p3->pyy = a*p3a->pyy + b*p3b->pyy;
p3->pyz = a*p3a->pyz + b*p3b->pyz;
p3->pzx = a*p3a->pzx + b*p3b->pzx;
p3->pzy = a*p3a->pzy + b*p3b->pzy;
p3->pzz = a*p3a->pzz + b*p3b->pzz;
p3->lxx = a*p3a->lxx + b*p3b->lxx;
p3->lyx = a*p3a->lyx + b*p3b->lyx;
p3->lyy = a*p3a->lyy + b*p3b->lyy;
p3->lzx = a*p3a->lzx + b*p3b->lzx;
p3->lzy = a*p3a->lzy + b*p3b->lzy;
p3->lzz = a*p3a->lzz + b*p3b->lzz;
if (clceta == 0) return 1;
p3->exx = a*p3a->exx + b*p3b->exx;
p3->eyx = a*p3a->eyx + b*p3b->eyx;
p3->eyy = a*p3a->eyy + b*p3b->eyy;
p3->ezx = a*p3a->ezx + b*p3b->ezx;
p3->ezy = a*p3a->ezy + b*p3b->ezy;
p3->ezz = a*p3a->ezz + b*p3b->ezz;
return 1;
}
else { /* Isotrope Modellierung */
a = 1 - pl->az;
b = pl->az;
p2a = (MD2REC*) AryPtr(Pma, is, js, k-1); if (!p2a) eX(25);
p2b = (MD2REC*) AryPtr(Pma, is, js, k ); if (!p2b) eX(26);
p3->tau = a*p2a->tau + b*p2b->tau;
if (OnlyAdvec) return 1;
if (clceta < 0) return 1;
p3->pxx = a*p2a->pvv + b*p2b->pvv;
p3->pyy = p3->pxx;
p3->pzz = a*p2a->pww + b*p2b->pww;
p3->lxx = a*p2a->lvv + b*p2b->lvv;
p3->lyy = p3->lxx;
p3->lzz = a*p2a->lww + b*p2b->lww;
if (clceta == 0) return 1;
p3->exx = a*p2a->evv + b*p2b->evv;
p3->eyy = p3->exx;
p3->ezz = a*p2a->eww + b*p2b->eww;
return 1;
}
}
else { /* 1-dimensionale Meteorologie */
if ((k <= Nz) && (k > 0)) {
memset(p3, 0, sizeof(MD3REC)); /* Zuerst alles auf 0 setzen */
a = 1 - pl->az;
b = pl->az;
if (PrfMode == 3) { /* Anisotrope Modellierung */
p3b = (MD3REC*) AryPtr(Pma, k ); if (!p3b) eX(27);
p3a = (MD3REC*) AryPtr(Pma, k-1); if (!p3a) eX(28);
p3->vx = a*p3a->vx + b*p3b->vx;
p3->vy = a*p3a->vy + b*p3b->vy;
p3->vz = a*p3a->vz + b*p3b->vz;
p3->tau = a*p3a->tau + b*p3b->tau;
if (OnlyAdvec) return 1;
p3->wz = p3b->wz;
p3->ths = p3b->ths;
if (clceta < 0) return 1;
p3->pxx = a*p3a->pxx + b*p3b->pxx;
p3->pxy = a*p3a->pxy + b*p3b->pxy;
p3->pxz = a*p3a->pxz + b*p3b->pxz;
p3->pyx = a*p3a->pyx + b*p3b->pyx;
p3->pyy = a*p3a->pyy + b*p3b->pyy;
p3->pyz = a*p3a->pyz + b*p3b->pyz;
p3->pzx = a*p3a->pzx + b*p3b->pzx;
p3->pzy = a*p3a->pzy + b*p3b->pzy;
p3->pzz = a*p3a->pzz + b*p3b->pzz;
p3->lxx = a*p3a->lxx + b*p3b->lxx;
p3->lyx = a*p3a->lyx + b*p3b->lyx;
p3->lyy = a*p3a->lyy + b*p3b->lyy;
p3->lzx = a*p3a->lzx + b*p3b->lzx;
p3->lzy = a*p3a->lzy + b*p3b->lzy;
p3->lzz = a*p3a->lzz + b*p3b->lzz;
if (clceta == 0) return 1;
p3->exx = a*p3a->exx + b*p3b->exx;
p3->eyx = a*p3a->eyx + b*p3b->eyx;
p3->eyy = a*p3a->eyy + b*p3b->eyy;
p3->ezx = a*p3a->ezx + b*p3b->ezx;
p3->ezy = a*p3a->ezy + b*p3b->ezy;
p3->ezz = a*p3a->ezz + b*p3b->ezz;
return 1; }
else { /* Isotrope Modellierung */
p2b = (MD2REC*) AryPtr(Pma, k ); if (!p2b) eX(29);
p2a = (MD2REC*) AryPtr(Pma, k-1); if (!p2a) eX(30);
p3->vx = a*p2a->vx + b*p2b->vx;
p3->vy = a*p2a->vy + b*p2b->vy;
p3->vz = a*p2a->vz + b*p2b->vz;
p3->tau = a*p2a->tau + b*p2b->tau;
if (OnlyAdvec) return 1;
p3->wz = p2b->wz;
p3->ths = p2b->ths;
if (clceta < 0) return 1;
p3->pxx = a*p2a->pvv + b*p2b->pvv;
p3->pyy = p3->pxx;
p3->pzz = a*p2a->pww + b*p2b->pww;
p3->lxx = a*p2a->lvv + b*p2b->lvv;
p3->lyy = p3->lxx;
p3->lzz = a*p2a->lww + b*p2b->lww;
if (clceta == 0) return 1;
p3->exx = a*p2a->evv + b*p2b->evv;
p3->eyy = p3->exx;
p3->ezz = a*p2a->eww + b*p2b->eww;
return 1; }
}
else eG(31); /* k > Nz */
}
return 1;
eX_1: eX_2: eX_4: eX_5:
eX_10: eX_11: eX_12: eX_13: eX_14: eX_15: eX_16: eX_17:
eX_20: eX_21: eX_22: eX_23: eX_24: eX_25: eX_26: eX_27: eX_28: eX_29:
eX_30: eX_31:
eMSG(_loc_$_, pl->i, pl->j, pl->k, pl->is, pl->js, pl->ax, pl->ay, pl->az);
}
//================================================================= Try2Escape
//
static int Try2Escape( /* In einem Gebaeude erzeugte Partikel herausdraengen */
PTLREC *pa,
LOCREC *pl ) {
dP(WRK:Try2Escape);
double x = pa->x; // Anfangs-Position merken
double y = pa->y;
double h = pa->h;
double d = 1;
int ii[9] = { 0, 1, -1, 0, 0, 1, 1, -1, -1 };
int jj[9] = { 0, 0, 0, 1, -1, 1, -1, 1, -1 };
int n, i, k, kk;
for (n=1; n<=10; n++) { // 10 Distanzen pruefen
d = n*0.1*Dd;
for (kk=0; kk<=2; kk++) {
k = (kk == 2) ? -1 : kk;
pa->h = h + k*d;
for (i=(k==0); i<9; i++) {
pa->x = x + ii[i]*d;
pa->y = y + jj[i]*d;
ClcPtlLoc(NULL, pa, pl, WRK_USE_H); eG(2); // neu einordnen
if (!INSOLID(pl->i,pl->j,pl->k)) { // geschafft!
d = 0;
break;
}
} // for i
if (d == 0) break;
}
if (d == 0) break;
}
if (d > 0) eX(3);
return 0;
eX_2: eMSG(_cant_locate_); //-2004-12-23
eX_3: eMSG(_cant_shift_); //-2004-12-23
}
//----------------------------------------------------------------------------
// utilities for rain drop drift -2023-07-17
//----------------------------------------------------------------------------
/**
* Register deposited masses
*/
static int WetDepoRegister(
float xdrop, // x-coordinate of wet deposition
float ydrop, // y-coordinate of wet deposition
float* wdrops) // deposition values [Nc]
{
dP(WetDepoRegister);
int inet, net, ic, ibuf, idrop, jdrop, irec, nd, nx, ny;
float *buffer_field = NULL;
float *record;
NETREC* pndst;
//
// deposition inside calculation area?
if (GrdPprm->numnet == 0)
return 0;
if (xdrop=WetXmax || ydrop=WetYmax)
return 0;
//
// identify target grid: only check grids not yet processed
for (net=GrdPprm->net-1; net>0; net--) {
pndst = (NETREC*) AryPtr(GrdPtab, net);
if (xdrop >= pndst->xmin && xdrop < pndst->xmax
&& ydrop >= pndst->ymin && ydropymax)
break;
}
inet = net - 1;
if (inet < 0)
return -1;
ibuf = Gp*numnet + inet;
nx = floor(pndst->nx + 0.5); //-2023-07-31
ny = floor(pndst->ny + 0.5); //-2023-07-31
if (WetDepoList[ibuf] == NULL) {
nd = nx*ny*Nc*sizeof(float);
WetDepoList[ibuf] = ALLOC(nd);
if (!WetDepoList[ibuf]) eX(1);
}
idrop = floor((xdrop - pndst->xmin)/pndst->delta);
jdrop = floor((ydrop - pndst->ymin)/pndst->delta);
if (idrop < 0) idrop = 0; //-2023-07-31
if (idrop > nx-1) idrop = nx-1; //-2023-07-31
if (jdrop < 0) jdrop = 0; //-2023-07-31
if (jdrop > ny-1) jdrop = ny-1; //-2023-07-31
buffer_field = WetDepoList[ibuf];
irec = Nc*(jdrop + ny*idrop);
record = buffer_field + irec;
for (ic=0; ic record[ic] += wdrops[ic];
return 0;
eX_1: eMSG("internal error");
}
/*
* Addition of deposition values stored in WetDepoList to the values of
* the current dose array DOSarr.
*/
static int WetDepoAdd() {
dP(WetDepoAdd);
int inet = Pgp->net - 1;
if (inet < 0)
return -1;
int ibuf = Gp*numnet + inet;
float *pwd = (float*)WetDepoList[ibuf];
if (pwd == NULL)
return 0;
for (int i=0; i for (int j=0; j float* padd = pwd + Nc*(j + Ny*i);
float* pdos = AryPtr(Pda, i+1, j+1, -1, 0);
for (int ic=0; ic pdos[ic] += padd[ic];
padd[ic] = 0;
}
}
}
return 0;
}
static int WetDepoPrint(int inet, int gp) {
int ibuf = gp*numnet + inet;
float *dd = (float*)WetDepoList[ibuf];
NETREC* pn = (NETREC*)AryPtr(GrdPtab, inet+1);
int nx = pn->nx;
int ny = pn->ny;
if (dd == NULL)
return -1;
for (int j=ny-1; j>=0; j--) {
printf("%3d: ", j+1);
for (int i=0; i float d = dd[Nc*(j + ny*i)];
printf("%8.2e ", d);
}
printf("\n");
}
return 0;
}
//----------------------------------------------------------------------------
// end utilities for rain drop drift
//----------------------------------------------------------------------------
//===================================================================== RunPtl
//
static int RunPtl(
void )
{
dP(RunPtl);
long rc;
static int totals = 0;
static int steps = 0;
PTLREC *pp, *pa, *pe;
#define a (*pa)
#define e (*pe)
LOCREC l, ll;
MD3REC m, mm;
PPMREC *pr;
SRCREC *ps;
float dg, gsum, r1, r2, r3, vsed, vdep, *pf, *pw, dafu;
float wdep, v2, tau, vred, hred, dzp, afufac;
float *wetbuf; //-2023-07-17
float vxdrop, vydrop; //-2023-07-17
float src_rf; //-2024-01-17
float wdeps[60];
double tp;
int handle, retry;
int icmp, jcmp, cmpnum, cmpoff, nzdos, np, kk, k, srcind;
int imin, imax, jmin, jmax; //-2023-07-17
int deposition, discard, created, wdep_unknown;
int washout; //-2011-11-21
int predict, predicting; //-2001-09-22
vLOG(5)("WRK:RunPtl() starting");
if (WetDrift) { //-2023-07-17
// transfer wet deposition that entered into the current grid by drop drift
// from buffer WetDepoList to this dose array DOSarr.
WetDepoAdd(); // applies Pgp, Pda, numnet
}
wetbuf = ALLOC(Nc*sizeof(float));
imin = Pda->bound[0].low;
imax = Pda->bound[0].hgh;
jmin = Pda->bound[1].low;
jmax = Pda->bound[1].hgh;
//
pa = ALLOC(PtlRecSize); if (!pa) eX(30);
pe = ALLOC(PtlRecSize); if (!pe) eX(30);
memset(&l, 0, sizeof(l)); /* Strukturen zu 0 initialisieren */
memset(&m, 0, sizeof(m));
memset(&e, 0, PtlRecSize);
nzdos = Pda->bound[2].hgh;
handle = PtlStart(Ppa); eG(1);
NumRun = 0;
predict = (NoPredict) ? 0 : (bDdmet != 0); //-2001-09-22
for (np=0; ; ) { /*+++++++++++++++++ Fuer alle Teilchen: ++++++++++++++++*/
pp = PtlNext(handle); /* Neues Teilchen anfordern */
if (pp == NULL) break; /* aufhoeren, wenn keins mehr da ist */
totals++;
if (pp->flag & SRC_REMOVE) /* naechstes nehmen, wenn aussortiert */
continue;
if (!Valid) { // Keine gueltige Meteorologie:
pp->flag |= SRC_REMOVE; // Alle vorhandenen Partikel werden
NumRemove++; // geloescht.
continue;
}
memcpy(&a, pp, PtlRecSize); /* umspeichern */
if (a.t >= T2) continue; /* naechstes nehmen, wenn abgearbeitet */
if (a.t < T1) eX(2);
tp = a.t; /* Zeit separat merken */
created = (a.flag & SRC_CREATED) != 0; /* neu erzeugt, falls Flag gesetzt */ //-2003-02-21
a.flag = SRC_EXIST;
WrkRndSetSeed(a.rnd); // Zufallszahlen initialisieren //-2003-02-21
cmpoff = a.offcmp; /* Komponenten-Offset im Dosis-Feld */
cmpnum = a.numcmp; /* Zahl der Komponenten */
srcind = a.srcind - 1; /* Index der zugehoerigen Quelle */
ps = (SRCREC*) AryPtr(Psa, srcind); // Pointer auf die zugehoerige Quelle
if (!ps) eX(13);
afufac = ps->fr; // Faktor bei Reflektion
src_rf = ps->rf; // fuer VDI 3673/1 Schwergas -2024-01-17
pr = (PPMREC*) Ppp->start + cmpoff;
deposition = 0;
washout = 0; //-2011-11-21
for (icmp=0; icmp if (pr[icmp].wdep > 0) deposition = 1;
if (pr[icmp].rwsh > 0) washout = 1; //-2011-11-21
}
if (Pwa) deposition = 1;
wdep_unknown = 0;
if (a.vg > 0) { /* Falls individuelles Vsed */
vsed = a.vg;
deposition = 1;
wdep_unknown = 1;
}
else vsed = pr->vsed; /* Sedimentations-Geschwindigkeit */
vred = pr->vred; /* Reduktion der vert. Startgeschwindigk.*/
hred = pr->hred; /* Reduktion der hor. Startgeschwindigk. */
np++; /* Teilchenzaehler hochzaehlen ... */
NumRun = np; /* ... und extern zur Verfuegung stellen */
ClcPtlLoc(NULL, &a, &l, WRK_USE_H); eG(3); /* Teilchen einordnen */
if (a.flag & SRC_EOUT) { /* Falls ausserhalb des Netzes, ... */
if (created) a.flag |= SRC_CREATED; //-2005-01-23
pp->flag = a.flag; /* ... vermerken und naechstes nehmen */
continue; }
if (INSOLID(l.i,l.j,l.k)) {
if (created) {
Try2Escape(&a, &l); eG(14); //-2004-10-10
}
else {
Try2Escape(&a, &l); eG(15); //-2004-12-09
}
}
kk = (a.z > Zref) ? -1 : created; /* Oberhalb Zref ? */
ClcLocMod(&l, &m, kk); eG(4); /* Felder interpolieren */
if (created) { /* Bei einem neuen Teilchen ... */
r1 = WrkRnd(); /* ... Geschwindigkeiten initalisieren */
r2 = WrkRnd();
r3 = WrkRnd();
vred = pr->vred; /* Reduktion der vert. Startgeschwindigk.*/
hred = pr->hred; /* Reduktion der hor. Startgeschwindigk. */
a.u = hred*m.exx*r1;
a.v = hred*(m.eyx*r1 + m.eyy*r2);
a.w = vred*(m.ezx*r1 + m.ezy*r2 + m.ezz*r3);
if (!ExactTau) m.tau *= 0.5 + WrkRnd01(); //-2001-09-23
v2 = m.vx*m.vx + m.vy*m.vy;
ClcAfu(ps, v2, Stab, &a.afuhgt, &a.afutsc); eG(5);
if (a.afutsc > 0) {
float vs;
vs = a.afuhgt/a.afutsc;
ClcAfuRxyz(ps, &vs, &a.afurx, &a.afury, &a.afurz);
a.afuhgt = vs*a.afutsc;
}
else {
a.afurx = 0; a.afury = 0; a.afurz = 1; }
}
retry = 0; /* Zaehler fuer Schrittversuche setzen */
memcpy(&e, &a, PtlRecSize); /* Alles uebernehmen (Massen g[] !) */
predicting = predict;
while (tp < T2) { /*---------- Zeitschritte durchfuehren -----------*/
tau = m.tau;
if (predicting) { //-2001-09-22
predicting = 0;
e.x = a.x + tau*m.vx; // Nur horizontale ...
e.y = a.y + tau*m.vy; // ... Advektion durchfuehren
e.z = a.z;
ll = l;
while (1) {
ClcPtlLoc(&a, &e, &ll, WRK_USE_AZ); eG(6); // Position berechnen
if (e.flag & SRC_EOUT)
break;
e.z += tau*(m.vz - vsed); // Vertikale Advektion
ClcPtlLoc(&a, &e, &ll, 0); eG(9); // Wieder Position berechnen
if (e.flag & SRC_EOUT)
break;
rc = ClcLocMod(&ll, &mm, -1); eG(7); // vx, vy, vz berechnen
if (rc <= 0) // in verbotenem Gebiet
break;
m.vx = 0.5*(m.vx + mm.vx);
m.vy = 0.5*(m.vy + mm.vy);
m.vz = 0.5*(m.vz + mm.vz);
break;
}
memcpy(&e, &a, PtlRecSize); // Position zuruecksetzen
continue;
}
steps++;
if (TrcMode) trcPtl(np,&e,tp,m);
CountStep++;
if (!ShowProgress && CountStep>=CountMax) {
CountStep = 0;
CountPhase++;
if (CountPhase >= 4) CountPhase = 0;
printf("%c\b", PhaseSymbols[CountPhase]); fflush(stdout);
}
e.refl = 0;
if (a.z > Zref) { /* Oberhalb Zref ... */
e.u = 0; e.v = 0; e.w = 0; } /* ... ohne Diffusion rechnen */
else {
r1 = WrkRnd();
r2 = WrkRnd();
r3 = WrkRnd();
e.u = m.pxx*a.u + m.pxy*a.v + m.pxz*a.w
+ m.wx + m.lxx*r1;
e.v = m.pyx*a.u + m.pyy*a.v + m.pyz*a.w
+ m.wy + m.lyx*r1 + m.lyy*r2;
e.w = m.pzx*a.u + m.pzy*a.v + m.pzz*a.w
+ m.wz + m.lzx*r1 + m.lzy*r2 + m.lzz*r3;
}
e.x = a.x + tau*m.vx; /* Horizontale ... */
e.y = a.y + tau*m.vy; /* ... Advektion durchfuehren */
e.z = a.z;
vxdrop = m.vx; //-2023-07-17
vydrop = m.vy; //-2023-07-17
dzp = tau*(m.vz + e.w - vsed); /* Vertikalen Transport merken */
if (retry > MinRetry && vsed > 0)
dzp += tau*vsed; //-2005-02-24
if (e.afuhgt > 0) { /* Abgasfahnen-Ueberhoehung */
if (e.afurz < 0. && afufac > 0.) {
/* Schwergas nach VDI 3783/1 */
if ((e.afutsc <= m.tau) || (ABS(e.afuhgt) < 0.01)) {
dafu = e.afuhgt;
dzp += dafu*e.afurz; /* Rest aufaddieren */
e.y += dafu*e.afury;
e.x += dafu*e.afurx;
e.afuhgt = 0;
}
else { /* Teil aufaddieren */
// src_rf: used to define fh*Lc
if (src_rf <= 0.) eX(56);
float red = pow(e.h/src_rf, 3.); //-2024-01-17
if (red > 1.) red = 1.;
dafu = e.afuhgt*tau/e.afutsc;
dzp += dafu*e.afurz*red; // reduced shift
e.y += dafu*e.afury;
e.x += dafu*e.afurx;
e.afuhgt -= dafu; // no reduction to guarantee user-defined decay time
}
}
else {
/* Standardverfahren */
if ((e.afutsc <= m.tau) || (ABS(e.afuhgt) < 1.00)) {
dafu = e.afuhgt;
dzp += dafu*e.afurz; /* Rest aufaddieren */
e.y += dafu*e.afury;
e.x += dafu*e.afurx;
e.afuhgt = 0;
}
else { /* Teil aufaddieren */
dafu = e.afuhgt*tau/e.afutsc;
dzp += dafu*e.afurz; /* im z-System */
e.y += dafu*e.afury;
e.x += dafu*e.afurx;
e.afuhgt -= dafu;
}
}
}
if (retry == 0) { /* Beim ersten Durchlauf ... */
ll = l; /* alte gueltige Position merken */
e.flag = SRC_EXIST; /* ... und Flags zuruecksetzen. */
ClcPtlLoc(&a, &e, &l, WRK_USE_AZ); eG(6); /* Position berechnen */
if (e.flag & SRC_EOUT) { /* Rechengebiet verlassen, */
a.flag = e.flag; /* vermerken ... */
break; } /* ... und aufhoeren */
if (INSOLID(l.i,l.j,l.k)) { /* Innerhalb eines Koerpers ? */
e.x = a.x; e.y = a.y; /* Ort zuruecksetzen */
l = ll; /* Alte Position uebernehmen */
m.vx = 0; m.vy = 0; /* Advektion ausschalten */
}
}
e.x += tau*e.u; /* turbulenter Transport */
e.y += tau*e.v;
e.z += dzp;
ClcPtlLoc(&a, &e, &l, 0); eG(9); /* Wieder Position berechnen */
if (e.flag & SRC_EOUT) { /* Rechengebiet verlassen, */
a.flag = e.flag; /* vermerken ... */
break; } /* ... und aufhoeren */
kk = (e.z > Zref) ? -1 : 0;
if (INSOLID(l.i,l.j,l.k)) rc = 0; /* Innerhalb eines Koerpers ? */
else {
rc = ClcLocMod(&l, &m, kk); eG(7); /* Felder interpolieren */
}
if (rc == 0) { /* In verbotenem Gebiet */
retry++; /* Anzahl Versuche pruefen */
if (retry > MaxRetry) {
a.flag |= SRC_REMOVE;
fprintf(MsgFile,
"*** %5d: %6.2f (%5.3lf,%5.3lf,%5.3lf) (%5.3f,%5.3f,%5.3f)"
" F(%5.3f,%5.3f,%5.3f)\n",
np, m.tau, a.x, a.y, a.z, a.u, a.v, a.w, m.vx, m.vy, m.vz);
break; }
if (ll.k > l.k) e.refl |= WRK_REFB; /* Kommt von oben : Deposition ! */
if (retry < MinRetry) { /* Bei den ersten Malen */
a.u=-a.u; a.v=-a.v; a.w=-a.w; } /* ... Teilchen umdrehen, */
else { /* sonst */
a.u = 0; a.v = 0; a.w = 0; /* Eigengeschwindigkeit loeschen */
m.wx = 0; m.wy = 0; m.wz = 0; /* Driftgeschwindigkeit loeschen */
m.vx = 0; m.vy = 0; m.vz = 0; } /* Windgeschwindigkeit loeschen */
continue; /* ... und noch einmal versuchen */
}
//-2011-10-13
if (Pxfr) { // Chemical reactions
XFRREC *prec;
int ir, ixfr;
for (icmp=0; icmp ir = icmp + cmpoff; // get the row index
ixfr = Ixfr[ir]; // get the first list index
if (ixfr < 0) // no reactions for this substance
continue;
gsum = a.g[2*icmp]; // get the old mass
while (ixfr < Nxfr) { // scan the reaction list
prec = Pxfr + ixfr; // get a reaction record
if (prec->ir != ir) // check the row index
break;
jcmp = prec->ic - cmpoff; // index of reaction partner
if (jcmp < 0 || jcmp >= cmpnum) eX(8);
gsum += tau*a.g[2*jcmp]*prec->val; // add contribution
ixfr++;
}
if (gsum < 0) eX(16);
e.g[2*icmp] = gsum; // ... new mass
}
}
//
if (e.refl & (WRK_REFB | WRK_REFT)) { /* Teilchen wurde reflektiert */
if (retry == 0) { /* Falls am Boden reflektiert, */
e.u = -e.u; /* alle Komponenten umkehren. */
e.v = -e.v;
e.w = -e.w /* - 2*m.wz */;
}
if (e.afuhgt > 0) {
if (afufac >= 0)
e.afuhgt *= afufac;
else {
e.afuhgt *= -afufac;
e.afurz *= -1;
}
}
if (deposition && e.refl&WRK_REFB && l.i>0 && l.j>0) {
k = 0;
pf = AryPtr(Pda, l.i, l.j, k, cmpoff); if (!pf) eX(20);
if (Pwa) {
pw = AryPtr(Pwa, l.i, l.j, cmpoff); if (!pw) eX(21);
}
if (wdep_unknown) {
ClcLocMod(&l, &m, 1); eG(7); /* ezz berechnen */
for (icmp=0; icmp if (Pwa) vdep = 0.8*m.ezz*pw[icmp]/(2 - pw[icmp]);
else vdep = pr[icmp].vdep;
PpmClcWdep(m.ezz, vdep+vsed, vsed, wdeps+icmp);
}
wdep_unknown = 0;
}
for (icmp=0; icmp if (e.vg > 0) wdep = wdeps[icmp];
else wdep = (Pwa) ? pw[icmp] : pr[icmp].wdep;
if (wdep < 0) eX(55); //-2013-07-08
if (wdep > 1) wdep = 1;
dg = e.g[2*icmp]*wdep; /* abgelagerte Masse ... */
e.g[2*icmp] -= dg; /* fehlt beim Teilchen ... */
pf[icmp] += dg; } /* und liegt am Boden. */
}
}
if (l.i>0 && l.j>0 && l.k>0 && l.k<=nzdos) {
pf = AryPtr(Pda, l.i, l.j, l.k, cmpoff); if (!pf) eX(12);
for (icmp=0; icmp pf[icmp] += e.g[2*icmp]*tau; /* ... Dosis aufaddieren */
}
if (washout && l.i>0 && l.j>0) { //---------------------- wet deposition
int outside=0; //-2023-07-17
float xdrop; //-2023-07-17
float ydrop; //-2023-07-17
if (WetDrift) {
// record deposited mass with drop drift -2023-07-17
float hdm; // height above ground
float xdm; // centre point x
float ydm; // centre point y
int idrop; // i-index
int jdrop; // j-index
hdm = 0.5*(a.z + e.z);
xdm = 0.5*(a.x + e.x);
ydm = 0.5*(a.y + e.y);
if (bDdmet) {
float h = 0;
GrdBottom(xdm, ydm, &h);
hdm -= h;
}
xdrop = xdm + WRK_DRIFT_VRED*(vxdrop + a.u)*(hdm/vdrop);
ydrop = ydm + WRK_DRIFT_VRED*(vydrop + a.v)*(hdm/vdrop);
idrop = 1 + floor(xdrop - Xmin)/Dd;
jdrop = 1 + floor(ydrop - Ymin)/Dd;
if (idropimax || jdropjmax) {
// store in an outer net
outside = 1;
for (int i=0; i wetbuf[i] = 0;
pf = wetbuf + cmpoff;
}
else {
pf = AryPtr(Pda, idrop, jdrop, -1, cmpoff); if (!pf) eX(10);
}
} // end DRIFT
else {
pf = AryPtr(Pda, l.i, l.j, -1, cmpoff); if (!pf) eX(10);
}
for (icmp=0; icmp float g = e.g[2*icmp];
float fac; //-2011-10-17
fac = pr[icmp].rwsh*tau;
if (fac > 1) fac = 1; //-2011-10-17
dg = g*fac; //-2010-10-17
g -= dg;
e.g[2*icmp] = g;
pf[icmp] += dg;
}
if (outside) {
WetDepoRegister(xdrop, ydrop, wetbuf); //-2023-07-17
}
} // ---------------------------------------------- end of wet deposition
discard = 1; // Noch weiter brauchbar?
for (icmp=0; icmp if (pr[icmp].mmin == 0 //-2012-04-09
|| e.g[2*icmp] > pr[icmp].mmin*e.g[2*icmp+1])
discard = 0; // Nicht abgereicherte weiter verwenden
}
if (discard) { // Leeres Teilchen als nicht ...
a.flag |= SRC_REMOVE; // ... mehr verwendbar markieren.
Ndsc++;
break; } /* und abbrechen. */
retry = 0; /* Versuchs-Zaehler zuruecksetzen */
tp += tau; /* Zeit weitersetzen */
memcpy(&a, &e, PtlRecSize); /* Anfangswerte setzen */
a.flag = SRC_EXIST; /* Flags zuruecksetzen */
a.refl = 0;
predicting = predict;
} /*-------------------------------- Ende der Zeitschleife ---------*/
a.rnd = WrkRndGetSeed(); /* Zufallszahl merken */
a.t = tp; /* Endzeit eintragen */ //-2019-03-14
memcpy(pp, &a, PtlRecSize); /* Werte in der Tabelle eintragen */
if (TrcMode) trcPtl(np,&a,tp,m);
} /*++++++++++++++++++++++++++++++++++ Teilchen fertig +++++++++++++++*/
PtlEnd(handle);
FREE(pa);
FREE(pe);
FREE(wetbuf); //-2023-07-17
vLOG(5)("WRK:RunPtl() returning, n=%d", NumRun);
return NumRun;
eX_30:
eMSG(_no_memory_);
eX_1:
eMSG(_no_handle_);
eX_2:
if (MsgFile != NULL)
fprintf(MsgFile, _current_interval_$$_, T1, T2);
prnPtl(MsgFile, "PTL", np+1, &a, (double)T1); //-2004-01-12
eMSG(_forgotten_particle_);
eX_3: eX_6: eX_9:
prnPtl(MsgFile, "PTL", np, &a, (double)T1); //-2004-01-12
eMSG(_invalid_position_);
eX_4: eX_7:
eMSG(_cant_interpolate_);
eX_5:
eMSG(_no_plume_rise_);
eX_8:
eMSG(_no_transfer_data_);
eX_10:
eX_12:
eX_20:
eMSG(_index_$$$_, l.i, l.j, k);
eX_21:
eMSG(_index_$$_, l.i, l.j);
eX_13:
eMSG(_invalid_source_count_$_, srcind);
eX_14:
eMSG(_inside_building_$$$_, pp->x, pp->y, pp->h);
eX_15:
prnPtl(stdout, "PTL", np, &a, (double)T1);
eMSG(_moved_inside_);
eX_16:
eMSG(_chemical_reactions_too_fast_);
eX_55:
eMSG(_negative_wdep_);
eX_56:
eMSG(_invalid_rf_);
}
#undef a
#undef e
/*================================================================= SetSimParm
*/
static void show_progress( long t ) { //-2001-07-10
long p;
p = (1000.0*(t-MI.start))/(MI.end-MI.start);
if (p != Progress && myid == 0) { //bp
//bp printf(_progress_$_, 0.1*p); //!!!!!!!!!!!!!!!!!!!!!!!!!
// printf("Proc %d %s: Done: %5.1f %%\n", myid ,MsgDate(), 0.1*p); //bp
if (start_sec == 0) start_sec = time(NULL); //bp
long diff = time(NULL) - start_sec; //bp
printf("Done: %d sec %4.1f %%\n", diff, 0.1*p); //bp
fflush(stdout);
Progress = p;
}
}
static long SetSimParm(
void )
{
dP(SetSimParm);
int modlen, i, j, k;
long igp, ima, ipi, iga, isa, ida, ipp, ixa, iwa, usage, ibn;
ARYDSC *pa;
char name[256], t1s[40], t2s[40]; //-2004-11-26
TXTSTR hdr = { NULL, 0 };
strcpy(t1s, TmString(&T1));
strcpy(t2s, TmString(&T2));
vLOG(5)("WRK: SetSimParm() [%s,%s] ", t1s, t2s);
TrcFile = MsgFile;
strcpy(name, NmsName(StdIdent));
igp = IDENT(GRDpar, 0, Gl, Gi);
pa = TmnAttach(igp, NULL, NULL, 0, NULL); eG(1);
if (!pa) eX(2);
Pgp = (GRDPARM*) pa->start;
Nx = Pgp->nx;
Ny = Pgp->ny;
Nz = Pgp->nz;
if (MI.flags & FLG_GAMMA) Nzm = MIN(Pgp->nzmap, Nz);
else Nzm = MIN(Pgp->nzdos, Nz);
Dd = Pgp->dd;
Xmin = Pgp->xmin;
Gx = Nx*Dd;
Xmax = Xmin + Gx;
Ymin = Pgp->ymin;
Gy = Ny*Dd;
Ymax = Ymin + Gy;
bPerx = (0 != (Pgp->bd & GRD_PERX));
if (bPerx) {
Xbmin = Xmin;
Xbmax = Xmax; }
else {
Xbmin = Xmin - Pgp->border;
Xbmax = Xmax + Pgp->border; }
bPery = (0 != (Pgp->bd & GRD_PERY));
if (bPery) {
Ybmin = Ymin;
Ybmax = Ymax; }
else {
Ybmin = Ymin - Pgp->border;
Ybmax = Ymax + Pgp->border; }
DefMode = Pgp->defmode;
PrfMode = Pgp->prfmode;
switch (PrfMode) {
case 2: modlen = sizeof(MD2REC);
break;
case 3: modlen = sizeof(MD3REC);
break;
default: eX(20);
}
Zscl = Pgp->zscl;
//-2005-01-12
if (Zscl != 0) eX(90);
Zref = Pgp->hmax;
iga = IDENT(GRDarr, 0, 0, 0);
pa = TmnAttach(iga, NULL, NULL, 0, NULL); eG(4);
if (!pa) eX(5);
Sk = (float*) pa->start;
TmnDetach(iga, NULL, NULL, 0, NULL); eG(6);
ipi = IDENT(PRMpar, 0, 0, 0);
PrmT2 = T1;
pa = TmnAttach(ipi, &T1, &PrmT2, 0, NULL); eG(7);
if (!pa) eX(8);
if (PrmT2 < T2) T2 = PrmT2;
if (T2 <= T1) eX(41);
if (ShowProgress) show_progress(T1);
precep = BlmPprm->Precep; //-2023-07-17
vdrop = (precep <= 0) ? 0 : WRK_DRIFT_VDROP*pow(precep, WRK_DRIFT_VDROPEXP); //-2023-07-17
Ppi = (PRMINF*) pa->start;
Nc = Ppi->sumcmp;
Flags = Ppi->flags;
TrcMode = (0 != (Flags & FLG_TRACING)); //-2002-01-05 lj
OnlyAdvec = (0 != (Flags & FLG_ONLYADVEC));
NoPredict = (0 != (Flags & FLG_NOPREDICT));
ExactTau = (0 != (Flags & FLG_EXACTTAU));
WetDrift = (0 != (Flags & FLG_WETDRIFT)); //-2023-07-17
TmnDetach(ipi, NULL, NULL, 0, NULL); eG(9);
ModT2 = T1;
ima = IDENT(MODarr, 0, Gl, Gi);
Pma = TmnAttach(ima, &T1, &ModT2, 0, &hdr); eG(10); //-2011-06-29
if (!Pma) eX(11);
bDdmet= (Pma->numdm == 3);
if (bDdmet) {
AryAssert(Pma, modlen, 3, 0, Nx, 0, Ny, 0, Nz); eG(22);
DmnDatAssert("Delta|Dd|Delt", hdr.s, Dd); eG(23); //-2011-06-29
DmnDatAssert("Xmin", hdr.s, Xmin); eG(24); //-2011-06-29
DmnDatAssert("Ymin", hdr.s, Ymin); eG(25); //-2011-06-29
ibn = IDENT(GRDarr, GRD_ISLD, Gl, Gi);
Pbn = TmnAttach(ibn, NULL, NULL, TMN_MODIFY, NULL); if (!Pbn) eX(50);
for (i=1; i<=Nx; i++)
for (j=1; j<=Ny; j++)
for (k=1; k<=Nz; k++)
if (((MD3REC*)AryPtrX(Pma,i,j,k-1))->vz <= WND_VOLOUT)
*(long*)AryPtrX(Pbn,i,j,k) |= GRD_VOLOUT;
}
else {
AryAssert(Pma, modlen, 1, 0, Nz); eG(26);
Pbn = NULL;
}
TmnInfo(ima, NULL, NULL, &usage, NULL, NULL); //-2001-06-29
Valid = (0 == (usage & TMN_INVALID)); //-2001-06-29
if (!Valid && DropT1 != T1) { //-2011-07-04
InterDrop += (T2 - T1);
AverDrop += (T2 - T1);
TotalDrop += (T2 - T1);
DropT1 = T1;
}
Lm = 0;
DmnGetFloat(hdr.s, "LM", "%f", &Lm, 1);
Z0 = 0.5;
DmnGetFloat(hdr.s, "Z0", "%f", &Z0, 1);
D0 = 6*Z0;
DmnGetFloat(hdr.s, "D0", "%f", &D0, 1); //-2018-10-04
Us = 0.3;
DmnGetFloat(hdr.s, "US", "%f", &Us, 1);
DmnGetFloat(hdr.s, "KL|KM", "%f", &Stab, 1);
Hm = 0;
DmnGetFloat(hdr.s, "HM", "%f", &Hm, 1);
Ta = 10.;
DmnGetFloat(hdr.s, "TA", "%f", &Ta, 1); //-2018-10-04
Rh = 70.;
DmnGetFloat(hdr.s, "RH", "%f", &Rh, 1); //-2018-10-04
if (bDdmet) {
if ( (Pgp->sscl == Sk[Nz])
&& (Pgp->bd & GRD_REFZ)
&& (Zref == Zscl)
&& (Zscl > 0));
else Zref = 9999;
}
else {
if (Zref <= 0) {
if (Pgp->bd & GRD_REFZ) Zref = Hm;
if (Zref <= 0) Zref = 9999;
}
else {
if (Zrefbd&GRD_REFZ)) Zref = Hm;
}
if (Zref <= Pgp->zmax) eX(51);
}
TmnDetach(igp, NULL, NULL, 0, NULL); eG(3);
if (ModT2 < T2) T2 = ModT2;
if (T2 <= T1) eX(42);
PpmT2 = T1;
ipp = IDENT(PPMarr, 0, Gl, Gi);
Ppp = TmnAttach(ipp, &T1, &PpmT2, 0, NULL); eG(12);
if (!Ppp) eX(13);
if (PpmT2 < T2) T2 = PpmT2;
if (T2 <= T1) eX(43);
if (bDdmet) {
WdpT2 = T1;
iwa = IDENT(WDParr, 0, Gl, Gi);
Pwa = TmnAttach(iwa, &T1, &WdpT2, 0, &hdr); eG(16);
if (!Pwa) eX(17);
if (WdpT2 < T2) T2 = WdpT2;
if (T2 <= T1) eX(44);
AryAssert(Pwa, sizeof(float), 3, 1, Nx, 1, Ny, 0, Nc-1); eG(27);
DmnDatAssert("Delta|Delt|Dd", hdr.s, Dd); eG(28); //-2011-06-29
DmnDatAssert("Xmin", hdr.s, Xmin); eG(29); //-2011-06-29
DmnDatAssert("Ymin", hdr.s, Ymin); eG(30); //-2011-06-29
}
else Pwa = NULL;
if (Flags & FLG_CHEM) {
float *pf;
XfrT2 = T1;
ixa = IDENT(CHMarr, 0, 0, 0); //-2001-04-25
Pxa = TmnAttach(ixa, &T1, &XfrT2, 0, NULL); eG(31);
if (!Pxa) eX(32);
if (XfrT2 < T2) T2 = XfrT2;
if (T2 <= T1) eX(45);
AryAssert(Pxa, sizeof(float), 2, 0, Nc-1, 0, Nc-1); eG(33);
// -2011-10-13
pf = (float*)Pxa->start;
Nxfr = 0;
for (i=0; i for (j=0; j if (*pf)
Nxfr++;
pf++;
}
if (Nxfr > 0) {
XFRREC *pxfr;
int nr=0;
Pxfr = ALLOC(Nxfr*sizeof(XFRREC));
Ixfr = ALLOC(Nc*sizeof(int));
pxfr = Pxfr;
pf = (float*)Pxa->start;
for (i=0; i Ixfr[i] = -1;
for (j=0; j if (*pf) {
pxfr->ir = i;
pxfr->ic = j;
pxfr->val = *pf;
if (Ixfr[i] < 0)
Ixfr[i] = nr;
pxfr++;
nr++;
}
pf++;
}
}
}
}
else {
Pxa = NULL;
Pxfr = NULL;
Nxfr = 0;
}
SrcT2 = T1;
isa = IDENT(SRCtab, 0, 0, 0);
Psa = TmnAttach(isa, &T1, &SrcT2, 0, NULL); eG(34);
if (!Psa) eX(35);
if (SrcT2 < T2) T2 = SrcT2;
if (T2 <= T1) eX(46);
//
if (Flags & FLG_PLURIS) { //-2018-10-04
int nsrc;
SRCREC *psrc;
nsrc = Psa->bound[0].hgh + 1;
psrc = (SRCREC*)Psa->start;
for (i=0; i int rc = ClcPlr(psrc + i);
vLOG(5)("WRK: ClcPlr('%s') = %d\n", psrc[i].name, rc);
}
}
//
ida = IDENT(DOSarr, Gp, Gl, Gi);
strcpy(name, NmsName(ida));
Pda = NULL;
if (TmnInfo(ida, &DosT1, &DosT2, &usage, NULL, NULL)) {
if (usage & TMN_DEFINED) {
if (DosT1>T1 || DosT2 Pda = TmnAttach(ida, NULL, NULL, TMN_MODIFY, NULL); eG(36);
if (!Pda) eX(39);
}
}
if (!Pda) {
Pda = TmnCreate(ida, sizeof(float),
4, 1, Nx, 1, Ny, -1, Nzm, 0, Nc-1); eG(37);
DosT1 = T1;
DosT2 = T2;
}
TxtClr(&hdr);
if (TrcMode) trcInit(Gp, Gl, Gi, T1, T2);
NumRemove = 0;
Ndsc = 0;
// -2023-07-17
if (GrdPtab != NULL) {
NETREC *pnet = (NETREC*)GrdPtab->start; // largest grid
WetXmin = pnet->xmin;
WetXmax = pnet->xmax;
WetYmin = pnet->ymin;
WetYmax = pnet->ymax;
}
else {
WetXmin = GrdPprm->xmin;
WetXmax = GrdPprm->xmax;
WetYmin = GrdPprm->ymin;
WetYmax = GrdPprm->ymax;
}
numnet = GrdPprm->numnet; //-2023-07-31
//
vLOG(5)("WRK:SetSimParm() returning");
return 0;
eX_1: eX_2: eX_3:
eMSG(_no_grid_data_);
eX_4: eX_5: eX_6:
eMSG(_no_vertical_grid_);
eX_7: eX_8: eX_9:
eMSG(_no_general_parameters_);
eX_10: eX_11:
strcpy(name, NmsName(ima));
eMSG(_no_model_field_$_, name);
eX_12: eX_13:
strcpy(name, NmsName(ipp));
eMSG(_no_particle_parameters_$_, name);
eX_16: eX_17:
strcpy(name, NmsName(iwa));
eMSG(_no_depo_prob_$_, name);
eX_20:
eMSG(_invalid_mode_$_, PrfMode);
eX_22: eX_23: eX_24: eX_25: eX_26:
strcpy(name, NmsName(ima));
eMSG(_improper_model_$_, name);
eX_28: eX_29: eX_30:
if (hdr.s) nMSG("hdr=%s<<<", hdr.s);
eX_27:
strcpy(name, NmsName(iwa));
eMSG(_improper_structure_$_, name);
eX_31: eX_32:
strcpy(name, NmsName(ixa));
eMSG(_no_transfer_$_, name);
eX_33:
strcpy(name, NmsName(iwa));
eMSG(_improper_transfer_$_, name);
eX_34: eX_35:
eMSG(_cant_get_plume_rise_);
eX_36: eX_37: eX_39:
eMSG(_no_dose_$_, name);
eX_38:
strcpy(t1s, TmString(&T1));
strcpy(t2s, TmString(&T2));
eMSG(_dose_$$$_parm_$$_, name, t1s, TmString(&DosT2), t1s, t2s);
eX_41: eX_42: eX_43: eX_44: eX_45: eX_46:
strcpy(t1s, TmString(&T1));
strcpy(t2s, TmString(&T2));
eMSG(_invalid_interval_$$_, t1s, t2s);
eX_50:
eMSG(_no_boundaries_);
eX_51:
eMSG(_reflection_$_$$$_, Zref, Gl, Gi, Pgp->zmax);
eX_90:
eMSG(_zscl_sscl_0_);
}
//================================================================ FreeSimParm
//
static long FreeSimParm( /* Simulations-Parameter freigeben. */
void )
{
dP(FreeSimParm);
long mode;
TmnDetach(IDENT(MODarr, 0, Gl, Gi), NULL, NULL, 0, NULL); eG(41);
Pma = NULL;
TmnDetach(IDENT(PPMarr, 0, Gl, Gi), NULL, NULL, 0, NULL); eG(42);
Ppp = NULL;
TmnDetach(IDENT(SRCtab, 0, 0, 0), NULL, NULL, 0, NULL); eG(43);
Psa = NULL;
if (Pwa) {
TmnDetach(IDENT(WDParr, 0, Gl, Gi), NULL, NULL, 0, NULL); eG(44);
Pwa = NULL;
}
if (Pxa) {
TmnDetach(IDENT(CHMarr, 0, 0, 0), NULL, NULL, 0, NULL); eG(45);
Pxa = NULL;
if (Pxfr) {
FREE(Pxfr); //-2011-10-13
Pxfr = NULL;
}
if (Ixfr) {
FREE(Ixfr); //-2011-10-13
Ixfr = NULL;
}
}
mode = TMN_MODIFY | ((Valid) ? TMN_SETVALID : TMN_INVALID); //-2001-06-29
TmnDetach(IDENT(DOSarr, Gp, Gl, Gi), &DosT1, &T2, mode, NULL); eG(46);
Pda = NULL;
if (Pbn) {
TmnDetach(IDENT(GRDarr, GRD_ISLD, Gl, Gi), NULL, NULL, TMN_MODIFY, NULL); eG(47);
Pbn = NULL;
}
if (NumRemove)
vLOG(4)("%d particles from inside bodies removed", NumRemove);
if (Ndsc)
vLOG(4)("%d particles discarded", Ndsc);
return 0;
eX_41: eX_42: eX_43: eX_44: eX_45: eX_46: eX_47:
eMSG(_cant_detach_);
}
//================================================================== WrkHeader
//
#define CAT(a,b) {TxtPrintf(&txt,(a),(b)); TxtCat(&hdr,txt.s);}
char *WrkHeader( /* the header (global storage) */
long id, /* identification */
long *pt1, /* start of the validity time */
long *pt2 ) /* end of validity time */
{
dP(WrkHeader);
char name[256], t1s[40], t2s[40], dts[40], rdat[40], format[80], vldf[40]; //-2008-12-04 //-2011-04-13 //-2018-10-04
int gl, gi, iga, igp, k, ict, nz, nc, seq;
char atype;
float *sk;
ARYDSC *pa, *pct, *pga;
GRDPARM *pgp;
CMPREC *pcmp;
enum DATA_TYPE dt;
TXTSTR txt = { NULL, 0 };
TXTSTR hdr = { NULL, 0 };
strcpy(name, NmsName(id));
dt = XTR_DTYPE(id);
gl = XTR_LEVEL(id);
gi = XTR_GRIDN(id);
switch (dt) {
case DOSarr: atype = 'D';
break;
case SUMarr: atype = 'D';
break;
case DOStab: atype = 'D';
break;
case CONtab: atype = 'C';
break;
case GAMarr: atype = 'G';
break;
default: atype = '?';
}
igp = IDENT(GRDpar, 0, gl, gi);
pa = TmnAttach(igp, NULL, NULL, 0, NULL); eG(1);
if (!pa) eX(2);
pgp = (GRDPARM*) pa->start;
iga = IDENT(GRDarr, 0, 0, 0);
pga = TmnAttach(iga, NULL, NULL, 0, NULL); eG(3);
if (!pga) eX(4);
sk = (float*) pga->start;
nz = pga->bound[0].hgh + 1;
ict = IDENT(CMPtab, 0, 0, 0);
pct = TmnAttach(ict, NULL, NULL, 0, NULL); eG(5);
if (!pct) eX(6);
pcmp = pct->start;
nc = pct->bound[0].hgh + 1;
strcpy(t1s, TmString(pt1));
strcpy(t2s, TmString(pt2));
strcpy(dts, TmString(&MI.cycle));
strcpy(format, "Con%10.2e");
strcpy(vldf, "V");
if (MI.flags & FLG_SCAT) {
strcat(format, "Dev%(*100)5.1f"); //-2011-06-29
strcat(vldf, "V");
}
seq = MI.cycind;
if (MI.average > 1) seq = 1 + (seq-1)/MI.average;
//
// dmna header //-2011-06-29
TxtCpy(&hdr, "\n");
// sprintf(s, "TALWRK_%d.%d.%s", StdVersion, StdRelease, StdPatch);
CAT("prgm \"%s\"\n", TalLabel); //-2011-12-05
CAT("idnt \"%s\"\n", MI.label);
CAT("artp \"%c\"\n", atype);
CAT("axes \"%s\"\n", "xyzs");
CAT("vldf \"%s\"\n", vldf);
CAT("file \"%s\"\n", name);
CAT("form \"%s\"\n", format);
CAT("t1 \"%s\"\n", t1s);
CAT("t2 \"%s\"\n", t2s);
CAT("dt \"%s\"\n", dts);
if (TmGetDate(MI.rbrk) < TM_MAX_DATE) {
TmFormatDate(rdat, "%24D", MI.rbrk);
CAT("rdat \"%s\"\n", rdat);
}
CAT("index %d\n", seq);
CAT("groups %d\n", MI.numgrp);
CAT("refx %d\n", pgp->refx);
CAT("refy %d\n", pgp->refy);
if (pgp->refx > 0 && *pgp->ggcs) //-2008-12-04
CAT("ggcs \"%s\"\n", pgp->ggcs); //-2008-12-04 //-2008-12-04
CAT("xmin %s\n", TxtForm(pgp->xmin, 6));
CAT("ymin %s\n", TxtForm(pgp->ymin, 6));
CAT("delta %s\n", TxtForm(pgp->dd, 6));
CAT("zscl %s\n", TxtForm(pgp->zscl, 6));
CAT("sscl %s\n", TxtForm(pgp->sscl, 6));
TxtCat(&hdr, "sk ");
for (k=0; k CAT(" %s", TxtForm(sk[k], 6));
TxtCat(&hdr, "\n");
TxtCat(&hdr, "name ");
for (k=0; k CAT(" \"%s\"", pcmp[k].name);
TxtCat(&hdr, "\n");
TxtCat(&hdr, "unit ");
for (k=0; k CAT(" \"%s\"", pcmp[k].unit);
TxtCat(&hdr, "\n");
TxtCat(&hdr, "refc ");
for (k=0; k CAT(" %10.2e", pcmp[k].refc);
TxtCat(&hdr, "\n");
TxtCat(&hdr, "refd ");
for (k=0; k CAT(" %10.2e", pcmp[k].refd);
TxtCat(&hdr, "\n");
CAT("valid %1.6f\n", WrkValid(VALID_AVER)); //-2011-07-04
TmnAttach(id, NULL, NULL, TMN_MODIFY, NULL); eG(7); //-2001-06-29
TmnDetach(id, NULL, NULL, TMN_MODIFY|TMN_SETVALID, NULL); eG(13);
TmnDetach(igp, NULL, NULL, 0, NULL); eG(10);
TmnDetach(iga, NULL, NULL, 0, NULL); eG(11);
TmnDetach(ict, NULL, NULL, 0, NULL); eG(12);
TxtClr(&txt);
return hdr.s;
eX_1: eX_2: eX_3: eX_4: eX_5: eX_6: eX_7:
nMSG(_cant_attach_);
return NULL;
eX_10: eX_11: eX_12: eX_13:
nMSG(_cant_detach_);
return NULL;
}
#undef CAT
//=================================================================== WrkValid
float WrkValid( int type ) { //-2011-07-04
dP(WrkValid);
float valid;
if (type == VALID_TOTAL) {
if (TotalTime <= 0) eX(1);
valid = (TotalTime - TotalDrop)/((float)TotalTime);
}
else if (type == VALID_AVER) {
if (AverTime <= 0) eX(2);
valid = (AverTime - AverDrop)/((float)AverTime);
}
else if (type == VALID_INTER) {
if (InterTime != MI.cycle) eX(3);
valid = (InterTime - InterDrop)/((float)InterTime);
}
else eX(4);
if (valid < 0) eX(5);
if (valid > 1) eX(6);
return valid;
eX_1: eX_2: eX_3: eX_4: eX_5: eX_6:
eMSG(_error_calculating_valid_);
return 0;
}
//==================================================================== WrkInit
//
long WrkInit( /* initialize server */
long flags, /* action flags */
char *istr ) /* server options */
{
dP(WrkInit);
long iddos, mask;
char *jstr, *ps, s[200];
int vrb, dsp;
if (StdStatus & STD_INIT) return 0;
if (istr) {
jstr = istr;
ps = strstr(istr, "-v");
if (ps) sscanf(ps+2, "%d", &StdLogLevel);
ps = strstr(istr, "-p");
if (ps) sscanf(ps+2, "%d", &ShowProgress);
ps = strstr(istr, "-y");
if (ps) sscanf(ps+2, "%d", &StdDspLevel);
ps = strstr(istr, "-d");
if (ps) strcpy(AltName, ps+2);
}
else jstr = "";
vLOG(3)("WRK_%d.%d.%s (%08lx,%s)", StdVersion, StdRelease, StdPatch, flags, jstr);
StdStatus |= flags;
iddos = IDENT(DOSarr, 0, 0, 0);
mask = ~(NMS_GROUP|NMS_GRIDN|NMS_LEVEL);
TmnCreator(iddos, mask, TMN_UNIQUE, "", WrkServer, WrkHeader); eG(1);
vrb = StdLogLevel;
dsp = StdDspLevel;
sprintf(s, " GRD -v%d -y%d -d%s", vrb, dsp, AltName);
GrdInit(flags, s); eG(2);
sprintf(s, " PRM -v%d -y%d -d%s", vrb, dsp, AltName);
PrmInit(flags, s); eG(3);
sprintf(s, " PTL -v%d -y%d -d%s", vrb, dsp, AltName);
PtlInit(flags, s); eG(4);
sprintf(s, " SRC -v%d -y%d -d%s", vrb, dsp, AltName);
SrcInit(flags, s); eG(5);
sprintf(s, " MOD -v%d -y%d -d%s", vrb, dsp, AltName);
ModInit(flags, s); eG(6);
sprintf(s, " PPM -v%d -y%d -d%s", vrb, dsp, AltName);
PpmInit(flags, s); eG(7);
if (MI.flags & FLG_PLURIS) {
vLOG(2)("PLURIS %s (bf=%1.1f,stdw=%d)", str_version(), PlrBreakFactor, PlrStDownwash);
}
StdStatus |= STD_INIT;
return 0;
eX_1:
eMSG(_creator_$_, NmsName(iddos));
eX_2: eX_3: eX_4: eX_5: eX_6: eX_7:
eMSG(_cant_init_servers_);
}
//================================================================= WrkServer
//
long WrkServer( /* server routine for DOSarr */
char *s ) /* calling option */
{
dP(WrkServer);
char t1s[40], t2s[40], name[256]; //-2004-11-26
int ida, ipa, mask;
if (StdArg(s)) return 0;
if (*s) {
switch (s[1]) {
case 'd': strcpy(AltName, s+2);
break;
case 'J': sscanf(s+2, "%d", &TrcMode);
break;
case 'n': if (s[2]) { //-2011-07-04
sscanf(s+2, "%d,%d", &AverTime, &AverDrop);
if (AverTime < 0) {
AverTime = 0;
AverDrop = 0;
TotalTime = 0;
TotalDrop = 0;
}
}
else {
TotalTime += MI.cycle;
AverTime += MI.cycle;
}
DropT1 = -1;
InterTime = MI.cycle;
InterDrop = 0;
break;
case 'o': if (strstr(s+2, "nostdw"))
PlrStDownwash = 0;
else if (strstr(s+2, "logpluris"))
LogPluris = 1;
break;
default: eX(17);
}
return 0;
}
if (!StdIdent) eX(20);
if ((StdStatus & STD_INIT) == 0) {
WrkInit(0, ""); eG(20);
}
Gl = XTR_LEVEL(StdIdent);
Gi = XTR_GRIDN(StdIdent);
Gp = XTR_GROUP(StdIdent);
if (StdStatus & STD_TIME) T1 = StdTime;
else eX(1);
strcpy(t1s, TmString(&T1));
ida = IDENT(DOSarr, 0, 0, 0);
mask = ~(NMS_GROUP|NMS_GRIDN|NMS_LEVEL);
TmnCreator(ida, mask, 0, "", NULL, NULL); eG(2);
PtlT2 = T1;
ipa = IDENT(PTLarr, Gp, Gl, Gi);
Ppa = TmnAttach(ipa, &T1, &PtlT2, TMN_MODIFY, NULL); eG(10);
if (!Ppa) eX(11);
T2 = PtlT2; if (T2 <= T1) eX(3);
Nptl = PtlCount(ipa, SRC_EXIST, Ppa); eG(12);
SetSimParm(); eG(13);
strcpy(t2s, TmString(&T2));
strcpy(name, NmsName(StdIdent));
vDSP(2)(_moving_$$$_, name, t1s, t2s);
RunPtl(); eG(14);
vDSP(2)("");
vLOG(4)("WRK: %d particles moved for %s [%s,%s]",
NumRun, name, t1s, t2s);
FreeSimParm(); eG(15);
TmnDetach(ipa, &T1, &PtlT2, TMN_MODIFY, NULL); eG(16);
TmnCreator(ida, mask, TMN_UNIQUE, "", WrkServer, WrkHeader); eG(21);
return 0;
eX_20:
eMSG(_missing_data_);
eX_1:
eMSG(_missing_time_);
eX_2: eX_21:
strcpy(name, NmsName(ida));
eMSG(_cant_change_creator_$_, name);
eX_3:
strcpy(name, NmsName(ipa));
strcpy(t2s, TmString(&T2));
eMSG(_invalid_times_$$$_, name, t1s, t2s);
eX_10: eX_11: eX_12:
strcpy(name, NmsName(ipa));
eMSG(_no_particles_$_, name);
eX_13: eX_15:
eMSG(_no_parameters_);
eX_14:
strcpy(name, NmsName(ida));
strcpy(t2s, TmString(&T2));
eMSG(_cant_build_$$$_, name, t1s, t2s);
eX_16:
strcpy(name, NmsName(ipa));
strcpy(t1s, TmString(&T2));
strcpy(t2s, TmString(&PtlT2));
eMSG(_cant_update_$$$_, name, t1s, t2s);
eX_17:
eMSG("unknown option %s!", s+1);
}
/*==============================================================================
*
* history:
*
* 2002-06-21 lj 0.13.0 final test version
* 2002-09-24 lj 1.0.0 final release candidate
* 2002-12-10 lj 1.0.4 source parameter Tt
* 2003-02-21 lj 1.1.2 internal random number generator
* 2003-06-16 lj 1.1.7 variabel vsed
* 2004-01-12 lj 1.1.14 error message more precise
* 2004-10-01 uj 2.1.0 version number upgraded
* 2004-10-10 uj 2.1.1 Try2Escape
* 2004-11-26 lj 2.1.7 string length for names = 256
* 2004-12-09 uj 2.1.8 shift particle outside body also if not created
* 2005-01-12 uj 2.1.9 check if Zscl = 0
* 2005-01-23 uj 2.1.12 re-mark particle created outside an inner grid
* 2005-02-24 uj 2.1.15 disable sedimentation for retry > MinRetry
* 2005-03-17 uj 2.2.0 version number upgrade
* 2006-10-26 lj 2.3.0 external strings
* 2008-08-04 lj 2.4.3 writes "valid" using 6 decimals
* 2008-12-04 lj 2.4.5 writes "refx", "refy", "ggcs", "rdat"
* 2011-04-13 uj 2.4.9 handle title with 256 chars in header
* 2011-06-29 uj 2.5.0 adjusted to dmna output
* 2011-07-04 uj handling of "valid" extended
* 2011-10-13 uj 2.5.2 optimization for chemical reactions
* 2011-11-21 lj 2.6.0 wet deposition
* 2013-07-08 uj 2.6.8 check for wdep smaller 0
* 2018-10-04 uj 3.0.0 pluris
* 2019-03-14 uj 3.0.3 end time of a particle as double
* 2021-06-23 uj 3.1.1 modification dense gas coupling for VDI 3783/1
* 2021-07-03 uj linux adjustment
* 2022-10-04 uj 3.1.5 ClcPlr: check if source is outside grid(s)
* 2023-07-17 uj 3.2.0 drop drift for wet deposition
* 2023-07-31 uj 3.2.1 drop drift corrections
* 2024-01-17 uj 3.3.0 BESMAX consistent calculation
* 2024-01-17 updated for heavy gas VDI 3783/1 E
*
*============================================================================*/
*
* Particle driver for AUSTAL
* ==========================
*
* Copyright (C) Umweltbundesamt, Dessau-Roßlau, Germany, 2002-2024
* Copyright (C) Janicke Consulting, 88662 Überlingen, Germany, 2002-2024
* Email: info@austal.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.
*
* last change: 2024-01-17 uj
*
* ===========================================================================*/
#include
#ifdef __linux__ //-2021-07-03
#include
#endif
#define STDMYMAIN WrkMain
#include "IBJmsg.h"
#include "IBJary.h"
#include "IBJtxt.h"
#include "IBJdmn.h"
#include "IBJstd.h"
static char *eMODn = "TalWrk";
//=========================================================================
STDPGM(tstwrk, WrkServer, 3, 3, 0);
//=========================================================================
#include "genutl.h"
//#include "genio.h"
#include "TalWnd.h"
#include "TalPrf.h"
#include "TalGrd.h"
#include "TalMat.h"
#include "TalMod.h"
#include "TalNms.h"
//#include "TalPlm.h"
#include "TalPpm.h"
#include "TalPrm.h"
#include "TalPtl.h"
#include "TalSrc.h"
#include "TalTmn.h"
#include "TalDef.h" //-2011-12-05
#include "TalWrk.h"
#include "TalWrk.nls"
#include "pluris.h" //-2018-10-04
#include "TalBlm.h" //-2018-10-04
#define LOC 0x02000000L
#ifndef PI
#define PI 3.1415926
#endif
#define WRK_REFB 0x0001
#define WRK_REFT 0x0002
#define WRK_USE_H 0x0100
#define WRK_USE_AZ 0x0200
#define WRK_DRIFT_VDROP 3.6 //-2023-07-17
#define WRK_DRIFT_VDROPEXP 0.16 //-2023-07-17
#define WRK_DRIFT_VRED 0.8 //-2023-07-17
#define INSOLID(a,b,c) ((Pbn) && *(long*)AryPtrX(Pbn,(a),(b),(c))&GRD_VOLOUT)
void** WetDepoList; // needed in TALdos() //-2023-07-17
typedef struct {
int i, j, k, is, js, ks;
float ax, ay, az;
} LOCREC;
typedef struct xfrrec { int ir, ic; float val; //-2011-10-13
} XFRREC;
/*======================== Funktions-Prototypen ============================*/
static int ClcAfu( // Berechnung der Abgasfahnen-Ueberhoehung
SRCREC *ps, // Pointer auf die Quelle
float v2, // Geschwindigkeitsquadrat
float stab, // Stabilitaetsklasse oder Temperatur-Gradient
float *pu, // Ergebnis: Ueberhoehung.
float *pt ); // Ergebnis: Anstiegszeit.
static int ClcPtlLoc( // Teilchen im Netz einordnen, evtl. reflektieren
PTLREC *pa, // Pointer auf alte Teilchen-Koordinaten
PTLREC *pp, // Pointer auf neue Teilchen-Koordinaten
LOCREC *pl, // Pointer auf Index-Werte
int h2z ); // h-Wert nach z-Wert umwandeln
static int ClcLocMod( // Feld am Ort des Teilchens berechnen
LOCREC *pl, // Pointer auf Index-Werte
MD3REC *p3, // Pointer auf Feld-Werte
int clceta ); // <=0, wenn Eta nicht interpoliert werden soll
static int RunPtl( // Teilchen-Trajektorien berechnen.
void );
static int Try2Escape( // In einem Gebaeude erzeugte Partikel herausdraengen
PTLREC *pa,
LOCREC *pl );
/*=================== Alle Groessen nur intern fuer Worker ====================*/
static long start_sec = 0; //bp
static long Flags = 0;
static int Valid = 1;
static int Nptl = 0;
static int Ndsc = 0;
static int NumRun = 0;
static int NumRemove = 0;
static int MaxRetry = 50;
static int MinRetry = 6;
static int Nx, Ny, Nz, Nzm, PrfMode, DefMode, TrcMode, Nc;
static int OnlyAdvec, NoPredict, ExactTau, WetDrift; //-2023-07-17
static long T1, T2, DosT1, DosT2;
static long PrmT2, PpmT2, WdpT2, SrcT2, XfrT2, PtlT2, ModT2;
static float Dd, Xmin, Xmax, Gx, Ymin, Ymax, Gy, Xbmin, Xbmax, Ybmin, Ybmax;
static float Zref, Zscl, Stab, *Sk, Hm, Lm, Z0, D0, Us; //-2018-10-04
static float Ta, Rh; //-2018-10-04
static int LogPluris; //-2018-10-04
static int Gl, Gi, Gp;
static char AltName[256];
static FILE *TrcFile;
static GRDPARM *Pgp; /* Grid-Parameter */
static PRMINF *Ppi; /* General Parameter */
static ARYDSC *Pma; /* Model-Array */
static ARYDSC *Ppp; /* Particle-Parameter */
static ARYDSC *Ppa; /* Particle-Array */
static ARYDSC *Pda; /* Dose-Array */
static ARYDSC *Pwa; /* Wdep-Array */
static ARYDSC *Pxa; /* Xfr-Array (Chem) */
static ARYDSC *Psa; /* Source-Array */
static ARYDSC *Pbn; /* Boundaries-Array */
static int bPerx, bPery, bDdmet;
static int TotalTime, TotalDrop, AverTime, AverDrop, InterTime, InterDrop, DropT1;
static char PhaseSymbols[] = "|/-\\";
static int CountStep, CountPhase, CountMax=50000;
static int ShowProgress, Progress=-1; //-2001-07-10
static int Nxfr; //-2011-10-13
static XFRREC *Pxfr;
static int *Ixfr;
static float precep; //-2023-07-17
static float vdrop; //-2023-07-17
static int numnet; //-2023-07-17
static double WetXmin, WetXmax, WetYmin, WetYmax; //-2023-07-17
//-------------------------------------------------------------
static unsigned long WrkRndSeed = 123456;
static int WrkRndIset = 0;
static unsigned long WrkRndULong() {
WrkRndSeed = 1664525*WrkRndSeed + 1013904223;
return WrkRndSeed;
}
static unsigned long WrkRndGetSeed() {
return WrkRndSeed;
}
static void WrkRndSetSeed(unsigned long seed) {
if (seed != 0) WrkRndSeed = seed;
else WrkRndSeed = 123456;
WrkRndIset = 0;
}
static float WrkRndFloat() {
unsigned long r = WrkRndULong();
float f = ((r>>8) & 0x00ffffff)/((float)0x01000000);
return f;
}
#define WrkRnd01() WrkRndFloat() // random numbers
#define WrkRnd() WrkRndNrm() // use normal distribution
//================================================================== WrkRndNrm
//
static float WrkRndNrm() // normal distribution, average=0, variance=1 -2008-12-11
{
static float gset;
float fac, r, v1, v2;
if (WrkRndIset == 0) {
do {
r = WrkRndFloat();
v1 = 2.0*r - 1.0;
r = WrkRndFloat();
v2 = 2.0*r - 1.0;
r = v1*v1 + v2*v2;
} while (r >= 1.0);
fac = sqrt( -2.0*log(r)/r );
gset = v1*fac;
WrkRndIset = 1;
return v2*fac;
}
else {
WrkRndIset = 0;
return gset;
}
}
static double WrkClock( // Get the clock in seconds
double base) // reference base
{ // RETURN clock - base
double clk;
#ifdef __linux__
static long tv_sec = 0;
static long tv_usec = 0;
struct timeval tval;
gettimeofday(&tval, NULL);
clk = tval.tv_sec + tval.tv_usec/1000000.0;
#else
clk = clock()/(double)CLOCKS_PER_SEC;
#endif
return clk - base;
}
/*--------------------------------------------------------------------------*/
#define a (*pa)
static void prnPtl( FILE *prn, char *txt, int np, PTLREC *pa, double tp )
{
int i;
fprintf(prn, "%s %4d: t=%1.0lf flg=%02x rfl=%02x src=%ld off=%d rnd=%ld tp=%1.1lf\n", //-2021-08-05
txt, np, a.t, a.flag, a.refl, a.srcind, a.offcmp, a.rnd, tp );
fprintf(prn, " x=(%11.3lf,%11.3lf,%11.3lf/%11.3lf) v=(%7.3f,%7.3f,%7.3f) \n",
a.x, a.y, a.z, a.h, a.u, a.v, a.w );
fprintf(prn, " afu=(%10.2e,%10.2e) \n", a.afuhgt, a.afutsc );
fprintf(prn, " g=(");
for (i=0; i
}
static void trcPtl( int np, PTLREC *pa, double tp, MD3REC m )
{
int i;
fprintf(TrcFile,
"# %3d %ld %d %02x %02x: %8.2lf | %8.2lf %8.2lf %7.3lf | "
"%6.3f %6.3f %6.3f | %6.3f %6.3f %6.3f |",
np, a.srcind, a.offcmp, a.flag, a.refl, tp, a.x, a.y, a.z,
m.vx, m.vy, m.vz, a.u, a.v, a.w );
for (i=0; i
}
#undef a
static void trcInit( int gp, int nl, int ni, long t1, long t2 )
{
fprintf( TrcFile,
"#### %d %d %d [%ld,%ld]\n", gp, nl, ni, t1, t2 );
}
/*===================================================================== ClcAfu
*/
static int ClcAfu( /* Berechnung der Abgasfahnen-Ueberhoeung */
SRCREC *ps, /* Pointer auf die Quelle */
float v2, /* Geschwindigkeitsquadrat */
float stab, /* Stabilitaetsklasse oder Temperatur-Gradient */
float *pu, /* Ergebnis: Ueberhoeung. */
float *pt ) /* Ergebnis: Anstiegszeit. */
{
dP(ClcAfu);
//float uq, dh, xe;
//int m, rc;
if (!ps) eX(1);
if (ps->ts >= 0) {
*pt = ps->ts;
*pu = ps->ts*ps->vq;
}
else if (ps->hrise >= 0) { //-2018-10-04
*pt = ps->trise;
*pu = ps->hrise;
}
/*
else {
m = PLM_VDI3782_3;
uq = sqrt(v2);
rc = TalPlm(m, &dh, &xe, ps->qq, ps->h, uq, stab,
ps->vq, ps->dq, ps->rhy, ps->lwc, ps->tmp); //-2002-12-10
if (rc < 0) {
if (rc == -1 && (PlmMsg)) eX(3); //-2002-09-23
eX(2);
}
*pu = dh;
if (uq < 1.0) uq = 1.0;
if (ps->qq>0 || ps->tmp>10) *pt = 0.4*xe/uq; //-2002-12-10
else if (ps->vq > 0) *pt = dh/ps->vq;
else *pt = 0;
}
*/
return 0;
eX_1:
eMSG(_no_source_data_);
/*
eX_3:
vMsg("%s", PlmMsg);
PlmMsg = NULL;
eX_2:
eMSG(_cant_calculate_plume_rise_);
* */
}
/*================================================================= ClcAfuRxyz
*/
static int ClcAfuRxyz( SRCREC *ps, float *pvs, float *prx, float *pry, float *prz )
{
dP(ClcAfuRxyz);
float sif, cof, sig, cog, dq;
float rx, ry, rz, ax, ay, az, r, vs, fq, gq, sh, sv, sl;
*prx = 0;
*pry = 0;
*prz = 1;
if (!ps) eX(1);
if (ps->tmp>0 && !(Flags & FLG_PLURIS)) //-2018-10-04
return 0;
fq = ps->fq;
sh = ps->sh;
sv = ps->sv;
sl = ps->sl;
if (fq==0 && sh<=0 && sv<=0 && sl<=0)
return 0;
vs = *pvs;
dq = ps->wq;
cof = cos(fq*PI/180);
sif = sin(fq*PI/180);
if (ps->gq > -999) {
gq = ps->gq; // gq: meterological orientation
sig = cos(gq*PI/180); // sig: mathematical orientation
cog = sin(gq*PI/180);
}
else {
gq = (ps->aq >= 0) ? 90-dq : dq;
cog = ps->co;
sig = ps->si;
}
rz = vs*cof; //-99-08-15
ry = vs*sif*sig;
rx = vs*sif*cog;
if (sh > 0) { //-99-08-05
ax = sig;
ay = -cog;
r = sh*WrkRndNrm();
rx += r*ax;
ry += r*ay;
}
if (sv > 0) {
ax = -cof*cog;
ay = -cof*sig;
az = sif;
r = sv*WrkRndNrm();
rx += r*ax;
ry += r*ay;
rz += r*az;
}
if (sl > 0) {
ax = sif*cog;
ay = sif*sig;
az = cof;
r = sl*WrkRndNrm();
rx += r*ax;
ry += r*ay;
rz += r*az;
}
r = sqrt(rx*rx + ry*ry + rz*rz);
*pvs = r;
if (r > 0) {
*prx = rx/r;
*pry = ry/r;
*prz = rz/r;
}
return 1;
eX_1:
eMSG(_no_source_data_);
}
/*======================================================================= ClcPlr
*/ //-2018-10-04
static int ClcPlr( // Call PLURIS to get the plume rise for this source
SRCREC *ps) // the source
{
dP(ClcPlr);
PLRSRC src;
PLRPRM prm;
PLRRSL rsl;
PLRPVL *ppvl=NULL;
PR3REC *p0, *p1;
ARYDSC *pprf=NULL;
char name[512];
//
int se = 100000;
int sn = 1000;
float sd = 0.01;
int adjust_step = 200;
float min_us = 0.05;
//
int k, k0, k1, nh, rc;
float h, vx, vy, v, d, tmp;
double *hhprf = NULL; // layer height at source position
long iprf=0, last_valid;
double c_0, c_1, c_2;
//
// explicit plume rise, nothing to do
if (ps->ts >= 0)
return 0;
//
// last calculated plume rise still valid, nothing to do
if (ps->valid >= fmin(SrcT2, ModT2)) //-2018-10-04
return 1;
//
// initialize to zero plume rise
ps->trise = 0;
ps->hrise = 0;
ps->sh = 0;
ps->sl = 0;
ps->sv = 0;
rsl.vr = 0; //-2024-01-17
rsl.tr = 0; //-2024-01-17
rsl.dr = 0; //-2024-01-17
//
// no valid meteorology
if (!Valid)
return 0;
//
// no volume flow, zero plume rise
last_valid = ps->valid;
if (last_valid < 0) last_valid = 0;
if (ps->vq <= 0) {
ps->valid = fmin(SrcT2, ModT2);
if (LogPluris) {
char t1s[40], t2s[40];
strcpy(t1s, TmString(&last_valid));
strcpy(t2s, TmString(&(ps->valid)));
vLOG(1)("WRK: [%s,%s] pluris('%s') = %s", t1s, t2s, ps->name, str_rsl(rsl));
}
return 0;
}
//
// prepare plume rise calculation with PLURIS
//
// check source parameters
if (ps->dq <= 0) eX(10);
if (ps->vq <= 0) eX(11);
if (ps->fq != 0) eX(12);
if (ps->rt < 0 || ps->rt > 1) eX(14);
if (ps->lwc < 0 || ps->lwc > 1) eX(21);
if (ps->vwc < 0 || ps->vwc > 1) eX(24);
if (ps->rhy < 0 || ps->rhy > 100) eX(22);
if (ps->rhy > 0 && ps->rhy <= 1) eX(15);
if (Rh < 0 || Rh > 100) eX(23);
if (Rh > 0 && Rh <= 1.0) eX(16);
if (Ta < 0 || Ta > 50) eX(17); //-2024-01-17
if (ps->rf < 0 || ps->rf > 1) eX(25); //-2024-01-17
//
c_0 = WrkClock(0);
vLOG(5)("WRK: ClcPlr('%s')", ps->name);
//
// profile heights Sk[] above ground (possibly reduced)
if (Gl > 1) { // nested grids
NETREC *pn = (NETREC*)GrdPtab->start; // outermost grid
iprf = IDENT(PRFarr, 0, (int)(pn->level), (int)(pn->index));
}
else { // single grid
iprf = IDENT(PRFarr, 0, Gl, Gi);
}
pprf = TmnAttach(iprf, NULL, NULL, 0, NULL); eG(106);
if (!pprf) eX(107);
k = (bDdmet) ? 2 : 0;
nh = pprf->bound[k].hgh - pprf->bound[k].low + 1;
TmnDetach(iprf, NULL, NULL, 0, NULL); eG(108);
pprf = NULL;
ppvl = (PLRPVL*)ALLOC(nh*sizeof(PLRPVL)); if (!ppvl) eX(3);
//
// source parameters
tmp = ps->tmp;
if (tmp < Ta)
tmp = Ta;
src.dm = ps->dq;
src.hb = ps->h;
src.zq = ps->wld;
src.lq = ps->lwc;
src.sq = ps->vwc;
src.rq = ps->rhy * 0.01;
src.te = tmp;
src.ve = ps->vq;
src.us = Us;
src.lm = Lm;
src.ta = Ta;
src.rf = ps->rf; //-2024-01-17
//
// calculation parameters
prm.sc = src.dm; // scale = source diameter
prm.sd = sd; // scaled step width
prm.se = se; // scaled maximum distance (ss)
prm.sn = sn; // number of log output lines
prm.adjust_step = adjust_step; // >0: dynamic step rr/adjust_step
prm.min_us = min_us; // minimum u* for break condition
prm.const_ust = 0;
prm.bf = PLRBREAKFACTOR;
if (MI.flags & FLG_BESMAX) { //-2024-01-17
prm.const_ust = 1;
prm.bf = PLRBREAKFACTOR_BESMAX;
}
vLOG(4)("WRK: ClcPlr, sc=%1.3f, sd=%1.4f, se=%1.0f, sn=%d, adjust_step=%d, min_us=%1.1f, bf=%1.2f, const_ust=%d\n",
prm.sc, prm.sd, prm.se, prm.sn, prm.adjust_step, prm.min_us, prm.bf, prm.const_ust);
//
// insert profile values
iprf = 0;
//
// 1D meteorology
if (!bDdmet) {
float a=0;
int nk;
//
if (Gl > 1) { // nested grids
NETREC *pn = (NETREC*)GrdPtab->start; // outermost grid
iprf = IDENT(PRFarr, 0, (int)(pn->level), (int)(pn->index));
}
else { // current (single) grid
iprf = IDENT(PRFarr, 0, Gl, Gi);
}
pprf = TmnAttach(iprf, NULL, NULL, 0, NULL); eG(6);
if (!pprf) eX(7);
nk = pprf->bound[0].hgh - pprf->bound[0].low + 1;
for (k=0; k
//
// for 1D meteorology, use vertical grid points Sk
// for defining the vertical profiles for PLURIS
h = Sk[k];
//
// find the vertical interval [k0, k1] containing h
for (k1=1; k1
if (p1->h >= h)
break;
}
if (k1 >= nk) eX(5);
k0 = k1 - 1;
p0 = (PR3REC*)AryPtr(pprf, k0); if (!p0) eX(101);
a = (h - p0->h)/(p1->h - p0->h);
vx = p0->vx + a*(p1->vx - p0->vx);
vy = p0->vy + a*(p1->vy - p0->vy);
v = sqrt(vx*vx + vy*vy);
d = (v == 0) ? 0 : atan2(-vx, -vy); // clockwise from north
d *= 180/PI;
if (d < 0) d += 360.;
ppvl[k].ha = h;
ppvl[k].ua = v;
ppvl[k].da = d;
}
TmnDetach(iprf, NULL, NULL, 0, NULL); eG(102);
pprf = NULL;
}
//
// 3D meteorology
else {
int in, nesting;
int ni, nj, nk; // number of points
double xmin, xmax, ymin, ymax, dd;
double xq, yq, ax, ay;
double zbot, ztop;
int iq, jq; // point index
int knext; // next valid k-value of ppvl
double hnext, hk; // next required height of PVL-array
int grid_found = 0;
NETREC *pn;
PR3REC *p000, *p100, *p010, *p110;
PR3REC *p001, *p101, *p011, *p111;
//
xq = ps->x;
yq = ps->y;
nesting = (GrdPtab != NULL);
hnext = 0;
//
in = GrdPprm->numnet;
knext = 0;
while(1) { // step through all grids, starting with the finest grid
in--;
vLOG(5)("WRK: ClcPlr, grid index = %d", in);
if (nesting) {
if (in < 0) {
nh = knext;
vLOG(5)("WRK: ClcPlr, grid restricts profiles to height %1.1f", Sk[nh-1]);
break;
}
pn = ((NETREC*)GrdPtab->start) + in;
xmin = pn->xmin;
xmax = pn->xmax;
ymin = pn->ymin;
ymax = pn->ymax;
if (xq < xmin || xq > xmax || yq < ymin || yq > ymax)
continue; // try the next coarser grid
grid_found = 1;
dd = pn->delta;
iprf = IDENT(PRFarr, 0, (int)(pn->level), (int)(pn->index));
}
else { // no nesting, use the current (single) grid
if (in < -1) {
nh = knext;
vLOG(5)("WRK: ClcPlr, grid restricts profiles to height %1.1f", Sk[nh-1]);
break;
}
xmin = GrdPprm->xmin;
xmax = GrdPprm->xmax;
ymin = GrdPprm->ymin;
ymax = GrdPprm->ymax;
if (xq < xmin || xq > xmax || yq < ymin || yq > ymax)
break;
grid_found = 1;
dd = GrdPprm->dd;
iprf = IDENT(PRFarr, 0, 0, 0);
}
strcpy(name, NmsName(iprf));
vLOG(5)("WRK: ClcPlr, using profile %s", name);
pprf = TmnAttach(iprf, NULL, NULL, 0, NULL); eG(103);
if (pprf == NULL) eX(104);
ni = pprf->bound[0].hgh - pprf->bound[0].low + 1;
nj = pprf->bound[1].hgh - pprf->bound[1].low + 1;
nk = pprf->bound[2].hgh - pprf->bound[2].low + 1;
// find the grid cell containing the source position
iq = (int)floor((xq - xmin)/dd);
if (iq == ni-1)
iq = ni-2;
jq = (int)floor((yq - ymin)/dd);
if (jq == nj-1)
jq = nj-2;
// just in case, should not happen
if (iq < 0 || iq > ni-1 || jq < 0 || jq > nj-1) eX(201); //-2022-10-04
ax = (xq - xmin - iq*dd)/dd;
ay = (yq - ymin - jq*dd)/dd;
p000 = AryPtr(pprf, iq, jq, 0);
p100 = AryPtr(pprf, iq+1, jq, 0);
p010 = AryPtr(pprf, iq, jq+1, 0);
p110 = AryPtr(pprf, iq+1, jq+1, 0);
zbot = (1-ax)*(1-ay)*p000->h + ax*(1-ay)*p100->h + (1-ax)*ay*p010->h
+ ax*ay*p110->h;
p001 = AryPtr(pprf, iq, jq, nk-1);
p101 = AryPtr(pprf, iq+1, jq, nk-1);
p011 = AryPtr(pprf, iq, jq+1, nk-1);
p111 = AryPtr(pprf, iq+1, jq+1, nk-1);
ztop = (1-ax)*(1-ay)*p001->h + ax*(1-ay)*p101->h + (1-ax)*ay*p011->h
+ ax*ay*p111->h;
vLOG(5)("WRK: ClcPlr, zbot=%e, ztop=%e, hnext=%e", zbot, ztop, hnext);
if (hnext > ztop-zbot) { // profile not suited
if (nesting) {
TmnDetach(iprf, NULL, NULL, 0, NULL); eG(105);
pprf = NULL;
continue;
}
else
break;
}
//
// Calculate the height above ground of the profile layers at the
// source position.
//
hhprf = ALLOC(nk*sizeof(double));
hhprf[0] = 0;
hhprf[nk-1] = ztop - zbot;
for (k=1; k
p100 = AryPtr(pprf, iq+1, jq, k);
p010 = AryPtr(pprf, iq, jq+1, k);
p110 = AryPtr(pprf, iq+1, jq+1, k);
hhprf[k] = (1-ax)*(1-ay)*p000->h + ax*(1-ay)*p100->h + (1-ax)*ay*p010->h
+ ax*ay*p110->h - zbot;
}
//
// Interpolate profile values at the heights of the PVL array.
//
for (k=knext; k
// for 3D meteorology, use the mid points of the vertical
// grid points Sk for defining the vertical profiles for PLURIS;
// this is consistent with the definition of the Arakawa-C wind
// components as area averages.
hk = (k == 0) ? 0 : 0.5*(Sk[k]+Sk[k-1]); //-2018-10-04
if (hk == 0) {
k1 = 1;
}
else {
for(k1=1; k1
break;
}
if (k1 >= nk) {
knext = k;
hnext = hk;
break;
}
k0 = k1 - 1;
p110 = AryPtr(pprf, iq+1, jq+1, k0);
if (p110->vz <= WND_VOLOUT) {
v = -1;
d = 0;
}
else {
p101 = AryPtr(pprf, iq+1, jq, k1);
p011 = AryPtr(pprf, iq, jq+1, k1);
p111 = AryPtr(pprf, iq+1, jq+1, k1);
vx = (1-ax)*p011->vx + ax*p111->vx; //-2018-10-04
vy = (1-ay)*p101->vy + ay*p111->vy; //-2018-10-04
v = sqrt(vx*vx + vy*vy);
d = (v == 0) ? 0 : atan2(-vx, -vy); // clockwise from north
d *= 180/PI;
if (d < 0) d += 360.;
}
ppvl[k].ha = hk;
ppvl[k].ua = v;
ppvl[k].da = d;
} // for k
FREE(hhprf);
TmnDetach(iprf, NULL, NULL, 0, NULL); eG(102);
pprf = NULL;
if (k >= nh)
break;
} // while
if (!grid_found) eX(200); //-2022-10-04
for (k=nh-1; k>=0; k--) { // insert wind speed in region with body cells
if (ppvl[k].ua >= 0)
continue;
ppvl[k].ua = ppvl[k+1].ua;
ppvl[k].da = ppvl[k+1].da;
}
} // 3d meteorology
//
// set meteorological parameters beside wind
for (k=0; k
ppvl[k].ca = 1;
ppvl[k].ra = Rh*0.01;
}
//
// set wind direction at the ground if required
if (ppvl[0].ua == 0 && ppvl[0].ra == 0) {
ppvl[0].da = ppvl[1].da;
}
//
// run pluris
PlrVerbose = StdLogLevel - 3 + 3*LogPluris;
PlrMessage = MsgFile;
c_1 = WrkClock(c_0);
rc = run_pluris(&src, nh, ppvl, &prm, &rsl); if (rc < 0) eX(4); //-2024-01-17
c_2 = WrkClock(c_0);
//
// rise parameters
if (MI.flags & FLG_BESMAX) { //-2024-01-17
rsl.tr = 0.0;
}
ps->hrise = rsl.dr;
ps->trise = rsl.tr;
//
// set source induced turbulence
ps->sh = rsl.vr * ps->rt;
ps->sl = ps->sh;
ps->sv = ps->sh;
//
if (MI.flags & FLG_BESMAX) { //-2024-01-17
ps->sh = 0.;
ps->sl = 0.;
ps->sv = 0.;
}
ps->valid = fmin(SrcT2, ModT2); //-2018-10-04
if (LogPluris) {
char t1s[40], t2s[40];
strcpy(t1s, TmString(&last_valid));
strcpy(t2s, TmString(&(ps->valid)));
vLOG(1)("WRK: [%s,%s] pluris('%s') = %s Set[rt=%5.3f, sh|sl|sv=%7.3f] Time[%1.4f]",
t1s, t2s, ps->name, str_rsl(rsl), ps->rt, ps->sh, c_2-c_1);
}
FREE(ppvl);
ppvl = NULL;
return 0;
eX_3:
eMSG("can't create metorological profile for source '%s'!", ps->name);
eX_4:
eMSG("PLURIS failed for source '%s' (%d)!", ps->name, rc);
eX_5:
eMSG("vertical extension of meteorological profile '%s' not sufficient!", name);
eX_106:
eMSG("can't get outer meteorological profile '%s'", name);
eX_6:
eMSG("can't get outer meteorological profile '%s'", name);
eX_107:
eMSG("missing outer meteorological profile '%s'", name);
eX_7:
eMSG("missing outer meteorological profile '%s'", name);
eX_10:
eMSG("plume rise calculation requires diameter > 0 for source '%s'!", ps->name);
eX_11:
eMSG("plume rise calculation requires exit velocity > 0 for source '%s'!", ps->name);
eX_12:
eMSG("Fq must be 0 (only plume rise for vertical exit implemented) for source '%s'!", ps->name);
eX_14:
eMSG("Rt (%1.1f) not in range [0,1] for source '%s'!", ps->rt, ps->name);
eX_15:
eMSG("Rh must be specified in percent (0 or in range [1,100]) for source '%s'!", ps->name);
eX_16:
eMSG("Rh of ambient air must be specified in percent (0 or in range [1,100])!");
eX_17:
eMSG("Ta of ambient air must be in range [0,50])!");
eX_21:
eMSG("Lw (%1.3f) not in range [0,1] for source '%s'!", ps->lwc, ps->name);
eX_22:
eMSG("Rh (%1.1f) not in range [0,100] for source '%s'!", ps->rhy, ps->name);
eX_23:
eMSG("Ambient Rh (%1.1f) not in range [0,100]!", Rh);
eX_24:
eMSG("Vw (%1.3f) not in range [0,1] for source '%s'!", ps->vwc, ps->name);
eX_25:
eMSG("Rf (%1.3f) not in range [0,1] for source '%s'!", ps->rf, ps->name);
eX_200:
eMSG("Source '%s' outside grid(s)!", ps->name);
eX_100: eX_101: eX_102: eX_103: eX_104: eX_105: eX_108: eX_201:
eMSG("internal error");
}
//===================================================================ClcPtlLoc
//
static int ClcPtlLoc( /* Teilchen im Netz einordnen, evtl. reflektieren */
PTLREC *pa, /* Pointer auf alte Teilchen-Koordinaten */
PTLREC *pp, /* Pointer auf neue Teilchen-Koordinaten */
LOCREC *pl, /* Pointer auf Index-Werte */
int h2z ) /* h-Wert nach z-Wert umwandeln */
{
// dP(ClcPtlLoc);
int i, j, k, flg, m;
float r, x, y, z, dz, xi, eta, zeta;
float zg;
flg = GRD_REFBOT;
z = (pa) ? pa->z : 0; //-2002-04-10
if ((pa) && (Zref) && (z<=Zref)) flg |= GRD_REFTOP; //-2002-04-10
zeta = 0;
if (h2z == WRK_USE_H) { z = pp->h; flg |= GRD_PARM_IS_H; }
else {
z = pp->z;
if (h2z == WRK_USE_AZ) flg |= GRD_USE_ZETA;
zeta = pl->az;
}
x = pp->x;
i = -1;
if (x < Xmin) {
if (bPerx) {
x = Xmax - fmod(Xmin-x, Gx); //-2001-03-23 lj
pp->x = x;
}
else i = 0;
}
else if (x > Xmax) {
if (bPerx) {
x = Xmin + fmod(x-Xmax, Gx); //-2001-03-23 lj
pp->x = x;
}
else i = 0;
}
y = pp->y;
j = -1;
if (y < Ymin) {
if (bPery) {
y = Ymax - fmod(Ymin-y, Gy); //-2001-03-23 lj
pp->y = y;
}
else j = 0;
}
else if (y > Ymax) {
if (bPery) {
y = Ymin + fmod(y-Ymax, Gy); //-2001-03-23 lj
pp->y = y;
}
else j = 0;
}
xi = (x - Xmin)/Dd;
eta = (y - Ymin)/Dd;
pl->i = 0;
pl->j = 0;
k = pl->k;
pl->k = 0;
if (bDdmet) {
if (i==0 || j==0) { pp->flag |= SRC_GOUT; return 0; }
m = GrdLocate(Pma, xi, eta, &z, flg, Zref, &i, &j, &k, &zg, &zeta);
if (m & GRD_REFTOP) pp->refl |= WRK_REFT;
if (m & GRD_REFBOT) pp->refl |= WRK_REFB;
if (!(m & GRD_INSIDE)) pp->flag |= SRC_GOUT;
}
else {
zg = 0;
if (i==0 || j==0) {
if (x<=Xbmin || x>=Xbmax) pp->flag |= SRC_XOUT;
if (y<=Ybmin || y>=Ybmax) pp->flag |= SRC_YOUT;
}
else {
i = 1 + xi; if (i > Nx) i = Nx;
j = 1 + eta; if (j > Ny) j = Ny;
}
if (z < 0) { z = -z; pp->refl |= WRK_REFB; }
if (flg&GRD_REFTOP && z>Zref) { z = 2*Zref-z; pp->refl |= WRK_REFT; }
for (k=1; k<=Nz; k++)
if (z <= Sk[k]) { //-2002-04-10
dz = Sk[k] - Sk[k-1];
zeta = (z - Sk[k-1])/dz;
break;
}
if (k > Nz) {
k = 0; pp->flag |= SRC_ZOUT; }
}
if (i < 0) i = 0;
if (j < 0) j = 0;
pl->i = i;
pl->j = j;
pl->k = k;
pp->z = z;
if ((i) && (j)) pp->h = z - zg;
if (pp->flag & SRC_GOUT) return 0;
pl->ax = xi - (i-1);
pl->ay = eta - (j-1);
pl->az = zeta;
if (bDdmet) {
r = WrkRnd01();
if (r > pl->ax) pl->is = i-1;
else pl->is = i;
r = WrkRnd01();
if (r > pl->ay) pl->js = j-1;
else pl->js = j;
pl->ks = k;
}
else {
pl->is = i;
pl->js = j;
pl->ks = k;
}
return 1;
}
//================================================================== ClcLocMod
//
static int ClcLocMod( /* Feld am Ort des Teilchens berechnen */
LOCREC *pl, /* Pointer auf Index-Werte */
MD3REC *p3, /* Pointer auf Feld-Werte */
int clceta ) /* <=0, wenn Eta nicht interpoliert werden soll */
{
dP(ClcLocMod);
MD2REC *p2a, *p2b;
MD3REC *p3a, *p3b;
MD3REC *p000, *p001, *p010, *p100, *p011, *p101, *p110, *p111;
float a000, a001, a010, a100, a011, a101, a110, a111;
float ax, bx, ay, by, az, bz;
int i, j, k, is, js;
float a, b;
i = pl->i;
j = pl->j;
k = pl->k;
is = pl->is;
js = pl->js;
if (k == 0) eX(1);
if (bDdmet) {
if (i==0 || j==0) eX(2);
p3a = (MD3REC*) AryPtr(Pma, i, j, k-1); if (!p3a) eX(4);
if (p3a->vz <= WND_VOLOUT) return 0; /* Innerhalb eines Koerpers */
p3b = (MD3REC*) AryPtr(Pma, i, j, k ); if (!p3b) eX(5);
if (p3b->ths <= WND_VOLOUT) return 0; /* Innerhalb eines Koerpers */
if (p3b->wz <= WND_VOLOUT) return 0; /* Innerhalb eines Koerpers */
memset(p3, 0, sizeof(MD3REC)); /* Zuerst alles auf 0 setzen */
if (!OnlyAdvec) {
p3->ths = p3b->ths; /* Gradient der potentiellen Temperatur */
p3->wx = p3b->wx; /* Kompensation, x-Komponente */
p3->wy = p3b->wy; /* Kompensation, y-Komponente */
p3->wz = p3b->wz; /* Kompensation, z-Komponente */
}
if (DefMode & GRD_3LIN) { /* Geschwindigkeiten 3*linear */
p000 = (MD3REC*) AryPtr(Pma, i-1, j-1, k-1); if (!p000) eX(10);
p001 = (MD3REC*) AryPtr(Pma, i-1, j-1, k ); if (!p001) eX(11);
p010 = (MD3REC*) AryPtr(Pma, i-1, j , k-1); if (!p010) eX(12);
p100 = (MD3REC*) AryPtr(Pma, i , j-1, k-1); if (!p100) eX(13);
p011 = (MD3REC*) AryPtr(Pma, i-1, j , k ); if (!p011) eX(14);
p101 = (MD3REC*) AryPtr(Pma, i , j-1, k ); if (!p101) eX(15);
p110 = (MD3REC*) AryPtr(Pma, i , j , k-1); if (!p110) eX(16);
p111 = (MD3REC*) AryPtr(Pma, i , j , k ); if (!p111) eX(17);
bx = pl->ax;
ax = 1.0 - bx;
by = pl->ay;
ay = 1.0 - by;
bz = pl->az;
az = 1.0 - bz;
a000 = ax*ay*az;
a001 = ax*ay*bz;
a010 = ax*by*az;
a100 = bx*ay*az;
a011 = ax*by*bz;
a101 = bx*ay*bz;
a110 = bx*by*az;
a111 = bx*by*bz;
p3->vx = a000*p000->vx + a001*p001->vx + a010*p010->vx + a100*p100->vx +
a011*p011->vx + a101*p101->vx + a110*p110->vx + a111*p111->vx;
p3->vy = a000*p000->vy + a001*p001->vy + a010*p010->vy + a100*p100->vy +
a011*p011->vy + a101*p101->vy + a110*p110->vy + a111*p111->vy;
p3->vz = a000*p000->vz + a001*p001->vz + a010*p010->vz + a100*p100->vz +
a011*p011->vz + a101*p101->vz + a110*p110->vz + a111*p111->vz;
}
else {
p3a = (MD3REC*) AryPtr(Pma, i-1, j, k); if (!p3a) eX(20);
p3->vx = (1-pl->ax)*p3a->vx + pl->ax*p3b->vx; /* Vx interpolieren */
p3a = (MD3REC*) AryPtr(Pma, i, j-1, k); if (!p3a) eX(21);
p3->vy = (1-pl->ay)*p3a->vy + pl->ay*p3b->vy; /* Vy interpolieren */
p3a = (MD3REC*) AryPtr(Pma, i, j, k-1); if (!p3a) eX(22);
p3->vz = (1-pl->az)*p3a->vz + pl->az*p3b->vz; /* Vz interpolieren */
}
if (PrfMode == 3) { /* Anisotrope Modellierung */
a = 1 - pl->az;
b = pl->az;
p3a = (MD3REC*) AryPtr(Pma, is, js, k-1); if (!p3a) eX(23);
p3b = (MD3REC*) AryPtr(Pma, is, js, k ); if (!p3b) eX(24);
p3->tau = a*p3a->tau + b*p3b->tau;
if (OnlyAdvec) return 1;
if (clceta < 0) return 1;
p3->pxx = a*p3a->pxx + b*p3b->pxx;
p3->pxy = a*p3a->pxy + b*p3b->pxy;
p3->pxz = a*p3a->pxz + b*p3b->pxz;
p3->pyx = a*p3a->pyx + b*p3b->pyx;
p3->pyy = a*p3a->pyy + b*p3b->pyy;
p3->pyz = a*p3a->pyz + b*p3b->pyz;
p3->pzx = a*p3a->pzx + b*p3b->pzx;
p3->pzy = a*p3a->pzy + b*p3b->pzy;
p3->pzz = a*p3a->pzz + b*p3b->pzz;
p3->lxx = a*p3a->lxx + b*p3b->lxx;
p3->lyx = a*p3a->lyx + b*p3b->lyx;
p3->lyy = a*p3a->lyy + b*p3b->lyy;
p3->lzx = a*p3a->lzx + b*p3b->lzx;
p3->lzy = a*p3a->lzy + b*p3b->lzy;
p3->lzz = a*p3a->lzz + b*p3b->lzz;
if (clceta == 0) return 1;
p3->exx = a*p3a->exx + b*p3b->exx;
p3->eyx = a*p3a->eyx + b*p3b->eyx;
p3->eyy = a*p3a->eyy + b*p3b->eyy;
p3->ezx = a*p3a->ezx + b*p3b->ezx;
p3->ezy = a*p3a->ezy + b*p3b->ezy;
p3->ezz = a*p3a->ezz + b*p3b->ezz;
return 1;
}
else { /* Isotrope Modellierung */
a = 1 - pl->az;
b = pl->az;
p2a = (MD2REC*) AryPtr(Pma, is, js, k-1); if (!p2a) eX(25);
p2b = (MD2REC*) AryPtr(Pma, is, js, k ); if (!p2b) eX(26);
p3->tau = a*p2a->tau + b*p2b->tau;
if (OnlyAdvec) return 1;
if (clceta < 0) return 1;
p3->pxx = a*p2a->pvv + b*p2b->pvv;
p3->pyy = p3->pxx;
p3->pzz = a*p2a->pww + b*p2b->pww;
p3->lxx = a*p2a->lvv + b*p2b->lvv;
p3->lyy = p3->lxx;
p3->lzz = a*p2a->lww + b*p2b->lww;
if (clceta == 0) return 1;
p3->exx = a*p2a->evv + b*p2b->evv;
p3->eyy = p3->exx;
p3->ezz = a*p2a->eww + b*p2b->eww;
return 1;
}
}
else { /* 1-dimensionale Meteorologie */
if ((k <= Nz) && (k > 0)) {
memset(p3, 0, sizeof(MD3REC)); /* Zuerst alles auf 0 setzen */
a = 1 - pl->az;
b = pl->az;
if (PrfMode == 3) { /* Anisotrope Modellierung */
p3b = (MD3REC*) AryPtr(Pma, k ); if (!p3b) eX(27);
p3a = (MD3REC*) AryPtr(Pma, k-1); if (!p3a) eX(28);
p3->vx = a*p3a->vx + b*p3b->vx;
p3->vy = a*p3a->vy + b*p3b->vy;
p3->vz = a*p3a->vz + b*p3b->vz;
p3->tau = a*p3a->tau + b*p3b->tau;
if (OnlyAdvec) return 1;
p3->wz = p3b->wz;
p3->ths = p3b->ths;
if (clceta < 0) return 1;
p3->pxx = a*p3a->pxx + b*p3b->pxx;
p3->pxy = a*p3a->pxy + b*p3b->pxy;
p3->pxz = a*p3a->pxz + b*p3b->pxz;
p3->pyx = a*p3a->pyx + b*p3b->pyx;
p3->pyy = a*p3a->pyy + b*p3b->pyy;
p3->pyz = a*p3a->pyz + b*p3b->pyz;
p3->pzx = a*p3a->pzx + b*p3b->pzx;
p3->pzy = a*p3a->pzy + b*p3b->pzy;
p3->pzz = a*p3a->pzz + b*p3b->pzz;
p3->lxx = a*p3a->lxx + b*p3b->lxx;
p3->lyx = a*p3a->lyx + b*p3b->lyx;
p3->lyy = a*p3a->lyy + b*p3b->lyy;
p3->lzx = a*p3a->lzx + b*p3b->lzx;
p3->lzy = a*p3a->lzy + b*p3b->lzy;
p3->lzz = a*p3a->lzz + b*p3b->lzz;
if (clceta == 0) return 1;
p3->exx = a*p3a->exx + b*p3b->exx;
p3->eyx = a*p3a->eyx + b*p3b->eyx;
p3->eyy = a*p3a->eyy + b*p3b->eyy;
p3->ezx = a*p3a->ezx + b*p3b->ezx;
p3->ezy = a*p3a->ezy + b*p3b->ezy;
p3->ezz = a*p3a->ezz + b*p3b->ezz;
return 1; }
else { /* Isotrope Modellierung */
p2b = (MD2REC*) AryPtr(Pma, k ); if (!p2b) eX(29);
p2a = (MD2REC*) AryPtr(Pma, k-1); if (!p2a) eX(30);
p3->vx = a*p2a->vx + b*p2b->vx;
p3->vy = a*p2a->vy + b*p2b->vy;
p3->vz = a*p2a->vz + b*p2b->vz;
p3->tau = a*p2a->tau + b*p2b->tau;
if (OnlyAdvec) return 1;
p3->wz = p2b->wz;
p3->ths = p2b->ths;
if (clceta < 0) return 1;
p3->pxx = a*p2a->pvv + b*p2b->pvv;
p3->pyy = p3->pxx;
p3->pzz = a*p2a->pww + b*p2b->pww;
p3->lxx = a*p2a->lvv + b*p2b->lvv;
p3->lyy = p3->lxx;
p3->lzz = a*p2a->lww + b*p2b->lww;
if (clceta == 0) return 1;
p3->exx = a*p2a->evv + b*p2b->evv;
p3->eyy = p3->exx;
p3->ezz = a*p2a->eww + b*p2b->eww;
return 1; }
}
else eG(31); /* k > Nz */
}
return 1;
eX_1: eX_2: eX_4: eX_5:
eX_10: eX_11: eX_12: eX_13: eX_14: eX_15: eX_16: eX_17:
eX_20: eX_21: eX_22: eX_23: eX_24: eX_25: eX_26: eX_27: eX_28: eX_29:
eX_30: eX_31:
eMSG(_loc_$_, pl->i, pl->j, pl->k, pl->is, pl->js, pl->ax, pl->ay, pl->az);
}
//================================================================= Try2Escape
//
static int Try2Escape( /* In einem Gebaeude erzeugte Partikel herausdraengen */
PTLREC *pa,
LOCREC *pl ) {
dP(WRK:Try2Escape);
double x = pa->x; // Anfangs-Position merken
double y = pa->y;
double h = pa->h;
double d = 1;
int ii[9] = { 0, 1, -1, 0, 0, 1, 1, -1, -1 };
int jj[9] = { 0, 0, 0, 1, -1, 1, -1, 1, -1 };
int n, i, k, kk;
for (n=1; n<=10; n++) { // 10 Distanzen pruefen
d = n*0.1*Dd;
for (kk=0; kk<=2; kk++) {
k = (kk == 2) ? -1 : kk;
pa->h = h + k*d;
for (i=(k==0); i<9; i++) {
pa->x = x + ii[i]*d;
pa->y = y + jj[i]*d;
ClcPtlLoc(NULL, pa, pl, WRK_USE_H); eG(2); // neu einordnen
if (!INSOLID(pl->i,pl->j,pl->k)) { // geschafft!
d = 0;
break;
}
} // for i
if (d == 0) break;
}
if (d == 0) break;
}
if (d > 0) eX(3);
return 0;
eX_2: eMSG(_cant_locate_); //-2004-12-23
eX_3: eMSG(_cant_shift_); //-2004-12-23
}
//----------------------------------------------------------------------------
// utilities for rain drop drift -2023-07-17
//----------------------------------------------------------------------------
/**
* Register deposited masses
*/
static int WetDepoRegister(
float xdrop, // x-coordinate of wet deposition
float ydrop, // y-coordinate of wet deposition
float* wdrops) // deposition values [Nc]
{
dP(WetDepoRegister);
int inet, net, ic, ibuf, idrop, jdrop, irec, nd, nx, ny;
float *buffer_field = NULL;
float *record;
NETREC* pndst;
//
// deposition inside calculation area?
if (GrdPprm->numnet == 0)
return 0;
if (xdrop
return 0;
//
// identify target grid: only check grids not yet processed
for (net=GrdPprm->net-1; net>0; net--) {
pndst = (NETREC*) AryPtr(GrdPtab, net);
if (xdrop >= pndst->xmin && xdrop < pndst->xmax
&& ydrop >= pndst->ymin && ydrop
break;
}
inet = net - 1;
if (inet < 0)
return -1;
ibuf = Gp*numnet + inet;
nx = floor(pndst->nx + 0.5); //-2023-07-31
ny = floor(pndst->ny + 0.5); //-2023-07-31
if (WetDepoList[ibuf] == NULL) {
nd = nx*ny*Nc*sizeof(float);
WetDepoList[ibuf] = ALLOC(nd);
if (!WetDepoList[ibuf]) eX(1);
}
idrop = floor((xdrop - pndst->xmin)/pndst->delta);
jdrop = floor((ydrop - pndst->ymin)/pndst->delta);
if (idrop < 0) idrop = 0; //-2023-07-31
if (idrop > nx-1) idrop = nx-1; //-2023-07-31
if (jdrop < 0) jdrop = 0; //-2023-07-31
if (jdrop > ny-1) jdrop = ny-1; //-2023-07-31
buffer_field = WetDepoList[ibuf];
irec = Nc*(jdrop + ny*idrop);
record = buffer_field + irec;
for (ic=0; ic
return 0;
eX_1: eMSG("internal error");
}
/*
* Addition of deposition values stored in WetDepoList to the values of
* the current dose array DOSarr.
*/
static int WetDepoAdd() {
dP(WetDepoAdd);
int inet = Pgp->net - 1;
if (inet < 0)
return -1;
int ibuf = Gp*numnet + inet;
float *pwd = (float*)WetDepoList[ibuf];
if (pwd == NULL)
return 0;
for (int i=0; i
float* pdos = AryPtr(Pda, i+1, j+1, -1, 0);
for (int ic=0; ic
padd[ic] = 0;
}
}
}
return 0;
}
static int WetDepoPrint(int inet, int gp) {
int ibuf = gp*numnet + inet;
float *dd = (float*)WetDepoList[ibuf];
NETREC* pn = (NETREC*)AryPtr(GrdPtab, inet+1);
int nx = pn->nx;
int ny = pn->ny;
if (dd == NULL)
return -1;
for (int j=ny-1; j>=0; j--) {
printf("%3d: ", j+1);
for (int i=0; i
printf("%8.2e ", d);
}
printf("\n");
}
return 0;
}
//----------------------------------------------------------------------------
// end utilities for rain drop drift
//----------------------------------------------------------------------------
//===================================================================== RunPtl
//
static int RunPtl(
void )
{
dP(RunPtl);
long rc;
static int totals = 0;
static int steps = 0;
PTLREC *pp, *pa, *pe;
#define a (*pa)
#define e (*pe)
LOCREC l, ll;
MD3REC m, mm;
PPMREC *pr;
SRCREC *ps;
float dg, gsum, r1, r2, r3, vsed, vdep, *pf, *pw, dafu;
float wdep, v2, tau, vred, hred, dzp, afufac;
float *wetbuf; //-2023-07-17
float vxdrop, vydrop; //-2023-07-17
float src_rf; //-2024-01-17
float wdeps[60];
double tp;
int handle, retry;
int icmp, jcmp, cmpnum, cmpoff, nzdos, np, kk, k, srcind;
int imin, imax, jmin, jmax; //-2023-07-17
int deposition, discard, created, wdep_unknown;
int washout; //-2011-11-21
int predict, predicting; //-2001-09-22
vLOG(5)("WRK:RunPtl() starting");
if (WetDrift) { //-2023-07-17
// transfer wet deposition that entered into the current grid by drop drift
// from buffer WetDepoList to this dose array DOSarr.
WetDepoAdd(); // applies Pgp, Pda, numnet
}
wetbuf = ALLOC(Nc*sizeof(float));
imin = Pda->bound[0].low;
imax = Pda->bound[0].hgh;
jmin = Pda->bound[1].low;
jmax = Pda->bound[1].hgh;
//
pa = ALLOC(PtlRecSize); if (!pa) eX(30);
pe = ALLOC(PtlRecSize); if (!pe) eX(30);
memset(&l, 0, sizeof(l)); /* Strukturen zu 0 initialisieren */
memset(&m, 0, sizeof(m));
memset(&e, 0, PtlRecSize);
nzdos = Pda->bound[2].hgh;
handle = PtlStart(Ppa); eG(1);
NumRun = 0;
predict = (NoPredict) ? 0 : (bDdmet != 0); //-2001-09-22
for (np=0; ; ) { /*+++++++++++++++++ Fuer alle Teilchen: ++++++++++++++++*/
pp = PtlNext(handle); /* Neues Teilchen anfordern */
if (pp == NULL) break; /* aufhoeren, wenn keins mehr da ist */
totals++;
if (pp->flag & SRC_REMOVE) /* naechstes nehmen, wenn aussortiert */
continue;
if (!Valid) { // Keine gueltige Meteorologie:
pp->flag |= SRC_REMOVE; // Alle vorhandenen Partikel werden
NumRemove++; // geloescht.
continue;
}
memcpy(&a, pp, PtlRecSize); /* umspeichern */
if (a.t >= T2) continue; /* naechstes nehmen, wenn abgearbeitet */
if (a.t < T1) eX(2);
tp = a.t; /* Zeit separat merken */
created = (a.flag & SRC_CREATED) != 0; /* neu erzeugt, falls Flag gesetzt */ //-2003-02-21
a.flag = SRC_EXIST;
WrkRndSetSeed(a.rnd); // Zufallszahlen initialisieren //-2003-02-21
cmpoff = a.offcmp; /* Komponenten-Offset im Dosis-Feld */
cmpnum = a.numcmp; /* Zahl der Komponenten */
srcind = a.srcind - 1; /* Index der zugehoerigen Quelle */
ps = (SRCREC*) AryPtr(Psa, srcind); // Pointer auf die zugehoerige Quelle
if (!ps) eX(13);
afufac = ps->fr; // Faktor bei Reflektion
src_rf = ps->rf; // fuer VDI 3673/1 Schwergas -2024-01-17
pr = (PPMREC*) Ppp->start + cmpoff;
deposition = 0;
washout = 0; //-2011-11-21
for (icmp=0; icmp
if (pr[icmp].rwsh > 0) washout = 1; //-2011-11-21
}
if (Pwa) deposition = 1;
wdep_unknown = 0;
if (a.vg > 0) { /* Falls individuelles Vsed */
vsed = a.vg;
deposition = 1;
wdep_unknown = 1;
}
else vsed = pr->vsed; /* Sedimentations-Geschwindigkeit */
vred = pr->vred; /* Reduktion der vert. Startgeschwindigk.*/
hred = pr->hred; /* Reduktion der hor. Startgeschwindigk. */
np++; /* Teilchenzaehler hochzaehlen ... */
NumRun = np; /* ... und extern zur Verfuegung stellen */
ClcPtlLoc(NULL, &a, &l, WRK_USE_H); eG(3); /* Teilchen einordnen */
if (a.flag & SRC_EOUT) { /* Falls ausserhalb des Netzes, ... */
if (created) a.flag |= SRC_CREATED; //-2005-01-23
pp->flag = a.flag; /* ... vermerken und naechstes nehmen */
continue; }
if (INSOLID(l.i,l.j,l.k)) {
if (created) {
Try2Escape(&a, &l); eG(14); //-2004-10-10
}
else {
Try2Escape(&a, &l); eG(15); //-2004-12-09
}
}
kk = (a.z > Zref) ? -1 : created; /* Oberhalb Zref ? */
ClcLocMod(&l, &m, kk); eG(4); /* Felder interpolieren */
if (created) { /* Bei einem neuen Teilchen ... */
r1 = WrkRnd(); /* ... Geschwindigkeiten initalisieren */
r2 = WrkRnd();
r3 = WrkRnd();
vred = pr->vred; /* Reduktion der vert. Startgeschwindigk.*/
hred = pr->hred; /* Reduktion der hor. Startgeschwindigk. */
a.u = hred*m.exx*r1;
a.v = hred*(m.eyx*r1 + m.eyy*r2);
a.w = vred*(m.ezx*r1 + m.ezy*r2 + m.ezz*r3);
if (!ExactTau) m.tau *= 0.5 + WrkRnd01(); //-2001-09-23
v2 = m.vx*m.vx + m.vy*m.vy;
ClcAfu(ps, v2, Stab, &a.afuhgt, &a.afutsc); eG(5);
if (a.afutsc > 0) {
float vs;
vs = a.afuhgt/a.afutsc;
ClcAfuRxyz(ps, &vs, &a.afurx, &a.afury, &a.afurz);
a.afuhgt = vs*a.afutsc;
}
else {
a.afurx = 0; a.afury = 0; a.afurz = 1; }
}
retry = 0; /* Zaehler fuer Schrittversuche setzen */
memcpy(&e, &a, PtlRecSize); /* Alles uebernehmen (Massen g[] !) */
predicting = predict;
while (tp < T2) { /*---------- Zeitschritte durchfuehren -----------*/
tau = m.tau;
if (predicting) { //-2001-09-22
predicting = 0;
e.x = a.x + tau*m.vx; // Nur horizontale ...
e.y = a.y + tau*m.vy; // ... Advektion durchfuehren
e.z = a.z;
ll = l;
while (1) {
ClcPtlLoc(&a, &e, &ll, WRK_USE_AZ); eG(6); // Position berechnen
if (e.flag & SRC_EOUT)
break;
e.z += tau*(m.vz - vsed); // Vertikale Advektion
ClcPtlLoc(&a, &e, &ll, 0); eG(9); // Wieder Position berechnen
if (e.flag & SRC_EOUT)
break;
rc = ClcLocMod(&ll, &mm, -1); eG(7); // vx, vy, vz berechnen
if (rc <= 0) // in verbotenem Gebiet
break;
m.vx = 0.5*(m.vx + mm.vx);
m.vy = 0.5*(m.vy + mm.vy);
m.vz = 0.5*(m.vz + mm.vz);
break;
}
memcpy(&e, &a, PtlRecSize); // Position zuruecksetzen
continue;
}
steps++;
if (TrcMode) trcPtl(np,&e,tp,m);
CountStep++;
if (!ShowProgress && CountStep>=CountMax) {
CountStep = 0;
CountPhase++;
if (CountPhase >= 4) CountPhase = 0;
printf("%c\b", PhaseSymbols[CountPhase]); fflush(stdout);
}
e.refl = 0;
if (a.z > Zref) { /* Oberhalb Zref ... */
e.u = 0; e.v = 0; e.w = 0; } /* ... ohne Diffusion rechnen */
else {
r1 = WrkRnd();
r2 = WrkRnd();
r3 = WrkRnd();
e.u = m.pxx*a.u + m.pxy*a.v + m.pxz*a.w
+ m.wx + m.lxx*r1;
e.v = m.pyx*a.u + m.pyy*a.v + m.pyz*a.w
+ m.wy + m.lyx*r1 + m.lyy*r2;
e.w = m.pzx*a.u + m.pzy*a.v + m.pzz*a.w
+ m.wz + m.lzx*r1 + m.lzy*r2 + m.lzz*r3;
}
e.x = a.x + tau*m.vx; /* Horizontale ... */
e.y = a.y + tau*m.vy; /* ... Advektion durchfuehren */
e.z = a.z;
vxdrop = m.vx; //-2023-07-17
vydrop = m.vy; //-2023-07-17
dzp = tau*(m.vz + e.w - vsed); /* Vertikalen Transport merken */
if (retry > MinRetry && vsed > 0)
dzp += tau*vsed; //-2005-02-24
if (e.afuhgt > 0) { /* Abgasfahnen-Ueberhoehung */
if (e.afurz < 0. && afufac > 0.) {
/* Schwergas nach VDI 3783/1 */
if ((e.afutsc <= m.tau) || (ABS(e.afuhgt) < 0.01)) {
dafu = e.afuhgt;
dzp += dafu*e.afurz; /* Rest aufaddieren */
e.y += dafu*e.afury;
e.x += dafu*e.afurx;
e.afuhgt = 0;
}
else { /* Teil aufaddieren */
// src_rf: used to define fh*Lc
if (src_rf <= 0.) eX(56);
float red = pow(e.h/src_rf, 3.); //-2024-01-17
if (red > 1.) red = 1.;
dafu = e.afuhgt*tau/e.afutsc;
dzp += dafu*e.afurz*red; // reduced shift
e.y += dafu*e.afury;
e.x += dafu*e.afurx;
e.afuhgt -= dafu; // no reduction to guarantee user-defined decay time
}
}
else {
/* Standardverfahren */
if ((e.afutsc <= m.tau) || (ABS(e.afuhgt) < 1.00)) {
dafu = e.afuhgt;
dzp += dafu*e.afurz; /* Rest aufaddieren */
e.y += dafu*e.afury;
e.x += dafu*e.afurx;
e.afuhgt = 0;
}
else { /* Teil aufaddieren */
dafu = e.afuhgt*tau/e.afutsc;
dzp += dafu*e.afurz; /* im z-System */
e.y += dafu*e.afury;
e.x += dafu*e.afurx;
e.afuhgt -= dafu;
}
}
}
if (retry == 0) { /* Beim ersten Durchlauf ... */
ll = l; /* alte gueltige Position merken */
e.flag = SRC_EXIST; /* ... und Flags zuruecksetzen. */
ClcPtlLoc(&a, &e, &l, WRK_USE_AZ); eG(6); /* Position berechnen */
if (e.flag & SRC_EOUT) { /* Rechengebiet verlassen, */
a.flag = e.flag; /* vermerken ... */
break; } /* ... und aufhoeren */
if (INSOLID(l.i,l.j,l.k)) { /* Innerhalb eines Koerpers ? */
e.x = a.x; e.y = a.y; /* Ort zuruecksetzen */
l = ll; /* Alte Position uebernehmen */
m.vx = 0; m.vy = 0; /* Advektion ausschalten */
}
}
e.x += tau*e.u; /* turbulenter Transport */
e.y += tau*e.v;
e.z += dzp;
ClcPtlLoc(&a, &e, &l, 0); eG(9); /* Wieder Position berechnen */
if (e.flag & SRC_EOUT) { /* Rechengebiet verlassen, */
a.flag = e.flag; /* vermerken ... */
break; } /* ... und aufhoeren */
kk = (e.z > Zref) ? -1 : 0;
if (INSOLID(l.i,l.j,l.k)) rc = 0; /* Innerhalb eines Koerpers ? */
else {
rc = ClcLocMod(&l, &m, kk); eG(7); /* Felder interpolieren */
}
if (rc == 0) { /* In verbotenem Gebiet */
retry++; /* Anzahl Versuche pruefen */
if (retry > MaxRetry) {
a.flag |= SRC_REMOVE;
fprintf(MsgFile,
"*** %5d: %6.2f (%5.3lf,%5.3lf,%5.3lf) (%5.3f,%5.3f,%5.3f)"
" F(%5.3f,%5.3f,%5.3f)\n",
np, m.tau, a.x, a.y, a.z, a.u, a.v, a.w, m.vx, m.vy, m.vz);
break; }
if (ll.k > l.k) e.refl |= WRK_REFB; /* Kommt von oben : Deposition ! */
if (retry < MinRetry) { /* Bei den ersten Malen */
a.u=-a.u; a.v=-a.v; a.w=-a.w; } /* ... Teilchen umdrehen, */
else { /* sonst */
a.u = 0; a.v = 0; a.w = 0; /* Eigengeschwindigkeit loeschen */
m.wx = 0; m.wy = 0; m.wz = 0; /* Driftgeschwindigkeit loeschen */
m.vx = 0; m.vy = 0; m.vz = 0; } /* Windgeschwindigkeit loeschen */
continue; /* ... und noch einmal versuchen */
}
//-2011-10-13
if (Pxfr) { // Chemical reactions
XFRREC *prec;
int ir, ixfr;
for (icmp=0; icmp
ixfr = Ixfr[ir]; // get the first list index
if (ixfr < 0) // no reactions for this substance
continue;
gsum = a.g[2*icmp]; // get the old mass
while (ixfr < Nxfr) { // scan the reaction list
prec = Pxfr + ixfr; // get a reaction record
if (prec->ir != ir) // check the row index
break;
jcmp = prec->ic - cmpoff; // index of reaction partner
if (jcmp < 0 || jcmp >= cmpnum) eX(8);
gsum += tau*a.g[2*jcmp]*prec->val; // add contribution
ixfr++;
}
if (gsum < 0) eX(16);
e.g[2*icmp] = gsum; // ... new mass
}
}
//
if (e.refl & (WRK_REFB | WRK_REFT)) { /* Teilchen wurde reflektiert */
if (retry == 0) { /* Falls am Boden reflektiert, */
e.u = -e.u; /* alle Komponenten umkehren. */
e.v = -e.v;
e.w = -e.w /* - 2*m.wz */;
}
if (e.afuhgt > 0) {
if (afufac >= 0)
e.afuhgt *= afufac;
else {
e.afuhgt *= -afufac;
e.afurz *= -1;
}
}
if (deposition && e.refl&WRK_REFB && l.i>0 && l.j>0) {
k = 0;
pf = AryPtr(Pda, l.i, l.j, k, cmpoff); if (!pf) eX(20);
if (Pwa) {
pw = AryPtr(Pwa, l.i, l.j, cmpoff); if (!pw) eX(21);
}
if (wdep_unknown) {
ClcLocMod(&l, &m, 1); eG(7); /* ezz berechnen */
for (icmp=0; icmp
else vdep = pr[icmp].vdep;
PpmClcWdep(m.ezz, vdep+vsed, vsed, wdeps+icmp);
}
wdep_unknown = 0;
}
for (icmp=0; icmp
else wdep = (Pwa) ? pw[icmp] : pr[icmp].wdep;
if (wdep < 0) eX(55); //-2013-07-08
if (wdep > 1) wdep = 1;
dg = e.g[2*icmp]*wdep; /* abgelagerte Masse ... */
e.g[2*icmp] -= dg; /* fehlt beim Teilchen ... */
pf[icmp] += dg; } /* und liegt am Boden. */
}
}
if (l.i>0 && l.j>0 && l.k>0 && l.k<=nzdos) {
pf = AryPtr(Pda, l.i, l.j, l.k, cmpoff); if (!pf) eX(12);
for (icmp=0; icmp
}
if (washout && l.i>0 && l.j>0) { //---------------------- wet deposition
int outside=0; //-2023-07-17
float xdrop; //-2023-07-17
float ydrop; //-2023-07-17
if (WetDrift) {
// record deposited mass with drop drift -2023-07-17
float hdm; // height above ground
float xdm; // centre point x
float ydm; // centre point y
int idrop; // i-index
int jdrop; // j-index
hdm = 0.5*(a.z + e.z);
xdm = 0.5*(a.x + e.x);
ydm = 0.5*(a.y + e.y);
if (bDdmet) {
float h = 0;
GrdBottom(xdm, ydm, &h);
hdm -= h;
}
xdrop = xdm + WRK_DRIFT_VRED*(vxdrop + a.u)*(hdm/vdrop);
ydrop = ydm + WRK_DRIFT_VRED*(vydrop + a.v)*(hdm/vdrop);
idrop = 1 + floor(xdrop - Xmin)/Dd;
jdrop = 1 + floor(ydrop - Ymin)/Dd;
if (idrop
// store in an outer net
outside = 1;
for (int i=0; i
pf = wetbuf + cmpoff;
}
else {
pf = AryPtr(Pda, idrop, jdrop, -1, cmpoff); if (!pf) eX(10);
}
} // end DRIFT
else {
pf = AryPtr(Pda, l.i, l.j, -1, cmpoff); if (!pf) eX(10);
}
for (icmp=0; icmp
float fac; //-2011-10-17
fac = pr[icmp].rwsh*tau;
if (fac > 1) fac = 1; //-2011-10-17
dg = g*fac; //-2010-10-17
g -= dg;
e.g[2*icmp] = g;
pf[icmp] += dg;
}
if (outside) {
WetDepoRegister(xdrop, ydrop, wetbuf); //-2023-07-17
}
} // ---------------------------------------------- end of wet deposition
discard = 1; // Noch weiter brauchbar?
for (icmp=0; icmp
|| e.g[2*icmp] > pr[icmp].mmin*e.g[2*icmp+1])
discard = 0; // Nicht abgereicherte weiter verwenden
}
if (discard) { // Leeres Teilchen als nicht ...
a.flag |= SRC_REMOVE; // ... mehr verwendbar markieren.
Ndsc++;
break; } /* und abbrechen. */
retry = 0; /* Versuchs-Zaehler zuruecksetzen */
tp += tau; /* Zeit weitersetzen */
memcpy(&a, &e, PtlRecSize); /* Anfangswerte setzen */
a.flag = SRC_EXIST; /* Flags zuruecksetzen */
a.refl = 0;
predicting = predict;
} /*-------------------------------- Ende der Zeitschleife ---------*/
a.rnd = WrkRndGetSeed(); /* Zufallszahl merken */
a.t = tp; /* Endzeit eintragen */ //-2019-03-14
memcpy(pp, &a, PtlRecSize); /* Werte in der Tabelle eintragen */
if (TrcMode) trcPtl(np,&a,tp,m);
} /*++++++++++++++++++++++++++++++++++ Teilchen fertig +++++++++++++++*/
PtlEnd(handle);
FREE(pa);
FREE(pe);
FREE(wetbuf); //-2023-07-17
vLOG(5)("WRK:RunPtl() returning, n=%d", NumRun);
return NumRun;
eX_30:
eMSG(_no_memory_);
eX_1:
eMSG(_no_handle_);
eX_2:
if (MsgFile != NULL)
fprintf(MsgFile, _current_interval_$$_, T1, T2);
prnPtl(MsgFile, "PTL", np+1, &a, (double)T1); //-2004-01-12
eMSG(_forgotten_particle_);
eX_3: eX_6: eX_9:
prnPtl(MsgFile, "PTL", np, &a, (double)T1); //-2004-01-12
eMSG(_invalid_position_);
eX_4: eX_7:
eMSG(_cant_interpolate_);
eX_5:
eMSG(_no_plume_rise_);
eX_8:
eMSG(_no_transfer_data_);
eX_10:
eX_12:
eX_20:
eMSG(_index_$$$_, l.i, l.j, k);
eX_21:
eMSG(_index_$$_, l.i, l.j);
eX_13:
eMSG(_invalid_source_count_$_, srcind);
eX_14:
eMSG(_inside_building_$$$_, pp->x, pp->y, pp->h);
eX_15:
prnPtl(stdout, "PTL", np, &a, (double)T1);
eMSG(_moved_inside_);
eX_16:
eMSG(_chemical_reactions_too_fast_);
eX_55:
eMSG(_negative_wdep_);
eX_56:
eMSG(_invalid_rf_);
}
#undef a
#undef e
/*================================================================= SetSimParm
*/
static void show_progress( long t ) { //-2001-07-10
long p;
p = (1000.0*(t-MI.start))/(MI.end-MI.start);
if (p != Progress && myid == 0) { //bp
//bp printf(_progress_$_, 0.1*p); //!!!!!!!!!!!!!!!!!!!!!!!!!
// printf("Proc %d %s: Done: %5.1f %%\n", myid ,MsgDate(), 0.1*p); //bp
if (start_sec == 0) start_sec = time(NULL); //bp
long diff = time(NULL) - start_sec; //bp
printf("Done: %d sec %4.1f %%\n", diff, 0.1*p); //bp
fflush(stdout);
Progress = p;
}
}
static long SetSimParm(
void )
{
dP(SetSimParm);
int modlen, i, j, k;
long igp, ima, ipi, iga, isa, ida, ipp, ixa, iwa, usage, ibn;
ARYDSC *pa;
char name[256], t1s[40], t2s[40]; //-2004-11-26
TXTSTR hdr = { NULL, 0 };
strcpy(t1s, TmString(&T1));
strcpy(t2s, TmString(&T2));
vLOG(5)("WRK: SetSimParm() [%s,%s] ", t1s, t2s);
TrcFile = MsgFile;
strcpy(name, NmsName(StdIdent));
igp = IDENT(GRDpar, 0, Gl, Gi);
pa = TmnAttach(igp, NULL, NULL, 0, NULL); eG(1);
if (!pa) eX(2);
Pgp = (GRDPARM*) pa->start;
Nx = Pgp->nx;
Ny = Pgp->ny;
Nz = Pgp->nz;
if (MI.flags & FLG_GAMMA) Nzm = MIN(Pgp->nzmap, Nz);
else Nzm = MIN(Pgp->nzdos, Nz);
Dd = Pgp->dd;
Xmin = Pgp->xmin;
Gx = Nx*Dd;
Xmax = Xmin + Gx;
Ymin = Pgp->ymin;
Gy = Ny*Dd;
Ymax = Ymin + Gy;
bPerx = (0 != (Pgp->bd & GRD_PERX));
if (bPerx) {
Xbmin = Xmin;
Xbmax = Xmax; }
else {
Xbmin = Xmin - Pgp->border;
Xbmax = Xmax + Pgp->border; }
bPery = (0 != (Pgp->bd & GRD_PERY));
if (bPery) {
Ybmin = Ymin;
Ybmax = Ymax; }
else {
Ybmin = Ymin - Pgp->border;
Ybmax = Ymax + Pgp->border; }
DefMode = Pgp->defmode;
PrfMode = Pgp->prfmode;
switch (PrfMode) {
case 2: modlen = sizeof(MD2REC);
break;
case 3: modlen = sizeof(MD3REC);
break;
default: eX(20);
}
Zscl = Pgp->zscl;
//-2005-01-12
if (Zscl != 0) eX(90);
Zref = Pgp->hmax;
iga = IDENT(GRDarr, 0, 0, 0);
pa = TmnAttach(iga, NULL, NULL, 0, NULL); eG(4);
if (!pa) eX(5);
Sk = (float*) pa->start;
TmnDetach(iga, NULL, NULL, 0, NULL); eG(6);
ipi = IDENT(PRMpar, 0, 0, 0);
PrmT2 = T1;
pa = TmnAttach(ipi, &T1, &PrmT2, 0, NULL); eG(7);
if (!pa) eX(8);
if (PrmT2 < T2) T2 = PrmT2;
if (T2 <= T1) eX(41);
if (ShowProgress) show_progress(T1);
precep = BlmPprm->Precep; //-2023-07-17
vdrop = (precep <= 0) ? 0 : WRK_DRIFT_VDROP*pow(precep, WRK_DRIFT_VDROPEXP); //-2023-07-17
Ppi = (PRMINF*) pa->start;
Nc = Ppi->sumcmp;
Flags = Ppi->flags;
TrcMode = (0 != (Flags & FLG_TRACING)); //-2002-01-05 lj
OnlyAdvec = (0 != (Flags & FLG_ONLYADVEC));
NoPredict = (0 != (Flags & FLG_NOPREDICT));
ExactTau = (0 != (Flags & FLG_EXACTTAU));
WetDrift = (0 != (Flags & FLG_WETDRIFT)); //-2023-07-17
TmnDetach(ipi, NULL, NULL, 0, NULL); eG(9);
ModT2 = T1;
ima = IDENT(MODarr, 0, Gl, Gi);
Pma = TmnAttach(ima, &T1, &ModT2, 0, &hdr); eG(10); //-2011-06-29
if (!Pma) eX(11);
bDdmet= (Pma->numdm == 3);
if (bDdmet) {
AryAssert(Pma, modlen, 3, 0, Nx, 0, Ny, 0, Nz); eG(22);
DmnDatAssert("Delta|Dd|Delt", hdr.s, Dd); eG(23); //-2011-06-29
DmnDatAssert("Xmin", hdr.s, Xmin); eG(24); //-2011-06-29
DmnDatAssert("Ymin", hdr.s, Ymin); eG(25); //-2011-06-29
ibn = IDENT(GRDarr, GRD_ISLD, Gl, Gi);
Pbn = TmnAttach(ibn, NULL, NULL, TMN_MODIFY, NULL); if (!Pbn) eX(50);
for (i=1; i<=Nx; i++)
for (j=1; j<=Ny; j++)
for (k=1; k<=Nz; k++)
if (((MD3REC*)AryPtrX(Pma,i,j,k-1))->vz <= WND_VOLOUT)
*(long*)AryPtrX(Pbn,i,j,k) |= GRD_VOLOUT;
}
else {
AryAssert(Pma, modlen, 1, 0, Nz); eG(26);
Pbn = NULL;
}
TmnInfo(ima, NULL, NULL, &usage, NULL, NULL); //-2001-06-29
Valid = (0 == (usage & TMN_INVALID)); //-2001-06-29
if (!Valid && DropT1 != T1) { //-2011-07-04
InterDrop += (T2 - T1);
AverDrop += (T2 - T1);
TotalDrop += (T2 - T1);
DropT1 = T1;
}
Lm = 0;
DmnGetFloat(hdr.s, "LM", "%f", &Lm, 1);
Z0 = 0.5;
DmnGetFloat(hdr.s, "Z0", "%f", &Z0, 1);
D0 = 6*Z0;
DmnGetFloat(hdr.s, "D0", "%f", &D0, 1); //-2018-10-04
Us = 0.3;
DmnGetFloat(hdr.s, "US", "%f", &Us, 1);
DmnGetFloat(hdr.s, "KL|KM", "%f", &Stab, 1);
Hm = 0;
DmnGetFloat(hdr.s, "HM", "%f", &Hm, 1);
Ta = 10.;
DmnGetFloat(hdr.s, "TA", "%f", &Ta, 1); //-2018-10-04
Rh = 70.;
DmnGetFloat(hdr.s, "RH", "%f", &Rh, 1); //-2018-10-04
if (bDdmet) {
if ( (Pgp->sscl == Sk[Nz])
&& (Pgp->bd & GRD_REFZ)
&& (Zref == Zscl)
&& (Zscl > 0));
else Zref = 9999;
}
else {
if (Zref <= 0) {
if (Pgp->bd & GRD_REFZ) Zref = Hm;
if (Zref <= 0) Zref = 9999;
}
else {
if (Zref
}
if (Zref <= Pgp->zmax) eX(51);
}
TmnDetach(igp, NULL, NULL, 0, NULL); eG(3);
if (ModT2 < T2) T2 = ModT2;
if (T2 <= T1) eX(42);
PpmT2 = T1;
ipp = IDENT(PPMarr, 0, Gl, Gi);
Ppp = TmnAttach(ipp, &T1, &PpmT2, 0, NULL); eG(12);
if (!Ppp) eX(13);
if (PpmT2 < T2) T2 = PpmT2;
if (T2 <= T1) eX(43);
if (bDdmet) {
WdpT2 = T1;
iwa = IDENT(WDParr, 0, Gl, Gi);
Pwa = TmnAttach(iwa, &T1, &WdpT2, 0, &hdr); eG(16);
if (!Pwa) eX(17);
if (WdpT2 < T2) T2 = WdpT2;
if (T2 <= T1) eX(44);
AryAssert(Pwa, sizeof(float), 3, 1, Nx, 1, Ny, 0, Nc-1); eG(27);
DmnDatAssert("Delta|Delt|Dd", hdr.s, Dd); eG(28); //-2011-06-29
DmnDatAssert("Xmin", hdr.s, Xmin); eG(29); //-2011-06-29
DmnDatAssert("Ymin", hdr.s, Ymin); eG(30); //-2011-06-29
}
else Pwa = NULL;
if (Flags & FLG_CHEM) {
float *pf;
XfrT2 = T1;
ixa = IDENT(CHMarr, 0, 0, 0); //-2001-04-25
Pxa = TmnAttach(ixa, &T1, &XfrT2, 0, NULL); eG(31);
if (!Pxa) eX(32);
if (XfrT2 < T2) T2 = XfrT2;
if (T2 <= T1) eX(45);
AryAssert(Pxa, sizeof(float), 2, 0, Nc-1, 0, Nc-1); eG(33);
// -2011-10-13
pf = (float*)Pxa->start;
Nxfr = 0;
for (i=0; i
Nxfr++;
pf++;
}
if (Nxfr > 0) {
XFRREC *pxfr;
int nr=0;
Pxfr = ALLOC(Nxfr*sizeof(XFRREC));
Ixfr = ALLOC(Nc*sizeof(int));
pxfr = Pxfr;
pf = (float*)Pxa->start;
for (i=0; i
for (j=0; j
pxfr->ir = i;
pxfr->ic = j;
pxfr->val = *pf;
if (Ixfr[i] < 0)
Ixfr[i] = nr;
pxfr++;
nr++;
}
pf++;
}
}
}
}
else {
Pxa = NULL;
Pxfr = NULL;
Nxfr = 0;
}
SrcT2 = T1;
isa = IDENT(SRCtab, 0, 0, 0);
Psa = TmnAttach(isa, &T1, &SrcT2, 0, NULL); eG(34);
if (!Psa) eX(35);
if (SrcT2 < T2) T2 = SrcT2;
if (T2 <= T1) eX(46);
//
if (Flags & FLG_PLURIS) { //-2018-10-04
int nsrc;
SRCREC *psrc;
nsrc = Psa->bound[0].hgh + 1;
psrc = (SRCREC*)Psa->start;
for (i=0; i
vLOG(5)("WRK: ClcPlr('%s') = %d\n", psrc[i].name, rc);
}
}
//
ida = IDENT(DOSarr, Gp, Gl, Gi);
strcpy(name, NmsName(ida));
Pda = NULL;
if (TmnInfo(ida, &DosT1, &DosT2, &usage, NULL, NULL)) {
if (usage & TMN_DEFINED) {
if (DosT1>T1 || DosT2
if (!Pda) eX(39);
}
}
if (!Pda) {
Pda = TmnCreate(ida, sizeof(float),
4, 1, Nx, 1, Ny, -1, Nzm, 0, Nc-1); eG(37);
DosT1 = T1;
DosT2 = T2;
}
TxtClr(&hdr);
if (TrcMode) trcInit(Gp, Gl, Gi, T1, T2);
NumRemove = 0;
Ndsc = 0;
// -2023-07-17
if (GrdPtab != NULL) {
NETREC *pnet = (NETREC*)GrdPtab->start; // largest grid
WetXmin = pnet->xmin;
WetXmax = pnet->xmax;
WetYmin = pnet->ymin;
WetYmax = pnet->ymax;
}
else {
WetXmin = GrdPprm->xmin;
WetXmax = GrdPprm->xmax;
WetYmin = GrdPprm->ymin;
WetYmax = GrdPprm->ymax;
}
numnet = GrdPprm->numnet; //-2023-07-31
//
vLOG(5)("WRK:SetSimParm() returning");
return 0;
eX_1: eX_2: eX_3:
eMSG(_no_grid_data_);
eX_4: eX_5: eX_6:
eMSG(_no_vertical_grid_);
eX_7: eX_8: eX_9:
eMSG(_no_general_parameters_);
eX_10: eX_11:
strcpy(name, NmsName(ima));
eMSG(_no_model_field_$_, name);
eX_12: eX_13:
strcpy(name, NmsName(ipp));
eMSG(_no_particle_parameters_$_, name);
eX_16: eX_17:
strcpy(name, NmsName(iwa));
eMSG(_no_depo_prob_$_, name);
eX_20:
eMSG(_invalid_mode_$_, PrfMode);
eX_22: eX_23: eX_24: eX_25: eX_26:
strcpy(name, NmsName(ima));
eMSG(_improper_model_$_, name);
eX_28: eX_29: eX_30:
if (hdr.s) nMSG("hdr=%s<<<", hdr.s);
eX_27:
strcpy(name, NmsName(iwa));
eMSG(_improper_structure_$_, name);
eX_31: eX_32:
strcpy(name, NmsName(ixa));
eMSG(_no_transfer_$_, name);
eX_33:
strcpy(name, NmsName(iwa));
eMSG(_improper_transfer_$_, name);
eX_34: eX_35:
eMSG(_cant_get_plume_rise_);
eX_36: eX_37: eX_39:
eMSG(_no_dose_$_, name);
eX_38:
strcpy(t1s, TmString(&T1));
strcpy(t2s, TmString(&T2));
eMSG(_dose_$$$_parm_$$_, name, t1s, TmString(&DosT2), t1s, t2s);
eX_41: eX_42: eX_43: eX_44: eX_45: eX_46:
strcpy(t1s, TmString(&T1));
strcpy(t2s, TmString(&T2));
eMSG(_invalid_interval_$$_, t1s, t2s);
eX_50:
eMSG(_no_boundaries_);
eX_51:
eMSG(_reflection_$_$$$_, Zref, Gl, Gi, Pgp->zmax);
eX_90:
eMSG(_zscl_sscl_0_);
}
//================================================================ FreeSimParm
//
static long FreeSimParm( /* Simulations-Parameter freigeben. */
void )
{
dP(FreeSimParm);
long mode;
TmnDetach(IDENT(MODarr, 0, Gl, Gi), NULL, NULL, 0, NULL); eG(41);
Pma = NULL;
TmnDetach(IDENT(PPMarr, 0, Gl, Gi), NULL, NULL, 0, NULL); eG(42);
Ppp = NULL;
TmnDetach(IDENT(SRCtab, 0, 0, 0), NULL, NULL, 0, NULL); eG(43);
Psa = NULL;
if (Pwa) {
TmnDetach(IDENT(WDParr, 0, Gl, Gi), NULL, NULL, 0, NULL); eG(44);
Pwa = NULL;
}
if (Pxa) {
TmnDetach(IDENT(CHMarr, 0, 0, 0), NULL, NULL, 0, NULL); eG(45);
Pxa = NULL;
if (Pxfr) {
FREE(Pxfr); //-2011-10-13
Pxfr = NULL;
}
if (Ixfr) {
FREE(Ixfr); //-2011-10-13
Ixfr = NULL;
}
}
mode = TMN_MODIFY | ((Valid) ? TMN_SETVALID : TMN_INVALID); //-2001-06-29
TmnDetach(IDENT(DOSarr, Gp, Gl, Gi), &DosT1, &T2, mode, NULL); eG(46);
Pda = NULL;
if (Pbn) {
TmnDetach(IDENT(GRDarr, GRD_ISLD, Gl, Gi), NULL, NULL, TMN_MODIFY, NULL); eG(47);
Pbn = NULL;
}
if (NumRemove)
vLOG(4)("%d particles from inside bodies removed", NumRemove);
if (Ndsc)
vLOG(4)("%d particles discarded", Ndsc);
return 0;
eX_41: eX_42: eX_43: eX_44: eX_45: eX_46: eX_47:
eMSG(_cant_detach_);
}
//================================================================== WrkHeader
//
#define CAT(a,b) {TxtPrintf(&txt,(a),(b)); TxtCat(&hdr,txt.s);}
char *WrkHeader( /* the header (global storage) */
long id, /* identification */
long *pt1, /* start of the validity time */
long *pt2 ) /* end of validity time */
{
dP(WrkHeader);
char name[256], t1s[40], t2s[40], dts[40], rdat[40], format[80], vldf[40]; //-2008-12-04 //-2011-04-13 //-2018-10-04
int gl, gi, iga, igp, k, ict, nz, nc, seq;
char atype;
float *sk;
ARYDSC *pa, *pct, *pga;
GRDPARM *pgp;
CMPREC *pcmp;
enum DATA_TYPE dt;
TXTSTR txt = { NULL, 0 };
TXTSTR hdr = { NULL, 0 };
strcpy(name, NmsName(id));
dt = XTR_DTYPE(id);
gl = XTR_LEVEL(id);
gi = XTR_GRIDN(id);
switch (dt) {
case DOSarr: atype = 'D';
break;
case SUMarr: atype = 'D';
break;
case DOStab: atype = 'D';
break;
case CONtab: atype = 'C';
break;
case GAMarr: atype = 'G';
break;
default: atype = '?';
}
igp = IDENT(GRDpar, 0, gl, gi);
pa = TmnAttach(igp, NULL, NULL, 0, NULL); eG(1);
if (!pa) eX(2);
pgp = (GRDPARM*) pa->start;
iga = IDENT(GRDarr, 0, 0, 0);
pga = TmnAttach(iga, NULL, NULL, 0, NULL); eG(3);
if (!pga) eX(4);
sk = (float*) pga->start;
nz = pga->bound[0].hgh + 1;
ict = IDENT(CMPtab, 0, 0, 0);
pct = TmnAttach(ict, NULL, NULL, 0, NULL); eG(5);
if (!pct) eX(6);
pcmp = pct->start;
nc = pct->bound[0].hgh + 1;
strcpy(t1s, TmString(pt1));
strcpy(t2s, TmString(pt2));
strcpy(dts, TmString(&MI.cycle));
strcpy(format, "Con%10.2e");
strcpy(vldf, "V");
if (MI.flags & FLG_SCAT) {
strcat(format, "Dev%(*100)5.1f"); //-2011-06-29
strcat(vldf, "V");
}
seq = MI.cycind;
if (MI.average > 1) seq = 1 + (seq-1)/MI.average;
//
// dmna header //-2011-06-29
TxtCpy(&hdr, "\n");
// sprintf(s, "TALWRK_%d.%d.%s", StdVersion, StdRelease, StdPatch);
CAT("prgm \"%s\"\n", TalLabel); //-2011-12-05
CAT("idnt \"%s\"\n", MI.label);
CAT("artp \"%c\"\n", atype);
CAT("axes \"%s\"\n", "xyzs");
CAT("vldf \"%s\"\n", vldf);
CAT("file \"%s\"\n", name);
CAT("form \"%s\"\n", format);
CAT("t1 \"%s\"\n", t1s);
CAT("t2 \"%s\"\n", t2s);
CAT("dt \"%s\"\n", dts);
if (TmGetDate(MI.rbrk) < TM_MAX_DATE) {
TmFormatDate(rdat, "%24D", MI.rbrk);
CAT("rdat \"%s\"\n", rdat);
}
CAT("index %d\n", seq);
CAT("groups %d\n", MI.numgrp);
CAT("refx %d\n", pgp->refx);
CAT("refy %d\n", pgp->refy);
if (pgp->refx > 0 && *pgp->ggcs) //-2008-12-04
CAT("ggcs \"%s\"\n", pgp->ggcs); //-2008-12-04 //-2008-12-04
CAT("xmin %s\n", TxtForm(pgp->xmin, 6));
CAT("ymin %s\n", TxtForm(pgp->ymin, 6));
CAT("delta %s\n", TxtForm(pgp->dd, 6));
CAT("zscl %s\n", TxtForm(pgp->zscl, 6));
CAT("sscl %s\n", TxtForm(pgp->sscl, 6));
TxtCat(&hdr, "sk ");
for (k=0; k
TxtCat(&hdr, "\n");
TxtCat(&hdr, "name ");
for (k=0; k
TxtCat(&hdr, "\n");
TxtCat(&hdr, "unit ");
for (k=0; k
TxtCat(&hdr, "\n");
TxtCat(&hdr, "refc ");
for (k=0; k
TxtCat(&hdr, "\n");
TxtCat(&hdr, "refd ");
for (k=0; k
TxtCat(&hdr, "\n");
CAT("valid %1.6f\n", WrkValid(VALID_AVER)); //-2011-07-04
TmnAttach(id, NULL, NULL, TMN_MODIFY, NULL); eG(7); //-2001-06-29
TmnDetach(id, NULL, NULL, TMN_MODIFY|TMN_SETVALID, NULL); eG(13);
TmnDetach(igp, NULL, NULL, 0, NULL); eG(10);
TmnDetach(iga, NULL, NULL, 0, NULL); eG(11);
TmnDetach(ict, NULL, NULL, 0, NULL); eG(12);
TxtClr(&txt);
return hdr.s;
eX_1: eX_2: eX_3: eX_4: eX_5: eX_6: eX_7:
nMSG(_cant_attach_);
return NULL;
eX_10: eX_11: eX_12: eX_13:
nMSG(_cant_detach_);
return NULL;
}
#undef CAT
//=================================================================== WrkValid
float WrkValid( int type ) { //-2011-07-04
dP(WrkValid);
float valid;
if (type == VALID_TOTAL) {
if (TotalTime <= 0) eX(1);
valid = (TotalTime - TotalDrop)/((float)TotalTime);
}
else if (type == VALID_AVER) {
if (AverTime <= 0) eX(2);
valid = (AverTime - AverDrop)/((float)AverTime);
}
else if (type == VALID_INTER) {
if (InterTime != MI.cycle) eX(3);
valid = (InterTime - InterDrop)/((float)InterTime);
}
else eX(4);
if (valid < 0) eX(5);
if (valid > 1) eX(6);
return valid;
eX_1: eX_2: eX_3: eX_4: eX_5: eX_6:
eMSG(_error_calculating_valid_);
return 0;
}
//==================================================================== WrkInit
//
long WrkInit( /* initialize server */
long flags, /* action flags */
char *istr ) /* server options */
{
dP(WrkInit);
long iddos, mask;
char *jstr, *ps, s[200];
int vrb, dsp;
if (StdStatus & STD_INIT) return 0;
if (istr) {
jstr = istr;
ps = strstr(istr, "-v");
if (ps) sscanf(ps+2, "%d", &StdLogLevel);
ps = strstr(istr, "-p");
if (ps) sscanf(ps+2, "%d", &ShowProgress);
ps = strstr(istr, "-y");
if (ps) sscanf(ps+2, "%d", &StdDspLevel);
ps = strstr(istr, "-d");
if (ps) strcpy(AltName, ps+2);
}
else jstr = "";
vLOG(3)("WRK_%d.%d.%s (%08lx,%s)", StdVersion, StdRelease, StdPatch, flags, jstr);
StdStatus |= flags;
iddos = IDENT(DOSarr, 0, 0, 0);
mask = ~(NMS_GROUP|NMS_GRIDN|NMS_LEVEL);
TmnCreator(iddos, mask, TMN_UNIQUE, "", WrkServer, WrkHeader); eG(1);
vrb = StdLogLevel;
dsp = StdDspLevel;
sprintf(s, " GRD -v%d -y%d -d%s", vrb, dsp, AltName);
GrdInit(flags, s); eG(2);
sprintf(s, " PRM -v%d -y%d -d%s", vrb, dsp, AltName);
PrmInit(flags, s); eG(3);
sprintf(s, " PTL -v%d -y%d -d%s", vrb, dsp, AltName);
PtlInit(flags, s); eG(4);
sprintf(s, " SRC -v%d -y%d -d%s", vrb, dsp, AltName);
SrcInit(flags, s); eG(5);
sprintf(s, " MOD -v%d -y%d -d%s", vrb, dsp, AltName);
ModInit(flags, s); eG(6);
sprintf(s, " PPM -v%d -y%d -d%s", vrb, dsp, AltName);
PpmInit(flags, s); eG(7);
if (MI.flags & FLG_PLURIS) {
vLOG(2)("PLURIS %s (bf=%1.1f,stdw=%d)", str_version(), PlrBreakFactor, PlrStDownwash);
}
StdStatus |= STD_INIT;
return 0;
eX_1:
eMSG(_creator_$_, NmsName(iddos));
eX_2: eX_3: eX_4: eX_5: eX_6: eX_7:
eMSG(_cant_init_servers_);
}
//================================================================= WrkServer
//
long WrkServer( /* server routine for DOSarr */
char *s ) /* calling option */
{
dP(WrkServer);
char t1s[40], t2s[40], name[256]; //-2004-11-26
int ida, ipa, mask;
if (StdArg(s)) return 0;
if (*s) {
switch (s[1]) {
case 'd': strcpy(AltName, s+2);
break;
case 'J': sscanf(s+2, "%d", &TrcMode);
break;
case 'n': if (s[2]) { //-2011-07-04
sscanf(s+2, "%d,%d", &AverTime, &AverDrop);
if (AverTime < 0) {
AverTime = 0;
AverDrop = 0;
TotalTime = 0;
TotalDrop = 0;
}
}
else {
TotalTime += MI.cycle;
AverTime += MI.cycle;
}
DropT1 = -1;
InterTime = MI.cycle;
InterDrop = 0;
break;
case 'o': if (strstr(s+2, "nostdw"))
PlrStDownwash = 0;
else if (strstr(s+2, "logpluris"))
LogPluris = 1;
break;
default: eX(17);
}
return 0;
}
if (!StdIdent) eX(20);
if ((StdStatus & STD_INIT) == 0) {
WrkInit(0, ""); eG(20);
}
Gl = XTR_LEVEL(StdIdent);
Gi = XTR_GRIDN(StdIdent);
Gp = XTR_GROUP(StdIdent);
if (StdStatus & STD_TIME) T1 = StdTime;
else eX(1);
strcpy(t1s, TmString(&T1));
ida = IDENT(DOSarr, 0, 0, 0);
mask = ~(NMS_GROUP|NMS_GRIDN|NMS_LEVEL);
TmnCreator(ida, mask, 0, "", NULL, NULL); eG(2);
PtlT2 = T1;
ipa = IDENT(PTLarr, Gp, Gl, Gi);
Ppa = TmnAttach(ipa, &T1, &PtlT2, TMN_MODIFY, NULL); eG(10);
if (!Ppa) eX(11);
T2 = PtlT2; if (T2 <= T1) eX(3);
Nptl = PtlCount(ipa, SRC_EXIST, Ppa); eG(12);
SetSimParm(); eG(13);
strcpy(t2s, TmString(&T2));
strcpy(name, NmsName(StdIdent));
vDSP(2)(_moving_$$$_, name, t1s, t2s);
RunPtl(); eG(14);
vDSP(2)("");
vLOG(4)("WRK: %d particles moved for %s [%s,%s]",
NumRun, name, t1s, t2s);
FreeSimParm(); eG(15);
TmnDetach(ipa, &T1, &PtlT2, TMN_MODIFY, NULL); eG(16);
TmnCreator(ida, mask, TMN_UNIQUE, "", WrkServer, WrkHeader); eG(21);
return 0;
eX_20:
eMSG(_missing_data_);
eX_1:
eMSG(_missing_time_);
eX_2: eX_21:
strcpy(name, NmsName(ida));
eMSG(_cant_change_creator_$_, name);
eX_3:
strcpy(name, NmsName(ipa));
strcpy(t2s, TmString(&T2));
eMSG(_invalid_times_$$$_, name, t1s, t2s);
eX_10: eX_11: eX_12:
strcpy(name, NmsName(ipa));
eMSG(_no_particles_$_, name);
eX_13: eX_15:
eMSG(_no_parameters_);
eX_14:
strcpy(name, NmsName(ida));
strcpy(t2s, TmString(&T2));
eMSG(_cant_build_$$$_, name, t1s, t2s);
eX_16:
strcpy(name, NmsName(ipa));
strcpy(t1s, TmString(&T2));
strcpy(t2s, TmString(&PtlT2));
eMSG(_cant_update_$$$_, name, t1s, t2s);
eX_17:
eMSG("unknown option %s!", s+1);
}
/*==============================================================================
*
* history:
*
* 2002-06-21 lj 0.13.0 final test version
* 2002-09-24 lj 1.0.0 final release candidate
* 2002-12-10 lj 1.0.4 source parameter Tt
* 2003-02-21 lj 1.1.2 internal random number generator
* 2003-06-16 lj 1.1.7 variabel vsed
* 2004-01-12 lj 1.1.14 error message more precise
* 2004-10-01 uj 2.1.0 version number upgraded
* 2004-10-10 uj 2.1.1 Try2Escape
* 2004-11-26 lj 2.1.7 string length for names = 256
* 2004-12-09 uj 2.1.8 shift particle outside body also if not created
* 2005-01-12 uj 2.1.9 check if Zscl = 0
* 2005-01-23 uj 2.1.12 re-mark particle created outside an inner grid
* 2005-02-24 uj 2.1.15 disable sedimentation for retry > MinRetry
* 2005-03-17 uj 2.2.0 version number upgrade
* 2006-10-26 lj 2.3.0 external strings
* 2008-08-04 lj 2.4.3 writes "valid" using 6 decimals
* 2008-12-04 lj 2.4.5 writes "refx", "refy", "ggcs", "rdat"
* 2011-04-13 uj 2.4.9 handle title with 256 chars in header
* 2011-06-29 uj 2.5.0 adjusted to dmna output
* 2011-07-04 uj handling of "valid" extended
* 2011-10-13 uj 2.5.2 optimization for chemical reactions
* 2011-11-21 lj 2.6.0 wet deposition
* 2013-07-08 uj 2.6.8 check for wdep smaller 0
* 2018-10-04 uj 3.0.0 pluris
* 2019-03-14 uj 3.0.3 end time of a particle as double
* 2021-06-23 uj 3.1.1 modification dense gas coupling for VDI 3783/1
* 2021-07-03 uj linux adjustment
* 2022-10-04 uj 3.1.5 ClcPlr: check if source is outside grid(s)
* 2023-07-17 uj 3.2.0 drop drift for wet deposition
* 2023-07-31 uj 3.2.1 drop drift corrections
* 2024-01-17 uj 3.3.0 BESMAX consistent calculation
* 2024-01-17 updated for heavy gas VDI 3783/1 E
*
*============================================================================*/