Changeset 1111


Ignore:
Timestamp:
Mar 8, 2013 11:54:10 PM (9 years ago)
Author:
raasch
Message:

New:
---

GPU porting of pres, swap_timelevel. Adjustments of openACC directives.
Further porting of poisfft, which now runs completely on GPU without any
host/device data transfer for serial an parallel runs (but parallel runs
require data transfer before and after the MPI transpositions).
GPU-porting of tridiagonal solver:
tridiagonal routines split into extermal subroutines (instead using CONTAINS),
no distinction between parallel/non-parallel in poisfft and tridia any more,
tridia routines moved to end of file because of probable bug in PGI compiler
(otherwise "invalid device function" is indicated during runtime).
(cuda_fft_interfaces, fft_xy, flow_statistics, init_3d_model, palm, poisfft, pres, prognostic_equations, swap_timelevel, time_integration, transpose)
output of accelerator board information. (header)

optimization of tridia routines: constant elements and coefficients of tri are
stored in seperate arrays ddzuw and tric, last dimension of tri reduced from 5 to 2,
(init_grid, init_3d_model, modules, palm, poisfft)

poisfft_init is now called internally from poisfft,
(Makefile, Makefile_check, init_pegrid, poisfft, poisfft_hybrid)

CPU-time per grid point and timestep is output to CPU_MEASURES file
(cpu_statistics, modules, time_integration)

Changed:


resorting from/to array work changed, work now has 4 dimensions instead of 1 (transpose)
array diss allocated only if required (init_3d_model)

pressure boundary condition "Neumann+inhomo" removed from the code
(check_parameters, header, poisfft, poisfft_hybrid, pres)

Errors:


bugfix: dependency added for cuda_fft_interfaces (Makefile)
bugfix: CUDA fft plans adjusted for domain decomposition (before they always
used total domain) (fft_xy)

