!  **************************************************************
 
! Code converted using TO_F90 by Alan Miller
! Date: 2005-11-20  Time: 21:16:14
 
!  * LPMOD                                                      *
!  *                                                            *
!  * THIS PROGRAM GIVES SOLUTIONS FOR AN UNSTEADY, NONUNIFORM   *
!  * FLOW, MOBILE-BED MATHEMATICAL MODEL.  ITS INTENDED         *
!  * APPLICATION IS MODELING OF ADJUSTMENTS IN RIVER LONGITUD-  *
!  * INAL PROFILE AND OTHER STREAM CHANNEL PARAMETERS.          *
!  * SI UNITS ARE USED THROUGHOUT (ANY ODD UNITS ARE NOTED).    *
!  *      THIS CODE WAS WRITTEN SPECIFICALLY TO RUN ON VAX      *
!  * MAINFRAMES, BUT THE ONLY MODIFICATIONS NEEDED FOR OTHER    *
!  * SYSTEMS SHOULD BE CHANGES IN INPUT/OUTPUT COMMANDS.        *
!  *      THIS VERSION IS MODIFIED TO SIMULATE CRUSTAL UPLIFT.  *
!  **************************************************************
 
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')
 
! READ IN CONTROL VARIABLES FROM "CONTROL" FILE:
 
!   RUN    - I.D. NUMBER (FOR RECORDS PURPOSES)
!   L      - NUMBER OF DATA POINTS
!   NTS    - NUMBER OF TIME STEPS
!   DELTY  - LENGTH OF TIME STEP (IN YEARS)
!   SEDMOV - VALUE CONTROLLING PRESENCE (1.0) OR ABSENCE (0.0) OF
!            SEDIMENT TRANSPORT
!   HTOP   - OPTION TO ALLOW NO BED ELEVATION CHANGE AT THE FIVE
!            MOST UPSTREAM POINTS (HTOP=1) OR LET THEM CHANGE (0)
!            (OCCASIONALLY USEFUL FOR DAMPING NUMERICAL ERRORS)
!   JPRINT - NUMBER OF TIME STEPS TO PROCESS PER PRINTING OF
!            OUTPUT DATA (ALLOWS REDUCTION OF DATA COLLECTED FOR
!            RUNS WITH MANY TIME STEPS)...SET AT 1 IF WISH DATA
!            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
 
! OTHER CONTROL VARIABLES ("MANUAL" INPUT AREA):
 
!   FRICT  - PRESENCE (1.0) OR ABSENCE (0.0) OF FLOW RESISTANCE
!   ITMAX  - MAXIMUM NUMBER OF ITERATIONS TO TRY AT ANY TIME
!   TH     - THETA, TIME WEIGHTING FACTOR FOR SOLUTION
!   CONCRI - CONVERGENCE CRITERION FOR ADEQUATE APPROXIMATION OF
!            DEPENDENT VARIABLE VALUES
frict = 1.0D0
itmax = 20
th = .550D0
concri = .0001D0
 
! GLOBAL VARIABLES AND PHYSICAL CONSTANTS (INPUT AREA):
 
!   QTOP  - WATER DISCHARGE FROM UPSTREAM AT X = 0
!   QSTOP - SEDIMENT DISCHARGE FROM UPSTREAM AT X = 0
!   BET   - BETA, LATERAL INFLOW DOWNSTREAM MOMENTUM FACTOR
!   G     - ACCEL. OF GRAVITY
!   ALP   - ALPHA, (1 - POROSITY) OF CHANNEL SEDIMENTS
!   RHO   - WATER DENSITY
!   RHOS  - SEDIMENT DENSITY
qtop = 5.0D0
qstop = 0.001D0
bet = 0.0D0
g = 9.80D0
alp = 0.7D0
rho = 1000.0D0
rhos = 2650.0D0
 
! CALCULATE GENERAL CONVENIENCE VARIABLES
 
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
 
! INITIAL CONDITIONS, POSITION-SPECIFIC VALUES, FROM "INPUT" FILE
 
