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

Last change on this file since 2036 was 2001, checked in by knoop, 8 years ago

last commit documented

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