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

Last change on this file since 1889 was 1889, checked in by suehring, 6 years ago

last commit documented

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