!!****f* setups/jet_injection/init_block
!!
!! NAME
!!
!!  init_block
!!
!!
!! SYNOPSIS
!!
!!  init_block(block_no)
!!
!!  init_block(integer)
!!
!!
!! DESCRIPTION
!!
!!  
!!
!!  Initializes the fluid data (density, pressure, velocity, etc.)  for the
!!  problem of jet injection through a nozzle in the left X boundary of the
!!  computational domain. 
!!
!!
!!  By Michal Hanasz (MH), based on the windtunnel 'init_block.F90'
!!
!!  Initializes fluid data (density, pressure, velocity, etc.) for
!!  a specified block.  
!!
!!
!!
!! ARGUMENTS
!!
!!  block_no        the number of the block to initialize
!!
!!
!! PARAMETERS
!!
!!  p_a       Initial ambient pressure
!!
!!  p_j       Jet pressure
!!
!!  rho_a       Initial ambient density
!!
!!  rho_j       Jet density
!!
!!  jet_vel        Inflow velocity
!!
!!  rj         Jet radius
!!
!!  mjb        Steepness parameter of the jet boundary 
!!             See user_bnd for the following lines:
!!
!!          solnData(idens,i,j,k) = rho_a - (rho_a-rho_j)/cosh((y(j)/rj)**mjb)
!!          solnData(ivelx,i,j,k) = jet_vel/cosh((y(j)/rj)**mjb)
!!          solnData(ipres,i,j,k) = p_a - (p_a-p_j)/cosh((y(j)/rj)**mjb)
!!
!!
!!
!!  gamma           the adiabatic index
!!
!!  smallp          the low pressure cut-off value
!!
!!***

subroutine init_block (block_no)

!===============================================================================

  use logfile, ONLY: stamp_logfile

  use runtime_parameters, ONLY: get_parm_from_context, GLOBAL_PARM_CONTEXT

  use dBase, ONLY: nguard, nxb, nyb, nzb, k2d, k3d, ionmax, &
       dBasePropertyInteger, &
       dBaseKeyNumber, &
       dBaseSpecies, &
       dBasePutData

  implicit none

  integer :: block_no
  integer :: i, j, k, n

  integer, save :: MyPE, MasterPE
  
  integer, save :: idens, iener, ivelx, ively, ivelz, ipres, igamc, igame, &
       inuc_begin

  integer, save :: iPt

  real :: rho_zone, velx_zone, vely_zone, velz_zone, pres_zone, &
       ener_zone, ekin_zone

  real, save :: gamma
  real, save :: smallp, smallx

! Jet parameters (MH) ------------

      integer, save :: mjb
      real, save    :: p_a, p_j, rho_j, rho_a, rj, jet_vel

!---------------------------------

  logical, save :: firstCall = .TRUE.


!===============================================================================

  if (firstCall) then

     MyPE = dBasePropertyInteger('MyProcessor')
     MasterPE = dBasePropertyInteger('MasterProcessor')

     idens = dBaseKeyNumber('dens')
     iener = dBaseKeyNumber('ener')

     ivelx = dBaseKeyNumber('velx')
     ively = dBaseKeyNumber('vely')
     ivelz = dBaseKeyNumber('velz')

     ipres = dBaseKeyNumber('pres')
     igamc = dBaseKeyNumber('gamc')
     igame = dBaseKeyNumber('game')

     inuc_begin = dBaseSpecies(1)

     iPt = dBaseKeyNumber('Point')

     call get_parm_from_context(GLOBAL_PARM_CONTEXT, 'gamma', gamma)

     call get_parm_from_context(GLOBAL_PARM_CONTEXT, 'smallp', smallp)
     call get_parm_from_context(GLOBAL_PARM_CONTEXT, 'smallx', smallx)

! (MH)
     call get_parm_from_context(GLOBAL_PARM_CONTEXT, 'p_a', p_a)
     call get_parm_from_context(GLOBAL_PARM_CONTEXT, 'rho_a', rho_a)
     call get_parm_from_context(GLOBAL_PARM_CONTEXT, 'rho_j', rho_j)
     call get_parm_from_context(GLOBAL_PARM_CONTEXT, 'jet_vel', jet_vel)
     call get_parm_from_context(GLOBAL_PARM_CONTEXT, 'rj', rj)
     call get_parm_from_context(GLOBAL_PARM_CONTEXT, 'mjb', mjb)
! ----

     if (MyPE == MasterPE) then

        call stamp_logfile("initializing for JET INJECTION", 'run_init')
        write (*,*) "flash:  initializing for JET INJECTION"

     endif

     firstCall = .FALSE.

  endif



! Loop over cells in the block.  In this problem the initial
! conditions are spatially uniform.

  do k = 1, 2*nguard*k3d + nzb
     do j = 1, 2*nguard*k2d + nyb
        do i = 1, 2*nguard + nxb
           
           rho_zone = rho_a 
           pres_zone = p_a 

           velx_zone = 0.0   
           vely_zone = 0.0
           velz_zone = 0.0


! Compute the gas energy and set the gamma-values needed for
! the equation of state.
           ekin_zone = 0.5 * (velx_zone**2 + vely_zone**2 + velz_zone**2)

           ener_zone = pres_zone / (gamma-1.)
           ener_zone = ener_zone / rho_zone
           ener_zone = ener_zone + ekin_zone
           ener_zone = max(ener_zone, smallp)


           call dBasePutData(inuc_begin,iPt, i,j,k, block_no, &
     &                       1.0-(ionmax-1)*smallx, 1)
           do n = 2, ionmax
              call dBasePutData(inuc_begin-1+n,iPt, i,j,k, block_no, smallx, 1)
           enddo


! store the variables in the current zone via the database put methods
           call dBasePutData(idens, iPt, i,j,k, block_no, rho_zone, 1)
           call dBasePutData(ipres, iPt, i,j,k, block_no, pres_zone, 1)
           call dBasePutData(iener, iPt, i,j,k, block_no, ener_zone, 1)
           
           call dBasePutData(igamc, iPt, i,j,k, block_no, gamma, 1)
           call dBasePutData(igame, iPt, i,j,k, block_no, gamma, 1)
           
           call dBasePutData(ivelx, iPt, i,j,k, block_no, velx_zone, 1)
           call dBasePutData(ively, iPt, i,j,k, block_no, vely_zone, 1)
           call dBasePutData(ivelz, iPt, i,j,k, block_no, velz_zone, 1)

 
        enddo
     enddo
  enddo



  return
end subroutine init_block