/* -------------------------------------------------------------
   file: channel.c
   ------------------------------------------------------------- */
/* -------------------------------------------------------------
   Battelle Memorial Institute
   Pacific Northwest Laboratory
   ------------------------------------------------------------- */
/* -------------------------------------------------------------
   Created October 24, 1995 by  William A Perkins
   Last Change: Tue Dec  3 09:02:25 2002 by Scott Waichler <waichler@tuff.pnl.gov>
   ------------------------------------------------------------- */


static const char* RCS_ID = "$Id: channel.c,v 1.6 2001/05/15 17:30:27 waichler Exp $";

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <errno.h>
#include <malloc.h>
#include "errorhandler.h"
#include "channel.h"
#include "tableio.h"
#include "boolean.h"
#include "settings.h"
/* #include "channels.h" */

/* for test msw */
#define TEST_MAIN 0
/* end test */


extern float route_segment_powerlaw(float inFlux, float dt, float startDepth, float width, float length, float slope, 
                    float roughness, float *storage);


extern float route_segment( int seg_id, float ivol, float dT, float *out_1 );

/* -------------------------------------------------------------
   -------------- ChannelClass Functions -----------------------
   ------------------------------------------------------------- */

/* -------------------------------------------------------------
   alloc_channel_class
   ------------------------------------------------------------- */
static ChannelClass *
alloc_channel_class(void) 
{
  ChannelClass *p;
  
  if ((p = (ChannelClass *)malloc(sizeof(ChannelClass))) == NULL) {
    error_handler(ERRHDL_ERROR, "alloc_channel_class: malloc failed: %s", 
                  strerror(errno));
    return NULL;
  }

  p->id = 0;
  p->width = 0.0;
  p->bank_height = 0.0;
  p->friction = 0.0;
  p->crown = CHAN_OUTSLOPED;
  p->next = (ChannelClass *)NULL;

  return p;
}

/* -------------------------------------------------------------
   channel_free_classes
   ------------------------------------------------------------- */
void
channel_free_classes(ChannelClass *head)
{
  if (head->next != NULL) {
    channel_free_classes(head->next);
  }
  free(head);
}


/* -------------------------------------------------------------
   find_channel_class
   ------------------------------------------------------------- */
static ChannelClass *
find_channel_class(ChannelClass *list, ClassID id) 
{
  while (list != (ChannelClass *)NULL) {
    if (list->id == id) break;
    list = list->next;
  }
  return list;
}


/* -------------------------------------------------------------
   channel_read_classes
   This function opens and reads the specified file and returns a
   linked list of ChannelClass structs.  If anything goes wrong, NULL
   is returned and any ChannelClass structs are destroyed.
   ------------------------------------------------------------- */
