c
c      gets velocity and pressure data for transonic duct solution 
c      for comparison with experimental data
c
      parameter (id=161,jd=97,kd=2,nbl=10,nvmax=10,nst=20)
c
      dimension x(id,jd,kd,3),q(id,jd,kd,nvmax),
     .idim(nbl),jdim(nbl),kdim(nbl),xprofile(nst),vel(nst,jd),
     .y(nst,jd),ptop(id),pbot(id),vel_exp(nst,jd),y_exp(nst,jd),
     .xbot_exp(id),xtop_exp(id),ptop_exp(id),pbot_exp(id),
     .x_exp(nst),nyval(nst)
c
      character*60 gfile,qfile
      character*80 string
c
      gamma = 1.4
c
c     dimensionalization values
c
      xref   = 0.1444 !ft
      astar  = 1029.8 !ft/sec
      pstar  = 1488.7 !lb/ft**2
      ptotal = 2819.5 !lb/ft**2
c
c     locations for experimental velocity profiles
c
      nsta = 4
      xprofile(1) = 2.882
      xprofile(2) = 4.611
      xprofile(3) = 6.340
      xprofile(4) = 7.493
c
c     write(6,*) 'name of plot3d grid file to read from? '
c     read(5,'(a60)') gfile
c     write(6,*) 'name of plot3d solution file to read from? '
c     read(5,'(a60)') qfile
c
      gfile='plot3dg.bin'
      qfile='plot3dq.bin'
c
      open(unit=3,file=gfile,form='unformatted',status='old')
      open(unit=4,file=qfile,form='unformatted',status='old')
c
c     read grid file to get number of blocks
c
      read(3) ngd
      if (ngd.gt.nbl) then
         write(6,*)'stopping...must set parameter nmax to at least ',
     .   nbl
         stop
      end if
c
c     loop over all blocks
c
      do 100 n=1,ngd
c
c     read grid file
c
      if (n.eq.1) read(3) (idim(nn),jdim(nn),nn=1,ngd)
      kdim(n)=1
      if (n.eq.1) then
         idmax = id
         jdmax = jd
         kdmax = kd
         do 5 nn=1,ngd
         if (idim(nn).gt.idmax) idmax = idim(nn)
         if (jdim(nn).gt.jdmax) jdmax = jdim(nn)
         if (kdim(nn).gt.kdmax) kdmax = kdim(nn)
 5       continue
         if (idmax.gt.id .or. jdmax.gt.jd .or. kdmax.gt.kd) then
            if (idmax.gt.id) then
               write(6,*)'stopping: must set parameter id to ',
     .                   'at least ',idmax
            end if
            if (jdmax.gt.jd) then
               write(6,*)'stopping: must set parameter jd to ',
     .                   'at least ',jdmax
            end if
            if (kdmax.gt.kd) then
               write(6,*)'stopping: must set parameter kd to ',
     .                   'at least ',kdmax
            end if
            stop
         end if
      end if
      ni=idim(n)
      nj=jdim(n)
      nk=kdim(n)
      read(3) (((x(i,j,k,1),i=1,ni),j=1,nj),k=1,nk),
     .        (((x(i,j,k,2),i=1,ni),j=1,nj),k=1,nk),
     .        (((idum,i=1,ni),j=1,nj),k=1,nk)
c
c     read solution file
c
      if (n.eq.1) read(4) ngd
      if (n.eq.1) read(4) (idim(nn),jdim(nn),nn=1,ngd)
      read(4) fsmach,alpha,re,time
      ni=idim(n)
      nj=jdim(n)
      nk=kdim(n)
      read(4) ((((q(i,j,k,l),i=1,ni),j=1,nj),k=1,nk),l=1,4)
c
c     convert plot3d quantities to primitive
c
      do 110 j=1,jdim(n)
      do 110 k=1,kdim(n)
      do 110 i=1,idim(n)
      rho  = q(i,j,k,1)
      u    = q(i,j,k,2)/rho
      v    = q(i,j,k,3)/rho
      vsq  = u*u + v*v
      p    = (gamma-1.)*(q(i,j,k,4)-.5*rho*vsq)
      csq  = gamma*p/rho
      dmsq = vsq/csq
      mach = sqrt(dmsq)
      p0  = p*(1.+.5*(gamma-1.)*dmsq)**(gamma/(gamma-1.))
      q(i,j,k,2) = u
      q(i,j,k,3) = v
      q(i,j,k,4) = p
      q(i,j,k,5) = mach
      q(i,j,k,6) = p0
 110  continue
