SUBROUTINE QS_Vrijn04 (hd,u,hs,tp,wavang,
     *            d10,d50,d90,dss,rhos,
     *            te,sa,sus,bed,susw,bedw,stot)


c    Both parameters should be the same as in the calling routine     
C     parameter     (nvert=100,nfrac=20)
      parameter     (nvert=100,nfrac=100)
C
C
      dimension zz(nvert),uca(nvert),ucwa(nvert),caf(nvert),
     *          cat(nvert),ssca(nvert),sswa(nvert)
c
      character line*80
      real bslx,bsly,type, pclay, dt10,dt50,dt90,dtss,sus,bed,susw,bedw
      dimension dsi(nfrac),psi(nfrac),wsi(nfrac),dsteri(nfrac),
     *          dshi(nfrac),egfi(nfrac),taucri(nfrac),tacr1i(nfrac),
     *          cai(nfrac),tcwi(nfrac),betai(nfrac),fc1i(nfrac),
     *          fw1i(nfrac)
      external fwl2004,tolx2004
      logical log
      integer ns

c Calibration coefficient used in original TR2004 code
      FCR=1. 
c      sus =1.0
c      bed =1.0
c      susw=1.0
c      bedw=1.0
      
C cilia
      bslx = 0.       !x-component bed slope
      bsly = 0.       !y-component bed slope
      v = 0.          !v-component current
      type = 3.       !environment (coastal sea)
      ns = 1          !number of fractions
      pclay = 0       !fraction mud
      dt10 = d10*1.e-6
      dt90 = d90*1.e-6
      dt50 = d50*1.e-6
      dtss = dss*1.e-6
c cilia

c======================================================================
c Interface to TRANSPOR1993/2000/2004 Conventions
c======================================================================
      pi=4.*atan(1.)
      deg2rad=pi/180.
      wavang   = wavang*deg2rad
      curang   = atan2(v,u)
      vr       = sqrt (u*u+v*v)
      ur       = 0.
      if (wavang.lt.0.) wavang=2.*pi+wavang
      wavcart  = 1.5*pi - wavang
      if (wavcart.lt.0.) wavcart=2.*pi+wavcart
      ywav=sin(wavcart)
      xwav=cos(wavcart)
      phi      = curang - wavcart
      phi      = 180.*phi/pi
      if (phi.lt.0.) phi=360.+phi
c======================================================================
c End of Interface to TRANSPOR1993/2000/2004 Conventions
c======================================================================
c
c   if ns (fractions ) > 1 calculate d10,d50,d90 etc
c   fraction 0 to 8 micron is clay (Pclay)
c   fraction 8 to 32 micron is fine silt
c   fraction 32 to 62 micron is coarse silt;
c   fraction 62 to 125 micron is fine sand;
c   fraction 125 to 500 micron is medium sand
c   fraction 500 to 2000 micron is coarse sand
c
c
c     fch1= cohesive sediment factor for sediments smaller 
c           than about 62 microns
c     pclay can be estimated 
c     with Pclay=1.-((d50-dclay)/(dsand-dclay))**0.1
c
      dclay=0.000008
      dsilt=0.000032
      dsand=0.000062
      dgravel=0.002
c
c
c     sediments smaller than 0.5dsilt=16 um can be assumed 
c     to be transported as flocs/aggregates??
c
      
c      
      if(abs(vr).lt.0.0001)vr=0.001
c      
      if(ns.gt.1)then
        call calc2004(dsi,psi,dt10,dt50,dt90,dm,egfi,pclay,ns)
      else
        psi(1)=1.
        dsi(1)=dt50
        dm=dt50
        egfi(1)=1.
      endif
c
      fch2=dt50/(1.5*dsand)
      if(fch2.ge.1.)fch2=1.
      if(fch2.lt.0.3)fch2=0.3
      Fsilt=dsand/dt50
      if(Fsilt.le.1.)Fsilt=1.
C
C
C  INITIALIZING PARAMETERS
C


      call init2004L(hd,hs,ub,tp,phi,nn,tanif,g,pi,cl,rho,rhos,rnu,
     *          del,dsh,d50h,wsi,wsb,rkap,dster,dsterm,
     *          dsteri,dt50,dsi,dm,ns,te,dshi,sa,bslx,bsly,dtss)

C
C   WAVE PARAMETERS
C
      IF(HD.LE.0.)STOP 'WATER DEPTH .LE. 0'

	
	IF (TP.GE.1.) THEN
         xx=.1
         yy= 2.*tp*tp
         uugg=VR
         VR=0.
         call zeroin2004(xx,yy,log,fwl2004,tolx2004,tp,vr,phi,g,pi,hd)
             if(.not.log)then
C                stop 'zeroin2004 no wave length'
                                sbc = 0.
                                ssc = 0.
                                sbw = 0
                                ssw = 0
                                goto 330
             else
                yy=4.0*xx
                xx=.2*xx
                VR=uugg
                call zeroin2004(xx,yy,log,fwl2004,tolx2004,tp,vr,phi,
     &                    g,pi,hd)
                  if(.not.log)then
C                               stop 'zeroin2004 no wave length'
                                        sbc = 0.
                                        ssc = 0.
                                        sbw = 0
                                        ssw = 0
                                        goto 330
                  else
                    rls=xx
                  endif
             endif
C
         arg=2.*pi*hd/rls
         tp1=tp/(1.-(VR*tp*cos(phi))/rls)
             IF (ARG.GT.50.) THEN
                ABW=0.
                FW=0.
                FW1=0.
                UBW=0.
             ELSE
                ABW=HS/(2.*SINH(ARG))
                FW=0.
                FW1=0.
                UBW  = 2.*PI/TP1*ABW
             ENDIF
       ELSE
C      If TP less then 1 s; TP is set to 1 s in subroutine init2004.
         ABW=0.
         RLS=0.
         FW=0.
         FW1=0.
         UBW=0.
         tp1=tp
       ENDIF

 

c
c   VELOCITY VECTOR DUE TO MAIN CURRENT AND WAVE-RELATED RETURN CURRENT
c
      htrough=hd*(0.95-0.35*(hs/hd))
      if(abs(ur-9.).lt.1.e-6)
     *ur=-0.125*g**0.5*hs**2./(hd**0.5*htrough)
c
c   ur is defined with respect to the wave dir.; ur is always negative;
c   which means against the wave dir.; if phi=180 (waves against current)
c   then ur cos(phi) will be pos.; thus ur in same direction as vr
c
      Vrr=(VR**2+ur**2+2.*vr*ur*cos(phi))**0.5
C


c   WAVE VELOCITY ASYMMETRY ACCORDING TO ISOBE-HORIKAWA
c    (modified by Grasmeijer April 2003)
c
       rr=-0.4*hs/hd+1.0 
	 
       umax=rr*2.*ubw
       t1=tp1*(g/hd)**0.5
       u1=umax/(g*hd)**0.5
       a11=-0.0049*(t1)**2-0.069*(t1)+0.2911
       raIH=-5.25-6.1*tanh(a11*u1-1.76)
       if(raIH.le.0.5)raIH=0.5
       bs=((bslx)**2.+(bsly)**2.)**0.5 
       rmax=-2.5*hd/rls+0.85
       if(rmax.ge.0.75)rmax=0.75
       if(rmax.le.0.62)rmax=0.62
       ubwfor=umax*(0.5+(rmax-0.5)*tanh((raIH-0.5)/(rmax-0.5)))
       ubwback=umax-ubwfor
  	 
	 
