/*
 * SUMMARY:      Aggregate.c - calculate basin-wide values
 * 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 Nov 27 11:48:33 2002 by Scott Waichler <waichler@tuff.pnl.gov>
 * DESCRIPTION:  Calculate the average values for the different fluxes and
 *               state variables over the basin.
 * DESCRIP-END.
 * FUNCTIONS:    Aggregate()
 * COMMENTS:
 */

#ifndef lint
static char vcid[] = "$Id: Aggregate.c,v 1.3 2001/04/04 16:43:45 waichler Exp $";
#endif /* lint */

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

/*****************************************************************************
  Aggregate()
  
  Calculate the average values for the different fluxes and state variables
  over the basin.  Only the runoff is calculated as a total volume instead
  of an average.  In the current implementation the local radiation
  elements are not stored for the entire area.  Therefore these components
  are aggregated in AggregateRadiation() inside MassEnergyBalance().

  The aggregated values are set to zero in the function RestAggregate,
  which is executed at the beginning of each time step.
*****************************************************************************/
void Aggregate(int SpatialExtent, MAPSIZE *Map, OPTIONSTRUCT *Options, TOPOPIX **TopoMap, 
	       LAYER *Soil, LAYER *Veg, VEGPIX **VegMap, EVAPPIX **Evap, 
	       PRECIPPIX **Precip, RADCLASSPIX **RadMap, SNOWPIX **Snow, 
	       SOILPIX **SoilMap, AGGREGATED *Total, VEGTABLE *VType, 
	       ROADSTRUCT **Network, GWPIX **Groundwater, FRZPIX **FrozenSoil)
{
  int NPixels;                  /* Number of pixels in the basin */
  int NSoilL;                   /* Number of soil layers for current pixel */
  int NVegL;                    /* Number of vegetation layers for current
                                   pixel */
  int i;                        /* counter */
  int j;                        /* counter */
  int x;
  int y;
  float DeepDepth;		/* depth to bottom of lowest rooting zone */
  double Area = 1;          /* relative grid cell area  */
  double TotalArea = 0;     /* accumulated area m^2 */
  double AverageArea = 1;

  NPixels = 0;

  /* set up 
  for (y = 0; y < Map->NY; y++) 
    for (x = 0; x < Map->NX; x++) 
      if (InBasin(TopoMap[y][x].Mask)) {
        TotalArea += TopoMap[y][x].Area;
        NPixels++;
      }
*/
  TotalArea = Total->BasinArea;
  NPixels = Total->ActiveCells;

  AverageArea = TotalArea / (float)NPixels;

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

        /* use relative area to avoid roundoff errors when accumulating */
        Area = TopoMap[y][x].Area / AverageArea;

        NSoilL = Soil->NLayers[SoilMap[y][x].Soil-1];
        NVegL = Veg->NLayers[VegMap[y][x].Veg-1];

        /* aggregate the evaporation data */
        Total->Evap.ETot += Evap[y][x].ETot * Area;
        for (i = 0; i < NVegL; i++) {
          Total->Evap.EPot[i] += Evap[y][x].EPot[i] * Area;
          Total->Evap.EAct[i] += Evap[y][x].EAct[i] * Area;
          Total->Evap.EInt[i] += Evap[y][x].EInt[i] * Area;
        }
        Total->Evap.EPot[Veg->MaxLayers] += Evap[y][x].EPot[NVegL] * Area;
        Total->Evap.EAct[Veg->MaxLayers] += Evap[y][x].EAct[NVegL] * Area;
        for (i = 0; i < NVegL; i++) {
          for (j = 0; j < NSoilL; j++) {
            Total->Evap.ESoil[i][j] += Evap[y][x].ESoil[i][j] * Area;
          }
        }
        Total->Evap.EvapSoil += Evap[y][x].EvapSoil * Area;
        
        /* aggregate precipitation data */
        Total->Precip.Precip += Precip[y][x].Precip * Area;
        for (i = 0; i < NVegL; i++){ 
          Total->Precip.IntRain[i] += Precip[y][x].IntRain[i] * Area;
          Total->Precip.IntSnow[i] += Precip[y][x].IntSnow[i] * Area;
          Total->CanopyWater += Precip[y][x].IntRain[i] * Area +
            Precip[y][x].IntSnow[i] * Area;
        }

        /* aggregate radiation data */
	if (Options->MM5 == TRUE) {
	  Total->RadClass.Beam = NOT_APPLICABLE;
	  Total->RadClass.Diffuse = NOT_APPLICABLE;
	}
	else {
	  Total->RadClass.Beam += RadMap[y][x].Beam * Area;
	  Total->RadClass.Diffuse += RadMap[y][x].Diffuse * Area;
	}

        /* aggregate snow data */
        if (Snow[y][x].HasSnow)
          Total->Snow.HasSnow = TRUE;
        Total->Snow.Swq += Snow[y][x].Swq * Area;
        /* Total->Snow.Melt += Snow[y][x].Melt; */
        Total->Snow.Melt += Snow[y][x].Outflow * Area;
        Total->Snow.PackWater += Snow[y][x].PackWater * Area;
        Total->Snow.TPack += Snow[y][x].TPack * Area;
        Total->Snow.SurfWater += Snow[y][x].SurfWater * Area;
        Total->Snow.TSurf += Snow[y][x].TSurf * Area;
        Total->Snow.ColdContent += Snow[y][x].ColdContent * Area;
        Total->Snow.Albedo += Snow[y][x].Albedo * Area;
        Total->Snow.Depth += Snow[y][x].Depth * Area;
        Total->Snow.VaporMassFlux += Snow[y][x].VaporMassFlux * Area;
        Total->Snow.CanopyVaporMassFlux += Snow[y][x].CanopyVaporMassFlux * Area;
        
        /* aggregate soil moisture data */
        Total->Soil.Depth += SoilMap[y][x].Depth;
		DeepDepth = 0.0;
        for (i = 0; i < NSoilL; i++) {
          Total->Soil.Moist[i] += SoilMap[y][x].Moist[i] * Area;
          /* assert(SoilMap[y][x].Moist[i] >= 0.0);  */
          Total->Soil.Perc[i] += SoilMap[y][x].Perc[i] * Area;
          Total->Soil.Temp[i] += SoilMap[y][x].Temp[i] * Area;
          Total->SoilWater +=  SoilMap[y][x].Moist[i] * 
            VType[VegMap[y][x].Veg-1].RootDepth[i] * Network[y][x].Adjust[i] * Area;
          DeepDepth += VType[VegMap[y][x].Veg-1].RootDepth[i];
        }
        Total->Soil.Moist[Soil->MaxLayers] += SoilMap[y][x].Moist[NSoilL] * Area;
        Total->SoilWater += SoilMap[y][x].Moist[NSoilL] *
          (SoilMap[y][x].Depth - DeepDepth) 
          * Network[y][x].Adjust[NSoilL] * Area;
        
        Total->Soil.TableDepth += SoilMap[y][x].TableDepth * Area;
        if (SoilMap[y][x].TableDepth <= 0)
          (Total->Saturated)++;
        Total->Soil.WaterLevel += SoilMap[y][x].WaterLevel * Area;
        Total->Soil.SatFlow += SoilMap[y][x].SatFlow * Area;
        Total->Soil.TSurf += SoilMap[y][x].TSurf * Area;
        Total->Soil.Qnet += SoilMap[y][x].Qnet * Area;
        Total->Soil.Qs += SoilMap[y][x].Qs * Area;
        Total->Soil.Qe += SoilMap[y][x].Qe * Area;
        Total->Soil.Qg += SoilMap[y][x].Qg * Area;
        Total->Soil.Qst += SoilMap[y][x].Qst * Area;
        Total->Soil.GwRecharge += SoilMap[y][x].GwRecharge * Area;
       /*  SoilMap[y][x].GwRecharge = 0.0; */
        if(Options->FrozenSoil)
          Total->frozenSoilHeatFlux += FrozenSoil[y][x].heatFlux * Area;
        
        if( Options->Groundwater ){
          Total->GwDepth += Groundwater[y][x].depth * Area;
          Total->GwDeepLoss += Groundwater[y][x].deepLoss;  /* already in m^3 */
        }
        else{
          Total->GwDepth = 0;
          Total->GwDeepLoss = 0;
        }
        
        Total->Soil.QuickVol += SoilMap[y][x].QuickVol * Area;

		/* For point mode simulations.
           Special treatment here and in UnsaturatedFlow.c, MassEnergyBalance.c, soilmoisture.h  */
		if (SpatialExtent == POINT) { 
			SoilMap[y][x].ChannelInt += SoilMap[y][x].Runoff;
			SoilMap[y][x].Runoff = 0.0;
		}

        Total->Runoff += SoilMap[y][x].Runoff * Area;
        Total->ChannelInt += SoilMap[y][x].ChannelInt * Area;
        SoilMap[y][x].ChannelInt = 0.0;
        Total->RoadInt += SoilMap[y][x].RoadInt * Area;
        SoilMap[y][x].RoadInt = 0.0;
      }
    }
  }

  /* calculate average values for all quantities except the surface flow */
  
  /* average evaporation data */
  Total->Evap.ETot /= NPixels;
  for (i = 0; i < Veg->MaxLayers + 1; i++) {
    Total->Evap.EPot[i] /= NPixels;
    Total->Evap.EAct[i] /= NPixels;
  }
  for (i = 0; i < Veg->MaxLayers; i++) 
    Total->Evap.EInt[i] /= NPixels;
  for (i = 0; i < Veg->MaxLayers; i++) {
    for (j = 0; j < Soil->MaxLayers; j++) {
      Total->Evap.ESoil[i][j] /= NPixels;
    }
  }
  Total->Evap.EvapSoil /=NPixels;;
  
  /* average precipitation data */
  Total->Precip.Precip /= NPixels;
  for (i = 0; i < Veg->MaxLayers; i++) {
    Total->Precip.IntRain[i] /= NPixels;
    Total->Precip.IntSnow[i] /= NPixels;
  }
  Total->CanopyWater /= NPixels;

  /* temperature data */
  Total->MeanTair /= NPixels;
  
  /* average radiation data */
  for (i = 0; i < Veg->MaxLayers + 1; i++) {
    Total->Rad.NetShort[i] /= NPixels;
    Total->Rad.LongIn[i] /= NPixels;
    Total->Rad.LongOut[i] /= NPixels;
  }
  Total->Rad.PixelNetShort /= NPixels;
  Total->Rad.PixelLongIn /= NPixels;
  Total->Rad.PixelLongOut /= NPixels;

  Total->RadClass.Beam /= NPixels;
  Total->RadClass.Diffuse /= NPixels;
  
  /* average snow data */
  Total->Snow.Swq /= NPixels;
  Total->Snow.Melt /= NPixels;
  Total->Snow.PackWater /= NPixels;
  Total->Snow.TPack /= NPixels;
  Total->Snow.SurfWater /= NPixels;
  Total->Snow.TSurf /= NPixels;
  Total->Snow.ColdContent /= NPixels;
  Total->Snow.Albedo /= NPixels;
  Total->Snow.Depth /= NPixels;
  Total->Snow.VaporMassFlux /= NPixels;
  Total->Snow.CanopyVaporMassFlux /= NPixels;  

  /* average soil moisture data */
  Total->Soil.Depth /= NPixels;
  for (i = 0; i < Soil->MaxLayers; i++) {
    Total->Soil.Moist[i] /= NPixels;
    Total->Soil.Perc[i] /= NPixels;
    Total->Soil.Temp[i] /= NPixels;
  }
  Total->Soil.Moist[Soil->MaxLayers] /= NPixels;
  Total->Soil.TableDepth /= NPixels;
  Total->Soil.WaterLevel /= NPixels;
  Total->Soil.SatFlow /= NPixels;
  Total->Soil.TSurf /= NPixels;
  Total->Soil.Qnet /= NPixels;
  Total->Soil.Qs /= NPixels;
  Total->Soil.Qe /= NPixels;
  Total->Soil.Qg /= NPixels;
  Total->Soil.Qst /= NPixels;
  Total->SoilWater /= NPixels;
  Total->Runoff /= NPixels;
  Total->ChannelInt /= NPixels;
  Total->RoadInt /= NPixels;
  Total->CulvertReturnFlow /= NPixels;
  Total->CulvertToChannel /= NPixels;
  Total->RunoffToChannel /= NPixels;
  Total->Soil.GwRecharge /= NPixels;
  Total->GwDepth /= NPixels;
  Total->frozenSoilHeatFlux /= NPixels; 
  Total->Soil.QuickVol /= NPixels;
  Total->GwToChannel /= TotalArea;
  Total->channelLoss /= TotalArea;
  Total->GwDeepLoss /= TotalArea;
}