ChannelClass *
channel_read_classes(const char *file)
{
  ChannelClass *head = NULL, *current = NULL;
  static const int fields = 6;
  int done;
  int err = 0;
  static char *crown_words[4] = { 
    "OUTSLOPED", "CROWNED", "INSLOPED", NULL 
  };
  
  static TableField class_fields[6] = {
    { "ID", TABLE_INTEGER, TRUE, FALSE, {0}, "", NULL },
    { "Channel Width", TABLE_REAL, TRUE, FALSE, {0.0}, "", NULL },
    { "Bank (stream) or Cut Height (road)", TABLE_REAL, TRUE, FALSE, {0.0}, "",
      NULL },  
    { "Friction Coefficient (Manning's n)", TABLE_REAL, TRUE, FALSE, {0.0}, "",
      NULL }, 
    { "Road Crown Type", TABLE_WORD, FALSE, FALSE, {0}, "", crown_words }
  };
  
  error_handler(ERRHDL_STATUS, 
                "channel_read_classes: reading file \"%s\"",
                file);

  if (table_open(file) != 0) {
    error_handler(ERRHDL_ERROR, 
                  "channel_read_classes: unable to open file \"%s\": %s",
                  file, strerror(errno));
    return NULL;
  }

  done = FALSE;
  while (!done) {
    int i;
    
    done = (table_get_fields(fields, class_fields) < 0);
    if (done) {
      for (i = 0; i < fields; i++) {
        if (class_fields[i].read) break;
      }
      if (i >= fields) continue;
    }

    if (head == NULL) {
      head = alloc_channel_class();
      current = head;
    } else {
      current->next = alloc_channel_class();
      current = current->next;
    }
    
    for (i = 0; i < fields; i++) {
      if (class_fields[i].read) {
        switch (i) {
        case 0:
          current->id = class_fields[i].value.integer;
          if (current->id <= 0) {
            error_handler(ERRHDL_ERROR,
                          "%s: class %d: class id invalid",
                          file, current->id);
            err++;
          }
          break;
        case 1:
          current->width = class_fields[i].value.real;
          break;
        case 2:
          current->bank_height = class_fields[i].value.real;

/* printf("bh1=%.3f ", current->bank_height); */ /* this comes out ok */

          break;
        case 3:
          current->friction = class_fields[i].value.real;
          break;
        case 4:
          switch (class_fields[i].value.integer) {
          case -1:
            error_handler(ERRHDL_ERROR, 
                          "channel_read_classes: %s: unknown road crown type: %s",
                          file, class_fields[i].field);
            err++;
            break;
          case 0:
            current->crown = CHAN_OUTSLOPED;
            break;
          case 1:
            current->crown = CHAN_CROWNED;
            break;
          case 2:
            current->crown = CHAN_INSLOPED;
            break;
          default:
            error_handler(ERRHDL_FATAL, 
                          "channel_read_classes: this should not happen");
          }
          break;
        default:
          error_handler(ERRHDL_FATAL, 
                        "channel_read_classes: this should not happen either");
        }
      }
    }
  }
  error_handler(ERRHDL_STATUS, 
                "channel_read_classes: %s: %d errors, %d warnings",
                file, table_errors, table_warnings);

  table_close();

  error_handler(ERRHDL_STATUS, 
                "channel_read_classes: done reading file \"%s\"",
                file);
  
  if (table_errors) {
    error_handler(ERRHDL_ERROR,
                  "channel_read_classes: %s: too many errors",
                  file);
    channel_free_classes(head);
    head = NULL;
  }

  return(head);
}

/* -------------------------------------------------------------
   ---------------------- Channel Functions --------------------
   ------------------------------------------------------------- */

/* -------------------------------------------------------------
   alloc_channel_segment
   ------------------------------------------------------------- */
static Channel *
alloc_channel_segment(void) 
{
  Channel *seg;
  if ((seg = (Channel *)malloc(sizeof(Channel))) == NULL) {
    error_handler(ERRHDL_ERROR, "alloc_channel_segment: malloc failed: %s",
                  strerror(errno));
    return NULL;
  }
  seg->id = 0;
  seg->order = 0;
  seg->record_name = NULL;
  seg->record = FALSE;
  seg->length = 0.0;
  seg->slope = 0.0;
  seg->channel_class = NULL;
  seg->lateral_inflow = 0.0;
  seg->quick_flow = 0.0;
  seg->last_inflow = 0.0;
  seg->last_outflow = 0.0;
  seg->last_storage = 0;
  seg->storage = 0;
  seg->inflow = 0.0;
  seg->outflow = 0.0;
  seg->last_lateral_inflow = 0.0;
  seg->outlet = NULL;
  seg->numCells = 0;
  seg->recharge = 0.0;
  seg->next = NULL;
  
  return seg;
}

/* -------------------------------------------------------------
   channel_find_segment
   A simple linear search of the channel network to find a segment
   with the given id
   ------------------------------------------------------------- */
Channel *
channel_find_segment(Channel *head, SegmentID id)
{
  for (; head != NULL; head = head->next) {
    if (head->id == id) break;
  }
  if (head == NULL) {
    error_handler(ERRHDL_WARNING, 
                  "channel_find_segment: unable to find segment %d", 
                  id); 
  } else {
    error_handler(ERRHDL_DEBUG,
                  "channel_find_segment: found segment %d", 
                  id); 
  }
  
  return head;
}

/* -------------------------------------------------------------
   channel_routing_parameters
   ------------------------------------------------------------- */
