C    <%Z%	%Y% %M%	%I% %G%>
      subroutine VDIFFUSE(idim,km,im,nx, ndet,extmode,
     &                    H, U, V, T, S, X , det, B,
     &                    Tdeep, Sdeep, Xdeep, detdeep,
     &                    ieast,iwest,inorth, isouth,
     &                    Area_h, Bu, Kappa, Nu, Ri, Ks, dtime)
      implicit none

c==> Arguments

      integer idim
      integer km
      integer im
      integer nx,ndet
      logical extmode
      real    H( idim, km)
      real    U( idim, km)
      real    V( idim, km)
      real    T( idim, km)
      real    S( idim, km)
      real    X( idim, km, 2, nx )
      real    det( idim, km, 2, ndet )
      real    B( idim, km)
      real    Tdeep
      real    Sdeep
      real    Xdeep( nx )
      real    detdeep( ndet )
      integer ieast( idim )
      integer iwest( idim )
      integer inorth( idim )
      integer isouth( idim )
      real    Area_h( idim )
      integer Bu( idim )
      real    hmax1
      real    dtime
      
c==> Output arrays for diagnostics

      real    Kappa( idim, km)
      real    Nu( idim, km )
      real    Ri( idim, km )
      real    Ks( idim, km )
      



c==> Local Parameters
      
      real MINRI            ! Richardson number at convection 
      real KAPPA_CONVECT    ! Kappa when convecting
      real NU_CONVECT       ! Viscosity when convecting
      real NU_ML            ! Nu at base of mixed layer (if isolated )
      real NU_BOTT          ! Nu at bottom of ocean (if extmode )
      integer NITS          ! Number of iterations
      
      parameter ( MINRI = 1.E-03 )
      real MINN2
      parameter ( MINN2 = 1.E-10 )
            
      parameter ( KAPPA_CONVECT = 1.E+03   )
      parameter ( NU_CONVECT    = 1.E-02   )
      
      
      parameter ( NU_ML         = 1.E-06   )
      parameter ( NU_BOTT       = 1.E-04   )
      
!     parameter ( NITS = 4 )
      parameter ( NITS = 1 )
      
 


c==> Local Variables

      real E1(idim, km)
      real A1(idim, km)
      real N2(idim, km)
      real Hu(idim, km)
      logical mix(idim, km)
      real    Area_u_r( idim )
      real Es(idim, km)
      real As(idim, km)
      real Alphs(idim)

      real Alpha(idim)
      real Wk1(idim)
      real Wk2(idim)
  
! Insitu conversion    
      real Z(idim,km)
      real Tbottom(idim) ! Abyssal in situ temperature
      real    Haa(idim)
      real Riu
      integer num_new
      integer nremain
      real Haas
      integer i,k,m,ntimes,k2
      real    dt2df,fact
      real    dh
      real    hmn
      real    hmx
      real    func1
      real    nu_bottom

 
c==> Functions    

      real FKAPPA
      real FKAPPA_SALT
      real FNU
      real INSITU
      real POTEMP
      real BYNCY
      
CFPP$ EXPAND ( FKAPPA ) R
CFPP$ EXPAND ( FNU ) R
CFPP$ EXPAND ( INSITU ) R
CFPP$ EXPAND ( POTEMP ) R
CFPP$ EXPAND ( BYNCY ) R

c---------------------------------------------------------------------
      if ( extmode ) then
        k2 = km-1
        do i=1,idim
         Kappa(i,km) = 0.0
         Ks(i,km)    = 0.0
         Nu(i,km)    = NU_BOTT
         Ri(i,km)    = 0.0
         mix(i,km)   = .FALSE.
        enddo
      else 
        k2 = km
      endif


      do k=1,k2
      do i=1,idim
       Ri(i,k) = 0.0
       N2(i,k) = 1.0
       Nu(i,k) = 0.0
       Kappa(i,k) = 0.0
       Ks(i,k)    = 0.0
       mix(i,k)= .FALSE.
      enddo
      enddo
      
      
      do i=1,idim
       Wk1(i) = 0.
      enddo
      do i=1,im
       Wk1(i) = Area_h(i) + Area_h(ieast(i))
      enddo
      do i=1,im
       Area_u_r(i) = Bu(i) / ( Wk1(i) + Wk1(inorth(i)) + 1.E-36 )
      enddo
      

      dt2df = 2.*dtime/NITS
      

