MOM6
Functions/Subroutines | Variables
regrid_edge_values Module Reference

Detailed Description

Edge value estimation for high-order resconstruction.

Functions/Subroutines

subroutine, public bound_edge_values (N, h, u, edge_val, h_neglect)
 Bound edge values by neighboring cell averages. More...
 
subroutine, public average_discontinuous_edge_values (N, edge_val)
 Replace discontinuous collocated edge values with their average. More...
 
subroutine, public check_discontinuous_edge_values (N, u, edge_val)
 Check discontinuous edge values and replace them with their average if not monotonic. More...
 
subroutine, public edge_values_explicit_h2 (N, h, u, edge_val, h_neglect)
 Compute h2 edge values (explicit second order accurate) in the same units as h. More...
 
subroutine, public edge_values_explicit_h4 (N, h, u, edge_val, h_neglect, answers_2018)
 Compute h4 edge values (explicit fourth order accurate) in the same units as h. More...
 
subroutine, public edge_values_implicit_h4 (N, h, u, edge_val, h_neglect, answers_2018)
 Compute ih4 edge values (implicit fourth order accurate) in the same units as h. More...
 
subroutine, public edge_values_implicit_h6 (N, h, u, edge_val, h_neglect, answers_2018)
 Compute ih6 edge values (implicit sixth order accurate) in the same units as h. More...
 

Variables

real, parameter hneglect_edge_dflt = 1.e-10
 The default value for cut-off minimum thickness for sum(h) in edge value inversions. More...
 
real, parameter hneglect_dflt = 1.e-30
 The default value for cut-off minimum thickness for sum(h) in other calculations. More...
 
real, parameter hminfrac = 1.e-5
 A minimum fraction for min(h)/sum(h) More...
 

Function/Subroutine Documentation

◆ average_discontinuous_edge_values()

subroutine, public regrid_edge_values::average_discontinuous_edge_values ( integer, intent(in)  N,
real, dimension(:,:), intent(inout)  edge_val 
)

Replace discontinuous collocated edge values with their average.

For each interior edge, check whether the edge values are discontinuous. If so, compute the average and replace the edge values by the average.

Parameters
[in]nNumber of cells
[in,out]edge_valEdge values that may be modified the second index size is 2.

Definition at line 141 of file regrid_edge_values.F90.

141  integer, intent(in) :: N !< Number of cells
142  real, dimension(:,:), intent(inout) :: edge_val !< Edge values that may be modified
143  !! the second index size is 2.
144  ! Local variables
145  integer :: k ! loop index
146  real :: u0_minus ! left value at given edge
147  real :: u0_plus ! right value at given edge
148  real :: u0_avg ! avg value at given edge
149 
150  ! Loop on interior edges
151  do k = 1,n-1
152 
153  ! Edge value on the left of the edge
154  u0_minus = edge_val(k,2)
155 
156  ! Edge value on the right of the edge
157  u0_plus = edge_val(k+1,1)
158 
159  if ( u0_minus /= u0_plus ) then
160  u0_avg = 0.5 * ( u0_minus + u0_plus )
161  edge_val(k,2) = u0_avg
162  edge_val(k+1,1) = u0_avg
163  endif
164 
165  enddo ! end loop on interior edges
166 

Referenced by p1m_functions::p1m_interpolation(), and p3m_functions::p3m_limiter().

Here is the caller graph for this function:

◆ bound_edge_values()

subroutine, public regrid_edge_values::bound_edge_values ( integer, intent(in)  N,
real, dimension(:), intent(in)  h,
real, dimension(:), intent(in)  u,
real, dimension(:,:), intent(inout)  edge_val,
real, intent(in), optional  h_neglect 
)

Bound edge values by neighboring cell averages.

In this routine, we loop on all cells to bound their left and right edge values by the cell averages. That is, the left edge value must lie between the left cell average and the central cell average. A similar reasoning applies to the right edge values.

Both boundary edge values are set equal to the boundary cell averages. Any extrapolation scheme is applied after this routine has been called. Therefore, boundary cells are treated as if they were local extrama.

Parameters
[in]nNumber of cells
[in]hcell widths (size N) [H]
[in]ucell average properties (size N) in arbitrary units [A]
[in,out]edge_valPotentially modified edge values [A]
[in]h_neglectA negligibly small width [H]

Definition at line 48 of file regrid_edge_values.F90.

48  integer, intent(in) :: N !< Number of cells
49  real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
50  real, dimension(:), intent(in) :: u !< cell average properties (size N) in arbitrary units [A]
51  real, dimension(:,:), intent(inout) :: edge_val !< Potentially modified edge values [A]
52  real, optional, intent(in) :: h_neglect !< A negligibly small width [H]
53  ! Local variables
54  integer :: k ! loop index
55  integer :: k0, k1, k2
56  real :: h_l, h_c, h_r ! Layer thicknesses [H]
57  real :: u_l, u_c, u_r ! Cell average properties [A]
58  real :: u0_l, u0_r ! Edge values of properties [A]
59  real :: sigma_l, sigma_c, sigma_r ! left, center and right
60  ! van Leer slopes [A H-1]
61  real :: slope ! retained PLM slope [A H-1]
62  real :: hNeglect ! A negligible thickness [H].
63 
64  hneglect = hneglect_dflt ; if (present(h_neglect)) hneglect = h_neglect
65 
66  ! Loop on cells to bound edge value
67  do k = 1,n
68 
69  ! For the sake of bounding boundary edge values, the left neighbor
70  ! of the left boundary cell is assumed to be the same as the left
71  ! boundary cell and the right neighbor of the right boundary cell
72  ! is assumed to be the same as the right boundary cell. This
73  ! effectively makes boundary cells look like extrema.
74  if ( k == 1 ) then
75  k0 = 1
76  k1 = 1
77  k2 = 2
78  elseif ( k == n ) then
79  k0 = n-1
80  k1 = n
81  k2 = n
82  else
83  k0 = k-1
84  k1 = k
85  k2 = k+1
86  endif
87 
88  ! All cells can now be treated equally
89  h_l = h(k0)
90  h_c = h(k1)
91  h_r = h(k2)
92 
93  u_l = u(k0)
94  u_c = u(k1)
95  u_r = u(k2)
96 
97  u0_l = edge_val(k,1)
98  u0_r = edge_val(k,2)
99 
100  sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + hneglect )
101  sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r + hneglect )
102  sigma_r = 2.0 * ( u_r - u_c ) / ( h_c + hneglect )
103 
104  if ( (sigma_l * sigma_r) > 0.0 ) then
105  slope = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c )
106  else
107  slope = 0.0
108  endif
109 
110  ! The limiter must be used in the local coordinate system to each cell.
111  ! Hence, we must multiply the slope by h1. The multiplication by 0.5 is
112  ! simply a way to make it useable in the limiter (cfr White and Adcroft
113  ! JCP 2008 Eqs 19 and 20)
114  slope = slope * h_c * 0.5
115 
116  if ( (u_l-u0_l)*(u0_l-u_c) < 0.0 ) then
117  u0_l = u_c - sign( min( abs(slope), abs(u0_l-u_c) ), slope )
118  endif
119 
120  if ( (u_r-u0_r)*(u0_r-u_c) < 0.0 ) then
121  u0_r = u_c + sign( min( abs(slope), abs(u0_r-u_c) ), slope )
122  endif
123 
124  ! Finally bound by neighboring cell means in case of round off
125  u0_l = max( min( u0_l, max(u_l, u_c) ), min(u_l, u_c) )
126  u0_r = max( min( u0_r, max(u_r, u_c) ), min(u_r, u_c) )
127 
128  ! Store edge values
129  edge_val(k,1) = u0_l
130  edge_val(k,2) = u0_r
131 
132  enddo ! loop on interior edges
133 

