Index: palm/trunk/DOC/app/chapter_4.1.html
===================================================================
--- palm/trunk/DOC/app/chapter_4.1.html	(revision 56)
+++ palm/trunk/DOC/app/chapter_4.1.html	(revision 57)
@@ -11404,5 +11404,5 @@
 
 
-    <tr>
+    <tr><td style="vertical-align: top;"><a name="pt_reference"></a><span style="font-weight: bold;">pt_reference</span></td><td style="vertical-align: top;">R</td><td style="vertical-align: top;"><span style="font-style: italic;">use horizontal average as refrence</span></td><td style="vertical-align: top;">Reference temperature to be used in all buoyancy terms (in K).<br><br>By default, the instantaneous horizontal average over the total model domain is used.</td></tr><tr>
 
 
Index: palm/trunk/DOC/app/chapter_4.6.html
===================================================================
--- palm/trunk/DOC/app/chapter_4.6.html	(revision 56)
+++ palm/trunk/DOC/app/chapter_4.6.html	(revision 57)
@@ -16681,5 +16681,5 @@
 
 
-    <tr>
+    <tr><td><a style="font-weight: bold;" href="chapter_4.1.html#pt_reference">pt_reference</a></td><td>I</td><td>R</td><td><span style="font-style: italic;">use horizontal average as reference</span></td><td>Reference temperature to be used in all buoyancy terms (in K).</td></tr><tr>
 
 
Index: palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 56)
+++ palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 57)
@@ -1,8 +1,14 @@
 New:
 ---
+
+particle reflection from vertical walls implemented, particle SGS model adjusted to walls
 
 Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.
 
 new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files
+
+new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).
+
+new user interface user_advec_particles
 
 new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays
@@ -12,5 +18,7 @@
 samples added to the user interface which show how to add user-define time series quantities.
 
-Makefile, check_open, check_parameters, header, init_3d_model, modules, netcdf, parin, user_interface
+unit 9 opened for debug output (file DEBUG_<pe#>)
+
+Makefile, advec_particles, buoyancy, check_open, check_parameters, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, user_interface, write_var_list
 
 New: wall_fluxes
@@ -19,4 +27,5 @@
 Changed:
 -------
++age_m in particle_type
 
 Move call of user_actions( 'after_integration' ) below increment of times
@@ -27,5 +36,7 @@
 Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.
 
-check_parameters, data_output_ts, flow_statistics, init_3d_model, modules, parin, time_integration
+__vtk directives removed from main program.
+
+check_parameters, data_output_ts, flow_statistics, init_particles, init_3d_model, modules, palm, parin, time_integration
 
 
@@ -35,7 +46,9 @@
 Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model
 
+Bugfix in sample for reading user defined data from restart file (user_init)
+
 in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
 # case of .mod files
 
 Makefile
-init_1d_model
+init_1d_model, user_interface
Index: palm/trunk/SOURCE/advec_particles.f90
===================================================================
--- palm/trunk/SOURCE/advec_particles.f90	(revision 56)
+++ palm/trunk/SOURCE/advec_particles.f90	(revision 57)
@@ -4,4 +4,7 @@
 ! Actual revisions:
 ! -----------------
+! Particle reflection at vertical walls (by Jin Zhang), regarding vertical
+! walls in the SGS model,
+! + user interface user_advec_particles
 ! TEST: PRINT statements on unit 9 (commented out)
 !
@@ -37,116 +40,4 @@
 ! sub-timesteps (if sgs velocities are included), +3d-arrays de_dx/y/z,
 ! izuf renamed iran, output of particle time series
-!
-! Revision 1.28  2006/02/23 09:41:21  raasch
-! Only a fraction of the particles may have tails, therefore output for unit
-! 85 and 90 has changed,
-! sorting of particles in gridbox order after each timestep, optimization of
-! the droplet growth by collision,
-! 0.5 replaced by 0.4999999999 in setting the random fluctuations,
-! psl, psr, pdx, etc. are now 1d arrays, i.e. they depend on the particle group,
-! improve particle release at PE boundaries, nt_anz renamed
-! current_timestep_number, idum (particle_type) renamed tail_id,
-! error number argument for handle_netcdf_error,
-!
-! Revision 1.27  2005/10/20 14:04:15  raasch
-! droplet growth by collision completed (still has to be optimized for speed),
-! new routine collision_efficiency (see end of this file),
-! 2*r replaced by r in the exponential terms (momentum equation for particles),
-! number of particles really used is additionally output on the netcdf particle
-! data file,
-! test print statement removed
-!
-! Revision 1.26  2005/06/26 19:50:14  raasch
-! Cloud droplet features added, radius is used instead of diameter,
-! particle type initial data assignment extended,
-! output on unit 80 for non-parallel runs
-!
-! Revision 1.25  2005/05/18 14:53:11  raasch
-! Extensions for NetCDF output
-!
-! Revision 1.24  2005/04/23 08:24:58  raasch
-! Revised calculation of output time counters regarding a possible decrease of
-! the output time interval in case of restart runs
-!
-! Revision 1.23  2005/03/26 14:20:05  raasch
-! calculation of vertical particle velocity (with inertia) corrected, exp_arg
-! had a wrong sign
-!
-! Revision 1.22  2004/04/30 07:48:57  raasch
-! Particle type initial data assignment extended
-!
-! Revision 1.21  2004/01/28 15:00:48  raasch
-! Output of particle infos in subroutine allocate_prt_memory on demand only
-!
-! Revision 1.20  2003/10/29 08:38:52  raasch
-! Module random_function_mod is used, modifications for new particle group
-! feature, version number is written on binary file
-!
-! Revision 1.19  2003/03/16 09:19:32  raasch
-! Two underscores (_) are placed in front of all define-strings
-!
-! Revision 1.18  2003/03/14 13:36:31  raasch
-! Error in reflection boundary condition removed: particle speed must also
-! be inverted
-!
-! Revision 1.17  2003/03/04 11:23:20  raasch
-! Error in particle inertia part removed (exp_arg must not contain the timestep)
-! Particle velocities are also stored in particles in case of zero density ratio
-!
-! Revision 1.16  2002/09/12 12:55:29  raasch
-! Transport of particles with inertia implemented, write density_ratio to
-! restart file
-!
-! Revision 1.15  2002/06/11 12:11:56  raasch
-! Declaration of u_int_l, u_int_u, v_int_l, v_int_u added
-!
-! Revision 1.14  2002/05/02  18:46:34  18:46:34  raasch (Siegfried Raasch)
-! Horizontal velocity components for particle advection are now interpolated
-! between horizontal grid levels
-! 
-! Revision 1.13  2002/04/24  18:46:11  18:46:11  raasch (Siegfried Raasch)
-! Temporary particle-arrays trlp, trnp, trrp and trsp are now allocatable arrays
-! 
-! Revision 1.12  2002/04/16 07:44:07  raasch
-! Additional boundary conditions added, particle data for later analysis can be
-! written on file
-!
-! Revision 1.11  2001/11/09 13:38:20  raasch
-! Colour information for particle tails is stored in array
-! particle_tail_coordinates, adding information to the particle tails
-! moved after call of user_particle_attributes
-!
-! Revision 1.10  2001/08/21 08:19:48  raasch
-! Storage of particle tails, precision of vertical advection improved,
-! particle attributes (e.g. size and color) can be defined by the user
-!
-! Revision 1.9  2001/07/12 12:01:27  raasch
-! Random start positions of initial particles possible
-!
-! Revision 1.8  2001/03/29 17:27:10  raasch
-! Translation of remaining German identifiers (variables, subroutines, etc.)
-!
-! Revision 1.7  2001/02/13 10:24:00  raasch
-! Particle advection possible in case of galilei transformation
-!
-! Revision 1.6  2001/01/25 06:51:04  raasch
-! Module test_variables removed, writing of particle informations is optional
-! now
-!
-! Revision 1.5  2001/01/02 17:15:29  raasch
-! Unit 80 is immediately closed after usage every time. Particle positions
-! are written on unit 90 instead of 81 (subroutine write_particles).
-!
-! Revision 1.4  2000/12/28 12:43:54  raasch
-! Packing of particles-array is done explicitly (not by PACK-function, which
-! seems to have problems when a large number of PEs is used)
-! Missing translations included,
-! code is used only optionally (cpp-directives are added)
-!
-! Revision 1.3  2000/01/20 08:51:37  letzel
-! All comments translated into English
-!
-! Revision 1.2  1999/11/25 16:41:48  raasch
-! Indexfehler im nichtparallelen Teil korrigiert
 !
 ! Revision 1.1  1999/11/25 16:16:06  raasch
@@ -693,6 +584,47 @@
           DO  j = nys, nyn
              DO  k = nzb, nzt+1
-                de_dx(k,j,i) = sgs_wfu_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx
-                de_dy(k,j,i) = sgs_wfv_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy
+
+                IF ( k <= nzb_s_inner(j,i-1)  .AND.  &
+                     k  > nzb_s_inner(j,i)    .AND.  &
+                     k  > nzb_s_inner(j,i+1) )  THEN
+                   de_dx(k,j,i) = 2.0 * sgs_wfu_part * &
+                                  ( e(k,j,i+1) - e(k,j,i) ) * ddx
+                ELSEIF ( k  > nzb_s_inner(j,i-1)  .AND.  &
+                         k  > nzb_s_inner(j,i)    .AND.  &
+                         k <= nzb_s_inner(j,i+1) )  THEN
+                   de_dx(k,j,i) = 2.0 * sgs_wfu_part * &
+                                  ( e(k,j,i) - e(k,j,i-1) ) * ddx
+                ELSEIF ( k < nzb_s_inner(j,i)  .AND.  k < nzb_s_inner(j,i+1) ) &
+                THEN
+                   de_dx(k,j,i) = 0.0
+                ELSEIF ( k < nzb_s_inner(j,i-1)  .AND.  k < nzb_s_inner(j,i) ) &
+                THEN
+                   de_dx(k,j,i) = 0.0
+                ELSE
+                   de_dx(k,j,i) = sgs_wfu_part * &
+                                  ( e(k,j,i+1) - e(k,j,i-1) ) * ddx
+                ENDIF
+
+                IF ( k <= nzb_s_inner(j-1,i)  .AND.  &
+                     k  > nzb_s_inner(j,i)    .AND.  &
+                     k  > nzb_s_inner(j+1,i) )  THEN
+                   de_dy(k,j,i) = 2.0 * sgs_wfv_part * &
+                                  ( e(k,j+1,i) - e(k,j,i) ) * ddy
+                ELSEIF ( k  > nzb_s_inner(j-1,i)  .AND.  &
+                         k  > nzb_s_inner(j,i)    .AND.  &
+                         k <= nzb_s_inner(j+1,i) )  THEN
+                   de_dy(k,j,i) = 2.0 * sgs_wfv_part * &
+                                  ( e(k,j,i) - e(k,j-1,i) ) * ddy
+                ELSEIF ( k < nzb_s_inner(j,i)  .AND.  k < nzb_s_inner(j+1,i) ) &
+                THEN
+                   de_dy(k,j,i) = 0.0
+                ELSEIF ( k < nzb_s_inner(j-1,i)  .AND.  k < nzb_s_inner(j,i) ) &
+                THEN
+                   de_dy(k,j,i) = 0.0
+                ELSE
+                   de_dy(k,j,i) = sgs_wfv_part * &
+                                  ( e(k,j+1,i) - e(k,j-1,i) ) * ddy
+                ENDIF
+
              ENDDO
           ENDDO
@@ -703,16 +635,18 @@
        DO  i = nxl, nxr
           DO  j = nys, nyn
-             DO  k = nzb+2, nzt-1
+
+             DO  k = nzb_s_inner(j,i)+2, nzt-1
                 de_dz(k,j,i) = 2.0 * sgs_wfw_part * &
                                ( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1)-zu(k-1) )
              ENDDO
