      REAL*8 U(1000),HE(1000),H(1000),B(1000),W(1000),QS(1000),CT(1000) 
      REAL*8 DEEP(1000),S(1000),Q(1000),LOC(1000),QSI(1000),ED(1000)
      REAL*8 FR(1000),DELX,AR,RHO,INCR,QI(1000),D(1000),SETV(1000)   
      REAL*8 BASE,QTOP,QSTOP,C,INC,SGTOP,SINC,DCTDS,CON,QSG,R,CTINC   
      REAL*8 N,T,RHOS,VISK,LOCB,LOCT,LM,XR,CTG,MU,RD,FGR
      INTEGER ITR(1000),ITMAX,L,X,XX,LX,RUN        
        OPEN (UNIT=1,FILE='OUTPUT.DAT',STATUS='UNKNOWN')
	OPEN (UNIT=22,FILE='GR1.DAT',STATUS='UNKNOWN')
	OPEN (UNIT=23,FILE='SEED.DAT',STATUS='UNKNOWN')
C ************************************************************
        RUN = 53
        N = .030D0
        T = 15.0D0
        RHO = 1000.0D0
        RHOS = 2650.0D0
        G = 9.81D0
        MU = .01139D0
        VISK = .00000115D0
        LOCT = 10.0D0
        LOCB = 100.0D0
        L = 91
C ************************************************************
        BASE = 0.0D0
        ITMAX=50
        INCR=.0000001D0
	CON=.00000000005D0
	WRITE (1,4) RUN
   4    FORMAT (' RUN NUMBER ',I4,/)
      WRITE (1,5) N,T,VISK,MU,RHOS
  5   FORMAT ('  N = ',F6.4,'   T = ',F5.1,' DEGR. C   KIN. VISC. = ',
     1  F11.9,'   VISC. = ',F11.9,' POISE    SED. DENSITY = ',F7.1)
      WRITE (1,6) L
  6   FORMAT ('  NUMBER DATA PTS. = ',I4,/)
       LM = L - 1.0D0
       DELX = (LOCB-LOCT)/LM
       DO 15 X=1,L
         XR = X
         LOC(X) = LOCT + DELX*XR - DELX
C ************************************************************
         Q(X) = 100.0D0*(LOC(X)/LOCB)**1.2
         CT(X) = 1200.0D0
         D(X) = 0.0004D0*EXP(-0.003*LOC(X))
         W(X) = 50.0D0*(Q(X)/100.0D0)**.7
C ************************************************************
         RD = D(X)*50.0D0
         FGR = (G*RHO*(RHOS-RHO)*RD**2)/10000.0D0
         SETV(X) = (-3.0D0*MU + (9.0D0*(MU**2)
     1     + FGR*(.015476D0 + .19841*RD))**(0.5))
     2     /((RHO/1000.0D0)*(.011607D0 + .14881D0*RD))
C        WRITE (1,7) RD, FGR, SETV(X)
C  7     FORMAT (' RD, FGR, SETV: ',3F16.8)
         SETV(X) = .01D0*1.473D0*SETV(X)**.842D0
C        WRITE (1,8) LOC(X),SETV(X)
C  8     FORMAT (' LOC, SETV: ',2F16.8,/)
         QS(X) = CT(X)*Q(X)*RHO/(1000000.0D0*RHOS)
  15   CONTINUE
      SGTOP = .001
      DO 300 I=1,10
        CALL STRANS (CTG,SGTOP,Q(1),W(1),N,D(1),SETV(1),VISK)
        IF (CTG.GT.CT(1)) SGTOP = SGTOP/2.0D0
        IF (CTG.LT.CT(1)) SGTOP = SGTOP*2.0D0
  300 CONTINUE
      DO 100 X=1,L            
C                   
C  NEWTON ITERATION SOLVING FOR ENERGY SLOPE      
C                   
        IF (X.EQ.1) S(X) = SGTOP        
        IF (X.GT.1) S(X) = S(X-1)                 
        ITR(X) = 0            
   20   ITR(X) = ITR(X) + 1             
        CALL STRANS(CTG,S(X),Q(X),W(X),N,D(X),SETV(X),VISK)      
        R = CTG - CT(X)                 
        AR = DABS(R)          
C       WRITE (1,21) X,ITR(X),S(X),CTG            
C  21   FORMAT (' ',I4,'  ITER ',I2,'  S =',F12.10,'  CT =',F16.5)   
        IF (AR.LT.CON) GO TO 100        
        IF (ITR(X).GT.ITMAX) GO TO 26             
        INC = INCR            
        GO TO 24              
   22   INC = INC*10.0D0                
   24   SINC = S(X) + INC               
