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

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

Adjustments according new topography and surface-modelling concept implemented

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