C     ---------------------
      SUBROUTINE WFAULT1(H)
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 GROUNDWATER FLOW IN A TWO-
C     DIMENSIONAL FAULT PLANE WITHIN THE THREE-DIMENSIONAL FINITE-
C     ELEMENT MESH, INCLUDING THE EXCHANGE OF WATER BETWEEN THE FAULT
C     PLANE AND THE ADJACENT THREE-DIMENSIONAL GROUNDWATER SYSTEM.
C     THE FIRST ENTRY POINT IS CALLED ONCE BYTHE MAIN PROGRAM.
C     THE SECOND ENTRY POINT IS CALLED ONCE FOR EACH TIMESTEP ITERATION.
C     THE THIRD ENTRY POINT IS CALLED 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     MESHFILE    SWITCH FOR OUTPUT OF MESH AS .3DM FILE FOR GMS.
C     XCNAME      COMPARTMENT NAME.
C     H()         HEAD VALUE AT NODE.
C     HL()        HEAD VALUE AT NODE.
C     TRAN()      TRANSMISSIVITY OF FAULT ELEMENT NORMAL TO FAULT PLANE.
C     STOR()      STORATIVITY OF FAULT ELEMENT.
C     LEAK()      LEAKANCE ASSIGNED TO ELEMENTS WITHIN GROUP.
C     XNFN        NUMBER OF FAULT NODES IN FNT.
C     XE()        X-COORDINATE FOR FAULT NODE.
C     YE()        Y-COORDINATE FOR FAULT NODE.
C     ZE()        Z-COORDINATE FOR FAULT NODE.
C     AE()        ELEMENTAL CONDUCTANCE MATRIX.
C     BE()        ELEMENTAL STORAGE MATRIX.
C     AREA()      AREA OF ELEMENT
C     IMESH       SWITCH FOR OUTPUT OF 2DM MESH.
C     NFN         NUMBER OF FAULT NODES.
C     NFE         NUMBER OF FAULT ELEMENTS.
C     NFL         NUMBER OF FAULT LINKS.
C     IECHO1      SWITCH FOR ECHO OF FAULT INPUTS.
C     IECHO2      SWITCH FOR ECHO OF FAULT INITIAL HEADS.
C     IECHO3      SWITCH FOR ECHO OF TIME-STEP FAULT HEADS.
C     IECHO4      SWITCH FOR ECHO OF TIME-STEP FAULT LINK FLUXES.
C     FACX        SCALING FACTOR FOR X-COORDINATES FOR FAULT NODES.
C     FACY        SCALING FACTOR FOR Y-COORDINATES.
C     FACZ        SCALING FACTOR FOR Z-COORDINATES.
C     X()         X-COORDINATE FOR FAULT NODE.
C     Y()         Y-COORDINATE FOR FAULT NODE.
C     Z()         Z-COORDINATE FOR FAULT NODE.
C     NFC         NUMBER OF FAULT NODE COLUMNS
C     NMAX        MAXIMUM NUMBER OF NODES IN A COLUMN.
C     ZIMN()      MINIMUM DISTANCE BETWEEN COLLAPSED NODES.
C     ZMAX()      MAXIMUM ELEVATION OF TOP NODE.
C     COL()       COLUMN NODE NUMBER.
C     IN2()       NODE IN FAULT ELEMENT.
C     MAT()       IDENTIFER FOR FAULT MATERAIAL.
C     IN2(,1)     NODE IN FAULT PLANE.
C     IN2(,2)     NODE IN ADJACENT GROUNDWATER SYSTEM.
C     NGROUP      NUMBER OF FAULT GROUPS.
C     FACT        SCALING FACTOR FOR FAULT TRANSMISSIVITY.
C     FACS        SCALING FACTOR FOR FAULT STORAGE COEFFICIENT.
C     LL1         FIRST ELEMENT NUMBER WITHIN RANGE FOR ELEMENT GROUP.
C     LL2         LAST ELEMENT NUMBER WITHIN RANGE FOR ELEMENT GROUP.
C     XTRAN()     TRANSMISSIVITY OF FAULT ELEMENT NORMAL TO FAULT PLANE.
C     XSTORE()    STORATIVITY OF FAULT ELEMENT.
C     FACL        SCALING FACTOR FOR LEAKANCEASSIGNED TO FAULT LINKS.
C     LL1         FIRST ELEMENT NUMBER WITHIN RANGE FOR ELEMENT GROUP.
C     LL2         LAST ELEMENT NUMBER WITHIN RANGE FOR ELEMENT GROUP.
C     XLEAK       LEAKANCE ASSIGNED TO ELEMENTS WITHIN GROUP.
C     ISTART      SWITCH FOR UNIFORM INITAL HEADS.
C     XH0         UNIFORM INITAL HEAD VALUE.
C     ISAVE       SWITCH FOR OVERWRITING INITIAL HEADS.
C     IWHEN       STEP FOR OVERWRITING INITIAL HEADS.
C
C     INTERNAL ARRAYS.
      CHARACTER*20 XCNAME
      CHARACTER*80 INFILE,OUTFILE,MESHFILE,OUTEXT
      REAL*8 H(MAXNN),HL(MAXNN)
      REAL*8 TRAN(MAXFLE),STOR(MAXFLE),LEAK(MAXFLL)
      INTEGER*4 MAT(MAXFLE),XNFN
      LOGICAL FOUND,DONE
      REAL*8 XE(4),YE(4),ZE(4)
      REAL*8 AE(4,4),BE(4,4),AREA(MAXFLE)
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     IRFLT       INTERNAL READ FOR FLT UNIT NUMBER.
C     IWFLT       INTERNAL WRITE FOR FLT UNIT NUMBER.
C     Z0          ELEVATION OF COLUMN NODE.
C     COL         COLUMN NODE NUMBER.
C     NLAY        NUMBER OF COLUMN LAYERS.
C     NCOL        NUMBER OF COLUMN NODES.
C     ZMIN        MINIMUM DISTANCE BETWEEN COLLAPSED NODES.
C     ZMAX        MAXIMUM ELEVATION OF TOP NODE.
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 COL(MAXCOL,MAXLAY),NLAY(MAXCOL)
      REAL*8 ZMIN(MAXCOL),ZMAX(MAXCOL),Z0(MAXNN)
      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 /ADJUST/Z0,COL,NLAY,NCOL,ZMIN,ZMAX
      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     OPEN INPUT AND OUTPUT FILES.
      WRITE(*,900)
      OUTEXT='.FLT.OUT'
      INFILE=FNAME(11)
      CALL OPENFILES(RNAME1,INFILE,OUTFILE,OUTEXT,IRFLT,IWFLT)
      OUTEXT='.FLT.PLT'
      INFILE=FNAME(11)
      CALL OPENFILES(RNAME1,INFILE,OUTFILE,OUTEXT,0,IWFLTP)
      OUTEXT='.FNT.OUT'
      INFILE=FNAME(12)
      CALL OPENFILES(RNAME1,INFILE,OUTFILE,OUTEXT,IRFNT,IWFNT)
      WRITE(IWFLT,901)
      WRITE(IWFLT,902) INFILE,OUTFILE
      READ(IRFLT,903) IMESH
      IF (IMESH.NE.0) THEN
          OUTEXT='.FLT.2DM'
          CALL OPENFILES(RNAME1,INFILE,MESHFILE,OUTEXT,0,IW2DM)
          WRITE(IWFLT,904) MESHFILE
      ENDIF
