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

Last change on this file since 2100 was 2100, checked in by suehring, 5 years ago

Prevent extremely large SGS-velocities by limiting TKE in stochastic terms.

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