C     ------------------
      SUBROUTINE WRIVER1
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 SIMULATE RIVER-AQUIFER INTERATIONS. THE
C     SIMULATION IS PERFORMED BY ROUTING INFLOWS THROUGH RIVER NETWORKS
C     THAT CONSIST OF A MAIN RIVER CHANNEL AND TRIBUTARY CHANNELS, WHILE
C     ACCOUNTING FOR RIVER FLOW DEPLETIONS AND ACCRETIONS OWING TO THE
C     EXCHANGE OF WATER BETWEEN THE RIVER AND THE 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     OUTFILE1    OUTPUT OF PLOT FILE.
C     OUTFILE2    OUTPUT OF GRAPH FILE.
C     OUTEXT      EXTENSION FOR OUTPUT FILE.
C     XCNAME      COMPARTMENT NAME.
C     IECHO2      SWITCH FOR ECHO OF STEP FLUXES.
C     SNAME()     NAME ASSOCIATED WITH RIVER REACH.
C     QNAME()     NAME ASSOCIATED WITH INFLOW.
C     GNAME()     NAME ASSOCIATED WITH HYDROGRAPH SITE.
C     NRN()       NUMBER OF LINKS IN RIVER REACH.
C     NRC()       NUMBER OF RIVER REACHES IN COMPARTMENT.
C     NTBXC       NUMBER OF CROSS-SECTION TABLE ASSOCIATED WITH
C                 COMPARTMENT.
C     NTBQC       NUMBER OF INFLOW TABLES ASSOCIATED WITH COMPARTMENT.
C     NGPHC       NUMBER OF HYDROGRAPH SITES ASSOCIATED WITH COMPARTMENT
C     ORDER()     ORDER OR REACHES.
C     LAST()      FLAG TO INDICATE A DOWNSTREAM REACH IN NETWORK.
C     FACK        SCALING FACTOR FOR HYDRAULIC CONDUCTIVITY OF CHANNEL
C                 BET FOR LINKS WITHIN REACH.
C     LENGTH()    CHANNEL LENGTH FOR LINK.
C     TABQ()      INDICATE DISCHARGE TABLE NUMBER.
C     TABQC()     INFLOW TABLE ASSIGNED TO LINK.
C     JOIN()      RIVER REACH NUMBER JOINING AT LINK.
C     CJOIN()     COMPARTMENT NAME FOR RIVER REACH JOININ AT LINK.
C     ICASE()     INDICATE CONNECTION TYPE BETWEEN STREAM AND AQUIFER.
C     TABQC       INDICATE DISCHARGE TABLE NUMBER.
C     CJOIN()     INDICATE COMPARTMENT FOR TRIBUTARY REACH.
C     ZBED()      ELEVATION OF CHANNEL BED AT NODE.
C     BBED()      THICKNESS OF CHANNEL BED FOR LINK.
C     KBED()      VERTICAL HYDRAULIC CONDUCTIVITY OF CHANNEL BED FOR
C                 LINK.
C     LEAK        CHANNEL LEAKANCE.
C     RNODE()     STREAM NODE.
C     RNODEC()    NODE NUMBER FOR RIVER REACH LINK.
C     TABX        INDICATES CROSS-SECTION TABLE.
C     TABXC()     CROSS-SECTION TABLE ASSIGNED TO LINK.
C     QXTB()      DISCHARGE ASSOCIATED WITH POINT IN CROSS-SECTION TABLE
C                 ASSIGNED RECHARGE IS ALWAYS POSITIVE.
C     WXTB()      STREAMFLOW WIDTH, CORRESPONDING TO DISCHARGE,
C                 ASSOCIATED WITH POINT IN TABLE.
C     DXTB()      STREAMFLOW AVERAGE DEPTH, CORRESPONDING TO DISCHARGE,
C                 ASSOCIATED WITH POINT IN TABLE.
C     TTB()       ELAPSED TIME ASSOCIATED WITH POINT IN FLOW TABLE.
C     QTB()       INFLOW ASSOCIATED WITH POINT IN TABLE. INFLOW TO
C                 CHANNEL IS POSITIVE AND DIVERSION FROM CHANNEL IS
C                 NEGATIVE.
C     NPTX()      NUMBER OF POINTS IN CROSS-SECTION TABLE.
C     NPTQ()      NUMBER OF POINTS IN INFLOW TABLE.
C     QFACT       SCALING FACTOR FOR DISCHARGES IN INFLOW TABLE.
C     IGET()      SWITCH FOR INTERPOLATION METHOD.
C     GREACH()    RIVER REACH ASSOCIATED WITH HYDROGRAPH SITE.
C     GLINK()     LINK WITH REACH ASSOCIATED WITH HYDROGRAPHY SITE.
C     GREACHC     STREAMFLOW FOR HYDROGRAPHY POINT.
C     GTIME()     TIME FOR HYDROGRAPHY POINT.
C     QGRAPH()    STREAMFLOW AT HYDROGRAPHY POINT.
C     Q1()        UPSTREAM LINK INFLOW.
C     Q2()        TRIBUTARY LINK INFLOW.
C     Q3()        ASSIGNED LINK INFLOW.
C     Q4()        RECHARGE FROM LINK.
C     Q5()        DOWNSTREAM LINK OUTFLOW.
C     XTTB()      DISCHARGE TABLE FOR PASSING TO TABLE FUNCTION.
C     XQTB()      DISCHARGE TABLE FOR PASSING TO TABLE FUNCTION.
C     XQXTB()     CROSS-SECTION TABLE FOR PASSING TO TABLE FUNCTION.
C     XWXTB()     CROSS-SECTION TABLE FOR PASSING TO TABLE FUNCTION.
C     XDXTB()     CROSS-SECTION TABLE FOR PASSING TO TABLE FUNCTION.
C     H()         HEAD VALUE AT NODE.
C     CLH()       LEFT HANDED LEAKANCE COEFFICIENT.
C     CRH()       RIGHT HANDED LEAKANCE COEFFICIENT.
C     IECHO1      SWITCH FOR ECHO OF RIVER INPUTS.
C     IECHO2      SWITCH FOR OUTPUT OF RIVER DISCHARGES AND STREAM
C                 AQUIFER EXCHANGES.
C
C     INTERNAL ARRAYS.
      CHARACTER*80 INFILE,OUTFILE1,OUTFILE2,OUTEXT
      CHARACTER*20 SNAME(MAXRCH),QNAME(MAXTBQ),GNAME(MAXSWG),XCNAME
      INTEGER*4 NRN(MAXRCH),NRC(MAXCEL),ORDER(MAXRCH),LAST(MAXRCH)
      REAL*8 LENGTH(MAXRCH,MAXLNK)
      INTEGER*4 TABQ(MAXRCH,MAXLNK),JOIN(MAXRCH,MAXLNK),
     +    ICASE(MAXRCH,MAXLNK),TABQC
      CHARACTER*20 CJOIN(MAXRCH,MAXLNK)
      REAL*8 ZBED(MAXRCH,MAXLNK),BBED(MAXRCH,MAXLNK),
     +    KBED(MAXRCH,MAXLNK),LEAK
      INTEGER*4 TABX(MAXRCH,MAXLNK),RNODE(MAXRCH,MAXLNK),RNODEC,TABXC
      REAL*8 QXTB(MAXTBX,MAXPTX),WXTB(MAXTBX,MAXPTX),
     +    DXTB(MAXTBX,MAXPTX),TTB(MAXTBQ,MAXPTQ),QTB(MAXTBQ,MAXPTQ)
      INTEGER*4 NPTX(MAXTBX),NPTQ(MAXTBQ),IGET(MAXTBQ)
      INTEGER*4 GREACH(MAXSWG),GLINK(MAXSWG),GREACHC
      REAL*8 GTIME(MAXSTP),QGRAPH(MAXSWG,MAXSTP)
      LOGICAL FOUND,DONE
      REAL*8 Q1(MAXRCH,MAXLNK),Q2(MAXRCH,MAXLNK),Q3(MAXRCH,MAXLNK),
     +    Q4(MAXRCH,MAXLNK),Q5(MAXRCH,MAXLNK)
      REAL*8 XTTB(MAXPTQ),XQTB(MAXPTQ),XQXTB(MAXPTX),XWXTB(MAXTBX),
     +    XDXTB(MAXTBX)
      REAL*8 H(MAXNN),CLH(MAXRCH,MAXLNK),CRH(MAXRCH,MAXLNK)
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     IRRIV       INTERNAL READ FOR RIV UNIT NUMBER.
C     IWRIV       INTERNAL WRITE FOR RIV 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)
      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='.RIV.OUT'
      INFILE=FNAME(10)
      CALL OPENFILES(RNAME1,INFILE,OUTFILE1,OUTEXT,IRRIV,IWRIV)
      OUTEXT='.RIV.GPH'
      CALL OPENFILES(RNAME1,INFILE,OUTFILE2,OUTEXT,0,IWRIVG)
      WRITE(IWRIV,902)
      WRITE(IWRIV,904) INFILE,OUTFILE1,OUTFILE2
