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

Last change on this file since 2699 was 2698, checked in by suehring, 7 years ago

Particle reflections at downward-facing walls; revision of particle speed interpolations at walls; bugfixes in get_topography_index and in date constants

  • Property svn:keywords set to Id
File size: 51.8 KB
Line 
1!> @file lpm_advec.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
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! Particle interpolations at walls in case of SGS velocities revised and not
23! required parts are removed. (responsible Philipp Thiele)
24!
25! Bugfix in get_topography_top_index
26!
27! Former revisions:
28! -----------------
29! $Id: lpm_advec.f90 2698 2017-12-14 18:46:24Z kanani $
30! Removed indices ilog and jlog which are no longer needed since particle box
31! locations are identical to scalar boxes and topography.
32!
33! 2628 2017-11-20 12:40:38Z raasch
34! bugfix in logarithmic interpolation of v-component (usws was used by mistake)
35!
36! 2606 2017-11-10 10:36:31Z schwenkel
37! Changed particle box locations: center of particle box now coincides
38! with scalar grid point of same index.
39! Renamed module and subroutines: lpm_pack_arrays_mod -> lpm_pack_and_sort_mod
40! lpm_pack_all_arrays -> lpm_sort_in_subboxes, lpm_pack_arrays -> lpm_pack
41! lpm_sort -> lpm_sort_timeloop_done
42!
43! 2417 2017-09-06 15:22:27Z suehring
44! Particle loops adapted for sub-box structure, i.e. for each sub-box the
45! particle loop runs from start_index up to end_index instead from 1 to
46! number_of_particles. This way, it is possible to skip unnecessary
47! computations for particles that already completed the LES timestep.
48!
49! 2318 2017-07-20 17:27:44Z suehring
50! Get topography top index via Function call
51!
52! 2317 2017-07-20 17:27:19Z suehring
53!
54! 2232 2017-05-30 17:47:52Z suehring
55! Adjustments to new topography and surface concept
56!
57! 2100 2017-01-05 16:40:16Z suehring
58! Prevent extremely large SGS-velocities in regions where TKE is zero, e.g.
59! at the begin of simulations and/or in non-turbulent regions.
60!
61! 2000 2016-08-20 18:09:15Z knoop
62! Forced header and separation lines into 80 columns
63!
64! 1936 2016-06-13 13:37:44Z suehring
65! Formatting adjustments
66!
67! 1929 2016-06-09 16:25:25Z suehring
68! Put stochastic equation in an extra subroutine.
69! Set flag for stochastic equation to communicate whether a particle is near
70! topography. This case, memory and drift term are disabled in the Weil equation.
71!
72! Enable vertical logarithmic interpolation also above topography. This case,
73! set a lower limit for the friction velocity, as it can become very small
74! in narrow street canyons, leading to too large particle speeds.
75!
76! 1888 2016-04-21 12:20:49Z suehring
77! Bugfix concerning logarithmic interpolation of particle speed
78!
79! 1822 2016-04-07 07:49:42Z hoffmann
80! Random velocity fluctuations for particles added. Terminal fall velocity
81! for droplets is calculated from a parameterization (which is better than
82! the previous, physically correct calculation, which demands a very short
83! time step that is not used in the model).
84!
85! Unused variables deleted.
86!
87! 1691 2015-10-26 16:17:44Z maronga
88! Renamed prandtl_layer to constant_flux_layer.
89!
90! 1685 2015-10-08 07:32:13Z raasch
91! TKE check for negative values (so far, only zero value was checked)
92! offset_ocean_nzt_m1 removed
93!
94! 1682 2015-10-07 23:56:08Z knoop
95! Code annotations made doxygen readable
96!
97! 1583 2015-04-15 12:16:27Z suehring
98! Bugfix: particle advection within Prandtl-layer in case of Galilei
99! transformation.
100!
101! 1369 2014-04-24 05:57:38Z raasch
102! usage of module interfaces removed
103!
104! 1359 2014-04-11 17:15:14Z hoffmann
105! New particle structure integrated.
106! Kind definition added to all floating point numbers.
107!
108! 1322 2014-03-20 16:38:49Z raasch
109! REAL constants defined as wp_kind
110!
111! 1320 2014-03-20 08:40:49Z raasch
112! ONLY-attribute added to USE-statements,
113! kind-parameters added to all INTEGER and REAL declaration statements,
114! kinds are defined in new module kinds,
115! revision history before 2012 removed,
116! comment fields (!:) to be used for variable explanations added to
117! all variable declaration statements
118!
119! 1314 2014-03-14 18:25:17Z suehring
120! Vertical logarithmic interpolation of horizontal particle speed for particles
121! between roughness height and first vertical grid level.
122!
123! 1036 2012-10-22 13:43:42Z raasch
124! code put under GPL (PALM 3.9)
125!
126! 849 2012-03-15 10:35:09Z raasch
127! initial revision (former part of advec_particles)
128!
129!
130! Description:
131! ------------
132!> Calculation of new particle positions due to advection using a simple Euler
133!> scheme. Particles may feel inertia effects. SGS transport can be included
134!> using the stochastic model of Weil et al. (2004, JAS, 61, 2877-2887).
135!------------------------------------------------------------------------------!
136 SUBROUTINE lpm_advec (ip,jp,kp)
137 
138
139    USE arrays_3d,                                                             &
140        ONLY:  de_dx, de_dy, de_dz, diss, e, km, u, v, w, zu, zw
141
142    USE cpulog
143
144    USE pegrid
145
146    USE control_parameters,                                                    &
147        ONLY:  atmos_ocean_sign, cloud_droplets, constant_flux_layer, dt_3d,   &
148               dt_3d_reached_l, dz, g, kappa, topography, u_gtrans, v_gtrans
149
150    USE grid_variables,                                                        &
151        ONLY:  ddx, dx, ddy, dy
152       
153    USE indices,                                                               &
154        ONLY:  nzb, nzt, wall_flags_0
155       
156    USE kinds
157   
158    USE particle_attributes,                                                   &
159        ONLY:  block_offset, c_0, dt_min_part, grid_particles,                 &
160               iran_part, log_z_z0, number_of_particles, number_of_sublayers,  &
161               particles, particle_groups, offset_ocean_nzt, sgs_wf_part,      &
162               use_sgs_for_particles, vertical_particle_advection, z0_av_global
163       
164    USE statistics,                                                            &
165        ONLY:  hom
166
167    USE surface_mod,                                                           &
168        ONLY:  get_topography_top_index_ji, surf_def_h, surf_lsm_h, surf_usm_h
169
170    IMPLICIT NONE
171
172    LOGICAL ::  subbox_at_wall !< flag to see if the current subgridbox is adjacent to a wall
173
174    INTEGER(iwp) ::  agp                         !< loop variable
175    INTEGER(iwp) ::  gp_outside_of_building(1:8) !< number of grid points used for particle interpolation in case of topography
176    INTEGER(iwp) ::  i                           !< index variable along x
177    INTEGER(iwp) ::  ip                          !< index variable along x
178    INTEGER(iwp) ::  j                           !< index variable along y
179    INTEGER(iwp) ::  jp                          !< index variable along y
180    INTEGER(iwp) ::  k                           !< index variable along z
181    INTEGER(iwp) ::  k_wall                      !< vertical index of topography top
182    INTEGER(iwp) ::  kp                          !< index variable along z
183    INTEGER(iwp) ::  kw                          !< index variable along z
184    INTEGER(iwp) ::  n                           !< loop variable over all particles in a grid box
185    INTEGER(iwp) ::  nb                          !< block number particles are sorted in
186    INTEGER(iwp) ::  num_gp                      !< number of adjacent grid points inside topography
187    INTEGER(iwp) ::  surf_start                  !< Index on surface data-type for current grid box
188
189    INTEGER(iwp), DIMENSION(0:7) ::  start_index !< start particle index for current block
190    INTEGER(iwp), DIMENSION(0:7) ::  end_index   !< start particle index for current block
191
192    REAL(wp) ::  aa                 !< dummy argument for horizontal particle interpolation
193    REAL(wp) ::  bb                 !< dummy argument for horizontal particle interpolation
194    REAL(wp) ::  cc                 !< dummy argument for horizontal particle interpolation
195    REAL(wp) ::  d_sum              !< dummy argument for horizontal particle interpolation in case of topography
196    REAL(wp) ::  d_z_p_z0           !< inverse of interpolation length for logarithmic interpolation
197    REAL(wp) ::  dd                 !< dummy argument for horizontal particle interpolation
198    REAL(wp) ::  de_dx_int_l        !< x/y-interpolated TKE gradient (x) at particle position at lower vertical level
199    REAL(wp) ::  de_dx_int_u        !< x/y-interpolated TKE gradient (x) at particle position at upper vertical level
200    REAL(wp) ::  de_dy_int_l        !< x/y-interpolated TKE gradient (y) at particle position at lower vertical level
201    REAL(wp) ::  de_dy_int_u        !< x/y-interpolated TKE gradient (y) at particle position at upper vertical level
202    REAL(wp) ::  de_dt              !< temporal derivative of TKE experienced by the particle
203    REAL(wp) ::  de_dt_min          !< lower level for temporal TKE derivative
204    REAL(wp) ::  de_dz_int_l        !< x/y-interpolated TKE gradient (z) at particle position at lower vertical level
205    REAL(wp) ::  de_dz_int_u        !< x/y-interpolated TKE gradient (z) at particle position at upper vertical level
206    REAL(wp) ::  diameter           !< diamter of droplet
207    REAL(wp) ::  diss_int_l         !< x/y-interpolated dissipation at particle position at lower vertical level
208    REAL(wp) ::  diss_int_u         !< x/y-interpolated dissipation at particle position at upper vertical level
209    REAL(wp) ::  dt_gap             !< remaining time until particle time integration reaches LES time
210    REAL(wp) ::  dt_particle_m      !< previous particle time step
211    REAL(wp) ::  e_int_l            !< x/y-interpolated TKE at particle position at lower vertical level
212    REAL(wp) ::  e_int_u            !< x/y-interpolated TKE at particle position at upper vertical level
213    REAL(wp) ::  e_mean_int         !< horizontal mean TKE at particle height
214    REAL(wp) ::  exp_arg            !<
215    REAL(wp) ::  exp_term           !<
216    REAL(wp) ::  gg                 !< dummy argument for horizontal particle interpolation
217    REAL(wp) ::  height_p           !< dummy argument for logarithmic interpolation
218    REAL(wp) ::  lagr_timescale     !< Lagrangian timescale
219    REAL(wp) ::  location(1:30,1:3) !< wall locations
220    REAL(wp) ::  log_z_z0_int       !< logarithmus used for surface_layer interpolation
221    REAL(wp) ::  random_gauss       !<
222    REAL(wp) ::  RL                 !< Lagrangian autocorrelation coefficient
223    REAL(wp) ::  rg1                !< Gaussian distributed random number
224    REAL(wp) ::  rg2                !< Gaussian distributed random number
225    REAL(wp) ::  rg3                !< Gaussian distributed random number
226    REAL(wp) ::  sigma              !< velocity standard deviation
227    REAL(wp) ::  u_int_l            !< x/y-interpolated u-component at particle position at lower vertical level
228    REAL(wp) ::  u_int_u            !< x/y-interpolated u-component at particle position at upper vertical level
229    REAL(wp) ::  us_int             !< friction velocity at particle grid box
230    REAL(wp) ::  usws_int           !< surface momentum flux (u component) at particle grid box
231    REAL(wp) ::  v_int_l            !< x/y-interpolated v-component at particle position at lower vertical level
232    REAL(wp) ::  v_int_u            !< x/y-interpolated v-component at particle position at upper vertical level
233    REAL(wp) ::  vsws_int           !< surface momentum flux (u component) at particle grid box
234    REAL(wp) ::  vv_int             !<
235    REAL(wp) ::  w_int_l            !< x/y-interpolated w-component at particle position at lower vertical level
236    REAL(wp) ::  w_int_u            !< x/y-interpolated w-component at particle position at upper vertical level
237    REAL(wp) ::  w_s                !< terminal velocity of droplets
238    REAL(wp) ::  x                  !< dummy argument for horizontal particle interpolation
239    REAL(wp) ::  y                  !< dummy argument for horizontal particle interpolation
240    REAL(wp) ::  z_p                !< surface layer height (0.5 dz)
241
242    REAL(wp), PARAMETER ::  a_rog = 9.65_wp      !< parameter for fall velocity
243    REAL(wp), PARAMETER ::  b_rog = 10.43_wp     !< parameter for fall velocity
244    REAL(wp), PARAMETER ::  c_rog = 0.6_wp       !< parameter for fall velocity
245    REAL(wp), PARAMETER ::  k_cap_rog = 4.0_wp   !< parameter for fall velocity
246    REAL(wp), PARAMETER ::  k_low_rog = 12.0_wp  !< parameter for fall velocity
247    REAL(wp), PARAMETER ::  d0_rog = 0.745_wp    !< separation diameter
248
249    REAL(wp), DIMENSION(1:30) ::  d_gp_pl !< dummy argument for particle interpolation scheme in case of topography
250    REAL(wp), DIMENSION(1:30) ::  de_dxi  !< horizontal TKE gradient along x at adjacent wall
251    REAL(wp), DIMENSION(1:30) ::  de_dyi  !< horizontal TKE gradient along y at adjacent wall
252    REAL(wp), DIMENSION(1:30) ::  de_dzi  !< horizontal TKE gradient along z at adjacent wall
253    REAL(wp), DIMENSION(1:30) ::  dissi   !< dissipation at adjacent wall
254    REAL(wp), DIMENSION(1:30) ::  ei      !< TKE at adjacent wall
255
256    REAL(wp), DIMENSION(number_of_particles) ::  term_1_2     !< flag to communicate whether a particle is near topography or not
257    REAL(wp), DIMENSION(number_of_particles) ::  dens_ratio   !<
258    REAL(wp), DIMENSION(number_of_particles) ::  de_dx_int    !< horizontal TKE gradient along x at particle position
259    REAL(wp), DIMENSION(number_of_particles) ::  de_dy_int    !< horizontal TKE gradient along y at particle position
260    REAL(wp), DIMENSION(number_of_particles) ::  de_dz_int    !< horizontal TKE gradient along z at particle position
261    REAL(wp), DIMENSION(number_of_particles) ::  diss_int     !< dissipation at particle position
262    REAL(wp), DIMENSION(number_of_particles) ::  dt_particle  !< particle time step
263    REAL(wp), DIMENSION(number_of_particles) ::  e_int        !< TKE at particle position
264    REAL(wp), DIMENSION(number_of_particles) ::  fs_int       !< weighting factor for subgrid-scale particle speed
265    REAL(wp), DIMENSION(number_of_particles) ::  u_int        !< u-component of particle speed
266    REAL(wp), DIMENSION(number_of_particles) ::  v_int        !< v-component of particle speed
267    REAL(wp), DIMENSION(number_of_particles) ::  w_int        !< w-component of particle speed
268    REAL(wp), DIMENSION(number_of_particles) ::  xv           !< x-position
269    REAL(wp), DIMENSION(number_of_particles) ::  yv           !< y-position
270    REAL(wp), DIMENSION(number_of_particles) ::  zv           !< z-position
271
272    REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !< vector of Gaussian distributed random numbers
273
274    CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' )
275
276!
277!-- Determine height of Prandtl layer and distance between Prandtl-layer
278!-- height and horizontal mean roughness height, which are required for
279!-- vertical logarithmic interpolation of horizontal particle speeds
280!-- (for particles below first vertical grid level).
281    z_p      = zu(nzb+1) - zw(nzb)
282    d_z_p_z0 = 1.0_wp / ( z_p - z0_av_global )
283
284    start_index = grid_particles(kp,jp,ip)%start_index
285    end_index   = grid_particles(kp,jp,ip)%end_index
286
287    xv = particles(1:number_of_particles)%x
288    yv = particles(1:number_of_particles)%y
289    zv = particles(1:number_of_particles)%z
290
291    DO  nb = 0, 7
292!
293!--    Interpolate u velocity-component       
294       i = ip
295       j = jp + block_offset(nb)%j_off
296       k = kp + block_offset(nb)%k_off
297
298       DO  n = start_index(nb), end_index(nb)
299!
300!--       Interpolation of the u velocity component onto particle position. 
301!--       Particles are interpolation bi-linearly in the horizontal and a
302!--       linearly in the vertical. An exception is made for particles below
303!--       the first vertical grid level in case of a prandtl layer. In this
304!--       case the horizontal particle velocity components are determined using
305!--       Monin-Obukhov relations (if branch).
306!--       First, check if particle is located below first vertical grid level
307!--       above topography (Prandtl-layer height)
308!--       Determine vertical index of topography top
309          k_wall = get_topography_top_index_ji( jp, ip, 's' )
310
311          IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
312!
313!--          Resolved-scale horizontal particle velocity is zero below z0.
314             IF ( zv(n) - zw(k_wall) < z0_av_global )  THEN
315                u_int(n) = 0.0_wp
316             ELSE
317!
318!--             Determine the sublayer. Further used as index.
319                height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
320                                     * REAL( number_of_sublayers, KIND=wp )    &
321                                     * d_z_p_z0 
322!
323!--             Calculate LOG(z/z0) for exact particle height. Therefore,   
324!--             interpolate linearly between precalculated logarithm.
325                log_z_z0_int = log_z_z0(INT(height_p))                         &
326                                 + ( height_p - INT(height_p) )                &
327                                 * ( log_z_z0(INT(height_p)+1)                 &
328                                      - log_z_z0(INT(height_p))                &
329                                   ) 
330!
331!--             Get friction velocity and momentum flux from new surface data
332!--             types.
333                IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
334                     surf_def_h(0)%end_index(jp,ip) )  THEN
335                   surf_start = surf_def_h(0)%start_index(jp,ip)
336!--                Limit friction velocity. In narrow canyons or holes the
337!--                friction velocity can become very small, resulting in a too
338!--                large particle speed.
339                   us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 
340                   usws_int  = surf_def_h(0)%usws(surf_start)
341                ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                  &
342                         surf_lsm_h%end_index(jp,ip) )  THEN
343                   surf_start = surf_lsm_h%start_index(jp,ip)
344                   us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 
345                   usws_int  = surf_lsm_h%usws(surf_start)
346                ELSEIF ( surf_usm_h%start_index(jp,ip) <=                  &
347                         surf_usm_h%end_index(jp,ip) )  THEN
348                   surf_start = surf_usm_h%start_index(jp,ip)
349                   us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 
350                   usws_int  = surf_usm_h%usws(surf_start)
351                ENDIF
352
353!
354!--             Neutral solution is applied for all situations, e.g. also for
355!--             unstable and stable situations. Even though this is not exact
356!--             this saves a lot of CPU time since several calls of intrinsic
357!--             FORTRAN procedures (LOG, ATAN) are avoided, This is justified
358!--             as sensitivity studies revealed no significant effect of
359!--             using the neutral solution also for un/stable situations.
360                u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp )           & 
361                            * log_z_z0_int - u_gtrans
362               
363             ENDIF
364!
365!--       Particle above the first grid level. Bi-linear interpolation in the
366!--       horizontal and linear interpolation in the vertical direction.
367          ELSE
368
369             x  = xv(n) + ( 0.5_wp - i ) * dx
370             y  = yv(n) - j * dy
371             aa = x**2          + y**2
372             bb = ( dx - x )**2 + y**2
373             cc = x**2          + ( dy - y )**2
374             dd = ( dx - x )**2 + ( dy - y )**2
375             gg = aa + bb + cc + dd
376
377             u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
378                         + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) *            &
379                         u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
380
381             IF ( k == nzt )  THEN
382                u_int(n) = u_int_l
383             ELSE
384                u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1)  &
385                            + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) *           &
386                            u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
387                u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dz *                  &
388                           ( u_int_u - u_int_l )
389             ENDIF
390
391          ENDIF
392
393       ENDDO
394!
395!--    Same procedure for interpolation of the v velocity-component
396       i = ip + block_offset(nb)%i_off
397       j = jp
398       k = kp + block_offset(nb)%k_off
399
400       DO  n = start_index(nb), end_index(nb)
401
402!
403!--       Determine vertical index of topography top
404          k_wall = get_topography_top_index_ji( jp,ip, 's' )
405
406          IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
407             IF ( zv(n) - zw(k_wall) < z0_av_global )  THEN
408!
409!--             Resolved-scale horizontal particle velocity is zero below z0.
410                v_int(n) = 0.0_wp
411             ELSE       
412!
413!--             Determine the sublayer. Further used as index. Please note,
414!--             logarithmus can not be reused from above, as in in case of
415!--             topography particle on u-grid can be above surface-layer height,
416!--             whereas it can be below on v-grid.
417                height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
418                                  * REAL( number_of_sublayers, KIND=wp )       &
419                                  * d_z_p_z0 
420!
421!--             Calculate LOG(z/z0) for exact particle height. Therefore,   
422!--             interpolate linearly between precalculated logarithm.
423                log_z_z0_int = log_z_z0(INT(height_p))                         &
424                                 + ( height_p - INT(height_p) )                &
425                                 * ( log_z_z0(INT(height_p)+1)                 &
426                                      - log_z_z0(INT(height_p))                &
427                                   ) 
428!
429!--             Get friction velocity and momentum flux from new surface data
430!--             types.
431                IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
432                     surf_def_h(0)%end_index(jp,ip) )  THEN
433                   surf_start = surf_def_h(0)%start_index(jp,ip)
434!--                Limit friction velocity. In narrow canyons or holes the
435!--                friction velocity can become very small, resulting in a too
436!--                large particle speed.
437                   us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 
438                   vsws_int  = surf_def_h(0)%vsws(surf_start)
439                ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                  &
440                         surf_lsm_h%end_index(jp,ip) )  THEN
441                   surf_start = surf_lsm_h%start_index(jp,ip)
442                   us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 
443                   vsws_int  = surf_lsm_h%vsws(surf_start)
444                ELSEIF ( surf_usm_h%start_index(jp,ip) <=                  &
445                         surf_usm_h%end_index(jp,ip) )  THEN
446                   surf_start = surf_usm_h%start_index(jp,ip)
447                   us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 
448                   vsws_int  = surf_usm_h%vsws(surf_start)
449                ENDIF 
450!
451!--             Neutral solution is applied for all situations, e.g. also for
452!--             unstable and stable situations. Even though this is not exact
453!--             this saves a lot of CPU time since several calls of intrinsic
454!--             FORTRAN procedures (LOG, ATAN) are avoided, This is justified
455!--             as sensitivity studies revealed no significant effect of
456!--             using the neutral solution also for un/stable situations.
457                v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp )           &
458                         * log_z_z0_int - v_gtrans
459
460             ENDIF
461
462          ELSE
463             x  = xv(n) - i * dx
464             y  = yv(n) + ( 0.5_wp - j ) * dy
465             aa = x**2          + y**2
466             bb = ( dx - x )**2 + y**2
467             cc = x**2          + ( dy - y )**2
468             dd = ( dx - x )**2 + ( dy - y )**2
469             gg = aa + bb + cc + dd
470
471             v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
472                       + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
473                       ) / ( 3.0_wp * gg ) - v_gtrans
474
475             IF ( k == nzt )  THEN
476                v_int(n) = v_int_l
477             ELSE
478                v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
479                          + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
480                          ) / ( 3.0_wp * gg ) - v_gtrans
481                v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dz * &
482                                  ( v_int_u - v_int_l )
483             ENDIF
484
485          ENDIF
486
487       ENDDO
488!
489!--    Same procedure for interpolation of the w velocity-component
490       i = ip + block_offset(nb)%i_off
491       j = jp + block_offset(nb)%j_off
492       k = kp - 1
493
494       DO  n = start_index(nb), end_index(nb)
495
496          IF ( vertical_particle_advection(particles(n)%group) )  THEN
497
498             x  = xv(n) - i * dx
499             y  = yv(n) - j * dy
500             aa = x**2          + y**2
501             bb = ( dx - x )**2 + y**2
502             cc = x**2          + ( dy - y )**2
503             dd = ( dx - x )**2 + ( dy - y )**2
504             gg = aa + bb + cc + dd
505
506             w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
507                       + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
508                       ) / ( 3.0_wp * gg )
509
510             IF ( k == nzt )  THEN
511                w_int(n) = w_int_l
512             ELSE
513                w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + &
514                            ( gg-bb ) * w(k+1,j,i+1) + &
515                            ( gg-cc ) * w(k+1,j+1,i) + &
516                            ( gg-dd ) * w(k+1,j+1,i+1) &
517                          ) / ( 3.0_wp * gg )
518                w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dz * &
519                           ( w_int_u - w_int_l )
520             ENDIF
521
522          ELSE
523
524             w_int(n) = 0.0_wp
525
526          ENDIF
527
528       ENDDO
529
530    ENDDO
531
532!-- Interpolate and calculate quantities needed for calculating the SGS
533!-- velocities
534    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
535   
536       DO  nb = 0,7
537         
538          subbox_at_wall = .FALSE.         
539!
540!--       In case of topography check if subbox is adjacent to a wall
541          IF ( .NOT. topography == 'flat' ) THEN
542             i = ip + MERGE( -1_iwp , 1_iwp, BTEST( nb, 2 ) )
543             j = jp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 1 ) )
544             k = kp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 0 ) )
545             IF ( .NOT. BTEST(wall_flags_0(k,  jp, ip), 0) .OR.                &
546                  .NOT. BTEST(wall_flags_0(kp, j,  ip), 0) .OR.                &
547                  .NOT. BTEST(wall_flags_0(kp, jp, i ), 0) )                   &
548             THEN
549                subbox_at_wall = .TRUE.
550             ENDIF
551          ENDIF
552          IF ( subbox_at_wall ) THEN
553             e_int(start_index(nb):end_index(nb))     = e(kp,jp,ip) 
554             diss_int(start_index(nb):end_index(nb))  = diss(kp,jp,ip)
555             de_dx_int(start_index(nb):end_index(nb)) = de_dx(kp,jp,ip)
556             de_dy_int(start_index(nb):end_index(nb)) = de_dy(kp,jp,ip)
557             de_dz_int(start_index(nb):end_index(nb)) = de_dz(kp,jp,ip)
558!
559!--          Set flag for stochastic equation.
560             term_1_2(start_index(nb):end_index(nb)) = 0.0_wp             
561          ELSE
562             i = ip + block_offset(nb)%i_off
563             j = jp + block_offset(nb)%j_off
564             k = kp + block_offset(nb)%k_off
565
566             DO  n = start_index(nb), end_index(nb)
567!
568!--             Interpolate TKE
569                x  = xv(n) - i * dx
570                y  = yv(n) - j * dy
571                aa = x**2          + y**2
572                bb = ( dx - x )**2 + y**2
573                cc = x**2          + ( dy - y )**2
574                dd = ( dx - x )**2 + ( dy - y )**2
575                gg = aa + bb + cc + dd
576
577                e_int_l = ( ( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1)   &
578                          + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
579                          ) / ( 3.0_wp * gg )
580
581                IF ( k+1 == nzt+1 )  THEN
582                   e_int(n) = e_int_l
583                ELSE
584                   e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
585                               ( gg - bb ) * e(k+1,j,i+1) + &
586                               ( gg - cc ) * e(k+1,j+1,i) + &
587                               ( gg - dd ) * e(k+1,j+1,i+1) &
588                            ) / ( 3.0_wp * gg )
589                   e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
590                                     ( e_int_u - e_int_l )
591                ENDIF
592!
593!--             Needed to avoid NaN particle velocities (this might not be
594!--             required any more)
595                IF ( e_int(n) <= 0.0_wp )  THEN
596                   e_int(n) = 1.0E-20_wp
597                ENDIF
598!
599!--             Interpolate the TKE gradient along x (adopt incides i,j,k and
600!--             all position variables from above (TKE))
601                de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
602                                ( gg - bb ) * de_dx(k,j,i+1) + &
603                                ( gg - cc ) * de_dx(k,j+1,i) + &
604                                ( gg - dd ) * de_dx(k,j+1,i+1) &
605                               ) / ( 3.0_wp * gg )
606
607                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
608                   de_dx_int(n) = de_dx_int_l
609                ELSE
610                   de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
611                                   ( gg - bb ) * de_dx(k+1,j,i+1) + &
612                                   ( gg - cc ) * de_dx(k+1,j+1,i) + &
613                                   ( gg - dd ) * de_dx(k+1,j+1,i+1) &
614                                  ) / ( 3.0_wp * gg )
615                   de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / dz * &
616                                              ( de_dx_int_u - de_dx_int_l )
617                ENDIF
618!
619!--             Interpolate the TKE gradient along y
620                de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
621                                ( gg - bb ) * de_dy(k,j,i+1) + &
622                                ( gg - cc ) * de_dy(k,j+1,i) + &
623                                ( gg - dd ) * de_dy(k,j+1,i+1) &
624                               ) / ( 3.0_wp * gg )
625                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
626                   de_dy_int(n) = de_dy_int_l
627                ELSE
628                   de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
629                                   ( gg - bb ) * de_dy(k+1,j,i+1) + &
630                                   ( gg - cc ) * de_dy(k+1,j+1,i) + &
631                                   ( gg - dd ) * de_dy(k+1,j+1,i+1) &
632                                  ) / ( 3.0_wp * gg )
633                      de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * &
634                                                 ( de_dy_int_u - de_dy_int_l )
635                ENDIF
636
637!
638!--             Interpolate the TKE gradient along z
639                IF ( zv(n) < 0.5_wp * dz )  THEN
640                   de_dz_int(n) = 0.0_wp
641                ELSE
642                   de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
643                                   ( gg - bb ) * de_dz(k,j,i+1) + &
644                                   ( gg - cc ) * de_dz(k,j+1,i) + &
645                                   ( gg - dd ) * de_dz(k,j+1,i+1) &
646                                  ) / ( 3.0_wp * gg )
647
648                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
649                      de_dz_int(n) = de_dz_int_l
650                   ELSE
651                      de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
652                                      ( gg - bb ) * de_dz(k+1,j,i+1) + &
653                                      ( gg - cc ) * de_dz(k+1,j+1,i) + &
654                                      ( gg - dd ) * de_dz(k+1,j+1,i+1) &
655                                     ) / ( 3.0_wp * gg )
656                      de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) / dz * &
657                                                 ( de_dz_int_u - de_dz_int_l )
658                   ENDIF
659                ENDIF
660
661!
662!--             Interpolate the dissipation of TKE
663                diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
664                               ( gg - bb ) * diss(k,j,i+1) + &
665                               ( gg - cc ) * diss(k,j+1,i) + &
666                               ( gg - dd ) * diss(k,j+1,i+1) &
667                               ) / ( 3.0_wp * gg )
668
669                IF ( k == nzt )  THEN
670                   diss_int(n) = diss_int_l
671                ELSE
672                   diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
673                                  ( gg - bb ) * diss(k+1,j,i+1) + &
674                                  ( gg - cc ) * diss(k+1,j+1,i) + &
675                                  ( gg - dd ) * diss(k+1,j+1,i+1) &
676                                 ) / ( 3.0_wp * gg )
677                   diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz * &
678                                            ( diss_int_u - diss_int_l )
679                ENDIF
680
681!
682!--             Set flag for stochastic equation.
683                term_1_2(n) = 1.0_wp
684             ENDDO
685          ENDIF
686       ENDDO
687
688       DO nb = 0,7
689          i = ip + block_offset(nb)%i_off
690          j = jp + block_offset(nb)%j_off
691          k = kp + block_offset(nb)%k_off
692
693          DO  n = start_index(nb), end_index(nb)
694!
695!--          Vertical interpolation of the horizontally averaged SGS TKE and
696!--          resolved-scale velocity variances and use the interpolated values
697!--          to calculate the coefficient fs, which is a measure of the ratio
698!--          of the subgrid-scale turbulent kinetic energy to the total amount
699!--          of turbulent kinetic energy.
700             IF ( k == 0 )  THEN
701                e_mean_int = hom(0,1,8,0)
702             ELSE
703                e_mean_int = hom(k,1,8,0) +                                    &
704                                           ( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
705                                           ( zu(k+1) - zu(k) ) *               &
706                                           ( zv(n) - zu(k) )
707             ENDIF
708
709             kw = kp - 1
710
711             IF ( k == 0 )  THEN
712                aa  = hom(k+1,1,30,0)  * ( zv(n) / &
713                                         ( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
714                bb  = hom(k+1,1,31,0)  * ( zv(n) / &
715                                         ( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
716                cc  = hom(kw+1,1,32,0) * ( zv(n) / &
717                                         ( 1.0_wp * ( zw(kw+1) - zw(kw) ) ) )
718             ELSE
719                aa  = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) *    &
720                           ( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
721                bb  = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) *    &
722                           ( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
723                cc  = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) *   &
724                           ( ( zv(n) - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
725             ENDIF
726
727             vv_int = ( 1.0_wp / 3.0_wp ) * ( aa + bb + cc )
728!
729!--          Needed to avoid NaN particle velocities. The value of 1.0 is just
730!--          an educated guess for the given case.
731             IF ( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int == 0.0_wp )  THEN
732                fs_int(n) = 1.0_wp
733             ELSE
734                fs_int(n) = ( 2.0_wp / 3.0_wp ) * e_mean_int /                 &
735                            ( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int )
736             ENDIF
737
738          ENDDO
739       ENDDO
740
741       DO  nb = 0, 7
742          DO  n = start_index(nb), end_index(nb)
743             rg(n,1) = random_gauss( iran_part, 5.0_wp )
744             rg(n,2) = random_gauss( iran_part, 5.0_wp )
745             rg(n,3) = random_gauss( iran_part, 5.0_wp )
746          ENDDO
747       ENDDO
748
749       DO  nb = 0, 7
750          DO  n = start_index(nb), end_index(nb)
751
752!
753!--          Calculate the Lagrangian timescale according to Weil et al. (2004).
754             lagr_timescale = ( 4.0_wp * e_int(n) + 1E-20_wp ) / &
755                              ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp )
756
757!
758!--          Calculate the next particle timestep. dt_gap is the time needed to
759!--          complete the current LES timestep.
760             dt_gap = dt_3d - particles(n)%dt_sum
761             dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale, dt_gap )
762
763!
764!--          The particle timestep should not be too small in order to prevent
765!--          the number of particle timesteps of getting too large
766             IF ( dt_particle(n) < dt_min_part  .AND.  dt_min_part < dt_gap )  THEN
767                dt_particle(n) = dt_min_part
768             ENDIF
769
770!
771!--          Calculate the SGS velocity components
772             IF ( particles(n)%age == 0.0_wp )  THEN
773!
774!--             For new particles the SGS components are derived from the SGS
775!--             TKE. Limit the Gaussian random number to the interval
776!--             [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities
777!--             from becoming unrealistically large.
778                particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf_part * e_int(n)     &
779                                          + 1E-20_wp ) *                       &
780                                     ( rg(n,1) - 1.0_wp )
781                particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf_part * e_int(n)     &
782                                          + 1E-20_wp ) *                       &
783                                     ( rg(n,2) - 1.0_wp )
784                particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf_part * e_int(n)     &
785                                          + 1E-20_wp ) *                       &
786                                     ( rg(n,3) - 1.0_wp )
787
788             ELSE
789!
790!--             Restriction of the size of the new timestep: compared to the
791!--             previous timestep the increase must not exceed 200%. First,
792!--             check if age > age_m, in order to prevent that particles get zero
793!--             timestep.
794                dt_particle_m = MERGE( dt_particle(n),                         &
795                                       particles(n)%age - particles(n)%age_m,  &
796                                       particles(n)%age - particles(n)%age_m < &
797                                       1E-8_wp )
798                IF ( dt_particle(n) > 2.0_wp * dt_particle_m )  THEN
799                   dt_particle(n) = 2.0_wp * dt_particle_m
800                ENDIF
801
802!
803!--             For old particles the SGS components are correlated with the
804!--             values from the previous timestep. Random numbers have also to
805!--             be limited (see above).
806!--             As negative values for the subgrid TKE are not allowed, the
807!--             change of the subgrid TKE with time cannot be smaller than
808!--             -e_int(n)/dt_particle. This value is used as a lower boundary
809!--             value for the change of TKE
810                de_dt_min = - e_int(n) / dt_particle(n)
811
812                de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m
813
814                IF ( de_dt < de_dt_min )  THEN
815                   de_dt = de_dt_min
816                ENDIF
817
818                CALL weil_stochastic_eq(particles(n)%rvar1, fs_int(n), e_int(n),& 
819                                        de_dx_int(n), de_dt, diss_int(n),       &
820                                        dt_particle(n), rg(n,1), term_1_2(n) )
821
822                CALL weil_stochastic_eq(particles(n)%rvar2, fs_int(n), e_int(n),& 
823                                        de_dy_int(n), de_dt, diss_int(n),       &
824                                        dt_particle(n), rg(n,2), term_1_2(n) )
825
826                CALL weil_stochastic_eq(particles(n)%rvar3, fs_int(n), e_int(n),& 
827                                        de_dz_int(n), de_dt, diss_int(n),       &
828                                        dt_particle(n), rg(n,3), term_1_2(n) )
829
830             ENDIF
831
832             u_int(n) = u_int(n) + particles(n)%rvar1
833             v_int(n) = v_int(n) + particles(n)%rvar2
834             w_int(n) = w_int(n) + particles(n)%rvar3
835!
836!--          Store the SGS TKE of the current timelevel which is needed for
837!--          for calculating the SGS particle velocities at the next timestep
838             particles(n)%e_m = e_int(n)
839          ENDDO
840       ENDDO
841
842    ELSE
843!
844!--    If no SGS velocities are used, only the particle timestep has to
845!--    be set
846       dt_particle = dt_3d
847
848    ENDIF
849
850    dens_ratio = particle_groups(particles(1:number_of_particles)%group)%density_ratio
851
852    IF ( ANY( dens_ratio == 0.0_wp ) )  THEN
853       DO  nb = 0, 7
854          DO  n = start_index(nb), end_index(nb)
855
856!
857!--          Particle advection
858             IF ( dens_ratio(n) == 0.0_wp )  THEN
859!
860!--             Pure passive transport (without particle inertia)
861                particles(n)%x = xv(n) + u_int(n) * dt_particle(n)
862                particles(n)%y = yv(n) + v_int(n) * dt_particle(n)
863                particles(n)%z = zv(n) + w_int(n) * dt_particle(n)
864
865                particles(n)%speed_x = u_int(n)
866                particles(n)%speed_y = v_int(n)
867                particles(n)%speed_z = w_int(n)
868
869             ELSE
870!
871!--             Transport of particles with inertia
872                particles(n)%x = particles(n)%x + particles(n)%speed_x * &
873                                                  dt_particle(n)
874                particles(n)%y = particles(n)%y + particles(n)%speed_y * &
875                                                  dt_particle(n)
876                particles(n)%z = particles(n)%z + particles(n)%speed_z * &
877                                                  dt_particle(n)
878
879!
880!--             Update of the particle velocity
881                IF ( cloud_droplets )  THEN
882!
883!--                Terminal velocity is computed for vertical direction (Rogers et
884!--                al., 1993, J. Appl. Meteorol.)
885                   diameter = particles(n)%radius * 2000.0_wp !diameter in mm
886                   IF ( diameter <= d0_rog )  THEN
887                      w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
888                   ELSE
889                      w_s = a_rog - b_rog * EXP( -c_rog * diameter )
890                   ENDIF
891
892!
893!--                If selected, add random velocities following Soelch and Kaercher
894!--                (2010, Q. J. R. Meteorol. Soc.)
895                   IF ( use_sgs_for_particles )  THEN
896                      lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
897                      RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
898                      sigma          = SQRT( e(kp,jp,ip) )
899
900                      rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
901                      rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
902                      rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
903
904                      particles(n)%rvar1 = RL * particles(n)%rvar1 +              &
905                                           SQRT( 1.0_wp - RL**2 ) * sigma * rg1
906                      particles(n)%rvar2 = RL * particles(n)%rvar2 +              &
907                                           SQRT( 1.0_wp - RL**2 ) * sigma * rg2
908                      particles(n)%rvar3 = RL * particles(n)%rvar3 +              &
909                                           SQRT( 1.0_wp - RL**2 ) * sigma * rg3
910
911                      particles(n)%speed_x = u_int(n) + particles(n)%rvar1
912                      particles(n)%speed_y = v_int(n) + particles(n)%rvar2
913                      particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
914                   ELSE
915                      particles(n)%speed_x = u_int(n)
916                      particles(n)%speed_y = v_int(n)
917                      particles(n)%speed_z = w_int(n) - w_s
918                   ENDIF
919
920                ELSE
921
922                   IF ( use_sgs_for_particles )  THEN
923                      exp_arg  = particle_groups(particles(n)%group)%exp_arg
924                      exp_term = EXP( -exp_arg * dt_particle(n) )
925                   ELSE
926                      exp_arg  = particle_groups(particles(n)%group)%exp_arg
927                      exp_term = particle_groups(particles(n)%group)%exp_term
928                   ENDIF
929                   particles(n)%speed_x = particles(n)%speed_x * exp_term +         &
930                                          u_int(n) * ( 1.0_wp - exp_term )
931                   particles(n)%speed_y = particles(n)%speed_y * exp_term +         &
932                                          v_int(n) * ( 1.0_wp - exp_term )
933                   particles(n)%speed_z = particles(n)%speed_z * exp_term +         &
934                                          ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * &
935                                          g / exp_arg ) * ( 1.0_wp - exp_term )
936                ENDIF
937
938             ENDIF
939          ENDDO
940       ENDDO
941   
942    ELSE
943
944       DO  nb = 0, 7
945          DO  n = start_index(nb), end_index(nb)
946!
947!--          Transport of particles with inertia
948             particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)
949             particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)
950             particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)
951!
952!--          Update of the particle velocity
953             IF ( cloud_droplets )  THEN
954!
955!--             Terminal velocity is computed for vertical direction (Rogers et al.,
956!--             1993, J. Appl. Meteorol.)
957                diameter = particles(n)%radius * 2000.0_wp !diameter in mm
958                IF ( diameter <= d0_rog )  THEN
959                   w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
960                ELSE
961                   w_s = a_rog - b_rog * EXP( -c_rog * diameter )
962                ENDIF
963
964!
965!--             If selected, add random velocities following Soelch and Kaercher
966!--             (2010, Q. J. R. Meteorol. Soc.)
967                IF ( use_sgs_for_particles )  THEN
968                    lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
969                    RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
970                    sigma          = SQRT( e(kp,jp,ip) )
971
972                    rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
973                    rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
974                    rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
975
976                    particles(n)%rvar1 = RL * particles(n)%rvar1 +                &
977                                         SQRT( 1.0_wp - RL**2 ) * sigma * rg1
978                    particles(n)%rvar2 = RL * particles(n)%rvar2 +                &
979                                         SQRT( 1.0_wp - RL**2 ) * sigma * rg2
980                    particles(n)%rvar3 = RL * particles(n)%rvar3 +                &
981                                         SQRT( 1.0_wp - RL**2 ) * sigma * rg3
982
983                    particles(n)%speed_x = u_int(n) + particles(n)%rvar1
984                    particles(n)%speed_y = v_int(n) + particles(n)%rvar2
985                    particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
986                ELSE
987                    particles(n)%speed_x = u_int(n)
988                    particles(n)%speed_y = v_int(n)
989                    particles(n)%speed_z = w_int(n) - w_s
990                ENDIF
991
992             ELSE
993
994                IF ( use_sgs_for_particles )  THEN
995                   exp_arg  = particle_groups(particles(n)%group)%exp_arg
996                   exp_term = EXP( -exp_arg * dt_particle(n) )
997                ELSE
998                   exp_arg  = particle_groups(particles(n)%group)%exp_arg
999                   exp_term = particle_groups(particles(n)%group)%exp_term
1000                ENDIF
1001                particles(n)%speed_x = particles(n)%speed_x * exp_term +             &
1002                                       u_int(n) * ( 1.0_wp - exp_term )
1003                particles(n)%speed_y = particles(n)%speed_y * exp_term +             &
1004                                       v_int(n) * ( 1.0_wp - exp_term )
1005                particles(n)%speed_z = particles(n)%speed_z * exp_term +             &
1006                                       ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &
1007                                       exp_arg ) * ( 1.0_wp - exp_term )
1008             ENDIF
1009          ENDDO
1010       ENDDO
1011
1012    ENDIF
1013
1014!
1015!-- Store the old age of the particle ( needed to prevent that a
1016!-- particle crosses several PEs during one timestep, and for the
1017!-- evaluation of the subgrid particle velocity fluctuations )
1018    particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age
1019
1020    DO  nb = 0, 7
1021       DO  n = start_index(nb), end_index(nb)
1022!
1023!--       Increment the particle age and the total time that the particle
1024!--       has advanced within the particle timestep procedure
1025          particles(n)%age    = particles(n)%age    + dt_particle(n)
1026          particles(n)%dt_sum = particles(n)%dt_sum + dt_particle(n)
1027
1028!
1029!--       Check whether there is still a particle that has not yet completed
1030!--       the total LES timestep
1031          IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp )  THEN
1032             dt_3d_reached_l = .FALSE.
1033          ENDIF
1034
1035       ENDDO
1036    ENDDO
1037
1038    CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
1039
1040
1041 END SUBROUTINE lpm_advec
1042
1043! Description:
1044! ------------
1045!> Calculation of subgrid-scale particle speed using the stochastic model
1046!> of Weil et al. (2004, JAS, 61, 2877-2887).
1047!------------------------------------------------------------------------------!
1048 SUBROUTINE weil_stochastic_eq( v_sgs, fs_n, e_n, dedxi_n, dedt_n, diss_n,     &
1049                                dt_n, rg_n, fac )
1050
1051    USE kinds
1052
1053    USE particle_attributes,                                                   &
1054        ONLY:  c_0, sgs_wf_part
1055
1056    IMPLICIT NONE
1057
1058    REAL(wp) ::  a1      !< dummy argument
1059    REAL(wp) ::  dedt_n  !< time derivative of TKE at particle position
1060    REAL(wp) ::  dedxi_n !< horizontal derivative of TKE at particle position
1061    REAL(wp) ::  diss_n  !< dissipation at particle position
1062    REAL(wp) ::  dt_n    !< particle timestep
1063    REAL(wp) ::  e_n     !< TKE at particle position
1064    REAL(wp) ::  fac     !< flag to identify adjacent topography
1065    REAL(wp) ::  fs_n    !< weighting factor to prevent that subgrid-scale particle speed becomes too large
1066    REAL(wp) ::  sgs_w   !< constant (1/3)
1067    REAL(wp) ::  rg_n    !< random number
1068    REAL(wp) ::  term1   !< memory term
1069    REAL(wp) ::  term2   !< drift correction term
1070    REAL(wp) ::  term3   !< random term
1071    REAL(wp) ::  v_sgs   !< subgrid-scale velocity component
1072
1073!-- At first, limit TKE to a small non-zero number, in order to prevent
1074!-- the occurrence of extremely large SGS-velocities in case TKE is zero,
1075!-- (could occur at the simulation begin).
1076    e_n = MAX( e_n, 1E-20_wp )
1077!
1078!-- Please note, terms 1 and 2 (drift and memory term, respectively) are
1079!-- multiplied by a flag to switch of both terms near topography.
1080!-- This is necessary, as both terms may cause a subgrid-scale velocity build up
1081!-- if particles are trapped in regions with very small TKE, e.g. in narrow street
1082!-- canyons resolved by only a few grid points. Hence, term 1 and term 2 are
1083!-- disabled if one of the adjacent grid points belongs to topography.
1084!-- Moreover, in this case, the  previous subgrid-scale component is also set
1085!-- to zero.
1086
1087    a1 = fs_n * c_0 * diss_n
1088!
1089!-- Memory term
1090    term1 = - a1 * v_sgs * dt_n / ( 4.0_wp * sgs_wf_part * e_n + 1E-20_wp )    &
1091                 * fac
1092!
1093!-- Drift correction term
1094    term2 = ( ( dedt_n * v_sgs / e_n ) + dedxi_n ) * 0.5_wp * dt_n              &
1095                 * fac
1096!
1097!-- Random term
1098    term3 = SQRT( MAX( a1, 1E-20 ) ) * ( rg_n - 1.0_wp ) * SQRT( dt_n )
1099!
1100!-- In cese one of the adjacent grid-boxes belongs to topograhy, the previous
1101!-- subgrid-scale velocity component is set to zero, in order to prevent a
1102!-- velocity build-up.
1103!-- This case, set also previous subgrid-scale component to zero.
1104    v_sgs = v_sgs * fac + term1 + term2 + term3
1105
1106 END SUBROUTINE weil_stochastic_eq
Note: See TracBrowser for help on using the repository browser.