void
channel_routing_parameters(Channel *network, int deltat)
{
  /*   float ck; */
  float y;
  Channel *segment;

  for (segment = network; segment != NULL; segment = segment->next) {

    y = segment->channel_class->bank_height*0.75;
  
                                /* must keep 0 <= X <= 0.5 */  

    /* Next line commented out by Bart Nijssen Wed Feb  3 18:57:09 1999 */
    /*     segment->X = 0.5*(1 - 0.6*y/(segment->slope*segment->length)); */

                                /* must keep 0 <= X <= 0.5 */  
    /* below portion commented out by westrick, shallow slopes and short segs */
    /* will yield very small y values */
    /*    if( segment->X < 0.0  ) { */  	    
    /*  segment->X = 0.;*/
    /*  y = (5.0 / 3.0 )*segment->slope*segment->length;*/	
    /* y which forces X = 0 */    
    /*  error_handler(ERRHDL_WARNING, */
    /*                "channel_routing_parameters: segment %d: X = %f", */
    /*                segment->id, segment->X); */
    /* } */

                                /* negative outflow can result unless 
                                   K < DT/2X so don't let K exceed that value */ 

    /* Next part commented out by Bart Nijssen Wed Feb  3 18:56:02 1999, since
       it is currently not being used.  Note that in the next few lines deltat
       is still in hours, NOT in seconds */
    /*     ck = 6000*sqrt(segment->slope)*pow(y,
	   2.0/3.0)/segment->class->friction; */ 
    /*     segment->K = segment->length/ck;	 */
    
    /*     if(segment->K > (deltat/(2.0*segment->X) ) ) { */
    /*       segment->K = deltat / (2.0*segment->X); */
    /*       error_handler(ERRHDL_WARNING, */
    /*                     "channel_routing_parameters: segment %d: K = %f", */
    /*                     segment->id, segment->K); */
    /*     } */

  /*  for new routing scheme */ 
   segment->K = sqrt(segment->slope)*pow((double) y, (double) (2.0/3.0)) / 
     (segment->channel_class->friction*segment->length); 
   segment->X = exp(-segment->K*deltat);
  }
  
  return;
}
  

/* -------------------------------------------------------------
   channel_read_network
   ------------------------------------------------------------- */
Channel *
channel_read_network(const char *file, ChannelClass *class_list)
{
  Channel *head = NULL, *current = NULL;
  int err = 0;
  int done;
  static const int fields = 8;
  static char *save_words[2] = {
    "SAVE", "\0"
  };
  static TableField chan_fields[8] = {
    { "ID", TABLE_INTEGER, TRUE, FALSE, {0}, "", NULL },
    { "Order", TABLE_INTEGER, TRUE, FALSE, {0}, "", NULL },
    { "Slope", TABLE_REAL, TRUE, FALSE, {0.0}, "", NULL },
    { "Length", TABLE_REAL, TRUE, FALSE, {0.0}, "", NULL }, 
    { "Class", TABLE_INTEGER, TRUE, FALSE, {0}, "", NULL },
    { "Outlet ID", TABLE_INTEGER, FALSE, FALSE, {0}, "", NULL }, 
    { "Save Flag", TABLE_WORD, FALSE, FALSE, {0}, "", save_words },
    { "Save Name", TABLE_STRING, FALSE, FALSE, {0.0}, "", NULL }
  };
  
  error_handler(ERRHDL_STATUS, 
                "channel_read_network: reading file \"%s\"",
                file);

  if (table_open(file) != 0) {
    error_handler(ERRHDL_ERROR, 
                  "channel_read_network: unable to open file \"%s\": %s",
                  file, strerror(errno));
    return NULL;
  }

  done = FALSE;
  while (!done) {
    int i;
    
    done = (table_get_fields(fields, chan_fields) < 0);

    if (done) {
      for (i = 0; i < fields; i++) {
        if (chan_fields[i].read) break;
      }
      if (i >= fields) continue;
    }
      
    if (head == NULL) {
      head = alloc_channel_segment();
      current = head;
    } else {
      current->next = alloc_channel_segment();
      current = current->next;
    }
    
      
    for (i = 0; i < fields; i++) {
      if (chan_fields[i].read) {
        switch (i) {
        case 0:
          current->id = chan_fields[i].value.integer;
          if (current->id <= 0) {
            error_handler(ERRHDL_ERROR,
                          "%s: segment %d: channel id invalid",
                          file, current->id);
            err++;
          }
        break;
        case 1:
          if (chan_fields[i].value.integer > 0) {
            current->order = chan_fields[i].value.integer;
          } else {
            error_handler(ERRHDL_ERROR,
                          "%s: segment %d: channel order (%d) invalid",
                          file, current->id, chan_fields[i].value.integer);
            err++;
          }
          break;
        case 2:
          if(chan_fields[i].value.real > 0) {
            current->slope = chan_fields[i].value.real;
          } else {
            error_handler(ERRHDL_ERROR,
                          "%s: segment %d: channel slope (%f) invalid",
                          file, current->id, chan_fields[i].value.real);
            err++;
          }
          break;
        case 3:
          if(chan_fields[i].value.real > 0) {
            current->length = chan_fields[i].value.real;
          } else {
            error_handler(ERRHDL_ERROR,
                          "%s: segment %d: channel length (%f) invalid",
                          file, current->id, chan_fields[i].value.real);
            err++;
          }
        break;
        case 4:
          if ((current->channel_class = 
               find_channel_class(class_list,
                                  chan_fields[i].value.integer)
               ) == NULL) {
            error_handler(ERRHDL_ERROR,
                          "%s: segment %d: channel class %d not found",
                          file, current->id, 
                          chan_fields[i].value.integer);
            err++;
          }
          break;
        case 5:
          current->outlet = (Channel *)chan_fields[i].value.integer;
          break;
        case 6:
          current->record = TRUE;
          break;
        case 7:
          current->record_name = (char *)strdup(chan_fields[i].field);
          break;
        default:
          error_handler(ERRHDL_FATAL,
                        "channel_read_network: what is this field %d?",
                        i);
          break;
        }
      }
    }
  }

  table_close();

                                /* find segment outlet segments, if
                                   specified */

  for (current = head; current != NULL; current = current->next) {
    int outid = (int) current->outlet;

    if (outid != 0) {
      current->outlet = channel_find_segment(head, outid);
      if (current->outlet == NULL) {
        error_handler(ERRHDL_ERROR,
                      "%s: cannot find outlet (%d) for segment %d",
                      file, outid, current->id);
        err++;
      }
    }
  }
  
  table_errors += err;

  error_handler(ERRHDL_STATUS,
                "channel_read_network: %s: %d errors, %d warnings",
                file, table_errors, table_warnings);

  
  
  if (table_errors) {
    error_handler(ERRHDL_ERROR,
                  "channel_read_network: %s: too many errors",
                  file);
    channel_free_network(head);
    head = NULL;
  }

  return(head);
}

