| 1 | SUBROUTINE intp3d ( x, y, z, f, kt, ft, nx, ny, nz, xt, yt, zt,d1, d2, d3 )
|
|---|
| 2 |
|
|---|
| 3 | implicit none
|
|---|
| 4 | !c
|
|---|
| 5 | !c---------------------------------------------------------------------
|
|---|
| 6 | !c
|
|---|
| 7 | !c purpose: interpolation of a function of three variables in an
|
|---|
| 8 | !c equidistant(!!!) table.
|
|---|
| 9 | !c
|
|---|
| 10 | !c method: 8-point Lagrange linear interpolation formula
|
|---|
| 11 | !c
|
|---|
| 12 | !c x input vector of first variable
|
|---|
| 13 | !c y input vector of second variable
|
|---|
| 14 | !c z input vector of third variable
|
|---|
| 15 | !c
|
|---|
| 16 | !c f output vector of interpolated function values
|
|---|
| 17 | !c
|
|---|
| 18 | !c kt vector length of input and output vectors
|
|---|
| 19 | !c
|
|---|
| 20 | !c ft 3d array of tabulated function values
|
|---|
| 21 | !c nx x-dimension of table
|
|---|
| 22 | !c ny y-dimension of table
|
|---|
| 23 | !c nz z-dimension of table
|
|---|
| 24 | !c xt vector of x-coordinates of table
|
|---|
| 25 | !c yt vector of y-coordinates of table
|
|---|
| 26 | !c zt vector of z-coordinates of table
|
|---|
| 27 | !c
|
|---|
| 28 | !c d1 centered derivative of ft with respect to x
|
|---|
| 29 | !c d2 centered derivative of ft with respect to y
|
|---|
| 30 | !c d3 centered derivative of ft with respect to z
|
|---|
| 31 | !c Note that d? only make sense when intp3d is called with kt=1
|
|---|
| 32 | !c---------------------------------------------------------------------
|
|---|
| 33 | !c
|
|---|
| 34 | !c
|
|---|
| 35 |
|
|---|
| 36 | integer :: kt,nx,ny,nz
|
|---|
| 37 | real :: x(kt),y(kt),z(kt),f(kt)
|
|---|
| 38 | real :: xt(nx),yt(ny),zt(nz)
|
|---|
| 39 | real :: ft(nx,ny,nz)
|
|---|
| 40 | real :: d1,d2,d3
|
|---|
| 41 |
|
|---|
| 42 | integer, PARAMETER :: ktx = 400
|
|---|
| 43 | real :: fh(ktx,8), delx(ktx), dely(ktx), delz(ktx), &
|
|---|
| 44 | & a1(ktx), a2(ktx), a3(ktx), a4(ktx), &
|
|---|
| 45 | & a5(ktx), a6(ktx), a7(ktx), a8(ktx)
|
|---|
| 46 |
|
|---|
| 47 | real :: dx,dy,dz,dxi,dyi,dzi,dxyi,dxzi,dyzi,dxyzi
|
|---|
| 48 | integer :: n,ix,iy,iz
|
|---|
| 49 |
|
|---|
| 50 |
|
|---|
| 51 | IF (kt .GT. ktx) STOP '***KTX**'
|
|---|
| 52 |
|
|---|
| 53 | !c------ determine spacing parameters of (equidistant!!!) table
|
|---|
| 54 |
|
|---|
| 55 | dx = (xt(nx) - xt(1)) / FLOAT(nx-1)
|
|---|
| 56 | dy = (yt(ny) - yt(1)) / FLOAT(ny-1)
|
|---|
| 57 | dz = (zt(nz) - zt(1)) / FLOAT(nz-1)
|
|---|
| 58 |
|
|---|
| 59 | dxi = 1. / dx
|
|---|
| 60 | dyi = 1. / dy
|
|---|
| 61 | dzi = 1. / dz
|
|---|
| 62 |
|
|---|
| 63 | dxyi = dxi * dyi
|
|---|
| 64 | dxzi = dxi * dzi
|
|---|
| 65 | dyzi = dyi * dzi
|
|---|
| 66 |
|
|---|
| 67 | dxyzi = dxi * dyi * dzi
|
|---|
| 68 |
|
|---|
| 69 |
|
|---|
| 70 | !c------- loop over all points to be interpolated
|
|---|
| 71 |
|
|---|
| 72 | DO n = 1, kt
|
|---|
| 73 |
|
|---|
| 74 | !c------- determine location in (equidistant!!!) table
|
|---|
| 75 |
|
|---|
| 76 | ix = 2 + INT( (x(n) - xt(1) - 1.e-10) * dxi )
|
|---|
| 77 | iy = 2 + INT( (y(n) - yt(1) - 1.e-10) * dyi )
|
|---|
| 78 | iz = 2 + INT( (z(n) - zt(1) - 1.e-10) * dzi )
|
|---|
| 79 |
|
|---|
| 80 | ix = MAX( 2, MIN( ix, nx ) )
|
|---|
| 81 | iy = MAX( 2, MIN( iy, ny ) )
|
|---|
| 82 | iz = MAX( 2, MIN( iz, nz ) )
|
|---|
| 83 |
|
|---|
| 84 | !c write(*,*) iy-1,iy,iy+1
|
|---|
| 85 |
|
|---|
| 86 | !c------- set-up auxiliary arrays for Lagrange interpolation
|
|---|
| 87 |
|
|---|
| 88 | delx(n) = xt(ix) - x(n)
|
|---|
| 89 | dely(n) = yt(iy) - y(n)
|
|---|
| 90 | delz(n) = zt(iz) - z(n)
|
|---|
| 91 |
|
|---|
| 92 | fh(n,1) = ft(ix , iy , iz )
|
|---|
| 93 | fh(n,2) = ft(ix-1, iy , iz )
|
|---|
| 94 | fh(n,3) = ft(ix , iy-1, iz )
|
|---|
| 95 | fh(n,4) = ft(ix , iy , iz-1)
|
|---|
| 96 | fh(n,5) = ft(ix-1, iy-1, iz )
|
|---|
| 97 | fh(n,6) = ft(ix-1, iy , iz-1)
|
|---|
| 98 | fh(n,7) = ft(ix , iy-1, iz-1)
|
|---|
| 99 | fh(n,8) = ft(ix-1, iy-1, iz-1)
|
|---|
| 100 |
|
|---|
| 101 | !c------ set up coefficients of the interpolation polynomial and
|
|---|
| 102 | !c evaluate function values
|
|---|
| 103 |
|
|---|
| 104 | a1(n) = fh(n,1)
|
|---|
| 105 | a2(n) = dxi * ( fh(n,2) - fh(n,1) )
|
|---|
| 106 | a3(n) = dyi * ( fh(n,3) - fh(n,1) )
|
|---|
| 107 | a4(n) = dzi * ( fh(n,4) - fh(n,1) )
|
|---|
| 108 | a5(n) = dxyi * ( fh(n,5) - fh(n,2) - fh(n,3) + fh(n,1) )
|
|---|
| 109 | a6(n) = dxzi * ( fh(n,6) - fh(n,2) - fh(n,4) + fh(n,1) )
|
|---|
| 110 | a7(n) = dyzi * ( fh(n,7) - fh(n,3) - fh(n,4) + fh(n,1) )
|
|---|
| 111 | a8(n) = dxyzi * ( fh(n,8) - fh(n,1) + fh(n,2) + fh(n,3) + &
|
|---|
| 112 | & fh(n,4) - fh(n,5) - fh(n,6) - fh(n,7) )
|
|---|
| 113 |
|
|---|
| 114 | d1 = -a2(n)
|
|---|
| 115 | d2 = -a3(n)
|
|---|
| 116 | d3 = -a4(n)
|
|---|
| 117 | f(n) = a1(n) + a2(n) * delx(n) &
|
|---|
| 118 | & + a3(n) * dely(n) &
|
|---|
| 119 | & + a4(n) * delz(n) &
|
|---|
| 120 | & + a5(n) * delx(n) * dely(n) &
|
|---|
| 121 | & + a6(n) * delx(n) * delz(n) &
|
|---|
| 122 | & + a7(n) * dely(n) * delz(n) &
|
|---|
| 123 | & + a8(n) * delx(n) * dely(n) * delz(n)
|
|---|
| 124 |
|
|---|
| 125 | ENDDO
|
|---|
| 126 |
|
|---|
| 127 | RETURN
|
|---|
| 128 | END SUBROUTINE intp3d
|
|---|
| 129 |
|
|---|
