#include "StorageDeficit_VadoseZone.h"

//Note: Catchment average soil moisture deficit (StorageDeficit_VadoseZone_m, AveSMD) determines water storage in each Topographic Index bin
//Note: Saturation excess runoff is an outcome of AveSMD and local SMD values, based on variable source area and topographic index (TI) theory
//Note: The catchment Bulk Area uses local SMD values to determine saturation excess in each TI bin, StorageDeficit_VadoseZone_TI_m
//Note: The Green Infrastructure (GI) Area uses AveSMD via Depth_to_Groundwater_Table_m = StorageDeficit_VadoseZone_m / Porosity_Drainable_m3pm3
//Note: Green Infrastructure (GI) uses Depth_to_Groundwater_Table_m in Groundwater_StormwaterDevice.cpp to determine saturation by groundwater

//Note: Theory of updating AveSMD (m) with drainage from vadose zone and runoff from saturated zone based on Wolock (1993).
//Note: Soil volumetric water has 3 zones ranging from wilting point, 0wp to field capacity, 0fc, to saturation, 0sat
//Note: 1) evapotranspiration dynamically contains 0wp to 0sat; 0sat-0fc (gravitational water) releases to lower vadose or saturated zone
//Note: 2) vadose zone dynamically contains 0sat to 0fc; 0sat-0fc drains to saturated zone; it is ground surface to water table
//Note: 3) saturated zone generates subsurface runoff; it can rise into vadose and evapotranspiration zones to fill gravitational voids
//Note: The evapotranspiration and vadose zones are conceptually situated alongside each other, each filling a different volumetric component
//Note: The evapotranspiration zone extends to the limit of evapotranspiration, only processing capillary held water
//Note: The vadose zone is represented by average soil moisture deficit, AveSMD, StorageDeficit_VadoseZone_TI_m, or StorageDeficit_VadoseZone_m
//Note: For ts=0, AveSMD (m) defined in TerrainProcessor::getAveSoilMoistureDeficit using Q0_mph, T0_m2ph, mM, nN with Exponential or Power Decay
//Note: For ts>0, StorageDeficit_VadoseZone_m & StorageDeficit_VadoseZone_TI_m defined in Flux_LateralGroundWater::calculate as above
//Note: StorageDeficit_FieldCapacity_SoilEvapZone_Max_m (m) = Evapotranspiration_Depth_m * (Soil_FieldCapacity_m3pm3 - Soil_WiltingPoint_m3pm3)
//Note: StorageDeficit_VadoseZone_m (m) = (Soil_SaturationPoint_m3pm3 - Soil_FieldCapacity_m3pm3) * Depth_SoilColumn_Ave_m
//Reference: Wolock, D. M. (1993). Simulating the variable-source-area concept of watershed hydrology with TOPMODEL (USGS Water-Resources Investigation Report 93-4124). 
//Note: Neither the Topmodel code of Beven nor or Wolock (1993), compute return flow, despite presentation of Eq 21 in Wolock (1993); it is impractical to separate its source as aquifer or precipitation
//Note: Wolock (1993) Eq 21 mentions return flow is not simulated when StorageDeficit_VadoseZone_TI_m < 0, given saturation may represent impounded water, not laterally redistributed groundwater breaching surface
//Note: TOPMODEL code is archived at end of this .cpp file

