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

      IMPLICIT NONE

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

      INTEGER ib,ie,jb,je,kb,ke

      INTEGER l

      REAL uht(0:i1,0:j1,0:k1)
      REAL vht(0:i1,0:j1,0:k1)
      REAL wht(0:i1,0:j1,0:k1)

      REAL up(0:i1,0:j1,0:k1)
      REAL vp(0:i1,0:j1,0:k1)
      REAL wp(0:i1,0:j1,0:k1)

      REAL abr(6,0:i1,0:j1,0:k1)
      REAL magabr(0:i1,0:j1,0:k1)
      REAL aht(6,0:i1,0:j1,0:k1)
      REAL magaht(0:i1,0:j1,0:k1)
      REAL bbr(0:i1,0:j1,0:k1)
      REAL bht(0:i1,0:j1,0:k1)
      REAL mkldkl(0:i1,0:j1,0:k1)
      REAL lijdij(0:i1,0:j1,0:k1)

      REAL im,jm,km
      REAL tmp,sum1,sum2
c--------------------------------------------------------------------

c... first test-filter velocity field ...

      CALL filter(2,unew,uht,ru,i1,j1,k1,ib,ie,jb,je,kb,ke)
      CALL filter(2,vnew,vht,rp,i1,j1,k1,ib,ie,jb,je,kb,ke)
      CALL filter(2,wnew,wht,rp,i1,j1,k1,ib,ie,jb,je,kb,ke)





c      do i=1,imax,4
c         do j=1,jmax,4
c            do k=1,2

c               print*,1,i,j,k,unew(i,j,k),uht(i,j,k)
c               print*,2,i,j,k,vnew(i,j,k),vht(i,j,k)
c               print*,3,i,j,k,wnew(i,j,k),wht(i,j,k)

c            end do
c         end do
c      end do
c      stop
               




c... now set the boundary values for the test-filtered velocity field ...

      CALL bound(uht,vht,wht,tnew,ynew,rhonew,time,ijet,dt,zp,zw,0
     $     ,ru,rp)

c... compute d^bar_ij and |d^bar_ij| ...

      CALL deviatoric(abr,magabr,unew,vnew,wnew,ru,rp,dphi,zw,zp,i1,j1
     $     ,k1,ib,ie,jb,je,kb,ke)





c      print*,ijet
c      do i=1,imax
c         do j=1,jmax,4
c            write(*,20) i,j,abr(1,i,j,1),abr(2,i,j,1),abr(3,i,j,1)
c            write(*,20) i,j,abr(4,i,j,1),abr(5,i,j,1),abr(6,i,j,1)
c            write(*,21) i,j,magabr(i,j,1)
c 20         format(2(i4),3(f12.5))
c 21         format(2(i4),f12.5)
c         end do
c      end do
c      stop










c... compute d^bar^hat_ij and |d^bar^hat_ij| ...

      CALL deviatoric(aht,magaht,uht,vht,wht,ru,rp,dphi,zw,zp,i1,j1,k1
     $     ,ib,ie,jb,je,kb,ke)



c      do i=1,imax
c         do j=1,jmax,4
c            print*,i,j,aht(1,i,j,2),aht(2,i,j,2),aht(3,i,j,2)
c            print*,i,j,aht(4,i,j,2),aht(5,i,j,2),aht(6,i,j,2)
c            write(*,10) i,j,aht(2,i,j,1),abr(2,i,j,1),aht(3,i,j,1),abr(3
c     $           ,i,j,1)
c            write(*,10) i,j,aht(5,i,j,1),abr(5,i,j,1),aht(6,i,j,1),abr(5
c     $           ,i,j,1)
c 10         format(2(i4),4(f12.5))
c         end do
c      end do
c      stop





c... initialize M_kl * d^bar_kl and L_ij * d^bar_ij to zero ...

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

               mkldkl(i,j,k) = 0.0
               lijdij(i,j,k) = 0.0

            END DO
         END DO
      END DO

