diff --git a/sources/algebra.F90 b/sources/algebra.F90
index 7251e91..97bde8b 100644
--- a/sources/algebra.F90
+++ b/sources/algebra.F90
@@ -51,7 +51,7 @@ module algebra
   public :: quadratic, quadratic_normalized
   public :: cubic, cubic_normalized
   public :: quartic
-  public :: tridiag, invert
+  public :: tridiag, pentadiag, invert
 
 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 !
@@ -1148,6 +1148,87 @@ module algebra
 !
 !===============================================================================
 !
+! subroutine PENTADIAG:
+! --------------------
+!
+!   Subroutine solves the system of pentadiagonal linear equations using
+!   PTRANS-I algorithm.
+!
+!   Arguments:
+!
+!     n    - the size of the matrix;
+!     e(:) - the vector of the outer lower off-diagonal values;
+!     c(:) - the vector of the inner lower off-diagonal values;
+!     d(:) - the vector of the diagonal values;
+!     a(:) - the vector of the inner upper off-diagonal values;
+!     b(:) - the vector of the outer upper off-diagonal values;
+!     y(:) - the vector of the right side values;
+!     x(:) - the solution vector;
+!
+!   References:
+!
+!     [1] Askar, S. S., and Karawia, A. A.,
+!         "On Solving Pentadiagonal Linear Systems via Transformations",
+!         Mathematical Problems in Engineering,
+!         2015, Volume 2015, Article ID 232456, 9 pages
+!         http://dx.doi.org/10.1155/2015/232456
+!
+!===============================================================================
+!
+  subroutine pentadiag(n, e, c, d, a, b, y, x)
+
+    implicit none
+
+    integer                   , intent(in)  :: n
+    real(kind=8), dimension(n), intent(in)  :: e, c, d, a, b, y
+    real(kind=8), dimension(n), intent(out) :: x
+
+    integer                    :: i
+    real(kind=8), dimension(n) :: al, bt, gm, mu, z
+!
+!-------------------------------------------------------------------------------
+!
+    mu(1) = d(1)
+    al(1) = a(1) / mu(1)
+    bt(1) = b(1) / mu(1)
+    z (1) = y(1) / mu(1)
+
+    gm(2) = c(2)
+    mu(2) = d(2) - al(1) * gm(2)
+    al(2) = (a(2) - bt(1) * gm(2)) / mu(2)
+    bt(2) = b(2) / mu(2)
+    z (2) = (y(2) - z(1) * gm(2)) / mu(2)
+
+    do i = 3, n-2
+      gm(i) = c(i) - al(i-2) * e(i)
+      mu(i) = d(i) - bt(i-2) * e(i) - al(i-1) * gm(i)
+      al(i) = (a(i) - bt(i-1) * gm(i)) / mu(i)
+      bt(i) = b(i) / mu(i)
+      z (i) = (y(i) - z(i-2) * e(i) - z(i-1) * gm(i)) / mu(i)
+    end do
+
+    gm(n-1) = c(n-1) - al(n-3) * e(n-1)
+    mu(n-1) = d(n-1) - bt(n-3) * e(n-1) - al(n-2) * gm(n-1)
+    al(n-1) = (a(n-1) - bt(n-2) * gm(n-1)) / mu(n-1)
+
+    gm(n) = c(n) - al(n-2) * e(n)
+    mu(n) = d(n) - bt(n-3) * e(n) - al(n-1) * gm(n)
+
+    z(n-1) = (y(n-1) - z(n-2) * e(n-1) - z(n-2) * gm(n-1)) / mu(n-1)
+    z(n  ) = (y(n) - z(n-1) * e(n) - z(n-1) * gm(n)) / mu(n)
+
+    x(n  ) = z(n)
+    x(n-1) = z(n-1) - al(n-1) * x(n)
+    do i = n - 2, 1, -1
+      x(i) = z(i) - al(i) * x(i+1) - bt(i) * x(i+2)
+    end do
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine pentadiag
+!
+!===============================================================================
+!
 ! subroutine INVERT:
 ! -----------------
 !
diff --git a/sources/interpolations.F90 b/sources/interpolations.F90
index cf8b6e6..9a8d6e4 100644
--- a/sources/interpolations.F90
+++ b/sources/interpolations.F90
@@ -114,6 +114,10 @@ module interpolations
 !
   real(kind=8), save :: ppm_const  = 1.25d+00
 
+! weight toward central scheme for low-dissipation methods
+!
+  real(kind=8), save :: cweight = 0.0d+00
+
 ! Gaussian process reconstruction coefficients vector
 !
   real(kind=8), dimension(:)    , allocatable, save :: cgp
@@ -125,6 +129,29 @@ module interpolations
   logical     , save :: positivity = .false.
   logical     , save :: clip       = .false.
 
+! interpolation coefficients
+!
+  real(kind=8), dimension(10), save :: ce9a
+  real(kind=8), dimension( 8), save :: ce7a
+  real(kind=8), dimension( 6), save :: ce5a
+  real(kind=8), dimension( 3), save :: di5a, di7a
+  real(kind=8), dimension( 5), save :: di9a
+  real(kind=8), dimension( 4), save :: ci5a
+  real(kind=8), dimension( 6), save :: ci7a, ci9a
+
+  real(kind=8), dimension(9), parameter ::                                     &
+        ce9 = [ 4.000d+00,-4.100d+01, 1.990d+02,-6.410d+02, 1.879d+03,         &
+                1.375d+03,-3.050d+02, 5.500d+01,-5.000d+00 ] / 2.520d+03
+  real(kind=8), dimension(7), parameter ::                                     &
+        ce7 = [-3.000d+00, 2.500d+01,-1.010d+02, 3.190d+02,                    &
+                           2.140d+02,-3.800d+01, 4.000d+00 ] / 4.200d+02
+  real(kind=8), dimension(5), parameter ::                                     &
+        ce5 = [ 2.0d+00,-1.3d+01, 4.7d+01, 2.7d+01,-3.0d+00 ] / 6.0d+01
+  real(kind=8), dimension(3), parameter ::                                     &
+        ce3 = [-1.0d+00, 5.0d+00, 2.0d+00 ] / 6.0d+00
+  real(kind=8), dimension(2), parameter ::                                     &
+        ce2 = [ 5.0d-01, 5.0d-01 ]
+
 ! by default everything is private
 !
   private
@@ -225,12 +252,17 @@ module interpolations
     call get_parameter("kappa"               , kappa          )
     call get_parameter("kbeta"               , kbeta          )
     call get_parameter("ppm_const"           , ppm_const      )
+    call get_parameter("central_weight"      , cweight        )
     call get_parameter("cfl"                 , cfl            )
 
 ! calculate κ = (1 - ν) / ν
 !
     kappa = min(kappa, (1.0d+00 - cfl) / cfl)
 
+! correct central weight
+!
+    cweight = max(0.0d+00, min(1.0d+00, cweight))
+
 ! check ngp
 !
     if (mod(ngp,2) == 0 .or. ngp < 3) then
@@ -325,48 +357,24 @@ module interpolations
       reconstruct_states => reconstruct_mp9
       order              = 9
       nghosts            = max(nghosts, 6)
-    case ("mp5ld", "MP5LD")
-      name_rec           =  "5th order Low-Dissipation Monotonicity Preserving"
-      interfaces         => interfaces_dir
-      reconstruct_states => reconstruct_mp5ld
-      order              = 5
-      nghosts            = max(nghosts, 4)
-    case ("mp7ld", "MP7LD")
-      name_rec           =  "7th order Low-Dissipation Monotonicity Preserving"
-      interfaces         => interfaces_dir
-      reconstruct_states => reconstruct_mp7ld
-      order              = 7
-      nghosts            = max(nghosts, 4)
-    case ("mp9ld", "MP9LD")
-      name_rec           =  "9th order Low-Dissipation Monotonicity Preserving"
-      interfaces         => interfaces_dir
-      reconstruct_states => reconstruct_mp9ld
-      order              = 9
-      nghosts            = max(nghosts, 6)
     case ("crmp5", "CRMP5")
       name_rec           =  "5th order Compact Monotonicity Preserving"
       interfaces         => interfaces_dir
       reconstruct_states => reconstruct_crmp5
       order              = 5
       nghosts            = max(nghosts, 4)
-    case ("crmp5l", "crmp5ld", "CRMP5L", "CRMP5LD")
-      name_rec           =  "5th order Low-Dissipation Compact Monotonicity Preserving"
-      interfaces         => interfaces_dir
-      reconstruct_states => reconstruct_crmp5ld
-      order              = 5
-      nghosts            = max(nghosts, 4)
     case ("crmp7", "CRMP7")
       name_rec           =  "7th order Compact Monotonicity Preserving"
       interfaces         => interfaces_dir
       reconstruct_states => reconstruct_crmp7
       order              = 7
       nghosts            = max(nghosts, 4)
-    case ("crmp7l", "crmp7ld", "CRMP7L", "CRMP7LD")
-      name_rec           =  "7th order Low-Dissipation Compact Monotonicity Preserving"
+    case ("crmp9", "CRMP9")
+      name_rec           =  "9th order Compact Monotonicity Preserving"
       interfaces         => interfaces_dir
