V3FIT
ezspline_setup.f90
1 !/////
2 ! R8 !
3 !/////
4 subroutine ezspline_setup1_r8(spline_o, f, ier, exact_dim)
5  use ezspline_obj
6  implicit none
7  type(EZspline1_r8) spline_o
8  real(ezspline_r8), dimension(:), intent(in) :: f
9  ! ier:
10  ! 0=ok
11  ! 98=some error occurred in EZspline_setup
12  integer, intent(out) :: ier
13 
14  logical, intent(in), OPTIONAL :: exact_dim ! set T to require exact
15  ! dimensioning match between f and spline_o%fspl; default is F
16 
17  logical :: iexact
18  integer ifail
19 
20  integer :: ipx
21  integer iper, imsg, itol, inum, in1
22  real(ezspline_r8) ztol, df1, df2
23 
24  !-------------------------
25  iexact=.false.
26  if(present(exact_dim)) iexact = exact_dim
27 
28  if( .not.ezspline_allocated(spline_o) ) then
29  ier = 98
30  return
31  endif
32 
33  in1 = size(spline_o%fspl,2)
34 
35  ier = 57
36  if(size(f,1).lt.in1) return
37 
38  if(iexact) then
39  ier = 58
40  if(size(f,1).gt.in1) return
41  endif
42  ier = 0
43 
44  !
45  ! recompute min/max in case user changed the grid manually
46  spline_o%x1max = maxval(spline_o%x1)
47  spline_o%x1min = minval(spline_o%x1)
48  !
49  ! set xpkg, useful for cloud and array interpolation
50  imsg=0
51  itol=0 ! range tolerance option
52  ztol=5.e-7_ezspline_r8 ! range tolerance, if itol is set
53  iper=0
54  if(spline_o%ibctype1(1)==-1 .OR. spline_o%ibctype1(2)==-1) iper=1
55  call r8genxpkg(spline_o%n1, spline_o%x1(1), spline_o%x1pkg(1,1),&
56  & iper,imsg,itol,ztol, spline_o%klookup1 ,ifail)
57  if(ifail/=0) ier=27
58 
59  spline_o%isReady = 0
60 
61  spline_o%fspl(1, 1:in1) = &
62  & f(1:in1)
63 
64  if (spline_o%isHermite == 0 .and. spline_o%isLinear == 0) then
65 
66  call r8mkspline( &
67  & spline_o%x1(1), spline_o%n1, &
68  & spline_o%fspl(1,1), &
69  & spline_o%ibctype1(1), spline_o%bcval1min, &
70  & spline_o%ibctype1(2), spline_o%bcval1max, &
71  & spline_o%ilin1, ifail)
72 
73  if(ifail /= 0) then
74  ier = 98
75  else
76  spline_o%isReady = 1
77  endif
78 
79  else if (spline_o%isLinear == 1) then
80 
81  spline_o%ilin1=0
82  if(in1.gt.2) then
83  if(spline_o%x1pkg(3,4).eq.0.0_ezspline_r8) spline_o%ilin1=1 ! evenly spaced grid
84  else
85  spline_o%ilin1=1
86  end if
87 
88  spline_o%isReady = 1 ! no coefficient setup necessary
89 
90  else
91  !
92  ! Hermite polynomial coefficient setup based on Akima's method
93  ! (df/dx..etc not required)
94  !
95  ipx = 0
96  if (spline_o%ibctype1(1)==-1 .or. spline_o%ibctype1(2)==-1) then
97  ipx=1
98  else if (spline_o%ibctype1(1)<-1 .or. spline_o%ibctype1(1)>1 .or. &
99  spline_o%ibctype1(2)<-1 .or. spline_o%ibctype1(2)>1 ) then
100  ipx=99 ! an error...
101  else if (spline_o%ibctype1(1)==1 .or. spline_o%ibctype1(2)==1) then
102  ipx=2
103  if(spline_o%ibctype1(1)==1) then
104  spline_o%fspl(2,1)=spline_o%bcval1min
105  else
106  ! hand implemented default BC
107  df1=(spline_o%fspl(1,2)-spline_o%fspl(1,1))/ &
108  (spline_o%x1(2)-spline_o%x1(1))
109  df2=(spline_o%fspl(1,3)-spline_o%fspl(1,2))/ &
110  (spline_o%x1(3)-spline_o%x1(2))
111  spline_o%fspl(2,1)=(3*df1-df2)/2
112  endif
113  inum=spline_o%n1
114  if(spline_o%ibctype1(2)==1) then
115  spline_o%fspl(2,inum)=spline_o%bcval1max
116  else
117  ! hand implemented default BC
118  df1=(spline_o%fspl(1,inum)-spline_o%fspl(1,inum-1))/ &
119  (spline_o%x1(inum)-spline_o%x1(inum-1))
120  df2=(spline_o%fspl(1,inum-1)-spline_o%fspl(1,inum-2))/ &
121  (spline_o%x1(inum-1)-spline_o%x1(inum-2))
122  spline_o%fspl(2,inum)=(3*df1-df2)/2
123  endif
124  endif
125  ifail = 0
126  call r8akherm1p(spline_o%x1(1), spline_o%n1, &
127  & spline_o%fspl(1,1), &
128  & spline_o%ilin1, &
129  & ipx, ifail)
130 
131  if (ifail /=0 ) then
132  ier = 91
133  else
134  spline_o%isReady = 1
135  endif
136 
137  endif
138 
139 
140 end subroutine ezspline_setup1_r8
141 
142 
143 subroutine ezspline_setup2_r8(spline_o, f, ier, exact_dim)
144  use ezspline_obj
145  implicit none
146  type(EZspline2_r8) spline_o
147  real(ezspline_r8), dimension(:,:), intent(in) :: f
148  ! ier:
149  ! 0=ok
150  ! 98=some error occurred in EZspline_setup
151  integer, intent(out) :: ier
152 
153  logical, intent(in), OPTIONAL :: exact_dim ! set T to require exact
154  ! dimensioning match between f and spline_o%fspl; default is F
155 
156  logical :: iexact
157  integer ifail
158 
159  integer :: ipx, ipy
160  integer iper, imsg, itol, inum, ii, jj, in0, in1, in2
161  real(ezspline_r8) ztol, df1, df2
162 
163  !-------------------------
164  iexact=.false.
165  if(present(exact_dim)) iexact = exact_dim
166 
167  if( .not.ezspline_allocated(spline_o) ) then
168  ier = 98
169  return
170  endif
171 
172  in0 = size(spline_o%fspl,1)
173  in1 = size(spline_o%fspl,2)
174  in2 = size(spline_o%fspl,3)
175 
176  ier = 57
177  if(size(f,1).lt.in1) return
178  if(size(f,2).lt.in2) return
179 
180  if(iexact) then
181  ier = 58
182  if(size(f,1).gt.in1) return
183  if(size(f,2).gt.in2) return
184  endif
185  ier = 0
186 
187  !
188  ! recompute min/max in case user changed the grid manually
189  spline_o%x1max = maxval(spline_o%x1)
190  spline_o%x2max = maxval(spline_o%x2)
191  spline_o%x1min = minval(spline_o%x1)
192  spline_o%x2min = minval(spline_o%x2)
193 
194  !
195  ! set xpkg, useful for cloud and array interpolation
196  imsg=0
197  itol=0 ! range tolerance option
198  ztol=5.e-7_ezspline_r8 ! range tolerance, if itol is set
199  iper=0
200  if(spline_o%ibctype1(1)==-1 .OR. spline_o%ibctype1(2)==-1) iper=1
201  call r8genxpkg(spline_o%n1,spline_o%x1(1),spline_o%x1pkg(1,1),&
202  & iper,imsg,itol,ztol,spline_o%klookup1,ifail)
203  if(ifail/=0) ier=27
204  iper=0
205  if(spline_o%ibctype2(1)==-1 .OR. spline_o%ibctype2(2)==-1) iper=1
206  call r8genxpkg(spline_o%n2,spline_o%x2(1),spline_o%x2pkg(1,1),&
207  & iper,imsg,itol,ztol,spline_o%klookup2,ifail)
208  if(ifail/=0) ier=27
209 
210  spline_o%isReady = 0
211 
212  spline_o%fspl(1, 1:in1, 1:in2) = &
213  & f(1:in1, 1:in2)
214 
215  ! this fixes a VMS f90 compiler optimizer problem:
216  if(ztol.eq.-1.2345d30) &
217  write(6,*) 'spline_o%fspl(1,1,1) = ', spline_o%fspl(1,1,1)
218 
219  if (spline_o%isHybrid == 1) then
220 
221  call r8mkintrp2d( &
222  & spline_o%x1(1), spline_o%n1, &
223  & spline_o%x2(1), spline_o%n2, &
224  & spline_o%hspline, spline_o%fspl(1,1,1), &
225  & in0,in1,in2, &
226  & spline_o%ibctype1(1), spline_o%bcval1min(1), &
227  & spline_o%ibctype1(2), spline_o%bcval1max(1), &
228  & spline_o%ibctype2(1), spline_o%bcval2min(1), &
229  & spline_o%ibctype2(2), spline_o%bcval2max(1), &
230  & ifail)
231 
232  spline_o%ilin1 = 0 ! Hybrid does not compute this
233  spline_o%ilin2 = 0 ! Hybrid does not compute this
234 
235  if(ifail /= 0) then
236  ier = 98
237  else
238  spline_o%isReady = 1
239  endif
240 
241  else if (spline_o%isHermite == 0 .and. spline_o%isLinear == 0) then
242 
243  call r8mkbicub( &
244  & spline_o%x1(1), spline_o%n1, &
245  & spline_o%x2(1), spline_o%n2, &
246  & spline_o%fspl(1,1,1), spline_o%n1, &
247  & spline_o%ibctype1(1), spline_o%bcval1min(1), &
248  & spline_o%ibctype1(2), spline_o%bcval1max(1), &
249  & spline_o%ibctype2(1), spline_o%bcval2min(1), &
250  & spline_o%ibctype2(2), spline_o%bcval2max(1), &
251  & spline_o%ilin1, spline_o%ilin2, ifail)
252 
253  if(ifail /= 0) then
254  ier = 98
255  else
256  spline_o%isReady = 1
257  endif
258 
259  else if (spline_o%isLinear == 1) then
260 
261  spline_o%ilin1=0
262  if(in1.gt.2) then
263  if(spline_o%x1pkg(3,4).eq.0.0_ezspline_r8) spline_o%ilin1=1 ! evenly spaced grid
264  else
265  spline_o%ilin1=1
266  end if
267 
268  spline_o%ilin2=0
269  if(in2.gt.2) then
270  if(spline_o%x2pkg(3,4).eq.0.0_ezspline_r8) spline_o%ilin2=1 ! evenly spaced grid
271  else
272  spline_o%ilin2=1
273  end if
274 
275  spline_o%isReady = 1 ! no coefficient setup necessary
276 
277  else
278  !
279  ! Hermite polynomial coefficient setup based on Akima's method
280  ! (df/dx..etc not required)
281  !
282  ipx = 0
283  if (spline_o%ibctype1(1)==-1 .or. spline_o%ibctype1(2)==-1) then
284  ipx=1
285  else if (spline_o%ibctype1(1)<-1 .or. spline_o%ibctype1(1)>1 .or. &
286  spline_o%ibctype1(2)<-1 .or. spline_o%ibctype1(2)>1 ) then
287  ipx=99 ! an error...
288  else if (spline_o%ibctype1(1)==1 .or. spline_o%ibctype1(2)==1) then
289  ipx=2
290  do jj=1,spline_o%n2
291  if(spline_o%ibctype1(1)==1) then
292  spline_o%fspl(2,1,jj)=spline_o%bcval1min(jj)
293  else
294  ! hand implemented default BC
295  df1=(spline_o%fspl(1,2,jj)-spline_o%fspl(1,1,jj))/ &
296  (spline_o%x1(2)-spline_o%x1(1))
297  df2=(spline_o%fspl(1,3,jj)-spline_o%fspl(1,2,jj))/ &
298  (spline_o%x1(3)-spline_o%x1(2))
299  spline_o%fspl(2,1,jj)=(3*df1-df2)/2
300  endif
301  inum=spline_o%n1
302  if(spline_o%ibctype1(2)==1) then
303  spline_o%fspl(2,inum,jj)=spline_o%bcval1max(jj)
304  else
305  ! hand implemented default BC
306  df1=(spline_o%fspl(1,inum,jj)-spline_o%fspl(1,inum-1,jj))/ &
307  (spline_o%x1(inum)-spline_o%x1(inum-1))
308  df2=(spline_o%fspl(1,inum-1,jj)-spline_o%fspl(1,inum-2,jj))/ &
309  (spline_o%x1(inum-1)-spline_o%x1(inum-2))
310  spline_o%fspl(2,inum,jj)=(3*df1-df2)/2
311  endif
312  enddo
313  endif
314  ipy = 0
315  if (spline_o%ibctype2(1)==-1 .or. spline_o%ibctype2(2)==-1) then
316  ipy=1
317  else if (spline_o%ibctype2(1)<-1 .or. spline_o%ibctype2(1)>1 .or. &
318  spline_o%ibctype2(2)<-1 .or. spline_o%ibctype2(2)>1 ) then
319  ipy=99 ! an error...
320  else if (spline_o%ibctype2(1)==1 .or. spline_o%ibctype2(2)==1) then
321  ipy=2
322  do ii=1,spline_o%n1
323  if(spline_o%ibctype2(1)==1) then
324  spline_o%fspl(3,ii,1)=spline_o%bcval2min(ii)
325  else
326  ! hand implemented default BC
327  df1=(spline_o%fspl(1,ii,2)-spline_o%fspl(1,ii,1))/ &
328  (spline_o%x2(2)-spline_o%x2(1))
329  df2=(spline_o%fspl(1,ii,3)-spline_o%fspl(1,ii,2))/ &
330  (spline_o%x2(3)-spline_o%x2(2))
331  spline_o%fspl(3,ii,1)=(3*df1-df2)/2
332  endif
333  inum=spline_o%n2
334  if(spline_o%ibctype2(2)==1) then
335  spline_o%fspl(3,ii,inum)=spline_o%bcval2max(ii)
336  else
337  ! hand implemented default BC
338  df1=(spline_o%fspl(1,ii,inum)-spline_o%fspl(1,ii,inum-1))/ &
339  (spline_o%x2(inum)-spline_o%x2(inum-1))
340  df2=(spline_o%fspl(1,ii,inum-1)-spline_o%fspl(1,ii,inum-2))/ &
341  (spline_o%x2(inum-1)-spline_o%x2(inum-2))
342  spline_o%fspl(3,ii,inum)=(3*df1-df2)/2
343  endif
344  enddo
345  endif
346  ifail = 0
347  call r8akherm2p(spline_o%x1(1), spline_o%n1, &
348  & spline_o%x2(1), spline_o%n2, &
349  & spline_o%fspl(1,1,1), spline_o%n1, &
350  & spline_o%ilin1, spline_o%ilin2, &
351  & ipx,ipy, ifail)
352 
353  if (ifail /=0 ) then
354  ier = 91
355  else
356  spline_o%isReady = 1
357  endif
358 
359  endif
360 
361 
362 end subroutine ezspline_setup2_r8
363 
364 
365 
366 subroutine ezspline_setup3_r8(spline_o, f, ier, exact_dim)
367  use ezspline_obj
368  implicit none
369  type(EZspline3_r8) spline_o
370  real(ezspline_r8), dimension(:,:,:), intent(in) :: f
371  ! ier:
372  ! 0=ok
373  ! 98=some error occurred in EZspline_setup
374  integer, intent(out) :: ier
375 
376  logical, intent(in), OPTIONAL :: exact_dim ! set T to require exact
377  ! dimensioning match between f and spline_o%fspl; default is F
378 
379  logical :: iexact
380  integer ifail
381 
382  integer :: ipx, ipy, ipz
383  integer iper, imsg, itol, inum, ii, jj, in0, in1, in2, in3
384  real(ezspline_r8) ztol, df1, df2
385 
386  !-------------------------
387  iexact=.false.
388  if(present(exact_dim)) iexact = exact_dim
389 
390  if( .not.ezspline_allocated(spline_o) ) then
391  ier = 98
392  return
393  endif
394 
395  in0 = size(spline_o%fspl,1)
396  in1 = size(spline_o%fspl,2)
397  in2 = size(spline_o%fspl,3)
398  in3 = size(spline_o%fspl,4)
399 
400  ier = 57
401  if(size(f,1).lt.in1) return
402  if(size(f,2).lt.in2) return
403  if(size(f,3).lt.in3) return
404 
405  if(iexact) then
406  ier = 58
407  if(size(f,1).gt.in1) return
408  if(size(f,2).gt.in2) return
409  if(size(f,3).gt.in3) return
410  endif
411  ier = 0
412 
413  !
414  ! recompute min/max in case user changed the grid manually
415  spline_o%x1max = maxval(spline_o%x1)
416  spline_o%x2max = maxval(spline_o%x2)
417  spline_o%x3max = maxval(spline_o%x3)
418  spline_o%x1min = minval(spline_o%x1)
419  spline_o%x2min = minval(spline_o%x2)
420  spline_o%x3min = minval(spline_o%x3)
421 
422  !
423  ! set xpkg, useful for cloud and array interpolation
424  imsg=0
425  itol=0 ! range tolerance option
426  ztol=5.e-7_ezspline_r8 ! range tolerance, if itol is set
427  iper=0
428  if(spline_o%ibctype1(1)==-1 .OR. spline_o%ibctype1(2)==-1) iper=1
429  call r8genxpkg(spline_o%n1,spline_o%x1(1),spline_o%x1pkg(1,1),&
430  & iper,imsg,itol,ztol,spline_o%klookup1,ifail)
431  if(ifail/=0) ier=27
432  iper=0
433  if(spline_o%ibctype2(1)==-1 .OR. spline_o%ibctype2(2)==-1) iper=1
434  call r8genxpkg(spline_o%n2,spline_o%x2(1),spline_o%x2pkg(1,1),&
435  & iper,imsg,itol,ztol,spline_o%klookup2,ifail)
436  if(ifail/=0) ier=27
437  iper=0
438  if(spline_o%ibctype3(1)==-1 .OR. spline_o%ibctype3(2)==-1) iper=1
439  call r8genxpkg(spline_o%n3,spline_o%x3(1),spline_o%x3pkg(1,1),&
440  & iper,imsg,itol,ztol,spline_o%klookup3,ifail)
441  if(ifail/=0) ier=27
442 
443  spline_o%isReady = 0
444 
445  spline_o%fspl(1, 1:in1, 1:in2, 1:in3) = &
446  & f(1:in1, 1:in2, 1:in3)
447 
448  if (spline_o%isHybrid == 1) then
449 
450  call r8mkintrp3d( &
451  & spline_o%x1(1), spline_o%n1, &
452  & spline_o%x2(1), spline_o%n2, &
453  & spline_o%x3(1), spline_o%n3, &
454  & spline_o%hspline, spline_o%fspl(1,1,1,1), &
455  & in0,in1,in2,in3, &
456  & spline_o%ibctype1(1), spline_o%bcval1min(1,1), &
457  & spline_o%ibctype1(2), spline_o%bcval1max(1,1), &
458  & spline_o%ibctype2(1), spline_o%bcval2min(1,1), &
459  & spline_o%ibctype2(2), spline_o%bcval2max(1,1), &
460  & spline_o%ibctype3(1), spline_o%bcval3min(1,1), &
461  & spline_o%ibctype3(2), spline_o%bcval3max(1,1), &
462  & ifail)
463 
464  spline_o%ilin1 = 0 ! Hybrid does not compute this
465  spline_o%ilin2 = 0 ! Hybrid does not compute this
466  spline_o%ilin3 = 0 ! Hybrid does not compute this
467 
468  if(ifail /= 0) then
469  ier = 98
470  else
471  spline_o%isReady = 1
472  endif
473 
474  else if (spline_o%isHermite == 0 .and. spline_o%isLinear == 0) then
475 
476  call r8mktricub( &
477  & spline_o%x1(1), spline_o%n1, &
478  & spline_o%x2(1), spline_o%n2, &
479  & spline_o%x3(1), spline_o%n3, &
480  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
481  & spline_o%ibctype1(1), spline_o%bcval1min(1,1), &
482  & spline_o%ibctype1(2), spline_o%bcval1max(1,1), spline_o%n2, &
483  & spline_o%ibctype2(1), spline_o%bcval2min(1,1), &
484  & spline_o%ibctype2(2), spline_o%bcval2max(1,1), spline_o%n1, &
485  & spline_o%ibctype3(1), spline_o%bcval3min(1,1), &
486  & spline_o%ibctype3(2), spline_o%bcval3max(1,1), spline_o%n1, &
487  & spline_o%ilin1, spline_o%ilin2, spline_o%ilin3, ifail)
488 
489  if(ifail /= 0) then
490  ier = 98
491  else
492  spline_o%isReady = 1
493  endif
494 
495  else if (spline_o%isLinear == 1) then
496 
497  spline_o%ilin1=0
498  if(in1.gt.2) then
499  if(spline_o%x1pkg(3,4).eq.0.0_ezspline_r8) spline_o%ilin1=1 ! evenly spaced grid
500  else
501  spline_o%ilin1=1
502  end if
503 
504  spline_o%ilin2=0
505  if(in2.gt.2) then
506  if(spline_o%x2pkg(3,4).eq.0.0_ezspline_r8) spline_o%ilin2=1 ! evenly spaced grid
507  else
508  spline_o%ilin2=1
509  end if
510 
511  spline_o%ilin3=0
512  if(in3.gt.2) then
513  if(spline_o%x3pkg(3,4).eq.0.0_ezspline_r8) spline_o%ilin3=1 ! evenly spaced grid
514  else
515  spline_o%ilin3=1
516  end if
517 
518  spline_o%isReady = 1 ! no coefficient setup necessary
519 
520  else
521  !
522  ! Hermite polynomial coefficient setup based on Akima's method
523  ! (df/dx..etc not required)
524  !
525  ipx = 0
526  if (spline_o%ibctype1(1)==-1 .or. spline_o%ibctype1(2)==-1) then
527  ipx=1
528  else if (spline_o%ibctype1(1)<-1 .or. spline_o%ibctype1(1)>1 .or. &
529  spline_o%ibctype1(2)<-1 .or. spline_o%ibctype1(2)>1 ) then
530  ipx=99 ! an error...
531  else if (spline_o%ibctype1(1)==1 .or. spline_o%ibctype1(2)==1) then
532  ipx=2
533  do jj=1,spline_o%n3
534  do ii=1,spline_o%n2
535  if(spline_o%ibctype1(1)==1) then
536  spline_o%fspl(2,1,ii,jj)=spline_o%bcval1min(ii,jj)
537  else
538  ! hand implemented default BC
539  df1=(spline_o%fspl(1,2,ii,jj)-spline_o%fspl(1,1,ii,jj))/ &
540  (spline_o%x1(2)-spline_o%x1(1))
541  df2=(spline_o%fspl(1,3,ii,jj)-spline_o%fspl(1,2,ii,jj))/ &
542  (spline_o%x1(3)-spline_o%x1(2))
543  spline_o%fspl(2,1,ii,jj)=(3*df1-df2)/2
544  endif
545  inum=spline_o%n1
546  if(spline_o%ibctype1(2)==1) then
547  spline_o%fspl(2,inum,ii,jj)=spline_o%bcval1max(ii,jj)
548  else
549  ! hand implemented default BC
550  df1=(spline_o%fspl(1,inum,ii,jj)-spline_o%fspl(1,inum-1,ii,jj))/ &
551  (spline_o%x1(inum)-spline_o%x1(inum-1))
552  df2=(spline_o%fspl(1,inum-1,ii,jj)-spline_o%fspl(1,inum-2,ii,jj))/ &
553  (spline_o%x1(inum-1)-spline_o%x1(inum-2))
554  spline_o%fspl(2,inum,ii,jj)=(3*df1-df2)/2
555  endif
556  enddo
557  enddo
558  endif
559  ipy = 0
560  if (spline_o%ibctype2(1)==-1 .or. spline_o%ibctype2(2)==-1) then
561  ipy=1
562  else if (spline_o%ibctype2(1)<-1 .or. spline_o%ibctype2(1)>1 .or. &
563  spline_o%ibctype2(2)<-1 .or. spline_o%ibctype2(2)>1 ) then
564  ipy=99 ! an error...
565  else if (spline_o%ibctype2(1)==1 .or. spline_o%ibctype2(2)==1) then
566  ipy=2
567  do jj=1,spline_o%n3
568  do ii=1,spline_o%n1
569  if(spline_o%ibctype2(1)==1) then
570  spline_o%fspl(3,ii,1,jj)=spline_o%bcval2min(ii,jj)
571  else
572  ! hand implemented default BC
573  df1=(spline_o%fspl(1,ii,2,jj)-spline_o%fspl(1,ii,1,jj))/ &
574  (spline_o%x2(2)-spline_o%x2(1))
575  df2=(spline_o%fspl(1,ii,3,jj)-spline_o%fspl(1,ii,2,jj))/ &
576  (spline_o%x2(3)-spline_o%x2(2))
577  spline_o%fspl(3,ii,1,jj)=(3*df1-df2)/2
578  endif
579  inum=spline_o%n2
580  if(spline_o%ibctype2(2)==1) then
581  spline_o%fspl(3,ii,inum,jj)=spline_o%bcval2max(ii,jj)
582  else
583  ! hand implemented default BC
584  df1=(spline_o%fspl(1,ii,inum,jj)-spline_o%fspl(1,ii,inum-1,jj))/ &
585  (spline_o%x2(inum)-spline_o%x2(inum-1))
586  df2=(spline_o%fspl(1,ii,inum-1,jj)-spline_o%fspl(1,ii,inum-2,jj))/ &
587  (spline_o%x2(inum-1)-spline_o%x2(inum-2))
588  spline_o%fspl(3,ii,inum,jj)=(3*df1-df2)/2
589  endif
590  enddo
591  enddo
592  endif
593  ipz = 0
594  if (spline_o%ibctype3(1)==-1 .or. spline_o%ibctype3(2)==-1) then
595  ipz=1
596  else if (spline_o%ibctype3(1)<-1 .or. spline_o%ibctype3(1)>1 .or. &
597  spline_o%ibctype3(2)<-1 .or. spline_o%ibctype3(2)>1 ) then
598  ipz=99 ! an error...
599  else if (spline_o%ibctype3(1)==1 .or. spline_o%ibctype3(2)==1) then
600  ipz=2
601  do jj=1,spline_o%n2
602  do ii=1,spline_o%n1
603  if(spline_o%ibctype3(1)==1) then
604  spline_o%fspl(4,ii,jj,1)=spline_o%bcval3min(ii,jj)
605  else
606  ! hand implemented default BC
607  df1=(spline_o%fspl(1,ii,jj,2)-spline_o%fspl(1,ii,jj,1))/ &
608  (spline_o%x3(2)-spline_o%x3(1))
609  df2=(spline_o%fspl(1,ii,jj,3)-spline_o%fspl(1,ii,jj,2))/ &
610  (spline_o%x3(3)-spline_o%x3(2))
611  spline_o%fspl(4,ii,jj,1)=(3*df1-df2)/2
612  endif
613  inum=spline_o%n3
614  if(spline_o%ibctype3(2)==1) then
615  spline_o%fspl(4,ii,jj,inum)=spline_o%bcval3max(ii,jj)
616  else
617  ! hand implemented default BC
618  df1=(spline_o%fspl(1,ii,jj,inum)-spline_o%fspl(1,ii,jj,inum-1))/ &
619  (spline_o%x3(inum)-spline_o%x3(inum-1))
620  df2=(spline_o%fspl(1,ii,jj,inum-1)-spline_o%fspl(1,ii,jj,inum-2))/ &
621  (spline_o%x3(inum-1)-spline_o%x3(inum-2))
622  spline_o%fspl(4,ii,jj,inum)=(3*df1-df2)/2
623  endif
624  enddo
625  enddo
626  endif
627  ifail = 0
628  call r8akherm3p(spline_o%x1(1), spline_o%n1, &
629  & spline_o%x2(1), spline_o%n2, &
630  & spline_o%x3(1), spline_o%n3, &
631  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
632  & spline_o%ilin1, spline_o%ilin2, spline_o%ilin3, &
633  & ipx,ipy,ipz, ifail)
634 
635  if (ifail /=0 ) then
636  ier = 91
637  else
638  spline_o%isReady = 1
639  endif
640 
641  endif
642 
643 
644 end subroutine ezspline_setup3_r8
645 !/////
646 ! R4 !
647 !/////
648 subroutine ezspline_setup1_r4(spline_o, f, ier, exact_dim)
649  use ezspline_obj
650  implicit none
651  type(EZspline1_r4) spline_o
652  real(ezspline_r4), dimension(:), intent(in) :: f
653  ! ier:
654  ! 0=ok
655  ! 98=some error occurred in EZspline_setup
656  integer, intent(out) :: ier
657 
658  logical, intent(in), OPTIONAL :: exact_dim ! set T to require exact
659  ! dimensioning match between f and spline_o%fspl; default is F
660 
661  logical :: iexact
662  integer ifail
663 
664  integer :: ipx
665  integer iper, imsg, itol, inum, in1
666  real(ezspline_r4) ztol, df1, df2
667 
668  !-------------------------
669  iexact=.false.
670  if(present(exact_dim)) iexact = exact_dim
671 
672  if( .not.ezspline_allocated(spline_o) ) then
673  ier = 98
674  return
675  endif
676 
677  in1 = size(spline_o%fspl,2)
678 
679  ier = 57
680  if(size(f,1).lt.in1) return
681 
682  if(iexact) then
683  ier = 58
684  if(size(f,1).gt.in1) return
685  endif
686  ier = 0
687 
688  !
689  ! recompute min/max in case user changed the grid manually
690  spline_o%x1max = maxval(spline_o%x1)
691  spline_o%x1min = minval(spline_o%x1)
692 
693  !
694  ! set xpkg, useful for cloud and array interpolation
695  imsg=0
696  itol=0 ! range tolerance option
697  ztol=5.e-7_ezspline_r4 ! range tolerance, if itol is set
698  iper=0
699  if(spline_o%ibctype1(1)==-1 .OR. spline_o%ibctype1(2)==-1) iper=1
700  call genxpkg(spline_o%n1, spline_o%x1(1), spline_o%x1pkg(1,1),&
701  & iper,imsg,itol,ztol,spline_o%klookup1,ifail)
702  if(ifail/=0) ier=27
703 
704  spline_o%isReady = 0
705 
706  spline_o%fspl(1, 1:in1) = &
707  & f(1:in1)
708 
709  if (spline_o%isHermite == 0 .and. spline_o%isLinear == 0) then
710 
711  call mkspline( &
712  & spline_o%x1(1), spline_o%n1, &
713  & spline_o%fspl(1,1), &
714  & spline_o%ibctype1(1), spline_o%bcval1min, &
715  & spline_o%ibctype1(2), spline_o%bcval1max, &
716  & spline_o%ilin1, ifail)
717 
718  if(ifail /= 0) then
719  ier = 98
720  else
721  spline_o%isReady = 1
722  endif
723 
724  else if (spline_o%isLinear == 1) then
725 
726  spline_o%ilin1=0
727  if(in1.gt.2) then
728  if(spline_o%x1pkg(3,4).eq.0.0_ezspline_r4) spline_o%ilin1=1 ! evenly spaced grid
729  else
730  spline_o%ilin1=1
731  end if
732 
733  spline_o%isReady = 1 ! no coefficient setup necessary
734 
735  else
736  !
737  ! Hermite polynomial coefficient setup based on Akima's method
738  ! (df/dx..etc not required)
739  !
740  ipx = 0
741  if (spline_o%ibctype1(1)==-1 .or. spline_o%ibctype1(2)==-1) then
742  ipx=1
743  else if (spline_o%ibctype1(1)<-1 .or. spline_o%ibctype1(1)>1 .or. &
744  spline_o%ibctype1(2)<-1 .or. spline_o%ibctype1(2)>1 ) then
745  ipx=99 ! an error...
746  else if (spline_o%ibctype1(1)==1 .or. spline_o%ibctype1(2)==1) then
747  ipx=2
748  if(spline_o%ibctype1(1)==1) then
749  spline_o%fspl(2,1)=spline_o%bcval1min
750  else
751  ! hand implemented default BC
752  df1=(spline_o%fspl(1,2)-spline_o%fspl(1,1))/ &
753  (spline_o%x1(2)-spline_o%x1(1))
754  df2=(spline_o%fspl(1,3)-spline_o%fspl(1,2))/ &
755  (spline_o%x1(3)-spline_o%x1(2))
756  spline_o%fspl(2,1)=(3*df1-df2)/2
757  endif
758  inum=spline_o%n1
759  if(spline_o%ibctype1(2)==1) then
760  spline_o%fspl(2,inum)=spline_o%bcval1max
761  else
762  ! hand implemented default BC
763  df1=(spline_o%fspl(1,inum)-spline_o%fspl(1,inum-1))/ &
764  (spline_o%x1(inum)-spline_o%x1(inum-1))
765  df2=(spline_o%fspl(1,inum-1)-spline_o%fspl(1,inum-2))/ &
766  (spline_o%x1(inum-1)-spline_o%x1(inum-2))
767  spline_o%fspl(2,inum)=(3*df1-df2)/2
768  endif
769  endif
770  ifail = 0
771  call akherm1p(spline_o%x1(1), spline_o%n1, &
772  & spline_o%fspl(1,1), &
773  & spline_o%ilin1, &
774  & ipx, ifail)
775 
776  if (ifail /=0 ) then
777  ier = 91
778  else
779  spline_o%isReady = 1
780  endif
781 
782  endif
783 
784 
785 end subroutine ezspline_setup1_r4
786 
787 
788 subroutine ezspline_setup2_r4(spline_o, f, ier, exact_dim)
789  use ezspline_obj
790  implicit none
791  type(EZspline2_r4) spline_o
792  real(ezspline_r4), dimension(:,:), intent(in) :: f
793  ! ier:
794  ! 0=ok
795  ! 98=some error occurred in EZspline_setup
796  integer, intent(out) :: ier
797 
798  logical, intent(in), OPTIONAL :: exact_dim ! set T to require exact
799  ! dimensioning match between f and spline_o%fspl; default is F
800 
801  logical :: iexact
802  integer ifail
803 
804  integer :: ipx, ipy
805  integer iper, imsg, itol, inum, ii, jj, in0, in1, in2
806  real(ezspline_r4) ztol, df1, df2
807 
808  !-------------------------
809  iexact=.false.
810  if(present(exact_dim)) iexact = exact_dim
811 
812  if( .not.ezspline_allocated(spline_o) ) then
813  ier = 98
814  return
815  endif
816 
817  in0 = size(spline_o%fspl,1)
818  in1 = size(spline_o%fspl,2)
819  in2 = size(spline_o%fspl,3)
820 
821  ier = 57
822  if(size(f,1).lt.in1) return
823  if(size(f,2).lt.in2) return
824 
825  if(iexact) then
826  ier = 58
827  if(size(f,1).gt.in1) return
828  if(size(f,2).gt.in2) return
829  endif
830  ier = 0
831 
832  !
833  ! recompute min/max in case user changed the grid manually
834  spline_o%x1max = maxval(spline_o%x1)
835  spline_o%x2max = maxval(spline_o%x2)
836  spline_o%x1min = minval(spline_o%x1)
837  spline_o%x2min = minval(spline_o%x2)
838  !
839  ! set xpkg, useful for cloud and array interpolation
840  imsg=0
841  itol=0 ! range tolerance option
842  ztol=5.e-7_ezspline_r4 ! range tolerance, if itol is set
843  iper=0
844  if(spline_o%ibctype1(1)==-1 .OR. spline_o%ibctype1(2)==-1) iper=1
845  call genxpkg(spline_o%n1,spline_o%x1(1),spline_o%x1pkg(1,1),&
846  & iper,imsg,itol,ztol,spline_o%klookup1,ifail)
847  if(ifail/=0) ier=27
848  iper=0
849  if(spline_o%ibctype2(1)==-1 .OR. spline_o%ibctype2(2)==-1) iper=1
850  call genxpkg(spline_o%n2,spline_o%x2(1),spline_o%x2pkg(1,1),&
851  & iper,imsg,itol,ztol,spline_o%klookup2,ifail)
852  if(ifail/=0) ier=27
853 
854  spline_o%isReady = 0
855 
856  spline_o%fspl(1, 1:in1, 1:in2) = &
857  & f(1:in1, 1:in2)
858 
859  if (spline_o%isHybrid == 1) then
860 
861  call mkintrp2d( &
862  & spline_o%x1(1), spline_o%n1, &
863  & spline_o%x2(1), spline_o%n2, &
864  & spline_o%hspline, spline_o%fspl(1,1,1), &
865  & in0,in1,in2, &
866  & spline_o%ibctype1(1), spline_o%bcval1min(1), &
867  & spline_o%ibctype1(2), spline_o%bcval1max(1), &
868  & spline_o%ibctype2(1), spline_o%bcval2min(1), &
869  & spline_o%ibctype2(2), spline_o%bcval2max(1), &
870  & ifail)
871 
872  spline_o%ilin1 = 0 ! Hybrid does not compute this
873  spline_o%ilin2 = 0 ! Hybrid does not compute this
874 
875  if(ifail /= 0) then
876  ier = 98
877  else
878  spline_o%isReady = 1
879  endif
880 
881  else if (spline_o%isHermite == 0 .and. spline_o%isLinear == 0) then
882 
883  call mkbicub( &
884  & spline_o%x1(1), spline_o%n1, &
885  & spline_o%x2(1), spline_o%n2, &
886  & spline_o%fspl(1,1,1), spline_o%n1, &
887  & spline_o%ibctype1(1), spline_o%bcval1min(1), &
888  & spline_o%ibctype1(2), spline_o%bcval1max(1), &
889  & spline_o%ibctype2(1), spline_o%bcval2min(1), &
890  & spline_o%ibctype2(2), spline_o%bcval2max(1), &
891  & spline_o%ilin1, spline_o%ilin2, ifail)
892 
893  if(ifail /= 0) then
894  ier = 98
895  else
896  spline_o%isReady = 1
897  endif
898 
899  else if (spline_o%isLinear == 1) then
900 
901  spline_o%ilin1=0
902  if(in1.gt.2) then
903  if(spline_o%x1pkg(3,4).eq.0.0_ezspline_r4) spline_o%ilin1=1 ! evenly spaced grid
904  else
905  spline_o%ilin1=1
906  end if
907 
908  spline_o%ilin2=0
909  if(in2.gt.2) then
910  if(spline_o%x2pkg(3,4).eq.0.0_ezspline_r4) spline_o%ilin2=1 ! evenly spaced grid
911  else
912  spline_o%ilin2=1
913  end if
914 
915  spline_o%isReady = 1 ! no coefficient setup necessary
916 
917  else
918  !
919  ! Hermite polynomial coefficient setup based on Akima's method
920  ! (df/dx..etc not required)
921  !
922  ipx = 0
923  if (spline_o%ibctype1(1)==-1 .or. spline_o%ibctype1(2)==-1) then
924  ipx=1
925  else if (spline_o%ibctype1(1)<-1 .or. spline_o%ibctype1(1)>1 .or. &
926  spline_o%ibctype1(2)<-1 .or. spline_o%ibctype1(2)>1 ) then
927  ipx=99 ! an error...
928  else if (spline_o%ibctype1(1)==1 .or. spline_o%ibctype1(2)==1) then
929  ipx=2
930  do jj=1,spline_o%n2
931  if(spline_o%ibctype1(1)==1) then
932  spline_o%fspl(2,1,jj)=spline_o%bcval1min(jj)
933  else
934  ! hand implemented default BC
935  df1=(spline_o%fspl(1,2,jj)-spline_o%fspl(1,1,jj))/ &
936  (spline_o%x1(2)-spline_o%x1(1))
937  df2=(spline_o%fspl(1,3,jj)-spline_o%fspl(1,2,jj))/ &
938  (spline_o%x1(3)-spline_o%x1(2))
939  spline_o%fspl(2,1,jj)=(3*df1-df2)/2
940  endif
941  inum=spline_o%n1
942  if(spline_o%ibctype1(2)==1) then
943  spline_o%fspl(2,inum,jj)=spline_o%bcval1max(jj)
944  else
945  ! hand implemented default BC
946  df1=(spline_o%fspl(1,inum,jj)-spline_o%fspl(1,inum-1,jj))/ &
947  (spline_o%x1(inum)-spline_o%x1(inum-1))
948  df2=(spline_o%fspl(1,inum-1,jj)-spline_o%fspl(1,inum-2,jj))/ &
949  (spline_o%x1(inum-1)-spline_o%x1(inum-2))
950  spline_o%fspl(2,inum,jj)=(3*df1-df2)/2
951  endif
952  enddo
953  endif
954  ipy = 0
955  if (spline_o%ibctype2(1)==-1 .or. spline_o%ibctype2(2)==-1) then
956  ipy=1
957  else if (spline_o%ibctype2(1)<-1 .or. spline_o%ibctype2(1)>1 .or. &
958  spline_o%ibctype2(2)<-1 .or. spline_o%ibctype2(2)>1 ) then
959  ipy=99 ! an error...
960  else if (spline_o%ibctype2(1)==1 .or. spline_o%ibctype2(2)==1) then
961  ipy=2
962  do ii=1,spline_o%n1
963  if(spline_o%ibctype2(1)==1) then
964  spline_o%fspl(3,ii,1)=spline_o%bcval2min(ii)
965  else
966  ! hand implemented default BC
967  df1=(spline_o%fspl(1,ii,2)-spline_o%fspl(1,ii,1))/ &
968  (spline_o%x2(2)-spline_o%x2(1))
969  df2=(spline_o%fspl(1,ii,3)-spline_o%fspl(1,ii,2))/ &
970  (spline_o%x2(3)-spline_o%x2(2))
971  spline_o%fspl(3,ii,1)=(3*df1-df2)/2
972  endif
973  inum=spline_o%n2
974  if(spline_o%ibctype2(2)==1) then
975  spline_o%fspl(3,ii,inum)=spline_o%bcval2max(ii)
976  else
977  ! hand implemented default BC
978  df1=(spline_o%fspl(1,ii,inum)-spline_o%fspl(1,ii,inum-1))/ &
979  (spline_o%x2(inum)-spline_o%x2(inum-1))
980  df2=(spline_o%fspl(1,ii,inum-1)-spline_o%fspl(1,ii,inum-2))/ &
981  (spline_o%x2(inum-1)-spline_o%x2(inum-2))
982  spline_o%fspl(3,ii,inum)=(3*df1-df2)/2
983  endif
984  enddo
985  endif
986  ifail = 0
987  call akherm2p(spline_o%x1(1), spline_o%n1, &
988  & spline_o%x2(1), spline_o%n2, &
989  & spline_o%fspl(1,1,1), spline_o%n1, &
990  & spline_o%ilin1, spline_o%ilin2, &
991  & ipx,ipy, ifail)
992 
993  if (ifail /=0 ) then
994  ier = 91
995  else
996  spline_o%isReady = 1
997  endif
998 
999  endif
1000 
1001 
1002 end subroutine ezspline_setup2_r4
1003 
1004 
1005 
1006 subroutine ezspline_setup3_r4(spline_o, f, ier, exact_dim)
1007  use ezspline_obj
1008  implicit none
1009  type(EZspline3_r4) spline_o
1010  real(ezspline_r4), dimension(:,:,:), intent(in) :: f
1011  ! ier:
1012  ! 0=ok
1013  ! 98=some error occurred in EZspline_setup
1014  integer, intent(out) :: ier
1015 
1016  logical, intent(in), OPTIONAL :: exact_dim ! set T to require exact
1017  ! dimensioning match between f and spline_o%fspl; default is F
1018 
1019  logical :: iexact
1020  integer ifail
1021 
1022  integer :: ipx, ipy, ipz
1023  integer iper, imsg, itol, inum, ii, jj, in0, in1, in2, in3
1024  real(ezspline_r4) ztol, df1, df2
1025 
1026  !-------------------------
1027  iexact=.false.
1028  if(present(exact_dim)) iexact = exact_dim
1029 
1030  if( .not.ezspline_allocated(spline_o) ) then
1031  ier = 98
1032  return
1033  endif
1034 
1035  in0 = size(spline_o%fspl,1)
1036  in1 = size(spline_o%fspl,2)
1037  in2 = size(spline_o%fspl,3)
1038  in3 = size(spline_o%fspl,4)
1039 
1040  ier = 57
1041  if(size(f,1).lt.in1) return
1042  if(size(f,2).lt.in2) return
1043  if(size(f,3).lt.in3) return
1044 
1045  if(iexact) then
1046  ier = 58
1047  if(size(f,1).gt.in1) return
1048  if(size(f,2).gt.in2) return
1049  if(size(f,3).gt.in3) return
1050  endif
1051  ier = 0
1052 
1053  !
1054  ! recompute min/max in case user changed the grid manually
1055  spline_o%x1max = maxval(spline_o%x1)
1056  spline_o%x2max = maxval(spline_o%x2)
1057  spline_o%x3max = maxval(spline_o%x3)
1058  spline_o%x1min = minval(spline_o%x1)
1059  spline_o%x2min = minval(spline_o%x2)
1060  spline_o%x3min = minval(spline_o%x3)
1061 
1062  !
1063  ! set xpkg, useful for cloud and array interpolation
1064  imsg=0
1065  itol=0 ! range tolerance option
1066  ztol=5.e-7_ezspline_r4 ! range tolerance, if itol is set
1067  iper=0
1068  if(spline_o%ibctype1(1)==-1 .OR. spline_o%ibctype1(2)==-1) iper=1
1069  call genxpkg(spline_o%n1,spline_o%x1(1),spline_o%x1pkg(1,1),&
1070  & iper,imsg,itol,ztol,spline_o%klookup1,ifail)
1071  if(ifail/=0) ier=27
1072  iper=0
1073  if(spline_o%ibctype2(1)==-1 .OR. spline_o%ibctype2(2)==-1) iper=1
1074  call genxpkg(spline_o%n2,spline_o%x2(1),spline_o%x2pkg(1,1),&
1075  & iper,imsg,itol,ztol,spline_o%klookup2,ifail)
1076  if(ifail/=0) ier=27
1077  iper=0
1078  if(spline_o%ibctype3(1)==-1 .OR. spline_o%ibctype3(2)==-1) iper=1
1079  call genxpkg(spline_o%n3,spline_o%x3(1),spline_o%x3pkg(1,1),&
1080  & iper,imsg,itol,ztol,spline_o%klookup3,ifail)
1081  if(ifail/=0) ier=27
1082 
1083  spline_o%isReady = 0
1084 
1085  spline_o%fspl(1, 1:in1, 1:in2, 1:in3) = &
1086  & f(1:in1, 1:in2, 1:in3)
1087 
1088  if (spline_o%isHybrid == 1) then
1089 
1090  call mkintrp3d( &
1091  & spline_o%x1(1), spline_o%n1, &
1092  & spline_o%x2(1), spline_o%n2, &
1093  & spline_o%x3(1), spline_o%n3, &
1094  & spline_o%hspline, spline_o%fspl(1,1,1,1), &
1095  & in0,in1,in2,in3, &
1096  & spline_o%ibctype1(1), spline_o%bcval1min(1,1), &
1097  & spline_o%ibctype1(2), spline_o%bcval1max(1,1), &
1098  & spline_o%ibctype2(1), spline_o%bcval2min(1,1), &
1099  & spline_o%ibctype2(2), spline_o%bcval2max(1,1), &
1100  & spline_o%ibctype3(1), spline_o%bcval3min(1,1), &
1101  & spline_o%ibctype3(2), spline_o%bcval3max(1,1), &
1102  & ifail)
1103 
1104  spline_o%ilin1 = 0 ! Hybrid does not compute this
1105  spline_o%ilin2 = 0 ! Hybrid does not compute this
1106  spline_o%ilin3 = 0 ! Hybrid does not compute this
1107 
1108  if(ifail /= 0) then
1109  ier = 98
1110  else
1111  spline_o%isReady = 1
1112  endif
1113 
1114  else if (spline_o%isHermite == 0 .and. spline_o%isLinear == 0) then
1115 
1116  call mktricub( &
1117  & spline_o%x1(1), spline_o%n1, &
1118  & spline_o%x2(1), spline_o%n2, &
1119  & spline_o%x3(1), spline_o%n3, &
1120  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
1121  & spline_o%ibctype1(1), spline_o%bcval1min(1,1), &
1122  & spline_o%ibctype1(2), spline_o%bcval1max(1,1), spline_o%n2, &
1123  & spline_o%ibctype2(1), spline_o%bcval2min(1,1), &
1124  & spline_o%ibctype2(2), spline_o%bcval2max(1,1), spline_o%n1, &
1125  & spline_o%ibctype3(1), spline_o%bcval3min(1,1), &
1126  & spline_o%ibctype3(2), spline_o%bcval3max(1,1), spline_o%n1, &
1127  & spline_o%ilin1, spline_o%ilin2, spline_o%ilin3, ifail)
1128 
1129  if(ifail /= 0) then
1130  ier = 98
1131  else
1132  spline_o%isReady = 1
1133  endif
1134 
1135  else if (spline_o%isLinear == 1) then
1136 
1137  spline_o%ilin1=0
1138  if(in1.gt.2) then
1139  if(spline_o%x1pkg(3,4).eq.0.0_ezspline_r4) spline_o%ilin1=1 ! evenly spaced grid
1140  else
1141  spline_o%ilin1=1
1142  end if
1143 
1144  spline_o%ilin2=0
1145  if(in2.gt.2) then
1146  if(spline_o%x2pkg(3,4).eq.0.0_ezspline_r4) spline_o%ilin2=1 ! evenly spaced grid
1147  else
1148  spline_o%ilin2=1
1149  end if
1150 
1151  spline_o%ilin3=0
1152  if(in3.gt.2) then
1153  if(spline_o%x3pkg(3,4).eq.0.0_ezspline_r4) spline_o%ilin3=1 ! evenly spaced grid
1154  else
1155  spline_o%ilin3=1
1156  end if
1157 
1158  spline_o%isReady = 1 ! no coefficient setup necessary
1159 
1160  else
1161  !
1162  ! Hermite polynomial coefficient setup based on Akima's method
1163  ! (df/dx..etc not required)
1164  !
1165  ipx = 0
1166  if (spline_o%ibctype1(1)==-1 .or. spline_o%ibctype1(2)==-1) then
1167  ipx=1
1168  else if (spline_o%ibctype1(1)<-1 .or. spline_o%ibctype1(1)>1 .or. &
1169  spline_o%ibctype1(2)<-1 .or. spline_o%ibctype1(2)>1 ) then
1170  ipx=99 ! an error...
1171  else if (spline_o%ibctype1(1)==1 .or. spline_o%ibctype1(2)==1) then
1172  ipx=2
1173  do jj=1,spline_o%n3
1174  do ii=1,spline_o%n2
1175  if(spline_o%ibctype1(1)==1) then
1176  spline_o%fspl(2,1,ii,jj)=spline_o%bcval1min(ii,jj)
1177  else
1178  ! hand implemented default BC
1179  df1=(spline_o%fspl(1,2,ii,jj)-spline_o%fspl(1,1,ii,jj))/ &
1180  (spline_o%x1(2)-spline_o%x1(1))
1181  df2=(spline_o%fspl(1,3,ii,jj)-spline_o%fspl(1,2,ii,jj))/ &
1182  (spline_o%x1(3)-spline_o%x1(2))
1183  spline_o%fspl(2,1,ii,jj)=(3*df1-df2)/2
1184  endif
1185  inum=spline_o%n1
1186  if(spline_o%ibctype1(2)==1) then
1187  spline_o%fspl(2,inum,ii,jj)=spline_o%bcval1max(ii,jj)
1188  else
1189  ! hand implemented default BC
1190  df1=(spline_o%fspl(1,inum,ii,jj)-spline_o%fspl(1,inum-1,ii,jj))/ &
1191  (spline_o%x1(inum)-spline_o%x1(inum-1))
1192  df2=(spline_o%fspl(1,inum-1,ii,jj)-spline_o%fspl(1,inum-2,ii,jj))/ &
1193  (spline_o%x1(inum-1)-spline_o%x1(inum-2))
1194  spline_o%fspl(2,inum,ii,jj)=(3*df1-df2)/2
1195  endif
1196  enddo
1197  enddo
1198  endif
1199  ipy = 0
1200  if (spline_o%ibctype2(1)==-1 .or. spline_o%ibctype2(2)==-1) then
1201  ipy=1
1202  else if (spline_o%ibctype2(1)<-1 .or. spline_o%ibctype2(1)>1 .or. &
1203  spline_o%ibctype2(2)<-1 .or. spline_o%ibctype2(2)>1 ) then
1204  ipy=99 ! an error...
1205  else if (spline_o%ibctype2(1)==1 .or. spline_o%ibctype2(2)==1) then
1206  ipy=2
1207  do jj=1,spline_o%n3
1208  do ii=1,spline_o%n1
1209  if(spline_o%ibctype2(1)==1) then
1210  spline_o%fspl(3,ii,1,jj)=spline_o%bcval2min(ii,jj)
1211  else
1212  ! hand implemented default BC
1213  df1=(spline_o%fspl(1,ii,2,jj)-spline_o%fspl(1,ii,1,jj))/ &
1214  (spline_o%x2(2)-spline_o%x2(1))
1215  df2=(spline_o%fspl(1,ii,3,jj)-spline_o%fspl(1,ii,2,jj))/ &
1216  (spline_o%x2(3)-spline_o%x2(2))
1217  spline_o%fspl(3,ii,1,jj)=(3*df1-df2)/2
1218  endif
1219  inum=spline_o%n2
1220  if(spline_o%ibctype2(2)==1) then
1221  spline_o%fspl(3,ii,inum,jj)=spline_o%bcval2max(ii,jj)
1222  else
1223  ! hand implemented default BC
1224  df1=(spline_o%fspl(1,ii,inum,jj)-spline_o%fspl(1,ii,inum-1,jj))/ &
1225  (spline_o%x2(inum)-spline_o%x2(inum-1))
1226  df2=(spline_o%fspl(1,ii,inum-1,jj)-spline_o%fspl(1,ii,inum-2,jj))/ &
1227  (spline_o%x2(inum-1)-spline_o%x2(inum-2))
1228  spline_o%fspl(3,ii,inum,jj)=(3*df1-df2)/2
1229  endif
1230  enddo
1231  enddo
1232  endif
1233  ipz = 0
1234  if (spline_o%ibctype3(1)==-1 .or. spline_o%ibctype3(2)==-1) then
1235  ipz=1
1236  else if (spline_o%ibctype3(1)<-1 .or. spline_o%ibctype3(1)>1 .or. &
1237  spline_o%ibctype3(2)<-1 .or. spline_o%ibctype3(2)>1 ) then
1238  ipz=99 ! an error...
1239  else if (spline_o%ibctype3(1)==1 .or. spline_o%ibctype3(2)==1) then
1240  ipz=2
1241  do jj=1,spline_o%n2
1242  do ii=1,spline_o%n1
1243  if(spline_o%ibctype3(1)==1) then
1244  spline_o%fspl(4,ii,jj,1)=spline_o%bcval3min(ii,jj)
1245  else
1246  ! hand implemented default BC
1247  df1=(spline_o%fspl(1,ii,jj,2)-spline_o%fspl(1,ii,jj,1))/ &
1248  (spline_o%x3(2)-spline_o%x3(1))
1249  df2=(spline_o%fspl(1,ii,jj,3)-spline_o%fspl(1,ii,jj,2))/ &
1250  (spline_o%x3(3)-spline_o%x3(2))
1251  spline_o%fspl(4,ii,jj,1)=(3*df1-df2)/2
1252  endif
1253  inum=spline_o%n3
1254  if(spline_o%ibctype3(2)==1) then
1255  spline_o%fspl(4,ii,jj,inum)=spline_o%bcval3max(ii,jj)
1256  else
1257  ! hand implemented default BC
1258  df1=(spline_o%fspl(1,ii,jj,inum)-spline_o%fspl(1,ii,jj,inum-1))/ &
1259  (spline_o%x3(inum)-spline_o%x3(inum-1))
1260  df2=(spline_o%fspl(1,ii,jj,inum-1)-spline_o%fspl(1,ii,jj,inum-2))/ &
1261  (spline_o%x3(inum-1)-spline_o%x3(inum-2))
1262  spline_o%fspl(4,ii,jj,inum)=(3*df1-df2)/2
1263  endif
1264  enddo
1265  enddo
1266  endif
1267  ifail = 0
1268  call akherm3p(spline_o%x1(1), spline_o%n1, &
1269  & spline_o%x2(1), spline_o%n2, &
1270  & spline_o%x3(1), spline_o%n3, &
1271  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
1272  & spline_o%ilin1, spline_o%ilin2, spline_o%ilin3, &
1273  & ipx,ipy,ipz, ifail)
1274 
1275  if (ifail /=0 ) then
1276  ier = 91
1277  else
1278  spline_o%isReady = 1
1279  endif
1280 
1281  endif
1282 
1283 
1284 end subroutine ezspline_setup3_r4
ezspline_obj::ezspline_allocated
Definition: ezspline_obj.f90:19
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17