/* -------------------------------------------------------------
   channel_route_segment
   ------------------------------------------------------------- */
static int
channel_route_segment(Channel *segment, int deltat, int ChannelRouting) 
{
  float inflow, outflow, lateral_inflow, storage;
  float roughness, width, startDepth;
  float K = segment->K;
  float X = segment->X;
  int err = 0;

  /* change masses to rates */
  inflow = segment->inflow/deltat;
  lateral_inflow = segment->lateral_inflow/deltat;


  /* linear routing -------------------------------------------------------- */
  if( ChannelRouting == LINEAR ){
    storage = ( (inflow+lateral_inflow) / K) + 
                (segment->storage - (inflow+lateral_inflow)/K)*X;
    if( storage < 0.0 )
      storage = 0.0;
    outflow = (inflow+lateral_inflow) - (storage - segment->storage)/deltat;
  
    segment->outflow = outflow*deltat;
    segment->storage = storage;
  } else if ( ChannelRouting == POWERLAW ) { /* nonlinear routing with power law method */
      roughness = segment->channel_class->friction;
      width = segment->channel_class->width;
      startDepth = segment->storage / (segment->length * width);

      segment->outflow = route_segment_powerlaw(inflow + lateral_inflow, (float)deltat, startDepth, 
                                     width, segment->length, 
                                     segment->slope, roughness, &(segment->storage));
  } else if ( ChannelRouting == FTABLE ) { /* nonlinear routing with FTABLE method */
      segment->outflow = route_segment( segment->id, (segment->inflow + segment->lateral_inflow), 
                                        (float)deltat, &(segment->recharge) );
  } else printf("Error in channel routing option--must be LINEAR, POWERLAW, or FTABLE.");

  /* end routing --------------------------------------------------------------------*/
 
  if (segment->outlet != NULL) 
    segment->outlet->inflow += segment->outflow;

  return(err);


  /* for old M-C code ----------------------------------------------------------
  float C1, C2, C3, C4; */


  /* The following lines are currently not used and have therefore been
     commented out, Bart Nijssen, Wed Feb  3 19:06:33 1999.  Note that in this
     sections the units for C1, C2, C3, and C4 need to be checked before using
     the code.  The units of deltat have been changed from hours to seconds, but 
     that might not yet be reflected correctly in the next couple of lines */
  
  /* last_inflow = segment->last_inflow/deltat; */
  /* last_outflow = segment->last_outflow/deltat; */
  /* last_lateral_inflow = segment->last_lateral_inflow/deltat; */

  /*   C1 = deltat - 2*K*X; */
  /*   C2 = deltat + 2*K*X; */
  /*   C3 = 2*K*(1 - X) - deltat; */
  /*   C4 = lateral_inflow*deltat; */
  
  /*   outflow = ( C1*inflow + C2*last_inflow + C3*last_outflow )/  */
  /*     (2*K*(1 - X) + deltat); */
  
  /*   storage = K*(X*inflow + (1 - X)*outflow); */


}


