Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 4307)
+++ palm/trunk/SOURCE/Makefile	(revision 4309)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# Add dependency to parallel random generator for synthetic turbulence generator
+# 
+# 4286 2019-10-30 16:01:14Z resler
 # delete boundary_conds, added missing dependencies
 # 
@@ -1077,4 +1080,5 @@
 	nesting_offl_mod.o \
 	pmc_interface_mod.o \
+	random_generator_parallel_mod.o \
 	surface_mod.o
 temperton_fft_mod.o: \
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 4307)
+++ palm/trunk/SOURCE/header.f90	(revision 4309)
@@ -20,9 +20,12 @@
 ! Current revisions:
 ! -----------------
-! replaced recycling_yshift by y_shift
+! 
 ! 
 ! Former revisions:
 ! -----------------
 ! $Id$
+! replaced recycling_yshift by y_shift
+! 
+! 4301 2019-11-22 12:09:09Z oliver.maas
 !
 ! 4297 2019-11-21 10:37:50Z oliver.maas
Index: palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4307)
+++ palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4309)
@@ -27,4 +27,7 @@
 ! -----------------
 ! $Id$
+! Typo
+! 
+! 4302 2019-11-22 13:15:56Z suehring
 ! Omit tall canopy mapped on top of buildings
 ! 
@@ -1115,5 +1118,5 @@
                          WRITE( message_string, * )                            &
                                           'Resolved plant-canopy is ' //       &
-                                          'is defined on top of a building ' //&
+                                          'defined on top of a building ' //   &
                                           '- lad is omitted at this grid '   //&
                                           'point: (i,j) = ', i, j
Index: palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
===================================================================
--- palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90	(revision 4307)
+++ palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90	(revision 4309)
@@ -25,4 +25,13 @@
 ! -----------------
 ! $Id$
+! Computation of velocity seeds optimized. This implies that random numbers 
+! are computed now using the parallel random number generator. Random numbers 
+! are now only computed and normalized locally, while distributed over all  
+! mpi ranks afterwards, instead of computing random numbers on a global array. 
+! Further, the number of calls for the time-consuming velocity-seed generation 
+! is reduced - now the left and right, as well as the north and south boundary 
+! share the same velocity-seed matrices.
+! 
+! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
 ! 
@@ -154,4 +163,5 @@
                neutral,                                                        &
                num_mean_inflow_profiles,                                       &
+               random_generator,                                               &
                rans_mode,                                                      &
                restart_string,                                                 &
@@ -206,4 +216,10 @@
     USE pmc_interface,                                                         &
         ONLY : rans_mode_parent
+        
+   USE random_generator_parallel,                                              &
+        ONLY:  random_dummy,                                                   &
+               random_number_parallel,                                         &
+               random_seed_parallel,                                           &
+               seq_random_array
 
     USE transpose_indices,                                                     &
@@ -227,8 +243,4 @@
     INTEGER(iwp) ::  id_stg_right       !< right lateral boundary core id in case of turbulence generator
     INTEGER(iwp) ::  id_stg_south       !< south lateral boundary core id in case of turbulence generator
-    INTEGER(iwp) ::  stg_type_xz        !< MPI type for full z range
-    INTEGER(iwp) ::  stg_type_xz_small  !< MPI type for small z range
-    INTEGER(iwp) ::  stg_type_yz        !< MPI type for full z range
-    INTEGER(iwp) ::  stg_type_yz_small  !< MPI type for small z range
     INTEGER(iwp) ::  merg               !< maximum length scale (in gp)
     INTEGER(iwp) ::  mergp              !< merg + nbgp
@@ -237,4 +249,8 @@
     INTEGER(iwp) ::  nzb_y_stg          !< lower bound of z coordinate (required for transposing z on PEs along y)
     INTEGER(iwp) ::  nzt_y_stg          !< upper bound of z coordinate (required for transposing z on PEs along y)
