/*
 * SUMMARY:      InitTerrainMaps() - Initialize terrain coverages
 * 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: Thu Dec 12 10:59:19 2002 by Scott Waichler <waichler@tuff.pnl.gov>
 * DESCRIPTION:  Initialize terrain coverages
 * DESCRIP-END.
 * FUNCTIONS:    InitTerrainMaps()
 *               InitTopoMap()
 *               InitSoilMap()
 *               InitVegMap()
 * COMMENTS:     
 */

#ifndef lint
static char vcid[] = "$Id: InitTerrainMaps.c,v 1.2 2001/03/31 00:10:28 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 "getinit.h"
#include "sizeofnt.h"
#include "slopeaspect.h"
#include "varid.h"

extern void sortElev( int numCells, float *elev, int **row_col ); 

/*****************************************************************************
  InitTerrainMaps()
*****************************************************************************/
void InitTerrainMaps(LISTPTR Input, OPTIONSTRUCT *Options, MAPSIZE *Map, 
		     LAYER *Soil, TOPOPIX ***TopoMap, SOILPIX ***SoilMap, 
		     VEGPIX ***VegMap)
{
  printf("Initializing terrain maps\n");

  InitTopoMap(Input, Options, Map, TopoMap);
  InitSoilMap(Input, Map, Soil, *TopoMap, SoilMap);
  InitVegMap(Input, Map, VegMap);
}

