V3FIT
•All Classes Namespaces Files Functions Variables Enumerations Macros Pages
ezspline_cinterp.f90
1 subroutine ezspline_cinterp(spline_o, k, p1, p2, p3, f, ier)
2  ! list of coordinate triplets
3  ! this version of interp that has been inlined for better performance on the
4  ! crays (pletzer@pppl.gov) Thu Jun 15 08:52:46 PDT 2000
5  ! NOTE:
6  ! This routine does not require the -dp switch on Crays since
7  ! all constant have been made kind=r8.
8 
9  ! DMC Apr 2007 -- only the tricubic SPLINE evaluation has been inlined;
10  ! if a different method interpolation is used, the standard f77-style
11  ! pspline subroutine is called.
12 
13  use ezspline_obj
14  implicit none
15  type(EZspline3_r8) spline_o
16  integer, intent(in) :: k
17  real(ezspline_r8), intent(in) :: p1(k), p2(k), p3(k) ! location arrays
18  real(ezspline_r8), intent(out):: f(k) ! interpolant array
19  integer, intent(out) :: ier
20 
21  integer :: ifail
22  integer, parameter :: ict(10)=(/1,0,0,0,0,0,0,0,0,0/)
23  integer:: iwarn=0
24 
25  integer ix(k) ! zone indices {j}
26  REAL(EZSPLINE_R8) dxn(k) ! normalized displacements w/in zones
27  REAL(EZSPLINE_R8) hx(k) ! h(j) vector
28  REAL(EZSPLINE_R8) hxi(k) ! 1/h(j) vector
29 !
30  integer iy(k) ! zone indices {j}
31  REAL(EZSPLINE_R8) dyn(k) ! normalized displacements w/in zones
32  REAL(EZSPLINE_R8) hy(k) ! h(j) vector
33  REAL(EZSPLINE_R8) hyi(k) ! 1/h(j) vector
34 !
35  integer iz(k) ! zone indices {j}
36  REAL(EZSPLINE_R8) dzn(k) ! normalized displacements w/in zones
37  REAL(EZSPLINE_R8) hz(k) ! h(j) vector
38  REAL(EZSPLINE_R8) hzi(k) ! 1/h(j) vector
39  integer v
40 !
41  REAL(EZSPLINE_R8) , parameter :: sixth=0.166666666666666667_ezspline_r8
42  REAL(EZSPLINE_R8) , parameter :: one=1.0_ezspline_r8
43  REAL(EZSPLINE_R8) , parameter :: three=3.0_ezspline_r8
44 
45  integer j1, j2 , j3
46  REAL(EZSPLINE_R8) z36th,z216th,xp,xpi,xp2,xpi2,cx,cxi,hx2,cxd,cxdi,yp, &
47  & ypi,yp2,ypi2,cy,cyi,hy2,cyd,cydi,zp,zpi,zp2,zpi2,cz,czi,hz2, &
48  & czd,czdi,somme
49 
50  ier = 0
51  ifail=0
52  if(spline_o%isReady /= 1) then
53  ier = 94
54  return
55  endif
56 
57  if (spline_o%isLinear == 1) then
58 
59  call r8vecpc3(ict, k, p1, p2, p3, k, f, &
60  & spline_o%n1, spline_o%x1pkg(1,1), &
61  & spline_o%n2, spline_o%x2pkg(1,1), &
62  & spline_o%n3, spline_o%x3pkg(1,1), &
63  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
64  & iwarn,ifail)
65 
66  else if (spline_o%isHybrid == 1) then
67 
68  call r8vecintrp3d(ict, k, p1, p2, p3, k, f, &
69  & spline_o%n1, spline_o%x1pkg(1,1), &
70  & spline_o%n2, spline_o%x2pkg(1,1), &
71  & spline_o%n3, spline_o%x3pkg(1,1), &
72  & spline_o%hspline, &
73  & spline_o%fspl(1,1,1,1), size(spline_o%fspl,1), &
74  & size(spline_o%fspl,2), size(spline_o%fspl,3), &
75  & size(spline_o%fspl,4), iwarn, ifail)
76 
77  else if (spline_o%isHermite ==0) then
78  !
79 !!$ call r8vectricub(ict, k, p1, p2, p3, k, f, &
80 !!$ & spline_o%n1,spline_o%x1pkg(1,1), &
81 !!$ & spline_o%n2,spline_o%x2pkg(1,1), &
82 !!$ & spline_o%n3,spline_o%x3pkg(1,1), &
83 !!$ & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
84 !!$ & iwarn, ifail)
85 
86  call r8xlookup(k,p1,spline_o%n1,spline_o%x1pkg(1,1),2,ix,dxn,hx,hxi,ifail)
87  call r8xlookup(k,p2,spline_o%n2,spline_o%x2pkg(1,1),2,iy,dyn,hy,hyi,ifail)
88  call r8xlookup(k,p3,spline_o%n3,spline_o%x3pkg(1,1),2,iz,dzn,hz,hzi,ifail)
89 
90 !!$ call r8fvtricub(ict,k,k,f,ix,iy,iz,dxn,dyn,dzn, &
91 !!$ & hx,hxi,hy,hyi,hz,hzi,spline_o%fspl(1,1,1,1), &
92 !!$ & spline_o%n1, spline_o%n2, spline_o%n3)
93 
94 
95  z36th=sixth*sixth
96  z216th=sixth*sixth*sixth
97 !
98 ! prepare useful parameters...
99 !
100  do v=1,k
101  j1=ix(v)
102  j2=iy(v)
103  j3=iz(v)
104 !
105 ! ...in x direction
106 !
107  xp=dxn(v)
108  xpi=one-xp
109  xp2=xp*xp
110  xpi2=xpi*xpi
111 !
112 !!$ if((ict(1).eq.1).or.(ict(3).eq.1).or.(ict(4).eq.1).or. &
113 !!$ & (ict(6).eq.1).or.(ict(7).eq.1).or.(ict(10).eq.1)) then
114  cx=xp*(xp2-one)
115  cxi=xpi*(xpi2-one)
116  hx2=hx(v)*hx(v)
117 !!$ endif
118 !!$ if((ict(2).eq.1).or.(ict(8).eq.1).or.(ict(9).eq.1)) then
119 !!$ cxd=THREE*xp2-ONE
120 !!$ cxdi=-THREE*xpi2+ONE
121 !!$ endif
122 !
123 ! ...and in y direction
124 !
125  yp=dyn(v)
126  ypi=one-yp
127  yp2=yp*yp
128  ypi2=ypi*ypi
129 !
130 !!$ if((ict(1).eq.1).or.(ict(2).eq.1).or.(ict(4).eq.1).or. &
131 !!$ & (ict(5).eq.1).or.(ict(7).eq.1).or.(ict(9).eq.1)) then
132  cy=yp*(yp2-one)
133  cyi=ypi*(ypi2-one)
134  hy2=hy(v)*hy(v)
135 !!$ endif
136 !!$ if((ict(3).eq.1).or.(ict(8).eq.1).or.(ict(10).eq.1)) then
137 !!$ cyd=THREE*yp2-ONE
138 !!$ cydi=-THREE*ypi2+ONE
139 !!$ endif
140 !
141 ! ...and in z direction
142 !
143  zp=dzn(v)
144  zpi=one-zp
145  zp2=zp*zp
146  zpi2=zpi*zpi
147 !
148 !!$ if((ict(1).eq.1).or.(ict(2).eq.1).or.(ict(3).eq.1).or. &
149 !!$ & (ict(5).eq.1).or.(ict(6).eq.1).or.(ict(8).eq.1)) then
150  cz=zp*(zp2-one)
151  czi=zpi*(zpi2-one)
152  hz2=hz(v)*hz(v)
153 !!$ endif
154 !!$ if((ict(4).eq.1).or.(ict(9).eq.1).or.(ict(10).eq.1)) then
155 !!$ czd=THREE*zp2-ONE
156 !!$ czdi=-THREE*zpi2+ONE
157 !!$ endif
158 !
159 ! get desired values:
160 !
161  if(ict(1).eq.1) then
162 !
163 ! function value:
164 !
165  somme=( &
166  & zpi*( &
167  & xpi*(ypi*spline_o%fspl(1,j1,j2,j3) +yp*spline_o%fspl(1,j1,j2+1,j3))+ &
168  & xp*(ypi*spline_o%fspl(1,j1+1,j2,j3)+yp*spline_o%fspl(1,j1+1,j2+1,j3))) &
169  & +zp*( &
170  & xpi*(ypi*spline_o%fspl(1,j1,j2,j3+1) +yp*spline_o%fspl(1,j1,j2+1,j3+1))+ &
171  & xp*(ypi*spline_o%fspl(1,j1+1,j2,j3+1)+yp*spline_o%fspl(1,j1+1,j2+1,j3+1))))
172 !
173  somme=somme+sixth*hx2*( &
174  & zpi*( &
175  & cxi*(ypi*spline_o%fspl(2,j1,j2,j3) +yp*spline_o%fspl(2,j1,j2+1,j3))+ &
176  & cx*(ypi*spline_o%fspl(2,j1+1,j2,j3)+yp*spline_o%fspl(2,j1+1,j2+1,j3))) &
177  & +zp*( &
178  & cxi*(ypi*spline_o%fspl(2,j1,j2,j3+1) +yp*spline_o%fspl(2,j1,j2+1,j3+1))+ &
179  & cx*(ypi*spline_o%fspl(2,j1+1,j2,j3+1)+yp*spline_o%fspl(2,j1+1,j2+1,j3+1))))
180 !
181  somme=somme+sixth*hy2*( &
182  & zpi*( &
183  & xpi*(cyi*spline_o%fspl(3,j1,j2,j3) +cy*spline_o%fspl(3,j1,j2+1,j3))+ &
184  & xp*(cyi*spline_o%fspl(3,j1+1,j2,j3)+cy*spline_o%fspl(3,j1+1,j2+1,j3))) &
185  & +zp*( &
186  & xpi*(cyi*spline_o%fspl(3,j1,j2,j3+1) +cy*spline_o%fspl(3,j1,j2+1,j3+1))+ &
187  & xp*(cyi*spline_o%fspl(3,j1+1,j2,j3+1)+cy*spline_o%fspl(3,j1+1,j2+1,j3+1))))
188 !
189  somme=somme+sixth*hz2*( &
190  & czi*( &
191  & xpi*(ypi*spline_o%fspl(4,j1,j2,j3) +yp*spline_o%fspl(4,j1,j2+1,j3))+ &
192  & xp*(ypi*spline_o%fspl(4,j1+1,j2,j3)+yp*spline_o%fspl(4,j1+1,j2+1,j3))) &
193  & +cz*( &
194  & xpi*(ypi*spline_o%fspl(4,j1,j2,j3+1) +yp*spline_o%fspl(4,j1,j2+1,j3+1))+ &
195  & xp*(ypi*spline_o%fspl(4,j1+1,j2,j3+1)+yp*spline_o%fspl(4,j1+1,j2+1,j3+1))))
196 !
197  somme=somme+z36th*hx2*hy2*( &
198  & zpi*( &
199  & cxi*(cyi*spline_o%fspl(5,j1,j2,j3) +cy*spline_o%fspl(5,j1,j2+1,j3))+ &
200  & cx*(cyi*spline_o%fspl(5,j1+1,j2,j3)+cy*spline_o%fspl(5,j1+1,j2+1,j3))) &
201  & +zp*( &
202  & cxi*(cyi*spline_o%fspl(5,j1,j2,j3+1) +cy*spline_o%fspl(5,j1,j2+1,j3+1))+ &
203  & cx*(cyi*spline_o%fspl(5,j1+1,j2,j3+1)+cy*spline_o%fspl(5,j1+1,j2+1,j3+1))))
204 !
205  somme=somme+z36th*hx2*hz2*( &
206  & czi*( &
207  & cxi*(ypi*spline_o%fspl(6,j1,j2,j3) +yp*spline_o%fspl(6,j1,j2+1,j3))+ &
208  & cx*(ypi*spline_o%fspl(6,j1+1,j2,j3)+yp*spline_o%fspl(6,j1+1,j2+1,j3))) &
209  & +cz*( &
210  & cxi*(ypi*spline_o%fspl(6,j1,j2,j3+1) +yp*spline_o%fspl(6,j1,j2+1,j3+1))+ &
211  & cx*(ypi*spline_o%fspl(6,j1+1,j2,j3+1)+yp*spline_o%fspl(6,j1+1,j2+1,j3+1))))
212 !
213  somme=somme+z36th*hy2*hz2*( &
214  & czi*( &
215  & xpi*(cyi*spline_o%fspl(7,j1,j2,j3) +cy*spline_o%fspl(7,j1,j2+1,j3))+ &
216  & xp*(cyi*spline_o%fspl(7,j1+1,j2,j3)+cy*spline_o%fspl(7,j1+1,j2+1,j3))) &
217  & +cz*( &
218  & xpi*(cyi*spline_o%fspl(7,j1,j2,j3+1) +cy*spline_o%fspl(7,j1,j2+1,j3+1))+ &
219  & xp*(cyi*spline_o%fspl(7,j1+1,j2,j3+1)+cy*spline_o%fspl(7,j1+1,j2+1,j3+1))))
220 !
221  somme=somme+z216th*hx2*hy2*hz2*( &
222  & czi*( &
223  & cxi*(cyi*spline_o%fspl(8,j1,j2,j3) +cy*spline_o%fspl(8,j1,j2+1,j3))+ &
224  & cx*(cyi*spline_o%fspl(8,j1+1,j2,j3)+cy*spline_o%fspl(8,j1+1,j2+1,j3))) &
225  & +cz*( &
226  & cxi*(cyi*spline_o%fspl(8,j1,j2,j3+1) +cy*spline_o%fspl(8,j1,j2+1,j3+1))+ &
227  & cx*(cyi*spline_o%fspl(8,j1+1,j2,j3+1)+cy*spline_o%fspl(8,j1+1,j2+1,j3+1))))
228 !
229  f(v)=somme
230  else
231  print*,'ERROR ict(1) must be 1 (interpolation)'
232  endif
233 !
234 !
235  enddo ! vector loop
236 
237 
238  if(ifail /= 0) ier = 97
239 
240  else
241 
242  call r8vecherm3(ict, k, p1, p2, p3, k, f, &
243  & spline_o%n1, spline_o%x1pkg(1,1), &
244  & spline_o%n2, spline_o%x2pkg(1,1), &
245  & spline_o%n3, spline_o%x3pkg(1,1), &
246  & spline_o%fspl(1,1,1,1), spline_o%n1, spline_o%n2, &
247  & iwarn,ifail)
248 
249  endif
250 
251  if(ifail /= 0) ier = 97
252 
253 end subroutine
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17