+    INTEGER(iwp) ::  stg_type_xz        !< MPI type for full z range
+    INTEGER(iwp) ::  stg_type_xz_small  !< MPI type for small z range
+    INTEGER(iwp) ::  stg_type_yz        !< MPI type for full z range
+    INTEGER(iwp) ::  stg_type_yz_small  !< MPI type for small z range
 
     INTEGER(iwp), DIMENSION(3) ::  nr_non_topo_xz = 0 !< number of non-topography grid points at xz cross-sections,
@@ -256,5 +272,4 @@
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nwy            !< length scale of w in y direction (in gp)
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nwz            !< length scale of w in z direction (in gp)
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  seed           !< seed of random number for rn-generator
 
     REAL(wp) ::  blend                    !< value to create gradually and smooth blending of Reynolds stress and length 
@@ -476,4 +491,11 @@
           CALL message( 'stg_check_parameters', 'PA0039', 1, 2, 0, 6, 0 )
        ENDIF
+!
+!--    Synthetic turbulence generator requires the parallel random generator
+       IF ( random_generator /= 'random-parallel' )  THEN
+          message_string = 'Using synthetic turbulence generator ' //          &
+                           'requires random_generator = random-parallel.'
+          CALL message( 'stg_check_parameters', 'PA0421', 1, 2, 0, 6, 0 )
+       ENDIF
 
     ENDIF
@@ -540,7 +562,5 @@
     INTEGER(iwp) :: newtype                  !< dummy MPI type
     INTEGER(iwp) :: realsize                 !< size of REAL variables
-    INTEGER(iwp) :: nseed                    !< dimension of random number seed
-    INTEGER(iwp) :: startseed = 1234567890   !< start seed for random number generator
-    
+
     INTEGER(iwp), DIMENSION(3) ::  nr_non_topo_xz_l = 0 !< number of non-topography grid points at xz-cross-section on subdomain
     INTEGER(iwp), DIMENSION(3) ::  nr_non_topo_yz_l = 0 !< number of non-topography grid points at yz-cross-section on subdomain
@@ -648,13 +668,4 @@
 
 #endif
-!
-!-- Define seed of random number
-    CALL RANDOM_SEED()
-    CALL RANDOM_SEED( size=nseed )
-    ALLOCATE( seed(1:nseed) )
-    DO  j = 1, nseed
-       seed(j) = startseed + j
-    ENDDO
-    CALL RANDOM_SEED(put=seed)
 !
 !-- Allocate required arrays.
@@ -1209,61 +1220,74 @@
 !-- Initial value of fu, fv, fw
     IF ( time_since_reference_point == 0.0_wp .AND. .NOT. velocity_seed_initialized )  THEN
-       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_left )
-       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_left )
-       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_left )
-
-       IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent      &
-                                .AND.  .NOT.  rans_mode ) )  THEN
-!
-!--       Generate turbulence at right boundary
-          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_right )
-          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_right )
-          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_right )
-!
-!--       Generate turbulence at north boundary
-          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_north )
-          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_north )
-          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_north )
-!
-!--       Generate turbulence at south boundary
-          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_south )
-          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_south )
-          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_south )
+!
+!--    Generate turbulence at the left and right boundary. Random numbers
+!--    for the yz-planes at the left/right boundary are generated by the
+!--    left-sided mpi ranks only. After random numbers are calculated, they 
+!--    are distributed to all other mpi ranks in the model, so that the
+!--    velocity seed matrices are available on all mpi ranks (i.e. also on the
+!--    right-sided boundary mpi ranks). In case of offline nesting, this implies, 
+!--    that the left- and the right-sided lateral boundary have the same initial
+!--    seeds. 
+!--    Note, in case of inflow from the right only, only turbulence at the left 
+!--    boundary is required.
+       IF ( .NOT. ( nesting_offline  .OR.                                      &
+                  ( child_domain .AND.  rans_mode_parent                       &
+                                 .AND.  .NOT.  rans_mode ) ) )  THEN
+
+          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_left )
+          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_left )
+          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_left )
+       ELSE
+          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz,               &
+                                     id_stg_left, id_stg_right )
+          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz,               &
+                                     id_stg_left, id_stg_right )
+          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz,               &
+                                     id_stg_left, id_stg_right )
+!
+!--       Generate turbulence at the south and north boundary. Random numbers
+!--       for the xz-planes at the south/north boundary are generated by the
+!--       south-sided mpi ranks only. Please see also comment above.
+          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz,               &
+                                     id_stg_south, id_stg_north )
+          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz,               &
+                                     id_stg_south, id_stg_north )
+          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz,               &
+                                     id_stg_south, id_stg_north )
        ENDIF
        velocity_seed_initialized = .TRUE.
     ENDIF
