6 implicit none ;
private
19 integer,
intent(in) :: n
20 real,
dimension(:),
intent(in) :: h
21 real,
dimension(:),
intent(in) :: u
22 real,
dimension(:,:),
intent(inout) :: ppoly_e
24 real,
dimension(:,:),
intent(inout) :: ppoly_coef
26 real,
optional,
intent(in) :: h_neglect
33 real :: h_l, h_c, h_r, h_cn
34 real :: sigma_l, sigma_c, sigma_r
38 real :: u_min, u_max, e_l, e_r, edge
39 real :: almost_one, almost_two
40 real,
dimension(N) :: slp, mslp
43 hneglect =
hneglect_dflt ;
if (
present(h_neglect)) hneglect = h_neglect
45 almost_one = 1. - epsilon(slope)
46 almost_two = 2. * almost_one
52 u_l = u(k-1) ; u_c = u(k) ; u_r = u(k+1)
55 h_l = h(k-1) ; h_c = h(k) ; h_r = h(k+1)
56 h_cn = max( h_c, hneglect )
72 sigma_c = 2.0 * ( u_r - u_l ) * ( h_c / ( h_l + 2.0*h_c + h_r + hneglect) )
74 if ( (sigma_l * sigma_r) > 0.0 )
then
77 u_min = min( u_l, u_c, u_r )
78 u_max = max( u_l, u_c, u_r )
79 slope = sign( min( abs(sigma_c), 2.*min( u_c - u_min, u_max - u_c ) ), sigma_c )
87 u_min = min( u_l, u_c, u_r )
88 u_max = max( u_l, u_c, u_r )
89 if (u_c - 0.5*abs(slope) < u_min .or. u_c + 0.5*abs(slope) > u_max)
then
90 slope = slope * almost_one
94 if (abs(slope) < 1.e-140) slope = 0.
109 ppoly_e(k,1) = u_c - 0.5 * slope
110 ppoly_e(k,2) = u_c + 0.5 * slope
122 u_l = u(k-1) ; u_c = u(k) ; u_r = u(k+1)
125 mslp(k) = abs(slp(k))
126 u_min = min(e_r, u_c)
127 u_max = max(e_r, u_c)
128 edge = u_c - 0.5 * slp(k)
129 if ( ( edge - e_r ) * ( u_c - edge ) < 0. )
then
130 edge = 0.5 * ( edge + e_r )
131 mslp(k) = min( mslp(k), abs( edge - u_c ) * almost_two )
133 edge = u_c + 0.5 * slp(k)
134 if ( ( edge - u_c ) * ( e_l - edge ) < 0. )
then
135 edge = 0.5 * ( edge + e_l )
136 mslp(k) = min( mslp(k), abs( edge - u_c ) * almost_two )
154 ppoly_coef(1,1) = u(1)
157 slope = sign( mslp(k), slp(k) )
158 u_l = u(k) - 0.5 * slope
159 u_r = u(k) + 0.5 * slope
162 u_min = min( u(k-1), u(k) )
163 u_max = max( u(k-1), u(k) )
164 if (u_l<u_min .or. u_l>u_max)
then
165 write(0,*)
'u(k-1)=',u(k-1),
'u(k)=',u(k),
'slp=',slp(k),
'u_l=',u_l
166 stop
'Left edge out of bounds'
168 u_min = min( u(k+1), u(k) )
169 u_max = max( u(k+1), u(k) )
170 if (u_r<u_min .or. u_r>u_max)
then
171 write(0,*)
'u(k)=',u(k),
'u(k+1)=',u(k+1),
'slp=',slp(k),
'u_r=',u_r
172 stop
'Right edge out of bounds'
177 ppoly_coef(k,1) = u_l
178 ppoly_coef(k,2) = ( u_r - u_l )
181 edge = ppoly_coef(k,2) + ppoly_coef(k,1)
182 e_r = u(k+1) - 0.5 * sign( mslp(k+1), slp(k+1) )
183 if ( (edge-u(k))*(e_r-edge)<0.)
then
184 ppoly_coef(k,2) = ppoly_coef(k,2) * almost_one
189 ppoly_coef(n,1) = u(n)
206 integer,
intent(in) :: n
207 real,
dimension(:),
intent(in) :: h
208 real,
dimension(:),
intent(in) :: u
209 real,
dimension(:,:),
intent(inout) :: ppoly_e
211 real,
dimension(:,:),
intent(inout) :: ppoly_coef
213 real,
optional,
intent(in) :: h_neglect
223 hneglect =
hneglect_dflt ;
if (
present(h_neglect)) hneglect = h_neglect
235 ppoly_e(1,2) = (u0*h1 + u1*h0) / (h0 + h1)
239 slope = 2.0 * ( ppoly_e(1,2) - u0 )
241 ppoly_e(1,1) = u0 - 0.5 * slope
242 ppoly_e(1,2) = u0 + 0.5 * slope
244 ppoly_coef(1,1) = ppoly_e(1,1)
245 ppoly_coef(1,2) = ppoly_e(1,2) - ppoly_e(1,1)
250 h0 = h(n-1) + hneglect
257 ppoly_e(n,1) = (u0*h1 + u1*h0) / (h0 + h1)
261 slope = 2.0 * ( u1 - ppoly_e(n,1) )
263 ppoly_e(n,1) = u1 - 0.5 * slope
264 ppoly_e(n,2) = u1 + 0.5 * slope
266 ppoly_coef(n,1) = ppoly_e(n,1)
267 ppoly_coef(n,2) = ppoly_e(n,2) - ppoly_e(n,1)