! put everything onto w-nodes, follow original version
subroutine sgs_stag(jt,txx,txy,txz,tyy,tyz,tzz)
use types,only:rprec
use param,only:ld,ld_big,nx,nx2,ny,ny2,nz,dx,dy,dz,z_i,u_star,&
     model,Co,nnn,vonk,cs_count,molec,nu_molec,dt,sgs
use sim_param,only: u,v,w,dudx,dudy,dudz,dvdx,dvdy,dvdz,dwdx,dwdy,dwdz
use sgsmodule,only:u_lag,v_lag,w_lag,dyn_count,Cs_opt2
use bottombc,only:zo
implicit none
integer,intent(in)::jt
real(kind=rprec),dimension(nz)::l,ziko,zz
real(kind=rprec),dimension(ld,ny,nz),intent(inout)::&
     txx,txy,txz,tyy,tyz,tzz
real(kind=rprec),dimension(ld,ny,nz)::S11,S12,S13,S22,S23,S33
real(kind=rprec),dimension(ld_big,ny2,nz)::Cs_opt2_big,&
     S11_big,S12_big,S13_big,S22_big,S23_big,S33_big,S_big,&
     txx_big,txy_big,txz_big,tyy_big,tyz_big,tzz_big
real(kind=rprec)::ux, uy, uz, vx, vy, vz, wx, wy, wz
real(kind=rprec),dimension(ld,ny) :: txzp, tyzp
real(kind=rprec) :: delta, nu, const
integer::jx,jy,jz

! save the wall level so that we make sure we don't mess with it.
! Wall level stress is computed from law of the wall
txzp(:,:)=txz(:,:,1)
tyzp(:,:)=tyz(:,:,1)

! calculate Sij on w-nodes
! calculate |S| on w-nodes

! at level z=dz/2.  (Where dudz and dvdz are on UVP nodes)
do jy=1,ny
do jx=1,nx              
   ux=dudx(jx,jy,1)  ! uvp-node
   uy=dudy(jx,jy,1)  ! uvp-node
   uz=dudz(jx,jy,1)  ! uvp-node
   vx=dvdx(jx,jy,1)  ! uvp-node
   vy=dvdy(jx,jy,1)  ! uvp-node
   vz=dvdz(jx,jy,1)  ! uvp-node 
! special case
   wx=0.5_rprec*(dwdx(jx,jy,1)+dwdx(jx,jy,2))  ! uvp-node
   wy=0.5_rprec*(dwdy(jx,jy,1)+dwdy(jx,jy,2))  ! uvp-node
   wz=dwdz(jx,jy,1)  ! uvp-node
   S11(jx,jy,1)=ux          ! uvp-node
   S12(jx,jy,1)=0.5_rprec*(uy+vx) ! uvp-node
! taken care of with wall stress routine
   S13(jx,jy,1)=0.5_rprec*(uz+wx) ! uvp
   S22(jx,jy,1)=vy          ! uvp-node
! taken care of with wall stress routine 
   S23(jx,jy,1)=0.5_rprec*(vz+wy) ! uvp
   S33(jx,jy,1)=wz          ! uvp-node
end do
end do

! calculate derivatives/strain on w-nodes
do jz=2,nz-1  
do jy=1,ny
do jx=1,nx              
   ux=0.5_rprec*(dudx(jx,jy,jz) + dudx(jx,jy,jz-1))  ! w-node
   uy=0.5_rprec*(dudy(jx,jy,jz) + dudy(jx,jy,jz-1))  ! w-node
   uz=dudz(jx,jy,jz)  ! w-node
   vx=0.5_rprec*(dvdx(jx,jy,jz) + dvdx(jx,jy,jz-1))  ! w-node
   vy=0.5_rprec*(dvdy(jx,jy,jz) + dvdy(jx,jy,jz-1))  ! w-node
   vz=dvdz(jx,jy,jz)  ! w-node
   wx=dwdx(jx,jy,jz)  ! w-node
   wy=dwdy(jx,jy,jz)  ! w-node
   wz=0.5_rprec*(dwdz(jx,jy,jz) + dwdz(jx,jy,jz-1))  ! w-node
   S11(jx,jy,jz)=ux          ! w-node
   S12(jx,jy,jz)=0.5_rprec*(uy+vx) ! w-node
   S13(jx,jy,jz)=0.5_rprec*(uz+wx) ! w-node
   S22(jx,jy,jz)=vy          ! w-node
   S23(jx,jy,jz)=0.5_rprec*(vz+wy) ! w-node
   S33(jx,jy,jz)=wz          ! w-node
