diff --git a/sources/boundaries.F90 b/sources/boundaries.F90
index 649db61..82c7708 100644
--- a/sources/boundaries.F90
+++ b/sources/boundaries.F90
@@ -45,6 +45,32 @@ module boundaries
       real(kind=8), dimension(:)      , intent(in)    :: x, y, z
       real(kind=8), dimension(:,:,:,:), intent(inout) :: qn
     end subroutine
+
+#if NDIMS == 3
+    subroutine block_face_prolong_iface(dir, pos, qn, qb, status)
+      implicit none
+      integer                         , intent(in)    :: dir
+      integer     , dimension(3)      , intent(in)    :: pos
+      real(kind=8), dimension(:,:,:,:), intent(in)    :: qn
+      real(kind=8), dimension(:,:,:,:), intent(inout) :: qb
+      integer                         , intent(out)   :: status
+    end subroutine
+#endif /* NDIMS == 3 */
+    subroutine block_edge_prolong_iface(dir, pos, qn, qb, status)
+      implicit none
+      integer                         , intent(in)    :: dir
+      integer     , dimension(3)      , intent(in)    :: pos
+      real(kind=8), dimension(:,:,:,:), intent(in)    :: qn
+      real(kind=8), dimension(:,:,:,:), intent(inout) :: qb
+      integer                         , intent(out)   :: status
+    end subroutine
+    subroutine block_corner_prolong_iface(pos, qn, qb, status)
+      implicit none
+      integer     , dimension(3)      , intent(in)    :: pos
+      real(kind=8), dimension(:,:,:,:), intent(in)    :: qn
+      real(kind=8), dimension(:,:,:,:), intent(inout) :: qb
+      integer                         , intent(out)   :: status
+    end subroutine
   end interface
 
 ! pointers to custom boundary conditions, one per direction handling both sides
@@ -53,6 +79,14 @@ module boundaries
                                              custom_boundary_x => null(), &
                                              custom_boundary_y => null(), &
                                              custom_boundary_z => null()
+#if NDIMS == 3
+  procedure(block_face_prolong_iface), pointer, save :: &
+                               block_face_prolong => block_face_prolong_tvd
+#endif /* NDIMS == 3 */
+  procedure(block_edge_prolong_iface), pointer, save :: &
+                               block_edge_prolong => block_edge_prolong_tvd
+  procedure(block_corner_prolong_iface), pointer, save :: &
+                               block_corner_prolong => block_corner_prolong_tvd
 
 ! supported boundary types
 !
@@ -108,12 +142,13 @@ module boundaries
 !
   subroutine initialize_boundaries(status)
 
-    use coordinates, only : periodic
-    use helpers    , only : print_message
+    use coordinates   , only : periodic
+    use helpers       , only : print_message
+    use interpolations, only : prolongation_method
 #ifdef MPI
-    use mpitools   , only : npmax
+    use mpitools      , only : npmax
 #endif /* MPI */
-    use parameters , only : get_parameter
+    use parameters    , only : get_parameter
 
     implicit none
 
@@ -129,6 +164,21 @@ module boundaries
 !
     status = 0
 
+    select case(trim(prolongation_method))
+    case ("ppm", "PPM")
+#if NDIMS == 3
+      block_face_prolong   => block_face_prolong_ppm
+#endif /* NDIMS == 3 */
+      block_edge_prolong   => block_edge_prolong_ppm
+      block_corner_prolong => block_corner_prolong_ppm
+    case default
+#if NDIMS == 3
+      block_face_prolong   => block_face_prolong_tvd
+#endif /* NDIMS == 3 */
+      block_edge_prolong   => block_edge_prolong_tvd
+      block_corner_prolong => block_corner_prolong_tvd
+    end select
+
     do n = 1, NDIMS
       do s = 1, 2
 
@@ -3793,11 +3843,12 @@ module boundaries
 !
 !===============================================================================
 !
-! subroutine BLOCK_FACE_PROLONG:
-! -----------------------------
+! subroutine BLOCK_FACE_PROLONG_TVD:
+! ---------------------------------
 !
-!   Subroutine prolongs the region of qn specified by the face position
-!   and inserts it in the corresponding face boundary region of qb.
+!   This subroutine uses Total Variation Diminishing Method (TVD) to prolong
+!   a region of qn, specified by the face position, and inserts it into
+!   the corresponding face boundary region of qb.
 !
 !   Arguments:
 !
@@ -3809,7 +3860,7 @@ module boundaries
 !
 !===============================================================================
 !
-  subroutine block_face_prolong(dir, pos, qn, qb, status)
+  subroutine block_face_prolong_tvd(dir, pos, qn, qb, status)
 
     use coordinates   , only : nb, ne, ncells_half, nghosts_half
     use equations     , only : nv, positive
@@ -3835,7 +3886,7 @@ module boundaries
     logical                , save :: first = .true.
     integer, dimension(2,2), save :: nlims, tlims
 
-    character(len=*), parameter :: loc = 'BOUNDARIES::block_face_prolong()'
+    character(len=*), parameter :: loc = 'BOUNDARIES::block_face_prolong_tvd()'
 
 !-------------------------------------------------------------------------------
 !
@@ -3943,7 +3994,190 @@ module boundaries
 
 !-------------------------------------------------------------------------------
 !
