/*
 * SUMMARY:      InitNewMonth.c - Initialize new time periods
 * 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: Wed May  2 13:06:15 2001 by Scott Waichler <markw@r101483.pnl.gov>
 * DESCRIPTION:  Initialization functions that have to be executed at the
 *               beginning of certain timesteps
 * DESCRIP-END.
 * FUNCTIONS:    InitNewMonth()
 *               InitNewDay()
 *               InitNewStep()
 * COMMENTS:     
 */

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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "settings.h"
#include "data.h"
#include "DHSVMerror.h"
#include "functions.h"
#include "constants.h"
#include "init.h"
#include "fifobin.h"
#include "fileio.h"
#include "rad.h"
#include "slopeaspect.h"
#include "sizeofnt.h"
#include "inittime.h"
#include "readweights.h"

/*****************************************************************************
  InitNewMonth()
  At the start of a new month, read the new radiation files 
  (diffuse and direct beam), and potentially a new LAI value.
*****************************************************************************/
void InitNewMonth(TIMESTRUCT *Time, OPTIONSTRUCT *Options, MAPSIZE *Map, TOPOPIX **TopoMap, 
                  float **PrismMap, int RadType, RADCLASSPIX** RadMap, RADTABLE *RadClass, 
                  INPUTFILES *InFiles, int NVegs, VEGTABLE *VType,
                  int NStats, METLOCATION *Stat, char *Path,
                  unshort ****MetWeights, WEIGHT_INFO **WeightInfo, WEIGHT_FILES WeightFiles)
{
  const char *Routine = "InitNewMonth";
  char FileName[MAXSTRING+1];
  char str[MAXSTRING+1];
  FILES TableFile;
  int i;
  int j;
  int y,x;
  int NumberType;
  float *Array = NULL; 
  uchar **BasinMask;

  if (DEBUG)
    printf("Initializing new month\n");

  if( Options->Interpolation == READ ) {
    printf("Reading new met station weights from file\n");

    if (!(BasinMask = (uchar **) calloc(Map->NY, sizeof(uchar *))))
      ReportError((char *) Routine, 1);
    for (y = 0; y < Map->NY; y++) {
      if (!(BasinMask[y] = (uchar *) calloc(Map->NX, sizeof(uchar))))
	ReportError((char *) Routine, 1);
    }
    
    for (y = 0; y < Map->NY; y++)
      for (x = 0; x < Map->NX; x++) 
	    BasinMask[y][x] = TopoMap[y][x].Mask;

    ReadWeights(Map->NX, Map->NY, BasinMask, MetWeights, WeightInfo, 
        Time->Current.Month, WeightFiles.Path);

    free( BasinMask );
  }


  /* If PRISM precipitation fields are being used to interpolate the 
     observed precipitation fields, then read in the new months field */

  if (Options->Prism == TRUE) {
    printf("reading in new PRISM field for month %d \n",Time->Current.Month);

    sprintf(str, ".%02d", Time->Current.Month);
    sprintf(FileName, "%s%s%s%s", Path,"PrismMap", str,".bin");


    if (!(Array = (float *) calloc(Map->NY * Map->NX, sizeof(float))))
    ReportError((char *) Routine, 1);    
    NumberType = NC_FLOAT;
  
    Read2DMatrixBin(FileName, Array, NumberType, Map->NY, Map->NX, 0);

    printf("done reading matrix \n");

  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) {
      PrismMap[y][x] = Array[y*Map->NX + x]; 
    }
  }
  
  free(Array);
   
  }

  /* If the short wave radiation is calculated with the aid of IPW generated
     radiation images then a new radiation image and table needs to be read
     at the beginning of each month */

  if (RadType == IPW) {

    sprintf(str, ".%02d", Time->Current.Month);
    sprintf(TableFile.FileName, "%s%s%s", InFiles->RadTablePath, str, fileext);
    OpenFile(&(TableFile.FilePtr), TableFile.FileName, "r", FALSE);
    
    ReadRadTable(&TableFile, Time->NDaySteps, RadClass);

    fclose(TableFile.FilePtr);

    sprintf(str, ".%02d", Time->Current.Month); 
    sprintf(FileName, "%s%s%s", InFiles->RadMapPath, str, fileext);
    ReadRadMap(Map, RadClass->End, Time->NDaySteps, NStats, FileName, 
               TopoMap, RadMap, Stat);
  }

}

