      SUBROUTINE smagor(ib,ie,jb,je,kb,ke)

      IMPLICIT NONE

      INCLUDE 'param.txt'
      INCLUDE 'common.txt'

      INTEGER ib,ie,jb,je,kb,ke,ip,im,jp,jm,kp,km

      REAL athird
c-----------------------------------------------------------------------

c *** divergence of velocity ***

      DO k=kb,ke
         
         kp=k+1
         km=k-1

         DO j=jb,je

            jp=j+1
            jm=j-1

            DO i=ib,ie

               ip=i+1
               im=i-1

               div(i,j,k) = 
     1    ( ru(i)*unew(i,j,k) - ru(im)*unew(im,j,k) ) / 
     1              (rp(i)*(ru(i)-ru(im)))
     2 +  ( vnew(i,j,k) - vnew(i,jm,k) )              / 
     2              (rp(i)*dphi)        
     3 +  ( wnew(i,j,k) - wnew(i,j,km) )              /
     3              (zw(k) - zw(km))


               div(i,j,k) = 0.0


            END DO
         END DO
      END DO

c *** deviatoric strain rate tensor ***

      athird = 1.0/3.0

      DO k=kb,ke

         kp=k+1
         km=k-1

         DO j=jb,je

            jp=j+1
            jm=j-1

c *** Exclude i=1: considered i=1 apart at the end of the 300 loop

            DO i= 2,ie

               ip=i+1
               im=i-1

               drr(i,j,k) = ( unew(i,j,k)-unew(im,j,k) ) / (ru(i)-ru(im
     $              )) - athird*div(i,j,k)             

               doo(i,j,k) = ( (vnew(i,j,k)-vnew(i,jm,k)) / dphi +(unew(i
     $              ,j,k)+unew(im,j,k)) / 2.0 ) / rp(i) - athird*div(i,j
     $              ,k)

               dzz(i,j,k) = ( wnew(i,j,k) - wnew(i,j,km) ) / (zw(k)
     $              -zw(km)) - athird*div(i,j,k)


c *** r-theta (theta-r) strain ***


               dro(i,j,k) = 0.5 * 0.25 * (
     & ( ru(i) * ( vnew(ip,j,k)/rp(ip)  - vnew(i,j,k)/rp(i)   ) /
     &           ( Rp(ip) - Rp(i) )                             +
     &           ( unew(i,jp,k)         - unew(i,j,k)         ) /
     &           ( ru(i) * dphi )                               )
     &      +
     & ( ru(i) * ( vnew(ip,jm,k)/rp(ip) - vnew(i,jm,k)/rp(i)  ) /
     &           ( rp(ip) - rp(i) )                             +
     &           ( unew(i,j,k)          - unew(i,jm,k)        ) /
     &           ( ru(i) * dphi )                               )
     &      +
     & ( ru(im)* ( vnew(i,j,k)/rp(i)    - vnew(im,j,k)/rp(im) ) /
     &           ( rp(i) - rp(im) )                             +
     &           ( unew(im,jp,k)        - unew(im,j,k)        ) /
     &           ( ru(im) * dphi )                              )
     &      +
     & ( ru(im)* ( vnew(i,jm,k)/rp(i)   - vnew(im,jm,k)/rp(im)) /
     &           ( rp(i) - rp(im) )                             +
     &           ( unew(im,j,k)         - unew(im,jm,k)       ) /
     &           ( ru(im) * dphi )                              )
     &                       )

