regrid_edge_values.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Edge value estimation for high-order reconstruction
7
8use mom_error_handler, only : mom_error, fatal
12
13implicit none ; private
14
15! -----------------------------------------------------------------------------
16! The following routines are visible to the outside world
17! -----------------------------------------------------------------------------
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.
27real, parameter :: hminfrac = 1.e-5 !< A minimum fraction for min(h)/sum(h) [nondim]
28
29contains
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.
41subroutine 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 use_2018_answers = .true. ; if (present(answer_date)) use_2018_answers = (answer_date < 20190101)
57
58 ! Loop on cells to bound edge value
59 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 km1 = max(1,k-1) ; kp1 = min(k+1,n)
66
67 slope_x_h = 0.0
68 if (use_2018_answers) then
69 sigma_l = 2.0 * ( u(k) - u(km1) ) / ( h(k) + h_neglect )
70 sigma_c = 2.0 * ( u(kp1) - u(km1) ) / ( h(km1) + 2.0*h(k) + h(kp1) + h_neglect )
71 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 if ( (sigma_l * sigma_r) > 0.0 ) &
76 slope_x_h = 0.5 * h(k) * sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c )
77 elseif ( ((h(km1) + h(kp1)) + 2.0*h(k)) > 0.0 ) then
78 sigma_l = ( u(k) - u(km1) )
79 sigma_c = ( u(kp1) - u(km1) ) * ( h(k) / ((h(km1) + h(kp1)) + 2.0*h(k)) )
80 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 if ( (sigma_l * sigma_r) > 0.0 ) &
85 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 if ( (u(km1)-edge_val(k,1)) * (edge_val(k,1)-u(k)) < 0.0 ) then
90 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 if ( (u(kp1)-edge_val(k,2)) * (edge_val(k,2)-u(k)) < 0.0 ) then
94 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 edge_val(k,1) = max( min( edge_val(k,1), max(u(km1), u(k)) ), min(u(km1), u(k)) )
99 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
103end 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.
109subroutine 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 do k = 1,n-1
119 ! Compare edge values on the right and left sides of the edge
120 if ( edge_val(k,2) /= edge_val(k+1,1) ) then
121 u0_avg = 0.5 * ( edge_val(k,2) + edge_val(k+1,1) )
122 edge_val(k,2) = u0_avg
123 edge_val(k+1,1) = u0_avg
124 endif
125
126 enddo ! end loop on interior edges
127
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.
134subroutine 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 do k = 1,n-1
144 if ( (edge_val(k+1,1) - edge_val(k,2)) * (u(k+1) - u(k)) < 0.0 ) then
145 u0_avg = 0.5 * ( edge_val(k,2) + edge_val(k+1,1) )
146 u0_avg = max( min( u0_avg, max(u(k), u(k+1)) ), min(u(k), u(k+1)) )
147 edge_val(k,2) = u0_avg
148 edge_val(k+1,1) = u0_avg
149 endif
150 enddo ! end loop on interior edges
151
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.
168subroutine 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 edge_val(1,1) = u(1)
180 edge_val(n,2) = u(n)
181
182 do k = 2,n
183 ! Compute left edge value
184 if (h(k-1) + h(k) == 0.0) then ! Avoid singularities when h0+h1=0
185 edge_val(k,1) = 0.5 * (u(k-1) + u(k))
186 else
187 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 edge_val(k-1,2) = edge_val(k,1)
192 enddo
193
194end 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.
215subroutine 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 use_2018_answers = .true. ; if (present(answer_date)) use_2018_answers = (answer_date < 20190101)
248
249 ! Loop on interior cells
250 do i = 3,n-1
251
252 h0 = h(i-2)
253 h1 = h(i-1)
254 h2 = h(i)
255 h3 = h(i+1)
256
257 ! Avoid singularities when consecutive pairs of h vanish
258 if (h0+h1==0.0 .or. h1+h2==0.0 .or. h2+h3==0.0) then
259 if (use_2018_answers) then
260 h_min = hminfrac*max( h_neglect, h0+h1+h2+h3 )
261 else
262 h_min = hminfrac*max( h_neglect, (h0+h1)+(h2+h3) )
263 endif
264 h0 = max( h_min, h(i-2) )
265 h1 = max( h_min, h(i-1) )
266 h2 = max( h_min, h(i) )
267 h3 = max( h_min, h(i+1) )
268 endif
269
270 if (use_2018_answers) then
271 f1 = (h0+h1) * (h2+h3) / (h1+h2)
272 f2 = h2 * u(i-1) + h1 * u(i)
273 f3 = 1.0 / (h0+h1+h2) + 1.0 / (h1+h2+h3)
274 et1 = f1 * f2 * f3
275 et2 = ( h2 * (h2+h3) / ( (h0+h1+h2)*(h0+h1) ) ) * &
276 ((h0+2.0*h1) * u(i-1) - h1 * u(i-2))
277 et3 = ( h1 * (h0+h1) / ( (h1+h2+h3)*(h2+h3) ) ) * &
278 ((2.0*h2+h3) * u(i) - h2 * u(i+1))
279 edge_val(i,1) = (et1 + et2 + et3) / ( h0 + h1 + h2 + h3)
280 else
281 i_h12 = 1.0 / (h1+h2)
282 i_den_et2 = 1.0 / ( ((h0+h1)+h2)*(h0+h1) ) ; i_h012 = (h0+h1) * i_den_et2
283 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 ( 1.0 + (h2 * i_h123 + (h2+h3) * i_h012) ) * i_h12 * (h1*(h0+h1)) * u(i)
287 et2 = ( h1 * (h2*(h2+h3)) * i_den_et2 ) * (u(i-1)-u(i-2))
288 et3 = ( h2 * (h1*(h0+h1)) * i_den_et3 ) * (u(i) - u(i+1))
289 edge_val(i,1) = (et1 + (et2 + et3)) / ((h0 + h1) + (h2 + h3))
290 endif
291 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 if (use_2018_answers) then
297 h_min = max( h_neglect, hminfrac*sum(h(1:4)) )
298 x(1) = 0.0
299 do i = 1,4
300 dx = max(h_min, h(i) )
301 x(i+1) = x(i) + dx
302 do j = 1,4 ; a(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j) ; enddo
303 b(i) = u(i) * dx
304 enddo
305
306 call solve_linear_system( a, b, c, 4 )
307
308 ! Set the edge values of the first cell
309 edge_val(1,1) = evaluation_polynomial( c, 4, x(1) )
310 edge_val(1,2) = evaluation_polynomial( c, 4, x(2) )
311 else ! Use expressions with less sensitivity to roundoff
312 do i=1,4 ; dz(i) = max(h_neglect, h(i) ) ; u_tmp(i) = u(i) ; enddo
313 call end_value_h4(dz, u_tmp, c)
314
315 ! Set the edge values of the first cell
316 edge_val(1,1) = c(1)
317 edge_val(1,2) = c(1) + dz(1)*(c(2) + dz(1)*(c(3) + dz(1)*c(4)))
318 endif
319 edge_val(2,1) = edge_val(1,2)
320
321 ! Determine two edge values of the last cell
322 if (use_2018_answers) then
323 h_min = max( h_neglect, hminfrac*sum(h(n-3:n)) )
324
325 x(1) = 0.0
326 do i = 1,4
327 dx = max(h_min, h(n-4+i) )
328 x(i+1) = x(i) + dx
329 do j = 1,4 ; a(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j) ; enddo
330 b(i) = u(n-4+i) * dx
331 enddo
332
333 call solve_linear_system( a, b, c, 4 )
334
335 ! Set the last and second to last edge values
336 edge_val(n,2) = evaluation_polynomial( c, 4, x(5) )
337 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 do i=1,4 ; dz(i) = max(h_neglect, h(n+1-i) ) ; u_tmp(i) = u(n+1-i) ; enddo
342 call end_value_h4(dz, u_tmp, c)
343
344 ! Set the last and second to last edge values
345 edge_val(n,2) = c(1)
346 edge_val(n,1) = c(1) + dz(1)*(c(2) + dz(1)*(c(3) + dz(1)*c(4)))
347 endif
348 edge_val(n-1,2) = edge_val(n,1)
349
350end 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.
368subroutine 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 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 real :: au(n) ! Scalar field difference across each cell [A]
383 real :: al(n), ar(n) ! Scalar field at the left and right edges of a cell [A]
384 real :: h112(n+1), h122(n+1) ! Combinations of thicknesses [H ~> m or kg m-2]
385 real :: i_h12(n+1) ! Inverses of combinations of thickesses [H-1 ~> m-1 or m2 kg-1]
386 real :: h2_h123(n) ! A ratio of a layer thickness of the sum of 3 adjacent thicknesses [nondim]
387 real :: i_h0123(n) ! Inverse of the sum of 4 adjacent thicknesses [H-1 ~> m-1 or m2 kg-1]
388 real :: h01_h112(n+1) ! A ratio of sums of adjacent thicknesses [nondim], 2/3 in the limit of uniform thicknesses.
389 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 do k=1,n ; dp(k) = max(h(k), h_neglect) ; enddo
394
395 !compute grid metrics
396 do k=2,n
397 h112(k) = 2.*dp(k-1) + dp(k)
398 h122(k) = dp(k-1) + 2.*dp(k)
399 i_h12(k) = 1.0 / (dp(k-1) + dp(k))
400 enddo !k
401 do k=2,n-1
402 h2_h123(k) = dp(k) / (dp(k) + (dp(k-1)+dp(k+1)))
403 enddo
404 do k=3,n-1
405 i_h0123(k) = 1.0 / ((dp(k-2) + dp(k-1)) + (dp(k) + dp(k+1)))
406
407 h01_h112(k) = (dp(k-2) + dp(k-1)) / (2.0*dp(k-1) + dp(k))
408 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 au(1) = 0.
413 do k=2,n-1
414 slk = u(k )-u(k-1)
415 srk = u(k+1)-u(k)
416 if (slk*srk > 0.) then
417 sck = h2_h123(k)*( h112(k)*srk*i_h12(k+1) + h122(k+1)*slk*i_h12(k) )
418 au(k) = sign(min(abs(2.0*slk), abs(sck), abs(2.0*srk)), sck)
419 else
420 au(k) = 0.
421 endif
422 enddo !k
423 au(n) = 0.
424
425 !Compute "first guess" edge values: Colella, Eq. (1.6)
426 al(1) = u(1) ! 1st layer PCM
427 ar(1) = u(1) ! 1st layer PCM
428 al(2) = u(1) ! 1st layer PCM
429 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 + (dp(k)*au(k-1)*h23_h122(k) - dp(k-1)*au(k)*h01_h112(k)) )
435 ar(k-1) = al(k)
436 enddo !k
437 ar(n-1) = u(n) ! last layer PCM
438 al(n) = u(n) ! last layer PCM
439 ar(n) = u(n) ! last layer PCM
440
441 !Set coefficients
442 do k=1,n
443 edge_val(k,1) = al(k)
444 edge_val(k,2) = ar(k)
445 enddo !k
446
447end 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.
475subroutine 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 real, dimension(N+1) :: tri_l, & ! tridiagonal system (lower diagonal) [nondim]
506 tri_d, & ! tridiagonal system (middle diagonal) [nondim]
507 tri_c, & ! tridiagonal system central value [nondim], with tri_d = tri_c+tri_l+tri_u
508 tri_u, & ! tridiagonal system (upper diagonal) [nondim]
509 tri_b, & ! tridiagonal system (right hand side) [A]
510 tri_x ! tridiagonal system (solution vector) [A]
511 logical :: use_2018_answers ! If true use older, less accurate expressions.
512
513 use_2018_answers = .true. ; if (present(answer_date)) use_2018_answers = (answer_date < 20190101)
514
515 ! Loop on cells (except last one)
516 do i = 1,n-1
517 if (use_2018_answers) then
518 ! Get cell widths
519 h0 = h(i)
520 h1 = h(i+1)
521 ! Avoid singularities when h0+h1=0
522 if (h0+h1==0.) then
523 h0 = h_neglect
524 h1 = h_neglect
525 endif
526
527 ! Auxiliary calculations
528 h0ph1_2 = (h0 + h1)**2
529 h0ph1_4 = h0ph1_2**2
530 h0h1 = h0 * h1
531 h0_2 = h0 * h0
532 h1_2 = h1 * h1
533
534 ! Coefficients
535 alpha = h1_2 / h0ph1_2
536 beta = h0_2 / h0ph1_2
537 a = 2.0 * h1_2 * ( h1_2 + 2.0 * h0_2 + 3.0 * h0h1 ) / h0ph1_4
538 b = 2.0 * h0_2 * ( h0_2 + 2.0 * h1_2 + 3.0 * h0h1 ) / h0ph1_4
539
540 tri_d(i+1) = 1.0
541 else ! Use expressions with less sensitivity to roundoff
542 ! Get cell widths
543 h0 = max(h(i), h_neglect)
544 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 if (abs(h0) < 1.0e-12*abs(h1)) h0 = 1.0e-12*h1
548 if (abs(h1) < 1.0e-12*abs(h0)) h1 = 1.0e-12*h0
549 i_h2 = 1.0 / ((h0 + h1)**2)
550 alpha = (h1 * h1) * i_h2
551 beta = (h0 * h0) * i_h2
552 abmix = (h0 * h1) * i_h2
553 a = 2.0 * alpha * ( alpha + 2.0 * beta + 3.0 * abmix )
554 b = 2.0 * beta * ( beta + 2.0 * alpha + 3.0 * abmix )
555
556 tri_c(i+1) = 2.0*abmix ! = 1.0 - alpha - beta
557 endif
558
559 tri_l(i+1) = alpha
560 tri_u(i+1) = beta
561
562 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 if (use_2018_answers) then
568 h_min = max( h_neglect, hminfrac*sum(h(1:4)) )
569 x(1) = 0.0
570 do i = 1,4
571 dx = max(h_min, h(i) )
572 x(i+1) = x(i) + dx
573 do j = 1,4 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo
574 bsys(i) = u(i) * dx
575 enddo
576
577 call solve_linear_system( asys, bsys, csys, 4 )
578
579 tri_b(1) = evaluation_polynomial( csys, 4, x(1) ) ! Set the first edge value
580 tri_d(1) = 1.0
581 else ! Use expressions with less sensitivity to roundoff
582 do i=1,4 ; dz(i) = max(h_neglect, h(i) ) ; u_tmp(i) = u(i) ; enddo
583 call end_value_h4(dz, u_tmp, csys)
584
585 tri_b(1) = csys(1) ! Set the first edge value.
586 tri_c(1) = 1.0
587 endif
588 tri_u(1) = 0.0 ! tri_l(1) = 0.0
589
590 ! Boundary conditions: set the last boundary value
591 if (use_2018_answers) then
592 h_min = max( h_neglect, hminfrac*sum(h(n-3:n)) )
593 x(1) = 0.0
594 do i=1,4
595 dx = max(h_min, h(n-4+i) )
596 x(i+1) = x(i) + dx
597 do j = 1,4 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo
598 bsys(i) = u(n-4+i) * dx
599 enddo
600
601 call solve_linear_system( asys, bsys, csys, 4 )
602
603 ! Set the last edge value
604 tri_b(n+1) = evaluation_polynomial( csys, 4, x(5) )
605 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 do i=1,4 ; dz(i) = max(h_neglect, h(n+1-i) ) ; u_tmp(i) = u(n+1-i) ; enddo
611 call end_value_h4(dz, u_tmp, csys)
612
613 tri_b(n+1) = csys(1) ! Set the last edge value
614 tri_c(n+1) = 1.0
615 endif
616 tri_l(n+1) = 0.0 ! tri_u(N+1) = 0.0
617
618 ! Solve tridiagonal system and assign edge values
619 if (use_2018_answers) then
620 call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, n+1 )
621 else
622 call solve_diag_dominant_tridiag( tri_l, tri_c, tri_u, tri_b, tri_x, n+1 )
623 endif
624
625 edge_val(1,1) = tri_x(1)
626 do i=2,n
627 edge_val(i,1) = tri_x(i)
628 edge_val(i-1,2) = tri_x(i)
629 enddo
630 edge_val(n,2) = tri_x(n+1)
631
632end 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.
636subroutine 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 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 if ((h2+h3) < min_frac*h1) h3 = min_frac*h1 - h2
681 if ((h3+h4) < min_frac*h1) h4 = min_frac*h1 - h3
682
683 h12 = h1+h2 ; h23 = h2+h3 ; h34 = h3+h4
684 h123 = h12 + h3 ; h234 = h2 + h34 ; h1234 = h12 + h34
685 ! Find 3 reciprocals with a single division for efficiency.
686 i_denb3 = 1.0 / (h123 * h12 * h23)
687 i_h12 = (h123 * h23) * i_denb3
688 i_h23 = (h12 * h123) * i_denb3
689 i_h123 = (h12 * h23) * i_denb3
690 i_denom = 1.0 / ( h1234 * (h234 * h34) )
691 i_h34 = (h1234 * h234) * i_denom
692 i_h234 = (h1234 * h34) * i_denom
693 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 wt(1,1) = -h1 * (i_h1234 + i_h123 + i_h12) ! > -3
713 wt(2,1) = h1 * h12 * ( i_h234 * i_h1234 + i_h23 * (i_h234 + i_h123) ) ! < (h1/h234) + (h1/h23)*(2+(h1/h234))
714 wt(3,1) = -h1 * h12 * h123 * i_denom ! > -(h1/h34)*(1+(h1/h234))
715
716 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 (h1+h12) * ( i_h1234*i_h234 + i_h23 * (i_h234 + i_h123) ) )
719 wt(3,2) = 2.0 * ((h1+h12) * h123 + h1*h12 ) * i_denom ! < (2+(6*h1/h234)) / h34
720
721 wt(1,3) = -3.0 * i_h12 * i_h123* ( 1.0 + i_h1234 * ((h1+h12)+h123) ) ! > -12 / (h12*h123)
722 wt(2,3) = 3.0 * i_h23 * ( i_h123 + i_h1234 * ((h1+h12)+h123) * (i_h123 + i_h234) ) ! < 12 / (h23^2)
723 wt(3,3) = -3.0 * ((h1+h12)+h123) * i_denom ! > -9 / (h234*h23)
724
725 wt(1,4) = 4.0 * i_h1234 * i_h123 * i_h12 ! Wt*h1^3 < 4
726 wt(2,4) = -4.0 * i_h1234 * (i_h23 * (i_h123 + i_h234)) ! Wt*h1^3 > -4* (h1/h23)*(1+h1/h234)
727 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 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 csys(2) = ((wt(1,2) * (u(2)-u(1))) + (wt(2,2) * (u(3)-u(2)))) + (wt(3,2) * (u(4)-u(3)))
731 csys(3) = ((wt(1,3) * (u(2)-u(1))) + (wt(2,3) * (u(3)-u(2)))) + (wt(3,3) * (u(4)-u(3)))
732 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
749end 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.
781subroutine 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 real, dimension(N+1) :: tri_l, & ! tridiagonal system (lower diagonal) [nondim]
813 tri_d, & ! tridiagonal system (middle diagonal) [nondim]
814 tri_c, & ! tridiagonal system central value [nondim], with tri_d = tri_c+tri_l+tri_u
815 tri_u, & ! tridiagonal system (upper diagonal) [nondim]
816 tri_b, & ! tridiagonal system (right hand side) [A H-1]
817 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 hneglect3 = h_neglect**3
822 use_2018_answers = .true. ; if (present(answer_date)) use_2018_answers = (answer_date < 20190101)
823
824 ! Loop on cells (except last one)
825 do i = 1,n-1
826
827 if (use_2018_answers) then
828 ! Get cell widths
829 h0 = h(i)
830 h1 = h(i+1)
831
832 ! Auxiliary calculations
833 h0h1 = h0 * h1
834 h0_2 = h0 * h0
835 h1_2 = h1 * h1
836 h0_3 = h0_2 * h0
837 h1_3 = h1_2 * h1
838
839 d = 4.0 * h0h1 * ( h0 + h1 ) + h1_3 + h0_3
840
841 ! Coefficients
842 alpha = h1 * (h0_2 + h0h1 - h1_2) / ( d + hneglect3 )
843 beta = h0 * (h1_2 + h0h1 - h0_2) / ( d + hneglect3 )
844 a = -12.0 * h0h1 / ( d + hneglect3 )
845 b = -a
846
847 tri_l(i+1) = alpha
848 tri_d(i+1) = 1.0
849 tri_u(i+1) = beta
850
851 tri_b(i+1) = a * u(i) + b * u(i+1)
852 else
853 ! Get cell widths
854 h0 = max(h(i), h_neglect)
855 h1 = max(h(i+1), h_neglect)
856
857 i_h = 1.0 / (h0 + h1)
858 h0 = h0 * i_h ; h1 = h1 * i_h
859
860 h0h1 = h0 * h1 ; h0_2 = h0 * h0 ; h1_2 = h1 * h1
861 h0_3 = h0_2 * h0 ; h1_3 = h1_2 * h1
862
863 ! Set the tridiagonal coefficients
864 i_d = 1.0 / (4.0 * h0h1 * ( h0 + h1 ) + h1_3 + h0_3) ! = 1 / ((h0 + h1)**3 + h0*h1*(h0 + h1))
865 tri_l(i+1) = (h1 * ((h0_2 + h0h1) - h1_2)) * i_d
866 ! tri_d(i+1) = 1.0
867 tri_c(i+1) = 2.0 * ((h0_2 + h1_2) * (h0 + h1)) * i_d
868 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 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 if (use_2018_answers) then
882 x(1) = 0.0
883 do i = 1,4
884 dx = h(i)
885 x(i+1) = x(i) + dx
886 do j = 1,4 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo
887 bsys(i) = u(i) * dx
888 enddo
889
890 call solve_linear_system( asys, bsys, csys, 4 )
891
892 dsys(1) = csys(2) ; dsys(2) = 2.0 * csys(3) ; dsys(3) = 3.0 * csys(4)
893 tri_b(1) = evaluation_polynomial( dsys, 3, x(1) ) ! Set the first edge slope
894 tri_d(1) = 1.0
895 else ! Use expressions with less sensitivity to roundoff
896 do i=1,4 ; dz(i) = max(h_neglect, h(i) ) ; u_tmp(i) = u(i) ; enddo
897 call end_value_h4(dz, u_tmp, csys)
898
899 ! Set the first edge slope
900 tri_b(1) = csys(2)
901 tri_c(1) = 1.0
902 endif
903 tri_u(1) = 0.0 ! tri_l(1) = 0.0
904
905 ! Boundary conditions: set the last edge slope
906 if (use_2018_answers) then
907 x(1) = 0.0
908 do i = 1,4
909 dx = h(n-4+i)
910 x(i+1) = x(i) + dx
911 do j = 1,4 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo
912 bsys(i) = u(n-4+i) * dx
913 enddo
914
915 call solve_linear_system( asys, bsys, csys, 4 )
916
917 dsys(1) = csys(2) ; dsys(2) = 2.0 * csys(3) ; dsys(3) = 3.0 * csys(4)
918 ! Set the last edge slope
919 tri_b(n+1) = evaluation_polynomial( dsys, 3, x(5) )
920 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 do i=1,4 ; dz(i) = max(h_neglect, h(n+1-i) ) ; u_tmp(i) = u(n+1-i) ; enddo
925
926 call end_value_h4(dz, u_tmp, csys)
927
928 ! Set the last edge slope
929 tri_b(n+1) = -csys(2)
930 tri_c(n+1) = 1.0
931 endif
932 tri_l(n+1) = 0.0 ! tri_u(N+1) = 0.0
933
934 ! Solve tridiagonal system and assign edge slopes
935 if (use_2018_answers) then
936 call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, n+1 )
937 else
938 call solve_diag_dominant_tridiag( tri_l, tri_c, tri_u, tri_b, tri_x, n+1 )
939 endif
940
941 do i = 2,n
942 edge_slopes(i,1) = tri_x(i)
943 edge_slopes(i-1,2) = tri_x(i)
944 enddo
945 edge_slopes(1,1) = tri_x(1)
946 edge_slopes(n,2) = tri_x(n+1)
947
948end subroutine edge_slopes_implicit_h3
949
950
951!------------------------------------------------------------------------------
952!> Compute ih5 edge slopes (implicit fifth order accurate)
953subroutine 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 real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal) [nondim]
1017 tri_d, & ! trid. system (middle diagonal) [nondim]
1018 tri_u, & ! trid. system (upper diagonal) [nondim]
1019 tri_b, & ! trid. system (rhs) [A H-1]
1020 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 do k = 2,n-2
1026 ! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes.
1027 hmin = max(h_neglect, h_min_frac*((h(k-1) + h(k)) + (h(k+1) + h(k+2))))
1028 h0 = max(h(k-1), hmin) ; h1 = max(h(k), hmin)
1029 h2 = max(h(k+1), hmin) ; h3 = max(h(k+2), hmin)
1030
1031 ! Auxiliary calculations
1032 h1_2 = h1 * h1 ; h1_3 = h1_2 * h1 ; h1_4 = h1_2 * h1_2 ; h1_5 = h1_3 * h1_2
1033 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 asys(1:6,1) = (/ 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 /)
1040 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 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 -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 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 (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 bsys(1:6) = (/ 0.0, -2.0, 0.0, 0.0, 0.0, 0.0 /)
1052
1053 call linear_solver( 6, asys, bsys, csys )
1054
1055 alpha = csys(1)
1056 beta = csys(2)
1057
1058 tri_l(k+1) = alpha
1059 tri_d(k+1) = 1.0
1060 tri_u(k+1) = beta
1061 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 hmin = max(h_neglect, h_min_frac*((h(1) + h(2)) + (h(3) + h(4))))
1069 h0 = max(h(1), hmin) ; h1 = max(h(2), hmin)
1070 h2 = max(h(3), hmin) ; h3 = max(h(4), hmin)
1071
1072 ! Auxiliary calculations
1073 h1_2 = h1 * h1 ; h1_3 = h1_2 * h1 ; h1_4 = h1_2 * h1_2 ; h1_5 = h1_3 * h1_2
1074 h2_2 = h2 * h2 ; h2_3 = h2_2 * h2 ; h2_4 = h2_2 * h2_2 ; h2_5 = h2_3 * h2_2
1075 h01 = h0 + h1 ; h01_2 = h01 * h01
1076
1077 ! Compute matrix entries
1078 asys(1:6,1) = (/ 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 /)
1079 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 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 -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 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 (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 bsys(1:6) = (/ 0.0, -2.0, -6.0*h1, 12.0*h1_2, -20.0*h1_3, 30.0*h1_4 /)
1091
1092 call linear_solver( 6, asys, bsys, csys )
1093
1094 alpha = csys(1)
1095 beta = csys(2)
1096
1097 tri_l(2) = alpha
1098 tri_d(2) = 1.0
1099 tri_u(2) = beta
1100 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 x(1) = 0.0
1104 do i = 1,6
1105 dx = h(i)
1106 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 xavg * (xavg**4 + c5_6*xavg**2*dx**2 + 0.0625*dx**4) /)
1110 bsys(i) = u(i)
1111 x(i+1) = x(i) + dx
1112 enddo
1113
1114 call linear_solver( 6, asys, bsys, csys )
1115
1116 tri_d(1) = 0.0
1117 tri_d(1) = 1.0
1118 tri_u(1) = 0.0
1119 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 hmin = max(h_neglect, h_min_frac*((h(n-3) + h(n-2)) + (h(n-1) + h(n))))
1125 h0 = max(h(n-3), hmin) ; h1 = max(h(n-2), hmin)
1126 h2 = max(h(n-1), hmin) ; h3 = max(h(n), hmin)
1127
1128 ! Auxiliary calculations
1129 h1_2 = h1 * h1 ; h1_3 = h1_2 * h1 ; h1_4 = h1_2 * h1_2 ; h1_5 = h1_3 * h1_2
1130 h2_2 = h2 * h2 ; h2_3 = h2_2 * h2 ; h2_4 = h2_2 * h2_2 ; h2_5 = h2_3 * h2_2
1131
1132 h23 = h2 + h3 ; h23_2 = h23 * h23
1133
1134 ! Compute matrix entries
1135 asys(1:6,1) = (/ 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 /)
1136 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 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 -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 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 (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 bsys(1:6) = (/ 0.0, -2.0, 6.0*h2, 12.0*h2_2, 20.0*h2_3, 30.0*h2_4 /)
1148
1149 call linear_solver( 6, asys, bsys, csys )
1150
1151 alpha = csys(1)
1152 beta = csys(2)
1153
1154 tri_l(n) = alpha
1155 tri_d(n) = 1.0
1156 tri_u(n) = beta
1157 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 x(1) = 0.0
1161 do i = 1,6
1162 dx = h(n+1-i)
1163 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 xavg * (xavg**4 + c5_6*xavg**2*dx**2 + 0.0625*dx**4) /)
1167 bsys(i) = u(n+1-i)
1168 x(i+1) = x(i) + dx
1169 enddo
1170
1171 call linear_solver( 6, asys, bsys, csys )
1172
1173 tri_l(n+1) = 0.0
1174 tri_d(n+1) = 1.0
1175 tri_u(n+1) = 0.0
1176 tri_b(n+1) = -csys(2)
1177
1178 ! Solve tridiagonal system and assign edge values
1179 call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, n+1 )
1180
1181 do i = 2,n
1182 edge_slopes(i,1) = tri_x(i)
1183 edge_slopes(i-1,2) = tri_x(i)
1184 enddo
1185 edge_slopes(1,1) = tri_x(1)
1186 edge_slopes(n,2) = tri_x(n+1)
1187
1188end 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.
1226subroutine 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 real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal) [nondim]
1256 tri_d, & ! trid. system (middle diagonal) [nondim]
1257 tri_u, & ! trid. system (upper diagonal) [nondim]
1258 tri_b, & ! trid. system (rhs) [A]
1259 tri_x ! trid. system (unknowns vector) [A]
1260 integer :: i, k ! loop indexes
1261
1262 ! Loop on interior cells
1263 do k = 2,n-2
1264 ! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes.
1265 hmin = max(h_neglect, hminfrac*((h(k-1) + h(k)) + (h(k+1) + h(k+2))))
1266 h0 = max(h(k-1), hmin) ; h1 = max(h(k), hmin)
1267 h2 = max(h(k+1), hmin) ; h3 = max(h(k+2), hmin)
1268
1269 ! Auxiliary calculations
1270 h1_2 = h1 * h1 ; h1_3 = h1_2 * h1 ; h1_4 = h1_2 * h1_2 ; h1_5 = h1_3 * h1_2
1271 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 asys(1:6,1) = (/ 1.0, 1.0, -1.0, -1.0, -1.0, -1.0 /)
1275 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 -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 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 -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 -(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 bsys(1:6) = (/ -1.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
1287
1288 call linear_solver( 6, asys, bsys, csys )
1289
1290 alpha = csys(1)
1291 beta = csys(2)
1292
1293 tri_l(k+1) = alpha
1294 tri_d(k+1) = 1.0
1295 tri_u(k+1) = beta
1296 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 hmin = max(h_neglect, hminfrac*((h(1) + h(2)) + (h(3) + h(4))))
1304 h0 = max(h(1), hmin) ; h1 = max(h(2), hmin)
1305 h2 = max(h(3), hmin) ; h3 = max(h(4), hmin)
1306
1307 ! Auxiliary calculations
1308 h1_2 = h1 * h1 ; h1_3 = h1_2 * h1 ; h1_4 = h1_2 * h1_2 ; h1_5 = h1_3 * h1_2
1309 h2_2 = h2 * h2 ; h2_3 = h2_2 * h2 ; h2_4 = h2_2 * h2_2 ; h2_5 = h2_3 * h2_2
1310 h01 = h0 + h1 ; h01_2 = h01 * h01 ; h01_3 = h01 * h01_2
1311
1312 ! Compute matrix entries
1313 asys(1:6,1) = (/ 1.0, 1.0, -1.0, -1.0, -1.0, -1.0 /)
1314 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 -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 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 -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 -(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 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 call linear_solver( 6, asys, bsys, csys )
1328
1329 alpha = csys(1)
1330 beta = csys(2)
1331
1332 tri_l(2) = alpha
1333 tri_d(2) = 1.0
1334 tri_u(2) = beta
1335 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 hmin = max( h_neglect, hminfrac*((h(1)+h(2)) + (h(5)+h(6)) + (h(3)+h(4))) )
1339 x(1) = 0.0
1340 do i = 1,6
1341 dx = max( hmin, h(i) )
1342 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 xavg * (xavg**4 + c5_6*xavg**2*dx**2 + 0.0625*dx**4) /)
1346 bsys(i) = u(i)
1347 x(i+1) = x(i) + dx
1348 enddo
1349
1350 call linear_solver( 6, asys, bsys, csys )
1351
1352 tri_l(1) = 0.0
1353 tri_d(1) = 1.0
1354 tri_u(1) = 0.0
1355 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 hmin = max(h_neglect, hminfrac*((h(n-3) + h(n-2)) + (h(n-1) + h(n))))
1361 h0 = max(h(n-3), hmin) ; h1 = max(h(n-2), hmin)
1362 h2 = max(h(n-1), hmin) ; h3 = max(h(n), hmin)
1363
1364 ! Auxiliary calculations
1365 h1_2 = h1 * h1 ; h1_3 = h1_2 * h1 ; h1_4 = h1_2 * h1_2 ; h1_5 = h1_3 * h1_2
1366 h2_2 = h2 * h2 ; h2_3 = h2_2 * h2 ; h2_4 = h2_2 * h2_2 ; h2_5 = h2_3 * h2_2
1367 h23 = h2 + h3 ; h23_2 = h23 * h23 ; h23_3 = h23 * h23_2
1368
1369 ! Compute matrix entries
1370 asys(1:6,1) = (/ 1.0, 1.0, -1.0, -1.0, -1.0, -1.0 /)
1371 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 -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 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 -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 -(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 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 call linear_solver( 6, asys, bsys, csys )
1385
1386 alpha = csys(1)
1387 beta = csys(2)
1388
1389 tri_l(n) = alpha
1390 tri_d(n) = 1.0
1391 tri_u(n) = beta
1392 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 hmin = max( h_neglect, hminfrac*(h(n-3) + h(n-2)) + ((h(n-1) + h(n)) + (h(n-5) + h(n-4))) )
1396 x(1) = 0.0
1397 do i = 1,6
1398 dx = max( hmin, h(n+1-i) )
1399 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 xavg * (xavg**4 + c5_6*xavg**2*dx**2 + 0.0625*dx**4) /)
1403 bsys(i) = u(n+1-i)
1404 x(i+1) = x(i) + dx
1405 enddo
1406
1407 call linear_solver( 6, asys, bsys, csys )
1408
1409 tri_l(n+1) = 0.0
1410 tri_d(n+1) = 1.0
1411 tri_u(n+1) = 0.0
1412 tri_b(n+1) = csys(1)
1413
1414 ! Solve tridiagonal system and assign edge values
1415 call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, n+1 )
1416
1417 do i = 2,n
1418 edge_val(i,1) = tri_x(i)
1419 edge_val(i-1,2) = tri_x(i)
1420 enddo
1421 edge_val(1,1) = tri_x(1)
1422 edge_val(n,2) = tri_x(n+1)
1423
1424end 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.
1428subroutine 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 tolerance = 1.0e-12 ; if (present(tol)) tolerance = tol
1443
1444 sum = 0.0 ; sum_mag = max(0.0,mag)
1445 do i=1,n
1446 sum = sum + a(i) * c(i)
1447 sum_mag = sum_mag + abs(a(i) * c(i))
1448 enddo
1449
1450 if (abs(sum - r) > tolerance * (sum_mag + abs(r))) then
1451 write(mesg2, '(", Fractional error = ", es12.4,", sum = ", es12.4)') (sum - r) / (sum_mag + abs(r)), sum
1452 call mom_error(fatal, "Failed line test: "//trim(msg)//trim(mesg2))
1453 endif
1454
1455end subroutine test_line
1456
1457end module regrid_edge_values