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

Last change on this file since 1930 was 1930, checked in by suehring, 6 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 72.6 KB
Line 
1!> @file lpm_advec.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
21
22!
23! Former revisions:
24! -----------------
25! $Id: lpm_advec.f90 1930 2016-06-09 16:32:12Z suehring $
26!
27! 1929 2016-06-09 16:25:25Z suehring
28! Put stochastic equation in an extra subroutine.
29! Set flag for stochastic equation to communicate whether a particle is near
30! topography. This case, memory and drift term are disabled in the Weil equation.
31!
32! Enable vertical logarithmic interpolation also above topography. This case,
33! set a lower limit for the friction velocity, as it can become very small
34! in narrow street canyons, leading to too large particle speeds.
35!
36! 1888 2016-04-21 12:20:49Z suehring
37! Bugfix concerning logarithmic interpolation of particle speed
38!
39! 1822 2016-04-07 07:49:42Z hoffmann
40! Random velocity fluctuations for particles added. Terminal fall velocity
41! for droplets is calculated from a parameterization (which is better than
42! the previous, physically correct calculation, which demands a very short
43! time step that is not used in the model).
44!
45! Unused variables deleted.
46!
47! 1691 2015-10-26 16:17:44Z maronga
48! Renamed prandtl_layer to constant_flux_layer.
49!
50! 1685 2015-10-08 07:32:13Z raasch
51! TKE check for negative values (so far, only zero value was checked)
52! offset_ocean_nzt_m1 removed
53!
54! 1682 2015-10-07 23:56:08Z knoop
55! Code annotations made doxygen readable
56!
57! 1583 2015-04-15 12:16:27Z suehring
58! Bugfix: particle advection within Prandtl-layer in case of Galilei
59! transformation.
60!
61! 1369 2014-04-24 05:57:38Z raasch
62! usage of module interfaces removed
63!
64! 1359 2014-04-11 17:15:14Z hoffmann
65! New particle structure integrated.
66! Kind definition added to all floating point numbers.
67!
68! 1322 2014-03-20 16:38:49Z raasch
69! REAL constants defined as wp_kind
70!
71! 1320 2014-03-20 08:40:49Z raasch
72! ONLY-attribute added to USE-statements,
73! kind-parameters added to all INTEGER and REAL declaration statements,
74! kinds are defined in new module kinds,
75! revision history before 2012 removed,
76! comment fields (!:) to be used for variable explanations added to
77! all variable declaration statements
78!
79! 1314 2014-03-14 18:25:17Z suehring
80! Vertical logarithmic interpolation of horizontal particle speed for particles
81! between roughness height and first vertical grid level.
82!
83! 1036 2012-10-22 13:43:42Z raasch
84! code put under GPL (PALM 3.9)
85!
86! 849 2012-03-15 10:35:09Z raasch
87! initial revision (former part of advec_particles)
88!
89!
90! Description:
91! ------------
92!> Calculation of new particle positions due to advection using a simple Euler
93!> scheme. Particles may feel inertia effects. SGS transport can be included
94!> using the stochastic model of Weil et al. (2004, JAS, 61, 2877-2887).
95!------------------------------------------------------------------------------!
96 SUBROUTINE lpm_advec (ip,jp,kp)
97 
98
99    USE arrays_3d,                                                             &
100        ONLY:  de_dx, de_dy, de_dz, diss, e, km, u, us, usws, v, vsws, w, zu, zw
101
102    USE cpulog
103
104    USE pegrid
105
106    USE control_parameters,                                                    &
107        ONLY:  atmos_ocean_sign, cloud_droplets, constant_flux_layer, dt_3d,   &
108               dt_3d_reached_l, dz, g, kappa, topography, u_gtrans, v_gtrans
109
110    USE grid_variables,                                                        &
111        ONLY:  ddx, dx, ddy, dy
112       
113    USE indices,                                                               &
114        ONLY:  nzb, nzb_s_inner, nzt
115       
116    USE kinds
117   
118    USE particle_attributes,                                                   &
119        ONLY:  block_offset, c_0, dt_min_part, grid_particles,                 &
120               iran_part, log_z_z0, number_of_particles, number_of_sublayers,  &
121               particles, particle_groups, offset_ocean_nzt, sgs_wf_part,      &
122               use_sgs_for_particles, vertical_particle_advection, z0_av_global
123       
124    USE statistics,                                                            &
125        ONLY:  hom
126
127    IMPLICIT NONE
128
129    INTEGER(iwp) ::  agp                         !< loop variable
130    INTEGER(iwp) ::  gp_outside_of_building(1:8) !< number of grid points used for particle interpolation in case of topography
131    INTEGER(iwp) ::  i                           !< index variable along x
132    INTEGER(iwp) ::  ip                          !< index variable along x
133    INTEGER(iwp) ::  ilog                        !< index variable along x
134    INTEGER(iwp) ::  j                           !< index variable along y
135    INTEGER(iwp) ::  jp                          !< index variable along y
136    INTEGER(iwp) ::  jlog                        !< index variable along y
137    INTEGER(iwp) ::  k                           !< index variable along z
138    INTEGER(iwp) ::  kp                          !< index variable along z
139    INTEGER(iwp) ::  kw                          !< index variable along z
140    INTEGER(iwp) ::  n                           !< loop variable over all particles in a grid box
141    INTEGER(iwp) ::  nb                          !< block number particles are sorted in
142    INTEGER(iwp) ::  num_gp                      !< number of adjacent grid points inside topography
143
144    INTEGER(iwp), DIMENSION(0:7) ::  start_index !< start particle index for current block
145    INTEGER(iwp), DIMENSION(0:7) ::  end_index   !< start particle index for current block
146
147    REAL(wp) ::  aa                 !< dummy argument for horizontal particle interpolation
148    REAL(wp) ::  bb                 !< dummy argument for horizontal particle interpolation
149    REAL(wp) ::  cc                 !< dummy argument for horizontal particle interpolation
150    REAL(wp) ::  d_sum              !< dummy argument for horizontal particle interpolation in case of topography
151    REAL(wp) ::  d_z_p_z0           !< inverse of interpolation length for logarithmic interpolation
152    REAL(wp) ::  dd                 !< dummy argument for horizontal particle interpolation
153    REAL(wp) ::  de_dx_int_l        !< x/y-interpolated TKE gradient (x) at particle position at lower vertical level
154    REAL(wp) ::  de_dx_int_u        !< x/y-interpolated TKE gradient (x) at particle position at upper vertical level
155    REAL(wp) ::  de_dy_int_l        !< x/y-interpolated TKE gradient (y) at particle position at lower vertical level
156    REAL(wp) ::  de_dy_int_u        !< x/y-interpolated TKE gradient (y) at particle position at upper vertical level
157    REAL(wp) ::  de_dt              !< temporal derivative of TKE experienced by the particle
158    REAL(wp) ::  de_dt_min          !< lower level for temporal TKE derivative
159    REAL(wp) ::  de_dz_int_l        !< x/y-interpolated TKE gradient (z) at particle position at lower vertical level
160    REAL(wp) ::  de_dz_int_u        !< x/y-interpolated TKE gradient (z) at particle position at upper vertical level
161    REAL(wp) ::  diameter           !< diamter of droplet
162    REAL(wp) ::  diss_int_l         !< x/y-interpolated dissipation at particle position at lower vertical level
163    REAL(wp) ::  diss_int_u         !< x/y-interpolated dissipation at particle position at upper vertical level
164    REAL(wp) ::  dt_gap             !< remaining time until particle time integration reaches LES time
165    REAL(wp) ::  dt_particle_m      !< previous particle time step
166    REAL(wp) ::  e_int_l            !< x/y-interpolated TKE at particle position at lower vertical level
167    REAL(wp) ::  e_int_u            !< x/y-interpolated TKE at particle position at upper vertical level
168    REAL(wp) ::  e_mean_int         !< horizontal mean TKE at particle height
169    REAL(wp) ::  exp_arg            !<
170    REAL(wp) ::  exp_term           !<
171    REAL(wp) ::  gg                 !< dummy argument for horizontal particle interpolation
172    REAL(wp) ::  height_p           !< dummy argument for logarithmic interpolation
173    REAL(wp) ::  lagr_timescale     !< Lagrangian timescale
174    REAL(wp) ::  location(1:30,1:3) !< wall locations
175    REAL(wp) ::  log_z_z0_int       !< logarithmus used for surface_layer interpolation
176    REAL(wp) ::  random_gauss       !<
177    REAL(wp) ::  RL                 !< Lagrangian autocorrelation coefficient
178    REAL(wp) ::  rg1                !< Gaussian distributed random number
179    REAL(wp) ::  rg2                !< Gaussian distributed random number
180    REAL(wp) ::  rg3                !< Gaussian distributed random number
181    REAL(wp) ::  sigma              !< velocity standard deviation
182    REAL(wp) ::  u_int_l            !< x/y-interpolated u-component at particle position at lower vertical level
183    REAL(wp) ::  u_int_u            !< x/y-interpolated u-component at particle position at upper vertical level
184    REAL(wp) ::  us_int             !< friction velocity at particle grid box
185    REAL(wp) ::  v_int_l            !< x/y-interpolated v-component at particle position at lower vertical level
186    REAL(wp) ::  v_int_u            !< x/y-interpolated v-component at particle position at upper vertical level
187    REAL(wp) ::  vv_int             !<
188    REAL(wp) ::  w_int_l            !< x/y-interpolated w-component at particle position at lower vertical level
189    REAL(wp) ::  w_int_u            !< x/y-interpolated w-component at particle position at upper vertical level
190    REAL(wp) ::  w_s                !< terminal velocity of droplets
191    REAL(wp) ::  x                  !< dummy argument for horizontal particle interpolation
192    REAL(wp) ::  y                  !< dummy argument for horizontal particle interpolation
193    REAL(wp) ::  z_p                !< surface layer height (0.5 dz)
194
195    REAL(wp), PARAMETER ::  a_rog = 9.65_wp      !< parameter for fall velocity
196    REAL(wp), PARAMETER ::  b_rog = 10.43_wp     !< parameter for fall velocity
197    REAL(wp), PARAMETER ::  c_rog = 0.6_wp       !< parameter for fall velocity
198    REAL(wp), PARAMETER ::  k_cap_rog = 4.0_wp   !< parameter for fall velocity
199    REAL(wp), PARAMETER ::  k_low_rog = 12.0_wp  !< parameter for fall velocity
200    REAL(wp), PARAMETER ::  d0_rog = 0.745_wp    !< separation diameter
201
202    REAL(wp), DIMENSION(1:30) ::  d_gp_pl !< dummy argument for particle interpolation scheme in case of topography
203    REAL(wp), DIMENSION(1:30) ::  de_dxi  !< horizontal TKE gradient along x at adjacent wall
204    REAL(wp), DIMENSION(1:30) ::  de_dyi  !< horizontal TKE gradient along y at adjacent wall
205    REAL(wp), DIMENSION(1:30) ::  de_dzi  !< horizontal TKE gradient along z at adjacent wall
206    REAL(wp), DIMENSION(1:30) ::  dissi   !< dissipation at adjacent wall
207    REAL(wp), DIMENSION(1:30) ::  ei      !< TKE at adjacent wall
208
209    REAL(wp), DIMENSION(number_of_particles) ::  term_1_2     !< flag to communicate whether a particle is near topography or not
210    REAL(wp), DIMENSION(number_of_particles) ::  dens_ratio   !<
211    REAL(wp), DIMENSION(number_of_particles) ::  de_dx_int    !< horizontal TKE gradient along x at particle position
212    REAL(wp), DIMENSION(number_of_particles) ::  de_dy_int    !< horizontal TKE gradient along y at particle position
213    REAL(wp), DIMENSION(number_of_particles) ::  de_dz_int    !< horizontal TKE gradient along z at particle position
214    REAL(wp), DIMENSION(number_of_particles) ::  diss_int     !< dissipation at particle position
215    REAL(wp), DIMENSION(number_of_particles) ::  dt_particle  !< particle time step
216    REAL(wp), DIMENSION(number_of_particles) ::  e_int        !< TKE at particle position
217    REAL(wp), DIMENSION(number_of_particles) ::  fs_int       !< weighting factor for subgrid-scale particle speed
218    REAL(wp), DIMENSION(number_of_particles) ::  u_int        !< u-component of particle speed
219    REAL(wp), DIMENSION(number_of_particles) ::  v_int        !< v-component of particle speed
220    REAL(wp), DIMENSION(number_of_particles) ::  w_int        !< w-component of particle speed
221    REAL(wp), DIMENSION(number_of_particles) ::  xv           !< x-position
222    REAL(wp), DIMENSION(number_of_particles) ::  yv           !< y-position
223    REAL(wp), DIMENSION(number_of_particles) ::  zv           !< z-position
224
225    REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !< vector of Gaussian distributed random numbers
226
227    CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' )
228
229!
230!-- Determine height of Prandtl layer and distance between Prandtl-layer
231!-- height and horizontal mean roughness height, which are required for
232!-- vertical logarithmic interpolation of horizontal particle speeds
233!-- (for particles below first vertical grid level).
234    z_p      = zu(nzb+1) - zw(nzb)
235    d_z_p_z0 = 1.0_wp / ( z_p - z0_av_global )
236
237    start_index = grid_particles(kp,jp,ip)%start_index
238    end_index   = grid_particles(kp,jp,ip)%end_index
239
240    xv = particles(1:number_of_particles)%x
241    yv = particles(1:number_of_particles)%y
242    zv = particles(1:number_of_particles)%z
243
244    DO  nb = 0, 7
245
246       i = ip
247       j = jp + block_offset(nb)%j_off
248       k = kp + block_offset(nb)%k_off
249
250
251!
252!--    Interpolate u velocity-component
253       DO  n = start_index(nb), end_index(nb)
254!
255!--       Interpolation of the u velocity component onto particle position. 
256!--       Particles are interpolation bi-linearly in the horizontal and a
257!--       linearly in the vertical. An exception is made for particles below
258!--       the first vertical grid level in case of a prandtl layer. In this
259!--       case the horizontal particle velocity components are determined using
260!--       Monin-Obukhov relations (if branch).
261!--       First, check if particle is located below first vertical grid level
262!--       (Prandtl-layer height)
263          ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx
264          jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy
265
266          IF ( constant_flux_layer .AND. zv(n) - zw(nzb_s_inner(jlog,ilog)) < z_p )  THEN
267!
268!--          Resolved-scale horizontal particle velocity is zero below z0.
269             IF ( zv(n) - zw(nzb_s_inner(jlog,ilog)) < z0_av_global )  THEN
270                u_int(n) = 0.0_wp
271             ELSE
272!
273!--             Determine the sublayer. Further used as index.
274                height_p     = ( zv(n) - zw(nzb_s_inner(jlog,ilog)) - z0_av_global )            &
275                                     * REAL( number_of_sublayers, KIND=wp ) &
276                                     * d_z_p_z0
277!
278!--             Calculate LOG(z/z0) for exact particle height. Therefore,   
279!--             interpolate linearly between precalculated logarithm.
280                log_z_z0_int = log_z_z0(INT(height_p))                         &
281                                 + ( height_p - INT(height_p) )                &
282                                 * ( log_z_z0(INT(height_p)+1)                 &
283                                      - log_z_z0(INT(height_p))                &
284                                   ) 
285!
286!--             Limit friction velocity. In narrow canyons or holes the
287!--             friction velocity can become very small, resulting in a too
288!--             large particle speed.
289                us_int   = MAX( 0.5_wp * ( us(jlog,ilog) + us(jlog,ilog-1) ),  &
290                                0.01_wp ) 
291!
292!--             Neutral solution is applied for all situations, e.g. also for
293!--             unstable and stable situations. Even though this is not exact
294!--             this saves a lot of CPU time since several calls of intrinsic
295!--             FORTRAN procedures (LOG, ATAN) are avoided, This is justified
296!--             as sensitivity studies revealed no significant effect of
297!--             using the neutral solution also for un/stable situations.
298                u_int(n) = -usws(jlog,ilog) / ( us_int * kappa + 1E-10_wp )          & 
299                            * log_z_z0_int - u_gtrans
300               
301             ENDIF
302!
303!--       Particle above the first grid level. Bi-linear interpolation in the
304!--       horizontal and linear interpolation in the vertical direction.
305          ELSE
306
307             = xv(n) + ( 0.5_wp - i ) * dx
308             y  = yv(n) - j * dy
309             aa = x**2          + y**2
310             bb = ( dx - x )**2 + y**2
311             cc = x**2          + ( dy - y )**2
312             dd = ( dx - x )**2 + ( dy - y )**2
313             gg = aa + bb + cc + dd
314
315             u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
316                         + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) *            &
317                         u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
318
319             IF ( k == nzt )  THEN
320                u_int(n) = u_int_l
321             ELSE
322                u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1)  &
323                            + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) *           &
324                            u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
325                u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dz *                  &
326                           ( u_int_u - u_int_l )
327             ENDIF
328
329          ENDIF
330
331       ENDDO
332
333       i = ip + block_offset(nb)%i_off
334       j = jp
335       k = kp + block_offset(nb)%k_off
336!
337!--    Same procedure for interpolation of the v velocity-component
338       DO  n = start_index(nb), end_index(nb)
339
340          ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx
341          jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy
342          IF ( constant_flux_layer .AND. zv(n) - zw(nzb_s_inner(jlog,ilog)) < z_p )  THEN
343
344             IF ( zv(n) - zw(nzb_s_inner(jlog,ilog)) < z0_av_global )  THEN
345!
346!--             Resolved-scale horizontal particle velocity is zero below z0.
347                v_int(n) = 0.0_wp
348             ELSE       
349!
350!--             Determine the sublayer. Further used as index. Please note,
351!--             logarithmus can not be reused from above, as in in case of
352!--             topography particle on u-grid can be above surface-layer height,
353!--             whereas it can be below on v-grid.
354                height_p     = ( zv(n) - zw(nzb_s_inner(jlog,ilog)) - z0_av_global ) &
355                                  * REAL( number_of_sublayers, KIND=wp )    &
356                                  * d_z_p_z0
357!
358!--             Calculate LOG(z/z0) for exact particle height. Therefore,   
359!--             interpolate linearly between precalculated logarithm.
360                log_z_z0_int = log_z_z0(INT(height_p))                         &
361                                 + ( height_p - INT(height_p) )                &
362                                 * ( log_z_z0(INT(height_p)+1)                 &
363                                      - log_z_z0(INT(height_p))                &
364                                   ) 
365!
366!--             Limit friction velocity. In narrow canyons or holes the
367!--             friction velocity can become very small, resulting in a too
368!--             large particle speed.
369                us_int   = MAX( 0.5_wp * ( us(jlog,ilog) + us(jlog-1,ilog) ), &
370                                0.01_wp )   
371!
372!--             Neutral solution is applied for all situations, e.g. also for
373!--             unstable and stable situations. Even though this is not exact
374!--             this saves a lot of CPU time since several calls of intrinsic
375!--             FORTRAN procedures (LOG, ATAN) are avoided, This is justified
376!--             as sensitivity studies revealed no significant effect of
377!--             using the neutral solution also for un/stable situations.
378                v_int(n) = -vsws(jlog,ilog) / ( us_int * kappa + 1E-10_wp )          &
379                         * log_z_z0_int - v_gtrans
380
381             ENDIF
382
383          ELSE
384             = xv(n) - i * dx
385             y  = yv(n) + ( 0.5_wp - j ) * dy
386             aa = x**2          + y**2
387             bb = ( dx - x )**2 + y**2
388             cc = x**2          + ( dy - y )**2
389             dd = ( dx - x )**2 + ( dy - y )**2
390             gg = aa + bb + cc + dd
391
392             v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
393                       + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
394                       ) / ( 3.0_wp * gg ) - v_gtrans
395
396             IF ( k == nzt )  THEN
397                v_int(n) = v_int_l
398             ELSE
399                v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
400                          + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
401                          ) / ( 3.0_wp * gg ) - v_gtrans
402                v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dz * &
403                                  ( v_int_u - v_int_l )
404             ENDIF
405
406          ENDIF
407
408       ENDDO
409
410       i = ip + block_offset(nb)%i_off
411       j = jp + block_offset(nb)%j_off
412       k = kp - 1
413!
414!--    Same procedure for interpolation of the w velocity-component
415       DO  n = start_index(nb), end_index(nb)
416
417          IF ( vertical_particle_advection(particles(n)%group) )  THEN
418
419             = xv(n) - i * dx
420             y  = yv(n) - j * dy
421             aa = x**2          + y**2
422             bb = ( dx - x )**2 + y**2
423             cc = x**2          + ( dy - y )**2
424             dd = ( dx - x )**2 + ( dy - y )**2
425             gg = aa + bb + cc + dd
426
427             w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
428                       + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
429                       ) / ( 3.0_wp * gg )
430
431             IF ( k == nzt )  THEN
432                w_int(n) = w_int_l
433             ELSE
434                w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + &
435                            ( gg-bb ) * w(k+1,j,i+1) + &
436                            ( gg-cc ) * w(k+1,j+1,i) + &
437                            ( gg-dd ) * w(k+1,j+1,i+1) &
438                          ) / ( 3.0_wp * gg )
439                w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dz * &
440                           ( w_int_u - w_int_l )
441             ENDIF
442
443          ELSE
444
445             w_int(n) = 0.0_wp
446
447          ENDIF
448
449       ENDDO
450
451    ENDDO
452
453!-- Interpolate and calculate quantities needed for calculating the SGS
454!-- velocities
455    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
456
457       IF ( topography == 'flat' )  THEN
458
459          DO  nb = 0,7
460
461             i = ip + block_offset(nb)%i_off
462             j = jp + block_offset(nb)%j_off
463             k = kp + block_offset(nb)%k_off
464
465             DO  n = start_index(nb), end_index(nb)
466!
467!--             Interpolate TKE
468                x  = xv(n) - i * dx
469                y  = yv(n) - j * dy
470                aa = x**2          + y**2
471                bb = ( dx - x )**2 + y**2
472                cc = x**2          + ( dy - y )**2
473                dd = ( dx - x )**2 + ( dy - y )**2
474                gg = aa + bb + cc + dd
475
476                e_int_l = ( ( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1)   &
477                          + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
478                          ) / ( 3.0_wp * gg )
479
480                IF ( k+1 == nzt+1 )  THEN
481                   e_int(n) = e_int_l
482                ELSE
483                   e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
484                               ( gg - bb ) * e(k+1,j,i+1) + &
485                               ( gg - cc ) * e(k+1,j+1,i) + &
486                               ( gg - dd ) * e(k+1,j+1,i+1) &
487                            ) / ( 3.0_wp * gg )
488                   e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
489                                     ( e_int_u - e_int_l )
490                ENDIF
491!
492!--             Needed to avoid NaN particle velocities (this might not be
493!--             required any more)
494                IF ( e_int(n) <= 0.0_wp )  THEN
495                   e_int(n) = 1.0E-20_wp
496                ENDIF
497!
498!--             Interpolate the TKE gradient along x (adopt incides i,j,k and
499!--             all position variables from above (TKE))
500                de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
501                                ( gg - bb ) * de_dx(k,j,i+1) + &
502                                ( gg - cc ) * de_dx(k,j+1,i) + &
503                                ( gg - dd ) * de_dx(k,j+1,i+1) &
504                               ) / ( 3.0_wp * gg )
505
506                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
507                   de_dx_int(n) = de_dx_int_l
508                ELSE
509                   de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
510                                   ( gg - bb ) * de_dx(k+1,j,i+1) + &
511                                   ( gg - cc ) * de_dx(k+1,j+1,i) + &
512                                   ( gg - dd ) * de_dx(k+1,j+1,i+1) &
513                                  ) / ( 3.0_wp * gg )
514                   de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / dz * &
515                                              ( de_dx_int_u - de_dx_int_l )
516                ENDIF
517!
518!--             Interpolate the TKE gradient along y
519                de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
520                                ( gg - bb ) * de_dy(k,j,i+1) + &
521                                ( gg - cc ) * de_dy(k,j+1,i) + &
522                                ( gg - dd ) * de_dy(k,j+1,i+1) &
523                               ) / ( 3.0_wp * gg )
524                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
525                   de_dy_int(n) = de_dy_int_l
526                ELSE
527                   de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
528                                  ( gg - bb ) * de_dy(k+1,j,i+1) + &
529                                  ( gg - cc ) * de_dy(k+1,j+1,i) + &
530                                  ( gg - dd ) * de_dy(k+1,j+1,i+1) &
531                                  ) / ( 3.0_wp * gg )
532                   de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * &
533                                              ( de_dy_int_u - de_dy_int_l )
534                ENDIF
535
536!
537!--             Interpolate the TKE gradient along z
538                IF ( zv(n) < 0.5_wp * dz )  THEN
539                   de_dz_int(n) = 0.0_wp
540                ELSE
541                   de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
542                                   ( gg - bb ) * de_dz(k,j,i+1) + &
543                                   ( gg - cc ) * de_dz(k,j+1,i) + &
544                                   ( gg - dd ) * de_dz(k,j+1,i+1) &
545                                  ) / ( 3.0_wp * gg )
546
547                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
548                      de_dz_int(n) = de_dz_int_l
549                   ELSE
550                      de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
551                                      ( gg - bb ) * de_dz(k+1,j,i+1) + &
552                                      ( gg - cc ) * de_dz(k+1,j+1,i) + &
553                                      ( gg - dd ) * de_dz(k+1,j+1,i+1) &
554                                     ) / ( 3.0_wp * gg )
555                      de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) / dz * &
556                                                 ( de_dz_int_u - de_dz_int_l )
557                   ENDIF
558                ENDIF
559
560!
561!--             Interpolate the dissipation of TKE
562                diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
563                               ( gg - bb ) * diss(k,j,i+1) + &
564                               ( gg - cc ) * diss(k,j+1,i) + &
565                               ( gg - dd ) * diss(k,j+1,i+1) &
566                              ) / ( 3.0_wp * gg )
567
568                IF ( k == nzt )  THEN
569                   diss_int(n) = diss_int_l
570                ELSE
571                   diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
572                                  ( gg - bb ) * diss(k+1,j,i+1) + &
573                                  ( gg - cc ) * diss(k+1,j+1,i) + &
574                                  ( gg - dd ) * diss(k+1,j+1,i+1) &
575                                 ) / ( 3.0_wp * gg )
576                   diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz * &
577                                           ( diss_int_u - diss_int_l )
578                ENDIF
579
580!
581!--             Set flag for stochastic equation.
582                term_1_2(n) = 1.0_wp
583
584             ENDDO
585          ENDDO
586
587       ELSE ! non-flat topography, e.g., buildings
588
589          DO  n = 1, number_of_particles
590             i = particles(n)%x * ddx
591             j = particles(n)%y * ddy
592             k = ( zv(n) + 0.5_wp * dz * atmos_ocean_sign ) / dz  &
593                 + offset_ocean_nzt                      ! only exact if eq.dist
594!
595!--          In case that there are buildings it has to be determined
596!--          how many of the gridpoints defining the particle box are
597!--          situated within a building
598!--          gp_outside_of_building(1): i,j,k
599!--          gp_outside_of_building(2): i,j+1,k
600!--          gp_outside_of_building(3): i,j,k+1
601!--          gp_outside_of_building(4): i,j+1,k+1
602!--          gp_outside_of_building(5): i+1,j,k
603!--          gp_outside_of_building(6): i+1,j+1,k
604!--          gp_outside_of_building(7): i+1,j,k+1
605!--          gp_outside_of_building(8): i+1,j+1,k+1
606
607             gp_outside_of_building = 0
608             location = 0.0_wp
609             num_gp = 0
610
611             IF ( k > nzb_s_inner(j,i)  .OR.  nzb_s_inner(j,i) == 0 )  THEN
612                num_gp = num_gp + 1
613                gp_outside_of_building(1) = 1
614                location(num_gp,1) = i * dx
615                location(num_gp,2) = j * dy
616                location(num_gp,3) = k * dz - 0.5_wp * dz
617                ei(num_gp)     = e(k,j,i)
618                dissi(num_gp)  = diss(k,j,i)
619                de_dxi(num_gp) = de_dx(k,j,i)
620                de_dyi(num_gp) = de_dy(k,j,i)
621                de_dzi(num_gp) = de_dz(k,j,i)
622             ENDIF
623             IF ( k > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 )  THEN
624                num_gp = num_gp + 1
625                gp_outside_of_building(2) = 1
626                location(num_gp,1) = i * dx
627                location(num_gp,2) = (j+1) * dy
628                location(num_gp,3) = k * dz - 0.5_wp * dz
629                ei(num_gp)     = e(k,j+1,i)
630                dissi(num_gp)  = diss(k,j+1,i)
631                de_dxi(num_gp) = de_dx(k,j+1,i)
632                de_dyi(num_gp) = de_dy(k,j+1,i)
633                de_dzi(num_gp) = de_dz(k,j+1,i)
634             ENDIF
635
636             IF ( k+1 > nzb_s_inner(j,i)  .OR.  nzb_s_inner(j,i) == 0 )  THEN
637                num_gp = num_gp + 1
638                gp_outside_of_building(3) = 1
639                location(num_gp,1) = i * dx
640                location(num_gp,2) = j * dy
641                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
642                ei(num_gp)     = e(k+1,j,i)
643                dissi(num_gp)  = diss(k+1,j,i)
644                de_dxi(num_gp) = de_dx(k+1,j,i)
645                de_dyi(num_gp) = de_dy(k+1,j,i)
646                de_dzi(num_gp) = de_dz(k+1,j,i)
647             ENDIF
648
649             IF ( k+1 > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 )  THEN
650                num_gp = num_gp + 1
651                gp_outside_of_building(4) = 1
652                location(num_gp,1) = i * dx
653                location(num_gp,2) = (j+1) * dy
654                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
655                ei(num_gp)     = e(k+1,j+1,i)
656                dissi(num_gp)  = diss(k+1,j+1,i)
657                de_dxi(num_gp) = de_dx(k+1,j+1,i)
658                de_dyi(num_gp) = de_dy(k+1,j+1,i)
659                de_dzi(num_gp) = de_dz(k+1,j+1,i)
660             ENDIF
661
662             IF ( k > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 )  THEN
663                num_gp = num_gp + 1
664                gp_outside_of_building(5) = 1
665                location(num_gp,1) = (i+1) * dx
666                location(num_gp,2) = j * dy
667                location(num_gp,3) = k * dz - 0.5_wp * dz
668                ei(num_gp)     = e(k,j,i+1)
669                dissi(num_gp)  = diss(k,j,i+1)
670                de_dxi(num_gp) = de_dx(k,j,i+1)
671                de_dyi(num_gp) = de_dy(k,j,i+1)
672                de_dzi(num_gp) = de_dz(k,j,i+1)
673             ENDIF
674
675             IF ( k > nzb_s_inner(j+1,i+1)  .OR.  nzb_s_inner(j+1,i+1) == 0 )  THEN
676                num_gp = num_gp + 1
677                gp_outside_of_building(6) = 1
678                location(num_gp,1) = (i+1) * dx
679                location(num_gp,2) = (j+1) * dy
680                location(num_gp,3) = k * dz - 0.5_wp * dz
681                ei(num_gp)     = e(k,j+1,i+1)
682                dissi(num_gp)  = diss(k,j+1,i+1)
683                de_dxi(num_gp) = de_dx(k,j+1,i+1)
684                de_dyi(num_gp) = de_dy(k,j+1,i+1)
685                de_dzi(num_gp) = de_dz(k,j+1,i+1)
686             ENDIF
687
688             IF ( k+1 > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 )  THEN
689                num_gp = num_gp + 1
690                gp_outside_of_building(7) = 1
691                location(num_gp,1) = (i+1) * dx
692                location(num_gp,2) = j * dy
693                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
694                ei(num_gp)     = e(k+1,j,i+1)
695                dissi(num_gp)  = diss(k+1,j,i+1)
696                de_dxi(num_gp) = de_dx(k+1,j,i+1)
697                de_dyi(num_gp) = de_dy(k+1,j,i+1)
698                de_dzi(num_gp) = de_dz(k+1,j,i+1)
699             ENDIF
700
701             IF ( k+1 > nzb_s_inner(j+1,i+1)  .OR.  nzb_s_inner(j+1,i+1) == 0)  THEN
702                num_gp = num_gp + 1
703                gp_outside_of_building(8) = 1
704                location(num_gp,1) = (i+1) * dx
705                location(num_gp,2) = (j+1) * dy
706                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
707                ei(num_gp)     = e(k+1,j+1,i+1)
708                dissi(num_gp)  = diss(k+1,j+1,i+1)
709                de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
710                de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
711                de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
712             ENDIF
713!
714!--          If all gridpoints are situated outside of a building, then the
715!--          ordinary interpolation scheme can be used.
716             IF ( num_gp == 8 )  THEN
717
718                = particles(n)%x - i * dx
719                y  = particles(n)%y - j * dy
720                aa = x**2          + y**2
721                bb = ( dx - x )**2 + y**2
722                cc = x**2          + ( dy - y )**2
723                dd = ( dx - x )**2 + ( dy - y )**2
724                gg = aa + bb + cc + dd
725 
726                e_int_l = ( ( gg - aa ) * e(k,j,i)   + ( gg - bb ) * e(k,j,i+1)   &
727                          + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) &
728                          ) / ( 3.0_wp * gg )
729   
730                IF ( k == nzt )  THEN
731                   e_int(n) = e_int_l
732                ELSE
733                   e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
734                               ( gg - bb ) * e(k+1,j,i+1) + &
735                               ( gg - cc ) * e(k+1,j+1,i) + &
736                               ( gg - dd ) * e(k+1,j+1,i+1) &
737                             ) / ( 3.0_wp * gg )
738                   e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
739                                       ( e_int_u - e_int_l )
740                ENDIF
741!
742!--             Needed to avoid NaN particle velocities (this might not be
743!--             required any more)
744                IF ( e_int(n) <= 0.0_wp )  THEN
745                   e_int(n) = 1.0E-20_wp
746                ENDIF
747!
748!--             Interpolate the TKE gradient along x (adopt incides i,j,k
749!--             and all position variables from above (TKE))
750                de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
751                                ( gg - bb ) * de_dx(k,j,i+1) + &
752                                ( gg - cc ) * de_dx(k,j+1,i) + &
753                                ( gg - dd ) * de_dx(k,j+1,i+1) &
754                              ) / ( 3.0_wp * gg )
755
756                IF ( ( k == nzt )  .OR.  ( k == nzb ) )  THEN
757                   de_dx_int(n) = de_dx_int_l
758                ELSE
759                   de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
760                                   ( gg - bb ) * de_dx(k+1,j,i+1) + &
761                                   ( gg - cc ) * de_dx(k+1,j+1,i) + &
762                                   ( gg - dd ) * de_dx(k+1,j+1,i+1) &
763                                 ) / ( 3.0_wp * gg )
764                   de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / &
765                                           dz * ( de_dx_int_u - de_dx_int_l )
766                ENDIF
767
768!
769!--             Interpolate the TKE gradient along y
770                de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
771                                ( gg - bb ) * de_dy(k,j,i+1) + &
772                                ( gg - cc ) * de_dy(k,j+1,i) + &
773                                ( gg - dd ) * de_dy(k,j+1,i+1) &
774                              ) / ( 3.0_wp * gg )
775                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
776                   de_dy_int(n) = de_dy_int_l
777                ELSE
778                   de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
779                                   ( gg - bb ) * de_dy(k+1,j,i+1) + &
780                                   ( gg - cc ) * de_dy(k+1,j+1,i) + &
781                                   ( gg - dd ) * de_dy(k+1,j+1,i+1) &
782                                 ) / ( 3.0_wp * gg )
783                   de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / &
784                                           dz * ( de_dy_int_u - de_dy_int_l )
785                ENDIF
786
787!
788!--             Interpolate the TKE gradient along z
789                IF ( zv(n) < 0.5_wp * dz )  THEN
790                   de_dz_int(n) = 0.0_wp
791                ELSE
792                   de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
793                                   ( gg - bb ) * de_dz(k,j,i+1) + &
794                                   ( gg - cc ) * de_dz(k,j+1,i) + &
795                                   ( gg - dd ) * de_dz(k,j+1,i+1) &
796                                 ) / ( 3.0_wp * gg )
797 
798                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
799                      de_dz_int(n) = de_dz_int_l
800                   ELSE
801                      de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
802                                      ( gg - bb ) * de_dz(k+1,j,i+1) + &
803                                      ( gg - cc ) * de_dz(k+1,j+1,i) + &
804                                      ( gg - dd ) * de_dz(k+1,j+1,i+1) &
805                                    ) / ( 3.0_wp * gg )
806                      de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) /&
807                                           dz * ( de_dz_int_u - de_dz_int_l )
808                   ENDIF
809                ENDIF
810
811!
812!--             Interpolate the dissipation of TKE
813                diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
814                               ( gg - bb ) * diss(k,j,i+1) + &
815                               ( gg - cc ) * diss(k,j+1,i) + &
816                               ( gg - dd ) * diss(k,j+1,i+1) &
817                             ) / ( 3.0_wp * gg )
818 
819                IF ( k == nzt )  THEN
820                   diss_int(n) = diss_int_l
821                ELSE
822                   diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
823                                  ( gg - bb ) * diss(k+1,j,i+1) + &
824                                  ( gg - cc ) * diss(k+1,j+1,i) + &
825                                  ( gg - dd ) * diss(k+1,j+1,i+1) &
826                                ) / ( 3.0_wp * gg )
827                   diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *&
828                                           ( diss_int_u - diss_int_l )
829                ENDIF
830!
831!--             Set flag for stochastic equation.
832                term_1_2(n) = 1.0_wp
833 
834             ELSE
835 
836!
837!--             If wall between gridpoint 1 and gridpoint 5, then
838!--             Neumann boundary condition has to be applied
839                IF ( gp_outside_of_building(1) == 1  .AND. &
840                     gp_outside_of_building(5) == 0 )  THEN
841                   num_gp = num_gp + 1
842                   location(num_gp,1) = i * dx + 0.5_wp * dx
843                   location(num_gp,2) = j * dy
844                   location(num_gp,3) = k * dz - 0.5_wp * dz
845                   ei(num_gp)     = e(k,j,i)
846                   dissi(num_gp)  = diss(k,j,i)
847                   de_dxi(num_gp) = 0.0_wp
848                   de_dyi(num_gp) = de_dy(k,j,i)
849                   de_dzi(num_gp) = de_dz(k,j,i)
850                ENDIF
851
852                IF ( gp_outside_of_building(5) == 1  .AND. &
853                     gp_outside_of_building(1) == 0 )  THEN
854                   num_gp = num_gp + 1
855                   location(num_gp,1) = i * dx + 0.5_wp * dx
856                   location(num_gp,2) = j * dy
857                   location(num_gp,3) = k * dz - 0.5_wp * dz
858                   ei(num_gp)     = e(k,j,i+1)
859                   dissi(num_gp)  = diss(k,j,i+1)
860                   de_dxi(num_gp) = 0.0_wp
861                   de_dyi(num_gp) = de_dy(k,j,i+1)
862                   de_dzi(num_gp) = de_dz(k,j,i+1)
863                ENDIF
864
865!
866!--             If wall between gridpoint 5 and gridpoint 6, then
867!--             then Neumann boundary condition has to be applied
868                IF ( gp_outside_of_building(5) == 1  .AND. &
869                     gp_outside_of_building(6) == 0 )  THEN
870                   num_gp = num_gp + 1
871                   location(num_gp,1) = (i+1) * dx
872                   location(num_gp,2) = j * dy + 0.5_wp * dy
873                   location(num_gp,3) = k * dz - 0.5_wp * dz
874                   ei(num_gp)     = e(k,j,i+1)
875                   dissi(num_gp)  = diss(k,j,i+1)
876                   de_dxi(num_gp) = de_dx(k,j,i+1)
877                   de_dyi(num_gp) = 0.0_wp
878                   de_dzi(num_gp) = de_dz(k,j,i+1)
879                ENDIF
880
881                IF ( gp_outside_of_building(6) == 1  .AND. &
882                     gp_outside_of_building(5) == 0 )  THEN
883                   num_gp = num_gp + 1
884                   location(num_gp,1) = (i+1) * dx
885                   location(num_gp,2) = j * dy + 0.5_wp * dy
886                   location(num_gp,3) = k * dz - 0.5_wp * dz
887                   ei(num_gp)     = e(k,j+1,i+1)
888                   dissi(num_gp)  = diss(k,j+1,i+1)
889                   de_dxi(num_gp) = de_dx(k,j+1,i+1)
890                   de_dyi(num_gp) = 0.0_wp
891                   de_dzi(num_gp) = de_dz(k,j+1,i+1)
892                ENDIF
893
894!
895!--             If wall between gridpoint 2 and gridpoint 6, then
896!--             Neumann boundary condition has to be applied
897                IF ( gp_outside_of_building(2) == 1  .AND. &
898                     gp_outside_of_building(6) == 0 )  THEN
899                   num_gp = num_gp + 1
900                   location(num_gp,1) = i * dx + 0.5_wp * dx
901                   location(num_gp,2) = (j+1) * dy
902                   location(num_gp,3) = k * dz - 0.5_wp * dz
903                   ei(num_gp)     = e(k,j+1,i)
904                   dissi(num_gp)  = diss(k,j+1,i)
905                   de_dxi(num_gp) = 0.0_wp
906                   de_dyi(num_gp) = de_dy(k,j+1,i)
907                   de_dzi(num_gp) = de_dz(k,j+1,i)
908                ENDIF
909
910                IF ( gp_outside_of_building(6) == 1  .AND. &
911                     gp_outside_of_building(2) == 0 )  THEN
912                   num_gp = num_gp + 1
913                   location(num_gp,1) = i * dx + 0.5_wp * dx
914                   location(num_gp,2) = (j+1) * dy
915                   location(num_gp,3) = k * dz - 0.5_wp * dz
916                   ei(num_gp)     = e(k,j+1,i+1)
917                   dissi(num_gp)  = diss(k,j+1,i+1)
918                   de_dxi(num_gp) = 0.0_wp
919                   de_dyi(num_gp) = de_dy(k,j+1,i+1)
920                   de_dzi(num_gp) = de_dz(k,j+1,i+1)
921                ENDIF
922
923!
924!--             If wall between gridpoint 1 and gridpoint 2, then
925!--             Neumann boundary condition has to be applied
926                IF ( gp_outside_of_building(1) == 1  .AND. &
927                     gp_outside_of_building(2) == 0 )  THEN
928                   num_gp = num_gp + 1
929                   location(num_gp,1) = i * dx
930                   location(num_gp,2) = j * dy + 0.5_wp * dy
931                   location(num_gp,3) = k * dz - 0.5_wp * dz
932                   ei(num_gp)     = e(k,j,i)
933                   dissi(num_gp)  = diss(k,j,i)
934                   de_dxi(num_gp) = de_dx(k,j,i)
935                   de_dyi(num_gp) = 0.0_wp
936                   de_dzi(num_gp) = de_dz(k,j,i) 
937                ENDIF
938
939                IF ( gp_outside_of_building(2) == 1  .AND. &
940                     gp_outside_of_building(1) == 0 )  THEN
941                   num_gp = num_gp + 1
942                   location(num_gp,1) = i * dx
943                   location(num_gp,2) = j * dy + 0.5_wp * dy
944                   location(num_gp,3) = k * dz - 0.5_wp * dz
945                   ei(num_gp)     = e(k,j+1,i)
946                   dissi(num_gp)  = diss(k,j+1,i)
947                   de_dxi(num_gp) = de_dx(k,j+1,i)
948                   de_dyi(num_gp) = 0.0_wp
949                   de_dzi(num_gp) = de_dz(k,j+1,i)
950                ENDIF
951
952!
953!--             If wall between gridpoint 3 and gridpoint 7, then
954!--             Neumann boundary condition has to be applied
955                IF ( gp_outside_of_building(3) == 1  .AND. &
956                     gp_outside_of_building(7) == 0 )  THEN
957                   num_gp = num_gp + 1
958                   location(num_gp,1) = i * dx + 0.5_wp * dx
959                   location(num_gp,2) = j * dy
960                   location(num_gp,3) = k * dz + 0.5_wp * dz
961                   ei(num_gp)     = e(k+1,j,i)
962                   dissi(num_gp)  = diss(k+1,j,i)
963                   de_dxi(num_gp) = 0.0_wp
964                   de_dyi(num_gp) = de_dy(k+1,j,i)
965                   de_dzi(num_gp) = de_dz(k+1,j,i) 
966                ENDIF
967
968                IF ( gp_outside_of_building(7) == 1  .AND. &
969                     gp_outside_of_building(3) == 0 )  THEN
970                   num_gp = num_gp + 1
971                   location(num_gp,1) = i * dx + 0.5_wp * dx
972                   location(num_gp,2) = j * dy
973                   location(num_gp,3) = k * dz + 0.5_wp * dz
974                   ei(num_gp)     = e(k+1,j,i+1)
975                   dissi(num_gp)  = diss(k+1,j,i+1)
976                   de_dxi(num_gp) = 0.0_wp
977                   de_dyi(num_gp) = de_dy(k+1,j,i+1)
978                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
979                ENDIF
980
981!
982!--             If wall between gridpoint 7 and gridpoint 8, then
983!--             Neumann boundary condition has to be applied
984                IF ( gp_outside_of_building(7) == 1  .AND. &
985                     gp_outside_of_building(8) == 0 )  THEN
986                   num_gp = num_gp + 1
987                   location(num_gp,1) = (i+1) * dx
988                   location(num_gp,2) = j * dy + 0.5_wp * dy
989                   location(num_gp,3) = k * dz + 0.5_wp * dz
990                   ei(num_gp)     = e(k+1,j,i+1)
991                   dissi(num_gp)  = diss(k+1,j,i+1)
992                   de_dxi(num_gp) = de_dx(k+1,j,i+1)
993                   de_dyi(num_gp) = 0.0_wp
994                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
995                ENDIF
996
997                IF ( gp_outside_of_building(8) == 1  .AND. &
998                     gp_outside_of_building(7) == 0 )  THEN
999                   num_gp = num_gp + 1
1000                   location(num_gp,1) = (i+1) * dx
1001                   location(num_gp,2) = j * dy + 0.5_wp * dy
1002                   location(num_gp,3) = k * dz + 0.5_wp * 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) = 0.0_wp
1007                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
1008                ENDIF
1009
1010!
1011!--             If wall between gridpoint 4 and gridpoint 8, then
1012!--             Neumann boundary condition has to be applied
1013                IF ( gp_outside_of_building(4) == 1  .AND. &
1014                     gp_outside_of_building(8) == 0 )  THEN
1015                   num_gp = num_gp + 1
1016                   location(num_gp,1) = i * dx + 0.5_wp * dx
1017                   location(num_gp,2) = (j+1) * dy
1018                   location(num_gp,3) = k * dz + 0.5_wp * dz
1019                   ei(num_gp)     = e(k+1,j+1,i)
1020                   dissi(num_gp)  = diss(k+1,j+1,i)
1021                   de_dxi(num_gp) = 0.0_wp
1022                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
1023                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
1024                ENDIF
1025
1026                IF ( gp_outside_of_building(8) == 1  .AND. &
1027                     gp_outside_of_building(4) == 0 )  THEN
1028                   num_gp = num_gp + 1
1029                   location(num_gp,1) = i * dx + 0.5_wp * dx
1030                   location(num_gp,2) = (j+1) * dy
1031                   location(num_gp,3) = k * dz + 0.5_wp * dz
1032                   ei(num_gp)     = e(k+1,j+1,i+1)
1033                   dissi(num_gp)  = diss(k+1,j+1,i+1)
1034                   de_dxi(num_gp) = 0.0_wp
1035                   de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
1036                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
1037                ENDIF
1038
1039!
1040!--             If wall between gridpoint 3 and gridpoint 4, then
1041!--             Neumann boundary condition has to be applied
1042                IF ( gp_outside_of_building(3) == 1  .AND. &
1043                     gp_outside_of_building(4) == 0 )  THEN
1044                   num_gp = num_gp + 1
1045                   location(num_gp,1) = i * dx
1046                   location(num_gp,2) = j * dy + 0.5_wp * dy
1047                   location(num_gp,3) = k * dz + 0.5_wp * dz
1048                   ei(num_gp)     = e(k+1,j,i)
1049                   dissi(num_gp)  = diss(k+1,j,i)
1050                   de_dxi(num_gp) = de_dx(k+1,j,i)
1051                   de_dyi(num_gp) = 0.0_wp
1052                   de_dzi(num_gp) = de_dz(k+1,j,i)
1053                ENDIF
1054
1055                IF ( gp_outside_of_building(4) == 1  .AND. &
1056                     gp_outside_of_building(3) == 0 )  THEN
1057                   num_gp = num_gp + 1
1058                   location(num_gp,1) = i * dx
1059                   location(num_gp,2) = j * dy + 0.5_wp * dy
1060                   location(num_gp,3) = k * dz + 0.5_wp * dz
1061                   ei(num_gp)     = e(k+1,j+1,i)
1062                   dissi(num_gp)  = diss(k+1,j+1,i)
1063                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
1064                   de_dyi(num_gp) = 0.0_wp
1065                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
1066                ENDIF
1067
1068!
1069!--             If wall between gridpoint 1 and gridpoint 3, then
1070!--             Neumann boundary condition has to be applied
1071!--             (only one case as only building beneath is possible)
1072                IF ( gp_outside_of_building(1) == 0  .AND. &
1073                     gp_outside_of_building(3) == 1 )  THEN
1074                   num_gp = num_gp + 1
1075                   location(num_gp,1) = i * dx
1076                   location(num_gp,2) = j * dy
1077                   location(num_gp,3) = k * dz
1078                   ei(num_gp)     = e(k+1,j,i)
1079                   dissi(num_gp)  = diss(k+1,j,i)
1080                   de_dxi(num_gp) = de_dx(k+1,j,i)
1081                   de_dyi(num_gp) = de_dy(k+1,j,i)
1082                   de_dzi(num_gp) = 0.0_wp
1083                ENDIF
1084
1085!
1086!--             If wall between gridpoint 5 and gridpoint 7, then
1087!--             Neumann boundary condition has to be applied
1088!--             (only one case as only building beneath is possible)
1089                IF ( gp_outside_of_building(5) == 0  .AND. &
1090                     gp_outside_of_building(7) == 1 )  THEN
1091                   num_gp = num_gp + 1
1092                   location(num_gp,1) = (i+1) * dx
1093                   location(num_gp,2) = j * dy
1094                   location(num_gp,3) = k * dz
1095                   ei(num_gp)     = e(k+1,j,i+1)
1096                   dissi(num_gp)  = diss(k+1,j,i+1)
1097                   de_dxi(num_gp) = de_dx(k+1,j,i+1)
1098                   de_dyi(num_gp) = de_dy(k+1,j,i+1)
1099                   de_dzi(num_gp) = 0.0_wp
1100                ENDIF
1101
1102!
1103!--             If wall between gridpoint 2 and gridpoint 4, then
1104!--             Neumann boundary condition has to be applied
1105!--             (only one case as only building beneath is possible)
1106                IF ( gp_outside_of_building(2) == 0  .AND. &
1107                     gp_outside_of_building(4) == 1 )  THEN
1108                   num_gp = num_gp + 1
1109                   location(num_gp,1) = i * dx
1110                   location(num_gp,2) = (j+1) * dy
1111                   location(num_gp,3) = k * dz
1112                   ei(num_gp)     = e(k+1,j+1,i)
1113                   dissi(num_gp)  = diss(k+1,j+1,i)
1114                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
1115                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
1116                   de_dzi(num_gp) = 0.0_wp
1117                ENDIF
1118
1119!
1120!--             If wall between gridpoint 6 and gridpoint 8, then
1121!--             Neumann boundary condition has to be applied
1122!--             (only one case as only building beneath is possible)
1123                IF ( gp_outside_of_building(6) == 0  .AND. &
1124                     gp_outside_of_building(8) == 1 )  THEN
1125                   num_gp = num_gp + 1
1126                   location(num_gp,1) = (i+1) * dx
1127                   location(num_gp,2) = (j+1) * dy
1128                   location(num_gp,3) = k * dz
1129                   ei(num_gp)     = e(k+1,j+1,i+1)
1130                   dissi(num_gp)  = diss(k+1,j+1,i+1)
1131                   de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
1132                   de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
1133                   de_dzi(num_gp) = 0.0_wp
1134                ENDIF
1135 
1136!
1137!--             Carry out the interpolation
1138                IF ( num_gp == 1 )  THEN
1139!
1140!--                If only one of the gridpoints is situated outside of the
1141!--                building, it follows that the values at the particle
1142!--                location are the same as the gridpoint values
1143                   e_int(n)     = ei(num_gp)   
1144                   diss_int(n)  = dissi(num_gp) 
1145                   de_dx_int(n) = de_dxi(num_gp)
1146                   de_dy_int(n) = de_dyi(num_gp)
1147                   de_dz_int(n) = de_dzi(num_gp)
1148!
1149!--                Set flag for stochastic equation which disables calculation
1150!--                of drift and memory term near topography.
1151                   term_1_2(n) = 0.0_wp
1152                ELSE IF ( num_gp > 1 )  THEN
1153 
1154                   d_sum = 0.0_wp
1155!
1156!--                Evaluation of the distances between the gridpoints
1157!--                contributing to the interpolated values, and the particle
1158!--                location
1159                   DO  agp = 1, num_gp
1160                      d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2  &
1161                                   + ( particles(n)%y-location(agp,2) )**2  &
1162                                   + ( zv(n)-location(agp,3) )**2
1163                      d_sum        = d_sum + d_gp_pl(agp)
1164                   ENDDO
1165 
1166!
1167!--                Finally the interpolation can be carried out
1168                   e_int(n)     = 0.0_wp 
1169                   diss_int(n)  = 0.0_wp 
1170                   de_dx_int(n) = 0.0_wp 
1171                   de_dy_int(n) = 0.0_wp 
1172                   de_dz_int(n) = 0.0_wp 
1173                   DO  agp = 1, num_gp
1174                      e_int(n)     = e_int(n)     + ( d_sum - d_gp_pl(agp) ) * &
1175                                             ei(agp) / ( (num_gp-1) * d_sum )
1176                      diss_int(n)  = diss_int(n)  + ( d_sum - d_gp_pl(agp) ) * &
1177                                          dissi(agp) / ( (num_gp-1) * d_sum )
1178                      de_dx_int(n) = de_dx_int(n) + ( d_sum - d_gp_pl(agp) ) * &
1179                                         de_dxi(agp) / ( (num_gp-1) * d_sum )
1180                      de_dy_int(n) = de_dy_int(n) + ( d_sum - d_gp_pl(agp) ) * &
1181                                         de_dyi(agp) / ( (num_gp-1) * d_sum )
1182                      de_dz_int(n) = de_dz_int(n) + ( d_sum - d_gp_pl(agp) ) * &
1183                                         de_dzi(agp) / ( (num_gp-1) * d_sum )
1184                   ENDDO
1185 
1186                ENDIF
1187                e_int(n)     = MAX( 1E-20_wp, e_int(n)     )
1188                diss_int(n)  = MAX( 1E-20_wp, diss_int(n)  )
1189                de_dx_int(n) = MAX( 1E-20_wp, de_dx_int(n) )
1190                de_dy_int(n) = MAX( 1E-20_wp, de_dy_int(n) ) 
1191                de_dz_int(n) = MAX( 1E-20_wp, de_dz_int(n) ) 
1192!
1193!--             Set flag for stochastic equation which disables calculation
1194!--             of drift and memory term near topography.
1195                term_1_2(n) = 0.0_wp
1196             ENDIF
1197          ENDDO
1198       ENDIF
1199
1200       DO nb = 0,7
1201          i = ip + block_offset(nb)%i_off
1202          j = jp + block_offset(nb)%j_off
1203          k = kp + block_offset(nb)%k_off
1204
1205          DO  n = start_index(nb), end_index(nb)
1206!
1207!--          Vertical interpolation of the horizontally averaged SGS TKE and
1208!--          resolved-scale velocity variances and use the interpolated values
1209!--          to calculate the coefficient fs, which is a measure of the ratio
1210!--          of the subgrid-scale turbulent kinetic energy to the total amount
1211!--          of turbulent kinetic energy.
1212             IF ( k == 0 )  THEN
1213                e_mean_int = hom(0,1,8,0)
1214             ELSE
1215                e_mean_int = hom(k,1,8,0) +                                    &
1216                                           ( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
1217                                           ( zu(k+1) - zu(k) ) *               &
1218                                           ( zv(n) - zu(k) )
1219             ENDIF
1220
1221             kw = kp - 1
1222
1223             IF ( k == 0 )  THEN
1224                aa  = hom(k+1,1,30,0)  * ( zv(n) / &
1225                                         ( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
1226                bb  = hom(k+1,1,31,0)  * ( zv(n) / &
1227                                         ( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
1228                cc  = hom(kw+1,1,32,0) * ( zv(n) / &
1229                                         ( 1.0_wp * ( zw(kw+1) - zw(kw) ) ) )
1230             ELSE
1231                aa  = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) *    &
1232                           ( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
1233                bb  = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) *    &
1234                           ( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
1235                cc  = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) *   &
1236                           ( ( zv(n) - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
1237             ENDIF
1238
1239             vv_int = ( 1.0_wp / 3.0_wp ) * ( aa + bb + cc )
1240!
1241!--          Needed to avoid NaN particle velocities. The value of 1.0 is just
1242!--          an educated guess for the given case.
1243             IF ( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int == 0.0_wp )  THEN
1244                fs_int(n) = 1.0_wp
1245             ELSE
1246                fs_int(n) = ( 2.0_wp / 3.0_wp ) * e_mean_int /                 &
1247                            ( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int )
1248             ENDIF
1249
1250          ENDDO
1251       ENDDO
1252
1253       DO  n = 1, number_of_particles
1254
1255         rg(n,1) = random_gauss( iran_part, 5.0_wp )
1256         rg(n,2) = random_gauss( iran_part, 5.0_wp )
1257         rg(n,3) = random_gauss( iran_part, 5.0_wp )
1258
1259       ENDDO
1260
1261       DO  n = 1, number_of_particles
1262!
1263!--       Calculate the Lagrangian timescale according to Weil et al. (2004).
1264          lagr_timescale = ( 4.0_wp * e_int(n) + 1E-20_wp ) / &
1265                           ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp )
1266
1267!
1268!--       Calculate the next particle timestep. dt_gap is the time needed to
1269!--       complete the current LES timestep.
1270          dt_gap = dt_3d - particles(n)%dt_sum
1271          dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale, dt_gap )
1272
1273!
1274!--       The particle timestep should not be too small in order to prevent
1275!--       the number of particle timesteps of getting too large
1276          IF ( dt_particle(n) < dt_min_part  .AND.  dt_min_part < dt_gap )  THEN
1277             dt_particle(n) = dt_min_part
1278          ENDIF
1279
1280!
1281!--       Calculate the SGS velocity components
1282          IF ( particles(n)%age == 0.0_wp )  THEN
1283!
1284!--          For new particles the SGS components are derived from the SGS
1285!--          TKE. Limit the Gaussian random number to the interval
1286!--          [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities
1287!--          from becoming unrealistically large.
1288             particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) *   &
1289                                  ( rg(n,1) - 1.0_wp )
1290             particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) *   &
1291                                  ( rg(n,2) - 1.0_wp )
1292             particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) *   &
1293                                  ( rg(n,3) - 1.0_wp )
1294
1295          ELSE
1296!
1297!--          Restriction of the size of the new timestep: compared to the
1298!--          previous timestep the increase must not exceed 200%
1299
1300             dt_particle_m = particles(n)%age - particles(n)%age_m
1301             IF ( dt_particle(n) > 2.0_wp * dt_particle_m )  THEN
1302                dt_particle(n) = 2.0_wp * dt_particle_m
1303             ENDIF
1304
1305!
1306!--          For old particles the SGS components are correlated with the
1307!--          values from the previous timestep. Random numbers have also to
1308!--          be limited (see above).
1309!--          As negative values for the subgrid TKE are not allowed, the
1310!--          change of the subgrid TKE with time cannot be smaller than
1311!--          -e_int(n)/dt_particle. This value is used as a lower boundary
1312!--          value for the change of TKE
1313
1314             de_dt_min = - e_int(n) / dt_particle(n)
1315
1316             de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m
1317
1318             IF ( de_dt < de_dt_min )  THEN
1319                de_dt = de_dt_min
1320             ENDIF
1321
1322             CALL weil_stochastic_eq(particles(n)%rvar1, fs_int(n), e_int(n),  & 
1323                                     de_dx_int(n), de_dt, diss_int(n),         &
1324                                     dt_particle(n), rg(n,1), term_1_2(n) )
1325
1326             CALL weil_stochastic_eq(particles(n)%rvar2, fs_int(n), e_int(n),  & 
1327                                     de_dy_int(n), de_dt, diss_int(n),         &
1328                                     dt_particle(n), rg(n,2), term_1_2(n) )
1329
1330             CALL weil_stochastic_eq(particles(n)%rvar3, fs_int(n), e_int(n),  & 
1331                                     de_dz_int(n), de_dt, diss_int(n),         &
1332                                     dt_particle(n), rg(n,3), term_1_2(n) )
1333
1334          ENDIF
1335
1336          u_int(n) = u_int(n) + particles(n)%rvar1
1337          v_int(n) = v_int(n) + particles(n)%rvar2
1338          w_int(n) = w_int(n) + particles(n)%rvar3
1339!
1340!--       Store the SGS TKE of the current timelevel which is needed for
1341!--       for calculating the SGS particle velocities at the next timestep
1342          particles(n)%e_m = e_int(n)
1343       ENDDO
1344
1345    ELSE
1346!
1347!--    If no SGS velocities are used, only the particle timestep has to
1348!--    be set
1349       dt_particle = dt_3d
1350
1351    ENDIF
1352!
1353!-- Store the old age of the particle ( needed to prevent that a
1354!-- particle crosses several PEs during one timestep, and for the
1355!-- evaluation of the subgrid particle velocity fluctuations )
1356    particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age
1357
1358    dens_ratio = particle_groups(particles(1:number_of_particles)%group)%density_ratio
1359
1360    IF ( ANY( dens_ratio == 0.0_wp ) )  THEN
1361       DO  n = 1, number_of_particles
1362
1363!
1364!--       Particle advection
1365          IF ( dens_ratio(n) == 0.0_wp )  THEN
1366!
1367!--          Pure passive transport (without particle inertia)
1368             particles(n)%x = xv(n) + u_int(n) * dt_particle(n)
1369             particles(n)%y = yv(n) + v_int(n) * dt_particle(n)
1370             particles(n)%z = zv(n) + w_int(n) * dt_particle(n)
1371
1372             particles(n)%speed_x = u_int(n)
1373             particles(n)%speed_y = v_int(n)
1374             particles(n)%speed_z = w_int(n)
1375
1376          ELSE
1377!
1378!--          Transport of particles with inertia
1379             particles(n)%x = particles(n)%x + particles(n)%speed_x * &
1380                                               dt_particle(n)
1381             particles(n)%y = particles(n)%y + particles(n)%speed_y * &
1382                                               dt_particle(n)
1383             particles(n)%z = particles(n)%z + particles(n)%speed_z * &
1384                                               dt_particle(n)
1385
1386!
1387!--          Update of the particle velocity
1388             IF ( cloud_droplets )  THEN
1389!
1390!--             Terminal velocity is computed for vertical direction (Rogers et
1391!--             al., 1993, J. Appl. Meteorol.)
1392                diameter = particles(n)%radius * 2000.0_wp !diameter in mm
1393                IF ( diameter <= d0_rog )  THEN
1394                   w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
1395                ELSE
1396                   w_s = a_rog - b_rog * EXP( -c_rog * diameter )
1397                ENDIF
1398
1399!
1400!--             If selected, add random velocities following Soelch and Kaercher
1401!--             (2010, Q. J. R. Meteorol. Soc.)
1402                IF ( use_sgs_for_particles )  THEN
1403                   lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
1404                   RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
1405                   sigma          = SQRT( e(kp,jp,ip) )
1406
1407                   rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
1408                   rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
1409                   rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
1410
1411                   particles(n)%rvar1 = RL * particles(n)%rvar1 +              &
1412                                        SQRT( 1.0_wp - RL**2 ) * sigma * rg1
1413                   particles(n)%rvar2 = RL * particles(n)%rvar2 +              &
1414                                        SQRT( 1.0_wp - RL**2 ) * sigma * rg2
1415                   particles(n)%rvar3 = RL * particles(n)%rvar3 +              &
1416                                        SQRT( 1.0_wp - RL**2 ) * sigma * rg3
1417
1418                   particles(n)%speed_x = u_int(n) + particles(n)%rvar1
1419                   particles(n)%speed_y = v_int(n) + particles(n)%rvar2
1420                   particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
1421                ELSE
1422                   particles(n)%speed_x = u_int(n)
1423                   particles(n)%speed_y = v_int(n)
1424                   particles(n)%speed_z = w_int(n) - w_s
1425                ENDIF
1426
1427             ELSE
1428
1429                IF ( use_sgs_for_particles )  THEN
1430                   exp_arg  = particle_groups(particles(n)%group)%exp_arg
1431                   exp_term = EXP( -exp_arg * dt_particle(n) )
1432                ELSE
1433                   exp_arg  = particle_groups(particles(n)%group)%exp_arg
1434                   exp_term = particle_groups(particles(n)%group)%exp_term
1435                ENDIF
1436                particles(n)%speed_x = particles(n)%speed_x * exp_term +         &
1437                                       u_int(n) * ( 1.0_wp - exp_term )
1438                particles(n)%speed_y = particles(n)%speed_y * exp_term +         &
1439                                       v_int(n) * ( 1.0_wp - exp_term )
1440                particles(n)%speed_z = particles(n)%speed_z * exp_term +         &
1441                                       ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * &
1442                                       g / exp_arg ) * ( 1.0_wp - exp_term )
1443             ENDIF
1444
1445          ENDIF
1446
1447       ENDDO
1448   
1449    ELSE
1450
1451       DO  n = 1, number_of_particles
1452
1453!--       Transport of particles with inertia
1454          particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)
1455          particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)
1456          particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)
1457!
1458!--       Update of the particle velocity
1459          IF ( cloud_droplets )  THEN
1460!
1461!--          Terminal velocity is computed for vertical direction (Rogers et al.,
1462!--          1993, J. Appl. Meteorol.)
1463             diameter = particles(n)%radius * 2000.0_wp !diameter in mm
1464             IF ( diameter <= d0_rog )  THEN
1465                w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
1466             ELSE
1467                w_s = a_rog - b_rog * EXP( -c_rog * diameter )
1468             ENDIF
1469
1470!
1471!--          If selected, add random velocities following Soelch and Kaercher
1472!--          (2010, Q. J. R. Meteorol. Soc.)
1473             IF ( use_sgs_for_particles )  THEN
1474                 lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
1475                 RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
1476                 sigma          = SQRT( e(kp,jp,ip) )
1477
1478                 rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
1479                 rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
1480                 rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
1481
1482                 particles(n)%rvar1 = RL * particles(n)%rvar1 +                &
1483                                      SQRT( 1.0_wp - RL**2 ) * sigma * rg1
1484                 particles(n)%rvar2 = RL * particles(n)%rvar2 +                &
1485                                      SQRT( 1.0_wp - RL**2 ) * sigma * rg2
1486                 particles(n)%rvar3 = RL * particles(n)%rvar3 +                &
1487                                      SQRT( 1.0_wp - RL**2 ) * sigma * rg3
1488
1489                 particles(n)%speed_x = u_int(n) + particles(n)%rvar1
1490                 particles(n)%speed_y = v_int(n) + particles(n)%rvar2
1491                 particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
1492             ELSE
1493                 particles(n)%speed_x = u_int(n)
1494                 particles(n)%speed_y = v_int(n)
1495                 particles(n)%speed_z = w_int(n) - w_s
1496             ENDIF
1497
1498          ELSE
1499
1500             IF ( use_sgs_for_particles )  THEN
1501                exp_arg  = particle_groups(particles(n)%group)%exp_arg
1502                exp_term = EXP( -exp_arg * dt_particle(n) )
1503             ELSE
1504                exp_arg  = particle_groups(particles(n)%group)%exp_arg
1505                exp_term = particle_groups(particles(n)%group)%exp_term
1506             ENDIF
1507             particles(n)%speed_x = particles(n)%speed_x * exp_term +             &
1508                                    u_int(n) * ( 1.0_wp - exp_term )
1509             particles(n)%speed_y = particles(n)%speed_y * exp_term +             &
1510                                    v_int(n) * ( 1.0_wp - exp_term )
1511             particles(n)%speed_z = particles(n)%speed_z * exp_term +             &
1512                                    ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &
1513                                    exp_arg ) * ( 1.0_wp - exp_term )
1514          ENDIF
1515
1516       ENDDO
1517
1518    ENDIF
1519
1520    DO  n = 1, number_of_particles
1521!
1522!--    Increment the particle age and the total time that the particle
1523!--    has advanced within the particle timestep procedure
1524       particles(n)%age    = particles(n)%age    + dt_particle(n)
1525       particles(n)%dt_sum = particles(n)%dt_sum + dt_particle(n)
1526
1527!
1528!--    Check whether there is still a particle that has not yet completed
1529!--    the total LES timestep
1530       IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp )  THEN
1531          dt_3d_reached_l = .FALSE.
1532       ENDIF
1533
1534    ENDDO
1535
1536    CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
1537
1538
1539 END SUBROUTINE lpm_advec
1540
1541! Description:
1542! ------------
1543!> Calculation of subgrid-scale particle speed using the stochastic model
1544!> of Weil et al. (2004, JAS, 61, 2877-2887).
1545!------------------------------------------------------------------------------!
1546 SUBROUTINE weil_stochastic_eq( v_sgs, fs_n, e_n, dedxi_n, dedt_n, diss_n,     &
1547                                dt_n, rg_n, fac )
1548
1549    USE kinds
1550
1551    USE particle_attributes,                                                   &
1552        ONLY:  c_0, sgs_wf_part
1553
1554    IMPLICIT NONE
1555
1556    REAL(wp) ::  a1      !< dummy argument
1557    REAL(wp) ::  dedt_n  !< time derivative of TKE at particle position
1558    REAL(wp) ::  dedxi_n !< horizontal derivative of TKE at particle position
1559    REAL(wp) ::  diss_n  !< dissipation at particle position
1560    REAL(wp) ::  dt_n    !< particle timestep
1561    REAL(wp) ::  e_n     !< TKE at particle position
1562    REAL(wp) ::  fac     !< flag to identify adjacent topography
1563    REAL(wp) ::  fs_n    !< weighting factor to prevent that subgrid-scale particle speed becomes too large
1564    REAL(wp) ::  sgs_w   !< constant (1/3)
1565    REAL(wp) ::  rg_n    !< random number
1566    REAL(wp) ::  term1   !< memory term
1567    REAL(wp) ::  term2   !< drift correction term
1568    REAL(wp) ::  term3   !< random term
1569    REAL(wp) ::  v_sgs   !< subgrid-scale velocity component
1570
1571!
1572!-- Please note, terms 1 and 2 (drift and memory term, respectively) are
1573!-- multiplied by a flag to switch of both terms near topography.
1574!-- This is necessary, as both terms may cause a subgrid-scale velocity build up
1575!-- if particles are trapped in regions with very small TKE, e.g. in narrow street
1576!-- canyons resolved by only a few grid points. Hence, term 1 and term 2 are
1577!-- disabled if one of the adjacent grid points belongs to topography.
1578!-- Moreover, in this case, the  previous subgrid-scale component is also set
1579!-- to zero.
1580
1581    a1 = fs_n * c_0 * diss_n
1582!
1583!-- Memory term
1584    term1 = - a1 * v_sgs * dt_n / ( 4.0_wp * sgs_wf_part * e_n + 1E-20_wp )    &
1585                 * fac
1586!
1587!-- Drift correction term
1588    term2 = ( ( dedt_n * v_sgs / e_n ) + dedxi_n ) * 0.5_wp * dt_n              &
1589                 * fac
1590!
1591!-- Random term
1592    term3 = SQRT( MAX( a1, 1E-20 ) ) * ( rg_n - 1.0_wp ) * SQRT( dt_n )
1593!
1594!-- In cese one of the adjacent grid-boxes belongs to topograhy, the previous
1595!-- subgrid-scale velocity component is set to zero, in order to prevent a
1596!-- velocity build-up.
1597
1598!-- This case, set also previous subgrid-scale component to zero.
1599    v_sgs = v_sgs * fac + term1 + term2 + term3
1600
1601 END SUBROUTINE weil_stochastic_eq
Note: See TracBrowser for help on using the repository browser.