c       ubwfor=
c       ubwback=
c
c       ubwfor=..... can be used if measured
c       ubwback=.....can be used if measured
c
C    REPRESENTATIVE PEAK ORBITAL VELOCITY FOR REFERENCE CONCENTRATION
c
      UBWr=(0.5*ubwfor**3.0+0.5*UBwback**3.0)**(1./3.)

c	write(44,'(12F8.2)') rr, umax,u1,a11,raIH, bs, rmax, 
c     *    ubwfor,ubwback,UBWr, hs, hd
 
C
c    BED ROUGHNESS HEIGHTS (RC and RW)
C


      Uhulp=UBWr+abs(ur)
      Uwc=((Uhulp)**2.+(VR)**2.)**0.5
      RMob=Uwc**2./(del*g*dt50)
c     current-related roughness due to ripples
      RCr=0.
      fcoarse=(0.25*dgravel/dt50)**1.5
      if(fcoarse.ge.1.)fcoarse=1.
      if(RMob.le.50.0)RCr=150*fcoarse*dt50
      if(RMob.gt.50.0.and.RMob.le.250.0)
     *RCr=(-0.65*RMob+182.5)*fcoarse*dt50
      if(RMob.gt.250.0)RCr=20.0*fcoarse*dt50
c      RCr=fcoarse*((85.-65.*tanh(0.015*(RMob-150.)))*dt50)
      if(dt50.le.dsilt)RCr=20.0*dsilt
      if(RCr.le.dt90)RCr=dt90
      if(RCr.ge.0.1)RCr=0.1
C     current-related roughness due to megaripples in rivers, 
c     estuaries and coastal seas and represents the larger
c     scale bed roughness 
      RCmr=0.
      if(RMob.le.50.0)RCmr=0.0002*fch2*Rmob*hd
      if(RMob.gt.50.0.and.RMob.le.550.0)
     *RCmr=fch2*(-0.00002*RMob+0.011)*hd
      if(RMob.gt.550.0)RCmr=0.0
      if(RCmr.gt.0.2)RCmr=0.2
      if(RCmr.lt.0.02.and.dt50.ge.1.5*dsand)RCmr=0.02
      if(RCmr.lt.0.02.and.dt50.lt.1.5*dsand)
     *RCmr=200.*(dt50/(1.5*dsand))*dt50
      if(dt50.lt.1.5*dsand)
     *RCmr=200.*(dt50/(1.5*dsand))*dt50 
      if(dt50.le.dsilt)RCmr=0.
C     current-related roughness due to dunes in rivers
       RCd=0.
      if(TYPE-1.0.le.0.000001)then
      if(RMob.le.100.0.and.Hs.le.0.01)
     *RCd=0.0004*fch2*RMob*hd
      if(RMob.gt.100.0.and.RMob.le.600.0.
     *and.Hs.le.0.01)RCd=fch2*(-0.00008*RMob+0.048)*hd
      if(RMob.gt.600.0.and.Hs.le.0.01)RCd=0.02
      if(dt50.lt.1.5*dsand)
     *RCd=200.*(dt50/(1.5*dsand))*dt50 
      if(dt50.le.dsilt)Rcd=0.
      if(rcd.gt.1.0)RCd=1.
      endif




c
c     wave-related roughness due to ripples
       rwr=0.
       if(RMob.le.50.0)RWr=150*fcoarse*dt50
       if(RMob.gt.50.0.and.RMob.le.250.0)
     *rwr=(-0.65*RMob+182.5)*fcoarse*dt50
       if(RMob.gt.250.0)RWr=20.0*fcoarse*dt50
c       rwr=fcoarse*((85.-65.*tanh(0.015*(RMob-150.)))*dt50)
       if(dt50.le.dsilt)RWr=20.0*dsilt
       if(RWr.ge.0.1)RWr=0.1
c     
       RC=(RCr**2+RCmr**2+RCd**2)**0.5
       RW=RWr
       if(RC.lt.0.001)RC=0.001




c
c
C    STREAMING AT EDGE OF WAVE BOUNDARY LAYER
C
       uspar=((UBWfor+UBWback)/2.)**2./(RLS/TP1)
       RRR=TP1*((UBWfor+UBWback)/2./(2.*pi))/RW
       if(RRR.le.1.)ustream=-1.0*uspar
       if(RRR.gt.1..and.RRR.le.100.)
     *ustream=(-1.+0.875*Alog10(RRR))*uspar
       if(RRR.gt.100.)ustream=0.75*uspar
       ub=ustream
       



c      
c
c     FACTOR FOR ASYMMETRY-RELATED SUSPENDED TRANSPORT
c
      asym=0.
      if(abw.gt.0.)asym=ubwfor/(ubwfor+ubwback)
      UBF3=(UBWFOR)**3.
      UBF4=(UBWFOR)**4.
      UBB3=(UBWBACK)**3.
      UBB4=(UBWBACK)**4.
      FAKTOR=0.
       Phasef=1.
      PP=RWr/(wsb*Tp1)
      PPCR=0.1
      Phasef=-tanh(100.*(PP-PPCR))
      IF(ABS(UBF3+UBB3).GT.0.0001)
     *FAKTOR=0.1*Phasef*(((UBF4-UBB4)/(UBF3+UBB3))+UB)

 
c
c    SUSPENDED SEDIMENT SIZE(DSUS), if ns=1
c     
c     
c      suspended size cannot be smaller than dt50=16 um 
c      disintegration of bed flocs; minimum floc size=16 um
c
       Dsus=(1+0.0006*(dt50/dt10-1.0)*(RMob-550))*dt50
       if(RMob.ge.550.0)Dsus=1.0*dt50
       if(Dsus.le.0.5*(dt10+dt50))Dsus=0.5*(dt10+dt50)
       if(dt50.le.dsilt)Dsus=dt50
       if(dsus.le.0.5*dsilt)dsus=0.5*dsilt
c      lowest value is dsus=16 um
       dshh=0.01*g*del*Dsus**3./rnu/rnu
       if(Dsus.lt.1.5*dsand)wss=(del*g*Dsus*Dsus)/(18.*rnu)
       if(Dsus.ge.1.5*dsand.and.Dsus.lt.0.5*dgravel)
     * wss=(10.*rnu/Dsus)*((1.+dshh)**0.5-1.)
       if(Dsus.ge.0.5*dgravel)wss=1.1*(del*g*Dsus)**0.5
       if(ns.eq.1)wsi(1)=wss

      