end do
end do
end do

! at top
do jy=1,ny
do jx=1,nx              
   ux=dudx(jx,jy,nz) ! uvp-node
   uy=dudy(jx,jy,nz) ! uvp-node
   uz=0.5_rprec*(dudz(jx,jy,nz) + dudz(jx,jy,nz-1)) ! uvp-node
   vx=dvdx(jx,jy,nz) ! uvp-node
   vy=dvdy(jx,jy,nz) ! uvp-node
   vz=0.5_rprec*(dvdz(jx,jy,nz) + dvdz(jx,jy,nz-1))! uvp-node
   wx=0.5_rprec*(dwdx(jx,jy,nz) + dwdx(jx,jy,nz-1)) ! uvp-node
   wy=0.5_rprec*(dwdy(jx,jy,nz) + dwdy(jx,jy,nz-1))! uvp-node
   wz=dwdz(jx,jy,nz) ! uvp-node
! note: there's probably something wrong with the stuff on the w-nodes here
! only very small, though
   S11(jx,jy,nz)=ux          ! uvp-node
   S12(jx,jy,nz)=0.5_rprec*(uy+vx) ! uvp-node
   S13(jx,jy,nz)=0.5_rprec*(uz+wx) ! uvp-node
   S22(jx,jy,nz)=vy          ! uvp-node
   S23(jx,jy,nz)=0.5_rprec*(vz+wy) ! uvp-node
   S33(jx,jy,nz)=wz          ! uvp-node
end do
end do

!          do jz=1,nz
!          print*,jz,sum(S12(:,:,jz)/nx/ny),sum(S13(:,:,jz)/nx/ny)
!          end do

!	This small part computes the average velocity during cs_count times steps
!	This is used with the lagrangian model only
if (model == 4 .OR. model==5) then
   u_lag = u_lag+u
   v_lag = v_lag+v
   w_lag = w_lag+w
end if


if (sgs)then!ref01
! For traditional Smagorinsky
   if(model == 1) then!ref02
      delta=(dx*dy*dz)**(1._rprec/3._rprec) ! nondimensional
      do jz=1,nz!ref03
! Define parameters (Co and nn) for wallfunction
         Cs_opt2(:,:,jz)=Co**2  ! constant coefficient
         if (jz == 1) then!ref04
            zz(jz)=(real(jz)-0.5_rprec)*dz
         else  ! w-nodes
            zz(jz)=(jz-1)*dz
         end if! end ref04
! z's nondimensional, l here is on w-nodes, except at top & bottom
       l(jz)=(Co**(nnn)*(vonk*zz(jz))**(-nnn)+(delta)**(-nnn))**(-1._rprec/nnn)
      end do!end ref03
   else ! for dynamic procedures  !cont ref02
      l=(dx*dy*dz)**(1._rprec/3._rprec)  ! l is constant
      if ((jt == 1) .and. (Cs_opt2(1,1,1) < 1.e-10)) then!ref05
         dyn_count = cs_count
         print *,'CS_opt2 initialiazed'  
!     initialize Cs_opt2 data with curve fitting of prev. results
!	    do jz=1,nz
! may wanna change this

         Cs_opt2(:,:,jz) = (0.019_rprec*log((jz-0.5_rprec)/nz) + 0.11_rprec)**2
!	Cs_opt2(:,:,:) = 0.03 

