#PCRGLOB-WB: run oldcalc -f pcrglobwb_v1.1.txt 365 10
#Rens van Beek, Dept. Physical Geography, Utrecht University, 9/11/2005
#Basic water balance model, thus omitting routing and water temperature
#-with 2 surface covers (tall and short vegetation) and 3 layers of soil compartment
#-includes surface runoff and soil evaporation based on the improved Arno scheme
# of Hagemann and Gates (2003)
#-climate input based on CRU data downscaled with ERA-40 daily surface fields
#-change made to script regarding fractioning of soil evaporation and transpiration
# and reports on fraction of EACT over EPOT included on 20/08/2008
# NOTE: WS... for reducing bare soil evaporation is obsolete as a result of changes
# Updates made to root fractions (evapotranspiration parameters dependent on relative
# root fractions RFW rather than on absolute fractions RFRAC) and CR2 limited to parameterized
# storage capacity of second layer rather than on porosity in root zone (17/11/2008)
# Capillary rise is restricted to the fractional area influenced by the groundwater, CRFRAC
# 20/04/2011: Soil moisture scheme has been changed to allow for full saturation of the soil layers
# by including percolation over the contact between the second to third layer to the infiltration;
# both percolation and capillary rise over the first and second layer has been made dependent on
# the square root of k(theta_eff) to include the effect of less pervious layers;
# to limit recursive drying and wetting of the first layer by percolation to and capillary rise from
# the second, drainage from the first layer is restricted to water in excess of field capacity
# under wet conditions (water in excess of field capacity) and capillary rise from the second to
# first layer is based on the gradient rather than on the difference in relative degree of saturation
# alone and is capped to maximum field capacity; see #* for comments on the changes
# Modification 28/06/2012: minor issue with precipitation partitioning following interception solved
# Modification 17/12/2012: water balance calculation modified to return monthly totals

# 3 Aug 2015: Edwin modified the following (for the simplification of debugging purpose)
# - FRACWAT = 0.0 and KC_WATSTACK = 0.0 (for the simplification of debugging process)
# - Reinstating B_ORO (so that it's consistent with Van Beek et al, 2011)
# - WMAX = SC1 + SC2 (instead of reading from file)
# - CR2_L[TYPE] = min(S3,CR2_L[TYPE]) (instead of CR2_L[TYPE]= min(VEGFRAC[TYPE]*S3,CR2_L[TYPE] as we don't have to reduce it with VEGFRAC) 
# - Fixing BCH_50 (now it uses BCH_50[TYPE] instead of BCH_50)
# - Constrain that 0.0 <= SATFRAC_L =<= 1.0
# - Evaporation from intercepted water <= T_p[TYPE]



binding

# Maps & TSS - input
######################################################################################################
######################################################################################################


# General
######################################################################################################

  Duration = scalar(1);               # timestep in days
  LANDMASK = maps\catclone.map;       # clone map representing landmask of earth surface
  CELLAREA = maps\cellarea30.map;     # surface (m2) of cell covered by total land surface


# Surface water
######################################################################################################

  LDD = maps\lddsound_30min.map;      # local drainage direction map

# FRACWAT  = maps\glwd130m_fracw.map; # fraction of cell area covered by fresh water
  FRACWAT  = scalar(0);               # set to zero for debugging purpose

# KC_WATSTACK = maps\kc_wat;		  # composite crop factor for channels and wetlands or lakes
  KC_WATSTACK = scalar(0.0);          


# Groundwater
######################################################################################################

  KQ3            = maps\globalalpha.map;            # recession coefficient for store 3 (day-1): drainage
  SPECYIELD3     = maps\specificyield.map;          # specific yield for aquifer
  DZS3INFLUENCED = scalar(5.0);                     # additional zone influenced by groundwater

  DZREL0001      = maps\hydro1k_dzrel0001.map;		# maps of relative elevation above floodplain, in percent
  DZREL0005      = maps\hydro1k_dzrel0005.map;
  DZREL0010      = maps\hydro1k_dzrel0010.map;
  DZREL0020      = maps\hydro1k_dzrel0020.map;
  DZREL0030      = maps\hydro1k_dzrel0030.map;
  DZREL0040      = maps\hydro1k_dzrel0040.map;
  DZREL0050      = maps\hydro1k_dzrel0050.map;
  DZREL0060      = maps\hydro1k_dzrel0060.map;
  DZREL0070      = maps\hydro1k_dzrel0070.map;
  DZREL0080      = maps\hydro1k_dzrel0080.map;
  DZREL0090      = maps\hydro1k_dzrel0090.map;
  DZREL0100      = maps\hydro1k_dzrel0100.map;


# Meteo
######################################################################################################
# Climatic input - values in m and degC per day

 PRPSTACK        = maps\meteo\pr;                   # map stack of daily precipitation rates (m/day)
 TASTACK         = maps\meteo\ta;                   # map stack of average daily temperature (degC)
 EVAPSTACK       = maps\meteo\e0p;                  # map stack of daily reference potential evapotranspiration rates (m/day)



# Vegetation parameters per cover type
######################################################################################################

 COVERTYPE = [
   SHORT = sv,
   TALL  = tv];							                       # array of cover type: 1) short, 2) tall
 COVERTABLE = maps\param_permafrost_edwinVersion.tbl;	       # table with parameterization per cover type

 VEGFRAC[COVERTYPE] = index(COVERTABLE);       # subdivision in cover type
                                                                   
 CF_SHORTSTACK   = maps\cover_fraction\cv_s;   # fractional vegetation cover (-) per vegetation type
 CF_TALLSTACK    = maps\cover_fraction\cv_t;			

 SMAX_SHORTSTACK = maps\interception_capacity_input\smax_s;    # interception storage (m) per vegetation type
 SMAX_TALLSTACK  = maps\interception_capacity_input\smax_t;

 KC_SHORTSTACK = maps\crop_coefficient\kc_s;   # crop factor (-) per vegetation type
 KC_TALLSTACK  = maps\crop_coefficient\kc_t;


# Minimum bare soil crop factor (-)     
######################################################################################################

  KCMIN = scalar(0.2);                         # used to separate soil evaporation from transpiration


# Snow routine parameters: constants    
######################################################################################################

  TT    = scalar(0.0);     # threshold temperature for freezing/thawing (degC)
  CFMAX = scalar(0.0025);  # degree-day factor (m??C-1?d-1) was: 0.0055
  SFCF  = scalar(1.00);    # snowfall correction factor (-)
  CWH   = scalar(0.10);    # water holding capacity snow cover (-)
  CFR   = scalar(0.05);    # refreezing coefficient (-)


# Topographical parameters 
######################################################################################################

 LSLOPE   = maps\globalbcat.map;               # slope length (m)
 TANSLOPE = maps\globalgradslope.map;          # gradient of slope (m/m)
 B_ORO    =	maps\globalboro.map;               # shape coefficient related to orography


