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

Last change on this file since 1818 was 1818, checked in by maronga, 6 years ago

last commit documented / copyright update

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