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

Last change on this file since 1360 was 1360, checked in by hoffmann, 8 years ago

last commit documented

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