C     ----------------
      SUBROUTINE ELEMS
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 INPUTS THE ELEMENT INCIDENCES FOR THE WEDGE
C     ELEMENTS, CREATES THE ELEMENT INCIDENCES FOR THE TETRAHEDRAL
C     ELEMENTS, AND OUTPUTS A .3DM FILE FOR USE WITH GMS. THE .3DM
C     FILE REPRESENTS THE AGGREGATION OF THE INDIVIDUAL COMPARTMENT
C     MESHES. THE .3DM FILE CAN BE USED IN GMS TO VIEW THE MESH AS
C     A THREE-DIMENSIONAL OBJECT WITH VARIOUS OPTIONS. THE
C     SUBROUTINE IS CALLED ONCE BY THE MAIN PROGRAM.
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     ITET1()     TETRAHEDRAL NODAL ORDER FOR WEDGE WITH ONE
C                 TETRAHEDRON.
C     ITET2()     TETRAHEDRAL NODAL ORDER FOR WEDGE WITH TWO
C                 TETRAHEDRONS.
C     ITET3()     TETRAHEDRAL NODAL ORDER FOR WEDGE WITH THREE
C                 TETRAHEDRONS.
C     MAT()       MATERIAL IDENTIFIER FOR ELEMENT MATERIAL.
C     IRELE       INTERNAL READ FOR ELE UNIT NUMBER.
C     IWELE       INTERNAL WRITE FOR ELE UNIT NUMBER.
C     IMESH       SWITCH FOR OUTPUT OF GMS .3DM FILE
C     IECHO1      SWITCH FOR OUTPUT OF WEDGE ELEMENTS.
C     IECHO2      SWITCH FOR OUTPUT OF TETRAHEDRAL ELEMENTS.
C
      CHARACTER*80 INFILE,OUTFILE,MESHFILE,OUTEXT
      CHARACTER*20 XCNAME
      INTEGER*4 ITET1(4),ITET2(2,4),ITET3(3,4),MAT(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     X()         X-COORDINATE OF A NODE.
C     Y()         Y-COORDINATE OF A NODE.
C     Z()         Z-COORDINATE OF A NODE.
C     IEL()       PARENT WEDGE ELEMENT FOR TETRAHEDRAL ELEMENT.
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 IN THE X-DIRECTION.
C     QY()        DARCY GROUNDWATER VELOCITY IN THE Y-DIRECTION.
C     QZ()        DARCY GROUNDWATER VELOCITY IN THE Z-DIRECTION.
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)
      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     OPEN INPUT AND OUTPUT FILES.
      WRITE(*,900)
      OUTEXT='.ELE.OUT'
      INFILE=FNAME(3)
      CALL OPENFILES(RNAME1,INFILE,OUTFILE,OUTEXT,IRELE,IWELE)
      WRITE(IWELE,902)
      WRITE(IWELE,904) INFILE,OUTFILE
      READ(IRELE,906) IMESH
      IF (IMESH.NE.0) THEN
          OUTEXT='.ELE.3DM'
          CALL OPENFILES(RNAME1,INFILE,MESHFILE,OUTEXT,0,IW3DM)
          WRITE(IWELE,908) MESHFILE
      ENDIF
C
C     LOOP OVER COMPARTMENTS
      NE=0
      L=0
      ILAST=0
      CELL_LOOP:DO ICELL=1,NCELL
C
C         INPUT COMPARTMENT ELEMENT PARAMETERS
          READ(IRELE,910) XCNAME
          READ(IRELE,912) NEC(ICELL),IECHO1,IECHO2
          WRITE(IWELE,914) NCELL,ICELL,XCNAME,NEC(ICELL),IECHO1,IECHO2,
     +        IMESH
          NE=NE+NEC(ICELL)
C
C         CHECK COMPARTMENT NAME.
          IF (XCNAME.NE.CNAME(ICELL)) THEN
              WRITE(*,916)
              WRITE(*,918)
              STOP
          ENDIF
C
C         CHECK DIMENSIONING.
          IF (NE.GT.MAXNE) THEN
              WRITE(*,920)
              WRITE(*,918)
              STOP
          ENDIF
