c
c   On input:
c   
c        cc          volume fractions
c        S11         arbitrary nx*ny array;
c        S22              "
c        S12              "
c        S11t, S22t   work arrays of size nx*ny.
c        sigma     surface tension coefficient.
c        
c   
c   On output
c   
c        S11, S22, S12 stresses Sij*tau/h 
c        in momentum.f  at locations described on figure.
c   
      SUBROUTINE stresses(S11,S22,S12,cc,S11t,S22t,sigma,tau,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny
      DOUBLE PRECISION S11(nx,ny),S22(nx,ny),S12(nx,ny),cc(nx,ny)
      DOUBLE PRECISION S11t(nx,ny),S22t(nx,ny)
      DOUBLE PRECISION sigma,tau
      DOUBLE PRECISION sigmag25,h,sigmag,sigh,mag,dcx,dcy
      INTEGER i,j
c
      h = 1.d0/(nx-2)
      sigh =.5d0*tau**2*sigma/h**3
      do j=2,ny
         do i=2,nx
            dcx= cc(i,j) + cc(i,j-1) - cc(i-1,j) - cc(i-1,j-1)
            dcy= cc(i,j) + cc(i-1,j) - cc(i,j-1) - cc(i-1,j-1)
            mag = dcx**2 + dcy**2 + 1.d-50
            sigmag = sigh/dsqrt(mag)
            sigmag25 = 0.25d0*sigmag
            S11t(i,j) = - sigmag25*dcy*dcy
            S22t(i,j) = - sigmag25*dcx*dcx
            S12(i,j) =  sigmag*dcx*dcy
         enddo
      enddo
      do j=2,ny-1
         do i=2,nx-1
            S11(i,j) = S11t(i,j)+S11t(i+1,j)+S11t(i,j+1)+S11t(i+1,j+1)
            S22(i,j) = S22t(i,j)+S22t(i+1,j)+S22t(i,j+1)+S22t(i+1,j+1)
         enddo
      enddo
      return
      end