Index: /palm/trunk/SCRIPTS/.palm.iofiles
===================================================================
--- /palm/trunk/SCRIPTS/.palm.iofiles	(revision 4429)
+++ /palm/trunk/SCRIPTS/.palm.iofiles	(revision 4430)
@@ -67,5 +67,5 @@
 DATA_1D_SP_NETCDF*         out:tr   *         $base_data/$run_identifier/OUTPUT         _sp        nc
 DATA_1D_TS_NETCDF*         out:tr   *         $base_data/$run_identifier/OUTPUT         _ts        nc
-DATA_1D_PTS_NETCDF         out:tr   *         $base_data/$run_identifier/OUTPUT         _pts       nc
+DATA_1D_PTS_NETCDF*        out:tr   *         $base_data/$run_identifier/OUTPUT         _pts       nc
 DATA_2D_XY_NETCDF*         out:tr   *         $base_data/$run_identifier/OUTPUT         _xy        nc
 DATA_2D_XY_AV_NETCDF*      out:tr   *         $base_data/$run_identifier/OUTPUT         _av_xy     nc
Index: /palm/trunk/SOURCE/lagrangian_particle_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/lagrangian_particle_model_mod.f90	(revision 4429)
+++ /palm/trunk/SOURCE/lagrangian_particle_model_mod.f90	(revision 4430)
@@ -25,4 +25,12 @@
 ! -----------------
 ! $Id$
+! - Bugfix in logarithmic interpolation of near-ground particle speed (density
+!   was not considered).
+! - Revise CFL-check when SGS particle speeds are considered.
+! - In nested case with SGS particle speeds in the child domain, do not give 
+!   warning that particles are on domain boundaries. At the end of the particle
+!   time integration these will be transferred to the parent domain anyhow.
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Introduction of wall_flags_total_0, which currently sets bits based on static
 ! topography information used in wall_flags_static_0
@@ -126,7 +134,9 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  de_dx, de_dy, de_dz, dzw, zu, zw,  ql_c, ql_v, ql_vp, hyp,      &
-               pt, q, exner, ql, diss, e, u, v, w, km, ql_1, ql_2, pt_p, q_p,  &
-               d_exner
+        ONLY:  de_dx, de_dy, de_dz,                                            &
+               d_exner,                                                        &
+               drho_air_zw,                                                    &
+               dzw, zu, zw,  ql_c, ql_v, ql_vp, hyp,                           &
+               pt, q, exner, ql, diss, e, u, v, w, km, ql_1, ql_2, pt_p, q_p
  
     USE averaging,                                                             &
@@ -139,4 +149,5 @@
     USE control_parameters,                                                    &
         ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
+               child_domain,                                                   &
                cloud_droplets, constant_flux_layer, current_timestep_number,   &
                dt_3d, dt_3d_reached, first_call_lpm, humidity,                 &
@@ -830,15 +841,14 @@
 !
 !--    Liquid water content, change in liquid water content
-       ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                      &
+       ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                         &
                   ql_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-!
 !--    Real volume of particles (with weighting), volume of particles