!   TIME  - TIME FROM BEGINNING OF SIMULATION, IN YEARS
!   LOC   - LOCATIONS OF DATA POINTS (KM DOWNSTREAM), IN ORDER
!           FROM UPSTREAM TO DOWNSTREAM
!   U     - FLOW VELOCITY
!   H     - WATER SURFACE ELEVATION
!   B     - AVG. BED SURFACE ELEVATION
!   W     - CHANNEL WIDTH
!   QI    - LATERAL INFLOW OF WATER (CU. METERS PER SEC PER KM)
!   D     - SEDIMENT DIAMETER (METERS)
!   TAUC  - CRITICAL SHEAR STRESS OR OTHER CHARACTERISTIC OF
!           SEDIMENT USED IN "SEDMOV" RELATIONSHIP
!   QSI   - LATERAL INFLOW OF SEDIMENT (SAME UNITS AS QI)
!   E     - CROSS-SECTION CHANGE FACTOR (INIT. SET AT ZERO)
!   BCHA  - PREVIOUS TOTAL BED ELEVATION CHANGE (INIT. SET ZERO)
!   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)
 
! ADDITIONAL POSITION-SPECIFIC VALUES, ADJUSTED INFLOW UNITS
 
DO  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 /= 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)
END DO
 
! OUTPUT CONTROL VALUES TO HARD COPY FILE
 
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, ' YEARS',/)
WRITE (1,219) t,time
WRITE (1,40)
WRITE (1,50)
DO  x=1,l
  WRITE (1,60) loc(x),delx(x),qi(x),qsi(x),d(x),tauc(x)
END DO
40  FORMAT ('0','  LOCATION      DELX       WATER INFLOW    ',  &
    '   SED INFLOW      SED SIZE    CRIT VALUE')
50  FORMAT (' ','     KM          M           CU.M/S        ',  &
    '     CU.M/S           M                  ',/)
60  FORMAT (2F12.5,2F18.14,2F12.7)
j = jprint
 
! BEGIN TIME STEPS
 
DO  t=1,nts
  iter=0
! CALCULATE VALUES (QS, ETC) FOR OLD TIME STEP
! (DONE ONCE EACH TIME STEP)
  DO  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)
!   (SEDIMENT TRANSPORT)
    CALL strans(qs(x),u(x),deep(x),ch(x),rho,rhos,d(x),tauc(x))
    IF (x == 1) qs(1) = qstop/w(1)
    qs(x) = sedmov*qs(x)
    wqs(x) = qs(x)*w(x)
    qso(x) = wqs(x)*3600.0D0
 
! OUTPUT DATA FOR OLD TIME STEP TO GRAPHICS FILE, HARD COPY FILE
 
    IF (j /= jprint) CYCLE
    WRITE (24,270) time,loc(x),u(x),qs(x),b(x),bcha(x)
  END DO
  time = time + delt/31557600.0D0
  IF (j /= jprint) GO TO 700
  WRITE (1,106)
  WRITE (1,107)
  DO  x=1,l
    WRITE (1,108) loc(x),u(x),h(x),b(x),w(x),deep(x),delb(x),  &
        bcha(x),qso(x),qsi(x)
  END DO
  106  FORMAT (' ','**  LOCATION    VELOCITY   WATER ELEV   ',  &
      'BED ELEV     WIDTH       DEPTH    BED CHA RATE  ',  &
      'EROS/DEP  SED DISCHARGE   SED INFLOW')
  107  FORMAT (' ','       KM         M/S          M        ',  &
      '   M           M           M         MM/DAY     ',  &
      '   M        CU.M/HR        CU.M/S/M',/)
  108  FORMAT (' ','* ',9F12.5,f18.14)
 
! CALCULATE VALUES (QS, ETC) FOR "NEW TIME STEP"
! (DONE ONCE EACH ITERATION)
 
  700   iter = iter + 1
  DO  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)
!   (SEDIMENT TRANSPORT & PARTIAL DERIVATIVES)
    CALL strans(qsn(x),un(x),deepn(x),chn(x),rho,rhos, d(x),tauc(x))
    IF (x == 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, 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, d(x),tauc(x))
    qsinc = sedmov*qsinc
    dqsddn(x) = (qsinc-qsn(x))/.001D0
    IF (x == 1) dqsdun(x) = 0.0D0
    IF (x == 1) dqsddn(x) = 0.0D0
  END DO
  con = 1
 
