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

Last change on this file since 1350 was 1323, checked in by raasch, 11 years ago

last commit documented

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