#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


binding

#Maps & TSS - input
#-General
 Duration= $2;						#timestep in days
 LANDMASK= $3;						#clone map representing landmask of earth surface
 CELLAREA= $4;						#surface (m2) of cell covered by total land surface

#-Surface water
 LDD= $5;													#local drainage direction map
 FRACWAT= $6;							#fraction of cell area covered by fresh water [m2/m2]
 KC_WATSTACK= $7;					#composite crop factor for channels and wetlands or lakes
 
#-Groundwater
 KQ3= $8;									#recession coefficient for store 3 (day-1): drainage
 SPECYIELD3= $9;					#specific yield for aquifer
 DZS3INFLUENCED= $10;			#additional zone influenced by groundwater
 DZREL0001= $91;					#maps of relative elevation above floodplain, in percent
 DZREL0005= $92;
 DZREL0010= $93;
 DZREL0020= $94;
 DZREL0030= $95;
 DZREL0040= $96;
 DZREL0050= $97;
 DZREL0060= $98;
 DZREL0070= $99;
 DZREL0080= $100;
 DZREL0090= $101;
 DZREL0100= $102;

#-Meteo
#-Precipitation [m/day], conversion constant and factor
 PRPSTACK= $11;					 	# input value
 PRPCONSTANT= $12;				# conversion constant, to get values as required
 PRPFACTOR= $13;					# conversion factor, to get values as required
#-Precipitation anomaly, multiplicative
 PRPANOSTACK= $14;			 	# input value
 PRPANOCONSTANT= $15;			# conversion constant, to get values as required
 PRPANOFACTOR= $16;				# conversion factor, to get values as required
#-Temperature [degC], conversion constant and factor 
 TASTACK= $17;
 TACONSTANT= $18;
 TAFACTOR= $19;
#-Temperature anomaly , additive
 TAANOSTACK= $20;
 TAANOCONSTANT= $21;
 TAANOFACTOR= $22;
#-Reference potential evapotranspiration [m/day], conversion constant and factor
 EVAPSTACK= $23;
 EVAPCONSTANT= $24;
 EVAPFACTOR= $25;
#-Reference potential evapotranspiration anomaly, multiplicative
 EVAPANOSTACK= $26;
 EVAPANOCONSTANT= $27;
 EVAPANOFACTOR= $28;

#-Vegetation parameters per cover type, interception and fractions read in per type (maximum 10 values)
 COVERTYPE= [
   $31= $41,
   $32= $42];							#array of cover type, e.g.: 1) short, 2) tall
#-fractional vegetation cover (m2/m2), interception storage (m per unit area)
# and crop factor per vegetation type
 CFSTACK_01= $51;								#-cover fraction
 CFSTACK_02= $52;
 CFSTACK_03= $53;
 CFSTACK_04= $54;
 CFSTACK_05= $55;
 CFSTACK_06= $56;
 CFSTACK_07= $57;
 CFSTACK_08= $58;
 CFSTACK_09= $59;
 CFSTACK_10= $60;
 SMAXSTACK_01= $61;							#-interception, maximum canopy storage
 SMAXSTACK_02= $62;
 SMAXSTACK_03= $63;
 SMAXSTACK_04= $64;
 SMAXSTACK_05= $65;
 SMAXSTACK_06= $66;
 SMAXSTACK_07= $67;
 SMAXSTACK_08= $68;
 SMAXSTACK_09= $69;
 SMAXSTACK_10= $70;
 KCSTACK_01= $71;								#-crop factor
 KCSTACK_02= $72;
 KCSTACK_03= $73;
 KCSTACK_04= $74;
 KCSTACK_05= $75;
 KCSTACK_06= $76;
 KCSTACK_07= $77;
 KCSTACK_08= $78;
 KCSTACK_09= $79;
 KCSTACK_10= $80;
#-table with parameterization per cover type; all variables read as maps
 COVERTABLE= $81;				#table with parameterization per cover type
#-parameter values: vegetation fraction
 VEGFRAC[COVERTYPE]= index(COVERTABLE);	#subdivision in cover type
