portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 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 | |
| 6 | module PLM_functions | |
| 7 | ||
| 8 | implicit none ; private | |
| 9 | ||
| 10 | public PLM_boundary_extrapolation | |
| 11 | public PLM_extrapolate_slope | |
| 12 | public PLM_monotonized_slope | |
| 13 | public PLM_reconstruction | |
| 14 | public PLM_slope_wa | |
| 15 | public PLM_slope_cw | |
| 16 | ||
| 17 | contains | |
| 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. | |
| 22 | 50119172 | real 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 | 50119172 | sigma_r = u_r - u_c |
| 37 | 50119172 | sigma_l = u_c - u_l |
| 38 | ||
| 39 | ! Quasi-second order difference | |
| 40 | 50119172 | 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 | 50119172 | u_min = min( u_l, u_c, u_r ) |
| 44 | 50119172 | u_max = max( u_l, u_c, u_r ) |
| 45 | 50119172 | 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 | 41614192 | 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 | 8504980 | 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 | 50119172 | if (u_c - 0.5*abs(PLM_slope_wa) < u_min .or. u_c + 0.5*abs(PLM_slope_wa) > u_max) then |
| 57 | 14664 | 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 | 50119172 | if (abs(PLM_slope_wa) < 1.E-140) PLM_slope_wa = 0. |
| 64 | ||
| 65 | 50119172 | end 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]. | |
| 69 | 0 | real 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 | 0 | h_cn = h_c + h_neglect |
| 84 | ||
| 85 | ! Side differences | |
| 86 | 0 | sigma_r = u_r - u_c |
| 87 | 0 | 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 | 0 | + ( 2.*h_r + h_c ) / ( h_l + h_cn ) * sigma_l ) |
| 96 | ||
| 97 | ! Limit slope so that reconstructions are bounded by neighbors | |
| 98 | 0 | u_min = min( u_l, u_c, u_r ) |
| 99 | 0 | u_max = max( u_l, u_c, u_r ) |
| 100 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | if (u_c - 0.5*abs(PLM_slope_cw) < u_min .or. u_c + 0.5*abs(PLM_slope_cw) > u_max) then |
| 112 | 0 | 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 | 0 | if (abs(PLM_slope_cw) < 1.E-140) PLM_slope_cw = 0. |
| 119 | ||
| 120 | 0 | end 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]. | |
| 124 | 50119172 | real 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 | 50119172 | almost_two = 2. * ( 1. - epsilon(s_c) ) |
| 137 | ||
| 138 | ! Edge values of neighbors abutting this cell | |
| 139 | 50119172 | e_r = u_l + 0.5*s_l |
| 140 | 50119172 | e_l = u_r - 0.5*s_r |
| 141 | 50119172 | slp = abs(s_c) |
| 142 | ||
| 143 | ! Check that left edge is between right edge of cell to the left and this cell mean | |
| 144 | 50119172 | edge = u_c - 0.5 * s_c |
| 145 | 50119172 | if ( ( edge - e_r ) * ( u_c - edge ) < 0. ) then |
| 146 | 21148604 | edge = 0.5 * ( edge + e_r ) |
| 147 | 21148604 | 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 | 50119172 | edge = u_c + 0.5 * s_c |
| 152 | 50119172 | if ( ( edge - u_c ) * ( e_l - edge ) < 0. ) then |
| 153 | 21148604 | edge = 0.5 * ( edge + e_l ) |
| 154 | 21148604 | slp = min( slp, abs( edge - u_c ) * almost_two ) |
| 155 | endif | |
| 156 | ||
| 157 | 50119172 | PLM_monotonized_slope = sign( slp, s_c ) |
| 158 | ||
| 159 | 50119172 | end 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. | |
| 164 | 726144 | real 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 | 726144 | hl = h_l + h_neglect |
| 176 | 726144 | hc = h_c + h_neglect |
| 177 | ||
| 178 | ! The h2 scheme is used to compute the left edge value | |
| 179 | 726144 | left_edge = (u_l*hc + u_c*hl) / (hl + hc) |
| 180 | ||
| 181 | 726144 | PLM_extrapolate_slope = 2.0 * ( u_c - left_edge ) |
| 182 | ||
| 183 | 726144 | end 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. | |
| 190 | 323492 | subroutine 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 | 646984 | real, dimension(N) :: slp ! The first guess at the normalized tracer slopes [A] |
| 210 | 646984 | real, dimension(N) :: mslp ! The monotonized normalized tracer slopes [A] |
| 211 | ||
| 212 | 323492 | almost_one = 1. - epsilon(slope) |
| 213 | ||
| 214 | ! Loop on interior cells | |
| 215 | 23938408 | do k = 2,N-1 |
| 216 | 23938408 | 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 | 323492 | slp(1) = 0. |
| 221 | 323492 | slp(N) = 0. |
| 222 | ||
| 223 | ! This loop adjusts the slope so that edge values are monotonic. | |
| 224 | 23938408 | do K = 2, N-1 |
| 225 | 23938408 | 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 | 323492 | mslp(1) = 0. |
| 228 | 323492 | mslp(N) = 0. |
| 229 | ||
| 230 | ! Store and return edge values and polynomial coefficients. | |
| 231 | 323492 | edge_values(1,1) = u(1) |
| 232 | 323492 | edge_values(1,2) = u(1) |
| 233 | 323492 | ppoly_coef(1,1) = u(1) |
| 234 | 323492 | ppoly_coef(1,2) = 0. |
| 235 | 23938408 | do k = 2, N-1 |
| 236 | 23614916 | slope = mslp(k) |
| 237 | 23614916 | u_l = u(k) - 0.5 * slope ! Left edge value of cell k |
| 238 | 23614916 | u_r = u(k) + 0.5 * slope ! Right edge value of cell k |
| 239 | ||
| 240 | 23614916 | edge_values(k,1) = u_l |
| 241 | 23614916 | edge_values(k,2) = u_r |
| 242 | 23614916 | ppoly_coef(k,1) = u_l |
| 243 | 23614916 | 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 | 23614916 | edge = ppoly_coef(k,2) + ppoly_coef(k,1) |
| 247 | 23614916 | e_r = u(k+1) - 0.5 * sign( mslp(k+1), slp(k+1) ) |
| 248 | 23938408 | if ( (edge-u(k))*(e_r-edge)<0.) then |
| 249 | 3881 | ppoly_coef(k,2) = ppoly_coef(k,2) * almost_one |
| 250 | endif | |
| 251 | enddo | |
| 252 | 323492 | edge_values(N,1) = u(N) |
| 253 | 323492 | edge_values(N,2) = u(N) |
| 254 | 323492 | ppoly_coef(N,1) = u(N) |
| 255 | 323492 | ppoly_coef(N,2) = 0. |
| 256 | ||
| 257 | 323492 | end 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. | |
| 269 | 0 | subroutine 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 | 0 | slope = - PLM_extrapolate_slope( h(2), h(1), h_neglect, u(2), u(1) ) |
| 285 | ||
| 286 | 0 | edge_values(1,1) = u(1) - 0.5 * slope |
| 287 | 0 | edge_values(1,2) = u(1) + 0.5 * slope |
| 288 | ||
| 289 | 0 | ppoly_coef(1,1) = edge_values(1,1) |
| 290 | 0 | ppoly_coef(1,2) = edge_values(1,2) - edge_values(1,1) |
| 291 | ||
| 292 | ! Extrapolate from N-1 to N to estimate slope | |
| 293 | 0 | slope = PLM_extrapolate_slope( h(N-1), h(N), h_neglect, u(N-1), u(N) ) |
| 294 | ||
| 295 | 0 | edge_values(N,1) = u(N) - 0.5 * slope |
| 296 | 0 | edge_values(N,2) = u(N) + 0.5 * slope |
| 297 | ||
| 298 | 0 | ppoly_coef(N,1) = edge_values(N,1) |
| 299 | 0 | ppoly_coef(N,2) = edge_values(N,2) - edge_values(N,1) |
| 300 | ||
| 301 | 0 | end 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 | ||
| 311 | end module PLM_functions |