SUBROUTINE PROFTR
     1 (F,WFSURF,DT2,ttt,srho,sink,svel)
	parameter(pi=3.14159)
c
C     ---------------------------------------------------------------
c     use for sediment tracer
C     ---------------------------------------------------------------
C
      INCLUDE 'comblk.h'

	real fluxin,fluxout
      DIMENSION TTT(IM,JM,KB),WFSURF(IM,JM),sink(im,jm,kb)
      DIMENSION F(IM,JM,KB),DH(IM,JM) 
      dimension x3(im,jm),x4(im,jm),hs(im,jm)
      dimension TEMPDP(13),TEMPIS(13),srho(im,jm)

	DATA NTEMP/13/
	DATA TEMPDP/0.0002,0.0004,0.0008,0.002,0.004,0.008,0.02
     1   ,0.04,0.06,0.2,0.4,0.6,1.0/
	DATA TEMPIS/6.3e-3,7.6e-3,9.e-3,1.2e-2,1.5e-2,1.8e-2,
     1   2.7e-2,3.6e-2,4.3e-2,8.e-2,1.3e-1,1.7e-1,2.3e-1/

c      EQUIVALENCE (TPS,DH)
C
      UMOLPR=UMOL
	srho0=10.0

c
C     ----------------------------------------------------------------
c     Resuspension constant as a function of sediment density and a
c     velocity scale, srho and svelo, see notes:
C     ----------------------------------------------------------------
c
c      srho = 1.
c      svel = 1.e-4
c      svel = 1.e-5
c      rm   = srho * svel
c
C     ----------------------------------------------------------------
c     critical bottom shear stress for erosion ttte
c     critical bottom shear stress for deposition tttd 
C     ----------------------------------------------------------------
c
      ttte=0.02
      tttd=0.02
c

c
C     ----------------------------------------------------------------
c     compute bottom shear stress at k=kbm1 at tracer point: METHOD 2 
C     ----------------------------------------------------------------
c

         DO J=2,JMM1
         DO I=1,IMM1

         rrho=(rho(i,j,kbm1) + 1.025)*1000.0

         tttx=rrho*0.5*(wubot(i,j)+wubot(i+1,j))         
         ttty=rrho*0.5*(wvbot(i,j)+wvbot(i+1,j))
         ttt(i,j,kbm1)=sqrt(tttx**2.+ttty**2.)

         end do
         end do

c
C     ******************************************************************
C                                                                      *
C     THE FOLLOWING SECTION SOLVES THE EQUATION  DTI2*(KH*F')'-F=-FB   *
C                                                                      *
C     ******************************************************************
c
      DO 10 J=1,JM
      DO 10 I=1,IM
         DH(I,J)=H(I,J)+ETF(I,J)
   10 CONTINUE

      DO 20 K=2,KBM1
      DO 20 J=1,JM
      DO 20 I=1,IM
         A(I,J,K-1)=-DT2*(KH(I,J,K)+UMOLPR)/(DZ(K-1)*DZZ(K-1)*DH(I,J)
     1             *DH(I,J))
         C(I,J,K)=-DT2*(KH(I,J,K)+UMOLPR)/(DZ(K)*DZZ(K-1)*DH(I,J)
     1           *DH(I,J))
   20 CONTINUE
c
C     ---------------------------------------------------------------
c     k=1, surface condition
C     ---------------------------------------------------------------
c
      DO 500 J=1,JM
      DO 500 I=1,IM
         EE(I,J,1)=A(I,J,1)/(A(I,J,1)-1.0)
      GG(I,J,1)=-DT2*WFSURF(I,J)/(-DZ(1)*DH(I,J))-F(I,J,1)
      GG(I,J,1)=GG(I,J,1)/(A(I,J,1)-1.0)
c         GG(I,J,1)=-F(I,J,1)/(A(I,J,1)-1.0)
  500 CONTINUE
c
C     --------------------------------------------------------------
c     k=2-kbm2, subsurface condition
C     --------------------------------------------------------------
c
      DO 101 K=2,KBM2
      DO 101 J=1,JM
      DO 101 I=1,IM
         GG(I,J,K)=1./(A(I,J,K)+C(I,J,K)*(1.-EE(I,J,K-1))-1.)
         EE(I,J,K)=A(I,J,K)*GG(I,J,K)
         GG(I,J,K)=(C(I,J,K)*GG(I,J,K-1)-F(I,J,K))*GG(I,J,K)
  101 CONTINUE
c
C     --------------------------------------------------------------
C     k=kbm1, bottom condition/BOTTOM ADIABATIC B.C.
C     --------------------------------------------------------------

      DO 102 J=2,JMM1
      DO 102 I=1,IM
c
C     ----------------------------------------------------------------
c        volumn of bottom box
C     ----------------------------------------------------------------
c
         vol=dx(i,j)*dy(i,j)*dh(i,j)*dz(kbm1)*fsm(i,j)
c
C     ----------------------------------------------------------------
c        flux into water column if shear is larger than critical shear
c        for resuspension/erosion, (-) for influx
C     ----------------------------------------------------------------
c
         if(ttt(i,j,kbm1).gt.ttte) then
            fluxin = srho0*svel*(-1.)*(ttt(i,j,kbm1)/ttte - 1.0)
            ttt(i,j,1) = -dt2*fluxin
	 srho(i,j)=srho(i,j)-ttt(i,j,1)
            x3(i,j)=-ttt(i,j,1)/(DZ(kbm1)*DH(I,J))
         else
            fluxin = 0.E0
             ttt(i,j,1) = 0.0E0
            x3(i,j)=0.
         end if
c
C     ----------------------------------------------------------------
c        flux out of water column if shear is less than critical shear         
c        for deposition
C     ----------------------------------------------------------------
c

         if(ttt(i,j,kbm1).lt.tttd) then

        fluxout=F(I,J,KBM1)*sink(i,j,kbm1)*(1.0-ttt(i,j,kbm1)/tttd)
             ttt(i,j,2) = -dt2*fluxout
       	     srho(i,j)=srho(i,j)+ttt(i,j,2)
             x4(i,j)=-ttt(i,j,2)/(DZ(kbm1)*DH(I,J)) 
             F(i,j,kbm1) = F(i,j,kbm1) + x4(i,j)

         else
             fluxout=0.E0
             ttt(i,j,2) = 0.0E0
            x4(i,j)=0.
         end if
c
         F(I,J,KBM1)=((C(I,J,KBM1)*GG(I,J,KBM2)-F(I,J,KBM1)+x3(i,j))
     +              /(C(I,J,KBM1)*(1.-EE(I,J,KBM2))-1.))
c
         x3(i,j) = vol * x3(i,j) * (-1.)
         x4(i,j) = vol * x4(i,j) * (-1.)
         
  102 CONTINUE
c
C     --------------------------------------------------------------
c     compute new flux
C     --------------------------------------------------------------
c
      DO 105 K=2,KBM1
         KI=KB-K
         DO 105 J=2,JMM1
         DO 105 I=1,IM
         F(I,J,KI)=(EE(I,J,KI)*F(I,J,KI+1)+GG(I,J,KI))
  105 CONTINUE

      RETURN
      END