/*
 * SUMMARY:      InitNetwork.c - Initialize road/channel work
 * USAGE:        
 *
 * AUTHOR:       DHSVM Project (Bart Nijssen)
 * ORG:          University of Washington, Department of Civil Engineering
 * E-MAIL:       nijssen@u.washington.edu
 * ORIG-DATE:    27-Aug-1996 at 18:34:01
 * LAST-MOD: Tue Dec  3 09:02:12 2002 by Scott Waichler <waichler@tuff.pnl.gov>
 * DESCRIPTION:  Initialize road/channel work.  Memory is allocated, and the
 *               necessary adjustments for the soil profile are calculated
 * DESCRIP-END.
 * FUNCTIONS:    InitNetwork()
 * COMMENTS:     
 */

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

#include <stdio.h>
#include <stdlib.h>
#include "constants.h"
#include "settings.h"
#include "data.h"
#include "DHSVMerror.h"
#include "functions.h"
#include "soilmoisture.h"
#include "DHSVMChannel.h"

/*****************************************************************************
  Function name: InitNetwork()

  Purpose      : Initialize road/channel work.  Memory is allocated, and the
                 necessary adjustments for the soil profile are calculated

  Required     :
    int HasChannels       - TRUE if Network, FALSE if NO_Network
    int NY                - Number of rows (North-South) in area
    int NX                - Number of columns (East-West) in area
    float DX              - Pixel width (m)
    SOILPIX **SoilMap     - Structure with information on soil type for each
                            pixel 
    VEGPIX **VegMap       - Structure with information on vegetation type for
                            each pixel, including number of rooting layers
    VEGTABLE *VType       - Lookup table with information in vegetation type,
                            including rooting depths (soil layer thickness)
                            for each vegetation type 

  Returns      : void

  Modifies     :
    ROADSTRUCT ***Network - Address of two dimensional array to contain
                            network information

  Comments     : It is assumed that we are dealing with square pixels,
                 i.e. DX is equal to DY
*****************************************************************************/
void InitNetwork(int HasNetwork, int NY, int NX, float DX, TOPOPIX **TopoMap,
                 SOILPIX **SoilMap, VEGPIX **VegMap, VEGTABLE *VType,
                 ROADSTRUCT ***Network, CHANNEL *ChannelData)
{
  const char *Routine = "InitNetwork";
  int i;                        /* counter */
  int x;                        /* column counter */
  int y;                        /* row counter */
  
  /* Allocate memory for network structure */
  
  if (!(*Network = (ROADSTRUCT **) calloc(NY, sizeof(ROADSTRUCT *))))
    ReportError((char *) Routine, 1);
  
  for (y = 0; y < NY; y++) {
    if (!((*Network)[y] = (ROADSTRUCT *) calloc(NX, sizeof(ROADSTRUCT))))
      ReportError((char *) Routine, 1);
  }
  
  for (y = 0; y < NY; y++) {
    for (x = 0; x < NX; x++) {
      if (InBasin(TopoMap[y][x].Mask)) {
        if (!((*Network)[y][x].Adjust = 
              (float *) calloc((VType[VegMap[y][x].Veg-1].NSoilLayers+1),
                               sizeof(float))))
          ReportError((char *) Routine, 1);
        
        if (!((*Network)[y][x].PercWidth = 
              (float *) calloc((VType[VegMap[y][x].Veg-1].NSoilLayers+1),
                               sizeof(float))))
          ReportError((char *) Routine, 1);
      }
    }
  }
  
  /* If a road/channel Network is imposed on the area, read the Network
     information, and calculate the storage adjustment factors */
  if (HasNetwork) {
    for (y = 0; y < NY; y++) {
      for (x = 0; x < NX; x++) {
        if (InBasin(TopoMap[y][x].Mask)) {

/* printf("bh4=%.3f ", ChannelData->stream_class->bank_height); */ /* this is ok */

          ChannelCut(y, x, ChannelData, &((*Network)[y][x]));

	      AdjustStorage(VType[VegMap[y][x].Veg-1].NSoilLayers,
                        SoilMap[y][x].Depth,
                        VType[VegMap[y][x].Veg-1].RootDepth,
                        (*Network)[y][x].Width, DX,
                        (*Network)[y][x].BankHeight,
                        (*Network)[y][x].PercWidth,
                        (*Network)[y][x].Adjust,
                        &((*Network)[y][x].CutBankZone));
          if (channel_grid_has_channel(ChannelData->road_map, x, y)) {
            (*Network)[y][x].fraction = 
              ChannelFraction(&(TopoMap[y][x]), ChannelData->road_map[x][y]);
          }
        }
      }
    }
  }
  /* if no road/channel Network is imposed, set the adjustment factors to the
     values they have in the absence of an imposed network */
  else {
    for (y = 0; y < NY; y++) {
      for (x = 0; x < NX; x++) {
        if (InBasin(TopoMap[y][x].Mask)) {
          for (i = 0; i <= VType[VegMap[y][x].Veg-1].NSoilLayers; i++) {
            (*Network)[y][x].Adjust[i] = 1.0;
            (*Network)[y][x].PercWidth[i] = 1.0;
            (*Network)[y][x].CutBankZone = NO_CUT;
          }
        }
      }
    }       
  }
}
