MOM6
Functions/Subroutines | Variables
regrid_edge_slopes Module Reference

Detailed Description

Routines that estimate edge slopes to be used in high-order reconstruction schemes.

Functions/Subroutines

subroutine, public edge_slopes_implicit_h3 (N, h, u, edge_slopes, h_neglect)
 Compute ih4 edge slopes (implicit third order accurate) in the same units as h. More...
 
subroutine, public edge_slopes_implicit_h5 (N, h, u, edge_slopes, h_neglect)
 Compute ih5 edge values (implicit fifth order accurate) More...
 

Variables

real, parameter hneglect_dflt = 1.E-30
 Default negligible cell thickness.
 

Function/Subroutine Documentation

◆ edge_slopes_implicit_h3()

subroutine, public regrid_edge_slopes::edge_slopes_implicit_h3 ( integer, intent(in)  N,
real, dimension(:), intent(in)  h,
real, dimension(:), intent(in)  u,
real, dimension(:,:), intent(inout)  edge_slopes,
real, intent(in), optional  h_neglect 
)

Compute ih4 edge slopes (implicit third order accurate) in the same units as h.

Compute edge slopes based on third-order implicit estimates. Note that the estimates are fourth-order accurate on uniform grids

Third-order implicit estimates of edge slopes are based on a two-cell stencil. A tridiagonal system is set up and is based on expressing the edge slopes in terms of neighboring cell averages. The generic relationship is

\[ \alpha u'_{i-1/2} + u'_{i+1/2} + \beta u'_{i+3/2} = a \bar{u}_i + b \bar{u}_{i+1} \]

and the stencil looks like this 

     i     i+1

..–o---—o---—o–.. i-1/2 i+1/2 i+3/2

In this routine, the coefficients \(\alpha\), \(\beta\), a and b are computed, the tridiagonal system is built, boundary conditions are prescribed and the system is solved to yield edge-slope estimates.

There are N+1 unknowns and we are able to write N-1 equations. The boundary conditions close the system.

Parameters
[in]nNumber of cells
[in]hcell widths (size N)
[in]ucell average properties (size N)
[in,out]edge_slopesReturned edge slopes, with the same units as u divided by the units of h.
[in]h_neglectA negligibly small width

Definition at line 50 of file regrid_edge_slopes.F90.

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 

◆ edge_slopes_implicit_h5()

subroutine, public regrid_edge_slopes::edge_slopes_implicit_h5 ( integer, intent(in)  N,
real, dimension(:), intent(in)  h,
real, dimension(:), intent(in)  u,
real, dimension(:,:), intent(inout)  edge_slopes,
real, intent(in), optional  h_neglect 
)

Compute ih5 edge values (implicit fifth order accurate)

Parameters
[in]nNumber of cells
[in]hcell widths (size N)
[in]ucell average properties (size N)
[in,out]edge_slopesReturned edge slopes, with the same units as u divided by the units of h.
[in]h_neglectA negligibly small width in the same units as h.

Definition at line 177 of file regrid_edge_slopes.F90.

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