!!****f* source/Simulation/SimulationMain/tomas/Simulation_initBlock
!!
!! NAME
!!
!!  Simulation_initBlock
!!
!!
!! SYNOPSIS
!!
!!  Simulation_initBlock(integer(IN) :: blockId)
!!
!!
!!
!! DESCRIPTION
!!  This routine applies initial conditions of a specific simulation
!!  to the specified block.
!!
!! 
!! ARGUMENTS
!!
!!  blockId -         the number of the block to update
!!
!! PARAMETERS
!!
!!
!! SEE ALSO
!!
!!  Eos_wrapped
!!***

subroutine Simulation_initBlock(blockId)
  use Simulation_data, ONLY : sim_dens, sim_temp, sim_energy, sim_Rexp
  use Grid_interface, ONLY : Grid_getBlkIndexLimits, Grid_getCellCoords, &
    Grid_getDeltas , Grid_getBlkPtr, Grid_releaseBlkPtr

  use PhysicalConstants_interface, ONLY : PhysicalConstants_get
  use Eos_interface, ONLY : Eos_wrapped

    
  ! to get definitions of macros LOW, HIGH, IAXIS, etc.
#include "constants.h"
#include "Flash.h"
  ! INPUT is the block number
  integer,intent(IN) ::  blockId
  
  ! VARIABLE definitions
  real,allocatable,dimension(:) :: xCoord,yCoord,zCoord
  integer,dimension(2,MDIM)     :: blkLimits,blkLimitsGC
  integer :: sizeX,sizeY,sizeZ
  real    :: del(MDIM)
  logical :: gcell = .true.
  integer :: ii, jj, kk
  real, POINTER, DIMENSION(:,:,:,:) :: solnData
  real    :: thigh, boltz, mH, pi, r2, mu = 0.609


  call PhysicalConstants_get('Boltzmann', boltz)
  call PhysicalConstants_get('proton mass', mH)
  call PhysicalConstants_get('pi', pi)

  gam = 1.66667
  thigh =  3 * (gam-1.0) * mu * mH * sim_energy &
        / (4 * pi * sim_dens * boltz * sim_Rexp**3)

  print *, thigh


  ! get the coordinate information for the current block from the database
  ! and allocate coordinate arrays
  call Grid_getBlkIndexLimits(blockId,blkLimits,blkLimitsGC)

  sizeX = blkLimitsGC(HIGH,IAXIS) - blkLimitsGC(LOW,IAXIS) + 1
  allocate(xCoord(sizeX),stat=istat)
  sizeY = blkLimitsGC(HIGH,JAXIS) - blkLimitsGC(LOW,JAXIS) + 1
  allocate(yCoord(sizeY),stat=istat)
  sizeZ = blkLimitsGC(HIGH,KAXIS) - blkLimitsGC(LOW,KAXIS) + 1
  allocate(zCoord(sizeZ),stat=istat)
  
  ! fill the coordinate arrays with values
  if (NDIM == 3) call Grid_getCellCoords&
                      (KAXIS, blockId, CENTER, gcell, zCoord, sizeZ)
  if (NDIM >= 2) call Grid_getCellCoords&
                      (JAXIS, blockId, CENTER,gcell, yCoord, sizeY)
  call Grid_getCellCoords(IAXIS, blockId, CENTER, gcell, xCoord, sizeX)
  
  ! eventually get the grid cell size
  call Grid_getDeltas(blockId,del)
  call Grid_getBlkPtr(blockId, solnData, CENTER)
  do ii = blkLimits(LOW, IAXIS), blkLimits(HIGH, IAXIS)
    do jj = blkLimits(LOW, JAXIS), blkLimits(HIGH, JAXIS)
      do kk = blkLimits(LOW, KAXIS), blkLimits(HIGH, KAXIS)

    if (xCoord(ii) .gt. 1.1e+19) then
            solnData(DENS_VAR, ii, jj, kk) = sim_dens*1000
    elseif (xCoord(ii) .gt. 0.9e+19) then
            solnData(DENS_VAR, ii, jj, kk) = sim_dens*100
    elseif (xCoord(ii) .gt. 0.7e+19) then
            solnData(DENS_VAR, ii, jj, kk) = sim_dens*10
    else
        solnData(DENS_VAR, ii, jj, kk) = sim_dens
    endif

        solnData(VELX_VAR, ii, jj, kk) = 0.0
        solnData(VELY_VAR, ii, jj, kk) = 0.0
        solnData(VELZ_VAR, ii, jj, kk) = 0.0
        r2 = xCoord(ii)**2 + yCoord(jj)**2 + zCoord(kk)**2
        if (r2 .lt. sim_Rexp**2) then
          solnData(TEMP_VAR, ii, jj, kk) = thigh
        else
          solnData(TEMP_VAR, ii, jj, kk) = sim_temp
        endif
        solnData(EINT_VAR, ii, jj, kk) = solnData(DENS_VAR, ii, jj, kk) * boltz * solnData(TEMP_VAR, ii, jj, kk) &
                                         / (mu * mH * (gam-1))
        solnData(ENER_VAR, ii, jj, kk) = solnData(DENS_VAR, ii, jj, kk) * boltz * solnData(TEMP_VAR, ii, jj, kk) &
                                         / (mu * mH * (gam-1))
        solnData(PRES_VAR, ii, jj, kk) = solnData(DENS_VAR, ii, jj, kk) * boltz * solnData(TEMP_VAR, ii, jj, kk) &
                                         / (mu * mH)
        !print *, solnData(:,ii,jj,kk)

        solnData(IHP_SPEC, ii, jj, kk) = 0.0
        solnData(IHA_SPEC, ii, jj, kk) = 1.0
        solnData(IHP_SPEC, ii, jj, kk) = 0.0
        solnData(ICO_SPEC, ii, jj, kk) = 0.0
        solnData(ICP_SPEC, ii, jj, kk) = 1.0

      enddo
    enddo
  enddo

  call Eos_wrapped(MODE_DENS_TEMP, blkLimits, blockID)
  call Grid_releaseBlkPtr(blockID,solnData, CENTER)


  deallocate(xCoord)
  deallocate(yCoord)
  deallocate(zCoord)



end subroutine Simulation_initBlock