C
C     INPUT FAULT DIMENSIONS.
      READ(IRFLT,906) NFN,NFE,NFL,IECHO1,IECHO2,IECHO3,IECHO4
      WRITE(IWFLT,908) NFN,NFE,NFL,IECHO1,IECHO2,IECHO3,IECHO4,
     +    IMESH
      NN5=NN3+NFN
      NE4=NE3+NFE
      NE5=NE4+NFL
C
C     CHECK NUMBER OF FAULT NODES.
      IF (NN5.GT.MAXNN) THEN
          WRITE(*,910)
          WRITE(*,912)
          STOP
      ENDIF
C
C     CHECK NUMBER OF FAULT ELEMENTS.
      IF (NE5.GT.MAXNE2) THEN
          WRITE(*,916)
          WRITE(*,912)
          STOP
      ENDIF
      IF (NFE.GT.MAXFLE) THEN
          WRITE(*,918)
          WRITE(*,912)
          STOP
      ENDIF
C
C     CHECK NUMBER OF LINK ELEMENTS.
      IF (NFL.GT.MAXFLL) THEN
          WRITE(*,919)
          WRITE(*,912)
          STOP
      ENDIF
C
C     INPUT FAULT NODES.
      READ(IRFLT,920) FACX,FACY,FACZ
      IF (IECHO1.NE.0) WRITE(IWFLT,922) FACX,FACY,FACZ
      I=NN3
      DO IFN=1,NFN
          I=I+1
          READ(IRFLT,924) X(I),Y(I),Z(I)
          X(I)=X(I)*FACX
          Y(I)=Y(I)*FACY
          Z(I)=Z(I)*FACZ
          Z0(I)=Z(I)
          IF (IECHO1.NE.0) WRITE(IWFLT,926) IFN,I,X(I),Y(I),Z(I)
      END DO