c      
c         
C   WAVE-RELATED FRICTION FACTORS
C
C     the "MAX (abrw,1.53056)" makes sure that fw1i(i) can not be larger than 0.3
C     the "MAX (fw1ii,17.9484)" makes sure that fw1i(i) can not be larger than 0.05
      abrw = ABW/RW
      abrw = MAX (abrw,1.53056)
      abdt90 = ABW/dt90
      abdt90 = MAX (abdt90,17.9484)
      IF(ABW.GT.0.)FW=EXP(-6.+5.2*(abrw)**(-0.19))
      IF(ABW.GT.0.)FW1=EXP(-6.+5.2*(abdt90)**(-0.19))
c      IF(FW.GT.0.3 )FW=0.3     !solved with "abrw = MAX (abrw,1.53056)"
c      IF(FW1.GT.0.05)FW1=0.05  !solved with "abdt90 = MAX (abdt90,17.9484)"
C
C   CURRENT-RELATED FRICTION FACTORS
c
      CC=18.*ALOG10(12.*HD/RC)
      CC1=18.*ALOG10(12.*HD/dt90)
      FC=0.24*ALOG10(12.*HD/RC)**(-2.0)
      FC1=0.24*ALOG10(12.*HD/dt90)**(-2.0)
C
C   APPARENT ROUGHNESS      
c     
      GAMMA=0.
      h2=phi
      if(h2.gt.pi)h2=2.*pi-phi
      gamma=0.8+h2-0.3*h2*h2
      uratio=UBWr/vrr
      if(uratio.gt.5.)uratio=5.
      RA=EXP(GAMMA*uratio)*RC
      RCC=10.*RC
      IF(RA.GE.RCC)RA=RCC
      If(RA.GE.0.5*HD)RA=0.5*HD
      CCA=18.*ALOG(12.*HD/RA)
      FCA=0.24*ALOG10(12.*HD/RA)**(-2.0)
      TAUCA=0.125*rho*FCA*vrr*vrr
C
C     CRITICAL SHEAR STRESS ACCORDING TO SHIELDS AND KOMAR-MILLER
C
      call fdster2004(dster,taucr,taucr1,thetcr,pclay,g,dt50,rhos,
     &  rho,fcr)
      if(ns.eq.1)then
       taucrm=taucr
       taucri(1)=taucr
       tacr1i(1)=taucr1
       dsteri(1)=dster
      endif
      umcr=5.75*(del*g*dt50)**0.5*(thetcr)**0.5*log10(4.*hd/dt90)
      ubwcr=0.
      if(tp.gt.0.)then
       if(dt50.le.0.0005)then
             ubwcr=(0.12*del*g*dt50**0.5*tp**0.5)**0.667
       else
             ubwcr=(1.09*del*g*dt50**0.75*tp**0.25)**0.57
       endif
      endif
      if(ns.gt.1)then
      call fdster2004(dsterm,taucrm,tacr1m,dum,pclay,g,dm,rhos,rho,fcr)
c
      do 10 i=1,ns
       call fdster2004(dsteri(i),taucri(i),tacr1i(i),dum,pclay,g,dsi(i),
     *             rhos,rho,fcr)
   10 continue
      endif
C
C  REFERENCE CONCENTRATION CA
C
      aa1=0.5*RCr
      aa2=0.5*RWr
      a=max(aa1,aa2)
      if(a.le.0.01)a=0.01
      if (a.le.rc/30.)a=rc/30.+max(aa1,aa2)
      UST=G**0.5*ABS(vrr)/CC
c
c   wave boundary layer and mixing layer thickness     
c
          DELm=0.
          DELw=0.
          IF(ABW.GT.0.)THEN
             DELW=0.36*ABW*(ABW/RW)**(-0.25)
             DELm=2.*DELw
          ENDIF
          DELm=max(DELm,rc)
          IF(DELm.le.Ra/29.9)DELm=ra/29.9
          If(DELm.ge.0.2)DELm=0.2
          If(DELM.le.0.05)DELm=0.05
c       
      h3=(ALOG(30.*DELM/RA)/ALOG(30.*DELM/RC))**2
      h4=((-1.+alog(30.*hd/rc))/(-1.+alog(30.*hd/ra)))**2 
      alfaw=h3*h4
      if(alfaw.lt.0.)alfaw=0.
      if(alfaw.gt.1.)alfaw=1.
      RMUC=FC1/FC
      if(rmuc.gt.1.)rmuc=1.
      RMUW=0.7/Dster
         if(Dster.ge.5.0)RMUW=0.14
         if(Dster.le.2.0)RMUW=0.35
      TAUC=0.125*rho*FC*vrr*vrr
      TAUW=0.25*rho*FW*UBWr*UBWr
      tauce=rmuc*tauc
      TAUCEF=RMUC*ALFAW*TAUC
      Tc=(tauce-taucr1)/taucr
      Tc=max(0.0001,Tc)
      TAUWEF=RMUW*TAUW
c      if(TAUWEF.le.0.25*rho*FW1*UBWR**2.)TAUWEF=0.25*rho*FW1*UBWR**2.
       TAUCWE=TAUCEF+TAUWEF
      RRR1=0.8+0.2*((taucwe/taucr1-0.8)/1.2)
      if(RRR1.le.0.8)RRR1=0.8
      if(RRR1.ge.1.)RRR1=1.
      Tcw=(TAUCWE-RRR1*TAUCR1)/TAUCR
      Tcw=MAX(.0001,Tcw)
      Thetacw=taucwe/((rhos-rho)*g*dt50)
      CA=0.015*Fsilt*(1.-pclay)*dt50/A*DSTER**(-.3)*Tcw**1.5
      ca=min(ca,0.05)
c
      do 20 i=1,ns
c
c     LOOP OVER FRACTIONS
c     (sum of sand fractions and pclay is equal to 1)
c
c      ks,grain=di
      fc1i(i)=0.24*alog10(12.*hd/dsi(i))**(-2.)
      fw1ii=abw/dsi(i)
C     the "MAX (fw1ii,17.9484)" makes sure that fw1i(i) can not be larger than 0.05
      fw1ii= MAX (fw1ii,17.9484)
      if(fw1ii.gt.0.)fw1ii=fw1ii**(-0.19)
      fw1i(i)=exp(-6.+5.2*fw1ii)
c      IF(fw1i(i).GT.0.05)fw1i(i)=0.05  !solved with "fw1ii= MAX (fw1ii,17.9484)"
c
      RMUWi=0.7/dsteri(i)
         if(dsteri(i).ge.5.0)RMUWi=0.14
         if(dsteri(i).le.2.0)RMUWi=0.35
c         
c      RMUWi is efficiency factor acting on fractions
c
       taucee=0.125*alfaw*rho*fc1i(i)*vrr*vrr
c      taucee=0.125*alfaw*rho*fc1*vrr*vrr
       tauwee=0.25*RMUWi*fw*UBWr*UBWr*rho
       tacwee=taucee+tauwee
c
c     dt50 method for taucr
c
      fsilti=dsand/dsi(i)
      if(fsilti.le.1.)fsilti=1.
      Rlamda=(dsi(i)/dt50)**0.25
      Tcwi(i)=Rlamda*((TACWEE-RRR1*TAUCR1*(dsi(i)/dt50)*
     *(egfi(i)))/(TAUCR*(dsi(i)/dt50)))
