C     -----------------
      SUBROUTINE NODES
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 NODAL COORDINATES FOR THE MESH, THE
C     NODAL COLUMNS TO BE USED IN MESH DEFORMATIONS, AND THE PARAMETRES
C     CONTROLLING THE DEFORMATION. THE SUBROUTINE PERFORMS THE
C     COORDINATE ADJUSTMENTS REQUIRED TO DEFORM THE MESH, AND IT
C     OUTPUTS THE ADJUSTMENTS. THE FIRST ENTRY POINT IS CALLED ONCE BY
C     THE MAIN PROGRAM. THE SECOND ENTRY POINT IS CALL ONCE FOR EACH
C     TIMESTEP ITERATION. THE THIRD ENTRY POINT IS CALL ONCE FOR EACH
C     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     XCNAME      COMPARTMENT NAME.
C     H()         COMPUTED HEAD FOR TIMESTEP INTERATION
C     ZC()        DISPLACEMENT OF NODE WITH RESPECT TO THE INTIAL
C                 CONFIGURATION OF MESH.
C     IRXYZ       INTERNAL READ FOR XYZ UNIT NUMBER.
C     IWXYZ       INTERNAL WRITE FOR XYZ UNIT NUMBER.
C     FACX        FACTOR FOR SCALING X-COORDINATE.
C     FACY        FACTOR FOR SCALING Y-COORDINATE.
C     FACZ        FACTOR FOR SCALING Z-COORDINATE.
C     IECHO1      SWITCH FOR ECHO ON COORDINATE INPUTS.
C
      CHARACTER*80 INFILE,OUTFILE,OUTEXT
      CHARACTER*20 XCNAME
      REAL*8 ZC(MAXNN),H(MAXNN)
      LOGICAL DONE
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     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     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     F()         RIGHTHAND VECTOR FOR SYSTEM OF EQUATIONS.
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     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     COL()       COLUMN NODE NUMBER.
C     NLAY()      NUMBER OF NODE LAYERS IN A COLUMN, AS REPRESENTED IN
C                 THEMODEL DOMAIN FROM THE TOP OF THE GRID TO THE BOTTOM
C     ZMIN()      MINIMUM SPACING BETWEEN COLLAPSING NODES.
C     ZMAX()      MAXIMUM ELEVATION OF TOP NODE.
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)
      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 FOR NODES.
      WRITE(*,902)
      OUTEXT='.XYZ.OUT'
      INFILE=FNAME(2)
      CALL OPENFILES(RNAME1,INFILE,OUTFILE,OUTEXT,IRXYZ,IWXYZ)
      WRITE(IWXYZ,903)
      WRITE(IWXYZ,904) INFILE,OUTFILE
C
C     LOOP OVER COMPARTMENTS TO INPUT COORDINATES.
      ILAST=0
      NN=0
      CELL_LOOP1:DO ICELL=1,NCELL
C
C         INPUT COMPARTMENT NODE PARAMETERS.
          READ(IRXYZ,905) XCNAME
          READ(IRXYZ,907) NNC(ICELL),FACX,FACY,FACZ,IECHO1
          WRITE(IWXYZ,908) NCELL,ICELL,XCNAME,NNC(ICELL),
     +        FACX,FACY,FACZ,IECHO1
          NN=NN+NNC(ICELL)
C
C         CHECK COMPARTMENT NAME IN LIST AT EXPECTED LOCATION.
          IF (XCNAME.NE.CNAME(ICELL)) THEN
              WRITE(*,911) ICELL,XCNAME,CNAME(ICELL)
              WRITE(*,912)
              STOP
          ENDIF
C
C         CHECK DIMENSIONING.
          IF (NN.GT.MAXNN) THEN
              WRITE(*,910) ICELL,NN,MAXNN
              WRITE(*,912)
              STOP
          ENDIF
C
C         INPUT COMPARTMENT NODE COORDINATES.
          IF (IECHO1.EQ.1) WRITE(IWXYZ,914)
          I=ILAST
          DO II=1,NNC(ICELL)
              I=I+1
              READ(IRXYZ,915) 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.EQ.1) WRITE(IWXYZ,916) II,I,X(I),Y(I),Z(I)
          END DO
          ILAST=ILAST+NNC(ICELL)
      END DO CELL_LOOP1
C
C     LOOP OVER COMPARTMENTS TO INPUT COLLAPSING MESH.
      ILAST=0
      JLAST=0
      NCOL=0
      CELL_LOOP2:DO ICELL=1,NCELL
