regrid_interp Module Reference

Detailed Description

Vertical interpolation for regridding.

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

Functions/Subroutines
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, edge_values, 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...

Variables
integer, parameter	interpolation_p1m_h2 = 0
	O(h^2)
integer, parameter	interpolation_p1m_h4 = 1
	O(h^2)
integer, parameter	interpolation_p1m_ih4 = 2
	O(h^2)
integer, parameter	interpolation_plm = 3
	O(h^2)
integer, parameter	interpolation_ppm_h4 = 4
	O(h^3)
integer, parameter	interpolation_ppm_ih4 = 5
	O(h^3)
integer, parameter	interpolation_p3m_ih4ih3 = 6
	O(h^4)
integer, parameter	interpolation_p3m_ih6ih5 = 7
	O(h^4)
integer, parameter	interpolation_pqm_ih4ih3 = 8
	O(h^4)
integer, parameter	interpolation_pqm_ih6ih5 = 9
	O(h^5)
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 [nondim].
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.
real, parameter, public	nr_tolerance = 1e-12
	Tolerance for Newton-Raphson iterations (stop when increment falls below this) [nondim].
integer, parameter	degree_1 = 1
	Interpolant degrees.
integer, parameter	degree_2 = 2
	Interpolant degrees.
integer, parameter	degree_3 = 3
	Interpolant degrees.
integer, parameter	degree_4 = 4
	Interpolant degrees.
integer, parameter, public	degree_max = 5
	Interpolant degrees.

Function/Subroutine Documentation

build_and_interpolate_grid()

subroutine, public regrid_interp::build_and_interpolate_grid	(	type(interp_cs_type), intent(in)	CS,
		real, dimension(n0), intent(in)	densities,
		integer, intent(in)	n0,
		real, dimension(n0), intent(in)	h0,
		real, dimension(n0+1), intent(in)	x0,
		real, dimension(n1+1), intent(in)	target_values,
		integer, intent(in)	n1,
		real, dimension(n1), intent(inout)	h1,
		real, dimension(n1+1), 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]	cs	A control structure for regrid_interp
[in]	n0	The number of points on the input grid
[in]	n1	The number of points on the output grid
[in]	densities	Input cell densities [R ~> kg m-3]
[in]	target_values	Target values of interfaces [R ~> kg m-3]
[in]	h0	Initial cell widths [H]
[in]	x0	Source interface positions [H]
[in,out]	h1	Output cell widths [H]
[in,out]	x1	Target interface positions [H]
[in]	h_neglect	A negligibly small width for the purpose of cell reconstructions [H] in the same units as h0.
[in]	h_neglect_edge	A negligibly small width for the purpose of edge value calculations [H] in the same units as h0.

Definition at line 309 of file regrid_interp.F90.

  type(interp_CS_type),  intent(in)    :: CS  !< A control structure for regrid_interp
  integer,               intent(in)    :: n0  !< The number of points on the input grid
  integer,               intent(in)    :: n1  !< The number of points on the output grid
  real, dimension(n0),   intent(in)    :: densities !< Input cell densities [R ~> kg m-3]
  real, dimension(n1+1), intent(in)    :: target_values !< Target values of interfaces [R ~> kg m-3]
  real, dimension(n0),   intent(in)    :: h0  !< Initial cell widths [H]
  real, dimension(n0+1), intent(in)    :: x0  !< Source interface positions [H]
  real, dimension(n1),   intent(inout) :: h1  !< Output cell widths [H]
  real, dimension(n1+1), intent(inout) :: x1  !< Target interface positions [H]
  real,        optional, intent(in)    :: h_neglect !< A negligibly small width for the
                                           !! purpose of cell reconstructions [H]
                                           !! in the same units as h0.
  real,        optional, intent(in)    :: h_neglect_edge !< A negligibly small width
                                           !! for the purpose of edge value calculations [H]
                                           !! in the same units as h0.
 
  real, dimension(n0,2) :: ppoly0_E   ! Polynomial edge values [R ~> kg m-3]
  real, dimension(n0,2) :: ppoly0_S   ! Polynomial edge slopes [R H-1]
  real, dimension(n0,DEGREE_MAX+1) :: ppoly0_C  ! Polynomial interpolant coeficients on the local 0-1 grid [R ~> kg m-3]
  integer :: degree
 
  call regridding_set_ppolys(cs, densities, n0, h0, ppoly0_e, ppoly0_s, ppoly0_c, &
       degree, h_neglect, h_neglect_edge)
  call interpolate_grid(n0, h0, x0, ppoly0_e, ppoly0_c, target_values, degree, &
       n1, h1, x1, answers_2018=cs%answers_2018)