References hneglect_dflt.

Referenced by p1m_functions::p1m_interpolation(), p3m_functions::p3m_limiter(), ppm_functions::ppm_limiter_standard(), and pqm_functions::pqm_limiter().

Here is the caller graph for this function:

◆ check_discontinuous_edge_values()

subroutine, public regrid_edge_values::check_discontinuous_edge_values ( integer, intent(in)  N,
real, dimension(:), intent(in)  u,
real, dimension(:,:), intent(inout)  edge_val 
)

Check discontinuous edge values and replace them with their average if not monotonic.

For each interior edge, check whether the edge values are discontinuous. If so and if they are not monotonic, replace each edge value by their average.

Parameters
[in]nNumber of cells
[in]ucell averages (size N) in arbitrary units [A]
[in,out]edge_valCell edge values [A].

Definition at line 174 of file regrid_edge_values.F90.

174  integer, intent(in) :: N !< Number of cells
175  real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A]
176  real, dimension(:,:), intent(inout) :: edge_val !< Cell edge values [A].
177  ! Local variables
178  integer :: k ! loop index
179  real :: u0_minus ! left value at given edge [A]
180  real :: u0_plus ! right value at given edge [A]
181  real :: um_minus ! left cell average [A]
182  real :: um_plus ! right cell average [A]
183  real :: u0_avg ! avg value at given edge [A]
184 
185  ! Loop on interior cells
186  do k = 1,n-1
187 
188  ! Edge value on the left of the edge
189  u0_minus = edge_val(k,2)
190 
191  ! Edge value on the right of the edge
192  u0_plus = edge_val(k+1,1)
193 
194  ! Left cell average
195  um_minus = u(k)
196 
197  ! Right cell average
198  um_plus = u(k+1)
199 
200  if ( (u0_plus - u0_minus)*(um_plus - um_minus) < 0.0 ) then
201  u0_avg = 0.5 * ( u0_minus + u0_plus )
202  u0_avg = max( min( u0_avg, max(um_minus, um_plus) ), min(um_minus, um_plus) )
203  edge_val(k,2) = u0_avg
204  edge_val(k+1,1) = u0_avg
205  endif
206 
207  enddo ! end loop on interior edges
208 

Referenced by ppm_functions::ppm_limiter_standard(), and pqm_functions::pqm_limiter().

Here is the caller graph for this function:

◆ edge_values_explicit_h2()

subroutine, public regrid_edge_values::edge_values_explicit_h2 ( integer, intent(in)  N,
real, dimension(:), intent(in)  h,
real, dimension(:), intent(in)  u,
real, dimension(:,:), intent(inout)  edge_val,
real, intent(in), optional  h_neglect 
)

Compute h2 edge values (explicit second order accurate) in the same units as h.

Parameters
[in]nNumber of cells
[in]hcell widths (size N) [H]
[in]ucell average properties (size N) in arbitrary units [A]
[in,out]edge_valReturned edge values [A]; the second index size is 2.
[in]h_neglectA negligibly small width [H]

Definition at line 226 of file regrid_edge_values.F90.

226  integer, intent(in) :: N !< Number of cells
227  real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
228  real, dimension(:), intent(in) :: u !< cell average properties (size N) in arbitrary units [A]
229  real, dimension(:,:), intent(inout) :: edge_val !< Returned edge values [A]; the second index size is 2.
230  real, optional, intent(in) :: h_neglect !< A negligibly small width [H]
231  ! Local variables
232  integer :: k ! loop index
233  real :: h0, h1 ! cell widths [H]
234  real :: u0, u1 ! cell averages [A]
235  real :: hNeglect ! A negligible thickness [H]
236 
237  hneglect = hneglect_edge_dflt ; if (present(h_neglect)) hneglect = h_neglect
238 
239  ! Loop on interior cells
240  do k = 2,n
241 
242  h0 = h(k-1)
243  h1 = h(k)
244 
245  ! Avoid singularities when h0+h1=0
246  if (h0+h1==0.) then
247  h0 = hneglect
248  h1 = hneglect
249  endif
250 
251  u0 = u(k-1)
252  u1 = u(k)
253 
254  ! Compute left edge value
255  edge_val(k,1) = ( u0*h1 + u1*h0 ) / ( h0 + h1 )
256 
257  ! Left edge value of the current cell is equal to right edge
258  ! value of left cell
259  edge_val(k-1,2) = edge_val(k,1)
260 
261  enddo ! end loop on interior cells
262 
263  ! Boundary edge values are simply equal to the boundary cell averages
264  edge_val(1,1) = u(1)
265  edge_val(n,2) = u(n)
266 

References hneglect_edge_dflt.

◆ edge_values_explicit_h4()

