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

Last change on this file since 1685 was 1685, checked in by raasch, 6 years ago

bugfix concerning vertical index calculation for particles in case of ocean runs

  • Property svn:keywords set to Id
File size: 61.4 KB
Line 
1!> @file lpm_advec.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2014 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
21! TKE check for negative values (so far, only zero value was checked)
22! offset_ocean_nzt_m1 removed
23!
24! Former revisions:
25! -----------------
26! $Id: lpm_advec.f90 1685 2015-10-08 07:32:13Z raasch $
27!
28! 1682 2015-10-07 23:56:08Z knoop
29! Code annotations made doxygen readable
30!
31! 1583 2015-04-15 12:16:27Z suehring
32! Bugfix: particle advection within Prandtl-layer in case of Galilei
33! transformation.
34!
35! 1369 2014-04-24 05:57:38Z raasch
36! usage of module interfaces removed
37!
38! 1359 2014-04-11 17:15:14Z hoffmann
39! New particle structure integrated.
40! Kind definition added to all floating point numbers.
41!
42! 1322 2014-03-20 16:38:49Z raasch
43! REAL constants defined as wp_kind
44!
45! 1320 2014-03-20 08:40:49Z raasch
46! ONLY-attribute added to USE-statements,
47! kind-parameters added to all INTEGER and REAL declaration statements,
48! kinds are defined in new module kinds,
49! revision history before 2012 removed,
50! comment fields (!:) to be used for variable explanations added to
51! all variable declaration statements
52!
53! 1314 2014-03-14 18:25:17Z suehring
54! Vertical logarithmic interpolation of horizontal particle speed for particles
55! between roughness height and first vertical grid level.
56!
57! 1036 2012-10-22 13:43:42Z raasch
58! code put under GPL (PALM 3.9)
59!
60! 849 2012-03-15 10:35:09Z raasch
61! initial revision (former part of advec_particles)
62!
63!
64! Description:
65! ------------
66!> Calculation of new particle positions due to advection using a simple Euler
67!> scheme. Particles may feel inertia effects. SGS transport can be included
68!> using the stochastic model of Weil et al. (2004, JAS, 61, 2877-2887).
69!------------------------------------------------------------------------------!
70 SUBROUTINE lpm_advec (ip,jp,kp)
71 
72
73    USE arrays_3d,                                                             &
74        ONLY:  de_dx, de_dy, de_dz, diss, e, u, us, usws, v, vsws, w, z0, zu,  &
75               zw
76
77    USE cpulog
78
79    USE pegrid
80
81    USE control_parameters,                                                    &
82        ONLY:  atmos_ocean_sign, cloud_droplets, dt_3d, dt_3d_reached_l, dz,   &
83               g, kappa, molecular_viscosity, prandtl_layer, topography,       &
84               u_gtrans, v_gtrans, simulated_time
85
86    USE grid_variables,                                                        &
87        ONLY:  ddx, dx, ddy, dy
88       
89    USE indices,                                                               &
90        ONLY:  nzb, nzb_s_inner, nzt
91       
92    USE kinds
93   
94    USE particle_attributes,                                                   &
95        ONLY:  block_offset, c_0, density_ratio, dt_min_part, grid_particles,  &
96               iran_part, log_z_z0, number_of_particles, number_of_sublayers,  &
97               particles, particle_groups, offset_ocean_nzt, sgs_wfu_part,     &
98               sgs_wfv_part, sgs_wfw_part, use_sgs_for_particles,              &
99               vertical_particle_advection, z0_av_global
100       
101    USE statistics,                                                            &
102        ONLY:  hom
103
104    IMPLICIT NONE
105
106    INTEGER(iwp) ::  agp                         !<
107    INTEGER(iwp) ::  gp_outside_of_building(1:8) !<
108    INTEGER(iwp) ::  i                           !<
109    INTEGER(iwp) ::  ip                          !<
110    INTEGER(iwp) ::  j                           !<
111    INTEGER(iwp) ::  jp                          !<
112    INTEGER(iwp) ::  k                           !<
113    INTEGER(iwp) ::  kp                          !<
114    INTEGER(iwp) ::  kw                          !<
115    INTEGER(iwp) ::  n                           !<
116    INTEGER(iwp) ::  nb                          !<
117    INTEGER(iwp) ::  num_gp                      !<
118
119    INTEGER(iwp), DIMENSION(0:7) ::  start_index !<
120    INTEGER(iwp), DIMENSION(0:7) ::  end_index   !<
121
122    REAL(wp) ::  aa                 !<
123    REAL(wp) ::  bb                 !<
124    REAL(wp) ::  cc                 !<
125    REAL(wp) ::  d_sum              !<
126    REAL(wp) ::  d_z_p_z0           !<
127    REAL(wp) ::  dd                 !<
128    REAL(wp) ::  de_dx_int_l        !<
129    REAL(wp) ::  de_dx_int_u        !<
130    REAL(wp) ::  de_dy_int_l        !<
131    REAL(wp) ::  de_dy_int_u        !<
132    REAL(wp) ::  de_dt              !<
133    REAL(wp) ::  de_dt_min          !<
134    REAL(wp) ::  de_dz_int_l        !<
135    REAL(wp) ::  de_dz_int_u        !<
136    REAL(wp) ::  diss_int_l         !<
137    REAL(wp) ::  diss_int_u         !<
138    REAL(wp) ::  dt_gap             !<
139    REAL(wp) ::  dt_particle_m      !<
140    REAL(wp) ::  e_int_l            !<
141    REAL(wp) ::  e_int_u            !<
142    REAL(wp) ::  e_mean_int         !<
143    REAL(wp) ::  exp_arg            !<
144    REAL(wp) ::  exp_term           !<
145    REAL(wp) ::  gg                 !<
146    REAL(wp) ::  height_int         !<
147    REAL(wp) ::  height_p           !<
148    REAL(wp) ::  lagr_timescale     !<
149    REAL(wp) ::  location(1:30,1:3) !<
150    REAL(wp) ::  random_gauss       !<
151    REAL(wp) ::  u_int_l            !<
152    REAL(wp) ::  u_int_u            !<
153    REAL(wp) ::  us_int             !<
154    REAL(wp) ::  v_int_l            !<
155    REAL(wp) ::  v_int_u            !<
156    REAL(wp) ::  vv_int             !<
157    REAL(wp) ::  w_int_l            !<
158    REAL(wp) ::  w_int_u            !<
159    REAL(wp) ::  x                  !<
160    REAL(wp) ::  y                  !<
161    REAL(wp) ::  z_p                !<   
162
163    REAL(wp), DIMENSION(1:30) ::  d_gp_pl !<
164    REAL(wp), DIMENSION(1:30) ::  de_dxi  !<
165    REAL(wp), DIMENSION(1:30) ::  de_dyi  !<
166    REAL(wp), DIMENSION(1:30) ::  de_dzi  !<
167    REAL(wp), DIMENSION(1:30) ::  dissi   !<
168    REAL(wp), DIMENSION(1:30) ::  ei      !<
169
170    REAL(wp), DIMENSION(number_of_particles) ::  dens_ratio   !<
171    REAL(wp), DIMENSION(number_of_particles) ::  de_dx_int    !<
172    REAL(wp), DIMENSION(number_of_particles) ::  de_dy_int    !<
173    REAL(wp), DIMENSION(number_of_particles) ::  de_dz_int    !<
174    REAL(wp), DIMENSION(number_of_particles) ::  diss_int     !<
175    REAL(wp), DIMENSION(number_of_particles) ::  dt_particle  !<
176    REAL(wp), DIMENSION(number_of_particles) ::  e_int        !<
177    REAL(wp), DIMENSION(number_of_particles) ::  fs_int       !<
178    REAL(wp), DIMENSION(number_of_particles) ::  log_z_z0_int !<
179    REAL(wp), DIMENSION(number_of_particles) ::  u_int        !<
180    REAL(wp), DIMENSION(number_of_particles) ::  v_int        !<
181    REAL(wp), DIMENSION(number_of_particles) ::  w_int        !<
182    REAL(wp), DIMENSION(number_of_particles) ::  xv           !<
183    REAL(wp), DIMENSION(number_of_particles) ::  yv           !<
184    REAL(wp), DIMENSION(number_of_particles) ::  zv           !<
185
186    REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !<
187
188    CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' )
189
190!
191!-- Determine height of Prandtl layer and distance between Prandtl-layer
192!-- height and horizontal mean roughness height, which are required for
193!-- vertical logarithmic interpolation of horizontal particle speeds
194!-- (for particles below first vertical grid level).
195    z_p      = zu(nzb+1) - zw(nzb)
196    d_z_p_z0 = 1.0_wp / ( z_p - z0_av_global )
197
198    start_index = grid_particles(kp,jp,ip)%start_index
199    end_index   = grid_particles(kp,jp,ip)%end_index
200
201    xv = particles(1:number_of_particles)%x
202    yv = particles(1:number_of_particles)%y
203    zv = particles(1:number_of_particles)%z
204
205    DO  nb = 0, 7
206
207       i = ip
208       j = jp + block_offset(nb)%j_off
209       k = kp + block_offset(nb)%k_off
210
211!
212!--    Interpolate u velocity-component
213       DO  n = start_index(nb), end_index(nb)
214!
215!--       Interpolation of the u velocity component onto particle position. 
216!--       Particles are interpolation bi-linearly in the horizontal and a
217!--       linearly in the vertical. An exception is made for particles below
218!--       the first vertical grid level in case of a prandtl layer. In this
219!--       case the horizontal particle velocity components are determined using
220!--       Monin-Obukhov relations (if branch).
221!--       First, check if particle is located below first vertical grid level
222!--       (Prandtl-layer height)
223          IF ( prandtl_layer  .AND.  particles(n)%z < z_p )  THEN
224!
225!--          Resolved-scale horizontal particle velocity is zero below z0.
226             IF ( particles(n)%z < z0_av_global )  THEN
227                u_int(n) = 0.0_wp
228             ELSE
229!
230!--             Determine the sublayer. Further used as index.
231                height_p     = ( particles(n)%z - z0_av_global )            &
232                                     * REAL( number_of_sublayers, KIND=wp ) &
233                                     * d_z_p_z0
234!
235!--             Calculate LOG(z/z0) for exact particle height. Therefore,   
236!--             interpolate linearly between precalculated logarithm.
237                log_z_z0_int(n) = log_z_z0(INT(height_p))                      &
238                                 + ( height_p - INT(height_p) )                &
239                                 * ( log_z_z0(INT(height_p)+1)                 &
240                                      - log_z_z0(INT(height_p))                &
241                                   ) 
242!
243!--             Neutral solution is applied for all situations, e.g. also for
244!--             unstable and stable situations. Even though this is not exact
245!--             this saves a lot of CPU time since several calls of intrinsic
246!--             FORTRAN procedures (LOG, ATAN) are avoided, This is justified
247!--             as sensitivity studies revealed no significant effect of
248!--             using the neutral solution also for un/stable situations.
249!--             Calculated left and bottom index on u grid.
250                us_int = 0.5_wp * ( us(j,i) + us(j,i-1) ) 
251
252                u_int  = -usws(j,i) / ( us_int * kappa + 1E-10_wp )           & 
253                         * log_z_z0_int(n) - u_gtrans
254
255             ENDIF
256!
257!--       Particle above the first grid level. Bi-linear interpolation in the
258!--       horizontal and linear interpolation in the vertical direction.
259          ELSE
260
261             = xv(n) + ( 0.5_wp - i ) * dx
262             y  = yv(n) - j * dy
263             aa = x**2          + y**2
264             bb = ( dx - x )**2 + y**2
265             cc = x**2          + ( dy - y )**2
266             dd = ( dx - x )**2 + ( dy - y )**2
267             gg = aa + bb + cc + dd
268
269             u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
270                         + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) *            &
271                         u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
272
273             IF ( k == nzt )  THEN
274                u_int(n) = u_int_l
275             ELSE
276                u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1)  &
277                            + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) *           &
278                            u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
279                u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dz *                  &
280                           ( u_int_u - u_int_l )
281             ENDIF
282          ENDIF
283
284       ENDDO
285
286       i = ip + block_offset(nb)%i_off
287       j = jp
288       k = kp + block_offset(nb)%k_off
289!
290!--    Same procedure for interpolation of the v velocity-component
291       DO  n = start_index(nb), end_index(nb)
292
293          IF ( prandtl_layer  .AND.  particles(n)%z < z_p )  THEN
294
295             IF ( particles(n)%z < z0_av_global )  THEN
296!
297!--             Resolved-scale horizontal particle velocity is zero below z0.
298                v_int(n) = 0.0_wp
299             ELSE       
300!
301!--             Neutral solution is applied for all situations, e.g. also for
302!--             unstable and stable situations. Even though this is not exact
303!--             this saves a lot of CPU time since several calls of intrinsic
304!--             FORTRAN procedures (LOG, ATAN) are avoided, This is justified
305!--             as sensitivity studies revealed no significant effect of
306!--             using the neutral solution also for un/stable situations.
307!--             Calculated left and bottom index on v grid.
308                us_int = 0.5_wp * ( us(j,i) + us(j-1,i) )   
309
310                v_int  = -vsws(j,i) / ( us_int * kappa + 1E-10_wp )             &
311                         * log_z_z0_int(n) - v_gtrans
312             ENDIF
313          ELSE
314             = xv(n) - i * dx
315             y  = yv(n) + ( 0.5_wp - j ) * dy
316             aa = x**2          + y**2
317             bb = ( dx - x )**2 + y**2
318             cc = x**2          + ( dy - y )**2
319             dd = ( dx - x )**2 + ( dy - y )**2
320             gg = aa + bb + cc + dd
321
322             v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
323                       + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
324                       ) / ( 3.0_wp * gg ) - v_gtrans
325
326             IF ( k == nzt )  THEN
327                v_int(n) = v_int_l
328             ELSE
329                v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
330                          + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
331                          ) / ( 3.0_wp * gg ) - v_gtrans
332                v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dz * &
333                                  ( v_int_u - v_int_l )
334             ENDIF
335          ENDIF
336
337       ENDDO
338
339       i = ip + block_offset(nb)%i_off
340       j = jp + block_offset(nb)%j_off
341       k = kp-1
342!
343!--    Same procedure for interpolation of the w velocity-component
344       DO  n = start_index(nb), end_index(nb)
345
346          IF ( vertical_particle_advection(particles(n)%group) )  THEN
347
348             = xv(n) - i * dx
349             y  = yv(n) - j * dy
350             aa = x**2          + y**2
351             bb = ( dx - x )**2 + y**2
352             cc = x**2          + ( dy - y )**2
353             dd = ( dx - x )**2 + ( dy - y )**2
354             gg = aa + bb + cc + dd
355
356             w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
357                       + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
358                       ) / ( 3.0_wp * gg )
359
360             IF ( k == nzt )  THEN
361                w_int(n) = w_int_l
362             ELSE
363                w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + &
364                            ( gg-bb ) * w(k+1,j,i+1) + &
365                            ( gg-cc ) * w(k+1,j+1,i) + &
366                            ( gg-dd ) * w(k+1,j+1,i+1) &
367                          ) / ( 3.0_wp * gg )
368                w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dz * &
369                           ( w_int_u - w_int_l )
370             ENDIF
371
372          ELSE
373
374             w_int(n) = 0.0_wp
375
376          ENDIF
377
378       ENDDO
379
380    ENDDO
381
382!-- Interpolate and calculate quantities needed for calculating the SGS
383!-- velocities
384    IF ( use_sgs_for_particles )  THEN
385
386       IF ( topography == 'flat' )  THEN
387
388          DO  nb = 0,7
389
390             i = ip + block_offset(nb)%i_off
391             j = jp + block_offset(nb)%j_off
392             k = kp + block_offset(nb)%k_off
393
394             DO  n = start_index(nb), end_index(nb)
395!
396!--             Interpolate TKE
397                x  = xv(n) - i * dx
398                y  = yv(n) - j * dy
399                aa = x**2          + y**2
400                bb = ( dx - x )**2 + y**2
401                cc = x**2          + ( dy - y )**2
402                dd = ( dx - x )**2 + ( dy - y )**2
403                gg = aa + bb + cc + dd
404
405                e_int_l = ( ( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1)   &
406                          + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
407                          ) / ( 3.0_wp * gg )
408
409                IF ( k+1 == nzt+1 )  THEN
410                   e_int(n) = e_int_l
411                ELSE
412                   e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
413                               ( gg - bb ) * e(k+1,j,i+1) + &
414                               ( gg - cc ) * e(k+1,j+1,i) + &
415                               ( gg - dd ) * e(k+1,j+1,i+1) &
416                            ) / ( 3.0_wp * gg )
417                   e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
418                                     ( e_int_u - e_int_l )
419                ENDIF
420!
421!--             Needed to avoid NaN particle velocities (this might not be
422!--             required any more)
423                IF ( e_int(n) <= 0.0_wp )  THEN
424                   e_int(n) = 1.0E-20_wp
425                ENDIF
426!
427!--             Interpolate the TKE gradient along x (adopt incides i,j,k and
428!--             all position variables from above (TKE))
429                de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
430                                ( gg - bb ) * de_dx(k,j,i+1) + &
431                                ( gg - cc ) * de_dx(k,j+1,i) + &
432                                ( gg - dd ) * de_dx(k,j+1,i+1) &
433                               ) / ( 3.0_wp * gg )
434
435                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
436                   de_dx_int(n) = de_dx_int_l
437                ELSE
438                   de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
439                                   ( gg - bb ) * de_dx(k+1,j,i+1) + &
440                                   ( gg - cc ) * de_dx(k+1,j+1,i) + &
441                                   ( gg - dd ) * de_dx(k+1,j+1,i+1) &
442                                  ) / ( 3.0_wp * gg )
443                   de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / dz * &
444                                              ( de_dx_int_u - de_dx_int_l )
445                ENDIF
446!
447!--             Interpolate the TKE gradient along y
448                de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
449                                ( gg - bb ) * de_dy(k,j,i+1) + &
450                                ( gg - cc ) * de_dy(k,j+1,i) + &
451                                ( gg - dd ) * de_dy(k,j+1,i+1) &
452                               ) / ( 3.0_wp * gg )
453                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
454                   de_dy_int(n) = de_dy_int_l
455                ELSE
456                   de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
457                                  ( gg - bb ) * de_dy(k+1,j,i+1) + &
458                                  ( gg - cc ) * de_dy(k+1,j+1,i) + &
459                                  ( gg - dd ) * de_dy(k+1,j+1,i+1) &
460                                  ) / ( 3.0_wp * gg )
461                   de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * &
462                                              ( de_dy_int_u - de_dy_int_l )
463                ENDIF
464
465!
466!--             Interpolate the TKE gradient along z
467                IF ( zv(n) < 0.5_wp * dz )  THEN
468                   de_dz_int(n) = 0.0_wp
469                ELSE
470                   de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
471                                   ( gg - bb ) * de_dz(k,j,i+1) + &
472                                   ( gg - cc ) * de_dz(k,j+1,i) + &
473                                   ( gg - dd ) * de_dz(k,j+1,i+1) &
474                                  ) / ( 3.0_wp * gg )
475
476                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
477                      de_dz_int(n) = de_dz_int_l
478                   ELSE
479                      de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
480                                      ( gg - bb ) * de_dz(k+1,j,i+1) + &
481                                      ( gg - cc ) * de_dz(k+1,j+1,i) + &
482                                      ( gg - dd ) * de_dz(k+1,j+1,i+1) &
483                                     ) / ( 3.0_wp * gg )
484                      de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) / dz * &
485                                                 ( de_dz_int_u - de_dz_int_l )
486                   ENDIF
487                ENDIF
488
489!
490!--             Interpolate the dissipation of TKE
491                diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
492                               ( gg - bb ) * diss(k,j,i+1) + &
493                               ( gg - cc ) * diss(k,j+1,i) + &
494                               ( gg - dd ) * diss(k,j+1,i+1) &
495                              ) / ( 3.0_wp * gg )
496
497                IF ( k == nzt )  THEN
498                   diss_int(n) = diss_int_l
499                ELSE
500                   diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
501                                  ( gg - bb ) * diss(k+1,j,i+1) + &
502                                  ( gg - cc ) * diss(k+1,j+1,i) + &
503                                  ( gg - dd ) * diss(k+1,j+1,i+1) &
504                                 ) / ( 3.0_wp * gg )
505                   diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz * &
506                                           ( diss_int_u - diss_int_l )
507                ENDIF
508
509             ENDDO
510          ENDDO
511
512       ELSE ! non-flat topography, e.g., buildings
513
514          DO  n = 1, number_of_particles
515
516             i = particles(n)%x * ddx
517             j = particles(n)%y * ddy
518             k = ( zv(n) + 0.5_wp * dz * atmos_ocean_sign ) / dz  &
519                 + offset_ocean_nzt                      ! only exact if eq.dist
520!
521!--          In case that there are buildings it has to be determined
522!--          how many of the gridpoints defining the particle box are
523!--          situated within a building
524!--          gp_outside_of_building(1): i,j,k
525!--          gp_outside_of_building(2): i,j+1,k
526!--          gp_outside_of_building(3): i,j,k+1
527!--          gp_outside_of_building(4): i,j+1,k+1
528!--          gp_outside_of_building(5): i+1,j,k
529!--          gp_outside_of_building(6): i+1,j+1,k
530!--          gp_outside_of_building(7): i+1,j,k+1
531!--          gp_outside_of_building(8): i+1,j+1,k+1
532
533             gp_outside_of_building = 0
534             location = 0.0_wp
535             num_gp = 0
536
537             IF ( k > nzb_s_inner(j,i)  .OR.  nzb_s_inner(j,i) == 0 )  THEN
538                num_gp = num_gp + 1
539                gp_outside_of_building(1) = 1
540                location(num_gp,1) = i * dx
541                location(num_gp,2) = j * dy
542                location(num_gp,3) = k * dz - 0.5_wp * dz
543                ei(num_gp)     = e(k,j,i)
544                dissi(num_gp)  = diss(k,j,i)
545                de_dxi(num_gp) = de_dx(k,j,i)
546                de_dyi(num_gp) = de_dy(k,j,i)
547                de_dzi(num_gp) = de_dz(k,j,i)
548             ENDIF
549
550             IF ( k > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 ) &
551             THEN
552                num_gp = num_gp + 1
553                gp_outside_of_building(2) = 1
554                location(num_gp,1) = i * dx
555                location(num_gp,2) = (j+1) * dy
556                location(num_gp,3) = k * dz - 0.5_wp * dz
557                ei(num_gp)     = e(k,j+1,i)
558                dissi(num_gp)  = diss(k,j+1,i)
559                de_dxi(num_gp) = de_dx(k,j+1,i)
560                de_dyi(num_gp) = de_dy(k,j+1,i)
561                de_dzi(num_gp) = de_dz(k,j+1,i)
562             ENDIF
563
564             IF ( k+1 > nzb_s_inner(j,i)  .OR.  nzb_s_inner(j,i) == 0 )  THEN
565                num_gp = num_gp + 1
566                gp_outside_of_building(3) = 1
567                location(num_gp,1) = i * dx
568                location(num_gp,2) = j * dy
569                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
570                ei(num_gp)     = e(k+1,j,i)
571                dissi(num_gp)  = diss(k+1,j,i)
572                de_dxi(num_gp) = de_dx(k+1,j,i)
573                de_dyi(num_gp) = de_dy(k+1,j,i)
574                de_dzi(num_gp) = de_dz(k+1,j,i)
575             ENDIF
576
577             IF ( k+1 > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 ) &
578             THEN
579                num_gp = num_gp + 1
580                gp_outside_of_building(4) = 1
581                location(num_gp,1) = i * dx
582                location(num_gp,2) = (j+1) * dy
583                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
584                ei(num_gp)     = e(k+1,j+1,i)
585                dissi(num_gp)  = diss(k+1,j+1,i)
586                de_dxi(num_gp) = de_dx(k+1,j+1,i)
587                de_dyi(num_gp) = de_dy(k+1,j+1,i)
588                de_dzi(num_gp) = de_dz(k+1,j+1,i)
589             ENDIF
590
591             IF ( k > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 ) &
592             THEN
593                num_gp = num_gp + 1
594                gp_outside_of_building(5) = 1
595                location(num_gp,1) = (i+1) * dx
596                location(num_gp,2) = j * dy
597                location(num_gp,3) = k * dz - 0.5_wp * dz
598                ei(num_gp)     = e(k,j,i+1)
599                dissi(num_gp)  = diss(k,j,i+1)
600                de_dxi(num_gp) = de_dx(k,j,i+1)
601                de_dyi(num_gp) = de_dy(k,j,i+1)
602                de_dzi(num_gp) = de_dz(k,j,i+1)
603             ENDIF
604
605             IF ( k > nzb_s_inner(j+1,i+1)  .OR.  nzb_s_inner(j+1,i+1) == 0 ) &
606             THEN
607                num_gp = num_gp + 1
608                gp_outside_of_building(6) = 1
609                location(num_gp,1) = (i+1) * dx
610                location(num_gp,2) = (j+1) * dy
611                location(num_gp,3) = k * dz - 0.5_wp * dz
612                ei(num_gp)     = e(k,j+1,i+1)
613                dissi(num_gp)  = diss(k,j+1,i+1)
614                de_dxi(num_gp) = de_dx(k,j+1,i+1)
615                de_dyi(num_gp) = de_dy(k,j+1,i+1)
616                de_dzi(num_gp) = de_dz(k,j+1,i+1)
617             ENDIF
618
619             IF ( k+1 > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 ) &
620             THEN
621                num_gp = num_gp + 1
622                gp_outside_of_building(7) = 1
623                location(num_gp,1) = (i+1) * dx
624                location(num_gp,2) = j * dy
625                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
626                ei(num_gp)     = e(k+1,j,i+1)
627                dissi(num_gp)  = diss(k+1,j,i+1)
628                de_dxi(num_gp) = de_dx(k+1,j,i+1)
629                de_dyi(num_gp) = de_dy(k+1,j,i+1)
630                de_dzi(num_gp) = de_dz(k+1,j,i+1)
631             ENDIF
632
633             IF ( k+1 > nzb_s_inner(j+1,i+1)  .OR.  nzb_s_inner(j+1,i+1) == 0)&
634             THEN
635                num_gp = num_gp + 1
636                gp_outside_of_building(8) = 1
637                location(num_gp,1) = (i+1) * dx
638                location(num_gp,2) = (j+1) * dy
639                location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
640                ei(num_gp)     = e(k+1,j+1,i+1)
641                dissi(num_gp)  = diss(k+1,j+1,i+1)
642                de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
643                de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
644                de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
645             ENDIF
646
647!
648!--          If all gridpoints are situated outside of a building, then the
649!--          ordinary interpolation scheme can be used.
650             IF ( num_gp == 8 )  THEN
651
652                = particles(n)%x - i * dx
653                y  = particles(n)%y - j * dy
654                aa = x**2          + y**2
655                bb = ( dx - x )**2 + y**2
656                cc = x**2          + ( dy - y )**2
657                dd = ( dx - x )**2 + ( dy - y )**2
658                gg = aa + bb + cc + dd
659
660                e_int_l = ( ( gg - aa ) * e(k,j,i)   + ( gg - bb ) * e(k,j,i+1)   &
661                          + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) &
662                          ) / ( 3.0_wp * gg )
663
664                IF ( k == nzt )  THEN
665                   e_int(n) = e_int_l
666                ELSE
667                   e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
668                               ( gg - bb ) * e(k+1,j,i+1) + &
669                               ( gg - cc ) * e(k+1,j+1,i) + &
670                               ( gg - dd ) * e(k+1,j+1,i+1) &
671                             ) / ( 3.0_wp * gg )
672                   e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
673                                     ( e_int_u - e_int_l )
674                ENDIF
675!
676!--             Needed to avoid NaN particle velocities (this might not be
677!--             required any more)
678                IF ( e_int(n) <= 0.0_wp )  THEN
679                   e_int(n) = 1.0E-20_wp
680                ENDIF
681!
682!--             Interpolate the TKE gradient along x (adopt incides i,j,k
683!--             and all position variables from above (TKE))
684                de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
685                                ( gg - bb ) * de_dx(k,j,i+1) + &
686                                ( gg - cc ) * de_dx(k,j+1,i) + &
687                                ( gg - dd ) * de_dx(k,j+1,i+1) &
688                              ) / ( 3.0_wp * gg )
689
690                IF ( ( k == nzt )  .OR.  ( k == nzb ) )  THEN
691                   de_dx_int(n) = de_dx_int_l
692                ELSE
693                   de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
694                                   ( gg - bb ) * de_dx(k+1,j,i+1) + &
695                                   ( gg - cc ) * de_dx(k+1,j+1,i) + &
696                                   ( gg - dd ) * de_dx(k+1,j+1,i+1) &
697                                 ) / ( 3.0_wp * gg )
698                   de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / &
699                                           dz * ( de_dx_int_u - de_dx_int_l )
700                ENDIF
701
702!
703!--             Interpolate the TKE gradient along y
704                de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
705                                ( gg - bb ) * de_dy(k,j,i+1) + &
706                                ( gg - cc ) * de_dy(k,j+1,i) + &
707                                ( gg - dd ) * de_dy(k,j+1,i+1) &
708                              ) / ( 3.0_wp * gg )
709                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
710                   de_dy_int(n) = de_dy_int_l
711                ELSE
712                   de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
713                                   ( gg - bb ) * de_dy(k+1,j,i+1) + &
714                                   ( gg - cc ) * de_dy(k+1,j+1,i) + &
715                                   ( gg - dd ) * de_dy(k+1,j+1,i+1) &
716                                 ) / ( 3.0_wp * gg )
717                   de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / &
718                                           dz * ( de_dy_int_u - de_dy_int_l )
719                ENDIF
720
721!
722!--             Interpolate the TKE gradient along z
723                IF ( zv(n) < 0.5_wp * dz )  THEN
724                   de_dz_int(n) = 0.0_wp
725                ELSE
726                   de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
727                                   ( gg - bb ) * de_dz(k,j,i+1) + &
728                                   ( gg - cc ) * de_dz(k,j+1,i) + &
729                                   ( gg - dd ) * de_dz(k,j+1,i+1) &
730                                 ) / ( 3.0_wp * gg )
731
732                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
733                      de_dz_int(n) = de_dz_int_l
734                   ELSE
735                      de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
736                                      ( gg - bb ) * de_dz(k+1,j,i+1) + &
737                                      ( gg - cc ) * de_dz(k+1,j+1,i) + &
738                                      ( gg - dd ) * de_dz(k+1,j+1,i+1) &
739                                    ) / ( 3.0_wp * gg )
740                      de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) /&
741                                           dz * ( de_dz_int_u - de_dz_int_l )
742                   ENDIF
743                ENDIF
744
745!
746!--             Interpolate the dissipation of TKE
747                diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
748                               ( gg - bb ) * diss(k,j,i+1) + &
749                               ( gg - cc ) * diss(k,j+1,i) + &
750                               ( gg - dd ) * diss(k,j+1,i+1) &
751                             ) / ( 3.0_wp * gg )
752
753                IF ( k == nzt )  THEN
754                   diss_int(n) = diss_int_l
755                ELSE
756                   diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
757                                  ( gg - bb ) * diss(k+1,j,i+1) + &
758                                  ( gg - cc ) * diss(k+1,j+1,i) + &
759                                  ( gg - dd ) * diss(k+1,j+1,i+1) &
760                                ) / ( 3.0_wp * gg )
761                   diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *&
762                                           ( diss_int_u - diss_int_l )
763                ENDIF
764
765             ELSE
766
767!
768!--             If wall between gridpoint 1 and gridpoint 5, then
769!--             Neumann boundary condition has to be applied
770                IF ( gp_outside_of_building(1) == 1  .AND. &
771                     gp_outside_of_building(5) == 0 )  THEN
772                   num_gp = num_gp + 1
773                   location(num_gp,1) = i * dx + 0.5_wp * dx
774                   location(num_gp,2) = j * dy
775                   location(num_gp,3) = k * dz - 0.5_wp * dz
776                   ei(num_gp)     = e(k,j,i)
777                   dissi(num_gp)  = diss(k,j,i)
778                   de_dxi(num_gp) = 0.0_wp
779                   de_dyi(num_gp) = de_dy(k,j,i)
780                   de_dzi(num_gp) = de_dz(k,j,i)
781                ENDIF
782
783                IF ( gp_outside_of_building(5) == 1  .AND. &
784                   gp_outside_of_building(1) == 0 )  THEN
785                   num_gp = num_gp + 1
786                   location(num_gp,1) = i * dx + 0.5_wp * dx
787                   location(num_gp,2) = j * dy
788                   location(num_gp,3) = k * dz - 0.5_wp * dz
789                   ei(num_gp)     = e(k,j,i+1)
790                   dissi(num_gp)  = diss(k,j,i+1)
791                   de_dxi(num_gp) = 0.0_wp
792                   de_dyi(num_gp) = de_dy(k,j,i+1)
793                   de_dzi(num_gp) = de_dz(k,j,i+1)
794                ENDIF
795
796!
797!--             If wall between gridpoint 5 and gridpoint 6, then
798!--             then Neumann boundary condition has to be applied
799                IF ( gp_outside_of_building(5) == 1  .AND. &
800                     gp_outside_of_building(6) == 0 )  THEN
801                   num_gp = num_gp + 1
802                   location(num_gp,1) = (i+1) * dx
803                   location(num_gp,2) = j * dy + 0.5_wp * dy
804                   location(num_gp,3) = k * dz - 0.5_wp * dz
805                   ei(num_gp)     = e(k,j,i+1)
806                   dissi(num_gp)  = diss(k,j,i+1)
807                   de_dxi(num_gp) = de_dx(k,j,i+1)
808                   de_dyi(num_gp) = 0.0_wp
809                   de_dzi(num_gp) = de_dz(k,j,i+1)
810                ENDIF
811
812                IF ( gp_outside_of_building(6) == 1  .AND. &
813                     gp_outside_of_building(5) == 0 )  THEN
814                   num_gp = num_gp + 1
815                   location(num_gp,1) = (i+1) * dx
816                   location(num_gp,2) = j * dy + 0.5_wp * dy
817                   location(num_gp,3) = k * dz - 0.5_wp * dz
818                   ei(num_gp)     = e(k,j+1,i+1)
819                   dissi(num_gp)  = diss(k,j+1,i+1)
820                   de_dxi(num_gp) = de_dx(k,j+1,i+1)
821                   de_dyi(num_gp) = 0.0_wp
822                   de_dzi(num_gp) = de_dz(k,j+1,i+1)
823                ENDIF
824
825!
826!--             If wall between gridpoint 2 and gridpoint 6, then
827!--             Neumann boundary condition has to be applied
828                IF ( gp_outside_of_building(2) == 1  .AND. &
829                     gp_outside_of_building(6) == 0 )  THEN
830                   num_gp = num_gp + 1
831                   location(num_gp,1) = i * dx + 0.5_wp * dx
832                   location(num_gp,2) = (j+1) * dy
833                   location(num_gp,3) = k * dz - 0.5_wp * dz
834                   ei(num_gp)     = e(k,j+1,i)
835                   dissi(num_gp)  = diss(k,j+1,i)
836                   de_dxi(num_gp) = 0.0_wp
837                   de_dyi(num_gp) = de_dy(k,j+1,i)
838                   de_dzi(num_gp) = de_dz(k,j+1,i)
839                ENDIF
840
841                IF ( gp_outside_of_building(6) == 1  .AND. &
842                     gp_outside_of_building(2) == 0 )  THEN
843                   num_gp = num_gp + 1
844                   location(num_gp,1) = i * dx + 0.5_wp * dx
845                   location(num_gp,2) = (j+1) * dy
846                   location(num_gp,3) = k * dz - 0.5_wp * dz
847                   ei(num_gp)     = e(k,j+1,i+1)
848                   dissi(num_gp)  = diss(k,j+1,i+1)
849                   de_dxi(num_gp) = 0.0_wp
850                   de_dyi(num_gp) = de_dy(k,j+1,i+1)
851                   de_dzi(num_gp) = de_dz(k,j+1,i+1)
852                ENDIF
853
854!
855!--             If wall between gridpoint 1 and gridpoint 2, then
856!--             Neumann boundary condition has to be applied
857                IF ( gp_outside_of_building(1) == 1  .AND. &
858                     gp_outside_of_building(2) == 0 )  THEN
859                   num_gp = num_gp + 1
860                   location(num_gp,1) = i * dx
861                   location(num_gp,2) = j * dy + 0.5_wp * dy
862                   location(num_gp,3) = k * dz - 0.5_wp * dz
863                   ei(num_gp)     = e(k,j,i)
864                   dissi(num_gp)  = diss(k,j,i)
865                   de_dxi(num_gp) = de_dx(k,j,i)
866                   de_dyi(num_gp) = 0.0_wp
867                   de_dzi(num_gp) = de_dz(k,j,i)
868                ENDIF
869
870                IF ( gp_outside_of_building(2) == 1  .AND. &
871                     gp_outside_of_building(1) == 0 )  THEN
872                   num_gp = num_gp + 1
873                   location(num_gp,1) = i * dx
874                   location(num_gp,2) = j * dy + 0.5_wp * dy
875                   location(num_gp,3) = k * dz - 0.5_wp * dz
876                   ei(num_gp)     = e(k,j+1,i)
877                   dissi(num_gp)  = diss(k,j+1,i)
878                   de_dxi(num_gp) = de_dx(k,j+1,i)
879                   de_dyi(num_gp) = 0.0_wp
880                   de_dzi(num_gp) = de_dz(k,j+1,i)
881                ENDIF
882
883!
884!--             If wall between gridpoint 3 and gridpoint 7, then
885!--             Neumann boundary condition has to be applied
886                IF ( gp_outside_of_building(3) == 1  .AND. &
887                     gp_outside_of_building(7) == 0 )  THEN
888                   num_gp = num_gp + 1
889                   location(num_gp,1) = i * dx + 0.5_wp * dx
890                   location(num_gp,2) = j * dy
891                   location(num_gp,3) = k * dz + 0.5_wp * dz
892                   ei(num_gp)     = e(k+1,j,i)
893                   dissi(num_gp)  = diss(k+1,j,i)
894                   de_dxi(num_gp) = 0.0_wp
895                   de_dyi(num_gp) = de_dy(k+1,j,i)
896                   de_dzi(num_gp) = de_dz(k+1,j,i) 
897                ENDIF
898
899                IF ( gp_outside_of_building(7) == 1  .AND. &
900                     gp_outside_of_building(3) == 0 )  THEN
901                   num_gp = num_gp + 1
902                   location(num_gp,1) = i * dx + 0.5_wp * dx
903                   location(num_gp,2) = j * dy
904                   location(num_gp,3) = k * dz + 0.5_wp * dz
905                   ei(num_gp)     = e(k+1,j,i+1)
906                   dissi(num_gp)  = diss(k+1,j,i+1)
907                   de_dxi(num_gp) = 0.0_wp
908                   de_dyi(num_gp) = de_dy(k+1,j,i+1)
909                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
910                ENDIF
911
912!
913!--             If wall between gridpoint 7 and gridpoint 8, then
914!--             Neumann boundary condition has to be applied
915                IF ( gp_outside_of_building(7) == 1  .AND. &
916                     gp_outside_of_building(8) == 0 )  THEN
917                   num_gp = num_gp + 1
918                   location(num_gp,1) = (i+1) * dx
919                   location(num_gp,2) = j * dy + 0.5_wp * dy
920                   location(num_gp,3) = k * dz + 0.5_wp * dz
921                   ei(num_gp)     = e(k+1,j,i+1)
922                   dissi(num_gp)  = diss(k+1,j,i+1)
923                   de_dxi(num_gp) = de_dx(k+1,j,i+1)
924                   de_dyi(num_gp) = 0.0_wp
925                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
926                ENDIF
927
928                IF ( gp_outside_of_building(8) == 1  .AND. &
929                     gp_outside_of_building(7) == 0 )  THEN
930                   num_gp = num_gp + 1
931                   location(num_gp,1) = (i+1) * dx
932                   location(num_gp,2) = j * dy + 0.5_wp * dy
933                   location(num_gp,3) = k * dz + 0.5_wp * dz
934                   ei(num_gp)     = e(k+1,j+1,i+1)
935                   dissi(num_gp)  = diss(k+1,j+1,i+1)
936                   de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
937                   de_dyi(num_gp) = 0.0_wp
938                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
939                ENDIF
940
941!
942!--             If wall between gridpoint 4 and gridpoint 8, then
943!--             Neumann boundary condition has to be applied
944                IF ( gp_outside_of_building(4) == 1  .AND. &
945                     gp_outside_of_building(8) == 0 )  THEN
946                   num_gp = num_gp + 1
947                   location(num_gp,1) = i * dx + 0.5_wp * dx
948                   location(num_gp,2) = (j+1) * dy
949                   location(num_gp,3) = k * dz + 0.5_wp * dz
950                   ei(num_gp)     = e(k+1,j+1,i)
951                   dissi(num_gp)  = diss(k+1,j+1,i)
952                   de_dxi(num_gp) = 0.0_wp
953                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
954                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
955                ENDIF
956
957                IF ( gp_outside_of_building(8) == 1  .AND. &
958                     gp_outside_of_building(4) == 0 )  THEN
959                   num_gp = num_gp + 1
960                   location(num_gp,1) = i * dx + 0.5_wp * dx
961                   location(num_gp,2) = (j+1) * dy
962                   location(num_gp,3) = k * dz + 0.5_wp * dz
963                   ei(num_gp)     = e(k+1,j+1,i+1)
964                   dissi(num_gp)  = diss(k+1,j+1,i+1)
965                   de_dxi(num_gp) = 0.0_wp
966                   de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
967                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
968                ENDIF
969
970!
971!--             If wall between gridpoint 3 and gridpoint 4, then
972!--             Neumann boundary condition has to be applied
973                IF ( gp_outside_of_building(3) == 1  .AND. &
974                     gp_outside_of_building(4) == 0 )  THEN
975                   num_gp = num_gp + 1
976                   location(num_gp,1) = i * dx
977                   location(num_gp,2) = j * dy + 0.5_wp * dy
978                   location(num_gp,3) = k * dz + 0.5_wp * dz
979                   ei(num_gp)     = e(k+1,j,i)
980                   dissi(num_gp)  = diss(k+1,j,i)
981                   de_dxi(num_gp) = de_dx(k+1,j,i)
982                   de_dyi(num_gp) = 0.0_wp
983                   de_dzi(num_gp) = de_dz(k+1,j,i)
984                ENDIF
985
986                IF ( gp_outside_of_building(4) == 1  .AND. &
987                     gp_outside_of_building(3) == 0 )  THEN
988                   num_gp = num_gp + 1
989                   location(num_gp,1) = i * dx
990                   location(num_gp,2) = j * dy + 0.5_wp * dy
991                   location(num_gp,3) = k * dz + 0.5_wp * dz
992                   ei(num_gp)     = e(k+1,j+1,i)
993                   dissi(num_gp)  = diss(k+1,j+1,i)
994                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
995                   de_dyi(num_gp) = 0.0_wp
996                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
997                ENDIF
998
999!
1000!--             If wall between gridpoint 1 and gridpoint 3, then
1001!--             Neumann boundary condition has to be applied
1002!--             (only one case as only building beneath is possible)
1003                IF ( gp_outside_of_building(1) == 0  .AND. &
1004                     gp_outside_of_building(3) == 1 )  THEN
1005                   num_gp = num_gp + 1
1006                   location(num_gp,1) = i * dx
1007                   location(num_gp,2) = j * dy
1008                   location(num_gp,3) = k * dz
1009                   ei(num_gp)     = e(k+1,j,i)
1010                   dissi(num_gp)  = diss(k+1,j,i)
1011                   de_dxi(num_gp) = de_dx(k+1,j,i)
1012                   de_dyi(num_gp) = de_dy(k+1,j,i)
1013                   de_dzi(num_gp) = 0.0_wp
1014                ENDIF
1015
1016!
1017!--             If wall between gridpoint 5 and gridpoint 7, then
1018!--             Neumann boundary condition has to be applied
1019!--             (only one case as only building beneath is possible)
1020                IF ( gp_outside_of_building(5) == 0  .AND. &
1021                     gp_outside_of_building(7) == 1 )  THEN
1022                   num_gp = num_gp + 1
1023                   location(num_gp,1) = (i+1) * dx
1024                   location(num_gp,2) = j * dy
1025                   location(num_gp,3) = k * dz
1026                   ei(num_gp)     = e(k+1,j,i+1)
1027                   dissi(num_gp)  = diss(k+1,j,i+1)
1028                   de_dxi(num_gp) = de_dx(k+1,j,i+1)
1029                   de_dyi(num_gp) = de_dy(k+1,j,i+1)
1030                   de_dzi(num_gp) = 0.0_wp
1031                ENDIF
1032
1033!
1034!--             If wall between gridpoint 2 and gridpoint 4, then
1035!--             Neumann boundary condition has to be applied
1036!--             (only one case as only building beneath is possible)
1037                IF ( gp_outside_of_building(2) == 0  .AND. &
1038                     gp_outside_of_building(4) == 1 )  THEN
1039                   num_gp = num_gp + 1
1040                   location(num_gp,1) = i * dx
1041                   location(num_gp,2) = (j+1) * dy
1042                   location(num_gp,3) = k * dz
1043                   ei(num_gp)     = e(k+1,j+1,i)
1044                   dissi(num_gp)  = diss(k+1,j+1,i)
1045                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
1046                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
1047                   de_dzi(num_gp) = 0.0_wp
1048                ENDIF
1049
1050!
1051!--             If wall between gridpoint 6 and gridpoint 8, then
1052!--             Neumann boundary condition has to be applied
1053!--             (only one case as only building beneath is possible)
1054                IF ( gp_outside_of_building(6) == 0  .AND. &
1055                     gp_outside_of_building(8) == 1 )  THEN
1056                   num_gp = num_gp + 1
1057                   location(num_gp,1) = (i+1) * dx
1058                   location(num_gp,2) = (j+1) * dy
1059                   location(num_gp,3) = k * dz
1060                   ei(num_gp)     = e(k+1,j+1,i+1)
1061                   dissi(num_gp)  = diss(k+1,j+1,i+1)
1062                   de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
1063                   de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
1064                   de_dzi(num_gp) = 0.0_wp
1065                ENDIF
1066
1067!
1068!--             Carry out the interpolation
1069                IF ( num_gp == 1 )  THEN
1070!
1071!--                If only one of the gridpoints is situated outside of the
1072!--                building, it follows that the values at the particle
1073!--                location are the same as the gridpoint values
1074                   e_int(n)     = ei(num_gp)
1075                   diss_int(n)  = dissi(num_gp)
1076                   de_dx_int(n) = de_dxi(num_gp)
1077                   de_dy_int(n) = de_dyi(num_gp)
1078                   de_dz_int(n) = de_dzi(num_gp)
1079                ELSE IF ( num_gp > 1 )  THEN
1080
1081                   d_sum = 0.0_wp
1082!
1083!--                Evaluation of the distances between the gridpoints
1084!--                contributing to the interpolated values, and the particle
1085!--                location
1086                   DO  agp = 1, num_gp
1087                      d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2  &
1088                                   + ( particles(n)%y-location(agp,2) )**2  &
1089                                   + ( zv(n)-location(agp,3) )**2
1090                      d_sum        = d_sum + d_gp_pl(agp)
1091                   ENDDO
1092
1093!
1094!--                Finally the interpolation can be carried out
1095                   e_int(n)     = 0.0_wp 
1096                   diss_int(n)  = 0.0_wp 
1097                   de_dx_int(n) = 0.0_wp 
1098                   de_dy_int(n) = 0.0_wp 
1099                   de_dz_int(n) = 0.0_wp 
1100                   DO  agp = 1, num_gp
1101                      e_int(n)     = e_int(n)     + ( d_sum - d_gp_pl(agp) ) * &
1102                                             ei(agp) / ( (num_gp-1) * d_sum )
1103                      diss_int(n)  = diss_int(n)  + ( d_sum - d_gp_pl(agp) ) * &
1104                                          dissi(agp) / ( (num_gp-1) * d_sum )
1105                      de_dx_int(n) = de_dx_int(n) + ( d_sum - d_gp_pl(agp) ) * &
1106                                         de_dxi(agp) / ( (num_gp-1) * d_sum )
1107                      de_dy_int(n) = de_dy_int(n) + ( d_sum - d_gp_pl(agp) ) * &
1108                                         de_dyi(agp) / ( (num_gp-1) * d_sum )
1109                      de_dz_int(n) = de_dz_int(n) + ( d_sum - d_gp_pl(agp) ) * &
1110                                         de_dzi(agp) / ( (num_gp-1) * d_sum )
1111                   ENDDO
1112
1113                ENDIF
1114
1115             ENDIF
1116          ENDDO
1117       ENDIF
1118
1119       DO nb = 0,7
1120          i = ip + block_offset(nb)%i_off
1121          j = jp + block_offset(nb)%j_off
1122          k = kp + block_offset(nb)%k_off
1123
1124          DO  n = start_index(nb), end_index(nb)
1125!
1126!--          Vertical interpolation of the horizontally averaged SGS TKE and
1127!--          resolved-scale velocity variances and use the interpolated values
1128!--          to calculate the coefficient fs, which is a measure of the ratio
1129!--          of the subgrid-scale turbulent kinetic energy to the total amount
1130!--          of turbulent kinetic energy.
1131             IF ( k == 0 )  THEN
1132                e_mean_int = hom(0,1,8,0)
1133             ELSE
1134                e_mean_int = hom(k,1,8,0) +                                    &
1135                                           ( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
1136                                           ( zu(k+1) - zu(k) ) *               &
1137                                           ( zv(n) - zu(k) )
1138             ENDIF
1139
1140             kw = kp - 1
1141
1142             IF ( k == 0 )  THEN
1143                aa  = hom(k+1,1,30,0)  * ( zv(n) / &
1144                                         ( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
1145                bb  = hom(k+1,1,31,0)  * ( zv(n) / &
1146                                         ( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
1147                cc  = hom(kw+1,1,32,0) * ( zv(n) / &
1148                                         ( 1.0_wp * ( zw(kw+1) - zw(kw) ) ) )
1149             ELSE
1150                aa  = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) *    &
1151                           ( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
1152                bb  = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) *    &
1153                           ( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
1154                cc  = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) *   &
1155                           ( ( zv(n) - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
1156             ENDIF
1157
1158             vv_int = ( 1.0_wp / 3.0_wp ) * ( aa + bb + cc )
1159!
1160!--          Needed to avoid NaN particle velocities. The value of 1.0 is just
1161!--          an educated guess for the given case.
1162             IF ( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int == 0.0_wp )  THEN
1163                fs_int(n) = 1.0_wp
1164             ELSE
1165                fs_int(n) = ( 2.0_wp / 3.0_wp ) * e_mean_int /                 &
1166                            ( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int )
1167             ENDIF
1168
1169          ENDDO
1170       ENDDO
1171
1172       DO  n = 1, number_of_particles
1173
1174         rg(n,1) = random_gauss( iran_part, 5.0_wp )
1175         rg(n,2) = random_gauss( iran_part, 5.0_wp )
1176         rg(n,3) = random_gauss( iran_part, 5.0_wp )
1177
1178       ENDDO
1179
1180       DO  n = 1, number_of_particles
1181!
1182!--       Calculate the Lagrangian timescale according to Weil et al. (2004).
1183          lagr_timescale = ( 4.0_wp * e_int(n) ) / &
1184                           ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) )
1185
1186!
1187!--       Calculate the next particle timestep. dt_gap is the time needed to
1188!--       complete the current LES timestep.
1189          dt_gap = dt_3d - particles(n)%dt_sum
1190          dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale, dt_gap )
1191
1192!
1193!--       The particle timestep should not be too small in order to prevent
1194!--       the number of particle timesteps of getting too large
1195          IF ( dt_particle(n) < dt_min_part  .AND.  dt_min_part < dt_gap )  THEN
1196             dt_particle(n) = dt_min_part
1197          ENDIF
1198
1199!
1200!--       Calculate the SGS velocity components
1201          IF ( particles(n)%age == 0.0_wp )  THEN
1202!
1203!--          For new particles the SGS components are derived from the SGS
1204!--          TKE. Limit the Gaussian random number to the interval
1205!--          [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities
1206!--          from becoming unrealistically large.
1207             particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wfu_part * e_int(n) ) *   &
1208                                  ( rg(n,1) - 1.0_wp )
1209             particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wfv_part * e_int(n) ) *   &
1210                                  ( rg(n,2) - 1.0_wp )
1211             particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wfw_part * e_int(n) ) *   &
1212                                  ( rg(n,3) - 1.0_wp )
1213
1214          ELSE
1215!
1216!--          Restriction of the size of the new timestep: compared to the
1217!--          previous timestep the increase must not exceed 200%
1218
1219             dt_particle_m = particles(n)%age - particles(n)%age_m
1220             IF ( dt_particle(n) > 2.0_wp * dt_particle_m )  THEN
1221                dt_particle(n) = 2.0_wp * dt_particle_m
1222             ENDIF
1223
1224!
1225!--          For old particles the SGS components are correlated with the
1226!--          values from the previous timestep. Random numbers have also to
1227!--          be limited (see above).
1228!--          As negative values for the subgrid TKE are not allowed, the
1229!--          change of the subgrid TKE with time cannot be smaller than
1230!--          -e_int(n)/dt_particle. This value is used as a lower boundary
1231!--          value for the change of TKE
1232
1233             de_dt_min = - e_int(n) / dt_particle(n)
1234
1235             de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m
1236
1237             IF ( de_dt < de_dt_min )  THEN
1238                de_dt = de_dt_min
1239             ENDIF
1240
1241             particles(n)%rvar1 = particles(n)%rvar1 - fs_int(n) * c_0 *       &
1242                           diss_int(n) * particles(n)%rvar1 * dt_particle(n) / &
1243                           ( 4.0_wp * sgs_wfu_part * e_int(n) ) +              &
1244                           ( 2.0_wp * sgs_wfu_part * de_dt *                   &
1245                             particles(n)%rvar1 /                              &
1246                             ( 2.0_wp * sgs_wfu_part * e_int(n) ) +            &
1247                             de_dx_int(n)                                      &
1248                           ) * dt_particle(n) / 2.0_wp          +              &
1249                           SQRT( fs_int(n) * c_0 * diss_int(n) ) *             &
1250                           ( rg(n,1) - 1.0_wp ) *                              &
1251                           SQRT( dt_particle(n) )
1252
1253             particles(n)%rvar2 = particles(n)%rvar2 - fs_int(n) * c_0 *       &
1254                           diss_int(n) * particles(n)%rvar2 * dt_particle(n) / &
1255                           ( 4.0_wp * sgs_wfv_part * e_int(n) ) +              &
1256                           ( 2.0_wp * sgs_wfv_part * de_dt *                   &
1257                             particles(n)%rvar2 /                              &
1258                             ( 2.0_wp * sgs_wfv_part * e_int(n) ) +            &
1259                             de_dy_int(n)                                      &
1260                           ) * dt_particle(n) / 2.0_wp          +              &
1261                           SQRT( fs_int(n) * c_0 * diss_int(n) ) *             &
1262                           ( rg(n,2) - 1.0_wp ) *                              &
1263                           SQRT( dt_particle(n) )
1264
1265             particles(n)%rvar3 = particles(n)%rvar3 - fs_int(n) * c_0 *       &
1266                           diss_int(n) * particles(n)%rvar3 * dt_particle(n) / &
1267                           ( 4.0_wp * sgs_wfw_part * e_int(n) ) +              &
1268                           ( 2.0_wp * sgs_wfw_part * de_dt *                   &
1269                             particles(n)%rvar3 /                              &
1270                             ( 2.0_wp * sgs_wfw_part * e_int(n) ) +            &
1271                             de_dz_int(n)                                      &
1272                           ) * dt_particle(n) / 2.0_wp          +              &
1273                           SQRT( fs_int(n) * c_0 * diss_int(n) ) *             &
1274                           ( rg(n,3) - 1.0_wp ) *                              &
1275                           SQRT( dt_particle(n) )
1276
1277          ENDIF
1278          u_int(n) = u_int(n) + particles(n)%rvar1
1279          v_int(n) = v_int(n) + particles(n)%rvar2
1280          w_int(n) = w_int(n) + particles(n)%rvar3
1281
1282!
1283!--       Store the SGS TKE of the current timelevel which is needed for
1284!--       for calculating the SGS particle velocities at the next timestep
1285          particles(n)%e_m = e_int(n)
1286       ENDDO
1287
1288    ELSE
1289!
1290!--    If no SGS velocities are used, only the particle timestep has to
1291!--    be set
1292       dt_particle = dt_3d
1293
1294    ENDIF
1295!
1296!-- Store the old age of the particle ( needed to prevent that a
1297!-- particle crosses several PEs during one timestep, and for the
1298!-- evaluation of the subgrid particle velocity fluctuations )
1299    particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age
1300
1301    dens_ratio = particle_groups(particles(1:number_of_particles)%group)%density_ratio
1302
1303    IF ( ANY( dens_ratio == 0.0_wp ) )  THEN
1304       DO  n = 1, number_of_particles
1305
1306!
1307!--       Particle advection
1308          IF ( dens_ratio(n) == 0.0_wp )  THEN
1309!
1310!--          Pure passive transport (without particle inertia)
1311             particles(n)%x = xv(n) + u_int(n) * dt_particle(n)
1312             particles(n)%y = yv(n) + v_int(n) * dt_particle(n)
1313             particles(n)%z = zv(n) + w_int(n) * dt_particle(n)
1314
1315             particles(n)%speed_x = u_int(n)
1316             particles(n)%speed_y = v_int(n)
1317             particles(n)%speed_z = w_int(n)
1318
1319          ELSE
1320!
1321!--          Transport of particles with inertia
1322             particles(n)%x = particles(n)%x + particles(n)%speed_x * &
1323                                               dt_particle(n)
1324             particles(n)%y = particles(n)%y + particles(n)%speed_y * &
1325                                               dt_particle(n)
1326             particles(n)%z = particles(n)%z + particles(n)%speed_z * &
1327                                               dt_particle(n)
1328
1329!
1330!--          Update of the particle velocity
1331             IF ( cloud_droplets )  THEN
1332                exp_arg  = 4.5_wp * dens_ratio(n) * molecular_viscosity /      &
1333                           ( particles(n)%radius )**2 *                        &
1334                           ( 1.0_wp + 0.15_wp * ( 2.0_wp * particles(n)%radius &
1335                             * SQRT( ( u_int(n) - particles(n)%speed_x )**2 +  &
1336                                     ( v_int(n) - particles(n)%speed_y )**2 +  &
1337                                     ( w_int(n) - particles(n)%speed_z )**2 )  &
1338                                          / molecular_viscosity )**0.687_wp    &
1339                           )
1340
1341                exp_term = EXP( -exp_arg * dt_particle(n) )
1342             ELSEIF ( use_sgs_for_particles )  THEN
1343                exp_arg  = particle_groups(particles(n)%group)%exp_arg
1344                exp_term = EXP( -exp_arg * dt_particle(n) )
1345             ELSE
1346                exp_arg  = particle_groups(particles(n)%group)%exp_arg
1347                exp_term = particle_groups(particles(n)%group)%exp_term
1348             ENDIF
1349             particles(n)%speed_x = particles(n)%speed_x * exp_term +          &
1350                                    u_int(n) * ( 1.0_wp - exp_term )
1351             particles(n)%speed_y = particles(n)%speed_y * exp_term +          &
1352                                    v_int(n) * ( 1.0_wp - exp_term )
1353             particles(n)%speed_z = particles(n)%speed_z * exp_term +          &
1354                                    ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) *  &
1355                                    g / exp_arg ) * ( 1.0_wp - exp_term )
1356          ENDIF
1357
1358       ENDDO
1359   
1360    ELSE
1361
1362       DO  n = 1, number_of_particles
1363
1364!--       Transport of particles with inertia
1365          particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)
1366          particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)
1367          particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)
1368!
1369!--       Update of the particle velocity
1370          IF ( cloud_droplets )  THEN
1371
1372             exp_arg  = 4.5_wp * dens_ratio(n) * molecular_viscosity /         &
1373                        ( particles(n)%radius )**2 *                           &
1374                        ( 1.0_wp + 0.15_wp * ( 2.0_wp * particles(n)%radius *  &
1375                          SQRT( ( u_int(n) - particles(n)%speed_x )**2 +       &
1376                                ( v_int(n) - particles(n)%speed_y )**2 +       &
1377                                ( w_int(n) - particles(n)%speed_z )**2 ) /     &
1378                                         molecular_viscosity )**0.687_wp       &
1379                        )
1380
1381             exp_term = EXP( -exp_arg * dt_particle(n) )
1382          ELSEIF ( use_sgs_for_particles )  THEN
1383             exp_arg  = particle_groups(particles(n)%group)%exp_arg
1384             exp_term = EXP( -exp_arg * dt_particle(n) )
1385          ELSE
1386             exp_arg  = particle_groups(particles(n)%group)%exp_arg
1387             exp_term = particle_groups(particles(n)%group)%exp_term
1388          ENDIF
1389          particles(n)%speed_x = particles(n)%speed_x * exp_term +             &
1390                                 u_int(n) * ( 1.0_wp - exp_term )
1391          particles(n)%speed_y = particles(n)%speed_y * exp_term +             &
1392                                 v_int(n) * ( 1.0_wp - exp_term )
1393          particles(n)%speed_z = particles(n)%speed_z * exp_term +             &
1394                                 ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &
1395                                 exp_arg ) * ( 1.0_wp - exp_term )
1396       ENDDO
1397
1398    ENDIF
1399
1400    DO  n = 1, number_of_particles
1401!
1402!--    Increment the particle age and the total time that the particle
1403!--    has advanced within the particle timestep procedure
1404       particles(n)%age    = particles(n)%age    + dt_particle(n)
1405       particles(n)%dt_sum = particles(n)%dt_sum + dt_particle(n)
1406
1407!
1408!--    Check whether there is still a particle that has not yet completed
1409!--    the total LES timestep
1410       IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp )  THEN
1411          dt_3d_reached_l = .FALSE.
1412       ENDIF
1413
1414    ENDDO
1415
1416    CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
1417
1418 END SUBROUTINE lpm_advec
Note: See TracBrowser for help on using the repository browser.