! this is the w-node version
subroutine lagrange_Ssim(S11,S12,S13,S22,S23,S33,jt)
! standard dynamic model to calculate the Smagorinsky coefficient
! this is done layer-by-layer to save memory
! everything is done to be on uv-nodes
! -note: we need to calculate |S| here, too.
! stuff is done on uv-nodes
! can save more mem if necessary. mem requirement ~ n^2, not n^3
!TSuse sgsmodule,only:F_LM,F_MM,F_QN,F_NN,u_lag,v_lag,w_lag,initest
use types,only:rprec
use param,only:ld,lh,nx,ny,nz,dx,dy,dz,dt,cs_count,use_bldg
use sim_param,only:u,v,w
use sgsmodule,only:F_LM,F_MM,Beta,initest,Cs_opt2
use test_filtermodule
!TSC
use immersedbc
implicit none
!TSC
integer::px,py,lx,ly,lz,i
integer :: jx,jy,jz,jt
real(kind=rprec), dimension(ld,ny,nz+1) :: w_nod
real(kind=rprec) :: u_rms,v_rms,uw_rms,w_rms,u_av,v_av,w_av
real(kind=rprec), dimension(nx,ny) :: u_pr,v_pr,w_pr
real(kind=rprec) :: fractus

real(kind=rprec), dimension(ld,ny,nz) :: S11,S12,S13,S22,S23,S33
real(kind=rprec), dimension(ld,ny) :: L11,L12,L13,L22,L23,L33
real(kind=rprec), dimension(ld,ny) :: M11,M12,M13,M22,M23,M33
real(kind=rprec), dimension(ld,ny) :: fourbeta

real(kind=rprec), dimension(nz) :: LMvert,MMvert

real(kind=rprec), dimension(ld,ny) :: LM,MM,Tn,epsi,dumfac

real(kind=rprec), dimension(ld,ny) :: S_bar,S11_bar,S12_bar,&
     S13_bar,S22_bar,S23_bar,S33_bar,S_S11_bar, S_S12_bar,&
     S_S13_bar, S_S22_bar, S_S23_bar, S_S33_bar

real(kind=rprec), dimension(ld,ny) :: u_bar,v_bar,w_bar
!     hat is tilda, bar, hat, three filterings
real(kind=rprec), dimension(ld,ny) :: S
real(kind=rprec) :: delta,const
real(kind=rprec) :: lagran_dt,opftdelta,powcoeff
!TSC
delta = 2._rprec*(dx*dy*dz)**(1._rprec/3._rprec)
opftdelta = 1.5_rprec*delta
powcoeff = -1._rprec/8._rprec
!TS FIX THE BUG IN COMPAQ COMPILER (USE INVERSE)
!opftdelta = 2._rprec/3_rprec/delta
!powcoeff = 1._rprec/8._rprec
!TSC
if (use_bldg) then
      do i=1,n_bldg
         px=bldg_pts(1,i)
         py=bldg_pts(2,i)
         lx=bldg_pts(3,i)
         ly=bldg_pts(4,i)
         lz=bldg_pts(5,i)
         print*,'u_bar',lx,u(px+2,py-3:py+ly+3,10)
         do jx=px+1,px+lx-1
         do jy=py+1,py+ly-1
         do jz=1,1+lz-1
         u(jx,jy,jz)=.5_rprec*(u(px,jy,jz)+&
              (jx-px)*(u(px+lx,jy,jz)-u(px,jy,jz))/lx)+&
              .5_rprec*(u(jx,py,jz)+&
              (jy-py)*(u(jx,py+ly,jz)-u(jx,py,jz))/ly)
         v(jx,jy,jz)=.5_rprec*(u(px,jy,jz)+&
              (jx-px)*(u(px+lx,jy,jz)-u(px,jy,jz))/lx)+&
              .5_rprec*(u(jx,py,jz)+&
              (jy-py)*(u(jx,py+ly,jz)-u(jx,py,jz))/ly)
         w(jx,jy,jz)=.5_rprec*(u(px,jy,jz)+&
              (jx-px)*(u(px+lx,jy,jz)-u(px,jy,jz))/lx)+&
              .5_rprec*(u(jx,py,jz)+&
              (jy-py)*(u(jx,py+ly,jz)-u(jx,py,jz))/ly)
         enddo
         enddo
         enddo
         print*,'u_aft',u(px+2,py-3:py+ly+3,10)
      end do
   end if