#-parameter values: snow module
 TT[COVERTYPE]= index(COVERTABLE);  				#threshold temperature for freezing/thawing (?C)
 CFMAX[COVERTYPE]= index(COVERTABLE);				#degree-day factor (m??C-1?d-1) was: 0.0055
 SFCF[COVERTYPE]= index(COVERTABLE);				#snowfall correction factor (-)
 CWH[COVERTYPE]= index(COVERTABLE);					#water holding capacity snow cover (-)
 CFR[COVERTYPE]= index(COVERTABLE);					#refreezing coefficient (-)
#-topographical parameters per cover type
 LSLOPE[COVERTYPE]= index(COVERTABLE);			#slope length (m)
 TANSLOPE[COVERTYPE]= index(COVERTABLE);		#gradient of slope (m/m)
#-soil parameters per cover type
 PSI_FC[COVERTYPE]= index(COVERTABLE);			#matric suction at field capacity (m)
 PSI_50[COVERTYPE]= index(COVERTABLE);			#matric suction at which transpiration is halved (m)
 BCH_FACTOR[COVERTYPE]= index(COVERTABLE);	#factor for kr-relationship of Clapp & Hornberger (1978; default 2)
 BCH_ADD[COVERTYPE]= index(COVERTABLE);			#addition for kr-relationship of Clapp & Hornberger (1978; default 3)
 BCF[COVERTYPE]= index(COVERTABLE);					#beta coefficient of the improved Arno scheme
 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
#-soil parameters per cover type, constants 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
#-vegetation parameters per crop type
 KCMIN[COVERTYPE]= index(COVERTABLE);				#used to separate soil evaporation from transpiration
#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)
#Initial storages for local variables, cell averages reported by default
 S3_INI= $82;																#initial storage in lower store (m)

#Output: only maps are reported
#-representative for land cover dependent variables: note cover type is not needed to report local fluxes!
#-initial values
 SC_L= $89\sc;					#snow cover (m)
 SCF_L= $89\scf;				#liquid water stored in snow cover (m)
 INTS_L= $89\ints;			#interception storage (m)
 Q2_L= $89\q2;					#runoff from second store (flow) (m) 
 S1_L= $89\s1;					#storage in upper store (m)
 S2_L= $89\s2;					#storage in second store (m)
#-optional output: see matches below for cell-based values
 PRP_L= $89\prec;
 TA_L= $89\temp;
 ET_p= $89\epot;
 THETA1_L= $89\th1;
 THETA2_L= $89\th2;
 THETAS_L= $89\ths;
 W= $89\stors;
 SATFRAC_L= $89\satf;
 EACT_L= $89\eact;
 ES_p= $89\espot;
 ES_a= $89\esact;
 T1_p= $89\t1pot;
 T1_a= $89\t1act;
 T2_p= $89\t2pot;
 T2_a= $89\t2act;
 T_p= $89\tpot;
 T_a= $89\tact;
 P0_L= $89\p0;
 P1_L= $89\p1;
 P2_L= $89\p2;
 Q1_L= $89\q1;
 CR1_L= $89\cr1;
 CR2_L= $89\cr2;

#-representative for total cell: meteo
 PRPTOT= $89\prec;			#total precipitiation (m/day)
 TAVG= $89\temp;				#average air temperature (degC)
 EPOT= $89\epot;				#total potential evapotranspiration (m)
#-representative for total cell: states
 SC= $89\snowcov;				#snow cover (m)
 SCF= $89\snowliq;			#liquid water stored in snow cover (m)
 INTS= $89\intstor;			#interception storage (m)
 THETA1= $89\theta1x;		#volumetric moisture content in the first soil (m3/m3)
 THETA2= $89\theta2x;		#volumetric moisture content in the second soil (m3/m3)
 THETAS= $89\theta;			#volumetric moisture content in soil column (m3/m3)
 S1= $89\stor1x;				#water storage in upper store (m)
 S2= $89\stor2x;				#water storage in second store (m)
 SS= $89\stors;					#total water storage in soil (m)
 SATFRAC= $89\satf;			#fraction saturated area (-)
 CRFRAC= $89\crfrac;		#fraction of soil surface influenced by capillary rise