C
C     INPUT PARAMETERS FOR COLLAPSING MESH.
      READ(IRFLT,1100) NFC,NMAX
      WRITE(IWFLT,1102) NFC,NMAX
      IF (NFC.NE.0) THEN
          ICOL=NCOL
          DO IFC=1,NFC
C
C             INPUT NODE COLUMN.
              ICOL=ICOL+1
              READ(IRFLT,1104) ZMIN(ICOL),ZMAX(ICOL),(COL(ICOL,ILAY),
     +            ILAY=1,NMAX)
C
C             CHECK THAT ZMAX EAUAL TO OR HIGHER THAN TOP NODE.
                  II=COL(ICOL,1)+NN3
                  IF (ZMAX(ICOL).LT.Z(II)) THEN
                      WRITE(*,1110)
                      WRITE(*,912)
                      STOP
                  ENDIF
C
C             LOOP OVER LAYERS.
              DO ILAY=1,NMAX
C
C                 IDENTIFY NUMBER OF NODES IN COLUMN.
                  IF (COL(ICOL,ILAY).EQ.0) THEN
                      NLAY(ICOL)=ILAY-1
                      EXIT
                  ENDIF
                  COL(ICOL,ILAY)=COL(ICOL,ILAY)+NN3
                  IF (ILAY.EQ.NMAX) THEN
                      NLAY(ICOL)=ILAY
                      EXIT
                  ENDIF
              END DO
C
C             OUTPUT NODE COLUMN.
              WRITE(IWFLT,1106) IFC,ZMIN(ICOL),ZMAX(ICOL),
     +            (COL(ICOL,ILAY)-NN3,ILAY=1,NLAY(ICOL))
              WRITE(IWFLT,1108) IFC,(COL(ICOL,ILAY),ILAY=1,NLAY(ICOL))
          END DO
          NCOL=ICOL
      ENDIF
C
C     INPUT FAULT ELEMENTS.
      IF (IECHO1.NE.0) WRITE(IWFLT,928)
      L2=NE3
      DO IFE=1,NFE
          L2=L2+1
          READ(IRFLT,930) (IN2(L2,J),J=1,3),MAT(IFE)
          IF (IECHO1.NE.0) WRITE(IWFLT,932) IFE,L2,(IN2(L2,J),J=1,3),
     +        (IN2(L2,J)+NN3,J=1,3),MAT(IFE)
          DO J=1,3
              IF (IN2(L2,J).GT.NFN) THEN
                  WRITE(*,934)
                  WRITE(*,912)
                  STOP
              ENDIF
              IN2(L2,J)=IN2(L2,J)+NN3
          END DO
      END DO
C
C     INPUT FALUT LINKS.
C         IN2(*,1) IS A NODE IN A TWO-DIMENSIONAL FALUT ELEMENT.
C         IN2(*,2) IS A NODE IN A THREE-DIMENSIONAL AQUIFER ELEMENT.
      IF (IECHO1.EQ.1) WRITE(IWFLT,938)
      L2=NE4
      DO IFL=1,NFL
C
C         INPUT LOCAL ELEMENT INCIDENCES.
          L2=L2+1
          READ(IRFLT,940) (IN2(L2,J),J=1,2),XCNAME
