C   *MIDAS : MODEL INVESTIGATING DENSITY AND SIZE SORTING A
C   *NONUNIFORM, QUASI-UNSTEADY, MOVABLE BED, SINGLE CHANNEL 
C   *FLOW MODEL FOR HETEROGENEOUS SIZE-DENSITY MIXTURES
C   *
C   *PROGRAMMED BY:
C   *RUDY SLINGERLAND AND KOEN VOGEL, PENN STATE UNIVERSITY 
C   *JOHN S. BRIDGE, SUNY BINGHAMTON, NY,
C   * AND A. VAN NIEKERK, WITS UNIVERSITY, JOHANNESBURG, 
C   *IN COOPERATION WITH THE CHAMBER OF MINES OF SOUTH AFRICA
C   *
C   * VERS. 5.2 most recent update  January, 1996 by Ruth Robinson
C   *
C   *AMAIN90.F
C   *
C   *******************************************************************
C
      IMPLICIT NONE
      INCLUDE 'common.h'
C----------------------------------------------------------------------
C   INTERNAL VARIABLES
C----------------------------------------------------------------------
      REAL*8 AWP(20,3),MANN1,P7INT
      INTEGER*4 IRSTRT(20),NTSTAR
      INTEGER*4 NRUN,IRUN,I,J,JN,NTFIN
      CHARACTER*24 F1(20),F2(20),F3(20)
      CHARACTER*4 MINSORT
      CHARACTER*10 FILENAME
C----------------------------------------------------------------------
C**********************************************************************
C----------------------------------------------------------------------
C     FILE MANAGEMENT PORTION : WHICH FILES ARE OPENED AS INPUT
C     OR OUTPUT FILES
C----------------------------------------------------------------------
      NRUN=1
      IRUN=0
      SCOUNT=1
        OPEN(3,FILE='./files')
139   READ(3,1,END=148) F1(NRUN)
      READ(3,1) F2(NRUN)
      IF (F1(NRUN).EQ.F2(NRUN)) THEN
         IRSTRT(NRUN)=1
         READ(3,1) F2(NRUN)
      ENDIF
      READ(3,1) F3(NRUN)
      NRUN=NRUN+1
      GOTO 139