C
C     LOOP OVER COMPARTMENTS.
      ILAST=0
      JLAST=0
      KLAST=0
      LLAST=0
      IRCH=0
      ITBX=0
      ITBQ=0
      IGPH=0
      NRCH=0
      NTBX=0
      NTBQ=0
      NJOIN=0
      NGPH=0
      READ(IRRIV,908) IECHO2
      CELL_LOOP:DO ICELL=1,NCELL
C
C         INPUT NUMBER OF REACHES, SECTION TABLES, AND INFLOW TABLES.
          READ(IRRIV,906) XCNAME
          READ(IRRIV,908) NRC(ICELL),NTBXC,NTBQC,NGPHC,IECHO1
          NTBQ=NTBQ+NTBQC
          WRITE(IWRIV,910) NCELL,ICELL,XCNAME,NRC(ICELL),NTBXC,NTBQC,
     +        NGPHC,IECHO1,IECHO2
          IF (NRC(ICELL).EQ.0) THEN
              ILAST=ILAST+NNC(ICELL)
              CYCLE
          END IF
          NRCH=NRCH+NRC(ICELL)
          NTBX=NTBX+NTBXC
          NTBQ=NTBQ+NTBQC
          NGPH=NGPH+NGPHC
C
C         CHECK COMPARTMENT NAME.
          IF (XCNAME.NE.CNAME(ICELL)) THEN
              WRITE(*,912)
              WRITE(*,914)
              STOP
          ENDIF
C
C         CHECK NUMBER OF REACHES
          IF (NRCH.GT.MAXRCH) THEN
              WRITE(*,916)
              WRITE(*,914)
              STOP
          ENDIF