-             de_dz(nzb,j,i)   = 0.0
-             de_dz(nzb+1,j,i) = 2.0 * sgs_wfw_part *              &
-                                ( e(nzb+2,j,i) - e(nzb+1,j,i) ) / &
-                                ( zu(nzb+2) - zu(nzb+1) )
+
+             k = nzb_s_inner(j,i)
+             de_dz(nzb:k,j,i)   = 0.0
+             de_dz(k+1,j,i) = 2.0 * sgs_wfw_part * ( e(k+2,j,i) - e(k+1,j,i) ) &
+                                                 / ( zu(k+2) - zu(k+1) )
              de_dz(nzt,j,i)   = 0.0
              de_dz(nzt+1,j,i) = 0.0
           ENDDO
-       ENDDO              
+       ENDDO
 
 !
@@ -967,107 +901,706 @@
              k = ( particles(n)%z + 0.5 * dz ) / dz  ! only exact if eq.dist
 
-             x  = particles(n)%x - i * dx
-             y  = particles(n)%y - j * dy
-             aa = x**2          + y**2
-             bb = ( dx - x )**2 + y**2
-             cc = x**2          + ( dy - y )**2
-             dd = ( dx - x )**2 + ( dy - y )**2
-             gg = aa + bb + cc + dd
-
-             e_int_l = ( ( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1)   &
-                       + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
-                       ) / ( 3.0 * gg )
-
-             IF ( k+1 == nzt+1 )  THEN
-                e_int = e_int_l
+             IF ( topography == 'flat' )  THEN        
+
+                x  = particles(n)%x - i * dx
+                y  = particles(n)%y - j * dy
+                aa = x**2          + y**2
+                bb = ( dx - x )**2 + y**2
+                cc = x**2          + ( dy - y )**2
+                dd = ( dx - x )**2 + ( dy - y )**2
+                gg = aa + bb + cc + dd
+
+                e_int_l = ( ( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1)   &
+                          + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
+                          ) / ( 3.0 * gg )
+
+                IF ( k+1 == nzt+1 )  THEN
+                   e_int = e_int_l
+                ELSE
+                   e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
+                               ( gg - bb ) * e(k+1,j,i+1) + &
+                               ( gg - cc ) * e(k+1,j+1,i) + &
+                               ( gg - dd ) * e(k+1,j+1,i+1) &
+                             ) / ( 3.0 * gg )
+                   e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
+                                     ( e_int_u - e_int_l )
+                ENDIF
+
+!
+!--             Interpolate the TKE gradient along x (adopt incides i,j,k and
+!--             all position variables from above (TKE))
+                de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
+                                ( gg - bb ) * de_dx(k,j,i+1) + &
+                                ( gg - cc ) * de_dx(k,j+1,i) + &
+                                ( gg - dd ) * de_dx(k,j+1,i+1) &
+                              ) / ( 3.0 * gg )
+
+                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
+                   de_dx_int = de_dx_int_l
+                ELSE
+                   de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
+                                   ( gg - bb ) * de_dx(k+1,j,i+1) + &
+                                   ( gg - cc ) * de_dx(k+1,j+1,i) + &
+                                   ( gg - dd ) * de_dx(k+1,j+1,i+1) &
+                                 ) / ( 3.0 * gg )
+                   de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * &
+                                             ( de_dx_int_u - de_dx_int_l )
+                ENDIF
+
+!
+!--             Interpolate the TKE gradient along y
+                de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
+                                ( gg - bb ) * de_dy(k,j,i+1) + &
+                                ( gg - cc ) * de_dy(k,j+1,i) + &
+                                ( gg - dd ) * de_dy(k,j+1,i+1) &
+                              ) / ( 3.0 * gg )
+                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
+                   de_dy_int = de_dy_int_l
+                ELSE
+                   de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
+                                   ( gg - bb ) * de_dy(k+1,j,i+1) + &
+                                   ( gg - cc ) * de_dy(k+1,j+1,i) + &
+                                   ( gg - dd ) * de_dy(k+1,j+1,i+1) &
+                                 ) / ( 3.0 * gg )
+                   de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * &
+                                             ( de_dy_int_u - de_dy_int_l )
+                ENDIF
+
+!
+!--             Interpolate the TKE gradient along z
+                IF ( particles(n)%z < 0.5 * dz ) THEN
+                   de_dz_int = 0.0
+                ELSE
+                   de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
+                                   ( gg - bb ) * de_dz(k,j,i+1) + &
+                                   ( gg - cc ) * de_dz(k,j+1,i) + &
+                                   ( gg - dd ) * de_dz(k,j+1,i+1) &
+                                 ) / ( 3.0 * gg )
+
+                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
+                      de_dz_int = de_dz_int_l
+                   ELSE
+                      de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
+                                      ( gg - bb ) * de_dz(k+1,j,i+1) + &
+                                      ( gg - cc ) * de_dz(k+1,j+1,i) + &
+                                      ( gg - dd ) * de_dz(k+1,j+1,i+1) &
+                                    ) / ( 3.0 * gg )
+                      de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * &
+                                                ( de_dz_int_u - de_dz_int_l )
+                   ENDIF
+                ENDIF
+
+!
+!--             Interpolate the dissipation of TKE
+                diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
+                               ( gg - bb ) * diss(k,j,i+1) + &
+                               ( gg - cc ) * diss(k,j+1,i) + &
+                               ( gg - dd ) * diss(k,j+1,i+1) &
+                              ) / ( 3.0 * gg )
+
+                IF ( k+1 == nzt+1 )  THEN
+                   diss_int = diss_int_l
+                ELSE
+                   diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
+                                  ( gg - bb ) * diss(k+1,j,i+1) + &
+                                  ( gg - cc ) * diss(k+1,j+1,i) + &
+                                  ( gg - dd ) * diss(k+1,j+1,i+1) &
+                                 ) / ( 3.0 * gg )
+                   diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * &
+                                           ( diss_int_u - diss_int_l )
+                ENDIF
+
              ELSE