c... compute M_kl * d^bar_kl ...

      DO l=1,6

         DO i=1,imax
            DO j=0,j1
               DO k=1,kmax
                  bbr(i,j,k) = magabr(i,j,k)*abr(l,i,j,k)
               END DO
            END DO
         END DO

         CALL filter(2,bbr,bht,rp,i1,j1,k1,ib+1,ie-1,jb,je,kb+1,ke-1)





c      if(l .eq. 3)then
c      do i=1,imax
c         do j=1,jmax,4
c            write(*,30) i,j,bbr(i,j,2),bht(i,j,2)
c 30         format(2(i4),2(f12.5))
c         end do
c      end do
c      stop
c      end if







         tmp = 1.0
         IF(l .GT. 3) tmp = 2.0
         DO i=2,imax-1
            DO j=1,jmax
               DO k=2,kmax-1
               
                  mkldkl(i,j,k) = mkldkl(i,j,k) + tmp*((delht(i,k)
     $                 /delbr(i,k))*magaht(i,j,k)*aht(l,i,j,k)*abr(l,i,j
     $                 ,k) - bht(i,j,k)*aht(l,i,j,k))

               END DO
            END DO
         END DO






c      if(l .eq. 3)then
c      do i=2,imax-1,4
c         do j=1,jmax,4
c            do k=2,10
c               write(*,30) i,j,k,mkldkl(i,j,k)
c 30            format(3(i4),1(f12.5))
c            end do
c         end do
c      end do
c      stop
c      end if






         
      END DO

c... compute L_ij * d^bar_ij ...

c... interpolate velocities onto pressure nodes ...

      DO i=1,imax
         DO j=1,jmax
            DO k=1,kmax

               up(i,j,k) = 0.5*(unew(i,j,k) + unew(i-1,j,k))

               vp(i,j,k) = 0.5*(vnew(i,j,k) + vnew(i,j-1,k))

               wp(i,j,k) = 0.5*(wnew(i,j,k) + wnew(i,j,k-1))

            END DO
         END DO
      END DO

c... apply boundary conditions in theta direction ...

      DO i=1,imax
         DO k=1,kmax

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

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

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

         END DO
      END DO

c... test filter velocity fields ...

      CALL filter(2,up,uht,rp,i1,j1,k1,ib+1,ie-1,jb,je,kb+1,ke-1)
      CALL filter(2,vp,vht,rp,i1,j1,k1,ib+1,ie-1,jb,je,kb+1,ke-1)
      CALL filter(2,wp,wht,rp,i1,j1,k1,ib+1,ie-1,jb,je,kb+1,ke-1)





c      do i=2,imax-1,6
c         do j=1,jmax,6
c            do k=1,kmax-1
c               
c               write(*,50) i,j,k,up(i,j,k),uht(i,j,k),vp(i,j,k),vht(i,j
c     $              ,k),wp(i,j,k),wht(i,j,k)
c 50            format(3(i4),6(f12.5))
c
c            end do
c         end do
c      end do
c      stop















c... compute u^bar^hat_i * u^bar^hat_j ...

      DO i=1,imax
         DO j=0,j1
            DO k=1,kmax

               aht(1,i,j,k) = uht(i,j,k)*uht(i,j,k)
               aht(2,i,j,k) = vht(i,j,k)*vht(i,j,k)
               aht(3,i,j,k) = wht(i,j,k)*wht(i,j,k)
               aht(4,i,j,k) = uht(i,j,k)*vht(i,j,k)
               aht(5,i,j,k) = wht(i,j,k)*vht(i,j,k)
               aht(6,i,j,k) = uht(i,j,k)*wht(i,j,k)

            END DO
         END DO
      END DO

c... compute u^bar_ij * u^bar_ij ...                

