8 implicit none ;
private
29 integer,
intent(in) :: n
30 real,
dimension(:),
intent(in) :: h
31 real,
dimension(:),
intent(in) :: u
32 real,
dimension(:,:),
intent(inout) :: ppoly_e
33 real,
dimension(:,:),
intent(inout) :: ppoly_s
34 real,
dimension(:,:),
intent(inout) :: ppoly_coef
35 real,
optional,
intent(in) :: h_neglect
43 call p3m_limiter( n, h, u, ppoly_e, ppoly_s, ppoly_coef, h_neglect )
60 subroutine p3m_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect )
61 integer,
intent(in) :: N
62 real,
dimension(:),
intent(in) :: h
63 real,
dimension(:),
intent(in) :: u
64 real,
dimension(:,:),
intent(inout) :: ppoly_E
65 real,
dimension(:,:),
intent(inout) :: ppoly_S
66 real,
dimension(:,:),
intent(inout) :: ppoly_coef
67 real,
optional,
intent(in) :: h_neglect
76 real :: sigma_l, sigma_c, sigma_r
81 hneglect =
hneglect_dflt ;
if (
present(h_neglect)) hneglect = h_neglect
126 sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + hneglect )
127 sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r + hneglect )
128 sigma_r = 2.0 * ( u_r - u_c ) / ( h_c + hneglect )
130 if ( (sigma_l * sigma_r) > 0.0 )
then
131 slope = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c )
137 if ( abs(u1_l*h_c) < epsilon(u_c)*abs(u_c) ) u1_l = 0.0
138 if ( abs(u1_r*h_c) < epsilon(u_c)*abs(u_c) ) u1_r = 0.0
142 if ( abs(u1_l) > abs(sigma_l) )
then
146 if ( abs(u1_r) > abs(sigma_r) )
then
159 if ( .not.monotonic )
then
160 call monotonize_cubic( h_c, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r )
189 h_neglect, h_neglect_edge )
190 integer,
intent(in) :: n
191 real,
dimension(:),
intent(in) :: h
192 real,
dimension(:),
intent(in) :: u
193 real,
dimension(:,:),
intent(inout) :: ppoly_e
194 real,
dimension(:,:),
intent(inout) :: ppoly_s
195 real,
dimension(:,:),
intent(inout) :: ppoly_coef
196 real,
optional,
intent(in) :: h_neglect
198 real,
optional,
intent(in) :: h_neglect_edge
209 real :: hneglect, hneglect_edge
211 hneglect =
hneglect_dflt ;
if (
present(h_neglect)) hneglect = h_neglect
213 if (
present(h_neglect_edge)) hneglect_edge = h_neglect_edge
218 h0 = h(i0) + hneglect_edge
219 h1 = h(i1) + hneglect_edge
229 slope = 2.0 * ( u1 - u0 ) / ( h0 + hneglect )
230 if ( abs(u1_r) > abs(slope) )
then
241 u0_l = 3.0 * u0 + 0.5 * h0*u1_r - 2.0 * u0_r
242 u1_l = ( - 6.0 * u0 - 2.0 * h0*u1_r + 6.0 * u0_r) / ( h0 + hneglect )
247 if ( (u0_r-u0_l) * slope < 0.0 )
then
263 if ( .not.monotonic )
then
276 h0 = h(i0) + hneglect_edge
277 h1 = h(i1) + hneglect_edge
286 u1_l = (b + 2*c + 3*d) / ( h0 + hneglect )
289 slope = 2.0 * ( u1 - u0 ) / ( h1 + hneglect )
290 if ( abs(u1_l) > abs(slope) )
then
301 u0_r = 3.0 * u1 - 0.5 * h1*u1_l - 2.0 * u0_l
302 u1_r = ( 6.0 * u1 - 2.0 * h1*u1_l - 6.0 * u0_l) / ( h1 + hneglect )
307 if ( (u0_r-u0_l) * slope < 0.0 )
then
322 if ( .not.monotonic )
then
342 real,
dimension(:),
intent(in) :: h
343 integer,
intent(in) :: k
344 real,
dimension(:,:),
intent(in) :: ppoly_E
345 real,
dimension(:,:),
intent(in) :: ppoly_S
346 real,
dimension(:,:),
intent(inout) :: ppoly_coef
352 real :: a0, a1, a2, a3
359 u1_l = ppoly_s(k,1) * h_c
360 u1_r = ppoly_s(k,2) * h_c
364 a2 = 3.0 * ( u0_r - u0_l ) - u1_r - 2.0 * u1_l
365 a3 = u1_r + u1_l + 2.0 * ( u0_l - u0_r )
384 real,
dimension(:,:),
intent(in) :: ppoly_coef
385 integer,
intent(in) :: k
390 b = 2.0 * ppoly_coef(k,3)
391 c = 3.0 * ppoly_coef(k,4)
394 if (b*b - 4.0*a*c <= 0.0)
then
396 elseif (a * (a + (b + c)) < 0.0)
then
398 elseif (b * (b + 2.0*c) < 0.0)
then
433 subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r )
434 real,
intent(in) :: h
435 real,
intent(in) :: u0_l
436 real,
intent(in) :: u0_r
437 real,
intent(in) :: sigma_l
438 real,
intent(in) :: sigma_r
439 real,
intent(in) :: slope
440 real,
intent(inout) :: u1_l
441 real,
intent(inout) :: u1_r
444 logical :: inflexion_l
445 logical :: inflexion_r
453 inflexion_l = .false.
454 inflexion_r = .false.
458 if ( u1_l*slope <= 0.0 )
then
462 if ( u1_r*slope <= 0.0 )
then
469 a2 = 3.0 * ( u0_r - u0_l ) - h*(u1_r + 2.0*u1_l)
470 a3 = h*(u1_r + u1_l) + 2.0*(u0_l - u0_r)
475 if ( a3 /= 0.0 )
then
477 xi_ip = - a2 / (3.0 * a3)
480 if ( (xi_ip >= 0.0) .AND. (xi_ip <= 1.0) )
then
491 slope_ip = a1 + 2.0*a2*xi_ip + 3.0*a3*xi_ip*xi_ip
494 if ( slope_ip*slope < 0.0 )
then
495 if ( abs(sigma_l) < abs(sigma_r) )
then
508 if ( inflexion_l )
then
510 u1_l_tmp = 1.5*(u0_r-u0_l)/h - 0.5*u1_r
511 u1_r_tmp = 3.0*(u0_r-u0_l)/h - 2.0*u1_l
513 if ( (u1_l_tmp*slope < 0.0) .AND. (u1_r_tmp*slope < 0.0) )
then
516 u1_r = 3.0 * (u0_r - u0_l) / h
518 elseif (u1_l_tmp*slope < 0.0)
then
521 u1_l = 1.5*(u0_r - u0_l)/h - 0.5*u1_r
523 elseif (u1_r_tmp*slope < 0.0)
then
526 u1_r = 3.0*(u0_r - u0_l)/h - 2.0*u1_l
538 if ( inflexion_r )
then
540 u1_l_tmp = 3.0*(u0_r-u0_l)/h - 2.0*u1_r
541 u1_r_tmp = 1.5*(u0_r-u0_l)/h - 0.5*u1_l
543 if ( (u1_l_tmp*slope < 0.0) .AND. (u1_r_tmp*slope < 0.0) )
then
545 u1_l = 3.0 * (u0_r - u0_l) / h
548 elseif (u1_l_tmp*slope < 0.0)
then
551 u1_l = 3.0*(u0_r - u0_l)/h - 2.0*u1_r
553 elseif (u1_r_tmp*slope < 0.0)
then
556 u1_r = 1.5*(u0_r - u0_l)/h - 0.5*u1_l
568 if ( abs(u1_l*h) < epsilon(u0_l) * (abs(u0_l) + abs(u0_r)) ) u1_l = 0.0
569 if ( abs(u1_r*h) < epsilon(u0_l) * (abs(u0_l) + abs(u0_r)) ) u1_r = 0.0