C
C         CHECK NUMBER OF CROSS SECTION TABLES.
          IF (NTBX.GT.MAXTBX) THEN
              WRITE(*,918)
              WRITE(*,914)
              STOP
          ENDIF
C
C         CHECK NUMBER OF INFLOW TABLES.
          IF (NTBQ.GT.MAXTBQ) THEN
              WRITE(*,920)
              WRITE(*,914)
              STOP
          ENDIF
C
C         INPUT RIVER REACH PARAMETERS.
          REACH_LOOP1:DO II=1,NRC(ICELL)
              IRCH=IRCH+1
              READ(IRRIV,922) NRN(IRCH),LAST(IRCH),FACK,
     +            SNAME(IRCH)
              IF (IECHO1.NE.0) WRITE(IWRIV,924) II,IRCH,NRN(IRCH),
     +            LAST(IRCH),SNAME(IRCH),FACK
C
C             CHECK NUMERER OF LINKS PER REACH.
              IF (NRN(IRCH).GT.MAXLNK) THEN
                  WRITE(*,925)
                  WRITE(*,914)
                  STOP
              ENDIF
C
C             LOOP OVER LINKS.
              LINK_LOOP1:DO ILNK=1,NRN(IRCH)
                  READ(IRRIV,926) RNODEC,LENGTH(IRCH,ILNK),
     +                ZBED(IRCH,ILNK),KBED(IRCH,ILNK),BBED(IRCH,ILNK),
     +                TABXC,TABQC,JOIN(IRCH,ILNK),
     +                CJOIN(IRCH,ILNK)
                  KBED(IRCH,ILNK)=KBED(IRCH,ILNK)*FACK
                  RNODE(IRCH,ILNK)=RNODEC+ILAST
                  TABX(IRCH,ILNK)=TABXC+JLAST
                  IF (TABQC.NE.0) THEN
                      TABQ(IRCH,ILNK)=TABQC+KLAST
                  ELSE
                      TABQ(IRCH,ILNK)=0
                  ENDIF
                  IF (JOIN(IRCH,ILNK).NE.0) THEN
                      NJOIN=NJOIN+1
                  ELSE
                      CJOIN(IRCH,ILNK)='None'
                  ENDIF
                  IF (IECHO1.NE.0) WRITE(IWRIV,928) ILNK,RNODEC,
     +                RNODE(IRCH,ILNK),LENGTH(IRCH,ILNK),
     +                ZBED(IRCH,ILNK),KBED(IRCH,ILNK),BBED(IRCH,ILNK),
     +                TABXC,TABX(IRCH,ILNK),TABQC,TABQ(IRCH,ILNK),
     +                JOIN(IRCH,ILNK),CJOIN(IRCH,ILNK)
C
C                 CHECK JOINING COMPARTMENT.
                  IF (JOIN(IRCH,ILNK).NE.0) THEN
                      FOUND=.FALSE.
                      DO JCELL=1,NCELL
                          IF (CJOIN(IRCH,ILNK).EQ.CNAME(JCELL)) THEN
                             FOUND=.TRUE.
                             EXIT
                          ENDIF
                      END DO
                      IF (.NOT.FOUND) THEN
                          WRITE(*,929)
                          WRITE(*,914)
                          STOP
                      ENDIF
                  ENDIF
C
C                 CHECK SECTION TABLE WITHIN RANGE.
                  IF (TABXC.GT.NTBXC) THEN
                      WRITE(*,932)
                      WRITE(*,914)
                      STOP
                  ENDIF
C
C                 CHECK INFLOW TABLE WITHIN RANGE.
                  IF (TABQC.GT.NTBQC) THEN
                      WRITE(*,934)
                      WRITE(*,914)
                      STOP
                  ENDIF
C
C                 CHECK NODE NUMBER WITHIN RANGE.
                  IF (RNODEC.GT.NNC(ICELL)) THEN
                      WRITE(*,936)
                      WRITE(*,914)
                      STOP
                  ENDIF
              END DO LINK_LOOP1
          END DO REACH_LOOP1
C
C         INPUT CHANNEL GEOMETRY TABLES.
          SECTION_LOOP:DO II=1,NTBXC
              ITBX=ITBX+1
C
C             INPUT TABLE SIZE.
              READ(IRRIV,938) NPTX(ITBX)
              IF (IECHO1.EQ.1) WRITE(IWRIV,940) II,ITBX,NPTX(ITBX)
C
C             CHECK TABLE SIZE.
              IF (NPTX(ITBX).GT.MAXPTX) THEN
                  WRITE(*,942)
                  WRITE(*,914)
                  STOP
              ENDIF
C
C             LOOP OVER TABLE POINTS.
              DO IPTX=1,NPTX(ITBX)
                  READ(IRRIV,944) QXTB(ITBX,IPTX),WXTB(ITBX,IPTX),
     +                DXTB(ITBX,IPTX)
                  IF (IECHO1.EQ.1) WRITE(IWRIV,946) IPTX,
     +                QXTB(ITBX,IPTX),WXTB(ITBX,IPTX),DXTB(ITBX,IPTX)
              END DO
          END DO SECTION_LOOP
C
C         INPUT LATERAL INFLOW TABLES.
          INFLOW_LOOP:DO II=1,NTBQC
              ITBQ=ITBQ+1
C
C             INPUT TABLE SIZE.
              READ(IRRIV,948) NPTQ(ITBQ),QFACT,IGET(ITBQ),QNAME(ITBQ)
              IF (IECHO1.EQ.1) WRITE(IWRIV,950) II,ITBQ,
     +            QNAME(ITBQ),NPTQ(ITBQ),QFACT,IGET(ITBQ)
