Changeset 4651 for palm/trunk/SOURCE/fft_xy_mod.f90
- Timestamp:
- Aug 27, 2020 7:17:45 AM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/fft_xy_mod.f90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/fft_xy_mod.f90
r4646 r4651 24 24 ! ----------------- 25 25 ! $Id$ 26 ! preprocessor branch for ibm removed 27 ! 28 ! 4646 2020-08-24 16:02:40Z raasch 26 29 ! file re-formatted to follow the PALM coding standard 27 30 ! … … 99 102 REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: f_vec_x 100 103 101 #if defined( __ibm ) 102 INTEGER(iwp), PARAMETER :: nau1 = 20000 !< 103 INTEGER(iwp), PARAMETER :: nau2 = 22000 !< 104 ! 105 !-- The following working arrays contain tables and have to be "save" and shared in OpenMP sense 106 REAL(wp), DIMENSION(nau1), SAVE :: aux1 !< 107 REAL(wp), DIMENSION(nau1), SAVE :: auy1 !< 108 REAL(wp), DIMENSION(nau1), SAVE :: aux3 !< 109 REAL(wp), DIMENSION(nau1), SAVE :: auy3 !< 110 111 #elif defined( __nec_fft ) 104 #if defined( __nec_fft ) 112 105 INTEGER(iwp), SAVE :: nz1 !< 113 106 … … 189 182 !-- The following temporary working arrays have to be on stack or private 190 183 !-- in OpenMP sense 191 #if defined( __ibm ) 192 REAL(wp), DIMENSION(nau2) :: aux2 !< 193 REAL(wp), DIMENSION(nau2) :: auy2 !< 194 REAL(wp), DIMENSION(nau2) :: aux4 !< 195 REAL(wp), DIMENSION(nau2) :: auy4 !< 196 REAL(wp), DIMENSION(0:nx+2) :: workx !< 197 REAL(wp), DIMENSION(0:ny+2) :: worky !< 198 #elif defined( __nec_fft ) 184 #if defined( __nec_fft ) 199 185 REAL(wp), DIMENSION(0:nx+3,nz+1) :: work_x !< 200 186 REAL(wp), DIMENSION(0:ny+3,nz+1) :: work_y !< … … 228 214 sqr_dnx = SQRT( dnx ) 229 215 sqr_dny = SQRT( dny ) 230 #if defined( __ibm ) 231 ! 232 !-- Initialize tables for fft along x 233 CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, aux2, nau2 ) 234 CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, aux4, nau2 ) 235 ! 236 !-- Initialize tables for fft along y 237 CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, auy2, nau2 ) 238 CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, auy4, nau2 ) 239 #elif defined( __nec_fft ) 216 217 #if defined( __nec_fft ) 240 218 message_string = 'fft method "' // TRIM( fft_method) // '" currently does not work on NEC' 241 219 CALL message( 'fft_init', 'PA0187', 1, 2, 0, 6, 0 ) … … 350 328 REAL(wp), DIMENSION(nys_x:nyn_x,nzb_x:nzt_x,0:nx), OPTIONAL :: ar_inv !< 351 329 352 #if defined( __ibm ) 353 REAL(wp), DIMENSION(nau2) :: aux2 !< 354 REAL(wp), DIMENSION(nau2) :: aux4 !< 355 #elif defined( __nec_fft ) 330 #if defined( __nec_fft ) 356 331 REAL(wp), DIMENSION(6*(nx+1)) :: work2 !< 357 332 #elif defined( __cuda_fft ) … … 643 618 ELSEIF ( fft_method == 'system-specific' ) THEN 644 619 645 #if defined( __ibm ) 646 IF ( forward_fft ) THEN 647 648 !$OMP PARALLEL PRIVATE ( work, i, j, k ) 649 !$OMP DO 650 DO k = nzb_x, nzt_x 651 DO j = nys_x, nyn_x 652 653 CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, aux2, nau2 ) 654 655 DO i = 0, (nx+1)/2 656 ar(i,j,k) = work(2*i) 657 ENDDO 658 DO i = 1, (nx+1)/2 - 1 659 ar(nx+1-i,j,k) = work(2*i+1) 660 ENDDO 661 662 ENDDO 663 ENDDO 664 !$OMP END PARALLEL 665 666 ELSE 667 668 !$OMP PARALLEL PRIVATE ( work, i, j, k ) 669 !$OMP DO 670 DO k = nzb_x, nzt_x 671 DO j = nys_x, nyn_x 672 673 DO i = 0, (nx+1)/2 674 work(2*i) = ar(i,j,k) 675 ENDDO 676 DO i = 1, (nx+1)/2 - 1 677 work(2*i+1) = ar(nx+1-i,j,k) 678 ENDDO 679 work(1) = 0.0_wp 680 work(nx+2) = 0.0_wp 681 682 CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, aux4, nau2 ) 683 684 DO i = 0, nx 685 ar(i,j,k) = work(i) 686 ENDDO 687 688 ENDDO 689 ENDDO 690 !$OMP END PARALLEL 691 692 ENDIF 693 694 #elif defined( __nec_fft ) 620 #if defined( __nec_fft ) 695 621 696 622 IF ( forward_fft ) THEN … … 822 748 COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: cwork !< 823 749 824 #if defined( __ibm ) 825 REAL(wp), DIMENSION(nau2) :: aux2 !< 826 REAL(wp), DIMENSION(nau2) :: aux4 !< 827 #elif defined( __nec_fft ) 750 #if defined( __nec_fft ) 828 751 REAL(wp), DIMENSION(6*(nx+1)) :: work2 !< 829 752 #endif … … 940 863 ELSEIF ( fft_method == 'system-specific' ) THEN 941 864 942 #if defined( __ibm ) 943 IF ( forward_fft ) THEN 944 945 CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, aux2, nau2 ) 946 947 DO i = 0, (nx+1)/2 948 ar(i) = work(2*i) 949 ENDDO 950 DO i = 1, (nx+1)/2 - 1 951 ar(nx+1-i) = work(2*i+1) 952 ENDDO 953 954 ELSE 955 956 DO i = 0, (nx+1)/2 957 work(2*i) = ar(i) 958 ENDDO 959 DO i = 1, (nx+1)/2 - 1 960 work(2*i+1) = ar(nx+1-i) 961 ENDDO 962 work(1) = 0.0_wp 963 work(nx+2) = 0.0_wp 964 965 CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, aux4, nau2 ) 966 967 DO i = 0, nx 968 ar(i) = work(i) 969 ENDDO 970 971 ENDIF 972 #elif defined( __nec_fft ) 865 #if defined( __nec_fft ) 973 866 IF ( forward_fft ) THEN 974 867 … … 1061 954 COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: cwork !< 1062 955 1063 #if defined( __ibm ) 1064 REAL(wp), DIMENSION(nau2) :: auy2 !< 1065 REAL(wp), DIMENSION(nau2) :: auy4 !< 1066 #elif defined( __nec_fft ) 956 #if defined( __nec_fft ) 1067 957 REAL(wp), DIMENSION(6*(ny+1)) :: work2 !< 1068 958 #elif defined( __cuda_fft ) … … 1347 1237 ELSEIF ( fft_method == 'system-specific' ) THEN 1348 1238 1349 #if defined( __ibm ) 1350 IF ( forward_fft) THEN 1351 1352 !$OMP PARALLEL PRIVATE ( work, i, j, k ) 1353 !$OMP DO 1354 DO k = nzb_y, nzt_y 1355 DO i = nxl_y_l, nxr_y_l 1356 1357 CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, auy2, nau2 ) 1358 1359 DO j = 0, (ny+1)/2 1360 ar_tr(j,i,k) = work(2*j) 1361 ENDDO 1362 DO j = 1, (ny+1)/2 - 1 1363 ar_tr(ny+1-j,i,k) = work(2*j+1) 1364 ENDDO 1365 1366 ENDDO 1367 ENDDO 1368 !$OMP END PARALLEL 1369 1370 ELSE 1371 1372 !$OMP PARALLEL PRIVATE ( work, i, j, k ) 1373 !$OMP DO 1374 DO k = nzb_y, nzt_y 1375 DO i = nxl_y_l, nxr_y_l 1376 1377 DO j = 0, (ny+1)/2 1378 work(2*j) = ar_tr(j,i,k) 1379 ENDDO 1380 DO j = 1, (ny+1)/2 - 1 1381 work(2*j+1) = ar_tr(ny+1-j,i,k) 1382 ENDDO 1383 work(1) = 0.0_wp 1384 work(ny+2) = 0.0_wp 1385 1386 CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, auy4, nau2 ) 1387 1388 DO j = 0, ny 1389 ar(j,i,k) = work(j) 1390 ENDDO 1391 1392 ENDDO 1393 ENDDO 1394 !$OMP END PARALLEL 1395 1396 ENDIF 1397 #elif defined( __nec_fft ) 1239 #if defined( __nec_fft ) 1398 1240 IF ( forward_fft ) THEN 1399 1241 … … 1523 1365 COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: cwork !< 1524 1366 1525 #if defined( __ibm ) 1526 REAL(wp), DIMENSION(nau2) :: auy2 !< 1527 REAL(wp), DIMENSION(nau2) :: auy4 !< 1528 #elif defined( __nec_fft ) 1367 #if defined( __nec_fft ) 1529 1368 REAL(wp), DIMENSION(6*(ny+1)) :: work2 !< 1530 1369 #endif … … 1643 1482 ELSEIF ( fft_method == 'system-specific' ) THEN 1644 1483 1645 #if defined( __ibm ) 1646 IF ( forward_fft ) THEN 1647 1648 CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, auy2, nau2 ) 1649 1650 DO j = 0, (ny+1)/2 1651 ar(j) = work(2*j) 1652 ENDDO 1653 DO j = 1, (ny+1)/2 - 1 1654 ar(ny+1-j) = work(2*j+1) 1655 ENDDO 1656 1657 ELSE 1658 1659 DO j = 0, (ny+1)/2 1660 work(2*j) = ar(j) 1661 ENDDO 1662 DO j = 1, (ny+1)/2 - 1 1663 work(2*j+1) = ar(ny+1-j) 1664 ENDDO 1665 work(1) = 0.0_wp 1666 work(ny+2) = 0.0_wp 1667 1668 CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, auy4, nau2 ) 1669 1670 DO j = 0, ny 1671 ar(j) = work(j) 1672 ENDDO 1673 1674 ENDIF 1675 #elif defined( __nec_fft ) 1484 #if defined( __nec_fft ) 1676 1485 IF ( forward_fft ) THEN 1677 1486
Note: See TracChangeset
for help on using the changeset viewer.