!	    end do
! make sure that the "node conventions" in these dynamic models
! match with those done here
      else if((model==2.OR.model==3).AND.mod(jt,dyn_count)==0) then!cont ref05
         if (model == 2) then  ! Standard dynamic model !ref06
            ziko = Cs_opt2(5,5,:)
            call std_dynamic(ziko,S11,S12,S13,S22,S23,S33)
            do jx=1,ld
            do jy =1,ny
               Cs_opt2(jx,jy,:)=ziko(:)
            end do
            end do
         end if!end ref06
         if (model == 3) then  ! Scale dependent dynamic model
            ziko = Cs_opt2(5,5,:)
            call scaledep_dynamic(ziko,S11,S12,S13,S22,S23,S33,jt)
            do jx=1,ld
            do jy =1,ny
               Cs_opt2(jx,jy,:)=ziko(:)
            end do
            end do
         end if
!	else if (model==4.AND.(mod(jt_count,dyn_count)==0)) then ! Lagrangian
      else if (model==4.AND.(mod(jt,cs_count)==0)) then ! LagSS	cont ref5
         call lagrange_Ssim(S11,S12,S13,S22,S23,S33,jt)
      else if (model==5.AND.(mod(jt,cs_count)==0)) then ! LagSD	cont ref5
         call lagrange_Sdep(S11,S12,S13,S22,S23,S33,jt)
      end if  ! for if jt == 1 and 2 else ifs		end ref5
   end if !(for model =1 or else)		end ref02
end if ! for sgs-if		end ref01
!	jt_count = jt_count+1

call dealias1(S11,S11_big) ! w
call dealias1(S12,S12_big) ! w
call dealias1(S13,S13_big) ! w
call dealias1(S22,S22_big) ! w
call dealias1(S23,S23_big) ! w
call dealias1(S33,S33_big) ! w
call dealias1(Cs_opt2,Cs_opt2_big) ! w

! define |S| on w-nodes
do jz=1,nz
do jy=1,ny2
do jx=1,nx2
   S_big(jx,jy,jz)=sqrt(2._rprec * (S11_big(jx,jy,jz)**2 + &
        S22_big(jx,jy,jz)**2+S33_big(jx,jy,jz)**2&
        + 2._rprec * (S12_big(jx,jy,jz)**2 + S13_big(jx,jy,jz)**2 +&
        S23_big(jx,jy,jz)**2)))
end do
end do
end do

! evaluate nu_t= c_s^2 l^2 |S|
! bottom
      
do jy=1,ny2
do jx=1,nx2
   nu = 0.
   const = 0.
   if (sgs) then
      const = Cs_opt2_big(jx,jy,1)*l(1)**2
   end if
   if (molec) then
      nu = (nu_molec/(u_star*z_i)) ! dim/less 
   end if
! everything on right node here
   txx_big(jx,jy,1) = -2._rprec*(const*S_big(jx,jy,1) + nu)*S11_big(jx,jy,1)
   txy_big(jx,jy,1) = -2._rprec*(const*S_big(jx,jy,1) + nu)*S12_big(jx,jy,1)
   tyy_big(jx,jy,1) = -2._rprec*(const*S_big(jx,jy,1) + nu)*S22_big(jx,jy,1)
   tzz_big(jx,jy,1) = -2._rprec*(const*S_big(jx,jy,1) + nu)*S33_big(jx,jy,1)
   
   if ((const*dt/(dx**2)) .GT. 0.5_rprec) then
      print*, (const*dt/(dx**2))
   end if

