program inj

      use avbp_constants
      use avbp_solution
      use avbp_thermo
      use avbp_coor
      use avbp_rundat

      implicit none

      integer i, ix, iy, k,
     +        choice
      double precision xx, yy, ut, vt,
     +                 gmm1, r, eps,
     +                 u, v, P, rho, T,
     +                 ec, e,rhoe,
     +                 YF, YO

      eps  = 1.0d-5

c -------------------------------------------------
c Reading data input file
c -------------------------------------------------
      open(1,file='init_inj.choices')
        read(1,*) solutfile
        read(1,*) coorfile
        read(1,*) choice
        read(1,*) P
        read(1,*) T
        read(1,*) u
        read(1,*) v
      close(1)

c -------------------------------------------------
c Reading AVBP input files
c -------------------------------------------------
      call readrundat
      call readthermo
      call readcoor

c -------------------------------------------------
c Allocation of memory
c -------------------------------------------------
      sol_niter    = 0
      sol_nnode    = coor_nnode
      sol_neq      = coor_neq
      sol_neqs     = 2
      sol_neqt     = 0
      sol_ndim     = coor_ndim
      sol_nreac    = 0
      sol_nadd     = 0
      sol_neq2pf   = 0
      sol_iles     = rundat_iles
      sol_ichem    = rundat_ichem
      sol_dtsum    = 0.0d0
      sol_cp       = thermo_cp
      sol_fmach    = sol_cp

      call createsolutionmemory

c -------------------------------------------------
c Beginnig initialisation
c -------------------------------------------------
      r    = rgas / thermo_mean_mole_mass
      gmm1 = (thermo_cp / (thermo_cp-r)) - 1

      sol_gmm1     = gmm1

      rho = P/(r*T)

      ix = 1
      iy = 2

c -------------------------------------------------
c First choice
c -------------------------------------------------
      if (choice.eq.0) then
        ut = 0.0d0
        vt = 0.0d0

        do i=1, sol_nnode
          xx = coor((ix-1)*sol_nnode+i)
          yy = coor((iy-1)*sol_nnode+i)

          rhoe = p/gmm1+half*rho*(ut*ut+vt*vt)

c NS variables
          k=1
          w((k-1)*sol_nnode+i) = rho
          k=2
          w((k-1)*sol_nnode+i) = rho * ut
          k=3
          w((k-1)*sol_nnode+i) = rho * vt
          k=4
          w((k-1)*sol_nnode+i) = rhoe

c Mass fractions
            if (dabs(yy).ge.0.055) then

            k=1
            w_mass((k-1)*sol_nnode+i) = rho * 1.0d0
            k=2
            w_mass((k-1)*sol_nnode+i) = 0.0d0

            else

            k=1
            w_mass((k-1)*sol_nnode+i) = 0.0d0
            k=2
            w_mass((k-1)*sol_nnode+i) = rho * 1.0d0
 

            endif
        enddo
 

c -------------------------------------------------
c Second choice
c -------------------------------------------------
      elseif (choice.eq.1) then

        do i=1, sol_nnode
          xx = coor((ix-1)*sol_nnode+i)
          yy = coor((iy-1)*sol_nnode+i)
c NS variables

c Mass fractions

        enddo

      endif

      winf(1)=rho
      winf(2)=rho
      winf(3)=rho
      winf(4)=rho

c -------------------------------------------------
c Writing solution... and that's it !!!
c -------------------------------------------------
      call writesolution
      call writesolutioninfo
      call destroysolutionmemory
 

      end