-       ALLOCATE ( ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                      &
-                     ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                         &
+                  ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     ENDIF
 
 
-    ALLOCATE( u_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
-              v_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
+    ALLOCATE( u_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
+              v_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
               w_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 !
@@ -2338,5 +2348,4 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE particles_from_parent_to_child
-    IMPLICIT NONE
 
     CALL pmcp_c_get_particle_from_parent                         ! Child actions
@@ -2344,4 +2353,5 @@
 
     RETURN
+
  END SUBROUTINE particles_from_parent_to_child
 
@@ -2353,5 +2363,4 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE particles_from_child_to_parent
-    IMPLICIT NONE
 
     CALL pmcp_c_send_particle_to_parent                         ! Child actions
@@ -2359,4 +2368,5 @@
 
     RETURN
+
  END SUBROUTINE particles_from_child_to_parent
  
@@ -3315,5 +3325,5 @@
     REAL(wp), DIMENSION(number_of_particles) ::  zv             !< z-position
 
-    REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !< vector of Gaussian distributed random numbers
+    REAL(wp), DIMENSION(number_of_particles, 3) ::  rg          !< vector of Gaussian distributed random numbers
 
     CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' )
@@ -3461,19 +3471,19 @@
 !
 !--                Determine the sublayer. Further used as index.
-                   height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
-                                        * REAL( number_of_sublayers, KIND=wp )    &
+                   height_p = ( zv(n) - zw(k_wall) - z0_av_global )            &
+                                        * REAL( number_of_sublayers, KIND=wp ) &
                                         * d_z_p_z0
 !
 !--                Calculate LOG(z/z0) for exact particle height. Therefore,
 !--                interpolate linearly between precalculated logarithm.
-                   log_z_z0_int = log_z_z0(INT(height_p))                         &
-                                    + ( height_p - INT(height_p) )                &
-                                    * ( log_z_z0(INT(height_p)+1)                 &
-                                         - log_z_z0(INT(height_p))                &
+                   log_z_z0_int = log_z_z0(INT(height_p))                      &
+                                    + ( height_p - INT(height_p) )             &
+                                    * ( log_z_z0(INT(height_p)+1)              &
+                                         - log_z_z0(INT(height_p))             &
                                       )
 !
 !--                Get friction velocity and momentum flux from new surface data
 !--                types.
-                   IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
+                   IF ( surf_def_h(0)%start_index(jp,ip) <=                    &
                         surf_def_h(0)%end_index(jp,ip) )  THEN
                       surf_start = surf_def_h(0)%start_index(jp,ip)
@@ -3482,15 +3492,18 @@
 !--                   large particle speed.
                       us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
-                      usws_int  = surf_def_h(0)%usws(surf_start)
-                   ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                  &
+                      usws_int  = surf_def_h(0)%usws(surf_start)               &
+                                * drho_air_zw(k_wall)
+                   ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                   &
                             surf_lsm_h%end_index(jp,ip) )  THEN
                       surf_start = surf_lsm_h%start_index(jp,ip)
                       us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
-                      usws_int  = surf_lsm_h%usws(surf_start)
-                   ELSEIF ( surf_usm_h%start_index(jp,ip) <=                  &
+                      usws_int  = surf_lsm_h%usws(surf_start)                  &
+                                * drho_air_zw(k_wall)
+                   ELSEIF ( surf_usm_h%start_index(jp,ip) <=                   &
                             surf_usm_h%end_index(jp,ip) )  THEN
                       surf_start = surf_usm_h%start_index(jp,ip)
                       us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp )
-                      usws_int  = surf_usm_h%usws(surf_start)
+                      usws_int  = surf_usm_h%usws(surf_start)                  &
+                                * drho_air_zw(k_wall)
                    ENDIF
 !
@@ -3501,5 +3514,5 @@
 !--                as sensitivity studies revealed no significant effect of
 !--                using the neutral solution also for un/stable situations.
-                   u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp )           &
+                   u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp )        &
                                * log_z_z0_int - u_gtrans
                 ENDIF
@@ -3553,19 +3566,19 @@
 !--                topography particle on u-grid can be above surface-layer height,
 !--                whereas it can be below on v-grid.
-                   height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
-                                     * REAL( number_of_sublayers, KIND=wp )       &
+                   height_p = ( zv(n) - zw(k_wall) - z0_av_global )            &
+                                     * REAL( number_of_sublayers, KIND=wp )    &
                                      * d_z_p_z0
 !
 !--                Calculate LOG(z/z0) for exact particle height. Therefore,
 !--                interpolate linearly between precalculated logarithm.
-                   log_z_z0_int = log_z_z0(INT(height_p))                         &
-                                    + ( height_p - INT(height_p) )                &
-                                    * ( log_z_z0(INT(height_p)+1)                 &
-                                         - log_z_z0(INT(height_p))                &
+                   log_z_z0_int = log_z_z0(INT(height_p))                      &
+                                    + ( height_p - INT(height_p) )             &
+                                    * ( log_z_z0(INT(height_p)+1)              &
+                                         - log_z_z0(INT(height_p))             &
                                       )
 !
 !--                Get friction velocity and momentum flux from new surface data
 !--                types.
-                   IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
+                   IF ( surf_def_h(0)%start_index(jp,ip) <=                    &
                         surf_def_h(0)%end_index(jp,ip) )  THEN
                       surf_start = surf_def_h(0)%start_index(jp,ip)
@@ -3574,15 +3587,18 @@
 !--                   large particle speed.
                       us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
-                      vsws_int  = surf_def_h(0)%vsws(surf_start)
-                   ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                  &
+                      vsws_int  = surf_def_h(0)%vsws(surf_start)               &
+                                * drho_air_zw(k_wall)
+                   ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                   &
                             surf_lsm_h%end_index(jp,ip) )  THEN
                       surf_start = surf_lsm_h%start_index(jp,ip)
                       us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
-                      vsws_int  = surf_lsm_h%vsws(surf_start)
-                   ELSEIF ( surf_usm_h%start_index(jp,ip) <=                  &
+                      vsws_int  = surf_lsm_h%vsws(surf_start)                  &
+                                * drho_air_zw(k_wall)
+                   ELSEIF ( surf_usm_h%start_index(jp,ip) <=                   &
                             surf_usm_h%end_index(jp,ip) )  THEN
                       surf_start = surf_usm_h%start_index(jp,ip)
                       us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp )
-                      vsws_int  = surf_usm_h%vsws(surf_start)
+                      vsws_int  = surf_usm_h%vsws(surf_start)                  &
+                                * drho_air_zw(k_wall)
                    ENDIF
 !
@@ -3593,5 +3609,5 @@
 !--                as sensitivity studies revealed no significant effect of
 !--                using the neutral solution also for un/stable situations.
-                   v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp )           &
+                   v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp )        &
                             * log_z_z0_int - v_gtrans
 
