/*
 * SUMMARY:      MainDHSVM.c - Distributed Hydrology-Soil-Vegetation Model
 * USAGE:        DHSVM
 *
 * AUTHOR:       Bart Nijssen
 * ORG:          University of Washington, Department of Civil Engineering
 * E-MAIL:       nijssen@u.washington.edu
 * ORIG-DATE:    Apr-96
 * Last Change: Thu Oct  2 15:55:33 2003 by Scott Waichler <waichler@snow.pnl.gov>
 * DESCRIPTION:  Main routine to drive DHSVM, the Distributed 
 *               Hydrology-Soil-Vegetation Model  
 * DESCRIP-END.cd
 * FUNCTIONS:    main()
 * COMMENTS:
 */

#ifndef lint
static char vcid[] = "$Id: MainDHSVM.cpp,v 1.6 2001/06/01 19:29:02 markw Exp $";
#endif /* lint */

/******************************************************************************/
/*				    INCLUDES                                  */
/******************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

extern "C" {
#include "settings.h"
#include "constants.h"
#include "data.h"
#include "DHSVMerror.h"
#include "functions.h"
#include "fileio.h"
#include "getinit.h"
#include "DHSVMChannel.h"
#include "elevationband.h"
#include "calcsolar.h"
#include "inittime.h"
#include "groundwater.h"
#include "frozensoil.h"
}
#include "ftablesdhsvm.h"

/******************************************************************************/
/*				GLOBAL VARIABLES                              */
/******************************************************************************/

/* global function pointers */
void (*CreateMapFile)(char *FileName, ...);
int (*Read2DMatrix)(char *FileName, void *Matrix, int NumberType, int NY,
 		    int NX, int NDataSet, ...);
int (*Write2DMatrix)(char *FileName, void *Matrix, int NumberType, int NY, 
		     int NX, ...); 

/* global strings */
char *version = "$Id: MainDHSVM.cpp,v 1.6 2001/06/01 19:29:02 markw Exp $";		/* store version string */
char commandline[BUFSIZE+1] = ""; /* store command line */
char fileext[BUFSIZ+1] = "";	/* file extension */