cB      Tcwi(i)=((TACWEE-TAUCR1*(dsi(i)/dt50)*
cB     *(egfi(i)))/(TAUCR*(dsi(i)/dt50)))
cC      Tcwi(i)=(TACWEE-TAUCRi(i)*egfi(i))/TAUCRi(i)      
cD      Tcwi(i)=(TACWEE-TAUCRi(i))/TAUCRi(i) 
      Tcwi(i)=MAX(.0001,Tcwi(i))
      Cai(i)=0.015*fsilti*(1.-pclay)*Dsi(i)/A*DSTERi(i)**(-.3)*
     *Tcwi(i)**1.5 
      Cai(i)=min(cai(i),0.05)
c
c     END LOOP OVER FRACTIONS
c
   20 continue
c
C     VELOCITY AND CONCENTRATION PROFILES
c
      sbx=0.
      sby=0.
      ta1som=0.
      tau1x=0.
      tau1y=0.
c
c   grid points in z direction
c
        do 206 ij=1,nn
          zz (ij)=a*(hd/a)**(float(ij-1)/float(nn-1))
  206   continue
c
c  velocities
c
      udel=0.
      udelw=0.
      uc=0.
      ucw=0.
c
c   Velocity profile in current direction UC: Logarithmic in two layers;
C
c   velocity profile in wave direction UCW;
c   Logarithmic in lower half of depth and power profile in upper half
c   (velocity curving back to zero at surface)
c
        alf1=0.5
        alf2=1.
        zaa=ra/30.
        c10=-1.+alog(hd/zaa)
        c20=alog(alf1*hd/zaa)
        c30=-alf1+alf1*alog(alf1*hd/zaa)
        if(abs(alf1-1.).lt.1.e-6.and.abs(alf2-1.).lt.1.e-6)fac1=1.
        if(alf1.lt.1.)fac1=c10/(c30+0.375*c20)
c    ur1 is velocity at mid depth in wave direction
        ur1=c20/c10* ur*fac1
      do 3005 i=1,nn
        z=zz(i)
        HULP30=-1.+ALOG(30.*HD/RA)
        UDELc=VR*ALOG(30.*DELM/RA)/HULP30
        UDELw=uR*ALOG(30.*DELM/RA)/HULP30
        IF(DELM.GT.0..and.z.lt.delm)THEN
            UC=UDELc*ALOG(30.*z/RC)/ALOG(30.*DELM/RC)
            ucw=fac1*UDELw*ALOG(30.*z/RC)/ALOG(30.*DELM/RC)
        else
           IF(z.GE.DELM)UC=VR*ALOG(30.*z/RA)/HULP30
           IF(z.GE.DELM)UCw=fac1*uR*ALOG(30.*z/RA)/HULP30
        endif
        if(z.le.rc/30.)then
          uc=0.
          ucw=0.
        endif
        uca(i)=uc
        ucwa(i)=ucw
        if(z.gt.alf1*hd)then
          ucw=ur1*(1.-((z-alf1*hd)/((alf2-alf1)*hd))**3)
        endif
        if(z.gt.alf2*hd)ucw=0.
        ucwa(i)=ucw
 3005 continue
c
c
      do 299 ii=1,nvert
          caf (ii)=0.
          cat (ii)=0.
          ssca(ii)=0.
          sswa(ii)=0.
  299 continue
c
c    LOOP OVER FRACTIONS
c
      do 305 is=1 ,ns
c
c    MIXING COEFFICIENTS
c
        betai(is)=1.+2.*(wsi(is)/ust)**2.
        if(betai(is).ge.1.5)betai(is)=1.5
        gambr=1.
        if(hs/hd.gt.0.4)gambr=1.+((HS/hd)-0.4)**0.5
         Ds=2.*gambr*DELW
         if(DS.le.0.05)DS=0.05
         if(DS.ge.0.2)DS=0.2
        ustw=(tauw/rho)**0.5
        Betaw=1.+2.*(wsi(is)/ustw)**2
        if(Betaw.ge.1.5)Betaw=1.5
        Betaa=1.+2.*(wsb/ust)**2
        if(betaa.ge.1.5)betaa=1.5
        Fdamp=(250/Rmob)**0.5
        if(Fdamp.ge.1.)Fdamp=1.
        if(Fdamp.le.0.1)Fdamp=0.1
        EBW=0.018*Fdamp*Betaw*gambr*DS*UBWr
        IF (TP.GE.1.) THEN
          EMAXW=0.035*gambr*HD*HS/TP 
        ELSE
          EMAXW=0.
        ENDIF
        IF(EMAXW.LE.EBW)EMAXW=EBW
        IF(EMAXW.GE.0.05)EMAXW=0.05
        EMAXC=0.25*RKAP*UST*HD*BETAi(is)
       JTAL= 4
       if(ns.gt.1)ca=cai(is)
       dym=ca/(nn*jtal)
       DXM = HD/(NN*jtal)
       DYX = DYM/DXM
C
C  COMPUTATION CONCENTRATION PROFILE DC/DY OR DC/DX
C
      C=CA
      Z=A
      ws=wsi(is)
      caf(1)=ca*rhos*psi(is)
      cat(1)=caf(1)+cat(1)
      itt=2
C
C  INTEGRATION FROM Z=A TO SURFACE
C
      Y=CA
      XEND=A
      z=a
      IT=2
      c=ca
C
  100 CONTINUE
      XOLD = XEND
      YOLD = Y
      IF(Z.LE.DS)ESW=EBW
      IF(Z.GT.DS.AND.Z.LE.0.5*HD)ESW=EBW+(EMAXW-EBW)*((Z-DS)/
     *(0.5*HD-DS))
      IF(Z.GE.0.5*HD)ESW=EMAXW
      IF(Z.GE.0.5*HD)ESC=EMAXC
      IF(Z.LT.0.5*HD)ESC=EMAXC-EMAXC*(1.-2.0*Z/HD)**2.
      ES=ESW**2+ESC**2
      if(es.gt.0.)es=sqrt(es)
      fcc=0.
c
c    Derivative dc/dz
c
      IF (C.GT.1.E-9) THEN
        IF(Z.GE.A)then
c        
c         fi=dampings term; fch2=extra damping due to fine silt 
c         and clay<32 microns
c         
c     cmaxs=maximum bed concentration in case of sandy bottom (=0.65)
      cmaxs=0.65
      cmax=(dt50/dsand)*cmaxs
      if(cmax.le.0.05)cmax=0.05
      if(cmax.ge.cmaxs)cmax=cmaxs
      fi=fch2*(1.0+((c/cmaxs)**0.8)-(2.0*(c/cmaxs)**0.4))
      if(fi.lt.0.01)fi=0.01
c
c   fhs=hindered settling term
c
      fhs=(1.-0.65*c/cmax)**5.0
      if(fhs.le.0.01)fhs=0.01
c
c   ffloc=flocculation term; if d<dsilt
c   salmax is limit of salinity; flocculation is fully active if sa>salmax
c
      Ffloc0=1.
      Efloc=0.
