13 use plm_functions,
only : plm_reconstruction, plm_boundary_extrapolation
14 use ppm_functions,
only : ppm_reconstruction, ppm_boundary_extrapolation
15 use pqm_functions,
only : pqm_reconstruction, pqm_boundary_extrapolation_v1
17 use p1m_functions,
only : p1m_interpolation, p1m_boundary_extrapolation
18 use p3m_functions,
only : p3m_interpolation, p3m_boundary_extrapolation
20 implicit none ;
private
28 integer :: interpolation_scheme = -1
32 logical :: boundary_extrapolation
35 public regridding_set_ppolys, interpolate_grid
36 public build_and_interpolate_grid
37 public set_interp_scheme, set_interp_extrap
40 integer,
parameter :: interpolation_p1m_h2 = 0
41 integer,
parameter :: interpolation_p1m_h4 = 1
42 integer,
parameter :: interpolation_p1m_ih4 = 2
43 integer,
parameter :: interpolation_plm = 3
44 integer,
parameter :: interpolation_ppm_h4 = 4
45 integer,
parameter :: interpolation_ppm_ih4 = 5
46 integer,
parameter :: interpolation_p3m_ih4ih3 = 6
47 integer,
parameter :: interpolation_p3m_ih6ih5 = 7
48 integer,
parameter :: interpolation_pqm_ih4ih3 = 8
49 integer,
parameter :: interpolation_pqm_ih6ih5 = 9
52 integer,
parameter :: degree_1 = 1, degree_2 = 2, degree_3 = 3, degree_4 = 4
53 integer,
public,
parameter :: degree_max = 5
59 real,
public,
parameter :: nr_offset = 1e-6
63 integer,
public,
parameter :: nr_iterations = 8
65 real,
public,
parameter :: nr_tolerance = 1e-12
75 subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, &
76 ppoly0_coefs, degree, h_neglect, h_neglect_edge)
78 real,
dimension(:),
intent(in) :: densities
79 integer,
intent(in) :: n0
80 real,
dimension(:),
intent(in) :: h0
81 real,
dimension(:,:),
intent(inout) :: ppoly0_e
82 real,
dimension(:,:),
intent(inout) :: ppoly0_s
83 real,
dimension(:,:),
intent(inout) :: ppoly0_coefs
84 integer,
intent(inout) :: degree
85 real,
optional,
intent(in) :: h_neglect
88 real,
optional,
intent(in) :: h_neglect_edge
92 logical :: extrapolate
97 ppoly0_coefs(:,:) = 0.0
99 extrapolate = cs%boundary_extrapolation
102 select case (cs%interpolation_scheme)
104 case ( interpolation_p1m_h2 )
106 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
107 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
108 if (extrapolate)
then
109 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
112 case ( interpolation_p1m_h4 )
115 call edge_values_explicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge )
117 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
119 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
120 if (extrapolate)
then
121 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
124 case ( interpolation_p1m_ih4 )
127 call edge_values_implicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge )
129 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
131 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
132 if (extrapolate)
then
133 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
136 case ( interpolation_plm )
138 call plm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
139 if (extrapolate)
then
140 call plm_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
143 case ( interpolation_ppm_h4 )
146 call edge_values_explicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge )
147 call ppm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
148 if (extrapolate)
then
149 call ppm_boundary_extrapolation( n0, h0, densities, ppoly0_e, &
150 ppoly0_coefs, h_neglect )
154 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
155 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
156 if (extrapolate)
then
157 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
161 case ( interpolation_ppm_ih4 )
164 call edge_values_implicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge )
165 call ppm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
166 if (extrapolate)
then
167 call ppm_boundary_extrapolation( n0, h0, densities, ppoly0_e, &
168 ppoly0_coefs, h_neglect )
172 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
173 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
174 if (extrapolate)
then
175 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
179 case ( interpolation_p3m_ih4ih3 )
182 call edge_values_implicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge )
183 call edge_slopes_implicit_h3( n0, h0, densities, ppoly0_s, h_neglect )
184 call p3m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_s, &
185 ppoly0_coefs, h_neglect )
186 if (extrapolate)
then
187 call p3m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_s, &
188 ppoly0_coefs, h_neglect, h_neglect_edge )
192 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
193 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
194 if (extrapolate)
then
195 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
199 case ( interpolation_p3m_ih6ih5 )
202 call edge_values_implicit_h6( n0, h0, densities, ppoly0_e, h_neglect_edge )
203 call edge_slopes_implicit_h5( n0, h0, densities, ppoly0_s, h_neglect )
204 call p3m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_s, &
205 ppoly0_coefs, h_neglect )
206 if (extrapolate)
then
207 call p3m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_s, &
208 ppoly0_coefs, h_neglect, h_neglect_edge )
212 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
213 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
214 if (extrapolate)
then
215 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
219 case ( interpolation_pqm_ih4ih3 )
222 call edge_values_implicit_h4( n0, h0, densities, ppoly0_e, h_neglect_edge )
223 call edge_slopes_implicit_h3( n0, h0, densities, ppoly0_s, h_neglect )
224 call pqm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_s, &
225 ppoly0_coefs, h_neglect )
226 if (extrapolate)
then
227 call pqm_boundary_extrapolation_v1( n0, h0, densities, ppoly0_e, ppoly0_s, &
228 ppoly0_coefs, h_neglect )
232 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
233 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
234 if (extrapolate)
then
235 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
239 case ( interpolation_pqm_ih6ih5 )
242 call edge_values_implicit_h6( n0, h0, densities, ppoly0_e, h_neglect_edge )
243 call edge_slopes_implicit_h5( n0, h0, densities, ppoly0_s, h_neglect )
244 call pqm_reconstruction( n0, h0, densities, ppoly0_e, ppoly0_s, &
245 ppoly0_coefs, h_neglect )
246 if (extrapolate)
then
247 call pqm_boundary_extrapolation_v1( n0, h0, densities, ppoly0_e, ppoly0_s, &
248 ppoly0_coefs, h_neglect )
252 call edge_values_explicit_h2( n0, h0, densities, ppoly0_e, h_neglect_edge )
253 call p1m_interpolation( n0, h0, densities, ppoly0_e, ppoly0_coefs, h_neglect )
254 if (extrapolate)
then
255 call p1m_boundary_extrapolation( n0, h0, densities, ppoly0_e, ppoly0_coefs )
259 end subroutine regridding_set_ppolys
266 subroutine interpolate_grid( n0, h0, x0, ppoly0_E, ppoly0_coefs, &
267 target_values, degree, n1, h1, x1 )
268 integer,
intent(in) :: n0
269 real,
dimension(:),
intent(in) :: h0
270 real,
dimension(:),
intent(in) :: x0
271 real,
dimension(:,:),
intent(in) :: ppoly0_e
272 real,
dimension(:,:),
intent(in) :: ppoly0_coefs
273 real,
dimension(:),
intent(in) :: target_values
274 integer,
intent(in) :: degree
275 integer,
intent(in) :: n1
276 real,
dimension(:),
intent(inout) :: h1
277 real,
dimension(:),
intent(inout) :: x1
290 x1(k) = get_polynomial_coordinate( n0, h0, x0, ppoly0_e, ppoly0_coefs, t, degree )
291 h1(k-1) = x1(k) - x1(k-1)
293 h1(n1) = x1(n1+1) - x1(n1)
295 end subroutine interpolate_grid
298 subroutine build_and_interpolate_grid(CS, densities, n0, h0, x0, target_values, &
299 n1, h1, x1, h_neglect, h_neglect_edge)
301 real,
dimension(:),
intent(in) :: densities
302 real,
dimension(:),
intent(in) :: target_values
303 integer,
intent(in) :: n0
304 real,
dimension(:),
intent(in) :: h0
305 real,
dimension(:),
intent(in) :: x0
306 integer,
intent(in) :: n1
307 real,
dimension(:),
intent(inout) :: h1
308 real,
dimension(:),
intent(inout) :: x1
309 real,
optional,
intent(in) :: h_neglect
312 real,
optional,
intent(in) :: h_neglect_edge
316 real,
dimension(n0,2) :: ppoly0_e, ppoly0_s
317 real,
dimension(n0,DEGREE_MAX+1) :: ppoly0_c
320 call regridding_set_ppolys(cs, densities, n0, h0, ppoly0_e, ppoly0_s, ppoly0_c, &
321 degree, h_neglect, h_neglect_edge)
322 call interpolate_grid(n0, h0, x0, ppoly0_e, ppoly0_c, target_values, degree, &
324 end subroutine build_and_interpolate_grid
342 function get_polynomial_coordinate( N, h, x_g, ppoly_E, ppoly_coefs, &
343 target_value, degree )
result ( x_tgt )
345 integer,
intent(in) :: n
346 real,
dimension(:),
intent(in) :: h
347 real,
dimension(:),
intent(in) :: x_g
348 real,
dimension(:,:),
intent(in) :: ppoly_e
349 real,
dimension(:,:),
intent(in) :: ppoly_coefs
350 real,
intent(in) :: target_value
351 integer,
intent(in) :: degree
358 real,
dimension(DEGREE_MAX) :: a
373 if ( target_value <= ppoly_e(1,1) )
then
381 if ( ( target_value >= ppoly_e(k-1,2) ) .AND. &
382 ( target_value <= ppoly_e(k,1) ) )
then
392 if ( target_value >= ppoly_e(n,2) )
then
403 if ( ( target_value > ppoly_e(k,1) ) .AND. &
404 ( target_value < ppoly_e(k,2) ) )
then
415 if ( k_found == -1 )
then
416 write(*,*) target_value, ppoly_e(1,1), ppoly_e(n,2)
417 write(*,*)
'Could not find target coordinate in ' //&
418 '"get_polynomial_coordinate". This is caused by an '//&
419 'inconsistent interpolant (perhaps not monotonically '//&
421 call mom_error( fatal,
'Aborting execution' )
428 a(i) = ppoly_coefs(k_found,i)
439 if ( ( iter > nr_iterations ) .OR. &
440 ( abs(delta) < nr_tolerance ) )
then
444 numerator = a(1) + a(2)*xi0 + a(3)*xi0*xi0 + a(4)*xi0*xi0*xi0 + &
445 a(5)*xi0*xi0*xi0*xi0 - target_value
447 denominator = a(2) + 2*a(3)*xi0 + 3*a(4)*xi0*xi0 + 4*a(5)*xi0*xi0*xi0
449 delta = -numerator / denominator
458 if ( xi0 < 0.0 )
then
461 if ( grad == 0.0 ) xi0 = xi0 + eps
464 if ( xi0 > 1.0 )
then
466 grad = a(2) + 2*a(3) + 3*a(4) + 4*a(5)
467 if ( grad == 0.0 ) xi0 = xi0 - eps
473 x_tgt = x_g(k_found) + xi0 * h(k_found)
474 end function get_polynomial_coordinate
477 integer function interpolation_scheme(interp_scheme)
478 character(len=*),
intent(in) :: interp_scheme
482 select case ( uppercase(trim(interp_scheme)) )
483 case (
"P1M_H2"); interpolation_scheme = interpolation_p1m_h2
484 case (
"P1M_H4"); interpolation_scheme = interpolation_p1m_h4
485 case (
"P1M_IH2"); interpolation_scheme = interpolation_p1m_ih4
486 case (
"PLM"); interpolation_scheme = interpolation_plm
487 case (
"PPM_H4"); interpolation_scheme = interpolation_ppm_h4
488 case (
"PPM_IH4"); interpolation_scheme = interpolation_ppm_ih4
489 case (
"P3M_IH4IH3"); interpolation_scheme = interpolation_p3m_ih4ih3
490 case (
"P3M_IH6IH5"); interpolation_scheme = interpolation_p3m_ih6ih5
491 case (
"PQM_IH4IH3"); interpolation_scheme = interpolation_pqm_ih4ih3
492 case (
"PQM_IH6IH5"); interpolation_scheme = interpolation_pqm_ih6ih5
493 case default ;
call mom_error(fatal,
"regrid_interp: "//&
494 "Unrecognized choice for INTERPOLATION_SCHEME ("//trim(interp_scheme)//
").")
496 end function interpolation_scheme
499 subroutine set_interp_scheme(CS, interp_scheme)
501 character(len=*),
intent(in) :: interp_scheme
505 cs%interpolation_scheme = interpolation_scheme(interp_scheme)
506 end subroutine set_interp_scheme
509 subroutine set_interp_extrap(CS, extrap)
511 logical,
intent(in) :: extrap
514 cs%boundary_extrapolation = extrap
515 end subroutine set_interp_extrap