C
C         INPUT PARAMETERS FOR COLLAPSING MESH.
          READ(IRXYZ,905) XCNAME
          READ(IRXYZ,930) NC,NMAX,IECHO2
          NCOL=NCOL+NC
C
C         CHECK COMPARTMENT NAME IN LIST AT EXPECTED LOCATION.
          IF (XCNAME.NE.CNAME(ICELL)) THEN
              WRITE(*,911) ICELL,XCNAME,CNAME(ICELL)
              WRITE(*,912)
              STOP
          ENDIF
C
C         CHECK CUMULATIVE NUMBER OF COLUMNS.
          IF (NCOL.GT.MAXCOL) THEN
              WRITE(*,960) ICELL,NCOL,MAXCOL
              WRITE(*,912)
              STOP
          ENDIF
C
C         CHECK NUMBER OF LAYERS.
          IF (NMAX.GT.MAXLAY) THEN
              WRITE(*,962) ICELL,NMAX,MAXLAY
              WRITE(*,912)
              STOP
          ENDIF
C
C         INPUT NODES FOR COLUMNS.
          IF (IECHO2.EQ.1) WRITE(IWXYZ,932) NCELL,ICELL,XCNAME,NC,NMAX
          IF (NCOL.NE.0) THEN
              J=JLAST
              DO ICOL=1,NC
C
C                 INPUT NODE COLUMN.
                  J=J+1
                  READ(IRXYZ,934) ZMIN(J),ZMAX(J),(COL(J,ILAY),
     +                ILAY=1,NMAX)
C
C                 CHECK THAT ZMAX EQUAL OR HIGHER THAN TOP NODE
                  II=COL(J,1)+ILAST
                  IF (ZMAX(J).LT.Z(II)) THEN
                      WRITE(*,966) ICELL,ICOL,II,ZMAX(J),Z(II)
                      WRITE(*,912)
                      STOP
                  ENDIF
C
C                 LOOP OVER LAYERS.
                  DO ILAY=1,NMAX
C
C                     IDENTIFY NUMBER OF NODES IN COLUMN.
                      IF (COL(J,ILAY).EQ.0) THEN
                          NLAY(J)=ILAY-1
                          EXIT
                      ENDIF
                      COL(J,ILAY)=COL(J,ILAY)+ILAST
                      IF (ILAY.EQ.NMAX) THEN
                          NLAY(J)=ILAY
                          EXIT
                      ENDIF
C
C                     ADJUST TO GLOBAL NODE NUMBERING.
                  END DO
C
C                 OUTPUT NODE COLUMN.
                  IF (IECHO2.EQ.1) THEN
                      WRITE(IWXYZ,936) ICOL,ZMIN(J),ZMAX(J),
     +                    (COL(J,ILAY)-ILAST,ILAY=1,NLAY(J))
                      WRITE(IWXYZ,938) ICOL,
     +                    (COL(J,ILAY),ILAY=1,NLAY(J))
                  ENDIF
              END DO
          ENDIF
          ILAST=ILAST+NNC(ICELL)
          JLAST=JLAST+NC
      END DO CELL_LOOP2
C
C     CLOSE FILES AND EXIT NODES.
      NN3=NN
      NN5=NN
      WRITE(IWXYZ,918) MAX(NN,NN3,NN5),MAX(NE,NE2,NE3,NE4,NE5)
      WRITE(*,920)
      CLOSE(IRXYZ)
      IF (NCOL.EQ.0) CLOSE(IWXYZ)
      ICHECK=1
      RETURN
C
C     FORMATS.
902   FORMAT(1X,'*START READING IN NODES')
903   FORMAT(1X,'**START READING IN NODES**')
904   FORMAT(/1X,'FILE NAMES:'/
     +    1X,'INPUT FILE..............................',A80/
     +    1X,'OUTPUT FILE.............................',A80)
905   FORMAT(A20)
907   FORMAT(I10,3F10.0,I10)
908   FORMAT(/1X,'COORDINATES PARAMETERS:'/
     +    1X,'NUMBER OF COMPARTMENTS..................',I10/
     +    1X,'COMPARTMENT NUMBER......................',I10/
     +    1X,'COMPARTMENT NAME........................',A20/
     +    1X,'NUMBER OF NODES IN COMPARTMENT..........',I10/
     +    1X,'SCALING FACTOR FOR X....................',1PE10.2/
     +    1X,'SCALING FACTOR FOR Y....................',1PE10.2/
     +    1X,'SCALING FACTOR FOR Z....................',1PE10.2/
     +    1X,'SWITCH FOR ECHO COORDINATES.............',I10)