c
c     interpolate velocity profiles to experimental x-stations
c
      do 200 ns=1,nsta
         iloc = -1
         do 210 i=1,idim(n)-1
         if(x(i,1,1,1) .lt. xprofile(ns) .and. 
     .      x(i+1,1,1,1) .ge. xprofile(ns)) iloc = i
 210     continue
         if (iloc.eq.-1) then
            write(6,*)'stopping...cannot find x location in ',
     .      'grid for data station ',n
            stop
         end if
         x2 = x(iloc+1,1,1,1)
         x1 = x(iloc,1,1,1)
         do 220 j=1,jdim(n)
         vel(ns,j) = q(iloc,j,1,2) + (xprofile(ns)-x1)/(x2-x1)*
     .               (q(iloc+1,j,1,2)-q(iloc,j,1,2))
         y(ns,j)   = x(iloc,j,1,2) + (xprofile(ns)-x1)/(x2-x1)*
     .               (x(iloc+1,j,1,2)-x(iloc,j,1,2))
c        make velocity dimensional
         vel(ns,j) = vel(ns,j)*astar*0.3048 !0.3048 converts ft/s to m/s
 220     continue 
c        make y dimensionless wrt local duct height
         ymax = y(ns,jdim(n))
         do j=1,jdim(n)
            y(ns,j) = y(ns,j)/ymax
         end do
 200  continue
c
c     top and bottom wall pressures (normalized with total pressure)
c     q(4) = p  normalized so that q(4) = 1/gamma at cfl3d reference
c     conditions
c
      do 230 i=1,idim(n)
      ptop(i) = gamma*q(i,jdim(n),1,4)*pstar/ptotal
      pbot(i) = gamma*q(i,1,1,4)*pstar/ptotal
 230  continue
c
 100  continue
c
c     read in exp. data files (tecplot), and write out to new files so that
c     the computational data can be appended to it.
c
      write(6,*)'input 0 for weak-shock case, 1 for strong-shock case'
      read(5,*) icase
      if (icase.eq.0) then
         open(unit=10,file='sajwvel.data',form='formatted',
     .   status='old')
         open(unit=11,file='sajwpres.data',form='formatted',
     .   status='old')
      else
         open(unit=10,file='sajsvel.data',form='formatted',
     .   status='old')
         open(unit=11,file='sajspres.data',form='formatted',
     .   status='old')
      end if
c
      if (icase.eq.0) then
         open(unit=8,file='u_vs_y_ws.dat',form='formatted',
     .   status='unknown')
         open(unit=9,file='p_vs_x_ws.dat',form='formatted',
     .   status='unknown')
      else
         open(unit=8,file='u_vs_y_ss.dat',form='formatted',
     .   status='unknown')
         open(unit=9,file='p_vs_x_ss.dat',form='formatted',
     .   status='unknown')
      end if
c
      do nn=1,9999
         read(10,'(a80)',end=998) string
         write(8,'(a80)') string
      end do
 998  continue
      do nn=1,9999
         read(11,'(a80)',end=999) string
         write(9,'(a80)') string
      end do
 999  continue
c
c     write velocity profile data file (tecplot format)
c
      write(8,'(''Zone T="Computation", F=Point'')')
      do j=1,jdim(1)
         write(8,'(i3,f8.4,f8.2,f8.4,f8.2,f8.4,f8.2,f8.4,f8.2)') 
     .   j,(y(ns,j),vel(ns,j),ns=1,nsta)
      end do
c
c     write wall pressure data file (tecplot format)
c
      write(9,'(''Zone T="Computation", F=Point'')')
      do i=1,idim(1)
         write(9,*) x(i,jdim(1),1,1),pbot(i),x(i,jdim(1),1,1),ptop(i)
      end do
c
      write(6,*)'done'
      if (icase.eq.0) then
         write(6,*) 'tecplot file for velocity: u_vs_y_ws.dat'
         write(6,*) 'tecplot file for pressure: p_vs_x_ws.dat'
      else
         write(6,*) 'tecplot file for velocity: u_vs_y_ss.dat'
         write(6,*) 'tecplot file for pressure: p_vs_x_ss.dat'
      end if

      stop
      end