subroutine, public regrid_edge_values::edge_values_explicit_h4 ( integer, intent(in)  N,
real, dimension(:), intent(in)  h,
real, dimension(:), intent(in)  u,
real, dimension(:,:), intent(inout)  edge_val,
real, intent(in), optional  h_neglect,
logical, intent(in), optional  answers_2018 
)

Compute h4 edge values (explicit fourth order accurate) in the same units as h.

Compute edge values based on fourth-order explicit estimates. These estimates are based on a cubic interpolant spanning four cells and evaluated at the location of the middle edge. An interpolant spanning cells i-2, i-1, i and i+1 is evaluated at edge i-1/2. The estimate for each edge is unique. 

  i-2    i-1     i     i+1

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

The first two edge values are estimated by evaluating the first available cubic interpolant, i.e., the interpolant spanning cells 1, 2, 3 and 4. Similarly, the last two edge values are estimated by evaluating the last available interpolant.

For this fourth-order scheme, at least four cells must exist.

Parameters
[in]nNumber of cells
[in]hcell widths (size N) [H]
[in]ucell average properties (size N) in arbitrary units [A]
[in,out]edge_valReturned edge values [A]; the second index size is 2.
[in]h_neglectA negligibly small width [H]
[in]answers_2018If true use older, less acccurate expressions.

Definition at line 289 of file regrid_edge_values.F90.

289  integer, intent(in) :: N !< Number of cells
290  real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
291  real, dimension(:), intent(in) :: u !< cell average properties (size N) in arbitrary units [A]
292  real, dimension(:,:), intent(inout) :: edge_val !< Returned edge values [A]; the second index size is 2.
293  real, optional, intent(in) :: h_neglect !< A negligibly small width [H]
294  logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions.
295 
296  ! Local variables
297  integer :: i, j
298  real :: u0, u1, u2, u3 ! temporary properties [A]
299  real :: h0, h1, h2, h3 ! temporary thicknesses [H]
300  real :: f1, f2, f3 ! auxiliary variables with various units
301  real :: e ! edge value
302  real, dimension(5) :: x ! Coordinate system with 0 at edges [H]
303  real, parameter :: C1_12 = 1.0 / 12.0
304  real :: dx, xavg ! Differences and averages of successive values of x [same units as h]
305  real, dimension(4,4) :: A ! values near the boundaries
306  real, dimension(4) :: B, C
307  real :: hNeglect ! A negligible thickness in the same units as h.
308  logical :: use_2018_answers ! If true use older, less acccurate expressions.
309 
310  use_2018_answers = .true. ; if (present(answers_2018)) use_2018_answers = answers_2018
311  hneglect = hneglect_edge_dflt ; if (present(h_neglect)) hneglect = h_neglect
312 
313  ! Loop on interior cells
314  do i = 3,n-1
315 
316  h0 = h(i-2)
317  h1 = h(i-1)
318  h2 = h(i)
319  h3 = h(i+1)
320 
321  ! Avoid singularities when consecutive pairs of h vanish
322  if (h0+h1==0. .or. h1+h2==0. .or. h2+h3==0.) then
323  f1 = max( hneglect, h0+h1+h2+h3 )
324  h0 = max( hminfrac*f1, h(i-2) )
325  h1 = max( hminfrac*f1, h(i-1) )
326  h2 = max( hminfrac*f1, h(i) )
327  h3 = max( hminfrac*f1, h(i+1) )
328  endif
329 
330  u0 = u(i-2)
331  u1 = u(i-1)
332  u2 = u(i)
333  u3 = u(i+1)
334 
335  f1 = (h0+h1) * (h2+h3) / (h1+h2)
336  f2 = u1 * h2 + u2 * h1
337  f3 = 1.0 / (h0+h1+h2) + 1.0 / (h1+h2+h3)
338 
339  e = f1 * f2 * f3
340 
341  f1 = h2 * (h2+h3) / ( (h0+h1+h2)*(h0+h1) )
342  f2 = u1*(h0+2.0*h1) - u0*h1
343 
344  e = e + f1*f2
345 
346  f1 = h1 * (h0+h1) / ( (h1+h2+h3)*(h2+h3) )
347  f2 = u2*(2.0*h2+h3) - u3*h2
348 
349  e = e + f1*f2
350 
351  e = e / ( h0 + h1 + h2 + h3)
352 
353  edge_val(i,1) = e
354  edge_val(i-1,2) = e
355 
356 #ifdef __DO_SAFETY_CHECKS__
357  if (e /= e) then
358  write(0,*) 'NaN in explicit_edge_h4 at k=',i
359  write(0,*) 'u0-u3=',u0,u1,u2,u3
360  write(0,*) 'h0-h3=',h0,h1,h2,h3
361  write(0,*) 'f1-f3=',f1,f2,f3
362  stop 'Nan during edge_values_explicit_h4'
363  endif
364 #endif
365 
366  enddo ! end loop on interior cells
367 
368  ! Determine first two edge values
369  f1 = max( hneglect, hminfrac*sum(h(1:4)) )
370  x(1) = 0.0
371  do i = 2,5
372  x(i) = x(i-1) + max(f1, h(i-1))
373  enddo
374 
375  do i = 1,4
376  dx = max(f1, h(i) )
377  if (use_2018_answers) then
378  do j = 1,4 ; a(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j) ; enddo
379  else ! Use expressions with less sensitivity to roundoff
380  xavg = 0.5 * (x(i+1) + x(i))
381  a(i,1) = dx
382  a(i,2) = dx * xavg
383  a(i,3) = dx * (xavg**2 + c1_12*dx**2)
384  a(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2)
385  endif
386 
387  b(i) = u(i) * dx
388 
389  enddo
390 
391  call solve_linear_system( a, b, c, 4 )
392 
393  ! First edge value
394  edge_val(1,1) = evaluation_polynomial( c, 4, x(1) )
395 
396  ! Second edge value
397  edge_val(1,2) = evaluation_polynomial( c, 4, x(2) )
398  edge_val(2,1) = edge_val(1,2)
399 
400 #ifdef __DO_SAFETY_CHECKS__
401  if (edge_val(1,1) /= edge_val(1,1) .or. edge_val(1,2) /= edge_val(1,2)) then
402  write(0,*) 'NaN in explicit_edge_h4 at k=',1
403  write(0,*) 'A=',a
404  write(0,*) 'B=',b
405  write(0,*) 'C=',c
406  write(0,*) 'h(1:4)=',h(1:4)
407  write(0,*) 'x=',x
408  stop 'Nan during edge_values_explicit_h4'
409  endif
410 #endif
411 
412  ! Determine last two edge values
413  f1 = max( hneglect, hminfrac*sum(h(n-3:n)) )
414  x(1) = 0.0
415  do i = 2,5
416  x(i) = x(i-1) + max(f1, h(n-5+i))
417  enddo
418 
419  do i = 1,4
420  dx = max(f1, h(n-4+i) )
421  if (use_2018_answers) then
422  do j = 1,4 ; a(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j) ; enddo
423  else ! Use expressions with less sensitivity to roundoff
424  xavg = 0.5 * (x(i+1) + x(i))
425  a(i,1) = dx
426  a(i,2) = dx * xavg
427  a(i,3) = dx * (xavg**2 + c1_12*dx**2)
428  a(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2)
429  endif
430 
431  b(i) = u(n-4+i) * dx
432 
433  enddo
434 
435  call solve_linear_system( a, b, c, 4 )
436 
437  ! Last edge value
438  edge_val(n,2) = evaluation_polynomial( c, 4, x(5) )
439 
440  ! Second to last edge value
441  edge_val(n,1) = evaluation_polynomial( c, 4, x(4) )
442  edge_val(n-1,2) = edge_val(n,1)
443 
444 #ifdef __DO_SAFETY_CHECKS__
445  if (edge_val(n,1) /= edge_val(n,1) .or. edge_val(n,2) /= edge_val(n,2)) then
446  write(0,*) 'NaN in explicit_edge_h4 at k=',n
447  write(0,*) 'A='
448  do i = 1,4
449  do j = 1,4
450  a(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j)
451  enddo
452  write(0,*) a(i,:)
453  b(i) = u(n-4+i) * ( h(n-4+i) )
454  enddo
455  write(0,*) 'B=',b
456  write(0,*) 'C=',c
457  write(0,*) 'h(:N)=',h(n-3:n)
458  write(0,*) 'x=',x
459  stop 'Nan during edge_values_explicit_h4'
460  endif
461 #endif
462 