# Soil parameters per cover type and layer (i)
######################################################################################################
 
 THETASAT1[COVERTYPE] = index(COVERTABLE);          # saturated volumetric moisture content (m3.m-3)
 THETASAT2[COVERTYPE] = index(COVERTABLE);          # first and second layer
 THETARES1[COVERTYPE] = index(COVERTABLE);          # residual volumetric moisture content (m3.m-3)
 THETARES2[COVERTYPE] = index(COVERTABLE);          # first and second layer
 KS1[COVERTYPE]       = index(COVERTABLE);          # saturated hydraulic conductivity (m.day-1)
 KS2[COVERTYPE]       = index(COVERTABLE);          # first and second layer
 PSI_A1[COVERTYPE]    = index(COVERTABLE);          # air entry value (m) according to SWRC of Clapp & Hornberger (1978)
 PSI_A2[COVERTYPE]    = index(COVERTABLE);          # first and second layer
 BCH1[COVERTYPE]      = index(COVERTABLE);          # Pore size distribution parameter according to Clapp and Hornberger (1978)
 BCH2[COVERTYPE]      = index(COVERTABLE);          # first and second layer
 Z1[COVERTYPE]        = index(COVERTABLE);          # depth of first store
 Z2[COVERTYPE]        = index(COVERTABLE);          # depth of second store
 RFRAC1[COVERTYPE]    = index(COVERTABLE);		    # root fraction per soil layer
 RFRAC2[COVERTYPE]    = index(COVERTABLE);		    # first and second layer
                                                    
#SC1[COVERTYPE]       = index(COVERTABLE);          # total storage per layer and for the entire soil profile
#SC2[COVERTYPE]       = index(COVERTABLE);          
#WMAX[COVERTYPE]      = index(COVERTABLE);

 MINFRAC[COVERTYPE]   = index(COVERTABLE);		    # ratio of min soil depth over average soil depth
 MAXFRAC[COVERTYPE]   = index(COVERTABLE);		    # ratio of max soil depth over average soil depth
 P2_IMP[COVERTYPE]    = index(COVERTABLE);		    # fractional area where percolation to groundwater store is impeded


# overall
######################################################################################################
  PSI_FC  = scalar(1.00);                           # matric suction at field capacity (m)
  PSI_50  = scalar(3.33);                           # matric suction at which transpiration is halved (m)
  BCH_ADD = scalar(3);                              # addition for kr-relationship of Clapp & Hornberger (1978; default 3)


# STATIONS        = maps\grdc_stations.map;		    # ID of gauging stations
# MONITORSTATIONS = maps\hol_climatelocs.map;		# locations to monitor soil behaviour, one per climate zone  #**


# Water consumption: abstractions in m3 per day from the open water and groundwater store
  POTABSTR_S0 = scalar(0);                          # from open water surface (m3/d)
  POTABSTR_S3 = scalar(0);                          # from groundwater store (m3/d)


# Initial storages for local variables, cell averages reported by default
##########################################################################################################################
  SC_INI[COVERTYPE]   = index(COVERTABLE);                            # initial snow cover (m)
  SCF_INI[COVERTYPE]  = index(COVERTABLE);                            # initial liquid water stored in snow cover  (m)
  INTS_INI[COVERTYPE] = index(COVERTABLE);                            # initial interception storage (m)
  S1_INI[COVERTYPE]   = index(COVERTABLE);                            # initial storage in upper store (m)
  S2_INI[COVERTYPE]   = index(COVERTABLE);                            # initial storage in second store (m)
  Q2_INI[COVERTYPE]   = index(COVERTABLE);                            # initial drainage from second store (m)
  S3_INI              = maps\initials\storGroundwater_initial.map;    # initial storage in lower store (m)



# Maps & TSS - output
##########################################################################################################################

 # forcing variables (re-written) 

 prp_input = oldcalc_results\pr;                # precipitation (m/day)  
 ta_input  = oldcalc_results\ta;                # temperature (degree C)
 e0p_input = oldcalc_results\e0p;               # reference potential evaporation (m/day) 


 SC       = oldcalc_results\snowcov;			# snow cover (m)
 SC_L     = oldcalc_results\sc;				    # idem, local value
 SCF      = oldcalc_results\snowliq;			# liquid water stored in snow cover (m)
 SCF_L    = oldcalc_results\scf;				# idem, local value
 SCTSS    = oldcalc_results\snowtot.tss;		# total snow storage in equivalent water height (m)
                                                 
 INTS     = oldcalc_results\intstor;			# interception storage (m)
 INTS_L   = oldcalc_results\ints;			    # idem, local value
 INTSTSS  = oldcalc_results\ints.tss;		    # as above, timeseries
          
 EACT     = oldcalc_results\eact;				# actual evapotranspiration (m)
 EACTTSS  = oldcalc_results\eact.tss;			# as above, timeseries
 ETPOT    = oldcalc_results\etpot;				# total potential evapotranspiration (m)
 EFRAC    = oldcalc_results\efrac;				# fraction of actual over potential evapotranspiration (-)
 ESPOT    = oldcalc_results\espot;				# potential soil evaporation (m)
 ESACT    = oldcalc_results\esact;				# actual soil evaporation (m)
 EWAT     = oldcalc_results\ewat;				# potential evaporation over open water
          
 T1POT    = oldcalc_results\t1pot;              # potential transpiration drawn from first soil layer (m)
 T1ACT    = oldcalc_results\tactUpp;            # idem, actual
          
 T2POT    = oldcalc_results\t2pot;              # potential transpiration drawn from second soil layer (m)
 T2ACT    = oldcalc_results\tactLow;            # idem, actual
          
 T_ACT    = oldcalc_results\tact;               # total transpiration from the entire layers

 INT_EVAP = oldcalc_results\int_evap;           # evaporation from intercepted water  
 SCF_EVAP = oldcalc_results\scf_evap;           # evaporation from snow free water

 SATFRAC = oldcalc_results\satf;				# fraction saturated area (-)
 WACT    = oldcalc_results\wact;				# actual water storage within root zone (m)

 CRFRAC  = oldcalc_results\crfrac;			    # fraction of soil surface influenced by capillary rise

 P0     = oldcalc_results\infl;                 # infiltration

 P1     = oldcalc_results\p1x;					# percolation from first layer
 P2     = oldcalc_results\p2x;                  # percolation from second layer
 CR1    = oldcalc_results\cr1x;					# capillary rise from first layer
 CR2    = oldcalc_results\cr2x;                 # capilary rise from second layer

 Q1     = oldcalc_results\qDr;					# direct runoff (m)
 Q1TSS  = oldcalc_results\q1.tss;				# idem, timeseries
 Q1S    = oldcalc_results\q1s;					# direct runoff attributable to snow melt (m)
 Q1STSS = oldcalc_results\q1snow.tss;		    # idem, timeseries

 Q2     = oldcalc_results\qSf;					# runoff from second store (flow) (m)
 Q2_L   = oldcalc_results\q2;					# idem, local value

 Q2TSS  = oldcalc_results\q2.tss;				# idem, timeseries

 Q3     = oldcalc_results\qBf;					# runoff from lower store (base flow) (m)
 Q3TSS  = oldcalc_results\q3.tss;				# idem, timeseries

 QLOC   = oldcalc_results\qLoc;				    # specific runoff (m)

 S1     = oldcalc_results\sUpp; 				# storage in upper store (m)
 S1_L   = oldcalc_results\s1;					# idem, local value
 S1TSS  = oldcalc_results\stor1.tss;			# as above, timeseries

 S2     = oldcalc_results\sLow; 				# storage in second store (m)
 S2_L   = oldcalc_results\s2;					# idem, local value
 S2TSS  = oldcalc_results\stor2.tss;			# as above, timeseries

 R3     = oldcalc_results\rch;					# recharge to third store
 R3AVG  = oldcalc_results\r3_avg.map;			# average recharge to the third, groundwater store (m)

 S3     = oldcalc_results\sGw;			        # storage in lower store (m)
 S3TSS  = oldcalc_results\stor3.tss;	        # as above, timeseries

 QW     = oldcalc_results\qw;					# change in storage of freshwater surface

 
 # channel discharge has been retained as check but is not used
 QAVG        = oldcalc_results\qaverage.map;	# channel discharge - annual average (m/s)
 QCHANNEL    = oldcalc_results\qc;				# channel discharge - daily value
 QCHANNELTSS = oldcalc_results\qchannel.tss;	# timeseries of total discharge (m/s)


 # water balance components and budget check included here

 # - monthly water balance fluxes  in m waterslice per unit surface area of pertinent type (water, land) per month

 PRECMON  = oldcalc_results\precmon;	# total precipitation over unit area
 ETPMON   = oldcalc_results\etpmon;	    # total potential evapotranspiration over unit land area (crop specific evapotranspiration)

 QWMON    = oldcalc_results\qwmon;		# change in storage of freshwater surface over water surface
 Q1MON    = oldcalc_results\q1mon;		# runoff over land surface from first, top store
 Q2MON    = oldcalc_results\q2mon;		# runoff over land surface from second store
 Q3MON    = oldcalc_results\q3mon;		# runoff over land surface from third, groundwater store
 EWATMON  = oldcalc_results\ewatmon;	# total potential evaporation over water surface
 EACTMON  = oldcalc_results\eactmon;	# total actual evapotranspiration over land surface
 ESACTMON = oldcalc_results\esactmon;	# actual evaporation from the soil surface
 TACTMON  = oldcalc_results\tactmon;	# actual transpiration from the soil


 # - budget check
 
 PTOTTSS   = oldcalc_results\ptot.tss;			   # total cumulative rainfall and initial accumulated storage (km3)
 ETOT      = oldcalc_results\etot.map;			   # total cumulative actual evapotranspiration (m)
 ETOTTSS   = oldcalc_results\etot.tss;			   # idem, as timeseries (km3)
 QTOTTSS   = oldcalc_results\qtot.tss;			   # total cumulative discharge (km3)
 STOTTSS   = oldcalc_results\stot.tss;			   # total accumulated storage (km3)
 STOT_ACT  = oldcalc_results\stot.map;		       # total active storage (km3)
 STOT_AVG  = oldcalc_results\stotavg.map;		   # average soil storage (m)
 MBE       = oldcalc_results\mbe.map;			   # absolute local mass balance error (m)
 MBR       = oldcalc_results\mbcheck.map;		   # relative, local mass balance error (-)
 MBRTSS    = oldcalc_results\mbcheck.tss;		   # idem, as time series at (selected) outlets
 MBRMONTSS = oldcalc_results\mbcheck_selected.tss; # idem, as time series at (selected) monitoring stations #**

