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

Last change on this file since 2610 was 2610, checked in by raasch, 6 years ago

use of wildcards in file connection statements enabled, bugfix in logarithmic interpolation of v-component for particles (usws was used by mistake)

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