c 
      if(dsi(is).lt.dsand)Efloc=(dsand/dsi(is))-1.
      if(Efloc.ge.3.)Efloc=3.
      Fhulp=4.0+alog10(2.*c/cmax)
      if(Fhulp.le.1.)Fhulp=1.
      if(dsi(is).lt.dsand)Ffloc0=Fhulp**Efloc
      if(ffloc0.le.1.)Ffloc0=1.
      if(ffloc0.ge.10.)Ffloc0=10.
      salmax=1.
      if(sa.lt.salmax)Ffloc=(ffloc0-1)*sa/salmax+1.
      if(sa.ge.salmax)Ffloc=ffloc0
      if(sa.le.0.)Ffloc=1.
c      
      fcc=-WS*C*fhs*ffloc/(ES*fi)
        endif
      ENDIF
      YPRIME=fcc
      FF = 1./Y*YPRIME
      IF (-YPRIME .GT. DYX) THEN
c
c  Integration along c-axis
c
         Y = YOLD-DYM
         IF (Y .LT. 2./3.*YOLD) Y = 2./3.*YOLD
         XEND = XOLD+ALOG(Y/YOLD)/FF
      ELSE
c
c  Integration vertical z-axis
c
         XEND =min( XOLD+DXM,zz(itt)+1.e-6)
         Y = EXP(ALOG(YOLD)+(XEND-XOLD)*FF)
      ENDIF
      C=Y
      Z=XEND
c
c     If Tcw<0 then Tcw=0.0001 yielding roughly ca=0.2 10-6 kg/m3
c     climit=1 10-9 equivalent with 2.65 10-6 kg/m3; if c<climit 
c     then fcc=0. Computed c will be kept constant from that point upward;
c     c is little bit smaller than climit (inherent to numerical integration
c      procedure )
c
C
C  RESULTS IN ARRAYS
C
      ittt=itt
      do 207 it=itt,nn
      if( z.ge.zz(it ))then
         caf (IT)  =Y*RHOS*psi(is)
         cat (IT)  =cat (it)  +caf (it)
         ittt=it+1
         goto 1001
      ENDIF
  207 continue
 1001 continue
      itt=ittt
      if(itt.le.nn)GOTO 100
c
c   DEPTH INTEGRATION OF SUSPENDED TRANSPORT
c
      SSC=0.
      sscis=0.
      ssw=0.
      sswis=0.
      ssasym=0.
      ssasis=0.
      uz1=ucwa(1)
      do 201 i=2,nn
           dz=zz(i)-zz(i-1)
           uc=uca(i)
           uz2=ucwa(i)
c    Depth-integration of current-related transport in current dir. (ssc)
           ssc=ssc+(cat(i)*uc+cat(i-1)*uca(i-1))/2.*dz
           ssca(i)=ssc
           sscis=sscis+(caf(i)*uc+caf(i-1)*uca(i-1))/2.*dz
           z=zz(i)
c    Depth-integration of oscillating susp. transport in wave dir.(ssasym)
           if(z.le.3.*DS)then
              ssasym=ssasym+(cat(i)+cat(i-1))/2.*dz
              ssasis=ssasis+(caf(i)+caf(i-1))/2.*dz
           endif
C    Depth-integration of current-related transport in wave direction (ssw)
           torm=(uz1*cat(i-1)+uz2*cat(i))/2.*dz
           turm=(uz1*caf(i-1)+uz2*caf(i))/2.*dz
           ssw=ssw+torm
           sswa(i)=ssw
           sswis=sswis+turm
           uz1=uz2
  201 continue
c
c   SUMMATION OF CURRENT-RELATED AND WAVE-RELATED 
c   TRANSPORT IN WAVE DIRECTION
c
      sswasym=faktor*ssasym
      ssw=ssw+faktor*ssasym
      sswis=sswis+faktor*ssasis
c
c   ssw is pos. in wave direction and neg. against wave dir.
c   ssasym is always pos. (in wave dir.)
c
c   BED LOAD TRANSPORT (instantaneous and averaging)
c
c      vrdelm=vr*alog(30.*delm/ra)/(-1.+alog(30.*hd/ra))
c      urdelm=ur*alog(30.*delm/ra)/(-1.+alog(30.*hd/ra))*fac1
       VRa=UDELc*ALOG(30.*a/RC)/ALOG(30.*DELM/RC)
       URa=fac1*UDELw*ALOG(30.*a/RC)/ALOG(30.*DELM/RC)      
      IF(DELM.LE.RA/30.)THEN
        VRa=0.
        URa=0.
      endif
c      vrdlmr=(vrdelm**2+urdelm**2+2.*vrdelm*urdelm*cos(phi))**0.5
       vrar=(VRa**2+URa**2+2.*VRa*URa*cos(phi))**0.5
c  vrar is absolute value of VRa and URa
      ubtot=ub+URa
      rmfor=(ubwfor+ubtot)**2/(del*g*dt50)
      rmback=(-ubwback+ubtot)**2/(del*g*dt50)
      ntime=201
      dtt=tp1/(ntime)
      ubwf=ubwfor
      ubwb=ubwback
      tfor=0.
      if(abw.gt.0.)tfor=ubwback/(ubwfor+ubwback)*tp1
      tback=tp1-tfor
      profact=((-1.+alog(30.*hd/rc))/alog(30.*a/rc))**2.
      acw1=abs(vrr)/(abs(ubw)+abs(vrr))
C      acw1=abs(vrar)/(abs(ubw)+abs(vrar))
      acw=acw1**1
      sbxis=0.
      sbyis=0.
      ta1xis=0.
      ta1yis=0.
      nnn=0
c      
c      ks,grain=di 
      fc11=0.25*fc1i(is)*profact
      fcw1=acw*fc11+(1.-acw)*fw1i(is)
c
c      fc11=0.25*fc1*profact
c      fcw1=acw*fc11+(1.-acw)*fw1
c
c     fc11 can also be represented as 2(Kappa/ln(30(a/dt90)))**2
c     which can bed derived by assuming a log profile
c
  150 continue
       if(ns.eq.1)then
         fc11=0.25*fc1*profact
         fcw1=acw*fc11+(1.-acw)*fw1
       endif
      time=nnn*dtt
      udt=0.
      Pangle=0.
      plead1=cos(Pangle/180.*pi)
      plead2=sin(Pangle/180.*pi)
c According to Nielsen the phase lead angle is: Pangle=40 degrees
c Phase lead angle (Pangle)=0 yields no effect
c
      if(abw.gt.0.)then
         if(time.lt.tfor)then
         udt=ubwfor*sin(pi*time/tfor)
         udt2=ubwfor*sin(pi*(time+dtt)/tfor)
         udt1=ubwfor*sin(pi*(time-dtt)/tfor)
         dudt=Tp1/(2.*pi)*(udt2-udt1)/(2.*dtt)
         udt=plead1*udt+dudt*plead2
      else
         udt=-ubwback*sin(pi*(time-tfor)/tback)
         udt2=-ubwback*sin(pi*(time+dtt-tfor)/tback)
         udt1=-ubwback*sin(pi*(time-dtt-tfor)/tback)
         dudt=Tp1/(2.*pi)*(udt2-udt1)/(2.*dtt)
         udt=plead1*udt+dudt*plead2
      endif
      endif
      udtx=udt*cos(phi)
      udty=udt*sin(phi)
      ubtotx=ubtot*cos(phi)
      ubtoty=ubtot*sin(phi)
