program min_rho_p_loc

c     Purpose:  Computes the minimum density, pressure and x,y,z location
c               of the minimum pressure from plot3d xyz and q files.

      implicit none

      character*7
     &     filex, fileq	

      integer
     &     icount, il, jl, kl, i, j, k, imin, jmin, kmin, n, nfn, iloc

      real
     &     ani, rho, ru, rv, rw, re, p, pdiff, rmin, pmin, gamma,
     &     gm1, tmp1, tmp2, tmp3, tim, yloc, ymin, ysum

      real, dimension(:), allocatable:: f1, xyz

      read(5,'(a)') filex
      read(5,'(a)') fileq
c      print *, filex, fileq 

      open( 1, file=filex, form='unformatted', status='old', err=100 )
      open( 2, file=fileq, form='unformatted', status='old', err=101 )

c-----Read grid file
         
      nfn = 3
      
      read(1) il, jl, kl
      
      allocate( xyz(1:il*jl*kl*nfn) )
      
      read(1) (xyz(i), i=1,il*jl*kl*nfn)
      
      close(1)
      
c-----Read solution file      
      
      nfn = 5

      read(2) il, jl, kl
      read(2) tmp1, tmp2, tmp3, tim

      allocate( f1(1:il*jl*kl*nfn) )

      read(2) (f1(i), i=1,il*jl*kl*nfn)
      
      close(2)

      rmin = 1.0e+10
      pmin = 1.0e+10
      yloc = 1.0e+10
      ymin = 1.0e+10
      gamma = 1.4
      gm1 = gamma - 1.0

      do k = 1,kl
        do j = 1,jl
          do i = 1,il
            n = 1
            iloc = (n - 1) * il * jl * kl + (k - 1) * il * jl +
     &             (j - 1) * il + i
            rho = f1(iloc)
            rmin = min(rho, rmin)

            n = 2
            iloc = (n - 1) * il * jl * kl + (k - 1) * il * jl +
     &             (j - 1) * il + i
            ru = f1(iloc)
	    
            n = 3
            iloc = (n - 1) * il * jl * kl + (k - 1) * il * jl +
     &             (j - 1) * il + i
            rv = f1(iloc)

            n = 4
            iloc = (n - 1) * il * jl * kl + (k - 1) * il * jl +
     &             (j - 1) * il + i
            rw = f1(iloc)

            n = 5
            iloc = (n - 1) * il * jl * kl + (k - 1) * il * jl +
     &             (j - 1) * il + i
            re = f1(iloc)

            p = gm1 * (re - 0.5 * (ru*ru + rv*rv + rw*rw) / rho)
            pmin = min(p * gamma, pmin) ! p/p_inf not p/(r_inf*a_inf^2)
	   
          end do
        end do
      end do

c-----Calculate location of min p
      
      icount = 0
      ysum = 0.0
      
      do k = 1,kl
        do j = 1,jl
          do i = 1,il
            n = 1
            iloc = (n - 1) * il * jl * kl + (k - 1) * il * jl +
     &             (j - 1) * il + i
            rho = f1(iloc)
            rmin = min(rho, rmin)

            n = 2
            iloc = (n - 1) * il * jl * kl + (k - 1) * il * jl +
     &             (j - 1) * il + i
            ru = f1(iloc)
	    yloc = xyz(iloc)
	    
            n = 3
            iloc = (n - 1) * il * jl * kl + (k - 1) * il * jl +
     &             (j - 1) * il + i
            rv = f1(iloc)

            n = 4
            iloc = (n - 1) * il * jl * kl + (k - 1) * il * jl +
     &             (j - 1) * il + i
            rw = f1(iloc)

            n = 5
            iloc = (n - 1) * il * jl * kl + (k - 1) * il * jl +
     &             (j - 1) * il + i
            re = f1(iloc)

            p = gamma*(gm1 * (re - 0.5 * (ru*ru + rv*rv + rw*rw) / rho))
	    

             pdiff = p-pmin
	    	    
	    if (pdiff .le. 0.0000001) then
c	      ymin = yloc
c	      print *, p, pmin, pdiff, ymin
             ysum = ysum + yloc
	      icount = icount + 1
	    end if      
	   
          end do
        end do
      end do
      
      ymin = ysum/real(icount)
      
      print*, tim, rmin, pmin, ymin

      deallocate( f1, xyz )

      stop

 100  print*, ' Error opening file: ', filex
 101  print*, ' Error opening file: ', fileq
      stop

      end