/*
 * SUMMARY:      MakeLocalMetData.c - Generates meteorological conditions
 * USAGE:        Part of DHSVM
 *
 * AUTHOR:       Bart Nijssen
 * ORG:          University of Washington, Department of Civil Engineering
 * E-MAIL:       nijssen@u.washington.edu
 * ORIG-DATE:    Apr-96
 * Last Change: Tue Dec 17 08:18:35 2002 by Scott Waichler <waichler@tuff.pnl.gov>
 * DESCRIPTION:  Generates meteorological conditions for each individual cell
 * DESCRIP-END.
 * FUNCTIONS:    MakeLocalMetData()
 * COMMENTS:     
 */

#ifndef lint
static char vcid[] = "$Id: MakeLocalMetData.c,v 1.4 2001/05/15 17:30:26 waichler Exp $";
#endif /* lint */

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "settings.h"
#include "data.h"
#include "snow.h"
#include "DHSVMerror.h"
#include "functions.h"
#include "constants.h"
#include "rad.h"
#include "svp.h"

/*****************************************************************************
  Function name: MakeLocalMetData()

  Purpose      : Generates meteorological for each individual cell

  Required     :
    int y 
    int x
    MAPSIZE Map
    int DayStep
    unsigned char PrecipType
    int NStats
    METLOCATION *Stat
    unshort *MetWeights
    float LocalElev
    RADCLASSPIX *RadMap 
    PRECIPPIX *PrecipMap
    MAPSIZE Radar
    RADARPIX **RadarMap

  Returns      :
    PIXMET LocalMet

  Modifies     :

  Comments     :
    Reference: Shuttleworth, W.J., Evaporation,  In: Maidment, D. R. (ed.),
               Handbook of hydrology,  1993, McGraw-Hill, New York, etc..
*****************************************************************************/
PIXMET MakeLocalMetData(int y, int x, MAPSIZE *Map, int DayStep, 
                        OPTIONSTRUCT *Options, int NStats, 
                        METLOCATION *Stat, unshort *MetWeights, 
                        float LocalElev, RADCLASSPIX *RadMap, 
                        PRECIPPIX *PrecipMap, MAPSIZE *Radar, 
                        RADARPIX **RadarMap, float **PrismMap, SNOWPIX *LocalSnow,
                        SNOWTABLE *SnowAlbedo, float ***MM5Input,
                        float ***WindModel, float **PrecipLapseMap, 
      MET_MAP_PIX ***MetMap, int NGraphics, int Month)  
{
  float BeamMultiplier;         /* Multiplier for beam radiation to include
                                   atmospheric effects*/
  float CurrentWeight;    /* weight for current station */
  float DiffuseMultiplier;      /* Multiplier for diffuse radiation to include
                                   atmospheric effects*/
  float ScaleWind = 1;    /* Wind to be scaled by model factors if 
                                   WindSource == MODEL */
  float Temp;     /* Temporary variable */
  float WeightSum;    /* sum of the weights */
  int i;      /* counter */
  int RadarX;                   /* X coordinate of radar map coordinate */
  int RadarY;                   /* Y coordinate of radar map coordinate */
  float TempLapseRate;
  int WindDirection = 0;  /* Direction of model wind */
  PIXMET LocalMet;    /* local met data */

  LocalMet.Tair = 0.0;
  LocalMet.Rh   = 0.0;
  LocalMet.Wind = 0.0;
  LocalMet.Sin  = 0.0;
  LocalMet.SinBeam    = 0.0;
  LocalMet.SinDiffuse = 0.0;
  LocalMet.Lin  = 0.0;
  TempLapseRate = 0.0;


  if (Options->MM5 == TRUE) {
    LocalMet.Tair   = MM5Input[MM5_temperature-1][y][x];
    LocalMet.Rh     = MM5Input[MM5_humidity-1]   [y][x];
    LocalMet.Wind   = MM5Input[MM5_wind-1]       [y][x];
    LocalMet.Sin    = MM5Input[MM5_shortwave-1]  [y][x];
    LocalMet.Lin    = MM5Input[MM5_longwave-1]   [y][x];
    LocalMet.Press  = MM5Input[MM5_pressure-1]   [y][x];
    PrecipMap->Precip = MM5Input[MM5_precip-1]   [y][x];
  }
  else /*MM5 is false */{  
    WeightSum = 0.0;
    for (i = 0; i < NStats; i++) {
      WeightSum += (float) MetWeights[i];
      if (Options->WindSource == MODEL && Stat[i].IsWindModelLocation) {
        ScaleWind = Stat[i].Data.Wind;
        WindDirection = Stat[i].Data.WindDirection;
      }
    }

#ifndef SHOW_MET_ONLY
    for (i = 0; i < NStats; i++) {
      CurrentWeight  = ((float)MetWeights[i])/WeightSum;
      LocalMet.Tair  += CurrentWeight * 
                      LapseT(Stat[i].Data.Tair, 
                        Stat[i].Elev, LocalElev, Stat[i].NTempLapse, 
                        Stat[i].Data.TempLapse[0], Stat[i].Data.TempLapse[1],
                        Stat[i].ElevBreak[0], Options->PrecipLapse);
                        
                        
                        /*LapseT(Stat[i].Data.Tair, Stat[i].Elev, LocalElev, 
                               Stat[i].Data.TempLapse);*/
      LocalMet.Rh   += CurrentWeight * Stat[i].Data.Rh;
      if (Options->WindSource == STATION)
        LocalMet.Wind += CurrentWeight * Stat[i].Data.Wind;
      LocalMet.Lin  += CurrentWeight * Stat[i].Data.Lin;
      TempLapseRate += CurrentWeight * Stat[i].Data.TempLapse[0];
    }
    if (Options->WindSource == MODEL)
      LocalMet.Wind = ScaleWind * WindModel[WindDirection-1][y][x];
    
    if (RadMap->Beam > 0.0 || RadMap->Diffuse > 0.0) {
      
      BeamMultiplier = 0.0;
      DiffuseMultiplier = 0.0;
      for (i = 0; i < NStats; i++) {
        CurrentWeight = ((float) MetWeights[i])/WeightSum;
        BeamMultiplier += CurrentWeight * Stat[i].Data.BeamRatio;
        DiffuseMultiplier += CurrentWeight * Stat[i].Data.DiffuseRatio; 
      }
      
      /*     printf("%g\t%g\t%g\t%g\t", RadMap->Beam, RadMap->Diffuse,  */
      /*            BeamMultiplier, DiffuseMultiplier); */
      RadMap->Beam *= BeamMultiplier;
      RadMap->Diffuse *= DiffuseMultiplier;
      /*     printf("%g\t%g\n", RadMap->Beam, RadMap->Diffuse); */
    }
    
    LocalMet.Sin = RadMap->Beam + RadMap->Diffuse;
    LocalMet.SinBeam = RadMap->Beam;
    LocalMet.SinDiffuse = RadMap->Diffuse;

    
    if (Options->PrecipType == RADAR) {
      RadarY = (int)((y+Radar->OffsetY)*Map->DY/Radar->DY);
      RadarX = (int)((x-Radar->OffsetX)*Map->DX/Radar->DX);
      PrecipMap->Precip = RadarMap[RadarY][RadarX].Precip;
    }
    /* WORK IN PROGRESS, taken from old DHSVM version */
    /* Air pressure */
    /* In rare cases - i.e. when the lapse rate has a different sign for 
       different met stations - you can end up with a TemplapseRate of 0.0
       This will result in a crash, so a check was put in (Jul 28, 1997 - Bart
       Nijssen).  It is somewhat awkward to interpolate lapse rates anyway, so
       a better way of doing this would be welcome */
    if (TempLapseRate != 0.0) {
      Temp = 9.8067/(TempLapseRate * 287.0);
      LocalMet.Press = 101300. * pow(((288.0-TempLapseRate*LocalElev)/288.0), 
            Temp);
    }
    else
      LocalMet.Press = 101300.;
#endif
  }

  if(Options->QPF==TRUE || Options->MM5==FALSE){
    if (Options->PrecipType == STATION && Options->Prism == FALSE) {
      PrecipMap->Precip = 0.0;
      for (i = 0; i < NStats; i++) {
        
        if( Options->Interpolation == READ ) /* WeightSum can exceed 1.0 if reading weights */
          CurrentWeight = ((float) MetWeights[i]) / (float)PRISM_SCALE_FACTOR;
        else
          CurrentWeight = ((float) MetWeights[i]) / WeightSum;

        if (Options->PrecipLapse == MAP)
          PrecipMap->Precip += CurrentWeight * 
          LapsePrecip(Stat[i].Data.Precip, 
                      Stat[i].Elev, LocalElev, Stat[i].NPrecipLapse, 
                      Stat[i].Data.PrecipLapse[0], Stat[i].Data.PrecipLapse[1],
                      Stat[i].ElevBreak[0], Options->PrecipLapse);
        else
          PrecipMap->Precip += CurrentWeight *
          LapsePrecip(Stat[i].Data.Precip, 
                      Stat[i].Elev, LocalElev, Stat[i].NPrecipLapse, 
                      Stat[i].Data.PrecipLapse[0], Stat[i].Data.PrecipLapse[1],
                      Stat[i].ElevBreak[0], Options->PrecipLapse);
      }
    }
    else if (Options->PrecipType == STATION && Options->Prism == TRUE) {
      PrecipMap->Precip = 0.0;
      for (i = 0; i < NStats; i++) {
        CurrentWeight = ((float) MetWeights[i])/WeightSum;
    /* this is the real prism interpolation*/
    /* note that X = position from left  boundary, ie # of columns*/
    /* note that Y = position from upper boundary, ie # of rows   */
        if(Options->Outside==FALSE)
          PrecipMap->Precip += 
          CurrentWeight *Stat[i].Data.Precip/
          PrismMap[Stat[i].Loc.N][Stat[i].Loc.E]*PrismMap[y][x];
        else
          PrecipMap->Precip += 
          CurrentWeight *Stat[i].Data.Precip/
          Stat[i].PrismPrecip[Month-1]*PrismMap[y][x];
      }
    }
  }

  /* due to the nature of the interpolation scheme in DHSVM and the */
  /* interpolation scheme to handle the mess of different formats of met stations */
  /* in the PRISM project */
  /* relative humidities can be quite low when precip is occuring */
  /* at times this will results in PET being greater than precip */
  /* allow an option in DHSVM to override RH if Precip is occuring */
#ifndef SHOW_MET_ONLY
  if (Options->Rhoverride==TRUE){
    if (PrecipMap->Precip > 0.0) LocalMet.Rh=100.0;
  }

  PrecipMap->CumPrecip += PrecipMap->Precip;  /*  SRW 12/05/02.  */

  /* Separate precipitation into rainfall and snowfall */
  if (PrecipMap->Precip > 0.0 && LocalMet.Tair < MAX_SNOW_TEMP) {
    if (LocalMet.Tair > MIN_RAIN_TEMP) 
      PrecipMap->SnowFall = PrecipMap->Precip * 
        (MAX_SNOW_TEMP - LocalMet.Tair)/(MAX_SNOW_TEMP - MIN_RAIN_TEMP);
    else
      PrecipMap->SnowFall = PrecipMap->Precip;
  }
  else 
    PrecipMap->SnowFall = 0.0;
  
  PrecipMap->RainFall = PrecipMap->Precip - PrecipMap->SnowFall;
  
  /* Local heat of vaporization, Eq. 4.2.1, Shuttleworth (1993) */
  LocalMet.Lv = 2501000 - 2361 * LocalMet.Tair;
  
  /* Psychrometric constant */
  LocalMet.Gamma = CP * LocalMet.Press/(EPS * LocalMet.Lv);
  
  /* Saturated vapor pressure, Eq. 4.2.2, Shuttleworth (1993) */
  LocalMet.Es = SatVaporPressure(LocalMet.Tair);
  
  /* Slope of vapor pressure curve, Eq. 4.2.3, Shuttleworth (1993) */
  LocalMet.Slope = 4098.0 * LocalMet.Es / 
    ((237.3 + LocalMet.Tair)*(237.3 + LocalMet.Tair));
  
  /* Actual vapor pressure */
  LocalMet.Eact = LocalMet.Es*(LocalMet.Rh/100.); 
  
  /* Vapor pressure deficit */
  LocalMet.Vpd = LocalMet.Es - LocalMet.Eact;
  
  /* Air density, Eq. 4.2.4 Shuttleworth (1993) */
  LocalMet.AirDens = 0.003486 * LocalMet.Press / (275 + LocalMet.Tair);
  
  if (LocalSnow->HasSnow) {
    if (PrecipMap->SnowFall > 0.0)
      LocalSnow->LastSnow = 0;
    else
      LocalSnow->LastSnow++;
    LocalSnow->Albedo = CalcSnowAlbedo(LocalSnow->TSurf, LocalSnow->LastSnow, 
                                       SnowAlbedo);
  }
  else 
    LocalSnow->LastSnow = 0;

#endif
  /* SRW 12/16/02. Enforce minimum wind speed of 0.01 m/s to avoid divide by zero in
     aero resistance calculations.  */
  LocalMet.Wind = (LocalMet.Wind > 0.01) ? LocalMet.Wind : 0.01;

  if (NGraphics>0){

    (*MetMap)[y][x].accum_precip=(*MetMap)[y][x].accum_precip+PrecipMap->Precip;
    (*MetMap)[y][x].air_temp=LocalMet.Tair;
    (*MetMap)[y][x].wind_speed=LocalMet.Wind;
    (*MetMap)[y][x].humidity=LocalMet.Rh;
  }

  
  return LocalMet;
}