C
C             CHECK FOR VALID INTERPOLATION METHOD.
              IF (IGET(ITBQ).LT.1 .OR. IGET(ITBQ).GT.2) THEN
                  WRITE(*,951)
                  WRITE(*,914)
                  STOP
              ENDIF
C
C             CHECK TABLE SIZE.
              IF (NPTQ(ITBQ).GT.MAXPTQ) THEN
                  WRITE(*,952)
                  WRITE(*,914)
                  STOP
              ENDIF
C
C             LOOP OVER TABLE POINTS.
              DO IPTQ=1,NPTQ(ITBQ)
                  READ(IRRIV,954) TTB(ITBQ,IPTQ),QTB(ITBQ,IPTQ)
                  QTB(ITBQ,IPTQ)=QTB(ITBQ,IPTQ)*QFACT
                  IF (IECHO1.EQ.1) WRITE(IWRIV,956) IPTQ,
     +                TTB(ITBQ,IPTQ),QTB(ITBQ,IPTQ)
              END DO
          END DO INFLOW_LOOP
C
C         INPUT HYDROGRAPH LOCATIONS.
          IF (NGPHC.GT.0) THEN
              IF (IECHO1.NE.0) WRITE(IWRIV,958) NGPHC
              DO II=1,NGPHC
                  IGPH=IGPH+1
                  READ(IRRIV,960) GREACHC,GLINK(IGPH),GNAME(IGPH)
                  GREACH(IGPH)=GREACHC+LLAST
                  IF (IECHO1.NE.0) WRITE(IWRIV,962) GREACHC,
     +                GREACH(IGPH),GLINK(IGPH),GNAME(IGPH)
              END DO
          ENDIF
          ILAST=ILAST+NNC(ICELL)
          JLAST=JLAST+NTBXC
          KLAST=KLAST+NTBQC
          LLAST=LLAST+NRC(ICELL)
      END DO CELL_LOOP
C
C     UPDATE JOINING REACHES FROM LOCAL TO GLOBAL NUMBER.
      IF (IECHO1.NE.0 .AND. NJOIN.NE.0) WRITE(IWRIV,964) NRCH,NJOIN,
     +    NTBX,NTBQ
      DO IRCH=1,NRCH
          IF (NJOIN.EQ.0) EXIT
          DO ILNK=1,NRN(IRCH)
              IF (JOIN(IRCH,ILNK).NE.0) THEN
C
C                 FIND COMPARTMENT.
                  JLAST=0
                  DO ICELL=1,NCELL
                      IF (CJOIN(IRCH,ILNK).EQ.CNAME(ICELL)) THEN
C
C                         UPDATE TO GLOBAL REACH NUMBER.
                          JOINC=JOIN(IRCH,ILNK)
                          JOIN(IRCH,ILNK)=JOINC+JLAST
                          WRITE(IWRIV, 966) IRCH,ILNK,
     +                        JOIN(IRCH,ILNK),JOINC,CJOIN(IRCH,ILNK)
                          EXIT
                      ENDIF
                      JLAST=JLAST+NRC(ICELL)
                  END DO
              ENDIF
          END DO
      END DO
C
C     IDENTIFY UPSTREAM TO DOWNSTREAM PROCESSING ORDER.
      DO IRCH=1,NRCH
          ORDER(IRCH)=IRCH
      END DO
      DO
          DONE=.TRUE.
          ORDER_LOOP:DO II=1,NRCH
              IRCH=ORDER(II)
              DO ILNK=1,NRN(IRCH)
                  JRCH=JOIN(IRCH,ILNK)
                  IF (JRCH.EQ.0) CYCLE
                  DO JJ=1,NRCH
                      IF (ORDER(JJ).EQ.JRCH) EXIT
                  END DO
                  IF (II.LT.JJ) THEN
                      ORDER(II)=JRCH
                      ORDER(JJ)=IRCH
                      DONE=.FALSE.
                      EXIT ORDER_LOOP
                  ENDIF
              END DO
          END DO ORDER_LOOP
          IF (DONE) EXIT
      END DO
      WRITE(IWRIV,968)
      DO II=1,NRCH
          IRCH=ORDER(II)
          WRITE(IWRIV,970) II,IRCH,SNAME(IRCH)
      END DO
C
C     CLOSE FILE AND EXIT WRIVER.
      WRITE(*,972)
      WRITE(IWRIV,974)
      CLOSE (IRRIV)
      RETURN
C
C     FORMATS:
900   FORMAT(1X,'*START READING IN RIVER')
902   FORMAT(1X,'**START READING IN RIVER**')
904   FORMAT(/1X,'FILE NAMES:'/
     +    1X,'INPUT FILE..............................',A80/
     +    1X,'PLOT OUTPUT FILE........................',A80/
     +    1X,'GRAPH OUTPUT FILE.......................',A80)
906   FORMAT(A20)
908   FORMAT(6I10)
910   FORMAT(/1X,'STREAM NETWORK PARAMETERS:'/
     +    1X,'NUMBER OF COMPARTMENTS..................',I10/
     +    1X,'COMPARTMENT NUMBER......................',I10/
     +    1X,'COMPARTMENT NAME........................',A20/
     +    1X,'NUMBER OF CHANNEL REACHES...............',I10/
     +    1X,'NUMBER OF GEOMETRY TABLES...............',I10/
     +    1X,'NUMBER OF INFLOW TABLES.................',I10/
     +    1X,'NUMBER OF STREAMFLOW HYDROGRAPHS........',I10/
     +    1X,'SWITCH FOR ECHO OF INPUTS...............',I10/
     +    1X,'SWITCH FOR ECHO OF SUMULATION...........',I10)
