subroutine initial()
use types,only:rprec
use param,only:initu,interp,dns_bc,use_bldg,nz
use sim_param,only:path,u,v,w,RHSx,RHSy,RHSz
use sgsmodule,only:Cs_opt2_avg
use immersedbc,only:fx,fy,fz,u_des,v_des,w_des,n_bldg,bldg_pts
use io,only:mean_u,mean_u2,mean_v,mean_v2,mean_dudx
implicit none
integer::i
Cs_opt2_avg=0._rprec
fx=0._rprec;fy=0._rprec;fz=0._rprec
u_des=0._rprec;v_des=0._rprec;w_des=0._rprec
mean_u=0._rprec;mean_u2=0._rprec;mean_v=0._rprec;mean_v2=0._rprec;mean_dudx=0._rprec

open(51,file=path//'output/check_ke.out')
open(121,file=path//'vel.out',form='unformatted')

if (initu) then
   if (interp) then
      print*, 'Using interpolated field'
      open (31,file=path//'vel_new.out',form='unformatted')
      read(31) u,v,w
      close(31)
   else
      print *,'Reading initial velocity field from file'
      read(121) u,v,w,RHSx,RHSy,RHSz
   end if
else
   if (dns_bc) then
      print*, 'Creating initial velocity field with DNS BCs'
      call ic_dns()
   else
      print*, 'Creating initial velocity field'
      call ic()
   end if
end if

! bldg stuff
if (use_bldg) then
   open(1,file=path//'bldg.dat')
   read(1,*) n_bldg
   allocate(bldg_pts(5,n_bldg))
   do i=1,n_bldg
      read(1,*) bldg_pts(1:5,i)
      if(bldg_pts(5,i).ge.nz)then
         write(*,*)"lz should be less than nz"
         stop
      end if
   end do
end if
end subroutine initial
