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

Last change on this file since 2749 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

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