P1M_functions.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Linear interpolation functions
6module p1m_functions
7
9
10implicit none ; private
11
12! The following routines are visible to the outside world
14
15contains
16
17!> Linearly interpolate between edge values
18!!
19!! The resulting piecewise interpolant is stored in 'ppoly'.
20!! See 'ppoly.F90' for a definition of this structure.
21!!
22!! The edge values MUST have been estimated prior to calling this routine.
23!!
24!! The estimated edge values must be limited to ensure monotonicity of the
25!! interpolant. We also make sure that edge values are NOT discontinuous.
26!!
27!! It is assumed that the size of the array 'u' is equal to the number of cells
28!! defining 'grid' and 'ppoly'. No consistency check is performed here.
29subroutine p1m_interpolation( N, h, u, edge_values, ppoly_coef, h_neglect, answer_date )
30 integer, intent(in) :: n !< Number of cells
31 real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
32 real, dimension(:), intent(in) :: u !< cell average properties (size N) [A]
33 real, dimension(:,:), intent(inout) :: edge_values !< Potentially modified edge values [A]
34 real, dimension(:,:), intent(inout) :: ppoly_coef !< Potentially modified
35 !! piecewise polynomial coefficients, mainly [A]
36 real, intent(in) :: h_neglect !< A negligibly small width [H]
37 integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use
38
39 ! Local variables
40 integer :: k ! loop index
41 real :: u0_l, u0_r ! edge values (left and right) [A]
42
43 ! Bound edge values (routine found in 'edge_values.F90')
44 call bound_edge_values( n, h, u, edge_values, h_neglect, answer_date=answer_date )
45
46 ! Systematically average discontinuous edge values (routine found in
47 ! 'edge_values.F90')
48 call average_discontinuous_edge_values( n, edge_values )
49
50 ! Loop on interior cells to build interpolants
51 do k = 1,n
52
53 u0_l = edge_values(k,1)
54 u0_r = edge_values(k,2)
55
56 ppoly_coef(k,1) = u0_l
57 ppoly_coef(k,2) = u0_r - u0_l
58
59 enddo ! end loop on interior cells
60
61end subroutine p1m_interpolation
62
63!> Interpolation by linear polynomials within boundary cells
64!!
65!! The left and right edge values in the left and right boundary cells,
66!! respectively, are estimated using a linear extrapolation within the cells.
67!!
68!! It is assumed that the size of the array 'u' is equal to the number of cells
69!! defining 'grid' and 'ppoly'. No consistency check is performed here.
70subroutine p1m_boundary_extrapolation( N, h, u, edge_values, ppoly_coef )
71 ! Arguments
72 integer, intent(in) :: n !< Number of cells
73 real, dimension(:), intent(in) :: h !< cell widths (size N) [H]
74 real, dimension(:), intent(in) :: u !< cell averages (size N) [A]
75 real, dimension(:,:), intent(inout) :: edge_values !< edge values of piecewise polynomials [A]
76 real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly [A]
77
78 ! Local variables
79 real :: u0, u1 ! cell averages [A]
80 real :: h0, h1 ! corresponding cell widths [H]
81 real :: slope ! retained PLM slope [A]
82 real :: u0_l, u0_r ! edge values [A]
83
84 ! -----------------------------------------
85 ! Left edge value in the left boundary cell
86 ! -----------------------------------------
87 h0 = h(1)
88 h1 = h(2)
89
90 u0 = u(1)
91 u1 = u(2)
92
93 ! The standard PLM slope is computed as a first estimate for the
94 ! interpolation within the cell
95 slope = 2.0 * ( u1 - u0 )
96
97 ! The right edge value is then computed and we check whether this
98 ! right edge value is consistent: it cannot be larger than the edge
99 ! value in the neighboring cell if the data set is increasing.
100 ! If the right value is found to too large, the slope is further limited
101 ! by using the edge value in the neighboring cell.
102 u0_r = u0 + 0.5 * slope
103
104 if ( (u1 - u0) * (edge_values(2,1) - u0_r) < 0.0 ) then
105 slope = 2.0 * ( edge_values(2,1) - u0 )
106 endif
107
108 ! Using the limited slope, the left edge value is reevaluated and
109 ! the interpolant coefficients recomputed
110 if ( h0 /= 0.0 ) then
111 edge_values(1,1) = u0 - 0.5 * slope
112 else
113 edge_values(1,1) = u0
114 endif
115
116 ppoly_coef(1,1) = edge_values(1,1)
117 ppoly_coef(1,2) = edge_values(1,2) - edge_values(1,1)
118
119 ! ------------------------------------------
120 ! Right edge value in the left boundary cell
121 ! ------------------------------------------
122 h0 = h(n-1)
123 h1 = h(n)
124
125 u0 = u(n-1)
126 u1 = u(n)
127
128 slope = 2.0 * ( u1 - u0 )
129
130 u0_l = u1 - 0.5 * slope
131
132 if ( (u1 - u0) * (u0_l - edge_values(n-1,2)) < 0.0 ) then
133 slope = 2.0 * ( u1 - edge_values(n-1,2) )
134 endif
135
136 if ( h1 /= 0.0 ) then
137 edge_values(n,2) = u1 + 0.5 * slope
138 else
139 edge_values(n,2) = u1
140 endif
141
142 ppoly_coef(n,1) = edge_values(n,1)
143 ppoly_coef(n,2) = edge_values(n,2) - edge_values(n,1)
144
145end subroutine p1m_boundary_extrapolation
146
147!> \namespace p1m_functions
148!!
149!! Date of creation: 2008.06.09
150!! L. White
151!!
152!! This module contains p1m (linear) interpolation routines.
153!!
154!! p1m interpolation is performed by estimating the edge values and
155!! linearly interpolating between them.
156!
157!! Once the edge values are estimated, the limiting process takes care of
158!! ensuring that (1) edge values are bounded by neighboring cell averages
159!! and (2) discontinuous edge values are averaged in order to provide a
160!! fully continuous interpolant throughout the domain. This last step is
161!! essential for the regridding problem to yield a unique solution.
162!! Also, a routine is provided that takes care of linear extrapolation
163!! within the boundary cells.
164
165end module p1m_functions