C  **************************************************************
C  * LPMOD                                                      *
C  *                                                            *
C  * THIS PROGRAM GIVES SOLUTIONS FOR AN UNSTEADY, NONUNIFORM   *      
C  * FLOW, MOBILE-BED MATHEMATICAL MODEL.  ITS INTENDED         *     
C  * APPLICATION IS MODELING OF ADJUSTMENTS IN RIVER LONGITUD-  *     
C  * INAL PROFILE AND OTHER STREAM CHANNEL PARAMETERS.          *     
C  * SI UNITS ARE USED THROUGHOUT (ANY ODD UNITS ARE NOTED).    *
C  *      THIS CODE WAS WRITTEN SPECIFICALLY TO RUN ON VAX      *
C  * MAINFRAMES, BUT THE ONLY MODIFICATIONS NEEDED FOR OTHER    *
C  * SYSTEMS SHOULD BE CHANGES IN INPUT/OUTPUT COMMANDS.        *
C  *      THIS VERSION IS MODIFIED TO SIMULATE CRUSTAL UPLIFT.  *
C  **************************************************************     
C 
      REAL*8 U(200),H(200),B(200),W(200),QI(200),D(200),E(200)
      REAL*8 QSI(200),HC(200),BCHA(200),BHOLD(200),TAUC(200)
      REAL*8 QS(200),QSO(200),DELX(200),LOC(200),SLOPE(200) 
      REAL*8 DELB(200),DEEP(200),UAUN(200),CH(200),CHN(200)         
      REAL*8 WD(200),WUD(200),WU2D(200),WHCMB(200),WQS(200) 
      REAL*8 UN(200),HN(200),BN(200),WN(200),QIN(200),QSIN(200)       
      REAL*8 QSN(200),U2DN(200),WAUN(200),DU(200),WUAU(200),UP(200)
      REAL*8 DQSDUN(200),DQSDDN(200),DEEPN(200),HCMBN(200),QSIC(200)
      REAL*8 WDN(200),WUDN(200),WU2DN(200),WHCMBN(200),WQSN(200)      
      REAL*8 WUAUN(200),WU2N(200),DH(200),DB(200),DW(200)   
      REAL*8 UDN(200),WUN(200),QSG(200),GC2(200),GC2N(200) 
      REAL*8 TH,THM,TH2,THM2,NEW,D2T,ADH,ADB,ADW,ADU,DELT,OLDB
      REAL*8 ALP,EC,S,RHO,RHOS,SEDMOV,G,QTOP,QSTOP,CONCRI   
      REAL*8 UINC,DINC,QSINC,TIME,Q,DELTY,BET
      REAL*8 JJ(4,200),KK(4,200),LL(4,200),MM(4,200),RR(4,200)
      REAL*8 A1,A2,A3,A4,B1,B2,B3,B4,C1,C2,C3,C4  
      REAL*8 D1,D2,D3,D4,E1,E2,E3,E4,F1,F2,F3,F4  
      REAL*8 G1,G2,G3,G4,H1,H2,H3,H4,R1,R2,R3,R4  
      REAL*8 J1,J2,J3,J4,K1,K2,K3,K4,L1,L2,L3,L4  
      REAL*8 M1,M2,M3,M4,EE
      INTEGER I,J,RUN,IMAX,T,NTS,X,L,XM,XX,LM,ITMAX,ITER,CON
      INTEGER HTOP  
       OPEN (UNIT=1,FILE='OUTPUT.DAT',STATUS='UNKNOWN')
       OPEN (UNIT=21,FILE='INPUT.DAT',STATUS='OLD')
       OPEN (UNIT=22,FILE='G1.DAT',STATUS='UNKNOWN')
       OPEN (UNIT=23,FILE='CONTIN.DAT',STATUS='UNKNOWN')
       OPEN (UNIT=24,FILE='G2.DAT',STATUS='UNKNOWN')
       OPEN (UNIT=25,FILE='CONTROL.DAT',STATUS='OLD')
C         
C READ IN CONTROL VARIABLES FROM "CONTROL" FILE:   
C
C   RUN    - I.D. NUMBER (FOR RECORDS PURPOSES)
C   L      - NUMBER OF DATA POINTS
C   NTS    - NUMBER OF TIME STEPS
C   DELTY  - LENGTH OF TIME STEP (IN YEARS)
C   SEDMOV - VALUE CONTROLLING PRESENCE (1.0) OR ABSENCE (0.0) OF
C            SEDIMENT TRANSPORT
C   HTOP   - OPTION TO ALLOW NO BED ELEVATION CHANGE AT THE FIVE
C            MOST UPSTREAM POINTS (HTOP=1) OR LET THEM CHANGE (0)
C            (OCCASIONALLY USEFUL FOR DAMPING NUMERICAL ERRORS)
C   JPRINT - NUMBER OF TIME STEPS TO PROCESS PER PRINTING OF 
C            OUTPUT DATA (ALLOWS REDUCTION OF DATA COLLECTED FOR 
C            RUNS WITH MANY TIME STEPS)...SET AT 1 IF WISH DATA
C            FOR EVERY TIME STEP
      READ (25,*) RUN
      READ (25,*) L
      READ (25,*) NTS
      READ (25,*) DELTY
      READ (25,*) SEDMOV
      READ (25,*) HTOP
      READ (25,*) JPRINT
C
C OTHER CONTROL VARIABLES ("MANUAL" INPUT AREA):
C
C   FRICT  - PRESENCE (1.0) OR ABSENCE (0.0) OF FLOW RESISTANCE
C   ITMAX  - MAXIMUM NUMBER OF ITERATIONS TO TRY AT ANY TIME
C   TH     - THETA, TIME WEIGHTING FACTOR FOR SOLUTION
C   CONCRI - CONVERGENCE CRITERION FOR ADEQUATE APPROXIMATION OF
C            DEPENDENT VARIABLE VALUES
      FRICT = 1.0D0 
      ITMAX = 20    
      TH = .550D0       
      CONCRI = .0001D0        
