C     Last change:  TJD  30 Dec 2005    4:53 pm
C     --------------------------
      SUBROUTINE WFLOW1(IMATRIX)
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 THE CENTRAL MODULE OF THE GROUNDWATER-FLOW
C     COMPONENT WITHIN FEMFLOW3D AND CONTROLS THE DATA INPUT AND
C     EXECUTES THE PRINCIPAL CALCULATINS RELATING TO THE IMPLEMENTATION
C     FINITE ELEMENT METHOD. THE FIRST ENTRY POINT IS CALLED ONCE BY
C     THE MAIN PROGRAM. THE SECOND ENTRY POINT IS CALLED ONCE FOR
C     EACH TIMESTEP ITERATION.
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     H0()        HEADS AT NODES.
C     H()         HEADS AT NODES.
C     HL()        HEADS AT NODES.
C     HLIT()      HEADS AT NODES.
C
C     INTERNAL ARRAYS.
      REAL*8 H0(MAXNN),H(MAXNN),HL(MAXNN),HLIT(MAXNN)
      REAL*8 KXX(MAXNE),KYX(MAXNE),KYY(MAXNE),KZX(MAXNE),KZY(MAXNE),
     +    KZZ(MAXNE),SS(MAXNE)
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
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)
      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
      SAVE
C
C     BEGIN READING FLOW INPUTS.
      CALL WPARAMS1
      CALL WSTART1(H)
      CALL WOUTPUT1(H)
      IF (FNAME(7).NE.'NOTUSED') CALL WCHEAD1
      IF (FNAME(8).NE.'NOTUSED') CALL WFLUX1
      IF (FNAME(9).NE.'NOTUSED') CALL WEVAP1
      IF (FNAME(10).NE.'NOTUSED') CALL WRIVER1
      IF (FNAME(11).NE.'NOTUSED') CALL WFAULT1(H)
      CALL PACK1(IMATRIX)
      RETURN
C
C     -------------------------------------------------------
      ENTRY WFLOW2(ISOLVE,NITER,CONV,ISS,DFMIN1,SHAPE1,RFACT)
C     -------------------------------------------------------
C
C     INITIALIZATIONS.
      DO II=1,NN5
          HL(II)=H(II)
          HLIT(II)=H(II)
      END DO
      IF (ISTEP.EQ.ISTEP1 .AND. ISS.EQ.1) THEN
          FSS=0.0
      ELSE
          FSS=1.0
      ENDIF
C
C     BEGIN TIME-STEP ITERATION LOOP.
      DO ITER=1,NITER
C
C         CONSTRUCT MATRICIES FOR SYSTEM OF EQUATIONS.
          DO IJA=1,NJA
              A(IJA)=0.0
              B(IJA)=0.0
          END DO
          DO II=1,NN5
              F(II)=0.0
          END DO
          CALL NODES2(H)
          CALL WPARAMS2(H,ISS)
          IF (FNAME(7).NE.'NOTUSED') CALL WCHEAD2(H)
          IF (FNAME(8).NE.'NOTUSED') CALL WFLUX2(H)
          IF (FNAME(9).NE.'NOTUSED') CALL WEVAP2(H)
          IF (FNAME(10).NE.'NOTUSED') CALL WRIVER2(H)
          IF (FNAME(11).NE.'NOTUSED') CALL WFAULT2(H,ISS)
C
C         MULTIPLY MAXTRIX [B] BY MAXTRIX [HL].
          DO II=1,NN5
              IJA1=IA(II)
              IJA2=IA(II+1)-1
              SROW=0.0
              DO IJA=IJA1,IJA2
                  JJ=JA(IJA)
                  SROW=SROW+FSS*B(IJA)*HL(JJ)
              END DO
              F(II)=F(II)+SROW/DELT
          END DO
C
C         ADD MATRICIES [A] AND [B].
          DO IJA=1,NJA
              A(IJA)=A(IJA)+FSS*B(IJA)/DELT
          END DO
C
C         SOLVE SYSTEM OF EQUATIONS.
          IF (ISOLVE.EQ.1) THEN
              CALL SRELAX2(IA,JA,A,H,F,NN5,ISTEP,JITER,1)
          ELSEIF (ISOLVE.EQ.2) THEN
              CALL SITPACK2(IA,JA,A,H,F,NN5,ISTEP,JITER,1)
          ELSEIF (ISOLVE.EQ.3) THEN
              CALL SAMG2(IA,JA,A,H,F,NN5,ISTEP,JITER,1)
          ENDIF