c *** theta-z (z-theta) strain ***


               doz(i,j,k) = 0.5 * 0.25 * (
     & (  ( wnew(i,jp,k)  - wnew(i,j,k)   ) / ( rp(i) * dphi )  +
     &    ( vnew(i,j,kp)  - vnew(i,j,k)   ) / ( zp(kp)-zp(k) )  )
     &      +
     & (  ( wnew(i,j,k)   - wnew(i,jm,k)  ) / ( rp(i) * dphi )  +
     &    ( vnew(i,jm,kp) - vnew(i,jm,k)  ) / ( zp(kp)-zp(k) )  )
     &      +
     & (  ( wnew(i,jp,km) - wnew(i,j,km)  ) / ( rp(i) * dphi )  +
     &    ( vnew(i,j,k)   - vnew(i,j,km)  ) / ( zp(k)-zp(km) )  )
     &      +
     & (  ( wnew(i,j,km)  - wnew(i,jm,km) ) / ( rp(i) * dphi )  +
     &    ( vnew(i,jm,k)  - vnew(i,jm,km) ) / ( zp(k)-zp(km) )  )
     &                       )


c *** r-z (z-r) strain ***


               drz(i,j,k) = 0.5 * 0.25 * (
     & (  ( wnew(ip,j,k)  - wnew(i,j,k)   ) / ( rp(ip)-rp(i) )  +
     &    ( unew(i,j,kp)  - unew(i,j,k)   ) / ( zp(kp)-zp(k) )  )
     &      +
     & (  ( wnew(ip,j,km) - wnew(i,j,km)  ) / ( rp(ip)-rp(i) )  +
     &    ( unew(i,j,k)   - unew(i,j,km)  ) / ( zp(k)-zp(km) )  )
     &      +
     & (  ( wnew(i,j,k)   - wnew(im,j,k)  ) / ( rp(i)-rp(im) )  +
     &    ( unew(im,j,kp) - unew(im,j,k)  ) / ( zp(kp)-zp(k) )  )
     &      +
     & (  ( wnew(i,j,km)  - wnew(im,j,km) ) / ( rp(i)-rp(im) )  +
     &    ( unew(im,j,k)  - unew(im,j,km) ) / ( zp(k)-zp(km) )  )
     &                       )

            END DO

c *** Consider i=1 as a special case because some stresses become zero
c *** at the centerline:

      im = 0
      i  = 1
      ip = 2

            drr(i,j,k) = ( unew(i,j,k)-unew(im,j,k) ) / (ru(i)-ru(im
     $         )) - athird*div(i,j,k)             

            doo(i,j,k) = ( (vnew(i,j,k)-vnew(i,jm,k)) / dphi +(unew(i
     $           ,j,k)+unew(im,j,k)) / 2.0 ) / rp(i) - athird*div(i,j
     $           ,k)

            dzz(i,j,k) = ( wnew(i,j,k) - wnew(i,j,km) ) / (zw(k)
     $           -zw(km)) - athird*div(i,j,k)


c *** r-theta (theta-r) strain ***


            dro(i,j,k) = 0.5 * 0.25 * (
     & ( ru(i) * ( vnew(ip,j,k)/rp(ip)  - vnew(i,j,k)/rp(i)   ) /
     &           ( rp(ip) - rp(i) )                             +
     &           ( unew(i,jp,k)         - unew(i,j,k)         ) /
     &           ( ru(i) * dphi )                               )
     &      +
     & ( ru(i) * ( vnew(ip,jm,k)/rp(ip) - vnew(i,jm,k)/rp(i)  ) /
     &           ( rp(ip) - rp(i) )                             +
     &           ( unew(i,j,k)          - unew(i,jm,k)        ) /
     &           ( ru(i) * dphi )                               )
cR=0 &      +
cR=0 & ( ru(im)* ( vnew(i,j,k)/rp(i)    - vnew(im,j,k)/rp(im) ) /
cR=0 &           ( rp(i) - rp(im) )                             +
cR=0 &           ( unew(im,jp,k)        - unew(im,j,k)        ) /
cR=0 &           ( ru(im) * dphi )                              )
cR=0 &      +
cR=0 & ( ru(im)* ( vnew(i,jm,k)/rp(i)   - vnew(im,jm,k)/rp(im)) /
cR=0 &           ( rp(i) - rp(im) )                             +
cR=0 &           ( unew(im,j,k)         - unew(im,jm,k)       ) /
cR=0 &           ( ru(im) * dphi )                              )
     &                       )


