MOM6
Data Types | Variables
regrid_interp Module Reference

Detailed Description

Vertical interpolation for regridding.

Data Types

type  interp_cs_type
 Control structure for regrid_interp module. More...
 

Variables

integer, parameter interpolation_p1m_h2 = 0
 O(h^2) More...
 
integer, parameter interpolation_p1m_h4 = 1
 O(h^2) More...
 
integer, parameter interpolation_p1m_ih4 = 2
 O(h^2) More...
 
integer, parameter interpolation_plm = 3
 O(h^2) More...
 
integer, parameter interpolation_ppm_h4 = 4
 O(h^3) More...
 
integer, parameter interpolation_ppm_ih4 = 5
 O(h^3) More...
 
integer, parameter interpolation_p3m_ih4ih3 = 6
 O(h^4) More...
 
integer, parameter interpolation_p3m_ih6ih5 = 7
 O(h^4) More...
 
integer, parameter interpolation_pqm_ih4ih3 = 8
 O(h^4) More...
 
integer, parameter interpolation_pqm_ih6ih5 = 9
 O(h^5) More...
 
integer, parameter degree_1 = 1
 Interpolant degrees. More...
 
integer, parameter degree_2 = 2
 Interpolant degrees. More...
 
integer, parameter degree_3 = 3
 Interpolant degrees. More...
 
integer, parameter degree_4 = 4
 Interpolant degrees. More...
 
integer, parameter, public degree_max = 5
 Interpolant degrees. More...
 
real, parameter, public nr_offset = 1e-6
 When the N-R algorithm produces an estimate that lies outside [0,1], the estimate is set to be equal to the boundary location, 0 or 1, plus or minus an offset, respectively, when the derivative is zero at the boundary. More...
 
integer, parameter, public nr_iterations = 8
 Maximum number of Newton-Raphson iterations. Newton-Raphson iterations are used to build the new grid by finding the coordinates associated with target densities and interpolations of degree larger than 1. More...
 
real, parameter, public nr_tolerance = 1e-12
 Tolerance for Newton-Raphson iterations (stop when increment falls below this) More...
 
subroutine, public regridding_set_ppolys (CS, densities, n0, h0, ppoly0_E, ppoly0_S, ppoly0_coefs, degree, h_neglect, h_neglect_edge)
 Builds an interpolated profile for the densities within each grid cell. More...
 
subroutine, public interpolate_grid (n0, h0, x0, ppoly0_E, ppoly0_coefs, target_values, degree, n1, h1, x1, answers_2018)
 Given target values (e.g., density), build new grid based on polynomial. More...
 
subroutine, public build_and_interpolate_grid (CS, densities, n0, h0, x0, target_values, n1, h1, x1, h_neglect, h_neglect_edge)
 Build a grid by interpolating for target values. More...
 
real function get_polynomial_coordinate (N, h, x_g, ppoly_E, ppoly_coefs, target_value, degree, answers_2018)
 Given a target value, find corresponding coordinate for given polynomial. More...
 
integer function interpolation_scheme (interp_scheme)
 Numeric value of interpolation_scheme corresponding to scheme name. More...
 
subroutine, public set_interp_scheme (CS, interp_scheme)
 Store the interpolation_scheme value in the interp_CS based on the input string. More...
 
subroutine, public set_interp_extrap (CS, extrap)
 Store the boundary_extrapolation value in the interp_CS. More...
 

Function/Subroutine Documentation

◆ build_and_interpolate_grid()

subroutine, public regrid_interp::build_and_interpolate_grid ( type(interp_cs_type), intent(in)  CS,
real, dimension(:), intent(in)  densities,
integer, intent(in)  n0,
real, dimension(:), intent(in)  h0,
real, dimension(:), intent(in)  x0,
real, dimension(:), intent(in)  target_values,
integer, intent(in)  n1,
real, dimension(:), intent(inout)  h1,
real, dimension(:), intent(inout)  x1,
real, intent(in), optional  h_neglect,
real, intent(in), optional  h_neglect_edge 
)

Build a grid by interpolating for target values.

Parameters
[in]csA control structure for regrid_interp
[in]densitiesInput cell densities [kg m-3]
[in]target_valuesTarget values of interfaces
[in]n0The number of points on the input grid
[in]h0Initial cell widths
[in]x0Source interface positions
[in]n1The number of points on the output grid
[in,out]h1Output cell widths
[in,out]x1Target interface positions
[in]h_neglectA negligibly small width for the purpose of cell reconstructions in the same units as h0.
[in]h_neglect_edgeA negligibly small width for the purpose of edge value calculations in the same units as h0.

Definition at line 307 of file regrid_interp.F90.

