InitialY4.f90

Go to the documentation of this file.

program impingingjetaxisym

  use modexamplescommon
  use modbc
  use modparsearguments
  use modstdio
  use modplot3d_io
  implicit none

  integer(ik) :: i, j, k
  character(len=CMD_ARG_LENGTH), dimension(:), allocatable :: RawArguments
  type(t_cmd_opt), dimension(12) :: Options
  integer(ik) :: Nx, Ny
  integer(ik) :: OutputFormat
  real(rk) :: U, V
  real(rk),dimension(:,:), allocatable :: xfac, yfac
  integer(ik) :: i1, j1, k1, i2, j2, k2, IBcheck
  real(rk) :: rhoinv
  character(120) :: ExecPath,GridPath,SolPath
  character :: cmd*100
  integer(ik), dimension(2) :: nPoints
  real(rk), dimension(2) :: Length
  real(rk), dimension(:,:,:), allocatable :: X
  real(rk), dimension(:,:,:,:,:), pointer :: XTemp, QTemp, QauxTemp
  integer(ik), dimension(:,:), allocatable :: IB
  integer(ik), dimension(:,:,:,:), pointer :: SSIB
  real(rk), dimension(:,:,:), allocatable :: Q, Qaux,massFractions
  type(t_bc), dimension(:),allocatable :: BCs
  real(rk), parameter :: Gamma = 1.4_rk
  real (rk), parameter :: GasConstant = 8314.4621_rk !J/(K kmol)
  real(rk) :: YO2, YH2, YH, YH2O, YHO2, nitrogenOxygenMoleRatio
  real(rk) :: MwRef, Mwa 
  real(rk) :: MachFactor,GL_csq
  real(rk) :: xtilde, ytilde!, Vtemp
  real(rk) :: Density, DensInv, Lref
  real(rk) :: Ly,Lx
  real(rk), dimension(2) :: Velocity
  integer(ik), dimension(2) :: nSponge
  real(rk), dimension(:), allocatable :: Mwcoeff, MwcoeffInv
  Integer (ik) :: VolRateSpecified = 0
  real(rk)::exp_x,exp_y
  integer(ik)::useDimensional 

  real(rk), dimension(:), allocatable :: thermo_mwts
  real(rk), dimension(:), allocatable :: thermo_h0,thermo_cp0,thermo_Y
  real(rk) :: thermo_rho, thermo_p, thermo_e, thermo_h, thermo_T, thermo_rhoe
  real(rk) :: densRef, Rair 

  Integer (ik), parameter :: nspecies=8
  Integer (ik), parameter :: iH2 = 1
  Integer (ik), parameter :: iO2 = 2
  Integer (ik), parameter :: iH = 3
  Integer (ik), parameter :: iH2O = 4
  Integer (ik), parameter :: iHO2 = 5
  Integer (ik), parameter :: iH2O2 = 6
  Integer (ik), parameter :: iO = 7
  Integer (ik), parameter :: iOH = 8
  Integer (ik) :: iN2 =9 
  Integer (ik), parameter :: GridStretch = 1
  Integer (ik), parameter :: includeElectrode = 1

  integer(ik)::iplug,jplug
  integer(ik)::jSlottedJetStart,jSlottedJetEnd
  real(rk)::AmbientPressure,AmbientDensity
  real(rk),dimension(:,:),allocatable::P ,T
  integer(ik)::jSpongeInlet
  integer(ik)::jm
  real(rk)::ystartJet,yendJet,ymid,y
  integer(ik),parameter::Nrefine=1
  real(rk),parameter::ReferenceLength=1.0
  real(rk)::xplug,xpos
  integer(ik),parameter::depthSlottedJet=20
  real(rk)::Mjet
  real(rk):: CPR 
  real(rk):: CDR
  real(rk):: Po
  real(rk):: rhoo
  real(rk):: Pjet
  real(rk):: rhojet

  real(rk)::plugRadius
  real(rk)::slotWidth
  real(rk)::dx,dy 
  real(rk)::machCrossStream

  real(rk)::Tref,rhoref,cref,rhoYO2Ref,rhoH2ref,DensityRatio
  real(rk)::Lp,dummyVal
  character(len=1000)::filename,dummyline
  integer::s,nskip
  logical::AIR_INTO_AIR, HYDROGEN_INTO_AIR
  integer::jetStart
  real(rk)::uvel,vvel

  real(rk)::tmp,tmp2,pressure,temperature
  real(rk)::MwcoeffRho

  real(rk)::densityLIBREF,temperatureLIBREF,velocityLIBREF
  real(rk)::densityLIB,vvelLIB,uvelLIB,temperatureLIB,YOLIB,YO2LIB,YO2LIBREF

  integer :: SSptsX, SSptsY, SSNdim, SSngrid, SSftype, SSNvar
  Character(LEN=2) :: SSprec, SSgf, SSvf, SSiblocal
  real(rk) :: SStau(4)
  integer,pointer,dimension(:,:) :: SSND
  integer,allocatable,dimension(:,:) :: ND

  real(rk)::unstartLoc

  real(rk):: MachNumberPreShock
  real(rk):: DensityPreShock
  real(rk):: TemperaturePreShock
  real(rk):: SoundSpeedPreShock

  real(rk):: MachNumberPostShock
  real(rk):: DensityPostShock
  real(rk):: TemperaturePostShock
  real(rk):: SoundSpeedPostShock

  real(rk)::VerticalOffSetFromJet
  real(rk)::HorizontalOffSetFromCenterline
  integer::jSlotJetCenter,iPlaceLIB,jPlaceLIB

  !inputs from Gyu Sub [labeled GS] 
  real(rk),parameter::airPoGS=167.1_rk*1.e3 !kPa
  real(rk),parameter::airToGs=298.0_rk !Kelvin
  real(rk),parameter::fuelPoGS=386.6_rk*1.e3 !kPa
  real(rk),parameter::fuelToGS=298.0_rk !Kelvin
  real(rk),parameter::airMdotGS=56.38_rk !g/s
  real(rk),parameter::fuelMdotGS=1.842_rk !g/s 
  real(rk)::airRoGS,fuelRoGs,airCoGS,fuelCoGs !derived below from Gyu Sub data
  real(rk),parameter::testSectionMachNumber=2.7_rk 
  real(rk)::testSectionDensity,testSectionPressure,testSectionTemperature,testSectionSoundSpeed,testSectionVelocity
  real(rk),parameter::fuelMachNumber=1.0_rk
  real(rk)::fuelDensity,fuelPressure,fuelTemperature,fuelSoundSpeed,fuelVelocity

  !to nondimensionalize the pressure and temperature
  real(rk),parameter::nondimT=298.*(gamma-1.0)
  real(rk),parameter::nondimP=101.325 !kPa
  real(rk),parameter::nondimR=1.2041 !kg/m^3 !what was in the input deck
  real(rk),parameter::nondimVel=343. !m/s

  airrogs = airpogs  / airtogs / 259.8 !indiviaul gas constant for oxygen [259.8]
  fuelrogs = fuelpogs  / fueltogs / 4124. !individual gas constant for hydrogen [4124]
  aircogs = sqrt(gamma * airpogs  / airrogs)
  fuelcogs = sqrt(gamma * fuelpogs  / fuelrogs)

  print *,"**************************************************************"
  print *,"Print the Gyu Sub inputs"
  print *,"airPo [Pa] = ",airpogs
  print *,"airTo [K] = ",airtogs
  print *,"airRo [kg/m^3] (derived) = ",airrogs
  print *,"fuelPo [Pa] = ",fuelpogs
  print *,"fuelTo [K] = ",fueltogs
  print *,"fuelRo [kg/m^3] (derived) = ",fuelrogs
  print *,"**************************************************************"

  testsectiondensity=get_isentropic_density(testsectionmachnumber,airrogs,gamma)
  testsectionpressure=get_isentropic_pressure(testsectionmachnumber,airpogs,gamma) 
  testsectiontemperature=get_isentropic_temperature(testsectionmachnumber,airtogs,gamma)
  testsectionsoundspeed=get_isentropic_soundspeed(testsectionmachnumber,aircogs,gamma)
  testsectionvelocity=testsectionmachnumber*testsectionsoundspeed

  print *,"**************************************************************"
  print *,"Print isentropic values in test section (M=2.7)"
  print *,"testSectionP [Pa] (derived) = ",testsectionpressure
  print *,"testSectionT [K] (derived) = ",testsectiontemperature
  print *,"testSectionR [kg/m^3] (derived) = ",testsectiondensity
  print *,"testSectionVel [m/s] (derived) = ",testsectionvelocity
  print *,"Recheck Mach = ",testsectionvelocity/testsectionsoundspeed 
  print *,"**************************************************************"

  fueldensity=get_isentropic_density(fuelmachnumber,fuelrogs,gamma)
  fuelpressure=get_isentropic_pressure(fuelmachnumber,fuelpogs,gamma)
  fueltemperature=get_isentropic_temperature(fuelmachnumber,fueltogs,gamma)
  fuelsoundspeed=get_isentropic_soundspeed(fuelmachnumber,fuelcogs,gamma)
  fuelvelocity=fuelmachnumber*fuelsoundspeed

  print *,"**************************************************************"
  print *,"Print isentropic values exiting fuel (M=1.0)"
  print *,"fuelP [Pa] (derived) = ",fuelpressure
  print *,"fuelT [K] (derived) = ",fueltemperature
  print *,"fuelR [kg/m^3] (derived) = ",fueldensity
  print *,"fuelVel [m/s] (derived) = ",fuelvelocity
  print *,"Recheck Mach = ",fuelvelocity/fuelsoundspeed
  print *,"**************************************************************"

  print *,"**************************************************************"
  print *,"Momentum Ratio of Test Section to Jet"
  print*,"J=",fueldensity*fuelvelocity**2/testsectiondensity/testsectionvelocity**2
  print *,"**************************************************************"

  print *,"**************************************************************"
  print *,"Stagnation jet pressure to ambient test section pressure"
  print*,"PR=",fuelpogs/testsectionpressure
  print *,"**************************************************************"

  print *,"**************************************************************"
  print *,"jet pressure to ambient test section pressure"
  print*,"PR_jet=",fuelpressure/testsectionpressure
  print *,"**************************************************************"

  print *,"**************************************************************"
  print *,"Reynolds number of the test section"
  print*,"Re=0.008 * testSectionVel * testSectionRho/9.0865272E-6=",0.008*testsectionvelocity*testsectiondensity/9.0865272e-6
  print *,"**************************************************************"

  print *,"**************************************************************"
  print *,"Reynolds number of the jet exit"
  print*,"Re=0.000195 * fueVel * fuelRho/7.7724307E-6=",0.000195*fuelvelocity*fueldensity/7.7724307e-6
  print *,"**************************************************************"

  !9.0865272E-6
  
  !get location of executable
  call get_command_argument(0, cmd)
  execpath(1:len_trim(cmd))=cmd(1:len_trim(cmd))
  !get command line arguments
  allocate(rawarguments(command_argument_count()))
  do i = 1, size(rawarguments)
    call get_command_argument(i, rawarguments(i))
  end do

  options(1) = t_cmd_opt_("format", "f", opt_value_type_string) !HDF5 or PLOT3D I/O
  options(2) = t_cmd_opt_("nptsx", "i", opt_value_type_integer) !number of points in x direction
  options(3) = t_cmd_opt_("nptsy", "j", opt_value_type_integer) !number of points in x direction
  options(4) = t_cmd_opt_("ly", "H", opt_value_type_real) !height of domain
  options(5) = t_cmd_opt_("lx", "W", opt_value_type_real) !width of domain
  options(6) = t_cmd_opt_("plug radius ", "r", opt_value_type_real) !plug radius
  options(7) = t_cmd_opt_("slot width ", "w", opt_value_type_real) !slotted jet width
  options(8) = t_cmd_opt_("cross stream Mach", "M", opt_value_type_real) !slotted jet width
  options(9) = t_cmd_opt_("dimensional or nondimensional", "D", opt_value_type_integer)!slotted jet width 
  options(10) = t_cmd_opt_("pressure ratio", "p",opt_value_type_real)!slotted jet width
  options(11) = t_cmd_opt_("x", "x",opt_value_type_real)!slottedjet width
  options(12) = t_cmd_opt_("y", "y",opt_value_type_real)!slottedjet width


  call parsearguments(rawarguments, options=options)
  
  !I/O Stuff
  if (options(1)%is_present) then
    if (options(1)%value == "plot3d") then
      outputformat = output_format_plot3d
    else if (options(1)%value == "hdf5") then
      outputformat = output_format_hdf5
    else
      write (error_unit, '(3a)') "ERROR: Unrecognized grid format '", trim(options(1)%value), "'."
      stop 1
    end if
  else
    outputformat = defaultoutputformat
  end if

  if (options(2)%is_present) then
    read (options(2)%value, *) nx
  else
    nx = 500*nrefine 
  end if

  if (options(3)%is_present) then
    read (options(3)%value, *) ny
  else
    ny = 1000*nrefine 
  end if

  if (options(4)%is_present) then
    read (options(4)%value, *) ly
  else
    ly = 56.0_rk*referencelength 
  end if

  if (options(5)%is_present) then
    read (options(5)%value, *) lx
  else
    lx = 14.0_rk*referencelength  
  end if
  
  if (options(6)%is_present) then
    read (options(6)%value, *) plugradius
  else
    plugradius = 4.0_rk*referencelength
  end if

  if (options(7)%is_present) then
    read (options(7)%value, *) slotwidth
  else
    slotwidth = 5.0_rk*referencelength
  end if

  if (options(8)%is_present) then
    read (options(8)%value, *) machcrossstream
  else
    machcrossstream = 0._rk
  end if

  if (options(10)%is_present) then
    read (options(10)%value, *) densityratio
  else
    densityratio=1.
  end if

