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

Last change on this file since 2254 was 2233, checked in by suehring, 7 years ago

last commit documented

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