NEW

AUSTAL2000 (modifiziert)

TalMtm.c

// ======================================================== TalMtm.c

//

// Calculation of daily averages for AUSTAL2000

// ============================================

//

// Copyright (C) Umweltbundesamt, 14191 Berlin, Germany, 2002-2007

// Copyright (C) Janicke Consulting, 26427 Dunum, Germany, 2002-2007

// email: info@austal2000.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:  2007-02-03 uj

//

//========================================================================

//

//  Berechnung des maximalen Tagesmittels aus einzelnen Tagesmittelwerten

//  Nur Tagesmittel aus mindestens 12 Stundenmitteln werden verwendet

//

//=========================================================================



static char *TalMtmVersion = "2.3.5";

static char *eMODn = "TalMtm";

static int CHECK = 0;



#include 

#include 

#include 

#include 

#include 

#include 



#include "IBJmsg.h"

#include "IBJary.h"

#include "IBJtxt.h"

#include "IBJdmn.h"

#include "TalUtl.h"

#include "TalInp.h"

#include "TalCfg.h"    //-2006-10-31

#include "TalMtm.h"

#include "TalMtm.nls"



#define MAXCMP 30



static ARYDSC CurDos, VecDos, VecDay, MaxCon, MaxDev, MaxDay, CncDsc, DevDsc;

static ARYDSC SumDos, SumCon, SumDev, SumDep, SumDps, SupCon, SupDev, SupDay;

static long Folge=1;

static long Anzahl, Dt, Ia, ElmSz, Nsum, Nval;

static char Eingabe[256];

static char DosName[80] = "d%04lda00.arr";

static char GrdName[80] = "";

static char Path[256];

static int Nx, Ny, Nz, Nc, Nh, Ns, Ng;

static int Maxima, Deposition, Concentration, AverageOnly;

static int Scatter=0, WriteHeader=1, Kref=-999, Kmax=-999;

static TXTSTR UsrHeader = { NULL, 0 };

static TXTSTR SysHeader = { NULL, 0 };

static int CmpIndex[8], CmpIsGas[8];

static int CmpDayNum, CmpHourNum;

static char CmpName[8], CmpFamily[8];

static char CmpConForm[32], CmpDayForm[32], CmpDepForm[32], CmpHourForm[32];

static char CmpConUnit[32], CmpDepUnit[32];

static float CmpConFac, CmpDepFac;

static float CmpConRef, CmpDayRef, CmpDepRef, CmpHourRef;

static char **Names;

static float Sk[101], Dd;

static double T1, T2;

static float SumValid, MinValid=0.5, RefDays=365;

static char locale[256]="C", locl[256]="";



//================================================================= WriteDmna

//