-    
-!
-!-- New set of fu, fv, fw
+!
+!-- New set of fu, fv, fw. Note, for inflow from the left side only, velocity
+!-- seeds are only required at the left boundary, while in case of offline
+!-- nesting or RANS-LES nesting velocity seeds are required also at the 
+!-- right, south and north boundaries. 
     IF ( stg_call )  THEN
-!      
-!--    Generate turbulence at left boundary (required in all applications, 
-!--    i.e. offline nesting and turbulent inflow from the left)
-       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left )
-       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
-       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
-       
-       IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent      &
-                             .AND.  .NOT.  rans_mode ) )  THEN
-!      
-!--       Generate turbulence at right boundary
-          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
-          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
-          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
-!      
-!--       Generate turbulence at north boundary
-          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
-          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
-          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
-!      
-!--       Generate turbulence at south boundary
-          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
-          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
-          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
+       IF ( .NOT. ( nesting_offline  .OR.                                      &
+                  ( child_domain .AND.  rans_mode_parent                       &
+                                 .AND.  .NOT.  rans_mode ) ) )  THEN
+          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left )
+          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
+          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
+
+       ELSE
+          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz,               &
+                                     id_stg_left, id_stg_right )
+          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz,               &
+                                     id_stg_left, id_stg_right )
+          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz,               &
+                                     id_stg_left, id_stg_right )
+
+          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz,               &
+                                     id_stg_south, id_stg_north )
+          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz,               &
+                                     id_stg_south, id_stg_north )
+          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz,               &
+                                     id_stg_south, id_stg_north )
        ENDIF
     ENDIF
     
 !
-!-- Turbulence generation at left and or right boundary
+!-- Turbulence generation at left and/or right boundary
     IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
 !
@@ -1313,8 +1337,8 @@
                 dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp ) *          &
                                     MERGE( 1.0_wp, 0.0_wp,                     &
-                                           BTEST( wall_flags_0(k,j,i), 1 ) )   
+                                           BTEST( wall_flags_0(k,j,i), 1 ) )
              ENDDO
           ENDDO
-          
+
           IF ( myidx == id_stg_left  )  i = nxl-1
           IF ( myidx == id_stg_right )  i = nxr+1
@@ -1341,5 +1365,4 @@
           ENDDO
        ENDIF
-
 !
 !--    Mass flux correction following Kim et al. (2013)
@@ -1514,5 +1537,5 @@
                                          + a22(k) * fv_xz(k,i) ), 3.0_wp ) *   &
                                     MERGE( 1.0_wp, 0.0_wp,                     &
-                                           BTEST( wall_flags_0(k,j,i), 2 ) )   
+                                           BTEST( wall_flags_0(k,j,i), 2 ) )
              ENDDO
           ENDDO
@@ -1648,14 +1671,18 @@
 !> parts are collected to form the full array.
 !------------------------------------------------------------------------------!
- SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id )
+ SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id_left, id_right )
+ 
+ USE pegrid
 
     IMPLICIT NONE
 