!  FORMULATE MODEL EQUATIONS AND BEGIN SOLUTION FOR THE
!  ITERATION, PROGRESSING FROM X=1 TO X=L
 
  DO  x=1,l
    IF (x /= 1) GO TO 130
 
! BOUNDARY CONDITIONS AT X=1:
!   HYDRAULIC BOUNDARY CONDITION (DISCHARGE)
!   .....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
!   SEDIMENT-RELATED BOUNDARY CONDITIONS (NO CHANGE THROUGH TIME)
!   .....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
!   .....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 == l) GO TO 140
 
! GENERAL EQUATIONS:
!   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)  &
        -bet*th2*qin(x)
    f1 = wun(x)/d2t-th*wu2n(x)/delx(x)-g*(th2*(wdn(x+1)+wdn(x))  &
        +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)  &
        +uaun(x)*gc2n(x)*th2
    j1=wdn(x+1)/d2t+2.0D0*th*wudn(x+1)/delx(x)  &
        +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)  &
        +wdn(x))+thm2*(wd(x+1)+wd(x)))*th/delx(x) +wn(x+1)*g*th2*slope(x)
    l1 = -wun(x+1)/d2t-th*wu2n(x+1)/delx(x) -g*th2*wn(x+1)*slope(x)
    m1 = udn(x+1)/d2t+th*u2dn(x+1)/delx(x)  &
        +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  &
        +th*(wu2dn(x+1)-wu2dn(x))/delx(x) +thm*(wu2d(x+1)-wu2d(x))/delx(x)  &
        +g*(th2*(wdn(x+1)+wdn(x)) +thm2*(wd(x+1)+wd(x)))*slope(x)  &
        +th2*(gc2n(x+1)*wuaun(x+1)+gc2n(x)*wuaun(x))  &
        +thm2*(gc2(x+1)*wuau(x+1)+gc2(x)*wuau(x))  &
        -bet*(th2*(qin(x)*un(x)+qin(x+1)*un(x+1))  &
        +thm2*(qi(x)*u(x)+qi(x+1)*u(x+1)))
    r1 = -r1
    IF (x == 1) GO TO 170
    140    CONTINUE
!   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  &
        + th*(wudn(x)-wudn(x-1))/delx(x-1) + thm*(wud(x)-wud(x-1))/delx(x-1)  &
        - th2*(qin(x)+qin(x-1))-thm2*(qi(x)+qi(x-1))
    r2 = -r2
    IF (x == l) GO TO 150
!   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  &
        -th*(wqsn(x+1)-wqsn(x-1))/(delx(x-1)+delx(x))  &
        -thm*(wqs(x+1)-wqs(x-1))/(delx(x-1)+delx(x))  &
        + (th*qsin(x) + thm*qsi(x))*sedmov
    r3 = -r3
    IF (htop /= 1) GO TO 145
!   OPTION TO HOLD TOP 5 POINTS WITH NO BED ELEV. CHANGE
    IF (x > 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
!   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 /= l) GO TO 160
    150    CONTINUE
 
! 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
!   HYDRAULIC BOUNDARY CONDITION (WATER SURFACE ELEVATION)
!   .....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
!   SEDIMENT-RELATED BOUNDARY COND. (BACKWARD-DIFFERENCE APPROX.)
!   .....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 - th*(wqsn(l)-wqsn(lm))/delx(x-1)  &
        - thm*(wqs(l)-wqs(lm))/delx(x-1) + (th*qsin(l) + thm*qsi(l))*sedmov
    r3 = -r3
!   .....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
 
