!        1         2         3         4         5         6         7
!23456789012345678901234567890123456789012345678901234567890123456789012
!
! 7/4/2015
! readif.f : code to read from vel.000000 and background.state
!            outputs one xz-plane of u', v', w', and theta'
!     NOTE : - this run had a domain of 40 km horizontal by 160 km
!              vertical.  The wavelength of the packet was 20 km.
!            - The dynamic viscosity is uniformly 1.5e-5 since an
!              isothermal atmosphere is used.  They kinimatic viscosity
!              is then just dynamic / density

      integer nx, ny, nz, nvarv, lx, ly, lz, nheadb
      integer i_rec_in, mean_recl, vel_recl
      character(72) str1, velfile, bsfile
      real(8),allocatable,dimension (:,:,:,:) :: volin
      real(8),allocatable,dimension (:,:    ) :: volu, volv, volw, volt
      real(8),allocatable,dimension (:,:    ) :: volup,volvp,volwp,voltp
      real(4),allocatable,dimension (:      ) :: meanu,meanv,meanw,meant
      real(4),allocatable,dimension (:      ) :: zb, rb, pb, tb, hb, nb         

! Declarations         
      nx    =  384
      ny    =    1  ! ny=1 for a 2d run, but there are 3 xz-planes written
      nz    = 1536  !      to vel.000000.  I believe they are identical.

      lx    = nx + 2
      ly    = ny + 2 
      lz    = nz + 2 

      nvarv =    8  ! vars 1-4 are u, v, w, t.

! Allocate arrays
      allocate( volin(lx,ly,lz,nvarv)                                  )
      allocate( volu(lx,lz) , volv(lx,lz) , volw(lx,lz) , volt(lx,lz)  )
      allocate( volup(lx,lz), volvp(lx,lz), volwp(lx,lz), voltp(lx,lz) )
      allocate( meanu(lz)   , meanv(lz)   , meanw(lz)   , meant(lz)    )
      allocate( zb(lz),  rb(lz),  pb(lz),  tb(lz),  hb(lz),  nb(lz)    )

! Open and read from files
!     vel.000000 is binary
!     background.state is ASCII

! file names
      velfile = 'vel.000000'
      bsfile  = 'background.state'

! vel.000000 is double precision
      i_rec_in = lx  * ly  * lz  * nvarv * 2      

      nheadb = 8

      close(100)
      close(200)
      open(unit=100,file=velfile,form='unformatted',access='direct',
     &       recl=i_rec_in, action='read')
      open(unit=200,file=bsfile,status='old')

! direct access read volin from vel.000000 (data stored as [x,y,z,var])
      read(100,rec=1) volin(:,:,:,:)    

! background.state : read past headers
      do k = 1, nheadb
        read(200,*) str1
      enddo
! background.state : read data (altitude, density, pressure, temp, p.temp, stability)
      do k = 1, lz
        read(200,*) zb(k), rb(k), pb(k), tb(k), hb(k), nb(k)
      enddo
! close files
      close(100)      
      close(200)      

! process (strip the density off and average)
      do k = 1, lz
        volu(:,k) = volin(:,2,k,1) / rb(k)
        volv(:,k) = volin(:,2,k,2) / rb(k) 
        volw(:,k) = volin(:,2,k,3) / rb(k) 
        volt(:,k) = volin(:,2,k,4) / rb(k) 
      
        meanu(k) = 0.
        meanv(k) = 0.
        meanw(k) = 0.
        meant(k) = 0.
        do i = 1, lx
          meanu(k) = meanu(k) + volu(i,k)/lx
          meanv(k) = meanv(k) + volv(i,k)/lx
          meanw(k) = meanw(k) + volw(i,k)/lx
          meant(k) = meant(k) + volt(i,k)/lx
        enddo
        do i = 1, lx
          volup(i,k) = volu(i,k) - meanu(k) 
          volvp(i,k) = volv(i,k) - meanv(k) 
          volwp(i,k) = volw(i,k) - meanw(k) 
          voltp(i,k) = volt(i,k) - meant(k) 
        enddo

      enddo


! write perturbation quantity files 
      close(301)
      close(302)
      close(303)
      close(304)

      vel_recl = lx * lz * 2

      open(unit=301,file='up.dat',form='unformatted',access='direct',
     &       recl=vel_recl, action='write')
      open(unit=302,file='vp.dat',form='unformatted',access='direct',
     &       recl=vel_recl, action='write')
      open(unit=303,file='wp.dat',form='unformatted',access='direct',
     &       recl=vel_recl, action='write')
      open(unit=304,file='tp.dat',form='unformatted',access='direct',
     &       recl=vel_recl, action='write')
      write(301,rec=1) volup(:,:)
      write(302,rec=1) volvp(:,:)
      write(303,rec=1) volwp(:,:)
      write(304,rec=1) voltp(:,:)
      close(301)
      close(302)
      close(303)
      close(304)


! write means
      close(401)
      close(402)
      close(403)
      close(404)

      mean_recl =      lz * 2

      open(unit=401,file='meanu.dat',form='unformatted',access='direct',
     &       recl=mean_recl, action='write')
      open(unit=402,file='meanv.dat',form='unformatted',access='direct',
     &       recl=mean_recl, action='write')
      open(unit=403,file='meanw.dat',form='unformatted',access='direct',
     &       recl=mean_recl, action='write')
      open(unit=404,file='meant.dat',form='unformatted',access='direct',
     &       recl=mean_recl, action='write')
      write(401,rec=1) meanu(:)
      write(402,rec=1) meanv(:)
      write(403,rec=1) meanw(:)
      write(404,rec=1) meant(:)
      close(401)
      close(402)
      close(403)
      close(404)


! END
      end