C         
C GLOBAL VARIABLES AND PHYSICAL CONSTANTS (INPUT AREA):      
C 
C   QTOP  - WATER DISCHARGE FROM UPSTREAM AT X = 0
C   QSTOP - SEDIMENT DISCHARGE FROM UPSTREAM AT X = 0
C   BET   - BETA, LATERAL INFLOW DOWNSTREAM MOMENTUM FACTOR
C   G     - ACCEL. OF GRAVITY
C   ALP   - ALPHA, (1 - POROSITY) OF CHANNEL SEDIMENTS
C   RHO   - WATER DENSITY
C   RHOS  - SEDIMENT DENSITY
      QTOP = 5.0D0
      QSTOP = 0.001D0
      BET = 0.0D0  
      G = 9.80D0 
      ALP = 0.7D0    
      RHO = 1000.0D0         
      RHOS = 2650.0D0         
C         
C CALCULATE GENERAL CONVENIENCE VARIABLES         
C         
      DELT = DELTY*31557600.0D0         
      LM = L-1      
      D2T = DELT*2.0D0        
      THM = 1.0D0 - TH        
      TH2 = TH/2.0D0
      THM2 = THM/2.0D0
      T = 0         
C         
C INITIAL CONDITIONS, POSITION-SPECIFIC VALUES, FROM "INPUT" FILE      
C         
C   TIME  - TIME FROM BEGINNING OF SIMULATION, IN YEARS
C   LOC   - LOCATIONS OF DATA POINTS (KM DOWNSTREAM), IN ORDER
C           FROM UPSTREAM TO DOWNSTREAM
C   U     - FLOW VELOCITY
C   H     - WATER SURFACE ELEVATION
C   B     - AVG. BED SURFACE ELEVATION
C   W     - CHANNEL WIDTH
C   QI    - LATERAL INFLOW OF WATER (CU. METERS PER SEC PER KM)
C   D     - SEDIMENT DIAMETER (METERS)
C   TAUC  - CRITICAL SHEAR STRESS OR OTHER CHARACTERISTIC OF
C           SEDIMENT USED IN "SEDMOV" RELATIONSHIP
C   QSI   - LATERAL INFLOW OF SEDIMENT (SAME UNITS AS QI)
C   E     - CROSS-SECTION CHANGE FACTOR (INIT. SET AT ZERO)
C   BCHA  - PREVIOUS TOTAL BED ELEVATION CHANGE (INIT. SET ZERO)
C   UP    - CRUSTAL UPLIFT (METERS PER YEAR)
      READ (21,*) TIME
      READ (21,*) (LOC(X),X=1,L)
      READ (21,*) (U(X),X=1,L)         
      READ (21,*) (H(X),X=1,L)      
      READ (21,*) (B(X),X=1,L)
      READ (21,*) (W(X),X=1,L)   
      READ (21,*) (QI(X),X=1,L)         
      READ (21,*) (D(X),X=1,L)       
      READ (21,*) (TAUC(X),X=1,L)       
      READ (21,*) (QSIC(X),X=1,L)                    
      READ (21,*) (E(X),X=1,L)
      READ (21,*) (BCHA(X),X=1,L)       
      READ (21,*) (UP(X),X=1,L)
C         
C ADDITIONAL POSITION-SPECIFIC VALUES, ADJUSTED INFLOW UNITS
C         
       DO 80 X=1,L  
        HC(X) = H(X) - BCHA(X)        
        DELB(X) = 0.0         
        UN(X) = U(X)
        HN(X) = H(X)
        BN(X) = B(X)
        BHOLD(X) = B(X)- BCHA(X)        
        WN(X) = W(X)
        IF (X.NE.L) DELX(X) = (LOC(X+1)-LOC(X))*1000.0D0
        QI(X) = QI(X)/1000.0D0
        QIN(X) = QI(X)        
        QSI(X) = QSIC(X)/1000.0D0 + UP(X)*W(X)*ALP/31557600.0D0
        QSIN(X) = QSI(X)      
 80   CONTINUE      
C         
C OUTPUT CONTROL VALUES TO HARD COPY FILE  
C         
      WRITE (1,4) RUN
 4    FORMAT (' RUN ',I7)
      WRITE (1,6) L
 6    FORMAT (' NUMBER OF DATA POINTS: ',I4)
      WRITE (1,5) NTS,DELTY
 5    FORMAT (' NO. TIME STEPS: ',I4,'    TIME STEP: ',F16.8,
     1  ' YEARS',/)
      WRITE (1,219) T,TIME    
         WRITE (1,40)         
         WRITE (1,50)         
       DO 70 X=1,L  
         WRITE (1,60) LOC(X),DELX(X),QI(X),QSI(X),D(X),TAUC(X)        
  70   CONTINUE     
  40  FORMAT ('0','  LOCATION      DELX       WATER INFLOW    ',   
     1 '   SED INFLOW      SED SIZE    CRIT VALUE')  
  50  FORMAT (' ','     KM          M           CU.M/S        ',
     1 '     CU.M/S           M                  ',/)
  60  FORMAT (2F12.5,2F18.14,2F12.7)    
      J = JPRINT
C         
C BEGIN TIME STEPS  
C         
      DO 400 T=1,NTS
       ITER=0       
