C     -----------------
      SUBROUTINE WEVAP1
C     -----------------
C
C     SUBROUTINE AUTHOR:
C     TIM DURBIN
C     TIMOTHY J. DURBIN, INC.
C     5330 PRIMROSE DRIVE, SUITE 228
C     FAIR OAKS, CA 95628
C     (916) 536-2314
C
C     PURPOSE:
C     THIS SUBROUTINE IS USED TO SIMULATE THE DISCHARGE OF GROUNDWATER
C     FROM A SHALLOW WATER TABLE AND IS ASSUMED TO BE LINERALY RELATED
C     TO THE DEPTH BELOW THE LAND OR UNTIL THE MAXIMUM EXTINCTION DEPTH
C     IS REACHED. THE FIRST ENTRY POINT IS CALLED ONCE BYTHE MAIN
C     PROGRAM. THE SECOND ENTRY POINT IS CALLED ONCE FOR EACH TIMESTEP
C     ITERATION. THE THIRD ENTRY POINT IS CALLED ONCE FOR EACH TIMESTEP.
C
      INCLUDE 'DIMEN.INC'
      INCLUDE 'WUNITS.INC'
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
C
C     INTERNAL SCALARS AND ARRAYS.
C     INFILE      NAME OF INPUT FILE.
C     OUTFILE     OUTPUT FILE NAME.
C     OUTEXT      EXTENSION FOR OUTPUT FILE.
C     XCNAME      COMPARTMENT NAME.
C     ETNAME()    GROUNDWATER-EVAPOTRANSPRIATION RATE NAME ASSOCIATED
C                 WITH NODE.
C     ETNODE()    ET NODE NUMBER.
C     ETTAB()     GROUNDWATER-EVAPOTRANSPRIATION TABLE ASSOCIATED WITH
C                 NODE.
C     ETAREA()    AREA ASSOCIATED WITH GROUNDWATER-EVAPOTRANSPRIATION
C                 NODE.
C     LAND()      LAND-SURFACE ELEVATION AT NODE.
C     DELH0()     EXTINCTION DEPTH FOR VEGETATION.
C     ETMAX()     MAXIMUM EVAPOTRANSPIRATION RATE FOR VEGETATION
C     CRH()       RIGHT HANDED LEAKANCE COEFFICIENT.
C     CLH()       LEFT HANDED LEAKANCE COEFFICIENT.
C     QET()       FLUX RATE.
C     TTB()       ELAPSED TIME ASSOCIATED WITH POINT IN TABLE.
C     FTB()       TRANSIENT-STATE SCALING FACTOR ASSOCIATED WITH POINT
C                 IN TABLE.
C     NPT()       NUMBER OF POINTS IN SPECIFIED-FLUX TABLE.
C     IGET()      SWITCH FOR INTERPOLATION METHOD.
C     XTT()       TIME ENTRIES IN TABLE FOR PASSING TO TABLE FUNCTION.
C     XFT()       SPECIFIED-ET ENTRIES IN TABLE FOR PASSING TO TABLE
C                 FUNCTION.
C     H()         HEAD VALUE AT NODE.
C     IECHO1      SWITCH FOR ECHO OF GROUNDWATER-EVAPOTRANSPRIATION
C                 INPUTS.
C     IECHO2      SWITCH FOR OUTPUT OF TIME STEP
C                 GROUNDWATER-EVAPOTRANSPRIATION FLUXES.
C     NETNC       NUMBER OF GROUNDWATER-EVAPOTRANSPRIATION NODES IN
C                 COMPARTMENTS.
C     NTAB        NUMBER OF GROUNDWATER-EVAPOTRANSPRIATION TABLES
C                 ASSOCIATED WITH COMPARTMENTS.
C
C     INTERNAL ARRAYS.
      CHARACTER*80 INFILE,OUTFILE,OUTEXT
      CHARACTER*20 XCNAME
      CHARACTER*20 ETNAME(MAXETN)
      INTEGER*4 ETNODE(MAXETN),ETTAB(MAXETN)
      REAL*8 ETAREA(MAXETN),LAND(MAXETN),DELH0(MAXETN),ETMAX(MAXETN),
     +   CRH(MAXETN),CLH(MAXETN),QET(MAXETN)
      REAL*8 TTB(MAXTBE,MAXPTE),FTB(MAXTBE,MAXPTE)
      INTEGER*4 NPT(MAXTBE),IGET(MAXTBE)
      REAL*8 XTT(MAXPTE),XFT(MAXPTE)
      REAL*8 H(MAXNN)