#-representative for total cell: fluxes
 EWAT= $89\ewat;				#potential evaporation over open water
 EACT= $89\eact;				#actual evapotranspiration (m)
 ESPOT= $89\espot;			#potential soil evaporation (m)
 ESACT= $89\esact;			#actual soil evaporation (m)
 T1POT= $89\t1pot;			#potential transpiration drawn from first soil layer (m)
 T1ACT= $89\t1act;			#idem, actual
 T2POT= $89\t2pot;			#potential transpiration drawn from second soil layer (m)
 T2ACT= $89\t2act;			#idem, actual
 TPOT= $89\tpot;				#potential transpiration drawn from total soil layer (m)
 TACT= $89\tact;				#idem, actual
 P0= $89\p0x;						#infiltration to first layer (m/day)
 P1= $89\p1x;						#percolation from first layer (m/day)
 P2= $89\p2x;           #percolation from second layer (m/day)
 CR1= $89\cr1x;					#capillary rise to first layer (m/day)
 CR2= $89\cr2x;         #capilary rise to second layer (m/day)
 Q1= $89\q1x;						#direct runoff (m)
 Q2=  $89\q2x;					#runoff from second store (flow) (m)
 Q3= $89\q3x;						#runoff from lower store (base flow) (m)
 QLOC= $89\qloc;				#specific runoff (m)
 R3= $89\r3x;						#recharge to third store
 S3= $89\stor3x;				#storage in lower store (m)
 QWAT= $89\qwat;				#change in storage of freshwater surface
 QAVG= $89\qaverage.map;#channel discharge - annual average (m/s)
#budget check
 MBE= $89\mbe.map;			#absolute local mass balance error (m)
 MBR= $89\mbcheck.map;	#relative, local mass balance error (-)

areamap
 LANDMASK;

timer
 1 $1 1;
 rep1= endtime;					#times to report TSS and initial maps
 rep2= 1+1..endtime;			#times to report soil states and fluxes


initial

#Patching LDD
 LDD= if(LANDMASK,cover(LDD,5));