areamap
 LANDMASK;

timer
 1 $12 1;
 rep1 = endtime;					                # times to report TSS and initial maps
 rep2 = 1+1..endtime;				                # times to report soil states and fluxes (default: daily)
#rep2 = $15+$15..endtime;				            # times to report soil states and fluxes
 rep3 = $1,$2,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12;	    # times to report monthly updates of water balance components

initial

#Initialization, dependent on COVERTYPE
#Initial storages
 SC= scalar(0);
 SCF= scalar(0);
 INTS= scalar(0);
 S1= scalar(0);
 S2= scalar(0);
 S3= S3_INI;

 foreach TYPE in COVERTYPE{

   #Initial storages
   #-states and fluxes per vegetation type
   SC_L[TYPE]= SC_INI[TYPE];			#initial snow cover (m)
   SCF_L[TYPE]= SCF_INI[TYPE];			#initial liquid water stored in snow cover  (m)
   INTS_L[TYPE]= INTS_INI[TYPE];		#initial interception storage (m)
   S1_L[TYPE]= S1_INI[TYPE];			#initial storage in upper store (m)
   S2_L[TYPE]= S2_INI[TYPE];			#initial storage in second store (m)
   Q2_L[TYPE]= Q2_INI[TYPE];			#initial drainage from second store (m)

   #-total storages
   SC= SC+VEGFRAC[TYPE]*SC_L[TYPE];
   SCF= SCF+VEGFRAC[TYPE]*SCF_L[TYPE];

   INTS= INTS+VEGFRAC[TYPE]*INTS_L[TYPE];

   S1= S1+VEGFRAC[TYPE]*S1_L[TYPE];
   S2= S2+VEGFRAC[TYPE]*S2_L[TYPE];

   #-soil parameters

   BCB1[TYPE]= 2*BCH1[TYPE]+BCH_ADD;		#Campbell's (1974) coefficient to calculate the relative,
   BCB2[TYPE]= 2*BCH2[TYPE]+BCH_ADD;		#unsaturated hydraulic conductivity

						#degree of saturation at field capacity for first and second layer
						#and corresponding unstaturated hydraulic conductivity
   THEFF1_FC[TYPE]= (PSI_FC/PSI_A1[TYPE])**(-1/BCH1[TYPE]);
   KTHEFF1_FC[TYPE]= max(0,THEFF1_FC[TYPE]**BCB1[TYPE]*KS1[TYPE]);
   THEFF2_FC[TYPE]= (PSI_FC/PSI_A2[TYPE])**(-1/BCH2[TYPE]); 
   KTHEFF2_FC[TYPE]= max(0,THEFF2_FC[TYPE]**BCB2[TYPE]*KS2[TYPE]);

   						#centroid-lag for constant drainable pore space (day),
						#rewritten to fractions given the duration of one timestep
   TCL[TYPE]= Duration*(2*KS2[TYPE]*TANSLOPE)/(LSLOPE*(1-THEFF2_FC[TYPE])*(THETASAT2[TYPE]-THETARES2[TYPE]));
   TCL[TYPE]= max(0.0, TCL[TYPE]);

   #-storage parameters related to Improved Arno Scheme

   SC1[TYPE] = Z1[TYPE] * (THETASAT1[TYPE] - THETARES1[TYPE]);
   SC2[TYPE] = Z2[TYPE] * (THETASAT2[TYPE] - THETARES2[TYPE]);

   # storage in top soil limiting bare soil evapotranspiration
   WSMAX[TYPE] = SC1[TYPE];
   WMAX[TYPE]  = SC1[TYPE] + SC2[TYPE];
   
   # minimum storage capacity set to fraction of mean total
   WMIN[TYPE]  = MINFRAC[TYPE]*WMAX[TYPE];
   WSMIN[TYPE] = MINFRAC[TYPE]*WSMAX[TYPE];
   
   # range in storage capacity
   WRANGE[TYPE]= WMAX[TYPE]-WMIN[TYPE];
   WSRANGE[TYPE]= WSMAX[TYPE]-WSMIN[TYPE];
   
   # shape coefficient for distribution of maximum storage capacity, maximum set to fraction of total
   BCF[TYPE] = max(0.001,(MAXFRAC[TYPE]-1)/(1-MINFRAC[TYPE])+B_ORO-0.01);

   # weighed root fractions
   RFW1[TYPE] = if(RFRAC1[TYPE]+RFRAC2[TYPE] > 0,
		min(1.0,RFRAC1[TYPE]/(RFRAC1[TYPE]+RFRAC2[TYPE])),0.0);
   RFW2[TYPE] = if(RFRAC1[TYPE]+RFRAC2[TYPE] > 0.0,
		min(1.0,RFRAC2[TYPE]/(RFRAC1[TYPE]+RFRAC2[TYPE])),0.0);

   # average degree of saturation at which actual transpiration is halved

   THEFF_50[TYPE] = if( (SC1[TYPE]*RFW1[TYPE]+SC2[TYPE]*RFW2[TYPE]) > 0,  
     (SC1[TYPE]*RFW1[TYPE]*(PSI_50/PSI_A1[TYPE])**(-1/BCH1[TYPE])+
     SC2[TYPE]*RFW2[TYPE]*(PSI_50/PSI_A2[TYPE])**(-1/BCH2[TYPE]))/
     (SC1[TYPE]*RFW1[TYPE]+SC2[TYPE]*RFW2[TYPE]),0.5);
   
   BCH_50[TYPE] = if( (SC1[TYPE]*RFW1[TYPE]+SC2[TYPE]*RFW2[TYPE]) > 0,(SC1[TYPE]*RFW1[TYPE]*BCH1[TYPE]+SC2[TYPE]*RFW2[TYPE]*BCH2[TYPE])/
    (SC1[TYPE]*RFW1[TYPE]+SC2[TYPE]*RFW2[TYPE]),0.5*(BCH1[TYPE]+BCH2[TYPE]));

}

