Index: /palm/trunk/SOURCE/lagrangian_particle_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/lagrangian_particle_model_mod.f90	(revision 4544)
+++ /palm/trunk/SOURCE/lagrangian_particle_model_mod.f90	(revision 4545)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Using parallel random generator, thus preventing dependency of PE number
+! 
+! 4535 2020-05-15 12:07:23Z raasch
 ! bugfix for restart data format query
 ! 
@@ -212,6 +215,11 @@
                netcdf_handle_error
 
-    USE random_function_mod,                                                   &
-        ONLY:  random_function
+    USE random_generator_parallel,                                             &
+        ONLY:  init_parallel_random_generator,                                 &
+               random_dummy,                                                   &
+               random_number_parallel,                                         &
+               random_number_parallel_gauss,                                   &
+               random_seed_parallel,                                           &
+               id_random_array
 
     USE restart_data_mpi_io_mod,                                                                   &
@@ -273,4 +281,6 @@
     INTEGER(iwp) ::  trnp_count_sum                               !< parameter for particle exchange of PEs
     INTEGER(iwp) ::  trnp_count_recv_sum                          !< parameter for particle exchange of PEs
+
+    INTEGER(isp), DIMENSION(:,:,:), ALLOCATABLE ::  seq_random_array_particle   !< sequence of random array for particle
 
     LOGICAL ::  lagrangian_particle_model = .FALSE.       !< namelist parameter (see documentation) 
@@ -1244,4 +1254,22 @@
 !--    different on the different PEs.
        iran_part = iran_part + myid
+
+!
+!--    Initialize parallel random number sequence seed for particles
+!--    This is done individually, as thus particle random numbers does
+!--    not affect random numbers used for the flow field.
+       ALLOCATE ( seq_random_array_particle(5,nys:nyn,nxl:nxr) )
+       seq_random_array_particle = 0
+
+!--    Initializing with random_seed_parallel for every vertical
+!--    gridpoint column.
+       random_dummy = 0
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             CALL random_seed_parallel (random_sequence=id_random_array(j, i))
+             CALL random_number_parallel (random_dummy)
+             CALL random_seed_parallel (get=seq_random_array_particle(:, j, i))
+          ENDDO
+       ENDDO
 !
 !--    Create the particle set, and set the initial particles
@@ -1507,4 +1535,7 @@
        DO  ip = nxl, nxr
           DO  jp = nys, nyn
+!
+!--          Put the random seeds at grid point jp, ip
+             CALL random_seed_parallel( put=seq_random_array_particle(:,jp,ip) )
              DO  kp = nzb+1, nzt
                 number_of_particles = prt_count(kp,jp,ip)
@@ -1519,5 +1550,6 @@
                 DO  n = local_start(kp,jp,ip), number_of_particles
                    IF ( psl(particles(n)%group) /= psr(particles(n)%group) )  THEN
-                      rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
+                      CALL random_number_parallel( random_dummy )
+                      rand_contr = ( random_dummy - 0.5_wp ) *                 &
                                      pdx(particles(n)%group)
                       particles(n)%x = particles(n)%x +                        &
@@ -1527,5 +1559,6 @@
                    ENDIF
                    IF ( pss(particles(n)%group) /= psn(particles(n)%group) )  THEN
-                      rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
+                      CALL random_number_parallel( random_dummy )
+                      rand_contr = ( random_dummy - 0.5_wp ) *                 &
                                      pdy(particles(n)%group)
                       particles(n)%y = particles(n)%y +                        &
@@ -1535,5 +1568,6 @@
                    ENDIF
                    IF ( psb(particles(n)%group) /= pst(particles(n)%group) )  THEN
-                      rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
+                      CALL random_number_parallel( random_dummy )
+                      rand_contr = ( random_dummy - 0.5_wp ) *                 &
                                      pdz(particles(n)%group)
                       particles(n)%z = particles(n)%z +                        &
@@ -1572,4 +1606,7 @@
                 ENDDO
              ENDDO
+!
+!--       Get the new random seeds from last call at grid point jp, ip
+          CALL random_seed_parallel( get=seq_random_array_particle(:,jp,ip) )
           ENDDO
        ENDDO
@@ -1699,4 +1736,7 @@
     DO  ip = nxl, nxr
        DO  jp = nys, nyn
+!
+!--       Put the random seeds at grid point jp, ip
+          CALL random_seed_parallel( put=seq_random_array_particle(:,jp,ip) )
           DO  kp = nzb+1, nzt
 
@@ -1730,7 +1770,9 @@
 !--             to increase or decrease the number of simulated aerosols
                 particles(n)%weight_factor = particles(n)%weight_factor * aero_weight