-  end subroutine block_face_prolong
+  end subroutine block_face_prolong_tvd
+!
+!===============================================================================
+!
+! subroutine BLOCK_FACE_PROLONG_PPM:
+! ---------------------------------
+!
+!   This subroutine uses Colella and Woodward's Piecewise Parabolic Method (PPM)
+!   to prolong a region of qn, specified by the face position, and inserts it
+!   into the corresponding face boundary region of qb.
+!
+!   Arguments:
+!
+!     dir    - the face direction;
+!     pos    - the face position;
+!     qn     - the array to prolong from;
+!     qb     - the array of the face boundary region;
+!     status - the call status;
+!
+!   References:
+!
+!     [1] Colella, P., Woodward, P. R.,
+!         "The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations",
+!         Journal of Computational Physics, 1984, vol. 54, pp. 174-201,
+!         https://doi.org/10.1016/0021-9991(84)90143-8
+!
+!===============================================================================
+!
+  subroutine block_face_prolong_ppm(dir, pos, qn, qb, status)
+
+    use coordinates   , only : nb, ne, ncells_half, nghosts_half
+    use equations     , only : nv, positive
+    use helpers       , only : print_message
+    use interpolations, only : limiter_prol
+
+    implicit none
+
+    integer                         , intent(in)    :: dir
+    integer     , dimension(3)      , intent(in)    :: pos
+    real(kind=8), dimension(:,:,:,:), intent(in)    :: qn
+    real(kind=8), dimension(:,:,:,:), intent(inout) :: qb
+    integer                         , intent(out)   :: status
+
+    integer      :: p, n
+    integer      :: i, il, iu, is, it, im, ip
+    integer      :: j, jl, ju, js, jt, jm, jp
+    integer      :: k, kl, ku, ks, kt, km, kp
+    real(kind=8) :: qc, dq1, dq2, dq3, dq4
+
+    real(kind=8), dimension(3) :: ql, qr, dq
+
+    logical                , save :: first = .true.
+    integer, dimension(2,2), save :: nlims, tlims
+
+    character(len=*), parameter :: loc = 'BOUNDARIES::block_face_prolong_ppm()'
+
+!-------------------------------------------------------------------------------
+!
+    if (first) then
+      nlims(1,1) = ne - nghosts_half + 1
+      nlims(1,2) = ne
+      nlims(2,1) = nb
+      nlims(2,2) = nb + nghosts_half - 1
+
+      tlims(1,1) = nb
+      tlims(1,2) = nb + ncells_half + nghosts_half - 1
+      tlims(2,1) = ne - ncells_half - nghosts_half + 1
+      tlims(2,2) = ne
+
+      first = .false.
+    end if
+
+    status = 0
+
+    select case(dir)
+    case(1)
+      il = nlims(pos(1),1)
+      iu = nlims(pos(1),2)
+      jl = tlims(pos(2),1)
+      ju = tlims(pos(2),2)
+      kl = tlims(pos(3),1)
+      ku = tlims(pos(3),2)
+    case(2)
+      il = tlims(pos(1),1)
+      iu = tlims(pos(1),2)
+      jl = nlims(pos(2),1)
+      ju = nlims(pos(2),2)
+      kl = tlims(pos(3),1)
+      ku = tlims(pos(3),2)
+    case default
+      il = tlims(pos(1),1)
+      iu = tlims(pos(1),2)
+      jl = tlims(pos(2),1)
+      ju = tlims(pos(2),2)
+      kl = nlims(pos(3),1)
+      ku = nlims(pos(3),2)
+    end select
+
+    do k = kl, ku
+      km = k - 1
+      kp = k + 1
+      ks = 2 * (k - kl) + 1
+      kt = ks + 1
+      do j = jl, ju
+        jm = j - 1
+        jp = j + 1
+        js = 2 * (j - jl) + 1
+        jt = js + 1
+        do i = il, iu
+          im = i - 1
+          ip = i + 1
+          is = 2 * (i - il) + 1
+          it = is + 1
+
+          do p = 1, nv
+
+            qc    = qn(p,i,j,k)
+
+            ql(1) = (7.0d+00 * (qn(p,i-1,j,k) + qn(p,i  ,j,k)) &
+                             - (qn(p,i-2,j,k) + qn(p,i+1,j,k))) / 1.2d+01
+            qr(1) = (7.0d+00 * (qn(p,i  ,j,k) + qn(p,i+1,j,k)) &
+                             - (qn(p,i-1,j,k) + qn(p,i+2,j,k))) / 1.2d+01
+            ql(2) = (7.0d+00 * (qn(p,i,j-1,k) + qn(p,i,j  ,k)) &
+                             - (qn(p,i,j-2,k) + qn(p,i,j+1,k))) / 1.2d+01
+            qr(2) = (7.0d+00 * (qn(p,i,j  ,k) + qn(p,i,j+1,k)) &
+                             - (qn(p,i,j-1,k) + qn(p,i,j+2,k))) / 1.2d+01
+#if NDIMS == 3
+            ql(3) = (7.0d+00 * (qn(p,i,j,k-1) + qn(p,i,j,k  )) &
+                             - (qn(p,i,j,k-2) + qn(p,i,j,k+1))) / 1.2d+01
+            qr(3) = (7.0d+00 * (qn(p,i,j,k  ) + qn(p,i,j,k+1)) &
+                             - (qn(p,i,j,k-1) + qn(p,i,j,k+2))) / 1.2d+01
+#endif /* NDIMS == 3 */
+
+            do n = 1, NDIMS
+              if ((qr(n) - qc) * (ql(n) - qc) > 0.0d+00) then
+                ql(n) = qc
+                qr(n) = qc
+              else
+                if (abs(qr(n) - qc) >= 2.0d+00 * abs(ql(n) - qc)) then
+                  qr(n) = qc - 2.0d+00 * (ql(n) - qc)
+                end if
+                if (abs(ql(n) - qc) >= 2.0d+00 * abs(qr(n) - qc)) then
+                  ql(n) = qc - 2.0d+00 * (qr(n) - qc)
+                end if
+              end if
+            end do
+
+            dq(:) = 2.5d-01 * (qr(:) - ql(:))
+
+            if (positive(p)) then
+              if (qn(p,i,j,k) > 0.0d+00) then
+                dq1 = sum(abs(dq(1:NDIMS)))
+                if (qn(p,i,j,k) <= dq1) &
+                  dq(:) = (5.0d-01 * qn(p,i,j,k) / dq1) * dq(:)
+              else
+                call print_message(loc, "Positive variable is not positive!")
+                status = 1
+                return
+              end if
+            end if
+
+            dq1 = dq(1) + dq(2) + dq(3)
+            dq2 = dq(1) - dq(2) - dq(3)
+            dq3 = dq(1) - dq(2) + dq(3)
+            dq4 = dq(1) + dq(2) - dq(3)
+
+            qb(p,is,js,ks) = qn(p,i,j,k) - dq1
+            qb(p,it,js,ks) = qn(p,i,j,k) + dq2
+            qb(p,is,jt,ks) = qn(p,i,j,k) - dq3
+            qb(p,it,jt,ks) = qn(p,i,j,k) + dq4
+            qb(p,is,js,kt) = qn(p,i,j,k) - dq4
+            qb(p,it,js,kt) = qn(p,i,j,k) + dq3
+            qb(p,is,jt,kt) = qn(p,i,j,k) - dq2
+            qb(p,it,jt,kt) = qn(p,i,j,k) + dq1
+
+          end do
+
+        end do
+      end do
+    end do
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine block_face_prolong_ppm
 #endif /* NDIMS == 3 */
 !
 !===============================================================================
