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