! REMOVE A-B-C-D VALUES FROM MATRIX
 
    e1 = e1 - a1*jj(1,x-1) - b1*jj(2,x-1) - c1*jj(3,x-1) - d1*jj(4,x-1)
    f1 = f1 - a1*kk(1,x-1) - b1*kk(2,x-1) - c1*kk(3,x-1) - d1*kk(4,x-1)
    g1 = g1 - a1*ll(1,x-1) - b1*ll(2,x-1) - c1*ll(3,x-1) - d1*ll(4,x-1)
    h1 = h1 - a1*mm(1,x-1) - b1*mm(2,x-1) - c1*mm(3,x-1) - d1*mm(4,x-1)
    r1 = r1 - a1*rr(1,x-1) - b1*rr(2,x-1) - c1*rr(3,x-1) - d1*rr(4,x-1)
    e2 = e2 - a2*jj(1,x-1) - b2*jj(2,x-1) - c2*jj(3,x-1) - d2*jj(4,x-1)
    f2 = f2 - a2*kk(1,x-1) - b2*kk(2,x-1) - c2*kk(3,x-1) - d2*kk(4,x-1)
    g2 = g2 - a2*ll(1,x-1) - b2*ll(2,x-1) - c2*ll(3,x-1) - d2*ll(4,x-1)
    h2 = h2 - a2*mm(1,x-1) - b2*mm(2,x-1) - c2*mm(3,x-1) - d2*mm(4,x-1)
    r2 = r2 - a2*rr(1,x-1) - b2*rr(2,x-1) - c2*rr(3,x-1) - d2*rr(4,x-1)
    e3 = e3 - a3*jj(1,x-1) - b3*jj(2,x-1) - c3*jj(3,x-1) - d3*jj(4,x-1)
    f3 = f3 - a3*kk(1,x-1) - b3*kk(2,x-1) - c3*kk(3,x-1) - d3*kk(4,x-1)
    g3 = g3 - a3*ll(1,x-1) - b3*ll(2,x-1) - c3*ll(3,x-1) - d3*ll(4,x-1)
    h3 = h3 - a3*mm(1,x-1) - b3*mm(2,x-1) - c3*mm(3,x-1) - d3*mm(4,x-1)
    r3 = r3 - a3*rr(1,x-1) - b3*rr(2,x-1) - c3*rr(3,x-1) - d3*rr(4,x-1)
    e4 = e4 - a4*jj(1,x-1) - b4*jj(2,x-1) - c4*jj(3,x-1) - d4*jj(4,x-1)
    f4 = f4 - a4*kk(1,x-1) - b4*kk(2,x-1) - c4*kk(3,x-1) - d4*kk(4,x-1)
    g4 = g4 - a4*ll(1,x-1) - b4*ll(2,x-1) - c4*ll(3,x-1) - d4*ll(4,x-1)
    h4 = h4 - a4*mm(1,x-1) - b4*mm(2,x-1) - c4*mm(3,x-1) - d4*mm(4,x-1)
    r4 = r4 - a4*rr(1,x-1) - b4*rr(2,x-1) - c4*rr(3,x-1) - d4*rr(4,x-1)
    170    CONTINUE
 
! PERFORM GAUSSIAN ELIMINATION ON E-F-G-H PART OF MATRIX
 
    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
  END DO
  du(l) = rr(1,l)
  dh(l) = rr(2,l)
  db(l) = rr(3,l)
  dw(l) = rr(4,l)
!   (CHECK FOR CONVERGENCE)
  adu = DABS(du(l))
  adh = DABS(dh(l))
  adb = DABS(db(l))
  adw = DABS(dw(l))
  IF (adu > concri) con = 0
  IF (adh > concri) con = 0
  IF (adb > concri) con = 0
  IF (adw > concri) con = 0
 
! COMPLETE THE SOLUTION FOR THE ITERATION,
! PROGRESSING FROM X=L TO X=1
 
  DO  xx=1,lm
    x = l-xx
    du(x) = rr(1,x) - du(x+1)*jj(1,x) - dh(x+1)*kk(1,x)  &
        - 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)  &
        - 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)  &
        - 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)  &
        - db(x+1)*ll(4,x) - dw(x+1)*mm(4,x)
!   (CHECK FOR CONVERGENCE)
    adu = DABS(du(x))
    adh = DABS(dh(x))
    adb = DABS(db(x))
    adw = DABS(dw(x))
    IF (adu > concri) con = 0
    IF (adh > concri) con = 0
    IF (adb > concri) con = 0
    IF (adw > concri) con = 0
  END DO
  IF (con == 1) GO TO 500
  IF (iter == itmax) GO TO 450
 
!  NO CONVERGENCE:
 
  DO  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)
  END DO
  GO TO 700
 