-                e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
-                            ( gg - bb ) * e(k+1,j,i+1) + &
-                            ( gg - cc ) * e(k+1,j+1,i) + &
-                            ( gg - dd ) * e(k+1,j+1,i+1) &
-                          ) / ( 3.0 * gg )
-                e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
-                                  ( e_int_u - e_int_l )
-             ENDIF
-
-!
-!--          Interpolate the TKE gradient along x (adopt incides i,j,k and all
-!--          position variables from above (TKE))
-             de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
-                             ( gg - bb ) * de_dx(k,j,i+1) + &
-                             ( gg - cc ) * de_dx(k,j+1,i) + &
-                             ( gg - dd ) * de_dx(k,j+1,i+1) &
-                           ) / ( 3.0 * gg )
-
-             IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
-                de_dx_int = de_dx_int_l
-             ELSE
-                de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
-                                ( gg - bb ) * de_dx(k+1,j,i+1) + &
-                                ( gg - cc ) * de_dx(k+1,j+1,i) + &
-                                ( gg - dd ) * de_dx(k+1,j+1,i+1) &
-                              ) / ( 3.0 * gg )
-                de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * &
-                                          ( de_dx_int_u - de_dx_int_l )
-             ENDIF
-
-!
-!--          Interpolate the TKE gradient along y
-             de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
-                             ( gg - bb ) * de_dy(k,j,i+1) + &
-                             ( gg - cc ) * de_dy(k,j+1,i) + &
-                             ( gg - dd ) * de_dy(k,j+1,i+1) &
-                           ) / ( 3.0 * gg )
-             IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
-                de_dy_int = de_dy_int_l
-             ELSE
-                de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
-                                ( gg - bb ) * de_dy(k+1,j,i+1) + &
-                                ( gg - cc ) * de_dy(k+1,j+1,i) + &
-                                ( gg - dd ) * de_dy(k+1,j+1,i+1) &
-                              ) / ( 3.0 * gg )
-                de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * &
-                                          ( de_dy_int_u - de_dy_int_l )
-             ENDIF
-
-!
-!--          Interpolate the TKE gradient along z
-             de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
-                             ( gg - bb ) * de_dz(k,j,i+1) + &
-                             ( gg - cc ) * de_dz(k,j+1,i) + &
-                             ( gg - dd ) * de_dz(k,j+1,i+1) &
-                           ) / ( 3.0 * gg )
-
-             IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
-                de_dz_int = de_dz_int_l
-             ELSE
-                de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
-                                ( gg - bb ) * de_dz(k+1,j,i+1) + &
-                                ( gg - cc ) * de_dz(k+1,j+1,i) + &
-                                ( gg - dd ) * de_dz(k+1,j+1,i+1) &
-                              ) / ( 3.0 * gg )
-                de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * &
-                                          ( de_dz_int_u - de_dz_int_l )
-             ENDIF
-
-!
-!--          Interpolate the dissipation of TKE
-             diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
-                            ( gg - bb ) * diss(k,j,i+1) + &
-                            ( gg - cc ) * diss(k,j+1,i) + &
-                            ( gg - dd ) * diss(k,j+1,i+1) &
-                           ) / ( 3.0 * gg )
-
-             IF ( k+1 == nzt+1 )  THEN
-                diss_int = diss_int_l
-             ELSE
-                diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
-                               ( gg - bb ) * diss(k+1,j,i+1) + &
-                               ( gg - cc ) * diss(k+1,j+1,i) + &
-                               ( gg - dd ) * diss(k+1,j+1,i+1) &
-                              ) / ( 3.0 * gg )
-                diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * &
-                                        ( diss_int_u - diss_int_l )
-             ENDIF
+
+!
+!--             In case that there are buildings it has to be determined 
+!--             how many of the gridpoints defining the particle box are 
+!--             situated within a building
+!--             gp_outside_of_building(1): i,j,k
+!--             gp_outside_of_building(2): i,j+1,k
+!--             gp_outside_of_building(3): i,j,k+1
+!--             gp_outside_of_building(4): i,j+1,k+1
+!--             gp_outside_of_building(5): i+1,j,k
+!--             gp_outside_of_building(6): i+1,j+1,k
+!--             gp_outside_of_building(7): i+1,j,k+1
+!--             gp_outside_of_building(8): i+1,j+1,k+1
+           
+                gp_outside_of_building = 0
+                location = 0.0
+                num_gp = 0
+
+                IF ( k > nzb_s_inner(j,i)  .OR.  nzb_s_inner(j,i) == 0 )  THEN 
+                   num_gp = num_gp + 1
+                   gp_outside_of_building(1) = 1
+                   location(num_gp,1) = i * dx
+                   location(num_gp,2) = j * dx
+                   location(num_gp,3) = k * dz - 0.5 * dz
+                   ei(num_gp)     = e(k,j,i)
+                   dissi(num_gp)  = diss(k,j,i)
+                   de_dxi(num_gp) = de_dx(k,j,i)
+                   de_dyi(num_gp) = de_dy(k,j,i)
+                   de_dzi(num_gp) = de_dz(k,j,i)
+                ENDIF
+
+                IF ( k > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 ) &
+                THEN
+                   num_gp = num_gp + 1
+                   gp_outside_of_building(2) = 1
+                   location(num_gp,1) = i * dx
+                   location(num_gp,2) = (j+1) * dx
+                   location(num_gp,3) = k * dz - 0.5 * dz
+                   ei(num_gp)     = e(k,j+1,i)
+                   dissi(num_gp)  = diss(k,j+1,i)        
+                   de_dxi(num_gp) = de_dx(k,j+1,i)
+                   de_dyi(num_gp) = de_dy(k,j+1,i)
+                   de_dzi(num_gp) = de_dz(k,j+1,i)
+                ENDIF
+
+                IF ( k+1 > nzb_s_inner(j,i)  .OR.  nzb_s_inner(j,i) == 0 )  THEN
+                   num_gp = num_gp + 1
+                   gp_outside_of_building(3) = 1
+                   location(num_gp,1) = i * dx
+                   location(num_gp,2) = j * dy
+                   location(num_gp,3) = (k+1) * dz - 0.5 * dz
+                   ei(num_gp)     = e(k+1,j,i)
+                   dissi(num_gp)  = diss(k+1,j,i)        
+                   de_dxi(num_gp) = de_dx(k+1,j,i)
+                   de_dyi(num_gp) = de_dy(k+1,j,i)
+                   de_dzi(num_gp) = de_dz(k+1,j,i)
+                ENDIF
+
+                IF ( k+1 > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 ) &
+                THEN
+                   num_gp = num_gp + 1
+                   gp_outside_of_building(4) = 1
+                   location(num_gp,1) = i * dx
+                   location(num_gp,2) = (j+1) * dy
+                   location(num_gp,3) = (k+1) * dz - 0.5 * dz
+                   ei(num_gp)     = e(k+1,j+1,i)
+                   dissi(num_gp)  = diss(k+1,j+1,i)        
+                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
+                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
+                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
+                ENDIF
+ 
+                IF ( k > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 ) &
+                THEN
+                   num_gp = num_gp + 1
+                   gp_outside_of_building(5) = 1
+                   location(num_gp,1) = (i+1) * dx
+                   location(num_gp,2) = j * dy
+                   location(num_gp,3) = k * dz - 0.5 * dz
+                   ei(num_gp)     = e(k,j,i+1)
+                   dissi(num_gp)  = diss(k,j,i+1) 
+                   de_dxi(num_gp) = de_dx(k,j,i+1)
+                   de_dyi(num_gp) = de_dy(k,j,i+1)
+                   de_dzi(num_gp) = de_dz(k,j,i+1)
+                ENDIF
+
+                IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) &
+                THEN
+                   num_gp = num_gp + 1
+                   gp_outside_of_building(6) = 1
+                   location(num_gp,1) = (i+1) * dx
+                   location(num_gp,2) = (j+1) * dy
+                   location(num_gp,3) = k * dz - 0.5 * dz
+                   ei(num_gp)     = e(k,j+1,i+1)
+                   dissi(num_gp)  = diss(k,j+1,i+1) 
+                   de_dxi(num_gp) = de_dx(k,j+1,i+1)
+                   de_dyi(num_gp) = de_dy(k,j+1,i+1)
+                   de_dzi(num_gp) = de_dz(k,j+1,i+1)
+                ENDIF
+
+                IF ( k+1 > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 ) &
+                THEN
+                   num_gp = num_gp + 1
+                   gp_outside_of_building(7) = 1
+                   location(num_gp,1) = (i+1) * dx
+                   location(num_gp,2) = j * dy
+                   location(num_gp,3) = (k+1) * dz - 0.5 * dz
+                   ei(num_gp)     = e(k+1,j,i+1)
+                   dissi(num_gp)  = diss(k+1,j,i+1) 
+                   de_dxi(num_gp) = de_dx(k+1,j,i+1)
+                   de_dyi(num_gp) = de_dy(k+1,j,i+1)
+                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
+                ENDIF
+
+                IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0)&
+                THEN
+                   num_gp = num_gp + 1
+                   gp_outside_of_building(8) = 1
+                   location(num_gp,1) = (i+1) * dx
+                   location(num_gp,2) = (j+1) * dy
+                   location(num_gp,3) = (k+1) * dz - 0.5 * dz
+                   ei(num_gp)     = e(k+1,j+1,i+1)
+                   dissi(num_gp)  = diss(k+1,j+1,i+1) 
+                   de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
+                   de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
+                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
+                ENDIF
+
+!
+!--             If all gridpoints are situated outside of a building, then the 
+!--             ordinary interpolation scheme can be used.
+                IF ( num_gp == 8 )  THEN
+
+                   x  = particles(n)%x - i * dx
+                   y  = particles(n)%y - j * dy
+                   aa = x**2          + y**2
+                   bb = ( dx - x )**2 + y**2
+                   cc = x**2          + ( dy - y )**2
+                   dd = ( dx - x )**2 + ( dy - y )**2
+                   gg = aa + bb + cc + dd
+
+                   e_int_l = (( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1) &
+                            + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1)&
+                             ) / ( 3.0 * gg )
+
+                   IF ( k+1 == nzt+1 )  THEN
+                      e_int = e_int_l
+                   ELSE
+                      e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
+                                  ( gg - bb ) * e(k+1,j,i+1) + &
+                                  ( gg - cc ) * e(k+1,j+1,i) + &
+                                  ( gg - dd ) * e(k+1,j+1,i+1) &
+                                ) / ( 3.0 * gg )
+                      e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
+                                        ( e_int_u - e_int_l )
+                   ENDIF
+
+!
+!--                Interpolate the TKE gradient along x (adopt incides i,j,k
+!--                and all position variables from above (TKE))
+                   de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
+                                   ( gg - bb ) * de_dx(k,j,i+1) + &
+                                   ( gg - cc ) * de_dx(k,j+1,i) + &
+                                   ( gg - dd ) * de_dx(k,j+1,i+1) &
+                                 ) / ( 3.0 * gg )
+
+                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
+                      de_dx_int = de_dx_int_l
+                   ELSE
+                      de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
+                                      ( gg - bb ) * de_dx(k+1,j,i+1) + &
+                                      ( gg - cc ) * de_dx(k+1,j+1,i) + &
+                                      ( gg - dd ) * de_dx(k+1,j+1,i+1) &
+                                    ) / ( 3.0 * gg )
+                      de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / &
+                                              dz * ( de_dx_int_u - de_dx_int_l )
+                   ENDIF
+
+!
+!--                Interpolate the TKE gradient along y
+                   de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
+                                   ( gg - bb ) * de_dy(k,j,i+1) + &
+                                   ( gg - cc ) * de_dy(k,j+1,i) + &
+                                   ( gg - dd ) * de_dy(k,j+1,i+1) &
+                                 ) / ( 3.0 * gg )
+                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
+                      de_dy_int = de_dy_int_l
+                   ELSE
+                      de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
+                                      ( gg - bb ) * de_dy(k+1,j,i+1) + &
+                                      ( gg - cc ) * de_dy(k+1,j+1,i) + &
+                                      ( gg - dd ) * de_dy(k+1,j+1,i+1) &
+                                    ) / ( 3.0 * gg )
+                      de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / &
+                                              dz * ( de_dy_int_u - de_dy_int_l )
+                   ENDIF
+
+!
+!--                Interpolate the TKE gradient along z
+                   IF ( particles(n)%z < 0.5 * dz )  THEN
+                      de_dz_int = 0.0
+                   ELSE
+                      de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
+                                      ( gg - bb ) * de_dz(k,j,i+1) + &
+                                      ( gg - cc ) * de_dz(k,j+1,i) + &
+                                      ( gg - dd ) * de_dz(k,j+1,i+1) &
+                                    ) / ( 3.0 * gg )
+
+                      IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
+                         de_dz_int = de_dz_int_l
+                      ELSE
+                         de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
+                                         ( gg - bb ) * de_dz(k+1,j,i+1) + &
+                                         ( gg - cc ) * de_dz(k+1,j+1,i) + &
+                                         ( gg - dd ) * de_dz(k+1,j+1,i+1) &
+                                       ) / ( 3.0 * gg )
+                         de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) /&
+                                              dz * ( de_dz_int_u - de_dz_int_l )
+                      ENDIF
+                   ENDIF
+
+!
+!--                Interpolate the dissipation of TKE
+                   diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
+                                  ( gg - bb ) * diss(k,j,i+1) + &
+                                  ( gg - cc ) * diss(k,j+1,i) + &
+                                  ( gg - dd ) * diss(k,j+1,i+1) &
+                                ) / ( 3.0 * gg )
+
+                   IF ( k+1 == nzt+1 )  THEN
+                      diss_int = diss_int_l
+                   ELSE
+                      diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
+                                     ( gg - bb ) * diss(k+1,j,i+1) + &
+                                     ( gg - cc ) * diss(k+1,j+1,i) + &
+                                     ( gg - dd ) * diss(k+1,j+1,i+1) &
+                                   ) / ( 3.0 * gg )
+                      diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz *&
+                                              ( diss_int_u - diss_int_l )
+                   ENDIF
+
+                ELSE
+
+!
+!--                If wall between gridpoint 1 and gridpoint 5, then
+!--                Neumann boundary condition has to be applied
+                   IF ( gp_outside_of_building(1) == 1  .AND. &
+                        gp_outside_of_building(5) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx + 0.5 * dx
+                      location(num_gp,2) = j * dy
+                      location(num_gp,3) = k * dz - 0.5 * dz
+                      ei(num_gp)     = e(k,j,i)
+                      dissi(num_gp)  = diss(k,j,i)
+                      de_dxi(num_gp) = 0.0
+                      de_dyi(num_gp) = de_dy(k,j,i)
+                      de_dzi(num_gp) = de_dz(k,j,i)                     
+                   ENDIF 
+   
+                   IF ( gp_outside_of_building(5) == 1  .AND. &
+                      gp_outside_of_building(1) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx + 0.5 * dx
+                      location(num_gp,2) = j * dy
+                      location(num_gp,3) = k * dz - 0.5 * dz
+                      ei(num_gp)     = e(k,j,i+1)
+                      dissi(num_gp)  = diss(k,j,i+1)
+                      de_dxi(num_gp) = 0.0
+                      de_dyi(num_gp) = de_dy(k,j,i+1)
+                      de_dzi(num_gp) = de_dz(k,j,i+1)
+                   ENDIF
+
+!
+!--                If wall between gridpoint 5 and gridpoint 6, then
+!--                then Neumann boundary condition has to be applied
+                   IF ( gp_outside_of_building(5) == 1  .AND. &
+                        gp_outside_of_building(6) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = (i+1) * dx
+                      location(num_gp,2) = j * dy + 0.5 * dy
+                      location(num_gp,3) = k * dz - 0.5 * dz
+                      ei(num_gp)     = e(k,j,i+1)
+                      dissi(num_gp)  = diss(k,j,i+1)
+                      de_dxi(num_gp) = de_dx(k,j,i+1)
+                      de_dyi(num_gp) = 0.0
+                      de_dzi(num_gp) = de_dz(k,j,i+1)
+                   ENDIF
+
+                   IF ( gp_outside_of_building(6) == 1  .AND. &
+                        gp_outside_of_building(5) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = (i+1) * dx
+                      location(num_gp,2) = j * dy + 0.5 * dy
+                      location(num_gp,3) = k * dz - 0.5 * dz
+                      ei(num_gp)     = e(k,j+1,i+1)
+                      dissi(num_gp)  = diss(k,j+1,i+1)
+                      de_dxi(num_gp) = de_dx(k,j+1,i+1)
+                      de_dyi(num_gp) = 0.0
+                      de_dzi(num_gp) = de_dz(k,j+1,i+1)
+                   ENDIF
+
+!
+!--                If wall between gridpoint 2 and gridpoint 6, then
+!--                Neumann boundary condition has to be applied
+                   IF ( gp_outside_of_building(2) == 1  .AND. &
+                        gp_outside_of_building(6) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx + 0.5 * dx
+                      location(num_gp,2) = (j+1) * dy
+                      location(num_gp,3) = k * dz - 0.5 * dz
+                      ei(num_gp)     = e(k,j+1,i)
+                      dissi(num_gp)  = diss(k,j+1,i)
+                      de_dxi(num_gp) = 0.0
+                      de_dyi(num_gp) = de_dy(k,j+1,i)
+                      de_dzi(num_gp) = de_dz(k,j+1,i)                     
+                   ENDIF 
+   
+                   IF ( gp_outside_of_building(6) == 1  .AND. &
+                        gp_outside_of_building(2) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx + 0.5 * dx
+                      location(num_gp,2) = (j+1) * dy
+                      location(num_gp,3) = k * dz - 0.5 * dz
+                      ei(num_gp)     = e(k,j+1,i+1)
+                      dissi(num_gp)  = diss(k,j+1,i+1)
+                      de_dxi(num_gp) = 0.0
+                      de_dyi(num_gp) = de_dy(k,j+1,i+1)
+                      de_dzi(num_gp) = de_dz(k,j+1,i+1)
+                   ENDIF
+
+!
+!--                If wall between gridpoint 1 and gridpoint 2, then
+!--                Neumann boundary condition has to be applied
+                   IF ( gp_outside_of_building(1) == 1  .AND. &
+                        gp_outside_of_building(2) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx
+                      location(num_gp,2) = j * dy + 0.5 * dy
+                      location(num_gp,3) = k * dz - 0.5 * dz
+                      ei(num_gp)     = e(k,j,i)
+                      dissi(num_gp)  = diss(k,j,i)
+                      de_dxi(num_gp) = de_dx(k,j,i)
+                      de_dyi(num_gp) = 0.0
+                      de_dzi(num_gp) = de_dz(k,j,i)
+                   ENDIF
+
+                   IF ( gp_outside_of_building(2) == 1  .AND. &
+                        gp_outside_of_building(1) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx
+                      location(num_gp,2) = j * dy + 0.5 * dy
+                      location(num_gp,3) = k * dz - 0.5 * dz
+                      ei(num_gp)     = e(k,j+1,i)
+                      dissi(num_gp)  = diss(k,j+1,i)
+                      de_dxi(num_gp) = de_dx(k,j+1,i)
+                      de_dyi(num_gp) = 0.0
+                      de_dzi(num_gp) = de_dz(k,j+1,i)
+                   ENDIF
+
+!
+!--                If wall between gridpoint 3 and gridpoint 7, then
+!--                Neumann boundary condition has to be applied
+                   IF ( gp_outside_of_building(3) == 1  .AND. &
+                        gp_outside_of_building(7) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx + 0.5 * dx
+                      location(num_gp,2) = j * dy
+                      location(num_gp,3) = k * dz + 0.5 * dz
+                      ei(num_gp)     = e(k+1,j,i)
+                      dissi(num_gp)  = diss(k+1,j,i)
+                      de_dxi(num_gp) = 0.0
+                      de_dyi(num_gp) = de_dy(k+1,j,i)
+                      de_dzi(num_gp) = de_dz(k+1,j,i)                     
+                   ENDIF 
+   
+                   IF ( gp_outside_of_building(7) == 1  .AND. &
+                        gp_outside_of_building(3) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx + 0.5 * dx
+                      location(num_gp,2) = j * dy
+                      location(num_gp,3) = k * dz + 0.5 * dz
+                      ei(num_gp)     = e(k+1,j,i+1)
+                      dissi(num_gp)  = diss(k+1,j,i+1)
+                      de_dxi(num_gp) = 0.0
+                      de_dyi(num_gp) = de_dy(k+1,j,i+1)
+                      de_dzi(num_gp) = de_dz(k+1,j,i+1)
+                   ENDIF
+
+!
+!--                If wall between gridpoint 7 and gridpoint 8, then
+!--                Neumann boundary condition has to be applied
+                   IF ( gp_outside_of_building(7) == 1  .AND. &
+                        gp_outside_of_building(8) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = (i+1) * dx
+                      location(num_gp,2) = j * dy + 0.5 * dy
+                      location(num_gp,3) = k * dz + 0.5 * dz
+                      ei(num_gp)     = e(k+1,j,i+1)
+                      dissi(num_gp)  = diss(k+1,j,i+1)
+                      de_dxi(num_gp) = de_dx(k+1,j,i+1)
+                      de_dyi(num_gp) = 0.0
+                      de_dzi(num_gp) = de_dz(k+1,j,i+1)
+                   ENDIF
+
+                   IF ( gp_outside_of_building(8) == 1  .AND. &
+                        gp_outside_of_building(7) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = (i+1) * dx
+                      location(num_gp,2) = j * dy + 0.5 * dy
+                      location(num_gp,3) = k * dz + 0.5 * dz
+                      ei(num_gp)     = e(k+1,j+1,i+1)
+                      dissi(num_gp)  = diss(k+1,j+1,i+1)
+                      de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
+                      de_dyi(num_gp) = 0.0
+                      de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
+                   ENDIF
+
+!
+!--                If wall between gridpoint 4 and gridpoint 8, then
+!--                Neumann boundary condition has to be applied
+                   IF ( gp_outside_of_building(4) == 1  .AND. &
+                        gp_outside_of_building(8) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx + 0.5 * dx
+                      location(num_gp,2) = (j+1) * dy
+                      location(num_gp,3) = k * dz + 0.5 * dz
+                      ei(num_gp)     = e(k+1,j+1,i)
+                      dissi(num_gp)  = diss(k+1,j+1,i)
+                      de_dxi(num_gp) = 0.0
+                      de_dyi(num_gp) = de_dy(k+1,j+1,i)
+                      de_dzi(num_gp) = de_dz(k+1,j+1,i)                     
+                   ENDIF 
+   
+                   IF ( gp_outside_of_building(8) == 1  .AND. &
+                        gp_outside_of_building(4) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx + 0.5 * dx
+                      location(num_gp,2) = (j+1) * dy
+                      location(num_gp,3) = k * dz + 0.5 * dz
+                      ei(num_gp)     = e(k+1,j+1,i+1)
+                      dissi(num_gp)  = diss(k+1,j+1,i+1)
+                      de_dxi(num_gp) = 0.0
+                      de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
+                      de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
+                   ENDIF
+
+!
+!--                If wall between gridpoint 3 and gridpoint 4, then
+!--                Neumann boundary condition has to be applied
+                   IF ( gp_outside_of_building(3) == 1  .AND. &
+                        gp_outside_of_building(4) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx
+                      location(num_gp,2) = j * dy + 0.5 * dy
+                      location(num_gp,3) = k * dz + 0.5 * dz
+                      ei(num_gp)     = e(k+1,j,i)
+                      dissi(num_gp)  = diss(k+1,j,i)
+                      de_dxi(num_gp) = de_dx(k+1,j,i)
+                      de_dyi(num_gp) = 0.0
+                      de_dzi(num_gp) = de_dz(k+1,j,i)
+                   ENDIF
+
+                   IF ( gp_outside_of_building(4) == 1  .AND. &
+                        gp_outside_of_building(3) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx
+                      location(num_gp,2) = j * dy + 0.5 * dy
+                      location(num_gp,3) = k * dz + 0.5 * dz
+                      ei(num_gp)     = e(k+1,j+1,i)
+                      dissi(num_gp)  = diss(k+1,j+1,i)
+                      de_dxi(num_gp) = de_dx(k+1,j+1,i)
+                      de_dyi(num_gp) = 0.0
+                      de_dzi(num_gp) = de_dz(k+1,j+1,i)
+                   ENDIF
+
+!
+!--                If wall between gridpoint 1 and gridpoint 3, then
+!--                Neumann boundary condition has to be applied
+!--                (only one case as only building beneath is possible)
+                   IF ( gp_outside_of_building(1) == 1  .AND. &
+                        gp_outside_of_building(3) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx
+                      location(num_gp,2) = j * dy
+                      location(num_gp,3) = k * dz
+                      ei(num_gp)     = e(k+1,j,i)
+                      dissi(num_gp)  = diss(k+1,j,i)
+                      de_dxi(num_gp) = de_dx(k+1,j,i)
+                      de_dyi(num_gp) = de_dy(k+1,j,i)
+                      de_dzi(num_gp) = 0.0
+                   ENDIF
+ 
+!
+!--                If wall between gridpoint 5 and gridpoint 7, then
+!--                Neumann boundary condition has to be applied
+!--                (only one case as only building beneath is possible)
+                   IF ( gp_outside_of_building(5) == 1  .AND. &
+                        gp_outside_of_building(7) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = (i+1) * dx
+                      location(num_gp,2) = j * dy
+                      location(num_gp,3) = k * dz
+                      ei(num_gp)     = e(k+1,j,i+1)
+                      dissi(num_gp)  = diss(k+1,j,i+1)
+                      de_dxi(num_gp) = de_dx(k+1,j,i+1)
+                      de_dyi(num_gp) = de_dy(k+1,j,i+1)
+                      de_dzi(num_gp) = 0.0
+                   ENDIF
+
+!
+!--                If wall between gridpoint 2 and gridpoint 4, then
+!--                Neumann boundary condition has to be applied
+!--                (only one case as only building beneath is possible)
+                   IF ( gp_outside_of_building(2) == 1  .AND. &
+                        gp_outside_of_building(4) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = i * dx
+                      location(num_gp,2) = (j+1) * dy
+                      location(num_gp,3) = k * dz
+                      ei(num_gp)     = e(k+1,j+1,i)
+                      dissi(num_gp)  = diss(k+1,j+1,i)
+                      de_dxi(num_gp) = de_dx(k+1,j+1,i)
+                      de_dyi(num_gp) = de_dy(k+1,j+1,i)
+                      de_dzi(num_gp) = 0.0
+                   ENDIF
+
+!
+!--                If wall between gridpoint 6 and gridpoint 8, then
+!--                Neumann boundary condition has to be applied
+!--                (only one case as only building beneath is possible)
+                   IF ( gp_outside_of_building(6) == 1  .AND. &
+                        gp_outside_of_building(8) == 0 )  THEN
+                      num_gp = num_gp + 1
+                      location(num_gp,1) = (i+1) * dx
+                      location(num_gp,2) = (j+1) * dy
+                      location(num_gp,3) = k * dz
+                      ei(num_gp)     = e(k+1,j+1,i+1)
+                      dissi(num_gp)  = diss(k+1,j+1,i+1)
+                      de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
+                      de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
+                      de_dzi(num_gp) = 0.0
+                   ENDIF
+
+!
+!--                Carry out the interpolation
+                   IF ( num_gp == 1 )  THEN
+!
+!--                   If only one of the gridpoints is situated outside of the 
+!--                   building, it follows that the values at the particle 
+!--                   location are the same as the gridpoint values
+                      e_int = ei(num_gp)
+
+                   ELSE IF ( num_gp > 1 )  THEN
+
+                      d_sum = 0.0
+!
+!--                   Evaluation of the distances between the gridpoints 
+!--                   contributing to the interpolated values, and the particle 
+!--                   location
+                      DO agp = 1, num_gp
+                         d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2  &
+                                      + ( particles(n)%y-location(agp,2) )**2  &
+                                      + ( particles(n)%z-location(agp,3) )**2
+                         d_sum        = d_sum + d_gp_pl(agp)
+                      ENDDO
+
+!
+!--                   Finally the interpolation can be carried out
+                      e_int     = 0.0
+                      diss_int  = 0.0
+                      de_dx_int = 0.0
+                      de_dy_int = 0.0
+                      de_dz_int = 0.0
+                      DO agp = 1, num_gp
+                         e_int     = e_int     + ( d_sum - d_gp_pl(agp) ) * &
+                                                ei(agp) / ( (num_gp-1) * d_sum )
+                         diss_int  = diss_int  + ( d_sum - d_gp_pl(agp) ) * &
+                                             dissi(agp) / ( (num_gp-1) * d_sum )
+                         de_dx_int = de_dx_int + ( d_sum - d_gp_pl(agp) ) * &
+                                            de_dxi(agp) / ( (num_gp-1) * d_sum )
+                         de_dy_int = de_dy_int + ( d_sum - d_gp_pl(agp) ) * &
+                                            de_dyi(agp) / ( (num_gp-1) * d_sum )
+                         de_dz_int = de_dz_int + ( d_sum - d_gp_pl(agp) ) * &
+                                            de_dzi(agp) / ( (num_gp-1) * d_sum )
+                      ENDDO
+
+                   ENDIF
+
+                ENDIF
+
+             ENDIF  
 
 !
@@ -1148,13 +1681,36 @@
 
              ELSE
+
+!
+!--             Restriction of the size of the new timestep: compared to the 
+!--             previous timestep the increase must not exceed 200%
+
+                dt_particle_old = particles(n)%age - particles(n)%age_m
+                IF ( dt_particle > 2.0 * dt_particle_old )  THEN 
+                   dt_particle = 2.0 * dt_particle_old
+                ENDIF
+
 !
 !--             For old particles the SGS components are correlated with the
 !--             values from the previous timestep. Random numbers have also to
 !--             be limited (see above).
+!--             As negative values for the subgrid TKE are not allowed, the
+!--             change of the subgrid TKE with time cannot be smaller than
+!--             -e_int/dt_particle. This value is used as a lower boundary
+!--             value for the change of TKE 
+
+                de_dt_min = - e_int / dt_particle
+
+                de_dt = ( e_int - particles(n)%e_m ) / dt_particle_old
+
+                IF ( de_dt < de_dt_min )  THEN
+                   de_dt = de_dt_min
+                ENDIF
+
                 particles(n)%speed_x_sgs = particles(n)%speed_x_sgs -          &
                        fs_int * c_0 * diss_int * particles(n)%speed_x_sgs *    &
                        dt_particle / ( 4.0 * sgs_wfu_part * e_int ) +          &
-                       ( ( 2.0 * sgs_wfu_part * ( e_int - particles(n)%e_m ) / &
-                         dt_particle ) * particles(n)%speed_x_sgs /            &
+                       ( 2.0 * sgs_wfu_part * de_dt *                          &
+                               particles(n)%speed_x_sgs /                      &
                          ( 2.0 * sgs_wfu_part * e_int ) + de_dx_int            &
                        ) * dt_particle / 2.0                        +          &
@@ -1166,6 +1722,6 @@
                        fs_int * c_0 * diss_int * particles(n)%speed_y_sgs *    &
                        dt_particle / ( 4.0 * sgs_wfv_part * e_int ) +          &
-                       ( ( 2.0 * sgs_wfv_part * ( e_int - particles(n)%e_m ) / &
-                         dt_particle ) * particles(n)%speed_y_sgs /            &
+                       ( 2.0 * sgs_wfv_part * de_dt *                          &
+                               particles(n)%speed_y_sgs /                      &
                          ( 2.0 * sgs_wfv_part * e_int ) + de_dy_int            &
                        ) * dt_particle / 2.0                        +          &
@@ -1177,6 +1733,6 @@
                        fs_int * c_0 * diss_int * particles(n)%speed_z_sgs *    &
                        dt_particle / ( 4.0 * sgs_wfw_part * e_int ) +          &
-                       ( ( 2.0 * sgs_wfw_part * ( e_int - particles(n)%e_m ) / &
-                         dt_particle ) * particles(n)%speed_z_sgs /            &
+                       ( 2.0 * sgs_wfw_part * de_dt *                          &
+                               particles(n)%speed_z_sgs /                      &
                          ( 2.0 * sgs_wfw_part * e_int ) + de_dz_int            &
                        ) * dt_particle / 2.0                        +          &
@@ -1203,4 +1759,11 @@
 
           ENDIF
+
+!
+!--       Remember the old age of the particle ( needed to prevent that a
+!--       particle crosses several PEs during one timestep and for the
+!--       evaluation of the subgrid particle velocity fluctuations )
+          particles(n)%age_m = particles(n)%age
+
 
 !
@@ -1270,4 +1833,689 @@
 
        ENDDO   ! advection loop
+
+
+!
+!--    Particle reflection from walls
+       CALL cpu_log( log_point_s(48), 'advec_part_refle', 'start' )
+
+       DO  n = 1, number_of_particles
+
+!-------------------------------------------------------
+      i2 =  (particles(n)%x + 0.5*dx) * ddx
+      j2 =  (particles(n)%y + 0.5*dy) * ddy
+      k2 =  particles(n)%z / dz + 1
+      
+      particles_x = particles(n)%x
+      particles_y = particles(n)%y
+      particles_z = particles(n)%z
+
+ 
+    IF (k2<=nzb_s_inner(j2,i2).and.nzb_s_inner(j2,i2)/=0) THEN
+
+   old_positions_x = particles(n)%x - particles(n)%speed_x * dt_particle
+   old_positions_y = particles(n)%y - particles(n)%speed_y * dt_particle
+   old_positions_z = particles(n)%z - particles(n)%speed_z * dt_particle
+   i1= (old_positions_x +0.5*dx) * ddx
+   j1 =(old_positions_y +0.5*dy) * ddy
+   k1 = old_positions_z / dz
+    
+
+!-- CASE 1
+
+   IF ( particles(n)%x  >  old_positions_x .and. &
+         particles(n)%y  >  old_positions_y )THEN
+ 
+     
+
+
+     t_index = 1
+   
+     Do i = i1, i2
+     xline = i*dx+0.5*dx
+     
+     t(t_index) =(xline-old_positions_x)/(particles(n)%x-old_positions_x)
+     t_index = t_index + 1
+     ENDDO
+
+     Do j = j1, j2
+     yline = j*dy+0.5*dy
+     
+     t(t_index) =(yline-old_positions_y)/(particles(n)%y-old_positions_y)
+     t_index = t_index + 1
+     ENDDO
+
+     IF ( particles(n)%z < old_positions_z ) THEN
+     Do k = k1, k2 , -1
+     zline = k*dz
+     t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z)
+     t_index = t_index + 1
+     ENDDO
+     ENDIF
+     t_index_number = t_index-1
+!-- sorting t
+                   inc = 1
+                   jr  = 1
+                   DO WHILE ( inc <= t_index_number )
+                      inc = 3 * inc + 1
+                   ENDDO
+
+                   DO WHILE ( inc > 1 )
+                      inc = inc / 3
+                      DO  ir = inc+1, t_index_number
+                         tmp_t = t(ir)
+                         jr = ir
+                         
+                         DO WHILE ( t(jr-inc) > tmp_t )
+                            t(jr) = t(jr-inc)
+                            jr = jr - inc
+                            IF ( jr <= inc )  THEN
+                           EXIT
+                            ENDIF
+                         ENDDO
+                        t(jr) = tmp_t
+                       
+             ENDDO
+        ENDDO
+    
+
+   
+   
+    Do t_index = 1, t_index_number
+  
+    positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x)
+    positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y)
+    positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z)
+
+     
+    
+     i3 = (positions_x + 0.5*dx) * ddx   
+     j3 = (positions_y + 0.5*dy) * ddy
+     k3 =  positions_z / dz
+
+     i5 = positions_x * ddx
+     j5 = positions_y * ddy
+     k5 = positions_z / dz
+
+     
+ 
+
+
+    IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN
+             
+         IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. &
+          k3 == nzb_s_inner(j5,i3) ) THEN
+      
+          particles(n)%z = 2*positions_z - particles_z
+          particles(n)%speed_z = - particles(n)%speed_z
+          
+           IF ( use_sgs_for_particles )  THEN
+             particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+           ENDIF
+          GOTO 999
+       ENDIF
+    ENDIF
+   
+IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN
+             
+         IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. &
+          k3 == nzb_s_inner(j3,i5) ) THEN
+      
+          particles(n)%z = 2*positions_z - particles_z
+          particles(n)%speed_z = - particles(n)%speed_z
+            
+            IF ( use_sgs_for_particles )  THEN
+              particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+            ENDIF
+           
+          GOTO 999
+       ENDIF
+    ENDIF
+  
+     IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN
+
+       IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. &
+          k3 == nzb_s_inner(j3,i3) ) THEN
+      
+          particles(n)%z = 2*positions_z - particles_z
+          particles(n)%speed_z = - particles(n)%speed_z
+           
+           IF ( use_sgs_for_particles )  THEN
+          particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+           ENDIF
+           
+          GOTO 999
+       ENDIF
+    
+       IF ( positions_y == j3*dy-0.5*dy .and. &
+         positions_z < nzb_s_inner(j3,i3)*dz) THEN
+              
+          particles(n)%y = 2*positions_y - particles_y
+          particles(n)%speed_y = - particles(n)%speed_y
+          
+          IF ( use_sgs_for_particles )  THEN
+            particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs
+          ENDIF
+   
+          GOTO 999
+       ENDIF
+
+  
+       IF ( positions_x == i3*dx-0.5*dx.and. &
+            positions_z < nzb_s_inner(j3,i3)*dz ) THEN
+              
+          particles(n)%x = 2*positions_x - particles_x
+          particles(n)%speed_x = - particles(n)%speed_x
+       
+        IF ( use_sgs_for_particles )  THEN 
+          particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs
+          GOTO 999
+        ENDIF  
+        
+        ENDIF
+     
+     ENDIF 
+    ENDDO
+   ENDIF
+    
+!-- CASE 2
+  
+     IF ( particles(n)%x  >  old_positions_x .and. &
+         particles(n)%y  <  old_positions_y )THEN
+ 
+
+     t_index = 1
+   
+     Do i = i1, i2
+     xline = i*dx+0.5*dx
+     
+     t(t_index) =(xline-old_positions_x)/(particles(n)%x-old_positions_x)
+     t_index = t_index + 1
+     ENDDO
+
+     Do j = j1, j2,-1
+     yline = j*dy-0.5*dy
+     
+     t(t_index) =(old_positions_y-yline)/(old_positions_y-particles(n)%y)
+     t_index = t_index + 1
+     ENDDO
+
+     IF ( particles(n)%z < old_positions_z ) THEN
+     Do k = k1, k2 , -1
+     zline = k*dz
+     t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z)
+     t_index = t_index + 1
+     ENDDO
+     ENDIF
+     t_index_number = t_index-1
+
+!-- sorting t
+                   inc = 1
+                   jr  = 1
+                   DO WHILE ( inc <= t_index_number )
+                      inc = 3 * inc + 1
+                   ENDDO
+
+                   DO WHILE ( inc > 1 )
+                      inc = inc / 3
+                      DO  ir = inc+1, t_index_number
+                         tmp_t = t(ir)
+                         jr = ir
+                         
+                         DO WHILE ( t(jr-inc) > tmp_t )
+                            t(jr) = t(jr-inc)
+                            jr = jr - inc
+                            IF ( jr <= inc )  THEN
+                           EXIT
+                            ENDIF
+                         ENDDO
+                        t(jr) = tmp_t
+                       
+             ENDDO
+        ENDDO
+
+     
+  
+     
+    
+   
+    Do t_index = 1, t_index_number
+  
+    positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x)
+    positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y)
+    positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z)
+
+     
+    
+     i3 = (positions_x + 0.5*dx) * ddx   
+     j3 = (positions_y + 0.5*dy) * ddy
+     k3 =  positions_z / dz
+
+     i5 = positions_x * ddx
+     j5 = positions_y * ddy
+     k5 = positions_z / dz
+    
+
+ 
+    IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN
+             
+         IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. &
+          k3 == nzb_s_inner(j3,i5) ) THEN
+      
+          particles(n)%z = 2*positions_z - particles_z
+          particles(n)%speed_z = - particles(n)%speed_z
+           
+           IF ( use_sgs_for_particles )  THEN
+            particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+           ENDIF
+       
+          GOTO 999
+       ENDIF
+    ENDIF
+     
+     IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN
+
+       IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. &
+          k3 == nzb_s_inner(j3,i3) ) THEN
+      
+          particles(n)%z = 2*positions_z - particles_z
+          particles(n)%speed_z = - particles(n)%speed_z
+          
+         IF ( use_sgs_for_particles )  THEN
+          particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+         ENDIF 
+         
+         GOTO 999
+       ENDIF
+  
+        
+       IF ( positions_x == i3*dx-0.5*dx.and.&
+       positions_z < nzb_s_inner(j3,i3)*dz ) THEN
+              
+          particles(n)%x = 2*positions_x - particles_x
+          particles(n)%speed_x = - particles(n)%speed_x
+          
+          IF ( use_sgs_for_particles )  THEN
+           particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs
+          ENDIF
+          
+          GOTO 999
+       ENDIF
+     
+     ENDIF
+
+     
+     IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN
+
+       IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. &
+          k3 == nzb_s_inner(j5,i3) ) THEN
+      
+          particles(n)%z = 2*positions_z - particles_z
+          particles(n)%speed_z = - particles(n)%speed_z
+           IF ( use_sgs_for_particles )  THEN
+            particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+           ENDIF
+
+          GOTO 999
+       ENDIF
+  
+       IF ( positions_y == j5 * dy + 0.5*dy.and.&
+        positions_z < nzb_s_inner(j5,i3)*dz ) THEN
+      
+          particles(n)%y = 2*positions_y - particles_y
+          particles(n)%speed_y = - particles(n)%speed_y
+           
+          IF ( use_sgs_for_particles )  THEN
+            particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs
+          ENDIF
+          GOTO 999
+       ENDIF
+
+     ENDIF 
+    ENDDO
+   ENDIF
+
+
+!-- CASE 3
+
+
+
+        IF ( particles(n)%x < old_positions_x .and. &
+         particles(n)%y > old_positions_y )THEN
+ 
+
+     t_index = 1
+   
+     Do i = i1, i2,-1
+     xline = i*dx-0.5*dx
+     
+     t(t_index) =(old_positions_x-xline)/(old_positions_x-particles(n)%x)
+     t_index = t_index + 1
+     ENDDO
+
+     Do j = j1, j2
+     yline = j*dy+0.5*dy
+     
+     t(t_index) =(yline-old_positions_y)/(particles(n)%y-old_positions_y)
+     t_index = t_index + 1
+     ENDDO
+     
+     IF ( particles(n)%z < old_positions_z ) THEN
+     Do k = k1, k2 , -1
+     zline = k*dz
+     t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z)
+     t_index = t_index + 1
+     ENDDO
+     ENDIF 
+    
+     t_index_number = t_index-1
+
+!-- sorting t
+                   inc = 1
+                   jr  = 1
+                   DO WHILE ( inc <= t_index_number )
+                      inc = 3 * inc + 1
+                   ENDDO
+
+                   DO WHILE ( inc > 1 )
+                      inc = inc / 3
+                      DO  ir = inc+1, t_index_number
+                         tmp_t = t(ir)
+                         jr = ir
+                         
+                         DO WHILE ( t(jr-inc) > tmp_t )
+                            t(jr) = t(jr-inc)
+                            jr = jr - inc
+                            IF ( jr <= inc )  THEN
+                           EXIT
+                            ENDIF
+                         ENDDO
+                        t(jr) = tmp_t
+                       
+             ENDDO
+        ENDDO
+
+    
+   
+    Do t_index = 1, t_index_number
+  
+    positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x)
+    positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y)
+    positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z)
+
+     
+    
+     i3 = (positions_x + 0.5*dx) * ddx   
+     j3 = (positions_y + 0.5*dy) * ddy
+     k3 =  positions_z / dz
+
+     i5 = positions_x * ddx
+     j5 = positions_y * ddy
+     k5 = positions_z / dz
+
+
+
+     IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN
+             
+         IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. &
+          k3 == nzb_s_inner(j5,i3) ) THEN
+      
+          particles(n)%z = 2*positions_z - particles_z
+          particles(n)%speed_z = - particles(n)%speed_z
+         
+         IF ( use_sgs_for_particles )  THEN  
+          particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+         ENDIF
+
+          GOTO 999
+       ENDIF
+    ENDIF
+     
+
+
+      IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN
+
+       
+      IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. &
+          k3 == nzb_s_inner(j3,i3) ) THEN
+      
+          particles(n)%z = 2*positions_z - particles_z
+          particles(n)%speed_z = - particles(n)%speed_z
+          
+         IF ( use_sgs_for_particles )  THEN
+          particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+         ENDIF
+
+          GOTO 999
+       ENDIF
+
+          IF ( positions_y == j3*dy-0.5*dy .and. &
+            positions_z < nzb_s_inner(j3,i3)*dz) THEN
+              
+          particles(n)%y = 2*positions_y - particles_y
+          particles(n)%speed_y = - particles(n)%speed_y
+          
+           IF ( use_sgs_for_particles )  THEN
+              particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs
+           ENDIF
+
+          GOTO 999
+       ENDIF  
+
+     ENDIF
+      
+        
+     IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN
+
+       IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. &
+          k3 == nzb_s_inner(j3,i5) ) THEN
+      
+          particles(n)%z = 2*positions_z - particles_z
+          particles(n)%speed_z = - particles(n)%speed_z
+          
+         IF ( use_sgs_for_particles )  THEN
+          particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+         ENDIF 
+
+
+          GOTO 999
+       ENDIF
+  
+       IF ( positions_x == i5 * dx + 0.5*dx .and. &
+               positions_z < nzb_s_inner(j3,i5)*dz ) THEN
+      
+          particles(n)%x= 2*positions_x - particles_x
+          particles(n)%speed_x = - particles(n)%speed_x
+          
+         IF ( use_sgs_for_particles )  THEN
+           particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs
+         ENDIF  
+        
+          GOTO 999
+       ENDIF
+
+     ENDIF 
+    ENDDO
+   ENDIF
+
+
+!-- CASE 4
+     
+
+       IF ( particles(n)%x < old_positions_x .and. &
+         particles(n)%y  < old_positions_y )THEN
+ 
+
+     t_index = 1
+   
+     Do i = i1, i2,-1
+     xline = i*dx-0.5*dx
+     
+     t(t_index) =(old_positions_x-xline)/(old_positions_x-particles(n)%x)
+     t_index = t_index + 1
+     
+     
+     
+     ENDDO
+
+     Do j = j1, j2,-1
+     yline = j*dy-0.5*dy
+     
+     t(t_index) =(old_positions_y-yline)/(old_positions_y-particles(n)%y)
+     t_index = t_index + 1
+     
+     
+     ENDDO
+
+     IF ( particles(n)%z < old_positions_z ) THEN
+     Do k = k1, k2 , -1
+     zline = k*dz
+     t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z)
+     t_index = t_index + 1
+     
+
+     ENDDO
+     ENDIF
+
+     t_index_number = t_index-1
+
+!-- sorting t
+                   inc = 1
+                   jr  = 1
+                   DO WHILE ( inc <= t_index_number )
+                      inc = 3 * inc + 1
+                   ENDDO
+
+                   DO WHILE ( inc > 1 )
+                      inc = inc / 3
+                      DO  ir = inc+1, t_index_number
+                         tmp_t = t(ir)
+                         jr = ir
+                         
+                         DO WHILE ( t(jr-inc) > tmp_t )
+                            t(jr) = t(jr-inc)
+                            jr = jr - inc
+                            IF ( jr <= inc )  THEN
+                           EXIT
+                            ENDIF
+                         ENDDO
+                        t(jr) = tmp_t
+                       
+             ENDDO
+        ENDDO
+    
+    
+    Do t_index = 1, t_index_number
+  
+    positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x)
+    positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y)
+    positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z)
+
+     
+    
+     i3 = (positions_x + 0.5*dx) * ddx   
+     j3 = (positions_y + 0.5*dy) * ddy
+     k3 =  positions_z / dz
+
+     i5 = positions_x * ddx
+     j5 = positions_y * ddy
+     k5 = positions_z / dz
+
+
+
+
+     IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN
+
+       IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. &
+          k3 == nzb_s_inner(j3,i3) ) THEN
+      
+          particles(n)%z = 2*positions_z - particles_z
+          particles(n)%speed_z = - particles(n)%speed_z
+         
+          IF ( use_sgs_for_particles )  THEN
+           particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+          ENDIF 
+         
+
+          GOTO 999
+       ENDIF
+     ENDIF
+           
+        
+     IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN
+
+       IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. &
+          k3 == nzb_s_inner(j3,i5) ) THEN
+      
+          particles(n)%z = 2*positions_z - particles_z
+          particles(n)%speed_z = - particles(n)%speed_z
+         
+         IF ( use_sgs_for_particles )  THEN 
+          particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+         ENDIF  
+        
+          GOTO 999
+       ENDIF
+  
+       IF ( positions_x == i5 * dx + 0.5*dx .and. nzb_s_inner(j3,i5)/=0.and.&
+                 positions_z < nzb_s_inner(j3,i5)*dz) THEN
+      
+          particles(n)%x= 2*positions_x - particles_x
+          particles(n)%speed_x = - particles(n)%speed_x
+         
+         IF ( use_sgs_for_particles )  THEN  
+           particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs
+         ENDIF 
+          GOTO 999
+       ENDIF
+      
+    ENDIF
+    
+     IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0) THEN
+       
+
+       IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. &
+          k5 == nzb_s_inner(j5,i3) ) THEN
+      
+          particles(n)%z = 2*positions_z - particles_z
+          particles(n)%speed_z = - particles(n)%speed_z
+          
+          IF ( use_sgs_for_particles )  THEN 
+           particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
+          ENDIF
+          
+          GOTO 999
+       ENDIF
+
+          IF ( positions_y == j5 * dy + 0.5*dy .and. nzb_s_inner(j5,i3)/=0.and.&
+                    positions_z < nzb_s_inner(j5,i3)*dz ) THEN
+      
+          particles(n)%y= 2*positions_y - particles_y
+          particles(n)%speed_y = - particles(n)%speed_y
+          
+         IF ( use_sgs_for_particles )  THEN
+           particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs
+         ENDIF  
+
+          GOTO 999
+       ENDIF
+
+
+      ENDIF 
+    ENDDO
+   ENDIF
+
+     999 CONTINUE
+    
+   ENDIF
+
+!---------------------------------------------------
+
+       ENDDO   ! particle reflection from walls
+
+       CALL cpu_log( log_point_s(48), 'advec_part_refle', 'stop' )
+
+!
+!--    User-defined actions after the evaluation of the new particle position
+       CALL user_advec_particles
 
 !