#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]= BCH_FACTOR[TYPE]*BCH1[TYPE]+BCH_ADD[TYPE];		#Campbell's (1974) coefficient to calculate the relative,
   BCB2[TYPE]= BCH_FACTOR[TYPE]*BCH2[TYPE]+BCH_ADD[TYPE];		#unsaturated hydraulic conductivity
						#degree of saturation at field capacity for first and second layer
						#and corresponding unstaturated hydraulic conductivity
   THEFF1_FC[TYPE]= (PSI_FC[TYPE]/PSI_A1[TYPE])**(-1/BCH1[TYPE]);
   KTHEFF1_FC[TYPE]= max(0,THEFF1_FC[TYPE]**BCB1[TYPE]*KS1[TYPE]);
   THEFF2_FC[TYPE]= (PSI_FC[TYPE]/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[TYPE])/(LSLOPE[TYPE]*(1-THEFF2_FC[TYPE])*(THETASAT2[TYPE]-THETARES2[TYPE]));
   #-storage parameters related to Improved Arno Scheme
   #-total active water storage capacity over first and second layer and total soil profile (m)
   SC1[TYPE]= Z1[TYPE]*(THETASAT1[TYPE]-THETARES1[TYPE]);
   SC2[TYPE]= Z2[TYPE]*(THETASAT2[TYPE]-THETARES2[TYPE]); 
   WMAX[TYPE]= SC1[TYPE]+SC2[TYPE];
   # storage in top soil limiting bare soil evapotranspiration
   WSMAX[TYPE]= SC1[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];
   #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(RFW1[TYPE]+RFW2[TYPE] > 0,  
     (SC1[TYPE]*RFW1[TYPE]*(PSI_50[TYPE]/PSI_A1[TYPE])**(-1/BCH1[TYPE])+
     SC2[TYPE]*RFW2[TYPE]*(PSI_50[TYPE]/PSI_A2[TYPE])**(-1/BCH2[TYPE]))/
     (SC1[TYPE]*RFW1[TYPE]+SC2[TYPE]*RFW2[TYPE]),0.5);
   BCH_50[TYPE]= if(RFW1[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 fluxes for average discharge and budget check
 QCUM= scalar(0);
 STOT_CUM= scalar(0);
 SLOCINI= (1-FRACWAT)*(S1+S2+S3+INTS+SC+SCF);
 PTOT= scalar(0);
 ETOT= scalar(0);
 QTOT= scalar(0);
 
dynamic
#----------------
# Meteo
#----------------
#Meteorological input as total/average per time step
# PRPORG	(m):	total precipitation from original forcing
# TA			(?C):	average temperature
# EVAP		(m):	reference potential evapotranspiration
 PRPORG= PRPCONSTANT+PRPFACTOR*
  timeinput(PRPSTACK);
 PRPANO= PRPANOCONSTANT+PRPANOFACTOR*
  timeinputsparse(PRPANOSTACK);
 PRPORG= max(0.,PRPANO*PRPORG)*Duration*timeslice();
 TA= TACONSTANT+TAFACTOR*timeinput(TASTACK);
 TAANO= TAANOCONSTANT+TAANOFACTOR*
  timeinputsparse(TAANOSTACK);
 TA= TAANO+TA;
 TA= TA+20; #!
 EVAP= EVAPCONSTANT+EVAPFACTOR*
  timeinput(EVAPSTACK);
 EVAPANO= EVAPANOCONSTANT+EVAPANOFACTOR*
  timeinputsparse(EVAPANOSTACK);
 EVAP= 0*max(0.,EVAPANO*EVAP)*Duration*timeslice(); #!
 
#----------------
# 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);

#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);
 PRPTOT= scalar(0);
 TAVG= scalar(0); 
 SC= scalar(0);
 SCF= scalar(0);
 INTS= scalar(0);
 EPOT= scalar(0);
 ESPOT= scalar(0);
 ESACT= scalar(0);
 T1POT= scalar(0);
 T1ACT= scalar(0);
 T2POT= scalar(0);
 T2ACT= scalar(0);
 TPOT= scalar(0);
 TACT= scalar(0);
 EACT= scalar(0);
 SATFRAC= scalar(0);
 WACT= scalar(0);
 THETA1= scalar(0);
 THETA2= scalar(0);
 THETAS= scalar(0);
 S1= scalar(0);
 S2= scalar(0);
 SS= 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);

 foreach TYPE in COVERTYPE{
  #-Dynamic vegetation properties: update vegetation counter and set the following variables
  #TYPECOUNT  : vegetation counter
  #CFRAC   (-):	vegetation cover fraction
  #KC	     (-):	crop factor
  #INTCMAX (m):	maximum interception storage capacity
   TYPECOUNT= TYPECOUNT+1;
   CFRAC[TYPE]= timeinputsparse(CFSTACK_01);
   CFRAC[TYPE]= if(TYPECOUNT ==  2, timeinputsparse(CFSTACK_02), CFRAC[TYPE]);
   CFRAC[TYPE]= if(TYPECOUNT ==  3, timeinputsparse(CFSTACK_03), CFRAC[TYPE]);
   CFRAC[TYPE]= if(TYPECOUNT ==  4, timeinputsparse(CFSTACK_04), CFRAC[TYPE]);
   CFRAC[TYPE]= if(TYPECOUNT ==  5, timeinputsparse(CFSTACK_05), CFRAC[TYPE]);
   CFRAC[TYPE]= if(TYPECOUNT ==  6, timeinputsparse(CFSTACK_06), CFRAC[TYPE]);
   CFRAC[TYPE]= if(TYPECOUNT ==  7, timeinputsparse(CFSTACK_07), CFRAC[TYPE]);
   CFRAC[TYPE]= if(TYPECOUNT ==  8, timeinputsparse(CFSTACK_08), CFRAC[TYPE]);
   CFRAC[TYPE]= if(TYPECOUNT ==  9, timeinputsparse(CFSTACK_09), CFRAC[TYPE]);
   CFRAC[TYPE]= if(TYPECOUNT == 10, timeinputsparse(CFSTACK_10), CFRAC[TYPE]);
   KC[TYPE]= timeinputsparse(KCSTACK_01);
   KC[TYPE]= if(TYPECOUNT ==  2, timeinputsparse(KCSTACK_02), KC[TYPE]);
   KC[TYPE]= if(TYPECOUNT ==  3, timeinputsparse(KCSTACK_03), KC[TYPE]);
   KC[TYPE]= if(TYPECOUNT ==  4, timeinputsparse(KCSTACK_04), KC[TYPE]);
   KC[TYPE]= if(TYPECOUNT ==  5, timeinputsparse(KCSTACK_05), KC[TYPE]);
   KC[TYPE]= if(TYPECOUNT ==  6, timeinputsparse(KCSTACK_06), KC[TYPE]);
   KC[TYPE]= if(TYPECOUNT ==  7, timeinputsparse(KCSTACK_07), KC[TYPE]);
   KC[TYPE]= if(TYPECOUNT ==  8, timeinputsparse(KCSTACK_08), KC[TYPE]);
   KC[TYPE]= if(TYPECOUNT ==  9, timeinputsparse(KCSTACK_09), KC[TYPE]);
   KC[TYPE]= if(TYPECOUNT == 10, timeinputsparse(KCSTACK_10), KC[TYPE]);
   INTCMAX[TYPE]= timeinputsparse(SMAXSTACK_01);
   INTCMAX[TYPE]= if(TYPECOUNT ==  2, timeinputsparse(SMAXSTACK_02), INTCMAX[TYPE]);
   INTCMAX[TYPE]= if(TYPECOUNT ==  3, timeinputsparse(SMAXSTACK_03), INTCMAX[TYPE]);
   INTCMAX[TYPE]= if(TYPECOUNT ==  4, timeinputsparse(SMAXSTACK_04), INTCMAX[TYPE]);
   INTCMAX[TYPE]= if(TYPECOUNT ==  5, timeinputsparse(SMAXSTACK_05), INTCMAX[TYPE]);
   INTCMAX[TYPE]= if(TYPECOUNT ==  6, timeinputsparse(SMAXSTACK_06), INTCMAX[TYPE]);
   INTCMAX[TYPE]= if(TYPECOUNT ==  7, timeinputsparse(SMAXSTACK_07), INTCMAX[TYPE]);
   INTCMAX[TYPE]= if(TYPECOUNT ==  8, timeinputsparse(SMAXSTACK_08), INTCMAX[TYPE]);
   INTCMAX[TYPE]= if(TYPECOUNT ==  9, timeinputsparse(SMAXSTACK_09), INTCMAX[TYPE]);
   INTCMAX[TYPE]= if(TYPECOUNT == 10, timeinputsparse(SMAXSTACK_10), INTCMAX[TYPE]);
   
   INTCMAX[TYPE]= scalar(0); #-smax #!
	
	#-dynamic vegetation types set, continue with water balance

  #-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]= KC[TYPE]*EVAP;	
   ES_p[TYPE]= KCMIN[TYPE]*EVAP;
   T_p[TYPE]= KC[TYPE]*EVAP-ES_p[TYPE];

  #-Temperature: mean value passed
  TA_L[TYPE]= TA;
	#-Precipitation: partitioning rain and snow and converting evapotransiration
  # PRP_L		(m):	liquid precipitation
  # SNOW_L	(m):	snow in water equivalent
   SNOW_L= if(TA_L[TYPE]<TT[TYPE],PRPORG,0);
	 PRP_L[TYPE]= max(0,PRPORG-SNOW_L);
	 SNOW_L= SFCF[TYPE]*SNOW_L;
	 PRP_L[TYPE]= PRP_L[TYPE]+SNOW_L;

  #-Interception
  #ICC		(m):	equivalent interception storage capacity
  #INTS		(m):	interception storage
  #PTF_L	(m):	precipitation passing the canopy
   ICC[TYPE]= CFRAC[TYPE]*INTCMAX[TYPE];
   PTF_L= (1-CFRAC[TYPE])*PRP_L[TYPE]+max(CFRAC[TYPE]*PRP_L[TYPE]+INTS_L[TYPE]-ICC[TYPE],0);
   INTS_L[TYPE]= max(0,INTS_L[TYPE]+PRP_L[TYPE]-PTF_L);
   SNOW_L= min(PTF_L,SNOW_L*if(PRP_L[TYPE]>0,PTF_L/PRP_L[TYPE],0)); #**
   PTF_L= max(0,PTF_L-SNOW_L); #**
   EACT_L[TYPE]= min(INTS_L[TYPE],(T_p[TYPE]*if(ICC[TYPE]>0,INTS_L[TYPE]/ICC[TYPE],0)**(2/3)));
   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_L[TYPE]<=TT[TYPE],CFR[TYPE]*SCF_L[TYPE],-min(SC_L[TYPE],max(TA_L[TYPE]-TT[TYPE],0)*CFMAX[TYPE]*Duration*timeslice()));
   SC_L[TYPE]= SC_L[TYPE]+DSC[TYPE]+SNOW_L;
   SCF_L[TYPE]= SCF_L[TYPE]-DSC[TYPE]+PTF_L;
   Pn= max(0,SCF_L[TYPE]-CWH[TYPE]*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]);
   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];

  #-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_L= max(0,S1_L[TYPE]/SC1[TYPE]);
   THEFF2_L= max(0,S2_L[TYPE]/SC2[TYPE]);
   PSI1= PSI_A1[TYPE]*max(0.01,THEFF1_L)**-BCH1[TYPE]; #*
   PSI2= PSI_A2[TYPE]*max(0.01,THEFF2_L)**-BCH2[TYPE]; #*
   GRAD= max(0,2*(PSI1-PSI2)/(Z1[TYPE]+Z2[TYPE])-1); #*
   KTHEFF1= max(0,THEFF1_L**BCB1[TYPE]*KS1[TYPE]);
   KTHEFF2= max(0,THEFF2_L**BCB2[TYPE]*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.0);
   WACT_L= (BCF[TYPE]+1)*WMAX[TYPE]-BCF[TYPE]*WMIN[TYPE]-(BCF[TYPE]+1)*WRANGE[TYPE]*WFRACB;
   Q1_L[TYPE]= 0*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))); #* #!
   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]);
   T1_a[TYPE]= FRACTA[TYPE]*T1_p[TYPE];
   T2_a[TYPE]= FRACTA[TYPE]*T2_p[TYPE];
   T_a[TYPE]= T1_a[TYPE]+T2_a[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_L > THEFF1_FC[TYPE],min(max(0,THEFF1_L-THEFF1_FC[TYPE])*SC1[TYPE],
     P1_L[TYPE]),P1_L[TYPE])+max(0,P0_L[TYPE]-(SC1[TYPE]-S1_L[TYPE])); #*
   CR1_L[TYPE]= 0*min(max(0,THEFF1_FC[TYPE]-THEFF1_L)*SC1[TYPE],
     KTHVERT*GRAD*Duration*timeslice()); #* #!
   CR2_L[TYPE]= 0*0.5*(SATFRAC_L[TYPE]+CRFRAC)*min((1-THEFF2_L)*sqrt(KS2[TYPE]*KTHEFF2)*Duration*timeslice(),
     max(0,THEFF2_FC[TYPE]-THEFF2_L)*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]+T1_a[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];
   T1_a[TYPE]= ADJUST*T1_a[TYPE];
   P1_L[TYPE]= ADJUST*P1_L[TYPE];
   # second layer
   ADJUST= T2_a[TYPE]+P2_L[TYPE]+Q2_L[TYPE];
   ADJUST= if(ADJUST>0,min(1,max(S2_L[TYPE]+P1_L[TYPE],0)/ADJUST),0);
   T2_a[TYPE]= ADJUST*T2_a[TYPE];
   P2_L[TYPE]= ADJUST*P2_L[TYPE];
   Q2_L[TYPE]= ADJUST*Q2_L[TYPE];
   CR2_L[TYPE]= min(VEGFRAC[TYPE]*S3,CR2_L[TYPE]);
   CR1_L[TYPE]= min(max(0,S2_L[TYPE]+P1_L[TYPE]-(T2_a[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]+T2_a[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]+T1_a[TYPE]+ES_a[TYPE])); #*
   S2_L[TYPE]= min(S2_L[TYPE],SC2[TYPE]);
   Q1_L[TYPE]= Q1_L[TYPE]+max(0,S1_L[TYPE]-SC1[TYPE]);
   S1_L[TYPE]= min(S1_L[TYPE],SC1[TYPE]);
  #-total actual evapotranspiration
   EACT_L[TYPE]= EACT_L[TYPE]+ES_a[TYPE]+T1_a[TYPE]+T2_a[TYPE];

	#-update volumetric moisture content (m3/m3)
	# THETA1_L	:	volumetric moisture content in the first soil layer (m3/m3)
	# THETA2_L	:	volumetric moisture content in the second soil layer (m3/m3)
	# THETAS_L	:	volumetric moisture content in the total soil layer (m3/m3)
	# W					:	total water content in the soil layer (m)
   THETA1_L[TYPE]= THETARES1[TYPE]+(THETASAT1[TYPE]-THETARES1[TYPE])*max(0,S1_L[TYPE]/SC1[TYPE]);
   THETA2_L[TYPE]= THETARES1[TYPE]+(THETASAT1[TYPE]-THETARES1[TYPE])*max(0,S2_L[TYPE]/SC2[TYPE]);
   THETAS_L[TYPE]= (Z1[TYPE]*THETA1_L[TYPE]+Z2[TYPE]*THETA2_L[TYPE])/max(0.001,Z1[TYPE]+Z2[TYPE]);
   W[TYPE]= max(0,S1_L[TYPE]+S2_L[TYPE]);
     
  #-adding local fluxes and states relative to vegetation fractions
  # as a function of RTN to correct for vegetation presence
	 PRPTOT= PRPTOT+VEGFRAC[TYPE]*PRP_L[TYPE];
	 TAVG= TAVG+VEGFRAC[TYPE]*TA_L[TYPE];
   EPOT= EPOT+VEGFRAC[TYPE]*ET_p[TYPE];	 
   SC= SC+VEGFRAC[TYPE]*SC_L[TYPE];
   SCF= SCF+VEGFRAC[TYPE]*SCF_L[TYPE];
   INTS= INTS+VEGFRAC[TYPE]*INTS_L[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]*T1_p[TYPE];
   T1ACT= T1ACT+VEGFRAC[TYPE]*T1_a[TYPE];
   T2POT= T2POT+VEGFRAC[TYPE]*T2_p[TYPE];
   T2ACT= T2ACT+VEGFRAC[TYPE]*T2_a[TYPE];
   TPOT= T1POT+T2POT;
   TACT= T1ACT+T2ACT;
   SATFRAC= SATFRAC+VEGFRAC[TYPE]*SATFRAC_L[TYPE];
   WACT= WACT+VEGFRAC[TYPE]*WACT_L;
   THETA1= THETA1+VEGFRAC[TYPE]*THETA1_L[TYPE];
   THETA2= THETA2+VEGFRAC[TYPE]*THETA2_L[TYPE];
   THETAS= THETAS+VEGFRAC[TYPE]*THETAS_L[TYPE];
   S1= S1+VEGFRAC[TYPE]*S1_L[TYPE];
   S2= S2+VEGFRAC[TYPE]*S2_L[TYPE];
   SS= SS+VEGFRAC[TYPE]*W[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];
   
  #-report values per land cover type
  #-initial and possible output
   report (rep2) SC_L[TYPE]= SC_L[TYPE];
   report (rep2)  SCF_L[TYPE]= SCF_L[TYPE];
   report (rep2)  INTS_L[TYPE]= INTS_L[TYPE];
   report (rep2) S1_L[TYPE]= S1_L[TYPE];
   report (rep2)  S2_L[TYPE]= S2_L[TYPE];
  #-possible output only
   report (rep2) PRP_L[TYPE]= PRP_L[TYPE];
   report (rep2) TA_L[TYPE]= TA_L[TYPE];
   report (rep2) ET_p[TYPE]= ET_p[TYPE];
   report (rep2) THETA1_L[TYPE]= THETA1_L[TYPE];
   report (rep2) THETA2_L[TYPE]= THETA2_L[TYPE];
   report (rep2) THETAS_L[TYPE]= THETAS_L[TYPE];
   report (rep2) W[TYPE]= W[TYPE];
   report (rep2) SATFRAC_L[TYPE]= SATFRAC_L[TYPE];
   report (rep2) ES_p[TYPE]= ES_p[TYPE];
   report (rep2) ES_a[TYPE]= ES_a[TYPE];
   report (rep2) T1_p[TYPE]= T1_p[TYPE];
   report (rep2) T1_a[TYPE]= T1_a[TYPE];
   report (rep2) T2_p[TYPE]= T2_p[TYPE];
   report (rep2) T2_a[TYPE]= T2_a[TYPE];
   report (rep2) T_p[TYPE]= T_p[TYPE];
   report (rep2) T_a[TYPE]= T_a[TYPE];
   report (rep2) P0_L[TYPE]= P0_L[TYPE];
   report (rep2) P1_L[TYPE]= P1_L[TYPE];
   report (rep2) P2_L[TYPE]= P2_L[TYPE];
   report (rep2) CR1_L[TYPE]= CR1_L[TYPE];
   report (rep2) CR2_L[TYPE]= CR2_L[TYPE];
   report (rep2) Q1_L[TYPE]= Q1_L[TYPE];   
   report (rep2) Q2_L[TYPE]= Q2_L[TYPE];
  
  }#-end of land cover type loop
 