C       WRITE (1,25) SINC,INC           
C  25   FORMAT (' SINC ',F15.10,6X,'INC ',F15.10)           
        CALL STRANS(CTINC,SINC,Q(X),W(X),N,D(X),SETV(X),VISK)    
        DCTDS = (CTINC-CTG)/(SINC-S(X))           
        IF (DCTDS.EQ.0.0D0) GO TO 22              
        S(X) = S(X) - R/DCTDS          
        IF (S(X).LE.0.0D0) S(X) = .0000001D0        
        GO TO 20              
   26   WRITE (1,27) LOC(X)
   27    FORMAT (' MAX. ITERATIONS AT LOCATION ',F10.3)
        WRITE (6,27) LOC(X)
	STOP
  100    CONTINUE
	DO 200 XX=1,L
	X=L-XX+1
	U(X)=N**(-.6)*(Q(X)/W(X))**.4*S(X)**.3
	DEEP(X)=Q(X)/(U(X)*W(X))
	IF (X.EQ.L) HE(L)=BASE+U(L)**2/19.6D0
	IF (X.LT.L) HE(X)=HE(X+1)+DELX*1000.0D0*(S(X)+S(X+1))/2.0D0
	H(X)=HE(X)-U(X)**2/19.6D0
	B(X)=H(X)-DEEP(X)
	FR(X)=U(X)/(DSQRT(G*DEEP(X)))
  200    CONTINUE
	WRITE (1,30)
 30    FORMAT (' ',7X,'LOC ',5X,' CT ',8X,'U',7X,           
     1       'Q ',9X,'QS',10X,'H',10X,'B',8X,               
     2       'DEPTH',6X,'SLOPE',8X,'W',9X,'FR',/)           
       DO 500 X=1,L           
        TIME = 0.0            
      IF (X.EQ.1) GO TO 32
      IF (X.EQ.L) GO TO 34
      QI(X) = (Q(X+1)-Q(X-1))/(DELX*2.0D0)
      QSI(X) = (QS(X+1)-QS(X-1))/(DELX*2.0D0)
      GO TO 36
  32  QI(1) = (Q(2)-Q(1))/DELX
      QSI(1) = (QS(2)-QS(1))/DELX
      GO TO 36
  34  QI(L) = (Q(L)-Q(L-1))/DELX
      QSI(L) = (QS(L)-QS(L-1))/DELX
  36  CONTINUE
	WRITE (1,40) LOC(X),CT(X),U(X),Q(X),QS(X),H(X),B(X),
     1    DEEP(X),S(X),W(X),FR(X)       
 40     FORMAT (F11.5,F11.3,2F11.5,F11.7,3F11.5,F11.7,2F11.5)         
	WRITE (22,270) RUN,LOC(X),H(X),B(X)
  270   FORMAT (I12,3F12.4)
 500   CONTINUE               
         WRITE (1,501)
 501   FORMAT ('0',4X,'LOC',6X,'ITR',5X,'HE',8X,'QI',8X,'QSI',12X,    
     1    'D',11X,'SETV',/)             
       DO 510 X=1,L           
         WRITE (1,502) LOC(X),ITR(X),HE(X),QI(X),QSI(X),D(X),SETV(X)
 502     FORMAT (F11.5,2X,I3,F12.5,F12.7,F12.9,2F12.7)      
 510   CONTINUE               
C                   
C OUTPUT TO TAPE DATA SET FOR USE IN PROGRAM "LPMOD" AS INITIAL       
C CONDITION FOR RUNS OF CHANNEL ADJUSTMENT OVER TIME        
C                   
      WRITE (23,409) TIME               
 409  FORMAT (F18.6)          
      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) (SETV(X),X=1,L)              
      WRITE (23,411) (QSI(X),X=1,L)                
      WRITE (23,410) (LOC(X),X=1,L)               
      WRITE (23,410) (ED(X),X=1,L)
      WRITE (23,410) (ED(X),X=1,L)                
 410  FORMAT (5F13.7)                   
  411    FORMAT (5F13.10)
 	CLOSE (UNIT=22)
	CLOSE (UNIT=23)
       WRITE (6,70)
  70   FORMAT ('0 .......CALCULATIONS COMPLETED........',/)
       WRITE (6,72)
  72   FORMAT (' RESULTS LISTED IN THE FILE OUTPUT.DAT',/)
      STOP
      END                     
       SUBROUTINE STRANS(CT,S,Q,W,N,D,SETV,VISK)
       REAL*8 QS,S,Q,W,N,D,SETV,USTAR,CRIVEL,VDIR,CT,V 
       VDIR=1.0D0
      IF (Q.LT.0.0D0) THEN
       VDIR=-1.0D0
       Q=-Q
       ENDIF
       IF (S.LT.0.0D0) S=-S
       USTAR=9.8D0**.5*(Q*N/W)**.3*(S)**.35
      IF (USTAR*D/VISK.GT.70.0D0) THEN
        CRIVEL=2.05D0*SETV
      ELSE
        CRIVEL=SETV*(2.5D0/(DLOG10(USTAR*D/VISK)-0.06D0)
     1          +0.66D0)
      ENDIF
      V=N**(-0.6)*(Q/W)**.4*(S)**.3
       IF (V.LT.CRIVEL) GO TO 600
       CT=(5.435D0-(0.286D0*DLOG10(SETV*D/VISK))
     2    -0.457D0*DLOG10(USTAR/SETV)+(1.799D0-0.409D0*DLOG10(SETV*D
     3    /VISK)-.314D0*DLOG10(USTAR/SETV))*DLOG10(V*S
     4    /SETV-CRIVEL*S/SETV))
        CT = 10.0D0**(CT)
       GO TO 800
  600  CT=0.0D0
 800   CONTINUE
C      QS=Q*2.65D0**(-1)*CT*VISK
C      QS=VDIR*QS
       END
