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

Last change on this file since 3655 was 3655, checked in by knoop, 3 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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