!LIB POSITIONS RELATIVE TO THE JET
if (options(11)%is_present) then
    read (options(11)%value, *) verticaloffsetfromjet
  else
    verticaloffsetfromjet=-4.00
  end if

if (options(12)%is_present) then
    read (options(12)%value, *) horizontaloffsetfromcenterline
  else
    horizontaloffsetfromcenterline=7.43
  end if

  !grid creation
  npoints = [nx, ny]
  length = [lx, ly]

  dx=length(1)/real(npoints(1)-1)
  dy=length(2)/real(npoints(2)-1)

  allocate(x(2,npoints(1),npoints(2)))
  allocate(ib(npoints(1),npoints(2)))
  ib(:,:) = 1

  !simple uniform mesh stretching
  do j = 1, npoints(2)
    do i = 1, npoints(1)
      xtilde = real(i-1, kind=rk)/real(npoints(1)-1, kind=rk)
      ytilde = real(j-1, kind=rk)/real(npoints(2)-1, kind=rk)
      u = xtilde 
      v = ytilde
      x(:,i,j) = (length * [u,v])- [0.0_rk,0.0_rk]
    end do
  end do

  !allocate Q vars
  allocate(q(4,npoints(1),npoints(2)))
  allocate(qaux(nspecies,npoints(1),npoints(2)))
  allocate(massfractions(nspecies,npoints(1),npoints(2)))
  allocate(p(npoints(1),npoints(2)))
  allocate(t(npoints(1),npoints(2)))

  !plug (i.e. the 'sting') radius is 4 mm
  !start the jet exit at 5\times 4 mm
  jslottedjetstart=int(5.0*plugradius/dy)
  jslottedjetend=jslottedjetstart+int(slotwidth/dy)
  print *,"Number of points across slot=",abs(jslottedjetstart-jslottedjetend)

  iplug=int(plugradius/dx)
  ystartjet=x(2,iplug,jslottedjetstart)
  yendjet=x(2,iplug,jslottedjetend)

  !iblank the plug
  do j = 1, npoints(2)
  do i = 1, iplug-1
  ib(i,j)=0.
  end do
  end do
 
  jetstart=7*iplug/8
  do j = jslottedjetstart,jslottedjetend
  do i = jetstart, iplug
  ib(i,j)=1.
  end do
  end do

  !copied directly from jet in cross flow initialization per model0 specification
  in2 = nspecies + 1
  allocate(mwcoeff(nspecies+1), mwcoeffinv(nspecies+1))
  mwcoeff(ih2) = 2.0;
  mwcoeff(io2) = 32.0;
  if  (nspecies>=3) mwcoeff(ih) = 1.0_rk
  if  (nspecies>=4) mwcoeff(ih2o) = 18.0_rk
  if  (nspecies>=5) mwcoeff(iho2) = 33.0_rk
  if  (nspecies>=6) mwcoeff(ih2o2) = 34.0_rk
  if  (nspecies>=7) mwcoeff(io) = 16.0_rk
  if  (nspecies>=8) mwcoeff(ioh) = 17.0_rk
  mwcoeff(in2) = 28.0_rk
  nitrogenoxygenmoleratio = 3.76_rk
  yo2 = 1.0_rk/(1.0_rk+7.0_rk/8.0_rk*nitrogenoxygenmoleratio);
  mwa =  (mwcoeff(io2) +  mwcoeff(in2)*nitrogenoxygenmoleratio)/(1.0_rk+nitrogenoxygenmoleratio)
  rair = gasconstant / mwa
  cref = sqrt(gamma*rair*tref)
  cref = cref / machfactor
  do i=1,nspecies+1
    mwcoeff(i) = 1.0_rk / mwcoeff(i)
  enddo
  mwref = yo2*mwcoeff(io2) + (1.0_rk-yo2)*mwcoeff(in2)
  mwcoeffrho=1.0_rk/mwcoeff(in2)/mwref
  mwcoeff = mwcoeff/mwref
  do i=1,nspecies
    mwcoeff(i) = mwcoeff(i)-mwcoeff(in2)
  enddo

  !some gut checks to make sure `luca code' makes physical sense
  print *,"**************************************************************"
  print *,"some gut checks to make sure `luca code' makes physical sense"
  print *,"**************************************************************"
  massfractions(:,:,:)=0.0_rk !pure nitrogen
  massfractions(in2,:,:)=1.0_rk !pure nitrogen 
  print*,"density of nitrogen[kg/m^3]=",1.0_rk/(mwcoeff(in2)+sum(mwcoeff(1:nspecies)*massfractions(1:nspecies,1,1)))*nondimr
  massfractions(:,:,:)=0.0_rk !pure nitrogen
  massfractions(ih2,:,:)=1.0_rk !pure hydrogen 
  print*,"density of hydrogen[kg/m^3]=",1.0_rk/(mwcoeff(in2)+sum(mwcoeff(1:nspecies)*massfractions(1:nspecies,1,1)))*nondimr
  massfractions(:,:,:)=0.0_rk 
  massfractions(io2,:,:)=yo2               
  massfractions(in2,:,:)=1.0-yo2 
  print *,"density of air[kg/m^3]=",1.0_rk/(mwcoeff(in2) +sum(mwcoeff(1:nspecies)*massfractions(1:nspecies,1,1)))*nondimr
  print *,"**************************************************************"

  print *,"**************************************************************"
  print *,"ad-hoc boundary layer of delta_99 = slot width"
  print *,"The functional form=tanh(2.646652412/(slotwidth)*(xpos-*xplug))"
  print *,"**************************************************************"
  cref=1.0
  massfractions(:,:,:)=0._rk
  massfractions(io2,:,:)=1.0_rk !YO2
  xplug=x(1,iplug,1)

  !density(l0)*temperature(l0) / gamma

  do i=1,npoints(1)
  do j=1,npoints(2) 

   if (i .ge. iplug) then
   xpos=x(1,i,j)
   q(1,i,j)=testsectiondensity/nondimr !1.0_rk/(Mwcoeff(iN2) + sum(Mwcoeff(1:nspecies)*massFractions(1:nspecies,i,j)))
   q(3,i,j)=q(1,i,j)*testsectiondensity/nondimvel*tanh(2.646652412/(slotwidth)*(xpos-xplug))
   q(4,i,j)=q(1,i,j)*testsectiontemperature/nondimt/gamma + 0.5_rk * q(3,i,j)**2/q(1,i,j)
   qaux(1:nspecies,i,j)=q(1,i,j) * massfractions(1:nspecies,i,j)
   else
   q(1,i,j)=testsectiondensity/nondimr !1.0_rk/(Mwcoeff(iN2) + sum(Mwcoeff(1:nspecies)*massFractions(1:nspecies,i,j)))
   q(2,i,j)=0. 
   q(3,i,j)=0. 
   q(4,i,j)=q(1,i,j)*testsectiontemperature/nondimt/gamma + 0.5_rk * q(3,i,j)**2/q(1,i,j)
   qaux(1:nspecies,i,j)=q(1,i,j) * massfractions(1:nspecies,i,j)
   end if

  end do
  end do


