subroutine dns_stress(txx,txy,txz,tyy,tyz,tzz)
! using the 'sgs' sign convention for stress, so there is a - sign
use types,only:rprec
use param,only:ld,ld_big,nx,ny,nz,z_i,u_star,nu_molec
use sim_param,only:dudx,dudy,dudz,dvdx,dvdy,dvdz,dwdx,dwdy,dwdz
implicit none
real(kind=rprec),dimension(ld,ny,nz),intent(out)::txx,txy,txz,tyy, tyz,tzz
real(kind=rprec)::S11,S12,S13,S22,S23,S33
real(kind=rprec),dimension(ld,ny)::txzp,tyzp
real(kind=rprec),dimension(nz)::zz
real(kind=rprec)::delta,nu
integer::jx,jy,jz

! non-dimensional molecular viscosity
nu=nu_molec/(z_i*u_star)
! uvp-nodes at height Delta_z/2
do jy=1,ny
do jx=1,nx
   S11=dudx(jx,jy,1)
   S12=0.5_rprec*(dudy(jx,jy,1)+dvdx(jx,jy,1))
   S22=dvdy(jx,jy,1)
   S33=dwdz(jx,jy,1)
   txx(jx,jy,1)=-2.*nu*S11
   txy(jx,jy,1)=-2.*nu*S12
   tyy(jx,jy,1)=-2.*nu*S22
   tzz(jx,jy,1)=-2.*nu*S33
end do
end do

! uvp-nodes
do jz=2,nz-2
do jy=1,ny
do jx=1,nx
   S11=dudx(jx,jy,jz)
   S12=0.5_rprec*(dudy(jx,jy,jz)+dvdx(jx,jy,jz))
   S22=dvdy(jx,jy,jz)
   S33=dwdz(jx,jy,jz)
   txx(jx,jy,jz)=-2.*nu*S11
   txy(jx,jy,jz)=-2.*nu*S12
   tyy(jx,jy,jz)=-2.*nu*S22
   tzz(jx,jy,jz)=-2.*nu*S33
end do
end do
end do

! uvp-nodes near top (don't need top, since we impose the vel. there)
do jy=1,ny
do jx=1,nx
   S11=dudx(jx,jy,nz-1)
   S12=0.5_rprec*(dudy(jx,jy,nz-1)+dvdx(jx,jy,nz-1))
   S22=dvdy(jx,jy,nz-1)
   S33=dwdz(jx,jy,nz-1)
   txx(jx,jy,nz-1)=-2._rprec*nu*S11
   txy(jx,jy,nz-1)=-2._rprec*nu*S12
   tyy(jx,jy,nz-1)=-2._rprec*nu*S22
   tzz(jx,jy,nz-1)=-2._rprec*nu*S33
end do
end do

! w-nodes
! leave the wall level alone: taken care of with wall stress
do jz=2,nz-1
do jy=1,ny
do jx=1,nx
   S13=0.5_rprec*(dudz(jx,jy,jz)+dwdx(jx,jy,jz))
   S23=0.5_rprec*(dvdz(jx,jy,jz)+dwdy(jx,jy,jz))
   txz(jx,jy,jz)=-2._rprec*nu*S13
   tyz(jx,jy,jz)=-2._rprec*nu*S23
end do
end do
end do
txz(:,:,nz)=0._rprec
tyz(:,:,nz)=0._rprec
end subroutine dns_stress