C CALCULATE VALUES (QS, ETC) FOR OLD TIME STEP    
C (DONE ONCE EACH TIME STEP)  
       DO 100 X=1,L 
        DEEP(X) = H(X) - B(X) 
        WD(X) = W(X)*DEEP(X)  
        WUD(X) = WD(X)*U(X)   
        WU2D(X) = WUD(X)*U(X) 
        WHCMB(X) = W(X)*(HC(X)-B(X))    
        WUAU(X) = W(X)*U(X)*DABS(U(X))
        CALL XSECT(E(X),W(X),LOC(X),HC(X),TIME)
        CALL RESIST(CH(X),DEEP(X),D(X))  
        GC2(X) = G/(CH(X)**2)
C   (SEDIMENT TRANSPORT)
        CALL STRANS(QS(X),U(X),DEEP(X),CH(X),RHO,RHOS,D(X),TAUC(X))
        IF (X.EQ.1) QS(1) = QSTOP/W(1)
        QS(X) = SEDMOV*QS(X)
        WQS(X) = QS(X)*W(X)
        QSO(X) = WQS(X)*3600.0D0
C         
C OUTPUT DATA FOR OLD TIME STEP TO GRAPHICS FILE, HARD COPY FILE  
C
      IF (J.NE.JPRINT) GO TO 100
       WRITE (24,270) TIME,LOC(X),U(X),QS(X),B(X),BCHA(X)
 100  CONTINUE     
      TIME = TIME + DELT/31557600.0D0  
      IF (J.NE.JPRINT) GO TO 700 
       WRITE (1,106)
       WRITE (1,107)
      DO 110 X=1,L  
       WRITE (1,108) LOC(X),U(X),H(X),B(X),W(X),DEEP(X),DELB(X),      
     1   BCHA(X),QSO(X),QSI(X)  
 110  CONTINUE      
 106  FORMAT (' ','**  LOCATION    VELOCITY   WATER ELEV   ',
     1  'BED ELEV     WIDTH       DEPTH    BED CHA RATE  ',
     2  'EROS/DEP  SED DISCHARGE   SED INFLOW')   
 107  FORMAT (' ','       KM         M/S          M        ',
     1  '   M           M           M         MM/DAY     ',
     2  '   M        CU.M/HR        CU.M/S/M',/)        
 108  FORMAT (' ','* ',9F12.5,F18.14)           
C         
C CALCULATE VALUES (QS, ETC) FOR "NEW TIME STEP"  
C (DONE ONCE EACH ITERATION)  
C         
 700   ITER = ITER + 1
       DO 800 X=1,L
        DEEPN(X) = (HN(X)-BN(X))        
        HCMBN(X) = (HC(X)-BN(X))        
        UDN(X) = UN(X)*DEEPN(X)         
        U2DN(X) = UN(X)*UDN(X)
        WDN(X) = WN(X)*DEEPN(X)         
        WUDN(X) = UN(X)*WDN(X)
        WU2DN(X) = UN(X)*WUDN(X)        
        WHCMBN(X) = WN(X)*HCMBN(X)      
        WAUN(X) = WN(X)*DABS(UN(X))     
        WUAUN(X) = WAUN(X)*UN(X)        
        WUN(X) =WN(X)*UN(X)
        WU2N(X) = WUN(X)*UN(X)
        UAUN(X) = UN(X)*DABS(UN(X))
        CALL RESIST(CHN(X),DEEPN(X),D(X))     
        GC2N(X) = G/(CHN(X)**2)
C   (SEDIMENT TRANSPORT & PARTIAL DERIVATIVES)
        CALL STRANS(QSN(X),UN(X),DEEPN(X),CHN(X),RHO,RHOS,
     1       D(X),TAUC(X))
        IF (X.EQ.1) QSN(1) = QSTOP/WN(1)
        QSN(X) = SEDMOV*QSN(X)
        WQSN(X)=QSN(X)*WN(X)
        UINC = UN(X) + .001D0
        CALL STRANS(QSINC,UINC,DEEPN(X),CHN(X),RHO,RHOS,
     1       D(X),TAUC(X))
        QSINC = SEDMOV*QSINC
        DQSDUN(X) = (QSINC-QSN(X))/.001D0
        DINC = DEEPN(X) + .001D0
        CALL STRANS(QSINC,UN(X),DINC,CHN(X),RHO,RHOS,
     1       D(X),TAUC(X))
        QSINC = SEDMOV*QSINC
        DQSDDN(X) = (QSINC-QSN(X))/.001D0
        IF (X.EQ.1) DQSDUN(X) = 0.0D0
        IF (X.EQ.1) DQSDDN(X) = 0.0D0
 800   CONTINUE     
       CON = 1      
C         
C  FORMULATE MODEL EQUATIONS AND BEGIN SOLUTION FOR THE 
C  ITERATION, PROGRESSING FROM X=1 TO X=L    
C         
       DO 200 X=1,L 
        IF (X.NE.1) GO TO 130 
C  
C BOUNDARY CONDITIONS AT X=1:         
C   HYDRAULIC BOUNDARY CONDITION (DISCHARGE)      
C   .....F2 MUST BE NON-ZERO  
        A2 = 0.0D0  
        B2 = 0.0D0  
        C2 = 0.0D0  
        D2 = 0.0D0  
        E2 = WDN(1) 
        F2 = WUN(1) 
        G2 = -WUN(1)
        H2 = UDN(1) 
        J2 = 0.0D0  
        K2 = 0.0D0  
        L2 = 0.0D0  
        M2 = 0.0D0  
        R2 = WUDN(1) - QTOP   
        R2 = -R2    
