C     -----------------
      SUBROUTINE WFLUX1
C     -----------------
C
C     SUBROUTINE AUTHOR:
C     TIM DURBIN
C     TIMOTHY J. DURBIN, INC.
C     5330 PRIMROSE DRIVE, SUITE 228
C     FAIR OAKS, CA 95628
C     (916) 536-2314
C
C     PURPOSE:
C     THIS SUBROUTINE IS USED TO REPRESENT SPECIFIED-FLUX BOUNDARY
C     CONDITIONS AND INTERNAL SOURCES AND SINKS, WHICH CAN VARY WITH
C     TIME AS SPECIFIED BY THE USER IN A TABLE. THE FIRST ENTRY POINT IS
C     CALLED ONCE BYTHE MAIN PROGRAM. THE SECOND ENTRY POINT IS CALLED
C     ONCE FOR EACH TIMESTEP ITERATION. THE THIRD ENTRY POINT IS CALLED
C     ONCE FOR EACH TIMESTEP.
C
      INCLUDE 'DIMEN.INC'
      INCLUDE 'WUNITS.INC'
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
C
C     INTERNAL SCALARS AND ARRAYS.
C     INFILE      NAME OF INPUT FILE.
C     OUTFILE     OUTPUT FILE NAME.
C     OUTEXT      EXTENSION FOR OUTPUT FILE.
C     XCNAME      COMPARTMENT NAME.
C     WNAME()     WELL NAME.
C     FXNAME()    SPECIFIED-FLUX NAME ASSIGNED TO NODE.
C     AE()        ELEMENT CONDUCTANCE MATRIX.
C     FXTAB()     SPECIFIED-FLUX TABLE ASSIGNED TO NODE.
C     NPT()       NUMBER OF POINTS IN SPECIFIED-FLUX TABLE.
C     IGET()      SWITCH FOR INTERPOLATION METHOD.
C     FXRATE()    STEADY-STATE FLUX RATE AT NODE.
C     LEAK()      LEAKANCE OF WELL LINK.
C     WLRATE()    WELL FLUX RATE.
C     WLTAB()     FLUX RATE TABLE.
C     TTB()       ELAPSED TIME ASSOCIATED WITH POINT IN TABLE.
C     QTB()       TRANSIENT-STATE SCALING FACTOR ASSOCIATED WITH POINT
C                 IN TABLE.
C     XTT()       TIME ENTRIES IN TABLE FOR PASSING TO TABLE FUNCTION.
C     XQT()       SPECIFIED-FLUX ENTRIES IN TABLE FOR PASSING TO TABLE
C                 FUNCTION.
C     H()         HEAD VALUE AT NODE.
C     IECHO1      SWITCH FOR ECHO OF SPECIFIED-FLUX INPUTS.
C     IECHO2      SWITCH FOR ECHO OF STEP FLUXES.
C     NFXNC       NUMBER OF SPECIFIED-FLUX NODES IN COMPARTMENT.
C     NTABC       NUMBER OF SPECIFIED-FLUX TABLES IN COMPARTMENT.
C     NWELC       NUMBER OF WELLS.
C     IN2()       WELL NODE.
C
C     INTERNAL ARRAYS.
      CHARACTER*80 INFILE,OUTFILE,OUTEXT
      CHARACTER*20 XCNAME,WNAME(MAXWEL)
      CHARACTER*20 FXNAME(MAXFLX)
      REAL*8 AE(4,4)
      INTEGER*4 FXTAB(MAXFLX),NPT(MAXTBF),IGET(MAXTBF)
      REAL*8 FXRATE(MAXFLX)
      REAL*8 LEAK(MAXWEL,MAXWND),WRATE(MAXWEL),HMIN(MAXWEL)
      REAL*8 WLRATE(MAXWEL)
      INTEGER*4 WLTAB(MAXWEL)
      REAL*8 TTB(MAXTBF,MAXPTF),QTB(MAXTBF,MAXPTF)
      REAL*8 XTT(MAXPTF),XQT(MAXPTF)
      REAL*8 H(MAXNN)