-      reconstruct_states => reconstruct_crmp7ld
-      order              = 7
-      nghosts            = max(nghosts, 4)
+      reconstruct_states => reconstruct_crmp9
+      order              = 9
+      nghosts            = max(nghosts, 6)
     case ("ocmp5", "OCMP5")
       name_rec           =  "5th order Optimized Compact Monotonicity Preserving"
       interfaces         => interfaces_dir
@@ -527,6 +535,39 @@ module interpolations
 !
     ng = nghosts
 
+! calculate coefficients of the low-dissipation methods
+!
+    ce5a = [ -           cweight +2.000d+00,  5.00d+00 * cweight -1.300d+01,   &
+             -1.00d+01 * cweight +4.700d+01,  1.00d+01 * cweight +2.700d+01,   &
+             -5.00d+00 * cweight -3.000d+00,             cweight ] / 6.00d+01
+    ce7a = [  3.00d+00 * cweight -6.000d+00, -2.10d+01 * cweight +5.000d+01,   &
+              6.30d+01 * cweight -2.020d+02, -1.05d+02 * cweight +6.380d+02,   &
+              1.05d+02 * cweight +4.280d+02, -6.30d+01 * cweight -7.600d+01,   &
+              2.10d+01 * cweight +8.000d+00, -3.00d+00 * cweight ] / 8.40d+02
+    ce9a = [ -2.00d+00 * cweight +4.000d+00,  1.80d+01 * cweight -4.100d+01,   &
+             -7.20d+01 * cweight +1.990d+02,  1.68d+02 * cweight -6.410d+02,   &
+             -2.52d+02 * cweight +1.879d+03,  2.52d+02 * cweight +1.375d+03,   &
+             -1.68d+02 * cweight -3.050d+02,  7.20d+01 * cweight +5.500d+01,   &
+             -1.80d+01 * cweight -5.000d+00,  2.00d+00 * cweight ] / 2.52d+03
+
+    di5a = [ -(cweight -3.0d+00) / 6.0d+00, 1.0d+00,                           &
+              (cweight +1.0d+00) / 6.0d+00 ]
+    ci5a = [ - cweight +2.0d+00, -9.0d+00 * cweight +3.8d+01,                  &
+               9.0d+00 * cweight +2.0d+01,  cweight ] / 3.6d+01
+
+    di7a = [ -(cweight -4.0d+00) / 8.0d+00, 1.0d+00,                           &
+              (cweight +2.0d+00) / 8.0d+00 ]
+    ci7a = [   cweight -2.0d+00, -1.50d+01 * cweight +3.8d+01,                 &
+              -8.0d+01 * cweight +4.78d+02, 8.0d+01 * cweight +3.18d+02,       &
+               1.5d+01 * cweight +8.00d+00, -cweight ] / 4.8d+02
+
+    di9a = [ -2.0d+00 * cweight +5.0d+00, -1.0d+01 * cweight +4.0d+01,         &
+              6.0d+01, 1.0d+01 * cweight +2.0d+01,                             &
+              2.0d+00 * cweight +1.0d+00 ] / 6.0d+01
+    ci9a = [ -3.00d+00 * cweight +6.000d+00,-1.75d+02 * cweight +4.810d+02,    &
+             -3.00d+02 * cweight +1.881d+03, 3.00d+02 * cweight +1.281d+03,    &
+              1.75d+02 * cweight +1.310d+02, 3.00d+00 * cweight ] / 1.8d+03
+
 #ifdef PROFILE
 ! stop accounting time for module initialization/finalization
 !
@@ -625,6 +666,7 @@ module interpolations
 !-------------------------------------------------------------------------------
 !
     if (verbose) then
+
       call print_section(verbose, "Interpolations")
       call print_parameter(verbose, "reconstruction"      , name_rec )
       call print_parameter(verbose, "TVD limiter"         , name_tlim)
@@ -645,6 +687,7 @@ module interpolations
       else
         call print_parameter(verbose, "clip extrema"        , "off"    )
       end if
+      call print_parameter(verbose, "central weight", cweight)
 
     end if
 
@@ -3794,8 +3837,8 @@ module interpolations
 ! subroutine RECONSTRUCT_MP5:
 ! --------------------------
 !
-!   Subroutine reconstructs the interface states using the fifth order
-!   Monotonicity Preserving (MP) method.
+!   Subroutine reconstructs the interface states using the 5th order
+!   Monotonicity Preserving (MP) method with adjustable dissipation.
 !
 !   Arguments are described in subroutine reconstruct().
 !
@@ -3807,105 +3850,66 @@ module interpolations
 !         Journal on Computational Physics,
 !         1997, vol. 136, pp. 83-99,
 !         http://dx.doi.org/10.1006/jcph.1997.5745
-!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
-!         "A 5th order monotonicity-preserving upwind compact difference
-!          scheme",
-!         Science China Physics, Mechanics and Astronomy,
-!         Volume 54, Issue 3, pp. 511-522,
-!         http://dx.doi.org/10.1007/s11433-010-4220-x
+!     [2] Jian Fang, Yufeng Yao, Zhaorui Li, Lipeng Lu,
+!         "Investigation of low-dissipation monotonicity-preserving scheme
+!          for direct numerical simulation of compressible turbulent flows",
+!         Computers & Fluids,
+!         2014, vol. 104, pp. 55-72,
+!         http://dx.doi.org/10.1016/j.compfluid.2014.07.024
 !
 !===============================================================================
 !
   subroutine reconstruct_mp5(h, fc, fl, fr)
 
-! local variables are not implicit by default
-!
     implicit none
 
-! subroutine arguments
-!
     real(kind=8)              , intent(in)  :: h
     real(kind=8), dimension(:), intent(in)  :: fc
     real(kind=8), dimension(:), intent(out) :: fl, fr
 
-! local variables
-!
     integer :: n, i
 
-! local arrays for derivatives
-!
     real(kind=8), dimension(size(fc)) :: fi
     real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
-    real(kind=8), dimension(5), parameter ::                                   &
-                                ce5 = (/ 2.0d+00,-1.3d+01, 4.7d+01             &
-                                                , 2.7d+01,-3.0d+00 /) / 6.0d+01
-    real(kind=8), dimension(3), parameter :: &
-                                ce3 = (/-1.0d+00, 5.0d+00, 2.0d+00 /) / 6.0d+00
-    real(kind=8), dimension(2), parameter :: ce2 = (/ 5.0d-01, 5.0d-01 /)
 !
 !-------------------------------------------------------------------------------
-!
-! get the input vector length
 !
     n = size(fc)
 
 !! === left-side interpolation ===
 !!
-! reconstruct the interface state using the 5th order interpolation
-!
-    do i = 3, n - 2
-      u(i) = sum(ce5(:) * fc(i-2:i+2))
+    do i = 3, n - 3
+      u(i) = sum(ce5a(:) * fc(i-2:i+3))
     end do
 
-! interpolate the interface state of the ghost zones using the interpolations
-! of lower orders
-!
-    u(  1) = sum(ce2(:) * fc(  1:  2))
-    u(  2) = sum(ce3(:) * fc(  1:  3))
-    u(n-1) = sum(ce3(:) * fc(n-2:  n))
-    u(n  ) =              fc(n      )
+    u(  1) = sum(ce2(:) * fc(  1:2))
+    u(  2) = sum(ce3(:) * fc(  1:3))
+    u(n-2) = sum(ce5(:) * fc(n-4:n))
+    u(n-1) = sum(ce3(:) * fc(n-2:n))
+    u(n  ) =              fc(n    )
 
-! apply the monotonicity preserving limiting
-!
     call mp_limiting(fc(:), u(:))
 
-! copy the interpolation to the respective vector
-!
     fl(1:n) = u(1:n)
 
 !! === right-side interpolation ===
 !!
-! invert the cell-centered value vector
-!
     fi(1:n) = fc(n:1:-1)
 
-! reconstruct the interface state using the 5th order interpolation
-!
-    do i = 3, n - 2
-      u(i) = sum(ce5(:) * fi(i-2:i+2))
+    do i = 3, n - 3
+      u(i) = sum(ce5a(:) * fi(i-2:i+3))
     end do
 
-! interpolate the interface state of the ghost zones using the interpolations
-! of lower orders
-!
-    u(  1) = sum(ce2(:) * fi(  1:  2))
-    u(  2) = sum(ce3(:) * fi(  1:  3))
-    u(n-1) = sum(ce3(:) * fi(n-2:  n))
-    u(n  ) =              fi(n      )
+    u(  1) = sum(ce2(:) * fi(  1:2))
+    u(  2) = sum(ce3(:) * fi(  1:3))
+    u(n-2) = sum(ce5(:) * fi(n-4:n))
+    u(n-1) = sum(ce3(:) * fi(n-2:n))
+    u(n  ) =              fi(n    )
 
-! apply the monotonicity preserving limiting
-!
     call mp_limiting(fi(:), u(:))
 
-! copy the interpolation to the respective vector
-!
     fr(1:n-1) = u(n-1:1:-1)
 
-! update the interpolation of the first and last points
-!
     i     = n - 1
     fl(1) = 0.5d+00 * (fc(1) + fc(2))
     fr(i) = 0.5d+00 * (fc(i) + fc(n))
@@ -3921,125 +3925,72 @@ module interpolations
 ! subroutine RECONSTRUCT_MP7:
 ! --------------------------
 !