C
C     COMMON BLOCKS.
C     FNAME()     NAMES OF INPUT FILES.
C     RNAME1      SIMULATION NAME FOR GROUNDWATER-FLOW FILES.
C     RNAME2      SIMULATION NAME FOR SOLUTE-TRANSPORT FILES.
C     NCELL       NUMBER OF COMPARTMENTS.
C     CNAME()     NAMES FOR COMPARTMENTS.
C     NNC()       NUMBER OF NODES IN A COMPARTMENT.
C     NEC()       NUMBER OF WEDGE ELEMENTS WITHIN A COMPARTMENT.
C     IN()        ELEMENT INCIDENCE FOR A WEDGE ELEMENT.
C     IN2()       ELEMENT INCIDENCES FOR A TETRAHEDRAL ELEMENT,
C                     WELL LINK, FAULT ELEMENT, OR FAULT LINK.
C     IEL()       PARENT WEDGE ELEMENT FOR TETRAHEDRAL ELEMENT.
C     X()         X-COORDINATE OF A NODE.
C     Y()         Y-COORDINATE OF A NODE.
C     Z()         Z-COORDINATE OF A NODE.
C     IA()        VECTOR INDICATING POSITIONS OF ROW ENDS IN
C                     VECTOR A() OR B().
C     JA()        VECTOR INDICATION COLUMN NUMBER ASSOCIATED WITH
C                     ENTRIES IN VECTORS A() OR B().
C     A()         VECTOR OF NON-ZERO COEFFICIENTS FOR CONDUCTANCE
C                     MATRIX.
C     B()         VECTOR OF NON-ZERO COEFFICIENTS FOR CAPACITANCE
C                     MATRIX.
C     F()         RIGHTHAND VECTOR FOR SYSTEM OF EQUATIONS.
C     QX()        DARCY GROUNDWATER VELOCITY IN THE X-DIRECTION.
C     QY()        DARCY GROUNDWATER VELOCITY IN THE Y-DIRECTION.
C     QZ()        DARCY GROUNDWATER VELOCITY IN THE Z-DIRECTION.
C     NE          GLOBAL NUMBER OF WEDGE ELEMENTS.
C     NE2         GLOBAL NUMBER OF TETRAHEDRAL ELEMENTS.
C     NE3         GLOBAL NUMBER OF TETRAHEDRAL ELEMENS AND WELL LINKS.
C     NE4         GLOBAL NUMBER OF TETHEDRALAL ELEMENTS, WELL LINKS,
C                     AND FAULT ELEMENTS.
C     NE5         GLOBAL NUMBER OF TETHEDRALAL ELEMENTS, WELL LINKS,
C                     FAULT ELEMENTS, AND FAULT LINKS.
C     NN          GLOBAL NUMBER OF MESH NODES.
C     NN3         GLOBAL NUMBER OF MESH NODES AND WELL NODES.
C     NN5         GLOBAL NUMBER OF MESH NODES, WELL NODES, AND
C                     FAULT NODES.
C     ISTEP       TIME STEP NUMBER WITHIN FLOW SIMULATION.
C     ISTEP1      STARTING STEP NUMBER FOR SIMULATION.
C     ISTEP2      ENDING STEP FOR SIMULATION.
C     DELT        LENGTH OF TIME STEP.
C     STIME       INITIAL ELAPSED TIME FOR SIMULATION.
C     IPRN        FLAG TO INDICATED PRINT OUTPUT FOR TIME STEP.
C     IPLT        FLAG TO INDICATE PLOT OUTPUT FOR TIME STEP.
C     NTSTEP      NUMBER OF TRANSPORT STEPS WITHIN FLOW STEP.
C     ITSTEP      TIME STEP NUMBER WITHIN TRANSPORT SIMULATION.
C     NIA         NUMBER OF ROWS REPRESENTED IN VECTORS A() OR B().
C     NJA         NUMBER OF NON-ZERO COEFFICIENTS IN VECTORS A() OR B().
C     ITRAN       SWITCH FOR TRANSPORT SIMULATION.
C     IRETN       INTERNAL READ FOR ET UNIT NUMBER.
C     IWETN       INTERNAL WRITE FOR ET UNIT NUMBER.
C
C     COMMON BLOCKS.
      CHARACTER*80 FNAME(30)
      CHARACTER*12 RNAME1,RNAME2
      CHARACTER*20 CNAME(MAXCEL)
      INTEGER*4 NNC(MAXCEL),NEC(MAXCEL)
      INTEGER*4 IN(MAXNE,6),IN2(MAXNE2,4),IEL(MAXNE2)
      REAL*8 X(MAXNN),Y(MAXNN),Z(MAXNN)
      INTEGER*4 IA(MAXNN+1),JA(MAXNON)
      REAL*8 A(MAXNON),B(MAXNON),F(MAXNN)
      REAL*8 QX(MAXNE2),QY(MAXNE2),QZ(MAXNE2)
      COMMON /NAMES/ FNAME,RNAME1,RNAME2
      COMMON /CELLS/ NCELL,CNAME,NNC,NEC
      COMMON /MESH/ X,Y,Z,IN,IN2,IEL,NE,NE2,NE3,NE4,NE5,NN,NN3,NN5
      COMMON /STEPS/ ISTEP,ISTEP1,ISTEP2,DELT,STIME,IPRN,IPLT,NTSTEP,
     +    ITSTEP
      COMMON /MATRIX/ A,B,IA,JA,NIA,NJA,F
      COMMON /SHARE/ QX,QY,QZ,ITRAN
      SAVE