@@ -4052,11 +4286,12 @@ module boundaries
 !
 !===============================================================================
 !
-! subroutine BLOCK_EDGE_PROLONG:
-! -----------------------------
+! subroutine BLOCK_EDGE_PROLONG_TVD:
+! ---------------------------------
 !
-!   Subroutine prolongss the region of qn specified by the edge position
-!   and inserts it in the corresponding edge boundary region of qb.
+!   This subroutine uses Total Variation Diminishing Method (TVD) to prolong
+!   a region of qn, specified by the edge position, and inserts it into
+!   the corresponding edge boundary region of qb.
 !
 !   Arguments:
 !
@@ -4068,7 +4303,7 @@ module boundaries
 !
 !===============================================================================
 !
-  subroutine block_edge_prolong(dir, pos, qn, qb, status)
+  subroutine block_edge_prolong_tvd(dir, pos, qn, qb, status)
 
     use coordinates   , only : nb, ne, ncells_half, nghosts_half
     use equations     , only : nv, positive
@@ -4101,7 +4336,7 @@ module boundaries
     logical                , save :: first = .true.
     integer, dimension(2,2), save :: nlims, tlims
 
-    character(len=*), parameter :: loc = 'BOUNDARIES::block_edge_prolong()'
+    character(len=*), parameter :: loc = 'BOUNDARIES::block_edge_prolong_tvd()'
 
 !-------------------------------------------------------------------------------
 !
@@ -4227,7 +4462,213 @@ module boundaries
 
 !-------------------------------------------------------------------------------
 !
