1subroutine trispline(f,x,y,z,n1,n2,n3,t1,t2,t3,fx,fy,fxy)
2!subroutine trispline(f,x,y,z,n1,n2,n3,t1,t2,t3,fx,fy,fz,fxy,fxz,fyz,fxyz)
3!     constructs interpolation coefficients for a periodic
4!     function f evaluated on a
5!     three dimensional grid (x,y,z) with size n1 x n2 x n3.
6!     t1,t2,t3 are the periods.
7implicit none
8integer n1,n2,n3
9double precision f(n1,n2,n3),x(n1),y(n2),z(n3),t1,t2,t3
10double precision fx(n1,n2,n3),fy(n1,n2,n3),fxy(n1,n2,n3)   !, &
11!&                fz(n1,n2,n3),fxz(n1,n2,n3),fyz(n1,n2,n3), &
12!&                fxyz(n1,n2,n3)
13double precision, allocatable :: coefs(:)
14integer ii,jj,kk
15
16allocate (coefs(n1))
17do jj=1,n2
18do kk=1,n3
19  call spline(f(:,jj,kk),x,n1,t1,coefs)
20  call dspline(f(:,jj,kk),x,coefs,n1,fx(:,jj,kk))
21enddo
22enddo
23
24deallocate (coefs)
25allocate (coefs(n2))
26do ii=1,n1
27do kk=1,n3
28  call spline(f(ii,:,kk),y,n2,t2,coefs)
29  call dspline(f(ii,:,kk),y,coefs,n2,fy(ii,:,kk))
30  call spline(fx(ii,:,kk),y,n2,t2,coefs)
31  call dspline(fx(ii,:,kk),y,coefs,n2,fxy(ii,:,kk))
32enddo
33enddo
34
35!deallocate (coefs)
36!allocate (coefs(n3))
37!do ii=1,n1
38!do jj=1,n2
39!  call spline(f(ii,jj,:),z,n3,t3,coefs)
40!  call dspline(f(ii,jj,:),z,coefs,n3,fz(ii,jj,:))
41!  call spline(fx(ii,jj,:),z,n3,t3,coefs)
42!  call dspline(fx(ii,jj,:),z,coefs,n3,fxz(ii,jj,:))
43!  call spline(fy(ii,jj,:),z,n3,t3,coefs)
44!  call dspline(fy(ii,jj,:),z,coefs,n3,fyz(ii,jj,:))
45!  call spline(fxy(ii,jj,:),z,n3,t3,coefs)
46!  call dspline(fxy(ii,jj,:),z,coefs,n3,fxyz(ii,jj,:))
47!enddo
48!enddo
49
50return
51end subroutine trispline
52
53!***************************************************************************
54
55subroutine xyzspline(f,x,y,z,n1,n2,n3,t1,t2,t3,fx,fy,fz,fxy,fxz,fyz,fxyz)
56!     constructs interpolation coefficients for a periodic
57!     function f evaluated on a
58!     three dimensional grid (x,y,z) with size n1 x n2 x n3.
59!     t1,t2,t3 are the periods.
60implicit none
61integer n1,n2,n3
62double precision f(n1,n2,n3),x(n1),y(n2),z(n3),t1,t2,t3
63double precision fx(n1,n2,n3),fy(n1,n2,n3),fxy(n1,n2,n3), &
64&                fz(n1,n2,n3),fxz(n1,n2,n3),fyz(n1,n2,n3), &
65&                fxyz(n1,n2,n3)
66double precision, allocatable :: coefs(:)
67integer ii,jj,kk
68
69allocate (coefs(n1))
70do jj=1,n2
71do kk=1,n3
72  call spline(f(:,jj,kk),x,n1,t1,coefs)
73  call dspline(f(:,jj,kk),x,coefs,n1,fx(:,jj,kk))
74enddo
75enddo
76
77deallocate (coefs)
78allocate (coefs(n2))
79do ii=1,n1
80do kk=1,n3
81  call spline(f(ii,:,kk),y,n2,t2,coefs)
82  call dspline(f(ii,:,kk),y,coefs,n2,fy(ii,:,kk))
83  call spline(fx(ii,:,kk),y,n2,t2,coefs)
84  call dspline(fx(ii,:,kk),y,coefs,n2,fxy(ii,:,kk))
85enddo
86enddo
87
88deallocate (coefs)
89allocate (coefs(n3))
90do ii=1,n1
91do jj=1,n2
92  call spline(f(ii,jj,:),z,n3,t3,coefs)
93  call dspline(f(ii,jj,:),z,coefs,n3,fz(ii,jj,:))
94  call spline(fx(ii,jj,:),z,n3,t3,coefs)
95  call dspline(fx(ii,jj,:),z,coefs,n3,fxz(ii,jj,:))
96  call spline(fy(ii,jj,:),z,n3,t3,coefs)
97  call dspline(fy(ii,jj,:),z,coefs,n3,fyz(ii,jj,:))
98  call spline(fxy(ii,jj,:),z,n3,t3,coefs)
99  call dspline(fxy(ii,jj,:),z,coefs,n3,fxyz(ii,jj,:))
100enddo
101enddo
102
103return
104end subroutine xyzspline
105
106!***************************************************************************
107
108subroutine spl23(f,fx,fy,fxy,x,y,z,n1,n2,n3,t1,t2,t3,xv,yv,nz,zv,fv)
109implicit none
110integer :: n1,n2,n3,nz
111double precision :: f(n1,n2,n3),x(n1),y(n2),z(n3),t1,t2,t3,xv,yv,zv(nz)
112double precision :: fv(nz)
113double precision :: fx(n1,n2,n3),fy(n1,n2,n3),fxy(n1,n2,n3),zvals(n3)
114double precision, allocatable :: coefs(:)
115double precision :: veca(4),vecb(4),dmat(4,4)
116integer :: ii,jj,kk
117integer :: ihi,ilo,jhi,jlo
118double precision :: dx,dxg,dy,dyg
119
120if (xv.lt.x(1)) then
121  ihi=1
122  ilo=n1
123  dxg=x(ihi)-x(ilo)+t1
124  dx=(xv-x(ilo)+t1)/dxg
125else
126  ihi=n1
127  ilo=1
128  do
129    if (ihi-ilo.le.1) exit
130    ii=(ihi+ilo)/2
131    if (x(ii).gt.xv) then
132      ihi=ii
133    else
134      ilo=ii
135    endif
136  enddo
137  dxg=x(ihi)-x(ilo)
138  dx=(xv-x(ilo))/dxg
139endif
140
141if (yv.lt.y(1)) then
142  jhi=1
143  jlo=n2
144  dyg=y(jhi)-y(jlo)+t2
145  dy=(yv-y(jlo)+t2)/dyg
146else
147  jhi=n2
148  jlo=1
149  do
150    if (jhi-jlo.le.1) exit
151    jj=(jhi+jlo)/2
152    if (y(jj).gt.yv) then
153      jhi=jj
154    else
155      jlo=jj
156    endif
157  enddo
158  dyg=y(jhi)-y(jlo)
159  dy=(yv-y(jlo))/dyg
160endif
161
162!if (zv.lt.z(1)) then
163!  khi=1
164!  klo=n3
165!  dzg=z(khi)-z(klo)+t3
166!  dz=(zv-z(klo)+t3)/dzg
167!else
168!  khi=n3
169!  klo=1
170!  do
171!    if (khi-klo.le.1) exit
172!    kk=(khi+klo)/2
173!    if (z(kk).gt.zv) then
174!      khi=kk
175!    else
176!      klo=kk
177!    endif
178!  enddo
179!  dzg=z(khi)-z(klo)
180!  dz=(zv-z(klo))/dzg
181!endif
182
183do kk=1,n3
184!  allocate (coefs(n1))
185!  do jj=1,n2
186!    call spline(f(:,jj,kk),x,n1,t1,coefs)
187!    call dspline(f(:,jj,kk),x,coefs,n1,fx(:,jj,kk))
188!  enddo
189!  deallocate (coefs)
190!
191!  allocate (coefs(n2))
192!  do ii=1,n1
193!    call spline(f(ii,:,kk),y,n2,t2,coefs)
194!    call dspline(f(ii,:,kk),y,coefs,n2,fy(ii,:,kk))
195!    call spline(fx(ii,:,kk),y,n2,t2,coefs)
196!    call dspline(fx(ii,:,kk),y,coefs,n2,fxy(ii,:,kk))
197!  enddo
198!  deallocate (coefs)
199!
200  veca(1)=(1.d0-dx)**2*(1.d0+2.d0*dx)
201  veca(2)=dx**2*(3.d0-2.d0*dx)
202  veca(3)=(1.d0-dx)**2*dx*dxg
203  veca(4)=(dx-1.d0)*dx**2*dxg
204
205  vecb(1)=(1.d0-dy)**2*(1.d0+2.d0*dy)
206  vecb(2)=dy**2*(3.d0-2.d0*dy)
207  vecb(3)=(1.d0-dy)**2*dy*dyg
208  vecb(4)=(dy-1.d0)*dy**2*dyg
209
210  dmat(1,1)=f(ilo,jlo,kk)
211  dmat(1,2)=f(ilo,jhi,kk)
212  dmat(1,3)=fy(ilo,jlo,kk)
213  dmat(1,4)=fy(ilo,jhi,kk)
214  dmat(2,1)=f(ihi,jlo,kk)
215  dmat(2,2)=f(ihi,jhi,kk)
216  dmat(2,3)=fy(ihi,jlo,kk)
217  dmat(2,4)=fy(ihi,jhi,kk)
218  dmat(3,1)=fx(ilo,jlo,kk)
219  dmat(3,2)=fx(ilo,jhi,kk)
220  dmat(3,3)=fxy(ilo,jlo,kk)
221  dmat(3,4)=fxy(ilo,jhi,kk)
222  dmat(4,1)=fx(ihi,jlo,kk)
223  dmat(4,2)=fx(ihi,jhi,kk)
224  dmat(4,3)=fxy(ihi,jlo,kk)
225  dmat(4,4)=fxy(ihi,jhi,kk)
226
227  zvals(kk)=0.d0
228  do ii=1,4
229  do jj=1,4
230    zvals(kk)=zvals(kk)+dmat(jj,ii)*vecb(ii)*veca(jj)
231  enddo
232  enddo
233enddo
234
235allocate (coefs(n3))
236call spline(zvals,z,n3,t3,coefs)
237do kk=1,nz
238  call evalspline(zvals,z,t3,coefs,n3,zv(kk),fv(kk))
239enddo
240deallocate (coefs)
241
242return
243end subroutine spl23
244
245