get_polynomial_coordinate()

real function regrid_interp::get_polynomial_coordinate ( integer, intent(in)  N, real, dimension(n), intent(in)  h, real, dimension(n+1), intent(in)  x_g, real, dimension(n,2), intent(in)  edge_values, real, dimension(n,degree_max+1), 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]	n	Number of grid cells
[in]	h	Grid cell thicknesses [H]
[in]	x_g	Grid interface locations [H]
[in]	edge_values	Edge values of interpolating polynomials [A]
[in]	ppoly_coefs	Coefficients of interpolating polynomials [A]
[in]	target_value	Target value to find position for [A]
[in]	degree	Degree of the interpolating polynomials
[in]	answers_2018	If true use older, less acccurate expressions.

Returns

The position of x_g at which target_value is found [H]

Definition at line 354 of file regrid_interp.F90.

  ! Arguments
  integer,              intent(in) :: N            !< Number of grid cells
  real, dimension(N),   intent(in) :: h            !< Grid cell thicknesses [H]
  real, dimension(N+1), intent(in) :: x_g          !< Grid interface locations [H]
  real, dimension(N,2), intent(in) :: edge_values  !< Edge values of interpolating polynomials [A]
  real, dimension(N,DEGREE_MAX+1), intent(in) :: ppoly_coefs  !< Coefficients of interpolating polynomials [A]
  real,                 intent(in) :: target_value !< Target value to find position for [A]
  integer,              intent(in) :: degree       !< Degree of the interpolating polynomials
  logical,    optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions.
  real                             :: x_tgt        !< The position of x_g at which target_value is found [H]
 
  ! Local variables
  real                        :: xi0         ! normalized target coordinate [nondim]
  real, dimension(DEGREE_MAX) :: a           ! polynomial coefficients [A]
  real                        :: numerator
  real                        :: denominator
  real                        :: delta       ! Newton-Raphson increment [nondim]