void StorageDeficit_VadoseZone::calculate(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
    double ratio_FolderArea_to_CatchmentArea = beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2") / input->SimulationNumericalParams["CatchmentArea_m2"];
    //Catchment_PerviousArea_AffectingStorageDeficit_m2 (m2) contains entire catchment pervious area draining to vadose zone
    //Note: Catchment_PerviousArea_AffectingStorageDeficit_m2 influences the catchment variable, StorageDeficit_VadoseZone_m
    double Catchment_PerviousArea_AffectingStorageDeficit_m2 = input->RepoDict["Catchment_PerviousArea_AffectingStorageDeficit_m2"];

    double Folder_Pervious_Area_m2 = beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2") * beC->by_key(MapPixel_ID, DataFolder_ID, "PerviousCover_frac");
    double Drainage_VadoseZone_m = 0.0;
    double Runoff_Subsurface_m = 0.0;  

    //If Type is BulkArea then 
    if (beC->by_key_str(MapPixel_ID, DataFolder_ID, "Type") == "BulkArea") {
        //Drainage_VadoseZone_m (m) is catchment pervious area depth, arriving from DrainageToSaturationZone.cpp as pervious cover depth, then 
        //Note: ... multiplied by Folder_Pervious_Area_m2 (m2) to get volume (m3), then divided by Catchment_PerviousArea_AffectingStorageDeficit_m2 (m2) to get catchment pervious area depth
        Drainage_VadoseZone_m = input->SafeDivide((beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_VadoseZone_m") * Folder_Pervious_Area_m2), Catchment_PerviousArea_AffectingStorageDeficit_m2);
    }
    //Else If Type is GI then implement GI drainage
    else if (input->RepoDict["Flag_GI_Simulated"] == 1 && beC->by_key_str(MapPixel_ID, DataFolder_ID, "Type") != "BulkArea" && Folder_Pervious_Area_m2 > 0) {
        //Drainage_VadoseZone_m (m) is catchment pervious area depth, arriving from Flux_Manager_StormwaterDevice.cpp as GI_Drainage_VadoseZone_m3 (m3), then
        //Note: ... divided by Catchment_PerviousArea_AffectingStorageDeficit_m2 (m2) to get catchment pervious area depth
        Drainage_VadoseZone_m = input->SafeDivide(beC->by_key(MapPixel_ID, DataFolder_ID, "GI_Drainage_VadoseZone_m3"), Catchment_PerviousArea_AffectingStorageDeficit_m2);
    }
    //Runoff_Subsurface_m (m) is catchment pervious area depth, arriving from RunoffSubsurface.cpp as pervious cover depth, then
    //Note: ... multiplied by Folder_Pervious_Area_m2 (m2) to get volume (m3), then divided by Catchment_PerviousArea_AffectingStorageDeficit_m2 (m2) to get catchment pervious area depth
    Runoff_Subsurface_m = input->SafeDivide((beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Subsurface_m") * Folder_Pervious_Area_m2), Catchment_PerviousArea_AffectingStorageDeficit_m2);
    //StorageDeficit_VadoseZone_m (m), pervious area depth, is decreased by Drainage_VadoseZone_m (m) and increased by Runoff_Subsurface_m (m) from prior time step 
    //Note: Groundwater table is lifted by Drainage_VadoseZone_m (m) and lowered by Runoff_Subsurface_m (m) 
    //Note: Terms Drainage_VadoseZone_m and Runoff_Subsurface_m are prior time step, based on sequence of functions in SimulationCoordinator::RunStatisticalModel
    //Note: StorageDeficit_VadoseZone_m iteratively updated while passing through For Loop of DataDrawers, and it replaces StorageDeficit_VadoseZone_TS_Prior_m at start of time step
    input->RepoDict["StorageDeficit_VadoseZone_m"] = input->RepoDict["StorageDeficit_VadoseZone_m"] - Drainage_VadoseZone_m + Runoff_Subsurface_m;
    //TOPMODEL Code: SBAR = SBAR - QUZ + QB, where SBAR is AveSMD (m), QUZ is Drainage_VadoseZone_m, and QB is Runoff_Subsurface_m

    //Depth_to_Groundwater_Table_m (m), scaled to pervious area, is quotient of StorageDeficit_VadoseZone_m (m) and Porosity_Drainable_m3pm3 (m3/m3)
    //Note: Eq 16 Wolock (1993) gives Z = SMD/ne, where SMD (m) = soil moisture deficit, Z (m) = depth to groundwater table, ne = drainable porosity
    //Note: Green Infrastructure (GI) uses Depth_to_Groundwater_Table_m in Groundwater_StormwaterDevice.cpp to determine saturation by groundwater
    input->InputXml["Depth_to_Groundwater_Table_m"] = input->RepoDict["StorageDeficit_VadoseZone_m"] / input->InputXml["Porosity_Drainable_m3pm3"];
}

/*  USGS TOPMODEL Code
C****************************************************************
CCTE TOPMODEL code developed by Keith Beven and used by USGS D.M.Wolock https ://pubs.er.usgs.gov/publication/wri934124
CCTE Wolock, D.M. (1993).Simulating the variable - source - area concept of watershed hydrology with TOPMODEL, USGS WRI Report 93 - 4124.
CCTE The algorithm requires the topographic index (TI) bins (ia) are sorted from highest (wettest) to lowest TI values
CCTE The algorithm requires when ia = 0, the TI bin TI_Value[0] should be the upper limit and its TI_Area_frac[0] = 0
CCTE The algorithm presumes for each ia, the area of catchment between discrete TI values TI_Value[ia] and TI_Value[ia+1] equals TI_Area_frac[ia]
CCTE The algorithm presumes for each ia, the area of catchment bounding TI_Value[ia] equals (TI_Area_frac[ia] + TI_Area_frac[ia+1])/2
CCTE The algorithm in C++ has TI values for ia = 0 to ia = [Size_TI_Value_Vector - 1], yet requires use of TI_Area_frac[Size_TI_Value_Vector] = 0 in the above fractional area equation

CCTE Note infiltration routine updates DF, which updates P in soil moisture routine, which updates SRZand QUZ
CCTE Note deficit is equivalent to capacity, when 0 there is no more storage.
CTE Routine loops through all IA, the TI_Value value within topographic index histogramand adds water to root zone, drains unsaturated zone, identify surface runoff, removes evapotranspiration, removes surface flow from fully saturated area, compute saturated subsurface drainage
CTE a.obtain P, precipitation, and its rate, RINT = P / DT
CTE b.compute DF, incremental infiltration, from function EXPINF
CTE c.compute REX, infiltration excess runoff, REX = P - DF
CTE d.update P, water available for root zone, P = P - REX(same as P = DF)
CTE 0. Loop IA
CTE 1. compute ACF, ACF = 0.5 * (AC(IA) + AC(IA + 1))
CTE 2. reset UZ and EX(IA) to 0, UZ = 0 and EX(IA) = 0
CTE 3. compute SD(IA), local storage deficit, SD(IA) = SBAR + SZM * (TL - ST(IA)), and IF(SD(IA).LT.0.)SD(IA) = 0
CTE 4. compute SRZ(IA), root zone storage deficit, SRZ(IA) = SRZ(IA) - P
CTE 5. constrain SRZ(IA) to 0; if IF(SRZ(IA).LT.0.)THEN SUZ(IA) = SUZ(IA) - SRZ(IA) and SRZ(IA) = 0.
CTE 6. constrain SUZ(IA) to SD, unsaturated zone storage (not deficit); IF(SUZ(IA).GT.SD(IA)) THEN EX(IA) = SUZ(IA) - SD(IA) and SUZ(IA) = SD(IA)
CTE 6a.note this updated EX(IA), saturation excess flow (to outlet)
CTE 6b.note this updated SUZ(IA), unsaturated zone storage (to subsurface)
CTE 7. constrain SD(IA) to 0; IF(SD(IA).GT.0.)THEN UZ = SUZ(IA) / (SD(IA) * TD * DT)
CTE 8. constrain UZ to 0, unsaturated zone flow; IF(UZ.GT.SUZ(IA)) THEN UZ = SUZ(IA)
CTE 9. compute SUZ(IA) = SUZ(IA) - UZ, adding constraint IF(SUZ(IA).LT.0.0000001)THEN SUZ(IA) = 0
CTE 10. compute QUZ = QUZ + UZ * ACF, unsaturated zone runoff
CTE 11. compute EA = EP * (1 - SRZ(IA) / SRMAX), actual evapotranspiration
CTE 12. constrain EA to available water; IF(EA.GT.SRMAX - SRZ(IA)) THEN EA = SRMAX - SRZ(IA)
CTE 13. update SRZ(IA), SRZ(IA) = SRZ(IA) + EA and sum EA to SUMAE = SUMAE + EA * ACF
CTE 14. compute OF, saturation excess flow, only for IF(IA.GT.1)THEN IB = IA - 1 and IF(EX(IA).GT.0.)THEN OF = AC(IA) * (EX(IB) + EX(IA)) / 2
CTE 14a.ELSE IF(EX(IB).GT.0.)THEN ACF = ACF * EX(IB) / (EX(IB) - EX(IA)) and OF = ACF * EX(IB) / 2,
CTE 14c.compute QOF = QOF + OF, total overland flow
CTE 14c.compute ACM = ACM + ACF and IF(ACM.GT.ACMAX)ACMAX = ACM
CTE 15. End IA Loop
CTE 16. add REX, infiltration excess runoff, to QOF
CTE 17. compute QB, saturated zone drainage, QB = SZQ * EXP(-SBAR / SZM)
CTE 18. update SBAR, average storage deficit, getting bigger with QB flows, smaller with QUZ flows, SBAR = SBAR - QUZ + QB

CTE IT is time step, DT is time step increment
CTE P is rainfall available for infiltration after surface controls, updated P = P - REX
CTE RINT is rainfall rate, P / DT
CTE DF is volumetric increment of infiltration, determined with infiltration routine
CTE CUMF is cumulative infiltration, which is reset to 0 when ponding ends
CTE IROF is flag for ponding, IROF = 1 when ponded
CTE IHROF is hour(s) of ponding
CTE TP is time to ponding
CTE IA, the TI_Value number of 30 topographic index values, where IA = 1 is wettest topographic index value
CTE ACF is fractional area of ST(IA), initialized as TL = 0 then summed loop IA = 1 to NAC, ACF = 0.5 * (AC(IA) + AC(IA + 1));
CTE UZ is unsaturated flow to subsurface zone for specific SD(IA)
CTE EX(IA) is saturation excess runoff depth for each IA
CTE SD(IA) is local storage deficit, which constrains SUZ unsaturated zone storage
CTE SRZ(IA) when less than zero will contribute the lower level SUZ unsaturated zone storage(not a deficit)
CTE SUZ(IA) is unsaturated zone storage, which receives excess water from SRZ(IA), and generates saturation excess runoff if > SD(IA)
CTE REX, is infiltration excess runoff depth, REX = P - DF
CTE NAC is number of bins for topopgraphic index values, typically 30.
CTE ST(IA) is the TI_Value value, ST(IA) = LN(A / TANB)
CTE AC(IA) is fractional area of TI_Value value, either read in as a fraction or initially an area with AC(1) = AC(1) / tarea
CTE AC(IA) continued; AC(1) = 0 and AC(NAC + 1) = 0, which is needed to integrate the entire area and compute weighted average of ST
CTE TL is aerial weighted average ST, initialized as TL = 0 then summed loop J = 2 to NAC, TL = TL + AC(J) * (ST(J) + ST(J - 1)) / 2
CTE ACMAX is maximum area
CTE ACM is fractional area with saturated soils, and saturation excess
CTE P has been set to infiltrated water in this calculation above REX = P - DF, P = P - REX
CTE SRZ is root zone storage deficit, and it is reduced by the infiltrated P, and has maximum of SRMAX
CTE QUZ is total unsaturated zone flow to subsurface
CTE TD is number of time steps to delay the release of QUZ water from the SUZ, multiplied by the DT time step, in SUZ / (SD * TD * DT)
CTE EA is actual evapotranspiration from EP potential
CTE QOF (m/dt) is overland flow from saturated areas
CTE SBAR is average soil storage deficit
CTE SZM is scaling parameter m
CTE SZQ (m/dt) is maximum subsurface flow with all soil saturated
CTE SR0 IS INITIAL ROOT ZONE STORAGE DEFICIT BELOW FIELD CAPACITY; Q0 IS THE INITIAL DISCHARGE FOR THIS SUBCATCHMENT


C  TOPMODEL DEMONSTRATION PROGRAM VERSION 95.02
C
C
C  Compiled using Lahey Fortran77and Grafmatic Graphics
C
C  This version by Keith Beven 1985
C  Revised for distribution 1993, 1995
C
C * ***************************************************************
C  This program is distributed freely with only two conditions.
C
C  1. In any use for commercial or paid consultancy purposes a
C     suitable royalty agreement must be negotiated with Lancaster
C     University(Contact Keith Beven)
C
C  2. In any publication arising from use for research purposes the
C     source of the program should be properly acknowledgedand a
C     pre - print of the publication sent to Keith Beven at the address
C     below.
C
C  All rights retained 1993, 1995
C  Keith Beven
C  Centre for Research on Environmental Systemsand Statistics
C  Institute of Environmentaland Biological Sciences
C  Lancaster University, Lancaster LA1 4YQ, UK
C
C  Tel : (+44) 1524 593892  Fax : (+44) 1524 593985
C  Email : K.Beven@UK.AC.LANCASTER
C
C * ***************************************************************
C
C  SIMPLE SUBCATCHMENT VERSION OF TOPMODEL
C
C  This program allows single or multiple subcatchment calculations
C  but with single average rainfall and potential evapotranspiration
C  inputs to the whole catchment.Subcatchment discharges are routed
C  to the catchment outlet using a linear routing algorithm with
C  constant main channel velocityand internal subcatchment
C  routing velocity.The program requires ln(a / tanB) distributions
C  for each subcatchment.These may be calculated using the
C  GRIDATB program which requires raster elevation data as input.
C  It is recommended that those data should be 50 m resolution or
C  better.
C
C  NOTE that TOPMODEL is not intended to be a traditional model
C  package but is more a collection of concepts that can be used
C * ***where appropriate * ***.It is up to the user to verify that
C  the assumptions are appropriate(see discussion in
    C  Beven et al.(1994).This version of the model  will be
    C  best suited to catchments with shallow soils and moderate
    C  topography which do not suffer from excessively long dry
    C  periods.Ideally predicted contributing areas should be
    C  checked against what actually happens in the catchment.
    C
    C  It includes infiltration excess calculationsand parameters
    C  based on the exponential conductivity Green - Ampt model of
    C  Beven(HSJ, 1984) but if infiltration excess does occur it
    C  does so over whole area of a subcatchment.Spatial variability
    C  in conductivities can however be handled by specifying
    C  Ko parameter values for different subcatchments, even if they
    C  have the same ln(a / tanB) and routing parameters, ie.to
    C  represent different parts of the area.
    C
    C  Note that time step calculations are explicit ie.SBAR
    C  at start of time step is used to determine contributing area.
    C  Thus with long(daily) time steps contributing area depends on
    C  initial value together with any volume filling effect of daily
    C  inputs.Also baseflow at start of time step is used to update
    C  SBAR at end of time step
    C
    C  Current program limits are :
C          Number of time steps = 2500
C          Number of subcatchments = 10
C          Number of ln(a / tanB) increments = 30
C          Number of subcatchment routing ordinates = 10
C          Number of time delay histogram ordinates = 20
C          Size of subcatchment pixel maps = 100 x 100
C
C  Limits are mostly set in Common blocks in file TMCOMMON.FOR
C * ****************************************************************
C
C  This version uses five files as follows :
C       Channel 4 "TOPMOD.DAT" contains run and file information
C       Channel 7 < INPUTS$ > contains rainfall, pe and qobs data
C       Channel 8 < SUBCAT$ > contains subcatchment data
C       Channel 9 < PARAMS$ > contains parameter data
C       Channel 10 < OUTPUT$ > is output file

C  In addition
C       Channel 12 < MAPFILE$ > is used to read subcatchment ln(a / tanB)
C                   maps if IMAP = 1
C
C
C * ****************************************************************
C
INCLUDE TMCOMMON.FOR
CHARACTER * 15 INPUTS$, SUBCAT$, PARAMS$, OUTPUT$
OPEN(4, FILE = "TOPMOD.RUN", STATUS = "OLD")
READ(4, "(A)")TITLE
READ(4, "(A)")INPUTS$
READ(4, "(A)")SUBCAT$
READ(4, "(A)")PARAMS$
READ(4, "(A)")OUTPUT$
OPEN(7, FILE = INPUTS$, STATUS = "OLD")
OPEN(8, FILE = SUBCAT$, STATUS = "OLD")
OPEN(9, FILE = PARAMS$, STATUS = "OLD")
OPEN(10, FILE = OUTPUT$)
WRITE(10, 1001)TITLE
1001 FORMAT(1x, A)
Write(6, 1002)title
1002 Format(///1x,'TOPMODEL Version: TMOD95.02'////
    11x, 'This run :' / 1x, A//////////
    11x, 'Centre for Research on Environmental Systems and Statistics' /
    21x, 'Lancaster University, Lancaster LA1 4YQ, UK')
Write(6, 602)
602 format(/ 1x, '                       Press return to continue' / )
Read(5, *)
C
C  READ IN DT and RAINFALL, PE, QOBS INPUTS
CALL INPUTS
C
C  READ IN SUBCATCHMENT TOPOGRAPHIC DATA
READ(8, *)NSC, IMAP, IOUT
C
C  OPEN PARAMETER FILE

C  START LOOP ON SUBCATCHMENTS
DO 10 ISC = 1, NSC
If(iout.ge.2)Write(10, 600)ISC
600 Format(1x, 'Starting Subcatchment', I6)
C
C  INITIALISATION FOR THIS SUBCATCHMENT
CALL TREAD
CALL INIT
C
C  RUN MODEL FOR THIS SUBCATCHMENT INCLUDING LINEAR ROUTING CALCULATIONS
CALL TOPMOD

C
C  END LOOP ON SUBCATCHMENTS
C
10 CONTINUE
C  CALL RESULTS ROUTINE : if IRUN = 0 on return stop
CALL RESULTS
c  IRUN Disabled at present
CLOSE(5)
CLOSE(7)
CLOSE(8)
CLOSE(9)
CLOSE(10)
STOP
END
C
C
SUBROUTINE TOPMOD
C
INCLUDE TMCOMMON.FOR
DIMENSION EX(30)
C
C * ****************************************************************
C
C  THIS ROUTINE RUNS TOPMODEL FOR ONE SUBCATCHMENT, INCLUDING THE
C  LINEAR CHANNEL ROUTING CALCULATIONS.
C
C  The calculations are made for areal subdivisions based on the
C  NAC ln(a / tanB) subdivisions.The saturation deficit for each
C  subdivision is calculated from SBAR at the start of each time
C  step.
C
C  Each increment also has a root zone storage(SRZ) deficit which
C  is 0 at 'field capcacity' and becomes more positive as the soil
C  dries out;and an unsaturated zone storage(SUZ) which is zero at
C  field capacityand becomes more positive as storage increases.
C  SUZ has an upper limit of the local saturation deficit SD.
C  The local contributing area is where SD - SUZ is less than or
C  equal to zero.
C
C  REMEMBER SBAR, SD AND SRZ VALUES ARE DEFICITS; SUZ IS A STORAGE.
C
******************************************************************
C IROF is infiltration excess runoff
IROF = 0
REX = 0.
CUMF = 0.
ACMAX = 0.
SUMP = 0.
SUMAE = 0.
SUMQ = 0.
C
C  Initialise contributing area counts
IHROF = 0
do 5 ia = 1, nac
    5 ihour(ia) = 0
    C
    C  START LOOP ON TIME STEPS
    If(IOUT.ge.2)Write(10, 101)
    101 format(1x, '  it       p        ep       q(it)       quz',
        1'      q       sbar       qof')
    C
    DO 10 IT = 1, NSTEP
    QOF = 0.
    QUZ = 0.
    C
    EP = PE(IT)
    P = R(IT)
    SUMP = SUMP + P
    C
    C  SKIP INFILTRATION EXCESS CALCULATIONS IF INFEX = 0
    IF(INFEX.EQ.1) THEN


    C * ***************************************************************
    C  INFILTRATION EXCESS CALCULATIONS USING EXPINF ROUTINE BASED ON
    C  GREEN - AMPT INFILTRATION IN A SOIL WITH CONDUCTIVITY DECLINING
    C  EXPONENTIALLY WITH DEPTH(REF.BEVEN, HSJ, 1984)
    C
    C  NOTE THAT IF INFILTRATION EXCESS DOES OCCUR IT WILL DO SO OVER
    C  THE WHOLE SUBCATCHMENT BECAUSE OF HOMOGENEOUS SOIL ASSUMPTION
    C
    C  ALL PARAMETERS AND VARIABLES ON INPUT MUST BE IN M / H
    C
    C  THIS SECTION CAN BE OMITTED WITHOUT PROBLEM
    C * ***********************************************************8 * **
    IF(P.GT.0.)THEN
    C
    C  Adjust Rainfall rate from m / time step to m / h
    RINT = P / DT
    CALL EXPINF(IROF, IT, RINT, DF, CUMF)
    C  DF is volumetric increment of infiltration and is returned in m / DT
    REX = P - DF
    P = P - REX
    If(IROF.EQ.1)IHROF = IHROF + 1
    ELSE
    REX = 0.
    IROF = 0
    CUMF = 0.
    ENDIF
    ENDIF
    C * ***************************************************************
    C
    C P IS RAINFALL AVAILABLE FOR INFILTRATION AFTER SURFACE CONTROL
    C   CALCULATION
    C
    ACM = 0.
    C  START LOOP ON A / TANB INCREMENTS
    DO 30 IA = 1, NAC
    ACF = 0.5 * (AC(IA) + AC(IA + 1))
    UZ = 0.
    EX(IA) = 0.
    C
    C  CALCULATE LOCAL STORAGE DEFICIT
    SD(IA) = SBAR + SZM * (TL - ST(IA))
    IF(SD(IA).LT.0.)SD(IA) = 0.
    C
    CTE Note infiltration routine updates DF, which updates P in soil moisture routine
    C  ROOT ZONE CALCULATIONS
    CTE P is rainfall available for infiltration after surface controls, and becomes volumetric increment of infiltration
    CTE DF is volumetric increment of infiltration
    CTE P has been set to infiltrated water REX = P - DF, P = P - REX
    CTE SRZ is root zone storage deficit, and it is reduced by the infiltrated P
    CTE SRZ when less than zero will contribute the lower level SUZ unsaturated zone storage (not a deficit)
    SRZ(IA) = SRZ(IA) - P
    IF(SRZ(IA).LT.0.)THEN
    SUZ(IA) = SUZ(IA) - SRZ(IA)
    SRZ(IA) = 0.
    ENDIF
    C
    C  UZ CALCULATIONS
    IF(SUZ(IA).GT.SD(IA))THEN
    EX(IA) = SUZ(IA) - SD(IA)
    SUZ(IA) = SD(IA)
    ENDIF
    C
    C  CALCULATE DRAINAGE FROM SUZ
    IF(SD(IA).GT.0.)THEN
    UZ = SUZ(IA) / (SD(IA) * TD * DT)
    IF(UZ.GT.SUZ(IA))UZ = SUZ(IA)
    SUZ(IA) = SUZ(IA) - UZ
    IF(SUZ(IA).LT.0.0000001)SUZ(IA) = 0.
    QUZ = QUZ + UZ * ACF
    ENDIF
    C
    C * **************************************************************
    C  CALCULATE EVAPOTRANSPIRATION FROM ROOT ZONE DEFICIT
    C
    EA = 0.
    IF(EP.GT.0.)THEN
    EA = EP * (1 - SRZ(IA) / SRMAX)
    IF(EA.GT.SRMAX - SRZ(IA))EA = SRMAX - SRZ(IA)
    SRZ(IA) = SRZ(IA) + EA
    ENDIF
    SUMAE = SUMAE + EA * ACF
    SAE = SAE + EA * ACF
    C
    C * **************************************************************
    C
    C
    C  CALCULATION OF FLOW FROM FULLY SATURATED AREA
    C  This section assumes that a / tanB values are ordered from high to low
    C
    OF = 0.
    IF(IA.GT.1)THEN
    IB = IA - 1
    IF(EX(IA).GT.0.)THEN
    c  Both limits are saturated
    OF = AC(IA) * (EX(IB) + EX(IA)) / 2
    ACM = ACM + ACF
    ihour(ib) = ihour(ib) + 1
    ELSE
    c  Check if lower limit saturated(higher a / tanB value)
    IF(EX(IB).GT.0.)THEN
    ACF = ACF * EX(IB) / (EX(IB) - EX(IA))
    OF = ACF * EX(IB) / 2
    ACM = ACM + ACF
    ihour(ib) = ihour(ib) + 1
    ENDIF
    ENDIF
    ENDIF
    QOF = QOF + OF
    C
    C  Set contributing area plotting array
    CA(IT) = ACM
    IF(ACM.GT.ACMAX)ACMAX = ACM
    C
    C  END OF A / TANB LOOP
    30 CONTINUE

    C
    C  ADD INFILTRATION EXCESS
    QOF = QOF + REX

    IF(IROF.EQ.1)ACMAX = 1.
    C
    C  CALCULATE SATURATED ZONE DRAINAGE
    QB = SZQ * EXP(-SBAR / SZM)
    SBAR = SBAR - QUZ + QB
    QOUT = QB + QOF
    SUMQ = SUMQ + QOUT
    C
    C  CHANNEL ROUTING CALCULATIONS
    C  allow for time delay to catchment outlet ND as well as
    C  internal routing array
    DO 40 IR = 1, NR
    IN = IT + ND + IR - 1
    IF(IN.GT.NSTEP)GO TO 10
    Q(IN) = Q(IN) + QOUT * AR(IR)
    40 CONTINUE
    C
    If(IOUT.ge.2) write(10, 100)it, p, ep, q(it), quz, qb, sbar, qof
    100 format(1x, i4, 7e10.3)

    C  END OF TIME STEP LOOP
    10 CONTINUE
    C
    C  CALCULATE BALANCE TERMS
    SUMRZ = 0.
    SUMUZ = 0.
    DO 50 IA = 1, NAC
    ACF = 0.5 * (AC(IA) + AC(IA + 1))
    SUMRZ = SUMRZ + SRZ(IA) * ACF
    SUMUZ = SUMUZ + SUZ(IA) * ACF
    50 CONTINUE
    BAL = BAL + SBAR + SUMP - SUMAE - SUMQ + SUMRZ - SUMUZ
    Write(10, 650)SUBCAT, SUMP, SUMAE, SUMQ, SUMRZ, SUMUZ, SBAR, BAL
    WRITE(6, 650)SUBCAT, SUMP, SUMAE, SUMQ, SUMRZ, SUMUZ, SBAR, BAL
    650 FORMAT(1X, 'Water Balance for Subcatchment : ', A /
        11x, '   SUMP      SUMAE       SUMQ      SUMRZ   ',
        2   '  SUMUZ      SBAR        BAL' / 7e11.4)
    If(IOUT.ge.1)WRITE(10, 651)ACMAX
    651 FORMAT(1X, 'Maximum contributing area ', e12.5)
    RETURN
    END
    *
    C * **************************************************************
    *
    SUBROUTINE INPUTS

    *
    INCLUDE TMCOMMON.FOR
    *
    *This subroutine must read in rainfall, peand observed
    * discharges for T = 1, NSTEP with time step DT hours
    *
    READ(7, *)NSTEP, DT
    READ(7, *)(R(I), PE(I), QOBS(I), I = 1, NSTEP)
    CLOSE(7)
    DO 10 IT = 1, NSTEP
    10 Q(IT) = 0.
    RETURN
    END
    C
    C * *************************************************************
    C
    SUBROUTINE TREAD
    C
    INCLUDE TMCOMMON.FOR
    C
    READ(8, "(A)")subcat
    Write(10, 1010)subcat
    1010 Format(1x, 'Subcatchment : ', A)
    READ(8, *)NAC, AREA
    * NAC IS NUMBER OF A / TANB ORDINATES
    * AREA IS SUBCATCHMENT AREA AS PROPORTION OF TOTAL CATCHMENT
    READ(8, *)(AC(J), ST(J), J = 1, NAC)
    * AC IS DISTRIBUTION OF AREA WITH LN(A / TANB)
    * ST IS LN(A / TANB) VALUE
    tarea = ac(1)
    do 10 j = 2, nac
        tarea = tarea + ac(j)
        10  continue

        *
        * CALCULATE AREAL INTEGRAL OF LN(A / TANB)
        * NB.a / tanB values should be ordered from high to low with ST(1)
        * as an upper limit such that AC(1) should be zero, with AC(2) representing
        * the area between ST(1) and ST(2)
        TL = 0.
        AC(1) = AC(1) / tarea
        SUMAC = AC(1)
        DO 11 J = 2, NAC
        AC(J) = AC(J) / tarea
        SUMAC = SUMAC + AC(J)
        TL = TL + AC(J) * (ST(J) + ST(J - 1)) / 2
        11 CONTINUE
        AC(NAC + 1) = 0.
        *
        *READ CHANNEL NETWORK DATA
        READ(8, *)NCH
        READ(8, *)(ACH(J), D(J), J = 1, NCH)
        * ACH IS CUMULATIVE DISTRIBUTION OF AREA WITH DISTANCE D
        * FROM OUTLET.D(1) is distance from subcatchment outlet
        * ACH(1) = 0.
        *
        If(IOUT.ge.1)Write(10, 600)TL, SUMAC
        600 Format(1x, 'TL = ', f8.2, / 'SUMAC = ', f8.2)
        RETURN
        END
        *
        C * **************************************************************
        *
        *
        SUBROUTINE INIT
        DIMENSION TCH(10)
        INCLUDE TMCOMMON.FOR
        *
        *READ PARAMETER DATA
        READ(9, "(A)")SUBCAT
        READ(9, *)SZM, T0, TD, CHV, RV, SRMAX, Q0, SR0, INFEX, XK0, HF, DTH
        *
        * Convert parameters to m / time step DT
        * with exception of XK0 which must stay in m / h
        * Q0 is already in m / time step
        * T0 is input as Ln(To)
        RVDT = RV * DT
        CHVDT = CHV * DT
        T0DT = T0 + ALOG(DT)

        * Calculate SZQ parameter
        SZQ = EXP(T0DT - TL)
        * *
        *CONVERT DISTANCE / AREA FORM TO TIME DELAY HISTOGRAM ORDINATES
        *
        TCH(1) = D(1) / CHVDT
        DO 15 J = 2, NCH
        TCH(J) = TCH(1) + (D(J) - D(1)) / RVDT
        15 CONTINUE
        NR = INT(TCH(NCH))
        IF(FLOAT(NR).LT.TCH(NCH))NR = NR + 1
        ND = INT(TCH(1))
        NR = NR - ND
        DO 20 IR = 1, NR
        TIME = ND + IR
        IF(TIME.GT.TCH(NCH))THEN
        AR(IR) = 1.0
        ELSE
        DO 21 J = 2, NCH
        IF(TIME.LE.TCH(J))THEN
        AR(IR) = ACH(J - 1) + (ACH(J) - ACH(J - 1)) * (TIME - TCH(J - 1)) /
        1      (TCH(J) - TCH(J - 1))
        GOTO 20
        ENDIF
        21    CONTINUE
        ENDIF
        20 CONTINUE
        A1 = AR(1)
        SUMAR = AR(1)
        AR(1) = AR(1) * AREA
        IF(NR.GT.1)THEN
        DO 22 IR = 2, NR
        A2 = AR(IR)
        AR(IR) = A2 - A1
        A1 = A2
        SUMAR = SUMAR + AR(IR)
        AR(IR) = AR(IR) * AREA
        22    CONTINUE
        ENDIF
        If(IOUT.ge.1)write(10, 603)szq
        603 format(1x, 'SZQ  ', e12.5)
        If(IOUT.ge.1)WRITE(10, 604)TCH(NCH), SUMAR, (AR(IR), IR = 1, NR)
        604 FORMAT(1X, 'SUBCATCHMENT ROUTING DATA' /
            1  1X, 'Maximum Routing Delay  ', E12.5 /
            2  1X, 'Sum of histogram ordinates  ', f10.4 / (1X, 5E12.5))
        *
        *INITIALISE SRZ AND Q0 VALUES HERE
        * SR0 IS INITIAL ROOT ZONE STORAGE DEFICIT BELOW FIELD CAPACITY
        * Q0 IS THE INITIAL DISCHARGE FOR THIS SUBCATCHMENT
        *
        *INITIALISE STORES
        DO 25 IA = 1, NAC
        SUZ(IA) = 0.
        25 SRZ(IA) = SR0
        SBAR = -SZM * ALOG(Q0 / SZQ)
        c
        c  Reinitialise discharge array
        SUM = 0.
        DO 29 I = 1, ND
        29 Q(I) = Q(I) + Q0 * AREA
        DO 30 I = 1, NR
        SUM = SUM + AR(I)
        IN = ND + I
        30 Q(IN) = Q(IN) + Q0 * (AREA - SUM)
        *
        *Initialise water balance.BAL is positive for storage
        BAL = -SBAR - SR0
        If(IOUT.ge.1)Write(10, 605)BAL, SBAR, SR0
        605 Format(1x, 'Initial Balance BAL ', e12.5 /
            1       1x, 'Initial SBAR        ', e12.5 /
            2       1x, 'Initial SR0         ', e12.5)
        *
        RETURN
        END
        C
        C * *************************************************************
        C
        SUBROUTINE EXPINF(IROF, IT, RINT, DF, CUMF)
        C
        INCLUDE TMCOMMON.FOR
        DOUBLE PRECISION CONST, SUM, FC, FUNC, CD, SZF, XKF
        DATA E / 0.00001 /
        C * ************************************************************
        C
        C  SUBROUTINE TO CALCULATE INFILTRATION EXCESS RUNOFF USING THE
        C  EXPONENTIAL GREEN - AMPT MODEL.
        C
        C * *************************************************************
        C
        C
        C  Note that HF and DTH only appear in product CD
        CD = HF * DTH
        SZF = 1. / SZM
        XKF = XK0
        IF(IROF.EQ.1)GO TO 10
        C  PONDING HAS ALREADY OCCURRED - GO TO EXCESS CALCULATION
        C
        IF(CUMF.EQ.0.)GOTO 7
        C  FIRST TIME STEP, OVERFLOW IF CUMF = 0, GO DIRECT TO F2 CALCULATION
        C  INITIAL ESTIMATE OF TIME TO PONDING
        F1 = CUMF
        R2 = -XKF * SZF * (CD + F1) / (1 - EXP(SZF * F1))
        IF(R2.LT.RINT)THEN
        C  PONDING STARTS AT BEGINNING OF TIME STEP
        TP = (IT - 1.) * DT
        IROF = 1
        F = CUMF
        GO TO 8
        ENDIF
        7 F2 = CUMF + DT * RINT
        IF(F2.EQ.0.)GO TO 20
        R2 = -XKF * SZF * (CD + F2) / (1 - EXP(SZF * F2))
        IF(R2.GT.RINT)GO TO 20
        F = CUMF + R2 * DT
        DO 9 I = 1, 20
        R2 = -XKF * SZF * (CD + F) / (1 - EXP(SZF * F))
        IF(R2.GT.RINT)THEN
        F1 = F
        F = (F2 + F) * 0.5
        IF(ABS(F - F1).LT.E)GO TO 11
        ELSE
        F2 = F
        F = (F1 + F) * 0.5
        IF(ABS(F - F2).LT.E)GO TO 11
        ENDIF
        9 CONTINUE
        WRITE(6, 600)
        600 FORMAT(1X, 'MAXIMUM NO OF ITERATIONS EXCEEDED')
        11 CONTINUE
        TP = (IT - 1) * DT + (F - CUMF) / RINT
        IF(TP.GT.IT * DT)GO TO 20
        C
        C  SET UP DEFINITE INTEGRAL CONSTANT USING FP
        C
        8 CONST = 0
        FAC = 1
        FC = (F + CD)
        DO 12 J = 1, 10
        FAC = FAC * J
        ADD = (FC * SZF) * *J / (J * FAC)
        CONST = CONST + ADD
        12 CONTINUE
        CONST = DLOG(FC) - (DLOG(FC) + CONST) / DEXP(SZF * CD)
        C
        IROF = 1
        F = F + 0.5 * RINT * (IT * DT - TP)
        10 CONTINUE
        C
        C  NEWTON - RAPHSON SOLUTION FOR F(T)
        DO 14 I = 1, 20
        C
        C  CALCULATE SUM OF SERIES TERMS
        FC = (F + CD)
        SUM = 0.
        FAC = 1.
        DO 13 J = 1, 10
        FAC = FAC * J
        ADD = (FC * SZF) * *J / (J * FAC)
        SUM = SUM + ADD
        13 CONTINUE
        FUNC = -(DLOG(FC) - (DLOG(FC) + SUM) / DEXP(SZF * CD) - CONST) / (XKF * SZF)
        1 - (IT * DT - TP)
        DFUNC = (EXP(SZF * F) - 1) / (XKF * SZF * FC)
        DF = -FUNC / DFUNC
        F = F + DF
        IF(ABS(DF).LE.E)GO TO 15
        14 CONTINUE
        WRITE(6, 600)
        15 CONTINUE
        IF(F.LT.CUMF + RINT)THEN
        DF = F - CUMF
        CUMF = F
        C  SET UP INITIAL ESTIMATE FOR NEXT TIME STEP
        F = F + DF
        RETURN
        ENDIF
        20 CONTINUE
        C  THERE IS NO PONDING IN THIS TIME STEP
        IROF = 0
        DF = RINT * DT
        CUMF = CUMF + DF
        RETURN
        END
        C
        C
        C * **************************************************************
        C
        SUBROUTINE RESULTS
        INCLUDE TMCOMMON.FOR
        C
        C  OBJECTIVE FUNCTION CALCULATIONS
        F1 = 0.
        F2 = 0.
        SUMQ = 0.
        SSQ = 0.
        DO 60 IT = 1, NSTEP
        SUMQ = SUMQ + QOBS(IT)
        SSQ = SSQ + QOBS(IT) * QOBS(IT)
        F1 = F1 + (Q(IT) - QOBS(IT)) * *2
        F2 = F2 + ABS(Q(IT) - QOBS(IT))
        60 CONTINUE
        QBAR = SUMQ / NSTEP
        VARQ = (SSQ / NSTEP - QBAR * QBAR)
        VARE = F1 / NSTEP
        E = 1 - VARE / VARQ
        C
        C  add objective function values to output file
        write(6, 621)f1, e, f2, qbar, varq, vare
        write(10, 621)f1, e, f2, qbar, varq, vare
        621 format(//1x,'Objective function values'/
            1 1x, 'F1 ', e12.5, '   E ', f12.5, '   F2 'e12.5//
            2 1x, 'Mean Obs Q ', e12.5, '    Variance Obs Q ', e12.5 /
            3 '     Error Variance', e12.5)
        C
        C
        C
        RETURN
        END
        C
        C * *************************************************************
        C

*/