!This adds in the LIB on top of the provided Q solution
if (.false.) then 

   !read grid from example directory
   !GridPath = trim(ExecPath)//trim(TemplateFlnm)
   ssprec = ""; ssgf = ""; ssvf = ""; ssiblocal = ""
   print *,"made it here"
   print *,"sizeof(SSND)",sizeof(ssnd)
   allocate(ssnd(1,3))
   print *,ssnd
   deallocate(ssnd)
   print *,"going to read the mesh"
   gridpath=trim("grid.xyz")
   solpath=trim("solution.q")
   call read_grid_size(ssndim, ssngrid, ssnd, gridpath, 1, ssiblocal)
   Call read_grid(ssndim, ssngrid, ssnd, xtemp, ssib,ssprec,ssgf,ssvf,ssiblocal,gridpath, 1)
   print*,'SS grid X extent is',minval(xtemp(:,:,:,:,1)),maxval(xtemp(:,:,:,:,1))
   print*,'SS grid Y extent is',minval(xtemp(:,:,:,:,2)),maxval(xtemp(:,:,:,:,2))
   Call read_soln_header(ssndim, ssngrid, ssnd,solpath, sstau)
   print*,'Read soln, size is', ssndim,ssngrid,ssnd,sstau
   Call read_soln(ssndim, ssngrid, ssnd, qtemp, sstau, ssprec, ssgf, ssvf,solpath, 1)
   print*,'Read soln, size is', shape(qtemp)

   do j=1,npoints(2)
   do i=1,npoints(1)
    q(1,i,j) = qtemp(1,i,j,1,1) 
    q(2,i,j) = qtemp(1,i,j,1,2) 
    q(3,i,j) = qtemp(1,i,j,1,3) 
    q(4,i,j) = qtemp(1,i,j,1,5) 
   end do
   end do

  jslotjetcenter=int((jslottedjetstart+jslottedjetend)/2)
  verticaloffsetfromjet=verticaloffsetfromjet+xtemp(1,1,jslotjetcenter,1,2)
  iplacelib=0
  jplacelib=0
  do j=1,npoints(2)
  do i=1,npoints(1)
    if (xtemp(1,i,j,1,1).le.horizontaloffsetfromcenterline) iplacelib=i
    if (xtemp(1,i,j,1,2).le.verticaloffsetfromjet)jplacelib=j 
  end do
  end do

  print *,"i,j of LIB",iplacelib,jplacelib
  print*,xtemp(1,iplacelib,jplacelib,1,1),xtemp(1,iplacelib,jplacelib,1,2)-xtemp(1,iplacelib,jslotjetcenter,1,2)
 
  massfractions(:,:,:)=0. 
  filename='128x128.dat'
  print *,"Reading lib data file: ",trim(filename)
  nskip=18
  OPEN(1,file=trim(filename))
  do s=1,nskip
  READ(1,*) dummyline
  end do
  do j=jplacelib-64,jplacelib+63
   do i=iplacelib-64,iplacelib+63
     xpos=x(1,i,j)
     read(1,*),dummyval,dummyval,dummyval,&
               dummyval,& !Qaux(7,i,j),&
               dummyval,& !Qaux(2,i,j),&
               uvellib,&! Q(2,i,j),&
               vvellib,&!Q(3,i,j),&
               dummyval,& !ignore the w-velocity
               temperaturelib,& !Q(4,i,j),&
               densitylib,&  !Q(1,i,j)
               yolib,&
               yo2lib

      if (temperaturelib .ge. 300.) then

      yo2libref=0.232917502522468567 
      densitylibref=1.22
      temperaturelibref=288.
      velocitylibref=343. 

      q(4,i,j)=1.89*temperaturelib/temperaturelibref*tref/gamma*q(1,i,j) + 0.5_rk *q(3,i,j)**2/q(1,i,j)+ 0.5_rk * q(2,i,j)**2/q(1,i,j)

      end if
    end do
  end do
  close(1)