! compute Z as the pressure of the center of each layer (not depth)
               
       do i=1,im
           Z(i,1) = 0.5 * H(i,1)
       enddo
       do k=2,km
        do i=1,im
         Z(i,k) = Z(i,k-1) + 0.5*( H(i,k) + H(i,k-1) )
        enddo
       enddo
       
       do i=1,im
        Tbottom(i)=INSITU(Tdeep,Sdeep,Z(i,km)+H(i,km)*.5)
       enddo
      
      nremain=NITS
      DO WHILE ( nremain .GT. 0 )
      

       num_new = 0
       do k=1,k2
        call VHRICH (idim, km, im, idim-im, k, Ri(1,k), N2(1,k),
     &               H,U,V,B,
     &               iwest,isouth )
        
          !  If we mixed on the last iteration, keep the diffusivity
          !  high, regardless of the present Richardson number.
          !  Also, update the array which flags high mixing (convection)
        
        do i=1,im
          ri(i,k)= amax1(Ri(i,k), MINRI)
          N2(i,k)= amax1(N2(i,k), MINN2)
          Kappa(i,k) = FKAPPA(Ri(i,k),N2(i,k))
          Ks(i,k)    = Kappa(i,k)
          Nu(i,k)    = FNU(Ri(i,k),N2(i,k))
        enddo ! i
       enddo ! k

         ! If we have not uncovered any new convection, just go
         ! ahead and do all the remaining iterations with the
         ! present values of the diffusivities.  This will save 
         ! time.  To further save time, we could gather all the
         ! points for which some_new applies, and re-do them,
         ! but that will await a later release
         
                
       if ( num_new .EQ.  0 ) then
        dt2df=dt2df*nremain
        nremain=0
       endif

       
       do k=1,km-1
        do i=1,im
         A1(i,k) = Kappa(i,k)*dt2df/(H(i,k) + H(i,k+1))
         As(i,k) = Ks(i,k)   *dt2df/(H(i,k) + H(i,k+1))
        enddo
       enddo
       k = km
       if ( extmode ) then
        do i=1,im
         A1(i,km) = 0.
         As(i,km) = 0.
        enddo
       else
        do i=1,im
         A1(i,km) = Kappa(i,km)*dt2df/H(i,km)
         As(i,km) = Ks(i,km)   *dt2df/H(i,km)
        enddo
       endif ! extmode


! Convert T to Insitu temp 
       do k=1,km
         do i=1,im
           T(i,k) = INSITU( T(i,k) , S(i,k) , Z(i,k) )
         enddo
       enddo       

       k=1
       do i=1,im
         Alpha(i) = H(i,k)*A1(i,k) / (H(i,k)+A1(i,k))
         Alphs(i) = H(i,k)*As(i,k) / (H(i,k)+As(i,k))
         E1(i,k)  =        A1(i,k) / (H(i,k)+A1(i,k)) 
         Es(i,k)  =        As(i,k) / (H(i,k)+As(i,k)) 
         T(i,k)   = T(i,k)* H(i,k) / (H(i,k)+A1(i,k))
         S(i,k)   = S(i,k)* H(i,k) / (H(i,k)+As(i,k))
        enddo
        do m=1,nx
         do i=1,im
          X(i,k,1,m) = X(i,k,1,m)*H(i,k)/(H(i,k)+A1(i,k))
         enddo
        enddo
        do m=1,ndet
         do i=1,im
          det(i,k,1,m) = det(i,k,1,m)*H(i,k)/(H(i,k)+A1(i,k))
         enddo
        enddo
        
       do k=2,km
       do i=1,im
         Haa(i)   = 1./( H(i,k) + A1(i,k) + Alpha(i) )
         Haas     = 1./( H(i,k) + As(i,k) + Alphs(i) )
         Alpha(i) = A1(i,k)*(H(i,k)+Alpha(i))*Haa(i)
         Alphs(i) = As(i,k)*(H(i,k)+Alphs(i))*Haas
         E1(i,k)  = A1(i,k)*Haa(i)
         Es(i,k)  = As(i,k)*Haas
         T(i,k)   = (T(i,k)*H(i,k) + A1(i,k-1)*T(i,k-1))*Haa(i)
         S(i,k)   = (S(i,k)*H(i,k) + As(i,k-1)*S(i,k-1))*Haas
        enddo
        
        do m=1,nx
         do i=1,im
           X(i,k,1,m) = 
     &          (X(i,k,1,m)*H(i,k) + A1(i,k-1)*X(i,k-1,1,m))*Haa(i)
         enddo
        enddo
        do m=1,ndet
         do i=1,im
           det(i,k,1,m) = 
     &          (det(i,k,1,m)*H(i,k) + A1(i,k-1)*det(i,k-1,1,m))*Haa(i)
         enddo
        enddo

       enddo !k

       if ( .NOT. extmode ) then
       do i=1,im
         T(i,km) = E1(i,km)*Tbottom(i) + T(i,km)
         S(i,km) = Es(i,km)*Sdeep      + S(i,km)
       enddo
       do m=1,nx
        do i=1,im
         X(i,km,1,m) = E1(i,km)*Xdeep(m) + X(i,km,1,m)
        enddo
       enddo
       do m=1,ndet
        do i=1,im
         det(i,km,1,m) = E1(i,km)*detdeep(m) + det(i,km,1,m)
        enddo
       enddo
       endif ! .NOT. extmode

       do k=km-1,1,-1
        do m=1,nx
         do i=1,im
           X(i,k,1,m) = E1(i,k)*X(i,k+1,1,m) + X(i,k,1,m)
         enddo
        enddo     
       enddo
       do k=km-1,1,-1
        do m=1,ndet
         do i=1,im
           det(i,k,1,m) = E1(i,k)*det(i,k+1,1,m) + det(i,k,1,m)
         enddo
        enddo     
       enddo
       do k=km-1,1,-1
         do i=1,im
           T(i,k) = E1(i,k)*T(i,k+1) + T(i,k)
           S(i,k) = Es(i,k)*S(i,k+1) + S(i,k)
         enddo
       enddo

       do k=1,km
        do i=1,im
          T(i,k) = POTEMP( T(i,k), S(i,k) , Z(i,k) )
        enddo
       enddo



