diff --git a/sources/interpolations.F90 b/sources/interpolations.F90 index a611165..4a17082 100644 --- a/sources/interpolations.F90 +++ b/sources/interpolations.F90 @@ -42,7 +42,8 @@ module interpolations real(kind=8), dimension(:,:,:) , intent(in) :: q real(kind=8), dimension(:,:,:,:,:), intent(out) :: qi end subroutine - subroutine reconstruct_iface(h, fc, fl, fr) + subroutine reconstruct_iface(positive, h, fc, fl, fr) + logical , intent(in) :: positive real(kind=8) , intent(in) :: h real(kind=8), dimension(:), intent(in) :: fc real(kind=8), dimension(:), intent(out) :: fl, fr @@ -76,7 +77,6 @@ module interpolations ! monotonicity preserving reconstruction coefficients ! real(kind=8), save :: kappa = 1.0d+00 - real(kind=8), save :: kbeta = 1.0d+00 ! number of ghost zones (required for compact schemes) ! @@ -218,7 +218,6 @@ module interpolations call get_parameter("eps" , eps ) call get_parameter("limo3_rad" , rad ) 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 ) @@ -978,20 +977,20 @@ module interpolations ! ! Arguments: ! -! positive - the variable positivity flag; -! h - the spatial step; -! q - the variable array; -! qi - the array of reconstructed interfaces (2 in each direction); +! p - the variable positivity flag; +! h - the spatial step; +! q - the variable array; +! qi - the array of reconstructed interfaces (2 in each direction); ! !=============================================================================== ! - subroutine interfaces_dir(positive, h, q, qi) + subroutine interfaces_dir(p, h, q, qi) use coordinates, only : nb, ne, nbl, neu implicit none - logical , intent(in) :: positive + logical , intent(in) :: p real(kind=8), dimension(:) , intent(in) :: h real(kind=8), dimension(:,:,:) , intent(in) :: q real(kind=8), dimension(:,:,:,:,:), intent(out) :: qi @@ -1025,23 +1024,23 @@ module interpolations do k = nbl, neu #endif /* NDIMS == 3 */ do j = nbl, neu - call reconstruct(h(1), q(:,j,k), qi(:,j,k,1,1), qi(:,j,k,2,1)) + call reconstruct(p, h(1), q(:,j,k), qi(:,j,k,1,1), qi(:,j,k,2,1)) end do ! j = nbl, neu do i = nbl, neu - call reconstruct(h(2), q(i,:,k), qi(i,:,k,1,2), qi(i,:,k,2,2)) + call reconstruct(p, h(2), q(i,:,k), qi(i,:,k,1,2), qi(i,:,k,2,2)) end do ! i = nbl, neu #if NDIMS == 3 end do ! k = nbl, neu do j = nbl, neu do i = nbl, neu - call reconstruct(h(3), q(i,j,:), qi(i,j,:,1,3), qi(i,j,:,2,3)) + call reconstruct(p, h(3), q(i,j,:), qi(i,j,:,1,3), qi(i,j,:,2,3)) end do ! i = nbl, neu end do ! j = nbl, neu #endif /* NDIMS == 3 */ ! make sure the interface states are positive for positive variables ! - if (positive) then + if (p) then #if NDIMS == 3 do k = nbl, neu @@ -1442,6 +1441,7 @@ module interpolations ! ! Arguments: ! +! p - the flag indicating if the reconstructed variable is positivite; ! h - the spatial step; this is required for some reconstruction methods; ! f - the input vector of cell averaged values; ! fl - the left side state reconstructed for location (i+1/2); @@ -1449,10 +1449,11 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct(h, f, fl, fr) + subroutine reconstruct(p, h, f, fl, fr) implicit none + logical , intent(in) :: p real(kind=8) , intent(in) :: h real(kind=8), dimension(:), intent(in) :: f real(kind=8), dimension(:), intent(out) :: fl, fr @@ -1461,7 +1462,7 @@ module interpolations ! ! reconstruct the states using the selected subroutine ! - call reconstruct_states(h, f(:), fl(:), fr(:)) + call reconstruct_states(p, h, f(:), fl(:), fr(:)) ! correct the reconstruction near extrema by clipping them in order to improve ! the stability of scheme @@ -1485,10 +1486,11 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_tvd(h, fc, fl, fr) + subroutine reconstruct_tvd(p, h, fc, fl, fr) 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 @@ -1553,10 +1555,11 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_weno3(h, fc, fl, fr) + subroutine reconstruct_weno3(p, h, fc, fl, fr) 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 @@ -1675,10 +1678,11 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_limo3(h, fc, fl, fr) + subroutine reconstruct_limo3(p, h, fc, fl, fr) 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 @@ -1818,10 +1822,11 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_ppm(h, fc, fl, fr) + subroutine reconstruct_ppm(p, h, fc, fl, fr) 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 @@ -2002,10 +2007,11 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_weno5z(h, fc, fl, fr) + subroutine reconstruct_weno5z(p, h, fc, fl, fr) 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 @@ -2159,10 +2165,11 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_weno5yc(h, fc, fl, fr) + subroutine reconstruct_weno5yc(p, h, fc, fl, fr) 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 @@ -2317,10 +2324,11 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_weno5ns(h, fc, fl, fr) + subroutine reconstruct_weno5ns(p, h, fc, fl, fr) 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 @@ -2503,12 +2511,13 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_crweno5z(h, fc, fl, fr) + subroutine reconstruct_crweno5z(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 @@ -2878,12 +2887,13 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_crweno5yc(h, fc, fl, fr) + subroutine reconstruct_crweno5yc(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 @@ -3256,12 +3266,13 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_crweno5ns(h, fc, fl, fr) + subroutine reconstruct_crweno5ns(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 @@ -3647,10 +3658,11 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_mp5(h, fc, fl, fr) + subroutine reconstruct_mp5(p, h, fc, fl, fr) 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 @@ -3676,7 +3688,7 @@ module interpolations u(n-1) = sum(ce3(:) * fc(n-2:n)) u(n ) = fc(n ) - call mp_limiting(fc(:), u(:)) + call mp_limiting(p, n, fc(:), u(:)) fl(1:n) = u(1:n) @@ -3694,7 +3706,7 @@ module interpolations u(n-1) = sum(ce3(:) * fi(n-2:n)) u(n ) = fi(n ) - call mp_limiting(fi(:), u(:)) + call mp_limiting(p, n, fi(:), u(:)) fr(1:n-1) = u(n-1:1:-1) @@ -3724,10 +3736,11 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_mp7(h, fc, fl, fr) + subroutine reconstruct_mp7(p, h, fc, fl, fr) 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 @@ -3755,7 +3768,7 @@ module interpolations u(n-1) = sum(ce3(:) * fc(n-2:n)) u(n ) = fc(n ) - call mp_limiting(fc(:), u(:)) + call mp_limiting(p, n, fc(:), u(:)) fl(1:n) = u(1:n) @@ -3775,7 +3788,7 @@ module interpolations u(n-1) = sum(ce3(:) * fi(n-2:n)) u(n ) = fi(n ) - call mp_limiting(fi(:), u(:)) + call mp_limiting(p, n, fi(:), u(:)) fr(1:n-1) = u(n-1:1:-1) @@ -3805,10 +3818,11 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_mp9(h, fc, fl, fr) + subroutine reconstruct_mp9(p, h, fc, fl, fr) 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 @@ -3838,7 +3852,7 @@ module interpolations u(n-1) = sum(ce3(:) * fc(n-2:n)) u(n ) = fc(n ) - call mp_limiting(fc(:), u(:)) + call mp_limiting(p, n, fc(:), u(:)) fl(1:n) = u(1:n) @@ -3860,7 +3874,7 @@ module interpolations u(n-1) = sum(ce3(:) * fi(n-2:n)) u(n ) = fi(n ) - call mp_limiting(fi(:), u(:)) + call mp_limiting(p, n, fi(:), u(:)) fr(1:n-1) = u(n-1:1:-1) @@ -3902,12 +3916,13 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_crmp5(h, fc, fl, fr) + subroutine reconstruct_crmp5(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 @@ -3958,7 +3973,7 @@ module interpolations call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n)) - call mp_limiting(fc(:), u(:)) + call mp_limiting(p, n, fc(:), u(:)) fl(1:n) = u(1:n) @@ -3983,7 +3998,7 @@ module interpolations call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n)) - call mp_limiting(fi(:), u(:)) + call mp_limiting(p, n, fi(:), u(:)) fr(1:n-1) = u(n-1:1:-1) @@ -4014,12 +4029,13 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_crmp7(h, fc, fl, fr) + subroutine reconstruct_crmp7(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 @@ -4072,7 +4088,7 @@ module interpolations call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n)) - call mp_limiting(fc(:), u(:)) + call mp_limiting(p, n, fc(:), u(:)) fl(1:n) = u(1:n) @@ -4099,7 +4115,7 @@ module interpolations call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n)) - call mp_limiting(fi(:), u(:)) + call mp_limiting(p, n, fi(:), u(:)) fr(1:n-1) = u(n-1:1:-1) @@ -4130,12 +4146,13 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_crmp9(h, fc, fl, fr) + subroutine reconstruct_crmp9(p, h, fc, fl, fr) use algebra, only : pentadiag 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 @@ -4196,7 +4213,7 @@ module interpolations call pentadiag(n, e(:), c(:), d(:), a(:), b(:), r(:), u(:)) - call mp_limiting(fc(:), u(:)) + call mp_limiting(p, n, fc(:), u(:)) fl(1:n) = u(1:n) @@ -4225,7 +4242,7 @@ module interpolations call pentadiag(n, e(:), c(:), d(:), a(:), b(:), r(:), u(:)) - call mp_limiting(fi(:), u(:)) + call mp_limiting(p, n, fi(:), u(:)) fr(1:n-1) = u(n-1:1:-1) @@ -4260,12 +4277,13 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_ocmp5(h, fc, fl, fr) + subroutine reconstruct_ocmp5(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 @@ -4327,7 +4345,7 @@ module interpolations call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n)) - call mp_limiting(fc(:), u(:)) + call mp_limiting(p, n, fc(:), u(:)) fl(1:n) = u(1:n) @@ -4352,7 +4370,7 @@ module interpolations call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n)) - call mp_limiting(fi(:), u(:)) + call mp_limiting(p, n, fi(:), u(:)) fr(1:n-1) = u(n-1:1:-1) @@ -4387,12 +4405,13 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_ocmp7(h, fc, fl, fr) + subroutine reconstruct_ocmp7(p, h, fc, fl, fr) use algebra, only : pentadiag 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 @@ -4466,7 +4485,7 @@ module interpolations call pentadiag(n, e(:), c(:), d(:), a(:), b(:), r(:), u(:)) - call mp_limiting(fc(:), u(:)) + call mp_limiting(p, n, fc(:), u(:)) fl(1:n) = u(1:n) @@ -4493,7 +4512,7 @@ module interpolations call pentadiag(n, e(:), c(:), d(:), a(:), b(:), r(:), u(:)) - call mp_limiting(fi(:), u(:)) + call mp_limiting(p, n, fi(:), u(:)) fr(1:n-1) = u(n-1:1:-1) @@ -4528,12 +4547,13 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_ocmp9(h, fc, fl, fr) + subroutine reconstruct_ocmp9(p, h, fc, fl, fr) use algebra, only : pentadiag 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 @@ -4611,7 +4631,7 @@ module interpolations call pentadiag(n, e(:), c(:), d(:), a(:), b(:), r(:), u(:)) - call mp_limiting(fc(:), u(:)) + call mp_limiting(p, n, fc(:), u(:)) fl(1:n) = u(1:n) @@ -4640,7 +4660,7 @@ module interpolations call pentadiag(n, e(:), c(:), d(:), a(:), b(:), r(:), u(:)) - call mp_limiting(fi(:), u(:)) + call mp_limiting(p, n, fi(:), u(:)) fr(1:n-1) = u(n-1:1:-1) @@ -4773,10 +4793,11 @@ module interpolations ! !=============================================================================== ! - subroutine reconstruct_gp(h, fc, fl, fr) + subroutine reconstruct_gp(p, h, fc, fl, fr) 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 @@ -4814,7 +4835,7 @@ module interpolations ! apply the monotonicity preserving limiting ! - call mp_limiting(fc(:), u(:)) + call mp_limiting(p, n, fc(:), u(:)) ! copy the interpolation to the respective vector ! @@ -4848,7 +4869,7 @@ module interpolations ! apply the monotonicity preserving limiting ! - call mp_limiting(fi(:), u(:)) + call mp_limiting(p, n, fi(:), u(:)) ! copy the interpolation to the respective vector ! @@ -5325,93 +5346,118 @@ module interpolations ! subroutine MP_LIMITING: ! ---------------------- ! -! Subroutine applies the monotonicity preserving (MP) limiter to a vector of -! high order reconstructed interface values. +! Subroutine applies the monotonicity preserving (MP) limitation to +! the interface states reconstructed using a high order upwind +! interpolation. ! ! Arguments: ! -! qc - the vector of cell-centered values; -! qi - the vector of interface values obtained from the high order -! interpolation as input and its monotonicity limited values as output; +! p - the logical flag indicating if the variable is positive; +! n - the length of vector u and q; +! u - the vector of centered cell-integrated values; +! q - the vector of interface values obtained from the high order +! interpolation to which the limitation is applied; ! ! 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 +! Journal on Computational Physics, 1997, 136, 83-99, +! https://doi.org/10.1006/jcph.1997.5745 +! [2] Myeong-Hwan Ahn1, Duck-Joo Lee, +! "Modified Monotonicity Preserving Constraints for High-Resolution +! Optimized Compact Scheme", +! Journal of Scientific Computing, 2020, 83:34, +! https://doi.org/10.1007/s10915-020-01221-0 ! !=============================================================================== ! - subroutine mp_limiting(qc, qi) + subroutine mp_limiting(p, n, u, q) implicit none - real(kind=8), dimension(:), intent(in) :: qc - real(kind=8), dimension(:), intent(inout) :: qi + logical , intent(in) :: p + integer , intent(in) :: n + real(kind=8), dimension(n), intent(in) :: u + real(kind=8), dimension(n), intent(inout) :: q - integer :: n, i, im1, ip1, ip2 - real(kind=8) :: dq, ds, dc0, dc4, dm1, dp1, dml, dmr + integer :: i, im2, im1, ip1, ip2 + real(kind=8) :: du, dm, dc, dp, dc4, dmm, dmp, di, de, bt real(kind=8) :: qlc, qmd, qmp, qmn, qmx, qul - real(kind=8), dimension(0:size(qc)+2) :: dm + real(kind=8), dimension(n) :: d + real(kind=8), dimension(2) :: umn, umx !------------------------------------------------------------------------------- ! - n = size(qc) + d(1 ) = 0.0d+00 + d(2:n) = u(2:n) - u(1:n-1) - dm(0 ) = 0.0d+00 - dm(1 ) = 0.0d+00 - dm(2:n) = qc(2:n) - qc(1:n-1) - dm(n+1) = 0.0d+00 - dm(n+2) = 0.0d+00 + im2 = 1 + im1 = 1 + ip1 = 2 + ip2 = 3 - do i = 1, n - 1 + do i = 1, n - ip1 = i + 1 + du = kappa * d(i) - if (dm(i) * dm(ip1) >= 0.0d+00) then - dq = kappa * dm(i) - else - dq = kbeta * dm(i) - end if + qmp = u(i) + minmod(d(ip1), du) - qmp = qc(i) + minmod(dm(ip1), dq) - ds = (qi(i) - qc(i)) * (qi(i) - qmp) + if ((q(i) - u(i)) * (q(i) - qmp) > eps) then - if (ds > eps) then + dm = d(i ) - d(im1) + dc = d(ip1) - d(i ) + dp = d(ip2) - d(ip1) + dc4 = 4.0d+00 * dc - im1 = i - 1 - ip2 = i + 2 + dmp = minmod4(dc4 - dp, 4.0d+00 * dp - dc, dc, dp) + dmm = minmod4(dc4 - dm, 4.0d+00 * dm - dc, dc, dm) - dm1 = dm(i ) - dm(im1) - dc0 = dm(ip1) - dm(i ) - dp1 = dm(ip2) - dm(ip1) - dc4 = 4.0d+00 * dc0 + qul = u(i) + du + qmd = u(i) + 5.0d-01 * (d(ip1) - dmp ) + qlc = u(i) + 5.0d-01 * (d(i ) + 8.0d+00 * dmm / 3.0d+00) - dml = 0.5d+00 * minmod4(dc4 - dm1, 4.0d+00 * dm1 - dc0, dc0, dm1) - dmr = 0.5d+00 * minmod4(dc4 - dp1, 4.0d+00 * dp1 - dc0, dc0, dp1) + umn(1) = min(u(im1), u(ip1)) + umn(2) = min(u(im2), u(ip2)) + umx(1) = max(u(im1), u(ip1)) + umx(2) = max(u(im2), u(ip2)) - qmd = qc(i) + 0.5d+00 * dm(ip1) - dmr - qul = qc(i) + dq - qlc = qc(i) + 0.5d+00 * dq + dml + if ((u(i) > umx(1) .and. umx(1) > umx(2) .and. & + u(im2) <= u(ip1) .and. u(im1) >= u(ip2)) .or. & + (u(i) < umn(1) .and. umn(1) < umn(2) .and. & + u(im2) >= u(ip1) .and. u(im1) <= u(ip2))) then - 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)) + bt = 2.0d+00 * u(i) + di = (u(im1) + u(ip1)) - bt + de = (u(im2) + u(ip2)) - bt + bt = di / de + if (bt > 3.0d-01) qlc = u(im2) - qi(i) = median(qi(i), qmn, qmx) + end if + + qmx = max(min(u(i), u(ip1), qmd), min(u(i), qul, qlc)) + qmn = min(max(u(i), u(ip1), qmd), max(u(i), qul, qlc)) + + if (p .and. qmn <= 0.0d+00) then + if (d(i) <= 0.0d+00 .and. d(ip1) >= 0.0d+00) then + qmn = u(i) * u(ip1) / (u(i) + u(ip1)) + else + qmn = max(eps, min(u(i), u(ip1))) + end if + end if + + q(i) = median(q(i), qmn, qmx) end if - end do ! i = 1, n - 1 + im2 = im1 + im1 = i + ip1 = ip2 + ip2 = min(ip2 + 1, n) + + end do !------------------------------------------------------------------------------- !