# cumulative of monthly water balance
  PRECMON  = scalar(0);
  ETPMON   = scalar(0);
  QWMON    = scalar(0);
  Q1MON    = scalar(0);
  Q2MON    = scalar(0);
  Q3MON    = scalar(0);
  EWATMON  = scalar(0);
  EACTMON  = scalar(0);
  ESACTMON = scalar(0);
  TACTMON  = scalar(0);

# cumulative fluxes for average discharge and budget check
  QCUM     = scalar(0);
  R3CUM    = scalar(0);
  STOT_CUM = scalar(0);
  SLOCINI  = (1-FRACWAT)*(S1+S2+S3+INTS+SC+SCF);
  PTOT     = scalar(0);
  ETOT     = scalar(0);
  QTOT     = scalar(0);

# check on meteo generated
 PTEST= scalar(0);
 TTEST= scalar(0);
 ETEST= scalar(0);

dynamic

#----------------
# Meteo
#----------------
#Meteorological input as total/average per time step
# TA		(?C):	average temperature
# EVAP		(m):	reference potential evapotranspiration
# EWAT		(m):	as above, imposed on water surface
# PRP		(m):	liquid precipitation
# SNOW		(m):	snow in water equivalent
# PRPTOT	(m):	total precipitation


# precipitation

 PRPTOT = max(0.,timeinput(PRPSTACK))*Duration*timeslice();

report (rep2) prp_input = PRPTOT;                 # precipitation (m/day)  

# temperature

 TA= timeinput(TASTACK);

report (rep2) ta_input  = TA;                     # temperature (degree C)

# partitioning rain and snow and converting evapotransiration

 SNOW = if(TA<TT,PRPTOT,0);                        

 PRP    = PRPTOT-SNOW;                            # not needed
 SNOW   = SFCF*SNOW;
 PRPTOT = PRP+SNOW;

# reference potential evaporation

 EVAP = max(0.,timeinputsparse(EVAPSTACK))*Duration*timeslice(); 

report (rep2) e0p_input = EVAP;                   # reference potential evaporation (m/day) 


#report (rep2) EWAT = timeinputsparse(KC_WATSTACK)*EVAP;
#report (rep2) EWAT = KC_WATSTACK*EVAP;
               EWAT = KC_WATSTACK*EVAP;
#--------------------
# Fresh water surface
#--------------------
#QW		(m):	local change in storage of fresh water surface (can be negative)
#report (rep2) QW = if(LANDMASK,PRPTOT-(EWAT+POTABSTR_S0));
               QW = if(LANDMASK,PRPTOT-(EWAT+POTABSTR_S0));


#----------------
# Land surface
#----------------

#Current approximate height of groundwater table and corresponding reach of cell under
# influence of capillary rise
DZS3= S3/SPECYIELD3+DZS3INFLUENCED;
CRFRAC= min(1.0,1.0-(DZREL0100-DZS3)*0.1/max(1e-3,DZREL0100-DZREL0090));
CRFRAC= if(DZS3<DZREL0090,0.9-(DZREL0090-DZS3)*0.1/max(1e-3,DZREL0090-DZREL0080),CRFRAC);
CRFRAC= if(DZS3<DZREL0080,0.8-(DZREL0080-DZS3)*0.1/max(1e-3,DZREL0080-DZREL0070),CRFRAC);
CRFRAC= if(DZS3<DZREL0070,0.7-(DZREL0070-DZS3)*0.1/max(1e-3,DZREL0070-DZREL0060),CRFRAC);
CRFRAC= if(DZS3<DZREL0060,0.6-(DZREL0060-DZS3)*0.1/max(1e-3,DZREL0060-DZREL0050),CRFRAC);
CRFRAC= if(DZS3<DZREL0050,0.5-(DZREL0050-DZS3)*0.1/max(1e-3,DZREL0050-DZREL0040),CRFRAC);
CRFRAC= if(DZS3<DZREL0040,0.4-(DZREL0040-DZS3)*0.1/max(1e-3,DZREL0040-DZREL0030),CRFRAC);
CRFRAC= if(DZS3<DZREL0030,0.3-(DZREL0030-DZS3)*0.1/max(1e-3,DZREL0030-DZREL0020),CRFRAC);
CRFRAC= if(DZS3<DZREL0020,0.2-(DZREL0020-DZS3)*0.1/max(1e-3,DZREL0020-DZREL0010),CRFRAC);
CRFRAC= if(DZS3<DZREL0010,0.1-(DZREL0010-DZS3)*0.05/max(1e-3,DZREL0010-DZREL0005),CRFRAC);
CRFRAC= if(DZS3<DZREL0005,0.05-(DZREL0005-DZS3)*0.04/max(1e-3,DZREL0005-DZREL0001),CRFRAC);
CRFRAC= if(DZS3<DZREL0001,0.01-(DZREL0001-DZS3)*0.01/max(1e-3,DZREL0001),CRFRAC);
CRFRAC= if(FRACWAT<1,max(0,CRFRAC-FRACWAT)/(1-FRACWAT),0);

CRFRAC=max(0.0, CRFRAC);
CRFRAC=min(1.0, CRFRAC);

