      PROGRAM MAIN
C     MAIN PROGRAM OF WATERSHED RESOURCES MANAGEMENT MODEL - A
C     CONTINUOUS SIMULATION OF AGRICULTURAL WATERSHED HYDROLOGY.
C     THIS PROGRAM CALLS THE FOLLOWING 1ST ORDER MODULES:
C     INITIALIZATION ROUTINE, TIMING ROUTINE, RAINFALL EVENT
C     ROUTINE, EVAPOTRANSPIRATION EVENT ROUTINE, AND AN OPTIONAL
C     REPORT GENERATOR .........................................
      IMPLICIT REAL(K)
C     DIRECT PROGRAM INPUT: NUMBER OF SOIL TYPES AS ARRAY SIZE IN
C     THE FOLLOWING SOIL PARAMETERS ..........................
      REAL BDENS(3),ORGMA(3),CLAY(3),SAND(3),SILT(3),WTCF(3),
     *CEC(3),DWBA(3),PSIZE(3),FC(3),WP(3),SAT(3),STAB(3),
     *WDCLAY(3),MG(3),FE(3),AL(3),COND(3),CA(3),NA(3),            
     *CECNH4(3),VFSAND(3),SPSURF(3),CACO3(3),SAR(3), 
     *EFFPOR(3),PORFAC(3),TOTPOR(3),CFCORR(3),EACORR(3),VOLCF(3),
     *WFCP(3),KONDS(3),RSWR(3),WISCHM(3)  
C     DIRECT PROGRAM INPUT: NUMBER OF VEGETATION TYPES AS ARRAY   
C     SIZE IN THE FOLLOWING VEGETATION PARAMETERS ..............
      REAL CISCAP(8),KC(8),LAI(8)
C     DIRECT PROGRAM INPUT: NUMBER OF LOCATIONS IN SPACIAL        
C     DISTRIBUTION, AS ARRAY SIZE IN THE FOLLOWING SYSTEM AND
C     WATERSHED PARAMETERS (WITH THE EXCEPTION OF POLY, LCULV, 
C     NCULV, LSPILL AND ETK; TO BE LEFT AS IS) ..................
      REAL CANFAC(73),BACAN(73),CAN(73),OPEN(73),BAOPEN(73), 
     *COVER(73),GI(73),DELTX2(73),MANN1(73),MANN2(73),
     *PCAN(73),CANHT(73),GRCE(73),PRESD(73),KTC(73),   
     *DW(73),KSTO(73),DSUB(73),DSURF(73),MANN3(73),ASW(73), 
     *KRILL(73),TAWC(73),KRILLY(73),FSFAC(73),FSANG(73),   
     *CANEF(73),CWIDTH(73),POLY(4),FRI(73),RETAIN(73), 
     *GIM1(73),GSC1(73),KOND(73),KONDC(73),KONDO(73),  
     *KONDE(73),PFC(73),SCCOR(73),FGCORR(73),LCULV,NCULV,NOUT,
     *LSPILL,FSLOPE(73),CSLOPE(73)
      INTEGER SOIL(73),VEG(73),DUAL(73),STRUC1(73),
     *STRUC2(73),STRUC3(73),ETK,NELEMV(73),NELEMC(73),
     *NELEMT(73),NELEMW(73),TEST(73)
C     DIRECT PROGRAM INPUT: NUMBER OF LOCATIONS IN SPACIAL        
C     DISTRIBUTION AS ARRAY SIZE IN THE FOLLOWING SYSTEM AND      
C     STATE VARIABLES, WITH THE EXCEPTION OF STIMES WITH SIZE EQUAL
C     TO NUMBER OF EVENTS - RAINFALL & ET PEROIDS CONSECUTIVELY 
C     INTERSPACED (STARTING WITH RAINFALL) - PLUS 1. 
      REAL CIS(73),FILT(73),CFILT(73), 
     *FLOWV(73),FLOWC(73),FLOWT(73),FLOWW(73),QVLAT(73), 
     *QCLAT(73),QTLAT(73),QWLAT(73),QGUP(73),QVUP(73),
     *QCUP(73),QTUP(73),QWUP(73),SEDLDV(73),SEDLDC(73), 
     *SEDLDT(73),SEDLDW(73),SUBS(73),SURFS(73),QSUB(73), 
     *RFLOW(73),RSTORE(73),GVUP(73),GCUP(73),GTUP(73),   
     *GWUP(73),RCOMP(73),POND(73),STIMES(609)

C	Newly Introduced variables 
	  REAL AREA1(73),AREA2(73),RK(73),RB(73),CAF,W,Z(73),Lm,LATRAD
C     NOTE: THE ABOVE ARRAY SIZE SPECIFICATION SHOULD BE APPLIED TO
C     THE APPROPRIATE PARAMETERS & VARIABLES IN THE FOLLOWING     
C     PROGRAM MODULES: INITIA, TIMER, RAINFL, EVAPT AND BASEFL.   
C     ALSO TO BE APPLIED IN DATA STATEMENTS IN EVAPT.
	  
      CHARACTER EVENT1*20,NAME*15,TITLE*70  
	NAME='output.dat'   
	CALL DELOUT()	
      OPEN (1,FILE='FOR001.DAT',STATUS='OLD')
	OPEN (2,FILE='NEWINP.DAT',STATUS='OLD')
      OPEN (3,FILE='FOR003.DAT',STATUS='OLD')
      OPEN (5,FILE='FOR005.DAT',STATUS='OLD')	
      OPEN (6,FILE=NAME,STATUS='NEW',POSITION='APPEND')
	OPEN (7,FILE='z.dat',STATUS='OLD')
C	WRITE (*,*) 'Enter catchment type: '
C	WRITE (*,*) '1 for forested and 2 for grassed: '
C	READ (*,*) CATCH
C     WATERSHED NAME AND LOCATION INPUT, AND OUTPUT FILE 
C     HEADINGS...........................................
      READ (1,*) TITLE
      WRITE (*,FMT='(14X,"CONTINUOUS SIMULATION OF WATERSHED", 
     *" HYDROLOGY",/,5X,A64,/,5X,"TIME(mins)   RUNOFF", 
     *" RATE(L/s)  SED. LOAD(kg/s)",
     *"   SED. CONC.(ppm)")') TITLE	 
C     INITIALIZATION OF SYSTEM AND STATE VARIABLES, COMPUTATION AND
C     ASSIGNMENT OF PARAMETER DATA ............................... 
      CALL INITIA(NSOIL,BDENS,ORGMA,CLAY,SAND,SILT,WTCF,CEC, 
     *DWBA,PSIZE,FC,WP,SAT,STAB,WDCLAY,MG,FE,AL,COND,CA,NA,CECNH4,
     *VFSAND,SPSURF,CACO3,TOTPOR,VOLCF,CFCORR,EACORR,EFFPOR,WFCP,
     *RSWR,KONDS,PORFAC,WISCHM,NVEG,CISCAP,KC,LAI,NGRID,NSEG,KP,
     *BETA,DELTT,RNWL,RMWL,TSLOPE,WSLOPE,BFLOW,POND,
     *MONITO,NOUT,POLY,NEVENT,PRIN,ENDS,KRET,GIM,GSC,SCULV,XAREA,
     *LCULV,EII,EOI,DCULV,KLOSSE,KLOSSC,KORIF,KSPILL,LSPILL,KB,
     *DELTX1,NCULV,TWIDTH,WWIDTH,BACAN,CAN,OPEN,BAOPEN,COVER,GI,
     *PFC,SCCOR,ASW,DELTX2,MANN1,MANN2,MANN3,FSLOPE,CWIDTH,CSLOPE,
     *CANHT,PRESD,KTC,DW,DSUB,DSURF,STRUC1,STRUC2,STRUC3,DUAL,
     *NELEMV,NELEMC,NELEMT,NELEMW,SOIL,VEG,GIM1,GSC1,FGCORR,KOND,
     *CANFAC,KONDC,KONDO,KONDE,KRILL,TAWC,KRILLY,FSANG,FSFAC,
     *PCAN,CANEF,GRCE,KSTO,FRI,RETAIN,MOUT,RI,ETK,CIS,CFILT,FILT,
     *FLOWV,FLOWC,FLOWT,FLOWW,QVLAT,QCLAT,QTLAT,QWLAT,RFLOW,QSUB,
     *RSTORE,RCOMP,QVUP,QCUP,QTUP,QWUP,QGUP,SEDLDV,SEDLDC,SEDLDT,
     *SEDLDW,GVUP,GCUP,GTUP,GWUP,SUBS,SURFS,STIME,NEXTE,TEST,
     *STIMES,SAR,Z,CAF,W,AREA1,AREA2,RK,RB,Lm,LATRAD)	
  1   IF(STIMES(NEXTE) .LT. ENDS)THEN
C       TIMING ROUTINE DETERMINES EVENT TYPE AND PROVIDES EVENT   
C       TIME ADVANCE ..........................................
        CALL TIMER(NEXTE,STIMES,DELTT,NSIM,ETT,EVENT1)
        IF(EVENT1 .EQ. 'RAINFALL')THEN
C         RAINFALL EVENT ROUTINE UPDATES SYSTEM STATE .........
          CALL RAINFL(NSIM,STIME,DELTT,NGRID,VEG,SOIL,CIS,CISCAP,
     *    CFILT,FILT,KONDE,EFFPOR,ASW,WFCP,DUAL,QVUP,DELTX1,
     *    QGUP,DSUB,SAT,WP,SUBS,RAIN,RETAIN,FLOWV,QVLAT,MANN1,
     *    FSLOPE,BETA,FLOWC,QCLAT,FLOWW,DELTX2,MANN2,CSLOPE,QCUP,
     *    CWIDTH,DWBA,DW,KTC,KRILL,TAWC,GVUP,KRILLY,FSFAC,CANEF,
     *    GRCE,SEDLDV,GCUP,PSIZE,SEDLDC,DSURF,KSTO,QSUB,STRUC1,
     *    STRUC2,STRUC3,RFLOW,RCOMP,RSTORE,POLY,RNWL,RMWL,KSPILL,
     *    LSPILL,DCULV,EII,EOI,XAREA,KLOSSE,KLOSSC,LCULV,SCULV,   
     *    KORIF,NCULV,FLOWT,QTLAT,MANN3,TSLOPE,QTUP,TWIDTH,GTUP,RI,
     *    QWLAT,WSLOPE,QWUP,WWIDTH,GWUP,MOUT,MONITO,NELEMV,NELEMC,
     *    NELEMT,NELEMW,SEDLDT,SEDLDW,SURFS,NAME,POND,TEST,
     *    BFLOW,NOUT,FC,NSEG,FLOW,AREA1,AREA2,RK,RB,CAF)
        ELSE
C         ET EVENT ROUTINE UPDATES SYSTEM STATE ............... 
          CALL EVAPT(STIME,DELTT,ETK,NGRID,VEG,SOIL,FILT,CFILT,
     *    KONDE,WFCP,EFFPOR,ASW,DUAL,DSUB,SAT,WP,SUBS,
     *    QVUP,FLOWV,QVLAT,MANN1,MANN2,FSLOPE,BETA,FLOWC,QCLAT,
     *    FLOWW,DELTX1,DELTX2,CSLOPE,QCUP,CWIDTH,DWBA,DW,KTC,KRILL,
     *    TAWC,GVUP,KRILLY,FSFAC,CANEF,GRCE,SEDLDV,GCUP,PSIZE,
     *    SEDLDC,DSURF,FC,QGUP,KSTO,QSUB,STRUC1,STRUC2,STRUC3,
     *    RFLOW,RCOMP,RSTORE,POLY,RNWL,RMWL,KSPILL,LSPILL,DCULV,
     *    EII,EOI,XAREA,KLOSSE,KLOSSC,LCULV,SCULV,KORIF,NCULV,
     *    FLOWT,QTLAT,MANN3,TWIDTH,WWIDTH,TSLOPE,WSLOPE,QTUP,GTUP, 
     *    RI,QWLAT,QWUP,GWUP,GI,PCAN,LAI,KC,KP,KB,CIS,MOUT,MONITO, 
     *    NELEMV,NELEMC,NELEMT,NELEMW,SEDLDT,SEDLDW,ETT,SURFS,
     *    TOTPOR,NOUT,BFLOW,CLAY,VFSAND,CACO3,SAR,
     *    SPSURF,WDCLAY,SAND,NAME,POND,TEST,RETAIN,NSEG,FLOW,
     *	W,Z,Lm,LATRAD,CAF)     
        ENDIF
      GO TO 1
      ENDIF 
C     SIMULATION OUTPUT FILE PROCESSING ........................
C     CALL REPORT(                                               )
      END
C     ............................................................
C     SUBROUTINE CALLED IN MAIN PROGRAM (1ST ORDER ROUTINE): 
C     INITIALIZATION ROUTINE ...................................
      SUBROUTINE INITIA(NSOIL,BDENS,ORGMA,CLAY,SAND,SILT,WTCF,CEC, 
     *DWBA,PSIZE,FC,WP,SAT,STAB,WDCLAY,MG,FE,AL,COND,CA,NA,CECNH4,
     *VFSAND,SPSURF,CACO3,TOTPOR,VOLCF,CFCORR,EACORR,EFFPOR,WFCP,
     *RSWR,KONDS,PORFAC,WISCHM,NVEG,CISCAP,KC,LAI,NGRID,NSEG,KP,
     *BETA,DELTT,RNWL,RMWL,TSLOPE,WSLOPE,BFLOW,POND,
     *MONITO,NOUT,POLY,NEVENT,PRIN,ENDS,KRET,GIM,GSC,SCULV,XAREA,
     *LCULV,EII,EOI,DCULV,KLOSSE,KLOSSC,KORIF,KSPILL,LSPILL,KB,
     *DELTX1,NCULV,TWIDTH,WWIDTH,BACAN,CAN,OPEN,BAOPEN,COVER,GI,
     *PFC,SCCOR,ASW,DELTX2,MANN1,MANN2,MANN3,FSLOPE,CWIDTH,CSLOPE,
     *CANHT,PRESD,KTC,DW,DSUB,DSURF,STRUC1,STRUC2,STRUC3,DUAL,
     *NELEMV,NELEMC,NELEMT,NELEMW,SOIL,VEG,GIM1,GSC1,FGCORR,KOND,
     *CANFAC,KONDC,KONDO,KONDE,KRILL,TAWC,KRILLY,FSANG,FSFAC,
     *PCAN,CANEF,GRCE,KSTO,FRI,RETAIN,MOUT,RI,ETK,CIS,CFILT,FILT,
     *FLOWV,FLOWC,FLOWT,FLOWW,QVLAT,QCLAT,QTLAT,QWLAT,RFLOW,QSUB,
     *RSTORE,RCOMP,QVUP,QCUP,QTUP,QWUP,QGUP,SEDLDV,SEDLDC,SEDLDT,
     *SEDLDW,GVUP,GCUP,GTUP,GWUP,SUBS,SURFS,STIME,NEXTE,TEST,
     *STIMES,SAR,Z,CAF,W,AREA1,AREA2,RK,RB,Lm,LATRAD)
	 
	 !DEC$ ATTRIBUTES DLLEXPORT, STDCALL :: INITIA
	 !DEC$ ATTRIBUTES ALIAS:'INITIA' :: INITIA
     !DEC$ ATTRIBUTES REFERENCE :: NSOIL,BDENS,ORGMA,CLAY,SAND,SILT,WTCF,CEC, 
     *DWBA,PSIZE,FC,WP,SAT,STAB,WDCLAY,MG,FE,AL,COND,CA,NA,CECNH4,
     *VFSAND,SPSURF,CACO3,TOTPOR,VOLCF,CFCORR,EACORR,EFFPOR,WFCP,
     *RSWR,KONDS,PORFAC,WISCHM,NVEG,CISCAP,KC,LAI,NGRID,NSEG,KP,
     *BETA,DELTT,RNWL,RMWL,TSLOPE,WSLOPE,BFLOW,POND,
     *MONITO,NOUT,POLY,NEVENT,PRIN,ENDS,KRET,GIM,GSC,SCULV,XAREA,
     *LCULV,EII,EOI,DCULV,KLOSSE,KLOSSC,KORIF,KSPILL,LSPILL,KB,
     *DELTX1,NCULV,TWIDTH,WWIDTH,BACAN,CAN,OPEN,BAOPEN,COVER,GI,
     *PFC,SCCOR,ASW,DELTX2,MANN1,MANN2,MANN3,FSLOPE,CWIDTH,CSLOPE,
     *CANHT,PRESD,KTC,DW,DSUB,DSURF,STRUC1,STRUC2,STRUC3,DUAL,
     *NELEMV,NELEMC,NELEMT,NELEMW,SOIL,VEG,GIM1,GSC1,FGCORR,KOND,
     *CANFAC,KONDC,KONDO,KONDE,KRILL,TAWC,KRILLY,FSANG,FSFAC,
     *PCAN,CANEF,GRCE,KSTO,FRI,RETAIN,MOUT,RI,ETK,CIS,CFILT,FILT,
     *FLOWV,FLOWC,FLOWT,FLOWW,QVLAT,QCLAT,QTLAT,QWLAT,RFLOW,QSUB,
     *RSTORE,RCOMP,QVUP,QCUP,QTUP,QWUP,QGUP,SEDLDV,SEDLDC,SEDLDT,
     *SEDLDW,GVUP,GCUP,GTUP,GWUP,SUBS,SURFS,STIME,NEXTE,TEST,
     *STIMES,SAR,Z,CAF,W,AREA1,AREA2,RK,RB,Lm,LATRAD 
	 
	 
	 
	 
      IMPLICIT REAL(K)