C
C     COMMON BLOCKS.
C     FNAME()     NAMES OF INPUT FILES.
C     RNAME1      SIMULATION NAME FOR GROUNDWATER-FLOW FILES.
C     RNAME2      SIMULATION NAME FOR SOLUTE-TRANSPORT FILES.
C     NCELL       NUMBER OF COMPARTMENTS.
C     CNAME()     NAMES FOR COMPARTMENTS.
C     NNC()       NUMBER OF NODES IN A COMPARTMENT.
C     NEC()       NUMBER OF WEDGE ELEMENTS WITHIN A COMPARTMENT.
C     IN()        ELEMENT INCIDENCE FOR A WEDGE ELEMENT.
C     IN2()       ELEMENT INCIDENCES FOR A TETRAHEDRAL ELEMENT,
C                     WELL LINK, FAULT ELEMENT, OR FAULT LINK.
C     IEL()       PARENT WEDGE ELEMENT FOR TETRAHEDRAL ELEMENT.
C     X()         X-COORDINATE OF A NODE.
C     Y()         Y-COORDINATE OF A NODE.
C     Z()         Z-COORDINATE OF A NODE.
C     IA()        VECTOR INDICATING POSITIONS OF ROW ENDS IN
C                     VECTOR A() OR B().
C     JA()        VECTOR INDICATION COLUMN NUMBER ASSOCIATED WITH
C                     ENTRIES IN VECTORS A() OR B().
C     A()         VECTOR OF NON-ZERO COEFFICIENTS FOR CONDUCTANCE
C                     MATRIX.
C     B()         VECTOR OF NON-ZERO COEFFICIENTS FOR CAPACITANCE
C                     MATRIX.
C     F()         RIGHTHAND VECTOR FOR SYSTEM OF EQUATIONS.
C     QX()        DARCY GROUNDWATER VELOCITY IN THE X-DIRECTION.
C     QY()        DARCY GROUNDWATER VELOCITY IN THE Y-DIRECTION.
C     QZ()        DARCY GROUNDWATER VELOCITY IN THE Z-DIRECTION.
C     NE          GLOBAL NUMBER OF WEDGE ELEMENTS.
C     NE2         GLOBAL NUMBER OF TETRAHEDRAL ELEMENTS.
C     NE3         GLOBAL NUMBER OF TETRAHEDRAL ELEMENS AND WELL LINKS.
C     NE4         GLOBAL NUMBER OF TETHEDRALAL ELEMENTS, WELL LINKS,
C                     AND FAULT ELEMENTS.
C     NE5         GLOBAL NUMBER OF TETHEDRALAL ELEMENTS, WELL LINKS,
C                     FAULT ELEMENTS, AND FAULT LINKS.
C     NN          GLOBAL NUMBER OF MESH NODES.
C     NN3         GLOBAL NUMBER OF MESH NODES AND WELL NODES.
C     NN5         GLOBAL NUMBER OF MESH NODES, WELL NODES, AND
C                     FAULT NODES.
C     ISTEP       TIME STEP NUMBER WITHIN FLOW SIMULATION.
C     ISTEP1      STARTING STEP NUMBER FOR SIMULATION.
C     ISTEP2      ENDING STEP FOR SIMULATION.
C     DELT        LENGTH OF TIME STEP.
C     STIME       INITIAL ELAPSED TIME FOR SIMULATION.
C     IPRN        FLAG TO INDICATED PRINT OUTPUT FOR TIME STEP.
C     IPLT        FLAG TO INDICATE PLOT OUTPUT FOR TIME STEP.
C     NTSTEP      NUMBER OF TRANSPORT STEPS WITHIN FLOW STEP.
C     ITSTEP      TIME STEP NUMBER WITHIN TRANSPORT SIMULATION.
C     NIA         NUMBER OF ROWS REPRESENTED IN VECTORS A() OR B().
C     NJA         NUMBER OF NON-ZERO COEFFICIENTS IN VECTORS A() OR B().
C     ITRAN       SWITCH FOR TRANSPORT SIMULATION.
C     QFLX        FLUX RATE.
C     FXNODE      NODE NUMBER FOR SPECIFIED FLUX.
C     QWEL        WELL FLUX.
C     QLNK        LINK FLUX.
C     NWLN        NUMBER OF WELL NODES.
C     IRFLX       INTERNAL READ FOR FLX UNIT NUMBER.
C     IWFLX       INTERNAL WRITE FOR FLX UNIT NUMBER.
C
C     COMMON BLOCKS.
      CHARACTER*80 FNAME(30)
      CHARACTER*12 RNAME1,RNAME2
      CHARACTER*20 CNAME(MAXCEL)
      INTEGER*4 NNC(MAXCEL),NEC(MAXCEL)
      INTEGER*4 IN(MAXNE,6),IN2(MAXNE2,4),IEL(MAXNE2)
      REAL*8 X(MAXNN),Y(MAXNN),Z(MAXNN)
      INTEGER*4 IA(MAXNN+1),JA(MAXNON)
      REAL*8 A(MAXNON),B(MAXNON),F(MAXNN)
      REAL*8 QX(MAXNE2),QY(MAXNE2),QZ(MAXNE2)
      INTEGER*4 FXNODE(MAXFLX),NWLN(MAXWEL)
      REAL*8 QFLX(MAXFLX),QLNK(MAXWEL,MAXWND),QWEL(MAXWEL)
      COMMON /NAMES/ FNAME,RNAME1,RNAME2
      COMMON /CELLS/ NCELL,CNAME,NNC,NEC
      COMMON /MESH/ X,Y,Z,IN,IN2,IEL,NE,NE2,NE3,NE4,NE5,NN,NN3,NN5
      COMMON /STEPS/ ISTEP,ISTEP1,ISTEP2,DELT,STIME,IPRN,IPLT,NTSTEP,
     +    ITSTEP
      COMMON /MATRIX/ A,B,IA,JA,NIA,NJA,F
      COMMON /SHARE/ QX,QY,QZ,ITRAN
      COMMON /FLUX/ QFLX,FXNODE,QWEL,QLNK,NWLN
      SAVE