-    INTEGER(iwp) :: id          !< core ids at respective boundaries
-    INTEGER(iwp) :: j           !< loop index in y-direction
-    INTEGER(iwp) :: jj          !< loop index in y-direction
-    INTEGER(iwp) :: k           !< loop index in z-direction
-    INTEGER(iwp) :: kk          !< loop index in z-direction
-    INTEGER(iwp) :: send_count  !< send count for MPI_GATHERV
+    INTEGER(iwp)           :: i           !< grid index x-direction
+    INTEGER(iwp)           :: id_left     !< core ids at respective boundaries
+    INTEGER(iwp), OPTIONAL :: id_right    !< core ids at respective boundaries
+    INTEGER(iwp)           :: j           !< loop index in y-direction
+    INTEGER(iwp)           :: jj          !< loop index in y-direction
+    INTEGER(iwp)           :: k           !< loop index in z-direction
+    INTEGER(iwp)           :: kk          !< loop index in z-direction
+    INTEGER(iwp)           :: send_count  !< send count for MPI_GATHERV
 
     INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y    !< length scale in y-direction
@@ -1670,45 +1697,110 @@
     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_y     !< filter function in y-direction
     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter function in z-direction
+    
     REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nysg:nyng) :: f_n_l   !<  local velocity seed
-    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng)     :: f_n     !<  velocity seed
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE        :: rand_it !<  random number
-
-
-!
-!-- Generate random numbers using a seed generated in stg_init.
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng)             :: f_n     !<  velocity seed
+    
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it   !< global array of random numbers
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it_l !< local array of random numbers
+
+!
+!-- Generate random numbers using the parallel random generator.
 !-- The set of random numbers are modified to have an average of 0 and
-!-- unit variance.
+!-- unit variance. Note, at the end the random number array must be defined
+!-- globally in order to compute the correlation matrices. However,
+!-- the calculation and normalization of random numbers is done locally,
+!-- while the result is later distributed to all processes. Further, 
+!-- please note, a set of random numbers is only calculated for the 
+!-- left boundary, while the right boundary uses the same random numbers
+!-- and thus also computes the same correlation matrix.
     ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:ny+mergp) )
+    rand_it   = 0.0_wp
 
     rand_av        = 0.0_wp
     rand_sigma_inv = 0.0_wp
     nyz_inv        = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( ny+1 ), KIND=wp )
-
-    DO  j = 0, ny
-       DO  k = nzb, nzt+1
-          CALL RANDOM_NUMBER( rand_it(k,j) )
-          rand_av = rand_av + rand_it(k,j)
+!
+!-- Compute and normalize random numbers only on left-boundary ranks.
+    IF ( myidx == id_stg_left )  THEN
+!
+!--    Allocate array for local set of random numbers
+       ALLOCATE( rand_it_l(nzb:nzt+1,nys:nyn) )
+       rand_it_l = 0.0_wp
+
+       i = nxl
+       DO  j = nys, nyn
+!
+!--       Put the random seeds at grid point j,i
+          CALL random_seed_parallel( put=seq_random_array(:, j, i) )
+          DO  k = nzb, nzt+1
+             CALL random_number_parallel( random_dummy )
+             rand_it_l(k,j) = random_dummy
+             rand_av        = rand_av + rand_it_l(k,j)
+          ENDDO
+!
+!--       Get the new random seeds from last call at grid point j,i
+          CALL random_seed_parallel( get=seq_random_array(:, j, i) )
        ENDDO
-    ENDDO
-
-    rand_av = rand_av * nyz_inv
-
-    DO  j = 0, ny
-       DO  k = nzb, nzt+1
-          rand_it(k,j)   = rand_it(k,j) - rand_av
-          rand_sigma_inv = rand_sigma_inv + rand_it(k,j) ** 2
+#if defined( __parallel )
+!
+!--    Sum-up the local averages of the random numbers
+       CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_av, 1, MPI_REAL,                 &
+                           MPI_SUM, comm1dy, ierr )
+#endif
+       rand_av = rand_av * nyz_inv
+!
+!--    Obtain zero mean
+       rand_it_l = rand_it_l - rand_av
+!
+!--    Now, compute the variance
+       DO  j = nys, nyn
+          DO  k = nzb, nzt+1
+             rand_sigma_inv = rand_sigma_inv + rand_it_l(k,j)**2
+          ENDDO
        ENDDO
