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 | !> Edge value estimation for high-order reconstruction | |
| 6 | module regrid_edge_values | |
| 7 | ||
| 8 | use MOM_error_handler, only : MOM_error, FATAL | |
| 9 | use regrid_solvers, only : solve_linear_system, linear_solver | |
| 10 | use regrid_solvers, only : solve_tridiagonal_system, solve_diag_dominant_tridiag | |
| 11 | use polynomial_functions, only : evaluation_polynomial | |
| 12 | ||
| 13 | implicit none ; private | |
| 14 | ||
| 15 | ! ----------------------------------------------------------------------------- | |
| 16 | ! The following routines are visible to the outside world | |
| 17 | ! ----------------------------------------------------------------------------- | |
| 18 | public bound_edge_values, average_discontinuous_edge_values, check_discontinuous_edge_values | |
| 19 | public edge_values_explicit_h2, edge_values_explicit_h4, edge_values_explicit_h4cw | |
| 20 | public edge_values_implicit_h4, edge_values_implicit_h6 | |
| 21 | public edge_slopes_implicit_h3, edge_slopes_implicit_h5 | |
| 22 | ||
| 23 | ! The following parameter is used to avoid singular matrices for boundary | |
| 24 | ! extrapolation. It is needed only in the case where thicknesses vanish | |
| 25 | ! to a small enough values such that the eigenvalues of the matrix can not | |
| 26 | ! be separated. | |
| 27 | real, parameter :: hMinFrac = 1.e-5 !< A minimum fraction for min(h)/sum(h) [nondim] | |
| 28 | ||
| 29 | contains | |
| 30 | ||
| 31 | !> Bound edge values by neighboring cell averages | |
| 32 | !! | |
| 33 | !! In this routine, we loop on all cells to bound their left and right | |
| 34 | !! edge values by the cell averages. That is, the left edge value must lie | |
| 35 | !! between the left cell average and the central cell average. A similar | |
| 36 | !! reasoning applies to the right edge values. | |
| 37 | !! | |
| 38 | !! Both boundary edge values are set equal to the boundary cell averages. | |
| 39 | !! Any extrapolation scheme is applied after this routine has been called. | |
| 40 | !! Therefore, boundary cells are treated as if they were local extrema. | |
| 41 | 303504 | subroutine bound_edge_values( N, h, u, edge_val, h_neglect, answer_date ) |
| 42 | integer, intent(in) :: N !< Number of cells | |
| 43 | real, dimension(N), intent(in) :: h !< cell widths [H] | |
| 44 | real, dimension(N), intent(in) :: u !< cell average properties in arbitrary units [A] | |
| 45 | real, dimension(N,2), intent(inout) :: edge_val !< Potentially modified edge values [A]; the | |
| 46 | !! second index is for the two edges of each cell. | |
| 47 | real, intent(in) :: h_neglect !< A negligibly small width [H] | |
| 48 | integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use | |
| 49 | ||
| 50 | ! Local variables | |
| 51 | real :: sigma_l, sigma_c, sigma_r ! left, center and right van Leer slopes [A H-1] or [A] | |
| 52 | real :: slope_x_h ! retained PLM slope times half grid step [A] | |
| 53 | logical :: use_2018_answers ! If true use older, less accurate expressions. | |
| 54 | integer :: k, km1, kp1 ! Loop index and the values to either side. | |
| 55 | ||
| 56 | 303504 | use_2018_answers = .true. ; if (present(answer_date)) use_2018_answers = (answer_date < 20190101) |
| 57 | ||
| 58 | ! Loop on cells to bound edge value | |
| 59 | 23066304 | do k = 1,N |
| 60 | ||
| 61 | ! For the sake of bounding boundary edge values, the left neighbor of the left boundary cell | |
| 62 | ! is assumed to be the same as the left boundary cell and the right neighbor of the right | |
| 63 | ! boundary cell is assumed to be the same as the right boundary cell. This effectively makes | |
| 64 | ! boundary cells look like extrema. | |
| 65 | 22762800 | km1 = max(1,k-1) ; kp1 = min(k+1,N) |
| 66 | ||
| 67 | 22762800 | slope_x_h = 0.0 |
| 68 | 22762800 | if (use_2018_answers) then |
| 69 | 0 | sigma_l = 2.0 * ( u(k) - u(km1) ) / ( h(k) + h_neglect ) |
| 70 | 0 | sigma_c = 2.0 * ( u(kp1) - u(km1) ) / ( h(km1) + 2.0*h(k) + h(kp1) + h_neglect ) |
| 71 | 0 | sigma_r = 2.0 * ( u(kp1) - u(k) ) / ( h(k) + h_neglect ) |
| 72 | ||
| 73 | ! The limiter is used in the local coordinate system to each cell, so for convenience store | |
| 74 | ! the slope times a half grid spacing. (See White and Adcroft JCP 2008 Eqs 19 and 20) | |
| 75 | 0 | if ( (sigma_l * sigma_r) > 0.0 ) & |
| 76 | 0 | slope_x_h = 0.5 * h(k) * sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c ) |
| 77 | 22762800 | elseif ( ((h(km1) + h(kp1)) + 2.0*h(k)) > 0.0 ) then |
| 78 | 22762800 | sigma_l = ( u(k) - u(km1) ) |
| 79 | 22762800 | sigma_c = ( u(kp1) - u(km1) ) * ( h(k) / ((h(km1) + h(kp1)) + 2.0*h(k)) ) |
| 80 | 22762800 | sigma_r = ( u(kp1) - u(k) ) |
| 81 | ||
| 82 | ! The limiter is used in the local coordinate system to each cell, so for convenience store | |
| 83 | ! the slope times a half grid spacing. (See White and Adcroft JCP 2008 Eqs 19 and 20) | |
| 84 | 22762800 | if ( (sigma_l * sigma_r) > 0.0 ) & |
| 85 | 17807876 | slope_x_h = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c ) |
| 86 | endif | |
| 87 | ||
| 88 | ! Limit the edge values | |
| 89 | 22762800 | if ( (u(km1)-edge_val(k,1)) * (edge_val(k,1)-u(k)) < 0.0 ) then |
| 90 | 5112335 | edge_val(k,1) = u(k) - sign( min( abs(slope_x_h), abs(edge_val(k,1)-u(k)) ), slope_x_h ) |
| 91 | endif | |
| 92 | ||
| 93 | 22762800 | if ( (u(kp1)-edge_val(k,2)) * (edge_val(k,2)-u(k)) < 0.0 ) then |
| 94 | 5078345 | edge_val(k,2) = u(k) + sign( min( abs(slope_x_h), abs(edge_val(k,2)-u(k)) ), slope_x_h ) |
| 95 | endif | |
| 96 | ||
| 97 | ! Finally bound by neighboring cell means in case of roundoff | |
| 98 | 22762800 | edge_val(k,1) = max( min( edge_val(k,1), max(u(km1), u(k)) ), min(u(km1), u(k)) ) |
| 99 | 23066304 | edge_val(k,2) = max( min( edge_val(k,2), max(u(kp1), u(k)) ), min(u(kp1), u(k)) ) |
| 100 | ||
| 101 | enddo ! loop on interior edges | |
| 102 | ||
| 103 | 303504 | end subroutine bound_edge_values |
| 104 | ||
| 105 | !> Replace discontinuous collocated edge values with their average | |
| 106 | !! | |
| 107 | !! For each interior edge, check whether the edge values are discontinuous. | |
| 108 | !! If so, compute the average and replace the edge values by the average. | |
| 109 | 0 | subroutine average_discontinuous_edge_values( N, edge_val ) |
| 110 | integer, intent(in) :: N !< Number of cells | |
| 111 | real, dimension(N,2), intent(inout) :: edge_val !< Edge values that may be modified [A]; the | |
| 112 | !! second index is for the two edges of each cell. | |
| 113 | ! Local variables | |
| 114 | integer :: k ! loop index | |
| 115 | real :: u0_avg ! avg value at given edge [A] | |
| 116 | ||
| 117 | ! Loop on interior edges | |
| 118 | 0 | do k = 1,N-1 |
| 119 | ! Compare edge values on the right and left sides of the edge | |
| 120 | 0 | if ( edge_val(k,2) /= edge_val(k+1,1) ) then |
| 121 | 0 | u0_avg = 0.5 * ( edge_val(k,2) + edge_val(k+1,1) ) |
| 122 | 0 | edge_val(k,2) = u0_avg |
| 123 | 0 | edge_val(k+1,1) = u0_avg |
| 124 | endif | |
| 125 | ||
| 126 | enddo ! end loop on interior edges | |
| 127 | ||
| 128 | 0 | end subroutine average_discontinuous_edge_values |
| 129 | ||
| 130 | !> Check discontinuous edge values and replace them with their average if not monotonic | |
| 131 | !! | |
| 132 | !! For each interior edge, check whether the edge values are discontinuous. | |
| 133 | !! If so and if they are not monotonic, replace each edge value by their average. | |
| 134 | 303504 | subroutine check_discontinuous_edge_values( N, u, edge_val ) |
| 135 | integer, intent(in) :: N !< Number of cells | |
| 136 | real, dimension(N), intent(in) :: u !< cell averages in arbitrary units [A] | |
| 137 | real, dimension(N,2), intent(inout) :: edge_val !< Cell edge values [A]; the | |
| 138 | !! second index is for the two edges of each cell. | |
| 139 | ! Local variables | |
| 140 | integer :: k ! loop index | |
| 141 | real :: u0_avg ! avg value at given edge [A] | |
| 142 | ||
| 143 | 22762800 | do k = 1,N-1 |
| 144 | 22762800 | if ( (edge_val(k+1,1) - edge_val(k,2)) * (u(k+1) - u(k)) < 0.0 ) then |
| 145 | 900210 | u0_avg = 0.5 * ( edge_val(k,2) + edge_val(k+1,1) ) |
| 146 | 900210 | u0_avg = max( min( u0_avg, max(u(k), u(k+1)) ), min(u(k), u(k+1)) ) |
| 147 | 900210 | edge_val(k,2) = u0_avg |
| 148 | 900210 | edge_val(k+1,1) = u0_avg |
| 149 | endif | |
| 150 | enddo ! end loop on interior edges | |
| 151 | ||
| 152 | 303504 | end subroutine check_discontinuous_edge_values |
| 153 | ||
| 154 | ||
| 155 | !> Compute h2 edge values (explicit second order accurate) | |
| 156 | !! in the same units as u. | |
| 157 | ! | |
| 158 | !! Compute edge values based on second-order explicit estimates. | |
| 159 | !! These estimates are based on a straight line spanning two cells and evaluated | |
| 160 | !! at the location of the middle edge. An interpolant spanning cells | |
| 161 | !! k-1 and k is evaluated at edge k-1/2. The estimate for each edge is unique. | |
| 162 | !! | |
| 163 | !! k-1 k | |
| 164 | !! ..--o------o------o--.. | |
| 165 | !! k-1/2 | |
| 166 | !! | |
| 167 | !! Boundary edge values are set to be equal to the boundary cell averages. | |
| 168 | 0 | subroutine edge_values_explicit_h2( N, h, u, edge_val ) |
| 169 | integer, intent(in) :: N !< Number of cells | |
| 170 | real, dimension(N), intent(in) :: h !< cell widths [H] | |
| 171 | real, dimension(N), intent(in) :: u !< cell average properties in arbitrary units [A] | |
| 172 | real, dimension(N,2), intent(inout) :: edge_val !< Returned edge values [A]; the | |
| 173 | !! second index is for the two edges of each cell. | |
| 174 | ||
| 175 | ! Local variables | |
| 176 | integer :: k ! loop index | |
| 177 | ||
| 178 | ! Boundary edge values are simply equal to the boundary cell averages | |
| 179 | 0 | edge_val(1,1) = u(1) |
| 180 | 0 | edge_val(N,2) = u(N) |
| 181 | ||
| 182 | 0 | do k = 2,N |
| 183 | ! Compute left edge value | |
| 184 | 0 | if (h(k-1) + h(k) == 0.0) then ! Avoid singularities when h0+h1=0 |
| 185 | 0 | edge_val(k,1) = 0.5 * (u(k-1) + u(k)) |
| 186 | else | |
| 187 | 0 | edge_val(k,1) = ( u(k-1)*h(k) + u(k)*h(k-1) ) / ( h(k-1) + h(k) ) |
| 188 | endif | |
| 189 | ||
| 190 | ! Left edge value of the current cell is equal to right edge value of left cell | |
| 191 | 0 | edge_val(k-1,2) = edge_val(k,1) |
| 192 | enddo | |
| 193 | ||
| 194 | 0 | end subroutine edge_values_explicit_h2 |
| 195 | ||
| 196 | !> Compute h4 edge values (explicit fourth order accurate) | |
| 197 | !! in the same units as u. | |
| 198 | !! | |
| 199 | !! Compute edge values based on fourth-order explicit estimates. | |
| 200 | !! These estimates are based on a cubic interpolant spanning four cells | |
| 201 | !! and evaluated at the location of the middle edge. An interpolant spanning | |
| 202 | !! cells i-2, i-1, i and i+1 is evaluated at edge i-1/2. The estimate for | |
| 203 | !! each edge is unique. | |
| 204 | !! | |
| 205 | !! i-2 i-1 i i+1 | |
| 206 | !! ..--o------o------o------o------o--.. | |
| 207 | !! i-1/2 | |
| 208 | !! | |
| 209 | !! The first two edge values are estimated by evaluating the first available | |
| 210 | !! cubic interpolant, i.e., the interpolant spanning cells 1, 2, 3 and 4. | |
| 211 | !! Similarly, the last two edge values are estimated by evaluating the last | |
| 212 | !! available interpolant. | |
| 213 | !! | |
| 214 | !! For this fourth-order scheme, at least four cells must exist. | |
| 215 | 0 | subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect, answer_date ) |
| 216 | integer, intent(in) :: N !< Number of cells | |
| 217 | real, dimension(N), intent(in) :: h !< cell widths [H] | |
| 218 | real, dimension(N), intent(in) :: u !< cell average properties in arbitrary units [A] | |
| 219 | real, dimension(N,2), intent(inout) :: edge_val !< Returned edge values [A]; the second index | |
| 220 | !! is for the two edges of each cell. | |
| 221 | real, intent(in) :: h_neglect !< A negligibly small width [H] | |
| 222 | integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use | |
| 223 | ||
| 224 | ! Local variables | |
| 225 | real :: h0, h1, h2, h3 ! temporary thicknesses [H] | |
| 226 | real :: h_min ! A minimal cell width [H] | |
| 227 | real :: f1 ! An auxiliary variable [H] | |
| 228 | real :: f2 ! An auxiliary variable [A H] | |
| 229 | real :: f3 ! An auxiliary variable [H-1] | |
| 230 | real :: et1, et2, et3 ! terms the expression for edge values [A H] | |
| 231 | real :: I_h12 ! The inverse of the sum of the two central thicknesses [H-1] | |
| 232 | real :: I_h012, I_h123 ! Inverses of sums of three successive thicknesses [H-1] | |
| 233 | real :: I_den_et2, I_den_et3 ! Inverses of denominators in edge value terms [H-2] | |
| 234 | real, dimension(5) :: x ! Coordinate system with 0 at edges [H] | |
| 235 | real, dimension(4) :: dz ! A temporary array of limited layer thicknesses [H] | |
| 236 | real, dimension(4) :: u_tmp ! A temporary array of cell average properties [A] | |
| 237 | real, parameter :: C1_12 = 1.0 / 12.0 ! A rational constant [nondim] | |
| 238 | real :: dx ! Difference of successive values of x [H] | |
| 239 | real, dimension(4,4) :: A ! Differences in successive positions raised to various powers, | |
| 240 | ! in units that vary with the second (j) index as [H^j] | |
| 241 | real, dimension(4) :: B ! The right hand side of the system to solve for C [A H] | |
| 242 | real, dimension(4) :: C ! The coefficients of a fit polynomial in units that vary | |
| 243 | ! with the index (j) as [A H^(j-1)] | |
| 244 | integer :: i, j | |
| 245 | logical :: use_2018_answers ! If true use older, less accurate expressions. | |
| 246 | ||
| 247 | 0 | use_2018_answers = .true. ; if (present(answer_date)) use_2018_answers = (answer_date < 20190101) |
| 248 | ||
| 249 | ! Loop on interior cells | |
| 250 | 0 | do i = 3,N-1 |
| 251 | ||
| 252 | 0 | h0 = h(i-2) |
| 253 | 0 | h1 = h(i-1) |
| 254 | 0 | h2 = h(i) |
| 255 | 0 | h3 = h(i+1) |
| 256 | ||
| 257 | ! Avoid singularities when consecutive pairs of h vanish | |
| 258 | 0 | if (h0+h1==0.0 .or. h1+h2==0.0 .or. h2+h3==0.0) then |
| 259 | 0 | if (use_2018_answers) then |
| 260 | 0 | h_min = hMinFrac*max( h_neglect, h0+h1+h2+h3 ) |
| 261 | else | |
| 262 | 0 | h_min = hMinFrac*max( h_neglect, (h0+h1)+(h2+h3) ) |
| 263 | endif | |
| 264 | 0 | h0 = max( h_min, h(i-2) ) |
| 265 | 0 | h1 = max( h_min, h(i-1) ) |
| 266 | 0 | h2 = max( h_min, h(i) ) |
| 267 | 0 | h3 = max( h_min, h(i+1) ) |
| 268 | endif | |
| 269 | ||
| 270 | 0 | if (use_2018_answers) then |
| 271 | 0 | f1 = (h0+h1) * (h2+h3) / (h1+h2) |
| 272 | 0 | f2 = h2 * u(i-1) + h1 * u(i) |
| 273 | 0 | f3 = 1.0 / (h0+h1+h2) + 1.0 / (h1+h2+h3) |
| 274 | 0 | et1 = f1 * f2 * f3 |
| 275 | et2 = ( h2 * (h2+h3) / ( (h0+h1+h2)*(h0+h1) ) ) * & | |
| 276 | 0 | ((h0+2.0*h1) * u(i-1) - h1 * u(i-2)) |
| 277 | et3 = ( h1 * (h0+h1) / ( (h1+h2+h3)*(h2+h3) ) ) * & | |
| 278 | 0 | ((2.0*h2+h3) * u(i) - h2 * u(i+1)) |
| 279 | 0 | edge_val(i,1) = (et1 + et2 + et3) / ( h0 + h1 + h2 + h3) |
| 280 | else | |
| 281 | 0 | I_h12 = 1.0 / (h1+h2) |
| 282 | 0 | I_den_et2 = 1.0 / ( ((h0+h1)+h2)*(h0+h1) ) ; I_h012 = (h0+h1) * I_den_et2 |
| 283 | 0 | I_den_et3 = 1.0 / ( (h1+(h2+h3))*(h2+h3) ) ; I_h123 = (h2+h3) * I_den_et3 |
| 284 | ||
| 285 | et1 = ( 1.0 + (h1 * I_h012 + (h0+h1) * I_h123) ) * I_h12 * (h2*(h2+h3)) * u(i-1) + & | |
| 286 | 0 | ( 1.0 + (h2 * I_h123 + (h2+h3) * I_h012) ) * I_h12 * (h1*(h0+h1)) * u(i) |
| 287 | 0 | et2 = ( h1 * (h2*(h2+h3)) * I_den_et2 ) * (u(i-1)-u(i-2)) |
| 288 | 0 | et3 = ( h2 * (h1*(h0+h1)) * I_den_et3 ) * (u(i) - u(i+1)) |
| 289 | 0 | edge_val(i,1) = (et1 + (et2 + et3)) / ((h0 + h1) + (h2 + h3)) |
| 290 | endif | |
| 291 | 0 | edge_val(i-1,2) = edge_val(i,1) |
| 292 | ||
| 293 | enddo ! end loop on interior cells | |
| 294 | ||
| 295 | ! Determine first two edge values | |
| 296 | 0 | if (use_2018_answers) then |
| 297 | 0 | h_min = max( h_neglect, hMinFrac*sum(h(1:4)) ) |
| 298 | 0 | x(1) = 0.0 |
| 299 | 0 | do i = 1,4 |
| 300 | 0 | dx = max(h_min, h(i) ) |
| 301 | 0 | x(i+1) = x(i) + dx |
| 302 | 0 | do j = 1,4 ; A(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j) ; enddo |
| 303 | 0 | B(i) = u(i) * dx |
| 304 | enddo | |
| 305 | ||
| 306 | 0 | call solve_linear_system( A, B, C, 4 ) |
| 307 | ||
| 308 | ! Set the edge values of the first cell | |
| 309 | 0 | edge_val(1,1) = evaluation_polynomial( C, 4, x(1) ) |
| 310 | 0 | edge_val(1,2) = evaluation_polynomial( C, 4, x(2) ) |
| 311 | else ! Use expressions with less sensitivity to roundoff | |
| 312 | 0 | do i=1,4 ; dz(i) = max(h_neglect, h(i) ) ; u_tmp(i) = u(i) ; enddo |
| 313 | 0 | call end_value_h4(dz, u_tmp, C) |
| 314 | ||
| 315 | ! Set the edge values of the first cell | |
| 316 | 0 | edge_val(1,1) = C(1) |
| 317 | 0 | edge_val(1,2) = C(1) + dz(1)*(C(2) + dz(1)*(C(3) + dz(1)*C(4))) |
| 318 | endif | |
| 319 | 0 | edge_val(2,1) = edge_val(1,2) |
| 320 | ||
| 321 | ! Determine two edge values of the last cell | |
| 322 | 0 | if (use_2018_answers) then |
| 323 | 0 | h_min = max( h_neglect, hMinFrac*sum(h(N-3:N)) ) |
| 324 | ||
| 325 | 0 | x(1) = 0.0 |
| 326 | 0 | do i = 1,4 |
| 327 | 0 | dx = max(h_min, h(N-4+i) ) |
| 328 | 0 | x(i+1) = x(i) + dx |
| 329 | 0 | do j = 1,4 ; A(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j) ; enddo |
| 330 | 0 | B(i) = u(N-4+i) * dx |
| 331 | enddo | |
| 332 | ||
| 333 | 0 | call solve_linear_system( A, B, C, 4 ) |
| 334 | ||
| 335 | ! Set the last and second to last edge values | |
| 336 | 0 | edge_val(N,2) = evaluation_polynomial( C, 4, x(5) ) |
| 337 | 0 | edge_val(N,1) = evaluation_polynomial( C, 4, x(4) ) |
| 338 | else | |
| 339 | ! Use expressions with less sensitivity to roundoff, including using a coordinate | |
| 340 | ! system that sets the origin at the last interface in the domain. | |
| 341 | 0 | do i=1,4 ; dz(i) = max(h_neglect, h(N+1-i) ) ; u_tmp(i) = u(N+1-i) ; enddo |
| 342 | 0 | call end_value_h4(dz, u_tmp, C) |
| 343 | ||
| 344 | ! Set the last and second to last edge values | |
| 345 | 0 | edge_val(N,2) = C(1) |
| 346 | 0 | edge_val(N,1) = C(1) + dz(1)*(C(2) + dz(1)*(C(3) + dz(1)*C(4))) |
| 347 | endif | |
| 348 | 0 | edge_val(N-1,2) = edge_val(N,1) |
| 349 | ||
| 350 | 0 | end subroutine edge_values_explicit_h4 |
| 351 | ||
| 352 | !> Compute h4 edge values (explicit fourth order accurate) | |
| 353 | !! in the same units as u. | |
| 354 | !! | |
| 355 | !! From (Colella & Woodward, JCP, 1984) and based on hybgen_ppm_coefs. | |
| 356 | !! | |
| 357 | !! Compute edge values based on fourth-order explicit estimates. | |
| 358 | !! These estimates are based on a cubic interpolant spanning four cells | |
| 359 | !! and evaluated at the location of the middle edge. An interpolant spanning | |
| 360 | !! cells i-2, i-1, i and i+1 is evaluated at edge i-1/2. The estimate for | |
| 361 | !! each edge is unique. | |
| 362 | !! | |
| 363 | !! i-2 i-1 i i+1 | |
| 364 | !! ..--o------o------o------o------o--.. | |
| 365 | !! i-1/2 | |
| 366 | !! | |
| 367 | !! For this fourth-order scheme, at least four cells must exist. | |
| 368 | 0 | subroutine edge_values_explicit_h4cw( N, h, u, edge_val, h_neglect ) |
| 369 | integer, intent(in) :: N !< Number of cells | |
| 370 | real, dimension(N), intent(in) :: h !< cell widths [H] | |
| 371 | real, dimension(N), intent(in) :: u !< cell average properties in arbitrary units [A] | |
| 372 | real, dimension(N,2), intent(inout) :: edge_val !< Returned edge values [A]; the second index | |
| 373 | !! is for the two edges of each cell. | |
| 374 | real, intent(in) :: h_neglect !< A negligibly small width [H] | |
| 375 | ||
| 376 | ! Local variables | |
| 377 | 0 | real :: dp(N) ! Input grid layer thicknesses, but with a minimum thickness [H ~> m or kg m-2] |
| 378 | real :: da ! Difference between the unlimited scalar edge value estimates [A] | |
| 379 | real :: a6 ! Scalar field differences that are proportional to the curvature [A] | |
| 380 | real :: slk, srk ! Differences between adjacent cell averages of scalars [A] | |
| 381 | real :: sck ! Scalar differences across a cell [A] | |
| 382 | 0 | real :: au(N) ! Scalar field difference across each cell [A] |
| 383 | 0 | real :: al(N), ar(N) ! Scalar field at the left and right edges of a cell [A] |
| 384 | 0 | real :: h112(N+1), h122(N+1) ! Combinations of thicknesses [H ~> m or kg m-2] |
| 385 | 0 | real :: I_h12(N+1) ! Inverses of combinations of thickesses [H-1 ~> m-1 or m2 kg-1] |
| 386 | 0 | real :: h2_h123(N) ! A ratio of a layer thickness of the sum of 3 adjacent thicknesses [nondim] |
| 387 | 0 | real :: I_h0123(N) ! Inverse of the sum of 4 adjacent thicknesses [H-1 ~> m-1 or m2 kg-1] |
| 388 | 0 | real :: h01_h112(N+1) ! A ratio of sums of adjacent thicknesses [nondim], 2/3 in the limit of uniform thicknesses. |
| 389 | 0 | real :: h23_h122(N+1) ! A ratio of sums of adjacent thicknesses [nondim], 2/3 in the limit of uniform thicknesses. |
| 390 | integer :: k | |
| 391 | ||
| 392 | ! Set the thicknesses for very thin layers to some minimum value. | |
| 393 | 0 | do k=1,N ; dp(k) = max(h(k), h_neglect) ; enddo |
| 394 | ||
| 395 | !compute grid metrics | |
| 396 | 0 | do k=2,N |
| 397 | 0 | h112(K) = 2.*dp(k-1) + dp(k) |
| 398 | 0 | h122(K) = dp(k-1) + 2.*dp(k) |
| 399 | 0 | I_h12(K) = 1.0 / (dp(k-1) + dp(k)) |
| 400 | enddo !k | |
| 401 | 0 | do k=2,N-1 |
| 402 | 0 | h2_h123(k) = dp(k) / (dp(k) + (dp(k-1)+dp(k+1))) |
| 403 | enddo | |
| 404 | 0 | do K=3,N-1 |
| 405 | 0 | I_h0123(K) = 1.0 / ((dp(k-2) + dp(k-1)) + (dp(k) + dp(k+1))) |
| 406 | ||
| 407 | 0 | h01_h112(K) = (dp(k-2) + dp(k-1)) / (2.0*dp(k-1) + dp(k)) |
| 408 | 0 | h23_h122(K) = (dp(k) + dp(k+1)) / (dp(k-1) + 2.0*dp(k)) |
| 409 | enddo | |
| 410 | ||
| 411 | !Compute average slopes: Colella, Eq. (1.8) | |
| 412 | 0 | au(1) = 0. |
| 413 | 0 | do k=2,N-1 |
| 414 | 0 | slk = u(k )-u(k-1) |
| 415 | 0 | srk = u(k+1)-u(k) |
| 416 | 0 | if (slk*srk > 0.) then |
| 417 | 0 | sck = h2_h123(k)*( h112(K)*srk*I_h12(K+1) + h122(K+1)*slk*I_h12(K) ) |
| 418 | 0 | au(k) = sign(min(abs(2.0*slk), abs(sck), abs(2.0*srk)), sck) |
| 419 | else | |
| 420 | 0 | au(k) = 0. |
| 421 | endif | |
| 422 | enddo !k | |
| 423 | 0 | au(N) = 0. |
| 424 | ||
| 425 | !Compute "first guess" edge values: Colella, Eq. (1.6) | |
| 426 | 0 | al(1) = u(1) ! 1st layer PCM |
| 427 | 0 | ar(1) = u(1) ! 1st layer PCM |
| 428 | 0 | al(2) = u(1) ! 1st layer PCM |
| 429 | 0 | do K=3,N-1 |
| 430 | ! This is a 4th order explicit edge value estimate. | |
| 431 | al(k) = (dp(k)*u(k-1) + dp(k-1)*u(k)) * I_h12(K) & | |
| 432 | + I_h0123(K)*( 2.*dp(k)*dp(k-1)*I_h12(K)*(u(k)-u(k-1)) * & | |
| 433 | ( h01_h112(K) - h23_h122(K) ) & | |
| 434 | 0 | + (dp(k)*au(k-1)*h23_h122(K) - dp(k-1)*au(k)*h01_h112(K)) ) |
| 435 | 0 | ar(k-1) = al(k) |
| 436 | enddo !k | |
| 437 | 0 | ar(N-1) = u(N) ! last layer PCM |
| 438 | 0 | al(N) = u(N) ! last layer PCM |
| 439 | 0 | ar(N) = u(N) ! last layer PCM |
| 440 | ||
| 441 | !Set coefficients | |
| 442 | 0 | do k=1,N |
| 443 | 0 | edge_val(k,1) = al(k) |
| 444 | 0 | edge_val(k,2) = ar(k) |
| 445 | enddo !k | |
| 446 | ||
| 447 | 0 | end subroutine edge_values_explicit_h4cw |
| 448 | ||
| 449 | !> Compute ih4 edge values (implicit fourth order accurate) | |
| 450 | !! in the same units as u. | |
| 451 | !! | |
| 452 | !! Compute edge values based on fourth-order implicit estimates. | |
| 453 | !! | |
| 454 | !! Fourth-order implicit estimates of edge values are based on a two-cell | |
| 455 | !! stencil. A tridiagonal system is set up and is based on expressing the | |
| 456 | !! edge values in terms of neighboring cell averages. The generic | |
| 457 | !! relationship is | |
| 458 | !! | |
| 459 | !! \f[ | |
| 460 | !! \alpha u_{i-1/2} + u_{i+1/2} + \beta u_{i+3/2} = a \bar{u}_i + b \bar{u}_{i+1} | |
| 461 | !! \f] | |
| 462 | !! | |
| 463 | !! and the stencil looks like this | |
| 464 | !! | |
| 465 | !! i i+1 | |
| 466 | !! ..--o------o------o--.. | |
| 467 | !! i-1/2 i+1/2 i+3/2 | |
| 468 | !! | |
| 469 | !! In this routine, the coefficients \f$\alpha\f$, \f$\beta\f$, \f$a\f$ and \f$b\f$ are | |
| 470 | !! computed, the tridiagonal system is built, boundary conditions are prescribed and | |
| 471 | !! the system is solved to yield edge-value estimates. | |
| 472 | !! | |
| 473 | !! There are N+1 unknowns and we are able to write N-1 equations. The | |
| 474 | !! boundary conditions close the system. | |
| 475 | 303504 | subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect, answer_date ) |
| 476 | integer, intent(in) :: N !< Number of cells | |
| 477 | real, dimension(N), intent(in) :: h !< cell widths [H] | |
| 478 | real, dimension(N), intent(in) :: u !< cell average properties in arbitrary units [A] | |
| 479 | real, dimension(N,2), intent(inout) :: edge_val !< Returned edge values [A]; the second index | |
| 480 | !! is for the two edges of each cell. | |
| 481 | real, intent(in) :: h_neglect !< A negligibly small width [H] | |
| 482 | integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use | |
| 483 | ||
| 484 | ! Local variables | |
| 485 | integer :: i, j ! loop indexes | |
| 486 | real :: h0, h1 ! cell widths [H] | |
| 487 | real :: h_min ! A minimal cell width [H] | |
| 488 | real :: h0_2, h1_2, h0h1 ! Squares or products of thicknesses [H2] | |
| 489 | real :: h0ph1_2 ! The square of a sum of thicknesses [H2] | |
| 490 | real :: h0ph1_4 ! The fourth power of a sum of thicknesses [H4] | |
| 491 | real :: alpha, beta ! stencil coefficients [nondim] | |
| 492 | real :: I_h2, abmix ! stencil coefficients [nondim] | |
| 493 | real :: a, b ! Combinations of stencil coefficients [nondim] | |
| 494 | real, dimension(5) :: x ! Coordinate system with 0 at edges [H] | |
| 495 | real, parameter :: C1_12 = 1.0 / 12.0 ! A rational constant [nondim] | |
| 496 | real, parameter :: C1_3 = 1.0 / 3.0 ! A rational constant [nondim] | |
| 497 | real, dimension(4) :: dz ! A temporary array of limited layer thicknesses [H] | |
| 498 | real, dimension(4) :: u_tmp ! A temporary array of cell average properties [A] | |
| 499 | real :: dx ! Differences and averages of successive values of x [H] | |
| 500 | real, dimension(4,4) :: Asys ! Differences in successive positions raised to various powers, | |
| 501 | ! in units that vary with the second (j) index as [H^j] | |
| 502 | real, dimension(4) :: Bsys ! The right hand side of the system to solve for C [A H] | |
| 503 | real, dimension(4) :: Csys ! The coefficients of a fit polynomial in units that vary | |
| 504 | ! with the index (j) as [A H^(j-1)] | |
| 505 | 607008 | real, dimension(N+1) :: tri_l, & ! tridiagonal system (lower diagonal) [nondim] |
| 506 | 607008 | tri_d, & ! tridiagonal system (middle diagonal) [nondim] |
| 507 | 607008 | tri_c, & ! tridiagonal system central value [nondim], with tri_d = tri_c+tri_l+tri_u |
| 508 | 607008 | tri_u, & ! tridiagonal system (upper diagonal) [nondim] |
| 509 | 607008 | tri_b, & ! tridiagonal system (right hand side) [A] |
| 510 | 607008 | tri_x ! tridiagonal system (solution vector) [A] |
| 511 | logical :: use_2018_answers ! If true use older, less accurate expressions. | |
| 512 | ||
| 513 | 303504 | use_2018_answers = .true. ; if (present(answer_date)) use_2018_answers = (answer_date < 20190101) |
| 514 | ||
| 515 | ! Loop on cells (except last one) | |
| 516 | 22762800 | do i = 1,N-1 |
| 517 | 22459296 | if (use_2018_answers) then |
| 518 | ! Get cell widths | |
| 519 | 0 | h0 = h(i) |
| 520 | 0 | h1 = h(i+1) |
| 521 | ! Avoid singularities when h0+h1=0 | |
| 522 | 0 | if (h0+h1==0.) then |
| 523 | 0 | h0 = h_neglect |
| 524 | 0 | h1 = h_neglect |
| 525 | endif | |
| 526 | ||
| 527 | ! Auxiliary calculations | |
| 528 | 0 | h0ph1_2 = (h0 + h1)**2 |
| 529 | 0 | h0ph1_4 = h0ph1_2**2 |
| 530 | 0 | h0h1 = h0 * h1 |
| 531 | 0 | h0_2 = h0 * h0 |
| 532 | 0 | h1_2 = h1 * h1 |
| 533 | ||
| 534 | ! Coefficients | |
| 535 | 0 | alpha = h1_2 / h0ph1_2 |
| 536 | 0 | beta = h0_2 / h0ph1_2 |
| 537 | 0 | a = 2.0 * h1_2 * ( h1_2 + 2.0 * h0_2 + 3.0 * h0h1 ) / h0ph1_4 |
| 538 | 0 | b = 2.0 * h0_2 * ( h0_2 + 2.0 * h1_2 + 3.0 * h0h1 ) / h0ph1_4 |
| 539 | ||
| 540 | 0 | tri_d(i+1) = 1.0 |
| 541 | else ! Use expressions with less sensitivity to roundoff | |
| 542 | ! Get cell widths | |
| 543 | 22459296 | h0 = max(h(i), h_neglect) |
| 544 | 22459296 | h1 = max(h(i+1), h_neglect) |
| 545 | ! The 1e-12 here attempts to balance truncation errors from the differences of | |
| 546 | ! large numbers against errors from approximating thin layers as non-vanishing. | |
| 547 | 22459296 | if (abs(h0) < 1.0e-12*abs(h1)) h0 = 1.0e-12*h1 |
| 548 | 22459296 | if (abs(h1) < 1.0e-12*abs(h0)) h1 = 1.0e-12*h0 |
| 549 | 22459296 | I_h2 = 1.0 / ((h0 + h1)**2) |
| 550 | 22459296 | alpha = (h1 * h1) * I_h2 |
| 551 | 22459296 | beta = (h0 * h0) * I_h2 |
| 552 | 22459296 | abmix = (h0 * h1) * I_h2 |
| 553 | 22459296 | a = 2.0 * alpha * ( alpha + 2.0 * beta + 3.0 * abmix ) |
| 554 | 22459296 | b = 2.0 * beta * ( beta + 2.0 * alpha + 3.0 * abmix ) |
| 555 | ||
| 556 | 22459296 | tri_c(i+1) = 2.0*abmix ! = 1.0 - alpha - beta |
| 557 | endif | |
| 558 | ||
| 559 | 22459296 | tri_l(i+1) = alpha |
| 560 | 22459296 | tri_u(i+1) = beta |
| 561 | ||
| 562 | 22762800 | tri_b(i+1) = a * u(i) + b * u(i+1) |
| 563 | ||
| 564 | enddo ! end loop on cells | |
| 565 | ||
| 566 | ! Boundary conditions: set the first boundary value | |
| 567 | 303504 | if (use_2018_answers) then |
| 568 | 0 | h_min = max( h_neglect, hMinFrac*sum(h(1:4)) ) |
| 569 | 0 | x(1) = 0.0 |
| 570 | 0 | do i = 1,4 |
| 571 | 0 | dx = max(h_min, h(i) ) |
| 572 | 0 | x(i+1) = x(i) + dx |
| 573 | 0 | do j = 1,4 ; Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo |
| 574 | 0 | Bsys(i) = u(i) * dx |
| 575 | enddo | |
| 576 | ||
| 577 | 0 | call solve_linear_system( Asys, Bsys, Csys, 4 ) |
| 578 | ||
| 579 | 0 | tri_b(1) = evaluation_polynomial( Csys, 4, x(1) ) ! Set the first edge value |
| 580 | 0 | tri_d(1) = 1.0 |
| 581 | else ! Use expressions with less sensitivity to roundoff | |
| 582 | 1517520 | do i=1,4 ; dz(i) = max(h_neglect, h(i) ) ; u_tmp(i) = u(i) ; enddo |
| 583 | 303504 | call end_value_h4(dz, u_tmp, Csys) |
| 584 | ||
| 585 | 303504 | tri_b(1) = Csys(1) ! Set the first edge value. |
| 586 | 303504 | tri_c(1) = 1.0 |
| 587 | endif | |
| 588 | 303504 | tri_u(1) = 0.0 ! tri_l(1) = 0.0 |
| 589 | ||
| 590 | ! Boundary conditions: set the last boundary value | |
| 591 | 303504 | if (use_2018_answers) then |
| 592 | 0 | h_min = max( h_neglect, hMinFrac*sum(h(N-3:N)) ) |
| 593 | 0 | x(1) = 0.0 |
| 594 | 0 | do i=1,4 |
| 595 | 0 | dx = max(h_min, h(N-4+i) ) |
| 596 | 0 | x(i+1) = x(i) + dx |
| 597 | 0 | do j = 1,4 ; Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo |
| 598 | 0 | Bsys(i) = u(N-4+i) * dx |
| 599 | enddo | |
| 600 | ||
| 601 | 0 | call solve_linear_system( Asys, Bsys, Csys, 4 ) |
| 602 | ||
| 603 | ! Set the last edge value | |
| 604 | 0 | tri_b(N+1) = evaluation_polynomial( Csys, 4, x(5) ) |
| 605 | 0 | tri_d(N+1) = 1.0 |
| 606 | ||
| 607 | else | |
| 608 | ! Use expressions with less sensitivity to roundoff, including using a coordinate | |
| 609 | ! system that sets the origin at the last interface in the domain. | |
| 610 | 1517520 | do i=1,4 ; dz(i) = max(h_neglect, h(N+1-i) ) ; u_tmp(i) = u(N+1-i) ; enddo |
| 611 | 303504 | call end_value_h4(dz, u_tmp, Csys) |
| 612 | ||
| 613 | 303504 | tri_b(N+1) = Csys(1) ! Set the last edge value |
| 614 | 303504 | tri_c(N+1) = 1.0 |
| 615 | endif | |
| 616 | 303504 | tri_l(N+1) = 0.0 ! tri_u(N+1) = 0.0 |
| 617 | ||
| 618 | ! Solve tridiagonal system and assign edge values | |
| 619 | 303504 | if (use_2018_answers) then |
| 620 | 0 | call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, N+1 ) |
| 621 | else | |
| 622 | 303504 | call solve_diag_dominant_tridiag( tri_l, tri_c, tri_u, tri_b, tri_x, N+1 ) |
| 623 | endif | |
| 624 | ||
| 625 | 303504 | edge_val(1,1) = tri_x(1) |
| 626 | 22762800 | do i=2,N |
| 627 | 22459296 | edge_val(i,1) = tri_x(i) |
| 628 | 22762800 | edge_val(i-1,2) = tri_x(i) |
| 629 | enddo | |
| 630 | 303504 | edge_val(N,2) = tri_x(N+1) |
| 631 | ||
| 632 | 303504 | end subroutine edge_values_implicit_h4 |
| 633 | ||
| 634 | !> Determine a one-sided 4th order polynomial fit of u to the data points for the purposes of specifying | |
| 635 | !! edge values, as described in the appendix of White and Adcroft JCP 2008. | |
| 636 | 607008 | subroutine end_value_h4(dz, u, Csys) |
| 637 | real, dimension(4), intent(in) :: dz !< The thicknesses of 4 layers, starting at the edge [H]. | |
| 638 | !! The values of dz must be positive. | |
| 639 | real, dimension(4), intent(in) :: u !< The average properties of 4 layers, starting at the edge [A] | |
| 640 | real, dimension(4), intent(out) :: Csys !< The four coefficients of a 4th order polynomial fit | |
| 641 | !! of u as a function of z [A H-(n-1)] | |
| 642 | ||
| 643 | ! Local variables | |
| 644 | real :: Wt(3,4) ! The weights of successive u differences in the 4 closed form expressions. | |
| 645 | ! The units of Wt vary with the second index as [H-(n-1)]. | |
| 646 | real :: h1, h2, h3, h4 ! Copies of the layer thicknesses [H] | |
| 647 | real :: h12, h23, h34 ! Sums of two successive thicknesses [H] | |
| 648 | real :: h123, h234 ! Sums of three successive thicknesses [H] | |
| 649 | real :: h1234 ! Sums of all four thicknesses [H] | |
| 650 | ! real :: I_h1 ! The inverse of the a thickness [H-1] | |
| 651 | real :: I_h12, I_h23, I_h34 ! The inverses of sums of two thicknesses [H-1] | |
| 652 | real :: I_h123, I_h234 ! The inverse of the sum of three thicknesses [H-1] | |
| 653 | real :: I_h1234 ! The inverse of the sum of all four thicknesses [H-1] | |
| 654 | real :: I_denom ! The inverse of the denominator some expressions [H-3] | |
| 655 | real :: I_denB3 ! The inverse of the product of three sums of thicknesses [H-3] | |
| 656 | real :: min_frac = 1.0e-6 ! The square of min_frac should be much larger than roundoff [nondim] | |
| 657 | real, parameter :: C1_3 = 1.0 / 3.0 ! A rational parameter [nondim] | |
| 658 | ||
| 659 | ! These are only used for code verification | |
| 660 | ! real, dimension(4) :: Atest ! The coefficients of an expression that is being tested. | |
| 661 | ! real :: zavg, u_mag, c_mag | |
| 662 | ! character(len=128) :: mesg | |
| 663 | ! real, parameter :: C1_12 = 1.0 / 12.0 | |
| 664 | ||
| 665 | ! if ((dz(1) == dz(2)) .and. (dz(1) == dz(3)) .and. (dz(1) == dz(4))) then | |
| 666 | ! ! There are simple closed-form expressions in this case | |
| 667 | ! I_h1 = 0.0 ; if (dz(1) > 0.0) I_h1 = 1.0 / dz(1) | |
| 668 | ! Csys(1) = u(1) + (-13.0 * (u(2)-u(1)) + 10.0 * (u(3)-u(2)) - 3.0 * (u(4)-u(3))) * (0.25*C1_3) | |
| 669 | ! Csys(2) = (35.0 * (u(2)-u(1)) - 34.0 * (u(3)-u(2)) + 11.0 * (u(4)-u(3))) * (0.25*C1_3 * I_h1) | |
| 670 | ! Csys(3) = (-5.0 * (u(2)-u(1)) + 8.0 * (u(3)-u(2)) - 3.0 * (u(4)-u(3))) * (0.25 * I_h1**2) | |
| 671 | ! Csys(4) = ((u(2)-u(1)) - 2.0 * (u(3)-u(2)) + (u(4)-u(3))) * (0.5*C1_3) | |
| 672 | ! else | |
| 673 | ||
| 674 | ! Express the coefficients as sums of the differences between properties of successive layers. | |
| 675 | ||
| 676 | 607008 | h1 = dz(1) ; h2 = dz(2) ; h3 = dz(3) ; h4 = dz(4) |
| 677 | ! Some of the weights used below are proportional to (h1/(h2+h3))**2 or (h1/(h2+h3))*(h2/(h3+h4)) | |
| 678 | ! so h2 and h3 should be adjusted to ensure that these ratios are not so large that property | |
| 679 | ! differences at the level of roundoff are amplified to be of order 1. | |
| 680 | 607008 | if ((h2+h3) < min_frac*h1) h3 = min_frac*h1 - h2 |
| 681 | 607008 | if ((h3+h4) < min_frac*h1) h4 = min_frac*h1 - h3 |
| 682 | ||
| 683 | 607008 | h12 = h1+h2 ; h23 = h2+h3 ; h34 = h3+h4 |
| 684 | 607008 | h123 = h12 + h3 ; h234 = h2 + h34 ; h1234 = h12 + h34 |
| 685 | ! Find 3 reciprocals with a single division for efficiency. | |
| 686 | 607008 | I_denB3 = 1.0 / (h123 * h12 * h23) |
| 687 | 607008 | I_h12 = (h123 * h23) * I_denB3 |
| 688 | 607008 | I_h23 = (h12 * h123) * I_denB3 |
| 689 | 607008 | I_h123 = (h12 * h23) * I_denB3 |
| 690 | 607008 | I_denom = 1.0 / ( h1234 * (h234 * h34) ) |
| 691 | 607008 | I_h34 = (h1234 * h234) * I_denom |
| 692 | 607008 | I_h234 = (h1234 * h34) * I_denom |
| 693 | 607008 | I_h1234 = (h234 * h34) * I_denom |
| 694 | ||
| 695 | ! Calculation coefficients in the four equations | |
| 696 | ||
| 697 | ! The expressions for Csys(3) and Csys(4) come from reducing the 4x4 matrix problem into the following 2x2 | |
| 698 | ! matrix problem, then manipulating the analytic solution to avoid any subtraction and simplifying. | |
| 699 | ! (C1_3 * h123 * h23) * Csys(3) + (0.25 * h123 * h23 * (h3 + 2.0*h2 + 3.0*h1)) * Csys(4) = | |
| 700 | ! (u(3)-u(1)) - (u(2)-u(1)) * (h12 + h23) * I_h12 | |
| 701 | ! (C1_3 * ((h23 + h34) * h1234 + h23 * h3)) * Csys(3) + | |
| 702 | ! (0.25 * ((h1234 + h123 + h12 + h1) * h23 * h3 + (h1234 + h12 + h1) * (h23 + h34) * h1234)) * Csys(4) = | |
| 703 | ! (u(4)-u(1)) - (u(2)-u(1)) * (h123 + h234) * I_h12 | |
| 704 | ! The final expressions for Csys(1) and Csys(2) were derived by algebraically manipulating the following expressions: | |
| 705 | ! Csys(1) = (C1_3 * h1 * h12 * Csys(3) + 0.25 * h1 * h12 * (2.0*h1+h2) * Csys(4)) - & | |
| 706 | ! (h1*I_h12)*(u(2)-u(1)) + u(1) | |
| 707 | ! Csys(2) = (-2.0*C1_3 * (2.0*h1+h2) * Csys(3) - 0.5 * (h1**2 + h12 * (2.0*h1+h2)) * Csys(4)) + & | |
| 708 | ! 2.0*I_h12 * (u(2)-u(1)) | |
| 709 | ! These expressions are typically evaluated at x=0 and x=h1, so it is important that these are well behaved | |
| 710 | ! for these values, suggesting that h1/h23 and h1/h34 should not be allowed to be too large. | |
| 711 | ||
| 712 | 607008 | Wt(1,1) = -h1 * (I_h1234 + I_h123 + I_h12) ! > -3 |
| 713 | 607008 | Wt(2,1) = h1 * h12 * ( I_h234 * I_h1234 + I_h23 * (I_h234 + I_h123) ) ! < (h1/h234) + (h1/h23)*(2+(h1/h234)) |
| 714 | 607008 | Wt(3,1) = -h1 * h12 * h123 * I_denom ! > -(h1/h34)*(1+(h1/h234)) |
| 715 | ||
| 716 | 607008 | Wt(1,2) = 2.0 * (I_h12*(1.0 + (h1+h12) * (I_h1234 + I_h123)) + h1 * I_h1234*I_h123) ! < 10/h12 |
| 717 | Wt(2,2) = -2.0 * ((h1 * h12 * I_h1234) * (I_h23 * (I_h234 + I_h123)) + & ! > -(10+6*(h1/h234))/h23 | |
| 718 | 607008 | (h1+h12) * ( I_h1234*I_h234 + I_h23 * (I_h234 + I_h123) ) ) |
| 719 | 607008 | Wt(3,2) = 2.0 * ((h1+h12) * h123 + h1*h12 ) * I_denom ! < (2+(6*h1/h234)) / h34 |
| 720 | ||
| 721 | 607008 | Wt(1,3) = -3.0 * I_h12 * I_h123* ( 1.0 + I_h1234 * ((h1+h12)+h123) ) ! > -12 / (h12*h123) |
| 722 | 607008 | Wt(2,3) = 3.0 * I_h23 * ( I_h123 + I_h1234 * ((h1+h12)+h123) * (I_h123 + I_h234) ) ! < 12 / (h23^2) |
| 723 | 607008 | Wt(3,3) = -3.0 * ((h1+h12)+h123) * I_denom ! > -9 / (h234*h23) |
| 724 | ||
| 725 | 607008 | Wt(1,4) = 4.0 * I_h1234 * I_h123 * I_h12 ! Wt*h1^3 < 4 |
| 726 | 607008 | Wt(2,4) = -4.0 * I_h1234 * (I_h23 * (I_h123 + I_h234)) ! Wt*h1^3 > -4* (h1/h23)*(1+h1/h234) |
| 727 | 607008 | Wt(3,4) = 4.0 * I_denom ! = 4.0*I_h1234 * I_h234 * I_h34 ! Wt*h1^3 < 4 * (h1/h234)*(h1/h34) |
| 728 | ||
| 729 | 607008 | Csys(1) = ((u(1) + (Wt(1,1) * (u(2)-u(1)))) + (Wt(2,1) * (u(3)-u(2)))) + (Wt(3,1) * (u(4)-u(3))) |
| 730 | 607008 | Csys(2) = ((Wt(1,2) * (u(2)-u(1))) + (Wt(2,2) * (u(3)-u(2)))) + (Wt(3,2) * (u(4)-u(3))) |
| 731 | 607008 | Csys(3) = ((Wt(1,3) * (u(2)-u(1))) + (Wt(2,3) * (u(3)-u(2)))) + (Wt(3,3) * (u(4)-u(3))) |
| 732 | 607008 | Csys(4) = ((Wt(1,4) * (u(2)-u(1))) + (Wt(2,4) * (u(3)-u(2)))) + (Wt(3,4) * (u(4)-u(3))) |
| 733 | ||
| 734 | ! endif ! End of non-uniform layer thickness branch. | |
| 735 | ||
| 736 | ! To verify that these answers are correct, uncomment the following: | |
| 737 | ! u_mag = 0.0 ; do i=1,4 ; u_mag = max(u_mag, abs(u(i))) ; enddo | |
| 738 | ! do i = 1,4 | |
| 739 | ! if (i==1) then ; zavg = 0.5*dz(i) ; else ; zavg = zavg + 0.5*(dz(i-1)+dz(i)) ; endif | |
| 740 | ! Atest(1) = 1.0 | |
| 741 | ! Atest(2) = zavg ! = ( (z(i+1)**2) - (z(i)**2) ) / (2*dz(i)) | |
| 742 | ! Atest(3) = (zavg**2 + 0.25*C1_3*dz(i)**2) ! = ( (z(i+1)**3) - (z(i)**3) ) / (3*dz(i)) | |
| 743 | ! Atest(4) = zavg * (zavg**2 + 0.25*dz(i)**2) ! = ( (z(i+1)**4) - (z(i)**4) ) / (4*dz(i)) | |
| 744 | ! c_mag = 1.0 ; do k=0,3 ; do j=1,3 ; c_mag = c_mag + abs(Wt(j,k+1) * zavg**k) ; enddo ; enddo | |
| 745 | ! write(mesg, '("end_value_h4 line ", i2, " c_mag = ", es10.2, " u_mag = ", es10.2)') i, c_mag, u_mag | |
| 746 | ! call test_line(mesg, 4, Atest, Csys, u(i), u_mag*c_mag, tol=1.0e-15) | |
| 747 | ! enddo | |
| 748 | ||
| 749 | 607008 | end subroutine end_value_h4 |
| 750 | ||
| 751 | ||
| 752 | !------------------------------------------------------------------------------ | |
| 753 | !> Compute ih3 edge slopes (implicit third order accurate) | |
| 754 | !! in the same units as h. | |
| 755 | !! | |
| 756 | !! Compute edge slopes based on third-order implicit estimates. Note that | |
| 757 | !! the estimates are fourth-order accurate on uniform grids | |
| 758 | !! | |
| 759 | !! Third-order implicit estimates of edge slopes are based on a two-cell | |
| 760 | !! stencil. A tridiagonal system is set up and is based on expressing the | |
| 761 | !! edge slopes in terms of neighboring cell averages. The generic | |
| 762 | !! relationship is | |
| 763 | !! | |
| 764 | !! \f[ | |
| 765 | !! \alpha u'_{i-1/2} + u'_{i+1/2} + \beta u'_{i+3/2} = | |
| 766 | !! a \bar{u}_i + b \bar{u}_{i+1} | |
| 767 | !! \f] | |
| 768 | !! | |
| 769 | !! and the stencil looks like this | |
| 770 | !! | |
| 771 | !! i i+1 | |
| 772 | !! ..--o------o------o--.. | |
| 773 | !! i-1/2 i+1/2 i+3/2 | |
| 774 | !! | |
| 775 | !! In this routine, the coefficients \f$\alpha\f$, \f$\beta\f$, a and b are computed, | |
| 776 | !! the tridiagonal system is built, boundary conditions are prescribed and | |
| 777 | !! the system is solved to yield edge-slope estimates. | |
| 778 | !! | |
| 779 | !! There are N+1 unknowns and we are able to write N-1 equations. The | |
| 780 | !! boundary conditions close the system. | |
| 781 | 0 | subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect, answer_date ) |
| 782 | integer, intent(in) :: N !< Number of cells | |
| 783 | real, dimension(N), intent(in) :: h !< cell widths [H] | |
| 784 | real, dimension(N), intent(in) :: u !< cell average properties in arbitrary units [A] | |
| 785 | real, dimension(N,2), intent(inout) :: edge_slopes !< Returned edge slopes [A H-1]; the | |
| 786 | !! second index is for the two edges of each cell. | |
| 787 | real, intent(in) :: h_neglect !< A negligibly small width [H] | |
| 788 | integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use | |
| 789 | ||
| 790 | ! Local variables | |
| 791 | integer :: i, j ! loop indexes | |
| 792 | real :: h0, h1 ! cell widths [H or nondim] | |
| 793 | real :: h0_2, h1_2, h0h1 ! products of cell widths [H2 or nondim] | |
| 794 | real :: h0_3, h1_3 ! products of three cell widths [H3 or nondim] | |
| 795 | real :: d ! A temporary variable [H3] | |
| 796 | real :: I_d ! A temporary variable [nondim] | |
| 797 | real :: I_h ! Inverses of thicknesses [H-1] | |
| 798 | real :: alpha, beta ! stencil coefficients [nondim] | |
| 799 | real :: a, b ! weights of cells [H-1] | |
| 800 | real, parameter :: C1_12 = 1.0 / 12.0 ! A rational parameter [nondim] | |
| 801 | real, dimension(4) :: dz ! A temporary array of limited layer thicknesses [H] | |
| 802 | real, dimension(4) :: u_tmp ! A temporary array of cell average properties [A] | |
| 803 | real, dimension(5) :: x ! Coordinate system with 0 at edges [H] | |
| 804 | real :: dx ! Differences and averages of successive values of x [H] | |
| 805 | real, dimension(4,4) :: Asys ! Differences in successive positions raised to various powers, | |
| 806 | ! in units that vary with the second (j) index as [H^j] | |
| 807 | real, dimension(4) :: Bsys ! The right hand side of the system to solve for C [A H] | |
| 808 | real, dimension(4) :: Csys ! The coefficients of a fit polynomial in units that vary with the | |
| 809 | ! index (j) as [A H^(j-1)] | |
| 810 | real, dimension(3) :: Dsys ! The coefficients of the first derivative of the fit polynomial | |
| 811 | ! in units that vary with the index (j) as [A H^(j-2)] | |
| 812 | 0 | real, dimension(N+1) :: tri_l, & ! tridiagonal system (lower diagonal) [nondim] |
| 813 | 0 | tri_d, & ! tridiagonal system (middle diagonal) [nondim] |
| 814 | 0 | tri_c, & ! tridiagonal system central value [nondim], with tri_d = tri_c+tri_l+tri_u |
| 815 | 0 | tri_u, & ! tridiagonal system (upper diagonal) [nondim] |
| 816 | 0 | tri_b, & ! tridiagonal system (right hand side) [A H-1] |
| 817 | 0 | tri_x ! tridiagonal system (solution vector) [A H-1] |
| 818 | real :: hNeglect3 ! hNeglect^3 [H3]. | |
| 819 | logical :: use_2018_answers ! If true use older, less accurate expressions. | |
| 820 | ||
| 821 | 0 | hNeglect3 = h_neglect**3 |
| 822 | 0 | use_2018_answers = .true. ; if (present(answer_date)) use_2018_answers = (answer_date < 20190101) |
| 823 | ||
| 824 | ! Loop on cells (except last one) | |
| 825 | 0 | do i = 1,N-1 |
| 826 | ||
| 827 | 0 | if (use_2018_answers) then |
| 828 | ! Get cell widths | |
| 829 | 0 | h0 = h(i) |
| 830 | 0 | h1 = h(i+1) |
| 831 | ||
| 832 | ! Auxiliary calculations | |
| 833 | 0 | h0h1 = h0 * h1 |
| 834 | 0 | h0_2 = h0 * h0 |
| 835 | 0 | h1_2 = h1 * h1 |
| 836 | 0 | h0_3 = h0_2 * h0 |
| 837 | 0 | h1_3 = h1_2 * h1 |
| 838 | ||
| 839 | 0 | d = 4.0 * h0h1 * ( h0 + h1 ) + h1_3 + h0_3 |
| 840 | ||
| 841 | ! Coefficients | |
| 842 | 0 | alpha = h1 * (h0_2 + h0h1 - h1_2) / ( d + hNeglect3 ) |
| 843 | 0 | beta = h0 * (h1_2 + h0h1 - h0_2) / ( d + hNeglect3 ) |
| 844 | 0 | a = -12.0 * h0h1 / ( d + hNeglect3 ) |
| 845 | 0 | b = -a |
| 846 | ||
| 847 | 0 | tri_l(i+1) = alpha |
| 848 | 0 | tri_d(i+1) = 1.0 |
| 849 | 0 | tri_u(i+1) = beta |
| 850 | ||
| 851 | 0 | tri_b(i+1) = a * u(i) + b * u(i+1) |
| 852 | else | |
| 853 | ! Get cell widths | |
| 854 | 0 | h0 = max(h(i), h_neglect) |
| 855 | 0 | h1 = max(h(i+1), h_neglect) |
| 856 | ||
| 857 | 0 | I_h = 1.0 / (h0 + h1) |
| 858 | 0 | h0 = h0 * I_h ; h1 = h1 * I_h |
| 859 | ||
| 860 | 0 | h0h1 = h0 * h1 ; h0_2 = h0 * h0 ; h1_2 = h1 * h1 |
| 861 | 0 | h0_3 = h0_2 * h0 ; h1_3 = h1_2 * h1 |
| 862 | ||
| 863 | ! Set the tridiagonal coefficients | |
| 864 | 0 | I_d = 1.0 / (4.0 * h0h1 * ( h0 + h1 ) + h1_3 + h0_3) ! = 1 / ((h0 + h1)**3 + h0*h1*(h0 + h1)) |
| 865 | 0 | tri_l(i+1) = (h1 * ((h0_2 + h0h1) - h1_2)) * I_d |
| 866 | ! tri_d(i+1) = 1.0 | |
| 867 | 0 | tri_c(i+1) = 2.0 * ((h0_2 + h1_2) * (h0 + h1)) * I_d |
| 868 | 0 | tri_u(i+1) = (h0 * ((h1_2 + h0h1) - h0_2)) * I_d |
| 869 | ! The following expressions have been simplified using the nondimensionalization above: | |
| 870 | ! I_d = 1.0 / (1.0 + h0h1) | |
| 871 | ! tri_l(i+1) = (h0h1 - h1_3) * I_d | |
| 872 | ! tri_c(i+1) = 2.0 * (h0_2 + h1_2) * I_d | |
| 873 | ! tri_u(i+1) = (h0h1 - h0_3) * I_d | |
| 874 | ||
| 875 | 0 | tri_b(i+1) = 12.0 * (h0h1 * I_d) * ((u(i+1) - u(i)) * I_h) |
| 876 | endif | |
| 877 | ||
| 878 | enddo ! end loop on cells | |
| 879 | ||
| 880 | ! Boundary conditions: set the first edge slope | |
| 881 | 0 | if (use_2018_answers) then |
| 882 | 0 | x(1) = 0.0 |
| 883 | 0 | do i = 1,4 |
| 884 | 0 | dx = h(i) |
| 885 | 0 | x(i+1) = x(i) + dx |
| 886 | 0 | do j = 1,4 ; Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo |
| 887 | 0 | Bsys(i) = u(i) * dx |
| 888 | enddo | |
| 889 | ||
| 890 | 0 | call solve_linear_system( Asys, Bsys, Csys, 4 ) |
| 891 | ||
| 892 | 0 | Dsys(1) = Csys(2) ; Dsys(2) = 2.0 * Csys(3) ; Dsys(3) = 3.0 * Csys(4) |
| 893 | 0 | tri_b(1) = evaluation_polynomial( Dsys, 3, x(1) ) ! Set the first edge slope |
| 894 | 0 | tri_d(1) = 1.0 |
| 895 | else ! Use expressions with less sensitivity to roundoff | |
| 896 | 0 | do i=1,4 ; dz(i) = max(h_neglect, h(i) ) ; u_tmp(i) = u(i) ; enddo |
| 897 | 0 | call end_value_h4(dz, u_tmp, Csys) |
| 898 | ||
| 899 | ! Set the first edge slope | |
| 900 | 0 | tri_b(1) = Csys(2) |
| 901 | 0 | tri_c(1) = 1.0 |
| 902 | endif | |
| 903 | 0 | tri_u(1) = 0.0 ! tri_l(1) = 0.0 |
| 904 | ||
| 905 | ! Boundary conditions: set the last edge slope | |
| 906 | 0 | if (use_2018_answers) then |
| 907 | 0 | x(1) = 0.0 |
| 908 | 0 | do i = 1,4 |
| 909 | 0 | dx = h(N-4+i) |
| 910 | 0 | x(i+1) = x(i) + dx |
| 911 | 0 | do j = 1,4 ; Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo |
| 912 | 0 | Bsys(i) = u(N-4+i) * dx |
| 913 | enddo | |
| 914 | ||
| 915 | 0 | call solve_linear_system( Asys, Bsys, Csys, 4 ) |
| 916 | ||
| 917 | 0 | Dsys(1) = Csys(2) ; Dsys(2) = 2.0 * Csys(3) ; Dsys(3) = 3.0 * Csys(4) |
| 918 | ! Set the last edge slope | |
| 919 | 0 | tri_b(N+1) = evaluation_polynomial( Dsys, 3, x(5) ) |
| 920 | 0 | tri_d(N+1) = 1.0 |
| 921 | else | |
| 922 | ! Use expressions with less sensitivity to roundoff, including using a coordinate | |
| 923 | ! system that sets the origin at the last interface in the domain. | |
| 924 | 0 | do i=1,4 ; dz(i) = max(h_neglect, h(N+1-i) ) ; u_tmp(i) = u(N+1-i) ; enddo |
| 925 | ||
| 926 | 0 | call end_value_h4(dz, u_tmp, Csys) |
| 927 | ||
| 928 | ! Set the last edge slope | |
| 929 | 0 | tri_b(N+1) = -Csys(2) |
| 930 | 0 | tri_c(N+1) = 1.0 |
| 931 | endif | |
| 932 | 0 | tri_l(N+1) = 0.0 ! tri_u(N+1) = 0.0 |
| 933 | ||
| 934 | ! Solve tridiagonal system and assign edge slopes | |
| 935 | 0 | if (use_2018_answers) then |
| 936 | 0 | call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, N+1 ) |
| 937 | else | |
| 938 | 0 | call solve_diag_dominant_tridiag( tri_l, tri_c, tri_u, tri_b, tri_x, N+1 ) |
| 939 | endif | |
| 940 | ||
| 941 | 0 | do i = 2,N |
| 942 | 0 | edge_slopes(i,1) = tri_x(i) |
| 943 | 0 | edge_slopes(i-1,2) = tri_x(i) |
| 944 | enddo | |
| 945 | 0 | edge_slopes(1,1) = tri_x(1) |
| 946 | 0 | edge_slopes(N,2) = tri_x(N+1) |
| 947 | ||
| 948 | 0 | end subroutine edge_slopes_implicit_h3 |
| 949 | ||
| 950 | ||
| 951 | !------------------------------------------------------------------------------ | |
| 952 | !> Compute ih5 edge slopes (implicit fifth order accurate) | |
| 953 | 0 | subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect, answer_date ) |
| 954 | integer, intent(in) :: N !< Number of cells | |
| 955 | real, dimension(N), intent(in) :: h !< cell widths [H] | |
| 956 | real, dimension(N), intent(in) :: u !< cell average properties in arbitrary units [A] | |
| 957 | real, dimension(N,2), intent(inout) :: edge_slopes !< Returned edge slopes [A H-1]; the | |
| 958 | !! second index is for the two edges of each cell. | |
| 959 | real, intent(in) :: h_neglect !< A negligibly small width [H] | |
| 960 | integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use | |
| 961 | ||
| 962 | ! ----------------------------------------------------------------------------- | |
| 963 | ! Fifth-order implicit estimates of edge slopes are based on a four-cell, | |
| 964 | ! three-edge stencil. A tridiagonal system is set up and is based on | |
| 965 | ! expressing the edge slopes in terms of neighboring cell averages. | |
| 966 | ! | |
| 967 | ! The generic relationship is | |
| 968 | ! | |
| 969 | ! \alpha u'_{i-1/2} + u'_{i+1/2} + \beta u'_{i+3/2} = | |
| 970 | ! a \bar{u}_{i-1} + b \bar{u}_i + c \bar{u}_{i+1} + d \bar{u}_{i+2} | |
| 971 | ! | |
| 972 | ! and the stencil looks like this | |
| 973 | ! | |
| 974 | ! i-1 i i+1 i+2 | |
| 975 | ! ..--o------o------o------o------o--.. | |
| 976 | ! i-1/2 i+1/2 i+3/2 | |
| 977 | ! | |
| 978 | ! In this routine, the coefficients \alpha, \beta, a, b, c and d are | |
| 979 | ! computed, the tridiagonal system is built, boundary conditions are | |
| 980 | ! prescribed and the system is solved to yield edge-value estimates. | |
| 981 | ! | |
| 982 | ! Note that the centered stencil only applies to edges 3 to N-1 (edges are | |
| 983 | ! numbered 1 to n+1), which yields N-3 equations for N+1 unknowns. Two other | |
| 984 | ! equations are written by using a right-biased stencil for edge 2 and a | |
| 985 | ! left-biased stencil for edge N. The prescription of boundary conditions | |
| 986 | ! (using sixth-order polynomials) closes the system. | |
| 987 | ! | |
| 988 | ! CAUTION: For each edge, in order to determine the coefficients of the | |
| 989 | ! implicit expression, a 6x6 linear system is solved. This may | |
| 990 | ! become computationally expensive if regridding is carried out | |
| 991 | ! often. Figuring out closed-form expressions for these coefficients | |
| 992 | ! on nonuniform meshes turned out to be intractable. | |
| 993 | ! ----------------------------------------------------------------------------- | |
| 994 | ||
| 995 | ! Local variables | |
| 996 | real :: h0, h1, h2, h3 ! cell widths [H] | |
| 997 | real :: hMin ! The minimum thickness used in these calculations [H] | |
| 998 | real :: h01, h01_2 ! Summed thicknesses to various powers [H^n ~> m^n or kg^n m-2n] | |
| 999 | real :: h23, h23_2 ! Summed thicknesses to various powers [H^n ~> m^n or kg^n m-2n] | |
| 1000 | real :: h1_2, h2_2 ! Squares of thicknesses [H2] | |
| 1001 | real :: h1_3, h2_3 ! Cubes of thicknesses [H3] | |
| 1002 | real :: h1_4, h2_4 ! Fourth powers of thicknesses [H4] | |
| 1003 | real :: h1_5, h2_5 ! Fifth powers of thicknesses [H5] | |
| 1004 | real :: alpha, beta ! stencil coefficients [nondim] | |
| 1005 | real, dimension(7) :: x ! Coordinate system with 0 at edges in the same units as h [H] | |
| 1006 | real, parameter :: C1_12 = 1.0 / 12.0 ! A rational parameter [nondim] | |
| 1007 | real, parameter :: C5_6 = 5.0 / 6.0 ! A rational parameter [nondim] | |
| 1008 | real :: dx ! Differences between successive values of x in the same units as h [H] | |
| 1009 | real :: xavg ! Average of successive values of x in the same units as h [H] | |
| 1010 | real, dimension(6,6) :: Asys ! The matrix that is being inverted for a solution, | |
| 1011 | ! in units that might vary with the second (j) index as [H^j] | |
| 1012 | real, dimension(6) :: Bsys ! The right hand side of the system to solve for C in various | |
| 1013 | ! units that sometimes vary with the intex (j) as [H^(j-1)] or [H^j] | |
| 1014 | ! or might be [A] | |
| 1015 | real, dimension(6) :: Csys ! The solution to a matrix equation usually [nondim] in this routine. | |
| 1016 | 0 | real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal) [nondim] |
| 1017 | 0 | tri_d, & ! trid. system (middle diagonal) [nondim] |
| 1018 | 0 | tri_u, & ! trid. system (upper diagonal) [nondim] |
| 1019 | 0 | tri_b, & ! trid. system (rhs) [A H-1] |
| 1020 | 0 | tri_x ! trid. system (unknowns vector) [A H-1] |
| 1021 | real :: h_Min_Frac = 1.0e-4 ! A minimum fractional thickness [nondim] | |
| 1022 | integer :: i, k ! loop indexes | |
| 1023 | ||
| 1024 | ! Loop on cells (except the first and last ones) | |
| 1025 | 0 | do k = 2,N-2 |
| 1026 | ! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes. | |
| 1027 | 0 | hMin = max(h_neglect, h_Min_Frac*((h(k-1) + h(k)) + (h(k+1) + h(k+2)))) |
| 1028 | 0 | h0 = max(h(k-1), hMin) ; h1 = max(h(k), hMin) |
| 1029 | 0 | h2 = max(h(k+1), hMin) ; h3 = max(h(k+2), hMin) |
| 1030 | ||
| 1031 | ! Auxiliary calculations | |
| 1032 | 0 | h1_2 = h1 * h1 ; h1_3 = h1_2 * h1 ; h1_4 = h1_2 * h1_2 ; h1_5 = h1_3 * h1_2 |
| 1033 | 0 | h2_2 = h2 * h2 ; h2_3 = h2_2 * h2 ; h2_4 = h2_2 * h2_2 ; h2_5 = h2_3 * h2_2 |
| 1034 | ||
| 1035 | ! Compute matrix entries as described in Eq. (52) of White and Adcroft (2009). The last 4 rows are | |
| 1036 | ! Asys(1:6,n) = (/ -n*(n-1)*(-h1)**(n-2), -n*(n-1)*h1**(n-2), (-1)**(n-1) * ((h0+h1)**n - h0**n) / h0, & | |
| 1037 | ! (-h1)**(n-1), h2**(n-1), ((h2+h3)**n - h2**n) / h3 /) | |
| 1038 | ||
| 1039 | 0 | Asys(1:6,1) = (/ 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 /) |
| 1040 | 0 | Asys(1:6,2) = (/ 2.0, 2.0, (2.0*h1 + h0), h1, -h2, -(2.0*h2 + h3) /) |
| 1041 | Asys(1:6,3) = (/ 6.0*h1, -6.0* h2, (3.0*h1_2 + h0*(3.0*h1 + h0)), & | |
| 1042 | 0 | h1_2, h2_2, (3.0*h2_2 + h3*(3.0*h2 + h3)) /) |
| 1043 | Asys(1:6,4) = (/ -12.0*h1_2, -12.0*h2_2, -(4.0*h1_3 + h0*(6.0*h1_2 + h0*(4.0*h1 + h0))), & | |
| 1044 | 0 | -h1_3, h2_3, (4.0*h2_3 + h3*(6.0*h2_2 + h3*(4.0*h2 + h3))) /) |
| 1045 | Asys(1:6,5) = (/ 20.0*h1_3, -20.0*h2_3, (5.0*h1_4 + h0*(10.0*h1_3 + h0*(10.0*h1_2 + h0*(5.0*h1 + h0)))), & | |
| 1046 | 0 | h1_4, h2_4, (5.0*h2_4 + h3*(10.0*h2_3 + h3*(10.0*h2_2 + h3*(5.0*h2 + h3)))) /) |
| 1047 | Asys(1:6,6) = (/ -30.0*h1_4, -30.0*h2_4, & | |
| 1048 | -(6.0*h1_5 + h0*(15.0*h1_4 + h0*(20.0*h1_3 + h0*(15.0*h1_2 + h0*(6.0*h1 + h0))))), & | |
| 1049 | -h1_5, h2_5, & | |
| 1050 | 0 | (6.0*h2_5 + h3*(15.0*h2_4 + h3*(20.0*h2_3 + h3*(15.0*h2_2 + h3*(6.0*h2 + h3))))) /) |
| 1051 | 0 | Bsys(1:6) = (/ 0.0, -2.0, 0.0, 0.0, 0.0, 0.0 /) |
| 1052 | ||
| 1053 | 0 | call linear_solver( 6, Asys, Bsys, Csys ) |
| 1054 | ||
| 1055 | 0 | alpha = Csys(1) |
| 1056 | 0 | beta = Csys(2) |
| 1057 | ||
| 1058 | 0 | tri_l(k+1) = alpha |
| 1059 | 0 | tri_d(k+1) = 1.0 |
| 1060 | 0 | tri_u(k+1) = beta |
| 1061 | 0 | tri_b(k+1) = Csys(3) * u(k-1) + Csys(4) * u(k) + Csys(5) * u(k+1) + Csys(6) * u(k+2) |
| 1062 | ||
| 1063 | enddo ! end loop on cells | |
| 1064 | ||
| 1065 | ! Use a right-biased stencil for the second row, as described in Eq. (53) of White and Adcroft (2009). | |
| 1066 | ||
| 1067 | ! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes. | |
| 1068 | 0 | hMin = max(h_neglect, h_Min_Frac*((h(1) + h(2)) + (h(3) + h(4)))) |
| 1069 | 0 | h0 = max(h(1), hMin) ; h1 = max(h(2), hMin) |
| 1070 | 0 | h2 = max(h(3), hMin) ; h3 = max(h(4), hMin) |
| 1071 | ||
| 1072 | ! Auxiliary calculations | |
| 1073 | 0 | h1_2 = h1 * h1 ; h1_3 = h1_2 * h1 ; h1_4 = h1_2 * h1_2 ; h1_5 = h1_3 * h1_2 |
| 1074 | 0 | h2_2 = h2 * h2 ; h2_3 = h2_2 * h2 ; h2_4 = h2_2 * h2_2 ; h2_5 = h2_3 * h2_2 |
| 1075 | 0 | h01 = h0 + h1 ; h01_2 = h01 * h01 |
| 1076 | ||
| 1077 | ! Compute matrix entries | |
| 1078 | 0 | Asys(1:6,1) = (/ 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 /) |
| 1079 | 0 | Asys(1:6,2) = (/ 2.0, 2.0, (2.0*h1 + h0), h1, -h2, -(2.0*h2 + h3) /) |
| 1080 | Asys(1:6,3) = (/ 6.0*h01, 0.0, (3.0*h1_2 + h0*(3.0*h1 + h0)), & | |
| 1081 | 0 | h1_2, h2_2, (3.0*h2_2 + h3*(3.0*h2 + h3)) /) |
| 1082 | Asys(1:6,4) = (/ -12.0*h01_2, 0.0, -(4.0*h1_3 + h0*(6.0*h1_2 + h0*(4.0*h1 + h0))), & | |
| 1083 | 0 | -h1_3, h2_3, (4.0*h2_3 + h3*(6.0*h2_2 + h3*(4.0*h2 + h3))) /) |
| 1084 | Asys(1:6,5) = (/ 20.0*(h01*h01_2), 0.0, (5.0*h1_4 + h0*(10.0*h1_3 + h0*(10.0*h1_2 + h0*(5.0*h1 + h0)))), & | |
| 1085 | 0 | h1_4, h2_4, (5.0*h2_4 + h3*(10.0*h2_3 + h3*(10.0*h2_2 + h3*(5.0*h2 + h3)))) /) |
| 1086 | Asys(1:6,6) = (/ -30.0*(h01_2*h01_2), 0.0, & | |
| 1087 | -(6.0*h1_5 + h0*(15.0*h1_4 + h0*(20.0*h1_3 + h0*(15.0*h1_2 + h0*(6.0*h1 + h0))))), & | |
| 1088 | -h1_5, h2_5, & | |
| 1089 | 0 | (6.0*h2_5 + h3*(15.0*h2_4 + h3*(20.0*h2_3 + h3*(15.0*h2_2 + h3*(6.0*h2 + h3))))) /) |
| 1090 | 0 | Bsys(1:6) = (/ 0.0, -2.0, -6.0*h1, 12.0*h1_2, -20.0*h1_3, 30.0*h1_4 /) |
| 1091 | ||
| 1092 | 0 | call linear_solver( 6, Asys, Bsys, Csys ) |
| 1093 | ||
| 1094 | 0 | alpha = Csys(1) |
| 1095 | 0 | beta = Csys(2) |
| 1096 | ||
| 1097 | 0 | tri_l(2) = alpha |
| 1098 | 0 | tri_d(2) = 1.0 |
| 1099 | 0 | tri_u(2) = beta |
| 1100 | 0 | tri_b(2) = Csys(3) * u(1) + Csys(4) * u(2) + Csys(5) * u(3) + Csys(6) * u(4) |
| 1101 | ||
| 1102 | ! Boundary conditions: left boundary | |
| 1103 | 0 | x(1) = 0.0 |
| 1104 | 0 | do i = 1,6 |
| 1105 | 0 | dx = h(i) |
| 1106 | 0 | xavg = x(i) + 0.5 * dx |
| 1107 | Asys(1:6,i) = (/ 1.0, xavg, (xavg**2 + C1_12*dx**2), xavg * (xavg**2 + 0.25*dx**2), & | |
| 1108 | (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4), & | |
| 1109 | 0 | xavg * (xavg**4 + C5_6*xavg**2*dx**2 + 0.0625*dx**4) /) |
| 1110 | 0 | Bsys(i) = u(i) |
| 1111 | 0 | x(i+1) = x(i) + dx |
| 1112 | enddo | |
| 1113 | ||
| 1114 | 0 | call linear_solver( 6, Asys, Bsys, Csys ) |
| 1115 | ||
| 1116 | 0 | tri_d(1) = 0.0 |
| 1117 | 0 | tri_d(1) = 1.0 |
| 1118 | 0 | tri_u(1) = 0.0 |
| 1119 | 0 | tri_b(1) = Csys(2) ! first edge value |
| 1120 | ||
| 1121 | ! Use a left-biased stencil for the second to last row, as described in Eq. (54) of White and Adcroft (2009). | |
| 1122 | ||
| 1123 | ! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes. | |
| 1124 | 0 | hMin = max(h_neglect, h_Min_Frac*((h(N-3) + h(N-2)) + (h(N-1) + h(N)))) |
| 1125 | 0 | h0 = max(h(N-3), hMin) ; h1 = max(h(N-2), hMin) |
| 1126 | 0 | h2 = max(h(N-1), hMin) ; h3 = max(h(N), hMin) |
| 1127 | ||
| 1128 | ! Auxiliary calculations | |
| 1129 | 0 | h1_2 = h1 * h1 ; h1_3 = h1_2 * h1 ; h1_4 = h1_2 * h1_2 ; h1_5 = h1_3 * h1_2 |
| 1130 | 0 | h2_2 = h2 * h2 ; h2_3 = h2_2 * h2 ; h2_4 = h2_2 * h2_2 ; h2_5 = h2_3 * h2_2 |
| 1131 | ||
| 1132 | 0 | h23 = h2 + h3 ; h23_2 = h23 * h23 |
| 1133 | ||
| 1134 | ! Compute matrix entries | |
| 1135 | 0 | Asys(1:6,1) = (/ 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 /) |
| 1136 | 0 | Asys(1:6,2) = (/ 2.0, 2.0, (2.0*h1 + h0), h1, -h2, -(2.0*h2 + h3) /) |
| 1137 | Asys(1:6,3) = (/ 0.0, -6.0*h23, (3.0*h1_2 + h0*(3.0*h1 + h0)), & | |
| 1138 | 0 | h1_2, h2_2, (3.0*h2_2 + h3*(3.0*h2 + h3)) /) |
| 1139 | Asys(1:6,4) = (/ 0.0, -12.0*h23_2, -(4.0*h1_3 + h0*(6.0*h1_2 + h0*(4.0*h1 + h0))), & | |
| 1140 | 0 | -h1_3, h2_3, (4.0*h2_3 + h3*(6.0*h2_2 + h3*(4.0*h2 + h3))) /) |
| 1141 | Asys(1:6,5) = (/ 0.0, -20.0*(h23*h23_2), (5.0*h1_4 + h0*(10.0*h1_3 + h0*(10.0*h1_2 + h0*(5.0*h1 + h0)))), & | |
| 1142 | 0 | h1_4, h2_4, (5.0*h2_4 + h3*(10.0*h2_3 + h3*(10.0*h2_2 + h3*(5.0*h2 + h3)))) /) |
| 1143 | Asys(1:6,6) = (/ 0.0, -30.0*(h23_2*h23_2), & | |
| 1144 | -(6.0*h1_5 + h0*(15.0*h1_4 + h0*(20.0*h1_3 + h0*(15.0*h1_2 + h0*(6.0*h1 + h0))))), & | |
| 1145 | -h1_5, h2_5, & | |
| 1146 | 0 | (6.0*h2_5 + h3*(15.0*h2_4 + h3*(20.0*h2_3 + h3*(15.0*h2_2 + h3*(6.0*h2 + h3))))) /) |
| 1147 | 0 | Bsys(1:6) = (/ 0.0, -2.0, 6.0*h2, 12.0*h2_2, 20.0*h2_3, 30.0*h2_4 /) |
| 1148 | ||
| 1149 | 0 | call linear_solver( 6, Asys, Bsys, Csys ) |
| 1150 | ||
| 1151 | 0 | alpha = Csys(1) |
| 1152 | 0 | beta = Csys(2) |
| 1153 | ||
| 1154 | 0 | tri_l(N) = alpha |
| 1155 | 0 | tri_d(N) = 1.0 |
| 1156 | 0 | tri_u(N) = beta |
| 1157 | 0 | tri_b(N) = Csys(3) * u(N-3) + Csys(4) * u(N-2) + Csys(5) * u(N-1) + Csys(6) * u(N) |
| 1158 | ||
| 1159 | ! Boundary conditions: right boundary | |
| 1160 | 0 | x(1) = 0.0 |
| 1161 | 0 | do i = 1,6 |
| 1162 | 0 | dx = h(N+1-i) |
| 1163 | 0 | xavg = x(i) + 0.5*dx |
| 1164 | Asys(1:6,i) = (/ 1.0, xavg, (xavg**2 + C1_12*dx**2), xavg * (xavg**2 + 0.25*dx**2), & | |
| 1165 | (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4), & | |
| 1166 | 0 | xavg * (xavg**4 + C5_6*xavg**2*dx**2 + 0.0625*dx**4) /) |
| 1167 | 0 | Bsys(i) = u(N+1-i) |
| 1168 | 0 | x(i+1) = x(i) + dx |
| 1169 | enddo | |
| 1170 | ||
| 1171 | 0 | call linear_solver( 6, Asys, Bsys, Csys ) |
| 1172 | ||
| 1173 | 0 | tri_l(N+1) = 0.0 |
| 1174 | 0 | tri_d(N+1) = 1.0 |
| 1175 | 0 | tri_u(N+1) = 0.0 |
| 1176 | 0 | tri_b(N+1) = -Csys(2) |
| 1177 | ||
| 1178 | ! Solve tridiagonal system and assign edge values | |
| 1179 | 0 | call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, N+1 ) |
| 1180 | ||
| 1181 | 0 | do i = 2,N |
| 1182 | 0 | edge_slopes(i,1) = tri_x(i) |
| 1183 | 0 | edge_slopes(i-1,2) = tri_x(i) |
| 1184 | enddo | |
| 1185 | 0 | edge_slopes(1,1) = tri_x(1) |
| 1186 | 0 | edge_slopes(N,2) = tri_x(N+1) |
| 1187 | ||
| 1188 | 0 | end subroutine edge_slopes_implicit_h5 |
| 1189 | ||
| 1190 | ||
| 1191 | !> Compute ih6 edge values (implicit sixth order accurate) in the same units as u. | |
| 1192 | !! | |
| 1193 | !! Sixth-order implicit estimates of edge values are based on a four-cell, | |
| 1194 | !! three-edge stencil. A tridiagonal system is set up and is based on | |
| 1195 | !! expressing the edge values in terms of neighboring cell averages. | |
| 1196 | !! | |
| 1197 | !! The generic relationship is | |
| 1198 | !! | |
| 1199 | !! \f[ | |
| 1200 | !! \alpha u_{i-1/2} + u_{i+1/2} + \beta u_{i+3/2} = | |
| 1201 | !! a \bar{u}_{i-1} + b \bar{u}_i + c \bar{u}_{i+1} + d \bar{u}_{i+2} | |
| 1202 | !! \f] | |
| 1203 | !! | |
| 1204 | !! and the stencil looks like this | |
| 1205 | !! | |
| 1206 | !! i-1 i i+1 i+2 | |
| 1207 | !! ..--o------o------o------o------o--.. | |
| 1208 | !! i-1/2 i+1/2 i+3/2 | |
| 1209 | !! | |
| 1210 | !! In this routine, the coefficients \f$\alpha\f$, \f$\beta\f$, a, b, c and d are | |
| 1211 | !! computed, the tridiagonal system is built, boundary conditions are prescribed and | |
| 1212 | !! the system is solved to yield edge-value estimates. This scheme is described in detail | |
| 1213 | !! by White and Adcroft, 2009, J. Comp. Phys, https://doi.org/10.1016/j.jcp.2008.04.026 | |
| 1214 | !! | |
| 1215 | !! Note that the centered stencil only applies to edges 3 to N-1 (edges are | |
| 1216 | !! numbered 1 to n+1), which yields N-3 equations for N+1 unknowns. Two other | |
| 1217 | !! equations are written by using a right-biased stencil for edge 2 and a | |
| 1218 | !! left-biased stencil for edge N. The prescription of boundary conditions | |
| 1219 | !! (using sixth-order polynomials) closes the system. | |
| 1220 | !! | |
| 1221 | !! CAUTION: For each edge, in order to determine the coefficients of the | |
| 1222 | !! implicit expression, a 6x6 linear system is solved. This may | |
| 1223 | !! become computationally expensive if regridding is carried out | |
| 1224 | !! often. Figuring out closed-form expressions for these coefficients | |
| 1225 | !! on nonuniform meshes turned out to be intractable. | |
| 1226 | 0 | subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect, answer_date ) |
| 1227 | integer, intent(in) :: N !< Number of cells | |
| 1228 | real, dimension(N), intent(in) :: h !< cell widths [H] | |
| 1229 | real, dimension(N), intent(in) :: u !< cell average properties (size N) in arbitrary units [A] | |
| 1230 | real, dimension(N,2), intent(inout) :: edge_val !< Returned edge values [A]; the second index | |
| 1231 | !! is for the two edges of each cell. | |
| 1232 | real, intent(in) :: h_neglect !< A negligibly small width [H] | |
| 1233 | integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use | |
| 1234 | ||
| 1235 | ! Local variables | |
| 1236 | real :: h0, h1, h2, h3 ! cell widths [H] | |
| 1237 | real :: hMin ! The minimum thickness used in these calculations [H] | |
| 1238 | real :: h01, h01_2, h01_3 ! Summed thicknesses to various powers [H^n ~> m^n or kg^n m-2n] | |
| 1239 | real :: h23, h23_2, h23_3 ! Summed thicknesses to various powers [H^n ~> m^n or kg^n m-2n] | |
| 1240 | real :: h1_2, h2_2, h1_3, h2_3 ! Cell widths raised to the 2nd and 3rd powers [H2] or [H3] | |
| 1241 | real :: h1_4, h2_4, h1_5, h2_5 ! Cell widths raised to the 4th and 5th powers [H4] or [H5] | |
| 1242 | real :: alpha, beta ! stencil coefficients [nondim] | |
| 1243 | real, dimension(7) :: x ! Coordinate system with 0 at edges in the same units as h [H] | |
| 1244 | real, parameter :: C1_12 = 1.0 / 12.0 ! A rational parameter [nondim] | |
| 1245 | real, parameter :: C5_6 = 5.0 / 6.0 ! A rational parameter [nondim] | |
| 1246 | real :: dx, xavg ! Differences and averages of successive values of x [H] | |
| 1247 | real, dimension(6,6) :: Asys ! The matrix that is being inverted for a solution, | |
| 1248 | ! in units that might vary with the second (j) index as [H^j] | |
| 1249 | real, dimension(6) :: Bsys ! The right hand side of the system to solve for C in various | |
| 1250 | ! units that sometimes vary with the intex (j) as [H^(j-1)] or [H^j] | |
| 1251 | ! or might be [A] | |
| 1252 | real, dimension(6) :: Csys ! The solution to a matrix equation, which might be [nondim] or the | |
| 1253 | ! coefficients of a fit polynomial in units that vary with the | |
| 1254 | ! index (j) as [A H^(j-1)] | |
| 1255 | 0 | real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal) [nondim] |
| 1256 | 0 | tri_d, & ! trid. system (middle diagonal) [nondim] |
| 1257 | 0 | tri_u, & ! trid. system (upper diagonal) [nondim] |
| 1258 | 0 | tri_b, & ! trid. system (rhs) [A] |
| 1259 | 0 | tri_x ! trid. system (unknowns vector) [A] |
| 1260 | integer :: i, k ! loop indexes | |
| 1261 | ||
| 1262 | ! Loop on interior cells | |
| 1263 | 0 | do k = 2,N-2 |
| 1264 | ! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes. | |
| 1265 | 0 | hMin = max(h_neglect, hMinFrac*((h(k-1) + h(k)) + (h(k+1) + h(k+2)))) |
| 1266 | 0 | h0 = max(h(k-1), hMin) ; h1 = max(h(k), hMin) |
| 1267 | 0 | h2 = max(h(k+1), hMin) ; h3 = max(h(k+2), hMin) |
| 1268 | ||
| 1269 | ! Auxiliary calculations | |
| 1270 | 0 | h1_2 = h1 * h1 ; h1_3 = h1_2 * h1 ; h1_4 = h1_2 * h1_2 ; h1_5 = h1_3 * h1_2 |
| 1271 | 0 | h2_2 = h2 * h2 ; h2_3 = h2_2 * h2 ; h2_4 = h2_2 * h2_2 ; h2_5 = h2_3 * h2_2 |
| 1272 | ||
| 1273 | ! Compute matrix entries as described in Eq. (48) of White and Adcroft (2009) | |
| 1274 | 0 | Asys(1:6,1) = (/ 1.0, 1.0, -1.0, -1.0, -1.0, -1.0 /) |
| 1275 | 0 | Asys(1:6,2) = (/ -2.0*h1, 2.0*h2, (2.0*h1 + h0), h1, -h2, -(2.0*h2 + h3) /) |
| 1276 | Asys(1:6,3) = (/ 3.0*h1_2, 3.0*h2_2, -(3.0*h1_2 + h0*(3.0*h1 + h0)), & ! = -((h0+h1)**3 - h1**3) / h0 | |
| 1277 | 0 | -h1_2, -h2_2, -(3.0*h2_2 + h3*(3.0*h2 + h3)) /) ! = -((h2+h3)**3 - h2**3) / h3 |
| 1278 | Asys(1:6,4) = (/ -4.0*h1_3, 4.0*h2_3, (4.0*h1_3 + h0*(6.0*h1_2 + h0*(4.0*h1 + h0))), & | |
| 1279 | 0 | h1_3, -h2_3, -(4.0*h2_3 + h3*(6.0*h2_2 + h3*(4.0*h2 + h3))) /) |
| 1280 | Asys(1:6,5) = (/ 5.0*h1_4, 5.0*h2_4, -(5.0*h1_4 + h0*(10.0*h1_3 + h0*(10.0*h1_2 + h0*(5.0*h1 + h0)))), & | |
| 1281 | 0 | -h1_4, -h2_4, -(5.0*h2_4 + h3*(10.0*h2_3 + h3*(10.0*h2_2 + h3*(5.0*h2 + h3)))) /) |
| 1282 | Asys(1:6,6) = (/ -6.0*h1_5, 6.0*h2_5, & | |
| 1283 | (6.0*h1_5 + h0*(15.0*h1_4 + h0*(20.0*h1_3 + h0*(15.0*h1_2 + h0*(6.0*h1 + h0))))), & | |
| 1284 | h1_5, -h2_5, & | |
| 1285 | 0 | -(6.0*h2_5 + h3*(15.0*h2_4 + h3*(20.0*h2_3 + h3*(15.0*h2_2 + h3*(6.0*h2 + h3))))) /) |
| 1286 | 0 | Bsys(1:6) = (/ -1.0, 0.0, 0.0, 0.0, 0.0, 0.0 /) |
| 1287 | ||
| 1288 | 0 | call linear_solver( 6, Asys, Bsys, Csys ) |
| 1289 | ||
| 1290 | 0 | alpha = Csys(1) |
| 1291 | 0 | beta = Csys(2) |
| 1292 | ||
| 1293 | 0 | tri_l(k+1) = alpha |
| 1294 | 0 | tri_d(k+1) = 1.0 |
| 1295 | 0 | tri_u(k+1) = beta |
| 1296 | 0 | tri_b(k+1) = Csys(3) * u(k-1) + Csys(4) * u(k) + Csys(5) * u(k+1) + Csys(6) * u(k+2) |
| 1297 | ||
| 1298 | enddo ! end loop on cells | |
| 1299 | ||
| 1300 | ! Use a right-biased stencil for the second row, as described in Eq. (49) of White and Adcroft (2009). | |
| 1301 | ||
| 1302 | ! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes. | |
| 1303 | 0 | hMin = max(h_neglect, hMinFrac*((h(1) + h(2)) + (h(3) + h(4)))) |
| 1304 | 0 | h0 = max(h(1), hMin) ; h1 = max(h(2), hMin) |
| 1305 | 0 | h2 = max(h(3), hMin) ; h3 = max(h(4), hMin) |
| 1306 | ||
| 1307 | ! Auxiliary calculations | |
| 1308 | 0 | h1_2 = h1 * h1 ; h1_3 = h1_2 * h1 ; h1_4 = h1_2 * h1_2 ; h1_5 = h1_3 * h1_2 |
| 1309 | 0 | h2_2 = h2 * h2 ; h2_3 = h2_2 * h2 ; h2_4 = h2_2 * h2_2 ; h2_5 = h2_3 * h2_2 |
| 1310 | 0 | h01 = h0 + h1 ; h01_2 = h01 * h01 ; h01_3 = h01 * h01_2 |
| 1311 | ||
| 1312 | ! Compute matrix entries | |
| 1313 | 0 | Asys(1:6,1) = (/ 1.0, 1.0, -1.0, -1.0, -1.0, -1.0 /) |
| 1314 | 0 | Asys(1:6,2) = (/ -2.0*h01, 0.0, (2.0*h1 + h0), h1, -h2, -(2.0*h2 + h3) /) |
| 1315 | Asys(1:6,3) = (/ 3.0*h01_2, 0.0, -(3.0*h1_2 + h0*(3.0*h1 + h0)), & | |
| 1316 | 0 | -h1_2, -h2_2, -(3.0*h2_2 + h3*(3.0*h2 + h3)) /) |
| 1317 | Asys(1:6,4) = (/ -4.0*h01_3, 0.0, (4.0*h1_3 + h0*(6.0*h1_2 + h0*(4.0*h1 + h0))), & | |
| 1318 | 0 | h1_3, -h2_3, -(4.0*h2_3 + h3*(6.0*h2_2 + h3*(4.0*h2 + h3))) /) |
| 1319 | Asys(1:6,5) = (/ 5.0*(h01_2*h01_2), 0.0, -(5.0*h1_4 + h0*(10.0*h1_3 + h0*(10.0*h1_2 + h0*(5.0*h1 + h0)))), & | |
| 1320 | 0 | -h1_4, -h2_4, -(5.0*h2_4 + h3*(10.0*h2_3 + h3*(10.0*h2_2 + h3*(5.0*h2 + h3)))) /) |
| 1321 | Asys(1:6,6) = (/ -6.0*(h01_3*h01_2), 0.0, & | |
| 1322 | (6.0*h1_5 + h0*(15.0*h1_4 + h0*(20.0*h1_3 + h0*(15.0*h1_2 + h0*(6.0*h1 + h0))))), & | |
| 1323 | h1_5, - h2_5, & | |
| 1324 | 0 | -(6.0*h2_5 + h3*(15.0*h2_4 + h3*(20.0*h2_3 + h3*(15.0*h2_2 + h3*(6.0*h2 + h3))))) /) |
| 1325 | 0 | Bsys(1:6) = (/ -1.0, 2.0*h1, -3.0*h1_2, 4.0*h1_3, -5.0*h1_4, 6.0*h1_5 /) |
| 1326 | ||
| 1327 | 0 | call linear_solver( 6, Asys, Bsys, Csys ) |
| 1328 | ||
| 1329 | 0 | alpha = Csys(1) |
| 1330 | 0 | beta = Csys(2) |
| 1331 | ||
| 1332 | 0 | tri_l(2) = alpha |
| 1333 | 0 | tri_d(2) = 1.0 |
| 1334 | 0 | tri_u(2) = beta |
| 1335 | 0 | tri_b(2) = Csys(3) * u(1) + Csys(4) * u(2) + Csys(5) * u(3) + Csys(6) * u(4) |
| 1336 | ||
| 1337 | ! Boundary conditions: left boundary | |
| 1338 | 0 | hMin = max( h_neglect, hMinFrac*((h(1)+h(2)) + (h(5)+h(6)) + (h(3)+h(4))) ) |
| 1339 | 0 | x(1) = 0.0 |
| 1340 | 0 | do i = 1,6 |
| 1341 | 0 | dx = max( hMin, h(i) ) |
| 1342 | 0 | xavg = x(i) + 0.5*dx |
| 1343 | Asys(1:6,i) = (/ 1.0, xavg, (xavg**2 + C1_12*dx**2), xavg * (xavg**2 + 0.25*dx**2), & | |
| 1344 | (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4), & | |
| 1345 | 0 | xavg * (xavg**4 + C5_6*xavg**2*dx**2 + 0.0625*dx**4) /) |
| 1346 | 0 | Bsys(i) = u(i) |
| 1347 | 0 | x(i+1) = x(i) + dx |
| 1348 | enddo | |
| 1349 | ||
| 1350 | 0 | call linear_solver( 6, Asys, Bsys, Csys ) |
| 1351 | ||
| 1352 | 0 | tri_l(1) = 0.0 |
| 1353 | 0 | tri_d(1) = 1.0 |
| 1354 | 0 | tri_u(1) = 0.0 |
| 1355 | 0 | tri_b(1) = evaluation_polynomial( Csys, 6, x(1) ) ! first edge value |
| 1356 | ||
| 1357 | ! Use a left-biased stencil for the second to last row, as described in Eq. (50) of White and Adcroft (2009). | |
| 1358 | ||
| 1359 | ! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes. | |
| 1360 | 0 | hMin = max(h_neglect, hMinFrac*((h(N-3) + h(N-2)) + (h(N-1) + h(N)))) |
| 1361 | 0 | h0 = max(h(N-3), hMin) ; h1 = max(h(N-2), hMin) |
| 1362 | 0 | h2 = max(h(N-1), hMin) ; h3 = max(h(N), hMin) |
| 1363 | ||
| 1364 | ! Auxiliary calculations | |
| 1365 | 0 | h1_2 = h1 * h1 ; h1_3 = h1_2 * h1 ; h1_4 = h1_2 * h1_2 ; h1_5 = h1_3 * h1_2 |
| 1366 | 0 | h2_2 = h2 * h2 ; h2_3 = h2_2 * h2 ; h2_4 = h2_2 * h2_2 ; h2_5 = h2_3 * h2_2 |
| 1367 | 0 | h23 = h2 + h3 ; h23_2 = h23 * h23 ; h23_3 = h23 * h23_2 |
| 1368 | ||
| 1369 | ! Compute matrix entries | |
| 1370 | 0 | Asys(1:6,1) = (/ 1.0, 1.0, -1.0, -1.0, -1.0, -1.0 /) |
| 1371 | 0 | Asys(1:6,2) = (/ 0.0, 2.0*h23, (2.0*h1 + h0), h1, -h2, -(2.0*h2 + h3) /) |
| 1372 | Asys(1:6,3) = (/ 0.0, 3.0*h23_2, -(3.0*h1_2 + h0*(3.0*h1 + h0)), & | |
| 1373 | 0 | -h1_2, -h2_2, -(3.0*h2_2 + h3*(3.0*h2 + h3)) /) |
| 1374 | Asys(1:6,4) = (/ 0.0, 4.0*h23_3, (4.0*h1_3 + h0*(6.0*h1_2 + h0*(4.0*h1 + h0))), & | |
| 1375 | 0 | h1_3, -h2_3, -(4.0*h2_3 + h3*(6.0*h2_2 + h3*(4.0*h2 + h3))) /) |
| 1376 | Asys(1:6,5) = (/ 0.0, 5.0*(h23_2*h23_2), -(5.0*h1_4 + h0*(10.0*h1_3 + h0*(10.0*h1_2 + h0*(5.0*h1 + h0)))), & | |
| 1377 | 0 | -h1_4, -h2_4, -(5.0*h2_4 + h3*(10.0*h2_3 + h3*(10.0*h2_2 + h3*(5.0*h2 + h3)))) /) |
| 1378 | Asys(1:6,6) = (/ 0.0, 6.0*(h23_3*h23_2), & | |
| 1379 | (6.0*h1_5 + h0*(15.0*h1_4 + h0*(20.0*h1_3 + h0*(15.0*h1_2 + h0*(6.0*h1 + h0))))), & | |
| 1380 | h1_5, -h2_5, & | |
| 1381 | 0 | -(6.0*h2_5 + h3*(15.0*h2_4 + h3*(20.0*h2_3 + h3*(15.0*h2_2 + h3*(6.0*h2 + h3))))) /) |
| 1382 | 0 | Bsys(1:6) = (/ -1.0, -2.0*h2, -3.0*h2_2, -4.0*h2_3, -5.0*h2_4, -6.0*h2_5 /) |
| 1383 | ||
| 1384 | 0 | call linear_solver( 6, Asys, Bsys, Csys ) |
| 1385 | ||
| 1386 | 0 | alpha = Csys(1) |
| 1387 | 0 | beta = Csys(2) |
| 1388 | ||
| 1389 | 0 | tri_l(N) = alpha |
| 1390 | 0 | tri_d(N) = 1.0 |
| 1391 | 0 | tri_u(N) = beta |
| 1392 | 0 | tri_b(N) = Csys(3) * u(N-3) + Csys(4) * u(N-2) + Csys(5) * u(N-1) + Csys(6) * u(N) |
| 1393 | ||
| 1394 | ! Boundary conditions: right boundary | |
| 1395 | 0 | hMin = max( h_neglect, hMinFrac*(h(N-3) + h(N-2)) + ((h(N-1) + h(N)) + (h(N-5) + h(N-4))) ) |
| 1396 | 0 | x(1) = 0.0 |
| 1397 | 0 | do i = 1,6 |
| 1398 | 0 | dx = max( hMin, h(N+1-i) ) |
| 1399 | 0 | xavg = x(i) + 0.5 * dx |
| 1400 | Asys(1:6,i) = (/ 1.0, xavg, (xavg**2 + C1_12*dx**2), xavg * (xavg**2 + 0.25*dx**2), & | |
| 1401 | (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4), & | |
| 1402 | 0 | xavg * (xavg**4 + C5_6*xavg**2*dx**2 + 0.0625*dx**4) /) |
| 1403 | 0 | Bsys(i) = u(N+1-i) |
| 1404 | 0 | x(i+1) = x(i) + dx |
| 1405 | enddo | |
| 1406 | ||
| 1407 | 0 | call linear_solver( 6, Asys, Bsys, Csys ) |
| 1408 | ||
| 1409 | 0 | tri_l(N+1) = 0.0 |
| 1410 | 0 | tri_d(N+1) = 1.0 |
| 1411 | 0 | tri_u(N+1) = 0.0 |
| 1412 | 0 | tri_b(N+1) = Csys(1) |
| 1413 | ||
| 1414 | ! Solve tridiagonal system and assign edge values | |
| 1415 | 0 | call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, N+1 ) |
| 1416 | ||
| 1417 | 0 | do i = 2,N |
| 1418 | 0 | edge_val(i,1) = tri_x(i) |
| 1419 | 0 | edge_val(i-1,2) = tri_x(i) |
| 1420 | enddo | |
| 1421 | 0 | edge_val(1,1) = tri_x(1) |
| 1422 | 0 | edge_val(N,2) = tri_x(N+1) |
| 1423 | ||
| 1424 | 0 | end subroutine edge_values_implicit_h6 |
| 1425 | ||
| 1426 | ||
| 1427 | !> Test that A*C = R to within a tolerance, issuing a fatal error with an explanatory message if they do not. | |
| 1428 | 0 | subroutine test_line(msg, N, A, C, R, mag, tol) |
| 1429 | character(len=*), intent(in) :: msg !< An identifying message for this test | |
| 1430 | integer, intent(in) :: N !< The number of points in the system | |
| 1431 | real, dimension(4), intent(in) :: A !< One of the two vectors being multiplied in arbitrary units [A] | |
| 1432 | real, dimension(4), intent(in) :: C !< One of the two vectors being multiplied in arbitrary units [B] | |
| 1433 | real, intent(in) :: R !< The expected solution of the equation [A B] | |
| 1434 | real, intent(in) :: mag !< The magnitude of leading order terms in this line [A B] | |
| 1435 | real, optional, intent(in) :: tol !< The fractional tolerance for the sums [nondim] | |
| 1436 | ||
| 1437 | real :: sum, sum_mag ! The sum of the products and their magnitude in arbitrary units [A B] | |
| 1438 | real :: tolerance ! The fractional tolerance for the sums [nondim] | |
| 1439 | character(len=128) :: mesg2 | |
| 1440 | integer :: i | |
| 1441 | ||
| 1442 | 0 | tolerance = 1.0e-12 ; if (present(tol)) tolerance = tol |
| 1443 | ||
| 1444 | 0 | sum = 0.0 ; sum_mag = max(0.0,mag) |
| 1445 | 0 | do i=1,N |
| 1446 | 0 | sum = sum + A(i) * C(i) |
| 1447 | 0 | sum_mag = sum_mag + abs(A(i) * C(i)) |
| 1448 | enddo | |
| 1449 | ||
| 1450 | 0 | if (abs(sum - R) > tolerance * (sum_mag + abs(R))) then |
| 1451 | 0 | write(mesg2, '(", Fractional error = ", es12.4,", sum = ", es12.4)') (sum - R) / (sum_mag + abs(R)), sum |
| 1452 | 0 | call MOM_error(FATAL, "Failed line test: "//trim(msg)//trim(mesg2)) |
| 1453 | endif | |
| 1454 | ||
| 1455 | 0 | end subroutine test_line |
| 1456 | ||
| 1457 | end module regrid_edge_values |