C     SOIL TYPE PARAMETERS ......................................
      REAL BDENS(3),ORGMA(3),CLAY(3),SAND(3),SILT(3),WTCF(3),
     *CEC(3),DWBA(3),PSIZE(3),FC(3),WP(3),SAT(3),STAB(3),
     *WDCLAY(3),MG(3),FE(3),AL(3),COND(3),CA(3),NA(3),            
     *CECNH4(3),VFSAND(3),SPSURF(3),CACO3(3),SAR(3), 
     *EFFPOR(3),PORFAC(3),TOTPOR(3),CFCORR(3),
     *EACORR(3),VOLCF(3),WFCP(3),KONDS(3),RSWR(3),WISCHM(3)       
C     VEGETATION TYPE PARAMETERS .................................
      REAL CISCAP(8),KC(8),LAI(8)
C     SYSTEM AND WATERSHED PARAMETERS ............................ 
      REAL CANFAC(73),BACAN(73),CAN(73),OPEN(73),BAOPEN(73), 
     *COVER(73),GI(73),DELTX2(73),MANN1(73),MANN2(73),
     *PCAN(73),CANHT(73),GRCE(73),PRESD(73),KTC(73),   
     *DW(73),KSTO(73),DSUB(73),DSURF(73),MANN3(73),ASW(73), 
     *KRILL(73),TAWC(73),KRILLY(73),FSFAC(73),FSANG(73),   
     *CANEF(73),CWIDTH(73),POLY(4),FRI(73),RETAIN(73), 
     *GIM1(73),GSC1(73),KOND(73),KONDC(73),KONDO(73),  
     *KONDE(73),PFC(73),SCCOR(73),FGCORR(73),LCULV,NCULV,NOUT,    
     *LSPILL,FSLOPE(73),CSLOPE(73),W,AREA1(73),AREA2(73),Z(73),
     *RK(73),RB(73),Lm,LATRAD
      INTEGER SOIL(73),VEG(73),DUAL(73),STRUC1(73),
     *STRUC2(73),STRUC3(73),ETK,NELEMV(73),NELEMC(73),
     *NELEMT(73),NELEMW(73),TEST(73)      
C     SYSTEM AND STATE VARIABLES.................................
      REAL CIS(73),FILT(73),CFILT(73), 
     *FLOWV(73),FLOWC(73),FLOWT(73),FLOWW(73),QVLAT(73), 
     *QCLAT(73),QTLAT(73),QWLAT(73),QGUP(73),QVUP(73),
     *QCUP(73),QTUP(73),QWUP(73),SEDLDV(73),SEDLDC(73), 
     *SEDLDT(73),SEDLDW(73),SUBS(73),SURFS(73),QSUB(73), 
     *RFLOW(73),RSTORE(73),GVUP(73),GCUP(73),GTUP(73),   
     *GWUP(73),RCOMP(73),POND(73),STIMES(609)
C     INPUT: SOILS DATA ................................ 
C     NOTE - FOLLOWING UNITS ARE REQUIRED: BDENS, g/cc;           
C     ORGMA/CLAY/SAND/SILT/VFSAND/WTCF/WDCLAY/FE/AL/CACO3, %;
C     STAB/FC/WP/SAT, %; CA/MG/NA/CECNH4, meq/100g; CEC/SAR,    
C     fraction; SPSURF, mgEGME/g; PSIZE, mm; DWBA, fraction; COND, 
C     mmhos	 
	W = 2.84
	  
      READ (1,*)
      READ (1,*) NSOIL	
      DO 2 L=1,NSOIL
         READ (1,*)
         READ (1,*) BDENS(L),ORGMA(L),CLAY(L),SAND(L),SILT(L) 
         READ (1,*) WTCF(L),CEC(L),DWBA(L),PSIZE(L),FC(L)
         READ (1,*) WP(L),SAT(L),STAB(L),WDCLAY(L),MG(L)
         READ (1,*) FE(L),AL(L),COND(L),CA(L),NA(L)               
         READ (1,*) CECNH4(L),VFSAND(L),SPSURF(L),CACO3(L),SAR(L) 
  2   CONTINUE
C     SOILS DATA BY COMPUTATION ..............................
C     NOTE - UNITS: TOTPOR/CFCORR/EACORR/EFFPOR/RSWR/PORFAC,      
C     fraction (vol.); VOLCF, %; WFCP, cm; KONDS, cm/hr; 
C     WISCHM, none
      DO 3 L=1,NSOIL
         TOTPOR(L)=(2.65-BDENS(L))/2.65
         VOLCF(L)=ABS((WTCF(L)/2.65)*100.)/(100.-                 
     *   (WTCF(L)/BDENS(L))+(WTCF(L)/2.65))
         CFCORR(L)=(100.-VOLCF(L))/100.
         EACORR(L)=ABS(1.-(3.8+0.19*(CLAY(L)**2.)/1000.-          
     *   0.337*SAND(L)+0.126*CEC(L)*CLAY(L)+
     *   ORGMA(L)*(SAND(L)/200.)**2.)/100.)
         EFFPOR(L)=TOTPOR(L)*CFCORR(L)*EACORR(L)
         X1=5.34+0.185*((CLAY(L)/100.)**2.)-2.484*TOTPOR(L)    
         X2=-0.0021*((CLAY(L)/100.)**2.)-
     *       0.0436*(SAND(L)/100.)*TOTPOR(L)       
         X3=-0.6175*(CLAY(L)/100.)*TOTPOR(L)
         X4=0.00144*((SAND(L)/100.)**2.)*(TOTPOR(L)**2.)          
         X5=-0.00855*((CLAY(L)/100.)**2.)*(TOTPOR(L)**2.)         
         X6=-0.000013*(CLAY(L)/100.)*((SAND(L)/100.)**2.)
         X7=0.009*TOTPOR(L)*((CLAY(L)/100.)**2.)                  
         X8=-0.00073*TOTPOR(L)*((SAND(L)/100.)**2.)
         X9=0.000005*(SAND(L)/100.)*((CLAY(L)/100.)**2.)
         X0=0.5003*(CLAY(L)/100.)*(TOTPOR(L)**2.)
         XTEMP=X1+X2+X3+X4+X5+X6+X7+X8+X9+X0
         WFCP(L)=EXP(XTEMP)
         RSWR(L)=0.2+0.1*ORGMA(L)+0.25*CLAY(L)*(CEC(L)**0.45)*    
     *   (BDENS(L)/100.)*EACORR(L)*CFCORR(L)
         Y1=-0.17+0.181*CLAY(L)                          
         Y2=-0.00000069*(SAND(L)**2.)*(CLAY(L)**2.)              
         Y3=-0.00000041*(SAND(L)**2.)*(SILT(L)**2.)
         Y4=0.000118*(SAND(L)**2.)*(BDENS(L)**2.)
         Y5=0.00069*(CLAY(L)**2.)*(BDENS(L)**2.)
         Y6=0.000049*(SAND(L)**2.)*CLAY(L)
         Y7=-0.000085*SILT(L)*(CLAY(L)**2.)
         CONST=Y1+Y2+Y3+Y4+Y5+Y6+Y7
         ABC=((EFFPOR(L)-RSWR(L))**3.)/((1.-EFFPOR(L))**2.)
         IF(RSWR(L) .EQ. 0.)THEN
           XYZ=0.
         ELSE
           XYZ=(BDENS(L)**2.)/RSWR(L)*0.35*(CONST**2.)/1000.
         ENDIF
         KONDS(L)=ABC*XYZ
         YTEMP=(0.96-0.032*SAND(L)+0.04*CLAY(L)-0.032*BDENS(L))
         PORFAC(L)=EXP(YTEMP)
         WISCHM(L)=(SILT(L)+VFSAND(L))*(100.-CLAY(L))
  3   CONTINUE 
C     INPUT: VEGETATION DATA .................................
C     NOTE - REQUIRED UNITS: CISCAP, mm; KC/LAI, none   
      READ (1,*)
      READ (1,*) NVEG
      DO 4 L=1,NVEG
         READ (1,*)
         READ (1,*) CISCAP(L),KC(L),LAI(L) 
  4   CONTINUE 
C     INPUT: WATERSHED DATA, SYSTEM PARAMETERS, CONSERVATION AND  
C     HYDRAULIC STRUCTURES' DATA ................................
C     NOTE (1): MULTIPLE HYDRAULIC STRUCTURES OF EACH TYPE -      
C     DETENTION POND OR CULVERT - ARE ASSUMED UNIFORM AND TREATED 
C     TO THE SAME DATA SET. ALSO, VEGETAL RETARDANCE INDEX IS     
C     ASSUMED UNIFORM FOR GRASS WATERWAYS, AS IS SLOPE FOR TERRACES
C     OR WATERWAYS. 
C     NOTE (2) - UNITS: DSURF/DSUB, cm; PRIN/DELTT/ENDS, mins;    
C     KB/KP/ALPHA/BETA/ALPHAR/BETAR/POLY/NGRID/NSEG/MONITO/NOUT/
C     NEVENT/KLOSSE/KLOSSC/KORIF/KSPILL/PFC/
C     STRUC*/DUAL/NELEM*/SOIL/VEG, none; XAREA, sq.m; DCULV/ 
C     RNWL/RMWL/GIM/GIM1/LCULV/EII/EOI/LSPILL/DELTX*/*WIDTH/CANHT, 
C     m; GSC/GSC1, per sq.m; TSLOPE/WSLOPE/FSLOPE/CSLOPE/SCULV/
C     GI/SCCOR/MANN*/PRESD/KTC/DW/NCULV, fraction; KRET, mm;      
C     BACAN/CAN/OPEN/BAOPEN/COVER/ASW, %; BFLOW, m3/s  
	READ (2,*)
      READ (1,*)
      READ (1,*) NGRID,NSEG,BETA,DELTT,RNWL,TWIDTH
      READ (1,*) RMWL,TSLOPE,WSLOPE,BFLOW,MONITO,NOUT,GIM,GSC     
      READ (1,*) POLY(1),POLY(2),POLY(3),POLY(4),PRIN,ENDS,KB,KP
      READ (1,*) LCULV,EII,EOI,DCULV,KLOSSE,NEVENT,KSTORE    
      READ (1,*) KLOSSC,KORIF,KSPILL,LSPILL,DELTX1,NCULV,KRET,	
     *WWIDTH 
	READ (2,*) CAF,Lm,LATRAD

      IF(LCULV .NE. 0.)THEN
        SCULV=(EII-EOI)/LCULV
      ENDIF
      XAREA=3.14*(DCULV**2.)/4.
      DO 5 L=1,NGRID
         READ (1,*)
         READ (1,*) BACAN(L),CAN(L),OPEN(L),BAOPEN(L),COVER(L)
         READ (1,*) GI(L),PFC(L),SCCOR(L),ASW(L),DELTX2(L)
         READ (1,*) MANN1(L),MANN2(L),MANN3(L),FSLOPE(L),CWIDTH(L) 
         READ (1,*) CSLOPE(L),CANHT(L),PRESD(L),KTC(L),DW(L)      
         READ (1,*) DSUB(L),DSURF(L),STRUC1(L),STRUC2(L),STRUC3(L) 
         READ (1,*) DUAL(L),NELEMV(L),NELEMC(L),NELEMT(L),NELEMW(L)
         READ (1,*) SOIL(L),VEG(L),GIM1(L),GSC1(L)  	  
	   READ (1,*) AREA1(L),AREA2(L),RK(L),RB(L)
	   READ (7,*) Z(L)
  5   CONTINUE
C     WATERSHED DATA BY COMPUTATION ...............................
C     UNITS: KSTO, days; FGCORR/CANFAC/FRI/RI/MOUT/ETK, none;     
C     KOND/KONDC/KONDO/KONDE, cm/hr; KRILL, s/m; TAWC, Pa; KRILLY, 
C     Kg s m-4; PCAN/CANEF/GRCE, fraction; FSANG, radians; 
C     RETAIN, mm
      DO 6 L=1,NGRID
         IF((PFC(L) .GE. 1.) .AND. (PFC(L) .NE. 52.6316))THEN
           FGCORR(L)=0.1
         ELSE
           FGCORR(L)=2.-0.019*PFC(L)
         ENDIF
         L1=SOIL(L)
         KOND(L)=KONDS(L1)*FGCORR(L)*SCCOR(L)
         CANFAC(L)=1.+COVER(L)/100.  
         KONDC(L)=KOND(L)*CANFAC(L)*((BACAN(L)/CAN(L))*SCCOR(L)+
     *   PORFAC(L1)*(1.-BACAN(L)/CAN(L)))
         KONDO(L)=KOND(L)*((BAOPEN(L)/OPEN(L))*SCCOR(L)+
     *   PORFAC(L1)*(1.-BAOPEN(L)/OPEN(L)))
         KONDE(L)=(CAN(L)*KONDC(L))/100.+(OPEN(L)*KONDO(L))/100.
         KRILL(L)=(196+(0.015*(WISCHM(L1)**0.2)*((WISCHM(L1)**0.8)-
     *   3500))-(32.7/(CECNH4(L1)**0.4))+
     *  ((0.16/(SAR(L1)**0.75))*(1./AL(L1)-(50.*(SAR(L1)**1.75))))+
     *   (7.98*((ORGMA(L1)**1.5)+3.5)/(ORGMA(L1)**0.5)))/1000
         IF(KRILL(L) .LT. 0.)THEN
           KRILL(L)=0.
         ENDIF 
         IF(CLAY(L1) .GT. 30.)THEN
           TAWC(L)=-0.5-284.*ASW(L)/100.*(ASW(L)/100.-0.3)        
         ELSE
           TAWC(L)=-2.85-(8.87/((VFSAND(L1)+0.1)**0.2))-          
     *     (0.16*CACO3(L1))+(3.65*SAR(L1))+
     *     (3.79*(SPSURF(L1)**0.2))+
     *     (28.1*((WDCLAY(L1)/CLAY(L1))**0.8)*1./(SAND(L1)**0.3))
         ENDIF
         IF((CLAY(L1) .GT. 35.) .AND. (STAB(L1) .GT. 17.) .AND.   
     *   (STAB(L1) .LT. 19.))THEN 
           KRILLY(L)=2670000.
         ELSE IF(CLAY(L1) .GT. 35.)THEN
           KRILLY(L)=(2.67-0.115*(ALOG((18.-                     
     *     STAB(L1))**2.)))*1000000.
         ELSE
           KRILLY(L)=(-2.92-2.71*(WDCLAY(L1)/CLAY(L1))-0.51*MG(L1)+
     *     0.10*WDCLAY(L1)+4.19*((CLAY(L1)/(FE(L1)+AL(L1)))**0.16)+
     *     1.24*COND(L1))*1000000.
         ENDIF
         IF(KRILLY(L) .LT. 0.)THEN
           KRILLY(L)=0.
         ENDIF 
         FSANG(L)=ASIN(FSLOPE(L)) 
         FSFAC(L)=1.05-0.85*EXP(-4.*SIN(FSANG(L)))
         PCAN(L)=CAN(L)/100.
         CANEF(L)=1.-PCAN(L)*EXP(-0.34*CANHT(L))
         GRCE(L)=EXP(-2.5*PRESD(L))
         KSTO(L)=(1000.*CLAY(L1)/(CLAY(L1)+EXP(10.047-           
     *   0.148*CLAY(L1))))*KSTORE
         FRI(L)=1.63*(GIM1(L)**0.29)*(GSC1(L)**0.12)
         RETAIN(L)=KRET*FRI(L)/(2.5*FSLOPE(L))
  6   CONTINUE 
      MOUT=INT((PRIN/DELTT)+0.5)                                  
      RI=1.63*(GIM**0.29)*(GSC**0.12)                   
      ETK=INT(30./DELTT)
