diff --git a/sources/interpolations.F90 b/sources/interpolations.F90 index c04bad5..ab77d7a 100644 --- a/sources/interpolations.F90 +++ b/sources/interpolations.F90 @@ -138,15 +138,6 @@ module interpolations real(kind=8), dimension(2), parameter :: & ce2 = [ 5.0d-01, 5.0d-01 ] -! coefficients for the optimized compact schemes -! - logical :: oc_scheme = .false. - logical :: oc_stable = .false. - real(kind=8), dimension(2) :: oc_a = 0.0d+00 - real(kind=8), dimension(3) :: di5 - real(kind=8), dimension(5) :: di7, di9, ci5, ci7 - real(kind=8), dimension(7) :: ci9 - ! by default everything is private ! private @@ -231,15 +222,6 @@ module interpolations call get_parameter("central_weight" , cweight ) call get_parameter("cfl" , cfl ) - stmp = 'minimum' - call get_parameter("ocmp_scheme_mode", stmp) - select case(trim(stmp)) - case ("stable", "STABLE", "stability", "STABILITY") - oc_stable = .true. - case default - oc_stable = .false. - end select - ! calculate κ = (1 - ν) / ν ! kappa = min(kappa, (1.0d+00 - cfl) / cfl) @@ -365,60 +347,34 @@ module interpolations reconstruct_states => reconstruct_ocmp5 order = 5 nghosts = max(nghosts, 4) - oc_scheme = .true. - if (oc_stable) then - oc_a = [ 4.98486328125000044409d-01, 2.50756835937499977796d-01 ] - else - oc_a = [ 5.01779599466318559919d-01, 2.54377898581479855444d-01 ] - end if - di5 = [ oc_a(1), 1.0d+00, oc_a(2) ] - ci5 = [- 3.0d+00 * oc_a(1) - 3.0d+00 * oc_a(2) + 2.0d+00, & - 2.7d+01 * oc_a(1) + 1.7d+01 * oc_a(2) - 1.3d+01, & - 4.7d+01 * oc_a(1) - 4.3d+01 * oc_a(2) + 4.7d+01, & - - 1.3d+01 * oc_a(1) + 7.7d+01 * oc_a(2) + 2.7d+01, & - 2.0d+00 * oc_a(1) + 1.2d+01 * oc_a(2) - 3.0d+00 ] / 6.0d+01 - case ("ocmp7", "OCMP7") - name_rec = "7th order Optimized Compact Monotonicity Preserving" + case ("ocmp7t", "OCMP7t") + name_rec = "7th order Tridiagonal Optimized " // & + "Compact Monotonicity Preserving" interfaces => interfaces_dir - reconstruct_states => reconstruct_ocmp7 + reconstruct_states => reconstruct_ocmp7t order = 7 nghosts = max(nghosts, 6) - oc_scheme = .true. - if (oc_stable) then - oc_a = [ 6.64390055338541563046d-01, 3.34471638997395903647d-01 ] - else - oc_a = [ 6.71157762607103580699d-01, 3.40895359923885754583d-01 ] - end if - di7 = [ (3.0d+00 * oc_a(1) + 2.0d+00 * oc_a(2) - 2.0d+00) / 8.0d+00, & - oc_a(1), 1.0d+00, oc_a(2), & - ( oc_a(1) + 6.0d+00 * oc_a(2) - 2.0d+00) / 4.0d+01 ] - ci7 = [ 1.80d+01 * oc_a(1) + 1.80d+01 * oc_a(2) - 1.60d+01, & - 5.43d+02 * oc_a(1) + 2.68d+02 * oc_a(2) - 2.91d+02, & - 3.43d+02 * oc_a(1) - 3.32d+02 * oc_a(2) + 5.09d+02, & - -1.07d+02 * oc_a(1) + 5.68d+02 * oc_a(2) + 3.09d+02, & - 4.30d+01 * oc_a(1) + 3.18d+02 * oc_a(2) - 9.10d+01 ] / 6.0d+02 - case ("ocmp9", "OCMP9") - name_rec = "9th order Optimized Compact Monotonicity Preserving" + case ("ocmp9t", "OCMP9t") + name_rec = "9th order Tridiagonal Optimized " // & + "Compact Monotonicity Preserving" interfaces => interfaces_dir - reconstruct_states => reconstruct_ocmp9 + reconstruct_states => reconstruct_ocmp9t + order = 7 + nghosts = max(nghosts, 6) + case ("ocmp7", "OCMP7", "ocmp7p", "OCMP7P") + name_rec = "7th order Pentadiagonal Optimized " // & + "Compact Monotonicity Preserving" + interfaces => interfaces_dir + reconstruct_states => reconstruct_ocmp7p + order = 7 + nghosts = max(nghosts, 6) + case ("ocmp9", "OCMP9", "ocmp9p", "OCMP9P") + name_rec = "9th order Pentadiagonal Optimized " // & + "Compact Monotonicity Preserving" + interfaces => interfaces_dir + reconstruct_states => reconstruct_ocmp9p order = 9 nghosts = max(nghosts, 6) - oc_scheme = .true. - if (oc_stable) then - oc_a = [ 6.64322884877522779057d-01, 4.01406269073486310361d-01 ] - else - oc_a = [ 6.70794426721399439373d-01, 4.09731006175920509094d-01 ] - end if - di9 = [ (9.0d+00 * oc_a(1) + 5.0d+00 * oc_a(2) - 6.0d+00) / 2.0d+01, & - oc_a(1), 1.0d+00, oc_a(2), & - ( oc_a(1) + 5.0d+00 * oc_a(2) - 2.0d+00) / 2.0d+01 ] - ci9 = [-1.000d+01 * oc_a(1) -1.000d+01 * oc_a(2) +1.000d+01, & - 2.420d+02 * oc_a(1) +2.000d+02 * oc_a(2) -2.140d+02, & - 4.085d+03 * oc_a(1) +1.635d+03 * oc_a(2) -2.146d+03, & - 2.335d+03 * oc_a(1) -1.865d+03 * oc_a(2) +3.454d+03, & - -8.150d+02 * oc_a(1) +3.385d+03 * oc_a(2) +2.404d+03, & - 4.450d+02 * oc_a(1) +2.895d+03 * oc_a(2) -9.560d+02, & - 1.800d+01 * oc_a(1) +6.000d+01 * oc_a(2) -3.200d+01 ] / 4.2d+03 case ("gp", "GP") write(stmp, '(f16.1)') sgp write(name_rec, '("Gaussian Process (",i1,"-point, δ=",a,")")') ngp & @@ -709,13 +665,6 @@ module interpolations call print_section(verbose, "Interpolations") call print_parameter(verbose, "reconstruction" , name_rec ) - if (oc_scheme) then - if (oc_stable) then - call print_parameter(verbose, "scheme optimized for", "stability") - else - call print_parameter(verbose, "scheme optimized for", "minimum error") - end if - end if call print_parameter(verbose, "TVD limiter" , name_tlim) call print_parameter(verbose, "prolongation limiter", name_plim) if (mlp) then @@ -4050,7 +3999,7 @@ module interpolations do i = ng, n - ng + 1 r(i) = sum(ci5a(:) * fi(i-1:i+2)) - end do ! i = ng, n - ng + 1 + end do r( 1) = sum(ce2(:) * fi( 1: 2)) r( 2) = sum(ce3(:) * fi( 1: 3)) @@ -4165,7 +4114,7 @@ module interpolations do i = ng, n - ng + 1 r(i) = sum(ci7a(:) * fi(i-2:i+3)) - end do ! i = ng, n - ng + 1 + end do r( 1) = sum(ce2(:) * fi( 1: 2)) r( 2) = sum(ce3(:) * fi( 1: 3)) @@ -4290,7 +4239,7 @@ module interpolations do i = ng, n - ng + 1 r(i) = sum(ci9a(:) * fi(i-2:i+3)) - end do ! i = ng, n - ng + 1 + end do r( 1) = sum(ce2(:) * fi( 1: 2)) r( 2) = sum(ce3(:) * fi( 1: 3)) @@ -4326,11 +4275,16 @@ module interpolations !=============================================================================== ! ! subroutine RECONSTRUCT_OCMP5: -! ----------------------------- +! ---------------------------- ! ! Subroutine reconstructs the interface states using the 5th 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 6.4e-12. +! ! Arguments are described in subroutine reconstruct(). ! ! References: @@ -4362,6 +4316,16 @@ module interpolations real(kind=8), dimension(size(fc)) :: r real(kind=8), dimension(size(fc)) :: u + real(kind=8), parameter :: a1 = 5.0148583220385203994203557d-01 + real(kind=8), parameter :: a2 = 2.5452055302783787784904423d-01 + real(kind=8), dimension(3), parameter :: di5 = [ a1, 1.0d+00, a2 ] + real(kind=8), dimension(5), parameter :: & + ci5 = [- 3.0d+00 * a1 - 3.0d+00 * a2 + 2.0d+00, & + 2.7d+01 * a1 + 1.7d+01 * a2 - 1.3d+01, & + 4.7d+01 * a1 - 4.3d+01 * a2 + 4.7d+01, & + - 1.3d+01 * a1 + 7.7d+01 * a2 + 2.7d+01, & + 2.0d+00 * a1 + 1.2d+01 * a2 - 3.0d+00 ] / 6.0d+01 + !------------------------------------------------------------------------------- ! n = size(fc) @@ -4411,7 +4375,7 @@ module interpolations do i = ng, n - ng + 1 r(i) = sum(ci5(:) * fi(i-2:i+2)) - end do ! i = ng, n - ng + 1 + end do r( 1) = sum(ce2(:) * fi( 1: 2)) r( 2) = sum(ce3(:) * fi( 1: 3)) @@ -4431,8 +4395,8 @@ module interpolations fr(1:n-1) = u(n-1:1:-1) i = n - 1 - fl(1) = 0.5d+00 * (fc(1) + fc(2)) - fr(i) = 0.5d+00 * (fc(i) + fc(n)) + 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) @@ -4442,26 +4406,286 @@ module interpolations ! !=============================================================================== ! -! subroutine RECONSTRUCT_OCMP7: +! subroutine RECONSTRUCT_OCMP7T: ! ----------------------------- ! ! Subroutine reconstructs the interface states using the 7th 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.0e-11. +! ! Arguments are described in subroutine reconstruct(). ! -! References: -! -! [1] 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 -! !=============================================================================== ! - subroutine reconstruct_ocmp7(p, h, fc, fl, fr) + subroutine reconstruct_ocmp7t(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.0119544102147968726471782d-01 + real(kind=8), parameter :: a2 = 3.0649545768561641522180678d-01 + real(kind=8), dimension(3), parameter :: di7 = [ a1, 1.0d+00, a2 ] + real(kind=8), dimension(7), parameter :: & + ci7 = [ 4.00d+00 * a1 + 4.00d+00 * a2 - 3.00d+00, & + - 3.80d+01 * a1 - 3.10d+01 * a2 + 2.50d+01, & + 2.14d+02 * a1 + 1.09d+02 * a2 - 1.01d+02, & + 3.19d+02 * a1 - 2.41d+02 * a2 + 3.19d+02, & + - 1.01d+02 * a1 + 4.59d+02 * a2 + 2.14d+02, & + 2.50d+01 * a1 + 1.30d+02 * a2 - 3.80d+01, & + - 3.00d+00 * a1 - 1.00d+01 * a2 + 4.00d+00 ] / 4.2d+02 + +!------------------------------------------------------------------------------- +! + 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) = di7(1) + b(i) = di7(2) + c(i) = di7(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(ci7(:) * fc(i-3:i+3)) + end do + + 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)) + end do + do i = n - ng, n - 3 + r(i) = sum(ce7(:) * fc(i-3:i+3)) + end do + 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(ci7(:) * fi(i-3:i+3)) + end do + + 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)) + end do + do i = n - ng, n - 3 + r(i) = sum(ce7(:) * fi(i-3:i+3)) + end do + 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_ocmp7t +! +!=============================================================================== +! +! 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: +! ----------------------------- +! +! Subroutine reconstructs the interface states using the 7th order Optimized +! Compact Reconstruction Monotonicity Preserving (CRMP) method. +! +! The method uses an implicit penta-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.38e-11. +! +! Arguments are described in subroutine reconstruct(). +! +!=============================================================================== +! + subroutine reconstruct_ocmp7p(p, h, fc, fl, fr) use algebra, only : pentadiag @@ -4479,12 +4703,23 @@ module interpolations real(kind=8), dimension(size(fc)) :: r real(kind=8), dimension(size(fc)) :: u + real(kind=8), parameter :: a1 = 6.7100477315625665203995935d-01 + real(kind=8), parameter :: a2 = 3.4097048071278782662422942d-01 + real(kind=8), dimension(5), parameter :: & + di7 = [ (3.0d+00 * a1 + 2.0d+00 * a2 - 2.0d+00) / 8.0d+00, & + a1, 1.0d+00, a2, & + ( a1 + 6.0d+00 * a2 - 2.0d+00) / 4.0d+01 ] + real(kind=8), dimension(5), parameter :: & + ci7 = [ 1.80d+01 * a1 + 1.80d+01 * a2 - 1.60d+01, & + 5.43d+02 * a1 + 2.68d+02 * a2 - 2.91d+02, & + 3.43d+02 * a1 - 3.32d+02 * a2 + 5.09d+02, & + -1.07d+02 * a1 + 5.68d+02 * a2 + 3.09d+02, & + 4.30d+01 * a1 + 3.18d+02 * a2 - 9.10d+01 ] / 6.0d+02 + !------------------------------------------------------------------------------- ! n = size(fc) -! prepare the diagonals of the tridiagonal matrix -! do i = 1, ng e(i) = 0.0d+00 c(i) = 0.0d+00 @@ -4538,7 +4773,7 @@ module interpolations do i = ng, n - ng + 1 r(i) = sum(ci7(:) * fi(i-2:i+2)) - end do ! i = ng, n - ng + 1 + end do r( 1) = sum(ce2(:) * fi( 1: 2)) r( 2) = sum(ce3(:) * fi( 1: 3)) @@ -4560,37 +4795,33 @@ module interpolations fr(1:n-1) = u(n-1:1:-1) i = n - 1 - fl(1) = 0.5d+00 * (fc(1) + fc(2)) - fr(i) = 0.5d+00 * (fc(i) + fc(n)) + 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_ocmp7 + end subroutine reconstruct_ocmp7p ! !=============================================================================== ! -! subroutine RECONSTRUCT_OCMP9: +! subroutine RECONSTRUCT_OCMP9P: ! ----------------------------- ! ! Subroutine reconstructs the interface states using the 9th order Optimized ! Compact Reconstruction Monotonicity Preserving (CRMP) method. ! +! The method uses an implicit penta-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.54e-11. +! ! Arguments are described in subroutine reconstruct(). ! -! References: -! -! [1] 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 -! !=============================================================================== ! - subroutine reconstruct_ocmp9(p, h, fc, fl, fr) + subroutine reconstruct_ocmp9p(p, h, fc, fl, fr) use algebra, only : pentadiag @@ -4608,12 +4839,25 @@ module interpolations real(kind=8), dimension(size(fc)) :: r real(kind=8), dimension(size(fc)) :: u + real(kind=8), parameter :: a1 = 6.7028117788580347100541103d-01 + real(kind=8), parameter :: a2 = 4.1002518667894543684869155d-01 + real(kind=8), dimension(5), parameter :: & + di9 = [ (9.0d+00 * a1 + 5.0d+00 * a2 - 6.0d+00) / 2.0d+01, & + a1, 1.0d+00, a2, & + ( a1 + 5.0d+00 * a2 - 2.0d+00) / 2.0d+01 ] + real(kind=8), dimension(7), parameter :: & + ci9 = [-1.000d+01 * a1 -1.000d+01 * a2 +1.000d+01, & + 2.420d+02 * a1 +2.000d+02 * a2 -2.140d+02, & + 4.085d+03 * a1 +1.635d+03 * a2 -2.146d+03, & + 2.335d+03 * a1 -1.865d+03 * a2 +3.454d+03, & + -8.150d+02 * a1 +3.385d+03 * a2 +2.404d+03, & + 4.450d+02 * a1 +2.895d+03 * a2 -9.560d+02, & + 1.800d+01 * a1 +6.000d+01 * a2 -3.200d+01 ] / 4.2d+03 + !------------------------------------------------------------------------------- ! n = size(fc) -! prepare the diagonals of the tridiagonal matrix -! do i = 1, ng e(i) = 0.0d+00 c(i) = 0.0d+00 @@ -4669,7 +4913,7 @@ module interpolations do i = ng, n - ng + 1 r(i) = sum(ci9(:) * fi(i-3:i+3)) - end do ! i = ng, n - ng + 1 + end do r( 1) = sum(ce2(:) * fi( 1: 2)) r( 2) = sum(ce3(:) * fi( 1: 3)) @@ -4693,14 +4937,14 @@ module interpolations fr(1:n-1) = u(n-1:1:-1) i = n - 1 - fl(1) = 0.5d+00 * (fc(1) + fc(2)) - fr(i) = 0.5d+00 * (fc(i) + fc(n)) + 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_ocmp9 + end subroutine reconstruct_ocmp9p ! !=============================================================================== !