module mesh
  implicit none

  !== domain size ==!

  real(8), parameter :: xmin = -0.5d0
  real(8), parameter :: xmax = +0.5d0
  real(8), parameter :: ymin =  0.0d0
  real(8), parameter :: ymax =  1.0d0
  real(8), parameter :: zmin =  0.0d0
  real(8), parameter :: zmax =  1.0d0

  !=== total resolution ===!

  integer, parameter :: nxg = 32
  integer, parameter :: nyg = 32
  integer, parameter :: nzg = 32

  !== information on local mesh ===!

  ! the number of conservative variables
  integer, parameter :: nd  = 8

  ! ghost cell
  integer, parameter :: ng  = 8
  integer, parameter :: ng2 = ng * 2

  ! domain size and cell size
  real(8), parameter :: Lx = xmax - xmin
  real(8), parameter :: Ly = ymax - ymin
  real(8), parameter :: Lz = zmax - zmin
  real(8), parameter :: dx = Lx / dble(nxg)
  real(8), parameter :: dy = Ly / dble(nyg)
  real(8), parameter :: dz = Lz / dble(nzg)

  integer :: nx, ny, nz
  
  ! coordinate systems
  integer, allocatable :: dnx(:), dny(:), dnz(:)
  real(8), allocatable :: xc(:), xb(:)
  real(8), allocatable :: yc(:), yb(:)
  real(8), allocatable :: zc(:), zb(:)

contains


  subroutine init_mesh
    use parallel
    integer :: i, j, k, ierr

    !== local mesh numbers ==!

    allocate(dnx(cart3d%dims(1)), dny(cart3d%dims(2)), dnz(cart3d%dims(3)))

    dnx(:) = nxg / cart3d%dims(1)
    dny(:) = nyg / cart3d%dims(2)
    dnz(:) = nzg / cart3d%dims(3)

    dnx(1:mod(nxg,cart3d%dims(1))) = dnx(1:mod(nxg,cart3d%dims(1))) + 1
    dny(1:mod(nyg,cart3d%dims(2))) = dny(1:mod(nyg,cart3d%dims(2))) + 1
    dnz(1:mod(nzg,cart3d%dims(3))) = dnz(1:mod(nzg,cart3d%dims(3))) + 1

    nx = dnx(cart3d%coords(1)+1) + ng2
    ny = dny(cart3d%coords(2)+1) + ng2
    nz = dnz(cart3d%coords(3)+1) + ng2

    allocate (xb(0:nx), yb(0:ny), zb(0:nz))
    allocate (xc(1:nx), yc(1:ny), zc(1:nz))

    !== x-coordinate system ==!

    do i = 0, nx
       xb(i) = xmin + dx * (dble(i-ng) + dble(sum(dnx(1:cart3d%coords(1)))))
    end do
    do i = 1, nx
       xc(i) = 0.5d0 * (xb(i) + xb(i-1))
    end do

    !== y-coordinate system ==!
    do j = 0, ny
       yb(j) = ymin + dy * (dble(j-ng) + dble(sum(dny(1:cart3d%coords(2)))))
    end do
    do j = 1, ny
       yc(j) = 0.5d0 * (yb(j) + yb(j-1))
    end do

    !== z-coorinate system ==!
    do k = 0, nz
       zb(k) = zmin + dz * (dble(k-ng) + dble(sum(dnz(1:cart3d%coords(3)))))
    end do
    do k = 1, nz
       zc(k) = 0.5d0 * (zb(k) + zb(k-1))
    end do

    return
  end subroutine init_mesh


end module mesh