SUBROUTINE VDF2LPF7FM(KITER,KSTP,KPER,IGRID) C ****************************************************************** C ADD LEAKAGE CORRECTION AND STORAGE TO HCOF AND RHS, AND CALCULATE C CONDUCTANCE AS REQUIRED. C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ USE GLOBAL, ONLY:NCOL,NROW,NLAY,IBOUND,BOTM,NBOTM,DELR,DELC, 1 LBOTM,CV,HNEW,RHS,HCOF,HOLD,ISSFLG,IOUT USE GWFBASMODULE,ONLY:DELT USE GWFLPFMODULE,ONLY:LAYTYP,SC1,SC2,NOVFC USE VDFMODULE, ONLY: DENSEREF,PS,ELEV,HSALT,MFNADVFD C ------------------------------------------------------------------ C C1------SET POINTERS TO DATA, GET STEADY-STATE FLAG FOR STRESS PERIOD, C1------DEFINE CONSTANT. CALL SGWF2LPF7PNT(IGRID) ISS=ISSFLG(KPER) ONE=1. C C2------FOR EACH LAYER: IF CONVERTIBLE, CALCULATE CONDUCTANCES. DO 100 K=1,NLAY KK=K IF(LAYTYP(K).NE.0) 1 CALL SGWF2LPF7HCOND(KK,KITER,KSTP,KPER) !,HSALT) 100 CONTINUE DO 101 K=1,NLAY KK=K IF(K.NE.NLAY) THEN IF(LAYTYP(K).NE.0 .OR. LAYTYP(K+1).NE.0) 1 CALL SGWF2LPF7VCOND(KK) END IF 101 CONTINUE C C3------IF THE STRESS PERIOD IS TRANSIENT, ADD STORAGE TO HCOF AND RHS IF(ISS.EQ.0) THEN TLED=ONE/DELT DO 200 K=1,NLAY C C4------SEE IF THIS LAYER IS CONVERTIBLE OR NON-CONVERTIBLE. IF(LAYTYP(K).EQ.0) THEN C5------NON-CONVERTIBLE LAYER, SO USE PRIMARY STORAGE DO 140 I=1,NROW DO 140 J=1,NCOL IF(IBOUND(J,I,K).LE.0) GO TO 140 RHO=SC1(J,I,K)*TLED C SEAWAT: CONSERVE MASS C HCOF(J,I,K)=HCOF(J,I,K)-RHO C RHS(J,I,K)=RHS(J,I,K)-RHO*HOLD(J,I,K) HCOF(J,I,K)=HCOF(J,I,K)-RHO*PS(J,I,K) RHS(J,I,K)=RHS(J,I,K)-RHO*HOLD(J,I,K)*PS(J,I,K) 140 CONTINUE ELSE C C5------A CONVERTIBLE LAYER, SO CHECK OLD AND NEW HEADS TO DETERMINE C5------WHEN TO USE PRIMARY AND SECONDARY STORAGE DO 180 I=1,NROW DO 180 J=1,NCOL C C5A-----IF THE CELL IS EXTERNAL THEN SKIP IT. IF(IBOUND(J,I,K).LE.0) GO TO 180 TP=BOTM(J,I,LBOTM(K)-1) RHO2=SC2(J,I,LAYTYP(K))*TLED RHO1=SC1(J,I,K)*TLED C C5B-----FIND STORAGE FACTOR AT START OF TIME STEP. SOLD=RHO2 C--SEAWAT: USE SALTHEADS FOR COMPARISON C IF(HOLD(J,I,K).GT.TP) SOLD=RHO1 HTMP=HOLD(J,I,K) IF(SALTHEAD(HTMP,PS(J,I,K),ELEV(J,I,K)).GT.TP) SOLD=RHO1 C C5C-----FIND STORAGE FACTOR AT END OF TIME STEP. C HTMP=HNEW(J,I,K) HTMP=HSALT(J,I,K) SNEW=RHO2 C--SEAWAT: USE SALTHEAD FOR COMPARISON C IF(HTMP.GT.TP) SNEW=RHO1 C IF(SALTHEAD(HTMP,PS(J,I,K),ELEV(J,I,K)).GT.TP) SNEW=RHO1 IF(HTMP.GT.TP) SNEW=RHO1 C C5D-----ADD STORAGE TERMS TO RHS AND HCOF. C--SEAWAT:CONSERVE MASS HCOF(J,I,K)=HCOF(J,I,K)-SNEW*PS(J,I,K) RHS(J,I,K)=RHS(J,I,K) - SOLD*(HOLD(J,I,K)-TP)*PS(J,I,K) & - SNEW*TP*PS(J,I,K) C 180 CONTINUE END IF C 200 CONTINUE END IF C C6------FOR EACH LAYER DETERMINE IF CORRECTION TERMS ARE NEEDED FOR C6------FLOW DOWN INTO PARTIALLY SATURATED LAYERS. DO 300 K=1,NLAY C C7------SEE IF CORRECTION IS NEEDED FOR LEAKAGE FROM ABOVE. IF(LAYTYP(K).NE.0 .AND. K.NE.1) THEN C C7A-----FOR EACH CELL MAKE THE CORRECTION IF NEEDED. DO 220 I=1,NROW DO 220 J=1,NCOL C C7B-----IF THE CELL IS EXTERNAL(IBOUND<=0) THEN SKIP IT. IF(IBOUND(J,I,K).LE.0) GO TO 220 C C7C-----IF HEAD IS ABOVE TOP THEN CORRECTION NOT NEEDED TOP=BOTM(J,I,LBOTM(K)-1) IF(HSALT(J,I,K).GE.TOP) GO TO 220 C C7D-----WITH HEAD BELOW TOP ADD CORRECTION TERMS TO RHS. C--SEAWAT: USE VD FORM, CONSERVE MASS C RHS(J,I,K)=RHS(J,I,K) + PS(J,I,K-1)*CV(J,I,K-1)* C & (ELEV(J,I,K)-HNFE+PS(J,I,K-1)/DENSEREF* C & (TOP-ELEV(J,I,K))) C C--SEAWAT: THIS SECTION REWRITTEN 9/24/07 USING CLEARER CODING C--FIRST CALCULATE Q_MCALC, WHICH IS THE MASS FLUX WITH THE STANDARD DARCY EQUATION DIS1=ELEV(J,I,K-1)-BOTM(J,I,K-1) DIS2=BOTM(J,I,K-1)-ELEV(J,I,K) AVGDENS=(DIS1*PS(J,I,K-1)+DIS2*PS(J,I,K))/(DIS1+DIS2) Q_MCALC=CV(J,I,K-1)*(HNEW(J,I,K-1)-HNEW(J,I,K)+ + (AVGDENS-DENSEREF)/DENSEREF*(ELEV(J,I,K-1)-ELEV(J,I,K))) IF(MFNADVFD.EQ.2) THEN Q_MCALC=Q_MCALC*AVGDENS ELSE IF (Q_MCALC.GT.0.) THEN Q_MCALC=Q_MCALC*PS(J,I,K-1) ELSE Q_MCALC=Q_MCALC*PS(J,I,K) ENDIF ENDIF C C--CALCULATE Q_MACTUAL, WHICH IS THE TRUE MASS FLUX Q_MACTUAL=PS(J,I,K-1)/DENSEREF*CV(J,I,K-1)* + (SALTHEAD(HNEW(J,I,K-1),PS(J,I,K-1),ELEV(J,I,K-1)) + -BOTM(J,I,LBOTM(K)-1)) Q_MACTUAL=Q_MACTUAL*PS(J,I,K-1) C C--CALCULATE Q_MCORRECT, WHICH IS THE CORRECTION TERM THAT IS ADDED TO THE RHS Q_MCORRECT=Q_MCALC-Q_MACTUAL RHS(J,I,K)=RHS(J,I,K)+Q_MCORRECT 220 CONTINUE END IF C C8------SEE IF THIS LAYER MAY NEED CORRECTION FOR LEAKAGE TO BELOW. IF(K.EQ.NLAY) GO TO 300 IF(LAYTYP(K+1).NE.0) THEN C C8A-----FOR EACH CELL MAKE THE CORRECTION IF NEEDED. DO 280 I=1,NROW DO 280 J=1,NCOL C C8B-----IF CELL IS EXTERNAL (IBOUND<=0) THEN SKIP IT. IF(IBOUND(J,I,K).LE.0) GO TO 280 C C8C-----IF HEAD IN THE LOWER CELL IS LESS THAN TOP ADD CORRECTION C8C-----TERM TO RHS. HTMP=HNEW(J,I,K+1) TOP=BOTM(J,I,LBOTM(K+1)-1) IF(SALTHEAD(HTMP,PS(J,I,K+1),ELEV(J,I,K+1)).GE. & BOTM(J,I,LBOTM(K+1)-1)) GOTO 280 C C--FIRST CALCULATE Q_MCALC, WHICH IS THE MASS FLUX WITH THE STANDARD DARCY EQUATION DIS1=BOTM(J,I,K)-ELEV(J,I,K+1) DIS2=ELEV(J,I,K)-BOTM(J,I,K) AVGDENS=(DIS1*PS(J,I,K+1)+DIS2*PS(J,I,K))/(DIS1+DIS2) Q_MCALC=CV(J,I,K)*(HNEW(J,I,K+1)-HNEW(J,I,K)+ + (AVGDENS-DENSEREF)/DENSEREF*(ELEV(J,I,K+1)-ELEV(J,I,K))) IF(MFNADVFD.EQ.2) THEN Q_MCALC=Q_MCALC*AVGDENS ELSE IF (Q_MCALC.GT.0.) THEN Q_MCALC=Q_MCALC*PS(J,I,K+1) ELSE Q_MCALC=Q_MCALC*PS(J,I,K) ENDIF ENDIF C C--CALCULATE Q_MACTUAL, WHICH IS THE TRUE OR ACTUAL MASS FLUX Q_MACTUAL=PS(J,I,K)/DENSEREF*CV(J,I,K)* + (BOTM(J,I,K)-SALTHEAD(HNEW(J,I,K),PS(J,I,K),ELEV(J,I,K))) Q_MACTUAL=Q_MACTUAL*PS(J,I,K) C C--CALCULATE Q_MCORRECT, WHICH IS THE CORRECTION TERM THAT IS ADDED TO THE RHS Q_MCORRECT=Q_MCALC-Q_MACTUAL RHS(J,I,K)=RHS(J,I,K)+Q_MCORRECT 280 CONTINUE END IF C 300 CONTINUE C C9------RETURN RETURN END SUBROUTINE VDF2LPF7BDS(KSTP,KPER,IGRID) C ****************************************************************** C COMPUTE STORAGE BUDGET FLOW TERM FOR LPF. C--SEAWAT: MODIFIED TO CONSERVE MASS C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ USE GLOBAL, ONLY:NCOL,NROW,NLAY,ISSFLG,IBOUND,HNEW,HOLD, 1 BUFF,BOTM,LBOTM,IOUT USE GWFBASMODULE,ONLY:MSUM,ICBCFL,VBVL,VBNM,DELT,PERTIM,TOTIM USE GWFLPFMODULE,ONLY:ILPFCB,LAYTYP,SC1,SC2 USE VDFMODULE, ONLY: PS,ELEV,HSALT CHARACTER*16 TEXT DOUBLE PRECISION STOIN,STOUT,SSTRG C--SEAWAT: ADD HTMP,HSING DOUBLE PRECISION HTMP,HSING C DATA TEXT /' STORAGE'/ C ------------------------------------------------------------------ C CALL SGWF2LPF7PNT(IGRID) C C1------INITIALIZE BUDGET ACCUMULATORS AND 1/DELT. ISS=ISSFLG(KPER) ZERO=0. STOIN=ZERO STOUT=ZERO C2------IF STEADY STATE, STORAGE TERM IS ZERO IF(ISS.NE.0) GOTO 400 ONE=1. TLED=ONE/DELT C C3------IF CELL-BY-CELL FLOWS WILL BE SAVED, SET FLAG IBD. IBD=0 IF(ILPFCB.GT.0) IBD=ICBCFL C C4------CLEAR BUFFER. DO 210 K=1,NLAY DO 210 I=1,NROW DO 210 J=1,NCOL BUFF(J,I,K)=ZERO 210 CONTINUE C C5------LOOP THROUGH EVERY CELL IN THE GRID. KT=0 DO 300 K=1,NLAY LC=LAYTYP(K) IF(LC.NE.0) KT=KT+1 DO 300 I=1,NROW DO 300 J=1,NCOL C C6------SKIP NO-FLOW AND CONSTANT-HEAD CELLS. IF(IBOUND(J,I,K).LE.0) GO TO 300 HSING=HNEW(J,I,K) C C7-----CHECK LAYER TYPE TO SEE IF ONE STORAGE CAPACITY OR TWO. IF(LC.EQ.0) GO TO 285 C C7A----TWO STORAGE CAPACITIES. TP=BOTM(J,I,LBOTM(K)-1) RHO2=SC2(J,I,KT)*TLED RHO1=SC1(J,I,K)*TLED SOLD=RHO2 C--SEAWAT: USE SALTHEAD FOR COMPARISON IF(SALTHEAD(HOLD(J,I,K),PS(J,I,K),ELEV(J,I,K)).GT.TP) SOLD=RHO1 SNEW=RHO2 C--SEAWAT: USE SALTHEAD FOR COMPARISON IF(SALTHEAD(HSING,PS(J,I,K),ELEV(J,I,K)).GT.TP) SNEW=RHO1 C--SEAWAT: CONSERVE MASS STRG=PS(J,I,K)*SOLD*(HOLD(J,I,K)-TP) + PS(J,I,K)*SNEW*TP + - PS(J,I,K)*SNEW*HSING GO TO 288 C C7B----ONE STORAGE CAPACITY. 285 RHO=SC1(J,I,K)*TLED C--SEAWAT: CONSERVE MASS STRG=PS(J,I,K)*(RHO*HOLD(J,I,K) - RHO*HSING) C C8-----STORE CELL-BY-CELL FLOW IN BUFFER AND ADD TO ACCUMULATORS. C--SEAWAT: CONVERT BACK TO VOLUME FOR CELL BY CELL FLOW C 288 BUFF(J,I,K)=STRG 288 BUFF(J,I,K)=STRG/PS(J,I,K) SSTRG=STRG IF(STRG.LT.ZERO) THEN STOUT=STOUT-SSTRG ELSE STOIN=STOIN+SSTRG END IF C 300 CONTINUE C C9-----IF IBD FLAG IS SET RECORD THE CONTENTS OF THE BUFFER. IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT, 1 ILPFCB,BUFF,NCOL,NROW,NLAY,IOUT) IF(IBD.EQ.2) CALL UBDSV1(KSTP,KPER,TEXT,ILPFCB, 1 BUFF,NCOL,NROW,NLAY,IOUT,DELT,PERTIM,TOTIM,IBOUND) C C10-----ADD TOTAL RATES AND VOLUMES TO VBVL & PUT TITLE IN VBNM. 400 CONTINUE SIN=STOIN SOUT=STOUT VBVL(1,MSUM)=VBVL(1,MSUM)+SIN*DELT VBVL(2,MSUM)=VBVL(2,MSUM)+SOUT*DELT VBVL(3,MSUM)=SIN VBVL(4,MSUM)=SOUT VBNM(MSUM)=TEXT MSUM=MSUM+1 C C11----RETURN. RETURN END SUBROUTINE VDF2LPF7BDCH(KSTP,KPER,IGRID) C ****************************************************************** C COMPUTE FLOW FROM CONSTANT-HEAD CELLS C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ USE GLOBAL, ONLY:NCOL,NROW,NLAY,IBOUND,HNEW,BUFF,CR,CC,CV, 1 BOTM,LBOTM,IOUT,IUNIT,DELR,DELC USE M_MF2005_IU, ONLY : IUANI USE GWFBASMODULE,ONLY:MSUM,VBVL,VBNM,DELT,PERTIM,TOTIM,ICBCFL, 1 ICHFLG USE GWFLPFMODULE,ONLY:ILPFCB,LAYTYP,NOVFC USE VDFMODULE, ONLY: MFNADVFD,IWTABLE,DENSEREF,PS,ELEV,HSALT CHARACTER*16 TEXT DOUBLE PRECISION HD,CHIN,CHOUT, & XX1,XX2,XX3,XX4,XX5,XX6, & H1,HS1,H2,HS2 C DATA TEXT /' CONSTANT HEAD'/ C ------------------------------------------------------------------ CALL SGWF2LPF7PNT(IGRID) C C1------SET IBD TO INDICATE IF CELL-BY-CELL BUDGET VALUES WILL BE SAVED. IBD=0 IF(ILPFCB.LT.0 .AND. ICBCFL.NE.0) IBD=-1 IF(ILPFCB.GT.0) IBD=ICBCFL C C2------CLEAR BUDGET ACCUMULATORS. ZERO=0. CHIN=ZERO CHOUT=ZERO IBDLBL=0 C C3------CLEAR BUFFER. DO 5 K=1,NLAY DO 5 I=1,NROW DO 5 J=1,NCOL BUFF(J,I,K)=ZERO 5 CONTINUE !## compute addition of diagonal fluxes from constant heads if(iunit(iuani).gt.0)then call gwf2ani8bd_chd(IGRID) DO K=1,NLAY DO I=1,NROW DO J=1,NCOL IF(buff(j,i,k).LT.ZERO) THEN CHOUT=CHOUT-denseref*buff(j,i,k) ELSE CHIN=CHIN+denseref*buff(j,i,k) END IF ENDDO ENDDO ENDDO ENDIF C C3A-----IF SAVING CELL-BY-CELL FLOW IN A LIST, COUNT CONSTANT-HEAD C3A-----CELLS AND WRITE HEADER RECORDS. IF(IBD.EQ.2) THEN NCH=0 DO 7 K=1,NLAY DO 7 I=1,NROW DO 7 J=1,NCOL IF(IBOUND(J,I,K).LT.0) NCH=NCH+1 7 CONTINUE CALL UBDSV2(KSTP,KPER,TEXT,ILPFCB,NCOL,NROW,NLAY, 1 NCH,IOUT,DELT,PERTIM,TOTIM,IBOUND) END IF C C4------LOOP THROUGH EACH CELL AND CALCULATE FLOW INTO MODEL FROM EACH C4------CONSTANT-HEAD CELL. DO 200 K=1,NLAY DO 200 I=1,NROW DO 200 J=1,NCOL C C5------IF CELL IS NOT CONSTANT HEAD SKIP IT & GO ON TO NEXT CELL. IF (IBOUND(J,I,K).GE.0)GO TO 200 C C6------CLEAR VALUES FOR FLOW RATE THROUGH EACH FACE OF CELL. X1=ZERO X2=ZERO X3=ZERO X4=ZERO X5=ZERO X6=ZERO CHCH1=ZERO CHCH2=ZERO CHCH3=ZERO CHCH4=ZERO CHCH5=ZERO CHCH6=ZERO C C--SEAWAT: INTRODUCE THE TERM FOR STORING VOLUMETRIC FLOW FLOW=ZERO C C7------CALCULATE FLOW THROUGH THE LEFT FACE. C7------COMMENTS A-C APPEAR ONLY IN THE SECTION HEADED BY COMMENT 7, C7------BUT THEY APPLY IN A SIMILAR MANNER TO SECTIONS 8-12. C C7A-----IF THERE IS NO FLOW TO CALCULATE THROUGH THIS FACE, THEN GO ON C7A-----TO NEXT FACE. NO FLOW OCCURS AT THE EDGE OF THE GRID, TO AN C7A-----ADJACENT NO-FLOW CELL, OR TO AN ADJACENT CONSTANT-HEAD CELL. IF(J.EQ.1) GO TO 30 IF(IBOUND(J-1,I,K).EQ.0) GO TO 30 IF(IBOUND(J-1,I,K).LT.0 .AND. ICHFLG.EQ.0) GO TO 30 C C7B-----CALCULATE FLOW THROUGH THIS FACE INTO THE ADJACENT CELL. H1=HNEW(J,I,K) HS1=HSALT(J,I,K) H2=HNEW(J-1,I,K) HS2=HSALT(J-1,I,K) Z1=ELEV(J,I,K) Z2=ELEV(J-1,I,K) PS1=PS(J,I,K) PS2=PS(J-1,I,K) TOP1=BOTM(J,I,K-1) TOP2=BOTM(J-1,I,K-1) BOT1=BOTM(J,I,K) BOT2=BOTM(J-1,I,K) C--SEAWAT: CALL VDWTABLE IF WATER TABLE CORRECTIONS ARE ACTIVE IF(IWTABLE.EQ.1.AND.LAYTYP(K).NE.0) & CALL VDWTABLE(H1,HS1,PS1,Z1,TOP1,BOT1, & H2,HS2,PS2,Z2,TOP2,BOT2) HDIFF=H2-H1 DIS1=DELR(J-1)/2 DIS2=DELR(J)/2 AVGDENS=(DIS1*PS(J-1,I,K)+DIS2*PS(J,I,K))/(DIS1+DIS2) D1=CR(J-1,I,K)*(AVGDENS-DENSEREF)/DENSEREF*(Z2-Z1) FLOWDIR=CR(J-1,I,K)*HDIFF+D1 FLOW=FLOW-FLOWDIR IF(MFNADVFD.EQ.2) THEN CHCH1=-FLOWDIR*AVGDENS ELSE IF(FLOWDIR.GE.0.) THEN CHCH1=-FLOWDIR*PS(J-1,I,K) ELSE CHCH1=-FLOWDIR*PS(J,I,K) ENDIF ENDIF IF(IBOUND(J-1,I,K).LT.0) GO TO 30 X1=CHCH1 XX1=X1 C C7C-----ACCUMULATE POSITIVE AND NEGATIVE FLOW. IF (X1) 10,30,20 10 CHOUT=CHOUT-XX1 GO TO 30 20 CHIN=CHIN+XX1 C C8------CALCULATE FLOW THROUGH THE RIGHT FACE. 30 IF(J.EQ.NCOL) GO TO 60 IF(IBOUND(J+1,I,K).EQ.0) GO TO 60 IF(IBOUND(J+1,I,K).LT.0 .AND. ICHFLG.EQ.0) GO TO 60 H1=HNEW(J,I,K) HS1=HSALT(J,I,K) H2=HNEW(J+1,I,K) HS2=HSALT(J+1,I,K) Z1=ELEV(J,I,K) Z2=ELEV(J+1,I,K) PS1=PS(J,I,K) PS2=PS(J+1,I,K) TOP1=BOTM(J,I,K-1) TOP2=BOTM(J+1,I,K-1) BOT1=BOTM(J,I,K) BOT2=BOTM(J+1,I,K) C--SEAWAT: CALL VDWTABLE IF WATER TABLE CORRECTIONS ARE ACTIVE IF(IWTABLE.EQ.1.AND.LAYTYP(K).NE.0) & CALL VDWTABLE(H1,HS1,PS1,Z1,TOP1,BOT1, & H2,HS2,PS2,Z2,TOP2,BOT2) HDIFF=H2-H1 DIS1=DELR(J+1)/2 DIS2=DELR(J)/2 AVGDENS=(DIS1*PS(J+1,I,K)+DIS2*PS(J,I,K))/(DIS1+DIS2) D2=CR(J,I,K)*(AVGDENS-DENSEREF)/DENSEREF*(Z2-Z1) FLOWDIR=CR(J,I,K)*HDIFF+D2 FLOW=FLOW-FLOWDIR IF(MFNADVFD.EQ.2) THEN CHCH2=-FLOWDIR*AVGDENS ELSE IF(FLOWDIR.GE.0.) THEN CHCH2=-FLOWDIR*PS(J+1,I,K) ELSE CHCH2=-FLOWDIR*PS(J,I,K) ENDIF ENDIF IF(IBOUND(J+1,I,K).LT.0) GO TO 60 X2=CHCH2 XX2=X2 IF(X2)40,60,50 40 CHOUT=CHOUT-XX2 GO TO 60 50 CHIN=CHIN+XX2 C C9------CALCULATE FLOW THROUGH THE BACK FACE. 60 IF(I.EQ.1) GO TO 90 IF (IBOUND(J,I-1,K).EQ.0) GO TO 90 IF (IBOUND(J,I-1,K).LT.0 .AND. ICHFLG.EQ.0) GO TO 90 H1=HNEW(J,I,K) HS1=HSALT(J,I,K) H2=HNEW(J,I-1,K) HS2=HSALT(J,I-1,K) Z1=ELEV(J,I,K) Z2=ELEV(J,I-1,K) PS1=PS(J,I,K) PS2=PS(J,I-1,K) TOP1=BOTM(J,I,K-1) TOP2=BOTM(J,I-1,K-1) BOT1=BOTM(J,I,K) BOT2=BOTM(J,I-1,K) C--SEAWAT: CALL VDWTABLE IF WATER TABLE CORRECTIONS ARE ACTIVE IF(IWTABLE.EQ.1.AND.LAYTYP(K).NE.0) & CALL VDWTABLE(H1,HS1,PS1,Z1,TOP1,BOT1, & H2,HS2,PS2,Z2,TOP2,BOT2) HDIFF=H2-H1 DIS1=DELC(I-1)/2 DIS2=DELC(I)/2 AVGDENS=(DIS1*PS(J,I-1,K)+DIS2*PS(J,I,K))/(DIS1+DIS2) D3=CC(J,I-1,K)*(AVGDENS-DENSEREF)/DENSEREF*(Z2-Z1) FLOWDIR=CC(J,I-1,K)*HDIFF+D3 FLOW=FLOW-FLOWDIR IF(MFNADVFD.EQ.2) THEN CHCH3=-FLOWDIR*AVGDENS ELSE IF(FLOWDIR.GE.0.) THEN CHCH3=-FLOWDIR*PS(J,I-1,K) ELSE CHCH3=-FLOWDIR*PS(J,I,K) ENDIF ENDIF IF(IBOUND(J,I-1,K).LT.0) GO TO 90 X3=CHCH3 XX3=X3 IF(X3) 70,90,80 70 CHOUT=CHOUT-XX3 GO TO 90 80 CHIN=CHIN+XX3 C C10-----CALCULATE FLOW THROUGH THE FRONT FACE. 90 IF(I.EQ.NROW) GO TO 120 IF(IBOUND(J,I+1,K).EQ.0) GO TO 120 IF(IBOUND(J,I+1,K).LT.0 .AND. ICHFLG.EQ.0) GO TO 120 H1=HNEW(J,I,K) HS1=HSALT(J,I,K) H2=HNEW(J,I+1,K) HS2=HSALT(J,I+1,K) Z1=ELEV(J,I,K) Z2=ELEV(J,I+1,K) PS1=PS(J,I,K) PS2=PS(J,I+1,K) TOP1=BOTM(J,I,K-1) TOP2=BOTM(J,I+1,K-1) BOT1=BOTM(J,I,K) BOT2=BOTM(J,I+1,K) IF(IWTABLE.EQ.1.AND.LAYTYP(K).NE.0) & CALL VDWTABLE(H1,HS1,PS1,Z1,TOP1,BOT1, & H2,HS2,PS2,Z2,TOP2,BOT2) HDIFF=H2-H1 DIS1=DELC(I+1)/2 DIS2=DELC(I)/2 AVGDENS=(DIS1*PS(J,I+1,K)+DIS2*PS(J,I,K))/(DIS1+DIS2) D4=CC(J,I,K)*(AVGDENS-DENSEREF)/DENSEREF*(Z2-Z1) FLOWDIR=CC(J,I,K)*HDIFF+D4 FLOW=FLOW-FLOWDIR IF(MFNADVFD.EQ.2) THEN CHCH4=-FLOWDIR*AVGDENS ELSE IF(FLOWDIR.GE.0.) THEN CHCH4=-FLOWDIR*PS(J,I+1,K) ELSE CHCH4=-FLOWDIR*PS(J,I,K) ENDIF ENDIF IF(IBOUND(J,I+1,K).LT.0) GO TO 120 X4=CHCH4 XX4=X4 IF (X4) 100,120,110 100 CHOUT=CHOUT-XX4 GO TO 120 110 CHIN=CHIN+XX4 C C11-----CALCULATE FLOW THROUGH THE UPPER FACE. 120 IF(K.EQ.1) GO TO 150 IF (IBOUND(J,I,K-1).EQ.0) GO TO 150 IF (IBOUND(J,I,K-1).LT.0 .AND. ICHFLG.EQ.0) GO TO 150 HD=HNEW(J,I,K) IF(LAYTYP(K).EQ.0) GO TO 122 TMP=HD TOP=BOTM(J,I,LBOTM(K)-1) IF(TMP.LT.TOP) HD=TOP 122 HDIFF=HNEW(J,I,K-1)-HD DIS1=ELEV(J,I,K-1)-BOTM(J,I,K-1) DIS2=BOTM(J,I,K-1)-ELEV(J,I,K) AVGDENS=(DIS1*PS(J,I,K-1)+DIS2*PS(J,I,K))/(DIS1+DIS2) D5=CV(J,I,K-1)*(AVGDENS-DENSEREF)/DENSEREF*(ELEV(J,I,K-1)- & ELEV(J,I,K)) FLOWDIR=CV(J,I,K-1)*HDIFF+D5 C--SEAWAT: CHECK AND CORRECT FOR DEWATERED CASE IF(LAYTYP(K).GT.0) THEN HS2=SALTHEAD(HNEW(J,I,K),PS(J,I,K),ELEV(J,I,K)) BOT1=BOTM(J,I,LBOTM(K)-1) IF(HS2.LT.BOT1) THEN HS1=SALTHEAD(HNEW(J,I,K-1),PS(J,I,K-1),ELEV(J,I,K-1)) FLOWDIR=PS(J,I,K-1)/DENSEREF*CV(J,I,K-1)*(HS1-BOT1) ENDIF ENDIF C--SEAWAT: END DEWATERED CORRECTION FLOW=FLOW-FLOWDIR IF(MFNADVFD.EQ.2) THEN CHCH5=-FLOWDIR*AVGDENS ELSE IF(FLOWDIR.GE.0.) THEN CHCH5=-FLOWDIR*PS(J,I,K-1) ELSE CHCH5=-FLOWDIR*PS(J,I,K) ENDIF ENDIF IF(IBOUND(J,I,K-1).LT.0) GO TO 150 X5=CHCH5 XX5=X5 IF(X5) 130,150,140 130 CHOUT=CHOUT-XX5 GO TO 150 140 CHIN=CHIN+XX5 C C12-----CALCULATE FLOW THROUGH THE LOWER FACE. 150 IF(K.EQ.NLAY) GO TO 180 IF(IBOUND(J,I,K+1).EQ.0) GO TO 180 IF(IBOUND(J,I,K+1).LT.0 .AND. ICHFLG.EQ.0) GO TO 180 HD=HNEW(J,I,K+1) IF(LAYTYP(K+1).EQ.0) GO TO 152 TMP=HD TOP=BOTM(J,I,LBOTM(K+1)-1) IF(TMP.LT.TOP) HD=TOP 152 HDIFF=HD-HNEW(J,I,K) DIS1=BOTM(J,I,K)-ELEV(J,I,K+1) DIS2=ELEV(J,I,K)-BOTM(J,I,K) AVGDENS=(DIS1*PS(J,I,K+1)+DIS2*PS(J,I,K))/(DIS1+DIS2) D6=CV(J,I,K)*(AVGDENS-DENSEREF)/DENSEREF*(ELEV(J,I,K+1)- & ELEV(J,I,K)) FLOWDIR=CV(J,I,K)*HDIFF+D6 C--SEAWAT: CHECK AND CORRECT FOR DEWATERED CASE IF(LAYTYP(K+1).GT.0) THEN HS2=SALTHEAD(HNEW(J,I,K+1),PS(J,I,K+1),ELEV(J,I,K+1)) BOT1=BOTM(J,I,LBOTM(K+1)-1) IF(HS2.LT.BOT1) THEN HS1=SALTHEAD(HNEW(J,I,K),PS(J,I,K),ELEV(J,I,K)) FLOWDIR=PS(J,I,K)/DENSEREF*CV(J,I,K)*(BOT1-HS1) ENDIF ENDIF C--SEAWAT: END DEWATERED CORRECTION FLOW=FLOW-FLOWDIR IF(MFNADVFD.EQ.2) THEN CHCH6=-FLOWDIR*AVGDENS ELSE IF(FLOWDIR.GE.0.) THEN CHCH6=-FLOWDIR*PS(J,I,K+1) ELSE CHCH6=-FLOWDIR*PS(J,I,K) ENDIF ENDIF IF(IBOUND(J,I,K+1).LT.0) GO TO 180 X6=CHCH6 XX6=X6 IF(X6) 160,180,170 160 CHOUT=CHOUT-XX6 GO TO 180 170 CHIN=CHIN+XX6 C C13-----SUM THE FLOWS THROUGH SIX FACES OF CONSTANT HEAD CELL, AND C13-----STORE SUM IN BUFFER. 180 RATE=CHCH1+CHCH2+CHCH3+CHCH4+CHCH5+CHCH6 if(iunit(iuani).gt.0)then BUFF(J,I,K)=denseref*buff(j,i,k)+RATE else BUFF(J,I,K)=RATE endif C C14-----PRINT THE FLOW FOR THE CELL IF REQUESTED. IF(IBD.LT.0) THEN IF(IBDLBL.EQ.0) WRITE(IOUT,899) TEXT,KPER,KSTP 899 FORMAT(1X,/1X,A,' PERIOD',I4,' STEP',I3) C--SEAWAT: WRITE FLOW RATHER THAN RATE WRITE(IOUT,900) K,I,J,FLOW 900 FORMAT(1X,'LAYER',I3,' ROW',I5,' COL',I5, 1 ' RATE',1PG15.6) IBDLBL=1 END IF C C--SEAWAT: SAVE VOLUMETRIX FLUX (FLOW) RATHER THAN MASS FLUX (RATE) RATE=FLOW BUFF(J,I,K)=RATE C15-----IF SAVING CELL-BY-CELL FLOW IN LIST, WRITE FLOW FOR CELL. IF(IBD.EQ.2) CALL UBDSVA(ILPFCB,NCOL,NROW,J,I,K,RATE,IBOUND,NLAY) 200 CONTINUE C C16-----IF SAVING CELL-BY-CELL FLOW IN 3-D ARRAY, WRITE THE ARRAY. IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT, 1 ILPFCB,BUFF,NCOL,NROW,NLAY,IOUT) C C17-----SAVE TOTAL CONSTANT HEAD FLOWS AND VOLUMES IN VBVL TABLE C17-----FOR INCLUSION IN BUDGET. PUT LABELS IN VBNM TABLE. CIN=CHIN COUT=CHOUT VBVL(1,MSUM)=VBVL(1,MSUM)+CIN*DELT VBVL(2,MSUM)=VBVL(2,MSUM)+COUT*DELT VBVL(3,MSUM)=CIN VBVL(4,MSUM)=COUT VBNM(MSUM)=TEXT MSUM=MSUM+1 C C18-----RETURN. RETURN END SUBROUTINE VDF2LPF7BDADJ(KSTP,KPER,IDIR,IBDRET, 1 IC1,IC2,IR1,IR2,IL1,IL2,IGRID) C ****************************************************************** C COMPUTE FLOW BETWEEN ADJACENT CELLS IN A SUBREGION OF THE GRID C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ USE GLOBAL, ONLY:NCOL,NROW,NLAY,IBOUND,HNEW,BUFF,CR,CC,CV, 1 BOTM,LBOTM,IOUT,IUNIT,DELR,DELC USE M_MF2005_IU, ONLY : IUANI USE GWFBASMODULE,ONLY:ICBCFL,DELT,PERTIM,TOTIM,ICHFLG USE GWFLPFMODULE,ONLY:ILPFCB,LAYTYP,NOVFC USE VDFMODULE, ONLY: IWTABLE,DENSEREF,PS,ELEV,HSALT C CHARACTER*16 TEXT(3) DOUBLE PRECISION HD,H1,HS1,H2,HS2 C DATA TEXT(1),TEXT(2),TEXT(3) 1 /'FLOW RIGHT FACE ','FLOW FRONT FACE ','FLOW LOWER FACE '/ C ------------------------------------------------------------------ C CALL SGWF2LPF7PNT(IGRID) C C1------IF CELL-BY-CELL FLOWS WILL BE SAVED IN A FILE, SET FLAG IBD. C1------RETURN IF FLOWS ARE NOT BEING SAVED OR RETURNED. ZERO=0. IBD=0 IF(ILPFCB.GT.0) IBD=ICBCFL IF(IBD.EQ.0 .AND. IBDRET.EQ.0) RETURN C C2------SET THE SUBREGION EQUAL TO THE ENTIRE GRID IF VALUES ARE BEING C2------SAVED IN A FILE. IF(IBD.NE.0) THEN K1=1 K2=NLAY I1=1 I2=NROW J1=1 J2=NCOL END IF C C3------TEST FOR DIRECTION OF CALCULATION; IF NOT ACROSS COLUMNS, GO TO C3------STEP 4. IF ONLY 1 COLUMN, RETURN. IF(IDIR.NE.1) GO TO 405 IF(NCOL.EQ.1) RETURN C C3A-----CALCULATE FLOW ACROSS COLUMNS (THROUGH RIGHT FACE). IF NOT C3A-----SAVING IN A FILE, SET THE SUBREGION. CLEAR THE BUFFER. IF(IBD.EQ.0) THEN K1=IL1 K2=IL2 I1=IR1 I2=IR2 J1=IC1-1 IF(J1.LT.1) J1=1 J2=IC2 END IF DO 310 K=K1,K2 DO 310 I=I1,I2 DO 310 J=J1,J2 BUFF(J,I,K)=ZERO 310 CONTINUE C C3B-----FOR EACH CELL CALCULATE FLOW THRU RIGHT FACE & STORE IN BUFFER. IF(J2.EQ.NCOL) J2=J2-1 DO 400 K=K1,K2 DO 400 I=I1,I2 DO 400 J=J1,J2 IF(ICHFLG.EQ.0) THEN IF((IBOUND(J,I,K).LE.0) .AND. (IBOUND(J+1,I,K).LE.0)) GO TO 400 ELSE IF((IBOUND(J,I,K).EQ.0) .OR. (IBOUND(J+1,I,K).EQ.0)) GO TO 400 END IF C--SEAWAT: WATER TABLE CORRECTIONS H1=HNEW(J,I,K) HS1=HSALT(J,I,K) H2=HNEW(J+1,I,K) HS2=HSALT(J+1,I,K) Z1=ELEV(J,I,K) Z2=ELEV(J+1,I,K) PS1=PS(J,I,K) PS2=PS(J+1,I,K) TOP1=BOTM(J,I,K-1) TOP2=BOTM(J+1,I,K-1) BOT1=BOTM(J,I,K) BOT2=BOTM(J+1,I,K) IF(IWTABLE.EQ.1.AND.LAYTYP(K).NE.0) & CALL VDWTABLE(H1,HS1,PS1,Z1,TOP1,BOT1, & H2,HS2,PS2,Z2,TOP2,BOT2) C--SEAWAT: END WATER TABLE CORRECTION PROGRAMMING DIS1=DELR(J+1)/2 DIS2=DELR(J)/2 AVGDENS=(DIS1*PS(J+1,I,K)+DIS2*PS(J,I,K))/(DIS1+DIS2) D=CR(J,I,K)*(AVGDENS-DENSEREF)/DENSEREF*(Z1-Z2) HDIFF=H1-H2 BUFF(J,I,K)=HDIFF*CR(J,I,K)+D 400 CONTINUE if(IUNIT(IUANI).gt.0) call gwf2ani7bd(igrid,idir) ! DLT C C3C-----RECORD CONTENTS OF BUFFER AND RETURN. IF(IBD.EQ.1) 1 CALL UBUDSV(KSTP,KPER,TEXT(1),ILPFCB,BUFF,NCOL,NROW,NLAY,IOUT) IF(IBD.EQ.2) CALL UBDSV1(KSTP,KPER,TEXT(1),ILPFCB,BUFF,NCOL,NROW, 1 NLAY,IOUT,DELT,PERTIM,TOTIM,IBOUND) RETURN C C4------TEST FOR DIRECTION OF CALCULATION; IF NOT ACROSS ROWS, GO TO C4------STEP 5. IF ONLY 1 ROW, RETURN. 405 IF(IDIR.NE.2) GO TO 505 IF(NROW.EQ.1) RETURN C C4A-----CALCULATE FLOW ACROSS ROWS (THROUGH FRONT FACE). IF NOT SAVING C4A-----IN A FILE, SET THE SUBREGION. CLEAR THE BUFFER. IF(IBD.EQ.0) THEN K1=IL1 K2=IL2 I1=IR1-1 IF(I1.LT.1) I1=1 I2=IR2 J1=IC1 J2=IC2 END IF DO 410 K=K1,K2 DO 410 I=I1,I2 DO 410 J=J1,J2 BUFF(J,I,K)=ZERO 410 CONTINUE C C4B-----FOR EACH CELL CALCULATE FLOW THRU FRONT FACE & STORE IN BUFFER. IF(I2.EQ.NROW) I2=I2-1 DO 500 K=K1,K2 DO 500 I=I1,I2 DO 500 J=J1,J2 IF(ICHFLG.EQ.0) THEN IF((IBOUND(J,I,K).LE.0) .AND. (IBOUND(J,I+1,K).LE.0)) GO TO 500 ELSE IF((IBOUND(J,I,K).EQ.0) .OR. (IBOUND(J,I+1,K).EQ.0)) GO TO 500 END IF C--SEAWAT WATER TABLE CORRECTIONS H1=HNEW(J,I,K) HS1=HSALT(J,I,K) H2=HNEW(J,I+1,K) HS2=HSALT(J,I+1,K) Z1=ELEV(J,I,K) Z2=ELEV(J,I+1,K) PS1=PS(J,I,K) PS2=PS(J,I+1,K) TOP1=BOTM(J,I,K-1) TOP2=BOTM(J,I+1,K-1) BOT1=BOTM(J,I,K) BOT2=BOTM(J,I+1,K) IF(IWTABLE.EQ.1.AND.LAYTYP(K).NE.0) & CALL VDWTABLE(H1,HS1,PS1,Z1,TOP1,BOT1, & H2,HS2,PS2,Z2,TOP2,BOT2) C--SEAWAT: END WATER TABLE CORRECTION PROGRAMMING DIS1=DELC(I+1)/2 DIS2=DELC(I)/2 AVGDENS=(DIS1*PS(J,I+1,K)+DIS2*PS(J,I,K))/(DIS1+DIS2) D=CC(J,I,K)*(AVGDENS-DENSEREF)/DENSEREF*(Z1-Z2) HDIFF=H1-H2 BUFF(J,I,K)=HDIFF*CC(J,I,K)+D 500 CONTINUE if(IUNIT(IUANI).gt.0) call gwf2ani7bd(igrid,idir) ! DLT C C4C-----RECORD CONTENTS OF BUFFER AND RETURN. IF(IBD.EQ.1) 1 CALL UBUDSV(KSTP,KPER,TEXT(2),ILPFCB,BUFF,NCOL,NROW,NLAY,IOUT) IF(IBD.EQ.2) CALL UBDSV1(KSTP,KPER,TEXT(2),ILPFCB,BUFF,NCOL,NROW, 1 NLAY,IOUT,DELT,PERTIM,TOTIM,IBOUND) RETURN C C5------DIRECTION OF CALCULATION IS ACROSS LAYERS BY ELIMINATION. IF C5------ONLY 1 LAYER, RETURN. 505 IF(NLAY.EQ.1) RETURN C C5A-----CALCULATE FLOW ACROSS LAYERS (THROUGH LOWER FACE). IF NOT C5A-----SAVING IN A FILE, SET THE SUBREGION. CLEAR THE BUFFER. IF(IBD.EQ.0) THEN K1=IL1-1 IF(K1.LT.1) K1=1 K2=IL2 I1=IR1 I2=IR2 J1=IC1 J2=IC2 END IF DO 510 K=K1,K2 DO 510 I=I1,I2 DO 510 J=J1,J2 BUFF(J,I,K)=ZERO 510 CONTINUE C C5B-----FOR EACH CELL CALCULATE FLOW THRU LOWER FACE & STORE IN BUFFER. IF(K2.EQ.NLAY) K2=K2-1 DO 600 K=1,K2 IF(K.LT.K1) GO TO 600 DO 590 I=I1,I2 DO 590 J=J1,J2 IF(ICHFLG.EQ.0) THEN IF((IBOUND(J,I,K).LE.0) .AND. (IBOUND(J,I,K+1).LE.0)) GO TO 590 ELSE IF((IBOUND(J,I,K).EQ.0) .OR. (IBOUND(J,I,K+1).EQ.0)) GO TO 590 END IF HD=HNEW(J,I,K+1) IF(NOVFC.NE.0 .OR. LAYTYP(K+1).EQ.0) GO TO 580 TMP=HD TOP=BOTM(J,I,LBOTM(K+1)-1) IF(TMP.LT.TOP) HD=TOP 580 HDIFF=HNEW(J,I,K)-HD DIS1=BOTM(J,I,K)-ELEV(J,I,K+1) DIS2=ELEV(J,I,K)-BOTM(J,I,K) AVGDENS=(DIS1*PS(J,I,K+1)+DIS2*PS(J,I,K))/(DIS1+DIS2) D=CV(J,I,K)*(AVGDENS-DENSEREF)/DENSEREF* & (ELEV(J,I,K)-ELEV(J,I,K+1)) BUFF(J,I,K)=HDIFF*CV(J,I,K)+D C--SEAWAT: CHECK AND RECALCULATE FOR DEWATERED CASE IF(LAYTYP(K+1).GT.0) THEN HS2=SALTHEAD(HNEW(J,I,K+1),PS(J,I,K+1),ELEV(J,I,K+1)) BOT1=BOTM(J,I,LBOTM(K+1)-1) IF(HS2.LT.BOT1) THEN HS1=SALTHEAD(HNEW(J,I,K),PS(J,I,K),ELEV(J,I,K)) BUFF(J,I,K)=PS(J,I,K)/DENSEREF*CV(J,I,K)*(HS1-BOT1) ENDIF ENDIF C--SEAWAT: END DEWATERED CORRECTION 590 CONTINUE 600 CONTINUE C C5C-----RECORD CONTENTS OF BUFFER AND RETURN. IF(IBD.EQ.1) 1 CALL UBUDSV(KSTP,KPER,TEXT(3),ILPFCB,BUFF,NCOL,NROW,NLAY,IOUT) IF(IBD.EQ.2) CALL UBDSV1(KSTP,KPER,TEXT(3),ILPFCB,BUFF,NCOL,NROW, 1 NLAY,IOUT,DELT,PERTIM,TOTIM,IBOUND) RETURN END