C     INPUT: ANTECEDENT VALUES OF STATE VARIABLES, AND EVENTS
C     LIST; EVENTS' DATA (RAIN/EPAN) INPUT BY PIECE-MEAL DURING   
C     EXECUTION .................................................
C     NOTE (1): DETENTION POND INITIAL STORAGE (m3) MUST BE INPUT 
C     FOR EACH LOCATION OF POND, EVEN IF ZERO. ONLY ONE LOCATION OF
C     POND IS PERMITTED IN A WATERSHED ELEMENT - AT THE LOWER END 
C     OF A CHANNEL SEGMENT  
C     NOTE (2) - UNITS: CIS/SUBS/SURFS/EPAN/POND, mm; CFILT, cm;  
C     FILT, cm/hr; FLOW*/RFLOW/FLOW/Q*UP, m3/s; RSTORE, m3;       
C     G*UP/SEDLD*, kg/s; GVUP/SEDLDV, kg/s/m; STIME/STIMES, mins; 
C     RAIN, mm/hr 
      DO 7 L=1,NGRID
         CIS(L)=0.
         CFILT(L)=0.
         FILT(L)=0.
         FLOWV(L)=0.
         QVLAT(L)=0.
         FLOWC(L)=BFLOW
         QCLAT(L)=0.
         FLOWT(L)=0.
         QTLAT(L)=0.
         FLOWW(L)=0.
         QWLAT(L)=0.
         RFLOW(L)=0.         
         QSUB(L)=0.
         IF(STRUC3(L) .EQ. 1)THEN
           READ (1,*)
           READ (1,*) RSTORE(L)
         ELSE
           RSTORE(L)=0.
         ENDIF
         STORE=POLY(1)+(POLY(2)*(RNWL-0.05))+
     *   (POLY(3)*(RNWL-0.05)**2)+(POLY(4)*(RNWL-0.05)**3)
         IF(RSTORE(L) .GT. STORE)THEN
           RSTORE(L)=STORE
         ENDIF 
         RCOMP(L)=2.*RSTORE(L)/(DELTT*60.)
         QVUP(L)=0.   
         QCUP(L)=0.
         QTUP(L)=0.
         QWUP(L)=0.
         QGUP(L)=0.
         SEDLDV(L)=0.
         SEDLDC(L)=0.
         SEDLDT(L)=0.
         SEDLDW(L)=0.
         GVUP(L)=0.
         GCUP(L)=0.
         GTUP(L)=0.
         GWUP(L)=0.
         POND(L)=0.
         TEST(L)=0
         L1=SOIL(L)
         ASW(L)=ASW(L)-WP(L1) 
         IF(ASW(L)/100. .GT. EFFPOR(L1))THEN
           ASW(L)=EFFPOR(L1)*100.
         ENDIF  
         SUBS(L)=DSUB(L)*10.*ASW(L)/100.
         SURFS(L)=DSURF(L)*10.*ASW(L)/100.
  7   CONTINUE 
      STIME=0.
      NEXTE=1
      NEVENT=NEVENT+1
      DO 8 L=1,NEVENT
         IF(L .EQ. NEVENT)THEN
           STIMES(L)=ENDS
         ELSE
           READ (1,*)
           READ (1,*) STIMES(L) 
         ENDIF
  8   CONTINUE 	
C     OUTPUT HYDROGRAPH/SEDIMENTGRAPH INITIAL VALUES ........
        CONC=0. 
      WRITE (6,FMT='(5X,I8,2X,F13.2,5X,F13.4,5X,I8)') INT(STIME), 
     *FLOWC(MONITO)*1000,SEDLDC(MONITO),INT(CONC*1000)
      WRITE (*,FMT='(5X,I8,2X,F13.2,5X,F13.4,5X,I8)') INT(STIME), 
     *FLOWC(MONITO)*1000,SEDLDC(MONITO),INT(CONC*1000)
      CLOSE (6)
      RETURN
      END
C     ..........................................................
C     SUBROUTINE CALLED IN MAIN PROGRAM (1ST ORDER ROUTINE): 
C     TIMING ROUTINE ...........................................
      SUBROUTINE TIMER(NEXTE,STIMES,DELTT,NSIM,ETT,EVENT1)
	  
	  !DEC$ ATTRIBUTES DLLEXPORT, STDCALL :: TIMER
	  !DEC$ ATTRIBUTES ALIAS:'TIMER' :: TIMER
      !DEC$ ATTRIBUTES REFERENCE :: NEXTE,STIMES,DELTT,NSIM,ETT,EVENT1
	  
	  
      REAL STIMES(609)
      CHARACTER EVENT1*20
      IEVENT=NEXTE
      NEXTE=NEXTE+1
      L=INT(IEVENT/2)
      IF((2*L) .EQ. IEVENT)THEN
        EVENT1='EVAPOTRANSPIRATION'
        ETT=STIMES(NEXTE)-STIMES(IEVENT)
      ELSE
        EVENT1='RAINFALL'
        NSIM=INT(((STIMES(NEXTE)-STIMES(IEVENT))/DELTT)+0.5)
      ENDIF        
      RETURN        
      END
C     ..........................................................
C     SUBROUTINE CALLED IN MAIN PROGRAM (1ST ORDER ROUTINE): 
C     RAINFALL EVENT ROUTINE ...................................
      SUBROUTINE RAINFL(NSIM,STIME,DELTT,NGRID,VEG,SOIL,CIS,CISCAP,
     *CFILT,FILT,KONDE,EFFPOR,ASW,WFCP,DUAL,QVUP,DELTX1,
     *QGUP,DSUB,SAT,WP,SUBS,RAIN,RETAIN,FLOWV,QVLAT,MANN1,
     *FSLOPE,BETA,FLOWC,QCLAT,FLOWW,DELTX2,MANN2,CSLOPE,QCUP,
     *CWIDTH,DWBA,DW,KTC,KRILL,TAWC,GVUP,KRILLY,FSFAC,CANEF,
     *GRCE,SEDLDV,GCUP,PSIZE,SEDLDC,DSURF,KSTO,QSUB,STRUC1,
     *STRUC2,STRUC3,RFLOW,RCOMP,RSTORE,POLY,RNWL,RMWL,KSPILL,
     *LSPILL,DCULV,EII,EOI,XAREA,KLOSSE,KLOSSC,LCULV,SCULV,   
     *KORIF,NCULV,FLOWT,QTLAT,MANN3,TSLOPE,QTUP,TWIDTH,GTUP,RI,
     *QWLAT,WSLOPE,QWUP,WWIDTH,GWUP,MOUT,MONITO,NELEMV,NELEMC,
     *NELEMT,NELEMW,SEDLDT,SEDLDW,SURFS,NAME,POND,TEST,
     *BFLOW,NOUT,FC,NSEG,FLOW,AREA1,AREA2,RK,RB,CAF)
	 
	 
	 !DEC$ ATTRIBUTES DLLEXPORT, STDCALL :: RAINFL
	  !DEC$ ATTRIBUTES ALIAS:'RAINFL' :: RAINFL
      !DEC$ ATTRIBUTES REFERENCE :: NSIM,STIME,DELTT,NGRID,VEG,SOIL,CIS,CISCAP,
     *CFILT,FILT,KONDE,EFFPOR,ASW,WFCP,DUAL,QVUP,DELTX1,
     *QGUP,DSUB,SAT,WP,SUBS,RAIN,RETAIN,FLOWV,QVLAT,MANN1,
     *FSLOPE,BETA,FLOWC,QCLAT,FLOWW,DELTX2,MANN2,CSLOPE,QCUP,
     *CWIDTH,DWBA,DW,KTC,KRILL,TAWC,GVUP,KRILLY,FSFAC,CANEF,
     *GRCE,SEDLDV,GCUP,PSIZE,SEDLDC,DSURF,KSTO,QSUB,STRUC1,
     *STRUC2,STRUC3,RFLOW,RCOMP,RSTORE,POLY,RNWL,RMWL,KSPILL,
     *LSPILL,DCULV,EII,EOI,XAREA,KLOSSE,KLOSSC,LCULV,SCULV,   
     *KORIF,NCULV,FLOWT,QTLAT,MANN3,TSLOPE,QTUP,TWIDTH,GTUP,RI,
     *QWLAT,WSLOPE,QWUP,WWIDTH,GWUP,MOUT,MONITO,NELEMV,NELEMC,
     *NELEMT,NELEMW,SEDLDT,SEDLDW,SURFS,NAME,POND,TEST,
     *BFLOW,NOUT,FC,NSEG,FLOW,AREA1,AREA2,RK,RB,CAF
	 
	 
      IMPLICIT REAL(K)
C     SOIL TYPE PARAMETERS ....................................
      REAL DWBA(3),PSIZE(3),WP(3),SAT(3),EFFPOR(3),WFCP(3),FC(3)  
C     VEGETATION TYPE PARAMETER ...............................
      REAL CISCAP(8)
C     SYSTEM AND WATERSHED PARAMETERS ........................... 
      REAL DELTX2(73),MANN1(73),MANN2(73),GRCE(73),KTC(73),  
     *DW(73),KSTO(73),DSUB(73),DSURF(73),MANN3(73),ASW(73), 
     *KRILL(73),TAWC(73),KRILLY(73),FSFAC(73),CANEF(73),     
     *CWIDTH(73),POLY(4),RETAIN(73),KONDE(73),LCULV,NCULV,NOUT,   
     *LSPILL,FSLOPE(73),CSLOPE(73),AREA1(73),AREA2(73),RK(73),RB(73)
      INTEGER SOIL(73),VEG(73),DUAL(73),STRUC1(73),STRUC2(73),
     *STRUC3(73),NELEMV(73),NELEMC(73),NELEMT(73),NELEMW(73), 
     *TEST(73),DUMMY
C     SYSTEM AND STATE VARIABLES.................................
      REAL CIS(73),FILT(73),CFILT(73), 
     *FLOWV(73),FLOWC(73),FLOWT(73),FLOWW(73),QVLAT(73), 
     *QCLAT(73),QTLAT(73),QWLAT(73),QGUP(73),QVUP(73),
     *QCUP(73),QTUP(73),QWUP(73),SEDLDV(73),SEDLDC(73), 
     *SEDLDT(73),SEDLDW(73),SUBS(73),SURFS(73),QSUB(73), 
     *POND(73),RFLOW(73),RSTORE(73),GVUP(73),GCUP(73),   
     *GTUP(73),GWUP(73),RCOMP(73),CAF,FRFDG,EFWC
      CHARACTER NAME*15
      OPEN (6,FILE=NAME,POSITION='APPEND')
      READ (3,*)	
      DO 14 J=1,NSIM
        STIME=STIME+DELTT
         READ (3,FMT='(7X,F17.2)',END=100,ERR=100) RAIN
         DO 12 I=1,NGRID    
            L=VEG(I)
            L1=SOIL(I)  
            IF(CIS(I) .NE. CISCAP(L))THEN
              CFILT(I)=(RAIN*DELTT/60.-CISCAP(L)+CIS(I))/10.
              IF(CFILT(I) .GT. 0.)THEN
                CIS(I)=CISCAP(L)
                RAIN=CFILT(I)*10./(DELTT/60.)
                FILT(I)=KONDE(I)*(1.+((EFFPOR(L1)-                
     *          ASW(I)/100.)*WFCP(L1)/CFILT(I)))                  
                DUMMY=1
              ELSE 
C               CIS(I)=CIS(I)+(RAIN*DELTT/60.) -OLD
C	---------------------NEW EQUATION ---------------------			  
C			  CIS(I)=CIS(I)+(DELTT*((AREA1(I)*AREA2(I)*(RAIN
C     *		  -(PET*CIS(I))/CISCAP(L)))-(RK(I)*EXP(RB(I)*(CIS(I)-
C     *		  CISCAP(L)))))) 
C	CCC is contant characteristic of the canopy
C	EFWC is the evaporation from wet canopy, 
C	FRFDG is the fraction of rain falling directly to the ground
			  FRFDG=0.26
			  EFWC=0.15
			  CIS(I)=CIS(I)+DELTT*((1-FRFDG)*RAIN-(CIS(I)/
     *		  CISCAP(L)*EFWC)-RK(I)*EXP(CIS(I)/RB(I)))			  