912   FORMAT(1X,'*ERROR* COMPARTMENT NAME NOT IN LIST')
914   FORMAT(1X,'*ERROR* RUNTEMINATED FOR ERROR IN WRIVER')
916   FORMAT(1X,'*ERROR* NUMBER OF REACHES EXCEEDS DIMENSION')
918   FORMAT(1X,'*ERROR* NUMBER OF SECTION TABLES EXCEEDS DIMENSION')
920   FORMAT(1X,'*ERROR* NUMBER OF INFLOW TABLES EXCEEDS DIMENSION')
922   FORMAT(2I10,F10.0,1X,A20)
924   FORMAT(/1X'REACH PARAMETERS:'/
     +    1X,'LOCAL REACH NUMBER......................',I10/
     +    1X,'GLOBAL REACH NUMBER.....................',I10/
     +    1X,'NUMBER OF LINKS.........................',I10/
     +    1X,'FLAG ENDING REACH.......................',I10/
     +    1X,'REACH NAME..............................',A20/
     +    1X,'CONDUCTIVITY SCALING FACTOR.............',1PE10.2//
     +    1X,75X,'LOCAL',4X,'GLOBAL',5X,'LOCAL',4X,'GLOBAL'/
     +    1X,15X,'LOCAL',4X,'GLOBAL',43X,'SECTION',3X,'SECTION',
     +        4X,'INFLOW',4X,'INFLOW',3X,'JOINING',1X,'JOINING'/
     +    1X,6X,'LINK',6X,'NODE',6X,'NODE',4X,'LENGTH',4X,'HEIGHT',
     +        6X,'PERM',1X,'THICKNESS',5X,'TABLE',5X,'TABLE',
     +        5X,'REACH',5X,'REACH',5X,'REACH',1X,'NAME')
925   FORMAT(1X,'*ERROR* NUMBER OF LINKS EXCEEDS DIMENSION')
926   FORMAT(I10,4F10.0,3I10,1X,A20)
928   FORMAT(1X,3I10,2F10.2,1PE10.2,0PF10.2,5I10,1X,A20)
929   FORMAT(1X,'*ERROR* COMPARTMENT NAME FOR REACH NOT IN LIST')
930   FORMAT(1X,'*ERROR* NUMBER OF NODES EXCEEDS DIMENSION')
932   FORMAT(1X,'*ERROR* TABLE EXCEEDS NUMBER OF LOCAL SECTION TABLES')
934   FORMAT(1X,'*ERROR* TABLE EXCEEDS NUMBER OF LOCAL INFLOW TABLES')
936   FORMAT(1X,'*ERROR* NODE NUMBER EXCEEDS NUMBER IN MESH')
938   FORMAT(I10)
940   FORMAT(/1X,'GEOMETRY TABLE PARAMETERS:'/
     +    1X,'LOCAL TABLE NUMBER......................',I10/
     +    1X,'GLOBAL TABLE NUMBER.....................',I10/
     +    1X,'NUMBER OF TABLE POINTS..................',I10//
     +    1X,5X,'POINT',1X,'DISCHARGE',5X,'WIDTH',5X,'DEPTH')
942   FORMAT(1X,'*ERROR* NUMBER OF SECTION POINTS EXCEEDS DIMENSION')
944   FORMAT(3F10.0)
946   FORMAT(1X,I10,3F10.1)
948   FORMAT(I10,F10.0,I10,1X,A20)
950   FORMAT(/1X,'INFLOW TABLE PARAMETERS:'/
     +    1X,'LOCAL TABLE NUMBER......................',I10/
     +    1X,'GLOBAL TABLE NUMBER.....................',I10/
     +    1X,'INFLOW NAME.............................',A20/
     +    1X,'NUMBER OF TABLE POINTS..................',I10/
     +    1X,'SCALING FACTOR FOR INFLOWS..............',1PE10.2/
     +    1X,'SWITCH FOR INTERPOLATION METHOD.........',I10//
     +    1X,5X,'POINT',6X,'TIME',1X,'DISCHARGE')
951   FORMAT(1X,'*ERROR* INTPERLATION METHOD NOT VALID')
952   FORMAT(1X,'*ERROR* NUMBER OF INFLOW POINTS EXCEEDS DIMENSION')
954   FORMAT(2F10.0)
956   FORMAT(1X,I10,2F10.1)
958   FORMAT(/1X,'HYDROGRAPH PARAMETERS:'/
     +    1X,'NUMBER OF HYDROGRAPHS...................',I10//
     +    1X,5X,'LOCAL',4X,'GLOBAL',11X,'SITE'/
     +    1X,5X,'REACH',5X,'REACH',6X,'LINK',1X,'NAME'/)
960   FORMAT(2I10,1X,A20)
962   FORMAT(1X,3I10,X,A20)
964   FORMAT(/1X,'NETWORK CONNECTIVITY PARAMETERS:'/
     +    1X,'NUMBER OF REACHES.......................',I10/
     +    1X,'NUMBER OF INTERSECTIONS.................',I10/
     +    1X,'GLOBAL NUMBER OF SECTION TABLES.........',I10/
     +    1X,'GLOBAL NUMBER OF INFLOW TABLES..........',I10//
     +    1X,'UPDATED JOINING REACHES:'/
     +    1X,4X,'GLOBAL',14X,'GLOBAL',5X,'LOCAL'/
     +    1X,5X,'REACH',6X,'LINK',6X,'JOIN',6X,'JOIN',1X,'COMPARTMENT')
