!TS 
! highly modified version with minimal output
! for use with fftw
! put in wisdom!
program main
use types,only:rprec
use param
use sim_param
use io,only:openfiles,avgslice,output_loop,output_final
use fft
use immersedbc
use sgsmodule,only:initest,Cs_opt2
use test_filtermodule
use topbc,only:setsponge,sponge
use bottombc,only:zo,patch,patchnum,num_patch,patches,avgpatch
implicit none
integer,parameter::wbase=1
integer :: jx,jy,jz,jt
real(kind=rprec),dimension(ld,ny,nz)::cx,cy,cz,&
     txx,txy,txz,tyy,tyz,tzz,divtx,divty,divtz
equivalence (txx,divtx),(tyy,divty),(tzz,divtz)
real(kind=rprec) rmsdivvel,ke,kestor,testu   
real(kind=rprec),dimension(nz)::u_ndim
real(kind=rprec),dimension(nx,ny):: T_s,q_s2d,&
     psi_m,psi_h,phi_h,phi_m,u_avg
real(kind=rprec)::u_avg_(num_patch)
real(kind=rprec),dimension(ld,ny,nz)::S_hat,Pr_0
real(kind=rprec)::const,tt

initest=1
call initial()
! formulate the fft plans--may want to use FFTW_USE_WISDOM
! initialize the kx,ky arrays
call init_fft()
! Open output files      
call openfiles()
! initialize test filter
!TSC standard dynamic or Lagrangian
   call test_filter_init(2._rprec*filter_size,G_test)
if (model == 3 .OR. model == 5) then  ! scale dependent dynamic
   call test_filter_init(4._rprec*filter_size,G_test_test)
end if
! sets roughness in the patches
! sc note: make sure that either patchnum is initialized (to 0) before
! calling, or that the patches routine initializes it to 0
call patches()

if (ubc == 1) call setsponge()
print*,'Number of timesteps', nsteps
print*,'dt and Cs= ',dt, cs
print*,'Nx, Ny, Nz= ',nx,ny,nz
print*,'Number of patches:',num_patch
!      print*,'sampling stats every ',c_count,' timesteps'
!      print*,'writing stats every ',p_count,' timesteps'
if (molec) print*, 'molecular viscosity (dimensional) ', nu_molec

! time Loop
do jt=1,nsteps   
! advance total time
   tt=tt+dt
! save previous time's right-hand-sides for Adams-Bashforth Integration
!  (In subroutine "STEP" use first order time advancement on first time 
!  step).
   if (jt > 1) then    
      RHSx_f=RHSx
      RHSy_f=RHSy
      RHSz_f=RHSz
   else ! for first timestep and for neutral conditions
      psi_m=0._rprec
      psi_h=0._rprec
      phi_m=1._rprec
      phi_h=1._rprec
   end if
! kill oddballs and calculate horizontal derivatives
   call filt_da(u,dudx,dudy)
   call filt_da(v,dvdx,dvdy)
   call filt_da(w,dwdx,dwdy)
! calculate average u at the wall over each patch
! don't need for dns case
   call avgpatch(u_avg)
!  finite differences      
   call ddz_uv(dudz,u)
   call ddz_uv(dvdz,v)
   call ddz_w(dwdz,w)   

! calculate wall stress 
   if (dns_bc) then
      call wallstress_dns(txz,tyz)
      call derivwall_dns()
   else
! "impose" wall stress
      call wallstress(txz,tyz)
! calculate derivatives at wall 
      call derivwall()
   end if
   
! compute convective term in physical space
   call convec(cx,cy,cz,u,v,w,dudy,dudz,dvdx,dvdz,dwdx,dwdy,jt)  
! compute turbulent viscosity (const.)
   if (dns_bc .and. molec) then
      call dns_stress(txx,txy,txz,tyy,tyz,tzz)
   else 
      call sgs_stag(jt,txx,txy,txz,tyy,tyz,tzz)
   end if
   if (use_bldg) then
      call bldg_stress(txx,txy,txz,tyy,tyz,tzz)
   end if
! compute divergence of SGS shear stresses     
! note: the divt's and the diagonal elements of t are equivalenced!
   call divstress_uv(divtx, txx, txy, txz)
   call divstress_uv(divty, txy, tyy, tyz)
   call divstress_w(divtz, txz, tyz, tzz)
! Compute preliminary RHS matrices for pressure calculation
   RHSx = -cx - divtx
   RHSy = -cy - divty
   RHSz = -cz - divtz
!        do jz=1,nz
!          do jy=1,ny
!            do jx=1,nx
!              RHSx(jx,jy,jz)=-cx(jx,jy,jz)-divtx(jx,jy,jz)
!	      RHSy(jx,jy,jz)=-cy(jx,jy,jz)-divty(jx,jy,jz)
!	      RHSz(jx,jy,jz)=-cz(jx,jy,jz)-divtz(jx,jy,jz) 
!!     &	         -sponge(jz)*w(jx,jy,jz)+beta(jx,jy,jz)  !stop
!            end do
!          end do
!        end do

! solve Poisson equation for pressure.
! do we ever need p itself, or only its gradient?  -> probably don't need
! store p
   if (use_bldg) then
      call press_stag_bldg(p,dpdx,dpdy,divtz)
   else
      call press_stag(p,dpdx,dpdy,divtz)     
   end if
! Calculate pressure gradients here
! careful! p is not dimensioned the same as the others!
   const=1._rprec/dz
   do jz=1,nz
   do jy=1,ny
   do jx=1,nx
      dpdz(jx,jy,jz)=const*(p(jx,jy,jz)-p(jx,jy,jz-1))
   end do
   end do
   end do

   call forcing(jt)
   call energy (ke)
   call rmsdiv(rmsdivvel,dudx,dvdy,dwdz)
   write(51,*) jt,ke
   if (mod(jt,wbase) == 0) then
      write(6,7777) float(jt), dt,rmsdivvel,ke
   end if
7777 format ('tt,dt,NU,diui,ke:',5(1x,e9.4))
   call output_loop(jt,txx,txy,txz,tyy,tyz,tzz)
! end time loop
end do
if (output)call output_final(jt)
end program main

