MOM6
regrid_edge_slopes.F90
Go to the documentation of this file.
1 !> Routines that estimate edge slopes to be used in
2 !! high-order reconstruction schemes.
3 module regrid_edge_slopes
4 
5 ! This file is part of MOM6. See LICENSE.md for the license.
6 
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 public edge_slopes_implicit_h3
13 public edge_slopes_implicit_h5
14 
15 ! Specifying a dimensional parameter value, as is done here, is a terrible idea.
16 real, parameter :: hneglect_dflt = 1.e-30 !< Default negligible cell thickness
17 
18 contains
19 
20 !------------------------------------------------------------------------------
21 !> Compute ih4 edge slopes (implicit third order accurate)
22 !! in the same units as h.
23 !!
24 !! Compute edge slopes based on third-order implicit estimates. Note that
25 !! the estimates are fourth-order accurate on uniform grids
26 !!
27 !! Third-order implicit estimates of edge slopes are based on a two-cell
28 !! stencil. A tridiagonal system is set up and is based on expressing the
29 !! edge slopes in terms of neighboring cell averages. The generic
30 !! relationship is
31 !!
32 !! \f[
33 !! \alpha u'_{i-1/2} + u'_{i+1/2} + \beta u'_{i+3/2} =
34 !! a \bar{u}_i + b \bar{u}_{i+1}
35 !! \f]
36 !!
37 !! and the stencil looks like this
38 !!
39 !! i i+1
40 !! ..--o------o------o--..
41 !! i-1/2 i+1/2 i+3/2
42 !!
43 !! In this routine, the coefficients \f$\alpha\f$, \f$\beta\f$, a and b are computed,
44 !! the tridiagonal system is built, boundary conditions are prescribed and
45 !! the system is solved to yield edge-slope estimates.
46 !!
47 !! There are N+1 unknowns and we are able to write N-1 equations. The
48 !! boundary conditions close the system.
49 subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect, answers_2018 )
50  integer, intent(in) :: n !< Number of cells
51  real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
52  real, dimension(:), intent(in) :: u !< cell average properties (size N) in arbitrary units [A]
53  real, dimension(:,:), intent(inout) :: edge_slopes !< Returned edge slopes [A H-1]
54  real, optional, intent(in) :: h_neglect !< A negligibly small width
55  logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions.
56  ! Local variables
57  integer :: i, j ! loop indexes
58  real :: h0, h1 ! cell widths [H]
59  real :: h0_2, h1_2, h0h1 ! products of cell widths [H2]
60  real :: h0_3, h1_3 ! products of three cell widths [H3]
61  real :: d ! A demporary variable [H3]
62  real :: alpha, beta ! stencil coefficients [nondim]
63  real :: a, b ! weights of cells [H-1]
64  real, parameter :: c1_12 = 1.0 / 12.0
65  real, dimension(5) :: x ! Coordinate system with 0 at edges [H]
66  real :: dx, xavg ! Differences and averages of successive values of x [H]
67  real, dimension(4,4) :: asys ! matrix used to find boundary conditions
68  real, dimension(4) :: bsys, csys
69  real, dimension(3) :: dsys
70  real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal) [nondim]
71  tri_d, & ! trid. system (middle diagonal) [nondim]
72  tri_u, & ! trid. system (upper diagonal) [nondim]
73  tri_b, & ! trid. system (unknowns vector) [A H-1]
74  tri_x ! trid. system (rhs) [A H-1]
75  real :: hneglect ! A negligible thickness [H].
76  real :: hneglect3 ! hNeglect^3 [H3].
77  logical :: use_2018_answers ! If true use older, less acccurate expressions.
78 
79  hneglect = hneglect_dflt ; if (present(h_neglect)) hneglect = h_neglect
80  hneglect3 = hneglect**3
81  use_2018_answers = .true. ; if (present(answers_2018)) use_2018_answers = answers_2018
82 
83  ! Loop on cells (except last one)
84  do i = 1,n-1
85 
86  ! Get cell widths
87  h0 = h(i)
88  h1 = h(i+1)
89 
90  ! Auxiliary calculations
91  h0h1 = h0 * h1
92  h0_2 = h0 * h0
93  h1_2 = h1 * h1
94  h0_3 = h0_2 * h0
95  h1_3 = h1_2 * h1
96 
97  d = 4.0 * h0h1 * ( h0 + h1 ) + h1_3 + h0_3
98 
99  ! Coefficients
100  alpha = h1 * (h0_2 + h0h1 - h1_2) / ( d + hneglect3 )
101  beta = h0 * (h1_2 + h0h1 - h0_2) / ( d + hneglect3 )
102  a = -12.0 * h0h1 / ( d + hneglect3 )
103  b = -a
104 
105  tri_l(i+1) = alpha
106  tri_d(i+1) = 1.0
107  tri_u(i+1) = beta
108 
109  tri_b(i+1) = a * u(i) + b * u(i+1)
110 
111  enddo ! end loop on cells
112 
113  ! Boundary conditions: left boundary
114  x(1) = 0.0
115  do i = 2,5
116  x(i) = x(i-1) + h(i-1)
117  enddo
118 
119  do i = 1,4
120  dx = h(i)
121  if (use_2018_answers) then
122  do j = 1,4 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo
123  else ! Use expressions with less sensitivity to roundoff
124  xavg = 0.5 * (x(i+1) + x(i))
125  asys(i,1) = dx
126  asys(i,2) = dx * xavg
127  asys(i,3) = dx * (xavg**2 + c1_12*dx**2)
128  asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2)
129  endif
130 
131  bsys(i) = u(i) * dx
132 
133  enddo
134 
135  call solve_linear_system( asys, bsys, csys, 4 )
136 
137  dsys(1) = csys(2)
138  dsys(2) = 2.0 * csys(3)
139  dsys(3) = 3.0 * csys(4)
140 
141  tri_d(1) = 1.0
142  tri_u(1) = 0.0
143  tri_b(1) = evaluation_polynomial( dsys, 3, x(1) ) ! first edge slope
144 
145  ! Boundary conditions: right boundary
146  x(1) = 0.0
147  do i = 2,5
148  x(i) = x(i-1) + h(n-5+i)
149  enddo
150 
151  do i = 1,4
152  dx = h(n-4+i)
153  if (use_2018_answers) then
154  do j = 1,4 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo
155  else ! Use expressions with less sensitivity to roundoff
156  xavg = 0.5 * (x(i+1) + x(i))
157  asys(i,1) = dx
158  asys(i,2) = dx * xavg
159  asys(i,3) = dx * (xavg**2 + c1_12*dx**2)
160  asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2)
161  endif
162  bsys(i) = u(n-4+i) * dx
163 
164  enddo
165 
166  call solve_linear_system( asys, bsys, csys, 4 )
167 
168  dsys(1) = csys(2)
169  dsys(2) = 2.0 * csys(3)
170  dsys(3) = 3.0 * csys(4)
171 
172  tri_l(n+1) = 0.0
173  tri_d(n+1) = 1.0
174  tri_b(n+1) = evaluation_polynomial( dsys, 3, x(5) ) ! last edge slope
175 
176  ! Solve tridiagonal system and assign edge values
177  call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, n+1 )
178 
179  do i = 2,n
180  edge_slopes(i,1) = tri_x(i)
181  edge_slopes(i-1,2) = tri_x(i)
182  enddo
183  edge_slopes(1,1) = tri_x(1)
184  edge_slopes(n,2) = tri_x(n+1)
185 
186 end subroutine edge_slopes_implicit_h3
187 
188 
189 !------------------------------------------------------------------------------
190 !> Compute ih5 edge values (implicit fifth order accurate)
191 subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect, answers_2018 )
192  integer, intent(in) :: n !< Number of cells
193  real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
194  real, dimension(:), intent(in) :: u !< cell average properties (size N) in arbitrary units [A]
195  real, dimension(:,:), intent(inout) :: edge_slopes !< Returned edge slopes [A H-1]
196  real, optional, intent(in) :: h_neglect !< A negligibly small width [H]
197  logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions.
198 ! -----------------------------------------------------------------------------
199 ! Fifth-order implicit estimates of edge values are based on a four-cell,
200 ! three-edge stencil. A tridiagonal system is set up and is based on
201 ! expressing the edge slopes in terms of neighboring cell averages.
202 !
203 ! The generic relationship is
204 !
205 ! \alpha u'_{i-1/2} + u'_{i+1/2} + \beta u'_{i+3/2} =
206 ! a \bar{u}_{i-1} + b \bar{u}_i + c \bar{u}_{i+1} + d \bar{u}_{i+2}
207 !
208 ! and the stencil looks like this
209 !
210 ! i-1 i i+1 i+2
211 ! ..--o------o------o------o------o--..
212 ! i-1/2 i+1/2 i+3/2
213 !
214 ! In this routine, the coefficients \alpha, \beta, a, b, c and d are
215 ! computed, the tridiagonal system is built, boundary conditions are
216 ! prescribed and the system is solved to yield edge-value estimates.
217 !
218 ! Note that the centered stencil only applies to edges 3 to N-1 (edges are
219 ! numbered 1 to n+1), which yields N-3 equations for N+1 unknowns. Two other
220 ! equations are written by using a right-biased stencil for edge 2 and a
221 ! left-biased stencil for edge N. The prescription of boundary conditions
222 ! (using sixth-order polynomials) closes the system.
223 !
224 ! CAUTION: For each edge, in order to determine the coefficients of the
225 ! implicit expression, a 6x6 linear system is solved. This may
226 ! become computationally expensive if regridding is carried out
227 ! often. Figuring out closed-form expressions for these coefficients
228 ! on nonuniform meshes turned out to be intractable.
229 ! -----------------------------------------------------------------------------
230 
231  ! Local variables
232  integer :: i, j, k ! loop indexes
233  real :: h0, h1, h2, h3 ! cell widths
234  real :: g, g_2, g_3 ! the following are
235  real :: g_4, g_5, g_6 ! auxiliary variables
236  real :: d2, d3, d4, d5, d6 ! to set up the systems
237  real :: n2, n3, n4, n5, n6 ! used to compute the
238  real :: h1_2, h2_2 ! the coefficients of the
239  real :: h1_3, h2_3 ! tridiagonal system
240  real :: h1_4, h2_4 ! ...
241  real :: h1_5, h2_5 ! ...
242  real :: h1_6, h2_6 ! ...
243  real :: h0ph1, h0ph1_2 ! ...
244  real :: h0ph1_3, h0ph1_4 ! ...
245  real :: h2ph3, h2ph3_2 ! ...
246  real :: h2ph3_3, h2ph3_4 ! ...
247  real :: alpha, beta ! stencil coefficients
248  real :: a, b, c, d ! "
249  real, dimension(7) :: x ! Coordinate system with 0 at edges [same units as h]
250  real, parameter :: c1_12 = 1.0 / 12.0
251  real, parameter :: c5_6 = 5.0 / 6.0
252  real :: dx, xavg ! Differences and averages of successive values of x [same units as h]
253  real, dimension(6,6) :: asys ! matrix used to find boundary conditions
254  real, dimension(6) :: bsys, csys ! ...
255  real, dimension(5) :: dsys ! derivative
256  real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal)
257  tri_d, & ! trid. system (middle diagonal)
258  tri_u, & ! trid. system (upper diagonal)
259  tri_b, & ! trid. system (unknowns vector)
260  tri_x ! trid. system (rhs)
261  real :: hneglect ! A negligible thickness in the same units as h.
262  logical :: use_2018_answers ! If true use older, less acccurate expressions.
263 
264  hneglect = hneglect_dflt ; if (present(h_neglect)) hneglect = h_neglect
265  use_2018_answers = .true. ; if (present(answers_2018)) use_2018_answers = answers_2018
266 
267  ! Loop on cells (except last one)
268  do k = 2,n-2
269 
270  ! Cell widths
271  h0 = h(k-1)
272  h1 = h(k+0)
273  h2 = h(k+1)
274  h3 = h(k+2)
275 
276  ! Auxiliary calculations
277  h1_2 = h1 * h1
278  h1_3 = h1_2 * h1
279  h1_4 = h1_2 * h1_2
280  h1_5 = h1_3 * h1_2
281  h1_6 = h1_3 * h1_3
282 
283  h2_2 = h2 * h2
284  h2_3 = h2_2 * h2
285  h2_4 = h2_2 * h2_2
286  h2_5 = h2_3 * h2_2
287  h2_6 = h2_3 * h2_3
288 
289  g = h0 + h1
290  g_2 = g * g
291  g_3 = g * g_2
292  g_4 = g_2 * g_2
293  g_5 = g_4 * g
294  g_6 = g_3 * g_3
295 
296  d2 = ( h1_2 - g_2 ) / ( h0 + hneglect )
297  d3 = ( h1_3 - g_3 ) / ( h0 + hneglect )
298  d4 = ( h1_4 - g_4 ) / ( h0 + hneglect )
299  d5 = ( h1_5 - g_5 ) / ( h0 + hneglect )
300  d6 = ( h1_6 - g_6 ) / ( h0 + hneglect )
301 
302  g = h2 + h3
303  g_2 = g * g
304  g_3 = g * g_2
305  g_4 = g_2 * g_2
306  g_5 = g_4 * g
307  g_6 = g_3 * g_3
308 
309  n2 = ( g_2 - h2_2 ) / ( h3 + hneglect )
310  n3 = ( g_3 - h2_3 ) / ( h3 + hneglect )
311  n4 = ( g_4 - h2_4 ) / ( h3 + hneglect )
312  n5 = ( g_5 - h2_5 ) / ( h3 + hneglect )
313  n6 = ( g_6 - h2_6 ) / ( h3 + hneglect )
314 
315  ! Compute matrix entries
316  asys(1,1) = 0.0
317  asys(1,2) = 0.0
318  asys(1,3) = 1.0
319  asys(1,4) = 1.0
320  asys(1,5) = 1.0
321  asys(1,6) = 1.0
322 
323  asys(2,1) = 1.0
324  asys(2,2) = 1.0
325  asys(2,3) = -0.5 * d2
326  asys(2,4) = 0.5 * h1
327  asys(2,5) = -0.5 * h2
328  asys(2,6) = -0.5 * n2
329 
330  asys(3,1) = h1
331  asys(3,2) = - h2
332  asys(3,3) = - d3 / 6.0
333  asys(3,4) = h1_2 / 6.0
334  asys(3,5) = h2_2 / 6.0
335  asys(3,6) = n3 / 6.0
336 
337  asys(4,1) = - h1_2 / 2.0
338  asys(4,2) = - h2_2 / 2.0
339  asys(4,3) = d4 / 24.0
340  asys(4,4) = - h1_3 / 24.0
341  asys(4,5) = h2_3 / 24.0
342  asys(4,6) = n4 / 24.0
343 
344  asys(5,1) = h1_3 / 6.0
345  asys(5,2) = - h2_3 / 6.0
346  asys(5,3) = - d5 / 120.0
347  asys(5,4) = h1_4 / 120.0
348  asys(5,5) = h2_4 / 120.0
349  asys(5,6) = n5 / 120.0
350 
351  asys(6,1) = - h1_4 / 24.0
352  asys(6,2) = - h2_4 / 24.0
353  asys(6,3) = d6 / 720.0
354  asys(6,4) = - h1_5 / 720.0
355  asys(6,5) = h2_5 / 720.0
356  asys(6,6) = n6 / 720.0
357 
358  bsys(:) = (/ 0.0, -1.0, 0.0, 0.0, 0.0, 0.0 /)
359 
360  call solve_linear_system( asys, bsys, csys, 6 )
361 
362  alpha = csys(1)
363  beta = csys(2)
364  a = csys(3)
365  b = csys(4)
366  c = csys(5)
367  d = csys(6)
368 
369  tri_l(k+1) = alpha
370  tri_d(k+1) = 1.0
371  tri_u(k+1) = beta
372  tri_b(k+1) = a * u(k-1) + b * u(k) + c * u(k+1) + d * u(k+2)
373 
374  enddo ! end loop on cells
375 
376  ! Use a right-biased stencil for the second row
377 
378  ! Cell widths
379  h0 = h(1)
380  h1 = h(2)
381  h2 = h(3)
382  h3 = h(4)
383 
384  ! Auxiliary calculations
385  h1_2 = h1 * h1
386  h1_3 = h1_2 * h1
387  h1_4 = h1_2 * h1_2
388  h1_5 = h1_3 * h1_2
389  h1_6 = h1_3 * h1_3
390 
391  h2_2 = h2 * h2
392  h2_3 = h2_2 * h2
393  h2_4 = h2_2 * h2_2
394  h2_5 = h2_3 * h2_2
395  h2_6 = h2_3 * h2_3
396 
397  g = h0 + h1
398  g_2 = g * g
399  g_3 = g * g_2
400  g_4 = g_2 * g_2
401  g_5 = g_4 * g
402  g_6 = g_3 * g_3
403 
404  h0ph1 = h0 + h1
405  h0ph1_2 = h0ph1 * h0ph1
406  h0ph1_3 = h0ph1_2 * h0ph1
407  h0ph1_4 = h0ph1_2 * h0ph1_2
408 
409  d2 = ( h1_2 - g_2 ) / ( h0 + hneglect )
410  d3 = ( h1_3 - g_3 ) / ( h0 + hneglect )
411  d4 = ( h1_4 - g_4 ) / ( h0 + hneglect )
412  d5 = ( h1_5 - g_5 ) / ( h0 + hneglect )
413  d6 = ( h1_6 - g_6 ) / ( h0 + hneglect )
414 
415  g = h2 + h3
416  g_2 = g * g
417  g_3 = g * g_2
418  g_4 = g_2 * g_2
419  g_5 = g_4 * g
420  g_6 = g_3 * g_3
421 
422  n2 = ( g_2 - h2_2 ) / ( h3 + hneglect )
423  n3 = ( g_3 - h2_3 ) / ( h3 + hneglect )
424  n4 = ( g_4 - h2_4 ) / ( h3 + hneglect )
425  n5 = ( g_5 - h2_5 ) / ( h3 + hneglect )
426  n6 = ( g_6 - h2_6 ) / ( h3 + hneglect )
427 
428  ! Compute matrix entries
429  asys(1,1) = 0.0
430  asys(1,2) = 0.0
431  asys(1,3) = 1.0
432  asys(1,4) = 1.0
433  asys(1,5) = 1.0
434  asys(1,6) = 1.0
435 
436  asys(2,1) = 1.0
437  asys(2,2) = 1.0
438  asys(2,3) = -0.5 * d2
439  asys(2,4) = 0.5 * h1
440  asys(2,5) = -0.5 * h2
441  asys(2,6) = -0.5 * n2
442 
443  asys(3,1) = h0ph1
444  asys(3,2) = 0.0
445  asys(3,3) = - d3 / 6.0
446  asys(3,4) = h1_2 / 6.0
447  asys(3,5) = h2_2 / 6.0
448  asys(3,6) = n3 / 6.0
449 
450  asys(4,1) = - h0ph1_2 / 2.0
451  asys(4,2) = 0.0
452  asys(4,3) = d4 / 24.0
453  asys(4,4) = - h1_3 / 24.0
454  asys(4,5) = h2_3 / 24.0
455  asys(4,6) = n4 / 24.0
456 
457  asys(5,1) = h0ph1_3 / 6.0
458  asys(5,2) = 0.0
459  asys(5,3) = - d5 / 120.0
460  asys(5,4) = h1_4 / 120.0
461  asys(5,5) = h2_4 / 120.0
462  asys(5,6) = n5 / 120.0
463 
464  asys(6,1) = - h0ph1_4 / 24.0
465  asys(6,2) = 0.0
466  asys(6,3) = d6 / 720.0
467  asys(6,4) = - h1_5 / 720.0
468  asys(6,5) = h2_5 / 720.0
469  asys(6,6) = n6 / 720.0
470 
471  bsys(:) = (/ 0.0, -1.0, -h1, h1_2/2.0, -h1_3/6.0, h1_4/24.0 /)
472 
473  call solve_linear_system( asys, bsys, csys, 6 )
474 
475  alpha = csys(1)
476  beta = csys(2)
477  a = csys(3)
478  b = csys(4)
479  c = csys(5)
480  d = csys(6)
481 
482  tri_l(2) = alpha
483  tri_d(2) = 1.0
484  tri_u(2) = beta
485  tri_b(2) = a * u(1) + b * u(2) + c * u(3) + d * u(4)
486 
487  ! Boundary conditions: left boundary
488  x(1) = 0.0
489  do i = 2,7
490  x(i) = x(i-1) + h(i-1)
491  enddo
492 
493  do i = 1,6
494 
495  dx = h(i)
496  if (use_2018_answers) then
497  do j = 1,6 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo
498  else ! Use expressions with less sensitivity to roundoff
499  xavg = 0.5 * (x(i+1) + x(i))
500  asys(i,1) = dx
501  asys(i,2) = dx * xavg
502  asys(i,3) = dx * (xavg**2 + c1_12*dx**2)
503  asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2)
504  asys(i,5) = dx * (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4)
505  asys(i,6) = dx * xavg * (xavg**4 + c5_6*xavg**2*dx**2 + 0.0625*dx**4)
506  endif
507 
508  bsys(i) = u(i) * dx
509 
510  enddo
511 
512  call solve_linear_system( asys, bsys, csys, 6 )
513 
514  dsys(1) = csys(2)
515  dsys(2) = 2.0 * csys(3)
516  dsys(3) = 3.0 * csys(4)
517  dsys(4) = 4.0 * csys(5)
518  dsys(5) = 5.0 * csys(6)
519 
520  tri_d(1) = 0.0
521  tri_d(1) = 1.0
522  tri_u(1) = 0.0
523  tri_b(1) = evaluation_polynomial( dsys, 5, x(1) ) ! first edge value
524 
525  ! Use a left-biased stencil for the second to last row
526 
527  ! Cell widths
528  h0 = h(n-3)
529  h1 = h(n-2)
530  h2 = h(n-1)
531  h3 = h(n)
532 
533  ! Auxiliary calculations
534  h1_2 = h1 * h1
535  h1_3 = h1_2 * h1
536  h1_4 = h1_2 * h1_2
537  h1_5 = h1_3 * h1_2
538  h1_6 = h1_3 * h1_3
539 
540  h2_2 = h2 * h2
541  h2_3 = h2_2 * h2
542  h2_4 = h2_2 * h2_2
543  h2_5 = h2_3 * h2_2
544  h2_6 = h2_3 * h2_3
545 
546  g = h0 + h1
547  g_2 = g * g
548  g_3 = g * g_2
549  g_4 = g_2 * g_2
550  g_5 = g_4 * g
551  g_6 = g_3 * g_3
552 
553  h2ph3 = h2 + h3
554  h2ph3_2 = h2ph3 * h2ph3
555  h2ph3_3 = h2ph3_2 * h2ph3
556  h2ph3_4 = h2ph3_2 * h2ph3_2
557 
558  d2 = ( h1_2 - g_2 ) / ( h0 + hneglect )
559  d3 = ( h1_3 - g_3 ) / ( h0 + hneglect )
560  d4 = ( h1_4 - g_4 ) / ( h0 + hneglect )
561  d5 = ( h1_5 - g_5 ) / ( h0 + hneglect )
562  d6 = ( h1_6 - g_6 ) / ( h0 + hneglect )
563 
564  g = h2 + h3
565  g_2 = g * g
566  g_3 = g * g_2
567  g_4 = g_2 * g_2
568  g_5 = g_4 * g
569  g_6 = g_3 * g_3
570 
571  n2 = ( g_2 - h2_2 ) / ( h3 + hneglect )
572  n3 = ( g_3 - h2_3 ) / ( h3 + hneglect )
573  n4 = ( g_4 - h2_4 ) / ( h3 + hneglect )
574  n5 = ( g_5 - h2_5 ) / ( h3 + hneglect )
575  n6 = ( g_6 - h2_6 ) / ( h3 + hneglect )
576 
577  ! Compute matrix entries
578  asys(1,1) = 0.0
579  asys(1,2) = 0.0
580  asys(1,3) = 1.0
581  asys(1,4) = 1.0
582  asys(1,5) = 1.0
583  asys(1,6) = 1.0
584 
585  asys(2,1) = 1.0
586  asys(2,2) = 1.0
587  asys(2,3) = -0.5 * d2
588  asys(2,4) = 0.5 * h1
589  asys(2,5) = -0.5 * h2
590  asys(2,6) = -0.5 * n2
591 
592  asys(3,1) = 0.0
593  asys(3,2) = - h2ph3
594  asys(3,3) = - d3 / 6.0
595  asys(3,4) = h1_2 / 6.0
596  asys(3,5) = h2_2 / 6.0
597  asys(3,6) = n3 / 6.0
598 
599  asys(4,1) = 0.0
600  asys(4,2) = - h2ph3_2 / 2.0
601  asys(4,3) = d4 / 24.0
602  asys(4,4) = - h1_3 / 24.0
603  asys(4,5) = h2_3 / 24.0
604  asys(4,6) = n4 / 24.0
605 
606  asys(5,1) = 0.0
607  asys(5,2) = - h2ph3_3 / 6.0
608  asys(5,3) = - d5 / 120.0
609  asys(5,4) = h1_4 / 120.0
610  asys(5,5) = h2_4 / 120.0
611  asys(5,6) = n5 / 120.0
612 
613  asys(6,1) = 0.0
614  asys(6,2) = - h2ph3_4 / 24.0
615  asys(6,3) = d6 / 720.0
616  asys(6,4) = - h1_5 / 720.0
617  asys(6,5) = h2_5 / 720.0
618  asys(6,6) = n6 / 720.0
619 
620  bsys(:) = (/ 0.0, -1.0, h2, h2_2/2.0, h2_3/6.0, h2_4/24.0 /)
621 
622  call solve_linear_system( asys, bsys, csys, 6 )
623 
624  alpha = csys(1)
625  beta = csys(2)
626  a = csys(3)
627  b = csys(4)
628  c = csys(5)
629  d = csys(6)
630 
631  tri_l(n) = alpha
632  tri_d(n) = 1.0
633  tri_u(n) = beta
634  tri_b(n) = a * u(n-3) + b * u(n-2) + c * u(n-1) + d * u(n)
635 
636  ! Boundary conditions: right boundary
637  x(1) = 0.0
638  do i = 2,7
639  x(i) = x(i-1) + h(n-7+i)
640  enddo
641 
642  do i = 1,6
643  dx = h(n-6+i)
644  if (use_2018_answers) then
645  do j = 1,6 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo
646  else ! Use expressions with less sensitivity to roundoff
647  xavg = 0.5 * (x(i+1) + x(i))
648  asys(i,1) = dx
649  asys(i,2) = dx * xavg
650  asys(i,3) = dx * (xavg**2 + c1_12*dx**2)
651  asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2)
652  asys(i,5) = dx * (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4)
653  asys(i,6) = dx * xavg * (xavg**4 + c5_6*xavg**2*dx**2 + 0.0625*dx**4)
654  endif
655  bsys(i) = u(n-6+i) * dx
656  enddo
657 
658  call solve_linear_system( asys, bsys, csys, 6 )
659 
660  dsys(1) = csys(2)
661  dsys(2) = 2.0 * csys(3)
662  dsys(3) = 3.0 * csys(4)
663  dsys(4) = 4.0 * csys(5)
664  dsys(5) = 5.0 * csys(6)
665 
666  tri_l(n+1) = 0.0
667  tri_d(n+1) = 1.0
668  tri_u(n+1) = 0.0
669  tri_b(n+1) = evaluation_polynomial( dsys, 5, x(7) ) ! last edge value
670 
671  ! Solve tridiagonal system and assign edge values
672  call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, n+1 )
673 
674  do i = 2,n
675  edge_slopes(i,1) = tri_x(i)
676  edge_slopes(i-1,2) = tri_x(i)
677  enddo
678  edge_slopes(1,1) = tri_x(1)
679  edge_slopes(n,2) = tri_x(n+1)
680 
681 end subroutine edge_slopes_implicit_h5
682 
683 end module regrid_edge_slopes
regrid_edge_slopes::edge_slopes_implicit_h3
subroutine, public edge_slopes_implicit_h3(N, h, u, edge_slopes, h_neglect, answers_2018)
Compute ih4 edge slopes (implicit third order accurate) in the same units as h.
Definition: regrid_edge_slopes.F90:50
polynomial_functions::evaluation_polynomial
real function, public evaluation_polynomial(coeff, ncoef, x)
Pointwise evaluation of a polynomial at x.
Definition: polynomial_functions.F90:20
regrid_solvers
Solvers of linear systems.
Definition: regrid_solvers.F90:2
regrid_edge_slopes::edge_slopes_implicit_h5
subroutine, public edge_slopes_implicit_h5(N, h, u, edge_slopes, h_neglect, answers_2018)
Compute ih5 edge values (implicit fifth order accurate)
Definition: regrid_edge_slopes.F90:192
regrid_edge_slopes::hneglect_dflt
real, parameter hneglect_dflt
Default negligible cell thickness.
Definition: regrid_edge_slopes.F90:16
regrid_solvers::solve_linear_system
subroutine, public solve_linear_system(A, B, X, system_size)
Solve the linear system AX = B by Gaussian elimination.
Definition: regrid_solvers.F90:20
polynomial_functions
Polynomial functions.
Definition: polynomial_functions.F90:2
regrid_solvers::solve_tridiagonal_system
subroutine, public solve_tridiagonal_system(Al, Ad, Au, B, X, system_size)
Solve the tridiagonal system AX = B.
Definition: regrid_solvers.F90:114
regrid_edge_slopes
Routines that estimate edge slopes to be used in high-order reconstruction schemes.
Definition: regrid_edge_slopes.F90:3