C
C         ADJUST LOCAL FAULT NODE NUMBER TO GLOBAL NUMBER.
          IN2(L2,1)=IN2(L2,1)+NN3
C
C         ADJUST LOCAL AQUIFER NODE NUMBER TO GLOBAL NUMBER.
          FOUND=.FALSE.
          ILAST=0
          DO ICELL=1,NCELL
              IF (CNAME(ICELL).EQ.XCNAME) THEN
                  FOUND=.TRUE.
                  EXIT
              ENDIF
              ILAST=ILAST+NNC(ICELL)
          END DO
          IF (FOUND) THEN
              IN2(L2,2)=IN2(L2,2)+ILAST
          ELSE
              WRITE(*,944)
              WRITE(*,912)
              STOP
          ENDIF
          IF (IECHO1.NE.0) WRITE(IWFLT,942) IFL,L2,IN2(L2,1)-NN3,
     +        IN2(L2,2)-ILAST,(IN2(L2,J),J=1,2),XCNAME
      END DO
C
C     INPUT FAULT HYDRAULIC PROPERTIES
      READ(IRFLT,948) NGROUP,FACT,FACS
      IF (IECHO1.NE.0) WRITE(IWFLT,950) NGROUP,FACT,FACS
      DO IGROUP=1,NGROUP
          READ(IRFLT,952) LL1,LL2,XTRAN,XSTOR
          XTRAN=XTRAN*FACT
          XSTOR=XSTOR*FACS
          L1=LL1+NE3
          L2=LL2+NE3
          IF (IECHO1.NE.0) WRITE(IWFLT,954) LL1,LL2,L1,L2,XTRAN,XSTOR
          IF (LL1.GT.NFE .OR. LL2.GT.NFE) THEN
              WRITE(*,956)
              WRITE(*,912)
              STOP
          ENDIF
          DO IFE=LL1,LL2
              TRAN(IFE)=XTRAN
              STOR(IFE)=XSTOR
          END DO
      END DO
C
C     INPUT LINK HYDRUALIC PROPERTIES.
      READ(IRFLT,980) NGROUP,FACL
      IF (IECHO1.NE.0) WRITE(IWFLT,982) NGROUP,FACL
      DO IGROUP=1,NGROUP
          READ(IRFLT,983) LL1,LL2,XLEAK
          XLEAK=XLEAK*FACL
          L1=LL1+NE4
          L2=LL2+NE4
          IF (IECHO1.NE.0) WRITE(IWFLT,984) LL1,LL2,L1,L2,XLEAK
          IF (LL1.GT.NFL .OR. LL2.GT.NFL) THEN
              WRITE(*,985)
              WRITE(*,912)
              STOP
          ENDIF
          DO IFL=LL1,LL2
              LEAK(IFL)=XLEAK
          END DO
      END DO
C
C     OUTPUT 2DM FILE.
      IF (IMESH.NE.0) THEN
          WRITE(IW2DM,986)
          L2=NE3
          DO IFE=1,NFE
              L2=L2+1
              WRITE(IW2DM,987) IFE,(IN2(L2,J)-NN3,J=1,3),MAT(IFE)
          END DO
          I=NN3
          DO IFN=1,NFN
              I=I+1
              WRITE(IW2DM,988) IFN,X(I),Y(I),Z(I)
          END DO
      ENDIF
C
C     INPUT FAULT INITIAL HEADS
      READ(IRFNT,972) XNFN,ISTART,XH0,ISAVE,IWHEN
      IF (IECHO1.NE.0) WRITE(IWFNT,974) XNFN,ISTART,XH0,ISAVE,IWHEN
      I1=NN3+1
      I2=NN3+NFN
      IF (ISTART.EQ.1) THEN
          DO I=I1,I2
              H(I)=XH0
          END DO
      ELSE
          READ(IRFNT,976) (H(I),I=I1,I2)
      ENDIF
C
C     OUTPUT INITIAL HEADS
      IF (IECHO1.EQ.1) THEN
          DONE=.FALSE.
          I1=NN3-4
          DO WHILE (.NOT.DONE)
              I1=I1+5
              I2=I1+4
              IF (I2.GE.NN5) THEN
                  I2=NN5
                  DONE=.TRUE.
              ENDIF
              WRITE(IWFNT,978) I1-NN3,I2-NN3,I1,I2,(H(I),I=I1,I2)
          END DO
      ENDIF