#--------------------------
#Overall fluxes third store
#--------------------------
#Third reservoir
#R3		(m):	groundwater recharge
#S3		(m):	storage in third store, updated with current fluxes
#Q3		(m):	discharge from third reservoir, based on storage previous timestep
 R3= P2-CR2;
 S3= max(0,S3+P2-CR2);
 Q3= min(S3,KQ3*S3*Duration*timeslice());
 S3= max(0,S3-Q3);

#--------------------
# Fresh water surface
#--------------------
# EWAT		(m):	potential evapotranspiration imposed on water surface
# QWAT	  (m):	local change in storage of fresh water surface (can be negative)
 EWAT= timeinputsparse(KC_WATSTACK)*EVAP;
 QWAT= if(LANDMASK,PRPTOT-EWAT);

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

#------------
#Budget check
#------------
#PTOT		(m):		total accumulated precipitation
#ETOT		(m):		total accumulated evapotranspiration
#QTOT		(m):		total accumulated local discharge
#INTOT, OUTTOT	(km3):		total incoming and outgoing water volumes per catchment
#SLOC, SLOCINI	(m):		local storage at any timestep and initially
#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;
 ETOT= ETOT+(FRACWAT*EWAT+(1-FRACWAT)*EACT);
 QTOT= QTOT+((1-FRACWAT)*QLOC+FRACWAT*QWAT);
 SLOC= (1-FRACWAT)*(S1+S2+S3+INTS+SC+SCF);
 MBE= PTOT+SLOCINI-(ETOT+QTOT+SLOC);
 STOT_ACT= if(LANDMASK,1E-9,1E-9)*maptotal((1-FRACWAT)*CELLAREA*(S1+S2+S3+INTS));
 INTOT= catchmenttotal(1E-9*CELLAREA*(SLOCINI+PTOT),LDD);
 OUTTOT= catchmenttotal(1E-9*CELLAREA*(ETOT+QTOT+SLOC),LDD);
 MBR= 1-if(INTOT>0,(INTOT-OUTTOT)/INTOT,0);

