MOM6
Functions/Subroutines | Variables
regrid_edge_slopes Module Reference

Functions/Subroutines

subroutine, public edge_slopes_implicit_h3 (N, h, u, edge_slopes)
 
subroutine, public edge_slopes_implicit_h5 (N, h, u, edge_slopes)
 

Variables

real, parameter h_neglect = 1.E-30
 

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 
)

Definition at line 34 of file regrid_edge_slopes.F90.

References polynomial_functions::evaluation_polynomial(), h_neglect, regrid_solvers::solve_linear_system(), and regrid_solvers::solve_tridiagonal_system().

Referenced by mom_remapping::build_reconstructions_1d(), and regrid_interp::regridding_set_ppolys().

34 ! -----------------------------------------------------------------------------
35 ! Compute edge slopes based on third-order implicit estimates. Note that
36 ! the estimates are fourth-order accurate on uniform grids
37 !
38 ! Third-order implicit estimates of edge slopes are based on a two-cell
39 ! stencil. A tridiagonal system is set up and is based on expressing the
40 ! edge slopes in terms of neighboring cell averages. The generic
41 ! relationship is
42 !
43 ! \alpha u'_{i-1/2} + u'_{i+1/2} + \beta u'_{i+3/2} =
44 ! a \bar{u}_i + b \bar{u}_{i+1}
45 !
46 ! and the stencil looks like this
47 !
48 ! i i+1
49 ! ..--o------o------o--..
50 ! i-1/2 i+1/2 i+3/2
51 !
52 ! In this routine, the coefficients \alpha, \beta, a and b are computed,
53 ! the tridiagonal system is built, boundary conditions are prescribed and
54 ! the system is solved to yield edge-slope estimates.
55 !
56 ! There are N+1 unknowns and we are able to write N-1 equations. The
57 ! boundary conditions close the system.
58 ! -----------------------------------------------------------------------------
59 
60  ! Arguments
61  integer, intent(in) :: n ! Number of cells
62  real, dimension(:), intent(in) :: h ! cell averages (size N)
63  real, dimension(:), intent(in) :: u ! cell averages (size N)
64  real, dimension(:,:), intent(inout) :: edge_slopes
65 
66  ! Local variables
67  integer :: i, j ! loop indexes
68  real :: h0, h1 ! cell widths
69  real :: h0_2, h1_2, h0h1
70  real :: h0_3, h1_3
71  real :: d
72  real :: alpha, beta ! stencil coefficients
73  real :: a, b
74  real, dimension(5) :: x ! system used to enforce
75  real, dimension(4,4) :: asys ! boundary conditions
76  real, dimension(4) :: bsys, csys
77  real, dimension(3) :: dsys
78  real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal)
79  tri_d, & ! trid. system (middle diagonal)
80  tri_u, & ! trid. system (upper diagonal)
81  tri_b, & ! trid. system (unknowns vector)
82  tri_x ! trid. system (rhs)
83 
84  ! Loop on cells (except last one)
85  do i = 1,n-1
86 
87  ! Get cell widths
88  h0 = h(i)
89  h1 = h(i+1)
90 
91  ! Auxiliary calculations
92  h0h1 = h0 * h1
93  h0_2 = h0 * h0
94  h1_2 = h1 * h1
95  h0_3 = h0_2 * h0
96  h1_3 = h1_2 * h1
97 
98  d = 4.0 * h0h1 * ( h0 + h1 ) + h1_3 + h0_3
99 
100  ! Coefficients
101  alpha = h1 * (h0_2 + h0h1 - h1_2) / ( d + h_neglect )
102  beta = h0 * (h1_2 + h0h1 - h0_2) / ( d + h_neglect )
103  a = -12.0 * h0h1 / ( d + h_neglect )
104  b = -a
105 
106  tri_l(i+1) = alpha
107  tri_d(i+1) = 1.0
108  tri_u(i+1) = beta
109 
110  tri_b(i+1) = a * u(i) + b * u(i+1)
111 
112  end do ! end loop on cells
113 
114  ! Boundary conditions: left boundary
115  x(1) = 0.0
116  do i = 2,5
117  x(i) = x(i-1) + h(i-1)
118  end do
119 
120  do i = 1,4
121 
122  do j = 1,4
123  asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j
124  end do
125 
126  bsys(i) = u(i) * ( h(i) )
127 
128  end do
129 
130  call solve_linear_system( asys, bsys, csys, 4 )
131 
132  dsys(1) = csys(2)
133  dsys(2) = 2.0 * csys(3)
134  dsys(3) = 3.0 * csys(4)
135 
136  tri_d(1) = 1.0
137  tri_u(1) = 0.0
138  tri_b(1) = evaluation_polynomial( dsys, 3, x(1) ) ! first edge slope
139 
140  ! Boundary conditions: right boundary
141  x(1) = 0.0
142  do i = 2,5
143  x(i) = x(i-1) + h(n-5+i)
144  end do
145 
146  do i = 1,4
147 
148  do j = 1,4
149  asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j
150  end do
151 
152  bsys(i) = u(n-4+i) * ( h(n-4+i) )
153 
154  end do
155 
156  call solve_linear_system( asys, bsys, csys, 4 )
157 
158  dsys(1) = csys(2)
159  dsys(2) = 2.0 * csys(3)
160  dsys(3) = 3.0 * csys(4)
161 
162  tri_l(n+1) = 0.0
163  tri_d(n+1) = 1.0
164  tri_b(n+1) = evaluation_polynomial( dsys, 3, x(5) ) ! last edge slope
165 
166  ! Solve tridiagonal system and assign edge values
167  call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, n+1 )
168 
169  do i = 2,n
170  edge_slopes(i,1) = tri_x(i)
171  edge_slopes(i-1,2) = tri_x(i)
172  end do
173  edge_slopes(1,1) = tri_x(1)
174  edge_slopes(n,2) = tri_x(n+1)
175 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 
)

