C     ------------------
      SUBROUTINE WCHEAD1
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 ASSIGNES CONSTANT HEADS TO NODES AND CALCULATES
C     FLOW BASED UPON PARAMETERS AND TRANSIENT TABLES AS DESCRIBED IN
C     THE FEMFLOW3D MANUAL. THE FIRST ENTRY POINT IS CALLED ONCE BY
C     THE MAIN PROGRAM. THE SECOND ENTRY POINT IS CALLED
C     ONCE FOR EACH TIMESTEP ITERATION. THE THIRD ENTRY POINT IS CALLED
C     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     XCNAME      COMPARTMENT NAME.
C     CHNAME()    NAME ASSOCIATED WITH SPECIFIED-HEAD NODE.
C     LEAK()      LEAKANCE FOR SPECIFIED-HEAD NODE.
C     CRH()       RIGHT HANDED LEAKANCE COEFFICIENT.
C     CLH()       LEFT HANDED LEAKANCE COEFFICIENT.
C     CHTAB()     TABLE ASSOCIATED WITH NODE.
C     ITYPE()     SWITCH FOR FLUX CONDITIONS ALLOWED.
C     NPT()       NUMBER OF POINTS IN SPECIFIED-HEAD TABLE.
C     IGET()      SWITCH FOR INTERPOLATION MEATHOD.
C     TTB()       ELAPSED TIME ASSOCIATED WITH POINT IN TABLE.
C     HTB()       TRANSIENT-STATE HEAD ASSOCIATED WITH POINT IN TABLE.
C     XTT()       TIME ENTRIES IN TABLE FOR PASSING TO TABLE FUNCTION.
C     XHT()       SPECIFIED-HEAD ENTRIES IN TABLE FOR PASSING TO TABLE
C                 FUNCTION.
C     H()         HEAD VALUE AT NODE.
C
C     INTERNAL ARRAYS.
      CHARACTER*20 XCNAME
      CHARACTER*80 INFILE,OUTFILE,OUTEXT
      CHARACTER*20 CHNAME(MAXCHN)
      REAL*8 CHBASE(MAXCHN),LEAK(MAXCHN),CRH(MAXCHN),CLH(MAXCHN),
     +    HCHD(MAXCHN)
      INTEGER*4 CHTAB(MAXCHN),ITYPE(MAXCHN)
      INTEGER*4 NPT(MAXTBC),IGET(MAXTBC)
      REAL*8 TTB(MAXTBC,MAXPTC),HTB(MAXTBC,MAXPTC),XTT(MAXPTC),
     +    XHT(MAXPTC)
      REAL*8 H(MAXNN)
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     QCH()       SPECIFIED-HEAD FLUX.
C     CHNODE()    SPECIFIED-HEAD NODE.
C     IRCHD       INTERNAL READ FOR CHD UNIT NUMBER.
C     IWCHD       INTERNAL WRITE FOR CHD UNIT NUMBER.
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 CHNODE(MAXCHN)
      REAL*8 QCH(MAXCHN)
      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
      COMMON /CHEAD/ QCH,CHNODE
      SAVE
C
C     OPEN INPUT AND OUTPUT FILES.
      WRITE(*,900)
      INFILE=FNAME(7)
      OUTEXT='.CHD.OUT'
      CALL OPENFILES(RNAME1,INFILE,OUTFILE,OUTEXT,IRCHD,IWCHD)
      WRITE(IWCHD,902)
      WRITE(IWCHD,903) INFILE,OUTFILE
C
C     LOOP OVER COMPARTMENTS
      ICHN=0
      ITAB=0
      NCHN=0
      NTABLE=0
      ILAST=0
      JLAST=0
      READ(IRCHD,906) IECHO2
      DO ICELL=1,NCELL