#Water balance of first and second store,
#split in two horizontal compartments for (1) short and (2) tall vegetation (COVERTYPE)
#to simulate effect of rooting depth over two vertical layers with maximum storage SC1 and SC2
#-initialization of states and fluxes for the present time step, local variables identified by '_L'
 TYPECOUNT= scalar(0);
 SC= scalar(0);
 SCF= scalar(0);
 INTS= scalar(0);
 ETPOT= scalar(0);
 ESPOT= scalar(0);
 ESACT= scalar(0);
 T1POT= scalar(0);
 T1ACT= scalar(0);
 T2POT= scalar(0);
 T2ACT= scalar(0);
 EACT= scalar(0);
 SATFRAC= scalar(0);
 WACT= scalar(0);
 S1= scalar(0);
 S2= scalar(0);
 P0= scalar(0);
 P1= scalar(0);
 P2= scalar(0);
 CR1= scalar(0);
 CR2= scalar(0);
 Es= scalar(0);
 Ta= scalar(0);
 Q1= scalar(0);
 Q1S= scalar(0);
 Q2= scalar(0);

 INT_EVAP = scalar(0);
 SCF_EVAP = scalar(0);

 foreach TYPE in COVERTYPE{
  #-Vegetation cover fraction

  #CFRAC	(-):	vegetation cover fraction
   TYPECOUNT= TYPECOUNT+1;
   CFRAC[TYPE]= if(TYPECOUNT== 1,timeinputsparse(CF_SHORTSTACK),timeinputsparse(CF_TALLSTACK));
  #-FAO crop factor for transpiration

  #KC		(-):	crop factor
   KC[TYPE]= if(TYPECOUNT== 1,timeinputsparse(KC_SHORTSTACK),timeinputsparse(KC_TALLSTACK));

  #-Potential bare soil evaporation and transpiration
  #ET_p		(m):	total potential evapotranspiration
  #ES_p		(m):	bare soil evaporation
  #T_p		(m):	transpiration
  #EACT		(-):	direct evapotranspiration

   ET_p[TYPE] = max(KCMIN, KC[TYPE]*EVAP);	
   ES_p[TYPE] = KCMIN*EVAP;
   T_p[TYPE]  = max(0.0, KC[TYPE]*EVAP-ES_p[TYPE]);

  #-Interception
  #INTCMAX	(m):	maximum interception storage capacity
  #ICC		(m):	equivalent interception storage capacity
  #INTS		(m):	interception storage
  #PRP		(m):	precipitation passing the canopy

   INTCMAX[TYPE] = if(TYPECOUNT== 1,timeinputsparse(SMAX_SHORTSTACK),timeinputsparse(SMAX_TALLSTACK));
   ICC[TYPE]     = CFRAC[TYPE]*INTCMAX[TYPE];
   PRP           = (1-CFRAC[TYPE])*PRPTOT+max(CFRAC[TYPE]*PRPTOT+INTS_L[TYPE]-ICC[TYPE],0);

   INTS_L[TYPE]  = max(0,INTS_L[TYPE]+PRPTOT-PRP);

   SNOW          = min(PRP,SNOW*min(1.0,if(PRPTOT>0,PRP/PRPTOT,0))); #**

   PRP           = max(0,PRP-SNOW); #**
   EACT_L[TYPE]  = min(INTS_L[TYPE],(T_p[TYPE]*if(ICC[TYPE]>0,INTS_L[TYPE]/ICC[TYPE],0)**(2/3)));

   EACT_L[TYPE]  = min(T_p[TYPE], EACT_L[TYPE]);
   
   INT_EVAP_L = EACT_L[TYPE];   # evaporation from intercepted water

   report (rep1) INTS_L[TYPE] = INTS_L[TYPE]-EACT_L[TYPE];
   T_p[TYPE]                  = max(0,T_p[TYPE]-EACT_L[TYPE]);

  #-Snow accumulation and melt
  #SC		(m):	snow cover
  #SCF		(m):	free water stored in snow cover
  #DSC		(m):	change in snow cover, - melt, + gain in snow or refreezing (CFR)
  #Pn		(m):	net liquid water transferred to soil
  #DSCR		(m):	relative contribution of snow melt to net liquid water transfer
  #ES_a		(m):	actual bare soil evaporation,
  #			here used to subtract any evaporation from liquid phase of snow cover
   
   DSC[TYPE] = if(TA<=TT,CFR*SCF_L[TYPE],-min(SC_L[TYPE],max(TA-TT,0)*CFMAX*Duration*timeslice()));
   SC_L[TYPE] = max(0.0, SC_L[TYPE] + SNOW + DSC[TYPE]);

   SCF_L[TYPE] = SCF_L[TYPE] - DSC[TYPE] + PRP;

   Pn = max(0,SCF_L[TYPE]- CWH*SC_L[TYPE]);
   
   DSCR = if(Pn>0,max(-DSC[TYPE],0)/Pn,0);
   
   SCF_L[TYPE]= max(0,SCF_L[TYPE]-Pn);
   ES_a[TYPE]= min(SCF_L[TYPE],ES_p[TYPE]);
   report (rep1) SCF_L[TYPE]= SCF_L[TYPE]-ES_a[TYPE];
   ES_p[TYPE]= max(0,ES_p[TYPE]-ES_a[TYPE]);
   EACT_L[TYPE]= EACT_L[TYPE]+ES_a[TYPE];

   SCF_EVAP_L = ES_a[TYPE];   # evaporation from snow free water

  #-Direct runoff and infiltration based on improved Arno scheme
  # partial runoff when not entirely saturated (condition 1), else complete saturation excess
  #BCF		(-):	b coefficient of soil water storage capacity distribution
  #WMIN, WMAX	(m):	root zone water storage capacity, area-averaged values
  #W		(m):	actual water storage in root zone
  #WRANGE, DW,
  #WFRAC	(m):	computation steps to ease runoff calculation, resp.
  #			(WMAX-WMIN), (WMAX-W) and DW/WRANGE
  #			WFRAC capped at 1
  #WFRACB	(nd):	DW/WRANGE raised to the power (1/(b+1))
  #SATFRAC	(-):	fractional saturated area
  #WACT		(m):	actual water storage within rootzone
  #THEFF(i)	(-):	effective degree of saturation
  #
  #-Saturated and unsaturated hydraulic conductivity, matric suction and gradient
  #*PSI(i)	(m):	matric suction in layer(i)
  #*GRAD		(-):	gradient for capillary rise
  #KS(i)	(m/d):	saturated hydraulic conductivity
  #KTHEFF(i)	(m/d):	unsaturated hydraulic conductivity
  #*KTHVERT	(m/d):	idem, exchange between layers, capped at field capacity
  #BCH(i)	(-):	pore size distribution factor of Clapp and Hornberger (1978)
  #BCB(i)	(-):	Campbell's (1974) coefficient to calculate the relative,
  #			unsaturated hydraulic conductivity
  #Pn		(m):	net liquid precipitation, reduced if WMIN not exceeded
  #Q1		(m):	direct or surface runoff
  #Q1S		(-):	direct or surface runoff directly attributable to snow melt
  #P0		(m):	infiltration
  #P(i)		(m):	percolation from layer(i) to layer(i+1)
  #*
  #* Note that here the scaling of the unsaturated hydraulic conductivity by that at field capacity
  #* has been placed which earlier was included in the section where the soil hydrological fluxes (percolation)
  #* were calculated. In order to allow for full saturation, the percolation from the second to the third
  #* layer is already calculated here and added to the infiltration.
  #* To ensure that scaling of k(thetaeff) by that value at field capacity for percolation and capillary rise of the first
  #* layer does not affect the other results, the additional variable k(thetaeff) in the vertical direction over 
  #* the two soil layers is introduced. The capillary rise into this layer is no longer dependent on moisture but 
  #* primarily on the gradient. To this end, the new variables PSI(i), GRAD are introduced.
  #*
   THEFF1= max(0,S1_L[TYPE]/SC1[TYPE]);
   THEFF2= max(0,S2_L[TYPE]/SC2[TYPE]);

   THEFF1= min(1.0, THEFF1);
   THEFF2= min(1.0, THEFF2);

   PSI1= PSI_A1[TYPE]*max(0.01,THEFF1)**-BCH1[TYPE]; #*
   PSI2= PSI_A2[TYPE]*max(0.01,THEFF2)**-BCH2[TYPE]; #*

   GRAD= max(0,2*(PSI1-PSI2)/(Z1[TYPE]+Z2[TYPE])-1); #*

   KTHEFF1= max(0,THEFF1**BCB1[TYPE]*KS1[TYPE]);
   KTHEFF2= max(0,THEFF2**BCB2[TYPE]*KS2[TYPE]);

   KTHEFF1= min(KTHEFF1, KS1[TYPE]);
   KTHEFF2= min(KTHEFF2, KS2[TYPE]);

   KTHVERT= min(sqrt(KTHEFF1*KTHEFF2),
     (KTHEFF1*KTHEFF2*KTHEFF1_FC[TYPE]*KTHEFF2_FC[TYPE])**0.25); #*
 
   P2_L[TYPE]= min(KTHEFF2,sqrt(KTHEFF2*KTHEFF2_FC[TYPE]))*Duration*timeslice(); #*
 
   W[TYPE]= max(0,S1_L[TYPE]+S2_L[TYPE]);
   P0_L[TYPE]= Pn;
   Pn= W[TYPE]+Pn;
   Pn= Pn-max(WMIN[TYPE],W[TYPE]);
   W[TYPE]= if(Pn<0,WMIN[TYPE]+Pn,max(W[TYPE],WMIN[TYPE]));
   Pn= max(0,Pn);
   DW= max(0,WMAX[TYPE]-W[TYPE]);
   WFRAC= min(1,DW/WRANGE[TYPE]);
   WFRACB= WFRAC**(1/(1+BCF[TYPE]));

   SATFRAC_L[TYPE] = if(WFRACB>0,1-WFRACB**BCF[TYPE],1);
   
   SATFRAC_L[TYPE] = min(1.0, SATFRAC_L[TYPE]);
   SATFRAC_L[TYPE] = max(0.0, SATFRAC_L[TYPE]);
   
   WACT_L= (BCF[TYPE]+1)*WMAX[TYPE]-BCF[TYPE]*WMIN[TYPE]-(BCF[TYPE]+1)*WRANGE[TYPE]*WFRACB;
   Q1_L[TYPE]= max(0,Pn-(WMAX[TYPE]+P2_L[TYPE]-W[TYPE])+
     if(Pn>=(BCF[TYPE]+1)*WRANGE[TYPE]*WFRACB, 0,
     WRANGE[TYPE]*(WFRACB-Pn/((BCF[TYPE]+1)*WRANGE[TYPE]))**(BCF[TYPE]+1))); #*
   
   Q1_L[TYPE] = min(Q1_L[TYPE], P0_L[TYPE]);
   
   Q1S_L[TYPE]= min(1,DSCR)*Q1_L[TYPE];
   P0_L[TYPE]= P0_L[TYPE]-Q1_L[TYPE];

  #-Actual bare soil evaporation and transpiration based on the remainder of the potential
  # and limited to the available moisture content; top soil for ES, entire root zone for T
  #RFW(i)	(-):	root fraction per layer, corrected to 100%
  #WF(i)	(-):	weighing factor for fractioning transpiration per layer,
  #					based on total available moisture storage in soil, or else RFW
  #ES_p		(m):	potential bare soil evaporation
  #ES_a		(m):	actual bare soil evaporation
  #T_p(i)	(m):	potential transpiration per layer
  #T_a(i)	(m):	actual transpiration per layer
  #FRACTA	(-):	fraction of actual over potential transpiration

   WF1= if((S1_L[TYPE]+S2_L[TYPE])>0,RFW1[TYPE]*S1_L[TYPE]/
    max(1e-9,RFW1[TYPE]*S1_L[TYPE]+RFW2[TYPE]*S2_L[TYPE]),RFW1[TYPE]);
   WF2= if((S1_L[TYPE]+S2_L[TYPE])>0,RFW2[TYPE]*S2_L[TYPE]/
    max(1e-9,RFW1[TYPE]*S1_L[TYPE]+RFW2[TYPE]*S2_L[TYPE]),RFW2[TYPE]);

   FRACTA[TYPE]= (WMAX[TYPE]+BCF[TYPE]*WRANGE[TYPE]*(1-(1+BCF[TYPE])/BCF[TYPE]*WFRACB))/
     (WMAX[TYPE]+BCF[TYPE]*WRANGE[TYPE]*(1-WFRACB));
   FRACTA[TYPE]= (1-SATFRAC_L[TYPE])/(1+(max(0.01,FRACTA[TYPE])/THEFF_50[TYPE])**(-3*BCH_50[TYPE]));

   T1_p[TYPE]= WF1*T_p[TYPE];
   T2_p[TYPE]= max(0,T_p[TYPE]-T1_p[TYPE]);

   T_a1[TYPE]= FRACTA[TYPE]*T1_p[TYPE];
   T_a2[TYPE]= FRACTA[TYPE]*T2_p[TYPE];

   T_a[TYPE]= T_a1[TYPE]+T_a2[TYPE];

  #-actual bare soil evaporation
   ES_a[TYPE]= SATFRAC_L[TYPE]*min(ES_p[TYPE],KS1[TYPE]*Duration*timeslice())+
     (1-SATFRAC_L[TYPE])*min(ES_p[TYPE],KTHEFF1*Duration*timeslice());

  #-Percolation, subsurface storm flow and capillary rise
  #P(i)		(m):	percolation from layer(i) to layer(i+1)
  #CR(i)	(m):	capillary rise into layer(i) #*
  #Q2		(m):	lateral drainage from second store:
  #			dependent on net recharge to saturated wedge
  #			and centroid lag (Sloan and Moore, 1984)
  #			simplified by considering drainable pore space only
  #RQ2		(m):	recharge adding to or drawing from saturated wedge
  #-fluxes
  #*
  #* Note that here the sequence of calculation has been reversed. First, the direct runoff
  #* and the infiltration are modified by the amount of infiltration that is in excess of
  #* the infiltration capacity. Next, the percolation from the first into the second layer is
  #* computed as the flux due to the vertical unsaturated hydraulic conductivity over the two layers
  #* this flux can only drain the layer to field capacity when the latyer is wet; infiltration in
  #* excess to the storage capacity of layer 1 is passed on to the second layer. If the second layer
  #* becomes saturated, the excess water is passed back to the first layer by means of the percolation
  #* which may become negative in extreme cases (gain in storage in layer 1 due to return flow).
  #* This was already included in the model but the return flow not assigned to a particular flux.
  #* The capillary rise from the second to the first layer has been modified to consider the
  #* gradient and the vertical unsaturated hydraulic conductivity.  
  #*
   Q1_L[TYPE]= Q1_L[TYPE]+max(0,P0_L[TYPE]-KS1[TYPE]*Duration*timeslice());
   P0_L[TYPE]= min(P0_L[TYPE],KS1[TYPE]*Duration*timeslice());
   P1_L[TYPE]= KTHVERT*Duration*timeslice(); #*
   P1_L[TYPE]= if(THEFF1 > THEFF1_FC[TYPE],min(max(0,THEFF1-THEFF1_FC[TYPE])*SC1[TYPE],
     P1_L[TYPE]),P1_L[TYPE])+max(0,P0_L[TYPE]-(SC1[TYPE]-S1_L[TYPE])); #*
   CR1_L[TYPE]= min(max(0,THEFF1_FC[TYPE]-THEFF1)*SC1[TYPE],
     KTHVERT*GRAD*Duration*timeslice()); #*
   CR2_L[TYPE]= 0.5*(SATFRAC_L[TYPE]+CRFRAC)*min((1-THEFF2)*sqrt(KS2[TYPE]*KTHEFF2)*Duration*timeslice(),
     max(0,THEFF2_FC[TYPE]-THEFF2)*SC2[TYPE]);
   
   RQ2=  P2_IMP[TYPE]*(P1_L[TYPE]+CR2_L[TYPE]-(P2_L[TYPE]+CR1_L[TYPE]));
   Q2_L[TYPE]= max(TCL[TYPE]*RQ2+(1-TCL[TYPE])*Q2_L[TYPE],0);   
   #-water balance: scaled fluxes and new states
   # first layer 
   ADJUST= ES_a[TYPE]+T_a1[TYPE]+P1_L[TYPE];
   ADJUST= if(ADJUST>0,min(1,(max(0,S1_L[TYPE]+P0_L[TYPE]))/ADJUST),0);
   ES_a[TYPE]= ADJUST*ES_a[TYPE];
   T_a1[TYPE]= ADJUST*T_a1[TYPE];
   P1_L[TYPE]= ADJUST*P1_L[TYPE];
   # second layer
   ADJUST= T_a2[TYPE]+P2_L[TYPE]+Q2_L[TYPE];
   ADJUST= if(ADJUST>0,min(1,max(S2_L[TYPE]+P1_L[TYPE],0)/ADJUST),0);
   T_a2[TYPE]= ADJUST*T_a2[TYPE];
   P2_L[TYPE]= ADJUST*P2_L[TYPE];
   report (rep1) Q2_L[TYPE]= ADJUST*Q2_L[TYPE];

  #CR2_L[TYPE] = min(VEGFRAC[TYPE]*S3,CR2_L[TYPE]);
   CR2_L[TYPE] = min(S3,CR2_L[TYPE]);                                    # Edwin modified this line (as we don't have to reduce CR2_L with VEGFRAC)

   CR1_L[TYPE] = min(max(0,S2_L[TYPE]+P1_L[TYPE]-(T_a2[TYPE]+P2_L[TYPE]+Q2_L[TYPE])),CR1_L[TYPE]);
   S2_L[TYPE]  = max(0,S2_L[TYPE]+P1_L[TYPE]+CR2_L[TYPE]-
    (P2_L[TYPE]+Q2_L[TYPE]+CR1_L[TYPE]+T_a2[TYPE]));

   P1_L[TYPE] = P1_L[TYPE]-max(0,S2_L[TYPE]-SC2[TYPE]); #*
   S1_L[TYPE] = max(0,S1_L[TYPE]+P0_L[TYPE]+CR1_L[TYPE]-
    (P1_L[TYPE]+T_a1[TYPE]+ES_a[TYPE])); #*

   report (rep1) S2_L[TYPE]= min(S2_L[TYPE],SC2[TYPE]);

   Q1_L[TYPE] = Q1_L[TYPE]+max(0,S1_L[TYPE]-SC1[TYPE]);
   report (rep1) S1_L[TYPE]= min(S1_L[TYPE],SC1[TYPE]);

  #-total actual evapotranspiration
   EACT_L[TYPE]= EACT_L[TYPE]+ES_a[TYPE]+T_a1[TYPE]+T_a2[TYPE];

  #-adding local fluxes and states relative to vegetation fractions
  # as a function of RTN to correct for vegetation presence
   SC= SC+VEGFRAC[TYPE]*SC_L[TYPE];
   SCF= SCF+VEGFRAC[TYPE]*SCF_L[TYPE];
   INTS= INTS+VEGFRAC[TYPE]*INTS_L[TYPE];
   ETPOT= ETPOT+VEGFRAC[TYPE]*ET_p[TYPE];
   EACT= EACT+VEGFRAC[TYPE]*EACT_L[TYPE];
   ESPOT= ESPOT+VEGFRAC[TYPE]*ES_p[TYPE];
   ESACT= ESACT+VEGFRAC[TYPE]*ES_a[TYPE];
   T1POT= T1POT+VEGFRAC[TYPE]*WF1*T_p[TYPE];
   T1ACT= T1ACT+VEGFRAC[TYPE]*T_a1[TYPE];
   T2POT= T2POT+VEGFRAC[TYPE]*WF2*T_p[TYPE];
   T2ACT= T2ACT+VEGFRAC[TYPE]*T_a2[TYPE];
   SATFRAC= SATFRAC+VEGFRAC[TYPE]*SATFRAC_L[TYPE];
   WACT= WACT+VEGFRAC[TYPE]*WACT_L;
   S1= S1+VEGFRAC[TYPE]*S1_L[TYPE];
   S2= S2+VEGFRAC[TYPE]*S2_L[TYPE];
   P0= P0+VEGFRAC[TYPE]*P0_L[TYPE];
   P1= P1+VEGFRAC[TYPE]*P1_L[TYPE];
   P2= P2+VEGFRAC[TYPE]*P2_L[TYPE];
   CR1= CR1+VEGFRAC[TYPE]*CR1_L[TYPE];
   CR2= CR2+VEGFRAC[TYPE]*CR2_L[TYPE];
   Q1= Q1+VEGFRAC[TYPE]*Q1_L[TYPE];
   Q1S= Q1S+VEGFRAC[TYPE]*Q1S_L[TYPE];
   Q2= Q2+VEGFRAC[TYPE]*Q2_L[TYPE];

   INT_EVAP = INT_EVAP + VEGFRAC[TYPE]*INT_EVAP_L;
   SCF_EVAP = SCF_EVAP + VEGFRAC[TYPE]*SCF_EVAP_L;

   }

