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

Last change on this file since 1574 was 1370, checked in by raasch, 11 years ago

last commit documented

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