end if
  
  call writeplot3dsolution("initial.q", ndims=2, npoints=npoints, q=q)
  call writeplot3dfunc("initial.aux.q", ndims=2, npoints=npoints, f=qaux)

  !initializing the target data for the sonic hydrogen jet 
  !underexapnded with the pressure from the input deck
  massfractions(:,:,:)=0._rk
  massfractions(ih2,:,:)=1.0_rk
  do j=jslottedjetstart,jslottedjetend
  do i=jetstart,iplug
   y=x(2,i,j)
   xpos=x(1,i,j)
   q(1,i,j)=fueldensity/nondimr !DensityRatio/(Mwcoeff(iN2) +sum(Mwcoeff(1:nspecies)*massFractions(1:nspecies,i,j)))
   q(2,i,j)=q(1,i,j)*fuelvelocity/nondimvel!*(tanh(10.*(y-ystartjet))-tanh(10.*(y-yendJet))-1.)
   q(3,i,j)=0._rk
   q(4,i,j)=q(1,i,j)*fueltemperature/nondimt/gamma+0.5*(q(2,i,j)**2)/q(1,i,j)
   qaux(1:nspecies,i,j)=q(1,i,j)*massfractions(1:nspecies,i,j)  
  end do
  end do 

  call writeplot3dsolution("initial.target.q", ndims=2, npoints=npoints, q=q)
  call writeplot3dfunc("initial.target.aux.q", ndims=2, npoints=npoints, f=qaux)

  print *,"**************************************************************"
  print *,"The ambient pressure to jet pressure ratio="!,!P 
  print *,"**************************************************************"

  allocate(bcs(11))
  bcs(1) =  t_bc_(1, bc_type_sat_slip_adiabatic_wall, -1, [-1,-1], [1,-1])
  bcs(2) =  t_bc_(1, bc_type_sat_far_field, 2, [iplug,-1], [1,1])
  bcs(3) =  t_bc_(1, bc_type_sat_far_field, -2, [iplug,-1], [-1,-1])
  bcs(4) =  t_bc_(1, bc_type_sat_slip_adiabatic_wall,1,[iplug,iplug],[1,jslottedjetstart])
  bcs(5) =  t_bc_(1, bc_type_sat_slip_adiabatic_wall,1,[iplug,iplug],[jslottedjetend,-1])
  bcs(6) =  t_bc_(1, bc_type_sponge, -2, [iplug,-1],[npoints(2)-int(plugradius/dy),-1])
  bcs(7) =  t_bc_(1, bc_type_sponge, 2, [iplug,-1], [1,int(plugradius/dy)])
  bcs(8) =  t_bc_(1, bc_type_sat_far_field,1,[jetstart,jetstart],[jslottedjetstart,jslottedjetend])
  bcs(9) =  t_bc_(1, bc_type_sat_slip_adiabatic_wall,2,[jetstart,iplug],[jslottedjetstart,jslottedjetstart])
  bcs(10) = t_bc_(1, bc_type_sat_slip_adiabatic_wall,-2,[jetstart,iplug],[jslottedjetend,jslottedjetend])
  bcs(11) = t_bc_(1, bc_type_sponge, 1,[jetstart,iplug],[jslottedjetstart,jslottedjetend])

  call writeplot3dgrid("grid.xyz", ndims=2, npoints=npoints, x=x, iblank = ib)
  call writebcs(bcs, "bc.dat")