/*****************************************************************************
  InitTopoMap()
*****************************************************************************/
void InitTopoMap(LISTPTR Input, OPTIONSTRUCT *Options, MAPSIZE *Map, 
		 TOPOPIX ***TopoMap)
{
  const char *Routine = "InitTopoMap";
  char VarName[BUFSIZE+1];	/* Variable name */
  char FlowDirFileName[BUFSIZE+1];	/* File with flow direction */
  int i;			/* Counter */
  int x;			/* Counter */
  int y, n;			/* Counter */
  int NumberType;		/* Number type of data set */
  unsigned char *Mask = NULL;   /* Basin mask */
  float *Elev;                  /* Surface elevation */
  short int **Coords;                  /* Coordinates of cells within basin (y,x) */
  int k; /* Counter */
  int **index;

  STRINIENTRY StrEnv[] = {
    {"TERRAIN", "DEM FILE"        , ""  , ""},
    {"TERRAIN", "BASIN MASK FILE", ""  , ""},
    {"TERRAIN", "FLOW DIRECTION FILE", ""  , ""},
    {NULL       , NULL            , ""  , NULL}
  };
  
  /* Process the [TERRAIN] section in the input file */

  /* don't look for flow direction file name if we are not reading 
     flow directions */
  if( Options->ReadFlowDir != TRUE )StrEnv[2].SectionName = NULL;

  /* Read the key-entry pairs from the input file */
  for (i = 0; StrEnv[i].SectionName; i++) {
    GetInitString(StrEnv[i].SectionName, StrEnv[i].KeyName, StrEnv[i].Default,
                  StrEnv[i].VarStr, (unsigned long) BUFSIZE, Input);
    if (IsEmptyStr(StrEnv[i].VarStr))
      ReportError(StrEnv[i].KeyName, 51);
  }


  /* assign name to flow direction file if required */
  if( Options->ReadFlowDir == TRUE ) {
    if (IsEmptyStr(StrEnv[flowdir_file].VarStr))
      ReportError(StrEnv[flowdir_file].KeyName, 51);
    else
      strcpy(FlowDirFileName, StrEnv[flowdir_file].VarStr);
  }
  else
      strcpy(FlowDirFileName, "none");

  /* Read the elevation data from the DEM dataset */
  GetVarName(001, 0, VarName);
  GetVarNumberType(001, &NumberType);
  if (!(Elev = (float *) calloc(Map->NX * Map->NY, 
				SizeOfNumberType(NumberType)))) 
    ReportError((char *) Routine, 1);
  Read2DMatrix(StrEnv[demfile].VarStr, Elev, NumberType, Map->NY, Map->NX, 0,
	       VarName);
  
  /* Read the mask */
  GetVarName(002, 0, VarName);
  GetVarNumberType(002, &NumberType);
  if (!(Mask = (unsigned char *) calloc(Map->NX * Map->NY, 
					SizeOfNumberType(NumberType)))) 
    ReportError((char *) Routine, 1);    
  Read2DMatrix(StrEnv[maskfile].VarStr, Mask, NumberType, Map->NY, Map->NX, 0,
	       VarName);

  /* Assign the attributes to the correct map pixel */
  if (!(*TopoMap = (TOPOPIX **) calloc(Map->NY, sizeof(TOPOPIX *))))
    ReportError((char *) Routine, 1);
  for (y = 0; y < Map->NY; y++) {
    if (!((*TopoMap)[y] = (TOPOPIX *) calloc(Map->NX, sizeof(TOPOPIX))))
      ReportError((char *) Routine, 1);
  }
  for (y = 0, i = 0; y < Map->NY; y++) 
    for (x = 0; x < Map->NX; x++, i++) 
      (*TopoMap)[y][x].Dem  = Elev[i]; 
  free(Elev);

  for (y = 0, i = 0; y < Map->NY; y++)
    for (x = 0; x < Map->NX; x++, i++) 
      (*TopoMap)[y][x].Mask   = Mask[i];
  free(Mask);

  /* set cell area to Map.DX*MAP.DY as default for raster, square grid, 
     cell-to-cell routing this may be reset in InitElevBand()
     initialize remaining variables  */
  Map->NumCells = 0;
  for (y = 0; y < Map->NY; y++)
    for (x = 0; x < Map->NX; x++) {
      (*TopoMap)[y][x].Area = Map->DX * Map->DY;
      (*TopoMap)[y][x].TotalDir = 0;
      (*TopoMap)[y][x].Travel = 0;
      (*TopoMap)[y][x].Grad = 0.0;
      (*TopoMap)[y][x].Slope = 0.0;
      (*TopoMap)[y][x].Aspect = 0.0;
      (*TopoMap)[y][x].FlowGrad = 0.0;
      for (n = 0; n < NDIRS; n++)(*TopoMap)[y][x].Dir[n] = 0;

/* Begin modifications to create ordered cell coordinates.  SRW 10/02. */

      /* Count the number of cells in the basin.  Need this to allocate memory for
         the new, smaller Elev[] and Coords[][].  */
      if (InBasin((*TopoMap)[y][x].Mask)) Map->NumCells++;
    }

  /* Create arrays to hold elevations of only those cells
     within the basin; and y,x of those cells.  Elev is now smaller. */
  GetVarNumberType(001, &NumberType);
  if (!(Elev = (float *) calloc(Map->NumCells, SizeOfNumberType(NumberType)))) 
    ReportError((char *) Routine, 1);

  if (!(Coords = (short int **) calloc(Map->NumCells, sizeof(short int *))))
    ReportError((char *) Routine, 1);
  for (y = 0; y < Map->NumCells; y++) {
    if (!((Coords)[y] = (short int *) calloc(2, sizeof(short int))))
      ReportError((char *) Routine, 1);
  }

  if (!(Map->OrderedCellsRowCol = (int **) calloc(Map->NumCells, sizeof(int *))))
    ReportError((char *) Routine, 1);
  for (y = 0; y < Map->NumCells; y++) {
    if (!((Map->OrderedCellsRowCol)[y] = (int *) calloc(2, sizeof(int))))
      ReportError((char *) Routine, 1);
  }

  k = 0;
  for (y = 0; y < Map->NY; y++)
    for (x = 0; x < Map->NX; x++) {
      /* Save the elevation, y, and x in temporary arrays. */
      if (InBasin((*TopoMap)[y][x].Mask)) {
        Elev[k] = (*TopoMap)[y][x].Dem;
        Coords[k][0] = y;
        Coords[k][1] = x;
        k++;
      }
    }

  /* Sort Elev in descending order--get the indices back.  */
  if( !( index = (int **)calloc( Map->NumCells, sizeof(int *) ) ) )
		fprintf(stderr, "Could not allocate memory in first stage\n");
  if( !( index[0] = (int *)calloc( Map->NumCells * 2, sizeof(int) ) ) )
		fprintf(stderr, "Could not allocate memory in second stage\n");

  for( i = 1; i < Map->NumCells; i++)
		index[i] = index[0] + i * 2;
  for( i = 0; i < Map->NumCells; i++ ){
    index[i][0] = i;
    index[i][1] = (Map->NumCells - 1) - i;
  }

  sortElev( Map->NumCells, Elev, index);

  i = 0;
  for( k = (Map->NumCells-1); k >= 0; k-- ){
    Map->OrderedCellsRowCol[i][0] = Coords[ index[k][0] ][0];
    Map->OrderedCellsRowCol[i][1] = Coords[ index[k][0] ][1];
    i++;
  }

  free(Elev);
  free(Coords);
/* End of modifications to create ordered cell coordinates.  SRW 10/02 */


  /* Calculate slope, aspect, magnitude of subsurface flow gradient, and 
     fraction of flow flowing in each direction based on the land surface 
     slope.  The slope is also used in the inline radiation calculations. */
  if (Options->FlowGradient == TOPOGRAPHY || Options->RadType == INLINE || 
      Options->MM5 == TRUE) 
    ElevationSlopeAspect(Map, *TopoMap, Options->ReadFlowDir, FlowDirFileName);

  /* After calculating the slopes and aspects for all the points, reset the mask 
     if the model is to be run in point mode */

  if (Options->Extent == POINT) {
    for (y = 0; y < Map->NY; y++) 
      for (x = 0; x < Map->NX; x++)
        (*TopoMap)[y][x].Mask = OUTSIDEBASIN;
    (*TopoMap)[Options->PointY][Options->PointX].Mask = (1 != OUTSIDEBASIN);
  }

}

