Index: palm/trunk/SOURCE/lpm_init.f90
===================================================================
--- palm/trunk/SOURCE/lpm_init.f90	(revision 2221)
+++ palm/trunk/SOURCE/lpm_init.f90	(revision 2223)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Add check for particle release at model top
 ! 
 ! Former revisions:
@@ -618,104 +618,109 @@
           DO WHILE ( pos_z <= pst(i) )
 
-             pos_y = pss(i)
-
-             DO WHILE ( pos_y <= psn(i) )
-
-                IF ( pos_y >= ( nys - 0.5_wp ) * dy  .AND.  &
-                     pos_y <  ( nyn + 0.5_wp ) * dy )  THEN
-
-                   pos_x = psl(i)
-
-            xloop: DO WHILE ( pos_x <= psr(i) )
-
-                      IF ( pos_x >= ( nxl - 0.5_wp ) * dx  .AND.  &
-                           pos_x <  ( nxr + 0.5_wp ) * dx )  THEN
-
-                         DO  j = 1, particles_per_point
-
-                            n = n + 1
-                            tmp_particle%x             = pos_x
-                            tmp_particle%y             = pos_y
-                            tmp_particle%z             = pos_z
-                            tmp_particle%age           = 0.0_wp
-                            tmp_particle%age_m         = 0.0_wp
-                            tmp_particle%dt_sum        = 0.0_wp
-                            tmp_particle%user          = 0.0_wp !unused, free for the user
-                            tmp_particle%e_m           = 0.0_wp
-                            IF ( curvature_solution_effects )  THEN
-!
-!--                            Initial values (internal timesteps, derivative)
-!--                            for Rosenbrock method
-                               tmp_particle%rvar1      = 1.0E-6_wp     !last Rosenbrock timestep
-                               tmp_particle%rvar2      = 0.1E-6_wp     !dry aerosol radius
-                               tmp_particle%rvar3      = -9999999.9_wp !unused in this configuration
-                            ELSE
-!
-!--                            Initial values for SGS velocities
-                               tmp_particle%rvar1      = 0.0_wp
-                               tmp_particle%rvar2      = 0.0_wp
-                               tmp_particle%rvar3      = 0.0_wp
-                            ENDIF
-                            tmp_particle%speed_x       = 0.0_wp
-                            tmp_particle%speed_y       = 0.0_wp
-                            tmp_particle%speed_z       = 0.0_wp
-                            tmp_particle%origin_x      = pos_x
-                            tmp_particle%origin_y      = pos_y
-                            tmp_particle%origin_z      = pos_z
-                            tmp_particle%radius        = particle_groups(i)%radius
-                            tmp_particle%weight_factor = initial_weighting_factor
-                            tmp_particle%class         = 1
-                            tmp_particle%group         = i
-                            tmp_particle%id1           = 0 
-                            tmp_particle%id2           = 0 
-                            tmp_particle%particle_mask = .TRUE.
-                            tmp_particle%block_nr      = -1
-!
-!--                         Determine the grid indices of the particle position
-                            ip = ( tmp_particle%x + 0.5_wp * dx ) * ddx
-                            jp = ( tmp_particle%y + 0.5_wp * dy ) * ddy
-                            kp = tmp_particle%z / dz + 1 + offset_ocean_nzt
-
-                            IF ( seed_follows_topography )  THEN
-!
-!--                            Particle height is given relative to topography
-                               kp = kp + nzb_w_inner(jp,ip)
-                               tmp_particle%z = tmp_particle%z +               &
+             IF ( pos_z >= 0.0_wp  .AND.  pos_z < zw(nzt) )  THEN
+
+
+                pos_y = pss(i)
+
+                DO WHILE ( pos_y <= psn(i) )
+
+                   IF ( pos_y >= ( nys - 0.5_wp ) * dy  .AND.                  &
+                        pos_y <  ( nyn + 0.5_wp ) * dy )  THEN
+
+                      pos_x = psl(i)
+
+               xloop: DO WHILE ( pos_x <= psr(i) )
+
+                         IF ( pos_x >= ( nxl - 0.5_wp ) * dx  .AND.            & 
+                              pos_x <  ( nxr + 0.5_wp ) * dx )  THEN
+
+                            DO  j = 1, particles_per_point
+
+                               n = n + 1
+                               tmp_particle%x             = pos_x
+                               tmp_particle%y             = pos_y
+                               tmp_particle%z             = pos_z
+                               tmp_particle%age           = 0.0_wp
+                               tmp_particle%age_m         = 0.0_wp
+                               tmp_particle%dt_sum        = 0.0_wp
+                               tmp_particle%user          = 0.0_wp !unused, free for the user
+                               tmp_particle%e_m           = 0.0_wp
+                               IF ( curvature_solution_effects )  THEN
+!
+!--                               Initial values (internal timesteps, derivative)
+!--                               for Rosenbrock method
+                                  tmp_particle%rvar1      = 1.0E-6_wp     !last Rosenbrock timestep
+                                  tmp_particle%rvar2      = 0.1E-6_wp     !dry aerosol radius
+                                  tmp_particle%rvar3      = -9999999.9_wp !unused in this configuration
+                               ELSE
+!
+!--                               Initial values for SGS velocities
+                                  tmp_particle%rvar1      = 0.0_wp
+                                  tmp_particle%rvar2      = 0.0_wp
+                                  tmp_particle%rvar3      = 0.0_wp
+                               ENDIF
+                               tmp_particle%speed_x       = 0.0_wp
+                               tmp_particle%speed_y       = 0.0_wp
+                               tmp_particle%speed_z       = 0.0_wp
+                               tmp_particle%origin_x      = pos_x
+                               tmp_particle%origin_y      = pos_y
+                               tmp_particle%origin_z      = pos_z
+                               tmp_particle%radius        = particle_groups(i)%radius
+                               tmp_particle%weight_factor = initial_weighting_factor
+                               tmp_particle%class         = 1
+                               tmp_particle%group         = i
+                               tmp_particle%id1           = 0 
+                               tmp_particle%id2           = 0 
+                               tmp_particle%particle_mask = .TRUE.
+                               tmp_particle%block_nr      = -1
+!
+!--                            Determine the grid indices of the particle position
+                              ip = ( tmp_particle%x + 0.5_wp * dx ) * ddx
+                               jp = ( tmp_particle%y + 0.5_wp * dy ) * ddy
+                               kp = tmp_particle%z / dz + 1 + offset_ocean_nzt
+
+                               IF ( seed_follows_topography )  THEN
+!
+!--                               Particle height is given relative to topography
+                                  kp = kp + nzb_w_inner(jp,ip)
+                                  tmp_particle%z = tmp_particle%z +            &
                                                          zw(nzb_w_inner(jp,ip))