C   SEDIMENT-RELATED BOUNDARY CONDITIONS (NO CHANGE THROUGH TIME)     
C   .....G3 MUST BE NON-ZERO  
        A3 = 0.0D0  
        B3 = 0.0D0  
        C3 = 0.0D0  
        D3 = 0.0D0  
        E3 = 0.0D0  
        F3 = 0.0D0  
        G3 = 1.0D0  
        H3 = 0.0D0  
        J3 = 0.0D0  
        K3 = 0.0D0  
        L3 = -1.0D0 
        M3 = 0.0D0  
        R3 = BN(1) - B(1) - BN(2) + B(2)
        R3 = -R3 
C   .....H4 MUST BE NON-ZERO  
        A4 = 0.0D0  
        B4 = 0.0D0  
        C4 = 0.0D0  
        D4 = 0.0D0  
        E4 = 0.0D0  
        F4 = 0.0D0  
        G4 = E(1)      
        H4 = 1.0D0  
        J4 = 0.0D0  
        K4 = 0.0D0  
        L4 = 0.0D0  
        M4 = 0.0D0  
        R4 = E(1)*BN(1) - E(1)*B(1) + WN(1) - W(1)      
        R4 = -R4    
 130    CONTINUE    
        IF (X.EQ.L) GO TO 140 
C         
C GENERAL EQUATIONS:
C   CONSERVATION OF MOMENTUM FOR CLEAR WATER      
       SLOPE(X) = (TH*(HN(X+1)-HN(X))+THM*(H(X+1)-H(X)))/DELX(X)
       A1 = 0.0D0  
       B1 = 0.0D0  
       C1 = 0.0D0  
       D1 = 0.0D0  
       E1 = WDN(X)/D2T-2.0D0*TH*WUDN(X)/DELX(X)+GC2N(X)*TH*WAUN(X)
     1     -BET*TH2*QIN(X)
       F1 = WUN(X)/D2T-TH*WU2N(X)/DELX(X)-G*(TH2*(WDN(X+1)+WDN(X))    
     1     +THM2*(WD(X+1)+WD(X)))*TH/DELX(X)+WN(X)*G*TH2*SLOPE(X)     
       G1 = -WUN(X)/D2T+TH*WU2N(X)/DELX(X)-G*TH2*WN(X)*SLOPE(X)       
       H1 = UDN(X)/D2T-TH*U2DN(X)/DELX(X)+G*TH2*DEEPN(X)*SLOPE(X)     
     1     +UAUN(X)*GC2N(X)*TH2   
       J1=WDN(X+1)/D2T+2.0D0*TH*WUDN(X+1)/DELX(X)
     1     +GC2N(X+1)*TH*WAUN(X+1)-BET*TH2*QIN(X+1)
       K1 = WUN(X+1)/D2T+TH*WU2N(X+1)/DELX(X)+G*(TH2*(WDN(X+1)
     1     +WDN(X))+THM2*(WD(X+1)+WD(X)))*TH/DELX(X)
     2     +WN(X+1)*G*TH2*SLOPE(X)   
       L1 = -WUN(X+1)/D2T-TH*WU2N(X+1)/DELX(X)
     1     -G*TH2*WN(X+1)*SLOPE(X) 
       M1 = UDN(X+1)/D2T+TH*U2DN(X+1)/DELX(X)
     1     +G*TH2*DEEPN(X+1)*SLOPE(X)+UAUN(X+1)*GC2N(X+1)*TH2 
       R1 = (WUDN(X)+WUDN(X+1)-WUD(X)-WUD(X+1))/D2T         
     1     +TH*(WU2DN(X+1)-WU2DN(X))/DELX(X)      
     2     +THM*(WU2D(X+1)-WU2D(X))/DELX(X)
     3     +G*(TH2*(WDN(X+1)+WDN(X))
     4     +THM2*(WD(X+1)+WD(X)))*SLOPE(X)   
     5     +TH2*(GC2N(X+1)*WUAUN(X+1)+GC2N(X)*WUAUN(X))
     6     +THM2*(GC2(X+1)*WUAU(X+1)+GC2(X)*WUAU(X))
     7     -BET*(TH2*(QIN(X)*UN(X)+QIN(X+1)*UN(X+1))
     8     +THM2*(QI(X)*U(X)+QI(X+1)*U(X+1)))
        R1 = -R1        
        IF (X.EQ.1) GO TO 170 
 140    CONTINUE    
C   CONTINUITY OF WATER       
       A2 = -TH*WDN(X-1)/DELX(X-1)      
       B2 = WN(X-1)/D2T - TH*WUN(X-1)/DELX(X-1)   
       C2 = -WN(X-1)/D2T +TH*WUN(X-1)/DELX(X-1)   
       D2 = DEEPN(X-1)/D2T - TH*UDN(X-1)/DELX(X-1)
       E2 = TH*WDN(X)/DELX(X-1)         
       F2 = WN(X)/D2T + TH*WUN(X)/DELX(X-1)       
       G2 = -WN(X)/D2T - TH*WUN(X)/DELX(X-1)      
       H2 = DEEPN(X)/D2T + TH*UDN(X)/DELX(X-1)    
       J2 = 0.0D0   
       K2 = 0.0D0   
       L2 = 0.0D0   
       M2 = 0.0D0   
       R2 = (WDN(X)+WDN(X-1)-WD(X)-WD(X-1))/D2T   
     1     + TH*(WUDN(X)-WUDN(X-1))/DELX(X-1)     
     2     + THM*(WUD(X)-WUD(X-1))/DELX(X-1)      
     3     - TH2*(QIN(X)+QIN(X-1))-THM2*(QI(X)+QI(X-1))     
        R2 = -R2    
        IF (X.EQ.L) GO TO 150 