@@ -1500,8 +2748,8 @@
           trlp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
                                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
-                                0, 0, 0, 0 )
+                                0.0, 0, 0, 0, 0 )
           trrp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
                                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
-                                0, 0, 0, 0 )
+                                0.0, 0, 0, 0, 0 )
 
           IF ( use_particle_tails )  THEN
@@ -1933,8 +3181,8 @@
           trsp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
                                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
-                                0, 0, 0, 0 )
+                                0.0, 0, 0, 0, 0 )
           trnp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
                                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
-                                0, 0, 0, 0 )
+                                0.0, 0, 0, 0, 0 )
 
           IF ( use_particle_tails )  THEN
@@ -2460,4 +3708,20 @@
 
           nn = particles(n)%tail_id
+
+!
+!--       Stop if particles have moved further than the length of one 
+!--       PE subdomain
+          IF ( ABS(particles(n)%speed_x) >                                  & 
+               ((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m)  .OR. &
+               ABS(particles(n)%speed_y) >                                  &
+               ((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) )  THEN
+
+             PRINT*, '+++ advec_particles: particle too fast.  n = ', n
+#if defined( __parallel )
+             CALL MPI_ABORT( comm2d, 9999, ierr )
+#else
+             CALL local_stop
+#endif
+          ENDIF
 
           IF ( particles(n)%age > particle_maximum_age  .AND.  &
Index: palm/trunk/SOURCE/buoyancy.f90
===================================================================
--- palm/trunk/SOURCE/buoyancy.f90	(revision 56)
+++ palm/trunk/SOURCE/buoyancy.f90	(revision 57)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Reference temperature pt_reference can be used.
 !
 ! Former revisions:
@@ -59,14 +59,27 @@
 !
 !--       Normal case: horizontal surface
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb_s_inner(j,i)+1, nzt-1
-                    tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                    &
+          IF ( use_pt_reference )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                      tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                 &
+                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / pt_reference + &
+                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / pt_reference   &
+                                                            )
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                      tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                  &
                         ( theta(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
                         ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
-                                                          )
-                ENDDO
-             ENDDO
-          ENDDO
+                                                            )
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
 
        ELSE
@@ -136,10 +149,19 @@
 !
 !--       Normal case: horizontal surface
-          DO  k = nzb_s_inner(j,i)+1, nzt-1
-              tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                          &
+          IF ( use_pt_reference )  THEN
+             DO  k = nzb_s_inner(j,i)+1, nzt-1
+                 tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                      &
+                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / pt_reference + &
+                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / pt_reference   &
+                                                       )
+             ENDDO
+          ELSE
+             DO  k = nzb_s_inner(j,i)+1, nzt-1
+                 tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                       &
                         ( theta(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
                         ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
-                                                    )
-          ENDDO
+                                                       )
+             ENDDO
+          ENDIF
 
        ELSE
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 56)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 57)
@@ -5,5 +5,6 @@
 ! -----------------
 ! "by_user" allowed as initializing action, -data_output_ts,
