Changeset 1216 for palm/trunk/SOURCE/fft_xy.f90
- Timestamp:
- Aug 26, 2013 9:31:42 AM (11 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/fft_xy.f90 (modified) (27 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/fft_xy.f90
r1211 r1216 286 286 287 287 288 SUBROUTINE fft_x( ar, direction )288 SUBROUTINE fft_x( ar, direction, ar_2d ) 289 289 290 290 !----------------------------------------------------------------------! … … 321 321 COMPLEX(dpk), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) :: ar_tmp 322 322 #endif 323 REAL, DIMENSION(0:nx,nys_x:nyn_x), OPTIONAL :: ar_2d 323 324 REAL, DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) :: ar 324 325 … … 454 455 CALL FFTW_EXECUTE_DFT_R2C( plan_xf, x_in, x_out ) 455 456 456 DO i = 0, (nx+1)/2 457 ar(i,j,k) = REAL( x_out(i) ) / ( nx+1 ) 458 ENDDO 459 DO i = 1, (nx+1)/2 - 1 460 ar(nx+1-i,j,k) = AIMAG( x_out(i) ) / ( nx+1 ) 461 ENDDO 457 IF ( PRESENT( ar_2d ) ) THEN 458 459 DO i = 0, (nx+1)/2 460 ar_2d(i,j) = REAL( x_out(i) ) / ( nx+1 ) 461 ENDDO 462 DO i = 1, (nx+1)/2 - 1 463 ar_2d(nx+1-i,j) = AIMAG( x_out(i) ) / ( nx+1 ) 464 ENDDO 465 466 ELSE 467 468 DO i = 0, (nx+1)/2 469 ar(i,j,k) = REAL( x_out(i) ) / ( nx+1 ) 470 ENDDO 471 DO i = 1, (nx+1)/2 - 1 472 ar(nx+1-i,j,k) = AIMAG( x_out(i) ) / ( nx+1 ) 473 ENDDO 474 475 ENDIF 462 476 463 477 ENDDO … … 465 479 !$OMP END PARALLEL 466 480 467 ELSE481 ELSE 468 482 !$OMP PARALLEL PRIVATE ( work, i, j, k ) 469 483 !$OMP DO … … 471 485 DO j = nys_x, nyn_x 472 486 473 x_out(0) = CMPLX( ar(0,j,k), 0.0 ) 474 DO i = 1, (nx+1)/2 - 1 475 x_out(i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) ) 476 ENDDO 477 x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 ) 487 IF ( PRESENT( ar_2d ) ) THEN 488 489 x_out(0) = CMPLX( ar_2d(0,j), 0.0 ) 490 DO i = 1, (nx+1)/2 - 1 491 x_out(i) = CMPLX( ar_2d(i,j), ar_2d(nx+1-i,j) ) 492 ENDDO 493 x_out((nx+1)/2) = CMPLX( ar_2d((nx+1)/2,j), 0.0 ) 494 495 ELSE 496 497 x_out(0) = CMPLX( ar(0,j,k), 0.0 ) 498 DO i = 1, (nx+1)/2 - 1 499 x_out(i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) ) 500 ENDDO 501 x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 ) 502 503 ENDIF 478 504 479 505 CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in) … … 484 510 !$OMP END PARALLEL 485 511 486 ENDIF512 ENDIF 487 513 #endif 488 514 … … 766 792 ENDIF 767 793 794 ELSEIF ( fft_method == 'fftw' ) THEN 795 796 #if defined( __fftw ) 797 IF ( forward_fft ) THEN 798 799 x_in(0:nx) = ar(0:nx) 800 CALL FFTW_EXECUTE_DFT_R2C( plan_xf, x_in, x_out ) 801 802 DO i = 0, (nx+1)/2 803 ar(i) = REAL( x_out(i) ) / ( nx+1 ) 804 ENDDO 805 DO i = 1, (nx+1)/2 - 1 806 ar(nx+1-i) = AIMAG( x_out(i) ) / ( nx+1 ) 807 ENDDO 808 809 ELSE 810 811 x_out(0) = CMPLX( ar(0), 0.0 ) 812 DO i = 1, (nx+1)/2 - 1 813 x_out(i) = CMPLX( ar(i), ar(nx+1-i) ) 814 ENDDO 815 x_out((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 ) 816 817 CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in) 818 ar(0:nx) = x_in(0:nx) 819 820 ENDIF 821 #endif 822 768 823 ELSEIF ( fft_method == 'system-specific' ) THEN 769 824 … … 843 898 END SUBROUTINE fft_x_1d 844 899 845 SUBROUTINE fft_y( ar, direction ) 900 SUBROUTINE fft_y( ar, direction, ar_tr, nxl_y_bound, nxr_y_bound, nxl_y_l, & 901 nxr_y_l ) 846 902 847 903 !----------------------------------------------------------------------! … … 853 909 ! fft_y uses internal algorithms (Singleton or Temperton) or ! 854 910 ! system-specific routines, if they are available ! 911 ! ! 912 ! direction: 'forward' or 'backward' ! 913 ! ar, ar_tr: 3D data arrays ! 914 ! forward: ar: before ar_tr: after transformation ! 915 ! backward: ar_tr: before ar: after transfosition ! 916 ! ! 917 ! In case of non-overlapping transposition/transformation: ! 918 ! nxl_y_bound = nxl_y_l = nxl_y ! 919 ! nxr_y_bound = nxr_y_l = nxr_y ! 920 ! ! 921 ! In case of overlapping transposition/transformation ! 922 ! - nxl_y_bound and nxr_y_bound have the original values of ! 923 ! nxl_y, nxr_y. ar_tr is dimensioned using these values. ! 924 ! - nxl_y_l = nxr_y_r. ar is dimensioned with these values, so that ! 925 ! transformation is carried out for a 2D-plane only. ! 855 926 !----------------------------------------------------------------------! 856 927 … … 864 935 CHARACTER (LEN=*) :: direction 865 936 INTEGER :: i, j, jshape(1), k 937 INTEGER :: nxl_y_bound, nxl_y_l, nxr_y_bound, nxr_y_l 866 938 867 939 LOGICAL :: forward_fft … … 878 950 COMPLEX(dpk), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) :: ar_tmp 879 951 #endif 880 REAL, DIMENSION(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) :: ar 952 REAL, DIMENSION(0:ny,nxl_y_l:nxr_y_l,nzb_y:nzt_y) :: ar 953 REAL, DIMENSION(0:ny,nxl_y_bound:nxr_y_bound,nzb_y:nzt_y) :: ar_tr 881 954 882 955 IF ( direction == 'forward' ) THEN … … 898 971 !$OMP DO 899 972 DO k = nzb_y, nzt_y 900 DO i = nxl_y , nxr_y973 DO i = nxl_y_l, nxr_y_l 901 974 902 975 DO j = 0, ny … … 908 981 909 982 DO j = 0, (ny+1)/2 910 ar (j,i,k) = REAL( cwork(j) )983 ar_tr(j,i,k) = REAL( cwork(j) ) 911 984 ENDDO 912 985 DO j = 1, (ny+1)/2 - 1 913 ar (ny+1-j,i,k) = -AIMAG( cwork(j) )986 ar_tr(ny+1-j,i,k) = -AIMAG( cwork(j) ) 914 987 ENDDO 915 988 … … 923 996 !$OMP DO 924 997 DO k = nzb_y, nzt_y 925 DO i = nxl_y , nxr_y926 927 cwork(0) = CMPLX( ar (0,i,k), 0.0 )998 DO i = nxl_y_l, nxr_y_l 999 1000 cwork(0) = CMPLX( ar_tr(0,i,k), 0.0 ) 928 1001 DO j = 1, (ny+1)/2 - 1 929 cwork(j) = CMPLX( ar (j,i,k), -ar(ny+1-j,i,k) )930 cwork(ny+1-j) = CMPLX( ar (j,i,k), ar(ny+1-j,i,k) )931 ENDDO 932 cwork((ny+1)/2) = CMPLX( ar ((ny+1)/2,i,k), 0.0 )1002 cwork(j) = CMPLX( ar_tr(j,i,k), -ar_tr(ny+1-j,i,k) ) 1003 cwork(ny+1-j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k) ) 1004 ENDDO 1005 cwork((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0 ) 933 1006 934 1007 jshape = SHAPE( cwork ) … … 957 1030 !$OMP DO 958 1031 DO k = nzb_y, nzt_y 959 DO i = nxl_y , nxr_y1032 DO i = nxl_y_l, nxr_y_l 960 1033 961 1034 work(0:ny) = ar(0:ny,i,k) … … 963 1036 964 1037 DO j = 0, (ny+1)/2 965 ar (j,i,k) = work(2*j)1038 ar_tr(j,i,k) = work(2*j) 966 1039 ENDDO 967 1040 DO j = 1, (ny+1)/2 - 1 968 ar (ny+1-j,i,k) = work(2*j+1)1041 ar_tr(ny+1-j,i,k) = work(2*j+1) 969 1042 ENDDO 970 1043 … … 978 1051 !$OMP DO 979 1052 DO k = nzb_y, nzt_y 980 DO i = nxl_y , nxr_y1053 DO i = nxl_y_l, nxr_y_l 981 1054 982 1055 DO j = 0, (ny+1)/2 983 work(2*j) = ar (j,i,k)1056 work(2*j) = ar_tr(j,i,k) 984 1057 ENDDO 985 1058 DO j = 1, (ny+1)/2 - 1 986 work(2*j+1) = ar (ny+1-j,i,k)1059 work(2*j+1) = ar_tr(ny+1-j,i,k) 987 1060 ENDDO 988 1061 work(1) = 0.0 … … 1006 1079 !$OMP DO 1007 1080 DO k = nzb_y, nzt_y 1008 DO i = nxl_y , nxr_y1081 DO i = nxl_y_l, nxr_y_l 1009 1082 1010 1083 y_in(0:ny) = ar(0:ny,i,k) … … 1012 1085 1013 1086 DO j = 0, (ny+1)/2 1014 ar (j,i,k) = REAL( y_out(j) ) /(ny+1)1087 ar_tr(j,i,k) = REAL( y_out(j) ) / (ny+1) 1015 1088 ENDDO 1016 1089 DO j = 1, (ny+1)/2 - 1 1017 ar (ny+1-j,i,k) = AIMAG( y_out(j) ) /(ny+1)1090 ar_tr(ny+1-j,i,k) = AIMAG( y_out(j) ) / (ny+1) 1018 1091 ENDDO 1019 1092 … … 1027 1100 !$OMP DO 1028 1101 DO k = nzb_y, nzt_y 1029 DO i = nxl_y , nxr_y1030 1031 y_out(0) = CMPLX( ar (0,i,k), 0.0 )1102 DO i = nxl_y_l, nxr_y_l 1103 1104 y_out(0) = CMPLX( ar_tr(0,i,k), 0.0 ) 1032 1105 DO j = 1, (ny+1)/2 - 1 1033 y_out(j) = CMPLX( ar (j,i,k), ar(ny+1-j,i,k) )1034 ENDDO 1035 y_out((ny+1)/2) = CMPLX( ar ((ny+1)/2,i,k), 0.0 )1106 y_out(j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k) ) 1107 ENDDO 1108 y_out((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0 ) 1036 1109 1037 1110 CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in ) … … 1053 1126 !$OMP DO 1054 1127 DO k = nzb_y, nzt_y 1055 DO i = nxl_y , nxr_y1128 DO i = nxl_y_l, nxr_y_l 1056 1129 1057 1130 CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, & … … 1059 1132 1060 1133 DO j = 0, (ny+1)/2 1061 ar (j,i,k) = work(2*j)1134 ar_tr(j,i,k) = work(2*j) 1062 1135 ENDDO 1063 1136 DO j = 1, (ny+1)/2 - 1 1064 ar (ny+1-j,i,k) = work(2*j+1)1137 ar_tr(ny+1-j,i,k) = work(2*j+1) 1065 1138 ENDDO 1066 1139 … … 1074 1147 !$OMP DO 1075 1148 DO k = nzb_y, nzt_y 1076 DO i = nxl_y , nxr_y1149 DO i = nxl_y_l, nxr_y_l 1077 1150 1078 1151 DO j = 0, (ny+1)/2 1079 work(2*j) = ar (j,i,k)1152 work(2*j) = ar_tr(j,i,k) 1080 1153 ENDDO 1081 1154 DO j = 1, (ny+1)/2 - 1 1082 work(2*j+1) = ar (ny+1-j,i,k)1155 work(2*j+1) = ar_tr(ny+1-j,i,k) 1083 1156 ENDDO 1084 1157 work(1) = 0.0 … … 1103 1176 !$OMP DO 1104 1177 DO k = nzb_y, nzt_y 1105 DO i = nxl_y , nxr_y1178 DO i = nxl_y_l, nxr_y_l 1106 1179 1107 1180 work(0:ny) = ar(0:ny,i,k) … … 1110 1183 1111 1184 DO j = 0, (ny+1)/2 1112 ar (j,i,k) = work(2*j)1185 ar_tr(j,i,k) = work(2*j) 1113 1186 ENDDO 1114 1187 DO j = 1, (ny+1)/2 - 1 1115 ar (ny+1-j,i,k) = work(2*j+1)1188 ar_tr(ny+1-j,i,k) = work(2*j+1) 1116 1189 ENDDO 1117 1190 … … 1125 1198 !$OMP DO 1126 1199 DO k = nzb_y, nzt_y 1127 DO i = nxl_y , nxr_y1200 DO i = nxl_y_l, nxr_y_l 1128 1201 1129 1202 DO j = 0, (ny+1)/2 1130 work(2*j) = ar (j,i,k)1203 work(2*j) = ar_tr(j,i,k) 1131 1204 ENDDO 1132 1205 DO j = 1, (ny+1)/2 - 1 1133 work(2*j+1) = ar (ny+1-j,i,k)1206 work(2*j+1) = ar_tr(ny+1-j,i,k) 1134 1207 ENDDO 1135 1208 work(1) = 0.0 … … 1322 1395 1323 1396 ENDIF 1397 1398 ELSEIF ( fft_method == 'fftw' ) THEN 1399 1400 #if defined( __fftw ) 1401 IF ( forward_fft ) THEN 1402 1403 y_in(0:ny) = ar(0:ny) 1404 CALL FFTW_EXECUTE_DFT_R2C( plan_yf, y_in, y_out ) 1405 1406 DO j = 0, (ny+1)/2 1407 ar(j) = REAL( y_out(j) ) / (ny+1) 1408 ENDDO 1409 DO j = 1, (ny+1)/2 - 1 1410 ar(ny+1-j) = AIMAG( y_out(j) ) / (ny+1) 1411 ENDDO 1412 1413 ELSE 1414 1415 y_out(0) = CMPLX( ar(0), 0.0 ) 1416 DO j = 1, (ny+1)/2 - 1 1417 y_out(j) = CMPLX( ar(j), ar(ny+1-j) ) 1418 ENDDO 1419 y_out((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 ) 1420 1421 CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in ) 1422 ar(0:ny) = y_in(0:ny) 1423 1424 ENDIF 1425 #endif 1324 1426 1325 1427 ELSEIF ( fft_method == 'system-specific' ) THEN
Note: See TracChangeset
for help on using the changeset viewer.