C   CONTINUITY OF SEDIMENT    
       A3 = TH*WN(X-1)*DQSDUN(X-1)/(DELX(X-1)+DELX(X))      
       B3 = TH*WN(X-1)*DQSDDN(X-1)/(DELX(X-1)+DELX(X))      
       C3 = -TH*WN(X-1)*DQSDDN(X-1)/(DELX(X-1)+DELX(X))     
       D3 = TH*QSN(X-1)/(DELX(X-1)+DELX(X))       
       E3 = 0.0D0   
       F3 = 0.0D0   
       G3 = -ALP*WN(X)/DELT   
       H3 = ALP*HCMBN(X)/DELT 
       J3 = -TH*WN(X+1)*DQSDUN(X+1)/(DELX(X-1)+DELX(X))     
       K3 = -TH*WN(X+1)*DQSDDN(X+1)/(DELX(X-1)+DELX(X))     
       L3 = TH*WN(X+1)*DQSDDN(X+1)/(DELX(X-1)+DELX(X))      
       M3 = -TH*QSN(X+1)/(DELX(X-1)+DELX(X))      
       R3 = ALP*(WHCMBN(X)-WHCMB(X))/DELT         
     1     -TH*(WQSN(X+1)-WQSN(X-1))/(DELX(X-1)+DELX(X))    
     2     -THM*(WQS(X+1)-WQS(X-1))/(DELX(X-1)+DELX(X))     
     3     + (TH*QSIN(X) + THM*QSI(X))*SEDMOV     
       R3 = -R3     
       IF (HTOP.NE.1) GO TO 145
C   OPTION TO HOLD TOP 5 POINTS WITH NO BED ELEV. CHANGE
       IF (X.GT.5) GO TO 145
       A3 = 0.0D0
       B3 = 0.0D0
       C3 = 0.0D0
       D3 = 0.0D0
       E3 = 0.0D0
       F3 = 0.0D0
       G3 = 1.0D0
       H3 = 0.0D0
       J3 = 0.0D0
       K3 = 0.0D0
       L3 = 0.0D0
       M3 = 0.0D0
       R3 = BN(X) - B(X)
       R3 = -R3
 145   CONTINUE         
C   CHANGE OF CROSS-SECTION GEOMETRY  
       A4 = 0.0D0   
       B4 = 0.0D0   
       C4 = 0.0D0   
       D4 = 0.0D0   
       E4 = 0.0D0   
       F4 = 0.0D0   
       G4 = E(X)       
       H4 = 1.0D0   
       J4 = 0.0D0   
       K4 = 0.0D0   
       L4 = 0.0D0   
       M4 = 0.0D0   
       R4 = E(X)*BN(X) - E(X)*B(X) + WN(X) - W(X)       
       R4 = -R4     
        IF (X.NE.L) GO TO 160 
 150    CONTINUE    
C         
C BOUNDARY CONDITIONS AT X=L:
        A1 = A2     
        B1 = B2     
        C1 = C2     
        D1 = D2     
        E1 = E2     
        F1 = F2     
        G1 = G2     
        H1 = H2     
        J1 = J2     
        K1 = K2     
        L1 = L2     
        M1 = M2     
        R1 = R2              
C   HYDRAULIC BOUNDARY CONDITION (WATER SURFACE ELEVATION)  
C   .....F2 MUST BE NON-ZERO  
        A2 = 0.0D0  
        B2 = 0.0D0  
        C2 = 0.0D0  
        D2 = 0.0D0  
        E2 = 0.0D0  
        F2 = 1.0D0  
        G2 = 0.0D0  
        H2 = 0.0D0  
        J2 = 0.0D0  
        K2 = 0.0D0  
        L2 = 0.0D0  
        M2 = 0.0D0  
        R2 = HN(L) -H(L)      
        R2 = -R2    
C   SEDIMENT-RELATED BOUNDARY COND. (BACKWARD-DIFFERENCE APPROX.)
C   .....G3 MUST BE NON-ZERO  
       A3 = TH*WN(LM)*DQSDUN(LM)/DELX(X-1)        
       B3 = TH*WN(LM)*DQSDDN(LM)/DELX(X-1)        
       C3 = -TH*WN(LM)*DQSDDN(LM)/DELX(X-1)       
       D3 = TH*QSN(LM)/DELX(X-1)        
       E3 = -TH*WN(L)*DQSDUN(L)/DELX(X-1)         
       F3 = -TH*WN(L)*DQSDDN(L)/DELX(X-1)         
       G3 = -ALP*WN(L)/DELT + TH*WN(L)*DQSDDN(L)/DELX(X-1)  
       H3 = ALP*HCMBN(L)/DELT - TH*QSN(L)/DELX(X-1)         
       J3 = 0.0D0   
       K3 = 0.0D0   
       L3 = 0.0D0   
       M3 = 0.0D0   
       R3 = ALP*(WHCMBN(L)-WHCMB(L))/DELT         
     1     - TH*(WQSN(L)-WQSN(LM))/DELX(X-1)      
     2     - THM*(WQS(L)-WQS(LM))/DELX(X-1)       
     3     + (TH*QSIN(L) + THM*QSI(L))*SEDMOV     
       R3 = -R3     
C   .....H4 MUST BE NON-ZERO  
       A4 = 0.0D0   
       B4 = 0.0D0   
       C4 = 0.0D0   
       D4 = 0.0D0   
       E4 = 0.0D0   
       F4 = 0.0D0   
       G4 = E(L)       
       H4 = 1.0D0   
       J4 = 0.0D0   
       K4 = 0.0D0   
       L4 = 0.0D0   
       M4 = 0.0D0   
       R4 = E(L)*BN(L) - E(L)*B(L) + WN(L) - W(L)       
       R4 = -R4     
 160    CONTINUE    
