/*
 * SUMMARY:      InitTables.c - Initialize lookup tables
 * 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-MOD: Tue May  1 15:38:36 2001 by Scott Waichler <markw@r101483.pnl.gov>
 * DESCRIPTION:  Initialize lookup tables
 * DESCRIP-END.
 * FUNCTIONS:    InitTables() 
 *               InitSoilTable() 
 *               InitVegTable() 
 *               InitSnowTable()
 * COMMENTS:     
 */

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

#include <ctype.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "settings.h"
#include "data.h"
#include "DHSVMerror.h"
#include "functions.h"
#include "Calendar.h"
#include "constants.h"
#include "aerodynamic.h"
#include "calckh.h"
#include "fileio.h"
#include "getinit.h"

/*******************************************************************************/
/*				  InitTables()                                 */
/*******************************************************************************/
void InitTables(int StepsPerDay, LISTPTR Input, SOILTABLE **SType, LAYER *Soil,
		VEGTABLE **VType, LAYER *Veg, SNOWTABLE **SnowAlbedo, int infiltrationRate,
     int hydraulicRouting)
{
  printf("Initializing tables\n");

  if ((Soil->NTypes = InitSoilTable(SType, Input, Soil, 
    infiltrationRate, hydraulicRouting)) == 0)
    ReportError("Input Options File", 8);

  if ((Veg->NTypes = InitVegTable(VType, Input, Veg)) == 0)
    ReportError("Input Options File", 8); 

  InitSnowTable(SnowAlbedo, StepsPerDay);
  InitSatVaporTable();
}

