Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module interp_1d_lib_sg
21 : implicit none
22 :
23 : contains
24 :
25 :
26 : ! this routine is a simply wrapper for making an interpolant and then using it.
27 0 : subroutine interpolate_vector_sg( &
28 0 : n_old, x_old, n_new, x_new, v_old, v_new, interp_vec_sg, nwork, work1, str, ierr)
29 : integer, intent(in) :: n_old, n_new
30 : real, intent(in) :: x_old(:) ! (n_old)
31 : real, intent(in) :: v_old(:) ! (n_old)
32 : real, intent(in) :: x_new(:) ! (n_new)
33 : real, intent(inout) :: v_new(:) ! (n_new)
34 : interface
35 : subroutine interp_vec_sg(x, nx, f1, nwork, work1, str, ierr) ! make cubic interpolant
36 : ! e.g., interp_pm, interp_m3a, interp_m3b, or interp_m3q
37 : implicit none
38 : integer, intent(in) :: nx ! length of x vector
39 : real, intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
40 : real, intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
41 : integer, intent(in) :: nwork
42 : real, intent(inout), pointer :: work1(:)
43 : character (len=*) :: str ! for debugging
44 : integer, intent(out) :: ierr
45 : end subroutine interp_vec_sg
46 : end interface
47 : integer, intent(in) :: nwork
48 : real, intent(inout), pointer :: work1(:) ! =(:,:) ! (n_old, nwork)
49 : character (len=*) :: str ! for debugging
50 : integer, intent(out) :: ierr
51 0 : real, pointer :: f1(:), f(:,:)
52 : integer :: k
53 0 : ierr = 0
54 0 : allocate(f1(4*n_old), stat=ierr)
55 0 : if (ierr /= 0) return
56 0 : f(1:4,1:n_old) => f1(1:4*n_old)
57 0 : forall (k=1:n_old) f(1,k) = v_old(k)
58 0 : call interp_vec_sg(x_old, n_old, f1, nwork, work1, str, ierr) ! make interpolant
59 0 : if (ierr /= 0) then
60 0 : deallocate(f1)
61 0 : return
62 : end if
63 0 : call interp_values_sg(x_old, n_old, f1, n_new, x_new, v_new, ierr)
64 0 : deallocate(f1)
65 0 : end subroutine interpolate_vector_sg
66 :
67 :
68 : ! this routine is a simply wrapper for making an interpolant with interp_pm and then using it.
69 1 : subroutine interpolate_vector_pm_sg( &
70 1 : n_old, x_old, n_new, x_new, v_old, v_new, work1, str, ierr)
71 : use interp_1d_def, only: pm_work_size
72 : integer, intent(in) :: n_old, n_new
73 : real, intent(in) :: x_old(:) ! (n_old)
74 : real, intent(in) :: v_old(:) ! (n_old)
75 : real, intent(in) :: x_new(:) ! (n_new)
76 : real, intent(inout) :: v_new(:) ! (n_new)
77 : real, intent(inout), pointer :: work1(:) ! =(:,:) ! (n_old, nwork)
78 : character (len=*) :: str ! for debugging
79 : integer, intent(out) :: ierr
80 1 : real, pointer :: f1(:), f(:,:)
81 : integer :: k
82 1 : ierr = 0
83 1 : allocate(f1(4*n_old), stat=ierr)
84 1 : if (ierr /= 0) return
85 1 : f(1:4,1:n_old) => f1(1:4*n_old)
86 4 : forall (k=1:n_old) f(1,k) = v_old(k)
87 1 : call interp_pm_sg(x_old, n_old, f1, pm_work_size, work1, str, ierr) ! make interpolant
88 1 : if (ierr /= 0) then
89 0 : deallocate(f1)
90 0 : return
91 : end if
92 1 : call interp_values_sg(x_old, n_old, f1, n_new, x_new, v_new, ierr)
93 1 : deallocate(f1)
94 1 : end subroutine interpolate_vector_pm_sg
95 :
96 :
97 0 : subroutine interp_4_to_1_sg( &
98 : pdqm1, pdq00, pdqp1, ndq00, pfm1, pf00, pfp1, pfp2, nf00, str, ierr)
99 : ! 4 points in, 1 point out
100 : ! piecewise monotonic cubic interpolation
101 : use interp_1d_def, only: pm_work_size
102 : real, intent(in) :: pdqm1, pdq00, pdqp1 ! spacing between input points
103 : real, intent(in) :: ndq00
104 : real, intent(in) :: pfm1, pf00, pfp1, pfp2 ! values at input points
105 : real, intent(out) :: nf00 ! new value at ndq00
106 : character (len=*) :: str ! for debugging
107 : integer, intent(out) :: ierr
108 : integer, parameter :: n_old=4, n_new=1
109 0 : real :: x_old(n_old), v_old(n_old), x_new(n_new), v_new(n_new)
110 0 : real, target :: work1_ary(n_old*pm_work_size)
111 : real, pointer :: work1(:)
112 0 : work1 => work1_ary
113 : ierr = 0
114 0 : x_old(1) = 0.0
115 0 : x_old(2) = pdqm1
116 0 : x_old(3) = pdqm1+pdq00
117 0 : x_old(4) = pdqm1+pdq00+pdqp1
118 0 : v_old(1) = pfm1
119 0 : v_old(2) = pf00
120 0 : v_old(3) = pfp1
121 0 : v_old(4) = pfp2
122 0 : x_new(1) = ndq00
123 : call interpolate_vector_pm_sg( &
124 0 : n_old, x_old, n_new, x_new, v_old, v_new, work1, str, ierr)
125 0 : nf00 = v_new(1)
126 0 : end subroutine interp_4_to_1_sg
127 :
128 :
129 0 : subroutine interp_3_to_1_sg(pdqm1, pdq00, ndqm1, pfm1, pf00, pfp1, nf00, str, ierr)
130 : ! 3 points in, 1 point out
131 : ! piecewise monotonic quadratic interpolation
132 : use interp_1d_def, only: pm_work_size
133 : real, intent(in) :: pdqm1, pdq00 ! spacing between input points
134 : real, intent(in) :: ndqm1 ! new spacing to nk
135 : real, intent(in) :: pfm1, pf00, pfp1 ! values at pkm1, pk, pkp1
136 : real, intent(out) :: nf00 ! new value at ndq00
137 : character (len=*) :: str ! for debugging
138 : integer, intent(out) :: ierr
139 : integer, parameter :: n_old=3, n_new=1
140 0 : real :: x_old(n_old), v_old(n_old), x_new(n_new), v_new(n_new)
141 0 : real, target :: work1_ary(n_old*pm_work_size)
142 : real, pointer :: work1(:)
143 0 : work1 => work1_ary
144 : ierr = 0
145 0 : x_old(1) = 0.0
146 0 : x_old(2) = pdqm1
147 0 : x_old(3) = pdqm1+pdq00
148 0 : v_old(1) = pfm1
149 0 : v_old(2) = pf00
150 0 : v_old(3) = pfp1
151 0 : x_new(1) = ndqm1
152 : call interpolate_vector_pm_sg( &
153 0 : n_old, x_old, n_new, x_new, v_old, v_new, work1, str, ierr)
154 0 : nf00 = v_new(1)
155 0 : end subroutine interp_3_to_1_sg
156 :
157 :
158 1 : subroutine interp_3_to_2_sg(pdqm1, pdq00, ndqm1, ndq00, pfm1, pf00, pfp1, nf00, nfp1, str, ierr)
159 : ! 3 points in, 2 points out
160 : ! piecewise monotonic quadratic interpolation
161 : use interp_1d_def, only: pm_work_size
162 : real, intent(in) :: pdqm1, pdq00 ! previous spacing to pk and pkp1
163 : real, intent(in) :: ndqm1, ndq00 ! new spacing to nk and nk+1
164 : real, intent(in) :: pfm1, pf00, pfp1 ! previous values at pkm1, pk, pkp1
165 : real, intent(out) :: nf00, nfp1 ! new values at nk, nk+1
166 : character (len=*) :: str ! for debugging
167 : integer, intent(out) :: ierr
168 : integer, parameter :: n_old=3, n_new=2
169 14 : real :: x_old(n_old), v_old(n_old), x_new(n_new), v_new(n_new)
170 10 : real, target :: work1_ary(n_old*pm_work_size)
171 : real, pointer :: work1(:)
172 1 : work1 => work1_ary
173 : ierr = 0
174 1 : x_old(1) = 0d0
175 1 : x_old(2) = pdqm1
176 1 : x_old(3) = pdqm1+pdq00
177 1 : v_old(1) = pfm1
178 1 : v_old(2) = pf00
179 1 : v_old(3) = pfp1
180 1 : x_new(1) = ndqm1
181 1 : x_new(2) = ndqm1+ndq00
182 : call interpolate_vector_pm_sg( &
183 1 : n_old, x_old, n_new, x_new, v_old, v_new, work1, str, ierr)
184 1 : nf00 = v_new(1)
185 1 : nfp1 = v_new(2)
186 1 : end subroutine interp_3_to_2_sg
187 :
188 :
189 : ! general routines
190 :
191 : ! these routines use previously created interpolant information (f)
192 : ! the interpolant can come from either the piecewise monotonic routines, or
193 : ! from the monotonicity preserving routines -- they use the same format for f.
194 :
195 15 : subroutine interp_values_sg(init_x, nx, f1, nv, x, vals, ierr)
196 : use interp_1d_def
197 : use interp_1d_misc_sg
198 : real, intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
199 : integer, intent(in) :: nx ! length of init_x vector
200 : real, intent(in), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
201 : integer, intent(in) :: nv ! length of new x vector and vals vector
202 : real, intent(in) :: x(:) ! (nv) ! locations where want interpolated values
203 : ! strictly monotonic in same way as init_x
204 : ! values out of range of init_x's are clipped to boundaries of init_x's
205 : real, intent(inout) :: vals(:) ! (nv)
206 : integer, intent(out) :: ierr ! 0 means AOK
207 15 : call do_interp_values_sg(init_x, nx, f1, nv, x, vals, ierr)
208 15 : end subroutine interp_values_sg
209 :
210 :
211 0 : subroutine interp_value_sg(init_x, nx, f1, xval, val, ierr)
212 : use interp_1d_def
213 : use interp_1d_misc_sg
214 : real, intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
215 : integer, intent(in) :: nx ! length of init_x vector
216 : real, intent(in), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
217 : real, intent(in) :: xval ! location where want interpolated values
218 : real, intent(out) :: val
219 : integer, intent(out) :: ierr ! 0 means AOK
220 : integer, parameter :: nv = 1
221 0 : real :: x(nv), vals(nv)
222 0 : x(1) = xval
223 0 : call do_interp_values_sg(init_x, nx, f1, nv, x, vals, ierr)
224 0 : val = vals(1)
225 0 : end subroutine interp_value_sg
226 :
227 :
228 0 : subroutine interp_values_and_slopes_sg(init_x, nx, f1, nv, x, vals, slopes, ierr)
229 : use interp_1d_def
230 : use interp_1d_misc_sg
231 : real, intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
232 : integer, intent(in) :: nx ! length of init_x vector
233 : real, intent(in), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
234 : integer, intent(in) :: nv ! length of new x vector and vals and slopes vectors
235 : real, intent(in) :: x(:) ! (nv) ! locations where want interpolated values
236 : ! strictly monotonic in same way as init_x
237 : ! values out of range of init_x's are clipped to boundaries of init_x's
238 : real, intent(inout) :: vals(:) ! (nv)
239 : real, intent(inout) :: slopes(:) ! (nv)
240 : integer, intent(out) :: ierr ! 0 means AOK
241 0 : call do_interp_values_and_slopes_sg(init_x, nx, f1, nv, x, vals, slopes, ierr)
242 0 : end subroutine interp_values_and_slopes_sg
243 :
244 :
245 0 : subroutine interp_value_and_slope_sg(init_x, nx, f1, xval, val, slope, ierr)
246 : use interp_1d_def
247 : use interp_1d_misc_sg
248 : real, intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
249 : integer, intent(in) :: nx ! length of init_x vector
250 : real, intent(in), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
251 : real, intent(in) :: xval ! location where want interpolated values
252 : real, intent(out) :: val, slope
253 : integer, intent(out) :: ierr ! 0 means AOK
254 : integer, parameter :: nv = 1
255 0 : real :: x(nv), vals(nv), slopes(nv)
256 0 : x(1) = xval
257 0 : call do_interp_values_and_slopes_sg(init_x, nx, f1, nv, x, vals, slopes, ierr)
258 0 : val = vals(1)
259 0 : slope = slopes(1)
260 0 : end subroutine interp_value_and_slope_sg
261 :
262 :
263 2 : subroutine integrate_values_sg(init_x, nx, f1, nv, x, vals, ierr)
264 : use interp_1d_def
265 : use interp_1d_misc_sg
266 : real, intent(in) :: init_x(:) ! (nx) ! junction points, strictly increasing
267 : integer, intent(in) :: nx ! length of init_x vector
268 : real, intent(in), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
269 : integer, intent(in) :: nv ! length of new x vector and vals vector
270 : real, intent(in) :: x(:) ! (nv)
271 : ! strictly monotonic in same way as init_x
272 : ! NOTE: no extrapolation allowed -- x's must be within range of init_x's
273 : real, intent(inout) :: vals(:) ! (nv)
274 : ! for i > 1, vals(i) = integral of interpolating poly from x(i-1) to x(i)
275 : ! vals(1) = 0
276 : integer, intent(out) :: ierr ! 0 means AOK
277 :
278 2 : call do_integrate_values_sg(init_x, nx, f1, nv, x, vals, ierr)
279 :
280 2 : end subroutine integrate_values_sg
281 :
282 :
283 : ! piecewise monotonic routines
284 :
285 : ! the following produce piecewise monotonic interpolants rather than monotonicity preserving
286 : ! this stricter limit never introduces interpolated values exceeding the given values,
287 : ! even in places where the given values are not monotonic.
288 : ! the downside is reduced accuracy on smooth data compared to the mp routines.
289 :
290 :
291 : ! Steffen, M., "A simple method for monotonic interpolation in one dimension",
292 : ! Astron. Astrophys., (239) 1990, 443-450.
293 :
294 :
295 46 : subroutine interp_pm_sg(x, nx, f1, nwork, work1, str, ierr)
296 : ! make piecewise monotonic cubic interpolant
297 : use interp_1d_def
298 : use interp_1d_pm_sg
299 : integer, intent(in) :: nx ! length of x vector (>= 2)
300 : real, intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
301 : real, intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
302 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
303 : real, intent(inout), pointer :: work1(:) ! =(nx,nwork)
304 : character (len=*) :: str ! for debugging
305 : integer, intent(out) :: ierr
306 46 : call mk_pmcub_sg(x, nx, f1, .false., nwork, work1, str, ierr)
307 46 : end subroutine interp_pm_sg
308 :
309 :
310 0 : subroutine interp_pm_slopes_only_sg(x, nx, f1, nwork, work1, str, ierr)
311 : ! identical to interp_pm, but only calculates slopes and stores them in f(2,:)
312 : ! this is a little faster for the special case in which you just want the slopes at x
313 : use interp_1d_def
314 : use interp_1d_pm_sg
315 : integer, intent(in) :: nx ! length of x vector (>= 2)
316 : real, intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
317 : real, intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
318 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
319 : real, intent(inout), pointer :: work1(:) ! =(nx,nwork)
320 : character (len=*) :: str ! for debugging
321 : integer, intent(out) :: ierr
322 0 : call mk_pmcub_sg(x, nx, f1, .true., nwork, work1, str, ierr)
323 0 : end subroutine interp_pm_slopes_only_sg
324 :
325 :
326 1 : subroutine interp_4pt_pm_sg(x, y, a)
327 : ! returns coefficients for monotonic cubic interpolation from x(2) to x(3)
328 : real, intent(in) :: x(4) ! junction points, strictly monotonic
329 : real, intent(in) :: y(4) ! data values at x's
330 : real, intent(inout) :: a(3) ! coefficients
331 1 : real :: h1, h2, h3, s1, s2, s3, p2, p3, as2, ss2, yp2, yp3
332 : ! for x(2) <= x <= x(3) and dx = x-x(2),
333 : ! y(x) = y(2) + dx*(a(1) + dx*(a(2) + dx*a(3)))
334 1 : h1 = x(2)-x(1)
335 1 : h2 = x(3)-x(2)
336 1 : h3 = x(4)-x(3)
337 1 : s1 = (y(2)-y(1))/h1
338 1 : s2 = (y(3)-y(2))/h2
339 1 : s3 = (y(4)-y(3))/h3
340 1 : p2 = (s1*h2+s2*h1)/(h1+h2)
341 1 : p3 = (s2*h3+s3*h2)/(h2+h3)
342 1 : as2 = abs(s2)
343 1 : ss2 = sign(1.0, s2)
344 1 : yp2 = (sign(1.0, s1)+ss2)*min(abs(s1), as2, 0.5*abs(p2))
345 1 : yp3 = (ss2+sign(1.0, s3))*min(as2, abs(s3), 0.5*abs(p3))
346 1 : a(1) = yp2
347 1 : a(2) = (3*s2-2*yp2-yp3)/h2
348 1 : a(3) = (yp2+yp3-2*s2)/(h2*h2)
349 1 : end subroutine interp_4pt_pm_sg
350 :
351 :
352 0 : subroutine interp_pm_on_uniform_grid_sg(dx, nx, f1, nwork, work1, str, ierr)
353 : ! make piecewise monotonic cubic interpolant on uniformly spaced mesh
354 : use interp_1d_def
355 : use interp_1d_pm_sg
356 : real, intent(in) :: dx ! grid spacing
357 : integer, intent(in) :: nx ! length of vector (>= 2)
358 : real, intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
359 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
360 : real, intent(inout), pointer :: work1(:) ! =(nx, nwork)
361 : character (len=*) :: str ! for debugging
362 : integer, intent(out) :: ierr
363 0 : call mk_pmcub_uniform_sg(dx, nx, f1, .false., nwork, work1, str, ierr)
364 0 : end subroutine interp_pm_on_uniform_grid_sg
365 :
366 :
367 :
368 : ! monotonicity preserving routines
369 :
370 : ! Huynh, H.T., "Accurate Monotone Cubic Interpolation", SIAM J Numer. Anal. (30) 1993, 57-100.
371 :
372 : ! Suresh, A, and H.T. Huynh, "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
373 : ! Time Stepping", JCP (136) 1997, 83-99.
374 :
375 :
376 0 : subroutine interp_m3_sg(x, nx, f1, which, nwork, work1, str, ierr)
377 : ! make monotonicity preserving cubic interpolant on arbitrarily spaced grid
378 : use interp_1d_def
379 : use interp_1d_mp_sg
380 : integer, intent(in) :: nx ! length of x vector (>= 4)
381 : real, intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
382 : real, intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
383 : integer, intent(in) :: which ! average, quartic, or super_bee
384 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
385 : real, intent(inout), pointer :: work1(:) ! =(nx, nwork)
386 : character (len=*) :: str ! for debugging
387 : integer, intent(out) :: ierr
388 0 : call m3_sg(x, nx, f1, which, .false., nwork, work1, str, ierr)
389 0 : end subroutine interp_m3_sg
390 :
391 :
392 0 : subroutine interp_m3a_sg(x, nx, f1, nwork, work1, str, ierr)
393 : ! make monotonicity preserving cubic interpolant on arbitrarily spaced grid
394 : use interp_1d_def
395 : use interp_1d_mp_sg
396 : integer, intent(in) :: nx ! length of x vector (>= 4)
397 : real, intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
398 : real, intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
399 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
400 : real, intent(inout), pointer :: work1(:) ! =(nx, nwork)
401 : character (len=*) :: str ! for debugging
402 : integer, intent(out) :: ierr
403 0 : call m3_sg(x, nx, f1, average, .false., nwork, work1, str, ierr)
404 0 : end subroutine interp_m3a_sg
405 :
406 :
407 0 : subroutine interp_m3b_sg(x, nx, f1, nwork, work1, str, ierr)
408 : ! make monotonicity preserving cubic interpolant on arbitrarily spaced grid
409 : use interp_1d_def
410 : use interp_1d_mp_sg
411 : integer, intent(in) :: nx ! length of x vector (>= 4)
412 : real, intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
413 : real, intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
414 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
415 : real, intent(inout), pointer :: work1(:) ! =(nx, nwork)
416 : character (len=*) :: str ! for debugging
417 : integer, intent(out) :: ierr
418 0 : call m3_sg(x, nx, f1, super_bee, .false., nwork, work1, str, ierr)
419 0 : end subroutine interp_m3b_sg
420 :
421 :
422 4 : subroutine interp_m3q_sg(x, nx, f1, nwork, work1, str, ierr)
423 : ! make monotonicity preserving cubic interpolant on arbitrarily spaced grid
424 : use interp_1d_def
425 : use interp_1d_mp_sg
426 : integer, intent(in) :: nx ! length of x vector (>= 4)
427 : real, intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
428 : real, intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
429 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
430 : real, intent(inout), pointer :: work1(:) ! =(nx, nwork)
431 : character (len=*) :: str ! for debugging
432 : integer, intent(out) :: ierr
433 4 : call m3_sg(x, nx, f1, quartic, .false., nwork, work1, str, ierr)
434 4 : end subroutine interp_m3q_sg
435 :
436 :
437 0 : subroutine interp_m3_on_uniform_grid_sg(dx, nx, f1, which, nwork, work1, str, ierr)
438 : ! make monotonicity preserving cubic interpolant on uniformly spaced grid
439 : use interp_1d_def
440 : use interp_1d_mp_sg
441 : real, intent(in) :: dx ! the grid spacing
442 : integer, intent(in) :: nx ! length of x vector (>= 4)
443 : real, intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
444 : integer, intent(in) :: which ! average, quartic, or super_bee
445 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
446 : real, intent(inout), pointer :: work1(:) ! =(nx, nwork)
447 : character (len=*) :: str ! for debugging
448 : integer, intent(out) :: ierr
449 0 : call m3_on_uniform_grid_sg(dx, nx, f1, which, .false., nwork, work1, str, ierr)
450 0 : end subroutine interp_m3_on_uniform_grid_sg
451 :
452 :
453 0 : subroutine interp_m3a_on_uniform_grid_sg(dx, nx, f1, nwork, work1, str, ierr)
454 : ! make monotonicity preserving cubic interpolant on uniformly spaced grid
455 : use interp_1d_def
456 : use interp_1d_mp_sg
457 : real, intent(in) :: dx ! the grid spacing
458 : integer, intent(in) :: nx ! length of x vector (>= 4)
459 : real, intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
460 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
461 : real, intent(inout), pointer :: work1(:) ! =(nx, nwork)
462 : character (len=*) :: str ! for debugging
463 : integer, intent(out) :: ierr
464 0 : call m3_on_uniform_grid_sg(dx, nx, f1, average, .false., nwork, work1, str, ierr)
465 0 : end subroutine interp_m3a_on_uniform_grid_sg
466 :
467 :
468 0 : subroutine interp_m3b_on_uniform_grid_sg(dx, nx, f1, nwork, work1, str, ierr)
469 : ! make monotonicity preserving cubic interpolant on uniformly spaced grid
470 : use interp_1d_def
471 : use interp_1d_mp_sg
472 : real, intent(in) :: dx ! the grid spacing
473 : integer, intent(in) :: nx ! length of x vector (>= 4)
474 : real, intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
475 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
476 : real, intent(inout), pointer :: work1(:) ! =(nx, nwork)
477 : character (len=*) :: str ! for debugging
478 : integer, intent(out) :: ierr
479 0 : call m3_on_uniform_grid_sg(dx, nx, f1, super_bee, .false., nwork, work1, str, ierr)
480 0 : end subroutine interp_m3b_on_uniform_grid_sg
481 :
482 :
483 0 : subroutine interp_m3q_on_uniform_grid_sg(dx, nx, f1, nwork, work1, str, ierr)
484 : ! make monotonicity preserving cubic interpolant on uniformly spaced grid
485 : use interp_1d_def
486 : use interp_1d_mp_sg
487 : real, intent(in) :: dx ! the grid spacing
488 : integer, intent(in) :: nx ! length of x vector (>= 4)
489 : real, intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
490 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
491 : real, intent(inout), pointer :: work1(:) ! =(nx, nwork)
492 : character (len=*) :: str ! for debugging
493 : integer, intent(out) :: ierr
494 0 : call m3_on_uniform_grid_sg(dx, nx, f1, quartic, .false., nwork, work1, str, ierr)
495 0 : end subroutine interp_m3q_on_uniform_grid_sg
496 :
497 :
498 : end module interp_1d_lib_sg
|