contains

!consider having total (stagnation properties as inputs)
real function get_isentropic_pressure(M,Po,gamma)
real(rk)::m
real(rk)::Po
real(rk)::gamma
get_isentropic_pressure=0.
get_isentropic_pressure=(1.+(gamma-1.)/2.*m**2)**(-gamma/(gamma-1.))*po
end function
real function get_isentropic_density(M,rhoo,gamma)
real(rk)::M
real(rk)::rhoo
real(rk)::gamma
get_isentropic_density=0.
get_isentropic_density=(1.+(gamma-1.)/2.*m**2)**(-1./(gamma-1.))*rhoo
end function
real function get_isentropic_temperature(M,To,gamma)
real(rk)::M
real(rk)::To
real(rk)::gamma
get_isentropic_temperature=0.
get_isentropic_temperature=(1.+(gamma-1.)/2.*m**2)**(-1)*to
end function
real function get_isentropic_soundspeed(M,co,gamma)
real(rk)::M
real(rk)::co
real(rk)::gamma
get_isentropic_soundspeed=0.
get_isentropic_soundspeed=(1.+(gamma-1.)/2.*m**2)**(-1./2.)*co
end function
 
end program impingingjetaxisym

Subroutine p3d_read_soln_header(NDIM, ngrid, ND, fname, tau)

    ! ... function call variables
    Integer :: NDIM, ngrid
    Integer, Dimension(:,:), Pointer :: ND
    Character(LEN=*) :: fname
    Real(KIND=8) :: tau(4)

    ! ... local variables
    Integer :: ibuf, I, J

    ! ... open the file
    Open (unit=10, file=trim(fname), form='unformatted', status='old')

    ! ... read the number of grids, exit if error
    Read (10) ibuf
    print *,ibuf

    ! ... read the size of each grid, throw value away
    Read (10) ((ibuf,j=1,ndim),i=1,ngrid)
    print *,ndim

    ! ... read the header
    Read (10) tau
    print *,tau

    ! ... close and exit
    Close (10)

  End Subroutine p3d_read_soln_header

  Subroutine p3d_read_soln(NDIM, ngrid, ND, X, tau, prec, gf, vf, file, OF)

    Implicit None

    Real(KIND=8), Dimension(:,:,:,:,:), Pointer :: X
    Real(KIND=8) :: tau(4)
    Character(LEN=2) :: prec, gf, vf, ib
    Character(LEN=80) :: file
    Integer :: ngrid, I, J, K, L, M, Nx, Ny, Nz, NDIM, ftype
    Integer, Dimension(:,:), Pointer :: ND
    Real(KIND=4), Dimension(:,:,:,:), Pointer :: fX
    Real(KIND=4) :: ftau(4)
    Logical :: gf_exists
    Integer :: OF

    Inquire(file=trim(file),exist=gf_exists)
    If (.NOT.gf_exists) Then
      Write (*,'(A,A,A)') 'File ', file(1:len_trim(file)), ' is not readable'
      stop
    End If

    prec = ""; gf = ""; vf = ""; ib = ""
    Call p3d_detect(len(trim(file)),trim(file),prec,ndim,gf,vf,ib,ftype)

    Open (unit=10, file=trim(file), form='unformatted', status='old')

    !  Read number of grids
    If (gf(1:1) .eq. 's') Then
      ngrid = 1
    Else
      Read (10) ngrid
    End If


    !  Read grid sizes
    Allocate(nd(ngrid,3))
    nd(:,3) = 1
    Read (10) ((nd(i,j),j=1,ndim),i=1,ngrid)

    !  Allocate soln memory
    nx = maxval(nd(:,1))
    ny = maxval(nd(:,2))
    If (ndim .EQ. 3) Then
      nz = maxval(nd(:,3))
    Else
      nz = 1
    End If
    If (associated(x)) deallocate(x)
    Allocate(x(ngrid,nx,ny,nz,ndim+2))
    If (prec(1:1) .eq. 's') Allocate(fx(nx,ny,nz,ndim+2))
    x = 0d0

    !  Read soln values
    If (vf(1:1) .eq. 'w') Then
      If (prec(1:1) .eq. 'd') Then
        Do l = 1, ngrid
          Read (10) tau(1), tau(2), tau(3), tau(4)
          Read (10)((((x(l,i,j,k,m),i=1,nd(l,1)),j=1,nd(l,2)),k=1,nd(l,3)),m=1,ndim+2)
        End Do
      Else
        Do l = 1, ngrid
          Read (10) ftau(1), ftau(2), ftau(3), ftau(4)
          Read (10)((((fx(i,j,k,m),i=1,nd(l,1)),j=1,nd(l,2)),k=1,nd(l,3)),m=1,ndim+2)
          x(l,:,:,:,:) = dble(fx(:,:,:,:))
          tau = dble(ftau)
        End Do
      End If

    Else
      If (prec(1:1) .eq. 'd') Then
        Do l = 1, ngrid
          Read (10) tau(1), tau(2), tau(3), tau(4)
          Do k = 1, nd(l,3)
            Read (10) (((x(l,i,j,k,m),i=1,nd(l,1)),j=1,nd(l,2)),m=1,ndim+2)
          End Do
        End Do
      Else
        Do l = 1, ngrid
          Read (10) ftau(1), ftau(2), ftau(3), ftau(4)
          Do k = 1, nd(l,3)
            Read (10) (((fx(i,j,k,m),i=1,nd(l,1)),j=1,nd(l,2)),m=1,ndim+2)
          End Do
          x(l,:,:,:,:) = dble(fx(:,:,:,:))
          tau = dble(ftau)
        End Do
      End If
    End If

    Close (10)

    If (prec(1:1) .eq. 's')  Deallocate(fx)

    Return
  End Subroutine p3d_read_soln

  Subroutine p3d_read_grid(NDIM, ngrid, ND, X, IBLANK, prec, gf, vf, ib, file,OF)

    Implicit None

    Real(KIND=8), Dimension(:,:,:,:,:), Pointer :: X
    Integer, Dimension(:,:,:,:), Pointer :: IBLANK
    Character(LEN=2) :: prec, gf, vf, ib
    Character(LEN=80) :: file
    Integer :: ngrid, I, J, K, L, Nx, Ny, Nz, NDIM, M, ftype
    Integer, Dimension(:,:), Pointer :: ND
    Integer :: OMap(3), NMap(3), ND_copy(3)
    Real(KIND=4), Dimension(:,:,:,:), Pointer :: fX
    Logical :: gf_exists
    Integer :: OF

    Inquire(file=trim(file),exist=gf_exists)
    If (.NOT.gf_exists) Then
      Write (*,'(A,A,A)') 'File ', file(1:len_trim(file)), ' is not readable'
      stop
    End If

    prec = ""; gf = ""; vf = ""; ib = ""
    Call p3d_detect(len(trim(file)),trim(file),prec,ndim,gf,vf,ib,ftype)

    Open (unit=10, file=trim(file), form='unformatted', status='old')

    !  Read number of grids
    If (gf(1:1) .eq. 's') Then
      ngrid = 1
    Else
      Read (10) ngrid
    End If
    Allocate(nd(ngrid,3))

    !  Read grid sizes
    nd(:,3) = 1
    Read (10) ((nd(i,j),j=1,ndim),i=1,ngrid)

    !  Allocate grid memory
    nx = maxval(nd(:,1))
    ny = maxval(nd(:,2))
    nz = maxval(nd(:,3))
    Allocate(x(ngrid,nx,ny,nz,ndim))
    If (prec(1:1) .eq. 's') Allocate(fx(nx,ny,nz,ndim))
    Allocate(iblank(ngrid,nx,ny,nz))

    !  Read grid values
    If (vf(1:1) .eq. 'w') Then
      If (prec(1:1) .eq. 'd') Then
        If (ib(1:1) .eq. 'n') Then
          Do l = 1, ngrid
            Read (10)((((x(l,i,j,k,m),i=1,nd(l,1)),j=1,nd(l,2)),k=1,nd(l,3)),m=1,ndim)
          End Do
        Else
          Do l = 1, ngrid
            Read (10)((((x(l,i,j,k,m),i=1,nd(l,1)),j=1,nd(l,2)),k=1,nd(l,3)),m=1,ndim), &
                 (((iblank(l,i,j,k),i=1,nd(l,1)),j=1,nd(l,2)),k=1,nd(l,3))
          End Do
        End If
      Else
        If (ib(1:1) .eq. 'n') Then
          Do l = 1, ngrid
            Read (10)((((fx(i,j,k,m),i=1,nd(l,1)),j=1,nd(l,2)),k=1,nd(l,3)),m=1,ndim)
            x(l,:,:,:,:) = dble(fx(:,:,:,:))
          End Do
        Else
          Do l = 1, ngrid
            Read (10)((((fx(i,j,k,m),i=1,nd(l,1)),j=1,nd(l,2)),k=1,nd(l,3)),m=1,ndim), &
                 (((iblank(l,i,j,k),i=1,nd(l,1)),j=1,nd(l,2)),k=1,nd(l,3))
            x(l,:,:,:,:) = dble(fx(:,:,:,:))
          End Do
        End If
      End If

    Else
      If (ndim .EQ. 3) Then
        If (prec(1:1) .eq. 'd') Then
          If (ib(1:1) .eq. 'n') Then
            Do l = 1, ngrid
              Do k = 1, nd(l,3)
                Read (10) ((x(l,i,j,k,1),i=1,nd(l,1)),j=1,nd(l,2)), &
                     ((x(l,i,j,k,2),i=1,nd(l,1)),j=1,nd(l,2)), &
                     ((x(l,i,j,k,3),i=1,nd(l,1)),j=1,nd(l,2))
              End Do
            End Do
          Else
            Do l = 1, ngrid
              Do k = 1, nd(l,3)
                Read (10) ((x(l,i,j,k,1),i=1,nd(l,1)),j=1,nd(l,2)), &
                     ((x(l,i,j,k,2),i=1,nd(l,1)),j=1,nd(l,2)), &
                     ((x(l,i,j,k,3),i=1,nd(l,1)),j=1,nd(l,2)), &
                     ((iblank(l,i,j,k),i=1,nd(l,1)),j=1,nd(l,2))
              End Do
            End Do
          End If
        Else
          If (ib(1:1) .eq. 'n') Then
            Do l = 1, ngrid
              Do k = 1, nd(l,3)
                Read (10) ((fx(i,j,k,1),i=1,nd(l,1)),j=1,nd(l,2)), &
                     ((fx(i,j,k,2),i=1,nd(l,1)),j=1,nd(l,2)), &
                     ((fx(i,j,k,3),i=1,nd(l,1)),j=1,nd(l,2))
              End Do
              x(l,:,:,:,:) = dble(fx(:,:,:,:))
            End Do
          Else
            Do l = 1, ngrid
              Do k = 1, nd(l,3)
                Read (10) ((fx(i,j,k,1),i=1,nd(l,1)),j=1,nd(l,2)), &
                     ((fx(i,j,k,2),i=1,nd(l,1)),j=1,nd(l,2)), &
                     ((fx(i,j,k,3),i=1,nd(l,1)),j=1,nd(l,2)), &
                     ((iblank(l,i,j,k),i=1,nd(l,1)),j=1,nd(l,2))
              End Do
              x(l,:,:,:,:) = dble(fx(:,:,:,:))
            End Do
          End If
        End If
      Else
        If (prec(1:1) .eq. 'd') Then
          If (ib(1:1) .eq. 'n') Then
            Do l = 1, ngrid
              Do j = 1, nd(l,2)
                Read (10) (x(l,i,j,1,1),i=1,nd(l,1)), &
                     (x(l,i,j,1,2),i=1,nd(l,1))
              End Do
            End Do
          Else
            Do l = 1, ngrid
              Do j = 1, nd(l,2)
                Read (10) (x(l,i,j,1,1),i=1,nd(l,1)), &
                     (x(l,i,j,1,2),i=1,nd(l,1)), &
                     (iblank(l,i,j,1),i=1,nd(l,1))
              End Do
            End Do
          End If
        Else
          If (ib(1:1) .eq. 'n') Then
            Do l = 1, ngrid
              Do k = j, nd(l,2)
                Read (10) (fx(i,j,1,1),i=1,nd(l,1)), &
                     (fx(i,j,1,2),i=1,nd(l,1))
              End Do
              x(l,:,:,:,:) = dble(fx(:,:,:,:))
            End Do
          Else
            Do l = 1, ngrid
              Do k = j, nd(l,2)
                Read (10) (fx(i,j,1,1),i=1,nd(l,1)), &
                     (fx(i,j,1,2),i=1,nd(l,1)), &
                     (iblank(l,i,j,1),i=1,nd(l,1))
              End Do
              x(l,:,:,:,:) = dble(fx(:,:,:,:))
            End Do
          End If
        End If
      End If
    End If

    Close (10)

    If (prec(1:1) .eq. 's')  Deallocate(fx)
    If (ib(1:1) .eq. 'n') iblank = 1

    Return
  End Subroutine p3d_read_grid