c... r-r component ...

      DO i=1,imax
         DO j=0,j1
            DO k=1,kmax

               bbr(i,j,k) = up(i,j,k)*up(i,j,k)

            END DO
         END DO
      END DO
      CALL filter(2,bbr,bht,rp,i1,j1,k1,ib+1,ie-1,jb,je,kb+1,ke-1)
      DO i=2,imax-1
         DO j=1,jmax
            DO k=2,kmax-1
               
               lijdij(i,j,k) = lijdij(i,j,k) + (bht(i,j,k) - aht(1,i,j,k
     $              ))*abr(1,i,j,k)

            END DO
         END DO
      END DO

c... theta-theta component ...

      DO i=1,imax
         DO j=0,j1
            DO k=1,kmax

               bbr(i,j,k) = vp(i,j,k)*vp(i,j,k)

            END DO
         END DO
      END DO
      CALL filter(2,bbr,bht,rp,i1,j1,k1,ib+1,ie-1,jb,je,kb+1,ke-1)
      DO i=2,imax-1
         DO j=1,jmax
            DO k=2,kmax-1
               
               lijdij(i,j,k) = lijdij(i,j,k) + (bht(i,j,k) - aht(2,i,j,k
     $              ))*abr(2,i,j,k)

            END DO
         END DO
      END DO

c... z-z component ...

      DO i=1,imax
         DO j=0,j1
            DO k=1,kmax

               bbr(i,j,k) = wp(i,j,k)*wp(i,j,k)

            END DO
         END DO
      END DO
      CALL filter(2,bbr,bht,rp,i1,j1,k1,ib+1,ie-1,jb,je,kb+1,ke-1)
      DO i=2,imax-1
         DO j=1,jmax
            DO k=2,kmax-1
               
               lijdij(i,j,k) = lijdij(i,j,k) + (bht(i,j,k) - aht(3,i,j,k
     $              ))*abr(3,i,j,k)

            END DO
         END DO
      END DO

c... r-theta component ...

      DO i=1,imax
         DO j=0,j1
            DO k=1,kmax

               bbr(i,j,k) = up(i,j,k)*vp(i,j,k)

            END DO
         END DO
      END DO
      CALL filter(2,bbr,bht,rp,i1,j1,k1,ib+1,ie-1,jb,je,kb+1,ke-1)
      DO i=2,imax-1
         DO j=1,jmax
            DO k=2,kmax-1
               
               lijdij(i,j,k) = lijdij(i,j,k) + 2.0*(bht(i,j,k) - aht(4,i
     $              ,j,k))*abr(4,i,j,k)

            END DO
         END DO
      END DO
               
c... z-theta component ...

      DO i=1,imax
         DO j=0,j1
            DO k=1,kmax

               bbr(i,j,k) = wp(i,j,k)*vp(i,j,k)

            END DO
         END DO
      END DO
      CALL filter(2,bbr,bht,rp,i1,j1,k1,ib+1,ie-1,jb,je,kb+1,ke-1)
      DO i=2,imax-1
         DO j=1,jmax
            DO k=2,kmax-1
               
               lijdij(i,j,k) = lijdij(i,j,k) + 2.0*(bht(i,j,k) - aht(5,i
     $              ,j,k))*abr(5,i,j,k)

            END DO
         END DO
      END DO
               
c... r-z component ...

      DO i=1,imax
         DO j=0,j1
            DO k=1,kmax

               bbr(i,j,k) = up(i,j,k)*wp(i,j,k)

            END DO
         END DO
      END DO
      CALL filter(2,bbr,bht,rp,i1,j1,k1,ib+1,ie-1,jb,je,kb+1,ke-1)
      DO i=2,imax-1
         DO j=1,jmax
            DO k=2,kmax-1
               
               lijdij(i,j,k) = lijdij(i,j,k) + 2.0*(bht(i,j,k) - aht(6,i
     $              ,j,k))*abr(6,i,j,k)

            END DO
         END DO
      END DO
               