C
C         INPUT LOCAL SPECIFIED-HEAD PARAMETERS.
          READ(IRCHD,904) XCNAME
          READ(IRCHD,906) NCHNC,NTABC,IECHO1
          WRITE(IWCHD,908) NCELL,ICELL,XCNAME,NCHNC,NTABC,IECHO1,IECHO2
          NCHN=NCHN+NCHNC
          NTABLE=NTABLE+NTABC
C
C         CHECK COMPARTMENT NAME.
          IF(XCNAME.NE.CNAME(ICELL))THEN
              WRITE(*,910)
              WRITE(*,912)
              STOP
          END IF
C
C         CHECK DIMENSIONING FOR NUMBER OF SPECIFIED HEADS.
          IF (NCHN.GT.MAXCHN) THEN
              WRITE(*,914)
              WRITE(*,912)
              STOP
          ENDIF
C
C         CHECK DIMENSIONING FOR NUMBER OF TABLES
          IF (NTABLE.GT.MAXTBC) THEN
              WRITE(*,916)
              WRITE(*,912)
              STOP
          ENDIF
C
C         INPUT LOCAL PARAMETERS FOR SPECIFIED-HEAD NODES.
          WRITE(IWCHD,918)
          DO II=1,NCHNC
C
C             INPUT NODAL PARAMETERS.
              ICHN=ICHN+1
              READ(IRCHD,920) CHNODE(ICHN),CHBASE(ICHN),LEAK(ICHN),
     +            CHTAB(ICHN),ITYPE(ICHN),CHNAME(ICHN)
              IF(CHTAB(ICHN).NE.0) THEN
                  ITAB=CHTAB(ICHN)+JLAST
              ELSE
                  ITAB=0
              ENDIF
              IF (IECHO1.EQ.1) WRITE(IWCHD,922) CHNODE(ICHN),
     +            CHNODE(ICHN)+ILAST,CHBASE(ICHN),LEAK(ICHN),
     +            CHTAB(ICHN),ITAB,ITYPE(ICHN),CHNAME(ICHN)
              CHNODE(ICHN)=CHNODE(ICHN)+ILAST
              CHTAB(ICHN)=ITAB
C
C             CHECK DIMENSIONING FOR NODAL NUMBERS.
              IF (CHNODE(ICHN).GT.MAXNN) THEN
                  WRITE(*,925)
                  WRITE(*,912)
                  STOP
              ENDIF
C
C             CHECK TABLE NUMBER AGAINST NUMBER OF TABLES.
              IF (CHTAB(ICHN).GT.NTABLE) THEN
                  WRITE(*,926)
                  WRITE(*,912)
                  STOP
              ENDIF
          END DO
C
C         INPUT LOCAL SPECIFIED-HEAD TABLES.
          IF (NTABC.NE.0) THEN
              DO J=1,NTABC
C
C                 INPUT NUMBER OF POINTS IN LOCAL TABLE
                  ITAB=J+JLAST
                  READ(IRCHD,928) NPT(ITAB),IGET(ITAB)
                  IF (IECHO1.NE.0) WRITE(IWCHD,930) J,ITAB,NPT(ITAB),
     +                IGET(ITAB)
C
C                 CHECKING FOR VALID INTERPOLATION METHOD.
                  IF (IGET(ITAB).LT.1 .OR. IGET(ITAB).GT.2) THEN
                      WRITE(*,931)
                      WRITE(*,912)
                      STOP
                  ENDIF
C
C                 CHECK DIMENSIONING FOR NUMBER OF POINTS IN TABLE.
                  IF (NPT(ITAB).GT.MAXPTC) THEN
                      WRITE(*,932)
                      WRITE(*,912)
                      STOP
                  ENDIF
C
C                 INPUT TABLE ENTRIES.
                  DO IPT=1,NPT(ITAB)
                      READ(IRCHD,934) TTB(ITAB,IPT),HTB(ITAB,IPT)
                      IF (IECHO1.NE.0) WRITE(IWCHD,936) IPT,
     +                    TTB(ITAB,IPT),HTB(ITAB,IPT)
                  END DO
              END DO
          ENDIF
          ILAST=ILAST+NNC(ICELL)
          JLAST=NTABLE
      END DO