-!   Subroutine reconstructs the interface states using the seventh order
-!   Monotonicity Preserving (MP) method.
+!   Subroutine reconstructs the interface states using the 7th order
+!   Monotonicity Preserving (MP) method with adjustable dissipation.
 !
 !   Arguments are described in subroutine reconstruct().
 !
 !   References:
 !
-!     [1] Suresh, A. & Huynh, H. T.,
-!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
-!          Time Stepping"
-!         Journal on Computational Physics,
-!         1997, vol. 136, pp. 83-99,
-!         http://dx.doi.org/10.1006/jcph.1997.5745
-!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
-!         "A 5th order monotonicity-preserving upwind compact difference
-!          scheme",
-!         Science China Physics, Mechanics and Astronomy,
-!         Volume 54, Issue 3, pp. 511-522,
-!         http://dx.doi.org/10.1007/s11433-010-4220-x
+!     - See RECONSTRUCT_MP5
 !
 !===============================================================================
 !
   subroutine reconstruct_mp7(h, fc, fl, fr)
 
-! local variables are not implicit by default
-!
     implicit none
 
-! subroutine arguments
-!
     real(kind=8)              , intent(in)  :: h
     real(kind=8), dimension(:), intent(in)  :: fc
     real(kind=8), dimension(:), intent(out) :: fl, fr
 
-! local variables
-!
     integer :: n, i
 
-! local arrays for derivatives
-!
     real(kind=8), dimension(size(fc)) :: fi
     real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
-    real(kind=8), dimension(7), parameter ::                                   &
-                                ce7 = (/-3.0d+00, 2.5d+01,-1.01d+02, 3.19d+02  &
-                                      , 2.14d+02,-3.8d+01, 4.0d+00 /) / 4.2d+02
-    real(kind=8), dimension(5), parameter ::                                   &
-                                ce5 = (/ 2.0d+00,-1.3d+01, 4.7d+01             &
-                                                , 2.7d+01,-3.0d+00 /) / 6.0d+01
-    real(kind=8), dimension(3), parameter :: &
-                                ce3 = (/-1.0d+00, 5.0d+00, 2.0d+00 /) / 6.0d+00
-    real(kind=8), dimension(2), parameter :: ce2 = (/ 5.0d-01, 5.0d-01 /)
 !
 !-------------------------------------------------------------------------------
-!
-! get the input vector length
 !
     n = size(fc)
 
 !! === left-side interpolation ===
 !!
-! reconstruct the interface state using the 5th order interpolation
-!
-    do i = 4, n - 3
-      u(i) = sum(ce7(:) * fc(i-3:i+3))
+    do i = 4, n - 4
+      u(i) = sum(ce7a(:) * fc(i-3:i+4))
     end do
 
-! interpolate the interface state of the ghost zones using the interpolations
-! of lower orders
-!
-    u(  1) = sum(ce2(:) * fc(  1:  2))
-    u(  2) = sum(ce3(:) * fc(  1:  3))
-    u(  3) = sum(ce5(:) * fc(  1:  5))
-    u(n-2) = sum(ce5(:) * fc(n-4:  n))
-    u(n-1) = sum(ce3(:) * fc(n-2:  n))
-    u(n  ) =              fc(n      )
+    u(  1) = sum(ce2(:) * fc(  1:2))
+    u(  2) = sum(ce3(:) * fc(  1:3))
+    u(  3) = sum(ce5(:) * fc(  1:5))
+    u(n-3) = sum(ce7(:) * fc(n-6:n))
+    u(n-2) = sum(ce5(:) * fc(n-4:n))
+    u(n-1) = sum(ce3(:) * fc(n-2:n))
+    u(n  ) =              fc(n    )
 
-! apply the monotonicity preserving limiting
-!
     call mp_limiting(fc(:), u(:))
 
-! copy the interpolation to the respective vector
-!
     fl(1:n) = u(1:n)
 
 !! === right-side interpolation ===
 !!
-! invert the cell-centered value vector
-!
     fi(1:n) = fc(n:1:-1)
 
-! reconstruct the interface state using the 5th order interpolation
-!
-    do i = 4, n - 3
-      u(i) = sum(ce7(:) * fi(i-3:i+3))
+    do i = 4, n - 4
+      u(i) = sum(ce7a(:) * fi(i-3:i+4))
     end do
 
-! interpolate the interface state of the ghost zones using the interpolations
-! of lower orders
-!
-    u(  1) = sum(ce2(:) * fi(  1:  2))
-    u(  2) = sum(ce3(:) * fi(  1:  3))
-    u(  3) = sum(ce5(:) * fi(  1:  5))
-    u(n-2) = sum(ce5(:) * fi(n-4:  n))
-    u(n-1) = sum(ce3(:) * fi(n-2:  n))
-    u(n  ) =              fi(n      )
+    u(  1) = sum(ce2(:) * fi(  1:2))
+    u(  2) = sum(ce3(:) * fi(  1:3))
+    u(  3) = sum(ce5(:) * fi(  1:5))
+    u(n-3) = sum(ce7(:) * fi(n-6:n))
+    u(n-2) = sum(ce5(:) * fi(n-4:n))
+    u(n-1) = sum(ce3(:) * fi(n-2:n))
+    u(n  ) =              fi(n    )
 
-! apply the monotonicity preserving limiting
-!
     call mp_limiting(fi(:), u(:))
 
-! copy the interpolation to the respective vector
-!
     fr(1:n-1) = u(n-1:1:-1)
 
-! update the interpolation of the first and last points
-!
     i     = n - 1
     fl(1) = 0.5d+00 * (fc(1) + fc(2))
     fr(i) = 0.5d+00 * (fc(i) + fc(n))
@@ -4055,131 +4006,76 @@ module interpolations
 ! subroutine RECONSTRUCT_MP9:
 ! --------------------------
 !
-!   Subroutine reconstructs the interface states using the ninth order
-!   Monotonicity Preserving (MP) method.
+!   Subroutine reconstructs the interface states using the 9th order
+!   Monotonicity Preserving (MP) method with adjustable dissipation.
 !
 !   Arguments are described in subroutine reconstruct().
 !
 !   References:
 !
-!     [1] Suresh, A. & Huynh, H. T.,
-!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
-!          Time Stepping"
-!         Journal on Computational Physics,
-!         1997, vol. 136, pp. 83-99,
-!         http://dx.doi.org/10.1006/jcph.1997.5745
-!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
-!         "A 5th order monotonicity-preserving upwind compact difference
-!          scheme",
-!         Science China Physics, Mechanics and Astronomy,
-!         Volume 54, Issue 3, pp. 511-522,
-!         http://dx.doi.org/10.1007/s11433-010-4220-x
+!     - See RECONSTRUCT_MP5LD
 !
 !===============================================================================
 !
   subroutine reconstruct_mp9(h, fc, fl, fr)
 
-! local variables are not implicit by default
-!
     implicit none
 
-! subroutine arguments
-!
     real(kind=8)              , intent(in)  :: h
     real(kind=8), dimension(:), intent(in)  :: fc
     real(kind=8), dimension(:), intent(out) :: fl, fr
 
-! local variables
-!
     integer :: n, i
 
-! local arrays for derivatives
-!
     real(kind=8), dimension(size(fc)) :: fi
     real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
-    real(kind=8), dimension(9), parameter ::                                   &
-             ce9 = (/ 4.000d+00,-4.100d+01, 1.990d+02,-6.410d+02, 1.879d+03,   &
-                      1.375d+03,-3.050d+02, 5.500d+01,-5.000d+00 /) / 2.520d+03
-    real(kind=8), dimension(7), parameter ::                                   &
-             ce7 = (/-3.000d+00, 2.500d+01,-1.010d+02, 3.190d+02, 2.140d+02,   &
-                     -3.800d+01, 4.000d+00 /) / 4.200d+02
-    real(kind=8), dimension(5), parameter ::                                   &
-             ce5 = (/ 2.0d+00,-1.3d+01, 4.7d+01, 2.7d+01,-3.0d+00 /) / 6.0d+01
-    real(kind=8), dimension(3), parameter ::                                   &
-             ce3 = (/-1.0d+00, 5.0d+00, 2.0d+00 /) / 6.0d+00
-    real(kind=8), dimension(2), parameter :: ce2 = (/ 5.0d-01, 5.0d-01 /)
 !
 !-------------------------------------------------------------------------------
-!
-! get the input vector length
 !
     n = size(fc)
 
 !! === left-side interpolation ===
 !!
-! reconstruct the interface state using the 9th order interpolation
-!
-    do i = 5, n - 4
-      u(i) = sum(ce9(:) * fc(i-4:i+4))
+    do i = 5, n - 5
+      u(i) = sum(ce9a(:) * fc(i-4:i+5))
     end do
 
-! interpolate the interface state of the ghost zones using the interpolations
-! of lower orders
-!
-    u(  1) = sum(ce2(:) * fc(  1:  2))
-    u(  2) = sum(ce3(:) * fc(  1:  3))
-    u(  3) = sum(ce5(:) * fc(  1:  5))
-    u(  4) = sum(ce7(:) * fc(  1:  7))
-    u(n-3) = sum(ce7(:) * fc(n-6:  n))
-    u(n-2) = sum(ce5(:) * fc(n-4:  n))
-    u(n-1) = sum(ce3(:) * fc(n-2:  n))
-    u(n  ) =              fc(n      )
+    u(  1) = sum(ce2(:) * fc(  1:2))
+    u(  2) = sum(ce3(:) * fc(  1:3))
+    u(  3) = sum(ce5(:) * fc(  1:5))
+    u(  4) = sum(ce7(:) * fc(  1:7))
+    u(n-4) = sum(ce9(:) * fc(n-8:n))
+    u(n-3) = sum(ce7(:) * fc(n-6:n))
+    u(n-2) = sum(ce5(:) * fc(n-4:n))
+    u(n-1) = sum(ce3(:) * fc(n-2:n))
+    u(n  ) =              fc(n    )
 
