8 implicit none ;
private
21 integer,
intent(in) :: n
22 real,
dimension(:),
intent(in) :: h
23 real,
dimension(:),
intent(in) :: u
24 real,
dimension(:,:),
intent(inout) :: ppoly_e
26 real,
dimension(:,:),
intent(inout) :: ppoly_s
28 real,
dimension(:,:),
intent(inout) :: ppoly_coef
30 real,
optional,
intent(in) :: h_neglect
42 call pqm_limiter( n, h, u, ppoly_e, ppoly_s, h_neglect )
57 c = 30.0 * u(k) - 12.0*u0_r - 18.0*u0_l + 1.5*h_c*(u1_r - 3.0*u1_l)
58 d = -60.0 * u(k) + h_c *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
59 e = 30.0 * u(k) + 2.5*h_c*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
78 subroutine pqm_limiter( N, h, u, ppoly_E, ppoly_S, h_neglect )
79 integer,
intent(in) :: N
80 real,
dimension(:),
intent(in) :: h
81 real,
dimension(:),
intent(in) :: u
82 real,
dimension(:,:),
intent(inout) :: ppoly_E
84 real,
dimension(:,:),
intent(inout) :: ppoly_S
86 real,
optional,
intent(in) :: h_neglect
91 integer :: inflexion_l
92 integer :: inflexion_r
97 real :: sigma_l, sigma_c, sigma_r
100 real :: alpha1, alpha2, alpha3
101 real :: rho, sqrt_rho
102 real :: gradient1, gradient2
106 hneglect =
hneglect_dflt ;
if (
present(h_neglect)) hneglect = h_neglect
138 sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + hneglect )
139 sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r + hneglect )
140 sigma_r = 2.0 * ( u_r - u_c ) / ( h_c + hneglect )
142 if ( (sigma_l * sigma_r) > 0.0 )
then
143 slope = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c )
150 if ( u1_l*slope <= 0.0 ) u1_l = slope
151 if ( u1_r*slope <= 0.0 ) u1_r = slope
154 if ( (u0_r - u_c) * (u_c - u0_l) <= 0.0)
then
167 if ( (inflexion_l == 0) .AND. (inflexion_r == 0) )
then
171 c = 30.0 * u(k) - 12.0*u0_r - 18.0*u0_l + 1.5*h_c*(u1_r - 3.0*u1_l)
172 d = -60.0 * u(k) + h_c *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
173 e = 30.0 * u(k) + 2.5*h_c*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
181 rho = alpha2 * alpha2 - 4.0 * alpha1 * alpha3
184 if (( alpha1 /= 0.0 ) .and. ( rho >= 0.0 ))
then
186 sqrt_rho = sqrt( rho )
188 x1 = 0.5 * ( - alpha2 - sqrt_rho ) / alpha1
189 x2 = 0.5 * ( - alpha2 + sqrt_rho ) / alpha1
192 if ( (x1 >= 0.0) .AND. (x1 <= 1.0) .AND. &
193 (x2 >= 0.0) .AND. (x2 <= 1.0) )
then
195 gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
196 gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b
199 if ( (gradient1 * slope < 0.0) .OR. &
200 (gradient2 * slope < 0.0) )
then
203 if ( abs(sigma_l) < abs(sigma_r) )
then
212 elseif ( (x1 >= 0.0) .AND. (x1 <= 1.0) )
then
214 gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
217 if ( gradient1 * slope < 0.0 )
then
220 if ( abs(sigma_l) < abs(sigma_r) )
then
228 elseif ( (x2 >= 0.0) .AND. (x2 <= 1.0) )
then
230 gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b
233 if ( gradient2 * slope < 0.0 )
then
236 if ( abs(sigma_l) < abs(sigma_r) )
then
249 if (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 ))
then
251 x1 = - alpha3 / alpha2
252 if ( (x1 >= 0.0) .AND. (x1 <= 1.0) )
then
254 gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
257 if ( gradient1 * slope < 0.0 )
then
260 if ( abs(sigma_l) < abs(sigma_r) )
then
274 if ( inflexion_l == 1 )
then
278 u1_l = ( 10.0 * u_c - 2.0 * u0_r - 8.0 * u0_l ) / (3.0*h_c + hneglect )
279 u1_r = ( -10.0 * u_c + 6.0 * u0_r + 4.0 * u0_l ) / ( h_c + hneglect )
285 if ( u1_l * slope < 0.0 )
then
288 u0_r = 5.0 * u_c - 4.0 * u0_l
289 u1_r = 20.0 * (u_c - u0_l) / ( h_c + hneglect )
291 elseif ( u1_r * slope < 0.0 )
then
294 u0_l = (5.0*u_c - 3.0*u0_r) / 2.0
295 u1_l = 10.0 * (-u_c + u0_r) / (3.0 * h_c + hneglect)
299 elseif ( inflexion_r == 1 )
then
303 u1_r = ( -10.0 * u_c + 8.0 * u0_r + 2.0 * u0_l ) / (3.0 * h_c + hneglect)
304 u1_l = ( 10.0 * u_c - 4.0 * u0_r - 6.0 * u0_l ) / (h_c + hneglect)
310 if ( u1_l * slope < 0.0 )
then
313 u0_r = ( 5.0 * u_c - 3.0 * u0_l ) / 2.0
314 u1_r = 10.0 * (u_c - u0_l) / (3.0 * h_c + hneglect)
316 elseif ( u1_r * slope < 0.0 )
then
319 u0_l = 5.0 * u_c - 4.0 * u0_r
320 u1_l = 20.0 * ( -u_c + u0_r ) / (h_c + hneglect)
359 integer,
intent(in) :: n
360 real,
dimension(:),
intent(in) :: h
361 real,
dimension(:),
intent(in) :: u
362 real,
dimension(:,:),
intent(inout) :: ppoly_e
364 real,
dimension(:,:),
intent(inout) :: ppoly_coef
370 real :: a, b, c, d, e
391 slope = 2.0 * ( u1 - u0 )
392 if ( abs(u1_r) > abs(slope) )
then
403 u0_l = 3.0 * u0 + 0.5 * u1_r - 2.0 * u0_r
406 exp1 = (u0_r - u0_l) * (u0 - 0.5*(u0_l+u0_r))
407 exp2 = (u0_r - u0_l) * (u0_r - u0_l) / 6.0
409 if ( exp1 > exp2 )
then
410 u0_l = 3.0 * u0 - 2.0 * u0_r
413 if ( exp1 < -exp2 )
then
414 u0_r = 3.0 * u0 - 2.0 * u0_l
421 b = 6.0 * u0 - 4.0 * u0_l - 2.0 * u0_r
422 c = 3.0 * ( u0_r + u0_l - 2.0 * u0 )
428 ppoly_coef(i0,4) = 0.0
429 ppoly_coef(i0,5) = 0.0
445 u1_l = (b + 2*c + 3*d + 4*e)
446 u1_l = u1_l * (h1/h0)
449 slope = 2.0 * ( u1 - u0 )
450 if ( abs(u1_l) > abs(slope) )
then
461 u0_r = 3.0 * u1 - 0.5 * u1_l - 2.0 * u0_l
464 exp1 = (u0_r - u0_l) * (u1 - 0.5*(u0_l+u0_r))
465 exp2 = (u0_r - u0_l) * (u0_r - u0_l) / 6.0
467 if ( exp1 > exp2 )
then
468 u0_l = 3.0 * u1 - 2.0 * u0_r
471 if ( exp1 < -exp2 )
then
472 u0_r = 3.0 * u1 - 2.0 * u0_l
479 b = 6.0 * u1 - 4.0 * u0_l - 2.0 * u0_r
480 c = 3.0 * ( u0_r + u0_l - 2.0 * u1 )
486 ppoly_coef(i1,4) = 0.0
487 ppoly_coef(i1,5) = 0.0
508 integer,
intent(in) :: n
509 real,
dimension(:),
intent(in) :: h
510 real,
dimension(:),
intent(in) :: u
511 real,
dimension(:,:),
intent(inout) :: ppoly_e
513 real,
dimension(:,:),
intent(inout) :: ppoly_s
515 real,
dimension(:,:),
intent(inout) :: ppoly_coef
517 real,
optional,
intent(in) :: h_neglect
522 integer :: inflexion_l
523 integer :: inflexion_r
526 real :: a, b, c, d, e
532 real :: alpha1, alpha2, alpha3
533 real :: rho, sqrt_rho
534 real :: gradient1, gradient2
538 hneglect =
hneglect_dflt ;
if (
present(h_neglect)) hneglect = h_neglect
551 slope = 2.0 * ( u1 - u0 ) / ( ( h0 + h1 ) + hneglect )
560 u1_r = b / (h1 + hneglect)
565 beta = 2.0 * ( u0_r - um ) / ( (h0 + hneglect)*u1_r) - 1.0
569 br = u0_r + beta*u0_r - um
570 ar = um + beta*um - br
576 u_plm = um - 0.5 * slope
582 if ( abs(um-u0_l) < abs(um-u_plm) )
then
583 u1_l = 2.0 * ( br - ar*beta)
584 u1_l = u1_l / (h0 + hneglect)
587 u1_l = slope / (h0 + hneglect)
595 c = 30.0 * um - 12.0*u0_r - 18.0*u0_l + 1.5*h0*(u1_r - 3.0*u1_l)
596 d = -60.0 * um + h0 *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
597 e = 30.0 * um + 2.5*h0*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
603 rho = alpha2 * alpha2 - 4.0 * alpha1 * alpha3
607 if (( alpha1 /= 0.0 ) .and. ( rho >= 0.0 ))
then
609 sqrt_rho = sqrt( rho )
611 x1 = 0.5 * ( - alpha2 - sqrt_rho ) / alpha1
612 if ( (x1 > 0.0) .and. (x1 < 1.0) )
then
613 gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
614 if ( gradient1 * slope < 0.0 )
then
619 x2 = 0.5 * ( - alpha2 + sqrt_rho ) / alpha1
620 if ( (x2 > 0.0) .and. (x2 < 1.0) )
then
621 gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b
622 if ( gradient2 * slope < 0.0 )
then
629 if (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 ))
then
631 x1 = - alpha3 / alpha2
632 if ( (x1 >= 0.0) .and. (x1 <= 1.0) )
then
633 gradient1 = 3.0 * d * (x1**2) + 2.0 * c * x1 + b
634 if ( gradient1 * slope < 0.0 )
then
641 if ( inflexion_l == 1 )
then
645 u1_l = ( 10.0 * um - 2.0 * u0_r - 8.0 * u0_l ) / (3.0*h0 + hneglect)
646 u1_r = ( -10.0 * um + 6.0 * u0_r + 4.0 * u0_l ) / (h0 + hneglect)
652 if ( u1_l * slope < 0.0 )
then
655 u0_r = 5.0 * um - 4.0 * u0_l
656 u1_r = 20.0 * (um - u0_l) / ( h0 + hneglect )
658 elseif ( u1_r * slope < 0.0 )
then
661 u0_l = (5.0*um - 3.0*u0_r) / 2.0
662 u1_l = 10.0 * (-um + u0_r) / (3.0 * h0 + hneglect )
676 c = 30.0 * um - 12.0*u0_r - 18.0*u0_l + 1.5*h0*(u1_r - 3.0*u1_l)
677 d = -60.0 * um + h0 *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
678 e = 30.0 * um + 2.5*h0*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
698 slope = 2.0 * ( u1 - u0 ) / ( h0 + h1 )
708 u0_l = a + b + c + d + e
709 u1_l = (b + 2*c + 3*d + 4*e) / h0
713 if (um-u0_l /= 0.)
then
714 beta = 0.5*h1*u1_l / (um-u0_l) - 1.0
718 br = beta*um + um - u0_l
722 if (1+beta /= 0.)
then
723 u0_r = (ar + 2*br + beta*br ) / ((1+beta)*(1+beta))
725 u0_r = um + 0.5 * slope
729 u_plm = um + 0.5 * slope
735 if ( abs(um-u0_r) < abs(um-u_plm) )
then
736 u1_r = 2.0 * ( br - ar*beta ) / ( (1+beta)*(1+beta)*(1+beta) )
748 c = 30.0 * um - 12.0*u0_r - 18.0*u0_l + 1.5*h1*(u1_r - 3.0*u1_l)
749 d = -60.0 * um + h1*(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
750 e = 30.0 * um + 2.5*h1*(u1_r - u1_l) - 15.0*(u0_l + u0_r)
756 rho = alpha2 * alpha2 - 4.0 * alpha1 * alpha3
760 if (( alpha1 /= 0.0 ) .and. ( rho >= 0.0 ))
then
762 sqrt_rho = sqrt( rho )
764 x1 = 0.5 * ( - alpha2 - sqrt_rho ) / alpha1
765 if ( (x1 > 0.0) .and. (x1 < 1.0) )
then
766 gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
767 if ( gradient1 * slope < 0.0 )
then
772 x2 = 0.5 * ( - alpha2 + sqrt_rho ) / alpha1
773 if ( (x2 > 0.0) .and. (x2 < 1.0) )
then
774 gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b
775 if ( gradient2 * slope < 0.0 )
then
782 if (( alpha1 == 0.0 ) .and. ( alpha2 /= 0.0 ))
then
784 x1 = - alpha3 / alpha2
785 if ( (x1 >= 0.0) .and. (x1 <= 1.0) )
then
786 gradient1 = 3.0 * d * (x1**2) + 2.0 * c * x1 + b
787 if ( gradient1 * slope < 0.0 )
then
794 if ( inflexion_r == 1 )
then
798 u1_r = ( -10.0 * um + 8.0 * u0_r + 2.0 * u0_l ) / (3.0 * h1)
799 u1_l = ( 10.0 * um - 4.0 * u0_r - 6.0 * u0_l ) / h1
805 if ( u1_l * slope < 0.0 )
then
808 u0_r = ( 5.0 * um - 3.0 * u0_l ) / 2.0
809 u1_r = 10.0 * (um - u0_l) / (3.0 * h1)
811 elseif ( u1_r * slope < 0.0 )
then
814 u0_l = 5.0 * um - 4.0 * u0_r
815 u1_l = 20.0 * ( -um + u0_r ) / h1
829 c = 30.0 * um - 12.0*u0_r - 18.0*u0_l + 1.5*h1*(u1_r - 3.0*u1_l)
830 d = -60.0 * um + h1 *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l
831 e = 30.0 * um + 2.5*h1*(u1_r - u1_l) - 15.0*(u0_l + u0_r)