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

Last change on this file since 849 was 849, checked in by raasch, 10 years ago

Changed:


Original routine advec_particles split into several new subroutines and renamed
lpm.
init_particles renamed lpm_init
user_advec_particles renamed user_lpm_advec,
particle_boundary_conds renamed lpm_boundary_conds,
set_particle_attributes renamed lpm_set_attributes,
user_init_particles renamed user_lpm_init,
user_particle_attributes renamed user_lpm_set_attributes
(Makefile, lpm_droplet_collision, lpm_droplet_condensation, init_3d_model, modules, palm, read_var_list, time_integration, write_var_list, deleted: advec_particles, init_particles, particle_boundary_conds, set_particle_attributes, user_advec_particles, user_init_particles, user_particle_attributes, new: lpm, lpm_advec, lpm_boundary_conds, lpm_calc_liquid_water_content, lpm_data_output_particles, lpm_droplet_collision, lpm_drollet_condensation, lpm_exchange_horiz, lpm_extend_particle_array, lpm_extend_tails, lpm_extend_tail_array, lpm_init, lpm_init_sgs_tke, lpm_pack_arrays, lpm_read_restart_file, lpm_release_set, lpm_set_attributes, lpm_sort_arrays, lpm_write_exchange_statistics, lpm_write_restart_file, user_lpm_advec, user_lpm_init, user_lpm_set_attributes

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