307  type(interp_CS_type), intent(in) :: CS !< A control structure for regrid_interp
308  real, dimension(:), intent(in) :: densities !< Input cell densities [kg m-3]
309  real, dimension(:), intent(in) :: target_values !< Target values of interfaces
310  integer, intent(in) :: n0 !< The number of points on the input grid
311  real, dimension(:), intent(in) :: h0 !< Initial cell widths
312  real, dimension(:), intent(in) :: x0 !< Source interface positions
313  integer, intent(in) :: n1 !< The number of points on the output grid
314  real, dimension(:), intent(inout) :: h1 !< Output cell widths
315  real, dimension(:), intent(inout) :: x1 !< Target interface positions
316  real, optional, intent(in) :: h_neglect !< A negligibly small width for the
317  !! purpose of cell reconstructions
318  !! in the same units as h0.
319  real, optional, intent(in) :: h_neglect_edge !< A negligibly small width
320  !! for the purpose of edge value calculations
321  !! in the same units as h0.
322 
323  real, dimension(n0,2) :: ppoly0_E, ppoly0_S
324  real, dimension(n0,DEGREE_MAX+1) :: ppoly0_C
325  integer :: degree
326 
327  call regridding_set_ppolys(cs, densities, n0, h0, ppoly0_e, ppoly0_s, ppoly0_c, &
328  degree, h_neglect, h_neglect_edge)
329  call interpolate_grid(n0, h0, x0, ppoly0_e, ppoly0_c, target_values, degree, &
330  n1, h1, x1, answers_2018=cs%answers_2018)

References interpolate_grid(), and regridding_set_ppolys().

Referenced by coord_hycom::build_hycom1_column(), coord_rho::build_rho_column(), and coord_rho::build_rho_column_iteratively().

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

◆ get_polynomial_coordinate()

real function regrid_interp::get_polynomial_coordinate ( integer, intent(in)  N,
real, dimension(:), intent(in)  h,
real, dimension(:), intent(in)  x_g,
real, dimension(:,:), intent(in)  ppoly_E,
real, dimension(:,:), intent(in)  ppoly_coefs,
real, intent(in)  target_value,
integer, intent(in)  degree,
logical, intent(in), optional  answers_2018 
)
private

Given a target value, find corresponding coordinate for given polynomial.

Here, 'ppoly' is assumed to be a piecewise discontinuous polynomial of degree 'degree' throughout the domain defined by 'grid'. A target value is given and we need to determine the corresponding grid coordinate to define the new grid.

If the target value is out of range, the grid coordinate is simply set to be equal to one of the boundary coordinates, which results in vanished layers near the boundaries.

IT IS ASSUMED THAT THE PIECEWISE POLYNOMIAL IS MONOTONICALLY INCREASING. IF THIS IS NOT THE CASE, THE NEW GRID MAY BE ILL-DEFINED.

It is assumed that the number of cells defining 'grid' and 'ppoly' are the same.

Parameters
[in]nNumber of grid cells
[in]hGrid cell thicknesses
[in]x_gGrid interface locations
[in]ppoly_eEdge values of interpolating polynomials
[in]ppoly_coefsCoefficients of interpolating polynomials
[in]target_valueTarget value to find position for
[in]degreeDegree of the interpolating polynomials
[in]answers_2018If true use older, less acccurate expressions.
Returns
The position of x_g at which target_value is found.

Definition at line 351 of file regrid_interp.F90.

