//======================================================== IBJprf.c
//
// Boundary layer profile functions for IBJ-programs
// =================================================
//
// Copyright (C) Janicke Consulting, Dunum, Germany, 1998-2001
// email: info@janicke.de
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
//
// history:
//
// last change:
// 2001-07-28 uj  
//
//==================================================================

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <ctype.h>
#include <math.h>

#define  PRF_FU 2.5
#define  PRF_FV 2.0
#define  PRF_FW 1.3
#define  PRF_FFU 0.17
#define  PRF_FFV 0.17
#define  PRF_FFW 0.13
#define  PRF_PS  4.7
#define  PRF_PN  1.0
#define  PRF_PL  15.0
#define  PRF_KAPPA 0.4
#define  PRF_ALPHA 0.3
#define  PRF_FCOR  1.e-4
#define  PRF_TMPGRADIENT -0.0098
#define  PRF_GRAV 9.81
#define VDI_EX 1.0
#define VDI_C0 5.7
#define V_STANDARD 21
#define V_VDI      26

extern char *PrfVersion;

float PrfGetTmp(float z, float lm, float ust, float z0, float tmp, float theta0);
float PrfGetU(float z, float lm, float hm, float ust, float z0, float interval);
float PrfGetSigU(float z, float lm, float hm, float ust, float z0, float interval);
float PrfGetSigV(float z, float lm, float hm, float ust, float z0, float interval);
float PrfGetSigW(float z, float lm, float hm, float ust, float z0, float interval);
float PrfGetKU(float z, float lm, float hm, float ust, float z0, float interval);
float PrfGetKV(float z, float lm, float hm, float ust, float z0, float interval);
float PrfGetKW(float z, float lm, float hm, float ust, float z0, float interval);
float PrfGetHM(float, float);
void  PrfSetKappa(float kappa);
void  PrfSetD0(float d0);
void  PrfSetBeta(float beta);
void  PrfSetVersion(int version);

extern char *PrfVersion = "1.3";

float Ha, Ua, Z0, Ust, Lm, Hm, Zz[101], Interval, D0=0., Kappa=PRF_KAPPA;
float Beta = 1;
int   Nz = 100;
char  Name[512];
int   Version = 21;
// Version 2.1: Standard Version (cf. BZU 2, 2000)
// Version 2.6: Profiles of VDI 3783/8 (draft)
int   ImplVersions[] = { V_STANDARD, V_VDI };
int   NVersions = 2;
int   StdVersion = V_STANDARD;
int   DEBUG = 0;

//----------------------------------------------------------------- PrfSetKappa
void PrfSetKappa(float kappa) { Kappa = kappa; }

//-------------------------------------------------------------------- PrfSetD0
void PrfSetD0(float d0) { D0 = d0; }

//------------------------------------------------------------------ PrfSetBeta
void PrfSetBeta(float beta) { Beta = beta; }

//--------------------------------------------------------------- PrfSetVersion
void PrfSetVersion(int version) {
  int i;
  Version = version;
  for (i=0; i<NVersions; i++) if (ImplVersions[i] == Version) break;
  if (i >= NVersions) Version = StdVersion;
  if (DEBUG) printf("Version set to %d\n", Version);
}

//------------------------------------------------------------------- PrfGetTmp
float PrfGetTmp(float z, float lm, float ust, float z0, float mean, float tmp0) {
  float zeta, zeta0, psi, y, y0, theta, tst;
  zeta = (z+z0)/lm;
  zeta0 = z0/lm;
  tst = ust*ust*mean/(Kappa*PRF_GRAV*lm);
  if (lm > 0) {
    psi = -PRF_PS*(zeta - zeta0);
  }
  else {
    y = sqrt(1. - PRF_PL*zeta);
    y0 = sqrt(1. - PRF_PL*zeta0);
    psi = 2*log((1.+y)/(1.+y0));
  }
  theta = tmp0 + PRF_PN*tst/Kappa*(log(zeta/zeta0) - psi);
  return (theta + PRF_TMPGRADIENT*z);
}