/* -------------------------------------------------------------
   channel_route_network
   ------------------------------------------------------------- */
int
channel_route_network(Channel *net, int deltat, int ChannelRouting)
{
  unsigned int order;
  int order_count;
  int err = 0;
  Channel *current;
  
  for (order = 1; ; order += 1) {
    order_count = 0;
    current = net;
    while (current != NULL) {
      if (current->order == order) {
        err += channel_route_segment(current, deltat, ChannelRouting);
        order_count += 1;
      }
      
/* this next statement gives good values */
/* if(current->recharge > 0.0) printf("r3a=%f ", current->recharge); */

      current = current->next;

/* this next statement causes segmentation fault */
/* if(current->recharge > 0.0) printf("r3b=%f ", current->recharge); */

    }

    if (order_count == 0) break;
  }

  return (err);
}

/* -------------------------------------------------------------
   channel_step_initialize_network
   ------------------------------------------------------------- */
int 
channel_step_initialize_network(Channel *net)
{
  if (net != NULL) {
    net->last_inflow = net->inflow;
    net->inflow = 0.0;
    net->lateral_inflow = 0.0;
    net->quick_flow = 0.0;
    net->last_outflow = net->outflow;
    net->last_storage = net->storage;
    net->last_lateral_inflow = net->lateral_inflow;
    channel_step_initialize_network(net->next);
  }
  return (0);
}

/* -------------------------------------------------------------
   channel_save_outflow
   This routine saves the channel output
   ------------------------------------------------------------- */
int
channel_save_outflow(double time, Channel *net, FILE *out, FILE *out2)
{
  char buffer[16];
  sprintf(buffer, "%12.5g", time);
  return (channel_save_outflow_text(buffer, net, out, out2,0));
}
  
/* -------------------------------------------------------------
   channel_save_outflow_wtext
   Saves the channel outflow using a text string as the time field
   ------------------------------------------------------------- */
int
channel_save_outflow_text(char *tstring, Channel *net, FILE *out, FILE *out2, int flag)
{
  int err = 0;
  float total_outflow = 0.0;
  float total_lateral_inflow = 0.0;
  float total_storage = 0.0;
  float total_storage_change = 0.0;
  float total_error = 0.0;
  FILE *NewOutFile;
  char name[255];
  int NumSegs = 0;
  Channel *head_net;

  head_net = net;


  if(flag==1){ 
    for (; net != NULL; net = net->next) {
      if(net->record) NumSegs++;
    }
    fprintf(out2, " %d\n", NumSegs );    

    /* test */ net = head_net;

    fprintf(out2,"DATE ");
    for (; net != NULL; net = net->next) {
    total_lateral_inflow += net->lateral_inflow;
    if (net->outlet == NULL) {
      total_outflow += net->outflow;
    }
    if(net->record) fprintf(out2,"%s ",net->record_name);
    }
    fprintf(out2,"\n");

    /* test */ net = head_net;

  }

     if (fprintf(out2, "%15s ", tstring) == EOF) {
        error_handler(ERRHDL_ERROR, 
                      "channel_save_outflow: write error:%s",
                      strerror(errno));
        err++;
      }  

  for (; net != NULL; net = net->next) {
    total_lateral_inflow += net->lateral_inflow;
    if (net->outlet == NULL) {
      total_outflow += net->outflow;
    }
    total_storage += net->storage;
    total_storage_change += net->storage - net->last_storage;

/*     printf("r4=%.1f ", net->recharge); */

    if (net->record) {  
      if (fprintf(out, "%15s %10d %12.5g %12.5g %12.5g %12.5g",
                  tstring, net->id, net->inflow, net->lateral_inflow, 
                  net->outflow, net->storage - net->last_storage) == EOF) {
        error_handler(ERRHDL_ERROR, 
                      "channel_save_outflow: write error:%s",
                      strerror(errno));
        err++;
      }
      if (fprintf(out2, "%12.5g ", net->outflow) == EOF) {
        error_handler(ERRHDL_ERROR, 
                      "channel_save_outflow: write error:%s",
                      strerror(errno));
        err++;
      }
      if (net->record_name != NULL) {
        if (fprintf(out, "   \"%s\"\n", net->record_name) == EOF) {
          error_handler(ERRHDL_ERROR, 
                        "channel_save_outflow: write error:%s",
                        strerror(errno));
          err++;
        }

      } else {
        if (fprintf(out, "\n") == EOF) {
          error_handler(ERRHDL_ERROR, 
                        "channel_save_outflow: write error:%s",
                        strerror(errno));
          err++;
        }
      }
    }
  }
  total_error = total_storage_change - total_lateral_inflow + total_outflow;
  if (fprintf(out, "%15s %10d %12.5g %12.5g %12.5g %12.5g %12.5g \"Totals\"\n",
              tstring, 0, total_lateral_inflow, 
              total_outflow, total_storage, 
              total_storage_change, total_error) == EOF) {
    error_handler(ERRHDL_ERROR, 
                  "channel_save_outflow: write error:%s",
                  strerror(errno));
    err++;
  }
  fprintf(out2,"\n");
  
  return (err);
    }