C
C     OPEN INPUT AND OUTPUT FILES.
      WRITE(*,902)
      OUTEXT='.GET.OUT'
      INFILE=FNAME(9)
      CALL OPENFILES(RNAME1,INFILE,OUTFILE,OUTEXT,IRETN,IWETN)
      WRITE(IWETN,904)
      WRITE(IWETN,906) INFILE,OUTFILE
C
C     BEGIN LOOP OVER COMPARTMENTS.
      ILAST=0
      JLAST=0
      IETN=0
      NETN=0
      NTABLE=0
      READ(IRETN,910) IECHO2
      DO ICELL=1,NCELL
C
C         NUMBER OF ETNODES
          READ(IRETN,908) XCNAME
          READ(IRETN,910) NETNC,NTABC,IECHO1
          WRITE(IWETN,912) NCELL,ICELL,XCNAME,NETNC,NTABC,IECHO1,IECHO2
          NETN=NETN+NETNC
          NTABLE=NTABLE+NTABC
C
C         CHECK CELL NAME
          IF (XCNAME.NE.CNAME(ICELL)) THEN
              WRITE(*,914)
              WRITE(*, 916)
              STOP
          ENDIF
C
C         CHECK NUMBER OF ET NODES.
          IF (NETN.GT.MAXETN) THEN
              WRITE(*,918)
              WRITE(*,916)
              STOP
          ENDIF
C
C         CHECK NUMBER OF TABLES.
          IF (NTABLE.GT.MAXTBE) THEN
              WRITE(*,920)
              WRITE(*,916)
              STOP
          ENDIF
C
C         NODES, AREA, AND LAND SURFACE.
          IF (IECHO1.EQ.1) WRITE(IWETN,922)
          DO II=1,NETNC
              IETN=IETN+1
              READ(IRETN,924) ETNODE(IETN),ETAREA(IETN),LAND(IETN),
     +            DELH0(IETN),ETMAX(IETN),ETTAB(IETN),ETNAME(IETN)
C
C             CHECK NODE NUMBER.
              IF (ETNODE(IETN).GT.NNC(ICELL)) THEN
                  WRITE(*,926)
                  WRITE(*,916)
                  STOP
              ENDIF
              IF (ETTAB(IETN).NE.0) THEN
                  ITAB=ETTAB(IETN)+JLAST
              ELSE
                  ITAB=0
              ENDIF
              IF (IECHO1.EQ.1) WRITE(IWETN,928) ETNODE(IETN),
     +            ETNODE(IETN)+ILAST,ETAREA(IETN),LAND(IETN),
     +            DELH0(IETN),ETMAX(IETN),ETTAB(IETN),ITAB,ETNAME(IETN)
              ETNODE(IETN)=ETNODE(IETN)+ILAST
              ETTAB(IETN)=ITAB
          END DO