#---------------------------
#Overall states upper stores
#---------------------------
 report (rep2) INTS = INTS;
 report (rep2) P0   = P0;
 report (rep2) SC   = SC;
 report (rep2) SCF  = SCF;
 SATFRAC= SATFRAC;
#report (rep2)  WACT = WACT;
                WACT = WACT;
 report (rep2)  S1   = S1;
 report (rep2)  S2   = S2;

#---------------------------
#Overall fluxes upper stores
#---------------------------
 ETPOT= ETPOT;
 ESPOT= ESPOT;
 report (rep2) ESACT = ESACT;
 T1POT= T1POT;
 report (rep2) T1ACT = T1ACT;
 T2POT= T2POT;
 report (rep2) T2ACT = T2ACT;
 report (rep2) EACT  = EACT;
 report (rep2) T_ACT = T1ACT + T2ACT;
 
 report (rep2) INT_EVAP = INT_EVAP;
 report (rep2) SCF_EVAP = SCF_EVAP;
 
 EFRAC= if(ETPOT>0,EACT/ETPOT,1);
 P1= P1;
 report (rep2) Q1 = Q1;
 P2= P2;
 Q1S= Q1S;
 report (rep2) Q2= Q2;
 CR1= CR1;
 CR2= CR2;

#--------------------------
#Overall fluxes third store
#--------------------------
#Third reservoir
#R3		(m):	groundwater recharge
#R3AVG		(m):	average recharge
#S3		(m):	storage in third store, updated with current fluxes
#Q3		(m):	discharge from third reservoir, based on storage previous timestep
 report (rep2) R3 = P2-CR2;
 R3CUM= R3CUM+R3;
 report (rep1) R3AVG = R3CUM/time();
 ABSTR_S3= min(POTABSTR_S3,S3); 
 S3= max(0,S3+P2-(ABSTR_S3+CR2));
 report (rep2) Q3 = min(S3,KQ3*S3*Duration*timeslice());
 report (rep2) S3 = max(0,S3-Q3);

