C     -------------------------
      SUBROUTINE PACK1(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 USED TO CREATE COEFFICIENTS FOR THE MATRICES IN
C     THE LINEAR SYSTEM OF EQUATIONS. THE FIRST ENTRY POINT IS CALLED
C     ONCE BYTHE 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     INFILE      NAME OF INPUT FILE.
C     OUTFILE     OUTPUT FILE NAME.
C     OUTEXT      EXTENSION FOR OUTPUT FILE.
C
C     INTERNAL ARRAYS.
      LOGICAL FOUND
      LOGICAL DONE
      CHARACTER*80 INFILE,OUTFILE,OUTEXT
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     IRPAK       INTERNAL READ FOR PAK UNIT NUMBER.
C     IWPAK       INTERNAL WRITE FOR PAK UNIT NUMBER.
C
C     COMMON BLOCKS.
      CHARACTER*80 FNAME(30)
      CHARACTER*12 RNAME1,RNAME2,IJNAME
      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 OR OUTPUT FILE.
      IF (IMATRIX.EQ.1) THEN
          OUTEXT='.PAK'
          CALL OPENFILES(RNAME1,INFILE,OUTFILE,OUTEXT,0,IWPAK)
      ENDIF
      IF (IMATRIX.EQ.2) THEN
          INFILE=FNAME(13)
          CALL OPENFILES(RNAME1,INFILE,OUTFILE,OUTEXT,IRPAK,0)
      ENDIF
C
C     READ MATRIX CONFIGURATION SAVE FROM PREVIOUS SIMULATION.
      IF (IMATRIX.EQ.2) THEN
          READ(IRPAK,910) NIA,NJA
          DO I=1,NIA+1
              READ(IRPAK,912) IA(I)
          END DO
          DO J=1,NJA
              READ(IRPAK,914) JA(J)
          END DO
          WRITE(*,908) NJA
          CLOSE(IRPAK)
          RETURN
      ENDIF
C
C     INITIALIZE IA AND JA ASSUMING ONE NON-ZERO COEFFICIENT PER ROW.
      NIA=NN5+1
      NJA=NN5
      DO II=1,NJA
          IA(II)=II
          JA(II)=II
      END DO
      IA(NIA)=NJA+1
C
C     CONSTRUCT IA AND JA FROM ELEMENT INCIDENCES.
      DO L2=1,NE5
          IF (L2.GT.0 .AND. L2.LE.NE2) THEN
              NNODE=4
          ELSEIF (L2.GT.NE2 .AND. L2.LE.NE3) THEN
              NNODE=2
          ELSEIF (L2.GT.NE3 .AND. L2.LE.NE4) THEN
              NNODE=3
          ELSEIF (L2.GT.NE4 .AND. L2.LE.NE5) THEN
              NNODE=2
          ENDIF
          DO II=1,NNODE
              IROW=IN2(L2,II)
              DO JJ=1,NNODE
                  ICOL=IN2(L2,JJ)
C
C                 IDENTIFY STARTING POSITION OF ROW IN JA.
                  IJA1=IA(IROW)
                  IJA2=IA(IROW+1)-1
C
C                 CHECK FOR COLUMN ALREADY ASSIGNED TO ROW.
                  FOUND=.FALSE.
                  DO IJA=IJA1,IJA2
                      IF(JA(IJA).EQ.ICOL) THEN
                          FOUND=.TRUE.
                          EXIT
                      ENDIF
                  END DO
                  IF (.NOT.FOUND) THEN
C
C                     UPDATE IA AND INSERT NEW COLUMN IN JA.
                      DO I=IROW+1,NIA
                          IA(I)=IA(I)+1
                      END DO
                      DO IJA=NJA,IJA1,-1
                          JA(IJA+1)=JA(IJA)
                      END DO
                      NJA=NJA+1
C
C                     CHECK FOR DINENSIONING.
                      IF (NJA.GT.MAXNON) THEN
                          WRITE(*,904)
                          WRITE(*,906)
                          STOP
                      ENDIF
                      JA(IJA1)=ICOL
                  ENDIF
              END DO
          END DO
      END DO
C
C     SORT COLUMNS WITHIN JA IN ASSENDING ORDER.
      DO I=1,NIA-1
          IJA1=IA(I)
          IJA2=IA(I+1)-1
          DO
              DONE=.TRUE.
              DO IJA=IJA1+1,IJA2
                  IF (JA(IJA).LT.JA(IJA-1)) THEN
                      ITEMP=JA(IJA-1)
                      JA(IJA-1)=JA(IJA)
                      JA(IJA)=ITEMP
                      DONE=.FALSE.
                  ENDIF
              END DO
              IF (DONE) EXIT
          END DO
      END DO
C
C     OUTPUT NUMBER OF NON-ZERO COEFFICIENTS AND CLOSE FILE.
      WRITE(*,908) NJA
C
C     OUTPUT MATRIX CONFIGURATION TO FILE.
      IF (IMATRIX.EQ.1) THEN
          WRITE(IWPAK,910) NIA,NJA
          DO I=1,NIA+1
              WRITE(IWPAK,912) IA(I)
          END DO
          DO J=1,NJA
              WRITE(IWPAK,914) JA(J)
          END DO
          CLOSE(IWPAK)
      ENDIF
      RETURN
C
C     FORMATS:
904   FORMAT(1X,'*ERROR* NUMBER OF NON-ZEROS EXCEEDS DIMENSION MAXNON')
906   FORMAT(1X,'*ERROR* SIMULATION TERMINATED FOR ERROR IN PACK1')
908   FORMAT(1X,'*      NON-ZERO COEFFICIENTS',I10)
910   FORMAT(2I10)
912   FORMAT(I10)
914   FORMAT(I10)
C
C     --------------------------------
      ENTRY PACK2(JROW,JCOL,COEF,IVEC)
C     --------------------------------
C
      IJA1=IA(JROW)
      IJA2=IA(JROW+1)-1
      FOUND=.FALSE.
      DO IJA=IJA1,IJA2
          IF (JCOL.EQ.JA(IJA)) THEN
              IF (IVEC.EQ.1) A(IJA)=A(IJA)+COEF
              IF (IVEC.EQ.2) B(IJA)=B(IJA)+COEF
              FOUND=.TRUE.
              EXIT
          ENDIF
      END DO
      IF (.NOT.FOUND) THEN
          WRITE(*,916)
          WRITE(*,906)
          STOP
      ENDIF
      RETURN
C
C     FORMATS:
916   FORMAT(1X,'*ERROR* POSITION NOT FOUND IN PACK')
      END SUBROUTINE