-  end subroutine block_edge_prolong
+  end subroutine block_edge_prolong_tvd
+!
+!===============================================================================
+!
+! subroutine BLOCK_EDGE_PROLONG_PPM:
+! ---------------------------------
+!
+!   This subroutine uses Colella and Woodward's Piecewise Parabolic Method (PPM)
+!   to prolong a region of qn, specified by the edge position, and inserts it
+!   into the corresponding edge boundary region of qb.
+!
+!   Arguments:
+!
+!     dir    - the edge direction;
+!     pos    - the edge position;
+!     qn     - the array to prolong from;
+!     qb     - the array of the edge boundary region;
+!     status - the call status;
+!
+!   References:
+!
+!     [1] Colella, P., Woodward, P. R.,
+!         "The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations",
+!         Journal of Computational Physics, 1984, vol. 54, pp. 174-201,
+!         https://doi.org/10.1016/0021-9991(84)90143-8
+!
+!===============================================================================
+!
+  subroutine block_edge_prolong_ppm(dir, pos, qn, qb, status)
+
+    use coordinates   , only : nb, ne, ncells_half, nghosts_half
+    use equations     , only : nv, positive
+    use helpers       , only : print_message
+    use interpolations, only : limiter_prol
+
+    implicit none
+
+    integer                         , intent(in)    :: dir
+    integer     , dimension(3)      , intent(in)    :: pos
+    real(kind=8), dimension(:,:,:,:), intent(in)    :: qn
+    real(kind=8), dimension(:,:,:,:), intent(inout) :: qb
+    integer                         , intent(out)   :: status
+
+    integer      :: p, n
+    integer      :: i, il, iu, is, it, im, ip
+    integer      :: j, jl, ju, js, jt, jm, jp
+#if NDIMS == 3
+    integer      :: k, kl, ku, ks, kt, km, kp
+#else /* NDIMS == 3 */
+    integer      :: k
+#endif /* NDIMS == 3 */
+    real(kind=8) :: qc, dq1, dq2
+#if NDIMS == 3
+    real(kind=8) :: dq3, dq4
+#endif /* NDIMS == 3 */
+
+    real(kind=8), dimension(NDIMS) :: ql, qr, dq
+
+    logical                , save :: first = .true.
+    integer, dimension(2,2), save :: nlims, tlims
+
+    character(len=*), parameter :: loc = 'BOUNDARIES::block_edge_prolong_ppm()'
+
+!-------------------------------------------------------------------------------
+!
+    if (first) then
+      nlims(1,1) = ne - nghosts_half + 1
+      nlims(1,2) = ne
+      nlims(2,1) = nb
+      nlims(2,2) = nb + nghosts_half - 1
+
+      tlims(1,1) = nb
+      tlims(1,2) = nb + ncells_half + nghosts_half - 1
+      tlims(2,1) = ne - ncells_half - nghosts_half + 1
+      tlims(2,2) = ne
+
+      first = .false.
+    end if
+
+    status = 0
+
+    il = nlims(pos(1),1)
+    iu = nlims(pos(1),2)
+    jl = nlims(pos(2),1)
+    ju = nlims(pos(2),2)
+#if NDIMS == 3
+    kl = nlims(pos(3),1)
+    ku = nlims(pos(3),2)
+#endif /* NDIMS == 3 */
+
+    select case(dir)
+    case(1)
+      il = tlims(pos(1),1)
+      iu = tlims(pos(1),2)
+    case(2)
+      jl = tlims(pos(2),1)
+      ju = tlims(pos(2),2)
+#if NDIMS == 3
+    case(3)
+      kl = tlims(pos(3),1)
+      ku = tlims(pos(3),2)
+#endif /* NDIMS == 3 */
+    end select
+
+#if NDIMS == 3
+    do k = kl, ku
+      km = k - 1
+      kp = k + 1
+      ks = 2 * (k - kl) + 1
+      kt = ks + 1
+#else /* NDIMS == 3 */
+    k  = 1
+#endif /* NDIMS == 3 */
+      do j = jl, ju
+        jm = j - 1
+        jp = j + 1
+        js = 2 * (j - jl) + 1
+        jt = js + 1
+        do i = il, iu
+          im = i - 1
+          ip = i + 1
+          is = 2 * (i - il) + 1
+          it = is + 1
+
+          do p = 1, nv
+
+            qc    = qn(p,i,j,k)
+
+            ql(1) = (7.0d+00 * (qn(p,i-1,j,k) + qn(p,i  ,j,k)) &
+                             - (qn(p,i-2,j,k) + qn(p,i+1,j,k))) / 1.2d+01
+            qr(1) = (7.0d+00 * (qn(p,i  ,j,k) + qn(p,i+1,j,k)) &
+                             - (qn(p,i-1,j,k) + qn(p,i+2,j,k))) / 1.2d+01
+            ql(2) = (7.0d+00 * (qn(p,i,j-1,k) + qn(p,i,j  ,k)) &
+                             - (qn(p,i,j-2,k) + qn(p,i,j+1,k))) / 1.2d+01
+            qr(2) = (7.0d+00 * (qn(p,i,j  ,k) + qn(p,i,j+1,k)) &
+                             - (qn(p,i,j-1,k) + qn(p,i,j+2,k))) / 1.2d+01
+#if NDIMS == 3
+            ql(3) = (7.0d+00 * (qn(p,i,j,k-1) + qn(p,i,j,k  )) &
+                             - (qn(p,i,j,k-2) + qn(p,i,j,k+1))) / 1.2d+01
+            qr(3) = (7.0d+00 * (qn(p,i,j,k  ) + qn(p,i,j,k+1)) &
+                             - (qn(p,i,j,k-1) + qn(p,i,j,k+2))) / 1.2d+01
+#endif /* NDIMS == 3 */
+
+            do n = 1, NDIMS
+              if ((qr(n) - qc) * (ql(n) - qc) > 0.0d+00) then
+                ql(n) = qc
+                qr(n) = qc
+              else
+                if (abs(qr(n) - qc) >= 2.0d+00 * abs(ql(n) - qc)) then
+                  qr(n) = qc - 2.0d+00 * (ql(n) - qc)
+                end if
+                if (abs(ql(n) - qc) >= 2.0d+00 * abs(qr(n) - qc)) then
+                  ql(n) = qc - 2.0d+00 * (qr(n) - qc)
+                end if
+              end if
+            end do
+
+            dq(:) = 2.5d-01 * (qr(:) - ql(:))
+
+            if (positive(p)) then
+              if (qn(p,i,j,k) > 0.0d+00) then
+                dq1 = sum(abs(dq(1:NDIMS)))
+                if (qn(p,i,j,k) <= dq1) &
+                  dq(:) = (5.0d-01 * qn(p,i,j,k) / dq1) * dq(:)
+              else
+                call print_message(loc, "Positive variable is not positive!")
+                status = 1
+                return
+              end if
+            end if
+
+#if NDIMS == 2
+            dq1 = dq(1) + dq(2)
+            dq2 = dq(1) - dq(2)
+
+            qb(p,is,js,k ) = qn(p,i,j,k) - dq1
+            qb(p,it,js,k ) = qn(p,i,j,k) + dq2
+            qb(p,is,jt,k ) = qn(p,i,j,k) - dq2
+            qb(p,it,jt,k ) = qn(p,i,j,k) + dq1
+#endif /* NDIMS == 2 */
+#if NDIMS == 3
+            dq1 = dq(1) + dq(2) + dq(3)
+            dq2 = dq(1) - dq(2) - dq(3)
+            dq3 = dq(1) - dq(2) + dq(3)
+            dq4 = dq(1) + dq(2) - dq(3)
+
+            qb(p,is,js,ks) = qn(p,i,j,k) - dq1
+            qb(p,it,js,ks) = qn(p,i,j,k) + dq2
+            qb(p,is,jt,ks) = qn(p,i,j,k) - dq3
+            qb(p,it,jt,ks) = qn(p,i,j,k) + dq4
+            qb(p,is,js,kt) = qn(p,i,j,k) - dq4
+            qb(p,it,js,kt) = qn(p,i,j,k) + dq3
+            qb(p,is,jt,kt) = qn(p,i,j,k) - dq2
+            qb(p,it,jt,kt) = qn(p,i,j,k) + dq1
+#endif /* NDIMS == 3 */
+
+          end do
+
+        end do
+      end do
+#if NDIMS == 3
+    end do
+#endif /* NDIMS == 3 */
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine block_edge_prolong_ppm
 !
 !===============================================================================
 !
@@ -4322,11 +4763,12 @@ module boundaries
 !
 !===============================================================================
 !
-! subroutine BLOCK_CORNER_PROLONG:
-! -------------------------------
+! subroutine BLOCK_CORNER_PROLONG_TVD:
+! -----------------------------------
 !