C
C         TEST FOR TIME-STEP CLOSURE.
C         IDENTIFY MAXIMUM DIFFERENCE IN HEAD.
          DIF=0.0
          DO II=1,NN5
              DH=H(II)-HLIT(II)
              DIF=MAX(DIF,ABS(DH))
          END DO
C
C         COMPUTE UNDERRELAXATION FACTOR.
          IF (ITER.EQ.1) THEN
C
C             FIRST ITERATION.
              DFMIN=DFMIN1
              SHAPE=SHAPE1
              DF=DFMIN+(1.0-DFMIN)*EXP(-SHAPE*(DIF-CONV))
          ELSE
C
C             SUBSEQUENT ITERATIONS.
              IF (DIF.LT.DIFLAST) THEN
C
C                 CONVERGENCE FROM LAST ITERATION.
                  DF=DFMIN+(1.0-DFMIN)*EXP(-SHAPE*(DIF-CONV))
              ELSE
C
C                 DIVERGENCE FROM LAST ITERATION.
                  DDF=(1.0-DFMIN)*EXP(-SHAPE*(DIF-CONV))*RFACT
                  DFMIN=DFMIN*RFACT
                  IF (DDF.GT.1.0E-04) THEN
                      SHAPE=-LOG(DDF/(1.0-DFMIN))/(DIF-CONV)
                  ENDIF
                  DF=DFMIN+(1.0-DFMIN)*EXP(-SHAPE*(DIF-CONV))
              ENDIF
          ENDIF
          DIFLAST=DIF
          WRITE(*,906) ITER,DIF,DF,DFMIN,SHAPE
          IF (DIF.LE.CONV) EXIT
C
C         UPDATE HEADS.
          DO II=1,NN5
              DH=H(II)-HLIT(II)
              H(II)=HLIT(II)+DH*DF
              HLIT(II)=H(II)
          END DO
      END DO
      ITER=MIN(ITER,NITER)
C
C     CONSTRUCT WATER BUDGET.
      SUMCH=0.0
      SUMFX=0.0
      SUMWL=0.0
      SUMET=0.0
      SUMRV=0.0
      SUMLK=0.0
      SUMFS=0.0
      SUMPS=0.0
      CALL NODES3(ITER)
      CALL WPARAMS3(H,HL,SUMPS)
      IF (FNAME(7).NE.'NOTUSED') CALL WCHEAD3(H,SUMCH)
      IF (FNAME(8).NE.'NOTUSED') CALL WFLUX3(H,SUMWL,SUMFX)
      IF (FNAME(9).NE.'NOTUSED') CALL WEVAP3(H,SUMET)
      IF (FNAME(10).NE.'NOTUSED') CALL WRIVER3(H,SUMRV)
      IF (FNAME(11).NE.'NOTUSED') CALL WFAULT3(H,HL,SUMLK,SUMFS)
      CALL WOUTPUT2(H,SUMCH,SUMFX,SUMWL,SUMET,SUMRV,SUMLK,SUMPS,SUMFS,
     +    ITER,DIFF)
      CALL WSTART2(H)
C
C     CALCULATE DARCY VELOCITIES.
      CALL WPARAMS4(H)
      RETURN
C
C     FORMAT:
906   FORMAT(1X,' STEP ITERATION=',I4,' HEAD CHANGE=',1PE9.2,
     +    ' RELAXATION FACTOR=',1PE9.2/
     +    1X,42X,' MINIMUM DAMPING=  ',1PE9.2/
     +    1X,42X,' SHAPE FACTOR=     ',1PE9.2)
      END SUBROUTINE
      INCLUDE 'WOUTPUT.FOR'
      INCLUDE 'WPARAMS.FOR'
      INCLUDE 'WSTART.FOR'
      INCLUDE 'WCHEAD.FOR'
      INCLUDE 'WFLUX.FOR'
      INCLUDE 'WEVAP.FOR'
      INCLUDE 'WRIVER.FOR'
      INCLUDE 'WFAULT.FOR'
      INCLUDE 'WSHAPE.FOR'