Definition at line 183 of file regrid_edge_slopes.F90.

References polynomial_functions::evaluation_polynomial(), h_neglect, regrid_solvers::solve_linear_system(), and regrid_solvers::solve_tridiagonal_system().

Referenced by mom_remapping::build_reconstructions_1d(), and regrid_interp::regridding_set_ppolys().

183 ! -----------------------------------------------------------------------------
184 ! Fifth-order implicit estimates of edge values are based on a four-cell,
185 ! three-edge stencil. A tridiagonal system is set up and is based on
186 ! expressing the edge slopes in terms of neighboring cell averages.
187 !
188 ! The generic relationship is
189 !
190 ! \alpha u'_{i-1/2} + u'_{i+1/2} + \beta u'_{i+3/2} =
191 ! a \bar{u}_{i-1} + b \bar{u}_i + c \bar{u}_{i+1} + d \bar{u}_{i+2}
192 !
193 ! and the stencil looks like this
194 !
195 ! i-1 i i+1 i+2
196 ! ..--o------o------o------o------o--..
197 ! i-1/2 i+1/2 i+3/2
198 !
199 ! In this routine, the coefficients \alpha, \beta, a, b, c and d are
200 ! computed, the tridiagonal system is built, boundary conditions are
201 ! prescribed and the system is solved to yield edge-value estimates.
202 !
203 ! Note that the centered stencil only applies to edges 3 to N-1 (edges are
204 ! numbered 1 to n+1), which yields N-3 equations for N+1 unknowns. Two other
205 ! equations are written by using a right-biased stencil for edge 2 and a
206 ! left-biased stencil for edge N. The prescription of boundary conditions
207 ! (using sixth-order polynomials) closes the system.
208 !
209 ! CAUTION: For each edge, in order to determine the coefficients of the
210 ! implicit expression, a 6x6 linear system is solved. This may
211 ! become computationally expensive if regridding is carried out
212 ! often. Figuring out closed-form expressions for these coefficients
213 ! on nonuniform meshes turned out to be intractable.
214 ! -----------------------------------------------------------------------------
215 
216  ! Arguments
217  integer, intent(in) :: n ! Number of cells
218  real, dimension(:), intent(in) :: h ! cell averages (size N)
219  real, dimension(:), intent(in) :: u ! cell averages (size N)
220  real, dimension(:,:), intent(inout) :: edge_slopes
221 
222  ! Local variables
223  integer :: i, j, k ! loop indexes
224  real :: h0, h1, h2, h3 ! cell widths
225  real :: g, g_2, g_3 ! the following are
226  real :: g_4, g_5, g_6 ! auxiliary variables
227  real :: d2, d3, d4, d5, d6 ! to set up the systems
228  real :: n2, n3, n4, n5, n6 ! used to compute the
229  real :: h1_2, h2_2 ! the coefficients of the
230  real :: h1_3, h2_3 ! tridiagonal system
231  real :: h1_4, h2_4 ! ...
232  real :: h1_5, h2_5 ! ...
233  real :: h1_6, h2_6 ! ...
234  real :: h0ph1, h0ph1_2 ! ...
235  real :: h0ph1_3, h0ph1_4 ! ...
236  real :: h2ph3, h2ph3_2 ! ...
237  real :: h2ph3_3, h2ph3_4 ! ...
238  real :: alpha, beta ! stencil coefficients
239  real :: a, b, c, d ! "
240  real, dimension(7) :: x ! system used to enforce
241  real, dimension(6,6) :: asys ! boundary conditions
242  real, dimension(6) :: bsys, csys ! ...
243  real, dimension(5) :: dsys ! derivative
244  real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal)
245  tri_d, & ! trid. system (middle diagonal)
246  tri_u, & ! trid. system (upper diagonal)
247  tri_b, & ! trid. system (unknowns vector)
248  tri_x ! trid. system (rhs)
249 
250  ! Loop on cells (except last one)
251  do k = 2,n-2
252 
253  ! Cell widths
254  h0 = h(k-1)
255  h1 = h(k+0)
256  h2 = h(k+1)
257  h3 = h(k+2)
258 
259  ! Auxiliary calculations
260  h1_2 = h1 * h1
261  h1_3 = h1_2 * h1
262  h1_4 = h1_2 * h1_2
263  h1_5 = h1_3 * h1_2
264  h1_6 = h1_3 * h1_3
265 
266  h2_2 = h2 * h2
267  h2_3 = h2_2 * h2
268  h2_4 = h2_2 * h2_2
269  h2_5 = h2_3 * h2_2
270  h2_6 = h2_3 * h2_3
271 
272  g = h0 + h1
273  g_2 = g * g
274  g_3 = g * g_2
275  g_4 = g_2 * g_2
276  g_5 = g_4 * g
277  g_6 = g_3 * g_3
278 
279  d2 = ( h1_2 - g_2 ) / ( h0 + h_neglect )
280  d3 = ( h1_3 - g_3 ) / ( h0 + h_neglect )
281  d4 = ( h1_4 - g_4 ) / ( h0 + h_neglect )
282  d5 = ( h1_5 - g_5 ) / ( h0 + h_neglect )
283  d6 = ( h1_6 - g_6 ) / ( h0 + h_neglect )
284 
285  g = h2 + h3
286  g_2 = g * g
287  g_3 = g * g_2
288  g_4 = g_2 * g_2
289  g_5 = g_4 * g
290  g_6 = g_3 * g_3
291 
292  n2 = ( g_2 - h2_2 ) / ( h3 + h_neglect )
293  n3 = ( g_3 - h2_3 ) / ( h3 + h_neglect )
294  n4 = ( g_4 - h2_4 ) / ( h3 + h_neglect )
295  n5 = ( g_5 - h2_5 ) / ( h3 + h_neglect )
296  n6 = ( g_6 - h2_6 ) / ( h3 + h_neglect )
297 
298  ! Compute matrix entries
299  asys(1,1) = 0.0
300  asys(1,2) = 0.0
301  asys(1,3) = 1.0
302  asys(1,4) = 1.0
303  asys(1,5) = 1.0
304  asys(1,6) = 1.0
305 
306  asys(2,1) = 1.0
307  asys(2,2) = 1.0
308  asys(2,3) = -0.5 * d2
309  asys(2,4) = 0.5 * h1
310  asys(2,5) = -0.5 * h2
311  asys(2,6) = -0.5 * n2
312 
313  asys(3,1) = h1
314  asys(3,2) = - h2
315  asys(3,3) = - d3 / 6.0
316  asys(3,4) = h1_2 / 6.0
317  asys(3,5) = h2_2 / 6.0
318  asys(3,6) = n3 / 6.0
319 
320  asys(4,1) = - h1_2 / 2.0
321  asys(4,2) = - h2_2 / 2.0
322  asys(4,3) = d4 / 24.0
323  asys(4,4) = - h1_3 / 24.0
324  asys(4,5) = h2_3 / 24.0
325  asys(4,6) = n4 / 24.0
326 
327  asys(5,1) = h1_3 / 6.0
328  asys(5,2) = - h2_3 / 6.0
329  asys(5,3) = - d5 / 120.0
330  asys(5,4) = h1_4 / 120.0
331  asys(5,5) = h2_4 / 120.0
332  asys(5,6) = n5 / 120.0
333 
334  asys(6,1) = - h1_4 / 24.0
335  asys(6,2) = - h2_4 / 24.0
336  asys(6,3) = d6 / 720.0
337  asys(6,4) = - h1_5 / 720.0
338  asys(6,5) = h2_5 / 720.0
339  asys(6,6) = n6 / 720.0
340 
341  bsys(:) = (/ 0.0, -1.0, 0.0, 0.0, 0.0, 0.0 /)
342 
343  call solve_linear_system( asys, bsys, csys, 6 )
344 
345  alpha = csys(1)
346  beta = csys(2)
347  a = csys(3)
348  b = csys(4)
349  c = csys(5)
350  d = csys(6)
351 
352  tri_l(k+1) = alpha
353  tri_d(k+1) = 1.0
354  tri_u(k+1) = beta
355  tri_b(k+1) = a * u(k-1) + b * u(k) + c * u(k+1) + d * u(k+2)
356 
357  end do ! end loop on cells
358 
359  ! Use a right-biased stencil for the second row
360 
361  ! Cell widths
362  h0 = h(1)
363  h1 = h(2)
364  h2 = h(3)
365  h3 = h(4)
366 
367  ! Auxiliary calculations
368  h1_2 = h1 * h1
369  h1_3 = h1_2 * h1
370  h1_4 = h1_2 * h1_2
371  h1_5 = h1_3 * h1_2
372  h1_6 = h1_3 * h1_3
373 
374  h2_2 = h2 * h2
375  h2_3 = h2_2 * h2
376  h2_4 = h2_2 * h2_2
377  h2_5 = h2_3 * h2_2
378  h2_6 = h2_3 * h2_3
379 
380  g = h0 + h1
381  g_2 = g * g
382  g_3 = g * g_2
383  g_4 = g_2 * g_2
384  g_5 = g_4 * g
385  g_6 = g_3 * g_3
386 
387  h0ph1 = h0 + h1
388  h0ph1_2 = h0ph1 * h0ph1
389  h0ph1_3 = h0ph1_2 * h0ph1
390  h0ph1_4 = h0ph1_2 * h0ph1_2
391 
392  d2 = ( h1_2 - g_2 ) / ( h0 + h_neglect )
393  d3 = ( h1_3 - g_3 ) / ( h0 + h_neglect )
394  d4 = ( h1_4 - g_4 ) / ( h0 + h_neglect )
395  d5 = ( h1_5 - g_5 ) / ( h0 + h_neglect )
396  d6 = ( h1_6 - g_6 ) / ( h0 + h_neglect )
397 
398  g = h2 + h3
399  g_2 = g * g
400  g_3 = g * g_2
401  g_4 = g_2 * g_2
402  g_5 = g_4 * g
403  g_6 = g_3 * g_3
404 
405  n2 = ( g_2 - h2_2 ) / ( h3 + h_neglect )
406  n3 = ( g_3 - h2_3 ) / ( h3 + h_neglect )
407  n4 = ( g_4 - h2_4 ) / ( h3 + h_neglect )
408  n5 = ( g_5 - h2_5 ) / ( h3 + h_neglect )
409  n6 = ( g_6 - h2_6 ) / ( h3 + h_neglect )
410 
411  ! Compute matrix entries
412  asys(1,1) = 0.0
413  asys(1,2) = 0.0
414  asys(1,3) = 1.0
415  asys(1,4) = 1.0
416  asys(1,5) = 1.0
417  asys(1,6) = 1.0
418 
419  asys(2,1) = 1.0
420  asys(2,2) = 1.0
421  asys(2,3) = -0.5 * d2
422  asys(2,4) = 0.5 * h1
423  asys(2,5) = -0.5 * h2
424  asys(2,6) = -0.5 * n2
425 
426  asys(3,1) = h0ph1
427  asys(3,2) = 0.0
428  asys(3,3) = - d3 / 6.0
429  asys(3,4) = h1_2 / 6.0
430  asys(3,5) = h2_2 / 6.0
431  asys(3,6) = n3 / 6.0
432 
433  asys(4,1) = - h0ph1_2 / 2.0
434  asys(4,2) = 0.0
435  asys(4,3) = d4 / 24.0
436  asys(4,4) = - h1_3 / 24.0
437  asys(4,5) = h2_3 / 24.0
438  asys(4,6) = n4 / 24.0
439 
440  asys(5,1) = h0ph1_3 / 6.0
441  asys(5,2) = 0.0
442  asys(5,3) = - d5 / 120.0
443  asys(5,4) = h1_4 / 120.0
444  asys(5,5) = h2_4 / 120.0
445  asys(5,6) = n5 / 120.0
446 
447  asys(6,1) = - h0ph1_4 / 24.0
448  asys(6,2) = 0.0
449  asys(6,3) = d6 / 720.0
450  asys(6,4) = - h1_5 / 720.0
451  asys(6,5) = h2_5 / 720.0
452  asys(6,6) = n6 / 720.0
453 
454  bsys(:) = (/ 0.0, -1.0, -h1, h1_2/2.0, -h1_3/6.0, h1_4/24.0 /)
455 
456  call solve_linear_system( asys, bsys, csys, 6 )
457 
458  alpha = csys(1)
459  beta = csys(2)
460  a = csys(3)
461  b = csys(4)
462  c = csys(5)
463  d = csys(6)
464 
465  tri_l(2) = alpha
466  tri_d(2) = 1.0
467  tri_u(2) = beta
468  tri_b(2) = a * u(1) + b * u(2) + c * u(3) + d * u(4)
469 
470  ! Boundary conditions: left boundary
471  x(1) = 0.0
472  do i = 2,7
473  x(i) = x(i-1) + h(i-1)
474  end do
475 
476  do i = 1,6
477 
478  do j = 1,6
479  asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j
480  end do
481 
482  bsys(i) = u(i) * h(i)
483 
484  end do
485 
486  call solve_linear_system( asys, bsys, csys, 6 )
487 
488  dsys(1) = csys(2)
489  dsys(2) = 2.0 * csys(3)
490  dsys(3) = 3.0 * csys(4)
491  dsys(4) = 4.0 * csys(5)
492  dsys(5) = 5.0 * csys(6)
493 
494  tri_d(1) = 0.0
495  tri_d(1) = 1.0
496  tri_u(1) = 0.0
497  tri_b(1) = evaluation_polynomial( dsys, 5, x(1) ) ! first edge value
498 
499  ! Use a left-biased stencil for the second to last row
500 
501  ! Cell widths
502  h0 = h(n-3)
503  h1 = h(n-2)
504  h2 = h(n-1)
505  h3 = h(n)
506 
507  ! Auxiliary calculations
508  h1_2 = h1 * h1
509  h1_3 = h1_2 * h1
510  h1_4 = h1_2 * h1_2
511  h1_5 = h1_3 * h1_2
512  h1_6 = h1_3 * h1_3
513 
514  h2_2 = h2 * h2
515  h2_3 = h2_2 * h2
516  h2_4 = h2_2 * h2_2
517  h2_5 = h2_3 * h2_2
518  h2_6 = h2_3 * h2_3
519 
520  g = h0 + h1
521  g_2 = g * g
522  g_3 = g * g_2
523  g_4 = g_2 * g_2
524  g_5 = g_4 * g
525  g_6 = g_3 * g_3
526 
527  h2ph3 = h2 + h3
528  h2ph3_2 = h2ph3 * h2ph3
529  h2ph3_3 = h2ph3_2 * h2ph3
530  h2ph3_4 = h2ph3_2 * h2ph3_2
531 
532  d2 = ( h1_2 - g_2 ) / ( h0 + h_neglect )
533  d3 = ( h1_3 - g_3 ) / ( h0 + h_neglect )
534  d4 = ( h1_4 - g_4 ) / ( h0 + h_neglect )
535  d5 = ( h1_5 - g_5 ) / ( h0 + h_neglect )
536  d6 = ( h1_6 - g_6 ) / ( h0 + h_neglect )
537 
538  g = h2 + h3
539  g_2 = g * g
540  g_3 = g * g_2
541  g_4 = g_2 * g_2
542  g_5 = g_4 * g
543  g_6 = g_3 * g_3
544 
545  n2 = ( g_2 - h2_2 ) / ( h3 + h_neglect )
546  n3 = ( g_3 - h2_3 ) / ( h3 + h_neglect )
547  n4 = ( g_4 - h2_4 ) / ( h3 + h_neglect )
548  n5 = ( g_5 - h2_5 ) / ( h3 + h_neglect )
549  n6 = ( g_6 - h2_6 ) / ( h3 + h_neglect )
550 
551  ! Compute matrix entries
552  asys(1,1) = 0.0
553  asys(1,2) = 0.0
554  asys(1,3) = 1.0
555  asys(1,4) = 1.0
556  asys(1,5) = 1.0
557  asys(1,6) = 1.0
558 
559  asys(2,1) = 1.0
560  asys(2,2) = 1.0
561  asys(2,3) = -0.5 * d2
562  asys(2,4) = 0.5 * h1
563  asys(2,5) = -0.5 * h2
564  asys(2,6) = -0.5 * n2
565 
566  asys(3,1) = 0.0
567  asys(3,2) = - h2ph3
568  asys(3,3) = - d3 / 6.0
569  asys(3,4) = h1_2 / 6.0
570  asys(3,5) = h2_2 / 6.0
571  asys(3,6) = n3 / 6.0
572 
573  asys(4,1) = 0.0
574  asys(4,2) = - h2ph3_2 / 2.0
575  asys(4,3) = d4 / 24.0
576  asys(4,4) = - h1_3 / 24.0
577  asys(4,5) = h2_3 / 24.0
578  asys(4,6) = n4 / 24.0
579 
580  asys(5,1) = 0.0
581  asys(5,2) = - h2ph3_3 / 6.0
582  asys(5,3) = - d5 / 120.0
583  asys(5,4) = h1_4 / 120.0
584  asys(5,5) = h2_4 / 120.0
585  asys(5,6) = n5 / 120.0
586 
587  asys(6,1) = 0.0
588  asys(6,2) = - h2ph3_4 / 24.0
589  asys(6,3) = d6 / 720.0
590  asys(6,4) = - h1_5 / 720.0
591  asys(6,5) = h2_5 / 720.0
592  asys(6,6) = n6 / 720.0
593 
594  bsys(:) = (/ 0.0, -1.0, h2, h2_2/2.0, h2_3/6.0, h2_4/24.0 /)
595 
596  call solve_linear_system( asys, bsys, csys, 6 )
597 
598  alpha = csys(1)
599  beta = csys(2)
600  a = csys(3)
601  b = csys(4)
602  c = csys(5)
603  d = csys(6)
604 
605  tri_l(n) = alpha
606  tri_d(n) = 1.0
607  tri_u(n) = beta
608  tri_b(n) = a * u(n-3) + b * u(n-2) + c * u(n-1) + d * u(n)
609 
610  ! Boundary conditions: right boundary
611  x(1) = 0.0
612  do i = 2,7
613  x(i) = x(i-1) + h(n-7+i)
614  end do
615 
616  do i = 1,6
617 
618  do j = 1,6
619  asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j
620  end do
621 
622  bsys(i) = u(n-6+i) * h(n-6+i)
623 
624  end do
625 
626  call solve_linear_system( asys, bsys, csys, 6 )
627 
628  dsys(1) = csys(2)
629  dsys(2) = 2.0 * csys(3)
630  dsys(3) = 3.0 * csys(4)
631  dsys(4) = 4.0 * csys(5)
632  dsys(5) = 5.0 * csys(6)
633 
634  tri_l(n+1) = 0.0
635  tri_d(n+1) = 1.0
636  tri_u(n+1) = 0.0
637  tri_b(n+1) = evaluation_polynomial( dsys, 5, x(7) ) ! last edge value
638 
639  ! Solve tridiagonal system and assign edge values
640  call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, n+1 )
641 
642  do i = 2,n
643  edge_slopes(i,1) = tri_x(i)
644  edge_slopes(i-1,2) = tri_x(i)
645  end do
646  edge_slopes(1,1) = tri_x(1)
647  edge_slopes(n,2) = tri_x(n+1)
648 
Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ h_neglect

real, parameter regrid_edge_slopes::h_neglect = 1.E-30
private

Definition at line 25 of file regrid_edge_slopes.F90.

Referenced by edge_slopes_implicit_h3(), and edge_slopes_implicit_h5().

25 real, parameter :: h_neglect = 1.e-30