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

Last change on this file since 1322 was 1322, checked in by raasch, 8 years ago

REAL functions and a lot of REAL constants provided with KIND-attribute,
some small bugfixes

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