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

Last change on this file since 2564 was 2417, checked in by suehring, 7 years ago

Major bugfix in modeling SGS particle speeds.

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