diff --git a/AMSS_NCKU_source/fmisc.f90 b/AMSS_NCKU_source/fmisc.f90
index 64f3b43..fee5c45 100644
--- a/AMSS_NCKU_source/fmisc.f90
+++ b/AMSS_NCKU_source/fmisc.f90
@@ -1241,10 +1241,10 @@ end subroutine d2dump
     y = y + dy
   end do
 
-  return
-  end subroutine polint
-
-  subroutine polint0(xa, ya, y, ordn)
+  return
+  end subroutine polint
+
+  subroutine polint0(xa, ya, y, ordn)
   ! Lagrange interpolation at x=0, O(n) direct formula
   implicit none
   integer, intent(in) :: ordn
@@ -1264,14 +1264,54 @@ end subroutine d2dump
     enddo
     y = y + wj * ya(j)
   enddo
-
-  return
-  end subroutine polint0
-!------------------------------------------------------------------------------
-!
-! interpolation in 2 dimensions, follow yx order
-!
-!------------------------------------------------------------------------------
+
+  return
+  end subroutine polint0
+!------------------------------------------------------------------------------
+!
+! interpolation in 2 dimensions, follow yx order
+!
+!------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
+! Compute Lagrange interpolation basis weights for one target point.
+!------------------------------------------------------------------------------
+!DIR$ ATTRIBUTES FORCEINLINE :: polint_lagrange_weights
+  subroutine polint_lagrange_weights(xa, x, w, ordn)
+  implicit none
+
+  integer, intent(in) :: ordn
+  real*8, dimension(1:ordn), intent(in) :: xa
+  real*8, intent(in) :: x
+  real*8, dimension(1:ordn), intent(out) :: w
+
+  integer :: i, j
+  real*8 :: num, den, dx
+
+  do i = 1, ordn
+    num = 1.d0
+    den = 1.d0
+    do j = 1, ordn
+      if (j /= i) then
+        dx = xa(i) - xa(j)
+        if (dx == 0.0d0) then
+          write(*,*) 'failure in polint for point',x
+          write(*,*) 'with input points: ',xa
+          stop
+        end if
+        num = num * (x - xa(j))
+        den = den * dx
+      end if
+    end do
+    w(i) = num / den
+  end do
+
+  return
+  end subroutine polint_lagrange_weights
+!------------------------------------------------------------------------------
+!
+! interpolation in 2 dimensions, follow yx order
+!
+!------------------------------------------------------------------------------
   subroutine polin2(x1a,x2a,ya,x1,x2,y,dy,ordn)
   implicit none
 
@@ -1319,11 +1359,11 @@ end subroutine d2dump
   real*8, intent(in) :: x1,x2,x3
   real*8, intent(out) :: y,dy
 
-#ifdef POLINT_LEGACY_ORDER
-  integer  :: i,j,m,n
-  real*8, dimension(ordn,ordn) :: yatmp
-  real*8, dimension(ordn) :: ymtmp
-  real*8, dimension(ordn) :: yntmp
+#ifdef POLINT_LEGACY_ORDER
+  integer  :: i,j,m,n
+  real*8, dimension(ordn,ordn) :: yatmp
+  real*8, dimension(ordn) :: ymtmp
+  real*8, dimension(ordn) :: yntmp
   real*8, dimension(ordn) :: yqtmp
 
   m=size(x1a)
@@ -1333,29 +1373,36 @@ end subroutine d2dump
     yqtmp=ya(i,j,:)
     call polint(x3a,yqtmp,x3,yatmp(i,j),dy,ordn)
    end do
-    yntmp=yatmp(i,:)
-    call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
-  end do
-  call polint(x1a,ymtmp,x1,y,dy,ordn)
-#else
-  integer  :: j, k
-  real*8, dimension(ordn,ordn) :: yatmp
-  real*8, dimension(ordn) :: ymtmp
-  real*8 :: dy_temp
-
-  do k=1,ordn
-    do j=1,ordn
-      call polint(x1a, ya(:,j,k), x1, yatmp(j,k), dy_temp, ordn)
-    end do
-  end do
-  do k=1,ordn
-    call polint(x2a, yatmp(:,k), x2, ymtmp(k), dy_temp, ordn)
-  end do
-  call polint(x3a, ymtmp, x3, y, dy, ordn)
-#endif
-
-  return
-  end subroutine polin3
+    yntmp=yatmp(i,:)
+    call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
+  end do
+  call polint(x1a,ymtmp,x1,y,dy,ordn)
+#else
+  integer  :: i, j, k
+  real*8, dimension(ordn) :: w1, w2
+  real*8, dimension(ordn) :: ymtmp
+  real*8 :: yx_sum, x_sum
+
+  call polint_lagrange_weights(x1a, x1, w1, ordn)
+  call polint_lagrange_weights(x2a, x2, w2, ordn)
+
+  do k = 1, ordn
+    yx_sum = 0.d0
+    do j = 1, ordn
+      x_sum = 0.d0
+      do i = 1, ordn
+        x_sum = x_sum + w1(i) * ya(i,j,k)
+      end do
+      yx_sum = yx_sum + w2(j) * x_sum
+    end do
+    ymtmp(k) = yx_sum
+  end do
+
+  call polint(x3a, ymtmp, x3, y, dy, ordn)
+#endif
+
+  return
+  end subroutine polin3
 !--------------------------------------------------------------------------------------
 ! calculate L2norm
   subroutine l2normhelper(ex, X, Y, Z,xmin,ymin,zmin,xmax,ymax,zmax,&