V3FIT
ezspline_derivative.f90
1 !/////
2 ! R8 !
3 !/////
4 !
5 ! 1-D
6 !
7 
8 subroutine ezspline_derivative1_r8(spline_o, i1, p1, f, ier)
9  use ezspline_obj
10  implicit none
11  type(EZspline1_r8) spline_o
12  integer, intent(in) :: i1
13  real(ezspline_r8), intent(in) :: p1
14  real(ezspline_r8), intent(out) :: f
15  integer, intent(out) :: ier
16 
17  integer ifail
18 
19  integer ict_herm(2)
20  integer ict(3)
21 
22  ier = 0
23  ifail=0
24  if( .not.ezspline_allocated(spline_o) .or. spline_o%isReady /= 1) then
25  ier = 94
26  return
27  endif
28 
29  f = 0.0_ezspline_r8
30  if(i1 < 0) then
31  ier = 12
32  return
33  endif
34  if(i1 > 3) then
35  ier = 13
36  return
37  endif
38 
39  if(i1 >1 .AND. (spline_o%isHermite==1 .or. spline_o%isLinear==1)) then
40  ier = 24; if(spline_o%isLinear==1) ier = 46
41  return
42  endif
43 
44  if(spline_o%isHermite==1 .or. spline_o%isLinear==1) then
45 
46  if(i1 == 1) then
47  ict_herm = (/0, 1/)
48  else
49  ict_herm = (/1, 0/)
50  endif
51 
52  if(spline_o%isLinear == 0) then
53  call r8herm1ev(p1, &
54  & spline_o%x1(1), spline_o%n1, &
55  & spline_o%ilin1, &
56  & spline_o%fspl(1,1), &
57  & ict_herm, f, ifail)
58  else
59  call r8pc1ev(p1, &
60  & spline_o%x1(1), spline_o%n1, &
61  & spline_o%ilin1, &
62  & spline_o%fspl(1,1), &
63  & ict_herm, f, ifail)
64  endif
65 
66  else
67 
68  call ezmake_ict1(i1,ict)
69 
70  call r8evspline(p1, &
71  & spline_o%x1(1), spline_o%n1, &
72  & spline_o%ilin1, &
73  & spline_o%fspl(1,1), &
74  & ict, f, ifail)
75 
76  endif
77  if(ifail /= 0) ier = 96
78 
79 end subroutine ezspline_derivative1_r8
80 
81 subroutine ezspline_derivative1_array_r8(spline_o, i1, k1, p1, f, ier)
82  use ezspline_obj
83  implicit none
84  type(EZspline1_r8) spline_o
85  integer, intent(in) :: i1
86  integer, intent(in) :: k1
87  real(ezspline_r8), intent(in) :: p1(k1)
88  real(ezspline_r8), intent(out) :: f(k1)
89  integer, intent(out) :: ier
90 
91  integer ifail
92 
93  integer ict(3)
94  integer ict_herm(2)
95  integer:: iwarn=0
96 
97  ier = 0
98  ifail=0
99  if( .not.ezspline_allocated(spline_o) .or. spline_o%isReady /= 1) then
100  ier = 94
101  return
102  endif
103 
104  f = 0.0_ezspline_r8
105  if(i1 < 0) then
106  ier = 12
107  return
108  endif
109  if(i1 > 3) then
110  ier = 13
111  return
112  endif
113 
114  if(i1 >1 .AND. (spline_o%isHermite==1 .or. spline_o%isLinear==1)) then
115  ier = 24; if(spline_o%isLinear==1) ier = 46
116  return
117  endif
118 
119  if(spline_o%isHermite==1 .or. spline_o%isLinear==1) then
120 
121  if(i1 == 1) then
122  ict_herm = (/0, 1/)
123  else
124  ict_herm = (/1, 0/)
125  endif
126 
127  if(spline_o%isLinear == 0) then
128  call r8vecherm1(ict_herm, k1, p1, k1, f, &
129  & spline_o%n1, spline_o%x1pkg(1,1), &
130  & spline_o%fspl(1,1), &
131  & iwarn, ifail)
132  else
133  call r8vecpc1(ict_herm, k1, p1, k1, f, &
134  & spline_o%n1, spline_o%x1pkg(1,1), &
135  & spline_o%fspl(1,1), &
136  & iwarn, ifail)
137  endif
138 
139  else
140 
141  call ezmake_ict1(i1,ict)
142 
143  call r8vecspline(ict, k1, p1, k1, f, &
144  & spline_o%n1, spline_o%x1pkg(1,1), &
145  & spline_o%fspl(1,1), &
146  & iwarn, ifail)
147 
148  endif
149 
150  if(ifail /= 0) ier = 96
151 
152 end subroutine ezspline_derivative1_array_r8
153 
154 
155 !!
156 !! 2-D
157 !!
158 
159 
160 subroutine ezspline_derivative2_r8(spline_o, i1, i2, p1, p2, f, ier)
161  use ezspline_obj
162  implicit none
163  type(EZspline2_r8) spline_o
164  integer, intent(in) :: i1, i2
165  real(ezspline_r8), intent(in) :: p1, p2
166  real(ezspline_r8), intent(out) :: f
167  integer, intent(out) :: ier
168 
169  integer ifail
170 
171  integer ict(6)
172  integer ict_herm(4)
173 
174  ier = 0
175  ifail=0
176  if( .not.ezspline_allocated(spline_o) .or. spline_o%isReady /= 1) then
177  ier = 94
178  return
179  endif
180 
181  if(i1<0 .OR. i2<0) then
182  ier = 12
183  return
184  endif
185  if(max(i1,i2) > 3) then
186  ier = 13
187  return
188  endif
189 
190  if(max(i1,i2) > 1 .AND. (spline_o%isHermite==1 .or. spline_o%isLinear==1)) then
191  ier = 24; if(spline_o%isLinear==1) ier = 46
192  return
193  endif
194 
195  if(spline_o%isHybrid == 1) then
196 
197  call ezmake_ict2(i1,i2,ict)
198 
199  call r8evintrp2d(p1, p2, &
200  & spline_o%x1(1), spline_o%n1, &
201  & spline_o%x2(1), spline_o%n2, &
202  & spline_o%hspline, spline_o%fspl(1,1,1), &
203  & size(spline_o%fspl,1), size(spline_o%fspl,2), &
204  & size(spline_o%fspl,3), &
205  & ict, f, ifail)
206 
207  else if(spline_o%isHermite==1 .or. spline_o%isLinear==1) then
208  if(i1 == 1 .AND. i2 == 0) then
209  ict_herm = (/0, 1, 0, 0/)
210  else if(i1 == 0 .AND. i2 == 1) then
211  ict_herm = (/0, 0, 1, 0/)
212  else if(i1 == 1 .AND. i2 == 1) then
213  ict_herm = (/0, 0, 0, 1/)
214  else
215  ict_herm = (/1, 0, 0, 0/)
216  endif
217 
218  if(spline_o%isLinear == 0) then
219  call r8herm2ev(p1, p2, &
220  & spline_o%x1(1), spline_o%n1, &
221  & spline_o%x2(1), spline_o%n2, &
222  & spline_o%ilin1, spline_o%ilin2, &
223  & spline_o%fspl(1,1,1), spline_o%n1, &
224  & ict_herm, f, ifail)
225  else
226  call r8pc2ev(p1, p2, &
227  & spline_o%x1(1), spline_o%n1, &
228  & spline_o%x2(1), spline_o%n2, &
229  & spline_o%ilin1, spline_o%ilin2, &
230  & spline_o%fspl(1,1,1), spline_o%n1, &
231  & ict_herm, f, ifail)
232  endif
233 
234  else
235 
236  call ezmake_ict2(i1,i2,ict)
237 
238  call r8evbicub(p1, p2, &
239  & spline_o%x1(1), spline_o%n1, &
240  & spline_o%x2(1), spline_o%n2, &
241  & spline_o%ilin1, spline_o%ilin2, &
242  & spline_o%fspl(1,1,1), spline_o%n1, &
243  & ict, f, ifail)
244 
245  endif
246  if(ifail /= 0) ier = 96
247 
248 end subroutine ezspline_derivative2_r8
249 
250 
251 
252 subroutine ezspline_derivative2_array_r8(spline_o, i1, i2, &
253  & k1, k2, p1, p2, f, ier)
254  use ezspline_obj
255  implicit none
256  type(EZspline2_r8) spline_o
257  integer, intent(in) :: i1, i2, k1, k2
258  real(ezspline_r8), intent(in) :: p1(k1), p2(k2)
259  real(ezspline_r8), intent(out) :: f(k1,k2)
260  integer, intent(out) :: ier
261 
262  integer ifail
263 
264  integer ict(6)
265  integer ict_herm(4)
266 
267  real(ezspline_r8), dimension(:), allocatable :: p1_cloud, p2_cloud
268  integer k12
269  integer:: iwarn=0
270 
271  ier = 0
272  ifail=0
273  if( .not.ezspline_allocated(spline_o) .or. spline_o%isReady /= 1) then
274  ier = 94
275  return
276  endif
277 
278  f = 0.0_ezspline_r8
279  if(i1<0 .OR. i2<0) then
280  ier = 12
281  return
282  endif
283  if(max(i1,i2) > 3) then
284  ier = 13
285  return
286  endif
287 
288  if(max(i1,i2) > 1 .AND. (spline_o%isHermite==1 .or. spline_o%isLinear==1)) then
289  ier = 24; if(spline_o%isLinear==1) ier = 46
290  return
291  endif
292 
293  k12 = k1*k2
294  allocate(p1_cloud(k12), p2_cloud(k12), stat=ifail)
295  if(ifail /= 0) then
296  ier = 32
297  return
298  endif
299 
300  p1_cloud = reshape( &
301  & source=spread(source=p1, dim=2, ncopies=k2), &
302  & shape=(/k12/))
303  p2_cloud = reshape( &
304  & source=spread(source=p2, dim=1, ncopies=k1), &
305  & shape=(/k12/))
306 
307  if(spline_o%isHybrid == 1) then
308 
309  call ezmake_ict2(i1,i2,ict)
310 
311  call r8vecintrp2d(ict, k12, p1_cloud, p2_cloud, k12, f, &
312  & spline_o%n1, spline_o%x1pkg(1,1), &
313  & spline_o%n2, spline_o%x2pkg(1,1), &
314  & spline_o%hspline, spline_o%fspl(1,1,1), &
315  & size(spline_o%fspl,1), size(spline_o%fspl,2), &
316  & size(spline_o%fspl,3), &
317  & iwarn, ifail)
318 
319  else if(spline_o%isHermite==1 .or. spline_o%isLinear==1)then
320  if(i1 == 1 .AND. i2 == 0) then
321  ict_herm = (/0, 1, 0, 0/)
322  else if(i1 == 0 .AND. i2 == 1) then
323  ict_herm = (/0, 0, 1, 0/)
324  else if(i1 == 1 .AND. i2 == 1) then
325  ict_herm = (/0, 0, 0, 1/)
326  else
327  ict_herm = (/1, 0, 0, 0/)
328  endif
329 
330  if(spline_o%isLinear == 0) then
331  call r8vecherm2(ict_herm, k12, p1_cloud, p2_cloud, k12, f, &
332  & spline_o%n1, spline_o%x1pkg(1,1), &
333  & spline_o%n2, spline_o%x2pkg(1,1), &
334  & spline_o%fspl(1,1,1), spline_o%n1, &
335  & iwarn, ifail)
336  else
337  call r8vecpc2(ict_herm, k12, p1_cloud, p2_cloud, k12, f, &
338  & spline_o%n1, spline_o%x1pkg(1,1), &
339  & spline_o%n2, spline_o%x2pkg(1,1), &
340  & spline_o%fspl(1,1,1), spline_o%n1, &
341  & iwarn, ifail)
342  endif
343 
344  else
345 
346  call ezmake_ict2(i1,i2,ict) ! higher derivatives
347 
348  call r8vecbicub(ict, k12, p1_cloud, p2_cloud, k12, f, &
349  & spline_o%n1,spline_o%x1pkg(1,1), &
350  & spline_o%n2,spline_o%x2pkg(1,1), &
351  & spline_o%fspl(1,1,1), spline_o%n1, &
352  & iwarn, ifail)
353 
354  endif
355 
356 
357  if(ifail /= 0) ier = 96
358 
359  deallocate(p1_cloud, p2_cloud, stat=ifail)
360  if(ifail /= 0) then
361  ier = 33
362  return
363  endif
364 
365 
366 end subroutine ezspline_derivative2_array_r8
367 
368 
369 subroutine ezspline_derivative2_cloud_r8(spline_o, i1, i2, &
370  & k, p1, p2, f, ier)
371  use ezspline_obj
372  implicit none
373  type(EZspline2_r8) spline_o
374  integer, intent(in) :: i1, i2, k
375  real(ezspline_r8), intent(in) :: p1(k), p2(k)
376  real(ezspline_r8), intent(out) :: f(k)
377  integer, intent(out) :: ier
378 
379  integer ifail
380 
381  integer ict(6)
382  integer ict_herm(4)
383  integer:: iwarn=0
384 
385  ier = 0
386  ifail=0
387  if( .not.ezspline_allocated(spline_o) .or. spline_o%isReady /= 1) then
388  ier = 94
389  return
390  endif
391 
392  f = 0.0_ezspline_r8
393  if(i1<0 .OR. i2<0) then
394  ier = 12
395  return
396  endif
397  if(max(i1,i2) > 3) then
398  ier = 13
399  return
400  endif
401 
402  if(max(i1,i2) > 1 .AND. (spline_o%isHermite==1 .or. spline_o%isLinear==1)) then
403  ier = 24; if(spline_o%isLinear==1) ier = 46
404  return
405  endif
406 
407  if(spline_o%isHybrid == 1) then
408 
409  call ezmake_ict2(i1,i2,ict)
410 
411  call r8vecintrp2d(ict, k, p1, p2, k, f, &
412  & spline_o%n1, spline_o%x1pkg(1,1), &
413  & spline_o%n2, spline_o%x2pkg(1,1), &
414  & spline_o%hspline, spline_o%fspl(1,1,1), &
415  & size(spline_o%fspl,1), size(spline_o%fspl,2), &
416  & size(spline_o%fspl,3), &
417  & iwarn, ifail)
418 
419  else if(spline_o%isHermite==1 .or. spline_o%isLinear==1)then
420  if(i1 == 1 .AND. i2 == 0) then
421  ict_herm = (/0, 1, 0, 0/)
422  else if(i1 == 0 .AND. i2 == 1) then
423  ict_herm = (/0, 0, 1, 0/)
424  else if(i1 == 1 .AND. i2 == 1) then
425  ict_herm = (/0, 0, 0, 1/)
426  else
427  ict_herm = (/1, 0, 0, 0/)
428  endif
429 
430  if(spline_o%isLinear == 0) then
431  call r8vecherm2(ict_herm, k, p1, p2, &
432  & k, f, &
433  & spline_o%n1, spline_o%x1pkg(1,1), &
434  & spline_o%n2, spline_o%x2pkg(1,1), &
435  & spline_o%fspl(1,1,1), spline_o%n1, &
436  & iwarn, ifail)
437  else
438  call r8vecpc2(ict_herm, k, p1, p2, &
439  & k, f, &
440  & spline_o%n1, spline_o%x1pkg(1,1), &
441  & spline_o%n2, spline_o%x2pkg(1,1), &
442  & spline_o%fspl(1,1,1), spline_o%n1, &
443  & iwarn, ifail)
444  endif
445 
446  else
447 
448  call ezmake_ict2(i1,i2,ict)
449 
450  call r8vecbicub(ict, k, p1, p2, &
451  & k, f, &
452  & spline_o%n1, spline_o%x1pkg(1,1), &
453  & spline_o%n2, spline_o%x2pkg(1,1), &
454  & spline_o%fspl(1,1,1), spline_o%n1, &
455  & iwarn, ifail)
456 
457  endif
458 
459  if(ifail /= 0) ier = 96
460 
461 end subroutine ezspline_derivative2_cloud_r8
462 
463 !!!
464 !!! 3-D
465 !!!
466 
467 
468 subroutine ezspline_derivative3_r8(spline_o, i1, i2, i3, p1, p2, p3, f, ier)
469  use ezspline_obj
470  implicit none
471  type(EZspline3_r8) spline_o
472  integer, intent(in) :: i1, i2, i3
473  real(ezspline_r8), intent(in) :: p1, p2, p3
474  real(ezspline_r8), intent(out) :: f
475  integer, intent(out) :: ier
476 
477  integer ifail
478 
479  integer ict(10)
480  integer ict_herm(8)
481 
482  ier = 0
483  ifail=0
484  if( .not.ezspline_allocated(spline_o) .or. spline_o%isReady /= 1) then
485  ier = 94
486  return
487  endif
488 
489  f = 0.0_ezspline_r8
490  if(i1<0 .OR. i2<0 .OR. i3<0) then
491  ier = 12
492  return
493  endif
494  if(max(i1,max(i2,i3)) > 3) then
495  ier = 13
496  return
497  endif
498 
499  if(maxval((/i1,i2,i3/)) > 1 .AND. (spline_o%isHermite==1 .or. spline_o%isLinear==1)) then
500  ier = 24; if(spline_o%isLinear==1) ier = 46
501  return
502  endif
503 
504  if(spline_o%isHybrid == 1) then
505 
506  call ezmake_ict3(i1,i2,i3,ict)
507 
508  call r8evintrp3d(p1, p2, p3, &
509  & spline_o%x1(1), spline_o%n1, &
510  & spline_o%x2(1), spline_o%n2, &
511  & spline_o%x3(1), spline_o%n3, &
512  & spline_o%hspline, spline_o%fspl(1,1,1,1), &
513  & size(spline_o%fspl,1), size(spline_o%fspl,2), &
514  & size(spline_o%fspl,3), size(spline_o%fspl,4), &
515  & ict, f, ifail)
516 
517  else if(spline_o%isHermite==1 .or. spline_o%isLinear==1)then
518  if(i1 == 1 .AND. i2 == 0 .AND. i3 == 0) then
519  ict_herm = (/0, 1, 0, 0, 0, 0, 0, 0/)
520  else if(i1 == 0 .AND. i2 == 1 .AND. i3 == 0) then
521  ict_herm = (/0, 0, 1, 0, 0, 0, 0, 0/)
522  else if(i1 == 0 .AND. i2 == 0 .AND. i3 == 1) then
523  ict_herm = (/0, 0, 0, 1, 0, 0, 0, 0/)
524  else if(i1 == 1 .AND. i2 == 1 .AND. i3 == 0) then
525  ict_herm = (/0, 0, 0, 0, 1, 0, 0, 0/)
526  else if(i1 == 1 .AND. i2 == 0 .AND. i3 == 1) then
527  ict_herm = (/0, 0, 0, 0, 0, 1, 0, 0/)
528  else if(i1 == 0 .AND. i2 == 1 .AND. i3 == 1) then
529  ict_herm = (/0, 0, 0, 0, 0, 0, 1, 0/)
530  else if(i1 == 1 .AND. i2 == 1 .AND. i3 == 1) then
531  ict_herm = (/0, 0, 0, 0, 0, 0, 0, 1/)
532  else
533  ict_herm = (/1, 0, 0, 0, 0, 0, 0, 0/)
534  endif
535 
536  if(spline_o%isLinear == 0) then
537  call r8herm3ev(p1, p2, p3, &
538  & spline_o%x1(1), spline_o%n1, &
539  & spline_o%x2(1), spline_o%n2, &
540  & spline_o%x3(1), spline_o%n3, &
541  & spline_o%ilin1, spline_o%ilin2, spline_o%ilin3, &
542  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
543  & ict_herm, f, ifail)
544  else
545  call r8pc3ev(p1, p2, p3, &
546  & spline_o%x1(1), spline_o%n1, &
547  & spline_o%x2(1), spline_o%n2, &
548  & spline_o%x3(1), spline_o%n3, &
549  & spline_o%ilin1, spline_o%ilin2, spline_o%ilin3, &
550  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
551  & ict_herm, f, ifail)
552  endif
553 
554  else
555 
556  call ezmake_ict3(i1,i2,i3,ict)
557 
558  call r8evtricub(p1, p2, p3, &
559  & spline_o%x1(1), spline_o%n1, &
560  & spline_o%x2(1), spline_o%n2, &
561  & spline_o%x3(1), spline_o%n3, &
562  & spline_o%ilin1, spline_o%ilin2, spline_o%ilin3, &
563  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
564  & ict, f, ifail)
565 
566  endif
567  if(ifail /= 0) ier = 96
568 
569 end subroutine ezspline_derivative3_r8
570 
571 
572 subroutine ezspline_derivative3_array_r8(spline_o, i1, i2, i3, &
573  & k1, k2, k3, p1, p2, p3, f, ier)
574  use ezspline_obj
575  implicit none
576  type(EZspline3_r8) spline_o
577  integer, intent(in) :: i1, i2, i3, k1, k2, k3
578  real(ezspline_r8), intent(in) :: p1(k1), p2(k2), p3(k3)
579  real(ezspline_r8), intent(out) :: f(k1,k2,k3)
580  integer, intent(out) :: ier
581 
582  integer ifail
583 
584  integer ict(10)
585  integer ict_herm(8)
586 
587  real(ezspline_r8), dimension(:), allocatable :: p1_cloud, p2_cloud, p3_cloud
588  integer k123
589  integer:: iwarn=0
590 
591  ier = 0
592  ifail=0
593  if( .not.ezspline_allocated(spline_o) .or. spline_o%isReady /= 1) then
594  ier = 94
595  return
596  endif
597 
598  f = 0.0_ezspline_r8
599  if(i1<0 .OR. i2<0 .OR. i3<0) then
600  ier = 12
601  return
602  endif
603  if(max(i1,max(i2,i3)) > 3) then
604  ier = 13
605  return
606  endif
607 
608  if(maxval((/i1,i2,i3/)) > 1 .AND. (spline_o%isHermite==1 .or. spline_o%isLinear==1)) then
609  ier = 24; if(spline_o%isLinear==1) ier = 46
610  return
611  endif
612 
613  k123 = k1*k2*k3
614  allocate(p1_cloud(k123), p2_cloud(k123), p3_cloud(k123), stat=ifail)
615  if(ifail /= 0) then
616  ier = 32
617  return
618  endif
619 
620  p1_cloud = reshape(source=spread( &
621  & source=spread(source=p1, dim=2, ncopies=k2), &
622  & dim=3, ncopies=k3), shape=(/k123/))
623  p2_cloud = reshape(source=spread( &
624  & source=spread(source=p2, dim=1, ncopies=k1), &
625  & dim=3, ncopies=k3), shape=(/k123/))
626  p3_cloud = reshape(source=spread( &
627  & source=spread(source=p3, dim=1, ncopies=k1), &
628  & dim=2, ncopies=k2), shape=(/k123/))
629 
630  if(spline_o%isHybrid == 1) then
631 
632  call ezmake_ict3(i1,i2,i3,ict)
633 
634  call r8vecintrp3d(ict, k123, p1_cloud, p2_cloud, p3_cloud, k123, f, &
635  & spline_o%n1, spline_o%x1pkg(1,1), &
636  & spline_o%n2, spline_o%x2pkg(1,1), &
637  & spline_o%n3, spline_o%x3pkg(1,1), &
638  & spline_o%hspline, spline_o%fspl(1,1,1,1), &
639  & size(spline_o%fspl,1), size(spline_o%fspl,2), &
640  & size(spline_o%fspl,3), size(spline_o%fspl,4), &
641  & iwarn, ifail)
642 
643  else if(spline_o%isHermite==1 .or. spline_o%isLinear==1)then
644  if(i1 == 1 .AND. i2 == 0 .AND. i3 == 0) then
645  ict_herm = (/0, 1, 0, 0, 0, 0, 0, 0/)
646  else if(i1 == 0 .AND. i2 == 1 .AND. i3 == 0) then
647  ict_herm = (/0, 0, 1, 0, 0, 0, 0, 0/)
648  else if(i1 == 0 .AND. i2 == 0 .AND. i3 == 1) then
649  ict_herm = (/0, 0, 0, 1, 0, 0, 0, 0/)
650  else if(i1 == 1 .AND. i2 == 1 .AND. i3 == 0) then
651  ict_herm = (/0, 0, 0, 0, 1, 0, 0, 0/)
652  else if(i1 == 1 .AND. i2 == 0 .AND. i3 == 1) then
653  ict_herm = (/0, 0, 0, 0, 0, 1, 0, 0/)
654  else if(i1 == 0 .AND. i2 == 1 .AND. i3 == 1) then
655  ict_herm = (/0, 0, 0, 0, 0, 0, 1, 0/)
656  else if(i1 == 1 .AND. i2 == 1 .AND. i3 == 1) then
657  ict_herm = (/0, 0, 0, 0, 0, 0, 0, 1/)
658  else
659  ict_herm = (/1, 0, 0, 0, 0, 0, 0, 0/)
660  endif
661 
662  if(spline_o%isLinear == 0) then
663  call r8vecherm3(ict_herm, k123, p1_cloud, p2_cloud, p3_cloud, k123, &
664  & f, &
665  & spline_o%n1, spline_o%x1pkg(1,1), &
666  & spline_o%n2, spline_o%x2pkg(1,1), &
667  & spline_o%n3, spline_o%x3pkg(1,1), &
668  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
669  & iwarn, ifail)
670  else
671  call r8vecpc3(ict_herm, k123, p1_cloud, p2_cloud, p3_cloud, k123, &
672  & f, &
673  & spline_o%n1, spline_o%x1pkg(1,1), &
674  & spline_o%n2, spline_o%x2pkg(1,1), &
675  & spline_o%n3, spline_o%x3pkg(1,1), &
676  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
677  & iwarn, ifail)
678  endif
679 
680  else
681 
682  call ezmake_ict3(i1,i2,i3,ict)
683 
684  call r8vectricub(ict, k123, p1_cloud, p2_cloud, p3_cloud, k123, f, &
685  & spline_o%n1,spline_o%x1pkg(1,1), &
686  & spline_o%n2,spline_o%x2pkg(1,1), &
687  & spline_o%n3,spline_o%x3pkg(1,1), &
688  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
689  & iwarn, ifail)
690 
691  endif
692 
693 
694  if(ifail /= 0) ier = 96
695 
696  deallocate(p1_cloud, p2_cloud, p3_cloud, stat=ifail)
697  if(ifail /= 0) then
698  ier = 33
699  return
700  endif
701 
702 
703 end subroutine ezspline_derivative3_array_r8
704 
705 subroutine ezspline_derivative3_cloud_r8(spline_o, i1, i2, i3, &
706  & k, p1, p2, p3, f, ier)
707  use ezspline_obj
708  implicit none
709  type(EZspline3_r8) spline_o
710  integer, intent(in) :: i1, i2, i3, k
711  real(ezspline_r8), intent(in) :: p1(k), p2(k), p3(k)
712  real(ezspline_r8), intent(out) :: f(k)
713  integer, intent(out) :: ier
714 
715  integer ifail
716 
717  integer ict(10)
718  integer ict_herm(8)
719  integer:: iwarn=0
720 
721  ier = 0
722  ifail=0
723  if( .not.ezspline_allocated(spline_o) .or. spline_o%isReady /= 1) then
724  ier = 94
725  return
726  endif
727 
728  f = 0.0_ezspline_r8
729  if(i1<0 .OR. i2<0 .OR. i3<0) then
730  ier = 12
731  return
732  endif
733  if(max(i1,max(i2,i3)) > 3) then
734  ier = 13
735  return
736  endif
737 
738  if(maxval((/i1,i2,i3/)) > 1 .AND. (spline_o%isHermite==1 .or. spline_o%isLinear==1)) then
739  ier = 24; if(spline_o%isLinear==1) ier = 46
740  return
741  endif
742 
743  if(spline_o%isHybrid == 1) then
744 
745  call ezmake_ict3(i1,i2,i3,ict)
746 
747  call r8vecintrp3d(ict, k, p1, p2, p3, k, f, &
748  & spline_o%n1, spline_o%x1pkg(1,1), &
749  & spline_o%n2, spline_o%x2pkg(1,1), &
750  & spline_o%n3, spline_o%x3pkg(1,1), &
751  & spline_o%hspline, spline_o%fspl(1,1,1,1), &
752  & size(spline_o%fspl,1), size(spline_o%fspl,2), &
753  & size(spline_o%fspl,3), size(spline_o%fspl,4), &
754  & iwarn, ifail)
755 
756  else if(spline_o%isHermite==1 .or. spline_o%isLinear==1)then
757  if(i1 == 1 .AND. i2 == 0 .AND. i3 == 0) then
758  ict_herm = (/0, 1, 0, 0, 0, 0, 0, 0/)
759  else if(i1 == 0 .AND. i2 == 1 .AND. i3 == 0) then
760  ict_herm = (/0, 0, 1, 0, 0, 0, 0, 0/)
761  else if(i1 == 0 .AND. i2 == 0 .AND. i3 == 1) then
762  ict_herm = (/0, 0, 0, 1, 0, 0, 0, 0/)
763  else if(i1 == 1 .AND. i2 == 1 .AND. i3 == 0) then
764  ict_herm = (/0, 0, 0, 0, 1, 0, 0, 0/)
765  else if(i1 == 1 .AND. i2 == 0 .AND. i3 == 1) then
766  ict_herm = (/0, 0, 0, 0, 0, 1, 0, 0/)
767  else if(i1 == 0 .AND. i2 == 1 .AND. i3 == 1) then
768  ict_herm = (/0, 0, 0, 0, 0, 0, 1, 0/)
769  else if(i1 == 1 .AND. i2 == 1 .AND. i3 == 1) then
770  ict_herm = (/0, 0, 0, 0, 0, 0, 0, 1/)
771  else
772  ict_herm = (/1, 0, 0, 0, 0, 0, 0, 0/)
773  endif
774 
775  if(spline_o%isLinear == 0) then
776  call r8vecherm3(ict_herm, k, p1, p2, p3, &
777  & k, f, &
778  & spline_o%n1, spline_o%x1pkg(1,1), &
779  & spline_o%n2, spline_o%x2pkg(1,1), &
780  & spline_o%n3, spline_o%x3pkg(1,1), &
781  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
782  & iwarn, ifail)
783  else
784  call r8vecpc3(ict_herm, k, p1, p2, p3, &
785  & k, f, &
786  & spline_o%n1, spline_o%x1pkg(1,1), &
787  & spline_o%n2, spline_o%x2pkg(1,1), &
788  & spline_o%n3, spline_o%x3pkg(1,1), &
789  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
790  & iwarn, ifail)
791  endif
792 
793  else
794 
795  call ezmake_ict3(i1,i2,i3,ict)
796 
797  call r8vectricub(ict, k, p1, p2, p3, &
798  & k, f, &
799  & spline_o%n1, spline_o%x1pkg(1,1), &
800  & spline_o%n2, spline_o%x2pkg(1,1), &
801  & spline_o%n3, spline_o%x3pkg(1,1), &
802  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
803  & iwarn, ifail)
804 
805  endif
806 
807  if(ifail /= 0) ier = 96
808 
809 end subroutine ezspline_derivative3_cloud_r8
810 !/////
811 ! R4 !
812 !/////
813 !
814 ! 1-D
815 !
816 
817 subroutine ezspline_derivative1_r4(spline_o, i1, p1, f, ier)
818  use ezspline_obj
819  implicit none
820  type(EZspline1_r4) spline_o
821  integer, intent(in) :: i1
822  real(ezspline_r4), intent(in) :: p1
823  real(ezspline_r4), intent(out) :: f
824  integer, intent(out) :: ier
825 
826  integer ifail
827 
828  integer ict_herm(2)
829  integer ict(3)
830 
831  ier = 0
832  ifail=0
833  if( .not.ezspline_allocated(spline_o) .or. spline_o%isReady /= 1) then
834  ier = 94
835  return
836  endif
837 
838  f = 0.0_ezspline_r4
839  if(i1 < 0) then
840  ier = 12
841  return
842  endif
843  if(i1 > 3) then
844  ier = 13
845  return
846  endif
847 
848  if(i1 >1 .AND. (spline_o%isHermite==1 .or. spline_o%isLinear==1)) then
849  ier = 24; if(spline_o%isLinear==1) ier = 46
850  return
851  endif
852 
853  if(spline_o%isHermite==1 .or. spline_o%isLinear==1) then
854 
855  if(i1 == 1) then
856  ict_herm = (/0, 1/)
857  else
858  ict_herm = (/1, 0/)
859  endif
860 
861  if(spline_o%isLinear == 0) then
862  call herm1ev(p1, &
863  & spline_o%x1(1), spline_o%n1, &
864  & spline_o%ilin1, &
865  & spline_o%fspl(1,1), &
866  & ict_herm, f, ifail)
867  else
868  call pc1ev(p1, &
869  & spline_o%x1(1), spline_o%n1, &
870  & spline_o%ilin1, &
871  & spline_o%fspl(1,1), &
872  & ict_herm, f, ifail)
873  endif
874 
875  else
876 
877  call ezmake_ict1(i1,ict)
878 
879  call evspline(p1, &
880  & spline_o%x1(1), spline_o%n1, &
881  & spline_o%ilin1, &
882  & spline_o%fspl(1,1), &
883  & ict, f, ifail)
884 
885  endif
886  if(ifail /= 0) ier = 96
887 
888 end subroutine ezspline_derivative1_r4
889 
890 subroutine ezspline_derivative1_array_r4(spline_o, i1, k1, p1, f, ier)
891  use ezspline_obj
892  implicit none
893  type(EZspline1_r4) spline_o
894  integer, intent(in) :: i1
895  integer, intent(in) :: k1
896  real(ezspline_r4), intent(in) :: p1(k1)
897  real(ezspline_r4), intent(out) :: f(k1)
898  integer, intent(out) :: ier
899 
900  integer ifail
901 
902  integer ict(3)
903  integer ict_herm(2)
904  integer:: iwarn=0
905 
906  ier = 0
907  ifail=0
908  if( .not.ezspline_allocated(spline_o) .or. spline_o%isReady /= 1) then
909  ier = 94
910  return
911  endif
912 
913  f = 0.0_ezspline_r4
914  if(i1 < 0) then
915  ier = 12
916  return
917  endif
918  if(i1 > 3) then
919  ier = 13
920  return
921  endif
922 
923  if(i1 >1 .AND. (spline_o%isHermite==1 .or. spline_o%isLinear==1)) then
924  ier = 24; if(spline_o%isLinear==1) ier = 46
925  return
926  endif
927 
928  if(spline_o%isHermite==1 .or. spline_o%isLinear==1) then
929 
930  if(i1 == 1) then
931  ict_herm = (/0, 1/)
932  else
933  ict_herm = (/1, 0/)
934  endif
935 
936  if(spline_o%isLinear == 0) then
937  call vecherm1(ict_herm, k1, p1, k1, f, &
938  & spline_o%n1, spline_o%x1pkg(1,1), &
939  & spline_o%fspl(1,1), &
940  & iwarn, ifail)
941  else
942  call vecpc1(ict_herm, k1, p1, k1, f, &
943  & spline_o%n1, spline_o%x1pkg(1,1), &
944  & spline_o%fspl(1,1), &
945  & iwarn, ifail)
946  endif
947 
948  else
949 
950  call ezmake_ict1(i1,ict)
951 
952  call vecspline(ict, k1, p1, k1, f, &
953  & spline_o%n1, spline_o%x1pkg(1,1), &
954  & spline_o%fspl(1,1), &
955  & iwarn, ifail)
956 
957  endif
958 
959  if(ifail /= 0) ier = 96
960 
961 end subroutine ezspline_derivative1_array_r4
962 
963 
964 !!
965 !! 2-D
966 !!
967 
968 
969 subroutine ezspline_derivative2_r4(spline_o, i1, i2, p1, p2, f, ier)
970  use ezspline_obj
971  implicit none
972  type(EZspline2_r4) spline_o
973  integer, intent(in) :: i1, i2
974  real(ezspline_r4), intent(in) :: p1, p2
975  real(ezspline_r4), intent(out) :: f
976  integer, intent(out) :: ier
977 
978  integer ifail
979 
980  integer ict(6)
981  integer ict_herm(4)
982 
983  ier = 0
984  ifail=0
985  if( .not.ezspline_allocated(spline_o) .or. spline_o%isReady /= 1) then
986  ier = 94
987  return
988  endif
989 
990  if(i1<0 .OR. i2<0) then
991  ier = 12
992  return
993  endif
994  if(max(i1,i2) > 3) then
995  ier = 13
996  return
997  endif
998 
999  if(max(i1,i2) > 1 .AND. (spline_o%isHermite==1 .or. spline_o%isLinear==1)) then
1000  ier = 24; if(spline_o%isLinear==1) ier = 46
1001  return
1002  endif
1003 
1004  if(spline_o%isHybrid == 1) then
1005 
1006  call ezmake_ict2(i1,i2,ict)
1007 
1008  call evintrp2d(p1, p2, &
1009  & spline_o%x1(1), spline_o%n1, &
1010  & spline_o%x2(1), spline_o%n2, &
1011  & spline_o%hspline, spline_o%fspl(1,1,1), &
1012  & size(spline_o%fspl,1), size(spline_o%fspl,2), &
1013  & size(spline_o%fspl,3), &
1014  & ict, f, ifail)
1015 
1016  else if(spline_o%isHermite==1 .or. spline_o%isLinear==1) then
1017  if(i1 == 1 .AND. i2 == 0) then
1018  ict_herm = (/0, 1, 0, 0/)
1019  else if(i1 == 0 .AND. i2 == 1) then
1020  ict_herm = (/0, 0, 1, 0/)
1021  else if(i1 == 1 .AND. i2 == 1) then
1022  ict_herm = (/0, 0, 0, 1/)
1023  else
1024  ict_herm = (/1, 0, 0, 0/)
1025  endif
1026 
1027  if(spline_o%isLinear == 0) then
1028  call herm2ev(p1, p2, &
1029  & spline_o%x1(1), spline_o%n1, &
1030  & spline_o%x2(1), spline_o%n2, &
1031  & spline_o%ilin1, spline_o%ilin2, &
1032  & spline_o%fspl(1,1,1), spline_o%n1, &
1033  & ict_herm, f, ifail)
1034  else
1035  call pc2ev(p1, p2, &
1036  & spline_o%x1(1), spline_o%n1, &
1037  & spline_o%x2(1), spline_o%n2, &
1038  & spline_o%ilin1, spline_o%ilin2, &
1039  & spline_o%fspl(1,1,1), spline_o%n1, &
1040  & ict_herm, f, ifail)
1041  endif
1042 
1043  else
1044 
1045  call ezmake_ict2(i1,i2,ict)
1046 
1047  call evbicub(p1, p2, &
1048  & spline_o%x1(1), spline_o%n1, &
1049  & spline_o%x2(1), spline_o%n2, &
1050  & spline_o%ilin1, spline_o%ilin2, &
1051  & spline_o%fspl(1,1,1), spline_o%n1, &
1052  & ict, f, ifail)
1053 
1054  endif
1055  if(ifail /= 0) ier = 96
1056 
1057 end subroutine ezspline_derivative2_r4
1058 
1059 
1060 
1061 subroutine ezspline_derivative2_array_r4(spline_o, i1, i2, &
1062  & k1, k2, p1, p2, f, ier)
1063  use ezspline_obj
1064  implicit none
1065  type(EZspline2_r4) spline_o
1066  integer, intent(in) :: i1, i2, k1, k2
1067  real(ezspline_r4), intent(in) :: p1(k1), p2(k2)
1068  real(ezspline_r4), intent(out) :: f(k1,k2)
1069  integer, intent(out) :: ier
1070 
1071  integer ifail
1072 
1073  integer ict(6)
1074  integer ict_herm(4)
1075 
1076  real(ezspline_r4), dimension(:), allocatable :: p1_cloud, p2_cloud
1077  integer k12
1078  integer:: iwarn=0
1079 
1080  ier = 0
1081  ifail=0
1082  if( .not.ezspline_allocated(spline_o) .or. spline_o%isReady /= 1) then
1083  ier = 94
1084  return
1085  endif
1086 
1087  f = 0.0_ezspline_r4
1088  if(i1<0 .OR. i2<0) then
1089  ier = 12
1090  return
1091  endif
1092  if(max(i1,i2) > 3) then
1093  ier = 13
1094  return
1095  endif
1096 
1097  if(max(i1,i2) > 1 .AND. (spline_o%isHermite==1 .or. spline_o%isLinear==1)) then
1098  ier = 24; if(spline_o%isLinear==1) ier = 46
1099  return
1100  endif
1101 
1102  k12 = k1*k2
1103  allocate(p1_cloud(k12), p2_cloud(k12), stat=ifail)
1104  if(ifail /= 0) then
1105  ier = 32
1106  return
1107  endif
1108 
1109  p1_cloud = reshape( &
1110  & source=spread(source=p1, dim=2, ncopies=k2), &
1111  & shape=(/k12/))
1112  p2_cloud = reshape( &
1113  & source=spread(source=p2, dim=1, ncopies=k1), &
1114  & shape=(/k12/))
1115 
1116  if(spline_o%isHybrid == 1) then
1117 
1118  call ezmake_ict2(i1,i2,ict)
1119 
1120  call vecintrp2d(ict, k12, p1_cloud, p2_cloud, k12, f, &
1121  & spline_o%n1, spline_o%x1pkg(1,1), &
1122  & spline_o%n2, spline_o%x2pkg(1,1), &
1123  & spline_o%hspline, spline_o%fspl(1,1,1), &
1124  & size(spline_o%fspl,1), size(spline_o%fspl,2), &
1125  & size(spline_o%fspl,3), &
1126  & iwarn, ifail)
1127 
1128  else if(spline_o%isHermite==1 .or. spline_o%isLinear==1)then
1129  if(i1 == 1 .AND. i2 == 0) then
1130  ict_herm = (/0, 1, 0, 0/)
1131  else if(i1 == 0 .AND. i2 == 1) then
1132  ict_herm = (/0, 0, 1, 0/)
1133  else if(i1 == 1 .AND. i2 == 1) then
1134  ict_herm = (/0, 0, 0, 1/)
1135  else
1136  ict_herm = (/1, 0, 0, 0/)
1137  endif
1138 
1139  if(spline_o%isLinear == 0) then
1140  call vecherm2(ict_herm, k12, p1_cloud, p2_cloud, k12, f, &
1141  & spline_o%n1, spline_o%x1pkg(1,1), &
1142  & spline_o%n2, spline_o%x2pkg(1,1), &
1143  & spline_o%fspl(1,1,1), spline_o%n1, &
1144  & iwarn, ifail)
1145  else
1146  call vecpc2(ict_herm, k12, p1_cloud, p2_cloud, k12, f, &
1147  & spline_o%n1, spline_o%x1pkg(1,1), &
1148  & spline_o%n2, spline_o%x2pkg(1,1), &
1149  & spline_o%fspl(1,1,1), spline_o%n1, &
1150  & iwarn, ifail)
1151  endif
1152 
1153 
1154  else
1155 
1156  call ezmake_ict2(i1,i2,ict)
1157 
1158  call vecbicub(ict, k12, p1_cloud, p2_cloud, k12, f, &
1159  & spline_o%n1,spline_o%x1pkg(1,1), &
1160  & spline_o%n2,spline_o%x2pkg(1,1), &
1161  & spline_o%fspl(1,1,1), spline_o%n1, &
1162  & iwarn, ifail)
1163 
1164  endif
1165 
1166 
1167  if(ifail /= 0) ier = 96
1168 
1169  deallocate(p1_cloud, p2_cloud, stat=ifail)
1170  if(ifail /= 0) then
1171  ier = 33
1172  return
1173  endif
1174 
1175 
1176 end subroutine ezspline_derivative2_array_r4
1177 
1178 
1179 subroutine ezspline_derivative2_cloud_r4(spline_o, i1, i2, &
1180  & k, p1, p2, f, ier)
1181  use ezspline_obj
1182  implicit none
1183  type(EZspline2_r4) spline_o
1184  integer, intent(in) :: i1, i2, k
1185  real(ezspline_r4), intent(in) :: p1(k), p2(k)
1186  real(ezspline_r4), intent(out) :: f(k)
1187  integer, intent(out) :: ier
1188 
1189  integer ifail
1190 
1191  integer ict(6)
1192  integer ict_herm(4)
1193  integer:: iwarn=0
1194 
1195  ier = 0
1196  ifail=0
1197  if( .not.ezspline_allocated(spline_o) .or. spline_o%isReady /= 1) then
1198  ier = 94
1199  return
1200  endif
1201 
1202  f = 0.0_ezspline_r4
1203  if(i1<0 .OR. i2<0) then
1204  ier = 12
1205  return
1206  endif
1207  if(max(i1,i2) > 3) then
1208  ier = 13
1209  return
1210  endif
1211 
1212  if(max(i1,i2) > 1 .AND. (spline_o%isHermite==1 .or. spline_o%isLinear==1)) then
1213  ier = 24; if(spline_o%isLinear==1) ier = 46
1214  return
1215  endif
1216 
1217  if(spline_o%isHybrid == 1) then
1218 
1219  call ezmake_ict2(i1,i2,ict)
1220 
1221  call vecintrp2d(ict, k, p1, p2, k, f, &
1222  & spline_o%n1, spline_o%x1pkg(1,1), &
1223  & spline_o%n2, spline_o%x2pkg(1,1), &
1224  & spline_o%hspline, spline_o%fspl(1,1,1), &
1225  & size(spline_o%fspl,1), size(spline_o%fspl,2), &
1226  & size(spline_o%fspl,3), &
1227  & iwarn, ifail)
1228 
1229  else if(spline_o%isHermite==1 .or. spline_o%isLinear==1)then
1230  if(i1 == 1 .AND. i2 == 0) then
1231  ict_herm = (/0, 1, 0, 0/)
1232  else if(i1 == 0 .AND. i2 == 1) then
1233  ict_herm = (/0, 0, 1, 0/)
1234  else if(i1 == 1 .AND. i2 == 1) then
1235  ict_herm = (/0, 0, 0, 1/)
1236  else
1237  ict_herm = (/1, 0, 0, 0/)
1238  endif
1239 
1240  if(spline_o%isLinear == 0) then
1241  call vecherm2(ict_herm, k, p1, p2, &
1242  & k, f, &
1243  & spline_o%n1, spline_o%x1pkg(1,1), &
1244  & spline_o%n2, spline_o%x2pkg(1,1), &
1245  & spline_o%fspl(1,1,1), spline_o%n1, &
1246  & iwarn, ifail)
1247  else
1248  call vecpc2(ict_herm, k, p1, p2, &
1249  & k, f, &
1250  & spline_o%n1, spline_o%x1pkg(1,1), &
1251  & spline_o%n2, spline_o%x2pkg(1,1), &
1252  & spline_o%fspl(1,1,1), spline_o%n1, &
1253  & iwarn, ifail)
1254  endif
1255 
1256  else
1257 
1258  call ezmake_ict2(i1,i2,ict)
1259 
1260  call vecbicub(ict, k, p1, p2, &
1261  & k, f, &
1262  & spline_o%n1, spline_o%x1pkg(1,1), &
1263  & spline_o%n2, spline_o%x2pkg(1,1), &
1264  & spline_o%fspl(1,1,1), spline_o%n1, &
1265  & iwarn, ifail)
1266 
1267  endif
1268 
1269  if(ifail /= 0) ier = 96
1270 
1271 end subroutine ezspline_derivative2_cloud_r4
1272 
1273 !!!
1274 !!! 3-D
1275 !!!
1276 
1277 
1278 subroutine ezspline_derivative3_r4(spline_o, i1, i2, i3, p1, p2, p3, f, ier)
1279  use ezspline_obj
1280  implicit none
1281  type(EZspline3_r4) spline_o
1282  integer, intent(in) :: i1, i2, i3
1283  real(ezspline_r4), intent(in) :: p1, p2, p3
1284  real(ezspline_r4), intent(out) :: f
1285  integer, intent(out) :: ier
1286 
1287  integer ifail
1288 
1289  integer ict(10)
1290  integer ict_herm(8)
1291 
1292  ier = 0
1293  ifail=0
1294  if( .not.ezspline_allocated(spline_o) .or. spline_o%isReady /= 1) then
1295  ier = 94
1296  return
1297  endif
1298 
1299  f = 0.0_ezspline_r4
1300  if(i1<0 .OR. i2<0 .OR. i3<0) then
1301  ier = 12
1302  return
1303  endif
1304  if(max(i1,max(i2,i3)) > 3) then
1305  ier = 13
1306  return
1307  endif
1308 
1309  if(maxval((/i1,i2,i3/)) > 1 .AND. (spline_o%isHermite==1 .or. spline_o%isLinear==1)) then
1310  ier = 24; if(spline_o%isLinear==1) ier = 46
1311  return
1312  endif
1313 
1314  if(spline_o%isHybrid == 1) then
1315 
1316  call ezmake_ict3(i1,i2,i3,ict)
1317 
1318  call evintrp3d(p1, p2, p3, &
1319  & spline_o%x1(1), spline_o%n1, &
1320  & spline_o%x2(1), spline_o%n2, &
1321  & spline_o%x3(1), spline_o%n3, &
1322  & spline_o%hspline, spline_o%fspl(1,1,1,1), &
1323  & size(spline_o%fspl,1), size(spline_o%fspl,2), &
1324  & size(spline_o%fspl,3), size(spline_o%fspl,4), &
1325  & ict, f, ifail)
1326 
1327  else if(spline_o%isHermite==1 .or. spline_o%isLinear==1)then
1328  if(i1 == 1 .AND. i2 == 0 .AND. i3 == 0) then
1329  ict_herm = (/0, 1, 0, 0, 0, 0, 0, 0/)
1330  else if(i1 == 0 .AND. i2 == 1 .AND. i3 == 0) then
1331  ict_herm = (/0, 0, 1, 0, 0, 0, 0, 0/)
1332  else if(i1 == 0 .AND. i2 == 0 .AND. i3 == 1) then
1333  ict_herm = (/0, 0, 0, 1, 0, 0, 0, 0/)
1334  else if(i1 == 1 .AND. i2 == 1 .AND. i3 == 0) then
1335  ict_herm = (/0, 0, 0, 0, 1, 0, 0, 0/)
1336  else if(i1 == 1 .AND. i2 == 0 .AND. i3 == 1) then
1337  ict_herm = (/0, 0, 0, 0, 0, 1, 0, 0/)
1338  else if(i1 == 0 .AND. i2 == 1 .AND. i3 == 1) then
1339  ict_herm = (/0, 0, 0, 0, 0, 0, 1, 0/)
1340  else if(i1 == 1 .AND. i2 == 1 .AND. i3 == 1) then
1341  ict_herm = (/0, 0, 0, 0, 0, 0, 0, 1/)
1342  else
1343  ict_herm = (/1, 0, 0, 0, 0, 0, 0, 0/)
1344  endif
1345 
1346  if(spline_o%isLinear == 0) then
1347  call herm3ev(p1, p2, p3, &
1348  & spline_o%x1(1), spline_o%n1, &
1349  & spline_o%x2(1), spline_o%n2, &
1350  & spline_o%x3(1), spline_o%n3, &
1351  & spline_o%ilin1, spline_o%ilin2, spline_o%ilin3, &
1352  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
1353  & ict_herm, f, ifail)
1354  else
1355  call pc3ev(p1, p2, p3, &
1356  & spline_o%x1(1), spline_o%n1, &
1357  & spline_o%x2(1), spline_o%n2, &
1358  & spline_o%x3(1), spline_o%n3, &
1359  & spline_o%ilin1, spline_o%ilin2, spline_o%ilin3, &
1360  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
1361  & ict_herm, f, ifail)
1362  endif
1363 
1364  else
1365 
1366  call ezmake_ict3(i1,i2,i3,ict)
1367 
1368  call evtricub(p1, p2, p3, &
1369  & spline_o%x1(1), spline_o%n1, &
1370  & spline_o%x2(1), spline_o%n2, &
1371  & spline_o%x3(1), spline_o%n3, &
1372  & spline_o%ilin1, spline_o%ilin2, spline_o%ilin3, &
1373  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
1374  & ict, f, ifail)
1375 
1376 
1377  endif
1378  if(ifail /= 0) ier = 96
1379 
1380 end subroutine ezspline_derivative3_r4
1381 
1382 
1383 subroutine ezspline_derivative3_array_r4(spline_o, i1, i2, i3, &
1384  & k1, k2, k3, p1, p2, p3, f, ier)
1385  use ezspline_obj
1386  implicit none
1387  type(EZspline3_r4) spline_o
1388  integer, intent(in) :: i1, i2, i3, k1, k2, k3
1389  real(ezspline_r4), intent(in) :: p1(k1), p2(k2), p3(k3)
1390  real(ezspline_r4), intent(out) :: f(k1,k2,k3)
1391  integer, intent(out) :: ier
1392 
1393  integer ifail
1394 
1395  integer ict(10)
1396  integer ict_herm(8)
1397 
1398  real(ezspline_r4), dimension(:), allocatable :: p1_cloud, p2_cloud, p3_cloud
1399  integer k123
1400  integer:: iwarn=0
1401 
1402  ier = 0
1403  ifail=0
1404  if( .not.ezspline_allocated(spline_o) .or. spline_o%isReady /= 1) then
1405  ier = 94
1406  return
1407  endif
1408 
1409  f = 0.0_ezspline_r4
1410  if(i1<0 .OR. i2<0 .OR. i3<0) then
1411  ier = 12
1412  return
1413  endif
1414  if(max(i1,max(i2,i3)) > 3) then
1415  ier = 13
1416  return
1417  endif
1418 
1419  if(maxval((/i1,i2,i3/)) > 1 .AND. (spline_o%isHermite==1 .or. spline_o%isLinear==1)) then
1420  ier = 24; if(spline_o%isLinear==1) ier = 46
1421  return
1422  endif
1423 
1424  k123 = k1*k2*k3
1425  allocate(p1_cloud(k123), p2_cloud(k123), p3_cloud(k123), stat=ifail)
1426  if(ifail /= 0) then
1427  ier = 32
1428  return
1429  endif
1430 
1431  p1_cloud = reshape(source=spread( &
1432  & source=spread(source=p1, dim=2, ncopies=k2), &
1433  & dim=3, ncopies=k3), shape=(/k123/))
1434  p2_cloud = reshape(source=spread( &
1435  & source=spread(source=p2, dim=1, ncopies=k1), &
1436  & dim=3, ncopies=k3), shape=(/k123/))
1437  p3_cloud = reshape(source=spread( &
1438  & source=spread(source=p3, dim=1, ncopies=k1), &
1439  & dim=2, ncopies=k2), shape=(/k123/))
1440 
1441  if(spline_o%isHybrid == 1) then
1442 
1443  call ezmake_ict3(i1,i2,i3,ict)
1444 
1445  call vecintrp3d(ict, k123, p1_cloud, p2_cloud, p3_cloud, k123, f, &
1446  & spline_o%n1, spline_o%x1pkg(1,1), &
1447  & spline_o%n2, spline_o%x2pkg(1,1), &
1448  & spline_o%n3, spline_o%x3pkg(1,1), &
1449  & spline_o%hspline, spline_o%fspl(1,1,1,1), &
1450  & size(spline_o%fspl,1), size(spline_o%fspl,2), &
1451  & size(spline_o%fspl,3), size(spline_o%fspl,4), &
1452  & iwarn, ifail)
1453 
1454  else if(spline_o%isHermite==1 .or. spline_o%isLinear==1)then
1455  if(i1 == 1 .AND. i2 == 0 .AND. i3 == 0) then
1456  ict_herm = (/0, 1, 0, 0, 0, 0, 0, 0/)
1457  else if(i1 == 0 .AND. i2 == 1 .AND. i3 == 0) then
1458  ict_herm = (/0, 0, 1, 0, 0, 0, 0, 0/)
1459  else if(i1 == 0 .AND. i2 == 0 .AND. i3 == 1) then
1460  ict_herm = (/0, 0, 0, 1, 0, 0, 0, 0/)
1461  else if(i1 == 1 .AND. i2 == 1 .AND. i3 == 0) then
1462  ict_herm = (/0, 0, 0, 0, 1, 0, 0, 0/)
1463  else if(i1 == 1 .AND. i2 == 0 .AND. i3 == 1) then
1464  ict_herm = (/0, 0, 0, 0, 0, 1, 0, 0/)
1465  else if(i1 == 0 .AND. i2 == 1 .AND. i3 == 1) then
1466  ict_herm = (/0, 0, 0, 0, 0, 0, 1, 0/)
1467  else if(i1 == 1 .AND. i2 == 1 .AND. i3 == 1) then
1468  ict_herm = (/0, 0, 0, 0, 0, 0, 0, 1/)
1469  else
1470  ict_herm = (/1, 0, 0, 0, 0, 0, 0, 0/)
1471  endif
1472 
1473  if(spline_o%isLinear == 0) then
1474  call vecherm3(ict_herm, k123, p1_cloud, p2_cloud, p3_cloud, k123, f, &
1475  & spline_o%n1, spline_o%x1pkg(1,1), &
1476  & spline_o%n2, spline_o%x2pkg(1,1), &
1477  & spline_o%n3, spline_o%x3pkg(1,1), &
1478  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
1479  & iwarn, ifail)
1480  else
1481  call vecpc3(ict_herm, k123, p1_cloud, p2_cloud, p3_cloud, k123, f, &
1482  & spline_o%n1, spline_o%x1pkg(1,1), &
1483  & spline_o%n2, spline_o%x2pkg(1,1), &
1484  & spline_o%n3, spline_o%x3pkg(1,1), &
1485  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
1486  & iwarn, ifail)
1487  endif
1488 
1489  else
1490 
1491  call ezmake_ict3(i1,i2,i3,ict)
1492 
1493  call vectricub(ict, k123, p1_cloud, p2_cloud, p3_cloud, k123, f, &
1494  & spline_o%n1,spline_o%x1pkg(1,1), &
1495  & spline_o%n2,spline_o%x2pkg(1,1), &
1496  & spline_o%n3,spline_o%x3pkg(1,1), &
1497  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
1498  & iwarn, ifail)
1499 
1500  endif
1501 
1502 
1503  if(ifail /= 0) ier = 96
1504 
1505  deallocate(p1_cloud, p2_cloud, p3_cloud, stat=ifail)
1506  if(ifail /= 0) then
1507  ier = 33
1508  return
1509  endif
1510 
1511 
1512 end subroutine ezspline_derivative3_array_r4
1513 
1514 subroutine ezspline_derivative3_cloud_r4(spline_o, i1, i2, i3, &
1515  & k, p1, p2, p3, f, ier)
1516  use ezspline_obj
1517  implicit none
1518  type(EZspline3_r4) spline_o
1519  integer, intent(in) :: i1, i2, i3, k
1520  real(ezspline_r4), intent(in) :: p1(k), p2(k), p3(k)
1521  real(ezspline_r4), intent(out) :: f(k)
1522  integer, intent(out) :: ier
1523 
1524  integer ifail
1525 
1526  integer ict(10)
1527  integer ict_herm(8)
1528  integer:: iwarn=0
1529 
1530  ier = 0
1531  ifail=0
1532  if( .not.ezspline_allocated(spline_o) .or. spline_o%isReady /= 1) then
1533  ier = 94
1534  return
1535  endif
1536 
1537  f = 0.0_ezspline_r4
1538  if(i1<0 .OR. i2<0 .OR. i3<0) then
1539  ier = 12
1540  return
1541  endif
1542  if(max(i1,max(i2,i3)) > 3) then
1543  ier = 13
1544  return
1545  endif
1546 
1547  if(maxval((/i1,i2,i3/)) > 1 .AND. (spline_o%isHermite==1 .or. spline_o%isLinear==1)) then
1548  ier = 24; if(spline_o%isLinear==1) ier = 46
1549  return
1550  endif
1551 
1552  if(spline_o%isHybrid == 1) then
1553 
1554  call ezmake_ict3(i1,i2,i3,ict)
1555 
1556  call vecintrp3d(ict, k, p1, p2, p3, k, f, &
1557  & spline_o%n1, spline_o%x1pkg(1,1), &
1558  & spline_o%n2, spline_o%x2pkg(1,1), &
1559  & spline_o%n3, spline_o%x3pkg(1,1), &
1560  & spline_o%hspline, spline_o%fspl(1,1,1,1), &
1561  & size(spline_o%fspl,1), size(spline_o%fspl,2), &
1562  & size(spline_o%fspl,3), size(spline_o%fspl,4), &
1563  & iwarn, ifail)
1564 
1565  else if(spline_o%isHermite==1 .or. spline_o%isLinear==1)then
1566  if(i1 == 1 .AND. i2 == 0 .AND. i3 == 0) then
1567  ict_herm = (/0, 1, 0, 0, 0, 0, 0, 0/)
1568  else if(i1 == 0 .AND. i2 == 1 .AND. i3 == 0) then
1569  ict_herm = (/0, 0, 1, 0, 0, 0, 0, 0/)
1570  else if(i1 == 0 .AND. i2 == 0 .AND. i3 == 1) then
1571  ict_herm = (/0, 0, 0, 1, 0, 0, 0, 0/)
1572  else if(i1 == 1 .AND. i2 == 1 .AND. i3 == 0) then
1573  ict_herm = (/0, 0, 0, 0, 1, 0, 0, 0/)
1574  else if(i1 == 1 .AND. i2 == 0 .AND. i3 == 1) then
1575  ict_herm = (/0, 0, 0, 0, 0, 1, 0, 0/)
1576  else if(i1 == 0 .AND. i2 == 1 .AND. i3 == 1) then
1577  ict_herm = (/0, 0, 0, 0, 0, 0, 1, 0/)
1578  else if(i1 == 1 .AND. i2 == 1 .AND. i3 == 1) then
1579  ict_herm = (/0, 0, 0, 0, 0, 0, 0, 1/)
1580  else
1581  ict_herm = (/1, 0, 0, 0, 0, 0, 0, 0/)
1582  endif
1583 
1584  if(spline_o%isLinear == 0) then
1585  call vecherm3(ict_herm, k, p1, p2, p3, &
1586  & k, f, &
1587  & spline_o%n1, spline_o%x1pkg(1,1), &
1588  & spline_o%n2, spline_o%x2pkg(1,1), &
1589  & spline_o%n3, spline_o%x3pkg(1,1), &
1590  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
1591  & iwarn, ifail)
1592  else
1593  call vecpc3(ict_herm, k, p1, p2, p3, &
1594  & k, f, &
1595  & spline_o%n1, spline_o%x1pkg(1,1), &
1596  & spline_o%n2, spline_o%x2pkg(1,1), &
1597  & spline_o%n3, spline_o%x3pkg(1,1), &
1598  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
1599  & iwarn, ifail)
1600  endif
1601 
1602  else
1603 
1604  call ezmake_ict3(i1,i2,i3,ict)
1605 
1606  call vectricub(ict, k, p1, p2, p3, &
1607  & k, f, &
1608  & spline_o%n1, spline_o%x1pkg(1,1), &
1609  & spline_o%n2, spline_o%x2pkg(1,1), &
1610  & spline_o%n3, spline_o%x3pkg(1,1), &
1611  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
1612  & iwarn, ifail)
1613 
1614  endif
1615 
1616  if(ifail /= 0) ier = 96
1617 
1618 end subroutine ezspline_derivative3_cloud_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