@@ -3901,4 +3917,5 @@
                 ENDIF
              ENDIF
+
              rvar1_temp(n) = particles(n)%rvar1
              rvar2_temp(n) = particles(n)%rvar2
@@ -3918,5 +3935,4 @@
                 rvar3_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n)          &
                                           + 1E-20_wp ) * ( rg(n,3) - 1.0_wp )
-
              ELSE
 !
@@ -3966,47 +3982,42 @@
 !
 !--    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
+!--    criterion. If yes, limt the SGS particle speed to match the
+!--    CFL criterion. Note, a re-calculation of the SGS particle speed with
+!--    smaller timestep does not necessarily fulfill the CFL criterion as the 
+!--    new SGS speed can be even larger (due to the random term with scales with
+!--    the square-root of dt_particle, for small dt the random contribution increases).
+!--    Thus, we would need to re-calculate the SGS speeds as long as they would
+!--    fulfill the requirements, which could become computationally expensive, 
+!--    Hence, we just limit them.
        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) )
+             IF ( 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
+!
+!--             If total speed exceeds the allowed speed according to CFL 
+!--             criterion, limit the SGS speed to
+!--             dx_i / dt_particle - u_resolved_i, considering a safty factor.
+                rvar1_temp(n) = MERGE( rvar1_temp(n),                          &
+                                       0.9_wp *                                &
+                                       SIGN( dx / dt_particle(n)               &
+                                           - ABS( u_int(n) ), rvar1_temp(n) ), &
+                                       ABS( u_int(n) + rvar1_temp(n) ) <       &
+                                       ( dx / dt_particle(n) ) )
+                rvar2_temp(n) = MERGE( rvar2_temp(n),                          &
+                                       0.9_wp *                                &
+                                       SIGN( dy / dt_particle(n)               &
+                                           - ABS( v_int(n) ), rvar2_temp(n) ), &
+                                       ABS( v_int(n) + rvar2_temp(n) ) <       &
+                                       ( dy / dt_particle(n) ) )
+                rvar3_temp(n) = MERGE( rvar3_temp(n),                          &
+                                       0.9_wp *                                &
+                                       SIGN( zw(kp)-zw(kp-1) / dt_particle(n)  &
+                                           - ABS( w_int(n) ), rvar3_temp(n) ), &
+                                       ABS( w_int(n) + rvar3_temp(n) ) <       &
+                                       ( zw(kp)-zw(kp-1) / dt_particle(n) ) )
              ENDIF
-
 !
 !--          Update particle velocites 
@@ -4276,19 +4287,31 @@
        ENDIF
 !
-!--    Loop from particle start to particle end
+!--    Loop from particle start to particle end and increment the particle 
+!--    age and the total time that the particle has advanced within the 
+!--    particle timestep procedure.
        DO  n = particle_start, particle_end
-!
-!--       Increment the particle age and the total time that the particle
-!--       has advanced within the particle timestep procedure
           particles(n)%age    = particles(n)%age    + dt_particle(n)
           particles(n)%dt_sum = particles(n)%dt_sum + dt_particle(n)
-
-!
-!--       Check whether there is still a particle that has not yet completed
-!--       the total LES timestep
-          IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp )  THEN
+       ENDDO
+!
+!--    Particles that leave the child domain during the SGS-timestep loop
+!--    must not continue timestepping until they are transferred to the 
+!--    parent. Hence, set their dt_sum to dt.
+       IF ( child_domain  .AND.  use_sgs_for_particles )  THEN
+          DO  n = particle_start, particle_end
+             IF ( particles(n)%x < 0.0_wp         .OR.                         &
+                  particles(n)%y < 0.0_wp         .OR.                         &
+                  particles(n)%x > ( nx+1 ) * dx  .OR.                         &
+                  particles(n)%y < ( ny+1 ) * dy )  THEN
+                particles(n)%dt_sum = dt_3d
+             ENDIF
+          ENDDO
+       ENDIF
+!
+!--    Check whether there is still a particle that has not yet completed
+!--    the total LES timestep
+       DO  n = particle_start, particle_end
+          IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp )                     &
              dt_3d_reached_l = .FALSE.