148   NRUN=NRUN-1
149   IRUN=IRUN+1
        OPEN(4,FILE='./'//F1(IRUN))
        OPEN(9,FILE='./'//F2(IRUN))
      DO 154 KR=1,24
         IF (F3(IRUN)(KR:KR).EQ.' ') GOTO 155
154   CONTINUE
155   KR=KR-1
      FOUT=F3(IRUN)(1:KR)//'a'//CHAR(0)
      OPEN(8,FILE='./'//FOUT)
      FOUT=F3(IRUN)(1:KR)//'b'//CHAR(0)
      OPEN(7,FILE='./'//FOUT)
      FOUT=F3(IRUN)(1:KR)//'c'//CHAR(0)
      OPEN(12,FILE='./'//FOUT)
      FOUT=F3(IRUN)(1:KR)//'d'//CHAR(0)
      OPEN(13,FILE='./'//FOUT)
      FOUT=F3(IRUN)(1:KR)//'e'
C----------------------------------------------------------------------
C     READ INPUT PARAMETERS
C          >BED SHEAR STRESS DISTRIBUTION AND DERIVATIVES. IF OPTURB IS
C               1, THEN BED SHEAR STRESS IS GAUSSIAN DISTRIBUTED,
C               OTHERWISE IT IS LOG-GAUSSIAN DISTRIBUTED.
C          >SEDIMENT TRANSPORT PARAMETERS. IF OPTHIDE IS 1 THEN NO
C               GRAIN HIDING IS ASSUMED.
C          >HIDING FUNCTION OPTION. IF OPTFUNC IS 1 THEN THE KOMAR FUNC.
C               IS USED OTHERWISE THE EGIAZAROFF FUNCTION IS USED.
C          >SALTATION HEIGHT OPTION. IF OPTSALT IS 1 THEN THE BRIDGE
C               FUNCTION IS USED OTHERWISE THE EINSTEIN FUNCTION
C               IS USED.
C          >GRAIN SIZE AND DENSITY DISTRIBUTION CONTROL PARAMETERS
C          >GRAIN DENSITY VALUES
C          >WEIGHT PROPORTION OF AVAILABLE BED MATERIAL IN EACH
C               SIZE-DENSITY FRACTION. IF OPTDIST= 1, THEN WEIGHTS
C               ARE CALCULATED ASSUMING A LOG-GAUSSIAN DISTRIBUTION.
C          >INITIAL BOUNDARY CONDITION. IF OPT*INP IS 1 THEN CLEAR WATER
C               INFLOW IS ASSUMED. IF OPT*INP IS 2 THEN AN EQUILIBRIUM
C               CONDITION IS ASSUMED, IF OPT*INP IS 3 THEN A USER-
C               DEFINED CONCENTRATION PROFILE IS ASSUMED. IF OPTERO
C               IS 1 THEN EROSION OR DEPOSITION IS ALLOWED AT THE
C               HEAD OF THE REACH, IF OPTERO IS 2 THEN NO EROSION OR
C               DEPOSITION IS ALLOWED AT THE HEAD OF THE REACH.
C----------------------------------------------------------------------
      READ(4,2) N,NK,OPTURB,COVAR
      READ(4,3) RLONG,DELX,DELT
      READ(4,4) RHO,G,VISC,MANN1
      READ(4,5) A,B,VOLCON,THETAC,ACTIVE
      READ(4,6) OPTHIDE,OPTFUNC,OPTSALT,OPTLNODE
      READ(4,7) CON1,CON2,PHI,PSI,BETA
      READ(4,8) OPTBINP,OPTSINP,OPTERO,OPTEL,OPTWID,OPTSUSP,OPTABRAD,
     1          OPTTECT,OPTDISC
      M=INT((RLONG/DELX)+0.1D0)+2
      IF(OPTEL.NE.1.AND.IRSTRT(IRUN).NE.1)THEN
        READ(4,9) ELEVTOP,SLOPE,REGF,GRAD,BRK
        BREAK=INT((BREAK/DELX)+0.1D0)+2.
         ELSEIF(OPTEL.EQ.1)THEN
           READ(4,10)FILENAME
           OPEN(17,FILE='./'//FILENAME)
           DO JN=2,M
             READ(17,12) ELEVM(JN)
           ENDDO
           ELEVM(M+1)=ELEVM(M)-(ELEVM(M-1)-ELEVM(M))
           ELEVM(1)=ELEVM(2)+(ELEVM(2)-ELEVM(3))
           CLOSE(17)
      ENDIF
      IF(OPTWID.NE.1.AND.IRSTRT(IRUN).NE.1)THEN
        READ(4,26)WIDTOP,DELWID,REGG
          ELSEIF(OPTWID.EQ.1)THEN
             READ(4,11)FILENAME
             OPEN(17,FILE='./'//FILENAME)
             DO JN=2,M
               READ(17,12) WIDM(JN)
             ENDDO
             WIDM(M+1)=WIDM(M)-(WIDM(M-1)-WIDM(M))
             WIDM(1)=WIDM(2)+(WIDM(2)-WIDM(3))
             CLOSE(17)
      ENDIF
      IF(IRSTRT(NRUN).NE.1)THEN
        CALL ELWIDCALC
      ENDIF
      READ(4,28)TECRATE,ZEROLOC,EFOLDX
      ZEROLOC=INT((ZEROLOC/DELX)+0.1D0)+2
      EFOLDX=INT((EFOLDX/DELX)+0.1D0)+2
      READ(4,14) NTIM
      IF(OPTDISC.NE.1)THEN
        READ(4,15) QDSMOD,QFAC,DELQ,OPTTRIB,QSTMOD,THALFB,RENEW,
     1              QSSMODSTRT,QSTMODSTRT 
         ELSEIF(OPTDISC.EQ.1)THEN
           READ(4,30)FILENAME
           OPEN(17,FILE='./'//FILENAME)
           READ(17,*) QTCOUNTER
           READ(17,29)(HLEVEL(NT),QT(NT),DATETIME(NT),NT=1,QTCOUNTER)
           CLOSE(17)
      ENDIF
        READ(4,16) (PRINT(J),J=1,7)
        READ(4,17) XPRINT,PXINTERVAL,TPRINT,PTINTERVAL,PLOTPRINT,
     #TSTART,TFILE
        P7INT=PTINTERVAL*PRINT(7)
        TPRINT7=TPRINT
        READ(4,18) NSIGMA
      DO 228 J=1,NSIGMA
        READ(4,19)OPTDIST(J),SIGMA(J),NSIZE(J),KK1(J)
228   CONTINUE
      DELT=DELT*3600.0D0
      VISKIN=VISC/RHO
C----------------------------------------------------------------------
C     READ IN VARIABLES FOR SIZE DISTRIBUTIONS
C----------------------------------------------------------------------
      IF(IRSTRT(IRUN).NE.1)THEN
C----------------------------------------------------------------------
C     IF OPTDIST=1 THEN A PHINORMAL DISTRIBUTION IS GENERATED FOR
C     MINERAL J.  MEAN, DIMIN AND STDEV MUST BE IN PHI UNITS
C----------------------------------------------------------------------
         READ(4,20)FILENAME
          DO 263 J=1,NSIGMA
            READ(4,21)FILENAME
             IF(OPTDIST(J).EQ.1)THEN
               READ(4,22)MEAN(J),STDEV(J),DIMIN(J),PROPTJ(J),BLWTTO(J)
c                 READ(4,27) AWPT(I,J)
               CALL GSDIST(J)
               ELSE
C----------------------------------------------------------------------
C     IF OPTDIST= 2, THEN USE A USER-SPECIFIED SIZE DISTRIBUTION
C     WITH ALL PARAMETERS READ IN USING METER UNITS. IN THE CASE
C     OF SPATIALLY VARYING SEDIMENT DISCHARGE WITH A DIFFERENT OR SAME
C     SIZE DISTRIBUTION FROM THE HEAD INPUT, THE DISTRIBUTION IS
C     READ IN HERE (AWPT). IF THERE IS NO TRIBUTARY EFFECT, AWPT=0.0
C     SET THICKNESS AND MASS DUE TO TRIBUTARY INFLUENCE TO ZERO BEFORE
C     TRIB SUBROUTINE
C----------------------------------------------------------------------
                 DO 254 I=1,NSIZE(J)
                  AWPT(I,J)=0.0
                  READ(4,22)AWP(I,J),DIMID(I,J),DIFEND(I,J),AWPT(I,J)
254              CONTINUE
                 READ(4,22)DIMIN(J),PROPTJ(J),BLWTTO(J)
                  DO 261 JN=1,M
                  DO 260 I=1,NSIZE(J)
                     AWPM(I,J,JN)=AWP(I,J)*PROPTJ(J)/100.
                     AWPOM(I,J,JN)=AWPM(I,J,JN)
                     TMSPRTHK(I,J,JN)=0.0
                     TPRTHICK(I,J,JN)=0.0
260               CONTINUE
261               CONTINUE
             ENDIF
263       CONTINUE
          IF(OPTTRIB.EQ.1)THEN
               CALL TRIB
          ENDIF
C----------------------------------------------------------------------
C     DEFINE D50 OF AVAILABLE SEDIMENT
C----------------------------------------------------------------------
          CALL DISCR
          CALL DFIFTY
C----------------------------------------------------------------------
C     CALCULATE TIME-LOOP PARAMETERS
C----------------------------------------------------------------------
          NTSTAR=0
          NTFIN=NTIM
      ENDIF
      READ(4,23)RUN
      READ(4,24)RUNNAME
      OPEN(18,file='./plot.'//RUNNAME)
      OPEN(20,file='./dch.'//RUNNAME)
      OPEN(21,file='./trans.'//RUNNAME)
      OPEN(22,file='./gsdtr.'//RUNNAME)
      OPEN(23,file='./shearst.'//RUNNAME)
      OPEN(24,file='./hydrograph.'//RUNNAME)
      OPEN(25,file='./orighydro.'//RUNNAME)
        DO 285 J=0,NSIGMA-1
          MINSORT='min'//CHAR(48+J+1)
          OPEN(14+J,file='./'//MINSORT//'.'
     #//RUNNAME)
285     CONTINUE
C----------------------------------------------------------------------
C     IF RESTART FILE THEN SKIP PRELIMINARIES
C----------------------------------------------------------------------
      IF(IRSTRT(IRUN).EQ.1)THEN
        CALL READRSTRT
        NTSTAR=NT+1
        NTFIN=NT+NTIM
      ENDIF
C----------------------------------------------------------------------
C     CALCULATE VALUES OF CRITICAL BED SHEAR STRESS, CRITICAL
C     SHEAR VELOCITY, AND SETTLING VELOCITY FOR EACH GRAIN FRACTION
C----------------------------------------------------------------------
      CALL ENTRAIN
      CALL SETTLE
C---------------------------------------------------------------------
C     CALCULATE AVDENS-- THIS IS THE MASS/UNIT VOLUME IN THE
C     ORIGINAL ACTIVE LAYER; WHICH INCLUDES HEAD INPUT AND NOW 
C     THIS STATEMENT INCLUDES AN INPUT FROM TRIBUTARIES -AWPM 
C     IS A SUMMATION OF AWPM AND AWPT IF OPTTRIB=1(CALCULATED EARLIER)
C---------------------------------------------------------------------
      DO 306 JN=1,M
         AVDENS(JN)=0.0D0
          DO 306 J=1,NSIGMA
          DO 306 I=1,NSIZE(J)
            AVDENS(JN)=AVDENS(JN)+(AWPM(I,J,JN)*SIGMA(J))
306       CONTINUE
C----------------------------------------------------------------------
C     BEGIN COMPUTATIONS OF BEDLOAD TRANSPORT AND BED EROSION AND
C     DEPOSITION, STARTING AT TIMESTEP NT = 0
C----------------------------------------------------------------------
      DO 406 NT=NTSTAR,NTFIN
C----------------------------------------------------------------------
C     CLEAR ALL TRANSPORT ARRAYS
C----------------------------------------------------------------------
        DO 320 JN=1,M
          MANN(JN)=MANN1
320     CONTINUE
C----------------------------------------------------------------------
C  DETERMINE ALL THE NECESSARY CHANNEL DATA, OR READ VALUES. THIS
C  SUBROUTINE WILL BE CHANGED, DEPENDING ON THE SITUATION BEING MODELLED
C----------------------------------------------------------------------
        CALL FLDTA
        CALL SVELA
C----------------------------------------------------------------------
C  CALCULATE MASS FLUXES FOR SUSPENSION AND BEDLOAD. IF MORE OF A SIZE-
C  FRACTION IS ERODED THAN IS AVAILABLE IN THE BED, THEN
C----------------------------------------------------------------------
        IF(NT.EQ.TPRINT)THEN
          IF(PRINT(1).EQ.1)THEN
            CALL WRFLDTA
          ENDIF
          IF(PRINT(6).EQ.1)THEN
            CALL PLOTFILES
          ENDIF
        ENDIF
C----------------------------------------------------------------------
C       BEGIN COMPUTATION OF BEDLOAD TRANSPORT
C       CALCULATE BED SHEAR STRESS AT EVERY NODE
C----------------------------------------------------------------------
        CALL TURB
C----------------------------------------------------------------------
C     CALCULATE WHICH GRAIN FRACTIONS ARE TO BE INCORPORATED IN BEDLOAD
C     TRANSPORT AT EVERY NODE, EITHER WITH (352) OR W/O HIDING EFFECTS
C----------------------------------------------------------------------
        GOTO(350,352)OPTHIDE
350     CALL DRANGE
        GOTO 354
352     CALL ENTRAINH
        CALL DRANGEH
354     CALL BEDLOAD
        IF(NT.EQ.TPRINT)THEN
	  IF(NT.EQ.TPRINT7)THEN
            IF(PRINT(2).EQ.1)THEN
              CALL TOPRINT
            ENDIF
          ENDIF
        ENDIF
C----------------------------------------------------------------------
C  BEGIN COMPUTATION OF SUSPENDED LOAD CONCENTRATION PROFILES
C  (CALCULATE OR GET THE BED LAYER THICKNESS, DETERMINE DELY)
C----------------------------------------------------------------------
C  GENERATE LOGARITHMIC VELOCITY DISTRIBUTION
C----------------------------------------------------------------------
      CALL LOGDIST
C----------------------------------------------------------------------
C  COMPUTE VERTICAL DIFFUSIVITIY PROFILE
C----------------------------------------------------------------------
      CALL EXVERT
C----------------------------------------------------------------------
C   CALCULATE THE ADDITION OF MASS TO THE BED DUE TO TRIBUTARY EFFECT
C----------------------------------------------------------------------
      IF(OPTTRIB.EQ.1.AND.NT.NE.0)THEN
       CALL TRIB
      ENDIF
        IF(NT.EQ.0)THEN
          CALL SUSP1
        ENDIF
      CALL SUSP
C----------------------------------------------------------------------
C  CALCULATE EROSION/DEPOSITION AND NEW WEIGHT PROPORTIONS IN ACTIVE LAYER
C----------------------------------------------------------------------
      CALL ERODEP
C----------------------------------------------------------------------
C  CALCULATE THE BED ELEVATION CHANGE DUE TO SUBSIDENCE OR UPLIFT
C----------------------------------------------------------------------
       IF(OPTTECT.NE.0)THEN
         CALL TECTONICS
       ENDIF
C----------------------------------------------------------------------
C     DEFINE D50 OF AVAILABLE SEDIMENT
C----------------------------------------------------------------------
      CALL DFIFTY
      CALL TABLE
C----------------------------------------------------------------------
C  WRITE OUT RESULTS TO THE .c AND .d FILES
C----------------------------------------------------------------------
       IF(NT.EQ.TPRINT)THEN
           IF((PRINT(3).EQ.1).OR.(PRINT(4).EQ.1))THEN
             CALL RESULT
           ENDIF
       ENDIF
C----------------------------------------------------------------------
C  PRINT RESULTS FOR RESTART FILE
C----------------------------------------------------------------------
       IF(PRINT(5).EQ.1)THEN
          IF(NT.EQ.TSTART)THEN
            IF(NT.GT.0)THEN
              CALL WRITERSTRT
              TSTART=TSTART+TFILE
            ENDIF
          ENDIF
       ENDIF
      IF(NT.EQ.TPRINT)THEN
        TPRINT=TPRINT+PTINTERVAL
      ENDIF
      IF(NT.EQ.TPRINT7)THEN
	TPRINT7=TPRINT7+P7INT
      ENDIF
406   CONTINUE
      IF (IRUN.LT.NRUN) THEN
          GOTO 149
      ENDIF
      STOP
1      FORMAT(A24)
2      FORMAT(3(59X,I2,/),59X,F6.4)
3      FORMAT(2(59X,F10.0,/),59X,F10.6)
4      FORMAT(59X,F10.0,/,59X,F5.3,/,59X,F8.6/,59X,F8.6)
5      FORMAT(59X,F6.3,4(/,59X,F8.6))
6      FORMAT(3(59X,I1,/),59X,I1)
7      FORMAT(59X,F6.3,4(/,59X,F7.4))
8      FORMAT(8(59X,I1,/),59X,I1)
9      FORMAT(59X,F8.3,3(/,59X,F8.6),/,59X,F10.0,/,/,/)
10     FORMAT(7(/),59X,A10)
11     FORMAT(5(/),59X,A10)
12     FORMAT(F8.3)
13     FORMAT(F8.0)
14     FORMAT(59X,I10)
15     FORMAT(59X,I1,/,2(59X,F8.4,/),2(59X,I1,/),3(59X,F8.0,/),
     1        59X,I8,3(/))
16     FORMAT(/,6(59X,I1,/),59X,I1)
17     FORMAT(3(59X,I8,/,59X,I10,/),59X,I10)
18     FORMAT(59X,I1,5(/))
19     FORMAT(15X,I1,8X,F10.4,5X,I2,5X,F10.9)
20     FORMAT(A8,15(/))
21     FORMAT(A8,/)
22     FORMAT(5(3X,F10.8))
23     FORMAT(/,1X,A72)
24     FORMAT(59X,A6,5(/))
25     FORMAT(1X,'TIMESTEP   ',I9)
26     FORMAT(59X,F8.3,2(/,59X,F8.6),/,/,/)
27     FORMAT(3X,F10.8)
28     FORMAT(59X,F10.9,/,59X,F10.0,/,59X,F10.0)
29     FORMAT (2X,F6.0,1X,F7.0,2X,A14)
30     FORMAT(11(/),59X,A10)
32     FORMAT(I12)
       END
C****************************************************************
C***********************SUBROUTINE BEDLOAD************************
C   *CALCULATES THE BEDLOAD TRANSPORT RATES AND WEIGHTS PER UNIT AREA
C   *FOR EACH SIZE-DENSITY FRACTION
C
C                   BEDLOAD TRANSPORT OF DIFFERENT SIZE
C                   AND FRACTIONS IS PROPORTIONED ACCORDING
C                   TO THE VOLUMES IN THE BED.
C
C      SUBROUTINE: BEDLOAD.F
C
C      VERS. 4.0  Update 15 October, 1990
C      VERS. 4.1  Update October, 1995
C******************************************************************
      SUBROUTINE BEDLOAD
C------------------------------------------------------------------
C            IMPLICIT STATEMENT
C------------------------------------------------------------------
      IMPLICIT NONE
C------------------------------------------------------------------
C            DEFINE VARIABLES
C------------------------------------------------------------------
C                  EXTERNAL AND MODIFIED VARIABLES
C------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
C                  INTERNAL VARIABLES
C----------------------------------------------------------------------
      REAL*8 VOL,PROPVOL,TMP,MINISIF,CONV
      INTEGER*4 J,I,K,JN,IS,IF
C----------------------------------------------------------------------
C     INTRINSIC FUNCTION DECLARATION
C----------------------------------------------------------------------
C     INITIALIZATION
C**********************************************************************
C----------------------------------------------------------------------
C     CALCULATE SEDIMENT INFLUX AT THE HEAD (BOUNDARY CONDITION)
C----------------------------------------------------------------------
      IF(NT.EQ.0)THEN
        IF(OPTBINP.EQ.1)THEN
          DO 81 J=1,NSIGMA
          DO 81 I=1,NSIZE(J)
            BLWTM(I,J,1)=0.0D0
            BLWTAM(I,J,1)=0.0D0
81        CONTINUE
          ELSEIF(OPTBINP.EQ.3)THEN
            DO 87 J=1,NSIGMA
C-------RAR changed format for blwto from input.1 on 12/30/93----------
C-------RAR made more changes to bedload feed rate input 1/23/95-------
C----------------------------------------------------------------------
C     BLWTO IS TOTAL FEED RATE AT HEAD. THIS IS SEPARATED INTO SIZE
C     DENSITY FRACTIONS BY MULTIPLYING IT WITH AWPM TO GET BLWTM(I,J,1)
C----------------------------------------------------------------------
            IF(OPTDIST(J).EQ.1)THEN
               READ(4,1) BLWTO(1,J)
            ELSE
               READ(4,1)(BLWTO(I,J),I=1,NSIZE(J))
            ENDIF 
            DO 87 I=1,NSIZE(J)
            IF(OPTDIST(J).EQ.1)THEN
               BLWTM(I,J,1)=BLWTO(1,J)*AWPM(I,J,1)
            ELSE
               BLWTM(I,J,1)=BLWTO(I,J)
            ENDIF
87          CONTINUE
            BCOUNT=0
         ENDIF
      ENDIF
C----------------------------------------------------------------------
C     TEMPORALLY VARYING SEDIMENT FEED AT SOME STARTING TIME QSTMODSTRT
C----------------------------------------------------------------------
C-----RAR modified temp. varying sed. discharge algorithmn so that the
C-----starting time of discharge decrease can begin within a run and system
C-----can reach new equilibrium with new values---7/29/93
      IF(QSTMOD.EQ.1.AND.NT.GE.QSTMODSTRT)THEN
          IF(BCOUNT.GE.INT(RENEW*3600.0D0/DELT))THEN
            BCOUNT=0
          ENDIF
	      IF(BCOUNT*DELT.LE.THALFB*3600)THEN
            DO 100 J=1,NSIGMA
            DO 100 I=1,NSIZE(J)
             CONV=LOG(0.5)
C-----RAR changed following line 7/8/93 by removing DELT from denominator--
             BLWTM(I,J,1)=BLWTO(I,J)*AWPM(I,J,1)*EXP(CONV*((BCOUNT*DELT)
     #         /(THALFB*3600.0D0)))
100         CONTINUE
	      ELSEIF(BCOUNT*DELT.GT.THALFB*3600)THEN
              DO 105 J=1,NSIGMA
              DO 105 I=1,NSIZE(J)
	            BLWTM(I,J,1)=BLWTO(I,J)*AWPM(I,J,1)*0.5
105	          CONTINUE
	      ELSEIF(NT.LT.QSTMODSTRT)THEN
	       BLWTM(I,J,1)=BLWTO(I,J)*AWPM(I,J,1)
          ENDIF
      ENDIF
      BCOUNT=BCOUNT+1
C----------------------------------------------------------------------
C     START BEDLOAD CALC'S
C----------------------------------------------------------------------
      DO 201 JN=2,M
        VOL = 0.0D0
C----------------------------------------------------------------------
C      CALCULATE TOTAL 'VOLUME' = A PROPORTIONALITY FACTOR NECESSARY
C      TO CALCULATE THE PROPORTION OF VOLUME A CERTAIN AWPM REPRESENTS
C----------------------------------------------------------------------
       DO 123 J=1,NSIGMA
       DO 123 I=1,NSIZE(J)
         VOL=VOL+(AWPM(I,J,JN)/SIGMA(J))
C----------------------------------------------------------------------
C      SET ALL VARIABLES TO ZERO
C----------------------------------------------------------------------
         BLWTM(I,J,JN)=0.0D0
         BLWTAM(I,J,JN)=0.0D0
         SLWTM(I,J,JN)=0.0D0
         SLWTAM(I,J,JN)=0.0D0
         BEDCON(I,J,JN)=0.0D0
123   CONTINUE
C----------------------------------------------------------------------
C     BEDLOAD TRANSPORT RATE CALCULATIONS
C----------------------------------------------------------------------
      DO 200 K=1,NK
      DO 199 J=1,NSIGMA
         IS=IMIN(K,J,JN)
         IF=IMAX(K,J,JN)
         IF(IF.GT.NSIZE(J))GOTO 199
         IF(IS.GT.NSIZE(J))GOTO 156
C----------------------------------------------------------------------
C       IS=MINIMUM SIZE THAT WILL STAY IN THE MOVING BED LAYER
C       IF=MAXIMUM ENTRAINABLE SIZE
C----------------------------------------------------------------------
         DO 155 I=IS,IF
           IF(TO(K,JN).GT.TC(I,J,JN))THEN
             PROPVOL=AWPM(I,J,JN)/(SIGMA(J)*VOL)
             TMP=PROPT(K)*PROPVOL*10.0D0*(TO(K,JN)-TC(I,J,JN))
     #       *(SHVEL(K,JN)-SHVELC(I,J,JN))
C----------------------------------------------------------------------
C     CONVERSION OF "WET" TRANSPORT RATE TO "DRY" TRANSPORT RATE
C     BLWTM(I,J,JN)= THE DRY WEIGHT TRANSPORTED THROUGH A METER
C            WIDTH OF STREAM IN ONE SECOND AS BEDLOAD TRANSPORT
C----------------------------------------------------------------------
             TMP=TMP*SIGMA(J)/((SIGMA(J)-RHO)*G)
             BLWTM(I,J,JN)=BLWTM(I,J,JN)+TMP
C----------------------------------------------------------------------
C     BEDLOAD WEIGHT TRANSPORT PER UNIT AREA CALCULATIONS
C----------------------------------------------------------------------
             BLWTAM(I,J,JN)=BLWTAM(I,J,JN)+TMP/
     #       (A*(SHVEL(K,JN)-SHVELC(I,J,JN)))
           ENDIF
155      CONTINUE
156      MINISIF=MIN((IS-1),IF)
         DO 198 I=1,MINISIF
           IF(TO(K,JN).GT.TC(I,J,JN))THEN
             PROPVOL=AWPM(I,J,JN)/(SIGMA(J)*VOL)
C----------------------------------------------------------------------
C       CALCULATE THE IMMERSED WT FLUX.
C       [TMP]=KG/S3=N/MS
C----------------------------------------------------------------------
             TMP=PROPT(K)*PROPVOL*10.00D0*(TO(K,JN)-TC(I,J,JN))
     #       *(SHVEL(K,JN)-SHVELC(I,J,JN))
C----------------------------------------------------------------------
C      CALCULATE THE MASS FLUX PER UNIT WIDTH
C               [TMP]=KG/MS
C               [SLWTM]=KG/MS
C----------------------------------------------------------------------
             TMP=TMP*SIGMA(J)/((SIGMA(J)-RHO)*G)
             SLWTM(I,J,JN)=SLWTM(I,J,JN)+TMP
C----------------------------------------------------------------------
C     CALCULATE THE DRY MASS OF GRAINS MOVING OVER UNIT BED AREA
C     [TMP]=KG/M2 (THE TMP HERE REFERS TO WHAT IS BELOW THIS COMMENT
C---------------------------------------------------------------------
             TMP=TMP/(A*(SHVEL(K,JN)-SHVELC(I,J,JN)))
             SLWTAM(I,J,JN)=SLWTAM(I,J,JN)+TMP
C----------------------------------------------------------------------
C       IF OPTSALT.EQ.1 THEN USE BRIDGE FUNCTION TO CALCULATE
C       BEDCON=(KG/M3)=MASS PER UNIT VOLUME.
C    CALCULATE THE SUSPENDABLE AMOUNT PRESENT IN THE MOVING BED
C    LAYER:SLWTM [KG M-1 S-1] AND SLWTAM [KG M-2]
C       (SEE BRIDGE AND DOMINIC (1984) FOR YMAX EQUATION)
C----------------------------------------------------------------------
      IF(OPTSALT.EQ.1)THEN
          YMAX(JN)=(DIMID(I,J)*(2.53D0*(DSQRT((TO(K,JN)-TC(I,J,JN))
     #           /((SIGMA(J)-RHO)*G*DIMID(I,J))))+0.5D0))
       ELSE
C----------------------------------------------------------------------
C     ELSE OPTSALT.EQ.2 AND YMAX IS CALCULATED USING EINSTEIN(2*D50)
C----------------------------------------------------------------------
          YMAX(JN)=D50M(JN)*2.0D0
       ENDIF
       BEDCON(I,J,JN)=BEDCON(I,J,JN)+TMP/YMAX(JN)
       ENDIF
       YMAX(1)=YMAX(2)
198   CONTINUE
199   CONTINUE
200   CONTINUE
201   CONTINUE
      RETURN
1     FORMAT(1X,F5.0,2X,11(F5.0,2X))
c     FORMAT(3X,F4.3)
      END
C***************************SUBROUTINE DFIFTY******************
C       CALCULATES D50 OF THE SEDIMENT
C
C       SUBROUTINE DFIFTY.F
C
C       VERS. 4.0  Update 15 October, 1990
C**************************************************************
      SUBROUTINE DFIFTY
C--------------------------------------------------------------
C       IMPLICIT STATEMENT
C--------------------------------------------------------------
      IMPLICIT NONE
C--------------------------------------------------------------
C               EXTERNAL AND MODIFIED VARIABLES
C--------------------------------------------------------------
       INCLUDE 'common.h'
C---------------------------------------------------------------
C               INTERNAL VARIABLES
C---------------------------------------------------------------
      REAL*8 VOLTOT(60),VOL,SUM,LOWLIM,LOWDEN
      INTEGER*4 J,I,K,JN
C****************************************************************
C----------------------------------------------------------------
C     CALCULATE THE TOTAL VOLUME PRESENT IN EACH SEQUENTIAL SIZE
C     INTERVAL
C----------------------------------------------------------------
      DO 88 JN=2,M
      DO 42 K=1,NTOT
42      VOLTOT(K)=0.0D0
        DO 64 J=1,NSIGMA
         LOWLIM=DIMIN(J)
         LOWDEN=DIMIN(J)
         K=1
         DO WHILE((DENDM(K)-DIMIN(J)).LT.1.0D-9)
           K=K+1
         END DO
C----------------------------------------------------------------
C     CALCULATE THE TOTAL VOLUME BELONGING TO A CERTAIN DENDM(K)
C     SIZE INTERVAL
C----------------------------------------------------------------
         DO 63 I=1,NSIZE(J)
            DO WHILE(((DENDM(K)-DIFEND(I,J)).LE.1.0D-9)
     #      .AND.(K.LE.NTOT))
                  VOLTOT(K)=VOLTOT(K)+AWPM(I,J,JN)/SIGMA(J)*
     #               (DENDM(K)-LOWDEN)/(DIFEND(I,J)-LOWLIM)
                  LOWDEN=DENDM(K)
                  K=K+1
            END DO
            LOWLIM=DIFEND(I,J)
63       CONTINUE
64    CONTINUE
C----------------------------------------------------------------
C  CALCULATE D50
C----------------------------------------------------------------
      VOL=0.0D0
      DO 70 K=1,NTOT
70      VOL=VOL+VOLTOT(K)
      SUM=0.0D0
      DO 87 K=1,NTOT
         SUM=SUM+VOLTOT(K)
         IF(SUM.GE.0.5*VOL)THEN
C----------------------------------------------------------------
C  LINEAR INTERPOLATION BETWEEN SIZES WHOSE PERCENTILES ARE
C  RESPECTIVELY JUST LESS THAN AND JUST GREATER THAN 50%
C----------------------------------------------------------------
          IF(K.EQ.1)THEN
            D50M(JN)=DENDM(1)
          ELSE
           D50M(JN)=DENDM(K-1)+(DENDM(K)-DENDM(K-1))/(VOLTOT(K))
     #         *(0.5D0*VOL-(SUM-VOLTOT(K)))
          ENDIF
            GOTO 88
        ENDIF
87    CONTINUE
88    CONTINUE
      D50M(M+1)=D50M(M)
      D50M(1)=D50M(2)
      RETURN
      END
C******************************************************************
C       THIS SUBROUTINE DISCRETISES THE PARTICLE SIZE RANGES THROUGH
C       ALL SIZE RANGES OF ALL THE COMBINED DENSITY FRACTIONS FOR
C       CALCULATION OF THE NEW D50'S
C
C       DISCR.F
C
C      VERS. 4.0 Update 15 October, 1990
C*******************************************************************
      SUBROUTINE DISCR
      IMPLICIT NONE
C-------------------------------------------------------------------
C       DEFINE VARIABLES
C-------------------------------------------------------------------
C      EXTERNAL AND MODIFIED VARIABLES
C-------------------------------------------------------------------
       INCLUDE 'common.h'
C-------------------------------------------------------------------
C       INTERNAL VARIABLES
C-------------------------------------------------------------------
      REAL*8 TEMP
      INTEGER*4 J,I,K1,K2,NRP,K
C*******************************************************************
C-------------------------------------------------------------------
C     GET ALL THE SIZE RANGES IN THE DENDM ARRAY
C-------------------------------------------------------------------
      NTOT=0
      DO 41 J=1,NSIGMA
         DO 41 I=1,NSIZE(J)
              NTOT=NTOT+1
              DENDM(NTOT)=DIFEND(I,J)
41     CONTINUE
C-------------------------------------------------------------------
C         DO A BUBBLE SORT TO HAVE DENDM IN ASCENDING ORDER.
C-------------------------------------------------------------------
      DO 52 K1=1,NTOT-1
         DO 51 K2=K1+1,NTOT
            IF(DENDM(K1).LE.DENDM(K2))GOTO 51
            TEMP=DENDM(K1)
            DENDM(K1)=DENDM(K2)
            DENDM(K2)=TEMP
51       CONTINUE
52    CONTINUE
C-------------------------------------------------------------------
C          REMOVE REPEAT SIZE FRACTIONS
C-------------------------------------------------------------------
      NRP=0
      DO 66 K=1,NTOT-1
         IF((DENDM(K-NRP+1)-DENDM(K-NRP)).LT.0.1E-7)THEN
            NTOT=NTOT-1
            NRP=NRP+1
            DO 63 K1=K-NRP+1,NTOT
               DENDM(K1)=DENDM(K1+1)
63          CONTINUE
            DENDM(NTOT+1)=0.0D0
         ENDIF
66    CONTINUE
      RETURN
      END
C**********************SUBROUTINE DRANGE******************************
C*********************************************************************
C      *CALCULATES THE GRAIN SIZE FRACTIONS THAT CAN BE TRANSPORTED AS
C      *BEDLOAD FOR EACH SHEAR STRESS RANGE, NO GRAIN HIDING
C
C      SUBROUTINE: DRANGE
C
C      VERS. 4.0 Update 15 October, 1990
C********************************************************************
      SUBROUTINE DRANGE
C-------------------------------------------------------------------
C            IMPLICIT STATEMENT
C--------------------------------------------------------------------
      IMPLICIT NONE
C--------------------------------------------------------------------
C            DECLARATION STATEMENTS--EXTERNAL VARIABLES
C-------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
C            DECLARATION STATEMENTS--INTERNAL VARIABLES
C----------------------------------------------------------------------
      REAL*8 REG,AA,BB,CC,DIMSV,DLOG10,DIMDL,FD,DFD,DIN,DI,SETVEL
      INTEGER*4 I,K,J,JN,COUNT 
C*********************************************************************
C----------------------------------------------------------------------
C     CALCULATE THE GRAIN DIAMETER AT ENTRAINMENT THRESHOLD (DMAX)
C     FOR A GIVEN GRAIN DENSITY USING BRIDGE(1981) REGRESSION EQNS
C----------------------------------------------------------------------
      DO 152 JN=2,M
        DO 95 K=1,NK
          REG = SHVEL(K,JN)*D50M(JN)/VISKIN
          DO 94 J=1,NSIGMA
           IF(REG.LE.1.0D0) THEN
C----------------------------------------------------------------------
C     HYDRAULICALLY SMOOTH BED
C----------------------------------------------------------------------
            DMAX(K,J,JN) = ((TO(K,JN)**1.15D0)/(0.1D0*(SIGMA(J)-RHO)
     #*G*(VISKIN**0.3D0)*(RHO**0.15D0)))**1.43D0
           ELSEIF(REG.LE.60.0D0)THEN
C---------------------------------------------------------------------
C     HYDRAULICALLY TRANSITIONAL BED: SOLVE EQUATION FOR DMAX USING
C     QUADRATIC EQUATION
C---------------------------------------------------------------------
            AA=+0.16841D0
            BB=AA*2.0D0*DLOG(SHVEL(K,JN)/VISKIN)+0.09512D0
            CC=-DLOG(TO(K,JN)/(G*(SIGMA(J)-RHO)))-2.26D0-
     #      0.90488D0*DLOG(SHVEL(K,JN)/VISKIN)+AA*(DLOG
     #      (SHVEL(K,JN)/VISKIN))**2
            IF(BB**2-4.0D0*AA*CC .LE. 0.0D0)THEN
              DMAX(K,J,JN)=((TO(K,JN)**1.15D0)/(0.1D0*(SIGMA(J)-
     #        RHO)*G*(VISKIN**0.3D0)*(RHO**0.15D0)))**1.43D0
            ELSE
              DMAX(K,J,JN)=DEXP((-BB+DSQRT(BB**2-4.0D0*AA*CC))
     #        /(2.0D0*AA))
            ENDIF
            ELSE
C---------------------------------------------------------------------
C     HYDRAULICALLY ROUGH BED
C---------------------------------------------------------------------
               DMAX(K,J,JN)=(TO(K,JN)/(THETAC*(SIGMA(J)-RHO)*G))
            ENDIF             
94    CONTINUE
95    CONTINUE
C----------------------------------------------------------------------
C     CALCULATE THE GRAIN DIAMETER AT THE SUSPENSION THRESHOLD, DMIN
C----------------------------------------------------------------------
      DO 132 K=1,NK
       SETVEL=B*SHVEL(K,JN)
       DO 131 J=1,NSIGMA
        DIMSV=RHO*SETVEL*SETVEL*SETVEL/((SIGMA(J)-RHO)*G*VISKIN)
C----------------------------------------------------------------------
C     FIRST APPROXIMATION TO GRAIN SIZE, DI
C----------------------------------------------------------------------
        IF(DIMSV.LE.1.0D0)THEN
          DI=DSQRT(18.0D0*SETVEL*VISC/((SIGMA(J)-RHO)*G))
        ELSE
          DI=0.3D0*SETVEL*SETVEL*RHO/((SIGMA(J)-RHO)*G)
        ENDIF
C----------------------------------------------------------------------
C     SOLVE DIETRICH EQUATION FOR GRAIN SIZE USING NEWTON-RAPHSON
C----------------------------------------------------------------------
        COUNT=0
        DIMDL=DLOG10(DI*DI*DI*G*(SIGMA(J)-RHO)/(RHO*VISKIN*VISKIN))
117     FD=-3.767D0+1.929D0*DIMDL-0.0981D0*DIMDL**(2)-0.00575D0*
     #  DIMDL**(3)+0.00056D0*DIMDL**(4)-DLOG10(DIMSV)
        DFD=(1.929D0-0.196D0*DIMDL-0.01725D0*DIMDL**(2)+
     #  0.00224D0*DIMDL**(3))
        COUNT=COUNT+1
        DIN=DIMDL-FD/DFD
        IF(DABS(DIN-DIMDL).GT.0.000001D0.OR.COUNT.GT.100)THEN
          DIMDL=DIN
          GO TO 117
        ELSE
          DIMDL=10.0D0**(DIMDL)
          DMIN(K,J,JN)=(DIMDL*RHO*VISKIN**(2)/((SIGMA(J)-RHO)*G))
     #    **(0.333D0)
         ENDIF
131     CONTINUE
132   CONTINUE
C----------------------------------------------------------------------
C     DEFINE GRAIN SIZE ARRAY INDICES FOR DMIN AND DMAX
C----------------------------------------------------------------------
      DO 151 K=1,NK
       DO 150 J=1,NSIGMA
        DO 147 I=2,NSIZE(J)
         IF(DMIN(K,J,JN).GE.DIMID(I-1,J))THEN
         IF(DMIN(K,J,JN).LT.DIMID(I,J))THEN
           IMIN(K,J,JN)=I
         ELSE
           IMIN(K,J,JN)=I+1
         ENDIF
         ENDIF
         IF(DMAX(K,J,JN).GE.DIMID(I,J))IMAX(K,J,JN)=I
147      CONTINUE
       IF(DMAX(K,J,JN).LT.DIMID(1,J))IMAX(K,J,JN)=NSIZE(J)+100
       IF(DMIN(K,J,JN).LT.DIMID(1,J))IMIN(K,J,JN)=1
150    CONTINUE
151   CONTINUE
152   CONTINUE
      RETURN
      END
C*********************************************************************
C*************************SUBROUTINE DRANGEH**************************
C*********************************************************************
C      *CALCULATES THE GRAIN SIZE FRACTIONS THAT CAN BE TRANSPORTED AS
C       BEDLOAD FOR EACH SHEAR STRESS RANGE WITH HIDING
C
C      VERS. 4.0 Update 15 October, 1990
C*********************************************************************
      SUBROUTINE DRANGEH
C---------------------------------------------------------------------
      IMPLICIT NONE
C---------------------------------------------------------------------
C                  DECLARATION STATEMENTS-EXTERNAL VARIABLES  COM BLOCK
C---------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
C                  DECLARATION STATEMENTS--INTERNAL VARIABLES
C----------------------------------------------------------------------
      REAL*8 DI,SUR,X1,X2,FMID,F,RTBIS,DX,SETVEL,DIMSV,DIMDL,FD
      REAL*8 DFD,DIN,XACC,FUN1,D50
      INTEGER*4 J,K,I,JN,COUNT 
C**********************************************************************
C----------------------------------------------------------------------
C     DEFINE EGIAZAROFF FUNCTION
C----------------------------------------------------------------------
      FUN1(DI,D50,SUR,CON2)=DLOG10(CON2*DI/D50)-DSQRT(DI/SUR)
C**********************************************************************
      DO 167 JN=2,M
        D50=D50M(JN)
C----------------------------------------------------------------------
C     DETERMINE THE HIDING FUNCTION TO BE USED
C----------------------------------------------------------------------
        GOTO(67,75)OPTFUNC
C----------------------------------------------------------------------
C     KOMAR EQUATION 6 HIDING FUNCTION
C----------------------------------------------------------------------
67      DO 69 J=1,NSIGMA
        DO 69 K=1,NK
69        DMAX(K,J,JN) = (TO(K,JN)*D50**(CON1)/(THET50(JN)*
     #    (SIGMA(J)-RHO)*G))**(1.0D0/(1.0D0+CON1))
          GOTO 114
C----------------------------------------------------------------------
C     EGIAZAROFF HIDING FUNCTION
C----------------------------------------------------------------------
75        XACC = 0.0000001D0
          X1 = (2.7183)**2*D50/CON2
          X2 = 80.0D0*D50
          DO 113 J=1,NSIGMA
          DO 112 K=1,NK
            SUR = TO(K,JN)/((SIGMA(J)-RHO)*G*(DLOG10(CON2))**2.
     #      *THET50(JN))
C----------------------------------------------------------------------
C     USE BISECTION TO FIND THE ROOT OF THE EGIAZAROFF FUNCTION
C----------------------------------------------------------------------
            FMID = FUN1(X2,D50,SUR,CON2)
            F=FUN1(X1,D50,SUR,CON2)
            IF(F*FMID.GE.0.0D0) THEN
C----------------------------------------------------------------------
C     IF NO ROOT EXISTS FOR THIS TO, THEN SET DMAX TO A VALUE THAT
C     WILL MAKE IMAX = NSIZE(J)+100
C----------------------------------------------------------------------
               DMAX(K,J,JN) = DIFEND(1,J)-1.0D0
               GOTO 112
            ENDIF
            IF(F.LT.0.0D0)THEN
              RTBIS=X1
              DX=X2-X1
            ELSE
              RTBIS=X2
              DX=X1-X2
            ENDIF
            COUNT=0
103         DX=DX*0.5D0
            COUNT=COUNT+1
            DI=RTBIS+DX
            FMID=FUN1(DI,D50,SUR,CON2)
            IF(FMID.LT.0.0D0)RTBIS=DI
            IF((DABS(DX).LT.XACC).OR.(FMID.EQ.0.0D0).OR.(COUNT
     #.GT.100)) GOTO 111
              GOTO 103
111           DMAX(K,J,JN) = DI
112   CONTINUE
113   CONTINUE
C----------------------------------------------------------------------
C     CALCULATE THE GRAIN DIAMETER AT THE SUSPENSION THRESHOLD
C----------------------------------------------------------------------
114   DO 148 K=1,NK
       SETVEL=B*SHVEL(K,JN)
       DO 148 J=1,NSIGMA
         DIMSV=RHO*SETVEL*SETVEL*SETVEL/((SIGMA(J)-RHO)*G*VISKIN)
C----------------------------------------------------------------------
C     FIRST APPROXIMATION
C----------------------------------------------------------------------
         IF(DIMSV.LE.1.0D0)THEN
          DI=DSQRT(18.0D0*SETVEL*VISC/(SIGMA(J)-RHO)*G)
         ELSE
          DI=0.3D0*SETVEL*SETVEL*RHO/((SIGMA(J)-RHO)*G)
         ENDIF
C----------------------------------------------------------------------
C     SOLVE DIETRICH EQUATION FOR GRAIN SIZE USING NEWTON-RAPHSON
C----------------------------------------------------------------------
         COUNT=0
         DIMDL=DLOG10(DI*DI*DI*G*(SIGMA(J)-RHO)/(RHO*VISKIN*VISKIN))
134      FD=-3.767D0+1.929D0*DIMDL-0.098D0*DIMDL**(2)-0.00575D0*
     #   DIMDL**(3)+0.00056D0*DIMDL**(4)-DLOG10(DIMSV)
         COUNT=COUNT+1
         DFD=(1.929D0-0.196D0*DIMDL-0.01725D0*DIMDL**(2)+0.00224D0*
     #   DIMDL**(3))
         DIN=DIMDL-FD/DFD
         IF(DABS(DIN-DIMDL).GE.0.000001D0.OR.COUNT.GT.100)THEN
            DIMDL=DIN
            GO TO 134
         ELSE
            DIMDL=10.0D0**(DIMDL)
            DMIN(K,J,JN)=(DIMDL*RHO*VISKIN**(2)/((SIGMA(J)-RHO)*G))
     #      **(0.333D0)
         ENDIF
148   CONTINUE
C----------------------------------------------------------------------
C     DEFINE GRAIN SIZE ARRAY INDICES FOR DMIN AND DMAX
C----------------------------------------------------------------------
      DO 166 K=1,NK
        DO 166 J=1,NSIGMA
        DO 163 I=2,NSIZE(J)
         IF(DMIN(K,J,JN).GE.DIMID(I-1,J))THEN
         IF(DMIN(K,J,JN).LT.DIMID(I,J))THEN
           IMIN(K,J,JN)=I
         ELSE
           IMIN(K,J,JN)=I+1
         ENDIF
         ENDIF
         IF(DMAX(K,J,JN).GE.DIMID(I,J))IMAX(K,J,JN)=I
163   CONTINUE
      IF(DMAX(K,J,JN).LT.DIMID(1,J))IMAX(K,J,JN)=NSIZE(J)+100
      IF(DMIN(K,J,JN).LT.DIMID(1,J))IMIN(K,J,JN)=1
166   CONTINUE
167   CONTINUE
        RETURN
        END
C***********************************************************************
C*****************SUBROUTINE  ELWIDCALC*********************************
C***********************************************************************
C*   DESIGNED TO CALCULATE THE ELEVATIONS AND WIDTHS OF THE NODE POINTS
C*   FOLLOWING A USER SPECIFIED ALGORITHM
C*
C*   SUBROUTINE ELWIDCALC
C*
C*   VERS. 4.0  Update 15 October, 1990
C*   VERS. 4.2  Update June 1995 by Ruth Robinson
C***********************************************************************
      SUBROUTINE ELWIDCALC
      IMPLICIT NONE
C----------------------------------------------------------------------
C----------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
C----------------------------------------------------------------------
C     INTERNAL VARIABLES
C----------------------------------------------------------------------
      INTEGER*4 JN
C***********************************************************************
      DO 56 JN=M,1,-1
          IF(OPTWID.EQ.1)THEN
            GOTO 47
          ELSEIF(OPTWID.EQ.2)THEN
            WIDM(JN)=WIDTOP+DELWID*(JN-2)*DELX
          ELSEIF(OPTWID.EQ.3)THEN
            WIDM(JN)=WIDTOP*EXP((JN-2)*DELX*DELWID)
          ELSEIF(OPTWID.EQ.4)THEN
            WIDM(JN)=WIDTOP+LOG((JN-2)*DELX*DELWID)
          ELSEIF(OPTWID.EQ.5)THEN
            WIDM(JN)=WIDTOP*((FLOAT(JN-2)/FLOAT(M-2))**DELWID)
            IF(JN.EQ.2)WIDM(2)=WIDM(3)
          ENDIF
47        IF(OPTEL.EQ.1)THEN
            GOTO 56
          ELSEIF(OPTEL.EQ.2)THEN
            ELEVM(JN)=ELEVTOP-SLOPE*(JN-2)*DELX
c         ELSEIF(OPTEL.EQ.3)THEN
c            ELEVM(JN)=ELEVTOP-(EXP(DELX/1000.0*(JN-2)*SLOPE))*REGF*DELX
          ELSEIF(OPTEL.EQ.4)THEN
            ELEVM(JN)=ELEVTOP-LOG((JN-2)*DELX*SLOPE)
C---RAR added the next statement 11/24/93 to create a "pool" at downstream end---
          ELSEIF(OPTEL.EQ.5.AND.JN.GE.BREAK)THEN 
            ELEVM(BREAK)=ELEVTOP-SLOPE*(BREAK-2)*DELX
   	        ELEVM(JN)=ELEVM(BREAK)-GRAD*(JN-BREAK)*DELX
          ELSEIF(OPTEL.EQ.5.AND.JN.LT.BREAK)THEN
	        ELEVM(JN)=ELEVTOP-SLOPE*(JN-2)*DELX
          ELSEIF(OPTEL.EQ.6.AND.JN.GT.BREAK)THEN
            ELEVM(BREAK)=ELEVTOP-(EXP(DELX/1000.0*(BREAK-2)*SLOPE))
     1                *REGF*DELX
            ELEVM(JN)=ELEVM(BREAK)-(EXP(DELX/1000.0*(JN-BREAK)*GRAD))
     1                *REGF*DELX
           ELSEIF(OPTEL.EQ.6.AND.JN.LE.BREAK)THEN
	       ELEVM(JN)=ELEVTOP-(EXP(DELX/1000.0*(JN-2)*SLOPE))*REGF*DELX
	  ENDIF
56     CONTINUE
        WIDM(1)=WIDM(2)-(WIDM(3)-WIDM(2))
        WIDM(M+1)=WIDM(M)
        DO 61 JN=2,M
          ELEVOM(JN)=ELEVM(JN)
61      CONTINUE
	RETURN
        END
C**********************************************************************
C*******************************SUBROUTINE ENTRAIN*********************
C**********************************************************************
C    *CALCULATES THE CRITICAL SHEAR STRESS OF THE
C    *MEDIAN SIZE OF EACH SIZE-DENSITY FRACTION USING YALIN AND
C     KARAHAN
C
C      VERS. 4.0 Update 15 October, 1990
C**********************************************************************
      SUBROUTINE ENTRAIN
C----------------------------------------------------------------------
C      IMPLICIT STATEMENT
C----------------------------------------------------------------------
       IMPLICIT NONE
C----------------------------------------------------------------------
C      DECLARATION STATEMENTS--EXTERNAL VARIABLES       COM BLOCK
C----------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
C      DECLARATION STATEMENTS--INTERNAL VARIABLES
C---------------------------------------------------------------------
      REAL*8 TCI,FTC,DFTC,TCIN,REG,THETLN,RELN,TC1(20,3)
C----------------------INTEGERS--------------------------------------
        INTEGER*4 I,J
        INTEGER*4 JN,COUNT
C*********************************************************************
       DO 88 J=1,NSIGMA
        DO 88 I=1,NSIZE(J)
C----------------------------------------------------------------------
C     INITIALIZE APPROXIMATION TO CRITICAL BED SHEAR STRESS
C     USING KOMAR APPROXIMATION  EQUATIONS
C----------------------------------------------------------------------
          IF(DIMID(I,J).LE.0.0008D0)THEN
            TCI=0.1D0*17.7D0*(DIMID(I,J)*100.0D0)**0.59D0
          ELSE
            TCI=0.1D0*65.62D0*(DIMID(I,J)*100.0D0)**1.11D0
          ENDIF
          REG=SQRT(TCI/RHO)*DIMID(I,J)/VISKIN
          IF (REG.LE.1.0D0)THEN
C----------------------------------------------------------------------
C     VALUE FOR HYDRAULICALLY SMOOTH BED
C----------------------------------------------------------------------
              TC1(I,J)=0.1D0*(SIGMA(J)-RHO)*G*(RHO**0.15D0)*
     #        (VISKIN**0.3D0)
              TC1(I,J)=(TC1(I,J)*(DIMID(I,J)**0.7D0))**0.869565
           ELSEIF(REG.LE.60.0) THEN
C----------------------------------------------------------------------
C     VALUE FOR HYDRAULICALLY TRANSITIONAL BED
C----------------------------------------------------------------------
              COUNT=0
64            THETLN=DLOG(TCI/((SIGMA(J)-RHO)*G*DIMID(I,J)))
              COUNT=COUNT+1
              RELN=DLOG(SQRT(TCI/RHO)*DIMID(I,J)/VISKIN)
              FTC=-2.26D0-0.905D0*RELN+0.168D0*RELN*RELN-THETLN
              DFTC=(-1.45D0+0.168D0*RELN)/TCI
              TCIN=TCI-FTC/DFTC
              IF(COUNT.GT.100)GOTO 74
              IF(ABS(TCIN-TCI)-0.001D0) 74,74,72
72               TCI=TCIN
               GO TO 64
74              TC1(I,J)=TCIN
              ELSE
C----------------------------------------------------------------------
C     VALUE FOR HYDRAULICALLY ROUGH BED
C----------------------------------------------------------------------
                  TC1(I,J)=THETAC*G*(SIGMA(J)-RHO)*DIMID(I,J)
              ENDIF
C----------------------------------------------------------------------
C     VALUE OF CRITICAL SHEAR VELOCITY
C----------------------------------------------------------------------
       DO 87 JN=2,M
         TC(I,J,JN)=TC1(I,J)
         SHVELC(I,J,JN)=SQRT(TC1(I,J)/RHO)
87     CONTINUE
88    CONTINUE
      RETURN
      END
C**********************************************************************
C************************SUBROUTINE ENTRAINH***************************
C**********************************************************************
C      *CALCULATES THE CRITICAL SHIELDS THETA FOR THE MEDIAN SIZE OF THE
C       DISTRIBUTION AND THEN CALCULATES THE CRITICAL SHEAR STRESS FOR
C       THE I,J FRACTION USING A HIDING FUNCTION
C
C      SUBROUTINE ENTRAINH
C
C      VERS. 4.0  Update 15 October, 1990
C**********************************************************************
      SUBROUTINE ENTRAINH
C----------------------------------------------------------------------
C            IMPLICIT STATEMENT
C----------------------------------------------------------------------
       IMPLICIT NONE
C----------------------------------------------------------------------
C            DECLARATION STATEMENTS--EXTERNAL VARIABLES      COM. BLOCK
C----------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
C            DECLARATION STATEMENTS--INTERNAL VARIABLES
C----------------------------------------------------------------------
      REAL*8 TCI,REG,THETLN,RELN,D50
C----------------------------------------------------------------------
      INTEGER*4 I,J,JN
C**********************************************************************
C----------------------------------------------------------------------
C     INITIALIZE APPROXIMATION TO CRITICAL BED SHEAR STRESS
C----------------------------------------------------------------------
      DO 96 JN=2,M
         D50=D50M(JN)
         IF(D50.LE.0.0008D0)THEN
           TCI=0.1D0*17.7D0*(D50*100.0D0)**0.59D0
         ELSE
           TCI=0.1D0*65.62D0*(D50*100.0D0)**1.11D0
         ENDIF
         REG=DSQRT(TCI/RHO)*D50/VISKIN
         IF (REG.LE.1.0D0)THEN
C----------------------------------------------------------------------
C     VALUE FOR HYDRAULICALLY SMOOTH BED
C----------------------------------------------------------------------
            THETLN = DLOG(THETAC)
         ELSEIF(REG.LE.60.0) THEN
C----------------------------------------------------------------------
C     VALUE FOR HYDRAULICALLY TRANSITIONAL BED
C----------------------------------------------------------------------
            RELN=DLOG(REG)
            THETLN = -2.26D0-0.905D0*RELN+0.168*RELN*RELN
         ELSE
C----------------------------------------------------------------------
C     VALUE FOR HYDRAULICALLY ROUGH BED
C----------------------------------------------------------------------
            THETLN = DLOG(THETAC)
         ENDIF
         THET50(JN) = DEXP(THETLN)
C----------------------------------------------------------------------
C     CALCULATE THE TC FOR MIXTURES WITH HIDING,USING KOMAR EQN 6
C     OR EGIAZAROFF, DEPENDING UPON OPTFUNC
C----------------------------------------------------------------------
         DO 96 J=1,NSIGMA
         DO 96 I=1,NSIZE(J)
           GOTO(86,89)OPTFUNC
86         TC(I,J,JN)=THET50(JN)*G*(SIGMA(J)-RHO)*DIMID(I,J)
     #     **(1.0D0+CON1)*D50**(-CON1)
           GOTO 95
89         IF((DIMID(I,J)/D50).LE.1.0D0/CON2)THEN
               TC(I,J,JN)=100.0D0
           ELSE
               TC(I,J,JN)=(SIGMA(J)-RHO)*G*DIMID(I,J)*THET50(JN)*
     #        (DLOG10(CON2))**2./(DLOG10(CON2*DIMID(I,J)/D50))**2.
           ENDIF
95         SHVELC(I,J,JN)=SQRT(TC(I,J,JN)/RHO)
96       CONTINUE
97    CONTINUE
      RETURN
      END
C*******************************************************************
C**********************  SUBROUTINE ERODEP  *************************
C*******************************************************************
C  THIS SUBROUTINE CALCULATES THE TOTAL EROSION OR DEPOSITION OCCURRING
C  IN A TIME STEP AND CALCULATES THE NEW PARTICLE SIZE DISTRIBUTIONS IN
C  ACTIVE LAYER
C
C    SUBROUTINE ERODEP
C
C      VERS. 4.0  Update 15 October, 1990
C**********************************************************************
      SUBROUTINE ERODEP
C----------------------------------------------------------------------
C         IMPLICIT STATEMENT
C----------------------------------------------------------------------
      IMPLICIT NONE
C----------------------------------------------------------------------
C         DECLARATION STATEMENTS--EXTERNAL VARIABLES IN COMMON BLOCKS
C----------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
C            DECLARATION STATEMENTS--INTERNAL VARIABLES
C----------------------------------------------------------------------
      REAL*8 MSDEP,MSACTIVE,MSEROACT,TOTNEWMS,NEWMSACT(20,3),TEMPWEIT
      INTEGER*4 I,J,JN
C**********************************************************************
      DO 164 JN=2,M
C----------------------------------------------------------------------
C          INITIALIZE
C----------------------------------------------------------------------
C----------------------------------------------------------------------
C  CALCULATE THE TOTAL MASS/UNIT AREA IN THE ACTIVE LAYER
C  AT THE PREVIOUS TIMESTEP : MSACTIVE=KG/M2
C----------------------------------------------------------------------
      MSACTIVE=ACTIVE*VOLCON*NWAVDENS(JN)
C----------------------------------------------------------------------
C     MSPRTHK=MASS ADDED OR REMOVED OF EACH SIZE AND DENSITY FRACTION
C     UNITS= KG/M2
C     CALCULATE EROSION/DEPOSITION AND NEW BED COMPOSITION
C     MSDEP=SUM OF MASS/AREA DEPOSITED
C     MSPRTHK=MASS/AREA ERODED FOR MINERAL FRACTION IJ( >0 FOR EROSION)
C     MSACTIVE=TOTAL MASS/AREA IN THE ACTIVE LAYER
C----------------------------------------------------------------------
      MSDEP=0.0D0
      THICKM(JN)=0.0D0
      DO 70 J=1,NSIGMA
      DO 70 I=1,NSIZE(J)
         PRTHICK(I,J,JN)=MSPRTHK(I,J,JN)/(SIGMA(J)*VOLCON)
         THICKM(JN)=THICKM(JN)+PRTHICK(I,J,JN)
         IF(MSPRTHK(I,J,JN).LT.0.0D0)THEN
             MSDEP=MSDEP+MSPRTHK(I,J,JN)
         ENDIF
70    CONTINUE
      TOTNEWMS=0.0D0
C----------------------------------------------------------------------
C     REPLACE ERODED MATERIAL AT BASE OF ACTIVE LAYER
C----------------------------------------------------------------------
C     CALCULATE THE NEW LAYER THICKNESS AND MASS AFTER THE ERODIBLE 
C     FRACTIONS HAVE BEEN REMOVED.
C----------------------------------------------------------------------
C     IF THICKM(JN).LT.0 THEN THERE IS DEPOSITION,
C----------------------------------------------------------------------
      IF (THICKM(JN).LT.0.0D0) THEN
C----------------------------------------------------------------------
C     IF A THICKNESS GREATER THAN THE ACTIVE LAYER HAS BEEN DEPOSITED THEN
C     THE NEW ACTIVE LAYER TAKES ON THE SIZE DISTRIBUTION OF THE DEPOSITED
C     MATERIAL.
C----------------------------------------------------------------------
        DO 107 J=1,NSIGMA
        DO 107 I=1,NSIZE(J)
         IF(-THICKM(JN).GE.ACTIVE)THEN
           AWPM(I,J,JN)=MSPRTHK(I,J,JN)/MSDEP
         ELSE
C----------------------------------------------------------------------
C  CALCULATE THE NEW LAYER THICKNESS AND MASS AFTER THE ERODIBLE 
C  FRACTIONS HAVE BEEN REMOVED.
C----------------------------------------------------------------------
C     CALCULATE THE NEW MASS/AREA IN THE NEW ACTIVE LAYER
C     NEWMSACT(I,J)=KG/M2
C     TOTNEWMS=TOTAL NEW MASS/AREA ADDED OR REMOVED AT A NODE
C     TOTNEWMS=KG/M2
C     THICKM(JN)=TOTAL THICKNESS ERODED OR DEP. AT NODE (M)
C----------------------------------------------------------------------
         NEWMSACT(I,J)=(AWPM(I,J,JN)*MSACTIVE-MSPRTHK(I,J,JN))
           IF(NEWMSACT(I,J).LT.0.0)THEN
             NEWMSACT(I,J)=0.0
           ENDIF
           TOTNEWMS=TOTNEWMS+NEWMSACT(I,J)    
         ENDIF
107     CONTINUE
      ELSE
C----------------------------------------------------------------------
C         IF THICKM(JN).GT.0) THEN THERE IS EROSION
C----------------------------------------------------------------------
C     CALCULATE MASS PER UNIT AREA OF EACH GRAIN FRACTION REMAINING IN
C     ACTIVE LAYER FOLLOWING EROSION IN THIS TIME INCREMENT
C     ACTERO=REVISED MASS/AREA LEFT IN THE ACTIVE LAYER
C     THICK=THICKNESS ERODED FROM THE ACTIVE LAYER( >0 FOR EROSION )
C     NEWMSACT(I,J)=REVISED MASS/AREA OF IJ IN THE ACTIVE LAYER
C----------------------------------------------------------------------
      DO 143 J=1,NSIGMA
      DO 143 I=1,NSIZE(J)
C----------------------------------------------------------------------
C     CALCULATE THE NEW MASS/AREA IN THE NEW ACTIVE LAYER
C     NEWMSACT(I,J)=KG/M2
C     TOTNEWMS=TOTAL NEW MASS/AREA ADDED OR REMOVED AT A NODE
C     TOTNEWMS=KG/M2
C     THICKM(JN)=TOTAL THICKNESS ERODED OR DEP. AT NODE (M)
C     MSEROACT: IF EROSION OCCURS, THIS IS THE MASS/UNIT AREA OF THE 
C     PARENT MATERIAL
C----------------------------------------------------------------------
         NEWMSACT(I,J)=AWPM(I,J,JN)*MSACTIVE-MSPRTHK(I,J,JN)
         IF(NEWMSACT(I,J).LT.0.0)THEN
           NEWMSACT(I,J)=0.0
         ENDIF
         IF(ELEVOM(JN).GT.ELEVM(JN))THEN
           TEMPWEIT=AWPOM(I,J,JN)
           MSEROACT=ACTIVE*VOLCON*AVDENS(JN)
         ELSE
           MSEROACT=MSACTIVE
           TEMPWEIT=AWPM(I,J,JN)
         ENDIF
         NEWMSACT(I,J)=NEWMSACT(I,J)+TEMPWEIT*MSEROACT*
     #   THICKM(JN)/ACTIVE
         TOTNEWMS=TOTNEWMS+NEWMSACT(I,J)
143   CONTINUE
C----------------------------------------------------------------------
C RESET THE ELEVOM VALUE TO THE NEW BED ELEVATION
C----------------------------------------------------------------------
      IF(ELEVOM(JN).GE.ELEVM(JN))THEN
        ELEVOM(JN)=ELEVM(JN)-THICKM(JN)
      ENDIF
      ENDIF
      DO 159 J=1,NSIGMA
      DO 159 I=1,NSIZE(J)
        IF(-THICKM(JN).LE.ACTIVE)THEN
         AWPM(I,J,JN)=NEWMSACT(I,J)/TOTNEWMS
        ENDIF
        IF(AWPM(I,J,JN).LE.0.0D0)THEN
           AWPM(I,J,JN)=0.0D0
        ENDIF
159   CONTINUE
C----------------------------------------------------------------------
C       UPDATE THE ELEVATION OF THE STREAM PROFILE
C----------------------------------------------------------------------
       ELEVM(JN)=ELEVM(JN)-THICKM(JN)
164   CONTINUE
      ELEVM(M+1)=2.*ELEVM(M)-ELEVM(M-1)
      RETURN
      END
C***************************SUBROUTINE EXVERT***********************
C**********************************************************************
C       *COMPUTES VERTICAL DIFFUSIVITIES AT THE MIDPOINTS OF VERTICAL
C        INTERVALS.
C
C      EXVERT.F
C
C      VERS. 4.0   Update 15 October, 1990
C**********************************************************************
      SUBROUTINE EXVERT
C----------------------------------------------------------------------
C          DECLARATION STATEMENTS--EXTERNAL VARIABLES--NO COMMON BLOCKS
C----------------------------------------------------------------------
      IMPLICIT NONE
C----------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
C               DECLARATION STATEMENTS--INTERNAL VARIABLES
C----------------------------------------------------------------------
      REAL*8 Y
      INTEGER*4 L,JN
C**********************************************************************
      DO 36 JN=2,M
      DO 36 L=1,N-1
         Y=DCHM(JN)-DELY(JN)/2-(L-1)*DELY(JN)
         EPSY(L,JN)=0.4*SVELUM(JN)*Y*(1.0-Y/DCHM(JN))
36    CONTINUE
      RETURN
      END
C***************************SUBROUTINE FLDTA***************************
C**********************************************************************
C  THIS SUBROUTINE: 1) READS INFORMATION ABOUT THE CHANNEL HYDRAULICS
C                   2) CALCULATES THE STREAM PROFILE ON THE
C                      BASIS OF THE ENERGY EQUATION. A COMPLETE
C                      DERIVATION CAN BE FOUND ON PAGE 143 IN
C                      HENDERSON, F.M., 1966,      OPEN CHANNEL FLOW,
C                      MACMILLAN PUBLISHING CO., 520 P.
C
C      FLDTA.F
C
C      VERS. 5.0  Update December 1995 by Ruth Robinson
C**********************************************************************
      SUBROUTINE FLDTA
C----------------------------------------------------------------------
C            IMPLICIT STATEMENT
C----------------------------------------------------------------------
      IMPLICIT NONE
C----------------------------------------------------------------------
C          DECLARATION STATEMENTS--EXTERNAL VARIABLES      COMMON BLOCK
C----------------------------------------------------------------------
      INCLUDE 'common.h'
C----------------------------------------------------------------------
C            DECLARATION STATEMENTS--INTERNAL VARIABLES
C----------------------------------------------------------------------
      REAL*8 DCONT,AREA,P,V,SF1,H1,H2,H21,RCH,D2,A2,P2,V2,CONST,DELTAD
      REAL*8 RCH2,SF2,POS,FR,SLOP,YCRIT,CONV
      REAL*8 CONST2,DMIN1
      INTEGER*4 JN,K,END,STEP,NODE,SIGN,COUNT
      CHARACTER*1 DUMMY
**********************************************************************
      real*8 hl(100),qn2(100),dq,dh
      integer*4 year,month,day(1000),enter,interval
      integer*4 i,j,hour(1000),minute(1000),second(1000)
C**********************************************************************
C----------------------------------------------------------------------
C  DEFINE PARAMETERS
C  DCONT=DOWNSTREAM DISTANCE OF THE CONTROL SECTION (M)
C  QO= WATER DISCHARGE AT HEAD; DUR=DURATION OF QO
C  QDTMOD,THALFQ,QDTMODSTRT = PARAMETERS FOR TEMPORALLY MODIFIED DISCHARGE
C  OPTCONT=FLAG FOR METHOD OF DEPTH CALCULATION AT CONTROL SECTION
C----------------------------------------------------------------------
      IF(NT.EQ.0.AND.OPTDISC.EQ.0)THEN
        READ(9,3)DCONT
        ICONT=INT(DCONT/DELX+2.1D0)
        READ(9,3)ELEVM(ICONT)
        READ(9,2)DUMMY
      ENDIF
122   IF((NT.GE.NEXT.AND.OPTDISC.EQ.0).OR.
     1                                (NT.EQ.0.AND.OPTDISC.EQ.0))THEN
        READ(9,1)DINIT,QO,DUR,OPTCONT,QDTMOD,THALFQ,QDTMODSTRT
        QCOUNT=0
        DUR=DUR*3600.0D0
        IF(OPTCONT.EQ.2)THEN
          FLOWEL=DINIT+ELEVM(ICONT)
        ENDIF
        IF(NT.EQ.0)THEN
          NEXT=50e6
c       ELSE
c         NEXT=NEXT+INT(DUR/DELT+0.5)
        ENDIF
      ELSEIF(OPTDISC.EQ.1.AND.NT.EQ.0)THEN
C-------------------------------------------------------------------------
C  BEGIN READING DISCHARGE FILE. EACH TIMESTEP WILL HAVE A PARTICULAR QT
C  BY LINEARLY INTERPOLATING BETWEEN KNOWN DISCHARGES AND THEIR TIMES.
C  CONVERT HEAD WATER LEVEL FROM MM TO M. CONVERT DISCHARGE FROM LITRES/S
C  TO M^3/S
C  THIS ROUTINE LOOKS AT THE DATETIME HYDROGRAPH AND CONVERTS CHARACTER
C  TIME INTO NUMERICAL TIME AND PERFORMS THE INTERPOLATION AND PRINTS
C  HYDROGRAPH TO FILE
C-------------------------------------------------------------------------
        READ(9,3)DCONT
        ICONT=INT(DCONT/DELX+2.1D0)
        READ(9,3)ELEVM(ICONT)
        READ(9,2)DUMMY
        READ(9,1)DINIT,QO,DUR,OPTCONT,QDTMOD,THALFQ,QDTMODSTRT
        DUR=DUR*3600.0D0
        IF(OPTCONT.EQ.2)THEN
          FLOWEL=DINIT+ELEVM(ICONT)
        ENDIF
c-next temporarily hardwired (3/5/9) 
        NEXT=50e6 
        DO J=1,QTCOUNTER
          HLEVEL(J)=HLEVEL(J)/1000.
          QT(J)=QT(J)*0.001307
          year=(ichar(datetime(J)(1:1))-48)*10+(ichar(datetime(J)(2:2))
     1          -48)
          month=(ichar(datetime(J)(3:3))-48)*10+
     1           (ichar(datetime(J)(4:4))-48)
          day(j)=(ichar(datetime(J)(5:5))-48)*10+
     1           (ichar(datetime(J)(6:6))-48)
          hour(j)=(ichar(datetime(J)(9:9))-48)*10+
     1           (ichar(datetime(J)(10:10))-48)
          minute(j)=(ichar(datetime(J)(11:11))-48)*10+
     1              (ichar(datetime(J)(12:12))-48)
          second(j)=(ichar(datetime(J)(13:13))-48)*10+
     1              (ichar(datetime(J)(14:14))-48)
          WRITE(25,4) QT(J),YEAR+1900,MONTH,DAY(J),HOUR(J),MINUTE(J),
     1           SECOND(J)
        ENDDO
        interval=((day(2)*24.+hour(2)+minute(2)/60.+second(2)/3600.)-
     1  (day(1)*24.+hour(1)+minute(1)/60.+second(1)/3600.))/(delt/3600.)
        ENTER=INTERVAL
        QTCOUNTER=1
        CLOSE(25)
      ENDIF
C------------------------------------------------------------------------
C BEGIN INTERPOLATION
C------------------------------------------------------------------------
      IF(OPTDISC.EQ.1.AND.ENTER.EQ.INTERVAL)THEN
        DQ=(QT(QTCOUNTER+1)-QT(QTCOUNTER))/INTERVAL
        DH=(HLEVEL(QTCOUNTER+1)-HLEVEL(QTCOUNTER))/INTERVAL
        QN2(1)=QT(QTCOUNTER)
        HL(1)=HLEVEL(QTCOUNTER)
        WRITE(24,5) QN2(1),HL(1)
          DO J=1,INTERVAL-1
            QN2(J+1)=QT(QTCOUNTER)+DQ*J
            HL(J+1)=HLEVEL(QTCOUNTER)+DH*J
            WRITE(24,5) QN2(J+1),HL(J+1)
          ENDDO 
          QTCOUNTER=QTCOUNTER+1
          ENTER=2
          QN(2)=QN2(1) 
          GOTO 99
      ELSEIF(OPTDISC.EQ.1.AND.ENTER.NE.INTERVAL)THEN
        QN(2)=QN2(ENTER)
        ENTER=ENTER+1
      ENDIF
C----------------------------------------------------------------------
      ELEVM(1)=2.*ELEVM(2)-ELEVM(3)
C----------------------------------------------------------------------
C MODIFY QO IF IT IS A FUNCTION OF TIME (QDTMOD = 1), STARTING AT
C SPECIFIED TIME(TIME=QDTMODSTRT)
C----------------------------------------------------------------------
C-RAR modified temporally varying discharge algorithmn- it can begin varying
C-at some specified time during the run- 7/29/93
99     IF((QDTMOD.NE.1.OR.NT.LT.QDTMODSTRT).AND.OPTDISC.EQ.0)THEN
	     QN(2)=QO
       ENDIF
       IF((QDTMOD.EQ.1.AND.NT.GE.QDTMODSTRT).AND.OPTDISC.EQ.0)THEN
         IF(QCOUNT*DELT.LE.THALFQ*3600.0D0)THEN
           CONV=LOG(0.5)
           QN(2)=QO*EXP(CONV*(QCOUNT*DELT)/(THALFQ*3600.0D0))
            ELSEIF(QCOUNT*DELT.GT.THALFQ*3600.0D0)THEN
	      QN(2)=QO*0.5
         ENDIF
       ENDIF
       QCOUNT=QCOUNT+1
C----------------------------------------------------------------------
C      MODIFY QO IF IT IS A FUNCTION OF DISTANCE ( QDSMOD = 1 )
C----------------------------------------------------------------------
       IF(QDSMOD.EQ.1)THEN
         QN(ICONT)=(WIDM(ICONT)/DELQ)**(1./QFAC)
c        QN(ICONT)=QN(2)+QFAC*(ABS(ICONT+1)*DELX)
         ELSE
          QN(ICONT)=QN(2)
       ENDIF
C----------------------------------------------------------------------
C    DO UNIFORM FLOW CALCULATIONS IF OPTCONT=1
C    WARNING! BED SLOPE AT CONTROL SECTION MUST BE > 0 !!!
C----------------------------------------------------------------------
       IF(OPTCONT.EQ.1)THEN
          YCRIT=((QN(ICONT)/WIDM(ICONT))**2./9.81)**(1./3.)
C----------------------------------------------------------------------
C CALCULATE DINIT USING A NEWTON-RAPHSON SOLUTION FOR MANNING'S EQ.
C----------------------------------------------------------------------
          SLOP=(ELEVM(ICONT-1)-ELEVM(ICONT))/DELX
          IF(SLOP.LE.0)THEN
            SLOP=0.00001
          ENDIF
          DINIT=((QN(ICONT)*MANN(ICONT))/(WIDM(ICONT)*SQRT(SLOP)))
     #          **(3./5.)
          IF(DINIT.LT.YCRIT)THEN
            DINIT=YCRIT
          ENDIF
          IF((2.0*DINIT/WIDM(ICONT)).GT.0.00001)THEN
            COUNT=0
            CONST=QN(ICONT)*MANN(ICONT)/((WIDM(ICONT)**(5./3.))
     #            *SQRT(SLOP))
128         CONST2=(DINIT**(5./3.))/((WIDM(ICONT)+2.*DINIT)**(2./3.))
            DELTAD=ABS(CONST-CONST2)
            DELTAD=(DELTAD*(WIDM(ICONT)**(2./3.)))**(3./5.)
              IF((DELTAD.LT.1E-4).OR.(COUNT.EQ.20))GOTO 145
                  SIGN=((CONST-CONST2)/ABS(CONST-CONST2))
                  DINIT=((DINIT**(5./3.))+(SIGN*(DELTAD**(5./3.))))
     #                  **(3./5.)
                  COUNT=COUNT+1
                  GOTO 128
          ENDIF
       ELSEIF(OPTCONT.EQ.2)THEN
         DINIT=FLOWEL-ELEVM(ICONT)
       ENDIF
C----------------------------------------------------------------------
C       CALCULATE INITIAL CROSS-SECT AREA, THE WETTED
C       PERIMETER(P), THE MEAN VEL. (V), THE HYD RAD (RCH)
C       AND THE DEPTH OF THE WATER (DCHM).
C----------------------------------------------------------------------
145    AREA=DINIT*WIDM(ICONT)
       P=WIDM(ICONT)+2.0D0*DINIT
       V=QN(ICONT)/AREA
       RCH=AREA/P
       DCHM(ICONT)=DINIT
C----------------------------------------------------------------------
C      CALCULATE THE ENERGY SLOPE (SFL) AND THE HYDRAULIC HEAD
C      (H1), OR IF NOT AT CONTROL POINT TAKE THE PREVIOUS NODES
C      VALUES FOR THESE VALUES.
C      K=1 GO DOWNSTREAM FROM CONTROL SECTION; K=2 GO UPSTREAM
c----------------------------------------------------------------------
       DO 206 K=1,2
          IF(K.EQ.1)THEN
             END=M-1
             STEP=1
           ELSE
             END=3
             STEP=-1
           ENDIF
       DO 205 JN=ICONT,END,STEP
         COUNT=0
         IF(JN.EQ.ICONT)THEN
           SF1=((MANN(JN)*V)/(RCH**0.666667D0))**2.0D0
           H1=ELEVM(JN)+DINIT+V**2.0D0/(2.0D0*9.81D0)
           DMIN1=DINIT
         ELSE
           SF1=SF2
           H1=H2
         ENDIF
C----------------------------------------------------------------------
C     GUESS A NEW STREAM DEPTH AT THE NEXT NODE STARTING
C     WITH THE DEPTH AT THE PRESENT NODE. THE ONLY NEW DATA
C     USED IS THE STREAM WIDTH AT THE NEW NODE.
C----------------------------------------------------------------------
         D2=DMIN1
         NODE=JN+STEP
         IF(QDSMOD.EQ.1)THEN
C----------------------------------------------------------------------
C     QN IS A FUNCTION OF WIDTH FOLLOWING  RELATIONSHIP THAT
C     W=aQ^b. QFAC=b; DELQ=a BOTH CAN BE VARIED TO SUIT. RAR January 1995
C----------------------------------------------------------------------
           QN(NODE)=(WIDM(JN)/DELQ)**(1./QFAC)
         ELSEIF(QDSMOD.EQ.2)THEN
C----------------------------------------------------------------------
C     QN LINEARLY INCREASES DOWNSTREAM
C----------------------------------------------------------------------
           QN(ICONT)=QN(2)+QFAC*(ABS(ICONT+1)*DELX)
         ELSE
           QN(NODE)=QN(2)
         ENDIF
182      IF(D2.LE.1.0D-9)D2=1.0D-3
         A2=D2*WIDM(NODE)
         P2=WIDM(NODE)+2.*D2
         V2=QN(NODE)/A2
         RCH2=A2/P2
C----------------------------------------------------------------------
C     CALCULATE ENERGY BY ENERGY EQUATION AND BY ENERGY SLOPE
C     THEN RECALCULATE UNTIL THE DIFFERENCE IS SMALL
C----------------------------------------------------------------------
         SF2=((MANN(NODE)*V2)/(RCH2**0.666667D0))**2.0D0
         H2=ELEVM(NODE)+D2+V2**2.0D0/(2.0D0*9.81D0)
         H21=H1-STEP*0.5*DELX*(SF1+SF2)
         IF((ABS(H2-H21).GE.1E-6).AND.(COUNT.LE.50))THEN
            COUNT=COUNT+1
            D2=D2+(H21-H2)/(1-V2**2./(9.81*D2)+(3.*SF2*DELX)/(2.*RCH2))
            GOTO 182
         ELSE
            DCHM(NODE)=D2
         ENDIF
205      CONTINUE
206   CONTINUE
C----------------------------------------------------------------------
C      WRITE TO THE ARRAY STRDATA THE VALUES OF THE ENERGY SLOPE(SF)
C            THE HYDRAULIC RADIUS (RCH), THE POSITION(POS).
C----------------------------------------------------------------------
      DCHM(1)=DCHM(2)
      DCHM(M+1)=DCHM(M)
      QN(1)=QN(2)
C----------------------------------------------------------------------
C     STORE STREAM DATA IN ARRAYS
C----------------------------------------------------------------------
      DO 232 JN=1,M
        DCHM(JN)=INT((DCHM(JN)+0.0005D0)*1000.0D0)/1000.0D0
        AREA=DCHM(JN)*WIDM(JN)
        P=WIDM(JN)+2*DCHM(JN)
        V=QN(JN)/AREA
        RCH=AREA/P
        SF1=((MANN(JN)*V)/(RCH**0.666667D0))**2.0D0
        POS=FLOAT(JN-2)*DELX
        FR=V/(SQRT(G*DCHM(JN)))
        STRDATA(1,JN)=RCH
        STRDATA(2,JN)=FR
        STRDATA(3,JN)=SF1
        STRDATA(4,JN)=POS
        STRDATA(5,JN)=V
        STRDATA(6,JN)=QN(JN)
232   CONTINUE
      RETURN
1     FORMAT(3X,F8.3,2(3X,F8.0),2(3X,I2),2(3X,F8.0))
2     FORMAT(A1,9(/))
3     FORMAT(59X,F10.3)
4     FORMAT (2X,F10.0,5X,I4,2(1X,I2),3(1X,F6.0))
5     FORMAT(F10.4,5X,F10.4)
      END
C***********************SUBROUTINE GSDIST*************************
C*****************************************************************
C        *CALCULATES THE PHINORMAL GRAIN SIZE DISTRIBUTION FOR
C         EACH GRAIN DENSITY
C
C       GSDIST.F
C
C       VERS. 4.0  Update 15 October, 1990
C******************************************************************
       SUBROUTINE GSDIST(J)
C------------------------------------------------------------------
C        IMPLICIT STATEMENT
C------------------------------------------------------------------
       IMPLICIT NONE
C------------------------------------------------------------------
C        DECLARATION STATEMENTS--EXTERNAL VARIABLES--COMMON BLOCK
C-------------------------------------------------------------------
       INCLUDE 'common.h'
C---------------------------------------------------------------------
C        DECLARATION STATEMENTS--INTERNAL VARIABLES
C----------------------------------------------------------------------
       REAL*8 FO,FI,H,T,P,A1,A2,A3,A4,A5,H1,PISQR
       REAL*8 ALIMIT,BLIMIT,AWP(20,3),CONV,CLLENGTH,DISTMIN1,DIST
       INTEGER*4 I,J,JN,FLAG
C**********************************************************************
C     THIS SUBROUTINE CALCULATES THE ERROR FUNCTION USING THE
C     FOLLOWING ALGORITHM:
C     ERF(X)= 1-(A1*T+A2*T^2+A3*T^3+A4*T^4+A5*T^5)*EXP(-X^2)
C     T=1/(1+P*X)
C     PROBABILITIES ( I.E. AWP'S ) OF A CERTAIN GRAIN OCCURRING
C     IN A SIZE INTERVAL ARE CALCULATED BY SUBTRACTING ERF'S
C
C     CONV=CONVERSION FACTOR BETWEEN LOG BASE E AND LOG BASE 2
C     DETERMINE NECESSARY COEFFICIENTS                                -
C----------------------------------------------------------------------
      P=.2316419
      A1=.31938153
      A2=-0.356563782
      A3=1.781477937
      A4=-1.821255978
      A5=1.330274427
      CONV=0.693147181
      PISQR=SQRT(2.0*3.1415926536)
C----------------------------------------------------------------------
C     MEAN AND STANDARD DEVIATION SHOULD BE IN PHI UNITS
C     CLASSES RANGE FROM MEAN-3*STDEV TO MEAN+3*STDEV
C     ALIMIT=SMALLEST GRAIN SIZE (=LARGEST PHI SIZE )
C     BLIMIT=LARGEST GRAIN SIZE
C----------------------------------------------------------------------
      ALIMIT=MEAN(J)+3.0D0*STDEV(J)
      DIMIN(J)=EXP(-CONV*ALIMIT)/1000.0D0
      BLIMIT=MEAN(J)-3.0D0*STDEV(J)
      CLLENGTH=(BLIMIT-ALIMIT)/FLOAT(NSIZE(J))
      DISTMIN1=(ALIMIT-MEAN(J))/STDEV(J)
      H1=(EXP(-DISTMIN1*DISTMIN1/2.))/PISQR
      T=1/(1+P*DISTMIN1)
      H=((((A5*T+A4)*T+A3)*T+A2)*T+A1)*T
      FO=H*H1
      DO 99 I=1,NSIZE(J)
        FLAG=0
        DIST=((ALIMIT+I*CLLENGTH)-MEAN(J))/STDEV(J)
        DIMID(I,J)=EXP(-CONV*(ALIMIT+(I-0.5)*CLLENGTH))/1000.0D0
        DIFEND(I,J)=EXP(-CONV*(ALIMIT+I*CLLENGTH))/1000.0D0
        DIMID(I,J)=INT((DIMID(I,J)*1000000.)+0.5)/1000000.0D0
        DIFEND(I,J)=INT((DIFEND(I,J)*1000000.0D0)+0.5)/1000000.0D0
        H1=(EXP(-DIST*DIST/2.0D0))/PISQR
        IF(DIST.LT.0.)THEN 
          DIST=-DIST
          FLAG=1
        ENDIF
        T=1/(1+P*DIST)
        H=((((A5*T+A4)*T+A3)*T+A2)*T+A1)*T
        FI=H*H1
        IF(FLAG.EQ.1)THEN
           FI=1-FI
        ENDIF
        AWP(I,J)=FI-FO
        FO=FI
       DO 98 JN=1,M
        AWPOM(I,J,JN)=AWP(I,J)*PROPTJ(J)
        AWPM(I,J,JN)=AWP(I,J)*PROPTJ(J)
98    CONTINUE
99    CONTINUE
      RETURN
      END
C**********************************************************************
C**************************SUBROUTINE LOGDIST**************************
C**********************************************************************
C       *CALCULATES THE LOGRITHMIC VELOCITY PROFILE
C
C      LOGDIST.F
C
C      VERS. 4.0  Updated 15 October, 1990
C**********************************************************************
      SUBROUTINE LOGDIST
C----------------------------------------------------------------------
C      IMPLICIT STATEMENT
C----------------------------------------------------------------------
        IMPLICIT NONE
C----------------------------------------------------------------------
C      DECLARATION STATEMENTS--EXTERNAL VARIABLES--NO COMMON
C----------------------------------------------------------------------
       INCLUDE 'common.h'
C---------------------------------------------------------------------
C      DECLARATION STATEMENTS--INTERNAL VARIABLES
C---------------------------------------------------------------------
      REAL*8 Y,V
      INTEGER*4 JN,L
C**********************************************************************
      DO 55 JN=2,M
C----------------------------------------------------------------------
C   CONSIDER THE BED REFERENCE LEVEL TO BE AT EINSTEIN'S 2*D50
C   FOR THE PURPOSES OF CALCULATION.
C----------------------------------------------------------------------
      DELY(JN)=(DCHM(JN)-2.0*D50M(JN))/(N-1)
      YMAX(JN)=2.0*D50M(JN)
      DO 54 L=1,N
         V=STRDATA(5,JN)
         Y=DCHM(JN)-(L-1)*DELY(JN)
         IF(L.EQ.N)Y=YMAX(JN)
         VEL(L,JN)=(SVELUM(JN)/0.4)*(DLOG(Y/DCHM(JN))+1.)+V
C----------------------------------------------------------------------
C CHECK VALIDITY OF THIS ASSUMPTION
C----------------------------------------------------------------------
         IF(VEL(L,JN).LT.SVELUM(JN))VEL(N,JN)=SVELUM(JN)
54    CONTINUE
55    CONTINUE
      RETURN
      END
C*****************SUBROUTINE PLOTFILES*********************************
C**********************************************************************
C       WRITES DATA TO FILES TO BE USED IN GRAPHICS APPLICATIONS
C
C     PLOTFILES.F
C
C     VERS. 4.2  Update December, 1995
C**********************************************************************
      SUBROUTINE PLOTFILES
C----------------------------------------------------------------------
C      IMPLICIT STATEMENT
C----------------------------------------------------------------------
      IMPLICIT NONE
C----------------------------------------------------------------------
C      DECLARATION STATEMENTS--EXTERNAL VARIABLES--NO COMMON
C----------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
C     LOCAL VARIABLES
C----------------------------------------------------------------------
      REAL*8 POS
      REAL*8 TOTSBLWTA,TOTSCVELIA
      INTEGER*4 I,J,JN,NUMNODES
C**********************************************************************
      OPEN(19,FILE='tst')
      WRITE(19,3)NT
      CLOSE(19)
      IF(NT.EQ.0)THEN
        NUMNODES=INT((M-2)*DELX/PLOTPRINT)+1
        WRITE(18,*) NUMNODES
        WRITE(20,*) NUMNODES
        DO 59 J=1,NSIGMA
           WRITE(13+J,*) NUMNODES,NSIZE(J)
           DO 59 I=NSIZE(J),1,-1
             WRITE(13+J,*)DIFEND(I,J)*1000
59      CONTINUE
      ENDIF
      WRITE(20,*) NT
      WRITE(22,*) NT
      DO 65 J=1,NSIGMA
         WRITE(13+J,*) NT
65    CONTINUE
      TOTSBLWTA=0.0D0
      TOTSCVELIA=0.0D0
      WRITE(18,*)NT
      XPRINT=0
      DO 88 JN=2,M
      POS=(JN-2)*DELX
        IF(POS.EQ.XPRINT)THEN
          WRITE(18,2)POS,ELEVM(JN),D50M(JN)*1000
          WRITE(20,2)POS,ELEVM(JN),DCHM(JN)
          DO 86 J=1,NSIGMA
          WRITE(13+J,*)POS
          DO 86 I=NSIZE(J),1,-1
            WRITE(13+J,1)AWPM(I,J,JN)*100
86        CONTINUE
          ENDIF
	  XPRINT=XPRINT+PLOTPRINT
88     CONTINUE
       RETURN
1      FORMAT(1X,F7.3,2X,F4.2)
2      FORMAT(F10.2,4X,F10.4,4X,F9.4)
3      FORMAT(1X,'TIMESTEP   ',I9)
       END
C*******************************************************************************
C*    SUBROUTINE READRSTRT                                                     *
C*                                                                             *
C*    CREATED BY AVN AND KRV TO WRITE TO THE .e FILE ALL THE PARAMETERS        *
C*    NECESSARY TO DO A RESTART,  AT TIME NT=TFILE                             *
C*                                                                             *
C*    READRSTRT.F                                                              *
C*                                                                             *
C*    VERS. 4.0  Update 15 October, 1990                                       *
C*    MODIFIED VERS. 4.1 January, 1995  Ruth Robinson                          * 
C*    MODIFIED VERS. 4.2 December, 1995 Ruth Robinson                          * 
C*******************************************************************************
      SUBROUTINE READRSTRT
      IMPLICIT NONE
C-------------------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
      INTEGER*4 I,J,JN,K
C***********************************************************************
      READ(4,24)QDTMOD,THALFQ,QCOUNT
      READ(4,16)(WIDM(JN),JN=1,M+1)
      READ(4,16)(ELEVM(JN),JN=1,M+1)
      READ(4,16)(ELEVOM(JN),JN=1,M+1)
      READ(4,16)(D50M(JN),JN=1,M+1)
      DO 79 JN=1,M
      DO 79 J=1,NSIGMA
         READ(4,16)(AWPM(I,J,JN),I=1,NSIZE(J))
         READ(4,16)(AWPOM(I,J,JN),I=1,NSIZE(J))
79    CONTINUE
      READ(4,18) NT,NTOT,COUNTER
      READ(4,19)(DENDM(K),K=1,NTOT)
      READ(4,20)DINIT,QO,OPTCONT,ICONT,NEXT
      IF(OPTCONT.EQ.2)THEN
        READ(4,21)FLOWEL
      ENDIF
      READ(4,22)(DIMIN(J),J=1,NSIGMA)
      DO 92 J=1,NSIGMA
            READ(4,16) (AWPT(I,J),I=1,NSIZE(J))
            READ(4,16) (DIMID(I,J),I=1,NSIZE(J))
            READ(4,16) (DIFEND(I,J),I=1,NSIZE(J))
            READ(4,16) (BLWTM(I,J,1),I=1,NSIZE(J))
            READ(4,16) (BLWTT(I,J),I=1,NSIZE(J))
            READ(4,16) (BEDCON(I,J,1),I=1,NSIZE(J))
92    CONTINUE
16     FORMAT(100E14.6)
18     FORMAT(4(I9))
19     FORMAT(60E14.6)
20     FORMAT(2F8.3,I2,I5,I12)
21     FORMAT(F8.3)
22     FORMAT(9E14.6)
23     FORMAT(6E14.6)
24     FORMAT(I2,1X,F10.6,1X,I12)
       RETURN
       END
C*********************************************************************
C*********************** SUBROUTINE RESULT.F ************************
C*********************************************************************
C      *WRITES TO FILE RESULTA THE BEDLOAD TRANSPORT INFORMATION
C      *AFTER THE ACTIVE BED LAYER IS TAKEN INTO ACCOUNT
C
C      RESULT.F
C
C      VERS. 4.0  Update 15 October, 1990
C
C**********************************************************************
      SUBROUTINE RESULT
C----------------------------------------------------------------------
C            IMPLICIT STATEMENT
C----------------------------------------------------------------------
      IMPLICIT NONE
C----------------------------------------------------------------------
C          DECLARATION STATEMENTS--EXTERNAL VARIABLES      COMMON BLOCK
C----------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
C      DECLARATION STATEMENTS--EXTERNAL VARIABLES--NOT IN COMMON BLOCKS
C----------------------------------------------------------------------
      INTEGER*4 P
C----------------------------------------------------------------------
C            DECLARATION STATEMENTS--INTERNAL VARIABLES
C----------------------------------------------------------------------
      REAL*8 TOTSDBDX,TOTSBDT,TOTSXFLUX,TOTSTFLUX,TOTSBLWTA
      REAL*8 TOTBLWT,TOTSCVELIA,TOTSCVELI,TOTSPRTHICK
      REAL*8 TOTALTHICK,INSTTHIC
      INTEGER*4 I,J,JN
C*******************************************************************
C            IF AT FIRST TIMESTEP THEN OPEN THE FILES.
C---------------------------------------------------------------------
      IF(NT.EQ.0) THEN
          IF(PRINT(3).EQ.1)THEN
C----------------------------------------------------------------------
C                  WRITE HEADER TO SUMMARY DATA FILE (UNIT=12)
C----------------------------------------------------------------------
           WRITE(12,1) RUN,ICONT,DELX,DELT,PXINTERVAL,PTINTERVAL
1        FORMAT (1X, 'SEDIMENT TRANSPORT MODEL FOR HETEROGENEOUS SIZE',
     1                 ' - DENSITY MIXTURES',/
     1                 '                          VERSION 3.0',//
     1                 ' TITLE : ',A72,//
     1                 '------------------------------------------',/
     1                 ' THE NUMBER OF NODES (ICONT)  ',I4,/
     1                 ' THE DIST BETW. NODES(DELX)   ',F12.3, ' M',/
     1                 ' THE SIZE OF TIMESTEP (DELT)  ',F12.0, ' S',/
     1                 ' PRINT EVERY (PXINTERVAL)     ',I8,4X, ' M',/
     1                 ' PRINT EVERY (PTINTERVAL)     ',I8,4X,  'NT',/
     1                 '------------------------------------------',//)
          ENDIF
          IF(PRINT(4).EQ.1) THEN
C----------------------------------------------------------------------
C                  WRITE HEADER TO DETAILED DATA FILE (UNIT=13)
C----------------------------------------------------------------------
              WRITE(13,2) RUN,ICONT,DELX,DELT,PXINTERVAL,PTINTERVAL
2        FORMAT (1X, 'SEDIMENT TRANSPORT MODEL FOR HETEROGENEOUS SIZE',
     1                 ' - DENSITY MIXTURES',/
     1                 '                          VERSION 4.0',//
     1                 ' TITLE : ',A72,//
     1                 '------------------------------------------',/
     1                 ' THE NUMBER OF NODES (ICONT)  ',I4,/
     1                 ' THE DIST BETW. NODES(DELX)   ',F11.3, ' M',/
     1                 ' THE SIZE OF TIMESTEP (DELT)  ',F11.3, ' S',/
     1                 ' PRINT EVERY (PXINTERVAL)     ',I8,3X,   ' M',/
     1                 ' PRINT EVERY (PTINTERVAL)     ',I8,3X,   'NT',/
     1                 '------------------------------------------',//)
          ENDIF
      ENDIF
C----------------------------------------------------------------------
C            YOU HAVE NOW FINISHED WRITING THE HEADERS TO FILE 13 AND
C            FILE 12 IF THEY WERE TO BE WRITTEN TO.
C----------------------------------------------------------------------
      IF(PRINT(3).EQ.1) THEN
C----------------------------------------------------------------------
C                  WRITE MORE HEADER TO SUMMARY FILE
C----------------------------------------------------------------------
         WRITE(12,9)NT
         WRITE(12,3)
3                  FORMAT(/,'  POSITION     AVERAGE    ',
     1       '   AVERAGE      AVERAGE     AVERAGE  ',
     1       '    CUMULATIVE    CUMULATIVE    CUMULATIVE   DFIFTY  ')
         WRITE(12,4)
4               FORMAT('               BEDLOAD',
     1            '       SUSP. LOAD   BEDLOAD     SUSP. LOAD ',
     1            '  AMT.ERODED(-) AMT.ERODED(-) AMT.ERODED(-)       ')
         WRITE(12,5)
5               FORMAT('                SPATIAL',
     1            '       SPATIAL      MASS TRANS MASS TRANS ',
     1            '   OVER DELT     OVER TPRINT   OVER N*DELT     ')
         WRITE(12,6)
6               FORMAT('               GRADIENT',
     1            '      GRADIENT  ',
     1             '   RATE PER     RATE PER')
         WRITE(12,7)
7               FORMAT('              PER UNIT A',
     1            '    PER UNIT A    UNIT WIDTH    UNIT WIDTH')
         WRITE(12,8)
8                FORMAT('      (M)       (KG/M2S) ',
     1       '     (KG/M2S) ',
     1       '   (KG/MS)       (KG/MS)',
     1       '         (M)          (M)          (M)',
     1       '         (M)')
         WRITE(12,18)
      ENDIF
      IF(PRINT(4).EQ.1)THEN
	IF(NT.EQ.TPRINT7)THEN
            WRITE(13,19)
            WRITE(13,9)NT
9               FORMAT(1X,'AT TIMESTEP',I9,/)
        ENDIF
      ENDIF
      XPRINT=0
      DO 253 JN=2,M
	 IF((JN-2)*DELX.EQ.XPRINT)THEN
           IF(PRINT(4).EQ.1.AND.NT.EQ.TPRINT7)THEN
               WRITE(13,10)FLOAT(JN-2)*DELX
10               FORMAT(//'RESULTS AT ',F12.3,'METERS DOWN REACH')
            WRITE(13,11)
11                 FORMAT(/,'GRAIN SIZE AVERAGE   INSTANT',
     1  '   AVERAGE   INSTANT',
     1  '   AVERAGE   INSTANT   AVERAGE   INSTANT',
     1  '  CUMULATIVE CUMULATIVE CUMULATIVE')
C-----------------------------------------------------------------
            WRITE(13,12)
12              FORMAT('MIDPOINT   BEDLOAD   BEDLOAD',
     1            '  SUSP.LOAD SUSP.LOAD',
     1            '  BEDLOAD   BEDLOAD  SUSP.LOAD SUSP.LOAD',
     1            ' AMT.EROD(-)AMT.EROD(-)AMT.EROD(-)')
            WRITE(13,13)
13              FORMAT('           SPATIAL   SPATIAL',
     1            '   SPATIAL   SPATIAL',
     1            ' MASS TRANS MASS TRANS MASS TRANS MASS TRAN',
     1            'OVER DELT OVER PRINT OVER N*DELT')
            WRITE(13,14)
14              FORMAT('           GRADIENT  GRADIENT',
     1            ' GRADIENT  GRADIENT  ',
     1             'RATE PER   RATE PER   RATE PER  RATE PER',
     1       '    ')
            WRITE(13,15)
15              FORMAT('          /UNIT AREA /UNIT A  ',
     1            '/UNIT AREA /UNIT A  ',
     1            '/UNIT WIDTH /UNIT W   /UNIT W   /UNIT W    ')
            WRITE(13,16)
16               FORMAT('  (M)      (KG/M2S)  (KG/M2S) ',
     1            ' (KG/M2S)  (KG/M2S)  ',
     1            '(KG/MS)   (KG/MS)   (KG/MS)     (KG/MS)',
     1            '    (M)        (M)        (M)    ')
            WRITE(13,18)
C---------------------------------------------------------------------
C                  WRITE THE FLUX VALUES TO THE DETAILED FILE
C---------------------------------------------------------------------
              DO 193 J=1,NSIGMA
              DO 192 I=1,NSIZE(J)
                WRITE (13,17) DIMID(I,J),SDBDX(I,J,JN),SDBDT(I,J,JN),
     1          SXFLUX(I,J,JN),STFLUX(I,J,JN),SBLWTA(I,J,JN),
     2          SBLWT(I,J,JN),SCVELIA(I,J,JN),SCVELI(I,J,JN),
     3          -PRTHICK(I,J,JN),- SPRTHICK(I,J,JN),- TOTTHICK(I,J,JN)
192           CONTINUE
193           CONTINUE
           ENDIF
         ENDIF
17    FORMAT(E10.4,11(E10.4))
C----------------------------------------------------------------------
C            PREPARE TO SUMMARIZE THE FLUX VALUES FROM EACH SIZE
C            FRACTION. FIRST SET VALUES TO ZERO
C----------------------------------------------------------------------
         TOTSDBDX=0.0D0
         TOTSBDT=0.0D0
         TOTSXFLUX=0.0D0
         TOTSTFLUX=0.0D0
         TOTSBLWTA=0.0D0
         TOTBLWT=0.0D0
         TOTSCVELI=0.0D0
         TOTSCVELIA=0.0D0
         TOTSPRTHICK=0.0D0
         TOTALTHICK=0.0D0
         INSTTHIC=0.0D0
C----------------------------------------------------------------------
C            SUM THE FLUX VALUES OVER EACH SIZE FRACTION
C----------------------------------------------------------------------
         DO 229 J=1,NSIGMA
            DO 228 I=1,NSIZE(J)
               TOTSDBDX=TOTSDBDX+SDBDX(I,J,JN)
               TOTSBDT=TOTSBDT+SDBDT(I,J,JN)
               TOTSXFLUX=TOTSXFLUX+SXFLUX(I,J,JN)
               TOTSTFLUX=TOTSTFLUX+STFLUX(I,J,JN)
               TOTSBLWTA=TOTSBLWTA+SBLWTA(I,J,JN)
               TOTBLWT=TOTBLWT+SBLWT(I,J,JN)
               TOTSCVELIA=TOTSCVELIA+SCVELIA(I,J,JN)
               TOTSCVELI=TOTSCVELI+SCVELI(I,J,JN)
               TOTSPRTHICK=TOTSPRTHICK+SPRTHICK(I,J,JN)
               TOTALTHICK=TOTALTHICK+TOTTHICK(I,J,JN)
               INSTTHIC=INSTTHIC+PRTHICK(I,J,JN)
228          CONTINUE
229        CONTINUE
18    FORMAT(1X,'-------------------------------------------------------
     #----------------------------------------------------------------')
19              FORMAT(/,'************************************',
     1                  '*************************************',
     1                  '************************************',
     1                  '***************')
       IF(PRINT(4).EQ.1.AND.NT.EQ.TPRINT7.AND.(JN-2)*DELX.EQ.XPRINT)THEN
          WRITE(13,18)
          WRITE(13,21)  TOTSDBDX,TOTSBDT,TOTSXFLUX
     1     ,TOTSTFLUX,TOTSBLWTA,TOTBLWT,TOTSCVELIA,
     2      TOTSCVELI,- INSTTHIC,- TOTSPRTHICK,
     3      - TOTALTHICK
           ENDIF
      IF(PRINT(3).EQ.1.AND.(JN-2)*DELX.EQ.XPRINT)THEN
             WRITE(12,20) FLOAT(JN-2)*DELX,TOTSDBDX,
     1TOTSXFLUX,TOTSBLWTA,TOTSCVELIA,-INSTTHIC,
     2-TOTSPRTHICK,-TOTALTHICK,D50M(JN)
      ENDIF
      IF((JN-2)*DELX.EQ.XPRINT)THEN
      XPRINT=XPRINT+PXINTERVAL
      ENDIF
253   CONTINUE
C-----------------------------------------------------------------------
C     WRITE TO THE TRANSPORT FILES THE LAST NODE AVERAGE AND TOTAL 
C     BEDLOAD AND SUSPENDED LOAD TRANSPORT RATES 
C-----------------------------------------------------------------------
        DO 74 J=1,NSIGMA
        DO 73 I=1,NSIZE(J)
          WRITE(22,*)SBLWTA(I,J,M),SCVELIA(I,J,M)
73      CONTINUE
74      CONTINUE
      WRITE(21,*)NT,TOTSBLWTA,TOTSCVELIA
20        FORMAT(E9.4,1X,8(3X,E10.4))
21        FORMAT(3X,'TOTALS:',11(E10.4))
      RETURN
      END
C*********************SUBROUTINE SETTLE********************************
C**********************************************************************
C     *CALCULATES THE CONSTANT TERMINAL SETTLING VELOCITY OF EACH SIZE-
C      DENSITY FRACTION'S MEDIAN SIZE FROM DIETRICH'S EQUATION
C
C      SETTLE.F
C
C      VERS. 4.0  Update 15 October, 1990
C**********************************************************************
C*********************************************************************
      SUBROUTINE SETTLE
C----------------------------------------------------------------------
C      IMPLICIT STATEMENT
C----------------------------------------------------------------------
       IMPLICIT NONE
C----------------------------------------------------------------------
C      DECLARATION STATEMENTS--EXTERNAL VARIABLES--COMMON BLOCK
C----------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
C      DECLARATION STATEMENTS--INTERNAL VARIABLES
C----------------------------------------------------------------------
      REAL*8 DIMDL,RHS
      INTEGER*4 I,J
C**********************************************************************
      DO 47 J=1,NSIGMA
      DO 47 I=1,NSIZE(J)
        DIMDL=(SIGMA(J)-RHO)*G*(DIMID(I,J)**(3))/(RHO*VISKIN*VISKIN)
        DIMDL=DLOG10(DIMDL)
        RHS=-3.76715D0+1.92944D0*DIMDL-0.09815D0*DIMDL**(2)
        RHS=RHS-0.00575D0*DIMDL**(3)+0.00056D0*DIMDL**(4)
       SETV(I,J)=(((SIGMA(J)-RHO)*G*VISKIN*10.0D0**(RHS))/RHO)**(0.333)
47    CONTINUE
        RETURN
        END
C**********************************************************************
C**************************SUBROUTINE SUSP*****************************
C**********************************************************************
C      *THIS SUBROUTINE CALCULATES SUSPENDED LOAD TRANSPORT RATES
C      AND PERFORMS THE TRANSPORT REDUCTION
C
C      SUSP.F
C
C      VERS. 4.0  Update 15 October, 1990
C      VERS. 4.1  June 1994 by Ruth Robinson
C**********************************************************************
      SUBROUTINE SUSP
C----------------------------------------------------------------------
C            IMPLICIT STATEMENT
C----------------------------------------------------------------------
      IMPLICIT NONE
C----------------------------------------------------------------------
C          DECLARATION STATEMENTS--EXTERNAL VARIABLES IN COMMON BLOCKS
C----------------------------------------------------------------------
      INCLUDE 'common.h'
C----------------------------------------------------------------------
C----------------------------------------------------------------------
C          INTERNAL VARIABLES
C----------------------------------------------------------------------
      INTEGER*4 I,IN,J,JN
      REAL*8 B1,BCON,BWT,BWTA,SWT,SWTA,W
      REAL*8 MSACT(100)
      REAL*8 AV,Y,CA,Z,DY1,PHI1
      REAL*8 LOWLIM,HILIM,FACTOR
      INTEGER*4 IV,L,SIGNMSP,NEWSIGN,REDUC
C**********************************************************************
      DO 349 JN=2,M+1
C----------------------------------------------------------------------
C  CALCULATE THE AVERAGE DENSITY (AVDENS) (KG/M3) IN THE ACTIVE LAYER
C  THEN CALCULATE THE MASS/UNIT AREA IN THE ACTIVE LAYER
C  THIS IS MSACT(JN)=KG/M2
C----------------------------------------------------------------------
       IF(OPTTRIB.EQ.0)THEN
         NWAVDENS(JN)=0.0D0
         MSACT(JN)=0.0D0
         DO 83 J = 1,NSIGMA
         DO 83 I = 1,NSIZE(J)
           NWAVDENS(JN)=NWAVDENS(JN)+(AWPM(I,J,JN)*SIGMA(J))
83       CONTINUE
       ENDIF
         MSACT(JN)=ACTIVE*VOLCON*NWAVDENS(JN)
         DO 348 J=1,NSIGMA
         DO 348 I=1,NSIZE(J)
C----------------------------------------------------------------------
C  IF OPTSUSP IS > 0, CALCULATE PROFILES FOR SUSPENDED PARTICLES
C----------------------------------------------------------------------
c      IF(OPTSUSP.EQ.0)THEN 
c          DO 107 IV=1,N
c             CONC(IV)=0.0D0
c07        CONTINUE
c          CONIM(I,J,JN)=0.0D0
c          CVELM(I,J,JN)=0.0D0
c          GOTO 155
c      ENDIF 
C----------------------------------------------------------------------
C  IN JN=M+1 WE ARE AT THE LAST NODE AND CVELM & BLWTM ARE
C  SET AS BOUNDARY CONDITIONS
C-------------------------------------------------------------------
           REDUC=0
96         CA=BEDCON(I,J,JN)
           IF(JN.EQ.(M+1))GOTO 155
           W=SETV(I,J)
C-------------------------------------------------------------------
C  IF NO SUSPENDIBLE MATERIAL AVAILABLE SET SUSP LOAD=0.0
C-------------------------------------------------------------------
         IF(CA.EQ.0.0D0.AND.CONIM(I,J,JN-1).EQ.0.0D0)THEN
           DO 105 IV=1,N
              CONC(IV)=0.0D0
105        CONTINUE
           CONIM(I,J,JN)=0.0D0
           CVELM(I,J,JN)=0.0D0
           GOTO 155
         ELSE
C---------------------------------------------------------------------
C  CALCULATE THE CONCENTRATIONS AT THE GRID POINTS
C  IF OPTSUSP=1 THEN ANDRE VAN NIEKERK'S ALGORITHM
C  INCORPORATING LONGITUDINALLY ADVECTIVE SUSPENDED
C  LOAD FLUX IS USED
C---------------------------------------------------------------------
       IF(OPTSUSP.EQ.1)THEN
C---------------------------------------------------------------------
C  ANDRE VAN NIEKERK'S ALGORITHM SHOULD BE IN HERE
C  I OMITTED IT DUE TO COMPILATION DIFFICULTIES DUE TO THE LENGTH
C  OF  THE SUBROUTINE
C---------------------------------------------------------------------
         DY1=1.0
      ELSE
C----------------------------------------------------------------------
C  CALCULATE EQUILIBRIUM CONCENTRATIONS USING THE ROUSE EQUATIONS
C----------------------------------------------------------------------
        Z=W/(0.4*SVELUM(JN))
        AV=YMAX(JN)
        DO 128 L=1,N-1
          Y=DCHM(JN)-(L-1)*DELY(JN)
            IF(L.EQ.1)Y=DCHM(JN)-DELY(JN)/4
          CONC(L)=CA*(((DCHM(JN)-Y)*AV)/(Y*(DCHM(JN)-AV)))**Z
128       CONTINUE
          CONC(N)=CA
        ENDIF
        ENDIF
C----------------------------------------------------------------------
C  NOW THAT CONCENTRATIONS OF SUSPENDED MATERIAL AT VARIOUS
C  STREAM LEVELS ARE KNOWN (CONC(IV)), INTEGRATE THESE
C  TO GET THE SUSPENDED LOAD FLUX.
C  CALCULATE THE CONCENTRATION PROFILE INTEGRALS
C  BY CALCULATING THE CONCENTRATION INTEGRALS AT EACH GRID
C----------------------------------------------------------------------
      CONIM(I,J,JN)=0.0D0
      CVELM(I,J,JN)=0.0D0
      DO 150 IV=1,N-1
        CONIM(I,J,JN)=CONIM(I,J,JN)+(CONC(IV)+CONC(IV+1))*(DELY(JN)/2.)
        CVELM(I,J,JN)=CVELM(I,J,JN)+(CONC(IV)*VEL(IV,JN)+CONC(IV+1)*
     #VEL(IV+1,JN))*(DELY(JN)/2.)
150   CONTINUE
      CONIM(I,J,JN)=CONIM(I,J,JN)+(CONC(N)*YMAX(JN)/2.)
      IF(REDUC.EQ.0)THEN
       BLWTM(I,J,JN)=BLWTM(I,J,JN)+(CONC(N)*VEL(N,JN))*YMAX(JN)/2.
      ENDIF
155   IF(JN.EQ.2)THEN
        IF(OPTBINP.EQ.2)THEN
          BLWTM(I,J,1)=BLWTM(I,J,2)
        ENDIF
        IF(OPTSINP.EQ.2)THEN
          CVELM(I,J,1)=CVELM(I,J,2)
        ENDIF
      ENDIF
C----------------------------------------------------------------------
C  CONTINUITY PORTION OF MASSFLUX CALCS
C  CALL CONTINUITY EQUATION CALCULATION TO DETERMINE
C  NET MASS ADDED OR REMOVED AT EACH NODE
C  DETERMINE THE COEFFICIENTS IN THE FINITE DIFFERENCE FORMULATION OF
C  THE SEDIMENT CONTINUITY EQUATION. NOTE: A WEIGHTING FACTOR OF PHI
C  IS USED.
C----------------------------------------------------------------------
C            [XFLUX]=KG/M, [CVEL*]=KG/MS
C----------------------------------------------------------------------
       IN=JN-1
C--------------------------------------------------------------------
C  IF AT LAST NODE: SELECT BOUNDARY CONDITION USING OPTLNODE
C--------------------------------------------------------------------
       IF(IN.EQ.M)THEN
          PHI1=PHI
          FACTOR=1
C--------------------------------------------------------------------
C  IF OPTLNODE=0: CONDITION OF RIGID, UNERODABLE BED AT LAST
C  NODE: FLUX IN AT M = FLUX OUT AT M
C--------------------------------------------------------------------
          IF(REDUC.EQ.0)THEN
            IF(OPTLNODE.EQ.0)THEN
             DBDX(I,J,IN)=0.0
             XFLUX(I,J,IN)=0.0
             GOTO 251
C--------------------------------------------------------------------
C  IF OPTLNODE = 1 THEN USE A BACKWARD DIFFERENCE USING CALCULATED
C  TRANSPORT RATES, ON HALF THE DELX INTERVAL AT M
C--------------------------------------------------------------------
            ELSEIF(OPTLNODE.EQ.1)THEN
             PHI1=1.
             FACTOR=2.
C--------------------------------------------------------------------
C  IF OPTLNODE = 2 DO A CENTRAL DIFFERENCE UNDER THE ASSUMPTION
C  THAT AT NODE M+1 THE TRANSPORT RATES ARE = TO THOSE AT M
C--------------------------------------------------------------------
            ELSEIF(OPTLNODE.EQ.2)THEN
             BLWTM(I,J,M+1)=BLWTM(I,J,M)
             CVELM(I,J,M+1)=CVELM(I,J,M)
             PHI1=1.0
C--------------------------------------------------------------------
C  IF OPTLNODE=3 ASSUME TRANSPORT RATES AT M+1 ARE = 0
C--------------------------------------------------------------------
            ELSEIF(OPTLNODE.EQ.3)THEN
             BLWTM(I,J,M+1)=2.*BLWTM(I,J,M)-BLWTM(I,J,M-1)
             CVELM(I,J,M+1)=2.*CVELM(I,J,M)-CVELM(I,J,M-1)
            ELSEIF(OPTLNODE.EQ.4)THEN
             DBDX(I,J,IN)=2.0*DBDX(I,J,IN-1)-DBDX(I,J,IN-2)
             XFLUX(I,J,IN)=2.0*XFLUX(I,J,IN-1)-XFLUX(I,J,IN-2)
             GOTO 251
            ELSEIF(OPTLNODE.EQ.5)THEN
             DBDX(I,J,IN)=DBDX(I,J,IN-1)
             XFLUX(I,J,IN)=XFLUX(I,J,IN-1)
             GOTO 251
            ENDIF
            DBDX(I,J,IN)=((1.-PHI1)*((BLWTM(I,J,IN+1)*WIDM(IN))-
     1      (BLWTM(I,J,IN)*WIDM(IN)))+PHI1*FACTOR*((BLWTM(I,J,IN)
     2      *WIDM(IN))-(BLWTM(I,J,IN-1)*WIDM(IN-1))))*DELT/DELX
C
            XFLUX(I,J,IN)=((1.-PHI1)*((CVELM(I,J,IN+1)*WIDM(IN))-
     1      (CVELM(I,J,IN)*WIDM(IN)))+PHI1*FACTOR*((CVELM(I,J,IN)
     2      *WIDM(IN))-(CVELM(I,J,IN-1)*WIDM(IN-1))))*DELT/DELX
          ELSE
            MSPRTHK(I,J,IN)=TMSPRTHK(I,J,JN)+AWPM(I,J,IN)*MSACT(IN)
            PRTHICK(I,J,IN)=TPRTHICK(I,J,JN)+MSPRTHK(I,J,IN)/
     1      (SIGMA(J)*VOLCON)
            GOTO 345
          ENDIF
       ENDIF
C--------------------------------------------------------------------
C  IF AT FIRST NODE 
C--------------------------------------------------------------------
       IF(IN.EQ.1)THEN
            AWPM(I,J,1)=AWPM(I,J,2)
            MSACT(1)=MSACT(2)
            XFLUX(I,J,1)=((CVELM(I,J,IN+1)*WIDM(IN+1))-
     1    (CVELM(I,J,IN)*WIDM(IN)))*DELT/DELX
            DBDX(I,J,1)=((BLWTM(I,J,2)*WIDM(2))-
     1                (BLWTM(I,J,1)*WIDM(1)))*DELT/DELX
       ELSE 
            DBDX(I,J,IN)=((1.-PHI)*((BLWTM(I,J,IN+1)*WIDM(IN+1))-
     1     (BLWTM(I,J,IN)*WIDM(IN)))+PHI*((BLWTM(I,J,IN)*WIDM(IN))-
     2     (BLWTM(I,J,IN-1)*WIDM(IN-1))))*DELT/DELX
C
           XFLUX(I,J,IN)=((1.-PHI)*((CVELM(I,J,IN+1)*WIDM(IN+1))-
     1     (CVELM(I,J,IN)*WIDM(IN)))+PHI*((CVELM(I,J,IN)*WIDM(IN))-
     2     (CVELM(I,J,IN-1)*WIDM(IN-1))))*DELT/DELX
       ENDIF
251    TXFLUX(IN)=TXFLUX(IN)+XFLUX(I,J,IN)
       TDBDX(IN) = TDBDX(IN)+DBDX(I,J,IN)
C-----------------------------------------------------------------------
C     CALCULATE AMOUNT AND SIZE DENSITY DISTRIBUTION OF MATERIAL TO BE
C     ERODED OR DEPOSITED. THE FIRST TWO STATEMENTS NOW INCLUDE AN AMOUNT
C     THAT CORRESPONDS TO THE TRIBUTARY MATERIAL, CALCULATED IN TRIB. RAR
C     ADDED THIS 1/29/95
C-----------------------------------------------------------------------
      PRTHICK(I,J,IN)=TPRTHICK(I,J,JN)+(DBDX(I,J,IN)+XFLUX(I,J,IN))/
     #(VOLCON*SIGMA(J)*WIDM(IN))
      MSPRTHK(I,J,IN)=TMSPRTHK(I,J,JN)+(DBDX(I,J,IN)+XFLUX(I,J,IN))/
     #WIDM(IN)
      IF(ABS(MSPRTHK(I,J,IN)).LT.1E-7)GOTO 345
      IF((MSPRTHK(I,J,IN).GT.AWPM(I,J,IN)*MSACT(IN))
     #.OR.(REDUC.EQ.1))THEN
      IF(REDUC.EQ.0)THEN
        SIGNMSP=0
        IF(MSPRTHK(I,J,IN).NE.0.0)THEN
          SIGNMSP=MSPRTHK(I,J,IN)/ABS(MSPRTHK(I,J,IN))
          NEWSIGN=SIGNMSP
        ELSE
          SIGNMSP=1
          NEWSIGN=SIGNMSP 
        ENDIF
       ELSE
         IF(MSPRTHK(I,J,IN).NE.0.0)THEN
            NEWSIGN=MSPRTHK(I,J,IN)/ABS(MSPRTHK(I,J,IN))
         ELSE
            NEWSIGN=1
         ENDIF
       ENDIF
C-----------------------------------------------------------------------
C  IF AVAILABLE BED MATERIAL IS EXCEEDED THEN DO INTERVAL HALVING
C  TO DETERMINE THE ALLOWABLE TRANSPORT RATES AND MASS FLUXES
C  START BY HALFING THE FLUXES. REMEMBER BEDCON IS THE ACTIVE LAYER
C  THICKNESS.
C------------------------------------------------------------------------
      IF(REDUC.EQ.0)THEN
         B1=BLWTM(I,J,IN+1)+SLWTM(I,J,IN+1)
         BCON=BEDCON(I,J,IN+1)
         BWT=BLWTM(I,J,IN+1)
         BWTA=BLWTAM(I,J,IN+1)
         SWT=SLWTM(I,J,IN+1)
         SWTA=SLWTAM(I,J,IN+1)
         LOWLIM=0.0D0
         HILIM=1.0D0
         FACTOR=(HILIM-LOWLIM)/2.
         BEDCON(I,J,IN+1)=FACTOR*BCON
         BLWTM(I,J,IN+1)=FACTOR*BWT
         BLWTAM(I,J,IN+1)=FACTOR*BWTA
         SLWTM(I,J,IN+1)=FACTOR*SWT
         SLWTAM(I,J,IN+1)=FACTOR*SWTA
C----------------------------------------------------------------------
C  RECALCULATE THE SUSPENDED LOAD TRANSPORT RATES BASED
C  ON THE NEW BEDLOAD TRANSPORT RATE
C  SET REDUC=1 TO LET SUSP KNOW WE ARE IN THE MASS REDUCTION
C  LOOP
C----------------------------------------------------------------------
          REDUC=1
          GOTO 96
         ENDIF
C----------------------------------------------------------------------
C  IF MSPRTHK HAS CHANGED SIGN DUE TO HALVING OR IF YOU ARE
C  APPROACHING FROM BELOW,BUT HAVE NOT APPROACHED THE TOLERANCE 
C  THEN START HALVING THE DISTANCE UP TO PRTHICK.
C  SET THE LOW LIMIT TO THE PREVIOUS GUESS
C----------------------------------------------------------------------
      IF(((MSPRTHK(I,J,IN).NE.0.0).AND.(NEWSIGN.NE.SIGNMSP))
     1.OR.((AWPM(I,J,IN)*MSACT(IN)-MSPRTHK(I,J,IN)).GT.1.0E-12))THEN
               LOWLIM=FACTOR
               FACTOR=FACTOR+(HILIM-LOWLIM)/2.
               BEDCON(I,J,IN+1)=FACTOR*BCON
               BLWTM(I,J,IN+1)=FACTOR*BWT
               BLWTAM(I,J,IN+1)=FACTOR*BWTA
               SLWTM(I,J,IN+1)=FACTOR*SWT
               SLWTAM(I,J,IN+1)=FACTOR*SWTA
               GOTO 96
       ELSEIF (ABS(HILIM).LT.1.0D-12) THEN
            MSPRTHK(I,J,IN)=AWPM(I,J,IN)*MSACT(IN)
C----------------------------------------------------------------------
C  THE RIGHT SIDE OF THE IF CLAUSE IS:
C  AWP*AVERAGE DENSITY OF THE ACTIVE LAYER/DENSITY OF
C  SIZE FRACTION CONSIDERED
C----------------------------------------------------------------------
      ELSEIF((MSPRTHK(I,J,IN)-AWPM(I,J,IN)*MSACT(IN)).GT.1.0E-12)THEN
C----------------------------------------------------------------------
C  IF STILL TOO MUCH FLUX THEN START HALVING AGAIN
C  UNTIL YOU BEGIN TO APPROACH PRTHICK FROM BELOW
C  SET THE UPPER LIMIT TO THE PREVIOUS GUESS
C----------------------------------------------------------------------
               HILIM=FACTOR
               FACTOR=FACTOR-(HILIM-LOWLIM)/2.
               BEDCON(I,J,IN+1)=FACTOR*BCON
               BLWTM(I,J,IN+1)=FACTOR*BWT
               BLWTAM(I,J,IN+1)=FACTOR*BWTA
               SLWTM(I,J,IN+1)=FACTOR*SWT
               SLWTAM(I,J,IN+1)=FACTOR*SWTA
               GOTO 96
         ELSE
           MSPRTHK(I,J,JN)=MSACT(JN)*AWPM(I,J,JN)
         ENDIF
      ENDIF
345   DO 347 IV=1,N
         CONCM(I,J,IV)=CONC(IV)
347   CONTINUE
348   CONTINUE
349   CONTINUE
      RETURN
      END
C**************************SUBROUTINE SUSP1****************************
C**********************************************************************
C      *THIS SUBROUTINE CALCULATES SUSPENDED LOAD FEED RATES
C
C      SUSP1.F
C
C      VERS. 4.0   Update 15 October, 1990
C**********************************************************************
      SUBROUTINE SUSP1
C----------------------------------------------------------------------
C            IMPLICIT STATEMENT
C----------------------------------------------------------------------
      IMPLICIT NONE
C----------------------------------------------------------------------
C----------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
C       DECLARATION STATEMENTS--INTERNAL VARIABLES NOT IN COMMON BLOCKS
C----------------------------------------------------------------------
C----------------------------------------------------------------------
      INTEGER*4 I,J,IV,L,NN
      REAL*8 CA,W,Z,AV,Y
      CHARACTER*1 DUMMY
C**********************************************************************
C----------------------------------------------------------------------
C      THIS SUBROUTINE IS DESIGNED TO  CALCULATE THE BED CONCENTRATION 
C      FOR ALL SIZE-DENSITY FRACTIONS AT THE FIRST NODE (JN=1), THE
C      INPUT NODE. AS THIS NEVER CHANGES, IT ONLY NEEDS TO BE 
C      CALCULATED AT TIMESTEP 0
C      IF OPTSINP=2 HOWEVER, THE EQUILIBRIUM CASE, INPUT AT NODE 1
C      IS DYNAMIC
C----------------------------------------------------------------------
      IF(OPTSINP.EQ.3)THEN
        READ(4,2)DUMMY
        DO 63 J=1,NSIGMA
          READ(4,1)(BEDCON(I,J,1),I=1,NSIZE(J))
63      CONTINUE
      ENDIF
      DO 110 J=1,NSIGMA
      DO 110 I=1,NSIZE(J)
C----------------------------------------------------------------------
C            OPTSINP=1 : CLEAR WATER INPUT AT THE HEAD
C            IV IS A VERTICAL COUNTER
C       SINCE ALL VALUES ARE SET TO 0.0 AT THE START OF THE SUBROUTINE
C       THE INPUT IS CLEAR WATER UNLESS SUSP. CONC'S ARE CHANGED
C----------------------------------------------------------------------
        IF(OPTSINP.EQ.1)THEN
           BEDCON(I,J,1)=0.0D0
           CVELM(I,J,1)=0.0D0
           CONIM(I,J,1)=0.0D0
C----------------------------------------------------------------------
C       OPTSINP=2 : SUSPENDED LOAD CONCENTRATIONS AT THE HEAD ARE
C       IN EQUILIBRIUM WITH THE FLOW CONDITIONS AT THE HEAD
C       THIS IS DEALT WITH IN SUSP
C----------------------------------------------------------------------
C      IF OPTSINP=3 A USER DEFINED BED- AND SUSPENDED LOAD IS INPUT
C      AT THE HEAD. VALUES OF BEDCON ARE READ IN AT THE HEAD AND
C      AN EQUILIBRIUM PROFILE IS THEN CALCULATED
C----------------------------------------------------------------------
        ELSEIF(OPTSINP.EQ.3)THEN
           CA=BEDCON(I,J,1)
           W=SETV(I,J)
           Z=W/(0.4D0*SVELUM(2))
           AV=YMAX(2)
           NN=N-1
           DO 99 L=1,NN
C----------------------------------------------------------------------
C     CALCULTE EQUILIBRIUM CONCENTRATIONS USING ROUSE EQUATIONS
C----------------------------------------------------------------------
             Y=DCHM(2)-(L-1)*DELY(2)
             IF(L.EQ.1)Y=DCHM(2)-DELY(2)/4.0D0
             CONC(L)=CA*(((DCHM(2)-Y)*AV)/(Y*(DCHM(2)-AV)))**Z
99         CONTINUE
           CONC(N)=CA
           DO 106 IV=1,N-1
             CONIM(I,J,1)=CONIM(I,J,1)+(CONC(IV)+CONC(IV+1))
     #       *DELY(2)/2.0D0
             CVELM(I,J,1)=CVELM(I,J,1)+(CONC(IV)*VEL(IV,2)
     #       +CONC(IV+1)*VEL(IV+1,2))*DELY(2)/2.0D0
106        CONTINUE
           CONIM(I,J,1)=CONIM(I,J,1)+(CONC(N)*YMAX(2))/2.0D0
           BLWTM(I,J,1)=BLWTM(I,J,1)+(CONC(N)*VEL(N,2))*YMAX(2)/2.0D0
        ENDIF
110    CONTINUE
1      FORMAT(3X,12(F4.3,2X))
2      FORMAT(A1,/)
       RETURN
       END
C*****************************SVELA************************************
C**********************************************************************
C      *CALCULATES SHEAR VELOCITY AND SHEAR STRESS FOLLOWING EINSTEIN
C     
C      SVELA.F
C
C      VERS. 4.0  Update 15 October, 1990
C**********************************************************************
      SUBROUTINE SVELA
C----------------------------------------------------------------------
C            IMPLICIT STATEMENT
C----------------------------------------------------------------------
      IMPLICIT NONE
C----------------------------------------------------------------------
C            DECLARATION STATEMENTS--EXTERNAL VARIABLES
C----------------------------------------------------------------------
       INCLUDE 'common.h'
C---------------------------------------------------------------------- -
C            DECLARATION STATEMENT--INTERNAL VARIABLES
C----------------------------------------------------------------------
      REAL*8 AK,DR1,DR2,DR3,DR,F3,RCH,SF,V
      INTEGER*4 JN,COUNT
C**********************************************************************
C----------------------------------------------------------------------
C  CALCULATE INITIAL DATA (RECTANGULAR CHANNEL SECTION IS ASSUMED)
C----------------------------------------------------------------------
      DO 73 JN=2,M
        RCH=STRDATA(1,JN)
        V=STRDATA(5,JN)
        SF=STRDATA(3,JN)
C----------------------------------------------------------------------
C  CALCULATE REDUCED HYDRAULIC RADIUS
C----------------------------------------------------------------------
        COUNT=0
        AK=2.5*D50M(JN)
        DR1=2.0*RCH
        DR2=0.1*RCH
        DR3=0.5*(DR1+DR2)
52      F3=6.25+2.5*DLOG(DR3/AK)-V/SQRT(G*DR3*SF)
        COUNT=COUNT+1
        IF(ABS(F3).LE.0.001)GO TO 62
        IF(DR1.EQ.DR2.OR.COUNT.GT.50)GOTO 62
        IF(F3)57,62,59
57      DR2=DR3
        GO TO 60
59      DR1=DR3
60      DR3=0.5*(DR1+DR2)
        GO TO 52
62       DR=DR3
        IF(DR.GT.RCH)DR=RCH
C----------------------------------------------------------------------
C  CALCULATE SHEAR VELOCITY
C----------------------------------------------------------------------
        SVELM(JN)=SQRT(G*DR*SF)
        SVELUM(JN)=SQRT(G*RCH*SF)
C----------------------------------------------------------------------
C     CALCULATE THE BED SHEAR STRESS
C----------------------------------------------------------------------
        TOXM(JN)=RHO*G*SF*DR
73    CONTINUE
      RETURN
      END
C******************************SUBROUTINE TABLE************************
C**********************************************************************
C     TABULATES THE TRANSPORT RATE RESULTS
C
C     TABLE.F
C
C      VERS. 4.0  Update 15 October, 1990
C**********************************************************************
      SUBROUTINE TABLE
C----------------------------------------------------------------------
C            IMPLICIT STATEMENT
C----------------------------------------------------------------------
      IMPLICIT NONE
C----------------------------------------------------------------------
C           DECLARATION STATEMENTS--EXTERNAL VARIABLES IN COMMON BLOCKS
C----------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
C            DECLARATION STATEMENTS--EXTERNAL VARIABLES NO IN COMMON
C----------------------------------------------------------------------
C----------------------------------------------------------------------
C            DECLARATION STATEMENTS--      INTERNAL VARIABLES
C----------------------------------------------------------------------
      REAL*8 TIMEPRES,TIMEMIN1
      INTEGER*4 J,JN,I
C**********************************************************************
C----------------------------------------------------------------------
C          SDBDX=TIME AVERAGED MASS EXCHANGE BETWEEN BEDLOAD LAYER
C                  AND ACTIVE BED (KG/M2S)
C          SDBDT=INSTANTANEOUS MASS EXCHANGE RATE BETWEEN BEDLOAD LAYER
C                  AND ACTIVE BED (KG/M2S)
C          SXFLUX=TIME AVERAGED MASS EXCHANGE BETWEEN SUSPENDED LOAD
C                  AND BED  (KG/M2S)
C          STFLUX=INSTANTANEOUS MASS EXCHANGE RATE BETWEEN BEDLOAD
C                 AND ACTIVE BED
C          SBLWT=INSTANTANEOUS BEDLOAD FLUX (KG/MS)
C          SBLWTA=AVERAGE BEDLOAD FLUX
C          SCVEL=INSTANTANEOUS SUSPENDED LOAD FLUX
C          SCVELIA=AVERAGE SUSPENDED LOAD FLUX
C          SPRTHICK=THICKNESS ERODED OR DEPOSITED OVER THE TIMESTEP
C          TOTTHICK=TOTAL THICKNESSES ERODED OR DEPOSITED SINCE THE
C                START OF SIMULATION.
C----------------------------------------------------------------------
C----------------------------------------------------------------------
      DO 135 JN=2,M
      IF(NT.EQ.0) THEN
         COUNTER=1
         DO 91 J=1,NSIGMA
            DO 90 I=1,NSIZE(J)
               SDBDX(I,J,JN)=0.0D0
               SDBDT(I,J,JN)=0.0D0
               SXFLUX(I,J,JN)=0.0D0
               STFLUX(I,J,JN)=0.0D0
               SBLWT(I,J,JN)=0.0D0
               SCVELI(I,J,JN)=0.0D0
               SBLWTA(I,J,JN)=0.0D0
               SCVELIA(I,J,JN)=0.0D0
               SPRTHICK(I,J,JN)=0.0D0
               TOTTHICK(I,J,JN)=0.0D0
90          CONTINUE
91       CONTINUE
        ENDIF
C----------------------------------------------------------------------
C                  IF PREVIOUS TIMESTEP IS A TPRINT THEN CLEAR THE
C                  AVERAGING ARRAYS SO THAT A NEW TIME AVERAGE CAN
C                  BE DEVELOPED
C----------------------------------------------------------------------
      IF((NT-1).EQ.(TPRINT-PTINTERVAL)) THEN
          COUNTER=1
          DO 110 J=1,NSIGMA
             DO 109 I=1,NSIZE(J)
               SDBDX(I,J,JN)=0.0D0
               SXFLUX(I,J,JN)=0.0D0
               SBLWTA(I,J,JN)=0.0D0
               SCVELIA(I,J,JN)=0.0D0
               SPRTHICK(I,J,JN)=0.0D0
109          CONTINUE
110       CONTINUE
      ENDIF
C
      TIMEPRES=FLOAT(COUNTER)
      TIMEMIN1=FLOAT(COUNTER-1)
      DO 133 J=1,NSIGMA
         DO 132 I =1,NSIZE(J)
            SDBDX(I,J,JN)=(((SDBDX(I,J,JN)*TIMEMIN1
     1          *WIDM(JN)*DELT)+DBDX(I,J,JN))/(TIMEPRES*DELT))/WIDM(JN)
C
            SDBDT(I,J,JN)=(DBDX(I,J,JN)/DELT)/WIDM(JN)
            SXFLUX(I,J,JN)=(((SXFLUX(I,J,JN)*TIMEMIN1
     1      *WIDM(JN)*DELT)+XFLUX(I,J,JN))/(TIMEPRES*DELT))/WIDM(JN)
            STFLUX(I,J,JN)=(XFLUX(I,J,JN)/DELT)/WIDM(JN)
            SBLWT(I,J,JN)= BLWTM(I,J,JN)
            SBLWTA(I,J,JN)=((TIMEMIN1*(SBLWTA(I,J,JN)))+BLWTM(I,J,JN))
     1/TIMEPRES
            SCVELI(I,J,JN)=CVELM(I,J,JN)
            SCVELIA(I,J,JN)=((TIMEMIN1*SCVELIA(I,J,JN))+CVELM(I,J,JN))
     1/TIMEPRES
            SPRTHICK(I,J,JN)=SPRTHICK(I,J,JN)+PRTHICK(I,J,JN)
            TOTTHICK(I,J,JN)=TOTTHICK(I,J,JN)+PRTHICK(I,J,JN)
132            CONTINUE
133      CONTINUE
      COUNTER=COUNTER+1
135   CONTINUE
      RETURN
      END
C    ******************************************************************
C    *                                                                *
C    *SUBROUTINE TECTONICS                                            *
C    *                                                                *
C    *THIS SUBROUTINE ADJUSTS BED ELEVATION DUE TO SUBSIDENCE/UPLIFT  *
C    *                                                                *
C    *TECTONICS.F                                                     *
C    *                                                                *
C    *VERS. 1.0 updated 7TH MARCH, 1995                                       *
C    *WRITTEN BY RUTH ROBINSON                                        *
C    ******************************************************************
      SUBROUTINE TECTONICS
C----------------------------------------------------------------------
      IMPLICIT NONE
C----------------------------------------------------------------------
C        DECLARATION STATEMENTS - EXTERNAL VARIABLES - COMMON BLOCKS
C----------------------------------------------------------------------
      INCLUDE 'common.h'
C----------------------------------------------------------------------
C        DECLARATION STATEMENTS - INTERNAL VARIABLES
C----------------------------------------------------------------------
      INTEGER*4 JN
C----------------------------------------------------------------------
      DO JN=1,ZEROLOC-1
         IF(OPTTECT.EQ.1)THEN
C----------------------------------------------------------------------
C     FLEXURAL SUBSIDENCE IS LINEAR, GREATEST AT JN=0
C----------------------------------------------------------------------
	        ELEVM(JN)=ELEVM(JN)-(((ZEROLOC-JN)/(ZEROLOC-2))*TECRATE)
         ELSEIF(OPTTECT.EQ.2)THEN
C----------------------------------------------------------------------
C     FLEXURAL SUBSIDENCE IS AN EXPONENTIAL FUNCTION, GREATEST AT JN=0
C     FLEXURE HAS AN E-FOLDING LENGTH EFOLDX
C----------------------------------------------------------------------
	        ELEVM(JN)=ELEVM(JN)-TECRATE*EXP(-((JN-2)/EFOLDX))
         ELSEIF(OPTTECT.EQ.3)THEN
C----------------------------------------------------------------------
C     FLEXURAL SUBSIDENCE IS SIMILAR TO AN EXP FUNCTION BUT WITH MORE
C     SUBSIDENCE AT JN=0
C----------------------------------------------------------------------
	        ELEVM(JN)=ELEVM(JN)-((1/JN)*NT*TECRATE)
         ENDIF
         ELEVOM(JN)=ELEVM(JN)
      ENDDO
      RETURN
      END
C    ******************************************************************
C    *                                                                *
C    *SUBROUTINE TOPRINT                                              *
C    *                                                                *
C    *THIS SUBROUTINE PRINTS OUT ENTRAINMENT VALUES                   *
C    *                                                                *
C    * TOPRINT.F                                                      *
C    *                                                                *
C    * VERS. 4.0  Update 15 October, 1990                             *
C    ******************************************************************
      SUBROUTINE TOPRINT
C----------------------------------------------------------------------
C               IMPLICIT STATEMENT
C----------------------------------------------------------------------
      IMPLICIT NONE
C----------------------------------------------------------------------
C          DECLARATION STATEMENTS--EXTERNAL VARIABLES--    COMMON BLOCK
C----------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
C----------------------------------------------------------------------
C               DECLARATION STATEMENTS--INTERNAL VARIABLES
C----------------------------------------------------------------------
      REAL*8 POS
      INTEGER*4 J
      INTEGER*4 I
      INTEGER*4 JN,K,DUM
C**********************************************************************
      IF(NT.EQ.0)THEN
C----------------------------------------------------------------------
C     WRITE OUT INITIAL PARAMETERS
C----------------------------------------------------------------------
        WRITE(7,2)
2         FORMAT(' ',7X, 'SEDIMENT TRANSPORT MODEL FOR HETEROGENEOUS SIZ
     #E-DENSITY MIXTURES'/,30X,'VERSION 4.0'/)
        WRITE(7,1)RUN
        WRITE(7,10)
        GOTO(83,86)OPTURB
83      WRITE(7,3)
3       FORMAT(' ','BED SHEAR STRESS IS GAUSSIAN DISTRIBUTED WITH')
        GOTO 87
86      WRITE(7,4)
4       FORMAT(' ','BED SHEAR STRESS IS LOG-GAUSSIAN DISTRIBUTED WITH')
87      WRITE(7,5)COVAR
5       FORMAT(' ',' COEFFICIENT OF VARIATION',25X,F10.2)
        WRITE(7,6)RHO,G,VISC,VISKIN,VOLCON,D50M(1),A,B,THETAC,ACTIVE
     #,DELT
6         FORMAT(' ','FLUID DENSITY',37X,F10.1,'KG M-3'/1X,    'GRAVITAT
     #IONAL ACCELERATION',24X,F10.2,'M S-2'/1X,'MOLECULAR VISCOSITY',31
     #X,F10.5,'PA S'/1X,'KINEMATIC VISCOSITY',31X,F10.7,'M2 S-1'/1X, 'VO
     #LUME CONCENTRATION OF BED SEDIMENT',14X,F10.2/1X, 'MEDIAN GRAIN DI
     #AMETER OF AVAILABLE SEDIMENT',7X,F10.7,'M'/1X, 'A IN BEDLOAD EQUAT
     #ION',29X,F10.2/1X,'B IN SUSPENSION EQUATION',26X,F10.2/1X, 'DIMENS
     #IONLESS BED SHEAR STRESS FOR ROUGH BOUNDARY',1X,F10.3/1X, 'ACTIVE
     #LAYER THICKNESS',28X,F10.3,1X,'M'/1X,'TIME INCREMENT STEP (DELT)',
     #23X,F10.1,1X,'S')
      IF((COVAR.GT.0.35).AND.(OPTURB.EQ.1))THEN
         WRITE(7,7)
7        FORMAT(' ','WARNING!CHOSEN COEFFICIENT OF VARIATION CAUSES NEGA
     #TIVE SHEAR STRESSES'//' IN A NORMAL DISTRIBUTION.A LOGNORMAL DISTR
     #IBUTION SHOULD BE USED')
      ENDIF
      WRITE(7,10)
      ENDIF
C----------------------------------------------------------------------
C     WRITE OUT SHEAR STRESS VALUES AND THEIR DURATION, PLUS ARRAY
C     INDICES FOR EACH MINERAL SPECIES
C----------------------------------------------------------------------
      XPRINT=0
      WRITE(23,*)NT
      DO 151 JN=2,M
       POS=(JN-2)*DELX
       IF(POS.EQ.XPRINT)THEN
          WRITE(7,20)
          WRITE(7,8)(JN-2)*DELX,NT,TOXM(JN),D50M(JN)
          WRITE(7,20)
          WRITE(23,12) (JN-2)*DELX,TOXM(JN),D50M(JN),THET50(JN)
8         FORMAT(' '//10X,'METRES ALONG THE REACH = ',F10.2,' M',
     #/24X,'TIMESTEP = ',I8,/11X,'MEAN BED SHEAR STRESS = ',F10.2,' PA'
     #,/11X,'MEDIAN GRAIN DIAMETER = ',F10.8,1X,'M',/)
          DO 148 J=1,NSIGMA
            WRITE(7,9)J,SIGMA(J)
9           FORMAT(' ',//,27X,'MINERAL SPECIES #',I2,/20X,'MINERAL DENSI
     #TY',F10.1,'KG M-3',//1X,'BED SHEAR STRESS  MINIMUM  ARRAY   MAXIMU
     #M ARRAY',13X,'GRAIN SIZE AV-BLE   CRITICAL BED    CRITICAL    SETT
     #LING  ARRAY'/
     #3X,'MIDPOINT (PA)   SIZE(M)  INDEX   SIZE(M) INDEX  DURATION',3X,
     #'MIDPOINT  WT. PROP. SHEAR STRESS SHEAR VELOCITY VELOCITY NUMBER
     #'/
     #70X,'    (M)                  (PA)        (M S-1)      (M S-1)'/)
          DUM=MIN(NK,NSIZE(J))
          DO 137 K=1,DUM
            I=K
            WRITE(7,11)TO(K,JN),DMIN(K,J,JN),IMIN(K,J,JN),DMAX(K,J,JN),
     1IMAX(K,J,JN),PROPT(K),DIMID(I,J),AWPM(I,J,JN),TC(I,J,JN),
     2SHVELC(I,J,JN),SETV(I,J),I
            WRITE(23,13) DIMID(I,J),TC(I,J,JN)
137       CONTINUE
13    FORMAT(F10.7,4X,F9.3)
          DO 143 K=DUM+1,NK
            WRITE(7,11)TO(K,JN),DMIN(K,J,JN),IMIN(K,J,JN),DMAX(K,J,JN),
     #IMAX(K,J,JN),PROPT(K)
11          FORMAT(F10.2,7X,F10.7,I4,2X,F10.7,I4,2X,F10.7,3X,
     #2F10.6,F9.3,F13.3,6X,F9.3,2X,I3)
12          FORMAT(F10.2,3X,F10.2,3X,F10.8,3X,F10.2)
143       CONTINUE
          DO 147 I=DUM+1,NSIZE(J)
           WRITE(7,11)DIMID(I,J),AWPM(I,J,JN),TC(I,J,JN),
     #SHVELC(I,J,JN),SETV(I,J),I
147       CONTINUE
148     CONTINUE
        WRITE(7,20)
       XPRINT=XPRINT+PXINTERVAL
       ENDIF
151   CONTINUE
      RETURN
1     FORMAT(' ',A72//)
10    FORMAT(1X,'-------------------------------------------------------
     #---------')
20    FORMAT(1X,'-------------------------------------------------------
     #------------------------------------------------------------------
     #---')
      END
C***********************************************************************
C***********************************************************************
C      SUBROUTINE TRIB.F
C***********************************************************************
C      IN ADDING MATERIAL AT EACH NODE TO SIMULATE A TRIBUTARY EFFECT, IT
C      IS ASSUMED THAT ALL TRIBUTARY MATERIAL IS DEPOSITED AND THEREFORE
C      THE HEIGHT OF THE BED IS ADJUSTED WITH A REPRESENTATIVE VOLUME 
C      TERM. THE VOLUMETRIC PROPORTIONS OF SIZES, AND DENSITIES IN THE
C      BED ARE ALSO READJUSTED
C
C      TRIB.F
C
C      VERS. 1.0 WRITTEN 23 JANUARY, 1995 BY RUTH ROBINSON
C***********************************************************************
       SUBROUTINE TRIB
C-----------------------------------------------------------------------
C      IMPLICIT STATEMENT
C-----------------------------------------------------------------------
       IMPLICIT NONE
C-----------------------------------------------------------------------
       INCLUDE 'common.h'
C-----------------------------------------------------------------------
C      INTERNAL DECLARATION STATEMENTS
C-----------------------------------------------------------------------
       REAL*8 TMSDEP(100),TTHICKM,TMSFLUX,TMP
       INTEGER*4 I,J,JN
C----------------------------------------------------------------------
C     CALCULATE THE CHANGE IN AWPM OF BED MATERIAL DUE TO TRIBUTARY
C     SEDIMENT IF OPTTRIB=1. THE CUMULATIVE EFFECT OF ADDING AWPT 
C     AT EACH NODE IS CALCULATED BY SUMMING THE ADDITIONS USING TMP
C     WHICH REPRESENTS THE AWPM FROM THE PREVIOUS NODE. RAR ADDED
C     THIS 1/25/95
C----------------------------------------------------------------------
      IF(NT.EQ.0)THEN
        DO 300 J=1,NSIGMA
        DO 300 I=1,NSIZE(J)
         AWPT(I,J)=AWPT(I,J)*PROPTJ(J)/100.
300     CONTINUE
        DO 200 JN=1,M
          AVDENS(JN)=0.0D0
           DO 205 J=1,NSIGMA
           DO 205 I=1,NSIZE(J)
             AVDENS(JN)=AVDENS(JN)+(AWPM(I,J,JN)*SIGMA(J))
205        CONTINUE
           TMSDEP(JN)=0.0
           TTHICKM=0.0
            DO 210 J=1,NSIGMA
            DO 210 I=1,NSIZE(J)
             BLWTT(I,J)=(BLWTTO(J)/M)*AWPT(I,J)
             TMSFLUX=BLWTT(I,J)*WIDM(JN)*(DELT/DELX)
C-----------------------------------------------------------------------
C       IN SUSP, AND ERODEP, THICKM AND MSPRTHK AND PRTHICK ARE NEGATIVE
C       FOR DEPOSITION; THEY ARE TRANSLATED BELOW FOR THE TRIB. VALUES
C-----------------------------------------------------------------------
             TPRTHICK(I,J,JN)=-(TMSFLUX/(VOLCON*SIGMA(J)*WIDM(JN)))
             TMSPRTHK(I,J,JN)=-(TMSFLUX/WIDM(JN))
             TTHICKM=TTHICKM+TPRTHICK(I,J,JN)
             TMSDEP(JN)=TMSDEP(JN)+TMSPRTHK(I,J,JN)
C-----------------------------------------------------------------------
C     FIND THE RATIO OF MASS IN ACTIVE LAYER TO MASS COMING IN FROM
C     THE SIDES. THIS IS THE WEIGHTING FACTOR THAT IS USED TO REDIST-
C     RIBUTE THE NEW AWPM IN THE ACTIVE LAYER.
C----------------------------------------------------------------------- 
210         CONTINUE
            WEITFAC(JN)=ABS(TMSDEP(JN))/(ACTIVE*VOLCON*AVDENS(JN))
            ELEVM(JN)=ELEVM(JN)-TTHICKM
200       CONTINUE   
          DO 100 JN=1,M
           DO 105 J=1,NSIGMA
           DO 105 I=1,NSIZE(J)
             AWPM(I,J,JN)=AWPM(I,J,JN)*(1-WEITFAC(JN))+(AWPT(I,J)*
     1                   WEITFAC(JN))
             AWPOM(I,J,JN)=AWPM(I,J,JN)
105          CONTINUE
C-----------------------------------------------------------------------
C     ASSIGN WEITFAC THE VALUE OF TMSDEP FOR EACH NODE FOR LATER 
C     CALCULATION
C-------------------------------------------------------------------------
             WEITFAC(JN)=TMSDEP(JN)
100       CONTINUE
          ELEVM(M+1)=2.*ELEVM(M)-ELEVM(M-1)
      ELSE
      DO 90 JN=1,M
        NWAVDENS(JN)=0.0
         DO J=1,NSIGMA
         DO I=1,NSIZE(J)
           NWAVDENS(JN)=NWAVDENS(JN)+(AWPM(I,J,JN)*SIGMA(J))
         END DO
         END DO
         TMP=WEITFAC(JN)
         WEITFAC(JN)=ABS(TMP)/(ACTIVE*VOLCON*NWAVDENS(JN))
          DO 95 J=1,NSIGMA
          DO 95 I=1,NSIZE(J)
           AWPM(I,J,JN)=AWPM(I,J,JN)*(1-WEITFAC(JN))+(AWPT(I,J)*
     1                 WEITFAC(JN))
95        CONTINUE
90    CONTINUE 
      ENDIF
      RETURN
      END
C*****************************SUBROUTINE TURB*************************
C*********************************************************************
C      *CALCULATES THE GAUSSIAN OR LOG-GAUSSIAN DISTRIBUTION OF
C       INSTANTANEOUS SHEAR STRESSES ON THE BED, GIVEN A MEAN AND
C       COEFFICIENT OF VARIATION
C
C      TURB.F
C
C      VERS. 4.0  Update 15 October, 1990
C**********************************************************************
      SUBROUTINE TURB
C----------------------------------------------------------------------
C            IMPLICIT STATEMENT
C----------------------------------------------------------------------
      IMPLICIT NONE
C----------------------------------------------------------------------
C----------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
C----------------------------------------------------------------------
        REAL*8 TOMEAN
C----------------------------------------------------------------------
C            DECLARATION STATEMENTS--INTERNAL VARIABLES
C----------------------------------------------------------------------
      REAL*8 CONV,DISTMIN1,DIST,DELTO,DELTOH,TOLO
      REAL*8 FO,FI,H,T,P,A1,A2,A3,A4,A5,H1,PISQR,STD
      INTEGER*4 JN,K,FLAG
C**********************************************************************
C     DETERMINE NECESSARY COEFFICIENTS                                -
C----------------------------------------------------------------------
      P=.2316419
      A1=.31938153
      A2=-0.356563782
      A3=1.781477937
      A4=-1.821255978
      A5=1.330274427
      CONV=0.693147181
      PISQR=SQRT(2.0*3.1415926536)
C----------------------------------------------------------------------
C    INITIALIZE INCREMENT IN INST BED SHEAR STRESS AND MIDPOINT
C    VALUES OF BED SHEAR STRESS RANGES
C----------------------------------------------------------------------
      DO 96 JN=2,M
C----------------------------------------------------------------------
C       IF COVAR=0.333 THEN TO(K,JN) VARIES FROM 0 TO TOXM(JN)+3SIGMA
C       IF COVAR>0.333 THEN THE CALCULATED TOXM(JN) IS NO LONGER
C       EQUAL TO THE INPUT TOXM(JN) DUE TO TRUNCATION
C----------------------------------------------------------------------
        IF(OPTURB.EQ.1)THEN
          TOMEAN=TOXM(JN)
        ELSE
          TOMEAN=LOG(TOXM(JN))
        ENDIF
        STD=COVAR*TOMEAN
        DELTO=(STD*6.0D0)/FLOAT(NK)
        DELTOH=DELTO*0.5D0
        IF(COVAR.LE.0.333D0) THEN
          TO(1,JN)=TOMEAN-(STD*3.0D0)+DELTOH
        ELSE
          TO(1,JN)=DELTOH
        ENDIF
        TOLO=TO(1,JN)-DELTOH+0.00001
        DO 85 K=2,NK
85        TO(K,JN)=TO(K-1,JN)+DELTO
        IF(OPTURB.EQ.2)THEN
          DO 88 K=1,NK
88           TO(K,JN)=EXP(TO(K,JN))
        ENDIF
C----------------------------------------------------------------------
C    CORRESPONDING VALUE OF SHEAR VELOCITY
C----------------------------------------------------------------------
       DO 95 K=1,NK
        SHVEL(K,JN)=DSQRT(TO(K,JN)/RHO)
95     CONTINUE
96     CONTINUE
C----------------------------------------------------------------------
C     GAUSSIAN DISTRIBUTION OF BED SHEAR STRESS
C----------------------------------------------------------------------
      DISTMIN1=(TOLO-TOMEAN)/STD
      H1=(EXP(-DISTMIN1*DISTMIN1/2.))/PISQR
      DISTMIN1=-DISTMIN1
      T=1/(1+P*DISTMIN1)
      H=((((A5*T+A4)*T+A3)*T+A2)*T+A1)*T
      FO=H*H1
      DO 122 K=1,NK
        FLAG=0
        DIST=((TOLO+K*DELTO)-TOMEAN)/STD
        H1=(EXP(-DIST*DIST/2.))/PISQR
        IF(DIST.LT.0.)THEN 
          DIST=-DIST
          FLAG=1
        ENDIF
        T=1/(1+P*DIST)
        H=((((A5*T+A4)*T+A3)*T+A2)*T+A1)*T
        FI=1-H*H1
        IF(FLAG.EQ.1)THEN
           FI=1-FI
        ENDIF
        PROPT(K)=FI-FO
        FO=FI
122   CONTINUE
      RETURN
      END
C**************************SUBROUTINE WRFLDTA**************************
C**********************************************************************
C      *WRITES INFORMATION ABOUT THE CHANNEL FORM AND PROFILE TO OUTPUT
C       FILE.
C
C      WRFLDTA.F
C
C      VERS. 4.0  Update 15 October, 1990
C**********************************************************************
      SUBROUTINE WRFLDTA
C----------------------------------------------------------------------
C            IMPLICIT STATEMENT
C----------------------------------------------------------------------
      IMPLICIT NONE
C----------------------------------------------------------------------
C          DECLARATION STATEMENTS--EXTERNAL VARIABLES      COMMON BLOCK
C----------------------------------------------------------------------
C----------------------------------------------------------------------
      INCLUDE 'common.h'
C----------------------------------------------------------------------
C----------------------------------------------------------------------
C            DECLARATION STATEMENTS--INTERNAL VARIABLES
C----------------------------------------------------------------------
      REAL*8 POS,VV,SF,RCH,FR
      INTEGER*4 JN,COUNT
C*******************************************************************
C----------------------------------------------------------------------
C      DEFINE PARAMETERS
C----------------------------------------------------------------------
      IF(NT.EQ.0) THEN
            WRITE(8,1)RUN,ICONT,DELX,DELT,PXINTERVAL,PTINTERVAL
1        FORMAT (1X, 'SEDIMENT TRANSPORT MODEL FOR HETEROGENEOUS SIZE',
     1                  ' - DENSITY MIXTURES',/
     1                  '                          VERSION 4.0',//
     1                  ' TITLE : ',A72,//
     1                  '------------------------------------------',/
     1                  ' THE NUMBER OF NODES (ICONT)  ',I4,/
     1                  ' THE DIST BETW. NODES(DELX)   ',F10.3, ' M',/
     1                  ' THE SIZE OF TIMESTEP (DELT)  ',E10.4,'  S',/
     1                  ' PRINT EVERY (PXINTERVAL)      ',I10 ,'  M',/
     1                  ' PRINT EVERY (PTINTERVAL)      ',I10 , ' NT',/
     1                  '------------------------------------------',//)
      ENDIF
      WRITE(8,2)NT
2     FORMAT (1X,' AT TIMESTEP',I9//
     #  ' DISTANCE   BED ELEVATION  WATER DISCHARGE  FLOW DEPTH ',
     #  ' CHANNEL WIDTH  HYD.RAD   MEAN VEL.   FROUDE#  ',
     #  '  FRICTION SLOPE   ',
     #  ' T eff.        D50',/
     #  '      (M)            (M)         (M3/S)          (M)      ',
     #  '     (M)        (M)     (M/S)                  ',
     #  '                  (Pa)         (M)')
      COUNT=0
      XPRINT=0
      DO 97 JN=2,M
C---------------------------------------------------------------------
C            CHECK TO SEE IF ITS IN POSITION TO PRINT
C---------------------------------------------------------------------
         IF((JN-2)*DELX.EQ.XPRINT)THEN
            POS=STRDATA(4,JN)
            RCH=STRDATA(1,JN)
            SF=STRDATA(3,JN)
            VV=STRDATA(5,JN)
            FR=STRDATA(2,JN)
            WRITE(8,4)POS,ELEVM(JN),QN(JN),DCHM(JN),WIDM(JN),RCH, VV,
     #           FR,SF,TOXM(JN),D50M(JN)
C----------------------------------------------------------------------
C               WRITE TO PLOT FILE
C----------------------------------------------------------------------
            COUNT=COUNT+1
            XPRINT=XPRINT+PXINTERVAL
         ENDIF
97    CONTINUE
      WRITE(8,3)
3     FORMAT(' ','----------------------------------',
     #          '-------------------------------------',
     #         '------------------------------',
     #        '------------------------------------------')
4         FORMAT(E10.4,7X,F8.3,4X,F10.3,6X,F8.3,7X,F8.3,3X,F6.3,2X,
     #    E10.3,4X,F6.3,8X,E10.3,3X,E10.3,2X,E10.3)
      RETURN
      END
C*******************************************************************************
C*    SUBROUTINE WRITERSTRT                                                    *
C*                                                                             *
C*    CREATED BY AVN AND KRV TO WRITE TO THE .e FILE ALL THE PARAMETERS        *
C*    NECESSARY TO DO A RESTART,  AT TIME NT=TFILE                             *
C*                                                                             *
C*    WRITERSTRT.F                                                             *
C*                                                                             *
C*    VERS. 4.0  Update 15 October, 1990                                       *
C*    VERS. 4.1  Update June, 1994 by Ruth Robinson
C*******************************************************************************
      SUBROUTINE WRITERSTRT
      IMPLICIT NONE
C-------------------------------------------------------------------------------
C--    DECLARATION STATEMENT -- EXTERNAL VARIABLE  
C-------------------------------------------------------------------------------
       INCLUDE 'common.h'
C----------------------------------------------------------------------
      INTEGER*4 I,J,JN,K,LEN
      CHARACTER*24  F,FF
      CHARACTER*10 FILENAME
C***********************************************************************
      IF(ABS(NT).GE.10.)LEN=2
      IF(ABS(NT).GE.100)LEN=3
      IF(ABS(NT).GE.1000)LEN=4
      IF(ABS(NT).GE.10000)LEN=5
      IF(ABS(NT).GE.100000)LEN=6
      IF(ABS(NT).GE.1000000)LEN=7
      IF(ABS(NT).GE.10000000)LEN=8
      WRITE(UNIT=F,FMT='("(i",i1,",a1)")')LEN
      WRITE(UNIT=F, FMT=F) INT(NT), CHAR(0)
      FF=FOUT(1:KR)//'e'//F
      OPEN(2,FILE='./rst'//FF)
      XPRINT=0
      WRITE(2,1) N,NK,OPTURB,COVAR
      WRITE(2,2) RLONG,DELX,DELT/3600.
      WRITE(2,3) RHO,G,VISC,MANN(2)
      WRITE(2,4) A,B,VOLCON,THETAC,ACTIVE
      WRITE(2,5) OPTHIDE,OPTFUNC,OPTSALT,OPTLNODE
      WRITE(2,6) CON1,CON2,PHI,PSI,BETA
      WRITE(2,7) OPTBINP,OPTSINP,OPTERO,OPTEL,OPTWID,OPTSUSP,
     #OPTABRAD,OPTTECT,OPTDISC,TECRATE,(ZEROLOC-2)*DELX,(EFOLDX-2)*DELX
      WRITE(2,8) NTIM
       IF(OPTDISC.NE.1)THEN
          WRITE(2,9) QDSMOD,QFAC,DELQ,OPTTRIB,QSTMOD,THALFB,RENEW,
     # QSSMODSTRT,QSTMODSTRT
       ELSEIF(OPTDISC.EQ.1)THEN
          WRITE(2,25)FILENAME
       ENDIF
       WRITE(2,10) (PRINT(J),J=1,7)
       WRITE(2,11) XPRINT,PXINTERVAL,(NT+PTINTERVAL),PTINTERVAL,
     # PLOTPRINT,(NT+TFILE),TFILE
      WRITE(2,12) NSIGMA
      DO 119 J=1,NSIGMA
         WRITE(2,13)OPTDIST(J),SIGMA(J),NSIZE(J),KK1(J)
119   CONTINUE
      WRITE(2,14)RUN
      WRITE(2,15)RUNNAME
      WRITE(2,24)QDTMOD,THALFQ,QCOUNT
      WRITE(2,16)(WIDM(JN),JN=1,M+1)
      WRITE (2,16)(ELEVM(JN),JN=1,M+1)
      WRITE (2,16)(ELEVOM(JN),JN=1,M+1)
      WRITE (2,16)(D50M(JN),JN=1,M+1)
      DO 131 JN=1,M
      DO 131 J=1,NSIGMA
         WRITE (2,16)(AWPM(I,J,JN),I=1,NSIZE(J))
         WRITE (2,16)(AWPOM(I,J,JN),I=1,NSIZE(J))
131   CONTINUE
      WRITE (2,18) NT,NTOT,COUNTER
      WRITE (2,19)(DENDM(K),K=1,NTOT)
      WRITE (2,20)DINIT,QO,OPTCONT,ICONT,NEXT
      IF(OPTCONT.EQ.2)THEN
        WRITE(2,21)FLOWEL
      ENDIF
      WRITE (2,22)(DIMIN(J),J=1,NSIGMA)
      DO 144 J=1,NSIGMA
            WRITE (2,16) (AWPT(I,J),I=1,NSIZE(J))
            WRITE (2,16) (DIMID(I,J),I=1,NSIZE(J))
            WRITE (2,16) (DIFEND(I,J),I=1,NSIZE(J))
            WRITE (2,16) (BLWTM(I,J,1),I=1,NSIZE(J))
            WRITE (2,16) (BLWTT(I,J),I=1,NSIZE(J))
            WRITE (2,16) (BEDCON(I,J,1),I=1,NSIZE(J))
144   CONTINUE
      CLOSE(2)
1      FORMAT(3(59X,I2,/),59X,F6.4)
2      FORMAT(2(59X,F10.0,/),59X,F10.6)
3      FORMAT(59X,F10.0,/,59X,F5.3,/,59X,F8.6,/,59X,F8.6)
4      FORMAT(59X,F6.3,4(/,59X,F8.6))
5      FORMAT(3(59X,I1,/),59X,I1)
6      FORMAT(59X,F6.3,4(/,59X,F7.4))
7      FORMAT(9(59X,I1,/),59X,F10.9,/,59X,F10.0,/,59X,F10.0)
8      FORMAT(59X,I10)
9      FORMAT(59X,I1,/,2(59X,F8.4,/),2(59X,I1,/),3(59X,F8.0,/)
     1        ,59X,I8,3(/))
10     FORMAT(/,6(59X,I1,/),59X,I1)
11     FORMAT(3(59X,I8,/,59X,I10,/),59X,I10)
12     FORMAT(59X,I1,5(/))
13     FORMAT(15X,I1,8X,F10.4,5X,I2,5X,F10.9)
14     FORMAT(/,1X,A72)
15     FORMAT(59X,A6,5(/))
16     FORMAT(100E14.6)
18     FORMAT(4(I9))
19     FORMAT(60E14.6)
20     FORMAT(2(F8.3),I2,I5,I12)
21     FORMAT(F8.3)
22     FORMAT(9E14.6)
23     FORMAT(6E14.6)
24     FORMAT(I2,1X,F10.6,1X,I12)
25     FORMAT(59X,A10)
       RETURN
       END
