P3M_functions.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Cubic interpolation functions
6module p3m_functions
7
9
10implicit none ; private
11
14
15contains
16
17!> Set up a piecewise cubic interpolation from cell averages and estimated
18!! edge slopes and values
19!!
20!! Cubic interpolation between edges.
21!!
22!! The edge values and slopes MUST have been estimated prior to calling
23!! this routine.
24!!
25!! It is assumed that the size of the array 'u' is equal to the number of cells
26!! defining 'grid' and 'ppoly'. No consistency check is performed here.
27subroutine p3m_interpolation( N, h, u, edge_values, ppoly_S, ppoly_coef, h_neglect, answer_date )
28 integer, intent(in) :: n !< Number of cells
29 real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
30 real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A]
31 real, dimension(:,:), intent(inout) :: edge_values !< Edge value of polynomial [A]
32 real, dimension(:,:), intent(inout) :: ppoly_s !< Edge slope of polynomial [A H-1].
33 real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A]
34 real, intent(in) :: h_neglect !< A negligibly small width for the
35 !! purpose of cell reconstructions [H]
36 integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use
37
38 ! Call the limiter for p3m, which takes care of everything from
39 ! computing the coefficients of the cubic to monotonizing it.
40 ! This routine could be called directly instead of having to call
41 ! 'P3M_interpolation' first but we do that to provide an homogeneous
42 ! interface.
43 call p3m_limiter( n, h, u, edge_values, ppoly_s, ppoly_coef, h_neglect, &
44 answer_date=answer_date )
45
46end subroutine p3m_interpolation
47
48!> Adust a piecewise cubic reconstruction with a limiter that adjusts the edge
49!! values and slopes
50!!
51!! The p3m limiter operates as follows:
52!!
53!! 1. Edge values are bounded
54!! 2. Discontinuous edge values are systematically averaged
55!! 3. Loop on cells and do the following
56!! a. Build cubic curve
57!! b. Check if cubic curve is monotonic
58!! c. If not, monotonize cubic curve and rebuild it
59!!
60!! Step 3 of the monotonization process leaves all edge values unchanged.
61subroutine p3m_limiter( N, h, u, edge_values, ppoly_S, ppoly_coef, h_neglect, answer_date )
62 integer, intent(in) :: N !< Number of cells
63 real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
64 real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A]
65 real, dimension(:,:), intent(inout) :: edge_values !< Edge value of polynomial [A]
66 real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial [A H-1]
67 real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A]
68 real, intent(in) :: h_neglect !< A negligibly small width for
69 !! the purpose of cell reconstructions [H]
70 integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use
71
72 ! Local variables
73 integer :: k ! loop index
74 logical :: monotonic ! boolean indicating whether the cubic is monotonic
75 real :: u0_l, u0_r ! edge values [A]
76 real :: u1_l, u1_r ! edge slopes [A H-1]
77 real :: u_l, u_c, u_r ! left, center and right cell averages [A]
78 real :: h_l, h_c, h_r ! left, center and right cell widths [H]
79 real :: sigma_l, sigma_c, sigma_r ! left, center and right van Leer slopes [A H-1]
80 real :: slope ! retained PLM slope [A H-1]
81
82 ! 1. Bound edge values (boundary cells are assumed to be local extrema)
83 call bound_edge_values( n, h, u, edge_values, h_neglect, answer_date=answer_date )
84
85 ! 2. Systematically average discontinuous edge values
86 call average_discontinuous_edge_values( n, edge_values )
87
88
89 ! 3. Loop on cells and do the following
90 ! (a) Build cubic curve
91 ! (b) Check if cubic curve is monotonic
92 ! (c) If not, monotonize cubic curve and rebuild it
93 do k = 1,n
94
95 ! Get edge values, edge slopes and cell width
96 u0_l = edge_values(k,1)
97 u0_r = edge_values(k,2)
98 u1_l = ppoly_s(k,1)
99 u1_r = ppoly_s(k,2)
100
101 ! Get cell widths and cell averages (boundary cells are assumed to
102 ! be local extrema for the sake of slopes)
103 u_c = u(k)
104 h_c = h(k)
105
106 if ( k == 1 ) then
107 h_l = h(k)
108 u_l = u(k)
109 else
110 h_l = h(k-1)
111 u_l = u(k-1)
112 endif
113
114 if ( k == n ) then
115 h_r = h(k)
116 u_r = u(k)
117 else
118 h_r = h(k+1)
119 u_r = u(k+1)
120 endif
121
122 ! Compute limited slope
123 sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + h_neglect )
124 sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r + h_neglect )
125 sigma_r = 2.0 * ( u_r - u_c ) / ( h_c + h_neglect )
126
127 if ( (sigma_l * sigma_r) > 0.0 ) then
128 slope = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c )
129 else
130 slope = 0.0
131 endif
132
133 ! If the slopes are small, set them to zero to prevent asymmetric representation near extrema.
134 if ( abs(u1_l*h_c) < epsilon(u_c)*abs(u_c) ) u1_l = 0.0
135 if ( abs(u1_r*h_c) < epsilon(u_c)*abs(u_c) ) u1_r = 0.0
136
137 ! The edge slopes are limited from above by the respective
138 ! one-sided slopes
139 if ( abs(u1_l) > abs(sigma_l) ) then
140 u1_l = sigma_l
141 endif
142
143 if ( abs(u1_r) > abs(sigma_r) ) then
144 u1_r = sigma_r
145 endif
146
147 ! Build cubic interpolant (compute the coefficients)
148 call build_cubic_interpolant( h, k, edge_values, ppoly_s, ppoly_coef )
149
150 ! Check whether cubic is monotonic
151 monotonic = is_cubic_monotonic( ppoly_coef, k )
152
153 ! If cubic is not monotonic, monotonize it by modifiying the
154 ! edge slopes, store the new edge slopes and recompute the
155 ! cubic coefficients
156 if ( .not.monotonic ) then
157 call monotonize_cubic( h_c, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r )
158 endif
159
160 ! Store edge slopes
161 ppoly_s(k,1) = u1_l
162 ppoly_s(k,2) = u1_r
163
164 ! Recompute coefficients of cubic
165 call build_cubic_interpolant( h, k, edge_values, ppoly_s, ppoly_coef )
166
167 enddo ! loop on cells
168
169end subroutine p3m_limiter
170
171
172!> Calculate the edge values and slopes at boundary cells as part of building a
173!! piecewise cubic sub-grid scale profiles
174!!
175!! The following explanations apply to the left boundary cell. The same
176!! reasoning holds for the right boundary cell.
177!!
178!! A cubic needs to be built in the cell and requires four degrees of freedom,
179!! which are the edge values and slopes. The right edge values and slopes are
180!! taken to be that of the neighboring cell (i.e., the left edge value and slope
181!! of the neighboring cell). The left edge value and slope are determined by
182!! computing the parabola based on the cell average and the right edge value
183!! and slope. The resulting cubic is not necessarily monotonic and the slopes
184!! are subsequently modified to yield a monotonic cubic.
185subroutine p3m_boundary_extrapolation( N, h, u, edge_values, ppoly_S, ppoly_coef, &
186 h_neglect, h_neglect_edge )
187 integer, intent(in) :: n !< Number of cells
188 real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
189 real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A]
190 real, dimension(:,:), intent(inout) :: edge_values !< Edge value of polynomial [A]
191 real, dimension(:,:), intent(inout) :: ppoly_s !< Edge slope of polynomial [A H-1]
192 real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A]
193 real, intent(in) :: h_neglect !< A negligibly small width for the
194 !! purpose of cell reconstructions [H]
195 real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose
196 !! of finding edge values [H]. The default is h_neglect.
197 ! Local variables
198 integer :: i0, i1
199 logical :: monotonic ! boolean indicating whether the cubic is monotonic
200 real :: u0, u1 ! Values of u in two adjacent cells [A]
201 real :: h0, h1 ! Values of h in two adjacent cells, plus a smal increment [H]
202 real :: b, c, d ! Temporary variables [A]
203 real :: u0_l, u0_r ! Left and right edge values [A]
204 real :: u1_l, u1_r ! Left and right edge slopes [A H-1]
205 real :: slope ! The cell center slope [A H-1]
206 real :: hneglect_edge ! Negligibly small thickness [H]
207
208 hneglect_edge = h_neglect ; if (present(h_neglect_edge)) hneglect_edge = h_neglect_edge
209
210 ! ----- Left boundary -----
211 i0 = 1
212 i1 = 2
213 h0 = h(i0) + hneglect_edge
214 h1 = h(i1) + hneglect_edge
215 u0 = u(i0)
216 u1 = u(i1)
217
218 ! Compute the left edge slope in neighboring cell and express it in
219 ! the global coordinate system
220 b = ppoly_coef(i1,2)
221 u1_r = b / h1 ! derivative evaluated at xi = 0.0, expressed w.r.t. x
222
223 ! Limit the right slope by the PLM limited slope
224 slope = 2.0 * ( u1 - u0 ) / ( h0 + h_neglect )
225 if ( abs(u1_r) > abs(slope) ) then
226 u1_r = slope
227 endif
228
229 ! The right edge value in the boundary cell is taken to be the left
230 ! edge value in the neighboring cell
231 u0_r = edge_values(i1,1)
232
233 ! Given the right edge value and slope, we determine the left
234 ! edge value and slope by computing the parabola as determined by
235 ! the right edge value and slope and the boundary cell average
236 u0_l = 3.0 * u0 + 0.5 * h0*u1_r - 2.0 * u0_r
237 u1_l = ( - 6.0 * u0 - 2.0 * h0*u1_r + 6.0 * u0_r) / ( h0 + h_neglect )
238
239 ! Check whether the edge values are monotonic. For example, if the left edge
240 ! value is larger than the right edge value while the slope is positive, the
241 ! edge values are inconsistent and we need to modify the left edge value
242 if ( (u0_r-u0_l) * slope < 0.0 ) then
243 u0_l = u0_r
244 u1_l = 0.0
245 u1_r = 0.0
246 endif
247
248 ! Store edge values and slope, build cubic and check monotonicity
249 edge_values(i0,1) = u0_l
250 edge_values(i0,2) = u0_r
251 ppoly_s(i0,1) = u1_l
252 ppoly_s(i0,2) = u1_r
253
254 ! Store edge values and slope, build cubic and check monotonicity
255 call build_cubic_interpolant( h, i0, edge_values, ppoly_s, ppoly_coef )
256 monotonic = is_cubic_monotonic( ppoly_coef, i0 )
257
258 if ( .not.monotonic ) then
259 call monotonize_cubic( h0, u0_l, u0_r, 0.0, slope, slope, u1_l, u1_r )
260
261 ! Rebuild cubic after monotonization
262 ppoly_s(i0,1) = u1_l
263 ppoly_s(i0,2) = u1_r
264 call build_cubic_interpolant( h, i0, edge_values, ppoly_s, ppoly_coef )
265
266 endif
267
268 ! ----- Right boundary -----
269 i0 = n-1
270 i1 = n
271 h0 = h(i0) + hneglect_edge
272 h1 = h(i1) + hneglect_edge
273 u0 = u(i0)
274 u1 = u(i1)
275
276 ! Compute the right edge slope in neighboring cell and express it in
277 ! the global coordinate system
278 b = ppoly_coef(i0,2)
279 c = ppoly_coef(i0,3)
280 d = ppoly_coef(i0,4)
281 u1_l = (b + 2*c + 3*d) / ( h0 + h_neglect ) ! derivative evaluated at xi = 1.0
282
283 ! Limit the left slope by the PLM limited slope
284 slope = 2.0 * ( u1 - u0 ) / ( h1 + h_neglect )
285 if ( abs(u1_l) > abs(slope) ) then
286 u1_l = slope
287 endif
288
289 ! The left edge value in the boundary cell is taken to be the right
290 ! edge value in the neighboring cell
291 u0_l = edge_values(i0,2)
292
293 ! Given the left edge value and slope, we determine the right
294 ! edge value and slope by computing the parabola as determined by
295 ! the left edge value and slope and the boundary cell average
296 u0_r = 3.0 * u1 - 0.5 * h1*u1_l - 2.0 * u0_l
297 u1_r = ( 6.0 * u1 - 2.0 * h1*u1_l - 6.0 * u0_l) / ( h1 + h_neglect )
298
299 ! Check whether the edge values are monotonic. For example, if the right edge
300 ! value is smaller than the left edge value while the slope is positive, the
301 ! edge values are inconsistent and we need to modify the right edge value
302 if ( (u0_r-u0_l) * slope < 0.0 ) then
303 u0_r = u0_l
304 u1_l = 0.0
305 u1_r = 0.0
306 endif
307
308 ! Store edge values and slope, build cubic and check monotonicity
309 edge_values(i1,1) = u0_l
310 edge_values(i1,2) = u0_r
311 ppoly_s(i1,1) = u1_l
312 ppoly_s(i1,2) = u1_r
313
314 call build_cubic_interpolant( h, i1, edge_values, ppoly_s, ppoly_coef )
315 monotonic = is_cubic_monotonic( ppoly_coef, i1 )
316
317 if ( .not.monotonic ) then
318 call monotonize_cubic( h1, u0_l, u0_r, slope, 0.0, slope, u1_l, u1_r )
319
320 ! Rebuild cubic after monotonization
321 ppoly_s(i1,1) = u1_l
322 ppoly_s(i1,2) = u1_r
323 call build_cubic_interpolant( h, i1, edge_values, ppoly_s, ppoly_coef )
324
325 endif
326
327end subroutine p3m_boundary_extrapolation
328
329
330!> Build cubic interpolant in cell k
331!!
332!! Given edge values and edge slopes, compute coefficients of cubic in cell k.
333!!
334!! NOTE: edge values and slopes MUST have been properly calculated prior to
335!! calling this routine.
336subroutine build_cubic_interpolant( h, k, edge_values, ppoly_S, ppoly_coef )
337 real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
338 integer, intent(in) :: k !< The index of the cell to work on
339 real, dimension(:,:), intent(in) :: edge_values !< Edge value of polynomial in arbitrary units [A]
340 real, dimension(:,:), intent(in) :: ppoly_S !< Edge slope of polynomial [A H-1]
341 real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A]
342
343 ! Local variables
344 real :: u0_l, u0_r ! edge values [A]
345 real :: u1_l, u1_r ! edge slopes times the cell width [A]
346 real :: h_c ! cell width [H]
347 real :: a0, a1, a2, a3 ! cubic coefficients [A]
348
349 h_c = h(k)
350
351 u0_l = edge_values(k,1)
352 u0_r = edge_values(k,2)
353
354 u1_l = ppoly_s(k,1) * h_c
355 u1_r = ppoly_s(k,2) * h_c
356
357 a0 = u0_l
358 a1 = u1_l
359 a2 = 3.0 * ( u0_r - u0_l ) - u1_r - 2.0 * u1_l
360 a3 = u1_r + u1_l + 2.0 * ( u0_l - u0_r )
361
362 ppoly_coef(k,1) = a0
363 ppoly_coef(k,2) = a1
364 ppoly_coef(k,3) = a2
365 ppoly_coef(k,4) = a3
366
367end subroutine build_cubic_interpolant
368
369
370!> Check whether the cubic reconstruction in cell k is monotonic
371!!
372!! This function checks whether the cubic curve in cell k is monotonic.
373!! If so, returns 1. Otherwise, returns 0.
374!!
375!! The cubic is monotonic if the first derivative is single-signed in (0,1).
376!! Hence, we check whether the roots (if any) lie inside this interval. If there
377!! is no root or if both roots lie outside this interval, the cubic is monotonic.
378logical function is_cubic_monotonic( ppoly_coef, k )
379 real, dimension(:,:), intent(in) :: ppoly_coef !< Coefficients of cubic polynomial in arbitrary units [A]
380 integer, intent(in) :: k !< The index of the cell to work on
381 ! Local variables
382 real :: a, b, c ! Coefficients of the first derivative of the cubic [A]
383
384 a = ppoly_coef(k,2)
385 b = 2.0 * ppoly_coef(k,3)
386 c = 3.0 * ppoly_coef(k,4)
387
388 ! Look for real roots of the quadratic derivative equation, c*x**2 + b*x + a = 0, in (0, 1)
389 if (b*b - 4.0*a*c <= 0.0) then ! The cubic is monotonic everywhere.
390 is_cubic_monotonic = .true.
391 elseif (a * (a + (b + c)) < 0.0) then ! The derivative changes sign between the endpoints of (0, 1)
392 is_cubic_monotonic = .false.
393 elseif (b * (b + 2.0*c) < 0.0) then ! The second derivative changes sign inside of (0, 1)
394 is_cubic_monotonic = .false.
395 else
396 is_cubic_monotonic = .true.
397 endif
398
399end function is_cubic_monotonic
400
401!> Monotonize a cubic curve by modifying the edge slopes.
402!!
403!! This routine takes care of monotonizing a cubic on [0,1] by modifying the
404!! edge slopes. The edge values are NOT modified. The cubic is entirely
405!! determined by the four degrees of freedom u0_l, u0_r, u1_l and u1_r.
406!!
407!! u1_l and u1_r are the edge slopes expressed in the GLOBAL coordinate system.
408!!
409!! The monotonization occurs as follows.
410!
411!! 1. The edge slopes are set to 0 if they are inconsistent with the limited
412!! PLM slope
413!! 2. We check whether we can find an inflexion point in [0,1]. At most one
414!! inflexion point may exist.
415!! a. If there is no inflexion point, the cubic is monotonic.
416!! b. If there is one inflexion point and it lies outside [0,1], the
417!! cubic is monotonic.
418!! c. If there is one inflexion point and it lies in [0,1] and the slope
419!! at the location of the inflexion point is consistent, the cubic
420!! is monotonic.
421!! d. If the inflexion point lies in [0,1] but the slope is inconsistent,
422!! we go to (3) to shift the location of the inflexion point to the left
423!! or to the right. To the left when the 2nd-order left slope is smaller
424!! than the 2nd order right slope.
425!! 3. Edge slopes are modified to shift the inflexion point, either onto the left
426!! edge or onto the right edge.
427
428subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r )
429 real, intent(in) :: h !< cell width [H]
430 real, intent(in) :: u0_l !< left edge value in arbitrary units [A]
431 real, intent(in) :: u0_r !< right edge value [A]
432 real, intent(in) :: sigma_l !< left 2nd-order slopes [A H-1]
433 real, intent(in) :: sigma_r !< right 2nd-order slopes [A H-1]
434 real, intent(in) :: slope !< limited PLM slope [A H-1]
435 real, intent(inout) :: u1_l !< left edge slopes [A H-1]
436 real, intent(inout) :: u1_r !< right edge slopes [A H-1]
437 ! Local variables
438 logical :: found_ip
439 logical :: inflexion_l ! bool telling if inflex. pt must be on left
440 logical :: inflexion_r ! bool telling if inflex. pt must be on right
441 real :: a1, a2, a3 ! Temporary slopes times the cell width [A]
442 real :: u1_l_tmp ! trial left edge slope [A H-1]
443 real :: u1_r_tmp ! trial right edge slope [A H-1]
444 real :: xi_ip ! location of inflexion point in cell coordinates (0,1) [nondim]
445 real :: slope_ip ! slope at inflexion point times cell width [A]
446
447 found_ip = .false.
448 inflexion_l = .false.
449 inflexion_r = .false.
450
451 ! If the edge slopes are inconsistent w.r.t. the limited PLM slope,
452 ! set them to zero
453 if ( u1_l*slope <= 0.0 ) then
454 u1_l = 0.0
455 endif
456
457 if ( u1_r*slope <= 0.0 ) then
458 u1_r = 0.0
459 endif
460
461 ! Compute the location of the inflexion point, which is the root
462 ! of the second derivative
463 a1 = h * u1_l
464 a2 = 3.0 * ( u0_r - u0_l ) - h*(u1_r + 2.0*u1_l)
465 a3 = h*(u1_r + u1_l) + 2.0*(u0_l - u0_r)
466
467 ! There is a possible root (and inflexion point) only if a3 is nonzero.
468 ! When a3 is zero, the second derivative of the cubic is constant (the
469 ! cubic degenerates into a parabola) and no inflexion point exists.
470 if ( a3 /= 0.0 ) then
471 ! Location of inflexion point
472 xi_ip = - a2 / (3.0 * a3)
473
474 ! If the inflexion point lies in [0,1], change boolean value
475 if ( (xi_ip >= 0.0) .AND. (xi_ip <= 1.0) ) then
476 found_ip = .true.
477 endif
478 endif
479
480 ! When there is an inflexion point within [0,1], check the slope
481 ! to see if it is consistent with the limited PLM slope. If not,
482 ! decide on which side we want to collapse the inflexion point.
483 ! If the inflexion point lies on one of the edges, the cubic is
484 ! guaranteed to be monotonic
485 if ( found_ip ) then
486 slope_ip = a1 + 2.0*a2*xi_ip + 3.0*a3*xi_ip*xi_ip
487
488 ! Check whether slope is consistent
489 if ( slope_ip*slope < 0.0 ) then
490 if ( abs(sigma_l) < abs(sigma_r) ) then
491 inflexion_l = .true.
492 else
493 inflexion_r = .true.
494 endif
495 endif
496 endif ! found_ip
497
498 ! At this point, if the cubic is not monotonic, we know where the
499 ! inflexion point should lie. When the cubic is monotonic, both
500 ! 'inflexion_l' and 'inflexion_r' are false and nothing is to be done.
501
502 ! Move inflexion point on the left
503 if ( inflexion_l ) then
504
505 u1_l_tmp = 1.5*(u0_r-u0_l)/h - 0.5*u1_r
506 u1_r_tmp = 3.0*(u0_r-u0_l)/h - 2.0*u1_l
507
508 if ( (u1_l_tmp*slope < 0.0) .AND. (u1_r_tmp*slope < 0.0) ) then
509
510 u1_l = 0.0
511 u1_r = 3.0 * (u0_r - u0_l) / h
512
513 elseif (u1_l_tmp*slope < 0.0) then
514
515 u1_r = u1_r_tmp
516 u1_l = 1.5*(u0_r - u0_l)/h - 0.5*u1_r
517
518 elseif (u1_r_tmp*slope < 0.0) then
519
520 u1_l = u1_l_tmp
521 u1_r = 3.0*(u0_r - u0_l)/h - 2.0*u1_l
522
523 else
524
525 u1_l = u1_l_tmp
526 u1_r = u1_r_tmp
527
528 endif
529
530 endif ! end treating case with inflexion point on the left
531
532 ! Move inflexion point on the right
533 if ( inflexion_r ) then
534
535 u1_l_tmp = 3.0*(u0_r-u0_l)/h - 2.0*u1_r
536 u1_r_tmp = 1.5*(u0_r-u0_l)/h - 0.5*u1_l
537
538 if ( (u1_l_tmp*slope < 0.0) .AND. (u1_r_tmp*slope < 0.0) ) then
539
540 u1_l = 3.0 * (u0_r - u0_l) / h
541 u1_r = 0.0
542
543 elseif (u1_l_tmp*slope < 0.0) then
544
545 u1_r = u1_r_tmp
546 u1_l = 3.0*(u0_r - u0_l)/h - 2.0*u1_r
547
548 elseif (u1_r_tmp*slope < 0.0) then
549
550 u1_l = u1_l_tmp
551 u1_r = 1.5*(u0_r - u0_l)/h - 0.5*u1_l
552
553 else
554
555 u1_l = u1_l_tmp
556 u1_r = u1_r_tmp
557
558 endif
559
560 endif ! end treating case with inflexion point on the right
561
562 ! Zero out negligibly small slopes.
563 if ( abs(u1_l*h) < epsilon(u0_l) * (abs(u0_l) + abs(u0_r)) ) u1_l = 0.0
564 if ( abs(u1_r*h) < epsilon(u0_l) * (abs(u0_l) + abs(u0_r)) ) u1_r = 0.0
565
566end subroutine monotonize_cubic
567
568!> \namespace p3m_functions
569!!
570!! Date of creation: 2008.06.09
571!! L. White
572!!
573!! This module contains p3m interpolation routines.
574!!
575!! p3m interpolation is performed by estimating the edge values and slopes
576!! and constructing a cubic polynomial. We then make sure that the edge values
577!! are bounded and continuous and we then modify the slopes to get a monotonic
578!! cubic curve.
579
580end module p3m_functions