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

Last change on this file since 2969 was 2969, checked in by thiele, 6 years ago

Bugfix in Interpolation indices

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