-! apply the monotonicity preserving limiting
-!
     call mp_limiting(fc(:), u(:))
 
-! copy the interpolation to the respective vector
-!
     fl(1:n) = u(1:n)
 
 !! === right-side interpolation ===
 !!
-! invert the cell-centered value vector
-!
     fi(1:n) = fc(n:1:-1)
 
-! reconstruct the interface state using the 9th order interpolation
-!
-    do i = 5, n - 4
-      u(i) = sum(ce9(:) * fi(i-4:i+4))
+    do i = 5, n - 5
+      u(i) = sum(ce9a(:) * fi(i-4:i+5))
     end do
 
-! interpolate the interface state of the ghost zones using the interpolations
-! of lower orders
-!
-    u(  1) = sum(ce2(:) * fi(  1:  2))
-    u(  2) = sum(ce3(:) * fi(  1:  3))
-    u(  3) = sum(ce5(:) * fi(  1:  5))
-    u(  4) = sum(ce7(:) * fi(  1:  7))
-    u(n-3) = sum(ce7(:) * fi(n-6:  n))
-    u(n-2) = sum(ce5(:) * fi(n-4:  n))
-    u(n-1) = sum(ce3(:) * fi(n-2:  n))
-    u(n  ) =              fi(n      )
+    u(  1) = sum(ce2(:) * fi(  1:2))
+    u(  2) = sum(ce3(:) * fi(  1:3))
+    u(  3) = sum(ce5(:) * fi(  1:5))
+    u(  4) = sum(ce7(:) * fi(  1:7))
+    u(n-4) = sum(ce9(:) * fi(n-8:n))
+    u(n-3) = sum(ce7(:) * fi(n-6:n))
+    u(n-2) = sum(ce5(:) * fi(n-4:n))
+    u(n-1) = sum(ce3(:) * fi(n-2:n))
+    u(n  ) =              fi(n    )
 
-! apply the monotonicity preserving limiting
-!
     call mp_limiting(fi(:), u(:))
 
-! copy the interpolation to the respective vector
-!
     fr(1:n-1) = u(n-1:1:-1)
 
-! update the interpolation of the first and last points
-!
     i     = n - 1
     fl(1) = 0.5d+00 * (fc(1) + fc(2))
     fr(i) = 0.5d+00 * (fc(i) + fc(n))
@@ -4192,422 +4088,12 @@ module interpolations
 !
 !===============================================================================
 !
