! this is the w-node version
subroutine lagrange_Sdep(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
use types,only:rprec
use param,only:ld,lh,nx,ny,nz,dx,dy,dz,dt,cs_count
use sim_param,only:u,v,w
use sgsmodule,only:F_LM,F_MM,F_QN,F_NN,beta,initest,dyn_count,Cs_opt2
use test_filtermodule
implicit none
integer :: jx,jy,jz,jt
integer :: countmc_s(nz),countmc_b(nz)
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,nz) :: Cs_opt2_2delta
real(kind=rprec), dimension(ld,ny) :: L11,L12,L13,L22,L23,L33
real(kind=rprec), dimension(ld,ny) :: Q11,Q12,Q13,Q22,Q23,Q33
real(kind=rprec), dimension(ld,ny) :: M11,M12,M13,M22,M23,M33
real(kind=rprec), dimension(ld,ny) :: N11,N12,N13,N22,N23,N33
real(kind=rprec), dimension(ld,ny) :: N11_,N12_,N13_,N22_,N23_,N33_
real(kind=rprec), dimension(ld,ny) :: L11_,L12_,L13_,L22_,L23_,L33_
real(kind=rprec), dimension(ld,ny) :: Q11_,Q12_,Q13_,Q22_,Q23_,Q33_
real(kind=rprec), dimension(ld,ny) :: M11_,M12_,M13_,M22_,M23_,M33_

real(kind=rprec), dimension(nz) :: LMvert,MMvert,QNvert,NNvert
      
real(kind=rprec), dimension(ld,ny) :: LM,MM,QN,NN,Tn,epsi,dumfac
real(kind=rprec), dimension(ld,ny) :: fourbeta

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) :: S_hat,S11_hat,S12_hat,&
     S13_hat,S22_hat,S23_hat,S33_hat,S_S11_hat, S_S12_hat,&
     S_S13_hat, S_S22_hat, S_S23_hat, S_S33_hat
real(kind=rprec), dimension(ld,ny) :: u_bar,v_bar,w_bar
real(kind=rprec), dimension(ld,ny) :: u_hat,v_hat,w_hat
!     hat is tilda, bar, hat, three filterings
real(kind=rprec), dimension(ld,ny) :: S
real(kind=rprec) :: delta,const,const2
real(kind=rprec) :: lagran_dt,opftdelta,powcoeff

if (initest .ne. 2) then
   print*, 'Beta initialized'
   Beta (:,:,:)  = 1
end if

delta = filter_size*(dx*dy*dz)**(1._rprec/3._rprec)
opftdelta = 1.5_rprec*delta
powcoeff = -1.0_rprec/8.0_rprec


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

do jz =1,nz
   w_nod(:,:,jz) = (w(:,:,jz)+w(:,:,jz+1))/2._rprec
end do
w_nod(:,:,nz+1) = w(:,:,nz)
const = 2._rprec*delta**2
const2 = 2._rprec*(2._rprec*delta)**2

do jz = nz,1,-1
! 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
      L11(:,:) = 0.5_rprec*(u(:,:,jz) + u(:,:,jz-1))*&
           0.5_rprec*(u(:,:,jz) + u(:,:,jz-1))
      L12(:,:) = 0.5_rprec*(u(:,:,jz) + u(:,:,jz-1))*&
           0.5_rprec*(v(:,:,jz) + v(:,:,jz-1))
      L13(:,:) = 0.5_rprec*(u(:,:,jz) + u(:,:,jz-1))*w(:,:,jz)
      L22(:,:) = 0.5_rprec*(v(:,:,jz) + v(:,:,jz-1))*&
           0.5_rprec*(v(:,:,jz) + v(:,:,jz-1))
      L23(:,:) = 0.5_rprec*(v(:,:,jz) + v(:,:,jz-1))*w(:,:,jz)
      L33(:,:) = w(:,:,jz)*w(:,:,jz)
      u_bar(:,:) = 0.5_rprec*(u(:,:,jz) + u(:,:,jz-1))
      v_bar(:,:) = 0.5_rprec*(v(:,:,jz) + v(:,:,jz-1))
      w_bar(:,:) = w(:,:,jz)
   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)
   u_hat = u_bar
   v_hat = v_bar
   w_hat = w_bar
   Q11 = u_bar*u_bar
   Q12 = u_bar*v_bar
   Q13 = u_bar*w_bar
   Q22 = v_bar*v_bar
   Q23 = v_bar*w_bar
   Q33 = w_bar*w_bar

   call test_filter(u_hat,G_test_test)
   call test_filter(v_hat,G_test_test)
   call test_filter(w_hat,G_test_test)
   call test_filter(Q11,G_test_test)
   Q11 = Q11 - u_hat*u_hat
   call test_filter(Q12,G_test_test)
   Q12 = Q12 - u_hat*v_hat
   call test_filter(Q13,G_test_test)
   Q13 = Q13 - u_hat*w_hat
   call test_filter(Q22,G_test_test)
   Q22 = Q22 - v_hat*v_hat
   call test_filter(Q23,G_test_test)
   Q23 = Q23 - v_hat*w_hat
   call test_filter(Q33,G_test_test)
   Q33 = Q33 - w_hat*w_hat