//--------------------------------------------------------------------- PrfGetU
float PrfGetU(float z, float lm, float hm, float ust, float z0, float interval) {
  float zeta, zeta0, psi, x1, x2, x3, x, x0, u, fac=1.;
  if (D0 > 0.) {
    if (z > 2*D0) {
      fac = 1.;
      z = z - D0;
    }
    else {
      fac = 0.5*z/D0;
      z = D0;
    }
  }
  zeta = (z+z0)/lm;
  zeta0 = z0/lm;
  if (lm > 0.) {
    x1 = -5*zeta0;
    x2 = 7*log(2.) - log(zeta0) - 5*zeta0 - 4;
    x3 = 8*log(20) - 11.165 - log(2.) - log(zeta0) - 5*zeta0;
    if (zeta < 0.5)
      psi = log(zeta/zeta0) + 5*zeta + x1;
    else if (zeta <= 10.)
      psi = 8*log(zeta) + 4.25/zeta - 0.5/zeta/zeta + x2;
    else
      psi = 0.7585*zeta + x3;
    u = ust/Kappa*psi;
  }
  else {
    x = pow(1 - 15*zeta, 0.25);
    x0 = pow(1 - 15*zeta0, 0.25);
    x1 = (1+x)*(1+x)/(1+x0)/(1+x0) * (1+x*x)/(1+x0*x0);
    psi = log(x1) - 2*(atan(x) - atan(x0));
    u = ust/Kappa*(log(zeta/zeta0) - psi);
  }
  if (u < 0) u = 0.;
  return fac*u;
}

//------------------------------------------------------------------ PrfGetSigU
float PrfGetSigU(float z, float lm, float hm, float ust, float z0, float interval) {
  float av, sig, wst;
  if (D0 > 0.)
    if (z <= 2*D0) z = 0.5*z;
    else           z = z - D0;
  if (Version == V_VDI) {
    if (z < z0) z = z0;
    if (lm > 0) {
      sig = 2.4*ust*exp(-VDI_EX*z/hm);
    }
    else {
      wst = ust*pow(-hm/Kappa/lm, 1./3.);
      sig = pow(2.4*ust, 3.) + pow(0.59*wst, 3.);
      sig = pow(sig, 0.333333) * exp(-VDI_EX*z/hm);
    }
  }
  else {
    if (interval < 600.) interval = 600.;
    if (interval > 3600.) interval = 3600.;
    av = pow(interval/3600., 0.2);
    sig = PRF_FU*av*ust;
    if (lm < 0.)
      sig *= pow(1. - 0.064*hm/lm, 1./3.);
  }
  if (sig < 0.01*ust) sig = 0.01*ust;
  return sig;
}

//------------------------------------------------------------------ PrfGetSigV
float PrfGetSigV(float z, float lm, float hm, float ust, float z0, float interval) {
  float av, sig, wst;
  if (D0 > 0.)
    if (z <= 2*D0) z = 0.5*z;
    else           z = z - D0;
  if (Version == V_VDI) {
    if (z < z0) z = z0;
    if (lm > 0.) {
      sig = 1.8*ust*exp(-VDI_EX*z/hm);
    }
    else {
      wst = ust*pow(-hm/Kappa/lm, 1./3.);
      sig = pow(1.8*ust, 3.) + pow(0.59*wst, 3.);
      sig = pow(sig, 0.333333) * exp(-VDI_EX*z/hm);
    }
  }
  else {
    if (interval < 600.) interval = 600.;
    if (interval > 3600.) interval = 3600.;
    av = pow(interval/3600., 0.2);
    sig = PRF_FV*av*ust;
    if (lm < 0.)
      sig *= pow(1. - 0.064*hm/lm, 1./3.);
  }
  if (sig < 0.01*ust) sig = 0.01*ust;
  return sig;
}

//------------------------------------------------------------------ PrfGetSigW
float PrfGetSigW(float z, float lm, float hm, float ust, float z0, float interval) {
  float sig, mu, wst, a, f, sl, ss;
  if (D0 > 0.)
    if (z <= 2*D0) z = 0.5*z;
    else           z = z - D0;
  a = z/hm;
  if (Version == V_VDI) {
    if (z < z0) z = z0;
    if (lm > 0.) {
      sig = 1.3*ust*exp(-VDI_EX*a);
    }
    else {
      wst = wst = ust*pow(-hm/Kappa/lm, 0.333333);
      sig = pow(1.3*ust*exp(-VDI_EX*a), 3.);
      if (1-0.8*a > 0.)
        sig += pow(1.3*pow(a,0.333333)*(1-0.8*a)*wst, 3.);
      sig = pow(sig, 0.333333);
    }
  }
  else {
    sig = PRF_FW*ust;
    mu = exp(-2*PRF_ALPHA*a);
    if (lm < 0.)
      sig *= pow(mu*mu*mu + 2.10*(-z/lm)*exp(-3*a), 0.333333);
    else
      sig *= mu;
  }
  if (sig < 0.01*ust) sig = 0.01*ust;
  return sig;
}

