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

Last change on this file since 1036 was 1036, checked in by raasch, 9 years ago

code has been put under the GNU General Public License (v3)

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