static int WriteDmna( int day, float valid ) {                      //-2002-08-09

  dP(WriteDmna);

  int i, j, k, k1, k2, l, ii, n=0, odor;

  float fc, *pcurdos, *pcnc, *pdev;

  float dos, sqd, dev, cnc;

  char *pc, name[40], fname[256], gname[256], srefv[40], t[80], s[256];

  char src, dst, ot[40], of[40];

  float bs = TipOdorThreshold;

  if (CurDos.numdm != 4)                                            eX(10);

  if (valid <= 0)  valid = 1;

  k2 = CurDos.bound[2].hgh;

  if (Kref >= 0) {

    if (Kref > k2)                                                   eX(1);

    k1 = Kref;

    k2 = Kref;

  }

  else if (Kmax > 0) {

    if (Kmax > k2)                                                   eX(1);

    k1 = 1;

    k2 = Kmax;

  }

  else return 0;

  //

  strcpy(name, CmpFamily);                                    

  odor = !strncmp(name, "odor",4); //bp                              //-2004-06-10

  if (odor) {

    sprintf(ot, "\"F%4.2f\"", bs);

    sprintf(of, "\"frq%s", CmpConForm+4);

  }

  //

  if (k2 > 0) {

    AryCreate(&CncDsc, sizeof(float), 3, 1,Nx, 1,Ny, k1,k2);         eG(11);

    if (Scatter) {

      AryCreate(&DevDsc, sizeof(float), 3, 1,Nx, 1,Ny, k1,k2);       eG(12);

    }

  }

  else {

    AryCreate(&CncDsc, sizeof(float), 2, 1,Nx, 1,Ny);                eG(13);

    if (Scatter) {

      AryCreate(&DevDsc, sizeof(float), 2, 1,Nx, 1,Ny);              eG(14);

    }

  }

  for (k=k1; k<=k2; k++) {

    if (k > 0) fc = CmpConFac/((Sk[k]-Sk[k-1])*Dd*Dd*Dt);

    else       fc = CmpDepFac/(Dd*Dd*Dt);

    for (i=1; i<=Nx; i++) {

      for (j=1; j<=Ny; j++) {

        l = CmpIndex[0];

        if (l < 0)  continue;

        dos = 0;

        sqd = 0;

        for (ii=0; ii<8; ii++) {

          l = CmpIndex[ii];

          if (l < 0)  break;

          if (k>0 && !CmpIsGas[ii])  continue;

          pcurdos = AryPtr(&CurDos, i, j, k, l);  if (!pcurdos)     eX(15);

          if (Scatter)  sqd += pcurdos[1] + 2*dos*pcurdos[0]/Ng;

          dos += pcurdos[0];

        }

        cnc = fc*dos/valid;                                         //-2002-08-09

        pcnc = AryPtr(&CncDsc, i, j, k);  if (!pcnc)                eX(16);

        *pcnc = cnc;

        if (Scatter) {

          dev = (dos > 0) ? (Ng*sqd/(dos*dos)-1)/(Ng-1) : 0;

          pdev = AryPtr(&DevDsc, i, j, k);  if (!pdev)              eX(17);

          *pdev = (dev > 0) ? sqrt(dev) : 0; 

          if (odor)  *pdev *= *pcnc;                                //-2004-06-10:

        }

      } // for j

    } // for i

  } // for k

  DmnRplName(&UsrHeader, "FORMAT", "form");

  DmnRplName(&UsrHeader, "VALDEF", "vldf");

  DmnRplName(&UsrHeader, "ARRTYPE", "artp");

  sprintf(t, "\"%s\"", name);

  for (pc=t; (*pc); pc++)  *pc = toupper(*pc);

  DmnRplValues(&UsrHeader, "NAME", t);

  DmnRplValues(&UsrHeader, "vldf", "\"V\"");

  DmnRplValues(&UsrHeader, "axes", "\"xyz\"");

  DmnRplValues(&UsrHeader, "sequ", "k+,j-,i+");

  DmnRplValues(&UsrHeader, "refc", NULL);

  DmnRplValues(&UsrHeader, "refd", NULL);

  DmnRplValues(&UsrHeader, "einheit", NULL);

  //

  if (*locl) {

    char s[256];

    sprintf(s, "\"%s\"", locl);

    DmnRplValues(&UsrHeader, "locl", s);

  }

  else DmnRplValues(&UsrHeader, "locl", NULL);

  if (!strcmp(locl, "german"))  { src = '.';  dst = ','; }

  else                          { src = ',';  dst = '.'; }

  DmnRplChar(&UsrHeader, "xmin", src, dst); 

  DmnRplChar(&UsrHeader, "ymin", src, dst); 

  DmnRplChar(&UsrHeader, "valid", src, dst); 

  DmnRplChar(&UsrHeader, "sk|zk", src, dst); 

  DmnRplChar(&UsrHeader, "dd|delta", src, dst); 

  //

  if (k2 > 0) {

    sprintf(gname, "%s-%03d%s%s", name, day, CfgAddString, GrdName); //-2006-10-31

    DmnRplValues(&UsrHeader, "file", TxtQuote(gname));

    sprintf(fname, "%s%s-%03d%s%s", Path, name, day, CfgAddString, GrdName); //-2006-10-31

    DmnRplValues(&UsrHeader, "form", (odor) ? of : CmpConForm);

    DmnRplValues(&UsrHeader, "unit", CmpConUnit);

    sprintf(srefv, "%10.3e", CmpConRef);

    DmnRplValues(&UsrHeader, "refv", srefv);

    DmnRplValues(&UsrHeader, "artp", (odor) ? ot : "\"C\"");

    DmnWrite(fname, UsrHeader.s, &CncDsc);                          eG(18);

    vMsg(_file_$_written_, fname);

    if (Scatter) {

      sprintf(fname, "\"%s-%03d%s%s\"", name, day, CfgDevString, GrdName);  //-2006-10-31

      DmnRplValues(&UsrHeader, "file", fname);

      sprintf(fname, "%s%s-%03d%s%s", Path, name, day, CfgDevString, GrdName); //-2006-10-31

      if (odor) {                                                   //-2004-06-10

        DmnRplValues(&UsrHeader, "form", "\"sct%5.3f\"");

        DmnRplValues(&UsrHeader, "unit", CmpConUnit);

        DmnRplValues(&UsrHeader, "refv", "1.0");

        DmnRplValues(&UsrHeader, "artp", "\"E\"");

      }

      else {

        DmnRplValues(&UsrHeader, "form", "\"dev%(*100)5.1f\"");

        DmnRplValues(&UsrHeader, "unit", NULL);

        DmnRplValues(&UsrHeader, "refv", "0.03");

        DmnRplValues(&UsrHeader, "artp", "\"E\"");

      }

      DmnWrite(fname, UsrHeader.s, &DevDsc);                        eG(19);

      vMsg(_file_$_written_, fname);

    }

  }

  else {

    DmnRplValues(&UsrHeader, "sequ", "j-,i+");

    DmnRplValues(&UsrHeader, "axes", "\"xy\"");

    sprintf(gname, "%s-dep%03d%s%s", name, day, CfgAddString, GrdName); //-2006-10-31

    DmnRplValues(&UsrHeader, "file", TxtQuote(gname));

    sprintf(fname, "%s%s-dep%03d%s%s", Path, name, day, CfgAddString, GrdName);  //-2006-10-31

    DmnRplValues(&UsrHeader, "form", CmpDepForm);

    DmnRplValues(&UsrHeader, "unit", CmpDepUnit);

    sprintf(srefv, "%10.3e", CmpDepRef);

    DmnRplValues(&UsrHeader, "refv", srefv);

    DmnRplValues(&UsrHeader, "artp", "\"XD\"");                     //-2002-07-24

    DmnWrite(fname, UsrHeader.s, &CncDsc);                          eG(20);

    vMsg("TMT: Datei %s ausgeschrieben.", fname);

    if (Scatter) {

      sprintf(fname, "\"%s-dep%03d%s%s\"", name, day, CfgDevString, GrdName); //-2006-10-31

      DmnRplValues(&UsrHeader, "file", fname);

      sprintf(fname, "%s%s-dep%03d%s%s", Path, name, day, CfgDevString, GrdName); //-2006-10-31

      DmnRplValues(&UsrHeader, "form", "\"dev%(*100)5.1f\"");

      DmnRplValues(&UsrHeader, "unit", NULL);

      DmnRplValues(&UsrHeader, "refv", "0.03");

      DmnRplValues(&UsrHeader, "artp", "\"E\"");

      DmnWrite(fname, UsrHeader.s, &DevDsc);                        eG(21);

      vMsg(_file_$_written_, fname);

    }

  }

  AryFree(&CncDsc);

  AryFree(&DevDsc);

  return 0;

eX_1:

  eMSG(_invalid_kref_);

eX_10: eX_11: eX_12: eX_13: eX_14: eX_15: eX_16: eX_17: eX_18: eX_19:

eX_20: eX_21:

  eMSG(_internal_error_);

}