C
C         INPUT LOCAL FLUX TABLES.
          IF (NTABC.NE.0) THEN
              DO J=1,NTABC
C
C                 INPUT NUMBER OF POINTS IN LOCAL TABLE
                  ITAB=J+JLAST
                  READ(IRETN,930) NPT(ITAB),IGET(ITAB)
                  IF (IECHO1.NE.0) WRITE(IWETN,932) J,ITAB,NPT(ITAB),
     +                IGET(ITAB)
C
C                 CHECKING FOR VALID INTERPOLATION METHOD.
                  IF (IGET(ITAB).LT.1 .OR. IGET(ITAB).GT.2) THEN
                      WRITE(*,933)
                      WRITE(*,916)
                      STOP
                  ENDIF
C
C                 CHECK DIMENSIONING FOR NUMBER OF POINTS IN TABLE.
                  IF (NPT(ITAB).GT.MAXPTE) THEN
                      WRITE(*,934)
                      WRITE(*,916)
                      STOP
                  ENDIF
C
C                 INPUT TABLE ENTRIES.
                  DO IPT=1,NPT(ITAB)
                      READ(IRETN,936) TTB(ITAB,IPT),FTB(ITAB,IPT)
                      IF (IECHO1.NE.0) WRITE(IWETN,948) IPT,
     +                    TTB(ITAB,IPT),FTB(ITAB,IPT)
                  END DO
              END DO
          ENDIF
          ILAST=ILAST+NNC(ICELL)
          JLAST=NTABLE
      END DO
C
C     CLOSE FLIE AND EXIT WEVAP
      WRITE(*,940)
      WRITE(IWETN, 942)
      CLOSE (IRETN)

      RETURN

902   FORMAT(1X,'*START READING IN WEVAP')
904   FORMAT(1X,'**START READING IN WEVAP**')
906   FORMAT(/1X,'FILE NAMES:'/
     +    1X,'INPUT FILE..............................',A80/
     +    1X,'OUTPUT FILE.............................',A80)
908   FORMAT(A20)
910   FORMAT(4I10)
912   FORMAT(/1X,'EVAPOTRANSPIRATION PARAMETERS:'/
     +    1X,'NUMBER OF COMPARTMENTS..................',I10/
     +    1X,'COMPARTMENT NUMBER......................',I10/
     +    1X,'COMPARTMENT NAME........................',A20/
     +    1X,'NUMBER OF NODES.........................',I10/
     +    1X,'NUMBER OF TABLES........................',I10/
     +    1X,'SWITCH FOR ECHO OF INPUT DATA...........',I10/
     +    1X,'SWITCH FOR ECHO OF FLUXES...............',I10)
914   FORMAT(1X,'*ERROR* CELL NAME NOT IN LIST')
916   FORMAT(1X,'*ERROR* RUN TERMINATED FOR ERROR IN WEVAP')
918   FORMAT(1X,'*ERROR* NUMBER OF NODES EXCEEDS DIMENSION')
920   FORMAT(1X,'*ERROR* NUMBER OF TABLES EXCEEDS DIMENSION')
922   FORMAT(/1X,'NODAL PARAMETERS:'/
     +    1X,5X,'LOCAL',4X,'GLOBAL',45X,'LOCAL',4X,'GLOBAL'/
     +    1X,2(6X,'NODE'),6X,'AREA',6X,'LAND',5X,'DELH0',5X,'ETMAX',
     +    2(5X,'TABLE'),1X,'NAME')
924   FORMAT(I10,4F10.0,I10,1X,A20)
926   FORMAT(1X,'*ERROR* NODE NUMBER EXCEEDS NUMBER OF NODES')
928   FORMAT(1X,2I10,1PE10.2,2(0PF10.2),1PE10.2,2I10,1X,A20)
930   FORMAT(2I10)
932   FORMAT(/1X,'EVAPOTRANSPIRATION TABLE PARAMETERS:'/
     +    1X,'LOCAL TABLE NUMBER......................',I10/
     +    1X,'GLOBAL TABLE NUMBER.....................',I10/
     +    1X,'NUMBER OF ENTRIES IN TABLE..............',I10/
     +    1X,'SWITCH FOR INTERPOLATION METHOD.........',I10//
     +    1X,'TABLE ENTRIES:'/
     +    1X,5X,'POINT',6X,'TIME',6X,'RATE')