lagran_dt=dt*real(cs_count,kind=rprec)
call interpolag_Ssim(jt)
!! (this is the time step used in the lagrangian computations)
fractus= 1._rprec/(ny*nx)

const = 2._rprec*delta**2
do jz =1,nz-1
! w_nod is the value of w interpolated on uvp nodes
   w_nod(:,:,jz) = (w(:,:,jz)+w(:,:,jz+1))/2._rprec
end do
w_nod(:,:,nz) = w(:,:,nz-1)

Beta(:,:,1) = (1._rprec-0.65_rprec*Exp(-0.85_rprec*0.5_rprec*dz/delta))
Beta(ld,:,1) = 0._rprec
Beta(ld-1,:,1) = 0._rprec
do jz= 2,nz
   Beta(:,:,jz) = (1._rprec-0.65_rprec*Exp(-0.85_rprec*(jz-1)*dz/delta))
   Beta(ld,:,jz) = 0._rprec
   Beta(ld-1,:,jz) = 0._rprec
end do

do jz = 1,nz
! using L_ij as temp storage here
!!! watch the 0.25's:  recall w = c*z^2 close to wall, so get 0.25
   if (jz == 1) then  ! put on uvp-nodes
      L11(:,:) = u(:,:,1)*u(:,:,1)  ! uv-node
      L12(:,:) = u(:,:,1)*v(:,:,1)  ! uv-node
      L13(:,:) = u(:,:,1)*0.25_rprec*w(:,:,2)  ! parabolic interp.
      L22(:,:) = v(:,:,1)*v(:,:,1)  ! uv-node
      L23(:,:) = v(:,:,jz)*0.25_rprec*w(:,:,2)  ! uv-node
      L33(:,:) = (0.25_rprec*w(:,:,2))**2  ! uv-node
      u_bar(:,:) = u(:,:,1)
      v_bar(:,:) = v(:,:,1)
      w_bar(:,:) = 0.25_rprec*w(:,:,2)
   else  ! w-nodes
      u_bar(:,:) = 0.5_rprec*(u(:,:,jz) + u(:,:,jz-1))
      v_bar(:,:) = 0.5_rprec*(v(:,:,jz) + v(:,:,jz-1))
      w_bar(:,:) = w(:,:,jz)   
      L11(:,:)=u_bar*u_bar
      L12(:,:)=u_bar*v_bar
      L13(:,:)=u_bar*w_bar
      L23(:,:)=v_bar*w_bar
      L22(:,:)=v_bar*v_bar
      L33(:,:)=w_bar*w_bar
   end if
   call test_filter(u_bar,G_test)
   call test_filter(v_bar,G_test)
   call test_filter(w_bar,G_test)
   call test_filter(L11,G_test)  ! in-place filtering
   L11 = L11 - u_bar*u_bar
   call test_filter(L12,G_test)
   L12 = L12 - u_bar*v_bar
   call test_filter(L13,G_test)
   L13 = L13 - u_bar*w_bar
   call test_filter(L22,G_test)
   L22 = L22 - v_bar*v_bar
   call test_filter(L23,G_test)
   L23 = L23 - v_bar*w_bar
   call test_filter(L33,G_test)
   L33 = L33 - w_bar*w_bar

! would be neat to have some rotating storage or something so we don't 
!   have to do this (use pointers)
! calculate |S|
   S(:,:) = sqrt(2._rprec*(S11(:,:,jz)**2 + S22(:,:,jz)**2 + &
        S33(:,:,jz)**2 + 2._rprec*(S12(:,:,jz)**2 + &
        S13(:,:,jz)**2 + S23(:,:,jz)**2)))