C	---------------------END OF NEW EQUATION ---------------------	
                IF(CIS(I) .GT. CISCAP(L))THEN
                  CFILT(I)=CIS(I)-CISCAP(L)
                  CIS(I)=CISCAP(L)
                  FILT(I)=KONDE(I)*(1.+((EFFPOR(L1)-              
     *            ASW(I)/100.)*WFCP(L1)/CFILT(I)))                
                  DUMMY=1
                ELSE
                  CFILT(I)=0.
                  FILT(I)=0.
                ENDIF
              ENDIF
            ELSE
              DUMMY=0
            ENDIF
            IF((CIS(I) .EQ. CISCAP(L)) .AND. (TEST(I) .EQ. 1))THEN
              IF(FILT(I) .LE. RAIN/10.)THEN
                FILTX=FILT(I)
                CALL PONDED(CFILT(I),FILT(I),DELTT,KONDE(I),     
     *          WFCP(L1),EFFPOR(L1),ASW(I))
                FILTX=(FILTX+FILT(I))/2.
                EXRAIN=RAIN-(10.*FILTX)
                IF(EXRAIN .LT. 0.)THEN
                  EXRAIN=0.
                ENDIF 
                POND(I)=EXRAIN*DELTT/60.
                SURFS(I)=DSURF(I)*10.*(SAT(L1)-WP(L1))/100.
                QSUB(I)=(FILTX/100.)*DELTX1*(DSUB(I)/100.)/3600.
              ELSE
                FILTX=FILT(I)
                CFILT(I)=CFILT(I)+(RAIN*DELTT/60.)/10.
                FILT(I)=KONDE(I)*(1.+((EFFPOR(L1)-                
     *          ASW(I)/100.)*WFCP(L1)/CFILT(I)))
                FILTX=(FILTX+FILT(I))/2.
                EXRAIN=RAIN-(10.*FILTX)
                IF(EXRAIN .LT. 0.)THEN
                  EXRAIN=0.
                ENDIF 
                POND(I)=EXRAIN*DELTT/60.
                SURFS(I)=DSURF(I)*10.*(SAT(L1)-WP(L1))/100.
                QSUB(I)=(FILTX/100.)*DELTX1*(DSUB(I)/100.)/3600.
              ENDIF
            ELSE IF((CIS(I) .EQ. CISCAP(L)) .AND. 
     *        (TEST(I) .EQ. 0))THEN
              IF(DUMMY .EQ. 1)THEN
                CFILTX=CFILT(I)
              ELSE  
                CFILTX=CFILT(I)+(RAIN*DELTT/60.)/10.
              ENDIF 
              IF(CFILTX .GT. 0.)THEN  
                FILTX=KONDE(I)*(1.+((EFFPOR(L1)-                
     *          ASW(I)/100.)*WFCP(L1)/CFILTX))
                IF((FILTX .LE. RAIN/10.) .AND. (POND(I) .EQ. 0.)
     *            .AND. (FILTX .GT. 0.) .AND. (RAIN .GT. 0.))THEN 
                  PFILT=(EFFPOR(L1)-                              
     *            ASW(I)/100.)*WFCP(L1)/(RAIN/10./KONDE(I)-1.)
                  CFILT(I)=PFILT
                  FILT(I)=KONDE(I)*(1.+((EFFPOR(L1)-              
     *            ASW(I)/100.)*WFCP(L1)/CFILT(I)))
                  EXRAIN=RAIN-(10.*FILT(I))
                  IF(EXRAIN .LT. 0.)THEN
                    EXRAIN=0.
                  ENDIF 
                  POND(I)=EXRAIN*DELTT/60.
                  SURFS(I)=CFILT(I)*10.
                ELSE IF((FILTX .LE. RAIN/10.) .AND. 
     *            (POND(I) .GT. 0.) .AND. (FILTX .GT. 0.) .AND.
     *            (RAIN .GT. 0.))THEN 
                  FILTX=FILT(I)
                  CALL PONDED(CFILT(I),FILT(I),DELTT,KONDE(I),    
     *            WFCP(L1),EFFPOR(L1),ASW(I))
                  FILTX=(FILTX+FILT(I))/2.
                  EXRAIN=RAIN-(10.*FILTX)
                  IF(EXRAIN .LT. 0.)THEN
                    EXRAIN=0.
                  ENDIF 
                  POND(I)=POND(I)+EXRAIN*DELTT/60.
                  IF(POND(I) .GE. RETAIN(I))THEN
                    TEST(I)=1 
                    SURFS(I)=DSURF(I)*10.*(SAT(L1)-WP(L1))/100.
                  ENDIF
                ENDIF
              ELSE
                FILT(I)=0.
                EXRAIN=0.
              ENDIF
            ENDIF
            IF(DUAL(I) .EQ. 1)THEN     
              IF((NELEMC(I) .LT. NGRID) .AND. 
     *          (NELEMW(I) .NE. 0))THEN
                L3=NELEMW(I)
                FLOWW(I)=FLOWW(L3)
              ENDIF
              CALL SATROF(DELTT,QGUP(I),DELTX1,DSUB(I),SAT(L1),
     *        WP(L1),SUBS(I),RAIN,POND(I),QVUP(I),FILTX,FLOWV(I),
     *        QVLAT(I),MANN1(I),FSLOPE(I),BETA,FLOWC(I),QCLAT(I), 
     *        FLOWW(I),DELTX2(I),MANN2(I),CSLOPE(I),QCUP(I),
     *        CWIDTH(I),DWBA(L1),DW(I),FC(L1),KTC(I),KRILL(I),
     *        TAWC(I),GVUP(I),KRILLY(I),FSFAC(I),CANEF(I),GRCE(I),
     *        SEDLDV(I),GCUP(I),PSIZE(L1),SEDLDC(I),DEPTH,SURFS(I),
     *        DSURF(I),KSTO(I),TEST(I),QSUB(I),BFLOW,NOUT,CAF)
            ELSE 
              CALL RUNOFF(SAT(L1),WP(L1),DSUB(I),DELTT,FILTX,  
     *        QGUP(I),DELTX1,SUBS(I),KSTO(I),QSUB(I),POND(I),
     *        FLOWV(I),QVLAT(I),MANN1(I),FSLOPE(I),QVUP(I),     
     *        BETA,DWBA(L1),DW(I),KTC(I),KRILL(I),TAWC(I),
     *        RAIN,GVUP(I),KRILLY(I),FSFAC(I),CANEF(I),GRCE(I),
     *        SEDLDV(I),TEST(I),FC(L1),SURFS(I),DSURF(I),NOUT,CAF)	    
            ENDIF
            IF((STRUC1(I) .EQ. 1) .AND. (DUAL(I) .EQ. 1) .AND.
     *        (STRUC3(I) .EQ. 1))THEN                           
              CALL HSTRUC(FLOWC(I),RFLOW(I),RCOMP(I),RSTORE(I),
     *        POLY,RNWL,RMWL,KSPILL,LSPILL,SEDLDC(I),DELTT,FLOW)
            ELSE IF((STRUC1(I) .EQ. 1) .AND.
     *        (DUAL(I) .EQ. 1))THEN 
              CALL CULVT(DEPTH,DCULV,EII,EOI,XAREA,KLOSSE,KLOSSC,
     *        LCULV,SCULV,KORIF,NCULV,FLOWC(I))                 
            ELSE IF(STRUC1(I) .EQ. 1)THEN                       
              CALL STRUCT(FLOWT(I),QTLAT(I),FLOWV(I),DELTX1,
     *        MANN3(I),TSLOPE,DELTT,QTUP(I),BETA,TWIDTH,GTUP(I),
     *        SEDLDV(I),PSIZE(L1),TAWC(I),KRILL(I),SEDLDT(I),
     *        STRUC2(I),RI,FLOWW(I),QWLAT(I),DELTX2(I),WSLOPE,
     *        QWUP(I),WWIDTH,GWUP(I),SEDLDW(I))          
            ENDIF     
  12     CONTINUE           
         DO 13 I=1,NGRID
            L=NELEMV(I)         
            L1=NELEMC(I)
            L2=NELEMT(I)
            L3=NELEMW(I)
            IF(L .EQ. 0)THEN
              QVUP(I)=0.
              QGUP(I)=0.
              GVUP(I)=0.
            ELSE IF(L .GT. NGRID)THEN
              L=L-NGRID
              QVUP(I)=FLOWV(L)+FLOWV(L2) 
              QGUP(I)=QSUB(L)+QSUB(L2)
              GVUP(I)=SEDLDV(L)+SEDLDV(L2)  
            ELSE
              QVUP(I)=FLOWV(L)
              QGUP(I)=QSUB(L)
              GVUP(I)=SEDLDV(L)
            ENDIF    
            IF(L1 .EQ. 0)THEN
              IF(I .EQ. NSEG)THEN
                QCUP(I)=BFLOW
                GCUP(I)=0.
              ELSE
                QCUP(I)=0.
                GCUP(I)=0.
              ENDIF
            ELSE IF(L1 .GT. NGRID)THEN
              L1=L1-NGRID
              QCUP(I)=FLOWC(L1)+FLOWC(L3)-BFLOW                   
              GCUP(I)=SEDLDC(L1)+SEDLDC(L3)
            ELSE
              QCUP(I)=FLOWC(L1)
              GCUP(I)=SEDLDC(L1)
            ENDIF    
            IF(L2 .EQ. 0)THEN
              QTUP(I)=0.
              GTUP(I)=0.
            ELSE
              QTUP(I)=FLOWT(L2)
              GTUP(I)=SEDLDT(L2)
            ENDIF    
            IF(L3 .EQ. 0)THEN
              QWUP(I)=0.
              GWUP(I)=0.
            ELSE
              QWUP(I)=FLOWW(L3)
              GWUP(I)=SEDLDW(L3)
            ENDIF    
  13     CONTINUE           
         L5=INT(J/MOUT)
         IF((MOUT*L5) .EQ. J)THEN
           IF(FLOWC(MONITO) .GT. 0.)THEN 
             CONC=SEDLDC(MONITO)/FLOWC(MONITO)
           ELSE
             CONC=0.
           ENDIF 
           WRITE (6,FMT='(5X,I8,2X,F13.2,5X,F13.4,5X,I8)')        
     *     INT(STIME),FLOWC(MONITO)*1000,SEDLDC(MONITO),
     *     INT(CONC*1000)
           WRITE (*,FMT='(5X,I8,2X,F13.2,5X,F13.4,5X,I8)')        
     *     INT(STIME),FLOWC(MONITO)*1000,SEDLDC(MONITO),
     *     INT(CONC*1000)
         ENDIF              
  14  CONTINUE
      CLOSE (6)
      RETURN
 100  WRITE (*,*) 'RAINFALL DATA INPUT ERROR'
      STOP 100
      END
C     .............................................................
C     PONDED-INFILTRATION ROUTINE (2ND ORDER ROUTINE)..............
      SUBROUTINE PONDED(CFILT1,FILT1,DELTT1,KONDE1,WFCP1,     
     *EFFPO1,ASW1)
	 
	 !DEC$ ATTRIBUTES DLLEXPORT, STDCALL :: PONDED
	  !DEC$ ATTRIBUTES ALIAS:'PONDED' :: PONDED
      !DEC$ ATTRIBUTES REFERENCE :: CFILT1,FILT1,DELTT1,KONDE1,WFCP1,     
     *EFFPO1,ASW1
	 
	 
	 
      REAL KONDE1
      TEMP=CFILT1
      CFILT1=CFILT1+(KONDE1*DELTT1/60.)
  21  TEMP1=CFILT1
      XYZ=(TEMP1+(WFCP1*(EFFPO1-ASW1/100.)))/
     *(TEMP+(WFCP1*(EFFPO1-ASW1/100.)))
      IF(XYZ .GT. 0.)THEN   
        CFILT1=(KONDE1*DELTT1/60.)+TEMP+(WFCP1*(EFFPO1-           
     *  ASW1/100.))*ALOG(XYZ)
        IF((CFILT1 .GT. 0) .AND. 
     *    ((ABS(CFILT1-TEMP1)/TEMP1) .LE. 0.05))THEN
          FILT1=KONDE1*(1.+((EFFPO1-ASW1/100.)*WFCP1/CFILT1))  
        ELSE IF(CFILT1 .GT. 0.)THEN
          GO TO 21
        ENDIF
      ELSE IF(CFILT1 .GT. 0.)THEN
        FILT1=KONDE1*(1.+((EFFPO1-ASW1/100.)*WFCP1/CFILT1))  
      ENDIF
      RETURN
      END
C     ..........................................................
C     SUBROUTINE CALLED IN MAIN PROGRAM (1ST ORDER ROUTINE): 
C     EVAPOTRANSPIRATION EVENT ROUTINE .......................
      SUBROUTINE EVAPT(STIME,DELTT,ETK,NGRID,VEG,SOIL,FILT,CFILT,
     *KONDE,WFCP,EFFPOR,ASW,DUAL,DSUB,SAT,WP,SUBS,QVUP,
     *FLOWV,QVLAT,MANN1,MANN2,FSLOPE,BETA,FLOWC,QCLAT,FLOWW,DELTX1,
     *DELTX2,CSLOPE,QCUP,CWIDTH,DWBA,DW,KTC,KRILL,TAWC,GVUP,KRILLY,
     *FSFAC,CANEF,GRCE,SEDLDV,GCUP,PSIZE,SEDLDC,DSURF,FC,QGUP,KSTO,
     *QSUB,STRUC1,STRUC2,STRUC3,RFLOW,RCOMP,RSTORE,POLY,RNWL,RMWL,
     *KSPILL,LSPILL,DCULV,EII,EOI,XAREA,KLOSSE,KLOSSC,LCULV,SCULV,
     *KORIF,NCULV,FLOWT,QTLAT,MANN3,TWIDTH,WWIDTH,TSLOPE,WSLOPE,
     *QTUP,GTUP,RI,QWLAT,QWUP,GWUP,GI,PCAN,LAI,KC,KP,KB,CIS,MOUT,
     *MONITO,NELEMV,NELEMC,NELEMT,NELEMW,SEDLDT,SEDLDW,ETT,SURFS,
     *TOTPOR,NOUT,BFLOW,CLAY,VFSAND,CACO3,SAR,SPSURF,WDCLAY,SAND,
     *NAME,POND,TEST,RETAIN,NSEG,FLOW,W,Z,Lm,LATRAD,CAF)
	 
	 
	 
	 !DEC$ ATTRIBUTES DLLEXPORT, STDCALL :: EVAPT
	  !DEC$ ATTRIBUTES ALIAS:'EVAPT' :: EVAPT
      !DEC$ ATTRIBUTES REFERENCE :: STIME,DELTT,ETK,NGRID,VEG,SOIL,FILT,CFILT,
     *KONDE,WFCP,EFFPOR,ASW,DUAL,DSUB,SAT,WP,SUBS,QVUP,
     *FLOWV,QVLAT,MANN1,MANN2,FSLOPE,BETA,FLOWC,QCLAT,FLOWW,DELTX1,
     *DELTX2,CSLOPE,QCUP,CWIDTH,DWBA,DW,KTC,KRILL,TAWC,GVUP,KRILLY,
     *FSFAC,CANEF,GRCE,SEDLDV,GCUP,PSIZE,SEDLDC,DSURF,FC,QGUP,KSTO,
     *QSUB,STRUC1,STRUC2,STRUC3,RFLOW,RCOMP,RSTORE,POLY,RNWL,RMWL,
     *KSPILL,LSPILL,DCULV,EII,EOI,XAREA,KLOSSE,KLOSSC,LCULV,SCULV,
     *KORIF,NCULV,FLOWT,QTLAT,MANN3,TWIDTH,WWIDTH,TSLOPE,WSLOPE,
     *QTUP,GTUP,RI,QWLAT,QWUP,GWUP,GI,PCAN,LAI,KC,KP,KB,CIS,MOUT,
     *MONITO,NELEMV,NELEMC,NELEMT,NELEMW,SEDLDT,SEDLDW,ETT,SURFS,
     *TOTPOR,NOUT,BFLOW,CLAY,VFSAND,CACO3,SAR,SPSURF,WDCLAY,SAND,
     *NAME,POND,TEST,RETAIN,NSEG,FLOW,W,Z,Lm,LATRAD,CAF
	 
	 
      IMPLICIT REAL(K)
C     SOIL TYPE PARAMETERS .......................................
      REAL CLAY(3),SAND(3),DWBA(3),PSIZE(3),FC(3),WP(3),SAT(3),
     *WDCLAY(3),SAR(3),VFSAND(3),SPSURF(3),CACO3(3),              
     *EFFPOR(3),TOTPOR(3),WFCP(3)
C     VEGETATION TYPE PARAMETERS .................................
      REAL KC(8),LAI(8)
C     SYSTEM AND WATERSHED PARAMETERS ............................ 
      REAL GI(73),DELTX2(73),MANN1(73),MANN2(73),PCAN(73),   
     *GRCE(73),KTC(73),DW(73),KSTO(73),DSUB(73),DSURF(73),  
     *MANN3(73),ASW(73),KRILL(73),TAWC(73),KRILLY(73),       
     *FSFAC(73),CANEF(73),CWIDTH(73),POLY(4),RETAIN(73),
     *KONDE(73),LCULV,NCULV,NOUT,LSPILL,FSLOPE(73),CSLOPE(73)
     
      INTEGER SOIL(73),VEG(73),DUAL(73),STRUC1(73),STRUC2(73),
     *STRUC3(73),NELEMV(73),NELEMC(73),NELEMT(73),NELEMW(73),
     *TEST(73),ETK,COUNT1,COUNT2,ETFLAG,MARK,MARKER   
C     SYSTEM AND STATE VARIABLES ..................................
      REAL CIS(73),AET(73),FILT(73),CFILT(73), 
     *FLOWV(73),FLOWC(73),FLOWT(73),FLOWW(73),QVLAT(73), 
     *QCLAT(73),QTLAT(73),QWLAT(73),QGUP(73),QVUP(73),
     *QCUP(73),QTUP(73),QWUP(73),SEDLDV(73),SEDLDC(73), 
     *SEDLDT(73),SEDLDW(73),SUBS(73),SURFS(73),QSUB(73), 
     *POND(73),DELTF(73),RFLOW(73),RSTORE(73),GVUP(73),   
     *GCUP(73),GTUP(73),GWUP(73),RCOMP(73),AVPOR(73),RT(73)
C	New variables		N=daylight hours=DAYLHR, n=actual duration
C						of sunshine = ADS,d=Solar Declination SOLARD
C						dr=inverse relative distance earth-sun=IRDES
C						LAT = latitude,Ws=Sunset hour angle=SHA
C						Gsc= Solar constant
	REAL EA,ES,G,RN,RA,DELTA,Pa,Cp,GAMMA,LAMBDA,Zm,Zo,Sn,Ln,
     *RS,Patm,T,ROW,Rsp,KT,U2,Z(73),Rso,ALFA,Tmax,Tmin,Ts,So,RHmean,
     *PW,Eomax,Eomin,DAYLHR,ADS,WS,Rns,RnL,Gsc,SOLARD,LATRAD,IRDES,
     *SHA,W1,W2,W,Lz,Lm,SC,T1,B,CAF
	INTEGER J,MONTH,DAY,YEAR
      CHARACTER NAME*15
      DATA (AET(I),I=1,73) /73*0./
      DATA (DELTF(I),I=1,73) /73*0./
      DATA COUNT2,ETFLAG /0,0/
      TIME=STIME
      READ (5,*)
      OPEN (6,FILE=NAME,POSITION='APPEND')
      MARKER=INT(ETT/30)
      MARK=0
      DO 31 COUNT1=1,900000000
         STIME=STIME+DELTT
         COUNT2=COUNT1/ETK  
         IF(COUNT2 .EQ. (MARK+1))THEN
           ETFLAG=1