933   FORMAT(1X,'*ERROR* INTERPOLATION METHOD NOT VALID')
934   FORMAT(1X,'*ERROR* NUMBER OF POINTS EXCEEDS DIMENSION')
936   FORMAT(2F10.0)
938   FORMAT(1X,I10,2(1PE10.2))
940   FORMAT(1X,'*END   READING IN WEVAP')
942   FORMAT(/1X,'**END READING IN WEVAP**')
C
C     END OF DATA INPUT
C
C     ---------------
      ENTRY WEVAP2(H)
C     ---------------
C
C     INITIALIZE COEFFICIENTS.
      IF (NETN.EQ.0) RETURN
      DO I=1,NETN
          CRH(I)=0.0
          CLH(I)=0.0
      END DO
C
C     COMPUTE LEFT-HAND AND RIGHT-HAND COEFFICIENTS.
      DO IETN=1,NETN
          ITAB=ETTAB(IETN)
          INDE=ETNODE(IETN)
          IF (ITAB.NE.0) THEN
              DO IPT=1,NPT(ITAB)
                  XTT(IPT)=TTB(ITAB,IPT)
                  XFT(IPT)=FTB(ITAB,IPT)
              END DO
              XF=TABLE(NPT(ITAB),XTT,XFT,STIME-DELT,STIME,IGET(ITAB))
          ELSE
              XF=1.0
          ENDIF
          IF (H(INDE).LE.LAND(IETN)-DELH0(IETN)) THEN
C
C             DEPTH GREATER THAN EXTINCTION DEPTH.
              CRH(IETN)=0.0
              CLH(IETN)=0.0
          ELSE
C
C             DEPTH LESS THAN EXTINCTION DEPTH.
              COEF=ETAREA(IETN)*ETMAX(IETN)*XF/DELH0(IETN)
              CRH(IETN)=COEF*(LAND(IETN)-DELH0(IETN))
              CLH(IETN)=COEF
          ENDIF
C
C         PACK COEFFICIENTS INTO LEFTHAND MATRIX AND RIGHTHAND VECTOR.
          II=ETNODE(IETN)
          CALL PACK2(II,II,CLH(IETN),1)
          F(II)=F(II)+CRH(IETN)
      END DO
      RETURN
C
C     ---------------------
      ENTRY WEVAP3(H,SUMET)
C     ---------------------
C
C     COMPUTE ET FLUXES.
      SUMET=0.0
      IF (ISTEP.EQ.ISTEP1) WRITE(IWETN,944)
      DO IETN=1,NETN
          INDE=ETNODE(IETN)
          QET(IETN)=CRH(IETN)-CLH(IETN)*H(INDE)
          SUMET=SUMET+QET(IETN)
      END DO
C
C     DISPLAY RESULTS.
      IF (IPRN.EQ.1 .AND. IECHO2.EQ.1) THEN
          WRITE(IWETN,946) ISTEP,STIME,SUMET
          DO IETN=1,NETN
              INDE=ETNODE(IETN)
              WRITE(IWETN,948) ETNODE(IETN),LAND(IETN)-H(INDE),
     +            QET(IETN),ETNAME(IETN)
          END DO
      ENDIF
      IF (ISTEP.EQ.ISTEP2) THEN
          WRITE(IWETN,950)
          CLOSE (IWETN)
      ENDIF
      RETURN
C
C     FORMATS:
944   FORMAT(/1X,'**START SIMULATION IN WEVAP**')
946   FORMAT(/1X,'TIME-STEP SUMMARY FLUXES:'/
     +    1X,'TIME STEP...............................',I10/
     +    1X,'CUMULATIVE TIME.........................',1PE10.2/
     +    1X,'CUMULATIVE FLUX.........................',1PE10.2//
     +    1X,'NODE FLUXES:'/
     +    1X,6X,'NODE',5X,'DEPTH',3X,'OUTFLOW',1X,'NAME')
948   FORMAT(1X,I10,0PF10.2,1PE10.2,1X,A20)
950   FORMAT(/1X,'**END SIMULATION IN WEVAP**')
      END