!   real                        :: x           ! global target coordinate
  real                        :: eps         ! offset used to get away from boundaries [nondim]
  real                        :: grad        ! gradient during N-R iterations [A]
  integer :: i, k, iter  ! loop indices
  integer :: k_found     ! index of target cell
  character(len=200) :: mesg
  logical :: use_2018_answers  ! If true use older, less acccurate expressions.
 
  eps = nr_offset
  k_found = -1
  use_2018_answers = .true. ; if (present(answers_2018)) use_2018_answers = answers_2018
 
  ! If the target value is outside the range of all values, we
  ! force the target coordinate to be equal to the lowest or
  ! largest value, depending on which bound is overtaken
  if ( target_value <= edge_values(1,1) ) then
    x_tgt = x_g(1)
    return  ! return because there is no need to look further
  endif
 
  ! Since discontinuous edge values are allowed, we check whether the target
  ! value lies between two discontinuous edge values at interior interfaces
  do k = 2,n
    if ( ( target_value >= edge_values(k-1,2) ) .AND. ( target_value <= edge_values(k,1) ) ) then
      x_tgt = x_g(k)
      return   ! return because there is no need to look further
    endif
  enddo
 
  ! If the target value is outside the range of all values, we
  ! force the target coordinate to be equal to the lowest or
  ! largest value, depending on which bound is overtaken
  if ( target_value >= edge_values(n,2) ) then
    x_tgt = x_g(n+1)
    return  ! return because there is no need to look further
  endif
 
  ! At this point, we know that the target value is bounded and does not
  ! lie between discontinuous, monotonic edge values. Therefore,
  ! there is a unique solution. We loop on all cells and find which one
  ! contains the target value. The variable k_found holds the index value
  ! of the cell where the taregt value lies.
  do k = 1,n
    if ( ( target_value > edge_values(k,1) ) .AND. ( target_value < edge_values(k,2) ) ) then
      k_found = k
      exit
    endif
  enddo
 
  ! At this point, 'k_found' should be strictly positive. If not, this is
  ! a major failure because it means we could not find any target cell
  ! despite the fact that the target value lies between the extremes. It
  ! means there is a major problem with the interpolant. This needs to be
  ! reported.
  if ( k_found == -1 ) then
    write(mesg,*) 'Could not find target coordinate', target_value, 'in get_polynomial_coordinate. This is '//&
                  'caused by an inconsistent interpolant (perhaps not monotonically increasing):', &
                  target_value, edge_values(1,1), edge_values(n,2)
    call mom_error( fatal, mesg )
  endif
 
  ! Reset all polynomial coefficients to 0 and copy those pertaining to
  ! the found cell
  a(:) = 0.0
  do i = 1,degree+1
    a(i) = ppoly_coefs(k_found,i)
  enddo
 
  ! Guess the middle of the cell to start Newton-Raphson iterations
  xi0 = 0.5
 
  ! Newton-Raphson iterations
  do iter = 1,nr_iterations
 
    if (use_2018_answers) then
      numerator = a(1) + a(2)*xi0 + a(3)*xi0*xi0 + a(4)*xi0*xi0*xi0 + &
                  a(5)*xi0*xi0*xi0*xi0 - target_value
      denominator = a(2) + 2*a(3)*xi0 + 3*a(4)*xi0*xi0 + 4*a(5)*xi0*xi0*xi0
    else  ! These expressions are mathematicaly equivalent but more accurate.
      numerator = (a(1) - target_value) + xi0*(a(2) + xi0*(a(3) + xi0*(a(4) + a(5)*xi0)))
      denominator = a(2) + xi0*(2.*a(3) + xi0*(3.*a(4) + 4.*a(5)*xi0))
    endif
 
    delta = -numerator / denominator
 
    xi0 = xi0 + delta
 
    ! Check whether new estimate is out of bounds. If the new estimate is
    ! indeed out of bounds, we manually set it to be equal to the overtaken
    ! bound with a small offset towards the interior when the gradient of
    ! the function at the boundary is zero (in which case, the Newton-Raphson
    ! algorithm does not converge).
    if ( xi0 < 0.0 ) then
      xi0 = 0.0
      grad = a(2)
      if ( grad == 0.0 ) xi0 = xi0 + eps
    endif
 
    if ( xi0 > 1.0 ) then
      xi0 = 1.0
      if (use_2018_answers) then
        grad = a(2) + 2*a(3) + 3*a(4) + 4*a(5)
      else  ! These expressions are mathematicaly equivalent but more accurate.
        grad = a(2) + (2.*a(3) + (3.*a(4) + 4.*a(5)))
      endif
      if ( grad == 0.0 ) xi0 = xi0 - eps
    endif
 
    ! break if converged or too many iterations taken
    if ( abs(delta) < nr_tolerance ) exit
  enddo ! end Newton-Raphson iterations
 
  x_tgt = x_g(k_found) + xi0 * h(k_found)

interpolate_grid()

subroutine, public regrid_interp::interpolate_grid	(	integer, intent(in)	n0,
		real, dimension(n0), intent(in)	h0,
		real, dimension(n0+1), intent(in)	x0,
		real, dimension(n0,2), intent(in)	ppoly0_E,
		real, dimension(n0,degree_max+1), intent(in)	ppoly0_coefs,
		real, dimension(n1+1), intent(in)	target_values,
		integer, intent(in)	degree,
		integer, intent(in)	n1,
		real, dimension(n1), intent(inout)	h1,
		real, dimension(n1+1), 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]	n0	Number of points on source grid
