brook90.h

#ifndef  M1D_MOD_WATERCYCLE_BROOK90_H
#define  M1D_MOD_WATERCYCLE_BROOK90_H
 
#include  <string>
#include "watercyclemodule.h"
 
using std::string;
 
/*  WaterCycleBROOK90
 *  The class WaterCycleBROOK90 contains all Methods and data from the WaterCycleBROOK90 Module.
 *
 *
*/
class  WaterCycleBROOK90  :  public  WaterCycleModule {
    public:
            WaterCycleBROOK90(ModelCore & modelCore, unsigned int timeInterval) : WaterCycleModule(modelCore, timeInterval),
                                modelCore_(modelCore),
                                time_(modelCore.getLandscape().getTime()),
                                si_(modelCore.getInput().getSite()),
                                sipar_(modelCore.getInput().getParameter().getSitePar()),
                                sppar_(modelCore.getInput().getParameter().getSpeciesPar()),
                                ac_(modelCore.getSubstates().getAc()),
                                mc_(modelCore.getSubstates().getMc()),
                                sc_(modelCore.getSubstates().getSc()),
                                wc_(const_cast<WaterCycle &>(modelCore.getSubstates().getWc())),
                                ph_(modelCore.getSubstates().getPh()),
                                vs_(modelCore.getSubstates().getVs())
                    {
                        this->setName("WaterCycleBROOK90");
                    }
 
            ~WaterCycleBROOK90() {
            }
 
 
            void  initialize()
            {
                WaterCycleModule::initialize();
                GAMMA = 0.0662;  // GAMMA CONSTANTE PSYCHROMETRIQUE EN KPA PAR DEGRE K
                PI = 3.14159265; // number pi
                CVAIR    = 1297.0;       // volumetric heat capacity of air, J m-3 K-1), wikipedia; BROOK90: 1240.0
            }
 
            void  solve() {
                this->run_();
            }
 
 
private:
        ModelCore & modelCore_;
        Timer & time_;
        Site  & si_;
        SitePar &  sipar_;
        SpeciesPar &  sppar_;
 
        AirChemistry const & ac_;
        MicroClimate const & mc_;
        SoilChemistry const & sc_;
        WaterCycle & wc_;
        Physiology const & ph_;
        VegStructure const & vs_;;
 
        // declaration of module specific parameters
      static const double DT;          // time step for DFILE interval,  must be 1 d
      double GAMMA;  // GAMMA CONSTANTE PSYCHROMETRIQUE EN KPA PAR DEGRE K
      double PI; // number pi
      double CVAIR;       // volumetric heat capacity of air, J m-3 K-1), wikipedia; BROOK90: 1240.0
 
      void run_();
      // functions
      void SUNDS
           (double LAT,double SLOPE,int DOY,double L1,double L2, double DEC,
            double &I0HDAY,double &SLFDAY);
 
      void EQUIVSLP
           (double LAT,double SLOPE,double ASPECT,
            double &L1,double &L2);
 
      void SOILPAR
           (int NLAYER,double *THICK,double *THSAT,double *STONEF,double *THETAF,
             double *PSIF,double *BEXP,double *WETINF,double *PSIM,double *KF,double PSICR,
            double *PSIG,double *SWATMX,double *WETF,double *WETC,double *CHM,
             double *CHN,double *WETNES,double *SWATI,double *KSAT);
 
      void SOILVAR
           (int NLAYER,double *PSIG,double *PSIM,double *WETNES,double *THSAT,
             double *KF,double *BEXP,double *WETF,double *SWATI,
            double *PSITI,double *THETA,double *KK,double &SWAT);
 
      void INFPAR
           (double INFEXP,double IDEPTH,int NLAYER,double *THICK,
            int &ILAYER, double *INFRAC);
 
      void SRFPAR
           (double QDEPTH,int NLAYER,double *THETAF,double *THICK,double *STONEF,double *SWATMX,
            int &QLAYER,double &SWATQX,double &SWATQF);
 
      void RTDEN
           (double *ROOTDEN,int NLAYER,double RDEPTH,double *THICK,
            double *RELDEN);
 
      void AVAILEN
           (double SLRAD,double ALBEDO,double C1,double C2,double C3,double TA,
             double EA,double RATIO,double SHEAT,double CR,double LAI,double SAI,
            double &AA,double &ASUBS);
 
      void ESAT
           (double TA,
            double &ES, double &DELTA);
 
      void SWGRA
           (double UA,double ZA,double HEIGHT,double Z0,double DISP,double Z0C,
             double DISPC,double Z0GS,double LWIDTH,double RHOTP,double NN,
             double LAI,double SAI,
            double &RAA,double &RAC,double &RAS);
 