-! leapfrog with non-flat topography not allowed any more
+! leapfrog with non-flat topography not allowed any more, pt_reference is
+! checked
 !
 ! Former revisions:
@@ -574,4 +575,8 @@
 
 !
+!-- Reference temperature to be used in buoyancy terms
+    IF ( pt_reference /= 9999999.9 )  use_pt_reference = .TRUE.
+
+!
 !-- In case of a given slope, compute the relevant quantities
     IF ( alpha_surface /= 0.0 )  THEN
Index: palm/trunk/SOURCE/diffusion_e.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_e.f90	(revision 56)
+++ palm/trunk/SOURCE/diffusion_e.f90	(revision 57)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Reference temperature pt_reference can be used in buoyancy term
 !
 ! Former revisions:
@@ -82,6 +82,11 @@
                 dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
                 IF ( dpt_dz > 0.0 ) THEN
-                   l_stable = 0.76 * SQRT( e(k,j,i) ) / &
-                                     SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
+                   IF ( use_pt_reference )  THEN
+                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
+                                        SQRT( g / pt_reference * dpt_dz ) + 1E-5
+                   ELSE
+                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
+                                        SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
+                   ENDIF
                 ELSE
                    l_stable = l_grid(k)
@@ -195,6 +200,11 @@
           dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
           IF ( dpt_dz > 0.0 ) THEN
