1* 2* $Id$ 3* 4CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 5C 6C File: sint.f 7C 8C Library: FFTPACK 4.1 9C 10C Author: Paul N. Swarztrauber 11C National Center for Atmospheric Research 12C PO 3000, Boulder, Colorado 13C 14C Date: Wed Mar 29 18:31:13 MST 1995 15C 16C Description: Forward and backward, 1D real sin FFT 17C 18 SUBROUTINE SINT(N,X,WSAVE) 19 DOUBLE PRECISION X(*),WSAVE(*) 20 21 NP1 = N + 1 22 IW1 = N/2 + 1 23 IW2 = IW1 + NP1 24 IW3 = IW2 + NP1 25 CALL SINT1(N,X,WSAVE,WSAVE(IW1),WSAVE(IW2),WSAVE(IW3)) 26 RETURN 27 END 28CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 29C 30C File: sinti.f 31C 32C Library: FFTPACK 4.1 33C 34C Author: Paul N. Swarztrauber 35C National Center for Atmospheric Research 36C PO 3000, Boulder, Colorado 37C 38C Date: Wed Mar 29 18:31:13 MST 1995 39C 40C Description: Initialization routine for SINT 41C 42 SUBROUTINE SINTI(N,WSAVE) 43 DOUBLE PRECISION PI 44 DOUBLE PRECISION PIMACH 45 DOUBLE PRECISION DUM 46 DOUBLE PRECISION DT 47 DOUBLE PRECISION WSAVE(*) 48 49 PI = PIMACH(DUM) 50 IF (N.LE.1) RETURN 51 NS2 = N/2 52 NP1 = N + 1 53 DT = PI/DBLE(NP1) 54 DO 101 K = 1,NS2 55 WSAVE(K) = 2.D0*SIN(K*DT) 56 101 CONTINUE 57 CALL RFFTI(NP1,WSAVE(NS2+1)) 58 RETURN 59 END 60CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 61C 62C File: sint1.f 63C 64C Library: FFTPACK 4.1 65C 66C Author: Paul N. Swarztrauber 67C National Center for Atmospheric Research 68C PO 3000, Boulder, Colorado 69C 70C Date: Wed Mar 29 18:31:13 MST 1995 71C 72C Description: Lower-level auxiliary routine 73C 74 SUBROUTINE SINT1(N,WAR,WAS,XH,X,IFAC) 75 DOUBLE PRECISION SQRT3 76 DOUBLE PRECISION XHOLD 77 DOUBLE PRECISION T1 78 DOUBLE PRECISION T2 79 DOUBLE PRECISION WAR(*),WAS(*),X(*),XH(*) 80 integer*4 n4 81 INTEGER IFAC(*) 82 DATA SQRT3/1.73205080756888D0/ 83C 84C FFTPACK 5.0 auxiliary routine 85C 86 DO 100 I = 1,N 87 XH(I) = WAR(I) 88 WAR(I) = X(I) 89 100 CONTINUE 90#if 1 91 n4=n-2 92 IF (n4) 101,102,103 93#else 94 IF (N-2) 101,102,103 95#endif 96 101 XH(1) = XH(1) + XH(1) 97 GO TO 106 98 102 XHOLD = SQRT3* (XH(1)+XH(2)) 99 XH(2) = SQRT3* (XH(1)-XH(2)) 100 XH(1) = XHOLD 101 GO TO 106 102 103 NP1 = N + 1 103 NS2 = N/2 104 X(1) = 0.D0 105 DO 104 K = 1,NS2 106 KC = NP1 - K 107 T1 = XH(K) - XH(KC) 108 T2 = WAS(K)* (XH(K)+XH(KC)) 109 X(K+1) = T1 + T2 110 X(KC+1) = T2 - T1 111 104 CONTINUE 112 MODN = MOD(N,2) 113 IF (MODN.NE.0) X(NS2+2) = 4.D0*XH(NS2+1) 114 CALL RFFTF1(NP1,X,XH,WAR,IFAC) 115 XH(1) = .5D0*X(1) 116 DO 105 I = 3,N,2 117 XH(I-1) = -X(I) 118 XH(I) = XH(I-2) + X(I-1) 119 105 CONTINUE 120 IF (MODN.NE.0) GO TO 106 121 XH(N) = -X(N+1) 122 106 DO 107 I = 1,N 123 X(I) = WAR(I) 124 WAR(I) = XH(I) 125 107 CONTINUE 126 RETURN 127 END 128CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 129C 130C File: rffti.f 131C 132C Library: FFTPACK 4.1 133C 134C Author: Paul N. Swarztrauber 135C National Center for Atmospheric Research 136C PO 3000, Boulder, Colorado 137C 138C Date: Wed Mar 29 18:31:13 MST 1995 139C 140C Description: Initialization routine for RFFTB, RFFTF 141C 142 SUBROUTINE RFFTI(N,WSAVE) 143 DOUBLE PRECISION WSAVE(*) 144 145 IF (N.EQ.1) RETURN 146 CALL RFFTI1(N,WSAVE(N+1),WSAVE(2*N+1)) 147 RETURN 148 END 149CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 150C 151C File: rfftf1.f 152C 153C Library: FFTPACK 4.1 154C 155C Author: Paul N. Swarztrauber 156C National Center for Atmospheric Research 157C PO 3000, Boulder, Colorado 158C 159C Date: Wed Mar 29 18:31:13 MST 1995 160C 161C Description: Lower-level auxiliary routine 162C 163 SUBROUTINE RFFTF1(N,C,CH,WA,IFAC) 164 DOUBLE PRECISION CH(*),C(*),WA(*) 165 INTEGER IFAC(*) 166C 167C FFTPACK 5.0 auxiliary routine 168C 169 NF = IFAC(2) 170 NA = 1 171 L2 = N 172 IW = N 173 DO 111 K1 = 1,NF 174 KH = NF - K1 175 IP = IFAC(KH+3) 176 L1 = L2/IP 177 IDO = N/L2 178 IDL1 = IDO*L1 179 IW = IW - (IP-1)*IDO 180 NA = 1 - NA 181 IF (IP.NE.4) GO TO 102 182 IX2 = IW + IDO 183 IX3 = IX2 + IDO 184 IF (NA.NE.0) GO TO 101 185 CALL RADF4(IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3)) 186 GO TO 110 187 101 CALL RADF4(IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3)) 188 GO TO 110 189 102 IF (IP.NE.2) GO TO 104 190 IF (NA.NE.0) GO TO 103 191 CALL RADF2(IDO,L1,C,CH,WA(IW)) 192 GO TO 110 193 103 CALL RADF2(IDO,L1,CH,C,WA(IW)) 194 GO TO 110 195 104 IF (IP.NE.3) GO TO 106 196 IX2 = IW + IDO 197 IF (NA.NE.0) GO TO 105 198 CALL RADF3(IDO,L1,C,CH,WA(IW),WA(IX2)) 199 GO TO 110 200 105 CALL RADF3(IDO,L1,CH,C,WA(IW),WA(IX2)) 201 GO TO 110 202 106 IF (IP.NE.5) GO TO 108 203 IX2 = IW + IDO 204 IX3 = IX2 + IDO 205 IX4 = IX3 + IDO 206 IF (NA.NE.0) GO TO 107 207 CALL RADF5(IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4)) 208 GO TO 110 209 107 CALL RADF5(IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4)) 210 GO TO 110 211 108 IF (IDO.EQ.1) NA = 1 - NA 212 IF (NA.NE.0) GO TO 109 213 CALL RADFG(IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW)) 214 NA = 1 215 GO TO 110 216 109 CALL RADFG(IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW)) 217 NA = 0 218 110 L2 = L1 219 111 CONTINUE 220 IF (NA.EQ.1) RETURN 221 DO 112 I = 1,N 222 C(I) = CH(I) 223 112 CONTINUE 224 RETURN 225 END 226CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 227C 228C File: rffti1.f 229C 230C Library: FFTPACK 4.1 231C 232C Author: Paul N. Swarztrauber 233C National Center for Atmospheric Research 234C PO 3000, Boulder, Colorado 235C 236C Date: Wed Mar 29 18:31:13 MST 1995 237C 238C Description: Lower-level auxiliary routine 239C 240 SUBROUTINE RFFTI1(N,WA,IFAC) 241 DOUBLE PRECISION TPI 242 DOUBLE PRECISION PIMACH 243 DOUBLE PRECISION DUM 244 DOUBLE PRECISION ARGH 245 DOUBLE PRECISION ARGLD 246 DOUBLE PRECISION FI 247 DOUBLE PRECISION ARG 248 DOUBLE PRECISION WA(*) 249 INTEGER IFAC(*),NTRYH(4) 250 integer*4 n4 251 DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/ 252C 253C FFTPACK 5.0 auxiliary routine 254C 255 NL = N 256 NF = 0 257 NTRY = 0 258 J = 0 259 101 J = J + 1 260#if 1 261 n4=j-4 262 IF (n4) 102,102,103 263#else 264 IF (J-4) 102,102,103 265#endif 266 102 NTRY = NTRYH(J) 267 GO TO 104 268 103 NTRY = NTRY + 2 269 104 NQ = NL/NTRY 270 NR = NL - NTRY*NQ 271#if 1 272 n4=NR 273 IF (n4) 101,105,101 274#else 275 IF (NR) 101,105,101 276#endif 277 105 NF = NF + 1 278 IFAC(NF+2) = NTRY 279 NL = NQ 280 IF (NTRY.NE.2) GO TO 107 281 IF (NF.EQ.1) GO TO 107 282 DO 106 I = 2,NF 283 IB = NF - I + 2 284 IFAC(IB+2) = IFAC(IB+1) 285 106 CONTINUE 286 IFAC(3) = 2 287 107 IF (NL.NE.1) GO TO 104 288 IFAC(1) = N 289 IFAC(2) = NF 290 TPI = 2.0D0*PIMACH(DUM) 291 ARGH = TPI/DBLE(N) 292 IS = 0 293 NFM1 = NF - 1 294 L1 = 1 295 IF (NFM1.EQ.0) RETURN 296 DO 110 K1 = 1,NFM1 297 IP = IFAC(K1+2) 298 LD = 0 299 L2 = L1*IP 300 IDO = N/L2 301 IPM = IP - 1 302 DO 109 J = 1,IPM 303 LD = LD + L1 304 I = IS 305 ARGLD = DBLE(LD)*ARGH 306 FI = 0.D0 307 DO 108 II = 3,IDO,2 308 I = I + 2 309 FI = FI + 1.D0 310 ARG = FI*ARGLD 311 WA(I-1) = COS(ARG) 312 WA(I) = SIN(ARG) 313 108 CONTINUE 314 IS = IS + IDO 315 109 CONTINUE 316 L1 = L2 317 110 CONTINUE 318 RETURN 319 END 320CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 321C 322C File: radf2.f 323C 324C Library: FFTPACK 4.1 325C 326C Author: Paul N. Swarztrauber 327C National Center for Atmospheric Research 328C PO 3000, Boulder, Colorado 329C 330C Date: Wed Mar 29 18:31:13 MST 1995 331C 332C Description: Lower-level auxiliary routine 333C 334 SUBROUTINE RADF2(IDO,L1,CC,CH,WA1) 335 DOUBLE PRECISION TR2 336 DOUBLE PRECISION TI2 337 DOUBLE PRECISION CH(IDO,2,L1),CC(IDO,L1,2),WA1(*) 338 integer*4 n4 339C 340C FFTPACK 5.0 auxiliary routine 341C 342 DO 101 K = 1,L1 343 CH(1,1,K) = CC(1,K,1) + CC(1,K,2) 344 CH(IDO,2,K) = CC(1,K,1) - CC(1,K,2) 345 101 CONTINUE 346#if 1 347 n4=ido-2 348 IF (n4) 107,105,102 349#else 350 IF (IDO-2) 107,105,102 351#endif 352 102 IDP2 = IDO + 2 353 DO 104 K = 1,L1 354 DO 103 I = 3,IDO,2 355 IC = IDP2 - I 356 TR2 = WA1(I-2)*CC(I-1,K,2) + WA1(I-1)*CC(I,K,2) 357 TI2 = WA1(I-2)*CC(I,K,2) - WA1(I-1)*CC(I-1,K,2) 358 CH(I,1,K) = CC(I,K,1) + TI2 359 CH(IC,2,K) = TI2 - CC(I,K,1) 360 CH(I-1,1,K) = CC(I-1,K,1) + TR2 361 CH(IC-1,2,K) = CC(I-1,K,1) - TR2 362 103 CONTINUE 363 104 CONTINUE 364 IF (MOD(IDO,2).EQ.1) RETURN 365 105 DO 106 K = 1,L1 366 CH(1,2,K) = -CC(IDO,K,2) 367 CH(IDO,1,K) = CC(IDO,K,1) 368 106 CONTINUE 369 107 RETURN 370 END 371CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 372C 373C File: radf3.f 374C 375C Library: FFTPACK 4.1 376C 377C Author: Paul N. Swarztrauber 378C National Center for Atmospheric Research 379C PO 3000, Boulder, Colorado 380C 381C Date: Wed Mar 29 18:31:13 MST 1995 382C 383C Description: Lower-level auxiliary routine 384C 385 SUBROUTINE RADF3(IDO,L1,CC,CH,WA1,WA2) 386 DOUBLE PRECISION TAUR 387 DOUBLE PRECISION TAUI 388 DOUBLE PRECISION CR2 389 DOUBLE PRECISION DR2 390 DOUBLE PRECISION DI2 391 DOUBLE PRECISION DR3 392 DOUBLE PRECISION DI3 393 DOUBLE PRECISION CI2 394 DOUBLE PRECISION TR2 395 DOUBLE PRECISION TI2 396 DOUBLE PRECISION TR3 397 DOUBLE PRECISION TI3 398 DOUBLE PRECISION CH(IDO,3,L1),CC(IDO,L1,3),WA1(*),WA2(*) 399 DATA TAUR,TAUI/-.5D0,.866025403784439D0/ 400C 401C FFTPACK 5.0 auxiliary routine 402C 403 DO 101 K = 1,L1 404 CR2 = CC(1,K,2) + CC(1,K,3) 405 CH(1,1,K) = CC(1,K,1) + CR2 406 CH(1,3,K) = TAUI* (CC(1,K,3)-CC(1,K,2)) 407 CH(IDO,2,K) = CC(1,K,1) + TAUR*CR2 408 101 CONTINUE 409 IF (IDO.EQ.1) RETURN 410 IDP2 = IDO + 2 411 DO 103 K = 1,L1 412 DO 102 I = 3,IDO,2 413 IC = IDP2 - I 414 DR2 = WA1(I-2)*CC(I-1,K,2) + WA1(I-1)*CC(I,K,2) 415 DI2 = WA1(I-2)*CC(I,K,2) - WA1(I-1)*CC(I-1,K,2) 416 DR3 = WA2(I-2)*CC(I-1,K,3) + WA2(I-1)*CC(I,K,3) 417 DI3 = WA2(I-2)*CC(I,K,3) - WA2(I-1)*CC(I-1,K,3) 418 CR2 = DR2 + DR3 419 CI2 = DI2 + DI3 420 CH(I-1,1,K) = CC(I-1,K,1) + CR2 421 CH(I,1,K) = CC(I,K,1) + CI2 422 TR2 = CC(I-1,K,1) + TAUR*CR2 423 TI2 = CC(I,K,1) + TAUR*CI2 424 TR3 = TAUI* (DI2-DI3) 425 TI3 = TAUI* (DR3-DR2) 426 CH(I-1,3,K) = TR2 + TR3 427 CH(IC-1,2,K) = TR2 - TR3 428 CH(I,3,K) = TI2 + TI3 429 CH(IC,2,K) = TI3 - TI2 430 102 CONTINUE 431 103 CONTINUE 432 RETURN 433 END 434CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 435C 436C File: radf4.f 437C 438C Library: FFTPACK 4.1 439C 440C Author: Paul N. Swarztrauber 441C National Center for Atmospheric Research 442C PO 3000, Boulder, Colorado 443C 444C Date: Wed Mar 29 18:31:13 MST 1995 445C 446C Description: Lower-level auxiliary routine 447C 448 SUBROUTINE RADF4(IDO,L1,CC,CH,WA1,WA2,WA3) 449 DOUBLE PRECISION HSQT2 450 DOUBLE PRECISION TR1 451 DOUBLE PRECISION TR2 452 DOUBLE PRECISION CR2 453 DOUBLE PRECISION CI2 454 DOUBLE PRECISION CR3 455 DOUBLE PRECISION CI3 456 DOUBLE PRECISION CR4 457 DOUBLE PRECISION CI4 458 DOUBLE PRECISION TR4 459 DOUBLE PRECISION TI1 460 DOUBLE PRECISION TI4 461 DOUBLE PRECISION TI2 462 DOUBLE PRECISION TI3 463 DOUBLE PRECISION TR3 464 DOUBLE PRECISION CC(IDO,L1,4),CH(IDO,4,L1),WA1(*),WA2(*),WA3(*) 465 integer*4 n4 466 DATA HSQT2/.7071067811865475D0/ 467C 468C FFTPACK 5.0 auxiliary routine 469C 470 DO 101 K = 1,L1 471 TR1 = CC(1,K,2) + CC(1,K,4) 472 TR2 = CC(1,K,1) + CC(1,K,3) 473 CH(1,1,K) = TR1 + TR2 474 CH(IDO,4,K) = TR2 - TR1 475 CH(IDO,2,K) = CC(1,K,1) - CC(1,K,3) 476 CH(1,3,K) = CC(1,K,4) - CC(1,K,2) 477 101 CONTINUE 478#if 1 479 n4=ido-2 480 IF (n4) 107,105,102 481#else 482 IF (IDO-2) 107,105,102 483#endif 484 102 IDP2 = IDO + 2 485 DO 104 K = 1,L1 486 DO 103 I = 3,IDO,2 487 IC = IDP2 - I 488 CR2 = WA1(I-2)*CC(I-1,K,2) + WA1(I-1)*CC(I,K,2) 489 CI2 = WA1(I-2)*CC(I,K,2) - WA1(I-1)*CC(I-1,K,2) 490 CR3 = WA2(I-2)*CC(I-1,K,3) + WA2(I-1)*CC(I,K,3) 491 CI3 = WA2(I-2)*CC(I,K,3) - WA2(I-1)*CC(I-1,K,3) 492 CR4 = WA3(I-2)*CC(I-1,K,4) + WA3(I-1)*CC(I,K,4) 493 CI4 = WA3(I-2)*CC(I,K,4) - WA3(I-1)*CC(I-1,K,4) 494 TR1 = CR2 + CR4 495 TR4 = CR4 - CR2 496 TI1 = CI2 + CI4 497 TI4 = CI2 - CI4 498 TI2 = CC(I,K,1) + CI3 499 TI3 = CC(I,K,1) - CI3 500 TR2 = CC(I-1,K,1) + CR3 501 TR3 = CC(I-1,K,1) - CR3 502 CH(I-1,1,K) = TR1 + TR2 503 CH(IC-1,4,K) = TR2 - TR1 504 CH(I,1,K) = TI1 + TI2 505 CH(IC,4,K) = TI1 - TI2 506 CH(I-1,3,K) = TI4 + TR3 507 CH(IC-1,2,K) = TR3 - TI4 508 CH(I,3,K) = TR4 + TI3 509 CH(IC,2,K) = TR4 - TI3 510 103 CONTINUE 511 104 CONTINUE 512 IF (MOD(IDO,2).EQ.1) RETURN 513 105 CONTINUE 514 DO 106 K = 1,L1 515 TI1 = -HSQT2* (CC(IDO,K,2)+CC(IDO,K,4)) 516 TR1 = HSQT2* (CC(IDO,K,2)-CC(IDO,K,4)) 517 CH(IDO,1,K) = TR1 + CC(IDO,K,1) 518 CH(IDO,3,K) = CC(IDO,K,1) - TR1 519 CH(1,2,K) = TI1 - CC(IDO,K,3) 520 CH(1,4,K) = TI1 + CC(IDO,K,3) 521 106 CONTINUE 522 107 RETURN 523 END 524CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 525C 526C File: radf5.f 527C 528C Library: FFTPACK 4.1 529C 530C Author: Paul N. Swarztrauber 531C National Center for Atmospheric Research 532C PO 3000, Boulder, Colorado 533C 534C Date: Wed Mar 29 18:31:13 MST 1995 535C 536C Description: Lower-level auxiliary routine 537C 538 SUBROUTINE RADF5(IDO,L1,CC,CH,WA1,WA2,WA3,WA4) 539 DOUBLE PRECISION TR11 540 DOUBLE PRECISION TI11 541 DOUBLE PRECISION TR12 542 DOUBLE PRECISION TI12 543 DOUBLE PRECISION CR2 544 DOUBLE PRECISION CI5 545 DOUBLE PRECISION CR3 546 DOUBLE PRECISION CI4 547 DOUBLE PRECISION DR2 548 DOUBLE PRECISION DI2 549 DOUBLE PRECISION DR3 550 DOUBLE PRECISION DI3 551 DOUBLE PRECISION DR4 552 DOUBLE PRECISION DI4 553 DOUBLE PRECISION DR5 554 DOUBLE PRECISION DI5 555 DOUBLE PRECISION CR5 556 DOUBLE PRECISION CI2 557 DOUBLE PRECISION CR4 558 DOUBLE PRECISION CI3 559 DOUBLE PRECISION TR2 560 DOUBLE PRECISION TI2 561 DOUBLE PRECISION TR3 562 DOUBLE PRECISION TI3 563 DOUBLE PRECISION TR5 564 DOUBLE PRECISION TI5 565 DOUBLE PRECISION TR4 566 DOUBLE PRECISION TI4 567 DOUBLE PRECISION CC(IDO,L1,5),CH(IDO,5,L1),WA1(*),WA2(*),WA3(*), 568 + WA4(*) 569 DATA TR11,TI11,TR12,TI12/.309016994374947D0,.951056516295154D0, 570 + -.809016994374947D0,.587785252292473D0/ 571C 572C FFTPACK 5.0 auxiliary routine 573C 574 DO 101 K = 1,L1 575 CR2 = CC(1,K,5) + CC(1,K,2) 576 CI5 = CC(1,K,5) - CC(1,K,2) 577 CR3 = CC(1,K,4) + CC(1,K,3) 578 CI4 = CC(1,K,4) - CC(1,K,3) 579 CH(1,1,K) = CC(1,K,1) + CR2 + CR3 580 CH(IDO,2,K) = CC(1,K,1) + TR11*CR2 + TR12*CR3 581 CH(1,3,K) = TI11*CI5 + TI12*CI4 582 CH(IDO,4,K) = CC(1,K,1) + TR12*CR2 + TR11*CR3 583 CH(1,5,K) = TI12*CI5 - TI11*CI4 584 101 CONTINUE 585 IF (IDO.EQ.1) RETURN 586 IDP2 = IDO + 2 587 DO 103 K = 1,L1 588 DO 102 I = 3,IDO,2 589 IC = IDP2 - I 590 DR2 = WA1(I-2)*CC(I-1,K,2) + WA1(I-1)*CC(I,K,2) 591 DI2 = WA1(I-2)*CC(I,K,2) - WA1(I-1)*CC(I-1,K,2) 592 DR3 = WA2(I-2)*CC(I-1,K,3) + WA2(I-1)*CC(I,K,3) 593 DI3 = WA2(I-2)*CC(I,K,3) - WA2(I-1)*CC(I-1,K,3) 594 DR4 = WA3(I-2)*CC(I-1,K,4) + WA3(I-1)*CC(I,K,4) 595 DI4 = WA3(I-2)*CC(I,K,4) - WA3(I-1)*CC(I-1,K,4) 596 DR5 = WA4(I-2)*CC(I-1,K,5) + WA4(I-1)*CC(I,K,5) 597 DI5 = WA4(I-2)*CC(I,K,5) - WA4(I-1)*CC(I-1,K,5) 598 CR2 = DR2 + DR5 599 CI5 = DR5 - DR2 600 CR5 = DI2 - DI5 601 CI2 = DI2 + DI5 602 CR3 = DR3 + DR4 603 CI4 = DR4 - DR3 604 CR4 = DI3 - DI4 605 CI3 = DI3 + DI4 606 CH(I-1,1,K) = CC(I-1,K,1) + CR2 + CR3 607 CH(I,1,K) = CC(I,K,1) + CI2 + CI3 608 TR2 = CC(I-1,K,1) + TR11*CR2 + TR12*CR3 609 TI2 = CC(I,K,1) + TR11*CI2 + TR12*CI3 610 TR3 = CC(I-1,K,1) + TR12*CR2 + TR11*CR3 611 TI3 = CC(I,K,1) + TR12*CI2 + TR11*CI3 612 TR5 = TI11*CR5 + TI12*CR4 613 TI5 = TI11*CI5 + TI12*CI4 614 TR4 = TI12*CR5 - TI11*CR4 615 TI4 = TI12*CI5 - TI11*CI4 616 CH(I-1,3,K) = TR2 + TR5 617 CH(IC-1,2,K) = TR2 - TR5 618 CH(I,3,K) = TI2 + TI5 619 CH(IC,2,K) = TI5 - TI2 620 CH(I-1,5,K) = TR3 + TR4 621 CH(IC-1,4,K) = TR3 - TR4 622 CH(I,5,K) = TI3 + TI4 623 CH(IC,4,K) = TI4 - TI3 624 102 CONTINUE 625 103 CONTINUE 626 RETURN 627 END 628CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 629C 630C File: radfg.f 631C 632C Library: FFTPACK 4.1 633C 634C Author: Paul N. Swarztrauber 635C National Center for Atmospheric Research 636C PO 3000, Boulder, Colorado 637C 638C Date: Wed Mar 29 18:31:13 MST 1995 639C 640C Description: Lower-level auxiliary routine 641C 642 SUBROUTINE RADFG(IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) 643 DOUBLE PRECISION TPI 644 DOUBLE PRECISION PIMACH 645 DOUBLE PRECISION DUM 646 DOUBLE PRECISION ARG 647 DOUBLE PRECISION DCP 648 DOUBLE PRECISION DSP 649 DOUBLE PRECISION AR1 650 DOUBLE PRECISION AI1 651 DOUBLE PRECISION AR1H 652 DOUBLE PRECISION DC2 653 DOUBLE PRECISION DS2 654 DOUBLE PRECISION AR2 655 DOUBLE PRECISION AI2 656 DOUBLE PRECISION AR2H 657 DOUBLE PRECISION CH(IDO,L1,IP),CC(IDO,IP,L1),C1(IDO,L1,IP), 658 + C2(IDL1,IP),CH2(IDL1,IP),WA(*) 659C 660C FFTPACK 5.0 auxiliary routine 661C 662 TPI = 2.0D0*PIMACH(DUM) 663 ARG = TPI/DBLE(IP) 664 DCP = COS(ARG) 665 DSP = SIN(ARG) 666 IPPH = (IP+1)/2 667 IPP2 = IP + 2 668 IDP2 = IDO + 2 669 NBD = (IDO-1)/2 670 IF (IDO.EQ.1) GO TO 119 671 DO 101 IK = 1,IDL1 672 CH2(IK,1) = C2(IK,1) 673 101 CONTINUE 674 DO 103 J = 2,IP 675 DO 102 K = 1,L1 676 CH(1,K,J) = C1(1,K,J) 677 102 CONTINUE 678 103 CONTINUE 679 IF (NBD.GT.L1) GO TO 107 680 IS = -IDO 681 DO 106 J = 2,IP 682 IS = IS + IDO 683 IDIJ = IS 684 DO 105 I = 3,IDO,2 685 IDIJ = IDIJ + 2 686 DO 104 K = 1,L1 687 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J) + 688 + WA(IDIJ)*C1(I,K,J) 689 CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J) - 690 + WA(IDIJ)*C1(I-1,K,J) 691 104 CONTINUE 692 105 CONTINUE 693 106 CONTINUE 694 GO TO 111 695 107 IS = -IDO 696 DO 110 J = 2,IP 697 IS = IS + IDO 698 DO 109 K = 1,L1 699 IDIJ = IS 700 DO 108 I = 3,IDO,2 701 IDIJ = IDIJ + 2 702 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J) + 703 + WA(IDIJ)*C1(I,K,J) 704 CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J) - 705 + WA(IDIJ)*C1(I-1,K,J) 706 108 CONTINUE 707 109 CONTINUE 708 110 CONTINUE 709 111 IF (NBD.LT.L1) GO TO 115 710 DO 114 J = 2,IPPH 711 JC = IPP2 - J 712 DO 113 K = 1,L1 713 DO 112 I = 3,IDO,2 714 C1(I-1,K,J) = CH(I-1,K,J) + CH(I-1,K,JC) 715 C1(I-1,K,JC) = CH(I,K,J) - CH(I,K,JC) 716 C1(I,K,J) = CH(I,K,J) + CH(I,K,JC) 717 C1(I,K,JC) = CH(I-1,K,JC) - CH(I-1,K,J) 718 112 CONTINUE 719 113 CONTINUE 720 114 CONTINUE 721 GO TO 121 722 115 DO 118 J = 2,IPPH 723 JC = IPP2 - J 724 DO 117 I = 3,IDO,2 725 DO 116 K = 1,L1 726 C1(I-1,K,J) = CH(I-1,K,J) + CH(I-1,K,JC) 727 C1(I-1,K,JC) = CH(I,K,J) - CH(I,K,JC) 728 C1(I,K,J) = CH(I,K,J) + CH(I,K,JC) 729 C1(I,K,JC) = CH(I-1,K,JC) - CH(I-1,K,J) 730 116 CONTINUE 731 117 CONTINUE 732 118 CONTINUE 733 GO TO 121 734 119 DO 120 IK = 1,IDL1 735 C2(IK,1) = CH2(IK,1) 736 120 CONTINUE 737 121 DO 123 J = 2,IPPH 738 JC = IPP2 - J 739 DO 122 K = 1,L1 740 C1(1,K,J) = CH(1,K,J) + CH(1,K,JC) 741 C1(1,K,JC) = CH(1,K,JC) - CH(1,K,J) 742 122 CONTINUE 743 123 CONTINUE 744C 745 AR1 = 1.D0 746 AI1 = 0.D0 747 DO 127 L = 2,IPPH 748 LC = IPP2 - L 749 AR1H = DCP*AR1 - DSP*AI1 750 AI1 = DCP*AI1 + DSP*AR1 751 AR1 = AR1H 752 DO 124 IK = 1,IDL1 753 CH2(IK,L) = C2(IK,1) + AR1*C2(IK,2) 754 CH2(IK,LC) = AI1*C2(IK,IP) 755 124 CONTINUE 756 DC2 = AR1 757 DS2 = AI1 758 AR2 = AR1 759 AI2 = AI1 760 DO 126 J = 3,IPPH 761 JC = IPP2 - J 762 AR2H = DC2*AR2 - DS2*AI2 763 AI2 = DC2*AI2 + DS2*AR2 764 AR2 = AR2H 765 DO 125 IK = 1,IDL1 766 CH2(IK,L) = CH2(IK,L) + AR2*C2(IK,J) 767 CH2(IK,LC) = CH2(IK,LC) + AI2*C2(IK,JC) 768 125 CONTINUE 769 126 CONTINUE 770 127 CONTINUE 771 DO 129 J = 2,IPPH 772 DO 128 IK = 1,IDL1 773 CH2(IK,1) = CH2(IK,1) + C2(IK,J) 774 128 CONTINUE 775 129 CONTINUE 776C 777 IF (IDO.LT.L1) GO TO 132 778 DO 131 K = 1,L1 779 DO 130 I = 1,IDO 780 CC(I,1,K) = CH(I,K,1) 781 130 CONTINUE 782 131 CONTINUE 783 GO TO 135 784 132 DO 134 I = 1,IDO 785 DO 133 K = 1,L1 786 CC(I,1,K) = CH(I,K,1) 787 133 CONTINUE 788 134 CONTINUE 789 135 DO 137 J = 2,IPPH 790 JC = IPP2 - J 791 J2 = J + J 792 DO 136 K = 1,L1 793 CC(IDO,J2-2,K) = CH(1,K,J) 794 CC(1,J2-1,K) = CH(1,K,JC) 795 136 CONTINUE 796 137 CONTINUE 797 IF (IDO.EQ.1) RETURN 798 IF (NBD.LT.L1) GO TO 141 799 DO 140 J = 2,IPPH 800 JC = IPP2 - J 801 J2 = J + J 802 DO 139 K = 1,L1 803 DO 138 I = 3,IDO,2 804 IC = IDP2 - I 805 CC(I-1,J2-1,K) = CH(I-1,K,J) + CH(I-1,K,JC) 806 CC(IC-1,J2-2,K) = CH(I-1,K,J) - CH(I-1,K,JC) 807 CC(I,J2-1,K) = CH(I,K,J) + CH(I,K,JC) 808 CC(IC,J2-2,K) = CH(I,K,JC) - CH(I,K,J) 809 138 CONTINUE 810 139 CONTINUE 811 140 CONTINUE 812 RETURN 813 141 DO 144 J = 2,IPPH 814 JC = IPP2 - J 815 J2 = J + J 816 DO 143 I = 3,IDO,2 817 IC = IDP2 - I 818 DO 142 K = 1,L1 819 CC(I-1,J2-1,K) = CH(I-1,K,J) + CH(I-1,K,JC) 820 CC(IC-1,J2-2,K) = CH(I-1,K,J) - CH(I-1,K,JC) 821 CC(I,J2-1,K) = CH(I,K,J) + CH(I,K,JC) 822 CC(IC,J2-2,K) = CH(I,K,JC) - CH(I,K,J) 823 142 CONTINUE 824 143 CONTINUE 825 144 CONTINUE 826 RETURN 827 END 828