bc_to_gc Subroutine

  subroutine bc_to_gc(mg, id, nc, iv, nb, bc_type)
    type(mg_t), intent(inout) :: mg
    integer, intent(in)       :: id
    integer, intent(in)       :: nc
    integer, intent(in)       :: iv
    integer, intent(in)       :: nb
    !< Neighbor direction
    integer, intent(in)       :: bc_type
    !< Type of b.c.
    real(dp)                  :: c0, c1, c2, dr

    ! If we call the interior point x1, x2 and the ghost point x0, then a
    ! Dirichlet boundary value b can be imposed as:
    ! x0 = -x1 + 2*b
    ! A Neumann b.c. can be imposed as:
    ! x0 = x1 +/- dx * b
    ! A continuous boundary (same slope) as:
    ! x0 = 2 * x1 - x2
    ! Below, we set coefficients to handle these cases
    select case (bc_type)
    case (mg_bc_dirichlet)
       c0 = 2
       c1 = -1
       c2 = 0
    case (mg_bc_neumann)
       dr = mg%dr(mg_neighb_dim(nb), mg%boxes(id)%lvl)
       c0 = dr * mg_neighb_high_pm(nb) ! This gives a + or - sign
       c1 = 1
       c2 = 0
    case (mg_bc_continuous)
       c0 = 0
       c1 = 2
       c2 = -1
    case default
       error stop "bc_to_gc: unknown boundary condition"
    end select

    select case (nb)
    case (mg_neighb_lowx)
       mg%boxes(id)%cc(0, 1:nc, iv) = c0 * mg%boxes(id)%cc(0, 1:nc,&
           iv) + c1 * mg%boxes(id)%cc(1, 1:nc, iv) + c2 * mg%boxes(id)%cc(2,&
           1:nc, iv)
    case (mg_neighb_highx)
       mg%boxes(id)%cc(nc+1, 1:nc, iv) = c0 * mg%boxes(id)%cc(nc+1, 1:nc,&
           iv) + c1 * mg%boxes(id)%cc(nc, 1:nc,&
           iv) + c2 * mg%boxes(id)%cc(nc-1, 1:nc, iv)
    case (mg_neighb_lowy)
       mg%boxes(id)%cc(1:nc, 0, iv) = c0 * mg%boxes(id)%cc(1:nc, 0,&
           iv) + c1 * mg%boxes(id)%cc(1:nc, 1, iv) + c2 * mg%boxes(id)%cc(1:nc,&
           2, iv)
    case (mg_neighb_highy)
       mg%boxes(id)%cc(1:nc, nc+1, iv) = c0 * mg%boxes(id)%cc(1:nc, nc+1,&
           iv) + c1 * mg%boxes(id)%cc(1:nc, nc,&
           iv) + c2 * mg%boxes(id)%cc(1:nc, nc-1, iv)
    end select
  end subroutine bc_to_gc