Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 2627)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 2628)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Enabled particle advection with grid stretching -> Removed parameter check
+! 
+! 2575 2017-10-24 09:57:58Z maronga
 ! Renamed phi --> latitude
 ! 
@@ -880,12 +883,4 @@
 
 !
-!-- Check if vertical grid stretching is used together with particles
-    IF ( dz_stretch_level < 100000.0_wp .AND. particle_advection )  THEN
-       message_string = 'Vertical grid stretching is not allowed together ' // &
-                        'with particle advection.'
-       CALL message( 'check_parameters', 'PA0017', 1, 2, 0, 6, 0 )
-    ENDIF
-
-!
 !-- Check topography setting (check for illegal parameter combinations)
     IF ( topography /= 'flat' )  THEN
Index: /palm/trunk/SOURCE/lpm_advec.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_advec.f90	(revision 2627)
+++ /palm/trunk/SOURCE/lpm_advec.f90	(revision 2628)
@@ -20,5 +20,7 @@
 ! Current revisions:
 ! ------------------
-! 
+!  Removed indices ilog and jlog which are no longer needed since particle box
+! locations are identical to scalar boxes and topography.
+
 ! 
 ! Former revisions:
@@ -167,8 +169,6 @@
     INTEGER(iwp) ::  i                           !< index variable along x
     INTEGER(iwp) ::  ip                          !< index variable along x 
-    INTEGER(iwp) ::  ilog                        !< index variable along x 
     INTEGER(iwp) ::  j                           !< index variable along y 
     INTEGER(iwp) ::  jp                          !< index variable along y 
-    INTEGER(iwp) ::  jlog                        !< index variable along y 
     INTEGER(iwp) ::  k                           !< index variable along z 
     INTEGER(iwp) ::  k_wall                      !< vertical index of topography top 
@@ -299,9 +299,6 @@
 !--       First, check if particle is located below first vertical grid level 
 !--       above topography (Prandtl-layer height)
-          ilog = particles(n)%x * ddx
-          jlog = particles(n)%y * ddy
-!
 !--       Determine vertical index of topography top
-          k_wall = get_topography_top_index( jlog, ilog, 's' )
+          k_wall = get_topography_top_index( jp, ip, 's' )
 
           IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
@@ -327,7 +324,7 @@
 !--             Get friction velocity and momentum flux from new surface data 
 !--             types.
-                IF ( surf_def_h(0)%start_index(jlog,ilog) <=                   &
-                     surf_def_h(0)%end_index(jlog,ilog) )  THEN
-                   surf_start = surf_def_h(0)%start_index(jlog,ilog)
+                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)
 !--                Limit friction velocity. In narrow canyons or holes the 
 !--                friction velocity can become very small, resulting in a too
@@ -335,12 +332,12 @@
                    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(jlog,ilog) <=                  &
-                         surf_lsm_h%end_index(jlog,ilog) )  THEN
-                   surf_start = surf_lsm_h%start_index(jlog,ilog)
+                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(jlog,ilog) <=                  &
-                         surf_usm_h%end_index(jlog,ilog) )  THEN
-                   surf_start = surf_usm_h%start_index(jlog,ilog)
+                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)
@@ -396,9 +393,7 @@
        DO  n = start_index(nb), end_index(nb)
 
-          ilog = particles(n)%x * ddx
-          jlog = particles(n)%y * ddy
 !
 !--       Determine vertical index of topography top
-          k_wall = get_topography_top_index( jlog,ilog, 's' )
+          k_wall = get_topography_top_index( jp,ip, 's' )
 
           IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
@@ -427,7 +422,7 @@
 !--             Get friction velocity and momentum flux from new surface data 
 !--             types.
-                IF ( surf_def_h(0)%start_index(jlog,ilog) <=                   &
-                     surf_def_h(0)%end_index(jlog,ilog) )  THEN
-                   surf_start = surf_def_h(0)%start_index(jlog,ilog)
+                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)
 !--                Limit friction velocity. In narrow canyons or holes the 
 !--                friction velocity can become very small, resulting in a too
@@ -435,12 +430,12 @@
                    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(jlog,ilog) <=                  &