C
C     OPEN INPUT AND OUTPUT FILES.
      WRITE(*,902)
      OUTEXT='.FLX.OUT'
      INFILE=FNAME(8)
      CALL OPENFILES(RNAME1,INFILE,OUTFILE,OUTEXT,IRFLX,IWFLX)
      WRITE(IWFLX,904)
      WRITE(IWFLX,906) INFILE,OUTFILE
C
C     LOOP OVER COMPARTMENTS.
      IFXN=0
      ITAB=0
      IWEL=0
      NFXN=0
      NTAB=0
      NWEL=0
      L2=NE2
      ILAST=0
      JLAST=0
      NE3=NE2
      READ(IRFLX,910) IECHO2
      CELL_LOOP:DO ICELL=1,NCELL
C
C         INPUT DATASET PARAMETERS AND CHECK PARAMETERS.
          READ(IRFLX,908) XCNAME
          READ(IRFLX,910) NFXNC,NTABC,NWELC,IECHO1
          WRITE(IWFLX,912) NCELL,ICELL,XCNAME,NFXNC,NTABC,NWELC,
     +        IECHO1,IECHO2
          NFXN=NFXN+NFXNC
          NTAB=NTAB+NTABC
          NWEL=NWEL+NWELC
C
C         CHECK COMPARTMENT IN LIST
          IF (XCNAME.NE.CNAME(ICELL)) THEN
              WRITE(*,914)
              WRITE(*,916)
              STOP
          ENDIF
C
C         CHECK NUMBER OF NODES.
          IF (NFXN.GT.MAXFLX) THEN
              WRITE(*,918)
              WRITE(*,916)
              STOP
          ENDIF