C           READ (5,FMT='(7X,F24.2)') EPAN
		 READ (5,*)MONTH,DAY,YEAR,RS,Tmin,Tmax,RHmean,Patm,U2		  		 
           MARKER=MARKER-1
           MARK=MARK+1
         ELSE
           ETFLAG=0
         ENDIF
         DO 32 I=1,NGRID		  
            L=VEG(I)
            L1=SOIL(I)
            RAIN=0.
            POND(I)=0.
            IF(COUNT1 .EQ. 1)THEN
              RT(I)=0.5*RETAIN(I)
            ENDIF
            FILTX=FILT(I)
            CFILTX=CFILT(I)
            DELTF(I)=DELTF(I)-CFILT(I)
            CALL PONDED(CFILT(I),FILT(I),DELTT,KONDE(I),        
     *      WFCP(L1),EFFPOR(L1),ASW(I))
            DELTF(I)=DELTF(I)+CFILT(I)
            IF(RT(I) .GE. (10.*DELTF(I)*DELTT/60.))THEN
              RT(I)=RT(I)-(10.*DELTF(I)*DELTT/60.)      
            ELSE IF(RT(I) .GT. 0.)THEN
              CFILT(I)=CFILTX+RT(I)/10.                     
              FILT(I)=KONDE(I)*(1.+((EFFPOR(L1)-                
     *        ASW(I)/100.)*WFCP(L1)/CFILT(I)))
              RT(I)=0.
            ELSE
              CFILT(I)=CFILTX
              FILT(I)=0.
            ENDIF
            FILTX=(FILTX+FILT(I))/2.
            QSUB(I)=(FILTX/100.)*DELTX1*(DSUB(I)/100.)/3600.
            IF(DUAL(I) .EQ. 1)THEN     
              IF((NELEMC(I) .LT. NGRID) .AND. 
     *          (NELEMW(I) .NE. 0))THEN
                L3=NELEMW(I)
                FLOWW(I)=FLOWW(L3)
              ENDIF
              CALL SATROF(DELTT,QGUP(I),DELTX1,DSUB(I),SAT(L1),
     *        WP(L1),SUBS(I),RAIN,POND(I),QVUP(I),FILTX,FLOWV(I),
     *        QVLAT(I),MANN1(I),FSLOPE(I),BETA,FLOWC(I),QCLAT(I), 
     *        FLOWW(I),DELTX2(I),MANN2(I),CSLOPE(I),QCUP(I),
     *        CWIDTH(I),DWBA(L1),DW(I),FC(L1),KTC(I),KRILL(I),
     *        TAWC(I),GVUP(I),KRILLY(I),FSFAC(I),CANEF(I),GRCE(I), 
     *        SEDLDV(I),GCUP(I),PSIZE(L1),SEDLDC(I),DEPTH,
     *        SURFS(I),DSURF(I),KSTO(I),TEST(I),QSUB(I),BFLOW,NOUT,
     *		CAF)
            ELSE 
              CALL RUNOFF(SAT(L1),WP(L1),DSUB(I),
     *        DELTT,FILTX,QGUP(I),DELTX1,SUBS(I),KSTO(I),QSUB(I),
     *        POND(I),FLOWV(I),QVLAT(I),MANN1(I),FSLOPE(I),
     *        QVUP(I),BETA,DWBA(L1),DW(I),KTC(I),KRILL(I),
     *        TAWC(I),RAIN,GVUP(I),KRILLY(I),FSFAC(I),CANEF(I),
     *        GRCE(I),SEDLDV(I),TEST(I),FC(L1),SURFS(I),DSURF(I),
     *        NOUT,CAF) 
            ENDIF
            IF((STRUC1(I) .EQ. 1) .AND. (DUAL(I) .EQ. 1) .AND.
     *        (STRUC3(I) .EQ. 1))THEN                             
              CALL HSTRUC(FLOWC(I),RFLOW(I),RCOMP(I),RSTORE(I),
     *        POLY,RNWL,RMWL,KSPILL,LSPILL,SEDLDC(I),DELTT,FLOW)
            ELSE IF((STRUC1(I) .EQ. 1) .AND. (DUAL(I) .EQ. 1))THEN 
              CALL CULVT(DEPTH,DCULV,EII,EOI,XAREA,KLOSSE,KLOSSC, 
     *        LCULV,SCULV,KORIF,NCULV,FLOWC(I))                   
            ELSE IF(STRUC1(I) .EQ. 1)THEN
              CALL STRUCT(FLOWT(I),QTLAT(I),FLOWV(I),DELTX1,
     *        MANN3(I),TSLOPE,DELTT,QTUP(I),BETA,TWIDTH,GTUP(I),
     *        SEDLDV(I),PSIZE(L1),TAWC(I),KRILL(I),SEDLDT(I),
     *        STRUC2(I),RI,FLOWW(I),QWLAT(I),DELTX2(I),WSLOPE,
     *        QWUP(I),WWIDTH,GWUP(I),SEDLDW(I))
            ENDIF 
            IF(ETFLAG .EQ. 1)THEN
C              PET=(GI(I)*PCAN(I)*LAI(L)*KC(L)*KP*EPAN)+((1.-      
C    *         PCAN(I))*KB*KP*EPAN)     - OLD  
C	----------NEW EQUATIONS AND VARIABLES    -----------------			 
			 Gsc=0.0820
			 Zm=2 
			 Lz=75.			 
			 T=(Tmin+Tmax)/2.			 
			 T1=0.17
		     J=INT(275.*MONTH/9.-30.+DAY)-2.
			 IF(MONTH.LT. 3)THEN
				J=J+2
			 ENDIF
			 IF((4*INT(YEAR/4.)) .EQ. YEAR .AND. MONTH .GT. 2) THEN		
				J=J+1
	         ENDIF						 
			 IRDES=1+0.033*COS(3.141593/365*J)
			 SOLARD=0.409*SIN(2*3.141593/365*J-1.39)
			 B=2*3.141593*(J-81)/364.
			 SC=0.1645*SIN(2*B)-0.1255*COS(B)-0.025*SIN(B)
			 LATRAD=3.141593/180.*LATRAD
			 SHA=ACOS(-TAN(LATRAD)*TAN(SOLARD))		     
			 GAMMA=0.665*10**(-3)/Patm
			 DELTA=4098*(0.6108*EXP(17.27*T/(T+237.3)))/(T+237.3)**2.			 
			 EA=0.6108*EXP(17.27*T/(T+237.3))*RHmean/100.			 
			 W=3.141593/12*((2.5+0.06667*(Lz-Lm)+SC)-12)
			 W1=W-3.141593*T1/24.
			 W2=W+3.141593*T1/24.

			 RA=12.*60/3.141593*Gsc*IRDES*((W2-W1)*SIN(LATRAD)*			 
     *		 SIN(SOLARD)+COS(LATRAD)*COS(SOLARD)*(SIN(W2)-SIN(W1)))
		     Rso=((0.75+2*Z(I))/100000.)*RA			 
			 Rns=(1-0.23)*RS
	         RnL=(((4.903*10**(-9))*(Tmax+273.16)**4)+
     *		 ((4.903*10**(-9))*(Tmin+273.16)**4))/2.*
     *	     (0.34-0.14*EA**0.5)*(1.35*RS/Rso-0.35)
			 RN=Rns-RnL
			 G=0.1*RN			 			 			 
			 PET=(0.408*DELTA*(RN-G)+(GAMMA*37/(T+273))*U2*((0.6108*
     *		 EXP(17.27*T/(T+237.3)))-EA))/(DELTA+GAMMA*(1+0.34*U2))			 			 
C	---------------- END-------------------------------------               
              IF(PET .GE. CIS(I))THEN
                AET(I)=PET-CIS(I)
                CIS(I)=0.
              ELSE
                CIS(I)=CIS(I)-PET
                AET(I)=PET
              ENDIF
              IF(SURFS(I) .GE. PET)THEN
                AET(I)=PET
                SURFS(I)=SURFS(I)-PET
              ELSE                                                
                H=(SURFS(I)+SUBS(I))/    
     *          ((DSURF(I)+DSUB(I))*10.*(SAT(L1)-WP(L1))/100.)
                IF((H .GT. 0.) .AND. (H .LE. 1.))THEN
                  AVPOR(I)=TOTPOR(L1)*(1.-H)
                ELSE IF(H .GT. 1.)THEN
                  AVPOR(I)=0.
                ELSE
                  AVPOR(I)=TOTPOR(L1)
                ENDIF
C               AET(I)=PET*((TOTPOR(L1)-AVPOR(I))/TOTPOR(L1))**0.1
C	------------------------NEW EQUATION----------------------
			  AET(I)=KC(L)*PET
C	------------------------END OF NEW EQUATION--------------
                SURFS(I)=(SURFS(I)+SUBS(I)-AET(I))*DSURF(I)/
     *          (DSURF(I)+DSUB(I))
                SUBS(I)=(SURFS(I)+SUBS(I)-AET(I))*DSUB(I)/
     *          (DSURF(I)+DSUB(I))
                IF(SURFS(I) .LT. 0.)THEN
                  SURFS(I)=0.
                ENDIF
                IF(SUBS(I) .LT. 0.)THEN
                  SUBS(I)=0.
                ENDIF  
              ENDIF
            ENDIF  
  32    CONTINUE
        DO 33 I=1,NGRID         
           L=NELEMV(I)         
           L1=NELEMC(I)
           L2=NELEMT(I)
           L3=NELEMW(I)
           IF(L .EQ. 0)THEN
             QVUP(I)=0.
             QGUP(I)=0.
             GVUP(I)=0.
           ELSE IF(L .GT. NGRID)THEN
             L=L-NGRID
             QVUP(I)=FLOWV(L)+FLOWV(L2) 
             QGUP(I)=QSUB(L)+QSUB(L2)
             GVUP(I)=SEDLDV(L)+SEDLDV(L2)  
           ELSE
             QVUP(I)=FLOWV(L)
             QGUP(I)=QSUB(L)
             GVUP(I)=SEDLDV(L)
           ENDIF    
           IF(L1 .EQ. 0)THEN
             IF(I .EQ. NSEG)THEN
               QCUP(I)=BFLOW
               GCUP(I)=0.
             ELSE
               QCUP(I)=0.
               GCUP(I)=0.
             ENDIF
           ELSE IF(L1 .GT. NGRID)THEN
             L1=L1-NGRID
             QCUP(I)=FLOWC(L1)+FLOWC(L3)-BFLOW                    
             GCUP(I)=SEDLDC(L1)+SEDLDC(L3)
           ELSE
             QCUP(I)=FLOWC(L1)
             GCUP(I)=SEDLDC(L1)
           ENDIF    
           IF(L2 .EQ. 0)THEN
             QTUP(I)=0.
             GTUP(I)=0.
           ELSE
             QTUP(I)=FLOWT(L2)
             GTUP(I)=SEDLDT(L2)
           ENDIF    
           IF(L3 .EQ. 0)THEN
             QWUP(I)=0.
             GWUP(I)=0.
           ELSE
             QWUP(I)=FLOWW(L3)
             GWUP(I)=SEDLDW(L3)
           ENDIF    
  33    CONTINUE
        L6=INT(COUNT1/MOUT)
        IF((MOUT*L6) .EQ. COUNT1)THEN
          IF(FLOWC(MONITO) .GT. 0.)THEN 
            CONC=SEDLDC(MONITO)/FLOWC(MONITO)
          ELSE
            CONC=0.
          ENDIF 
          WRITE (6,FMT='(5X,I8,2X,F13.2,5X,F13.4,5X,I8)')         
     *    INT(STIME),FLOWC(MONITO)*1000,SEDLDC(MONITO),
     *    INT(CONC*1000)
          WRITE (*,FMT='(5X,I8,2X,F13.2,5X,F13.4,5X,I8)')         
     *    INT(STIME),FLOWC(MONITO)*1000,SEDLDC(MONITO),
     *    INT(CONC*1000)
        ENDIF              
        IF(((FLOWC(MONITO)-BFLOW) .LT. 0.00075) .OR.
     *    (MARKER .EQ. 0))THEN
          IF(MARKER .EQ. 0)THEN
            GO TO 36
          ENDIF
          STIME=TIME+(MARK*30)
          ETT=ETT-(MARK*30)
          CALL BASEFL(ETT,STIME,NGRID,VEG,SOIL,DUAL,     
     *    STRUC1,STRUC3,SEDLDT,FLOWW,SEDLDW,GI,PCAN,LAI,KC,KP, 
     *    KB,SURFS,SUBS,TOTPOR,DSURF,DSUB,BFLOW,FLOWC,MONITO, 
     *    SEDLDC,NELEMV,NELEMC,NELEMT,NELEMW,QVUP,QGUP,GVUP,FLOWV,
     *    QSUB,SEDLDV,QCUP,GCUP,QTUP,GTUP,FLOWT,QWUP,GWUP,QCLAT,
     *    SAT,FC,WP,DELTX1,DELTX2,MANN2,CSLOPE,BETA,CWIDTH,PSIZE,
     *    TAWC,KRILL,AET,KSTO,RFLOW,RCOMP,RSTORE,POLY,RNWL,RMWL,
     *    KSPILL,LSPILL,DCULV,EII,EOI,XAREA,KLOSSE,KLOSSC,LCULV,
     *    SCULV,KORIF,NCULV,DELTT,NSEG,FLOW,Z,Lm,LATRAD)             
  36      DO 35 I=1,NGRID         
             L=SOIL(I)         
             ASW(I)=((SURFS(I)+SUBS(I))/
     *       ((DSURF(I)+DSUB(I))*(SAT(L)-WP(L))/10.))*100.
             IF(ASW(I) .GT. EFFPOR(L))THEN
               ASW(I)=EFFPOR(L)
             ENDIF  
             IF(CLAY(L) .GT. 30.)THEN
               TAWC(I)=-0.5-284.*ASW(I)/100.*(ASW(I)/100.-0.3)   
             ELSE
               TAWC(I)=-2.85-(8.87/((VFSAND(L)+0.1)**0.2))-      
     *         (0.16*CACO3(L))+(3.65*SAR(L))+
     *         (3.79*(SPSURF(L)**0.2))+
     *         (28.1*((WDCLAY(L)/CLAY(L))**0.8)*1./(SAND(L)**0.3))
             ENDIF
             CFILT(I)=0.
             CIS(I)=0.
             FILT(I)=0.
             IF(FLOWV(I) .EQ. 0.)THEN
               TEST(I)=0
               QVLAT(I)=0.
             ENDIF
  35      CONTINUE 
          CLOSE (6)
          RETURN
        ENDIF  
  31  CONTINUE
      END