-
+!
+!--             create random numver with parallel number generator
+                CALL random_number_parallel( random_dummy )
                 IF ( particles(n)%weight_factor - FLOOR(particles(n)%weight_factor,KIND=wp) &
-                     > random_function( iran_part ) )  THEN
+                     > random_dummy )  THEN
                    particles(n)%weight_factor = FLOOR(particles(n)%weight_factor,KIND=wp) + 1.0_wp
                 ELSE
@@ -1776,4 +1818,7 @@
 
           ENDDO
+!
+!--    Get the new random seeds from last call at grid point j
+       CALL random_seed_parallel( get=seq_random_array_particle(:,jp,ip) )
        ENDDO
     ENDDO
@@ -2168,4 +2213,8 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
+!
+!--                   Put the random seeds at grid point j, i
+                      CALL random_seed_parallel( put=seq_random_array_particle(:,j,i) )
+
                       DO  k = nzb+1, nzt
 
@@ -2250,4 +2299,8 @@
 
                       ENDDO
+!
+!--                   Get the new random seeds from last call at grid point jp, ip
+                      CALL random_seed_parallel( get=seq_random_array_particle(:,j,i) )
+
                    ENDDO
                 ENDDO
@@ -3394,5 +3447,4 @@
     REAL(wp) ::  height_p           !< dummy argument for logarithmic interpolation
     REAL(wp) ::  log_z_z0_int       !< logarithmus used for surface_layer interpolation
-    REAL(wp) ::  random_gauss       !< Gaussian-distributed random number used for SGS particle advection
     REAL(wp) ::  RL                 !< Lagrangian autocorrelation coefficient
     REAL(wp) ::  rg1                !< Gaussian distributed random number
@@ -4010,7 +4062,10 @@
        DO  nb = 0, 7
           DO  n = start_index(nb), end_index(nb)
-             rg(n,1) = random_gauss( iran_part, 5.0_wp )
-             rg(n,2) = random_gauss( iran_part, 5.0_wp )
-             rg(n,3) = random_gauss( iran_part, 5.0_wp )
+             CALL random_number_parallel_gauss( random_dummy )
+             rg(n,1) = random_dummy
+             CALL random_number_parallel_gauss( random_dummy )
+             rg(n,2) = random_dummy
+             CALL random_number_parallel_gauss( random_dummy )
+             rg(n,3) = random_dummy
           ENDDO
        ENDDO
@@ -4054,9 +4109,9 @@
 !--             from becoming unrealistically large.
                 rvar1_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n)          &
-                                          + 1E-20_wp ) * ( rg(n,1) - 1.0_wp )
+                                          + 1E-20_wp ) * rg(n,1)
                 rvar2_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n)          &
-                                          + 1E-20_wp ) * ( rg(n,2) - 1.0_wp )
+                                          + 1E-20_wp ) * rg(n,2)
                 rvar3_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n)          &
-                                          + 1E-20_wp ) * ( rg(n,3) - 1.0_wp )
+                                          + 1E-20_wp ) * rg(n,3)
              ELSE
 !
@@ -4244,8 +4299,13 @@
                                              1.0E-20_wp ) )
                       sigma          = SQRT( e(kp,jp,ip) )
-
-                      rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
-                      rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
-                      rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
+!
+!--                   Calculate random component of particle sgs velocity using parallel
+!--                   random generator
+                      CALL random_number_parallel_gauss( random_dummy )
+                      rg1 = random_dummy
+                      CALL random_number_parallel_gauss( random_dummy )
+                      rg2 = random_dummy
+                      CALL random_number_parallel_gauss( random_dummy )
+                      rg3 = random_dummy
 
                       particles(n)%rvar1 = RL * particles(n)%rvar1 +              &
@@ -4340,7 +4400,13 @@
                     sigma          = SQRT( e(kp,jp,ip) )
 
-                    rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
-                    rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
-                    rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
+!
+!--                 Calculate random component of particle sgs velocity using parallel
+!--                 random generator
+                    CALL random_number_parallel_gauss( random_dummy )
+                    rg1 = random_dummy
+                    CALL random_number_parallel_gauss( random_dummy )
+                    rg2 = random_dummy
+                    CALL random_number_parallel_gauss( random_dummy )
+                    rg3 = random_dummy
 
                     particles(n)%rvar1 = RL * particles(n)%rvar1 +                &
@@ -4603,5 +4669,4 @@
     REAL(wp) ::  prt_y          !< current particle position (y)
     REAL(wp) ::  prt_z          !< current particle position (z)
-    REAL(wp) ::  ran_val        !< location of wall in z
     REAL(wp) ::  reset_top      !< location of wall in z
     REAL(wp) ::  t_old          !< previous reflection time
@@ -4682,7 +4747,6 @@
                 index_reset = MINLOC( prt_count(nzb+1:particles_top,j,i), DIM = 1 )
                 reset_top = zu(index_reset)