-                         surf_lsm_h%end_index(jlog,ilog) )  THEN
-                   surf_start = surf_lsm_h%start_index(jlog,ilog)
+                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(jlog,ilog) <=                  &
-                         surf_usm_h%end_index(jlog,ilog) )  THEN
-                   surf_start = surf_usm_h%start_index(jlog,ilog)
+                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)
Index: /palm/trunk/SOURCE/lpm_exchange_horiz.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_exchange_horiz.f90	(revision 2627)
+++ /palm/trunk/SOURCE/lpm_exchange_horiz.f90	(revision 2628)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Enabled particle advection with grid stretching. Furthermore, the CFL-
+! criterion is checked for every particle at every time step.
+! 
+! 2606 2017-11-10 10:36:31Z schwenkel
 ! Changed particle box locations: center of particle box now coincides 
 ! with scalar grid point of same index.
@@ -124,5 +128,8 @@
  
     USE, INTRINSIC ::  ISO_C_BINDING
-
+    
+    USE arrays_3d,                                                             &
+       ONLY:  zw
+    
     USE control_parameters,                                                    &
         ONLY:  dz, message_string, simulated_time
@@ -887,4 +894,14 @@
        jp = particle_array(n)%y * ddy
        kp = particle_array(n)%z / dz + 1 + offset_ocean_nzt
+!
+!--    In case of grid stretching a particle might be above or below the  
+!--    previously calculated particle grid box (indices).     
+       DO WHILE( zw(kp) < particle_array(n)%z )
+          kp = kp + 1
+       ENDDO
+
+       DO WHILE( zw(kp-1) > particle_array(n)%z )
+          kp = kp - 1
+       ENDDO
 
        IF ( ip >= nxl  .AND.  ip <= nxr  .AND.  jp >= nys  .AND.  jp <= nyn    &
@@ -1022,5 +1039,5 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE lpm_move_particle
-
+ 
     IMPLICIT NONE
 
@@ -1038,5 +1055,5 @@
 
     CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'start' )
-
+    CALL lpm_check_cfl
     DO  ip = nxl, nxr
        DO  jp = nys, nyn
@@ -1050,6 +1067,11 @@
                 i = particles_before_move(n)%x * ddx
                 j = particles_before_move(n)%y * ddy
-                k = particles_before_move(n)%z / dz + 1 + offset_ocean_nzt
-
+                k = kp
+!
+!--             Find correct vertical particle grid box (necessary in case of grid stretching)
+!--             Due to the CFL limitations only the neighbouring grid boxes are considered. 
+                IF( zw(k)   < particles_before_move(n)%z ) k = k + 1
+                IF( zw(k-1) > particles_before_move(n)%z ) k = k - 1 
+                
 !--             For lpm_exchange_horiz to work properly particles need to be moved to the outermost gridboxes
 !--             of the respective processor. If the particle index is inside the processor the following lines
@@ -1059,6 +1081,8 @@
                 j = MIN ( j , nyn )
                 j = MAX ( j , nys )
+                
                 k = MIN ( k , nzt )
                 k = MAX ( k , nzb+1 )
+                
 !
 !--             Check, if particle has moved to another grid cell.
@@ -1092,5 +1116,46 @@
 
  END SUBROUTINE lpm_move_particle
