Index: palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3298)
+++ palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3299)
@@ -25,4 +25,10 @@
 ! -----------------
 ! $Id$
+! - removed global array wall_flags_0_global, hence reduced accuracy of l_wall
+!   calculation
+! - removed maxloc call as this produced different results for different
+!   compiler options
+! 
+! 3294 2018-10-01 02:37:10Z raasch
 ! changes concerning modularization of ocean option
 ! 
@@ -1320,4 +1326,6 @@
 ! -----------------------------------------------------------------------------!
 !> Pre-computation of grid-dependent and near-wall mixing length.
+!> @todo consider walls in horizontal direction at a distance further than a
+!>       single grid point (RANS mode)
 !------------------------------------------------------------------------------!
  SUBROUTINE tcm_init_mixing_length
@@ -1362,5 +1370,4 @@
     INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k)
 
-    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: wall_flags_0_global !< wall_flags_0 of whole domain
     INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: wall_flags_dummy    !< dummy array required for MPI_ALLREDUCE command
 
@@ -1510,35 +1517,9 @@
 
 !
-!--    Gather topography information of whole domain
-       !> @todo reduce amount of data sent by MPI call
-       !>   By now, a whole global 3D-array is sent and received with
-       !>   MPI_ALLREDUCE although most of the array is 0. This can be
-       !>   drastically reduced if only the local subarray is sent and stored
-       !>   in a global array. For that, an MPI data type or subarray must be
-       !>   defined.
-       !>   2018-03-19, gronemeier
-       ALLOCATE( wall_flags_0_global(nzb:nzt+1,0:ny,0:nx) )
-
-#if defined ( __parallel )
-       ALLOCATE( wall_flags_dummy(nzb:nzt+1,0:ny,0:nx) )
-       wall_flags_dummy = 0
-       wall_flags_dummy(nzb:nzt+1,nys:nyn,nxl:nxr) =  &
-           wall_flags_0(nzb:nzt+1,nys:nyn,nxl:nxr)
-
-       CALL MPI_ALLREDUCE( wall_flags_dummy,                  &
-                           wall_flags_0_global,               &
-                           (nzt-nzb+2)*(ny+1)*(nx+1),         &
-                           MPI_INTEGER, MPI_SUM, comm2d, ierr )
-       DEALLOCATE( wall_flags_dummy )
-#else
-       wall_flags_0_global(nzb:nzt+1,nys:nyn,nxl:nxr) =  &
-              wall_flags_0(nzb:nzt+1,nys:nyn,nxl:nxr)
-#endif
-!
-!--    Get height level of highest topography
-       DO  i = 0, nx
-          DO  j = 0, ny
+!--    Get height level of highest topography within local subdomain
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              DO  k = nzb+1, nzt-1
-                IF ( .NOT. BTEST( wall_flags_0_global(k,j,i), 0 ) .AND.  &
+                IF ( .NOT. BTEST( wall_flags_0(k,j,i), 0 ) .AND.  &
                      k > k_max_topo )  &
                    k_max_topo = k
@@ -1614,9 +1595,17 @@
 !--                   First, limit size of data copied to vicinity by the domain
 !--                   border
-                      rad_i_l = MIN( rad_i, i )
-                      rad_i_r = MIN( rad_i, nx-i )
-
-                      rad_j_s = MIN( rad_j, j )
-                      rad_j_n = MIN( rad_j, ny-j )
+                      !> @note limit copied area to 1 grid point in hor. dir.
+                      !>   Ignore walls in horizontal direction which are
+                      !>   further away than a single grid point. This allows to
+                      !>   only search within local subdomain without the need
+                      !>   of global topography information.
+                      !>   The error made by this assumption are acceptable at
+                      !>   the moment.
+                      !>   2018-10-01, gronemeier
+                      rad_i_l = MIN( 1, rad_i, i )
+                      rad_i_r = MIN( 1, rad_i, nx-i )
+
+                      rad_j_s = MIN( 1, rad_j, j )
+                      rad_j_n = MIN( 1, rad_j, ny-j )
 
                       CALL copy_into_vicinity( k, j, i,           &
@@ -1624,74 +1613,13 @@
                                                -rad_j_s, rad_j_n, &
                                                -rad_i_l, rad_i_r  )
-!
-!--                   In case of cyclic boundaries, copy parts into vicinity
-!--                   where vicinity reaches over the domain borders.
-                      IF ( bc_lr_cyc )  THEN
-!
-!--                      Vicinity reaches over left domain boundary
-                         IF ( rad_i > rad_i_l )  THEN
-                            CALL copy_into_vicinity( k, j, nx+rad_i_l+1, &
-                                                     -rad_k_b, rad_k_t,  &
-                                                     -rad_j_s, rad_j_n,  &
-                                                     -rad_i, -rad_i_l-1  )
-!
-!--                         ...and over southern domain boundary
-                            IF ( bc_ns_cyc .AND. rad_j > rad_j_s )  &
-                               CALL copy_into_vicinity( k, ny+rad_j_s+1,    &
-                                                        nx+rad_i_l+1,       &
-                                                        -rad_k_b, rad_k_t,  &
-                                                        -rad_j, -rad_j_s-1, &
-                                                        -rad_i, -rad_i_l-1  )
-!
-!--                         ...and over northern domain boundary
-                            IF ( bc_ns_cyc .AND. rad_j > rad_j_n )  &
-                               CALL copy_into_vicinity( k, 0-rad_j_n-1,    &
-                                                        nx+rad_i_l+1,      &
-                                                        -rad_k_b, rad_k_t, &
-                                                         rad_j_n+1, rad_j, &
-                                                        -rad_i, -rad_i_l-1 )
-                         ENDIF
-!
-!--                      Vicinity reaches over right domain boundary
-                         IF ( rad_i > rad_i_r )  THEN
-                            CALL copy_into_vicinity( k, j, 0-rad_i_r-1, &
-                                                     -rad_k_b, rad_k_t, &
-                                                     -rad_j_s, rad_j_n, &
-                                                      rad_i_r+1, rad_i  )
-!
-!--                         ...and over southern domain boundary
-                            IF ( bc_ns_cyc .AND. rad_j > rad_j_s )  &
-                               CALL copy_into_vicinity( k, ny+rad_j_s+1,    &
-                                                        0-rad_i_r-1,        &
-                                                        -rad_k_b, rad_k_t,  &
-                                                        -rad_j, -rad_j_s-1, &
-                                                         rad_i_r+1, rad_i   )
-!
-!--                         ...and over northern domain boundary
-                            IF ( bc_ns_cyc .AND. rad_j > rad_j_n )  &
-                               CALL copy_into_vicinity( k, 0-rad_j_n-1,    &
-                                                        0-rad_i_r-1,       &
-                                                        -rad_k_b, rad_k_t, &
-                                                         rad_j_n+1, rad_j, &
-                                                         rad_i_r+1, rad_i  )
-                         ENDIF
-                      ENDIF
-
-                      IF ( bc_ns_cyc )  THEN
-!
-!--                      Vicinity reaches over southern domain boundary
-                         IF ( rad_j > rad_j_s )  &
-                            CALL copy_into_vicinity( k, ny+rad_j_s+1, i, &
-                                                     -rad_k_b, rad_k_t,  &
-                                                     -rad_j, -rad_j_s-1, &
-                                                     -rad_i_l, rad_i_r   )
-!
-!--                      Vicinity reaches over northern domain boundary
-                         IF ( rad_j > rad_j_n )  &
-                            CALL copy_into_vicinity( k, 0-rad_j_n-1, i, &
-                                                     -rad_k_b, rad_k_t, &
-                                                      rad_j_n+1, rad_j, &
-                                                      rad_i_l, rad_i_r  )
-                      ENDIF
+                      !> @note in case of cyclic boundaries, those parts of the
+                      !>   topography which lies beyond the domain borders but
+                      !>   still within the search radius still needs to be
+                      !>   copied into 'vicinity'. As the effective search
+                      !>   radius is limited to 1 at the moment, no further
+                      !>   copying is needed. Old implementation (prior to
+                      !>   2018-10-01) had this covered but used a global array.
+                      !>   2018-10-01, gronemeier
+
 !
 !--                   Search for walls only if there is any within vicinity
@@ -1791,6 +1719,4 @@
        ENDDO  !k loop
 
-       DEALLOCATE( wall_flags_0_global )
-
     ENDIF  !LES or RANS mode
 
@@ -1806,6 +1732,4 @@
 !> (pos_i/jj/kk), where (jj/kk) is the position of the maximum of 'array'
 !> closest to the origin (0/0) of 'array'.
-!> @todo this part of PALM does not reproduce the same results for optimized
-!>   and debug options for the compiler. This should be fixed
 !------------------------------------------------------------------------------!
     REAL(wp) FUNCTION shortest_distance( array, orientation, pos_i )
@@ -1817,6 +1741,10 @@
        INTEGER(iwp), INTENT(IN) :: pos_i     !< x position of the yz-plane 'array'
 
+       INTEGER(iwp) :: a                     !< loop index
+       INTEGER(iwp) :: b                     !< loop index
        INTEGER(iwp) :: jj                    !< loop index
 
+       INTEGER(KIND=1) :: maximum            !< maximum of array along z dimension
+
        INTEGER(iwp), DIMENSION(0:rad_j) :: loc_k !< location of closest wall along vertical dimension
 
@@ -1825,8 +1753,15 @@
 !
 !--    Get coordinate of first maximum along vertical dimension
-!--    at each y position of array.
-!--    Substract 1 because indices count from 1 instead of 0 by MAXLOC
-       loc_k = MAXLOC( array, DIM = 1) - 1
-
+!--    at each y position of array (similar to function maxloc but more stable).
+       DO  a = 0, rad_j
+          loc_k(a) = rad_k+1
+          maximum = MAXVAL( array(:,a) )
+          DO  b = 0, rad_k+1
+             IF ( array(b,a) == maxi )  THEN
+                loc_k(a) = b
+                EXIT
+             ENDIF
+          ENDDO
+       ENDDO
 !
 !--    Set distance to the default maximum value (=search radius)
@@ -1888,10 +1823,9 @@
        INTEGER(iwp) :: k   !< loop index
 
-
        DO  i = il, ir
           DO  j = js, jn
              DO  k = kb, kt
                 vicinity(k,j,i) = MERGE( 0, 1,               &
-                       BTEST( wall_flags_0_global(kp+k,jp+j,ip+i), 0 ) )
+                       BTEST( wall_flags_0(kp+k,jp+j,ip+i), 0 ) )
              ENDDO
           ENDDO
@@ -2536,5 +2470,5 @@
        IF ( intermediate_timestep_count == 2 )  diss_diff2(:,j,i) = tend(:,j,i) - diss_adve2(:,j,i) - diss_prod2(:,j,i)
        IF ( intermediate_timestep_count == 3 )  diss_diff3(:,j,i) = tend(:,j,i) - diss_adve3(:,j,i) - diss_prod3(:,j,i)
-       IF ( intermediate_timestep_count == 3 )  dummy3(:,j,i) = km(:,j,i)
+!        IF ( intermediate_timestep_count == 3 )  dummy3(:,j,i) = km(:,j,i)
 
        dum_dif = tend(:,j,i) - dum_adv - dum_pro                                !> @todo remove later