end do
end do
! middle
do jz=2,nz-2  ! following original code
do jy=1,ny2
do jx=1,nx2
   nu=0._rprec
   const=0._rprec
   if (sgs) then 
      const=0.5_rprec*(Cs_opt2_big(jx,jy,jz) + Cs_opt2_big(jx,jy,jz+1))*&
           0.5_rprec*(l(jz)**2 + l(jz+1)**2)
   end if
   if (molec) then
      nu=(nu_molec/(u_star*z_i)) ! dim/less
   end if
   txx_big(jx,jy,jz)=&
        -2._rprec*(const*0.5_rprec*(S_big(jx,jy,jz)*S11_big(jx,jy,jz) +& 
        S_big(jx,jy,jz+1)*S11_big(jx,jy,jz+1))+&
        nu*0.5_rprec*(S11_big(jx,jy,jz) + S11_big(jx,jy,jz+1)))
   txy_big(jx,jy,jz)=&
        -2._rprec*(const*0.5_rprec*(S_big(jx,jy,jz)*S12_big(jx,jy,jz) +& 
        S_big(jx,jy,jz+1)*S12_big(jx,jy,jz+1)) +&
        nu*0.5_rprec*(S12_big(jx,jy,jz) + S12_big(jx,jy,jz+1)))
   tyy_big(jx,jy,jz)=&
        -2._rprec*(const*0.5_rprec*(S_big(jx,jy,jz)*S22_big(jx,jy,jz) +& 
        S_big(jx,jy,jz+1)*S22_big(jx,jy,jz+1)) +&
        nu*0.5_rprec*(S22_big(jx,jy,jz) + S22_big(jx,jy,jz+1)))
   tzz_big(jx,jy,jz)=&
        -2._rprec*(const*0.5_rprec*(S_big(jx,jy,jz)*S33_big(jx,jy,jz) +& 
        S_big(jx,jy,jz+1)*S33_big(jx,jy,jz+1)) +&
        nu*0.5_rprec*(S33_big(jx,jy,jz) + S33_big(jx,jy,jz+1)))

   if ((const*dt/(dx**2)) .GT. 0.5_rprec) then
      print*, (const*dt/(dx**2))
   end if

end do
end do
end do

! like original code:
! this is slightly messed up, the levels nz, nz-1 weird
! watch the nz and nz-1 here!
do jy=1,ny2
do jx=1,nx2
   nu=0._rprec
   const=0._rprec
   if (sgs) const=Cs_opt2_big(jx,jy,nz)*l(nz)**2
   if (molec) nu=(nu_molec/(u_star*z_i)) ! make dimensionless
   txx_big(jx,jy,nz-1)=-2._rprec*(const*S_big(jx,jy,nz) + nu)*S11_big(jx,jy,nz)
   txy_big(jx,jy,nz-1)=-2._rprec*(const*S_big(jx,jy,nz) + nu)*S12_big(jx,jy,nz)
   tyy_big(jx,jy,nz-1)=-2._rprec*(const*S_big(jx,jy,nz) + nu)*S22_big(jx,jy,nz)
   tzz_big(jx,jy,nz-1)=-2._rprec*(const*S_big(jx,jy,nz) + nu)*S33_big(jx,jy,nz)

   if ((const*dt/(dx**2)) .GT. 0.5_rprec) then
      print*, (const*dt/(dx**2))
   end if

end do
end do

do jz=2,nz-1
do jy=1,ny2
do jx=1,nx2
   nu = 0._rprec
   const = 0._rprec
   if (sgs) const=Cs_opt2_big(jx,jy,jz)*l(jz)**2
   if (molec) nu=(nu_molec/(u_star*z_i)) ! dim/less
   txz_big(jx,jy,jz)=-2._rprec*(const*S_big(jx,jy,jz) + nu)*&
        S13_big(jx,jy,jz)
   tyz_big(jx,jy,jz)=-2._rprec*(const*S_big(jx,jy,jz) + nu)*&
        S23_big(jx,jy,jz)
end do
end do
end do
! what happens at the very top nodes? for txx, tyy, etc. 
! u(k)=u(k-1), etc? make unimportant?

call dealias2(txx,txx_big)
call dealias2(txy,txy_big)
call dealias2(txz,txz_big)
call dealias2(tyy,tyy_big)
call dealias2(tyz,tyz_big)
call dealias2(tzz,tzz_big)
!TSCcall dealias2(Cs_opt2,Cs_opt2_big)

! don't need to do this, unless we accidentally mess with the wall level
txz(:,:,1)=txzp(:,:)
tyz(:,:,1)=tyzp(:,:)

txz(:,:,nz)=0._rprec
tyz(:,:,nz)=0._rprec

if (mod(jt,10)==0) then
   write(97,5400) (SUM(SUM(txz(:,:,:),DIM=1),DIM=1)/nx/ny)
   print*,SUM(txz(:,:,1))
   print*,SUM(tyz(:,:,1))
end if

5400 format(2050(E14.5))
end subroutine sgs_stag