C     .............................................................
C     RUNOFF ROUTINE (2ND ORDER ROUTINE)...........................
      SUBROUTINE RUNOFF(SAT1,WP1,DSUB1,DELTT1,FILTX1,  
     *QGUP1,DELTX1,SUBS1,KSTO1,QSUB1,POND1,FLOWV1,QVLAT1,
     *MANN,SLOPE,QVUP1,BETA,DWBA1,DW1,KTC1,KRILL1,TAWC1,RAIN,GVUP1,
     *KRILL2,FSFAC1,CANEF1,GRCE1,SEDLD1,TEST1,FC1,SURFS1,DSURF1,
     *NOUT,CAF) 

		!DEC$ ATTRIBUTES DLLEXPORT, STDCALL :: RUNOFF
	  !DEC$ ATTRIBUTES ALIAS:'RUNOFF' :: RUNOFF
      !DEC$ ATTRIBUTES REFERENCE :: SAT1,WP1,DSUB1,DELTT1,FILTX1,  
     *QGUP1,DELTX1,SUBS1,KSTO1,QSUB1,POND1,FLOWV1,QVLAT1,
     *MANN,SLOPE,QVUP1,BETA,DWBA1,DW1,KTC1,KRILL1,TAWC1,RAIN,GVUP1,
     *KRILL2,FSFAC1,CANEF1,GRCE1,SEDLD1,TEST1,FC1,SURFS1,DSURF1,
     *NOUT,CAF

	 
      REAL KSTO1,MANN,KTC1,KRILL1,KRILL2,NOUT
      INTEGER TEST1
      TEMP=DSURF1*10.*(FC1-WP1)/100.
      TEMPS=DSUB1*10.*(SAT1-WP1)/100.
      DRAIN=DELTT1*60.*(FILTX1*10./3600 +(1000*QGUP1/(DELTX1**2.)))
      IF((DRAIN .GT. 0.) .AND. (SUBS1 .EQ. TEMPS))THEN
        QSUB1=(QSUB1+QGUP1)*(1-EXP(-DELTT1/KSTO1))       
      ELSE IF((DRAIN .GT. 0.) .AND. ((SUBS1+DRAIN) .GT. TEMPS))THEN
        QSUB1=(QSUB1+QGUP1)*(1-EXP(-DELTT1/KSTO1))       
        SUBS1=TEMPS
      ELSE
        QSUB1=QGUP1*(1-EXP(-DELTT1/KSTO1))
        SUBS1=SUBS1+DELTT1*60.*(FILTX1*10./3600.)
        SURFS1=TEMP
      ENDIF
      IF(TEST1 .EQ. 1)THEN
        Q=FLOWV1
        QVLAT1=(QVLAT1+(POND1*DELTX1/(1000.*NOUT))/(DELTT1*60.))/2.
        ALPHA=(MANN*(DELTX1**0.67)/(SLOPE**0.5))**0.6
        IF(Q .LT. 0)THEN
          Q=0.
        ENDIF
        IF(QVUP1 .LT. 0.)THEN
          QVUP1=0.
        ENDIF 
        IF((Q+QVUP1) .EQ. 0.)THEN
          B=1.
        ELSE
          B=BETA-1.
        ENDIF
        FLOWV1=((DELTT1*60./DELTX1)*QVUP1+                  
     *  ALPHA*BETA*Q*(((Q+QVUP1)/2.)**B)+
     *  DELTT1*60.*QVLAT1)/((DELTT1*60./DELTX1)+
     *  ALPHA*BETA*(((Q+QVUP1)/2.)**B))
        CONST1=(DELTT1*60./DELTX1)*QVUP1+
     *  (ALPHA*(Q**BETA))+ 
     *  DELTT1*60.*QVLAT1
        IF(FLOWV1 .LT. 0.)THEN
          FLOWV1=0.
        ENDIF
        CALL KMATIC(DELTT1,DELTX1,ALPHA,BETA,FLOWV1,CONST1) 
      ELSE
        FLOWV1=0.
      ENDIF
      TAW=(9.81*1000.*DWBA1/((8.*9.81)**(1./3.)*(DW1)**(2./3.)))*
     *(SLOPE*FLOWV1/DELTX1)**(2./3.)
      TC=KTC1*(TAW**1.5)
      DETCAP=KRILL1*(TAW-TAWC1)
      IF((TC .EQ. 0.) .OR. (DETCAP .LE. 0.) .OR.
     *  (GVUP1 .GT. TC))THEN
        RILL=0.
      ELSE
        RILL=DETCAP*(1.-GVUP1/TC)
      ENDIF
C      RILLY=KRILL2*(RAIN**2.)*FSFAC1*CANEF1*GRCE1
C	------------------------NEW EQUATION----------------------
C	CAF - cover adjustment factor
      RILLY=KRILL2*(RAIN*FLOWV1)*FSFAC1*CAF
C	------------------------END OF NEW EQUATION---------------
      DELTG=(RILL+RILLY)/DELTX1
      IF((TC-GVUP1) .GE. DELTG)THEN
        SEDLD1=GVUP1+DELTG
      ELSE
        SEDLD1=TC
      ENDIF
      RETURN
      END
C     .............................................................
C     KINEMATIC ROUTINE (4TH ORDER ROUTINE)........................
      SUBROUTINE KMATIC(DELTT1,DELTX1,ALPHA,BETA,FLOW,CONST1)
	  
	  !DEC$ ATTRIBUTES DLLEXPORT, STDCALL :: KMATIC
	  !DEC$ ATTRIBUTES ALIAS:'KMATIC' :: KMATIC
      !DEC$ ATTRIBUTES REFERENCE :: DELTT1,DELTX1,ALPHA,BETA,FLOW,CONST1
	  
	  
	  
  41  IF(FLOW .GT. 0.)THEN
        VALUE=DELTT1*60./DELTX1
        FFLOW=(VALUE*FLOW)+(ALPHA*(FLOW**BETA))-CONST1
        DFLOW=VALUE+(ALPHA*BETA*FLOW**(BETA-1))
        TEMP=FLOW
        FLOW=FLOW-(FFLOW/DFLOW)
      ENDIF
      IF((FLOW .GT. 0.) .AND. (ABS(FLOW-TEMP) .GT. 0.00005))THEN
        GO TO 41
      ELSE IF(FLOW .LT. 0.)THEN
        FLOW=TEMP
      ENDIF
      RETURN
      END
C     .............................................................
C     SATURATION RUNOFF ROUTINE (2ND ORDER ROUTINE)...............
      SUBROUTINE SATROF(DELTT1,QGUP1,DELTX1,DSUB1,SAT1,WP1,SUBS1, 
     *RAIN,POND1,QVUP1,FILTX1,FLOWV1,QVLAT1,MANN,SLOPE,BETA,FLOWC1,
     *QCLAT1,FLOWW1,DELTX,MANNN,SLOPE1,QCUP1,WIDTH1,DWBA1,DW1,FC1,
     *KTC1,KRILL1,TAWC1,GVUP1,KRILL2,FSFAC1,CANEF1,GRCE1,SEDLD1,
     *GCUP1,PSIZE1,SEDLD2,DEPTH,SURFS1,DSURF1,KSTO1,TEST1,QSUB1,
     *BFLOW,NOUT,CAF) 
	 
	 !DEC$ ATTRIBUTES DLLEXPORT, STDCALL :: SATROF
	  !DEC$ ATTRIBUTES ALIAS:'SATROF' :: SATROF
      !DEC$ ATTRIBUTES REFERENCE :: DELTT1,QGUP1,DELTX1,DSUB1,SAT1,WP1,SUBS1, 
     *RAIN,POND1,QVUP1,FILTX1,FLOWV1,QVLAT1,MANN,SLOPE,BETA,FLOWC1,
     *QCLAT1,FLOWW1,DELTX,MANNN,SLOPE1,QCUP1,WIDTH1,DWBA1,DW1,FC1,
     *KTC1,KRILL1,TAWC1,GVUP1,KRILL2,FSFAC1,CANEF1,GRCE1,SEDLD1,
     *GCUP1,PSIZE1,SEDLD2,DEPTH,SURFS1,DSURF1,KSTO1,TEST1,QSUB1,
     *BFLOW,NOUT,CAF
	 
	 
      REAL MANN,MANNN,KTC1,KRILL1,KRILL2,KSTO1,NOUT,CAF
      INTEGER TEST1
      TEMP=DSURF1*10.*(FC1-WP1)/100.
      TEMPS=DSUB1*10.*(SAT1-WP1)/100.
      DRAIN=DELTT1*60.*(FILTX1*10./3600 +(1000*QGUP1/(DELTX1**2.)))
      IF((DRAIN .GT. 0.) .AND. (SUBS1 .EQ. TEMPS))THEN
        QSUB1=(QSUB1+QGUP1)*(1-EXP(-DELTT1/KSTO1))       
      ELSE IF((DRAIN .GT. 0.) .AND. ((SUBS1+DRAIN) .GT. TEMPS))THEN
        QSUB1=(QSUB1+QGUP1)*(1-EXP(-DELTT1/KSTO1))       
        SUBS1=TEMPS
      ELSE
        QSUB1=QGUP1*(1-EXP(-DELTT1/KSTO1))
        SUBS1=SUBS1+DELTT1*60.*(FILTX1*10./3600.)
        SURFS1=TEMP
      ENDIF
      IF(TEST1 .EQ. 1)THEN
        Q=FLOWV1
        QVLAT1=(QVLAT1+(POND1*DELTX1/(1000.*NOUT))/(DELTT1*60.))/2.
        ALPHA=(MANN*(DELTX1**0.67)/(SLOPE**0.5))**0.6
        IF(Q .LT. 0)THEN
          Q=0.
        ENDIF
        IF(QVUP1 .LT. 0.)THEN
          QVUP1=0.
        ENDIF 
        IF((Q+QVUP1) .EQ. 0.)THEN
          B=1.
        ELSE
          B=BETA-1.
        ENDIF
        FLOWV1=((DELTT1*60./DELTX1)*QVUP1+               
     *  ALPHA*BETA*Q*(((Q+QVUP1)/2.)**B)+
     *  DELTT1*60.*QVLAT1)/((DELTT1*60./DELTX1)+
     *  ALPHA*BETA*(((Q+QVUP1)/2.)**B))
        CONST1=(DELTT1*60./DELTX1)*QVUP1+
     *  (ALPHA*(Q**BETA))+ 
     *  DELTT1*60.*QVLAT1
        IF(FLOWV1 .LT. 0.)THEN
          FLOWV1=0.
        ENDIF
        CALL KMATIC(DELTT1,DELTX1,ALPHA,BETA,FLOWV1,CONST1) 
      ELSE
        FLOWV1=0.
      ENDIF
      Q=FLOWC1
      QCLAT1=(QCLAT1+((FLOWV1+QSUB1+FLOWW1)/DELTX))/2.
      ALPHA=(MANNN*((2.5*WIDTH1)**0.67)/(SLOPE1**0.5))**0.6
      IF(Q .LT. BFLOW)THEN
        Q=BFLOW
      ENDIF
      IF(QCUP1 .LT. BFLOW)THEN
        QCUP1=BFLOW
      ENDIF 
      IF((Q+QCUP1) .EQ. 0.)THEN
        B=1.
      ELSE
        B=BETA-1.
      ENDIF
      FLOWC1=((DELTT1*60./DELTX)*QCUP1+                       
     *ALPHA*BETA*Q*(((Q+QCUP1)/2.)**B)+
     *DELTT1*60.*QCLAT1)/((DELTT1*60./DELTX)+
     *ALPHA*BETA*(((Q+QCUP1)/2.)**B))
      CONST1=(DELTT1*60./DELTX)*QCUP1+(ALPHA*(Q**BETA))+      
     *DELTT1*60.*QCLAT1
      IF(FLOWC1 .LT. BFLOW)THEN
        FLOWC1=BFLOW
      ENDIF
      CALL KMATIC(DELTT1,DELTX,ALPHA,BETA,FLOWC1,CONST1)   
      DEPTH=ALPHA*(FLOWC1**BETA)/WIDTH1
      TAW=(9.81*1000.*DWBA1/((8.*9.81)**(1./3.)*(DW1)**(2./3.)))*
     *(SLOPE*FLOWV1/DELTX1)**(2./3.)
      TC=KTC1*(TAW**1.5)
      DETCAP=KRILL1*(TAW-TAWC1)
      IF((TC .EQ. 0.) .OR. (DETCAP .LE. 0.) .OR.
     *  (GVUP1 .GT. TC))THEN
        RILL=0.
      ELSE
        RILL=DETCAP*(1.-GVUP1/TC)
      ENDIF
C      RILLY=KRILL2*(RAIN**2.)*FSFAC1*CANEF1*GRCE1
C	------------------------NEW EQUATION----------------------
       RILLY=KRILL2*(RAIN*FLOWV1)*FSFAC1*CAF
C	-------------------END OF NEW EQUATION--------------------
      DELTG=(RILL+RILLY)/DELTX1
      IF((TC-GVUP1) .GE. DELTG)THEN
        SEDLD1=GVUP1+DELTG
      ELSE	
        SEDLD1=TC
      ENDIF
      GCUP1=GCUP1+(DELTX1*SEDLD1)
      BSVEL=9.81*DEPTH*SLOPE1
      WIE=(BSVEL**2.)/(1.65*9.81*PSIZE1/1000.)
      EX=0.1*(WIE**(5./2.))
      IF((FLOWC1 .EQ. 0.) .OR. (DEPTH .EQ. 0. ))THEN
        BFFAC=0.
        FIE=0.
      ELSE
        BFFAC=8.*BSVEL/((FLOWC1/(DEPTH*WIDTH1))**2.)
        FIE=4.*EX/BFFAC
      ENDIF
      QS=FIE*(1000.*((1.65*9.81*PSIZE1)**(3./2.)))
      CTC=QS*WIDTH1/1000.
      TAWS=1000.*9.81*DEPTH*SLOPE1
      IF((CTC .GT. GCUP1) .AND. (TAWS .GT. TAWC1))THEN
        EROS=WIDTH1*KRILL1*(TAWS-TAWC1)*(1.-GCUP1/CTC)
        SEDLD2=(GCUP1+EROS)
      ELSE IF(CTC .LT. GCUP1)THEN  
        SEDLD2=CTC   
      ELSE
        SEDLD2=GCUP1
      ENDIF
      RETURN
      END
C     .............................................................
C     STRUCTURES ROUTINE (2ND ORDER ROUTINE).......................
      SUBROUTINE STRUCT(FLOWT1,QTLAT1,FLOWV1,DELTX1,MANNN,TSLOPE, 
     *DELTT1,QTUP1,BETA,TWIDTH,GTUP1,SEDLD1,PSIZE1,TAWC1,KRILL1,
     *SEDLD2,STRUC,RI,FLOWW1,QWLAT1,DELTX,WSLOPE,QWUP1,WWIDTH,
     *GWUP1,SEDLD3)
	 
	 
	 !DEC$ ATTRIBUTES DLLEXPORT, STDCALL :: STRUCT
      !DEC$ ATTRIBUTES ALIAS:'STRUCT' :: STRUCT
      !DEC$ ATTRIBUTES REFERENCE :: FLOWT1,QTLAT1,FLOWV1,DELTX1,MANNN,TSLOPE, 
     *DELTT1,QTUP1,BETA,TWIDTH,GTUP1,SEDLD1,PSIZE1,TAWC1,KRILL1,
     *SEDLD2,STRUC,RI,FLOWW1,QWLAT1,DELTX,WSLOPE,QWUP1,WWIDTH,
     *GWUP1,SEDLD3
	 
	 
      REAL MANNN,KRILL1
      INTEGER STRUC
      Q=FLOWT1
      QTLAT1=(QTLAT1+(FLOWV1/DELTX1))/2.
      ALPHA=(MANNN*((2.5*TWIDTH)**0.67)/(TSLOPE**0.5))**0.6
      IF(Q .LT. 0)THEN
        Q=0.
      ENDIF
      IF(QTUP1 .LT. 0.)THEN
        QTUP1=0.
      ENDIF 
      IF((Q+QTUP1) .EQ. 0.)THEN
        B=1.
      ELSE
        B=BETA-1.
      ENDIF
      FLOWT1=((DELTT1*60./DELTX1)*QTUP1+                       
     *ALPHA*BETA*Q*(((Q+QTUP1)/2.)**B)+
     *DELTT1*60.*QTLAT1)/((DELTT1*60./DELTX1)+
     *ALPHA*BETA*(((Q+QTUP1)/2.)**B))
      CONST1=(DELTT1*60./DELTX1)*QTUP1+(ALPHA*Q**BETA)+      
     *DELTT1*60.*QTLAT1
      IF(FLOWT1 .LT. 0.)THEN
        FLOWT1=0.
      ENDIF
      CALL KMATIC(DELTT1,DELTX1,ALPHA,BETA,FLOWT1,CONST1)  
      DEPTH1=ALPHA*(FLOWT1**BETA)/TWIDTH
      GTUP1=GTUP1+(DELTX1*SEDLD1)
      BSVEL=9.81*DEPTH1*TSLOPE
      WIE=(BSVEL**2.)/(1.65*9.81*PSIZE1/1000.)
      EX=0.1*(WIE**(5./2.))
      IF((FLOWT1 .EQ. 0.) .OR. (DEPTH1 .EQ. 0. ))THEN
        BFFAC=0.
        FIE=0.
      ELSE
        BFFAC=8.*BSVEL/((FLOWT1/(DEPTH1*TWIDTH))**2.)
        FIE=4.*EX/BFFAC
      ENDIF
      QS=FIE*(1000.*((1.65*9.81*PSIZE1)**(3./2.)))
      CTC=QS*TWIDTH/1000.
      TAWS=1000.*9.81*DEPTH1*TSLOPE
      IF((CTC .GT. GTUP1) .AND. (TAWS .GT. TAWC1))THEN
        EROS=TWIDTH*KRILL1*(TAWS-TAWC1)*(1.-GTUP1/CTC)
        SEDLD2=(GTUP1+EROS)
      ELSE IF(CTC .LT. GTUP1)THEN  
        SEDLD2=CTC   
      ELSE
        SEDLD2=GTUP1
      ENDIF
      IF(STRUC .EQ. 1)THEN
        CALL GWWAY(RI,FLOWW1,QWLAT1,FLOWT1,DELTX,WSLOPE,DELTT1,
     *  QWUP1,BETA,WWIDTH,GWUP1,SEDLD2,PSIZE1,SEDLD3)             
      ELSE
        FLOWW1=0.
        SEDLD3=0.
      ENDIF
      RETURN
      END