/*****************************************************************************
  InitSoilMap()
*****************************************************************************/
void InitSoilMap(LISTPTR Input, MAPSIZE *Map, LAYER *Soil, 
		 TOPOPIX **TopoMap, SOILPIX ***SoilMap)
{
  const char *Routine = "InitSoilMap";
  char VarName[BUFSIZE+1];	/* Variable name */
  int i;			/* counter */
  int x;			/* counter */
  int y;			/* counter */
  int NumberType;		/* number type */
  short int *Type;          /* Soil type */
  float *Depth;			/* Soil depth */
  STRINIENTRY StrEnv[] = {
    {"SOILS", "SOIL MAP FILE"  , ""  , ""},
    {"SOILS", "SOIL DEPTH FILE", ""  , ""},
    {NULL       , NULL         , ""  , NULL}
  };
  float SoilDepthMultiplier;     /* SRW 10-24-01 */
  char SectionName[] = "SOILS";  /* SRW 10-24-01 */
  char NewStr[BUFSIZE+1];        /* SRW 10-24-01 */

  /* Process the filenames in the [SOILS] section in the input file */

  /* Read the key-entry pairs from the input file */
  for (i = 0; StrEnv[i].SectionName; i++) {
    GetInitString(StrEnv[i].SectionName, StrEnv[i].KeyName, StrEnv[i].Default,
                  StrEnv[i].VarStr, (unsigned long) BUFSIZE, Input);
    if (IsEmptyStr(StrEnv[i].VarStr))
      ReportError(StrEnv[i].KeyName, 51);
  }

  /* Add this to give more flexibility in calibration, SRW 10-24-01.
  Get the multiplier to be applied to all soil depths in the binary file */
  GetInitString(SectionName, "SOIL DEPTH MULTIPLIER", "", NewStr, 
		(unsigned long) BUFSIZE, Input);
  if (!CopyFloat(&SoilDepthMultiplier, NewStr, 1))
    ReportError("SOIL DEPTH MULTIPLIER", 51);
  
  /* Read the soil type */
  GetVarName(003, 0, VarName);
  GetVarNumberType(003, &NumberType);
  if (!(Type = (short int *) calloc(Map->NX * Map->NY, 
                                        SizeOfNumberType(NumberType))))
    ReportError((char *) Routine, 1);
  Read2DMatrix(StrEnv[soiltype_file].VarStr, Type, NumberType, Map->NY, Map->NX, 
	       0, VarName); 

  /* Read the total soil depth  */
  GetVarName(004, 0, VarName);
  GetVarNumberType(004, &NumberType);
  if (!(Depth = (float *) calloc(Map->NX * Map->NY, 
                                 SizeOfNumberType(NumberType)))) 
    ReportError((char *) Routine, 1);
  Read2DMatrix(StrEnv[soildepth_file].VarStr, Depth, NumberType, Map->NY,
	       Map->NX, 0, VarName); 

  /* Assign the attributes to the correct map pixel */
  if (!(*SoilMap = (SOILPIX **) calloc(Map->NY, sizeof(SOILPIX *))))
    ReportError((char *) Routine, 1);
  for (y = 0; y < Map->NY; y++) {
    if (!((*SoilMap)[y] = (SOILPIX *) calloc(Map->NX, sizeof(SOILPIX))))
      ReportError((char *) Routine, 1);
  }
  for (y = 0, i = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++, i++) {
      if (((int)Type[i]) > Soil->NTypes) 
	ReportError(StrEnv[soiltype_file].VarStr, 32);
      (*SoilMap)[y][x].Soil = Type[i];
      (*SoilMap)[y][x].Depth = Depth[i] * SoilDepthMultiplier;  /* modified SRW 10-24-01 */
  
      /* allocate memory for the number of root layers, plus an additional 
	 layer below the deepest root layer */
      if (InBasin(TopoMap[y][x].Mask)) {
        if (!((*SoilMap)[y][x].Moist = 
              (float *) calloc((Soil->NLayers[Type[i]-1]+1), sizeof(float))))
          ReportError((char *) Routine, 1);
        if (!((*SoilMap)[y][x].Perc = 
              (float *) calloc(Soil->NLayers[Type[i]-1], sizeof(float))))
          ReportError((char *) Routine, 1);
        if (!((*SoilMap)[y][x].Temp = 
              (float *) calloc(Soil->NLayers[Type[i]-1], sizeof(float))))
          ReportError((char *) Routine, 1);
      }
      else {
        (*SoilMap)[y][x].Moist = NULL;
        (*SoilMap)[y][x].Perc  = NULL;
        (*SoilMap)[y][x].Temp  = NULL;
      }
    }
  }
  free(Type);
  free(Depth);
}

