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

Last change on this file since 1315 was 1315, checked in by suehring, 10 years ago

last commit documented

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