C
C     CLOSE INPUT FILE AND EXIT WFAULT.
      WRITE(*,989)
      WRITE(IWFLT,990)
      CLOSE (IRFLT)
      IF(ISAVE.EQ.0)THEN
          CLOSE(IRFNT)
          CLOSE(IWFNT)
      ELSE
          REWIND(IRFNT)
          CLOSE(IWFNT)
      ENDIF
      IF (IMESH.NE.0) CLOSE (IW2DM)
C
C     WRITE PLOT HEADER AND EXIT WOUTPUT.
      WRITE(IWFLTP,991)
      WRITE(IWFLTP,992)
      WRITE(IWFLTP,993)
      WRITE(IWFLTP,994) NFN
      WRITE(IWFLTP,995) NFE
      WRITE(IWFLTP,996)
      WRITE(IWFLTP,997) 0
      II=NN3
      DO I=1,NFN
          II=II+1
          WRITE(IWFLTP,998) H(II)
      END DO
      RETURN
C
C     FORMATS:
900   FORMAT(1X,'*START READING IN WFAULT')
901   FORMAT(1X,'**START READING IN WFAULT**')
902   FORMAT(/1X,'FILE NAMES:'/
     +    1X,'INPUT FILE..............................',A80/
     +    1X,'OUTPUT FILE.............................'A80)
903   FORMAT(I10)
904   FORMAT(1X,'MESH 2DM FILE............................',A80)
906   FORMAT(7I10)
908   FORMAT(/1X,'FAULT PARAMETERS:'/
     +    1X,'NUMBER OF FAULT NODES...................',I10/
     +    1X,'NUMBER OF FAULT ELEMENTS................',I10/
     +    1X,'NUMBER OF FALUT LINKS...................',I10/
     +    1x,'SWITCH FOR ECHO OF INPUTS...............',I10/
     +    1X,'SWITCH FOR OUTPUT OF FAULT FLUXES.......',I10/
     +    1X,'SWITCH FOR ECHO OF INITIAL HEADS........',I10/
     +    1X,'SWITCH FOR OUTPUT OF HEADS..............',I10/
     +    1X,'SWITCH FOR OUTPUT OF MESH FILE (2DM)....',I10)
910   FORMAT(1X,'*ERROR* NUMBER OF NODES EXCEEDS DIMENSION')
912   FORMAT(1X,'*ERROR* SIMULATION TERMINATED FOR ERROR IN WFALUT')
916   FORMAT(1X,'*ERROR* NUMBER OF ELEMENTS EXCEEDS DIMESION')
918   FORMAT(1X,'*ERROR* NUMBER OF FAULT ELEMENTS EXCEEDS DIMENSION')
919   FORMAT(1X,'*ERROR* NUMBER OF LINKS EXCEEDS DIMENSION')
920   FORMAT(3F10.0)
922   FORMAT(/1X,'FAULT NODES PARAMETERS:'/
     +    1X,'FACTOR FOR X COORDINATE.................',1PE10.2/
     +    1X,'FACTOR FOR Y COORDINATE.................',1PE10.2/
     +    1X,'FACTOR FOR Z COORDINATE.................',1PE10.2//
     +    1X,'NODE COORDINATES:'//
     +    1X,5X,'LOCAL',4X,'GLOBAL'/
     +    1X,6X,'NODE',6X,'NODE',9X,'X',9X,'Y',9X,'Z')
924   FORMAT(10X,3F10.0)
926   FORMAT(1X,2I10,3F10.1)
928   FORMAT(/1X,'FAULT ELEMENT LOCAL/GLOBAL INCIDENCES:'//
     +    1X,5X,'LOCAL',4X,'GLOBAL',3(5X,'LOCAL'),3(4X,'GLOBAL')/
     +    1X,6X,'ELEM',6X,'ELEM',6(6X,'NODE'),2X,'MATERIAL')
