Ignore:
Timestamp:
Aug 26, 2013 9:31:42 AM (8 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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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.