c *** theta-z (z-theta) strain ***


            doz(i,j,k) = 0.5 * 0.25 * (
     & (  ( wnew(i,jp,k)  - wnew(i,j,k)   ) / ( rp(i) * dphi )  +
     &    ( vnew(i,j,kp)  - vnew(i,j,k)   ) / ( zp(kp)-zp(k) )  )
     &      +
     & (  ( wnew(i,j,k)   - wnew(i,jm,k)  ) / ( rp(i) * dphi )  +
     &    ( vnew(i,jm,kp) - vnew(i,jm,k)  ) / ( zp(kp)-zp(k) )  )
     &      +
     & (  ( wnew(i,jp,km) - wnew(i,j,km)  ) / ( rp(i) * dphi )  +
     &    ( vnew(i,j,k)   - vnew(i,j,km)  ) / ( zp(k)-zp(km) )  )
     &      +
     & (  ( wnew(i,j,km)  - wnew(i,jm,km) ) / ( rp(i) * dphi )  +
     &    ( vnew(i,jm,k)  - vnew(i,jm,km) ) / ( zp(k)-zp(km) )  )
     &                       )

c *** r-z (z-r) strain ***


            drz(i,j,k) = 0.5 * 0.25 * (
     & (  ( wnew(ip,j,k)  - wnew(i,j,k)   ) / ( rp(ip)-rp(i) )  +
     &    ( unew(i,j,kp)  - unew(i,j,k)   ) / ( zp(kp)-zp(k) )  )
     &      +
     & (  ( wnew(ip,j,km) - wnew(i,j,km)  ) / ( rp(ip)-rp(i) )  +
     &    ( unew(i,j,k)   - unew(i,j,km)  ) / ( zp(k)-zp(km) )  )
cR=0 &      +
cR=0 & (  ( wnew(i,j,k)   - wnew(im,j,k)  ) / ( rp(i)-rp(im) )  +
cR=0 &    ( unew(im,j,kp) - unew(im,j,k)  ) / ( dz )            )
cR=0 &      +
cR=0 & (  ( wnew(i,j,km)  - wnew(im,j,km) ) / ( rp(i)-rp(im) )  +
cR=0 &    ( unew(im,j,k)  - unew(im,j,km) ) / ( dz )            )
     &                       )

         END DO

      END DO


c... compute |d_i,j| = sqrt(2*d_(i,j)*d_(i,j)) ... 
c... and mu_t = C_s * delta^2 * |d_i,j| = lmix(i,k)*magdij(i,j,k) 

      DO k=kb,ke
         DO j=jb,je
            DO i=ib,ie

               magdij(i,j,k) = SQRT( 2*drr(i,j,k)**2 + 2*doo(i,j,k)**2 +
     $              2*dzz(i,j,k)**2 + 4*dro(i,j,k)**2 + 4*doz(i,j,k)**2
     $              + 4*drz(i,j,k)**2 )
               
               mut(i,j,k) = lmix(i,k) * magdij(i,j,k)

            END DO
         END DO
      END DO

c*************************************************************
c
c     Set boundary boundary condition for K-closure factors.
c  
c     probably don't need to do this there because it's done
c     in the routine viscous, but what the hell.
c
c*************************************************************

      DO k=1,kmax

c *** Theta-direction ***

         DO i=1,imax

            mut(i,0,k)  = mut(i,jmax,k)
            mut(i,j1,k) = mut(i,1,k)

         END DO

c *** R-direction ***

         DO j=0,j1

            mut(0,j,k)  =   mut(1,j,k)
            mut(i1,j,k) = - mut(imax,j,k)

         END DO

      END DO

c *** Z-direction ***

      DO j=0,j1
         DO i=0,i1

            mut(i,j,0)  = mut(i,j,kmax)
            mut(i,j,k1) = mut(i,j,1)

         END DO
      END DO
      
      RETURN
      END
