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

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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