-                iran_part = iran_part + myid
-                ran_val = random_function( iran_part )
-                particles(n)%z       = reset_top *  ( 1.0  + ( ran_val / 10.0_wp) )
+                CALL random_number_parallel( random_dummy )
+                particles(n)%z       = reset_top *  ( 1.0  + ( random_dummy / 10.0_wp) )
                 particles(n)%speed_z = 0.0_wp
                 IF ( curvature_solution_effects )  THEN
@@ -5582,6 +5646,7 @@
 !--          Calculate the number of collections and consider multiple collections.
 !--          (Accordingly, p_crit will be 0.0, 1.0, 2.0, ...)
+             CALL random_number_parallel( random_dummy )
              IF ( collection_probability - FLOOR(collection_probability)    &
-                  > random_function( iran_part ) )  THEN
+                  > random_dummy )  THEN
                 collection_probability = FLOOR(collection_probability) + 1.0_wp
              ELSE
Index: /palm/trunk/SOURCE/random_generator_parallel_mod.f90
===================================================================
--- /palm/trunk/SOURCE/random_generator_parallel_mod.f90	(revision 4544)
+++ /palm/trunk/SOURCE/random_generator_parallel_mod.f90	(revision 4545)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Add generator (using parallel mode) returning gaussian distributed random
+! number
+! 
+! 4438 2020-03-03 20:49:28Z raasch
 ! - Rename variables to avoid confusion with subdomain grid indices
 ! - Some formatting adjustments
@@ -45,4 +49,6 @@
 !> ran_parallel returns a uniform random deviate between 0.0 and 1.0 
 !> (exclusive of the end point values).
+!> Moreover, it contains a routine returning a normally distributed random number
+!> with mean zero and unity standard deviation.
 !> Additionally it provides the generator with five integer for use as initial state space.
 !> The first tree integers (iran, jran, kran) are maintained as non negative values, 
@@ -95,4 +101,8 @@
    END INTERFACE
 
+   INTERFACE random_number_parallel_gauss
+      MODULE PROCEDURE gasdev_s
+   END INTERFACE
+
    INTERFACE random_seed_parallel
       MODULE PROCEDURE random_seed_parallel
@@ -114,5 +124,6 @@
 
    PUBLIC random_number_parallel, random_seed_parallel, random_dummy,          &
-          id_random_array, seq_random_array, init_parallel_random_generator
+          id_random_array, seq_random_array, init_parallel_random_generator,   &
+          random_number_parallel_gauss
 
  CONTAINS
@@ -248,4 +259,54 @@
 ! Description:
 ! ------------
+!> Returns in harvest a normally distributed deviate with zero mean and unit
+!> variance, using ran0_s as the source of uniform deviates. Following
+!> Numerical Recipes in Fortran90 (Press et al., 2nd edition, 1996, pp 1152ff).
+!> Note, instead of ran1_s, ran0_s is used.
+!------------------------------------------------------------------------------!
+   SUBROUTINE gasdev_s(harvest)
+
+      REAL(wp), INTENT(OUT) ::  harvest   !<
+
+      REAL(wp) ::  rsq      !<
+      REAL(wp) ::  v1       !<
+      REAL(wp) ::  v2       !<
+      REAL(wp), SAVE ::  g  !<
+
+      LOGICAL, SAVE ::  gaus_stored = .FALSE. !<
+!
+!--   We have an extra deviate handy, so return it, and unset the flag.
+      IF (gaus_stored)  THEN
+         harvest = g
+         gaus_stored = .FALSE.
+!
+!--   We donât have an extra deviate handy, so pick two uniform numbers in the
+!--   square extending from -1 to +1 in each direction
+      ELSE
+         DO
+            CALL ran0_s(v1)
+            CALL ran0_s(v2)
+            v1 = 2.0_wp * v1 - 1.0_wp
+            v2 = 2.0_wp * v2 - 1.0_wp
+!
+!--         see if they are in the unit circle
+            rsq = v1**2 + v2**2
+!
+!--         otherwise try again.
+            IF (rsq > 0.0  .AND.  rsq < 1.0) EXIT
+         ENDDO
+!
+!--      Now make the Box-Muller transformation to get two normal deviates.
+!--      Return one and save the other for next time. Set flag.
+         rsq = SQRT(-2.0_sp * LOG(rsq)/rsq)
+         harvest = v1 * rsq
+         g = v2 * rsq
+         gaus_stored = .TRUE.
+      ENDIF
+
+   END SUBROUTINE gasdev_s
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Initialize or reinitialize the random generator state space to vectors of size length. 
 !> The saved variable seq is hashed (via calls to the module routine ran_hash) 