c  ubtot= current velocity near bed in wave direction;
c  ubtot is pos. in wave dir. and neg. against wave dir.
c  x-direction in current direction
c  y-direction normal to current direction
c  wave direction under angle phi to current direction
      uxt=vra+udtx+ubtotx
      uyt=udty+ubtoty
      utvec=(uxt**2 +uyt**2 )
      if(utvec.gt.0.)utvec=sqrt(utvec)
      BSLS=(bslx*UXT/abs(UTVEC))+(bsly*UYT/abs(UTVEC))
      BSLN=(bslx*UYT/abs(UTVEC))+(bsly*UXT/abs(UTVEC))
      fsl1=1./(1.-bsls/tanif)
C     fsl1=Bagnold factor; can be better applied after time-averaging of bed load (see below fsl3)
C     BSLS and BSLN are positive for downsloping beds
C     fsl1=1./(1.+bsls/tanif)if BSLS is defined to be negative for downsloping bed
      angleB=atan(BSLS)
      angleG=atan(BSLN)
      angleP=atan(tanif) 
      fsl2=(1-angleB/angleP)**0.75*(1-angleG/angleP)**0.37
c      fsl2=sin(atan(tanif)+atan(bsls))/sin(atan(tanif))
c     BSLS is the tangens of the slope in direction of velocity vecor (s-direction)
      ta1som=ta1som+0.5*rho*fcw1*utvec*utvec
      tau1t=0.5*rho*fcw1*utvec*utvec
c  dt50 or di method for taucr
      taucr2=taucr1*fsl2
      RRR2=0.8+0.2*((tau1t/taucr1-0.8)/1.2)
      if(RRR2.le.0.8)RRR2=0.8
      if(RRR2.ge.1.)RRR2=1.
      Rlamda=(dsi(is)/dt50)**0.25
      TT=Rlamda*((tau1t-RRR2*taucr2*(dsi(is)/dt50)*egfi(is))/
     *(taucr*(dsi(is)/dt50)))
cB      TT=(tau1t-taucr2*(dsi(is)/dt50)*egfi(is))/
cB     *(taucr*(dsi(is)/dt50))
cC      TT=(tau1t-taucri(is)*egfi(is))/taucri(is)
cD      TT=(tau1t-taucri(is))/taucri(is)
      TT=max(0.0001,TT)
      uster1t=(tau1t/rho)**0.5
      fsilti=dsand/dsi(is)
      if(fsilti.le.1.)fsilti=1.
      sbt=0.5*Fsilti*(1.-pclay)*fsl1*dsi(is)*rhos*uster1t*TT**1.0/
     *(dsteri(is)**0.3)
c
c   bed l. tr. for rivers
c   sh=((rhos-rho)/rho*g)**0.5
c   if(TT.le.3.)then
c   sbt=0.053*sh*rhos*(dsi(is))**1.5*TT**2.1/(dsteri(is)**0.3)
c   else
c   sbt=0.1*sh*rhos*(dsi(is))**1.5*TT**1.5/(dsteri(is)**0.3)
c  endif
c
      sbxis=sbxis+uxt/utvec*sbt*psi(is)
      sbyis=sbyis+uyt/utvec*sbt*psi(is)
      ta1xis=ta1xis+uxt/utvec*tau1t*psi(is)
      ta1yis=ta1yis+uyt/utvec*tau1t*psi(is)
      nnn=nnn+1
      if(nnn.le.ntime)goto 150
      ta1xis=ta1xis/(ntime-1)
      ta1yis=ta1yis/(ntime-1)
      tau1x=tau1x+ta1xis
      tau1y=tau1y+ta1yis
      ta1net=(tau1x**2+tau1y**2)
      if(ta1net.gt.0.)ta1net=sqrt(ta1net)
      phi1t=(180./pi)*atan2(tau1y,tau1x)
      sbxis=sbxis/(ntime-1)
      sbx=sbx+sbxis
      sbyis=sbyis/(ntime-1)
      sby=sby+sbyis
c
c   SUSPENDED LOAD TRANSPORT VECTORS
c
      ssx=ssw*cos(phi)+ssc
      ssxis=sswis*cos(phi)+sscis
      ssy=ssw*sin(phi)
      ssyis=sswis*sin(phi)
c     ssx is pos. in current direction (which is pos. x-direction)
      ssvec=sqrt(ssx**2 +ssy**2 )
      ssvis=sqrt(ssxis**2 +ssyis**2 )
      phis1=(180./pi)*atan2(ssy,ssx)
      phis2=(180./pi)*phi-phis1
c
c   BED LOADTRANSPORT VECTORS
c
      sbvec=sqrt(sbx**2 +sby**2 )
      sbvis=sqrt(sbxis**2 +sbyis**2 )
c      fsl3=1.
c      if(abs(sbvec).gt.0.)then
c      tanb=bslx*sbx/abs(sbx)+bsly*sby/abs(sby)
c      fsl3=1./(1.+1.*tanb/tanif)**1.
c      if(fsl3.gt.2.)fsl3=2.
c      if(fsl3.lt.0.5)fsl3=0.5
c      sbvec=sbvec*fsl3
c      sbvis=sbvis*fsl3
c      endif
      sbc=sbx-sby/tan(phi)
      sbw=sby/sin(phi)
      phib1=0.
      if(.not.(abs(sby).lt.1.e-6.and.abs(sbx).lt.1.e-6))
     *phib1=(180./pi)*atan2(sby,sbx)
      phib2=(180./pi)*phi-phib1
c
c   TOTAL LOAD TRANSPORT VECTORS
c
      stx=ssx+sbx
      stxis=ssxis+sbxis
      sty=ssy+sby
      styis=ssyis+sbyis
      stvec=sqrt(stx**2 +sty**2 )
      stvis=sqrt(stxis**2 +styis**2 )
      phit1=(180./pi)*atan2(sty,stx)
      phit2=(180./pi)*phi-phit1
      stc=stx-sty/tan(phi)
      stw=sty/sin(phi)
c
c
c   BED FORMS
c
      fr=vr/sqrt(g*hd)
      if(fr.gt.1.)fr=1.
c
c
c     dd1,rrl1; dd5,rrl5; dd9,rrl9   =ripples
c     dd2,rrl2; dd6,rrl6; dd10,rrl10 =mega-ripples
c     dd3,rrl3; dd7,rrl7; dd11,rrl11 =dunes
c     dd4,rrl4; dd8,rrl8; dd12,rrl12 =sand waves
c      
c   Bed forms in rivers (small waves are assumed to have no effect)
c
c     ripples
      if(dster.lt.10.0.and.Tcw.lt.25.0)then
        dd1=100.*dt50
        rrl1=750.*dt50
      else
      dd1=0.
      rrl1=0.
      endif