References polynomial_functions::evaluation_polynomial(), hminfrac, hneglect_edge_dflt, and regrid_solvers::solve_linear_system().

Here is the call graph for this function:

◆ edge_values_implicit_h4()

subroutine, public regrid_edge_values::edge_values_implicit_h4 ( integer, intent(in)  N,
real, dimension(:), intent(in)  h,
real, dimension(:), intent(in)  u,
real, dimension(:,:), intent(inout)  edge_val,
real, intent(in), optional  h_neglect,
logical, intent(in), optional  answers_2018 
)

Compute ih4 edge values (implicit fourth order accurate) in the same units as h.

Compute edge values based on fourth-order implicit estimates.

Fourth-order implicit estimates of edge values are based on a two-cell stencil. A tridiagonal system is set up and is based on expressing the edge values 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-value 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_valReturned edge values [A]; the second index size is 2.
[in]h_neglectA negligibly small width [H]
[in]answers_2018If true use older, less acccurate expressions.

Definition at line 492 of file regrid_edge_values.F90.

492  integer, intent(in) :: N !< Number of cells
493  real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
494  real, dimension(:), intent(in) :: u !< cell average properties (size N) in arbitrary units [A]
495  real, dimension(:,:), intent(inout) :: edge_val !< Returned edge values [A]; the second index size is 2.
496  real, optional, intent(in) :: h_neglect !< A negligibly small width [H]
497  logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions.
498 
499  ! Local variables
500  integer :: i, j ! loop indexes
501  real :: h0, h1 ! cell widths [H]
502  real :: h0_2, h1_2, h0h1
503  real :: d2, d4
504  real :: alpha, beta ! stencil coefficients
505  real :: a, b
506  real, dimension(5) :: x ! Coordinate system with 0 at edges [H]
507  real, parameter :: C1_12 = 1.0 / 12.0
508  real :: dx, xavg ! Differences and averages of successive values of x [H]
509  real, dimension(4,4) :: Asys ! boundary conditions
510  real, dimension(4) :: Bsys, Csys
511  real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal)
512  tri_d, & ! trid. system (middle diagonal)
513  tri_u, & ! trid. system (upper diagonal)
514  tri_b, & ! trid. system (unknowns vector)
515  tri_x ! trid. system (rhs)
516  real :: hNeglect ! A negligible thickness [H]
517  logical :: use_2018_answers ! If true use older, less acccurate expressions.
518 
519  use_2018_answers = .true. ; if (present(answers_2018)) use_2018_answers = answers_2018
520  hneglect = hneglect_edge_dflt ; if (present(h_neglect)) hneglect = h_neglect
521 
522  ! Loop on cells (except last one)
523  do i = 1,n-1
524 
525  ! Get cell widths
526  h0 = h(i)
527  h1 = h(i+1)
528 
529  ! Avoid singularities when h0+h1=0
530  if (h0+h1==0.) then
531  h0 = hneglect
532  h1 = hneglect
533  endif
534 
535  ! Auxiliary calculations
536  d2 = (h0 + h1) ** 2
537  d4 = d2 ** 2
538  h0h1 = h0 * h1
539  h0_2 = h0 * h0
540  h1_2 = h1 * h1
541 
542  ! Coefficients
543  alpha = h1_2 / d2
544  beta = h0_2 / d2
545  a = 2.0 * h1_2 * ( h1_2 + 2.0 * h0_2 + 3.0 * h0h1 ) / d4
546  b = 2.0 * h0_2 * ( h0_2 + 2.0 * h1_2 + 3.0 * h0h1 ) / d4
547 
548  tri_l(i+1) = alpha
549  tri_d(i+1) = 1.0
550  tri_u(i+1) = beta
551 
552  tri_b(i+1) = a * u(i) + b * u(i+1)
553 
554  enddo ! end loop on cells
555 
556  ! Boundary conditions: left boundary
557  h0 = max( hneglect, hminfrac*sum(h(1:4)) )
558  x(1) = 0.0
559  do i = 2,5
560  x(i) = x(i-1) + max( h0, h(i-1) )
561  enddo
562 
563  do i = 1,4
564  dx = max(h0, h(i) )
565  if (use_2018_answers) then
566  do j = 1,4 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo
567  else ! Use expressions with less sensitivity to roundoff
568  xavg = 0.5 * (x(i+1) + x(i))
569  asys(i,1) = dx
570  asys(i,2) = dx * xavg
571  asys(i,3) = dx * (xavg**2 + c1_12*dx**2)
572  asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2)
573  endif
574 
575  bsys(i) = u(i) * dx
576 
577  enddo
578 
579  call solve_linear_system( asys, bsys, csys, 4 )
580 
581  tri_d(1) = 1.0
582  tri_u(1) = 0.0
583  tri_b(1) = evaluation_polynomial( csys, 4, x(1) ) ! first edge value
584 
585  ! Boundary conditions: right boundary
586  h0 = max( hneglect, hminfrac*sum(h(n-3:n)) )
587  x(1) = 0.0
588  do i = 2,5
589  x(i) = x(i-1) + max( h0, h(n-5+i) )
590  enddo
591 
592  do i = 1,4
593  dx = max(h0, h(n-4+i) )
594  if (use_2018_answers) then
595  do j = 1,4 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo
596  else ! Use expressions with less sensitivity to roundoff
597  xavg = 0.5 * (x(i+1) + x(i))
598  asys(i,1) = dx
599  asys(i,2) = dx * xavg
600  asys(i,3) = dx * (xavg**2 + c1_12*dx**2)
601  asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2)
602  endif
603  bsys(i) = u(n-4+i) * dx
604 
605  enddo
606 
607  call solve_linear_system( asys, bsys, csys, 4 )
608 
609  tri_l(n+1) = 0.0
610  tri_d(n+1) = 1.0
611  tri_b(n+1) = evaluation_polynomial( csys, 4, x(5) ) ! last edge value
612 
613  ! Solve tridiagonal system and assign edge values
614  call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, n+1 )
615 
616  do i = 2,n
617  edge_val(i,1) = tri_x(i)
618  edge_val(i-1,2) = tri_x(i)
619  enddo
620  edge_val(1,1) = tri_x(1)
621  edge_val(n,2) = tri_x(n+1)
622 

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

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

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