-    ENDDO
-
-    rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nyz_inv)
-
-    DO  j = 0, ny
-       DO  k = nzb, nzt+1
-          rand_it(k,j) = rand_it(k,j) * rand_sigma_inv
-       ENDDO
-    ENDDO
-
-!
-!-- Periodic fill of random number in space
+
+#if defined( __parallel )
+!
+!--    Sum-up the local quadratic averages of the random numbers
+       CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_sigma_inv, 1, MPI_REAL,          &
+                           MPI_SUM, comm1dy, ierr )
+#endif
+!
+!--    Compute standard deviation
+       IF ( rand_sigma_inv /= 0.0_wp )  THEN
+          rand_sigma_inv = 1.0_wp / SQRT( rand_sigma_inv * nyz_inv )
+       ELSE
+          rand_sigma_inv = 1.0_wp
+       ENDIF
+!
+!--    Normalize with standard deviation to obtain unit variance
+       rand_it_l = rand_it_l * rand_sigma_inv
+!
+!--    Copy local random numbers on the global array
+       rand_it(nzb:nzt+1,nys:nyn) = rand_it_l(nzb:nzt+1,nys:nyn)
+!
+!--    Deallocate local array
+       DEALLOCATE( rand_it_l )
+!
+!--    Now, distribute the final set of random numbers to all mpi ranks located
+!--    on the left boundary. Here, an allreduce with sum reduction is sufficient,
+!--    or, in the non-parallel case, nothing need to be done at all.
+#if defined( __parallel )
+       CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_it, SIZE( rand_it ), MPI_REAL,   &
+                           MPI_SUM, comm1dy, ierr )
+#endif
+    ENDIF
+!
+!-- Finally, distribute the set of random numbers (defined on the leftmost-
+!-- located mpi ranks) to all other mpi ranks. Here, a allreduce with sum 
+!-- option is sufficient, because rand_it is zero on all other mpi_ranks. 
+!-- Note, the reduce operation is only performed with communicator comm1dx, 
+!-- where only 1 rank within the communicator has non-zero random numbers.
+#if defined( __parallel )
+    CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_it, SIZE( rand_it ), MPI_REAL,      &
+                        MPI_SUM, comm1dx, ierr )
+#endif
+!
+!-- Periodic fill of random numbers in space
     DO  j = 0, ny
        DO  k = 1, mergp
@@ -1748,7 +1840,17 @@
 
 #if defined( __parallel )
-    CALL MPI_GATHERV( f_n_l(nzb_x_stg,nysg), send_count, stg_type_yz_small,    &
-                      f_n(nzb+1,nysg), recv_count_yz, displs_yz, stg_type_yz,  &
-                      id, comm1dx, ierr )
+!
+!-- Gather the velocity seed matrix on left boundary mpi ranks.
+    CALL MPI_GATHERV( f_n_l(nzb_x_stg,nysg), send_count, stg_type_yz_small,     &
+                      f_n(nzb+1,nysg), recv_count_yz, displs_yz, stg_type_yz,   &
+                      id_left, comm1dx, ierr )
+!
+!-- If required, gather the same velocity seed matrix on right boundary 
+!-- mpi ranks (in offline nesting for example).
+    IF ( PRESENT( id_right ) )  THEN
+       CALL MPI_GATHERV( f_n_l(nzb_x_stg,nysg), send_count, stg_type_yz_small,  &
+                         f_n(nzb+1,nysg), recv_count_yz, displs_yz, stg_type_yz,&
+                         id_right, comm1dx, ierr )
+    ENDIF
 #else
     f_n(nzb+1:nzt+1,nysg:nyng) = f_n_l(nzb_x_stg:nzt_x_stg+1,nysg:nyng)
@@ -1757,6 +1859,4 @@
 
  END SUBROUTINE stg_generate_seed_yz