-! subroutine RECONSTRUCT_MP5LD:
-! ----------------------------
-!
-!   Subroutine reconstructs the interface states using the fifth order
-!   low dissipation Monotonicity Preserving (MP) method.
-!
-!   Arguments are described in subroutine reconstruct().
-!
-!   References:
-!
-!     [1] Suresh, A. & Huynh, H. T.,
-!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
-!          Time Stepping"
-!         Journal on Computational Physics,
-!         1997, vol. 136, pp. 83-99,
-!         http://dx.doi.org/10.1006/jcph.1997.5745
-!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
-!         "A 5th order monotonicity-preserving upwind compact difference
-!          scheme",
-!         Science China Physics, Mechanics and Astronomy,
-!         Volume 54, Issue 3, pp. 511-522,
-!         http://dx.doi.org/10.1007/s11433-010-4220-x
-!
-!===============================================================================
-!
-  subroutine reconstruct_mp5ld(h, fc, fl, fr)
-
-! local variables are not implicit by default
-!
-    implicit none
-
-! subroutine arguments
-!
-    real(kind=8)              , intent(in)  :: h
-    real(kind=8), dimension(:), intent(in)  :: fc
-    real(kind=8), dimension(:), intent(out) :: fl, fr
-
-! local variables
-!
-    integer :: n, i
-
-! local arrays for derivatives
-!
-    real(kind=8), dimension(size(fc)) :: fi
-    real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
-    real(kind=8), dimension(6), parameter ::                                   &
-                              ce5 = [ 1.20d+01,-8.10d+01, 3.09d+02,            &
-                                      2.09d+02,-3.10d+01, 2.00d+00 ] / 4.2d+02
-    real(kind=8), dimension(3), parameter ::                                   &
-                              ce3 = [-1.00d+00, 5.00d+00, 2.00d+00 ] / 6.0d+00
-    real(kind=8), dimension(2), parameter ::                                   &
-                              ce2 = [ 5.0d-01, 5.0d-01 ]
-!
-!-------------------------------------------------------------------------------
-!
-! get the input vector length
-!
-    n = size(fc)
-
-!! === left-side interpolation ===
-!!
-! reconstruct the interface state using the 5th order interpolation
-!
-    do i = 3, n - 3
-      u(i) = sum(ce5(:) * fc(i-2:i+3))
-    end do
-
-! interpolate the interface state of the ghost zones using the interpolations
-! of lower orders
-!
-    u(  1) = sum(ce2(:) * fc(  1:  2))
-    u(  2) = sum(ce3(:) * fc(  1:  3))
-    u(n-2) = sum(ce3(:) * fc(n-3:n-1))
-    u(n-1) = sum(ce3(:) * fc(n-2:n  ))
-    u(n  ) =              fc(n      )
-
-! apply the monotonicity preserving limiting
-!
-    call mp_limiting(fc(:), u(:))
-
-! copy the interpolation to the respective vector
-!
-    fl(1:n) = u(1:n)
-
-!! === right-side interpolation ===
-!!
-! invert the cell-centered value vector
-!
-    fi(1:n) = fc(n:1:-1)
-
-! reconstruct the interface state using the 5th order interpolation
-!
-    do i = 3, n - 3
-      u(i) = sum(ce5(:) * fi(i-2:i+3))
-    end do
-
-! interpolate the interface state of the ghost zones using the interpolations
-! of lower orders
-!
-    u(  1) = sum(ce2(:) * fi(  1:  2))
-    u(  2) = sum(ce3(:) * fi(  1:  3))
-    u(n-2) = sum(ce3(:) * fi(n-3:n-1))
-    u(n-1) = sum(ce3(:) * fi(n-2:n  ))
-    u(n  ) =              fi(n      )
-
-! apply the monotonicity preserving limiting
-!
-    call mp_limiting(fi(:), u(:))
-
-! copy the interpolation to the respective vector
-!
-    fr(1:n-1) = u(n-1:1:-1)
-
-! update the interpolation of the first and last points
-!
-    i     = n - 1
-    fl(1) = 0.5d+00 * (fc(1) + fc(2))
-    fr(i) = 0.5d+00 * (fc(i) + fc(n))
-    fl(n) = fc(n)
-    fr(n) = fc(n)
-
-!-------------------------------------------------------------------------------
-!
-  end subroutine reconstruct_mp5ld
-!
-!===============================================================================
-!
-! subroutine RECONSTRUCT_MP7LD:
-! ----------------------------
-!
-!   Subroutine reconstructs the interface states using the seventh order
-!   low dissipation Monotonicity Preserving (MP) method.
-!
-!   Arguments are described in subroutine reconstruct().
-!
-!   References:
-!
-!     [1] Suresh, A. & Huynh, H. T.,
-!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
-!          Time Stepping"
-!         Journal on Computational Physics,
-!         1997, vol. 136, pp. 83-99,
-!         http://dx.doi.org/10.1006/jcph.1997.5745
-!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
-!         "A 5th order monotonicity-preserving upwind compact difference
-!          scheme",
-!         Science China Physics, Mechanics and Astronomy,
-!         Volume 54, Issue 3, pp. 511-522,
-!         http://dx.doi.org/10.1007/s11433-010-4220-x
-!
-!===============================================================================
-!
-  subroutine reconstruct_mp7ld(h, fc, fl, fr)
-
-! local variables are not implicit by default
-!
-    implicit none
-
-! subroutine arguments
-!
-    real(kind=8)              , intent(in)  :: h
-    real(kind=8), dimension(:), intent(in)  :: fc
-    real(kind=8), dimension(:), intent(out) :: fl, fr
-
-! local variables
-!
-    integer :: n, i
-
-! local arrays for derivatives
-!
-    real(kind=8), dimension(size(fc)) :: fi
-    real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
-    real(kind=8), dimension(8), parameter ::                                   &
-                    ce7 = [-8.00d+00, 6.80d+01,-2.82d+02, 9.22d+02,            &
-                            6.77d+02,-1.35d+02, 1.9d+01, -1.0d+00 ] / 1.26d+03
-    real(kind=8), dimension(5), parameter ::                                   &
-                    ce5 = [ 2.0d+00,-1.3d+01, 4.7d+01                          &
-                                            , 2.7d+01,-3.0d+00 ] / 6.0d+01
-    real(kind=8), dimension(3), parameter ::                                   &
-                    ce3 = [-1.0d+00, 5.0d+00, 2.0d+00 ] / 6.0d+00
-    real(kind=8), dimension(2), parameter ::                                   &
-                    ce2 = [ 5.0d-01, 5.0d-01 ]
-!
-!-------------------------------------------------------------------------------
-!
-! get the input vector length
-!
-    n = size(fc)
-
-!! === left-side interpolation ===
-!!
-! reconstruct the interface state using the 5th order interpolation
-!
-    do i = 4, n - 4
-      u(i) = sum(ce7(:) * fc(i-3:i+4))
-    end do
-
-! interpolate the interface state of the ghost zones using the interpolations
-! of lower orders
-!
-    u(  1) = sum(ce2(:) * fc(  1:  2))
-    u(  2) = sum(ce3(:) * fc(  1:  3))
-    u(  3) = sum(ce5(:) * fc(  1:  5))
-    u(n-3) = sum(ce5(:) * fc(n-5:n-1))
-    u(n-2) = sum(ce5(:) * fc(n-4:n  ))
-    u(n-1) = sum(ce3(:) * fc(n-2:n  ))
-    u(n  ) =              fc(n      )
-
-! apply the monotonicity preserving limiting
-!
-    call mp_limiting(fc(:), u(:))
-
-! copy the interpolation to the respective vector
-!
-    fl(1:n) = u(1:n)
-
-!! === right-side interpolation ===
-!!
-! invert the cell-centered value vector
-!
-    fi(1:n) = fc(n:1:-1)
-
-! reconstruct the interface state using the 5th order interpolation
-!
-    do i = 4, n - 4
-      u(i) = sum(ce7(:) * fi(i-3:i+4))
-    end do
-
-! interpolate the interface state of the ghost zones using the interpolations
-! of lower orders
-!
-    u(  1) = sum(ce2(:) * fi(  1:  2))
-    u(  2) = sum(ce3(:) * fi(  1:  3))
-    u(  3) = sum(ce5(:) * fi(  1:  5))
-    u(n-3) = sum(ce5(:) * fi(n-5:n-1))
-    u(n-2) = sum(ce5(:) * fi(n-4:n  ))
-    u(n-1) = sum(ce3(:) * fi(n-2:n  ))
-    u(n  ) =              fi(n      )
-
-! apply the monotonicity preserving limiting
-!
-    call mp_limiting(fi(:), u(:))
-
-! copy the interpolation to the respective vector
-!
-    fr(1:n-1) = u(n-1:1:-1)
-
-! update the interpolation of the first and last points
-!
-    i     = n - 1
-    fl(1) = 0.5d+00 * (fc(1) + fc(2))
-    fr(i) = 0.5d+00 * (fc(i) + fc(n))
-    fl(n) = fc(n)
-    fr(n) = fc(n)
-
-!-------------------------------------------------------------------------------
-!
-  end subroutine reconstruct_mp7ld
-!
-!===============================================================================
-!
-! subroutine RECONSTRUCT_MP9LD:
-! ----------------------------
-!
-!   Subroutine reconstructs the interface states using the ninth order
-!   low dissipation Monotonicity Preserving (MP) method.
-!
-!   Arguments are described in subroutine reconstruct().
-!
-!   References:
-!
-!     [1] Suresh, A. & Huynh, H. T.,
-!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
-!          Time Stepping"
-!         Journal on Computational Physics,
-!         1997, vol. 136, pp. 83-99,
-!         http://dx.doi.org/10.1006/jcph.1997.5745
-!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
-!         "A 5th order monotonicity-preserving upwind compact difference
-!          scheme",
-!         Science China Physics, Mechanics and Astronomy,
-!         Volume 54, Issue 3, pp. 511-522,
-!         http://dx.doi.org/10.1007/s11433-010-4220-x
-!
-!===============================================================================
-!
-  subroutine reconstruct_mp9ld(h, fc, fl, fr)
-
-! local variables are not implicit by default
-!
-    implicit none
-
-! subroutine arguments
-!
-    real(kind=8)              , intent(in)  :: h
-    real(kind=8), dimension(:), intent(in)  :: fc
-    real(kind=8), dimension(:), intent(out) :: fl, fr
-
-! local variables
-!
-    integer :: n, i
-
-! local arrays for derivatives
-!
-    real(kind=8), dimension(size(fc)) :: fi
-    real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
-    real(kind=8), dimension(10), parameter ::                                  &
-           ce9 = [ 4.0000d+01,-4.1500d+02, 2.0450d+03,-6.7150d+03,             &
-                   2.0165d+04, 1.5629d+04,-3.6910d+03, 7.4900d+02,             &
-                  -9.1000d+01, 4.0000d+00 ] / 2.7720d+04
-    real(kind=8), dimension(7), parameter ::                                   &
-           ce7 = [-3.000d+00, 2.500d+01,-1.010d+02, 3.190d+02, 2.140d+02,      &
-                  -3.800d+01, 4.000d+00 ] / 4.200d+02
-    real(kind=8), dimension(5), parameter ::                                   &
-           ce5 = [ 2.0d+00,-1.3d+01, 4.7d+01, 2.7d+01,-3.0d+00 ] / 6.0d+01
-    real(kind=8), dimension(3), parameter ::                                   &
-           ce3 = [-1.0d+00, 5.0d+00, 2.0d+00 ] / 6.0d+00
-    real(kind=8), dimension(2), parameter ::                                   &
-           ce2 = [ 5.0d-01, 5.0d-01 ]
-!
-!-------------------------------------------------------------------------------
-!
-! get the input vector length
-!
-    n = size(fc)
-
-!! === left-side interpolation ===
-!!
-! reconstruct the interface state using the 9th order interpolation
-!
-    do i = 5, n - 5
-      u(i) = sum(ce9(:) * fc(i-4:i+5))
-    end do
-
-! interpolate the interface state of the ghost zones using the interpolations
-! of lower orders
-!
-    u(  1) = sum(ce2(:) * fc(  1:  2))
-    u(  2) = sum(ce3(:) * fc(  1:  3))
-    u(  3) = sum(ce5(:) * fc(  1:  5))
-    u(  4) = sum(ce7(:) * fc(  1:  7))
-    u(n-4) = sum(ce7(:) * fc(n-7:n-1))
-    u(n-3) = sum(ce7(:) * fc(n-6:n  ))
-    u(n-2) = sum(ce5(:) * fc(n-4:n  ))
-    u(n-1) = sum(ce3(:) * fc(n-2:n  ))
-    u(n  ) =              fc(n      )
-
-! apply the monotonicity preserving limiting
-!
-    call mp_limiting(fc(:), u(:))
-
-! copy the interpolation to the respective vector
-!
-    fl(1:n) = u(1:n)
-
-!! === right-side interpolation ===
-!!
-! invert the cell-centered value vector
-!
-    fi(1:n) = fc(n:1:-1)
-
-! reconstruct the interface state using the 9th order interpolation
-!
-    do i = 5, n - 5
-      u(i) = sum(ce9(:) * fi(i-4:i+5))
-    end do
-
-! interpolate the interface state of the ghost zones using the interpolations
-! of lower orders
-!
-    u(  1) = sum(ce2(:) * fi(  1:  2))
-    u(  2) = sum(ce3(:) * fi(  1:  3))
-    u(  3) = sum(ce5(:) * fi(  1:  5))
-    u(  4) = sum(ce7(:) * fi(  1:  7))
-    u(n-4) = sum(ce7(:) * fi(n-7:n-1))
-    u(n-3) = sum(ce7(:) * fi(n-6:n  ))
-    u(n-2) = sum(ce5(:) * fi(n-4:n  ))
-    u(n-1) = sum(ce3(:) * fi(n-2:n  ))
-    u(n  ) =              fi(n      )
-
-! apply the monotonicity preserving limiting
-!
-    call mp_limiting(fi(:), u(:))
-
-! copy the interpolation to the respective vector
-!
-    fr(1:n-1) = u(n-1:1:-1)
-
-! update the interpolation of the first and last points
-!
-    i     = n - 1
-    fl(1) = 0.5d+00 * (fc(1) + fc(2))
-    fr(i) = 0.5d+00 * (fc(i) + fc(n))
-    fl(n) = fc(n)
-    fr(n) = fc(n)
-
-!-------------------------------------------------------------------------------
-!
-  end subroutine reconstruct_mp9ld
-!
-!===============================================================================
-!
 ! subroutine RECONSTRUCT_CRMP5:
 ! ----------------------------
 !
-!   Subroutine reconstructs the interface states using the fifth order
-!   Compact Reconstruction Monotonicity Preserving (CRMP) method.
+!   Subroutine reconstructs the interface states using the 5th order
+!   Compact Reconstruction Monotonicity Preserving (CRMP) method with
+!   adjustable dissipation.
 !
 !   Arguments are described in subroutine reconstruct().
 !
@@ -4619,72 +4105,45 @@ module interpolations
 !         Journal on Computational Physics,
 !         1997, vol. 136, pp. 83-99,
 !         http://dx.doi.org/10.1006/jcph.1997.5745
-!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
-!         "A 5th order monotonicity-preserving upwind compact difference
-!          scheme",
-!         Science China Physics, Mechanics and Astronomy,
-!         Volume 54, Issue 3, pp. 511-522,
-!         http://dx.doi.org/10.1007/s11433-010-4220-x
+!     [2] Jian Fang, Yufeng Yao, Zhaorui Li, Lipeng Lu,
+!         "Investigation of low-dissipation monotonicity-preserving scheme
+!          for direct numerical simulation of compressible turbulent flows",
+!         Computers & Fluids,
+!         2014, vol. 104, pp. 55-72,
+!         http://dx.doi.org/10.1016/j.compfluid.2014.07.024
 !
 !===============================================================================
 !
   subroutine reconstruct_crmp5(h, fc, fl, fr)
 
