| 1 | c********************************************************************************* |
|---|
| 2 | c********************************************************************************* |
|---|
| 3 | c |
|---|
| 4 | c This file contains the subroutines transform3d, five, six and seven |
|---|
| 5 | c These routines must be linked to the Temperton FFT package to run |
|---|
| 6 | c |
|---|
| 7 | c******************************************************************************* |
|---|
| 8 | c******************************************************************************* |
|---|
| 9 | c Modified to operate on physical space data distributed in |
|---|
| 10 | c horizontal "slabs" across processors and wavenumber space |
|---|
| 11 | c data distributed in columns, i.e. each processor stores a |
|---|
| 12 | c finite z range of P. space data and a finite kx range of |
|---|
| 13 | c wavenumber space data. 2 data swapping routines perform the |
|---|
| 14 | c interprocessor exchanges: slabs_to_columns and columns_to_slabs. |
|---|
| 15 | c Except for the number of transforms executed as a block, these |
|---|
| 16 | c transform routines are applied as usual. |
|---|
| 17 | |
|---|
| 18 | c Each processor performs the following operations on a portion |
|---|
| 19 | c of the globally distributed data: |
|---|
| 20 | |
|---|
| 21 | c Forward Transforms: real data arranged in "slabs" |
|---|
| 22 | c FFT in x, FFT in y, swap data and construct column, S/C transform in z |
|---|
| 23 | c Inverse Transforms: complex data arranged in "columns" |
|---|
| 24 | c IC/S transform in z, swap data and construct slab, IFFT in y, IFFT in x |
|---|
| 25 | c***************************************************************************** |
|---|
| 26 | |
|---|
| 27 | subroutine transform3d (x, x_cplx, nx, ny, nz, iswitch, tswitch, itrig, |
|---|
| 28 | * trigx, trigy, trigz, wn, ifax, ifay, ifaz, |
|---|
| 29 | * num_dims,work,scratch,xt,myid,numprocs,locnx,locnz, |
|---|
| 30 | * comm,twoslice,subslice) |
|---|
| 31 | |
|---|
| 32 | C--------------------------------------------------------------------- |
|---|
| 33 | C |
|---|
| 34 | C transform3d.f - This FORTRAN subroutine accepts data |
|---|
| 35 | c stored in a 3-dimensional array and performs |
|---|
| 36 | c one of the the following transforms: |
|---|
| 37 | c |
|---|
| 38 | c iswitch=1 |
|---|
| 39 | c transform real data to the complex frequency domain |
|---|
| 40 | c Fourier transforms are applied in the x and y directions |
|---|
| 41 | c sine or cosine transforms are applied in the z direction |
|---|
| 42 | c tswitch='s' for sine transform, 'c' for cosine |
|---|
| 43 | c |
|---|
| 44 | c iswitch=-1 |
|---|
| 45 | c transform complex frequency domain data to real physical space data |
|---|
| 46 | c Inverse Fourier transforms are applied in the x and y directions |
|---|
| 47 | c sine or cosine transforms are applied in the z direction |
|---|
| 48 | c N.B. sine/cosine transforms are thei own inverses |
|---|
| 49 | c tswitch='s' for sine transform, 'c' for cosine |
|---|
| 50 | c |
|---|
| 51 | c This subroutine makes use of the Temperton FFT package |
|---|
| 52 | c cfft99.f along with the subroutines five.f, six.f and |
|---|
| 53 | c seven.f |
|---|
| 54 | c |
|---|
| 55 | c when num_dims=3, the minimum values for nx,ny and nz are(8,8,8) |
|---|
| 56 | c when num_dims=2, the minimum values for nx,ny and nz are(8,1,8) |
|---|
| 57 | c and ny=0 so that nyp1=ny+1=1 |
|---|
| 58 | |
|---|
| 59 | c |
|---|
| 60 | c code written by Kraig Winters June 30, 1997 |
|---|
| 61 | c this is a modified version of the Bryan Johns/KW code |
|---|
| 62 | c extended to allow the cases FFS/FFC/FFF in 2 or 3 dimensions |
|---|
| 63 | c extended for MPI parallel execution KW 8/18/00 |
|---|
| 64 | c |
|---|
| 65 | |
|---|
| 66 | c |
|---|
| 67 | c inputs: |
|---|
| 68 | c |
|---|
| 69 | c nx,ny,nz integers defining (base) array sizes |
|---|
| 70 | c num_dims integer specifying 2 (x-z) or 3 (xyz) dimensional problem (z="up") |
|---|
| 71 | c x real array of size (nx+2,ny+1,nz+1) for iswitch=1 or |
|---|
| 72 | c complex array of size (nx/2+1,ny+1,nz+1) for iswitch=-1 |
|---|
| 73 | c by convention first index ==> x direction, second ==> y, third ==> z |
|---|
| 74 | c iswitch integer specifying transform direction +1==> forward transform |
|---|
| 75 | c tswitch character*1 specifying symmetry of z transform 'c'/'s'/'f' |
|---|
| 76 | c for cosine/sine/fourier transforms in the vertical direction |
|---|
| 77 | c itrig integer flag specifying whether trig setup calculations need to be |
|---|
| 78 | c done. itrig=1 specifies tables need to be computed. |
|---|
| 79 | c N.B. if itrig=1, tables will be computed and then itrig set to 0. |
|---|
| 80 | c |
|---|
| 81 | c The following arrays must be dimensioned and passed but values are |
|---|
| 82 | c determined within transform3d.f when itrig=1. These arrays should not be |
|---|
| 83 | c overwritten if subsequent calls with itrig=0 are made to transform3d.f. |
|---|
| 84 | c |
|---|
| 85 | c trigx real array of size 3*(nx+2)/2+1 contains trig tables for x transforms |
|---|
| 86 | c trigy real array of size 2*(ny+1) contains trig tables for y transforms |
|---|
| 87 | c trigz real array of size 2*nz contains trig tables for z transforms |
|---|
| 88 | c wn complex array of size nz/4 + 1 auxilliary trig info for sine/cosine |
|---|
| 89 | c ifax integer array of size 13 factorization tables for FFTs |
|---|
| 90 | c ifay integer array of size 13 factorization tables for FFTs |
|---|
| 91 | c ifaz integer array of size 13 factorization tables for FFTs |
|---|
| 92 | c |
|---|
| 93 | c Scratch space: |
|---|
| 94 | c work real scratch array dimensioned in variable_declarations.h |
|---|
| 95 | c xt complex array of size (nx+2,nz/2+1) |
|---|
| 96 | c |
|---|
| 97 | c Output: |
|---|
| 98 | c |
|---|
| 99 | c x_cplx complex array of frequency data for iswitch=1 |
|---|
| 100 | c x real array of physical space data for iswitch=-1 |
|---|
| 101 | c |
|---|
| 102 | c N.B. Global real data "x" has size (nx+2,ny+1,nz+1) |
|---|
| 103 | c of which (nx,ny,nz+1) |
|---|
| 104 | c elements define the periodic sequences on the spatial grid. |
|---|
| 105 | c Complex data "x_cplx" has global size (nx/2+1,ny+1,nz+1) |
|---|
| 106 | c Each processor stores a portion of the input and |
|---|
| 107 | c output arrays. See real and complex data maps. |
|---|
| 108 | c |
|---|
| 109 | c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ |
|---|
| 110 | implicit none |
|---|
| 111 | integer i,j,k,kend,myid,numprocs,locnx,locnz |
|---|
| 112 | integer comm,twoslice,subslice,ierr |
|---|
| 113 | integer num_dims,plane,npz,nx,ny,nz,nxp2,nyp1 |
|---|
| 114 | integer ifax(13),ifay(13),ifaz(13),iswitch,isign,itrig,nyplanes |
|---|
| 115 | integer lotx,loty,lotz,incx,incy,incz,jumpx,jumpy,jumpz |
|---|
| 116 | real trigx(3*(nx+2)/2+1),trigy(2*(ny+1)),trigz(2*nz),pi,xnorm |
|---|
| 117 | real x(nx+2,ny+1,locnz+1),work(1) |
|---|
| 118 | complex arg ,wn(nz/4+1),xt(nx+2,nz/2+1) |
|---|
| 119 | complex x_cplx(locnx,ny+1,nz+1) |
|---|
| 120 | complex scratch(locnx,ny+1,nz+1) |
|---|
| 121 | character*1 tswitch |
|---|
| 122 | c |
|---|
| 123 | c |
|---|
| 124 | pi = 4.*atan(1.) |
|---|
| 125 | nxp2 = nx + 2 |
|---|
| 126 | nyp1 = ny + 1 ! note nyp1=0+1=1 for 2d problem |
|---|
| 127 | npz = nz/2 |
|---|
| 128 | |
|---|
| 129 | c |
|---|
| 130 | c Initialize trig arrays for the ffts. This code executes only if itrig=1. |
|---|
| 131 | c After execution, the results are passed back to the calling routine and |
|---|
| 132 | c remain usable if the arrays are not overwritten. itrig is then set to 0 |
|---|
| 133 | c so that trig initialization is only performed once. |
|---|
| 134 | c |
|---|
| 135 | if( itrig .eq. 1 ) then |
|---|
| 136 | call fftfax (nx,ifax,trigx) |
|---|
| 137 | if( num_dims .eq. 3 ) call cftfax (ny,ifay,trigy) |
|---|
| 138 | if( tswitch .eq. 's' .or. tswitch .eq. 'c' ) then |
|---|
| 139 | call cftfax (npz,ifaz,trigz) ! for sin/cos transforms via 5&6&7 |
|---|
| 140 | do i = 1, nz/4+1 |
|---|
| 141 | arg = pi * (i-1)/(nz/2.0) |
|---|
| 142 | wn(i) = cos(arg) + (0.0,1.0)*sin(arg) |
|---|
| 143 | enddo |
|---|
| 144 | elseif( tswitch .eq. 'f' .or. tswitch .eq. 'F' ) then |
|---|
| 145 | call cftfax (nz,ifaz,trigz) ! for complex vertical fourier transforms |
|---|
| 146 | endif |
|---|
| 147 | if( myid .eq. 0 ) then |
|---|
| 148 | write(6,*) 'trig tables set up ' |
|---|
| 149 | endif |
|---|
| 150 | itrig = 0 ! flag indicating that setup calculations have been done |
|---|
| 151 | endif |
|---|
| 152 | |
|---|
| 153 | if (iswitch .le. 0) go to 60 ! i.e. do inverse transform |
|---|
| 154 | |
|---|
| 155 | |
|---|
| 156 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 157 | c Forward transforms: (real data to frequency components) |
|---|
| 158 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 159 | |
|---|
| 160 | c Number of horizontal planes in each slab |
|---|
| 161 | if( myid .eq. numprocs-1 ) then |
|---|
| 162 | kend=locnz+1 |
|---|
| 163 | else |
|---|
| 164 | kend=locnz |
|---|
| 165 | endif |
|---|
| 166 | |
|---|
| 167 | if( num_dims .eq. 3 ) then |
|---|
| 168 | c For 3d problems do x then y transforms for each xy plane: |
|---|
| 169 | c |
|---|
| 170 | c ! real to complex x transforms |
|---|
| 171 | isign = -1 ! forward transform flag for Temperton package |
|---|
| 172 | incx = 1 ! words between adjacent elements in a transform |
|---|
| 173 | jumpx = nxp2 ! real words between successive first elements |
|---|
| 174 | lotx = ny ! number of x transforms to do |
|---|
| 175 | do k = 1,kend ! Fourier transform in x for each horizontal plane |
|---|
| 176 | call fft991 (x(1,1,k),work,trigx,ifax,incx,jumpx,nx,lotx,isign) |
|---|
| 177 | enddo |
|---|
| 178 | c |
|---|
| 179 | c complex to complex transforms in the y direction |
|---|
| 180 | c |
|---|
| 181 | incy = nxp2/2 ! complex words between successive elements |
|---|
| 182 | jumpy = 1 ! complex words between adjacent elements in a transform |
|---|
| 183 | loty = nxp2/2 ! number of complex y transforms to do |
|---|
| 184 | do k = 1,kend ! Fourier transform in y |
|---|
| 185 | call cfft99 (x(1,1,k),work,trigy,ifay,incy,jumpy,ny,loty,isign) |
|---|
| 186 | enddo |
|---|
| 187 | elseif( num_dims .eq. 2 ) then ! do all x transforms in a single call |
|---|
| 188 | c ! no y transforms required |
|---|
| 189 | c ! Do real to complex transforms in x direction |
|---|
| 190 | isign = -1 ! forward transform flag for Temperton package |
|---|
| 191 | incx = 1 ! words between adjacent elements in a transform |
|---|
| 192 | jumpx = nxp2 ! real words between successive first elements |
|---|
| 193 | lotx = kend ! number of x transforms to do |
|---|
| 194 | do j = 1,1 ! Fourier transform in x for each j(xz) plane |
|---|
| 195 | call fft991 (x(1,j,1),work,trigx,ifax,incx,jumpx,nx,lotx,isign) |
|---|
| 196 | enddo |
|---|
| 197 | |
|---|
| 198 | endif |
|---|
| 199 | |
|---|
| 200 | c |
|---|
| 201 | c z transforms in xz plane -- determine number of planes |
|---|
| 202 | c |
|---|
| 203 | if( num_dims .eq. 3 ) then |
|---|
| 204 | nyplanes=ny |
|---|
| 205 | elseif( num_dims .eq. 2) then |
|---|
| 206 | nyplanes=1 |
|---|
| 207 | endif |
|---|
| 208 | |
|---|
| 209 | |
|---|
| 210 | c ***** EXCHANGE DATA BETWEEN PROCESSORS SO THAT EACH *************** |
|---|
| 211 | c PROCESSOR HOLDS A "COLUMN" OF DATA RATHER THAN A "SLAB" |
|---|
| 212 | call slabs_to_columns(x,x_cplx,nx,ny,nz,locnx,locnz,numprocs, |
|---|
| 213 | * myid,comm,twoslice,subslice,work) |
|---|
| 214 | c************************************************************************* |
|---|
| 215 | |
|---|
| 216 | |
|---|
| 217 | c Do sine/cosine/fourier transforms in z direction, sine and cosine transforms |
|---|
| 218 | c are their own inverses and use isign=+1 by convention |
|---|
| 219 | c |
|---|
| 220 | |
|---|
| 221 | if ((tswitch .eq. 'C') .or. (tswitch .eq. 'c')) then |
|---|
| 222 | isign = +1 |
|---|
| 223 | incz = 2*locnx ! words between adjacent elements in a transform |
|---|
| 224 | jumpz = 1 ! real words between successive first elements |
|---|
| 225 | lotz = 2*locnx ! number of z transforms to do |
|---|
| 226 | do j = 1,nyplanes ! for each x-z plane of data |
|---|
| 227 | plane = j ! cosine transform in z |
|---|
| 228 | call six (x_cplx, 2*locnx, ny, nz, plane, wn, xt, work, |
|---|
| 229 | * trigz, ifaz, incz, jumpz, npz, lotz, isign) |
|---|
| 230 | enddo |
|---|
| 231 | elseif((tswitch .eq. 'S') .or. (tswitch .eq. 's')) then |
|---|
| 232 | isign = +1 |
|---|
| 233 | incz = 2*locnx ! words between adjacent elements in a transform |
|---|
| 234 | jumpz = 1 ! real words between successive first elements |
|---|
| 235 | lotz = 2*locnx ! number of z transforms to do |
|---|
| 236 | do j = 1, nyplanes |
|---|
| 237 | plane = j ! sine transform in z |
|---|
| 238 | call seven (x_cplx, 2*locnx, ny, nz, plane, wn, xt, |
|---|
| 239 | * work, trigz, ifaz, incz, jumpz, npz, lotz, isign) |
|---|
| 240 | enddo |
|---|
| 241 | elseif((tswitch .eq. 'F') .or. (tswitch .eq. 'f')) then ! Fourier transform |
|---|
| 242 | isign=-1 |
|---|
| 243 | incz = locnx*(ny+1) ! words between adjacent elements in a transform |
|---|
| 244 | jumpz = 1 ! words between successive first elements |
|---|
| 245 | lotz = locnx ! number of z transforms to do |
|---|
| 246 | do j = 1, nyplanes ! Complex Fourier transform in z |
|---|
| 247 | call cfft99 (x_cplx(1,j,1),work,trigz,ifaz,incz,jumpz,nz,lotz,isign) |
|---|
| 248 | enddo |
|---|
| 249 | endif |
|---|
| 250 | |
|---|
| 251 | c |
|---|
| 252 | c Normalize the results |
|---|
| 253 | c |
|---|
| 254 | if( num_dims .eq. 3 ) then |
|---|
| 255 | if( tswitch .eq. 's' ) xnorm = float(nz)/(2.*ny) |
|---|
| 256 | if( tswitch .eq. 'c' ) xnorm = float(nz)/(2.*ny) |
|---|
| 257 | if( tswitch .eq. 'f' ) xnorm = 1./float(ny*nz) |
|---|
| 258 | elseif( num_dims .eq. 2 ) then |
|---|
| 259 | if( tswitch .eq. 's' ) xnorm = float(nz)/2. |
|---|
| 260 | if( tswitch .eq. 'c' ) xnorm = float(nz)/2. |
|---|
| 261 | if( tswitch .eq. 'f' ) xnorm = 1./float(nz) |
|---|
| 262 | endif |
|---|
| 263 | |
|---|
| 264 | do k = 1, nz+1 |
|---|
| 265 | do j = 1,nyp1 ! nyp1=1 for 2d problem |
|---|
| 266 | do i = 1,locnx |
|---|
| 267 | x_cplx(i,j,k) = x_cplx(i,j,k)*xnorm |
|---|
| 268 | enddo |
|---|
| 269 | enddo |
|---|
| 270 | enddo |
|---|
| 271 | |
|---|
| 272 | return |
|---|
| 273 | |
|---|
| 274 | |
|---|
| 275 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 276 | c Perform the inverse transform (frequency components to real data) |
|---|
| 277 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 278 | 60 continue |
|---|
| 279 | |
|---|
| 280 | c inverse z transforms in xz plane -- determine number of planes |
|---|
| 281 | |
|---|
| 282 | if( num_dims .eq. 3 ) then |
|---|
| 283 | nyplanes=ny |
|---|
| 284 | elseif( num_dims .eq. 2) then |
|---|
| 285 | nyplanes=1 |
|---|
| 286 | endif |
|---|
| 287 | |
|---|
| 288 | c make a copy of the input data so it can be left unaltered |
|---|
| 289 | do k=1,nz+1 |
|---|
| 290 | do j=1,nyplanes |
|---|
| 291 | do i=1,locnx |
|---|
| 292 | scratch(i,j,k)=x_cplx(i,j,k) |
|---|
| 293 | enddo |
|---|
| 294 | enddo |
|---|
| 295 | enddo |
|---|
| 296 | |
|---|
| 297 | c |
|---|
| 298 | c Do inverse sin/cos/fourier transforms in z direction, |
|---|
| 299 | c sine and cosine transforms |
|---|
| 300 | c are their own inverses and use isign=+1 by convention |
|---|
| 301 | c |
|---|
| 302 | if ((tswitch .eq. 'C') .or. (tswitch .eq. 'c')) then |
|---|
| 303 | isign = +1 |
|---|
| 304 | incz = 2*locnx ! words between adjacent elements in a transform |
|---|
| 305 | jumpz = 1 ! real words between successive first elements |
|---|
| 306 | lotz = 2*locnx ! number of z transforms to do |
|---|
| 307 | do j = 1,nyplanes ! for each x-z plane of frequency data |
|---|
| 308 | plane = j ! inverse cosine transform in z |
|---|
| 309 | call six (scratch, 2*locnx, ny, nz, plane, wn, xt, work, |
|---|
| 310 | * trigz, ifaz, incz, jumpz, npz, lotz, isign) |
|---|
| 311 | enddo |
|---|
| 312 | elseif ((tswitch .eq. 'S') .or. (tswitch .eq. 's')) then |
|---|
| 313 | isign = +1 |
|---|
| 314 | incz = 2*locnx ! words between adjacent elements in a transform |
|---|
| 315 | jumpz = 1 ! real words between successive first elements |
|---|
| 316 | lotz = 2*locnx ! number of z transforms to do |
|---|
| 317 | do j = 1,nyplanes |
|---|
| 318 | plane = j ! inverse sine transform in Z |
|---|
| 319 | call seven (scratch, 2*locnx, ny, nz, plane, wn, xt, |
|---|
| 320 | * work, trigz, ifaz, incz, jumpz, npz, lotz, isign) |
|---|
| 321 | enddo |
|---|
| 322 | elseif ((tswitch .eq. 'F') .or. (tswitch .eq. 'f')) then |
|---|
| 323 | isign = +1 |
|---|
| 324 | incz = locnx*(ny+1) ! words between adjacent elements in a transform |
|---|
| 325 | jumpz = 1 ! real words between successive first elements |
|---|
| 326 | lotz = locnx ! number of z transforms to do |
|---|
| 327 | do j = 1, nyplanes ! Complex Fourier transform in z |
|---|
| 328 | call cfft99 (scratch(1,j,1),work,trigz,ifaz,incz,jumpz,nz,lotz,isign) |
|---|
| 329 | enddo |
|---|
| 330 | endif |
|---|
| 331 | |
|---|
| 332 | c ***** EXCHANGE DATA BETWEEN PROCESSORS SO THAT EACH *************** |
|---|
| 333 | c PROCESSOR HOLDS A "SLAB" OF DATA RATHER THAN A "COLUMN" |
|---|
| 334 | call columns_to_slabs(x,scratch,nx,ny,nz,locnx,locnz,numprocs, |
|---|
| 335 | * myid,comm,twoslice,subslice,work) |
|---|
| 336 | c************************************************************************* |
|---|
| 337 | |
|---|
| 338 | c Number of horizontal planes in each slab |
|---|
| 339 | if( myid .eq. numprocs-1 ) then |
|---|
| 340 | kend=locnz+1 |
|---|
| 341 | else |
|---|
| 342 | kend=locnz |
|---|
| 343 | endif |
|---|
| 344 | |
|---|
| 345 | c zero out the k nyquist positions prior to inverse transforming |
|---|
| 346 | do i=nx+1,nx+2 |
|---|
| 347 | do j=1,nyplanes |
|---|
| 348 | do k=1,kend |
|---|
| 349 | x(i,j,k)=0.0 |
|---|
| 350 | enddo |
|---|
| 351 | enddo |
|---|
| 352 | enddo |
|---|
| 353 | |
|---|
| 354 | c inverse Fourier transform in y if num_dims=3 |
|---|
| 355 | c |
|---|
| 356 | if( num_dims .eq. 3 ) then ! inverse transform in y then x |
|---|
| 357 | isign = +1 |
|---|
| 358 | incy = nxp2/2 ! complex words between successive first elements |
|---|
| 359 | jumpy = 1 ! complex words between adjacent elements in a transform |
|---|
| 360 | loty = nxp2/2 ! number of complex y transforms to do |
|---|
| 361 | do k = 1,kend ! inverse Fourier transform in y for each horiz. plane |
|---|
| 362 | call cfft99 (x(1,1,k),work,trigy,ifay,incy,jumpy,ny,loty,isign) |
|---|
| 363 | enddo |
|---|
| 364 | c |
|---|
| 365 | c inverse Fourier transform in x |
|---|
| 366 | c |
|---|
| 367 | isign=+1 |
|---|
| 368 | incx = 1 ! real words between adjacent elements in a transform |
|---|
| 369 | jumpx = nxp2 ! real words between successive first elements |
|---|
| 370 | lotx = ny ! number of x transforms to do |
|---|
| 371 | do k = 1,kend ! inverse Fourier transform in x for each horiz. plane |
|---|
| 372 | call fft991 (x(1,1,k),work,trigx,ifax,incx,jumpx,nx,lotx,isign) |
|---|
| 373 | enddo |
|---|
| 374 | elseif( num_dims .eq. 2 ) then ! do all x transforms in a single call |
|---|
| 375 | c ! no y transforms required |
|---|
| 376 | c ! Do real to complex transforms in x direction |
|---|
| 377 | isign = +1 ! inverse transform flag for Temperton package |
|---|
| 378 | incx = 1 ! words between adjacent elements in a transform |
|---|
| 379 | jumpx = nxp2 ! real words between successive first elements |
|---|
| 380 | lotx = kend ! number of x transforms to do |
|---|
| 381 | do j = 1,1 ! Fourier transform in x for each j(xz) plane |
|---|
| 382 | call fft991 (x(1,j,1),work,trigx,ifax,incx,jumpx,nx,lotx,isign) |
|---|
| 383 | enddo |
|---|
| 384 | endif |
|---|
| 385 | c |
|---|
| 386 | return |
|---|
| 387 | end |
|---|
| 388 | |
|---|
| 389 | |
|---|
| 390 | subroutine five (c, nsets, n, wn, work, trig, ifa, |
|---|
| 391 | * inc, jump, np, lot, isign) |
|---|
| 392 | |
|---|
| 393 | C----------------------------------------------------------------------------- |
|---|
| 394 | C |
|---|
| 395 | C five.f - This FORTRAN subroutine computes the inverse Fourier Transform |
|---|
| 396 | C of a number of sets of conjugate even data. Input consists of |
|---|
| 397 | C nsets vectors of n+1 complex elements of conjugate even data. |
|---|
| 398 | C Output consists of nsets vectors of n complex elements containing |
|---|
| 399 | C the results of the inverse transform. Since the input data is |
|---|
| 400 | C conjugate even, the results from the transform are real. Actually, |
|---|
| 401 | C these "real" results are stored in alternating real and imaginary |
|---|
| 402 | C components of the output vectors. |
|---|
| 403 | C |
|---|
| 404 | C Note that the inverse transform in this routine utilizes |
|---|
| 405 | C the Temperton FFT package. |
|---|
| 406 | C |
|---|
| 407 | C Procedure Five: as discussed in |
|---|
| 408 | C Cooley's paper in the Journal of Sound Vibration (1970), 12(3), |
|---|
| 409 | C pg. 315-337. |
|---|
| 410 | C |
|---|
| 411 | c code written by Kraig Winters July 1 1990 |
|---|
| 412 | c this is a modified version of Bryan Johns/KW code |
|---|
| 413 | C |
|---|
| 414 | C input: c (complex array of nsets by n+1) - contains the sets of |
|---|
| 415 | C input conjugate even data |
|---|
| 416 | C nsets (integer) - number of sets of data to be inverse transformed |
|---|
| 417 | C n (integer) - number of elements resulting from the inverse |
|---|
| 418 | C transform |
|---|
| 419 | C wn (complex array of n/2+1) - contains the complex exponential |
|---|
| 420 | C tables required by this routine |
|---|
| 421 | C work (real array) - work array required by the Temperton package |
|---|
| 422 | C trig (real array of np) - contains trig tables required by the |
|---|
| 423 | C Temperton package |
|---|
| 424 | C ifa (integer array of 13) - contains factorization tables |
|---|
| 425 | C required by the Temperton package |
|---|
| 426 | C inc (integer) - indicates direction of inverse transform |
|---|
| 427 | C jump (integer) - distance to next data sets in the input array |
|---|
| 428 | C np (integer) - number of data in each data set |
|---|
| 429 | C lot (integer) - number of data sets to be inverse transformed |
|---|
| 430 | C isign (integer) - indicates forward or inverse transform |
|---|
| 431 | C (for this routine, should be inverse: isign = +1) |
|---|
| 432 | C |
|---|
| 433 | C output: c (complex array of nsets by n+1) - contains the results |
|---|
| 434 | C of the inverse transform |
|---|
| 435 | C (NOTE: only nsets by n elements are used in output) |
|---|
| 436 | C------------------------------------------------------------------------------ |
|---|
| 437 | |
|---|
| 438 | integer nsets, n, np |
|---|
| 439 | complex c(nsets,n+1), wn(n/2+1), a1, a2 |
|---|
| 440 | real work(1), trig(np) |
|---|
| 441 | integer ifa(13), i, k |
|---|
| 442 | integer inc, jump, lot, isign, npt, nb2, nb2p1 |
|---|
| 443 | |
|---|
| 444 | npt = n |
|---|
| 445 | nb2 = npt/2 |
|---|
| 446 | nb2p1 = nb2+1 |
|---|
| 447 | |
|---|
| 448 | c |
|---|
| 449 | c Compute the input to the inverse Fourier transform "c" |
|---|
| 450 | c |
|---|
| 451 | |
|---|
| 452 | do i = 1, nsets |
|---|
| 453 | a1 = c(i,1) + conjg(c(i,npt+1)) |
|---|
| 454 | a2 = (c(i,1) - conjg(c(i,npt+1))) * wn(1) |
|---|
| 455 | c(i,1) = a1 + (0.0,1.0)*a2 |
|---|
| 456 | do k = 2, nb2 |
|---|
| 457 | a1 = c(i,k) + conjg(c(i,npt+2-k)) |
|---|
| 458 | a2 = (c(i,k) - conjg(c(i,npt+2-k))) * wn(k) |
|---|
| 459 | c(i,k) = a1 + (0.0,1.0)*a2 |
|---|
| 460 | c(i,npt+2-k) = conjg(a1 - (0.0,1.0)*a2) |
|---|
| 461 | enddo |
|---|
| 462 | a1 = c(i,nb2p1) + conjg(c(i,nb2p1)) |
|---|
| 463 | a2 = (c(i,nb2p1) - conjg(c(i,nb2p1))) * wn(nb2p1) |
|---|
| 464 | c(i,nb2p1) = a1 + (0.0,1.0)*a2 |
|---|
| 465 | enddo |
|---|
| 466 | |
|---|
| 467 | c |
|---|
| 468 | c Inverse Fourier transform the array "c" |
|---|
| 469 | c |
|---|
| 470 | call cfft99 (c, work, trig, ifa, inc, jump, np, lot, isign) |
|---|
| 471 | |
|---|
| 472 | return |
|---|
| 473 | end |
|---|
| 474 | |
|---|
| 475 | |
|---|
| 476 | |
|---|
| 477 | subroutine six (x, nxp2, ny, nz, plane, wn, xt, work, |
|---|
| 478 | * trigz, ifaz, incz, jumpz, npz, lotz, isign) |
|---|
| 479 | |
|---|
| 480 | C------------------------------------------------------------------------------ |
|---|
| 481 | C |
|---|
| 482 | C six.f - This FORTRAN subroutine computes the Fourier cosine transform |
|---|
| 483 | C in the z direction of a number of sets of data in a particular |
|---|
| 484 | C x-z plane of 3-D data. |
|---|
| 485 | C |
|---|
| 486 | C Input consists of nxp2 vectors of data of length nz+1 in the plane. |
|---|
| 487 | C Output are npx2 vectors of length nz+1 containing the coefficients of |
|---|
| 488 | C the cosine transform for each vector of data in the plane. The |
|---|
| 489 | C coefficients are overwritten onto the input data. |
|---|
| 490 | C |
|---|
| 491 | C This routine is an implementation of |
|---|
| 492 | C Procedure Six as described in Cooley's paper in the Journal of |
|---|
| 493 | C Sound Vibration, (1970) 12(3), pg. 315-337. |
|---|
| 494 | C This procedure allows a cosine transform to be performed via |
|---|
| 495 | C a Fourier transform routine but using only 1/4 of the input data. |
|---|
| 496 | C |
|---|
| 497 | C This routine utilizes the Temperton FFT package to perform the |
|---|
| 498 | C required Fourier transforms. |
|---|
| 499 | C |
|---|
| 500 | C N.B.: This routine must be linked to five.f and the Temperton FFT package. |
|---|
| 501 | C |
|---|
| 502 | C |
|---|
| 503 | c code written by Kraig Winters July 1 1990 |
|---|
| 504 | c this is a modified version of Bryan Johns/KW code |
|---|
| 505 | C |
|---|
| 506 | C input: x (real array of nxp2 by ny+1 by nz+1) - contains input data |
|---|
| 507 | C nxp2 (integer) - number of vectors in the plane to be transformed |
|---|
| 508 | C ny (integer) - number of x-z planes in x array |
|---|
| 509 | C nz (integer) - length of each vector in the z direction |
|---|
| 510 | C plane (integer) - indicates which x-z plane in x to transform |
|---|
| 511 | C wn (complex array of nz/4+1) - contains complex exponential table |
|---|
| 512 | C xt (complex array of nxp2 by nz/2+1) - work array containing |
|---|
| 513 | C x tilda terms required as input to Procedure five |
|---|
| 514 | C work (real array of nxp2*(ny+1)*nz) - work array required by |
|---|
| 515 | C Temperton FFT package |
|---|
| 516 | C trigz (real array of npz) - contains trig tables in z direction |
|---|
| 517 | C required by Temperton package |
|---|
| 518 | C ifaz (integer array of 13) - contains factorization tables |
|---|
| 519 | C required by Temperton package |
|---|
| 520 | C incz (integer) - indicates direction of cosine transform in |
|---|
| 521 | C the array x |
|---|
| 522 | C jumpz (integer) - distance between vectors in the plane |
|---|
| 523 | C npz (integer) - number of data in each vector |
|---|
| 524 | C lotz (integer) - number of vectors of data in the plane |
|---|
| 525 | C isign (integer) - indicates forward or inverse transform |
|---|
| 526 | C (should be inverse : isign = +1) |
|---|
| 527 | C |
|---|
| 528 | C output: x - contains cosine coefficients of input data |
|---|
| 529 | C |
|---|
| 530 | C------------------------------------------------------------------------------- |
|---|
| 531 | |
|---|
| 532 | integer nxp2, ny, nz, plane, incz, jumpz, npz, lotz, isign |
|---|
| 533 | integer ifaz(13) |
|---|
| 534 | real x(nxp2, ny+1, nz+1), work(1), trigz(npz) |
|---|
| 535 | real a1, a2, sum, pi |
|---|
| 536 | complex xt(nxp2,nz/2+1) |
|---|
| 537 | complex wn(nz/4+1) |
|---|
| 538 | integer i,k,n |
|---|
| 539 | |
|---|
| 540 | n = nz/2 |
|---|
| 541 | pi = 4.*atan(1.) |
|---|
| 542 | |
|---|
| 543 | C+ |
|---|
| 544 | C Compute a complex conjugate even sequence for each real vector of data |
|---|
| 545 | C (in Cooley's paper, these are the x tilda / N terms) |
|---|
| 546 | C (NOTE: we only need to form the first nz/2+1 of them because of symmetry) |
|---|
| 547 | C- |
|---|
| 548 | do i = 1, nxp2 |
|---|
| 549 | xt(i,1) = cmplx(x(i,plane,1)/nz, 0.0) |
|---|
| 550 | do k = 2, n |
|---|
| 551 | xt(i,k) = cmplx(x(i,plane,2*k-1)/nz, |
|---|
| 552 | * (x(i,plane,2*k-2) - x(i,plane,2*k))/nz) |
|---|
| 553 | enddo |
|---|
| 554 | xt(i,n+1) = cmplx(x(i,plane,nz+1)/nz,0.0) |
|---|
| 555 | enddo |
|---|
| 556 | C+ |
|---|
| 557 | C Inverse transform the conjugate even sequence into a real sequence using |
|---|
| 558 | C Procedure Five as outlined in Cooley's paper |
|---|
| 559 | C- |
|---|
| 560 | call five (xt, nxp2, n, wn, work, trigz, ifaz, incz, |
|---|
| 561 | * jumpz, npz, lotz, isign) |
|---|
| 562 | |
|---|
| 563 | C+ |
|---|
| 564 | C Post-process the results to form the cosine coefficients |
|---|
| 565 | C (see Cooley's Paper for the details here) |
|---|
| 566 | C- |
|---|
| 567 | do 100 i = 1, nxp2 |
|---|
| 568 | a1 = real(xt(i,1)) |
|---|
| 569 | sum = 0.0 |
|---|
| 570 | do 20 k = 2, nz, 2 |
|---|
| 571 | sum = sum + x(i,plane,k) |
|---|
| 572 | 20 continue |
|---|
| 573 | a2 = 2.0*sum/nz |
|---|
| 574 | x(i,plane,1) = a1 + a2 |
|---|
| 575 | x(i,plane,nz+1) = a1 - a2 |
|---|
| 576 | do 30 k = 2, nz-2, 2 |
|---|
| 577 | x(i,plane,k) = 0.5 * ((aimag(xt(i,k/2)) + |
|---|
| 578 | & aimag(xt(i,(nz+2-k)/2))) - |
|---|
| 579 | & (aimag(xt(i,k/2)) - aimag(xt(i,(nz+2-k)/2))) / |
|---|
| 580 | & (2.0*sin((k-1)*pi/nz))) |
|---|
| 581 | x(i,plane,k+1) = 0.5 * ((real(xt(i,k/2+1)) + |
|---|
| 582 | & real(xt(i,(nz+1-k)/2+1))) - |
|---|
| 583 | & (real(xt(i,k/2+1)) - real(xt(i,(nz+1-k)/2+1))) / |
|---|
| 584 | & (2.0*sin(k*pi/nz))) |
|---|
| 585 | 30 continue |
|---|
| 586 | x(i,plane,nz) = 0.5 * ((aimag(xt(i,nz/2)) + aimag(xt(i,1))) - |
|---|
| 587 | & (aimag(xt(i,nz/2)) - aimag(xt(i,1))) / |
|---|
| 588 | & (2.0*sin((nz-1)*pi/nz))) |
|---|
| 589 | 100 continue |
|---|
| 590 | |
|---|
| 591 | return |
|---|
| 592 | end |
|---|
| 593 | |
|---|
| 594 | |
|---|
| 595 | subroutine seven (x, nxp2, ny, nz, plane, wn, xt, work, |
|---|
| 596 | & trigz, ifaz, incz, jumpz, npz, lotz, isign) |
|---|
| 597 | |
|---|
| 598 | C------------------------------------------------------------------------------ |
|---|
| 599 | C |
|---|
| 600 | C seven.f - This FORTRAN subroutine computes the Fourier sine transform |
|---|
| 601 | C in the z direction of a number of sets of data in a particular |
|---|
| 602 | C x-z plane of 3-D data. |
|---|
| 603 | C |
|---|
| 604 | C Input consists of nxp2 vectors of data of length nz+1 in the plane. |
|---|
| 605 | C Output are npx2 vectors of length nz+1 containing the coefficients of |
|---|
| 606 | C the sine transform for each vector of data in the plane. The |
|---|
| 607 | C coefficients are overwritten onto the input data. |
|---|
| 608 | C |
|---|
| 609 | C This routine is an implementation of |
|---|
| 610 | C Procedure Seven as described in Cooley's paper in the Journal of |
|---|
| 611 | C Sound Vibration, (1970) 12(3), pg. 315-337. |
|---|
| 612 | C This procedure allows a sine transform to be performed via |
|---|
| 613 | C a Fourier transform routine but using only 1/4 of the input data. |
|---|
| 614 | C |
|---|
| 615 | C This routine employs the Temperton FFT package to perform the |
|---|
| 616 | C required Fourier transforms. |
|---|
| 617 | C |
|---|
| 618 | C N.B. This routine must be linked to five.f and the Temperton FFT package. |
|---|
| 619 | C |
|---|
| 620 | c code written by Kraig Winters July 1 1990 |
|---|
| 621 | c this is a modified version of Bryan Johns/KW code |
|---|
| 622 | |
|---|
| 623 | C |
|---|
| 624 | C input: x (real array of nxp2 by ny+1 by nz+1) - contains input data |
|---|
| 625 | C nxp2 (integer) - number of vectors in the plane to be transformed |
|---|
| 626 | C ny (integer) - number of x-z planes in x array |
|---|
| 627 | C nz (integer) - length of each vector in the z direction |
|---|
| 628 | C plane (integer) - indicates which x-z plane in x to transform |
|---|
| 629 | C wn (complex array of nz/4+1) - contains complex exponential table |
|---|
| 630 | C xt (complex array of nxp2 by nz/2+1) - work array containing |
|---|
| 631 | C x tilda terms required as input to Procedure five |
|---|
| 632 | C work (real array of nxp2*(ny+1)*nz) - work array required by |
|---|
| 633 | C Temperton FFT package |
|---|
| 634 | C trigz (real array of npz) - contains trig tables in z direction |
|---|
| 635 | C required by Temperton package |
|---|
| 636 | C ifaz (integer array of 13) - contains factorization tables |
|---|
| 637 | C required by Temperton package |
|---|
| 638 | C incz (integer) - indicates direction of sine transform in |
|---|
| 639 | C the array xt |
|---|
| 640 | C jumpz (integer) - distance between vectors in the plane |
|---|
| 641 | C npz (integer) - number of data in each vector |
|---|
| 642 | C lotz (integer) - number of vectors of data in the plane |
|---|
| 643 | C isign (integer) - indicates forward or inverse transform |
|---|
| 644 | C (should be inverse : isign = +1) |
|---|
| 645 | C |
|---|
| 646 | C output: x - contains sine coefficients of input data |
|---|
| 647 | C |
|---|
| 648 | C---------------------------------------------------------------------------- |
|---|
| 649 | |
|---|
| 650 | integer nxp2, ny, nz, plane, incz, jumpz, npz, lotz, isign |
|---|
| 651 | integer ifaz(13) |
|---|
| 652 | real x(nxp2, ny+1, nz+1), work(1), trigz(npz) |
|---|
| 653 | real pi |
|---|
| 654 | complex xt(nxp2,nz/2+1) |
|---|
| 655 | complex wn(nz/4+1) |
|---|
| 656 | integer i,k,n |
|---|
| 657 | |
|---|
| 658 | n = nz/2 |
|---|
| 659 | pi = 4.*atan(1.) |
|---|
| 660 | |
|---|
| 661 | C+ |
|---|
| 662 | C Compute a complex conjugate even sequence for each real vector of data |
|---|
| 663 | C (in Cooley's paper, these are the x tilda / N terms) |
|---|
| 664 | C (NOTE: we only need to form the first nz/2+1 of them because of symmetry) |
|---|
| 665 | C- |
|---|
| 666 | do i = 1, nxp2 |
|---|
| 667 | xt(i,1) = cmplx(-2.0*x(i,plane,2)/nz, 0.0) |
|---|
| 668 | do k = 2, n |
|---|
| 669 | xt(i,k) = |
|---|
| 670 | * cmplx ((x(i,plane,2*k-2) - x(i,plane,2*k))/nz, |
|---|
| 671 | * -x(i,plane,2*k-1)/nz) |
|---|
| 672 | enddo |
|---|
| 673 | xt(i,n+1) = cmplx(2.0*x(i,plane,nz)/nz,0.0) |
|---|
| 674 | enddo |
|---|
| 675 | |
|---|
| 676 | C+ |
|---|
| 677 | C Inverse transform the conjugate even sequence into a real sequence using |
|---|
| 678 | C Procedure Five as outlined in Cooley's paper |
|---|
| 679 | C- |
|---|
| 680 | call five (xt, nxp2, n, wn, work, trigz, ifaz, incz, |
|---|
| 681 | * jumpz, npz, lotz, isign) |
|---|
| 682 | |
|---|
| 683 | C+ |
|---|
| 684 | C Post-process the results to form the sine coefficients |
|---|
| 685 | C (see Cooley's Paper for the details here) |
|---|
| 686 | C- |
|---|
| 687 | do 100 i = 1, nxp2 |
|---|
| 688 | x(i,plane,1) = 0.0 |
|---|
| 689 | x(i,plane,nz+1) = 0.0 |
|---|
| 690 | do 30 k = 2, nz-2, 2 |
|---|
| 691 | x(i,plane,k) = 0.5 * ((aimag(xt(i,k/2)) - |
|---|
| 692 | & aimag(xt(i,(nz+2-k)/2))) - |
|---|
| 693 | & (aimag(xt(i,k/2)) + aimag(xt(i,(nz+2-k)/2))) / |
|---|
| 694 | & (2.0*sin((k-1)*pi/nz))) |
|---|
| 695 | x(i,plane,k+1) = 0.5 * ((real(xt(i,k/2+1)) - |
|---|
| 696 | & real(xt(i,(nz+1-k)/2+1))) - |
|---|
| 697 | & (real(xt(i,k/2+1)) + real(xt(i,(nz+1-k)/2+1))) / |
|---|
| 698 | & (2.0*sin(k*pi/nz))) |
|---|
| 699 | 30 continue |
|---|
| 700 | x(i,plane,nz) = 0.5 * ((aimag(xt(i,nz/2)) - aimag(xt(i,1))) - |
|---|
| 701 | & (aimag(xt(i,nz/2)) + aimag(xt(i,1))) / |
|---|
| 702 | & (2.0*sin((nz-1)*pi/nz))) |
|---|
| 703 | 100 continue |
|---|
| 704 | |
|---|
| 705 | return |
|---|
| 706 | end |
|---|
| 707 | |
|---|