351  ! Arguments
352  integer, intent(in) :: N !< Number of grid cells
353  real, dimension(:), intent(in) :: h !< Grid cell thicknesses
354  real, dimension(:), intent(in) :: x_g !< Grid interface locations
355  real, dimension(:,:), intent(in) :: ppoly_E !< Edge values of interpolating polynomials
356  real, dimension(:,:), intent(in) :: ppoly_coefs !< Coefficients of interpolating polynomials
357  real, intent(in) :: target_value !< Target value to find position for
358  integer, intent(in) :: degree !< Degree of the interpolating polynomials
359  logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions.
360  real :: x_tgt !< The position of x_g at which target_value is found.
361  ! Local variables
362  integer :: i, k ! loop indices
363  integer :: k_found ! index of target cell
364  integer :: iter
365  real :: xi0 ! normalized target coordinate
366  real, dimension(DEGREE_MAX) :: a ! polynomial coefficients
367  real :: numerator
368  real :: denominator
369  real :: delta ! Newton-Raphson increment
370  real :: x ! global target coordinate
371  real :: eps ! offset used to get away from
372  ! boundaries
373  real :: grad ! gradient during N-R iterations
374  logical :: use_2018_answers ! If true use older, less acccurate expressions.
375 
376  eps = nr_offset
377  k_found = -1
378  use_2018_answers = .true. ; if (present(answers_2018)) use_2018_answers = answers_2018
379 
380  ! If the target value is outside the range of all values, we
381  ! force the target coordinate to be equal to the lowest or
382  ! largest value, depending on which bound is overtaken
383  if ( target_value <= ppoly_e(1,1) ) then
384  x_tgt = x_g(1)
385  return ! return because there is no need to look further
386  endif
387 
388  ! Since discontinuous edge values are allowed, we check whether the target
389  ! value lies between two discontinuous edge values at interior interfaces
390  do k = 2,n
391  if ( ( target_value >= ppoly_e(k-1,2) ) .AND. &
392  ( target_value <= ppoly_e(k,1) ) ) then
393  x_tgt = x_g(k)
394  return ! return because there is no need to look further
395  exit
396  endif
397  enddo
398 
399  ! If the target value is outside the range of all values, we
400  ! force the target coordinate to be equal to the lowest or
401  ! largest value, depending on which bound is overtaken
402  if ( target_value >= ppoly_e(n,2) ) then
403  x_tgt = x_g(n+1)
404  return ! return because there is no need to look further
405  endif
406 
407  ! At this point, we know that the target value is bounded and does not
408  ! lie between discontinuous, monotonic edge values. Therefore,
409  ! there is a unique solution. We loop on all cells and find which one
410  ! contains the target value. The variable k_found holds the index value
411  ! of the cell where the taregt value lies.
412  do k = 1,n
413  if ( ( target_value > ppoly_e(k,1) ) .AND. &
414  ( target_value < ppoly_e(k,2) ) ) then
415  k_found = k
416  exit
417  endif
418  enddo
419 
420  ! At this point, 'k_found' should be strictly positive. If not, this is
421  ! a major failure because it means we could not find any target cell
422  ! despite the fact that the target value lies between the extremes. It
423  ! means there is a major problem with the interpolant. This needs to be
424  ! reported.
425  if ( k_found == -1 ) then
426  write(*,*) target_value, ppoly_e(1,1), ppoly_e(n,2)
427  write(*,*) 'Could not find target coordinate in ' //&
428  '"get_polynomial_coordinate". This is caused by an '//&
429  'inconsistent interpolant (perhaps not monotonically '//&
430  'increasing)'
431  call mom_error( fatal, 'Aborting execution' )
432  endif
433 
434  ! Reset all polynomial coefficients to 0 and copy those pertaining to
435  ! the found cell
436  a(:) = 0.0
437  do i = 1,degree+1
438  a(i) = ppoly_coefs(k_found,i)
439  enddo
440 
441  ! Guess value to start Newton-Raphson iterations (middle of cell)
442  xi0 = 0.5
443  iter = 1
444  delta = 1e10
445 
446  ! Newton-Raphson iterations
447  do
448  ! break if converged or too many iterations taken
449  if ( ( iter > nr_iterations ) .OR. &
450  ( abs(delta) < nr_tolerance ) ) then
451  exit
452  endif
453 
454  if (use_2018_answers) then
455  numerator = a(1) + a(2)*xi0 + a(3)*xi0*xi0 + a(4)*xi0*xi0*xi0 + &
456  a(5)*xi0*xi0*xi0*xi0 - target_value
457  denominator = a(2) + 2*a(3)*xi0 + 3*a(4)*xi0*xi0 + 4*a(5)*xi0*xi0*xi0
458  else ! These expressions are mathematicaly equivalent but more accurate.
459  numerator = (a(1) - target_value) + xi0*(a(2) + xi0*(a(3) + xi0*(a(4) + a(5)*xi0)))
460  denominator = a(2) + xi0*(2.*a(3) + xi0*(3.*a(4) + 4.*a(5)*xi0))
461  endif
462 
463  delta = -numerator / denominator
464 
465  xi0 = xi0 + delta
466 
467  ! Check whether new estimate is out of bounds. If the new estimate is
468  ! indeed out of bounds, we manually set it to be equal to the overtaken
469  ! bound with a small offset towards the interior when the gradient of
470  ! the function at the boundary is zero (in which case, the Newton-Raphson
471  ! algorithm does not converge).
472  if ( xi0 < 0.0 ) then
473  xi0 = 0.0
474  grad = a(2)
475  if ( grad == 0.0 ) xi0 = xi0 + eps
476  endif
477 
478  if ( xi0 > 1.0 ) then
479  xi0 = 1.0
480  if (use_2018_answers) then
481  grad = a(2) + 2*a(3) + 3*a(4) + 4*a(5)
482  else ! These expressions are mathematicaly equivalent but more accurate.
483  grad = a(2) + (2.*a(3) + (3.*a(4) + 4.*a(5)))
484  endif
485  if ( grad == 0.0 ) xi0 = xi0 - eps
486  endif
487 
488  iter = iter + 1
489  enddo ! end Newton-Raphson iterations
490 
491  x_tgt = x_g(k_found) + xi0 * h(k_found)

