program nozzlesetup
      implicit none

! we use 60x2x2 cubes of our base grid
      integer, parameter :: nx=60, ny=2, nz=2
      integer, parameter :: nbase=1000, npart=nx*ny*nz*nbase
      real,    dimension(3,nbase) :: grid
      real,    dimension(3,npart) :: pos,vel
      real,    dimension(  npart) :: ein,mass
      integer, dimension(  npart) :: id
      double precision, dimension(6) :: masstab
      integer :: i,j,k,l
      real :: gamma,cs,rho0,eps0,p0,v0

! static properties of medium
      gamma = 5./3.
      mass  = 1.
      rho0  = 8000. ! 8000 particles of mass 1 per unit cube
      eps0  = 1.
      p0    = (gamma-1.)*eps0*rho0
      cs    = sqrt(gamma*p0/rho0)
      v0    = 0.1

! the base grid has coordinate values 0...1 in all three dimensions
      open(1,file='glass.txt')
      read(1,*) grid
      close(1)

! we scale the coordinates of our initial setup to 30x1x1
      do i=1,nx
        do j=1,ny
          do k=1,nz
            l = (ny*nz*(i-1)+nz*(j-1)+k-1)*nbase
            pos(1,l+1:l+nbase) = 0.5*grid(1,:) + 0.5*(i-1)
            pos(2,l+1:l+nbase) = 0.5*grid(2,:) + 0.5*(j-1)
            pos(3,l+1:l+nbase) = 0.5*grid(3,:) + 0.5*(k-1)
          enddo
        enddo
      enddo
      forall(i=1:npart) id(i)=i
      vel = 0.
      vel(1,:) = v0
      ein = eps0

! set up the nozzle
      do i=1,npart
        if(pos(1,i).gt.10.5 .and. pos(1,i).lt.14.5) then
          if(pos(3,i).lt.   (0.3*(pos(1,i)-10.5)/4.) .or. &
             pos(3,i).gt.1.-(0.3*(pos(1,i)-10.5)/4.)) then
            id(i) = 0
            vel(:,i) = 0. ! the nozzle does not move
          else
            vel(1,i) = v0/(1.-0.6*(pos(1,i)-10.5)/4.)
          endif
        elseif(pos(1,i).ge.14.5 .and. pos(1,i).lt.20.) then
          if(pos(3,i).lt.   0.3 .or. &
             pos(3,i).gt.1.-0.3) then
            id(i) = 0
            vel(:,i) = 0. ! the nozzle does not move
          else
            vel(1,i) = v0/0.4 ! higher velocity of flow through nozzle
          endif
        endif
      enddo

! set up the moving wall
      do i=1,npart
        if(pos(1,i).lt.0.5) then
          id(i) = 0
!         vel(1,i) = v0 ! already done above
        endif
      enddo

! increase the mass of the boundary particles
      do i=1,npart
        if(id(i).eq.0) mass(i) = 10.
      enddo

! dump coordinates to text file for debugging
      open(1,file='grid2.txt')
      write(1,'(i6,6f8.4)') (id(i),pos(:,i),vel(:,i),i=1,npart)
      close(1)

! write the initial conditions file
      masstab = 0.
      open(1,file='nozzle.ic',form='unformatted')
      write(1) 'HEAD',8+256
      write(1) [npart,0,0,0,0,0],masstab, &
        (0,i=1,(256-6*4-6*8)/4) ! pad to 256 bytes
      write(1) 'ID  ',8+4*npart
      write(1) id
      write(1) 'POS ',8+3*4*npart
      write(1) pos
      write(1) 'VEL ',8+3*4*npart
      write(1) vel
      write(1) 'U   ',8+4*npart
      write(1) ein
      write(1) 'MASS',8+4*npart
      write(1) mass
      close(1)


      end