program slabsetup
      implicit none

! we use 50 cubes of the base grid for the low-density region,
! and 8*50 reduced cubes of the base grid for the high-density part
      integer, parameter :: nl=50,nx=100, ny=2, nz=2
      integer, parameter :: nbase=1000,npart=nl*nbase+nx*ny*nz*nbase
      real,    dimension(3,nbase) :: grid
      real,    dimension(3,npart) :: pos,vel,mag
      real,    dimension(  npart) :: ein
      integer, dimension(  npart) :: id
      double precision :: masstab(6)
      integer :: i,j,k,l
      real :: gamma,mass,rho0h,rho0l,p0h,p0l,eps0h,eps0l

! initial conditions of medium
      gamma = 5./3.
      rho0h = 1.
      rho0l = 0.125
      p0h   = 1.
      p0l   = 0.1
      eps0h = p0h/(gamma-1.)/rho0h ! p = (gamma-1)*eps*rho
      eps0l = p0l/(gamma-1.)/rho0l
      mass  = rho0h/nbase/8 ! = rho0l/nbase

      forall(i=1:npart) id(i)=i
      vel = 0.

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

! set up the high-density part
      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)
            mag(1,l+1:l+nbase) = 0.75
            mag(2,l+1:l+nbase) = 1.0
            mag(3,l+1:l+nbase) = 0.0
            ein(l+1:l+nbase) = eps0h
          enddo
        enddo
      enddo

! set up the low-density part
      do i=1,nl
        l = nx*ny*nz*nbase + (i-1)*nbase
        pos(1,l+1:l+nbase) = grid(1,:) + 1.0*(i-1) + 50.
        pos(2,l+1:l+nbase) = grid(2,:)
        pos(3,l+1:l+nbase) = grid(3,:)
        mag(1,l+1:l+nbase) =  0.75
        mag(2,l+1:l+nbase) = -1.0
        mag(3,l+1:l+nbase) =  0.0
        ein(l+1:l+nbase) = eps0l
      enddo

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

! write the initial conditions file
      masstab = 0.
      masstab(1) = mass
      id = [(i,i=1,npart)]
      open(1,file='slab.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) 'BFLD',8+3*4*npart
      write(1) mag
      close(1)

      end