Changeset 1216 for palm/trunk/SOURCE


Ignore:
Timestamp:
Aug 26, 2013 9:31:42 AM (11 years ago)
Author:
raasch
Message:

overlapping execution of fft and transpositions (MPI_ALLTOALL), but real overlapping is not activated so far,
fftw implemented for 1D-decomposition
resorting of arrays moved to separate routines resort_for_...
bugfix in mbuild concerning Makefile_check

Location:
palm/trunk/SOURCE
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile_check

    r1213 r1216  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# +tridia_solver
    2323#
    2424# Former revisions:
     
    7171      local_system.f90 message.f90 modules.f90 package_parin.f90 parin.f90 \
    7272      poisfft.f90 random_function.f90 singleton.f90 \
    73       subsidence.f90 temperton_fft.f90 \
     73      subsidence.f90 temperton_fft.f90 tridia_solver.f90 \
    7474      user_3d_data_averaging.f90 user_actions.f90 \
    7575      user_additional_routines.f90 user_check_data_output.f90 \
     
    8282      user_lpm_init.f90 user_lpm_set_attributes.f90 user_module.f90 \
    8383      user_parin.f90 user_read_restart_data.f90 user_spectra.f90 \
    84       user_statistics.f90 \
    85 
     84      user_statistics.f90
    8685
    8786OBJS=$(SOURCES:.f90=.o)
     
    131130package_parin.o: modules.o
    132131parin.o: modules.o
    133 poisfft.o: modules.o fft_xy.o
     132poisfft.o: modules.o fft_xy.o tridia_solver.o
    134133random_function.o: modules.o
    135134singleton.o: singleton.f90
    136135subsidence.o: modules.o
    137136temperton_fft.o: modules.o
     137tridia_solver.o: modules.o
    138138user_3d_data_averaging.o: modules.o user_module.o
    139139user_actions.o: modules.o user_module.o
  • palm/trunk/SOURCE/calc_spectra.f90

    r1121 r1216  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! resorting of array moved to separate routine resort_for_zx,
     23! one argument removed from the transpose_..d routines
    2324!
    2425! Former revisions:
     
    121122#if defined( __parallel )
    122123          IF ( pdims(2) /= 1 )  THEN
    123              CALL transpose_zx( d, tend, d )
     124             CALL resort_for_zx( d, tend )
     125             CALL transpose_zx( tend, d )
    124126          ELSE
    125              CALL transpose_yxd( d, tend, d )
     127             CALL transpose_yxd( d, d )
    126128          ENDIF
    127129          CALL calc_spectra_x( d, pr, m )
     
    153155
    154156#if defined( __parallel )
    155           CALL transpose_zyd( d, tend, d )
     157          CALL transpose_zyd( d, d )
    156158          CALL calc_spectra_y( d, pr, m )
    157159#else
  • palm/trunk/SOURCE/check_parameters.f90

    r1215 r1216  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! check for transpose_compute_overlap (temporary)
    2323!
    2424! Former revisions:
     
    319319    LOGICAL ::  found, ldum
    320320    REAL    ::  gradient, remote = 0.0, simulation_time_since_reference
     321
     322!
     323!-- Check for overlap combinations, which are not realized yet
     324    IF ( transpose_compute_overlap )  THEN
     325       IF ( numprocs == 1 )  STOP '+++ transpose-compute-overlap not implemented for single PE runs'
     326#if defined( __openacc )
     327       STOP '+++ transpose-compute-overlap not implemented for GPU usage'
     328#endif
     329    ENDIF
    321330
    322331!
  • palm/trunk/SOURCE/fft_xy.f90

    r1211 r1216  
    286286
    287287
    288     SUBROUTINE fft_x( ar, direction )
     288    SUBROUTINE fft_x( ar, direction, ar_2d )
    289289
    290290!----------------------------------------------------------------------!
     
    321321       COMPLEX(dpk), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) ::  ar_tmp
    322322#endif
     323       REAL, DIMENSION(0:nx,nys_x:nyn_x), OPTIONAL   ::  ar_2d
    323324       REAL, DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) ::  ar
    324325
     
    454455                   CALL FFTW_EXECUTE_DFT_R2C( plan_xf, x_in, x_out )
    455456
    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
    462476
    463477                ENDDO
     
    465479             !$OMP END PARALLEL
    466480
    467          ELSE
     481          ELSE
    468482             !$OMP PARALLEL PRIVATE ( work, i, j, k )
    469483             !$OMP DO
     
    471485                DO  j = nys_x, nyn_x
    472486
    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
    478504
    479505                   CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in)
     
    484510             !$OMP END PARALLEL
    485511
    486          ENDIF
     512          ENDIF
    487513#endif
    488514
     
    766792          ENDIF
    767793
     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
    768823       ELSEIF ( fft_method == 'system-specific' )  THEN
    769824
     
    843898    END SUBROUTINE fft_x_1d
    844899
    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 )
    846902
    847903!----------------------------------------------------------------------!
     
    853909!      fft_y uses internal algorithms (Singleton or Temperton) or      !
    854910!           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.                 !
    855926!----------------------------------------------------------------------!
    856927
     
    864935       CHARACTER (LEN=*) ::  direction
    865936       INTEGER ::  i, j, jshape(1), k
     937       INTEGER ::  nxl_y_bound, nxl_y_l, nxr_y_bound, nxr_y_l
    866938
    867939       LOGICAL ::  forward_fft
     
    878950       COMPLEX(dpk), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) ::  ar_tmp
    879951#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
    881954
    882955       IF ( direction == 'forward' )  THEN
     
    898971             !$OMP DO
    899972             DO  k = nzb_y, nzt_y
    900                 DO  i = nxl_y, nxr_y
     973                DO  i = nxl_y_l, nxr_y_l
    901974
    902975                   DO  j = 0, ny
     
    908981
    909982                   DO  j = 0, (ny+1)/2
    910                       ar(j,i,k) = REAL( cwork(j) )
     983                      ar_tr(j,i,k) = REAL( cwork(j) )
    911984                   ENDDO
    912985                   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) )
    914987                   ENDDO
    915988
     
    923996             !$OMP DO
    924997             DO  k = nzb_y, nzt_y
    925                 DO  i = nxl_y, nxr_y
    926 
    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 )
    9281001                   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 )
    9331006
    9341007                   jshape = SHAPE( cwork )
     
    9571030             !$OMP DO
    9581031             DO  k = nzb_y, nzt_y
    959                 DO  i = nxl_y, nxr_y
     1032                DO  i = nxl_y_l, nxr_y_l
    9601033
    9611034                   work(0:ny) = ar(0:ny,i,k)
     
    9631036
    9641037                   DO  j = 0, (ny+1)/2
    965                       ar(j,i,k) = work(2*j)
     1038                      ar_tr(j,i,k) = work(2*j)
    9661039                   ENDDO
    9671040                   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)
    9691042                   ENDDO
    9701043
     
    9781051             !$OMP DO
    9791052             DO  k = nzb_y, nzt_y
    980                 DO  i = nxl_y, nxr_y
     1053                DO  i = nxl_y_l, nxr_y_l
    9811054
    9821055                   DO  j = 0, (ny+1)/2
    983                       work(2*j) = ar(j,i,k)
     1056                      work(2*j) = ar_tr(j,i,k)
    9841057                   ENDDO
    9851058                   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)
    9871060                   ENDDO
    9881061                   work(1)    = 0.0
     
    10061079             !$OMP DO
    10071080             DO  k = nzb_y, nzt_y
    1008                 DO  i = nxl_y, nxr_y
     1081                DO  i = nxl_y_l, nxr_y_l
    10091082
    10101083                   y_in(0:ny) = ar(0:ny,i,k)
     
    10121085
    10131086                   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)
    10151088                   ENDDO
    10161089                   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)
    10181091                   ENDDO
    10191092
     
    10271100             !$OMP DO
    10281101             DO  k = nzb_y, nzt_y
    1029                 DO  i = nxl_y, nxr_y
    1030 
    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 )
    10321105                   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 )
    10361109
    10371110                   CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
     
    10531126             !$OMP DO
    10541127             DO  k = nzb_y, nzt_y
    1055                 DO  i = nxl_y, nxr_y
     1128                DO  i = nxl_y_l, nxr_y_l
    10561129
    10571130                   CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
     
    10591132
    10601133                   DO  j = 0, (ny+1)/2
    1061                       ar(j,i,k) = work(2*j)
     1134                      ar_tr(j,i,k) = work(2*j)
    10621135                   ENDDO
    10631136                   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)
    10651138                   ENDDO
    10661139
     
    10741147             !$OMP DO
    10751148             DO  k = nzb_y, nzt_y
    1076                 DO  i = nxl_y, nxr_y
     1149                DO  i = nxl_y_l, nxr_y_l
    10771150
    10781151                   DO  j = 0, (ny+1)/2
    1079                       work(2*j) = ar(j,i,k)
     1152                      work(2*j) = ar_tr(j,i,k)
    10801153                   ENDDO
    10811154                   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)
    10831156                   ENDDO
    10841157                   work(1)    = 0.0
     
    11031176             !$OMP DO
    11041177             DO  k = nzb_y, nzt_y
    1105                 DO  i = nxl_y, nxr_y
     1178                DO  i = nxl_y_l, nxr_y_l
    11061179
    11071180                   work(0:ny) = ar(0:ny,i,k)
     
    11101183
    11111184                   DO  j = 0, (ny+1)/2
    1112                       ar(j,i,k) = work(2*j)
     1185                      ar_tr(j,i,k) = work(2*j)
    11131186                   ENDDO
    11141187                   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)
    11161189                   ENDDO
    11171190
     
    11251198             !$OMP DO
    11261199             DO  k = nzb_y, nzt_y
    1127                 DO  i = nxl_y, nxr_y
     1200                DO  i = nxl_y_l, nxr_y_l
    11281201
    11291202                   DO  j = 0, (ny+1)/2
    1130                       work(2*j) = ar(j,i,k)
     1203                      work(2*j) = ar_tr(j,i,k)
    11311204                   ENDDO
    11321205                   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)
    11341207                   ENDDO
    11351208                   work(1) = 0.0
     
    13221395
    13231396          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
    13241426
    13251427       ELSEIF ( fft_method == 'system-specific' )  THEN
  • palm/trunk/SOURCE/header.f90

    r1213 r1216  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! output for transpose_compute_overlap
    2323!
    2424! Former revisions:
     
    332332    IF ( psolver(1:7) == 'poisfft' )  THEN
    333333       WRITE ( io, 111 )  TRIM( fft_method )
     334       IF ( transpose_compute_overlap )  WRITE( io, 115 )
    334335    ELSEIF ( psolver == 'sor' )  THEN
    335336       WRITE ( io, 112 )  nsor_ini, nsor, omega_sor
     
    16531654113 FORMAT (' --> Momentum advection via Piascek-Williams-Scheme (Form C3)', &
    16541655                  ' or Upstream')
     1656115 FORMAT ('     FFT and transpositions are overlapping')
    16551657116 FORMAT (' --> Scalar advection via Piascek-Williams-Scheme (Form C3)', &
    16561658                  ' or Upstream')
  • palm/trunk/SOURCE/modules.f90

    r1213 r1216  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! +transpose_compute_overlap,
     23! several variables are now defined in the serial (non-parallel) case also
    2324!
    2425! Former revisions:
     
    776777                scalar_rayleigh_damping = .TRUE., sloping_surface = .FALSE., &
    777778                stop_dt = .FALSE., synchronous_exchange = .FALSE., &
    778                 terminate_run = .FALSE., turbulence = .FALSE., &
    779                 turbulent_inflow = .FALSE., use_cmax = .TRUE., &
    780                 use_initial_profile_as_reference = .FALSE., &
     779                terminate_run = .FALSE., transpose_compute_overlap = .FALSE., &
     780                turbulence = .FALSE., turbulent_inflow = .FALSE., &
     781                use_cmax = .TRUE., use_initial_profile_as_reference = .FALSE., &
    781782                use_prescribed_profile_data = .FALSE., &
    782783                use_single_reference_value = .FALSE., &
     
    14891490    CHARACTER(LEN=2) ::  send_receive = 'al'
    14901491    CHARACTER(LEN=5) ::  myid_char = ''
    1491     INTEGER          ::  acc_rank, id_inflow = 0, id_recycling = 0,      &
    1492                          myid = 0, num_acc_per_node = 0, req_count = 0,  &
    1493                          target_id, npex = -1, npey = -1, numprocs = 1,  &
    1494                          numprocs_previous_run = -1,                     &
    1495                          tasks_per_node = -9999, threads_per_task = 1
     1492    INTEGER          ::  acc_rank, comm1dx, comm1dy, comm2d, comm_inter,       &
     1493                         comm_palm, id_inflow = 0, id_recycling = 0, ierr,     &
     1494                         myid = 0, myidx = 0, myidy = 0, ndim = 2, ngp_a,      &
     1495                         ngp_o, ngp_xy, ngp_y, npex = -1, npey = -1,           &
     1496                         numprocs = 1, numprocs_previous_run = -1,             &
     1497                         num_acc_per_node = 0, pleft, pnorth, pright, psouth,  &
     1498                         req_count = 0, sendrecvcount_xy, sendrecvcount_yz,    &
     1499                         sendrecvcount_zx, sendrecvcount_zyd,                  &
     1500                         sendrecvcount_yxd, target_id, tasks_per_node = -9999, &
     1501                         threads_per_task = 1, type_x, type_x_int, type_xy,    &
     1502                         type_y, type_y_int
    14961503
    14971504    INTEGER          ::  pdims(2) = 1, req(100)
     
    15071514    CHARACTER (LEN=MPI_MAX_PORT_NAME) ::  port_name
    15081515#endif
    1509 
    1510     INTEGER ::  comm1dx, comm1dy, comm2d, comm_inter, comm_palm, ierr, myidx,  &
    1511                 myidy, ndim = 2, ngp_a, ngp_o, ngp_xy, ngp_y, pleft, pnorth,   &
    1512                 pright, psouth, sendrecvcount_xy, sendrecvcount_yz,            &
    1513                 sendrecvcount_zx, sendrecvcount_zyd, sendrecvcount_yxd,        &
    1514                 type_x, type_x_int, type_xy, type_y, type_y_int
    15151516
    15161517    INTEGER ::  ibuf(12), pcoord(2)
  • palm/trunk/SOURCE/parin.f90

    r1196 r1216  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! +transpose_compute_overlap in inipar
    2323!
    2424! Former revisions:
     
    256256             subs_vertical_gradient_level, surface_heatflux, surface_pressure, &
    257257             surface_scalarflux, surface_waterflux, &
    258              s_surface, &
    259              s_surface_initial_change, s_vertical_gradient, &
     258             s_surface, s_surface_initial_change, s_vertical_gradient, &
    260259             s_vertical_gradient_level, timestep_scheme, &
    261260             topography, topography_grid_convention, top_heatflux, &
    262261             top_momentumflux_u, top_momentumflux_v, top_salinityflux, &
    263              turbulence, turbulent_inflow, ug_surface, ug_vertical_gradient, &
     262             transpose_compute_overlap, turbulence, turbulent_inflow, &
     263             ug_surface, ug_vertical_gradient, &
    264264             ug_vertical_gradient_level, use_surface_fluxes, use_cmax, &
    265265             use_top_fluxes, use_ug_for_galilei_tr, use_upstream_for_tke, &
  • palm/trunk/SOURCE/poisfft.f90

    r1213 r1216  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! resorting of arrays moved to separate routines resort_for_...,
     23! one argument, used as temporary work array, removed from all transpose
     24! routines
     25! overlapping fft / transposition implemented
    2326!
    2427! Former revisions:
     
    221224
    222225#if ! defined ( __check )
    223     SUBROUTINE poisfft( ar, work )
    224 
    225        USE control_parameters,  ONLY : fft_method
     226    SUBROUTINE poisfft( ar, ar_inv_test )
     227
     228       USE control_parameters,  ONLY : fft_method, transpose_compute_overlap
    226229       USE cpulog
    227230       USE interfaces
     
    230233       IMPLICIT NONE
    231234
    232        REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar, work
     235       INTEGER ::  ind_even, ind_odd, ind_third, ii, iind, inew, jj, jind,  &
     236                   jnew, ki, kk, knew, n, nblk, nnx_y, nny_z, nnz_t, nnz_x, &
     237                   nxl_y_bound, nxr_y_bound
     238       INTEGER, DIMENSION(4) ::  isave
     239
     240       REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar
     241       REAL, DIMENSION(nys:nyn,nxl:nxr,1:nz) ::  ar_inv_test ! work array tend from pres
     242       !$acc declare create( ar_inv )
     243       REAL, DIMENSION(nys:nyn,nxl:nxr,1:nz) ::  ar_inv
     244
     245       REAL, DIMENSION(:,:,:),   ALLOCATABLE ::  f_in, f_inv, f_out_y, f_out_z
     246       REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  ar1
    233247
    234248
     
    244258!--       1d-domain-decomposition along x:
    245259!--       FFT along y and transposition y --> x
    246           CALL ffty_tr_yx( ar, work, ar )
     260          CALL ffty_tr_yx( ar, ar )
    247261
    248262!
     
    252266!
    253267!--       Transposition x --> y and backward FFT along y
    254           CALL tr_xy_ffty( ar, work, ar )
     268          CALL tr_xy_ffty( ar, ar )
    255269
    256270       ELSEIF ( pdims(1) == 1  .AND.  pdims(2) > 1 )  THEN
     
    259273!--       1d-domain-decomposition along y:
    260274!--       FFT along x and transposition x --> y
    261           CALL fftx_tr_xy( ar, work, ar )
     275          CALL fftx_tr_xy( ar, ar )
    262276
    263277!
     
    267281!
    268282!--       Transposition y --> x and backward FFT along x
    269           CALL tr_yx_fftx( ar, work, ar )
    270 
    271        ELSE
     283          CALL tr_yx_fftx( ar, ar )
     284
     285       ELSEIF ( .NOT. transpose_compute_overlap )  THEN
    272286
    273287!
     
    275289!--       Transposition z --> x
    276290          CALL cpu_log( log_point_s(5), 'transpo forward', 'start' )
    277           CALL transpose_zx( ar, work, ar )
     291          CALL resort_for_zx( ar, ar_inv )
     292          CALL transpose_zx( ar_inv, ar )
    278293          CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
    279294
     
    291306!--       Transposition x --> y
    292307          CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
    293           CALL transpose_xy( ar, work, ar )
     308          CALL resort_for_xy( ar, ar_inv )
     309          CALL transpose_xy( ar_inv, ar )
    294310          CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
    295311
     
    298314             !$acc update host( ar )
    299315          ENDIF
    300           CALL fft_y( ar, 'forward' )
     316          CALL fft_y( ar, 'forward', ar_tr = ar,                &
     317                      nxl_y_bound = nxl_y, nxr_y_bound = nxr_y, &
     318                      nxl_y_l = nxl_y, nxr_y_l = nxr_y )
    301319          IF ( fft_method /= 'system-specific' )  THEN
    302320             !$acc update device( ar )
     
    307325!--       Transposition y --> z
    308326          CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
    309           CALL transpose_yz( ar, work, ar )
     327          CALL resort_for_yz( ar, ar_inv )
     328          CALL transpose_yz( ar_inv, ar )
    310329          CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
    311330
     
    320339!--       Transposition z --> y
    321340          CALL cpu_log( log_point_s(8), 'transpo invers', 'start' )
    322           CALL transpose_zy( ar, work, ar )
     341          CALL transpose_zy( ar, ar_inv )
     342          CALL resort_for_zy( ar_inv, ar )
    323343          CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
    324344
     
    327347             !$acc update host( ar )
    328348          ENDIF
    329           CALL fft_y( ar, 'backward' )
     349          CALL fft_y( ar, 'backward', ar_tr = ar,               &
     350                      nxl_y_bound = nxl_y, nxr_y_bound = nxr_y, &
     351                      nxl_y_l = nxl_y, nxr_y_l = nxr_y )
    330352          IF ( fft_method /= 'system-specific' )  THEN
    331353             !$acc update device( ar )
     
    336358!--       Transposition y --> x
    337359          CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
    338           CALL transpose_yx( ar, work, ar )
     360          CALL transpose_yx( ar, ar_inv )
     361          CALL resort_for_yx( ar_inv, ar )
    339362          CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
    340363
     
    352375!--       Transposition x --> z
    353376          CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
    354           CALL transpose_xz( ar, work, ar )
     377          CALL transpose_xz( ar, ar_inv )
     378          CALL resort_for_xz( ar_inv, ar )
    355379          CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
    356380
     381       ELSE
     382
     383!
     384!--       2d-domain-decomposition or no decomposition (1 PE run) with
     385!--       overlapping transposition / fft
     386          ALLOCATE( f_out_y(0:ny,nxl_y:nxr_y,nzb_y:nzt_y), &
     387                    f_out_z(0:nx,nys_x:nyn_x,nzb_x:nzt_x) )
     388!
     389!--       Transposition z --> x + subsequent fft along x
     390          ALLOCATE( f_inv(nys:nyn,nxl:nxr,1:nz) )
     391          CALL resort_for_zx( ar, f_inv )
     392!
     393!--       Save original indices and gridpoint counter
     394          isave(1) = nz
     395          isave(2) = nzb_x
     396          isave(3) = nzt_x
     397          isave(4) = sendrecvcount_zx
     398!
     399!--       Set new indices for transformation
     400          nblk  = nz / pdims(1)
     401          nz    = pdims(1)
     402          nnz_x = 1
     403          nzb_x = 1 + myidx * nnz_x
     404          nzt_x = ( myidx + 1 ) * nnz_x
     405          sendrecvcount_zx = nnx * nny * nnz_x
     406
     407          ALLOCATE( ar1(0:nx,nys_x:nyn_x,nzb_x:nzt_x,2) )
     408          ALLOCATE( f_in(nys:nyn,nxl:nxr,1:nz) )
     409
     410          DO  kk = 1, nblk+1
     411             ind_odd  = MOD( kk,   2 ) + 1
     412             ind_even = MOD( kk+1, 2 ) + 1
     413!$OMP sections private(ki,knew,n)
     414!$OMP section
     415             IF ( kk <= nblk )  THEN
     416
     417                IF ( kk == 1 )  THEN
     418                   CALL cpu_log( log_point_s(5), 'transpo forward', 'start' )
     419                ELSE
     420                   CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
     421                ENDIF
     422
     423                DO  knew = 1, nz
     424                   ki = kk + nblk * ( knew - 1 )
     425                   f_in(:,:,knew) = f_inv(:,:,ki)
     426                ENDDO
     427
     428                CALL transpose_zx( f_in, ar1(:,:,:,ind_odd))
     429                CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
     430
     431             ENDIF
     432
     433!$OMP section
     434             IF ( kk >= 2 )  THEN
     435
     436                IF ( kk == 2 )  THEN
     437                   CALL cpu_log( log_point_s(4), 'fft_x', 'start' )
     438                ELSE
     439                   CALL cpu_log( log_point_s(4), 'fft_x', 'continue' )
     440                ENDIF
     441
     442                n = isave(2) + kk - 2
     443                CALL fft_x( ar1(:,:,:,ind_even), 'forward',  ar_2d = f_out_z(:,:,n))
     444                CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
     445
     446             ENDIF
     447!$OMP end sections
     448
     449          ENDDO
     450!
     451!--       Restore original indices/counters
     452          nz               = isave(1)
     453          nzb_x            = isave(2)
     454          nzt_x            = isave(3)
     455          sendrecvcount_zx = isave(4)
     456
     457          DEALLOCATE( ar1, f_in, f_inv )
     458
     459!
     460!--       Transposition x --> y + subsequent fft along y
     461          ALLOCATE( f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
     462          CALL resort_for_xy( f_out_z, f_inv )
     463!
     464!--       Save original indices and gridpoint counter
     465          isave(1) = nx
     466          isave(2) = nxl_y
     467          isave(3) = nxr_y
     468          isave(4) = sendrecvcount_xy
     469!
     470!--       Set new indices for transformation
     471          nblk  = ( ( nx+1 ) / pdims(2) ) - 1
     472          nx    = pdims(2)
     473          nnx_y = 1
     474          nxl_y = myidy * nnx_y
     475          nxr_y = ( myidy + 1 ) * nnx_y - 1
     476          sendrecvcount_xy = nnx_y * ( nyn_x-nys_x+1 ) * ( nzt_x-nzb_x+1 )
     477
     478          ALLOCATE( ar1(0:ny,nxl_y:nxr_y,nzb_y:nzt_y,2) )
     479          ALLOCATE( f_in(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
     480
     481          DO  ii = 0, nblk+1
     482             ind_odd  = MOD( ii+1, 2 ) + 1
     483             ind_even = MOD( ii+2, 2 ) + 1
     484!$OMP sections private(ki,knew,n)
     485!$OMP section
     486             IF ( ii <= nblk )  THEN
     487
     488                CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
     489
     490                DO  inew = 0, nx-1
     491                   iind = ii + ( nblk + 1 ) * inew
     492                   f_in(:,:,inew) = f_inv(:,:,iind)
     493                ENDDO
     494
     495                CALL transpose_xy( f_in, ar1(:,:,:,ind_odd) )
     496
     497                CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
     498
     499             ENDIF
     500
     501!$OMP section
     502             IF ( ii >= 1 )  THEN
     503
     504                IF ( ii == 1 )  THEN
     505                   CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
     506                ELSE
     507                   CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
     508                ENDIF
     509
     510                nxl_y_bound = isave(2)
     511                nxr_y_bound = isave(3)
     512                n           = isave(2) + ii - 1
     513!                CALL fft_y( ar1(:,:,:,ind_even), 'forward', ar_3d = f_out_y, &
     514!                            ni = n )
     515                CALL fft_y( ar1(:,:,:,ind_even), 'forward', ar_tr = f_out_y, &
     516                            nxl_y_bound = nxl_y_bound, nxr_y_bound = nxr_y_bound, &
     517                            nxl_y_l = n, nxr_y_l = n )
     518
     519                CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
     520
     521             ENDIF
     522!$OMP end sections
     523
     524          ENDDO
     525!
     526!--       Restore original indices/counters
     527          nx               = isave(1)
     528          nxl_y            = isave(2)
     529          nxr_y            = isave(3)
     530          sendrecvcount_xy = isave(4)
     531
     532          DEALLOCATE( ar1, f_in, f_inv )
     533
     534!
     535!--       Transposition y --> z + subsequent tridia + resort for z --> y
     536          ALLOCATE( f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) )
     537          CALL resort_for_yz( f_out_y, f_inv )
     538!
     539!--       Save original indices and gridpoint counter
     540          isave(1) = ny
     541          isave(2) = nys_z
     542          isave(3) = nyn_z
     543          isave(4) = sendrecvcount_yz
     544!
     545!--       Set new indices for transformation
     546          nblk             = ( ( ny+1 ) / pdims(1) ) - 1
     547          ny               = pdims(1)
     548          nny_z            = 1
     549          nys_z            = myidx * nny_z
     550          nyn_z            = ( myidx + 1 ) * nny_z - 1
     551          sendrecvcount_yz = ( nxr_y-nxl_y+1 ) * nny_z * ( nzt_y-nzb_y+1 )
     552
     553          ALLOCATE( ar1(nxl_z:nxr_z,nys_z:nyn_z,1:nz,3) )
     554          ALLOCATE( f_in(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) )
     555
     556          DO  jj = 0, nblk+2
     557             ind_odd   = MOD( jj+3, 3 ) + 1
     558             ind_even  = MOD( jj+2, 3 ) + 1
     559             ind_third = MOD( jj+1, 3 ) + 1
     560!$OMP sections private(ki,knew,n)
     561!$OMP section
     562             IF ( jj <= nblk )  THEN
     563!
     564!--             Forward Fourier Transformation
     565!--             Transposition y --> z
     566                CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
     567
     568                DO  jnew = 0, ny-1
     569                   jind = jj + ( nblk + 1 ) * jnew
     570                   f_in(:,:,jnew) =f_inv(:,:,jind)
     571                ENDDO
     572
     573                CALL transpose_yz( f_in, ar1(:,:,:,ind_odd) )
     574
     575                IF ( jj == nblk )  THEN
     576                   CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
     577                ELSE
     578                   CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
     579                ENDIF
     580
     581             ENDIF
     582
     583             IF ( jj >= 2 )  THEN
     584!
     585!--             Inverse Fourier Transformation
     586!--             Transposition z --> y
     587!--             Only one thread should call MPI routines, therefore forward and
     588!--             backward tranpose are in the same section
     589                IF ( jj == 2 )  THEN
     590                   CALL cpu_log( log_point_s(8), 'transpo invers', 'start' )
     591                ELSE
     592                   CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
     593                ENDIF
     594
     595                CALL transpose_zy( ar1(:,:,:,ind_third), f_in )
     596
     597                DO  jnew = 0, ny-1
     598                   jind = jj-2 + ( nblk + 1 ) * jnew
     599                   f_inv(:,:,jind) = f_in(:,:,jnew)
     600                ENDDO
     601
     602                CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
     603
     604             ENDIF
     605
     606!$OMP section
     607             IF ( jj >= 1  .AND.  jj <= nblk+1 )  THEN
     608!
     609!--             Solve the tridiagonal equation system along z
     610                CALL cpu_log( log_point_s(6), 'tridia', 'start' )
     611
     612                n = isave(2) + jj - 1
     613                CALL tridia_substi_overlap( ar1(:,:,:,ind_even), n )
     614
     615                CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
     616             ENDIF
     617!$OMP end sections
     618
     619          ENDDO
     620!
     621!--       Restore original indices/counters
     622          ny               = isave(1)
     623          nys_z            = isave(2)
     624          nyn_z            = isave(3)
     625          sendrecvcount_yz = isave(4)
     626
     627          CALL resort_for_zy( f_inv, f_out_y )
     628
     629          DEALLOCATE( ar1, f_in, f_inv )
     630
     631!
     632!--       fft along y backward + subsequent transposition y --> x
     633          ALLOCATE( f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
     634!
     635!--       Save original indices and gridpoint counter
     636          isave(1) = nx
     637          isave(2) = nxl_y
     638          isave(3) = nxr_y
     639          isave(4) = sendrecvcount_xy
     640!
     641!--       Set new indices for transformation
     642          nblk             = (( nx+1 ) / pdims(2) ) - 1
     643          nx               = pdims(2)
     644          nnx_y            = 1
     645          nxl_y            = myidy * nnx_y
     646          nxr_y            = ( myidy + 1 ) * nnx_y - 1
     647          sendrecvcount_xy = nnx_y * ( nyn_x-nys_x+1 ) * ( nzt_x-nzb_x+1 )
     648
     649          ALLOCATE( ar1(0:ny,nxl_y:nxr_y,nzb_y:nzt_y,2) )
     650          ALLOCATE( f_in(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
     651
     652          DO  ii = 0, nblk+1
     653             ind_odd  = MOD( ii+1, 2 ) + 1
     654             ind_even = MOD( ii+2, 2 ) + 1
     655!$OMP sections private(ki,knew,n)
     656!$OMP section
     657             IF ( ii <= nblk )  THEN
     658
     659                CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
     660
     661                n = isave(2) + ii
     662                nxl_y_bound = isave(2)
     663                nxr_y_bound = isave(3)
     664
     665!                CALL fft_y( ar1(:,:,:,ind_even), 'backward', ar_3d = f_out_y, &
     666!                            ni = n )
     667                CALL fft_y( ar1(:,:,:,ind_even), 'backward', ar_tr = f_out_y, &
     668                            nxl_y_bound = nxl_y_bound, nxr_y_bound = nxr_y_bound, &
     669                            nxl_y_l = n, nxr_y_l = n )
     670
     671                IF ( ii == nblk )  THEN
     672                   CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
     673                ELSE
     674                   CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
     675                ENDIF
     676
     677             ENDIF
     678
     679!$OMP section
     680             IF ( ii >= 1 )  THEN
     681
     682                CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
     683
     684                CALL transpose_yx( ar1(:,:,:,ind_odd), f_in )
     685
     686                DO  inew = 0, nx-1
     687                   iind = ii-1 + (nblk+1) * inew
     688                   f_inv(:,:,iind) = f_in(:,:,inew)
     689                ENDDO
     690
     691                CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
     692
     693             ENDIF
     694!$OMP end sections
     695
     696          ENDDO
     697!
     698!--       Restore original indices/counters
     699          nx               = isave(1)
     700          nxl_y            = isave(2)
     701          nxr_y            = isave(3)
     702          sendrecvcount_xy = isave(4)
     703
     704          CALL resort_for_yx( f_inv, f_out_z )
     705
     706          DEALLOCATE( ar1, f_in, f_inv )
     707
     708!
     709!--       fft along x backward + subsequent final transposition x --> z
     710          ALLOCATE( f_inv(nys:nyn,nxl:nxr,1:nz) )
     711!
     712!--       Save original indices and gridpoint counter
     713          isave(1) = nz
     714          isave(2) = nzb_x
     715          isave(3) = nzt_x
     716          isave(4) = sendrecvcount_zx
     717!
     718!--       Set new indices for transformation
     719          nblk             = nz / pdims(1)
     720          nz               = pdims(1)
     721          nnz_x            = 1
     722          nzb_x            = 1 + myidx * nnz_x
     723          nzt_x            = ( myidx + 1 ) * nnz_x
     724          sendrecvcount_zx = nnx * nny * nnz_x
     725
     726          ALLOCATE( ar1(0:nx,nys_x:nyn_x,nzb_x:nzt_x,2) )
     727          ALLOCATE( f_in(nys:nyn,nxl:nxr,1:nz) )
     728
     729          DO  kk = 1, nblk+1
     730             ind_odd  = MOD( kk,   2 ) + 1
     731             ind_even = MOD( kk+1, 2 ) + 1
     732!$OMP sections private(ki,knew,n)
     733!$OMP section
     734             IF ( kk <= nblk )  THEN
     735
     736                CALL cpu_log( log_point_s(4), 'fft_x', 'continue' )
     737
     738                n = isave(2) + kk - 1
     739                CALL fft_x( ar1(:,:,:,ind_even), 'backward', f_out_z(:,:,n))
     740
     741                IF ( kk == nblk )  THEN
     742                   CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
     743                ELSE
     744                   CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
     745                ENDIF
     746
     747             ENDIF
     748
     749!$OMP section
     750             IF ( kk >= 2 )  THEN
     751
     752                CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
     753
     754                CALL transpose_xz( ar1(:,:,:,ind_odd), f_in )
     755
     756                DO  knew = 1, nz
     757                   ki = kk-1 + nblk * (knew-1)
     758                   f_inv(:,:,ki) = f_in(:,:,knew)
     759                ENDDO
     760
     761                IF ( kk == nblk+1 )  THEN
     762                   CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
     763                ELSE
     764                   CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
     765                ENDIF
     766
     767             ENDIF
     768!$OMP end sections
     769
     770          ENDDO
     771!
     772!--       Restore original indices/counters
     773          nz               = isave(1)
     774          nzb_x            = isave(2)
     775          nzt_x            = isave(3)
     776          sendrecvcount_zx = isave(4)
     777
     778          CALL resort_for_xz( f_inv, ar )
     779
     780          DEALLOCATE( ar1, f_in, f_inv )
     781
    357782       ENDIF
    358783
     
    363788
    364789
    365     SUBROUTINE ffty_tr_yx( f_in, work, f_out )
     790    SUBROUTINE ffty_tr_yx( f_in, f_out )
    366791
    367792!------------------------------------------------------------------------------!
     
    489914
    490915
    491     SUBROUTINE tr_xy_ffty( f_in, work, f_out )
     916    SUBROUTINE tr_xy_ffty( f_in, f_out )
    492917
    493918!------------------------------------------------------------------------------!
     
    7381163
    7391164
    740     SUBROUTINE fftx_tr_xy( f_in, work, f_out )
     1165    SUBROUTINE fftx_tr_xy( f_in, f_out )
    7411166
    7421167!------------------------------------------------------------------------------!
     
    8481273
    8491274
    850     SUBROUTINE tr_yx_fftx( f_in, work, f_out )
     1275    SUBROUTINE tr_yx_fftx( f_in, f_out )
    8511276
    8521277!------------------------------------------------------------------------------!
  • palm/trunk/SOURCE/transpose.f90

    r1112 r1216  
    1  SUBROUTINE transpose_xy( f_in, work, f_out )
     1 SUBROUTINE resort_for_xy( f_in, f_inv )
    22
    33!--------------------------------------------------------------------------------!
     
    2020! Current revisions:
    2121! -----------------
    22 !
     22! re-sorting of the transposed / to be transposed arrays moved to separate
     23! routines resort_for_...
    2324!
    2425! Former revisions:
     
    6970! Initial revision
    7071!
    71 !
     72!------------------------------------------------------------------------------!
     73! Description:
     74! ------------
     75! Resorting data for the transposition from x to y. The transposition itself
     76! is carried out in transpose_xy
     77!------------------------------------------------------------------------------!
     78
     79     USE indices
     80     USE transpose_indices
     81
     82     IMPLICIT NONE
     83
     84     REAL ::  f_in(0:nx,nys_x:nyn_x,nzb_x:nzt_x)
     85     REAL ::  f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx)
     86
     87
     88     INTEGER ::  i, j, k
     89
     90!
     91!-- Rearrange indices of input array in order to make data to be send
     92!-- by MPI contiguous
     93    !$OMP  PARALLEL PRIVATE ( i, j, k )
     94    !$OMP  DO
     95    !$acc kernels present( f_in, f_inv )
     96    !$acc loop
     97     DO  i = 0, nx
     98         DO  k = nzb_x, nzt_x
     99             !$acc loop vector( 32 )
     100             DO  j = nys_x, nyn_x
     101                 f_inv(j,k,i) = f_in(i,j,k)
     102             ENDDO
     103         ENDDO
     104     ENDDO
     105     !$acc end kernels
     106     !$OMP  END PARALLEL
     107
     108 END SUBROUTINE resort_for_xy
     109
     110
     111 SUBROUTINE transpose_xy( f_inv, f_out )
     112
     113!------------------------------------------------------------------------------!
    72114! Description:
    73115! ------------
     
    87129    INTEGER ::  i, j, k, l, ys
    88130   
    89     REAL ::  f_in(0:nx,nys_x:nyn_x,nzb_x:nzt_x), f_out(0:ny,nxl_y:nxr_y,nzb_y:nzt_y)
     131    REAL ::  f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx), f_out(0:ny,nxl_y:nxr_y,nzb_y:nzt_y)
    90132
    91133    REAL, DIMENSION(nyn_x-nys_x+1,nzb_y:nzt_y,nxl_y:nxr_y,0:pdims(2)-1) ::  work
    92134
    93     !$acc declare create( f_inv )
    94     REAL ::  f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx)
    95 
    96 
    97 !
    98 !-- Rearrange indices of input array in order to make data to be send
    99 !-- by MPI contiguous
    100 !$OMP  PARALLEL PRIVATE ( i, j, k )
    101 !$OMP  DO
    102     !$acc kernels present( f_in )
    103     !$acc loop
    104     DO  i = 0, nx
    105        DO  k = nzb_x, nzt_x
    106           !$acc loop vector( 32 )
    107           DO  j = nys_x, nyn_x
    108              f_inv(j,k,i) = f_in(i,j,k)
    109           ENDDO
    110        ENDDO
    111     ENDDO
    112     !$acc end kernels
    113 !$OMP  END PARALLEL
    114135
    115136    IF ( numprocs /= 1 )  THEN
     
    124145                          work(1,nzb_y,nxl_y,0), sendrecvcount_xy, MPI_REAL, &
    125146                          comm1dy, ierr )
    126        !$acc update device( work )
    127147       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
    128148
     
    131151!$OMP  PARALLEL PRIVATE ( i, j, k, l, ys )
    132152!$OMP  DO
     153       !$acc data copyin( work )
    133154       DO  l = 0, pdims(2) - 1
    134155          ys = 0 + l * ( nyn_x - nys_x + 1 )
     
    145166          !$acc end kernels
    146167       ENDDO
     168       !$acc end data
    147169!$OMP  END PARALLEL
    148170#endif
     
    154176!$OMP  PARALLEL PRIVATE ( i, j, k )
    155177!$OMP  DO
    156        !$acc kernels present( f_out )
     178       !$acc kernels present( f_inv, f_out )
    157179       !$acc loop
    158180       DO  k = nzb_y, nzt_y
     
    172194
    173195
    174  SUBROUTINE transpose_xz( f_in, work, f_out )
     196 SUBROUTINE resort_for_xz( f_inv, f_out )
     197
     198!------------------------------------------------------------------------------!
     199! Description:
     200! ------------
     201! Resorting data after the transposition from x to z. The transposition itself
     202! is carried out in transpose_xz
     203!------------------------------------------------------------------------------!
     204
     205     USE indices
     206     USE transpose_indices
     207
     208     IMPLICIT NONE
     209
     210     REAL ::  f_inv(nys:nyn,nxl:nxr,1:nz)
     211     REAL ::  f_out(1:nz,nys:nyn,nxl:nxr)
     212
     213
     214     INTEGER ::  i, j, k
     215
     216!
     217!-- Rearrange indices of input array in order to make data to be send
     218!-- by MPI contiguous.
     219!-- In case of parallel fft/transposition, scattered store is faster in
     220!-- backward direction!!!
     221    !$OMP  PARALLEL PRIVATE ( i, j, k )
     222    !$OMP  DO
     223    !$acc kernels present( f_inv, f_out )
     224    !$acc loop
     225     DO  k = 1, nz
     226         DO  i = nxl, nxr
     227             !$acc loop vector( 32 )
     228             DO  j = nys, nyn
     229                 f_out(k,j,i) = f_inv(j,i,k)
     230             ENDDO
     231         ENDDO
     232     ENDDO
     233     !$acc end kernels
     234     !$OMP  END PARALLEL
     235
     236 END SUBROUTINE resort_for_xz
     237
     238
     239 SUBROUTINE transpose_xz( f_in, f_inv )
    175240
    176241!------------------------------------------------------------------------------!
     
    192257    INTEGER ::  i, j, k, l, xs
    193258   
    194     REAL ::  f_in(0:nx,nys_x:nyn_x,nzb_x:nzt_x), f_out(1:nz,nys:nyn,nxl:nxr)
     259    REAL ::  f_in(0:nx,nys_x:nyn_x,nzb_x:nzt_x), f_inv(nys:nyn,nxl:nxr,1:nz)
    195260
    196261    REAL, DIMENSION(nys_x:nyn_x,nnx,nzb_x:nzt_x,0:pdims(1)-1) ::  work
    197 
    198     !$acc declare create( f_inv )
    199     REAL ::  f_inv(nys:nyn,nxl:nxr,1:nz)
    200262
    201263
     
    210272!$OMP  PARALLEL PRIVATE ( i, j, k, l, xs )
    211273!$OMP  DO
     274       !$acc data copyout( work )
    212275       DO  l = 0, pdims(1) - 1
    213276          xs = 0 + l * nnx
     
    224287          !$acc end kernels
    225288       ENDDO
     289       !$acc end data
    226290!$OMP  END PARALLEL
    227291
     
    230294       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
    231295       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    232        !$acc update host( work )
    233296       CALL MPI_ALLTOALL( work(nys_x,1,nzb_x,0), sendrecvcount_zx, MPI_REAL, &
    234297                          f_inv(nys,nxl,1),      sendrecvcount_zx, MPI_REAL, &
     
    236299       !$acc update device( f_inv )
    237300       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
    238 
    239 !
    240 !--    Reorder transposed array in a way that the z index is in first position
    241 !$OMP  PARALLEL PRIVATE ( i, j, k )
    242 !$OMP  DO
    243        !$acc kernels present( f_out )
    244        !$acc loop
    245        DO  k = 1, nz
    246           DO  i = nxl, nxr
    247              !$acc loop vector( 32 )
    248              DO  j = nys, nyn
    249                 f_out(k,j,i) = f_inv(j,i,k)
    250              ENDDO
    251           ENDDO
    252        ENDDO
    253        !$acc end kernels
    254 !$OMP  END PARALLEL
    255301#endif
    256302
     
    261307!$OMP  PARALLEL PRIVATE ( i, j, k )
    262308!$OMP  DO
    263        !$acc kernels present( f_in )
     309       !$acc kernels present( f_in, f_inv )
    264310       !$acc loop
    265311       DO  i = nxl, nxr
     
    274320!$OMP  END PARALLEL
    275321
    276 !$OMP  PARALLEL PRIVATE ( i, j, k )
    277 !$OMP  DO
    278        !$acc kernels present( f_out )
    279        !$acc loop
    280        DO  k = 1, nz
    281           DO  i = nxl, nxr
    282              !$acc loop vector( 32 )
    283              DO  j = nys, nyn
    284                 f_out(k,j,i) = f_inv(j,i,k)
    285              ENDDO
    286           ENDDO
    287        ENDDO
    288        !$acc end kernels
    289 !$OMP  END PARALLEL
    290 
    291322    ENDIF
    292323
     
    294325
    295326
    296  SUBROUTINE transpose_yx( f_in, work, f_out )
     327 SUBROUTINE resort_for_yx( f_inv, f_out )
     328
     329!------------------------------------------------------------------------------!
     330! Description:
     331! ------------
     332! Resorting data after the transposition from y to x. The transposition itself
     333! is carried out in transpose_yx
     334!------------------------------------------------------------------------------!
     335
     336     USE indices
     337     USE transpose_indices
     338
     339     IMPLICIT NONE
     340
     341     REAL ::  f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx)
     342     REAL ::  f_out(0:nx,nys_x:nyn_x,nzb_x:nzt_x)
     343
     344
     345     INTEGER ::  i, j, k
     346
     347!
     348!-- Rearrange indices of input array in order to make data to be send
     349!-- by MPI contiguous
     350    !$OMP  PARALLEL PRIVATE ( i, j, k )
     351    !$OMP  DO
     352    !$acc kernels present( f_inv, f_out )
     353    !$acc loop
     354     DO  i = 0, nx
     355         DO  k = nzb_x, nzt_x
     356             !$acc loop vector( 32 )
     357             DO  j = nys_x, nyn_x
     358                 f_out(i,j,k) = f_inv(j,k,i)
     359             ENDDO
     360         ENDDO
     361     ENDDO
     362     !$acc end kernels
     363     !$OMP  END PARALLEL
     364
     365 END SUBROUTINE resort_for_yx
     366
     367
     368 SUBROUTINE transpose_yx( f_in, f_inv )
    297369
    298370!------------------------------------------------------------------------------!
     
    314386    INTEGER ::  i, j, k, l, ys
    315387   
    316     REAL ::  f_in(0:ny,nxl_y:nxr_y,nzb_y:nzt_y), f_out(0:nx,nys_x:nyn_x,nzb_x:nzt_x)
     388    REAL ::  f_in(0:ny,nxl_y:nxr_y,nzb_y:nzt_y), f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx)
    317389
    318390    REAL, DIMENSION(nyn_x-nys_x+1,nzb_y:nzt_y,nxl_y:nxr_y,0:pdims(2)-1) ::  work
    319 
    320     !$acc declare create( f_inv )
    321     REAL ::  f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx)
    322391
    323392
     
    329398!$OMP  PARALLEL PRIVATE ( i, j, k, l, ys )
    330399!$OMP  DO
     400       !$acc data copyout( work )
    331401       DO  l = 0, pdims(2) - 1
    332402          ys = 0 + l * ( nyn_x - nys_x + 1 )
     
    343413          !$acc end kernels
    344414       ENDDO
     415       !$acc end data
    345416!$OMP  END PARALLEL
    346417
     
    349420       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
    350421       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    351        !$acc update host( work )
    352422       CALL MPI_ALLTOALL( work(1,nzb_y,nxl_y,0), sendrecvcount_xy, MPI_REAL, &
    353423                          f_inv(nys_x,nzb_x,0),  sendrecvcount_xy, MPI_REAL, &
     
    363433!$OMP  PARALLEL PRIVATE ( i, j, k )
    364434!$OMP  DO
    365        !$acc kernels present( f_in )
     435       !$acc kernels present( f_in, f_inv )
    366436       !$acc loop
    367437       DO  i = nxl_y, nxr_y
     
    378448    ENDIF
    379449
    380 !
    381 !-- Reorder transposed array in a way that the x index is in first position
    382 !$OMP  PARALLEL PRIVATE ( i, j, k )
    383 !$OMP  DO
    384     !$acc kernels present( f_out )
    385     !$acc loop
    386     DO  i = 0, nx
    387        DO  k = nzb_x, nzt_x
    388           !$acc loop vector( 32 )
    389           DO  j = nys_x, nyn_x
    390              f_out(i,j,k) = f_inv(j,k,i)
    391           ENDDO
    392        ENDDO
    393     ENDDO
    394     !$acc end kernels
    395 !$OMP  END PARALLEL
    396 
    397450 END SUBROUTINE transpose_yx
    398451
    399452
    400  SUBROUTINE transpose_yxd( f_in, work, f_out )
     453 SUBROUTINE transpose_yxd( f_in, f_out )
    401454
    402455!------------------------------------------------------------------------------!
     
    466519
    467520
    468  SUBROUTINE transpose_yz( f_in, work, f_out )
     521 SUBROUTINE resort_for_yz( f_in, f_inv )
     522
     523!------------------------------------------------------------------------------!
     524! Description:
     525! ------------
     526! Resorting data for the transposition from y to z. The transposition itself
     527! is carried out in transpose_yz
     528!------------------------------------------------------------------------------!
     529
     530     USE indices
     531     USE transpose_indices
     532
     533     IMPLICIT NONE
     534
     535     REAL ::  f_in(0:ny,nxl_y:nxr_y,nzb_y:nzt_y)
     536     REAL ::  f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny)
     537
     538
     539     INTEGER ::  i, j, k
     540
     541!
     542!-- Rearrange indices of input array in order to make data to be send
     543!-- by MPI contiguous
     544    !$OMP  PARALLEL PRIVATE ( i, j, k )
     545    !$OMP  DO
     546    !$acc kernels present( f_in, f_inv )
     547    !$acc loop
     548     DO  j = 0, ny
     549         DO  k = nzb_y, nzt_y
     550             !$acc loop vector( 32 )
     551             DO  i = nxl_y, nxr_y
     552                 f_inv(i,k,j) = f_in(j,i,k)
     553             ENDDO
     554         ENDDO
     555     ENDDO
     556     !$acc end kernels
     557     !$OMP  END PARALLEL
     558
     559 END SUBROUTINE resort_for_yz
     560
     561
     562 SUBROUTINE transpose_yz( f_inv, f_out )
    469563
    470564!------------------------------------------------------------------------------!
     
    486580    INTEGER ::  i, j, k, l, zs
    487581   
    488     REAL ::  f_in(0:ny,nxl_y:nxr_y,nzb_y:nzt_y), f_out(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
     582    REAL ::  f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny), f_out(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
    489583
    490584    REAL, DIMENSION(nxl_z:nxr_z,nzt_y-nzb_y+1,nys_z:nyn_z,0:pdims(1)-1) ::  work
    491585
    492     !$acc declare create( f_inv )
    493     REAL ::  f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny)
    494 
    495 
    496 !
    497 !-- Rearrange indices of input array in order to make data to be send
    498 !-- by MPI contiguous
    499 !$OMP  PARALLEL PRIVATE ( i, j, k )
    500 !$OMP  DO
    501     !$acc kernels present( f_in )
    502     !$acc loop
    503     DO  j = 0, ny
    504        DO  k = nzb_y, nzt_y
    505           !$acc loop vector( 32 )
    506           DO  i = nxl_y, nxr_y
    507              f_inv(i,k,j) = f_in(j,i,k)
    508           ENDDO
    509        ENDDO
    510     ENDDO
    511     !$acc end kernels
    512 !$OMP  END PARALLEL
    513 
    514 !
    515 !-- Move data to different array, because memory location of work1 is
    516 !-- needed further below (work1 = work2).
     586
     587!
    517588!-- If the PE grid is one-dimensional along y, only local reordering
    518589!-- of the data is necessary and no transposition has to be done.
     
    521592!$OMP  PARALLEL PRIVATE ( i, j, k )
    522593!$OMP  DO
    523        !$acc kernels present( f_out )
     594       !$acc kernels present( f_inv, f_out )
    524595       !$acc loop
    525596       DO  j = 0, ny
     
    545616                          work(nxl_z,1,nys_z,0), sendrecvcount_yz, MPI_REAL, &
    546617                          comm1dx, ierr )
    547        !$acc update device( work )
    548618       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
    549619
     
    552622!$OMP  PARALLEL PRIVATE ( i, j, k, l, zs )
    553623!$OMP  DO
     624       !$acc data copyin( work )
    554625       DO  l = 0, pdims(1) - 1
    555626          zs = 1 + l * ( nzt_y - nzb_y + 1 )
    556           !$acc kernels present( f_out, work )
     627          !$acc kernels present( f_out )
    557628          !$acc loop
    558629          DO  j = nys_z, nyn_z
     
    566637          !$acc end kernels
    567638       ENDDO
     639       !$acc end data
    568640!$OMP  END PARALLEL
    569641#endif
     
    574646
    575647
    576  SUBROUTINE transpose_zx( f_in, work, f_out )
     648 SUBROUTINE resort_for_zx( f_in, f_inv )
     649
     650!------------------------------------------------------------------------------!
     651! Description:
     652! ------------
     653! Resorting data for the transposition from z to x. The transposition itself
     654! is carried out in transpose_zx
     655!------------------------------------------------------------------------------!
     656
     657     USE indices
     658     USE transpose_indices
     659
     660     IMPLICIT NONE
     661
     662     REAL ::  f_in(1:nz,nys:nyn,nxl:nxr)
     663     REAL ::  f_inv(nys:nyn,nxl:nxr,1:nz)
     664
     665
     666     INTEGER ::  i, j, k
     667
     668!
     669!-- Rearrange indices of input array in order to make data to be send
     670!-- by MPI contiguous
     671    !$OMP  PARALLEL PRIVATE ( i, j, k )
     672    !$OMP  DO
     673    !$acc kernels present( f_in, f_inv )
     674    !$acc loop
     675     DO  k = 1,nz
     676         DO  i = nxl, nxr
     677             !$acc loop vector( 32 )
     678             DO  j = nys, nyn
     679                 f_inv(j,i,k) = f_in(k,j,i)
     680             ENDDO
     681         ENDDO
     682     ENDDO
     683     !$acc end kernels
     684     !$OMP  END PARALLEL
     685
     686 END SUBROUTINE resort_for_zx
     687
     688
     689 SUBROUTINE transpose_zx( f_inv, f_out )
    577690
    578691!------------------------------------------------------------------------------!
     
    594707    INTEGER ::  i, j, k, l, xs
    595708   
    596     REAL ::  f_in(1:nz,nys:nyn,nxl:nxr), f_out(0:nx,nys_x:nyn_x,nzb_x:nzt_x)
     709    REAL ::  f_inv(nys:nyn,nxl:nxr,1:nz), f_out(0:nx,nys_x:nyn_x,nzb_x:nzt_x)
    597710
    598711    REAL, DIMENSION(nys_x:nyn_x,nnx,nzb_x:nzt_x,0:pdims(1)-1) ::  work
    599712
    600     !$acc declare create( f_inv )
    601     REAL ::  f_inv(nys:nyn,nxl:nxr,1:nz)
    602 
    603 
    604 !
    605 !-- Rearrange indices of input array in order to make data to be send
    606 !-- by MPI contiguous
    607 !$OMP  PARALLEL PRIVATE ( i, j, k )
    608 !$OMP  DO
    609     !$acc kernels present( f_in )
    610     !$acc loop
    611     DO  k = 1,nz
    612        DO  i = nxl, nxr
    613           !$acc loop vector( 32 )
    614           DO  j = nys, nyn
    615              f_inv(j,i,k) = f_in(k,j,i)
    616           ENDDO
    617        ENDDO
    618     ENDDO
    619     !$acc end kernels
    620 !$OMP  END PARALLEL
    621 
    622 !
    623 !-- Move data to different array, because memory location of work1 is
    624 !-- needed further below (work1 = work2).
     713
     714!
    625715!-- If the PE grid is one-dimensional along y, only local reordering
    626716!-- of the data is necessary and no transposition has to be done.
     
    629719!$OMP  PARALLEL PRIVATE ( i, j, k )
    630720!$OMP  DO
    631        !$acc kernels present( f_out )
     721       !$acc kernels present( f_inv, f_out )
    632722       !$acc loop
    633723       DO  k = 1, nz
     
    653743                          work(nys_x,1,nzb_x,0), sendrecvcount_zx, MPI_REAL, &
    654744                          comm1dx, ierr )
    655        !$acc update device( work )
    656745       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
    657746
     
    660749!$OMP  PARALLEL PRIVATE ( i, j, k, l, xs )
    661750!$OMP  DO
     751       !$acc data copyin( work )
    662752       DO  l = 0, pdims(1) - 1
    663753          xs = 0 + l * nnx
    664           !$acc kernels present( f_out, work )
     754          !$acc kernels present( f_out )
    665755          !$acc loop
    666756          DO  k = nzb_x, nzt_x
     
    674764          !$acc end kernels
    675765       ENDDO
     766       !$acc end data
    676767!$OMP  END PARALLEL
    677768#endif
     
    682773
    683774
    684  SUBROUTINE transpose_zy( f_in, work, f_out )
     775 SUBROUTINE resort_for_zy( f_inv, f_out )
     776
     777!------------------------------------------------------------------------------!
     778! Description:
     779! ------------
     780! Resorting data after the transposition from z to y. The transposition itself
     781! is carried out in transpose_zy
     782!------------------------------------------------------------------------------!
     783
     784     USE indices
     785     USE transpose_indices
     786
     787     IMPLICIT NONE
     788
     789     REAL ::  f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny)
     790     REAL ::  f_out(0:ny,nxl_y:nxr_y,nzb_y:nzt_y)
     791
     792
     793     INTEGER ::  i, j, k
     794
     795!
     796!-- Rearrange indices of input array in order to make data to be send
     797!-- by MPI contiguous
     798    !$OMP  PARALLEL PRIVATE ( i, j, k )
     799    !$OMP  DO
     800    !$acc kernels present( f_inv, f_out )
     801    !$acc loop
     802     DO  k = nzb_y, nzt_y
     803         DO  j = 0, ny
     804             !$acc loop vector( 32 )
     805             DO  i = nxl_y, nxr_y
     806                 f_out(j,i,k) = f_inv(i,k,j)
     807             ENDDO
     808         ENDDO
     809     ENDDO
     810     !$acc end kernels
     811     !$OMP  END PARALLEL
     812
     813 END SUBROUTINE resort_for_zy
     814
     815
     816 SUBROUTINE transpose_zy( f_in, f_inv )
    685817
    686818!------------------------------------------------------------------------------!
     
    702834    INTEGER ::  i, j, k, l, zs
    703835   
    704     REAL ::  f_in(nxl_z:nxr_z,nys_z:nyn_z,1:nz), f_out(0:ny,nxl_y:nxr_y,nzb_y:nzt_y)
     836    REAL ::  f_in(nxl_z:nxr_z,nys_z:nyn_z,1:nz), f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny)
    705837
    706838    REAL, DIMENSION(nxl_z:nxr_z,nzt_y-nzb_y+1,nys_z:nyn_z,0:pdims(1)-1) ::  work
    707 
    708     !$acc declare create( f_inv )
    709     REAL ::  f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny)
    710839
    711840
     
    720849!$OMP  PARALLEL PRIVATE ( i, j, k, l, zs )
    721850!$OMP  DO
     851       !$acc data copyout( work )
    722852       DO  l = 0, pdims(1) - 1
    723853          zs = 1 + l * ( nzt_y - nzb_y + 1 )
     
    734864          !$acc end kernels
    735865       ENDDO
     866       !$acc end data
    736867!$OMP  END PARALLEL
    737868
     
    740871       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
    741872       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    742        !$acc update host( work )
    743873       CALL MPI_ALLTOALL( work(nxl_z,1,nys_z,0), sendrecvcount_yz, MPI_REAL, &
    744874                          f_inv(nxl_y,nzb_y,0),  sendrecvcount_yz, MPI_REAL, &
     
    753883!$OMP  PARALLEL PRIVATE ( i, j, k )
    754884!$OMP  DO
    755        !$acc kernels present( f_in )
     885       !$acc kernels present( f_in, f_inv )
    756886       !$acc loop
    757887       DO  k = nzb_y, nzt_y
     
    768898    ENDIF
    769899
    770 !
    771 !-- Reorder transposed array in a way that the y index is in first position
    772 !$OMP  PARALLEL PRIVATE ( i, j, k )
    773 !$OMP  DO
    774     !$acc kernels present( f_out )
    775     !$acc loop
    776     DO  k = nzb_y, nzt_y
    777        DO  i = nxl_y, nxr_y
    778           !$acc loop vector( 32 )
    779           DO  j = 0, ny
    780              f_out(j,i,k) = f_inv(i,k,j)
    781           ENDDO
    782        ENDDO
    783     ENDDO
    784     !$acc end kernels
    785 !$OMP  END PARALLEL
    786 
    787900 END SUBROUTINE transpose_zy
    788901
    789902
    790  SUBROUTINE transpose_zyd( f_in, work, f_out )
     903 SUBROUTINE transpose_zyd( f_in, f_out )
    791904
    792905!------------------------------------------------------------------------------!
  • palm/trunk/SOURCE/tridia_solver.f90

    r1213 r1216  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! +tridia_substi_overlap for handling overlapping fft / transposition
    2323!
    2424! Former revisions:
     
    5757    END INTERFACE tridia_substi
    5858
    59     PUBLIC  tridia_substi, tridia_init, tridia_1dd
     59    INTERFACE tridia_substi_overlap
     60       MODULE PROCEDURE tridia_substi_overlap
     61    END INTERFACE tridia_substi_overlap
     62
     63    PUBLIC  tridia_substi, tridia_substi_overlap, tridia_init, tridia_1dd
    6064
    6165 CONTAINS
     
    260264
    261265
     266    SUBROUTINE tridia_substi_overlap( ar, jj )
     267
     268!------------------------------------------------------------------------------!
     269! Substitution (Forward and Backward) (Thomas algorithm)
     270!------------------------------------------------------------------------------!
     271
     272          USE arrays_3d,  ONLY: tri
     273          USE control_parameters
     274
     275          IMPLICIT NONE
     276
     277          INTEGER ::  i, j, jj, k
     278
     279          REAL    ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
     280
     281          !$acc declare create( ar1 )
     282          REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1)   ::  ar1
     283
     284!
     285!--       Forward substitution
     286          DO  k = 0, nz - 1
     287             !$acc kernels present( ar, tri )
     288             !$acc loop
     289             DO  j = nys_z, nyn_z
     290                DO  i = nxl_z, nxr_z
     291
     292                   IF ( k == 0 )  THEN
     293                      ar1(i,j,k) = ar(i,j,k+1)
     294                   ELSE
     295                      ar1(i,j,k) = ar(i,j,k+1) - tri(i,jj,k,2) * ar1(i,j,k-1)
     296                   ENDIF
     297
     298                ENDDO
     299             ENDDO
     300             !$acc end kernels
     301          ENDDO
     302
     303!
     304!--       Backward substitution
     305!--       Note, the 1.0E-20 in the denominator is due to avoid divisions
     306!--       by zero appearing if the pressure bc is set to neumann at the top of
     307!--       the model domain.
     308          DO  k = nz-1, 0, -1
     309             !$acc kernels present( ar, tri )
     310             !$acc loop
     311             DO  j = nys_z, nyn_z
     312                DO  i = nxl_z, nxr_z
     313
     314                   IF ( k == nz-1 )  THEN
     315                      ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20 )
     316                   ELSE
     317                      ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
     318                              / tri(i,jj,k,1)
     319                   ENDIF
     320                ENDDO
     321             ENDDO
     322             !$acc end kernels
     323          ENDDO
     324
     325!
     326!--       Indices i=0, j=0 correspond to horizontally averaged pressure.
     327!--       The respective values of ar should be zero at all k-levels if
     328!--       acceleration of horizontally averaged vertical velocity is zero.
     329          IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
     330             IF ( nys_z == 0  .AND.  nxl_z == 0 )  THEN
     331                !$acc kernels loop present( ar )
     332                DO  k = 1, nz
     333                   ar(nxl_z,nys_z,k) = 0.0
     334                ENDDO
     335             ENDIF
     336          ENDIF
     337
     338    END SUBROUTINE tridia_substi_overlap
     339
     340
    262341    SUBROUTINE split
    263342
Note: See TracChangeset for help on using the changeset viewer.