/* -------------------------------------------------------------
   file: DHSVMChannel.c
   ------------------------------------------------------------- */
/* -------------------------------------------------------------
   Battelle Memorial Institute
   Pacific Northwest Laboratory
   ------------------------------------------------------------- */
/* -------------------------------------------------------------
   Created August 30, 1996 by  William A Perkins
   Last Change: Tue Dec 10 09:12:33 2002 by Scott Waichler <waichler@tuff.pnl.gov>
   ------------------------------------------------------------- */

static const char* RCS_ID = "$Id: DHSVMChannel.c,v 1.4 2001/05/15 17:30:24 waichler Exp $\nBattelle PNNL";

#include <stdio.h>
#include <math.h>
#include <string.h>
#include "getinit.h"
#include "DHSVMChannel.h"
#include "DHSVMerror.h"
#include "functions.h"
#include "settings.h"
#include "errorhandler.h"
#include "fileio.h"
#include "elevationband.h"

extern 
   float soilPowerLawOutflow( float *store, float coef, float exp, float inputVol, 
                            float dt );

/* -----------------------------------------------------------------------------
   InitChannel
   Reads stream and road files and builds the networks.
   -------------------------------------------------------------------------- */
void
InitChannel(LISTPTR Input, MAPSIZE *Map, int deltat, CHANNEL *channel, 
            SOILPIX **SoilMap, TOPOPIX **TopoMap)
{
  int i,x, y;
  STRINIENTRY StrEnv[] = {
    {"ROUTING", "STREAM NETWORK FILE", "", ""},
    {"ROUTING", "STREAM MAP FILE"    , "", ""},
    {"ROUTING", "STREAM CLASS FILE"  , "", ""},
    {"ROUTING", "ROAD NETWORK FILE"  , "", "none"},
    {"ROUTING", "ROAD MAP FILE"      , "", "none"},
    {"ROUTING", "ROAD CLASS FILE"    , "", "none"},
    {NULL     , NULL                 , "", NULL}
  };

  printf("Initializing Road/Stream Networks\n");  

  /* Read the key-entry pairs from the ROUTING section in 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);
  }

  channel->stream_class = NULL;
  channel->road_class = NULL;
  channel->streams = NULL;
  channel->roads = NULL;
  channel->stream_map = NULL;
  channel->road_map = NULL;
  
  channel_init();
  channel_grid_init(Map->NX, Map->NY);
  
  if (strncmp(StrEnv[stream_class].VarStr, "none", 4)) {

    printf("\tReading Stream data\n");

    if ((channel->stream_class = 
         channel_read_classes(StrEnv[stream_class].VarStr)) == NULL) {
      ReportError(StrEnv[stream_class].VarStr, 5);
    }
    if ((channel->streams =
         channel_read_network(StrEnv[stream_network].VarStr, 
			      channel->stream_class)) == NULL) {
      ReportError(StrEnv[stream_network].VarStr, 5);
    }
    if ((channel->stream_map = 
         channel_grid_read_map(channel->streams, 
			       StrEnv[stream_map].VarStr, SoilMap)) == NULL) {
      ReportError(StrEnv[stream_map].VarStr, 5);
    }
    error_handler (ERRHDL_STATUS,
                   "InitChannel: computing stream network routing coefficients");
    channel_routing_parameters (channel->streams, (double)deltat);
  
    for (y = 0; y < Map->NY; y++) {
      for (x = 0; x < Map->NX; x++) {
        if (InBasin (TopoMap[y][x].Mask)) {
          if ( channel_grid_has_channel(channel->stream_map, x, y)) 
            channel_grid_count_cell(channel->stream_map, x, y);

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

        }
      }
    }


  }
  
  if (strncmp(StrEnv[road_class].VarStr, "none", 4)) {

    printf("\tReading Road data\n");

    if ((channel->road_class = 
         channel_read_classes(StrEnv[road_class].VarStr)) == NULL) {
      ReportError(StrEnv[road_class].VarStr, 5);
    }
    if ((channel->roads =
         channel_read_network(StrEnv[road_network].VarStr, 
			      channel->road_class)) == NULL) {
      ReportError(StrEnv[road_network].VarStr, 5);
    }
    if ((channel->road_map = 
         channel_grid_read_map(channel->roads, 
			       StrEnv[road_map].VarStr, SoilMap)) == NULL) {
      ReportError(StrEnv[road_map].VarStr, 5);
    }
    error_handler (ERRHDL_STATUS,
                   "InitChannel: computing road network routing coefficients");
    channel_routing_parameters (channel->roads, (double)deltat);
  }
}

/* -------------------------------------------------------------
   InitChannelDump
   ------------------------------------------------------------- */
void
InitChannelDump(CHANNEL *channel, char *DumpPath)
{
  char buffer[NAMESIZE];
  
  if (channel->streams != NULL) {
    sprintf(buffer, "%sStream.Flow", DumpPath);
    OpenFile(&(channel->streamout), buffer, "w", TRUE);
    sprintf(buffer, "%sStreamflow.Only", DumpPath);
    OpenFile(&(channel->streamflowout), buffer, "w", TRUE);
  }
  if (channel->roads != NULL) {
    sprintf(buffer, "%sRoad.Flow", DumpPath);
    OpenFile(&(channel->roadout), buffer, "w", TRUE);
    sprintf(buffer, "%sRoadflow.Only", DumpPath);
    OpenFile(&(channel->roadflowout), buffer, "w", TRUE);

  }
}

/* -------------------------------------------------------------
   ChannelCulvertFlow    
   computes outflow of channel/road network to a grid cell, if it
   contains a sink
   ------------------------------------------------------------- */
double
ChannelCulvertFlow(int y, int x, CHANNEL *ChannelData)
{
  if (channel_grid_has_channel(ChannelData->road_map, x, y)) { 
    return channel_grid_outflow(ChannelData->road_map, x, y);
  } else {
    return 0;
  }
}

/* -------------------------------------------------------------
   RouteChannel
   ------------------------------------------------------------- */
void
RouteChannel(CHANNEL *ChannelData, TIMESTRUCT *Time, MAPSIZE *Map, 
             TOPOPIX **TopoMap, SOILPIX **SoilMap, AGGREGATED *Total,
             OPTIONSTRUCT Options, SOILTABLE *SType  )
{
  int x, y;
  int flag;
  char buffer[32];
  float CulvertFlow, inputVol = 0, outputVol = 0;
  

                                /* give any surface water to roads w/o sinks */

  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) {
      if (InBasin(TopoMap[y][x].Mask)) {
        if (channel_grid_has_channel(ChannelData->road_map, x, y) &&  
           !channel_grid_has_sink(ChannelData->road_map, x, y) ) {  /* road w/o sink */
             SoilMap[y][x].RoadInt += SoilMap[y][x].Runoff;  
             channel_grid_inc_inflow(ChannelData->road_map, x, y, 
                                      SoilMap[y][x].Runoff * TopoMap[y][x].Area);
             SoilMap[y][x].Runoff = 0.0f;
        }
      }
    }
  }

                                /* route the road network and save
                                   results*/
  
  SPrintDate(&(Time->Current), buffer);
  flag=IsEqualTime(&(Time->Current),&(Time->Start));
  if (ChannelData->roads != NULL) {
    channel_route_network(ChannelData->roads, Time->Dt, Options.ChannelRouting);
    channel_save_outflow_text(buffer, ChannelData->roads, 
                              ChannelData->roadout,ChannelData->roadflowout,flag);
  }

                                /* add culvert outflow to surface water */
  Total->CulvertReturnFlow = 0.0; 
  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) {
      if (InBasin(TopoMap[y][x].Mask)) {
        CulvertFlow =  ChannelCulvertFlow(y, x, ChannelData);
        /* CulvertFlow /= Map->DX * Map->DY; */
        /* CulvertFlow = (CulvertFlow > 0.0) ? CulvertFlow : 0.0; */
        
        if (channel_grid_has_channel(ChannelData->stream_map, x, y)) {
          if( Options.Extent == ELEVBAND && Options.HydraulicRouting != LINEAR ){
            inputVol = (SoilMap[y][x].Runoff - SoilMap[y][x].startRunoff);
              
            outputVol = soilPowerLawOutflow( &(SoilMap[y][x].startRunoff), 
                                 SType[SoilMap[y][x].Soil-1].quickflowCoeff, 
                                 SType[SoilMap[y][x].Soil-1].quickflowExp, 
                                 inputVol, (float) Time->Dt );

            /* test --  
            SoilMap[y][x].startRunoff = 0;
            outputVol = SoilMap[y][x].Runoff;
            /* end test */

            /* keep track of quickflow seperately for geochemistry */
            channel_grid_inc_quickflow(ChannelData->stream_map, x, y, 
                                    outputVol * TopoMap[y][x].Area);

            channel_grid_inc_inflow(ChannelData->stream_map, x, y, 
                                  outputVol * TopoMap[y][x].Area + CulvertFlow);
            SoilMap[y][x].ChannelInt += outputVol;
            Total->CulvertToChannel += CulvertFlow / TopoMap[y][x].Area;
            Total->RunoffToChannel += outputVol;

            SoilMap[y][x].Runoff = SoilMap[y][x].startRunoff;
          } else {
            if( Options.Extent == ELEVBAND){
            /* keep track of quickflow seperately for geochemistry */
            channel_grid_inc_quickflow(ChannelData->stream_map, x, y, 
                                    SoilMap[y][x].Runoff * TopoMap[y][x].Area);
            }

            channel_grid_inc_inflow(ChannelData->stream_map, x, y, 
                                  SoilMap[y][x].Runoff * TopoMap[y][x].Area + CulvertFlow);
            SoilMap[y][x].ChannelInt += SoilMap[y][x].Runoff;
            Total->CulvertToChannel += CulvertFlow / TopoMap[y][x].Area;
            Total->RunoffToChannel += SoilMap[y][x].Runoff;
            SoilMap[y][x].Runoff = 0.0f;
          }
        }
        else {
          SoilMap[y][x].Runoff += CulvertFlow / TopoMap[y][x].Area;
          Total->CulvertReturnFlow += CulvertFlow / TopoMap[y][x].Area;
        } 
      } /* end if in basin */  
    } /* end x loop */
  } /* end y loop */
                               /* route stream channels */
  if (ChannelData->streams != NULL) {
    channel_route_network(ChannelData->streams, Time->Dt, Options.ChannelRouting);

    channel_save_outflow_text(buffer, ChannelData->streams, 
                              ChannelData->streamout, ChannelData->streamflowout,flag);
  }  


}

