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 | !> Provides functions used with the Piecewise-Parabolic-Method in the vertical ALE algorithm. | |
| 6 | module PPM_functions | |
| 7 | ||
| 8 | ! First version was created by Laurent White, June 2008. | |
| 9 | ! Substantially re-factored January 2016. | |
| 10 | ||
| 11 | !! @todo Re-factor PPM_boundary_extrapolation to give round-off safe and | |
| 12 | !! optimization independent results. | |
| 13 | ||
| 14 | use regrid_edge_values, only : bound_edge_values, check_discontinuous_edge_values | |
| 15 | ||
| 16 | implicit none ; private | |
| 17 | ||
| 18 | public PPM_reconstruction, PPM_boundary_extrapolation, PPM_monotonicity | |
| 19 | ||
| 20 | contains | |
| 21 | ||
| 22 | !> Builds quadratic polynomials coefficients from cell mean and edge values. | |
| 23 | 303504 | subroutine PPM_reconstruction( N, h, u, edge_values, ppoly_coef, h_neglect, answer_date) |
| 24 | integer, intent(in) :: N !< Number of cells | |
| 25 | real, dimension(N), intent(in) :: h !< Cell widths [H] | |
| 26 | real, dimension(N), intent(in) :: u !< Cell averages in arbitrary coordinates [A] | |
| 27 | real, dimension(N,2), intent(inout) :: edge_values !< Edge values [A] | |
| 28 | real, dimension(N,3), intent(inout) :: ppoly_coef !< Polynomial coefficients, mainly [A] | |
| 29 | real, intent(in) :: h_neglect !< A negligibly small width [H] | |
| 30 | integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use | |
| 31 | ||
| 32 | ! Local variables | |
| 33 | integer :: k ! Loop index | |
| 34 | real :: edge_l, edge_r ! Edge values (left and right) [A] | |
| 35 | ||
| 36 | ! PPM limiter | |
| 37 | 303504 | call PPM_limiter_standard( N, h, u, edge_values, h_neglect, answer_date=answer_date ) |
| 38 | ||
| 39 | ! Loop over all cells | |
| 40 | 23066304 | do k = 1,N |
| 41 | ||
| 42 | 22762800 | edge_l = edge_values(k,1) |
| 43 | 22762800 | edge_r = edge_values(k,2) |
| 44 | ||
| 45 | ! Store polynomial coefficients | |
| 46 | 22762800 | ppoly_coef(k,1) = edge_l |
| 47 | 22762800 | ppoly_coef(k,2) = 4.0 * ( u(k) - edge_l ) + 2.0 * ( u(k) - edge_r ) |
| 48 | 23066304 | ppoly_coef(k,3) = 3.0 * ( ( edge_r - u(k) ) + ( edge_l - u(k) ) ) |
| 49 | ||
| 50 | enddo | |
| 51 | ||
| 52 | 303504 | end subroutine PPM_reconstruction |
| 53 | ||
| 54 | !> Adjusts edge values using the standard PPM limiter (Colella & Woodward, JCP 1984) | |
| 55 | !! after first checking that the edge values are bounded by neighbors cell averages | |
| 56 | !! and that the edge values are monotonic between cell averages. | |
| 57 | 303504 | subroutine PPM_limiter_standard( N, h, u, edge_values, h_neglect, answer_date ) |
| 58 | integer, intent(in) :: N !< Number of cells | |
| 59 | real, dimension(:), intent(in) :: h !< cell widths (size N) [H] | |
| 60 | real, dimension(:), intent(in) :: u !< cell average properties (size N) [A] | |
| 61 | real, dimension(:,:), intent(inout) :: edge_values !< Potentially modified edge values [A] | |
| 62 | real, intent(in) :: h_neglect !< A negligibly small width [H] | |
| 63 | integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use | |
| 64 | ||
| 65 | ! Local variables | |
| 66 | integer :: k ! Loop index | |
| 67 | real :: u_l, u_c, u_r ! Cell averages (left, center and right) [A] | |
| 68 | real :: edge_l, edge_r ! Edge values (left and right) [A] | |
| 69 | real :: expr1, expr2 ! Temporary expressions [A2] | |
| 70 | ||
| 71 | ! Bound edge values | |
| 72 | 303504 | call bound_edge_values( N, h, u, edge_values, h_neglect, answer_date=answer_date ) |
| 73 | ||
| 74 | ! Make discontinuous edge values monotonic | |
| 75 | 303504 | call check_discontinuous_edge_values( N, u, edge_values ) |
| 76 | ||
| 77 | ! Loop on interior cells to apply the standard | |
| 78 | ! PPM limiter (Colella & Woodward, JCP 84) | |
| 79 | 22459296 | do k = 2,N-1 |
| 80 | ||
| 81 | ! Get cell averages | |
| 82 | 22155792 | u_l = u(k-1) |
| 83 | 22155792 | u_c = u(k) |
| 84 | 22155792 | u_r = u(k+1) |
| 85 | ||
| 86 | 22155792 | edge_l = edge_values(k,1) |
| 87 | 22155792 | edge_r = edge_values(k,2) |
| 88 | ||
| 89 | 22155792 | if ( (u_r - u_c)*(u_c - u_l) <= 0.0) then |
| 90 | ! Flatten extremum | |
| 91 | 4347916 | edge_l = u_c |
| 92 | 4347916 | edge_r = u_c |
| 93 | else | |
| 94 | 17807876 | expr1 = 3.0 * (edge_r - edge_l) * ( (u_c - edge_l) + (u_c - edge_r)) |
| 95 | 17807876 | expr2 = (edge_r - edge_l) * (edge_r - edge_l) |
| 96 | 17807876 | if ( expr1 > expr2 ) then |
| 97 | ! Place extremum at right edge of cell by adjusting left edge value | |
| 98 | 622466 | edge_l = u_c + 2.0 * ( u_c - edge_r ) |
| 99 | 622466 | edge_l = max( min( edge_l, max(u_l, u_c) ), min(u_l, u_c) ) ! In case of round off |
| 100 | 17185410 | elseif ( expr1 < -expr2 ) then |
| 101 | ! Place extremum at left edge of cell by adjusting right edge value | |
| 102 | 505101 | edge_r = u_c + 2.0 * ( u_c - edge_l ) |
| 103 | 505101 | edge_r = max( min( edge_r, max(u_r, u_c) ), min(u_r, u_c) ) ! In case of round off |
| 104 | endif | |
| 105 | endif | |
| 106 | ! This checks that the difference in edge values is representable | |
| 107 | ! and avoids overshoot problems due to round off. | |
| 108 | !### The 1.e-60 needs to have units of [A], so this dimensionally inconsistent. | |
| 109 | 22155792 | if ( abs( edge_r - edge_l )<max(1.e-60,epsilon(u_c)*abs(u_c)) ) then |
| 110 | 4500885 | edge_l = u_c |
| 111 | 4500885 | edge_r = u_c |
| 112 | endif | |
| 113 | ||
| 114 | 22155792 | edge_values(k,1) = edge_l |
| 115 | 22459296 | edge_values(k,2) = edge_r |
| 116 | ||
| 117 | enddo ! end loop on interior cells | |
| 118 | ||
| 119 | ! PCM within boundary cells | |
| 120 | 910512 | edge_values(1,:) = u(1) |
| 121 | 910512 | edge_values(N,:) = u(N) |
| 122 | ||
| 123 | 303504 | end subroutine PPM_limiter_standard |
| 124 | ||
| 125 | !> Adjusts edge values using the original monotonicity constraint (Colella & Woodward, JCP 1984) | |
| 126 | !! Based on hybgen_ppm_coefs | |
| 127 | 0 | subroutine PPM_monotonicity( N, u, edge_values ) |
| 128 | integer, intent(in) :: N !< Number of cells | |
| 129 | real, dimension(:), intent(in) :: u !< cell average properties (size N) [A] | |
| 130 | real, dimension(:,:), intent(inout) :: edge_values !< Potentially modified edge values [A] | |
| 131 | ||
| 132 | ! Local variables | |
| 133 | integer :: k ! Loop index | |
| 134 | real :: a6, da ! Normalized scalar curvature and slope [A] | |
| 135 | ||
| 136 | ! Loop on interior cells to impose monotonicity | |
| 137 | ! Eq. 1.10 of (Colella & Woodward, JCP 84) | |
| 138 | 0 | do k = 2,N-1 |
| 139 | 0 | if (((u(k+1)-u(k))*(u(k)-u(k-1)) <= 0.)) then !local extremum |
| 140 | 0 | edge_values(k,1) = u(k) |
| 141 | 0 | edge_values(k,2) = u(k) |
| 142 | else | |
| 143 | 0 | da = edge_values(k,2)-edge_values(k,1) |
| 144 | 0 | a6 = 6.0*u(k) - 3.0*(edge_values(k,1)+edge_values(k,2)) |
| 145 | 0 | if (da*a6 > da*da) then !peak in right half of zone |
| 146 | 0 | edge_values(k,1) = 3.0*u(k) - 2.0*edge_values(k,2) |
| 147 | 0 | elseif (da*a6 < -da*da) then !peak in left half of zone |
| 148 | 0 | edge_values(k,2) = 3.0*u(k) - 2.0*edge_values(k,1) |
| 149 | endif | |
| 150 | endif | |
| 151 | enddo ! end loop on interior cells | |
| 152 | ||
| 153 | 0 | end subroutine PPM_monotonicity |
| 154 | ||
| 155 | !------------------------------------------------------------------------------ | |
| 156 | !> Reconstruction by parabolas within boundary cells | |
| 157 | 0 | subroutine PPM_boundary_extrapolation( N, h, u, edge_values, ppoly_coef, h_neglect) |
| 158 | !------------------------------------------------------------------------------ | |
| 159 | ! Reconstruction by parabolas within boundary cells. | |
| 160 | ! | |
| 161 | ! The following explanations apply to the left boundary cell. The same | |
| 162 | ! reasoning holds for the right boundary cell. | |
| 163 | ! | |
| 164 | ! A parabola needs to be built in the cell and requires three degrees of | |
| 165 | ! freedom, which are the right edge value and slope and the cell average. | |
| 166 | ! The right edge values and slopes are taken to be that of the neighboring | |
| 167 | ! cell (i.e., the left edge value and slope of the neighboring cell). | |
| 168 | ! The resulting parabola is not necessarily monotonic and the traditional | |
| 169 | ! PPM limiter is used to modify one of the edge values in order to yield | |
| 170 | ! a monotonic parabola. | |
| 171 | ! | |
| 172 | ! N: number of cells in grid | |
| 173 | ! h: thicknesses of grid cells | |
| 174 | ! u: cell averages to use in constructing piecewise polynomials | |
| 175 | ! edge_values : edge values of piecewise polynomials | |
| 176 | ! ppoly_coef : coefficients of piecewise polynomials | |
| 177 | ! | |
| 178 | ! It is assumed that the size of the array 'u' is equal to the number of cells | |
| 179 | ! defining 'grid' and 'ppoly'. No consistency check is performed here. | |
| 180 | !------------------------------------------------------------------------------ | |
| 181 | ||
| 182 | ! Arguments | |
| 183 | integer, intent(in) :: N !< Number of cells | |
| 184 | real, dimension(:), intent(in) :: h !< cell widths (size N) [H] | |
| 185 | real, dimension(:), intent(in) :: u !< cell averages (size N) [A] | |
| 186 | real, dimension(:,:), intent(inout) :: edge_values !< edge values of piecewise polynomials [A] | |
| 187 | real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly [A] | |
| 188 | real, intent(in) :: h_neglect !< A negligibly small width for | |
| 189 | !! the purpose of cell reconstructions [H] | |
| 190 | ||
| 191 | ! Local variables | |
| 192 | integer :: i0, i1 | |
| 193 | real :: u0, u1 ! Average concentrations in the two neighboring cells [A] | |
| 194 | real :: h0, h1 ! Thicknesses of the two neighboring cells [H] | |
| 195 | real :: a, b, c ! An edge value, normalized slope and normalized curvature | |
| 196 | ! of a reconstructed distribution [A] | |
| 197 | real :: u0_l, u0_r ! Edge values of a neighboring cell [A] | |
| 198 | real :: u1_l, u1_r ! Neighboring cell slopes renormalized by the thickness of | |
| 199 | ! the cell being worked on [A] | |
| 200 | real :: slope ! The normalized slope [A] | |
| 201 | real :: exp1, exp2 ! Temporary expressions [A2] | |
| 202 | ||
| 203 | ! ----- Left boundary ----- | |
| 204 | 0 | i0 = 1 |
| 205 | 0 | i1 = 2 |
| 206 | 0 | h0 = h(i0) |
| 207 | 0 | h1 = h(i1) |
| 208 | 0 | u0 = u(i0) |
| 209 | 0 | u1 = u(i1) |
| 210 | ||
| 211 | ! Compute the left edge slope in neighboring cell and express it in | |
| 212 | ! the global coordinate system | |
| 213 | 0 | b = ppoly_coef(i1,2) |
| 214 | 0 | u1_r = b *((h0+h_neglect)/(h1+h_neglect)) ! derivative evaluated at xi = 0.0, |
| 215 | ! expressed w.r.t. xi (local coord. system) | |
| 216 | ||
| 217 | ! Limit the right slope by the PLM limited slope | |
| 218 | 0 | slope = 2.0 * ( u1 - u0 ) |
| 219 | 0 | if ( abs(u1_r) > abs(slope) ) then |
| 220 | 0 | u1_r = slope |
| 221 | endif | |
| 222 | ||
| 223 | ! The right edge value in the boundary cell is taken to be the left | |
| 224 | ! edge value in the neighboring cell | |
| 225 | 0 | u0_r = edge_values(i1,1) |
| 226 | ||
| 227 | ! Given the right edge value and slope, we determine the left | |
| 228 | ! edge value and slope by computing the parabola as determined by | |
| 229 | ! the right edge value and slope and the boundary cell average | |
| 230 | 0 | u0_l = 3.0 * u0 + 0.5 * u1_r - 2.0 * u0_r |
| 231 | ||
| 232 | ! Apply the traditional PPM limiter | |
| 233 | 0 | exp1 = (u0_r - u0_l) * (u0 - 0.5*(u0_l+u0_r)) |
| 234 | 0 | exp2 = (u0_r - u0_l) * (u0_r - u0_l) / 6.0 |
| 235 | ||
| 236 | 0 | if ( exp1 > exp2 ) then |
| 237 | 0 | u0_l = 3.0 * u0 - 2.0 * u0_r |
| 238 | endif | |
| 239 | ||
| 240 | 0 | if ( exp1 < -exp2 ) then |
| 241 | 0 | u0_r = 3.0 * u0 - 2.0 * u0_l |
| 242 | endif | |
| 243 | ||
| 244 | 0 | edge_values(i0,1) = u0_l |
| 245 | 0 | edge_values(i0,2) = u0_r |
| 246 | ||
| 247 | 0 | a = u0_l |
| 248 | 0 | b = 6.0 * u0 - 4.0 * u0_l - 2.0 * u0_r |
| 249 | 0 | c = 3.0 * ( u0_r + u0_l - 2.0 * u0 ) |
| 250 | ||
| 251 | 0 | ppoly_coef(i0,1) = a |
| 252 | 0 | ppoly_coef(i0,2) = b |
| 253 | 0 | ppoly_coef(i0,3) = c |
| 254 | ||
| 255 | ! ----- Right boundary ----- | |
| 256 | 0 | i0 = N-1 |
| 257 | 0 | i1 = N |
| 258 | 0 | h0 = h(i0) |
| 259 | 0 | h1 = h(i1) |
| 260 | 0 | u0 = u(i0) |
| 261 | 0 | u1 = u(i1) |
| 262 | ||
| 263 | ! Compute the right edge slope in neighboring cell and express it in | |
| 264 | ! the global coordinate system | |
| 265 | 0 | b = ppoly_coef(i0,2) |
| 266 | 0 | c = ppoly_coef(i0,3) |
| 267 | 0 | u1_l = (b + 2*c) ! derivative evaluated at xi = 1.0 |
| 268 | 0 | u1_l = u1_l * ((h1+h_neglect)/(h0+h_neglect)) |
| 269 | ||
| 270 | ! Limit the left slope by the PLM limited slope | |
| 271 | 0 | slope = 2.0 * ( u1 - u0 ) |
| 272 | 0 | if ( abs(u1_l) > abs(slope) ) then |
| 273 | 0 | u1_l = slope |
| 274 | endif | |
| 275 | ||
| 276 | ! The left edge value in the boundary cell is taken to be the right | |
| 277 | ! edge value in the neighboring cell | |
| 278 | 0 | u0_l = edge_values(i0,2) |
| 279 | ||
| 280 | ! Given the left edge value and slope, we determine the right | |
| 281 | ! edge value and slope by computing the parabola as determined by | |
| 282 | ! the left edge value and slope and the boundary cell average | |
| 283 | 0 | u0_r = 3.0 * u1 - 0.5 * u1_l - 2.0 * u0_l |
| 284 | ||
| 285 | ! Apply the traditional PPM limiter | |
| 286 | 0 | exp1 = (u0_r - u0_l) * (u1 - 0.5*(u0_l+u0_r)) |
| 287 | 0 | exp2 = (u0_r - u0_l) * (u0_r - u0_l) / 6.0 |
| 288 | ||
| 289 | 0 | if ( exp1 > exp2 ) then |
| 290 | 0 | u0_l = 3.0 * u1 - 2.0 * u0_r |
| 291 | endif | |
| 292 | ||
| 293 | 0 | if ( exp1 < -exp2 ) then |
| 294 | 0 | u0_r = 3.0 * u1 - 2.0 * u0_l |
| 295 | endif | |
| 296 | ||
| 297 | 0 | edge_values(i1,1) = u0_l |
| 298 | 0 | edge_values(i1,2) = u0_r |
| 299 | ||
| 300 | 0 | a = u0_l |
| 301 | 0 | b = 6.0 * u1 - 4.0 * u0_l - 2.0 * u0_r |
| 302 | 0 | c = 3.0 * ( u0_r + u0_l - 2.0 * u1 ) |
| 303 | ||
| 304 | 0 | ppoly_coef(i1,1) = a |
| 305 | 0 | ppoly_coef(i1,2) = b |
| 306 | 0 | ppoly_coef(i1,3) = c |
| 307 | ||
| 308 | 0 | end subroutine PPM_boundary_extrapolation |
| 309 | ||
| 310 | end module PPM_functions |