C     -----------------
      PROGRAM FEMFLOW3D
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 SIMULATE FLOW AND TRANSPORT IN THREE-
C     DIMENSIONAL GOURNDWATER SYSTEMS USING THE FINITE ELEMENT METHOD.
C     THE FIRST ENTRY POINT IS CALLED ONCE BY THE MAIN PROGRAM.
C
C     INTERNAL VARIABLES.
      INCLUDE 'DIMEN.INC'
      INCLUDE 'WUNITS.INC'
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
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     X()         COORDINATE OF A NODE.
C     Y()         COORDINATE OF A NODE.
C     Z()         COORDINATE OF A NODE.
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     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     A()         VECTOR OF NON-ZERO COEFFICIENTS FOR CONDUCTANCE
C                     MATRIX.
C     B()         VECTOR OF NON-ZERO COEFFICIENTS FOR CAPACITANCE
C                     MATRIX.
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     NIA         NUMBER OF ROWS REPRESENTED IN VECTORS A() OR B().
C     NJA         NUMBER OF NON-ZERO COEFFICIENTS IN VECTORS A() OR B().
C     F()         RIGHTHAND VECTOR FOR SYSTEM OF EQUATIONS.
C     QX()        DARCY GROUNDWATER VELOCITY.
C     QY()        DARCY GROUNDWATER VELOCITY.
C     QZ()        DARCY GROUNDWATER VELOCITY.
C     ITRAN       SWITCH FOR TRANSPORT SIMULATION.
C
      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),SW(MAXNE2),POR(MAXNE),
     +    C(MAXNN),DC(MAXNN),RHO(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
C
C     INPUT DATA FOR FLOW MODULE.
      WRITE(*,902)
      CALL SYSTEM_CLOCK(IC,ICR,ICM)
      CLOCK1=IC/ICR
C
      CALL SETFILES1
      CALL BASIC(NITER,CONV,ISOLVE1,ISOLVE2,ISS,IMATRIX,ICHECK,
     +    DFMIN,SHAPE,RFACT)
      IF (ITRAN.EQ.1) CALL SETFILES2
      CALL NODES
      CALL ELEMS
      CALL WFLOW1(IMATRIX)
C     IF (ITRAN.EQ.1) CALL CTRAN1
      WRITE(*,904)
      IF (ICHECK.EQ.1) THEN
          WRITE(*,905)
          STOP
      ENDIF
C
C     BEGIN TIME STEP LOOP.
      DO ISTEP=ISTEP1,ISTEP2
          IF (ISTEP.NE.ISTEP1) ISS=0
          WRITE(*,906) ISTEP
          CALL STEPS2
          CALL WFLOW2(ISOLVE1,NITER,CONV,ISS,DFMIN,SHAPE,RFACT)
          DO ITSTEP=1,NTSTEP
              WRITE(*,907) ITSTEP
C         CALL CTRAN2(ISOLVE2)
          END DO
      END DO
C
      CALL SYSTEM_CLOCK(IC,ICR,ICM)
      CLOCK2=IC/ICR
      WRITE(*,908) CLOCK2-CLOCK1
      STOP
C
C     FORMATS:
902   FORMAT(1X,'*START READING INPUTS')
904   FORMAT(1X,'*END   READING INPUTS'/
     +    1X,'*START SIMULATION')
905   FORMAT(1X,'*CHECKING ONLY')
906   FORMAT(1X,'*START FLOW TIME STEP     ',I10)
907   FORMAT(1X,'*START TRANSPORT TIME STEP',I10)
908   FORMAT(1X,'*SIMULATION ELAPSED TIME  ',F10.1)
909   FORMAT(/1X,'*SIMULATION ELAPSED TIME  ',F10.1)
      END PROGRAM
      INCLUDE 'SETFILES.FOR'
      INCLUDE 'OPENFILES.FOR'
      INCLUDE 'TABLE.FOR'
      INCLUDE 'BASIC.FOR'
      INCLUDE 'STEPS.FOR'
      INCLUDE 'NODES.FOR'
      INCLUDE 'ELEMS.FOR'
      INCLUDE 'PACK.FOR'
      INCLUDE 'WFLOW.FOR'
      INCLUDE 'SRELAX.FOR'
      INCLUDE 'SITPACK.FOR'
      INCLUDE 'SAMG.FOR'