-! include external procedures
-!
-    use algebra   , only : tridiag
+    use algebra, only : tridiag
 
-! local variables are not implicit by default
-!
     implicit none
 
-! subroutine arguments
-!
     real(kind=8)              , intent(in)  :: h
     real(kind=8), dimension(:), intent(in)  :: fc
     real(kind=8), dimension(:), intent(out) :: fl, fr
 
-! local variables
-!
     integer :: n, i
 
-! local arrays for derivatives
-!
     real(kind=8), dimension(size(fc)) :: fi
     real(kind=8), dimension(size(fc)) :: a, b, c
     real(kind=8), dimension(size(fc)) :: r
     real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
-    real(kind=8), dimension(3), parameter ::                                   &
-                                di  = (/ 3.0d-01, 6.0d-01, 1.0d-01 /)
-    real(kind=8), dimension(3), parameter ::                                   &
-                                ci5 = (/ 1.0d+00, 1.9d+01, 1.0d+01 /) / 3.0d+01
-    real(kind=8), dimension(5), parameter ::                                   &
-                                ce5 = (/ 2.0d+00,-1.3d+01, 4.7d+01             &
-                                                , 2.7d+01,-3.0d+00 /) / 6.0d+01
-    real(kind=8), dimension(3), parameter :: &
-                                ce3 = (/-1.0d+00, 5.0d+00, 2.0d+00 /) / 6.0d+00
-    real(kind=8), dimension(2), parameter :: ce2 = (/ 5.0d-01, 5.0d-01 /)
 !
 !-------------------------------------------------------------------------------
-!
-! get the input vector length
 !
     n = size(fc)
 
-! prepare the diagonals of the tridiagonal matrix
-!
     do i = 1, ng
       a(i) = 0.0d+00
       b(i) = 1.0d+00
       c(i) = 0.0d+00
     end do
     do i = ng + 1, n - ng - 1
-      a(i) = di(1)
-      b(i) = di(2)
-      c(i) = di(3)
+      a(i) = di5a(1)
+      b(i) = di5a(2)
+      c(i) = di5a(3)
     end do
     do i = n - ng, n
       a(i) = 0.0d+00
@@ -4694,14 +4153,10 @@ module interpolations
 
 !! === left-side interpolation ===
 !!
-! prepare the tridiagonal system coefficients for the interior part
-!
     do i = ng, n - ng + 1
-      r(i) = sum(ci5(:) * fc(i-1:i+1))
+      r(i) = sum(ci5a(:) * fc(i-1:i+2))
     end do
 
-! interpolate ghost zones using the explicit method
-!
     r(  1) = sum(ce2(:) * fc(  1:  2))
     r(  2) = sum(ce3(:) * fc(  1:  3))
     do i = 3, ng
@@ -4713,32 +4168,20 @@ module interpolations
     r(n-1) = sum(ce3(:) * fc(n-2:  n))
     r(n  ) =              fc(n      )
 
-! solve the tridiagonal system of equations
-!
     call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
 
-! apply the monotonicity preserving limiting
-!
     call mp_limiting(fc(:), u(:))
 
-! copy the interpolation to the respective vector
-!
     fl(1:n) = u(1:n)
 
 !! === right-side interpolation ===
 !!
-! invert the cell-centered value vector
-!
     fi(1:n) = fc(n:1:-1)
 
-! prepare the tridiagonal system coefficients for the interior part
-!
     do i = ng, n - ng + 1
-      r(i) = sum(ci5(:) * fi(i-1:i+1))
+      r(i) = sum(ci5a(:) * fi(i-1:i+2))
     end do ! i = ng, n - ng + 1
 
-! interpolate ghost zones using the explicit method
-!
     r(  1) = sum(ce2(:) * fi(  1:  2))
     r(  2) = sum(ce3(:) * fi(  1:  3))
     do i = 3, ng
@@ -4750,20 +4193,12 @@ module interpolations
     r(n-1) = sum(ce3(:) * fi(n-2:  n))
     r(n  ) =              fi(n      )
 
-! solve the tridiagonal system of equations
-!
     call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
 
-! apply the monotonicity preserving limiting
-!
     call mp_limiting(fi(:), u(:))
 
-! copy the interpolation to the respective vector
-!
     fr(1:n-1) = u(n-1:1:-1)
 
-! update the interpolation of the first and last points
-!
     i     = n - 1
     fl(1) = 0.5d+00 * (fc(1) + fc(2))
     fr(i) = 0.5d+00 * (fc(i) + fc(n))
@@ -4776,273 +4211,51 @@ module interpolations
 !
 !===============================================================================
 !
-! subroutine RECONSTRUCT_CRMP5LD:
-! ------------------------------
-!
-!   Subroutine reconstructs the interface states using the fifth order
-!   Low-Dissipation Compact Reconstruction Monotonicity Preserving (CRMP)
-!   method.
-!
-!   Arguments are described in subroutine reconstruct().
-!
-!   References:
-!
-!     [1] Suresh, A. & Huynh, H. T.,
-!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
-!          Time Stepping"
-!         Journal on Computational Physics,
-!         1997, vol. 136, pp. 83-99,
-!         http://dx.doi.org/10.1006/jcph.1997.5745
-!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
-!         "A 5th order monotonicity-preserving upwind compact difference
-!          scheme",
-!         Science China Physics, Mechanics and Astronomy,
-!         Volume 54, Issue 3, pp. 511-522,
-!         http://dx.doi.org/10.1007/s11433-010-4220-x
-!     [3] Ghosh, D. & Baeder, J.,
-!         "Compact Reconstruction Schemes With Weighted ENO Limiting For
-!          Hyperbolic Conservation Laws",
-!         SIAM Journal on Scientific Computing,
-!         2012, vol. 34, no. 3, pp. A1678-A1705,
-!         http://dx.doi.org/10.1137/110857659
-!
-!===============================================================================
-!
-  subroutine reconstruct_crmp5ld(h, fc, fl, fr)
-
-! include external procedures
-!
-    use algebra   , only : tridiag
-
-! local variables are not implicit by default
-!
-    implicit none
-
-! subroutine arguments
-!
-    real(kind=8)              , intent(in)  :: h
-    real(kind=8), dimension(:), intent(in)  :: fc
-    real(kind=8), dimension(:), intent(out) :: fl, fr
-
-! local variables
-!
-    integer :: n, i
-
-! local arrays for derivatives
-!
-    real(kind=8), dimension(size(fc)) :: fi
-    real(kind=8), dimension(size(fc)) :: a, b, c
-    real(kind=8), dimension(size(fc)) :: r
-    real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
-    real(kind=8), dimension(3), parameter ::                                   &
-                                di  = (/ 2.5d-01, 6.0d-01, 1.5d-01 /)
-    real(kind=8), dimension(4), parameter ::                                   &
-                                ci5 = (/ 3.0d+00, 6.7d+01, 4.9d+01             &
-                                                         , 1.0d+00 /) / 1.2d+02
-    real(kind=8), dimension(5), parameter ::                                   &
-                                ce5 = (/ 2.0d+00,-1.3d+01, 4.7d+01             &
-                                                , 2.7d+01,-3.0d+00 /) / 6.0d+01
-    real(kind=8), dimension(3), parameter :: &
-                                ce3 = (/-1.0d+00, 5.0d+00, 2.0d+00 /) / 6.0d+00
-    real(kind=8), dimension(2), parameter :: ce2 = (/ 5.0d-01, 5.0d-01 /)
-!
-!-------------------------------------------------------------------------------
-!
-! get the input vector length
-!
-    n = size(fc)
-
-! prepare the diagonals of the tridiagonal matrix
-!
-    do i = 1, ng
-      a(i) = 0.0d+00
-      b(i) = 1.0d+00
-      c(i) = 0.0d+00
-    end do
-    do i = ng + 1, n - ng - 1
-      a(i) = di(1)
-      b(i) = di(2)
-      c(i) = di(3)
-    end do
-    do i = n - ng, n
-      a(i) = 0.0d+00
-      b(i) = 1.0d+00
-      c(i) = 0.0d+00
-    end do
-
-!! === left-side interpolation ===
-!!
-! prepare the tridiagonal system coefficients for the interior part
-!
-    do i = ng, n - ng + 1
-      r(i) = sum(ci5(:) * fc(i-1:i+2))
-    end do
-
-! interpolate ghost zones using the explicit method
-!
-    r(  1) = sum(ce2(:) * fc(  1:  2))
-    r(  2) = sum(ce3(:) * fc(  1:  3))
-    do i = 3, ng
-      r(i) = sum(ce5(:) * fc(i-2:i+2))
-    end do
-    do i = n - ng, n - 2
-      r(i) = sum(ce5(:) * fc(i-2:i+2))
-    end do
-    r(n-1) = sum(ce3(:) * fc(n-2:  n))
-    r(n  ) =              fc(n      )
-
-! solve the tridiagonal system of equations
-!
-    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
-
-! apply the monotonicity preserving limiting
-!
-    call mp_limiting(fc(:), u(:))
-
-! copy the interpolation to the respective vector
-!
-    fl(1:n) = u(1:n)
-
-!! === right-side interpolation ===
-!!
-! invert the cell-centered value vector
-!
-    fi(1:n) = fc(n:1:-1)
-
-! prepare the tridiagonal system coefficients for the interior part
-!
-    do i = ng, n - ng + 1
-      r(i) = sum(ci5(:) * fi(i-1:i+2))
-    end do ! i = ng, n - ng + 1
-
-! interpolate ghost zones using the explicit method
-!
-    r(  1) = sum(ce2(:) * fi(  1:  2))
-    r(  2) = sum(ce3(:) * fi(  1:  3))
-    do i = 3, ng
-      r(i) = sum(ce5(:) * fi(i-2:i+2))
-    end do
-    do i = n - ng, n - 2
-      r(i) = sum(ce5(:) * fi(i-2:i+2))
-    end do
-    r(n-1) = sum(ce3(:) * fi(n-2:  n))
-    r(n  ) =              fi(n      )
-
-! solve the tridiagonal system of equations
-!
-    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
-
-! apply the monotonicity preserving limiting
-!
-    call mp_limiting(fi(:), u(:))
-
-! copy the interpolation to the respective vector
-!
-    fr(1:n-1) = u(n-1:1:-1)
-
-! update the interpolation of the first and last points
-!
-    i     = n - 1
-    fl(1) = 0.5d+00 * (fc(1) + fc(2))
-    fr(i) = 0.5d+00 * (fc(i) + fc(n))
-    fl(n) = fc(n)
-    fr(n) = fc(n)
-
-!-------------------------------------------------------------------------------
-!
-  end subroutine reconstruct_crmp5ld
-!
-!===============================================================================
-!
 ! subroutine RECONSTRUCT_CRMP7:
 ! ----------------------------
 !