C
C     CLOSE INPUT FILE AND EXIT WCHEAD
      WRITE(*,938)
      WRITE(IWCHD,940)
      CLOSE (IRCHD)
      RETURN
C
C     FORMATS:
900   FORMAT(1X,'*START READING IN CHEAD')
902   FORMAT(1X,'**START READING IN CHEAD**')
903   FORMAT(/1X,'FILE NAMES:'/
     +    1X,'INPUT FILE..............................',A80/
     +    1X,'OUTPUT FILE.............................',A80)
904   FORMAT(A20)
906   FORMAT(4I10)
908   FORMAT(/1X,'SPECIFIED-HEAD PARAMETERS:'/
     +    1X,'NUMBER OF COMPARTMENTS..................',I10/
     +    1X,'COMPARTMENT NUMBER......................',I10/
     +    1X,'COMPARTMENT NAME........................',A20/
     +    1X,'NUMBER OF NODES.........................',I10/
     +    1X,'NUMBER OF TABLES........................',I10/
     +    1X,'SWITCH FOR ECHO OF INPUTS...............',I10/
     +    1X,'SWITCH FOR ECHO OF STEP FLUXES..........',I10)
910   FORMAT(1X,'*ERROR* COMPARTMENT NAME NOT IN LIST')
912   FORMAT(1X,'*ERROR* RUN TERMINATED FOR ERROR IN CHEAD')
914   FORMAT(1X,'*ERROR* NUMBER OF NODES EXCEEDS DIMENSION')
916   FORMAT(1X,'*ERROR* NUMBER OF TABLES EXCEEDS DIMENSION')
918   FORMAT(/1X,'CONSTANT-HEAD NODES:'/
     +    1X,5X,'LOCAL',4X,'GLOBAL',25X,'LOCAL',4X,'GLOBAL'/
     +    1X,6X,'NODE',6X,'NODE',6X,'HEAD',2X,'LEAKANCE',5X,'TABLE',
     +    5X,'TABLE',6X,'TYPE',1X,'NAME')
920   FORMAT(I10,2F10.0,2I10,1X,A20)
922   FORMAT(1X,2I10,F10.2,1PE10.2,3I10,1X,A20)
925   FORMAT(1X,'*ERROR* NODE NUMBER EXCEEDS DIMENSION')
926   FORMAT(1X,'*ERROR* TABLE NUMBER EXCEEDS DIMENSION')
928   FORMAT(2I10)
930   FORMAT(/1X,'CONSTANT-HEAD TABLE:'/
     +    1X,'LOCAL TABLE NUMBER......................',I10/
     +    1X,'GLOBAL TABLE NUMBER.....................',I10/
     +    1X,'NUMBER OF ENTRIES IN TABLE..............',I10/
     +    1X,'SWITCH FOR INTERPOLATION METHOD.........',I10//
     +    1X,'TABLE ENTRIES:'/
     +    1X,5X,'POINT',6X,'TIME',6X,'HEAD')
931   FORMAT(1X,'*ERROE* INTERPOLATION MEHTOD NOT VALID')
932   FORMAT(1X,'*ERROR* NUMBER OF POINTS EXCEEDS DIMENSION MAXPT')
934   FORMAT(2F10.0)
936   FORMAT(1X,I10,2F10.2)
938   FORMAT(1X,'*END   READING IN CHEAD')
940   FORMAT(/1X,'**END READING IN CHEAD**')
C
C     ----------------
      ENTRY WCHEAD2(H)
C     ----------------
C
C     COMPUTE LEFT-HAND AND RIGHT-HAND COEFFICIENTS.
      IF (NCHN.EQ.0) RETURN
      DO ICHN=1,NCHN