References mom_error_handler::mom_error(), nr_iterations, nr_offset, and nr_tolerance.

Referenced by interpolate_grid().

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

◆ interpolate_grid()

subroutine, public regrid_interp::interpolate_grid ( integer, intent(in)  n0,
real, dimension(:), intent(in)  h0,
real, dimension(:), intent(in)  x0,
real, dimension(:,:), intent(in)  ppoly0_E,
real, dimension(:,:), intent(in)  ppoly0_coefs,
real, dimension(:), intent(in)  target_values,
integer, intent(in)  degree,
integer, intent(in)  n1,
real, dimension(:), intent(inout)  h1,
real, dimension(:), intent(inout)  x1,
logical, intent(in), optional  answers_2018 
)

Given target values (e.g., density), build new grid based on polynomial.

Given the grid 'grid0' and the piecewise polynomial interpolant 'ppoly0' (possibly discontinuous), the coordinates of the new grid 'grid1' are determined by finding the corresponding target interface densities.

Parameters
[in]n0Number of points on source grid
[in]h0Thicknesses of source grid cells
[in]x0Source interface positions
[in]ppoly0_eEdge values of interpolating polynomials
[in]ppoly0_coefsCoefficients of interpolating polynomials
[in]target_valuesTarget values of interfaces
[in]degreeDegree of interpolating polynomials
[in]n1Number of points on target grid
[in,out]h1Thicknesses of target grid cells
[in,out]x1Target interface positions
[in]answers_2018If true use older, less acccurate expressions.

Definition at line 271 of file regrid_interp.F90.

271  integer, intent(in) :: n0 !< Number of points on source grid
272  real, dimension(:), intent(in) :: h0 !< Thicknesses of source grid cells
273  real, dimension(:), intent(in) :: x0 !< Source interface positions
274  real, dimension(:,:), intent(in) :: ppoly0_E !< Edge values of interpolating polynomials
275  real, dimension(:,:), intent(in) :: ppoly0_coefs !< Coefficients of interpolating polynomials
276  real, dimension(:), intent(in) :: target_values !< Target values of interfaces
277  integer, intent(in) :: degree !< Degree of interpolating polynomials
278  integer, intent(in) :: n1 !< Number of points on target grid
279  real, dimension(:), intent(inout) :: h1 !< Thicknesses of target grid cells
280  real, dimension(:), intent(inout) :: x1 !< Target interface positions
281  logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions.
282 
283  ! Local variables
284  logical :: use_2018_answers ! If true use older, less acccurate expressions.
285  integer :: k ! loop index
286  real :: t ! current interface target density
287 
288  ! Make sure boundary coordinates of new grid coincide with boundary
289  ! coordinates of previous grid
290  x1(1) = x0(1)
291  x1(n1+1) = x0(n0+1)
292 
293  ! Find coordinates for interior target values
294  do k = 2,n1
295  t = target_values(k)
296  x1(k) = get_polynomial_coordinate( n0, h0, x0, ppoly0_e, ppoly0_coefs, t, degree, &
297  answers_2018=answers_2018 )
298  h1(k-1) = x1(k) - x1(k-1)
299  enddo
300  h1(n1) = x1(n1+1) - x1(n1)
301 

References get_polynomial_coordinate().

Referenced by build_and_interpolate_grid().

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

◆ interpolation_scheme()

integer function regrid_interp::interpolation_scheme ( character(len=*), intent(in)  interp_scheme)
private

Numeric value of interpolation_scheme corresponding to scheme name.

Parameters
[in]interp_schemeName of the interpolation scheme Valid values include "P1M_H2", "P1M_H4", "P1M_IH2", "PLM", "PPM_H4", "PPM_IH4", "P3M_IH4IH3", "P3M_IH6IH5", "PQM_IH4IH3", and "PQM_IH6IH5"

Definition at line 496 of file regrid_interp.F90.