-          ENDIF
-
        ENDDO
     ENDDO
@@ -7060,5 +7083,5 @@
                       IF ( i < 0 )  THEN
 !
-!--                   Apply boundary condition along x
+!--                      Apply boundary condition along x
                          IF ( ibc_par_lr == 0 )  THEN
 !
@@ -7440,4 +7463,5 @@
 
        IF ( trsp_count_recv > 0 )  CALL lpm_add_particles_to_gridcell(rvsp(1:trsp_count_recv))
+
 
        DEALLOCATE(rvsp)
@@ -7742,5 +7766,7 @@
              ENDIF
           ELSE
-             WRITE(0,'(a,8i7)') 'particle out of range ',myid,ip,jp,kp,nxl,nxr,nys,nyn
+             IF ( .NOT. child_domain )  THEN
+                WRITE(0,'(a,8i7)') 'particle out of range ',myid,ip,jp,kp,nxl,nxr,nys,nyn
+             ENDIF
           ENDIF
        ENDIF
@@ -7863,16 +7889,17 @@
 !--             Note, check for CFL does not work at first particle timestep 
 !--             when both, age and age_m are zero. 
-                IF ( particles(n)%age - particles(n)%age_m > 0.0_wp )  THEN  
-                   IF(ABS(particles(n)%speed_x) >                              &
-                      (dx/(particles(n)%age-particles(n)%age_m))  .OR.         &
-                      ABS(particles(n)%speed_y) >                              & 
-                      (dy/(particles(n)%age-particles(n)%age_m))  .OR.         &
-                      ABS(particles(n)%speed_z) >                              &
-                      ((zw(k)-zw(k-1))/(particles(n)%age-particles(n)%age_m))) &
-                   THEN
+                IF ( particles(n)%age - particles(n)%age_m > 0.0_wp )  THEN
+                   IF( ABS( particles(n)%speed_x ) >                           &
+                      ( dx / ( particles(n)%age - particles(n)%age_m) )  .OR.  &
+                       ABS( particles(n)%speed_y ) >                           & 
+                      ( dy / ( particles(n)%age - particles(n)%age_m) )  .OR.  &
+                       ABS( particles(n)%speed_z ) >                           &
+                      ( ( zw(k)-zw(k-1) )                                      &
+                      / ( particles(n)%age - particles(n)%age_m) ) )  THEN
                       WRITE( message_string, * )                               &
                       'Particle violated CFL-criterion: &particle with id ',   &
                       particles(n)%id, ' will be deleted!'   
                       CALL message( 'lpm_check_cfl', 'PA0475', 0, 1, -1, 6, 0 )
+
                       particles(n)%particle_mask= .FALSE.
                    ENDIF
