Index: palm/trunk/SOURCE/lpm_advec.f90
===================================================================
--- palm/trunk/SOURCE/lpm_advec.f90	(revision 2885)
+++ palm/trunk/SOURCE/lpm_advec.f90	(revision 2886)
@@ -25,4 +25,10 @@
 ! -----------------
 ! $Id$
+! Bugfix in passive particle SGS Model:
+! Sometimes the added SGS velocities would lead to a violation of the CFL 
+! criterion for single particles. For this a check was added after the 
+! calculation of SGS velocities.
+! 
+! 2718 2018-01-02 08:49:38Z maronga
 ! Corrected "Former revisions" section
 ! 
@@ -218,6 +224,6 @@
     REAL(wp) ::  diss_int_l         !< x/y-interpolated dissipation at particle position at lower vertical level
     REAL(wp) ::  diss_int_u         !< x/y-interpolated dissipation at particle position at upper vertical level
-    REAL(wp) ::  dt_gap             !< remaining time until particle time integration reaches LES time
     REAL(wp) ::  dt_particle_m      !< previous particle time step
+    REAL(wp) ::  dz_temp            !< 
     REAL(wp) ::  e_int_l            !< x/y-interpolated TKE at particle position at lower vertical level
     REAL(wp) ::  e_int_u            !< x/y-interpolated TKE at particle position at upper vertical level
@@ -227,5 +233,4 @@
     REAL(wp) ::  gg                 !< dummy argument for horizontal particle interpolation 
     REAL(wp) ::  height_p           !< dummy argument for logarithmic interpolation
-    REAL(wp) ::  lagr_timescale     !< Lagrangian timescale
     REAL(wp) ::  location(1:30,1:3) !< wall locations
     REAL(wp) ::  log_z_z0_int       !< logarithmus used for surface_layer interpolation
@@ -265,19 +270,24 @@
     REAL(wp), DIMENSION(1:30) ::  ei      !< TKE at adjacent wall
 
-    REAL(wp), DIMENSION(number_of_particles) ::  term_1_2     !< flag to communicate whether a particle is near topography or not
-    REAL(wp), DIMENSION(number_of_particles) ::  dens_ratio   !<
-    REAL(wp), DIMENSION(number_of_particles) ::  de_dx_int    !< horizontal TKE gradient along x at particle position
-    REAL(wp), DIMENSION(number_of_particles) ::  de_dy_int    !< horizontal TKE gradient along y at particle position
-    REAL(wp), DIMENSION(number_of_particles) ::  de_dz_int    !< horizontal TKE gradient along z at particle position
-    REAL(wp), DIMENSION(number_of_particles) ::  diss_int     !< dissipation at particle position
-    REAL(wp), DIMENSION(number_of_particles) ::  dt_particle  !< particle time step
-    REAL(wp), DIMENSION(number_of_particles) ::  e_int        !< TKE at particle position
-    REAL(wp), DIMENSION(number_of_particles) ::  fs_int       !< weighting factor for subgrid-scale particle speed
-    REAL(wp), DIMENSION(number_of_particles) ::  u_int        !< u-component of particle speed
-    REAL(wp), DIMENSION(number_of_particles) ::  v_int        !< v-component of particle speed 
-    REAL(wp), DIMENSION(number_of_particles) ::  w_int        !< w-component of particle speed
-    REAL(wp), DIMENSION(number_of_particles) ::  xv           !< x-position
-    REAL(wp), DIMENSION(number_of_particles) ::  yv           !< y-position
-    REAL(wp), DIMENSION(number_of_particles) ::  zv           !< z-position
+    REAL(wp), DIMENSION(number_of_particles) ::  term_1_2       !< flag to communicate whether a particle is near topography or not
+    REAL(wp), DIMENSION(number_of_particles) ::  dens_ratio     !<
+    REAL(wp), DIMENSION(number_of_particles) ::  de_dx_int      !< horizontal TKE gradient along x at particle position
+    REAL(wp), DIMENSION(number_of_particles) ::  de_dy_int      !< horizontal TKE gradient along y at particle position
+    REAL(wp), DIMENSION(number_of_particles) ::  de_dz_int      !< horizontal TKE gradient along z at particle position
+    REAL(wp), DIMENSION(number_of_particles) ::  diss_int       !< dissipation at particle position
+    REAL(wp), DIMENSION(number_of_particles) ::  dt_gap         !< remaining time until particle time integration reaches LES time
+    REAL(wp), DIMENSION(number_of_particles) ::  dt_particle    !< particle time step
+    REAL(wp), DIMENSION(number_of_particles) ::  e_int          !< TKE at particle position
+    REAL(wp), DIMENSION(number_of_particles) ::  fs_int         !< weighting factor for subgrid-scale particle speed
+    REAL(wp), DIMENSION(number_of_particles) ::  lagr_timescale !< Lagrangian timescale
+    REAL(wp), DIMENSION(number_of_particles) ::  rvar1_temp     !< 
+    REAL(wp), DIMENSION(number_of_particles) ::  rvar2_temp     !< 
+    REAL(wp), DIMENSION(number_of_particles) ::  rvar3_temp     !< 
+    REAL(wp), DIMENSION(number_of_particles) ::  u_int          !< u-component of particle speed
+    REAL(wp), DIMENSION(number_of_particles) ::  v_int          !< v-component of particle speed 
+    REAL(wp), DIMENSION(number_of_particles) ::  w_int          !< w-component of particle speed
+    REAL(wp), DIMENSION(number_of_particles) ::  xv             !< x-position
+    REAL(wp), DIMENSION(number_of_particles) ::  yv             !< y-position
+    REAL(wp), DIMENSION(number_of_particles) ::  zv             !< z-position
 
     REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !< vector of Gaussian distributed random numbers 