C         
C REMOVE A-B-C-D VALUES FROM MATRIX     
C         
        E1 = E1 - A1*JJ(1,X-1) - B1*JJ(2,X-1) - C1*JJ(3,X-1)
     1       - D1*JJ(4,X-1)   
        F1 = F1 - A1*KK(1,X-1) - B1*KK(2,X-1) - C1*KK(3,X-1)
     1       - D1*KK(4,X-1)   
        G1 = G1 - A1*LL(1,X-1) - B1*LL(2,X-1) - C1*LL(3,X-1)
     1       - D1*LL(4,X-1)   
        H1 = H1 - A1*MM(1,X-1) - B1*MM(2,X-1) - C1*MM(3,X-1)
     1       - D1*MM(4,X-1)   
        R1 = R1 - A1*RR(1,X-1) - B1*RR(2,X-1) - C1*RR(3,X-1)
     1       - D1*RR(4,X-1)   
        E2 = E2 - A2*JJ(1,X-1) - B2*JJ(2,X-1) - C2*JJ(3,X-1)
     1       - D2*JJ(4,X-1)   
        F2 = F2 - A2*KK(1,X-1) - B2*KK(2,X-1) - C2*KK(3,X-1)
     1       - D2*KK(4,X-1)   
        G2 = G2 - A2*LL(1,X-1) - B2*LL(2,X-1) - C2*LL(3,X-1)
     1       - D2*LL(4,X-1)   
        H2 = H2 - A2*MM(1,X-1) - B2*MM(2,X-1) - C2*MM(3,X-1)
     1       - D2*MM(4,X-1)   
        R2 = R2 - A2*RR(1,X-1) - B2*RR(2,X-1) - C2*RR(3,X-1)
     1       - D2*RR(4,X-1)   
        E3 = E3 - A3*JJ(1,X-1) - B3*JJ(2,X-1) - C3*JJ(3,X-1)
     1       - D3*JJ(4,X-1)   
        F3 = F3 - A3*KK(1,X-1) - B3*KK(2,X-1) - C3*KK(3,X-1)
     1       - D3*KK(4,X-1)   
        G3 = G3 - A3*LL(1,X-1) - B3*LL(2,X-1) - C3*LL(3,X-1)
     1       - D3*LL(4,X-1)   
        H3 = H3 - A3*MM(1,X-1) - B3*MM(2,X-1) - C3*MM(3,X-1)
     1       - D3*MM(4,X-1)   
        R3 = R3 - A3*RR(1,X-1) - B3*RR(2,X-1) - C3*RR(3,X-1)
     1       - D3*RR(4,X-1)   
        E4 = E4 - A4*JJ(1,X-1) - B4*JJ(2,X-1) - C4*JJ(3,X-1)
     1       - D4*JJ(4,X-1)   
        F4 = F4 - A4*KK(1,X-1) - B4*KK(2,X-1) - C4*KK(3,X-1)
     1       - D4*KK(4,X-1)   
        G4 = G4 - A4*LL(1,X-1) - B4*LL(2,X-1) - C4*LL(3,X-1)
     1       - D4*LL(4,X-1)   
        H4 = H4 - A4*MM(1,X-1) - B4*MM(2,X-1) - C4*MM(3,X-1)
     1       - D4*MM(4,X-1)   
        R4 = R4 - A4*RR(1,X-1) - B4*RR(2,X-1) - C4*RR(3,X-1)
     1       - D4*RR(4,X-1)   
 170    CONTINUE    
C         
C PERFORM GAUSSIAN ELIMINATION ON E-F-G-H PART OF MATRIX    
C         
        EE = E2/E1  
        F2 = F2 - F1*EE       
        G2 = G2 - G1*EE       
        H2 = H2 - H1*EE       
        J2 = J2 - J1*EE       
        K2 = K2 - K1*EE       
        L2 = L2 - L1*EE       
        M2 = M2 - M1*EE       
        R2 = R2 - R1*EE       
        EE = E3/E1  
        F3 = F3 - F1*EE       
        G3 = G3 - G1*EE       
        H3 = H3 - H1*EE       
        J3 = J3 - J1*EE       
        K3 = K3 - K1*EE       
        L3 = L3 - L1*EE       
        M3 = M3 - M1*EE       
        R3 = R3 - R1*EE      
        EE = E4/E1  
        F4 = F4 - F1*EE       
        G4 = G4 - G1*EE       
        H4 = H4 - H1*EE       
        J4 = J4 - J1*EE       
        K4 = K4 - K1*EE       
        L4 = L4 - L1*EE       
        M4 = M4 - M1*EE       
        R4 = R4 - R1*EE       
        EE = F3/F2  
        G3 = G3 - G2*EE       
        H3 = H3 - H2*EE       
        J3 = J3 - J2*EE       
        K3 = K3 - K2*EE       
        L3 = L3 - L2*EE       
        M3 = M3 - M2*EE       
        R3 = R3 - R2*EE       
        EE = F4/F2  
        G4 = G4 - G2*EE       
        H4 = H4 - H2*EE       
        J4 = J4 - J2*EE       
        K4 = K4 - K2*EE       
        L4 = L4 - L2*EE       
        M4 = M4 - M2*EE       
        R4 = R4 - R2*EE       
        EE = G4/G3  
        H4 = H4 - H3*EE       
        J4 = J4 - J3*EE       
        K4 = K4 - K3*EE       
        L4 = L4 - L3*EE       
        M4 = M4 - M3*EE       
        R4 = R4 - R3*EE       
        F1 = F1/E1  
        G1 = G1/E1  
        H1 = H1/E1  
        J1 = J1/E1  
        K1 = K1/E1  
        L1 = L1/E1  
        M1 = M1/E1  
        R1 = R1/E1  
        G2 = G2/F2  
        H2 = H2/F2  
        J2 = J2/F2  
        K2 = K2/F2  
        L2 = L2/F2  
        M2 = M2/F2  
        R2 = R2/F2  
        H3 = H3/G3  
        J3 = J3/G3  
        K3 = K3/G3  
        L3 = L3/G3  
        M3 = M3/G3  
        R3 = R3/G3  
        J4 = J4/H4  
        K4 = K4/H4  
        L4 = L4/H4  
        M4 = M4/H4  
        R4 = R4/H4  
        JJ(4,X) = J4
        KK(4,X) = K4
        LL(4,X) = L4
        MM(4,X) = M4
        RR(4,X) = R4
        JJ(3,X) = J3 - J4*H3  
        KK(3,X) = K3 - K4*H3  
        LL(3,X) = L3 - L4*H3  
        MM(3,X) = M3 - M4*H3  
        RR(3,X) = R3 - R4*H3  
        JJ(2,X) = J2 - J4*H2 - JJ(3,X)*G2         
        KK(2,X) = K2 - K4*H2 - KK(3,X)*G2         
        LL(2,X) = L2 - L4*H2 - LL(3,X)*G2         
        MM(2,X) = M2 - M4*H2 - MM(3,X)*G2         
        RR(2,X) = R2 - R4*H2 - RR(3,X)*G2         
        JJ(1,X) = J1 - J4*H1 - JJ(3,X)*G1 - JJ(2,X)*F1      
        KK(1,X) = K1 - K4*H1 - KK(3,X)*G1 - KK(2,X)*F1      
        LL(1,X) = L1 - L4*H1 - LL(3,X)*G1 - LL(2,X)*F1      
        MM(1,X) = M1 - M4*H1 - MM(3,X)*G1 - MM(2,X)*F1
        RR(1,X) = R1 - R4*H1 - RR(3,X)*G1 - RR(2,X)*F1
 200   CONTINUE     
       DU(L) = RR(1,L)        
       DH(L) = RR(2,L)        
       DB(L) = RR(3,L)        
       DW(L) = RR(4,L)        