//-------------------------------------------------------------------- PrfGetE
float PrfGetE(float z, float lm, float hm, float ust, float z0, float interval) {
  float e, wst, a, f1;
  if (z < z0) z = z0;
  if (lm > 0.) {
    e = 1. + 4*z/lm;
  }
  else {
    a = z/hm;
    wst = -hm/Kappa/lm;
    e = ((1-a)*(1-a) + 2.5*Kappa*a) + Kappa*z*wst/hm*(1.5 - 1.3*pow(a, 0.333333));
    if (e < 1.) e = 1.;
  }
  return (e*ust*ust*ust/Kappa/z);
}


//-------------------------------------------------------------------- PrfGetKW
float PrfGetKW(float z, float lm, float hm, float ust, float z0, float interval) {
  float k, sig, kmin, e, a, f, wst, sf, zz;
  zz = z;
  sig = PrfGetSigW(zz, lm, hm, ust, z0, interval);
  if (D0 > 0.)
    if (z <= 2*D0) z = 0.5*z;
    else           z = z - D0;
  a = z/hm;
  if (Version == V_VDI) { 
    e = PrfGetE(zz, lm, hm, ust, z0, interval);
    k = 2*sig*sig*sig*sig/VDI_C0/e;
  }
  else {
    k = Kappa*ust*z/PRF_PN;
    f = z/lm;
    if (lm > 0.)
      k *= 1./(1 + PRF_PS*f) * exp(-6*PRF_ALPHA*a);
    else {
      sf = sqrt(0.01*hm/(-lm));
      a = ((1 - a) > 0.) ? (1-a)*(1-a) : 0.;
      k *= pow(exp(-18*PRF_ALPHA*z/hm) + (3*PRF_PS*(-f) + 120*f*f)*a*a*a, 0.333333);
    }
  }
  kmin = Beta*sig*sig*z0/ust;
  if (k < kmin) k = kmin;
  return k;
}

//-------------------------------------------------------------------- PrfGetKU
float PrfGetKU(float z, float lm, float hm, float ust, float z0, float interval) {
  float av, u, sig, k, kmin, e, zz;
  zz = z;
  sig = PrfGetSigU(zz, lm, hm, ust, z0, interval);
  if (D0 > 0.)
    if (z <= 2*D0) z = 0.5*z;
    else           z = z - D0;
  if (Version == V_VDI) {
    e = PrfGetE(zz, lm, hm, ust, z0, interval);
    k = 2*sig*sig*sig*sig/VDI_C0/e;
  }
  else {
    if (interval < 600.) interval = 600.;
    if (interval > 3600.) interval = 3600.;
    av = pow(interval/3600., 0.2);
    u = PrfGetU(zz, lm, hm, ust, z0, interval);
    k = PRF_FFU*av*u/(10*PRF_FU*ust)*hm*sig;
  }
  kmin = Beta*sig*sig*z0/ust;
  if (k < kmin) k = kmin;
  return k;
}

//-------------------------------------------------------------------- PrfGetKV
float PrfGetKV(float z, float lm, float hm, float ust, float z0, float interval) {
  float av, u, sig, k, kmin, e, zz;
  zz = z;
  sig = PrfGetSigV(zz, lm, hm, ust, z0, interval);
  if (D0 > 0.)
    if (z <= 2*D0) z = 0.5*z;
    else           z = z - D0;
  if (Version == V_VDI) {
    e = PrfGetE(zz, lm, hm, ust, z0, interval);
    k = 2*sig*sig*sig*sig/VDI_C0/e;
  }
  else {
    if (interval < 600.) interval = 600.;
    if (interval > 3600.) interval = 3600.;
    av = pow(interval/3600., 0.2);
    u = PrfGetU(zz, lm, hm, ust, z0, interval);
    k = PRF_FFV*av*u/(10*PRF_FV*ust)*hm*sig;
  }
  kmin = Beta*sig*sig*z0/ust;
  if (k < kmin) k = kmin;
  return k;
}