[in]	n1	Number of points on target grid
[in]	h0	Thicknesses of source grid cells [H]
[in]	x0	Source interface positions [H]
[in]	ppoly0_e	Edge values of interpolating polynomials [A]
[in]	ppoly0_coefs	Coefficients of interpolating polynomials [A]
[in]	target_values	Target values of interfaces [A]
[in]	degree	Degree of interpolating polynomials
[in,out]	h1	Thicknesses of target grid cells [H]
[in,out]	x1	Target interface positions [H]
[in]	answers_2018	If true use older, less acccurate expressions.

Definition at line 272 of file regrid_interp.F90.

  integer,               intent(in)     :: n0            !< Number of points on source grid
  integer,               intent(in)     :: n1            !< Number of points on target grid
  real, dimension(n0),   intent(in)     :: h0            !< Thicknesses of source grid cells [H]
  real, dimension(n0+1), intent(in)     :: x0            !< Source interface positions [H]
  real, dimension(n0,2), intent(in)     :: ppoly0_E      !< Edge values of interpolating polynomials [A]
  real, dimension(n0,DEGREE_MAX+1), &
                          intent(in)    :: ppoly0_coefs  !< Coefficients of interpolating polynomials [A]
  real, dimension(n1+1),  intent(in)    :: target_values !< Target values of interfaces [A]
  integer,                intent(in)    :: degree        !< Degree of interpolating polynomials
  real, dimension(n1),    intent(inout) :: h1            !< Thicknesses of target grid cells [H]
  real, dimension(n1+1),  intent(inout) :: x1            !< Target interface positions [H]
  logical,      optional, intent(in)    :: answers_2018  !< If true use older, less acccurate expressions.
 
  ! Local variables
  logical   :: use_2018_answers  ! If true use older, less acccurate expressions.
  integer        :: k ! loop index
  real           :: t ! current interface target density
 
  ! Make sure boundary coordinates of new grid coincide with boundary
  ! coordinates of previous grid
  x1(1) = x0(1)
  x1(n1+1) = x0(n0+1)
 
  ! Find coordinates for interior target values
  do k = 2,n1
    t = target_values(k)
    x1(k) = get_polynomial_coordinate( n0, h0, x0, ppoly0_e, ppoly0_coefs, t, degree, &
                                        answers_2018=answers_2018 )
    h1(k-1) = x1(k) - x1(k-1)
  enddo
  h1(n1) = x1(n1+1) - x1(n1)
 

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_scheme