!  TOO MANY ITERATIONS:
 
  450   WRITE (1,430)
  WRITE (6,430)
  430   FORMAT (' ','REACHED MAX. ALLOWED ITERATIONS--NO CONVERGENCE',/)
  WRITE (1,106)
  WRITE (1,107)
  DO  x=1,l
    WRITE (1,108) loc(x),un(x),hn(x),bn(x),wn(x),deepn(x),delb(x),  &
        bcha(x),qso(x),qsi(x)
  END DO
  STOP
 
!  CONVERGENCE:
 
!   CALCULATE NEW VALUES AT NEW TIME STEP, EXTRAPOLATE INITIAL
!   GUESSES FOR NEXT TIME STEP, CALCULATE ADDITIONAL OUTPUT VALUES
  500   DO  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
  END DO
  WRITE (6,217) t
  217   FORMAT (' TIME STEP ',i3,' COMPLETED')
  j = j+1
  IF (j > jprint) j=1
  IF (j /= jprint) GO TO 261
!   (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 = ',  &
      f16.4,'  YEARS')
  WRITE (1,220) iter
  220   FORMAT (' ','AFTER ',i2,' ITERATIONS')
  261     CONTINUE
END DO
 
! END OF TIME STEPS
 
t = nts
 
! OUTPUT FINAL VALUES:
!  ...HARD COPY FILE
WRITE (1,106)
WRITE (1,107)
DO  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),  &
      bcha(x),qso(x),qsi(x)
!  ...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)
!  ...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))
!  ...FINAL VALUES TO CONTINUATION FILE "CONTIN", USED AS INPUT
!     TO ALLOW SIMULATION TO RESUME WHERE THIS RUN LETS OFF
  qi(x) = qi(x)*1000.0D0
END DO
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)
!  ...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
 
!  RELATIONSHIP FOR CHANNEL CROSS-SECTION ADJUSTMENT
!  (SIMPLE FACTOR PRODUCING GIVEN PERCENT CHANGE IN
!  WIDTH FOR EACH METER CHANGE IN BED ELEVATION)
 
SUBROUTINE xsect(e,w,loc,hc,time)
 
 
REAL*8, INTENT(OUT)                      :: e
REAL*8, INTENT(IN)                       :: w
REAL*8, INTENT(IN OUT)                   :: loc
REAL*8, INTENT(IN OUT)                   :: hc
REAL*8, INTENT(IN OUT)                   :: time
REAL*8  ec
IF (time /= 0.0D0) GO TO 11
ec = 0.0D0
e = w*ec
11  CONTINUE
END SUBROUTINE xsect
 
!  FLOW RESISTANCE RELATIONSHIP
!  (MANNING)
 
SUBROUTINE resist(ch,deep,d)
 
 
REAL*8, INTENT(OUT)                      :: ch
REAL*8, INTENT(IN OUT)                   :: deep
REAL*8, INTENT(IN OUT)                   :: d
REAL*8  n
n = 0.03D0
ch = (deep**0.166)/n
END SUBROUTINE resist
 
!  SEDIMENT TRANSPORT RELATIONSHIP
!  (ENGELUND BEDLOAD EQN.)
 
SUBROUTINE strans(qs,u,deep,ch,rho,rhos,d,tauc)
 
REAL*8, INTENT(OUT)                      :: qs
REAL*8, INTENT(IN)                       :: u
REAL*8, INTENT(IN OUT)                   :: deep
REAL*8, INTENT(IN)                       :: ch
REAL*8, INTENT(IN)                       :: rho
REAL*8, INTENT(IN OUT)                   :: rhos
REAL*8, INTENT(IN)                       :: d
REAL*8, INTENT(IN)                       :: 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 /= 0) udir = u/(DABS(u))
IF (u == 0) udir = 1
tau = u**2*kn2
thr1 = tau - tauc
IF (thr1 <= 0) thr1=0.0D0
IF (thr1 > 0)  thr1 = (1.0D0+(0.267D0/thr1)**4)**(-0.25)
thr2 = DSQRT(tau) -  0.7D0*DSQRT(tauc)
IF (thr2 <= 0) thr2=0.0D0
qs = udir*kn1*thr1*thr2
END SUBROUTINE strans