-!   Subroutine prolongs the region of qn specified by the corner position
-!   and inserts it in the corresponding corner boundary region of qb.
+!   This subroutine uses Total Variation Diminishing Method (TVD) to prolong
+!   a region of qn, specified by the corner position, and inserts it into
+!   the corresponding corner boundary region of qb.
 !
 !   Remarks:
 !
@@ -4342,7 +4784,7 @@ module boundaries
 !
 !===============================================================================
 !
-  subroutine block_corner_prolong(pos, qn, qb, status)
+  subroutine block_corner_prolong_tvd(pos, qn, qb, status)
 
     use coordinates   , only : nb, ne, nghosts_half
     use equations     , only : nv, positive
@@ -4374,7 +4816,7 @@ module boundaries
     logical                , save :: first = .true.
     integer, dimension(2,2), save :: nlims
 
-    character(len=*), parameter :: loc = 'BOUNDARIES::block_corner_prolong()'
+    character(len=*), parameter :: loc = 'BOUNDARIES::block_corner_prolong_tvd()'
 
 !-------------------------------------------------------------------------------
 !
@@ -4479,7 +4921,195 @@ module boundaries
 
 !-------------------------------------------------------------------------------
 !
-  end subroutine block_corner_prolong
+  end subroutine block_corner_prolong_tvd
+!
+!===============================================================================
+!
+! subroutine BLOCK_CORNER_PROLONG_PPM:
+! -----------------------------------
+!
+!   This subroutine uses Colella and Woodward's Piecewise Parabolic Method (PPM)
+!   to prolong a region of qn, specified by the corner position, and inserts it
+!   into the corresponding corner boundary region of qb.
+!
+!   Remarks:
+!
+!     This subroutine requires the ghost zone regions of the qn already
+!     updated in order to perform a consistent interpolation.
+!
+!   Arguments:
+!
+!     pos    - the corner position;
+!     qn     - the array to prolong from;
+!     qb     - the array of the corner boundary region;
+!     status - the call status;
+!
+!   References:
+!
+!     [1] Colella, P., Woodward, P. R.,
+!         "The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations",
+!         Journal of Computational Physics, 1984, vol. 54, pp. 174-201,
+!         https://doi.org/10.1016/0021-9991(84)90143-8
+!
+!===============================================================================
+!
+  subroutine block_corner_prolong_ppm(pos, qn, qb, status)
+
+    use coordinates   , only : nb, ne, nghosts_half
+    use equations     , only : nv, positive
+    use helpers       , only : print_message
+    use interpolations, only : limiter_prol
+
+    implicit none
+
+    integer     , dimension(3)      , intent(in)    :: pos
+    real(kind=8), dimension(:,:,:,:), intent(in)    :: qn
+    real(kind=8), dimension(:,:,:,:), intent(inout) :: qb
+    integer                         , intent(out)   :: status
+
+    integer      :: p, n
+    integer      :: i, il, iu, im, ip, is, it
+    integer      :: j, jl, ju, jm, jp, js, jt
+#if NDIMS == 3
+    integer      :: k, kl, ku, km, kp, ks, kt
+#else /* NDIMS == 3 */
+    integer      :: k
+#endif /* NDIMS == 3 */
+    real(kind=8) :: qc, dq1, dq2
+#if NDIMS == 3
+    real(kind=8) :: dq3, dq4
+#endif /* NDIMS == 3 */
+
+    real(kind=8), dimension(NDIMS) :: ql, qr, dq
+
+    logical                , save :: first = .true.
+    integer, dimension(2,2), save :: nlims
+
+    character(len=*), parameter :: loc = 'BOUNDARIES::block_corner_prolong_ppm()'
+
+!-------------------------------------------------------------------------------
+!
+    if (first) then
+      nlims(1,1) = ne - nghosts_half + 1
+      nlims(1,2) = ne
+      nlims(2,1) = nb
+      nlims(2,2) = nb + nghosts_half - 1
+
+      first = .false.
+    end if
+
+    status = 0
+
+    il = nlims(pos(1),1)
+    iu = nlims(pos(1),2)
+    jl = nlims(pos(2),1)
+    ju = nlims(pos(2),2)
+#if NDIMS == 3
+    kl = nlims(pos(3),1)
+    ku = nlims(pos(3),2)
+
+    do k = kl, ku
+      km = k - 1
+      kp = k + 1
+      ks = 2 * (k - kl) + 1
+      kt = ks + 1
+#else /* NDIMS == 3 */
+    k  = 1
+#endif /* NDIMS == 3 */
+      do j = jl, ju
+        jm = j - 1
+        jp = j + 1
+        js = 2 * (j - jl) + 1
+        jt = js + 1
+        do i = il, iu
+          im = i - 1
+          ip = i + 1
+          is = 2 * (i - il) + 1
+          it = is + 1
+
+          do p = 1, nv
+
+            qc    = qn(p,i,j,k)
+
+            ql(1) = (7.0d+00 * (qn(p,i-1,j,k) + qn(p,i  ,j,k)) &
+                             - (qn(p,i-2,j,k) + qn(p,i+1,j,k))) / 1.2d+01
+            qr(1) = (7.0d+00 * (qn(p,i  ,j,k) + qn(p,i+1,j,k)) &
+                             - (qn(p,i-1,j,k) + qn(p,i+2,j,k))) / 1.2d+01
+            ql(2) = (7.0d+00 * (qn(p,i,j-1,k) + qn(p,i,j  ,k)) &
+                             - (qn(p,i,j-2,k) + qn(p,i,j+1,k))) / 1.2d+01
+            qr(2) = (7.0d+00 * (qn(p,i,j  ,k) + qn(p,i,j+1,k)) &
+                             - (qn(p,i,j-1,k) + qn(p,i,j+2,k))) / 1.2d+01
+#if NDIMS == 3
+            ql(3) = (7.0d+00 * (qn(p,i,j,k-1) + qn(p,i,j,k  )) &
+                             - (qn(p,i,j,k-2) + qn(p,i,j,k+1))) / 1.2d+01
+            qr(3) = (7.0d+00 * (qn(p,i,j,k  ) + qn(p,i,j,k+1)) &
+                             - (qn(p,i,j,k-1) + qn(p,i,j,k+2))) / 1.2d+01
+#endif /* NDIMS == 3 */
+
+            do n = 1, NDIMS
+              if ((qr(n) - qc) * (ql(n) - qc) > 0.0d+00) then
+                ql(n) = qc
+                qr(n) = qc
+              else
+                if (abs(qr(n) - qc) >= 2.0d+00 * abs(ql(n) - qc)) then
+                  qr(n) = qc - 2.0d+00 * (ql(n) - qc)
+                end if
+                if (abs(ql(n) - qc) >= 2.0d+00 * abs(qr(n) - qc)) then
+                  ql(n) = qc - 2.0d+00 * (qr(n) - qc)
+                end if
+              end if
+            end do
+
+            dq(:) = 2.5d-01 * (qr(:) - ql(:))
+
+            if (positive(p)) then
+              if (qn(p,i,j,k) > 0.0d+00) then
+                dq1 = sum(abs(dq(1:NDIMS)))
+                if (qn(p,i,j,k) <= dq1) &
+                  dq(:) = (5.0d-01 * qn(p,i,j,k) / dq1) * dq(:)
+              else
+                call print_message(loc, "Positive variable is not positive!")
+                status = 1
+                return
+              end if
+            end if
+
+#if NDIMS == 2
+            dq1 = dq(1) + dq(2)
+            dq2 = dq(1) - dq(2)
+
+            qb(p,is,js,k ) = qn(p,i,j,k) - dq1
+            qb(p,it,js,k ) = qn(p,i,j,k) + dq2
+            qb(p,is,jt,k ) = qn(p,i,j,k) - dq2
+            qb(p,it,jt,k ) = qn(p,i,j,k) + dq1
+#endif /* NDIMS == 2 */
+#if NDIMS == 3
+            dq1 = dq(1) + dq(2) + dq(3)
+            dq2 = dq(1) - dq(2) - dq(3)
+            dq3 = dq(1) - dq(2) + dq(3)
+            dq4 = dq(1) + dq(2) - dq(3)
+
+            qb(p,is,js,ks) = qn(p,i,j,k) - dq1
+            qb(p,it,js,ks) = qn(p,i,j,k) + dq2
+            qb(p,is,jt,ks) = qn(p,i,j,k) - dq3
+            qb(p,it,jt,ks) = qn(p,i,j,k) + dq4
+            qb(p,is,js,kt) = qn(p,i,j,k) - dq4
+            qb(p,it,js,kt) = qn(p,i,j,k) + dq3
+            qb(p,is,jt,kt) = qn(p,i,j,k) - dq2
+            qb(p,it,jt,kt) = qn(p,i,j,k) + dq1
+#endif /* NDIMS == 3 */
+
+          end do
+
+        end do
+      end do
+#if NDIMS == 3
+    end do
+#endif /* NDIMS == 3 */
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine block_corner_prolong_ppm
 !
 !===============================================================================
 !