@@ -763,5 +773,5 @@
 !
 !--          Calculate the Lagrangian timescale according to Weil et al. (2004).
-             lagr_timescale = ( 4.0_wp * e_int(n) + 1E-20_wp ) / &
+             lagr_timescale(n) = ( 4.0_wp * e_int(n) + 1E-20_wp ) / &
                               ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp )
 
@@ -769,14 +779,17 @@
 !--          Calculate the next particle timestep. dt_gap is the time needed to
 !--          complete the current LES timestep.
-             dt_gap = dt_3d - particles(n)%dt_sum
-             dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale, dt_gap )
-
+             dt_gap(n) = dt_3d - particles(n)%dt_sum
+             dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale(n), dt_gap(n) )
+             particles(n)%aux1 = lagr_timescale(n)
+             particles(n)%aux2 = dt_gap(n)
 !
 !--          The particle timestep should not be too small in order to prevent
 !--          the number of particle timesteps of getting too large
-             IF ( dt_particle(n) < dt_min_part  .AND.  dt_min_part < dt_gap )  THEN
+             IF ( dt_particle(n) < dt_min_part  .AND.  dt_min_part < dt_gap(n) )  THEN
                 dt_particle(n) = dt_min_part
              ENDIF
-
+             rvar1_temp(n) = particles(n)%rvar1
+             rvar2_temp(n) = particles(n)%rvar2
+             rvar3_temp(n) = particles(n)%rvar3
 !
 !--          Calculate the SGS velocity components
@@ -787,13 +800,10 @@
 !--             [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities
 !--             from becoming unrealistically large.
-                particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf_part * e_int(n)     &
-                                          + 1E-20_wp ) *                       &
-                                     ( rg(n,1) - 1.0_wp )
-                particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf_part * e_int(n)     &
-                                          + 1E-20_wp ) *                       &
-                                     ( rg(n,2) - 1.0_wp )
-                particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf_part * e_int(n)     &
-                                          + 1E-20_wp ) *                       &
-                                     ( rg(n,3) - 1.0_wp )
+                rvar1_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n)          &
+                                          + 1E-20_wp ) * ( rg(n,1) - 1.0_wp )
+                rvar2_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n)          &
+                                          + 1E-20_wp ) * ( rg(n,2) - 1.0_wp )
+                rvar3_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n)          &
+                                          + 1E-20_wp ) * ( rg(n,3) - 1.0_wp )
 
              ELSE
@@ -811,5 +821,4 @@
                 ENDIF
 
-!
 !--             For old particles the SGS components are correlated with the
 !--             values from the previous timestep. Random numbers have also to
@@ -827,13 +836,13 @@
                 ENDIF
 