-
-
 
 
@@ -1770,11 +1870,13 @@
 !> parts are collected to form the full array.
 !------------------------------------------------------------------------------!
- SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id )
+ SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id_south, id_north )
 
     IMPLICIT NONE
 
-    INTEGER(iwp) :: id          !< core ids at respective boundaries
     INTEGER(iwp) :: i           !< loop index in x-direction
+    INTEGER(iwp) :: id_north    !< core ids at respective boundaries
+    INTEGER(iwp) :: id_south    !< core ids at respective boundaries
     INTEGER(iwp) :: ii          !< loop index in x-direction
+    INTEGER(iwp) :: j           !< grid index y-direction
     INTEGER(iwp) :: k           !< loop index in z-direction
     INTEGER(iwp) :: kk          !< loop index in z-direction
@@ -1792,43 +1894,103 @@
     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_x     !< filter function in x-direction
     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter function in z-direction
+    
     REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg) :: f_n_l   !<  local velocity seed
-    REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg)     :: f_n     !<  velocity seed
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE        :: rand_it !<  random number
-
-
-!
-!-- Generate random numbers using a seed generated in stg_init.
+    REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg)             :: f_n     !<  velocity seed
+
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it   !< global array of random numbers
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it_l !< local array of random numbers
+
+!
+!-- Generate random numbers using the parallel random generator.
 !-- The set of random numbers are modified to have an average of 0 and
-!-- unit variance.
+!-- unit variance. Note, at the end the random number array must be defined
+!-- globally in order to compute the correlation matrices. However,
+!-- the calculation and normalization of random numbers is done locally,
+!-- while the result is later distributed to all processes. Further, 
+!-- please note, a set of random numbers is only calculated for the 
+!-- left boundary, while the right boundary uses the same random numbers
+!-- and thus also computes the same correlation matrix.
     ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:nx+mergp) )
+    rand_it   = 0.0_wp
 
     rand_av        = 0.0_wp
     rand_sigma_inv = 0.0_wp
     nxz_inv        = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( nx+1 ), KIND=wp )
-
-    DO  i = 0, nx
-       DO  k = nzb, nzt+1
-          CALL RANDOM_NUMBER( rand_it(k,i) )
-          rand_av = rand_av + rand_it(k,i)
+!
+!-- Compute and normalize random numbers only on south-boundary ranks.
+    IF ( myidy == id_stg_south )  THEN
+!
+!--    Allocate array for local set of random numbers
+       ALLOCATE( rand_it_l(nzb:nzt+1,nxl:nxr) )
+       rand_it_l = 0.0_wp
+
+       j = nys
+       DO  i = nxl, nxr
+!
+!--       Put the random seeds at grid point j,i
+          CALL random_seed_parallel( put=seq_random_array(:, j, i) )
+          DO  k = nzb, nzt+1
+             CALL random_number_parallel( random_dummy )
+             rand_it_l(k,i) = random_dummy
+             rand_av = rand_av + rand_it_l(k,i)
+          ENDDO
+!
+!--       Get the new random seeds from last call at grid point j,i
+          CALL random_seed_parallel( get=seq_random_array(:, j, i) )
        ENDDO
-    ENDDO
-
-    rand_av = rand_av * nxz_inv
-
-    DO  i = 0, nx
-       DO  k = nzb, nzt+1
-          rand_it(k,i)   = rand_it(k,i) - rand_av
-          rand_sigma_inv = rand_sigma_inv + rand_it(k,i) ** 2
+#if defined( __parallel )
+       CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_av, 1, MPI_REAL,                 &
+                           MPI_SUM, comm1dx, ierr )
+#endif
+       rand_av = rand_av * nxz_inv
+!
+!--    Obtain zero mean
+       rand_it_l = rand_it_l - rand_av
+!
+!--    Now, compute the variance
+       DO  i = nxl, nxr
+          DO  k = nzb, nzt+1
+             rand_sigma_inv = rand_sigma_inv + rand_it_l(k,i)**2
+          ENDDO
        ENDDO
