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

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

last commit documented

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