source: palm/trunk/SOURCE/lpm_advec.f90 @ 2318

Last change on this file since 2318 was 2318, checked in by suehring, 4 years ago

get topograpyhy top index via function call

  • Property svn:keywords set to Id
File size: 77.0 KB
Line 
1!> @file lpm_advec.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: lpm_advec.f90 2318 2017-07-20 17:27:44Z suehring $
27! Get topography top index via Function call
28!
29! 2317 2017-07-20 17:27:19Z suehring
30!
31! 2232 2017-05-30 17:47:52Z suehring
32! Adjustments to new topography and surface concept
33!
34! 2100 2017-01-05 16:40:16Z suehring
35! Prevent extremely large SGS-velocities in regions where TKE is zero, e.g.
36! at the begin of simulations and/or in non-turbulent regions.
37!
38! 2000 2016-08-20 18:09:15Z knoop
39! Forced header and separation lines into 80 columns
40!
41! 1936 2016-06-13 13:37:44Z suehring
42! Formatting adjustments
43!
44! 1929 2016-06-09 16:25:25Z suehring
45! Put stochastic equation in an extra subroutine.
46! Set flag for stochastic equation to communicate whether a particle is near
47! topography. This case, memory and drift term are disabled in the Weil equation.
48!
49! Enable vertical logarithmic interpolation also above topography. This case,
50! set a lower limit for the friction velocity, as it can become very small
51! in narrow street canyons, leading to too large particle speeds.
52!
53! 1888 2016-04-21 12:20:49Z suehring
54! Bugfix concerning logarithmic interpolation of particle speed
55!
56! 1822 2016-04-07 07:49:42Z hoffmann
57! Random velocity fluctuations for particles added. Terminal fall velocity
58! for droplets is calculated from a parameterization (which is better than
59! the previous, physically correct calculation, which demands a very short
60! time step that is not used in the model).
61!
62! Unused variables deleted.
63!
64! 1691 2015-10-26 16:17:44Z maronga
65! Renamed prandtl_layer to constant_flux_layer.
66!
67! 1685 2015-10-08 07:32:13Z raasch
68! TKE check for negative values (so far, only zero value was checked)
69! offset_ocean_nzt_m1 removed
70!
71! 1682 2015-10-07 23:56:08Z knoop
72! Code annotations made doxygen readable
73!
74! 1583 2015-04-15 12:16:27Z suehring
75! Bugfix: particle advection within Prandtl-layer in case of Galilei
76! transformation.
77!
78! 1369 2014-04-24 05:57:38Z raasch
79! usage of module interfaces removed
80!
81! 1359 2014-04-11 17:15:14Z hoffmann
82! New particle structure integrated.
83! Kind definition added to all floating point numbers.
84!
85! 1322 2014-03-20 16:38:49Z raasch
86! REAL constants defined as wp_kind
87!
88! 1320 2014-03-20 08:40:49Z raasch
89! ONLY-attribute added to USE-statements,
90! kind-parameters added to all INTEGER and REAL declaration statements,
91! kinds are defined in new module kinds,
92! revision history before 2012 removed,
93! comment fields (!:) to be used for variable explanations added to
94! all variable declaration statements
95!
96! 1314 2014-03-14 18:25:17Z suehring
97! Vertical logarithmic interpolation of horizontal particle speed for particles
98! between roughness height and first vertical grid level.
99!
100! 1036 2012-10-22 13:43:42Z raasch
101! code put under GPL (PALM 3.9)
102!
103! 849 2012-03-15 10:35:09Z raasch
104! initial revision (former part of advec_particles)
105!
106!
107! Description:
108! ------------
109!> Calculation of new particle positions due to advection using a simple Euler
110!> scheme. Particles may feel inertia effects. SGS transport can be included
111!> using the stochastic model of Weil et al. (2004, JAS, 61, 2877-2887).
112!------------------------------------------------------------------------------!
113 SUBROUTINE lpm_advec (ip,jp,kp)
114 
115
116    USE arrays_3d,                                                             &
117        ONLY:  de_dx, de_dy, de_dz, diss, e, km, u, v, w, zu, zw
118
119    USE cpulog
120
121    USE pegrid
122
123    USE control_parameters,                                                    &
124        ONLY:  atmos_ocean_sign, cloud_droplets, constant_flux_layer, dt_3d,   &
125               dt_3d_reached_l, dz, g, kappa, topography, u_gtrans, v_gtrans
126
127    USE grid_variables,                                                        &
128        ONLY:  ddx, dx, ddy, dy
129       
130    USE indices,                                                               &
131        ONLY:  nzb, nzt
132       
133    USE kinds
134   
135    USE particle_attributes,                                                   &
136        ONLY:  block_offset, c_0, dt_min_part, grid_particles,                 &
137               iran_part, log_z_z0, number_of_particles, number_of_sublayers,  &
138               particles, particle_groups, offset_ocean_nzt, sgs_wf_part,      &
139               use_sgs_for_particles, vertical_particle_advection, z0_av_global
140       
141    USE statistics,                                                            &
142        ONLY:  hom
143
144    USE surface_mod,                                                           &
145        ONLY:  get_topography_top_index, surf_def_h, surf_lsm_h, surf_usm_h
146
147    IMPLICIT NONE
148
149    INTEGER(iwp) ::  agp                         !< loop variable
150    INTEGER(iwp) ::  gp_outside_of_building(1:8) !< number of grid points used for particle interpolation in case of topography
151    INTEGER(iwp) ::  i                           !< index variable along x
152    INTEGER(iwp) ::  ip                          !< index variable along x
153    INTEGER(iwp) ::  ilog                        !< index variable along x
154    INTEGER(iwp) ::  j                           !< index variable along y
155    INTEGER(iwp) ::  jp                          !< index variable along y
156    INTEGER(iwp) ::  jlog                        !< index variable along y
157    INTEGER(iwp) ::  k                           !< index variable along z
158    INTEGER(iwp) ::  k_wall                      !< vertical index of topography top
159    INTEGER(iwp) ::  kp                          !< index variable along z
160    INTEGER(iwp) ::  kw                          !< index variable along z
161    INTEGER(iwp) ::  n                           !< loop variable over all particles in a grid box
162    INTEGER(iwp) ::  nb                          !< block number particles are sorted in
163    INTEGER(iwp) ::  num_gp                      !< number of adjacent grid points inside topography
164    INTEGER(iwp) ::  surf_start                  !< Index on surface data-type for current grid box
165
166    INTEGER(iwp), DIMENSION(0:7) ::  start_index !< start particle index for current block
167    INTEGER(iwp), DIMENSION(0:7) ::  end_index   !< start particle index for current block
168
169    REAL(wp) ::  aa                 !< dummy argument for horizontal particle interpolation
170    REAL(wp) ::  bb                 !< dummy argument for horizontal particle interpolation
171    REAL(wp) ::  cc                 !< dummy argument for horizontal particle interpolation
172    REAL(wp) ::  d_sum              !< dummy argument for horizontal particle interpolation in case of topography
173    REAL(wp) ::  d_z_p_z0           !< inverse of interpolation length for logarithmic interpolation
174    REAL(wp) ::  dd                 !< dummy argument for horizontal particle interpolation
175    REAL(wp) ::  de_dx_int_l        !< x/y-interpolated TKE gradient (x) at particle position at lower vertical level
176    REAL(wp) ::  de_dx_int_u        !< x/y-interpolated TKE gradient (x) at particle position at upper vertical level
177    REAL(wp) ::  de_dy_int_l        !< x/y-interpolated TKE gradient (y) at particle position at lower vertical level
178    REAL(wp) ::  de_dy_int_u        !< x/y-interpolated TKE gradient (y) at particle position at upper vertical level
179    REAL(wp) ::  de_dt              !< temporal derivative of TKE experienced by the particle
180    REAL(wp) ::  de_dt_min          !< lower level for temporal TKE derivative
181    REAL(wp) ::  de_dz_int_l        !< x/y-interpolated TKE gradient (z) at particle position at lower vertical level
182    REAL(wp) ::  de_dz_int_u        !< x/y-interpolated TKE gradient (z) at particle position at upper vertical level
183    REAL(wp) ::  diameter           !< diamter of droplet
184    REAL(wp) ::  diss_int_l         !< x/y-interpolated dissipation at particle position at lower vertical level
185    REAL(wp) ::  diss_int_u         !< x/y-interpolated dissipation at particle position at upper vertical level
186    REAL(wp) ::  dt_gap             !< remaining time until particle time integration reaches LES time
187    REAL(wp) ::  dt_particle_m      !< previous particle time step
188    REAL(wp) ::  e_int_l            !< x/y-interpolated TKE at particle position at lower vertical level
189    REAL(wp) ::  e_int_u            !< x/y-interpolated TKE at particle position at upper vertical level
190    REAL(wp) ::  e_mean_int         !< horizontal mean TKE at particle height
191    REAL(wp) ::  exp_arg            !<
192    REAL(wp) ::  exp_term           !<
193    REAL(wp) ::  gg                 !< dummy argument for horizontal particle interpolation
194    REAL(wp) ::  height_p           !< dummy argument for logarithmic interpolation
195    REAL(wp) ::  lagr_timescale     !< Lagrangian timescale
196    REAL(wp) ::  location(1:30,1:3) !< wall locations
197    REAL(wp) ::  log_z_z0_int       !< logarithmus used for surface_layer interpolation
198    REAL(wp) ::  random_gauss       !<
199    REAL(wp) ::  RL                 !< Lagrangian autocorrelation coefficient
200    REAL(wp) ::  rg1                !< Gaussian distributed random number
201    REAL(wp) ::  rg2                !< Gaussian distributed random number
202    REAL(wp) ::  rg3                !< Gaussian distributed random number
203    REAL(wp) ::  sigma              !< velocity standard deviation
204    REAL(wp) ::  u_int_l            !< x/y-interpolated u-component at particle position at lower vertical level
205    REAL(wp) ::  u_int_u            !< x/y-interpolated u-component at particle position at upper vertical level
206    REAL(wp) ::  us_int             !< friction velocity at particle grid box
207    REAL(wp) ::  usws_int           !< surface momentum flux (u component) at particle grid box
208    REAL(wp) ::  v_int_l            !< x/y-interpolated v-component at particle position at lower vertical level
209    REAL(wp) ::  v_int_u            !< x/y-interpolated v-component at particle position at upper vertical level
210    REAL(wp) ::  vsws_int           !< surface momentum flux (u component) at particle grid box
211    REAL(wp) ::  vv_int             !<
212    REAL(wp) ::  w_int_l            !< x/y-interpolated w-component at particle position at lower vertical level
213    REAL(wp) ::  w_int_u            !< x/y-interpolated w-component at particle position at upper vertical level
214    REAL(wp) ::  w_s                !< terminal velocity of droplets
215    REAL(wp) ::  x                  !< dummy argument for horizontal particle interpolation
216    REAL(wp) ::  y                  !< dummy argument for horizontal particle interpolation
217    REAL(wp) ::  z_p                !< surface layer height (0.5 dz)
218
219    REAL(wp), PARAMETER ::  a_rog = 9.65_wp      !< parameter for fall velocity
220    REAL(wp), PARAMETER ::  b_rog = 10.43_wp     !< parameter for fall velocity
221    REAL(wp), PARAMETER ::  c_rog = 0.6_wp       !< parameter for fall velocity
222    REAL(wp), PARAMETER ::  k_cap_rog = 4.0_wp   !< parameter for fall velocity
223    REAL(wp), PARAMETER ::  k_low_rog = 12.0_wp  !< parameter for fall velocity
224    REAL(wp), PARAMETER ::  d0_rog = 0.745_wp    !< separation diameter
225
226    REAL(wp), DIMENSION(1:30) ::  d_gp_pl !< dummy argument for particle interpolation scheme in case of topography
227    REAL(wp), DIMENSION(1:30) ::  de_dxi  !< horizontal TKE gradient along x at adjacent wall
228    REAL(wp), DIMENSION(1:30) ::  de_dyi  !< horizontal TKE gradient along y at adjacent wall
229    REAL(wp), DIMENSION(1:30) ::  de_dzi  !< horizontal TKE gradient along z at adjacent wall
230    REAL(wp), DIMENSION(1:30) ::  dissi   !< dissipation at adjacent wall
231    REAL(wp), DIMENSION(1:30) ::  ei      !< TKE at adjacent wall
232
233    REAL(wp), DIMENSION(number_of_particles) ::  term_1_2     !< flag to communicate whether a particle is near topography or not
234    REAL(wp), DIMENSION(number_of_particles) ::  dens_ratio   !<
235    REAL(wp), DIMENSION(number_of_particles) ::  de_dx_int    !< horizontal TKE gradient along x at particle position
236    REAL(wp), DIMENSION(number_of_particles) ::  de_dy_int    !< horizontal TKE gradient along y at particle position
237    REAL(wp), DIMENSION(number_of_particles) ::  de_dz_int    !< horizontal TKE gradient along z at particle position
238    REAL(wp), DIMENSION(number_of_particles) ::  diss_int     !< dissipation at particle position
239    REAL(wp), DIMENSION(number_of_particles) ::  dt_particle  !< particle time step
240    REAL(wp), DIMENSION(number_of_particles) ::  e_int        !< TKE at particle position
241    REAL(wp), DIMENSION(number_of_particles) ::  fs_int       !< weighting factor for subgrid-scale particle speed
242    REAL(wp), DIMENSION(number_of_particles) ::  u_int        !< u-component of particle speed
243    REAL(wp), DIMENSION(number_of_particles) ::  v_int        !< v-component of particle speed
244    REAL(wp), DIMENSION(number_of_particles) ::  w_int        !< w-component of particle speed
245    REAL(wp), DIMENSION(number_of_particles) ::  xv           !< x-position
246    REAL(wp), DIMENSION(number_of_particles) ::  yv           !< y-position
247    REAL(wp), DIMENSION(number_of_particles) ::  zv           !< z-position
248
249    REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !< vector of Gaussian distributed random numbers
250
251    CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' )
252
253!
254!-- Determine height of Prandtl layer and distance between Prandtl-layer
255!-- height and horizontal mean roughness height, which are required for
256!-- vertical logarithmic interpolation of horizontal particle speeds
257!-- (for particles below first vertical grid level).
258    z_p      = zu(nzb+1) - zw(nzb)
259    d_z_p_z0 = 1.0_wp / ( z_p - z0_av_global )
260
261    start_index = grid_particles(kp,jp,ip)%start_index
262    end_index   = grid_particles(kp,jp,ip)%end_index
263
264    xv = particles(1:number_of_particles)%x
265    yv = particles(1:number_of_particles)%y
266    zv = particles(1:number_of_particles)%z
267
268    DO  nb = 0, 7
269
270       i = ip
271       j = jp + block_offset(nb)%j_off
272       k = kp + block_offset(nb)%k_off
273!
274!--    Interpolate u velocity-component
275       DO  n = start_index(nb), end_index(nb)
276!
277!--       Interpolation of the u velocity component onto particle position. 
278!--       Particles are interpolation bi-linearly in the horizontal and a
279!--       linearly in the vertical. An exception is made for particles below
280!--       the first vertical grid level in case of a prandtl layer. In this
281!--       case the horizontal particle velocity components are determined using
282!--       Monin-Obukhov relations (if branch).
283!--       First, check if particle is located below first vertical grid level
284!--       above topography (Prandtl-layer height)
285          ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx
286          jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy
287!
288!--       Determine vertical index of topography top
289          k_wall = get_topography_top_index( jlog,ilog, 's' )
290
291          IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
292!
293!--          Resolved-scale horizontal particle velocity is zero below z0.
294             IF ( zv(n) - zw(k_wall) < z0_av_global )  THEN
295                u_int(n) = 0.0_wp
296             ELSE
297!
298!--             Determine the sublayer. Further used as index.
299                height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
300                                     * REAL( number_of_sublayers, KIND=wp )    &
301                                     * d_z_p_z0
302!
303!--             Calculate LOG(z/z0) for exact particle height. Therefore,   
304!--             interpolate linearly between precalculated logarithm.
305                log_z_z0_int = log_z_z0(INT(height_p))                         &
306                                 + ( height_p - INT(height_p) )                &
307                                 * ( log_z_z0(INT(height_p)+1)                 &
308                                      - log_z_z0(INT(height_p))                &
309                                   ) 
310!
311!--             Get friction velocity and momentum flux from new surface data
312!--             types.
313                IF ( surf_def_h(0)%start_index(jlog,ilog) <=                   &
314                     surf_def_h(0)%end_index(jlog,ilog) )  THEN
315                   surf_start = surf_def_h(0)%start_index(jlog,ilog)
316!--                Limit friction velocity. In narrow canyons or holes the
317!--                friction velocity can become very small, resulting in a too
318!--                large particle speed.
319                   us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 
320                   usws_int  = surf_def_h(0)%usws(surf_start)
321                ELSEIF ( surf_lsm_h%start_index(jlog,ilog) <=                  &
322                         surf_lsm_h%end_index(jlog,ilog) )  THEN
323                   surf_start = surf_lsm_h%start_index(jlog,ilog)
324                   us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 
325                   usws_int  = surf_lsm_h%usws(surf_start)
326                ELSEIF ( surf_usm_h%start_index(jlog,ilog) <=                  &
327                         surf_usm_h%end_index(jlog,ilog) )  THEN
328                   surf_start = surf_usm_h%start_index(jlog,ilog)
329                   us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 
330                   usws_int  = surf_usm_h%usws(surf_start)
331                ENDIF
332
333!
334!--             Neutral solution is applied for all situations, e.g. also for
335!--             unstable and stable situations. Even though this is not exact
336!--             this saves a lot of CPU time since several calls of intrinsic
337!--             FORTRAN procedures (LOG, ATAN) are avoided, This is justified
338!--             as sensitivity studies revealed no significant effect of
339!--             using the neutral solution also for un/stable situations.
340                u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp )           & 
341                            * log_z_z0_int - u_gtrans
342               
343             ENDIF
344!
345!--       Particle above the first grid level. Bi-linear interpolation in the
346!--       horizontal and linear interpolation in the vertical direction.
347          ELSE
348
349             = xv(n) + ( 0.5_wp - i ) * dx
350             y  = yv(n) - j * dy
351             aa = x**2          + y**2
352             bb = ( dx - x )**2 + y**2
353             cc = x**2          + ( dy - y )**2
354             dd = ( dx - x )**2 + ( dy - y )**2
355             gg = aa + bb + cc + dd
356
357             u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
358                         + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) *            &
359                         u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
360
361             IF ( k == nzt )  THEN
362                u_int(n) = u_int_l
363             ELSE
364                u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1)  &
365                            + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) *           &
366                            u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
367                u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dz *                  &
368                           ( u_int_u - u_int_l )
369             ENDIF
370
371          ENDIF
372
373       ENDDO
374
375       i = ip + block_offset(nb)%i_off
376       j = jp
377       k = kp + block_offset(nb)%k_off
378!
379!--    Same procedure for interpolation of the v velocity-component
380       DO  n = start_index(nb), end_index(nb)
381
382          ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx
383          jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy
384!
385!--       Determine vertical index of topography top
386          k_wall = get_topography_top_index( jlog,ilog, 's' )
387
388          IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
389             IF ( zv(n) - zw(k_wall) < z0_av_global )  THEN
390!
391!--             Resolved-scale horizontal particle velocity is zero below z0.
392                v_int(n) = 0.0_wp
393             ELSE       
394!
395!--             Determine the sublayer. Further used as index. Please note,
396!--             logarithmus can not be reused from above, as in in case of
397!--             topography particle on u-grid can be above surface-layer height,
398!--             whereas it can be below on v-grid.
399                height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
400                                  * REAL( number_of_sublayers, KIND=wp )       &
401                                  * d_z_p_z0
402!
403!--             Calculate LOG(z/z0) for exact particle height. Therefore,   
404!--             interpolate linearly between precalculated logarithm.
405                log_z_z0_int = log_z_z0(INT(height_p))                         &
406                                 + ( height_p - INT(height_p) )                &
407                                 * ( log_z_z0(INT(height_p)+1)                 &
408                                      - log_z_z0(INT(height_p))                &
409                                   ) 
410!
411!--             Get friction velocity and momentum flux from new surface data
412!--             types.
413                IF ( surf_def_h(0)%start_index(jlog,ilog) <=                   &
414                     surf_def_h(0)%end_index(jlog,ilog) )  THEN
415                   surf_start = surf_def_h(0)%start_index(jlog,ilog)
416!--                Limit friction velocity. In narrow canyons or holes the
417!--                friction velocity can become very small, resulting in a too
418!--                large particle speed.
419                   us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 
420                   vsws_int  = surf_def_h(0)%usws(surf_start)
421                ELSEIF ( surf_lsm_h%start_index(jlog,ilog) <=                  &
422                         surf_lsm_h%end_index(jlog,ilog) )  THEN
423                   surf_start = surf_lsm_h%start_index(jlog,ilog)
424                   us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 
425                   vsws_int  = surf_lsm_h%usws(surf_start)
426                ELSEIF ( surf_usm_h%start_index(jlog,ilog) <=                  &
427                         surf_usm_h%end_index(jlog,ilog) )  THEN
428                   surf_start = surf_usm_h%start_index(jlog,ilog)
429                   us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 
430                   vsws_int  = surf_usm_h%usws(surf_start)
431                ENDIF
432!
433!--             Neutral solution is applied for all situations, e.g. also for
434!--             unstable and stable situations. Even though this is not exact
435!--             this saves a lot of CPU time since several calls of intrinsic
436!--             FORTRAN procedures (LOG, ATAN) are avoided, This is justified
437!--             as sensitivity studies revealed no significant effect of
438!--             using the neutral solution also for un/stable situations.
439                v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp )           &
440                         * log_z_z0_int - v_gtrans
441
442             ENDIF
443
444          ELSE
445             = xv(n) - i * dx
446             y  = yv(n) + ( 0.5_wp - j ) * dy
447             aa = x**2          + y**2
448             bb = ( dx - x )**2 + y**2
449             cc = x**2          + ( dy - y )**2
450             dd = ( dx - x )**2 + ( dy - y )**2
451             gg = aa + bb + cc + dd
452
453             v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
454                       + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
455                       ) / ( 3.0_wp * gg ) - v_gtrans
456
457             IF ( k == nzt )  THEN
458                v_int(n) = v_int_l
459             ELSE
460                v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
461                          + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
462                          ) / ( 3.0_wp * gg ) - v_gtrans
463                v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dz * &
464                                  ( v_int_u - v_int_l )
465             ENDIF
466
467          ENDIF
468
469       ENDDO
470
471       i = ip + block_offset(nb)%i_off
472       j = jp + block_offset(nb)%j_off
473       k = kp - 1
474!
475!--    Same procedure for interpolation of the w velocity-component
476       DO  n = start_index(nb), end_index(nb)
477
478          IF ( vertical_particle_advection(particles(n)%group) )  THEN
479
480             = xv(n) - i * dx
481             y  = yv(n) - j * dy
482             aa = x**2          + y**2
483             bb = ( dx - x )**2 + y**2
484             cc = x**2          + ( dy - y )**2
485             dd = ( dx - x )**2 + ( dy - y )**2
486             gg = aa + bb + cc + dd
487
488             w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
489                       + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
490                       ) / ( 3.0_wp * gg )
491
492             IF ( k == nzt )  THEN
493                w_int(n) = w_int_l
494             ELSE
495                w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + &
496                            ( gg-bb ) * w(k+1,j,i+1) + &
497                            ( gg-cc ) * w(k+1,j+1,i) + &
498                            ( gg-dd ) * w(k+1,j+1,i+1) &
499                          ) / ( 3.0_wp * gg )
500                w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dz * &
501                           ( w_int_u - w_int_l )
502             ENDIF
503
504          ELSE
505
506             w_int(n) = 0.0_wp
507
508          ENDIF
509
510       ENDDO
511
512    ENDDO
513
514!-- Interpolate and calculate quantities needed for calculating the SGS
515!-- velocities
516    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
517
518       IF ( topography == 'flat' )  THEN
519
520          DO  nb = 0,7
521
522             i = ip + block_offset(nb)%i_off
523             j = jp + block_offset(nb)%j_off
524             k = kp + block_offset(nb)%k_off
525
526             DO  n = start_index(nb), end_index(nb)
527!
528!--             Interpolate TKE
529                x  = xv(n) - i * dx
530                y  = yv(n) - j * dy
531                aa = x**2          + y**2
532                bb = ( dx - x )**2 + y**2
533                cc = x**2          + ( dy - y )**2
534                dd = ( dx - x )**2 + ( dy - y )**2
535                gg = aa + bb + cc + dd
536
537                e_int_l = ( ( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1)   &
538                          + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
539                          ) / ( 3.0_wp * gg )
540
541                IF ( k+1 == nzt+1 )  THEN
542                   e_int(n) = e_int_l
543                ELSE
544                   e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
545                               ( gg - bb ) * e(k+1,j,i+1) + &
546                               ( gg - cc ) * e(k+1,j+1,i) + &
547                               ( gg - dd ) * e(k+1,j+1,i+1) &
548                            ) / ( 3.0_wp * gg )
549                   e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
550                                     ( e_int_u - e_int_l )
551                ENDIF
552!
553!--             Needed to avoid NaN particle velocities (this might not be
554!--             required any more)
555                IF ( e_int(n) <= 0.0_wp )  THEN
556                   e_int(n) = 1.0E-20_wp
557                ENDIF
558!
559!--             Interpolate the TKE gradient along x (adopt incides i,j,k and
560!--             all position variables from above (TKE))
561                de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
562                                ( gg - bb ) * de_dx(k,j,i+1) + &
563                                ( gg - cc ) * de_dx(k,j+1,i) + &
564                                ( gg - dd ) * de_dx(k,j+1,i+1) &
565                               ) / ( 3.0_wp * gg )
566
567                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
568                   de_dx_int(n) = de_dx_int_l
569                ELSE
570                   de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
571                                   ( gg - bb ) * de_dx(k+1,j,i+1) + &
572                                   ( gg - cc ) * de_dx(k+1,j+1,i) + &
573                                   ( gg - dd ) * de_dx(k+1,j+1,i+1) &
574                                  ) / ( 3.0_wp * gg )
575                   de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / dz * &
576                                              ( de_dx_int_u - de_dx_int_l )
577                ENDIF
578!
579!--             Interpolate the TKE gradient along y
580                de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
581                                ( gg - bb ) * de_dy(k,j,i+1) + &
582                                ( gg - cc ) * de_dy(k,j+1,i) + &
583                                ( gg - dd ) * de_dy(k,j+1,i+1) &
584                               ) / ( 3.0_wp * gg )
585                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
586                   de_dy_int(n) = de_dy_int_l
587                ELSE
588                   de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
589                                  ( gg - bb ) * de_dy(k+1,j,i+1) + &
590                                  ( gg - cc ) * de_dy(k+1,j+1,i) + &
591                                  ( gg - dd ) * de_dy(k+1,j+1,i+1) &
592                                  ) / ( 3.0_wp * gg )
593                   de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * &
594                                              ( de_dy_int_u - de_dy_int_l )
595                ENDIF
596
597!
598!--             Interpolate the TKE gradient along z
599                IF ( zv(n) < 0.5_wp * dz )  THEN
600                   de_dz_int(n) = 0.0_wp
601                ELSE
602                   de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
603                                   ( gg - bb ) * de_dz(k,j,i+1) + &
604                                   ( gg - cc ) * de_dz(k,j+1,i) + &
605                                   ( gg - dd ) * de_dz(k,j+1,i+1) &
606                                  ) / ( 3.0_wp * gg )
607
608                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
609                      de_dz_int(n) = de_dz_int_l
610                   ELSE
611                      de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
612                                      ( gg - bb ) * de_dz(k+1,j,i+1) + &
613                                      ( gg - cc ) * de_dz(k+1,j+1,i) + &
614                                      ( gg - dd ) * de_dz(k+1,j+1,i+1) &
615                                     ) / ( 3.0_wp * gg )
616                      de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) / dz * &
617                                                 ( de_dz_int_u - de_dz_int_l )
618                   ENDIF
619                ENDIF
620
621!
622!--             Interpolate the dissipation of TKE
623                diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
624                               ( gg - bb ) * diss(k,j,i+1) + &
625                               ( gg - cc ) * diss(k,j+1,i) + &
626                               ( gg - dd ) * diss(k,j+1,i+1) &
627                              ) / ( 3.0_wp * gg )
628
629                IF ( k == nzt )  THEN
630                   diss_int(n) = diss_int_l
631                ELSE
632                   diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
633                                  ( gg - bb ) * diss(k+1,j,i+1) + &
634                                  ( gg - cc ) * diss(k+1,j+1,i) + &
635                                  ( gg - dd ) * diss(k+1,j+1,i+1) &
636                                 ) / ( 3.0_wp * gg )
637                   diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz * &
638                                           ( diss_int_u - diss_int_l )
639                ENDIF
640
641!
642!--             Set flag for stochastic equation.
643                term_1_2(n) = 1.0_wp
644
645             ENDDO
646          ENDDO
647
648       ELSE ! non-flat topography, e.g., buildings
649
650          DO  n = 1, number_of_particles
651             i = particles(n)%x * ddx
652             j = particles(n)%y * ddy
653             k = ( zv(n) + 0.5_wp * dz * atmos_ocean_sign ) / dz  &
654                 + offset_ocean_nzt                      ! only exact if eq.dist
655!
656!--          In case that there are buildings it has to be determined
657!--          how many of the gridpoints defining the particle box are
658!--          situated within a building
659!--          gp_outside_of_building(1): i,j,k
660!--          gp_outside_of_building(2): i,j+1,k
661!--          gp_outside_of_building(3): i,j,k+1
662!--          gp_outside_of_building(4): i,j+1,k+1
663!--          gp_outside_of_building(5): i+1,j,k
664!--          gp_outside_of_building(6): i+1,j+1,k
665!--          gp_outside_of_building(7): i+1,j,k+1
666!--          gp_outside_of_building(8): i+1,j+1,k+1
667
668             gp_outside_of_building = 0
669             location = 0.0_wp
670             num_gp = 0
671
672!
673!--          Determine vertical index of topography top at (j,i)
674             k_wall = get_topography_top_index( j, i, 's' )
675!
676!--          To do: Reconsider order of computations in order to avoid
677!--          unnecessary index calculations.
678             IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
679                num_gp = num_gp + 1
680                gp_outside_of_building(1) = 1
681                location(num_gp,1) = i * dx
682                location(num_gp,2) = j * dy
683                location(num_gp,3) = k * dz - 0.5_wp * dz
684                ei(num_gp)     = e(k,j,i)
685                dissi(num_gp)  = diss(k,j,i)
686                de_dxi(num_gp) = de_dx(k,j,i)
687                de_dyi(num_gp) = de_dy(k,j,i)
688                de_dzi(num_gp) = de_dz(k,j,i)
689             ENDIF
690
691!
692!--          Determine vertical index of topography top at (j+1,i)
693             k_wall = get_topography_top_index( j+1, i, 's' )
694
695             IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
696                num_gp = num_gp + 1
697                gp_outside_of_building(2) = 1
698                location(num_gp,1) = i * dx
699                location(num_gp,2) = (j+1) * dy
700                location(num_gp,3) = k * dz - 0.5_wp * dz
701                ei(num_gp)     = e(k,j+1,i)
702                dissi(num_gp)  = diss(k,j+1,i)
703                de_dxi(num_gp) = de_dx(k,j+1,i)
704                de_dyi(num_gp) = de_dy(k,j+1,i)
705                de_dzi(num_gp) = de_dz(k,j+1,i)
706             ENDIF
707
708!
709!--          Determine vertical index of topography top at (j,i)
710             k_wall = get_topography_top_index( j, i, 's' )
711
712             IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
713                num_gp = num_gp + 1
714                gp_outside_of_building(3) = 1
715                location(num_gp,1) = i * dx
716                location(num_gp,2) = j * dy
717                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
718                ei(num_gp)     = e(k+1,j,i)
719                dissi(num_gp)  = diss(k+1,j,i)
720                de_dxi(num_gp) = de_dx(k+1,j,i)
721                de_dyi(num_gp) = de_dy(k+1,j,i)
722                de_dzi(num_gp) = de_dz(k+1,j,i)
723             ENDIF
724
725!
726!--          Determine vertical index of topography top at (j+1,i)
727             k_wall = get_topography_top_index( j+1, i, 's' )
728             IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
729                num_gp = num_gp + 1
730                gp_outside_of_building(4) = 1
731                location(num_gp,1) = i * dx
732                location(num_gp,2) = (j+1) * dy
733                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
734                ei(num_gp)     = e(k+1,j+1,i)
735                dissi(num_gp)  = diss(k+1,j+1,i)
736                de_dxi(num_gp) = de_dx(k+1,j+1,i)
737                de_dyi(num_gp) = de_dy(k+1,j+1,i)
738                de_dzi(num_gp) = de_dz(k+1,j+1,i)
739             ENDIF
740
741!
742!--          Determine vertical index of topography top at (j,i+1)
743             k_wall = get_topography_top_index( j, i+1, 's' )
744             IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
745                num_gp = num_gp + 1
746                gp_outside_of_building(5) = 1
747                location(num_gp,1) = (i+1) * dx
748                location(num_gp,2) = j * dy
749                location(num_gp,3) = k * dz - 0.5_wp * dz
750                ei(num_gp)     = e(k,j,i+1)
751                dissi(num_gp)  = diss(k,j,i+1)
752                de_dxi(num_gp) = de_dx(k,j,i+1)
753                de_dyi(num_gp) = de_dy(k,j,i+1)
754                de_dzi(num_gp) = de_dz(k,j,i+1)
755             ENDIF
756
757!
758!--          Determine vertical index of topography top at (j+1,i+1)
759             k_wall = get_topography_top_index( j+1, i+1, 's' )
760
761             IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
762                num_gp = num_gp + 1
763                gp_outside_of_building(6) = 1
764                location(num_gp,1) = (i+1) * dx
765                location(num_gp,2) = (j+1) * dy
766                location(num_gp,3) = k * dz - 0.5_wp * dz
767                ei(num_gp)     = e(k,j+1,i+1)
768                dissi(num_gp)  = diss(k,j+1,i+1)
769                de_dxi(num_gp) = de_dx(k,j+1,i+1)
770                de_dyi(num_gp) = de_dy(k,j+1,i+1)
771                de_dzi(num_gp) = de_dz(k,j+1,i+1)
772             ENDIF
773
774!
775!--          Determine vertical index of topography top at (j,i+1)
776             k_wall = get_topography_top_index( j, i+1, 's' )
777
778             IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
779                num_gp = num_gp + 1
780                gp_outside_of_building(7) = 1
781                location(num_gp,1) = (i+1) * dx
782                location(num_gp,2) = j * dy
783                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
784                ei(num_gp)     = e(k+1,j,i+1)
785                dissi(num_gp)  = diss(k+1,j,i+1)
786                de_dxi(num_gp) = de_dx(k+1,j,i+1)
787                de_dyi(num_gp) = de_dy(k+1,j,i+1)
788                de_dzi(num_gp) = de_dz(k+1,j,i+1)
789             ENDIF
790
791!
792!--          Determine vertical index of topography top at (j+1,i+1)
793             k_wall = get_topography_top_index( j+1, i+1, 's' )
794
795             IF ( k+1 > k_wall  .OR.  k_wall == 0)  THEN
796                num_gp = num_gp + 1
797                gp_outside_of_building(8) = 1
798                location(num_gp,1) = (i+1) * dx
799                location(num_gp,2) = (j+1) * dy
800                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
801                ei(num_gp)     = e(k+1,j+1,i+1)
802                dissi(num_gp)  = diss(k+1,j+1,i+1)
803                de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
804                de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
805                de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
806             ENDIF
807!
808!--          If all gridpoints are situated outside of a building, then the
809!--          ordinary interpolation scheme can be used.
810             IF ( num_gp == 8 )  THEN
811
812                = particles(n)%x - i * dx
813                y  = particles(n)%y - j * dy
814                aa = x**2          + y**2
815                bb = ( dx - x )**2 + y**2
816                cc = x**2          + ( dy - y )**2
817                dd = ( dx - x )**2 + ( dy - y )**2
818                gg = aa + bb + cc + dd
819 
820                e_int_l = ( ( gg - aa ) * e(k,j,i)   + ( gg - bb ) * e(k,j,i+1)   &
821                          + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) &
822                          ) / ( 3.0_wp * gg )
823   
824                IF ( k == nzt )  THEN
825                   e_int(n) = e_int_l
826                ELSE
827                   e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
828                               ( gg - bb ) * e(k+1,j,i+1) + &
829                               ( gg - cc ) * e(k+1,j+1,i) + &
830                               ( gg - dd ) * e(k+1,j+1,i+1) &
831                             ) / ( 3.0_wp * gg )
832                   e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
833                                       ( e_int_u - e_int_l )
834                ENDIF
835!
836!--             Needed to avoid NaN particle velocities (this might not be
837!--             required any more)
838                IF ( e_int(n) <= 0.0_wp )  THEN
839                   e_int(n) = 1.0E-20_wp
840                ENDIF
841!
842!--             Interpolate the TKE gradient along x (adopt incides i,j,k
843!--             and all position variables from above (TKE))
844                de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
845                                ( gg - bb ) * de_dx(k,j,i+1) + &
846                                ( gg - cc ) * de_dx(k,j+1,i) + &
847                                ( gg - dd ) * de_dx(k,j+1,i+1) &
848                              ) / ( 3.0_wp * gg )
849
850                IF ( ( k == nzt )  .OR.  ( k == nzb ) )  THEN
851                   de_dx_int(n) = de_dx_int_l
852                ELSE
853                   de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
854                                   ( gg - bb ) * de_dx(k+1,j,i+1) + &
855                                   ( gg - cc ) * de_dx(k+1,j+1,i) + &
856                                   ( gg - dd ) * de_dx(k+1,j+1,i+1) &
857                                 ) / ( 3.0_wp * gg )
858                   de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / &
859                                           dz * ( de_dx_int_u - de_dx_int_l )
860                ENDIF
861
862!
863!--             Interpolate the TKE gradient along y
864                de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
865                                ( gg - bb ) * de_dy(k,j,i+1) + &
866                                ( gg - cc ) * de_dy(k,j+1,i) + &
867                                ( gg - dd ) * de_dy(k,j+1,i+1) &
868                              ) / ( 3.0_wp * gg )
869                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
870                   de_dy_int(n) = de_dy_int_l
871                ELSE
872                   de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
873                                   ( gg - bb ) * de_dy(k+1,j,i+1) + &
874                                   ( gg - cc ) * de_dy(k+1,j+1,i) + &
875                                   ( gg - dd ) * de_dy(k+1,j+1,i+1) &
876                                 ) / ( 3.0_wp * gg )
877                   de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / &
878                                           dz * ( de_dy_int_u - de_dy_int_l )
879                ENDIF
880
881!
882!--             Interpolate the TKE gradient along z
883                IF ( zv(n) < 0.5_wp * dz )  THEN
884                   de_dz_int(n) = 0.0_wp
885                ELSE
886                   de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
887                                   ( gg - bb ) * de_dz(k,j,i+1) + &
888                                   ( gg - cc ) * de_dz(k,j+1,i) + &
889                                   ( gg - dd ) * de_dz(k,j+1,i+1) &
890                                 ) / ( 3.0_wp * gg )
891 
892                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
893                      de_dz_int(n) = de_dz_int_l
894                   ELSE
895                      de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
896                                      ( gg - bb ) * de_dz(k+1,j,i+1) + &
897                                      ( gg - cc ) * de_dz(k+1,j+1,i) + &
898                                      ( gg - dd ) * de_dz(k+1,j+1,i+1) &
899                                    ) / ( 3.0_wp * gg )
900                      de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) /&
901                                           dz * ( de_dz_int_u - de_dz_int_l )
902                   ENDIF
903                ENDIF
904
905!
906!--             Interpolate the dissipation of TKE
907                diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
908                               ( gg - bb ) * diss(k,j,i+1) + &
909                               ( gg - cc ) * diss(k,j+1,i) + &
910                               ( gg - dd ) * diss(k,j+1,i+1) &
911                             ) / ( 3.0_wp * gg )
912 
913                IF ( k == nzt )  THEN
914                   diss_int(n) = diss_int_l
915                ELSE
916                   diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
917                                  ( gg - bb ) * diss(k+1,j,i+1) + &
918                                  ( gg - cc ) * diss(k+1,j+1,i) + &
919                                  ( gg - dd ) * diss(k+1,j+1,i+1) &
920                                ) / ( 3.0_wp * gg )
921                   diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *&
922                                           ( diss_int_u - diss_int_l )
923                ENDIF
924!
925!--             Set flag for stochastic equation.
926                term_1_2(n) = 1.0_wp
927 
928             ELSE
929 
930!
931!--             If wall between gridpoint 1 and gridpoint 5, then
932!--             Neumann boundary condition has to be applied
933                IF ( gp_outside_of_building(1) == 1  .AND. &
934                     gp_outside_of_building(5) == 0 )  THEN
935                   num_gp = num_gp + 1
936                   location(num_gp,1) = i * dx + 0.5_wp * dx
937                   location(num_gp,2) = j * dy
938                   location(num_gp,3) = k * dz - 0.5_wp * dz
939                   ei(num_gp)     = e(k,j,i)
940                   dissi(num_gp)  = diss(k,j,i)
941                   de_dxi(num_gp) = 0.0_wp
942                   de_dyi(num_gp) = de_dy(k,j,i)
943                   de_dzi(num_gp) = de_dz(k,j,i)
944                ENDIF
945
946                IF ( gp_outside_of_building(5) == 1  .AND. &
947                     gp_outside_of_building(1) == 0 )  THEN
948                   num_gp = num_gp + 1
949                   location(num_gp,1) = i * dx + 0.5_wp * dx
950                   location(num_gp,2) = j * dy
951                   location(num_gp,3) = k * dz - 0.5_wp * dz
952                   ei(num_gp)     = e(k,j,i+1)
953                   dissi(num_gp)  = diss(k,j,i+1)
954                   de_dxi(num_gp) = 0.0_wp
955                   de_dyi(num_gp) = de_dy(k,j,i+1)
956                   de_dzi(num_gp) = de_dz(k,j,i+1)
957                ENDIF
958
959!
960!--             If wall between gridpoint 5 and gridpoint 6, then
961!--             then Neumann boundary condition has to be applied
962                IF ( gp_outside_of_building(5) == 1  .AND. &
963                     gp_outside_of_building(6) == 0 )  THEN
964                   num_gp = num_gp + 1
965                   location(num_gp,1) = (i+1) * dx
966                   location(num_gp,2) = j * dy + 0.5_wp * dy
967                   location(num_gp,3) = k * dz - 0.5_wp * dz
968                   ei(num_gp)     = e(k,j,i+1)
969                   dissi(num_gp)  = diss(k,j,i+1)
970                   de_dxi(num_gp) = de_dx(k,j,i+1)
971                   de_dyi(num_gp) = 0.0_wp
972                   de_dzi(num_gp) = de_dz(k,j,i+1)
973                ENDIF
974
975                IF ( gp_outside_of_building(6) == 1  .AND. &
976                     gp_outside_of_building(5) == 0 )  THEN
977                   num_gp = num_gp + 1
978                   location(num_gp,1) = (i+1) * dx
979                   location(num_gp,2) = j * dy + 0.5_wp * dy
980                   location(num_gp,3) = k * dz - 0.5_wp * dz
981                   ei(num_gp)     = e(k,j+1,i+1)
982                   dissi(num_gp)  = diss(k,j+1,i+1)
983                   de_dxi(num_gp) = de_dx(k,j+1,i+1)
984                   de_dyi(num_gp) = 0.0_wp
985                   de_dzi(num_gp) = de_dz(k,j+1,i+1)
986                ENDIF
987
988!
989!--             If wall between gridpoint 2 and gridpoint 6, then
990!--             Neumann boundary condition has to be applied
991                IF ( gp_outside_of_building(2) == 1  .AND. &
992                     gp_outside_of_building(6) == 0 )  THEN
993                   num_gp = num_gp + 1
994                   location(num_gp,1) = i * dx + 0.5_wp * dx
995                   location(num_gp,2) = (j+1) * dy
996                   location(num_gp,3) = k * dz - 0.5_wp * dz
997                   ei(num_gp)     = e(k,j+1,i)
998                   dissi(num_gp)  = diss(k,j+1,i)
999                   de_dxi(num_gp) = 0.0_wp
1000                   de_dyi(num_gp) = de_dy(k,j+1,i)
1001                   de_dzi(num_gp) = de_dz(k,j+1,i)
1002                ENDIF
1003
1004                IF ( gp_outside_of_building(6) == 1  .AND. &
1005                     gp_outside_of_building(2) == 0 )  THEN
1006                   num_gp = num_gp + 1
1007                   location(num_gp,1) = i * dx + 0.5_wp * dx
1008                   location(num_gp,2) = (j+1) * dy
1009                   location(num_gp,3) = k * dz - 0.5_wp * dz
1010                   ei(num_gp)     = e(k,j+1,i+1)
1011                   dissi(num_gp)  = diss(k,j+1,i+1)
1012                   de_dxi(num_gp) = 0.0_wp
1013                   de_dyi(num_gp) = de_dy(k,j+1,i+1)
1014                   de_dzi(num_gp) = de_dz(k,j+1,i+1)
1015                ENDIF
1016
1017!
1018!--             If wall between gridpoint 1 and gridpoint 2, then
1019!--             Neumann boundary condition has to be applied
1020                IF ( gp_outside_of_building(1) == 1  .AND. &
1021                     gp_outside_of_building(2) == 0 )  THEN
1022                   num_gp = num_gp + 1
1023                   location(num_gp,1) = i * dx
1024                   location(num_gp,2) = j * dy + 0.5_wp * dy
1025                   location(num_gp,3) = k * dz - 0.5_wp * dz
1026                   ei(num_gp)     = e(k,j,i)
1027                   dissi(num_gp)  = diss(k,j,i)
1028                   de_dxi(num_gp) = de_dx(k,j,i)
1029                   de_dyi(num_gp) = 0.0_wp
1030                   de_dzi(num_gp) = de_dz(k,j,i) 
1031                ENDIF
1032
1033                IF ( gp_outside_of_building(2) == 1  .AND. &
1034                     gp_outside_of_building(1) == 0 )  THEN
1035                   num_gp = num_gp + 1
1036                   location(num_gp,1) = i * dx
1037                   location(num_gp,2) = j * dy + 0.5_wp * dy
1038                   location(num_gp,3) = k * dz - 0.5_wp * dz
1039                   ei(num_gp)     = e(k,j+1,i)
1040                   dissi(num_gp)  = diss(k,j+1,i)
1041                   de_dxi(num_gp) = de_dx(k,j+1,i)
1042                   de_dyi(num_gp) = 0.0_wp
1043                   de_dzi(num_gp) = de_dz(k,j+1,i)
1044                ENDIF
1045
1046!
1047!--             If wall between gridpoint 3 and gridpoint 7, then
1048!--             Neumann boundary condition has to be applied
1049                IF ( gp_outside_of_building(3) == 1  .AND. &
1050                     gp_outside_of_building(7) == 0 )  THEN
1051                   num_gp = num_gp + 1
1052                   location(num_gp,1) = i * dx + 0.5_wp * dx
1053                   location(num_gp,2) = j * dy
1054                   location(num_gp,3) = k * dz + 0.5_wp * dz
1055                   ei(num_gp)     = e(k+1,j,i)
1056                   dissi(num_gp)  = diss(k+1,j,i)
1057                   de_dxi(num_gp) = 0.0_wp
1058                   de_dyi(num_gp) = de_dy(k+1,j,i)
1059                   de_dzi(num_gp) = de_dz(k+1,j,i) 
1060                ENDIF
1061
1062                IF ( gp_outside_of_building(7) == 1  .AND. &
1063                     gp_outside_of_building(3) == 0 )  THEN
1064                   num_gp = num_gp + 1
1065                   location(num_gp,1) = i * dx + 0.5_wp * dx
1066                   location(num_gp,2) = j * dy
1067                   location(num_gp,3) = k * dz + 0.5_wp * dz
1068                   ei(num_gp)     = e(k+1,j,i+1)
1069                   dissi(num_gp)  = diss(k+1,j,i+1)
1070                   de_dxi(num_gp) = 0.0_wp
1071                   de_dyi(num_gp) = de_dy(k+1,j,i+1)
1072                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
1073                ENDIF
1074
1075!
1076!--             If wall between gridpoint 7 and gridpoint 8, then
1077!--             Neumann boundary condition has to be applied
1078                IF ( gp_outside_of_building(7) == 1  .AND. &
1079                     gp_outside_of_building(8) == 0 )  THEN
1080                   num_gp = num_gp + 1
1081                   location(num_gp,1) = (i+1) * dx
1082                   location(num_gp,2) = j * dy + 0.5_wp * dy
1083                   location(num_gp,3) = k * dz + 0.5_wp * dz
1084                   ei(num_gp)     = e(k+1,j,i+1)
1085                   dissi(num_gp)  = diss(k+1,j,i+1)
1086                   de_dxi(num_gp) = de_dx(k+1,j,i+1)
1087                   de_dyi(num_gp) = 0.0_wp
1088                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
1089                ENDIF
1090
1091                IF ( gp_outside_of_building(8) == 1  .AND. &
1092                     gp_outside_of_building(7) == 0 )  THEN
1093                   num_gp = num_gp + 1
1094                   location(num_gp,1) = (i+1) * dx
1095                   location(num_gp,2) = j * dy + 0.5_wp * dy
1096                   location(num_gp,3) = k * dz + 0.5_wp * dz
1097                   ei(num_gp)     = e(k+1,j+1,i+1)
1098                   dissi(num_gp)  = diss(k+1,j+1,i+1)
1099                   de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
1100                   de_dyi(num_gp) = 0.0_wp
1101                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
1102                ENDIF
1103
1104!
1105!--             If wall between gridpoint 4 and gridpoint 8, then
1106!--             Neumann boundary condition has to be applied
1107                IF ( gp_outside_of_building(4) == 1  .AND. &
1108                     gp_outside_of_building(8) == 0 )  THEN
1109                   num_gp = num_gp + 1
1110                   location(num_gp,1) = i * dx + 0.5_wp * dx
1111                   location(num_gp,2) = (j+1) * dy
1112                   location(num_gp,3) = k * dz + 0.5_wp * dz
1113                   ei(num_gp)     = e(k+1,j+1,i)
1114                   dissi(num_gp)  = diss(k+1,j+1,i)
1115                   de_dxi(num_gp) = 0.0_wp
1116                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
1117                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
1118                ENDIF
1119
1120                IF ( gp_outside_of_building(8) == 1  .AND. &
1121                     gp_outside_of_building(4) == 0 )  THEN
1122                   num_gp = num_gp + 1
1123                   location(num_gp,1) = i * dx + 0.5_wp * dx
1124                   location(num_gp,2) = (j+1) * dy
1125                   location(num_gp,3) = k * dz + 0.5_wp * dz
1126                   ei(num_gp)     = e(k+1,j+1,i+1)
1127                   dissi(num_gp)  = diss(k+1,j+1,i+1)
1128                   de_dxi(num_gp) = 0.0_wp
1129                   de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
1130                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
1131                ENDIF
1132
1133!
1134!--             If wall between gridpoint 3 and gridpoint 4, then
1135!--             Neumann boundary condition has to be applied
1136                IF ( gp_outside_of_building(3) == 1  .AND. &
1137                     gp_outside_of_building(4) == 0 )  THEN
1138                   num_gp = num_gp + 1
1139                   location(num_gp,1) = i * dx
1140                   location(num_gp,2) = j * dy + 0.5_wp * dy
1141                   location(num_gp,3) = k * dz + 0.5_wp * dz
1142                   ei(num_gp)     = e(k+1,j,i)
1143                   dissi(num_gp)  = diss(k+1,j,i)
1144                   de_dxi(num_gp) = de_dx(k+1,j,i)
1145                   de_dyi(num_gp) = 0.0_wp
1146                   de_dzi(num_gp) = de_dz(k+1,j,i)
1147                ENDIF
1148
1149                IF ( gp_outside_of_building(4) == 1  .AND. &
1150                     gp_outside_of_building(3) == 0 )  THEN
1151                   num_gp = num_gp + 1
1152                   location(num_gp,1) = i * dx
1153                   location(num_gp,2) = j * dy + 0.5_wp * dy
1154                   location(num_gp,3) = k * dz + 0.5_wp * dz
1155                   ei(num_gp)     = e(k+1,j+1,i)
1156                   dissi(num_gp)  = diss(k+1,j+1,i)
1157                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
1158                   de_dyi(num_gp) = 0.0_wp
1159                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
1160                ENDIF
1161
1162!
1163!--             If wall between gridpoint 1 and gridpoint 3, then
1164!--             Neumann boundary condition has to be applied
1165!--             (only one case as only building beneath is possible)
1166                IF ( gp_outside_of_building(1) == 0  .AND. &
1167                     gp_outside_of_building(3) == 1 )  THEN
1168                   num_gp = num_gp + 1
1169                   location(num_gp,1) = i * dx
1170                   location(num_gp,2) = j * dy
1171                   location(num_gp,3) = k * dz
1172                   ei(num_gp)     = e(k+1,j,i)
1173                   dissi(num_gp)  = diss(k+1,j,i)
1174                   de_dxi(num_gp) = de_dx(k+1,j,i)
1175                   de_dyi(num_gp) = de_dy(k+1,j,i)
1176                   de_dzi(num_gp) = 0.0_wp
1177                ENDIF
1178
1179!
1180!--             If wall between gridpoint 5 and gridpoint 7, then
1181!--             Neumann boundary condition has to be applied
1182!--             (only one case as only building beneath is possible)
1183                IF ( gp_outside_of_building(5) == 0  .AND. &
1184                     gp_outside_of_building(7) == 1 )  THEN
1185                   num_gp = num_gp + 1
1186                   location(num_gp,1) = (i+1) * dx
1187                   location(num_gp,2) = j * dy
1188                   location(num_gp,3) = k * dz
1189                   ei(num_gp)     = e(k+1,j,i+1)
1190                   dissi(num_gp)  = diss(k+1,j,i+1)
1191                   de_dxi(num_gp) = de_dx(k+1,j,i+1)
1192                   de_dyi(num_gp) = de_dy(k+1,j,i+1)
1193                   de_dzi(num_gp) = 0.0_wp
1194                ENDIF
1195
1196!
1197!--             If wall between gridpoint 2 and gridpoint 4, then
1198!--             Neumann boundary condition has to be applied
1199!--             (only one case as only building beneath is possible)
1200                IF ( gp_outside_of_building(2) == 0  .AND. &
1201                     gp_outside_of_building(4) == 1 )  THEN
1202                   num_gp = num_gp + 1
1203                   location(num_gp,1) = i * dx
1204                   location(num_gp,2) = (j+1) * dy
1205                   location(num_gp,3) = k * dz
1206                   ei(num_gp)     = e(k+1,j+1,i)
1207                   dissi(num_gp)  = diss(k+1,j+1,i)
1208                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
1209                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
1210                   de_dzi(num_gp) = 0.0_wp
1211                ENDIF
1212
1213!
1214!--             If wall between gridpoint 6 and gridpoint 8, then
1215!--             Neumann boundary condition has to be applied
1216!--             (only one case as only building beneath is possible)
1217                IF ( gp_outside_of_building(6) == 0  .AND. &
1218                     gp_outside_of_building(8) == 1 )  THEN
1219                   num_gp = num_gp + 1
1220                   location(num_gp,1) = (i+1) * dx
1221                   location(num_gp,2) = (j+1) * dy
1222                   location(num_gp,3) = k * dz
1223                   ei(num_gp)     = e(k+1,j+1,i+1)
1224                   dissi(num_gp)  = diss(k+1,j+1,i+1)
1225                   de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
1226                   de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
1227                   de_dzi(num_gp) = 0.0_wp
1228                ENDIF
1229 
1230!
1231!--             Carry out the interpolation
1232                IF ( num_gp == 1 )  THEN
1233!
1234!--                If only one of the gridpoints is situated outside of the
1235!--                building, it follows that the values at the particle
1236!--                location are the same as the gridpoint values
1237                   e_int(n)     = ei(num_gp)   
1238                   diss_int(n)  = dissi(num_gp) 
1239                   de_dx_int(n) = de_dxi(num_gp)
1240                   de_dy_int(n) = de_dyi(num_gp)
1241                   de_dz_int(n) = de_dzi(num_gp)
1242!
1243!--                Set flag for stochastic equation which disables calculation
1244!--                of drift and memory term near topography.
1245                   term_1_2(n) = 0.0_wp
1246                ELSE IF ( num_gp > 1 )  THEN
1247 
1248                   d_sum = 0.0_wp
1249!
1250!--                Evaluation of the distances between the gridpoints
1251!--                contributing to the interpolated values, and the particle
1252!--                location
1253                   DO  agp = 1, num_gp
1254                      d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2  &
1255                                   + ( particles(n)%y-location(agp,2) )**2  &
1256                                   + ( zv(n)-location(agp,3) )**2
1257                      d_sum        = d_sum + d_gp_pl(agp)
1258                   ENDDO
1259 
1260!
1261!--                Finally the interpolation can be carried out
1262                   e_int(n)     = 0.0_wp 
1263                   diss_int(n)  = 0.0_wp 
1264                   de_dx_int(n) = 0.0_wp 
1265                   de_dy_int(n) = 0.0_wp 
1266                   de_dz_int(n) = 0.0_wp 
1267                   DO  agp = 1, num_gp
1268                      e_int(n)     = e_int(n)     + ( d_sum - d_gp_pl(agp) ) * &
1269                                             ei(agp) / ( (num_gp-1) * d_sum )
1270                      diss_int(n)  = diss_int(n)  + ( d_sum - d_gp_pl(agp) ) * &
1271                                          dissi(agp) / ( (num_gp-1) * d_sum )
1272                      de_dx_int(n) = de_dx_int(n) + ( d_sum - d_gp_pl(agp) ) * &
1273                                         de_dxi(agp) / ( (num_gp-1) * d_sum )
1274                      de_dy_int(n) = de_dy_int(n) + ( d_sum - d_gp_pl(agp) ) * &
1275                                         de_dyi(agp) / ( (num_gp-1) * d_sum )
1276                      de_dz_int(n) = de_dz_int(n) + ( d_sum - d_gp_pl(agp) ) * &
1277                                         de_dzi(agp) / ( (num_gp-1) * d_sum )
1278                   ENDDO
1279 
1280                ENDIF
1281                e_int(n)     = MAX( 1E-20_wp, e_int(n)     )
1282                diss_int(n)  = MAX( 1E-20_wp, diss_int(n)  )
1283                de_dx_int(n) = MAX( 1E-20_wp, de_dx_int(n) )
1284                de_dy_int(n) = MAX( 1E-20_wp, de_dy_int(n) ) 
1285                de_dz_int(n) = MAX( 1E-20_wp, de_dz_int(n) ) 
1286!
1287!--             Set flag for stochastic equation which disables calculation
1288!--             of drift and memory term near topography.
1289                term_1_2(n) = 0.0_wp
1290             ENDIF
1291          ENDDO
1292       ENDIF
1293
1294       DO nb = 0,7
1295          i = ip + block_offset(nb)%i_off
1296          j = jp + block_offset(nb)%j_off
1297          k = kp + block_offset(nb)%k_off
1298
1299          DO  n = start_index(nb), end_index(nb)
1300!
1301!--          Vertical interpolation of the horizontally averaged SGS TKE and
1302!--          resolved-scale velocity variances and use the interpolated values
1303!--          to calculate the coefficient fs, which is a measure of the ratio
1304!--          of the subgrid-scale turbulent kinetic energy to the total amount
1305!--          of turbulent kinetic energy.
1306             IF ( k == 0 )  THEN
1307                e_mean_int = hom(0,1,8,0)
1308             ELSE
1309                e_mean_int = hom(k,1,8,0) +                                    &
1310                                           ( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
1311                                           ( zu(k+1) - zu(k) ) *               &
1312                                           ( zv(n) - zu(k) )
1313             ENDIF
1314
1315             kw = kp - 1
1316
1317             IF ( k == 0 )  THEN
1318                aa  = hom(k+1,1,30,0)  * ( zv(n) / &
1319                                         ( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
1320                bb  = hom(k+1,1,31,0)  * ( zv(n) / &
1321                                         ( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
1322                cc  = hom(kw+1,1,32,0) * ( zv(n) / &
1323                                         ( 1.0_wp * ( zw(kw+1) - zw(kw) ) ) )
1324             ELSE
1325                aa  = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) *    &
1326                           ( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
1327                bb  = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) *    &
1328                           ( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
1329                cc  = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) *   &
1330                           ( ( zv(n) - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
1331             ENDIF
1332
1333             vv_int = ( 1.0_wp / 3.0_wp ) * ( aa + bb + cc )
1334!
1335!--          Needed to avoid NaN particle velocities. The value of 1.0 is just
1336!--          an educated guess for the given case.
1337             IF ( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int == 0.0_wp )  THEN
1338                fs_int(n) = 1.0_wp
1339             ELSE
1340                fs_int(n) = ( 2.0_wp / 3.0_wp ) * e_mean_int /                 &
1341                            ( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int )
1342             ENDIF
1343
1344          ENDDO
1345       ENDDO
1346
1347       DO  n = 1, number_of_particles
1348
1349         rg(n,1) = random_gauss( iran_part, 5.0_wp )
1350         rg(n,2) = random_gauss( iran_part, 5.0_wp )
1351         rg(n,3) = random_gauss( iran_part, 5.0_wp )
1352
1353       ENDDO
1354
1355       DO  n = 1, number_of_particles
1356!
1357!--       Calculate the Lagrangian timescale according to Weil et al. (2004).
1358          lagr_timescale = ( 4.0_wp * e_int(n) + 1E-20_wp ) / &
1359                           ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp )
1360
1361!
1362!--       Calculate the next particle timestep. dt_gap is the time needed to
1363!--       complete the current LES timestep.
1364          dt_gap = dt_3d - particles(n)%dt_sum
1365          dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale, dt_gap )
1366
1367!
1368!--       The particle timestep should not be too small in order to prevent
1369!--       the number of particle timesteps of getting too large
1370          IF ( dt_particle(n) < dt_min_part  .AND.  dt_min_part < dt_gap )  THEN
1371             dt_particle(n) = dt_min_part
1372          ENDIF
1373
1374!
1375!--       Calculate the SGS velocity components
1376          IF ( particles(n)%age == 0.0_wp )  THEN
1377!
1378!--          For new particles the SGS components are derived from the SGS
1379!--          TKE. Limit the Gaussian random number to the interval
1380!--          [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities
1381!--          from becoming unrealistically large.
1382             particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) *   &
1383                                  ( rg(n,1) - 1.0_wp )
1384             particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) *   &
1385                                  ( rg(n,2) - 1.0_wp )
1386             particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) *   &
1387                                  ( rg(n,3) - 1.0_wp )
1388
1389          ELSE
1390!
1391!--          Restriction of the size of the new timestep: compared to the
1392!--          previous timestep the increase must not exceed 200%
1393
1394             dt_particle_m = particles(n)%age - particles(n)%age_m
1395             IF ( dt_particle(n) > 2.0_wp * dt_particle_m )  THEN
1396                dt_particle(n) = 2.0_wp * dt_particle_m
1397             ENDIF
1398
1399!
1400!--          For old particles the SGS components are correlated with the
1401!--          values from the previous timestep. Random numbers have also to
1402!--          be limited (see above).
1403!--          As negative values for the subgrid TKE are not allowed, the
1404!--          change of the subgrid TKE with time cannot be smaller than
1405!--          -e_int(n)/dt_particle. This value is used as a lower boundary
1406!--          value for the change of TKE
1407
1408             de_dt_min = - e_int(n) / dt_particle(n)
1409
1410             de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m
1411
1412             IF ( de_dt < de_dt_min )  THEN
1413                de_dt = de_dt_min
1414             ENDIF
1415
1416             CALL weil_stochastic_eq(particles(n)%rvar1, fs_int(n), e_int(n),  & 
1417                                     de_dx_int(n), de_dt, diss_int(n),         &
1418                                     dt_particle(n), rg(n,1), term_1_2(n) )
1419
1420             CALL weil_stochastic_eq(particles(n)%rvar2, fs_int(n), e_int(n),  & 
1421                                     de_dy_int(n), de_dt, diss_int(n),         &
1422                                     dt_particle(n), rg(n,2), term_1_2(n) )
1423
1424             CALL weil_stochastic_eq(particles(n)%rvar3, fs_int(n), e_int(n),  & 
1425                                     de_dz_int(n), de_dt, diss_int(n),         &
1426                                     dt_particle(n), rg(n,3), term_1_2(n) )
1427
1428          ENDIF
1429
1430          u_int(n) = u_int(n) + particles(n)%rvar1
1431          v_int(n) = v_int(n) + particles(n)%rvar2
1432          w_int(n) = w_int(n) + particles(n)%rvar3
1433!
1434!--       Store the SGS TKE of the current timelevel which is needed for
1435!--       for calculating the SGS particle velocities at the next timestep
1436          particles(n)%e_m = e_int(n)
1437       ENDDO
1438
1439    ELSE
1440!
1441!--    If no SGS velocities are used, only the particle timestep has to
1442!--    be set
1443       dt_particle = dt_3d
1444
1445    ENDIF
1446!
1447!-- Store the old age of the particle ( needed to prevent that a
1448!-- particle crosses several PEs during one timestep, and for the
1449!-- evaluation of the subgrid particle velocity fluctuations )
1450    particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age
1451
1452    dens_ratio = particle_groups(particles(1:number_of_particles)%group)%density_ratio
1453
1454    IF ( ANY( dens_ratio == 0.0_wp ) )  THEN
1455       DO  n = 1, number_of_particles
1456
1457!
1458!--       Particle advection
1459          IF ( dens_ratio(n) == 0.0_wp )  THEN
1460!
1461!--          Pure passive transport (without particle inertia)
1462             particles(n)%x = xv(n) + u_int(n) * dt_particle(n)
1463             particles(n)%y = yv(n) + v_int(n) * dt_particle(n)
1464             particles(n)%z = zv(n) + w_int(n) * dt_particle(n)
1465
1466             particles(n)%speed_x = u_int(n)
1467             particles(n)%speed_y = v_int(n)
1468             particles(n)%speed_z = w_int(n)
1469
1470          ELSE
1471!
1472!--          Transport of particles with inertia
1473             particles(n)%x = particles(n)%x + particles(n)%speed_x * &
1474                                               dt_particle(n)
1475             particles(n)%y = particles(n)%y + particles(n)%speed_y * &
1476                                               dt_particle(n)
1477             particles(n)%z = particles(n)%z + particles(n)%speed_z * &
1478                                               dt_particle(n)
1479
1480!
1481!--          Update of the particle velocity
1482             IF ( cloud_droplets )  THEN
1483!
1484!--             Terminal velocity is computed for vertical direction (Rogers et
1485!--             al., 1993, J. Appl. Meteorol.)
1486                diameter = particles(n)%radius * 2000.0_wp !diameter in mm
1487                IF ( diameter <= d0_rog )  THEN
1488                   w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
1489                ELSE
1490                   w_s = a_rog - b_rog * EXP( -c_rog * diameter )
1491                ENDIF
1492
1493!
1494!--             If selected, add random velocities following Soelch and Kaercher
1495!--             (2010, Q. J. R. Meteorol. Soc.)
1496                IF ( use_sgs_for_particles )  THEN
1497                   lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
1498                   RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
1499                   sigma          = SQRT( e(kp,jp,ip) )
1500
1501                   rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
1502                   rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
1503                   rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
1504
1505                   particles(n)%rvar1 = RL * particles(n)%rvar1 +              &
1506                                        SQRT( 1.0_wp - RL**2 ) * sigma * rg1
1507                   particles(n)%rvar2 = RL * particles(n)%rvar2 +              &
1508                                        SQRT( 1.0_wp - RL**2 ) * sigma * rg2
1509                   particles(n)%rvar3 = RL * particles(n)%rvar3 +              &
1510                                        SQRT( 1.0_wp - RL**2 ) * sigma * rg3
1511
1512                   particles(n)%speed_x = u_int(n) + particles(n)%rvar1
1513                   particles(n)%speed_y = v_int(n) + particles(n)%rvar2
1514                   particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
1515                ELSE
1516                   particles(n)%speed_x = u_int(n)
1517                   particles(n)%speed_y = v_int(n)
1518                   particles(n)%speed_z = w_int(n) - w_s
1519                ENDIF
1520
1521             ELSE
1522
1523                IF ( use_sgs_for_particles )  THEN
1524                   exp_arg  = particle_groups(particles(n)%group)%exp_arg
1525                   exp_term = EXP( -exp_arg * dt_particle(n) )
1526                ELSE
1527                   exp_arg  = particle_groups(particles(n)%group)%exp_arg
1528                   exp_term = particle_groups(particles(n)%group)%exp_term
1529                ENDIF
1530                particles(n)%speed_x = particles(n)%speed_x * exp_term +         &
1531                                       u_int(n) * ( 1.0_wp - exp_term )
1532                particles(n)%speed_y = particles(n)%speed_y * exp_term +         &
1533                                       v_int(n) * ( 1.0_wp - exp_term )
1534                particles(n)%speed_z = particles(n)%speed_z * exp_term +         &
1535                                       ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * &
1536                                       g / exp_arg ) * ( 1.0_wp - exp_term )
1537             ENDIF
1538
1539          ENDIF
1540
1541       ENDDO
1542   
1543    ELSE
1544
1545       DO  n = 1, number_of_particles
1546
1547!--       Transport of particles with inertia
1548          particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)
1549          particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)
1550          particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)
1551!
1552!--       Update of the particle velocity
1553          IF ( cloud_droplets )  THEN
1554!
1555!--          Terminal velocity is computed for vertical direction (Rogers et al.,
1556!--          1993, J. Appl. Meteorol.)
1557             diameter = particles(n)%radius * 2000.0_wp !diameter in mm
1558             IF ( diameter <= d0_rog )  THEN
1559                w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
1560             ELSE
1561                w_s = a_rog - b_rog * EXP( -c_rog * diameter )
1562             ENDIF
1563
1564!
1565!--          If selected, add random velocities following Soelch and Kaercher
1566!--          (2010, Q. J. R. Meteorol. Soc.)
1567             IF ( use_sgs_for_particles )  THEN
1568                 lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
1569                 RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
1570                 sigma          = SQRT( e(kp,jp,ip) )
1571
1572                 rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
1573                 rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
1574                 rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
1575
1576                 particles(n)%rvar1 = RL * particles(n)%rvar1 +                &
1577                                      SQRT( 1.0_wp - RL**2 ) * sigma * rg1
1578                 particles(n)%rvar2 = RL * particles(n)%rvar2 +                &
1579                                      SQRT( 1.0_wp - RL**2 ) * sigma * rg2
1580                 particles(n)%rvar3 = RL * particles(n)%rvar3 +                &
1581                                      SQRT( 1.0_wp - RL**2 ) * sigma * rg3
1582
1583                 particles(n)%speed_x = u_int(n) + particles(n)%rvar1
1584                 particles(n)%speed_y = v_int(n) + particles(n)%rvar2
1585                 particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
1586             ELSE
1587                 particles(n)%speed_x = u_int(n)
1588                 particles(n)%speed_y = v_int(n)
1589                 particles(n)%speed_z = w_int(n) - w_s
1590             ENDIF
1591
1592          ELSE
1593
1594             IF ( use_sgs_for_particles )  THEN
1595                exp_arg  = particle_groups(particles(n)%group)%exp_arg
1596                exp_term = EXP( -exp_arg * dt_particle(n) )
1597             ELSE
1598                exp_arg  = particle_groups(particles(n)%group)%exp_arg
1599                exp_term = particle_groups(particles(n)%group)%exp_term
1600             ENDIF
1601             particles(n)%speed_x = particles(n)%speed_x * exp_term +             &
1602                                    u_int(n) * ( 1.0_wp - exp_term )
1603             particles(n)%speed_y = particles(n)%speed_y * exp_term +             &
1604                                    v_int(n) * ( 1.0_wp - exp_term )
1605             particles(n)%speed_z = particles(n)%speed_z * exp_term +             &
1606                                    ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &
1607                                    exp_arg ) * ( 1.0_wp - exp_term )
1608          ENDIF
1609
1610       ENDDO
1611
1612    ENDIF
1613
1614    DO  n = 1, number_of_particles
1615!
1616!--    Increment the particle age and the total time that the particle
1617!--    has advanced within the particle timestep procedure
1618       particles(n)%age    = particles(n)%age    + dt_particle(n)
1619       particles(n)%dt_sum = particles(n)%dt_sum + dt_particle(n)
1620
1621!
1622!--    Check whether there is still a particle that has not yet completed
1623!--    the total LES timestep
1624       IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp )  THEN
1625          dt_3d_reached_l = .FALSE.
1626       ENDIF
1627
1628    ENDDO
1629
1630    CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
1631
1632
1633 END SUBROUTINE lpm_advec
1634
1635! Description:
1636! ------------
1637!> Calculation of subgrid-scale particle speed using the stochastic model
1638!> of Weil et al. (2004, JAS, 61, 2877-2887).
1639!------------------------------------------------------------------------------!
1640 SUBROUTINE weil_stochastic_eq( v_sgs, fs_n, e_n, dedxi_n, dedt_n, diss_n,     &
1641                                dt_n, rg_n, fac )
1642
1643    USE kinds
1644
1645    USE particle_attributes,                                                   &
1646        ONLY:  c_0, sgs_wf_part
1647
1648    IMPLICIT NONE
1649
1650    REAL(wp) ::  a1      !< dummy argument
1651    REAL(wp) ::  dedt_n  !< time derivative of TKE at particle position
1652    REAL(wp) ::  dedxi_n !< horizontal derivative of TKE at particle position
1653    REAL(wp) ::  diss_n  !< dissipation at particle position
1654    REAL(wp) ::  dt_n    !< particle timestep
1655    REAL(wp) ::  e_n     !< TKE at particle position
1656    REAL(wp) ::  fac     !< flag to identify adjacent topography
1657    REAL(wp) ::  fs_n    !< weighting factor to prevent that subgrid-scale particle speed becomes too large
1658    REAL(wp) ::  sgs_w   !< constant (1/3)
1659    REAL(wp) ::  rg_n    !< random number
1660    REAL(wp) ::  term1   !< memory term
1661    REAL(wp) ::  term2   !< drift correction term
1662    REAL(wp) ::  term3   !< random term
1663    REAL(wp) ::  v_sgs   !< subgrid-scale velocity component
1664
1665!-- At first, limit TKE to a small non-zero number, in order to prevent
1666!-- the occurrence of extremely large SGS-velocities in case TKE is zero,
1667!-- (could occur at the simulation begin).
1668    e_n = MAX( e_n, 1E-20_wp )
1669!
1670!-- Please note, terms 1 and 2 (drift and memory term, respectively) are
1671!-- multiplied by a flag to switch of both terms near topography.
1672!-- This is necessary, as both terms may cause a subgrid-scale velocity build up
1673!-- if particles are trapped in regions with very small TKE, e.g. in narrow street
1674!-- canyons resolved by only a few grid points. Hence, term 1 and term 2 are
1675!-- disabled if one of the adjacent grid points belongs to topography.
1676!-- Moreover, in this case, the  previous subgrid-scale component is also set
1677!-- to zero.
1678
1679    a1 = fs_n * c_0 * diss_n
1680!
1681!-- Memory term
1682    term1 = - a1 * v_sgs * dt_n / ( 4.0_wp * sgs_wf_part * e_n + 1E-20_wp )    &
1683                 * fac
1684!
1685!-- Drift correction term
1686    term2 = ( ( dedt_n * v_sgs / e_n ) + dedxi_n ) * 0.5_wp * dt_n              &
1687                 * fac
1688!
1689!-- Random term
1690    term3 = SQRT( MAX( a1, 1E-20 ) ) * ( rg_n - 1.0_wp ) * SQRT( dt_n )
1691!
1692!-- In cese one of the adjacent grid-boxes belongs to topograhy, the previous
1693!-- subgrid-scale velocity component is set to zero, in order to prevent a
1694!-- velocity build-up.
1695
1696!-- This case, set also previous subgrid-scale component to zero.
1697    v_sgs = v_sgs * fac + term1 + term2 + term3
1698
1699 END SUBROUTINE weil_stochastic_eq
Note: See TracBrowser for help on using the repository browser.