-    ENDDO
-
-    rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nxz_inv)
-
-    DO  i = 0, nx
-       DO  k = nzb, nzt+1
-          rand_it(k,i) = rand_it(k,i) * rand_sigma_inv
-       ENDDO
-    ENDDO
-
+
+#if defined( __parallel )
+       CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_sigma_inv, 1, MPI_REAL,          &
+                           MPI_SUM, comm1dx, ierr )
+#endif
+!
+!--    Compute standard deviation
+       IF ( rand_sigma_inv /= 0.0_wp )  THEN
+          rand_sigma_inv = 1.0_wp / SQRT( rand_sigma_inv * nxz_inv )
+       ELSE
+          rand_sigma_inv = 1.0_wp
+       ENDIF
+!
+!--    Normalize with standard deviation to obtain unit variance
+       rand_it_l = rand_it_l * rand_sigma_inv
+!
+!--    Copy local random numbers on the global array
+       rand_it(nzb:nzt+1,nxl:nxr) = rand_it_l(nzb:nzt+1,nxl:nxr)
+!
+!--    Deallocate local array
+       DEALLOCATE( rand_it_l )
+!
+!--    Now, distribute the final set of random numbers to all mpi ranks located 
+!--    on the south boundary. Here, an allreduce with sum reduction is sufficient.
+#if defined( __parallel )
+       CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_it, SIZE( rand_it ), MPI_REAL,   &
+                           MPI_SUM, comm1dx, ierr )
+#endif
+    ENDIF
+!
+!-- Finally, distribute the set of random numbers (defined on the  
+!-- southmost-located mpi ranks) to all other mpi ranks. Here, a allreduce 
+!-- with sum option is sufficient, because rand_it is zero on all other 
+!-- mpi_ranks. Note, the reduce operation is only performed with communicator 
+!-- comm1dy, where only 1 rank within the communicator has non zero random numbers.
+#if defined( __parallel )
+    CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_it, SIZE( rand_it ), MPI_REAL,      &
+                        MPI_SUM, comm1dy, ierr )
+#endif
 !
 !-- Periodic fill of random number in space
@@ -1845,5 +2007,4 @@
        ENDDO
     ENDDO
-
 !
 !-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
@@ -1870,13 +2031,18 @@
     IF ( nzt_y_stg == nzt )  send_count = send_count + 1
 
-
 #if defined( __parallel )
-    CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxlg), send_count, stg_type_xz_small,    &
-                      f_n(nzb+1,nxlg), recv_count_xz, displs_xz, stg_type_xz,  &
-                      id, comm1dy, ierr )
+!
+!-- Gather the processed velocity seed on south boundary mpi ranks.
+    CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxlg), send_count, stg_type_xz_small,   &
+                      f_n(nzb+1,nxlg), recv_count_xz, displs_xz, stg_type_xz, &
+                      id_south, comm1dy, ierr )
+!
+!-- Gather the processed velocity seed on north boundary mpi ranks.
+    CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxlg), send_count, stg_type_xz_small,   &
+                      f_n(nzb+1,nxlg), recv_count_xz, displs_xz, stg_type_xz, &
+                      id_north, comm1dy, ierr )
 #else
     f_n(nzb+1:nzt+1,nxlg:nxrg) = f_n_l(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg)
 #endif
-
 
  END SUBROUTINE stg_generate_seed_xz
Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4307)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4309)
@@ -23,11 +23,14 @@
 ! Current revisions:
 ! ------------------
+! 
+! 
+! Former revisions:
+! -----------------
+! $Id$
 ! - Bugfix, include m_liq into restarts
 ! - Remove unused arrays for liquid water and saturation moisture at vertical
 !   walls
 !
-! Former revisions:
-! -----------------
-! $Id$
+! 4305 2019-11-25 11:15:40Z suehring
 ! Revision of some indoor-model parameters
 ! 