C
C         INPUT LOCAL ELEMENT INCIDENCES FOR WEDGE ELEMENT.
          IF (IECHO1.NE.0) WRITE(IWELE,922)
          DO LL=1,NEC(ICELL)
              L=L+1
              IF (LL.EQ.1) I1=L
              IF (LL.EQ.NEC(ICELL)) I2=L
              READ(IRELE,924) (IN(L,J),J=1,6),MAT(L)
C
C             OUTPUT LOCAL ELEMENT INCIDENCES.
              IF (IECHO1.NE.0) WRITE(IWELE,926) LL,(IN(L,J),J=1,6),
     +            MAT(L)
              DO J=1,6
                  IF(IN(L,J).EQ.0) EXIT
                  IN(L,J)=IN(L,J)+ILAST
C
C                 CHECK DIMENSIONING.
                  IF (IN(L,J).GT.NN) THEN
                      WRITE(*,928)
                      WRITE(*,918)
                      STOP
                  ENDIF
              END DO
          END DO
          IF (IECHO1.NE.0) THEN
              WRITE(IWELE,930)
              DO II=I1,I2
                  WRITE(IWELE,932) II,(IN(II,J),J=1,6),MAT(II)
              END DO
          ENDIF
          ILAST=ILAST+NNC(ICELL)
      END DO CELL_LOOP
C
C     DEFINE SUBDIVIION OF GLOBAL WEDGE ELEMENT INTO TETRAHEDRONS.
      ITET3(1,1)=1
      ITET3(1,2)=2
      ITET3(1,3)=6
      ITET3(1,4)=3
      ITET3(2,1)=5
      ITET3(2,2)=6
      ITET3(2,3)=1
      ITET3(2,4)=2
      ITET3(3,1)=4
      ITET3(3,2)=5
      ITET3(3,3)=6
      ITET3(3,4)=1
      ITET2(1,1)=1
      ITET2(1,2)=5
      ITET2(1,3)=3
      ITET2(1,4)=2
      ITET2(2,1)=4
      ITET2(2,2)=5
      ITET2(2,3)=3
      ITET2(2,4)=1
      ITET1(1)=4
      ITET1(2)=2
      ITET1(3)=3
      ITET1(4)=1
C
C     GENERATE INCIDENCES FOR TETRAHEDRAL ELEMENTS.
      L2=0
      DO L=1,NE
C
C         THREE TETRAHEDRONS IN WEDGE.
          IF(IN(L,6).NE.0) THEN
              DO LT=1,3
                  L2=L2+1
                  IEL(L2)=L
                  DO I=1,4
                      II=ITET3(LT,I)
                      IN2(L2,I)=IN(L,II)
                  END DO
              END DO
              CYCLE
          ENDIF
C
C         TWO TETRAHEDRONS IN WEDGE.
          IF(IN(L,6).EQ.0.AND.IN(L,5).NE.0) THEN
              DO LT=1,2
                  L2=L2+1
                  IEL(L2)=L
                  DO I=1,4
                      II=ITET2(LT,I)
                      IN2(L2,I)=IN(L,II)
                  END DO
              END DO
              CYCLE
          ENDIF
C
C         ONE TETRAHEDRON IN WEDGE.
          IF(IN(L,6).EQ.0.AND.IN(L,5).EQ.0) THEN
              L2=L2+1
              IEL(L2)=L
              DO I=1,4
                  II=ITET1(I)
                  IN2(L2,I)=IN(L,II)
              END DO
              CYCLE
          ENDIF
      END DO
      NE2=L2
      IF (IECHO2.EQ.1) THEN
          WRITE(IWELE,936)
          DO L2=1,NE2
              WRITE(IWELE,938) L2,IEL(L2),(IN2(L2,J),J=1,4),
     +            MAT(IEL(L2))
          END DO
      ENDIF
C
C     OUTPUT 3DM FILE.
      IF (IMESH.NE.0) THEN
          WRITE(IW3DM,940)
          DO L=1,NE
C
C             IDENTIFY ELEMENT TYPE.
              IF (IN(L,6).EQ.0) THEN
C
C                 PINCHOUT ELEMENT (TETRAHEDRON OR PYRAMID).
                  IF (IN(L,5).EQ.0) THEN
C
C                     TWO ZERO-HEIGHT EDGES (TETRAHEDRON).
                      WRITE(IW3DM,942) L,IN(L,3),IN(L,1),IN(L,4),
     +                    IN(L,2),MAT(L)
                  ELSE
