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

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

changes from last commit documented

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