! S_ij already on w-nodes
   S11_bar(:,:) = S11(:,:,jz)  
   S12_bar(:,:) = S12(:,:,jz)  
   S13_bar(:,:) = S13(:,:,jz)  
   S22_bar(:,:) = S22(:,:,jz)  
   S23_bar(:,:) = S23(:,:,jz)  
   S33_bar(:,:) = S33(:,:,jz)  

   call test_filter(S11_bar,G_test)
   call test_filter(S12_bar,G_test)
   call test_filter(S13_bar,G_test)
   call test_filter(S22_bar,G_test)
   call test_filter(S23_bar,G_test)
   call test_filter(S33_bar,G_test)

   S_bar = sqrt(2._rprec*(S11_bar**2 + S22_bar**2 + S33_bar**2 + &
        2._rprec*(S12_bar**2 + S13_bar**2 + S23_bar**2)))

   S_S11_bar(:,:) = S(:,:)*S11(:,:,jz)
   S_S12_bar(:,:) = S(:,:)*S12(:,:,jz)
   S_S13_bar(:,:) = S(:,:)*S13(:,:,jz)
   S_S22_bar(:,:) = S(:,:)*S22(:,:,jz)
   S_S23_bar(:,:) = S(:,:)*S23(:,:,jz)
   S_S33_bar(:,:) = S(:,:)*S33(:,:,jz)

   call test_filter(S_S11_bar,G_test)
   call test_filter(S_S12_bar,G_test)
   call test_filter(S_S13_bar,G_test)
   call test_filter(S_S22_bar,G_test)
   call test_filter(S_S23_bar,G_test)
   call test_filter(S_S33_bar,G_test)     

! now put beta back into M_ij
   fourbeta=4._rprec*Beta(:,:,jz)
! it might be faster just to go ahead and not store M_ij--just code directly
   M11 = const*(S_S11_bar - fourbeta*S_bar*S11_bar)
   M12 = const*(S_S12_bar - fourbeta*S_bar*S12_bar)
   M13 = const*(S_S13_bar - fourbeta*S_bar*S13_bar)
   M22 = const*(S_S22_bar - fourbeta*S_bar*S22_bar)
   M23 = const*(S_S23_bar - fourbeta*S_bar*S23_bar)
   M33 = const*(S_S33_bar - fourbeta*S_bar*S33_bar)
   LM = L11*M11+L22*M22+L33*M33+2._rprec*(L12*M12+L13*M13+L23*M23)
   MM = M11**2+M22**2+M33**2+2._rprec*(M12**2+M13**2+M23**2)

! not an efficient initialization
   if (initest .ne. 2) then
      F_MM (:,:,jz) = MM
      F_LM (:,:,jz) = .025_rprec*MM
      F_MM(ld,:,jz)=1._rprec
      F_LM(ld,:,jz)=1._rprec
      F_MM(ld-1,:,jz)=1._rprec
      F_LM(ld-1,:,jz)=1._rprec
   end if
!TS FIX THE BUG IN COMPAQ COMPILER
   Tn =max(real(F_LM(:,:,jz)*F_MM(:,:,jz)),real(1E-32))
   Tn = opftdelta*(Tn**powcoeff)
   dumfac = lagran_dt/Tn
   epsi = dumfac / (1._rprec+dumfac)
