1      SUBROUTINE EVT5AL(ISUM,LENX,LCIEVT,LCEVTR,LCEXDP,LCSURF,
2     1                  NCOL,NROW,NEVTOP,IN,IOUT,IEVTCB,IFREFM)
3C
4C-----VERSION 0957 21FEB1996 EVT5AL
5C     ******************************************************************
6C     ALLOCATE ARRAY STORAGE FOR EVAPOTRANSPIRATION
7C     ******************************************************************
8C
9C        SPECIFICATIONS:
10C     ------------------------------------------------------------------
11C     ------------------------------------------------------------------
12C
13C1------IDENTIFY PACKAGE.
14      WRITE(IOUT,1)IN
15    1 FORMAT(1X,/1X,'EVT5 -- EVAPOTRANSPIRATION PACKAGE, VERSION 5,',
16     1     ' 9/1/93',' INPUT READ FROM UNIT',I3)
17C
18C2------READ ET OPTION (NEVTOP) AND UNIT OR FLAG FOR CELL-BY-CELL FLOW
19C2------TERMS (IEVTCB).
20      IF(IFREFM.EQ.0) THEN
21         READ(IN,'(2I10)') NEVTOP,IEVTCB
22      ELSE
23         READ(IN,*) NEVTOP,IEVTCB
24      END IF
25C
26C3------CHECK TO SEE THAT ET OPTION IS LEGAL.
27      IF(NEVTOP.GE.1.AND.NEVTOP.LE.2)GO TO 200
28C
29C3A-----IF ILLEGAL PRINT A MESSAGE & ABORT SIMULATION.
30      WRITE(IOUT,8)
31    8 FORMAT(1X,'ILLEGAL ET OPTION CODE. SIMULATION ABORTING')
32      STOP
33C
34C4------IF THE OPTION IS LEGAL THEN PRINT THE OPTION CODE.
35  200 IF(NEVTOP.EQ.1) WRITE(IOUT,201)
36  201 FORMAT(1X,'OPTION 1 -- EVAPOTRANSPIRATION FROM TOP LAYER')
37      IF(NEVTOP.EQ.2) WRITE(IOUT,202)
38  202 FORMAT(1X,'OPTION 2 -- EVAPOTRANSPIRATION FROM ONE SPECIFIED',
39     1      ' NODE IN EACH VERTICAL COLUMN')
40      IRK=ISUM
41C
42C5------IF CELL-BY-CELL FLOWS ARE TO BE SAVED, THEN PRINT UNIT NUMBER.
43      IF(IEVTCB.GT.0) WRITE(IOUT,203) IEVTCB
44  203 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE SAVED ON UNIT',I3)
45C
46C6------ALLOCATE SPACE FOR THE ARRAYS EVTR, EXDP AND SURF.
47      LCEVTR=ISUM
48      ISUM=ISUM+NCOL*NROW
49      LCEXDP=ISUM
50      ISUM=ISUM+NCOL*NROW
51      LCSURF=ISUM
52      ISUM=ISUM+NCOL*NROW
53C
54C7------IF OPTION 2 THEN ALLOCATE SPACE FOR THE INDICATOR ARRAY(IEVT)
55      LCIEVT=ISUM
56      IF(NEVTOP.NE.2)GO TO 300
57      ISUM=ISUM+NCOL*NROW
58C
59C8------CALCULATE & PRINT AMOUNT OF SPACE USED BY ET PACKAGE.
60  300 IRK=ISUM-IRK
61      WRITE(IOUT,4)IRK
62    4 FORMAT(1X,I10,' ELEMENTS OF X ARRAY ARE USED BY EVT')
63      ISUM1=ISUM-1
64      WRITE(IOUT,5)ISUM1,LENX
65    5 FORMAT(1X,I10,' ELEMENTS OF X ARRAY USED OUT OF ',I10)
66      IF(ISUM1.GT.LENX)WRITE(IOUT,6)
67    6 FORMAT(1X,'   ***X ARRAY MUST BE MADE LARGER***')
68C
69C9------RETURN.
70      RETURN
71      END
72      SUBROUTINE EVT5RP(NEVTOP,IEVT,EVTR,EXDP,SURF,DELR,DELC,
73     1                  NCOL,NROW,IN,IOUT,IFREFM)
74C
75C-----VERSION 1001 21FEB1996 EVT5RP
76C     ******************************************************************
77C     READ EVAPOTRANSPIRATION DATA
78C     ******************************************************************
79C
80C        SPECIFICATIONS:
81C     ------------------------------------------------------------------
82      CHARACTER*24 ANAME(4)
83      DIMENSION IEVT(NCOL,NROW),EVTR(NCOL,NROW),EXDP(NCOL,NROW),
84     1          SURF(NCOL,NROW),DELR(NCOL),DELC(NROW)
85C
86      DATA ANAME(1) /'          ET LAYER INDEX'/
87      DATA ANAME(2) /'              ET SURFACE'/
88      DATA ANAME(3) /' EVAPOTRANSPIRATION RATE'/
89      DATA ANAME(4) /'        EXTINCTION DEPTH'/
90C     ------------------------------------------------------------------
91C
92C1------READ FLAGS SHOWING WHETHER DATA IS TO BE REUSED.
93      IF(NEVTOP.EQ.2) THEN
94         IF(IFREFM.EQ.0) THEN
95            READ(IN,'(4I10)') INSURF,INEVTR,INEXDP,INIEVT
96         ELSE
97            READ(IN,*) INSURF,INEVTR,INEXDP,INIEVT
98         END IF
99      ELSE
100         IF(IFREFM.EQ.0) THEN
101            READ(IN,'(3I10)') INSURF,INEVTR,INEXDP
102         ELSE
103            READ(IN,*) INSURF,INEVTR,INEXDP
104         END IF
105      END IF
106C
107C2------TEST INSURF TO SEE WHERE SURFACE ELEVATION COMES FROM.
108      IF(INSURF.GE.0)GO TO 32
109C
110C2A------IF INSURF<0 THEN REUSE SURFACE ARRAY FROM LAST STRESS PERIOD
111      WRITE(IOUT,3)
112    3 FORMAT(1X,/1X,'REUSING SURF FROM LAST STRESS PERIOD')
113      GO TO 35
114C
115C3-------IF INSURF=>0 THEN CALL MODULE U2DREL TO READ SURFACE.
116   32 CALL U2DREL(SURF,ANAME(2),NROW,NCOL,0,IN,IOUT)
117C
118C4------TEST INEVTR TO SEE WHERE MAX ET RATE COMES FROM.
119   35 IF(INEVTR.GE.0)GO TO 37
120C
121C4A-----IF INEVTR<0 THEN REUSE MAX ET RATE.
122      WRITE(IOUT,4)
123    4 FORMAT(1X,/1X,'REUSING EVTR FROM LAST STRESS PERIOD')
124      GO TO 45
125C
126C5------IF INEVTR=>0 CALL MODULE U2DREL TO READ MAX ET RATE.
127   37 CALL U2DREL(EVTR,ANAME(3),NROW,NCOL,0,IN,IOUT)
128C
129C6------MULTIPLY MAX ET RATE BY CELL AREA TO GET VOLUMETRIC RATE
130      DO 40 IR=1,NROW
131      DO 40 IC=1,NCOL
132      EVTR(IC,IR)=EVTR(IC,IR)*DELR(IC)*DELC(IR)
133   40 CONTINUE
134C
135C7------TEST INEXDP TO SEE WHERE EXTINCTION DEPTH COMES FROM
136   45 IF(INEXDP.GE.0)GO TO 47
137C
138C7A------IF INEXDP<0 REUSE EXTINCTION DEPTH FROM LAST STRESS PERIOD
139      WRITE(IOUT,5)
140    5 FORMAT(1X,/1X,'REUSING EXDP FROM LAST STRESS PERIOD')
141      GO TO 48
142C
143C8-------IF INEXDP=>0 CALL MODULE U2DREL TO READ EXTINCTION DEPTH
144   47 CALL U2DREL(EXDP,ANAME(4),NROW,NCOL,0,IN,IOUT)
145C
146C9------IF OPTION(NEVTOP) IS 2 THEN WE NEED AN INDICATOR ARRAY.
147  48  IF(NEVTOP.NE.2)GO TO 50
148C
149C10------IF INIEVT<0 THEN REUSE LAYER INDICATOR ARRAY.
150      IF(INIEVT.GE.0)GO TO 49
151      WRITE(IOUT,2)
152    2 FORMAT(1X,/1X,'REUSING IEVT FROM LAST STRESS PERIOD')
153      GO TO 50
154C
155C11------IF INIEVT=>0 THEN CALL MODULE U2DINT TO READ INDICATOR ARRAY.
156   49 CALL U2DINT(IEVT,ANAME(1),NROW,NCOL,0,IN,IOUT)
157C
158C12-----RETURN
159   50 RETURN
160      END
161      SUBROUTINE EVT5FM(NEVTOP,IEVT,EVTR,EXDP,SURF,RHS,HCOF,
162     1                  IBOUND,HNEW,NCOL,NROW,NLAY)
163C
164C-----VERSION 1616 16JULY1992 EVT5FM
165C     ******************************************************************
166C        ADD EVAPOTRANSPIRATION TO RHS AND HCOF
167C     ******************************************************************
168C
169C        SPECIFICATIONS:
170C     ------------------------------------------------------------------
171      DOUBLE PRECISION HNEW,HH,SS,XX,DD
172      DIMENSION IEVT(NCOL,NROW),EVTR(NCOL,NROW),EXDP(NCOL,NROW),
173     1          SURF(NCOL,NROW),RHS(NCOL,NROW,NLAY),
174     2          HCOF(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY),
175     3          HNEW(NCOL,NROW,NLAY)
176C     ------------------------------------------------------------------
177C
178C1------PROCESS EACH HORIZONTAL CELL LOCATION
179      DO 10 IR=1,NROW
180      DO 10 IC=1,NCOL
181C
182C2------SET THE LAYER INDEX EQUAL TO 1      .
183      IL=1
184C
185C3------IF OPTION 2 IS SPECIFIED THEN GET LAYER INDEX FROM IEVT ARRAY
186      IF(NEVTOP.EQ.2)IL=IEVT(IC,IR)
187C
188C4------IF THE CELL IS EXTERNAL IGNORE IT.
189      IF(IBOUND(IC,IR,IL).LE.0)GO TO 10
190      C=EVTR(IC,IR)
191      S=SURF(IC,IR)
192      SS=S
193      HH=HNEW(IC,IR,IL)
194C
195C5------IF AQUIFER HEAD IS GREATER THAN OR EQUAL TO SURF, ET IS CONSTANT
196      IF(HH.LT.SS) GO TO 5
197C
198C5A-----SUBTRACT -EVTR FROM RHS
199      RHS(IC,IR,IL)=RHS(IC,IR,IL) + C
200      GO TO 10
201C
202C6------IF DEPTH TO WATER>=EXTINCTION DEPTH THEN ET IS 0
203    5 DD=SS-HH
204      X=EXDP(IC,IR)
205      XX=X
206      IF(DD.GE.XX)GO TO 10
207C
208C7------LINEAR RANGE. ADD ET TERMS TO BOTH RHS AND HCOF.
209      RHS(IC,IR,IL)=RHS(IC,IR,IL)+C-C*S/X
210      HCOF(IC,IR,IL)=HCOF(IC,IR,IL)-C/X
211   10 CONTINUE
212C
213C8------RETURN
214      RETURN
215      END
216      SUBROUTINE EVT5BD(NEVTOP,IEVT,EVTR,EXDP,SURF,IBOUND,HNEW,
217     1           NCOL,NROW,NLAY,DELT,VBVL,VBNM,MSUM,KSTP,KPER,
218     2           IEVTCB,ICBCFL,BUFF,IOUT,PERTIM,TOTIM)
219C-----VERSION 0829 18DEC1992 EVT5BD
220C     ******************************************************************
221C     CALCULATE VOLUMETRIC BUDGET FOR EVAPOTRANSPIRATION
222C     ******************************************************************
223C
224C        SPECIFICATIONS:
225C     ------------------------------------------------------------------
226      CHARACTER*16 VBNM(MSUM),TEXT
227      DOUBLE PRECISION HNEW,RATOUT,QQ,HH,SS,DD,XX,HHCOF,RRHS
228      DIMENSION IEVT(NCOL,NROW),EVTR(NCOL,NROW),EXDP(NCOL,NROW),
229     1          SURF(NCOL,NROW),IBOUND(NCOL,NROW,NLAY),
230     2          VBVL(4,MSUM),HNEW(NCOL,NROW,NLAY),BUFF(NCOL,NROW,NLAY)
231C
232      DATA TEXT /'              ET'/
233C     ------------------------------------------------------------------
234C
235C1------CLEAR THE RATE ACCUMULATOR.
236      ZERO=0.
237      RATOUT=ZERO
238C
239C2------SET CELL-BY-CELL BUDGET SAVE FLAG (IBD) AND CLEAR THE BUFFER.
240      IBD=0
241      IF(IEVTCB.GT.0) IBD=ICBCFL
242      DO 2 IL=1,NLAY
243      DO 2 IR=1,NROW
244      DO 2 IC=1,NCOL
245      BUFF(IC,IR,IL)=ZERO
246    2 CONTINUE
247C
248C3------PROCESS EACH HORIZONTAL CELL LOCATION.
249      DO 10 IR=1,NROW
250      DO 10 IC=1,NCOL
251C
252C4------SET THE LAYER INDEX EQUAL TO 1.
253      IL=1
254C
255C5------IF OPTION 2 IS SPECIFIED THEN GET LAYER INDEX FROM IEVT ARRAY.
256      IF(NEVTOP.EQ.2)IL=IEVT(IC,IR)
257C
258C6------IF CELL IS EXTERNAL THEN IGNORE IT.
259      IF(IBOUND(IC,IR,IL).LE.0)GO TO 10
260      C=EVTR(IC,IR)
261      S=SURF(IC,IR)
262      SS=S
263      HH=HNEW(IC,IR,IL)
264C
265C7------IF AQUIFER HEAD => SURF,SET Q=MAX ET RATE.
266      IF(HH.LT.SS) GO TO 7
267      QQ=-C
268      GO TO 9
269C
270C8------IF DEPTH=>EXTINCTION DEPTH, ET IS 0.
271    7 X=EXDP(IC,IR)
272      XX=X
273      DD=SS-HH
274      IF(DD.GE.XX)GO TO 10
275C
276C9------LINEAR RANGE.  Q=-EVTR*(HNEW-(SURF-EXDP))/EXDP, WHICH IS
277C9------FORMULATED AS Q= -HNEW*EVTR/EXDP + (EVTR*SURF/EXDP -EVTR).
278      HHCOF=-C/X
279      RRHS=(C*S/X)-C
280      QQ=HH*HHCOF+RRHS
281C
282C10-----ACCUMULATE TOTAL FLOW RATE.
283    9 Q=QQ
284      RATOUT=RATOUT-QQ
285C
286C11-----ADD Q TO BUFFER.
287      BUFF(IC,IR,IL)=Q
288   10 CONTINUE
289C
290C12-----IF CELL-BY-CELL FLOW TO BE SAVED, CALL APPROPRIATE UTILITY
291C12-----MODULE SAVE THEM.
292      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,IEVTCB,BUFF,NCOL,NROW,
293     1                          NLAY,IOUT)
294      IF(IBD.EQ.2) CALL UBDSV3(KSTP,KPER,TEXT,IEVTCB,BUFF,IEVT,NEVTOP,
295     1                   NCOL,NROW,NLAY,IOUT,DELT,PERTIM,TOTIM,IBOUND)
296C
297C13-----MOVE TOTAL ET RATE INTO VBVL FOR PRINTING BY BAS1OT.
298      ROUT=RATOUT
299      VBVL(3,MSUM)=ZERO
300      VBVL(4,MSUM)=ROUT
301C
302C14-----ADD ET(ET_RATE TIMES STEP LENGTH) TO VBVL.
303      VBVL(2,MSUM)=VBVL(2,MSUM)+ROUT*DELT
304C
305C15-----MOVE BUDGET TERM LABELS TO VBNM FOR PRINT BY MODULE BAS1OT.
306      VBNM(MSUM)=TEXT
307C
308C16-----INCREMENT BUDGET TERM COUNTER.
309      MSUM=MSUM+1
310C
311C17-----RETURN.
312      RETURN
313      END
314