regrid_interp.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!> Vertical interpolation for regridding
6module regrid_interp
7
8use mom_error_handler, only : mom_error, fatal
10
15
20
23
24implicit none ; private
25
26!> Control structure for regrid_interp module
27type, public :: interp_cs_type ; private
28
29 !> The following parameter is only relevant when used with the target
30 !! interface densities regridding scheme. It indicates which interpolation
31 !! to use to determine the grid.
32 integer :: interpolation_scheme = -1
33
34 !> Indicate whether high-order boundary extrapolation should be used within
35 !! boundary cells
36 logical :: boundary_extrapolation
37
38 !> The vintage of the expressions to use for regridding
39 integer :: answer_date = 99991231
40end type interp_cs_type
41
44
45! List of interpolation schemes
46integer, parameter :: interpolation_p1m_h2 = 0 !< O(h^2)
47integer, parameter :: interpolation_p1m_h4 = 1 !< O(h^2)
48integer, parameter :: interpolation_p1m_ih4 = 2 !< O(h^2)
49integer, parameter :: interpolation_plm = 3 !< O(h^2)
50integer, parameter :: interpolation_ppm_cw =10 !< O(h^3)
51integer, parameter :: interpolation_ppm_h4 = 4 !< O(h^3)
52integer, parameter :: interpolation_ppm_ih4 = 5 !< O(h^3)
53integer, parameter :: interpolation_p3m_ih4ih3 = 6 !< O(h^4)
54integer, parameter :: interpolation_p3m_ih6ih5 = 7 !< O(h^4)
55integer, parameter :: interpolation_pqm_ih4ih3 = 8 !< O(h^4)
56integer, parameter :: interpolation_pqm_ih6ih5 = 9 !< O(h^5)
57
58!>@{ Interpolant degrees
59integer, parameter :: degree_1 = 1, degree_2 = 2, degree_3 = 3, degree_4 = 4
60integer, public, parameter :: degree_max = 5
61!>@}
62
63!> When the N-R algorithm produces an estimate that lies outside [0,1], the
64!! estimate is set to be equal to the boundary location, 0 or 1, plus or minus
65!! an offset, respectively, when the derivative is zero at the boundary [nondim].
66real, public, parameter :: nr_offset = 1e-6
67!> Maximum number of Newton-Raphson iterations. Newton-Raphson iterations are
68!! used to build the new grid by finding the coordinates associated with
69!! target densities and interpolations of degree larger than 1.
70integer, public, parameter :: nr_iterations = 8
71!> Tolerance for Newton-Raphson iterations (stop when increment falls below this) [nondim]
72real, public, parameter :: nr_tolerance = 1e-12
73
74contains
75
76!> Builds an interpolated profile for the densities within each grid cell.
77!!
78!! It may happen that, given a high-order interpolator, the number of
79!! available layers is insufficient (e.g., there are two available layers for
80!! a third-order PPM ih4 scheme). In these cases, we resort to the simplest
81!! continuous linear scheme (P1M h2).
82subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, &
83 ppoly0_coefs, degree, h_neglect, h_neglect_edge)
84 type(interp_cs_type), intent(in) :: cs !< Interpolation control structure
85 integer, intent(in) :: n0 !< Number of cells on source grid
86 real, dimension(n0), intent(in) :: densities !< Actual cell densities [A]
87 real, dimension(n0), intent(in) :: h0 !< cell widths on source grid [H]
88 real, dimension(n0,2), intent(inout) :: ppoly0_e !< Edge value of polynomial [A]
89 real, dimension(n0,2), intent(inout) :: ppoly0_s !< Edge slope of polynomial [A H-1]
90 real, dimension(n0,DEGREE_MAX+1), intent(inout) :: ppoly0_coefs !< Coefficients of polynomial [A]
91 integer, intent(inout) :: degree !< The degree of the polynomials
92 real, intent(in) :: h_neglect !< A negligibly small width for the
93 !! purpose of cell reconstructions [H]
94 !! in the same units as h0.
95 real, optional, intent(in) :: h_neglect_edge !< A negligibly small width
96 !! for the purpose of edge value calculations [H]
97 !! in the same units as h0.
98 ! Local variables
99 real :: h_neg_edge ! A negligibly small width for the purpose of edge value
100 ! calculations in the same units as h0 [H]
101 logical :: extrapolate
102
103 h_neg_edge = h_neglect ; if (present(h_neglect_edge)) h_neg_edge = h_neglect_edge
104
105 ! Reset piecewise polynomials
106 ppoly0_e(:,:) = 0.0
107 ppoly0_s(:,:) = 0.0
108 ppoly0_coefs(:,:) = 0.0
109
110 extrapolate = cs%boundary_extrapolation
111
112 ! Compute the interpolated profile of the density field and build grid
113 select case (cs%interpolation_scheme)
114
115 case ( interpolation_p1m_h2 )
116 degree = degree_1
117 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
118 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
119 if (extrapolate) then
120 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
121 endif
122
123 case ( interpolation_p1m_h4 )
124 degree = degree_1
125 if ( n0 >= 4 ) then
126 call edge_values_explicit_h4( n0, h0, densities, ppoly0_e, h_neg_edge, answer_date=cs%answer_date )
127 else
128 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
129 endif
130 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
131 if (extrapolate) then
132 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
133 endif
134
135 case ( interpolation_p1m_ih4 )
136 degree = degree_1
137 if ( n0 >= 4 ) then
138 call edge_values_implicit_h4( n0, h0, densities, ppoly0_e, h_neg_edge, answer_date=cs%answer_date )
139 else
140 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
141 endif
142 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
143 if (extrapolate) then
144 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
145 endif
146
147 case ( interpolation_plm )
148 degree = degree_1
149 call plm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
150 if (extrapolate) then
151 call plm_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
152 endif
153
154 case ( interpolation_ppm_cw )
155 if ( n0 >= 4 ) then
156 degree = degree_2
157 call edge_values_explicit_h4cw( n0, h0, densities, ppoly0_e, h_neg_edge )
158 call ppm_monotonicity( n0, densities, ppoly0_e )
159 call ppm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
160 if (extrapolate) then
161 call ppm_boundary_extrapolation( n0, h0, densities, ppoly0_e, &
162 ppoly0_coefs, h_neglect )
163 endif
164 else
165 degree = degree_1
166 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
167 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
168 if (extrapolate) then
169 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
170 endif
171 endif
172
173 case ( interpolation_ppm_h4 )
174 if ( n0 >= 4 ) then
175 degree = degree_2
176 call edge_values_explicit_h4( n0, h0, densities, ppoly0_e, h_neg_edge, answer_date=cs%answer_date )
177 call ppm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
178 if (extrapolate) then
179 call ppm_boundary_extrapolation( n0, h0, densities, ppoly0_e, &
180 ppoly0_coefs, h_neglect )
181 endif
182 else
183 degree = degree_1
184 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
185 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
186 if (extrapolate) then
187 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
188 endif
189 endif
190
191 case ( interpolation_ppm_ih4 )
192 if ( n0 >= 4 ) then
193 degree = degree_2
194 call edge_values_implicit_h4( n0, h0, densities, ppoly0_e, h_neg_edge, answer_date=cs%answer_date )
195 call ppm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
196 if (extrapolate) then
197 call ppm_boundary_extrapolation( n0, h0, densities, ppoly0_e, &
198 ppoly0_coefs, h_neglect )
199 endif
200 else
201 degree = degree_1
202 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
203 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
204 if (extrapolate) then
205 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
206 endif
207 endif
208
209 case ( interpolation_p3m_ih4ih3 )
210 if ( n0 >= 4 ) then
211 degree = degree_3
212 call edge_values_implicit_h4( n0, h0, densities, ppoly0_e, h_neg_edge, answer_date=cs%answer_date )
213 call edge_slopes_implicit_h3( n0, h0, densities, ppoly0_s, h_neglect, answer_date=cs%answer_date )
214 call p3m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_s, &
215 ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
216 if (extrapolate) then
217 call p3m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_s, &
218 ppoly0_coefs, h_neglect, h_neg_edge )
219 endif
220 else
221 degree = degree_1
222 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
223 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
224 if (extrapolate) then
225 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
226 endif
227 endif
228
229 case ( interpolation_p3m_ih6ih5 )
230 if ( n0 >= 6 ) then
231 degree = degree_3
232 call edge_values_implicit_h6( n0, h0, densities, ppoly0_e, h_neg_edge, answer_date=cs%answer_date )
233 call edge_slopes_implicit_h5( n0, h0, densities, ppoly0_s, h_neglect, answer_date=cs%answer_date )
234 call p3m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_s, &
235 ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
236 if (extrapolate) then
237 call p3m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_s, &
238 ppoly0_coefs, h_neglect, h_neglect_edge )
239 endif
240 else
241 degree = degree_1
242 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
243 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
244 if (extrapolate) then
245 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
246 endif
247 endif
248
249 case ( interpolation_pqm_ih4ih3 )
250 if ( n0 >= 4 ) then
251 degree = degree_4
252 call edge_values_implicit_h4( n0, h0, densities, ppoly0_e, h_neg_edge, answer_date=cs%answer_date )
253 call edge_slopes_implicit_h3( n0, h0, densities, ppoly0_s, h_neglect, answer_date=cs%answer_date )
254 call pqm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_s, &
255 ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
256 if (extrapolate) then
257 call pqm_boundary_extrapolation_v1( n0, h0, densities, ppoly0_e, ppoly0_s, &
258 ppoly0_coefs, h_neglect )
259 endif
260 else
261 degree = degree_1
262 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
263 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
264 if (extrapolate) then
265 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
266 endif
267 endif
268
269 case ( interpolation_pqm_ih6ih5 )
270 if ( n0 >= 6 ) then
271 degree = degree_4
272 call edge_values_implicit_h6( n0, h0, densities, ppoly0_e, h_neg_edge, answer_date=cs%answer_date )
273 call edge_slopes_implicit_h5( n0, h0, densities, ppoly0_s, h_neglect, answer_date=cs%answer_date )
274 call pqm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_s, &
275 ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
276 if (extrapolate) then
277 call pqm_boundary_extrapolation_v1( n0, h0, densities, ppoly0_e, ppoly0_s, &
278 ppoly0_coefs, h_neglect )
279 endif
280 else
281 degree = degree_1
282 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e )
283 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect, answer_date=cs%answer_date )
284 if (extrapolate) then
285 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
286 endif
287 endif
288 end select
289
290end subroutine regridding_set_ppolys
291
292!> Given target values (e.g., density), build new grid based on polynomial
293!!
294!! Given the grid 'grid0' and the piecewise polynomial interpolant
295!! 'ppoly0' (possibly discontinuous), the coordinates of the new grid 'grid1'
296!! are determined by finding the corresponding target interface densities.
297subroutine interpolate_grid( n0, h0, x0, ppoly0_E, ppoly0_coefs, &
298 target_values, degree, n1, h1, x1, answer_date )
299 integer, intent(in) :: n0 !< Number of points on source grid
300 integer, intent(in) :: n1 !< Number of points on target grid
301 real, dimension(n0), intent(in) :: h0 !< Thicknesses of source grid cells [H]
302 real, dimension(n0+1), intent(in) :: x0 !< Source interface positions [H]
303 real, dimension(n0,2), intent(in) :: ppoly0_E !< Edge values of interpolating polynomials [A]
304 real, dimension(n0,DEGREE_MAX+1), &
305 intent(in) :: ppoly0_coefs !< Coefficients of interpolating polynomials [A]
306 real, dimension(n1+1), intent(in) :: target_values !< Target values of interfaces [A]
307 integer, intent(in) :: degree !< Degree of interpolating polynomials
308 real, dimension(n1), intent(inout) :: h1 !< Thicknesses of target grid cells [H]
309 real, dimension(n1+1), intent(inout) :: x1 !< Target interface positions [H]
310 integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use
311
312 ! Local variables
313 integer :: k ! loop index
314 real :: t ! current interface target density [A]
315
316 ! Make sure boundary coordinates of new grid coincide with boundary
317 ! coordinates of previous grid
318 x1(1) = x0(1)
319 x1(n1+1) = x0(n0+1)
320
321 ! Find coordinates for interior target values
322 do k = 2,n1
323 t = target_values(k)
324 x1(k) = get_polynomial_coordinate( n0, h0, x0, ppoly0_e, ppoly0_coefs, t, degree, &
325 answer_date=answer_date )
326 h1(k-1) = x1(k) - x1(k-1)
327 enddo
328 h1(n1) = x1(n1+1) - x1(n1)
329
330end subroutine interpolate_grid
331
332!> Build a grid by interpolating for target values
333subroutine build_and_interpolate_grid(CS, densities, n0, h0, x0, target_values, &
334 n1, h1, x1, h_neglect, h_neglect_edge)
335 type(interp_cs_type), intent(in) :: cs !< A control structure for regrid_interp
336 integer, intent(in) :: n0 !< The number of points on the input grid
337 integer, intent(in) :: n1 !< The number of points on the output grid
338 real, dimension(n0), intent(in) :: densities !< Input cell densities [R ~> kg m-3]
339 real, dimension(n1+1), intent(in) :: target_values !< Target values of interfaces [R ~> kg m-3]
340 real, dimension(n0), intent(in) :: h0 !< Initial cell widths usually in [H ~> m or kg m-2] or [Z ~> m]
341 real, dimension(n0+1), intent(in) :: x0 !< Source interface positions [H ~> m or kg m-2] or [Z ~> m]
342 real, dimension(n1), intent(inout) :: h1 !< Output cell widths [H ~> m or kg m-2] or [Z ~> m]
343 real, dimension(n1+1), intent(inout) :: x1 !< Target interface positions [H ~> m or kg m-2] or [Z ~> m]
344 real, intent(in) :: h_neglect !< A negligibly small width for the
345 !! purpose of cell reconstructions in the same
346 !! units as h0 [H ~> m or kg m-2] or [Z ~> m].
347 real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the
348 !! purpose of edge value calculations in the same
349 !! units as h0 [H ~> m or kg m-2] or [Z ~> m]
350
351 real, dimension(n0,2) :: ppoly0_e ! Polynomial edge values [R ~> kg m-3]
352 real, dimension(n0,2) :: ppoly0_s ! Polynomial edge slopes [R H-1 ~> kg m-4 or m-1] or [R Z-1 ~> kg m-4]
353 real, dimension(n0,DEGREE_MAX+1) :: ppoly0_c ! Polynomial interpolant coeficients on the local 0-1 grid [R ~> kg m-3]
354 integer :: degree
355
356 call regridding_set_ppolys(cs, densities, n0, h0, ppoly0_e, ppoly0_s, ppoly0_c, &
357 degree, h_neglect, h_neglect_edge)
358 call interpolate_grid(n0, h0, x0, ppoly0_e, ppoly0_c, target_values, degree, &
359 n1, h1, x1, answer_date=cs%answer_date)
360end subroutine build_and_interpolate_grid
361
362!> Given a target value, find corresponding coordinate for given polynomial
363!!
364!! Here, 'ppoly' is assumed to be a piecewise discontinuous polynomial of degree
365!! 'degree' throughout the domain defined by 'grid'. A target value is given
366!! and we need to determine the corresponding grid coordinate to define the
367!! new grid.
368!!
369!! If the target value is out of range, the grid coordinate is simply set to
370!! be equal to one of the boundary coordinates, which results in vanished layers
371!! near the boundaries.
372!!
373!! IT IS ASSUMED THAT THE PIECEWISE POLYNOMIAL IS MONOTONICALLY INCREASING.
374!! IF THIS IS NOT THE CASE, THE NEW GRID MAY BE ILL-DEFINED.
375!!
376!! It is assumed that the number of cells defining 'grid' and 'ppoly' are the
377!! same.
378function get_polynomial_coordinate( N, h, x_g, edge_values, ppoly_coefs, &
379 target_value, degree, answer_date ) result ( x_tgt )
380 ! Arguments
381 integer, intent(in) :: n !< Number of grid cells
382 real, dimension(N), intent(in) :: h !< Grid cell thicknesses [H]
383 real, dimension(N+1), intent(in) :: x_g !< Grid interface locations [H]
384 real, dimension(N,2), intent(in) :: edge_values !< Edge values of interpolating polynomials [A]
385 real, dimension(N,DEGREE_MAX+1), intent(in) :: ppoly_coefs !< Coefficients of interpolating polynomials [A]
386 real, intent(in) :: target_value !< Target value to find position for [A]
387 integer, intent(in) :: degree !< Degree of the interpolating polynomials
388 integer, intent(in) :: answer_date !< The vintage of the expressions to use
389 real :: x_tgt !< The position of x_g at which target_value is found [H]
390
391 ! Local variables
392 real :: xi0 ! normalized target coordinate [nondim]
393 real, dimension(DEGREE_MAX) :: a ! polynomial coefficients [A]
394 real :: numerator ! The numerator of an expression [A]
395 real :: denominator ! The denominator of an expression [A]
396 real :: delta ! Newton-Raphson increment [nondim]
397! real :: x ! global target coordinate [nondim]
398 real :: eps ! offset used to get away from boundaries [nondim]
399 real :: grad ! gradient during N-R iterations [A]
400 integer :: i, k, iter ! loop indices
401 integer :: k_found ! index of target cell
402 character(len=320) :: mesg
403 logical :: use_2018_answers ! If true use older, less accurate expressions.
404
405 eps = nr_offset
406 k_found = -1
407 use_2018_answers = (answer_date < 20190101)
408
409 ! If the target value is outside the range of all values, we
410 ! force the target coordinate to be equal to the lowest or
411 ! largest value, depending on which bound is overtaken
412 if ( target_value <= edge_values(1,1) ) then
413 x_tgt = x_g(1)
414 return ! return because there is no need to look further
415 endif
416
417 ! Since discontinuous edge values are allowed, we check whether the target
418 ! value lies between two discontinuous edge values at interior interfaces
419 do k = 2,n
420 if ( ( target_value >= edge_values(k-1,2) ) .AND. ( target_value <= edge_values(k,1) ) ) then
421 x_tgt = x_g(k)
422 return ! return because there is no need to look further
423 endif
424 enddo
425
426 ! If the target value is outside the range of all values, we
427 ! force the target coordinate to be equal to the lowest or
428 ! largest value, depending on which bound is overtaken
429 if ( target_value >= edge_values(n,2) ) then
430 x_tgt = x_g(n+1)
431 return ! return because there is no need to look further
432 endif
433
434 ! At this point, we know that the target value is bounded and does not
435 ! lie between discontinuous, monotonic edge values. Therefore,
436 ! there is a unique solution. We loop on all cells and find which one
437 ! contains the target value. The variable k_found holds the index value
438 ! of the cell where the taregt value lies.
439 do k = 1,n
440 if ( ( target_value > edge_values(k,1) ) .AND. ( target_value < edge_values(k,2) ) ) then
441 k_found = k
442 exit
443 endif
444 enddo
445
446 ! At this point, 'k_found' should be strictly positive. If not, this is
447 ! a major failure because it means we could not find any target cell
448 ! despite the fact that the target value lies between the extremes. It
449 ! means there is a major problem with the interpolant. This needs to be
450 ! reported.
451 if ( k_found == -1 ) then
452 write(mesg,*) 'Could not find target coordinate', target_value, 'in get_polynomial_coordinate. This is '//&
453 'caused by an inconsistent interpolant (perhaps not monotonically increasing):', &
454 target_value, edge_values(1,1), edge_values(n,2)
455 call mom_error( fatal, mesg )
456 endif
457
458 ! Reset all polynomial coefficients to 0 and copy those pertaining to
459 ! the found cell
460 a(:) = 0.0
461 do i = 1,degree+1
462 a(i) = ppoly_coefs(k_found,i)
463 enddo
464
465 ! Guess the middle of the cell to start Newton-Raphson iterations
466 xi0 = 0.5
467
468 ! Newton-Raphson iterations
469 do iter = 1,nr_iterations
470
471 if (use_2018_answers) then
472 numerator = a(1) + a(2)*xi0 + a(3)*xi0*xi0 + a(4)*xi0*xi0*xi0 + &
473 a(5)*xi0*xi0*xi0*xi0 - target_value
474 denominator = a(2) + 2*a(3)*xi0 + 3*a(4)*xi0*xi0 + 4*a(5)*xi0*xi0*xi0
475 else ! These expressions are mathematicaly equivalent but more accurate.
476 numerator = (a(1) - target_value) + xi0*(a(2) + xi0*(a(3) + xi0*(a(4) + a(5)*xi0)))
477 denominator = a(2) + xi0*(2.*a(3) + xi0*(3.*a(4) + 4.*a(5)*xi0))
478 endif
479
480 delta = -numerator / denominator
481
482 xi0 = xi0 + delta
483
484 ! Check whether new estimate is out of bounds. If the new estimate is
485 ! indeed out of bounds, we manually set it to be equal to the overtaken
486 ! bound with a small offset towards the interior when the gradient of
487 ! the function at the boundary is zero (in which case, the Newton-Raphson
488 ! algorithm does not converge).
489 if ( xi0 < 0.0 ) then
490 xi0 = 0.0
491 grad = a(2)
492 if ( grad == 0.0 ) xi0 = xi0 + eps
493 endif
494
495 if ( xi0 > 1.0 ) then
496 xi0 = 1.0
497 if (use_2018_answers) then
498 grad = a(2) + 2*a(3) + 3*a(4) + 4*a(5)
499 else ! These expressions are mathematicaly equivalent but more accurate.
500 grad = a(2) + (2.*a(3) + (3.*a(4) + 4.*a(5)))
501 endif
502 if ( grad == 0.0 ) xi0 = xi0 - eps
503 endif
504
505 ! break if converged or too many iterations taken
506 if ( abs(delta) < nr_tolerance ) exit
507 enddo ! end Newton-Raphson iterations
508
509 x_tgt = x_g(k_found) + xi0 * h(k_found)
510end function get_polynomial_coordinate
511
512!> Numeric value of interpolation_scheme corresponding to scheme name
513integer function interpolation_scheme(interp_scheme)
514 character(len=*), intent(in) :: interp_scheme !< Name of the interpolation scheme
515 !! Valid values include "P1M_H2", "P1M_H4", "P1M_IH2", "PLM", "PPM_CW", "PPM_H4",
516 !! "PPM_IH4", "P3M_IH4IH3", "P3M_IH6IH5", "PQM_IH4IH3", and "PQM_IH6IH5"
517
518 select case ( uppercase(trim(interp_scheme)) )
519 case ("P1M_H2"); interpolation_scheme = interpolation_p1m_h2
520 case ("P1M_H4"); interpolation_scheme = interpolation_p1m_h4
521 case ("P1M_IH2"); interpolation_scheme = interpolation_p1m_ih4
522 case ("PLM"); interpolation_scheme = interpolation_plm
523 case ("PPM_CW"); interpolation_scheme = interpolation_ppm_cw
524 case ("PPM_H4"); interpolation_scheme = interpolation_ppm_h4
525 case ("PPM_IH4"); interpolation_scheme = interpolation_ppm_ih4
526 case ("P3M_IH4IH3"); interpolation_scheme = interpolation_p3m_ih4ih3
527 case ("P3M_IH6IH5"); interpolation_scheme = interpolation_p3m_ih6ih5
528 case ("PQM_IH4IH3"); interpolation_scheme = interpolation_pqm_ih4ih3
529 case ("PQM_IH6IH5"); interpolation_scheme = interpolation_pqm_ih6ih5
530 case default ; call mom_error(fatal, "regrid_interp: "//&
531 "Unrecognized choice for INTERPOLATION_SCHEME ("//trim(interp_scheme)//").")
532 end select
533end function interpolation_scheme
534
535!> Store the interpolation_scheme value in the interp_CS based on the input string.
536subroutine set_interp_scheme(CS, interp_scheme)
537 type(interp_cs_type), intent(inout) :: cs !< A control structure for regrid_interp
538 character(len=*), intent(in) :: interp_scheme !< Name of the interpolation scheme
539 !! Valid values include "P1M_H2", "P1M_H4", "P1M_IH2", "PLM", "PPM_CW", "PPM_H4",
540 !! "PPM_IH4", "P3M_IH4IH3", "P3M_IH6IH5", "PQM_IH4IH3", and "PQM_IH6IH5"
541
542 cs%interpolation_scheme = interpolation_scheme(interp_scheme)
543end subroutine set_interp_scheme
544
545!> Store the boundary_extrapolation value in the interp_CS
546subroutine set_interp_extrap(CS, extrap)
547 type(interp_cs_type), intent(inout) :: cs !< A control structure for regrid_interp
548 logical, intent(in) :: extrap !< Indicate whether high-order boundary
549 !! extrapolation should be used in boundary cells
550
551 cs%boundary_extrapolation = extrap
552end subroutine set_interp_extrap
553
554!> Store the value of the answer_date in the interp_CS
555subroutine set_interp_answer_date(CS, answer_date)
556 type(interp_cs_type), intent(inout) :: cs !< A control structure for regrid_interp
557 integer, intent(in) :: answer_date !< An integer encoding the vintage of
558 !! the expressions to use for regridding
559
560 cs%answer_date = answer_date
561end subroutine set_interp_answer_date
562
563end module regrid_interp