-             l_stable = 0.76 * SQRT( e(k,j,i) ) / &
-                               SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
+             IF ( use_pt_reference )  THEN
+                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
+                                  SQRT( g / pt_reference * dpt_dz ) + 1E-5
+             ELSE
+                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
+                                  SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
+             ENDIF
           ELSE
              l_stable = l_grid(k)
Index: palm/trunk/SOURCE/diffusion_u.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_u.f90	(revision 56)
+++ palm/trunk/SOURCE/diffusion_u.f90	(revision 57)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! Wall functions now include diabatic conditions, call of routine wall_fluxes
+! Wall functions now include diabatic conditions, call of routine wall_fluxes,
+! z0 removed from argument list
 ! 
 ! Former revisions:
@@ -50,5 +51,5 @@
 ! Call for all grid points
 !------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w, z0 )
+    SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w )
 
        USE control_parameters
@@ -61,5 +62,4 @@
        REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
-       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
        REAL, DIMENSION(:,:),   POINTER ::  usws
@@ -207,5 +207,5 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
-                               v, w, z0 )
+                               v, w )
 
        USE control_parameters
@@ -218,5 +218,4 @@
        REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
-       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
        REAL, DIMENSION(nzb:nzt+1)      ::  usvs
Index: palm/trunk/SOURCE/diffusion_v.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_v.f90	(revision 56)
+++ palm/trunk/SOURCE/diffusion_v.f90	(revision 57)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! Wall functions now include diabatic conditions, call of routine wall_fluxes
+! Wall functions now include diabatic conditions, call of routine wall_fluxes,
+! z0 removed from argument list
 !
 ! Former revisions:
@@ -48,5 +49,5 @@
 ! Call for all grid points
 !------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w, z0 )
+    SUBROUTINE diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w )
 
        USE control_parameters
@@ -59,5 +60,4 @@
        REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1)
-       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
        REAL, DIMENSION(:,:),   POINTER ::  vsws
@@ -205,5 +205,5 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE diffusion_v_ij( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
-                               vsws, w, z0 )
+                               vsws, w )
 
        USE control_parameters
@@ -216,5 +216,4 @@
        REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1)
-       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
        REAL, DIMENSION(nzb:nzt+1)      ::  vsus
Index: palm/trunk/SOURCE/diffusion_w.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_w.f90	(revision 56)
+++ palm/trunk/SOURCE/diffusion_w.f90	(revision 57)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! Wall functions now include diabatic conditions, call of routine wall_fluxes
+! Wall functions now include diabatic conditions, call of routine wall_fluxes,
+! z0 removed from argument list
 !
 ! Former revisions:
@@ -46,5 +47,5 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, &
-                            w, z0 )
+                            w )
 
        USE control_parameters
@@ -59,5 +60,4 @@
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1), &
                    km_damp_y(nys-1:nyn+1)
-       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
@@ -194,5 +194,5 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
-                               tend, u, v, w, z0 )
+                               tend, u, v, w )
 
        USE control_parameters
@@ -207,5 +207,4 @@
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1), &
                    km_damp_y(nys-1:nyn+1)
-       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
        REAL, DIMENSION(nzb:nzt+1)      ::  wsus, wsvs
Index: palm/trunk/SOURCE/diffusivities.f90
===================================================================
--- palm/trunk/SOURCE/diffusivities.f90	(revision 56)
+++ palm/trunk/SOURCE/diffusivities.f90	(revision 57)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Reference temperature pt_reference can be used in buoyancy term
 ! 
 ! Former revisions:
@@ -90,6 +90,11 @@
              dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
              IF ( dpt_dz > 0.0 ) THEN
-                l_stable = 0.76 * sqrt_e(k) / &
-                                  SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
+                IF ( use_pt_reference )  THEN
+                   l_stable = 0.76 * sqrt_e(k) / &
+                                     SQRT( g / pt_reference * dpt_dz ) + 1E-5
+                ELSE
+                   l_stable = 0.76 * sqrt_e(k) / &
+                                     SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
+                ENDIF
              ELSE
                 l_stable = l_grid(k)
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 56)
+++ palm/trunk/SOURCE/header.f90	(revision 57)
@@ -861,10 +861,11 @@
 !-- Other quantities
     WRITE ( io, 411 )  g
+    IF ( use_pt_reference )  WRITE ( io, 412 )  pt_reference
 
 !
 !-- Cloud physics parameters
     IF ( cloud_physics ) THEN
-       WRITE ( io, 412 )
-       WRITE ( io, 413 ) surface_pressure, r_d, rho_surface, cp, l_v
+       WRITE ( io, 415 )
+       WRITE ( io, 416 ) surface_pressure, r_d, rho_surface, cp, l_v
     ENDIF
 
@@ -1326,7 +1327,8 @@
             '                            f*    = ',F9.6,' 1/s')
 411 FORMAT (/'    Gravity             :   g     = ',F4.1,' m/s**2')
-412 FORMAT (/'    Cloud physics parameters:'/ &
+412 FORMAT (/'    Reference temperature in buoyancy terms: ',F8.4,' K')
+415 FORMAT (/'    Cloud physics parameters:'/ &
              '    ------------------------'/)
-413 FORMAT ('        Surface pressure   :   p_0   = ',F7.2,' hPa'/      &
+416 FORMAT ('        Surface pressure   :   p_0   = ',F7.2,' hPa'/      &
             '        Gas constant       :   R     = ',F5.1,' J/(kg K)'/ &
             '        Density of air     :   rho_0 = ',F5.3,' kg/m**3'/  &
Index: palm/trunk/SOURCE/init_particles.f90
===================================================================
--- palm/trunk/SOURCE/init_particles.f90	(revision 56)
+++ palm/trunk/SOURCE/init_particles.f90	(revision 57)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! displacements for mpi_particle_type changed, age_m initialized
 !
 ! Former revisions:
@@ -61,7 +61,7 @@
     blocklengths(1)  = 18; blocklengths(2)  = 4; blocklengths(3)  =  1
 #if defined( __t3eb )
-    displacements(1) = 0; displacements(2) = 144; displacements(3) = 176
+    displacements(1) = 0; displacements(2) = 152; displacements(3) = 184
 #else
-    displacements(1) = 0; displacements(2) = 144; displacements(3) = 160
+    displacements(1) = 0; displacements(2) = 152; displacements(3) = 168
 #endif
     types(1) = MPI_REAL
@@ -192,5 +192,5 @@
        particles = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
                                   0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
-                                  0, 0, 0, 0 )
+                                  0.0, 0, 0, 0, 0 )
        particle_groups = particle_groups_type( 0.0, 0.0, 0.0, 0.0 )
 
@@ -267,4 +267,5 @@
                             particles(n)%z             = pos_z
                             particles(n)%age           = 0.0
+                            particles(n)%age_m         = 0.0
                             particles(n)%dt_sum        = 0.0
                             particles(n)%dvrp_psize    = dvrp_psize
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 56)
+++ palm/trunk/SOURCE/modules.f90	(revision 57)
@@ -7,5 +7,6 @@
 ! +array rif_wall
 ! +netcdf_64bit_3d, zu_s_inner, zw_w_inner, id_var_zusi_*, id_var_zwwi_*,
-! ts_value, u_nzb_p1_for_vfc, v_nzb_p1_for_vfc
+! ts_value, u_nzb_p1_for_vfc, v_nzb_p1_for_vfc, pt_reference, use_pt_reference,
+! +age_m in particle_type
 ! -data_output_ts, dots_n
 ! arrays dots_label and dots_unit now dimensioned with dots_max
@@ -296,6 +297,6 @@
                 sloping_surface = .FALSE., stop_dt = .FALSE., &
                 terminate_run = .FALSE., use_prior_plot1d_parameters = .FALSE.,&
-                use_surface_fluxes = .FALSE., use_top_fluxes = .FALSE., &
-                use_ug_for_galilei_tr = .TRUE., &
+                use_pt_reference = .FALSE., use_surface_fluxes = .FALSE., &
+                use_top_fluxes = .FALSE., use_ug_for_galilei_tr = .TRUE., &
                 use_upstream_for_tke = .FALSE., wall_adjustment = .TRUE.
 
@@ -330,10 +331,11 @@
              overshoot_limit_u = 0.0, overshoot_limit_v = 0.0, &
              overshoot_limit_w = 0.0, particle_maximum_age = 9999999.9, &
-             phi = 55.0, prandtl_number = 1.0, pt_slope_offset = 0.0, &
-             pt_surface = 300.0, pt_surface_initial_change = 0.0, &
-             q_surface = 0.0, q_surface_initial_change = 0.0, &
-             rayleigh_damping_factor = -1.0, rayleigh_damping_height = -1.0, &
-             residual_limit = 1.0E-4, restart_time = 9999999.9, rif_max = 1.0, &
-             rif_min = -5.0, roughness_length = 0.1, simulated_time = 0.0, &
+             phi = 55.0, prandtl_number = 1.0, pt_reference = 9999999.9, &
+             pt_slope_offset = 0.0, pt_surface = 300.0, &
+             pt_surface_initial_change = 0.0, q_surface = 0.0, &
+             q_surface_initial_change = 0.0, rayleigh_damping_factor = -1.0, &
+             rayleigh_damping_height = -1.0, residual_limit = 1.0E-4, &
+             restart_time = 9999999.9, rif_max = 1.0, rif_min = -5.0, &
+             roughness_length = 0.1, simulated_time = 0.0, &
              simulated_time_at_begin, sin_alpha_surface, &
              skip_time_data_output = 0.0, skip_time_data_output_av = 9999999.9,&