◆ edge_values_implicit_h6()

subroutine, public regrid_edge_values::edge_values_implicit_h6 ( integer, intent(in)  N,
real, dimension(:), intent(in)  h,
real, dimension(:), intent(in)  u,
real, dimension(:,:), intent(inout)  edge_val,
real, intent(in), optional  h_neglect,
logical, intent(in), optional  answers_2018 
)

Compute ih6 edge values (implicit sixth order accurate) in the same units as h.

Sixth-order implicit estimates of edge values are based on a four-cell, three-edge stencil. A tridiagonal system is set up and is based on expressing the edge values 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-1} + b \bar{u}_i + c \bar{u}_{i+1} + d \bar{u}_{i+2} \]

and the stencil looks like this 

    i-1     i     i+1    i+2

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

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

Note that the centered stencil only applies to edges 3 to N-1 (edges are numbered 1 to n+1), which yields N-3 equations for N+1 unknowns. Two other equations are written by using a right-biased stencil for edge 2 and a left-biased stencil for edge N. The prescription of boundary conditions (using sixth-order polynomials) closes the system.

CAUTION: For each edge, in order to determine the coefficients of the implicit expression, a 6x6 linear system is solved. This may become computationally expensive if regridding is carried out often. Figuring out closed-form expressions for these coefficients on nonuniform meshes turned out to be intractable.

Parameters
[in]nNumber of cells
[in]hcell widths (size N) [H]
[in]ucell average properties (size N) in arbitrary units [A]
[in,out]edge_valReturned edge values [A]; the second index size is 2.
[in]h_neglectA negligibly small width [H]
[in]answers_2018If true use older, less acccurate expressions.

Definition at line 661 of file regrid_edge_values.F90.