c     mega-ripples
      if(dster.lt.10.0.and.Tcw.ge.1.0
     *.and.Tcw.lt.10.0)then
        dd2=0.02*hd*(10.-Tcw)*(1.-exp(-0.1*Tcw))
        rrl2=0.5*hd
      else
      dd2=0.
      rrl2=0.
      endif
c     dunes
      if(Tcw.ge.1.0.and.Tcw.lt.25.0)then
        dd3=0.11*hd*(dt50/hd)**0.3*(25.-Tcw)*(1.-exp(-0.5*Tcw))
        rrl3=7.3*hd
      else
      dd3=0.
      rrl3=0.
      endif
c     sand waves
      if(Tcw.ge.15.0.and.Fr.lt.0.8.and.
     *hd.ge.1.0)then
        dd4=0.15*hd*(1.-fr**2.)*(1.-exp(-0.5*(Tcw-15.)))
        rrl4=10.*hd
      else
      dd4=0.
      rrl4=0.
      endif
c
c   Bed forms in estuaries (small waves are assumed to have no effect)
c
c     ripples
      if(dster.lt.10.0.and.Tcw.lt.25.0)then
        dd5=100.*dt50
        rrl5=750.*dt50
      else
      dd5=0.
      rrl5=0.
      endif
c     mega-ripples
      if(dster.lt.10..and.Tcw.ge.1.0
     *.and.Tcw.lt.10.0)then
        dd6=0.02*hd*(10.-Tcw)*(1.-exp(-0.1*Tcw))
        rrl6=0.5*hd 
      else
      dd6=0.
      rrl6=0.
      endif 
c     dunes
      if(Tcw.ge.1.0.and.Tcw.lt.25.0.and.
     *hd.ge.1.0)then 
        dd7=0.05*hd
        rrl7=2.*hd
        else
        dd7=0.
        rrl7=0.
      endif
c    sand waves
      if(Tcw.ge.15.0.and.FR.lt.0.8.and.
     *hd.ge.1.0)then
        dd8=0.05*hd
        rrl8=3.*hd
        else
        dd8=0.
        rrl8=0.
      endif
C
C   Bed forms in coastal seas
C
C    ripples
      if(hs.ge.0.01.and.rmob.lt.10.0)then
        dd9=0.22*abw
        rrl9=5.6*dd9
      endif
      if(hs.ge.0.01.and.rmob.gt.10.0.and.rmob.lt.250.0)then
        dd9=2.8e-13*(250.-rmob)**5.*abw
        rrl9=(1./(2.e-7*(250.-rmob)**2.5))*dd9
      endif
      if(hs.ge.0.01.and.rmob.ge.250.0)then
        dd9=0.
        rrl9=0.
      endif
c    mega-ripples
      if(hs.ge.0.01.and.rmob.lt.550.0.and.vr.ge.0.3.and.
     *hd.ge.1.0)then
        dd10=0.02*hd
        rrl10=0.5*hd
        else
        dd10=0.
        rrl10=0.
      endif
c    dunes (dd11;rrl11; no dunes in coastal seas)
      if(hs.ge.0.01.and.rmob.gt.0..and.vr.ge.0.3.and.
     *hd.ge.1.0)then
      dd11=0.
      rrl11=0.
      endif
c    sand waves
      if(hs.gt.0.01.and.vr.ge.0.5.and.hd.ge.10.0)then
        dd12=0.1*hd
        rrl12=10.*hd
        else
        dd12=0.
        rrl12=0.
      endif
C
  305 continue
  330 continue
c END LOOP OVER FRACTIONS
C
c======================================================================
c Interface to Cartesian Conventions
c======================================================================
      sbcx=u*sbc/vr*bed
      sbwx=xwav*sbw*bedw
      sscx=u*ssc/vr*sus
      sswx=xwav*ssw*susw
      sbx = sbcx + sbwx
      ssx = sscx + sswx
      stot = sbx + ssx

              
c======================================================================
c End of Interface 
c======================================================================
C      STOP 'NORMAL END TRANSPOR'
      END

      subroutine init2004L(hd,hs,ur,tp,phi,nn,tanif,g,pi,cl,rho,rhos,
     *                rnu,del,dsh,d50h,wsi,wsb,rkap,dster,dsterm,
     *                dsteri,dt50,dsi,dm,ns,te,dshi,sa,bslx,bsly,dtss)
      real*4 dtss
      dimension wsi(*),dsi(*),dsteri(*),dshi(*)
      if(abs(hs).lt.1.e-6)then
        ur=0.
        ub=0.
      endif
      if(abs(hs).lt.1.e-6.and.abs(tp).lt.1.e-6)tp=7.
      if(tp.ge.0..and.tp.lt.1.)tp=1.
      if(Hs.eq.0..and.phi.eq.0.)phi=90.
      if(abs(hs).lt.1.e-6)hs=hd/5000.
      hs=max(hs,.001)
      if(abs(phi).lt.1.e-6)phi=0.000001
      if(abs(phi-180.).lt.1.e-6)phi=179.99999
      if(abs(phi-360.).lt.1.e-6)phi=359.99999
      NN     = 50
      G      = 9.81
      PI     = 4.*ATAN(1.)
      phi    = phi/180.*pi
c
      dclay=0.000008
      dsilt=0.000032
      dsand=0.000062
      dgravel=0.002
C
C     bslx= tangens of bed slope in current direction 
C      (slope is pos. if velocity is upsloping, bed level to datum 
C       increases with positive x-coordinate; 
C       slope is negative if velocity is downsloping))
C
C     bsly= tangens of bed slope normal to current direction
C  
      bslx=0.
      bsly=0.
      TANIF=0.6
c     bslx and bsly are tangens values of bed slopes (positive for downsloping)
c     TANIF=tangens of angle of repose
      cl=(sa-0.03)/1.805
      rho=1000.+1.455*cl-0.0065*(te-4.+0.4*cl)**2.
c      rhos=2650.
      rnu=(4.e-5)/(20.+te)
      DEL=(RHOS-rho)/rho
      d50h=0.01*g*del*dt50**3./rnu/rnu
      if(dt50.le.1.5*dsand)wsb=(del*g*dt50*dt50)/(18.*rnu)
      if(dt50.gt.1.5*dsand.and.dt50.lt.0.5*dgravel)
     * wsb=(10.*rnu/dt50)*((1.+d50h)**0.5-1.)
      if(dt50.ge.0.5*dgravel)wsb=1.1*(del*g*dt50)**0.5
      do 10 i=1,ns
         if(ns.gt.1)dtss=dsi(i)
         if(dtss.le.0.5*dsilt)dtss=0.5*dsilt
         dsh=0.01*g*del*dtss**3./rnu/rnu
         if(dtss.lt.1.5*dsand)ws=(del*g*dtss*dtss)/(18.*rnu)
         if(dtss.ge.1.5*dsand.and.dtss.lt.0.5*dgravel)
     *   ws=(10.*rnu/dtss)*((1.+dsh)**0.5-1.)
         if(dtss.ge.0.5*dgravel)ws=1.1*(del*g*dtss)**0.5
         DSTERi(i)=dtss*(DEL*G/RNU**2)**(1./3.)
         wsi(i)=ws
         dshi(i)=dsh
   10 continue
      RKAP=.4
      DSTER=dt50*(DEL*G/RNU**2)**(1./3.)
      DSTERm=Dm*(DEL*G/RNU**2)**(1./3.)
      end