910   FORMAT(1X,'*ERROR* NUMBER OF NODES EXCEEDS DIMENSION MAXNN'/
     +    1X,'ICELL=    ',I10/
     +    1X,'NN=       ',I10/
     +    1X,'MAXNN=    ',I10)
911   FORMAT(1X,'*ERROR* COMPARTMENT NAME NOT IN LIST'/
     +    1X,'ICELL=    ',I10/
     +    1X,'XCNAME=   ',A20/
     +    1X,'CNAME=    ',A20)
912   FORMAT(1X,'*ERROR* SIMULATION TERMINATED FOR ERROR IN NODES')
913   FORMAT(A20)
914   FORMAT(/1X,'NODE COORDINATES:'/
     +    1X,5X,'LOCAL',4X,'GLOBAL'/
     +    1X,6X,'NODE',6X,'NODE',11X,'X',11X,'Y',11X,'Z')
915   FORMAT(10X,3F10.0)
916   FORMAT(1X,2I10,3F12.2)
918   FORMAT(/1X,'PROBLEM SIZE:'/
     +    1X,'NUMBER OF NODES.........................',I10/
     +    1X,'NUMNER OF ELEMENTS......................',I10//
     +    1X,'**END READING IN NODES**')
920   FORMAT(1X,'*END   READING IN NODES')
930   FORMAT(3I10)
932   FORMAT(/1X,'COLLAPSING MESH PARAMETERS:'/
     +    1X,'NUMBER OF COMPARTMENTS..................',I10/
     +    1X,'COMPARTMENT NUMBER......................',I10/
     +    1X,'COMPARTMENT NAME........................',A20/
     +    1X,'NUMBER OF NODE COLUMNS.................',I10/
     +    1X,'MAXIMUM NODES IN COLUMN................',I10//
     +    1X,4X,'COLUMN',6X,'ZMIN',6X,'ZMAX',5X,'NODES')
934   FORMAT(2F10.0,20I10)
936   FORMAT(1X,I10,2F10.1,20I10)
938   FORMAT(1X,I10,20X,20I10)
960   FORMAT(1X,'*ERROR* NUMBER OF COLUMNS EXCEEDS DIMENSION'/
     +    1X,'ICELL=    ',I10/
     +    1X,'NCOL=     ',I10/
     +    1X,'MAXCOL=   ',I10)
962   FORMAT(1X,'*ERROR* NUMBER OF NODES IN COLUMN EXCEEDS DIMENSION'
     +    1X,'ICELL=    ',I10/
     +    1X,'NMAX=     ',I10/
     +    1X,'MAXLAY=   ',I10)
966   FORMAT(1X,'*ERROR* ZMAX LOWER THAN TOP NODE'/
     +    1X,'ICELL=    ',I10/
     +    1X,'ICOL=     ',I10/
     +    1X,'II=       ',I10/
     +    1X,'ZMAX=     ',F10.1/
     +    1X,'Z=        ',F10.1)
C
C     ---------------
      ENTRY NODES2(H)
C     ---------------
C
C     BEGIN LOOP OVER COLUMNS.
      IF (NCOL.EQ.0) RETURN
      COL_LOOP:DO ICOL=1,NCOL
C
C         CHECK FOR INVERTED NODES PRIOR TO ADJUSTMENT.
          IF (ICHECK.EQ.1) THEN
              DO ILAY=2,NLAY(ICOL)
                  I1=COL(ICOL,ILAY-1)
                  I2=COL(ICOL,ILAY)
                  DELZ=Z(I1)-Z(I2)
                  IF (DELZ.LT.ZMIN(ICOL)) THEN
                      WRITE(*,968) DELZ,ICOL,I1,I2
                      WRITE(*,912)
                      STOP
                  ENDIF
              END DO
          ENDIF
C
C         IDENTIFY DIRECTION OF GRID ADJUSTMENT.
          I1=COL(ICOL,1)
          IF (H(I1).LT.Z(I1)) THEN
C
C             GRID TO BE CONTRACTED DOWNWARD.
C             ADJUST ALTITUDE OF TOP NODE.
              I1=COL(ICOL,1)
              Z(I1)=H(I1)
C
C             MOVE OTHER NODES DOWN TO MAINTAIN MINIMUM SPACING.
              DO ILAY=2,NLAY(ICOL)-1
                  I1=COL(ICOL,ILAY-1)
                  I2=COL(ICOL,ILAY)
                  DELZ=Z(I1)-Z(I2)
                  IF (DELZ.GE.ZMIN(ICOL)) CYCLE
                  Z(I2)=Z(I1)-ZMIN(ICOL)
              END DO