-                CALL weil_stochastic_eq(particles(n)%rvar1, fs_int(n), e_int(n),& 
+                CALL weil_stochastic_eq(rvar1_temp(n), fs_int(n), e_int(n),& 
                                         de_dx_int(n), de_dt, diss_int(n),       &
                                         dt_particle(n), rg(n,1), term_1_2(n) )
 
-                CALL weil_stochastic_eq(particles(n)%rvar2, fs_int(n), e_int(n),& 
+                CALL weil_stochastic_eq(rvar2_temp(n), fs_int(n), e_int(n),& 
                                         de_dy_int(n), de_dt, diss_int(n),       &
                                         dt_particle(n), rg(n,2), term_1_2(n) )
 
-                CALL weil_stochastic_eq(particles(n)%rvar3, fs_int(n), e_int(n),& 
+                CALL weil_stochastic_eq(rvar3_temp(n), fs_int(n), e_int(n),& 
                                         de_dz_int(n), de_dt, diss_int(n),       &
                                         dt_particle(n), rg(n,3), term_1_2(n) )
@@ -841,4 +850,56 @@
              ENDIF
 
+          ENDDO
+       ENDDO
+!
+!--    Check if the added SGS velocities result in a violation of the CFL-
+!--    criterion. If yes choose a smaller timestep based on the new velocities
+!--    and calculate SGS velocities again
+       dz_temp = zw(kp)-zw(kp-1)
+       
+       DO  nb = 0, 7
+          DO  n = start_index(nb), end_index(nb)
+             IF ( .NOT. particles(n)%age == 0.0_wp .AND.                       &
+                (ABS( u_int(n) + rvar1_temp(n) ) > (dx/dt_particle(n))  .OR.   &
+                 ABS( v_int(n) + rvar2_temp(n) ) > (dy/dt_particle(n))  .OR.   &
+                 ABS( w_int(n) + rvar3_temp(n) ) > (dz_temp/dt_particle(n)))) THEN
+                
+                dt_particle(n) = 0.9_wp * MIN(                                 &
+                                 ( dx / ABS( u_int(n) + rvar1_temp(n) ) ),     &
+                                 ( dy / ABS( v_int(n) + rvar2_temp(n) ) ),     &
+                                 ( dz_temp / ABS( w_int(n) + rvar3_temp(n) ) ) )
+
+!
+!--             Reset temporary SGS velocites to "current" ones
+                rvar1_temp(n) = particles(n)%rvar1
+                rvar2_temp(n) = particles(n)%rvar2
+                rvar3_temp(n) = particles(n)%rvar3
+                
+                de_dt_min = - e_int(n) / dt_particle(n)
+
+                de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m
+
+                IF ( de_dt < de_dt_min )  THEN
+                   de_dt = de_dt_min
+                ENDIF
+
+                CALL weil_stochastic_eq(rvar1_temp(n), fs_int(n), e_int(n),& 
+                                        de_dx_int(n), de_dt, diss_int(n),       &
+                                        dt_particle(n), rg(n,1), term_1_2(n) )
+
+                CALL weil_stochastic_eq(rvar2_temp(n), fs_int(n), e_int(n),& 
+                                        de_dy_int(n), de_dt, diss_int(n),       &
+                                        dt_particle(n), rg(n,2), term_1_2(n) )
+
+                CALL weil_stochastic_eq(rvar3_temp(n), fs_int(n), e_int(n),& 
+                                        de_dz_int(n), de_dt, diss_int(n),       &
+                                        dt_particle(n), rg(n,3), term_1_2(n) )
+             ENDIF                           
+             
+!
+!--          Update particle velocites 
+             particles(n)%rvar1 = rvar1_temp(n)
+             particles(n)%rvar2 = rvar2_temp(n)
+             particles(n)%rvar3 = rvar3_temp(n)
              u_int(n) = u_int(n) + particles(n)%rvar1
              v_int(n) = v_int(n) + particles(n)%rvar2
@@ -850,5 +911,5 @@
           ENDDO
        ENDDO
-
+       
     ELSE
 !
@@ -905,6 +966,6 @@
 !--                (2010, Q. J. R. Meteorol. Soc.)
                    IF ( use_sgs_for_particles )  THEN
-                      lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
-                      RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
+                      lagr_timescale(n) = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
+                      RL             = EXP( -1.0_wp * dt_3d / lagr_timescale(n) )
                       sigma          = SQRT( e(kp,jp,ip) )
 
@@ -977,6 +1038,6 @@
 !--             (2010, Q. J. R. Meteorol. Soc.)
                 IF ( use_sgs_for_particles )  THEN
-                    lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
-                    RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
+                    lagr_timescale(n) = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
+                    RL             = EXP( -1.0_wp * dt_3d / lagr_timescale(n) )
                     sigma          = SQRT( e(kp,jp,ip) )
 