@@ -828,7 +830,7 @@
     TYPE particle_type
        SEQUENCE
-       REAL    ::  age, dt_sum, dvrp_psize, e_m, origin_x, origin_y, origin_z, &
-                   radius, speed_x, speed_x_sgs, speed_y, speed_y_sgs,         &
-                   speed_z, speed_z_sgs, weight_factor, x, y, z
+       REAL    ::  age, age_m, dt_sum, dvrp_psize, e_m, origin_x, origin_y, &
+                   origin_z, radius, speed_x, speed_x_sgs, speed_y,         &
+                   speed_y_sgs, speed_z, speed_z_sgs, weight_factor, x, y, z
        INTEGER ::  color, group, tailpoints, tail_id
     END TYPE particle_type
Index: palm/trunk/SOURCE/palm.f90
===================================================================
--- palm/trunk/SOURCE/palm.f90	(revision 56)
+++ palm/trunk/SOURCE/palm.f90	(revision 57)
@@ -1,12 +1,8 @@
-#if defined( __vtk )
- SUBROUTINE palm
-#else
  PROGRAM palm
-#endif
 
 !------------------------------------------------------------------------------!
 ! Actual revisions:
 ! -----------------
-! TEST: open(9)
+! __vtk directives removed, open unit 9 for debug output
 !
 ! Former revisions:
@@ -99,4 +95,8 @@
 
 !
+!-- Open a file for debug output
+    OPEN( 9, FILE='DEBUG'//myid_char, FORM='FORMATTED' )
+
+!
 !-- Generate grid parameters
     CALL init_grid
@@ -106,5 +106,4 @@
     CALL check_parameters
 
-    OPEN( 9, FILE='DEBUG'//myid_char, FORM='FORMATTED' )
 !
 !-- Initialize all necessary variables
@@ -169,9 +168,4 @@
 #endif
 
-#if defined( __vtk )
- END SUBROUTINE palm
-#else
  END PROGRAM palm
-#endif
 
-
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 56)
+++ palm/trunk/SOURCE/parin.f90	(revision 57)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! +netcdf_64bit_3d in d3par, -data_output_ts
+! +netcdf_64bit_3d in d3par, +pt_reference in inipar, -data_output_ts
 !
 ! Former revisions:
@@ -65,15 +65,16 @@
                        overshoot_limit_pt, overshoot_limit_u, &
                        overshoot_limit_v, overshoot_limit_w, passive_scalar, &
-                       phi, prandtl_layer, precipitation, pt_surface, &
-                       pt_surface_initial_change, pt_vertical_gradient, &
-                       pt_vertical_gradient_level, q_surface, &
-                       q_surface_initial_change, q_vertical_gradient, &
-                       q_vertical_gradient_level, radiation, random_generator, &
-                       random_heatflux, rif_max, rif_min, roughness_length, &
-                       scalar_advec, statistic_regions, surface_heatflux, &
-                       surface_pressure, surface_scalarflux, surface_waterflux,&
-                       s_surface, s_surface_initial_change, &
-                       s_vertical_gradient, s_vertical_gradient_level, &
-                       top_heatflux, timestep_scheme, topography, ug_surface, &
+                       phi, prandtl_layer, precipitation, pt_reference, &
+                       pt_surface, pt_surface_initial_change, &
+                       pt_vertical_gradient, pt_vertical_gradient_level, &
+                       q_surface, q_surface_initial_change, &
+                       q_vertical_gradient, q_vertical_gradient_level, &
+                       radiation, random_generator, random_heatflux, rif_max, &
+                       rif_min, roughness_length, scalar_advec, &
+                       statistic_regions, surface_heatflux, surface_pressure, &
+                       surface_scalarflux, surface_waterflux, s_surface, &
+                       s_surface_initial_change, s_vertical_gradient, &
+                       s_vertical_gradient_level, top_heatflux, &
+                       timestep_scheme, topography, ug_surface, &
                        ug_vertical_gradient, ug_vertical_gradient_level, &
                        ups_limit_e, ups_limit_pt, ups_limit_u, ups_limit_v, &
Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 56)
+++ palm/trunk/SOURCE/production_e.f90	(revision 57)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! Wall functions now include diabatic conditions, call of routine wall_fluxes_e
+! Wall functions now include diabatic conditions, call of routine wall_fluxes_e,
+! reference temperature pt_reference can be used in buoyancy term
 !
 ! Former revisions:
@@ -359,22 +360,43 @@
           IF ( .NOT. moisture )  THEN
 
-             DO  j = nys, nyn
-
-                DO  k = nzb_diff_s_inner(j,i), nzt_diff
-                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
-                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
+             IF ( use_pt_reference )  THEN
+
+                DO  j = nys, nyn
+                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
+                      tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g/pt_reference * &
+                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
+                   ENDDO
+
+                   IF ( use_surface_fluxes )  THEN
+                      k = nzb_diff_s_inner(j,i)-1
+                      tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
+                   ENDIF
+
+                   IF ( use_top_fluxes )  THEN
+                      k = nzt
+                      tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
+                   ENDIF
                 ENDDO
 
-                IF ( use_surface_fluxes )  THEN
-                   k = nzb_diff_s_inner(j,i)-1
-                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
-                ENDIF
-
-                IF ( use_top_fluxes )  THEN
-                   k = nzt
-                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
-                ENDIF
-
-             ENDDO
+             ELSE
+
+                DO  j = nys, nyn
+                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
+                      tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
+                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
+                   ENDDO
+
+                   IF ( use_surface_fluxes )  THEN
+                      k = nzb_diff_s_inner(j,i)-1
+                      tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
+                   ENDIF
+
+                   IF ( use_top_fluxes )  THEN
+                      k = nzt
+                      tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
+                   ENDIF
+                ENDDO
+
+             ENDIF
 
           ELSE
@@ -736,18 +758,39 @@
        IF ( .NOT. moisture )  THEN
 
-          DO  k = nzb_diff_s_inner(j,i), nzt_diff
-             tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
-                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
-          ENDDO
-
-          IF ( use_surface_fluxes )  THEN
-             k = nzb_diff_s_inner(j,i)-1
-             tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
-          ENDIF
-
-          IF ( use_top_fluxes )  THEN
-             k = nzt
-             tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
-          ENDIF
+          IF ( use_pt_reference )  THEN
+
+             DO  k = nzb_diff_s_inner(j,i), nzt_diff
+                tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt_reference * &
+                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
+             ENDDO
+
+             IF ( use_surface_fluxes )  THEN
+                k = nzb_diff_s_inner(j,i)-1
+                tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
+             ENDIF
+
+             IF ( use_top_fluxes )  THEN
+                k = nzt
+                tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
+             ENDIF
+
+          ELSE
+
+             DO  k = nzb_diff_s_inner(j,i), nzt_diff
+                tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
+                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
+             ENDDO
+
+             IF ( use_surface_fluxes )  THEN
+                k = nzb_diff_s_inner(j,i)-1
+                tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
+             ENDIF
+
+             IF ( use_top_fluxes )  THEN
+                k = nzt
+                tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
+             ENDIF
+
+         ENDIF
 
        ELSE
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 56)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 57)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! z0 removed from arguments in calls of diffusion_u/v/w
 !
 ! Former revisions:
@@ -130,8 +130,8 @@
           IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
              CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, u_m,   &
-                               usws_m, v_m, w_m, z0 )
+                               usws_m, v_m, w_m )
           ELSE
              CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
-                               v, w, z0 )
+                               v, w )
           ENDIF
           CALL coriolis( i, j, 1 )
@@ -198,8 +198,8 @@
           IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
              CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, u_m, &
-                               v_m, vsws_m, w_m, z0 )
+                               v_m, vsws_m, w_m )
           ELSE
              CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v,  &
-                               vsws, w, z0 )
+                               vsws, w )
           ENDIF
           CALL coriolis( i, j, 2 )
@@ -265,8 +265,8 @@
           IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
              CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, km_damp_y,  &
-                               tend, u_m, v_m, w_m, z0 )
+                               tend, u_m, v_m, w_m )
           ELSE
              CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y,    &
-                               tend, u, v, w, z0 )
+                               tend, u, v, w )
           ENDIF
           CALL coriolis( i, j, 3 )
@@ -668,8 +668,8 @@
              THEN
                 CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend,  &
-                                  u_m, usws_m, v_m, w_m, z0 )
+                                  u_m, usws_m, v_m, w_m )
              ELSE
                 CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, &
-                                  usws, v, w, z0 )
+                                  usws, v, w )
              ENDIF
              CALL coriolis( i, j, 1 )
@@ -718,8 +718,8 @@
              THEN
                 CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend,     &
-                                  u_m, v_m, vsws_m, w_m, z0 )
+                                  u_m, v_m, vsws_m, w_m )
              ELSE
                 CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
-                                  vsws, w, z0 )
+                                  vsws, w )
              ENDIF
              CALL coriolis( i, j, 2 )
@@ -767,8 +767,8 @@
              THEN
                 CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x,          &
-                                  km_damp_y, tend, u_m, v_m, w_m, z0 )
+                                  km_damp_y, tend, u_m, v_m, w_m )
              ELSE
                 CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
-                                  tend, u, v, w, z0 )
+                                  tend, u, v, w )
              ENDIF
              CALL coriolis( i, j, 3 )
@@ -1041,7 +1041,7 @@
     IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
        CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, v_m, &
-                         w_m, z0 )
+                         w_m )
     ELSE
-       CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w, z0 )
+       CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w )
     ENDIF
     CALL coriolis( 1 )
@@ -1113,7 +1113,7 @@
     IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
        CALL diffusion_v( ddzu, ddzw, km_m, km_damp_x, tend, u_m, v_m, vsws_m, &
-                         w_m, z0 )
+                         w_m )
     ELSE
-       CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w, z0 )
+       CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w )
     ENDIF
     CALL coriolis( 2 )
@@ -1184,8 +1184,7 @@
     IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
        CALL diffusion_w( ddzu, ddzw, km_m, km_damp_x, km_damp_y, tend, u_m, &
-                         v_m, w_m, z0 )
+                         v_m, w_m )
     ELSE
-       CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w, &
-                         z0 )
+       CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w )
     ENDIF
     CALL coriolis( 3 )
Index: palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- palm/trunk/SOURCE/read_var_list.f90	(revision 56)
+++ palm/trunk/SOURCE/read_var_list.f90	(revision 57)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +pt_reference
 !
 ! Former revisions:
@@ -252,4 +252,6 @@
           CASE ( 'pt_init' )
              READ ( 13 )  pt_init
+          CASE ( 'pt_reference' )
+             READ ( 13 )  pt_reference
           CASE ( 'pt_surface' )
              READ ( 13 )  pt_surface
Index: palm/trunk/SOURCE/user_interface.f90
===================================================================
--- palm/trunk/SOURCE/user_interface.f90	(revision 56)
+++ palm/trunk/SOURCE/user_interface.f90	(revision 57)
@@ -8,4 +8,5 @@
 ! routine user_statistics now has one argument (sr),
 ! sample for generating time series quantities added
+! Bugfix in sample for reading user defined data from restart file (user_init)
 !
 ! Former revisions:
@@ -174,4 +175,7 @@
 !
 !          END SELECT
+!
+!          READ ( 13 )  field_chr
+!
 !       ENDDO
 !    ENDIF
Index: palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- palm/trunk/SOURCE/write_var_list.f90	(revision 56)
+++ palm/trunk/SOURCE/write_var_list.f90	(revision 57)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +pt_refrence
 !
 ! Former revisions:
@@ -218,4 +218,6 @@
     WRITE ( 14 )  'pt_init                       '
     WRITE ( 14 )  pt_init
+    WRITE ( 14 )  'pt_reference                  '
+    WRITE ( 14 )  pt_reference
     WRITE ( 14 )  'pt_surface                    '
     WRITE ( 14 )  pt_surface