930   FORMAT(10X,4I10)
932   FORMAT(1X,9I10)
934   FORMAT(1X,'*ERROR* NODE NUMBER EXCEEDS NUMBER OF NODES')
938   FORMAT(/1X,'LINK ELEMENT LOCAL/GLOBAL INCIDENCES:'//
     +    1X,5X,'LOCAL',4X,'GLOBAL',2(5X,'LOCAL'),2(4X,'GLOBAL')/
     +    1X,6X,'ELEM',6X,'ELEM',4(6X,'NODE'),1X,'COMPARTMENT')
940   FORMAT(10X,2I10,1X,A20)
942   FORMAT(1X,6I10,1X,A20)
944   FORMAT(1X,'*ERROR* COMPARTMENT NAME NOT IN LIST')
946   FORMAT(1X,20X,2I10)
948   FORMAT(I10,2F10.0)
950   FORMAT(/1X,'FAULT ELEMENT HYDRAULIC PROPERTIES:'/
     +    1X,'NUMBER OF GROUPS........................',I10/
     +    1X,'FACTOR FOR TRANSMISSIVITY...............',1PE10.2/
     +    1X,'FACTOR FOR STORATIVITY..................',1PE10.2//
     +    1X,5X,'LOCAL',5X,'LOCAL',4X,'GLOBAL',4X,'GLOBAL'/
     +    1X,4(6X,'ELEM'),6X,'TRAN',6X,'STOR',5X,'TABLE')
952   FORMAT(2I10,2F10.0,I10)
954   FORMAT(1X,4I10,2(1PE10.2),I10)
956   FORMAT(1X,'*ERROR* FAULT ELEMENT RANGE EXCEEDS DIMENSION')
958   FORMAT(/1X,'UNSATURATED PARAMETER TABLES:'/
     +    1X,'NUMBER OF TABLES........................',I10)
960   FORMAT(I10,F10.0)
962   FORMAT(/1X,'TABLE PARAMETERS:'/
     +    1X,'TABLE NUMBER............................',I10/
     +    1X,'NUMBER OF POINTS IN TABLE...............',I10/
     +    1X,'MAXIMUM INTEGRATION SPACING.............',F10.1//
     +    5X,'POINT',3X,'PRESSURE',8X,'KR')
964   FORMAT(1X,'*ERROR* NUMBER OF POINTS EXCEEDS DIMENSIONS')
965   FORMAT(1X,'*ERROR* MAXIMUM SPACING LESS THAN OR TO EQUAL ZERO')
968   FORMAT(3F10.0)
970   FORMAT(1X,I10,2(1PE10.2))
972   FORMAT(2I10,F10.0,2I10)
974   FORMAT(/1X,'FAULT INITIAL HEADS:'/
     +    1X,'NUMBER OF FAULT NODES...................',I10/
     +    1X,'SWICH FOR UNIFORM INITIAL HEADS.........',I10/
     +    1X,'UNIFORM INITIAL HEAD....................',F10.2/
     +    1X,'SWITCH FOR OVERWRITTING INITIAL HEADS...',I10/
     +    1X,'STEP FOR OVERWRITTING INITIAL HEADS.....',I10//
     +    1X,2(5X,'LOCAL'),2(4X,'GLOBAL')/
     +    1X,4(6X,'NODE'),5(6X,'HEAD'))
976   FORMAT(20X,5F10.0)
978   FORMAT(1X,4I10,5F10.2)
980   FORMAT(I10,F10.0)
982   FORMAT(/1X,'LINK ELEMENT HYDRAULIC PARAMETERS:'/
     +    1X,'NUMBER OF GROUPS........................',I10/
     +    1X,'FACTOR FOR LEAKANCE.....................',1PE10.2//
     +    1X,5X,'LOCAL',5X,'LOCAL',4X,'GLOBAL',4X,'GLOBAL'/
     +    1X,4(6X,'ELEM'),2X,'LEAKANCE')