-!   Subroutine reconstructs the interface states using the seventh order
-!   Compact Reconstruction Monotonicity Preserving (CRMP) method.
+!   Subroutine reconstructs the interface states using the 7th order
+!   Compact Reconstruction Monotonicity Preserving (CRMP) method with
+!   adjustable dissipation.
 !
 !   Arguments are described in subroutine reconstruct().
 !
 !   References:
 !
-!     [1] Suresh, A. & Huynh, H. T.,
-!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
-!          Time Stepping"
-!         Journal on Computational Physics,
-!         1997, vol. 136, pp. 83-99,
-!         http://dx.doi.org/10.1006/jcph.1997.5745
-!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
-!         "A 5th order monotonicity-preserving upwind compact difference
-!          scheme",
-!         Science China Physics, Mechanics and Astronomy,
-!         Volume 54, Issue 3, pp. 511-522,
-!         http://dx.doi.org/10.1007/s11433-010-4220-x
+!     - See RECONSTRUCT_CRMP5
 !
 !===============================================================================
 !
   subroutine reconstruct_crmp7(h, fc, fl, fr)
 
-! include external procedures
-!
-    use algebra   , only : tridiag
+    use algebra, only : tridiag
 
-! local variables are not implicit by default
-!
     implicit none
 
-! subroutine arguments
-!
     real(kind=8)              , intent(in)  :: h
     real(kind=8), dimension(:), intent(in)  :: fc
     real(kind=8), dimension(:), intent(out) :: fl, fr
 
-! local variables
-!
     integer :: n, i
 
-! local arrays for derivatives
-!
     real(kind=8), dimension(size(fc)) :: fi
     real(kind=8), dimension(size(fc)) :: a, b, c
     real(kind=8), dimension(size(fc)) :: r
     real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
-    real(kind=8), dimension(3), parameter ::                                   &
-                                di  = (/ 2.0d+00, 4.0d+00, 1.0d+00/) / 7.0d+00
-    real(kind=8), dimension(5), parameter ::                                   &
-                                ci  = (/-1.0d+00, 1.9d+01, 2.39d+02            &
-                                              , 1.59d+02, 4.0d+00 /) / 4.2d+02
-    real(kind=8), dimension(7), parameter ::                                   &
-                                ce7 = (/-3.0d+00, 2.5d+01,-1.01d+02, 3.19d+02  &
-                                      , 2.14d+02,-3.8d+01, 4.0d+00/) / 4.2d+02
-    real(kind=8), dimension(5), parameter ::                                   &
-                                ce5 = (/ 2.0d+00,-1.3d+01, 4.70d+01            &
-                                               , 2.70d+01,-3.0d+00/) / 6.0d+01
-    real(kind=8), dimension(3), parameter :: &
-                                ce3 = (/-1.0d+00, 5.0d+00, 2.0d+00/) / 6.0d+00
-    real(kind=8), dimension(2), parameter :: ce2 = (/ 5.0d-01, 5.0d-01 /)
 !
 !-------------------------------------------------------------------------------
-!
-! get the input vector length
 !
     n = size(fc)
 
-! prepare the diagonals of the tridiagonal matrix
-!
     do i = 1, ng
       a(i) = 0.0d+00
       b(i) = 1.0d+00
       c(i) = 0.0d+00
     end do
     do i = ng + 1, n - ng - 1
-      a(i) = di(1)
-      b(i) = di(2)
-      c(i) = di(3)
+      a(i) = di7a(1)
+      b(i) = di7a(2)
+      c(i) = di7a(3)
     end do
     do i = n - ng, n
       a(i) = 0.0d+00
@@ -5052,14 +4265,10 @@ module interpolations
 
 !! === left-side interpolation ===
 !!
-! prepare the tridiagonal system coefficients for the interior part
-!
     do i = ng, n - ng + 1
-      r(i) = sum( ci(:) * fc(i-2:i+2))
+      r(i) = sum(ci7a(:) * fc(i-2:i+3))
     end do
 
-! interpolate ghost zones using the explicit method
-!
     r(  1) = sum(ce2(:) * fc(  1:  2))
     r(  2) = sum(ce3(:) * fc(  1:  3))
     r(  3) = sum(ce5(:) * fc(  1:  5))
@@ -5073,32 +4282,20 @@ module interpolations
     r(n-1) = sum(ce3(:) * fc(n-2:  n))
     r(n  ) =              fc(n      )
 
-! solve the tridiagonal system of equations
-!
     call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
 
-! apply the monotonicity preserving limiting
-!
     call mp_limiting(fc(:), u(:))
 
-! copy the interpolation to the respective vector
-!
     fl(1:n) = u(1:n)
 
 !! === right-side interpolation ===
 !!
-! invert the cell-centered value vector
-!
     fi(1:n) = fc(n:1:-1)
 
-! prepare the tridiagonal system coefficients for the interior part
-!
     do i = ng, n - ng + 1
-      r(i) = sum( ci(:) * fi(i-2:i+2))
+      r(i) = sum(ci7a(:) * fi(i-2:i+3))
     end do ! i = ng, n - ng + 1
 
-! interpolate ghost zones using the explicit method
-!
     r(  1) = sum(ce2(:) * fi(  1:  2))
     r(  2) = sum(ce3(:) * fi(  1:  3))
     r(  3) = sum(ce5(:) * fi(  1:  5))
@@ -5112,20 +4309,12 @@ module interpolations
     r(n-1) = sum(ce3(:) * fi(n-2:  n))
     r(n  ) =              fi(n      )
 
-! solve the tridiagonal system of equations
-!
     call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
 
-! apply the monotonicity preserving limiting
-!
     call mp_limiting(fi(:), u(:))
 
-! copy the interpolation to the respective vector
-!
     fr(1:n-1) = u(n-1:1:-1)
 
-! update the interpolation of the first and last points
-!
     i     = n - 1
     fl(1) = 0.5d+00 * (fc(1) + fc(2))
     fr(i) = 0.5d+00 * (fc(i) + fc(n))
@@ -5138,177 +4327,120 @@ module interpolations
 !
 !===============================================================================
 !
-! subroutine RECONSTRUCT_CRMP7LD:
-! ------------------------------
+! subroutine RECONSTRUCT_CRMP9:
+! ----------------------------
 !
-!   Subroutine reconstructs the interface states using the seventh order
-!   low dissipation Compact Reconstruction Monotonicity Preserving (CRMP)
-!   method.
+!   Subroutine reconstructs the interface states using the 9th order
+!   Compact Reconstruction Monotonicity Preserving (CRMP) method with
+!   adjustable dissipation.
 !
 !   Arguments are described in subroutine reconstruct().
 !
 !   References:
 !
-!     [1] Suresh, A. & Huynh, H. T.,
-!         "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
-!          Time Stepping"
-!         Journal on Computational Physics,
-!         1997, vol. 136, pp. 83-99,
-!         http://dx.doi.org/10.1006/jcph.1997.5745
-!     [2] He, ZhiWei, Li, XinLiang, Fu, DeXun, & Ma, YanWen,
-!         "A 5th order monotonicity-preserving upwind compact difference
-!          scheme",
-!         Science China Physics, Mechanics and Astronomy,
-!         Volume 54, Issue 3, pp. 511-522,
-!         http://dx.doi.org/10.1007/s11433-010-4220-x
+!     - See RECONSTRUCT_CRMP5
 !
 !===============================================================================
 !
-  subroutine reconstruct_crmp7ld(h, fc, fl, fr)
+  subroutine reconstruct_crmp9(h, fc, fl, fr)
 