#------------------------
# Reports: overall values
#------------------------
#-Meteo
 report (rep2) PRPTOT= PRPTOT;
 report (rep2) EPOT= EPOT;
 report (rep2) TAVG= TAVG;
#-states: snow pack
 report (rep2) SC= SC;
 report (rep2) SCF= SCF;
#-states: soil moisture and saturated area
 report (rep2) S1= S1;
 report (rep2) S2= S2;
 report (rep2) SS= SS;
 report (rep2) THETA1= THETA1;
 report (rep2) THETA2= THETA2;
 report (rep2) THETAS= THETAS;
 report (rep2) SATFRAC= SATFRAC;
#-groundwater: storage, recharge and area influenced by capillary rise
 report (rep2) S3= S3;
 report (rep2) R3= R3;
 report (rep2) CRFRAC= CRFRAC;
#-fluxes: evapotranspiration
 report EWAT= EWAT;
 report EACT= EACT;
 report ESPOT= ESPOT;
 report ESACT= ESACT;
 report T1POT= T1POT;
 report T2POT= T2POT;
 report T1ACT= T1ACT;
 report T2ACT= T2ACT;
 report TPOT= TPOT;
 report TACT= TACT;
#-fluxes: percolation and capillary rise
 report (rep2) P0= P0;
 report (rep2) P1= P1;
 report (rep2) P2= P2;
 report (rep2) CR1= CR1;
 report (rep2) CR2= CR2;
#-runoff, discharge and surface water
 report (rep2) Q1= Q1;
 report (rep2) Q2= Q2;
 report (rep2) Q3= Q3;
 report (rep2) QLOC= QLOC;
 report (rep2) QWAT= QWAT;
 report (rep1) QAVG= QAVG;

#-budget checks
 report (rep1) MBE= MBE;
 report (rep1) MBR= MBR;

#-all output reported, end of file