#-----------------------------
#Channel storage and discharge
#-----------------------------
#QLOC		(m):		local discharge
#QCHANNEL	(m3.s-1):	channel discharge
#QAVG		(m3.s-1):	channel discharge, averaged per year
 report (rep2) QLOC= Q1+Q2+Q3;
 QLOC = (1-FRACWAT)*QLOC+FRACWAT*QW;
#report (rep2) QCHANNEL= max(0,catchmenttotal(QLOC*CELLAREA,LDD)/(3600*24*Duration*timeslice()));
               QCHANNEL= max(0,catchmenttotal(QLOC*CELLAREA,LDD)/(3600*24*Duration*timeslice()));
 QCUM= QCUM+QCHANNEL;
 report (rep1) QAVG= QCUM/time();

#------------
# Water balance and budget check
#------------
#-update of monthly water balance components and budget check
 report (rep3) PRECMON= PRECMON+PRPTOT;
 report (rep3) ETPMON= ETPMON+EVAP;
 report (rep3) QWMON= QWMON+QW;
 report (rep3) Q1MON= Q1MON+Q1;
 report (rep3) Q2MON= Q2MON+Q2;
 report (rep3) Q3MON= Q3MON+Q3;
 report (rep3) EWATMON= EWATMON+EWAT;
 report (rep3) EACTMON=  EACTMON+EACT;
 report (rep3) ESACTMON= ESACTMON+ESACT;
 report (rep3) TACTMON= TACTMON+T1ACT+T2ACT;
 TIMERESET= boolean(0);
 TIMERESET= if(time() ==  $1,boolean(1),TIMERESET);
 TIMERESET= if(time() ==  $2,boolean(1),TIMERESET);
 TIMERESET= if(time() ==  $3,boolean(1),TIMERESET);
 TIMERESET= if(time() ==  $4,boolean(1),TIMERESET);
 TIMERESET= if(time() ==  $5,boolean(1),TIMERESET);
 TIMERESET= if(time() ==  $6,boolean(1),TIMERESET);
 TIMERESET= if(time() ==  $7,boolean(1),TIMERESET);
 TIMERESET= if(time() ==  $8,boolean(1),TIMERESET);
 TIMERESET= if(time() ==  $9,boolean(1),TIMERESET);
 TIMERESET= if(time() == $10,boolean(1),TIMERESET);
 TIMERESET= if(time() == $11,boolean(1),TIMERESET);
 TIMERESET= if(time() == $12,boolean(1),TIMERESET);
  PRECMON= if(TIMERESET,scalar(0),PRECMON);
  ETPMON= if(TIMERESET,scalar(0),ETPMON);
  QWMON= if(TIMERESET,scalar(0),QWMON);
  Q1MON= if(TIMERESET,scalar(0),Q1MON);
  Q2MON= if(TIMERESET,scalar(0),Q2MON);
  Q3MON= if(TIMERESET,scalar(0),Q3MON);
  EWATMON= if(TIMERESET,scalar(0),EWATMON);
  EACTMON=  if(TIMERESET,scalar(0),EACTMON);
  ESACTMON= if(TIMERESET,scalar(0),ESACTMON);
  TACTMON= if(TIMERESET,scalar(0),TACTMON);	