966   FORMAT(1X,4I10,1X,A20)
968   FORMAT(/1X,'PROCESSING ORDER:'/
     +    1X,5X,'ORDER',5X,'REACH',21X,'REACH NAME')
970   FORMAT(1X,2I10,21X,A20)
972   FORMAT(1X,'*END   READING IN RIVER')
974   FORMAT(/1X,'**END READING IN RIVER**')
C
C     ----------------
      ENTRY WRIVER2(H)
C     ----------------
C
C     DEFINITION OF WATER-BUDGET QUANTITIES
C         Q1  INFLOW FROM UPSTREAM
C         Q2  INFLOW FROM TRIBUTARY REACH
C         Q3  INFLOW FROM LATERAL INFLOW
C         Q4  OUTFLOW THROUGH RIVER BED
C         Q5  OUTFLOW TO DOWNSTREAM
C
C     LOOP OVER REACHES AND LINKS WITHIN REACHES.
      REACH_LOOP2: DO JJ=1,NRCH
          IRCH=ORDER(JJ)
C
C         LOOP OVER LINKS WITHIN REACH.
          LINK_LOOP2: DO ILNK=1,NRN(IRCH)
C
C             SET VALUE FOR UPSTREAM INFLOW TO REACH.
              Q1(IRCH,ILNK)=0.0
              IF (ILNK.NE.1) THEN
                  Q1(IRCH,ILNK)=MAX(Q5(IRCH,ILNK-1),0.0)
              ENDIF
C
C             SET VALUE FOR INFLOWS FROM TRIBUTARY REACHES.
              Q2(IRCH,ILNK)=0.0
              IF (JOIN(IRCH,ILNK).NE.0) THEN
                  JRCH=JOIN(IRCH,ILNK)
                  JLNK=NRN(JRCH)
                  Q2(IRCH,ILNK)=MAX(Q5(JRCH,JLNK),0.0)
              ENDIF
C
C             SET VALUE FOR LATERAL INFLOW.
              Q3(IRCH,ILNK)=0.0
              IF (TABQ(IRCH,ILNK).NE.0) THEN
                  ITBQ=TABQ(IRCH,ILNK)
                  DO IPTQ=1,NPTQ(ITBQ)
                      XTTB(IPTQ)=TTB(ITBQ,IPTQ)
                      XQTB(IPTQ)=QTB(ITBQ,IPTQ)
                  END DO
                  Q3(IRCH,ILNK)=TABLE(NPTQ(ITBQ),XTTB,XQTB,STIME-DELT,
     +                STIME,IGET(ITBQ))
              ENDIF
C
C             IF LATERAL INFLOW IS DIVERSION, POSSIBLY ADJUST DIVERSION
C                 TO REMOVE NEGATIVE STREAMFLOW.
              QSUM=Q1(IRCH,ILNK)+Q2(IRCH,ILNK)+Q3(IRCH,ILNK)
              IF (QSUM.LT.0.0) THEN
                  QSUM=0.0
                  Q3(IRCH,ILNK)=-Q1(IRCH,ILNK)-Q2(IRCH,ILNK)
              END IF
C
C             SETUP CHANNEL GEOMETRY TABLE.
              ITBX=TABX(IRCH,ILNK)
              DO IPTX=1,NPTX(ITBX)
                  XQXTB(IPTX)=QXTB(ITBX,IPTX)
                  XWXTB(IPTX)=WXTB(ITBX,IPTX)
                  XDXTB(IPTX)=DXTB(ITBX,IPTX)
              END DO
C
C             INTERPOLATE FLOW WIDTH AND DEPTH.
              IF (QSUM.GT.0.0) THEN
                  W=TABLE(NPTX(ITBX),XQXTB,XWXTB,QSUM,QSUM,1)
                  D=TABLE(NPTX(ITBX),XQXTB,XDXTB,QSUM,QSUM,1)
              ELSE
                  W=0.0
                  D=0.0
              ENDIF
              WMAX=XWXTB(NPTX(ITBX))
              ZSURF=ZBED(IRCH,ILNK)+D
              ZBOT=ZBED(IRCH,ILNK)-BBED(IRCH,ILNK)
C
C             IDENTIFY MESH NODE AND LEAKANCE.
              II=RNODE(IRCH,ILNK)
              LEAK=KBED(IRCH,ILNK)/BBED(IRCH,ILNK)
C
C             IDENTIFY LEFTHAND AND RIGHTHAND COEFFICIENTS FOR CASE.
C
C             CASE 1: INFLOW IS ZERO AND GROUNDWATER LEVEL BELOW
C                 CHANNEL BED.
              IF (QSUM.LE.0.0 .AND. H(II).LE.ZSURF) THEN
                  ICASE(IRCH,ILNK)=1
                  CLH(IRCH,ILNK)=0.0
                  CRH(IRCH,ILNK)=0.0
C
C             CASE 2: INFLOW NOT ZERO, GROUNDWATER LEVEL BELOW STREAM
C                 SURFACE, AND HYDRAULIC DISCONNECT
              ELSEIF (QSUM.GT.0.0 .AND. H(II).LE.ZBOT) THEN
