← back to index

src/ALE/regrid_interp.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

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
9use MOM_string_functions, only : uppercase
10
11use regrid_edge_values, only : edge_values_explicit_h2, edge_values_explicit_h4
12use regrid_edge_values, only : edge_values_explicit_h4cw
13use regrid_edge_values, only : edge_values_implicit_h4, edge_values_implicit_h6
14use regrid_edge_values, only : edge_slopes_implicit_h3, edge_slopes_implicit_h5
15
16use PLM_functions, only : PLM_reconstruction, PLM_boundary_extrapolation
17use PPM_functions, only : PPM_reconstruction, PPM_boundary_extrapolation
18use PPM_functions, only : PPM_monotonicity
19use PQM_functions, only : PQM_reconstruction, PQM_boundary_extrapolation_v1
20
21use P1M_functions, only : P1M_interpolation, P1M_boundary_extrapolation
22use P3M_functions, only : P3M_interpolation, P3M_boundary_extrapolation
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
42public regridding_set_ppolys, build_and_interpolate_grid
43public set_interp_scheme, set_interp_extrap, set_interp_answer_date
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).
820subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, &
830 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
1030 h_neg_edge = h_neglect ; if (present(h_neglect_edge)) h_neg_edge = h_neglect_edge
104
105 ! Reset piecewise polynomials
1060 ppoly0_E(:,:) = 0.0
1070 ppoly0_S(:,:) = 0.0
1080 ppoly0_coefs(:,:) = 0.0
109
1100 extrapolate = CS%boundary_extrapolation
111
112 ! Compute the interpolated profile of the density field and build grid
1130 select case (CS%interpolation_scheme)
114
115 case ( INTERPOLATION_P1M_H2 )
1160 degree = DEGREE_1
1170 call edge_values_explicit_h2( n0, h0, densities, ppoly0_E )
1180 call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
1190 if (extrapolate) then
1200 call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs )
121 endif
122
123 case ( INTERPOLATION_P1M_H4 )
1240 degree = DEGREE_1
1250 if ( n0 >= 4 ) then
1260 call edge_values_explicit_h4( n0, h0, densities, ppoly0_E, h_neg_edge, answer_date=CS%answer_date )
127 else
1280 call edge_values_explicit_h2( n0, h0, densities, ppoly0_E )
129 endif
1300 call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
1310 if (extrapolate) then
1320 call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs )
133 endif
134
135 case ( INTERPOLATION_P1M_IH4 )
1360 degree = DEGREE_1
1370 if ( n0 >= 4 ) then
1380 call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neg_edge, answer_date=CS%answer_date )
139 else
1400 call edge_values_explicit_h2( n0, h0, densities, ppoly0_E )
141 endif
1420 call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
1430 if (extrapolate) then
1440 call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs )
145 endif
146
147 case ( INTERPOLATION_PLM )
1480 degree = DEGREE_1
1490 call PLM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect )
1500 if (extrapolate) then
1510 call PLM_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect )
152 endif
153
154 case ( INTERPOLATION_PPM_CW )
1550 if ( n0 >= 4 ) then
1560 degree = DEGREE_2
1570 call edge_values_explicit_h4cw( n0, h0, densities, ppoly0_E, h_neg_edge )
1580 call PPM_monotonicity( n0, densities, ppoly0_E )
1590 call PPM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
1600 if (extrapolate) then
161 call PPM_boundary_extrapolation( n0, h0, densities, ppoly0_E, &
1620 ppoly0_coefs, h_neglect )
163 endif
164 else
1650 degree = DEGREE_1
1660 call edge_values_explicit_h2( n0, h0, densities, ppoly0_E )
1670 call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
1680 if (extrapolate) then
1690 call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs )
170 endif
171 endif
172
173 case ( INTERPOLATION_PPM_H4 )
1740 if ( n0 >= 4 ) then
1750 degree = DEGREE_2
1760 call edge_values_explicit_h4( n0, h0, densities, ppoly0_E, h_neg_edge, answer_date=CS%answer_date )
1770 call PPM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
1780 if (extrapolate) then
179 call PPM_boundary_extrapolation( n0, h0, densities, ppoly0_E, &
1800 ppoly0_coefs, h_neglect )
181 endif
182 else
1830 degree = DEGREE_1
1840 call edge_values_explicit_h2( n0, h0, densities, ppoly0_E )
1850 call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
1860 if (extrapolate) then
1870 call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs )
188 endif
189 endif
190
191 case ( INTERPOLATION_PPM_IH4 )
1920 if ( n0 >= 4 ) then
1930 degree = DEGREE_2
1940 call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neg_edge, answer_date=CS%answer_date )
1950 call PPM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
1960 if (extrapolate) then
197 call PPM_boundary_extrapolation( n0, h0, densities, ppoly0_E, &
1980 ppoly0_coefs, h_neglect )
199 endif
200 else
2010 degree = DEGREE_1
2020 call edge_values_explicit_h2( n0, h0, densities, ppoly0_E )
2030 call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
2040 if (extrapolate) then
2050 call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs )
206 endif
207 endif
208
209 case ( INTERPOLATION_P3M_IH4IH3 )
2100 if ( n0 >= 4 ) then
2110 degree = DEGREE_3
2120 call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neg_edge, answer_date=CS%answer_date )
2130 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, &
2150 ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
2160 if (extrapolate) then
217 call P3M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_S, &
2180 ppoly0_coefs, h_neglect, h_neg_edge )
219 endif
220 else
2210 degree = DEGREE_1
2220 call edge_values_explicit_h2( n0, h0, densities, ppoly0_E )
2230 call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
2240 if (extrapolate) then
2250 call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs )
226 endif
227 endif
228
229 case ( INTERPOLATION_P3M_IH6IH5 )
2300 if ( n0 >= 6 ) then
2310 degree = DEGREE_3
2320 call edge_values_implicit_h6( n0, h0, densities, ppoly0_E, h_neg_edge, answer_date=CS%answer_date )
2330 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, &
2350 ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
2360 if (extrapolate) then
237 call P3M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_S, &
2380 ppoly0_coefs, h_neglect, h_neglect_edge )
239 endif
240 else
2410 degree = DEGREE_1
2420 call edge_values_explicit_h2( n0, h0, densities, ppoly0_E )
2430 call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
2440 if (extrapolate) then
2450 call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs )
246 endif
247 endif
248
249 case ( INTERPOLATION_PQM_IH4IH3 )
2500 if ( n0 >= 4 ) then
2510 degree = DEGREE_4
2520 call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neg_edge, answer_date=CS%answer_date )
2530 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, &
2550 ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
2560 if (extrapolate) then
257 call PQM_boundary_extrapolation_v1( n0, h0, densities, ppoly0_E, ppoly0_S, &
2580 ppoly0_coefs, h_neglect )
259 endif
260 else
2610 degree = DEGREE_1
2620 call edge_values_explicit_h2( n0, h0, densities, ppoly0_E )
2630 call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
2640 if (extrapolate) then
2650 call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs )
266 endif
267 endif
268
269 case ( INTERPOLATION_PQM_IH6IH5 )
2700 if ( n0 >= 6 ) then
2710 degree = DEGREE_4
2720 call edge_values_implicit_h6( n0, h0, densities, ppoly0_E, h_neg_edge, answer_date=CS%answer_date )
2730 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, &
2750 ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
2760 if (extrapolate) then
277 call PQM_boundary_extrapolation_v1( n0, h0, densities, ppoly0_E, ppoly0_S, &
2780 ppoly0_coefs, h_neglect )
279 endif
280 else
2810 degree = DEGREE_1
2820 call edge_values_explicit_h2( n0, h0, densities, ppoly0_E )
2830 call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date )
2840 if (extrapolate) then
2850 call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefs )
286 endif
287 endif
288 end select
289
2900end 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.
2970subroutine interpolate_grid( n0, h0, x0, ppoly0_E, ppoly0_coefs, &
2980 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
3180 x1(1) = x0(1)
3190 x1(n1+1) = x0(n0+1)
320
321 ! Find coordinates for interior target values
3220 do k = 2,n1
3230 t = target_values(k)
324 x1(k) = get_polynomial_coordinate ( n0, h0, x0, ppoly0_E, ppoly0_coefs, t, degree, &
3250 answer_date=answer_date )
3260 h1(k-1) = x1(k) - x1(k-1)
327 enddo
3280 h1(n1) = x1(n1+1) - x1(n1)
329
3300end subroutine interpolate_grid
331
332!> Build a grid by interpolating for target values
3330subroutine build_and_interpolate_grid(CS, densities, n0, h0, x0, target_values, &
3340 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
3510 real, dimension(n0,2) :: ppoly0_E ! Polynomial edge values [R ~> kg m-3]
3520 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]
3530 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, &
3570 degree, h_neglect, h_neglect_edge)
358 call interpolate_grid(n0, h0, x0, ppoly0_E, ppoly0_C, target_values, degree, &
3590 n1, h1, x1, answer_date=CS%answer_date)
3600end 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.
3780function 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
4050 eps = NR_OFFSET
4060 k_found = -1
4070 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
4120 if ( target_value <= edge_values(1,1) ) then
4130 x_tgt = x_g(1)
4140 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
4190 do k = 2,N
4200 if ( ( target_value >= edge_values(k-1,2) ) .AND. ( target_value <= edge_values(k,1) ) ) then
4210 x_tgt = x_g(k)
4220 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
4290 if ( target_value >= edge_values(N,2) ) then
4300 x_tgt = x_g(N+1)
4310 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.
4390 do k = 1,N
4400 if ( ( target_value > edge_values(k,1) ) .AND. ( target_value < edge_values(k,2) ) ) then
4410 k_found = k
4420 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.
4510 if ( k_found == -1 ) then
4520 write(mesg,*) 'Could not find target coordinate', target_value, 'in get_polynomial_coordinate. This is '//&
4530 'caused by an inconsistent interpolant (perhaps not monotonically increasing):', &
4540 target_value, edge_values(1,1), edge_values(N,2)
4550 call MOM_error( FATAL, mesg )
456 endif
457
458 ! Reset all polynomial coefficients to 0 and copy those pertaining to
459 ! the found cell
4600 a(:) = 0.0
4610 do i = 1,degree+1
4620 a(i) = ppoly_coefs(k_found,i)
463 enddo
464
465 ! Guess the middle of the cell to start Newton-Raphson iterations
4660 xi0 = 0.5
467
468 ! Newton-Raphson iterations
4690 do iter = 1,NR_ITERATIONS
470
4710 if (use_2018_answers) then
472 numerator = a(1) + a(2)*xi0 + a(3)*xi0*xi0 + a(4)*xi0*xi0*xi0 + &
4730 a(5)*xi0*xi0*xi0*xi0 - target_value
4740 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.
4760 numerator = (a(1) - target_value) + xi0*(a(2) + xi0*(a(3) + xi0*(a(4) + a(5)*xi0)))
4770 denominator = a(2) + xi0*(2.*a(3) + xi0*(3.*a(4) + 4.*a(5)*xi0))
478 endif
479
4800 delta = -numerator / denominator
481
4820 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).
4890 if ( xi0 < 0.0 ) then
4900 xi0 = 0.0
4910 grad = a(2)
4920 if ( grad == 0.0 ) xi0 = xi0 + eps
493 endif
494
4950 if ( xi0 > 1.0 ) then
4960 xi0 = 1.0
4970 if (use_2018_answers) then
4980 grad = a(2) + 2*a(3) + 3*a(4) + 4*a(5)
499 else ! These expressions are mathematicaly equivalent but more accurate.
5000 grad = a(2) + (2.*a(3) + (3.*a(4) + 4.*a(5)))
501 endif
5020 if ( grad == 0.0 ) xi0 = xi0 - eps
503 endif
504
505 ! break if converged or too many iterations taken
5060 if ( abs(delta) < NR_TOLERANCE ) exit
507 enddo ! end Newton-Raphson iterations
508
5090 x_tgt = x_g(k_found) + xi0 * h(k_found)
5100end function get_polynomial_coordinate
511
512!> Numeric value of interpolation_scheme corresponding to scheme name
5130integer 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
5180 select case ( uppercase(trim(interp_scheme)) )
5190 case ("P1M_H2"); interpolation_scheme = INTERPOLATION_P1M_H2
5200 case ("P1M_H4"); interpolation_scheme = INTERPOLATION_P1M_H4
5210 case ("P1M_IH2"); interpolation_scheme = INTERPOLATION_P1M_IH4
5220 case ("PLM"); interpolation_scheme = INTERPOLATION_PLM
5230 case ("PPM_CW"); interpolation_scheme = INTERPOLATION_PPM_CW
5240 case ("PPM_H4"); interpolation_scheme = INTERPOLATION_PPM_H4
5250 case ("PPM_IH4"); interpolation_scheme = INTERPOLATION_PPM_IH4
5260 case ("P3M_IH4IH3"); interpolation_scheme = INTERPOLATION_P3M_IH4IH3
5270 case ("P3M_IH6IH5"); interpolation_scheme = INTERPOLATION_P3M_IH6IH5
5280 case ("PQM_IH4IH3"); interpolation_scheme = INTERPOLATION_PQM_IH4IH3
5290 case ("PQM_IH6IH5"); interpolation_scheme = INTERPOLATION_PQM_IH6IH5
530 case default ; call MOM_error(FATAL, "regrid_interp: "//&
5310 "Unrecognized choice for INTERPOLATION_SCHEME ("//trim(interp_scheme)//").")
532 end select
5330end function interpolation_scheme
534
535!> Store the interpolation_scheme value in the interp_CS based on the input string.
5360subroutine 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
5420 CS%interpolation_scheme = interpolation_scheme(interp_scheme)
5430end subroutine set_interp_scheme
544
545!> Store the boundary_extrapolation value in the interp_CS
5462subroutine 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
5512 CS%boundary_extrapolation = extrap
5522end subroutine set_interp_extrap
553
554!> Store the value of the answer_date in the interp_CS
5550subroutine 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
5600 CS%answer_date = answer_date
5610end subroutine set_interp_answer_date
562
5630end module regrid_interp