-! include external procedures
-!
-    use algebra   , only : tridiag
+    use algebra, only : pentadiag
 
-! local variables are not implicit by default
-!
     implicit none
 
-! subroutine arguments
-!
     real(kind=8)              , intent(in)  :: h
     real(kind=8), dimension(:), intent(in)  :: fc
     real(kind=8), dimension(:), intent(out) :: fl, fr
 
-! local variables
-!
     integer :: n, i
 
-! local arrays for derivatives
-!
     real(kind=8), dimension(size(fc)) :: fi
-    real(kind=8), dimension(size(fc)) :: a, b, c
+    real(kind=8), dimension(size(fc)) :: e, c, d, a, b
     real(kind=8), dimension(size(fc)) :: r
     real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
-    real(kind=8), dimension(3), parameter ::                                   &
-                              di  = [ 5.0d+00, 1.2d+01, 4.0d+00 ] / 2.1d+01
-    real(kind=8), dimension(6), parameter ::                                   &
-                              ci  = [-2.00d+00, 4.20d+01, 6.37d+02,            &
-                                      5.57d+02, 2.70d+01,-1.00d+00 ] / 1.26d+03
-    real(kind=8), dimension(7), parameter ::                                   &
-                              ce7 = [-3.0d+00, 2.5d+01,-1.01d+02, 3.19d+02,    &
-                                      2.14d+02,-3.8d+01, 4.0d+00 ] / 4.2d+02
-    real(kind=8), dimension(5), parameter ::                                   &
-                              ce5 = [ 2.0d+00,-1.3d+01, 4.70d+01,              &
-                                               2.70d+01,-3.0d+00 ] / 6.0d+01
-    real(kind=8), dimension(3), parameter ::                                   &
-                              ce3 = [-1.0d+00, 5.0d+00, 2.0d+00 ] / 6.0d+00
-    real(kind=8), dimension(2), parameter ::                                   &
-                              ce2 = [ 5.0d-01, 5.0d-01 ]
 !
 !-------------------------------------------------------------------------------
-!
-! get the input vector length
 !
     n = size(fc)
 
-! prepare the diagonals of the tridiagonal matrix
-!
     do i = 1, ng
-      a(i) = 0.0d+00
-      b(i) = 1.0d+00
+      e(i) = 0.0d+00
       c(i) = 0.0d+00
+      d(i) = 1.0d+00
+      a(i) = 0.0d+00
+      b(i) = 0.0d+00
     end do
     do i = ng + 1, n - ng - 1
-      a(i) = di(1)
-      b(i) = di(2)
-      c(i) = di(3)
+      e(i) = di9a(1)
+      c(i) = di9a(2)
+      d(i) = di9a(3)
+      a(i) = di9a(4)
+      b(i) = di9a(5)
     end do
     do i = n - ng, n
-      a(i) = 0.0d+00
-      b(i) = 1.0d+00
+      e(i) = 0.0d+00
       c(i) = 0.0d+00
+      d(i) = 1.0d+00
+      a(i) = 0.0d+00
+      b(i) = 0.0d+00
     end do
 
 !! === left-side interpolation ===
 !!
-! prepare the tridiagonal system coefficients for the interior part
-!
     do i = ng, n - ng + 1
-      r(i) = sum( ci(:) * fc(i-2:i+3))
+      r(i) = sum(ci9a(:) * fc(i-2:i+3))
     end do
 
-! interpolate ghost zones using the explicit method
-!
     r(  1) = sum(ce2(:) * fc(  1:  2))
     r(  2) = sum(ce3(:) * fc(  1:  3))
     r(  3) = sum(ce5(:) * fc(  1:  5))
-    do i = 4, ng
-      r(i) = sum(ce7(:) * fc(i-3:i+3))
+    r(  4) = sum(ce7(:) * fc(  1:  7))
+    do i = 5, ng
+      r(i) = sum(ce9(:) * fc(i-4:i+4))
     end do
-    do i = n - ng, n - 3
-      r(i) = sum(ce7(:) * fc(i-3:i+3))
+    do i = n - ng, n - 4
+      r(i) = sum(ce9(:) * fc(i-4:i+4))
     end do
+    r(n-3) = sum(ce7(:) * fc(n-6:  n))
     r(n-2) = sum(ce5(:) * fc(n-4:  n))
     r(n-1) = sum(ce3(:) * fc(n-2:  n))
     r(n  ) =              fc(n      )
 
-! solve the tridiagonal system of equations
-!
-    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
+    call pentadiag(n, e(:), c(:), d(:), a(:), b(:), r(:), u(:))
 
-! apply the monotonicity preserving limiting
-!
     call mp_limiting(fc(:), u(:))
 
-! copy the interpolation to the respective vector
-!
     fl(1:n) = u(1:n)
 
 !! === right-side interpolation ===
 !!
-! invert the cell-centered value vector
-!
     fi(1:n) = fc(n:1:-1)
 
-! prepare the tridiagonal system coefficients for the interior part
-!
     do i = ng, n - ng + 1
-      r(i) = sum( ci(:) * fi(i-2:i+3))
+      r(i) = sum(ci9a(:) * fi(i-2:i+3))
     end do ! i = ng, n - ng + 1
 
-! interpolate ghost zones using the explicit method
-!
     r(  1) = sum(ce2(:) * fi(  1:  2))
     r(  2) = sum(ce3(:) * fi(  1:  3))
     r(  3) = sum(ce5(:) * fi(  1:  5))
-    do i = 4, ng
-      r(i) = sum(ce7(:) * fi(i-3:i+3))
+    r(  4) = sum(ce7(:) * fi(  1:  7))
+    do i = 5, ng
+      r(i) = sum(ce9(:) * fi(i-4:i+4))
     end do
-    do i = n - ng, n - 3
-      r(i) = sum(ce7(:) * fi(i-3:i+3))
+    do i = n - ng, n - 4
+      r(i) = sum(ce9(:) * fi(i-4:i+4))
     end do
+    r(n-3) = sum(ce7(:) * fi(n-6:  n))
     r(n-2) = sum(ce5(:) * fi(n-4:  n))
     r(n-1) = sum(ce3(:) * fi(n-2:  n))
     r(n  ) =              fi(n      )
 
-! solve the tridiagonal system of equations
-!
-    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
+    call pentadiag(n, e(:), c(:), d(:), a(:), b(:), r(:), u(:))
 
-! apply the monotonicity preserving limiting
-!
     call mp_limiting(fi(:), u(:))
 
-! copy the interpolation to the respective vector
-!
     fr(1:n-1) = u(n-1:1:-1)
 
-! update the interpolation of the first and last points
-!
     i     = n - 1
     fl(1) = 0.5d+00 * (fc(1) + fc(2))
     fr(i) = 0.5d+00 * (fc(i) + fc(n))
@@ -5317,7 +4449,7 @@ module interpolations
 
 !-------------------------------------------------------------------------------
 !
-  end subroutine reconstruct_crmp7ld
+  end subroutine reconstruct_crmp9
 !
 !===============================================================================
 !
@@ -5360,14 +4492,11 @@ module interpolations
     real(kind=8), dimension(3), parameter ::                                   &
             di  = [ 5.0163016d-01, 1.0d+00, 2.5394716d-01 ]
     real(kind=8), dimension(5), parameter ::                                   &
-            ci5 = [-4.44553d-03,8.101861d-02, 9.9428149d-01,                   &
-                                6.6721233d-01, 1.751043d-02 ]
-    real(kind=8), dimension(5), parameter ::                                   &
-            ce5 = [ 2.0d+00,-1.3d+01, 4.7d+01, 2.7d+01,-3.0d+00 ] / 6.0d+01
-    real(kind=8), dimension(3), parameter ::                                   &
-            ce3 = [-1.0d+00, 5.0d+00, 2.0d+00 ] / 6.0d+00
-    real(kind=8), dimension(2), parameter ::                                   &
-            ce2 = [ 5.0d-01, 5.0d-01 ]
+            ci5 = [- 3.0d+00 * di(1) - 3.0d+00 * di(3) + 2.0d+00,              &
+                     2.7d+01 * di(1) + 1.7d+01 * di(3) - 1.3d+01,              &
+                     4.7d+01 * di(1) - 4.3d+01 * di(3) + 4.7d+01,              &
+                   - 1.3d+01 * di(1) + 7.7d+01 * di(3) + 2.7d+01,              &
+                     2.0d+00 * di(1) + 1.2d+01 * di(3) - 3.0d+00 ] / 6.0d+01
 !
 !-------------------------------------------------------------------------------
 !
@@ -5625,15 +4754,6 @@ module interpolations
 !
     real(kind=8), dimension(size(fc)) :: fi
     real(kind=8), dimension(size(fc)) :: u
-
-! local parameters
-!
-    real(kind=8), dimension(5), parameter ::                                   &
-                                ce5 = (/ 2.0d+00,-1.3d+01, 4.70d+01            &
-                                               , 2.70d+01,-3.0d+00/) / 6.0d+01
-    real(kind=8), dimension(3), parameter :: &
-                                ce3 = (/-1.0d+00, 5.0d+00, 2.0d+00/) / 6.0d+00
-    real(kind=8), dimension(2), parameter :: ce2 = (/ 5.0d-01, 5.0d-01 /)
 !
 !-------------------------------------------------------------------------------
 !