Line data Source code
1 : ! FFT99F
2 : !
3 : ! PURPOSE PERFORMS MULTIPLE FAST FOURIER TRANSFORMS. THIS PACKAGE
4 : ! WILL PERFORM A NUMBER OF SIMULTANEOUS REAL/HALF-COMPLEX
5 : ! PERIODIC FOURIER TRANSFORMS OR CORRESPONDING INVERSE
6 : ! TRANSFORMS, I.E. GIVEN A SET OF REAL DATA VECTORS, THE
7 : ! PACKAGE RETURNS A SET OF 'HALF-COMPLEX' FOURIER
8 : ! COEFFICIENT VECTORS, OR VICE VERSA. THE LENGTH OF THE
9 : ! TRANSFORMS MUST BE AN EVEN NUMBER GREATER THAN 4 THAT HAS
10 : ! NO OTHER FACTORS EXCEPT POSSIBLY POWERS OF 2, 3, AND 5.
11 : ! THIS IS AN ALL FORTRAN VERSION OF THE CRAYLIB PACKAGE
12 : ! THAT IS MOSTLY WRITTEN IN CAL.
13 : !
14 : ! THE PACKAGE FFT99F CONTAINS SEVERAL USER-LEVEL ROUTINES:
15 : !
16 : ! SUBROUTINE SET99
17 : ! AN INITIALIZATION ROUTINE THAT MUST BE CALLED ONCE
18 : ! BEFORE A SEQUENCE OF CALLS TO THE FFT ROUTINES
19 : ! (PROVIDED THAT N IS NOT CHANGED).
20 : !
21 : ! SUBROUTINES FFT99 AND FFT991
22 : ! TWO FFT ROUTINES THAT RETURN SLIGHTLY DIFFERENT
23 : ! ARRANGEMENTS OF THE DATA IN GRIDPOINT SPACE.
24 : !
25 : !
26 : ! ACCESS THIS FORTRAN VERSION MAY BE ACCESSED WITH
27 : !
28 : ! *FORTRAN,P=XLIB,SN=FFT99F
29 : !
30 : ! TO ACCESS THE CRAY OBJECT CODE, CALLING THE USER ENTRY
31 : ! POINTS FROM A CRAY PROGRAM IS SUFFICIENT. THE SOURCE
32 : ! FORTRAN AND CAL CODE FOR THE CRAYLIB VERSION MAY BE
33 : ! ACCESSED USING
34 : !
35 : ! FETCH P=CRAYLIB,SN=FFT99
36 : ! FETCH P=CRAYLIB,SN=CAL99
37 : !
38 : ! USAGE LET N BE OF THE FORM 2**P * 3**Q * 5**R, WHERE P .GE. 1,
39 : ! Q .GE. 0, AND R .GE. 0. THEN A TYPICAL SEQUENCE OF
40 : ! CALLS TO TRANSFORM A GIVEN SET OF REAL VECTORS OF LENGTH
41 : ! N TO A SET OF 'HALF-COMPLEX' FOURIER COEFFICIENT VECTORS
42 : ! OF LENGTH N IS
43 : !
44 : ! DIMENSION IFAX(13),TRIGS(3*N/2+1),A(M*(N+2)),
45 : ! + WORK(M*(N+1))
46 : !
47 : ! CALL SET99 (TRIGS, IFAX, N)
48 : ! CALL FFT99 (A,WORK,TRIGS,IFAX,INC,JUMP,N,M,ISIGN)
49 : !
50 : ! SEE THE INDIVIDUAL WRITE-UPS FOR SET99, FFT99, AND
51 : ! FFT991 BELOW, FOR A DETAILED DESCRIPTION OF THE
52 : ! ARGUMENTS.
53 : !
54 : ! HISTORY THE PACKAGE WAS WRITTEN BY CLIVE TEMPERTON AT ECMWF IN
55 : ! NOVEMBER, 1978. IT WAS MODIFIED, DOCUMENTED, AND TESTED
56 : ! FOR NCAR BY RUSS REW IN SEPTEMBER, 1980.
57 : !
58 : !-----------------------------------------------------------------------
59 : !
60 : ! SUBROUTINE SET99 (TRIGS, IFAX, N)
61 : !
62 : ! PURPOSE A SET-UP ROUTINE FOR FFT99 AND FFT991. IT NEED ONLY BE
63 : ! CALLED ONCE BEFORE A SEQUENCE OF CALLS TO THE FFT
64 : ! ROUTINES (PROVIDED THAT N IS NOT CHANGED).
65 : !
66 : ! ARGUMENT IFAX(13),TRIGS(3*N/2+1)
67 : ! DIMENSIONS
68 : !
69 : ! ARGUMENTS
70 : !
71 : ! ON INPUT TRIGS
72 : ! A FLOATING POINT ARRAY OF DIMENSION 3*N/2 IF N/2 IS
73 : ! EVEN, OR 3*N/2+1 IF N/2 IS ODD.
74 : !
75 : ! IFAX
76 : ! AN INTEGER ARRAY. THE NUMBER OF ELEMENTS ACTUALLY USED
77 : ! WILL DEPEND ON THE FACTORIZATION OF N. DIMENSIONING
78 : ! IFAX FOR 13 SUFFICES FOR ALL N LESS THAN A MILLION.
79 : !
80 : ! N
81 : ! AN EVEN NUMBER GREATER THAN 4 THAT HAS NO PRIME FACTOR
82 : ! GREATER THAN 5. N IS THE LENGTH OF THE TRANSFORMS (SEE
83 : ! THE DOCUMENTATION FOR FFT99 AND FFT991 FOR THE
84 : ! DEFINITIONS OF THE TRANSFORMS).
85 : !
86 : ! ON OUTPUT IFAX
87 : ! CONTAINS THE FACTORIZATION OF N/2. IFAX(1) IS THE
88 : ! NUMBER OF FACTORS, AND THE FACTORS THEMSELVES ARE STORED
89 : ! IN IFAX(2),IFAX(3),... IF SET99 IS CALLED WITH N ODD,
90 : ! OR IF N HAS ANY PRIME FACTORS GREATER THAN 5, IFAX(1)
91 : ! IS SET TO -99.
92 : !
93 : ! TRIGS
94 : ! AN ARRAY OF TRIGONOMETRIC FUNCTION VALUES SUBSEQUENTLY
95 : ! USED BY THE FFT ROUTINES.
96 : !
97 : !-----------------------------------------------------------------------
98 : !
99 : ! SUBROUTINE FFT991 (A,WORK,TRIGS,IFAX,INC,JUMP,N,M,ISIGN)
100 : ! AND
101 : ! SUBROUTINE FFT99 (A,WORK,TRIGS,IFAX,INC,JUMP,N,M,ISIGN)
102 : !
103 : ! PURPOSE PERFORM A NUMBER OF SIMULTANEOUS REAL/HALF-COMPLEX
104 : ! PERIODIC FOURIER TRANSFORMS OR CORRESPONDING INVERSE
105 : ! TRANSFORMS, USING ORDINARY SPATIAL ORDER OF GRIDPOINT
106 : ! VALUES (FFT991) OR EXPLICIT CYCLIC CONTINUITY IN THE
107 : ! GRIDPOINT VALUES (FFT99). GIVEN A SET
108 : ! OF REAL DATA VECTORS, THE PACKAGE RETURNS A SET OF
109 : ! 'HALF-COMPLEX' FOURIER COEFFICIENT VECTORS, OR VICE
110 : ! VERSA. THE LENGTH OF THE TRANSFORMS MUST BE AN EVEN
111 : ! NUMBER THAT HAS NO OTHER FACTORS EXCEPT POSSIBLY POWERS
112 : ! OF 2, 3, AND 5. THESE VERSION OF FFT991 AND FFT99 ARE
113 : ! OPTIMIZED FOR USE ON THE CRAY-1.
114 : !
115 : ! ARGUMENT A(M*(N+2)), WORK(M*(N+1)), TRIGS(3*N/2+1), IFAX(13)
116 : ! DIMENSIONS
117 : !
118 : ! ARGUMENTS
119 : !
120 : ! ON INPUT A
121 : ! AN ARRAY OF LENGTH M*(N+2) CONTAINING THE INPUT DATA
122 : ! OR COEFFICIENT VECTORS. THIS ARRAY IS OVERWRITTEN BY
123 : ! THE RESULTS.
124 : !
125 : ! WORK
126 : ! A WORK ARRAY OF DIMENSION M*(N+1)
127 : !
128 : ! TRIGS
129 : ! AN ARRAY SET UP BY SET99, WHICH MUST BE CALLED FIRST.
130 : !
131 : ! IFAX
132 : ! AN ARRAY SET UP BY SET99, WHICH MUST BE CALLED FIRST.
133 : !
134 : ! INC
135 : ! THE INCREMENT (IN WORDS) BETWEEN SUCCESSIVE ELEMENTS OF
136 : ! EACH DATA OR COEFFICIENT VECTOR (E.G. INC=1 FOR
137 : ! CONSECUTIVELY STORED DATA).
138 : !
139 : ! JUMP
140 : ! THE INCREMENT (IN WORDS) BETWEEN THE FIRST ELEMENTS OF
141 : ! SUCCESSIVE DATA OR COEFFICIENT VECTORS. ON THE CRAY-1,
142 : ! TRY TO ARRANGE DATA SO THAT JUMP IS NOT A MULTIPLE OF 8
143 : ! (TO AVOID MEMORY BANK CONFLICTS). FOR CLARIFICATION OF
144 : ! INC AND JUMP, SEE THE EXAMPLES BELOW.
145 : !
146 : ! N
147 : ! THE LENGTH OF EACH TRANSFORM (SEE DEFINITION OF
148 : ! TRANSFORMS, BELOW).
149 : !
150 : ! M
151 : ! THE NUMBER OF TRANSFORMS TO BE DONE SIMULTANEOUSLY.
152 : !
153 : ! ISIGN
154 : ! = +1 FOR A TRANSFORM FROM FOURIER COEFFICIENTS TO
155 : ! GRIDPOINT VALUES.
156 : ! = -1 FOR A TRANSFORM FROM GRIDPOINT VALUES TO FOURIER
157 : ! COEFFICIENTS.
158 : !
159 : ! ON OUTPUT A
160 : ! IF ISIGN = +1, AND M COEFFICIENT VECTORS ARE SUPPLIED
161 : ! EACH CONTAINING THE SEQUENCE:
162 : !
163 : ! A(0),B(0),A(1),B(1),...,A(N/2),B(N/2) (N+2 VALUES)
164 : !
165 : ! THEN THE RESULT CONSISTS OF M DATA VECTORS EACH
166 : ! CONTAINING THE CORRESPONDING N+2 GRIDPOINT VALUES:
167 : !
168 : ! FOR FFT991, X(0), X(1), X(2),...,X(N-1),0,0.
169 : ! FOR FFT99, X(N-1),X(0),X(1),X(2),...,X(N-1),X(0).
170 : ! (EXPLICIT CYCLIC CONTINUITY)
171 : !
172 : ! WHEN ISIGN = +1, THE TRANSFORM IS DEFINED BY:
173 : ! X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N))
174 : ! WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K)
175 : ! AND I=SQRT (-1)
176 : !
177 : ! IF ISIGN = -1, AND M DATA VECTORS ARE SUPPLIED EACH
178 : ! CONTAINING A SEQUENCE OF GRIDPOINT VALUES X(J) AS
179 : ! DEFINED ABOVE, THEN THE RESULT CONSISTS OF M VECTORS
180 : ! EACH CONTAINING THE CORRESPONDING FOURIER COFFICIENTS
181 : ! A(K), B(K), 0 .LE. K .LE N/2.
182 : !
183 : ! WHEN ISIGN = -1, THE INVERSE TRANSFORM IS DEFINED BY:
184 : ! C(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*EXP(-2*I*J*K*PI/N))
185 : ! WHERE C(K)=A(K)+I*B(K) AND I=SQRT(-1)
186 : !
187 : ! A CALL WITH ISIGN=+1 FOLLOWED BY A CALL WITH ISIGN=-1
188 : ! (OR VICE VERSA) RETURNS THE ORIGINAL DATA.
189 : !
190 : ! NOTE: THE FACT THAT THE GRIDPOINT VALUES X(J) ARE REAL
191 : ! IMPLIES THAT B(0)=B(N/2)=0. FOR A CALL WITH ISIGN=+1,
192 : ! IT IS NOT ACTUALLY NECESSARY TO SUPPLY THESE ZEROS.
193 : !
194 : ! EXAMPLES GIVEN 19 DATA VECTORS EACH OF LENGTH 64 (+2 FOR EXPLICIT
195 : ! CYCLIC CONTINUITY), COMPUTE THE CORRESPONDING VECTORS OF
196 : ! FOURIER COEFFICIENTS. THE DATA MAY, FOR EXAMPLE, BE
197 : ! ARRANGED LIKE THIS:
198 : !
199 : ! FIRST DATA A(1)= . . . A(66)= A(70)
200 : ! VECTOR X(63) X(0) X(1) X(2) ... X(63) X(0) (4 EMPTY LOCATIONS)
201 : !
202 : ! SECOND DATA A(71)= . . . A(140)
203 : ! VECTOR X(63) X(0) X(1) X(2) ... X(63) X(0) (4 EMPTY LOCATIONS)
204 : !
205 : ! AND SO ON. HERE INC=1, JUMP=70, N=64, M=19, ISIGN=-1,
206 : ! AND FFT99 SHOULD BE USED (BECAUSE OF THE EXPLICIT CYCLIC
207 : ! CONTINUITY).
208 : !
209 : ! ALTERNATIVELY THE DATA MAY BE ARRANGED LIKE THIS:
210 : !
211 : ! FIRST SECOND LAST
212 : ! DATA DATA DATA
213 : ! VECTOR VECTOR VECTOR
214 : !
215 : ! A(1)= A(2)= A(19)=
216 : !
217 : ! X(63) X(63) . . . X(63)
218 : ! A(20)= X(0) X(0) . . . X(0)
219 : ! A(39)= X(1) X(1) . . . X(1)
220 : ! . . .
221 : ! . . .
222 : ! . . .
223 : !
224 : ! IN WHICH CASE WE HAVE INC=19, JUMP=1, AND THE REMAINING
225 : ! PARAMETERS ARE THE SAME AS BEFORE. IN EITHER CASE, EACH
226 : ! COEFFICIENT VECTOR OVERWRITES THE CORRESPONDING INPUT
227 : ! DATA VECTOR.
228 : !-----------------------------------------------------------------------
229 : !
230 : ! $Id$
231 : ! $Author$
232 :
233 : !================================================================================================
234 :
235 0 : SUBROUTINE FFT99(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN)
236 : !
237 : !-----------------------------------------------------------------------
238 : ! SUBROUTINE "FFT99" - MULTIPLE FAST REAL PERIODIC TRANSFORM
239 : ! CORRESPONDING TO OLD SCALAR ROUTINE FFT9
240 : ! PROCEDURE USED TO CONVERT TO HALF-LENGTH COMPLEX TRANSFORM
241 : ! IS GIVEN BY COOLEY, LEWIS AND WELCH (J. SOUND VIB., VOL. 12
242 : ! (1970), 315-337)
243 : !
244 : ! A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA
245 : ! WORK IS AN AREA OF SIZE (N+1)*LOT
246 : ! TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES
247 : ! IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N/2
248 : ! INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR'
249 : ! (E.G. INC=1 FOR CONSECUTIVELY STORED DATA)
250 : ! JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR
251 : ! N IS THE LENGTH OF THE DATA VECTORS
252 : ! LOT IS THE NUMBER OF DATA VECTORS
253 : ! ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT
254 : ! = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL
255 : !
256 : ! ORDERING OF COEFFICIENTS:
257 : ! A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2)
258 : ! WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED
259 : !
260 : ! ORDERING OF DATA:
261 : ! X(N-1),X(0),X(1),X(2),...,X(N),X(0)
262 : ! I.E. EXPLICIT CYCLIC CONTINUITY; (N+2) LOCATIONS REQUIRED
263 : !
264 : ! VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN
265 : ! PARALLEL
266 : !
267 : ! *** N.B. N IS ASSUMED TO BE AN EVEN NUMBER
268 : !
269 : ! DEFINITION OF TRANSFORMS:
270 : ! -------------------------
271 : !
272 : ! ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N))
273 : ! WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K)
274 : !
275 : ! ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N))
276 : ! B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N))
277 : !
278 : !-----------------------------------------------------------------------
279 : !
280 : use shr_kind_mod, only: r8 => shr_kind_r8
281 : implicit none
282 : !
283 : !------------------------------Arguments--------------------------------
284 : !
285 : integer IFAX(13), inc, jump, n, lot, isign
286 : real(r8) A(LOT* (N+2) ), WORK(LOT*(N+1)), TRIGS(3*N/2+1)
287 : !
288 : !---------------------------Local variables-----------------------------
289 : !
290 : integer i, j, k, l, m, ia, ib, la, ink, nh, nx, nfax
291 : integer ibase, jbase, igo
292 : !
293 : !-----------------------------------------------------------------------
294 : !
295 0 : NFAX=IFAX(1)
296 0 : NX=N+1
297 0 : NH=N/2
298 0 : INK=INC+INC
299 0 : IF (ISIGN.EQ.+1) GO TO 30
300 : !
301 : ! IF NECESSARY, TRANSFER DATA TO WORK AREA
302 0 : IGO=50
303 0 : IF (MOD(NFAX,2).EQ.1) GOTO 40
304 0 : IBASE=INC+1
305 0 : JBASE=1
306 0 : DO 20 L=1,LOT
307 0 : I=IBASE
308 0 : J=JBASE
309 0 : DO 10 M=1,N
310 0 : WORK(J)=A(I)
311 0 : I=I+INC
312 0 : J=J+1
313 0 : 10 CONTINUE
314 0 : IBASE=IBASE+JUMP
315 0 : JBASE=JBASE+NX
316 0 : 20 CONTINUE
317 : !
318 : IGO=60
319 0 : GO TO 40
320 : !
321 : ! PREPROCESSING (ISIGN=+1)
322 : ! ------------------------
323 : !
324 : 30 CONTINUE
325 0 : CALL FFT99A(A,WORK,TRIGS,INC,JUMP,N,LOT)
326 0 : IGO=60
327 : !
328 : ! COMPLEX TRANSFORM
329 : ! -----------------
330 : !
331 : 40 CONTINUE
332 0 : IA=INC+1
333 0 : LA=1
334 0 : DO 80 K=1,NFAX
335 0 : IF (IGO.EQ.60) GO TO 60
336 : 50 CONTINUE
337 0 : CALL VPASSM(A(IA),A(IA+INC),WORK(1),WORK(2),TRIGS, &
338 0 : INK,2,JUMP,NX,LOT,NH,IFAX(K+1),LA)
339 0 : IGO=60
340 0 : GO TO 70
341 : 60 CONTINUE
342 0 : CALL VPASSM(WORK(1),WORK(2),A(IA),A(IA+INC),TRIGS, &
343 0 : 2,INK,NX,JUMP,LOT,NH,IFAX(K+1),LA)
344 0 : IGO=50
345 : 70 CONTINUE
346 0 : LA=LA*IFAX(K+1)
347 0 : 80 CONTINUE
348 : !
349 0 : IF (ISIGN.EQ.-1) GO TO 130
350 : !
351 : ! IF NECESSARY, TRANSFER DATA FROM WORK AREA
352 0 : IF (MOD(NFAX,2).EQ.1) GO TO 110
353 0 : IBASE=1
354 0 : JBASE=IA
355 0 : DO 100 L=1,LOT
356 0 : I=IBASE
357 0 : J=JBASE
358 0 : DO 90 M=1,N
359 0 : A(J)=WORK(I)
360 0 : I=I+1
361 0 : J=J+INC
362 0 : 90 CONTINUE
363 0 : IBASE=IBASE+NX
364 0 : JBASE=JBASE+JUMP
365 0 : 100 CONTINUE
366 : !
367 : ! FILL IN CYCLIC BOUNDARY POINTS
368 : 110 CONTINUE
369 0 : IA=1
370 0 : IB=N*INC+1
371 0 : DO 120 L=1,LOT
372 0 : A(IA)=A(IB)
373 0 : A(IB+INC)=A(IA+INC)
374 0 : IA=IA+JUMP
375 0 : IB=IB+JUMP
376 0 : 120 CONTINUE
377 0 : GO TO 140
378 : !
379 : ! POSTPROCESSING (ISIGN=-1):
380 : ! --------------------------
381 : !
382 : 130 CONTINUE
383 0 : CALL FFT99B(WORK,A,TRIGS,INC,JUMP,N,LOT)
384 : !
385 : 140 CONTINUE
386 0 : RETURN
387 : END SUBROUTINE FFT99
388 :
389 : !================================================================================================
390 :
391 3625440 : SUBROUTINE FFT99A(A,WORK,TRIGS,INC,JUMP,N,LOT)
392 : !
393 : !-----------------------------------------------------------------------
394 : ! SUBROUTINE FFT99A - PREPROCESSING STEP FOR FFT99, ISIGN=+1
395 : ! (SPECTRAL TO GRIDPOINT TRANSFORM)
396 : !-----------------------------------------------------------------------
397 : !
398 : use shr_kind_mod, only: r8 => shr_kind_r8
399 : implicit none
400 : !
401 : !------------------------------Arguments--------------------------------
402 : !
403 : integer inc, jump, n, lot
404 : real(r8) A(*), WORK(*), TRIGS(*)
405 : !
406 : !---------------------------Local variables-----------------------------
407 : !
408 : integer iabase, ibbase, jabase, jbbase, ia, ib, ink
409 : integer ja, jb, k, l, nh, nx
410 : real(r8) c, s
411 : !
412 : !-----------------------------------------------------------------------
413 : !
414 3625440 : NH=N/2
415 3625440 : NX=N+1
416 3625440 : INK=INC+INC
417 : !
418 : ! A(0) AND A(N/2)
419 3625440 : IA=1
420 3625440 : IB=N*INC+1
421 3625440 : JA=1
422 3625440 : JB=2
423 17843280 : DO 10 L=1,LOT
424 14217840 : WORK(JA)=A(IA)+A(IB)
425 14217840 : WORK(JB)=A(IA)-A(IB)
426 14217840 : IA=IA+JUMP
427 14217840 : IB=IB+JUMP
428 14217840 : JA=JA+NX
429 14217840 : JB=JB+NX
430 3625440 : 10 CONTINUE
431 : !
432 : ! REMAINING WAVENUMBERS
433 3625440 : IABASE=2*INC+1
434 3625440 : IBBASE=(N-2)*INC+1
435 3625440 : JABASE=3
436 3625440 : JBBASE=N-1
437 : !
438 3625440 : DO 30 K=3,NH,2
439 126890400 : IA=IABASE
440 126890400 : IB=IBBASE
441 126890400 : JA=JABASE
442 126890400 : JB=JBBASE
443 126890400 : C=TRIGS(N+K)
444 126890400 : S=TRIGS(N+K+1)
445 624514800 : DO 20 L=1,LOT
446 0 : WORK(JA)=(A(IA)+A(IB))- &
447 497624400 : (S*(A(IA)-A(IB))+C*(A(IA+INC)+A(IB+INC)))
448 0 : WORK(JB)=(A(IA)+A(IB))+ &
449 497624400 : (S*(A(IA)-A(IB))+C*(A(IA+INC)+A(IB+INC)))
450 0 : WORK(JA+1)=(C*(A(IA)-A(IB))-S*(A(IA+INC)+A(IB+INC)))+ &
451 497624400 : (A(IA+INC)-A(IB+INC))
452 0 : WORK(JB+1)=(C*(A(IA)-A(IB))-S*(A(IA+INC)+A(IB+INC)))- &
453 497624400 : (A(IA+INC)-A(IB+INC))
454 497624400 : IA=IA+JUMP
455 497624400 : IB=IB+JUMP
456 497624400 : JA=JA+NX
457 497624400 : JB=JB+NX
458 126890400 : 20 CONTINUE
459 126890400 : IABASE=IABASE+INK
460 126890400 : IBBASE=IBBASE-INK
461 126890400 : JABASE=JABASE+2
462 126890400 : JBBASE=JBBASE-2
463 3625440 : 30 CONTINUE
464 : !
465 3625440 : IF (IABASE.NE.IBBASE) GO TO 50
466 : ! WAVENUMBER N/4 (IF IT EXISTS)
467 3625440 : IA=IABASE
468 3625440 : JA=JABASE
469 17843280 : DO 40 L=1,LOT
470 14217840 : WORK(JA)=2.0_r8*A(IA)
471 14217840 : WORK(JA+1)=-2.0_r8*A(IA+INC)
472 14217840 : IA=IA+JUMP
473 14217840 : JA=JA+NX
474 3625440 : 40 CONTINUE
475 : !
476 : 50 CONTINUE
477 3625440 : RETURN
478 : END SUBROUTINE FFT99A
479 :
480 : !================================================================================================
481 :
482 3625440 : SUBROUTINE FFT99B(WORK,A,TRIGS,INC,JUMP,N,LOT)
483 : !
484 : !-----------------------------------------------------------------------
485 : ! SUBROUTINE FFT99B - POSTPROCESSING STEP FOR FFT99, ISIGN=-1
486 : ! (GRIDPOINT TO SPECTRAL TRANSFORM)
487 : !-----------------------------------------------------------------------
488 : !
489 : use shr_kind_mod, only: r8 => shr_kind_r8
490 : implicit none
491 : !
492 : !------------------------------Arguments--------------------------------
493 : !
494 : integer inc, jump, n, lot
495 : real(r8) WORK(*), A(*), TRIGS(*)
496 : !
497 : !---------------------------Local variables-----------------------------
498 : !
499 : integer iabase, ibbase, jabase, jbbase, ia, ib, ink
500 : integer ja, jb, k, l, nh, nx
501 : real(r8) scale, s, c
502 : !
503 : !-----------------------------------------------------------------------
504 : !
505 3625440 : NH=N/2
506 3625440 : NX=N+1
507 3625440 : INK=INC+INC
508 : !
509 : ! A(0) AND A(N/2)
510 3625440 : SCALE=1.0_r8/real(N,r8)
511 3625440 : IA=1
512 3625440 : IB=2
513 3625440 : JA=1
514 3625440 : JB=N*INC+1
515 17843280 : DO 10 L=1,LOT
516 14217840 : A(JA)=SCALE*(WORK(IA)+WORK(IB))
517 14217840 : A(JB)=SCALE*(WORK(IA)-WORK(IB))
518 14217840 : A(JA+INC)=0.0_r8
519 14217840 : A(JB+INC)=0.0_r8
520 14217840 : IA=IA+NX
521 14217840 : IB=IB+NX
522 14217840 : JA=JA+JUMP
523 14217840 : JB=JB+JUMP
524 3625440 : 10 CONTINUE
525 : !
526 : ! REMAINING WAVENUMBERS
527 3625440 : SCALE=0.5_r8*SCALE
528 3625440 : IABASE=3
529 3625440 : IBBASE=N-1
530 3625440 : JABASE=2*INC+1
531 3625440 : JBBASE=(N-2)*INC+1
532 : !
533 3625440 : DO 30 K=3,NH,2
534 126890400 : IA=IABASE
535 126890400 : IB=IBBASE
536 126890400 : JA=JABASE
537 126890400 : JB=JBBASE
538 126890400 : C=TRIGS(N+K)
539 126890400 : S=TRIGS(N+K+1)
540 624514800 : DO 20 L=1,LOT
541 0 : A(JA)=SCALE*((WORK(IA)+WORK(IB)) &
542 497624400 : +(C*(WORK(IA+1)+WORK(IB+1))+S*(WORK(IA)-WORK(IB))))
543 0 : A(JB)=SCALE*((WORK(IA)+WORK(IB)) &
544 497624400 : -(C*(WORK(IA+1)+WORK(IB+1))+S*(WORK(IA)-WORK(IB))))
545 0 : A(JA+INC)=SCALE*((C*(WORK(IA)-WORK(IB))-S*(WORK(IA+1)+WORK(IB+1))) &
546 497624400 : +(WORK(IB+1)-WORK(IA+1)))
547 0 : A(JB+INC)=SCALE*((C*(WORK(IA)-WORK(IB))-S*(WORK(IA+1)+WORK(IB+1))) &
548 497624400 : -(WORK(IB+1)-WORK(IA+1)))
549 497624400 : IA=IA+NX
550 497624400 : IB=IB+NX
551 497624400 : JA=JA+JUMP
552 497624400 : JB=JB+JUMP
553 126890400 : 20 CONTINUE
554 126890400 : IABASE=IABASE+2
555 126890400 : IBBASE=IBBASE-2
556 126890400 : JABASE=JABASE+INK
557 126890400 : JBBASE=JBBASE-INK
558 3625440 : 30 CONTINUE
559 : !
560 3625440 : IF (IABASE.NE.IBBASE) GO TO 50
561 : ! WAVENUMBER N/4 (IF IT EXISTS)
562 3625440 : IA=IABASE
563 3625440 : JA=JABASE
564 3625440 : SCALE=2.0_r8*SCALE
565 17843280 : DO 40 L=1,LOT
566 14217840 : A(JA)=SCALE*WORK(IA)
567 14217840 : A(JA+INC)=-SCALE*WORK(IA+1)
568 14217840 : IA=IA+NX
569 14217840 : JA=JA+JUMP
570 3625440 : 40 CONTINUE
571 : !
572 : 50 CONTINUE
573 3625440 : RETURN
574 : END SUBROUTINE FFT99B
575 :
576 : !================================================================================================
577 :
578 7250880 : SUBROUTINE FFT991(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN)
579 : !
580 : !-----------------------------------------------------------------------
581 : ! SUBROUTINE "FFT991" - MULTIPLE REAL/HALF-COMPLEX PERIODIC
582 : ! FAST FOURIER TRANSFORM
583 : !
584 : ! SAME AS FFT99 EXCEPT THAT ORDERING OF DATA CORRESPONDS TO
585 : ! THAT IN MRFFT2
586 : !
587 : ! PROCEDURE USED TO CONVERT TO HALF-LENGTH COMPLEX TRANSFORM
588 : ! IS GIVEN BY COOLEY, LEWIS AND WELCH (J. SOUND VIB., VOL. 12
589 : ! (1970), 315-337)
590 : !
591 : ! A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA
592 : ! WORK IS AN AREA OF SIZE (N+1)*LOT
593 : ! TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES
594 : ! IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N/2
595 : ! INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR'
596 : ! (E.G. INC=1 FOR CONSECUTIVELY STORED DATA)
597 : ! JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR
598 : ! N IS THE LENGTH OF THE DATA VECTORS
599 : ! LOT IS THE NUMBER OF DATA VECTORS
600 : ! ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT
601 : ! = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL
602 : !
603 : ! ORDERING OF COEFFICIENTS:
604 : ! A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2)
605 : ! WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED
606 : !
607 : ! ORDERING OF DATA:
608 : ! X(0),X(1),X(2),...,X(N-1)
609 : !
610 : ! VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN
611 : ! PARALLEL
612 : !
613 : ! *** N.B. N IS ASSUMED TO BE AN EVEN NUMBER
614 : !
615 : ! DEFINITION OF TRANSFORMS:
616 : ! -------------------------
617 : !
618 : ! ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N))
619 : ! WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K)
620 : !
621 : ! ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N))
622 : ! B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N))
623 : !-----------------------------------------------------------------------
624 : !
625 : use shr_kind_mod, only: r8 => shr_kind_r8
626 : implicit none
627 : !
628 : !------------------------------Arguments--------------------------------
629 : !
630 : integer IFAX(13), inc, jump, n, lot, isign
631 : real(r8) A(*), WORK(*), TRIGS(*)
632 : !
633 : !---------------------------Local variables-----------------------------
634 : !
635 : integer ibase, jbase, i, j, ia, ib, igo, ink, k
636 : integer l, la, m, nh, nfax, nx
637 : !
638 : !-----------------------------------------------------------------------
639 : !
640 7250880 : NFAX=IFAX(1)
641 7250880 : NX=N+1
642 7250880 : NH=N/2
643 7250880 : INK=INC+INC
644 7250880 : IF (ISIGN.EQ.+1) GO TO 30
645 : !
646 : ! IF NECESSARY, TRANSFER DATA TO WORK AREA
647 3625440 : IGO=50
648 3625440 : IF (MOD(NFAX,2).EQ.1) GOTO 40
649 3625440 : IBASE=1
650 3625440 : JBASE=1
651 17843280 : DO 20 L=1,LOT
652 14217840 : I=IBASE
653 14217840 : J=JBASE
654 2061586800 : DO 10 M=1,N
655 2047368960 : WORK(J)=A(I)
656 2047368960 : I=I+INC
657 2047368960 : J=J+1
658 14217840 : 10 CONTINUE
659 14217840 : IBASE=IBASE+JUMP
660 14217840 : JBASE=JBASE+NX
661 3625440 : 20 CONTINUE
662 : !
663 : IGO=60
664 3625440 : GO TO 40
665 : !
666 : ! PREPROCESSING (ISIGN=+1)
667 : ! ------------------------
668 : !
669 : 30 CONTINUE
670 3625440 : CALL FFT99A(A,WORK,TRIGS,INC,JUMP,N,LOT)
671 3625440 : IGO=60
672 : !
673 : ! COMPLEX TRANSFORM
674 : ! -----------------
675 : !
676 : 40 CONTINUE
677 7250880 : IA=1
678 7250880 : LA=1
679 36254400 : DO 80 K=1,NFAX
680 29003520 : IF (IGO.EQ.60) GO TO 60
681 : 50 CONTINUE
682 0 : CALL VPASSM(A(IA),A(IA+INC),WORK(1),WORK(2),TRIGS, &
683 14501760 : INK,2,JUMP,NX,LOT,NH,IFAX(K+1),LA)
684 14501760 : IGO=60
685 29003520 : GO TO 70
686 : 60 CONTINUE
687 0 : CALL VPASSM(WORK(1),WORK(2),A(IA),A(IA+INC),TRIGS, &
688 14501760 : 2,INK,NX,JUMP,LOT,NH,IFAX(K+1),LA)
689 14501760 : IGO=50
690 : 70 CONTINUE
691 29003520 : LA=LA*IFAX(K+1)
692 7250880 : 80 CONTINUE
693 : !
694 7250880 : IF (ISIGN.EQ.-1) GO TO 130
695 : !
696 : ! IF NECESSARY, TRANSFER DATA FROM WORK AREA
697 3625440 : IF (MOD(NFAX,2).EQ.1) GO TO 110
698 3625440 : IBASE=1
699 3625440 : JBASE=1
700 17843280 : DO 100 L=1,LOT
701 14217840 : I=IBASE
702 14217840 : J=JBASE
703 2061586800 : DO 90 M=1,N
704 2047368960 : A(J)=WORK(I)
705 2047368960 : I=I+1
706 2047368960 : J=J+INC
707 14217840 : 90 CONTINUE
708 14217840 : IBASE=IBASE+NX
709 14217840 : JBASE=JBASE+JUMP
710 3625440 : 100 CONTINUE
711 : !
712 : ! FILL IN ZEROS AT END
713 : 110 CONTINUE
714 3625440 : IB=N*INC+1
715 17843280 : DO 120 L=1,LOT
716 14217840 : A(IB)=0.0_r8
717 14217840 : A(IB+INC)=0.0_r8
718 14217840 : IB=IB+JUMP
719 3625440 : 120 CONTINUE
720 3625440 : GO TO 140
721 : !
722 : ! POSTPROCESSING (ISIGN=-1):
723 : ! --------------------------
724 : !
725 : 130 CONTINUE
726 3625440 : CALL FFT99B(WORK,A,TRIGS,INC,JUMP,N,LOT)
727 : !
728 : 140 CONTINUE
729 7250880 : RETURN
730 : END SUBROUTINE FFT991
731 :
732 : !================================================================================================
733 :
734 768 : SUBROUTINE SET99 (TRIGS, IFAX, N)
735 : !
736 : !-----------------------------------------------------------------------
737 : !
738 : use shr_kind_mod, only: r8 => shr_kind_r8
739 : use cam_logfile, only: iulog
740 : implicit none
741 : !
742 : !------------------------------Arguments--------------------------------
743 : !
744 : integer n, IFAX(13)
745 : real(r8) TRIGS(*)
746 : !
747 : !---------------------------Local variables-----------------------------
748 : !
749 : integer i
750 : !
751 : ! MODE 3 IS USED FOR REAL/HALF-COMPLEX TRANSFORMS. IT IS POSSIBLE
752 : ! TO DO COMPLEX/COMPLEX TRANSFORMS WITH OTHER VALUES OF MODE, BUT
753 : ! DOCUMENTATION OF THE DETAILS WERE NOT AVAILABLE WHEN THIS ROUTINE
754 : ! WAS WRITTEN.
755 : !
756 : integer mode
757 : DATA MODE /3/
758 : !
759 : !-----------------------------------------------------------------------
760 : !
761 768 : CALL FAX (IFAX, N, MODE)
762 768 : I = IFAX(1)
763 768 : IF (IFAX(I+1) .GT. 5 .OR. N .LE. 4) IFAX(1) = -99
764 768 : IF (IFAX(1) .LE. 0 ) THEN
765 0 : write(iulog,*) ' SET99 -- INVALID N'
766 0 : STOP 'SET99'
767 : ENDIF
768 768 : CALL FFTRIG (TRIGS, N, MODE)
769 768 : RETURN
770 : END SUBROUTINE SET99
771 :
772 : !================================================================================================
773 :
774 768 : SUBROUTINE FAX(IFAX,N,MODE)
775 : !
776 : !-----------------------------------------------------------------------
777 : !
778 : use shr_kind_mod, only: r8 => shr_kind_r8
779 : implicit none
780 : !
781 : !------------------------------Arguments--------------------------------
782 : !
783 : integer IFAX(10), n, mode
784 : !
785 : !---------------------------Local variables-----------------------------
786 : !
787 : integer ii, nfax, inc, item, i, istop, l, k, nn
788 : !
789 : !-----------------------------------------------------------------------
790 : !
791 768 : NN=N
792 768 : IF (IABS(MODE).EQ.1) GO TO 10
793 768 : IF (IABS(MODE).EQ.8) GO TO 10
794 768 : NN=N/2
795 768 : IF ((NN+NN).EQ.N) GO TO 10
796 0 : IFAX(1)=-99
797 768 : RETURN
798 : 10 K=1
799 : ! TEST FOR FACTORS OF 4
800 1536 : 20 IF (MOD(NN,4).NE.0) GO TO 30
801 768 : K=K+1
802 768 : IFAX(K)=4
803 768 : NN=NN/4
804 768 : IF (NN.EQ.1) GO TO 80
805 768 : GO TO 20
806 : ! TEST FOR EXTRA FACTOR OF 2
807 768 : 30 IF (MOD(NN,2).NE.0) GO TO 40
808 768 : K=K+1
809 768 : IFAX(K)=2
810 768 : NN=NN/2
811 768 : IF (NN.EQ.1) GO TO 80
812 : ! TEST FOR FACTORS OF 3
813 1536 : 40 IF (MOD(NN,3).NE.0) GO TO 50
814 1536 : K=K+1
815 1536 : IFAX(K)=3
816 1536 : NN=NN/3
817 1536 : IF (NN.EQ.1) GO TO 80
818 0 : GO TO 40
819 : ! NOW FIND REMAINING FACTORS
820 : 50 L=5
821 : INC=2
822 : ! INC ALTERNATELY TAKES ON VALUES 2 AND 4
823 0 : 60 IF (MOD(NN,L).NE.0) GO TO 70
824 0 : K=K+1
825 0 : IFAX(K)=L
826 0 : NN=NN/L
827 0 : IF (NN.EQ.1) GO TO 80
828 0 : GO TO 60
829 0 : 70 L=L+INC
830 0 : INC=6-INC
831 768 : GO TO 60
832 768 : 80 IFAX(1)=K-1
833 : ! IFAX(1) CONTAINS NUMBER OF FACTORS
834 768 : NFAX=IFAX(1)
835 : ! SORT FACTORS INTO ASCENDING ORDER
836 768 : IF (NFAX.EQ.1) GO TO 110
837 3072 : DO 100 II=2,NFAX
838 2304 : ISTOP=NFAX+2-II
839 6912 : DO 90 I=2,ISTOP
840 4608 : IF (IFAX(I+1).GE.IFAX(I)) GO TO 90
841 2304 : ITEM=IFAX(I)
842 2304 : IFAX(I)=IFAX(I+1)
843 4608 : IFAX(I+1)=ITEM
844 2304 : 90 CONTINUE
845 768 : 100 CONTINUE
846 : 110 CONTINUE
847 : RETURN
848 : END SUBROUTINE FAX
849 :
850 : !================================================================================================
851 :
852 768 : SUBROUTINE FFTRIG(TRIGS,N,MODE)
853 : !
854 : !-----------------------------------------------------------------------
855 : !
856 : use shr_kind_mod, only: r8 => shr_kind_r8
857 : implicit none
858 : !
859 : !------------------------------Arguments--------------------------------
860 : !
861 : integer n, mode
862 : real(r8) TRIGS(*)
863 : !
864 : !---------------------------Local variables-----------------------------
865 : !
866 : integer i, l, la, nh, imode, nn
867 : real(r8) del, pi, angle
868 : !
869 : !-----------------------------------------------------------------------
870 : !
871 768 : PI=2.0_r8*ASIN(1.0_r8)
872 768 : IMODE=IABS(MODE)
873 768 : NN=N
874 768 : IF (IMODE.GT.1.AND.IMODE.LT.6) NN=N/2
875 768 : DEL=(PI+PI)/real(NN,r8)
876 768 : L=NN+NN
877 768 : DO 10 I=1,L,2
878 55296 : ANGLE=0.5_r8*real(I-1,r8)*DEL
879 55296 : TRIGS(I)=COS(ANGLE)
880 55296 : TRIGS(I+1)=SIN(ANGLE)
881 768 : 10 CONTINUE
882 768 : IF (IMODE.EQ.1) RETURN
883 768 : IF (IMODE.EQ.8) RETURN
884 768 : DEL=0.5_r8*DEL
885 768 : NH=(NN+1)/2
886 768 : L=NH+NH
887 768 : LA=NN+NN
888 768 : DO 20 I=1,L,2
889 27648 : ANGLE=0.5_r8*real(I-1,r8)*DEL
890 27648 : TRIGS(LA+I)=COS(ANGLE)
891 27648 : TRIGS(LA+I+1)=SIN(ANGLE)
892 768 : 20 CONTINUE
893 768 : IF (IMODE.LE.3) RETURN
894 0 : DEL=0.5_r8*DEL
895 0 : LA=LA+NN
896 0 : IF (MODE.EQ.5) GO TO 40
897 0 : DO 30 I=2,NN
898 0 : ANGLE=real(I-1,r8)*DEL
899 0 : TRIGS(LA+I)=2.0_r8*SIN(ANGLE)
900 0 : 30 CONTINUE
901 0 : RETURN
902 : 40 CONTINUE
903 0 : DEL=0.5_r8*DEL
904 0 : DO 50 I=2,N
905 0 : ANGLE=real(I-1,r8)*DEL
906 0 : TRIGS(LA+I)=SIN(ANGLE)
907 0 : 50 CONTINUE
908 : RETURN
909 : END SUBROUTINE FFTRIG
910 :
911 : !================================================================================================
912 :
913 29003520 : SUBROUTINE VPASSM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,LOT,N,IFAC,LA)
914 : !
915 : !-----------------------------------------------------------------------
916 : ! SUBROUTINE "VPASSM" - MULTIPLE VERSION OF "VPASSA"
917 : ! PERFORMS ONE PASS THROUGH DATA
918 : ! AS PART OF MULTIPLE COMPLEX FFT ROUTINE
919 : ! A IS FIRST REAL INPUT VECTOR
920 : ! B IS FIRST IMAGINARY INPUT VECTOR
921 : ! C IS FIRST REAL OUTPUT VECTOR
922 : ! D IS FIRST IMAGINARY OUTPUT VECTOR
923 : ! TRIGS IS PRECALCULATED TABLE OF SINES " COSINES
924 : ! INC1 IS ADDRESSING INCREMENT FOR A AND B
925 : ! INC2 IS ADDRESSING INCREMENT FOR C AND D
926 : ! INC3 IS ADDRESSING INCREMENT BETWEEN A"S & B"S
927 : ! INC4 IS ADDRESSING INCREMENT BETWEEN C"S & D"S
928 : ! LOT IS THE NUMBER OF VECTORS
929 : ! N IS LENGTH OF VECTORS
930 : ! IFAC IS CURRENT FACTOR OF N
931 : ! LA IS PRODUCT OF PREVIOUS FACTORS
932 : !-----------------------------------------------------------------------
933 : !
934 : use shr_kind_mod, only: r8 => shr_kind_r8
935 : implicit none
936 : !
937 : !------------------------------Arguments--------------------------------
938 : !
939 : integer inc1, inc2, inc3, inc4, lot, n, ifac, la
940 : real(r8) A(*),B(*),C(*),D(*),TRIGS(*)
941 : !
942 : !---------------------------Local variables-----------------------------
943 : !
944 : integer ie, je, ke, kd, ib, ja, ia, i, l, jb, igo, jink
945 : integer iink, m, jbase, ibase, jump, j, kc, jc, jd, id
946 : integer ic, k, la1, ijk, kb
947 :
948 : real(r8) s3, c3, s4, c4, c2, s2, s1, c1
949 :
950 : real(r8) sin36, cos36, sin72, cos72, sin60
951 : DATA SIN36/0.587785252292473_r8/,COS36/0.809016994374947_r8/, &
952 : SIN72/0.951056516295154_r8/,COS72/0.309016994374947_r8/, &
953 : SIN60/0.866025403784437_r8/
954 : !
955 : !-----------------------------------------------------------------------
956 : !
957 29003520 : M=N/IFAC
958 29003520 : IINK=M*INC1
959 29003520 : JINK=LA*INC2
960 29003520 : JUMP=(IFAC-1)*JINK
961 29003520 : IBASE=0
962 29003520 : JBASE=0
963 29003520 : IGO=IFAC-1
964 29003520 : IF (IGO.GT.4) RETURN
965 29003520 : GO TO (10,50,90,130),IGO
966 : !
967 : ! CODING FOR FACTOR 2
968 : !
969 7250880 : 10 IA=1
970 7250880 : JA=1
971 7250880 : IB=IA+IINK
972 7250880 : JB=JA+JINK
973 14501760 : DO 20 L=1,LA
974 7250880 : I=IBASE
975 7250880 : J=JBASE
976 35686560 : DO 15 IJK=1,LOT
977 28435680 : C(JA+J)=A(IA+I)+A(IB+I)
978 28435680 : D(JA+J)=B(IA+I)+B(IB+I)
979 28435680 : C(JB+J)=A(IA+I)-A(IB+I)
980 28435680 : D(JB+J)=B(IA+I)-B(IB+I)
981 28435680 : I=I+INC3
982 28435680 : J=J+INC4
983 7250880 : 15 CONTINUE
984 7250880 : IBASE=IBASE+INC1
985 7250880 : JBASE=JBASE+INC2
986 7250880 : 20 CONTINUE
987 7250880 : IF (LA.EQ.M) RETURN
988 7250880 : LA1=LA+1
989 7250880 : JBASE=JBASE+JUMP
990 7250880 : DO 40 K=LA1,M,LA
991 253780800 : KB=K+K-2
992 253780800 : C1=TRIGS(KB+1)
993 253780800 : S1=TRIGS(KB+2)
994 507561600 : DO 30 L=1,LA
995 253780800 : I=IBASE
996 253780800 : J=JBASE
997 1249029600 : DO 25 IJK=1,LOT
998 995248800 : C(JA+J)=A(IA+I)+A(IB+I)
999 995248800 : D(JA+J)=B(IA+I)+B(IB+I)
1000 995248800 : C(JB+J)=C1*(A(IA+I)-A(IB+I))-S1*(B(IA+I)-B(IB+I))
1001 995248800 : D(JB+J)=S1*(A(IA+I)-A(IB+I))+C1*(B(IA+I)-B(IB+I))
1002 995248800 : I=I+INC3
1003 995248800 : J=J+INC4
1004 253780800 : 25 CONTINUE
1005 253780800 : IBASE=IBASE+INC1
1006 253780800 : JBASE=JBASE+INC2
1007 253780800 : 30 CONTINUE
1008 253780800 : JBASE=JBASE+JUMP
1009 7250880 : 40 CONTINUE
1010 : RETURN
1011 : !
1012 : ! CODING FOR FACTOR 3
1013 : !
1014 14501760 : 50 IA=1
1015 14501760 : JA=1
1016 14501760 : IB=IA+IINK
1017 14501760 : JB=JA+JINK
1018 14501760 : IC=IB+IINK
1019 14501760 : JC=JB+JINK
1020 72508800 : DO 60 L=1,LA
1021 58007040 : I=IBASE
1022 58007040 : J=JBASE
1023 285492480 : DO 55 IJK=1,LOT
1024 227485440 : C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))
1025 227485440 : D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I))
1026 227485440 : C(JB+J)=(A(IA+I)-0.5_r8*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I)))
1027 227485440 : C(JC+J)=(A(IA+I)-0.5_r8*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I)))
1028 227485440 : D(JB+J)=(B(IA+I)-0.5_r8*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I)))
1029 227485440 : D(JC+J)=(B(IA+I)-0.5_r8*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I)))
1030 227485440 : I=I+INC3
1031 227485440 : J=J+INC4
1032 58007040 : 55 CONTINUE
1033 58007040 : IBASE=IBASE+INC1
1034 58007040 : JBASE=JBASE+INC2
1035 14501760 : 60 CONTINUE
1036 14501760 : IF (LA.EQ.M) RETURN
1037 14501760 : LA1=LA+1
1038 14501760 : JBASE=JBASE+JUMP
1039 14501760 : DO 80 K=LA1,M,LA
1040 101512320 : KB=K+K-2
1041 101512320 : KC=KB+KB
1042 101512320 : C1=TRIGS(KB+1)
1043 101512320 : S1=TRIGS(KB+2)
1044 101512320 : C2=TRIGS(KC+1)
1045 101512320 : S2=TRIGS(KC+2)
1046 391547520 : DO 70 L=1,LA
1047 290035200 : I=IBASE
1048 290035200 : J=JBASE
1049 1427462400 : DO 65 IJK=1,LOT
1050 1137427200 : C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))
1051 1137427200 : D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I))
1052 0 : C(JB+J)= &
1053 : C1*((A(IA+I)-0.5_r8*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I)))) &
1054 1137427200 : -S1*((B(IA+I)-0.5_r8*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))))
1055 : D(JB+J)= &
1056 : S1*((A(IA+I)-0.5_r8*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I)))) &
1057 1137427200 : +C1*((B(IA+I)-0.5_r8*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))))
1058 0 : C(JC+J)= &
1059 : C2*((A(IA+I)-0.5_r8*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I)))) &
1060 1137427200 : -S2*((B(IA+I)-0.5_r8*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))))
1061 : D(JC+J)= &
1062 : S2*((A(IA+I)-0.5_r8*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I)))) &
1063 1137427200 : +C2*((B(IA+I)-0.5_r8*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))))
1064 1137427200 : I=I+INC3
1065 1137427200 : J=J+INC4
1066 290035200 : 65 CONTINUE
1067 290035200 : IBASE=IBASE+INC1
1068 290035200 : JBASE=JBASE+INC2
1069 101512320 : 70 CONTINUE
1070 101512320 : JBASE=JBASE+JUMP
1071 14501760 : 80 CONTINUE
1072 : RETURN
1073 : !
1074 : ! CODING FOR FACTOR 4
1075 : !
1076 7250880 : 90 IA=1
1077 7250880 : JA=1
1078 7250880 : IB=IA+IINK
1079 7250880 : JB=JA+JINK
1080 7250880 : IC=IB+IINK
1081 7250880 : JC=JB+JINK
1082 7250880 : ID=IC+IINK
1083 7250880 : JD=JC+JINK
1084 137766720 : DO 100 L=1,LA
1085 130515840 : I=IBASE
1086 130515840 : J=JBASE
1087 642358080 : DO 95 IJK=1,LOT
1088 511842240 : C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))
1089 511842240 : C(JC+J)=(A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))
1090 511842240 : D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I))
1091 511842240 : D(JC+J)=(B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I))
1092 511842240 : C(JB+J)=(A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I))
1093 511842240 : C(JD+J)=(A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I))
1094 511842240 : D(JB+J)=(B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I))
1095 511842240 : D(JD+J)=(B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I))
1096 511842240 : I=I+INC3
1097 511842240 : J=J+INC4
1098 130515840 : 95 CONTINUE
1099 130515840 : IBASE=IBASE+INC1
1100 130515840 : JBASE=JBASE+INC2
1101 7250880 : 100 CONTINUE
1102 7250880 : IF (LA.EQ.M) RETURN
1103 0 : LA1=LA+1
1104 0 : JBASE=JBASE+JUMP
1105 0 : DO 120 K=LA1,M,LA
1106 0 : KB=K+K-2
1107 0 : KC=KB+KB
1108 0 : KD=KC+KB
1109 0 : C1=TRIGS(KB+1)
1110 0 : S1=TRIGS(KB+2)
1111 0 : C2=TRIGS(KC+1)
1112 0 : S2=TRIGS(KC+2)
1113 0 : C3=TRIGS(KD+1)
1114 0 : S3=TRIGS(KD+2)
1115 0 : DO 110 L=1,LA
1116 0 : I=IBASE
1117 0 : J=JBASE
1118 0 : DO 105 IJK=1,LOT
1119 0 : C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))
1120 0 : D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I))
1121 0 : C(JC+J)= &
1122 : C2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))) &
1123 0 : -S2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I)))
1124 : D(JC+J)= &
1125 : S2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))) &
1126 0 : +C2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I)))
1127 0 : C(JB+J)= &
1128 : C1*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I))) &
1129 0 : -S1*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I)))
1130 : D(JB+J)= &
1131 : S1*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I))) &
1132 0 : +C1*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I)))
1133 0 : C(JD+J)= &
1134 : C3*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I))) &
1135 0 : -S3*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I)))
1136 : D(JD+J)= &
1137 : S3*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I))) &
1138 0 : +C3*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I)))
1139 0 : I=I+INC3
1140 0 : J=J+INC4
1141 0 : 105 CONTINUE
1142 0 : IBASE=IBASE+INC1
1143 0 : JBASE=JBASE+INC2
1144 0 : 110 CONTINUE
1145 0 : JBASE=JBASE+JUMP
1146 0 : 120 CONTINUE
1147 : RETURN
1148 : !
1149 : ! CODING FOR FACTOR 5
1150 : !
1151 0 : 130 IA=1
1152 0 : JA=1
1153 0 : IB=IA+IINK
1154 0 : JB=JA+JINK
1155 0 : IC=IB+IINK
1156 0 : JC=JB+JINK
1157 0 : ID=IC+IINK
1158 0 : JD=JC+JINK
1159 0 : IE=ID+IINK
1160 0 : JE=JD+JINK
1161 0 : DO 140 L=1,LA
1162 0 : I=IBASE
1163 0 : J=JBASE
1164 0 : DO 135 IJK=1,LOT
1165 0 : C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))
1166 0 : D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I))
1167 0 : C(JB+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) &
1168 0 : -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))
1169 0 : C(JE+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) &
1170 0 : +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))
1171 : D(JB+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) &
1172 0 : +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))
1173 : D(JE+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) &
1174 0 : -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))
1175 0 : C(JC+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) &
1176 0 : -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))
1177 0 : C(JD+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) &
1178 0 : +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))
1179 : D(JC+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) &
1180 0 : +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))
1181 : D(JD+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) &
1182 0 : -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))
1183 0 : I=I+INC3
1184 0 : J=J+INC4
1185 0 : 135 CONTINUE
1186 0 : IBASE=IBASE+INC1
1187 0 : JBASE=JBASE+INC2
1188 0 : 140 CONTINUE
1189 0 : IF (LA.EQ.M) RETURN
1190 0 : LA1=LA+1
1191 0 : JBASE=JBASE+JUMP
1192 0 : DO 160 K=LA1,M,LA
1193 0 : KB=K+K-2
1194 0 : KC=KB+KB
1195 0 : KD=KC+KB
1196 0 : KE=KD+KB
1197 0 : C1=TRIGS(KB+1)
1198 0 : S1=TRIGS(KB+2)
1199 0 : C2=TRIGS(KC+1)
1200 0 : S2=TRIGS(KC+2)
1201 0 : C3=TRIGS(KD+1)
1202 0 : S3=TRIGS(KD+2)
1203 0 : C4=TRIGS(KE+1)
1204 0 : S4=TRIGS(KE+2)
1205 0 : DO 150 L=1,LA
1206 0 : I=IBASE
1207 0 : J=JBASE
1208 0 : DO 145 IJK=1,LOT
1209 0 : C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))
1210 0 : D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I))
1211 0 : C(JB+J)= &
1212 : C1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) &
1213 : -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))) &
1214 : -S1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) &
1215 0 : +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))
1216 : D(JB+J)= &
1217 : S1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) &
1218 : -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))) &
1219 : +C1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) &
1220 0 : +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))
1221 0 : C(JE+J)= &
1222 : C4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) &
1223 : +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))) &
1224 : -S4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) &
1225 0 : -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))
1226 : D(JE+J)= &
1227 : S4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) &
1228 : +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))) &
1229 : +C4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) &
1230 0 : -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))))
1231 0 : C(JC+J)= &
1232 : C2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) &
1233 : -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))) &
1234 : -S2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) &
1235 0 : +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))
1236 : D(JC+J)= &
1237 : S2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) &
1238 : -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))) &
1239 : +C2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) &
1240 0 : +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))
1241 0 : C(JD+J)= &
1242 : C3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) &
1243 : +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))) &
1244 : -S3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) &
1245 0 : -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))
1246 : D(JD+J)= &
1247 : S3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) &
1248 : +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))) &
1249 : +C3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) &
1250 0 : -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))))
1251 0 : I=I+INC3
1252 0 : J=J+INC4
1253 0 : 145 CONTINUE
1254 0 : IBASE=IBASE+INC1
1255 0 : JBASE=JBASE+INC2
1256 0 : 150 CONTINUE
1257 0 : JBASE=JBASE+JUMP
1258 0 : 160 CONTINUE
1259 : RETURN
1260 : END SUBROUTINE VPASSM
1261 :
1262 : !===============================================================================
1263 :
1264 :
|