661  integer, intent(in) :: N !< Number of cells
662  real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
663  real, dimension(:), intent(in) :: u !< cell average properties (size N) in arbitrary units [A]
664  real, dimension(:,:), intent(inout) :: edge_val !< Returned edge values [A]; the second index size is 2.
665  real, optional, intent(in) :: h_neglect !< A negligibly small width [H]
666  logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions.
667 
668  ! Local variables
669  integer :: i, j, k ! loop indexes
670  real :: h0, h1, h2, h3 ! cell widths [H]
671  real :: g, g_2, g_3 ! the following are
672  real :: g_4, g_5, g_6 ! auxiliary variables
673  real :: d2, d3, d4, d5, d6 ! to set up the systems
674  real :: n2, n3, n4, n5, n6 ! used to compute the
675  real :: h1_2, h2_2 ! the coefficients of the
676  real :: h1_3, h2_3 ! tridiagonal system
677  real :: h1_4, h2_4 ! ...
678  real :: h1_5, h2_5 ! ...
679  real :: h1_6, h2_6 ! ...
680  real :: h0ph1, h0ph1_2 ! ...
681  real :: h0ph1_3, h0ph1_4 ! ...
682  real :: h2ph3, h2ph3_2 ! ...
683  real :: h2ph3_3, h2ph3_4 ! ...
684  real :: h0ph1_5, h2ph3_5 ! ...
685  real :: alpha, beta ! stencil coefficients
686  real :: a, b, c, d ! "
687  real, dimension(7) :: x ! Coordinate system with 0 at edges [same units as h]
688  real, parameter :: C1_12 = 1.0 / 12.0
689  real, parameter :: C5_6 = 5.0 / 6.0
690  real :: dx, xavg ! Differences and averages of successive values of x [same units as h]
691  real, dimension(6,6) :: Asys ! boundary conditions
692  real, dimension(6) :: Bsys, Csys ! ...
693  real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal)
694  tri_d, & ! trid. system (middle diagonal)
695  tri_u, & ! trid. system (upper diagonal)
696  tri_b, & ! trid. system (unknowns vector)
697  tri_x ! trid. system (rhs)
698  real :: hNeglect ! A negligible thickness [H].
699  logical :: use_2018_answers ! If true use older, less acccurate expressions.
700 
701  use_2018_answers = .true. ; if (present(answers_2018)) use_2018_answers = answers_2018
702  hneglect = hneglect_edge_dflt ; if (present(h_neglect)) hneglect = h_neglect
703 
704  ! Loop on cells (except last one)
705  do k = 2,n-2
706 
707  ! Cell widths
708  h0 = h(k-1)
709  h1 = h(k+0)
710  h2 = h(k+1)
711  h3 = h(k+2)
712 
713  ! Avoid singularities when h0=0 or h3=0
714  if (h0*h3==0.) then
715  g = max( hneglect, h0+h1+h2+h3 )
716  h0 = max( hminfrac*g, h0 )
717  h1 = max( hminfrac*g, h1 )
718  h2 = max( hminfrac*g, h2 )
719  h3 = max( hminfrac*g, h3 )
720  endif
721 
722  ! Auxiliary calculations
723  h1_2 = h1 * h1
724  h1_3 = h1_2 * h1
725  h1_4 = h1_2 * h1_2
726  h1_5 = h1_3 * h1_2
727  h1_6 = h1_3 * h1_3
728 
729  h2_2 = h2 * h2
730  h2_3 = h2_2 * h2
731  h2_4 = h2_2 * h2_2
732  h2_5 = h2_3 * h2_2
733  h2_6 = h2_3 * h2_3
734 
735  g = h0 + h1
736  g_2 = g * g
737  g_3 = g * g_2
738  g_4 = g_2 * g_2
739  g_5 = g_4 * g
740  g_6 = g_3 * g_3
741 
742  d2 = ( h1_2 - g_2 ) / h0
743  d3 = ( h1_3 - g_3 ) / h0
744  d4 = ( h1_4 - g_4 ) / h0
745  d5 = ( h1_5 - g_5 ) / h0
746  d6 = ( h1_6 - g_6 ) / h0
747 
748  g = h2 + h3
749  g_2 = g * g
750  g_3 = g * g_2
751  g_4 = g_2 * g_2
752  g_5 = g_4 * g
753  g_6 = g_3 * g_3
754 
755  n2 = ( g_2 - h2_2 ) / h3
756  n3 = ( g_3 - h2_3 ) / h3
757  n4 = ( g_4 - h2_4 ) / h3
758  n5 = ( g_5 - h2_5 ) / h3
759  n6 = ( g_6 - h2_6 ) / h3
760 
761  ! Compute matrix entries
762  asys(1,1) = 1.0
763  asys(1,2) = 1.0
764  asys(1,3) = -1.0
765  asys(1,4) = -1.0
766  asys(1,5) = -1.0
767  asys(1,6) = -1.0
768 
769  asys(2,1) = - h1
770  asys(2,2) = h2
771  asys(2,3) = -0.5 * d2
772  asys(2,4) = 0.5 * h1
773  asys(2,5) = -0.5 * h2
774  asys(2,6) = -0.5 * n2
775 
776  asys(3,1) = 0.5 * h1_2
777  asys(3,2) = 0.5 * h2_2
778  asys(3,3) = d3 / 6.0
779  asys(3,4) = - h1_2 / 6.0
780  asys(3,5) = - h2_2 / 6.0
781  asys(3,6) = - n3 / 6.0
782 
783  asys(4,1) = - h1_3 / 6.0
784  asys(4,2) = h2_3 / 6.0
785  asys(4,3) = - d4 / 24.0
786  asys(4,4) = h1_3 / 24.0
787  asys(4,5) = - h2_3 / 24.0
788  asys(4,6) = - n4 / 24.0
789 
790  asys(5,1) = h1_4 / 24.0
791  asys(5,2) = h2_4 / 24.0
792  asys(5,3) = d5 / 120.0
793  asys(5,4) = - h1_4 / 120.0
794  asys(5,5) = - h2_4 / 120.0
795  asys(5,6) = - n5 / 120.0
796 
797  asys(6,1) = - h1_5 / 120.0
798  asys(6,2) = h2_5 / 120.0
799  asys(6,3) = - d6 / 720.0
800  asys(6,4) = h1_5 / 720.0
801  asys(6,5) = - h2_5 / 720.0
802  asys(6,6) = - n6 / 720.0
803 
804  bsys(:) = (/ -1.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
805 
806  call solve_linear_system( asys, bsys, csys, 6 )
807 
808  alpha = csys(1)
809  beta = csys(2)
810  a = csys(3)
811  b = csys(4)
812  c = csys(5)
813  d = csys(6)
814 
815  tri_l(k+1) = alpha
816  tri_d(k+1) = 1.0
817  tri_u(k+1) = beta
818  tri_b(k+1) = a * u(k-1) + b * u(k) + c * u(k+1) + d * u(k+2)
819 
820  enddo ! end loop on cells
821 
822  ! Use a right-biased stencil for the second row
823 
824  ! Cell widths
825  h0 = h(1)
826  h1 = h(2)
827  h2 = h(3)
828  h3 = h(4)
829 
830  ! Avoid singularities when h0=0 or h3=0
831  if (h0*h3==0.) then
832  g = max( hneglect, h0+h1+h2+h3 )
833  h0 = max( hminfrac*g, h0 )
834  h1 = max( hminfrac*g, h1 )
835  h2 = max( hminfrac*g, h2 )
836  h3 = max( hminfrac*g, h3 )
837  endif
838 
839  ! Auxiliary calculations
840  h1_2 = h1 * h1
841  h1_3 = h1_2 * h1
842  h1_4 = h1_2 * h1_2
843  h1_5 = h1_3 * h1_2
844  h1_6 = h1_3 * h1_3
845 
846  h2_2 = h2 * h2
847  h2_3 = h2_2 * h2
848  h2_4 = h2_2 * h2_2
849  h2_5 = h2_3 * h2_2
850  h2_6 = h2_3 * h2_3
851 
852  g = h0 + h1
853  g_2 = g * g
854  g_3 = g * g_2
855  g_4 = g_2 * g_2
856  g_5 = g_4 * g
857  g_6 = g_3 * g_3
858 
859  h0ph1 = h0 + h1
860  h0ph1_2 = h0ph1 * h0ph1
861  h0ph1_3 = h0ph1_2 * h0ph1
862  h0ph1_4 = h0ph1_2 * h0ph1_2
863  h0ph1_5 = h0ph1_3 * h0ph1_2
864 
865  d2 = ( h1_2 - g_2 ) / h0
866  d3 = ( h1_3 - g_3 ) / h0
867  d4 = ( h1_4 - g_4 ) / h0
868  d5 = ( h1_5 - g_5 ) / h0
869  d6 = ( h1_6 - g_6 ) / h0
870 
871  g = h2 + h3
872  g_2 = g * g
873  g_3 = g * g_2
874  g_4 = g_2 * g_2
875  g_5 = g_4 * g
876  g_6 = g_3 * g_3
877 
878  n2 = ( g_2 - h2_2 ) / h3
879  n3 = ( g_3 - h2_3 ) / h3
880  n4 = ( g_4 - h2_4 ) / h3
881  n5 = ( g_5 - h2_5 ) / h3
882  n6 = ( g_6 - h2_6 ) / h3
883 
884  ! Compute matrix entries
885  asys(1,1) = 1.0
886  asys(1,2) = 1.0
887  asys(1,3) = -1.0
888  asys(1,4) = -1.0
889  asys(1,5) = -1.0
890  asys(1,6) = -1.0
891 
892  asys(2,1) = - h0ph1
893  asys(2,2) = 0.0
894  asys(2,3) = -0.5 * d2
895  asys(2,4) = 0.5 * h1
896  asys(2,5) = -0.5 * h2
897  asys(2,6) = -0.5 * n2
898 
899  asys(3,1) = 0.5 * h0ph1_2
900  asys(3,2) = 0.0
901  asys(3,3) = d3 / 6.0
902  asys(3,4) = - h1_2 / 6.0
903  asys(3,5) = - h2_2 / 6.0
904  asys(3,6) = - n3 / 6.0
905 
906  asys(4,1) = - h0ph1_3 / 6.0
907  asys(4,2) = 0.0
908  asys(4,3) = - d4 / 24.0
909  asys(4,4) = h1_3 / 24.0
910  asys(4,5) = - h2_3 / 24.0
911  asys(4,6) = - n4 / 24.0
912 
913  asys(5,1) = h0ph1_4 / 24.0
914  asys(5,2) = 0.0
915  asys(5,3) = d5 / 120.0
916  asys(5,4) = - h1_4 / 120.0
917  asys(5,5) = - h2_4 / 120.0
918  asys(5,6) = - n5 / 120.0
919 
920  asys(6,1) = - h0ph1_5 / 120.0
921  asys(6,2) = 0.0
922  asys(6,3) = - d6 / 720.0
923  asys(6,4) = h1_5 / 720.0
924  asys(6,5) = - h2_5 / 720.0
925  asys(6,6) = - n6 / 720.0
926 
927  bsys(:) = (/ -1.0, h1, -0.5*h1_2, h1_3/6.0, -h1_4/24.0, h1_5/120.0 /)
928 
929  call solve_linear_system( asys, bsys, csys, 6 )
930 
931  alpha = csys(1)
932  beta = csys(2)
933  a = csys(3)
934  b = csys(4)
935  c = csys(5)
936  d = csys(6)
937 
938  tri_l(2) = alpha
939  tri_d(2) = 1.0
940  tri_u(2) = beta
941  tri_b(2) = a * u(1) + b * u(2) + c * u(3) + d * u(4)
942 
943  ! Boundary conditions: left boundary
944  g = max( hneglect, hminfrac*sum(h(1:6)) )
945  x(1) = 0.0
946  do i = 2,7
947  x(i) = x(i-1) + max( g, h(i-1) )
948  enddo
949 
950  do i = 1,6
951  dx = max( g, h(i) )
952  if (use_2018_answers) then
953  do j = 1,6 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo
954  else ! Use expressions with less sensitivity to roundoff
955  xavg = 0.5 * (x(i+1) + x(i))
956  asys(i,1) = dx
957  asys(i,2) = dx * xavg
958  asys(i,3) = dx * (xavg**2 + c1_12*dx**2)
959  asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2)
960  asys(i,5) = dx * (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4)
961  asys(i,6) = dx * xavg * (xavg**4 + c5_6*xavg**2*dx**2 + 0.0625*dx**4)
962  endif
963  bsys(i) = u(i) * dx
964 
965  enddo
966 
967  call solve_linear_system( asys, bsys, csys, 6 )
968 
969  tri_l(1) = 0.0
970  tri_d(1) = 1.0
971  tri_u(1) = 0.0
972  tri_b(1) = evaluation_polynomial( csys, 6, x(1) ) ! first edge value
973 
974  ! Use a left-biased stencil for the second to last row
975 
976  ! Cell widths
977  h0 = h(n-3)
978  h1 = h(n-2)
979  h2 = h(n-1)
980  h3 = h(n)
981 
982  ! Avoid singularities when h0=0 or h3=0
983  if (h0*h3==0.) then
984  g = max( hneglect, h0+h1+h2+h3 )
985  h0 = max( hminfrac*g, h0 )
986  h1 = max( hminfrac*g, h1 )
987  h2 = max( hminfrac*g, h2 )
988  h3 = max( hminfrac*g, h3 )
989  endif
990 
991  ! Auxiliary calculations
992  h1_2 = h1 * h1
993  h1_3 = h1_2 * h1
994  h1_4 = h1_2 * h1_2
995  h1_5 = h1_3 * h1_2
996  h1_6 = h1_3 * h1_3
997 
998  h2_2 = h2 * h2
999  h2_3 = h2_2 * h2
1000  h2_4 = h2_2 * h2_2
1001  h2_5 = h2_3 * h2_2
1002  h2_6 = h2_3 * h2_3
1003 
1004  g = h0 + h1
1005  g_2 = g * g
1006  g_3 = g * g_2
1007  g_4 = g_2 * g_2
1008  g_5 = g_4 * g
1009  g_6 = g_3 * g_3
1010 
1011  h2ph3 = h2 + h3
1012  h2ph3_2 = h2ph3 * h2ph3
1013  h2ph3_3 = h2ph3_2 * h2ph3
1014  h2ph3_4 = h2ph3_2 * h2ph3_2
1015  h2ph3_5 = h2ph3_3 * h2ph3_2
1016 
1017  d2 = ( h1_2 - g_2 ) / h0
1018  d3 = ( h1_3 - g_3 ) / h0
1019  d4 = ( h1_4 - g_4 ) / h0
1020  d5 = ( h1_5 - g_5 ) / h0
1021  d6 = ( h1_6 - g_6 ) / h0
1022 
1023  g = h2 + h3
1024  g_2 = g * g
1025  g_3 = g * g_2
1026  g_4 = g_2 * g_2
1027  g_5 = g_4 * g
1028  g_6 = g_3 * g_3
1029 
1030  n2 = ( g_2 - h2_2 ) / h3
1031  n3 = ( g_3 - h2_3 ) / h3
1032  n4 = ( g_4 - h2_4 ) / h3
1033  n5 = ( g_5 - h2_5 ) / h3
1034  n6 = ( g_6 - h2_6 ) / h3
1035 
1036  ! Compute matrix entries
1037  asys(1,1) = 1.0
1038  asys(1,2) = 1.0
1039  asys(1,3) = -1.0
1040  asys(1,4) = -1.0
1041  asys(1,5) = -1.0
1042  asys(1,6) = -1.0
1043 
1044  asys(2,1) = 0.0
1045  asys(2,2) = h2ph3
1046  asys(2,3) = -0.5 * d2
1047  asys(2,4) = 0.5 * h1
1048  asys(2,5) = -0.5 * h2
1049  asys(2,6) = -0.5 * n2
1050 
1051  asys(3,1) = 0.0
1052  asys(3,2) = 0.5 * h2ph3_2
1053  asys(3,3) = d3 / 6.0
1054  asys(3,4) = - h1_2 / 6.0
1055  asys(3,5) = - h2_2 / 6.0
1056  asys(3,6) = - n3 / 6.0
1057 
1058  asys(4,1) = 0.0
1059  asys(4,2) = h2ph3_3 / 6.0
1060  asys(4,3) = - d4 / 24.0
1061  asys(4,4) = h1_3 / 24.0
1062  asys(4,5) = - h2_3 / 24.0
1063  asys(4,6) = - n4 / 24.0
1064 
1065  asys(5,1) = 0.0
1066  asys(5,2) = h2ph3_4 / 24.0
1067  asys(5,3) = d5 / 120.0
1068  asys(5,4) = - h1_4 / 120.0
1069  asys(5,5) = - h2_4 / 120.0
1070  asys(5,6) = - n5 / 120.0
1071 
1072  asys(6,1) = 0.0
1073  asys(6,2) = h2ph3_5 / 120.0
1074  asys(6,3) = - d6 / 720.0
1075  asys(6,4) = h1_5 / 720.0
1076  asys(6,5) = - h2_5 / 720.0
1077  asys(6,6) = - n6 / 720.0
1078 
1079  bsys(:) = (/ -1.0, -h2, -0.5*h2_2, -h2_3/6.0, -h2_4/24.0, -h2_5/120.0 /)
1080 
1081  call solve_linear_system( asys, bsys, csys, 6 )
1082 
1083  alpha = csys(1)
1084  beta = csys(2)
1085  a = csys(3)
1086  b = csys(4)
1087  c = csys(5)
1088  d = csys(6)
1089 
1090  tri_l(n) = alpha
1091  tri_d(n) = 1.0
1092  tri_u(n) = beta
1093  tri_b(n) = a * u(n-3) + b * u(n-2) + c * u(n-1) + d * u(n)
1094 
1095  ! Boundary conditions: right boundary
1096  g = max( hneglect, hminfrac*sum(h(n-5:n)) )
1097  x(1) = 0.0
1098  do i = 2,7
1099  x(i) = x(i-1) + max( g, h(n-7+i) )
1100  enddo
1101 
1102  do i = 1,6
1103  dx = max( g, h(n-6+i) )
1104  if (use_2018_answers) then
1105  do j = 1,6 ; asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo
1106  else ! Use expressions with less sensitivity to roundoff
1107  xavg = 0.5 * (x(i+1) + x(i))
1108  asys(i,1) = dx
1109  asys(i,2) = dx * xavg
1110  asys(i,3) = dx * (xavg**2 + c1_12*dx**2)
1111  asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2)
1112  asys(i,5) = dx * (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4)
1113  asys(i,6) = dx * xavg * (xavg**4 + c5_6*xavg**2*dx**2 + 0.0625*dx**4)
1114  endif
1115  bsys(i) = u(n-6+i) * dx
1116 
1117  enddo
1118 
1119  call solve_linear_system( asys, bsys, csys, 6 )
1120 
1121  tri_l(n+1) = 0.0
1122  tri_d(n+1) = 1.0
1123  tri_u(n+1) = 0.0
1124  tri_b(n+1) = evaluation_polynomial( csys, 6, x(7) ) ! last edge value
1125 
1126  ! Solve tridiagonal system and assign edge values
1127  call solve_tridiagonal_system( tri_l, tri_d, tri_u, tri_b, tri_x, n+1 )
1128 
1129  do i = 2,n
1130  edge_val(i,1) = tri_x(i)
1131  edge_val(i-1,2) = tri_x(i)
1132  enddo
1133  edge_val(1,1) = tri_x(1)
1134  edge_val(n,2) = tri_x(n+1)
1135 