-
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check CFL-criterion for each particle. If one particle violated the 
+!> criterion the particle will be deleted and a warning message is given.
+!------------------------------------------------------------------------------!
+ SUBROUTINE lpm_check_cfl  
+       
+    IMPLICIT NONE
+    
+    INTEGER(iwp)  ::  i
+    INTEGER(iwp)  ::  j
+    INTEGER(iwp)  ::  k
+    INTEGER(iwp)  ::  n
+    
+    DO  i = nxl, nxr
+       DO  j = nys, nyn
+          DO  k = nzb+1, nzt
+             number_of_particles = prt_count(k,j,i)
+             IF ( number_of_particles <= 0 )  CYCLE
+             particles => grid_particles(k,j,i)%particles(1:number_of_particles)         
+             DO n = 1, number_of_particles
+   
+                IF(ABS(particles(n)%speed_x) >                                 &
+                   (dx/(particles(n)%age-particles(n)%age_m))  .OR.            &
+                   ABS(particles(n)%speed_y) >                                 & 
+                   (dx/(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', 'PA0500', 0, 1, 0, 6, 0 )
+                   particles(n)%particle_mask= .FALSE.
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDDO   
+
+ END SUBROUTINE lpm_check_cfl
+ 
 !------------------------------------------------------------------------------!
 ! Description:
Index: /palm/trunk/SOURCE/lpm_init.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_init.f90	(revision 2627)
+++ /palm/trunk/SOURCE/lpm_init.f90	(revision 2628)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Enabled particle advection with grid stretching.
+! 
+! 2608 2017-11-13 14:04:26Z schwenkel
 ! Calculation of magnus equation in external module (diagnostic_quantities_mod).
 ! 
@@ -624,9 +627,11 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE lpm_create_particle (phase)
-
+    
+    USE arrays_3d,                                                             &
+       ONLY:  zw
     USE lpm_exchange_horiz_mod,                                                &
         ONLY: lpm_exchange_horiz, lpm_move_particle, realloc_particles_array
 
-    USE lpm_pack_and_sort_mod,                                                   &
+    USE lpm_pack_and_sort_mod,                                                 &
         ONLY: lpm_sort_in_subboxes
 
@@ -742,5 +747,11 @@
                                ip = tmp_particle%x * ddx
                                jp = tmp_particle%y * ddy
-                               kp = tmp_particle%z / dz + 1 + offset_ocean_nzt
+                               kp = tmp_particle%z / dz + 1 + offset_ocean_nzt                               
+                               DO WHILE( zw(kp) < tmp_particle%z ) 
+                                  kp = kp + 1
+                               ENDDO
+                               DO WHILE( zw(kp-1) > tmp_particle%z )
+                                  kp = kp - 1
+                               ENDDO 
 !
 !--                            Determine surface level. Therefore, check for
@@ -929,5 +940,11 @@
                    i = particles(n)%x * ddx
                    j = particles(n)%y * ddy
-                   k =   particles(n)%z / dz + 1 + offset_ocean_nzt
+                   k = particles(n)%z / dz + 1 + offset_ocean_nzt
+                   DO WHILE( zw(k) < particles(n)%z )
+                      k = k + 1
+                   ENDDO
+                   DO WHILE( zw(k-1) > particles(n)%z )
+                      k = k - 1
+                   ENDDO
 !
 !--                Check if particle is within topography
Index: /palm/trunk/SOURCE/lpm_pack_arrays.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_pack_arrays.f90	(revision 2627)
+++ /palm/trunk/SOURCE/lpm_pack_arrays.f90	(revision 2628)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Enabled particle advection with grid stretching.
+! 
+! 2609 2017-11-14 14:14:44Z schwenkel
 ! Integrated subroutine pack_and_sort into lpm_sort_in_subboxes
 ! 
@@ -112,17 +115,19 @@
 
        USE cpulog,                                                              &
-           ONLY:  cpu_log, log_point_s
+          ONLY:  cpu_log, log_point_s
 
        USE indices,                                                             &
-           ONLY:  nxl, nxr, nys, nyn, nzb, nzt
+          ONLY:  nxl, nxr, nys, nyn, nzb, nzt
 
        USE kinds
 
        USE control_parameters,                                                  &
-           ONLY: dz
+          ONLY: dz
        
        USE grid_variables,                                                      &
-           ONLY: dx,dy,ddx, ddy
-          
+          ONLY: dx,dy,ddx, ddy
+           
+       USE arrays_3d,                                                           &
+          ONLY:  zu
        IMPLICIT NONE
 
@@ -167,9 +172,8 @@
                       i = (particles(n)%x + 0.5_wp * dx) * ddx
                       j = (particles(n)%y + 0.5_wp * dy) * ddy
-                      k = ( particles(n)%z+ 0.5_wp *dz ) / dz+1 + offset_ocean_nzt
                       
                       IF ( i == ip )  sort_index = sort_index + 4 
-                      IF ( j == jp )  sort_index = sort_index + 2
-                      IF ( k == kp )  sort_index = sort_index + 1
+                      IF ( j == jp )  sort_index = sort_index + 2                      
+                      IF ( zu(kp) > particles(n)%z ) sort_index = sort_index + 1
                       
                       sort_count(sort_index) = sort_count(sort_index) + 1
@@ -343,3 +347,3 @@
 
 
- END module lpm_pack_and_sort_mod
+ END MODULE lpm_pack_and_sort_mod