diff --git a/sources/interpolations.F90 b/sources/interpolations.F90
index 6b937bf..6f7b3a9 100644
--- a/sources/interpolations.F90
+++ b/sources/interpolations.F90
@@ -79,6 +79,7 @@ module interpolations
 
 ! method names
 !
+  character(len=4)         :: prolongation_method = "TVD"
   character(len=255), save :: name_rec  = ""
   character(len=255), save :: name_tlim = ""
   character(len=255), save :: name_plim = ""
@@ -169,7 +170,7 @@ module interpolations
   public :: print_interpolations
   public :: interfaces, reconstruct, limiter_prol
   public :: fix_positivity
-  public :: order
+  public :: order, prolongation_method
 
 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 !
@@ -242,6 +243,13 @@ module interpolations
     call get_parameter("ppm_const"           , ppm_const      )
     call get_parameter("central_weight"      , cweight        )
     call get_parameter("cfl"                 , cfl            )
+    call get_parameter("prolongation_method" , prolongation_method)
+    select case(trim(prolongation_method))
+    case ("ppm", "PPM")
+      prolongation_method = 'PPM'
+    case default
+      prolongation_method = 'TVD'
+    end select
 
 ! calculate κ = (1 - ν) / ν
 !
@@ -811,7 +819,9 @@ module interpolations
         call print_parameter(verbose, "WENO sensitivity", weno_sensitivity)
       end if
       call print_parameter(verbose, "TVD limiter"         , name_tlim)
-      call print_parameter(verbose, "prolongation limiter", name_plim)
+      call print_parameter(verbose, "prolongation method" , prolongation_method)
+      if (prolongation_method == 'TVD') &
+        call print_parameter(verbose, "prolongation limiter", name_plim)
       if (mlp) then
         call print_parameter(verbose, "MLP limiting"        , "on"     )
       else
diff --git a/sources/mesh.F90 b/sources/mesh.F90
index 2e3b143..9351a8c 100644
--- a/sources/mesh.F90
+++ b/sources/mesh.F90
@@ -46,6 +46,15 @@ module mesh
     end subroutine
   end interface
 
+  abstract interface
+    subroutine prolong_block_iface(pmeta, status)
+      use blocks, only : block_meta
+      implicit none
+      type(block_meta), pointer, intent(inout) :: pmeta
+      integer                  , intent(out)   :: status
+    end subroutine
+  end interface
+
 ! pointer to the problem domain setup subroutine
 !
   procedure(setup_domain_iface) , pointer, save :: setup_domain => null()
@@ -54,6 +63,11 @@ module mesh
 !
   procedure(setup_problem_iface), pointer, save :: setup_problem => null()
 
+! pointer to the block prolongation subroutine
+!
+  procedure(prolong_block_iface), pointer, save :: &
+                                            prolong_block => prolong_block_tvd
+
 ! the flag indicating that the block structure or distribution has changed
 !
   logical, save :: changed = .true.