C     .............................................................
C     GRASS WATERWAY ROUTINE (3RD ORDER ROUTINE)...................
      SUBROUTINE GWWAY(RI,FLOWW1,QWLAT1,FLOWT1,DELTX,WSLOPE,DELTT1,
     *QWUP1,BETA,WWIDTH,GWUP1,SEDLD2,PSIZE1,SEDLD3)
	 
	 
	 !DEC$ ATTRIBUTES DLLEXPORT, STDCALL :: GWWAY
	  !DEC$ ATTRIBUTES ALIAS:'GWWAY' :: GWWAY
      !DEC$ ATTRIBUTES REFERENCE :: RI,FLOWW1,QWLAT1,FLOWT1,DELTX,WSLOPE,DELTT1,
     *QWUP1,BETA,WWIDTH,GWUP1,SEDLD2,PSIZE1,SEDLD3
	 
	 
      REAL MANN
      IF((RI .GE. 1.) .AND. (RI .LE. 3.))THEN
        MANN=0.035
      ELSE IF((RI .GT. 3.) .AND. (RI .LE. 5.))THEN
        MANN=0.05
      ELSE IF((RI .GT. 5.) .AND. (RI .LE. 7.))THEN
        MANN=0.065
      ELSE IF((RI .GT. 7.) .AND. (RI .LE. 9.))THEN
        MANN=0.13
      ELSE
        MANN=0.2
      ENDIF
  51  Q=FLOWW1
      QWLAT1=(QWLAT1+(FLOWT1/DELTX))/2.
      ALPHA=(MANN*((2.5*WWIDTH)**0.67)/(WSLOPE**0.5))**0.6
      IF(Q .LT. 0)THEN
        Q=0.
      ENDIF
      IF(QWUP1 .LT. 0.)THEN
        QWUP1=0.
      ENDIF 
      IF((Q+QWUP1) .EQ. 0.)THEN
        B=1.
      ELSE
        B=BETA-1.
      ENDIF
      FLOWW1=((DELTT1*60./DELTX)*QWUP1+                       
     *ALPHA*BETA*Q*(((Q+QWUP1)/2.)**B)+
     *DELTT1*60.*QWLAT1)/((DELTT1*60./DELTX)+
     *ALPHA*BETA*(((Q+QWUP1)/2.)**B))
      CONST1=(DELTT1*60./DELTX)*QWUP1+(ALPHA*(Q**BETA))+      
     *DELTT1*60.*QWLAT1
      IF(FLOWW1 .LT. 0.)THEN
        FLOWW1=0.
      ENDIF
      CALL KMATIC(DELTT1,DELTX,ALPHA,BETA,FLOWW1,CONST1)   
      DEPTH1=ALPHA*(FLOWW1**BETA)/WWIDTH
      IF(DEPTH1 .EQ. 0.)THEN
        GO TO 52
      ENDIF
      VEL=(DEPTH1**(2./3.))*(WSLOPE**0.5)/MANN
      X=13.3-1.4*RI
      TEMP=1./(2.1+2.3*X+6.*ALOG(10.8*VEL*DEPTH1))
      IF((TEMP .GT. 0.) .AND. 
     *  (ABS((TEMP-MANN)/MANN) .GT. 0.05))THEN 
        MANN=TEMP
        GO TO 51
      ENDIF
  52  GWUP1=GWUP1+SEDLD2
      BSVEL=9.81*DEPTH1*WSLOPE
      WIE=(BSVEL**2.)/(1.65*9.81*PSIZE1/1000.)
      EX=0.1*(WIE**(5./2.))
      IF((FLOWW1 .EQ. 0.) .OR. (DEPTH1 .EQ. 0. ))THEN
        BFFAC=0.
        FIE=0.
      ELSE
        BFFAC=8.*BSVEL/((FLOWW1/(DEPTH1*WWIDTH))**2.)
        FIE=4.*EX/BFFAC
      ENDIF
      QS=FIE*(1000.*((1.65*9.81*PSIZE1)**(3./2.)))
      CTC=QS*WWIDTH/1000.
      IF(CTC .LT. GWUP1)THEN
        SEDLD3=CTC
      ELSE   
        SEDLD3=GWUP1   
      ENDIF
      RETURN                 
      END
C     .............................................................
C     HYDRAULIC STRUCTURES ROUTINE (3RD ORDER ROUTINE).............
      SUBROUTINE HSTRUC(FLOWC1,RFLOW1,RCOMP1,RSTOR,POLY,RNWL,RMWL,
     *KSPILL,LSPILL,SEDLD2,DELTT,FLOW) 
	
	!DEC$ ATTRIBUTES DLLEXPORT, STDCALL :: HSTRUC
	  !DEC$ ATTRIBUTES ALIAS:'HSTRUC' :: HSTRUC
      !DEC$ ATTRIBUTES REFERENCE :: FLOWC1,RFLOW1,RCOMP1,RSTOR,POLY,RNWL,RMWL,
     *KSPILL,LSPILL,SEDLD2,DELTT,FLOW

	 
      REAL POLY(4),KSPILL,LSPILL
      RFLOW2=FLOWC1
      TEMP1=RFLOW2+RFLOW1+RCOMP1
      A1=POLY(3)/POLY(4)
      A2=POLY(2)/POLY(4)
      A3=(POLY(1)-RSTOR)/POLY(4)
      Q=((3.*A2)-(A1**2))/9.
      R=((9.*A1*A2)-(27.*A3)-(2.*(A1**3)))/54.
      IF(((Q**3)+(R**2)) .GE. 0.)THEN
        C=((Q**3)+(R**2))**0.5
      ELSE
        C=0.
      ENDIF
      IF((R+C) .GE. 0.)THEN 
        S=(R+C)**(1./3.)
      ELSE
        S=0.
      ENDIF
      IF((R-C) .GE. 0.)THEN
        T=(R-C)**(1./3.)
      ELSE
        T=0.
      ENDIF
      RELEV=S+T-(A1/3.)       
      IF((RELEV .GT. RNWL) .AND. (RELEV .LT. RMWL))THEN
        RHEAD=RELEV-RNWL
        FLOWC1=KSPILL*((2.*9.81)**0.5)*LSPILL*(RHEAD**1.5)
        RSTOR=(TEMP1-FLOWC1)*DELTT*60./2.
        RCOMP1=TEMP1-(2.*FLOWC1)
        STORE=POLY(1)+(POLY(2)*RNWL)+(POLY(3)*RNWL**2)+
     *  (POLY(4)*RNWL**3)
        IF(RCOMP1 .LT. (2.*STORE/(DELTT*60.)))THEN
          RCOMP1=(2.*STORE/(DELTT*60.))
          RSTOR=STORE
          FLOWC1=(TEMP1-RCOMP1)/2.
        ENDIF
      ELSE IF(RELEV .GE. RMWL)THEN
        RHEAD=RMWL-RNWL   
        FLOWC1=KSPILL*((2.*9.81)**0.5)*LSPILL*(RHEAD**1.5)
        RSTOR=(TEMP1-FLOWC1)*DELTT*60./2.
        RCOMP1=TEMP1-(2.*FLOWC1)
        STORE=POLY(1)+(POLY(2)*RNWL)+(POLY(3)*RNWL**2)+
     *  (POLY(4)*RNWL**3)
        IF(RCOMP1 .LT. (2.*STORE/(DELTT*60.)))THEN
          RCOMP1=(2.*STORE/(DELTT*60.))
          RSTOR=STORE
          FLOWC1=(TEMP1-RCOMP1)/2.
        ENDIF
      ELSE
        IF(INT(RELEV+0.05) .EQ. INT(RNWL))THEN
          FLOWC1=FLOW
        ELSE
          FLOWC1=0.
        ENDIF
        RSTOR=(TEMP1-FLOWC1)*DELTT*60./2.
        IF(RSTOR .LE. 0.)THEN
          RSTOR=0.
          RCOMP1=0.
          FLOWC1=0.
        ELSE
          RCOMP1=TEMP1-(2.*FLOWC1)
          IF(RCOMP1 .LT. 0.)THEN
            RCOMP1=0.
            FLOWC1=0.
          ENDIF
        ENDIF
      ENDIF
      RFLOW1=RFLOW2
      FLOW=FLOWC1
      SEDLD2=0.
      RETURN
      END 
C     .............................................................
C     CULVERT ROUTINE (4TH ORDER ROUTINE)..........................
      SUBROUTINE CULVT(DEPTH1,DCULV,EII,EOI,XAREA,KLOSSE,KLOSSC,  
     *LCULV,SCULV,KORIF,NCULV,FLOWC1) 


		!DEC$ ATTRIBUTES DLLEXPORT, STDCALL :: CULVT
	  !DEC$ ATTRIBUTES ALIAS:'CULVT' :: CULVT
      !DEC$ ATTRIBUTES REFERENCE :: DEPTH1,DCULV,EII,EOI,XAREA,KLOSSE,KLOSSC,  
     *LCULV,SCULV,KORIF,NCULV,FLOWC1

	 
      REAL NCULV,LCULV,KLOSSE,KLOSSC,KORIF
      IF(DEPTH1 .GT. DCULV)THEN
        PHEAD=(EII+DEPTH1)-(EOI +0.6*DCULV)
        IF(PHEAD .LT. 0.)THEN
          PHEAD=0.
        ENDIF
        FLOWC1=(XAREA*(19.62*PHEAD)**0.5)/(1.+KLOSSE+KLOSSC*LCULV)
        SN=(KLOSSC*(FLOWC1/XAREA)**2.)/19.62      
        IF(SCULV .GT. SN)THEN
        FLOWC1=XAREA*KORIF*(19.62*(DEPTH1-0.5*DCULV))**0.5
        ENDIF 
      ELSE IF((DEPTH1 .GT. 0.) .AND. 
     *  (DCULV*DEPTH1 .GE. DEPTH1**2.))THEN
 61     A=((DCULV/2.)**2.)*ACOS(((DCULV/2.)-DEPTH1)/(DCULV/2.))
        B=-(((DCULV/2.)-DEPTH1)*(DCULV*DEPTH1-DEPTH1**2.)**0.5)
        XAREA1=A+B
        C=DCULV*ACOS(((DCULV/2.)-DEPTH1)/(DCULV/2.))
        IF((XAREA1 .LE. 0.) .OR. (C .LE. 0.))THEN
          HYDR=0.
        ELSE 
          HYDR=XAREA1/C
        ENDIF
        VEL=(HYDR**(2./3.))*(SCULV**0.5)/NCULV   
        STLOSS=(VEL**2.)/19.62       
        DEPTH1=DEPTH1-STLOSS 
        FLOWC1=XAREA1*VEL
        IF((DEPTH1 .GT. 0.) .AND. ((STLOSS/DEPTH1) .GT. 0.05))THEN 
          GO TO 61
        ENDIF
      ELSE
        FLOWC1=0. 
      ENDIF
      RETURN
      END
C     .............................................................
C     BASEFLOW ET ROUTINE (2ND ORDER ROUTINE)......................
      SUBROUTINE BASEFL(ETT,STIME,NGRID,VEG,SOIL,DUAL,     
     *STRUC1,STRUC3,SEDLDT,FLOWW,SEDLDW,GI,PCAN,LAI,KC,KP,  
     *KB,SURFS,SUBS,TOTPOR,DSURF,DSUB,BFLOW,FLOWC,MONITO,    
     *SEDLDC,NELEMV,NELEMC,NELEMT,NELEMW,QVUP,QGUP,GVUP,FLOWV,QSUB,
     *SEDLDV,QCUP,GCUP,QTUP,GTUP,FLOWT,QWUP,GWUP,
     *QCLAT,SAT,FC,WP,DELTX1,DELTX2,MANN2,CSLOPE,BETA,
     *CWIDTH,PSIZE,TAWC,KRILL,AET,KSTO,RFLOW,RCOMP,RSTORE,        
     *POLY,RNWL,RMWL,KSPILL,LSPILL,DCULV,EII,EOI,XAREA,KLOSSE,
     *KLOSSC,LCULV,SCULV,KORIF,NCULV,DELTT,NSEG,FLOW,Z,Lm,LATRAD) 
	 
	 
	 !DEC$ ATTRIBUTES DLLEXPORT, STDCALL :: BASEFL
	  !DEC$ ATTRIBUTES ALIAS:'BASEFL' :: BASEFL
      !DEC$ ATTRIBUTES REFERENCE :: ETT,STIME,NGRID,VEG,SOIL,DUAL,     
     *STRUC1,STRUC3,SEDLDT,FLOWW,SEDLDW,GI,PCAN,LAI,KC,KP,  
     *KB,SURFS,SUBS,TOTPOR,DSURF,DSUB,BFLOW,FLOWC,MONITO,    
     *SEDLDC,NELEMV,NELEMC,NELEMT,NELEMW,QVUP,QGUP,GVUP,FLOWV,QSUB,
     *SEDLDV,QCUP,GCUP,QTUP,GTUP,FLOWT,QWUP,GWUP,
     *QCLAT,SAT,FC,WP,DELTX1,DELTX2,MANN2,CSLOPE,BETA,
     *CWIDTH,PSIZE,TAWC,KRILL,AET,KSTO,RFLOW,RCOMP,RSTORE,        
     *POLY,RNWL,RMWL,KSPILL,LSPILL,DCULV,EII,EOI,XAREA,KLOSSE,
     *KLOSSC,LCULV,SCULV,KORIF,NCULV,DELTT,NSEG,FLOW,Z,Lm,LATRAD
	 
	 
	 
      IMPLICIT REAL(K)
C     SOIL TYPE PARAMETERS .......................................
      REAL PSIZE(3),FC(3),WP(3),SAT(3),TOTPOR(3)
C     VEGETATION TYPE PARAMETERS .................................
      REAL KC(8),LAI(8)
C     SYSTEM AND WATERSHED PARAMETERS .............................
      REAL GI(73),DELTX2(73),MANN2(73),PCAN(73),KSTO(73),
     *DSUB(73),DSURF(73),KRILL(73),TAWC(73),CWIDTH(73),      
     *POLY(4),LCULV,NCULV,LSPILL,CSLOPE(73)
      INTEGER SOIL(73),VEG(73),DUAL(73),STRUC1(73),STRUC3(73),
     *NELEMV(73),NELEMC(73),NELEMT(73),NELEMW(73)   
C     SYSTEM AND STATE VARIABLES ................................. 
      REAL AET(73),AVPOR(73),FLOWV(73),FLOWC(73),FLOWT(73),  
     *FLOWW(73),QCLAT(73),QGUP(73),QVUP(73),QCUP(73),
     *QTUP(73),QWUP(73),SEDLDV(73),SEDLDC(73),SEDLDT(73), 
     *SEDLDW(73),SUBS(73),SURFS(73),QSUB(73),RFLOW(73), 
     *RSTORE(73),GVUP(73),GCUP(73),GTUP(73),GWUP(73),
     *RCOMP(73)