#-budget check
#STOT_ACT	(km3)		total active storage (e.g, excluding snow cover) to decide whether equilibrium has been
#				achieved (also included STOT_AVG, average total storage in (m) for soil,
#				which excludes snow accumulation)
#MBE		(m):		local mass balance error
#MBR		(-):		total mass balance error per catchment, relative to total input
 PTOT= PTOT+PRPTOT;
 report (rep1) ETOT= ETOT+(FRACWAT*EWAT+(1-FRACWAT)*EACT);
 QTOT= QTOT+QLOC;
 SLOC= (1-FRACWAT)*(S1+S2+S3+INTS+SC+SCF);
 report (rep1) MBE= PTOT+SLOCINI-(ETOT+QTOT+SLOC);
 report (rep1) STOT_ACT= if(LANDMASK,1E-9,1E-9)*maptotal((1-FRACWAT)*CELLAREA*(S1+S2+S3+INTS));
 STOT_CUM= STOT_CUM+(S1+S2+S3);
 report (rep1) STOT_AVG= STOT_CUM/time();

 INTOT= catchmenttotal(1E-9*CELLAREA*(SLOCINI+PTOT),LDD);
 OUTTOT= catchmenttotal(1E-9*CELLAREA*(POTABSTR_S0+ABSTR_S3+ETOT+QTOT+SLOC),LDD);
 report (rep1) MBR= 1-if(INTOT>0,(INTOT-OUTTOT)/INTOT,0);
#report (rep1) MBRTSS= timeoutput(STATIONS,MBR);
#report (rep1) MBRMONTSS= timeoutput(MONITORSTATIONS,MBR);

#-report of on total volumes as timeseries per catchment
#report (rep1) PTOTTSS = timeoutput(STATIONS,catchmenttotal(1E-9*CELLAREA*(SLOCINI+PTOT),LDD));
#report (rep1) ETOTTSS = timeoutput(STATIONS,catchmenttotal(1E-9*CELLAREA*ETOT,LDD));
#report (rep1) QTOTTSS = timeoutput(STATIONS,catchmenttotal(1E-9*CELLAREA*QTOT,LDD));
#report (rep1) STOTTSS = timeoutput(STATIONS,catchmenttotal(1E-9*CELLAREA*SLOC,LDD));

#------------
#Reporting timeseries
#--------------------
#report (rep1) SCTSS       = timeoutput(MONITORSTATIONS,SC+SCF);
#report (rep1) INTSTSS     = timeoutput(MONITORSTATIONS,INTS); 
#report (rep1) EACTTSS     = timeoutput(MONITORSTATIONS,EACT); 
#report (rep1) S1TSS       = timeoutput(MONITORSTATIONS,S1); 
#report (rep1) S2TSS       = timeoutput(MONITORSTATIONS,S2); 
#report (rep1) S3TSS       = timeoutput(MONITORSTATIONS,S3); 
#report (rep1) Q1TSS       = timeoutput(MONITORSTATIONS,Q1); 
#report (rep1) Q1STSS      = timeoutput(MONITORSTATIONS,Q1S); 
#report (rep1) Q2TSS       = timeoutput(MONITORSTATIONS,Q2); 
#report (rep1) Q3TSS       = timeoutput(MONITORSTATIONS,Q3); 
#report (rep1) QCHANNELTSS = timeoutput(STATIONS,QCHANNEL);
 
#--------------------
#Temporary Reports
#--------------------
#report (rep1) results\relstorsoil.tss= timeoutput(MONITORSTATIONS,(S1+S2)/WMAX[SHORT]);
#report (rep1) results\psi1.tss= timeoutput(MONITORSTATIONS,PSI1);
#report (rep1) results\psi2.tss= timeoutput(MONITORSTATIONS,PSI2);
#report (rep1) results\grad.tss= timeoutput(MONITORSTATIONS,GRAD);
#report (rep1) results\theff1.tss= timeoutput(MONITORSTATIONS,S1/SC1[SHORT]);
#report (rep1) results\theff2.tss= timeoutput(MONITORSTATIONS,S2/SC2[SHORT]);
#report (rep1) results\cr1.tss= timeoutput(MONITORSTATIONS,CR1);
#report (rep1) results\cr2.tss= timeoutput(MONITORSTATIONS,CR2);
#report (rep2) results\th1x= S1/SC1[SHORT];
#report (rep2) results\th2x= S2/SC2[SHORT];