References polynomial_functions::evaluation_polynomial(), hminfrac, hneglect_edge_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

◆ hminfrac

real, parameter regrid_edge_values::hminfrac = 1.e-5
private

A minimum fraction for min(h)/sum(h)

Definition at line 33 of file regrid_edge_values.F90.

33 real, parameter :: hMinFrac = 1.e-5 !< A minimum fraction for min(h)/sum(h)

Referenced by edge_values_explicit_h4(), edge_values_implicit_h4(), and edge_values_implicit_h6().

◆ hneglect_dflt

real, parameter regrid_edge_values::hneglect_dflt = 1.e-30
private

The default value for cut-off minimum thickness for sum(h) in other calculations.

Definition at line 31 of file regrid_edge_values.F90.

31 real, parameter :: hNeglect_dflt = 1.e-30 !< The default value for cut-off minimum

Referenced by bound_edge_values().

◆ hneglect_edge_dflt

real, parameter regrid_edge_values::hneglect_edge_dflt = 1.e-10
private

The default value for cut-off minimum thickness for sum(h) in edge value inversions.

Definition at line 29 of file regrid_edge_values.F90.

29 real, parameter :: hNeglect_edge_dflt = 1.e-10 !< The default value for cut-off minimum

Referenced by edge_values_explicit_h2(), edge_values_explicit_h4(), edge_values_implicit_h4(), and edge_values_implicit_h6().