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

Last change on this file since 1826 was 1823, checked in by hoffmann, 8 years ago

last commit documented

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