c***********************************************************************
c     This subroutine calculates the part of the stress tensor T due to
c     surface tension. Temporary values are calculated in the center
c     of the cell shown below, then the appropriate averages are
c     performed. The equation for T is : T = - sigh*|dc|*(I - n@n) 
c        
c                            ON OUTPUT:
c
c     S11, S22, S33, S12, S13, S23:        Sij*tau/h 
c             
c***   
c
c                                       u               cc,p,S11
c             p(i-1) *------------------*------------------*
c                   /.                                    /|
c                  / .                S12                / |
c          v(i-1) *  .               *                v *  |
c                /   .                                 /   |
c               /    .       u(j-1)            p(j-1) /    |
c   p(i-1,j-1) *------------------*------------------*     |
c              |     .                 S13           |     |
c              |     + w(i-1)           +            |     * w
c              |     .                               |     |
c              |     .                               | S23 |
c     S23(i-1) |  +  .                               |  *  |
c              |     .                               |     |
c              |     .         S13(j-1)              |     |
c   w(i-1,j-1) *     .            *           w(j-1) *     |
c              |     .                               |     |
c              |     + . . . . . . . . .+ . . . . . .| . . * p(k-1) 
c              |    . p(i-1,k-1)      u(k-1)         |    /         
c              |   .                                 |   /         
c              |  + v(i-1,k-1)       +               |  * v(k-1)   
c              | .                  S12(k-1)         | /
c              |.                                    |/
c              *------------------*------------------* p(j-1,k-1)
c     p(i-1,j-1,k-1)          u(j-1,k-1)
c
c
c***********************************************************************
      SUBROUTINE stresses( cc, S11, S22, S33, S12, S13, S23,
     %                    S11t,S22t,S33t,S12t,S13t,S23t,mag)
c***
      INCLUDE 'undefined.h'
      INCLUDE 'dimension.h'
      INCLUDE 'physics.h'
      INTEGER i, j, k
      DOUBLE PRECISION   cc(nx,ny,nz), mag(nx,ny,nz)
      DOUBLE PRECISION  S11(nx,ny,nz), S22(nx,ny,nz), S33(nx,ny,nz)
      DOUBLE PRECISION  S12(nx,ny,nz), S13(nx,ny,nz), S23(nx,ny,nz)
      DOUBLE PRECISION S12t(nx,ny,nz),S13t(nx,ny,nz),S23t(nx,ny,nz)
      DOUBLE PRECISION S11t(nx,ny,nz),S22t(nx,ny,nz),S33t(nx,ny,nz)
      DOUBLE PRECISION sigh,MINGRAD,diff
      PARAMETER (MINGRAD=1.0d-50)
      INTRINSIC DSQRT,DABS
c***
      sigh = .25d0*tau**2*sigma/h**3
c***
c     calculate the components of the gradient of cc in the center
c     of the above cell. (Stored in S11,S22,S33 to save memory space)
c***
      do k=2,nz
         do j=2,ny
            do i=2,nx
               S11(i,j,k) = (cc(i,j,k) + cc(i,j-1,k) + cc(i,j,k-1)
     %                + cc(i,j-1,k-1)) - (cc(i-1,j,k) + cc(i-1,j-1,k)
     %                + cc(i-1,j,k-1) + cc (i-1,j-1,k-1))
               S22(i,j,k) = (cc(i,j,k) + cc(i-1,j,k) + cc(i,j,k-1)
     %                + cc(i-1,j,k-1)) - (cc(i,j-1,k) + cc(i-1,j-1,k)
     %                + cc(i,j-1,k-1) + cc(i-1,j-1,k-1))
               S33(i,j,k) = (cc(i,j,k) + cc(i,j-1,k) + cc(i-1,j,k)
     %                + cc(i-1,j-1,k)) - (cc(i,j,k-1) + cc(i,j-1,k-1)
     %                + cc(i-1,j,k-1) + cc(i-1,j-1,k-1))
            enddo 
         enddo 
      enddo 
c***
c     calculate the unit normal vector n and the modulus of grad cc as:
c     mag  = |dc| = MAX(|dc|,MINGRAD), to avoid division by zero.
c***
      do k=2,nz
         do j=2,ny
            do i=2,nx
               mag(i,j,k) = DSQRT(S11(i,j,k)**2 + S22(i,j,k)**2 +
     %              S33(i,j,k)**2)
               diff = MINGRAD - mag(i,j,k)
               mag(i,j,k) = mag(i,j,k) + 0.5d0*(DABS(diff) + diff)
               S11(i,j,k) = S11(i,j,k)/mag(i,j,k)
               S22(i,j,k) = S22(i,j,k)/mag(i,j,k)
               S33(i,j,k) = S33(i,j,k)/mag(i,j,k)
               mag(i,j,k) = sigh*mag(i,j,k)
            enddo 
         enddo 
      enddo 
c***
c     get temporary stress tensor coefficients: S11t, S22t, S33t,
c     S12t, S13t, S23t, from: T = - sigh*|dc|*(I - n@n) 
c***
      do k=2,nz
         do j=2,ny
            do i=2,nx
               S11t(i,j,k)  = - mag(i,j,k)*
     %              (1.d0 - S11(i,j,k)*S11(i,j,k))
               S22t(i,j,k)  = - mag(i,j,k)*
     %              (1.d0 - S22(i,j,k)*S22(i,j,k))
               S33t(i,j,k)  = - mag(i,j,k)*
     %              (1.d0 - S33(i,j,k)*S33(i,j,k))
               S12t(i,j,k) =  mag(i,j,k)*
     %                    S11(i,j,k)*S22(i,j,k)
               S13t(i,j,k) =  mag(i,j,k)*
     %                    S11(i,j,k)*S33(i,j,k)
               S23t(i,j,k) =  mag(i,j,k)*
     %                    S22(i,j,k)*S33(i,j,k)
            enddo
         enddo
      enddo
c***
c     stress tensor coefficients in the vertices of the cell
c     are given as averages of the surrounding 8 points
c***
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               S11(i,j,k) = 0.125d0*(S11t(i,j,k) + S11t(i+1,j,k)
     %              +S11t(i,j+1,k)+S11t(i+1,j+1,k)+S11t(i,j,k+1)
     %              +S11t(i+1,j,k+1)+S11t(i,j+1,k+1)+S11t(i+1,j+1,k+1))
               S22(i,j,k) = 0.125d0*(S22t(i,j,k) + S22t(i+1,j,k)
     %              +S22t(i,j+1,k)+S22t(i+1,j+1,k)+S22t(i,j,k+1)
     %              +S22t(i+1,j,k+1)+S22t(i,j+1,k+1)+S22t(i+1,j+1,k+1))
               S33(i,j,k) = 0.125d0*(S33t(i,j,k) + S33t(i+1,j,k)
     %              +S33t(i,j+1,k)+S33t(i+1,j+1,k)+S33t(i,j,k+1)
     %              +S33t(i+1,j,k+1)+S33t(i,j+1,k+1)+S33t(i+1,j+1,k+1))
            enddo
         enddo
      enddo
c***
c     stress tensor coefficients in the faces of the cell are
c     given as averages along the normal to the surface
c***
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               S12(i,j,k) =  0.5d0*(S12t(i,j,k) + S12t(i,j,k+1))
               S13(i,j,k) =  0.5d0*(S13t(i,j,k) + S13t(i,j+1,k))
               S23(i,j,k) =  0.5d0*(S23t(i,j,k) + S23t(i+1,j,k))
            enddo
         enddo
      enddo
c***
      return
      end