C	New variables		N=daylight hours=DAYLHR, n=actual duration
C						of sunshine = ADS,d=Solar Declination SOLARD
C						dr=inverse relative distance earth-sun=IRDES
C						LAT = latitude,Ws=Sunset hour angle=SHA
C						Gsc= Solar constant
	REAL EA,ES,G,RN,RA,DELTA,Pa,Cp,GAMMA,LAMBDA,Zm,Zo,Sn,Ln,
     *RS,Patm,T,ROW,Rsp,KT,U2,Z(73),Rso,ALFA,Tmax,Tmin,Ts,So,RHmean,
     *PW,Eomax,Eomin,DAYLHR,ADS,WS,Rns,RnL,Gsc,SOLARD,LATRAD,IRDES,
     *SHA,W1,W2,W,Lz,Lm,SC,T1,B
	INTEGER JAY,DAY,MONTH,YEAR
      NSIM=INT(ETT/30)
      IF(NSIM .LT. 1)THEN
        GO TO 66
      ENDIF	
      DO 62 J=1,NSIM
C     READ (5,FMT='(7X,F24.2)',END=300,ERR=300) EPAN
	READ(5,*,END=300,ERR=300)MONTH,DAY,YEAR,RS,Tmin,Tmax,RHmean,Patm,
     * U2	
         STIME=STIME+30. 
         DO 63 I=1,NGRID
            L=VEG(I)
            L1=SOIL(I)
            IF(DUAL(I) .EQ. 1)THEN     
              CALL ETTWO(SUBS(I),SURFS(I),FLOWC(I),QCLAT(I),
     *        SEDLDC(I),QGUP(I),DELTX1,DSUB(I),DSURF(I),FC(L1),
     *        WP(L1),DELTX2(I),MANN2(I),CSLOPE(I),QCUP(I),BETA,
     *        CWIDTH(I),PSIZE(L1),GCUP(I),TAWC(I),KRILL(I),DEPTH,
     *        KSTO(I),QSUB(I),DELTT,FLOWV(I),SEDLDV(I),BFLOW)
            ELSE
              CALL ETONE(SURFS(I),DSURF(I),FC(L1),WP(L1),
     *        QGUP(I),DELTX1,SUBS(I),DSUB(I),KSTO(I),QSUB(I),
     *        FLOWV(I),SEDLDV(I),DELTT)              
            ENDIF
            IF((STRUC1(I) .EQ. 1) .AND. (DUAL(I) .EQ. 1) .AND.
     *        (STRUC3(I) .EQ. 1))THEN
              CALL HSTRUC(FLOWC(I),RFLOW(I),RCOMP(I),RSTORE(I),
     *        POLY,RNWL,RMWL,KSPILL,LSPILL,SEDLDC(I),DELTT,FLOW)
            ELSE IF((STRUC1(I) .EQ. 1) .AND. (DUAL(I) .EQ. 1))THEN
              CALL CULVT(DEPTH,DCULV,EII,EOI,XAREA,KLOSSE,KLOSSC, 
     *        LCULV,SCULV,KORIF,NCULV,FLOWC(I))                   
            ELSE IF(STRUC1(I) .EQ. 1)THEN
              FLOWV(I)=0.
              SEDLDV(I)=0.
              FLOWT(I)=0.
              SEDLDT(I)=0
              FLOWW(I)=0.
              SEDLDW(I)=0.
            ENDIF
C            PET=(GI(I)*PCAN(I)*LAI(L)*KC(L)*KP*EPAN)+((1.-  
C     *      PCAN(I))*KB*KP*EPAN)
C	----------NEW EQUATIONS AND VARIABLES    -----------------
			 Gsc=0.0820
			 Zm=2
			 Lz=75.			 
			 T=(Tmin+Tmax)/2.				 		 
			 T1=0.17
			 JAY=INT((275.*MONTH/9.)-30.+DAY)-2
			 IF(MONTH.LT. 3)THEN
				JAY=JAY+2
			 ENDIF
			 IF((4*INT(YEAR/4.)) .EQ. YEAR .AND. MONTH .GT. 2) THEN
				JAY=JAY+1
	         ENDIF				 		 
			 IRDES=1+0.033*COS(3.141593/365*JAY)
			 SOLARD=0.409*SIN(2*3.141593/365*JAY-1.39)
			 B=2*3.141593*(JAY-81)/364.
			 SC=0.1645*SIN(2*B)-0.1255*COS(B)-0.025*SIN(B)
			 LATRAD=3.141593/180.*LATRAD
			 SHA=ACOS(-TAN(LATRAD)*TAN(SOLARD))		     
			 GAMMA=0.665*10**(-3)/Patm
			 DELTA=4098*(0.6108*EXP(17.27*T/(T+237.3)))/(T+237.3)**2.
			 EA=0.6108*EXP(17.27*T/(T+237.3))*RHmean/100.			 
			 W=3.141593/12*((DELTT+0.06667*(Lz-Lm)+SC)-12)
			 W1=W-3.141593*T1/24.
			 W2=W+3.141593*T1/24.

			 RA=12.*60/3.141593*Gsc*IRDES*((W2-W1)*SIN(LATRAD)*			 
     *		 SIN(SOLARD)+COS(LATRAD)*COS(SOLARD)*(SIN(W2)-SIN(W1)))
		     Rso=((0.75+2*Z(I))/100000.)*RA			 
			 Rns=(1-0.23)*RS
	         RnL=(((4.903*10**(-9))*(Tmax+273.16)**4)+
     *		 ((4.903*10**(-9))*(Tmin+273.16)**4))/2.*
     *	     (0.34-0.14*EA**0.5)*(1.35*RS/Rso-0.35)
			 RN=Rns-RnL
			 G=0.1*RN			 			 			 
			 PET=(0.408*DELTA*(RN-G)+(GAMMA*37/(T+273))*U2*((0.6108*
     *		 EXP(17.27*T/(T+237.3)))-EA))/(DELTA+GAMMA*(1+0.34*U2))				 
C	---------------- END-------------------------------------               
            IF(SURFS(I) .GE. PET)THEN
              AET(I)=PET
              SURFS(I)=SURFS(I)-PET
            ELSE                                                
              H=(SURFS(I)+SUBS(I))/    
     *        ((DSURF(I)+DSUB(I))*10.*(SAT(L1)-WP(L1))/100.)
              IF((H .GT. 0.) .AND. (H .LE. 1.))THEN
                AVPOR(I)=TOTPOR(L1)*(1.-H)
              ELSE IF(H .GT. 1.)THEN
                AVPOR(I)=0.
              ELSE
                AVPOR(I)=TOTPOR(L1)
              ENDIF
C             AET(I)=PET*((TOTPOR(L1)-AVPOR(I))/TOTPOR(L1))**0.1
C	------------------------NEW EQUATION----------------------
			AET(I)=KC(L)*PET
C	------------------------END OF NEW EQUATION---------------
              SURFS(I)=(SURFS(I)+SUBS(I)-AET(I))*DSURF(I)/
     *        (DSURF(I)+DSUB(I))
              SUBS(I)=(SURFS(I)+SUBS(I)-AET(I))*DSUB(I)/
     *        (DSURF(I)+DSUB(I))
              IF(SURFS(I) .LT. 0.)THEN
                SURFS(I)=0.
              ENDIF
              IF(SUBS(I) .LT. 0.)THEN
                SUBS(I)=0.
              ENDIF  
            ENDIF
  63     CONTINUE 
         DO 64 I=1,NGRID
            L=NELEMV(I)         
            L1=NELEMC(I)
            L2=NELEMT(I)
            L3=NELEMW(I)
            IF(L .EQ. 0)THEN
              QVUP(I)=0.
              QGUP(I)=0.
              GVUP(I)=0.
            ELSE IF(L .GT. NGRID)THEN
              L=L-NGRID
              QVUP(I)=FLOWV(L)+FLOWV(L2) 
              QGUP(I)=QSUB(L)+QSUB(L2)
              GVUP(I)=SEDLDV(L)+SEDLDV(L2)  
            ELSE
              QVUP(I)=FLOWV(L)
              QGUP(I)=QSUB(L)
              GVUP(I)=SEDLDV(L)
            ENDIF    
            IF(L1 .EQ. 0)THEN
              IF(I .EQ. NSEG)THEN
                QCUP(I)=BFLOW
                GCUP(I)=0.
              ELSE
                QCUP(I)=0.
                GCUP(I)=0.
              ENDIF
            ELSE IF(L1 .GT. NGRID)THEN
              L1=L1-NGRID
              QCUP(I)=FLOWC(L1)+FLOWC(L3)-BFLOW
              GCUP(I)=SEDLDC(L1)+SEDLDC(L3)
            ELSE
              QCUP(I)=FLOWC(L1)
              GCUP(I)=SEDLDC(L1)
            ENDIF    
            IF(L2 .EQ. 0)THEN
              QTUP(I)=0.
              GTUP(I)=0.
            ELSE
              QTUP(I)=FLOWT(L2)
              GTUP(I)=SEDLDT(L2)
            ENDIF    
            IF(L3 .EQ. 0)THEN
              QWUP(I)=0.
              GWUP(I)=0.
            ELSE
              QWUP(I)=FLOWW(L3)
              GWUP(I)=SEDLDW(L3)
            ENDIF    
  64     CONTINUE 
         IF((2*INT(J/2.)) .EQ. J)THEN
           IF(FLOWC(MONITO) .GT. 0.)THEN 
             CONC=SEDLDC(MONITO)/FLOWC(MONITO)
           ELSE
             CONC=0.
           ENDIF 
           WRITE (6,FMT='(5X,I8,2X,F13.2,5X,F13.4,5X,I8)')        
     *     INT(STIME),FLOWC(MONITO)*1000,SEDLDC(MONITO),
     *     INT(CONC*1000)
           WRITE (*,FMT='(5X,I8,2X,F13.2,5X,F13.4,5X,I8)')        
     *     INT(STIME),FLOWC(MONITO)*1000,SEDLDC(MONITO),
     *     INT(CONC*1000)
         ENDIF
  62  CONTINUE 
  66  RETURN
 300  WRITE (*,*) 'PAN EVAPORATION DATA INPUT ERROR'
      STOP 300
      END
C     ........................................................
C     ET ROUTINE I (2ND ORDER ROUTINE) ........................
      SUBROUTINE ETONE(SURFS1,DSURF1,FC1,WP1,QGUP1,DELTX1,   
     *SUBS1,DSUB1,KSTO1,QSUB1,FLOWV1,SEDLD1,DELTT)  

		
		!DEC$ ATTRIBUTES DLLEXPORT, STDCALL :: ETONE
	  !DEC$ ATTRIBUTES ALIAS:'ETONE' :: ETONE
      !DEC$ ATTRIBUTES REFERENCE :: SURFS1,DSURF1,FC1,WP1,QGUP1,DELTX1,   
     *SUBS1,DSUB1,KSTO1,QSUB1,FLOWV1,SEDLD1,DELTT
		
	 
      REAL KSTO1 
      TEMPS=DSURF1*10.*(FC1-WP1)/100.
      IF(SURFS1 .GT. TEMPS)THEN
        DRAIN=SURFS1-TEMPS
        SURFS1=TEMPS
      ELSE
        DRAIN=0.
      ENDIF
      TEMPS=DSUB1*10.*(FC1-WP1)/100.
      TT=30./DELTT
      IF(SUBS1 .GT. TEMPS)THEN
        QSUB1=(QGUP1+(DRAIN/1000.)*(DELTX1**2.)/1800.)*
     *  (1-EXP(-30./(KSTO1*TT)))
        SUBS1=TEMPS
      ELSE IF(QSUB1 .GT. 0.)THEN 
        QSUB1=QSUB1*EXP(-30./(KSTO1*TT))       
      ELSE
        QSUB1=0.
      ENDIF
      FLOWV1=0.
      SEDLD1=0.
      RETURN
      END
C     ........................................................
C     ET ROUTINE II (2ND ORDER ROUTINE) ........................
      SUBROUTINE ETTWO(SUBS1,SURFS1,FLOWC1,QCLAT1,SEDLD2,QGUP1,   
     *DELTX1,DSUB1,DSURF1,FC1,WP1,DELTX,MANNN,SLOPE1,QCUP1,BETA,
     *WIDTH1,PSIZE1,GCUP1,TAWC1,KRILL1,DEPTH1,KSTO1,QSUB1,DELTT,
     *FLOWV1,SEDLD1,BFLOW) 


	  !DEC$ ATTRIBUTES DLLEXPORT, STDCALL :: ETTWO
	  !DEC$ ATTRIBUTES ALIAS:'ETTWO' :: ETTWO
      !DEC$ ATTRIBUTES REFERENCE :: SUBS1,SURFS1,FLOWC1,QCLAT1,SEDLD2,QGUP1,   
     *DELTX1,DSUB1,DSURF1,FC1,WP1,DELTX,MANNN,SLOPE1,QCUP1,BETA,
     *WIDTH1,PSIZE1,GCUP1,TAWC1,KRILL1,DEPTH1,KSTO1,QSUB1,DELTT,
     *FLOWV1,SEDLD1,BFLOW
	 
	 
      REAL MANNN,KRILL1,KSTO1
      TEMPS=DSURF1*10.*(FC1-WP1)/100.
      IF(SURFS1 .GT. TEMPS)THEN
        DRAIN=SURFS1-TEMPS
        SURFS1=TEMPS
      ELSE
        DRAIN=0.
      ENDIF
      TEMPS=DSUB1*10.*(FC1-WP1)/100.
      TT=30./DELTT
      IF(SUBS1 .GT. TEMPS)THEN
        QSUB1=(QGUP1+(DRAIN/1000.)*(DELTX1**2.)/1800.)*
     *  (1-EXP(-30./(KSTO1*TT)))
        SUBS1=TEMPS
      ELSE IF(QSUB1 .GT. 0.)THEN 
        QSUB1=QSUB1*EXP(-30./(KSTO1*TT))       
      ELSE
        QSUB1=0.
      ENDIF
      FLOWV1=0.
      SEDLD1=0.
      Q=FLOWC1
      QCLAT1=(QCLAT1+(QSUB1/DELTX))/2.
      ALPHA=(MANNN*((2.5*WIDTH1)**0.67)/(SLOPE1**0.5))**0.6
      IF(Q .LT. BFLOW)THEN
        Q=BFLOW
      ENDIF
      IF(QCUP1 .LT. BFLOW)THEN
        QCUP1=BFLOW
      ENDIF 
      IF((Q+QCUP1) .EQ. 0.)THEN
        B=1.
      ELSE
        B=BETA-1.
      ENDIF
      FLOWC1=((1800./DELTX)*QCUP1+                       
     *ALPHA*BETA*Q*(((Q+QCUP1)/2.)**B)+
     *1800.*QCLAT1)/((1800./DELTX)+
     *ALPHA*BETA*(((Q+QCUP1)/2.)**B))
      CONST1=(1800./DELTX)*QCUP1+(ALPHA*(Q**BETA))+
     *1800.*QCLAT1
      IF(FLOWC1 .LT. BFLOW)THEN
        FLOWC1=BFLOW
      ENDIF
      DELTT1=30.
      CALL KMATIC(DELTT1,DELTX,ALPHA,BETA,FLOWC1,CONST1)
      DEPTH1=ALPHA*(FLOWC1**BETA)/WIDTH1
      BSVEL=9.81*DEPTH1*SLOPE1
      WIE=(BSVEL**2.)/(1.65*9.81*PSIZE1/1000.)
      EX=0.1*(WIE**(5./2.))
      IF((FLOWC1 .EQ. 0.) .OR. (DEPTH1 .EQ. 0. ))THEN
        BFFAC=0.
        FIE=0.
      ELSE
        BFFAC=8.*BSVEL/((FLOWC1/(DEPTH1*WIDTH1))**2.)
        FIE=4.*EX/BFFAC
      ENDIF
      QS=FIE*(1000.*((1.65*9.81*PSIZE1/1000.)**(3./2.)))
      CTC=QS*WIDTH1/1000.
      TAWS=1000.*9.81*DEPTH1*SLOPE1
      IF(CTC .EQ. 0.)THEN
         SEDLD2=0.
      ELSE IF((CTC .GT. GCUP1) .AND. (TAWS .GT. TAWC1))THEN
        EROS=WIDTH1*KRILL1*(TAWS-TAWC1)*(1.-GCUP1/CTC)
        SEDLD2=(GCUP1+EROS)
      ELSE IF(CTC .LT. GCUP1)THEN  
        SEDLD2=CTC   
      ELSE
        SEDLD2=GCUP1
      ENDIF
      RETURN
      END

	SUBROUTINE DELOUT()		
		OPEN(1,FILE='output.dat',STATUS='UNKNOWN')		
		CLOSE(1, STATUS='DELETE')				
	END