C
C         CHECK NUMBER OF TABLES.
          IF (NTAB.GT.MAXTBF) THEN
              WRITE(*,920)
              WRITE(*,916)
              STOP
          ENDIF
C
C         INPUT LOCAL PARAMETERS FOR SPECIFIED-FLUX NODES.
          IF (IECHO1.NE.0 .AND. NFXNC.NE.0) WRITE(IWFLX,922)
          NODE_LOOP1:DO II=1,NFXNC
C
C             INPUT NODAL PARAMETERS.
              IFXN=IFXN+1
              READ(IRFLX,924) FXNODE(IFXN),FXRATE(IFXN),FXTAB(IFXN),
     +            FXNAME(IFXN)
              IF(FXTAB(IFXN).NE.0) THEN
                  ITAB=FXTAB(IFXN)+JLAST
              ELSE
                  ITAB=0
              ENDIF
              IF (IECHO1.EQ.1) WRITE(IWFLX,926) FXNODE(IFXN),
     +            FXNODE(IFXN)+ILAST,FXRATE(IFXN),FXTAB(IFXN),ITAB,
     +            FXNAME(IFXN)
              FXNODE(IFXN)=FXNODE(IFXN)+ILAST
              FXTAB(IFXN)=ITAB
C
C             CHECK DIMENSIONING FOR NODAL NUMBERS.
              IF (FXNODE(IFXN).GT.MAXNN) THEN
                  WRITE(*,928)
                  WRITE(*,916)
                  STOP
              ENDIF
C
C             CHECK TABLE NUMBER AGAINST NUMBER OF TABLES.
              IF (FXTAB(IFXN).GT.NTAB) THEN
                  WRITE(*,930)
                  WRITE(*,916)
                  STOP
              ENDIF
          END DO NODE_LOOP1
C
C         INPUT LOCAL FLUX TABLES.
          TABLE_LOOP1:DO J=1,NTABC
C
C             INPUT NUMBER OF POINTS IN LOCAL TABLE
              ITAB=J+JLAST
              READ(IRFLX,932) NPT(ITAB),IGET(ITAB)
              IF (IECHO1.NE.0) WRITE(IWFLX,934) J,ITAB,NPT(ITAB),
     +            IGET(ITAB)
C
C             CHECKING FOR VALID INTERPOLATION METHOD.
              IF (IGET(ITAB).LT.1 .OR. IGET(ITAB).GT.2) THEN
                  WRITE(*,935)
                  WRITE(*,916)
                  STOP
              ENDIF
C
C             CHECK DIMENSIONING FOR NUMBER OF POINTS IN TABLE.
              IF (NPT(ITAB).GT.MAXPTF) THEN
                  WRITE(*,936)
                  WRITE(*,916)
                  STOP
              ENDIF
C
C             INPUT TABLE ENTRIES.
              DO IPT=1,NPT(ITAB)
                  READ(IRFLX,938) TTB(ITAB,IPT),QTB(ITAB,IPT)
                  IF (IECHO1.NE.0) WRITE(IWFLX,940) IPT,
     +                TTB(ITAB,IPT),QTB(ITAB,IPT)
              END DO
          END DO TABLE_LOOP1
C
C         INPUT WELLS.
          IF (NWELC.NE.0 .AND. IECHO1.NE.0) WRITE(IWFLX,942) NWELC
          WELL_LOOP1:DO II=1,NWELC
C
C             GENERAL SPECIFICATIONS.
              IWEL=IWEL+1
              READ(IRFLX,943) NWLN(IWEL),WLRATE(IWEL),WLTAB(IWEL),
     +            WNAME(IWEL)
              NE3=NE3+NWLN(IWEL)
              IF (IECHO1.NE.0) WRITE(IWFLX,944) IWEL,WNAME(IWEL),
     +            NWLN(IWEL),WLRATE(IWEL),WLTAB(IWEL),WLTAB(IWEL)+JLAST
              WLTAB(IWEL)=WLTAB(IWEL)+JLAST
