MODULE GWFSUBMODULE INTEGER,SAVE,POINTER ::IIBSCB,ITMIN,NNDB,NDB,NMZ,NN,ND2,IDSAVE REAL, SAVE,POINTER ::AC1,AC2 LOGICAL,SAVE,POINTER ::NDF,NNDF INTEGER,SAVE, DIMENSION(:), POINTER ::ISBOCF,ISBOCU LOGICAL,SAVE, DIMENSION(:,:), POINTER ::OCFLGS LOGICAL,SAVE, DIMENSION(:), POINTER ::OCLAY INTEGER,SAVE, DIMENSION(:), POINTER ::ILSYS INTEGER,SAVE, DIMENSION(:), POINTER ::NTSSUM INTEGER,SAVE, DIMENSION(:), POINTER ::LN INTEGER,SAVE, DIMENSION(:), POINTER ::LDN INTEGER,SAVE, DIMENSION(:), POINTER ::NZ REAL, SAVE, DIMENSION(:), POINTER ::RNB REAL, SAVE, DIMENSION(:), POINTER ::DH REAL, SAVE, DIMENSION(:), POINTER ::DHP REAL, SAVE, DIMENSION(:), POINTER ::DHC REAL, SAVE, DIMENSION(:), POINTER ::DZ REAL, SAVE, DIMENSION(:), POINTER ::HC REAL, SAVE, DIMENSION(:), POINTER ::SCE REAL, SAVE, DIMENSION(:), POINTER ::SCV REAL, SAVE, DIMENSION(:), POINTER ::DCOM REAL, SAVE, DIMENSION(:), POINTER ::A1 REAL, SAVE, DIMENSION(:), POINTER ::A2 REAL, SAVE, DIMENSION(:), POINTER ::BB REAL, SAVE, DIMENSION(:), POINTER ::SUB REAL, SAVE, DIMENSION(:,:), POINTER ::DP REAL, SAVE, DIMENSION(:,:), POINTER ::DVB TYPE GWFSUBTYPE INTEGER, POINTER ::IIBSCB,ITMIN,NNDB,NDB,NMZ,NN,ND2,IDSAVE REAL, POINTER ::AC1,AC2 LOGICAL, POINTER ::NDF,NNDF INTEGER, DIMENSION(:), POINTER ::ISBOCF,ISBOCU LOGICAL, DIMENSION(:,:), POINTER ::OCFLGS LOGICAL, DIMENSION(:), POINTER ::OCLAY INTEGER, DIMENSION(:), POINTER ::ILSYS INTEGER, DIMENSION(:), POINTER ::NTSSUM INTEGER, DIMENSION(:), POINTER ::LN INTEGER, DIMENSION(:), POINTER ::LDN INTEGER, DIMENSION(:), POINTER ::NZ REAL, DIMENSION(:), POINTER ::RNB REAL, DIMENSION(:), POINTER ::DH REAL, DIMENSION(:), POINTER ::DHP REAL, DIMENSION(:), POINTER ::DHC REAL, DIMENSION(:), POINTER ::DZ REAL, DIMENSION(:), POINTER ::HC REAL, DIMENSION(:), POINTER ::SCE REAL, DIMENSION(:), POINTER ::SCV REAL, DIMENSION(:), POINTER ::DCOM REAL, DIMENSION(:), POINTER ::A1 REAL, DIMENSION(:), POINTER ::A2 REAL, DIMENSION(:), POINTER ::BB REAL, DIMENSION(:), POINTER ::SUB REAL, DIMENSION(:,:), POINTER ::DP REAL, DIMENSION(:,:), POINTER ::DVB END TYPE TYPE(GWFSUBTYPE), SAVE ::GWFSUBDAT(10) END MODULE GWFSUBMODULE SUBROUTINE GWF2SUB7AR(IN,IGRID) C ****************************************************************** C ALLOCATE ARRAY STORAGE FOR SUBSIDENCE PACKAGE. C READ SUBSIDENCE PACKAGE DATA. C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ USE GLOBAL, ONLY:IOUT,NCOL,NROW,NLAY,ISSFLG,NPER,NSTP,HNEW, 1 DELR,DELC,BUFF USE GWFSUBMODULE,ONLY:IIBSCB,ITMIN,NNDB,NDB,NMZ,NN,ND2,IDSAVE, 1 AC1,AC2,NDF,NNDF,ISBOCF,ISBOCU, 2 OCFLGS,OCLAY,ILSYS,NTSSUM,LN,LDN,NZ,RNB, 3 DH,DHP,DHC,DZ,HC,SCE,SCV,DCOM,A1,A2,BB, 4 SUB,DP,DVB C DIMENSION IFL(13) CHARACTER*24 ANAME(10) CHARACTER*200 LINE DATA ANAME(1) /' PRECONSOLIDATION HEAD'/ DATA ANAME(2) /'ELASTIC INTERBED STORAGE'/ DATA ANAME(3) /' VIRGIN INTERBED STORAGE'/ DATA ANAME(4) /' STARTING COMPACTION'/ DATA ANAME(5) /' DELAY STARTING HEAD'/ DATA ANAME(6) /' DELAY PRECOLSOL. HEAD'/ DATA ANAME(7) /'DELAY INITIAL COMPACTION'/ DATA ANAME(8) /'DELAY INTERBED THICKNESS'/ DATA ANAME(9) /' MATERIAL ZONE INDICES'/ DATA ANAME(10)/'NUMBER OF BEDS IN SYSTEM'/ DIMENSION IBUFF(NCOL,NROW) C ------------------------------------------------------------------ ALLOCATE (IIBSCB,ITMIN,NNDB,NDB,NMZ,NN,ND2,IDSAVE) ALLOCATE (AC1,AC2) ALLOCATE (NDF,NNDF) ALLOCATE (ISBOCF(6),ISBOCU(6)) ZERO=0.0 C C1------IDENTIFY PACKAGE. WRITE(IOUT,1)IN 1 FORMAT(/,'SUB7 -- SUBSIDENCE PACKAGE, VERSION 7,', 1 ' 03/31/2006',' INPUT READ FROM UNIT',I3) C C2------CHECK TO SEE THAT SUBSIDENCE OPTION IS APPROPRIATE C2------IF INAPPROPRIATE PRINT A MESSAGE & STOP THE SIMULATION. C2------ALSO, SUM TO GET THE TOTAL NUMBER OF TIME STEPS IN THE C2------SIMULATION. C NSTPT=0 DO 12 NS=1,NPER NSTPT=NSTPT+NSTP(NS) IF(ISSFLG(NS).NE.0.AND.NS.GT.1) THEN WRITE(IOUT,10) 10 FORMAT(1X,'SUBSIDENCE CANNOT BE USED IN SIMULATIONS', 1 ' IN WHICH STRESS PERIODS OTHER THAN THE ',/,1X, 2 ' FIRST ARE STEADY-STATE. SIMULATION ABORTED.') CALL USTOP(' ') ENDIF 12 CONTINUE C C3------ALLOCATE SPACE FOR ARRAY NTSSUM, WHICH WILL CONTAIN THE TOTAL C3------NUMBER OF TIME STEPS PRIOR TO THE CURRENT TIME STEP. ALLOCATE(NTSSUM(NPER)) C C4------READ FLAG FOR STORING CELL-BY-CELL STORAGE CHANGES AND C4------FLAG FOR PRINTING AND STORING COMPACTION, SUBSIDENCE, AND C4------CRITICAL HEAD ARRAYS. CALL URDCOM(IN,IOUT,LINE) LLOC=1 CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IIBSCB,R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISUBOC,R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NNDB,R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NDB,R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NMZ,R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NN,R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,3,I,AC1,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,3,I,AC2,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ITMIN,R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IDSAVE,R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IDREST,R,IOUT,IN) IF(AC2.EQ.ZERO) AC2=1.0 NDF=.TRUE. NNDF=.TRUE. IF(NNDB.LT.1) THEN NNDF=.FALSE. NNDB=0 ENDIF IF(NDB.LT.1) THEN NDF=.FALSE. NDB=0 NMZ=0 NN=0 ENDIF WRITE(IOUT,50) NNDB,NDB,NMZ,NN 50 FORMAT(/,' NUMBER OF SYSTEMS OF NO-DELAYED INTERBEDS:', 1 I3,/,' NUMBER OF SYSTEMS OF DELAY INTERBEDS:', 2 I3,/,' NUMBER OF MATERIAL ZONES:', 3 I3,/,' NUMBER OF NODES IN EACH STRING:',I3) IF(IDSAVE.GT.0) THEN WRITE(IOUT,52) IDSAVE 52 FORMAT(' RESTART INFORMATION WILL BE SAVED ON UNIT ', I5, 1 ' FOR DELAY INTERBEDS') ELSE WRITE(IOUT,53) 53 FORMAT(' RESTART INFORMATION WILL NOT BE SAVED FOR DELAY', 1 ' INTERBEDS') ENDIF IF(IDREST.GT.0) THEN WRITE(IOUT,54) IDREST 54 FORMAT(' RESTART INFORMATION WILL BE READ FROM UNIT ', I5, 1 ' FOR DELAY INTERBEDS') ELSE WRITE(IOUT,55) 55 FORMAT(' RESTART INFORMATION WILL NOT BE READ FOR DELAY', 1 ' INTERBEDS') ENDIF C C4A-----ABORT IF NO LAYERS ARE SPECIFIED FOR INTERBED STORAGE IF(.NOT.NNDF.AND..NOT.NDF) THEN WRITE(IOUT,60) 60 FORMAT(1X,'NO LAYERS WITH INTERBED STORAGE OF EITHER TYPE ', 1 'WERE SPECIFIED IN INPUT.',/,1X,'SIMULATION ABORTED.') CALL USTOP(' ') ENDIF C4B-----ABORT IF NO PROPERTY ZONES ARE SPECIFIED IF(NDF.AND.NMZ.LT.1) THEN WRITE(IOUT,*) ' STOPPING-- At least one property zone must ', & 'be specified for delay beds.' CALL USTOP(' ') ENDIF C4C-----ABORT IF NOT ENOUGH NODES ARE SPECIFIED IF(NDF.AND.NN.LT.2) THEN WRITE(IOUT,*) ' STOPPING-- Number of nodes in strings for ', & 'delay beds (NN) should be at least 2.' CALL USTOP(' ') ENDIF C C5------IF CELL-BY-CELL TERMS TO BE SAVED THEN PRINT UNIT NUMBER. 70 IF(IIBSCB.GT.0) WRITE(IOUT,80) IIBSCB 80 FORMAT(1X,'CELL-BY-CELL FLOW TERMS WILL BE SAVED ON UNIT',I3) C C5A-----IF OUTPUT CONTROL FOR PRINTING ARRAYS IS SELECTED PRINT MESSAGE. IF(ISUBOC.GT.0) WRITE(IOUT,90) 90 FORMAT(1X,'OUTPUT CONTROL RECORDS FOR SUB PACKAGE WILL BE ', 1 'READ EACH TIME STEP.') C C6------READ IN MODEL LAYER NUMBERS FOR EACH SYSTEM OF INTERBEDS, C6------FOR LAYERS WITHOUT DELAY. IF(NNDF) THEN ALLOCATE(LN(NNDB)) WRITE(IOUT,100) NNDB 100 FORMAT(/,' MODEL LAYER ASSIGNMENTS FOR EACH OF',I3,' NO-DELAY', 1 ' SYSTEMS OF INTERBEDS:') CALL URDCOM(IN,IOUT,LINE) READ(LINE,*) (LN(N),N=1,NNDB) WRITE(IOUT,115) (LN(N),N=1,NNDB) 115 FORMAT(1X,25I4) DO 120 N=1,NNDB IF(LN(N).GE.1.AND.LN(N).LE.NLAY) GO TO 120 WRITE(IOUT,118) 118 FORMAT(/,' IMPROPER LAYER ASSIGNMENT FOR NO-DELAY SYSTEM OF ', 1 'INTERBEDS.',/,' ABORTING...') CALL USTOP(' ') 120 CONTINUE ELSE ALLOCATE(LN(1)) ENDIF C C7------READ IN MODEL LAYER NUMBERS FOR EACH SYSTEM OF INTERBEDS, C7------FOR LAYERS WITH DELAY. IF(NDF) THEN ALLOCATE(LDN(NDB)) WRITE(IOUT,135) NDB 135 FORMAT(/,' MODEL LAYER ASSIGNMENTS FOR EACH OF',I3,' DELAY', 1 ' SYSTEMS OF INTERBEDS:') CALL URDCOM(IN,IOUT,LINE) READ(LINE,*) (LDN(N),N=1,NDB) WRITE(IOUT,115) (LDN(N),N=1,NDB) DO 140 N=1,NDB IF(LDN(N).GE.1.AND.LDN(N).LE.NLAY) GO TO 140 WRITE(IOUT,138) 138 FORMAT(/,' IMPROPER LAYER ASSIGNMENT FOR DELAY SYSTEM OF ', 1 'INTERBEDS.',/,' ABORTING...') CALL USTOP(' ') 140 CONTINUE ELSE ALLOCATE(LDN(1)) ENDIF C C8------ALLOCATE SPACE FOR THE ARRAYS HC, SCE, SCV, AND SUB. NCR=NROW*NCOL NND1=NCR*NNDB ND1=NCR*NDB ND2=0 C C9-----READ IN ARRAY RNB TO SEE HOW MANY STRINGS OF NN CELLS ARE NEEDED. IF(NDF) THEN ALLOCATE(RNB(ND1)) NNSUM=0 DO 190 KQ=1,NDB LOC1 = 1+(KQ-1)*NCR LAYNUM=LDN(KQ) WRITE(IOUT,144) KQ 144 FORMAT(/,1X,' SYSTEM',I4,' OF DELAY BEDS:') C CALL U2DREL(RNB(LOC1),ANAME(10),NROW,NCOL,LAYNUM,IN,IOUT) CALL U2DREL(BUFF(:,:,1),ANAME(10),NROW,NCOL,LAYNUM,IN,IOUT) CALL GWF2SUB72D1D(BUFF(:,:,1),NCOL,NROW,RNB,ND1,LOC1) DO 180 N=1,NCR IF(RNB(LOC1+N-1).GE.1.0) NNSUM=NNSUM+1 180 CONTINUE 190 CONTINUE ND2=NN*NNSUM ELSE ALLOCATE(RNB(1)) ENDIF IF(ND2.LT.1.AND.NDF) THEN WRITE(IOUT,*) ' STOPPING-- Delay beds were not found in ', & ' array specifying numbers of delay beds (RNB).' CALL USTOP(' ') ENDIF C C10-----ALLOCATE MEMORY. ALLOCATE(OCFLGS(13,NSTPT)) ALLOCATE(OCLAY(NLAY)) IF(NNDF) THEN ALLOCATE(HC(NND1)) ALLOCATE(SCE(NND1)) ALLOCATE(SCV(NND1)) ALLOCATE(SUB(NND1)) ALLOCATE(ILSYS(NNDB)) ELSE ALLOCATE(HC(1)) ALLOCATE(SCE(1)) ALLOCATE(SCV(1)) ALLOCATE(SUB(1)) ALLOCATE(ILSYS(1)) ENDIF IF(NDF) THEN ALLOCATE(NZ(ND1)) ALLOCATE(DZ(ND1)) ALLOCATE(DCOM(ND1)) ALLOCATE(DHP(ND2)) ALLOCATE(DH(ND2)) ALLOCATE(DHC(ND2)) ALLOCATE(DP(NMZ,3)) ALLOCATE(DVB(NDB,4)) ALLOCATE(A1(NN)) ALLOCATE(A2(NN)) ALLOCATE(BB(NN)) ELSE ALLOCATE(NZ(1)) ALLOCATE(DZ(1)) ALLOCATE(DCOM(1)) ALLOCATE(DHP(1)) ALLOCATE(DH(1)) ALLOCATE(DHC(1)) ALLOCATE(DP(1,1)) ALLOCATE(DVB(1,1)) ALLOCATE(A1(1)) ALLOCATE(A2(1)) ALLOCATE(BB(1)) ENDIF C C11-----READ ARRAYS. NCR=NROW*NCOL ANNI=0.5/(FLOAT(NN)-.5) C C12-----READ RESTART RECORDS IF THIS SIMULATION CONTINUES FROM A C12-----PREVIOUS SIMULATION IF(NDF) THEN IF(IDREST.GT.0) THEN READ(IDREST) NND2 IF(NND2.EQ.ND2) THEN WRITE(IOUT,242) 242 FORMAT(' HEAD AND PRECONSOLIDATION HEAD FOR DELAY BEDS ARE', 1 ' BEING READ FROM RESTART RECORDS') READ(IDREST) (DH(N),N=1,ND2) READ(IDREST) (DHC(N),N=1,ND2) DO 250 N2=1,ND2 DHP(N2)=DH(N2) 250 CONTINUE ELSE WRITE(IOUT,252) 252 FORMAT(' HEAD AND PRECONSOLIDATION HEAD FOR DELAY BEDS ', 1 'CANNOT BE READ FROM RESTART RECORDS',/, 2 ' SIMULATION ABORTING') CALL USTOP(' ') ENDIF ENDIF ENDIF C C13-----READ IN ARRAYS FOR SYSTEMS OF NO-DELAY INTERBEDS. IF(NNDF) THEN DO 260 KQ=1,NNDB K=LN(KQ) LOC1=1+(KQ-1)*NCR WRITE(IOUT,256) KQ 256 FORMAT(/,1X,' SYSTEM',I4,' OF NO-DELAY BEDS:') C CALL U2DREL(HC(LOC1),ANAME(1),NROW,NCOL,K,IN,IOUT) CALL U2DREL(BUFF(:,:,1),ANAME(1),NROW,NCOL,K,IN,IOUT) CALL GWF2SUB72D1D(BUFF(:,:,1),NCOL,NROW,HC,NND1,LOC1) WRITE(IOUT,256) KQ C CALL U2DREL(SCE(LOC1),ANAME(2),NROW,NCOL,K,IN,IOUT) CALL U2DREL(BUFF(:,:,1),ANAME(2),NROW,NCOL,K,IN,IOUT) CALL GWF2SUB72D1D(BUFF(:,:,1),NCOL,NROW,SCE,NND1,LOC1) WRITE(IOUT,256) KQ C CALL U2DREL(SCV(LOC1),ANAME(3),NROW,NCOL,K,IN,IOUT) CALL U2DREL(BUFF(:,:,1),ANAME(3),NROW,NCOL,K,IN,IOUT) CALL GWF2SUB72D1D(BUFF(:,:,1),NCOL,NROW,SCV,NND1,LOC1) WRITE(IOUT,256) KQ C CALL U2DREL(SUB(LOC1),ANAME(4),NROW,NCOL,K,IN,IOUT) CALL U2DREL(BUFF(:,:,1),ANAME(4),NROW,NCOL,K,IN,IOUT) CALL GWF2SUB72D1D(BUFF(:,:,1),NCOL,NROW,SUB,NND1,LOC1) 260 CONTINUE C C14-----INITIALIZE ARRAYS FOR SYSTEMS OF NO-DELAY INTERBEDS. DO 280 KQ=1,NNDB K=LN(KQ) NQ=(KQ-1)*NCR C NK=(K-1)*NCR DO 270 IR=1,NROW NQR=NQ+(IR-1)*NCOL C NKR=NK+(IR-1)*NCOL DO 270 IC=1,NCOL LOC2=NQR+IC C LOC2H=NKR+IC C C15------MULTIPLY STORAGE BY AREA TO GET STORAGE CAPACITY. AREA=DELR(IC)*DELC(IR) SCE(LOC2)=SCE(LOC2)*AREA SCV(LOC2)=SCV(LOC2)*AREA C C16-----MAKE SURE THAT PRECONSOLIDATION HEAD VALUES C16-----ARE CONSISTANT WITH STARTING HEADS. IF(HC(LOC2).GT.HNEW(IC,IR,K)) HC(LOC2)=HNEW(IC,IR,K) 270 CONTINUE 280 CONTINUE ENDIF IF(NDF) THEN C C17-----READ IN TABLE OF MATERIAL PROPERTIES: K, Sse, Ssv FOR EACH C17-----OF NMZ ZONES. WRITE(IOUT,295) 295 FORMAT(/,' MATERIAL PROPERTIES OF INTERBEDS WITH DELAY PROPERTIES' 1 ,//,' ZONE HYDRAULIC ELASTIC INEL', 2 'ASTIC ',/,' NUMBER CONDUCTIVITY SPECIFIC STORA', 3 'GE SPECIFIC STORAGE ',/,' ',69('-')) DO 300 N=1,NMZ READ(IN,*) (DP(N,NP),NP=1,3) 300 CONTINUE WRITE(IOUT,305) (N,(DP(N,NP),NP=1,3),N=1,NMZ) 305 FORMAT(I5,4X,G15.5,5X,G15.5,5X,G15.5) LOC3=0 LOC4=0 DO 380 KQ=1,NDB K=LDN(KQ) LOC1=1+(KQ-1)*NCR C C18-----READ IN ARRAYS FOR SYSTEMS OF DELAY INTERBEDS. IF(IDREST.LE.0) THEN WRITE(IOUT,308) KQ 308 FORMAT(/,1X,' SYSTEM',I4,' OF DELAY BEDS:') CALL U2DREL(BUFF(:,:,1),ANAME(5),NROW,NCOL,K,IN,IOUT) N1=0 DO 320 IR=1,NROW DO 320 IC=1,NCOL N1=N1+1 LOC2=LOC1+N1-1 IF(RNB(LOC2).LT.1.0) GO TO 320 DO 315 N2=1,NN LOC3=LOC3+1 DHP(LOC3)=BUFF(IC,IR,1) DH(LOC3)=BUFF(IC,IR,1) 315 CONTINUE 320 CONTINUE WRITE(IOUT,308) KQ CALL U2DREL(BUFF(:,:,1),ANAME(6),NROW,NCOL,K,IN,IOUT) N1=0 DO 330 IR=1,NROW DO 330 IC=1,NCOL N1=N1+1 LOC2=LOC1+N1-1 IF(RNB(LOC2).LT.1.0) GO TO 330 DO 325 N2=1,NN LOC4=LOC4+1 DHC(LOC4)=BUFF(IC,IR,1) IF(DHC(LOC4).GT.DH(LOC4)) DHC(LOC4)=DH(LOC4) 325 CONTINUE 330 CONTINUE ENDIF WRITE(IOUT,308) KQ C CALL U2DREL(DCOM(LOC1),ANAME(7),NROW,NCOL,K,IN,IOUT) CALL U2DREL(BUFF(:,:,1),ANAME(7),NROW,NCOL,K,IN,IOUT) CALL GWF2SUB72D1D(BUFF(:,:,1),NCOL,NROW,DCOM,ND1,LOC1) WRITE(IOUT,308) KQ C CALL U2DREL(DZ(LOC1),ANAME(8),NROW,NCOL,K,IN,IOUT) CALL U2DREL(BUFF(:,:,1),ANAME(8),NROW,NCOL,K,IN,IOUT) CALL GWF2SUB72D1D(BUFF(:,:,1),NCOL,NROW,DZ,ND1,LOC1) WRITE(IOUT,308) KQ C CALL U2DINT(NZ(LOC1),ANAME(9),NROW,NCOL,K,IN,IOUT) CALL U2DINT(IBUFF,ANAME(9),NROW,NCOL,K,IN,IOUT) L=LOC1-1 DO 335 I=1,NROW DO 335 J=1,NCOL L=L+1 NZ(L)=IBUFF(J,I) 335 CONTINUE C C19-----INITIALIZE ARRAYS FOR SYSTEMS OF DELAY INTERBEDS. DO 360 NL=1,NCR LOC2=LOC1+NL-1 IF(RNB(LOC2).GE.1.0.AND.DZ(LOC2).LE.ZERO) THEN WRITE(IOUT,355) 355 FORMAT(' A VALUE OF ZERO WAS FOUND IN THE DZ ARRAY WHERE ', 1 'DELAY INTERBEDS OCCUR.',/,' MAKE SURE THAT', 2 ' DZ IS GREATER THAN 0.0 AT ALL CELLS WHERE RNB ',/, 3 ' IS 1.0 OR MORE. SIMULATION ABORTING') CALL USTOP(' ') ENDIF DZ(LOC2)=DZ(LOC2)*ANNI 360 CONTINUE DO 370 N=1,4 DVB(KQ,N)=ZERO 370 CONTINUE 380 CONTINUE ENDIF C C20-----SET ALL FLAGS FOR OUTPUT CONTROL TO "FALSE". DO 390 I=1,NSTPT DO 385 N=1,13 OCFLGS(N,I)=.FALSE. 385 CONTINUE 390 CONTINUE C C21-----READ FORMATS AND UNIT NUMBERS OUTPUT FLAGS. IF(ISUBOC.GT.0) THEN CALL URDCOM(IN,IOUT,LINE) LLOC=1 CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCF(1),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCU(1),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCF(2),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCU(2),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCF(3),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCU(3),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCF(4),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCU(4),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCF(5),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCU(5),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCF(6),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCU(6),R,IOUT,IN) WRITE(IOUT,410) (ISBOCF(N),ISBOCU(N),N=1,6) 410 FORMAT(/,' SUBSIDENCE PRINT FORMAT IS NUMBER',I4/ & ' UNIT FOR SAVING SUBSIDENCE IS',I4/ & ' COMPACTION BY LAYER PRINT FORMAT IS NUMBER',I4/ & ' UNIT FOR SAVING COMPACTION BY LAYER IS',I4/ & ' COMPACTION BY SYSTEM PRINT FORMAT IS NUMBER',I4/ & ' UNIT FOR SAVING COMPACTION BY SYSTEM IS',I4/ & ' VERTICAL DISPLACEMENT PRINT FORMAT IS NUMBER',I4/ & ' UNIT FOR SAVING VERTICAL DISPLACEMENT IS',I4/ & ' NO-DELAY CRITICAL HEAD PRINT FORMAT IS NUMBER',I4/ & ' UNIT FOR SAVING NO-DELAY CRITICAL HEAD IS',I4/ & ' DELAY CRITICAL HEAD PRINT FORMAT IS NUMBER',I4/ & ' UNIT FOR SAVING DELAY CRITICAL HEAD IS',I4) NTSSUM(1)=0 IF(NPER.GT.1) THEN DO 415 N=2,NPER NTSSUM(N)=NTSSUM(N-1)+NSTP(N-1) 415 CONTINUE ENDIF DO 450 NOCLIN=1,ISUBOC CALL URDCOM(IN,IOUT,LINE) LLOC=1 CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISP1,R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISP2,R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,JTS1,R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,JTS2,R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(1),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(2),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(3),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(4),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(5),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(6),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(7),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(8),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(9),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(10),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(11),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(12),R,IOUT,IN) CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(13),R,IOUT,IN) IF(ISP1.LT.1) ISP1=1 IF(ISP1.GT.NPER) ISP1=NPER IF(ISP2.LT.1) ISP2=1 IF(ISP2.GT.NPER) ISP2=NPER IF(ISP1.GT.ISP2) ISP1=ISP2 DO 440 I=ISP1,ISP2 J1=JTS1 J2=JTS2 IF(J1.LT.1) J1=1 IF(J1.GT.NSTP(I)) J1=NSTP(I) IF(J2.LT.1) J2=1 IF(J2.GT.NSTP(I)) J2=NSTP(I) IF(J1.GT.J2) J1=J2 DO 430 J=J1,J2 ILOC=NTSSUM(I)+J DO 420 N=1,13 IF(IFL(N).GT.0) OCFLGS(N,ILOC)=.TRUE. IF(IFL(N).EQ.0) OCFLGS(N,ILOC)=.FALSE. 420 CONTINUE 430 CONTINUE 440 CONTINUE 450 CONTINUE ENDIF C C22-----RETURN 500 CALL SGWF2SUB7PSV(IGRID) RETURN END SUBROUTINE GWF2SUB7ST(KPER,IGRID) C ****************************************************************** C SET PRECONSOLIDATION HEAD (HC AND DHC) EQUAL TO THE STEADY- C STATE HEAD IF HEAD IS LOWER THAN PRECONSOLIDATION HEAD. C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ USE GLOBAL, ONLY:HNEW,NCOL,NROW,ISSFLG USE GWFSUBMODULE ,ONLY: RNB,LN,LDN,HC,DHP,DH,DHC,NZ,DZ, 1 NDF,NNDF,NNDB,NDB,NN C C ------------------------------------------------------------------ CALL SGWF2SUB7PNT(IGRID) C C1------RETURN IF THIS IS NOT THE SECOND STRESS PERIOD OR IF THE FIRST C1------STRESS PERIOD WAS TRANSIENT. IF(KPER.NE.2) RETURN IF(ISSFLG(1).EQ.0) RETURN C2-----MAKE SURE THAT NO-DELAY PRECONSOLIDATION HEAD VALUES ARE CONSISTENT C2-----WITH STEADY-STATE HEADS. NCR=NROW*NCOL IF(NNDF) THEN DO 20 KQ=1,NNDB K=LN(KQ) NQ=(KQ-1)*NCR C NK=(K-1)*NCR DO 10 IR=1,NROW NQR=NQ+(IR-1)*NCOL C NKR=NK+(IR-1)*NCOL DO 10 IC=1,NCOL LOC2=NQR+IC C LOC2H=NKR+IC HTMP=HNEW(IC,IR,K) IF(HC(LOC2).GT.HTMP) HC(LOC2)=HTMP 10 CONTINUE 20 CONTINUE ENDIF C3-----MAKE SURE THAT DELAY PRECONSOLIDATION HEAD VALUES ARE CONSISTENT C3-----WITH STEADY-STATE HEADS. ALSO, SET HEAD (DH) AND HEAD FOR C3-----PREVIOUS TIME STEP (DHP) EQUAL TO THE AQUIFER HEAD (HNEW). IF(NDF) THEN LOC3=0 DO 40 KQ=1,NDB K=LDN(KQ) NQ=(KQ-1)*NCR C NK=(K-1)*NCR N1=0 DO 30 IR=1,NROW DO 30 IC=1,NCOL N1=N1+1 LOC2=NQ+N1 C LOC2H=NK+N1 HTMP=HNEW(IC,IR,K) IF(RNB(LOC2).LT.1.0) GO TO 30 DO 18 N2=1,NN LOC3=LOC3+1 IF(DHC(LOC3).GT.HTMP) DHC(LOC3)=HTMP DH(LOC3)=HTMP DHP(LOC3)=HTMP 18 CONTINUE 30 CONTINUE 40 CONTINUE ENDIF C4-----RETURN. RETURN END SUBROUTINE GWF2SUB7FM(KPER,KITER,ISIP,IGRID) C ****************************************************************** C ADD INTERBED STORAGE TERMS TO RHS AND HCOF C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ USE GLOBAL, ONLY: RHS,HCOF,HNEW,HOLD,IBOUND,DELR,DELC, 1 NCOL,NROW,ISSFLG USE GWFBASMODULE, ONLY: DELT USE SIPMODULE, ONLY: V,HCLOSE USE GWFSUBMODULE ,ONLY: RNB,LN,LDN,HC,SCE,SCV,DHP,DH,DHC,NZ,DZ,DP, 1 BB,NDF,NNDF,AC1,AC2,ITMIN,NN, 1 NDB,NNDB,NMZ LOGICAL ICHK C DATA ICHK/.FALSE./ C ------------------------------------------------------------------ CALL SGWF2SUB7PNT(IGRID) ZERO=0.0 C C0------SKIP CALCULATIONS IF THIS IS A STEADY-STATE STRESS PERIOD. IF(ISSFLG(KPER).EQ.1) RETURN C C1------INITIALIZE TLED=1./DELT NCR=NCOL*NROW C IF(NNDF) THEN C2------FIND LAYERS WITH INTERBED STORAGE DO 110 KQ=1,NNDB K=LN(KQ) LOCT=(KQ-1)*NCR N=0 DO 100 IR=1,NROW DO 100 IC=1,NCOL N=N+1 LOC2=LOCT+N IF(IBOUND(IC,IR,K).LE.0) GO TO 100 C C3------DETERMINE STORAGE CAPACITIES FOR CELL AT START AND END OF STEP RHO1=SCE(LOC2)*TLED RHO2=RHO1 HCTMP=HC(LOC2) IF(HNEW(IC,IR,K).LT.HCTMP) RHO2=SCV(LOC2)*TLED C C4------ADD APPROPRIATE TERMS TO RHS AND HCOF RHS(IC,IR,K)=RHS(IC,IR,K)-HCTMP*(RHO2-RHO1)-RHO1*HOLD(IC,IR,K) HCOF(IC,IR,K)=HCOF(IC,IR,K)-RHO2 100 CONTINUE 110 CONTINUE ENDIF IF(NDF) THEN LOC3=1-NN DO 420 KQ=1,NDB K=LDN(KQ) NQ=(KQ-1)*NCR DO 410 IR=1,NROW NQR=NQ+(IR-1)*NCOL DO 410 IC=1,NCOL LOC2=NQR+IC RNB2=RNB(LOC2) IF(RNB2.LT.1.0) GO TO 410 LOC3=LOC3+NN IF(IBOUND(IC,IR,K).LE.0) GO TO 410 IF(ISIP.GT.0.AND.ICHK) THEN VV=V(IC,IR,K) ELSE VV=ZERO ICHK=.TRUE. ENDIF NZONE=NZ(LOC2) DZZ=DZ(LOC2) CI=DP(NZONE,1)/DZZ IF(ISIP.GT.0) THEN IF(KITER.GT.ITMIN.AND.ABS(VV).LT.HCLOSE) GO TO 205 END IF C C5------ASSEMBLE COEFFICIENTS FOR DIRECT SOLUTION OF HEAD IN INTERBED HAQ=HNEW(IC,IR,K)+VV*AC1 SSE=DP(NZONE,2) SSV=DP(NZONE,3) NEND=LOC3+NN-1 CALL SGWF2SUB7A(HAQ,TLED,CI,SSE,SSV,DZZ, 1 DH(LOC3:NEND),DHP(LOC3:NEND),DHC(LOC3:NEND),NN) C C6------SOLVE FOR HEAD CHANGES IN STRING USING GAUSSIAN ELIMINATION. C6------ADD CHANGES TO HEAD VALUES TO GET HEAD AT CURRENT ITERATION. CALL SGWF2SUB7S(NN) DO 200 N=1,NN DH(LOC3+N-1)=DH(LOC3+N-1)+BB(N)*AC2 200 CONTINUE C C7------CALCULATE STORAGE CHANGE IN INTERBEDS 205 AREA=DELR(IC)*DELC(IR) STRGS=ZERO L1=LOC3 L2=LOC3+NN-1 DO 210 LOC4=L1,L2 HHOLD=DHP(LOC4) HHNEW=DH(LOC4) HHC=DHC(LOC4) C C8------GET STORAGE CAPACITIES AT BEGINNING AND END OF TIME STEP. SBGN=DP(NZ(LOC2),2) SEND=SBGN IF(HHNEW.LT.HHC) SEND=DP(NZ(LOC2),3) C C9------CALCULATE VOLUME CHANGE IN INTERBED STORAGE FOR TIME STEP. STR1=(HHC*(SEND-SBGN)+SBGN*HHOLD- 1 SEND*HHNEW)*DZ(LOC2)*RNB(LOC2)*2. IF(LOC4.EQ.L2) STR1=STR1*.5 STRGS=STRGS+STR1 210 CONTINUE RATES=STRGS*AREA*TLED C C10-----ADD APPROPRIATE TERMS TO RHS AND HCOF RHS(IC,IR,K)=RHS(IC,IR,K)-RATES 410 CONTINUE 420 CONTINUE ENDIF C C11-----RETURN RETURN END SUBROUTINE GWF2SUB7BD(KSTP,KPER,IGRID) C ****************************************************************** C CALCULATE VOLUMETRIC BUDGET FOR INTERBED STORAGE C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ USE GLOBAL, ONLY: IOUT,NCOL,NROW,NLAY,IBOUND,HNEW,HOLD, 1 BUFF,DELR,DELC,ISSFLG USE GWFBASMODULE, ONLY: VBVL,VBNM,MSUM,ICBCFL,DELT USE GWFSUBMODULE ,ONLY: RNB,LN,LDN,HC,SCE,SCV,SUB,DHP,DH,DHC, 1 NZ,DZ,DCOM,DP,DVB,NDF,NNDF,NN,ND2,NDB, 2 NNDB,NMZ,IIBSCB CHARACTER*16 TEXT(2) C DATA TEXT(1) /'INST. IB STORAGE'/ DATA TEXT(2) /'DELAY IB STORAGE'/ C ------------------------------------------------------------------ CALL SGWF2SUB7PNT(IGRID) ZERO=0.0 TLED=ZERO IF(ISSFLG(KPER).EQ.0) TLED=1./DELT C C1------INITIALIZE CELL-BY-CELL FLOW TERM FLAG (IBD) AND C1------SET IF CELL-BY-CELL FLOW TERMS ARE NEEDED. IBD=0 IF(ICBCFL.NE.0 .AND. IIBSCB.GT.0 ) IBD=1 NCR=NCOL*NROW C C2------RUN THROUGH EVERY CELL IN THE GRID WITH INTERBED STORAGE. IF(NNDF) THEN C C3-------CELL-BY-CELL FLOW TERMS ARE NEEDED SET IBD AND CLEAR BUFFER. IF(IBD.EQ.1) THEN DO 90 K=1,NLAY DO 90 IR=1,NROW DO 90 IC=1,NCOL BUFF(IC,IR,K)=ZERO 90 CONTINUE ENDIF STOIN=ZERO STOUT=ZERO C C4------IF THIS IS A STEADY-STATE STRESS PERIOD, SKIP CALCULATIONS IF(ISSFLG(KPER).EQ.1) GO TO 111 C C5------CALCULATE NO-DELAY INTERBED STORAGE CHANGE FOR EACH CELL DO 110 KQ=1,NNDB K=LN(KQ) NQ=(KQ-1)*NCR DO 100 IR=1,NROW NQR=NQ+(IR-1)*NCOL DO 100 IC=1,NCOL LOC2=NQR+IC C C6------CALCULATE FLOW FROM STORAGE (VARIABLE HEAD CELLS ONLY) IF(IBOUND(IC,IR,K).LE.0) GO TO 100 HHOLD=HOLD(IC,IR,K) HHNEW=HNEW(IC,IR,K) HHC=HC(LOC2) C C7------GET STORAGE CAPACITIES AT BEGINNING AND END OF TIME STEP. SBGN=SCE(LOC2) SEND=SBGN IF(HHNEW.LT.HHC) SEND=SCV(LOC2) C C8------CALCULATE VOLUME CHANGE IN INTERBED STORAGE FOR TIME STEP. STRG=HHC*(SEND-SBGN)+SBGN*HHOLD-SEND*HHNEW C C9------ACCUMULATE SUBSIDENCE ASSOCIATED WITH CHANGE IN STORAGE SUB(LOC2)=SUB(LOC2)+STRG/(DELR(IC)*DELC(IR)) C C10-----IF C-B-C FLOW TERMS ARE TO BE SAVED THEN ADD RATE TO BUFFER. IF(IBD.EQ.1) BUFF(IC,IR,K)=BUFF(IC,IR,K)+STRG*TLED C C11-----SEE IF FLOW IS INTO OR OUT OF STORAGE. IF(STRG.LE.ZERO) THEN STOUT=STOUT-STRG ELSE STOIN=STOIN+STRG ENDIF 100 CONTINUE 110 CONTINUE C C12-----IF C-B-C FLOW TERMS WILL BE SAVED CALL UBUDSV TO RECORD THEM. 111 IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT(1),IIBSCB,BUFF,NCOL, 1 NROW,NLAY,IOUT) C C13-----MOVE RATES,VOLUMES & LABELS INTO ARRAYS FOR PRINTING. VBVL(3,MSUM)=STOIN*TLED VBVL(4,MSUM)=STOUT*TLED VBVL(1,MSUM)=VBVL(1,MSUM)+STOIN VBVL(2,MSUM)=VBVL(2,MSUM)+STOUT VBNM(MSUM)=TEXT(1) C C14-----INCREMENT BUDGET TERM COUNTER MSUM=MSUM+1 C C15-----UPDATE PRECONSOLIDATION HEAD ARRAY DO 310 KQ=1,NNDB K=LN(KQ) LOCT=(KQ-1)*NCR N=0 DO 300 IR=1,NROW DO 300 IC=1,NCOL N=N+1 LOC2=LOCT+N IF(IBOUND(IC,IR,K).LE.0) GO TO 300 HHNEW=HNEW(IC,IR,K) IF(HHNEW.LT.HC(LOC2)) HC(LOC2)=HHNEW 300 CONTINUE 310 CONTINUE ENDIF IF(NDF) THEN IF(IBD.EQ.1) THEN DO 390 K=1,NLAY DO 390 IR=1,NROW DO 390 IC=1,NCOL BUFF(IC,IR,K)=ZERO 390 CONTINUE ENDIF RATIN=ZERO RATOUT=ZERO C C16-----IF THIS IS A STEADY-STATE STRESS PERIOD, SKIP CALCULATIONS IF(ISSFLG(KPER).EQ.1) GO TO 421 C C17-----CALCULATE NO-DELAY INTERBED STORAGE CHANGE FOR EACH CELL LOC3=1-NN DO 420 KQ=1,NDB K=LDN(KQ) NQ=(KQ-1)*NCR STRGT=ZERO RATBSM=ZERO DO 410 IR=1,NROW NQR=NQ+(IR-1)*NCOL DO 410 IC=1,NCOL AREA=DELR(IC)*DELC(IR) LOC2=NQR+IC RNB2=RNB(LOC2) IF(RNB2.LT.1.0) GO TO 410 LOC3=LOC3+NN IF(IBOUND(IC,IR,K).LE.0) GO TO 410 HD1=DH(LOC3) HHNEW=HNEW(IC,IR,K) C C18-----CALCULATE CONDUCTANCE BETWEEN AQUIFER AND FIRST CELL IN INTERBED C18-----ACCOUNTING FOR BOTH HALVES OF RNB(LOC2) BEDS IN SYSTEM COND=4.*RNB2*DP(NZ(LOC2),1)*AREA/DZ(LOC2) C C19-----CALCULATE THE FLOW RATE INTO THE CELL RATB=DBLE(COND)*(DBLE(HD1)-HNEW(IC,IR,K)) STRGS=ZERO L1=LOC3 L2=LOC3+NN-1 DO 401 LOC4=L1,L2 HHOLD=DHP(LOC4) HHNEW=DH(LOC4) HHC=DHC(LOC4) C C20-----GET STORAGE CAPACITIES AT BEGINNING AND END OF TIME STEP. SBGN=DP(NZ(LOC2),2) SEND=SBGN IF(HHNEW.LT.HHC) SEND=DP(NZ(LOC2),3) C C21-----CALCULATE VOLUME CHANGE IN INTERBED STORAGE FOR TIME STEP. STR1=(HHC*(SEND-SBGN)+SBGN*HHOLD- 1 SEND*HHNEW)*DZ(LOC2)*RNB(LOC2)*2. IF(LOC4.EQ.L2) STR1=STR1*.5 STRGS=STRGS+STR1 C C22-----ACCUMULATE SUBSIDENCE ASSOCIATED WITH CHANGE IN STORAGE DCOM(LOC2)=DCOM(LOC2)+STR1 401 CONTINUE STRGS=STRGS*AREA STRGT=STRGT+STRGS RATS=STRGS*TLED IF(IBD.EQ.1) BUFF(IC,IR,K)=BUFF(IC,IR,K)+RATS RATBSM=RATBSM-RATS IF(RATS.LE.ZERO) THEN RATOUT=RATOUT-RATS ELSE RATIN=RATIN+RATS ENDIF 410 CONTINUE DVB(KQ,1)=DVB(KQ,1)+STRGT DVB(KQ,2)=DVB(KQ,2)+RATBSM*DELT DVB(KQ,3)=STRGT*TLED DVB(KQ,4)=RATBSM 420 CONTINUE C C23-----IF C-B-C FLOW TERMS WILL BE SAVED CALL UBUDSV TO RECORD THEM. 421 IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT(2),IIBSCB,BUFF, 1 NCOL,NROW,NLAY,IOUT) C C24-----MOVE RATES,VOLUMES & LABELS INTO ARRAYS FOR PRINTING. VBVL(3,MSUM)=RATIN VBVL(4,MSUM)=RATOUT VBVL(1,MSUM)=VBVL(1,MSUM)+RATIN*DELT VBVL(2,MSUM)=VBVL(2,MSUM)+RATOUT*DELT VBNM(MSUM)=TEXT(2) C C25-----INCREMENT BUDGET TERM COUNTER MSUM=MSUM+1 C C26-----UPDATE PRECONSOLIDATION HEAD ARRAY, PREVIOUS HEAD ARRAY FOR C26-----SYSTEMS OF DELAY INTERBEDS. DO 500 N=1,ND2 IF(DH(N).LT.DHC(N)) DHC(N)=DH(N) DHP(N)=DH(N) 500 CONTINUE ENDIF C C27----RETURN RETURN END SUBROUTINE GWF2SUB7OT(KSTP,KPER,IN,IGRID) C ****************************************************************** C PRINT AND STORE SUBSIDENCE, COMPACTION AND CRITICAL HEAD. C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ USE GLOBAL, ONLY: IOUT,NCOL,NROW,NLAY,NSTP,BUFF,ISSFLG USE GWFBASMODULE, ONLY:PERTIM,TOTIM USE GWFSUBMODULE ,ONLY: LN,LDN,SUB,HC,RNB,DCOM,DHC,DVB,NDF,NNDF, & NTSSUM,OCFLGS,OCLAY,ILSYS,ISBOCF,ISBOCU, & NN,NNDB,NDB CHARACTER*16 TEXT(7) LOGICAL IBDPR DATA TEXT & /' SUBSIDENCE', 'LAYER COMPACTION', & 'NDSYS COMPACTION', ' DSYS COMPACTION', & ' Z DISPLACEMENT', 'ND CRITICAL HEAD', 2 ' D CRITICAL HEAD'/ C ------------------------------------------------------------------ CALL SGWF2SUB7PNT(IGRID) C ZERO=0.0 NCR=NCOL*NROW C1------INITIALIZE TIME STEP POINTER TO RETRIEVE FLAGS FOR PRINTING AND C1------SAVING ARRAYS. SET FLAG FOR PRINTING BUDGET FOR DELAY INTERBEDS. NNSTP=NTSSUM(KPER)+KSTP IBDPR=.FALSE. IF(KSTP.EQ.NSTP(KPER).OR.OCFLGS(13,NNSTP)) IBDPR=.TRUE. IF(ISSFLG(KPER).EQ.1) IBDPR=.FALSE. C C3------PRINT AND STORE SUBSIDENCE, FIRST, CLEAR OUT BUFF. IF(OCFLGS(1,NNSTP).OR.OCFLGS(2,NNSTP)) THEN DO 30 K=1,NLAY DO 30 IR=1,NROW DO 30 IC=1,NCOL BUFF(IC,IR,K)=ZERO 30 CONTINUE C C4-------SUM COMPACTION IN ALL LAYERS TO GET SUBSIDENCE. IF(NNDF) THEN DO 50 KQ=1,NNDB LOCT=(KQ-1)*NCR N=0 DO 40 IR=1,NROW DO 40 IC=1,NCOL N=N+1 LOC2=LOCT+N BUFF(IC,IR,1)=BUFF(IC,IR,1)+SUB(LOC2) 40 CONTINUE 50 CONTINUE ENDIF IF(NDF) THEN DO 70 KQ=1,NDB LOCT=(KQ-1)*NCR N=0 DO 60 IR=1,NROW DO 60 IC=1,NCOL N=N+1 LOC2=LOCT+N BUFF(IC,IR,1)=BUFF(IC,IR,1)+DCOM(LOC2) 60 CONTINUE 70 CONTINUE ENDIF C C5-------PRINT SUBSIDENCE. IF(OCFLGS(1,NNSTP)) THEN IF(ISBOCF(1).LT.0) CALL ULAPRS(BUFF(:,:,1),TEXT(1),KSTP,KPER, 1 NCOL,NROW,1,-ISBOCF(1),IOUT) IF(ISBOCF(1).GE.0) CALL ULAPRW(BUFF(:,:,1),TEXT(1),KSTP,KPER, 1 NCOL,NROW,1,ISBOCF(1),IOUT) ENDIF C C6-------STORE SUBSIDENCE. IF(OCFLGS(2,NNSTP)) THEN CALL ULASAV(BUFF(:,:,1),TEXT(1),KSTP,KPER,PERTIM,TOTIM,NCOL, 1 NROW,1,ISBOCU(1)) ENDIF ENDIF C C7------PRINT AND STORE COMPACTION FOR EACH SYSTEM OF INTERBEDS, C7------INCLUDING DELAY AND NO-DELAY INTERBEDS. IF(OCFLGS(5,NNSTP).OR.OCFLGS(6,NNSTP)) THEN IF(NNDF) THEN DO 80 KQ=1,NNDB K=LN(KQ) LOC2=(KQ-1)*NCR+1 IF(OCFLGS(5,NNSTP)) THEN WRITE(IOUT,76) KQ 76 FORMAT(/,1X,' SYSTEM',I4,' OF NO-DELAY BEDS:') NEND=LOC2+NCOL*NROW-1 IF(ISBOCF(3).LT.0) CALL ULAPRS(SUB(LOC2:NEND),TEXT(3),KSTP, 1 KPER,NCOL,NROW,K,-ISBOCF(3),IOUT) IF(ISBOCF(3).GE.0) CALL ULAPRW(SUB(LOC2:NEND),TEXT(3),KSTP, 1 KPER,NCOL,NROW,K,ISBOCF(3),IOUT) ENDIF IF(OCFLGS(6,NNSTP)) THEN CALL ULASAV(SUB(LOC2:NEND),TEXT(3),KSTP,KPER,PERTIM,TOTIM, 1 NCOL,NROW,KQ,ISBOCU(3)) ENDIF 80 CONTINUE ENDIF IF(NDF) THEN DO 90 KQ=1,NDB K=LDN(KQ) LOC2=(KQ-1)*NCR+1 IF(OCFLGS(5,NNSTP)) THEN WRITE(IOUT,82) KQ 82 FORMAT(/,1X,' SYSTEM',I4,' OF DELAY BEDS:') NEND=LOC2+NCOL*NROW-1 IF(ISBOCF(3).LT.0) CALL ULAPRS(DCOM(LOC2:NEND),TEXT(4),KSTP, 1 KPER,NCOL,NROW,K,-ISBOCF(3),IOUT) IF(ISBOCF(3).GE.0) CALL ULAPRW(DCOM(LOC2:NEND),TEXT(4),KSTP, 1 KPER,NCOL,NROW,K,ISBOCF(3),IOUT) ENDIF IF(OCFLGS(6,NNSTP)) THEN CALL ULASAV(DCOM(LOC2:NEND),TEXT(4),KSTP,KPER,PERTIM,TOTIM, 1 NCOL,NROW,KQ,ISBOCU(3)) ENDIF 90 CONTINUE ENDIF ENDIF C C8------SUM COMPACTION IN EACH LAYER IN THE BUFF ARRAY FOR SAVING C8------OR PRINTING COMPACTION OR VERTICAL DISPLACEMENT BY MODEL C8------LAYER. FIRST, CLEAR OUT BUFF. IF(OCFLGS(3,NNSTP).OR.OCFLGS(4,NNSTP).OR. & OCFLGS(7,NNSTP).OR.OCFLGS(8,NNSTP)) THEN DO 93 NL=1,NLAY OCLAY(NL)=.FALSE. 93 CONTINUE DO 95 K=1,NLAY DO 95 IR=1,NROW DO 95 IC=1,NCOL BUFF(IC,IR,K)=ZERO 95 CONTINUE C C9-------SUM NO-DELAY COMPACTION IN ALL MODEL LAYERS. IF(NNDF) THEN DO 105 KQ=1,NNDB K=LN(KQ) OCLAY(K)=.TRUE. LOCT=(KQ-1)*NCR N=0 DO 100 IR=1,NROW DO 100 IC=1,NCOL N=N+1 LOC2=LOCT+N BUFF(IC,IR,K)=BUFF(IC,IR,K)+SUB(LOC2) 100 CONTINUE 105 CONTINUE ENDIF C C10------SUM DELAY COMPACTION IN ALL MODEL LAYERS. IF(NDF) THEN DO 115 KQ=1,NDB K=LDN(KQ) OCLAY(K)=.TRUE. LOCT=(KQ-1)*NCR N=0 DO 110 IR=1,NROW DO 110 IC=1,NCOL N=N+1 LOC2=LOCT+N BUFF(IC,IR,K)=BUFF(IC,IR,K)+DCOM(LOC2) 110 CONTINUE 115 CONTINUE ENDIF ENDIF C C11-----PRINT COMPACTION BY LAYER. IF(OCFLGS(3,NNSTP)) THEN DO 120 KL=1,NLAY IF(.NOT.OCLAY(KL)) GO TO 120 KKL=KL IF(ISBOCF(2).LT.0) CALL ULAPRS(BUFF(:,:,KL),TEXT(2),KSTP,KPER, 1 NCOL,NROW,KKL,-ISBOCF(2),IOUT) IF(ISBOCF(2).GE.0) CALL ULAPRW(BUFF(:,:,KL),TEXT(2),KSTP,KPER, 1 NCOL,NROW,KKL,ISBOCF(2),IOUT) 120 CONTINUE ENDIF C C12-----STORE COMPACTION BY LAYER. IF(OCFLGS(4,NNSTP)) THEN DO 125 KL=1,NLAY IF(.NOT.OCLAY(KL)) GO TO 125 KKL=KL CALL ULASAV(BUFF(:,:,KL),TEXT(2),KSTP,KPER,PERTIM,TOTIM,NCOL, 1 NROW,KKL,ISBOCU(2)) 125 CONTINUE ENDIF C C13----CALCULATE VERTICAL DISPLACEMENT. IF(OCFLGS(7,NNSTP).OR.OCFLGS(8,NNSTP)) THEN NL1=NLAY-1 IF(NLAY.GT.1) THEN DO 140 KL=NL1,1,-1 LOCT3=(KL-1)*NCR N=0 DO 135 IR=1,NROW DO 135 IC=1,NCOL N=N+1 BUFF(IC,IR,KL)=BUFF(IC,IR,KL)+BUFF(IC,IR,KL+1) 135 CONTINUE 140 CONTINUE ENDIF C C14-----PRINT VERTICAL DISPLACEMENT FOR ALL MODEL LAYERS. IF(OCFLGS(7,NNSTP)) THEN DO 145 KL=1,NLAY KKL=KL IF(ISBOCF(4).LT.0) CALL ULAPRS(BUFF(:,:,KL),TEXT(5),KSTP,KPER, 1 NCOL,NROW,KKL,-ISBOCF(4),IOUT) IF(ISBOCF(4).GE.0) CALL ULAPRW(BUFF(:,:,KL),TEXT(5),KSTP,KPER, 1 NCOL,NROW,KKL,ISBOCF(4),IOUT) 145 CONTINUE ENDIF C C15-----SAVE VERTICAL DISPLACEMENT FOR ALL MODEL LAYERS. IF(OCFLGS(8,NNSTP)) THEN DO 150 KL=1,NLAY KKL=KL CALL ULASAV(BUFF(:,:,KL),TEXT(5),KSTP,KPER,PERTIM,TOTIM,NCOL, 1 NROW,KKL,ISBOCU(4)) 150 CONTINUE ENDIF ENDIF C C16-----PRINT CRITICAL HEAD FOR SYSTEMS OF NO-DELAY INTERBEDS. C16-----STORAGE. IF(OCFLGS(9,NNSTP).OR.OCFLGS(10,NNSTP)) THEN IF(NNDF) THEN DO 155 NL=1,NLAY OCLAY(NL)=.TRUE. 155 CONTINUE DO 160 KQ=1,NNDB K=LN(KQ) LOC2=(KQ-1)*NCR+1 IF(.NOT.OCLAY(K)) GO TO 160 NSYS=1 ILSYS(1)=KQ IF(KQ.LT.NNDB) THEN KS1=KQ+1 DO 152 KS=KS1,NNDB IF(LN(KS).EQ.K) THEN NSYS=NSYS+1 ILSYS(NSYS)=KS ENDIF 152 CONTINUE ENDIF OCLAY(K)=.FALSE. IF(OCFLGS(9,NNSTP)) THEN WRITE(IOUT,154) (ILSYS(NS),NS=1,NSYS) 154 FORMAT(/,1X,' SYSTEM OR SYSTEMS OF NO-DELAY BEDS:',20I3) NEND=LOC2+NCOL*NROW-1 IF(ISBOCF(5).LT.0) CALL ULAPRS(HC(LOC2:NEND),TEXT(6),KSTP,KPER, 1 NCOL,NROW,K,-ISBOCF(5),IOUT) IF(ISBOCF(5).GE.0) CALL ULAPRW(HC(LOC2:NEND),TEXT(6),KSTP,KPER, 1 NCOL,NROW,K,ISBOCF(5),IOUT) ENDIF IF(OCFLGS(10,NNSTP)) THEN CALL ULASAV(HC(LOC2:NEND),TEXT(6),KSTP,KPER,PERTIM,TOTIM,NCOL, 1 NROW,K,ISBOCU(5)) ENDIF 160 CONTINUE ENDIF ENDIF C C17-----PRINT CRITICAL HEAD FOR ALL SYSTEMS OF DELAY INTERBED. IF(OCFLGS(11,NNSTP).OR.OCFLGS(12,NNSTP)) THEN IF(NDF) THEN LOC4=0 DO 190 KQ=1,NDB K=LDN(KQ) LOCT=(KQ-1)*NCR N=0 DO 180 IR=1,NROW DO 180 IC=1,NCOL N=N+1 BUFF(IC,IR,1)=ZERO LOC2=LOCT+N IF(RNB(LOC2).LT.1.0) GO TO 180 LOC4=LOC4+NN BUFF(IC,IR,1)=DHC(LOC4) 180 CONTINUE IF(OCFLGS(11,NNSTP)) THEN WRITE(IOUT,82) KQ IF(ISBOCF(6).LT.0) CALL ULAPRS(BUFF(:,:,1),TEXT(7),KSTP,KPER, 1 NCOL,NROW,K,-ISBOCF(6),IOUT) IF(ISBOCF(6).GE.0) CALL ULAPRW(BUFF(:,:,1),TEXT(7),KSTP,KPER, 1 NCOL,NROW,K,ISBOCF(6),IOUT) ENDIF IF(OCFLGS(12,NNSTP)) THEN CALL ULASAV(BUFF(:,:,1),TEXT(7),KSTP,KPER,PERTIM,TOTIM,NCOL, 1 NROW,KQ,ISBOCU(6)) ENDIF 190 CONTINUE ENDIF ENDIF C C18-----PRINT VOLUMETRIC BUDGET FOR SYSTEMS OF DELAY INTERBEDS IF(NDF.AND.IBDPR) THEN WRITE(IOUT,230) KSTP,KPER SUMSV=ZERO SUMBV=ZERO SUMSR=ZERO SUMBR=ZERO DO 200 KQ=1,NDB DISCV=ZERO DISCR=ZERO SUMV=DVB(KQ,1)+DVB(KQ,2) SUMR=DVB(KQ,3)+DVB(KQ,4) SUMSV=SUMSV+DVB(KQ,1) SUMBV=SUMBV+DVB(KQ,2) SUMSR=SUMSR+DVB(KQ,3) SUMBR=SUMBR+DVB(KQ,4) IF(DVB(KQ,1).NE.ZERO) DISCV=100.*SUMV/DVB(KQ,1) IF(DVB(KQ,3).NE.ZERO) DISCR=100.*SUMR/DVB(KQ,3) WRITE(IOUT,240) KQ,DVB(KQ,1),DVB(KQ,2),SUMV,DISCV, 1 DVB(KQ,3),DVB(KQ,4),SUMR,DISCR 200 CONTINUE IF(NDB.GT.1) THEN DISCV=ZERO DISCR=ZERO SUMV=SUMSV+SUMBV SUMR=SUMSR+SUMBR IF(SUMSV.NE.ZERO) DISCV=100.*SUMV/SUMSV IF(SUMSR.NE.ZERO) DISCR=100.*SUMR/SUMSR WRITE(IOUT,250) SUMSV,SUMBV,SUMV,DISCV, 1 SUMSR,SUMBR,SUMR,DISCR ENDIF ENDIF C C19-----RETURN RETURN C C20-----FORMATS C 230 FORMAT(/,31X,'VOLUMETRIC BUDGET FOR SYSTEMS OF DELAY INTERBEDS', 1 /,42X,'AT END OF TIME STEP',I3,' IN ','STRESS PERIOD',I3, 2 //,' | C U M U L A T I V E ', 3 'V O L U M E S L**3 | R A T E S F O R T H I S', 4 ' T I M E S T E P L**3/T |',/,' SYSTEM | CHANGE IN', 5 ' BOUNDARY PERCENT | CHANGE IN ', 6 ' BOUNDARY PERCENT |',/,' NUMBER | ', 7 ' STORAGE FLOW SUM DISCREPANCY | ', 8 ' STORAGE FLOW SUM DISCREPANCY |',/,2X, 9 8('-'),'|',56('-'),'|',56('-'),'|') 240 FORMAT(I7,4G15.5,4G15.5) 250 FORMAT(2X,8('-'),'|',56('-'),'|',56('-'),'|',/,' TOTALS:', 1 4G15.5,4G15.5) C END SUBROUTINE SGWF2SUB7A(HAQ,TLED,CI,SSE,SSV,DZ,DH,DHP,DHC,NN) C ****************************************************************** C ASSEMBLE COEFFICIENTS FOR SOLVING FOR HEAD DISTRIBUTION C IN ONE STRING OF CELLS REPRESENTING ONE-HALF OF A DOUBLY C DRAINING INTERBED C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ USE GWFSUBMODULE ,ONLY:A1,A2,BB DIMENSION DH(NN),DHP(NN),DHC(NN) C ------------------------------------------------------------------ C C1------INITIALIZE C CI2=CI*2.0 DS=DZ*TLED NN1=NN-1 C2------SET COEFFICIENTS FOR CELL BORDERING AQUIFER SS=SSE A2(1)=CI HD=DH(1) HC=DHC(1) IF(HD.LT.HC) SS=SSV A1(1)=-3.*CI-DS*SS BB(1)=DS*(SSE*(HC-DHP(1))-HC*SS)-CI2*HAQ-A1(1)*HD C3------SET COEFFICIENTS FOR INTERIOR CELLS BB(2)=-CI*HD DO 10 N=2,NN1 SS=SSE A2(N)=CI HD=DH(N) HC=DHC(N) IF(HD.LT.HC) SS=SSV CHN=-CI*HD BB(N-1)=BB(N-1)+CHN BB(N+1)=CHN A1(N)=-CI2-DS*SS BB(N)=BB(N)+DS*(SSE*(HC-DHP(N))-HC*SS)-A1(N)*HD 10 CONTINUE C4------SET COEFFICIENTS FOR CELL BORDERING MIDPLANE OF INTERBED SS=SSE A2(NN)=CI HD=DH(NN) HC=DHC(NN) BB(NN1)=BB(NN1)-CI*HD IF(HD.LT.HC) SS=SSV A1(NN)=-CI-0.5*DS*SS BB(NN)=BB(NN)+DS*0.5*(SSE*(HC-DHP(NN))-HC*SS)-A1(NN)*HD C5------RETURN RETURN END SUBROUTINE SGWF2SUB7S(NN) C ****************************************************************** C SOLVE SYSTEM OF EQUATIONS WITH A SYMMETRICAL TRI-DIAGONAL C COEFFICIENT MATRIX C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ USE GWFSUBMODULE ,ONLY: A1,A2,BB C ------------------------------------------------------------------ C C1------TRIANGULARIZE LEFT-HAND SIDE MATRIX NN1=NN-1 DO 30 N=1,NN1 F=1./A1(N) C=A2(N)*F I=N+1 A1(I)=A1(I)-C*A2(N) A2(N)=C BB(I)=BB(I)-C*BB(N) BB(N)=BB(N)*F 30 CONTINUE BB(NN)=BB(NN)/A1(NN) C C2------BACK SUBSTITE FOR SOLUTION DO 40 N=NN1,1,-1 BB(N)=BB(N)-A2(N)*BB(N+1) 40 CONTINUE RETURN END SUBROUTINE GWF2SUB72D1D(BUFF,NCOL,NROW,D,ND,LOC) C ****************************************************************** C Move 2-D array into 1-D array C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ DIMENSION BUFF(NCOL,NROW),D(ND) C ------------------------------------------------------------------ L=LOC-1 DO 10 I=1,NROW DO 10 J=1,NCOL L=L+1 D(L)=BUFF(J,I) 10 CONTINUE RETURN END SUBROUTINE GWF2SUB7SV(IGRID) C ****************************************************************** C SAVE INTERBED STORAGE DATA FOR FUTURE RESTART C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ USE GWFSUBMODULE ,ONLY: DH,DHC,NDF,ND2,IDSAVE C ------------------------------------------------------------------ CALL SGWF2SUB7PNT(IGRID) C C1-----PROCESS IF SAVE OPTION SELECTED AND DELAY INTERBEDS EXIST IF(IDSAVE.GT.0.AND.NDF) THEN C2-----WRITE OUT NUMBER OF NODAL HEAD AND PRECONSOLIDATION HEAD VALUES C2-----THAT ARE BEING SAVED TO DISK WRITE(IDSAVE) ND2 C3-----WRITE ARRAYS FOR EACH SYSTEM OF INTERBEDS WRITE(IDSAVE) (DH(N),N=1,ND2) WRITE(IDSAVE) (DHC(N),N=1,ND2) ENDIF C4-----RETURN RETURN END SUBROUTINE GWF2SUB7DA(IGRID) C C ****************************************************************** C DEALLOCATE DYNAMIC STORAGE FOR SUB PACKAGE C ****************************************************************** C C SPECIFICATIONS: C ------------------------------------------------------------------ USE GWFSUBMODULE C ------------------------------------------------------------------ C DEALLOCATE (GWFSUBDAT(IGRID)%IIBSCB) DEALLOCATE (GWFSUBDAT(IGRID)%ITMIN) DEALLOCATE (GWFSUBDAT(IGRID)%NNDB) DEALLOCATE (GWFSUBDAT(IGRID)%NDB) DEALLOCATE (GWFSUBDAT(IGRID)%NMZ) DEALLOCATE (GWFSUBDAT(IGRID)%NN) DEALLOCATE (GWFSUBDAT(IGRID)%ND2) DEALLOCATE (GWFSUBDAT(IGRID)%IDSAVE) DEALLOCATE (GWFSUBDAT(IGRID)%AC1) DEALLOCATE (GWFSUBDAT(IGRID)%AC2) DEALLOCATE (GWFSUBDAT(IGRID)%NDF) DEALLOCATE (GWFSUBDAT(IGRID)%NNDF) DEALLOCATE (GWFSUBDAT(IGRID)%ISBOCF) DEALLOCATE (GWFSUBDAT(IGRID)%ISBOCU) DEALLOCATE (GWFSUBDAT(IGRID)%OCFLGS) DEALLOCATE (GWFSUBDAT(IGRID)%OCLAY) DEALLOCATE (GWFSUBDAT(IGRID)%ILSYS) DEALLOCATE (GWFSUBDAT(IGRID)%NTSSUM) DEALLOCATE (GWFSUBDAT(IGRID)%LN) DEALLOCATE (GWFSUBDAT(IGRID)%LDN) DEALLOCATE (GWFSUBDAT(IGRID)%NZ) DEALLOCATE (GWFSUBDAT(IGRID)%RNB) DEALLOCATE (GWFSUBDAT(IGRID)%DH) DEALLOCATE (GWFSUBDAT(IGRID)%DHP) DEALLOCATE (GWFSUBDAT(IGRID)%DHC) DEALLOCATE (GWFSUBDAT(IGRID)%DZ) DEALLOCATE (GWFSUBDAT(IGRID)%HC) DEALLOCATE (GWFSUBDAT(IGRID)%SCE) DEALLOCATE (GWFSUBDAT(IGRID)%SCV) DEALLOCATE (GWFSUBDAT(IGRID)%DCOM) DEALLOCATE (GWFSUBDAT(IGRID)%A1) DEALLOCATE (GWFSUBDAT(IGRID)%A2) DEALLOCATE (GWFSUBDAT(IGRID)%BB) DEALLOCATE (GWFSUBDAT(IGRID)%SUB) DEALLOCATE (GWFSUBDAT(IGRID)%DP) DEALLOCATE (GWFSUBDAT(IGRID)%DVB) C2-----RETURN RETURN END SUBROUTINE SGWF2SUB7PNT(IGRID) C Change SUB data to a different grid. USE GWFSUBMODULE C IIBSCB=>GWFSUBDAT(IGRID)%IIBSCB ITMIN=>GWFSUBDAT(IGRID)%ITMIN NNDB=>GWFSUBDAT(IGRID)%NNDB NDB=>GWFSUBDAT(IGRID)%NDB NMZ=>GWFSUBDAT(IGRID)%NMZ NN=>GWFSUBDAT(IGRID)%NN ND2=>GWFSUBDAT(IGRID)%ND2 IDSAVE=>GWFSUBDAT(IGRID)%IDSAVE AC1=>GWFSUBDAT(IGRID)%AC1 AC2=>GWFSUBDAT(IGRID)%AC2 NDF=>GWFSUBDAT(IGRID)%NDF NNDF=>GWFSUBDAT(IGRID)%NNDF ISBOCF=>GWFSUBDAT(IGRID)%ISBOCF ISBOCU=>GWFSUBDAT(IGRID)%ISBOCU OCFLGS=>GWFSUBDAT(IGRID)%OCFLGS OCLAY=>GWFSUBDAT(IGRID)%OCLAY ILSYS=>GWFSUBDAT(IGRID)%ILSYS NTSSUM=>GWFSUBDAT(IGRID)%NTSSUM LN=>GWFSUBDAT(IGRID)%LN LDN=>GWFSUBDAT(IGRID)%LDN NZ=>GWFSUBDAT(IGRID)%NZ RNB=>GWFSUBDAT(IGRID)%RNB DH=>GWFSUBDAT(IGRID)%DH DHP=>GWFSUBDAT(IGRID)%DHP DHC=>GWFSUBDAT(IGRID)%DHC DZ=>GWFSUBDAT(IGRID)%DZ HC=>GWFSUBDAT(IGRID)%HC SCE=>GWFSUBDAT(IGRID)%SCE SCV=>GWFSUBDAT(IGRID)%SCV DCOM=>GWFSUBDAT(IGRID)%DCOM A1=>GWFSUBDAT(IGRID)%A1 A2=>GWFSUBDAT(IGRID)%A2 BB=>GWFSUBDAT(IGRID)%BB SUB=>GWFSUBDAT(IGRID)%SUB DP=>GWFSUBDAT(IGRID)%DP DVB=>GWFSUBDAT(IGRID)%DVB C RETURN END SUBROUTINE SGWF2SUB7PSV(IGRID) C Save SUB data for a grid. USE GWFSUBMODULE C GWFSUBDAT(IGRID)%IIBSCB=>IIBSCB GWFSUBDAT(IGRID)%ITMIN=>ITMIN GWFSUBDAT(IGRID)%NNDB=>NNDB GWFSUBDAT(IGRID)%NDB=>NDB GWFSUBDAT(IGRID)%NMZ=>NMZ GWFSUBDAT(IGRID)%NN=>NN GWFSUBDAT(IGRID)%ND2=>ND2 GWFSUBDAT(IGRID)%IDSAVE=>IDSAVE GWFSUBDAT(IGRID)%AC1=>AC1 GWFSUBDAT(IGRID)%AC2=>AC2 GWFSUBDAT(IGRID)%NDF=>NDF GWFSUBDAT(IGRID)%NNDF=>NNDF GWFSUBDAT(IGRID)%ISBOCF=>ISBOCF GWFSUBDAT(IGRID)%ISBOCU=>ISBOCU GWFSUBDAT(IGRID)%OCFLGS=>OCFLGS GWFSUBDAT(IGRID)%OCLAY=>OCLAY GWFSUBDAT(IGRID)%ILSYS=>ILSYS GWFSUBDAT(IGRID)%NTSSUM=>NTSSUM GWFSUBDAT(IGRID)%LN=>LN GWFSUBDAT(IGRID)%LDN=>LDN GWFSUBDAT(IGRID)%NZ=>NZ GWFSUBDAT(IGRID)%RNB=>RNB GWFSUBDAT(IGRID)%DH=>DH GWFSUBDAT(IGRID)%DHP=>DHP GWFSUBDAT(IGRID)%DHC=>DHC GWFSUBDAT(IGRID)%DZ=>DZ GWFSUBDAT(IGRID)%HC=>HC GWFSUBDAT(IGRID)%SCE=>SCE GWFSUBDAT(IGRID)%SCV=>SCV GWFSUBDAT(IGRID)%DCOM=>DCOM GWFSUBDAT(IGRID)%A1=>A1 GWFSUBDAT(IGRID)%A2=>A2 GWFSUBDAT(IGRID)%BB=>BB GWFSUBDAT(IGRID)%SUB=>SUB GWFSUBDAT(IGRID)%DP=>DP GWFSUBDAT(IGRID)%DVB=>DVB C RETURN END