496  character(len=*), intent(in) :: interp_scheme !< Name of the interpolation scheme
497  !! Valid values include "P1M_H2", "P1M_H4", "P1M_IH2", "PLM", "PPM_H4",
498  !! "PPM_IH4", "P3M_IH4IH3", "P3M_IH6IH5", "PQM_IH4IH3", and "PQM_IH6IH5"
499 
500  select case ( uppercase(trim(interp_scheme)) )
501  case ("P1M_H2"); interpolation_scheme = interpolation_p1m_h2
502  case ("P1M_H4"); interpolation_scheme = interpolation_p1m_h4
503  case ("P1M_IH2"); interpolation_scheme = interpolation_p1m_ih4
504  case ("PLM"); interpolation_scheme = interpolation_plm
505  case ("PPM_H4"); interpolation_scheme = interpolation_ppm_h4
506  case ("PPM_IH4"); interpolation_scheme = interpolation_ppm_ih4
507  case ("P3M_IH4IH3"); interpolation_scheme = interpolation_p3m_ih4ih3
508  case ("P3M_IH6IH5"); interpolation_scheme = interpolation_p3m_ih6ih5
509  case ("PQM_IH4IH3"); interpolation_scheme = interpolation_pqm_ih4ih3
510  case ("PQM_IH6IH5"); interpolation_scheme = interpolation_pqm_ih6ih5
511  case default ; call mom_error(fatal, "regrid_interp: "//&
512  "Unrecognized choice for INTERPOLATION_SCHEME ("//trim(interp_scheme)//").")
513  end select

References interpolation_p1m_h2, interpolation_p1m_h4, interpolation_p1m_ih4, interpolation_p3m_ih4ih3, interpolation_p3m_ih6ih5, interpolation_plm, interpolation_ppm_h4, interpolation_ppm_ih4, interpolation_pqm_ih4ih3, interpolation_pqm_ih6ih5, mom_error_handler::mom_error(), and mom_string_functions::uppercase().

Referenced by set_interp_scheme().

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

◆ regridding_set_ppolys()

subroutine, public regrid_interp::regridding_set_ppolys ( type(interp_cs_type), intent(in)  CS,
real, dimension(:), intent(in)  densities,
integer, intent(in)  n0,
real, dimension(:), intent(in)  h0,
real, dimension(:,:), intent(inout)  ppoly0_E,
real, dimension(:,:), intent(inout)  ppoly0_S,
real, dimension(:,:), intent(inout)  ppoly0_coefs,
integer, intent(inout)  degree,
real, intent(in), optional  h_neglect,
real, intent(in), optional  h_neglect_edge 
)

Builds an interpolated profile for the densities within each grid cell.

It may happen that, given a high-order interpolator, the number of available layers is insufficient (e.g., there are two available layers for a third-order PPM ih4 scheme). In these cases, we resort to the simplest continuous linear scheme (P1M h2).

Parameters
[in]csInterpolation control structure
[in]densitiesActual cell densities
[in]n0Number of cells on source grid
[in]h0cell widths on source grid
[in,out]ppoly0_eEdge value of polynomial
[in,out]ppoly0_sEdge slope of polynomial
[in,out]ppoly0_coefsCoefficients of polynomial
[in,out]degreeThe degree of the polynomials
[in]h_neglectA negligibly small width for the purpose of cell reconstructions in the same units as h0.
[in]h_neglect_edgeA negligibly small width for the purpose of edge value calculations in the same units as h0.

Definition at line 80 of file regrid_interp.F90.

80  type(interp_CS_type),intent(in) :: CS !< Interpolation control structure
81  real, dimension(:), intent(in) :: densities !< Actual cell densities
82  integer, intent(in) :: n0 !< Number of cells on source grid
83  real, dimension(:), intent(in) :: h0 !< cell widths on source grid
84  real, dimension(:,:),intent(inout) :: ppoly0_E !< Edge value of polynomial
85  real, dimension(:,:),intent(inout) :: ppoly0_S !< Edge slope of polynomial
86  real, dimension(:,:),intent(inout) :: ppoly0_coefs !< Coefficients of polynomial
87  integer, intent(inout) :: degree !< The degree of the polynomials
88  real, optional, intent(in) :: h_neglect !< A negligibly small width for the
89  !! purpose of cell reconstructions
90  !! in the same units as h0.
91  real, optional, intent(in) :: h_neglect_edge !< A negligibly small width
92  !! for the purpose of edge value calculations
93  !! in the same units as h0.
94  ! Local variables
95  logical :: extrapolate
96 
97  ! Reset piecewise polynomials
98  ppoly0_e(:,:) = 0.0
99  ppoly0_s(:,:) = 0.0
100  ppoly0_coefs(:,:) = 0.0
101 
102  extrapolate = cs%boundary_extrapolation
103 
104  ! Compute the interpolated profile of the density field and build grid
105  select case (cs%interpolation_scheme)
106 
107  case ( interpolation_p1m_h2 )
108  degree = degree_1
109  call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
110  call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
111  if (extrapolate) then
112  call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
113  endif
114 
115  case ( interpolation_p1m_h4 )
116  degree = degree_1
117  if ( n0 >= 4 ) then
118  call edge_values_explicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge, answers_2018=cs%answers_2018 )
119  else
120  call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
121  endif
122  call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
123  if (extrapolate) then
124  call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
125  endif
126 
127  case ( interpolation_p1m_ih4 )
128  degree = degree_1
129  if ( n0 >= 4 ) then
130  call edge_values_implicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge, answers_2018=cs%answers_2018 )
131  else
132  call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
133  endif
134  call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
135  if (extrapolate) then
136  call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
137  endif
138 
139  case ( interpolation_plm )
140  degree = degree_1
141  call plm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
142  if (extrapolate) then
143  call plm_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
144  endif
145 
146  case ( interpolation_ppm_h4 )
147  if ( n0 >= 4 ) then
148  degree = degree_2
149  call edge_values_explicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge, answers_2018=cs%answers_2018 )
150  call ppm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
151  if (extrapolate) then
152  call ppm_boundary_extrapolation( n0, h0, densities, ppoly0_e, &
153  ppoly0_coefs, h_neglect )
154  endif
155  else
156  degree = degree_1
157  call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
158  call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
159  if (extrapolate) then
160  call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
161  endif
162  endif
163 
164  case ( interpolation_ppm_ih4 )
165  if ( n0 >= 4 ) then
166  degree = degree_2
167  call edge_values_implicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge, answers_2018=cs%answers_2018 )
168  call ppm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
169  if (extrapolate) then
170  call ppm_boundary_extrapolation( n0, h0, densities, ppoly0_e, &
171  ppoly0_coefs, h_neglect )
172  endif
173  else
174  degree = degree_1
175  call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
176  call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
177  if (extrapolate) then
178  call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
179  endif
180  endif
181 
182  case ( interpolation_p3m_ih4ih3 )
183  if ( n0 >= 4 ) then
184  degree = degree_3
185  call edge_values_implicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge, answers_2018=cs%answers_2018 )
186  call edge_slopes_implicit_h3( n0, h0, densities, ppoly0_s, h_neglect, answers_2018=cs%answers_2018 )
187  call p3m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_s, &
188  ppoly0_coefs, h_neglect )
189  if (extrapolate) then
190  call p3m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_s, &
191  ppoly0_coefs, h_neglect, h_neglect_edge )
192  endif
193  else
194  degree = degree_1
195  call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
196  call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
197  if (extrapolate) then
198  call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
199  endif
200  endif
201 
202  case ( interpolation_p3m_ih6ih5 )
203  if ( n0 >= 6 ) then
204  degree = degree_3
205  call edge_values_implicit_h6( n0, h0, densities, ppoly0_e, h_neglect_edge, answers_2018=cs%answers_2018 )
206  call edge_slopes_implicit_h5( n0, h0, densities, ppoly0_s, h_neglect, answers_2018=cs%answers_2018 )
207  call p3m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_s, &
208  ppoly0_coefs, h_neglect )
209  if (extrapolate) then
210  call p3m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_s, &
211  ppoly0_coefs, h_neglect, h_neglect_edge )
212  endif
213  else
214  degree = degree_1
215  call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
216  call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
217  if (extrapolate) then
218  call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
219  endif
220  endif
221 
222  case ( interpolation_pqm_ih4ih3 )
223  if ( n0 >= 4 ) then
224  degree = degree_4
225  call edge_values_implicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge, answers_2018=cs%answers_2018 )
226  call edge_slopes_implicit_h3( n0, h0, densities, ppoly0_s, h_neglect, answers_2018=cs%answers_2018 )
227  call pqm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_s, &
228  ppoly0_coefs, h_neglect )
229  if (extrapolate) then
230  call pqm_boundary_extrapolation_v1( n0, h0, densities, ppoly0_e, ppoly0_s, &
231  ppoly0_coefs, h_neglect )
232  endif
233  else
234  degree = degree_1
235  call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
236  call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
237  if (extrapolate) then
238  call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
239  endif
240  endif
241 
242  case ( interpolation_pqm_ih6ih5 )
243  if ( n0 >= 6 ) then
244  degree = degree_4
245  call edge_values_implicit_h6( n0, h0, densities, ppoly0_e, h_neglect_edge, answers_2018=cs%answers_2018 )
246  call edge_slopes_implicit_h5( n0, h0, densities, ppoly0_s, h_neglect, answers_2018=cs%answers_2018 )
247  call pqm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_s, &
248  ppoly0_coefs, h_neglect )
249  if (extrapolate) then
250  call pqm_boundary_extrapolation_v1( n0, h0, densities, ppoly0_e, ppoly0_s, &
251  ppoly0_coefs, h_neglect )
252  endif
253  else
254  degree = degree_1
255  call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
256  call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
257  if (extrapolate) then
258  call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
259  endif
260  endif
261  end select