-                               IF ( kp > nzt )  THEN
+                                  IF ( kp > nzt )  THEN
+                                     pos_x = pos_x + pdx(i)
+                                     CYCLE xloop
+                                  ENDIF
+                               ELSEIF ( .NOT. seed_follows_topography .AND.    &
+                                        tmp_particle%z <= zw(nzb_w_inner(jp,ip)) )  THEN
                                   pos_x = pos_x + pdx(i)
-                                  CYCLE xloop
+                                  CYCLE xloop                               
                                ENDIF
-                            ELSEIF ( .NOT. seed_follows_topography .AND.       &
-                                     tmp_particle%z <= zw(nzb_w_inner(jp,ip)) )  THEN
-                               pos_x = pos_x + pdx(i)
-                               CYCLE xloop                               
-                            ENDIF
-
-                            local_count(kp,jp,ip) = local_count(kp,jp,ip) + 1
-
-                            IF ( .NOT. first_stride )  THEN
-                               IF ( ip < nxl  .OR.  jp < nys  .OR.  kp < nzb+1 )  THEN
-                                  write(6,*) 'xl ',ip,jp,kp,nxl,nys,nzb+1
+
+                               local_count(kp,jp,ip) = local_count(kp,jp,ip) + 1
+
+                               IF ( .NOT. first_stride )  THEN
+                                  IF ( ip < nxl  .OR.  jp < nys  .OR.  kp < nzb+1 )  THEN
+                                     write(6,*) 'xl ',ip,jp,kp,nxl,nys,nzb+1
+                                  ENDIF
+                                  IF ( ip > nxr  .OR.  jp > nyn  .OR.  kp > nzt )  THEN
+                                     write(6,*) 'xu ',ip,jp,kp,nxr,nyn,nzt
+                                  ENDIF
+                                  grid_particles(kp,jp,ip)%particles(local_count(kp,jp,ip)) = tmp_particle
+
                                ENDIF
-                               IF ( ip > nxr  .OR.  jp > nyn  .OR.  kp > nzt )  THEN
-                                  write(6,*) 'xu ',ip,jp,kp,nxr,nyn,nzt
-                               ENDIF
-                               grid_particles(kp,jp,ip)%particles(local_count(kp,jp,ip)) = tmp_particle
-
-                            ENDIF
-                         ENDDO
-
-                      ENDIF
-
-                      pos_x = pos_x + pdx(i)
-
-                   ENDDO xloop
-
-                ENDIF
-
-                pos_y = pos_y + pdy(i)
-
-             ENDDO
+                            ENDDO
+
+                         ENDIF
+
+                         pos_x = pos_x + pdx(i)
+
+                      ENDDO xloop
+
+                   ENDIF
+
+                   pos_y = pos_y + pdy(i)
+
+                ENDDO
+
+             ENDIF
 
              pos_z = pos_z + pdz(i)
