PLM_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!> Piecewise linear reconstruction functions
6module plm_functions
7
8implicit none ; private
9
14public plm_slope_wa
15public plm_slope_cw
16
17contains
18
19!> Returns a limited PLM slope following White and Adcroft, 2008, in the same arbitrary
20!! units [A] as the input values.
21!! Note that this is not the same as the Colella and Woodward method.
22real elemental pure function plm_slope_wa(h_l, h_c, h_r, h_neglect, u_l, u_c, u_r)
23 real, intent(in) :: h_l !< Thickness of left cell in arbitrary grid thickness units [H]
24 real, intent(in) :: h_c !< Thickness of center cell in arbitrary grid thickness units [H]
25 real, intent(in) :: h_r !< Thickness of right cell in arbitrary grid thickness units [H]
26 real, intent(in) :: h_neglect !< A negligible thickness [H]
27 real, intent(in) :: u_l !< Value of left cell in arbitrary units [A]
28 real, intent(in) :: u_c !< Value of center cell in arbitrary units [A]
29 real, intent(in) :: u_r !< Value of right cell in arbitrary units [A]
30 ! Local variables
31 real :: sigma_l, sigma_c, sigma_r ! Left, central and right slope estimates as
32 ! differences across the cell [A]
33 real :: u_min, u_max ! Minimum and maximum value across cell [A]
34
35 ! Side differences
36 sigma_r = u_r - u_c
37 sigma_l = u_c - u_l
38
39 ! Quasi-second order difference
40 sigma_c = 2.0 * ( u_r - u_l ) * ( h_c / ( h_l + 2.0*h_c + h_r + h_neglect) )
41
42 ! Limit slope so that reconstructions are bounded by neighbors
43 u_min = min( u_l, u_c, u_r )
44 u_max = max( u_l, u_c, u_r )
45 if ( (sigma_l * sigma_r) > 0.0 ) then
46 ! This limits the slope so that the edge values are bounded by the
47 ! two cell averages spanning the edge.
48 plm_slope_wa = sign( min( abs(sigma_c), 2.*min( u_c - u_min, u_max - u_c ) ), sigma_c )
49 else
50 ! Extrema in the mean values require a PCM reconstruction avoid generating
51 ! larger extreme values.
52 plm_slope_wa = 0.0
53 endif
54
55 ! This block tests to see if roundoff causes edge values to be out of bounds
56 if (u_c - 0.5*abs(plm_slope_wa) < u_min .or. u_c + 0.5*abs(plm_slope_wa) > u_max) then
57 plm_slope_wa = plm_slope_wa * ( 1. - epsilon(plm_slope_wa) )
58 endif
59
60 ! An attempt to avoid inconsistency when the values become unrepresentable.
61 ! ### The following 1.E-140 is dimensionally inconsistent. A newer version of
62 ! PLM is progress that will avoid the need for such rounding.
63 if (abs(plm_slope_wa) < 1.e-140) plm_slope_wa = 0.
64
65end function plm_slope_wa
66
67!> Returns a limited PLM slope following Colella and Woodward 1984, in the same
68!! arbitrary units as the input values [A].
69real elemental pure function plm_slope_cw(h_l, h_c, h_r, h_neglect, u_l, u_c, u_r)
70 real, intent(in) :: h_l !< Thickness of left cell in arbitrary grid thickness units [H]
71 real, intent(in) :: h_c !< Thickness of center cell in arbitrary grid thickness units [H]
72 real, intent(in) :: h_r !< Thickness of right cell in arbitrary grid thickness units [H]
73 real, intent(in) :: h_neglect !< A negligible thickness [H]
74 real, intent(in) :: u_l !< Value of left cell in arbitrary units [A]
75 real, intent(in) :: u_c !< Value of center cell in arbitrary units [A]
76 real, intent(in) :: u_r !< Value of right cell in arbitrary units [A]
77 ! Local variables
78 real :: sigma_l, sigma_c, sigma_r ! Left, central and right slope estimates as
79 ! differences across the cell [A]
80 real :: u_min, u_max ! Minimum and maximum value across cell [A]
81 real :: h_cn ! Thickness of center cell [H]
82
83 h_cn = h_c + h_neglect
84
85 ! Side differences
86 sigma_r = u_r - u_c
87 sigma_l = u_c - u_l
88
89 ! This is the second order slope given by equation 1.7 of
90 ! Piecewise Parabolic Method, Colella and Woodward (1984),
91 ! http://dx.doi.org/10.1016/0021-991(84)90143-8.
92 ! For uniform resolution it simplifies to ( u_r - u_l )/2 .
93 sigma_c = ( h_c / ( h_cn + ( h_l + h_r ) ) ) * ( &
94 ( 2.*h_l + h_c ) / ( h_r + h_cn ) * sigma_r &
95 + ( 2.*h_r + h_c ) / ( h_l + h_cn ) * sigma_l )
96
97 ! Limit slope so that reconstructions are bounded by neighbors
98 u_min = min( u_l, u_c, u_r )
99 u_max = max( u_l, u_c, u_r )
100 if ( (sigma_l * sigma_r) > 0.0 ) then
101 ! This limits the slope so that the edge values are bounded by the
102 ! two cell averages spanning the edge.
103 plm_slope_cw = sign( min( abs(sigma_c), 2.*min( u_c - u_min, u_max - u_c ) ), sigma_c )
104 else
105 ! Extrema in the mean values require a PCM reconstruction avoid generating
106 ! larger extreme values.
107 plm_slope_cw = 0.0
108 endif
109
110 ! This block tests to see if roundoff causes edge values to be out of bounds
111 if (u_c - 0.5*abs(plm_slope_cw) < u_min .or. u_c + 0.5*abs(plm_slope_cw) > u_max) then
112 plm_slope_cw = plm_slope_cw * ( 1. - epsilon(plm_slope_cw) )
113 endif
114
115 ! An attempt to avoid inconsistency when the values become unrepresentable.
116 ! ### The following 1.E-140 is dimensionally inconsistent. A newer version of
117 ! PLM is progress that will avoid the need for such rounding.
118 if (abs(plm_slope_cw) < 1.e-140) plm_slope_cw = 0.
119
120end function plm_slope_cw
121
122!> Returns a limited PLM slope following Colella and Woodward 1984, in the same
123!! arbitrary units as the input values [A].
124real elemental pure function plm_monotonized_slope(u_l, u_c, u_r, s_l, s_c, s_r)
125 real, intent(in) :: u_l !< Value of left cell in arbitrary units [A]
126 real, intent(in) :: u_c !< Value of center cell in arbitrary units [A]
127 real, intent(in) :: u_r !< Value of right cell in arbitrary units [A]
128 real, intent(in) :: s_l !< PLM slope of left cell [A]
129 real, intent(in) :: s_c !< PLM slope of center cell [A]
130 real, intent(in) :: s_r !< PLM slope of right cell [A]
131 ! Local variables
132 real :: e_r, e_l, edge ! Right, left and temporary edge values [A]
133 real :: almost_two ! The number 2, almost [nondim]
134 real :: slp ! Magnitude of PLM central slope [A]
135
136 almost_two = 2. * ( 1. - epsilon(s_c) )
137
138 ! Edge values of neighbors abutting this cell
139 e_r = u_l + 0.5*s_l
140 e_l = u_r - 0.5*s_r
141 slp = abs(s_c)
142
143 ! Check that left edge is between right edge of cell to the left and this cell mean
144 edge = u_c - 0.5 * s_c
145 if ( ( edge - e_r ) * ( u_c - edge ) < 0. ) then
146 edge = 0.5 * ( edge + e_r )
147 slp = min( slp, abs( edge - u_c ) * almost_two )
148 endif
149
150 ! Check that right edge is between left edge of cell to the right and this cell mean
151 edge = u_c + 0.5 * s_c
152 if ( ( edge - u_c ) * ( e_l - edge ) < 0. ) then
153 edge = 0.5 * ( edge + e_l )
154 slp = min( slp, abs( edge - u_c ) * almost_two )
155 endif
156
157 plm_monotonized_slope = sign( slp, s_c )
158
159end function plm_monotonized_slope
160
161!> Returns a PLM slope using h2 extrapolation from a cell to the left, in the same
162!! arbitrary units as the input values [A].
163!! Use the negative to extrapolate from the cell to the right.
164real elemental pure function plm_extrapolate_slope(h_l, h_c, h_neglect, u_l, u_c)
165 real, intent(in) :: h_l !< Thickness of left cell in arbitrary grid thickness units [H]
166 real, intent(in) :: h_c !< Thickness of center cell in arbitrary grid thickness units [H]
167 real, intent(in) :: h_neglect !< A negligible thickness [H]
168 real, intent(in) :: u_l !< Value of left cell in arbitrary units [A]
169 real, intent(in) :: u_c !< Value of center cell in arbitrary units [A]
170 ! Local variables
171 real :: left_edge ! Left edge value [A]
172 real :: hl, hc ! Left and central cell thicknesses [H]
173
174 ! Avoid division by zero for vanished cells
175 hl = h_l + h_neglect
176 hc = h_c + h_neglect
177
178 ! The h2 scheme is used to compute the left edge value
179 left_edge = (u_l*hc + u_c*hl) / (hl + hc)
180
181 plm_extrapolate_slope = 2.0 * ( u_c - left_edge )
182
183end function plm_extrapolate_slope
184
185
186!> Reconstruction by linear polynomials within each cell
187!!
188!! It is assumed that the size of the array 'u' is equal to the number of cells
189!! defining 'grid' and 'ppoly'. No consistency check is performed here.
190subroutine plm_reconstruction( N, h, u, edge_values, ppoly_coef, h_neglect )
191 integer, intent(in) :: n !< Number of cells
192 real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
193 real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A]
194 real, dimension(:,:), intent(inout) :: edge_values !< edge values of piecewise polynomials,
195 !! with the same units as u [A].
196 real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly
197 !! with the same units as u [A].
198 real, intent(in) :: h_neglect !< A negligibly small width for
199 !! the purpose of cell reconstructions
200 !! in the same units as h [H]
201
202 ! Local variables
203 integer :: k ! loop index
204 real :: u_l, u_r ! left and right cell averages [A]
205 real :: slope ! retained PLM slope for a normalized cell width [A]
206 real :: e_r ! The edge value in the neighboring cell [A]
207 real :: edge ! The projected edge value in the cell [A]
208 real :: almost_one ! A value that is slightly smaller than 1 [nondim]
209 real, dimension(N) :: slp ! The first guess at the normalized tracer slopes [A]
210 real, dimension(N) :: mslp ! The monotonized normalized tracer slopes [A]
211
212 almost_one = 1. - epsilon(slope)
213
214 ! Loop on interior cells
215 do k = 2,n-1
216 slp(k) = plm_slope_wa(h(k-1), h(k), h(k+1), h_neglect, u(k-1), u(k), u(k+1))
217 enddo ! end loop on interior cells
218
219 ! Boundary cells use PCM. Extrapolation is handled after monotonization.
220 slp(1) = 0.
221 slp(n) = 0.
222
223 ! This loop adjusts the slope so that edge values are monotonic.
224 do k = 2, n-1
225 mslp(k) = plm_monotonized_slope( u(k-1), u(k), u(k+1), slp(k-1), slp(k), slp(k+1) )
226 enddo ! end loop on interior cells
227 mslp(1) = 0.
228 mslp(n) = 0.
229
230 ! Store and return edge values and polynomial coefficients.
231 edge_values(1,1) = u(1)
232 edge_values(1,2) = u(1)
233 ppoly_coef(1,1) = u(1)
234 ppoly_coef(1,2) = 0.
235 do k = 2, n-1
236 slope = mslp(k)
237 u_l = u(k) - 0.5 * slope ! Left edge value of cell k
238 u_r = u(k) + 0.5 * slope ! Right edge value of cell k
239
240 edge_values(k,1) = u_l
241 edge_values(k,2) = u_r
242 ppoly_coef(k,1) = u_l
243 ppoly_coef(k,2) = ( u_r - u_l )
244 ! Check to see if this evaluation of the polynomial at x=1 would be
245 ! monotonic w.r.t. the next cell's edge value. If not, scale back!
246 edge = ppoly_coef(k,2) + ppoly_coef(k,1)
247 e_r = u(k+1) - 0.5 * sign( mslp(k+1), slp(k+1) )
248 if ( (edge-u(k))*(e_r-edge)<0.) then
249 ppoly_coef(k,2) = ppoly_coef(k,2) * almost_one
250 endif
251 enddo
252 edge_values(n,1) = u(n)
253 edge_values(n,2) = u(n)
254 ppoly_coef(n,1) = u(n)
255 ppoly_coef(n,2) = 0.
256
257end subroutine plm_reconstruction
258
259
260!> Reconstruction by linear polynomials within boundary cells
261!!
262!! The left and right edge values in the left and right boundary cells,
263!! respectively, are estimated using a linear extrapolation within the cells.
264!!
265!! This extrapolation is EXACT when the underlying profile is linear.
266!!
267!! It is assumed that the size of the array 'u' is equal to the number of cells
268!! defining 'grid' and 'ppoly'. No consistency check is performed here.
269subroutine plm_boundary_extrapolation( N, h, u, edge_values, ppoly_coef, h_neglect )
270 integer, intent(in) :: n !< Number of cells
271 real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
272 real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A]
273 real, dimension(:,:), intent(inout) :: edge_values !< edge values of piecewise polynomials,
274 !! with the same units as u [A].
275 real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly
276 !! with the same units as u [A].
277 real, intent(in) :: h_neglect !< A negligibly small width for
278 !! the purpose of cell reconstructions
279 !! in the same units as h [H]
280 ! Local variables
281 real :: slope ! retained PLM slope for a normalized cell width [A]
282
283 ! Extrapolate from 2 to 1 to estimate slope
284 slope = - plm_extrapolate_slope( h(2), h(1), h_neglect, u(2), u(1) )
285
286 edge_values(1,1) = u(1) - 0.5 * slope
287 edge_values(1,2) = u(1) + 0.5 * slope
288
289 ppoly_coef(1,1) = edge_values(1,1)
290 ppoly_coef(1,2) = edge_values(1,2) - edge_values(1,1)
291
292 ! Extrapolate from N-1 to N to estimate slope
293 slope = plm_extrapolate_slope( h(n-1), h(n), h_neglect, u(n-1), u(n) )
294
295 edge_values(n,1) = u(n) - 0.5 * slope
296 edge_values(n,2) = u(n) + 0.5 * slope
297
298 ppoly_coef(n,1) = edge_values(n,1)
299 ppoly_coef(n,2) = edge_values(n,2) - edge_values(n,1)
300
301end subroutine plm_boundary_extrapolation
302
303!> \namespace plm_functions
304!!
305!! Date of creation: 2008.06.06
306!! L. White
307!!
308!! This module contains routines that handle one-dimensional finite volume
309!! reconstruction using the piecewise linear method (PLM).
310
311end module plm_functions