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

Last change on this file since 1929 was 1929, checked in by suehring, 5 years ago

several bugfixes in particle model and serial mode

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