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

Last change on this file since 1314 was 1314, checked in by suehring, 8 years ago

Vertical logarithmic interpolation of horizontal particle speed for particles
between roughness height and first vertical grid level.

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