C
C                 CALCULATE RECHARGE CAPACITY WITHIN LINK.
                  WETL=LENGTH(IRCH,ILNK)
                  RCAP=LEAK*WETL*W*(ZSURF-ZBOT)
C
C                 CALCULATE LEFTHAND AND RIGHTHAND COEFFICIENTS.
                  ICASE(IRCH,ILNK)=2
                  CLH(IRCH,ILNK)=0.0
                  CRH(IRCH,ILNK)=MIN(RCAP,QSUM)
C
C             CASE 3: INFLOW NOT ZERO, GROUNDWATER LEVEL BELOW STREAM
C                 SURFACE, NOT HYDRAULIC DISCONNECT.
              ELSEIF (QSUM.GT.0.0 .AND. H(II).LE.ZSURF .AND.
     +            H(II).GT.ZBOT) THEN
C
C                 CALCULATE RECHARGE CAPACITY OF LINK.
                  WETL=LENGTH(IRCH,ILNK)
                  RCAP=LEAK*WETL*W*(ZSURF-H(II))
C
C                 ADJUST WETTED LENGTH OF LINK.
                  IF (RCAP.GT.QSUM) THEN
                      WETL=WETL*QSUM/RCAP
                  ENDIF
C
C                 CALCULATE LEFTHAND AND RIGHTHAND COEFFICIENTS.
                  ICASE(IRCH,ILNK)=3
                  CLH(IRCH,ILNK)=LEAK*WETL*W
                  CRH(IRCH,ILNK)=LEAK*WETL*W*ZSURF
C
C             CASE 4: GROUNDWATER LEVEL ABOVE STREAM SURFACE.
              ELSEIF (H(II).GT.ZSURF) THEN
C
C                 CALCULATE LEFTHAND AND RIGHTHAND COEFFICIENTS.
                  ICASE(IRCH,ILNK)=4
                  WETL=LENGTH(IRCH,ILNK)
                  CLH(IRCH,ILNK)=LEAK*WETL*WMAX
                  CRH(IRCH,ILNK)=LEAK*WETL*WMAX*ZSURF
              ENDIF
C
C             CALCULATE RECHARGE AND LINK OUTFLOW
              Q4(IRCH,ILNK)=CRH(IRCH,ILNK)-CLH(IRCH,ILNK)*H(II)
              Q5(IRCH,ILNK)=QSUM-Q4(IRCH,ILNK)
C
C             PACK COEFFICIENTS INTO LEFTHAND MATRIX AND RIGHTHAND VECTOR.
              CALL PACK2(II,II,CLH(IRCH,ILNK),1)
              F(II)=F(II)+CRH(IRCH,ILNK)
          END DO LINK_LOOP2
      END DO REACH_LOOP2
      RETURN
C
C     ----------------------
      ENTRY WRIVER3(H,SUMRV)
C     ----------------------
C
C     INITIALIZE OUTPUTS.
      SUMQ3=0.0
      SUMQ4=0.0
      SUMQ5=0.0
      IF (ISTEP.EQ.ISTEP1) WRITE(IWRIV,976)
C
C     UPDATE BUDGETS FOR RIVER LINKS
      REACH_LOOP3:DO JJ=1,NRCH
          IRCH=ORDER(JJ)
          IF (IPRN.EQ.1 .AND. IECHO2.EQ.1) THEN
              WRITE(IWRIV,978) IRCH,SNAME(IRCH),ISTEP,STIME
          ENDIF
          LINK_LOOP3:DO ILNK=1,NRN(IRCH)
C
C             SET VALUE FOR UPSTREAM INFLOW TO REACH.
              Q1(IRCH,ILNK)=0.0
              IF (ILNK.NE.1) THEN
                  Q1(IRCH,ILNK)=MAX(Q5(IRCH,ILNK-1),0.0)
              ENDIF
C
C             SET VALUE FOR INFLOWS FROM TRIBUTARY REACHES.
              Q2(IRCH,ILNK)=0.0
              IF (JOIN(IRCH,ILNK).NE.0) THEN
                  JRCH=JOIN(IRCH,ILNK)
                  JLNK=NRN(JRCH)
                  Q2(IRCH,ILNK)=MAX(Q5(JRCH,JLNK),0.0)
              ENDIF
C
C             SET VALUE FOR LATERAL INFLOW.
              Q3(IRCH,ILNK)=0.0
              IF (TABQ(IRCH,ILNK).NE.0) THEN
                  ITBQ=TABQ(IRCH,ILNK)
                  DO IPTQ=1,NPTQ(ITBQ)
                      XTTB(IPTQ)=TTB(ITBQ,IPTQ)
                      XQTB(IPTQ)=QTB(ITBQ,IPTQ)
                  END DO
                  Q3(IRCH,ILNK)=TABLE(NPTQ(ITBQ),XTTB,XQTB,STIME-DELT,
     +                STIME,IGET(ITBQ))
              ENDIF
C
C             IF LATERAL INFLOW IS DIVERSION, POSSIBLY ADJUST DIVERSION
C                 TO REMOVE NEGATIVE STREAMFLOW.
              QSUM=Q1(IRCH,ILNK)+Q2(IRCH,ILNK)+Q3(IRCH,ILNK)
              IF (QSUM.LT.0.0) THEN
                  QSUM=0.0
                  Q3(IRCH,ILNK)=-Q1(IRCH,ILNK)-Q2(IRCH,ILNK)
              END IF