C
C             CHECK NUMBER OF WELL LINKS.
              IF (NWLN(IWEL).GT.MAXWND) THEN
                  WRITE(*,947)
                  WRITE(*,916)
                  STOP
              ENDIF
C
C             WELL LINKS.
              DO I=1,NWLN(IWEL)
                  L2=L2+1
                  READ(IRFLX,945) IN2(L2,1),LEAK(IWEL,I)
                  IN2(L2,2)=NN+IWEL
                  IF (IECHO1.NE.0) WRITE(IWFLX,946) IN2(L2,1),
     +                IN2(L2,1)+ILAST,IN2(L2,2),LEAK(IWEL,I)
                  IN2(L2,1)=IN2(L2,1)+ILAST
              END DO
          END DO WELL_LOOP1
          ILAST=ILAST+NNC(ICELL)
          JLAST=JLAST+NTABC
      END DO CELL_LOOP
      NN3=NN+NWEL
      NN5=NN3
      NE4=NE3
      NE5=NE3
C
C     CLOSE FILE AND EXIT WFLUX
      WRITE(*,948)
      WRITE(IWFLX,950) MAX(NN,NN3,NN5),MAX(NE,NE2,NE3,NE4,NE5)
      CLOSE (IRFLX)
      RETURN
C
C     FORMATS:
902   FORMAT(1X,'*START READING IN WFLUX')
904   FORMAT(1X,'**START READING IN WFLUX**')
906   FORMAT(/1X,'FILE NAMES:'/
     +    1X,'INPUT FILE..............................',A80/
     +    1X,'OUTPUT FILE.............................',A80)
908   FORMAT(A20)
910   FORMAT(5I10)
912   FORMAT(/1X,'FLUX PARAMETERS:'/
     +    1X,'NUMBER OF COMPARTMENTS..................',I10/
     +    1X,'COMPARTMENT NUMBER......................',I10/
     +    1X,'COMPARTMENT NAME........................',A20/
     +    1X,'NUMBER OF NODES.........................',I10/
     +    1X,'NUMBER OF TABLES........................',I10/
     +    1X,'NUMBER OF WELL LINKS....................',I10/
     +    1X,'SWITCH FOR ECHO OF INPUTS...............',I10/
     +    1X,'SWITCH FOR OUTPUT OF STEP FLUXES........',I10)
914   FORMAT(1X,'*ERROR* COMPARTMENT NOT IN LIST')
916   FORMAT(1X,'*ERROR* RUN TERMINATED FOR ERROR IN WFLUX')
918   FORMAT(1X,'*ERROR* NUMBER OF NODES EXCEEDS DIMENSION')
920   FORMAT(1X,'*ERROR* NUMBER OF TABLES EXCEEDS DIMENSION')
922   FORMAT(/1X,'FLUX DATASET:'/
     +    1X,5X,'LOCAL',4X,'GLOBAL',15X,'LOCAL',4X,'GLOBAL'/
     +    1X,6X,'NODE',6X,'NODE',6X,'RATE',5X,'TABLE',5X,'TABLE')
924   FORMAT(I10,F10.0,I10,1X,A20)
926   FORMAT(1X,2I10,1PE10.2,2I10,1X,A20)
928   FORMAT(1X,'*ERROR* NUMBER OF NODES EXCEEDS DIMENSION')
930   FORMAT(1X,'*ERROR* TABLE NUMBER EXCEENS NUMBER OF TABLES')
932   FORMAT(2I10)
934   FORMAT(/1X,'SPECIFIED FLUX TABLE:'/
     +    1X,'LOCAL TABLE NUMBER......................',I10/
     +    1X,'GLOBAL TABLE NUMBER.....................',I10/
     +    1X,'NUMBER OF ENTRIES IN TABLE..............',I10/
     +    1X,'SWITCH FOR INTERPOLATION METHOD.........',I10//
     +    1X,5X,'POINT',6X,'TIME',6X,'RATE')