C
C                     ONE ZERO-HEIGHT EDGE (PYRAMID).
                      WRITE(IW3DM,944) L,IN(L,1),IN(L,2),IN(L,5),
     +                    IN(L,4),IN(L,3),MAT(L)
                  ENDIF
              ELSE
C
C                 NOT PICHOUT ELEMENT (WEDGE).
                  WRITE(IW3DM,946) L,IN(L,4),IN(L,5),IN(L,6),
     +                IN(L,1),IN(L,2),IN(L,3),MAT(L)
              ENDIF
          END DO
          DO I=1,NN
              WRITE(IW3DM,948) I,X(I),Y(I),Z(I)
          END DO
      ENDIF
C
C     CLOSE INPUT FILE AND EXIT ELEMS
      NE3=NE2
      NE4=NE2
      NE5=NE2
      WRITE(IWELE,950) MAX(NN,NN3,NN5),MAX(NE,NE2,NE3,NE4,NE5)
      WRITE(*,952)
      CLOSE(IRELE)
      CLOSE(IWELE)
      IF (IMESH.NE.0) CLOSE (IW3DM)
      RETURN
C
C     FORMATS
900   FORMAT(1X,'*START READING IN ELEMS')
902   FORMAT(1X,'**START READING IN ELEMS**')
904   FORMAT(/1X,'FILE NAMES:'/
     +    1X,'INPUT FILE..............................',A80/
     +    1X,'OUTPUT FILE.............................',A80)
905   FORMAT(1X,'MATRIX OUTPUT FILE......................',A80)
906   FORMAT(I10)
907   FORMAT(1X,'MATRIX INPUT FILE.......................',A80)
908   FORMAT(1X,'MESH 3DM FILE...........................',A80)
910   FORMAT(A20)
912   FORMAT(3I10)
914   FORMAT(/1X,'ELEMENT PARAMETERS:'/
     +    1X,'NUMBER OF COMPARTMENTS..................',I10/
     +    1X,'COMPARTMENT NUMBER......................',I10/
     +    1X,'COMPARTMENT NAME........................',A20/
     +    1X,'NUMBER OF WEDGE ELEMENTS................',I10/
     +    1X,'SWITCH FOR ECHO OF WEDGE NODES..........',I10/
     +    1X,'SWITCH FOR ECHO OF TETRAHEDRAL NODES....',I10/
     +    1X,'SWITCH FOR OUTPUT OF MESH FILE (3DM)....',I10)
916   FORMAT(1X,'*ERROR* COMPARTMENT NAME NOT IN LIST')
918   FORMAT(1X,'*ERROR* SIMULATION TERMINATED FOR ERROR IN ELEMS')
920   FORMAT(1X,'*ERROR* NUMBER OF ELEMENTS EXCEEDS DIMENSION MAXNE')
922   FORMAT(/1X,'WEDGE ELEMENT LOCAL INCIDENCES:'//
     +    1X,6X,'ELEM',6(6X,'NODE'),2X,'MATERIAL')
924   FORMAT(10X,7I10)
926   FORMAT(1X,8I10)
928   FORMAT(1X,'*ERROR* NODE NUMBER EXCEEDS DIMENSION MAXNN')
930   FORMAT(/1X,'WEDGE ELEMENT GLOBAL INCIDENCES:'//
     +    1X,6X,'ELEM',6(6X,'NODE'),2X,'MATERIAL')
932   FORMAT(1X,8I10)
936   FORMAT(/1X,'TETRAHEDRAL ELEMENT GLOBAL INCIDENCES:'//
     +    1X,5X,'TETRA',5X,'WEDGE'/
     +    1X,2(6X,'ELEM'),4(6X,'NODE'),2X,'MATERIAL')
938   FORMAT(1X,7I10)
940   FORMAT('MESH3D')
942   FORMAT('E4T',7X,5I10,20X,I10)
944   FORMAT('E5P',7X,6I10,10X,I10)
946   FORMAT('E6W',7X,8I10)
948   FORMAT('ND',8X,I10,3F20.4)
950   FORMAT(/1X,'PROBLEM SIZE:'/
     +    1X,'NUMBER OF NODES.........................',I10/
     +    1X,'NUMNER OF ELEMENTS......................',I10//
     +    1X,'**END READING IN ELEMS**')
952   FORMAT(1X,'*END   READING IN ELEMS')
      END SUBROUTINE