C
      subroutine calc2004(dsi,psi,dt10,dt50,dt90,dm,egfi,pclay,ns)
      dimension dsi(100),psi(100),egfi(100),hulp(100)
      dm=0.
      dclay=0.000008
      do 10 i=1,ns
         hulp(i)=0.
         dm=dm+dsi(i)*psi(i)
         do 10 j=1,i
           hulp(i)=hulp(i)+psi(j)
   10 continue
      if(abs(hulp(ns)+pclay-1.).gt.1.e-3)then
         stop
      endif
c
c  dt50
c
      index=0
      do 20 i=1,ns-1
       if(hulp(i).le.0.5.and.0.5.le.hulp(i+1))then
         index=i
         goto 30
       endif
   20 continue
   30 continue
      if(index.eq.0)stop 'dt50 not found'
      dt50=dsi(index)+(0.5-hulp(i))/(hulp(i+1)-hulp(i))*
     *    (dsi(index+1)-dsi(index))
c
c  dt90
c
      index=0
      do 40 i=1,ns-1
       if(hulp(i).le.0.9.and.0.9.le.hulp(i+1))then
         index=i
         goto 50
       endif
   40 continue
   50 continue
      if(hulp(ns).lt.0.9)dt90=1.5*dsi(ns)
      if(hulp(ns).ge.0.9)dt90=dsi(index)+
     *(0.9-hulp(i))/(hulp(i+1)-hulp(i))*(dsi(index+1)-dsi(index))
c      if(index.eq.0)stop 'dt90 not found'
c      dt90=dsi(index)+(0.9-hulp(i))/(hulp(i+1)-hulp(i))*
c     *    (dsi(index+1)-dsi(index))
c
c  dt10
c
      index=0
      do 60 i=1,ns-1
       if(hulp(i).le.0.1.and.0.1.le.hulp(i+1))then
         index=i
         goto 70
       endif
   60 continue
   70 continue
c     smallest particle size is set to 8 um
      if(hulp(1).gt.0.1)dt10=dclay+0.1/hulp(1)*(dsi(1)-dclay)
      if(hulp(1).le.0.1)dt10=dsi(index)+
     *(0.1-hulp(i))/(hulp(i+1)-hulp(i))*(dsi(index+1)-dsi(index))
      if(dt10.le.dclay)dt10=dclay
c     if(index.eq.0)stop 'dt10 not found, use more fractions'
c      dt10=dsi(index)+(0.1-hulp(i))/(hulp(i+1)-hulp(i))*
c     *    (dsi(index+1)-dsi(index))
c
c egfi
c
      al19=log(19.)
      do 80 i=1,ns
      egfi(i)=(al19/(al19+log(dsi(i)/dt50)))**2
   80 continue
      end
C
C
C
      subroutine fdster2004(dster,taucr,taucr1,thetcr,pclay,g,dt50,rhos,
     *                  rho,FCR)
c     IF(DSTER.LE.4.)THETCR=.24/DSTER
      IF(DSTER.LE.4.)THETCR=0.115/(DSTER)**0.5
      IF(4. .LT.DSTER.AND.DSTER.LE.10.)THETCR=.14*DSTER**(-.64)
      IF(10..LT.DSTER.AND.DSTER.LE.20.)THETCR=.04*DSTER**(-.1 )
      IF(20..LT.DSTER.AND.DSTER.LE.150.)THETCR=.013*DSTER**(.29 )
      IF(DSTER.GT.150.)THETCR=.055
CC    Soulsby gives one single formula 
CC    THETCR=(0.24/DSTER)+0.055*(1.0-exp(-0.02*DSTER))
      dsand=0.000062
      dsilt=0.000032
      cmaxs=0.65
      fch1=(dsand/dt50)**1.5
      cmax=(dt50/dsand)*cmaxs
c     cmaxs=maximum bed concentration in case of sandy bottom (=0.65)
      if(cmax.lt.0.05)cmax=0.05
      if(cmax.gt.cmaxs)cmax=cmaxs
      fpack=cmax/cmaxs
      if(fch1.lt.1.)fch1=1.
      if(fpack.gt.1.)fpack=1.
      fclay=1.
      if(pclay.ge.0.)fclay=(1.+Pclay)**3.
cc      if(pclay.ge.0..and.dt50.ge.dsand)fclay=(1.+Pclay)**3.
      if(fclay.ge.2.)fclay=2.
      TAUCR=FCR*fpack*fch1*fclay*(RHOS-rho)*G*dt50*THETCR
      taucr1=taucr
      end
CC
      SUBROUTINE zeroin2004(X,Y,LOG,F,tolx2004,tp,vr,phi,g,pi,hd)
        LOGICAL LOG
      A=X
      FA=F(X,tp,vr,phi,g,pi,hd)
      B=Y
      X=Y
      FB=F(X,tp,vr,phi,g,pi,hd)
C INTERPOLATE
    1 C=A
      FC=FA
C EXTRAPOLATE:
    2 IF(ABS(FC).GE.ABS(FB))GOTO 19
      A=B
      FA=FB
      X=C
      B=C
      FB=FC
      C=A
      FC=FA
C END INTERCHANGE
   19 TOL=tolx2004(X)
      RM=(C+B)*.5
      IF(ABS(RM-B).LE.TOL)GOTO 200
      P=(B-A)*FB
      IF(P.LT.0.)GOTO 5
      Q=FA-FB
      GOTO 6
    5 Q=FB-FA
      P=-P
    6 A=B
      FA=FB
      IF(P.GT.ABS(Q)*TOL)GOTO 8
      IF(C.LE.B)GOTO 10
      B=B+TOL
      X=B
      GOTO 11
   10 B=B-TOL
      X=B
   11 GOTO 12
    8 IF(P.GE.(RM-B)*Q)GOTO 14
      B=P/Q+B
      X=B
      GOTO 12
   14 B=RM
      X=B
   12 FB=F(X,tp,vr,phi,g,pi,hd)
      IF(FC.LT.0.)GOTO 16
      IF(FB.LT.0.)GOTO 18
      GOTO 1
   18 GOTO 2
   16 IF(FB.GT.0.)GOTO 21
      GOTO 1
   21 GOTO 2
  200 Y=C
      IF(FC.LT.0.)GOTO 31
      IF(FB.GT.0.)GOTO 33
      LOG=.TRUE.
      GOTO 300
   33 LOG=.FALSE.
      GOTO 300
   31 IF(FB.LT.0.)GOTO 35
      LOG=.TRUE.
      GOTO 300
   35 LOG=.FALSE.
  300 RETURN
      END
CC
      function fwl2004(rls,tp,vr,phi,g,pi,hd)
C     common/niorez/tp,vr,phi,g,pi,hd
      fwl2004=(rls/tp-vr*cos(phi))**2-g*rls/2./pi*tanh(2.*pi*hd/rls)
      return
      end
CC
      function tolx2004(x)
      tolx2004=.001*abs(x)+1.e-5
      return
      end
C