C   (CHECK FOR CONVERGENCE)     
        ADU = DABS(DU(L))     
        ADH = DABS(DH(L))     
        ADB = DABS(DB(L))     
        ADW = DABS(DW(L))     
        IF (ADU.GT.CONCRI) CON = 0      
        IF (ADH.GT.CONCRI) CON = 0      
        IF (ADB.GT.CONCRI) CON = 0      
        IF (ADW.GT.CONCRI) CON = 0      
C         
C COMPLETE THE SOLUTION FOR THE ITERATION,
C PROGRESSING FROM X=L TO X=1  
C         
       DO 300 XX=1,LM         
        X = L-XX      
        DU(X) = RR(1,X) - DU(X+1)*JJ(1,X) - DH(X+1)*KK(1,X) 
     1         - DB(X+1)*LL(1,X) - DW(X+1)*MM(1,X)
        DH(X) = RR(2,X) - DU(X+1)*JJ(2,X) - DH(X+1)*KK(2,X) 
     1         - DB(X+1)*LL(2,X) - DW(X+1)*MM(2,X)
        DB(X) = RR(3,X) - DU(X+1)*JJ(3,X) - DH(X+1)*KK(3,X) 
     1         - DB(X+1)*LL(3,X) - DW(X+1)*MM(3,X)
        DW(X) = RR(4,X) - DU(X+1)*JJ(4,X) - DH(X+1)*KK(4,X) 
     1         - DB(X+1)*LL(4,X) - DW(X+1)*MM(4,X)
C   (CHECK FOR CONVERGENCE)     
        ADU = DABS(DU(X))     
        ADH = DABS(DH(X))     
        ADB = DABS(DB(X))     
        ADW = DABS(DW(X))     
        IF (ADU.GT.CONCRI) CON = 0      
        IF (ADH.GT.CONCRI) CON = 0      
        IF (ADB.GT.CONCRI) CON = 0      
        IF (ADW.GT.CONCRI) CON = 0           
 300   CONTINUE
       IF (CON.EQ.1) GO TO 500
       IF (ITER.EQ.ITMAX) GO TO 450     
C         
C  NO CONVERGENCE:  
C         
       DO 425 X=1,L 
        UN(X) = UN(X) + DU(X) 
        HN(X) = HN(X) + DH(X) 
        BN(X) = BN(X) + DB(X) 
        WN(X) = WN(X) + DW(X) 
 425   CONTINUE           
       GO TO 700    
C         
C  TOO MANY ITERATIONS:       
C         
 450   WRITE (1,430)
       WRITE (6,430)
 430   FORMAT (' ','REACHED MAX. ALLOWED ITERATIONS--NO CONVERGENCE',/)         
       WRITE (1,106)
       WRITE (1,107)
      DO 431 X=1,L  
       WRITE (1,108) LOC(X),UN(X),HN(X),BN(X),WN(X),DEEPN(X),DELB(X), 
     1   BCHA(X),QSO(X),QSI(X)  
 431  CONTINUE      
       STOP         
C         
C  CONVERGENCE:     
C         
C   CALCULATE NEW VALUES AT NEW TIME STEP, EXTRAPOLATE INITIAL        
C   GUESSES FOR NEXT TIME STEP, CALCULATE ADDITIONAL OUTPUT VALUES     
 500   DO 550 X=1,L 
        NEW = UN(X)+DU(X)     
        UN(X) = NEW+NEW-U(X)  
        U(X) = NEW  
        NEW = HN(X)+DH(X)     
        HN(X) = NEW+NEW-H(X)  
        H(X) = NEW  
        NEW = BN(X)+DB(X)     
        DELB(X) = NEW-B(X)    
        BN(X) = NEW+DELB(X)   
        DELB(X) = DELB(X)*86400000.0D0/DELT       
        B(X) = NEW  
        BCHA(X) = B(X) - BHOLD(X)       
        NEW = WN(X)+DW(X)     
        WN(X) = NEW+NEW-W(X)  
        W(X) = NEW  
 550   CONTINUE     
       WRITE (6,217) T
 217   FORMAT (' TIME STEP ',I3,' COMPLETED')
        J = J+1     
        IF (J.GT.JPRINT) J=1  
        IF (J.NE.JPRINT) GO TO 261      