      void SRSC
           (double RAD,double TA,double VPD,double LAI,double SAI,double GLMIN,
             double GLMAX,double R5,double CVPD,double RM,double CR,double TL,
             double T1,double T2,double TH,
            double &RSC);
 
      void SWPE
           (double AA,double ASUBS,double VPD,double RAA,double RAC,
             double RAS,double RSC,double RSS,double DELTA,
            double &PRATE,double &ERATE);
 
      void TBYLAYER
           (int J,double PTR,double DISPC,double *ALPHA,double *KK,double *RROOTI,
             double RXYLEM,double *PSITI,int NLAYER,double PSICR,int NOOUTF,
            double &ATR,double *ATRANI);
 
      void SWGE
           (double AA,double ASUBS,double VPD,double RAA,double RAS,double RSS,
             double DELTA,double ARATE,
            double &ERATE);
 
      void INTER
           (double RFAL,double PINT,double LAI,double SAI,double FRINTL,
             double FRINTS,double CINTRL,double CINTRS,double DTP,double INTR,
            double &RINT,double &IRVP);
 
      void INTER24
           (double RFAL,double PINT,double LAI,double SAI,double FRINTL,
             double FRINTS,double CINTRL,double CINTRS,double DURATN,double INTR,
            double &RINT,double &IRVP);
 
      void SNOWPACK
           (double RTHR,double STHR,double PSNVP,double SNOEN,double DTP,
             double TA,double MAXLQF,double GRDMLT,
            double &CC,double &SNOW,double &SNOWLQ,double &RSNO,double &SNVP,
             double &SMLT);
 
      void SRFLFR
           (double SWAT,double SWATQX,double QFPAR,double SWATQF,double QFFC,
            double &SAFRAC);
 
      void BYFLFR
           (int BYPAR,int NLAYER,double *WETNES,double *WETF,double QFFC,double QFPAR,
            double *BYFRAC);
 
      void DSLOP
           (double DSLOPE,double LENGTH,double THICK,double STONEF,double PSIM,double KK,
            double &DSFLI);
 
      void VERT
           (double KK,double KK1,double KSAT,double KSAT1,double THICK,double THICK1,
             double PSIT,double PSIT1,double STONE,double STONE1,
            double &VRFLI);
 
      void INFLOW
           (int NLAYER,double DTI,double *INFRAC,double *BYFRAC,double SLFL,
             double *VRFLI,double *DSFLI,double *TRANI,double SLVP,double *SWATMX,
             double *SWATI,double *THETA,double *THICK,double *STONEF,
            double *VV,double *INFLI,double *BYFLI,double *NTFLI);
 
      void ITER
           (int NLAYER,double DTI,double *DPSIDW,double *NTFLI,double *SWATMX,
             double *PSITI,double DSWMAX,double DPSIMX,
            double &DTINEW);
 
      void GWATER
           (double GWAT,double GSC,double GSP,double DT,double VRFLN,
            double &GWFL,double &SEEP);
 
      void CANOPY
           (int DOY,double MAXHT,double RELHT,/*double MAXLAI,*/double RELLAI,double SNOW,
             double SNODEN,double MXRTLN,double MXKPL,double CS,double DENSEF,
            double &HEIGHT,double &LAI,double &SAI,double &RTLEN,double &RPLANT);
 
      void ROUGH
           (double HEIGHT,double ZMINH,double LAI,double SAI,double CZS,double CZR,
             double HS,double HR,double LPC,double CS,
            double &Z0GS,double &Z0C,double &DISPC,double &Z0,double &DISP,double &ZA);
 
      void PLNTRES
           (int NLAYER,double *THICK,double *STONEF,double RTLEN,double *RELDEN,
             double RTRAD,double RPLANT,double FXYLEM,
            double &RXYLEM,double *RROOTI,double *ALPHA);
 
      void WEATHER
           (double TMAX,double TMIN,double DAYLEN,double I0HDAY,double EA,double UW,
             double ZA,double DISP,double Z0,double WNDRAT,double FETCH,double Z0W,
             double ZW,double SOLRAD,double TA,double GTRANS,
            double &SOLRADC,double &TADTM,double &TANTM,double &UA,
             double &UADTM,double &UANTM);
 
      void SNOFRAC
           (double TMAX,double TMIN,double RSTEMP,
            double &SNOFRC);
 
