Changeset 4559 for palm/trunk/SOURCE/temperton_fft_mod.f90
- Timestamp:
- Jun 11, 2020 8:51:48 AM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/temperton_fft_mod.f90 (modified) (163 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/temperton_fft_mod.f90
r4370 r4559 1 1 !> @file temperton_fft_mod.f90 2 !------------------------------------------------------------------------------! 2 !--------------------------------------------------------------------------------------------------! 3 ! This file is part of the PALM model system. 4 ! 5 ! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General 6 ! Public License as published by the Free Software Foundation, either version 3 of the License, or 7 ! (at your option) any later version. 8 ! 9 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the 10 ! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 11 ! Public License for more details. 12 ! 13 ! You should have received a copy of the GNU General Public License along with PALM. If not, see 14 ! <http://www.gnu.org/licenses/>. 15 ! 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 17 !--------------------------------------------------------------------------------------------------! 18 ! 3 19 ! 4 20 ! Current revisions: … … 9 25 ! ----------------- 10 26 ! $Id$ 27 ! File re-formatted to follow the PALM coding standard 28 ! 29 ! 4370 2020-01-10 14:00:44Z raasch 11 30 ! unused variables removed 12 ! 31 ! 13 32 ! 4366 2020-01-09 08:12:43Z raasch 14 33 ! vectorized routines added 15 ! 34 ! 16 35 ! 4182 2019-08-22 15:20:23Z scharf 17 36 ! Corrected "Former revisions" section 18 ! 37 ! 19 38 ! 3725 2019-02-07 10:11:02Z raasch 20 39 ! GOTO statements replaced, file re-formatted corresponding to coding standards … … 27 46 ! ------------ 28 47 !> Fast Fourier transformation developed by Clive Temperton, ECMWF. 29 !------------------------------------------------------------------------------ !48 !--------------------------------------------------------------------------------------------------! 30 49 MODULE temperton_fft 31 50 … … 57 76 CONTAINS 58 77 59 !------------------------------------------------------------------------------ !78 !--------------------------------------------------------------------------------------------------! 60 79 ! Description: 61 80 ! ------------ 62 81 !> Calls Fortran-versions of fft's. 63 !> 82 !> 64 83 !> Method: 65 !> 66 !> Subroutine 'fft991cy' - multiple fast real periodic transform 67 !> supersedes previous routine'fft991cy'.68 !> 69 !> Real transform of length n performed by removing redundant 70 !> operations from complex transform oflength n.71 !> 84 !> 85 !> Subroutine 'fft991cy' - multiple fast real periodic transform supersedes previous routine 86 !> 'fft991cy'. 87 !> 88 !> Real transform of length n performed by removing redundant operations from complex transform of 89 !> length n. 90 !> 72 91 !> a is the array containing input & output data. 73 92 !> work is an area of size (n+1)*min(lot,nfft). … … 81 100 !> isign = +1 for transform from spectral to gridpoint 82 101 !> = -1 for transform from gridpoint to spectral. 83 !> 102 !> 84 103 !> ordering of coefficients: 85 104 !> a(0),b(0),a(1),b(1),a(2),b(2),.,a(n/2),b(n/2) 86 105 !> where b(0)=b(n/2)=0; (n+2) locations required. 87 !> 106 !> 88 107 !> ordering of data: 89 108 !> x(0),x(1),x(2),.,x(n-1), 0 , 0 ; (n+2) locations required. 90 !> 91 !> Vectorization is achieved on cray by doing the transforms 92 !> in parallel. 93 !> 109 !> 110 !> Vectorization is achieved on cray by doing the transforms in parallel. 111 !> 94 112 !> n must be composed of factors 2,3 & 5 but does not have to be even. 95 !> 113 !> 96 114 !> definition of transforms: 97 !> 115 !> 98 116 !> isign=+1: x(j)=sum(k=0,.,n-1)(c(k)*exp(2*i*j*k*pi/n)) 99 117 !> where c(k)=a(k)+i*b(k) and c(n-k)=a(k)-i*b(k) 100 !> 118 !> 101 119 !> isign=-1: a(k)=(1/n)*sum(j=0,.,n-1)(x(j)*cos(2*j*k*pi/n)) 102 120 !> b(k)=-(1/n)*sum(j=0,.,n-1)(x(j)*sin(2*j*k*pi/n)) 103 !> 121 !> 104 122 !> calls fortran-versions of fft's !!! 105 123 !> dimension a(n),work(n),trigs(n),ifax(1) 106 !------------------------------------------------------------------------------ !124 !--------------------------------------------------------------------------------------------------! 107 125 SUBROUTINE fft991cy( a, work, trigs, ifax, inc, jump, n, lot, isign ) 108 126 … … 111 129 IMPLICIT NONE 112 130 113 INTEGER(iwp) :: inc !<114 INTEGER(iwp) :: isign !<115 INTEGER(iwp) :: jump !<116 INTEGER(iwp) :: lot !<117 INTEGER(iwp) :: n !<118 119 REAL(wp) :: a(*) !<120 REAL(wp) :: trigs(*) !<121 REAL(wp) :: work(*) !<122 INTEGER(iwp) :: i fax(*) !<123 124 INTEGER(iwp) :: i !<125 INTEGER(iwp) :: i a !<126 INTEGER(iwp) :: i base !<127 INTEGER(iwp) :: i err !<128 INTEGER(iwp) :: i fac !<129 INTEGER(iwp) :: i go !<130 INTEGER(iwp) :: ii !<131 INTEGER(iwp) :: istart !<132 INTEGER(iwp) :: ix !<133 INTEGER(iwp) :: iz !<134 INTEGER(iwp) :: j !<135 INTEGER(iwp) :: jbase !<136 INTEGER(iwp) :: jj !<137 INTEGER(iwp) :: k !<138 INTEGER(iwp) :: la !<139 INTEGER(iwp) :: n b !<140 INTEGER(iwp) :: nblox !< 141 INTEGER(iwp) :: nfax !<142 INTEGER(iwp) :: nvex !<143 INTEGER(iwp) :: nx !<131 INTEGER(iwp) :: inc !< 132 INTEGER(iwp) :: isign !< 133 INTEGER(iwp) :: jump !< 134 INTEGER(iwp) :: lot !< 135 INTEGER(iwp) :: n !< 136 137 INTEGER(iwp) :: i !< 138 INTEGER(iwp) :: ia !< 139 INTEGER(iwp) :: ibase !< 140 INTEGER(iwp) :: ierr !< 141 INTEGER(iwp) :: ifac !< 142 INTEGER(iwp) :: igo !< 143 INTEGER(iwp) :: ii !< 144 INTEGER(iwp) :: ifax(*) !< 145 INTEGER(iwp) :: istart !< 146 INTEGER(iwp) :: ix !< 147 INTEGER(iwp) :: iz !< 148 INTEGER(iwp) :: j !< 149 INTEGER(iwp) :: jbase !< 150 INTEGER(iwp) :: jj !< 151 INTEGER(iwp) :: k !< 152 INTEGER(iwp) :: la !< 153 INTEGER(iwp) :: nb !< 154 INTEGER(iwp) :: nblox !< 155 INTEGER(iwp) :: nfax !< 156 INTEGER(iwp) :: nvex !< 157 INTEGER(iwp) :: nx !< 158 159 REAL(wp) :: a(*) !< 160 REAL(wp) :: trigs(*) !< 161 REAL(wp) :: work(*) !< 144 162 145 163 … … 147 165 nfax = ifax(1) 148 166 nx = n + 1 149 IF ( MOD( n,2) == 1 ) nx = n167 IF ( MOD( n, 2) == 1 ) nx = n 150 168 nblox = 1 + (lot-1) / nfft 151 nvex = lot - ( nblox-1) * nfft169 nvex = lot - ( nblox - 1 ) * nfft 152 170 153 171 … … 160 178 161 179 ia = istart 162 i = istart180 i = istart 163 181 !DIR$ IVDEP 164 182 !CDIR NODEP … … 166 184 DO j = 1, nvex 167 185 a(i+inc) = 0.5_wp * a(i) 168 i = i + jump186 i = i + jump 169 187 ENDDO 170 188 IF ( MOD(n,2) /= 1 ) THEN 171 189 i = istart + n * inc 172 DO j = 1, nvex190 DO j = 1, nvex 173 191 a(i) = 0.5_wp * a(i) 174 i = i + jump192 i = i + jump 175 193 ENDDO 176 194 ENDIF … … 211 229 ! 212 230 !-- If necessary, copy results back to a 213 IF ( MOD( nfax,2) /= 0 ) THEN231 IF ( MOD( nfax, 2 ) /= 0 ) THEN 214 232 ibase = 1 215 233 jbase = ia … … 219 237 DO ii = 1, n 220 238 a(j) = work(i) 221 i = i + 1222 j = j + inc239 i = i + 1 240 j = j + inc 223 241 ENDDO 224 242 ibase = ibase + nx … … 228 246 ! 229 247 !-- Fill in zeros at end 230 ix = istart + n *inc248 ix = istart + n * inc 231 249 !DIR$ IVDEP 232 250 !CDIR NODEP 233 251 !OCL NOVREC 234 252 DO j = 1, nvex 235 a(ix) = 0.0_wp253 a(ix) = 0.0_wp 236 254 a(ix+inc) = 0.0_wp 237 ix = ix + jump255 ix = ix + jump 238 256 ENDDO 239 istart = istart + nvex *jump257 istart = istart + nvex * jump 240 258 nvex = nfft 241 259 … … 254 272 255 273 ifac = ifax(nfax+2-k) 256 la = la / ifac274 la = la / ifac 257 275 ierr = -1 258 276 IF ( igo /= -1 ) THEN … … 291 309 DO ii = 1, n 292 310 a(j) = work(i) 293 i = i + 1294 j = j + inc311 i = i + 1 312 j = j + inc 295 313 ENDDO 296 314 ibase = ibase + nx … … 305 323 !OCL NOVREC 306 324 DO j = 1, nvex 307 a(ix) = a(ix+inc)325 a(ix) = a(ix+inc) 308 326 a(ix+inc) = 0.0_wp 309 ix = ix + jump327 ix = ix + jump 310 328 ENDDO 311 329 IF ( MOD(n,2) /= 1 ) THEN … … 313 331 DO j = 1, nvex 314 332 a(iz) = 0.0_wp 315 iz = iz + jump333 iz = iz + jump 316 334 ENDDO 317 335 ENDIF … … 325 343 END SUBROUTINE fft991cy 326 344 327 !------------------------------------------------------------------------------ !345 !--------------------------------------------------------------------------------------------------! 328 346 ! Description: 329 347 ! ------------ 330 !> Performs one pass through data as part of 331 !> multiple real fft (fourier analysis) routine. 332 !> 348 !> Performs one pass through data as part of multiple real fft (fourier analysis) routine. 349 !> 333 350 !> Method: 334 !> 351 !> 335 352 !> a is first real input vector 336 353 !> equivalence b(1) with a(ifac*la*inc1+1) … … 351 368 !> 2 - ifac not catered for 352 369 !> 3 - ifac only catered for if la=n/ifac 353 !------------------------------------------------------------------------------ !370 !--------------------------------------------------------------------------------------------------! 354 371 SUBROUTINE qpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr ) 355 372 356 373 USE kinds 357 374 358 IMPLICIT NONE 359 360 INTEGER(iwp) :: ierr !< 361 INTEGER(iwp) :: ifac !< 362 INTEGER(iwp) :: inc1 !< 363 INTEGER(iwp) :: inc2 !< 364 INTEGER(iwp) :: inc3 !< 365 INTEGER(iwp) :: inc4 !< 366 INTEGER(iwp) :: la !< 367 INTEGER(iwp) :: lot !< 368 INTEGER(iwp) :: n !< 375 IMPLICIT NONE 376 377 INTEGER(iwp) :: ierr !< 378 INTEGER(iwp) :: ifac !< 379 INTEGER(iwp) :: inc1 !< 380 INTEGER(iwp) :: inc2 !< 381 INTEGER(iwp) :: inc3 !< 382 INTEGER(iwp) :: inc4 !< 383 INTEGER(iwp) :: la !< 384 INTEGER(iwp) :: lot !< 385 INTEGER(iwp) :: n !< 386 387 INTEGER(iwp) :: i !< 388 INTEGER(iwp) :: ia !< 389 INTEGER(iwp) :: ib !< 390 INTEGER(iwp) :: ibase !< 391 INTEGER(iwp) :: ic !< 392 INTEGER(iwp) :: id !< 393 INTEGER(iwp) :: ie !< 394 INTEGER(iwp) :: if !< 395 INTEGER(iwp) :: ig !< 396 INTEGER(iwp) :: igo !< 397 INTEGER(iwp) :: ih !< 398 INTEGER(iwp) :: iink !< 399 INTEGER(iwp) :: ijk !< 400 INTEGER(iwp) :: ijump !< 401 INTEGER(iwp) :: j !< 402 INTEGER(iwp) :: ja !< 403 INTEGER(iwp) :: jb !< 404 INTEGER(iwp) :: jbase !< 405 INTEGER(iwp) :: jc !< 406 INTEGER(iwp) :: jd !< 407 INTEGER(iwp) :: je !< 408 INTEGER(iwp) :: jf !< 409 INTEGER(iwp) :: jink !< 410 INTEGER(iwp) :: k !< 411 INTEGER(iwp) :: kb !< 412 INTEGER(iwp) :: kc !< 413 INTEGER(iwp) :: kd !< 414 INTEGER(iwp) :: ke !< 415 INTEGER(iwp) :: kf !< 416 INTEGER(iwp) :: kstop !< 417 INTEGER(iwp) :: l !< 418 INTEGER(iwp) :: m !< 369 419 370 420 ! 371 421 !-- Arrays are dimensioned with n 372 REAL(wp) :: a(*) !< 373 REAL(wp) :: b(*) !< 374 REAL(wp) :: c(*) !< 375 REAL(wp) :: d(*) !< 376 REAL(wp) :: trigs(*) !< 377 378 REAL(wp) :: a0 !< 379 REAL(wp) :: a1 !< 380 REAL(wp) :: a10 !< 381 REAL(wp) :: a11 !< 382 REAL(wp) :: a2 !< 383 REAL(wp) :: a20 !< 384 REAL(wp) :: a21 !< 385 REAL(wp) :: a3 !< 386 REAL(wp) :: a4 !< 387 REAL(wp) :: a5 !< 388 REAL(wp) :: a6 !< 389 REAL(wp) :: b0 !< 390 REAL(wp) :: b1 !< 391 REAL(wp) :: b10 !< 392 REAL(wp) :: b11 !< 393 REAL(wp) :: b2 !< 394 REAL(wp) :: b20 !< 395 REAL(wp) :: b21 !< 396 REAL(wp) :: b3 !< 397 REAL(wp) :: b4 !< 398 REAL(wp) :: b5 !< 399 REAL(wp) :: b6 !< 400 REAL(wp) :: c1 !< 401 REAL(wp) :: c2 !< 402 REAL(wp) :: c3 !< 403 REAL(wp) :: c4 !< 404 REAL(wp) :: c5 !< 405 REAL(wp) :: qrt5 !< 406 REAL(wp) :: s1 !< 407 REAL(wp) :: s2 !< 408 REAL(wp) :: s3 !< 409 REAL(wp) :: s4 !< 410 REAL(wp) :: s5 !< 411 REAL(wp) :: sin36 !< 412 REAL(wp) :: sin45 !< 413 REAL(wp) :: sin60 !< 414 REAL(wp) :: sin72 !< 415 REAL(wp) :: z !< 416 REAL(wp) :: zqrt5 !< 417 REAL(wp) :: zsin36 !< 418 REAL(wp) :: zsin45 !< 419 REAL(wp) :: zsin60 !< 420 REAL(wp) :: zsin72 !< 421 422 INTEGER(iwp) :: i !< 423 INTEGER(iwp) :: ia !< 424 INTEGER(iwp) :: ib !< 425 INTEGER(iwp) :: ibase !< 426 INTEGER(iwp) :: ic !< 427 INTEGER(iwp) :: id !< 428 INTEGER(iwp) :: ie !< 429 INTEGER(iwp) :: if !< 430 INTEGER(iwp) :: ig !< 431 INTEGER(iwp) :: igo !< 432 INTEGER(iwp) :: ih !< 433 INTEGER(iwp) :: iink !< 434 INTEGER(iwp) :: ijk !< 435 INTEGER(iwp) :: ijump !< 436 INTEGER(iwp) :: j !< 437 INTEGER(iwp) :: ja !< 438 INTEGER(iwp) :: jb !< 439 INTEGER(iwp) :: jbase !< 440 INTEGER(iwp) :: jc !< 441 INTEGER(iwp) :: jd !< 442 INTEGER(iwp) :: je !< 443 INTEGER(iwp) :: jf !< 444 INTEGER(iwp) :: jink !< 445 INTEGER(iwp) :: k !< 446 INTEGER(iwp) :: kb !< 447 INTEGER(iwp) :: kc !< 448 INTEGER(iwp) :: kd !< 449 INTEGER(iwp) :: ke !< 450 INTEGER(iwp) :: kf !< 451 INTEGER(iwp) :: kstop !< 452 INTEGER(iwp) :: l !< 453 INTEGER(iwp) :: m !< 422 REAL(wp) :: a(*) !< 423 REAL(wp) :: b(*) !< 424 REAL(wp) :: c(*) !< 425 REAL(wp) :: d(*) !< 426 REAL(wp) :: trigs(*) !< 427 428 REAL(wp) :: a0 !< 429 REAL(wp) :: a1 !< 430 REAL(wp) :: a10 !< 431 REAL(wp) :: a11 !< 432 REAL(wp) :: a2 !< 433 REAL(wp) :: a20 !< 434 REAL(wp) :: a21 !< 435 REAL(wp) :: a3 !< 436 REAL(wp) :: a4 !< 437 REAL(wp) :: a5 !< 438 REAL(wp) :: a6 !< 439 REAL(wp) :: b0 !< 440 REAL(wp) :: b1 !< 441 REAL(wp) :: b10 !< 442 REAL(wp) :: b11 !< 443 REAL(wp) :: b2 !< 444 REAL(wp) :: b20 !< 445 REAL(wp) :: b21 !< 446 REAL(wp) :: b3 !< 447 REAL(wp) :: b4 !< 448 REAL(wp) :: b5 !< 449 REAL(wp) :: b6 !< 450 REAL(wp) :: c1 !< 451 REAL(wp) :: c2 !< 452 REAL(wp) :: c3 !< 453 REAL(wp) :: c4 !< 454 REAL(wp) :: c5 !< 455 REAL(wp) :: qrt5 !< 456 REAL(wp) :: s1 !< 457 REAL(wp) :: s2 !< 458 REAL(wp) :: s3 !< 459 REAL(wp) :: s4 !< 460 REAL(wp) :: s5 !< 461 REAL(wp) :: sin36 !< 462 REAL(wp) :: sin45 !< 463 REAL(wp) :: sin60 !< 464 REAL(wp) :: sin72 !< 465 REAL(wp) :: z !< 466 REAL(wp) :: zqrt5 !< 467 REAL(wp) :: zsin36 !< 468 REAL(wp) :: zsin45 !< 469 REAL(wp) :: zsin60 !< 470 REAL(wp) :: zsin72 !< 454 471 455 472 … … 468 485 iink = la * inc1 469 486 jink = la * inc2 470 ijump = (ifac -1) * iink471 kstop = ( n -ifac ) / ( 2*ifac )487 ijump = (ifac - 1) * iink 488 kstop = ( n - ifac ) / ( 2 * ifac ) 472 489 473 490 ibase = 0 … … 490 507 ib = ia + iink 491 508 ja = 1 492 jb = ja + ( 2*m-la) * inc2509 jb = ja + ( 2 * m - la ) * inc2 493 510 494 511 IF ( la /= m ) THEN … … 510 527 ENDDO 511 528 512 ja = ja + jink513 jink = 2 * jink514 jb = jb - jink529 ja = ja + jink 530 jink = 2 * jink 531 jb = jb - jink 515 532 ibase = ibase + ijump 516 533 ijump = 2 * ijump + iink … … 530 547 !OCL NOVREC 531 548 DO ijk = 1, lot 532 c(ja+j) = a(ia+i) + ( c1*a(ib+i)+s1*b(ib+i))533 c(jb+j) = a(ia+i) - ( c1*a(ib+i)+s1*b(ib+i))534 d(ja+j) = ( c1*b(ib+i)-s1*a(ib+i)) + b(ia+i)535 d(jb+j) = ( c1*b(ib+i)-s1*a(ib+i)) - b(ia+i)549 c(ja+j) = a(ia+i) + ( c1 * a(ib+i) + s1 * b(ib+i) ) 550 c(jb+j) = a(ia+i) - ( c1 * a(ib+i) + s1 * b(ib+i) ) 551 d(ja+j) = ( c1 * b(ib+i) - s1 * a(ib+i) ) + b(ia+i) 552 d(jb+j) = ( c1 * b(ib+i) - s1 * a(ib+i) ) - b(ia+i) 536 553 i = i + inc3 537 554 j = j + inc4 … … 562 579 i = i + inc3 563 580 j = j + inc4 564 565 581 ENDDO 566 582 ibase = ibase + inc1 … … 570 586 ELSE 571 587 572 z = 1.0_wp /REAL(n,KIND=wp)588 z = 1.0_wp / REAL( n, KIND = wp ) 573 589 DO l = 1, la 574 590 i = ibase … … 578 594 !OCL NOVREC 579 595 DO ijk = 1, lot 580 c(ja+j) = z *(a(ia+i)+a(ib+i))581 c(jb+j) = z *(a(ia+i)-a(ib+i))596 c(ja+j) = z * ( a(ia+i) + a(ib+i) ) 597 c(jb+j) = z * ( a(ia+i) - a(ib+i) ) 582 598 i = i + inc3 583 599 j = j + inc4 … … 597 613 ic = ib + iink 598 614 ja = 1 599 jb = ja + ( 2*m-la) * inc2615 jb = ja + ( 2 * m - la) * inc2 600 616 jc = jb 601 617 … … 609 625 !OCL NOVREC 610 626 DO ijk = 1, lot 611 c(ja+j) = a(ia+i) + ( a(ib+i)+a(ic+i))612 c(jb+j) = a(ia+i) - 0.5_wp *(a(ib+i)+a(ic+i))613 d(jb+j) = sin60 *(a(ic+i)-a(ib+i))627 c(ja+j) = a(ia+i) + ( a(ib+i) + a(ic+i) ) 628 c(jb+j) = a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) 629 d(jb+j) = sin60 * ( a(ic+i) - a(ib+i) ) 614 630 i = i + inc3 615 631 j = j + inc4 … … 644 660 !OCL NOVREC 645 661 DO ijk = 1, lot 646 a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c2*a(ic+i)+s2*b(ic+i)) 647 b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c2*b(ic+i)-s2*a(ic+i)) 648 a2 = a(ia+i) - 0.5_wp*a1 649 b2 = b(ia+i) - 0.5_wp*b1 650 a3 = sin60*((c1*a(ib+i)+s1*b(ib+i))-(c2*a(ic+i)+s2*b(ic+i))) 651 b3 = sin60*((c1*b(ib+i)-s1*a(ib+i))-(c2*b(ic+i)-s2*a(ic+i))) 662 a1 = ( c1 * a(ib+i) + s1 * b(ib+i) ) + ( c2 * a(ic+i) + s2 * b(ic+i) ) 663 b1 = ( c1 * b(ib+i) - s1 * a(ib+i) ) + ( c2 * b(ic+i) - s2 * a(ic+i) ) 664 a2 = a(ia+i) - 0.5_wp * a1 665 b2 = b(ia+i) - 0.5_wp * b1 666 a3 = sin60 * ( ( c1 * a(ib+i) + s1 * b(ib+i) ) - ( c2 * a(ic+i) + s2 & 667 * b(ic+i) ) ) 668 b3 = sin60 * ( ( c1 * b(ib+i) - s1 * a(ib+i) ) - ( c2 * b(ic+i) - s2 & 669 * a(ic+i) ) ) 652 670 c(ja+j) = a(ia+i) + a1 653 671 d(ja+j) = b(ia+i) + b1 … … 655 673 d(jb+j) = b2 - a3 656 674 c(jc+j) = a2 - b3 657 d(jc+j) = - (b2+a3)675 d(jc+j) = - ( b2 + a3 ) 658 676 i = i + inc3 659 677 j = j + inc4 … … 680 698 !OCL NOVREC 681 699 DO ijk = 1, lot 682 c(ja+j) = a(ia+i) + 0.5_wp *(a(ib+i)-a(ic+i))683 d(ja+j) = -sin60 *(a(ib+i)+a(ic+i))684 c(jb+j) = a(ia+i) - ( a(ib+i)-a(ic+i))700 c(ja+j) = a(ia+i) + 0.5_wp * ( a(ib+i) - a(ic+i) ) 701 d(ja+j) = -sin60 * ( a(ib+i) + a(ic+i) ) 702 c(jb+j) = a(ia+i) - ( a(ib+i) - a(ic+i) ) 685 703 i = i + inc3 686 704 j = j + inc4 … … 692 710 ELSE 693 711 694 z = 1.0_wp / REAL( n, KIND =wp )695 zsin60 = z *sin60712 z = 1.0_wp / REAL( n, KIND = wp ) 713 zsin60 = z * sin60 696 714 DO l = 1, la 697 715 i = ibase … … 701 719 !OCL NOVREC 702 720 DO ijk = 1, lot 703 c(ja+j) = z *(a(ia+i)+(a(ib+i)+a(ic+i)))704 c(jb+j) = z *(a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))705 d(jb+j) = zsin60 *(a(ic+i)-a(ib+i))721 c(ja+j) = z * ( a(ia+i) + ( a (ib+i) + a(ic+i) ) ) 722 c(jb+j) = z * ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) 723 d(jb+j) = zsin60 * ( a(ic+i) -a(ib+i) ) 706 724 i = i + inc3 707 725 j = j + inc4 … … 722 740 id = ic + iink 723 741 ja = 1 724 jb = ja + ( 2*m-la) * inc2725 jc = jb + 2 *m*inc2742 jb = ja + ( 2 * m - la ) * inc2 743 jc = jb + 2 * m * inc2 726 744 jd = jb 727 745 … … 736 754 !OCL NOVREC 737 755 DO ijk = 1, lot 738 c(ja+j) = ( a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))739 c(jc+j) = ( a(ia+i)+a(ic+i)) - (a(ib+i)+a(id+i))756 c(ja+j) = ( a(ia+i) + a(ic+i) ) + ( a(ib+i) + a(id+i) ) 757 c(jc+j) = ( a(ia+i) + a(ic+i) ) - ( a(ib+i) + a(id+i) ) 740 758 c(jb+j) = a(ia+i) - a(ic+i) 741 759 d(jb+j) = a(id+i) - a(ib+i) … … 747 765 ENDDO 748 766 749 ja = ja + jink750 jink = 2 * jink751 jb = jb + jink752 jc = jc - jink753 jd = jd - jink767 ja = ja + jink 768 jink = 2 * jink 769 jb = jb + jink 770 jc = jc - jink 771 jd = jd - jink 754 772 ibase = ibase + ijump 755 773 ijump = 2 * ijump + iink … … 775 793 !OCL NOVREC 776 794 DO ijk = 1, lot 777 a0 = a(ia+i) + ( c2*a(ic+i)+s2*b(ic+i))778 a2 = a(ia+i) - ( c2*a(ic+i)+s2*b(ic+i))779 a1 = ( c1*a(ib+i)+s1*b(ib+i)) + (c3*a(id+i)+s3*b(id+i))780 a3 = ( c1*a(ib+i)+s1*b(ib+i)) - (c3*a(id+i)+s3*b(id+i))781 b0 = b(ia+i) + ( c2*b(ic+i)-s2*a(ic+i))782 b2 = b(ia+i) - ( c2*b(ic+i)-s2*a(ic+i))783 b1 = ( c1*b(ib+i)-s1*a(ib+i)) + (c3*b(id+i)-s3*a(id+i))784 b3 = ( c1*b(ib+i)-s1*a(ib+i)) - (c3*b(id+i)-s3*a(id+i))795 a0 = a(ia+i) + ( c2 * a(ic+i) + s2 * b(ic+i) ) 796 a2 = a(ia+i) - ( c2 * a(ic+i) + s2 * b(ic+i) ) 797 a1 = ( c1 * a(ib+i) + s1 * b(ib+i) ) + ( c3 * a(id+i) + s3 * b(id+i) ) 798 a3 = ( c1 * a(ib+i) + s1 * b(ib+i) ) - ( c3 * a(id+i) + s3 * b(id+i) ) 799 b0 = b(ia+i) + ( c2 * b(ic+i) - s2 * a(ic+i) ) 800 b2 = b(ia+i) - ( c2 * b(ic+i) - s2 * a(ic+i) ) 801 b1 = ( c1 * b(ib+i) - s1 * a(ib+i) ) + ( c3 * b(id+i) - s3 * a(id+i) ) 802 b3 = ( c1 * b(ib+i) - s1 * a(ib+i) ) - ( c3 * b(id+i) - s3 * a(id+i) ) 785 803 c(ja+j) = a0 + a1 786 804 c(jc+j) = a0 - a1 … … 790 808 c(jd+j) = a2 - b3 791 809 d(jb+j) = b2 - a3 792 d(jd+j) = - (b2+a3)810 d(jd+j) = - ( b2 + a3 ) 793 811 i = i + inc3 794 812 j = j + inc4 … … 817 835 !OCL NOVREC 818 836 DO ijk = 1, lot 819 c(ja+j) = a(ia+i) + sin45 *(a(ib+i)-a(id+i))820 c(jb+j) = a(ia+i) - sin45 *(a(ib+i)-a(id+i))821 d(ja+j) = - a(ic+i) - sin45*(a(ib+i)+a(id+i))822 d(jb+j) = a(ic+i) - sin45 *(a(ib+i)+a(id+i))837 c(ja+j) = a(ia+i) + sin45 * ( a(ib+i) - a(id+i) ) 838 c(jb+j) = a(ia+i) - sin45 * ( a(ib+i) - a(id+i) ) 839 d(ja+j) = - a(ic+i) - sin45 * ( a(ib+i) + a(id+i) ) 840 d(jb+j) = a(ic+i) - sin45 * ( a(ib+i) + a(id+i) ) 823 841 i = i + inc3 824 842 j = j + inc4 … … 830 848 ELSE 831 849 832 z = 1.0_wp / REAL( n, KIND =wp )850 z = 1.0_wp / REAL( n, KIND = wp ) 833 851 DO l = 1, la 834 852 i = ibase … … 838 856 !OCL NOVREC 839 857 DO ijk = 1, lot 840 c(ja+j) = z *((a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i)))841 c(jc+j) = z *((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i)))842 c(jb+j) = z *(a(ia+i)-a(ic+i))843 d(jb+j) = z *(a(id+i)-a(ib+i))858 c(ja+j) = z * ( ( a(ia+i) + a(ic+i) ) + ( a(ib+i) + a(id+i) ) ) 859 c(jc+j) = z * ( ( a(ia+i) + a(ic+i) ) - (a (ib+i) + a(id+i) ) ) 860 c(jb+j) = z * ( a(ia+i) - a(ic+i) ) 861 d(jb+j) = z * ( a(id+i) - a(ib+i) ) 844 862 i = i + inc3 845 863 j = j + inc4 … … 861 879 ie = id + iink 862 880 ja = 1 863 jb = ja + ( 2*m-la) * inc2864 jc = jb + 2 *m*inc2881 jb = ja + ( 2 * m - la ) * inc2 882 jc = jb + 2 * m * inc2 865 883 jd = jc 866 884 je = jb … … 879 897 a2 = a(ic+i) + a(id+i) 880 898 a4 = a(ic+i) - a(id+i) 881 a5 = a(ia+i) - 0.25_wp *(a1+a2)882 a6 = qrt5 *(a1-a2)883 c(ja+j) = a(ia+i) + ( a1+a2)899 a5 = a(ia+i) - 0.25_wp * ( a1 + a2 ) 900 a6 = qrt5 * ( a1 - a2 ) 901 c(ja+j) = a(ia+i) + ( a1 + a2 ) 884 902 c(jb+j) = a5 + a6 885 903 c(jc+j) = a5 - a6 886 d(jb+j) = - sin72*a3 - sin36*a4887 d(jc+j) = - sin36*a3 + sin72*a4904 d(jb+j) = - sin72 * a3 - sin36 * a4 905 d(jc+j) = - sin36 * a3 + sin72 * a4 888 906 i = i + inc3 889 907 j = j + inc4 … … 925 943 !OCL NOVREC 926 944 DO ijk = 1, lot 927 a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c4*a(ie+i)+s4*b(ie+i))928 a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c4*a(ie+i)+s4*b(ie+i))929 a2 = (c2*a(ic+i)+s2*b(ic+i)) + (c3*a(id+i)+s3*b(id+i))930 a4 = (c2*a(ic+i)+s2*b(ic+i)) - (c3*a(id+i)+s3*b(id+i))931 b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c4*b(ie+i)-s4*a(ie+i))932 b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c4*b(ie+i)-s4*a(ie+i))933 b2 = (c2*b(ic+i)-s2*a(ic+i)) + (c3*b(id+i)-s3*a(id+i))934 b4 = (c2*b(ic+i)-s2*a(ic+i)) - (c3*b(id+i)-s3*a(id+i))935 a5 = a(ia+i) - 0.25_wp*(a1+a2)936 a6 = qrt5*(a1-a2)937 b5 = b(ia+i) - 0.25_wp*(b1+b2)938 b6 = qrt5*(b1-b2)945 a1 = (c1 * a(ib+i) + s1 * b(ib+i) ) + ( c4 * a(ie+i) + s4 * b(ie+i) ) 946 a3 = (c1 * a(ib+i) + s1 * b(ib+i) ) - ( c4 * a(ie+i) + s4 * b(ie+i) ) 947 a2 = (c2 * a(ic+i) + s2 * b(ic+i) ) + ( c3 * a(id+i) + s3 * b(id+i) ) 948 a4 = (c2 * a(ic+i) + s2 * b(ic+i) ) - ( c3 * a(id+i) + s3 * b(id+i) ) 949 b1 = (c1 * b(ib+i) - s1 * a(ib+i) ) + ( c4 * b(ie+i) - s4 * a(ie+i) ) 950 b3 = (c1 * b(ib+i) - s1 * a(ib+i) ) - ( c4 * b(ie+i) - s4 * a(ie+i) ) 951 b2 = (c2 * b(ic+i) - s2 * a(ic+i) ) + ( c3 * b(id+i) - s3 * a(id+i) ) 952 b4 = (c2 * b(ic+i) - s2 * a(ic+i) ) - ( c3 * b(id+i) - s3 * a(id+i) ) 953 a5 = a(ia+i) - 0.25_wp * ( a1 + a2 ) 954 a6 = qrt5 * ( a1 - a2 ) 955 b5 = b(ia+i) - 0.25_wp * ( b1 + b2 ) 956 b6 = qrt5 * ( b1 - b2 ) 939 957 a10 = a5 + a6 940 958 a20 = a5 - a6 941 959 b10 = b5 + b6 942 960 b20 = b5 - b6 943 a11 = sin72 *b3 + sin36*b4944 a21 = sin36 *b3 - sin72*b4945 b11 = sin72 *a3 + sin36*a4946 b21 = sin36 *a3 - sin72*a4947 c(ja+j) = a(ia+i) + ( a1+a2)961 a11 = sin72 * b3 + sin36 * b4 962 a21 = sin36 * b3 - sin72 * b4 963 b11 = sin72 * a3 + sin36 * a4 964 b21 = sin36 * a3 - sin72 * a4 965 c(ja+j) = a(ia+i) + ( a1 + a2 ) 948 966 c(jb+j) = a10 + a11 949 967 c(je+j) = a10 - a11 950 968 c(jc+j) = a20 + a21 951 969 c(jd+j) = a20 - a21 952 d(ja+j) = b(ia+i) + ( b1+b2)970 d(ja+j) = b(ia+i) + ( b1 + b2 ) 953 971 d(jb+j) = b10 - b11 954 d(je+j) = - (b10+b11)972 d(je+j) = - ( b10 + b11 ) 955 973 d(jc+j) = b20 - b21 956 d(jd+j) = - (b20+b21)974 d(jd+j) = - ( b20 + b21 ) 957 975 i = i + inc3 958 976 j = j + inc4 … … 985 1003 a2 = a(ic+i) + a(id+i) 986 1004 a4 = a(ic+i) - a(id+i) 987 a5 = a(ia+i) + 0.25_wp *(a3-a4)988 a6 = qrt5 *(a3+a4)1005 a5 = a(ia+i) + 0.25_wp * ( a3 - a4 ) 1006 a6 = qrt5 * ( a3 + a4 ) 989 1007 c(ja+j) = a5 + a6 990 1008 c(jb+j) = a5 - a6 991 c(jc+j) = a(ia+i) - ( a3-a4)992 d(ja+j) = - sin36*a1 - sin72*a2993 d(jb+j) = - sin72*a1 + sin36*a21009 c(jc+j) = a(ia+i) - ( a3 - a4 ) 1010 d(ja+j) = - sin36 * a1 - sin72 * a2 1011 d(jb+j) = - sin72 * a1 + sin36 * a2 994 1012 i = i + inc3 995 1013 j = j + inc4 … … 1001 1019 ELSE 1002 1020 1003 z = 1.0_wp / REAL( n, KIND =wp )1021 z = 1.0_wp / REAL( n, KIND = wp ) 1004 1022 zqrt5 = z * qrt5 1005 1023 zsin36 = z * sin36 … … 1016 1034 a2 = a(ic+i) + a(id+i) 1017 1035 a4 = a(ic+i) - a(id+i) 1018 a5 = z *(a(ia+i)-0.25_wp*(a1+a2))1019 a6 = zqrt5 *(a1-a2)1020 c(ja+j) = z *(a(ia+i)+(a1+a2))1036 a5 = z * ( a (ia+i) - 0.25_wp * ( a1 + a2 ) ) 1037 a6 = zqrt5 * ( a1 - a2 ) 1038 c(ja+j) = z * ( a(ia+i) + (a1+a2) ) 1021 1039 c(jb+j) = a5 + a6 1022 1040 c(jc+j) = a5 - a6 1023 d(jb+j) = - zsin72*a3 - zsin36*a41024 d(jc+j) = - zsin36*a3 + zsin72*a41041 d(jb+j) = - zsin72 * a3 - zsin36 * a4 1042 d(jc+j) = - zsin36 * a3 + zsin72 * a4 1025 1043 i = i + inc3 1026 1044 j = j + inc4 … … 1043 1061 if = ie + iink 1044 1062 ja = 1 1045 jb = ja + ( 2*m-la) * inc21046 jc = jb + 2 *m*inc21047 jd = jc + 2 *m*inc21063 jb = ja + ( 2 * m - la ) * inc2 1064 jc = jb + 2 * m * inc2 1065 jd = jc + 2 * m * inc2 1048 1066 je = jc 1049 1067 jf = jb … … 1058 1076 !OCL NOVREC 1059 1077 DO ijk = 1, lot 1060 a11 = ( a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))1061 c(ja+j) = ( a(ia+i)+a(id+i)) + a111062 c(jc+j) = ( a(ia+i)+a(id+i)-0.5_wp*a11)1063 d(jc+j) = sin60 *((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))1064 a11 = ( a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))1065 c(jb+j) = ( a(ia+i)-a(id+i)) - 0.5_wp*a111066 d(jb+j) = sin60 *((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))1067 c(jd+j) = ( a(ia+i)-a(id+i)) + a111078 a11 = ( a (ic+i) + a(if+i) ) + ( a(ib+i) + a(ie+i) ) 1079 c(ja+j) = ( a(ia+i) + a(id+i) ) + a11 1080 c(jc+j) = ( a(ia+i) + a(id+i) - 0.5_wp * a11 ) 1081 d(jc+j) = sin60 * ( ( a(ic+i) + a(if+i) ) - ( a(ib+i) + a(ie+i) ) ) 1082 a11 = ( a (ic+i) - a(if+i) ) + ( a(ie+i) - a(ib+i) ) 1083 c(jb+j) = ( a(ia+i) - a(id+i) ) - 0.5_wp * a11 1084 d(jb+j) = sin60 * ( ( a(ie+i) - a(ib+i) ) - ( a(ic+i) - a(if+i) ) ) 1085 c(jd+j) = ( a (ia+i) - a(id+i) ) + a11 1068 1086 i = i + inc3 1069 1087 j = j + inc4 … … 1072 1090 jbase = jbase + inc2 1073 1091 ENDDO 1074 ja = ja + jink1075 jink = 2 * jink1076 jb = jb + jink1077 jc = jc + jink1078 jd = jd - jink1079 je = je - jink1080 jf = jf - jink1092 ja = ja + jink 1093 jink = 2 * jink 1094 jb = jb + jink 1095 jc = jc + jink 1096 jd = jd - jink 1097 je = je - jink 1098 jf = jf - jink 1081 1099 ibase = ibase + ijump 1082 1100 ijump = 2 * ijump + iink … … 1108 1126 !OCL NOVREC 1109 1127 DO ijk = 1, lot 1110 a1 = c1*a(ib+i) + s1*b(ib+i)1111 b1 = c1*b(ib+i) - s1*a(ib+i)1112 a2 = c2*a(ic+i) + s2*b(ic+i)1113 b2 = c2*b(ic+i) - s2*a(ic+i)1114 a3 = c3*a(id+i) + s3*b(id+i)1115 b3 = c3*b(id+i) - s3*a(id+i)1116 a4 = c4*a(ie+i) + s4*b(ie+i)1117 b4 = c4*b(ie+i) - s4*a(ie+i)1118 a5 = c5*a(if+i) + s5*b(if+i)1119 b5 = c5*b(if+i) - s5*a(if+i)1120 a11 = ( a2+a5) + (a1+a4)1121 a20 = ( a(ia+i)+a3) - 0.5_wp*a111122 a21 = sin60 *((a2+a5)-(a1+a4))1123 b11 = ( b2+b5) + (b1+b4)1124 b20 = ( b(ia+i)+b3) - 0.5_wp*b111125 b21 = sin60 *((b2+b5)-(b1+b4))1126 c(ja+j) = ( a(ia+i)+a3) + a111127 d(ja+j) = ( b(ia+i)+b3) + b111128 a1 = c1 * a(ib+i) + s1 * b(ib+i) 1129 b1 = c1 * b(ib+i) - s1 * a(ib+i) 1130 a2 = c2 * a(ic+i) + s2 * b(ic+i) 1131 b2 = c2 * b(ic+i) - s2 * a(ic+i) 1132 a3 = c3 * a(id+i) + s3 * b(id+i) 1133 b3 = c3 * b(id+i) - s3 * a(id+i) 1134 a4 = c4 * a(ie+i) + s4 * b(ie+i) 1135 b4 = c4 * b(ie+i) - s4 * a(ie+i) 1136 a5 = c5 * a(if+i) + s5 * b(if+i) 1137 b5 = c5 * b(if+i) - s5 * a(if+i) 1138 a11 = ( a2 + a5 ) + ( a1 + a4 ) 1139 a20 = ( a(ia+i) + a3 ) - 0.5_wp * a11 1140 a21 = sin60 * ( ( a2 + a5 ) - ( a1 + a4 ) ) 1141 b11 = ( b2 + b5 ) + ( b1 + b4 ) 1142 b20 = ( b(ia+i) + b3 ) - 0.5_wp * b11 1143 b21 = sin60 * ( ( b2 + b5 ) - ( b1 + b4 ) ) 1144 c(ja+j) = ( a(ia+i) + a3 ) + a11 1145 d(ja+j) = ( b(ia+i) + b3 ) + b11 1128 1146 c(jc+j) = a20 - b21 1129 1147 d(jc+j) = a21 + b20 1130 1148 c(je+j) = a20 + b21 1131 1149 d(je+j) = a21 - b20 1132 a11 = ( a2-a5) + (a4-a1)1133 a20 = ( a(ia+i)-a3) - 0.5_wp*a111134 a21 = sin60 *((a4-a1)-(a2-a5))1135 b11 = ( b5-b2) - (b4-b1)1136 b20 = ( b3-b(ia+i)) - 0.5_wp*b111137 b21 = sin60 *((b5-b2)+(b4-b1))1150 a11 = ( a2 - a5 ) + ( a4 - a1 ) 1151 a20 = ( a(ia+i) - a3 ) - 0.5_wp * a11 1152 a21 = sin60 * ( ( a4 - a1 ) - ( a2 - a5 ) ) 1153 b11 = ( b5 - b2 ) - ( b4 - b1 ) 1154 b20 = ( b3 - b(ia+i) ) - 0.5_wp * b11 1155 b21 = sin60 * ( ( b5 - b2 ) + ( b4 - b1 ) ) 1138 1156 c(jb+j) = a20 - b21 1139 1157 d(jb+j) = a21 - b20 1140 c(jd+j) = a11 + ( a(ia+i)-a3)1141 d(jd+j) = b11 + ( b3-b(ia+i))1158 c(jd+j) = a11 + ( a(ia+i) - a3 ) 1159 d(jd+j) = b11 + ( b3 - b(ia+i) ) 1142 1160 c(jf+j) = a20 + b21 1143 1161 d(jf+j) = a21 + b20 … … 1169 1187 !OCL NOVREC 1170 1188 DO ijk = 1, lot 1171 c(ja+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i)) 1172 d(ja+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i)) 1173 c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i)) 1174 d(jb+j) = a(id+i) - (a(ib+i)+a(if+i)) 1175 c(jc+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i)) 1176 d(jc+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i)) 1189 c(ja+j) = ( a (ia+i) + 0.5_wp * ( a(ic+i) - a(ie+i) ) ) + sin60 * & 1190 ( a(ib+i) - a(if+i) ) 1191 d(ja+j) = - ( a(id+i) + 0.5_wp * ( a(ib+i) + a(if+i) ) ) - sin60 * ( a(ic+i) + & 1192 a(ie+i) ) 1193 c(jb+j) = a(ia+i) - ( a(ic+i) - a(ie+i) ) 1194 d(jb+j) = a(id+i) - ( a(ib+i) + a(if+i) ) 1195 c(jc+j) = ( a(ia+i) + 0.5_wp * ( a(ic+i) - a(ie+i) ) ) - sin60 * ( a(ib+i) - & 1196 a(if+i) ) 1197 d(jc+j) = - ( a(id+i) + 0.5_wp * ( a(ib+i) + a(if+i) ) ) + sin60 * ( a(ic+i) + & 1198 a(ie+i) ) 1177 1199 i = i + inc3 1178 1200 j = j + inc4 … … 1184 1206 ELSE 1185 1207 1186 z = 1.0_wp /REAL(n,KIND=wp)1187 zsin60 = z *sin601208 z = 1.0_wp / REAL( n, KIND = wp ) 1209 zsin60 = z * sin60 1188 1210 DO l = 1, la 1189 1211 i = ibase … … 1193 1215 !OCL NOVREC 1194 1216 DO ijk = 1, lot 1195 a11 = ( a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))1196 c(ja+j) = z *((a(ia+i)+a(id+i))+a11)1197 c(jc+j) = z *((a(ia+i)+a(id+i))-0.5_wp*a11)1198 d(jc+j) = zsin60 *((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))1199 a11 = ( a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))1200 c(jb+j) = z *((a(ia+i)-a(id+i))-0.5_wp*a11)1201 d(jb+j) = zsin60 *((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))1202 c(jd+j) = z *((a(ia+i)-a(id+i))+a11)1217 a11 = ( a(ic+i) + a(if+i) ) + ( a(ib+i) + a(ie+i) ) 1218 c(ja+j) = z * ( ( a(ia+i) + a(id+i) ) + a11 ) 1219 c(jc+j) = z * ( ( a(ia+i) + a(id+i) ) - 0.5_wp * a11 ) 1220 d(jc+j) = zsin60 * ( ( a(ic+i) + a(if+i) ) - ( a(ib+i) + a(ie+i) ) ) 1221 a11 = ( a(ic+i) - a(if+i) ) + ( a(ie+i) - a(ib+i) ) 1222 c(jb+j) = z * ( ( a(ia+i) - a(id+i) ) - 0.5_wp * a11 ) 1223 d(jb+j) = zsin60 * ( ( a(ie+i) - a(ib+i) ) - ( a(ic+i) - a(if+i) ) ) 1224 c(jd+j) = z * ( ( a(ia+i) - a(id+i) ) + a11 ) 1203 1225 i = i + inc3 1204 1226 j = j + inc4 … … 1229 1251 ja = 1 1230 1252 jb = ja + la * inc2 1231 jc = jb + 2 *m*inc21232 jd = jc + 2 *m*inc21233 je = jd + 2 *m*inc21234 z = 1.0_wp / REAL( n, KIND =wp )1253 jc = jb + 2 * m * inc2 1254 jd = jc + 2 * m * inc2 1255 je = jd + 2 * m * inc2 1256 z = 1.0_wp / REAL( n, KIND = wp ) 1235 1257 zsin45 = z * SQRT( 0.5_wp ) 1236 1258 … … 1243 1265 !OCL NOVREC 1244 1266 DO ijk = 1, lot 1245 c(ja+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))+((a(id+i)+ a(ih+i))+(a(ib+i)+a(if+i)))) 1246 c(je+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))-((a(id+i)+ a(ih+i))+(a(ib+i)+a(if+i)))) 1247 c(jc+j) = z*((a(ia+i)+a(ie+i))-(a(ic+i)+a(ig+i))) 1248 d(jc+j) = z*((a(id+i)+a(ih+i))-(a(ib+i)+a(if+i))) 1249 c(jb+j) = z*(a(ia+i)-a(ie+i)) + zsin45*((a(ih+i)-a(id+i))-(a(if+i)-a(ib+i))) 1250 c(jd+j) = z*(a(ia+i)-a(ie+i)) - zsin45*((a(ih+i)-a(id+i))-(a(if+i)-a(ib+i))) 1251 d(jb+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) + z*(a(ig+i)-a(ic+i)) 1252 d(jd+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) - z*(a(ig+i)-a(ic+i)) 1267 c(ja+j) = z * ( ( ( a(ia+i) + a(ie+i) ) + ( a(ic+i) + a(ig+i) ) ) + ( ( a(id+i) & 1268 + a(ih+i) ) + ( a(ib+i) + a(if+i) ) ) ) 1269 c(je+j) = z * ( ( ( a(ia+i) + a(ie+i) ) + (a(ic+i) + a(ig+i) ) ) - ( ( a(id+i) & 1270 + a(ih+i) ) + ( a(ib+i) + a(if+i) ) ) ) 1271 c(jc+j) = z * ( ( a(ia+i) + a(ie+i) ) - ( a(ic+i) + a(ig+i) ) ) 1272 d(jc+j) = z * ( ( a(id+i) + a(ih+i) ) - ( a(ib+i) + a(if+i) ) ) 1273 c(jb+j) = z * ( a(ia+i) - a(ie+i) ) + zsin45 * ( ( a(ih+i) - a(id+i) ) - ( a(if+i) & 1274 - a(ib+i) ) ) 1275 c(jd+j) = z * ( a(ia+i) - a(ie+i) ) - zsin45 * ( ( a(ih+i) - a(id+i) ) - ( a(if+i) & 1276 - a(ib+i) ) ) 1277 d(jb+j) = zsin45 * ( ( a(ih+i) - a(id+i) ) + ( a(if+i) - a(ib+i) ) ) + z * & 1278 ( a(ig+i) - a(ic+i) ) 1279 d(jd+j) = zsin45 * ( ( a(ih+i) - a(id+i) ) + ( a(if+i) - a(ib+i) ) ) - z * & 1280 ( a(ig+i) - a(ic+i) ) 1253 1281 i = i + inc3 1254 1282 j = j + inc4 … … 1262 1290 END SUBROUTINE qpassm 1263 1291 1264 !------------------------------------------------------------------------------ !1292 !--------------------------------------------------------------------------------------------------! 1265 1293 ! Description: 1266 1294 ! ------------ 1267 1295 !> Same as qpassm, but for backward fft 1268 !------------------------------------------------------------------------------ !1296 !--------------------------------------------------------------------------------------------------! 1269 1297 SUBROUTINE rpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr ) 1270 1298 … … 1273 1301 IMPLICIT NONE 1274 1302 1275 INTEGER(iwp) :: ierr !<1276 INTEGER(iwp) :: ifac !<1277 INTEGER(iwp) :: inc1 !<1278 INTEGER(iwp) :: inc2 !<1279 INTEGER(iwp) :: inc3 !<1280 INTEGER(iwp) :: inc4 !<1281 INTEGER(iwp) :: la !<1282 INTEGER(iwp) :: lot !<1283 INTEGER(iwp) :: n !<1303 INTEGER(iwp) :: ierr !< 1304 INTEGER(iwp) :: ifac !< 1305 INTEGER(iwp) :: inc1 !< 1306 INTEGER(iwp) :: inc2 !< 1307 INTEGER(iwp) :: inc3 !< 1308 INTEGER(iwp) :: inc4 !< 1309 INTEGER(iwp) :: la !< 1310 INTEGER(iwp) :: lot !< 1311 INTEGER(iwp) :: n !< 1284 1312 ! 1285 1313 !-- Arrays are dimensioned with n 1286 REAL(wp) :: a(*) !<1287 REAL(wp) :: b(*) !<1288 REAL(wp) :: c(*) !<1289 REAL(wp) :: d(*) !<1290 REAL(wp) :: trigs(*) !<1291 1292 REAL(wp) :: c1 !<1293 REAL(wp) :: c2 !<1294 REAL(wp) :: c3 !<1295 REAL(wp) :: c4 !<1296 REAL(wp) :: c5 !<1297 REAL(wp) :: qqrt5 !<1298 REAL(wp) :: qrt5 !<1299 REAL(wp) :: s1 !<1300 REAL(wp) :: s2 !<1301 REAL(wp) :: s3 !<1302 REAL(wp) :: s4 !<1303 REAL(wp) :: s5 !<1304 REAL(wp) :: sin36 !<1305 REAL(wp) :: sin45 !<1306 REAL(wp) :: sin60 !<1307 REAL(wp) :: sin72 !<1308 REAL(wp) :: ssin36 !<1309 REAL(wp) :: ssin45 !<1310 REAL(wp) :: ssin60 !<1311 REAL(wp) :: ssin72 !<1312 1313 INTEGER(iwp) :: i !<1314 INTEGER(iwp) :: ia !<1315 INTEGER(iwp) :: ib !<1316 INTEGER(iwp) :: ibase !<1317 INTEGER(iwp) :: ic !<1318 INTEGER(iwp) :: id !< 1319 INTEGER(iwp) :: ie !<1320 INTEGER(iwp) :: if !<1321 INTEGER(iwp) :: igo !<1322 INTEGER(iwp) :: iink !<1323 INTEGER(iwp) :: ijk !<1324 INTEGER(iwp) :: j !< 1325 INTEGER(iwp) :: ja !<1326 INTEGER(iwp) :: jb !<1327 INTEGER(iwp) :: jbase !<1328 INTEGER(iwp) :: jc !<1329 INTEGER(iwp) :: jd !<1330 INTEGER(iwp) :: je !<1331 INTEGER(iwp) :: jf !<1332 INTEGER(iwp) :: jg !<1333 INTEGER(iwp) :: jh !<1334 INTEGER(iwp) :: jink !<1335 INTEGER(iwp) :: jump !<1336 INTEGER(iwp) :: k !<1337 INTEGER(iwp) :: kb !<1338 INTEGER(iwp) :: kc !<1339 INTEGER(iwp) :: kd !<1340 INTEGER(iwp) :: ke !<1341 INTEGER(iwp) :: kf !<1342 INTEGER(iwp) :: kstop !<1343 INTEGER(iwp) :: l !<1344 INTEGER(iwp) :: m !<1345 1346 REAL(wp) :: a10(nfft) !<1347 REAL(wp) :: a11(nfft) !<1348 REAL(wp) :: a20(nfft) !<1349 REAL(wp) :: a21(nfft) !<1350 REAL(wp) :: b10(nfft) !<1351 REAL(wp) :: b11(nfft) !<1352 REAL(wp) :: b20(nfft) !<1353 REAL(wp) :: b21(nfft) !<1314 INTEGER(iwp) :: i !< 1315 INTEGER(iwp) :: ia !< 1316 INTEGER(iwp) :: ib !< 1317 INTEGER(iwp) :: ibase !< 1318 INTEGER(iwp) :: ic !< 1319 INTEGER(iwp) :: id !< 1320 INTEGER(iwp) :: ie !< 1321 INTEGER(iwp) :: if !< 1322 INTEGER(iwp) :: igo !< 1323 INTEGER(iwp) :: iink !< 1324 INTEGER(iwp) :: ijk !< 1325 INTEGER(iwp) :: j !< 1326 INTEGER(iwp) :: ja !< 1327 INTEGER(iwp) :: jb !< 1328 INTEGER(iwp) :: jbase !< 1329 INTEGER(iwp) :: jc !< 1330 INTEGER(iwp) :: jd !< 1331 INTEGER(iwp) :: je !< 1332 INTEGER(iwp) :: jf !< 1333 INTEGER(iwp) :: jg !< 1334 INTEGER(iwp) :: jh !< 1335 INTEGER(iwp) :: jink !< 1336 INTEGER(iwp) :: jump !< 1337 INTEGER(iwp) :: k !< 1338 INTEGER(iwp) :: kb !< 1339 INTEGER(iwp) :: kc !< 1340 INTEGER(iwp) :: kd !< 1341 INTEGER(iwp) :: ke !< 1342 INTEGER(iwp) :: kf !< 1343 INTEGER(iwp) :: kstop !< 1344 INTEGER(iwp) :: l !< 1345 INTEGER(iwp) :: m !< 1346 1347 REAL(wp) :: a(*) !< 1348 REAL(wp) :: b(*) !< 1349 REAL(wp) :: c(*) !< 1350 REAL(wp) :: d(*) !< 1351 REAL(wp) :: trigs(*) !< 1352 1353 REAL(wp) :: c1 !< 1354 REAL(wp) :: c2 !< 1355 REAL(wp) :: c3 !< 1356 REAL(wp) :: c4 !< 1357 REAL(wp) :: c5 !< 1358 REAL(wp) :: qqrt5 !< 1359 REAL(wp) :: qrt5 !< 1360 REAL(wp) :: s1 !< 1361 REAL(wp) :: s2 !< 1362 REAL(wp) :: s3 !< 1363 REAL(wp) :: s4 !< 1364 REAL(wp) :: s5 !< 1365 REAL(wp) :: sin36 !< 1366 REAL(wp) :: sin45 !< 1367 REAL(wp) :: sin60 !< 1368 REAL(wp) :: sin72 !< 1369 REAL(wp) :: ssin36 !< 1370 REAL(wp) :: ssin45 !< 1371 REAL(wp) :: ssin60 !< 1372 REAL(wp) :: ssin72 !< 1373 1374 REAL(wp) :: a10(nfft) !< 1375 REAL(wp) :: a11(nfft) !< 1376 REAL(wp) :: a20(nfft) !< 1377 REAL(wp) :: a21(nfft) !< 1378 REAL(wp) :: b10(nfft) !< 1379 REAL(wp) :: b11(nfft) !< 1380 REAL(wp) :: b20(nfft) !< 1381 REAL(wp) :: b21(nfft) !< 1354 1382 1355 1383 … … 1363 1391 iink = la * inc1 1364 1392 jink = la * inc2 1365 jump = ( ifac-1) * jink1366 kstop = ( n-ifac) / (2*ifac)1393 jump = ( ifac - 1 ) * jink 1394 kstop = ( n - ifac ) / ( 2 * ifac ) 1367 1395 1368 1396 IF ( lot > nfft ) THEN … … 1386 1414 1387 1415 ia = 1 1388 ib = ia + ( 2*m-la) * inc11416 ib = ia + ( 2 * m - la ) * inc1 1389 1417 ja = 1 1390 1418 jb = ja + jink … … 1409 1437 1410 1438 ia = ia + iink 1411 iink = 2 * iink1439 iink = 2 * iink 1412 1440 ib = ib - iink 1413 1441 ibase = 0 … … 1431 1459 c(ja+j) = a(ia+i) + a(ib+i) 1432 1460 d(ja+j) = b(ia+i) - b(ib+i) 1433 c(jb+j) = c1 *(a(ia+i)-a(ib+i)) - s1*(b(ia+i)+b(ib+i))1434 d(jb+j) = s1 *(a(ia+i)-a(ib+i)) + c1*(b(ia+i)+b(ib+i))1461 c(jb+j) = c1 * ( a(ia+i) - a(ib+i) ) - s1 * ( b(ia+i) + b(ib+i) ) 1462 d(jb+j) = s1 * ( a(ia+i) - a(ib+i) ) + c1 * ( b(ia+i) + b(ib+i) ) 1435 1463 i = i + inc3 1436 1464 j = j + inc4 … … 1457 1485 DO ijk = 1, lot 1458 1486 c(ja+j) = a(ia+i) 1459 c(jb+j) = - b(ia+i)1487 c(jb+j) = - b(ia+i) 1460 1488 i = i + inc3 1461 1489 j = j + inc4 … … 1474 1502 !OCL NOVREC 1475 1503 DO ijk = 1, lot 1476 c(ja+j) = 2.0_wp *(a(ia+i)+a(ib+i))1477 c(jb+j) = 2.0_wp *(a(ia+i)-a(ib+i))1504 c(ja+j) = 2.0_wp * ( a(ia+i) + a(ib+i) ) 1505 c(jb+j) = 2.0_wp * ( a(ia+i) - a(ib+i) ) 1478 1506 i = i + inc3 1479 1507 j = j + inc4 … … 1490 1518 1491 1519 ia = 1 1492 ib = ia + ( 2*m-la) * inc11520 ib = ia + ( 2 * m - la ) * inc1 1493 1521 ic = ib 1494 1522 ja = 1 … … 1506 1534 DO ijk = 1, lot 1507 1535 c(ja+j) = a(ia+i) + a(ib+i) 1508 c(jb+j) = ( a(ia+i)-0.5_wp*a(ib+i)) - (sin60*(b(ib+i)))1509 c(jc+j) = ( a(ia+i)-0.5_wp*a(ib+i)) + (sin60*(b(ib+i)))1536 c(jb+j) = ( a(ia+i) - 0.5_wp * a(ib+i) ) - ( sin60 * ( b(ib+i) ) ) 1537 c(jc+j) = ( a(ia+i) - 0.5_wp * a(ib+i) ) + ( sin60 * ( b(ib+i) ) ) 1510 1538 i = i + inc3 1511 1539 j = j + inc4 … … 1538 1566 !OCL NOVREC 1539 1567 DO ijk = 1, lot 1540 c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i)) 1541 d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i)) 1542 c(jb+j) = c1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ b(ic+i)))) & 1543 - s1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- a(ic+i)))) 1544 d(jb+j) = s1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ b(ic+i)))) & 1545 + c1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- a(ic+i)))) 1546 c(jc+j) = c2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ b(ic+i)))) & 1547 - s2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- a(ic+i)))) 1548 d(jc+j) = s2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ b(ic+i)))) & 1549 + c2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- a(ic+i)))) 1568 c(ja+j) = a(ia+i) + ( a(ib+i) + a(ic+i) ) 1569 d(ja+j) = b(ia+i) + ( b(ib+i) - b(ic+i) ) 1570 c(jb+j) = c1 * ( ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) - ( sin60 & 1571 * ( b(ib+i) + b(ic+i) ) ) ) - s1 * ( ( b(ia+i) - 0.5_wp & 1572 * ( b(ib+i) - b(ic+i) ) ) + ( sin60 * ( a(ib+i) - a(ic+i) ) ) ) 1573 d(jb+j) = s1 * ( ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) - ( sin60 & 1574 * ( b(ib+i) + b(ic+i) ) ) ) + c1 * ( ( b(ia+i) - 0.5_wp & 1575 * ( b(ib+i) - b(ic+i) ) ) + ( sin60 * ( a(ib+i) - a(ic+i) ) ) ) 1576 c(jc+j) = c2 * ( ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) + ( sin60 & 1577 * ( b(ib+i) + b(ic+i) ) ) ) - s2 * ( ( b(ia+i) - 0.5_wp & 1578 * ( b(ib+i) - b(ic+i) ) ) - ( sin60 * ( a(ib+i) - a(ic+i) ) ) ) 1579 d(jc+j) = s2 * ( ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) + ( sin60 & 1580 * ( b(ib+i) + b(ic+i) ) ) ) + c2 * ( ( b(ia+i) - 0.5_wp & 1581 * ( b(ib+i) - b(ic+i) ) ) - ( sin60 * ( a(ib+i) - a(ic+i) ) ) ) 1550 1582 i = i + inc3 1551 1583 j = j + inc4 … … 1573 1605 DO ijk = 1, lot 1574 1606 c(ja+j) = a(ia+i) + a(ib+i) 1575 c(jb+j) = ( 0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))1576 c(jc+j) = - (0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))1607 c(jb+j) = ( 0.5_wp * a(ia+i) - a(ib+i) ) - ( sin60 * b(ia+i) ) 1608 c(jc+j) = - ( 0.5_wp * a(ia+i) - a(ib+i) ) - ( sin60 * b(ia+i) ) 1577 1609 i = i + inc3 1578 1610 j = j + inc4 … … 1592 1624 !OCL NOVREC 1593 1625 DO ijk = 1, lot 1594 c(ja+j) = 2.0_wp *(a(ia+i)+a(ib+i))1595 c(jb+j) = ( 2.0_wp*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i))1596 c(jc+j) = ( 2.0_wp*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i))1626 c(ja+j) = 2.0_wp * ( a(ia+i) + a(ib+i) ) 1627 c(jb+j) = ( 2.0_wp * a(ia+i) - a(ib+i) ) - ( ssin60 * b(ib+i) ) 1628 c(jc+j) = ( 2.0_wp * a(ia+i) - a(ib+i) ) + ( ssin60 * b(ib+i) ) 1597 1629 i = i + inc3 1598 1630 j = j + inc4 … … 1609 1641 1610 1642 ia = 1 1611 ib = ia + ( 2*m-la) * inc11612 ic = ib + 2 *m*inc11643 ib = ia + ( 2 * m - la ) * inc1 1644 ic = ib + 2 * m * inc1 1613 1645 id = ib 1614 1646 ja = 1 … … 1626 1658 !OCL NOVREC 1627 1659 DO ijk = 1, lot 1628 c(ja+j) = ( a(ia+i)+a(ic+i)) + a(ib+i)1629 c(jb+j) = ( a(ia+i)-a(ic+i)) - b(ib+i)1630 c(jc+j) = ( a(ia+i)+a(ic+i)) - a(ib+i)1631 c(jd+j) = ( a(ia+i)-a(ic+i)) + b(ib+i)1660 c(ja+j) = ( a(ia+i) + a(ic+i) ) + a(ib+i) 1661 c(jb+j) = ( a(ia+i) - a(ic+i) ) - b(ib+i) 1662 c(jc+j) = ( a(ia+i) + a(ic+i) ) - a(ib+i) 1663 c(jd+j) = ( a(ia+i) - a(ic+i) ) + b(ib+i) 1632 1664 i = i + inc3 1633 1665 j = j + inc4 … … 1636 1668 jbase = jbase + inc2 1637 1669 ENDDO 1638 ia = ia + iink1639 iink = 2 * iink1640 ib = ib + iink1641 ic = ic - iink1642 id = id - iink1670 ia = ia + iink 1671 iink = 2 * iink 1672 ib = ib + iink 1673 ic = ic - iink 1674 id = id - iink 1643 1675 jbase = jbase + jump 1644 jump = 2 * jump + jink1676 jump = 2 * jump + jink 1645 1677 1646 1678 IF ( ib /= ic ) THEN … … 1664 1696 !OCL NOVREC 1665 1697 DO ijk = 1, lot 1666 c(ja+j) = ( a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))1667 d(ja+j) = ( b(ia+i)-b(ic+i)) + (b(ib+i)-b(id+i))1668 c(jc+j) = c2 *((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) - s2*((b(ia+i)-b(ic+i))&1669 -(b(ib+i)-b(id+i)))1670 d(jc+j) = s2 *((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) + c2*((b(ia+i)-b(ic+i))&1671 -(b(ib+i)-b(id+i)))1672 c(jb+j) = c1 *((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) - s1*((b(ia+i)+b(ic+i))&1673 +(a(ib+i)-a(id+i)))1674 d(jb+j) = s1 *((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) + c1*((b(ia+i)+b(ic+i))&1675 +(a(ib+i)-a(id+i)))1676 c(jd+j) = c3 *((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) - s3*((b(ia+i)+b(ic+i))&1677 -(a(ib+i)-a(id+i)))1678 d(jd+j) = s3 *((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) + c3*((b(ia+i)+b(ic+i))&1679 -(a(ib+i)-a(id+i)))1698 c(ja+j) = ( a(ia+i) + a(ic+i) ) + ( a(ib+i) + a(id+i) ) 1699 d(ja+j) = ( b(ia+i) - b(ic+i) ) + ( b(ib+i) - b(id+i) ) 1700 c(jc+j) = c2 * ( ( a(ia+i) + a(ic+i) ) - ( a(ib+i) + a(id+i) ) ) - s2 & 1701 * ( ( b(ia+i) - b(ic+i) ) - ( b(ib+i) - b(id+i) ) ) 1702 d(jc+j) = s2 * ( ( a(ia+i) + a(ic+i) ) - ( a(ib+i) + a(id+i) ) ) + c2 & 1703 * ( ( b(ia+i) - b(ic+i) ) - ( b(ib+i) - b(id+i) ) ) 1704 c(jb+j) = c1 * ( ( a(ia+i) - a(ic+i) ) - ( b(ib+i) + b(id+i) ) ) - s1 & 1705 * ( ( b(ia+i) + b(ic+i) ) + ( a(ib+i) - a(id+i) ) ) 1706 d(jb+j) = s1 * ( ( a(ia+i) - a(ic+i) ) - ( b (ib+i) + b(id+i) ) ) + c1 & 1707 * ( ( b(ia+i) + b(ic+i) ) + ( a(ib+i) - a(id+i) ) ) 1708 c(jd+j) = c3 * ( ( a(ia+i) - a(ic+i) ) + ( b(ib+i) + b(id+i) ) ) - s3 & 1709 * ( ( b(ia+i) + b(ic+i) ) - ( a(ib+i) - a(id+i) ) ) 1710 d(jd+j) = s3 * ( ( a(ia+i) - a(ic+i) ) + ( b(ib+i) + b(id+i) ) ) + c3 & 1711 * ( ( b(ia+i) + b(ic+i) ) - ( a(ib+i) - a(id+i) ) ) 1680 1712 i = i + inc3 1681 1713 j = j + inc4 … … 1705 1737 DO ijk = 1, lot 1706 1738 c(ja+j) = a(ia+i) + a(ib+i) 1707 c(jb+j) = sin45 *((a(ia+i)-a(ib+i))-(b(ia+i)+b(ib+i)))1739 c(jb+j) = sin45 * ( ( a(ia+i) - a(ib+i) ) - ( b(ia+i) + b(ib+i) ) ) 1708 1740 c(jc+j) = b(ib+i) - b(ia+i) 1709 c(jd+j) = - sin45*((a(ia+i)-a(ib+i))+(b(ia+i)+b(ib+i)))1741 c(jd+j) = - sin45 * ( ( a(ia+i) - a(ib+i) ) + ( b(ia+i) + b(ib+i) ) ) 1710 1742 i = i + inc3 1711 1743 j = j + inc4 … … 1724 1756 !OCL NOVREC 1725 1757 DO ijk = 1, lot 1726 c(ja+j) = 2.0_wp *((a(ia+i)+a(ic+i))+a(ib+i))1727 c(jb+j) = 2.0_wp *((a(ia+i)-a(ic+i))-b(ib+i))1728 c(jc+j) = 2.0_wp *((a(ia+i)+a(ic+i))-a(ib+i))1729 c(jd+j) = 2.0_wp *((a(ia+i)-a(ic+i))+b(ib+i))1758 c(ja+j) = 2.0_wp * ( ( a(ia+i) + a(ic+i) ) + a(ib+i) ) 1759 c(jb+j) = 2.0_wp * ( ( a(ia+i) - a(ic+i) ) - b(ib+i) ) 1760 c(jc+j) = 2.0_wp * ( ( a(ia+i) + a(ic+i) ) - a(ib+i) ) 1761 c(jd+j) = 2.0_wp * ( ( a(ia+i) - a(ic+i) ) + b(ib+i) ) 1730 1762 i = i + inc3 1731 1763 j = j + inc4 … … 1742 1774 1743 1775 ia = 1 1744 ib = ia + ( 2*m-la) * inc11745 ic = ib + 2 *m*inc11776 ib = ia + ( 2 * m - la ) * inc1 1777 ic = ib + 2 * m * inc1 1746 1778 id = ic 1747 1779 ie = ib … … 1761 1793 !OCL NOVREC 1762 1794 DO ijk = 1, lot 1763 c(ja+j) = a(ia+i) + ( a(ib+i)+a(ic+i))1764 c(jb+j) = ( (a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) -&1765 ( sin72*b(ib+i)+sin36*b(ic+i))1766 c(jc+j) = ( (a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) -&1767 ( sin36*b(ib+i)-sin72*b(ic+i))1768 c(jd+j) = ( (a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) +&1769 ( sin36*b(ib+i)-sin72*b(ic+i))1770 c(je+j) = ( (a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) +&1771 ( sin72*b(ib+i)+sin36*b(ic+i))1795 c(ja+j) = a(ia+i) + ( a(ib+i) + a(ic+i) ) 1796 c(jb+j) = ( ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) + qrt5 * & 1797 ( a(ib+i) - a(ic+i) ) ) - ( sin72 * b(ib+i) + sin36 * b(ic+i) ) 1798 c(jc+j) = ( ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) - qrt5 * & 1799 ( a(ib+i) - a(ic+i) ) ) - ( sin36 * b(ib+i) - sin72 * b(ic+i) ) 1800 c(jd+j) = ( ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) - qrt5 * & 1801 ( a(ib+i) - a(ic+i) ) ) + ( sin36 * b(ib+i) - sin72 * b(ic+i) ) 1802 c(je+j) = ( ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) + qrt5 * & 1803 ( a(ib+i) - a(ic+i) ) ) + ( sin72 * b(ib+i) + sin36 * b(ic+i) ) 1772 1804 i = i + inc3 1773 1805 j = j + inc4 … … 1776 1808 jbase = jbase + inc2 1777 1809 ENDDO 1778 ia = ia + iink1779 iink = 2 * iink1780 ib = ib + iink1781 ic = ic + iink1782 id = id - iink1783 ie = ie - iink1810 ia = ia + iink 1811 iink = 2 * iink 1812 ib = ib + iink 1813 ic = ic + iink 1814 id = id - iink 1815 ie = ie - iink 1784 1816 jbase = jbase + jump 1785 jump = 2 * jump + jink1817 jump = 2 * jump + jink 1786 1818 1787 1819 IF ( ib /= id ) THEN … … 1808 1840 !OCL NOVREC 1809 1841 DO ijk = 1, lot 1810 a10(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + & 1811 qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i))) 1812 a20(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - & 1813 qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i))) 1814 b10(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + & 1815 qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i))) 1816 b20(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - & 1817 qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i))) 1818 a11(ijk) = sin72*(b(ib+i)+b(ie+i)) + sin36*(b(ic+i)+b(id+i)) 1819 a21(ijk) = sin36*(b(ib+i)+b(ie+i)) - sin72*(b(ic+i)+b(id+i)) 1820 b11(ijk) = sin72*(a(ib+i)-a(ie+i)) + sin36*(a(ic+i)-a(id+i)) 1821 b21(ijk) = sin36*(a(ib+i)-a(ie+i)) - sin72*(a(ic+i)-a(id+i)) 1822 c(ja+j) = a(ia+i) + ((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i))) 1823 d(ja+j) = b(ia+i) + ((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i))) 1824 c(jb+j) = c1*(a10(ijk)-a11(ijk)) - s1*(b10(ijk)+b11(ijk)) 1825 d(jb+j) = s1*(a10(ijk)-a11(ijk)) + c1*(b10(ijk)+b11(ijk)) 1826 c(je+j) = c4*(a10(ijk)+a11(ijk)) - s4*(b10(ijk)-b11(ijk)) 1827 d(je+j) = s4*(a10(ijk)+a11(ijk)) + c4*(b10(ijk)-b11(ijk)) 1828 c(jc+j) = c2*(a20(ijk)-a21(ijk)) - s2*(b20(ijk)+b21(ijk)) 1829 d(jc+j) = s2*(a20(ijk)-a21(ijk)) + c2*(b20(ijk)+b21(ijk)) 1830 c(jd+j) = c3*(a20(ijk)+a21(ijk)) - s3*(b20(ijk)-b21(ijk)) 1831 d(jd+j) = s3*(a20(ijk)+a21(ijk)) + c3*(b20(ijk)-b21(ijk)) 1832 1842 a10(ijk) = ( a(ia+i) - 0.25_wp * ( ( a(ib+i) + a(ie+i) ) + ( a(ic+i) & 1843 + a(id+i) ) ) ) + qrt5 * ( ( a(ib+i) + a(ie+i) ) - ( a(ic+i) & 1844 + a(id+i) ) ) 1845 a20(ijk) = ( a(ia+i) - 0.25_wp * ( ( a(ib+i) + a(ie+i) ) + ( a(ic+i) & 1846 + a(id+i) ) ) ) - qrt5 * ( ( a(ib+i) + a(ie+i) ) - ( a(ic+i) & 1847 + a(id+i) ) ) 1848 b10(ijk) = ( b(ia+i) - 0.25_wp * ( ( b(ib+i) - b(ie+i) ) + ( b(ic+i) & 1849 - b(id+i) ) ) ) + qrt5 * ( ( b(ib+i) - b(ie+i) ) - ( b(ic+i) & 1850 - b(id+i) ) ) 1851 b20(ijk) = ( b(ia+i) - 0.25_wp * ( ( b(ib+i) - b(ie+i) ) + ( b(ic+i) & 1852 - b(id+i) ) ) ) - qrt5 * ( ( b(ib+i) - b(ie+i) ) - ( b(ic+i) & 1853 - b(id+i) ) ) 1854 a11(ijk) = sin72 * ( b(ib+i) + b(ie+i) ) + sin36 * ( b(ic+i) + b(id+i) ) 1855 a21(ijk) = sin36 * ( b(ib+i) + b(ie+i) ) - sin72 * ( b(ic+i) + b(id+i) ) 1856 b11(ijk) = sin72 * ( a(ib+i) - a(ie+i) ) + sin36 * ( a(ic+i) - a(id+i) ) 1857 b21(ijk) = sin36 * ( a(ib+i) - a(ie+i) ) - sin72 * ( a(ic+i) - a(id+i) ) 1858 c(ja+j) = a(ia+i) + ( ( a(ib+i) + a(ie+i) ) + ( a(ic+i) + a(id+i) ) ) 1859 d(ja+j) = b(ia+i) + ( ( b(ib+i) - b(ie+i) ) + ( b(ic+i) - b(id+i) ) ) 1860 c(jb+j) = c1 * ( a10(ijk) - a11(ijk) ) - s1 * ( b10(ijk) + b11(ijk) ) 1861 d(jb+j) = s1 * ( a10(ijk) - a11(ijk) ) + c1 * ( b10(ijk) + b11(ijk) ) 1862 c(je+j) = c4 * ( a10(ijk) + a11(ijk) ) - s4 * ( b10(ijk) - b11(ijk) ) 1863 d(je+j) = s4 * ( a10(ijk) + a11(ijk) ) + c4 * ( b10(ijk) - b11(ijk) ) 1864 c(jc+j) = c2 * ( a20(ijk) - a21(ijk) ) - s2 * ( b20(ijk) + b21(ijk) ) 1865 d(jc+j) = s2 * ( a20(ijk) - a21(ijk) ) + c2 * ( b20(ijk) + b21(ijk) ) 1866 c(jd+j) = c3 * ( a20(ijk) + a21(ijk) ) - s3 * ( b20(ijk) - b21(ijk) ) 1867 d(jd+j) = s3 * ( a20(ijk) + a21(ijk) ) + c3 * ( b20(ijk) - b21(ijk) ) 1833 1868 i = i + inc3 1834 1869 j = j + inc4 … … 1857 1892 !OCL NOVREC 1858 1893 DO ijk = 1, lot 1859 c(ja+j) = ( a(ia+i)+a(ib+i)) + a(ic+i)1860 c(jb+j) = ( qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -&1861 (sin36*b(ia+i)+sin72*b(ib+i))1862 c(je+j) = - (qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -&1863 (sin36*b(ia+i)+sin72*b(ib+i))1864 c(jc+j) = ( qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -&1865 (sin72*b(ia+i)-sin36*b(ib+i))1866 c(jd+j) = - (qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -&1867 (sin72*b(ia+i)-sin36*b(ib+i))1894 c(ja+j) = ( a(ia+i) + a(ib+i) ) + a(ic+i) 1895 c(jb+j) = ( qrt5 * ( a(ia+i) - a(ib+i) ) + ( 0.25_wp * ( a(ia+i) + a(ib+i) ) & 1896 - a(ic+i) ) ) - ( sin36 * b(ia+i) + sin72 * b(ib+i) ) 1897 c(je+j) = - ( qrt5 * ( a(ia+i) - a(ib+i) ) + ( 0.25_wp * ( a(ia+i) + a(ib+i) ) & 1898 - a(ic+i) ) ) - ( sin36 * b(ia+i) + sin72 * b(ib+i) ) 1899 c(jc+j) = ( qrt5 * ( a(ia+i) - a(ib+i) ) - ( 0.25_wp * ( a(ia+i) + a(ib+i) ) & 1900 - a(ic+i) ) ) - ( sin72 * b(ia+i) - sin36 * b(ib+i) ) 1901 c(jd+j) = - ( qrt5 * ( a(ia+i) - a(ib+i) ) - ( 0.25_wp * ( a(ia+i) + a(ib+i) ) & 1902 - a(ic+i) ) ) - ( sin72 * b(ia+i) - sin36 * b(ib+i) ) 1868 1903 i = i + inc3 1869 1904 j = j + inc4 … … 1885 1920 !OCL NOVREC 1886 1921 DO ijk = 1, lot 1887 c(ja+j) = 2.0_wp *(a(ia+i)+(a(ib+i)+a(ic+i)))1888 c(jb+j) = ( 2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+i)))&1889 - (ssin72*b(ib+i)+ssin36*b(ic+i))1890 c(jc+j) = ( 2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+i)))&1891 - (ssin36*b(ib+i)-ssin72*b(ic+i))1892 c(jd+j) = ( 2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+i)))&1893 + (ssin36*b(ib+i)-ssin72*b(ic+i))1894 c(je+j) = ( 2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+i)))&1895 + (ssin72*b(ib+i)+ssin36*b(ic+i))1922 c(ja+j) = 2.0_wp * ( a(ia+i) + ( a(ib+i) + a(ic+i) ) ) 1923 c(jb+j) = ( 2.0_wp * ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) + qqrt5 & 1924 * ( a(ib+i) - a(ic+i) ) ) - ( ssin72 * b(ib+i) + ssin36 * b(ic+i) ) 1925 c(jc+j) = ( 2.0_wp * ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) - qqrt5 & 1926 * ( a(ib+i) - a(ic+i) ) ) - ( ssin36 * b(ib+i) - ssin72 * b(ic+i) ) 1927 c(jd+j) = ( 2.0_wp * ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) - qqrt5 & 1928 * ( a(ib+i) - a(ic+i) ) ) + ( ssin36 * b(ib+i) - ssin72 * b(ic+i) ) 1929 c(je+j) = ( 2.0_wp * ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) + qqrt5 & 1930 * ( a(ib+i) - a(ic+i) ) ) + (ssin72 * b(ib+i) + ssin36 * b(ic+i) ) 1896 1931 i = i + inc3 1897 1932 j = j + inc4 … … 1908 1943 1909 1944 ia = 1 1910 ib = ia + ( 2*m-la) * inc11911 ic = ib + 2 *m*inc11912 id = ic + 2 *m*inc11945 ib = ia + ( 2 * m - la ) * inc1 1946 ic = ib + 2 * m * inc1 1947 id = ic + 2 * m * inc1 1913 1948 ie = ic 1914 1949 if = ib … … 1929 1964 !OCL NOVREC 1930 1965 DO ijk = 1, lot 1931 c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i)) 1932 c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i)) 1933 c(jb+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) - (sin60*(b(ib+i)+b(ic+i))) 1934 c(jf+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) + (sin60*(b(ib+i)+b(ic+i))) 1935 c(jc+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) - (sin60*(b(ib+i)-b(ic+i))) 1936 c(je+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) + (sin60*(b(ib+i)-b(ic+i))) 1966 c(ja+j) = ( a (ia+i) + a(id+i) ) + ( a(ib+i) + a(ic+i) ) 1967 c(jd+j) = ( a (ia+i) - a(id+i) ) - ( a(ib+i) - a(ic+i) ) 1968 c(jb+j) = ( ( a(ia+i) - a(id+i) ) + 0.5_wp * ( a(ib+i) - a(ic+i) ) ) & 1969 - ( sin60 * ( b(ib+i) + b(ic+i) ) ) 1970 c(jf+j) = ( ( a(ia+i) - a(id+i) ) + 0.5_wp * ( a(ib+i) - a(ic+i) ) ) & 1971 + ( sin60 * b(ib+i) + b(ic+i) ) 1972 c(jc+j) = ( ( a(ia+i) + a(id+i) ) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) & 1973 - ( sin60 * b(ib+i) - b(ic+i) ) 1974 c(je+j) = ( ( a(ia+i) + a(id+i) ) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) & 1975 + ( sin60 * b(ib+i) - b(ic+i) ) 1937 1976 i = i + inc3 1938 1977 j = j + inc4 … … 1941 1980 jbase = jbase + inc2 1942 1981 ENDDO 1943 ia = ia + iink1944 iink = 2 * iink1945 ib = ib + iink1946 ic = ic + iink1947 id = id - iink1948 ie = ie - iink1949 if = if - iink1982 ia = ia + iink 1983 iink = 2 * iink 1984 ib = ib + iink 1985 ic = ic + iink 1986 id = id - iink 1987 ie = ie - iink 1988 if = if - iink 1950 1989 jbase = jbase + jump 1951 jump = 2 * jump + jink1990 jump = 2 * jump + jink 1952 1991 1953 1992 IF ( ic /= id ) THEN … … 1977 2016 !OCL NOVREC 1978 2017 DO ijk = 1, lot 1979 a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i)) 1980 a20(ijk) = (a(ia+i)+a(id+i)) - 0.5_wp*a11(ijk) 1981 a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i))) 1982 b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i)) 1983 b20(ijk) = (b(ia+i)-b(id+i)) - 0.5_wp*b11(ijk) 1984 b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i))) 1985 1986 c(ja+j) = (a(ia+i)+a(id+i)) + a11(ijk) 1987 d(ja+j) = (b(ia+i)-b(id+i)) + b11(ijk) 1988 c(jc+j) = c2*(a20(ijk)-b21(ijk)) - s2*(b20(ijk)+a21(ijk)) 1989 d(jc+j) = s2*(a20(ijk)-b21(ijk)) + c2*(b20(ijk)+a21(ijk)) 1990 c(je+j) = c4*(a20(ijk)+b21(ijk)) - s4*(b20(ijk)-a21(ijk)) 1991 d(je+j) = s4*(a20(ijk)+b21(ijk)) + c4*(b20(ijk)-a21(ijk)) 1992 1993 a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i)) 1994 b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i)) 1995 a20(ijk) = (a(ia+i)-a(id+i)) - 0.5_wp*a11(ijk) 1996 a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i))) 1997 b20(ijk) = (b(ia+i)+b(id+i)) + 0.5_wp*b11(ijk) 1998 b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i))) 1999 2000 c(jd+j) = c3*((a(ia+i)-a(id+i))+a11(ijk)) - s3*((b(ia+i)+b(id+i))-b11(ijk)) 2001 d(jd+j) = s3*((a(ia+i)-a(id+i))+a11(ijk)) + c3*((b(ia+i)+b(id+i))-b11(ijk)) 2002 c(jb+j) = c1*(a20(ijk)-b21(ijk)) - s1*(b20(ijk)-a21(ijk)) 2003 d(jb+j) = s1*(a20(ijk)-b21(ijk)) + c1*(b20(ijk)-a21(ijk)) 2004 c(jf+j) = c5*(a20(ijk)+b21(ijk)) - s5*(b20(ijk)+a21(ijk)) 2005 d(jf+j) = s5*(a20(ijk)+b21(ijk)) + c5*(b20(ijk)+a21(ijk)) 2018 a11(ijk) = ( a(ie+i) + a(ib+i) ) + ( a(ic+i) + a(if+i) ) 2019 a20(ijk) = ( a(ia+i) + a(id+i) ) - 0.5_wp * a11(ijk) 2020 a21(ijk) = sin60 * ( ( a(ie+i) + a(ib+i) ) - ( a(ic+i) + a(if+i) ) ) 2021 b11(ijk) = ( b(ib+i) - b(ie+i) ) + ( b(ic+i) - b(if+i) ) 2022 b20(ijk) = ( b(ia+i) - b(id+i) ) - 0.5_wp * b11(ijk) 2023 b21(ijk) = sin60 * ( ( b(ib+i) - b(ie+i) ) - ( b(ic+i) - b(if+i) ) ) 2024 2025 c(ja+j) = ( a(ia+i) + a(id+i) ) + a11(ijk) 2026 d(ja+j) = ( b(ia+i) - b(id+i) ) + b11(ijk) 2027 c(jc+j) = c2 * ( a20(ijk) - b21(ijk) ) - s2 * ( b20(ijk) + a21(ijk) ) 2028 d(jc+j) = s2 * ( a20(ijk) - b21(ijk) ) + c2 * ( b20(ijk) + a21(ijk) ) 2029 c(je+j) = c4 * ( a20(ijk) + b21(ijk) ) - s4 * ( b20(ijk) - a21(ijk) ) 2030 d(je+j) = s4 * ( a20(ijk) + b21(ijk) ) + c4 * ( b20(ijk) - a21(ijk) ) 2031 2032 a11(ijk) = ( a(ie+i) - a(ib+i) ) + ( a(ic+i) - a(if+i) ) 2033 b11(ijk) = ( b(ie+i) + b(ib+i) ) - ( b(ic+i) + b(if+i) ) 2034 a20(ijk) = ( a(ia+i) - a(id+i) ) - 0.5_wp * a11(ijk) 2035 a21(ijk) = sin60 * ( ( a(ie+i) - a(ib+i) ) - ( a(ic+i) - a(if+i) ) ) 2036 b20(ijk) = ( b(ia+i) + b(id+i) ) + 0.5_wp * b11(ijk) 2037 b21(ijk) = sin60 * ( ( b(ie+i) + b(ib+i) ) + ( b(ic+i) + b(if+i) ) ) 2038 2039 c(jd+j) = c3 * ( (a(ia+i) - a(id+i) ) + a11(ijk) ) - s3 & 2040 * ( (b(ia+i) + b(id+i) ) - b11(ijk) ) 2041 d(jd+j) = s3 * ( (a(ia+i) - a(id+i) ) + a11(ijk) ) + c3 & 2042 * ( (b(ia+i) + b(id+i) ) - b11(ijk) ) 2043 c(jb+j) = c1 * ( a20(ijk) - b21(ijk) ) - s1 * ( b20(ijk) - a21(ijk) ) 2044 d(jb+j) = s1 * ( a20(ijk) - b21(ijk) ) + c1 * ( b20(ijk) - a21(ijk) ) 2045 c(jf+j) = c5 * ( a20(ijk) + b21(ijk) ) - s5 * ( b20(ijk) + a21(ijk) ) 2046 d(jf+j) = s5 * ( a20(ijk) + b21(ijk) ) + c5 * ( b20(ijk) + a21(ijk) ) 2006 2047 2007 2048 i = i + inc3 … … 2032 2073 !OCL NOVREC 2033 2074 DO ijk = 1, lot 2034 c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i)) 2035 c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i)) 2036 c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i)) 2037 c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i)) 2038 c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i)) 2039 c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i)) 2075 c(ja+j) = a(ib+i) + ( a(ia+i) + a(ic+i) ) 2076 c(jd+j) = b(ib+i) - ( b(ia+i) + b(ic+i) ) 2077 c(jb+j) = ( sin60 * ( a(ia+i) - a(ic+i) ) ) - ( 0.5_wp * ( b(ia+i) + b(ic+i) ) & 2078 + b(ib+i) ) 2079 c(jf+j) = - ( sin60 * ( a(ia+i) - a(ic+i) ) ) - ( 0.5_wp * ( b(ia+i) & 2080 + b(ic+i) ) + b(ib+i) ) 2081 c(jc+j) = sin60 * ( b(ic+i) - b(ia+i) ) + ( 0.5_wp * ( a(ia+i) + a(ic+i) ) & 2082 - a(ib+i) ) 2083 c(je+j) = sin60 * ( b(ic+i) - b(ia+i) ) - ( 0.5_wp * ( a(ia+i) + a(ic+i) ) & 2084 - a(ib+i) ) 2040 2085 i = i + inc3 2041 2086 j = j + inc4 … … 2055 2100 !OCL NOVREC 2056 2101 DO ijk = 1, lot 2057 c(ja+j) = (2.0_wp*(a(ia+i)+a(id+i))) + (2.0_wp*(a(ib+i)+a(ic+i))) 2058 c(jd+j) = (2.0_wp*(a(ia+i)-a(id+i))) - (2.0_wp*(a(ib+i)-a(ic+i))) 2059 c(jb+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+i)+b(ic+i))) 2060 c(jf+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+i)+b(ic+i))) 2061 c(jc+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+i)-b(ic+i))) 2062 c(je+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+i)-b(ic+i))) 2102 c(ja+j) = ( 2.0_wp * ( a(ia+i) + a(id+i) ) ) + ( 2.0_wp * ( a(ib+i) + a(ic+i) ) ) 2103 c(jd+j) = ( 2.0_wp * ( a(ia+i) - a(id+i) ) ) - ( 2.0_wp * ( a(ib+i) - a(ic+i) ) ) 2104 c(jb+j) = ( 2.0_wp * ( a(ia+i) - a(id+i) ) + ( a(ib+i) - a(ic+i) ) ) & 2105 - ( ssin60 * ( b(ib+i) + b(ic+i) ) ) 2106 c(jf+j) = ( 2.0_wp * ( a(ia+i) - a(id+i) ) + ( a(ib+i) - a(ic+i) ) ) & 2107 + ( ssin60 * ( b(ib+i) + b(ic+i) ) ) 2108 c(jc+j) = ( 2.0_wp * ( a(ia+i) + a(id+i) ) - ( a(ib+i) + a(ic+i) ) ) & 2109 - ( ssin60 * ( b(ib+i) - b(ic+i) ) ) 2110 c(je+j) = ( 2.0_wp * ( a(ia+i) + a(id+i) ) - ( a(ib+i) + a(ic+i) ) ) & 2111 + ( ssin60 * ( b(ib+i) - b(ic+i) ) ) 2063 2112 i = i + inc3 2064 2113 j = j + inc4 … … 2080 2129 2081 2130 ia = 1 2082 ib = ia + la *inc12083 ic = ib + 2 *la*inc12084 id = ic + 2 *la*inc12085 ie = id + 2 *la*inc12131 ib = ia + la * inc1 2132 ic = ib + 2 * la * inc1 2133 id = ic + 2 * la * inc1 2134 ie = id + 2 * la * inc1 2086 2135 ja = 1 2087 2136 jb = ja + jink … … 2101 2150 !OCL NOVREC 2102 2151 DO ijk = 1, lot 2103 c(ja+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i))) 2104 c(je+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i))) 2105 c(jc+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i))) 2106 c(jg+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i))) 2107 c(jb+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45*((a(ib+i)-a(id+i))-(b(ib+i)+b(id+i))) 2108 c(jf+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+i))-(b(ib+i)+b(id+i))) 2109 c(jd+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+i))+(b(ib+i)+b(id+i))) 2110 c(jh+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+i))+(b(ib+i)+b(id+i))) 2152 c(ja+j) = 2.0_wp * ( ( ( a(ia+i) + a(ie+i) ) + a(ic+i) ) + ( a(ib+i) + a(id+i) ) ) 2153 c(je+j) = 2.0_wp * ( ( ( a(ia+i) + a(ie+i) ) + a(ic+i) ) - ( a(ib+i) + a(id+i) ) ) 2154 c(jc+j) = 2.0_wp * ( ( ( a(ia+i) + a(ie+i) ) - a(ic+i) ) - ( b(ib+i) - b(id+i) ) ) 2155 c(jg+j) = 2.0_wp * ( ( ( a(ia+i) + a(ie+i) ) - a(ic+i) ) + ( b(ib+i) - b(id+i) ) ) 2156 c(jb+j) = 2.0_wp * ( ( a(ia+i) - a(ie+i) ) - b(ic+i) ) + ssin45 * ( ( a(ib+i) & 2157 - a(id+i) ) - ( b(ib+i) + b(id+i) ) ) 2158 c(jf+j) = 2.0_wp * ( ( a(ia+i) - a(ie+i) ) - b(ic+i) ) - ssin45 * ( ( a(ib+i) & 2159 - a(id+i) ) - ( b(ib+i) + b(id+i) ) ) 2160 c(jd+j) = 2.0_wp * ( ( a(ia+i) - a(ie+i) ) + b(ic+i) ) - ssin45 * ( ( a(ib+i) & 2161 - a(id+i) ) + ( b(ib+i) + b(id+i) ) ) 2162 c(jh+j) = 2.0_wp * ( ( a(ia+i) - a(ie+i) ) + b(ic+i) ) + ssin45 * ( ( a(ib+i) & 2163 - a(id+i) ) + ( b(ib+i) + b(id+i) ) ) 2111 2164 i = i + inc3 2112 2165 j = j + inc4 … … 2120 2173 END SUBROUTINE rpassm 2121 2174 2122 !------------------------------------------------------------------------------ !2175 !--------------------------------------------------------------------------------------------------! 2123 2176 ! Description: 2124 2177 ! ------------ 2125 2178 !> Computes factors of n & trigonometric functins required by fft99 & fft991cy 2126 2179 !> Method: Dimension trigs(n),ifax(1),jfax(10),lfax(7) 2127 !> subroutine 'set99' - computes factors of n & trigonometric 2128 !> functins required by fft99 & fft991cy 2129 !------------------------------------------------------------------------------! 2180 !> Subroutine 'set99' - computes factors of n & trigonometric functins required by fft99 & fft991cy 2181 !--------------------------------------------------------------------------------------------------! 2130 2182 SUBROUTINE set99( trigs, ifax, n ) 2131 2183 … … 2138 2190 IMPLICIT NONE 2139 2191 2140 INTEGER(iwp) :: n !< 2141 2142 INTEGER(iwp) :: ifax(*) !< 2143 REAL(wp) :: trigs(*) !< 2144 2145 2146 REAL(wp) :: angle !< 2147 REAL(wp) :: del !< 2148 INTEGER(iwp) :: i !< 2149 INTEGER(iwp) :: ifac !< 2150 INTEGER(iwp) :: ixxx !< 2151 INTEGER(iwp) :: k !< 2152 INTEGER(iwp) :: l !< 2153 INTEGER(iwp) :: nfax !< 2154 INTEGER(iwp) :: nhl !< 2155 INTEGER(iwp) :: nil !< 2156 INTEGER(iwp) :: nu !< 2157 2158 INTEGER(iwp) :: jfax(10) !< 2159 INTEGER(iwp) :: lfax(7) !< 2192 INTEGER(iwp) :: ifax(*) !< 2193 2194 INTEGER(iwp) :: i !< 2195 INTEGER(iwp) :: ifac !< 2196 INTEGER(iwp) :: ixxx !< 2197 INTEGER(iwp) :: jfax(10) !< 2198 INTEGER(iwp) :: k !< 2199 INTEGER(iwp) :: l !< 2200 INTEGER(iwp) :: lfax(7) !< 2201 INTEGER(iwp) :: n !< 2202 INTEGER(iwp) :: nfax !< 2203 INTEGER(iwp) :: nhl !< 2204 INTEGER(iwp) :: nil !< 2205 INTEGER(iwp) :: nu !< 2206 2207 REAL(wp) :: angle !< 2208 REAL(wp) :: del !< 2209 2210 REAL(wp) :: trigs(*) !< 2160 2211 2161 2212 DATA lfax / 6, 8, 5, 4, 3, 2, 1 / … … 2164 2215 ixxx = 1 2165 2216 2166 del = 4.0_wp * ASIN( 1.0_wp ) / REAL( n, KIND =wp )2217 del = 4.0_wp * ASIN( 1.0_wp ) / REAL( n, KIND = wp ) 2167 2218 nil = 0 2168 nhl = ( n/2) - 12219 nhl = ( n / 2 ) - 1 2169 2220 DO k = nil, nhl 2170 angle = REAL( k, KIND =wp ) * del2221 angle = REAL( k, KIND = wp ) * del 2171 2222 trigs(2*k+1) = COS( angle ) 2172 2223 trigs(2*k+2) = SIN( angle ) … … 2198 2249 ifac = lfax(l) 2199 2250 IF (ifac < 2 ) THEN 2200 message_string = 'number of gridpoints along x or/and y ' // & 2201 'contain illegal factors' // & 2251 message_string = 'number of gridpoints along x or/and y contain illegal factors' // & 2202 2252 '&only factors 2, 3, 5 are allowed' 2203 2253 CALL message( 'temperton_fft', 'PA0311', 1, 2, 0, 6, 0 ) … … 2224 2274 IMPLICIT NONE 2225 2275 2226 REAL(wp),DIMENSION(:,:) :: a !< 2227 REAL(wp),DIMENSION(:) :: trigs !< 2228 REAL(wp),DIMENSION(:,:) :: work !< 2276 2277 INTEGER(iwp) :: i !< 2278 INTEGER(iwp) :: ia !< 2279 INTEGER(iwp) :: ibase !< 2280 INTEGER(iwp) :: ierr !< 2281 INTEGER(iwp) :: ifac !< 2282 INTEGER(iwp) :: igo !< 2283 INTEGER(iwp) :: ii !< 2284 INTEGER(iwp) :: inc !< 2285 INTEGER(iwp) :: isign !< 2286 INTEGER(iwp) :: istart !< 2287 INTEGER(iwp) :: ix !< 2288 INTEGER(iwp) :: iz !< 2289 INTEGER(iwp) :: j !< 2290 INTEGER(iwp) :: jbase !< 2291 INTEGER(iwp) :: jump !< 2292 INTEGER(iwp) :: k !< 2293 INTEGER(iwp) :: la !< 2294 INTEGER(iwp) :: lot !< 2295 INTEGER(iwp) :: n !< 2296 INTEGER(iwp) :: nfax !< 2297 INTEGER(iwp) :: nvex !< 2298 INTEGER(iwp) :: nx !< 2299 2229 2300 INTEGER(iwp),DIMENSION(:),INTENT(IN) :: ifax !< 2230 2301 2231 INTEGER(iwp) :: inc !< 2232 INTEGER(iwp) :: isign !< 2233 INTEGER(iwp) :: jump !< 2234 INTEGER(iwp) :: lot !< 2235 INTEGER(iwp) :: n !< 2236 2237 2238 INTEGER(iwp) :: i !< 2239 INTEGER(iwp) :: ia !< 2240 INTEGER(iwp) :: ibase !< 2241 INTEGER(iwp) :: ierr !< 2242 INTEGER(iwp) :: ifac !< 2243 INTEGER(iwp) :: igo !< 2244 INTEGER(iwp) :: ii !< 2245 INTEGER(iwp) :: istart !< 2246 INTEGER(iwp) :: ix !< 2247 INTEGER(iwp) :: iz !< 2248 INTEGER(iwp) :: j !< 2249 INTEGER(iwp) :: jbase !< 2250 INTEGER(iwp) :: k !< 2251 INTEGER(iwp) :: la !< 2252 INTEGER(iwp) :: nfax !< 2253 INTEGER(iwp) :: nvex !< 2254 INTEGER(iwp) :: nx !< 2302 REAL(wp),DIMENSION(:,:) :: a !< 2303 REAL(wp),DIMENSION(:) :: trigs !< 2304 REAL(wp),DIMENSION(:,:) :: work !< 2255 2305 2256 2306 … … 2262 2312 nfax = ifax(1) 2263 2313 nx = n + 1 2264 IF ( MOD( n,2) == 1 ) nx = n2314 IF ( MOD( n, 2 ) == 1 ) nx = n 2265 2315 nvex = 1 2266 2316 … … 2273 2323 i = istart 2274 2324 a(:,i+inc) = 0.5_wp * a(:,i) 2275 IF ( MOD( n,2) /= 1 ) THEN2325 IF ( MOD( n, 2 ) /= 1 ) THEN 2276 2326 i = istart + n * inc 2277 2327 a(:,i) = 0.5_wp * a(:,i) … … 2313 2363 ! 2314 2364 !-- If necessary, copy results back to a 2315 IF ( MOD( nfax,2) /= 0 ) THEN2365 IF ( MOD( nfax, 2 ) /= 0 ) THEN 2316 2366 ibase = 1 2317 2367 jbase = ia … … 2326 2376 ! 2327 2377 !-- Fill in zeros at end 2328 ix = istart + n *inc2378 ix = istart + n * inc 2329 2379 a(:,ix) = 0.0_wp 2330 2380 a(:,ix+inc) = 0.0_wp … … 2370 2420 ! 2371 2421 !-- If necessary, copy results back to a 2372 IF ( MOD( nfax,2) /= 0 ) THEN2422 IF ( MOD( nfax, 2 ) /= 0 ) THEN 2373 2423 ibase = 1 2374 2424 jbase = ia … … 2387 2437 a(:,ix+inc) = 0.0_wp 2388 2438 2389 IF ( MOD( n,2) /= 1 ) THEN2390 iz = istart + ( n+1) * inc2439 IF ( MOD( n, 2 ) /= 1 ) THEN 2440 iz = istart + ( n + 1 ) * inc 2391 2441 a(:,iz) = 0.0_wp 2392 2442 ENDIF … … 2396 2446 END SUBROUTINE fft991cy_vec 2397 2447 2398 !------------------------------------------------------------------------------ !2448 !--------------------------------------------------------------------------------------------------! 2399 2449 ! Description: 2400 2450 ! ------------ 2401 !> Performs one pass through data as part of 2402 !> multiple real fft (fourier analysis) routine. 2451 !> Performs one pass through data as part of multiple real fft (fourier analysis) routine. 2403 2452 !> 2404 2453 !> Method: … … 2422 2471 !> 2 - ifac not catered for 2423 2472 !> 3 - ifac only catered for if la=n/ifac 2424 !------------------------------------------------------------------------------ !2473 !--------------------------------------------------------------------------------------------------! 2425 2474 SUBROUTINE qpassm_vec( a, b, c, d, trigs, inc1, inc2, n, ifac, la, ierr ) 2426 2475 … … 2429 2478 IMPLICIT NONE 2430 2479 2431 INTEGER(iwp),INTENT(IN) :: ifac !< 2432 INTEGER(iwp),INTENT(IN) :: inc1 !< 2433 INTEGER(iwp),INTENT(IN) :: inc2 !< 2434 INTEGER(iwp),INTENT(IN) :: la !< 2435 INTEGER(iwp),INTENT(IN) :: n !< 2436 INTEGER(iwp),INTENT(OUT) :: ierr !< 2480 INTEGER(iwp),INTENT(IN) :: ifac !< 2481 INTEGER(iwp),INTENT(IN) :: inc1 !< 2482 INTEGER(iwp),INTENT(IN) :: inc2 !< 2483 INTEGER(iwp),INTENT(IN) :: la !< 2484 INTEGER(iwp),INTENT(IN) :: n !< 2485 2486 INTEGER(iwp),INTENT(OUT) :: ierr !< 2487 2488 INTEGER(iwp) :: i !< 2489 INTEGER(iwp) :: ia !< 2490 INTEGER(iwp) :: ib !< 2491 INTEGER(iwp) :: ibase !< 2492 INTEGER(iwp) :: ic !< 2493 INTEGER(iwp) :: id !< 2494 INTEGER(iwp) :: ie !< 2495 INTEGER(iwp) :: if !< 2496 INTEGER(iwp) :: ig !< 2497 INTEGER(iwp) :: igo !< 2498 INTEGER(iwp) :: ih !< 2499 INTEGER(iwp) :: iink !< 2500 INTEGER(iwp) :: ijump !< 2501 INTEGER(iwp) :: ipl !< loop index parallel loop 2502 INTEGER(iwp) :: j !< 2503 INTEGER(iwp) :: ja !< 2504 INTEGER(iwp) :: jb !< 2505 INTEGER(iwp) :: jbase !< 2506 INTEGER(iwp) :: jc !< 2507 INTEGER(iwp) :: jd !< 2508 INTEGER(iwp) :: je !< 2509 INTEGER(iwp) :: jf !< 2510 INTEGER(iwp) :: jink !< 2511 INTEGER(iwp) :: k !< 2512 INTEGER(iwp) :: kb !< 2513 INTEGER(iwp) :: kc !< 2514 INTEGER(iwp) :: kd !< 2515 INTEGER(iwp) :: ke !< 2516 INTEGER(iwp) :: kf !< 2517 INTEGER(iwp) :: kstop !< 2518 INTEGER(iwp) :: l !< 2519 INTEGER(iwp) :: m !< 2437 2520 2438 2521 ! 2439 2522 !-- Arrays are dimensioned with n 2440 REAL(wp),DIMENSION(:,:) :: a !< 2441 REAL(wp),DIMENSION(:,:) :: b !< 2442 REAL(wp),DIMENSION(:,:) :: c !< 2443 REAL(wp),DIMENSION(:,:) :: d !< 2444 REAL(wp),DIMENSION(:),INTENT(IN) :: trigs !< 2445 2446 REAL(wp) :: a0 !< 2447 REAL(wp) :: a1 !< 2448 REAL(wp) :: a10 !< 2449 REAL(wp) :: a11 !< 2450 REAL(wp) :: a2 !< 2451 REAL(wp) :: a20 !< 2452 REAL(wp) :: a21 !< 2453 REAL(wp) :: a3 !< 2454 REAL(wp) :: a4 !< 2455 REAL(wp) :: a5 !< 2456 REAL(wp) :: a6 !< 2457 REAL(wp) :: b0 !< 2458 REAL(wp) :: b1 !< 2459 REAL(wp) :: b10 !< 2460 REAL(wp) :: b11 !< 2461 REAL(wp) :: b2 !< 2462 REAL(wp) :: b20 !< 2463 REAL(wp) :: b21 !< 2464 REAL(wp) :: b3 !< 2465 REAL(wp) :: b4 !< 2466 REAL(wp) :: b5 !< 2467 REAL(wp) :: b6 !< 2468 REAL(wp) :: c1 !< 2469 REAL(wp) :: c2 !< 2470 REAL(wp) :: c3 !< 2471 REAL(wp) :: c4 !< 2472 REAL(wp) :: c5 !< 2473 REAL(wp) :: qrt5 !< 2474 REAL(wp) :: s1 !< 2475 REAL(wp) :: s2 !< 2476 REAL(wp) :: s3 !< 2477 REAL(wp) :: s4 !< 2478 REAL(wp) :: s5 !< 2479 REAL(wp) :: sin36 !< 2480 REAL(wp) :: sin45 !< 2481 REAL(wp) :: sin60 !< 2482 REAL(wp) :: sin72 !< 2483 REAL(wp) :: z !< 2484 REAL(wp) :: zqrt5 !< 2485 REAL(wp) :: zsin36 !< 2486 REAL(wp) :: zsin45 !< 2487 REAL(wp) :: zsin60 !< 2488 REAL(wp) :: zsin72 !< 2489 2490 INTEGER(iwp) :: i !< 2491 INTEGER(iwp) :: ia !< 2492 INTEGER(iwp) :: ib !< 2493 INTEGER(iwp) :: ibase !< 2494 INTEGER(iwp) :: ic !< 2495 INTEGER(iwp) :: id !< 2496 INTEGER(iwp) :: ie !< 2497 INTEGER(iwp) :: if !< 2498 INTEGER(iwp) :: ig !< 2499 INTEGER(iwp) :: igo !< 2500 INTEGER(iwp) :: ih !< 2501 INTEGER(iwp) :: iink !< 2502 INTEGER(iwp) :: ijump !< 2503 INTEGER(iwp) :: ipl !< loop index parallel loop 2504 INTEGER(iwp) :: j !< 2505 INTEGER(iwp) :: ja !< 2506 INTEGER(iwp) :: jb !< 2507 INTEGER(iwp) :: jbase !< 2508 INTEGER(iwp) :: jc !< 2509 INTEGER(iwp) :: jd !< 2510 INTEGER(iwp) :: je !< 2511 INTEGER(iwp) :: jf !< 2512 INTEGER(iwp) :: jink !< 2513 INTEGER(iwp) :: k !< 2514 INTEGER(iwp) :: kb !< 2515 INTEGER(iwp) :: kc !< 2516 INTEGER(iwp) :: kd !< 2517 INTEGER(iwp) :: ke !< 2518 INTEGER(iwp) :: kf !< 2519 INTEGER(iwp) :: kstop !< 2520 INTEGER(iwp) :: l !< 2521 INTEGER(iwp) :: m !< 2523 REAL(wp),DIMENSION(:,:) :: a !< 2524 REAL(wp),DIMENSION(:,:) :: b !< 2525 REAL(wp),DIMENSION(:,:) :: c !< 2526 REAL(wp),DIMENSION(:,:) :: d !< 2527 2528 REAL(wp),DIMENSION(:),INTENT(IN) :: trigs !< 2529 2530 REAL(wp) :: a0 !< 2531 REAL(wp) :: a1 !< 2532 REAL(wp) :: a10 !< 2533 REAL(wp) :: a11 !< 2534 REAL(wp) :: a2 !< 2535 REAL(wp) :: a20 !< 2536 REAL(wp) :: a21 !< 2537 REAL(wp) :: a3 !< 2538 REAL(wp) :: a4 !< 2539 REAL(wp) :: a5 !< 2540 REAL(wp) :: a6 !< 2541 REAL(wp) :: b0 !< 2542 REAL(wp) :: b1 !< 2543 REAL(wp) :: b10 !< 2544 REAL(wp) :: b11 !< 2545 REAL(wp) :: b2 !< 2546 REAL(wp) :: b20 !< 2547 REAL(wp) :: b21 !< 2548 REAL(wp) :: b3 !< 2549 REAL(wp) :: b4 !< 2550 REAL(wp) :: b5 !< 2551 REAL(wp) :: b6 !< 2552 REAL(wp) :: c1 !< 2553 REAL(wp) :: c2 !< 2554 REAL(wp) :: c3 !< 2555 REAL(wp) :: c4 !< 2556 REAL(wp) :: c5 !< 2557 REAL(wp) :: qrt5 !< 2558 REAL(wp) :: s1 !< 2559 REAL(wp) :: s2 !< 2560 REAL(wp) :: s3 !< 2561 REAL(wp) :: s4 !< 2562 REAL(wp) :: s5 !< 2563 REAL(wp) :: sin36 !< 2564 REAL(wp) :: sin45 !< 2565 REAL(wp) :: sin60 !< 2566 REAL(wp) :: sin72 !< 2567 REAL(wp) :: z !< 2568 REAL(wp) :: zqrt5 !< 2569 REAL(wp) :: zsin36 !< 2570 REAL(wp) :: zsin45 !< 2571 REAL(wp) :: zsin60 !< 2572 REAL(wp) :: zsin72 !< 2573 2522 2574 2523 2575 … … 2531 2583 iink = la * inc1 2532 2584 jink = la * inc2 2533 ijump = ( ifac-1) * iink2534 kstop = ( n -ifac ) / ( 2*ifac )2585 ijump = ( ifac - 1 ) * iink 2586 kstop = ( n - ifac ) / ( 2 * ifac ) 2535 2587 2536 2588 ibase = 0 … … 2553 2605 ib = ia + iink 2554 2606 ja = 1 2555 jb = ja + ( 2*m-la) * inc22607 jb = ja + ( 2 * m - la ) * inc2 2556 2608 2557 2609 IF ( la /= m ) THEN 2558 2610 2559 2611 DO l = 1, la 2560 i = ibase2561 j = jbase2612 i = ibase 2613 j = jbase 2562 2614 c(:,ja+j) = a(:,ia+i) + a(:,ib+i) 2563 2615 c(:,jb+j) = a(:,ia+i) - a(:,ib+i) 2564 ibase = ibase + inc12565 jbase = jbase + inc22566 ENDDO 2567 2568 ja = ja + jink2569 jink = 2 * jink2570 jb = jb - jink2616 ibase = ibase + inc1 2617 jbase = jbase + inc2 2618 ENDDO 2619 2620 ja = ja + jink 2621 jink = 2 * jink 2622 jb = jb - jink 2571 2623 ibase = ibase + ijump 2572 2624 ijump = 2 * ijump + iink … … 2582 2634 i = ibase 2583 2635 j = jbase 2584 c(:,ja+j) = a(:,ia+i) + ( c1*a(:,ib+i)+s1*b(:,ib+i))2585 c(:,jb+j) = a(:,ia+i) - ( c1*a(:,ib+i)+s1*b(:,ib+i))2586 d(:,ja+j) = ( c1*b(:,ib+i)-s1*a(:,ib+i)) + b(:,ia+i)2587 d(:,jb+j) = ( c1*b(:,ib+i)-s1*a(:,ib+i)) - b(:,ia+i)2636 c(:,ja+j) = a(:,ia+i) + ( c1 * a(:,ib+i) + s1 * b(:,ib+i) ) 2637 c(:,jb+j) = a(:,ia+i) - ( c1 * a(:,ib+i) + s1 * b(:,ib+i) ) 2638 d(:,ja+j) = ( c1 * b(:,ib+i) - s1 * a(:,ib+i) ) + b(:,ia+i) 2639 d(:,jb+j) = ( c1 * b(:,ib+i) - s1 * a(:,ib+i) ) - b(:,ia+i) 2588 2640 ibase = ibase + inc1 2589 2641 jbase = jbase + inc2 … … 2611 2663 ELSE 2612 2664 2613 z = 1.0_wp /REAL(n,KIND=wp)2614 DO l = 1, la 2615 i = ibase 2616 j = jbase 2617 c(:,ja+j) = z *(a(:,ia+i)+a(:,ib+i))2618 c(:,jb+j) = z *(a(:,ia+i)-a(:,ib+i))2665 z = 1.0_wp / REAL( n, KIND = wp ) 2666 DO l = 1, la 2667 i = ibase 2668 j = jbase 2669 c(:,ja+j) = z * ( a(:,ia+i) + a(:,ib+i) ) 2670 c(:,jb+j) = z * ( a(:,ia+i) - a(:,ib+i) ) 2619 2671 ibase = ibase + inc1 2620 2672 jbase = jbase + inc2 … … 2631 2683 ic = ib + iink 2632 2684 ja = 1 2633 jb = ja + ( 2*m-la) * inc22685 jb = ja + ( 2 * m - la ) * inc2 2634 2686 jc = jb 2635 2687 … … 2637 2689 2638 2690 DO l = 1, la 2639 i = ibase2640 j = jbase2641 c(:,ja+j) = a(:,ia+i) + ( a(:,ib+i)+a(:,ic+i))2642 c(:,jb+j) = a(:,ia+i) - 0.5_wp *(a(:,ib+i)+a(:,ic+i))2643 d(:,jb+j) = sin60 *(a(:,ic+i)-a(:,ib+i))2644 ibase = ibase + inc12645 jbase = jbase + inc22646 ENDDO 2647 2648 ja = ja + jink2649 jink = 2 * jink2650 jb = jb + jink2651 jc = jc - jink2691 i = ibase 2692 j = jbase 2693 c(:,ja+j) = a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) ) 2694 c(:,jb+j) = a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) 2695 d(:,jb+j) = sin60 * ( a(:,ic+i) - a(:,ib+i) ) 2696 ibase = ibase + inc1 2697 jbase = jbase + inc2 2698 ENDDO 2699 2700 ja = ja + jink 2701 jink = 2 * jink 2702 jb = jb + jink 2703 jc = jc - jink 2652 2704 ibase = ibase + ijump 2653 2705 ijump = 2 * ijump + iink … … 2667 2719 i = ibase 2668 2720 j = jbase 2669 DO ipl = 1, SIZE(a,1) 2670 a1 = (c1*a(ipl,ib+i)+s1*b(ipl,ib+i)) + (c2*a(ipl,ic+i)+s2*b(ipl,ic+i)) 2671 b1 = (c1*b(ipl,ib+i)-s1*a(ipl,ib+i)) + (c2*b(ipl,ic+i)-s2*a(ipl,ic+i)) 2672 a2 = a(ipl,ia+i) - 0.5_wp*a1 2673 b2 = b(ipl,ia+i) - 0.5_wp*b1 2674 a3 = sin60*((c1*a(ipl,ib+i)+s1*b(ipl,ib+i))-(c2*a(ipl,ic+i)+s2*b(ipl,ic+i))) 2675 b3 = sin60*((c1*b(ipl,ib+i)-s1*a(ipl,ib+i))-(c2*b(ipl,ic+i)-s2*a(ipl,ic+i))) 2721 DO ipl = 1, SIZE( a, 1 ) 2722 a1 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) ) & 2723 + ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) ) 2724 b1 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) ) & 2725 + ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) ) 2726 a2 = a(ipl,ia+i) - 0.5_wp * a1 2727 b2 = b(ipl,ia+i) - 0.5_wp * b1 2728 a3 = sin60 * ( ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) ) & 2729 - ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) ) ) 2730 b3 = sin60 * ( ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) ) & 2731 - ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) ) ) 2676 2732 c(ipl,ja+j) = a(ipl,ia+i) + a1 2677 2733 d(ipl,ja+j) = b(ipl,ia+i) + b1 … … 2679 2735 d(ipl,jb+j) = b2 - a3 2680 2736 c(ipl,jc+j) = a2 - b3 2681 d(ipl,jc+j) = - (b2+a3)2737 d(ipl,jc+j) = - ( b2 + a3 ) 2682 2738 ENDDO 2683 2739 ibase = ibase + inc1 … … 2685 2741 ENDDO 2686 2742 ibase = ibase + ijump 2687 ja = ja + jink2688 jb = jb + jink2689 jc = jc - jink2743 ja = ja + jink 2744 jb = jb + jink 2745 jc = jc - jink 2690 2746 ENDDO 2691 2747 … … 2698 2754 i = ibase 2699 2755 j = jbase 2700 c(:,ja+j) = a(:,ia+i) + 0.5_wp *(a(:,ib+i)-a(:,ic+i))2701 d(:,ja+j) = - sin60*(a(:,ib+i)+a(:,ic+i))2702 c(:,jb+j) = a(:,ia+i) - ( a(:,ib+i)-a(:,ic+i))2756 c(:,ja+j) = a(:,ia+i) + 0.5_wp * ( a(:,ib+i) - a(:,ic+i) ) 2757 d(:,ja+j) = - sin60 * ( a(:,ib+i) + a(:,ic+i) ) 2758 c(:,jb+j) = a(:,ia+i) - ( a(:,ib+i) - a(:,ic+i) ) 2703 2759 ibase = ibase + inc1 2704 2760 jbase = jbase + inc2 … … 2707 2763 ELSE 2708 2764 2709 z = 1.0_wp / REAL( n, KIND =wp )2710 zsin60 = z *sin602711 DO l = 1, la 2712 i = ibase 2713 j = jbase 2714 c(:,ja+j) = z *(a(:,ia+i)+(a(:,ib+i)+a(:,ic+i)))2715 c(:,jb+j) = z *(a(:,ia+i)-0.5_wp*(a(:,ib+i)+a(:,ic+i)))2716 d(:,jb+j) = zsin60 *(a(:,ic+i)-a(:,ib+i))2765 z = 1.0_wp / REAL( n, KIND = wp ) 2766 zsin60 = z * sin60 2767 DO l = 1, la 2768 i = ibase 2769 j = jbase 2770 c(:,ja+j) = z * ( a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) ) ) 2771 c(:,jb+j) = z * ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) ) 2772 d(:,jb+j) = zsin60 * ( a(:,ic+i) - a(:,ib+i) ) 2717 2773 ibase = ibase + inc1 2718 2774 jbase = jbase + inc2 … … 2730 2786 id = ic + iink 2731 2787 ja = 1 2732 jb = ja + ( 2*m-la) * inc22733 jc = jb + 2 *m*inc22788 jb = ja + ( 2 * m - la ) * inc2 2789 jc = jb + 2 * m * inc2 2734 2790 jd = jb 2735 2791 … … 2737 2793 2738 2794 DO l = 1, la 2739 2740 i = ibase 2741 j = jbase 2742 c(:,ja+j) = (a(:,ia+i)+a(:,ic+i)) + (a(:,ib+i)+a(:,id+i)) 2743 c(:,jc+j) = (a(:,ia+i)+a(:,ic+i)) - (a(:,ib+i)+a(:,id+i)) 2795 i = ibase 2796 j = jbase 2797 c(:,ja+j) = ( a(:,ia+i) + a(:,ic+i) ) + ( a(:,ib+i) + a(:,id+i) ) 2798 c(:,jc+j) = ( a(:,ia+i) + a(:,ic+i) ) - ( a(:,ib+i) + a(:,id+i) ) 2744 2799 c(:,jb+j) = a(:,ia+i) - a(:,ic+i) 2745 2800 d(:,jb+j) = a(:,id+i) - a(:,ib+i) … … 2748 2803 ENDDO 2749 2804 2750 ja = ja + jink2751 jink = 2 * jink2752 jb = jb + jink2753 jc = jc - jink2754 jd = jd - jink2805 ja = ja + jink 2806 jink = 2 * jink 2807 jb = jb + jink 2808 jc = jc - jink 2809 jd = jd - jink 2755 2810 ibase = ibase + ijump 2756 2811 ijump = 2 * ijump + iink … … 2772 2827 i = ibase 2773 2828 j = jbase 2774 DO ipl = 1, SIZE(a,1) 2775 a0 = a(ipl,ia+i) + (c2*a(ipl,ic+i)+s2*b(ipl,ic+i)) 2776 a2 = a(ipl,ia+i) - (c2*a(ipl,ic+i)+s2*b(ipl,ic+i)) 2777 a1 = (c1*a(ipl,ib+i)+s1*b(ipl,ib+i)) + (c3*a(ipl,id+i)+s3*b(ipl,id+i)) 2778 a3 = (c1*a(ipl,ib+i)+s1*b(ipl,ib+i)) - (c3*a(ipl,id+i)+s3*b(ipl,id+i)) 2779 b0 = b(ipl,ia+i) + (c2*b(ipl,ic+i)-s2*a(ipl,ic+i)) 2780 b2 = b(ipl,ia+i) - (c2*b(ipl,ic+i)-s2*a(ipl,ic+i)) 2781 b1 = (c1*b(ipl,ib+i)-s1*a(ipl,ib+i)) + (c3*b(ipl,id+i)-s3*a(ipl,id+i)) 2782 b3 = (c1*b(ipl,ib+i)-s1*a(ipl,ib+i)) - (c3*b(ipl,id+i)-s3*a(ipl,id+i)) 2829 DO ipl = 1, SIZE( a, 1 ) 2830 a0 = a(ipl,ia+i) + ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) ) 2831 a2 = a(ipl,ia+i) - ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) ) 2832 a1 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) ) & 2833 + ( c3 * a(ipl,id+i) + s3 * b(ipl,id+i) ) 2834 a3 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) ) & 2835 - ( c3 * a(ipl,id+i) + s3 * b(ipl,id+i) ) 2836 b0 = b(ipl,ia+i) + ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) ) 2837 b2 = b(ipl,ia+i) - ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) ) 2838 b1 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) ) & 2839 + ( c3 * b(ipl,id+i) - s3 * a(ipl,id+i) ) 2840 b3 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) ) & 2841 - ( c3 * b(ipl,id+i) - s3 * a(ipl,id+i) ) 2783 2842 c(ipl,ja+j) = a0 + a1 2784 2843 c(ipl,jc+j) = a0 - a1 … … 2788 2847 c(ipl,jd+j) = a2 - b3 2789 2848 d(ipl,jb+j) = b2 - a3 2790 d(ipl,jd+j) = - (b2+a3)2849 d(ipl,jd+j) = - ( b2 + a3 ) 2791 2850 ENDDO 2792 2851 ibase = ibase + inc1 … … 2809 2868 i = ibase 2810 2869 j = jbase 2811 c(:,ja+j) = a(:,ia+i) + sin45*(a(:,ib+i)-a(:,id+i))2812 c(:,jb+j) = a(:,ia+i) - sin45*(a(:,ib+i)-a(:,id+i))2813 d(:,ja+j) = - a(:,ic+i) - sin45*(a(:,ib+i)+a(:,id+i))2814 d(:,jb+j) = a(:,ic+i) - sin45*(a(:,ib+i)+a(:,id+i))2870 c(:,ja+j) = a(:,ia+i) + sin45 * ( a(:,ib+i) - a(:,id+i) ) 2871 c(:,jb+j) = a(:,ia+i) - sin45 * ( a(:,ib+i) - a(:,id+i) ) 2872 d(:,ja+j) = - a(:,ic+i) - sin45 * ( a(:,ib+i) + a(:,id+i) ) 2873 d(:,jb+j) = a(:,ic+i) - sin45 * ( a(:,ib+i) + a(:,id+i) ) 2815 2874 ibase = ibase + inc1 2816 2875 jbase = jbase + inc2 … … 2819 2878 ELSE 2820 2879 2821 z = 1.0_wp / REAL( n, KIND =wp )2822 DO l = 1, la 2823 i = ibase 2824 j = jbase 2825 c(:,ja+j) = z *((a(:,ia+i)+a(:,ic+i))+(a(:,ib+i)+a(:,id+i)))2826 c(:,jc+j) = z *((a(:,ia+i)+a(:,ic+i))-(a(:,ib+i)+a(:,id+i)))2827 c(:,jb+j) = z *(a(:,ia+i)-a(:,ic+i))2828 d(:,jb+j) = z *(a(:,id+i)-a(:,ib+i))2880 z = 1.0_wp / REAL( n, KIND = wp ) 2881 DO l = 1, la 2882 i = ibase 2883 j = jbase 2884 c(:,ja+j) = z * ( ( a(:,ia+i) + a(:,ic+i) ) + ( a(:,ib+i) + a(:,id+i) ) ) 2885 c(:,jc+j) = z * ( ( a(:,ia+i) + a(:,ic+i) ) - ( a(:,ib+i) + a(:,id+i) ) ) 2886 c(:,jb+j) = z * ( a(:,ia+i) - a(:,ic+i) ) 2887 d(:,jb+j) = z * ( a(:,id+i) - a(:,ib+i) ) 2829 2888 ibase = ibase + inc1 2830 2889 jbase = jbase + inc2 … … 2843 2902 ie = id + iink 2844 2903 ja = 1 2845 jb = ja + ( 2*m-la) * inc22846 jc = jb + 2 *m*inc22904 jb = ja + ( 2 * m - la ) * inc2 2905 jc = jb + 2 * m * inc2 2847 2906 jd = jc 2848 2907 je = jb … … 2853 2912 i = ibase 2854 2913 j = jbase 2855 DO ipl = 1, SIZE( a,1)2914 DO ipl = 1, SIZE( a, 1 ) 2856 2915 a1 = a(ipl,ib+i) + a(ipl,ie+i) 2857 2916 a3 = a(ipl,ib+i) - a(ipl,ie+i) 2858 2917 a2 = a(ipl,ic+i) + a(ipl,id+i) 2859 2918 a4 = a(ipl,ic+i) - a(ipl,id+i) 2860 a5 = a(ipl,ia+i) - 0.25_wp *(a1+a2)2861 a6 = qrt5 *(a1-a2)2862 c(ipl,ja+j) = a(ipl,ia+i) + ( a1+a2)2919 a5 = a(ipl,ia+i) - 0.25_wp * ( a1 + a2 ) 2920 a6 = qrt5 * ( a1 - a2 ) 2921 c(ipl,ja+j) = a(ipl,ia+i) + ( a1 + a2 ) 2863 2922 c(ipl,jb+j) = a5 + a6 2864 2923 c(ipl,jc+j) = a5 - a6 2865 d(ipl,jb+j) = - sin72*a3 - sin36*a42866 d(ipl,jc+j) = - sin36*a3 + sin72*a42867 ENDDO 2868 ibase = ibase + inc1 2869 jbase = jbase + inc2 2870 ENDDO 2871 2872 ja = ja + jink2873 jink = 2 * jink2874 jb = jb + jink2875 jc = jc + jink2876 jd = jd - jink2877 je = je - jink2924 d(ipl,jb+j) = - sin72 * a3 - sin36 * a4 2925 d(ipl,jc+j) = - sin36 * a3 + sin72 * a4 2926 ENDDO 2927 ibase = ibase + inc1 2928 jbase = jbase + inc2 2929 ENDDO 2930 2931 ja = ja + jink 2932 jink = 2 * jink 2933 jb = jb + jink 2934 jc = jc + jink 2935 jd = jd - jink 2936 je = je - jink 2878 2937 ibase = ibase + ijump 2879 2938 ijump = 2 * ijump + iink … … 2899 2958 j = jbase 2900 2959 DO ipl = 1, SIZE(a,1) 2901 a1 = (c1*a(ipl,ib+i)+s1*b(ipl,ib+i)) + (c4*a(ipl,ie+i)+s4*b(ipl,ie+i)) 2902 a3 = (c1*a(ipl,ib+i)+s1*b(ipl,ib+i)) - (c4*a(ipl,ie+i)+s4*b(ipl,ie+i)) 2903 a2 = (c2*a(ipl,ic+i)+s2*b(ipl,ic+i)) + (c3*a(ipl,id+i)+s3*b(ipl,id+i)) 2904 a4 = (c2*a(ipl,ic+i)+s2*b(ipl,ic+i)) - (c3*a(ipl,id+i)+s3*b(ipl,id+i)) 2905 b1 = (c1*b(ipl,ib+i)-s1*a(ipl,ib+i)) + (c4*b(ipl,ie+i)-s4*a(ipl,ie+i)) 2906 b3 = (c1*b(ipl,ib+i)-s1*a(ipl,ib+i)) - (c4*b(ipl,ie+i)-s4*a(ipl,ie+i)) 2907 b2 = (c2*b(ipl,ic+i)-s2*a(ipl,ic+i)) + (c3*b(ipl,id+i)-s3*a(ipl,id+i)) 2908 b4 = (c2*b(ipl,ic+i)-s2*a(ipl,ic+i)) - (c3*b(ipl,id+i)-s3*a(ipl,id+i)) 2909 a5 = a(ipl,ia+i) - 0.25_wp*(a1+a2) 2910 a6 = qrt5*(a1-a2) 2911 b5 = b(ipl,ia+i) - 0.25_wp*(b1+b2) 2912 b6 = qrt5*(b1-b2) 2960 a1 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) ) + ( c4 * a(ipl,ie+i) + s4 & 2961 * b(ipl,ie+i) ) 2962 a3 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) ) - ( c4 * a(ipl,ie+i) + s4 & 2963 * b(ipl,ie+i) ) 2964 a2 = ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) ) + ( c3 * a(ipl,id+i) + s3 & 2965 * b(ipl,id+i) ) 2966 a4 = ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) ) - ( c3 * a(ipl,id+i) + s3 & 2967 * b(ipl,id+i) ) 2968 2969 b1 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) ) + ( c4 * b(ipl,ie+i) - s4 & 2970 * a(ipl,ie+i) ) 2971 b3 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) ) - ( c4 * b(ipl,ie+i) - s4 & 2972 * a(ipl,ie+i) ) 2973 b2 = ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) ) + ( c3 * b(ipl,id+i) - s3 & 2974 * a(ipl,id+i) ) 2975 b4 = ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) ) - ( c3 * b(ipl,id+i) - s3 & 2976 * a(ipl,id+i) ) 2977 2978 a5 = a(ipl,ia+i) - 0.25_wp * ( a1 + a2 ) 2979 a6 = qrt5 * ( a1 - a2 ) 2980 b5 = b(ipl,ia+i) - 0.25_wp * ( b1 + b2 ) 2981 b6 = qrt5 * ( b1 - b2) 2982 2913 2983 a10 = a5 + a6 2914 2984 a20 = a5 - a6 2915 2985 b10 = b5 + b6 2916 2986 b20 = b5 - b6 2917 a11 = sin72*b3 + sin36*b4 2918 a21 = sin36*b3 - sin72*b4 2919 b11 = sin72*a3 + sin36*a4 2920 b21 = sin36*a3 - sin72*a4 2921 c(ipl,ja+j) = a(ipl,ia+i) + (a1+a2) 2987 a11 = sin72 * b3 + sin36 * b4 2988 a21 = sin36 * b3 - sin72 * b4 2989 2990 b11 = sin72 * a3 + sin36 * a4 2991 b21 = sin36 * a3 - sin72 * a4 2992 2993 c(ipl,ja+j) = a(ipl,ia+i) + ( a1 + a2 ) 2922 2994 c(ipl,jb+j) = a10 + a11 2923 2995 c(ipl,je+j) = a10 - a11 2924 2996 c(ipl,jc+j) = a20 + a21 2925 2997 c(ipl,jd+j) = a20 - a21 2926 d(ipl,ja+j) = b(ipl,ia+i) + (b1+b2) 2998 2999 d(ipl,ja+j) = b(ipl,ia+i) + ( b1 + b2 ) 2927 3000 d(ipl,jb+j) = b10 - b11 2928 d(ipl,je+j) = - (b10+b11)3001 d(ipl,je+j) = - ( b10 + b11 ) 2929 3002 d(ipl,jc+j) = b20 - b21 2930 d(ipl,jd+j) = - (b20+b21)3003 d(ipl,jd+j) = - ( b20 + b21 ) 2931 3004 ENDDO 2932 3005 ibase = ibase + inc1 … … 2934 3007 ENDDO 2935 3008 ibase = ibase + ijump 2936 ja = ja + jink2937 jb = jb + jink2938 jc = jc + jink2939 jd = jd - jink2940 je = je - jink3009 ja = ja + jink 3010 jb = jb + jink 3011 jc = jc + jink 3012 jd = jd - jink 3013 je = je - jink 2941 3014 ENDDO 2942 3015 … … 2949 3022 i = ibase 2950 3023 j = jbase 2951 DO ipl = 1, SIZE( a,1)3024 DO ipl = 1, SIZE( a, 1 ) 2952 3025 a1 = a(ipl,ib+i) + a(ipl,ie+i) 2953 3026 a3 = a(ipl,ib+i) - a(ipl,ie+i) 2954 3027 a2 = a(ipl,ic+i) + a(ipl,id+i) 2955 3028 a4 = a(ipl,ic+i) - a(ipl,id+i) 2956 a5 = a(ipl,ia+i) + 0.25_wp*(a3-a4) 2957 a6 = qrt5*(a3+a4) 3029 a5 = a(ipl,ia+i) + 0.25_wp * ( a3 - a4 ) 3030 a6 = qrt5 * ( a3 + a4 ) 3031 2958 3032 c(ipl,ja+j) = a5 + a6 2959 3033 c(ipl,jb+j) = a5 - a6 2960 c(ipl,jc+j) = a(ipl,ia+i) - (a3-a4) 2961 d(ipl,ja+j) = -sin36*a1 - sin72*a2 2962 d(ipl,jb+j) = -sin72*a1 + sin36*a2 3034 c(ipl,jc+j) = a(ipl,ia+i) - ( a3 - a4 ) 3035 3036 d(ipl,ja+j) = - sin36 * a1 - sin72 * a2 3037 d(ipl,jb+j) = - sin72 * a1 + sin36 * a2 2963 3038 ENDDO 2964 3039 ibase = ibase + inc1 … … 2968 3043 ELSE 2969 3044 2970 z = 1.0_wp / REAL( n, KIND =wp )3045 z = 1.0_wp / REAL( n, KIND = wp ) 2971 3046 zqrt5 = z * qrt5 2972 3047 zsin36 = z * sin36 … … 2975 3050 i = ibase 2976 3051 j = jbase 2977 DO ipl = 1, SIZE( a,1)3052 DO ipl = 1, SIZE( a, 1 ) 2978 3053 a1 = a(ipl,ib+i) + a(ipl,ie+i) 2979 3054 a3 = a(ipl,ib+i) - a(ipl,ie+i) 2980 3055 a2 = a(ipl,ic+i) + a(ipl,id+i) 2981 3056 a4 = a(ipl,ic+i) - a(ipl,id+i) 2982 a5 = z*(a(ipl,ia+i)-0.25_wp*(a1+a2)) 2983 a6 = zqrt5*(a1-a2) 2984 c(ipl,ja+j) = z*(a(ipl,ia+i)+(a1+a2)) 3057 a5 = z * ( a(ipl,ia+i) - 0.25_wp * ( a1 + a2 ) ) 3058 a6 = zqrt5 * ( a1 - a2 ) 3059 3060 c(ipl,ja+j) = z * ( a(ipl,ia+i) + ( a1 + a2 ) ) 2985 3061 c(ipl,jb+j) = a5 + a6 2986 3062 c(ipl,jc+j) = a5 - a6 2987 d(ipl,jb+j) = -zsin72*a3 - zsin36*a4 2988 d(ipl,jc+j) = -zsin36*a3 + zsin72*a4 3063 3064 d(ipl,jb+j) = - zsin72 * a3 - zsin36 * a4 3065 d(ipl,jc+j) = - zsin36 * a3 + zsin72 * a4 2989 3066 ENDDO 2990 3067 ibase = ibase + inc1 … … 3005 3082 if = ie + iink 3006 3083 ja = 1 3007 jb = ja + ( 2*m-la) * inc23008 jc = jb + 2 *m*inc23009 jd = jc + 2 *m*inc23084 jb = ja + ( 2 * m - la ) * inc2 3085 jc = jb + 2 * m * inc2 3086 jd = jc + 2 * m * inc2 3010 3087 je = jc 3011 3088 jf = jb … … 3016 3093 i = ibase 3017 3094 j = jbase 3018 DO ipl = 1, SIZE(a,1) 3019 a11 = (a(ipl,ic+i)+a(ipl,if+i)) + (a(ipl,ib+i)+a(ipl,ie+i)) 3020 c(ipl,ja+j) = (a(ipl,ia+i)+a(ipl,id+i)) + a11 3021 c(ipl,jc+j) = (a(ipl,ia+i)+a(ipl,id+i)-0.5_wp*a11) 3022 d(ipl,jc+j) = sin60*((a(ipl,ic+i)+a(ipl,if+i))-(a(ipl,ib+i)+a(ipl,ie+i))) 3023 a11 = (a(ipl,ic+i)-a(ipl,if+i)) + (a(ipl,ie+i)-a(ipl,ib+i)) 3024 c(ipl,jb+j) = (a(ipl,ia+i)-a(ipl,id+i)) - 0.5_wp*a11 3025 d(ipl,jb+j) = sin60*((a(ipl,ie+i)-a(ipl,ib+i))-(a(ipl,ic+i)-a(ipl,if+i))) 3026 c(ipl,jd+j) = (a(ipl,ia+i)-a(ipl,id+i)) + a11 3095 DO ipl = 1, SIZE( a, 1 ) 3096 a11 = ( a(ipl,ic+i) + a(ipl,if+i) ) + ( a(ipl,ib+i) + a(ipl,ie+i) ) 3097 c(ipl,ja+j) = ( a(ipl,ia+i) + a(ipl,id+i) ) + a11 3098 c(ipl,jc+j) = ( a(ipl,ia+i) + a(ipl,id+i) - 0.5_wp * a11 ) 3099 d(ipl,jc+j) = sin60 * ( ( a(ipl,ic+i) + a(ipl,if+i) ) - & 3100 ( a(ipl,ib+i) + a(ipl,ie+i) ) ) 3101 a11 = ( a(ipl,ic+i) - a(ipl,if+i) ) + ( a(ipl,ie+i) - a(ipl,ib+i) ) 3102 c(ipl,jb+j) = ( a(ipl,ia+i) - a(ipl,id+i) ) - 0.5_wp * a11 3103 d(ipl,jb+j) = sin60 * ( ( a(ipl,ie+i) - a(ipl,ib+i) ) - & 3104 ( a(ipl,ic+i) - a(ipl,if+i) ) ) 3105 c(ipl,jd+j) = ( a(ipl,ia+i) - a(ipl,id+i) ) + a11 3027 3106 END DO 3028 3107 ibase = ibase + inc1 … … 3061 3140 i = ibase 3062 3141 j = jbase 3063 DO ipl = 1, SIZE( a,1)3064 a1 = c1 *a(ipl,ib+i) + s1*b(ipl,ib+i)3065 b1 = c1 *b(ipl,ib+i) - s1*a(ipl,ib+i)3066 a2 = c2 *a(ipl,ic+i) + s2*b(ipl,ic+i)3067 b2 = c2 *b(ipl,ic+i) - s2*a(ipl,ic+i)3068 a3 = c3 *a(ipl,id+i) + s3*b(ipl,id+i)3069 b3 = c3 *b(ipl,id+i) - s3*a(ipl,id+i)3070 a4 = c4 *a(ipl,ie+i) + s4*b(ipl,ie+i)3071 b4 = c4 *b(ipl,ie+i) - s4*a(ipl,ie+i)3072 a5 = c5 *a(ipl,if+i) + s5*b(ipl,if+i)3073 b5 = c5 *b(ipl,if+i) - s5*a(ipl,if+i)3074 a11 = ( a2+a5) + (a1+a4)3075 a20 = ( a(ipl,ia+i)+a3) - 0.5_wp*a113076 a21 = sin60 *((a2+a5)-(a1+a4))3077 b11 = ( b2+b5) + (b1+b4)3078 b20 = ( b(ipl,ia+i)+b3) - 0.5_wp*b113079 b21 = sin60 *((b2+b5)-(b1+b4))3080 c(ipl,ja+j) = ( a(ipl,ia+i)+a3) + a113081 d(ipl,ja+j) = ( b(ipl,ia+i)+b3) + b113142 DO ipl = 1, SIZE( a, 1 ) 3143 a1 = c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) 3144 b1 = c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) 3145 a2 = c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) 3146 b2 = c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) 3147 a3 = c3 * a(ipl,id+i) + s3 * b(ipl,id+i) 3148 b3 = c3 * b(ipl,id+i) - s3 * a(ipl,id+i) 3149 a4 = c4 * a(ipl,ie+i) + s4 * b(ipl,ie+i) 3150 b4 = c4 * b(ipl,ie+i) - s4 * a(ipl,ie+i) 3151 a5 = c5 * a(ipl,if+i) + s5 * b(ipl,if+i) 3152 b5 = c5 * b(ipl,if+i) - s5 * a(ipl,if+i) 3153 a11 = ( a2 + a5 ) + ( a1 + a4 ) 3154 a20 = ( a(ipl,ia+i) + a3 ) - 0.5_wp * a11 3155 a21 = sin60 * ( ( a2 + a5 ) - ( a1 + a4 ) ) 3156 b11 = ( b2 + b5 ) + ( b1 + b4 ) 3157 b20 = ( b(ipl,ia+i) + b3 ) - 0.5_wp * b11 3158 b21 = sin60 * ( ( b2 + b5 ) - ( b1 + b4 ) ) 3159 c(ipl,ja+j) = ( a(ipl,ia+i) + a3 ) + a11 3160 d(ipl,ja+j) = ( b(ipl,ia+i) + b3 ) + b11 3082 3161 c(ipl,jc+j) = a20 - b21 3083 3162 d(ipl,jc+j) = a21 + b20 3084 3163 c(ipl,je+j) = a20 + b21 3085 3164 d(ipl,je+j) = a21 - b20 3086 a11 = ( a2-a5) + (a4-a1)3087 a20 = ( a(ipl,ia+i)-a3) - 0.5_wp*a113088 a21 = sin60 *((a4-a1)-(a2-a5))3089 b11 = ( b5-b2) - (b4-b1)3090 b20 = ( b3-b(ipl,ia+i)) - 0.5_wp*b113091 b21 = sin60 *((b5-b2)+(b4-b1))3165 a11 = ( a2 - a5 ) + ( a4 - a1 ) 3166 a20 = ( a(ipl,ia+i) - a3 ) - 0.5_wp * a11 3167 a21 = sin60 * ( ( a4 - a1 ) - ( a2 - a5 ) ) 3168 b11 = ( b5 - b2 ) - ( b4 - b1 ) 3169 b20 = ( b3 - b(ipl,ia+i) ) - 0.5_wp * b11 3170 b21 = sin60 * ( ( b5 - b2 ) + ( b4 - b1 ) ) 3092 3171 c(ipl,jb+j) = a20 - b21 3093 3172 d(ipl,jb+j) = a21 - b20 3094 c(ipl,jd+j) = a11 + ( a(ipl,ia+i)-a3)3095 d(ipl,jd+j) = b11 + ( b3-b(ipl,ia+i))3173 c(ipl,jd+j) = a11 + ( a(ipl,ia+i) - a3 ) 3174 d(ipl,jd+j) = b11 + ( b3 - b(ipl,ia+i) ) 3096 3175 c(ipl,jf+j) = a20 + b21 3097 3176 d(ipl,jf+j) = a21 + b20 … … 3117 3196 i = ibase 3118 3197 j = jbase 3119 c(:,ja+j) = (a(:,ia+i)+0.5_wp*(a(:,ic+i)-a(:,ie+i))) + sin60*(a(:,ib+i)-a(:,if+i)) 3120 d(:,ja+j) = -(a(:,id+i)+0.5_wp*(a(:,ib+i)+a(:,if+i))) - sin60*(a(:,ic+i)+a(:,ie+i)) 3121 c(:,jb+j) = a(:,ia+i) - (a(:,ic+i)-a(:,ie+i)) 3122 d(:,jb+j) = a(:,id+i) - (a(:,ib+i)+a(:,if+i)) 3123 c(:,jc+j) = (a(:,ia+i)+0.5_wp*(a(:,ic+i)-a(:,ie+i))) - sin60*(a(:,ib+i)-a(:,if+i)) 3124 d(:,jc+j) = -(a(:,id+i)+0.5_wp*(a(:,ib+i)+a(:,if+i))) + sin60*(a(:,ic+i)+a(:,ie+i)) 3198 c(:,ja+j) = ( a(:,ia+i) + 0.5_wp * ( a(:,ic+i) - a(:,ie+i) ) ) + sin60 & 3199 * ( a(:,ib+i) - a(:,if+i) ) 3200 d(:,ja+j) = - ( a(:,id+i) + 0.5_wp * ( a(:,ib+i) + a(:,if+i) ) ) - sin60 & 3201 * ( a(:,ic+i) + a(:,ie+i) ) 3202 c(:,jb+j) = a(:,ia+i) - ( a(:,ic+i) - a(:,ie+i) ) 3203 d(:,jb+j) = a(:,id+i) - ( a(:,ib+i) + a(:,if+i) ) 3204 c(:,jc+j) = ( a(:,ia+i) + 0.5_wp * ( a(:,ic+i) - a(:,ie+i) ) ) - sin60 & 3205 * ( a(:,ib+i) - a(:,if+i) ) 3206 d(:,jc+j) = - ( a(:,id+i) + 0.5_wp * ( a(:,ib+i) + a(:,if+i) ) ) + sin60 & 3207 * ( a(:,ic+i) + a(:,ie+i) ) 3125 3208 ibase = ibase + inc1 3126 3209 jbase = jbase + inc2 … … 3129 3212 ELSE 3130 3213 3131 z = 1.0_wp/REAL(n,KIND=wp) 3132 zsin60 = z*sin60 3133 DO l = 1, la 3134 i = ibase 3135 j = jbase 3136 DO ipl = 1, SIZE(a,1) 3137 a11 = (a(ipl,ic+i)+a(ipl,if+i)) + (a(ipl,ib+i)+a(ipl,ie+i)) 3138 c(ipl,ja+j) = z*((a(ipl,ia+i)+a(ipl,id+i))+a11) 3139 c(ipl,jc+j) = z*((a(ipl,ia+i)+a(ipl,id+i))-0.5_wp*a11) 3140 d(ipl,jc+j) = zsin60*((a(ipl,ic+i)+a(ipl,if+i))-(a(ipl,ib+i)+a(ipl,ie+i))) 3141 a11 = (a(ipl,ic+i)-a(ipl,if+i)) + (a(ipl,ie+i)-a(ipl,ib+i)) 3142 c(ipl,jb+j) = z*((a(ipl,ia+i)-a(ipl,id+i))-0.5_wp*a11) 3143 d(ipl,jb+j) = zsin60*((a(ipl,ie+i)-a(ipl,ib+i))-(a(ipl,ic+i)-a(ipl,if+i))) 3144 c(ipl,jd+j) = z*((a(ipl,ia+i)-a(ipl,id+i))+a11) 3214 z = 1.0_wp / REAL( n, KIND = wp ) 3215 zsin60 = z * sin60 3216 DO l = 1, la 3217 i = ibase 3218 j = jbase 3219 DO ipl = 1, SIZE( a, 1 ) 3220 a11 = ( a (ipl,ic+i) + a(ipl,if+i) ) + ( a(ipl,ib+i) + a(ipl,ie+i) ) 3221 c(ipl,ja+j) = z * ( ( a(ipl,ia+i) + a(ipl,id+i) ) + a11 ) 3222 c(ipl,jc+j) = z * ( ( a(ipl,ia+i) + a(ipl,id+i) ) - 0.5_wp * a11 ) 3223 d(ipl,jc+j) = zsin60 * ( ( a(ipl,ic+i) + a(ipl,if+i) ) & 3224 - ( a(ipl,ib+i) + a(ipl,ie+i) ) ) 3225 a11 = ( a(ipl,ic+i) - a(ipl,if+i) ) + ( a(ipl,ie+i) - a(ipl,ib+i) ) 3226 c(ipl,jb+j) = z * ( ( a(ipl,ia+i) - a(ipl,id+i) ) - 0.5_wp * a11 ) 3227 d(ipl,jb+j) = zsin60 * ( ( a(ipl,ie+i) - a(ipl,ib+i) ) & 3228 - ( a(ipl,ic+i) - a(ipl,if+i) ) ) 3229 c(ipl,jd+j) = z * ( ( a(ipl,ia+i) - a(ipl,id+i) ) + a11 ) 3145 3230 ENDDO 3146 3231 ibase = ibase + inc1 … … 3169 3254 ja = 1 3170 3255 jb = ja + la * inc2 3171 jc = jb + 2 *m*inc23172 jd = jc + 2 *m*inc23173 je = jd + 2 *m*inc23174 z = 1.0_wp / REAL( n, KIND =wp )3256 jc = jb + 2 * m * inc2 3257 jd = jc + 2 * m * inc2 3258 je = jd + 2 * m * inc2 3259 z = 1.0_wp / REAL( n, KIND = wp ) 3175 3260 zsin45 = z * SQRT( 0.5_wp ) 3176 3261 … … 3178 3263 i = ibase 3179 3264 j = jbase 3180 c(:,ja+j) = z*(((a(:,ia+i)+a(:,ie+i))+(a(:,ic+i)+a(:,ig+i)))+((a(:,id+i)+ a(:,ih+i))+(a(:,ib+i)+a(:,if+i)))) 3181 c(:,je+j) = z*(((a(:,ia+i)+a(:,ie+i))+(a(:,ic+i)+a(:,ig+i)))-((a(:,id+i)+ a(:,ih+i))+(a(:,ib+i)+a(:,if+i)))) 3182 c(:,jc+j) = z*((a(:,ia+i)+a(:,ie+i))-(a(:,ic+i)+a(:,ig+i))) 3183 d(:,jc+j) = z*((a(:,id+i)+a(:,ih+i))-(a(:,ib+i)+a(:,if+i))) 3184 c(:,jb+j) = z*(a(:,ia+i)-a(:,ie+i)) + zsin45*((a(:,ih+i)-a(:,id+i))-(a(:,if+i)-a(:,ib+i))) 3185 c(:,jd+j) = z*(a(:,ia+i)-a(:,ie+i)) - zsin45*((a(:,ih+i)-a(:,id+i))-(a(:,if+i)-a(:,ib+i))) 3186 d(:,jb+j) = zsin45*((a(:,ih+i)-a(:,id+i))+(a(:,if+i)-a(:,ib+i))) + z*(a(:,ig+i)-a(:,ic+i)) 3187 d(:,jd+j) = zsin45*((a(:,ih+i)-a(:,id+i))+(a(:,if+i)-a(:,ib+i))) - z*(a(:,ig+i)-a(:,ic+i)) 3265 c(:,ja+j) = z * ( ( ( a(:,ia+i) + a(:,ie+i) ) + ( a(:,ic+i) + a(:,ig+i) ) ) & 3266 + ( ( a(:,id+i) + a(:,ih+i) ) + ( a(:,ib+i) + a(:,if+i) ) ) ) 3267 c(:,je+j) = z * ( ( ( a(:,ia+i) + a(:,ie+i) ) + ( a(:,ic+i) + a(:,ig+i) ) ) & 3268 - ( ( a(:,id+i) + a(:,ih+i) ) + ( a(:,ib+i) + a(:,if+i) ) ) ) 3269 c(:,jc+j) = z * ( ( a(:,ia+i) + a(:,ie+i) ) - ( a(:,ic+i) + a(:,ig+i) ) ) 3270 d(:,jc+j) = z * ( ( a(:,id+i) + a(:,ih+i) ) - ( a(:,ib+i) + a(:,if+i) ) ) 3271 c(:,jb+j) = z * ( a(:,ia+i) - a(:,ie+i) ) + zsin45 * ( ( a(:,ih+i) - a(:,id+i) ) & 3272 - ( a(:,if+i) - a(:,ib+i) ) ) 3273 c(:,jd+j) = z * ( a(:,ia+i) - a(:,ie+i) ) - zsin45 * ( ( a(:,ih+i) - a(:,id+i) ) & 3274 - ( a(:,if+i) - a(:,ib+i) ) ) 3275 d(:,jb+j) = zsin45 * ( ( a(:,ih+i) - a(:,id+i) ) + ( a(:,if+i) - a(:,ib+i) ) ) & 3276 + z * ( a(:,ig+i) - a(:,ic+i) ) 3277 d(:,jd+j) = zsin45 * ( ( a(:,ih+i) - a(:,id+i) ) + ( a(:,if+i) - a(:,ib+i) ) ) & 3278 - z * ( a(:,ig+i) - a(:,ic+i) ) 3188 3279 ibase = ibase + inc1 3189 3280 jbase = jbase + inc2 … … 3194 3285 END SUBROUTINE qpassm_vec 3195 3286 3196 !------------------------------------------------------------------------------ !3287 !--------------------------------------------------------------------------------------------------! 3197 3288 ! Description: 3198 3289 ! ------------ 3199 3290 !> Same as qpassm, but for backward fft 3200 !------------------------------------------------------------------------------ !3291 !--------------------------------------------------------------------------------------------------! 3201 3292 SUBROUTINE rpassm_vec(a, b, c, d, trigs, inc1, inc2, n, ifac, la, ierr ) 3202 3293 … … 3205 3296 IMPLICIT NONE 3206 3297 3207 INTEGER(iwp) :: ierr !< 3208 INTEGER(iwp) :: ifac !< 3209 INTEGER(iwp) :: inc1 !< 3210 INTEGER(iwp) :: inc2 !< 3211 INTEGER(iwp) :: la !< 3212 INTEGER(iwp) :: n !< 3298 3299 INTEGER(iwp) :: i !< 3300 INTEGER(iwp) :: ia !< 3301 INTEGER(iwp) :: ib !< 3302 INTEGER(iwp) :: ibase !< 3303 INTEGER(iwp) :: ic !< 3304 INTEGER(iwp) :: id !< 3305 INTEGER(iwp) :: ie !< 3306 INTEGER(iwp) :: ierr !< 3307 INTEGER(iwp) :: ifac !< 3308 INTEGER(iwp) :: if !< 3309 INTEGER(iwp) :: igo !< 3310 INTEGER(iwp) :: iink !< 3311 INTEGER(iwp) :: inc1 !< 3312 INTEGER(iwp) :: inc2 !< 3313 INTEGER(iwp) :: ipl !< loop index parallel loop 3314 INTEGER(iwp) :: j !< 3315 INTEGER(iwp) :: ja !< 3316 INTEGER(iwp) :: jb !< 3317 INTEGER(iwp) :: jbase !< 3318 INTEGER(iwp) :: jc !< 3319 INTEGER(iwp) :: jd !< 3320 INTEGER(iwp) :: je !< 3321 INTEGER(iwp) :: jf !< 3322 INTEGER(iwp) :: jg !< 3323 INTEGER(iwp) :: jh !< 3324 INTEGER(iwp) :: jink !< 3325 INTEGER(iwp) :: jump !< 3326 INTEGER(iwp) :: k !< 3327 INTEGER(iwp) :: kb !< 3328 INTEGER(iwp) :: kc !< 3329 INTEGER(iwp) :: kd !< 3330 INTEGER(iwp) :: ke !< 3331 INTEGER(iwp) :: kf !< 3332 INTEGER(iwp) :: kstop !< 3333 INTEGER(iwp) :: l !< 3334 INTEGER(iwp) :: la !< 3335 INTEGER(iwp) :: m !< 3336 INTEGER(iwp) :: n !< 3213 3337 ! 3214 3338 !-- Arrays are dimensioned with n 3215 REAL(wp),DIMENSION(:,:) :: a !< 3216 REAL(wp),DIMENSION(:,:) :: b !< 3217 REAL(wp),DIMENSION(:,:) :: c !< 3218 REAL(wp),DIMENSION(:,:) :: d !< 3219 REAL(wp),DIMENSION(:),INTENT(IN) :: trigs !< 3220 3221 REAL(wp) :: c1 !< 3222 REAL(wp) :: c2 !< 3223 REAL(wp) :: c3 !< 3224 REAL(wp) :: c4 !< 3225 REAL(wp) :: c5 !< 3226 REAL(wp) :: qqrt5 !< 3227 REAL(wp) :: qrt5 !< 3228 REAL(wp) :: s1 !< 3229 REAL(wp) :: s2 !< 3230 REAL(wp) :: s3 !< 3231 REAL(wp) :: s4 !< 3232 REAL(wp) :: s5 !< 3233 REAL(wp) :: sin36 !< 3234 REAL(wp) :: sin45 !< 3235 REAL(wp) :: sin60 !< 3236 REAL(wp) :: sin72 !< 3237 REAL(wp) :: ssin36 !< 3238 REAL(wp) :: ssin45 !< 3239 REAL(wp) :: ssin60 !< 3240 REAL(wp) :: ssin72 !< 3241 3242 INTEGER(iwp) :: i !< 3243 INTEGER(iwp) :: ia !< 3244 INTEGER(iwp) :: ib !< 3245 INTEGER(iwp) :: ibase !< 3246 INTEGER(iwp) :: ic !< 3247 INTEGER(iwp) :: id !< 3248 INTEGER(iwp) :: ie !< 3249 INTEGER(iwp) :: if !< 3250 INTEGER(iwp) :: igo !< 3251 INTEGER(iwp) :: iink !< 3252 INTEGER(iwp) :: ipl !< loop index parallel loop 3253 INTEGER(iwp) :: j !< 3254 INTEGER(iwp) :: ja !< 3255 INTEGER(iwp) :: jb !< 3256 INTEGER(iwp) :: jbase !< 3257 INTEGER(iwp) :: jc !< 3258 INTEGER(iwp) :: jd !< 3259 INTEGER(iwp) :: je !< 3260 INTEGER(iwp) :: jf !< 3261 INTEGER(iwp) :: jg !< 3262 INTEGER(iwp) :: jh !< 3263 INTEGER(iwp) :: jink !< 3264 INTEGER(iwp) :: jump !< 3265 INTEGER(iwp) :: k !< 3266 INTEGER(iwp) :: kb !< 3267 INTEGER(iwp) :: kc !< 3268 INTEGER(iwp) :: kd !< 3269 INTEGER(iwp) :: ke !< 3270 INTEGER(iwp) :: kf !< 3271 INTEGER(iwp) :: kstop !< 3272 INTEGER(iwp) :: l !< 3273 INTEGER(iwp) :: m !< 3274 3275 REAL(wp) :: a10 !< 3276 REAL(wp) :: a11 !< 3277 REAL(wp) :: a20 !< 3278 REAL(wp) :: a21 !< 3279 REAL(wp) :: b10 !< 3280 REAL(wp) :: b11 !< 3281 REAL(wp) :: b20 !< 3282 REAL(wp) :: b21 !< 3339 REAL(wp),DIMENSION(:,:) :: a !< 3340 REAL(wp),DIMENSION(:,:) :: b !< 3341 REAL(wp),DIMENSION(:,:) :: c !< 3342 REAL(wp),DIMENSION(:,:) :: d !< 3343 REAL(wp),DIMENSION(:),INTENT(IN) :: trigs !< 3344 3345 REAL(wp) :: a10 !< 3346 REAL(wp) :: a11 !< 3347 REAL(wp) :: a20 !< 3348 REAL(wp) :: a21 !< 3349 REAL(wp) :: b10 !< 3350 REAL(wp) :: b11 !< 3351 REAL(wp) :: b20 !< 3352 REAL(wp) :: b21 !< 3353 REAL(wp) :: c1 !< 3354 REAL(wp) :: c2 !< 3355 REAL(wp) :: c3 !< 3356 REAL(wp) :: c4 !< 3357 REAL(wp) :: c5 !< 3358 REAL(wp) :: qqrt5 !< 3359 REAL(wp) :: qrt5 !< 3360 REAL(wp) :: s1 !< 3361 REAL(wp) :: s2 !< 3362 REAL(wp) :: s3 !< 3363 REAL(wp) :: s4 !< 3364 REAL(wp) :: s5 !< 3365 REAL(wp) :: sin36 !< 3366 REAL(wp) :: sin45 !< 3367 REAL(wp) :: sin60 !< 3368 REAL(wp) :: sin72 !< 3369 REAL(wp) :: ssin36 !< 3370 REAL(wp) :: ssin45 !< 3371 REAL(wp) :: ssin60 !< 3372 REAL(wp) :: ssin72 !< 3373 3283 3374 3284 3375 … … 3292 3383 iink = la * inc1 3293 3384 jink = la * inc2 3294 jump = ( ifac-1) * jink3295 kstop = ( n-ifac) / (2*ifac)3385 jump = ( ifac - 1 ) * jink 3386 kstop = ( n - ifac ) / ( 2 * ifac ) 3296 3387 3297 3388 ibase = 0 … … 3311 3402 3312 3403 ia = 1 3313 ib = ia + ( 2*m-la) * inc13404 ib = ia + ( 2 * m - la ) * inc1 3314 3405 ja = 1 3315 3406 jb = ja + jink … … 3345 3436 c(:,ja+j) = a(:,ia+i) + a(:,ib+i) 3346 3437 d(:,ja+j) = b(:,ia+i) - b(:,ib+i) 3347 c(:,jb+j) = c1 *(a(:,ia+i)-a(:,ib+i)) - s1*(b(:,ia+i)+b(:,ib+i))3348 d(:,jb+j) = s1 *(a(:,ia+i)-a(:,ib+i)) + c1*(b(:,ia+i)+b(:,ib+i))3438 c(:,jb+j) = c1 * ( a(:,ia+i) - a(:,ib+i) ) - s1 * ( b(:,ia+i) + b(:,ib+i) ) 3439 d(:,jb+j) = s1 * ( a(:,ia+i) - a(:,ib+i) ) + c1 * ( b(:,ia+i) + b(:,ib+i) ) 3349 3440 ibase = ibase + inc1 3350 3441 jbase = jbase + inc2 … … 3364 3455 j = jbase 3365 3456 c(:,ja+j) = a(:,ia+i) 3366 c(:,jb+j) = - b(:,ia+i)3457 c(:,jb+j) = - b(:,ia+i) 3367 3458 ibase = ibase + inc1 3368 3459 jbase = jbase + inc2 … … 3374 3465 i = ibase 3375 3466 j = jbase 3376 c(:,ja+j) = 2.0_wp *(a(:,ia+i)+a(:,ib+i))3377 c(:,jb+j) = 2.0_wp *(a(:,ia+i)-a(:,ib+i))3467 c(:,ja+j) = 2.0_wp * ( a(:,ia+i) + a(:,ib+i) ) 3468 c(:,jb+j) = 2.0_wp * ( a(:,ia+i) - a(:,ib+i) ) 3378 3469 ibase = ibase + inc1 3379 3470 jbase = jbase + inc2 … … 3387 3478 3388 3479 ia = 1 3389 ib = ia + ( 2*m-la) * inc13480 ib = ia + ( 2 * m - la ) * inc1 3390 3481 ic = ib 3391 3482 ja = 1 … … 3399 3490 j = jbase 3400 3491 c(:,ja+j) = a(:,ia+i) + a(:,ib+i) 3401 c(:,jb+j) = ( a(:,ia+i)-0.5_wp*a(:,ib+i)) - (sin60*(b(:,ib+i)))3402 c(:,jc+j) = ( a(:,ia+i)-0.5_wp*a(:,ib+i)) + (sin60*(b(:,ib+i)))3492 c(:,jb+j) = ( a(:,ia+i) - 0.5_wp * a(:,ib+i) ) - ( sin60 * ( b(:,ib+i) ) ) 3493 c(:,jc+j) = ( a(:,ia+i) - 0.5_wp * a(:,ib+i) ) + ( sin60 * ( b(:,ib+i) ) ) 3403 3494 ibase = ibase + inc1 3404 3495 jbase = jbase + inc2 … … 3424 3515 i = ibase 3425 3516 j = jbase 3426 c(:,ja+j) = a(:,ia+i) + (a(:,ib+i)+a(:,ic+i)) 3427 d(:,ja+j) = b(:,ia+i) + (b(:,ib+i)-b(:,ic+i)) 3428 c(:,jb+j) = c1*((a(:,ia+i)-0.5_wp*(a(:,ib+i)+a(:,ic+i)))-(sin60*(b(:,ib+i)+ b(:,ic+i)))) & 3429 - s1*((b(:,ia+i)-0.5_wp*(b(:,ib+i)-b(:,ic+i)))+(sin60*(a(:,ib+i)- a(:,ic+i)))) 3430 d(:,jb+j) = s1*((a(:,ia+i)-0.5_wp*(a(:,ib+i)+a(:,ic+i)))-(sin60*(b(:,ib+i)+ b(:,ic+i)))) & 3431 + c1*((b(:,ia+i)-0.5_wp*(b(:,ib+i)-b(:,ic+i)))+(sin60*(a(:,ib+i)- a(:,ic+i)))) 3432 c(:,jc+j) = c2*((a(:,ia+i)-0.5_wp*(a(:,ib+i)+a(:,ic+i)))+(sin60*(b(:,ib+i)+ b(:,ic+i)))) & 3433 - s2*((b(:,ia+i)-0.5_wp*(b(:,ib+i)-b(:,ic+i)))-(sin60*(a(:,ib+i)- a(:,ic+i)))) 3434 d(:,jc+j) = s2*((a(:,ia+i)-0.5_wp*(a(:,ib+i)+a(:,ic+i)))+(sin60*(b(:,ib+i)+ b(:,ic+i)))) & 3435 + c2*((b(:,ia+i)-0.5_wp*(b(:,ib+i)-b(:,ic+i)))-(sin60*(a(:,ib+i)- a(:,ic+i)))) 3517 c(:,ja+j) = a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) ) 3518 d(:,ja+j) = b(:,ia+i) + ( b(:,ib+i) - b(:,ic+i) ) 3519 c(:,jb+j) = c1 * ( ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) ) & 3520 - ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) ) - s1 * ( ( b(:,ia+i) & 3521 - 0.5_wp * ( b(:,ib+i) - b(:,ic+i) ) ) + ( sin60 * ( a(:,ib+i) & 3522 - a(:,ic+i) ) ) ) 3523 d(:,jb+j) = s1 * ( ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) ) & 3524 - ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) ) + c1 * ( ( b(:,ia+i) & 3525 - 0.5_wp * ( b(:,ib+i) - b(:,ic+i) ) ) + ( sin60 * ( a(:,ib+i) & 3526 - a(:,ic+i) ) ) ) 3527 c(:,jc+j) = c2 * ( ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) ) & 3528 + ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) ) - s2 * ( ( b(:,ia+i) & 3529 - 0.5_wp * ( b(:,ib+i) - b(:,ic+i) ) ) - ( sin60 * ( a(:,ib+i) & 3530 - a(:,ic+i) ) ) ) 3531 d(:,jc+j) = s2 * ( ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) ) & 3532 + ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) ) + c2 * ( ( b(:,ia+i) & 3533 - 0.5_wp * ( b(:,ib+i) - b(:,ic+i) ) ) - ( sin60 * ( a(:,ib+i) & 3534 - a(:,ic+i) ) ) ) 3436 3535 ibase = ibase + inc1 3437 3536 jbase = jbase + inc2 … … 3452 3551 j = jbase 3453 3552 c(:,ja+j) = a(:,ia+i) + a(:,ib+i) 3454 c(:,jb+j) = ( 0.5_wp*a(:,ia+i)-a(:,ib+i)) - (sin60*b(:,ia+i))3455 c(:,jc+j) = - (0.5_wp*a(:,ia+i)-a(:,ib+i)) - (sin60*b(:,ia+i))3553 c(:,jb+j) = ( 0.5_wp * a(:,ia+i) - a(:,ib+i) ) - ( sin60 * b(:,ia+i) ) 3554 c(:,jc+j) = - ( 0.5_wp * a(:,ia+i) - a(:,ib+i) ) - ( sin60 * b(:,ia+i) ) 3456 3555 ibase = ibase + inc1 3457 3556 jbase = jbase + inc2 … … 3464 3563 i = ibase 3465 3564 j = jbase 3466 c(:,ja+j) = 2.0_wp *(a(:,ia+i)+a(:,ib+i))3467 c(:,jb+j) = ( 2.0_wp*a(:,ia+i)-a(:,ib+i)) - (ssin60*b(:,ib+i))3468 c(:,jc+j) = ( 2.0_wp*a(:,ia+i)-a(:,ib+i)) + (ssin60*b(:,ib+i))3565 c(:,ja+j) = 2.0_wp * ( a(:,ia+i) + a(:,ib+i) ) 3566 c(:,jb+j) = ( 2.0_wp * a(:,ia+i) - a(:,ib+i) ) - ( ssin60 * b(:,ib+i) ) 3567 c(:,jc+j) = ( 2.0_wp * a(:,ia+i) - a(:,ib+i) ) + ( ssin60 * b(:,ib+i) ) 3469 3568 ibase = ibase + inc1 3470 3569 jbase = jbase + inc2 … … 3478 3577 3479 3578 ia = 1 3480 ib = ia + ( 2*m-la) * inc13481 ic = ib + 2 *m*inc13579 ib = ia + ( 2 * m - la ) * inc1 3580 ic = ib + 2 * m * inc1 3482 3581 id = ib 3483 3582 ja = 1 … … 3491 3590 i = ibase 3492 3591 j = jbase 3493 c(:,ja+j) = ( a(:,ia+i)+a(:,ic+i)) + a(:,ib+i)3494 c(:,jb+j) = ( a(:,ia+i)-a(:,ic+i)) - b(:,ib+i)3495 c(:,jc+j) = ( a(:,ia+i)+a(:,ic+i)) - a(:,ib+i)3496 c(:,jd+j) = ( a(:,ia+i)-a(:,ic+i)) + b(:,ib+i)3592 c(:,ja+j) = ( a(:,ia+i) + a(:,ic+i) ) + a(:,ib+i) 3593 c(:,jb+j) = ( a(:,ia+i) - a(:,ic+i) ) - b(:,ib+i) 3594 c(:,jc+j) = ( a(:,ia+i) + a(:,ic+i) ) - a(:,ib+i) 3595 c(:,jd+j) = ( a(:,ia+i) - a(:,ic+i) ) + b(:,ib+i) 3497 3596 ibase = ibase + inc1 3498 3597 jbase = jbase + inc2 … … 3522 3621 i = ibase 3523 3622 j = jbase 3524 c(:,ja+j) = ( a(:,ia+i)+a(:,ic+i)) + (a(:,ib+i)+a(:,id+i))3525 d(:,ja+j) = ( b(:,ia+i)-b(:,ic+i)) + (b(:,ib+i)-b(:,id+i))3526 c(:,jc+j) = c2 *((a(:,ia+i)+a(:,ic+i))-(a(:,ib+i)+a(:,id+i))) - s2*((b(:,ia+i)-b(:,ic+i))&3527 -(b(:,ib+i)-b(:,id+i)))3528 d(:,jc+j) = s2 *((a(:,ia+i)+a(:,ic+i))-(a(:,ib+i)+a(:,id+i))) + c2*((b(:,ia+i)-b(:,ic+i))&3529 -(b(:,ib+i)-b(:,id+i)))3530 c(:,jb+j) = c1 *((a(:,ia+i)-a(:,ic+i))-(b(:,ib+i)+b(:,id+i))) - s1*((b(:,ia+i)+b(:,ic+i))&3531 +(a(:,ib+i)-a(:,id+i)))3532 d(:,jb+j) = s1 *((a(:,ia+i)-a(:,ic+i))-(b(:,ib+i)+b(:,id+i))) + c1*((b(:,ia+i)+b(:,ic+i))&3533 +(a(:,ib+i)-a(:,id+i)))3534 c(:,jd+j) = c3 *((a(:,ia+i)-a(:,ic+i))+(b(:,ib+i)+b(:,id+i))) - s3*((b(:,ia+i)+b(:,ic+i))&3535 -(a(:,ib+i)-a(:,id+i)))3536 d(:,jd+j) = s3 *((a(:,ia+i)-a(:,ic+i))+(b(:,ib+i)+b(:,id+i))) + c3*((b(:,ia+i)+b(:,ic+i))&3537 -(a(:,ib+i)-a(:,id+i)))3623 c(:,ja+j) = ( a(:,ia+i) + a(:,ic+i) ) + ( a(:,ib+i) + a(:,id+i) ) 3624 d(:,ja+j) = ( b(:,ia+i) - b(:,ic+i) ) + ( b(:,ib+i) - b(:,id+i) ) 3625 c(:,jc+j) = c2 * ( ( a(:,ia+i) + a(:,ic+i) ) - ( a(:,ib+i) + a(:,id+i) ) ) & 3626 - s2 * ( ( b(:,ia+i) - b(:,ic+i) ) - ( b(:,ib+i) - b(:,id+i) ) ) 3627 d(:,jc+j) = s2 * ( ( a(:,ia+i) + a(:,ic+i) ) - ( a(:,ib+i) + a(:,id+i) ) ) & 3628 + c2 * ( ( b(:,ia+i) - b(:,ic+i) ) - ( b(:,ib+i) - b(:,id+i) ) ) 3629 c(:,jb+j) = c1 * ( ( a(:,ia+i) - a(:,ic+i) ) - ( b(:,ib+i) + b(:,id+i) ) ) & 3630 - s1 * ( ( b(:,ia+i) + b(:,ic+i) ) + ( a(:,ib+i) - a(:,id+i) ) ) 3631 d(:,jb+j) = s1 * ( ( a(:,ia+i) - a(:,ic+i) ) - ( b(:,ib+i) + b(:,id+i) ) ) & 3632 + c1 * ( ( b(:,ia+i) + b(:,ic+i) ) + ( a(:,ib+i) - a(:,id+i) ) ) 3633 c(:,jd+j) = c3 * ( ( a(:,ia+i) - a(:,ic+i) ) + ( b(:,ib+i) + b(:,id+i) ) ) & 3634 - s3 * ( ( b(:,ia+i) + b(:,ic+i) ) - ( a(:,ib+i) - a(:,id+i) ) ) 3635 d(:,jd+j) = s3 * ( ( a(:,ia+i) - a(:,ic+i) ) + ( b(:,ib+i) + b(:,id+i) ) ) & 3636 + c3 * ( ( b(:,ia+i) + b(:,ic+i) ) - ( a(:,ib+i) - a(:,id+i) ) ) 3538 3637 ibase = ibase + inc1 3539 3638 jbase = jbase + inc2 … … 3556 3655 j = jbase 3557 3656 c(:,ja+j) = a(:,ia+i) + a(:,ib+i) 3558 c(:,jb+j) = sin45 *((a(:,ia+i)-a(:,ib+i))-(b(:,ia+i)+b(:,ib+i)))3657 c(:,jb+j) = sin45 * ( ( a(:,ia+i) - a(:,ib+i) ) - ( b(:,ia+i) + b(:,ib+i) ) ) 3559 3658 c(:,jc+j) = b(:,ib+i) - b(:,ia+i) 3560 c(:,jd+j) = - sin45*((a(:,ia+i)-a(:,ib+i))+(b(:,ia+i)+b(:,ib+i)))3659 c(:,jd+j) = - sin45 * ( ( a(:,ia+i) - a(:,ib+i) ) + ( b(:,ia+i) + b(:,ib+i) ) ) 3561 3660 ibase = ibase + inc1 3562 3661 jbase = jbase + inc2 … … 3568 3667 i = ibase 3569 3668 j = jbase 3570 c(:,ja+j) = 2.0_wp *((a(:,ia+i)+a(:,ic+i))+a(:,ib+i))3571 c(:,jb+j) = 2.0_wp *((a(:,ia+i)-a(:,ic+i))-b(:,ib+i))3572 c(:,jc+j) = 2.0_wp *((a(:,ia+i)+a(:,ic+i))-a(:,ib+i))3573 c(:,jd+j) = 2.0_wp *((a(:,ia+i)-a(:,ic+i))+b(:,ib+i))3669 c(:,ja+j) = 2.0_wp * ( ( a(:,ia+i) + a(:,ic+i) ) + a(:,ib+i) ) 3670 c(:,jb+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ic+i) ) - b(:,ib+i) ) 3671 c(:,jc+j) = 2.0_wp * ( ( a(:,ia+i) + a(:,ic+i) ) - a(:,ib+i) ) 3672 c(:,jd+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ic+i) ) + b(:,ib+i) ) 3574 3673 ibase = ibase + inc1 3575 3674 jbase = jbase + inc2 … … 3583 3682 3584 3683 ia = 1 3585 ib = ia + ( 2*m-la) * inc13586 ic = ib + 2 *m*inc13684 ib = ia + ( 2 * m - la ) * inc1 3685 ic = ib + 2 * m * inc1 3587 3686 id = ic 3588 3687 ie = ib … … 3598 3697 i = ibase 3599 3698 j = jbase 3600 c(:,ja+j) = a(:,ia+i) + ( a(:,ib+i)+a(:,ic+i))3601 c(:,jb+j) = ( (a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))+qrt5*(a(:,ib+i)-a(:,ic+i))) -&3602 (sin72*b(:,ib+i)+sin36*b(:,ic+i))3603 c(:,jc+j) = ( (a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))-qrt5*(a(:,ib+i)-a(:,ic+i))) -&3604 (sin36*b(:,ib+i)-sin72*b(:,ic+i))3605 c(:,jd+j) = ( (a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))-qrt5*(a(:,ib+i)-a(:,ic+i))) +&3606 (sin36*b(:,ib+i)-sin72*b(:,ic+i))3607 c(:,je+j) = ( (a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))+qrt5*(a(:,ib+i)-a(:,ic+i))) +&3608 (sin72*b(:,ib+i)+sin36*b(:,ic+i))3699 c(:,ja+j) = a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) ) 3700 c(:,jb+j) = ( ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) + qrt5 & 3701 * ( a(:,ib+i) - a(:,ic+i) ) ) - ( sin72 * b(:,ib+i) + sin36 * b(:,ic+i) ) 3702 c(:,jc+j) = ( ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) - qrt5 & 3703 * ( a(:,ib+i) - a(:,ic+i) ) ) - ( sin36 * b(:,ib+i) - sin72 * b(:,ic+i) ) 3704 c(:,jd+j) = ( ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) - qrt5 & 3705 * ( a(:,ib+i) - a(:,ic+i) ) ) + ( sin36 * b(:,ib+i) - sin72 * b(:,ic+i) ) 3706 c(:,je+j) = ( ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) + qrt5 & 3707 * ( a(:,ib+i) - a(:,ic+i) ) ) + ( sin72 * b(:,ib+i) + sin36 * b(:,ic+i) ) 3609 3708 ibase = ibase + inc1 3610 3709 jbase = jbase + inc2 … … 3639 3738 j = jbase 3640 3739 !DIR$ IVDEP 3641 DO ipl = 1, SIZE(a,1) 3642 a10 = (a(ipl,ia+i)-0.25_wp*((a(ipl,ib+i)+a(ipl,ie+i))+(a(ipl,ic+i)+a(ipl,id+i)))) + & 3643 qrt5*((a(ipl,ib+i)+a(ipl,ie+i))-(a(ipl,ic+i)+a(ipl,id+i))) 3644 a20 = (a(ipl,ia+i)-0.25_wp*((a(ipl,ib+i)+a(ipl,ie+i))+(a(ipl,ic+i)+a(ipl,id+i)))) - & 3645 qrt5*((a(ipl,ib+i)+a(ipl,ie+i))-(a(ipl,ic+i)+a(ipl,id+i))) 3646 b10 = (b(ipl,ia+i)-0.25_wp*((b(ipl,ib+i)-b(ipl,ie+i))+(b(ipl,ic+i)-b(ipl,id+i)))) + & 3647 qrt5*((b(ipl,ib+i)-b(ipl,ie+i))-(b(ipl,ic+i)-b(ipl,id+i))) 3648 b20 = (b(ipl,ia+i)-0.25_wp*((b(ipl,ib+i)-b(ipl,ie+i))+(b(ipl,ic+i)-b(ipl,id+i)))) - & 3649 qrt5*((b(ipl,ib+i)-b(ipl,ie+i))-(b(ipl,ic+i)-b(ipl,id+i))) 3650 a11 = sin72*(b(ipl,ib+i)+b(ipl,ie+i)) + sin36*(b(ipl,ic+i)+b(ipl,id+i)) 3651 a21 = sin36*(b(ipl,ib+i)+b(ipl,ie+i)) - sin72*(b(ipl,ic+i)+b(ipl,id+i)) 3652 b11 = sin72*(a(ipl,ib+i)-a(ipl,ie+i)) + sin36*(a(ipl,ic+i)-a(ipl,id+i)) 3653 b21 = sin36*(a(ipl,ib+i)-a(ipl,ie+i)) - sin72*(a(ipl,ic+i)-a(ipl,id+i)) 3654 c(ipl,ja+j) = a(ipl,ia+i) + ((a(ipl,ib+i)+a(ipl,ie+i))+(a(ipl,ic+i)+a(ipl,id+i))) 3655 d(ipl,ja+j) = b(ipl,ia+i) + ((b(ipl,ib+i)-b(ipl,ie+i))+(b(ipl,ic+i)-b(ipl,id+i))) 3656 c(ipl,jb+j) = c1*(a10 -a11 ) - s1*(b10 +b11 ) 3657 d(ipl,jb+j) = s1*(a10 -a11 ) + c1*(b10 +b11 ) 3658 c(ipl,je+j) = c4*(a10 +a11 ) - s4*(b10 -b11 ) 3659 d(ipl,je+j) = s4*(a10 +a11 ) + c4*(b10 -b11 ) 3660 c(ipl,jc+j) = c2*(a20 -a21 ) - s2*(b20 +b21 ) 3661 d(ipl,jc+j) = s2*(a20 -a21 ) + c2*(b20 +b21 ) 3662 c(ipl,jd+j) = c3*(a20 +a21 ) - s3*(b20 -b21 ) 3663 d(ipl,jd+j) = s3*(a20 +a21 ) + c3*(b20 -b21 ) 3740 DO ipl = 1, SIZE( a, 1 ) 3741 a10 = ( a (ipl,ia+i) - 0.25_wp * ( ( a(ipl,ib+i) + a(ipl,ie+i) ) & 3742 + ( a(ipl,ic+i) + a(ipl,id+i) ) ) ) + qrt5 * ( ( a(ipl,ib+i) & 3743 + a(ipl,ie+i) ) - ( a(ipl,ic+i) + a(ipl,id+i) ) ) 3744 a20 = ( a(ipl,ia+i) - 0.25_wp * ( ( a(ipl,ib+i) + a(ipl,ie+i) ) & 3745 + ( a(ipl,ic+i) + a(ipl,id+i) ) ) ) - qrt5 * ( ( a(ipl,ib+i) & 3746 + a(ipl,ie+i) ) - ( a(ipl,ic+i) + a(ipl,id+i) ) ) 3747 b10 = ( b(ipl,ia+i) - 0.25_wp * ( ( b(ipl,ib+i) - b(ipl,ie+i) ) & 3748 + ( b(ipl,ic+i) - b(ipl,id+i) ) ) ) + qrt5 * ( ( b(ipl,ib+i) & 3749 - b(ipl,ie+i) ) - ( b(ipl,ic+i) - b(ipl,id+i) ) ) 3750 b20 = ( b(ipl,ia+i) - 0.25_wp * ( ( b(ipl,ib+i) - b(ipl,ie+i) ) & 3751 + ( b(ipl,ic+i) - b(ipl,id+i) ) ) ) - qrt5 * ( ( b(ipl,ib+i) & 3752 - b(ipl,ie+i) ) - ( b(ipl,ic+i) - b(ipl,id+i) ) ) 3753 a11 = sin72 * ( b(ipl,ib+i) + b(ipl,ie+i) ) + sin36 & 3754 * ( b(ipl,ic+i) + b(ipl,id+i) ) 3755 a21 = sin36 * ( b(ipl,ib+i) + b(ipl,ie+i) ) - sin72 & 3756 * ( b(ipl,ic+i) + b(ipl,id+i) ) 3757 b11 = sin72 * ( a(ipl,ib+i) - a(ipl,ie+i) ) + sin36 & 3758 * ( a(ipl,ic+i) - a(ipl,id+i) ) 3759 b21 = sin36 * ( a(ipl,ib+i) - a(ipl,ie+i) ) - sin72 & 3760 * ( a(ipl,ic+i) - a(ipl,id+i) ) 3761 c(ipl,ja+j) = a(ipl,ia+i) + ( ( a(ipl,ib+i) + a(ipl,ie+i) ) & 3762 + ( a(ipl,ic+i) + a(ipl,id+i) ) ) 3763 d(ipl,ja+j) = b(ipl,ia+i) + ( ( b(ipl,ib+i) - b(ipl,ie+i) ) & 3764 + ( b(ipl,ic+i) - b(ipl,id+i) ) ) 3765 c(ipl,jb+j) = c1 * ( a10 - a11 ) - s1 * ( b10 + b11 ) 3766 d(ipl,jb+j) = s1 * ( a10 - a11 ) + c1 * ( b10 + b11 ) 3767 c(ipl,je+j) = c4 * ( a10 + a11 ) - s4 * ( b10 - b11 ) 3768 d(ipl,je+j) = s4 * ( a10 + a11 ) + c4 * ( b10 - b11 ) 3769 c(ipl,jc+j) = c2 * ( a20 - a21 ) - s2 * ( b20 + b21 ) 3770 d(ipl,jc+j) = s2 * ( a20 - a21 ) + c2 * ( b20 + b21 ) 3771 c(ipl,jd+j) = c3 * ( a20 + a21 ) - s3 * ( b20 - b21 ) 3772 d(ipl,jd+j) = s3 * ( a20 + a21 ) + c3 * ( b20 - b21 ) 3664 3773 ENDDO 3665 3774 ibase = ibase + inc1 … … 3682 3791 i = ibase 3683 3792 j = jbase 3684 c(:,ja+j) = ( a(:,ia+i)+a(:,ib+i)) + a(:,ic+i)3685 c(:,jb+j) = ( qrt5*(a(:,ia+i)-a(:,ib+i))+(0.25_wp*(a(:,ia+i)+a(:,ib+i))-a(:,ic+i))) -&3686 (sin36*b(:,ia+i)+sin72*b(:,ib+i))3687 c(:,je+j) = - (qrt5*(a(:,ia+i)-a(:,ib+i))+(0.25_wp*(a(:,ia+i)+a(:,ib+i))-a(:,ic+i))) -&3688 (sin36*b(:,ia+i)+sin72*b(:,ib+i))3689 c(:,jc+j) = ( qrt5*(a(:,ia+i)-a(:,ib+i))-(0.25_wp*(a(:,ia+i)+a(:,ib+i))-a(:,ic+i))) -&3690 (sin72*b(:,ia+i)-sin36*b(:,ib+i))3691 c(:,jd+j) = - (qrt5*(a(:,ia+i)-a(:,ib+i))-(0.25_wp*(a(:,ia+i)+a(:,ib+i))-a(:,ic+i))) -&3692 (sin72*b(:,ia+i)-sin36*b(:,ib+i))3793 c(:,ja+j) = ( a(:,ia+i) + a(:,ib+i) ) + a(:,ic+i) 3794 c(:,jb+j) = ( qrt5 * ( a(:,ia+i) - a(:,ib+i) ) + ( 0.25_wp * ( a(:,ia+i) & 3795 + a(:,ib+i) ) - a(:,ic+i) ) ) - ( sin36 * b(:,ia+i) + sin72 * b(:,ib+i) ) 3796 c(:,je+j) = - ( qrt5 * ( a(:,ia+i) - a(:,ib+i) ) + ( 0.25_wp * ( a(:,ia+i) & 3797 + a(:,ib+i) ) - a(:,ic+i) ) ) - ( sin36 * b(:,ia+i) + sin72 * b(:,ib+i) ) 3798 c(:,jc+j) = ( qrt5 * ( a(:,ia+i) - a(:,ib+i) ) - ( 0.25_wp * ( a(:,ia+i) & 3799 + a(:,ib+i) ) - a(:,ic+i) ) ) - ( sin72 * b(:,ia+i) - sin36 * b(:,ib+i) ) 3800 c(:,jd+j) = - ( qrt5 * ( a(:,ia+i) - a(:,ib+i) ) - ( 0.25_wp * ( a(:,ia+i) & 3801 + a(:,ib+i) ) - a(:,ic+i) ) ) - ( sin72 * b(:,ia+i) - sin36 * b(:,ib+i) ) 3693 3802 ibase = ibase + inc1 3694 3803 jbase = jbase + inc2 … … 3703 3812 i = ibase 3704 3813 j = jbase 3705 c(:,ja+j) = 2.0_wp*(a(:,ia+i)+(a(:,ib+i)+a(:,ic+i))) 3706 c(:,jb+j) = (2.0_wp*(a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))+qqrt5*(a(:,ib+i)-a(:,ic+i))) & 3707 - (ssin72*b(:,ib+i)+ssin36*b(:,ic+i)) 3708 c(:,jc+j) = (2.0_wp*(a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))-qqrt5*(a(:,ib+i)-a(:,ic+i))) & 3709 - (ssin36*b(:,ib+i)-ssin72*b(:,ic+i)) 3710 c(:,jd+j) = (2.0_wp*(a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))-qqrt5*(a(:,ib+i)-a(:,ic+i))) & 3711 + (ssin36*b(:,ib+i)-ssin72*b(:,ic+i)) 3712 c(:,je+j) = (2.0_wp*(a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))+qqrt5*(a(:,ib+i)-a(:,ic+i))) & 3713 + (ssin72*b(:,ib+i)+ssin36*b(:,ic+i)) 3814 c(:,ja+j) = 2.0_wp * ( a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) ) ) 3815 c(:,jb+j) = ( 2.0_wp * ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) & 3816 + qqrt5 * ( a(:,ib+i) - a(:,ic+i) ) ) - ( ssin72 * b(:,ib+i) & 3817 + ssin36 * b(:,ic+i) ) 3818 c(:,jc+j) = ( 2.0_wp * ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) & 3819 - qqrt5 * ( a(:,ib+i) - a(:,ic+i) ) ) - ( ssin36 * b(:,ib+i) & 3820 - ssin72 * b(:,ic+i) ) 3821 c(:,jd+j) = ( 2.0_wp * ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) & 3822 - qqrt5 * ( a(:,ib+i) - a(:,ic+i) ) ) + ( ssin36 * b(:,ib+i) & 3823 - ssin72 * b(:,ic+i) ) 3824 c(:,je+j) = ( 2.0_wp * ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) & 3825 + qqrt5 * ( a(:,ib+i) - a(:,ic+i) ) ) + ( ssin72 * b(:,ib+i) & 3826 + ssin36 * b(:,ic+i) ) 3714 3827 ibase = ibase + inc1 3715 3828 jbase = jbase + inc2 … … 3723 3836 3724 3837 ia = 1 3725 ib = ia + ( 2*m-la) * inc13726 ic = ib + 2 *m*inc13727 id = ic + 2 *m*inc13838 ib = ia + ( 2 * m - la ) * inc1 3839 ic = ib + 2 * m * inc1 3840 id = ic + 2 * m * inc1 3728 3841 ie = ic 3729 3842 if = ib … … 3740 3853 i = ibase 3741 3854 j = jbase 3742 c(:,ja+j) = (a(:,ia+i)+a(:,id+i)) + (a(:,ib+i)+a(:,ic+i)) 3743 c(:,jd+j) = (a(:,ia+i)-a(:,id+i)) - (a(:,ib+i)-a(:,ic+i)) 3744 c(:,jb+j) = ((a(:,ia+i)-a(:,id+i))+0.5_wp*(a(:,ib+i)-a(:,ic+i))) - (sin60*(b(:,ib+i)+b(:,ic+i))) 3745 c(:,jf+j) = ((a(:,ia+i)-a(:,id+i))+0.5_wp*(a(:,ib+i)-a(:,ic+i))) + (sin60*(b(:,ib+i)+b(:,ic+i))) 3746 c(:,jc+j) = ((a(:,ia+i)+a(:,id+i))-0.5_wp*(a(:,ib+i)+a(:,ic+i))) - (sin60*(b(:,ib+i)-b(:,ic+i))) 3747 c(:,je+j) = ((a(:,ia+i)+a(:,id+i))-0.5_wp*(a(:,ib+i)+a(:,ic+i))) + (sin60*(b(:,ib+i)-b(:,ic+i))) 3855 c(:,ja+j) = ( a(:,ia+i) + a(:,id+i) ) + ( a(:,ib+i) + a(:,ic+i) ) 3856 c(:,jd+j) = ( a(:,ia+i) - a(:,id+i) ) - ( a(:,ib+i) - a(:,ic+i) ) 3857 c(:,jb+j) = ( ( a(:,ia+i) - a(:,id+i) ) + 0.5_wp * ( a(:,ib+i) - a(:,ic+i) ) ) & 3858 - ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) 3859 c(:,jf+j) = ( ( a(:,ia+i) - a(:,id+i) ) + 0.5_wp * ( a(:,ib+i) - a(:,ic+i) ) ) & 3860 + ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) 3861 c(:,jc+j) = ( ( a(:,ia+i) + a(:,id+i) ) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) ) & 3862 - ( sin60 * ( b(:,ib+i) - b(:,ic+i) ) ) 3863 c(:,je+j) = ( ( a(:,ia+i) + a(:,id+i) ) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) ) & 3864 + ( sin60 * ( b(:,ib+i) - b(:,ic+i) ) ) 3748 3865 ibase = ibase + inc1 3749 3866 jbase = jbase + inc2 … … 3781 3898 i = ibase 3782 3899 j = jbase 3783 DO ipl = 1, SIZE(a,1) 3784 a11 = (a(ipl,ie+i)+a(ipl,ib+i)) + (a(ipl,ic+i)+a(ipl,if+i)) 3785 a20 = (a(ipl,ia+i)+a(ipl,id+i)) - 0.5_wp*a11 3786 a21 = sin60*((a(ipl,ie+i)+a(ipl,ib+i))-(a(ipl,ic+i)+a(ipl,if+i))) 3787 b11 = (b(ipl,ib+i)-b(ipl,ie+i)) + (b(ipl,ic+i)-b(ipl,if+i)) 3788 b20 = (b(ipl,ia+i)-b(ipl,id+i)) - 0.5_wp*b11 3789 b21 = sin60*((b(ipl,ib+i)-b(ipl,ie+i))-(b(ipl,ic+i)-b(ipl,if+i))) 3790 3791 c(ipl,ja+j) = (a(ipl,ia+i)+a(ipl,id+i)) + a11 3792 d(ipl,ja+j) = (b(ipl,ia+i)-b(ipl,id+i)) + b11 3793 c(ipl,jc+j) = c2*(a20 -b21 ) - s2*(b20 +a21 ) 3794 d(ipl,jc+j) = s2*(a20 -b21 ) + c2*(b20 +a21 ) 3795 c(ipl,je+j) = c4*(a20 +b21 ) - s4*(b20 -a21 ) 3796 d(ipl,je+j) = s4*(a20 +b21 ) + c4*(b20 -a21 ) 3797 3798 a11 = (a(ipl,ie+i)-a(ipl,ib+i)) + (a(ipl,ic+i)-a(ipl,if+i)) 3799 b11 = (b(ipl,ie+i)+b(ipl,ib+i)) - (b(ipl,ic+i)+b(ipl,if+i)) 3800 a20 = (a(ipl,ia+i)-a(ipl,id+i)) - 0.5_wp*a11 3801 a21 = sin60*((a(ipl,ie+i)-a(ipl,ib+i))-(a(ipl,ic+i)-a(ipl,if+i))) 3802 b20 = (b(ipl,ia+i)+b(ipl,id+i)) + 0.5_wp*b11 3803 b21 = sin60*((b(ipl,ie+i)+b(ipl,ib+i))+(b(ipl,ic+i)+b(ipl,if+i))) 3804 3805 c(ipl,jd+j) = c3*((a(ipl,ia+i)-a(ipl,id+i))+a11 ) - s3*((b(ipl,ia+i)+b(ipl,id+i))-b11 ) 3806 d(ipl,jd+j) = s3*((a(ipl,ia+i)-a(ipl,id+i))+a11 ) + c3*((b(ipl,ia+i)+b(ipl,id+i))-b11 ) 3807 c(ipl,jb+j) = c1*(a20 -b21 ) - s1*(b20 -a21 ) 3808 d(ipl,jb+j) = s1*(a20 -b21 ) + c1*(b20 -a21 ) 3809 c(ipl,jf+j) = c5*(a20 +b21 ) - s5*(b20 +a21 ) 3810 d(ipl,jf+j) = s5*(a20 +b21 ) + c5*(b20 +a21 ) 3900 DO ipl = 1, SIZE( a, 1) 3901 a11 = ( a(ipl,ie+i) + a(ipl,ib+i) ) + ( a(ipl,ic+i) + a(ipl,if+i) ) 3902 a20 = ( a(ipl,ia+i) + a(ipl,id+i) ) - 0.5_wp * a11 3903 a21 = sin60 * ( ( a(ipl,ie+i) + a(ipl,ib+i) ) & 3904 - ( a(ipl,ic+i) + a(ipl,if+i) ) ) 3905 b11 = ( b(ipl,ib+i) - b(ipl,ie+i) ) + ( b(ipl,ic+i) - b(ipl,if+i) ) 3906 b20 = ( b(ipl,ia+i) - b(ipl,id+i) ) - 0.5_wp * b11 3907 b21 = sin60 * ( ( b(ipl,ib+i) - b(ipl,ie+i) ) & 3908 - ( b(ipl,ic+i) - b(ipl,if+i) ) ) 3909 3910 c(ipl,ja+j) = ( a(ipl,ia+i) + a(ipl,id+i) ) + a11 3911 d(ipl,ja+j) = ( b(ipl,ia+i) - b(ipl,id+i) ) + b11 3912 c(ipl,jc+j) = c2 * ( a20 - b21 ) - s2 * ( b20 + a21 ) 3913 d(ipl,jc+j) = s2 * ( a20 - b21 ) + c2 * ( b20 + a21 ) 3914 c(ipl,je+j) = c4 * ( a20 + b21 ) - s4 * ( b20 - a21 ) 3915 d(ipl,je+j) = s4 * ( a20 + b21 ) + c4 * ( b20 - a21 ) 3916 3917 a11 = ( a(ipl,ie+i) - a(ipl,ib+i) ) + ( a(ipl,ic+i) - a(ipl,if+i) ) 3918 b11 = ( b(ipl,ie+i) + b(ipl,ib+i) ) - ( b(ipl,ic+i) + b(ipl,if+i) ) 3919 a20 = ( a(ipl,ia+i) - a(ipl,id+i) ) - 0.5_wp * a11 3920 a21 = sin60 * ( ( a(ipl,ie+i) - a(ipl,ib+i) ) - ( a(ipl,ic+i) & 3921 - a(ipl,if+i) ) ) 3922 b20 = ( b(ipl,ia+i) + b(ipl,id+i) ) + 0.5_wp * b11 3923 b21 = sin60 * ( ( b(ipl,ie+i) + b(ipl,ib+i) ) + ( b(ipl,ic+i) & 3924 + b(ipl,if+i) ) ) 3925 3926 c(ipl,jd+j) = c3 * ( ( a(ipl,ia+i) - a(ipl,id+i) ) + a11 ) & 3927 - s3 * ( ( b(ipl,ia+i) + b(ipl,id+i) ) - b11 ) 3928 d(ipl,jd+j) = s3 * ( ( a(ipl,ia+i) - a(ipl,id+i) ) + a11 ) & 3929 + c3 * ( ( b(ipl,ia+i) + b(ipl,id+i) ) - b11 ) 3930 c(ipl,jb+j) = c1 * ( a20 - b21 ) - s1 * ( b20 - a21 ) 3931 d(ipl,jb+j) = s1 * ( a20 - b21 ) + c1 * ( b20 - a21 ) 3932 c(ipl,jf+j) = c5 * ( a20 + b21 ) - s5 * ( b20 + a21 ) 3933 d(ipl,jf+j) = s5 * ( a20 + b21 ) + c5 * ( b20 + a21 ) 3811 3934 3812 3935 ENDDO … … 3831 3954 i = ibase 3832 3955 j = jbase 3833 DO ipl = 1, SIZE(a,1) 3834 c(ipl,ja+j) = a(ipl,ib+i) + (a(ipl,ia+i)+a(ipl,ic+i)) 3835 c(ipl,jd+j) = b(ipl,ib+i) - (b(ipl,ia+i)+b(ipl,ic+i)) 3836 c(ipl,jb+j) = (sin60*(a(ipl,ia+i)-a(ipl,ic+i))) - (0.5_wp*(b(ipl,ia+i)+b(ipl,ic+i))+b(ipl,ib+i)) 3837 c(ipl,jf+j) = -(sin60*(a(ipl,ia+i)-a(ipl,ic+i))) - (0.5_wp*(b(ipl,ia+i)+b(ipl,ic+i))+b(ipl,ib+i)) 3838 c(ipl,jc+j) = sin60*(b(ipl,ic+i)-b(ipl,ia+i)) + (0.5_wp*(a(ipl,ia+i)+a(ipl,ic+i))-a(ipl,ib+i)) 3839 c(ipl,je+j) = sin60*(b(ipl,ic+i)-b(ipl,ia+i)) - (0.5_wp*(a(ipl,ia+i)+a(ipl,ic+i))-a(ipl,ib+i)) 3956 DO ipl = 1, SIZE( a ,1 ) 3957 c(ipl,ja+j) = a(ipl,ib+i) + ( a(ipl,ia+i) + a(ipl,ic+i) ) 3958 c(ipl,jd+j) = b(ipl,ib+i) - ( b(ipl,ia+i) + b(ipl,ic+i) ) 3959 c(ipl,jb+j) = ( sin60 * ( a(ipl,ia+i) - a(ipl,ic+i) ) ) & 3960 - ( 0.5_wp * ( b(ipl,ia+i) + b(ipl,ic+i) ) + b(ipl,ib+i) ) 3961 c(ipl,jf+j) = - ( sin60 * ( a(ipl,ia+i) - a(ipl,ic+i) ) ) & 3962 - ( 0.5_wp * ( b(ipl,ia+i) + b(ipl,ic+i) ) + b(ipl,ib+i) ) 3963 c(ipl,jc+j) = sin60 * ( b(ipl,ic+i) - b(ipl,ia+i) ) & 3964 + ( 0.5_wp * ( a(ipl,ia+i) + a(ipl,ic+i) ) - a(ipl,ib+i) ) 3965 c(ipl,je+j) = sin60 * ( b(ipl,ic+i) - b(ipl,ia+i) ) & 3966 - ( 0.5_wp * ( a(ipl,ia+i) + a(ipl,ic+i) ) - a(ipl,ib+i) ) 3840 3967 ENDDO 3841 3968 ibase = ibase + inc1 … … 3849 3976 i = ibase 3850 3977 j = jbase 3851 c(:,ja+j) = (2.0_wp*(a(:,ia+i)+a(:,id+i))) + (2.0_wp*(a(:,ib+i)+a(:,ic+i))) 3852 c(:,jd+j) = (2.0_wp*(a(:,ia+i)-a(:,id+i))) - (2.0_wp*(a(:,ib+i)-a(:,ic+i))) 3853 c(:,jb+j) = (2.0_wp*(a(:,ia+i)-a(:,id+i))+(a(:,ib+i)-a(:,ic+i))) - (ssin60*(b(:,ib+i)+b(:,ic+i))) 3854 c(:,jf+j) = (2.0_wp*(a(:,ia+i)-a(:,id+i))+(a(:,ib+i)-a(:,ic+i))) + (ssin60*(b(:,ib+i)+b(:,ic+i))) 3855 c(:,jc+j) = (2.0_wp*(a(:,ia+i)+a(:,id+i))-(a(:,ib+i)+a(:,ic+i))) - (ssin60*(b(:,ib+i)-b(:,ic+i))) 3856 c(:,je+j) = (2.0_wp*(a(:,ia+i)+a(:,id+i))-(a(:,ib+i)+a(:,ic+i))) + (ssin60*(b(:,ib+i)-b(:,ic+i))) 3978 c(:,ja+j) = ( 2.0_wp * ( a(:,ia+i) + a(:,id+i) ) ) & 3979 + ( 2.0_wp * ( a(:,ib+i) + a(:,ic+i) ) ) 3980 c(:,jd+j) = ( 2.0_wp * ( a(:,ia+i) - a(:,id+i) ) ) & 3981 - ( 2.0_wp * ( a(:,ib+i) - a(:,ic+i) ) ) 3982 c(:,jb+j) = ( 2.0_wp * ( a(:,ia+i) - a(:,id+i) ) + ( a(:,ib+i) - a(:,ic+i) ) ) & 3983 - ( ssin60 * ( b(:,ib+i) + b(:,ic+i) ) ) 3984 c(:,jf+j) = ( 2.0_wp * ( a(:,ia+i) - a(:,id+i) ) + ( a(:,ib+i) - a(:,ic+i) ) ) & 3985 + ( ssin60 * ( b(:,ib+i) + b(:,ic+i) ) ) 3986 c(:,jc+j) = ( 2.0_wp * ( a(:,ia+i) + a(:,id+i) ) - ( a(:,ib+i) + a(:,ic+i) ) ) & 3987 - ( ssin60 * ( b(:,ib+i) - b(:,ic+i) ) ) 3988 c(:,je+j) = ( 2.0_wp * ( a(:,ia+i) + a(:,id+i) ) - ( a(:,ib+i) + a(:,ic+i) ) ) & 3989 + ( ssin60 * ( b(:,ib+i) - b(:,ic+i) ) ) 3857 3990 ibase = ibase + inc1 3858 3991 jbase = jbase + inc2 … … 3871 4004 3872 4005 ia = 1 3873 ib = ia + la *inc13874 ic = ib + 2 *la*inc13875 id = ic + 2 *la*inc13876 ie = id + 2 *la*inc14006 ib = ia + la * inc1 4007 ic = ib + 2 * la * inc1 4008 id = ic + 2 * la * inc1 4009 ie = id + 2 * la * inc1 3877 4010 ja = 1 3878 4011 jb = ja + jink … … 3888 4021 i = ibase 3889 4022 j = jbase 3890 c(:,ja+j) = 2.0_wp*(((a(:,ia+i)+a(:,ie+i))+a(:,ic+i))+(a(:,ib+i)+a(:,id+i))) 3891 c(:,je+j) = 2.0_wp*(((a(:,ia+i)+a(:,ie+i))+a(:,ic+i))-(a(:,ib+i)+a(:,id+i))) 3892 c(:,jc+j) = 2.0_wp*(((a(:,ia+i)+a(:,ie+i))-a(:,ic+i))-(b(:,ib+i)-b(:,id+i))) 3893 c(:,jg+j) = 2.0_wp*(((a(:,ia+i)+a(:,ie+i))-a(:,ic+i))+(b(:,ib+i)-b(:,id+i))) 3894 c(:,jb+j) = 2.0_wp*((a(:,ia+i)-a(:,ie+i))-b(:,ic+i)) + ssin45*((a(:,ib+i)-a(:,id+i))-(b(:,ib+i)+b(:,id+i))) 3895 c(:,jf+j) = 2.0_wp*((a(:,ia+i)-a(:,ie+i))-b(:,ic+i)) - ssin45*((a(:,ib+i)-a(:,id+i))-(b(:,ib+i)+b(:,id+i))) 3896 c(:,jd+j) = 2.0_wp*((a(:,ia+i)-a(:,ie+i))+b(:,ic+i)) - ssin45*((a(:,ib+i)-a(:,id+i))+(b(:,ib+i)+b(:,id+i))) 3897 c(:,jh+j) = 2.0_wp*((a(:,ia+i)-a(:,ie+i))+b(:,ic+i)) + ssin45*((a(:,ib+i)-a(:,id+i))+(b(:,ib+i)+b(:,id+i))) 4023 c(:,ja+j) = 2.0_wp * ( ( ( a(:,ia+i) + a(:,ie+i) ) + a(:,ic+i) ) & 4024 + ( a(:,ib+i) + a(:,id+i) ) ) 4025 c(:,je+j) = 2.0_wp * ( ( ( a(:,ia+i) + a(:,ie+i) ) + a(:,ic+i) ) & 4026 - ( a(:,ib+i) + a(:,id+i) ) ) 4027 c(:,jc+j) = 2.0_wp * ( ( ( a(:,ia+i) + a(:,ie+i) ) - a(:,ic+i) ) & 4028 - ( b(:,ib+i) - b(:,id+i) ) ) 4029 c(:,jg+j) = 2.0_wp * ( ( ( a(:,ia+i) + a(:,ie+i) ) - a(:,ic+i) ) & 4030 + ( b(:,ib+i) - b(:,id+i) ) ) 4031 c(:,jb+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ie+i) ) - b(:,ic+i) ) + ssin45 & 4032 * ( (a(:,ib+i) - a(:,id+i) ) - ( b(:,ib+i) + b(:,id+i) ) ) 4033 c(:,jf+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ie+i) ) - b(:,ic+i) ) - ssin45 & 4034 * ( (a(:,ib+i) - a(:,id+i) ) - ( b(:,ib+i) + b(:,id+i) ) ) 4035 c(:,jd+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ie+i) ) + b(:,ic+i) ) - ssin45 & 4036 * ( (a(:,ib+i) - a(:,id+i) ) + ( b(:,ib+i) + b(:,id+i) ) ) 4037 c(:,jh+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ie+i) ) + b(:,ic+i) ) + ssin45 & 4038 * ( (a(:,ib+i) - a(:,id+i) ) + ( b(:,ib+i) + b(:,id+i) ) ) 3898 4039 ibase = ibase + inc1 3899 4040 jbase = jbase + inc2
Note: See TracChangeset
for help on using the changeset viewer.