References degree_1, degree_2, degree_3, degree_4, regrid_edge_slopes::edge_slopes_implicit_h3(), regrid_edge_slopes::edge_slopes_implicit_h5(), regrid_edge_values::edge_values_implicit_h4(), regrid_edge_values::edge_values_implicit_h6(), interpolation_p1m_h2, interpolation_p1m_h4, interpolation_p1m_ih4, interpolation_p3m_ih4ih3, interpolation_p3m_ih6ih5, interpolation_plm, interpolation_ppm_h4, interpolation_ppm_ih4, interpolation_pqm_ih4ih3, interpolation_pqm_ih6ih5, p1m_functions::p1m_boundary_extrapolation(), p1m_functions::p1m_interpolation(), p3m_functions::p3m_boundary_extrapolation(), p3m_functions::p3m_interpolation(), plm_functions::plm_boundary_extrapolation(), plm_functions::plm_reconstruction(), ppm_functions::ppm_boundary_extrapolation(), ppm_functions::ppm_reconstruction(), pqm_functions::pqm_boundary_extrapolation_v1(), and pqm_functions::pqm_reconstruction().

Referenced by build_and_interpolate_grid().

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

◆ set_interp_extrap()

subroutine, public regrid_interp::set_interp_extrap ( type(interp_cs_type), intent(inout)  CS,
logical, intent(in)  extrap 
)