@@ -96,8 +110,9 @@ module mesh
 !
   subroutine initialize_mesh(verbose, status)
 
-    use coordinates, only : ncells, nghosts
-    use refinement , only : initialize_refinement
+    use coordinates   , only : ncells, nghosts
+    use interpolations, only : prolongation_method
+    use refinement    , only : initialize_refinement
 
     implicit none
 
@@ -108,6 +123,13 @@ module mesh
 !
     status = 0
 
+    select case(trim(prolongation_method))
+    case ("ppm", "PPM")
+      prolong_block => prolong_block_ppm
+    case default
+      prolong_block => prolong_block_tvd
+    end select
+
     setup_domain => setup_domain_default
 
     pm(1:NDIMS) = 2 * (ncells + nghosts)
@@ -1774,13 +1796,11 @@ module mesh
 !
 !===============================================================================
 !
-! subroutine PROLONG_BLOCK:
-! ------------------------
+! subroutine PROLONG_BLOCK_TVD:
+! ----------------------------
 !
-!   Subroutine prolongs variables from a data blocks linked to the input
-!   meta block and copy the resulting array of variables to
-!   the newly created children data block.  The process of data restriction
-!   conserves stored variables.
+!   This subroutine prolongs conserved variables from a parent data block
+!   to its children using the TVD method with a selected prolongation limiter.
 !
 !   Arguments:
 !
@@ -1789,7 +1809,7 @@ module mesh
 !
 !===============================================================================
 !
-  subroutine prolong_block(pmeta, status)
+  subroutine prolong_block_tvd(pmeta, status)
 
     use blocks        , only : block_meta, block_data, nchildren
     use coordinates   , only : ni => ncells, nn => bcells
@@ -1836,7 +1856,7 @@ module mesh
 !$  integer :: omp_get_thread_num
 !$omp threadprivate(first, nt, tmp)
 
-    character(len=*), parameter :: loc = 'MESH::prolong_block()'
+    character(len=*), parameter :: loc = 'MESH::prolong_block_tvd()'
 
 !-------------------------------------------------------------------------------
 !
@@ -1972,7 +1992,233 @@ module mesh
 
 !-------------------------------------------------------------------------------
 !
