MODULE initial_conditions

  USE shared_data
  USE neutral

  IMPLICIT NONE

  PRIVATE

  PUBLIC:: set_initial_conditions

CONTAINS

  !---------------------------------------------------------------------------
  ! This function sets up the initial condition for the code
  ! The variables which must be set are
  ! Rho - density
  ! V{x, y, z} - Velocities in x, y, z
  ! B{x, y, z} - Magnetic fields in x, y, z
  ! Energy - Specific internal energy
  ! Since temperature is a more intuitive quantity than specific internal energy
  ! There is a helper function get_energy which converts temperature to energy
  ! The syntax for this function (which is in core/neutral.f90) is
  !
  ! CALL get_energy(density, temperature, equation_of_state, &
  !     output_energy) 
  !
  ! REAL(num) :: density - The density at point (ix, iy) on the grid
  ! REAL(num) :: temperature - The temperature at point (ix, iy) on the grid
  ! INTEGER :: equation_of_state - The code for the equation of state to use.
  !            The global equation of state for the code is eos_number
  ! REAL(num) :: output_energy - The specific internal energy returned by
  !              the routine
  !
  ! You may also need the neutral fraction. This can be calculated by a function
  ! call to  get_neutral(temperature, rho, z). This routine is in core/neutral.f90
  ! and requires the local temperature and mass density. For example to set
  ! xi_n to the neutral fraction use
  ! xi_n = get_neutral(temperature, rho, z)
  !---------------------------------------------------------------------------
  
  SUBROUTINE set_initial_conditions
 INTEGER :: ix, iy,iz
    REAL(num) ::  beta,a1,b1,c1,a,zs,hw

    vx = 0.0_num
    vy = 0.0_num
    vz = 0.0_num
    bx = 0.0_num
    by = 0.0_num
    grav=0.0_num
 !   rho=1.0_num
    bz=1.0_num
    
    a = 0.01_num
    zs = 5.0_num
    hw = 7.5_num
   
    
    beta=0.0_num
    
DO ix = -1, nx+2
DO iy = -1, ny+2
DO iz = -1, nz+2
rho(ix,iy,iz)=(1.0_num+15.0_num*exp(-((xc(ix)-5.0_num*pi)/(0.5_num*hw))**6 ) )
END DO  
END DO
END DO



DO ix = -1, nx+2
DO iy = -2, ny+2
DO iz = -1, nz+2
by(ix,iy,iz)=a*exp(-((zb(iz)-zs)/2.5_num)**2)
END DO  
END DO
END DO


 

DO ix = -2, nx+2
DO iy = -2, ny+2
DO iz = -2, nz+2
vy(ix,iy,iz)=-a*exp(-((zb(iz)-zs)/2.5_num)**2)/sqrt(1.0_num+15.0_num*exp(-((xb(ix)-5.0_num*pi)/(0.5_num*hw))**6 ) )
END DO  
END DO
END DO

energy=0.5_num*(beta*1.0_num) / (rho * (gamma-1.0_num))

END SUBROUTINE set_initial_conditions


END MODULE initial_conditions