c==>  Average viscosity  and layer thicknesses
C     Set Diffusion coefficients 
      
      do i = 1,idim
      Wk1(i)=0.
      Wk2(i)=0.
      enddo
        
      DO k=1,km
       do i=1,im
         Wk1(i) = ( Nu(ieast(i),k) + Nu(i,k) )*Area_h(i)
         Wk2(i) = (  H(ieast(i),k) +  H(i,k) )*Area_h(i)
       enddo ! i

       do i=1,im
          Nu(i,k)  = (Wk1(i)+Wk1(inorth(i)))*Area_u_r(i)
          Hu(i,k)  = (Wk2(i)+Wk2(inorth(i)))*Area_u_r(i) + ( 1.-Bu(i) )
       enddo ! i
         
      ENDDO ! k  Finish this loop to compute Hu at all k
     
     
       
      do k=1,km-1
       do i=1,im
         A1(i,k) = Nu(i,k)*dt2df/( Hu(i,k) + Hu(i,k+1)  )
       enddo
      enddo
      do i=1,im   ! Note that Nu(*,km) was set at top if extmode
        A1(i,km) = Nu(i,km)*dt2df/( Hu(i,km) )
      enddo



      k=1
      do i=1,im
         Alpha(i) = A1(i,k)*Hu(i,k) / (Hu(i,k)+A1(i,k))
         E1(i,k)  =         A1(i,k) / (Hu(i,k)+A1(i,k))
         U(i,k)   =  U(i,k)*Hu(i,k) / (Hu(i,k)+A1(i,k))
         V(i,k)   =  V(i,k)*Hu(i,k) / (Hu(i,k)+A1(i,k))
      enddo
      
       do k=2,km
        do i=1,im
         fact      = Bu(i)/(Hu(i,k) + A1(i,k) + Alpha(i))
         Alpha(i)  = A1(i,k)*(Hu(i,k)+Alpha(i))*fact
         E1(i,k)   = A1(i,k)*fact
         U(i,k) = (U(i,k)*Hu(i,k) + A1(i,k-1)*U(i,k-1))*fact
         V(i,k) = (V(i,k)*Hu(i,k) + A1(i,k-1)*V(i,k-1))*fact
       enddo
      enddo
      
      do k=km-1,1,-1
       do i=1,im
        U(i,k) = E1(i,k)*U(i,k+1) + U(i,k)
        V(i,k) = E1(i,k)*V(i,k+1) + V(i,k)
       enddo
      enddo
       
      

      do k=1,km
       do i=1,im
        B(i,k) = BYNCY( T(i,k), S(i,k), 0. )
       enddo
      end do
      
      nremain = nremain-1
      
      END DO ! while 

      end
