c Copyright (C) Stichting Deltares, 2005-2019. c c This file is part iMOD-WQ, which itself is part of iMOD. c c This program is free software: you can redistribute it and/or modify c it under the terms of the GNU General Public License as published by c the Free Software Foundation, either version 3 of the License, or c (at your option) any later version. c c This program is distributed in the hope that it will be useful, c but WITHOUT ANY WARRANTY; without even the implied warranty of c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the c GNU General Public License for more details. c c You should have received a copy of the GNU General Public License c along with this program. If not, see . c c Contact: imod.support@deltares.nl c Stichting Deltares c P.O. Box 177 c 2600 MH Delft, The Netherlands. c c iMOD-WQ is partly based on the following source codes: c 1. MT3DMS 5.3 source code of the University of Alabama. c The original MT3DMS 5.3 source code can be downloaded from the c university of Alabama website https://hydro.geo.ua.edu/mt3d/htm. c 2. SEAWAT v4 source code of the USGS. c The original SEAWAT v4 source code can be downloaded from the c USGS website https://www.usgs.gov/software/seawat-a-computer- c program-simulation-three-dimensional-variable-density-ground-water-flow c 3. RT3D2.5. The original RT3D2.5 source code can be downloaded from c http://tpclement.weebly.com/computer-tools.html. c For iMOD-WQ these source codes have been expanded c and extensively modified by Stichting Deltares. SUBROUTINE VDF2BCF7FM(KITER,KSTP,KPER,IGRID) C-----VERSION 24JAN2020 VDF2BCF6FM C ****************************************************************** C ADD LEAKAGE CORRECTION AND STORAGE TO HCOF AND RHS, AND CALCULATE C CONDUCTANCE AS REQUIRED C--SEAWAT: ADD UTIITY FOR VARIABLE DENSITY FLOW C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ USE VDFMODULE, ONLY: DENSEREF,PS,ELEV,HSALT,MFNADVFD USE GLOBAL, ONLY:IOUT,NCOL,NROW,NLAY,IBOUND,BOTM,NBOTM, 1 LBOTM,CV,HNEW,RHS,HCOF,HOLD,ISSFLG USE GWFBASMODULE,ONLY:DELT USE GWFBCFMODULE,ONLY:LAYCON,SC1,SC2 C--SEAWAT: ADD HTMP,HN,HNFE TO DOUBLE PRECISION DOUBLE PRECISION HTMP,HN,HNFE C ------------------------------------------------------------------ ISS=ISSFLG(KPER) KB=0 KT=0 ONE=1. IF(ISS.EQ.0) TLED=ONE/DELT C C1------FOR EACH LAYER: IF T VARIES CALCULATE HORIZONTAL CONDUCTANCES DO 100 K=1,NLAY KK=K C C1A-----IF LAYER TYPE IS NOT 1 OR 3 THEN SKIP THIS LAYER. IF(LAYCON(K).NE.3 .AND. LAYCON(K).NE.1) GO TO 100 KB=KB+1 C C1B-----FOR LAYER TYPES 1 & 3 CALL SGWF1BCF6H TO CALCULATE C1B-----HORIZONTAL CONDUCTANCES. CALL SVDF2BCF7H(HNEW,IBOUND,CR,CC,CV,HY,TRPY,DELR,DELC,BOTM,NBOTM, 1 KK,KB,KITER,KSTP,KPER,NCOL,NROW,NLAY,IOUT,WETDRY,IWDFLG, 2 CVWD,WETFCT,IWETIT,IHDWET,HDRY,BUFF,PS,ELEV,HSALT) 100 CONTINUE C C2------IF THE SIMULATION IS TRANSIENT ADD STORAGE TO HCOF AND RHS IF(ISS.NE.0) GO TO 201 KT=0 DO 200 K=1,NLAY C C3------SEE IF THIS LAYER IS CONVERTIBLE OR NON-CONVERTIBLE. IF(LAYCON(K).EQ.3 .OR. LAYCON(K).EQ.2) GO TO 150 C4------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 HCOF(J,I,K)=HCOF(J,I,K)-RHO*PS(J,I,K) C RHS(J,I,K)=RHS(J,I,K)-RHO*HOLD(J,I,K) RHS(J,I,K)=RHS(J,I,K)-RHO*HOLD(J,I,K)*PS(J,I,K) 140 CONTINUE GO TO 200 C C5------A CONVERTIBLE LAYER, SO CHECK OLD AND NEW HEADS TO DETERMINE C5------WHEN TO USE PRIMARY AND SECONDARY STORAGE 150 KT=KT+1 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 C TP=TOP(J,I,KT) TP=BOTM(J,I,LBOTM(K)-1) RHO2=SC2(J,I,KT)*TLED RHO1=SC1(J,I,K)*TLED C C5B-----FIND STORAGE FACTOR AT START OF TIME STEP. SOLD=RHO2 C--SEAWAT: MAKE COMPARISON WITH SALTHEAD C IF(HOLD(J,I,K).GT.TP) SOLD=RHO1 IF(SALTHEAD(HOLD(J,I,K),PS(J,I,K),ELEV(J,I,K)).GT.TP) SOLD=RHO1 C C5C-----FIND STORAGE FACTOR AT END OF TIME STEP. HTMP=HNEW(J,I,K) SNEW=RHO2 C--SEAWAT: MAKE COMP. WITH SALTHEAD C IF(HTMP.GT.TP) SNEW=RHO1 IF(SALTHEAD(HTMP,PS(J,I,K),ELEV(J,I,K)).GT.TP) SNEW=RHO1 C C5D-----ADD STORAGE TERMS TO RHS AND HCOF. C--SEAWAT: CONSERVE MASS C HCOF(J,I,K)=HCOF(J,I,K)-SNEW HCOF(J,I,K)=HCOF(J,I,K)-SNEW*PS(J,I,K) C RHS(J,I,K)=RHS(J,I,K) - SOLD*(HOLD(J,I,K)-TP) - SNEW*TP RHS(J,I,K)=RHS(J,I,K) - PS(J,I,K)*SOLD*(HOLD(J,I,K)-TP) & - PS(J,I,K)*SNEW*TP C 180 CONTINUE C 200 CONTINUE C C6------FOR EACH LAYER DETERMINE IF CORRECTION TERMS ARE NEEDED FOR C6------FLOW DOWN INTO PARTIALLY SATURATED LAYERS. 201 DO 300 K=1,NLAY C C7------SEE IF CORRECTION IS NEEDED FOR LEAKAGE FROM ABOVE. IF(LAYCON(K).NE.3 .AND. LAYCON(K).NE.2) GO TO 250 IF(K.EQ.1) GO TO 250 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 HTMP=HNEW(J,I,K) C C7C-----IF HEAD IS ABOVE TOP THEN CORRECTION NOT NEEDED C--SEAWAT: USE SALTHEAD C IF(HTMP.GE.BOTM(J,I,LBOTM(K)-1)) GO TO 220 IF(SALTHEAD(HTMP,PS(J,I,K),ELEV(J,I,K)).GE.BOTM(J,I,LBOTM(K)-1)) & GO TO 220 C C7D-----WITH HEAD BELOW TOP ADD CORRECTION TERMS TO RHS. C--SEAWAT: USE VD FORM OF DARCY LAW, CONSERVE MASS C RHS(J,I,K)=RHS(J,I,K) + CV(J,I,K-1)*(BOTM(J,I,LBOTM(K)-1)-HTMP) C RHS(J,I,K)=RHS(J,I,K) + PS(J,I,K-1)*CV(J,I,K-1)* C & (ELEV(J,I,K)-HTMP+PS(J,I,K-1)/DENSEREF* C & (BOTM(J,I,LBOTM(K)-1)-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 C C8------SEE IF THIS LAYER MAY NEED CORRECTION FOR LEAKAGE TO BELOW. 250 IF(K.EQ.NLAY) GO TO 300 IF(LAYCON(K+1).NE.3 .AND. LAYCON(K+1).NE.2) GO TO 300 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) C--SEAWAT: CONSERVE MASS C IF(HTMP.LT.BOTM(J,I,LBOTM(K+1)-1)) RHS(J,I,K)=RHS(J,I,K) C 1 - CV(J,I,K)*(BOTM(J,I,LBOTM(K+1)-1)-HTMP) C IF(SALTHEAD(HTMP,PS(J,I,K+1),ELEV(J,I,K+1)).LT. C & BOTM(J,I,LBOTM(K+1)-1)) C & RHS(J,I,K)=RHS(J,I,K) + PS(J,I,K)* C & CV(J,I,K)*(HTMP-ELEV(J,I,K+1)+PS(J,I,K)/DENSEREF C & *(ELEV(J,I,K+1)-BOTM(J,I,LBOTM(K+1)-1))) C C--SEAWAT: THIS SECTION REWRITTEN 9/24/07 USING CLEARER CODING 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 300 CONTINUE C C9------RETURN RETURN !rhs(45,24,3) END SUBROUTINE SVDF2BCF7H(HNEW,IBOUND,CR,CC,CV,HY,TRPY,DELR,DELC 1,BOTM,NBOTM,K,KB,KITER,KSTP,KPER,NCOL,NROW,NLAY,IOUT 2,WETDRY,IWDFLG,CVWD,WETFCT,IWETIT,IHDWET,HDRY,BUFF,PS,ELEV,HSALT) C-----VERSION 11JAN2000 SGWF1BCF6H C ****************************************************************** C COMPUTE CONDUCTANCE FOR ONE LAYER FROM SATURATED THICKNESS AND C HYDRAULIC CONDUCTIVITY C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ DOUBLE PRECISION HNEW,HD,BBOT,TTOP,HSALT C DIMENSION HNEW(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY) 1,CR(NCOL,NROW,NLAY), CC(NCOL,NROW,NLAY), CV(NCOL,NROW,NLAY) 2,HY(NCOL,NROW,NLAY), TRPY(NLAY), DELR(NCOL), DELC(NROW) 3,BOTM(NCOL,NROW,0:NBOTM),WETDRY(NCOL,NROW,NLAY) 4,CVWD(NCOL,NROW,NLAY),BUFF(NCOL,NROW,NLAY) C--SEAWAT: ADD THESE DIMENSION PS(NCOL,NROW,NLAY),ELEV(NCOL,NROW,NLAY), 1 HSALT(NCOL,NROW,NLAY) CHARACTER*3 ACNVRT DIMENSION ICNVRT(5),JCNVRT(5),ACNVRT(5) C COMMON /DISCOM/LBOTM(999),LAYCBD(999) COMMON /BCFCOM/LAYCON(999) COMMON /FLWAVG/LAYAVG(999) C ------------------------------------------------------------------ C C1------LOOP THROUGH EACH CELL IN LAYER AND CALCULATE TRANSMISSIVITY AT C1------EACH ACTIVE CELL. ZERO=0. NCNVRT=0 IHDCNV=0 ITFLG=1 IF(IWDFLG.NE.0) ITFLG=MOD(KITER,IWETIT) DO 200 I=1,NROW DO 200 J=1,NCOL C C2------IF CELL IS ACTIVE, THEN SKIP TO CODE THAT CALCULATES SATURATED C2------THICKNESS. IF(IBOUND(J,I,K).NE.0) GO TO 20 C C3------DETERMINE IF THE CELL CAN CONVERT BETWEEN CONFINED AND C3------UNCONFINED. IF NOT, SKIP TO CODE THAT SETS TRANSMISSIVITY TO 0. IF(ITFLG.NE.0) GO TO 6 IF(WETDRY(J,I,KB).EQ.ZERO)GO TO 6 WD=WETDRY(J,I,KB) IF(WD.LT.ZERO) WD=-WD C TURNON=BOT(J,I,KB)+WD TURNON=BOTM(J,I,LBOTM(K))+WD C C3A-----CHECK HEAD IN CELL BELOW TO SEE IF WETTING THRESHOLD HAS BEEN C3A-----REACHED. IF(K.EQ.NLAY)GO TO 2 C--SEAWAT: USE SALTHEAD HTMP=SALTHEAD(HNEW(J,I,K+1),PS(J,I,K+1),ELEV(J,I,K+1)) IF(IBOUND(J,I,K+1).GT.0.AND.HTMP.GE.TURNON)GO TO 9 C C3B-----CHECK HEAD IN ADJACENT HORIZONTAL CELLS TO SEE IF WETTING C3B-----THRESHOLD HAS BEEN REACHED. 2 IF(WETDRY(J,I,KB).LT.ZERO) GO TO 6 IF(J.EQ.1)GO TO 3 C--SEAWAT: USE SALTHEAD HTMP=SALTHEAD(HNEW(J-1,I,K),PS(J-1,I,K),ELEV(J-1,I,K)) IF(IBOUND(J-1,I,K).GT.0.AND.IBOUND(J-1,I,K).NE.30000.AND. 1 HTMP.GE.TURNON)GO TO 9 3 IF(J.EQ.NCOL)GO TO 4 C--SEAWAT: USE SALTHEAD HTMP=SALTHEAD(HNEW(J+1,I,K),PS(J+1,I,K),ELEV(J+1,I,K)) IF(IBOUND(J+1,I,K).GT.0.AND.HTMP.GE.TURNON)GO TO 9 4 IF(I.EQ.1)GO TO 5 C--SEAWAT: USE SALTHEAD HTMP=SALTHEAD(HNEW(J,I-1,K),PS(J,I-1,K),ELEV(J,I-1,K)) IF(IBOUND(J,I-1,K).GT.0.AND.IBOUND(J,I-1,K).NE.30000.AND. 1 HTMP.GE.TURNON)GO TO 9 5 IF(I.EQ.NROW)GO TO 6 C--SEAWAT: USE SALTHEAD HTMP=SALTHEAD(HNEW(J,I+1,K),PS(J,I+1,K),ELEV(J,I+1,K)) IF(IBOUND(J,I+1,K).GT.0.AND.HTMP.GE.TURNON)GO TO 9 C C3C-----CELL IS DRY AND STAYS DRY. SET TRANSMISSIVITY TO 0, SET C3C-----SATURATED THICKNESS (BUFF) TO 0, AND SKIP TO THE NEXT CELL. 6 CC(J,I,K)=ZERO IF(LAYAVG(K).EQ.30) BUFF(J,I,K)=ZERO GO TO 200 C C4------CELL BECOMES WET. SET INITIAL HEAD AND VERTICAL CONDUCTANCE. C 9 IF(IHDWET.NE.0) HNEW(J,I,K)=BOT(J,I,KB)+WETFCT*WD C IF(IHDWET.EQ.0) HNEW(J,I,K)=BOT(J,I,KB)+WETFCT*(HTMP-BOT(J,I,KB)) 9 IF(IHDWET.NE.0) HNEW(J,I,K)=BOTM(J,I,LBOTM(K))+WETFCT*WD IF(IHDWET.EQ.0) HNEW(J,I,K)=BOTM(J,I,LBOTM(K))+WETFCT* 1 (HTMP-BOTM(J,I,LBOTM(K))) IF(K.EQ.NLAY) GO TO 12 IF(IBOUND(J,I,K+1).NE.0) CV(J,I,K)= CVWD(J,I,K) 12 IF(K.EQ.1) GO TO 14 IF(IBOUND(J,I,K-1).NE.0) CV(J,I,K-1)= CVWD(J,I,K-1) 14 IBOUND(J,I,K)=30000 C--SEAWAT: RESET HSALT NOW THAT THE CELL HAS REWET HSALT(J,I,K)=SALTHEAD(HNEW(J,I,K),PS(J,I,K),ELEV(J,I,K)) C C4A-----PRINT MESSAGE SAYING CELL HAS BEEN CONVERTED TO WET. NCNVRT=NCNVRT+1 ICNVRT(NCNVRT)=I JCNVRT(NCNVRT)=J ACNVRT(NCNVRT)='WET' IF(NCNVRT.LT.5) GO TO 20 IF(IHDCNV.EQ.0) WRITE(IOUT,17) KITER,K,KSTP,KPER 17 FORMAT(1X,/1X,'CELL CONVERSIONS FOR ITER.=',I3,' LAYER=', 1 I3,' STEP=',I3,' PERIOD=',I3,' (ROW,COL)') IHDCNV=1 WRITE(IOUT,18) (ACNVRT(L),ICNVRT(L),JCNVRT(L),L=1,NCNVRT) 18 FORMAT(1X,3X,5(A,'(',I3,',',I3,') ')) NCNVRT=0 C C5------CALCULATE SATURATED THICKNESS. C--SEAWAT: USE SALTHEAD C 20 HD=HNEW(J,I,K) 20 HD=SALTHEAD(HNEW(J,I,K),PS(J,I,K),ELEV(J,I,K)) C BBOT=BOT(J,I,KB) BBOT=BOTM(J,I,LBOTM(K)) IF(LAYCON(K).EQ.1) GO TO 50 C TTOP=TOP(J,I,KT) TTOP=BOTM(J,I,LBOTM(K)-1) IF(BBOT.GT.TTOP) THEN WRITE(IOUT,35) K,I,J 35 FORMAT(1X,'Negative cell thickness at (Layer,row,col)', 1 I4,',',I4,',',I4) CALL USTOP(' ') END IF IF(HD.GT.TTOP) HD=TTOP 50 THCK=HD-BBOT C C6------CHECK TO SEE IF SATURATED THICKNESS IS GREATER THAN ZERO. IF(THCK.LE.ZERO) GO TO 100 C C6A-----IF SATURATED THICKNESS>0 THEN EITHER CALCULATE TRANSMISSIVITY C6A-----AS HYDRAULIC CONDUCTIVITY TIMES SATURATED THICKNESS OR STORE C6A-----K IN CC AND SATURATED THICKNESS IN BUFF. IF(LAYAVG(K).EQ.30) THEN CC(J,I,K)=HY(J,I,KB) BUFF(J,I,K)=THCK ELSE CC(J,I,K)=THCK*HY(J,I,KB) END IF GO TO 200 C C6B-----WHEN SATURATED THICKNESS < 0, PRINT A MESSAGE AND SET C6B-----TRANSMISSIVITY, IBOUND, AND VERTICAL CONDUCTANCE =0 100 NCNVRT=NCNVRT+1 ICNVRT(NCNVRT)=I JCNVRT(NCNVRT)=J ACNVRT(NCNVRT)='DRY' IF(NCNVRT.LT.5) GO TO 150 IF(IHDCNV.EQ.0) WRITE(IOUT,17) KITER,K,KSTP,KPER IHDCNV=1 WRITE(IOUT,18) (ACNVRT(L),ICNVRT(L),JCNVRT(L),L=1,NCNVRT) NCNVRT=0 150 HNEW(J,I,K)=HDRY CC(J,I,K)=ZERO IF(IBOUND(J,I,K).GE.0) GO TO 160 WRITE(IOUT,151) 151 FORMAT(1X,/1X,'CONSTANT-HEAD CELL WENT DRY', 1 ' -- SIMULATION ABORTED') WRITE(IOUT,152) K,I,J,KITER,KSTP,KPER 152 FORMAT(1X,'LAYER=',I2,' ROW=',I3,' COLUMN=',I3, 1 ' ITERATION=',I3,' TIME STEP=',I3,' STRESS PERIOD=',I3) CALL USTOP(' ') 160 IBOUND(J,I,K)=0 IF(K.LT.NLAY) CV(J,I,K)=ZERO IF(K.GT.1) CV(J,I,K-1)=ZERO 200 CONTINUE C C7------PRINT ANY REMAINING CELL CONVERSIONS NOT YET PRINTED IF(NCNVRT.EQ.0) GO TO 203 IF(IHDCNV.EQ.0) WRITE(IOUT,17) KITER,K,KSTP,KPER IHDCNV=1 WRITE(IOUT,18) (ACNVRT(L),ICNVRT(L),JCNVRT(L),L=1,NCNVRT) NCNVRT=0 C C8------CHANGE IBOUND VALUE FOR CELLS THAT CONVERTED TO WET THIS C8------ITERATION FROM 30000 to 1. 203 IF(IWDFLG.EQ.0) GO TO 210 DO 205 I=1,NROW DO 205 J=1,NCOL IF(IBOUND(J,I,K).EQ.30000) IBOUND(J,I,K)=1 205 CONTINUE C C9------COMPUTE HORIZONTAL BRANCH CONDUCTANCES FROM TRANSMISSIVITY. 210 IF(LAYAVG(K).EQ.0) THEN CALL SGWF2BCF7C(CR,CC,TRPY,DELR,DELC,K,NCOL,NROW,NLAY) ELSE IF(LAYAVG(K).EQ.10) THEN CALL SGWF2BCF7A(CR,CC,TRPY,DELR,DELC,K,NCOL,NROW,NLAY) ELSE IF(LAYAVG(K).EQ.20) THEN CALL SGWF2BCF7L(CR,CC,TRPY,DELR,DELC,K,NCOL,NROW,NLAY) ELSE CALL SGWF2BCF7U(CR,CC,TRPY,DELR,DELC,BUFF,K,NCOL,NROW,NLAY) END IF C C10-----RETURN. RETURN END SUBROUTINE VDF2BCF7BDS(KSTP,KPER,IGRID) C ****************************************************************** C COMPUTE STORAGE BUDGET FLOW TERM FOR BCF. 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 GWFBCFMODULE,ONLY:IBCFCB,LAYCON,SC1,SC2 USE VDFMODULE, ONLY: PS,ELEV,HSALT C CHARACTER*16 TEXT DOUBLE PRECISION STOIN,STOUT,SSTRG C--SEAWAT: ADD HTMP DOUBLE PRECISION HSING,HTMP C DATA TEXT /' STORAGE'/ C ------------------------------------------------------------------ CALL SGWF2BCF7PNT(IGRID) ISS=ISSFLG(KPER) C C1------INITIALIZE BUDGET ACCUMULATORS AND 1/DELT. 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(IBCFCB.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=LAYCON(K) IF(LC.EQ.3 .OR. LC.EQ.2) 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.NE.3 .AND. LC.NE.2) 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:COMPARE WITH SALTHEAD C IF(HOLD(J,I,K).GT.TP) SOLD=RHO1 IF(SALTHEAD(HOLD(J,I,K),PS(J,I,K),ELEV(J,I,K)).GT.TP) SOLD=RHO1 SNEW=RHO2 C--SEAWAT:COMPARE WITH SALTHEAD C IF(HSING.GT.TP) SNEW=RHO1 IF(SALTHEAD(HSING,PS(J,I,K),ELEV(J,I,K)).GT.TP) SNEW=RHO1 C--SEAWAT: CONSERVE MASS C STRG=SOLD*(HOLD(J,I,K)-TP) + SNEW*TP - SNEW*HSING 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 C STRG=RHO*HOLD(J,I,K) - RHO*HSING 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) 292,300,294 292 STOUT=STOUT-SSTRG GO TO 300 294 STOIN=STOIN+SSTRG 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 IBCFCB,BUFF,NCOL,NROW,NLAY,IOUT) IF(IBD.EQ.2) CALL UBDSV1(KSTP,KPER,TEXT,IBCFCB, 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 VDF2BCF7BDCH(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,DELC,DELR USE M_MF2005_IU, ONLY : IUANI USE GWFBASMODULE,ONLY:MSUM,VBVL,VBNM,DELT,PERTIM,TOTIM,ICBCFL, 1 ICHFLG USE GWFBCFMODULE,ONLY:IBCFCB,LAYCON USE VDFMODULE, ONLY: MFNADVFD,IWTABLE,DENSEREF,PS,ELEV,HSALT C CHARACTER*16 TEXT DOUBLE PRECISION HD,CHIN,CHOUT,XX1,XX2,XX3,XX4,XX5,XX6 C DATA TEXT /' CONSTANT HEAD'/ C ------------------------------------------------------------------ CALL SGWF2BCF7PNT(IGRID) C C1------SET IBD TO INDICATE IF CELL-BY-CELL BUDGET VALUES WILL BE SAVED. IBD=0 IF(IBCFCB.LT.0 .AND. ICBCFL.NE.0) IBD=-1 IF(IBCFCB.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,IBCFCB,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 LC=LAYCON(K) 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--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.LAYCON(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.LAYCON(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.LAYCON(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.LAYCON(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(LC.NE.3 .AND. LC.NE.2) GO TO 122 TMP=HD IF(TMP.LT.BOTM(J,I,LBOTM(K)-1)) HD=BOTM(J,I,LBOTM(K)-1) 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=1000.0; IF(DIS1+DIS2.NE.0.0)THEN AVGDENS=(DIS1*PS(J,I,K-1)+DIS2*PS(J,I,K))/(DIS1+DIS2) ENDIF 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(LAYCON(K).EQ.3 .OR. LAYCON(K).EQ.2) 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(LAYCON(K+1).NE.3 .AND. LAYCON(K+1).NE.2) GO TO 152 TMP=HD IF(TMP.LT.BOTM(J,I,LBOTM(K+1)-1)) HD=BOTM(J,I,LBOTM(K+1)-1) 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=1000.0; IF(DIS1+DIS2.NE.0.0)THEN AVGDENS=(DIS1*PS(J,I,K+1)+DIS2*PS(J,I,K))/(DIS1+DIS2) ENDIF 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--CHECK AND CORRECT FOR DEWATERED CASE IF(LAYCON(K+1).EQ.3 .OR. LAYCON(K+1).EQ.2) 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',I3,' STEP',I3) C--SEAWAT: WRITE FLOW RATHER THAN RATE WRITE(IOUT,900) K,I,J,FLOW 900 FORMAT(1X,'LAYER',I3,' ROW',I4,' COL',I4, 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(IBCFCB,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 IBCFCB,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 VDF2BCF7BDADJ(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,DELR,DELC USE GWFBASMODULE,ONLY:ICBCFL,DELT,PERTIM,TOTIM,ICHFLG USE GWFBCFMODULE,ONLY:IBCFCB,LAYCON, 1 seepage, seepagemv ! DLT use global, only: iunit ! DLT USE M_MF2005_IU 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 ------------------------------------------------------------------ CALL SGWF2BCF7PNT(IGRID) c c determine the seepage if (idir.eq.3 .and. nlay.gt.1) then ! DLT seepage = seepagemv ! DLT K=1 ! DLT DO 290 I=IR1,IR2 ! DLT DO 290 J=IC1,IC2 ! DLT IF(ICHFLG.EQ.0) THEN ! DLT IF((IBOUND(J,I,K).LE.0) .AND. (IBOUND(J,I,K+1).LE.0)) ! DLT 1 GO TO 290 ! DLT ELSE ! DLT IF((IBOUND(J,I,K).EQ.0) .OR. (IBOUND(J,I,K+1).EQ.0)) ! DLT 1 GO TO 290 ! DLT END IF ! DLT HD=HNEW(J,I,K+1) ! DLT IF(LAYCON(K+1).NE.3 .AND. LAYCON(K+1).NE.2) GO TO 280 ! DLT TMP=HD ! DLT IF(TMP.LT.BOTM(J,I,LBOTM(K+1)-1)) HD=BOTM(J,I,LBOTM(K+1)-1) ! DLT 280 HDIFF=HNEW(J,I,K)-HD ! DLT seepage(j,i) = hdiff*cv(j,i,k) ! DLT 290 CONTINUE ! DLT call sgwf2bcf7psv(igrid) ! DLT end if ! DLT 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. IBD=0 ZERO=0. IF(IBCFCB.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.LAYCON(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) C HDIFF=HNEW(J,I,K)-HNEW(J+1,I,K) HDIFF=H1-H2 BUFF(J,I,K)=CR(J,I,K)*HDIFF+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),IBCFCB,BUFF,NCOL,NROW,NLAY,IOUT) IF(IBD.EQ.2) CALL UBDSV1(KSTP,KPER,TEXT(1),IBCFCB,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.LAYCON(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) C HDIFF=HNEW(J,I,K)-HNEW(J,I+1,K) HDIFF=H1-H2 BUFF(J,I,K)=CC(J,I,K)*HDIFF+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),IBCFCB,BUFF,NCOL,NROW,NLAY,IOUT) IF(IBD.EQ.2) CALL UBDSV1(KSTP,KPER,TEXT(2),IBCFCB,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 continue ! DLT IF(NLAY.EQ.1) RETURN ! DLT 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(LAYCON(K+1).NE.3 .AND. LAYCON(K+1).NE.2) GO TO 580 TMP=HD IF(TMP.LT.BOTM(J,I,LBOTM(K+1)-1)) HD=BOTM(J,I,LBOTM(K+1)-1) 580 HDIFF=HNEW(J,I,K)-HD C 580 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=1000.0; IF(DIS1+DIS2.NE.0.0)THEN AVGDENS=(DIS1*PS(J,I,K+1)+DIS2*PS(J,I,K))/(DIS1+DIS2) ENDIF D=CV(J,I,K)*(AVGDENS-DENSEREF)/DENSEREF* & (ELEV(J,I,K)-ELEV(J,I,K+1)) BUFF(J,I,K)=CV(J,I,K)*HDIFF+D C--CHECK AND RECALCULATE FOR DEWATERED CASE IF(LAYCON(K+1).EQ.3 .OR. LAYCON(K+1).EQ.2) 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),IBCFCB,BUFF,NCOL,NROW,NLAY,IOUT) IF(IBD.EQ.2) CALL UBDSV1(KSTP,KPER,TEXT(3),IBCFCB,BUFF,NCOL,NROW, 1 NLAY,IOUT,DELT,PERTIM,TOTIM,IBOUND) RETURN END