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

Last change on this file since 1985 was 1937, checked in by suehring, 8 years ago

last commit documented

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