diff --git a/sources/interpolations.F90 b/sources/interpolations.F90
index 6c87fa1..3485b9c 100644
--- a/sources/interpolations.F90
+++ b/sources/interpolations.F90
@@ -354,6 +354,13 @@ module interpolations
       reconstruct_states => reconstruct_ocmp7t
       order              = 7
       nghosts            = max(nghosts, 6)
+    case ("ocmp9t", "OCMP9t")
+      name_rec           =  "9th order Tridiagonal Optimized " // &
+                            "Compact Monotonicity Preserving"
+      interfaces         => interfaces_dir
+      reconstruct_states => reconstruct_ocmp9t
+      order              = 7
+      nghosts            = max(nghosts, 6)
     case ("ocmp7", "OCMP7", "ocmp7p", "OCMP7P")
       name_rec           =  "7th order Pentadiagonal Optimized " // &
                             "Compact Monotonicity Preserving"
@@ -4528,6 +4535,141 @@ module interpolations
 !
 !===============================================================================
 !
+! subroutine RECONSTRUCT_OCMP9T:
+! -----------------------------
+!
+!   Subroutine reconstructs the interface states using the 9th order Optimized
+!   Compact Reconstruction Monotonicity Preserving (CRMP) method.
+!
+!   The method uses an implicit tri-diagonal solver and its coefficients are
+!   optimized for the maximum value of imaginary part of the scaled modified
+!   wave number (Im(κ)) to be 2e-6 and the integrated error of compact
+!   discretization E to be 1.32e-11.
+!
+!   Arguments are described in subroutine reconstruct().
+!
+!===============================================================================
+!
+  subroutine reconstruct_ocmp9t(p, h, fc, fl, fr)
+
+    use algebra, only : tridiag
+
+    implicit none
+
+    logical                   , intent(in)  :: p
+    real(kind=8)              , intent(in)  :: h
+    real(kind=8), dimension(:), intent(in)  :: fc
+    real(kind=8), dimension(:), intent(out) :: fl, fr
+
+    integer :: n, i
+
+    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
+
+    real(kind=8), parameter :: a1 = 5.0087697436175480833840357d-01
+    real(kind=8), parameter :: a2 = 3.4091832657841234002695365d-01
+    real(kind=8), dimension(3), parameter :: di9 = [ a1, 1.0d+00, a2 ]
+    real(kind=8), dimension(9), parameter :: &
+             ci9 = [ - 5.000d+00 * a1 - 5.000d+00 * a2 + 4.000d+00, &
+                       5.500d+01 * a1 + 4.900d+01 * a2 - 4.100d+01, &
+                     - 3.050d+02 * a1 - 2.210d+02 * a2 + 1.990d+02, &
+                       1.375d+03 * a1 + 6.190d+02 * a2 - 6.410d+02, &
+                       1.879d+03 * a1 - 1.271d+03 * a2 + 1.879d+03, &
+                     - 6.410d+02 * a1 + 2.509d+03 * a2 + 1.375d+03, &
+                       1.990d+02 * a1 + 9.550d+02 * a2 - 3.050d+02, &
+                     - 4.100d+01 * a1 - 1.250d+02 * a2 + 5.500d+01, &
+                       4.000d+00 * a1 + 1.000d+01 * a2 - 5.000d+00 ] / 2.52d+03
+
+!-------------------------------------------------------------------------------
+!
+    n = size(fc)
+
+    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) = di9(1)
+      b(i) = di9(2)
+      c(i) = di9(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 ===
+!!
+    do i = ng, n - ng + 1
+      r(i) = sum(ci9(:) * fc(i-4:i+4))
+    end do
+
+    r(  1) = sum(ce2(:) * fc(  1:  2))
+    r(  2) = sum(ce3(:) * fc(  1:  3))
+    r(  3) = sum(ce5(:) * fc(  1:  5))
+    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 - 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      )
+
+    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
+
+    call mp_limiting(p, n, fc(:), u(:))
+
+    fl(1:n) = u(1:n)
+
+!! === right-side interpolation ===
+!!
+    fi(1:n) = fc(n:1:-1)
+
+    do i = ng, n - ng + 1
+      r(i) = sum(ci9(:) * fi(i-4:i+4))
+    end do
+
+    r(  1) = sum(ce2(:) * fi(  1:  2))
+    r(  2) = sum(ce3(:) * fi(  1:  3))
+    r(  3) = sum(ce5(:) * fi(  1:  5))
+    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 - 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      )
+
+    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
+
+    call mp_limiting(p, n, fi(:), u(:))
+
+    fr(1:n-1) = u(n-1:1:-1)
+
+    i     = n - 1
+    fl(1) = 5.0d-01 * (fc(1) + fc(2))
+    fr(i) = 5.0d-01 * (fc(i) + fc(n))
+    fl(n) = fc(n)
+    fr(n) = fc(n)
+
+!-------------------------------------------------------------------------------
+!
+  end subroutine reconstruct_ocmp9t
+!
+!===============================================================================
+!
 ! subroutine RECONSTRUCT_OCMP7P:
 ! -----------------------------
 !