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