Location:
palm/trunk
Files:
20 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SCRIPTS/.mrun.config.imuk_gpu

    • Property svn:keywords set to Id
    r1016 r1111  
     1#$Id$
    12#column 1          column 2                                   column 3
    23#name of variable  value of variable (~ must not be used)     scope
     
    89%add_source_path   $base_directory/USER_CODE/$fname
    910%depository_path   $base_directory/MAKE_DEPOSITORY
    10 #%use_makefile      true
    11 #
    12 # Enter your own host below by adding another line containing in the second
    13 # column your hostname (as provided by the unix command "hostname") and in the
    14 # third column the host identifier. Depending on your operating system, the
    15 # first characters of the host identifier should be "lc" (Linux cluster), "ibm"
    16 # (IBM-AIX), or "nec" (NEC-SX), respectively.
    1711#
    1812%host_identifier   inferno      lcmuk
    1913#
    20 # version 27/09/2012
    21 #
     14# pure MPI version
    2215%remote_username   <replace by your IMUK username>               lcmuk parallel pgi
    2316%tmp_user_catalog  /localdata                                    lcmuk parallel pgi
     
    2922%lopts             -Mcray=pointer:-fastsse:-r8                   lcmuk parallel pgi
    3023#
     24# pure MPI version with debug options
     25%remote_username   <replace by your IMUK username>               lcmuk parallel pgidbg
     26%tmp_user_catalog  /localdata                                    lcmuk parallel pgidbg
     27%compiler_name     mpif90                                        lcmuk parallel pgidbg
     28%compiler_name_ser pgf90                                         lcmuk parallel pgidbg
     29%cpp_options       -Mpreprocess:-DMPI_REAL=MPI_DOUBLE_PRECISION:-DMPI_2REAL=MPI_2DOUBLE_PRECISION:-D__nopointer   lcmuk parallel pgidbg
     30%mopts             -j:4                                          lcmuk parallel pgidbg
     31%fopts             -Mcray=pointer:-O0:-C:-g:-Mbounds:-Mchkstk:-traceback:-r8   lcmuk parallel pgidbg
     32%lopts             -Mcray=pointer:-O0:-C:-g:-Mbounds:-Mchkstk:-traceback:-r8   lcmuk parallel pgidbg
     33#
     34# pure GPU version
     35%remote_username   <replace by your IMUK username>                       lcmuk pgigpu
     36%tmp_user_catalog  /localdata                                            lcmuk pgigpu
     37%compiler_name     pgf90                                                 lcmuk pgigpu
     38%compiler_name_ser pgf90                                                 lcmuk pgigpu
     39%cpp_options       -Mpreprocess:-D__nopointer:-D__openacc:-D__cuda_fft   lcmuk pgigpu
     40%mopts             -j:4                                                  lcmuk pgigpu
     41%fopts             -acc:-ta=nvidia,4.1:-Minfo=acc:-Mcray=pointer:-fastsse:-r8:-Mcuda    lcmuk pgigpu
     42%lopts             -acc:-ta=nvidia,4.1:-Minfo=acc:-Mcray=pointer:-fastsse:-r8:-Mcuda:-L/localdata/opt/pgi/linux86-64/2012/cuda/4.1/lib64:-lcufft    lcmuk pgigpu
     43#
     44# MPI+GPU
    3145%remote_username   <replace by your IMUK username>               lcmuk parallel pgigpu
    3246%tmp_user_catalog  /localdata                                    lcmuk parallel pgigpu
    3347%compiler_name     mpif90                                        lcmuk parallel pgigpu
    3448%compiler_name_ser pgf90                                         lcmuk parallel pgigpu
    35 %cpp_options       -Mpreprocess:-DMPI_REAL=MPI_DOUBLE_PRECISION:-DMPI_2REAL=MPI_2DOUBLE_PRECISION:-D__nopointer:-D__openacc   lcmuk parallel pgigpu
     49%cpp_options       -Mpreprocess:-DMPI_REAL=MPI_DOUBLE_PRECISION:-DMPI_2REAL=MPI_2DOUBLE_PRECISION:-D__nopointer:-D__openacc:-D__cuda_fft   lcmuk parallel pgigpu
    3650%mopts             -j:4                                          lcmuk parallel pgigpu
    37 %fopts             -acc:-ta=nvidia,4.1:-Minfo=acc:-Mcray=pointer:-fastsse:-r8        lcmuk parallel pgigpu
    38 %lopts             -acc:-ta=nvidia,4.1:-Minfo=acc:-Mcray=pointer:-fastsse:-r8        lcmuk parallel pgigpu
    39 
     51%fopts             -acc:-ta=nvidia,4.1:-Minfo=acc:-Mcray=pointer:-fastsse:-r8:-Mcuda    lcmuk parallel pgigpu
     52%lopts             -acc:-ta=nvidia,4.1:-Minfo=acc:-Mcray=pointer:-fastsse:-r8:-Mcuda:-L/localdata/opt/pgi/linux86-64/2012/cuda/4.1/lib64:-lcufft   lcmuk parallel pgigpu
     53#
    4054%write_binary                true                             restart
    4155#
  • palm/trunk/SOURCE/Makefile

    r1107 r1111  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# dependencies removed from init_pegrid
     23# bugfix: dependency added for cuda_fft_interfaces
    2324#
    2425# Former revisions:
     
    267268cpu_log.o: modules.o
    268269cpu_statistics.o: modules.o
    269 cuda_fft_interfaces.o: cuda_fft_interfaces.f90
     270cuda_fft_interfaces.o: cuda_fft_interfaces.f90 modules.o
    270271data_log.o: modules.o
    271272data_output_dvrp.o: modules.o
     
    303304init_masks.o: modules.o
    304305init_ocean.o: modules.o eqn_state_seawater.o
    305 init_pegrid.o: modules.o fft_xy.o poisfft.o poisfft_hybrid.o
     306init_pegrid.o: modules.o
    306307init_pt_anomaly.o: modules.o
    307308init_rankine.o: modules.o
  • palm/trunk/SOURCE/Makefile_check

    r1107 r1111  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# dependencies removed from init_pegrid
    2323#
    2424# Former revisions:
     
    134134init_grid.o: modules.o
    135135init_masks.o: modules.o
    136 init_pegrid.o: modules.o fft_xy.o poisfft.o poisfft_hybrid.o
     136init_pegrid.o: modules.o
    137137local_stop.o: modules.o
    138138message.o: modules.o
  • palm/trunk/SOURCE/check_parameters.f90

    r1104 r1111  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! ibc_p_b = 2 removed
    2323!
    2424! Former revisions:
     
    16011601    ELSEIF ( bc_p_b == 'neumann' )  THEN
    16021602       ibc_p_b = 1
    1603     ELSEIF ( bc_p_b == 'neumann+inhomo' )  THEN
    1604        ibc_p_b = 2
    16051603    ELSE
    16061604       message_string = 'unknown boundary condition: bc_p_b = "' // &
     
    16081606       CALL message( 'check_parameters', 'PA0059', 1, 2, 0, 6, 0 )
    16091607    ENDIF
    1610     IF ( ibc_p_b == 2  .AND.  .NOT. prandtl_layer )  THEN
    1611        message_string = 'boundary condition: bc_p_b = "' // TRIM( bc_p_b ) // &
    1612                         '" not allowed with prandtl_layer = .FALSE.'
    1613        CALL message( 'check_parameters', 'PA0060', 1, 2, 0, 6, 0 )
    1614     ENDIF
     1608
    16151609    IF ( bc_p_t == 'dirichlet' )  THEN
    16161610       ibc_p_t = 0
  • palm/trunk/SOURCE/cpu_statistics.f90

    r1093 r1111  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! output of grid point numbers and average CPU time per grid point and timestep
    2323!
    2424! Former revisions:
     
    6969!------------------------------------------------------------------------------!
    7070
     71    USE control_parameters
    7172    USE cpulog
     73    USE indices,  ONLY: nx, ny, nz
    7274    USE pegrid
    73     USE control_parameters
    7475
    7576    IMPLICIT NONE
    7677
    7778    INTEGER    ::  i, ii(1), iii, sender
     79    REAL       ::  average_cputime
    7880    REAL, SAVE ::  norm = 1.0
    7981    REAL, DIMENSION(:),   ALLOCATABLE ::  pe_max, pe_min, pe_rms, sum
     
    157159
    158160!
     161!--    Get total time in order to calculate CPU-time per gridpoint and timestep
     162       IF ( nr_timesteps_this_run /= 0 )  THEN
     163          average_cputime = log_point(1)%sum / REAL( (nx+1) * (ny+1) * nz ) / &
     164                            REAL( nr_timesteps_this_run ) * 1E6  ! in micro-sec
     165       ELSE
     166          average_cputime = -1.0
     167       ENDIF
     168
     169!
    159170!--    Write cpu-times sorted by size
    160171       CALL check_open( 18 )
    161172#if defined( __parallel )
    162        WRITE ( 18, 100 )  TRIM( run_description_header ),        &
    163                           numprocs * threads_per_task, pdims(1), pdims(2), &
    164                           threads_per_task
     173       WRITE ( 18, 100 )  TRIM( run_description_header ),                          &
     174                          numprocs * threads_per_task, pdims(1), pdims(2),         &
     175                          threads_per_task, nx+1, ny+1, nz, nr_timesteps_this_run, &
     176                          average_cputime
     177                         
    165178       IF ( num_acc_per_node /= 0 )  WRITE ( 18, 108 )  num_acc_per_node
    166179       WRITE ( 18, 110 )
    167180#else
    168        WRITE ( 18, 100 )  TRIM( run_description_header ),        &
    169                           numprocs * threads_per_task, 1, 1, &
    170                           threads_per_task
     181       WRITE ( 18, 100 )  TRIM( run_description_header ),                          &
     182                          numprocs * threads_per_task, 1, 1,                       &
     183                          threads_per_task, nx+1, ny+1, nz, nr_timesteps_this_run, &
     184                          average_cputime
     185
    171186       IF ( num_acc_per_node /= 0 )  WRITE ( 18, 109 )  num_acc_per_node
    172187       WRITE ( 18, 110 )
     
    308323
    309324100 FORMAT (A/11('-')//'CPU measures for ',I5,' PEs (',I5,'(x) * ',I5,'(y', &
    310             &') tasks *',I5,' threads):')
     325            &') tasks *',I5,' threads):'//                                  &
     326            'gridpoints (x/y/z): ',20X,I5,' * ',I5,' * ',I5/                &
     327            'nr of timesteps: ',22X,I6/                                     &
     328            'cpu time per grid point and timestep: ',5X,F8.5,' * 10**-6 s')
    311329
    312330101 FORMAT (/'special measures:'/ &
     
    320338106 FORMAT (/'Exchange of ghostpoints via MPI_ISEND/MPI_IRECV')
    321339107 FORMAT (//)
    322 108 FORMAT ('Accelerator boards per node: ',I2)
    323 109 FORMAT ('Accelerator boards: ',I2)
     340108 FORMAT ('Accelerator boards per node: ',14X,I2)
     341109 FORMAT ('Accelerator boards: ',23X,I2)
    324342110 FORMAT ('----------------------------------------------------------',   &
    325343            &'------------'//&
  • palm/trunk/SOURCE/cuda_fft_interfaces.f90

    r1107 r1111  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! idata and odata changed from 1d- to 3d-arrays
    2323!
    2424! Former revisions:
     
    8787
    8888          INTEGER(C_INT), value ::  plan
    89           COMPLEX(dpk), device  ::  idata(*)
    90           REAL(dpk), device     ::  odata(*)
     89          COMPLEX(dpk), device  ::  idata(:,:,:)
     90          REAL(dpk), device     ::  odata(:,:,:)
    9191
    9292       END SUBROUTINE CUFFTEXECZ2D
     
    103103
    104104          INTEGER(C_INT), value ::  plan
    105           REAL(dpk), device     ::  idata(*)
    106           COMPLEX(dpk), device  ::  odata(*)
     105          REAL(dpk), device     ::  idata(:,:,:)
     106          COMPLEX(dpk), device  ::  odata(:,:,:)
    107107
    108108       END SUBROUTINE CUFFTEXECD2Z
  • palm/trunk/SOURCE/fft_xy.f90

    r1107 r1111  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! further openACC statements added, CUDA branch completely runs on GPU
     23! bugfix: CUDA fft plans adjusted for domain decomposition (before they always
     24! used total domain)
    2325!
    2426! Former revisions:
     
    213215          total_points_x_transpo = (nx+1) * (nyn_x-nys_x+1) * (nzt_x-nzb_x+1)
    214216          total_points_y_transpo = (ny+1) * (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1)
    215           CALL CUFFTPLAN1D( plan_xf, nx+1, CUFFT_D2Z, (ny+1)*nz )
    216           CALL CUFFTPLAN1D( plan_xi, nx+1, CUFFT_Z2D, (ny+1)*nz )
    217           CALL CUFFTPLAN1D( plan_yf, ny+1, CUFFT_D2Z, (nx+1)*nz )
    218           CALL CUFFTPLAN1D( plan_yi, ny+1, CUFFT_Z2D, (nx+1)*nz )
     217          CALL CUFFTPLAN1D( plan_xf, nx+1, CUFFT_D2Z, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
     218          CALL CUFFTPLAN1D( plan_xi, nx+1, CUFFT_Z2D, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
     219          CALL CUFFTPLAN1D( plan_yf, ny+1, CUFFT_D2Z, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
     220          CALL CUFFTPLAN1D( plan_yi, ny+1, CUFFT_Z2D, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
    219221#else
    220222          message_string = 'no system-specific fft-call available'
     
    261263
    262264       CHARACTER (LEN=*) ::  direction
    263        INTEGER ::  i, ishape(1), j, k, m
     265       INTEGER ::  i, ishape(1), j, k
    264266
    265267       LOGICAL ::  forward_fft
     
    273275       REAL, DIMENSION(6*(nx+1)) ::  work2
    274276#elif defined( __cuda_fft )
    275        REAL(dpk), DEVICE, DIMENSION(:), ALLOCATABLE    ::  cuda_a_device
    276        COMPLEX(dpk), DEVICE, DIMENSION(:), ALLOCATABLE ::  cuda_b_device
    277        COMPLEX(dpk), DIMENSION(:), ALLOCATABLE         ::  cuda_host
     277       !$acc declare create( ar_tmp )
     278       COMPLEX(dpk), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) ::  ar_tmp
    278279#endif
    279280       REAL, DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) ::  ar
     
    502503#elif defined( __cuda_fft )
    503504
    504           ALLOCATE( cuda_a_device(0:total_points_x_transpo-1) )
    505           ALLOCATE( cuda_b_device(0:((nx+1)/2+1) * (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) - 1) )
    506           ALLOCATE( cuda_host(0:((nx+1)/2+1) * (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) - 1) )
    507 
    508           m = 0
    509 
    510505          IF ( forward_fft )  THEN
    511506
    512              cuda_a_device = ar(0:total_points_x_transpo-1,nys_x,nzb_x)
    513 
    514              CALL CUFFTEXECD2Z( plan_xf, cuda_a_device, cuda_b_device )
    515              cuda_host = cuda_b_device
    516 
     507             !$acc data present( ar )
     508             CALL CUFFTEXECD2Z( plan_xf, ar, ar_tmp )
     509
     510             !$acc kernels
     511             !$acc loop
    517512             DO  k = nzb_x, nzt_x
    518513                DO  j = nys_x, nyn_x
    519514
     515                   !$acc loop vector( 32 )
    520516                   DO  i = 0, (nx+1)/2
    521                       ar(i,j,k)      = REAL( cuda_host(m+i) )  * dnx
    522                    ENDDO
    523 
     517                      ar(i,j,k)      = REAL( ar_tmp(i,j,k) )  * dnx
     518                   ENDDO
     519
     520                   !$acc loop vector( 32 )
    524521                   DO  i = 1, (nx+1)/2 - 1
    525                       ar(nx+1-i,j,k) = AIMAG( cuda_host(m+i) ) * dnx
    526                    ENDDO
    527 
    528                    m = m + (nx+1)/2 + 1
    529 
    530                 ENDDO
    531              ENDDO
    532 
    533           ELSE
    534 
     522                      ar(nx+1-i,j,k) = AIMAG( ar_tmp(i,j,k) ) * dnx
     523                   ENDDO
     524
     525                ENDDO
     526             ENDDO
     527             !$acc end kernels
     528             !$acc end data
     529
     530          ELSE
     531
     532             !$acc data present( ar )
     533             !$acc kernels
     534             !$acc loop
    535535             DO  k = nzb_x, nzt_x
    536536                DO  j = nys_x, nyn_x
    537537
    538                    cuda_host(m) = CMPLX( ar(0,j,k), 0.0 )
    539 
     538                   ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0 )
     539
     540                   !$acc loop vector( 32 )
    540541                   DO  i = 1, (nx+1)/2 - 1
    541                       cuda_host(m+i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
    542                    ENDDO
    543                    cuda_host(m+(nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
    544 
    545                    m = m + (nx+1)/2 + 1
    546 
    547                 ENDDO
    548              ENDDO
    549 
    550              cuda_b_device = cuda_host
    551              CALL CUFFTEXECZ2D( plan_xi, cuda_b_device, cuda_a_device )
    552 
    553              ar(0:total_points_x_transpo-1,nys_x,nzb_x) = cuda_a_device
    554 
    555           ENDIF
    556 
    557           DEALLOCATE( cuda_a_device, cuda_b_device, cuda_host )
     542                      ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
     543                   ENDDO
     544                   ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
     545
     546                ENDDO
     547             ENDDO
     548             !$acc end kernels
     549
     550             CALL CUFFTEXECZ2D( plan_xi, ar_tmp, ar )
     551             !$acc end data
     552
     553          ENDIF
    558554
    559555#else
     
    775771
    776772       CHARACTER (LEN=*) ::  direction
    777        INTEGER ::  i, j, jshape(1), k, m
     773       INTEGER ::  i, j, jshape(1), k
    778774
    779775       LOGICAL ::  forward_fft
     
    787783       REAL, DIMENSION(6*(ny+1)) ::  work2
    788784#elif defined( __cuda_fft )
    789        REAL(dpk), DEVICE, DIMENSION(:), ALLOCATABLE    ::  cuda_a_device
    790        COMPLEX(dpk), DEVICE, DIMENSION(:), ALLOCATABLE ::  cuda_b_device
    791        COMPLEX(dpk), DIMENSION(:), ALLOCATABLE         ::  cuda_host
     785       !$acc declare create( ar_tmp )
     786       COMPLEX(dpk), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) ::  ar_tmp
    792787#endif
    793788       REAL, DIMENSION(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) ::  ar
     
    10131008#elif defined( __cuda_fft )
    10141009
    1015           ALLOCATE( cuda_a_device(0:total_points_y_transpo-1) )
    1016           ALLOCATE( cuda_b_device(0:((ny+1)/2+1) * (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) - 1) )
    1017           ALLOCATE( cuda_host(0:((ny+1)/2+1) * (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) - 1) )
    1018 
    1019           m = 0
    1020 
    10211010          IF ( forward_fft )  THEN
    10221011
    1023              cuda_a_device = ar(0:total_points_y_transpo-1,nxl_y,nzb_y)
    1024 
    1025              CALL CUFFTEXECD2Z( plan_yf, cuda_a_device, cuda_b_device )
    1026              cuda_host = cuda_b_device
    1027 
     1012             !$acc data present( ar )
     1013             CALL CUFFTEXECD2Z( plan_yf, ar, ar_tmp )
     1014
     1015             !$acc kernels
     1016             !$acc loop
    10281017             DO  k = nzb_y, nzt_y
    10291018                DO  i = nxl_y, nxr_y
    10301019
     1020                   !$acc loop vector( 32 )
    10311021                   DO  j = 0, (ny+1)/2
    1032                       ar(j,i,k)      = REAL( cuda_host(m+j) )  * dny
    1033                    ENDDO
    1034 
     1022                      ar(j,i,k)      = REAL( ar_tmp(j,i,k) )  * dny
     1023                   ENDDO
     1024
     1025                   !$acc loop vector( 32 )
    10351026                   DO  j = 1, (ny+1)/2 - 1
    1036                       ar(ny+1-j,i,k) = AIMAG( cuda_host(m+j) ) * dny
    1037                    ENDDO
    1038 
    1039                    m = m + (ny+1)/2 + 1
    1040 
    1041                 ENDDO
    1042              ENDDO
    1043 
    1044           ELSE
    1045 
     1027                      ar(ny+1-j,i,k) = AIMAG( ar_tmp(j,i,k) ) * dny
     1028                   ENDDO
     1029
     1030                ENDDO
     1031             ENDDO
     1032             !$acc end kernels
     1033             !$acc end data
     1034
     1035          ELSE
     1036
     1037             !$acc data present( ar )
     1038             !$acc kernels
     1039             !$acc loop
    10461040             DO  k = nzb_y, nzt_y
    10471041                DO  i = nxl_y, nxr_y
    10481042
    1049                    cuda_host(m) = CMPLX( ar(0,i,k), 0.0 )
    1050 
     1043                   ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0 )
     1044
     1045                   !$acc loop vector( 32 )
    10511046                   DO  j = 1, (ny+1)/2 - 1
    1052                       cuda_host(m+j) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
    1053                    ENDDO
    1054                    cuda_host(m+(ny+1)/2) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
    1055 
    1056                    m = m + (ny+1)/2 + 1
    1057 
    1058                 ENDDO
    1059              ENDDO
    1060 
    1061              cuda_b_device = cuda_host
    1062              CALL CUFFTEXECZ2D( plan_yi, cuda_b_device, cuda_a_device )
    1063 
    1064              ar(0:total_points_y_transpo-1,nxl_y,nzb_y) = cuda_a_device
    1065 
    1066           ENDIF
    1067 
    1068           DEALLOCATE( cuda_a_device, cuda_b_device, cuda_host )
     1047                      ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
     1048                   ENDDO
     1049                   ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
     1050
     1051                ENDDO
     1052             ENDDO
     1053             !$acc end kernels
     1054
     1055             CALL CUFFTEXECZ2D( plan_yi, ar_tmp, ar )
     1056             !$acc end data
     1057
     1058          ENDIF
    10691059
    10701060#else
  • palm/trunk/SOURCE/flow_statistics.f90

    r1054 r1111  
    2020! Current revisions:
    2121! -----------------
     22! openACC directive added
    2223!
    2324! Former revisions:
     
    182183
    183184    ENDIF
     185
     186    !$acc update host( km, kh, e, pt, qs, qsws, rif, shf, ts, u, v, w )
    184187
    185188!
  • palm/trunk/SOURCE/header.f90

    r1109 r1111  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! output of accelerator board information
     23! ibc_p_b = 2 removed
    2324!
    2425! Former revisions:
     
    291292                          threads_per_task, pdims(1), pdims(2), TRIM( char1 )
    292293    ENDIF
     294    IF ( num_acc_per_node /= 0 )  WRITE ( io, 117 )  num_acc_per_node   
    293295    IF ( ( host(1:3) == 'ibm'  .OR.  host(1:3) == 'nec'  .OR.    &
    294296           host(1:2) == 'lc'   .OR.  host(1:3) == 'dec' )  .AND. &
     
    305307       WRITE ( io, 108 )  maximum_parallel_io_streams
    306308    ENDIF
     309#else
     310    IF ( num_acc_per_node /= 0 )  WRITE ( io, 120 )  num_acc_per_node
    307311#endif
    308312    WRITE ( io, 99 )
     
    593597    ELSEIF ( ibc_p_b == 1 )  THEN
    594598       runten = 'p(0)     = p(1)   |'
    595     ELSE
    596        runten = 'p(0)     = p(1) +R|'
    597599    ENDIF
    598600    IF ( ibc_p_t == 0 )  THEN
     
    16131615            37X,'independent precursor runs'/             &
    16141616            37X,42('-'))
     1617117 FORMAT (' Accelerator boards / node:  ',I2)
    16151618#endif
    16161619110 FORMAT (/' Numerical Schemes:'/ &
     
    16271630            '     translation velocity = ',A/ &
    16281631            '     distance advected ',A,':  ',F8.3,' km(x)  ',F8.3,' km(y)')
     1632120 FORMAT (' Accelerator boards: ',8X,I2)
    16291633122 FORMAT (' --> Time differencing scheme: ',A)
    16301634123 FORMAT (' --> Rayleigh-Damping active, starts ',A,' z = ',F8.2,' m'/ &
     
    16801684             ' CPU-time used:       ',F9.3,' s     per timestep:               ', &
    16811685               '  ',F9.3,' s'/                                                    &
    1682              '                                   per second of simulated tim',    &
     1686             '                                      per second of simulated tim', &
    16831687               'e: ',F9.3,' s')
    16841688207 FORMAT ( ' Coupling start time: ',F9.3,' s')
  • palm/trunk/SOURCE/init_3d_model.f90

    r1093 r1111  
    2323! Current revisions:
    2424! ------------------
    25 !
     25! openACC directives added for pres
     26! array diss allocated only if required
    2627!
    2728! Former revisions:
     
    234235    USE random_function_mod
    235236    USE statistics
     237    USE transpose_indices
    236238
    237239    IMPLICIT NONE
     
    334336    ENDIF
    335337
     338!
     339!-- Array for storing constant coeffficients of the tridiagonal solver
     340    IF ( psolver == 'poisfft' )  THEN
     341       ALLOCATE( tric(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) )
     342    ENDIF
     343
    336344    IF ( humidity  .OR.  passive_scalar ) THEN
    337345!
     
    474482    IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence )  THEN
    475483       ALLOCATE ( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    476     ELSE
    477        ALLOCATE ( diss(2,2,2) )  ! required because diss is used as a
    478                                  ! formal parameter
    479484    ENDIF
    480485
     
    14271432       CALL disturb_field( nzb_v_inner, tend, v )
    14281433       n_sor = nsor_ini
     1434       !$acc data copy( d, ddzw, nzb_s_inner, tric, u, v, w, tend )
    14291435       CALL pres
     1436       !$acc end data
    14301437       n_sor = nsor
    14311438    ENDIF
  • palm/trunk/SOURCE/init_pegrid.f90

    r1093 r1111  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! initialization of poisfft moved to poisfft
    2323!
    2424! Former revisions:
     
    150150
    151151    USE control_parameters
    152     USE fft_xy
    153152    USE grid_variables
    154153    USE indices
    155154    USE pegrid
    156     USE poisfft_mod
    157     USE poisfft_hybrid_mod
    158155    USE statistics
    159156    USE transpose_indices
     
    11351132    ENDIF
    11361133
    1137     IF ( psolver == 'poisfft_hybrid' )  THEN
    1138        CALL poisfft_hybrid_ini
    1139     ELSEIF ( psolver == 'poisfft' )  THEN
    1140        CALL poisfft_init
    1141     ENDIF
    1142 
    11431134!
    11441135!-- Allocate wall flag arrays used in the multigrid solver
  • palm/trunk/SOURCE/modules.f90

    r1107 r1111  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! +tric, nr_timesteps_this_run
    2323!
    2424! Former revisions:
     
    407407
    408408    REAL, DIMENSION(:,:), ALLOCATABLE ::                                       &
    409           c_u, c_v, c_w, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q, diss_s_qr,  &
    410           diss_s_sa, diss_s_u, diss_s_v, diss_s_w, dzu_mg, dzw_mg, flux_s_e,   &
    411           flux_s_nr, flux_s_pt, flux_s_q, flux_s_qr, flux_s_sa, flux_s_u,      &
    412           flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg, mean_inflow_profiles, nrs,  &
    413           nrsws, nrswst, pt_slope_ref, qs, qsws, qswst, qswst_remote, qrs,     &
    414           qrsws, qrswst, rif, saswsb, saswst, shf, total_2d_a, total_2d_o, ts, &
    415           tswst, us, usws, uswst, vsws, vswst, z0, z0h
    416          
     409          c_u, c_v, c_w, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q,             &
     410          diss_s_qr, diss_s_sa, diss_s_u, diss_s_v, diss_s_w, dzu_mg, dzw_mg,  &
     411          flux_s_e, flux_s_nr, flux_s_pt, flux_s_q, flux_s_qr, flux_s_sa,      &
     412          flux_s_u, flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg,                   &
     413          mean_inflow_profiles, nrs, nrsws, nrswst, pt_slope_ref, qs, qsws,    &
     414          qswst, qswst_remote, qrs, qrsws, qrswst, rif, saswsb, saswst, shf,  &
     415          total_2d_a, total_2d_o, ts, tswst, us, usws, uswst, vsws, vswst, z0, &
     416          z0h
    417417
    418418    REAL, DIMENSION(:,:,:), ALLOCATABLE ::                                     &
     
    422422          flux_l_qr, flux_l_sa, flux_l_u, flux_l_v, flux_l_w, kh, km, lad_s,   &
    423423          lad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls, tend, tend_pt,    &
    424           tend_nr, tend_q, tend_qr, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n,  &
    425           v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s
     424          tend_nr, tend_q, tend_qr, tric, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l,   &
     425          v_m_n, v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s
    426426           
    427427
     
    684684                maximum_parallel_io_streams = -1, max_pr_user = 0, &
    685685                mgcycles = 0, mg_cycles = -1, mg_switch_to_pe0_level = 0, mid, &
    686                 netcdf_data_format = 2, ngsrb = 2, nsor = 20, &
    687                 nsor_ini = 100, n_sor, normalizing_region = 0, &
     686                netcdf_data_format = 2, ngsrb = 2, nr_timesteps_this_run = 0, &
     687                nsor = 20, nsor_ini = 100, n_sor, normalizing_region = 0, &
    688688                nz_do3d = -9999, pch_index = 0, prt_time_count = 0, &
    689689                recycling_plane, runnr = 0, &
  • palm/trunk/SOURCE/palm.f90

    r1093 r1111  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! openACC statements updated
    2323!
    2424! Former revisions:
     
    234234!-- Declare and initialize variables in the accelerator memory with their
    235235!-- host values
    236     !$acc  data copyin( diss, e, e_p, kh, km, pt, pt_p, q, ql, tend, te_m, tpt_m, tu_m, tv_m, tw_m, u, u_p, v, vpt, v_p, w, w_p )          &
    237     !$acc       copyin( ddzu, ddzw, dd2zu, l_grid, l_wall, ptdf_x, ptdf_y, pt_init, rdf, rdf_sc, ug, vg, zu, zw )   &
     236    !$acc  data copyin( d, diss, e, e_p, kh, km, pt, pt_p, q, ql, tend, te_m, tpt_m, tu_m, tv_m, tw_m, u, u_p, v, vpt, v_p, w, w_p )          &
     237    !$acc       copyin( tric, ddzu, ddzw, dd2zu, l_grid, l_wall, ptdf_x, ptdf_y, pt_init, rdf, rdf_sc, ug, vg, zu, zw )   &
    238238    !$acc       copyin( hom, qs, qsws, qswst, rif, rif_wall, shf, ts, tswst, us, usws, uswst, vsws, vswst, z0, z0h )      &
    239239    !$acc       copyin( fxm, fxp, fym, fyp, fwxm, fwxp, fwym, fwyp, nzb_diff_s_inner, nzb_diff_s_outer, nzb_diff_u )       &
  • palm/trunk/SOURCE/poisfft.f90

    r1107 r1111  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! further openACC porting of non-parallel (MPI) branch:
     23! tridiagonal routines split into extermal subroutines (instead using CONTAINS),
     24! no distinction between parallel/non-parallel in poisfft and tridia any more,
     25! tridia routines moved to end of file because of probable bug in PGI compiler
     26! (otherwise "invalid device function" is indicated during runtime),
     27! optimization of tridia routines: constant elements and coefficients of tri are
     28! stored in seperate arrays ddzuw and tric, last dimension of tri reduced from 5
     29! to 2,
     30! poisfft_init is now called internally from poisfft, maketri is called from
     31! poisfft_init,
     32! ibc_p_b = 2 removed
    2333!
    2434! Former revisions:
     
    157167    IMPLICIT NONE
    158168
     169    LOGICAL, SAVE ::  poisfft_initialized = .FALSE.
     170
     171    REAL, DIMENSION(:,:), ALLOCATABLE ::  ddzuw
     172
    159173    PRIVATE
    160174
     
    181195    SUBROUTINE poisfft_init
    182196
     197       USE arrays_3d,  ONLY:  ddzu_pres, ddzw
     198
     199       IMPLICIT NONE
     200
     201       INTEGER ::  k
     202
     203
    183204       CALL fft_init
    184205
     206       ALLOCATE( ddzuw(0:nz-1,3) )
     207
     208       DO  k = 0, nz-1
     209          ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1)
     210          ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1)
     211          ddzuw(k,3) = -1.0 * &
     212                       ( ddzu_pres(k+2) * ddzw(k+1) + ddzu_pres(k+1) * ddzw(k+1) )
     213       ENDDO
     214!
     215!--    Calculate constant coefficients of the tridiagonal matrix
     216#if ! defined ( __check )
     217       CALL maketri
     218#endif
     219
     220       poisfft_initialized = .TRUE.
     221
    185222    END SUBROUTINE poisfft_init
     223
    186224
    187225#if ! defined ( __check )
     
    199237       CALL cpu_log( log_point_s(3), 'poisfft', 'start' )
    200238
     239       IF ( .NOT. poisfft_initialized )  CALL poisfft_init
     240
    201241!
    202242!--    Two-dimensional Fourier Transformation in x- and y-direction.
    203 #if defined( __parallel )
    204        IF ( pdims(2) == 1 )  THEN
     243       IF ( pdims(2) == 1  .AND.  pdims(1) > 1 )  THEN
    205244
    206245!
     
    217256          CALL tr_xy_ffty( ar, work, ar )
    218257
    219        ELSEIF ( pdims(1) == 1 )  THEN
     258       ELSEIF ( pdims(1) == 1  .AND.  pdims(2) > 1 )  THEN
    220259
    221260!
     
    235274
    236275!
    237 !--       2d-domain-decomposition
     276!--       2d-domain-decomposition or no decomposition (1 PE run)
    238277!--       Transposition z --> x
    239278          CALL cpu_log( log_point_s(5), 'transpo forward', 'start' )
     
    296335       ENDIF
    297336
    298 #else
    299 
    300 !
    301 !--    Two-dimensional Fourier Transformation along x- and y-direction.
    302        CALL cpu_log( log_point_s(5), 'transpo forward', 'start' )
    303        !$acc data copyin( ar, work )
    304        CALL transpose_zx( ar, work, ar )
    305        !$acc update host( ar )
    306        CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
    307        
    308 
    309        CALL cpu_log( log_point_s(4), 'fft_x', 'start' )
    310        CALL fft_x( ar, 'forward' )
    311        CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
    312 
    313        CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
    314        CALL transpose_xy( ar, work, ar )
    315        CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
    316 
    317        CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
    318        CALL fft_y( ar, 'forward' )
    319        CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
    320 
    321 !
    322 !--    Solve the tridiagonal equation system along z
    323        CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
    324        CALL transpose_yz( ar, work, ar )
    325        CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
    326 
    327        CALL cpu_log( log_point_s(6), 'tridia', 'start' )
    328        CALL tridia( ar )
    329        CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
    330 
    331        CALL cpu_log( log_point_s(8), 'transpo invers', 'start' )
    332        CALL transpose_zy( ar, work, ar )
    333        CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
    334 
    335 !
    336 !--    Inverse Fourier Transformation.
    337        CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
    338        CALL fft_y( ar, 'backward' )
    339        CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
    340 
    341        CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
    342        CALL transpose_yx( ar, work, ar )
    343        CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
    344 
    345        CALL cpu_log( log_point_s(4), 'fft_x', 'continue' )
    346        CALL fft_x( ar, 'backward' )
    347        CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
    348 
    349        CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
    350        CALL transpose_xz( ar, work, ar )
    351        CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
    352 
    353        !$acc end data
    354 
    355 #endif
    356 
    357337       CALL cpu_log( log_point_s(3), 'poisfft', 'stop' )
    358338
     
    361341
    362342
    363     SUBROUTINE tridia( ar )
    364 
    365 !------------------------------------------------------------------------------!
    366 ! solves the linear system of equations:
    367 !
    368 ! -(4 pi^2(i^2/(dx^2*nnx^2)+j^2/(dy^2*nny^2))+
    369 !   1/(dzu(k)*dzw(k))+1/(dzu(k-1)*dzw(k)))*p(i,j,k)+
    370 ! 1/(dzu(k)*dzw(k))*p(i,j,k+1)+1/(dzu(k-1)*dzw(k))*p(i,j,k-1)=d(i,j,k)
    371 !
    372 ! by using the Thomas algorithm
    373 !------------------------------------------------------------------------------!
    374 
    375        USE arrays_3d
    376 
    377        IMPLICIT NONE
    378 
    379        INTEGER ::  i, j, k, nnyh
    380 
    381        REAL, DIMENSION(nxl_z:nxr_z,0:nz-1)   ::  ar1
    382        REAL, DIMENSION(5,nxl_z:nxr_z,0:nz-1) ::  tri
    383 
    384        REAL    ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
    385 
    386 
    387        nnyh = (ny+1) / 2
    388 
    389 !
    390 !--    Define constant elements of the tridiagonal matrix.
    391 !$OMP  PARALLEL PRIVATE ( k, i )
    392 !$OMP  DO
    393        DO  k = 0, nz-1
    394           DO  i = nxl_z, nxr_z
    395              tri(2,i,k) = ddzu_pres(k+1) * ddzw(k+1)
    396              tri(3,i,k) = ddzu_pres(k+2) * ddzw(k+1)
    397           ENDDO
    398        ENDDO
    399 !$OMP  END PARALLEL
    400 
    401 #if defined( __parallel )
    402 !
    403 !--    Repeat for all y-levels.
    404 !$OMP  PARALLEL FIRSTPRIVATE( tri ) PRIVATE ( ar1, j )
    405 !$OMP  DO
    406        DO  j = nys_z, nyn_z
    407           IF ( j <= nnyh )  THEN
    408              CALL maketri( j )
    409           ELSE
    410              CALL maketri( ny+1-j )
    411           ENDIF
    412           CALL split
    413           CALL substi( j )
    414        ENDDO
    415 !$OMP  END PARALLEL
    416 #else
    417 !
    418 !--    First y-level.
    419        CALL maketri( nys_z )
    420        CALL split
    421        CALL substi( 0 )
    422 
    423 !
    424 !--    Further y-levels.
    425        DO  j = 1, nnyh - 1
    426           CALL maketri( j )
    427           CALL split
    428           CALL substi( j )
    429           CALL substi( ny+1-j )
    430        ENDDO
    431        CALL maketri( nnyh )
    432        CALL split
    433        CALL substi( nnyh+nys )
    434 #endif
    435 
    436     CONTAINS
    437 
    438        SUBROUTINE maketri( j )
    439 
    440 !------------------------------------------------------------------------------!
    441 ! Computes the i- and j-dependent component of the matrix
    442 !------------------------------------------------------------------------------!
    443 
    444           USE arrays_3d
    445           USE constants
    446           USE control_parameters
    447           USE grid_variables
    448 
    449           IMPLICIT NONE
    450 
    451           INTEGER ::  i, j, k, nnxh
    452           REAL    ::  a, c
    453           REAL    ::  ll(nxl_z:nxr_z)
    454 
    455 
    456           nnxh = ( nx + 1 ) / 2
    457 
    458 !
    459 !--       Provide the tridiagonal matrix for solution of the Poisson equation in
    460 !--       Fourier space. The coefficients are computed following the method of
    461 !--       Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
    462 !--       Siano's original version by discretizing the Poisson equation,
    463 !--       before it is Fourier-transformed
    464 #if defined( __parallel )
    465           DO  i = nxl_z, nxr_z
    466              IF ( i >= 0 .AND. i <= nnxh )  THEN
    467                 ll(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
    468                                           REAL( nx+1 ) ) ) / ( dx * dx ) + &
    469                         2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
    470                                           REAL( ny+1 ) ) ) / ( dy * dy )
    471              ELSE
    472                 ll(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
    473                                           REAL( nx+1 ) ) ) / ( dx * dx ) + &
    474                         2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
    475                                           REAL( ny+1 ) ) ) / ( dy * dy )
    476              ENDIF
    477              DO  k = 0,nz-1
    478                 a = -1.0 * ddzu_pres(k+2) * ddzw(k+1)
    479                 c = -1.0 * ddzu_pres(k+1) * ddzw(k+1)
    480                 tri(1,i,k) = a + c - ll(i)
    481              ENDDO
    482           ENDDO
    483 #else
    484           DO  i = 0, nnxh
    485              ll(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / REAL( nx+1 ) ) ) / &
    486                            ( dx * dx ) + &
    487                      2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / REAL( ny+1 ) ) ) / &
    488                            ( dy * dy )
    489              DO  k = 0, nz-1
    490                 a = -1.0 * ddzu_pres(k+2) * ddzw(k+1)
    491                 c = -1.0 * ddzu_pres(k+1) * ddzw(k+1)
    492                 tri(1,i,k) = a + c - ll(i)
    493                 IF ( i >= 1 .and. i < nnxh )  THEN
    494                    tri(1,nx+1-i,k) = tri(1,i,k)
    495                 ENDIF
    496              ENDDO
    497           ENDDO
    498 #endif
    499           IF ( ibc_p_b == 1  .OR.  ibc_p_b == 2 )  THEN
    500              DO  i = nxl_z, nxr_z
    501                 tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
    502              ENDDO
    503           ENDIF
    504           IF ( ibc_p_t == 1 )  THEN
    505              DO  i = nxl_z, nxr_z
    506                 tri(1,i,nz-1) = tri(1,i,nz-1) + tri(3,i,nz-1)
    507              ENDDO
    508           ENDIF
    509 
    510        END SUBROUTINE maketri
    511 
    512 
    513        SUBROUTINE substi( j )
    514 
    515 !------------------------------------------------------------------------------!
    516 ! Substitution (Forward and Backward) (Thomas algorithm)
    517 !------------------------------------------------------------------------------!
    518 
    519           USE control_parameters
    520 
    521           IMPLICIT NONE
    522 
    523           INTEGER ::  i, j, k
    524 
    525 !
    526 !--       Forward substitution.
    527           DO  i = nxl_z, nxr_z
    528              ar1(i,0) = ar(i,j,1)
    529           ENDDO
    530           DO  k = 1, nz - 1
    531              DO  i = nxl_z, nxr_z
    532                 ar1(i,k) = ar(i,j,k+1) - tri(5,i,k) * ar1(i,k-1)
    533              ENDDO
    534           ENDDO
    535 
    536 !
    537 !--       Backward substitution
    538 !--       Note, the 1.0E-20 in the denominator is due to avoid divisions
    539 !--       by zero appearing if the pressure bc is set to neumann at the top of
    540 !--       the model domain.
    541           DO  i = nxl_z, nxr_z
    542              ar(i,j,nz) = ar1(i,nz-1) / ( tri(4,i,nz-1) + 1.0E-20 )
    543           ENDDO
    544           DO  k = nz-2, 0, -1
    545              DO  i = nxl_z, nxr_z
    546                 ar(i,j,k+1) = ( ar1(i,k) - tri(3,i,k) * ar(i,j,k+2) ) &
    547                               / tri(4,i,k)
    548              ENDDO
    549           ENDDO
    550 
    551 !
    552 !--       Indices i=0, j=0 correspond to horizontally averaged pressure.
    553 !--       The respective values of ar should be zero at all k-levels if
    554 !--       acceleration of horizontally averaged vertical velocity is zero.
    555           IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
    556              IF ( j == 0  .AND.  nxl_z == 0 )  THEN
    557                 DO  k = 1, nz
    558                    ar(nxl_z,j,k) = 0.0
    559                 ENDDO
    560              ENDIF
    561           ENDIF
    562 
    563        END SUBROUTINE substi
    564 
    565 
    566        SUBROUTINE split
    567 
    568 !------------------------------------------------------------------------------!
    569 ! Splitting of the tridiagonal matrix (Thomas algorithm)
    570 !------------------------------------------------------------------------------!
    571 
    572           IMPLICIT NONE
    573 
    574           INTEGER ::  i, k
    575 
    576 !
    577 !--       Splitting.
    578           DO  i = nxl_z, nxr_z
    579              tri(4,i,0) = tri(1,i,0)
    580           ENDDO
    581           DO  k = 1, nz-1
    582              DO  i = nxl_z, nxr_z
    583                 tri(5,i,k) = tri(2,i,k) / tri(4,i,k-1)
    584                 tri(4,i,k) = tri(1,i,k) - tri(3,i,k-1) * tri(5,i,k)
    585              ENDDO
    586           ENDDO
    587 
    588        END SUBROUTINE split
    589 
    590     END SUBROUTINE tridia
    591 
    592 
    593 #if defined( __parallel )
    594343    SUBROUTINE ffty_tr_yx( f_in, work, f_out )
    595344
     
    706455!
    707456!--    Transpose array
     457#if defined( __parallel )
    708458       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
    709459       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     
    712462                          comm1dx, ierr )
    713463       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
     464#endif
    714465
    715466    END SUBROUTINE ffty_tr_yx
     
    745496!
    746497!--    Transpose array
     498#if defined( __parallel )
    747499       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
    748500       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     
    751503                          comm1dx, ierr )
    752504       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
     505#endif
    753506
    754507!
     
    1061814!
    1062815!--    Transpose array
     816#if defined( __parallel )
    1063817       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
    1064818       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     
    1067821                          comm1dy, ierr )
    1068822       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
     823#endif
    1069824
    1070825    END SUBROUTINE fftx_tr_xy
     
    1096851!
    1097852!--    Transpose array
     853#if defined( __parallel )
    1098854       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
    1099855       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     
    1102858                          comm1dy, ierr )
    1103859       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
     860#endif
    1104861
    1105862!
     
    14261183             ENDDO
    14271184          ENDDO
    1428           IF ( ibc_p_b == 1  .OR.  ibc_p_b == 2 )  THEN
     1185          IF ( ibc_p_b == 1 )  THEN
    14291186             DO  i = 0, nx
    14301187                tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
     
    15301287    END SUBROUTINE tridia_1dd
    15311288
    1532 #endif
    1533 #endif
     1289
     1290    SUBROUTINE tridia( ar )
     1291
     1292!------------------------------------------------------------------------------!
     1293! solves the linear system of equations:
     1294!
     1295! -(4 pi^2(i^2/(dx^2*nnx^2)+j^2/(dy^2*nny^2))+
     1296!   1/(dzu(k)*dzw(k))+1/(dzu(k-1)*dzw(k)))*p(i,j,k)+
     1297! 1/(dzu(k)*dzw(k))*p(i,j,k+1)+1/(dzu(k-1)*dzw(k))*p(i,j,k-1)=d(i,j,k)
     1298!
     1299! by using the Thomas algorithm
     1300!------------------------------------------------------------------------------!
     1301
     1302       USE arrays_3d
     1303
     1304       IMPLICIT NONE
     1305
     1306       INTEGER ::  i, j, k
     1307
     1308       !$acc declare create( tri )
     1309       REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) ::  tri
     1310
     1311       REAL    ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
     1312
     1313
     1314       CALL split( tri )
     1315       CALL substi( ar, tri )
     1316
     1317    END SUBROUTINE tridia
     1318
     1319
     1320    SUBROUTINE maketri
     1321
     1322!------------------------------------------------------------------------------!
     1323! Computes the i- and j-dependent component of the matrix
     1324!------------------------------------------------------------------------------!
     1325
     1326          USE arrays_3d,  ONLY: tric
     1327          USE constants
     1328          USE control_parameters
     1329          USE grid_variables
     1330
     1331          IMPLICIT NONE
     1332
     1333          INTEGER ::  i, j, k, nnxh, nnyh
     1334
     1335          !$acc declare create( ll )
     1336          REAL    ::  ll(nxl_z:nxr_z,nys_z:nyn_z)
     1337
     1338
     1339          nnxh = ( nx + 1 ) / 2
     1340          nnyh = ( ny + 1 ) / 2
     1341
     1342!
     1343!--       Provide the constant coefficients of the tridiagonal matrix for solution
     1344!--       of the Poisson equation in Fourier space.
     1345!--       The coefficients are computed following the method of
     1346!--       Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
     1347!--       Siano's original version by discretizing the Poisson equation,
     1348!--       before it is Fourier-transformed.
     1349
     1350          !$acc kernels present( tric )
     1351          !$acc loop vector( 32 )
     1352          DO  j = nys_z, nyn_z
     1353             DO  i = nxl_z, nxr_z
     1354                IF ( j >= 0  .AND.  j <= nnyh )  THEN
     1355                   IF ( i >= 0  .AND.  i <= nnxh )  THEN
     1356                      ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
     1357                                                  REAL( nx+1 ) ) ) / ( dx * dx ) + &
     1358                                2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
     1359                                                  REAL( ny+1 ) ) ) / ( dy * dy )
     1360                   ELSE
     1361                      ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
     1362                                                  REAL( nx+1 ) ) ) / ( dx * dx ) + &
     1363                                2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
     1364                                                  REAL( ny+1 ) ) ) / ( dy * dy )
     1365                   ENDIF
     1366                ELSE
     1367                   IF ( i >= 0  .AND.  i <= nnxh )  THEN
     1368                      ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
     1369                                                  REAL( nx+1 ) ) ) / ( dx * dx ) + &
     1370                                2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
     1371                                                  REAL( ny+1 ) ) ) / ( dy * dy )
     1372                   ELSE
     1373                      ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
     1374                                                  REAL( nx+1 ) ) ) / ( dx * dx ) + &
     1375                                2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
     1376                                                  REAL( ny+1 ) ) ) / ( dy * dy )
     1377                   ENDIF
     1378                ENDIF
     1379             ENDDO
     1380          ENDDO
     1381
     1382          !$acc loop
     1383          DO  k = 0, nz-1
     1384             DO  j = nys_z, nyn_z
     1385                !$acc loop vector( 32 )
     1386                DO  i = nxl_z, nxr_z
     1387                   tric(i,j,k) = ddzuw(k,3) - ll(i,j)
     1388                ENDDO
     1389             ENDDO
     1390          ENDDO
     1391          !$acc end kernels
     1392
     1393          IF ( ibc_p_b == 1 )  THEN
     1394             !$acc kernels present( tric )
     1395             !$acc loop
     1396             DO  j = nys_z, nyn_z
     1397                DO  i = nxl_z, nxr_z
     1398                   tric(i,j,0) = tric(i,j,0) + ddzuw(0,1)
     1399                ENDDO
     1400             ENDDO
     1401             !$acc end kernels
     1402          ENDIF
     1403          IF ( ibc_p_t == 1 )  THEN
     1404             !$acc kernels present( tric )
     1405             !$acc loop
     1406             DO  j = nys_z, nyn_z
     1407                DO  i = nxl_z, nxr_z
     1408                   tric(i,j,nz-1) = tric(i,j,nz-1) + ddzuw(nz-1,2)
     1409                ENDDO
     1410             ENDDO
     1411             !$acc end kernels
     1412          ENDIF
     1413
     1414    END SUBROUTINE maketri
     1415
     1416
     1417    SUBROUTINE substi( ar, tri )
     1418
     1419!------------------------------------------------------------------------------!
     1420! Substitution (Forward and Backward) (Thomas algorithm)
     1421!------------------------------------------------------------------------------!
     1422
     1423          USE control_parameters
     1424
     1425          IMPLICIT NONE
     1426
     1427          INTEGER ::  i, j, k
     1428
     1429          REAL    ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
     1430          REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) ::  tri
     1431
     1432          !$acc declare create( ar1 )
     1433          REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1)   ::  ar1
     1434
     1435!
     1436!--       Forward substitution
     1437          DO  k = 0, nz - 1
     1438             !$acc kernels present( ar, tri )
     1439             !$acc loop
     1440             DO  j = nys_z, nyn_z
     1441                DO  i = nxl_z, nxr_z
     1442
     1443                   IF ( k == 0 )  THEN
     1444                      ar1(i,j,k) = ar(i,j,k+1)
     1445                   ELSE
     1446                      ar1(i,j,k) = ar(i,j,k+1) - tri(i,j,k,2) * ar1(i,j,k-1)
     1447                   ENDIF
     1448
     1449                ENDDO
     1450             ENDDO
     1451             !$acc end kernels
     1452          ENDDO
     1453
     1454!
     1455!--       Backward substitution
     1456!--       Note, the 1.0E-20 in the denominator is due to avoid divisions
     1457!--       by zero appearing if the pressure bc is set to neumann at the top of
     1458!--       the model domain.
     1459          DO  k = nz-1, 0, -1
     1460             !$acc kernels present( ar, tri )
     1461             !$acc loop
     1462             DO  j = nys_z, nyn_z
     1463                DO  i = nxl_z, nxr_z
     1464
     1465                   IF ( k == nz-1 )  THEN
     1466                      ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20 )
     1467                   ELSE
     1468                      ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
     1469                              / tri(i,j,k,1)
     1470                   ENDIF
     1471                ENDDO
     1472             ENDDO
     1473             !$acc end kernels
     1474          ENDDO
     1475
     1476!
     1477!--       Indices i=0, j=0 correspond to horizontally averaged pressure.
     1478!--       The respective values of ar should be zero at all k-levels if
     1479!--       acceleration of horizontally averaged vertical velocity is zero.
     1480          IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
     1481             IF ( nys_z == 0  .AND.  nxl_z == 0 )  THEN
     1482                !$acc kernels loop present( ar )
     1483                DO  k = 1, nz
     1484                   ar(nxl_z,nys_z,k) = 0.0
     1485                ENDDO
     1486             ENDIF
     1487          ENDIF
     1488
     1489    END SUBROUTINE substi
     1490
     1491
     1492    SUBROUTINE split( tri )
     1493
     1494!------------------------------------------------------------------------------!
     1495! Splitting of the tridiagonal matrix (Thomas algorithm)
     1496!------------------------------------------------------------------------------!
     1497
     1498          USE arrays_3d,  ONLY: tric
     1499
     1500          IMPLICIT NONE
     1501
     1502          INTEGER ::  i, j, k
     1503
     1504          REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) ::  tri
     1505
     1506!
     1507!--       Splitting
     1508          !$acc kernels present( tri, tric )
     1509          !$acc loop
     1510          DO  j = nys_z, nyn_z
     1511             !$acc loop vector( 32 )
     1512             DO  i = nxl_z, nxr_z
     1513                tri(i,j,0,1) = tric(i,j,0)
     1514             ENDDO
     1515          ENDDO
     1516          !$acc end kernels
     1517
     1518          DO  k = 1, nz-1
     1519             !$acc kernels present( tri, tric )
     1520             !$acc loop
     1521             DO  j = nys_z, nyn_z
     1522                !$acc loop vector( 32 )
     1523                DO  i = nxl_z, nxr_z
     1524                   tri(i,j,k,2) = ddzuw(k,1) / tri(i,j,k-1,1)
     1525                   tri(i,j,k,1) = tric(i,j,k) - ddzuw(k-1,2) * tri(i,j,k,2)
     1526                ENDDO
     1527             ENDDO
     1528             !$acc end kernels
     1529          ENDDO
     1530
     1531    END SUBROUTINE split
     1532
     1533#endif
     1534
    15341535 END MODULE poisfft_mod
  • palm/trunk/SOURCE/poisfft_hybrid.f90

    r1107 r1111  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! poisfft_hybrid_ini is now called internally from poisfft_hybrid,
     23! ibc_p_b = 2 removed
    2324!
    2425! Former revisions:
     
    115116                     tasks_per_logical_node = -1    ! default no cluster
    116117
     118    LOGICAL, SAVE ::  poisfft_initialized = .FALSE.
     119
    117120
    118121    PRIVATE
     
    283286       ENDIF
    284287
     288       poisfft_initialized = .TRUE.
     289
    285290    END SUBROUTINE poisfft_hybrid_ini
    286291
     
    294299
    295300       REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar
     301
     302       IF ( .NOT. poisfft_initialized )  CALL poisfft_hybrid_ini
    296303
    297304       IF ( host(1:3) == 'nec' )  THEN
     
    961968             ENDDO
    962969          ENDDO
    963           IF ( ibc_p_b == 1  .OR.  ibc_p_b == 2 )  THEN
     970          IF ( ibc_p_b == 1 )  THEN
    964971             DO  i = 0,nx
    965972                tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
  • palm/trunk/SOURCE/pres.f90

    r1093 r1111  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! openACC statements added,
     23! ibc_p_b = 2 removed
    2324!
    2425! Former revisions:
     
    304305    ELSE
    305306       !$OMP PARALLEL DO SCHEDULE( STATIC )
     307       !$acc kernels present( d )
     308       !$acc loop
    306309       DO  i = nxl, nxr
    307310          DO  j = nys, nyn
     311             !$acc loop vector(32)
    308312             DO  k = nzb+1, nzt
    309313                d(k,j,i) = 0.0
     
    311315          ENDDO
    312316       ENDDO
     317       !$acc end kernels
    313318    ENDIF
    314319
     
    328333          ENDDO
    329334!
    330 !--       Additional pressure boundary condition at the bottom boundary for
    331 !--       inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
    332 !--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
    333 !--       This condition must not be applied at the start of a run, because then
    334 !--       flow_statistics has not yet been called and thus sums = 0.
    335           IF ( ibc_p_b == 2  .AND.  sums(nzb+1,4) /= 0.0 )  THEN
    336              k = nzb_s_inner(j,i)
    337              d(k+1,j,i) = d(k+1,j,i) + (                                     &
    338                                          ( usws(j,i+1) - usws(j,i) ) * ddx   &
    339                                        + ( vsws(j+1,i) - vsws(j,i) ) * ddy   &
    340                                        - g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
    341                                          sums(k+1,4)                         &
    342                                        ) * ddzw(k+1) * ddt_3d * d_weight_pres
    343           ENDIF
    344 
    345 !
    346335!--       Compute possible PE-sum of divergences for flow_statistics
    347336          DO  k = nzb_s_inner(j,i)+1, nzt
     
    357346    !$OMP END PARALLEL
    358347#else
    359     IF ( ibc_p_b == 2 .AND. sums(nzb+1,4) /= 0.0 )  THEN
    360        !$OMP PARALLEL PRIVATE (i,j,k)
    361        !$OMP DO SCHEDULE( STATIC )
    362        DO  i = nxl, nxr
    363           DO  j = nys, nyn
    364              DO  k = nzb_s_inner(j,i)+1, nzt
    365              d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
    366                           ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
    367                           ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
    368                                                                 * d_weight_pres
    369              ENDDO
    370           ENDDO
    371 !
    372 !--       Additional pressure boundary condition at the bottom boundary for
    373 !--       inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
    374 !--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
    375 !--       This condition must not be applied at the start of a run, because then
    376 !--       flow_statistics has not yet been called and thus sums = 0.
    377           DO  j = nys, nyn
    378               k = nzb_s_inner(j,i)
    379               d(k+1,j,i) = d(k+1,j,i) + (                        &
    380                              ( usws(j,i+1) - usws(j,i) ) * ddx   &
    381                            + ( vsws(j+1,i) - vsws(j,i) ) * ddy   &
    382                            - g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
    383                              sums(k+1,4)                         &
    384                                         ) * ddzw(k+1) * ddt_3d   &
    385                                           * d_weight_pres
    386           ENDDO
    387        ENDDO
    388        !$OMP END PARALLEL
    389 
    390     ELSE
    391 
    392        !$OMP PARALLEL PRIVATE (i,j,k)
    393        !$OMP DO SCHEDULE( STATIC )
    394        DO  i = nxl, nxr
    395           DO  j = nys, nyn
    396              DO  k = nzb_s_inner(j,i)+1, nzt
     348
     349    !$OMP PARALLEL PRIVATE (i,j,k)
     350    !$OMP DO SCHEDULE( STATIC )
     351    !$acc kernels present( d, ddzw, nzb_s_inner, u, v, w )
     352    !$acc loop
     353    DO  i = nxl, nxr
     354       DO  j = nys, nyn
     355          !$acc loop vector(32)
     356          DO  k = 1, nzt
     357             IF ( k > nzb_s_inner(j,i) )  THEN
    397358                d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
    398                           ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
    399                           ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
    400                                                                 * d_weight_pres
    401              ENDDO
    402           ENDDO
    403        ENDDO
    404        !$OMP END PARALLEL
    405 
    406     ENDIF
     359                           ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
     360                           ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
     361                           * d_weight_pres
     362             ENDIF
     363          ENDDO
     364       ENDDO
     365    ENDDO
     366    !$acc end kernels
     367    !$OMP END PARALLEL
    407368
    408369!
     
    439400!--       Solver for 2d-decomposition
    440401          CALL poisfft( d, tend )
     402          !$acc update host( d )
    441403       ELSEIF ( psolver == 'poisfft_hybrid' )  THEN
    442404!
     
    466428!
    467429!--       Neumann (dp/dz = 0)
    468           !$OMP PARALLEL DO
    469           DO  i = nxlg, nxrg
    470              DO  j = nysg, nyng
    471                 tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
    472              ENDDO
    473           ENDDO
    474 
    475        ELSEIF ( ibc_p_b == 2 )  THEN
    476 !
    477 !--       Neumann condition for inhomogeneous surfaces,
    478 !--       here currently still in the form of a zero gradient. Actually
    479 !--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0 would have to be used for
    480 !--       the computation (cf. above: computation of divergences).
    481430          !$OMP PARALLEL DO
    482431          DO  i = nxlg, nxrg
     
    611560!-- Correction of the provisional velocities with the current perturbation
    612561!-- pressure just computed
     562    !$acc update host( u, v, w )
    613563    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc .OR. bc_ns_cyc ) )  THEN
    614564       volume_flow_l(1) = 0.0
     
    748698    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
    749699
     700    !$acc update device( u, v, w )
     701
    750702    CALL cpu_log( log_point(8), 'pres', 'stop' )
    751703   
  • palm/trunk/SOURCE/prognostic_equations.f90

    r1107 r1111  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! update directives for prognostic quantities removed
    2323!
    2424! Former revisions:
     
    15741574
    15751575    CALL cpu_log( log_point(5), 'u-equation', 'stop' )
    1576     !$acc update host( u_p )
    15771576
    15781577!
     
    16371636
    16381637    CALL cpu_log( log_point(6), 'v-equation', 'stop' )
    1639     !$acc update host( v_p )
    16401638
    16411639!
     
    17001698
    17011699    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
    1702     !$acc update host( w_p )
    17031700
    17041701
     
    17961793
    17971794       CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
    1798        !$acc update host( pt_p )
    17991795
    18001796    ENDIF
     
    20462042
    20472043       CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
    2048        !$acc update host( e_p )
    20492044
    20502045    ENDIF
  • palm/trunk/SOURCE/swap_timelevel.f90

    r1054 r1111  
    2020! Current revisions:
    2121! -----------------
     22! openACC directives added
    2223!
    2324! Former revisions:
     
    8283    CALL cpu_log( log_point(28), 'swap_timelevel (nop)', 'start' )
    8384
     85    !$acc kernels present( pt, pt_p, u, u_p, v, v_p, w, w_p )
    8486    u  = u_p
    8587    v  = v_p
    8688    w  = w_p
    8789    pt = pt_p
     90    !$acc end kernels
    8891    IF ( .NOT. constant_diffusion )  THEN
     92       !$acc kernels present( e, e_p )
    8993       e = e_p
     94       !$acc end kernels
    9095    ENDIF
    9196    IF ( ocean )  THEN
  • palm/trunk/SOURCE/time_integration.f90

    r1093 r1111  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! +internal timestep counter for cpu statistics added,
     23! openACC directives updated
    2324!
    2425! Former revisions:
     
    238239!--       Exchange of ghost points (lateral boundary conditions)
    239240          CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'start' )
     241          !$acc update host( e_p, pt_p, u_p, v_p, w_p )
    240242          CALL exchange_horiz( u_p, nbgp )
    241243          CALL exchange_horiz( v_p, nbgp )
     
    272274!
    273275!--       Swap the time levels in preparation for the next time step.
     276          !$acc update device( e_p, pt_p, u_p, v_p, w_p )
    274277          CALL swap_timelevel
    275278
     
    298301             time_disturb = time_disturb + dt_3d
    299302             IF ( time_disturb >= dt_disturb )  THEN
     303                !$acc update host( u, v )
    300304                IF ( hom(nzb+5,1,pr_palm,0) < disturbance_energy_limit )  THEN
    301305                   CALL disturb_field( nzb_u_inner, tend, u )
     
    310314                   dist_range = 0
    311315                ENDIF
     316                !$acc update device( u, v )
    312317                time_disturb = time_disturb - dt_disturb
    313318             ENDIF
     
    321326             CALL pres
    322327          ENDIF
    323 !
    324 !--       Update device memory for calculating diffusion quantities and for next
    325 !--       timestep
    326           !$acc update device( e, pt, u, v, w )
    327           !$acc update device( q )  if ( allocated( q ) )
    328328
    329329!
     
    351351                CALL prandtl_fluxes
    352352                CALL cpu_log( log_point(19), 'prandtl_fluxes', 'stop' )
    353 !
    354 !++             Statistics still require updates on host
    355                 !$acc update host( qs, qsws, rif, shf, ts )
    356353             ENDIF
    357354
     
    369366             ENDIF
    370367             CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
    371 !
    372 !++          Statistics still require update of diffusivities on host
    373              !$acc update host( kh, km )
    374368
    375369          ENDIF
     
    379373!
    380374!--    Increase simulation time and output times
     375       nr_timesteps_this_run      = nr_timesteps_this_run + 1
    381376       current_timestep_number    = current_timestep_number + 1
    382377       simulated_time             = simulated_time   + dt_3d
  • palm/trunk/SOURCE/transpose.f90

    r1107 r1111  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! openACC directives added,
     23! resorting data from/to work changed, work got 4 dimensions instead of 1
    2324!
    2425! Former revisions:
     
    8182    IMPLICIT NONE
    8283
    83     INTEGER ::  i, j, k, l, m, ys
     84    INTEGER ::  i, j, k, l, ys
    8485   
    85     REAL ::  f_in(0:nx,nys_x:nyn_x,nzb_x:nzt_x),   &
    86              f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx),  &
    87              f_out(0:ny,nxl_y:nxr_y,nzb_y:nzt_y),  &
    88              work(nnx*nny*nnz)
     86    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)
     87
     88    REAL, DIMENSION(nyn_x-nys_x+1,nzb_y:nzt_y,nxl_y:nxr_y,0:pdims(2)-1) ::  work
     89
     90    !$acc declare create( f_inv )
     91    REAL ::  f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx)
     92
    8993
    9094!
     
    9397!$OMP  PARALLEL PRIVATE ( i, j, k )
    9498!$OMP  DO
     99    !$acc kernels present( f_in )
     100    !$acc loop
    95101    DO  i = 0, nx
    96102       DO  k = nzb_x, nzt_x
     103          !$acc loop vector( 32 )
    97104          DO  j = nys_x, nyn_x
    98105             f_inv(j,k,i) = f_in(i,j,k)
     
    100107       ENDDO
    101108    ENDDO
     109    !$acc end kernels
    102110!$OMP  END PARALLEL
    103111
     
    109117       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
    110118       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    111        CALL MPI_ALLTOALL( f_inv(nys_x,nzb_x,0), sendrecvcount_xy, MPI_REAL, &
    112                           work(1),              sendrecvcount_xy, MPI_REAL, &
     119       !$acc update host( f_inv )
     120       CALL MPI_ALLTOALL( f_inv(nys_x,nzb_x,0),  sendrecvcount_xy, MPI_REAL, &
     121                          work(1,nzb_y,nxl_y,0), sendrecvcount_xy, MPI_REAL, &
    113122                          comm1dy, ierr )
     123       !$acc update device( work )
    114124       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
    115125
    116126!
    117127!--    Reorder transposed array
    118 !$OMP  PARALLEL PRIVATE ( i, j, k, l, m, ys )
     128!$OMP  PARALLEL PRIVATE ( i, j, k, l, ys )
    119129!$OMP  DO
    120130       DO  l = 0, pdims(2) - 1
    121           m  = l * ( nxr_y - nxl_y + 1 ) * ( nzt_y - nzb_y + 1 ) * &
    122                    ( nyn_x - nys_x + 1 )
    123131          ys = 0 + l * ( nyn_x - nys_x + 1 )
     132          !$acc kernels present( f_out, work )
     133          !$acc loop
    124134          DO  i = nxl_y, nxr_y
    125135             DO  k = nzb_y, nzt_y
     136                !$acc loop vector( 32 )
    126137                DO  j = ys, ys + nyn_x - nys_x
    127                    m = m + 1
    128                    f_out(j,i,k) = work(m)
     138                   f_out(j,i,k) = work(j-ys+1,k,i,l)
    129139                ENDDO
    130140             ENDDO
    131141          ENDDO
     142          !$acc end kernels
    132143       ENDDO
    133144!$OMP  END PARALLEL
     
    140151!$OMP  PARALLEL PRIVATE ( i, j, k )
    141152!$OMP  DO
     153       !$acc kernels present( f_out )
     154       !$acc loop
    142155       DO  k = nzb_y, nzt_y
    143156          DO  i = nxl_y, nxr_y
     157             !$acc loop vector( 32 )
    144158             DO  j = 0, ny
    145159                f_out(j,i,k) = f_inv(j,k,i)
     
    147161          ENDDO
    148162       ENDDO
     163       !$acc end kernels
    149164!$OMP  END PARALLEL
    150165
     
    172187    IMPLICIT NONE
    173188
    174     INTEGER ::  i, j, k, l, m, xs
     189    INTEGER ::  i, j, k, l, xs
    175190   
    176     REAL ::  f_in(0:nx,nys_x:nyn_x,nzb_x:nzt_x),  &
    177              f_inv(nys:nyn,nxl:nxr,1:nz),         &
    178              f_out(1:nz,nys:nyn,nxl:nxr),         &
    179              work(nnx*nny*nnz)
     191    REAL ::  f_in(0:nx,nys_x:nyn_x,nzb_x:nzt_x), f_out(1:nz,nys:nyn,nxl:nxr)
     192
     193    REAL, DIMENSION(nys_x:nyn_x,nnx,nzb_x:nzt_x,0:pdims(1)-1) ::  work
     194
     195    !$acc declare create( f_inv )
     196    REAL ::  f_inv(nys:nyn,nxl:nxr,1:nz)
    180197
    181198
     
    188205!
    189206!--    Reorder input array for transposition
    190 !$OMP  PARALLEL PRIVATE ( i, j, k, l, m, xs )
     207!$OMP  PARALLEL PRIVATE ( i, j, k, l, xs )
    191208!$OMP  DO
    192209       DO  l = 0, pdims(1) - 1
    193           m  = l * ( nzt_x - nzb_x + 1 ) * nnx * ( nyn_x - nys_x + 1 )
    194210          xs = 0 + l * nnx
     211          !$acc kernels present( f_in, work )
     212          !$acc loop
    195213          DO  k = nzb_x, nzt_x
    196214             DO  i = xs, xs + nnx - 1
     215                !$acc loop vector( 32 )
    197216                DO  j = nys_x, nyn_x
    198                    m = m + 1
    199                    work(m) = f_in(i,j,k)
     217                   work(j,i-xs+1,k,l) = f_in(i,j,k)
    200218                ENDDO
    201219             ENDDO
    202220          ENDDO
     221          !$acc end kernels
    203222       ENDDO
    204223!$OMP  END PARALLEL
     
    208227       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
    209228       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    210        CALL MPI_ALLTOALL( work(1),          sendrecvcount_zx, MPI_REAL, &
    211                           f_inv(nys,nxl,1), sendrecvcount_zx, MPI_REAL, &
     229       !$acc update host( work )
     230       CALL MPI_ALLTOALL( work(nys_x,1,nzb_x,0), sendrecvcount_zx, MPI_REAL, &
     231                          f_inv(nys,nxl,1),      sendrecvcount_zx, MPI_REAL, &
    212232                          comm1dx, ierr )
     233       !$acc update device( f_inv )
    213234       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
    214235
     
    217238!$OMP  PARALLEL PRIVATE ( i, j, k )
    218239!$OMP  DO
     240       !$acc kernels present( f_out )
     241       !$acc loop
    219242       DO  k = 1, nz
    220243          DO  i = nxl, nxr
     244             !$acc loop vector( 32 )
    221245             DO  j = nys, nyn
    222246                f_out(k,j,i) = f_inv(j,i,k)
     
    224248          ENDDO
    225249       ENDDO
     250       !$acc end kernels
    226251!$OMP  END PARALLEL
    227252#endif
     
    233258!$OMP  PARALLEL PRIVATE ( i, j, k )
    234259!$OMP  DO
     260       !$acc kernels present( f_in )
     261       !$acc loop
    235262       DO  i = nxl, nxr
    236263          DO  j = nys, nyn
     264             !$acc loop vector( 32 )
    237265             DO  k = 1, nz
    238266                f_inv(j,i,k) = f_in(i,j,k)
     
    240268          ENDDO
    241269       ENDDO
    242 !$OMP  END PARALLEL
    243 
    244 !$OMP  PARALLEL PRIVATE ( i, j, k )
    245 !$OMP  DO
     270       !$acc end kernels
     271!$OMP  END PARALLEL
     272
     273!$OMP  PARALLEL PRIVATE ( i, j, k )
     274!$OMP  DO
     275       !$acc kernels present( f_out )
     276       !$acc loop
    246277       DO  k = 1, nz
    247278          DO  i = nxl, nxr
     279             !$acc loop vector( 32 )
    248280             DO  j = nys, nyn
    249281                f_out(k,j,i) = f_inv(j,i,k)
     
    251283          ENDDO
    252284       ENDDO
     285       !$acc end kernels
    253286!$OMP  END PARALLEL
    254287
     
    276309    IMPLICIT NONE
    277310
    278     INTEGER ::  i, j, k, l, m, ys
     311    INTEGER ::  i, j, k, l, ys
    279312   
    280     REAL ::  f_in(0:ny,nxl_y:nxr_y,nzb_y:nzt_y),  &
    281              f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx), &
    282              f_out(0:nx,nys_x:nyn_x,nzb_x:nzt_x), &
    283              work(nnx*nny*nnz)
     313    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)
     314
     315    REAL, DIMENSION(nyn_x-nys_x+1,nzb_y:nzt_y,nxl_y:nxr_y,0:pdims(2)-1) ::  work
     316
     317    !$acc declare create( f_inv )
     318    REAL ::  f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx)
     319
    284320
    285321    IF ( numprocs /= 1 )  THEN
     
    288324!
    289325!--    Reorder input array for transposition
    290 !$OMP  PARALLEL PRIVATE ( i, j, k, l, m, ys )
     326!$OMP  PARALLEL PRIVATE ( i, j, k, l, ys )
    291327!$OMP  DO
    292328       DO  l = 0, pdims(2) - 1
    293           m  = l * ( nxr_y - nxl_y + 1 ) * ( nzt_y - nzb_y + 1 ) * &
    294                    ( nyn_x - nys_x + 1 )
    295329          ys = 0 + l * ( nyn_x - nys_x + 1 )
     330          !$acc kernels present( f_in, work )
     331          !$acc loop
    296332          DO  i = nxl_y, nxr_y
    297333             DO  k = nzb_y, nzt_y
     334                !$acc loop vector( 32 )
    298335                DO  j = ys, ys + nyn_x - nys_x
    299                    m = m + 1
    300                    work(m) = f_in(j,i,k)
     336                   work(j-ys+1,k,i,l) = f_in(j,i,k)
    301337                ENDDO
    302338             ENDDO
    303339          ENDDO
     340          !$acc end kernels
    304341       ENDDO
    305342!$OMP  END PARALLEL
     
    309346       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
    310347       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    311        CALL MPI_ALLTOALL( work(1),              sendrecvcount_xy, MPI_REAL, &
    312                           f_inv(nys_x,nzb_x,0), sendrecvcount_xy, MPI_REAL, &
     348       !$acc update host( work )
     349       CALL MPI_ALLTOALL( work(1,nzb_y,nxl_y,0), sendrecvcount_xy, MPI_REAL, &
     350                          f_inv(nys_x,nzb_x,0),  sendrecvcount_xy, MPI_REAL, &
    313351                          comm1dy, ierr )
     352       !$acc update device( f_inv )
    314353       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
    315354#endif
     
    321360!$OMP  PARALLEL PRIVATE ( i, j, k )
    322361!$OMP  DO
     362       !$acc kernels present( f_in )
     363       !$acc loop
    323364       DO  i = nxl_y, nxr_y
    324365          DO  k = nzb_y, nzt_y
     366             !$acc loop vector( 32 )
    325367             DO  j = 0, ny
    326368                f_inv(j,k,i) = f_in(j,i,k)
     
    328370          ENDDO
    329371       ENDDO
     372       !$acc end kernels
    330373!$OMP  END PARALLEL
    331374
     
    336379!$OMP  PARALLEL PRIVATE ( i, j, k )
    337380!$OMP  DO
     381    !$acc kernels present( f_out )
     382    !$acc loop
    338383    DO  i = 0, nx
    339384       DO  k = nzb_x, nzt_x
     385          !$acc loop vector( 32 )
    340386          DO  j = nys_x, nyn_x
    341387             f_out(i,j,k) = f_inv(j,k,i)
     
    343389       ENDDO
    344390    ENDDO
     391    !$acc end kernels
    345392!$OMP  END PARALLEL
    346393
     
    434481    IMPLICIT NONE
    435482
    436     INTEGER ::  i, j, k, l, m, zs
     483    INTEGER ::  i, j, k, l, zs
    437484   
    438     REAL ::  f_in(0:ny,nxl_y:nxr_y,nzb_y:nzt_y),  &
    439              f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny), &
    440              f_out(nxl_z:nxr_z,nys_z:nyn_z,1:nz), &
    441              work(nnx*nny*nnz)
     485    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)
     486
     487    REAL, DIMENSION(nxl_z:nxr_z,nzt_y-nzb_y+1,nys_z:nyn_z,0:pdims(1)-1) ::  work
     488
     489    !$acc declare create( f_inv )
     490    REAL ::  f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny)
     491
    442492
    443493!
     
    446496!$OMP  PARALLEL PRIVATE ( i, j, k )
    447497!$OMP  DO
     498    !$acc kernels present( f_in )
     499    !$acc loop
    448500    DO  j = 0, ny
    449501       DO  k = nzb_y, nzt_y
     502          !$acc loop vector( 32 )
    450503          DO  i = nxl_y, nxr_y
    451504             f_inv(i,k,j) = f_in(j,i,k)
     
    453506       ENDDO
    454507    ENDDO
     508    !$acc end kernels
    455509!$OMP  END PARALLEL
    456510
     
    464518!$OMP  PARALLEL PRIVATE ( i, j, k )
    465519!$OMP  DO
     520       !$acc kernels present( f_out )
     521       !$acc loop
    466522       DO  j = 0, ny
    467523          DO  k = nzb_y, nzt_y
     524             !$acc loop vector( 32 )
    468525             DO  i = nxl_y, nxr_y
    469526                f_out(i,j,k) = f_inv(i,k,j)
     
    471528          ENDDO
    472529       ENDDO
     530       !$acc end kernels
    473531!$OMP  END PARALLEL
    474532
     
    480538       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
    481539       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    482        CALL MPI_ALLTOALL( f_inv(nxl_y,nzb_y,0), sendrecvcount_yz, MPI_REAL, &
    483                           work(1),              sendrecvcount_yz, MPI_REAL, &
     540       !$acc update host( f_inv )
     541       CALL MPI_ALLTOALL( f_inv(nxl_y,nzb_y,0),  sendrecvcount_yz, MPI_REAL, &
     542                          work(nxl_z,1,nys_z,0), sendrecvcount_yz, MPI_REAL, &
    484543                          comm1dx, ierr )
     544       !$acc update device( work )
    485545       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
    486546
    487547!
    488548!--    Reorder transposed array
    489 !$OMP  PARALLEL PRIVATE ( i, j, k, l, m, zs )
     549!$OMP  PARALLEL PRIVATE ( i, j, k, l, zs )
    490550!$OMP  DO
    491551       DO  l = 0, pdims(1) - 1
    492           m  = l * ( nyn_z - nys_z + 1 ) * ( nzt_y - nzb_y + 1 ) * &
    493                    ( nxr_z - nxl_z + 1 )
    494552          zs = 1 + l * ( nzt_y - nzb_y + 1 )
     553          !$acc kernels present( f_out, work )
     554          !$acc loop
    495555          DO  j = nys_z, nyn_z
    496556             DO  k = zs, zs + nzt_y - nzb_y
     557                !$acc loop vector( 32 )
    497558                DO  i = nxl_z, nxr_z
    498                    m = m + 1
    499                    f_out(i,j,k) = work(m)
     559                   f_out(i,j,k) = work(i,k-zs+1,j,l)
    500560                ENDDO
    501561             ENDDO
    502562          ENDDO
     563          !$acc end kernels
    503564       ENDDO
    504565!$OMP  END PARALLEL
     
    528589    IMPLICIT NONE
    529590
    530     INTEGER ::  i, j, k, l, m, xs
     591    INTEGER ::  i, j, k, l, xs
    531592   
    532     REAL ::  f_in(1:nz,nys:nyn,nxl:nxr), f_out(0:nx,nys_x:nyn_x,nzb_x:nzt_x), &
    533              work(nnx*nny*nnz)
    534 
    535     !$acc declare create ( f_inv )
     593    REAL ::  f_in(1:nz,nys:nyn,nxl:nxr), f_out(0:nx,nys_x:nyn_x,nzb_x:nzt_x)
     594
     595    REAL, DIMENSION(nys_x:nyn_x,nnx,nzb_x:nzt_x,0:pdims(1)-1) ::  work
     596
     597    !$acc declare create( f_inv )
    536598    REAL ::  f_inv(nys:nyn,nxl:nxr,1:nz)
    537599
     
    552614       ENDDO
    553615    ENDDO
     616    !$acc end kernels
    554617!$OMP  END PARALLEL
    555618
     
    573636          ENDDO
    574637       ENDDO
     638       !$acc end kernels
    575639!$OMP  END PARALLEL
    576640
     
    582646       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
    583647       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    584        CALL MPI_ALLTOALL( f_inv(nys,nxl,1), sendrecvcount_zx, MPI_REAL, &
    585                           work(1),          sendrecvcount_zx, MPI_REAL, &
     648       !$acc update host( f_inv )
     649       CALL MPI_ALLTOALL( f_inv(nys,nxl,1),      sendrecvcount_zx, MPI_REAL, &
     650                          work(nys_x,1,nzb_x,0), sendrecvcount_zx, MPI_REAL, &
    586651                          comm1dx, ierr )
     652       !$acc update device( work )
    587653       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
    588654
    589655!
    590656!--    Reorder transposed array
    591 !$OMP  PARALLEL PRIVATE ( i, j, k, l, m, xs )
     657!$OMP  PARALLEL PRIVATE ( i, j, k, l, xs )
    592658!$OMP  DO
    593659       DO  l = 0, pdims(1) - 1
    594           m  = l * ( nzt_x - nzb_x + 1 ) * nnx * ( nyn_x - nys_x + 1 )
    595660          xs = 0 + l * nnx
     661          !$acc kernels present( f_out, work )
     662          !$acc loop
    596663          DO  k = nzb_x, nzt_x
    597664             DO  i = xs, xs + nnx - 1
     665                !$acc loop vector( 32 )
    598666                DO  j = nys_x, nyn_x
    599                    m = m + 1
    600                    f_out(i,j,k) = work(m)
     667                   f_out(i,j,k) = work(j,i-xs+1,k,l)
    601668                ENDDO
    602669             ENDDO
    603670          ENDDO
     671          !$acc end kernels
    604672       ENDDO
    605673!$OMP  END PARALLEL
     
    629697    IMPLICIT NONE
    630698
    631     INTEGER ::  i, j, k, l, m, zs
     699    INTEGER ::  i, j, k, l, zs
    632700   
    633     REAL ::  f_in(nxl_z:nxr_z,nys_z:nyn_z,1:nz),  &
    634              f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny), &
    635              f_out(0:ny,nxl_y:nxr_y,nzb_y:nzt_y), &
    636              work(nnx*nny*nnz)
     701    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)
     702
     703    REAL, DIMENSION(nxl_z:nxr_z,nzt_y-nzb_y+1,nys_z:nyn_z,0:pdims(1)-1) ::  work
     704
     705    !$acc declare create( f_inv )
     706    REAL ::  f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny)
     707
    637708
    638709!
     
    644715!
    645716!--    Reorder input array for transposition
    646 !$OMP  PARALLEL PRIVATE ( i, j, k, l, m, zs )
     717!$OMP  PARALLEL PRIVATE ( i, j, k, l, zs )
    647718!$OMP  DO
    648719       DO  l = 0, pdims(1) - 1
    649           m  = l * ( nyn_z - nys_z + 1 ) * ( nzt_y - nzb_y + 1 ) * &
    650                    ( nxr_z - nxl_z + 1 )
    651720          zs = 1 + l * ( nzt_y - nzb_y + 1 )
     721          !$acc kernels present( f_in, work )
     722          !$acc loop
    652723          DO  j = nys_z, nyn_z
    653724             DO  k = zs, zs + nzt_y - nzb_y
     725                !$acc loop vector( 32 )
    654726                DO  i = nxl_z, nxr_z
    655                    m = m + 1
    656                    work(m) = f_in(i,j,k)
     727                   work(i,k-zs+1,j,l) = f_in(i,j,k)
    657728                ENDDO
    658729             ENDDO
    659730          ENDDO
     731          !$acc end kernels
    660732       ENDDO
    661733!$OMP  END PARALLEL
     
    665737       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
    666738       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    667        CALL MPI_ALLTOALL( work(1),              sendrecvcount_yz, MPI_REAL, &
    668                           f_inv(nxl_y,nzb_y,0), sendrecvcount_yz, MPI_REAL, &
     739       !$acc update host( work )
     740       CALL MPI_ALLTOALL( work(nxl_z,1,nys_z,0), sendrecvcount_yz, MPI_REAL, &
     741                          f_inv(nxl_y,nzb_y,0),  sendrecvcount_yz, MPI_REAL, &
    669742                          comm1dx, ierr )
     743       !$acc update device( f_inv )
    670744       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
    671745#endif
     
    676750!$OMP  PARALLEL PRIVATE ( i, j, k )
    677751!$OMP  DO
     752       !$acc kernels present( f_in )
     753       !$acc loop
    678754       DO  k = nzb_y, nzt_y
    679755          DO  j = 0, ny
     756             !$acc loop vector( 32 )
    680757             DO  i = nxl_y, nxr_y
    681758                f_inv(i,k,j) = f_in(i,j,k)
     
    683760          ENDDO
    684761       ENDDO
     762       !$acc end kernels
    685763!$OMP  END PARALLEL
    686764
     
    691769!$OMP  PARALLEL PRIVATE ( i, j, k )
    692770!$OMP  DO
     771    !$acc kernels present( f_out )
     772    !$acc loop
    693773    DO  k = nzb_y, nzt_y
    694774       DO  i = nxl_y, nxr_y
     775          !$acc loop vector( 32 )
    695776          DO  j = 0, ny
    696777             f_out(j,i,k) = f_inv(i,k,j)
     
    698779       ENDDO
    699780    ENDDO
     781    !$acc end kernels
    700782!$OMP  END PARALLEL
    701783
Note: See TracChangeset for help on using the changeset viewer.