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

Last change on this file since 3377 was 3274, checked in by knoop, 6 years ago

Modularization of all bulk cloud physics code components

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