module io
use types,only:rprec
use param,only:ld,nx,ny,nz,jx_s
implicit none
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer,parameter::base=100000,nwrite=1000
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer,parameter::jx_pls=jx_s+4,jx_ple=nx/16*9,width=int(ny/5)
integer,parameter::jy_pls=ny/2-width,jy_ple=ny/2+width+1
real(kind=rprec),dimension(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)::&
     mean_u,mean_v,mean_u2,mean_v2,mean_dudx
real(kind=rprec),dimension(nx,ny,nz)::lam2
contains

subroutine openfiles()
use sim_param,only:path
! note: this ain't in the 'output' subdir
open(64,file=path//'output/velspec.out')
open(78,file=path//'output/cs2_lag_planar.out')
open(79,file=path//'output/cs2_lag.out')
open(60,file=path//'output/cs2_2Del_planar.out')
open(61,file=path//'output/cs2_2Del.out')
open(62,file=path//'output/Beta_planar.out')
open(63,file=path//'output/Beta.out')

open(65,file=path//'output/FLM_vertprof.out')
open(66,file=path//'output/FMM_vertprof.out')
open(67,file=path//'output/LM_vertprof.out')
open(68,file=path//'output/MM_vertprof.out')
open(69,file=path//'output/nondim_u_vertprof.out')

open(45,file=path//'output/FQN_vertprof.out')
open(46,file=path//'output/FNN_vertprof.out')
open(47,file=path//'output/QN_vertprof.out')
open(48,file=path//'output/NN_vertprof.out')
 
open(80,file=path//'output/explode.out')

open(81,file=path//'output/clip.out')

open(90,file=path//'output/LMplan.out')
open(91,file=path//'output/MMplan.out')
open(92,file=path//'output/FLMplan.out')
open(93,file=path//'output/FMMplan.out')
open(94,file=path//'output/LMoverMM.out')

open(70,file=path//'output/QNplan.out')
open(71,file=path//'output/NNplan.out')
open(72,file=path//'output/FQNplan.out')
open(73,file=path//'output/FNNplan.out')
open(74,file=path//'output/QNoverNN.out')
 
open(98,file=path//'output/cs2plan-smallt.out')
open(96,file=path//'output/rmsprofiles.out')
open(97,file=path//'output/tawprofiles.out')

open(95,file=path//'output/mean_u.out',form='unformatted')
write(95)jx_pls,jx_ple,jy_pls,jy_ple
open(195,file=path//'output/lambda2.out',form='unformatted')
write(195)nx,ny,nz
end subroutine openfiles

!subroutine initial
!end subroutine initial

subroutine output_loop(jt,txx,txy,txz,tyy,tyz,tzz)
use param,only:output,dt,c_count
use sim_param,only:path,u,v,w,dudz,dudx,p,&
     RHSx,RHSy,RHSz
use sgsmodule,only:Cs_opt2,Cs_opt2_avg
implicit none
integer,intent(in)::jt
real(kind=rprec),dimension(ld,ny,nz)::txx,txy,txz,tyy,tyz,tzz
real(kind=rprec),dimension(nz)::u_ndim
character(len=20)::req
character(len=24)::fname
if (mod(jt,100) == 0) then
   u_ndim=sum(sum(dudz,DIM=1),DIM=1)/ny/nx
   write(69,5130) real(jt*dt),real(u_ndim)
end if
5130   format(2100(E14.5))
if (output.and.(mod(jt,base) == 0)) then
! checkpoint
!TS   open(1,file=path//'request')
!TS   read(1,*) req
!TS   close(1)
!TS   if (req.eq.'save') then
!TS      rewind (21)
!TS      write(21) u,v,w,RHSx,RHSy,RHSz,Cs_opt2
!TS   end if
! write out vel.out along the way
   if (mod(jt,base) == 0) then
      write(fname,'(A12,i5.5,A4)')path//'output/vel',jt,'.out'
      open(1,file=fname,form='unformatted')
      write(1) u,v,w,RHSx,RHSy,RHSz,Cs_opt2
      close(1)
      write(fname,'(A12,i5.5,A4)')path//'output/t13',jt,'.out'
      open(1,file=fname,form='unformatted')
      write(1) txz
      close(1)
   end if
end if
If (mod(jt,c_count)==0) then
   call avgslice (txx,txz,tyy,tyz,tzz,p,jt)
end if
Cs_opt2_avg(:,:,:)=Cs_opt2_avg(:,:,:)+Cs_opt2(:,:,:)/nwrite
!TS
mean_u(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)=&
     mean_u(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)+&
     u(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)/nwrite
mean_u2(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)=&
     mean_u2(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)+&
     u(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)**2/nwrite
mean_v(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)=&
     mean_v(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)+&
     v(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)/nwrite
mean_v2(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)=&
     mean_v2(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)+&
     v(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)**2/nwrite
mean_dudx(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)=&
     mean_dudx(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)+&
     dudx(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)/nwrite
if (output) then
   if (mod(jt,nwrite)==0) then
      write(fname,'(A13,i5.5,A4)')path//'output/vel-',jt,'.out'
      open(21,file=fname,form='unformatted')
      write(21) real(u),real(v),real(w),real(Cs_opt2_avg),real(p)
      Cs_opt2_avg(:,:,:)=0._rprec
      close(21)
      write(95)real(mean_u),real(mean_u2),real(mean_v),real(mean_v2),real(mean_dudx)
      mean_u(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)=0_rprec
      mean_u2(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)=0_rprec
      mean_v(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)=0_rprec
      mean_v2(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)=0_rprec
      mean_dudx(jx_pls:jx_ple,jy_pls:jy_ple,2:nz-1)=0_rprec
   end if
   if(jt.gt.1000.and.mod(jt,nwrite/2).eq.0)then
      call lambda2()
      write(195)real(lam2)
   endif
end if
end subroutine output_loop

subroutine avgslice (txx,txz,tyy,tyz,tzz,p,jt)
use sim_param,only:path,u,v,w
use param,only:dz,p_count,c_count
use sgsmodule,only:Cs_opt2
implicit none
integer::i,j,k, jt
real(kind=rprec),dimension(ld,ny,nz)::txx,txz,tyy,tyz,tzz,p
real(kind=rprec),dimension(nx,nz)::ap,au,av,aw,p2,u2,v2,w2,auw,avw,acs
real(kind=rprec),dimension(nx,nz)::atxx,atxz,atyy,atyz,atzz
real(kind=rprec)::tu1,tv1,tw1,ttxx,ttxz,ttyy,ttyz,ttzz,&
     tu2,tv2,tw2,tp1,tp2,tuw,tvw,tCs,fr, arg1, arg2     

fr=(1._rprec/real(p_count,kind=rprec))*real(c_count,kind=rprec)
do k=1,Nz
do i=1,Nx
   tu1=0._rprec
   tv1=0._rprec
   tw1=0._rprec
   tp1=0._rprec
   ttxx=0._rprec
   ttxz=0._rprec
   ttyy=0._rprec
   ttyz=0._rprec
   ttzz=0._rprec
   tu2=0._rprec
   tv2=0._rprec
   tw2=0._rprec
   tp2=0._rprec
   tuw=0._rprec
   tvw=0._rprec
   tCs=0._rprec
   do j=1,Ny
      tu1=tu1+u(i,j,k)
      tv1=tv1+v(i,j,k)
      tw1=tw1+w(i,j,k)
      tp1=tp1+p(i,j,k)
      ttxx=ttxx+txx(i,j,k)
      ttxz=ttxz+txz(i,j,k)
      ttyy=ttyy+tyy(i,j,k)
      ttyz=ttyz+tyz(i,j,k)
      ttzz=ttzz+tzz(i,j,k)
      tu2=tu2+u(i,j,k)*u(i,j,k)
      tv2=tv2+v(i,j,k)*v(i,j,k)
      tw2=tw2+w(i,j,k)*w(i,j,k)
      tp2=tp2+p(i,j,k)*p(i,j,k)
      tCs=tCs+Cs_opt2(i,j,k)
      if (k.gt.1) then
         arg1=(u(i,j,k)+u(i,j,k-1))/2.
         arg2=(v(i,j,k)+v(i,j,k-1))/2.
      else
         arg1=0._rprec
         arg2=0._rprec
      end if
      tuw=tuw+w(i,j,k)*arg1
      tvw=tvw+w(i,j,k)*arg2
   end do
   au(i,k)=au(i,k)+(fr)*tu1/Ny
   av(i,k)=av(i,k)+(fr)*tv1/Ny
   aw(i,k)=aw(i,k)+(fr)*tw1/Ny
   ap(i,k)=ap(i,k)+(fr)*tp1/Ny
   u2(i,k)=u2(i,k)+(fr)*tu2/Ny
   v2(i,k)=v2(i,k)+(fr)*tv2/Ny
   w2(i,k)=w2(i,k)+(fr)*tw2/Ny
   atxx(i,k)=atxx(i,k)+(fr)*ttxx/Ny
   atxz(i,k)=atxz(i,k)+(fr)*ttxz/Ny
   atyy(i,k)=atyy(i,k)+(fr)*ttyy/Ny
   atyz(i,k)=atyz(i,k)+(fr)*ttyz/Ny
   atzz(i,k)=atzz(i,k)+(fr)*ttzz/Ny
   p2(i,k)=p2(i,k)+fr*tp2/Ny
   aCs(i,k)=aCs(i,k)+(fr)*tCs/Ny
   auw(i,k)=auw(i,k)+(fr)*tuw/Ny
   avw(i,k)=avw(i,k)+(fr)*tvw/Ny
end do
end do

if (mod(jt,p_count)==0.and.jt.gt.2000) then
   open(20,file=path//'output/aver_u.out')
   open(19,file=path//'output/aver_v.out')
   open(22,file=path//'output/aver_w.out')
   open(23,file=path//'output/aver_p.out')
!	open(24,file=path//'output/aver_u2.out')
!	open(25,file=path//'output/aver_v2.out')
!	open(26,file=path//'output/aver_w2.out')
   open(27,file=path//'output/aver_txx.out')
   open(28,file=path//'output/aver_txz.out')
   open(29,file=path//'output/aver_tyy.out')
   open(30,file=path//'output/aver_tyz.out')
   open(31,file=path//'output/aver_tzz.out')
   !	open(32,file=path//'output/aver_p2.out')
   open(33,file=path//'output/aver_uw.out')
   open(34,file=path//'output/aver_vw.out')
   open(35,file=path//'output/aver_Cs2.out')
   do k=1,nz
      write(20,5168) real(jt*dz),(real(au(i,k)),i=1,nx)
      write(19,5168) real(jt*dz),(real(av(i,k)),i=1,nx)
      write(22,5168) real(jt*dz),(real(aw(i,k)),i=1,nx)
      write(23,5168) real(jt*dz),(real(ap(i,k)),i=1,nx)
!	write(24,5168) jt*dz,(u2(i,k),i=1,nx)
!	write(25,5168) jt*dz,(v2(i,k),i=1,nx)
!	write(26,5168) jt*dz,(w2(i,k),i=1,nx)
      write(27,5168) real(jt*dz),(real(atxx(i,k)),i=1,nx)
      write(28,5168) real(jt*dz),(real(atxz(i,k)),i=1,nx)
      write(29,5168) real(jt*dz),(real(atyy(i,k)),i=1,nx)
      write(30,5168) real(jt*dz),(real(atyz(i,k)),i=1,nx)
      write(31,5168) real(jt*dz),(real(atzz(i,k)),i=1,nx)
!	write(32,5168) jt*dz,(p2(i,k),i=1,nx)
      write(33,5168) real(jt*dz),(real(auw(i,k)),i=1,nx)
      write(34,5168) real(jt*dz),(real(avw(i,k)),i=1,nx)
      write(35,5168) real(jt*dz),(real(aCs(i,k)),i=1,nx)
   end do

   au=0._rprec
   av=0._rprec
   aw=0._rprec
   ap=0._rprec
   u2=0._rprec
   v2=0._rprec
   w2=0._rprec
   atxx=0._rprec
   atxz=0._rprec
   atyy=0._rprec
   atyz=0._rprec
   atzz=0._rprec
   p2=0._rprec
   auw=0._rprec
   avw=0._rprec
   aCs=0._rprec   
end if
 5168     format(1400(E14.5))
end subroutine avgslice

subroutine output_final(jt)
use types,only:rprec
use param,only:nsteps
use sim_param,only:u,v,w,usp,vsp,wsp,RHSx,RHSy,RHSz
use sgsmodule,only:Cs_opt2
use sgsmodule,only:F_LM,F_MM
implicit none
integer,intent(in)::jt
integer::jx,jy,jz

rewind(121)
write(121) u,v,w,RHSx,RHSy,RHSz,Cs_opt2,F_LM,F_MM

usp=0._rprec
vsp=0._rprec
wsp=0._rprec
if (jt .GT. (nsteps-100)) then
   usp=usp+u
   vsp=vsp+v
   wsp=wsp+w
end if

usp=usp/100._rprec
vsp=vsp/100._rprec
wsp=wsp/100._rprec

rewind(64)
do jz=1,Nz
do jy=1,Ny
do jx=1,Nx
   write (64,5167) usp(jx,jy,jz),vsp(jx,jy,jz),wsp(jx,jy,jz)
end do
end do
end do
5167  format (110(1x,e10.4))
end subroutine output_final

subroutine lambda2()
use types,only:rprec
use sim_param,only:u,v,w,&
     dudx,dudy,dudz,dvdx,dvdy,dvdz,dwdx,dwdy,dwdz
use param,only:dx,dy,dz
!TSreal(kind=rprec),dimension(nx)::dpdx,dpdy,dpdz,aeu,awu,anu,asu,afu,abu
!TSreal(kind=rprec)::fne,fnn,fnf,S11,S22,S33,S12,O12,S13,O13,S23,O23
real(kind=rprec)::S11,S22,S33,S12,O12,S13,O13,S23,O23,&
     ux,uy,uz,vx,vy,vz,wx,wy,wz
integer::jx,jy,jz
! following used for eispack call...
integer::neis,nmeis,matzeis,ierreis,iv1eis(3)
double precision,dimension(3,3)::aeis,zeis
double precision,dimension(3)::wreis,wieis,fv1eis
double precision::ave

! assignments for eispack call
neis=3
nmeis=3
matzeis=0
ierreis=0
!TSfne=1._rprec/dx
!TSfnn=1._rprec/dy
!TSfnf=1._rprec/dz
!TSprint*,fne,fnn,fnf
lam2=0._rprec

! 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=ux          ! uvp-node
   S12=0.5_rprec*(uy+vx) ! uvp-node
! taken care of with wall stress routine
   S13=0.5_rprec*(uz+wx) ! uvp
   O12=0.5_rprec*(uy-vx) ! w-node
   O13=0.5_rprec*(uz-wx) ! w-node
   S22=vy          ! uvp-node
! taken care of with wall stress routine 
   S23=0.5_rprec*(vz+wy) ! uvp
   O23=0.5_rprec*(vz-wy) ! w-node
   S33=wz          ! uvp-node
   aeis(1,1)=s11*s11+s12*s12+s13*s13-O12*O12-O13*O13
   aeis(1,2)=s11*s12+s12*s22+s13*s23-O13*O23
   aeis(1,3)=s11*s13+s12*s23+s13*s33+O12*O23
   aeis(2,2)=s12*s12+s22*s22+s23*s23-O12*O12-O23*O23
   aeis(2,3)=s12*s13+s22*s23+s23*s33-O12*O13
   aeis(3,3)=s13*s13+s23*s23+s33*s33-O13*O13-O23*O23
   aeis(2,1)=aeis(1,2)
   aeis(3,1)=aeis(1,3)
   aeis(3,2)=aeis(2,3)
   call rg(nmeis,neis,aeis,wreis,wieis,matzeis,zeis,iv1eis,fv1eis,ierreis)
   if(wreis(1).ge.wreis(2).and.wreis(1).le.wreis(3)) then
      lam2(jx,jy,jz)=real(wreis(1),kind=rprec)
   elseif(wreis(2).ge.wreis(1).and.wreis(2).le.wreis(3)) then
      lam2(jx,jy,jz)=real(wreis(2),kind=rprec)
   else
      lam2(jx,jy,jz)=real(wreis(3),kind=rprec)
   endif
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=ux          ! w-node
   S12=0.5_rprec*(uy+vx) ! w-node
   S13=0.5_rprec*(uz+wx) ! w-node
   O12=0.5_rprec*(uy-vx) ! w-node
   O13=0.5_rprec*(uz-wx) ! w-node
   S22=vy          ! w-node
   S23=0.5_rprec*(vz+wy) ! w-node
   O23=0.5_rprec*(vz-wy) ! w-node
   S33=wz          ! w-node
   aeis(1,1)=s11*s11+s12*s12+s13*s13-O12*O12-O13*O13
   aeis(1,2)=s11*s12+s12*s22+s13*s23-O13*O23
   aeis(1,3)=s11*s13+s12*s23+s13*s33+O12*O23
   aeis(2,2)=s12*s12+s22*s22+s23*s23-O12*O12-O23*O23
   aeis(2,3)=s12*s13+s22*s23+s23*s33-O12*O13
   aeis(3,3)=s13*s13+s23*s23+s33*s33-O13*O13-O23*O23
   aeis(2,1)=aeis(1,2)
   aeis(3,1)=aeis(1,3)
   aeis(3,2)=aeis(2,3)
   call rg(nmeis,neis,aeis,wreis,wieis,matzeis,zeis,iv1eis,fv1eis,ierreis)
   if(wreis(1).ge.wreis(2).and.wreis(1).le.wreis(3)) then
      lam2(jx,jy,jz)=real(wreis(1),kind=rprec)
   elseif(wreis(2).ge.wreis(1).and.wreis(2).le.wreis(3)) then
      lam2(jx,jy,jz)=real(wreis(2),kind=rprec)
   else
      lam2(jx,jy,jz)=real(wreis(3),kind=rprec)
   endif
end do
end do
end do

!do k=2,nz-1
!do j=2,ny-1
!do i=2,nx-1
!   dpdx(i) = u(i+1,j,k) - u(i-1,j,k)
!   dpdy(i) = u(i,j+1,k) - u(i,j-1,k)
!   dpdz(i) = u(i,j,k+1) - u(i,j,k-1)
!   aeu(i)  = v(i+1,j,k) - v(i-1,j,k)
!   awu(i)  = v(i,j+1,k) - v(i,j-1,k)
!   anu(i)  = v(i,j,k+1) - v(i,j,k-1)
!   asu(i)  = w(i+1,j,k) - w(i-1,j,k)
!   afu(i)  = w(i,j+1,k) - w(i,j-1,k)
!   abu(i)  = w(i,j,k+1) - w(i,j,k-1)

!   S11=fne*dpdx(i)
!   S22=fnn*awu(i)
!   S33=fnf*abu(i)
!   S12=0.5_rprec*(fnn*dpdy(i)+fne*aeu(i))
!   O12=0.5_rprec*(fnn*dpdy(i)-fne*aeu(i))
!   S13=0.5_rprec*(fnf*dpdz(i)+fne*asu(i))
!   O13=0.5_rprec*(fnf*dpdz(i)-fne*asu(i))
!   S23=0.5_rprec*(fnf*anu(i)+fnn*afu(i))
!   O23=0.5_rprec*(fnf*anu(i)-fnn*afu(i))
!   aeis(1,1)=s11*s11+s12*s12+s13*s13-O12*O12-O13*O13
!   aeis(1,2)=s11*s12+s12*s22+s13*s23-O13*O23
!   aeis(1,3)=s11*s13+s12*s23+s13*s33+O12*O23
!   aeis(2,2)=s12*s12+s22*s22+s23*s23-O12*O12-O23*O23
!   aeis(2,3)=s12*s13+s22*s23+s23*s33-O12*O13
!   aeis(3,3)=s13*s13+s23*s23+s33*s33-O13*O13-O23*O23
!   aeis(2,1)=aeis(1,2)
!   aeis(3,1)=aeis(1,3)
!   aeis(3,2)=aeis(2,3)
!   call rg(nmeis,neis,aeis,wreis,wieis,matzeis,zeis,iv1eis,fv1eis,ierreis)
!   if(wreis(1).ge.wreis(2).and.wreis(1).le.wreis(3)) then
!      lam2(i,j,k)=real(wreis(1),kind=rprec)
!   elseif(wreis(2).ge.wreis(1).and.wreis(2).le.wreis(3)) then
!      lam2(i,j,k)=real(wreis(2),kind=rprec)
!   else
!      lam2(i,j,k)=real(wreis(3),kind=rprec)
!   endif
!enddo
!enddo
!enddo
print*,'minmax',minval(lam2),maxval(lam2)
end subroutine lambda2
end module io