C
C             SETUP CHANNEL GEOMETRY TABLE.
              ITBX=TABX(IRCH,ILNK)
              DO IPTX=1,NPTX(ITBX)
                  XQXTB(IPTX)=QXTB(ITBX,IPTX)
                  XWXTB(IPTX)=WXTB(ITBX,IPTX)
                  XDXTB(IPTX)=DXTB(ITBX,IPTX)
              END DO
C
C             INTERPOLATE FLOW WIDTH AND DEPTH.
              IF (QSUM.GT.0.0) THEN
                  W=TABLE(NPTX(ITBX),XQXTB,XWXTB,QSUM,QSUM,1)
                  D=TABLE(NPTX(ITBX),XQXTB,XDXTB,QSUM,QSUM,1)
              ELSE
                  W=0.0
                  D=0.0
              ENDIF
              WMAX=XWXTB(NPTX(ITBX))
              ZSURF=ZBED(IRCH,ILNK)+D
              ZBOT=ZBED(IRCH,ILNK)-BBED(IRCH,ILNK)
C
C             CALCULATE RECHARGE AND LINK OUTFLOW.
              II=RNODE(IRCH,ILNK)
              Q4(IRCH,ILNK)=CRH(IRCH,ILNK)-CLH(IRCH,ILNK)*H(II)
              Q5(IRCH,ILNK)=QSUM-Q4(IRCH,ILNK)
C
C             ACCUMULATE FLOWS AND OUTPUT FOR LINK.
              SUMQ3=SUMQ3+Q3(IRCH,ILNK)
              SUMQ4=SUMQ4+Q4(IRCH,ILNK)
              RESID=Q1(IRCH,ILNK)+Q2(IRCH,ILNK)+Q3(IRCH,ILNK)-
     +            Q4(IRCH,ILNK)-Q5(IRCH,ILNK)
              IF (IPRN.EQ.1 .AND. IECHO2.NE.0) WRITE(IWRIV,980)
     +            RNODE(IRCH,ILNK),ZSURF,H(II),W,D,
     +            Q1(IRCH,ILNK),Q2(IRCH,ILNK),Q3(IRCH,ILNK),
     +            Q4(IRCH,ILNK),Q5(IRCH,ILNK),ICASE(IRCH,ILNK),RESID
          END DO LINK_LOOP3
          IF (LAST(IRCH).EQ.1) SUMQ5=SUMQ5+Q5(IRCH,NRN(IRCH))
      END DO REACH_LOOP3
      SUMRV=SUMQ4
C
C     UPDATE HYDROGRAPH VECTORS.
      DO IGPH=1,NGPH
          IRCH=GREACH(IGPH)
          ILNK=GLINK(IGPH)
          QGRAPH(IGPH,ISTEP)=Q5(IRCH,ILNK)
          GTIME(ISTEP)=STIME
      END DO

C
C     COMPILE AND DISPLAY SUMMARY BUDGET.
      RESID=SUMQ3-SUMQ4-SUMQ5
      IF (IPRN.EQ.1) WRITE(IWRIV,982) ISTEP,STIME,SUMQ3,SUMQ4,SUMQ5,
     +    RESID
C
C     CLOSE FILES AND EXIT WRIVER.
      IF (ISTEP.EQ.ISTEP2) THEN
          DO IGPH=1,NGPH
              WRITE(IWRIVG,984) GNAME(IGPH)
              DO JSTEP=ISTEP1,ISTEP2
                  WRITE(IWRIVG,986) JSTEP,GTIME(JSTEP),
     +                QGRAPH(IGPH,JSTEP)
              END DO
          END DO
          WRITE(IWRIV,988)
          CLOSE (IWRIV)
          CLOSE (IWRIVG)
      ENDIF
      RETURN
C
C     FORMATS:
976   FORMAT(1X,'**START SIMULATION IN WRIVER**')
978   FORMAT(/1X,'INFLOW AND OUTFLOW FOR REACH:'/
     +    1X,'REACH NUMBER............................',I10/
     +    1X,'REACH NAME..............................',A20/
     +    1X,'TIME STEP...............................',I10/
     +    1X,'ELAPSED TIME............................',1PE10.2//
     +    1X,6X,'NODE',3X,'SURFACE',6X,'HEAD',5X,'WIDTH',5X,'DEPTH',
     +        2X,'UPSTREAM',1X,'TRIBUTARY',3X,'LATERAL',2X,'RECHARGE',
     +        3X,'OUTFLOW',6X,'CASE',2X,'RESIDUAL')
980   FORMAT(1X,I10,4F10.2,5(1PE10.2),I10,1PE10.2)
982   FORMAT(/1X,'SUMMARY RIVER BUDGET:'/
     +    1X,'TIME STEP...............................',I10/
     +    1X,'ELAPSED TIME............................',1PE10.2/
     +    1X,'LATERAL INFLOWS.........................',1PE10.2/
     +    1X,'NET GROUNDWATER RECHARGE................',1PE10.2/
     +    1X,'OUTFLOW FROM LAST REACHES...............',1PE10.2/
     +    1X,'RESIDUAL................................',1PE10.2)
984   FORMAT(1X,'SITE: ',A20/
     +    1X,6X,'STEP',6X,'TIME',1X,'DISCHARGE')
986   FORMAT(1X,I10,2(1PE10.2))
988   FORMAT(/1X,'**END SIMULATION IN WRIVER**')
      END SUBROUTINE
