hydrorain.c

/*-------------------------------------------------------------------------------------------
 * hydrorain.c
 *
 *	Author: Albert Kettner, March 2006
 *
 * 	Calculates the daily rain fall and rain derived runoff for each altitude bin. Also
 *	calculates the rain derived groundwater flow, evaporation and time lag.
 *
 *	Physically based precipitation/groundwater/evaporation routine based on:
 *	M. Sivapalan, J.K. Ruprecht, N.R. Viney, 1996.  Water and Salt Balance modelling to
 *	predict the effects of Land-Use changes in forested catchments. 
 *	1. Small Catchment water balance model. Hydrological Processes, V.10, P.393-411.
 *
 * Variable		Def.Location	Type	Units	Usage
 * --------		------------	----	-----	-----
 * Mannualin	HydroRain.c		double	m^3/a	annual mass input to rain routine
 * Mannualout	HydroRain.c		double	m^3/a	annual mass output from rain routine
 * Minput		HydroRain.c		double	m^3		daily mass input to rain routine
 * Mout			HydroRain.c		double	m^3		daily mass input to rain routine
 * Mwrap		HydroRain.c		double	m^3		mass input to rain routine from previous year
 * Pground		HydroRain.c		double	m/day	ground precipitation
 * Qi1, Qi2		HydroRain.c		double	m^3/s	temporary variables for Qinfiltration
 * Qinfiltexcess HydroRain.c	double	m^3/s	river discharge from infiltration excess
 * Qinfiltration HydroRain.c	double	m^3/s	= Pground-Qrain[ii]
 * Qrainlag		HydroRain.c		double	m^3/s	temporary array for routing
 * Qsaturationexcess HydroRain.c double	m^3/s	river discharge from saturation excess
 * Qsslag		HydroRain.c		double	m^3/s	temporary array for routing
 * Qtogw		HydroRain.c		double	m^3/s	Rain derived discharge to the GW pool
 * err			various			int		-		error flag, halts program
 * ii			various			int		-		temporary loop counter
 * infiltrate	HydroRain.c		double	m/day	infiltration capacity
 * jj			various			int		-		temporary loop counter
 * kk			various			int		-		temporary loop counter
 * ratio1		HydroRain.c		double	-		ratio of the fullness of the GW pool
 * ratio2		HydroRain.c		double	-		linear interpolation for infiltration
 * shldday[]	HydroRain.c		double	m^3/s	shoulder discharge array
 *
 *-------------------------------------------------------------------------------------------*/

#include "hydroclimate.h"
#include "hydroparams.h"
#include "hydrotimeser.h"
#include <stdlib.h>

int FLAindex[daysiy];
long *daysievent;

/*-------------------------------
 *  Daily suspended and bedload
 *-------------------------------*/
double Cs[daysiy], **Csoutlet, Qb[daysiy], **Qboutlet, Qs[daysiy],
  Qspsi[daysiy], Qsglacier[daysiy], **Qsoutlet;

/*---------------------------------------------------------
 *  Daily Temperature, Precipitation and Snow time series
 *---------------------------------------------------------*/
double rainarea[daysiy], Ecanopy[daysiy];
double Pdaily[daysiy], Tdaily[daysiy], **Snowelevday;

/*--------------------------------------------------------------------
 *  Arrays for each type of daily discharge (rain,snow,ice,total,GW)
 *--------------------------------------------------------------------*/
double Qrain[maxday], Qice[maxday], Qnival[maxday], Qsumtot[maxday], **Qsum,
  Qss[maxday];

/*--------------------------------------------------
 *  Arrays for carryover from one year to the next
 *--------------------------------------------------*/
double Qrainwrap[wrapday], Qicewrap[wrapday], Qnivalwrap[wrapday];
double Qsswrap[wrapday], *Snowcarry;

/*-------------------------------------------
 *  Arrays for the groundwater storage pool
 *-------------------------------------------*/
double gwstore[daysiy], Qicetogw[daysiy], Qnivaltogw[daysiy];
double Egw[daysiy], Qexceedgw[daysiy];


/*----------------------
 *  Start of HydroRain
 *----------------------*/