!TS	print*,'bobby'
   F_LM(:,:,jz)=(epsi*LM + (1._rprec-epsi)*F_LM(:,:,jz))
   F_MM(:,:,jz)=(epsi*MM + (1._rprec-epsi)*F_MM(:,:,jz))
   F_LM(:,:,jz)= max(real(1E-32),real(F_LM(:,:,jz)))
   Cs_opt2(:,:,jz) = F_LM(:,:,jz)/F_MM(:,:,jz)
   Cs_opt2(ld,:,jz) = 0._rprec
   Cs_opt2(ld-1,:,jz) = 0._rprec
   Cs_opt2(:,:,jz)= max(real(1E-32),real(Cs_opt2(:,:,jz)))
   if (mod(jt,1000) == 0) then
      if (jz==1) then
         do jy= 1,ny
            write(90,5168) jt*dt,float(jz),(LM(jx,jy),jx=1,nx)
            write(91,5168) jt*dt,float(jz),(MM(jx,jy),jx=1,nx)
            write(92,5168) jt*dt,float(jz),(F_LM(jx,jy,jz),jx=1,nx)
            write(93,5168) jt*dt,float(jz),(F_MM(jx,jy,jz),jx=1,nx)
            write(94,5168) jt*dt,float(jz),(LM(jx,jy)/MM(jx,jy),jx=1,nx)
         end do
      end if
      if (jz==nz/2) then
         do jy= 1,ny
            write(90,5168) jt*dt,float(jz),(LM(jx,jy),jx=1,nx)
            write(91,5168) jt*dt,float(jz),(MM(jx,jy),jx=1,nx)
            write(92,5168) jt*dt,float(jz),(F_LM(jx,jy,jz),jx=1,nx)
            write(93,5168) jt*dt,float(jz),(F_MM(jx,jy,jz),jx=1,nx)
            write(94,5168) jt*dt,float(jz),(LM(jx,jy)/MM(jx,jy),jx=1,nx)
         end do
      end if
      if (jz==nz-2) then
         do jy= 1,ny
            write(90,5168) jt*dt,float(jz),(LM(jx,jy),jx=1,nx)
            write(91,5168) jt*dt,float(jz),(MM(jx,jy),jx=1,nx)
            write(92,5168) jt*dt,float(jz),(F_LM(jx,jy,jz),jx=1,nx)
            write(93,5168) jt*dt,float(jz),(F_MM(jx,jy,jz),jx=1,nx)
            write(94,5168) jt*dt,float(jz),(LM(jx,jy)/MM(jx,jy),jx=1,nx)
         end do
      end if
   end if

   if (mod(jt,200) == 0) then
      LMvert(jz) = sum(sum(LM,DIM=1),DIM=1)/ny/nx
      MMvert(jz) = sum(sum(MM,DIM=1),DIM=1)/ny/nx
   end if
! Note: these are not the rms but the standard dev
! Square these values to get the rms
   if (mod(jt,10)==0 .AND. jt>1)then
 !	all the folowing is on UVP nodes
      u_av  = sum(sum(u(:,:,jz),DIM=1),DIM=1)*fractus
      v_av  = sum(sum(v(:,:,jz),DIM=1),DIM=1)*fractus
      w_av  = sum(sum(w_nod(:,:,jz),DIM=1),DIM=1)*fractus
      do jx=1,nx
         u_pr (jx,:) = u(jx,:,jz)-u_av
         v_pr (jx,:) = v(jx,:,jz)-v_av
         w_pr (jx,:) = w_nod(jx,:,jz)-w_av
      end do
      u_rms = sqrt(sum(sum((u_pr*u_pr),DIM=1),DIM=1)*fractus)
      v_rms = sqrt(sum(sum((v_pr*v_pr),DIM=1),DIM=1)*fractus)
      w_rms = sqrt(sum(sum((w_pr*w_pr),DIM=1),DIM=1)*fractus)
      uw_rms=- sum(sum((u_pr*w_pr),DIM=1),DIM=1)*fractus
      write(96,5200) jt*dz,u_av,u_rms,v_rms,w_rms,uw_rms
   end if
! this ends the main jz=1-nz loop          
end do
initest = 2
if (mod(jt,200) == 0) then
   write(79,5130) jt*dt,sum(sum(Cs_opt2,DIM=1),DIM=1)/ny/nx
   write(62,5130) jt*dt,sum(sum(Beta,DIM=1),DIM=1)/ny/nx
   write(65,5130) jt*dt,sum(sum(F_LM,DIM=1),DIM=1)/ny/nx
   write(66,5130) jt*dt,sum(sum(F_MM,DIM=1),DIM=1)/ny/nx
   write(67,5130) jt*dt,LMvert
   write(68,5130) jt*dt,MMvert
end if
if (mod(jt,1000) == 0) then
   do jy= 1,ny
      write(78,5168) jt*dt,1.0,(Cs_opt2(jx,jy,1),jx=1,nx)
   end do
   do jy= 1,ny
      write(78,5168) jt*dt,nz/2,(Cs_opt2(jx,jy,nz/2),jx=1,nx)
   end do
   do jy= 1,ny
      write(78,5168) jt*dt,nz-2,(Cs_opt2(jx,jy,nz-2),jx=1,nx)
   end do
end if
!CTS
print*,'CS_OPT2',Cs_opt2(30-3:30+8+3,15+2,10)
!CTS
print*,'CS_OPT2',maxval(CS_opt2)
 5168     format(1400(E14.5))
 5200     format(130(E14.5))
 5130     format(2100(E14.5))
end subroutine lagrange_Ssim