! 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)
   
   
   S11_hat = S11_bar
   S12_hat = S12_bar
   S13_hat = S13_bar
   S22_hat = S22_bar
   S23_hat = S23_bar
   S33_hat = S33_bar
   
   call test_filter(S11_hat,G_test_test)
   call test_filter(S12_hat,G_test_test)
   call test_filter(S13_hat,G_test_test)
   call test_filter(S22_hat,G_test_test)
   call test_filter(S23_hat,G_test_test)
   call test_filter(S33_hat,G_test_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_hat = sqrt(2._rprec*(S11_hat**2 + S22_hat**2 + S33_hat**2 + &
        2._rprec*(S12_hat**2 + S13_hat**2 + S23_hat**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)     

   S_S11_hat(:,:) = S_bar(:,:)*S11_bar(:,:)
   S_S12_hat(:,:) = S_bar(:,:)*S12_bar(:,:)
   S_S13_hat(:,:) = S_bar(:,:)*S13_bar(:,:)
   S_S22_hat(:,:) = S_bar(:,:)*S22_bar(:,:)
   S_S23_hat(:,:) = S_bar(:,:)*S23_bar(:,:)
   S_S33_hat(:,:) = S_bar(:,:)*S33_bar(:,:)
   
   call test_filter(S_S11_hat,G_test_test)
   call test_filter(S_S12_hat,G_test_test)
   call test_filter(S_S13_hat,G_test_test)
   call test_filter(S_S22_hat,G_test_test)
   call test_filter(S_S23_hat,G_test_test)
   call test_filter(S_S33_hat,G_test_test)  

! now put beta back into M_ij
   fourbeta=4._rprec*Beta(:,:,jz)
!	print*, fourbeta

! 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)
  
   N11 = const2*(S_S11_hat - fourbeta*S_hat*S11_hat)
   N12 = const2*(S_S12_hat - fourbeta*S_hat*S12_hat)
   N13 = const2*(S_S13_hat - fourbeta*S_hat*S13_hat)
   N22 = const2*(S_S22_hat - fourbeta*S_hat*S22_hat)
   N23 = const2*(S_S23_hat - fourbeta*S_hat*S23_hat)
   N33 = const2*(S_S33_hat - fourbeta*S_hat*S33_hat)

   if(jz==3) then
      N11_ = N11
      N12_ = N12
      N13_ = N13
      N23_ = N23
      N22_ = N22
      N33_ = N33
      L11_ = L11
      L12_ = L12
      L13_ = L13
      L23_ = L23
      L22_ = L22
      L33_ = L33
      Q11_ = Q11
      Q12_ = Q12
      Q13_ = Q13
      Q23_ = Q23
      Q22_ = Q22
      Q33_ = Q33
      M11_ = M11
      M12_ = M12
      M13_ = M13
      M23_ = M23
      M22_ = M22
      M33_ = M33
   else if ((jz==1) .OR. (jz==2)) then
      N11 = N11_
      N12 = N12_
      N13 = N13_
      N23 = N23_
      N22 = N22_
      N33 = N33_
      M11 = M11_
      M12 = M12_
      M13 = M13_
      M23 = M23_
      M22 = M22_
      M33 = M33_
      L11 = L11_
      L12 = L12_
      L13 = L13_
      L23 = L23_
      L22 = L22_
      L33 = L33_
      Q11 = Q11_
      Q12 = Q12_
      Q13 = Q13_
      Q23 = Q23_
      Q22 = Q22_
      Q33 = Q33_
   end if

   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)
   QN = Q11*N11 + Q22*N22 + Q33*N33 + 2._rprec*(Q12*N12 + Q13*N13 + Q23*N23)
   NN = N11**2 + N22**2 + N33**2 + 2._rprec*(N12**2 + N13**2 + N23**2)

! not an efficient initialization
   if (initest .ne. 2) then
!          print *,'F_MM and F_LM initialized' 
      F_MM (:,:,jz) = MM
      F_LM (:,:,jz) = .025_rprec*MM
      F_NN (:,:,jz) = NN
      F_QN (:,:,jz) = .025_rprec*NN
   end if

   Tn=opftdelta*((F_LM(:,:,jz)*F_MM(:,:,jz))**powcoeff)
   dumfac = lagran_dt/Tn
   epsi = dumfac / (1._rprec+dumfac)

   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

   Tn=opftdelta*2._rprec*((F_QN(:,:,jz)*F_NN(:,:,jz))**powcoeff)
   dumfac = lagran_dt/Tn
   epsi = dumfac / (1._rprec+dumfac)

   F_QN(:,:,jz)=(epsi*QN + (1._rprec-epsi)*F_QN(:,:,jz))
   F_NN(:,:,jz)=(epsi*NN + (1._rprec-epsi)*F_NN(:,:,jz))
   F_QN(:,:,jz)= max(real(1E-32),real(F_QN(:,:,jz)))
   Cs_opt2_2delta(:,:,jz) = F_QN(:,:,jz)/F_NN(:,:,jz)
   Cs_opt2_2delta(ld,:,jz) = 0._rprec
   Cs_opt2_2delta(ld-1,:,jz) = 0._rprec
   
   Beta(:,:,jz) = (Sum(F_QN(:,:,jz))/Sum(F_NN(:,:,jz)))/&
                      (Sum(F_LM(:,:,jz))/Sum(F_MM(:,:,jz)))
!	end if
   Beta(ld,:,jz) = 0._rprec
   Beta(ld-1,:,jz) = 0._rprec
  
   if (mod(jt,100) == 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)
            write(70,5168) jt*dt,float(jz),(QN(jx,jy),jx=1,nx)
            write(71,5168) jt*dt,float(jz),(NN(jx,jy),jx=1,nx)
            write(72,5168) jt*dt,float(jz),(F_QN(jx,jy,jz),jx=1,nx)
            write(73,5168) jt*dt,float(jz),(F_NN(jx,jy,jz),jx=1,nx)
            write(74,5168) jt*dt,float(jz),(QN(jx,jy)/NN(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)
            write(70,5168) jt*dt,float(jz),(QN(jx,jy),jx=1,nx)
            write(71,5168) jt*dt,float(jz),(NN(jx,jy),jx=1,nx)
            write(72,5168) jt*dt,float(jz),(F_QN(jx,jy,jz),jx=1,nx)
            write(73,5168) jt*dt,float(jz),(F_NN(jx,jy,jz),jx=1,nx)
            write(74,5168) jt*dt,float(jz),(QN(jx,jy)/NN(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)
            write(70,5168) jt*dt,float(jz),(QN(jx,jy),jx=1,nx)
            write(71,5168) jt*dt,float(jz),(NN(jx,jy),jx=1,nx)
            write(72,5168) jt*dt,float(jz),(F_QN(jx,jy,jz),jx=1,nx)
            write(73,5168) jt*dt,float(jz),(F_NN(jx,jy,jz),jx=1,nx)
            write(74,5168) jt*dt,float(jz),(QN(jx,jy)/NN(jx,jy),jx=1,nx)
         end do
      end if
   end if

   if (mod(jt,50) == 0) then
      LMvert(jz) = sum(sum(LM,DIM=1),DIM=1)/ny/nx
      MMvert(jz) = sum(sum(MM,DIM=1),DIM=1)/ny/nx
      QNvert(jz) = sum(sum(QN,DIM=1),DIM=1)/ny/nx
      NNvert(jz) = sum(sum(NN,DIM=1),DIM=1)/ny/nx
   end if
   
   if (mod(jt,10)==0 .AND. jt>1) then
      
      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 = sum(sum((u_pr*u_pr),DIM=1),DIM=1)*fractus
      v_rms = sum(sum((v_pr*v_pr),DIM=1),DIM=1)*fractus
      w_rms = 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
5168 format(1400(E14.5))
5200 format(130(E14.5))
   
   write(99,5300) jt*dz,u(4,16,16),v(16,16,16),w(16,16,16)

5300 format(82(E14.5))

! this ends the main jz=1-nz loop          
end do
Cs_opt2(:,:,2)=Cs_opt2(:,:,3)/2._rprec
Cs_opt2(:,:,1)=Cs_opt2(:,:,3)/4._rprec
Cs_opt2_2delta(:,:,2)=Cs_opt2_2delta(:,:,3)/2._rprec
Cs_opt2_2delta(:,:,1)=Cs_opt2_2delta(:,:,3)/4._rprec
!	countmc_b(:)=0
!	countmc_s(:)=0
!	do jx=1,nx
!	do jy=1,ny
!	do jz=1,nz
!	if (Beta(jx,jy,jz) .LT. (0.5d0*dfloat(jz)/dfloat(nz))) then
!	Beta(jx,jy,jz) = (0.5d0*dfloat(jz)/dfloat(nz))
!	countmc_s(jz)=countmc_s(jz)+1
!	if (Beta(jx,jy,jz) .GT. 1.2d0) then
!	Beta(jx,jy,jz) = 1.2d0
!	countmc_b(jz)=countmc_b(jz)+1
!	end if
!	end do
!	end do
!	end do

!	print*, sum(Beta(:,:,:)/nx/ny/nz)

!      Cs_opt2(:,:,nz) = Cs_opt2(:,:,nz-1)

initest = 2
if (mod(jt,50) == 0) then
   write(79,5130) jt*dt,sum(sum(Cs_opt2,DIM=1),DIM=1)/ny/nx
   write(61,5130) jt*dt,sum(sum(Cs_opt2_2delta,DIM=1),DIM=1)/ny/nx
   write(63,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
   
   write(45,5130) jt*dt,sum(sum(F_QN,DIM=1),DIM=1)/ny/nx
   write(46,5130) jt*dt,sum(sum(F_NN,DIM=1),DIM=1)/ny/nx
   write(47,5130) jt*dt,QNvert
   write(48,5130) jt*dt,NNvert
end if
!	 if (mod(jt,5) == 0 .AND. jt<51) then
!          do jy= 1,ny
!          write(98,5168) jt*dt,2,(Cs_opt2(jx,jy,2),jx=1,nx)
!         end do
!          do jy= 1,ny
!          write(98,5168) jt*dt,12,(Cs_opt2(jx,jy,12),jx=1,nx)
!          end do
!          do jy= 1,ny
!          write(98,5168) jt*dt,22,(Cs_opt2(jx,jy,22),jx=1,nx)
!          end do
!          end if
!	 if (mod(jt,50) == 0 .AND. (jt>51 .AND. jt<1000)) then
!         do jy= 1,ny
!          write(98,5168) jt*dt,2,(Cs_opt2(jx,jy,2),jx=1,nx)
!         end do
!          do jy= 1,ny
!          write(98,5168) jt*dt,12,(Cs_opt2(jx,jy,12),jx=1,nx)
!          end do
!          do jy= 1,ny
!          write(98,5168) jt*dt,22,(Cs_opt2(jx,jy,22),jx=1,nx)
!          end do
!          end if
if (mod(jt,100) == 0) then
   
   do jy= 1,ny
      write(78,5168) jt*dt,1.0,(Cs_opt2(jx,jy,1),jx=1,nx)
      write(60,5168) jt*dt,1.0,(Cs_opt2_2delta(jx,jy,1),jx=1,nx)
      write(62,5168) jt*dt,1.0,(Beta(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)
      write(60,5168) jt*dt,nz/2,(Cs_opt2_2delta(jx,jy,nz/2),jx=1,nx)
      write(62,5168) jt*dt,nz/2,(Beta(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)
      write(60,5168) jt*dt,nz-2,(Cs_opt2_2delta(jx,jy,nz-2),jx=1,nx)
      write(62,5168) jt*dt,nz-2,(Beta(jx,jy,nz-2),jx=1,nx)
   end do

end if
!        do jz = 1,nz
!	Cs_opt2(:,:,jz)=Cs_opt2_tempo(jz)
!	end do
5130 format(2100(E14.5))
end subroutine lagrange_Sdep
