[6] | 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 | |
---|