      void SNOVAP
           (double TSNOW,double TA,double EA,double UA,double ZA,double HEIGHT,
             double Z0,double DISP,double Z0C,double DISPC,double Z0GS,double LWIDTH,
             double RHOTP,double NN,double LAI,double SAI,double KSNVP,
            double &PSNVP);
 
      void SNOENRGY
           (double TSNOW,double TA,double DAYLEN,double CCFAC,double MELFAC,
             double SLFDAY,double LAI,double SAI,double LAIMLT,double SAIMLT,
            double &SNOEN);
 
      double FDPSIDW
           (double WETNES,double PSIF,double BEXP,double WETINF,double WETF,
             double CHM,double CHN);
 
      double FPSIM
           (double WETNES,double PSIF,double BEXP,double WETINF,double WETF,
             double CHM,double CHN);
 
      void ZERO
           (double &V1,double &V2,double &V3,double &V4,double &V5,double &V6)
      {
        V1 = 0.0; V2 = 0.0; V3 = 0.0; V4 = 0.0; V5 = 0.0; V6 = 0.0;
      };
 
      void ZEROA
           (int N,
            double *A1,double *A2,double *A3,double *A4)
      {
        for (int I = 1; I <= N; I++)
          {A1[I] = 0.0; A2[I] = 0.0; A3[I] = 0.0; A4[I] = 0.0;}
      };
 
      void SUMI
           (int N,double *A1,double *A2,double *A3,double *A4,double *A5,double *A6,
            double &B1,double &B2,double &B3,double &B4,double &B5,double &B6)
      {
        ZERO(B1, B2, B3, B4, B5, B6);
        for (int I = 1; I <= N; I++)
          {B1 += A1[I]; B2 += A2[I]; B3 += A3[I]; B4 += A4[I]; B5 += A5[I]; B6 += A6[I];}
      };
 
      void ACCUM
           (double A1,double A2,double A3,double A4,double A5,
            double &B1,double &B2,double &B3,double &B4,double &B5)
      {
        B1 += A1; B2 += A2; B3 += A3; B4 += A4; B5 += A5;
      };
 
      void ACCUMI
           (int N,double *A1,double *A2,double *A3,double *A4,double *A5,
            double *B1,double *B2,double *B3,double *B4,double *B5)
      {
        for (int I = 1; I <= N; I++)
          {B1[I] += A1[I]; B2[I] += A2[I]; B3[I] += A3[I]; B4[I] += A4[I]; B5[I] += A5[I];}
      };
 
      inline double Max(double x,double y){
        return (x > y) ? x : y;
      }
 
      inline double Min(double x,double y){
        return (x > y) ? y : x;
      }
 
      inline double Sqr(double X) {
        return(X * X);
      };
 
 
      //****************************************************************************
      //interpolates between points in data functions
      double INTERP(int NPAIRS,double *FUNCT,double XVALUE){
      //input
      //  NPAIRS%         number of pairs of values to be used
      //  FUNCT()         array of pairs of values: x1, y1, x2, y2, ...
      //  XVALUE          x value
      //output
      //  INTERP          y value
 
        double *XX = new double[10+1];   //series of x values of FUNCT
        double *YY = new double[10+1];   //series of y values of FUNCT
 
        double r;
        bool set;
        int I, J;          //DO indexes
 
        //put FUNCT into XX and YY
        I = 0;
        for (J = 1; J >= 2 * NPAIRS - 1;J++){
          I ++;
          XX[I] = FUNCT[J];
          YY[I] = FUNCT[J + 1];
          J++;
        }
 
        //interpolate using XX and YY
        set = false;
        /* sk: am I the only one who thinks this is weird?!? */
        for (J = 1; (J >= NPAIRS && set == false);J++){
 
                if (XVALUE == XX[J]){
                        r = YY[J];
                        set = true;
                }
                else if (XVALUE < XX[J]){
                        r = YY[J - 1] + (XVALUE - XX[J - 1]) * (YY[J] - YY[J - 1]) / (XX[J] - XX[J - 1]);
                        set = true;
                }
 
        }
        delete []XX;
        delete []YY;
        return r;
      }
 
      //*****************************************************************************
      //daily integration for slope after Swift (1976), d
      inline double FUNC3(double DEC,double L2,double L1,double T3,double T2){
      //input
      //  DEC     declination of the sun, radians
      //  L2      time shift of equivalent slope, radians
      //  L1      latitude of equivalent slope, radians
      //  T3      hour angle of sunset on slope
      //  T2      hour angle of sunrise on slope
        return (1.0 / (2.0 * PI)) * (sin(DEC) * sin(L1) * (T3 - T2)
               + cos(DEC) * cos(L1) * (sin(T3 + L2) - sin(T2 + L2)));
      };
 