Name 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 488 of file regrid_interp.F90.

  character(len=*), intent(in) :: interp_scheme !< Name 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"
 
  select case ( uppercase(trim(interp_scheme)) )
    case ("P1M_H2");     interpolation_scheme = interpolation_p1m_h2
    case ("P1M_H4");     interpolation_scheme = interpolation_p1m_h4
    case ("P1M_IH2");    interpolation_scheme = interpolation_p1m_ih4
    case ("PLM");        interpolation_scheme = interpolation_plm
    case ("PPM_H4");     interpolation_scheme = interpolation_ppm_h4
    case ("PPM_IH4");    interpolation_scheme = interpolation_ppm_ih4
    case ("P3M_IH4IH3"); interpolation_scheme = interpolation_p3m_ih4ih3
    case ("P3M_IH6IH5"); interpolation_scheme = interpolation_p3m_ih6ih5
    case ("PQM_IH4IH3"); interpolation_scheme = interpolation_pqm_ih4ih3
    case ("PQM_IH6IH5"); interpolation_scheme = interpolation_pqm_ih6ih5
    case default ; call mom_error(fatal, "regrid_interp: "//&
     "Unrecognized choice for INTERPOLATION_SCHEME ("//trim(interp_scheme)//").")
  end select

regridding_set_ppolys()

subroutine, public regrid_interp::regridding_set_ppolys	(	type(interp_cs_type), intent(in)	CS,
		real, dimension(n0), intent(in)	densities,
		integer, intent(in)	n0,
		real, dimension(n0), intent(in)	h0,
		real, dimension(n0,2), intent(inout)	ppoly0_E,
		real, dimension(n0,2), intent(inout)	ppoly0_S,
		real, dimension(n0,degree_max+1), 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]	cs	Interpolation control structure
[in]	n0	Number of cells on source grid
[in]	densities	Actual cell densities [A]
[in]	h0	cell widths on source grid [H]
[in,out]	ppoly0_e	Edge value of polynomial [A]
[in,out]	ppoly0_s	Edge slope of polynomial [A H-1]
[in,out]	ppoly0_coefs	Coefficients of polynomial [A]
[in,out]	degree	The degree of the polynomials
[in]	h_neglect	A negligibly small width for the purpose of cell reconstructions [H] in the same units as h0.
[in]	h_neglect_edge	A negligibly small width for the purpose of edge value calculations [H] in the same units as h0.

Definition at line 80 of file regrid_interp.F90.

  type(interp_CS_type),  intent(in)    :: CS !< Interpolation control structure
  integer,               intent(in)    :: n0 !< Number of cells on source grid
  real, dimension(n0),   intent(in)    :: densities !< Actual cell densities [A]
  real, dimension(n0),   intent(in)    :: h0 !< cell widths on source grid [H]
  real, dimension(n0,2), intent(inout) :: ppoly0_E  !< Edge value of polynomial [A]
  real, dimension(n0,2), intent(inout) :: ppoly0_S  !< Edge slope of polynomial [A H-1]
  real, dimension(n0,DEGREE_MAX+1), intent(inout) :: ppoly0_coefs !< Coefficients of polynomial [A]
  integer,               intent(inout) :: degree    !< The degree of the polynomials
  real,        optional, intent(in)    :: h_neglect !< A negligibly small width for the
                                             !! purpose of cell reconstructions [H]
                                             !! in the same units as h0.
  real,        optional, intent(in)    :: h_neglect_edge !< A negligibly small width
                                             !! for the purpose of edge value calculations [H]
                                             !! in the same units as h0.
  ! Local variables
  logical :: extrapolate
 
  ! Reset piecewise polynomials
  ppoly0_e(:,:) = 0.0
  ppoly0_s(:,:) = 0.0
  ppoly0_coefs(:,:) = 0.0
 
  extrapolate = cs%boundary_extrapolation
 
  ! Compute the interpolated profile of the density field and build grid
  select case (cs%interpolation_scheme)
 
    case ( interpolation_p1m_h2 )
      degree = degree_1
      call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
      call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answers_2018=cs%answers_2018 )
      if (extrapolate) then
        call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
      endif
 
    case ( interpolation_p1m_h4 )
      degree = degree_1
      if ( n0 >= 4 ) then
        call edge_values_explicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge, answers_2018=cs%answers_2018 )
      else
        call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
      endif
      call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answers_2018=cs%answers_2018 )
      if (extrapolate) then
        call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
      endif
 
    case ( interpolation_p1m_ih4 )
      degree = degree_1
      if ( n0 >= 4 ) then
        call edge_values_implicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge, answers_2018=cs%answers_2018 )
      else
        call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
      endif
      call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answers_2018=cs%answers_2018 )
      if (extrapolate) then
        call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
      endif
 
    case ( interpolation_plm )
      degree = degree_1
      call plm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
      if (extrapolate) then
        call plm_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
      endif
 
    case ( interpolation_ppm_h4 )
      if ( n0 >= 4 ) then
        degree = degree_2
        call edge_values_explicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge, answers_2018=cs%answers_2018 )
        call ppm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answers_2018=cs%answers_2018 )
        if (extrapolate) then
          call ppm_boundary_extrapolation( n0, h0, densities, ppoly0_e, &
                                           ppoly0_coefs, h_neglect )
        endif
      else
        degree = degree_1
        call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
        call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answers_2018=cs%answers_2018 )
        if (extrapolate) then
          call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
        endif
      endif
 
    case ( interpolation_ppm_ih4 )
      if ( n0 >= 4 ) then
        degree = degree_2
        call edge_values_implicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge, answers_2018=cs%answers_2018 )
        call ppm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answers_2018=cs%answers_2018 )
        if (extrapolate) then
          call ppm_boundary_extrapolation( n0, h0, densities, ppoly0_e, &
                                           ppoly0_coefs, h_neglect )
        endif
      else
        degree = degree_1
        call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
        call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answers_2018=cs%answers_2018 )
        if (extrapolate) then
          call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
        endif
      endif
 
    case ( interpolation_p3m_ih4ih3 )
      if ( n0 >= 4 ) then
        degree = degree_3
        call edge_values_implicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge, answers_2018=cs%answers_2018 )
        call edge_slopes_implicit_h3( n0, h0, densities, ppoly0_s, h_neglect, answers_2018=cs%answers_2018 )
        call p3m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_s, &
                                ppoly0_coefs, h_neglect, answers_2018=cs%answers_2018 )
        if (extrapolate) then
          call p3m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_s, &
                                           ppoly0_coefs, h_neglect, h_neglect_edge )
        endif
      else
        degree = degree_1
        call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
        call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answers_2018=cs%answers_2018 )
        if (extrapolate) then
          call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
        endif
      endif
 
    case ( interpolation_p3m_ih6ih5 )
      if ( n0 >= 6 ) then
        degree = degree_3
        call edge_values_implicit_h6( n0, h0, densities, ppoly0_e, h_neglect_edge, answers_2018=cs%answers_2018 )
        call edge_slopes_implicit_h5( n0, h0, densities, ppoly0_s, h_neglect, answers_2018=cs%answers_2018 )
        call p3m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_s, &
                                ppoly0_coefs, h_neglect, answers_2018=cs%answers_2018 )
        if (extrapolate) then
          call p3m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_s, &
                   ppoly0_coefs, h_neglect, h_neglect_edge )
        endif
      else
        degree = degree_1
        call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
        call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answers_2018=cs%answers_2018 )
        if (extrapolate) then
          call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
        endif
      endif
 
    case ( interpolation_pqm_ih4ih3 )
      if ( n0 >= 4 ) then
        degree = degree_4
        call edge_values_implicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge, answers_2018=cs%answers_2018 )
        call edge_slopes_implicit_h3( n0, h0, densities, ppoly0_s, h_neglect, answers_2018=cs%answers_2018 )
        call pqm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_s, &
                                 ppoly0_coefs, h_neglect, answers_2018=cs%answers_2018 )
        if (extrapolate) then
          call pqm_boundary_extrapolation_v1( n0, h0, densities, ppoly0_e, ppoly0_s, &
                                 ppoly0_coefs, h_neglect )
        endif
      else
        degree = degree_1
        call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
        call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answers_2018=cs%answers_2018 )
        if (extrapolate) then
          call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
        endif
      endif
 
    case ( interpolation_pqm_ih6ih5 )
      if ( n0 >= 6 ) then
        degree = degree_4
        call edge_values_implicit_h6( n0, h0, densities, ppoly0_e, h_neglect_edge, answers_2018=cs%answers_2018 )
        call edge_slopes_implicit_h5( n0, h0, densities, ppoly0_s, h_neglect, answers_2018=cs%answers_2018 )
        call pqm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_s, &
                                 ppoly0_coefs, h_neglect, answers_2018=cs%answers_2018 )
        if (extrapolate) then
          call pqm_boundary_extrapolation_v1( n0, h0, densities, ppoly0_e, ppoly0_s, &
                                 ppoly0_coefs, h_neglect )
        endif
      else
        degree = degree_1
        call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
        call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answers_2018=cs%answers_2018 )
        if (extrapolate) then
          call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
        endif
      endif
  end select
 

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]	cs	A control structure for regrid_interp
[in]	extrap	Indicate whether high-order boundary extrapolation should be used in boundary cells

Definition at line 520 of file regrid_interp.F90.

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

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]	cs	A control structure for regrid_interp
[in]	interp_scheme	Name 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 510 of file regrid_interp.F90.

  type(interp_CS_type), intent(inout) :: CS  !< A control structure for regrid_interp
  character(len=*),     intent(in) :: interp_scheme !< Name 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"
 
  cs%interpolation_scheme = interpolation_scheme(interp_scheme)