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

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

changes in LPM and bulk cloud microphysics

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