/*
 * SUMMARY:      InitModelState.c - Initialize the model state variables
 * 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: Mon Dec  2 13:40:23 2002 by Scott Waichler <waichler@tuff.pnl.gov>
 * DESCRIPTION:  Initialize the model state variables using initial conditions
 *               or a saved state from an earlier model run
 * DESCRIP-END.
 * FUNCTIONS:    InitModelState()
 * COMMENTS:     
 */

#ifndef lint
static char vcid[] = "$Id: InitModelState.c,v 1.2 2001/03/31 00:10:27 waichler Exp $";
#endif /* lint */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "settings.h"
#include "data.h"
#include "DHSVMerror.h"
#include "fileio.h"
#include "functions.h"
#include "constants.h"
#include "sizeofnt.h"
#include "soilmoisture.h"
#include "varid.h"

/*****************************************************************************
  Function name: InitModelState()

  Purpose      : Initialize the state of the model using initial conditions
                 or a saved state from an earlier model run

  Required     :

  Returns      : void

  Modifies     :

  Comments : 
    Initialize the model state, by reading the state variables from a series
    of files.  This allows restarts of the model from any timestep for which
    the model state is known.  These model states can be stored using the
    routine StoreModelState().  Timesteps at which to dump the model state
    can be specified in the file with dump information.

*****************************************************************************/
void InitModelState(DATE *Start, MAPSIZE *Map, OPTIONSTRUCT *Options,
		    PRECIPPIX **PrecipMap, SNOWPIX **SnowMap, 
		    SOILPIX **SoilMap, LAYER Soil, SOILTABLE *SType, 
		    VEGPIX **VegMap, LAYER Veg, VEGTABLE *VType, char *Path, 
		    SNOWTABLE *SnowAlbedo, TOPOPIX **TopoMap, 
		    ROADSTRUCT **Network, UNITHYDRINFO *HydrographInfo,
		    float *Hydrograph)
{
  const char *Routine = "InitModelState";
  char Str[NAMESIZE+1];
  char FileName[NAMESIZE+1];
  FILE *HydroStateFile;
  int i;			/* counter */
  int x;			/* counter */
  int y;			/* counter */
  int NSet;			/* Number of dataset to be read */
  int NSoil;			/* Number of soil layers for current pixel */
  int NVeg;			/* Number of veg layers for current pixel */
  void *Array;
  MAPDUMP DMap;			/* Dump Info */
  float TotalRootDepth = 0;
  float remove;

  printf("Restoring model state\n");

  /* Restore canopy interception */

  NSet = 0;
  
  if (DEBUG)
    printf("Restoring canopy conditions\n");

  sprintf(Str, "%02d.%02d.%02d.%02d.%02d.%02d", Start->Month, Start->Day,
	  Start->Year, Start->Hour, Start->Min, Start->Sec); 
  sprintf(FileName, "%sInterception.State.%s%s", Path, Str, fileext);

  DMap.ID = 202;
  DMap.Layer = 0;
  DMap.Resolution = MAP_OUTPUT;
  strcpy(DMap.FileName, "");
  GetVarAttr(&DMap);
  if (!(Array = calloc(Map->NY * Map->NX, SizeOfNumberType(DMap.NumberType))))
    ReportError((char *) Routine, 1);    
  
  for (i = 0; i < Veg.MaxLayers; i++) {
    
    DMap.ID = 202;
    DMap.Layer = i;
    DMap.Resolution = MAP_OUTPUT;
    strcpy(DMap.FileName, "");
    GetVarAttr(&DMap);
    
    Read2DMatrix(FileName, Array, DMap.NumberType, Map->NY, Map->NX, NSet++,
		 DMap.Name);
    
    for (y = 0; y < Map->NY; y++) {
      for (x = 0; x < Map->NX; x++) {
        if (InBasin(TopoMap[y][x].Mask)) {
          NVeg = Veg.NLayers[(VegMap[y][x].Veg-1)];
          if (i < NVeg)
            PrecipMap[y][x].IntRain[i] = ((float *)Array)[y*Map->NX + x]; 
        }
      }
    }
  }

  for (i = 0; i < Veg.MaxLayers; i++) {

    DMap.ID = 203;
    DMap.Layer = i;
    DMap.Resolution = MAP_OUTPUT;
    strcpy(DMap.FileName, "");
    GetVarAttr(&DMap);

    Read2DMatrix(FileName, Array, DMap.NumberType, Map->NY, Map->NX, NSet++,
		 DMap.Name); 
    
    for (y = 0; y < Map->NY; y++) {
      for (x = 0; x < Map->NX; x++) {
        if (InBasin(TopoMap[y][x].Mask)) {
          NVeg = Veg.NLayers[(VegMap[y][x].Veg-1)];
          if (i < NVeg)
            PrecipMap[y][x].IntSnow[i] = ((float *)Array)[y*Map->NX + x]; 
        }
      }
    }
  }
  
  DMap.ID = 204;
  DMap.Resolution = MAP_OUTPUT;
  strcpy(DMap.FileName, "");
  GetVarAttr(&DMap);
  Read2DMatrix(FileName, Array, DMap.NumberType, Map->NY, Map->NX, NSet++,
	       DMap.Name);
  
  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) {
      if (InBasin(TopoMap[y][x].Mask))
        PrecipMap[y][x].TempIntStorage = ((float *)Array)[y*Map->NX + x]; 
    }
  }
  
  free(Array);

  /* Restore snow pack conditions */

  NSet = 0;

  if (DEBUG)
    printf("Restoring snow pack conditions\n");

  sprintf(FileName,"%sSnow.State.%s%s", Path, Str, fileext);

  DMap.ID = 401;
  DMap.Resolution = MAP_OUTPUT;
  strcpy(DMap.FileName, "");
  GetVarAttr(&DMap);

  if (!(Array = (float *) calloc(Map->NY * Map->NX, 
                                 SizeOfNumberType(DMap.NumberType))))
    ReportError((char *) Routine, 1);    
  
  Read2DMatrix(FileName, Array, DMap.NumberType, Map->NY, Map->NX, NSet++,
	       DMap.Name);

  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) { 
      if (InBasin(TopoMap[y][x].Mask)) {
        SnowMap[y][x].HasSnow = (unsigned char)((float *)Array)[y*Map->NX + x]; 
      }
    }
  }

  DMap.ID = 403;
  DMap.Resolution = MAP_OUTPUT;
  strcpy(DMap.FileName, "");
  GetVarAttr(&DMap);

  Read2DMatrix(FileName, Array, DMap.NumberType, Map->NY, Map->NX, NSet++,
	       DMap.Name);

  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) {
      if (InBasin(TopoMap[y][x].Mask)) {
      SnowMap[y][x].LastSnow = (unsigned short)((float *)Array)[y*Map->NX + x]; 
      }
    }
  }

  DMap.ID = 404;
  DMap.Resolution = MAP_OUTPUT;
  strcpy(DMap.FileName, "");
  GetVarAttr(&DMap);

  Read2DMatrix(FileName, Array, DMap.NumberType, Map->NY, Map->NX, NSet++,
	       DMap.Name);
  
  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) { 
      if (InBasin(TopoMap[y][x].Mask)) {
        SnowMap[y][x].Swq = ((float *)Array)[y*Map->NX + x]; 
      }
    }
  }

  DMap.ID = 406;
  DMap.Resolution = MAP_OUTPUT;
  strcpy(DMap.FileName, "");
  GetVarAttr(&DMap);

  Read2DMatrix(FileName, Array, DMap.NumberType, Map->NY, Map->NX, NSet++,
	       DMap.Name);
  
  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) {
      if (InBasin(TopoMap[y][x].Mask)) {
        SnowMap[y][x].PackWater = ((float *)Array)[y*Map->NX + x]; 
      }
    }
  }

  DMap.ID = 407;
  DMap.Resolution = MAP_OUTPUT;
  strcpy(DMap.FileName, "");
  GetVarAttr(&DMap);

  Read2DMatrix(FileName, Array, DMap.NumberType, Map->NY, Map->NX, NSet++,
	       DMap.Name);
  
  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) { 
      if (InBasin(TopoMap[y][x].Mask)) {
      SnowMap[y][x].TPack = ((float *)Array)[y*Map->NX + x]; 
      }
    }
  }
  
  DMap.ID = 408;
  DMap.Resolution = MAP_OUTPUT;
  strcpy(DMap.FileName, "");
  GetVarAttr(&DMap);

  Read2DMatrix(FileName, Array, DMap.NumberType, Map->NY, Map->NX, NSet++,
	       DMap.Name);
  
  for (y = 0; y < Map->NY; y++) { 
    for (x = 0; x < Map->NX; x++) {
      if (InBasin(TopoMap[y][x].Mask)) {      
        SnowMap[y][x].SurfWater = ((float *)Array)[y*Map->NX + x]; 
      }
    }
  }

  DMap.ID = 409;
  DMap.Resolution = MAP_OUTPUT;
  strcpy(DMap.FileName, "");
  GetVarAttr(&DMap);

  Read2DMatrix(FileName, Array, DMap.NumberType, Map->NY, Map->NX, NSet++,
	       DMap.Name);
  
  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) {
      if (InBasin(TopoMap[y][x].Mask)) {
        SnowMap[y][x].TSurf = ((float *)Array)[y*Map->NX + x]; 
      }
    }
  }
  
  DMap.ID = 410;
  DMap.Resolution = MAP_OUTPUT;
  strcpy(DMap.FileName, "");
  GetVarAttr(&DMap);

  Read2DMatrix(FileName, Array, DMap.NumberType, Map->NY, Map->NX, NSet++,
	       DMap.Name);
  
  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) { 
      if (InBasin(TopoMap[y][x].Mask)) {
        SnowMap[y][x].ColdContent = ((float *)Array)[y*Map->NX + x]; 
      }
    }
  }

  free(Array);

  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) { 
      if (InBasin(TopoMap[y][x].Mask)) {
        if (SnowMap[y][x].HasSnow)
          SnowMap[y][x].Albedo = CalcSnowAlbedo(SnowMap[y][x].TSurf, 
                                                SnowMap[y][x].LastSnow,
                                                SnowAlbedo);
        else
          SnowMap[y][x].Albedo = 0;
      }
    }
  }
  
  /* Restore soil conditions */

  NSet = 0;

  if (DEBUG)
    printf("Restoring soil conditions\n");

  sprintf(FileName, "%sSoil.State.%s%s", Path, Str, fileext);
  
  DMap.ID = 501;
  DMap.Layer = 0;
  DMap.Resolution = MAP_OUTPUT;
  strcpy(DMap.FileName, "");
  GetVarAttr(&DMap);

  if (!(Array = (float *) calloc(Map->NY * Map->NX, 
                                 SizeOfNumberType(DMap.NumberType))))
    ReportError((char *) Routine, 1);    
  
  for (i = 0; i < Soil.MaxLayers+1; i++) {

    DMap.ID = 501;
    DMap.Layer = i;
    DMap.Resolution = MAP_OUTPUT;
    strcpy(DMap.FileName, "");
    GetVarAttr(&DMap);

    Read2DMatrix(FileName, Array, DMap.NumberType, Map->NY, Map->NX, NSet++,
		 DMap.Name);
  
    for (y = 0; y < Map->NY; y++) {
      for (x = 0; x < Map->NX; x++) {
        if (InBasin(TopoMap[y][x].Mask)) {
          NSoil = Soil.NLayers[(SoilMap[y][x].Soil-1)];
          if (i <= NSoil)
	    SoilMap[y][x].Moist[i] = 
	      ((float *)Array)[y*Map->NX + x];
	   
        }
      }
    }
  }

  DMap.ID = 505;
  DMap.Resolution = MAP_OUTPUT;
  strcpy(DMap.FileName, "");
  GetVarAttr(&DMap);

  Read2DMatrix(FileName, Array, DMap.NumberType, Map->NY, Map->NX, NSet++,
	       DMap.Name);

  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) {
      if (InBasin(TopoMap[y][x].Mask)) {
      SoilMap[y][x].TSurf = ((float *)Array)[y*Map->NX + x]; 
      }
    }
  }
  
  for (i = 0; i < Soil.MaxLayers; i++) {

    DMap.ID = 511;
    DMap.Layer = i;
    DMap.Resolution = MAP_OUTPUT;
    strcpy(DMap.FileName, "");
    GetVarAttr(&DMap);

    Read2DMatrix(FileName, Array, DMap.NumberType, Map->NY, Map->NX, NSet++,
		 DMap.Name); 
  
    for (y = 0; y < Map->NY; y++) {
      for (x = 0; x < Map->NX; x++) {
        if (InBasin(TopoMap[y][x].Mask)) {
          NSoil = Soil.NLayers[(SoilMap[y][x].Soil-1)];
          if (i < NSoil)
            SoilMap[y][x].Temp[i] = ((float *)Array)[y*Map->NX + x]; 
        }
      }
    }
  }
  
  DMap.ID = 510;
  DMap.Resolution = MAP_OUTPUT;
  strcpy(DMap.FileName, "");
  GetVarAttr(&DMap);

  Read2DMatrix(FileName, Array, DMap.NumberType, Map->NY, Map->NX, NSet++,
	       DMap.Name);
  
  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) { 
      if (InBasin(TopoMap[y][x].Mask)) {
        SoilMap[y][x].Qst = ((float *)Array)[y*Map->NX + x]; 
      }
    }
  }

  DMap.ID = 512;
  DMap.Resolution = MAP_OUTPUT;
  strcpy(DMap.FileName, "");
  GetVarAttr(&DMap);

  Read2DMatrix(FileName, Array, DMap.NumberType, Map->NY, Map->NX, NSet++,
	       DMap.Name); 
  
  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) {
      if (InBasin(TopoMap[y][x].Mask)) {
        SoilMap[y][x].Runoff = ((float *)Array)[y*Map->NX + x]; 
      }
    }
  }

  free(Array);

  /* Calculate the water table depth at each point based on the soil 
     moisture profile. Give an error message if the water ponds on the
     surface since that should not be allowed at this point */
  remove=0.0;

  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) {

      /* SatFlow needs to be initialized properly in the future.  For now it 
         will just be set to zero here */
      SoilMap[y][x].SatFlow = 0.0;
      if (InBasin(TopoMap[y][x].Mask)) {

        /* make sure depth at bottom of lowest root zone < soil depth */
        TotalRootDepth = 0;
        for (i = 0; i < (Soil.NLayers[SoilMap[y][x].Soil-1]); i++) 
           TotalRootDepth += VType[VegMap[y][x].Veg-1].RootDepth[i];
        if( TotalRootDepth > SoilMap[y][x].Depth ) {
          printf("Total depth of root zone > soil depth in cell y: %d x: %d\n", y, x );
          exit(1);
        }


        if ((SoilMap[y][x].TableDepth = 
             WaterTableDepth((Soil.NLayers[SoilMap[y][x].Soil-1]),
                             SoilMap[y][x].Depth,
                             VType[VegMap[y][x].Veg-1].RootDepth,
                             SType[SoilMap[y][x].Soil-1].Porosity,
                             SType[SoilMap[y][x].Soil-1].FCap,
                             Network[y][x].Adjust, SoilMap[y][x].Moist)) 
            < 0.0)  /*  ReportError((char *) Routine, 35);*/{
	      remove-=SoilMap[y][x].TableDepth*Map->DX*Map->DY;
	      SoilMap[y][x].TableDepth = 0.0;
	    }
        SoilMap[y][x].SatThickness = 0;
        SoilMap[y][x].CumGwRecharge = 0;
        SoilMap[y][x].CumChannelLoss = 0;
      }
      else {
        SoilMap[y][x].SatThickness = 0;
        SoilMap[y][x].TableDepth = 0;
        SoilMap[y][x].CumGwRecharge = 0;
        SoilMap[y][x].CumChannelLoss = 0;
      }
    }
  }
  if (remove>0.0){
  printf("WARNING:excess water in soil profile is %f m^3 \n",remove);
  printf("Expect possible large flood wave during first timesteps \n");
  }
  
  /* If the unit hydrograph is used for flow routing, initialize the unit 
     hydrograph array */

  if (Options->Extent == BASIN && Options->HasNetwork == FALSE) {
    sprintf(FileName, "%sHydrograph.State.%s", Path, Str);
    OpenFile(&HydroStateFile, FileName, "r", FALSE);
    for (i = 0; i < HydrographInfo->TotalWaveLength; i++)
      fscanf(HydroStateFile, "%f\n", &(Hydrograph[i]));
    fclose(HydroStateFile);
  }
}