983   FORMAT(2I10,F10.0)
984   FORMAT(1X,4I10,1PE10.2)
985   FORMAT(1X,'*ERROR* LINK ELEMENT RANGE EXCEED DIMENSION')
986   FORMAT('MESH2D')
987   FORMAT('E3T',7X,5I10)
988   FORMAT('ND',8X,I10,3F20.4)
989   FORMAT(1X,'*END   READING IN WFALUT')
990   FORMAT(/1X,'**END READING IN WFAULT**')
991   FORMAT('DATASET')
992   FORMAT('OBJTYPE mesh2d')
993   FORMAT('BEGSCL')
994   FORMAT('ND',I13)
995   FORMAT('NC',I13)
996   FORMAT('NAME "Scalar Heads"')
997   FORMAT('TS 0',I11)
998   FORMAT(F15.5)
1100  FORMAT(2I10)
1102  FORMAT(/1X,'COLLAPSING MESH PARAMETERS:'/
     +    1X,'NUMBER OF NODE COLUMNS.................',I10/
     +    1X,'MAXIMUM NODES IN COLUMN................',I10//
     +    1X,4X,'COLUMN',6X,'ZMIN',6X,'ZMAX',5X,'NODES')
1104  FORMAT(2F10.0,20I10)
1106  FORMAT(1X,I10,2F10.1,20I10)
1108  FORMAT(1X,I10,20X,20I10)
1110  FORMAT(1X,'*ERROR* ZMAX LOWER THAN TOP NODE')
C
C     --------------------
      ENTRY WFAULT2(H,ISS)
C     --------------------
C
C     TRIANGUALER ELEMENTS.
      L2=NE3
      DO IFE=1,NFE
          L2=L2+1
C
C         ASSIGN COORDINATES
          DO I=1,3
              II=IN2(L2,I)
              XE(I)=X(II)
              YE(I)=Y(II)
              ZE(I)=Z(II)
          END DO
C
C         VECTORS A AND B AND MAGNITUDE.
          A1=XE(2)-XE(1)
          A2=YE(2)-YE(1)
          A3=ZE(2)-ZE(1)
          B1=XE(3)-XE(1)
          B2=YE(3)-YE(1)
          B3=ZE(3)-ZE(1)
          AMAG=(A1*A1+A2*A2+A3*A3)**0.5
          BMAG=(B1*B1+B2*B2+B3*B3)**0.5
C
C         ANGLE BETWEEN VECTORS.
          ABDOT=A1*B1+A2*B2+A3*B3
          ANGLE=ACOS(ABDOT/(AMAG*BMAG))
C
C         LOCAL COORDINATES.
          XE(1)=0.0
          YE(1)=0.0
          XE(2)=AMAG
          YE(2)=0.0
          XE(3)=BMAG*COS(ANGLE)
          YE(3)=BMAG*SIN(ANGLE)
C
C         INTEGRATE OVER PLANE OF FAULT.
          IF (ISS.EQ.1 .AND. ISTEP.EQ.ISTEP1) THEN
              FSS=0.0
          ELSE
              FSS=1.0
          ENDIF
          CALL WSHAPE2(TRAN(IFE),STOR(IFE)*FSS,XE,YE,AE,BE,
     +        AREA(IFE),IFE)
          AREA(IFE)=AREA(IFE)*STOR(IFE)*FSS
          DO II=1,3
              IROW=IN2(L2,II)
              DO JJ=1,3
                  ICOL=IN2(L2,JJ)
                  CALL PACK2(IROW,ICOL,AE(II,JJ),1)
                  CALL PACK2(IROW,ICOL,BE(II,JJ),2)
              END DO
          END DO
      END DO