/******************************************************************************/
/*				      MAIN                                    */
/******************************************************************************/
int main(int argc, char **argv) 
{
  float *Hydrograph = NULL;
  float ***MM5Input  = NULL;
  float **PrecipLapseMap = NULL;
  float **PrismMap = NULL;
  float ***WindModel = NULL;

  /*  char tempchar[255];
      FILE *outputfile;*/
  /* i is needed for an ugly hack, and needs to be get rid of */

  int i;
  int x;			/* row counter */
  int y;			/* column counter */
  int NStats;			/* Number of meteorological stations */
  int ReadBinFiles = TRUE;  /*  use to see if it is necessary to read *.bin files when in 
                                Options.Extent = ELEVBAND  */
  unshort ***MetWeights = NULL;	/* 3D array with weights for interpolating 
				   meteorological variables between the 
				   stations */

  int NGraphics;                 /* number of graphics for X11*/
  int *which_graphics;           /* which graphics for X11*/
  int TimeLoopIsOver = 0;   /* flag for indicating after last simulation step (1=TRUE) */

  AGGREGATED Total = {             /* Total or average value of a 
                                      variable over the entire basin */
    { 0.0, NULL, NULL, NULL, NULL, 0.0 }, /* EVAPPIX */
    { 0.0, 0.0, 0.0, NULL, NULL, 0.0, 0.0 }, /* PRECIPPIX */
    { {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, 0.0, 0.0, 0.0 }, /* PIXRAD */
    { 0, 0.0, 0, 0.0 },         /* RADCLASSPIX */
    { 0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0,
      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, /* SNOWPIX */
    { 0, 0.0, NULL, NULL, NULL, 0.0, 0.0, 0.0, 0.0,
      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, /*SOILPIX */
    0.0, 0.0, 0.0, 0.0, 0.0, 0l, 0.0, 0.0, 0.0, 0.0, 0.0, 0, 0.0, 0.0, 0.0, 0.0, 0.0 
  };
  CHANNEL ChannelData = { NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };
  DUMPSTRUCT Dump;
  EVAPPIX **EvapMap = NULL;
  INPUTFILES InFiles;
  LAYER Soil;
  LAYER Veg;
  LISTPTR Input = NULL;		/* Linked list with input strings */
  MAPSIZE Map;			/* Size and location of model area */
  MAPSIZE Radar;		/* Size and location of area covered by 
				   precipitation radar*/
  METLOCATION *Stat = NULL; 
  OPTIONSTRUCT Options;         /* Structure with information which program
                                   options to follow */
  PIXMET LocalMet;		/* Meteorological conditions for current pixel
				 */
  PRECIPPIX **PrecipMap = NULL;
  RADARPIX **RadarMap = NULL;
  RADCLASSPIX **RadMap = NULL;
  RADTABLE RadClass;
  ROADSTRUCT **Network = NULL;	/* 2D Array with channel information for each
                                   pixel */ 
  SNOWPIX **SnowMap = NULL;
  MET_MAP_PIX **MetMap = NULL;
  SNOWTABLE *SnowAlbedo = NULL;
  SOILPIX **SoilMap = NULL;
  SOILTABLE *SType = NULL;
  SOLARGEOMETRY SolarGeo;       /* Geometry of Sun-Earth system (needed for
                                   INLINE radiation calculations */
  BANDSTRUCT *ElevBand;
  BANDINFO ElevBandInfo;

  TIMESTRUCT Time;
  TOPOPIX **TopoMap = NULL;
  UNITHYDR **UnitHydrograph = NULL;
  UNITHYDRINFO HydrographInfo;  /* Information about unit hydrograph */
  VEGPIX **VegMap = NULL;
  VEGTABLE *VType = NULL;
  WATERBALANCE Mass =           /* parameter for mass balance calculations */
  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

  WEIGHT_INFO **WeightInfo;
  WEIGHT_FILES WeightFiles;
  STATESTRUCT State;
  GWPIX **Groundwater = NULL;
  FRZPIX **FrozenSoil = NULL;


/*****************************************************************************
  Initialization Procedures 
*****************************************************************************/
  if (argc != 2) {
    fprintf(stderr, "\nUsage: %s inputfile\n\n", argv[0]);
    exit(EXIT_FAILURE);
  }

  sprintf(commandline, "%s %s", argv[0], argv[1]);
  strcpy(InFiles.Const, argv[1]);

  printf("\nRunning DHSVM version: %s\n", version);
  printf("\nSTARTING INITIALIZATION PROCEDURES\n\n");
 
  ReadInitFile(InFiles.Const, &Input);

  InitConstants(Input, &Options, &Map, &SolarGeo, &Time);

  InitFileIO(Options.FileFormat);
  InitTables(Time.NDaySteps, Input, &SType, &Soil, &VType, &Veg, &SnowAlbedo, 
             Options.Infiltration, Options.HydraulicRouting);

  if(Options.Extent == ELEVBAND) {
  	InitElevBand(Input, Map, &ElevBand, &ElevBandInfo);
    ReadBinFiles = ElevBandInfo.ReadBinFiles;
    if( ReadBinFiles == FALSE )
      ElevBandBuildGrids(Map, &TopoMap, &VegMap, &ElevBand, &ElevBandInfo, &SoilMap, 
                         &Soil); 
  }

  if( ReadBinFiles ) {
    InitTerrainMaps(Input, &Options, &Map, &Soil, &TopoMap, &SoilMap, 
		                &VegMap);
  }

  if(Options.Extent == ELEVBAND && ReadBinFiles) {
	  ElevBandUpdate(Map, &TopoMap, &VegMap, &ElevBand, &ElevBandInfo, &SoilMap);
  }

  if(Options.Groundwater)  
     InitGroundwater(Input, &Options, &Map, TopoMap, &Groundwater);

  if(Options.FrozenSoil)  
     InitFrozenSoil(Input, &Options, &Map, &FrozenSoil);

  if (Options.HasNetwork){
    if(Options.ChannelRouting == FTABLE ) InitFtables(Input);
    InitChannel(Input, &Map, Time.Dt, &ChannelData, SoilMap, TopoMap);
  }
  else if (Options.Extent != POINT)
    InitUnitHydrograph(Input, &Map, TopoMap, &UnitHydrograph, 
		       &Hydrograph, &HydrographInfo);
  InitNetwork(Options.HasNetwork, Map.NY, Map.NX, Map.DX, TopoMap,
              SoilMap, VegMap, VType, &Network, &ChannelData);
 
  InitMetSources(Input, &Options, &Map, Soil.MaxLayers, &Time,
		 &InFiles, &NStats, &Stat, &RadClass, &Radar);

  /* the following piece of code is for the UW PRISM project */
  /* for real-time verification of SWE at Snotel sites */
  /* Other users, set OPTION.SNOTEL to FALSE, or use TRUE with caution */

  if(Options.Snotel == TRUE) {
    printf("Warning: All met stations locations are being set to vegetation class 19\n");
    printf("To disable this feature set Snotel OPTION to FALSE\n");
    for (i=0; i< NStats; i++){
      printf("veg type for station %d is %d ",i,VegMap[Stat[i].Loc.N][Stat[i].Loc.E].Veg);
      printf("setting to veg type 19 (assumed bare class)\n");
      VegMap[Stat[i].Loc.N][Stat[i].Loc.E].Veg = 19;
    }
  }

  InitMetMaps(Time.NDaySteps, &Map, &Radar, &Options, InFiles.WindMapPath, 
	      InFiles.PrecipLapseFile, &PrecipLapseMap, &PrismMap, &EvapMap, &PrecipMap, 
	      &RadarMap, &RadMap, SoilMap, &Soil, VegMap, &Veg, TopoMap, 
	      &MM5Input, &WindModel); 
  InitInterpolationWeights(&Map, &Options, TopoMap, &MetWeights, Stat, NStats,
                            Input, &WeightInfo, &WeightFiles, Time.Start.Month);

  InitDump(Input, &Options, &Map, Soil.MaxLayers, Veg.MaxLayers, Time.Dt,
	   TopoMap, &Dump, &NGraphics, &which_graphics, &State); 

  if (Options.HasNetwork == TRUE) {
    InitChannelDump(&ChannelData, Dump.Path);
    if( Options.ReadState == TRUE )
      ReadChannelState(Dump.InitStatePath, &(Time.Start), ChannelData.streams);
    else
      ZeroChannelState(ChannelData.streams);

  }

  InitSnowMap(&Map, &SnowMap);
  InitAggregated(Veg.MaxLayers, Soil.MaxLayers, &Total, Map, TopoMap ); 

  if( Options.ReadState == TRUE ){
    InitModelState(&(Time.Start), &Map, &Options, PrecipMap, SnowMap, SoilMap,
		 Soil, SType, VegMap, Veg, VType, Dump.InitStatePath,
		 SnowAlbedo, TopoMap, Network, &HydrographInfo, Hydrograph);
    
    if(Options.Groundwater)
      ReadGroundwaterState(Dump.InitStatePath, &(Time.Start), &Map,  SoilMap, 
                           TopoMap, Groundwater);
      
  }
  else{
    DefaultModelState(&(Time.Start), &Map, &Options, PrecipMap, SnowMap, SoilMap,
		 Soil, SType, VegMap, Veg, VType, Dump.InitStatePath,
		 SnowAlbedo, TopoMap, Network, &HydrographInfo, Hydrograph, State);

     MakeState(&Map, &(Time.Start), &Options, &Dump, TopoMap, EvapMap,
	     PrecipMap, RadMap, SnowMap, MetMap, VegMap, &Veg, SoilMap, &Soil, &Total,
	     &HydrographInfo, ChannelData.streams, Hydrograph, TimeLoopIsOver,
       Groundwater, FrozenSoil);

  }

  InitNewMonth(&Time, &Options, &Map, TopoMap, PrismMap, Options.RadType, RadMap, &RadClass, 
               &InFiles, Veg.NTypes, VType, NStats, Stat, Dump.InitStatePath,
               &MetWeights, WeightInfo, WeightFiles );

  if(NGraphics>0){
  printf("Initialzing X11 display and graphics \n");
  InitXGraphics(argc,argv,Map.NY,Map.NX,NGraphics,&MetMap);  
  printf("done \n");
  }

  printf("There are %d cells in the basin mask.\n", Map.NumCells);

  /* Done with initialization, delete the list with input strings */
  DeleteList(Input);

  /* setup for mass balance calculations */
  Aggregate(Options.Extent, &Map, &Options, TopoMap, &Soil, &Veg, VegMap, EvapMap, PrecipMap, 
	    RadMap, SnowMap, SoilMap, &Total, VType, Network, Groundwater, FrozenSoil);
  Mass.StartWaterStorage = 
    Total.Runoff + Total.CanopyWater + Total.SoilWater + Total.Snow.Swq +
    Total.Soil.SatFlow + Total.GwDepth;
  Mass.OldWaterStorage = Mass.StartWaterStorage;

  DumpList( &(Dump.List), Soil.MaxLayers, Veg.MaxLayers );

/*****************************************************************************
  Perform Calculations 
*****************************************************************************/

  printf("\nSTARTING CALCULATIONS\n\n");

  while (Before(&(Time.Current), &(Time.End)) || 
	 IsEqualTime(&(Time.Current), &(Time.End))) {
    ResetAggregate(&Soil, &Veg, &Total);
    if (IsNewMonth(&(Time.Current), Time.Dt)) 
  InitNewMonth(&Time, &Options, &Map, TopoMap, PrismMap, Options.RadType, RadMap, &RadClass, 
               &InFiles, Veg.NTypes, VType, NStats, Stat, Dump.InitStatePath,
               &MetWeights, WeightInfo, WeightFiles );

    if (IsNewDay(Time.DayStep)) {
        PrintDate(&(Time.Current), stdout);
	      printf("\n");
      
      InitNewDay(Time.Current.JDay, &SolarGeo, Veg.NTypes, VType, Time.Current);
    
        if(Options.FrozenSoil)  
          UpdateSoilPhase(Time.NDaySteps, Time.Current.JDay, Map, SnowMap, FrozenSoil);    
    }
  
    /*    if(Time.Current.JDay==1 && Time.DayStep == 0){
      sprintf(tempchar,"annual.precip.%d",Time.Current.Year);
      if((outputfile = fopen(tempchar,"wt"))==NULL){
	printf("error opening annual precip output file \n");
	exit(-1);
      }
      for (y = 0; y < Map.NY; y++)
	for (x = 0; x < Map.NX; x++){
	fprintf(outputfile,"%f ",MetMap[y][x].accum_precip);
	MetMap[y][x].accum_precip=0.0;
	}
      fclose(outputfile);
    }
    */
    InitNewStep(&InFiles, &Map, &Time, Soil.MaxLayers, &Options, NStats, Stat, 
                InFiles.RadarFile, &Radar, &RadClass, RadarMap, &SolarGeo, 
                TopoMap, RadMap, SoilMap, MM5Input, WindModel);

    /* initialize channel/road networks for time step */
    
    if (Options.HasNetwork) {
      channel_step_initialize_network(ChannelData.streams);
      channel_step_initialize_network(ChannelData.roads);
    }
  
    for (y = 0; y < Map.NY; y++) {
      for (x = 0; x < Map.NX; x++) {
	if (InBasin(TopoMap[y][x].Mask))  {

   //    MetMap[y][x].accum_precip=0.0; 
      /* SRW 12/04/02 */

    if( Options.Interpolation == READ ) {
       TempFullWeightArray( &WeightFiles, WeightInfo[y][x], MetWeights[y][x], NStats ); 
          
       LocalMet = 
            MakeLocalMetData(y, x, &Map, Time.DayStep, &Options, NStats, Stat, 
                             WeightFiles.WeightVector, TopoMap[y][x].Dem, 
                             &(RadMap[y][x]), &(PrecipMap[y][x]), &Radar, 
                             RadarMap, PrismMap, &(SnowMap[y][x]), SnowAlbedo,
                             MM5Input, WindModel, PrecipLapseMap, &MetMap, 
			                       NGraphics,Time.Current.Month);
    }
    else{
          LocalMet = 
            MakeLocalMetData(y, x, &Map, Time.DayStep, &Options, NStats, Stat, 
                             MetWeights[y][x], TopoMap[y][x].Dem, 
                             &(RadMap[y][x]), &(PrecipMap[y][x]), &Radar, 
                             RadarMap, PrismMap, &(SnowMap[y][x]), SnowAlbedo,
                             MM5Input, WindModel, PrecipLapseMap, &MetMap, 
			                       NGraphics,Time.Current.Month);
    }

    
    if(Options.FrozenSoil) FrozenSoil[y][x].airTemp += LocalMet.Tair;  

    Total.MeanTair +=       /* for basin-mean */
      LocalMet.Tair * (Total.ActiveCells * TopoMap[y][x].Area / Total.BasinArea);  


#ifndef SHOW_MET_ONLY
	  
          for (i = 0; i < Soil.MaxLayers; i++) {
            if (Options.MM5 == TRUE) 
              SoilMap[y][x].Temp[i] = MM5Input[i+N_MM5_MAPS][y][x];
            else
              SoilMap[y][x].Temp[i] = Stat[0].Data.Tsoil[i];
          }
          
      
          MassEnergyBalance(y, x, Map.DX, Time.Dt, Options.HeatFlux,
                            Options.Extent, Veg.MaxLayers, &LocalMet, &(Network[y][x]), 
                            &(PrecipMap[y][x]), &(VType[VegMap[y][x].Veg-1]), 
                            &(VegMap[y][x]), &(SType[SoilMap[y][x].Soil-1]), 
                            &(SoilMap[y][x]), &(SnowMap[y][x]),
                            &(EvapMap[y][x]), &(Total.Rad), Options.Infiltration, 
                            &(Time.Current), FrozenSoil, Options.FrozenSoil); 

#endif	  
	}
      }
    } 
    
#ifndef SHOW_MET_ONLY    
 
    if (Options.Extent != POINT)  /* SRW 1-3-01.  To obtain correct mass balance in POINT mode. */
      RouteSubSurface( Time.Dt, &Map, TopoMap, VType, VegMap, Network,
                    SType, SoilMap, &ChannelData, &Total );

    if (Options.Groundwater)  
      RouteGroundwater( Time.Dt, &Map, TopoMap, SType, SoilMap, 
                   &ChannelData, Groundwater, &(Total.GwToChannel) );

    if (Options.HasNetwork)
    RouteChannel(&ChannelData, &Time, &Map, TopoMap, SoilMap, &Total, Options, SType ); 

    if (Options.Extent == BASIN)
      RouteSurface(&Map, &Time, TopoMap, SoilMap, Options.HasNetwork, 
                   UnitHydrograph, &HydrographInfo, Hydrograph,
                   &(Dump.Stream), Options.HydraulicRouting, SType, &ChannelData );  

#endif
 
    if(NGraphics>0)
    draw(&(Time.Current),IsEqualTime(&(Time.Current),&(Time.Start)),
	 Map.NX,Map.NY,NGraphics,which_graphics,
	 VType,SType,SnowMap,SoilMap,VegMap,TopoMap,PrecipMap,PrismMap,EvapMap,
	 RadMap,MetMap);

#ifndef SHOW_MET_ONLY
                 
    if( Options.Extent == ELEVBAND && ElevBandInfo.PrintOutflows == TRUE )
      PrintElevBand (&(Time.Current), ElevBandInfo, ChannelData.streams, ElevBand);
    
    Aggregate(Options.Extent, &Map, &Options, TopoMap, &Soil, &Veg, VegMap, EvapMap, PrecipMap, 
	      RadMap, SnowMap, SoilMap, &Total, VType, Network, Groundwater, FrozenSoil);
    
    MassBalance(&(Time.Current), &(Dump.Balance), &Total, &Mass);
    
    ExecDump(&Map, &(Time.Current), &Options, &Dump, TopoMap, EvapMap,
	     PrecipMap, RadMap, SnowMap, MetMap, VegMap, &Veg, SoilMap, &Soil, &Total,
	     &HydrographInfo, ChannelData.streams, Hydrograph, TimeLoopIsOver,
       Groundwater, FrozenSoil);
#endif    

    IncreaseTime(&Time);

  }

  TimeLoopIsOver = 1;
  
  ExecDump(&Map, &(Time.Current), &Options, &Dump, TopoMap, EvapMap, PrecipMap,
	   RadMap, SnowMap, MetMap, VegMap, &Veg, SoilMap, &Soil, &Total,
	   &HydrographInfo, ChannelData.streams, Hydrograph, TimeLoopIsOver,
     Groundwater, FrozenSoil);
  
  FinalMassBalance(&(Dump.Balance), &Total, &Mass);


/*****************************************************************************
  Cleanup
*****************************************************************************/
  
  printf("\nSTARTING CLEANUP\n\n");
  
  printf("\nEND OF MODEL RUN\n\n");

  return EXIT_SUCCESS;
 
}