//================================================================== ReadArray

//

static long ReadArray( void )

  {

  dP(ReadArray);

  long n;

  char name[256];

  if (CurDos.start)  AryFree(&CurDos);

  n = Folge + Ia;

  sprintf(name, DosName, n);

  sprintf(Eingabe, "%s%s", Path, name);

#ifndef MAIN

  MsgQuiet = 1;

#endif

  DmnRead(Eingabe, &UsrHeader, &SysHeader, &CurDos);

  MsgQuiet = 0;

  if (MsgCode < 0) {

    MsgCode = 0;

    printf("   \r");  fflush(stdout);

    return 0;

  }

  printf("%3d\r", n);  fflush(stdout);

  if (!ElmSz) {

    ElmSz = CurDos.elmsz;

    if (ElmSz > sizeof(float)) {

      Scatter = 1;

      if (ElmSz != 2*sizeof(float))                                     eX(2);

    }

    else  Scatter = 0;                                                  //-2002-08-09

  }

  else if (ElmSz > CurDos.elmsz)                                        eX(1);

  return 1;

eX_1: eX_2:

  eMSG(_inconsistent_data_);

}



/*================================================================ CompareArray

*/

static long CompareArray( void )

  {

  dP(CompareArray);

  int n, i, ii, i1, i2, j, j1, j2, k, k1, k2, l, l1, l2, m, dm;

  int nh, ni, ic;

  float *pcurdos, *pvecdos, *psumdos, *pfs, *pfd, dos, sqd;

  short day, *pvecday, *pis, *pid;

  char t1s[40], t2s[40], *_pt1s=NULL, *_pt2s=NULL;

  char name[40];

  double t1, t2;

  float valid = 1;

  int not_used = 0;

  DmnGetFloat(UsrHeader.s, "Valid", "%f", &valid, 1);

  if (valid < 0) valid = 0;                                        //-2006-07-13

  if (valid > 1) valid = 1;                                        //-2006-07-13

  if (valid < MinValid)  not_used = 1;                             //-2001-06-29

  else  SumValid += valid;        //-2002-08-09   //-2002-05-24    //-2001-06-29

  *t1s = 0;

  *t2s = 0;

  if (1 != DmnGetString(UsrHeader.s, "T1", &_pt1s, 1))                  eX(7);

  strcpy(t1s, _pt1s);  FREE(_pt1s);

  if (1 != DmnGetString(UsrHeader.s, "T2", &_pt2s, 1))                  eX(8);

  strcpy(t2s, _pt2s);  FREE(_pt2s);

  t1 = MsgDateVal(t1s);

  t2 = MsgDateVal(t2s);

  if (CurDos.numdm != 4)                                                eX(1);

  i1 = CurDos.bound[0].low;  i2 = CurDos.bound[0].hgh;

  j1 = CurDos.bound[1].low;  j2 = CurDos.bound[1].hgh;

  k1 = CurDos.bound[2].low;  k2 = CurDos.bound[2].hgh;

  l1 = CurDos.bound[3].low;  l2 = CurDos.bound[3].hgh;

  if (i1!=1 || j1!=1 || k1>0 || l1!=0 )                                 eX(2);

  dm = 1 + Scatter;

  k1 = 1 - Deposition;

  if (SumDos.start == NULL) {

    Nsum = 0;

    Nval = 0;

    T1 = t1;

    T2 = t2;

    Nx = i2;

    Ny = j2;

    Nz = k2;                                                    //-2002-03-25

    Nc = l2 + 1;

    Dt = MsgDateSeconds(t1, t2) + 0.5;                          //-2005-08-06

    if (Dt <= 0)                                                        eX(3);

    Names = ALLOC((Nc+1)*sizeof(char*));

    n = DmnGetString(UsrHeader.s, "name", Names, Nc);  if (n != Nc)     eX(9);

    n = DmnGetFloat(UsrHeader.s, "delta", "%f", &Dd, 1);  if (n != 1)   eX(14);

    n = DmnGetFloat(UsrHeader.s, "sk", "%f", Sk, 101);  if (n < Nz+1)   eX(15);

    n = DmnGetInt(UsrHeader.s, "gruppen", "%d", &Ng, 1); if (n != 1)    eX(30); //-2006-07-13

    //

    // get the indices for the requested species

    //

    strcpy(name, "gas.");

    strcat(name, CmpName);

    for (ii=0; ii<8; ii++) {

      CmpIndex[ii] = -1;

      CmpIsGas[ii] = 0;

    }

    ii = 0;

    for (ic=0; ic<6; ic++) {

      sprintf(name, "%s.%s%s", TipGruppen[ic], CmpName, TipGrpXten[ic]);

      for (l=0; l
        if (!strcmp(Names[l], name)) {

          CmpIndex[ii] = l;

          CmpIsGas[ii] = (ic < 3);

          ii++;

          break;

        }

      }

    }

    if (ii == 0)                                                    eX(16);

    for (n=0; ; n++) if (!Names[n])  break;  else FREE(Names[n]);   //-2001-11-02

    FREE(Names);  Names = NULL;

    Nh = CmpDayNum;

    if (Anzahl < 2)  Maxima = 0;

    //

    // create the arrays

    //

    if (Ng <= 4) {

      Scatter = 0;

      ElmSz = 4;

      dm = 1;                                                       //-2006-07-13

    }

    AryCreate(&SumDos, ElmSz, 3, 1,Nx, 1,Ny, k1,Nz);                eG(5);

    if (Maxima) {

      AryCreate(&VecDos, ElmSz, 4, 1,Nx, 1,Ny, 1,Nz, 0,Nh);         eG(4);

      AryCreate(&VecDay, sizeof(short), 4, 1,Nx, 1,Ny, 1,Nz, 0,Nh); eG(6);

    }

  }

  if (i2!=Nx || j2!=Ny || k2
  if ((long)(MsgDateSeconds(t1, t2) + 0.5) != Dt)                   eX(13); //-2005-08-06

  T2 = t2;

  day = Folge + Ia - 1;     // first value of Ia is 1 at this point

  if (Kref>=0 || Kmax>=0) {                                        //-2002-04-09

    WriteDmna(day, valid);                                         eG(28);  //-2002-08-09

  }

  for (k=k1; k<=Nz; k++) {

    for (i=1; i<=Nx; i++) {

      for (j=1; j<=Ny; j++) {

        l = CmpIndex[0];

        if (l < 0)  continue;

        nh = CmpDayNum;

        dos = 0;

        sqd = 0;

        for (ii=0; ii<8; ii++) {

          l = CmpIndex[ii];

          if (l < 0)  break;

          if (k>0 && !CmpIsGas[ii])  continue;

          pcurdos = AryPtr(&CurDos, i, j, k, l);  if (!pcurdos)     eX(23);

          if (Scatter)  sqd += pcurdos[1] + 2*dos*pcurdos[0]/Ng;

          dos += pcurdos[0];

        }

        if (!not_used) {                                    //-2002-05-24

          psumdos = AryPtr(&SumDos, i, j, k);  if (!psumdos)       eX(27);

          if (Scatter)  psumdos[1] += sqd + 2*dos*psumdos[0]/Ng;

          psumdos[0] += dos;

        }

        if (k>0 && Maxima && !not_used) {

          dos /= valid;                                     //-2002-08-09

          sqd /= valid*valid;                               //-2002-08-09

          pvecdos = AryPtr(&VecDos, i, j, k, 0);  if (!pvecdos)  eX(24);

          pvecday = AryPtr(&VecDay, i, j, k, 0);  if (!pvecday)  eX(25);

          for (pfs=pvecdos, m=0; m<=nh; m++, pfs+=dm)

            if (dos >= *pfs)  break;

          if (m <= nh) {

            pfd = pvecdos + dm*nh + Scatter;

            pfs = pfd - dm;

            pid = pvecday + nh;

            pis = pid - 1;

            ni = nh - m;

            for (ii=ni; ii>0; ii--) {

              *pfd-- = *pfs--;

              if (Scatter)  *pfd-- = *pfs--;

              *pid-- = *pis--;

            }

            pfd = pvecdos + dm*m;

            pid = pvecday + m;

            pfd[0] = dos;

            if (Scatter)  pfd[1] = sqd;

            pid[0] = day;                                         //-2001-06-09

          }  // if (m <= nh)

        }  // if (k > 0)

      }  // for j

    }  // for i

  }  // for k

  Nsum++;

  if (!not_used)  Nval++;

  return not_used;

eX_1:  eX_2:  eX_3:  eX_7: eX_8: eX_9: eX_11: eX_13: eX_14: eX_15:

  eMSG(_invalid_array_);

eX_4:  eX_5:  eX_6:

  eMSG(_cant_create_array_);

eX_16:

  eMSG(_no_results_for_$_, CmpName);

eX_23: eX_24: eX_25: eX_27: eX_28:

  eMSG(_internal_error_);

eX_30:

  eMSG(_gruppen_not_found_);

  }



/*================================================================= WriteResult

*/

static long WriteResult( void )

  {

  dP(WriteResult);

  char fname[256], gname[256], t[80], name[40], srefv[40], ot[40], of[40], oe[40];

  char *pc, src, dst;

  int i, j, k, k1, nh, odor;

  float *pvecdos, *pmaxcon, *pmaxdev, *psumdos, dos, sqd, dev, con;

  float *psupcon, *psupdev;

  float *psumcon, *psumdev, *psumdep, *psumdps;

  float xdep, xcon, fd, fc, valid;

  float bs = TipOdorThreshold;

  short *pvecday, *pmaxday, *psupday;

  sprintf(t, "%d", Folge+Ia-1);

  DmnRplValues(&UsrHeader, "FOLGE", t);

  nh = CmpDayNum;

  strcpy(name, CmpFamily);                                    //-2001-11-02

  odor = !strncmp(name, "odor",4); //bp                              //-2004-06-10

  if (odor) {

    sprintf(ot, "\"FHP%4.2f\"", bs);                          //-2004-12-09

    sprintf(oe, "\"EHP\"");                                   //-2004-12-09

    sprintf(of, "\"frq%s", CmpConForm+4);

  }

  xcon = CmpConFac;

  xdep = CmpDepFac;

  valid = (SumValid > 0) ? SumValid/Nsum : 1;                 //-2001-06-29

  if (Nval < 0.9*RefDays) {                                   //-2002-08-10

    nh = nh*Nval/RefDays;                                     //-2002-08-10

    if (Maxima) vMsg(_exceedings_changed_$_, nh);

  }

  k1 = (Deposition) ? 0 : 1;

  if (Concentration) {

    AryCreate(&SumCon, sizeof(float), 3, 1,Nx, 1,Ny, 1,Nz);   eG(20);

    if (Scatter) {

      AryCreate(&SumDev, sizeof(float), 3, 1,Nx, 1,Ny, 1,Nz); eG(21);

    }

  }

  if (Deposition) {

    AryCreate(&SumDep, sizeof(float), 2, 1,Nx, 1,Ny);         eG(22);

    if (Scatter) {

      AryCreate(&SumDps, sizeof(float), 2, 1,Nx, 1,Ny);       eG(23);

    }

  }

  if (Maxima) {

    AryCreate(&MaxCon, sizeof(float), 3, 1,Nx, 1,Ny, 1,Nz);   eG(4);

    AryCreate(&SupCon, sizeof(float), 3, 1,Nx, 1,Ny, 1,Nz);   eG(4);

    if (Scatter) {

      AryCreate(&SupDev, sizeof(float), 3, 1,Nx, 1,Ny, 1,Nz); eG(5);

      AryCreate(&MaxDev, sizeof(float), 3, 1,Nx, 1,Ny, 1,Nz); eG(5);

    }

    AryCreate(&MaxDay, sizeof(short), 3, 1,Nx, 1,Ny, 1,Nz);   eG(6);

    AryCreate(&SupDay, sizeof(short), 3, 1,Nx, 1,Ny, 1,Nz);   eG(6);

  }

  fd = xdep/(Dd*Dd*Dt);

  for (k=k1; k<=Nz; k++) {

    fc = xcon/((Sk[k]-Sk[k-1])*Dd*Dd*Dt);

    for (i=1; i<=Nx; i++) {

      for (j=1; j<=Ny; j++) {

        psumdos = AryPtr(&SumDos, i, j, k);  if (!psumdos)          eX(15);

        dos = psumdos[0];

        sqd = (Scatter) ? psumdos[1] : 0;

        if (k == 0) {

          psumdep = AryPtr(&SumDep, i, j);  if (!psumdep)           eX(16);

          *psumdep = fd*dos/(valid*Nsum);                           //-2001-06-29

          if (Scatter) {

            psumdps = AryPtr(&SumDps, i, j);  if (!psumdps)         eX(17);

            dev = (dos > 0) ? (Ng*sqd/(dos*dos)-1)/(Ng-1) : 0;

            *psumdps = (dev > 0) ? sqrt(dev) : 0;

          }

        }

        else {

          if (Concentration) {

            psumcon = AryPtr(&SumCon, i, j, k);  if (!psumcon)      eX(18);

            *psumcon = fc*dos/(valid*Nsum);                         //-2001-06-29

            if (Scatter) {

              psumdev = AryPtr(&SumDev, i, j, k);  if (!psumdev)    eX(19);

              dev = (dos > 0) ? (Ng*sqd/(dos*dos)-1)/(Ng-1) : 0;

              *psumdev = (dev > 0) ? sqrt(dev) : 0;

              if (odor)  *psumdev *= *psumcon;                      //-2004-06-10:

            }

          }

          if (Maxima) {

            pmaxcon = AryPtr(&MaxCon, i, j, k);  if (!pmaxcon)        eX(10);

            pmaxday = AryPtr(&MaxDay, i, j, k);  if (!pmaxday)        eX(11);

            pvecdos = AryPtr(&VecDos, i, j, k, nh);  if (!pvecdos)    eX(12);

            pvecday = AryPtr(&VecDay, i, j, k, nh);  if (!pvecday)    eX(13);

            dos = pvecdos[0];

            con = fc*dos;

            *pmaxcon = con;

            if (Scatter) {

              pmaxdev = AryPtr(&MaxDev, i, j, k);  if (!pmaxdev)      eX(14);

              sqd = pvecdos[1];

              dev = (dos > 0) ? (Ng*sqd/(dos*dos)-1)/(Ng-1) : 0;

              *pmaxdev = (dev > 0) ? sqrt(dev) : 0;

             }

            *pmaxday = *pvecday;

            if (nh >= 0) {                                              //-2002-09-28

              psupcon = AryPtr(&SupCon, i, j, k);  if (!psupcon)        eX(10);

              psupday = AryPtr(&SupDay, i, j, k);  if (!psupday)        eX(11);

              pvecdos = AryPtr(&VecDos, i, j, k, 0);  if (!pvecdos)     eX(12);

              pvecday = AryPtr(&VecDay, i, j, k, 0);  if (!pvecday)     eX(13);

              dos = pvecdos[0];

              con = fc*dos;

              *psupcon = con;

              if (Scatter) {

                psupdev = AryPtr(&SupDev, i, j, k);  if (!psupdev)      eX(14);

                sqd = pvecdos[1];

                dev = (dos > 0) ? (Ng*sqd/(dos*dos)-1)/(Ng-1) : 0;

                *psupdev = (dev > 0) ? sqrt(dev) : 0;

               }

              *psupday = *pvecday;

            } // if nh

          } // if maxima

        } // k>0

      } // for j

    } // for i

  } // for k

  DmnRplName(&UsrHeader, "FORMAT", "form");

  DmnRplName(&UsrHeader, "VALDEF", "vldf");

  DmnRplName(&UsrHeader, "ARRTYPE", "artp");

  sprintf(t, "\"%s\"", name);

  for (pc=t; (*pc); pc++)  *pc = toupper(*pc);

  DmnRplValues(&UsrHeader, "NAME", t);

  DmnRplValues(&UsrHeader, "vldf", "\"V\"");

  DmnRplValues(&UsrHeader, "axes", "\"xyz\"");

  DmnRplValues(&UsrHeader, "sequ", "k+,j-,i+");

  DmnRplValues(&UsrHeader, "refc", NULL);

  DmnRplValues(&UsrHeader, "refd", NULL);

  DmnRplValues(&UsrHeader, "einheit", NULL);

  DmnRplValues(&UsrHeader, "T1", MsgDateString(T1));

  DmnRplValues(&UsrHeader, "T2", MsgDateString(T2));

  if (*locl) {

    char s[256];

    sprintf(s, "\"%s\"", locl);

    DmnRplValues(&UsrHeader, "locl", s);

  }

  else DmnRplValues(&UsrHeader, "locl", NULL);

  if (!strcmp(locl, "german"))  { src = '.';  dst = ','; }

  else                          { src = ',';  dst = '.'; }

  DmnRplChar(&UsrHeader, "xmin", src, dst); 

  DmnRplChar(&UsrHeader, "ymin", src, dst); 

  DmnRplChar(&UsrHeader, "valid", src, dst); 

  DmnRplChar(&UsrHeader, "sk|zk", src, dst); 

  DmnRplChar(&UsrHeader, "dd|delta", src, dst); 

  sprintf(t, "%1.4f", valid);

  DmnRplValues(&UsrHeader, "valid", t);

  if (Concentration) {

    sprintf(gname, "%s-%s00%s%s", name, CfgYearString, CfgAddString, GrdName); //-2006-10-31

    DmnRplValues(&UsrHeader, "file", TxtQuote(gname));

    sprintf(fname, "%s%s-%s00%s%s", Path, name, CfgYearString, CfgAddString, GrdName); //-2006-10-31

    DmnRplValues(&UsrHeader, "form", (odor) ? of : CmpConForm);

    DmnRplValues(&UsrHeader, "unit", CmpConUnit);

    sprintf(srefv, "%10.3e", CmpConRef);

    DmnRplValues(&UsrHeader, "refv", (odor) ? "0" : srefv);  //-2004-12-09   

    DmnRplValues(&UsrHeader, "artp", (odor) ? ot : "\"C\"");

    DmnWrite(fname, UsrHeader.s, &SumCon);                          eG(31);

    vMsg("TMT: Datei %s ausgeschrieben.", fname);

    if (Scatter) {

      sprintf(fname, "\"%s-%s00%s%s\"", name, CfgYearString, CfgDevString, GrdName);  //-2006-10-31

      DmnRplValues(&UsrHeader, "file", fname);

      sprintf(fname, "%s%s-%s00%s%s", Path, name, CfgYearString, CfgDevString, GrdName); //-2006-10-31

      if (odor) {                                                   //-2004-06-10

        DmnRplValues(&UsrHeader, "form", "\"sct%5.3f\"");

        DmnRplValues(&UsrHeader, "unit", CmpConUnit);

        DmnRplValues(&UsrHeader, "refv", "0.0");       //-2004-12-09

        DmnRplValues(&UsrHeader, "artp", oe);          //-2004-12-09

      }

      else {

        DmnRplValues(&UsrHeader, "form", "\"dev%(*100)5.1f\"");

        DmnRplValues(&UsrHeader, "unit", NULL);

        DmnRplValues(&UsrHeader, "refv", "0.03");

        DmnRplValues(&UsrHeader, "artp", "\"E\"");

      }

      DmnWrite(fname, UsrHeader.s, &SumDev);                        eG(32);

      vMsg(_file_$_written_, fname);

    }

  }

  if (Maxima) {

    int ii, nn;

    ARYDSC *pa;

    DmnRplValues(&UsrHeader, "valid", NULL);                          //-2002-08-10

    for (ii=0; ii<2; ii++) {

      nn = (ii) ? 0: CmpDayNum;

      sprintf(gname, "%s-%s%02d%s%s", name, CfgDayString, nn, CfgAddString, GrdName); //-2006-10-31

      DmnRplValues(&UsrHeader, "file", TxtQuote(gname));

      sprintf(fname, "%s%s-%s%02d%s%s", Path, name, CfgDayString, nn, CfgAddString, GrdName); //-2006-10-31

      DmnRplValues(&UsrHeader, "form", CmpDayForm);

      DmnRplValues(&UsrHeader, "unit", CmpConUnit);

      sprintf(srefv, "%10.3e", CmpDayRef);

      DmnRplValues(&UsrHeader, "refv", srefv);

      sprintf(t, "\"CD%02d\"", (ii) ? 0 : nh);                        //-2002-08-10

      DmnRplValues(&UsrHeader, "artp", t);                            //-2002-07-24

      sprintf(t, "%d", (ii) ? 0 : nh);                                //-2002-08-10

      DmnRplValues(&UsrHeader, "exceed", t);

      pa = (ii) ? &SupCon : &MaxCon;                                  //-2001-07-11

      DmnWrite(fname, UsrHeader.s, pa);                               eG(33);

      vMsg(_file_$_written_, fname);

      if (Scatter) {

        sprintf(fname, "\"%s-%s%02d%s%s\"", name, CfgDayString, nn, CfgDevString, GrdName);  //-2006-10-31

        DmnRplValues(&UsrHeader, "file", fname);

        sprintf(fname, "%s%s-%s%02d%s%s", Path, name, CfgDayString, nn, CfgDevString, GrdName); //-2006-10-31

        DmnRplValues(&UsrHeader, "form", "\"dev%(*100)5.1f\"");

        DmnRplValues(&UsrHeader, "unit", NULL);

        DmnRplValues(&UsrHeader, "refv", "0.30");

        DmnRplValues(&UsrHeader, "artp", "\"E\"");

        pa = (ii) ? &SupDev : &MaxDev;                                //-2001-07-11

        DmnWrite(fname, UsrHeader.s, pa);                             eG(34);

        vMsg(_file_$_written_, fname);

      }

      sprintf(fname, "\"%s-%s%02d%s%s\"", name, CfgDayString, nn, CfgIndString, GrdName); //-2006-10-31

      DmnRplValues(&UsrHeader, "file", fname);

      sprintf(fname, "%s%s-%s%02d%s%s", Path, name, CfgDayString, nn, CfgIndString, GrdName); //-2006-10-31

      DmnRplValues(&UsrHeader, "form", "\"day%4hd\"");

      DmnRplValues(&UsrHeader, "unit", NULL);

      DmnRplValues(&UsrHeader, "refv", NULL);

      DmnRplValues(&UsrHeader, "artp", "\"N\"");

      pa = (ii) ? &SupDay : &MaxDay;                                  //-2001-07-11

      DmnWrite(fname, UsrHeader.s, pa);                               eG(35);

      vMsg(_file_$_written_, fname);

      if (nn == 0)  break;                                            //-2001-07-11

      nn = 0;                                                         //-2001-07-11

    }

  }

  if (Deposition) {

    sprintf(t, "%1.4f", valid);                                       //-2002-08-10

    DmnRplValues(&UsrHeader, "valid", t);                             //-2002-08-10

    DmnRplValues(&UsrHeader, "sequ", "j-,i+");

    DmnRplValues(&UsrHeader, "axes", "\"xy\"");

    sprintf(gname, "%s-dep%s%s", name, CfgAddString, GrdName);  //-2006-10-31

    DmnRplValues(&UsrHeader, "file", TxtQuote(gname));

    sprintf(fname, "%s%s-dep%s%s", Path, name, CfgAddString, GrdName); //-2006-10-31

    DmnRplValues(&UsrHeader, "form", CmpDepForm);

    DmnRplValues(&UsrHeader, "unit", CmpDepUnit);

    sprintf(srefv, "%10.3e", CmpDepRef);

    DmnRplValues(&UsrHeader, "refv", srefv);

    DmnRplValues(&UsrHeader, "artp", "\"XD\"");                       //-2002-07-24

    DmnWrite(fname, UsrHeader.s, &SumDep);                            eG(36);

    vMsg(_file_$_written_, fname);

    if (Scatter) {

      sprintf(fname, "\"%s-dep%s%s\"", name, CfgDevString, GrdName); //-2006-10-31

      DmnRplValues(&UsrHeader, "file", fname);

      sprintf(fname, "%s%s-dep%s%s", Path, name, CfgDevString, GrdName);  //-2006-10-31

      DmnRplValues(&UsrHeader, "form", "\"dev%(*100)5.1f\"");

      DmnRplValues(&UsrHeader, "unit", NULL);

      DmnRplValues(&UsrHeader, "refv", "0.03");

      DmnRplValues(&UsrHeader, "artp", "\"E\"");

      DmnWrite(fname, UsrHeader.s, &SumDps);                        eG(37);

      vMsg(_file_$_written_, fname);

    }

  }

  TxtClr(&UsrHeader);                                             //-2006-02-15

  AryFree(&SumCon);

  AryFree(&SumDev);

  AryFree(&SumDep);

  AryFree(&SumDps);

  AryFree(&MaxCon);

  AryFree(&MaxDev);

  AryFree(&MaxDay);

  AryFree(&SupCon);

  AryFree(&SupDev);

  AryFree(&SupDay);

  AryFree(&VecDos);

  AryFree(&VecDay);

  AryFree(&SumDos);

  AryFree(&CurDos);

  return 0;

eX_4: eX_5: eX_6:

eX_20: eX_21: eX_22: eX_23:

eX_31: eX_32: eX_33: eX_34: eX_35: eX_36: eX_37:

  eMSG(_cant_write_array_);

eX_10: eX_11: eX_12: eX_13: eX_14: eX_15: eX_16: eX_17: eX_18: eX_19:

  eMSG(_internal_error_);

  }



/*======================================================================= main

*/

static void define_component( char *s ) {

  char name[40], format[40], units[40], t[400], *pc, uc[40], un[40];

  int n, dy, dd, dh, dn, nh, nd;

  float fc, fn, ry, rd, rh, rn, vd;

  char locale[256]="C";

  strcpy(locale, setlocale(LC_NUMERIC, NULL));      //-2003-07-07

  setlocale(LC_NUMERIC, "C");

  strcpy(t, s);

  for (pc=t; (*pc); pc++)  if (*pc == ';')  *pc = ' ';

  n = sscanf(t, "%s %f %f %s %f %s %f %d %d %f %d %d %f %d %f %d",

    name, // name of the component

    &vd,  // deposition velocity

    // concentration:

    &fc,  // conversion factor

    uc,   // units used

    // deposition:

    &fn,  // conversion factor

    un,   // units used

    // yearly average:

    &ry,  // reference value

    &dy,  // number of digits after the decimal point

    // daily average:

    &nd,  // number of allowed exceedings (days)

    &rd,  // reference value for daily maximum

    &dd,  // number of digits after the decimal point

    // hourly average:

    &nh,  // number of allowed exceedings (hours)

    &rh,  // reference value for hourly maximum

    &dh,  // number of digits after the decimal point

    // deposition:

    &rn,  // reference value

    &dn   // number of digits after the decimal point

    );

  if (n < 16)  goto done;

  MsgLow(name);

  strcpy(CmpName, name);

  strcpy(CmpFamily, name);

  CmpConFac = fc;

  CmpConRef = ry;

  if (dy < 0)  strcpy(format, "\"con%10.3e\"");         //-2003-02-21

  else         sprintf(format, "\"con%%5.%df\"", dy);

  strcpy(CmpConForm, format);

  sprintf(units, "\"%s\"", uc);

  strcpy(CmpConUnit, units);

  CmpDayNum = nd;

  CmpDayRef = rd;

  if (rd <= 0)  nd = -1;

  if (dd < 0)  strcpy(format, "\"con%10.3e\"");         //-2003-02-21

  else         sprintf(format, "\"con%%5.%df\"", dd);

  strcpy(CmpDayForm, format);

  CmpHourNum = nh;

  CmpHourRef = rh;

  if (rh <= 0)  nh = -1;

  if (dh < 0)  strcpy(format, "\"con%10.3e\"");         //-2003-02-21

  else         sprintf(format, "\"con%%5.%df\"", dh);

  strcpy(CmpHourForm, format);

  CmpDepFac = fn;

  CmpDepRef = rn;

  if (dn < 0)  strcpy(format, "\"dep%10.3e\"");         //-2003-02-21

  else         sprintf(format, "\"dep%%5.%df\"", dn);

  strcpy(CmpDepForm, format);

  sprintf(units, "\"%s\"", un);

  strcpy(CmpDepUnit, units);

  if (CmpDepFac > 0)  Deposition = 1;

  if (CHECK) vMsg(

    "Name=%s, fc=%10.3e, ry=%10.3e, mc=%s, uc=%s, nd=%d, rd=%10.3e, md=%s, "

    "fn=%10.3e, rn=%10.3e, mn=%s, un=%s",

    CmpName, CmpConFac, CmpConRef, CmpConForm, CmpConUnit, CmpDayNum,

    CmpDayRef, CmpDayForm, CmpDepFac, CmpDepRef, CmpDepForm, CmpDepUnit);

done:

  setlocale(LC_NUMERIC, locale);

}



int TmtMain( char *s )

  {

  dP(TmtMain);

  int i, l;

  int ninvalid = 0;

  long rc;

  char *pc;

  if (CHECK) vMsg("TmtMain(%s)", s);

  if (*s) {

    pc = s;

    if (*pc == '-')

      switch (pc[1]) {

        case '0': Nh=0; Ns=0; Ng=0; Anzahl=0; Folge=1; Deposition=0;

                  strcpy(DosName, "d%04lda00.arr");

                  *GrdName = 0;

                  *CmpName = 0;

                  WriteHeader = 1;

                  Kref = -999;

                  Kmax = -999;

                  break;

        case 'A': strcpy(locl, pc+2);                                 //-2003-07-07

                  break;

        case 'a': AverageOnly = 1;                                    //-2007-02-03

                  break;

        case 'b': sscanf(pc+2, "%ld", &Folge);

                  break;

        case 'f': strcpy(DosName, pc+2);

                  break;

        case 'g': strcpy(GrdName, pc+2);

                  break;

        case 'k': Kref = 1;

                  sscanf(pc+2, "%d", &Kref);

                  break;

        case 'm': Kmax = 1;

                  sscanf(pc+2, "%d", &Kmax);

                  break;

        case 'n': sscanf(pc+2, "%ld", &Anzahl);

                  break;

        case 'q': MsgQuiet = 1;

                  sscanf(pc+2, "%ld", &MsgQuiet);

                  break;

        case 's': define_component(pc+2);

                  break;

        default:  vMsg(_unknown_option_$_, pc);

                  break;

      }

    else {

      strcpy(Path, pc);

      l = strlen(Path);

      for (i=0; i
      if ((l) && (Path[l-1]!='/') && (Path[l-1]!=':')) {

        Path[l] = '/';

        Path[l+1] = 0;

      }

    }

    return 0;

  }

  if (WriteHeader) vMsg(_evaluation_$_, CmpName);

  if (*locl) {                                                        //-2003-07-07

    strcpy(locale, setlocale(LC_NUMERIC, NULL));

    setlocale(LC_NUMERIC, locl);

  }

  Nh = 0;

  Ng = 0;

  Nsum = 0;

  Concentration = (CmpConRef > 0);

  Deposition    = (CmpDepRef > 0); 

  Maxima        = (CmpDayRef > 0 && !AverageOnly);          //-2007-02-03

  SumValid = 0;

  for (Ia=0; ; ) {

    rc = ReadArray();                                                   eG(1);

    if (rc==0 && Anzahl>0 && Ia
    if (rc == 0)  break;

    Ia++;

    rc = CompareArray();                                                eG(3);

    if (rc > 0)  ninvalid++;

    AryFree(&CurDos);

    if (Anzahl>0 && Ia>=Anzahl)  break;

    }

  if (Anzahl>1 && WriteHeader>0)

    vMsg(_$_daily_averages_$_, Ia, ninvalid);

  WriteHeader = 0;

  WriteResult();                                                        eG(4);

  if (*locl)  setlocale(LC_NUMERIC, locale);                            //-2003-07-07

  return 0;

eX_1:

  eMSG(_cant_read_$_, Eingabe);

eX_2:

  eMSG(_only_$_$_found_, Ia, Anzahl);

eX_3:

  eMSG(_error_maximum_);

eX_4:

  eMSG(_error_writing_);

  }



//=============================================================================

//

// history:

//

// 2002-06-21 lj 0.13.0  final test version

// 2002-07-24 lj 0.13.1  "artp" corrected

// 2002-08-10 lj 0.13.2  use of "valid" corrected, setting of "exceed", RefDays 

// 2002-09-24 lj 1.0.0   final release candidate

// 2002-09-28 lj 1.0.1   write highest values even if nh==0

// 2003-02-21 lj 1.1.1   optional scientific notation

// 2003-07-07 lj 1.1.8   localisation

// 2004-06-10 lj 2.0.0   odor

// 2004-12-09 uj 2.1.1   artp for odor: FHP, EHP

// 2005-03-17 uj 2.2.0   version number upgrade

// 2005-08-06 uj 2.2.1   rounding of double returned by MsgDateSeconds() 

// 2006-02-15 lj 2.2.9   TxtClr inserted in WriteResult

// 2006-07-13 uj 2.2.12  case gruppen<=4 corrected

// 2006-10-18 lj 2.3.0   externalization of strings

// 2006-10-31 lj 2.3.1   configuration data

// 2007-02-03 uj 2.3.5   option -a (AverageOnly)

//

//=============================================================================