/*****************************************************************************
  InitVegMap()
*****************************************************************************/
void InitVegMap(LISTPTR Input, MAPSIZE *Map, VEGPIX ***VegMap)
{
  const char *Routine = "InitVegMap";
  char VarName[BUFSIZE+1];
  char VegMapFileName[BUFSIZE+1];
  int i;			/* counter */
  int x;			/* counter */
  int y;			/* counter */
  int NumberType;		/* number type */
  short int *Type;          /* Vegetation type */

  /* Get the map filename from the [VEGETATION] section */
  GetInitString("VEGETATION", "VEGETATION MAP FILE", "", VegMapFileName,
		(unsigned long) BUFSIZE, Input);
  if (!VegMapFileName)
    ReportError("VEGETATION MAP FILE", 51);

  /* Read the vegetation type */
  GetVarName(005, 0, VarName);
  GetVarNumberType(005, &NumberType);
  if (!(Type = (short int *) calloc(Map->NX * Map->NY, 
                                        SizeOfNumberType(NumberType))))
    ReportError((char *) Routine, 1);
  Read2DMatrix(VegMapFileName, Type, NumberType, Map->NY, Map->NX, 0, VarName);

  /* Assign the attributes to the correct map pixel */
  if (!(*VegMap = (VEGPIX **) calloc(Map->NY, sizeof(VEGPIX *))))
    ReportError((char *) Routine, 1);
  for (y = 0; y < Map->NY; y++) {
    if (!((*VegMap)[y] = (VEGPIX *) calloc(Map->NX, sizeof(VEGPIX))))
      ReportError((char *) Routine, 1);
  }
  for (y = 0, i = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++, i++) {
      (*VegMap)[y][x].Veg = Type[i];
      (*VegMap)[y][x].Tcanopy = 0.0;
    }
  }
  free(Type);
}