/* -------------------------------------------------------------
   ChannelCut
   computes necessary parameters for cell storage adjustment from
   channel/road dimensions
   ------------------------------------------------------------- */
void
ChannelCut(int y, int x, CHANNEL *ChannelData, ROADSTRUCT *Network)
{
  float bank_height = 0.0;
  float cut_width = 0.0;

  if (channel_grid_has_channel(ChannelData->stream_map, x, y)) {
    bank_height = channel_grid_cell_bankht(ChannelData->stream_map, x, y);

/* printf("bh6=%.3f ", bank_height); */ /* bad value */

    cut_width = channel_grid_cell_width(ChannelData->stream_map, x, y);
  } else if (channel_grid_has_channel(ChannelData->road_map, x, y)) {
    bank_height = channel_grid_cell_bankht(ChannelData->road_map, x, y);
    cut_width = channel_grid_cell_width(ChannelData->road_map, x, y);
  }
  Network->Width = cut_width;
  Network->BankHeight = bank_height;
}

/* -------------------------------------------------------------
   ChannelFraction
   This computes the (sub)surface flow fraction for a road
   ------------------------------------------------------------- */
uchar
ChannelFraction(TOPOPIX *topo, ChannelMapRec *rds)
{
  float effective_width = 0;
  float total_width;
  float sine, cosine;
  ChannelMapRec *r;
  float fract = 0.0;
  
  if (rds == NULL) {
    return 0;
  }
  cosine = cos(topo->Aspect);
  sine = sin(topo->Aspect);
  total_width = topo->FlowGrad/topo->Slope;
  effective_width = 0.0;

  for (r = rds; r != NULL; r = r->next) {
    effective_width += r->length*sin(fabs(topo->Aspect - r->aspect));
  }
  fract = effective_width/total_width*255.0;
  fract = (fract > 255.0 ? 255.0 : floor(fract + 0.5));
  
  return (uchar)fract;
}


    