C
C         IDENTIFY NODE AND INTERPOLATE FROM TABLE.
          II=CHNODE(ICHN)
          ITAB=CHTAB(ICHN)
          IF (ITAB.NE.0) THEN
              DO IPT=1,NPT(ITAB)
                  XTT(IPT)=TTB(ITAB,IPT)
                  XHT(IPT)=HTB(ITAB,IPT)
              END DO
          ENDIF
C
C         ADD TEMPORAL COMPONENT OF SPECIFIED HEAD.
          IF (ITAB.NE.0) THEN
              HCHD(ICHN)=CHBASE(ICHN)+
     +            TABLE(NPT(ITAB),XTT,XHT,STIME-DELT,STIME,IGET(ITAB))
          ELSE
              HCHD(ICHN)=CHBASE(ICHN)
          ENDIF
          IF (ITYPE(ICHN).EQ.1) THEN
C
C             FLOW POSSIBLE IN TWO DIRECTIONS.
              CLH(ICHN)=LEAK(ICHN)
              CRH(ICHN)=LEAK(ICHN)*HCHD(ICHN)
          ELSE
C
C             FLOW POSSIBLE OUTWARD ONLY.
              IF (H(II).GT.HCHD(ICHN)) THEN
C
C                 FLOW OCCURS.
                  CLH(ICHN)=LEAK(ICHN)
                  CRH(ICHN)=LEAK(ICHN)*HCHD(ICHN)
              ELSE
C
C                 NO FLOW OCCURS.
                  CLH(ICHN)=0.0
                  CRH(ICHN)=0.0
              ENDIF
          ENDIF
C
C         PACK COEFFICIENTS INTO LEFTHAND MATRIX AND RIGHTHAND VECTOR.
          CALL PACK2(II,II,CLH(ICHN),1)
          F(II)=F(II)+CRH(ICHN)
      END DO
      RETURN
C
C     ----------------------
      ENTRY WCHEAD3(H,SUMCH)
C     ----------------------
C
C     INITIALIZE OUTPUTS..
      SUMCH=0.0
      IF (ISTEP.EQ.ISTEP1) WRITE(IWCHD,942)
C
C     COMPUTE NODAL FLUX AND SUM.
      DO ICHN=1,NCHN
          II=CHNODE(ICHN)
          QCH(ICHN)=CRH(ICHN)-CLH(ICHN)*H(II)
          SUMCH=SUMCH+QCH(ICHN)
      END DO
C
C     OUTPUT NODAL FLUX AND SUM.
      IF (IECHO2.EQ.1 .AND. IPRN.EQ.1) THEN
          WRITE(IWCHD,944) ISTEP,STIME,SUMCH
          DO ICHN=1,NCHN
              II=CHNODE(ICHN)
              WRITE(IWCHD,946) ICHN,CHNODE(ICHN),HCHD(ICHN),H(II),
     +            QCH(ICHN),ITYPE(ICHN),CHNAME(ICHN)
          END DO
      END IF
      IF (ISTEP.EQ.ISTEP2) THEN
          WRITE(IWCHD,948)
          CLOSE (IWCHD)
      ENDIF
      RETURN
C
C     FORMATS.
942   FORMAT(1X,'**START SIMULATION IN CHEAD**')
944   FORMAT(/1X,'SPECIFIED-HEAD NODAL FLUXES:'/
     +    1X,'TIME STEP...............................',I10/
     +    1X,'SIMULATION TIME.........................',1PE10.2/
     +    1X,'NET FLUX................................',1PE10.2//
     +    1X,5X,'POINT',6X,'NODE',5X,'CHEAD',6X,'HEAD',6X,'FLUX',
     +    6X,'TYPE',1X,'NAME')
946   FORMAT(1X,2I10,2F10.2,1PE10.2,I10,1X,A20)
948   FORMAT(/1X,'**END SIMULATION IN CHEAD**')
      END SUBROUTINE
