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

Last change on this file since 2214 was 2101, checked in by suehring, 7 years ago

last commit documented

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