935   FORMAT(1X,'*ERROR* INTERPOLATION METHOD NOT VALID')
936   FORMAT(1X,'*ERROR* NUMBER OF POINTS EXCEEDS DIMENSION')
938   FORMAT(2F10.0)
940   FORMAT(1X,I10,2(1PE10.2))
942   FORMAT(/1X,'WELL PARAMETERS:'/
     +    1X,'NUMBER OF WELLS.........................',I10)
943   FORMAT(I10,F10.0,I10,1X,A20)
944   FORMAT(/1X,'WELL SPECIFICATIONS:'/
     +    1X,'WELL NUMBER.............................',I10/
     +    1X,'WELL NAME...............................',A20/
     +    1X,'NUMBER OF WELL NODES....................',I10/
     +    1X,'FLUX RATE...............................',1PE10.2/
     +    1X,'LOCAL FLUX RATE TABLE...................',I10/
     +    1X,'GLOBAL FLUX RATE TABLE..................',I10//
     +    1X,5X,'LOCAL',4X,'GLOBAL',6X,'WELL'/
     +    1X,6X,'NODE',6X,'NODE',6X,'NODE',2X,'LEAKANCE')
945   FORMAT(I10,F10.0)
946   FORMAT(1X,3I10,1PE10.2)
947   FORMAT(1X,'*ERROR* NUMBER OF WELL LINKS ECXEEDS DIMENSION')
948   FORMAT(1X,'*END   READING IN WFLUX')
950   FORMAT(/1X,'PROBLEM SIZE:'/
     +    1X,'NUMBER OF NODES.........................',I10/
     +    1X,'NUMNER OF ELEMENTS......................',I10//
     +    1X,'**END READING IN WFLUX**')
C
C     ---------------
      ENTRY WFLUX2(H)
C     ---------------
C
C     COMPUTE HEAD-INDEPENDENT FLUXES FOR TIME STEP.
      NODE_LOOP2:DO IFXN=1,NFXN
C
C         INTERPOLATE FROM TABLE.
          ITAB=FXTAB(IFXN)
          IF (ITAB.NE.0) THEN
              DO IPT=1,NPT(ITAB)
                  XTT(IPT)=TTB(ITAB,IPT)
                  XQT(IPT)=QTB(ITAB,IPT)
              END DO
              QFLX(IFXN)=
     +            TABLE(NPT(ITAB),XTT,XQT,STIME-DELT,STIME,IGET(ITAB))*
     +            FXRATE(IFXN)
          ELSE
              QFLX(IFXN)=FXRATE(IFXN)
          ENDIF
C
C         ASSIGN FLUX TO LOAD VECTOR.
          II=FXNODE(IFXN)
          F(II)=F(II)+QFLX(IFXN)
      END DO NODE_LOOP2
C
C     COMPUTE HEAD-DEPENDENT FLUXES FOR TIME STEP.
      L2=NE2
      WELL_LOOP2:DO IWEL=1,NWEL
C
C         INTERPOLATE FROM TABLE.
          ITAB=WLTAB(IWEL)
          IF (ITAB.NE.0) THEN
              DO IPT=1,NPT(ITAB)
                  XTT(IPT)=TTB(ITAB,IPT)
                  XQT(IPT)=QTB(ITAB,IPT)
              END DO
              QWEL(IWEL)=
     +            TABLE(NPT(ITAB),XTT,XQT,STIME-DELT,STIME,IGET(ITAB))*
     +            WLRATE(IWEL)
          ELSE
              QWEL(IWEL)=WLRATE(IWEL)
          ENDIF
          II=NN+IWEL
          F(II)=F(II)+QWEL(IWEL)