C
C     LINK ELEMENTS.
      L2=NE4
      DO IFL=1,NFL
          L2=L2+1
          CALL WSHAPE3(LEAK(IFL),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
      RETURN
C
C     -------------------------------
      ENTRY WFAULT3(H,HL,SUMLK,SUMFS)
C     -------------------------------
C
C     INITIALIZE OUTPUTS.
      SUMLK=0.0
      IF (NFE.EQ.0) RETURN
      IF (ISTEP.EQ.ISTEP1) WRITE(IWFLT,1091)
C
C     COMPUTE LINK FLUXES.
      L2=NE4
      DO IFL=1,NFL
          L2=L2+1
          II=IN2(L2,1)
          JJ=IN2(L2,2)
          QLK=LEAK(IFL)*(H(II)-H(JJ))
          SUMLK=SUMLK+QLK
      END DO
C
C     COMPUTE STORAGE CHANGE WITHIN FAULTS.
      L2=NE3
      SUMFS=0.0
      DO IFE=1,NFE
          L2=L2+1
          DELH=0.0
          DO J=1,3
              II=IN2(L2,J)
              DELH=DELH+(H(II)-HL(II))/3.0
          END DO
          SUMFS=SUMFS+DELH*AREA(IFE)/DELT
      END DO
      RESID=-SUMFS-SUMLK
C
C     OUTPUT HEADS
      IF (IECHO1.EQ.1) THEN
          WRITE(IWFLT,1092) ISTEP,STIME,SUMLK,SUMFS,RESID
          DONE=.FALSE.
          I1=NN3-4
          DO WHILE (.NOT.DONE)
              I1=I1+5
              I2=I1+4
              IF (I2.GE.NN5) THEN
                  I2=NN5
                  DONE=.TRUE.
              ENDIF
              WRITE(IWFLT,1093) I1,I2,(H(I),I=I1,I2)
          END DO
      ENDIF
C
C     OVERWRITE FILE OF INITIAL HEADS.
      IF(ISAVE.NE.0.AND.ISTEP.EQ.IWHEN)THEN
          ISTART=0
          XHO=0.0
          WRITE(IRFNT,1098)XNFN,ISTART,XHO,ISAVE,IWHEN
                I1=NN3+1
                I2=NN5
          WRITE(IRFNT,1099) (H(I),I=I1,I2)
          CLOSE(IRFNT)
      ENDIF
C
C     OUTPUT LINK FLUXES AND SUM.
      IF (IECHO2.NE.0 .AND. IPRN.EQ.1) THEN
          WRITE(IWFLT,1094) ISTEP,STIME,SUMLK
          L2=NE4
          DO IFL=1,NFL
              L2=L2+1
              II=IN2(L2,1)
              JJ=IN2(L2,2)
              QLK=LEAK(IFL)*(H(II)-H(JJ))
              WRITE(IWFLT,1095) IFL,II,JJ,H(II),H(JJ),QLK
          END DO
      ENDIF
C
C     OUTPUT HEADS TO PLOT FILE.
      IF (IPLT.EQ.1) THEN
          WRITE(IWFLTP,997) ISTEP
          II=NN3
          DO I=1,NFN
              II=II+1
              WRITE(IWFLTP,998) H(II)
          END DO
      ENDIF
C
C     CLOSE OUTPUT FILE.
      IF (ISTEP.EQ.ISTEP2) THEN
          WRITE(IWFLT,1096)
          CLOSE (IWFLT)
      ENDIF
      RETURN
C
C     FORMATS:
1091  FORMAT(1X,'**START SIMULATION IN WFAULT**')
1092  FORMAT(/1X,'COMPUTED FAULT HEADS:'/
     +    1X,'TIME STEP...............................',I10/
     +    1X,'ELAPSED TIME............................',1PE10.2/
     +    1X,'DISCHARHE FROM FAULT....................',1PE10.2/
     +    1X,'STORAGE CHANGE WITHIN FAULT.............',1PE10.2/
     +    1X,'WATER BUDGET RESIDUAL FOR FAULT.........',1PE10.2//
     +    1X,9X,'GLOBAL',9X,'GLOBAL'/
     +    1X,2(11X,'NODE'),5(6X,'NODE'))
1093  FORMAT(6X,I10,5X,I10,5F10.3)
1094  FORMAT(/1X,'FAULT LINK FLUXES:'/
     +    1X,'TIME STEP...............................',I10/
     +    1X,'ELAPSED TIME............................',1PE10.2/
     +    1X,'DISCHARGE FROM FAULT....................',1PE10.2//
     +    1X,14X,'GLOBAL',4X,'GLOBAL'/
     +    1X,15X,'FAULT',6X,'MESH',10X,'FAULT',11X,'MESH'/
     +    1X,6X,'LINK',2(6X,'NODE'),2(11X,'HEAD'),6X,'FLUX')
1095  FORMAT(1X,3I10,2F15.5,1PE10.2)
1096  FORMAT(/1X,'**END SIMULATION IN WFAULT**')
1098  FORMAT(2I10,F10.1,2I10)
1099  FORMAT(20X,5F10.3)
      END SUBROUTINE
