diff --git a/sources/interpolations.F90 b/sources/interpolations.F90 index 1cf144c..56295b9 100644 --- a/sources/interpolations.F90 +++ b/sources/interpolations.F90 @@ -5657,12 +5657,6 @@ module interpolations ! Science China Physics, Mechanics and Astronomy, ! Volume 54, Issue 3, pp. 511-522, ! http://dx.doi.org/10.1007/s11433-010-4220-x -! [3] Myeong-Hwan Ahn, Duck-Joo Lee, -! "Modified Monotonicity Preserving Constraints for High-Resolution -! Optimized Compact Scheme", -! Journal of Scientific Computing, -! 2020, vol. 83, p. 34 -! https://doi.org/10.1007/s10915-020-01221-0 ! !=============================================================================== ! @@ -5679,9 +5673,8 @@ module interpolations ! local variables ! - logical :: test - integer :: n, i, im2, im1, ip1, ip2 - real(kind=8) :: dq, ds, dc0, dc4, dm1, dp1, dml, dmr, bt + integer :: n, i, im1, ip1, ip2 + real(kind=8) :: dq, ds, dc0, dc4, dm1, dp1, dml, dmr real(kind=8) :: qlc, qmd, qmp, qmn, qmx, qul ! local vectors @@ -5689,10 +5682,12 @@ module interpolations real(kind=8), dimension(0:size(qc)+2) :: dm ! !------------------------------------------------------------------------------- +! +! get the input vector size ! n = size(qc) -! 1st order derivatives +! calculate derivatives ! dm(0 ) = 0.0d+00 dm(1 ) = 0.0d+00 @@ -5700,7 +5695,7 @@ module interpolations dm(n+1) = 0.0d+00 dm(n+2) = 0.0d+00 -! check the monotonicity condition and apply limiting if necessary +! check monotonicity condition for all elements and apply limiting if required ! do i = 1, n - 1 @@ -5717,7 +5712,6 @@ module interpolations if (ds > eps) then - im2 = i - 2 im1 = i - 1 ip2 = i + 2 @@ -5733,17 +5727,6 @@ module interpolations qul = qc(i) + dq qlc = qc(i) + 0.5d+00 * dq + dml - test = qc(i) > max(qc(im1), qc(ip1)) .and. & - min(qc(im1), qc(ip1)) > max(qc(im2), qc(ip2)) - test = test .or. qc(i) < min(qc(im1), qc(ip1)) .and. & - max(qc(im1), qc(ip1)) < min(qc(im2), qc(ip2)) - if (test) then - if (qc(im2) <= qc(ip1) .and. qc(ip2) <= qc(im1)) then - bt = dc0 / (qc(ip2) + qc(im2) - 2.0d+00 * qc(i)) - if (bt >= 3.0d-01) qlc = qc(im2) - end if - end if - qmx = max(min(qc(i), qc(ip1), qmd), min(qc(i), qul, qlc)) qmn = min(max(qc(i), qc(ip1), qmd), max(qc(i), qul, qlc))