-  end subroutine prolong_block
+  end subroutine prolong_block_tvd
+!
+!===============================================================================
+!
+! subroutine PROLONG_BLOCK_PPM:
+! ----------------------------
+!
+!   This subroutine prolongs conserved variables from a parent data block to
+!   its children using Colella and Woodward's Piecewise Parabolic Method (PPM).
+!
+!   Arguments:
+!
+!     pmeta  - the input meta block;
+!     status - the subroutine call status: 0 for success, otherwise failure;
+!
+!   References:
+!
+!     [1] Colella, P., Woodward, P. R.,
+!         "The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations",
+!         Journal of Computational Physics, 1984, vol. 54, pp. 174-201,
+!         https://doi.org/10.1016/0021-9991(84)90143-8
+!
+!===============================================================================
+!
+  subroutine prolong_block_ppm(pmeta, status)
+
+    use blocks        , only : block_meta, block_data, nchildren
+    use coordinates   , only : ni => ncells, nn => bcells
+    use coordinates   , only : nh => nghosts_half
+    use coordinates   , only : nb, ne
+    use equations     , only : nv, positive
+    use helpers       , only : print_message
+    use interpolations, only : limiter_prol
+    use workspace     , only : resize_workspace, work, work_in_use
+
+    implicit none
+
+    type(block_meta), pointer, intent(inout) :: pmeta
+    integer                  , intent(out)   :: status
+
+    type(block_meta), pointer :: pchild
+    type(block_data), pointer :: pdata
+
+    logical, save :: first = .true.
+
+    integer      :: n, p, nl, nu
+    integer      :: i, im, ip
+    integer      :: j, jm, jp
+#if NDIMS == 3
+    integer      :: k, km, kp
+#else /* NDIMS == 3 */
+    integer      :: k
+#endif /* NDIMS == 3 */
+    real(kind=8) :: uc
+#if NDIMS == 3
+    real(kind=8) :: du1, du2, du3, du4
+#else /* NDIMS == 3 */
+    real(kind=8) :: du1, du2
+#endif /* NDIMS == 3 */
+
+    character(len=80) :: msg
+
+    integer     , dimension(NDIMS) :: l, u
+    real(kind=8), dimension(NDIMS) :: ul, ur, du
+
+    real(kind=8), dimension(:,:,:), pointer, save :: tmp
+
+    integer, save :: nt
+!$  integer :: omp_get_thread_num
+!$omp threadprivate(first, nt, tmp)
+
+    character(len=*), parameter :: loc = 'MESH::prolong_block_ppm()'
+
+!-------------------------------------------------------------------------------
+!
+    nt = 0
+!$  nt = omp_get_thread_num()
+    status = 0
+
+    if (first) then
+      n = product(pm(:))
+
+      call resize_workspace(n, status)
+      if (status /= 0) then
+        call print_message(loc, "Could not resize the workspace!")
+        go to 100
+      end if
+
+      tmp(1:pm(1),1:pm(2),1:pm(3)) => work(1:n,nt)
+
+      first = .false.
+    end if
+
+    if (work_in_use(nt)) &
+      call print_message(loc, "Workspace is being used right now! " //         &
+                              "Corruptions can occur!")
+
+    work_in_use(nt) = .true.
+
+#if NDIMS == 2
+    k  = 1
+#endif /* NDIMS == 2 */
+
+    nl = nb - nh
+    nu = ne + nh
+
+    pdata => pmeta%data
+
+    do n = 1, nv
+
+#if NDIMS == 3
+      do k = nl, nu
+        km   = k - 1
+        kp   = k + 1
+        l(3) = 2 * (k - nl) + 1
+        u(3) = l(3) + 1
+#endif /* NDIMS == 3 */
+        do j = nl, nu
+          jm   = j - 1
+          jp   = j + 1
+          l(2) = 2 * (j - nl) + 1
+          u(2) = l(2) + 1
+          do i = nl, nu
+            im   = i - 1
+            ip   = i + 1
+            l(1) = 2 * (i - nl) + 1
+            u(1) = l(1) + 1
+
+            uc    = pdata%u(n,i,j,k)
+
+            ul(1) = (7.0d+00 * (pdata%u(n,i-1,j,k) + pdata%u(n,i  ,j,k)) &
+                             - (pdata%u(n,i-2,j,k) + pdata%u(n,i+1,j,k))) / 1.2d+01
+            ur(1) = (7.0d+00 * (pdata%u(n,i  ,j,k) + pdata%u(n,i+1,j,k)) &
+                             - (pdata%u(n,i-1,j,k) + pdata%u(n,i+2,j,k))) / 1.2d+01
+            ul(2) = (7.0d+00 * (pdata%u(n,i,j-1,k) + pdata%u(n,i,j  ,k)) &
+                             - (pdata%u(n,i,j-2,k) + pdata%u(n,i,j+1,k))) / 1.2d+01
+            ur(2) = (7.0d+00 * (pdata%u(n,i,j  ,k) + pdata%u(n,i,j+1,k)) &
+                             - (pdata%u(n,i,j-1,k) + pdata%u(n,i,j+2,k))) / 1.2d+01
+#if NDIMS == 3
+            ul(3) = (7.0d+00 * (pdata%u(n,i,j,k-1) + pdata%u(n,i,j,k  )) &
+                             - (pdata%u(n,i,j,k-2) + pdata%u(n,i,j,k+1))) / 1.2d+01
+            ur(3) = (7.0d+00 * (pdata%u(n,i,j,k  ) + pdata%u(n,i,j,k+1)) &
+                             - (pdata%u(n,i,j,k-1) + pdata%u(n,i,j,k+2))) / 1.2d+01
+#endif /* NDIMS == 3 */
+
+            do p = 1, NDIMS
+              if ((ur(p) - uc) * (ul(p) - uc) > 0.0d+00) then
+                ul(p) = uc
+                ur(p) = uc
+              else
+                if (abs(ur(p) - uc) >= 2.0d+00 * abs(ul(p) - uc)) then
+                  ur(p) = uc - 2.0d+00 * (ul(p) - uc)
+                end if
+                if (abs(ul(p) - uc) >= 2.0d+00 * abs(ur(p) - uc)) then
+                  ul(p) = uc - 2.0d+00 * (ur(p) - uc)
+                end if
+              end if
+            end do
+
+            du(:) = 2.5d-01 * (ur(:) - ul(:))
+
+            if (positive(n) .and. pdata%u(n,i,j,k) < sum(abs(du(1:NDIMS)))) then
+              if (pdata%u(n,i,j,k) > 0.0d+00) then
+                do while (pdata%u(n,i,j,k) <= sum(abs(du(1:NDIMS))))
+                  du(:) = 0.5d+00 * du(:)
+                end do
+              else
+                write(msg,"(a,3i4,a)")                                         &
+                      "Positive variable is not positive at (", i, j, k, " )"
+                call print_message(loc, msg)
+                du(:) = 0.0d+00
+              end if
+            end if
+
+#if NDIMS == 2
+            du1 = du(1) + du(2)
+            du2 = du(1) - du(2)
+            tmp(l(1),l(2),k) = pdata%u(n,i,j,k) - du1
+            tmp(u(1),l(2),k) = pdata%u(n,i,j,k) + du2
+            tmp(l(1),u(2),k) = pdata%u(n,i,j,k) - du2
+            tmp(u(1),u(2),k) = pdata%u(n,i,j,k) + du1
+#endif /* NDIMS == 2 */
+
+#if NDIMS == 3
+            du1 = du(1) + du(2) + du(3)
+            du2 = du(1) - du(2) - du(3)
+            du3 = du(1) - du(2) + du(3)
+            du4 = du(1) + du(2) - du(3)
+            tmp(l(1),l(2),l(3)) = pdata%u(n,i,j,k) - du1
+            tmp(u(1),l(2),l(3)) = pdata%u(n,i,j,k) + du2
+            tmp(l(1),u(2),l(3)) = pdata%u(n,i,j,k) - du3
+            tmp(u(1),u(2),l(3)) = pdata%u(n,i,j,k) + du4
+            tmp(l(1),l(2),u(3)) = pdata%u(n,i,j,k) - du4
+            tmp(u(1),l(2),u(3)) = pdata%u(n,i,j,k) + du3
+            tmp(l(1),u(2),u(3)) = pdata%u(n,i,j,k) - du2
+            tmp(u(1),u(2),u(3)) = pdata%u(n,i,j,k) + du1
+#endif /* NDIMS == 3 */
+          end do
+        end do
+#if NDIMS == 3
+      end do
+#endif /* NDIMS == 3 */
+
+      do p = 1, nchildren
+
+        pchild => pmeta%child(p)%ptr
+
+        l(1:NDIMS) = pchild%pos(1:NDIMS) * ni + 1
+        u(1:NDIMS) =          l(1:NDIMS) + nn - 1
+
+#if NDIMS == 2
+        pchild%data%u(n,:,:,1) = tmp(l(1):u(1),l(2):u(2),k)
+#endif /* NDIMS == 2 */
+#if NDIMS == 3
+        pchild%data%u(n,:,:,:) = tmp(l(1):u(1),l(2):u(2),l(3):u(3))
+#endif /* NDIMS == 3 */
+
+      end do ! nchildren
+    end do ! n = 1, nv
+
+    work_in_use(nt) = .false.
+
+    100 continue
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine prolong_block_ppm
 !
 !===============================================================================
 !