c... average in theta direction ...

      DO i=2,imax-1
         DO k=2,kmax-1
            
            sum1 = 0.0
            sum2 = 0.0
            DO j=1,jmax

               sum1 = sum1 + lijdij(i,j,k)
               sum2 = sum2 + mkldkl(i,j,k)

            END DO
            lijdij(i,0,k) = sum1/jmax
            mkldkl(i,0,k) = sum2/jmax

         END DO
      END DO

      DO i=2,imax-1
         DO k=2,kmax-1



c            write(*,100) i,k,lijdij(i,0,k),mkldkl(i,0,k)
c 100        format(2(i4),2(f12.5))


            
            lmix(i,k) = -0.5*lijdij(i,0,k)/mkldkl(i,0,k)
            
         END DO
      END DO


c      stop





      DO k=2,kmax-1

         lmix(1,k) = lmix(2,k)
         lmix(imax,k) = lmix(imax-1,k)

      END DO
      lmix(1,1) = lmix(1,2)
      lmix(imax,kmax) = lmix(imax,kmax-1)

      DO i=2,imax-1

         lmix(i,1) = lmix(i,2)
         lmix(i,kmax) = lmix(i,kmax-1)

      END DO
      lmix(1,kmax) = lmix(1,kmax-1)
      lmix(imax,1) = lmix(imax-1,1)


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

               mut(i,j,k) = lmix(i,k) * magabr(i,j,k)

            END DO
         END DO
      END DO

c... set boundary values for mut ...


      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)
