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