Changeset 1216 for palm/trunk/SOURCE/tridia_solver.f90

Timestamp:

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

• palm/trunk/SOURCE/tridia_solver.f90 (modified) (3 diffs)

palm/trunk/SOURCE/tridia_solver.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-!
+! +tridia_substi_overlap for handling overlapping fft / transposition
 !
 ! Former revisions:
@@ -57 +57 @@
     END INTERFACE tridia_substi
 
-    PUBLIC  tridia_substi, tridia_init, tridia_1dd
+    INTERFACE tridia_substi_overlap
+       MODULE PROCEDURE tridia_substi_overlap
+    END INTERFACE tridia_substi_overlap
+
+    PUBLIC  tridia_substi, tridia_substi_overlap, tridia_init, tridia_1dd
 
  CONTAINS
@@ -260 +264 @@
 
 
+    SUBROUTINE tridia_substi_overlap( ar, jj )
+
+!------------------------------------------------------------------------------!
+! Substitution (Forward and Backward) (Thomas algorithm)
+!------------------------------------------------------------------------------!
+
+          USE arrays_3d,  ONLY: tri
+          USE control_parameters
+
+          IMPLICIT NONE
+
+          INTEGER ::  i, j, jj, k
+
+          REAL    ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
+
+          !$acc declare create( ar1 )
+          REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1)   ::  ar1
+
+!
+!--       Forward substitution
+          DO  k = 0, nz - 1
+             !$acc kernels present( ar, tri )
+             !$acc loop
+             DO  j = nys_z, nyn_z
+                DO  i = nxl_z, nxr_z
+
+                   IF ( k == 0 )  THEN
+                      ar1(i,j,k) = ar(i,j,k+1)
+                   ELSE
+                      ar1(i,j,k) = ar(i,j,k+1) - tri(i,jj,k,2) * ar1(i,j,k-1)
+                   ENDIF
+
+                ENDDO
+             ENDDO
+             !$acc end kernels
+          ENDDO
+
+!
+!--       Backward substitution
+!--       Note, the 1.0E-20 in the denominator is due to avoid divisions
+!--       by zero appearing if the pressure bc is set to neumann at the top of
+!--       the model domain.
+          DO  k = nz-1, 0, -1
+             !$acc kernels present( ar, tri )
+             !$acc loop
+             DO  j = nys_z, nyn_z
+                DO  i = nxl_z, nxr_z
+
+                   IF ( k == nz-1 )  THEN
+                      ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20 )
+                   ELSE
+                      ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
+                              / tri(i,jj,k,1)
+                   ENDIF
+                ENDDO
+             ENDDO
+             !$acc end kernels
+          ENDDO
+
+!
+!--       Indices i=0, j=0 correspond to horizontally averaged pressure.
+!--       The respective values of ar should be zero at all k-levels if
+!--       acceleration of horizontally averaged vertical velocity is zero.
+          IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
+             IF ( nys_z == 0  .AND.  nxl_z == 0 )  THEN
+                !$acc kernels loop present( ar )
+                DO  k = 1, nz
+                   ar(nxl_z,nys_z,k) = 0.0
+                ENDDO
+             ENDIF
+          ENDIF
+
+    END SUBROUTINE tridia_substi_overlap
+
+
     SUBROUTINE split