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