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, answers_2018)
 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, answers_2018)
 Compute ih5 edge values (implicit fifth order accurate) More...
 

Variables

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

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,
logical, intent(in), optional  answers_2018 
)

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) [H]
[in]ucell average properties (size N) in arbitrary units [A]
[in,out]edge_slopesReturned edge slopes [A H-1]
[in]h_neglectA negligibly small width
[in]answers_2018If true use older, less acccurate expressions.

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

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

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

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,
real, intent(in), optional  h_neglect,
logical, intent(in), optional  answers_2018 
)

Compute ih5 edge values (implicit fifth order accurate)

Parameters
[in]nNumber of cells
[in]hcell widths (size N) [H]
[in]ucell average properties (size N) in arbitrary units [A]
[in,out]edge_slopesReturned edge slopes [A H-1]
[in]h_neglectA negligibly small width [H]
[in]answers_2018If true use older, less acccurate expressions.

Definition at line 192 of file regrid_edge_slopes.F90.

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 

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

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

Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ hneglect_dflt

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

Default negligible cell thickness.

Definition at line 16 of file regrid_edge_slopes.F90.

16 real, parameter :: hNeglect_dflt = 1.e-30 !< Default negligible cell thickness

Referenced by edge_slopes_implicit_h3(), and edge_slopes_implicit_h5().