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

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

last commit documented

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