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

Last change on this file since 2606 was 2606, checked in by schwenkel, 4 years ago

Modified particle box location and further changes in particle model

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