//-------------------------------------------------------------------- PrfGetHM
float PrfGetHM(float lm, float ust) {
  float hm;
  if (lm > 0.) {
    hm = PRF_ALPHA*ust/PRF_FCOR;
    if (lm < ust/PRF_FCOR) hm *= sqrt(PRF_FCOR*lm/ust);
  }
  else {
    if (lm < -200) hm = 800.;
    else           hm = 1100.;
  }
  return hm;
}

#ifdef MAIN /*###############################################################*/
#define  ISOPTION(a)   ((a) == '-')

//-------------------------------------------------------------------- profiles
void profiles() {
  int i;
  float vx[40], su[40], sv[40], sw[40], ku[40], kv[40], kw[40], z;
  FILE *fp;

  if (Nz <= 0) {
    Nz = 30;
    Zz[0] = 0.;
    i = 1;
    while (i++ < 20) Zz[i] = Zz[i-1] + 10.;
    while (i++ < 30) Zz[i] = Zz[i-1] + 30.;
    while (i++ < 40) Zz[i] = Zz[i-1] + 50.;
  }
  if (Ust == 0.) {
    if (Ha <= 0.) Ha = 10.;
    if (Ua <= 0.) Ua = 4.;
  }
  if (Z0 <= 0.) Z0 = 1.0;
  if (Lm == 0.) Lm = 9999.;
  if (Ust <= 0.) Ust = Ua/PrfGetU(Ha, Lm, Hm, 1., Z0, Interval);

  if (Hm <= 0.) Hm = PrfGetHM(Lm, Ust);

  for (i=0; i<=Nz; i++) {
    z = Zz[i];
    vx[i] = PrfGetU(z, Lm, Hm, Ust, Z0, Interval);
    su[i] = PrfGetSigU(z, Lm, Hm, Ust, Z0, Interval);
    sv[i] = PrfGetSigV(z, Lm, Hm, Ust, Z0, Interval);
    sw[i] = PrfGetSigW(z, Lm, Hm, Ust, Z0, Interval);
    ku[i] = PrfGetKU(z, Lm, Hm, Ust, Z0, Interval);
    kv[i] = PrfGetKV(z, Lm, Hm, Ust, Z0, Interval);
    kw[i] = PrfGetKW(z, Lm, Hm, Ust, Z0, Interval);
  }

  remove(Name);
  fp = (strlen(Name) > 1) ? fopen(Name, "w") : stdout;
  fprintf(fp, "   ' Version=%d\n", Version);
  fprintf(fp, "   ' us=%8.4f lm=%8.2f hm=%8.2f\n", Ust, Lm, Hm);
  fprintf(fp, "   ' nz=%8.0d ha=%8.2f ua=%8.3f\n", Nz, Ha, Ua);
  fprintf(fp, "   ' z0=%8.4f d0=%8.2f in=%8.2f kp=%8.2f Beta=%5.3f\n",
    Z0, D0, Interval, Kappa, Beta);
  fprintf(fp, "zz  ");
  for (i=0; i<Nz; i++) fprintf(fp, " %11.3e", Zz[i]); fprintf(fp, "\n");
  fprintf(fp, "vx  ");
  for (i=0; i<Nz; i++) fprintf(fp, " %11.3e", vx[i]); fprintf(fp, "\n");
  fprintf(fp, "su  ");
  for (i=0; i<Nz; i++) fprintf(fp, " %11.3e", su[i]); fprintf(fp, "\n");
  fprintf(fp, "sv  ");
  for (i=0; i<Nz; i++) fprintf(fp, " %11.3e", sv[i]); fprintf(fp, "\n");
  fprintf(fp, "sw  ");
  for (i=0; i<Nz; i++) fprintf(fp, " %11.3e", sw[i]); fprintf(fp, "\n");
  fprintf(fp, "ku  ");
  for (i=0; i<Nz; i++) fprintf(fp, " %11.3e", ku[i]); fprintf(fp, "\n");
  fprintf(fp, "kv  ");
  for (i=0; i<Nz; i++) fprintf(fp, " %11.3e", kv[i]); fprintf(fp, "\n");
  fprintf(fp, "kw  ");
  for (i=0; i<Nz; i++) fprintf(fp, " %11.3e", kw[i]); fprintf(fp, "\n");
}