/********************************************************************************
  Function Name: InitSoilTable()

  Purpose      : Initialize the soil lookup table
                 Processes most of the following section in InFileName:
		 [SOILS]

  Required     :
    SOILTABLE **SType - Pointer to lookup table
    LISTPTR Input     - Pointer to linked list with input info
    LAYER *Soil       - Pointer to structure with soil layer information

  Returns      : Number of soil layers

  Modifies     : SoilTable and Soil

  Comments     :
********************************************************************************/
int InitSoilTable(SOILTABLE **SType, LISTPTR Input, LAYER *Soil, 
                  int infiltrationRate, int hydraulicRouting)
{
  const char *Routine = "InitSoilTable";
  int i;			/* counter */
  int j;                        /* counter */
  int NSoils;			/* Number of soil types */
  char KeyName[max_ground_melt+1][BUFSIZE+1];
  char *KeyStr[] = {
    "SOIL DESCRIPTION",
    "LATERAL CONDUCTIVITY",
    "EXPONENTIAL DECREASE",
    "MAXIMUM INFILTRATION",
    "SURFACE ALBEDO",
    "NUMBER OF SOIL LAYERS",
    "POROSITY",
    "PORE SIZE DISTRIBUTION",
    "BUBBLING PRESSURE",
    "FIELD CAPACITY",
    "WILTING POINT",
    "BULK DENSITY",
    "VERTICAL CONDUCTIVITY",
    "THERMAL CONDUCTIVITY",
    "THERMAL CAPACITY",
    "INFILTRATION COEFFICIENT",
    "INFILTRATION EXPONENT",
    "QUICK FLOW COEFFICIENT",
    "QUICK FLOW EXPONENT",
    "MAXIMUM SOIL RESISTANCE",
    "SOIL RESISTANCE COEFF",
    "GROUNDWATER CONDUCTIVITY",
    "GROUNDWATER CONDUCTIVITY LAT",
    "BASE LAYER CONDUCTIVITY",
    "MAXIMUM GROUND MELT"
  };
  char SectionName[] = "SOILS";
  char VarStr[max_ground_melt+1][BUFSIZE+1];

  /* Get the number of different soil types */
  GetInitString(SectionName, "NUMBER OF SOIL TYPES", "", VarStr[0], 
		(unsigned long) BUFSIZE, Input);
  if (!CopyInt(&NSoils, VarStr[0], 1))
    ReportError("NUMBER OF SOIL TYPES", 51);

  if (NSoils == 0)
    return NSoils;

  if (!(Soil->NLayers = (int *) calloc(NSoils, sizeof(int))))
    ReportError((char *) Routine, 1);

  if (!(*SType = (SOILTABLE *) calloc(NSoils, sizeof(SOILTABLE))))
    ReportError((char *) Routine, 1);
  
  /********** Read information and allocate memory for each soil type *********/

  Soil->MaxLayers = 0;

  for (i = 0; i < NSoils; i++) {

    /* Read the key-entry pairs from the input file */
    for (j = 0; j <= max_ground_melt; j++) {
      sprintf(KeyName[j], "%s %d", KeyStr[j], i+1);
      GetInitString(SectionName, KeyName[j], "", VarStr[j], 
		    (unsigned long) BUFSIZE, Input);
    }

    /* Assign the entries to the appropriate variables */
    if (IsEmptyStr(VarStr[soil_description]))
      ReportError(KeyName[soil_description], 51);
    strcpy((*SType)[i].Desc, VarStr[soil_description]);
    (*SType)[i].Index = i;

    if (!CopyFloat(&((*SType)[i].KsLat), VarStr[lateral_ks], 1))
      ReportError(KeyName[lateral_ks], 51);

    if (!CopyFloat(&((*SType)[i].KsLatExp), VarStr[exponent], 1))
      ReportError(KeyName[exponent], 51);

    if (!CopyFloat(&((*SType)[i].MaxInfiltrationRate), 
		    VarStr[max_infiltration], 1))
      ReportError(KeyName[max_infiltration], 51);

    if (!CopyFloat(&((*SType)[i].Albedo), VarStr[soil_albedo], 1))
      ReportError(KeyName[soil_albedo], 51);

    if( infiltrationRate == VARIABLE ) { /* infiltration rate */
      if (!CopyFloat(&((*SType)[i].infiltrationCoeff), VarStr[infilt_coeff], 1))
        ReportError(KeyName[infilt_coeff], 51);

      if (!CopyFloat(&((*SType)[i].infiltrationExp), VarStr[infilt_exp], 1))
        ReportError(KeyName[infilt_exp], 51);
    }

    if( hydraulicRouting == TRUE ) { /* quickflow routing */
      if (!CopyFloat(&((*SType)[i].quickflowCoeff), VarStr[quickflow_coeff], 1))
        ReportError(KeyName[quickflow_coeff], 51);

      if (!CopyFloat(&((*SType)[i].quickflowExp), VarStr[quickflow_exp], 1))
        ReportError(KeyName[quickflow_exp], 51);
    }
      
    if (!CopyFloat(&((*SType)[i].maxResistance), VarStr[maximum_soil_resistance], 1))
        ReportError(KeyName[maximum_soil_resistance], 51);

    if (!CopyFloat(&((*SType)[i].resistanceCoeff), VarStr[soil_resistance_coeff], 1))
        ReportError(KeyName[soil_resistance_coeff], 51);

    /* deep groundwater vertical saturated hydraulic conductivity */
    if (!CopyFloat(&((*SType)[i].groundwaterKs), VarStr[groundwater_conductivity], 1))
        ReportError(KeyName[groundwater_conductivity], 51);

    /* deep groundwater lateral saturated hydraulic conductivity */
    if (!CopyFloat(&((*SType)[i].groundwaterKsLat), VarStr[groundwater_conductivity_lat], 1))
        ReportError(KeyName[groundwater_conductivity_lat], 51);

    /* saturated hydraulic conductivity of layer below groundwater zone*/
    if (!CopyFloat(&((*SType)[i].baseLayerKs), VarStr[base_layer_conductivity], 1))
        ReportError(KeyName[base_layer_conductivity], 51);

    /* maximum rate of snowmelt from ground heat flux */
    if (!CopyFloat(&((*SType)[i].maxGroundMelt), VarStr[max_ground_melt], 1))
        ReportError(KeyName[max_ground_melt], 51);

     /* allocate memory for the soil layers */
    if (!CopyInt(&(*SType)[i].NLayers, VarStr[number_of_layers], 1))
      ReportError(KeyName[number_of_layers], 51);
    Soil->NLayers[i] = (*SType)[i].NLayers;
    if (Soil->NLayers[i] > Soil->MaxLayers)
      Soil->MaxLayers = Soil->NLayers[i];
    
    if (!((*SType)[i].Porosity = (float *) calloc((*SType)[i].NLayers, 
						  sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*SType)[i].PoreDist = (float *) calloc((*SType)[i].NLayers,
						  sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*SType)[i].Press = (float *) calloc((*SType)[i].NLayers, 
					       sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*SType)[i].FCap = (float *) calloc((*SType)[i].NLayers,
					      sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*SType)[i].WP = (float *) calloc((*SType)[i].NLayers,
					    sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*SType)[i].Dens = (float *) calloc((*SType)[i].NLayers,
					      sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*SType)[i].Ks = (float *) calloc((*SType)[i].NLayers,
					    sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*SType)[i].KhDry = (float *) calloc((*SType)[i].NLayers,
					       sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*SType)[i].KhSol = (float *) calloc((*SType)[i].NLayers,
					       sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*SType)[i].Ch = (float *) calloc((*SType)[i].NLayers,
					    sizeof(float))))
      ReportError((char *) Routine, 1);
  
    if (!CopyFloat((*SType)[i].Porosity, VarStr[porosity], 
		    (*SType)[i].NLayers))
      ReportError(KeyName[porosity], 51);

    if (!CopyFloat((*SType)[i].PoreDist, VarStr[pore_size], 
		    (*SType)[i].NLayers))
      ReportError(KeyName[pore_size], 51);

    if (!CopyFloat((*SType)[i].Press, VarStr[bubbling_pressure], 
		    (*SType)[i].NLayers))
      ReportError(KeyName[bubbling_pressure], 51);

    if (!CopyFloat((*SType)[i].FCap, VarStr[field_capacity], 
		    (*SType)[i].NLayers))
      ReportError(KeyName[field_capacity], 51);

    if (!CopyFloat((*SType)[i].WP, VarStr[wilting_point], 
		    (*SType)[i].NLayers))
      ReportError(KeyName[wilting_point], 51);

    if (!CopyFloat((*SType)[i].Dens, VarStr[bulk_density], 
		    (*SType)[i].NLayers))
      ReportError(KeyName[bulk_density], 51);

    if (!CopyFloat((*SType)[i].Ks, VarStr[vertical_ks], 
		    (*SType)[i].NLayers))
      ReportError(KeyName[vertical_ks], 51);

    if (!CopyFloat((*SType)[i].KhSol, VarStr[solids_thermal], 
		    (*SType)[i].NLayers))
      ReportError(KeyName[solids_thermal], 51);

    if (!CopyFloat((*SType)[i].Ch, VarStr[thermal_capacity], 
		    (*SType)[i].NLayers))
      ReportError(KeyName[thermal_capacity], 51);
  }

  for (i = 0; i < NSoils; i++)
    for (j = 0; j < (*SType)[i].NLayers; j++){
      (*SType)[i].KhDry[j] = CalcKhDry((*SType)[i].Dens[j]);
    if (((*SType)[i].Porosity[j] < (*SType)[i].FCap[j]) || ((*SType)[i].Porosity[j] < (*SType)[i].WP[j]) ||
       ((*SType)[i].FCap[j] < (*SType)[i].WP[j]))
       ReportError((*SType)[i].Desc, 11);
    }

  return NSoils;
}

