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

Last change on this file since 1310 was 1310, checked in by raasch, 7 years ago

update of GPL copyright

  • Property svn:keywords set to Id
File size: 47.7 KB
Line 
1 SUBROUTINE lpm_advec
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: lpm_advec.f90 1310 2014-03-14 08:01:56Z raasch $
27!
28! 1036 2012-10-22 13:43:42Z raasch
29! code put under GPL (PALM 3.9)
30!
31! 849 2012-03-15 10:35:09Z raasch
32! initial revision (former part of advec_particles)
33!
34!
35! Description:
36! ------------
37! Calculation of new particle positions due to advection using a simple Euler
38! scheme. Particles may feel inertia effects. SGS transport can be included
39! using the stochastic model of Weil et al. (2004, JAS, 61, 2877-2887).
40!------------------------------------------------------------------------------!
41
42    USE arrays_3d
43    USE control_parameters
44    USE grid_variables
45    USE indices
46    USE particle_attributes
47    USE statistics
48
49    IMPLICIT NONE
50
51    INTEGER ::  i, j, k, n
52
53    REAL ::  aa, bb, cc, dd, dens_ratio, exp_arg, exp_term, gg, u_int,  &
54             u_int_l, u_int_u, v_int, v_int_l, v_int_u, w_int, w_int_l, &
55             w_int_u, x, y
56
57
58    INTEGER ::  agp, kw, num_gp
59    INTEGER ::  gp_outside_of_building(1:8)
60
61    REAL ::  d_sum, de_dx_int, de_dx_int_l, de_dx_int_u, de_dy_int,       &
62             de_dy_int_l, de_dy_int_u, de_dt, de_dt_min, de_dz_int,       &
63             de_dz_int_l, de_dz_int_u, diss_int, diss_int_l, diss_int_u,  &
64             dt_gap, dt_particle, dt_particle_m, e_int, e_int_l, e_int_u, &
65             e_mean_int, fs_int, lagr_timescale, random_gauss, vv_int
66
67    REAL ::  location(1:30,1:3)
68    REAL, DIMENSION(1:30) ::  de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei
69
70
71    DO  n = 1, number_of_particles
72
73!
74!--    Move particle only if the LES timestep has not (approximately) been
75!--    reached
76       IF ( ( dt_3d - particles(n)%dt_sum ) < 1E-8 )  CYCLE
77
78!
79!--    Interpolate u velocity-component, determine left, front, bottom
80!--    index of u-array
81       i = ( particles(n)%x + 0.5 * dx ) * ddx
82       j =   particles(n)%y * ddy
83       k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
84           + offset_ocean_nzt                     ! only exact if equidistant
85
86!
87!--    Interpolation of the velocity components in the xy-plane
88       x  = particles(n)%x + ( 0.5 - i ) * dx
89       y  = particles(n)%y - j * dy
90       aa = x**2          + y**2
91       bb = ( dx - x )**2 + y**2
92       cc = x**2          + ( dy - y )**2
93       dd = ( dx - x )**2 + ( dy - y )**2
94       gg = aa + bb + cc + dd
95
96       u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
97                 + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) &
98                 ) / ( 3.0 * gg ) - u_gtrans
99       IF ( k+1 == nzt+1 )  THEN
100          u_int = u_int_l
101       ELSE
102          u_int_u = ( ( gg-aa ) * u(k+1,j,i)   + ( gg-bb ) * u(k+1,j,i+1)    &
103                    + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1)  &
104                    ) / ( 3.0 * gg ) - u_gtrans
105          u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * &
106                            ( u_int_u - u_int_l )
107       ENDIF
108
109!
110!--    Same procedure for interpolation of the v velocity-component (adopt
111!--    index k from u velocity-component)
112       i =   particles(n)%x * ddx
113       j = ( particles(n)%y + 0.5 * dy ) * ddy
114
115       x  = particles(n)%x - i * dx
116       y  = particles(n)%y + ( 0.5 - j ) * dy
117       aa = x**2          + y**2
118       bb = ( dx - x )**2 + y**2
119       cc = x**2          + ( dy - y )**2
120       dd = ( dx - x )**2 + ( dy - y )**2
121       gg = aa + bb + cc + dd
122
123       v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
124                 + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
125                 ) / ( 3.0 * gg ) - v_gtrans
126       IF ( k+1 == nzt+1 )  THEN
127          v_int = v_int_l
128       ELSE
129          v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)    &
130                    + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1)  &
131                    ) / ( 3.0 * gg ) - v_gtrans
132          v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * &
133                            ( v_int_u - v_int_l )
134       ENDIF
135
136!
137!--    Same procedure for interpolation of the w velocity-component (adopt
138!--    index i from v velocity-component)
139       IF ( vertical_particle_advection(particles(n)%group) )  THEN
140          j = particles(n)%y * ddy
141          k = particles(n)%z / dz + offset_ocean_nzt_m1
142
143          x  = particles(n)%x - i * dx
144          y  = particles(n)%y - j * dy
145          aa = x**2          + y**2
146          bb = ( dx - x )**2 + y**2
147          cc = x**2          + ( dy - y )**2
148          dd = ( dx - x )**2 + ( dy - y )**2
149          gg = aa + bb + cc + dd
150
151          w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)    &
152                    + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1)  &
153                       ) / ( 3.0 * gg )
154          IF ( k+1 == nzt+1 )  THEN
155             w_int = w_int_l
156          ELSE
157             w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + &
158                         ( gg-bb ) * w(k+1,j,i+1) + &
159                         ( gg-cc ) * w(k+1,j+1,i) + &
160                         ( gg-dd ) * w(k+1,j+1,i+1) &
161                        ) / ( 3.0 * gg )
162             w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * &
163                               ( w_int_u - w_int_l )
164          ENDIF
165       ELSE
166          w_int = 0.0
167       ENDIF
168
169!
170!--    Interpolate and calculate quantities needed for calculating the SGS
171!--    velocities
172       IF ( use_sgs_for_particles )  THEN
173!
174!--       Interpolate TKE
175          i = particles(n)%x * ddx
176          j = particles(n)%y * ddy
177          k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
178              + offset_ocean_nzt                      ! only exact if eq.dist
179
180          IF ( topography == 'flat' )  THEN
181
182             = particles(n)%x - i * dx
183             y  = particles(n)%y - j * dy
184             aa = x**2          + y**2
185             bb = ( dx - x )**2 + y**2
186             cc = x**2          + ( dy - y )**2
187             dd = ( dx - x )**2 + ( dy - y )**2
188             gg = aa + bb + cc + dd
189
190             e_int_l = ( ( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1)   &
191                       + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
192                       ) / ( 3.0 * gg )
193
194             IF ( k+1 == nzt+1 )  THEN
195                e_int = e_int_l
196             ELSE
197                e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
198                            ( gg - bb ) * e(k+1,j,i+1) + &
199                            ( gg - cc ) * e(k+1,j+1,i) + &
200                            ( gg - dd ) * e(k+1,j+1,i+1) &
201                          ) / ( 3.0 * gg )
202                e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
203                                  ( e_int_u - e_int_l )
204             ENDIF
205
206!
207!--          Interpolate the TKE gradient along x (adopt incides i,j,k and
208!--          all position variables from above (TKE))
209             de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
210                             ( gg - bb ) * de_dx(k,j,i+1) + &
211                             ( gg - cc ) * de_dx(k,j+1,i) + &
212                             ( gg - dd ) * de_dx(k,j+1,i+1) &
213                           ) / ( 3.0 * gg )
214
215             IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
216                de_dx_int = de_dx_int_l
217             ELSE
218                de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
219                                ( gg - bb ) * de_dx(k+1,j,i+1) + &
220                                ( gg - cc ) * de_dx(k+1,j+1,i) + &
221                                ( gg - dd ) * de_dx(k+1,j+1,i+1) &
222                              ) / ( 3.0 * gg )
223                de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * &
224                                          ( de_dx_int_u - de_dx_int_l )
225             ENDIF
226
227!
228!--          Interpolate the TKE gradient along y
229             de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
230                             ( gg - bb ) * de_dy(k,j,i+1) + &
231                             ( gg - cc ) * de_dy(k,j+1,i) + &
232                             ( gg - dd ) * de_dy(k,j+1,i+1) &
233                           ) / ( 3.0 * gg )
234             IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
235                de_dy_int = de_dy_int_l
236             ELSE
237                de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
238                                ( gg - bb ) * de_dy(k+1,j,i+1) + &
239                                ( gg - cc ) * de_dy(k+1,j+1,i) + &
240                                ( gg - dd ) * de_dy(k+1,j+1,i+1) &
241                              ) / ( 3.0 * gg )
242                de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * &
243                                          ( de_dy_int_u - de_dy_int_l )
244             ENDIF
245
246!
247!--          Interpolate the TKE gradient along z
248             IF ( particles(n)%z < 0.5 * dz ) THEN
249                de_dz_int = 0.0
250             ELSE
251                de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
252                                ( gg - bb ) * de_dz(k,j,i+1) + &
253                                ( gg - cc ) * de_dz(k,j+1,i) + &
254                                ( gg - dd ) * de_dz(k,j+1,i+1) &
255                              ) / ( 3.0 * gg )
256
257                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
258                   de_dz_int = de_dz_int_l
259                ELSE
260                   de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
261                                   ( gg - bb ) * de_dz(k+1,j,i+1) + &
262                                   ( gg - cc ) * de_dz(k+1,j+1,i) + &
263                                   ( gg - dd ) * de_dz(k+1,j+1,i+1) &
264                                 ) / ( 3.0 * gg )
265                   de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * &
266                                             ( de_dz_int_u - de_dz_int_l )
267                ENDIF
268             ENDIF
269
270!
271!--          Interpolate the dissipation of TKE
272             diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
273                            ( gg - bb ) * diss(k,j,i+1) + &
274                            ( gg - cc ) * diss(k,j+1,i) + &
275                            ( gg - dd ) * diss(k,j+1,i+1) &
276                           ) / ( 3.0 * gg )
277
278             IF ( k+1 == nzt+1 )  THEN
279                diss_int = diss_int_l
280             ELSE
281                diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
282                               ( gg - bb ) * diss(k+1,j,i+1) + &
283                               ( gg - cc ) * diss(k+1,j+1,i) + &
284                               ( gg - dd ) * diss(k+1,j+1,i+1) &
285                              ) / ( 3.0 * gg )
286                diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * &
287                                        ( diss_int_u - diss_int_l )
288             ENDIF
289
290          ELSE
291
292!
293!--          In case that there are buildings it has to be determined
294!--          how many of the gridpoints defining the particle box are
295!--          situated within a building
296!--          gp_outside_of_building(1): i,j,k
297!--          gp_outside_of_building(2): i,j+1,k
298!--          gp_outside_of_building(3): i,j,k+1
299!--          gp_outside_of_building(4): i,j+1,k+1
300!--          gp_outside_of_building(5): i+1,j,k
301!--          gp_outside_of_building(6): i+1,j+1,k
302!--          gp_outside_of_building(7): i+1,j,k+1
303!--          gp_outside_of_building(8): i+1,j+1,k+1
304
305             gp_outside_of_building = 0
306             location = 0.0
307             num_gp = 0
308
309             IF ( k > nzb_s_inner(j,i)  .OR.  nzb_s_inner(j,i) == 0 )  THEN
310                num_gp = num_gp + 1
311                gp_outside_of_building(1) = 1
312                location(num_gp,1) = i * dx
313                location(num_gp,2) = j * dy
314                location(num_gp,3) = k * dz - 0.5 * dz
315                ei(num_gp)     = e(k,j,i)
316                dissi(num_gp)  = diss(k,j,i)
317                de_dxi(num_gp) = de_dx(k,j,i)
318                de_dyi(num_gp) = de_dy(k,j,i)
319                de_dzi(num_gp) = de_dz(k,j,i)
320             ENDIF
321
322             IF ( k > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 ) &
323             THEN
324                num_gp = num_gp + 1
325                gp_outside_of_building(2) = 1
326                location(num_gp,1) = i * dx
327                location(num_gp,2) = (j+1) * dy
328                location(num_gp,3) = k * dz - 0.5 * dz
329                ei(num_gp)     = e(k,j+1,i)
330                dissi(num_gp)  = diss(k,j+1,i)
331                de_dxi(num_gp) = de_dx(k,j+1,i)
332                de_dyi(num_gp) = de_dy(k,j+1,i)
333                de_dzi(num_gp) = de_dz(k,j+1,i)
334             ENDIF
335
336             IF ( k+1 > nzb_s_inner(j,i)  .OR.  nzb_s_inner(j,i) == 0 )  THEN
337                num_gp = num_gp + 1
338                gp_outside_of_building(3) = 1
339                location(num_gp,1) = i * dx
340                location(num_gp,2) = j * dy
341                location(num_gp,3) = (k+1) * dz - 0.5 * dz
342                ei(num_gp)     = e(k+1,j,i)
343                dissi(num_gp)  = diss(k+1,j,i)
344                de_dxi(num_gp) = de_dx(k+1,j,i)
345                de_dyi(num_gp) = de_dy(k+1,j,i)
346                de_dzi(num_gp) = de_dz(k+1,j,i)
347             ENDIF
348
349             IF ( k+1 > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 ) &
350             THEN
351                num_gp = num_gp + 1
352                gp_outside_of_building(4) = 1
353                location(num_gp,1) = i * dx
354                location(num_gp,2) = (j+1) * dy
355                location(num_gp,3) = (k+1) * dz - 0.5 * dz
356                ei(num_gp)     = e(k+1,j+1,i)
357                dissi(num_gp)  = diss(k+1,j+1,i)
358                de_dxi(num_gp) = de_dx(k+1,j+1,i)
359                de_dyi(num_gp) = de_dy(k+1,j+1,i)
360                de_dzi(num_gp) = de_dz(k+1,j+1,i)
361             ENDIF
362
363             IF ( k > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 ) &
364             THEN
365                num_gp = num_gp + 1
366                gp_outside_of_building(5) = 1
367                location(num_gp,1) = (i+1) * dx
368                location(num_gp,2) = j * dy
369                location(num_gp,3) = k * dz - 0.5 * dz
370                ei(num_gp)     = e(k,j,i+1)
371                dissi(num_gp)  = diss(k,j,i+1)
372                de_dxi(num_gp) = de_dx(k,j,i+1)
373                de_dyi(num_gp) = de_dy(k,j,i+1)
374                de_dzi(num_gp) = de_dz(k,j,i+1)
375             ENDIF
376
377             IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) &
378             THEN
379                num_gp = num_gp + 1
380                gp_outside_of_building(6) = 1
381                location(num_gp,1) = (i+1) * dx
382                location(num_gp,2) = (j+1) * dy
383                location(num_gp,3) = k * dz - 0.5 * dz
384                ei(num_gp)     = e(k,j+1,i+1)
385                dissi(num_gp)  = diss(k,j+1,i+1)
386                de_dxi(num_gp) = de_dx(k,j+1,i+1)
387                de_dyi(num_gp) = de_dy(k,j+1,i+1)
388                de_dzi(num_gp) = de_dz(k,j+1,i+1)
389             ENDIF
390
391             IF ( k+1 > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 ) &
392             THEN
393                num_gp = num_gp + 1
394                gp_outside_of_building(7) = 1
395                location(num_gp,1) = (i+1) * dx
396                location(num_gp,2) = j * dy
397                location(num_gp,3) = (k+1) * dz - 0.5 * dz
398                ei(num_gp)     = e(k+1,j,i+1)
399                dissi(num_gp)  = diss(k+1,j,i+1)
400                de_dxi(num_gp) = de_dx(k+1,j,i+1)
401                de_dyi(num_gp) = de_dy(k+1,j,i+1)
402                de_dzi(num_gp) = de_dz(k+1,j,i+1)
403             ENDIF
404
405             IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0)&
406             THEN
407                num_gp = num_gp + 1
408                gp_outside_of_building(8) = 1
409                location(num_gp,1) = (i+1) * dx
410                location(num_gp,2) = (j+1) * dy
411                location(num_gp,3) = (k+1) * dz - 0.5 * dz
412                ei(num_gp)     = e(k+1,j+1,i+1)
413                dissi(num_gp)  = diss(k+1,j+1,i+1)
414                de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
415                de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
416                de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
417             ENDIF
418
419!
420!--          If all gridpoints are situated outside of a building, then the
421!--          ordinary interpolation scheme can be used.
422             IF ( num_gp == 8 )  THEN
423
424                = particles(n)%x - i * dx
425                y  = particles(n)%y - j * dy
426                aa = x**2          + y**2
427                bb = ( dx - x )**2 + y**2
428                cc = x**2          + ( dy - y )**2
429                dd = ( dx - x )**2 + ( dy - y )**2
430                gg = aa + bb + cc + dd
431
432                e_int_l = (( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1) &
433                         + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1)&
434                          ) / ( 3.0 * gg )
435
436                IF ( k+1 == nzt+1 )  THEN
437                   e_int = e_int_l
438                ELSE
439                   e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
440                               ( gg - bb ) * e(k+1,j,i+1) + &
441                               ( gg - cc ) * e(k+1,j+1,i) + &
442                               ( gg - dd ) * e(k+1,j+1,i+1) &
443                             ) / ( 3.0 * gg )
444                   e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
445                                     ( e_int_u - e_int_l )
446                ENDIF
447
448!
449!--             Interpolate the TKE gradient along x (adopt incides i,j,k
450!--             and 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 * gg )
456
457                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
458                   de_dx_int = 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 * gg )
465                   de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / &
466                                           dz * ( de_dx_int_u - de_dx_int_l )
467                ENDIF
468
469!
470!--             Interpolate the TKE gradient along y
471                de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
472                                ( gg - bb ) * de_dy(k,j,i+1) + &
473                                ( gg - cc ) * de_dy(k,j+1,i) + &
474                                ( gg - dd ) * de_dy(k,j+1,i+1) &
475                              ) / ( 3.0 * gg )
476                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
477                   de_dy_int = de_dy_int_l
478                ELSE
479                   de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
480                                   ( gg - bb ) * de_dy(k+1,j,i+1) + &
481                                   ( gg - cc ) * de_dy(k+1,j+1,i) + &
482                                   ( gg - dd ) * de_dy(k+1,j+1,i+1) &
483                                 ) / ( 3.0 * gg )
484                   de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / &
485                                           dz * ( de_dy_int_u - de_dy_int_l )
486                ENDIF
487
488!
489!--             Interpolate the TKE gradient along z
490                IF ( particles(n)%z < 0.5 * dz )  THEN
491                   de_dz_int = 0.0
492                ELSE
493                   de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
494                                   ( gg - bb ) * de_dz(k,j,i+1) + &
495                                   ( gg - cc ) * de_dz(k,j+1,i) + &
496                                   ( gg - dd ) * de_dz(k,j+1,i+1) &
497                                 ) / ( 3.0 * gg )
498
499                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
500                      de_dz_int = de_dz_int_l
501                   ELSE
502                      de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
503                                      ( gg - bb ) * de_dz(k+1,j,i+1) + &
504                                      ( gg - cc ) * de_dz(k+1,j+1,i) + &
505                                      ( gg - dd ) * de_dz(k+1,j+1,i+1) &
506                                    ) / ( 3.0 * gg )
507                      de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) /&
508                                           dz * ( de_dz_int_u - de_dz_int_l )
509                   ENDIF
510                ENDIF
511
512!
513!--             Interpolate the dissipation of TKE
514                diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
515                               ( gg - bb ) * diss(k,j,i+1) + &
516                               ( gg - cc ) * diss(k,j+1,i) + &
517                               ( gg - dd ) * diss(k,j+1,i+1) &
518                             ) / ( 3.0 * gg )
519
520                IF ( k+1 == nzt+1 )  THEN
521                   diss_int = diss_int_l
522                ELSE
523                   diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
524                                  ( gg - bb ) * diss(k+1,j,i+1) + &
525                                  ( gg - cc ) * diss(k+1,j+1,i) + &
526                                  ( gg - dd ) * diss(k+1,j+1,i+1) &
527                                ) / ( 3.0 * gg )
528                   diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz *&
529                                           ( diss_int_u - diss_int_l )
530                ENDIF
531
532             ELSE
533
534!
535!--             If wall between gridpoint 1 and gridpoint 5, then
536!--             Neumann boundary condition has to be applied
537                IF ( gp_outside_of_building(1) == 1  .AND. &
538                     gp_outside_of_building(5) == 0 )  THEN
539                   num_gp = num_gp + 1
540                   location(num_gp,1) = i * dx + 0.5 * dx
541                   location(num_gp,2) = j * dy
542                   location(num_gp,3) = k * dz - 0.5 * dz
543                   ei(num_gp)     = e(k,j,i)
544                   dissi(num_gp)  = diss(k,j,i)
545                   de_dxi(num_gp) = 0.0
546                   de_dyi(num_gp) = de_dy(k,j,i)
547                   de_dzi(num_gp) = de_dz(k,j,i)
548                ENDIF
549
550                IF ( gp_outside_of_building(5) == 1  .AND. &
551                   gp_outside_of_building(1) == 0 )  THEN
552                   num_gp = num_gp + 1
553                   location(num_gp,1) = i * dx + 0.5 * dx
554                   location(num_gp,2) = j * dy
555                   location(num_gp,3) = k * dz - 0.5 * dz
556                   ei(num_gp)     = e(k,j,i+1)
557                   dissi(num_gp)  = diss(k,j,i+1)
558                   de_dxi(num_gp) = 0.0
559                   de_dyi(num_gp) = de_dy(k,j,i+1)
560                   de_dzi(num_gp) = de_dz(k,j,i+1)
561                ENDIF
562
563!
564!--             If wall between gridpoint 5 and gridpoint 6, then
565!--             then Neumann boundary condition has to be applied
566                IF ( gp_outside_of_building(5) == 1  .AND. &
567                     gp_outside_of_building(6) == 0 )  THEN
568                   num_gp = num_gp + 1
569                   location(num_gp,1) = (i+1) * dx
570                   location(num_gp,2) = j * dy + 0.5 * dy
571                   location(num_gp,3) = k * dz - 0.5 * dz
572                   ei(num_gp)     = e(k,j,i+1)
573                   dissi(num_gp)  = diss(k,j,i+1)
574                   de_dxi(num_gp) = de_dx(k,j,i+1)
575                   de_dyi(num_gp) = 0.0
576                   de_dzi(num_gp) = de_dz(k,j,i+1)
577                ENDIF
578
579                IF ( gp_outside_of_building(6) == 1  .AND. &
580                     gp_outside_of_building(5) == 0 )  THEN
581                   num_gp = num_gp + 1
582                   location(num_gp,1) = (i+1) * dx
583                   location(num_gp,2) = j * dy + 0.5 * dy
584                   location(num_gp,3) = k * dz - 0.5 * dz
585                   ei(num_gp)     = e(k,j+1,i+1)
586                   dissi(num_gp)  = diss(k,j+1,i+1)
587                   de_dxi(num_gp) = de_dx(k,j+1,i+1)
588                   de_dyi(num_gp) = 0.0
589                   de_dzi(num_gp) = de_dz(k,j+1,i+1)
590                ENDIF
591
592!
593!--             If wall between gridpoint 2 and gridpoint 6, then
594!--             Neumann boundary condition has to be applied
595                IF ( gp_outside_of_building(2) == 1  .AND. &
596                     gp_outside_of_building(6) == 0 )  THEN
597                   num_gp = num_gp + 1
598                   location(num_gp,1) = i * dx + 0.5 * dx
599                   location(num_gp,2) = (j+1) * dy
600                   location(num_gp,3) = k * dz - 0.5 * dz
601                   ei(num_gp)     = e(k,j+1,i)
602                   dissi(num_gp)  = diss(k,j+1,i)
603                   de_dxi(num_gp) = 0.0
604                   de_dyi(num_gp) = de_dy(k,j+1,i)
605                   de_dzi(num_gp) = de_dz(k,j+1,i)
606                ENDIF
607
608                IF ( gp_outside_of_building(6) == 1  .AND. &
609                     gp_outside_of_building(2) == 0 )  THEN
610                   num_gp = num_gp + 1
611                   location(num_gp,1) = i * dx + 0.5 * dx
612                   location(num_gp,2) = (j+1) * dy
613                   location(num_gp,3) = k * dz - 0.5 * dz
614                   ei(num_gp)     = e(k,j+1,i+1)
615                   dissi(num_gp)  = diss(k,j+1,i+1)
616                   de_dxi(num_gp) = 0.0
617                   de_dyi(num_gp) = de_dy(k,j+1,i+1)
618                   de_dzi(num_gp) = de_dz(k,j+1,i+1)
619                ENDIF
620
621!
622!--             If wall between gridpoint 1 and gridpoint 2, then
623!--             Neumann boundary condition has to be applied
624                IF ( gp_outside_of_building(1) == 1  .AND. &
625                     gp_outside_of_building(2) == 0 )  THEN
626                   num_gp = num_gp + 1
627                   location(num_gp,1) = i * dx
628                   location(num_gp,2) = j * dy + 0.5 * dy
629                   location(num_gp,3) = k * dz - 0.5 * dz
630                   ei(num_gp)     = e(k,j,i)
631                   dissi(num_gp)  = diss(k,j,i)
632                   de_dxi(num_gp) = de_dx(k,j,i)
633                   de_dyi(num_gp) = 0.0
634                   de_dzi(num_gp) = de_dz(k,j,i)
635                ENDIF
636
637                IF ( gp_outside_of_building(2) == 1  .AND. &
638                     gp_outside_of_building(1) == 0 )  THEN
639                   num_gp = num_gp + 1
640                   location(num_gp,1) = i * dx
641                   location(num_gp,2) = j * dy + 0.5 * dy
642                   location(num_gp,3) = k * dz - 0.5 * dz
643                   ei(num_gp)     = e(k,j+1,i)
644                   dissi(num_gp)  = diss(k,j+1,i)
645                   de_dxi(num_gp) = de_dx(k,j+1,i)
646                   de_dyi(num_gp) = 0.0
647                   de_dzi(num_gp) = de_dz(k,j+1,i)
648                ENDIF
649
650!
651!--             If wall between gridpoint 3 and gridpoint 7, then
652!--             Neumann boundary condition has to be applied
653                IF ( gp_outside_of_building(3) == 1  .AND. &
654                     gp_outside_of_building(7) == 0 )  THEN
655                   num_gp = num_gp + 1
656                   location(num_gp,1) = i * dx + 0.5 * dx
657                   location(num_gp,2) = j * dy
658                   location(num_gp,3) = k * dz + 0.5 * dz
659                   ei(num_gp)     = e(k+1,j,i)
660                   dissi(num_gp)  = diss(k+1,j,i)
661                   de_dxi(num_gp) = 0.0
662                   de_dyi(num_gp) = de_dy(k+1,j,i)
663                   de_dzi(num_gp) = de_dz(k+1,j,i) 
664                ENDIF
665
666                IF ( gp_outside_of_building(7) == 1  .AND. &
667                     gp_outside_of_building(3) == 0 )  THEN
668                   num_gp = num_gp + 1
669                   location(num_gp,1) = i * dx + 0.5 * dx
670                   location(num_gp,2) = j * dy
671                   location(num_gp,3) = k * dz + 0.5 * dz
672                   ei(num_gp)     = e(k+1,j,i+1)
673                   dissi(num_gp)  = diss(k+1,j,i+1)
674                   de_dxi(num_gp) = 0.0
675                   de_dyi(num_gp) = de_dy(k+1,j,i+1)
676                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
677                ENDIF
678
679!
680!--             If wall between gridpoint 7 and gridpoint 8, then
681!--             Neumann boundary condition has to be applied
682                IF ( gp_outside_of_building(7) == 1  .AND. &
683                     gp_outside_of_building(8) == 0 )  THEN
684                   num_gp = num_gp + 1
685                   location(num_gp,1) = (i+1) * dx
686                   location(num_gp,2) = j * dy + 0.5 * dy
687                   location(num_gp,3) = k * dz + 0.5 * dz
688                   ei(num_gp)     = e(k+1,j,i+1)
689                   dissi(num_gp)  = diss(k+1,j,i+1)
690                   de_dxi(num_gp) = de_dx(k+1,j,i+1)
691                   de_dyi(num_gp) = 0.0
692                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
693                ENDIF
694
695                IF ( gp_outside_of_building(8) == 1  .AND. &
696                     gp_outside_of_building(7) == 0 )  THEN
697                   num_gp = num_gp + 1
698                   location(num_gp,1) = (i+1) * dx
699                   location(num_gp,2) = j * dy + 0.5 * dy
700                   location(num_gp,3) = k * dz + 0.5 * dz
701                   ei(num_gp)     = e(k+1,j+1,i+1)
702                   dissi(num_gp)  = diss(k+1,j+1,i+1)
703                   de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
704                   de_dyi(num_gp) = 0.0
705                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
706                ENDIF
707
708!
709!--             If wall between gridpoint 4 and gridpoint 8, then
710!--             Neumann boundary condition has to be applied
711                IF ( gp_outside_of_building(4) == 1  .AND. &
712                     gp_outside_of_building(8) == 0 )  THEN
713                   num_gp = num_gp + 1
714                   location(num_gp,1) = i * dx + 0.5 * dx
715                   location(num_gp,2) = (j+1) * dy
716                   location(num_gp,3) = k * dz + 0.5 * dz
717                   ei(num_gp)     = e(k+1,j+1,i)
718                   dissi(num_gp)  = diss(k+1,j+1,i)
719                   de_dxi(num_gp) = 0.0
720                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
721                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
722                ENDIF
723
724                IF ( gp_outside_of_building(8) == 1  .AND. &
725                     gp_outside_of_building(4) == 0 )  THEN
726                   num_gp = num_gp + 1
727                   location(num_gp,1) = i * dx + 0.5 * dx
728                   location(num_gp,2) = (j+1) * dy
729                   location(num_gp,3) = k * dz + 0.5 * dz
730                   ei(num_gp)     = e(k+1,j+1,i+1)
731                   dissi(num_gp)  = diss(k+1,j+1,i+1)
732                   de_dxi(num_gp) = 0.0
733                   de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
734                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
735                ENDIF
736
737!
738!--             If wall between gridpoint 3 and gridpoint 4, then
739!--             Neumann boundary condition has to be applied
740                IF ( gp_outside_of_building(3) == 1  .AND. &
741                     gp_outside_of_building(4) == 0 )  THEN
742                   num_gp = num_gp + 1
743                   location(num_gp,1) = i * dx
744                   location(num_gp,2) = j * dy + 0.5 * dy
745                   location(num_gp,3) = k * dz + 0.5 * dz
746                   ei(num_gp)     = e(k+1,j,i)
747                   dissi(num_gp)  = diss(k+1,j,i)
748                   de_dxi(num_gp) = de_dx(k+1,j,i)
749                   de_dyi(num_gp) = 0.0
750                   de_dzi(num_gp) = de_dz(k+1,j,i)
751                ENDIF
752
753                IF ( gp_outside_of_building(4) == 1  .AND. &
754                     gp_outside_of_building(3) == 0 )  THEN
755                   num_gp = num_gp + 1
756                   location(num_gp,1) = i * dx
757                   location(num_gp,2) = j * dy + 0.5 * dy
758                   location(num_gp,3) = k * dz + 0.5 * dz
759                   ei(num_gp)     = e(k+1,j+1,i)
760                   dissi(num_gp)  = diss(k+1,j+1,i)
761                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
762                   de_dyi(num_gp) = 0.0
763                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
764                ENDIF
765
766!
767!--             If wall between gridpoint 1 and gridpoint 3, then
768!--             Neumann boundary condition has to be applied
769!--             (only one case as only building beneath is possible)
770                IF ( gp_outside_of_building(1) == 0  .AND. &
771                     gp_outside_of_building(3) == 1 )  THEN
772                   num_gp = num_gp + 1
773                   location(num_gp,1) = i * dx
774                   location(num_gp,2) = j * dy
775                   location(num_gp,3) = k * dz
776                   ei(num_gp)     = e(k+1,j,i)
777                   dissi(num_gp)  = diss(k+1,j,i)
778                   de_dxi(num_gp) = de_dx(k+1,j,i)
779                   de_dyi(num_gp) = de_dy(k+1,j,i)
780                   de_dzi(num_gp) = 0.0
781                ENDIF
782
783!
784!--             If wall between gridpoint 5 and gridpoint 7, then
785!--             Neumann boundary condition has to be applied
786!--             (only one case as only building beneath is possible)
787                IF ( gp_outside_of_building(5) == 0  .AND. &
788                     gp_outside_of_building(7) == 1 )  THEN
789                   num_gp = num_gp + 1
790                   location(num_gp,1) = (i+1) * dx
791                   location(num_gp,2) = j * dy
792                   location(num_gp,3) = k * dz
793                   ei(num_gp)     = e(k+1,j,i+1)
794                   dissi(num_gp)  = diss(k+1,j,i+1)
795                   de_dxi(num_gp) = de_dx(k+1,j,i+1)
796                   de_dyi(num_gp) = de_dy(k+1,j,i+1)
797                   de_dzi(num_gp) = 0.0
798                ENDIF
799
800!
801!--             If wall between gridpoint 2 and gridpoint 4, then
802!--             Neumann boundary condition has to be applied
803!--             (only one case as only building beneath is possible)
804                IF ( gp_outside_of_building(2) == 0  .AND. &
805                     gp_outside_of_building(4) == 1 )  THEN
806                   num_gp = num_gp + 1
807                   location(num_gp,1) = i * dx
808                   location(num_gp,2) = (j+1) * dy
809                   location(num_gp,3) = k * dz
810                   ei(num_gp)     = e(k+1,j+1,i)
811                   dissi(num_gp)  = diss(k+1,j+1,i)
812                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
813                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
814                   de_dzi(num_gp) = 0.0
815                ENDIF
816
817!
818!--             If wall between gridpoint 6 and gridpoint 8, then
819!--             Neumann boundary condition has to be applied
820!--             (only one case as only building beneath is possible)
821                IF ( gp_outside_of_building(6) == 0  .AND. &
822                     gp_outside_of_building(8) == 1 )  THEN
823                   num_gp = num_gp + 1
824                   location(num_gp,1) = (i+1) * dx
825                   location(num_gp,2) = (j+1) * dy
826                   location(num_gp,3) = k * dz
827                   ei(num_gp)     = e(k+1,j+1,i+1)
828                   dissi(num_gp)  = diss(k+1,j+1,i+1)
829                   de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
830                   de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
831                   de_dzi(num_gp) = 0.0
832                ENDIF
833
834!
835!--             Carry out the interpolation
836                IF ( num_gp == 1 )  THEN
837!
838!--                If only one of the gridpoints is situated outside of the
839!--                building, it follows that the values at the particle
840!--                location are the same as the gridpoint values
841                   e_int     = ei(num_gp)
842                   diss_int  = dissi(num_gp)
843                   de_dx_int = de_dxi(num_gp)
844                   de_dy_int = de_dyi(num_gp)
845                   de_dz_int = de_dzi(num_gp)
846                ELSE IF ( num_gp > 1 )  THEN
847
848                   d_sum = 0.0
849!
850!--                Evaluation of the distances between the gridpoints
851!--                contributing to the interpolated values, and the particle
852!--                location
853                   DO  agp = 1, num_gp
854                      d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2  &
855                                   + ( particles(n)%y-location(agp,2) )**2  &
856                                   + ( particles(n)%z-location(agp,3) )**2
857                      d_sum        = d_sum + d_gp_pl(agp)
858                   ENDDO
859
860!
861!--                Finally the interpolation can be carried out
862                   e_int     = 0.0
863                   diss_int  = 0.0
864                   de_dx_int = 0.0
865                   de_dy_int = 0.0
866                   de_dz_int = 0.0
867                   DO  agp = 1, num_gp
868                      e_int     = e_int     + ( d_sum - d_gp_pl(agp) ) * &
869                                             ei(agp) / ( (num_gp-1) * d_sum )
870                      diss_int  = diss_int  + ( d_sum - d_gp_pl(agp) ) * &
871                                          dissi(agp) / ( (num_gp-1) * d_sum )
872                      de_dx_int = de_dx_int + ( d_sum - d_gp_pl(agp) ) * &
873                                         de_dxi(agp) / ( (num_gp-1) * d_sum )
874                      de_dy_int = de_dy_int + ( d_sum - d_gp_pl(agp) ) * &
875                                         de_dyi(agp) / ( (num_gp-1) * d_sum )
876                      de_dz_int = de_dz_int + ( d_sum - d_gp_pl(agp) ) * &
877                                         de_dzi(agp) / ( (num_gp-1) * d_sum )
878                   ENDDO
879
880                ENDIF
881
882             ENDIF
883
884          ENDIF 
885
886!
887!--       Vertically interpolate the horizontally averaged SGS TKE and
888!--       resolved-scale velocity variances and use the interpolated values
889!--       to calculate the coefficient fs, which is a measure of the ratio
890!--       of the subgrid-scale turbulent kinetic energy to the total amount
891!--       of turbulent kinetic energy.
892          IF ( k == 0 )  THEN
893             e_mean_int = hom(0,1,8,0)
894          ELSE
895             e_mean_int = hom(k,1,8,0) +                                    &
896                                        ( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
897                                        ( zu(k+1) - zu(k) ) *               &
898                                        ( particles(n)%z - zu(k) )
899          ENDIF
900
901          kw = particles(n)%z / dz
902
903          IF ( k == 0 )  THEN
904             aa  = hom(k+1,1,30,0)  * ( particles(n)%z / &
905                                      ( 0.5 * ( zu(k+1) - zu(k) ) ) )
906             bb  = hom(k+1,1,31,0)  * ( particles(n)%z / &
907                                      ( 0.5 * ( zu(k+1) - zu(k) ) ) )
908             cc  = hom(kw+1,1,32,0) * ( particles(n)%z / &
909                                      ( 1.0 * ( zw(kw+1) - zw(kw) ) ) )
910          ELSE
911             aa  = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) * &
912                        ( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) )
913             bb  = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) * &
914                        ( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) )
915             cc  = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) *&
916                        ( ( particles(n)%z - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
917          ENDIF
918
919          vv_int = ( 1.0 / 3.0 ) * ( aa + bb + cc )
920
921          fs_int = ( 2.0 / 3.0 ) * e_mean_int / &
922                      ( vv_int + ( 2.0 / 3.0 ) * e_mean_int )
923
924!
925!--       Calculate the Lagrangian timescale according to Weil et al. (2004).
926          lagr_timescale = ( 4.0 * e_int ) / &
927                           ( 3.0 * fs_int * c_0 * diss_int )
928
929!
930!--       Calculate the next particle timestep. dt_gap is the time needed to
931!--       complete the current LES timestep.
932          dt_gap = dt_3d - particles(n)%dt_sum
933          dt_particle = MIN( dt_3d, 0.025 * lagr_timescale, dt_gap )
934
935!
936!--       The particle timestep should not be too small in order to prevent
937!--       the number of particle timesteps of getting too large
938          IF ( dt_particle < dt_min_part   .AND.  dt_min_part < dt_gap ) &
939          THEN
940             dt_particle = dt_min_part
941          ENDIF
942
943!
944!--       Calculate the SGS velocity components
945          IF ( particles(n)%age == 0.0 )  THEN
946!
947!--          For new particles the SGS components are derived from the SGS
948!--          TKE. Limit the Gaussian random number to the interval
949!--          [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities
950!--          from becoming unrealistically large.
951             particles(n)%rvar1 = SQRT( 2.0 * sgs_wfu_part * e_int ) * &
952                                        ( random_gauss( iran_part, 5.0 ) - 1.0 )
953             particles(n)%rvar2 = SQRT( 2.0 * sgs_wfv_part * e_int ) * &
954                                        ( random_gauss( iran_part, 5.0 ) - 1.0 )
955             particles(n)%rvar3 = SQRT( 2.0 * sgs_wfw_part * e_int ) * &
956                                        ( random_gauss( iran_part, 5.0 ) - 1.0 )
957
958          ELSE
959
960!
961!--          Restriction of the size of the new timestep: compared to the
962!--          previous timestep the increase must not exceed 200%
963
964             dt_particle_m = particles(n)%age - particles(n)%age_m
965             IF ( dt_particle > 2.0 * dt_particle_m )  THEN
966                dt_particle = 2.0 * dt_particle_m
967             ENDIF
968
969!
970!--          For old particles the SGS components are correlated with the
971!--          values from the previous timestep. Random numbers have also to
972!--          be limited (see above).
973!--          As negative values for the subgrid TKE are not allowed, the
974!--          change of the subgrid TKE with time cannot be smaller than
975!--          -e_int/dt_particle. This value is used as a lower boundary
976!--          value for the change of TKE
977
978             de_dt_min = - e_int / dt_particle
979
980             de_dt = ( e_int - particles(n)%e_m ) / dt_particle_m
981
982             IF ( de_dt < de_dt_min )  THEN
983                de_dt = de_dt_min
984             ENDIF
985
986             particles(n)%rvar1 = particles(n)%rvar1 - fs_int * c_0 *          &
987                                  diss_int * particles(n)%rvar1 * dt_particle /&
988                                  ( 4.0 * sgs_wfu_part * e_int ) +             &
989                                  ( 2.0 * sgs_wfu_part * de_dt *               &
990                                    particles(n)%rvar1 /                       &
991                                    ( 2.0 * sgs_wfu_part * e_int ) + de_dx_int &
992                                  ) * dt_particle / 2.0          +             &
993                                  SQRT( fs_int * c_0 * diss_int ) *            &
994                                  ( random_gauss( iran_part, 5.0 ) - 1.0 ) *   &
995                                  SQRT( dt_particle )
996
997             particles(n)%rvar2 = particles(n)%rvar2 - fs_int * c_0 *          &
998                                  diss_int * particles(n)%rvar2 * dt_particle /&
999                                  ( 4.0 * sgs_wfv_part * e_int ) +             &
1000                                  ( 2.0 * sgs_wfv_part * de_dt *               &
1001                                    particles(n)%rvar2 /                       &
1002                                    ( 2.0 * sgs_wfv_part * e_int ) + de_dy_int &
1003                                  ) * dt_particle / 2.0          +             &
1004                                  SQRT( fs_int * c_0 * diss_int ) *            &
1005                                  ( random_gauss( iran_part, 5.0 ) - 1.0 ) *   &
1006                                  SQRT( dt_particle )
1007
1008             particles(n)%rvar3 = particles(n)%rvar3 - fs_int * c_0 *          &
1009                                  diss_int * particles(n)%rvar3 * dt_particle /&
1010                                  ( 4.0 * sgs_wfw_part * e_int ) +             &
1011                                  ( 2.0 * sgs_wfw_part * de_dt *               &
1012                                    particles(n)%rvar3 /                       &
1013                                    ( 2.0 * sgs_wfw_part * e_int ) + de_dz_int &
1014                                  ) * dt_particle / 2.0          +             &
1015                                  SQRT( fs_int * c_0 * diss_int ) *            &
1016                                  ( random_gauss( iran_part, 5.0 ) - 1.0 ) *   &
1017                                  SQRT( dt_particle )
1018
1019          ENDIF
1020
1021          u_int = u_int + particles(n)%rvar1
1022          v_int = v_int + particles(n)%rvar2
1023          w_int = w_int + particles(n)%rvar3
1024
1025!
1026!--       Store the SGS TKE of the current timelevel which is needed for
1027!--       for calculating the SGS particle velocities at the next timestep
1028          particles(n)%e_m = e_int
1029
1030       ELSE
1031!
1032!--       If no SGS velocities are used, only the particle timestep has to
1033!--       be set
1034          dt_particle = dt_3d
1035
1036       ENDIF
1037
1038!
1039!--    Store the old age of the particle ( needed to prevent that a
1040!--    particle crosses several PEs during one timestep, and for the
1041!--    evaluation of the subgrid particle velocity fluctuations )
1042       particles(n)%age_m = particles(n)%age
1043
1044
1045!
1046!--    Particle advection
1047       IF ( particle_groups(particles(n)%group)%density_ratio == 0.0 )  THEN
1048!
1049!--       Pure passive transport (without particle inertia)
1050          particles(n)%x = particles(n)%x + u_int * dt_particle
1051          particles(n)%y = particles(n)%y + v_int * dt_particle
1052          particles(n)%z = particles(n)%z + w_int * dt_particle
1053
1054          particles(n)%speed_x = u_int
1055          particles(n)%speed_y = v_int
1056          particles(n)%speed_z = w_int
1057
1058       ELSE
1059!
1060!--       Transport of particles with inertia
1061          particles(n)%x = particles(n)%x + particles(n)%speed_x * &
1062                                            dt_particle
1063          particles(n)%y = particles(n)%y + particles(n)%speed_y * &
1064                                            dt_particle
1065          particles(n)%z = particles(n)%z + particles(n)%speed_z * &
1066                                            dt_particle
1067
1068!
1069!--       Update of the particle velocity
1070          dens_ratio = particle_groups(particles(n)%group)%density_ratio
1071          IF ( cloud_droplets )  THEN
1072             exp_arg  = 4.5 * dens_ratio * molecular_viscosity /        &
1073                        ( particles(n)%radius )**2 *                    &
1074                        ( 1.0 + 0.15 * ( 2.0 * particles(n)%radius *    &
1075                          SQRT( ( u_int - particles(n)%speed_x )**2 +   &
1076                                ( v_int - particles(n)%speed_y )**2 +   &
1077                                ( w_int - particles(n)%speed_z )**2 ) / &
1078                                         molecular_viscosity )**0.687   &
1079                        )
1080             exp_term = EXP( -exp_arg * dt_particle )
1081          ELSEIF ( use_sgs_for_particles )  THEN
1082             exp_arg  = particle_groups(particles(n)%group)%exp_arg
1083             exp_term = EXP( -exp_arg * dt_particle )
1084          ELSE
1085             exp_arg  = particle_groups(particles(n)%group)%exp_arg
1086             exp_term = particle_groups(particles(n)%group)%exp_term
1087          ENDIF
1088          particles(n)%speed_x = particles(n)%speed_x * exp_term +             &
1089                                 u_int * ( 1.0 - exp_term )
1090          particles(n)%speed_y = particles(n)%speed_y * exp_term +             &
1091                                 v_int * ( 1.0 - exp_term )
1092          particles(n)%speed_z = particles(n)%speed_z * exp_term +             &
1093                                 ( w_int - ( 1.0 - dens_ratio ) * g / exp_arg )&
1094                                 * ( 1.0 - exp_term )
1095       ENDIF
1096
1097!
1098!--    Increment the particle age and the total time that the particle
1099!--    has advanced within the particle timestep procedure
1100       particles(n)%age    = particles(n)%age    + dt_particle
1101       particles(n)%dt_sum = particles(n)%dt_sum + dt_particle
1102
1103!
1104!--    Check whether there is still a particle that has not yet completed
1105!--    the total LES timestep
1106       IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8 )  THEN
1107          dt_3d_reached_l = .FALSE.
1108       ENDIF
1109
1110    ENDDO
1111
1112
1113 END SUBROUTINE lpm_advec
Note: See TracBrowser for help on using the repository browser.