//------------------------------------------------------------------------ Help
void Help() {
  printf("Calculation of meteorological standard profiles (LASAT Version 2.1)\n");
  printf("usage: profiles [options]\n");
  printf("options:\n");
  printf("  -d<number>   Displacement height (m, default: 0)\n");
  printf("  -i<number>   Interval (s, default: 3600)\n");
  printf("  -k<number>   Kappa (default: 0.4)\n");
  printf("  -l<number>   Monin-Obukhov length (m, default: 9999)\n");
  printf("  -m<number>   Height of the mixing layer (m)\n");
  printf("  -h<number>   Anemometer height (m, default: 10)\n");
  printf("  -u<number>   Wind velocity at anemometer height (m/s, default: 4)\n");
  printf("  -s<number>   Friction velocity (m/s)\n");
  printf("  -z<number>   Roughness length (m, default: 1)\n");
  printf("  -v<z1>,<z2>,...  Vertical grid (m)\n");
  printf("  -V<version>  Version (default 2.1: BZU 2 (2000), 2.6: VDI 3783/3 draft)\n");
  printf("  -f<file>     File name\n");
  printf("\n");
  exit(0);
}

//------------------------------------------------------------------------ main
void main( int argc, char *argv[] ) {
  int i, j;
  float z;
  char *pc;
  printf("IBJprf, version %s of %s, %s\n", PrfVersion, __DATE__, __TIME__);
  printf("1-dim. boundary layer model\n");
  printf("Copyright (C) Janicke Consulting, Dunum, Germany, 2000-2001\n");
  printf("IBJprf comes with ABSOLUTELY NO WARRANTY\n");
  printf("IBJprf is free software under the GNU PUBLIC LICENCE\n\n");
  Lm = 0.;
  Hm = 0.;
  Ua = 0.;
  Ha = 0.;
  Ust = 0.;
  Nz = 0.;
  Z0 = 0.;
  D0 = 0.;
  Kappa = PRF_KAPPA;
  Interval = 3600.;
  sprintf(Name, "");
  if (argc < 2) Help();
  for (i=1; i<argc; i++) {
    pc = argv[i];
    if (ISOPTION(*pc)) {
      switch (pc[1]) {
        case 'l': j=sscanf(pc+2, "%f", &Lm);
                  if (j != 1) Help();
                  break;
        case 'm': j=sscanf(pc+2, "%f", &Hm);
                  if (j != 1) Help();
                  break;
        case 'u': j=sscanf(pc+2, "%f", &Ua);
                  if (j != 1) Help();
                  break;
        case 'h': j=sscanf(pc+2, "%f", &Ha);
                  if (j != 1) Help();
                  break;
        case 'i': j=sscanf(pc+2, "%f", &Interval);
                  if (j != 1) Help();
                  break;
        case 'd': j=sscanf(pc+2, "%f", &D0);
                  if (j != 1) Help();
                  break;
        case 'k': j=sscanf(pc+2, "%f", &Kappa);
                  if (j != 1) Help();
                  break;
        case 's': j=sscanf(pc+2, "%f", &Ust);
                  if (j != 1) Help();
                  break;
        case 'z': j=sscanf(pc+2, "%f", &Z0);
                  if (j != 1) Help();
                  break;
        case 'f': j=sscanf(pc+2, "%s", &Name);
                  if (j != 1) Help();
                  break;
        case 'V': j=sscanf(pc+2, "%f", &z);
                  if (j != 1) Help();
                  Version = (int)(z*10);
                  PrfSetVersion(Version);
                  break;
        case 'v': j = 0;
                  pc += 2;
                  while (1 == sscanf(pc, "%f", &z)) {
                    while (*pc != ',') pc++;
                    pc++;
                    Zz[j++] = z;
                  }
                  Nz = j;
                  break;
      }
    }
    else Help();
  }
  profiles();
  printf("Program IBJprf finished.\n");
  exit(0);
}

#endif  /*###################################################################*/