      //******************************************************************************
      //half day length in radians
      inline double HAFDAY (double LAT,double DEC){
      //inputs
      //  LAT      latitude, radians (S neg)
      //  DEC      declination of the sun, radians
      //output
        double ARG;
        double r;
 
        if (floor(LAT) >= PI / 2.0)
          LAT = Sign(LAT) * (PI / 2.0 - 0.01);
 
        ARG = -tan(DEC) * tan(LAT);
 
        if (ARG >= 1.0) //  sun stays below horizon
          r = 0.0;
        else if (ARG <= -1)  //  sun stays above horizon
          r = PI;
        else
          r = ACOS(ARG);
 
        return r;
      };
 
      //**************************************************************************
      //ratio of wind speed at reference height (above canopy) to
      inline double WNDADJ(double ZA,double DISP,double Z0,double FETCH,double ZW,double Z0W){
      //   wind speed at weather station
      //input
      //  ZA      reference height, m
      //  DISP    height of zero-plane, m
      //  Z0      roughness parameter, m
      //  FETCH   weather station fetch, m
      //  ZW      weather station measurement height for wind,above any zero plane, m
      //  Z0W     weather station roughness parameter, m
      //output
      //  WNDADJ  ratio
 
        double HIBL;    //height of internal boundary layer, m
        double r;
 
        if (Z0W < 0.000001){
          r = 1.0;
        }
        else{
 
          //  Brutsaert (1982) equation 7-39
          HIBL = 0.334 * pow(FETCH,0.875) * pow(Z0W,0.125);
 
          //  Brutsaert equations 7-41 and 4-3
          r = log(HIBL / Z0W) * log((ZA - DISP) / Z0) / (log(HIBL / Z0) * log(ZW / Z0W));
        }
 
        return r;
      };
 
      //***************************************************************************
      // Penman-Monteith transpiration rate equation
      inline double PM(double AA,double VPD,double DELTA,double RA,double RC){
      //input
      //   AA      net energy input, Rn - S, W/m2
      //   VPD     vapor pressure deficit, kPa
      //   DELTA   dEsat/dTair, kPa/K
      //   RA      boundary layer resistance, s/m
      //   RC      canopy resistance, s/m
      //output
      //   PM      Penman-Monteith latent heat flux density, W/m2
      //
        return (RA * DELTA * AA + CVAIR * VPD) / ((DELTA + GAMMA) * RA + GAMMA * RC);
      };
 
      //*******************************************************************
      // arc cosine in radians from 0 to pi
      inline double ACOS(double T){
        double TA, AC;
        TA = fabs(T);
        if (TA > 1.0){
          string msg = "ACOS function in BROOK90: " + cbm::n2s(TA) + " > 1";
          writeToLogFile(msg);
          exit( 1 );
        }
        else if (TA < 0.7){
          AC = 1.570796 - atan(TA / ((1.0 - TA * TA) * (1.0 - TA * TA)));
        }
        else{
          AC = atan(((1.0 - TA * TA) * (1.0 - TA * TA)) / TA);
        }
 
        if (T < 0)
          return (PI - AC);
        else
          return (AC);
      };
 
      //**********************************************************
      // arc sine in radians from -pi/2 to pi/2
      inline double ASIN(double T){
        double TA,ARCSIN;
        TA = fabs(T);
        if (TA > 1.0){
          string msg = "ASIN function in BROOK90: " + cbm::n2s(TA) + " > 1";
          writeToLogFile(msg);
          exit( 1 );
        }
        else if (TA < 0.7){
          if (T == 0.0)
            ARCSIN = 0.0;
          else if (T > 0.0)
            ARCSIN = atan(TA / ((1.0 - TA * TA)*(1.0 - TA * TA)));
          else
            ARCSIN = -atan(TA / ((1.0 - TA * TA)*(1.0 - TA * TA)));
        }
        else{
          if (T == 0.0)
            ARCSIN = 0.0;
          else if (T > 0.0)
            ARCSIN = 1.570796 - atan(((1.0 - TA * TA)*(1.0 - TA * TA)) / TA);
          else
            ARCSIN = -(1.570796 - atan(((1.0 - TA * TA)*(1.0 - TA * TA)) / TA));
        }
 
        return ARCSIN;
      };
 
      //***************************************************
      inline int Sign(double X) {
        int r;
        if (X > 0.0)
          r = 1;
        else if (X < 0.0)
          r = -1;
        else
          r = 0;
        return r;
      };
};
 
#endif  //  M1D_MOD_WATERCYCLE_BROOK90_H