int
hydrorain ()
{

/*-------------------
 *  Local Variables
 *-------------------*/
  double infiltrate, Pground[daysiy];
  double Qsaturationexcess[daysiy], Qinfiltexcess[daysiy];
  double Qinfiltration[daysiy], Qi1, Qi2, Qtogw, shldday[maxday];
  double Mout, Mannualin, Mannualout;
  double Qrainlag[maxday], Qsslag[maxday];
  double ratio1, ratio2;
  double test;
  int ii, jj, kk, err;

/*------------------------
 *  Initialize Variables
 *------------------------*/
  err = 0;
  MPrain = 0.0;
  Mannualin = 0.0;
  Mannualout = 0.0;
  for (ii = 0; ii < maxday; ii++)
    {
      Qrainlag[ii] = 0.0;
      Qsslag[ii] = 0.0;
      shldday[ii] = 0.0;
    }

/*-------------------------
 *  Loop through the year
 *-------------------------*/
  for (ii = 0; ii < daysiy; ii++)
    {

   /*----------------------------------------------------------------------
    *  Find the area being rained on
    *  Temperature must be above freezing
    *    All rain (P with T>=0.0) is handled here.
    *    All snow/ice (P with T<0.0) is handled in HydroSnow/HydroGlacier
    *----------------------------------------------------------------------*/
      rainarea[ii] = 0.0;
      for (kk = 0; kk < nelevbins && kk < FLAindex[ii]; kk++)
        rainarea[ii] += areabins[kk];

   /*--------------------------------------------------------------
    *  Get the remaining GW pool from the previous day
    *  the first day of the year is taken care of in HydroTrend.c
    *--------------------------------------------------------------*/
      if (ii > 0)
        gwstore[ii] = gwstore[ii - 1];

   /*---------------------------------------------------------------------------
    *  Calculate groundwater discharge (subsurface storm flow) and evaporation
    *  Qss (m^3/s), Egw (m/day), gwstore (m^3)
    *  Only drain down to gwmin
    *
    *  Having this first builds in a one delay lag to the GW discharge
    *  This also helps prevent the Glacier and Snow from totally filling
    *  the pool, and leaves some room for Rain.
    *
    *  The lines E=max(E,0.0) and Qss=max(Qss,0.0) prevent rare errors due
    *  to very small numbers when gwstore is close to gwmin.
    *---------------------------------------------------------------------------*/
      Egw[ii] =
        mn (alphagwe[ep] *
            pow ((gwstore[ii] - gwmin[ep]) / (gwmax[ep] - gwmin[ep]),
                 betagwe[ep]), (gwstore[ii] - gwmin[ep]) / totalarea[ep]);
      Egw[ii] = mx (Egw[ii], 0.0);
      gwstore[ii] -= Egw[ii] * totalarea[ep];

      Qss[ii] =
        mn (alphass[ep] *
            pow ((gwstore[ii] - gwmin[ep]) / (gwmax[ep] - gwmin[ep]),
                 betass[ep]), (gwstore[ii] - gwmin[ep]) / dTOs);
      Qss[ii] = mx (Qss[ii], 0.0);
      gwstore[ii] -= Qss[ii] * dTOs;

   /*-------------------------------------------------------------------------------
    *  Add the other sources to the Groundwater pool
    *  allocate overflow to Qexceedgw and add to river discharge in HydroSumFlow.c
    *-------------------------------------------------------------------------------*/
      if (gwstore[ii] < gwmax[ep])
        Qtogw =
          mn (Qicetogw[ii] + Qnivaltogw[ii], (gwmax[ep] - gwstore[ii]) / dTOs);
      else
        Qtogw = 0.0;

      gwstore[ii] += Qtogw * dTOs;
      Qexceedgw[ii] = Qicetogw[ii] + Qnivaltogw[ii] - Qtogw;

   /*--------------------------------------
    *  Calculate commonly used quantities
    *--------------------------------------*/
      ratio1 = (gwstore[ii] - gwmin[ep]) / (gwmax[ep] - gwmin[ep]);
      ratio2 = (Ko[ep] - pcr) / (pmax - pcr);

   /*------------------------------------------------------------
    *  Calculate the Ground Precipitation (m/day)
    *  the remainder is evaporated from canopy interception
    *  keep track of evaporation (m/day) for mass balance check
    *------------------------------------------------------------*/

        /*--------------------------------------------------------
	 *  Check if there is a input file with evaporation data
	 *--------------------------------------------------------*/
      if (evaporationdatafile == 0)
        {
          Pground[ii] = mx (alphag[ep] + betag[ep] * Pdaily[ii], 0.0);
//              printf("%f, %f\n",Pground[ii], Pdaily[ii]);
        }
      Ecanopy[ii] = Pdaily[ii] - Pground[ii];
//   printf("%f, %f, %f\n",Ecanopy[ii], Pground[ii], betag[ep]);

   /*----------------------------------------------------
    *  Calculate Surface Runoff (m^3/s)
    *  Sum of saturation excess and infiltration excess
    *  remainder goes to groundwater pool
    *----------------------------------------------------*/

   /*-----------------------------------------------------	
    *  Calculate the Saturation Excess
    *  this is all the rain that falls on saturated soil
    *  and cannot go into the groundwater pool
    *-----------------------------------------------------*/
      if (gwstore[ii] > gwmin[ep])
        Qsaturationexcess[ii] =
          (Pground[ii] * rainarea[ii] / dTOs) * alphac * pow (ratio1, betac);
      else
        Qsaturationexcess[ii] = 0.0;

   /*--------------------------------------------------------------------
    *  Calculate Theoretical Infiltration Rate based on Precip	(m/day)
    *    a) below pcr, 100% infiltration
    *    b) above pmax, max hydr. cond.
    *    c) linear fit in between
    *--------------------------------------------------------------------*/
      if (Pground[ii] < pcr)
        infiltrate = Pground[ii];
      else if (Pground[ii] > pmax)
        infiltrate = Ko[ep];
      else
        infiltrate = Pground[ii] * ratio2 + pcr * (1 - ratio2);

   /*-------------------------------------------------------------------------
    *  Calculate infiltration "discharge" to the groundwater pool (m^3/s)
    *  1st) calculate the theoretical infiltration discharge
    *  2nd) calculate the discharge that fills the groundwater pool in a day
    *  take the minimum so as not to overfill the pool
    *-------------------------------------------------------------------------*/
      if (gwstore[ii] > gwmin[ep])
        Qi1 = infiltrate * (1.0 - ratio1) * rainarea[ii] / dTOs;
      else
        Qi1 = infiltrate * rainarea[ii] / dTOs;
      Qi2 = (gwmax[ep] - gwstore[ii]) / dTOs;

      Qinfiltration[ii] = mn (Qi1, Qi2);
      gwstore[ii] += Qinfiltration[ii] * dTOs;

   /*--------------------------------------------------
    *  Calculate the infiltration excess (m^3/s)
    *  this is the runoff that:
    *    a) is not due to saturation excess
    *    b) and can not get into the groundwater pool
    *--------------------------------------------------*/
      Qinfiltexcess[ii] = mx (0.0,
                              Pground[ii] * rainarea[ii] / dTOs -
                              Qsaturationexcess[ii] - Qinfiltration[ii]);

   /*-------------------------------------------------------
    *  Calculate Surface Runoff (m^3/s)
    *    Sum of saturation excess and infiltration excess
    *    remainder went to groundwater pool or evaporation
    *-------------------------------------------------------*/
      Qrain[ii] = Qsaturationexcess[ii] + Qinfiltexcess[ii];
//   printf("%f\n",Qrain[ii]);

   /*---------------------------------------------------------------------------
    *  Add the carryover from the previous year
    *
    *  NOTE: adding at this stage assumes the wrapped rain is due to the wave
    *  propagating down the river and that this water does not have additional
    *  contributions to the Evap or gw pools
    *---------------------------------------------------------------------------*/
      if (ii < wrapday)
        Qrain[ii] += Qrainwrap[ii];

   /*------------------------------------------------
    *  Add to the daily rain mass balance (m^3/day)
    *------------------------------------------------*/
      Minput = Pdaily[ii] * rainarea[ii];
      if (ii < maxday - daysiy)
        Minput += Qrainwrap[ii] * dTOs;
      Mout =
        Ecanopy[ii] * rainarea[ii] + (Qrain[ii] + Qinfiltration[ii]) * dTOs;

   /*-----------------------------------------------
    *  Add to the annual rain mass balance (m^3/a)
    *-----------------------------------------------*/
      MPrain += Pdaily[ii] * rainarea[ii];
      Mannualin += Minput;
      Mannualout += Mout;

   /*---------------------------------------------
    *  Check the daily rain mass balance (m^3/d)
    *---------------------------------------------*/
      if (Minput > 0. && fabs (Mout - Minput) / Minput > masscheck)
        {
          fprintf (stderr, "ERROR: in HydroRain.c: \n");
          fprintf (stderr, "   Daily Mass Balance Error: Mout != Minput \n\n");
          fprintf (stderr, "   fabs(Mout-Minput)/Minput > masscheck \n");
          fprintf (stderr, "   note: masscheck set in HydroParams.h \n");
          fprintf (stderr, "   masscheck = %e (m^3) \n", masscheck);
          fprintf (stderr, "   fabs(Mout-Minput)/Minput = %e (m^3) \n\n",
                   fabs (Mout - Minput) / Minput);
          fprintf (stderr, " \t Minput = Pdaily + Qrainwrap \n");
          fprintf (stderr, " \t Minput \t\t = %e (m^3) \n", Minput);
          fprintf (stderr, " \t Pdaily[ii]*rainarea[ii] \t = %e (m^3) \n",
                   Pdaily[ii] * rainarea[ii]);
          if (ii < maxday - daysiy)
            fprintf (stderr, " \t Qrainwrap[ii] \t\t = %e (m^3) \n\n",
                     Qrainwrap[ii] * dTOs);
          else
            fprintf (stderr, " \t Qrainwrap[ii] \t\t = 0.0 (m^3) \n\n");
          fprintf (stderr, " \t Mout = Ecanopy + Qrain + Qinfiltration \n");
          fprintf (stderr, " \t Mout \t\t = %e (m^3) \n", Mout);
          fprintf (stderr, " \t Ecanopy[ii]*rainarea[ii] \t\t = %e (m^3) \n",
                   Ecanopy[ii] * rainarea[ii]);
          fprintf (stderr, " \t Qrain[ii]*dTOs \t\t = %e (m^3) \n",
                   Qrain[ii] * dTOs);
          fprintf (stderr, " \t Qinfiltration[ii]*dTOs \t = %e (m^3) \n\n",
                   Qinfiltration[ii] * dTOs);
          fprintf (stderr, " \t Day \t d \t %d \n", ii + 1);
          fprintf (stderr, " \t Year \t d \t %d \n", yr + 1);
          fprintf (stderr, " \t Epoch \t - \t %d \n", ep + 1);
          exit (-1);
        }

   /*------------------------------------------------
    *  Check for over/under-filled groundwater pool
    *  return error
    *------------------------------------------------*/
      if (gwstore[ii] < gwmin[ep] / 1.001 || gwstore[ii] > gwmax[ep] * 1.001)
        {
          fprintf (stderr, "WARNING in HydroRain.c \n");
          fprintf (stderr, "   gwstore out of range. \n");
          fprintf (stderr, "   criteria:  gwmin < gwstore < gwmax \n");
          fprintf (stderr, "      gwmax   (m^3) = %e \n", gwmax[ep]);
          fprintf (stderr, "      gwstore (m^3) = %e \n", gwstore[ii]);
          fprintf (stderr, "      gwmin   (m^3) = %e \n", gwmin[ep]);
          fprintf (stderr, "      Qss[ii] (m^3) = %e \n", Qss[ii] * dTOs);
          fprintf (stderr, "      Egw[ii] (m^3) = %e \n",
                   Egw[ii] * totalarea[ep]);
          err = 1;
        }
    }                           /* end day loop */

/*------------------------------------------------------------
 *  Use shoulder on Rain events if the basin is large enough
 *------------------------------------------------------------*/
  if ((int) (basinlength[ep] / (avgvel[ep] * dTOs)) >= 3 || Rvol[ep] > 0.0)
    {

   /*------------------------------------------------------------
    *  Create the shoulder events (Murray's version of routing)
    *    there is one left (preceeding) day scaled as:
    *      shoulderleft*event
    *    the main event is scaled down to:
    *      shouldermain*event
    *    there are 1 or more right (following days) scaled to:
    *      shoulderright[]*event
    *
    *    1.0 = Sum(shoulderleft+shouldermain+shoulderright[])
    *------------------------------------------------------------*/
      ii = 0;
      if (Qrain[ii] > 0.0)
        {
          shldday[ii] += shoulderleft * Qrain[ii];
          for (jj = 0; jj < shouldern - 2; jj++)
            shldday[ii + jj + 1] += shoulderright[jj] * Qrain[ii];
          Qrain[ii] = shouldermain * Qrain[ii];
        }
      for (ii = 1; ii < daysiy; ii++)
        {
          if (Qrain[ii] > 0.0)
            {
              shldday[ii - 1] += shoulderleft * Qrain[ii];
              for (jj = 0; jj < shouldern - 2; jj++)
                shldday[ii + jj + 1] += shoulderright[jj] * Qrain[ii];
              Qrain[ii] = shouldermain * Qrain[ii];
            }
        }

   /*----------------------------------------------------------------
    *  Add the shoulder events and the main events to get the total
    *  rain derived discharge
    *----------------------------------------------------------------*/
      for (ii = 0; ii < maxday; ii++)
        Qrain[ii] += shldday[ii];
    }                           /* endif shoulder */

/*---------------------------------------------------------------------
 *  (Mark's version of routing)
 *  Remove any lagged discharge from the initial rain discharge event
 *---------------------------------------------------------------------*/
  for (ii = 0; ii < daysiy; ii++)
    {
      for (kk = 0; kk < nelevbins && kk < FLAindex[ii]; kk++)
        {
          if (rainarea[ii] > 0.0)
            {
              Qrainlag[ii + distbins[kk]] +=
                Qrain[ii] * areabins[kk] / rainarea[ii];
              Qsslag[ii + distbins[kk]] +=
                Qss[ii] * areabins[kk] / rainarea[ii];
              Qrain[ii] -= Qrain[ii] * areabins[kk] / rainarea[ii];
              Qss[ii] -= Qss[ii] * areabins[kk] / rainarea[ii];
            }
        }
    }

/*------------------------------------
 *  Add the lagged discharge back in
 *------------------------------------*/
  for (ii = 0; ii < maxday; ii++)
    {
      Qrain[ii] += Qrainlag[ii];
      Qss[ii] += Qsslag[ii];
    }

/*--------------------------------
 *  Add the GW wrapped discharge
 *--------------------------------*/
  for (ii = 0; ii < wrapday; ii++)
    Qss[ii] += Qsswrap[ii];

/*---------------------------------------
 *  Check the annual mass balance (m^3)
 *---------------------------------------*/
  if (fabs (Mannualout - Mannualin) / Mannualin > masscheck)
    {
      fprintf (stderr, "ERROR: in HydroRain.c: \n");
      fprintf (stderr, "   Annual Mass Balance Error: Mout != Minput \n\n");

      fprintf (stderr,
               "   fabs(Mannualout-Mannualin)/Mannualin > masscheck \n");
      fprintf (stderr, "   note: masscheck set in HydroParams.h \n");
      fprintf (stderr, "   masscheck = %e (m^3) \n", masscheck);
      fprintf (stderr,
               "   fabs(Mannualout-Mannualin)/Mannualin = %e (m^3) \n\n",
               fabs (Mannualout - Mannualin) / Mannualin);

      fprintf (stderr,
               " \t Mannualin = Sum(Pdaily+Qrainwrap) \t\t = %e (m^3) \n",
               Mannualin);
      fprintf (stderr,
               " \t Mannualout = Sum(Ecanopy+Qrain+Qinfiltration) \t = %e (m^3) \n\n",
               Mannualout);

      fprintf (stderr, " \t Year \t d \t %d \n", yr + 1);
      fprintf (stderr, " \t Epoch \t - \t %d \n\n", ep + 1);

      exit (-1);
    }
  return (err);
}                               /* end of HydroRain.c */