C   (OUTPUT HEADINGS FOR NEW TIME STEP TO HARD COPY FILE)
       WRITE (1,219) T,TIME   
 219   FORMAT ('0',8X,'BELOW VALUES AT T = ',I3,5X,'TIME = ',
     1   F16.4,'  YEARS')  
       WRITE (1,220) ITER     
 220   FORMAT (' ','AFTER ',I2,' ITERATIONS')     
 261     CONTINUE        
 400   CONTINUE     
C         
C END OF TIME STEPS 
C         
       T = NTS      
C         
C OUTPUT FINAL VALUES:
C  ...HARD COPY FILE
       WRITE (1,106)
       WRITE (1,107)       
       DO 280 X=1,L 
        DEEP(X) = H(X)-B(X)
        QSO(X) = QSN(X)*WN(X)*3600.0D0
       WRITE (1,108) LOC(X),U(X),H(X),B(X),W(X),DEEP(X),DELB(X),      
     1   BCHA(X),QSO(X),QSI(X)  
C  ...GRAPHICS FILE
       WRITE (24,270) TIME,LOC(X),U(X),QSN(X),B(X),BCHA(X)        
 270   FORMAT (F11.4,1X,F7.3,1X,F8.4,1X,F12.10,1X,F8.3,1X,F9.5)      
C  ...SECOND GRAPHICS FILE, RECORD OF FINAL VALUES    
        Q = DEEP(X)*U(X)*W(X)       
        WRITE (22,271) RUN,LOC(X),H(X),B(X),Q,D(X),TAUC(X)  
 271    FORMAT (2X,I3,4(2X,F8.3),2(2X,F8.6))
C  ...FINAL VALUES TO CONTINUATION FILE "CONTIN", USED AS INPUT      
C     TO ALLOW SIMULATION TO RESUME WHERE THIS RUN LETS OFF 
      QI(X) = QI(X)*1000.0D0
 280  CONTINUE
      WRITE (23,409) TIME     
 409  FORMAT (F18.6)
      WRITE (23,410) (LOC(X),X=1,L)
      WRITE (23,410) (U(X),X=1,L)       
      WRITE (23,410) (H(X),X=1,L)       
      WRITE (23,410) (B(X),X=1,L)       
      WRITE (23,410) (W(X),X=1,L)       
      WRITE (23,410) (QI(X),X=1,L)      
      WRITE (23,410) (D(X),X=1,L)       
      WRITE (23,410) (TAUC(X),X=1,L)    
      WRITE (23,411) (QSIC(X),X=1,L)     
      WRITE (23,410) (E(X),X=1,L)      
      WRITE (23,410) (BCHA(X),X=1,L)
      WRITE (23,410) (UP(X),X=1,L)    
 410  FORMAT (5F13.6)
 411  FORMAT (5F13.10)
      CLOSE (UNIT=1)
      CLOSE (UNIT=21)
      CLOSE (UNIT=22)
      CLOSE (UNIT=23)
      CLOSE (UNIT=24)
      CLOSE (UNIT=25)
C  ...NOTE TO USER AT TERMINAL
      WRITE (6,291)
 291  FORMAT ('0.......CALCULATIONS COMPLETED.........',/)
      WRITE (6,292)
 292  FORMAT ('     RESULTS LISTED IN FILE, OUTPUT.DAT',/)
      STOP
      END 
C
C  RELATIONSHIP FOR CHANNEL CROSS-SECTION ADJUSTMENT
C  (SIMPLE FACTOR PRODUCING GIVEN PERCENT CHANGE IN
C  WIDTH FOR EACH METER CHANGE IN BED ELEVATION)
C
      SUBROUTINE XSECT(E,W,LOC,HC,TIME)
      REAL*8 E,W,LOC,TIME,EC,HC
      IF (TIME.NE.0.0D0) GO TO 11
      EC = 0.0D0
      E = W*EC
  11  CONTINUE
      END
C
C  FLOW RESISTANCE RELATIONSHIP
C  (MANNING)
C
      SUBROUTINE RESIST(CH,DEEP,D)
      REAL*8 CH,DEEP,D,N
      N = 0.03D0
      CH = (DEEP**0.166)/N
      END
C
C  SEDIMENT TRANSPORT RELATIONSHIP
C  (ENGELUND BEDLOAD EQN.)
C
      SUBROUTINE STRANS(QS,U,DEEP,CH,RHO,RHOS,D,TAUC)
      REAL*8 QS,U,DEEP,CH,RHO,RHOS,D,TAUC
      REAL*8 KN1,KN2,UDIR,TAU,THR1,THR2,G
      G = 9.8D0
      KN1 = 5.0D0*D*DSQRT((RHOS-RHO)*G*D/RHO)    
      KN2 = RHO/((RHOS-RHO)*D*CH**2)      
      IF (U.NE.0) UDIR = U/(DABS(U))   
      IF (U.EQ.0) UDIR = 1         
      TAU = U**2*KN2  
      THR1 = TAU - TAUC  
      IF (THR1.LE.0) THR1=0.0D0       
      IF (THR1.GT.0)  THR1 = (1.0D0+(0.267D0/THR1)**4)**(-0.25)
      THR2 = DSQRT(TAU) -  0.7D0*DSQRT(TAUC) 
      IF (THR2.LE.0) THR2=0.0D0       
      QS = UDIR*KN1*THR1*THR2
      END