/* -------------------------------------------------------------
   channel_free_network
   ------------------------------------------------------------- */
void
channel_free_network(Channel *net) 
{
  if (net->next != NULL) {
    channel_free_network(net->next);
  }
  free(net);
}

/* -------------------------------------------------------------
   channel_init
   ------------------------------------------------------------- */
void
channel_init(void)
{
  /* do nothing */
  return;
}


/* -------------------------------------------------------------
   channel_done
   ------------------------------------------------------------- */
void
channel_done(void)
{
  /* do nothing */
  return;
}



#if TEST_MAIN

/* -------------------------------------------------------------
   interpolate
   ------------------------------------------------------------- */
static float
interpolate(int n, float *x, float *y, float x0)
{
  int i;
  if (x0 <= x[0]) {
    return ( (x0 - x[0])/(x[1] - x[0])*(y[1] - y[0]) + y[0] );
  }
  for (i = 0; i < n - 1; i++) {
    if (x0 < x[i + 1]) {
      return ( (x0 - x[i])/(x[i + 1] - x[i])*(y[i + 1] - y[i]) + y[i] );
    }
  }
  return ( (x0 - x[i - 1])/(x[i] - x[i - 1])*(y[i] - y[i - 1]) + y[i] );
}
  
    
/* -------------------------------------------------------------
   Main Program
   ------------------------------------------------------------- */
int
main(int argc, char **argv)
{
  static int interval = 3600;	/* timestep in seconds */
  static int timestep;
  static int endtime  = 144;
#define TIMES 6
  static float bndflow[TIMES] = { 0.0, 0.0, 300.0, 300.0, 0.0, 0.0 };
  static float bndtime[TIMES] = { 0.0, 12.0, 36.0, 48.0, 60.0, 1000.0 };

  float time;
  ChannelClass *class;
  Channel *simple = NULL, *current, *tail;

  error_handler_init(argv[0], NULL, ERRHDL_ERROR);
  channel_init();

                                /* read classes */

  if ((class = channel_read_classes("example_classes.dat")) == NULL) {
    error_handler(ERRHDL_FATAL, 
                  "example_classes.dat: trouble reading file");
  }

                                /* read a network */

  if ((simple = channel_read_network("example_network.dat", class)) == NULL) {
    error_handler(ERRHDL_FATAL,
                  "example_network.dat: trouble reading file");
  }
    
                                /* initialize flows */
    
  for (current = simple; current != NULL; current = current->next) {
    current->inflow = bndflow[0];
    current->outflow = bndflow[0];
    current->outlet = current->next;
    tail = current;
  }

                                /* time loop */

  for (timestep = 0; timestep <= endtime; timestep++) {
    float inflow = interpolate(TIMES, bndtime, bndflow, timestep)*interval;
    float outflow;
    
    channel_step_initialize_network(simple);
    simple->inflow = inflow;
    (void) channel_route_network(simple, interval, 1);
    outflow = tail->outflow/interval;
    channel_save_outflow(timestep*interval, simple, stdout);
  }
  
  channel_free_network(simple);
  channel_free_classes(class);
  channel_done();
  error_handler_done();
  exit(0);
}
#endif  


