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

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

last commit documented

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