/*****************************************************************************
  Function name: InitNewDay()

  Purpose      : Initialize the Earth-Sun geometry variables at the beginning
                 of each day

  Required     : 
    int DayOfYear           - day of year (January 1 = 1)
    SOLARGEOMETRY *SolarGeo - structure with information about Earth-Sun 
                              geometry

  Returns      : void

  Modifies     : 
    SOLARGEOMETRY *SolarGeo

  Comments     : To be excuted at the beginning of each new day
*****************************************************************************/
void InitNewDay(int DayOfYear, SOLARGEOMETRY *SolarGeo, int NVegs, VEGTABLE *VType, DATE Current)
{
  int DaysPerMonth[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
  float fract;
  int i, j, month;
  month = Current.Month;

  if( month == 2 ){
    if (IsLeapYear(Current.Year))
      DaysPerMonth[1] = 29;
    else
      DaysPerMonth[1] = 28;
  }

  /* update LAI */
  fract = (float) (Current.Day - 1) / DaysPerMonth[ month - 1 ];
  /* printf("fract %f\n", fract ); */
  
  
  for (i = 0; i < NVegs; i++) {
    for (j = 0; j < VType[i].NVegLayers; j++) {
      
      if( month == 12 ){
        VType[i].LAI[j] = (1.0 - fract)*VType[i].MonthLAI[j][11] 
          + fract*VType[i].MonthLAI[j][0];
      }
      else {
        VType[i].LAI[j] = (1.0 - fract)*VType[i].MonthLAI[j][month-1] 
          + fract*VType[i].MonthLAI[j][month];        
      }
      
      VType[i].MaxInt[j] = VType[i].LAI[j] * VType[i].Fract[j] *
        LAI_WATER_MULTIPLIER;
      /* printf("Change LAI: VegType %d Layer %d Month %d Day %d LAI %f\n", 
        i, j, Current.Month, Current.Day, VType[i].LAI[j]); */
    }
  }

  

  SolarDay(DayOfYear, SolarGeo->Longitude, SolarGeo->Latitude,
           SolarGeo->StandardMeridian, &(SolarGeo->NoonHour), 
           &(SolarGeo->Declination), &(SolarGeo->HalfDayLength),
           &(SolarGeo->Sunrise), &(SolarGeo->Sunset),
           &(SolarGeo->TimeAdjustment), &(SolarGeo->SunEarthDistance));
}

/*****************************************************************************
  Function name: InitNewStep()

  Purpose      : Initialize Earth-Sun geometry and meteorological data at the
                 beginning of each timestep

  Required     :
    MAPSIZE Map              - Structure with information about location
    TIMESTRUCT Time          - Structure with time information
    int RadType              - Method by which to calculate radiation at each
                               pixel 
    int PrecipType           - Type of precipitation input, RADAR, STATION or
                               OROGRAPHIC
    int FlowGradient         - Type of FlowGradient calculation
    int NStats               - Number of meteorological stations
    METLOCATION *Stat        - Structure with information about the
                               meteorological stations in or near the study
                               area 
    char *RadarFileName      - Name of file with radar images
    MAPSIZE Radar            - Structure with information about the
                               precipitation radar coverage
    RADCLASSPIX **RadMap     - Structure with radiation data for each pixel
    RADARPIX **RadarMap      - Structure with precipitation information for
                               each radar pixel
    RADTABLE RadClass        - Radiation levels for IPW generated radiation
                               classes 
    SOLARGEOMETRY *SolarGeo  - structure with information about Earth-Sun 
                               geometry
    SOILPIX **SoilMap        - structure with soil information
    float ***MM5Input        - MM5 input maps
    float ***WindModel       - Wind model maps
                           
  Returns      : void

  Modifies     :

  Comments     : To be executed at the beginning of each time step
*****************************************************************************/
void InitNewStep(INPUTFILES *InFiles, MAPSIZE *Map, TIMESTRUCT *Time, 
                 int NSoilLayers, OPTIONSTRUCT *Options, int NStats, 
                 METLOCATION *Stat, char *RadarFileName, MAPSIZE *Radar, 
                 RADTABLE *RadClass, RADARPIX **RadarMap, 
                 SOLARGEOMETRY *SolarGeo, TOPOPIX **TopoMap, 
                 RADCLASSPIX **RadMap, SOILPIX **SoilMap,
                 float ***MM5Input, float ***WindModel) 
{
  const char *Routine = "InitNewStep";
  float Slope;
  int i;                        /* counter */
  int j;                        /* counter */
  int x;                        /* counter */
  int y;                        /* counter */
  int NumberType;               /* number type in MM5 input */
  int Step;                     /* Step in the MM5 Input */
  float *Array = NULL;
  float test = 0.0;

  if (Options->MM5 == TRUE) {

    /* Read the data from the MM5 files */
    
    if (!(Array = (float *) calloc(Map->NY * Map->NX, sizeof(float))))
      ReportError((char *) Routine, 1);    
    NumberType = NC_FLOAT;

    Step = NumberOfSteps(&(Time->StartMM5), &(Time->Current), Time->Dt);
    
    Read2DMatrixBin(InFiles->MM5Temp, Array, NumberType, Map->NY, Map->NX, Step);
    for (y = 0; y < Map->NY; y++)
      for (x = 0; x < Map->NX; x++)
        MM5Input[MM5_temperature-1][y][x] = Array[y*Map->NX+x];

    Read2DMatrixBin(InFiles->MM5Humidity, Array, NumberType, Map->NY, Map->NX,
		 Step); 
    for (y = 0; y < Map->NY; y++)
      for (x = 0; x < Map->NX; x++)
        MM5Input[MM5_humidity-1][y][x] = Array[y*Map->NX+x];

    Read2DMatrixBin(InFiles->MM5Wind, Array, NumberType, Map->NY, Map->NX, Step);
    for (y = 0; y < Map->NY; y++)
      for (x = 0; x < Map->NX; x++)
        MM5Input[MM5_wind-1][y][x] = Array[y*Map->NX+x];

    Read2DMatrixBin(InFiles->MM5ShortWave, Array, NumberType, Map->NY, Map->NX,
		 Step); 
    for (y = 0; y < Map->NY; y++)
      for (x = 0; x < Map->NX; x++)
        MM5Input[MM5_shortwave-1][y][x] = Array[y*Map->NX+x];

    Read2DMatrixBin(InFiles->MM5LongWave, Array, NumberType, Map->NY, Map->NX,
		 Step); 
    for (y = 0; y < Map->NY; y++)
      for (x = 0; x < Map->NX; x++)
        MM5Input[MM5_longwave-1][y][x] = Array[y*Map->NX+x];

    Read2DMatrixBin(InFiles->MM5Pressure, Array, NumberType, Map->NY, Map->NX,
		 Step); 
    for (y = 0; y < Map->NY; y++)
      for (x = 0; x < Map->NX; x++)
        MM5Input[MM5_pressure-1][y][x] = Array[y*Map->NX+x];

    Read2DMatrixBin(InFiles->MM5Precipitation, Array, NumberType, Map->NY, Map->NX,
		 Step);
    for (y = 0; y < Map->NY; y++)
      for (x = 0; x < Map->NX; x++)
        MM5Input[MM5_precip-1][y][x] = Array[y*Map->NX+x];

    for (i = 0, j = MM5_precip; i < NSoilLayers; i++, j++) {
      Read2DMatrixBin(InFiles->MM5SoilTemp[i], Array, NumberType, Map->NY, Map->NX,
		   Step);
      for (y = 0; y < Map->NY; y++)
	for (x = 0; x < Map->NX; x++)
	  MM5Input[j][y][x] = Array[y*Map->NX+x];
    }
    free(Array);
  }
  else {
    /* Calculate variables related to the position of the sun above the
       horizon */
    SolarHour(SolarGeo->Latitude, 
	      (Time->DayStep+1)*((float)Time->Dt)/SECPHOUR, 
	      ((float)Time->Dt)/SECPHOUR, SolarGeo->NoonHour,
	      SolarGeo->Declination, SolarGeo->Sunrise, SolarGeo->Sunset,
	      SolarGeo->TimeAdjustment, SolarGeo->SunEarthDistance,
	      &(SolarGeo->SineSolarAltitude), &(SolarGeo->DayLight),
	      &(SolarGeo->SolarTimeStep), &(SolarGeo->SunMax),
	      &(SolarGeo->SolarAzimuth)); 
    
    /* Calculate the radiation received by each pixel corrected for topographic
       effects.  The resulting values are corrected for atmospheric effects in
       MakeLocalMetData() */
    
    if (Options->RadType == IPW) {
      for (y = 0; y < Map->NY; y++) {
        for (x = 0; x < Map->NX; x++) {
          if (InBasin(TopoMap[y][x].Mask)) {
            RadMap[y][x].Beam = 
              RadClass->Beam[Time->DayStep][(int)RadMap[y][x].BeamCode[Time->DayStep]];
            RadMap[y][x].Diffuse = 
              RadClass->Diffuse[Time->DayStep][(int)RadMap[y][x].DiffuseCode[Time->DayStep]];
          }
        }
      }
    }
    else if (Options->RadType == INLINE) {
      for (y = 0; y < Map->NY; y++) {
        for (x = 0; x < Map->NX; x++) {
          if (InBasin(TopoMap[y][x].Mask)) {
            
            Slope = atan(TopoMap[y][x].Slope);
            
            /* WORK IN PROGRESS:  Currently a constant albedo of ALBEDO (see
               rad.h) is assumed.  This should be replaced by some albedo that
               varies with each pixel */
           

            SolarAngle(SolarGeo->Latitude, 0.15, SolarGeo->Declination, 
                       TopoMap[y][x].Aspect, Slope, SolarGeo->SunMax, 
                       SolarGeo->SineSolarAltitude, SolarGeo->SolarTimeStep, 
  		               SolarGeo->DayLight, SolarGeo->SolarAzimuth,
		               ((float)Time->Dt)/3600., &(RadMap[y][x].Beam), 
                       &(RadMap[y][x].Diffuse));

            /*SolarAngle(SolarGeo->Latitude, ALBEDO, SolarGeo->Declination, 
                       TopoMap[y][x].Aspect, Slope, SolarGeo->SunMax, 
                       SolarGeo->SineSolarAltitude, SolarGeo->SolarTimeStep, 
  		               SolarGeo->DayLight, SolarGeo->SolarAzimuth,
		               ((float)Time->Dt)/SECPHOUR, &(RadMap[y][x].Beam), 
                       &(RadMap[y][x].Diffuse));*/
          }
        }
      }
    }
    else
      ReportError((char *) Routine, 43);
    
    /* also do the met stations */
    for (i = 0; i < NStats; i++) {
      y = Stat[i].Loc.N;
      x = Stat[i].Loc.E;
      
      if (Options->RadType == IPW) {
        RadMap[y][x].Beam = 
          RadClass->Beam[Time->DayStep][(int)Stat[i].Data.BeamCode[Time->DayStep]];
        RadMap[y][x].Diffuse = 
          RadClass->Diffuse[Time->DayStep][(int)Stat[i].Data.DiffuseCode[Time->DayStep]]; 
      }
      else if (Options->RadType == INLINE) {

	/* The slope for the stations themselves is set to zero, because
           radiation is normally measured (or generated) for a flat surface.
           This means that the radiation density for the pixel that contains the
           station is NOT necessarily the same as the radiation received by the
           station.  The two are ONLY identical if the station pixel is on a
           flat slope as well. */

	Slope = 0.0;
	if((Stat[i].Loc.N > Map->NY || Stat[i].Loc.N < 0 || 
	    Stat[i].Loc.E > Map->NX || Stat[i].Loc.E < 0)){ /*outside bounding box*/
	  if(SolarGeo->DayLight==TRUE){
	  Stat[i].Data.SinBeamMod = SolarGeo->SunMax;
	  Stat[i].Data.SinDiffuseMod = SolarGeo->SunMax;;
	  }
	  else
	    {
	  Stat[i].Data.SinBeamMod = 0.0;
	  Stat[i].Data.SinDiffuseMod = 1.0;
	    }
	}
	else{/*inside bounding box*/
        SolarAngle(SolarGeo->Latitude, ALBEDO, SolarGeo->Declination, 
                   TopoMap[y][x].Aspect, Slope, SolarGeo->SunMax, 
                   SolarGeo->SineSolarAltitude,
                   SolarGeo->SolarTimeStep, SolarGeo->DayLight,
                   SolarGeo->SolarAzimuth, ((float)Time->Dt)/SECPHOUR, 
		   &(RadMap[y][x].Beam), &(RadMap[y][x].Diffuse));
      Stat[i].Data.SinBeamMod = RadMap[y][x].Beam;
      Stat[i].Data.SinDiffuseMod = RadMap[y][x].Diffuse;
	}
      }  
    }
    
    
    /* if the flow gradient is based on the water table, recalculate the water
       table gradients and the flow directions */
    if (Options->FlowGradient == WATERTABLE) {
      /* Calculate the WaterLebel, i.e.e tha height of the water table above 
         some datum */
      for (y = 0; y < Map->NY; y++) {
        for (x = 0; x < Map->NX; x++) {
          if (InBasin(TopoMap[y][x].Mask)) {
            SoilMap[y][x].WaterLevel = 
              TopoMap[y][x].Dem - SoilMap[y][x].TableDepth;
          }
        }
      }
      HeadSlopeAspect(Map, TopoMap, SoilMap);
    }
    
    GetMetData(Options, Time, NSoilLayers, NStats, SolarGeo->SunMax, Stat, 
               Radar, RadarMap, RadarFileName); 

  }
}