/********************************************************************************
  Function Name: InitVegTable()

  Purpose      : Initialize the vegetation lookup table
                 Processes most of the following section in the input file:
		 [VEGETATION]

  Required     :
    VEGTABLE **VType - Pointer to lookup table
    LISTPTR Input    - Pointer to linked list with input info
    LAYER *Veg       - Pointer to structure with veg layer information

  Returns      : Number of vegetation types

  Modifies     : VegTable and Veg

  Comments     :
********************************************************************************/
int InitVegTable(VEGTABLE **VType, LISTPTR Input, LAYER *Veg)
{
  const char *Routine = "InitVegTable";
  int i;			/* Counter */
  int j;			/* Counter */
  int NVegs;			/* Number of vegetation types */
  char KeyName[understory_fraction+1][BUFSIZE+1];
  char *KeyStr[] = {
    "VEGETATION DESCRIPTION",
    "OVERSTORY PRESENT",
    "UNDERSTORY PRESENT",
    "FRACTIONAL COVERAGE",
    "TRUNK SPACE",
    "AERODYNAMIC ATTENUATION",
    "RADIATION ATTENUATION",
    "MAX SNOW INT CAPACITY",
    "SNOW INTERCEPTION EFF",
    "MASS RELEASE DRIP RATIO",
    "HEIGHT",
    "OVERSTORY LAI",
    "UNDERSTORY LAI",
    "MAXIMUM RESISTANCE",
    "MINIMUM RESISTANCE",
    "MOISTURE THRESHOLD",
    "VAPOR PRESSURE DEFICIT",
    "RPC",
    "ALBEDO",
    "NUMBER OF ROOT ZONES",
    "ROOT ZONE DEPTHS",
    "OVERSTORY ROOT FRACTION",
    "UNDERSTORY ROOT FRACTION"
  };
  char SectionName[] = "VEGETATION";
  char VarStr[understory_fraction+1][BUFSIZE+1];
  

  /* Get the number of different vegetation types */
  GetInitString(SectionName, "NUMBER OF VEGETATION TYPES", "", VarStr[0], 
		(unsigned long) BUFSIZE, Input);
  if (!CopyInt(&NVegs, VarStr[0], 1))
    ReportError("NUMBER OF VEGETATION TYPES", 51);

  if (NVegs == 0)
    return NVegs;
  
  if (!(Veg->NLayers = (int *) calloc(NVegs, sizeof(int))))
    ReportError((char *) Routine, 1);
  
  if (!(*VType = (VEGTABLE *) calloc(NVegs, sizeof(VEGTABLE))))
    ReportError((char *) Routine, 1);
  
  /******* Read information and allocate memory for each vegetation type ******/

  Veg->MaxLayers = 0;

  printf("\n\n Summary information about the overstory vegetation classes\n");
  printf("For each category value is percent transfer of quantity from above to below canopy\n\n");
  printf("        Descripition\tS Rad Summer\tS Rad Winter\tWind\tMax Snow Int (mm)\n");

  for (i = 0; i < NVegs; i++) {

    /* Read the key-entry pairs from the input file */
    for (j = 0; j <= understory_fraction; j++) {
      sprintf(KeyName[j], "%s %d", KeyStr[j], i+1);
      GetInitString(SectionName, KeyName[j], "", VarStr[j], 
		    (unsigned long) BUFSIZE, Input);
    }
    
    /* Assign the entries to the appropriate variables */
    if (IsEmptyStr(VarStr[veg_description]))
      ReportError(KeyName[veg_description], 51);
    strcpy((*VType)[i].Desc, VarStr[veg_description]);
    (*VType)[i].Index = i;
    
    (*VType)[i].NVegLayers = 0;
    
    if (strncmp(VarStr[overstory], "TRUE", 4) == 0) {
      (*VType)[i].OverStory = TRUE;
      ((*VType)[i].NVegLayers)++;
    }
    else if (strncmp(VarStr[overstory], "FALSE", 5) == 0)
      (*VType)[i].OverStory = FALSE;
    else
      ReportError(KeyName[overstory], 51);
    
    if (strncmp(VarStr[understory], "TRUE", 4) == 0) {
      (*VType)[i].UnderStory = TRUE;
      ((*VType)[i].NVegLayers)++;
    }
    else if (strncmp(VarStr[understory], "FALSE", 5) == 0)
      (*VType)[i].UnderStory = FALSE;
    else
      ReportError(KeyName[understory], 51);
    
    Veg->NLayers[i] = (*VType)[i].NVegLayers;
    if ((*VType)[i].NVegLayers > Veg->MaxLayers)
      Veg->MaxLayers = (*VType)[i].NVegLayers;

    if (!CopyInt(&(*VType)[i].NSoilLayers, VarStr[number_of_root_zones], 1))
      ReportError(KeyName[number_of_root_zones], 51);
    
    /* allocate memory for the vegetation layers */
    if (!((*VType)[i].Fract = (float *) calloc((*VType)[i].NVegLayers,
					       sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*VType)[i].Height = (float *) calloc((*VType)[i].NVegLayers,
						sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*VType)[i].RsMax = (float *) calloc((*VType)[i].NVegLayers,
					       sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*VType)[i].RsMin = (float *) calloc((*VType)[i].NVegLayers,
					       sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*VType)[i].MoistThres = (float *) calloc((*VType)[i].NVegLayers,
						    sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*VType)[i].VpdThres = (float *) calloc((*VType)[i].NVegLayers,
						  sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*VType)[i].Rpc = (float *) calloc((*VType)[i].NVegLayers,
					     sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*VType)[i].Albedo = (float *) calloc(((*VType)[i].NVegLayers+1),
						sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*VType)[i].MaxInt = (float *) calloc((*VType)[i].NVegLayers,
						sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*VType)[i].LAI = (float *) calloc((*VType)[i].NVegLayers,
					     sizeof(float))))
      ReportError((char *) Routine, 1);
    if (!((*VType)[i].RootFract = (float **) calloc((*VType)[i].NVegLayers,
						    sizeof(float *))))
      ReportError((char *) Routine, 1);    
    for (j = 0; j < (*VType)[i].NVegLayers; j++) {
      if (!((*VType)[i].RootFract[j] = 
	    (float *) calloc((*VType)[i].NSoilLayers, sizeof(float))))
	ReportError((char *) Routine, 1);
    }      

    if (!((*VType)[i].MonthLAI = (float **) calloc((*VType)[i].NVegLayers,
						    sizeof(float *))))
      ReportError((char *) Routine, 1);    
    for (j = 0; j < (*VType)[i].NVegLayers; j++) {
      if (!((*VType)[i].MonthLAI[j] = 
	    (float *) calloc(12, sizeof(float))))
	ReportError((char *) Routine, 1);
    }      

    
    if (!((*VType)[i].RootDepth = (float *) calloc((*VType)[i].NSoilLayers,
						   sizeof(float))))
      ReportError((char *) Routine, 1);

    /* assign the entries to the appropriate variables */
    if ((*VType)[i].OverStory == TRUE) {

      if (!CopyFloat(&((*VType)[i].Fract[0]), VarStr[fraction], 1))
	ReportError(KeyName[fraction], 51);

      if (!CopyFloat(&((*VType)[i].Trunk), VarStr[trunk_space], 1))
	ReportError(KeyName[trunk_space], 51);

      if (!CopyFloat(&((*VType)[i].Cn), VarStr[aerodynamic_att], 1))
	ReportError(KeyName[aerodynamic_att], 51);

      if (!CopyFloat(&((*VType)[i].Atten), VarStr[radiation_att], 1))
	ReportError(KeyName[radiation_att], 51);

      if (!CopyFloat(&((*VType)[i].MaxSnowInt), VarStr[snow_int_cap], 1))
	ReportError(KeyName[snow_int_cap], 51);

      if (!CopyFloat(&((*VType)[i].MDRatio), VarStr[mass_drip_ratio], 1))
	ReportError(KeyName[mass_drip_ratio], 51);

      if (!CopyFloat(&((*VType)[i].SnowIntEff), VarStr[snow_int_eff], 1))
    	ReportError(KeyName[snow_int_eff], 51);

      if (!CopyFloat((*VType)[i].RootFract[0], VarStr[overstory_fraction], 
		      (*VType)[i].NSoilLayers))
    	ReportError(KeyName[overstory_fraction], 51);

      if (!CopyFloat((*VType)[i].MonthLAI[0], VarStr[overstory_lai], 12) ) 
    	ReportError(KeyName[overstory_lai], 51);

      if ((*VType)[i].UnderStory == TRUE) {
	(*VType)[i].Fract[1] = 1.0;
	if (!CopyFloat((*VType)[i].RootFract[1], VarStr[understory_fraction],
			(*VType)[i].NSoilLayers))
	  ReportError(KeyName[understory_fraction], 51);
      }
    }
    else {
      if ((*VType)[i].UnderStory == TRUE) {
	(*VType)[i].Fract[0] = 1.0;
      	if (!CopyFloat((*VType)[i].RootFract[0], VarStr[understory_fraction],
			(*VType)[i].NSoilLayers))
	  ReportError(KeyName[understory_fraction], 51);
      }
      (*VType)[i].Trunk = NOT_APPLICABLE;
      (*VType)[i].Cn    = NOT_APPLICABLE;
      (*VType)[i].Atten = NOT_APPLICABLE;
    }
    
    if ((*VType)[i].UnderStory == TRUE) {
      if (!CopyFloat((*VType)[i].MonthLAI[(*VType)[i].NVegLayers -1], 
        VarStr[understory_lai], 12) ) 
    	ReportError(KeyName[understory_lai], 51);
    }

    if (!CopyFloat((*VType)[i].Height, VarStr[height], 
		    (*VType)[i].NVegLayers))
      ReportError(KeyName[height], 51);

    if (!CopyFloat((*VType)[i].RsMax, VarStr[max_resistance],
		    (*VType)[i].NVegLayers))
      ReportError(KeyName[max_resistance], 51);

    if (!CopyFloat((*VType)[i].RsMin, VarStr[min_resistance], 
		    (*VType)[i].NVegLayers))
      ReportError(KeyName[min_resistance], 51);

    if (!CopyFloat((*VType)[i].MoistThres, VarStr[moisture_threshold], 
		    (*VType)[i].NVegLayers))
      ReportError(KeyName[moisture_threshold], 51);

    if (!CopyFloat((*VType)[i].VpdThres, VarStr[vpd], (*VType)[i].NVegLayers))
      ReportError(KeyName[vpd], 51);

    if (!CopyFloat((*VType)[i].Rpc, VarStr[rpc], (*VType)[i].NVegLayers))
      ReportError(KeyName[rpc], 51);

    if (!CopyFloat((*VType)[i].Albedo, VarStr[veg_albedo], 
		    (*VType)[i].NVegLayers))
      ReportError(KeyName[veg_albedo], 51);

    if (!CopyFloat((*VType)[i].RootDepth, VarStr[root_zone_depth], 
		    (*VType)[i].NSoilLayers))
      ReportError(KeyName[root_zone_depth], 51);
    
    /* Calculate the wind speed profiles and the aerodynamical resistances
       for each layer.  The values are normalized for a reference height wind
       speed of 1 m/s, and are adjusted each timestep using actual reference 
       height wind speeds */
    
    CalcAerodynamic((*VType)[i].NVegLayers, (*VType)[i].OverStory, 
		    (*VType)[i].Cn, (*VType)[i].Height, (*VType)[i].Trunk, 
		    (*VType)[i].U, &((*VType)[i].USnow), (*VType)[i].Ra, 
		    &((*VType)[i].RaSnow));
    /* Do some basic quality control on the input data: */
    /* Provide errors if outside of range */
    /* k*LAI > 3.0 will yield radiation transmission less than 5% */

    

    /*
    if ( (*VType)[i].OverStory == TRUE ) {
    printf("%20s\t%9.2f\t%9.2f\t%9.2f%9.2f\n",(*VType)[i].Desc,
	   100.0*exp(-((*VType)[i].Atten * (*VType)[i].LAISummer[0])),
           100.0*exp(-((*VType)[i].Atten * (*VType)[i].LAIWinter[0])),
	   (*VType)[i].USnow*100.0,(*VType)[i].MaxSnowInt*1000.0 );
    /*    printf("Ra Overstory = %f Ra UnderStory = %f Ra Snow = %f \n",
	  (*VType)[i].Ra[0],(*VType)[i].Ra[1],(*VType)[i].RaSnow);
    if( (*VType)[i].Atten * (*VType)[i].LAISummer[0] > 3.0)
      ReportError((*VType)[i].Desc, 12);
    if( (*VType)[i].Atten * (*VType)[i].LAIWinter[0] > 3.0)
      ReportError((*VType)[i].Desc, 12);
    if( (*VType)[i].USnow < 0.05)
      ReportError((*VType)[i].Desc, 13);
    }
    */
      


  }
  printf("\n");
  return NVegs;
}