C
C         ASSIGN COEFFICIENTS.
          DO I=1,NWLN(IWEL)
              L2=L2+1
              CALL WSHAPE3(LEAK(IWEL,I),AE)
              DO II=1,2
                  IROW=IN2(L2,II)
                  DO JJ=1,2
                      ICOL=IN2(L2,JJ)
                      CALL PACK2(IROW,ICOL,AE(II,JJ),1)
                  END DO
              END DO
          END DO
      END DO WELL_LOOP2
      RETURN
C
C     ---------------------------
      ENTRY WFLUX3(H,SUMWL,SUMFX)
C     ---------------------------
C
C     OUTPUT SPECIFIED FLUXES.
      IF (ISTEP.EQ.ISTEP1) WRITE(IWFLX,952)
      SUMFX=0.0
      DO IFXN=1,NFXN
          SUMFX=SUMFX+QFLX(IFXN)
      END DO
      IF (IPRN.EQ.1 .AND. IECHO2.NE.0 .AND. NFXN.NE.0) THEN
          WRITE(IWFLX,954) ISTEP,STIME,SUMFX
          DO IFXN=1,NFXN
              WRITE(IWFLX,956) IFXN,FXNODE(IFXN),QFLX(IFXN),FXNAME(IFXN)
          END DO
      END IF
C
C     OUTPUT WELL FLUXES.
      SUMWL=0.0
      DO IWEL=1,NWEL
          SUMWL=SUMWL+QWEL(IWEL)
      END DO
      IF (IECHO2.NE.0 .AND. IPRN.EQ.1 .AND. NWEL.NE.0) THEN
          WRITE(IWFLX,958) ISTEP,STIME,SUMWL
          L2=NE2
          DO IWEL=1,NWEL
              II=NN+IWEL
              WRITE(IWFLX,960) IWEL,WNAME(IWEL),ISTEP,STIME,
     +            QWEL(IWEL),H(II)
              DO I=1,NWLN(IWEL)
                  L2=L2+1
                  II=IN2(L2,1)
                  JJ=IN2(L2,2)
                  QLNK(IWEL,I)=LEAK(IWEL,I)*(H(JJ)-H(II))
                  WRITE(IWFLX,962) I,II,QLNK(IWEL,I),H(II)
              END DO
          END DO
      ENDIF
      IF (ISTEP.EQ.ISTEP2) THEN
          WRITE(IWFLX,964)
          CLOSE (IWFLX)
      ENDIF
      RETURN
C
C     FORMATS:
952   FORMAT(1X,'**START SIMULATION IN WFLUX**')
954   FORMAT(/1X,'TIME-STEP SUMMARY FLUXES:'/
     +    1X,'TIME STEP...............................',I10/
     +    1X,'CUMULATIVE TIME.........................',1PE10.2/
     +    1X,'CUMULATIVE FLUX.........................',1PE10.2//
     +    1X,'NODE FLUXES:'/
     +    1X,5X,'POINT'6X,'NODE',6X,'FLUX')
956   FORMAT(1X,2I10,1PE10.2,1X,A20)
958   FORMAT(/1X,'TIME-STEP SUMMARY WELL FLUXES:'/
     +    1X,'TIME STEP...............................',I10/
     +    1X,'CUMULATIVE TIME.........................',1PE10.2/
     +    1X,'CUMULATIVE FLUX.........................',1PE10.2)
960   FORMAT(/1X,'WELL FLUXES:'/
     +    1X,'WELL NUMBER.............................',I10/
     +    1X,'WELL NAME...............................',A20/
     +    1X,'TIME STEP...............................',I10/
     +    1X,'CUMULATIVE TIME.........................',1PE10.2/
     +    1X,'WELL CUMULATIVE FLUX....................',1PE10.2/
     +    1X,'WELL HEAD...............................',0PF10.2//
     +    1X,5X,'POINT',6X,'NODE',6X,'FLUX',6X,'HEAD')
962   FORMAT(1X,2I10,1PE10.2,0PF10.2)
964   FORMAT(/1X,'**END SIMULATION IN WFLUX**')
      END SUBROUTINE