Store the boundary_extrapolation value in the interp_CS.

Parameters
[in,out]csA control structure for regrid_interp
[in]extrapIndicate whether high-order boundary extrapolation should be used in boundary cells

Definition at line 528 of file regrid_interp.F90.

528  type(interp_CS_type), intent(inout) :: CS !< A control structure for regrid_interp
529  logical, intent(in) :: extrap !< Indicate whether high-order boundary
530  !! extrapolation should be used in boundary cells
531 
532  cs%boundary_extrapolation = extrap

Referenced by mom_regridding::set_regrid_params().

Here is the caller graph for this function:

◆ set_interp_scheme()

subroutine, public regrid_interp::set_interp_scheme ( type(interp_cs_type), intent(inout)  CS,
character(len=*), intent(in)  interp_scheme 
)

Store the interpolation_scheme value in the interp_CS based on the input string.

Parameters
[in,out]csA control structure for regrid_interp
[in]interp_schemeName of the interpolation scheme Valid values include "P1M_H2", "P1M_H4", "P1M_IH2", "PLM", "PPM_H4", "PPM_IH4", "P3M_IH4IH3", "P3M_IH6IH5", "PQM_IH4IH3", and "PQM_IH6IH5"

Definition at line 518 of file regrid_interp.F90.

518  type(interp_CS_type), intent(inout) :: CS !< A control structure for regrid_interp
519  character(len=*), intent(in) :: interp_scheme !< Name of the interpolation scheme
520  !! Valid values include "P1M_H2", "P1M_H4", "P1M_IH2", "PLM", "PPM_H4",
521  !! "PPM_IH4", "P3M_IH4IH3", "P3M_IH6IH5", "PQM_IH4IH3", and "PQM_IH6IH5"
522 
523  cs%interpolation_scheme = interpolation_scheme(interp_scheme)

References interpolation_scheme().

Referenced by mom_regridding::set_regrid_params().

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

Variable Documentation

◆ degree_1

integer, parameter regrid_interp::degree_1 = 1
private

Interpolant degrees.

Definition at line 55 of file regrid_interp.F90.

55 integer, parameter :: DEGREE_1 = 1, degree_2 = 2, degree_3 = 3, degree_4 = 4

Referenced by regridding_set_ppolys().

◆ degree_2

integer, parameter regrid_interp::degree_2 = 2
private

Interpolant degrees.

Definition at line 55 of file regrid_interp.F90.

Referenced by regridding_set_ppolys().

◆ degree_3

integer, parameter regrid_interp::degree_3 = 3
private

Interpolant degrees.

Definition at line 55 of file regrid_interp.F90.

Referenced by regridding_set_ppolys().

◆ degree_4

integer, parameter regrid_interp::degree_4 = 4
private

Interpolant degrees.

Definition at line 55 of file regrid_interp.F90.

Referenced by regridding_set_ppolys().

◆ degree_max

integer, parameter, public regrid_interp::degree_max = 5

Interpolant degrees.

Definition at line 56 of file regrid_interp.F90.

56 integer, public, parameter :: DEGREE_MAX = 5

◆ interpolation_p1m_h2

integer, parameter regrid_interp::interpolation_p1m_h2 = 0

O(h^2)

Definition at line 43 of file regrid_interp.F90.

43 integer, parameter :: INTERPOLATION_P1M_H2 = 0 !< O(h^2)

