/*
 * SUMMARY:      RouteSurface.c - Route surface flow
 * 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: Wed Oct 30 11:00:49 2002 by Scott Waichler <waichler@tuff.pnl.gov>
 * DESCRIPTION:  Route surface flow
 * DESCRIP-END.
 * FUNCTIONS:    RouteSurface()
 * COMMENTS:     
 */

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

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include "settings.h"
#include "data.h"
#include "slopeaspect.h"
#include "DHSVMerror.h"
#include "functions.h"
#include "constants.h"

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


/*****************************************************************************
  RouteSurface()

  If the watertable calculated in WaterTableDepth() was negative, then water
  is ponding on the surface.  At the moment no permanent ponding of water is allowed 
  in DHSVM, and all of the "excess" water is routed to the outlet or stream channel.
  
  There are two basic considerations for moving surface water:  (I) How much water
  in a cell is to be moved off of that cell; (II) How frequently does water move from
  one cell to the next.  For (I), the user may choose either 1) all of the water in the
  cell; or 2) some of the water, according to a power law function.  For (II), the user
  may choose either a) water moves one cell per timestep; or b) water moves all the way
  down the hillslope, across all cells in a flowpath, as a cascade within the timestep.
  In summary:
  1a.  All water moved off the cell; moves one cell per timestep.
  1b.  All water moved off the cell; moves thru all cells per timestep.
  2a.  Some water moved off the cell; moves one cell per timestep.
  2b.  Some water moved off the cell; moves thru all cells per timestep.

  In the original DHSVM, water moved as 1a.  Subsequently, Mark added the capability
  to move water as 2a.  Comparing the two, 1a is the upper limit on water movement
  rate, while 2a slows it down.  For either option, 
  if the timestep or number of cells in a hillslope is too large, the user may 
  experience too-slow movement of water to the channels,
  and difficulty in reproducing peaks in observed hydrographs.

  The latest version here by Scott and Mark allows movement by
  1b or 2b.  In the case of 1b, all of the Runoff in the watershed will be moved down
  to the stream channel and be intercepted at the next timestep (linear cascade).
  In the case of 2b, some of the water on each cell will remain according to the power
  law function, and will ultimately infiltrate (power law cascade).  For both 1b and 2b,
  the movement of water in a cascade downslope to the receiving channel
  or basin pour point is ensured by processing the grid cells in descending order
  of elevation.  This results in moving water out of the highest cells first, then down through
  the lowest cells, all within the same timestep.  

  If a grid cell contains a channel, no movement out of that 
  cell is allowed, with the result that the water (Runoff) will 
  remain on the cell until RouteChannel() is called the next timestep.  A perhaps better
  solution would be to have the water enter the stream channel, and have the
  channel routing process all of the water intercepted from the subsurface and surface
  in the same timestep.  This would require moving the call to RouteChannel()
  after RouteSurface() in main().  It is not there now for some good reason, which is
  now forgotten.  It probably has to do with avoiding flow loops when road channels
  are involved.   
  
  Option 1b will probably rarely be used, as it results in 
  little difference in the output from the input (they would be exactly
  the same if channel interception took place after surface routing in the main program).
  Option 2b is a true kinematic wave with a delayed rise, peak,
  and gradual recession.  Numerical stability of the wave should be checked in practice.
  Remember that normally the kinematic wave is solved on the order of seconds or minutes,
  so some unusual values of the power law coefficients might be needed to slow down
  flow across the hillslope to realistic rates for a given dx, dt, and hillslope length.  
  Both 1b and 2b are a big departure from the original DHSVM.  Be careful! 

  Ultimately, none of these solutions seem satisfactory for matching the conceptual
  model of macropore flow and translatory flow dependent on moisture content.  
    

*****************************************************************************/
void RouteSurface(MAPSIZE *Map, TIMESTRUCT *Time, TOPOPIX **TopoMap, 
                  SOILPIX **SoilMap, int HasNetwork, 
                  UNITHYDR **UnitHydrograph, 
                  UNITHYDRINFO *HydrographInfo, float *Hydrograph, 
                  FILES *StreamFile, int HydraulicRouting, SOILTABLE *SType, 
                  CHANNEL *ChannelData )
{

  const char *Routine = "RouteSurface";
  int Lag;                      /* Lag time for hydrograph */
  int Step;
  float StreamFlow, testRunoff;
  int TravelTime;
  int WaveLength;

  int i, j, k;
  int x;
  int y;
  int n;
  float **surface;
  float total_in = 0.0, total_out = 0.0, inputVol = 0.0;
  float outflow;

  if (HasNetwork) {
    if ((surface = (float * *)malloc(Map->NY*sizeof(float *))) == NULL) {
      ReportError((char *) Routine, 1);
    }

    /* Decide to do either one cell-at-a-time routing (first block), or
       cascade routing (second block).  */

    if( HydraulicRouting != CASCADE_LINEAR && HydraulicRouting != CASCADE_POWERLAW) {  

      /* original DHSVM code:  one cell-at-a-time routing */
      for (y = 0; y < Map->NY; y++) {
        if ((surface[y] = (float *)malloc(Map->NX*sizeof(float))) == NULL) {
          ReportError((char *) Routine, 1);
        } else {
          for (x = 0; x < Map->NX; x++) {
            /*surface[y][x] = SoilMap[y][x].Runoff;*/
            if (InBasin(TopoMap[y][x].Mask)) {
              if( HydraulicRouting != LINEAR && SoilMap[y][x].Runoff > 0) {
                /* option 2a */
                inputVol = (SoilMap[y][x].Runoff - SoilMap[y][x].startRunoff);

                surface[y][x] = soilPowerLawOutflow( &(SoilMap[y][x].startRunoff), 
                                   SType[SoilMap[y][x].Soil-1].quickflowCoeff, 
                                   SType[SoilMap[y][x].Soil-1].quickflowExp, 
                                   inputVol, (float) Time->Dt );

                assert( surface[y][x] > 0 ); /* && surface[y][x] < testRunoff + inputVol);*/

                SoilMap[y][x].Runoff = SoilMap[y][x].startRunoff;
              } else {  
                /* constant velocity (v = dx/dt) routing */
                /* option 1a */
                surface[y][x] = SoilMap[y][x].Runoff;
                total_in += SoilMap[y][x].Runoff;
                SoilMap[y][x].Runoff = 0.0;
              }
            }
          }
        }
      }

      for (y = 0; y < Map->NY; y++) {
        for (x = 0; x < Map->NX; x++) {
          if (InBasin(TopoMap[y][x].Mask)) {
            for (n = 0; n < NDIRS; n++) {
              int xn = x + Map->xneighbor[n];
              int yn = y + Map->yneighbor[n];
              if (valid_cell(Map, xn, yn)) {
                SoilMap[yn][xn].Runoff += 
                  surface[y][x]*((float)TopoMap[y][x].Dir[n]/
                                 (float)TopoMap[y][x].TotalDir);
              }
            }
          }
        }
      }

      for (y = 0; y < Map->NY; y++) {
        for (x = 0; x < Map->NX; x++) {
          if (InBasin(TopoMap[y][x].Mask)) {
            total_out += SoilMap[y][x].Runoff;
          }
        }
      }

      if (DEBUG)
        printf("RouteSurface: in = %.4g, out = %.4g\n", total_in, total_out);

      for (y = 0; y < Map->NY; y++) {
        free(surface[y]);
      }
      free(surface);
  
    } else { /* use cascade routing--water moves across entire hillslope in timestep */

      /* Loop thru all of the cells in descending order of elevation */
      for (k = 0; k < Map->NumCells; k++) {
        y = Map->OrderedCellsRowCol[k][0];
        x = Map->OrderedCellsRowCol[k][1];
  
        outflow = 0.0;
  
        /* If a channel cell or no Runoff is available, skip this and move no water */
        if(!channel_grid_has_channel(ChannelData->stream_map, x, y) 
           && SoilMap[y][x].Runoff > 0) {
  
          if( HydraulicRouting == CASCADE_POWERLAW) {
            /* option 2b */
            inputVol = (SoilMap[y][x].Runoff - SoilMap[y][x].startRunoff);

            outflow = soilPowerLawOutflow( &(SoilMap[y][x].startRunoff), 
                                     SType[SoilMap[y][x].Soil-1].quickflowCoeff, 
                                     SType[SoilMap[y][x].Soil-1].quickflowExp, 
                                     inputVol, (float) Time->Dt );
            if(outflow < 0.0) outflow = 0.0;
            /* assert(outflow >= 0 );  *//* && surface[y][x] < testRunoff + inputVol);*/

            SoilMap[y][x].Runoff = SoilMap[y][x].startRunoff;

          } else if ( (HydraulicRouting == CASCADE_LINEAR) ||
                      (HydraulicRouting == CASCADE_LINEAR_LIMITED &&
                       SoilMap[y][x].startRunoff > 0) ) {
            /* option 1b */
            /* First condition: all of the water moves to the next cell */
            /* Second condition: all of the water moves to the next cell, but only
               because the cell was wet before starting RouteSurface(). */
            outflow = SoilMap[y][x].Runoff;
            SoilMap[y][x].Runoff = 0.0;
          }        
        } /*  end if not a channel cell, and Runoff > 0 */
  
        if(outflow > 0) {  
          for (n = 0; n < NDIRS; n++) {
            int xn = x + Map->xneighbor[n];
            int yn = y + Map->yneighbor[n];
            if (valid_cell(Map, xn, yn)) {
                SoilMap[yn][xn].Runoff += 
                    outflow*((float)TopoMap[y][x].Dir[n]/(float)TopoMap[y][x].TotalDir);
            } 
          } /* end loop thru possible flow directions */
        }

        total_in += SoilMap[y][x].Runoff;
        total_out += SoilMap[y][x].Runoff;

      } /* end loop thru ordered basin cells */

      if (DEBUG)
        printf("RouteSurface: in = %.4g, out = %.4g\n", total_in, total_out);

    } /* end else block to use cascade routing */
  
  /* MAKE SURE THIS WORKS WITH A TIMESTEP IN SECONDS */
  } else {  /* No network, so use unit hydrograph method */
                                   
    for (y = 0; y < Map->NY; y++) {
      for (x = 0; x < Map->NX; x++) {
        if (InBasin(TopoMap[y][x].Mask)) {
          TravelTime = (int) TopoMap[y][x].Travel;
          if (TravelTime != 0) {
            WaveLength = HydrographInfo->WaveLength[TravelTime-1];
            for (Step = 0; Step < WaveLength; Step++) {
              Lag = UnitHydrograph[TravelTime-1][Step].TimeStep;
              Hydrograph[Lag] += SoilMap[y][x].Runoff * 
                UnitHydrograph[TravelTime-1][Step].Fraction;
            }
            SoilMap[y][x].Runoff = 0.0;
          }
        }
      } 
    }
    
    StreamFlow = 0.0;
    for (i = 0; i < Time->Dt; i++)
      StreamFlow += (Hydrograph[i] * Map->DX * Map->DY)/Time->Dt;
    
    /* Advance Hydrograph */
    for (i = 0; i < Time->Dt; i++) {
      for (j = 0; j < HydrographInfo->TotalWaveLength-1; j++) {
        Hydrograph[j] = Hydrograph[j+1];
      }
    }
    
    /* Set the last elements of the hydrograph to zero */
    for (i = 0; i < Time->Dt; i++)
      Hydrograph[HydrographInfo->TotalWaveLength-(i+1)] = 0.0;

    PrintDate(&(Time->Current), StreamFile->FilePtr);
    fprintf(StreamFile->FilePtr, " %g\n", StreamFlow);
  } /* end else, no network */
}  /* end RouteSurface() */