C
C             CHECK IF NEXT-TO-LAST NODE TOO CLOSE TO BOTTOM NODE.
              I1=COL(ICOL,NLAY(ICOL)-1)
              I2=COL(ICOL,NLAY(ICOL))
              DELZ=Z(I1)-Z(I2)
              IF (DELZ.LT.ZMIN(ICOL)) THEN
C
C                 MOVE ALL NODES BACK UP TO MAINTAIN MINIMUM SPACING.
                  DO ILAY=1,NLAY(ICOL)-1
                      IBACK=NLAY(ICOL)-ILAY
                      I1=COL(ICOL,IBACK)
                      I2=COL(ICOL,IBACK+1)
                      Z(I1)=Z(I2)+ZMIN(ICOL)
                  END DO
              ENDIF
          ELSE
C
C             GRID TO BE EXPANDED UPWARD.
C             ADJUST ALTITUDE OF TOP NODE.
              I1=COL(ICOL,1)
              Z(I1)=H(I1)
              IF (Z(I1).GT.ZMAX(ICOL)) Z(I1)=ZMAX(ICOL)
C
C             CHECK FOR NODES BELOW INITIAL POSITION.
              DO ILAY=2,NLAY(ICOL)-1
                  I1=COL(ICOL,ILAY-1)
                  I2=COL(ICOL,ILAY)
                  IF (Z(I2).LT.Z0(I2)) THEN
                      Z(I2)=Z0(I2)
                      DELZ=Z(I1)-Z(I2)
                      IF (DELZ.LT.ZMIN(ICOL)) Z(I2)=Z(I1)-ZMIN(ICOL)
                  ENDIF
              END DO
          ENDIF
      END DO COL_LOOP
      ICHECK=0
      RETURN
C
C     FORMATS:
968   FORMAT(1X,'*ERROR* NODE ELEVATIONS INVERTED OR TOO CLOSE'/
     +    1X,'DELZ=     ',F10.1/
     +    1X,'ICOL=     ',I10/
     +    1X,'I1=       ',I10/
     +    1X,'I2=       ',I10)
C
C     ------------------
      ENTRY NODES3(ITER)
C     ------------------
C
C     DISPLAY GRID DISPLACEMENTS.
      IF (NCOL.EQ.0) RETURN
      DO I=1,NN5
          ZC(I)=Z(I)-Z0(I)
      END DO
C
C     OUTPUT HEADS TO PRINT FILE.
      IF (IPRN.EQ.1) THEN
          ZCMAX=-1.0E+20
          ZCMIN=1.0E+20
          DO I=1,NN5
              ZCMAX=MAX(ZCMAX,ZC(I))
              ZCMIN=MIN(ZCMIN,ZC(I))
          END DO
          WRITE(IWXYZ,942) ISTEP,ISTEP1,ISTEP2,STIME,ITER,NN3,NN5,
     +        ZCMIN,ZCMAX
          DONE=.FALSE.
          I1=-4
          DO WHILE (.NOT.DONE)
              I1=I1+5
              I2=I1+4
                  IF (I2.GE.NN5) THEN
                      I2=NN5
                      DONE=.TRUE.
                  ENDIF
                  WRITE(IWXYZ,944) I1,I2,(ZC(I),I=I1,I2)
          END DO
      ENDIF
      CALL FLUSH(IWXYZ)
      IF (ISTEP.EQ.ISTEP2) THEN
          CLOSE(IWXYZ)
      ENDIF
      RETURN
C
C     FORMATS:
942   FORMAT(/1X,'COMPUTED NODE DISPLACEMENTS:'/
     +    1X,'STEP NUMBER.............................',I10/
     +    1X,'STARTING STEP NUMBER....................',I10/
     +    1X,'ENDING STEP NUMBER......................',I10/
     +    1X,'ELAPSED TIME............................',1PE10.2/
     +    1X,'NUMBER OF STEP ITERATIONS...............',I10/
     +    1X,'NUMBER FOR LAST NODE IN MESH............',I10/
     +    1X,'NUMBER FPR LAST NODE IN FAULTS..........',I10/
     +    1X,'MINIMUM DISPLACEMENT....................',0PF10.2/
     +    1X,'MAXIMUM DISPLACEMENT....................',0PF10.2//
     +    1X,2(6X,'NODE'),5(3X,'DISPLACEMENT'))
944   FORMAT(1X,2I10,5F15.5)
      RETURN
      END SUBROUTINE
