14 implicit none ;
private
29 integer,
intent(in) :: n
30 real,
dimension(N),
intent(in) :: h
31 real,
dimension(N),
intent(in) :: u
32 real,
dimension(N,2),
intent(inout) :: ppoly_e
34 real,
dimension(N,3),
intent(inout) :: ppoly_coef
36 real,
optional,
intent(in) :: h_neglect
40 real :: edge_l, edge_r
52 ppoly_coef(k,1) = edge_l
53 ppoly_coef(k,2) = 4.0 * ( u(k) - edge_l ) + 2.0 * ( u(k) - edge_r )
54 ppoly_coef(k,3) = 3.0 * ( ( edge_r - u(k) ) + ( edge_l - u(k) ) )
64 integer,
intent(in) :: N
65 real,
dimension(:),
intent(in) :: h
66 real,
dimension(:),
intent(in) :: u
67 real,
dimension(:,:),
intent(inout) :: ppoly_E
69 real,
optional,
intent(in) :: h_neglect
74 real :: edge_l, edge_r
95 if ( (u_r - u_c)*(u_c - u_l) <= 0.0)
then
100 expr1 = 3.0 * (edge_r - edge_l) * ( (u_c - edge_l) + (u_c - edge_r))
101 expr2 = (edge_r - edge_l) * (edge_r - edge_l)
102 if ( expr1 > expr2 )
then
104 edge_l = u_c + 2.0 * ( u_c - edge_r )
105 edge_l = max( min( edge_l, max(u_l, u_c) ), min(u_l, u_c) )
106 elseif ( expr1 < -expr2 )
then
108 edge_r = u_c + 2.0 * ( u_c - edge_l )
109 edge_r = max( min( edge_r, max(u_r, u_c) ), min(u_r, u_c) )
114 if ( abs( edge_r - edge_l )<max(1.e-60,epsilon(u_c)*abs(u_c)) )
then
119 ppoly_e(k,1) = edge_l
120 ppoly_e(k,2) = edge_r
159 integer,
intent(in) :: n
160 real,
dimension(:),
intent(in) :: h
161 real,
dimension(:),
intent(in) :: u
162 real,
dimension(:,:),
intent(inout) :: ppoly_e
164 real,
dimension(:,:),
intent(inout) :: ppoly_coef
166 real,
optional,
intent(in) :: h_neglect
181 hneglect =
hneglect_dflt ;
if (
present(h_neglect)) hneglect = h_neglect
194 u1_r = b *((h0+hneglect)/(h1+hneglect))
198 slope = 2.0 * ( u1 - u0 )
199 if ( abs(u1_r) > abs(slope) )
then
210 u0_l = 3.0 * u0 + 0.5 * u1_r - 2.0 * u0_r
213 exp1 = (u0_r - u0_l) * (u0 - 0.5*(u0_l+u0_r))
214 exp2 = (u0_r - u0_l) * (u0_r - u0_l) / 6.0
216 if ( exp1 > exp2 )
then
217 u0_l = 3.0 * u0 - 2.0 * u0_r
220 if ( exp1 < -exp2 )
then
221 u0_r = 3.0 * u0 - 2.0 * u0_l
228 b = 6.0 * u0 - 4.0 * u0_l - 2.0 * u0_r
229 c = 3.0 * ( u0_r + u0_l - 2.0 * u0 )
248 u1_l = u1_l * ((h1+hneglect)/(h0+hneglect))
251 slope = 2.0 * ( u1 - u0 )
252 if ( abs(u1_l) > abs(slope) )
then
263 u0_r = 3.0 * u1 - 0.5 * u1_l - 2.0 * u0_l
266 exp1 = (u0_r - u0_l) * (u1 - 0.5*(u0_l+u0_r))
267 exp2 = (u0_r - u0_l) * (u0_r - u0_l) / 6.0
269 if ( exp1 > exp2 )
then
270 u0_l = 3.0 * u1 - 2.0 * u0_r
273 if ( exp1 < -exp2 )
then
274 u0_r = 3.0 * u1 - 2.0 * u0_l
281 b = 6.0 * u1 - 4.0 * u0_l - 2.0 * u0_r
282 c = 3.0 * ( u0_r + u0_l - 2.0 * u1 )