c            mut(i1,j,k) = - mut(imax,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


c--------------------------------------------------------------------
c--------------------------------------------------------------------
      SUBROUTINE filter(dim,fbr,fht,rf,i1,j1,k1,ib,ie,jb,je,kb,ke)

      IMPLICIT NONE

      INTEGER dim,i1,j1,k1,ib,ie,jb,je,kb,ke
      INTEGER i,j,k,im,ip,jm,jp,km,kp

      REAL fbr(0:i1,0:j1,0:k1)
      REAL fht(0:i1,0:j1,0:k1)
      REAL rf(0:i1)

      REAL a1,a2,a3,a4
c--------------------------------------------------------------------

c... filter in theta-z plane ...

      IF(dim .EQ. 2)THEN

         a1 = 1.0  / 9.0
         a2 = 4.0  / 9.0
         a3 = 16.0 / 9.0

         DO i=ib,ie
            DO j=jb,je
               
               jm = j-1
               jp = j+1

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

                  fht(i,j,k) = 0.25*(
     $        a1*fbr(i,jm,km) + a2*fbr(i,jm,k) + a1*fbr(i,jm,kp) 
     $      + a2*fbr(i,j,km)  + a3*fbr(i,j,k)  + a2*fbr(i,j,kp)
     $      + a1*fbr(i,jp,km) + a2*fbr(i,jp,k) + a1*fbr(i,jp,kp)
     $                               )

               END DO
            END DO
         END DO

      END IF

c... filter in all three directions ...

      IF(dim .EQ. 3)THEN 

         a1 = 1.0  / 27.0
         a2 = 4.0  / 27.0
         a3 = 16.0 / 27.0
         a4 = 64.0 / 27.0

         DO i=ib,ie

            im = i-1
            ip = i+1

            DO j=jb,je

               jm = j-1
               jp = j+1

               DO k=kb,ke

                  km = k-1
                  kp = k+1

                  fht(i,j,k) = 1.0/(8.0*rf(i))*(
     1      a1*rf(im)*fbr(im,jm,km) + a2*rf(i)*fbr(i,jm,km)
     1    + a1*rf(ip)*fbr(ip,jm,km)
     2    + a2*rf(im)*fbr(im,jm,k)  + a3*rf(i)*fbr(i,jm,k)
     2    + a2*rf(ip)*fbr(ip,jm,k)
     3    + a1*rf(im)*fbr(im,jm,kp) + a2*rf(i)*fbr(i,jm,kp)
     3    + a1*rf(ip)*fbr(ip,jm,kp)
     4    + a2*rf(im)*fbr(im,j,km)  + a3*rf(i)*fbr(i,j,km)
     4    + a2*rf(ip)*fbr(ip,j,km)
     5    + a3*rf(im)*fbr(im,j,k)   + a4*rf(i)*fbr(i,j,k)
     5    + a3*rf(ip)*fbr(ip,j,k)
     6    + a2*rf(im)*fbr(im,j,kp)  + a3*rf(i)*fbr(i,j,kp)
     6    + a2*rf(ip)*fbr(ip,j,kp)
     7    + a1*rf(im)*fbr(im,jp,km) + a2*rf(i)*fbr(i,jp,km)
     7    + a1*rf(ip)*fbr(ip,jp,km)
     8    + a2*rf(im)*fbr(im,jp,k)  + a3*rf(i)*fbr(i,jp,k)
     8    + a2*rf(ip)*fbr(ip,jp,k)
     9    + a1*rf(im)*fbr(im,jp,kp) + a2*rf(i)*fbr(i,jp,kp)
     9    + a1*rf(ip)*fbr(ip,jp,kp)
     $                                          )

               END DO
            END DO
         END DO

      END IF

      RETURN
      END

c---------------------------------------------------------------------
c---------------------------------------------------------------------
      SUBROUTINE deviatoric(dij,magdij,u,v,w,ru,rp,dphi,zw,zp,i1,j1,k1
     $     ,ib,ie,jb,je,kb,ke)

      IMPLICIT NONE

      INTEGER i1,j1,k1,ib,ie,jb,je,kb,ke
      INTEGER i,j,k,im,ip,jm,jp,km,kp

      INTEGER l

      REAL dij(6,0:i1,0:j1,0:k1)
      REAL magdij(0:i1,0:j1,0:k1)
      REAL u(0:i1,0:j1,0:k1)
      REAL v(0:i1,0:j1,0:k1)
      REAL w(0:i1,0:j1,0:k1)
      REAL ru(0:i1),rp(0:i1),dphi,zw(0:k1),zp(0:k1)

      REAL div(0:i1,0:j1,0:k1)
      
      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)*u(i,j,k) - ru(im)*u(im,j,k) ) / 
     1              (rp(i)*(ru(i)-ru(im)))
     2 +  ( v(i,j,k) - v(i,jm,k) )              / 
     2              (rp(i)*dphi)        
     3 +  ( w(i,j,k) - w(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 ***

c... d(1...)=drr, d(2,...)=doo, d(3,...)=dzz, 
c... d(4,..)=dro, d(5,...)=dzo, d(6,...)=drz


      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

               dij(1,i,j,k) = ( u(i,j,k)-u(im,j,k) ) / (ru(i)-ru(im
     $              )) - athird*div(i,j,k)             

               dij(2,i,j,k) = ( (v(i,j,k)-v(i,jm,k)) / dphi +(u(i
     $              ,j,k)+u(im,j,k)) / 2.0 ) / rp(i) - athird*div(i,j
     $              ,k)

               dij(3,i,j,k) = ( w(i,j,k) - w(i,j,km) ) / (zw(k)
     $              -zw(km)) - athird*div(i,j,k)


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


               dij(4,i,j,k) = 0.5 * 0.25 * (
     & ( ru(i) * ( v(ip,j,k)/rp(ip)  - v(i,j,k)/rp(i)   ) /
     &           ( Rp(ip) - Rp(i) )                             +
     &           ( u(i,jp,k)         - u(i,j,k)         ) /
     &           ( ru(i) * dphi )                               )
     &      +
     & ( ru(i) * ( v(ip,jm,k)/rp(ip) - v(i,jm,k)/rp(i)  ) /
     &           ( rp(ip) - rp(i) )                             +
     &           ( u(i,j,k)          - u(i,jm,k)        ) /
     &           ( ru(i) * dphi )                               )
     &      +
     & ( ru(im)* ( v(i,j,k)/rp(i)    - v(im,j,k)/rp(im) ) /
     &           ( rp(i) - rp(im) )                             +
     &           ( u(im,jp,k)        - u(im,j,k)        ) /
     &           ( ru(im) * dphi )                              )
     &      +
     & ( ru(im)* ( v(i,jm,k)/rp(i)   - v(im,jm,k)/rp(im)) /
     &           ( rp(i) - rp(im) )                             +
     &           ( u(im,j,k)         - u(im,jm,k)       ) /
     &           ( ru(im) * dphi )                              )
     &                       )


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


               dij(5,i,j,k) = 0.5 * 0.25 * (
     & (  ( w(i,jp,k)  - w(i,j,k)   ) / ( rp(i) * dphi )  +
     &    ( v(i,j,kp)  - v(i,j,k)   ) / ( zp(kp)-zp(k) )  )
     &      +
     & (  ( w(i,j,k)   - w(i,jm,k)  ) / ( rp(i) * dphi )  +
     &    ( v(i,jm,kp) - v(i,jm,k)  ) / ( zp(kp)-zp(k) )  )
     &      +
     & (  ( w(i,jp,km) - w(i,j,km)  ) / ( rp(i) * dphi )  +
     &    ( v(i,j,k)   - v(i,j,km)  ) / ( zp(k)-zp(km) )  )
     &      +
     & (  ( w(i,j,km)  - w(i,jm,km) ) / ( rp(i) * dphi )  +
     &    ( v(i,jm,k)  - v(i,jm,km) ) / ( zp(k)-zp(km) )  )
     &                       )


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


               dij(6,i,j,k) = 0.5 * 0.25 * (
     & (  ( w(ip,j,k)  - w(i,j,k)   ) / ( rp(ip)-rp(i) )  +
     &    ( u(i,j,kp)  - u(i,j,k)   ) / ( zp(kp)-zp(k) )  )
     &      +
     & (  ( w(ip,j,km) - w(i,j,km)  ) / ( rp(ip)-rp(i) )  +
     &    ( u(i,j,k)   - u(i,j,km)  ) / ( zp(k)-zp(km) )  )
     &      +
     & (  ( w(i,j,k)   - w(im,j,k)  ) / ( rp(i)-rp(im) )  +
     &    ( u(im,j,kp) - u(im,j,k)  ) / ( zp(kp)-zp(k) )  )
     &      +
     & (  ( w(i,j,km)  - w(im,j,km) ) / ( rp(i)-rp(im) )  +
     &    ( u(im,j,k)  - u(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

            dij(1,i,j,k) = ( u(i,j,k)-u(im,j,k) ) / (ru(i)-ru(im
     $         )) - athird*div(i,j,k)             

            dij(2,i,j,k) = ( (v(i,j,k)-v(i,jm,k)) / dphi +(u(i
     $           ,j,k)+u(im,j,k)) / 2.0 ) / rp(i) - athird*div(i,j
     $           ,k)

            dij(3,i,j,k) = ( w(i,j,k) - w(i,j,km) ) / (zw(k)
     $           -zw(km)) - athird*div(i,j,k)


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


            dij(4,i,j,k) = 0.5 * 0.25 * (
     & ( ru(i) * ( v(ip,j,k)/rp(ip)  - v(i,j,k)/rp(i)   ) /
     &           ( rp(ip) - rp(i) )                             +
     &           ( u(i,jp,k)         - u(i,j,k)         ) /
     &           ( ru(i) * dphi )                               )
     &      +
     & ( ru(i) * ( v(ip,jm,k)/rp(ip) - v(i,jm,k)/rp(i)  ) /
     &           ( rp(ip) - rp(i) )                             +
     &           ( u(i,j,k)          - u(i,jm,k)        ) /
     &           ( ru(i) * dphi )                               )
cR=0 &      +
cR=0 & ( ru(im)* ( v(i,j,k)/rp(i)    - v(im,j,k)/rp(im) ) /
cR=0 &           ( rp(i) - rp(im) )                             +
cR=0 &           ( u(im,jp,k)        - u(im,j,k)        ) /
cR=0 &           ( ru(im) * dphi )                              )
cR=0 &      +
cR=0 & ( ru(im)* ( v(i,jm,k)/rp(i)   - v(im,jm,k)/rp(im)) /
cR=0 &           ( rp(i) - rp(im) )                             +
cR=0 &           ( u(im,j,k)         - u(im,jm,k)       ) /
cR=0 &           ( ru(im) * dphi )                              )
     &                       )


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


            dij(5,i,j,k) = 0.5 * 0.25 * (
     & (  ( w(i,jp,k)  - w(i,j,k)   ) / ( rp(i) * dphi )  +
     &    ( v(i,j,kp)  - v(i,j,k)   ) / ( zp(kp)-zp(k) )  )
     &      +
     & (  ( w(i,j,k)   - w(i,jm,k)  ) / ( rp(i) * dphi )  +
     &    ( v(i,jm,kp) - v(i,jm,k)  ) / ( zp(kp)-zp(k) )  )
     &      +
     & (  ( w(i,jp,km) - w(i,j,km)  ) / ( rp(i) * dphi )  +
     &    ( v(i,j,k)   - v(i,j,km)  ) / ( zp(k)-zp(km) )  )
     &      +
     & (  ( w(i,j,km)  - w(i,jm,km) ) / ( rp(i) * dphi )  +
     &    ( v(i,jm,k)  - v(i,jm,km) ) / ( zp(k)-zp(km) )  )
     &                       )

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


            dij(6,i,j,k) = 0.5 * 0.25 * (
     & (  ( w(ip,j,k)  - w(i,j,k)   ) / ( rp(ip)-rp(i) )  +
     &    ( u(i,j,kp)  - u(i,j,k)   ) / ( zp(kp)-zp(k) )  )
     &      +
     & (  ( w(ip,j,km) - w(i,j,km)  ) / ( rp(ip)-rp(i) )  +
     &    ( u(i,j,k)   - u(i,j,km)  ) / ( zp(k)-zp(km) )  )
cR=0 &      +
cR=0 & (  ( w(i,j,k)   - w(im,j,k)  ) / ( rp(i)-rp(im) )  +
cR=0 &    ( u(im,j,kp) - u(im,j,k)  ) / ( dz )            )
cR=0 &      +
cR=0 & (  ( w(i,j,km)  - w(im,j,km) ) / ( rp(i)-rp(im) )  +
cR=0 &    ( u(im,j,k)  - u(im,j,km) ) / ( dz )            )
     &                       )

         END DO

      END DO

c... compute |d_i,j| = sqrt(2*d_(i,j)*d_(i,j)) ... 

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

               magdij(i,j,k) = SQRT( 2*dij(1,i,j,k)**2 + 2*dij(2,i,j,k)
     $              **2 +2*dij(3,i,j,k)**2 + 4*dij(4,i,j,k)**2 + 4*dij(5
     $              ,i,j,k)**2+ 4*dij(6,i,j,k)**2 )

            END DO
         END DO
      END DO

c... in the theta direction things are periodic so we can supply
c... the following boundary values for j=0 and j=j1 -- these 
c... values are used when we test filter |d^bar_ij| * d^bar_ij

      DO l=1,6
         DO i=ib,ie
            DO k=kb,ke
            
               dij(l,i,0,k)  = dij(l,i,je,k)
               dij(l,i,j1,k) = dij(l,i,1,k)

               magdij(i,0,k)  = magdij(i,je,k)
               magdij(i,j1,k) = magdij(i,1,k)

            END DO
         END DO
      END DO
      

      RETURN 
      END 
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------

      