Referenced by interpolation_scheme(), and regridding_set_ppolys().

◆ interpolation_p1m_h4

integer, parameter regrid_interp::interpolation_p1m_h4 = 1
private

O(h^2)

Definition at line 44 of file regrid_interp.F90.

44 integer, parameter :: INTERPOLATION_P1M_H4 = 1 !< O(h^2)

Referenced by interpolation_scheme(), and regridding_set_ppolys().

◆ interpolation_p1m_ih4

integer, parameter regrid_interp::interpolation_p1m_ih4 = 2
private

O(h^2)

Definition at line 45 of file regrid_interp.F90.

45 integer, parameter :: INTERPOLATION_P1M_IH4 = 2 !< O(h^2)

Referenced by interpolation_scheme(), and regridding_set_ppolys().

◆ interpolation_p3m_ih4ih3

integer, parameter regrid_interp::interpolation_p3m_ih4ih3 = 6
private

O(h^4)

Definition at line 49 of file regrid_interp.F90.

49 integer, parameter :: INTERPOLATION_P3M_IH4IH3 = 6 !< O(h^4)

Referenced by interpolation_scheme(), and regridding_set_ppolys().

◆ interpolation_p3m_ih6ih5

integer, parameter regrid_interp::interpolation_p3m_ih6ih5 = 7
private

O(h^4)

Definition at line 50 of file regrid_interp.F90.

50 integer, parameter :: INTERPOLATION_P3M_IH6IH5 = 7 !< O(h^4)

Referenced by interpolation_scheme(), and regridding_set_ppolys().

◆ interpolation_plm

integer, parameter regrid_interp::interpolation_plm = 3
private

O(h^2)

Definition at line 46 of file regrid_interp.F90.

46 integer, parameter :: INTERPOLATION_PLM = 3 !< O(h^2)

Referenced by interpolation_scheme(), and regridding_set_ppolys().

◆ interpolation_ppm_h4

integer, parameter regrid_interp::interpolation_ppm_h4 = 4
private

O(h^3)

Definition at line 47 of file regrid_interp.F90.

47 integer, parameter :: INTERPOLATION_PPM_H4 = 4 !< O(h^3)

Referenced by interpolation_scheme(), and regridding_set_ppolys().

◆ interpolation_ppm_ih4

integer, parameter regrid_interp::interpolation_ppm_ih4 = 5
private

O(h^3)

Definition at line 48 of file regrid_interp.F90.

48 integer, parameter :: INTERPOLATION_PPM_IH4 = 5 !< O(h^3)

Referenced by interpolation_scheme(), and regridding_set_ppolys().

◆ interpolation_pqm_ih4ih3

integer, parameter regrid_interp::interpolation_pqm_ih4ih3 = 8
private

O(h^4)

Definition at line 51 of file regrid_interp.F90.

51 integer, parameter :: INTERPOLATION_PQM_IH4IH3 = 8 !< O(h^4)

Referenced by interpolation_scheme(), and regridding_set_ppolys().

◆ interpolation_pqm_ih6ih5

integer, parameter regrid_interp::interpolation_pqm_ih6ih5 = 9
private

O(h^5)

Definition at line 52 of file regrid_interp.F90.

52 integer, parameter :: INTERPOLATION_PQM_IH6IH5 = 9 !< O(h^5)

Referenced by interpolation_scheme(), and regridding_set_ppolys().

◆ nr_iterations

integer, parameter, public regrid_interp::nr_iterations = 8

Maximum number of Newton-Raphson iterations. Newton-Raphson iterations are used to build the new grid by finding the coordinates associated with target densities and interpolations of degree larger than 1.

Definition at line 66 of file regrid_interp.F90.

66 integer, public, parameter :: NR_ITERATIONS = 8

Referenced by get_polynomial_coordinate(), and coord_slight::rho_interfaces_col().

◆ nr_offset

real, parameter, public regrid_interp::nr_offset = 1e-6

When the N-R algorithm produces an estimate that lies outside [0,1], the estimate is set to be equal to the boundary location, 0 or 1, plus or minus an offset, respectively, when the derivative is zero at the boundary.

Definition at line 62 of file regrid_interp.F90.

62 real, public, parameter :: NR_OFFSET = 1e-6

Referenced by get_polynomial_coordinate().

◆ nr_tolerance

real, parameter, public regrid_interp::nr_tolerance = 1e-12

Tolerance for Newton-Raphson iterations (stop when increment falls below this)

Definition at line 68 of file regrid_interp.F90.

68 real, public, parameter :: NR_TOLERANCE = 1e-12

Referenced by get_polynomial_coordinate(), and coord_slight::rho_interfaces_col().