/********************************************************************************
  InitSnowTable()

  Source:
  Laramie, R. L., and J. C. Schaake, Jr., Simulation of the continuous
      snowmelt process, Ralph M. Parsons Laboratory, Mass. Inst. of Technol.,
      1972

  Snow albedo is calculated as a function of the number of days since the
  last observed snow fall.  There are separete albedo curves for the freeze
  and thaw conditions.
********************************************************************************/
void InitSnowTable(SNOWTABLE **SnowAlbedo, int StepsPerDay)
{
  const char *Routine = "InitSnowTable";
  int i;
  
  if (!(*SnowAlbedo = (SNOWTABLE *) calloc((int)((DAYPYEAR+1)*StepsPerDay), 
                                           sizeof(SNOWTABLE))))
    ReportError((char *) Routine, 1);
  
  /* Laramie and Schaake (1972) */
  for (i = 0; i < ((DAYPYEAR+1)*StepsPerDay); i++) {
    (*SnowAlbedo)[i].Freeze = 
      0.85 * pow(0.92, pow(((float)i)/StepsPerDay, 0.58)); 
    if ((*SnowAlbedo)[i].Freeze<0.4) (*SnowAlbedo)[i].Freeze=0.4;
    (*SnowAlbedo)[i].Thaw   = 
      0.85 * pow(0.87, pow(((float)i)/StepsPerDay, 0.46));
    if ((*SnowAlbedo)[i].Thaw<0.4) (*SnowAlbedo)[i].Thaw=0.4;
  }
}
