Changeset 849 for palm/trunk/SOURCE/lpm.f90
- Timestamp:
- Mar 15, 2012 10:35:09 AM (12 years ago)
- File:
-
- 1 moved
-
palm/trunk/SOURCE/lpm.f90 (moved) (moved from palm/trunk/SOURCE/advec_particles.f90) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm.f90
r848 r849 1 SUBROUTINE advec_particles1 SUBROUTINE lpm 2 2 3 3 !------------------------------------------------------------------------------! 4 4 ! Current revisions: 5 5 ! ------------------ 6 ! 6 ! Original routine advec_particles split into several subroutines and renamed 7 ! lpm 7 8 ! 8 9 ! Former revisions: 9 10 ! ----------------- 10 11 ! $Id$ 11 !12 ! 835 2012-02-22 11:21:19Z raasch13 ! Bugfix: array diss can be used only in case of Wang kernel14 12 ! 15 13 ! 831 2012-02-22 00:29:39Z raasch … … 43 41 ! Splitting of parallel I/O (routine write_particles) 44 42 ! 45 ! 667 2010-12-23 12:06:00Z suehring/gryschka46 ! Declaration of de_dx, de_dy, de_dz adapted to additional47 ! ghost points. Furthermore the calls of exchange_horiz were modified.48 !49 ! 622 2010-12-10 08:08:13Z raasch50 ! optional barriers included in order to speed up collective operations51 !52 ! 519 2010-03-19 05:30:02Z raasch53 ! NetCDF4 output format allows size of particle array to be extended54 !55 ! 420 2010-01-13 15:10:53Z franke56 ! Own weighting factor for every cloud droplet is implemented and condensation57 ! and collision of cloud droplets are adjusted accordingly. +delta_v, -s_r3,58 ! -s_r459 ! Initialization of variables for the (sub-) timestep is moved to the beginning60 ! in order to enable deletion of cloud droplets due to collision processes.61 ! Collision efficiency for large cloud droplets has changed according to table62 ! of Rogers and Yau. (collision_efficiency)63 ! Bugfix: calculation of cloud droplet velocity64 ! Bugfix: transfer of particles at south/left edge when new position65 ! y>=(ny+0.5)*dy-1.e-12 or x>=(nx+0.5)*dx-1.e-12, very rare66 ! Bugfix: calculation of y (collision_efficiency)67 !68 ! 336 2009-06-10 11:19:35Z raasch69 ! Particle attributes are set with new routine set_particle_attributes.70 ! Vertical particle advection depends on the particle group.71 ! Output of NetCDF messages with aid of message handling routine.72 ! Output of messages replaced by message handling routine73 ! Bugfix: error in check, if particles moved further than one subdomain length.74 ! This check must not be applied for newly released particles75 ! Bugfix: several tail counters are initialized, particle_tail_coordinates is76 ! only written to file if its third index is > 077 !78 ! 212 2008-11-11 09:09:24Z raasch79 ! Bugfix in calculating k index in case of oceans runs (sort_particles)80 !81 ! 150 2008-02-29 08:19:58Z raasch82 ! Bottom boundary condition and vertical index calculations adjusted for83 ! ocean runs.84 !85 ! 119 2007-10-17 10:27:13Z raasch86 ! Sorting of particles is controlled by dt_sort_particles and moved from87 ! the SGS timestep loop after the end of this loop.88 ! Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail89 ! numbers along y90 ! Small bugfixes in the SGS part91 !92 ! 106 2007-08-16 14:30:26Z raasch93 ! remaining variables iran changed to iran_part94 !95 ! 95 2007-06-02 16:48:38Z raasch96 ! hydro_press renamed hyp97 !98 ! 75 2007-03-22 09:54:05Z raasch99 ! Particle reflection at vertical walls implemented in new subroutine100 ! particle_boundary_conds,101 ! vertical walls are regarded in the SGS model,102 ! + user_advec_particles, particles-package is now part of the defaut code,103 ! array arguments in sendrecv calls have to refer to first element (1) due to104 ! mpich (mpiI) interface requirements,105 ! 2nd+3rd argument removed from exchange horiz106 !107 ! 16 2007-02-15 13:16:47Z raasch108 ! Bugfix: wrong if-clause from revision 1.32109 !110 ! r4 | raasch | 2007-02-13 12:33:16 +0100 (Tue, 13 Feb 2007)111 ! RCS Log replace by Id keyword, revision history cleaned up112 !113 ! Revision 1.32 2007/02/11 12:48:20 raasch114 ! Allways the lower level k is used for interpolation115 ! Bugfix: new particles are released only if end_time_prel > simulated_time116 ! Bugfix: transfer of particles when x < -0.5*dx (0.0 before), etc.,117 ! index i,j used instead of cartesian (x,y) coordinate to check for118 ! transfer because this failed under very rare conditions119 ! Bugfix: calculation of number of particles with same radius as the current120 ! particle (cloud droplet code)121 !122 ! Revision 1.31 2006/08/17 09:21:01 raasch123 ! Two more compilation errors removed from the last revision124 !125 ! Revision 1.30 2006/08/17 09:11:17 raasch126 ! Two compilation errors removed from the last revision127 !128 ! Revision 1.29 2006/08/04 14:05:01 raasch129 ! Subgrid scale velocities are (optionally) included for calculating the130 ! particle advection, new counters trlp_count_sum, etc. for accumulating131 ! the number of particles exchanged between the subdomains during all132 ! sub-timesteps (if sgs velocities are included), +3d-arrays de_dx/y/z,133 ! izuf renamed iran, output of particle time series134 !135 43 ! Revision 1.1 1999/11/25 16:16:06 raasch 136 44 ! Initial revision … … 143 51 144 52 USE arrays_3d 145 USE cloud_parameters146 USE constants147 53 USE control_parameters 148 54 USE cpulog 149 USE grid_variables150 USE indices151 55 USE interfaces 152 USE lpm_collision_kernels_mod153 USE netcdf_control154 56 USE particle_attributes 155 57 USE pegrid 156 USE random_function_mod157 58 USE statistics 158 59 159 60 IMPLICIT NONE 160 61 161 INTEGER :: agp, deleted_particles, deleted_tails, eclass, i, ie, ii, inc, & 162 internal_timestep_count, is, j, jj, js, jtry, k, kk, kw, m, n, & 163 nc, nd, nn, num_gp, pse, psi, rclass_l, rclass_s, & 164 tlength, trlp_count, trlp_count_sum, trlp_count_recv, & 165 trlp_count_recv_sum, trlpt_count, trlpt_count_recv, & 166 trnp_count, trnp_count_sum, trnp_count_recv, & 167 trnp_count_recv_sum, trnpt_count, trnpt_count_recv, & 168 trrp_count, trrp_count_sum, trrp_count_recv, & 169 trrp_count_recv_sum, trrpt_count, trrpt_count_recv, & 170 trsp_count, trsp_count_sum, trsp_count_recv, & 171 trsp_count_recv_sum, trspt_count, trspt_count_recv 172 173 INTEGER :: gp_outside_of_building(1:8) 174 175 INTEGER, PARAMETER :: maxtry = 40 ! for Rosenbrock method 176 177 LOGICAL :: dt_3d_reached, dt_3d_reached_l, prt_position 178 179 REAL :: aa, afactor, arg, bb, cc, dd, ddenom, delta_r, delta_v, & 180 dens_ratio, de_dt, de_dt_min, de_dx_int, de_dx_int_l, & 181 de_dx_int_u, de_dy_int, de_dy_int_l, de_dy_int_u, de_dz_int, & 182 de_dz_int_l, de_dz_int_u, diss_int, diss_int_l, diss_int_u, & 183 distance, drdt, drdt_ini, drdt_m, dt_ros, dt_ros_last, & 184 dt_ros_next, dt_ros_sum, dt_ros_sum_ini, d2rdt2, d2rdtdr, & 185 dt_gap, dt_particle, dt_particle_m, d_sum, epsilon, e_a, e_int,& 186 e_int_l, e_int_u, e_mean_int, e_s, err_ros, errmax, exp_arg, & 187 exp_term, fs_int, gg, g1, g2, g3, g4, integral, & 188 lagr_timescale, lw_max, mean_r, new_r, p_int, pt_int, & 189 pt_int_l, pt_int_u, q_int, q_int_l, q_int_u, ql_int, ql_int_l, & 190 ql_int_u, random_gauss, r_ros, r_ros_ini, & 191 sigma, sl_r3, sl_r4, t_int, u_int, u_int_l, u_int_u,vv_int, & 192 v_int, v_int_l, v_int_u, w_int, w_int_l, w_int_u, x, y 193 ! 194 !-- Parameters for Rosenbrock method 195 REAL, PARAMETER :: a21 = 2.0, a31 = 48.0/25.0, a32 = 6.0/25.0, & 196 a2x = 1.0, a3x = 3.0/5.0, b1 = 19.0/9.0, b2 = 0.5, & 197 b3 = 25.0/108.0, b4 = 125.0/108.0, c21 = -8.0, & 198 c31 = 372.0/25.0, c32 = 12.0/5.0, & 199 c41 = -112.0/125.0, c42 = -54.0/125.0, & 200 c43 = -2.0/5.0, c1x = 0.5, c2x= -3.0/2.0, & 201 c3x = 121.0/50.0, c4x = 29.0/250.0, & 202 errcon = 0.1296, e1 = 17.0/54.0, e2 = 7.0/36.0, & 203 e3 = 0.0, e4 = 125.0/108.0, gam = 0.5, grow = 1.5, & 204 pgrow = -0.25, pshrnk = -1.0/3.0, shrnk = 0.5 205 206 REAL, DIMENSION(1:30) :: de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei 207 208 REAL :: location(1:30,1:3) 209 210 REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: de_dx, de_dy, de_dz 211 212 REAL, DIMENSION(:,:,:), ALLOCATABLE :: trlpt, trnpt, trrpt, trspt 213 214 TYPE(particle_type) :: tmp_particle 215 216 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: trlp, trnp, trrp, trsp 217 218 219 CALL cpu_log( log_point(25), 'advec_particles', 'start' ) 220 221 222 ! 223 !-- Write particle data on file for later analysis. 224 !-- This has to be done here (before particles are advected) in order 225 !-- to allow correct output in case of dt_write_particle_data = dt_prel = 226 !-- particle_maximum_age. Otherwise (if output is done at the end of this 227 !-- subroutine), the relevant particles would have been already deleted. 228 !-- The MOD function allows for changes in the output interval with restart 229 !-- runs. 230 !-- Attention: change version number for unit 85 (in routine check_open) 231 !-- whenever the output format for this unit is changed! 62 INTEGER :: m 63 64 65 CALL cpu_log( log_point(25), 'lpm', 'start' ) 66 67 ! 68 !-- Write particle data at current time on file. 69 !-- This has to be done here, before particles are further processed, 70 !-- because they may be deleted within this timestep (in case that 71 !-- dt_write_particle_data = dt_prel = particle_maximum_age). 232 72 time_write_particle_data = time_write_particle_data + dt_3d 233 73 IF ( time_write_particle_data >= dt_write_particle_data ) THEN 234 74 235 CALL cpu_log( log_point_s(40), 'advec_part_io', 'start' ) 236 CALL check_open( 85 ) 237 WRITE ( 85 ) simulated_time, maximum_number_of_particles, & 238 number_of_particles 239 WRITE ( 85 ) particles 240 WRITE ( 85 ) maximum_number_of_tailpoints, maximum_number_of_tails, & 241 number_of_tails 242 IF ( maximum_number_of_tails > 0 ) THEN 243 WRITE ( 85 ) particle_tail_coordinates 244 ENDIF 245 CALL close_file( 85 ) 246 247 IF ( netcdf_output ) CALL output_particles_netcdf 248 75 CALL lpm_data_output_particles 76 ! 77 !-- The MOD function allows for changes in the output interval with restart 78 !-- runs. 249 79 time_write_particle_data = MOD( time_write_particle_data, & 250 80 MAX( dt_write_particle_data, dt_3d ) ) 251 CALL cpu_log( log_point_s(40), 'advec_part_io', 'stop' ) 252 ENDIF 253 254 ! 255 !-- Initialize variables for the (sub-) timestep, i.e. for 256 !-- marking those particles to be deleted after the timestep 257 particle_mask = .TRUE. 258 deleted_particles = 0 259 trlp_count_recv = 0 260 trnp_count_recv = 0 261 trrp_count_recv = 0 262 trsp_count_recv = 0 263 trlpt_count_recv = 0 264 trnpt_count_recv = 0 265 trrpt_count_recv = 0 266 trspt_count_recv = 0 267 IF ( use_particle_tails ) THEN 268 tail_mask = .TRUE. 269 ENDIF 270 deleted_tails = 0 271 272 ! 273 !-- Calculate exponential term used in case of particle inertia for each 274 !-- of the particle groups 275 CALL cpu_log( log_point_s(41), 'advec_part_exp', 'start' ) 276 DO m = 1, number_of_particle_groups 277 IF ( particle_groups(m)%density_ratio /= 0.0 ) THEN 278 particle_groups(m)%exp_arg = & 279 4.5 * particle_groups(m)%density_ratio * & 280 molecular_viscosity / ( particle_groups(m)%radius )**2 281 particle_groups(m)%exp_term = EXP( -particle_groups(m)%exp_arg * & 282 dt_3d ) 283 ENDIF 284 ENDDO 285 CALL cpu_log( log_point_s(41), 'advec_part_exp', 'stop' ) 286 287 ! 288 !-- Particle (droplet) growth by condensation/evaporation and collision 289 IF ( cloud_droplets ) THEN 290 291 ! 292 !-- Reset summation arrays 293 ql_c = 0.0; ql_v = 0.0; ql_vp = 0.0 294 delta_r=0.0; delta_v=0.0 295 296 ! 297 !-- Particle growth by condensation/evaporation 298 CALL cpu_log( log_point_s(42), 'advec_part_cond', 'start' ) 299 DO n = 1, number_of_particles 300 ! 301 !-- Interpolate temperature and humidity. 302 !-- First determine left, south, and bottom index of the arrays. 303 i = particles(n)%x * ddx 304 j = particles(n)%y * ddy 305 k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz & 306 + offset_ocean_nzt ! only exact if equidistant 307 308 x = particles(n)%x - i * dx 309 y = particles(n)%y - j * dy 310 aa = x**2 + y**2 311 bb = ( dx - x )**2 + y**2 312 cc = x**2 + ( dy - y )**2 313 dd = ( dx - x )**2 + ( dy - y )**2 314 gg = aa + bb + cc + dd 315 316 pt_int_l = ( ( gg - aa ) * pt(k,j,i) + ( gg - bb ) * pt(k,j,i+1) & 317 + ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) & 318 ) / ( 3.0 * gg ) 319 320 pt_int_u = ( ( gg-aa ) * pt(k+1,j,i) + ( gg-bb ) * pt(k+1,j,i+1) & 321 + ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) & 322 ) / ( 3.0 * gg ) 323 324 pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * & 325 ( pt_int_u - pt_int_l ) 326 327 q_int_l = ( ( gg - aa ) * q(k,j,i) + ( gg - bb ) * q(k,j,i+1) & 328 + ( gg - cc ) * q(k,j+1,i) + ( gg - dd ) * q(k,j+1,i+1) & 329 ) / ( 3.0 * gg ) 330 331 q_int_u = ( ( gg-aa ) * q(k+1,j,i) + ( gg-bb ) * q(k+1,j,i+1) & 332 + ( gg-cc ) * q(k+1,j+1,i) + ( gg-dd ) * q(k+1,j+1,i+1) & 333 ) / ( 3.0 * gg ) 334 335 q_int = q_int_l + ( particles(n)%z - zu(k) ) / dz * & 336 ( q_int_u - q_int_l ) 337 338 ql_int_l = ( ( gg - aa ) * ql(k,j,i) + ( gg - bb ) * ql(k,j,i+1) & 339 + ( gg - cc ) * ql(k,j+1,i) + ( gg - dd ) * ql(k,j+1,i+1) & 340 ) / ( 3.0 * gg ) 341 342 ql_int_u = ( ( gg-aa ) * ql(k+1,j,i) + ( gg-bb ) * ql(k+1,j,i+1) & 343 + ( gg-cc ) * ql(k+1,j+1,i) + ( gg-dd ) * ql(k+1,j+1,i+1) & 344 ) / ( 3.0 * gg ) 345 346 ql_int = ql_int_l + ( particles(n)%z - zu(k) ) / dz * & 347 ( ql_int_u - ql_int_l ) 348 349 ! 350 !-- Calculate real temperature and saturation vapor pressure 351 p_int = hyp(k) + ( particles(n)%z - zu(k) ) / dz * ( hyp(k+1)-hyp(k) ) 352 t_int = pt_int * ( p_int / 100000.0 )**0.286 353 354 e_s = 611.0 * EXP( l_d_rv * ( 3.6609E-3 - 1.0 / t_int ) ) 355 356 ! 357 !-- Current vapor pressure 358 e_a = q_int * p_int / ( 0.378 * q_int + 0.622 ) 359 360 ! 361 !-- Thermal conductivity for water (from Rogers and Yau, Table 7.1), 362 !-- diffusivity for water vapor (after Hall und Pruppacher, 1976) 363 thermal_conductivity_l = 7.94048E-05 * t_int + 0.00227011 364 diff_coeff_l = 0.211E-4 * ( t_int / 273.15 )**1.94 * & 365 ( 101325.0 / p_int) 366 ! 367 !-- Change in radius by condensation/evaporation 368 IF ( particles(n)%radius >= 1.0E-6 .OR. & 369 .NOT. curvature_solution_effects ) THEN 370 ! 371 !-- Approximation for large radii, where curvature and solution 372 !-- effects can be neglected 373 arg = particles(n)%radius**2 + 2.0 * dt_3d * & 374 ( e_a / e_s - 1.0 ) / & 375 ( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / & 376 thermal_conductivity_l + & 377 rho_l * r_v * t_int / diff_coeff_l / e_s ) 378 IF ( arg < 1.0E-16 ) THEN 379 new_r = 1.0E-8 380 ELSE 381 new_r = SQRT( arg ) 382 ENDIF 383 ENDIF 384 385 IF ( curvature_solution_effects .AND. & 386 ( ( particles(n)%radius < 1.0E-6 ) .OR. ( new_r < 1.0E-6 ) ) ) & 387 THEN 388 ! 389 !-- Curvature and solutions effects are included in growth equation. 390 !-- Change in Radius is calculated with 4th-order Rosenbrock method 391 !-- for stiff o.d.e's with monitoring local truncation error to adjust 392 !-- stepsize (see Numerical recipes in FORTRAN, 2nd edition, p. 731). 393 !-- For larger radii the simple analytic method (see ELSE) gives 394 !-- almost the same results. 395 ! 396 !-- Surface tension after (Straka, 2009) 397 sigma = 0.0761 - 0.00155 * ( t_int - 273.15 ) 398 399 r_ros = particles(n)%radius 400 dt_ros_sum = 0.0 ! internal integrated time (s) 401 internal_timestep_count = 0 402 403 ddenom = 1.0 / ( rho_l * r_v * t_int / ( e_s * diff_coeff_l ) + & 404 ( l_v / ( r_v * t_int ) - 1.0 ) * & 405 rho_l * l_v / ( thermal_conductivity_l * t_int )& 406 ) 407 408 afactor = 2.0 * sigma / ( rho_l * r_v * t_int ) 409 410 IF ( particles(n)%rvar3 == -9999999.9 ) THEN 411 ! 412 !-- First particle timestep. Derivative has to be calculated. 413 drdt_m = ddenom / r_ros * ( e_a / e_s - 1.0 - & 414 afactor / r_ros + & 415 bfactor / r_ros**3 ) 416 ELSE 417 ! 418 !-- Take value from last PALM timestep 419 drdt_m = particles(n)%rvar3 420 ENDIF 421 ! 422 !-- Take internal timestep values from the end of last PALM timestep 423 dt_ros_last = particles(n)%rvar1 424 dt_ros_next = particles(n)%rvar2 425 ! 426 !-- Internal timestep must not be larger than PALM timestep 427 dt_ros = MIN( dt_ros_next, dt_3d ) 428 ! 429 !-- Integrate growth equation in time unless PALM timestep is reached 430 DO WHILE ( dt_ros_sum < dt_3d ) 431 432 internal_timestep_count = internal_timestep_count + 1 433 434 ! 435 !-- Derivative at starting value 436 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - afactor / r_ros + & 437 bfactor / r_ros**3 ) 438 drdt_ini = drdt 439 dt_ros_sum_ini = dt_ros_sum 440 r_ros_ini = r_ros 441 442 ! 443 !-- Calculate time derivative of dr/dt 444 d2rdt2 = ( drdt - drdt_m ) / dt_ros_last 445 446 ! 447 !-- Calculate radial derivative of dr/dt 448 d2rdtdr = ddenom * ( ( 1.0 - e_a/e_s ) / r_ros**2 + & 449 2.0 * afactor / r_ros**3 - & 450 4.0 * bfactor / r_ros**5 ) 451 ! 452 !-- Adjust stepsize unless required accuracy is reached 453 DO jtry = 1, maxtry+1 454 455 IF ( jtry == maxtry+1 ) THEN 456 message_string = 'maxtry > 40 in Rosenbrock method' 457 CALL message( 'advec_particles', 'PA0347', 2, 2, 0, 6, 0 ) 458 ENDIF 459 460 aa = 1.0 / ( gam * dt_ros ) - d2rdtdr 461 g1 = ( drdt_ini + dt_ros * c1x * d2rdt2 ) / aa 462 r_ros = r_ros_ini + a21 * g1 463 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - & 464 afactor / r_ros + & 465 bfactor / r_ros**3 ) 466 467 g2 = ( drdt + dt_ros * c2x * d2rdt2 + c21 * g1 / dt_ros )& 468 / aa 469 r_ros = r_ros_ini + a31 * g1 + a32 * g2 470 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - & 471 afactor / r_ros + & 472 bfactor / r_ros**3 ) 473 474 g3 = ( drdt + dt_ros * c3x * d2rdt2 + & 475 ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa 476 g4 = ( drdt + dt_ros * c4x * d2rdt2 + & 477 ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa 478 r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4 479 480 dt_ros_sum = dt_ros_sum_ini + dt_ros 481 482 IF ( dt_ros_sum == dt_ros_sum_ini ) THEN 483 message_string = 'zero stepsize in Rosenbrock method' 484 CALL message( 'advec_particles', 'PA0348', 2, 2, 0, 6, 0 ) 485 ENDIF 486 ! 487 !-- Calculate error 488 err_ros = e1*g1 + e2*g2 + e3*g3 + e4*g4 489 errmax = 0.0 490 errmax = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros 491 ! 492 !-- Leave loop if accuracy is sufficient, otherwise try again 493 !-- with a reduced stepsize 494 IF ( errmax <= 1.0 ) THEN 495 EXIT 496 ELSE 497 dt_ros = SIGN( & 498 MAX( ABS( 0.9 * dt_ros * errmax**pshrnk ), & 499 shrnk * ABS( dt_ros ) ), & 500 dt_ros & 501 ) 502 ENDIF 503 504 ENDDO ! loop for stepsize adjustment 505 506 ! 507 !-- Calculate next internal timestep 508 IF ( errmax > errcon ) THEN 509 dt_ros_next = 0.9 * dt_ros * errmax**pgrow 510 ELSE 511 dt_ros_next = grow * dt_ros 512 ENDIF 513 514 ! 515 !-- Estimated timestep is reduced if the PALM time step is exceeded 516 dt_ros_last = dt_ros 517 IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d ) THEN 518 dt_ros = dt_3d - dt_ros_sum 519 ELSE 520 dt_ros = dt_ros_next 521 ENDIF 522 523 drdt_m = drdt 524 525 ENDDO 526 ! 527 !-- Store derivative and internal timestep values for next PALM step 528 particles(n)%rvar1 = dt_ros_last 529 particles(n)%rvar2 = dt_ros_next 530 particles(n)%rvar3 = drdt_m 531 532 new_r = r_ros 533 ! 534 !-- Radius should not fall below 1E-8 because Rosenbrock method may 535 !-- lead to errors otherwise 536 new_r = MAX( new_r, 1.0E-8 ) 537 538 ENDIF 539 540 delta_r = new_r - particles(n)%radius 541 542 ! 543 !-- Sum up the change in volume of liquid water for the respective grid 544 !-- volume (this is needed later on for calculating the release of 545 !-- latent heat) 546 i = ( particles(n)%x + 0.5 * dx ) * ddx 547 j = ( particles(n)%y + 0.5 * dy ) * ddy 548 k = particles(n)%z / dz + 1 + offset_ocean_nzt_m1 549 ! only exact if equidistant 550 551 ql_c(k,j,i) = ql_c(k,j,i) + particles(n)%weight_factor * & 552 rho_l * 1.33333333 * pi * & 553 ( new_r**3 - particles(n)%radius**3 ) / & 554 ( rho_surface * dx * dy * dz ) 555 IF ( ql_c(k,j,i) > 100.0 ) THEN 556 WRITE( message_string, * ) 'k=',k,' j=',j,' i=',i, & 557 ' ql_c=',ql_c(k,j,i), ' &part(',n,')%wf=', & 558 particles(n)%weight_factor,' delta_r=',delta_r 559 CALL message( 'advec_particles', 'PA0143', 2, 2, -1, 6, 1 ) 560 ENDIF 561 562 ! 563 !-- Change the droplet radius 564 IF ( ( new_r - particles(n)%radius ) < 0.0 .AND. new_r < 0.0 ) & 565 THEN 566 WRITE( message_string, * ) '#1 k=',k,' j=',j,' i=',i, & 567 ' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int, & 568 ' &delta_r=',delta_r, & 569 ' particle_radius=',particles(n)%radius 570 CALL message( 'advec_particles', 'PA0144', 2, 2, -1, 6, 1 ) 571 ENDIF 572 573 ! 574 !-- Sum up the total volume of liquid water (needed below for 575 !-- re-calculating the weighting factors) 576 ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * new_r**3 577 578 particles(n)%radius = new_r 579 580 ! 581 !-- Determine radius class of the particle needed for collision 582 IF ( ( hall_kernel .OR. wang_kernel ) .AND. use_kernel_tables ) & 583 THEN 584 particles(n)%class = ( LOG( new_r ) - rclass_lbound ) / & 585 ( rclass_ubound - rclass_lbound ) * & 586 radius_classes 587 particles(n)%class = MIN( particles(n)%class, radius_classes ) 588 particles(n)%class = MAX( particles(n)%class, 1 ) 589 ENDIF 590 591 ENDDO 592 CALL cpu_log( log_point_s(42), 'advec_part_cond', 'stop' ) 593 594 ! 595 !-- Particle growth by collision 596 IF ( collision_kernel /= 'none' ) THEN 597 598 CALL cpu_log( log_point_s(43), 'advec_part_coll', 'start' ) 599 600 DO i = nxl, nxr 601 DO j = nys, nyn 602 DO k = nzb+1, nzt 603 ! 604 !-- Collision requires at least two particles in the box 605 IF ( prt_count(k,j,i) > 1 ) THEN 606 ! 607 !-- First, sort particles within the gridbox by their size, 608 !-- using Shell's method (see Numerical Recipes) 609 !-- NOTE: In case of using particle tails, the re-sorting of 610 !-- ---- tails would have to be included here! 611 psi = prt_start_index(k,j,i) - 1 612 inc = 1 613 DO WHILE ( inc <= prt_count(k,j,i) ) 614 inc = 3 * inc + 1 615 ENDDO 616 617 DO WHILE ( inc > 1 ) 618 inc = inc / 3 619 DO is = inc+1, prt_count(k,j,i) 620 tmp_particle = particles(psi+is) 621 js = is 622 DO WHILE ( particles(psi+js-inc)%radius > & 623 tmp_particle%radius ) 624 particles(psi+js) = particles(psi+js-inc) 625 js = js - inc 626 IF ( js <= inc ) EXIT 627 ENDDO 628 particles(psi+js) = tmp_particle 629 ENDDO 630 ENDDO 631 632 psi = prt_start_index(k,j,i) 633 pse = psi + prt_count(k,j,i)-1 634 635 ! 636 !-- Now apply the different kernels 637 IF ( ( hall_kernel .OR. wang_kernel ) .AND. & 638 use_kernel_tables ) THEN 639 ! 640 !-- Fast method with pre-calculated efficiencies for 641 !-- discrete radius- and dissipation-classes. 642 ! 643 !-- Determine dissipation class index of this gridbox 644 IF ( wang_kernel ) THEN 645 eclass = INT( diss(k,j,i) * 1.0E4 / 1000.0 * & 646 dissipation_classes ) + 1 647 epsilon = diss(k,j,i) 648 ELSE 649 epsilon = 0.0 650 ENDIF 651 IF ( hall_kernel .OR. epsilon * 1.0E4 < 0.001 ) THEN 652 eclass = 0 ! Hall kernel is used 653 ELSE 654 eclass = MIN( dissipation_classes, eclass ) 655 ENDIF 656 657 DO n = pse, psi+1, -1 658 659 integral = 0.0 660 lw_max = 0.0 661 rclass_l = particles(n)%class 662 ! 663 !-- Calculate growth of collector particle radius using all 664 !-- droplets smaller than current droplet 665 DO is = psi, n-1 666 667 rclass_s = particles(is)%class 668 integral = integral + & 669 ( particles(is)%radius**3 * & 670 ckernel(rclass_l,rclass_s,eclass) * & 671 particles(is)%weight_factor ) 672 ! 673 !-- Calculate volume of liquid water of the collected 674 !-- droplets which is the maximum liquid water available 675 !-- for droplet growth 676 lw_max = lw_max + ( particles(is)%radius**3 * & 677 particles(is)%weight_factor ) 678 679 ENDDO 680 681 ! 682 !-- Change in radius of collector droplet due to collision 683 delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * & 684 integral * dt_3d * ddx * ddy / dz 685 686 ! 687 !-- Change in volume of collector droplet due to collision 688 delta_v = particles(n)%weight_factor & 689 * ( ( particles(n)%radius + delta_r )**3 & 690 - particles(n)%radius**3 ) 691 692 IF ( lw_max < delta_v .AND. delta_v > 0.0 ) THEN 693 !-- replace by message call 694 PRINT*, 'Particle has grown to much because the', & 695 ' change of volume of particle is larger', & 696 ' than liquid water available!' 697 698 ELSEIF ( lw_max == delta_v .AND. delta_v > 0.0 ) THEN 699 !-- can this case really happen?? 700 DO is = psi, n-1 701 702 particles(is)%weight_factor = 0.0 703 particle_mask(is) = .FALSE. 704 deleted_particles = deleted_particles + 1 705 706 ENDDO 707 708 ELSEIF ( lw_max > delta_v .AND. delta_v > 0.0 ) THEN 709 ! 710 !-- Calculate new weighting factor of collected droplets 711 DO is = psi, n-1 712 713 rclass_s = particles(is)%class 714 particles(is)%weight_factor = & 715 particles(is)%weight_factor - & 716 ( ( ckernel(rclass_l,rclass_s,eclass) * particles(is)%weight_factor & 717 / integral ) * delta_v ) 718 719 IF ( particles(is)%weight_factor < 0.0 ) THEN 720 WRITE( message_string, * ) 'negative ', & 721 'weighting factor: ', & 722 particles(is)%weight_factor 723 CALL message( 'advec_particles', '', 2, 2,& 724 -1, 6, 1 ) 725 726 ELSEIF ( particles(is)%weight_factor == 0.0 ) & 727 THEN 728 729 particles(is)%weight_factor = 0.0 730 particle_mask(is) = .FALSE. 731 deleted_particles = deleted_particles + 1 732 733 ENDIF 734 735 ENDDO 736 737 ENDIF 738 739 particles(n)%radius = particles(n)%radius + delta_r 740 ql_vp(k,j,i) = ql_vp(k,j,i) + & 741 particles(n)%weight_factor & 742 * particles(n)%radius**3 743 744 ENDDO 745 746 ELSEIF ( ( hall_kernel .OR. wang_kernel ) .AND. & 747 .NOT. use_kernel_tables ) THEN 748 ! 749 !-- Collision efficiencies are calculated for every new 750 !-- grid box. First, allocate memory for kernel table. 751 !-- Third dimension is 1, because table is re-calculated for 752 !-- every new dissipation value. 753 ALLOCATE( ckernel(prt_start_index(k,j,i): & 754 prt_start_index(k,j,i)+prt_count(k,j,i)-1, & 755 prt_start_index(k,j,i): & 756 prt_start_index(k,j,i)+prt_count(k,j,i)-1,1:1) ) 757 ! 758 !-- Now calculate collision efficiencies for this box 759 CALL recalculate_kernel( i, j, k ) 760 761 DO n = pse, psi+1, -1 762 763 integral = 0.0 764 lw_max = 0.0 765 ! 766 !-- Calculate growth of collector particle radius using all 767 !-- droplets smaller than current droplet 768 DO is = psi, n-1 769 770 integral = integral + particles(is)%radius**3 * & 771 ckernel(n,is,1) * & 772 particles(is)%weight_factor 773 ! 774 !-- Calculate volume of liquid water of the collected 775 !-- droplets which is the maximum liquid water available 776 !-- for droplet growth 777 lw_max = lw_max + ( particles(is)%radius**3 * & 778 particles(is)%weight_factor ) 779 780 ENDDO 781 782 ! 783 !-- Change in radius of collector droplet due to collision 784 delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * & 785 integral * dt_3d * ddx * ddy / dz 786 787 ! 788 !-- Change in volume of collector droplet due to collision 789 delta_v = particles(n)%weight_factor & 790 * ( ( particles(n)%radius + delta_r )**3 & 791 - particles(n)%radius**3 ) 792 793 IF ( lw_max < delta_v .AND. delta_v > 0.0 ) THEN 794 !-- replace by message call 795 PRINT*, 'Particle has grown to much because the', & 796 ' change of volume of particle is larger', & 797 ' than liquid water available!' 798 799 ELSEIF ( lw_max == delta_v .AND. delta_v > 0.0 ) THEN 800 !-- can this case really happen?? 801 DO is = psi, n-1 802 803 particles(is)%weight_factor = 0.0 804 particle_mask(is) = .FALSE. 805 deleted_particles = deleted_particles + 1 806 807 ENDDO 808 809 ELSEIF ( lw_max > delta_v .AND. delta_v > 0.0 ) THEN 810 ! 811 !-- Calculate new weighting factor of collected droplets 812 DO is = psi, n-1 813 814 particles(is)%weight_factor = & 815 particles(is)%weight_factor - & 816 ( ckernel(n,is,1) / integral * delta_v * & 817 particles(is)%weight_factor ) 818 819 IF ( particles(is)%weight_factor < 0.0 ) THEN 820 WRITE( message_string, * ) 'negative ', & 821 'weighting factor: ', & 822 particles(is)%weight_factor 823 CALL message( 'advec_particles', '', 2, 2,& 824 -1, 6, 1 ) 825 826 ELSEIF ( particles(is)%weight_factor == 0.0 ) & 827 THEN 828 829 particles(is)%weight_factor = 0.0 830 particle_mask(is) = .FALSE. 831 deleted_particles = deleted_particles + 1 832 833 ENDIF 834 835 ENDDO 836 837 ENDIF 838 839 particles(n)%radius = particles(n)%radius + delta_r 840 ql_vp(k,j,i) = ql_vp(k,j,i) + & 841 particles(n)%weight_factor & 842 * particles(n)%radius**3 843 844 ENDDO 845 846 DEALLOCATE( ckernel ) 847 848 ELSEIF ( palm_kernel ) THEN 849 ! 850 !-- PALM collision kernel 851 ! 852 !-- Calculate the mean radius of all those particles which 853 !-- are of smaller size than the current particle and 854 !-- use this radius for calculating the collision efficiency 855 DO n = psi+prt_count(k,j,i)-1, psi+1, -1 856 857 sl_r3 = 0.0 858 sl_r4 = 0.0 859 860 DO is = n-1, psi, -1 861 IF ( particles(is)%radius < particles(n)%radius ) & 862 THEN 863 sl_r3 = sl_r3 + particles(is)%weight_factor & 864 *( particles(is)%radius**3 ) 865 sl_r4 = sl_r4 + particles(is)%weight_factor & 866 *( particles(is)%radius**4 ) 867 ENDIF 868 ENDDO 869 870 IF ( ( sl_r3 ) > 0.0 ) THEN 871 mean_r = ( sl_r4 ) / ( sl_r3 ) 872 873 CALL collision_efficiency( mean_r, & 874 particles(n)%radius, & 875 effective_coll_efficiency ) 876 877 ELSE 878 effective_coll_efficiency = 0.0 879 ENDIF 880 881 IF ( effective_coll_efficiency > 1.0 .OR. & 882 effective_coll_efficiency < 0.0 ) & 883 THEN 884 WRITE( message_string, * ) 'collision_efficien' , & 885 'cy out of range:' ,effective_coll_efficiency 886 CALL message( 'advec_particles', 'PA0145', 2, 2, & 887 -1, 6, 1 ) 888 ENDIF 889 890 ! 891 !-- Interpolation of ... 892 ii = particles(n)%x * ddx 893 jj = particles(n)%y * ddy 894 kk = ( particles(n)%z + 0.5 * dz ) / dz 895 896 x = particles(n)%x - ii * dx 897 y = particles(n)%y - jj * dy 898 aa = x**2 + y**2 899 bb = ( dx - x )**2 + y**2 900 cc = x**2 + ( dy - y )**2 901 dd = ( dx - x )**2 + ( dy - y )**2 902 gg = aa + bb + cc + dd 903 904 ql_int_l = ( (gg-aa) * ql(kk,jj,ii) + (gg-bb) * & 905 ql(kk,jj,ii+1) & 906 + (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) * & 907 ql(kk,jj+1,ii+1) & 908 ) / ( 3.0 * gg ) 909 910 ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii) + (gg-bb) * & 911 ql(kk+1,jj,ii+1) & 912 + (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) * & 913 ql(kk+1,jj+1,ii+1) & 914 ) / ( 3.0 * gg ) 915 916 ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *& 917 ( ql_int_u - ql_int_l ) 918 919 ! 920 !-- Interpolate u velocity-component 921 ii = ( particles(n)%x + 0.5 * dx ) * ddx 922 jj = particles(n)%y * ddy 923 kk = ( particles(n)%z + 0.5 * dz ) / dz ! only if eqist 924 925 IF ( ( particles(n)%z - zu(kk) ) > (0.5*dz) ) kk = kk+1 926 927 x = particles(n)%x + ( 0.5 - ii ) * dx 928 y = particles(n)%y - jj * dy 929 aa = x**2 + y**2 930 bb = ( dx - x )**2 + y**2 931 cc = x**2 + ( dy - y )**2 932 dd = ( dx - x )**2 + ( dy - y )**2 933 gg = aa + bb + cc + dd 934 935 u_int_l = ( (gg-aa) * u(kk,jj,ii) + (gg-bb) * & 936 u(kk,jj,ii+1) & 937 + (gg-cc) * u(kk,jj+1,ii) + (gg-dd) * & 938 u(kk,jj+1,ii+1) & 939 ) / ( 3.0 * gg ) - u_gtrans 940 IF ( kk+1 == nzt+1 ) THEN 941 u_int = u_int_l 942 ELSE 943 u_int_u = ( (gg-aa) * u(kk+1,jj,ii) + (gg-bb) * & 944 u(kk+1,jj,ii+1) & 945 + (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) * & 946 u(kk+1,jj+1,ii+1) & 947 ) / ( 3.0 * gg ) - u_gtrans 948 u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz & 949 * ( u_int_u - u_int_l ) 950 ENDIF 951 952 ! 953 !-- Same procedure for interpolation of the v velocity-com- 954 !-- ponent (adopt index k from u velocity-component) 955 ii = particles(n)%x * ddx 956 jj = ( particles(n)%y + 0.5 * dy ) * ddy 957 958 x = particles(n)%x - ii * dx 959 y = particles(n)%y + ( 0.5 - jj ) * dy 960 aa = x**2 + y**2 961 bb = ( dx - x )**2 + y**2 962 cc = x**2 + ( dy - y )**2 963 dd = ( dx - x )**2 + ( dy - y )**2 964 gg = aa + bb + cc + dd 965 966 v_int_l = ( ( gg-aa ) * v(kk,jj,ii) + ( gg-bb ) * & 967 v(kk,jj,ii+1) & 968 + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) * & 969 v(kk,jj+1,ii+1) & 970 ) / ( 3.0 * gg ) - v_gtrans 971 IF ( kk+1 == nzt+1 ) THEN 972 v_int = v_int_l 973 ELSE 974 v_int_u = ( (gg-aa) * v(kk+1,jj,ii) + (gg-bb) * & 975 v(kk+1,jj,ii+1) & 976 + (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) * & 977 v(kk+1,jj+1,ii+1) & 978 ) / ( 3.0 * gg ) - v_gtrans 979 v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz & 980 * ( v_int_u - v_int_l ) 981 ENDIF 982 983 ! 984 !-- Same procedure for interpolation of the w velocity-com- 985 !-- ponent (adopt index i from v velocity-component) 986 jj = particles(n)%y * ddy 987 kk = particles(n)%z / dz 988 989 x = particles(n)%x - ii * dx 990 y = particles(n)%y - jj * dy 991 aa = x**2 + y**2 992 bb = ( dx - x )**2 + y**2 993 cc = x**2 + ( dy - y )**2 994 dd = ( dx - x )**2 + ( dy - y )**2 995 gg = aa + bb + cc + dd 996 997 w_int_l = ( ( gg-aa ) * w(kk,jj,ii) + ( gg-bb ) * & 998 w(kk,jj,ii+1) & 999 + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) * & 1000 w(kk,jj+1,ii+1) & 1001 ) / ( 3.0 * gg ) 1002 IF ( kk+1 == nzt+1 ) THEN 1003 w_int = w_int_l 1004 ELSE 1005 w_int_u = ( (gg-aa) * w(kk+1,jj,ii) + (gg-bb) * & 1006 w(kk+1,jj,ii+1) & 1007 + (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) * & 1008 w(kk+1,jj+1,ii+1) & 1009 ) / ( 3.0 * gg ) 1010 w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz & 1011 * ( w_int_u - w_int_l ) 1012 ENDIF 1013 1014 ! 1015 !-- Change in radius due to collision 1016 delta_r = effective_coll_efficiency / 3.0 & 1017 * pi * sl_r3 * ddx * ddy / dz & 1018 * SQRT( ( u_int - particles(n)%speed_x )**2 & 1019 + ( v_int - particles(n)%speed_y )**2 & 1020 + ( w_int - particles(n)%speed_z )**2 & 1021 ) * dt_3d 1022 ! 1023 !-- Change in volume due to collision 1024 delta_v = particles(n)%weight_factor & 1025 * ( ( particles(n)%radius + delta_r )**3 & 1026 - particles(n)%radius**3 ) 1027 1028 ! 1029 !-- Check if collected particles provide enough LWC for 1030 !-- volume change of collector particle 1031 IF ( delta_v >= sl_r3 .AND. sl_r3 > 0.0 ) THEN 1032 1033 delta_r = ( ( sl_r3/particles(n)%weight_factor ) & 1034 + particles(n)%radius**3 )**( 1./3. ) & 1035 - particles(n)%radius 1036 1037 DO is = n-1, psi, -1 1038 IF ( particles(is)%radius < & 1039 particles(n)%radius ) THEN 1040 particles(is)%weight_factor = 0.0 1041 particle_mask(is) = .FALSE. 1042 deleted_particles = deleted_particles + 1 1043 ENDIF 1044 ENDDO 1045 1046 ELSE IF ( delta_v < sl_r3 .AND. sl_r3 > 0.0 ) THEN 1047 1048 DO is = n-1, psi, -1 1049 IF ( particles(is)%radius < particles(n)%radius & 1050 .AND. sl_r3 > 0.0 ) THEN 1051 particles(is)%weight_factor = & 1052 ( ( particles(is)%weight_factor & 1053 * ( particles(is)%radius**3 ) ) & 1054 - ( delta_v & 1055 * particles(is)%weight_factor & 1056 * ( particles(is)%radius**3 ) & 1057 / sl_r3 ) ) & 1058 / ( particles(is)%radius**3 ) 1059 1060 IF ( particles(is)%weight_factor < 0.0 ) THEN 1061 WRITE( message_string, * ) 'negative ', & 1062 'weighting factor: ', & 1063 particles(is)%weight_factor 1064 CALL message( 'advec_particles', '', 2, & 1065 2, -1, 6, 1 ) 1066 ENDIF 1067 ENDIF 1068 ENDDO 1069 ENDIF 1070 1071 particles(n)%radius = particles(n)%radius + delta_r 1072 ql_vp(k,j,i) = ql_vp(k,j,i) + & 1073 particles(n)%weight_factor * & 1074 ( particles(n)%radius**3 ) 1075 ENDDO 1076 1077 ENDIF ! collision kernel 1078 1079 ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor & 1080 * particles(psi)%radius**3 1081 1082 1083 ELSE IF ( prt_count(k,j,i) == 1 ) THEN 1084 1085 psi = prt_start_index(k,j,i) 1086 ql_vp(k,j,i) = particles(psi)%weight_factor * & 1087 particles(psi)%radius**3 1088 ENDIF 1089 1090 ! 1091 !-- Check if condensation of LWC was conserved during collision 1092 !-- process 1093 IF ( ql_v(k,j,i) /= 0.0 ) THEN 1094 IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001 .OR. & 1095 ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999 ) THEN 1096 WRITE( message_string, * ) 'LWC is not conserved during',& 1097 ' collision! ', & 1098 'LWC after condensation: ', & 1099 ql_v(k,j,i), & 1100 ' LWC after collision: ', & 1101 ql_vp(k,j,i) 1102 CALL message( 'advec_particles', '', 2, 2, -1, 6, & 1103 1 ) 1104 ENDIF 1105 ENDIF 1106 1107 ENDDO 1108 ENDDO 1109 ENDDO 1110 1111 ENDIF ! collision handling 1112 1113 CALL cpu_log( log_point_s(43), 'advec_part_coll', 'stop' ) 1114 1115 ENDIF ! cloud droplet handling 1116 1117 1118 ! 1119 !-- Particle advection. 1120 !-- In case of including the SGS velocities, the LES timestep has probably 1121 !-- to be split into several smaller timesteps because of the Lagrangian 1122 !-- timescale condition. Because the number of timesteps to be carried out is 1123 !-- not known at the beginning, these steps are carried out in an infinite loop 1124 !-- with exit condition. 1125 ! 1126 !-- If SGS velocities are used, gradients of the TKE have to be calculated and 1127 !-- boundary conditions have to be set first. Also, horizontally averaged 1128 !-- profiles of the SGS TKE and the resolved-scale velocity variances are 1129 !-- needed. 1130 IF ( use_sgs_for_particles ) THEN 1131 1132 ! 1133 !-- TKE gradient along x and y 1134 DO i = nxl, nxr 1135 DO j = nys, nyn 1136 DO k = nzb, nzt+1 1137 1138 IF ( k <= nzb_s_inner(j,i-1) .AND. & 1139 k > nzb_s_inner(j,i) .AND. & 1140 k > nzb_s_inner(j,i+1) ) THEN 1141 de_dx(k,j,i) = 2.0 * sgs_wfu_part * & 1142 ( e(k,j,i+1) - e(k,j,i) ) * ddx 1143 ELSEIF ( k > nzb_s_inner(j,i-1) .AND. & 1144 k > nzb_s_inner(j,i) .AND. & 1145 k <= nzb_s_inner(j,i+1) ) THEN 1146 de_dx(k,j,i) = 2.0 * sgs_wfu_part * & 1147 ( e(k,j,i) - e(k,j,i-1) ) * ddx 1148 ELSEIF ( k < nzb_s_inner(j,i) .AND. k < nzb_s_inner(j,i+1) ) & 1149 THEN 1150 de_dx(k,j,i) = 0.0 1151 ELSEIF ( k < nzb_s_inner(j,i-1) .AND. k < nzb_s_inner(j,i) ) & 1152 THEN 1153 de_dx(k,j,i) = 0.0 1154 ELSE 1155 de_dx(k,j,i) = sgs_wfu_part * & 1156 ( e(k,j,i+1) - e(k,j,i-1) ) * ddx 1157 ENDIF 1158 1159 IF ( k <= nzb_s_inner(j-1,i) .AND. & 1160 k > nzb_s_inner(j,i) .AND. & 1161 k > nzb_s_inner(j+1,i) ) THEN 1162 de_dy(k,j,i) = 2.0 * sgs_wfv_part * & 1163 ( e(k,j+1,i) - e(k,j,i) ) * ddy 1164 ELSEIF ( k > nzb_s_inner(j-1,i) .AND. & 1165 k > nzb_s_inner(j,i) .AND. & 1166 k <= nzb_s_inner(j+1,i) ) THEN 1167 de_dy(k,j,i) = 2.0 * sgs_wfv_part * & 1168 ( e(k,j,i) - e(k,j-1,i) ) * ddy 1169 ELSEIF ( k < nzb_s_inner(j,i) .AND. k < nzb_s_inner(j+1,i) ) & 1170 THEN 1171 de_dy(k,j,i) = 0.0 1172 ELSEIF ( k < nzb_s_inner(j-1,i) .AND. k < nzb_s_inner(j,i) ) & 1173 THEN 1174 de_dy(k,j,i) = 0.0 1175 ELSE 1176 de_dy(k,j,i) = sgs_wfv_part * & 1177 ( e(k,j+1,i) - e(k,j-1,i) ) * ddy 1178 ENDIF 1179 1180 ENDDO 1181 ENDDO 1182 ENDDO 1183 1184 ! 1185 !-- TKE gradient along z, including bottom and top boundary conditions 1186 DO i = nxl, nxr 1187 DO j = nys, nyn 1188 1189 DO k = nzb_s_inner(j,i)+2, nzt-1 1190 de_dz(k,j,i) = 2.0 * sgs_wfw_part * & 1191 ( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1)-zu(k-1) ) 1192 ENDDO 1193 1194 k = nzb_s_inner(j,i) 1195 de_dz(nzb:k,j,i) = 0.0 1196 de_dz(k+1,j,i) = 2.0 * sgs_wfw_part * ( e(k+2,j,i) - e(k+1,j,i) ) & 1197 / ( zu(k+2) - zu(k+1) ) 1198 de_dz(nzt,j,i) = 0.0 1199 de_dz(nzt+1,j,i) = 0.0 1200 ENDDO 1201 ENDDO 1202 1203 ! 1204 !-- Lateral boundary conditions 1205 CALL exchange_horiz( de_dx, nbgp ) 1206 CALL exchange_horiz( de_dy, nbgp ) 1207 CALL exchange_horiz( de_dz, nbgp ) 1208 CALL exchange_horiz( diss, nbgp ) 1209 1210 ! 1211 !-- Calculate the horizontally averaged profiles of SGS TKE and resolved 1212 !-- velocity variances (they may have been already calculated in routine 1213 !-- flow_statistics). 1214 IF ( .NOT. flow_statistics_called ) THEN 1215 ! 1216 !-- First calculate horizontally averaged profiles of the horizontal 1217 !-- velocities. 1218 sums_l(:,1,0) = 0.0 1219 sums_l(:,2,0) = 0.0 1220 1221 DO i = nxl, nxr 1222 DO j = nys, nyn 1223 DO k = nzb_s_outer(j,i), nzt+1 1224 sums_l(k,1,0) = sums_l(k,1,0) + u(k,j,i) 1225 sums_l(k,2,0) = sums_l(k,2,0) + v(k,j,i) 1226 ENDDO 1227 ENDDO 1228 ENDDO 1229 1230 #if defined( __parallel ) 1231 ! 1232 !-- Compute total sum from local sums 1233 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1234 CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, & 1235 MPI_REAL, MPI_SUM, comm2d, ierr ) 1236 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1237 CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, & 1238 MPI_REAL, MPI_SUM, comm2d, ierr ) 1239 #else 1240 sums(:,1) = sums_l(:,1,0) 1241 sums(:,2) = sums_l(:,2,0) 1242 #endif 1243 1244 ! 1245 !-- Final values are obtained by division by the total number of grid 1246 !-- points used for the summation. 1247 hom(:,1,1,0) = sums(:,1) / ngp_2dh_outer(:,0) ! u 1248 hom(:,1,2,0) = sums(:,2) / ngp_2dh_outer(:,0) ! v 1249 1250 ! 1251 !-- Now calculate the profiles of SGS TKE and the resolved-scale 1252 !-- velocity variances 1253 sums_l(:,8,0) = 0.0 1254 sums_l(:,30,0) = 0.0 1255 sums_l(:,31,0) = 0.0 1256 sums_l(:,32,0) = 0.0 1257 DO i = nxl, nxr 1258 DO j = nys, nyn 1259 DO k = nzb_s_outer(j,i), nzt+1 1260 sums_l(k,8,0) = sums_l(k,8,0) + e(k,j,i) 1261 sums_l(k,30,0) = sums_l(k,30,0) + & 1262 ( u(k,j,i) - hom(k,1,1,0) )**2 1263 sums_l(k,31,0) = sums_l(k,31,0) + & 1264 ( v(k,j,i) - hom(k,1,2,0) )**2 1265 sums_l(k,32,0) = sums_l(k,32,0) + w(k,j,i)**2 1266 ENDDO 1267 ENDDO 1268 ENDDO 1269 1270 #if defined( __parallel ) 1271 ! 1272 !-- Compute total sum from local sums 1273 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1274 CALL MPI_ALLREDUCE( sums_l(nzb,8,0), sums(nzb,8), nzt+2-nzb, & 1275 MPI_REAL, MPI_SUM, comm2d, ierr ) 1276 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1277 CALL MPI_ALLREDUCE( sums_l(nzb,30,0), sums(nzb,30), nzt+2-nzb, & 1278 MPI_REAL, MPI_SUM, comm2d, ierr ) 1279 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1280 CALL MPI_ALLREDUCE( sums_l(nzb,31,0), sums(nzb,31), nzt+2-nzb, & 1281 MPI_REAL, MPI_SUM, comm2d, ierr ) 1282 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1283 CALL MPI_ALLREDUCE( sums_l(nzb,32,0), sums(nzb,32), nzt+2-nzb, & 1284 MPI_REAL, MPI_SUM, comm2d, ierr ) 1285 1286 #else 1287 sums(:,8) = sums_l(:,8,0) 1288 sums(:,30) = sums_l(:,30,0) 1289 sums(:,31) = sums_l(:,31,0) 1290 sums(:,32) = sums_l(:,32,0) 1291 #endif 1292 1293 ! 1294 !-- Final values are obtained by division by the total number of grid 1295 !-- points used for the summation. 1296 hom(:,1,8,0) = sums(:,8) / ngp_2dh_outer(:,0) ! e 1297 hom(:,1,30,0) = sums(:,30) / ngp_2dh_outer(:,0) ! u*2 1298 hom(:,1,31,0) = sums(:,31) / ngp_2dh_outer(:,0) ! v*2 1299 hom(:,1,32,0) = sums(:,32) / ngp_2dh_outer(:,0) ! w*2 1300 1301 ENDIF 1302 1303 ENDIF 81 ENDIF 82 83 84 ! 85 !-- Initialize arrays for marking those particles/tails to be deleted after the 86 !-- (sub-) timestep 87 particle_mask = .TRUE. 88 deleted_particles = 0 89 90 IF ( use_particle_tails ) THEN 91 tail_mask = .TRUE. 92 ENDIF 93 deleted_tails = 0 94 1304 95 1305 96 ! … … 1317 108 trnp_count_recv_sum = 0 1318 109 110 111 ! 112 !-- Calculate exponential term used in case of particle inertia for each 113 !-- of the particle groups 114 DO m = 1, number_of_particle_groups 115 IF ( particle_groups(m)%density_ratio /= 0.0 ) THEN 116 particle_groups(m)%exp_arg = & 117 4.5 * particle_groups(m)%density_ratio * & 118 molecular_viscosity / ( particle_groups(m)%radius )**2 119 particle_groups(m)%exp_term = EXP( -particle_groups(m)%exp_arg * & 120 dt_3d ) 121 ENDIF 122 ENDDO 123 124 125 ! 126 !-- Particle (droplet) growth by condensation/evaporation and collision 127 IF ( cloud_droplets ) THEN 128 129 ! 130 !-- Reset summation arrays 131 ql_c = 0.0; ql_v = 0.0; ql_vp = 0.0 132 133 ! 134 !-- Droplet growth by condensation / evaporation 135 CALL lpm_droplet_condensation 136 137 ! 138 !-- Particle growth by collision 139 IF ( collision_kernel /= 'none' ) CALL lpm_droplet_collision 140 141 ENDIF 142 143 144 ! 145 !-- If particle advection includes SGS velocity components, calculate the 146 !-- required SGS quantities (i.e. gradients of the TKE, as well as horizontally 147 !-- averaged profiles of the SGS TKE and the resolved-scale velocity variances) 148 IF ( use_sgs_for_particles ) CALL lpm_init_sgs_tke 149 150 1319 151 ! 1320 152 !-- Initialize the variable storing the total time that a particle has advanced … … 1322 154 particles(1:number_of_particles)%dt_sum = 0.0 1323 155 1324 ! 1325 !-- Timestep loop. 156 157 ! 158 !-- Timestep loop for particle advection. 1326 159 !-- This loop has to be repeated until the advection time of every particle 1327 !-- (in the total domain!) has reached the LES timestep (dt_3d) 160 !-- (within the total domain!) has reached the LES timestep (dt_3d). 161 !-- In case of including the SGS velocities, the particle timestep may be 162 !-- smaller than the LES timestep (because of the Lagrangian timescale restric- 163 !-- tion) and particles may require to undergo several particle timesteps, 164 !-- before the LES timestep is reached. Because the number of these particle 165 !-- timesteps to be carried out is unknown at first, these steps are carried 166 !-- out in the following infinite loop with exit condition. 1328 167 DO 1329 168 1330 CALL cpu_log( log_point_s(44), ' advec_part_advec', 'start' )169 CALL cpu_log( log_point_s(44), 'lpm_advec', 'start' ) 1331 170 1332 171 ! … … 1337 176 dt_3d_reached_l = .TRUE. 1338 177 1339 DO n = 1, number_of_particles 1340 ! 1341 !-- Move particles only if the LES timestep has not (approximately) been 1342 !-- reached 1343 IF ( ( dt_3d - particles(n)%dt_sum ) < 1E-8 ) CYCLE 1344 1345 ! 1346 !-- Interpolate u velocity-component, determine left, front, bottom 1347 !-- index of u-array 1348 i = ( particles(n)%x + 0.5 * dx ) * ddx 1349 j = particles(n)%y * ddy 1350 k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz & 1351 + offset_ocean_nzt ! only exact if equidistant 1352 1353 ! 1354 !-- Interpolation of the velocity components in the xy-plane 1355 x = particles(n)%x + ( 0.5 - i ) * dx 1356 y = particles(n)%y - j * dy 1357 aa = x**2 + y**2 1358 bb = ( dx - x )**2 + y**2 1359 cc = x**2 + ( dy - y )**2 1360 dd = ( dx - x )**2 + ( dy - y )**2 1361 gg = aa + bb + cc + dd 1362 1363 u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) & 1364 + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) & 1365 ) / ( 3.0 * gg ) - u_gtrans 1366 IF ( k+1 == nzt+1 ) THEN 1367 u_int = u_int_l 1368 ELSE 1369 u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) & 1370 + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) & 1371 ) / ( 3.0 * gg ) - u_gtrans 1372 u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * & 1373 ( u_int_u - u_int_l ) 1374 ENDIF 1375 1376 ! 1377 !-- Same procedure for interpolation of the v velocity-component (adopt 1378 !-- index k from u velocity-component) 1379 i = particles(n)%x * ddx 1380 j = ( particles(n)%y + 0.5 * dy ) * ddy 1381 1382 x = particles(n)%x - i * dx 1383 y = particles(n)%y + ( 0.5 - j ) * dy 1384 aa = x**2 + y**2 1385 bb = ( dx - x )**2 + y**2 1386 cc = x**2 + ( dy - y )**2 1387 dd = ( dx - x )**2 + ( dy - y )**2 1388 gg = aa + bb + cc + dd 1389 1390 v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) & 1391 + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) & 1392 ) / ( 3.0 * gg ) - v_gtrans 1393 IF ( k+1 == nzt+1 ) THEN 1394 v_int = v_int_l 1395 ELSE 1396 v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) & 1397 + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) & 1398 ) / ( 3.0 * gg ) - v_gtrans 1399 v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * & 1400 ( v_int_u - v_int_l ) 1401 ENDIF 1402 1403 ! 1404 !-- Same procedure for interpolation of the w velocity-component (adopt 1405 !-- index i from v velocity-component) 1406 IF ( vertical_particle_advection(particles(n)%group) ) THEN 1407 j = particles(n)%y * ddy 1408 k = particles(n)%z / dz + offset_ocean_nzt_m1 1409 1410 x = particles(n)%x - i * dx 1411 y = particles(n)%y - j * dy 1412 aa = x**2 + y**2 1413 bb = ( dx - x )**2 + y**2 1414 cc = x**2 + ( dy - y )**2 1415 dd = ( dx - x )**2 + ( dy - y )**2 1416 gg = aa + bb + cc + dd 1417 1418 w_int_l = ( ( gg - aa ) * w(k,j,i) + ( gg - bb ) * w(k,j,i+1) & 1419 + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) & 1420 ) / ( 3.0 * gg ) 1421 IF ( k+1 == nzt+1 ) THEN 1422 w_int = w_int_l 1423 ELSE 1424 w_int_u = ( ( gg-aa ) * w(k+1,j,i) + & 1425 ( gg-bb ) * w(k+1,j,i+1) + & 1426 ( gg-cc ) * w(k+1,j+1,i) + & 1427 ( gg-dd ) * w(k+1,j+1,i+1) & 1428 ) / ( 3.0 * gg ) 1429 w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * & 1430 ( w_int_u - w_int_l ) 1431 ENDIF 1432 ELSE 1433 w_int = 0.0 1434 ENDIF 1435 1436 ! 1437 !-- Interpolate and calculate quantities needed for calculating the SGS 1438 !-- velocities 1439 IF ( use_sgs_for_particles ) THEN 1440 ! 1441 !-- Interpolate TKE 1442 i = particles(n)%x * ddx 1443 j = particles(n)%y * ddy 1444 k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz & 1445 + offset_ocean_nzt ! only exact if eq.dist 1446 1447 IF ( topography == 'flat' ) THEN 1448 1449 x = particles(n)%x - i * dx 1450 y = particles(n)%y - j * dy 1451 aa = x**2 + y**2 1452 bb = ( dx - x )**2 + y**2 1453 cc = x**2 + ( dy - y )**2 1454 dd = ( dx - x )**2 + ( dy - y )**2 1455 gg = aa + bb + cc + dd 1456 1457 e_int_l = ( ( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) & 1458 + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) & 1459 ) / ( 3.0 * gg ) 1460 1461 IF ( k+1 == nzt+1 ) THEN 1462 e_int = e_int_l 1463 ELSE 1464 e_int_u = ( ( gg - aa ) * e(k+1,j,i) + & 1465 ( gg - bb ) * e(k+1,j,i+1) + & 1466 ( gg - cc ) * e(k+1,j+1,i) + & 1467 ( gg - dd ) * e(k+1,j+1,i+1) & 1468 ) / ( 3.0 * gg ) 1469 e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * & 1470 ( e_int_u - e_int_l ) 1471 ENDIF 1472 1473 ! 1474 !-- Interpolate the TKE gradient along x (adopt incides i,j,k and 1475 !-- all position variables from above (TKE)) 1476 de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + & 1477 ( gg - bb ) * de_dx(k,j,i+1) + & 1478 ( gg - cc ) * de_dx(k,j+1,i) + & 1479 ( gg - dd ) * de_dx(k,j+1,i+1) & 1480 ) / ( 3.0 * gg ) 1481 1482 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1483 de_dx_int = de_dx_int_l 1484 ELSE 1485 de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + & 1486 ( gg - bb ) * de_dx(k+1,j,i+1) + & 1487 ( gg - cc ) * de_dx(k+1,j+1,i) + & 1488 ( gg - dd ) * de_dx(k+1,j+1,i+1) & 1489 ) / ( 3.0 * gg ) 1490 de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * & 1491 ( de_dx_int_u - de_dx_int_l ) 1492 ENDIF 1493 1494 ! 1495 !-- Interpolate the TKE gradient along y 1496 de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + & 1497 ( gg - bb ) * de_dy(k,j,i+1) + & 1498 ( gg - cc ) * de_dy(k,j+1,i) + & 1499 ( gg - dd ) * de_dy(k,j+1,i+1) & 1500 ) / ( 3.0 * gg ) 1501 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1502 de_dy_int = de_dy_int_l 1503 ELSE 1504 de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + & 1505 ( gg - bb ) * de_dy(k+1,j,i+1) + & 1506 ( gg - cc ) * de_dy(k+1,j+1,i) + & 1507 ( gg - dd ) * de_dy(k+1,j+1,i+1) & 1508 ) / ( 3.0 * gg ) 1509 de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * & 1510 ( de_dy_int_u - de_dy_int_l ) 1511 ENDIF 1512 1513 ! 1514 !-- Interpolate the TKE gradient along z 1515 IF ( particles(n)%z < 0.5 * dz ) THEN 1516 de_dz_int = 0.0 1517 ELSE 1518 de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + & 1519 ( gg - bb ) * de_dz(k,j,i+1) + & 1520 ( gg - cc ) * de_dz(k,j+1,i) + & 1521 ( gg - dd ) * de_dz(k,j+1,i+1) & 1522 ) / ( 3.0 * gg ) 1523 1524 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1525 de_dz_int = de_dz_int_l 1526 ELSE 1527 de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + & 1528 ( gg - bb ) * de_dz(k+1,j,i+1) + & 1529 ( gg - cc ) * de_dz(k+1,j+1,i) + & 1530 ( gg - dd ) * de_dz(k+1,j+1,i+1) & 1531 ) / ( 3.0 * gg ) 1532 de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * & 1533 ( de_dz_int_u - de_dz_int_l ) 1534 ENDIF 1535 ENDIF 1536 1537 ! 1538 !-- Interpolate the dissipation of TKE 1539 diss_int_l = ( ( gg - aa ) * diss(k,j,i) + & 1540 ( gg - bb ) * diss(k,j,i+1) + & 1541 ( gg - cc ) * diss(k,j+1,i) + & 1542 ( gg - dd ) * diss(k,j+1,i+1) & 1543 ) / ( 3.0 * gg ) 1544 1545 IF ( k+1 == nzt+1 ) THEN 1546 diss_int = diss_int_l 1547 ELSE 1548 diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + & 1549 ( gg - bb ) * diss(k+1,j,i+1) + & 1550 ( gg - cc ) * diss(k+1,j+1,i) + & 1551 ( gg - dd ) * diss(k+1,j+1,i+1) & 1552 ) / ( 3.0 * gg ) 1553 diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * & 1554 ( diss_int_u - diss_int_l ) 1555 ENDIF 1556 1557 ELSE 1558 1559 ! 1560 !-- In case that there are buildings it has to be determined 1561 !-- how many of the gridpoints defining the particle box are 1562 !-- situated within a building 1563 !-- gp_outside_of_building(1): i,j,k 1564 !-- gp_outside_of_building(2): i,j+1,k 1565 !-- gp_outside_of_building(3): i,j,k+1 1566 !-- gp_outside_of_building(4): i,j+1,k+1 1567 !-- gp_outside_of_building(5): i+1,j,k 1568 !-- gp_outside_of_building(6): i+1,j+1,k 1569 !-- gp_outside_of_building(7): i+1,j,k+1 1570 !-- gp_outside_of_building(8): i+1,j+1,k+1 1571 1572 gp_outside_of_building = 0 1573 location = 0.0 1574 num_gp = 0 1575 1576 IF ( k > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN 1577 num_gp = num_gp + 1 1578 gp_outside_of_building(1) = 1 1579 location(num_gp,1) = i * dx 1580 location(num_gp,2) = j * dy 1581 location(num_gp,3) = k * dz - 0.5 * dz 1582 ei(num_gp) = e(k,j,i) 1583 dissi(num_gp) = diss(k,j,i) 1584 de_dxi(num_gp) = de_dx(k,j,i) 1585 de_dyi(num_gp) = de_dy(k,j,i) 1586 de_dzi(num_gp) = de_dz(k,j,i) 1587 ENDIF 1588 1589 IF ( k > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) & 1590 THEN 1591 num_gp = num_gp + 1 1592 gp_outside_of_building(2) = 1 1593 location(num_gp,1) = i * dx 1594 location(num_gp,2) = (j+1) * dy 1595 location(num_gp,3) = k * dz - 0.5 * dz 1596 ei(num_gp) = e(k,j+1,i) 1597 dissi(num_gp) = diss(k,j+1,i) 1598 de_dxi(num_gp) = de_dx(k,j+1,i) 1599 de_dyi(num_gp) = de_dy(k,j+1,i) 1600 de_dzi(num_gp) = de_dz(k,j+1,i) 1601 ENDIF 1602 1603 IF ( k+1 > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN 1604 num_gp = num_gp + 1 1605 gp_outside_of_building(3) = 1 1606 location(num_gp,1) = i * dx 1607 location(num_gp,2) = j * dy 1608 location(num_gp,3) = (k+1) * dz - 0.5 * dz 1609 ei(num_gp) = e(k+1,j,i) 1610 dissi(num_gp) = diss(k+1,j,i) 1611 de_dxi(num_gp) = de_dx(k+1,j,i) 1612 de_dyi(num_gp) = de_dy(k+1,j,i) 1613 de_dzi(num_gp) = de_dz(k+1,j,i) 1614 ENDIF 1615 1616 IF ( k+1 > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) & 1617 THEN 1618 num_gp = num_gp + 1 1619 gp_outside_of_building(4) = 1 1620 location(num_gp,1) = i * dx 1621 location(num_gp,2) = (j+1) * dy 1622 location(num_gp,3) = (k+1) * dz - 0.5 * dz 1623 ei(num_gp) = e(k+1,j+1,i) 1624 dissi(num_gp) = diss(k+1,j+1,i) 1625 de_dxi(num_gp) = de_dx(k+1,j+1,i) 1626 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1627 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1628 ENDIF 1629 1630 IF ( k > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) & 1631 THEN 1632 num_gp = num_gp + 1 1633 gp_outside_of_building(5) = 1 1634 location(num_gp,1) = (i+1) * dx 1635 location(num_gp,2) = j * dy 1636 location(num_gp,3) = k * dz - 0.5 * dz 1637 ei(num_gp) = e(k,j,i+1) 1638 dissi(num_gp) = diss(k,j,i+1) 1639 de_dxi(num_gp) = de_dx(k,j,i+1) 1640 de_dyi(num_gp) = de_dy(k,j,i+1) 1641 de_dzi(num_gp) = de_dz(k,j,i+1) 1642 ENDIF 1643 1644 IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) & 1645 THEN 1646 num_gp = num_gp + 1 1647 gp_outside_of_building(6) = 1 1648 location(num_gp,1) = (i+1) * dx 1649 location(num_gp,2) = (j+1) * dy 1650 location(num_gp,3) = k * dz - 0.5 * dz 1651 ei(num_gp) = e(k,j+1,i+1) 1652 dissi(num_gp) = diss(k,j+1,i+1) 1653 de_dxi(num_gp) = de_dx(k,j+1,i+1) 1654 de_dyi(num_gp) = de_dy(k,j+1,i+1) 1655 de_dzi(num_gp) = de_dz(k,j+1,i+1) 1656 ENDIF 1657 1658 IF ( k+1 > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) & 1659 THEN 1660 num_gp = num_gp + 1 1661 gp_outside_of_building(7) = 1 1662 location(num_gp,1) = (i+1) * dx 1663 location(num_gp,2) = j * dy 1664 location(num_gp,3) = (k+1) * dz - 0.5 * dz 1665 ei(num_gp) = e(k+1,j,i+1) 1666 dissi(num_gp) = diss(k+1,j,i+1) 1667 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1668 de_dyi(num_gp) = de_dy(k+1,j,i+1) 1669 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1670 ENDIF 1671 1672 IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0)& 1673 THEN 1674 num_gp = num_gp + 1 1675 gp_outside_of_building(8) = 1 1676 location(num_gp,1) = (i+1) * dx 1677 location(num_gp,2) = (j+1) * dy 1678 location(num_gp,3) = (k+1) * dz - 0.5 * dz 1679 ei(num_gp) = e(k+1,j+1,i+1) 1680 dissi(num_gp) = diss(k+1,j+1,i+1) 1681 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 1682 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 1683 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 1684 ENDIF 1685 1686 ! 1687 !-- If all gridpoints are situated outside of a building, then the 1688 !-- ordinary interpolation scheme can be used. 1689 IF ( num_gp == 8 ) THEN 1690 1691 x = particles(n)%x - i * dx 1692 y = particles(n)%y - j * dy 1693 aa = x**2 + y**2 1694 bb = ( dx - x )**2 + y**2 1695 cc = x**2 + ( dy - y )**2 1696 dd = ( dx - x )**2 + ( dy - y )**2 1697 gg = aa + bb + cc + dd 1698 1699 e_int_l = (( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) & 1700 + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1)& 1701 ) / ( 3.0 * gg ) 1702 1703 IF ( k+1 == nzt+1 ) THEN 1704 e_int = e_int_l 1705 ELSE 1706 e_int_u = ( ( gg - aa ) * e(k+1,j,i) + & 1707 ( gg - bb ) * e(k+1,j,i+1) + & 1708 ( gg - cc ) * e(k+1,j+1,i) + & 1709 ( gg - dd ) * e(k+1,j+1,i+1) & 1710 ) / ( 3.0 * gg ) 1711 e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * & 1712 ( e_int_u - e_int_l ) 1713 ENDIF 1714 1715 ! 1716 !-- Interpolate the TKE gradient along x (adopt incides i,j,k 1717 !-- and all position variables from above (TKE)) 1718 de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + & 1719 ( gg - bb ) * de_dx(k,j,i+1) + & 1720 ( gg - cc ) * de_dx(k,j+1,i) + & 1721 ( gg - dd ) * de_dx(k,j+1,i+1) & 1722 ) / ( 3.0 * gg ) 1723 1724 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1725 de_dx_int = de_dx_int_l 1726 ELSE 1727 de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + & 1728 ( gg - bb ) * de_dx(k+1,j,i+1) + & 1729 ( gg - cc ) * de_dx(k+1,j+1,i) + & 1730 ( gg - dd ) * de_dx(k+1,j+1,i+1) & 1731 ) / ( 3.0 * gg ) 1732 de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / & 1733 dz * ( de_dx_int_u - de_dx_int_l ) 1734 ENDIF 1735 1736 ! 1737 !-- Interpolate the TKE gradient along y 1738 de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + & 1739 ( gg - bb ) * de_dy(k,j,i+1) + & 1740 ( gg - cc ) * de_dy(k,j+1,i) + & 1741 ( gg - dd ) * de_dy(k,j+1,i+1) & 1742 ) / ( 3.0 * gg ) 1743 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1744 de_dy_int = de_dy_int_l 1745 ELSE 1746 de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + & 1747 ( gg - bb ) * de_dy(k+1,j,i+1) + & 1748 ( gg - cc ) * de_dy(k+1,j+1,i) + & 1749 ( gg - dd ) * de_dy(k+1,j+1,i+1) & 1750 ) / ( 3.0 * gg ) 1751 de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / & 1752 dz * ( de_dy_int_u - de_dy_int_l ) 1753 ENDIF 1754 1755 ! 1756 !-- Interpolate the TKE gradient along z 1757 IF ( particles(n)%z < 0.5 * dz ) THEN 1758 de_dz_int = 0.0 1759 ELSE 1760 de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + & 1761 ( gg - bb ) * de_dz(k,j,i+1) + & 1762 ( gg - cc ) * de_dz(k,j+1,i) + & 1763 ( gg - dd ) * de_dz(k,j+1,i+1) & 1764 ) / ( 3.0 * gg ) 1765 1766 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 1767 de_dz_int = de_dz_int_l 1768 ELSE 1769 de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + & 1770 ( gg - bb ) * de_dz(k+1,j,i+1) + & 1771 ( gg - cc ) * de_dz(k+1,j+1,i) + & 1772 ( gg - dd ) * de_dz(k+1,j+1,i+1) & 1773 ) / ( 3.0 * gg ) 1774 de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) /& 1775 dz * ( de_dz_int_u - de_dz_int_l ) 1776 ENDIF 1777 ENDIF 1778 1779 ! 1780 !-- Interpolate the dissipation of TKE 1781 diss_int_l = ( ( gg - aa ) * diss(k,j,i) + & 1782 ( gg - bb ) * diss(k,j,i+1) + & 1783 ( gg - cc ) * diss(k,j+1,i) + & 1784 ( gg - dd ) * diss(k,j+1,i+1) & 1785 ) / ( 3.0 * gg ) 1786 1787 IF ( k+1 == nzt+1 ) THEN 1788 diss_int = diss_int_l 1789 ELSE 1790 diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + & 1791 ( gg - bb ) * diss(k+1,j,i+1) + & 1792 ( gg - cc ) * diss(k+1,j+1,i) + & 1793 ( gg - dd ) * diss(k+1,j+1,i+1) & 1794 ) / ( 3.0 * gg ) 1795 diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz *& 1796 ( diss_int_u - diss_int_l ) 1797 ENDIF 1798 1799 ELSE 1800 1801 ! 1802 !-- If wall between gridpoint 1 and gridpoint 5, then 1803 !-- Neumann boundary condition has to be applied 1804 IF ( gp_outside_of_building(1) == 1 .AND. & 1805 gp_outside_of_building(5) == 0 ) THEN 1806 num_gp = num_gp + 1 1807 location(num_gp,1) = i * dx + 0.5 * dx 1808 location(num_gp,2) = j * dy 1809 location(num_gp,3) = k * dz - 0.5 * dz 1810 ei(num_gp) = e(k,j,i) 1811 dissi(num_gp) = diss(k,j,i) 1812 de_dxi(num_gp) = 0.0 1813 de_dyi(num_gp) = de_dy(k,j,i) 1814 de_dzi(num_gp) = de_dz(k,j,i) 1815 ENDIF 1816 1817 IF ( gp_outside_of_building(5) == 1 .AND. & 1818 gp_outside_of_building(1) == 0 ) THEN 1819 num_gp = num_gp + 1 1820 location(num_gp,1) = i * dx + 0.5 * dx 1821 location(num_gp,2) = j * dy 1822 location(num_gp,3) = k * dz - 0.5 * dz 1823 ei(num_gp) = e(k,j,i+1) 1824 dissi(num_gp) = diss(k,j,i+1) 1825 de_dxi(num_gp) = 0.0 1826 de_dyi(num_gp) = de_dy(k,j,i+1) 1827 de_dzi(num_gp) = de_dz(k,j,i+1) 1828 ENDIF 1829 1830 ! 1831 !-- If wall between gridpoint 5 and gridpoint 6, then 1832 !-- then Neumann boundary condition has to be applied 1833 IF ( gp_outside_of_building(5) == 1 .AND. & 1834 gp_outside_of_building(6) == 0 ) THEN 1835 num_gp = num_gp + 1 1836 location(num_gp,1) = (i+1) * dx 1837 location(num_gp,2) = j * dy + 0.5 * dy 1838 location(num_gp,3) = k * dz - 0.5 * dz 1839 ei(num_gp) = e(k,j,i+1) 1840 dissi(num_gp) = diss(k,j,i+1) 1841 de_dxi(num_gp) = de_dx(k,j,i+1) 1842 de_dyi(num_gp) = 0.0 1843 de_dzi(num_gp) = de_dz(k,j,i+1) 1844 ENDIF 1845 1846 IF ( gp_outside_of_building(6) == 1 .AND. & 1847 gp_outside_of_building(5) == 0 ) THEN 1848 num_gp = num_gp + 1 1849 location(num_gp,1) = (i+1) * dx 1850 location(num_gp,2) = j * dy + 0.5 * dy 1851 location(num_gp,3) = k * dz - 0.5 * dz 1852 ei(num_gp) = e(k,j+1,i+1) 1853 dissi(num_gp) = diss(k,j+1,i+1) 1854 de_dxi(num_gp) = de_dx(k,j+1,i+1) 1855 de_dyi(num_gp) = 0.0 1856 de_dzi(num_gp) = de_dz(k,j+1,i+1) 1857 ENDIF 1858 1859 ! 1860 !-- If wall between gridpoint 2 and gridpoint 6, then 1861 !-- Neumann boundary condition has to be applied 1862 IF ( gp_outside_of_building(2) == 1 .AND. & 1863 gp_outside_of_building(6) == 0 ) THEN 1864 num_gp = num_gp + 1 1865 location(num_gp,1) = i * dx + 0.5 * dx 1866 location(num_gp,2) = (j+1) * dy 1867 location(num_gp,3) = k * dz - 0.5 * dz 1868 ei(num_gp) = e(k,j+1,i) 1869 dissi(num_gp) = diss(k,j+1,i) 1870 de_dxi(num_gp) = 0.0 1871 de_dyi(num_gp) = de_dy(k,j+1,i) 1872 de_dzi(num_gp) = de_dz(k,j+1,i) 1873 ENDIF 1874 1875 IF ( gp_outside_of_building(6) == 1 .AND. & 1876 gp_outside_of_building(2) == 0 ) THEN 1877 num_gp = num_gp + 1 1878 location(num_gp,1) = i * dx + 0.5 * dx 1879 location(num_gp,2) = (j+1) * dy 1880 location(num_gp,3) = k * dz - 0.5 * dz 1881 ei(num_gp) = e(k,j+1,i+1) 1882 dissi(num_gp) = diss(k,j+1,i+1) 1883 de_dxi(num_gp) = 0.0 1884 de_dyi(num_gp) = de_dy(k,j+1,i+1) 1885 de_dzi(num_gp) = de_dz(k,j+1,i+1) 1886 ENDIF 1887 1888 ! 1889 !-- If wall between gridpoint 1 and gridpoint 2, then 1890 !-- Neumann boundary condition has to be applied 1891 IF ( gp_outside_of_building(1) == 1 .AND. & 1892 gp_outside_of_building(2) == 0 ) THEN 1893 num_gp = num_gp + 1 1894 location(num_gp,1) = i * dx 1895 location(num_gp,2) = j * dy + 0.5 * dy 1896 location(num_gp,3) = k * dz - 0.5 * dz 1897 ei(num_gp) = e(k,j,i) 1898 dissi(num_gp) = diss(k,j,i) 1899 de_dxi(num_gp) = de_dx(k,j,i) 1900 de_dyi(num_gp) = 0.0 1901 de_dzi(num_gp) = de_dz(k,j,i) 1902 ENDIF 1903 1904 IF ( gp_outside_of_building(2) == 1 .AND. & 1905 gp_outside_of_building(1) == 0 ) THEN 1906 num_gp = num_gp + 1 1907 location(num_gp,1) = i * dx 1908 location(num_gp,2) = j * dy + 0.5 * dy 1909 location(num_gp,3) = k * dz - 0.5 * dz 1910 ei(num_gp) = e(k,j+1,i) 1911 dissi(num_gp) = diss(k,j+1,i) 1912 de_dxi(num_gp) = de_dx(k,j+1,i) 1913 de_dyi(num_gp) = 0.0 1914 de_dzi(num_gp) = de_dz(k,j+1,i) 1915 ENDIF 1916 1917 ! 1918 !-- If wall between gridpoint 3 and gridpoint 7, then 1919 !-- Neumann boundary condition has to be applied 1920 IF ( gp_outside_of_building(3) == 1 .AND. & 1921 gp_outside_of_building(7) == 0 ) THEN 1922 num_gp = num_gp + 1 1923 location(num_gp,1) = i * dx + 0.5 * dx 1924 location(num_gp,2) = j * dy 1925 location(num_gp,3) = k * dz + 0.5 * dz 1926 ei(num_gp) = e(k+1,j,i) 1927 dissi(num_gp) = diss(k+1,j,i) 1928 de_dxi(num_gp) = 0.0 1929 de_dyi(num_gp) = de_dy(k+1,j,i) 1930 de_dzi(num_gp) = de_dz(k+1,j,i) 1931 ENDIF 1932 1933 IF ( gp_outside_of_building(7) == 1 .AND. & 1934 gp_outside_of_building(3) == 0 ) THEN 1935 num_gp = num_gp + 1 1936 location(num_gp,1) = i * dx + 0.5 * dx 1937 location(num_gp,2) = j * dy 1938 location(num_gp,3) = k * dz + 0.5 * dz 1939 ei(num_gp) = e(k+1,j,i+1) 1940 dissi(num_gp) = diss(k+1,j,i+1) 1941 de_dxi(num_gp) = 0.0 1942 de_dyi(num_gp) = de_dy(k+1,j,i+1) 1943 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1944 ENDIF 1945 1946 ! 1947 !-- If wall between gridpoint 7 and gridpoint 8, then 1948 !-- Neumann boundary condition has to be applied 1949 IF ( gp_outside_of_building(7) == 1 .AND. & 1950 gp_outside_of_building(8) == 0 ) THEN 1951 num_gp = num_gp + 1 1952 location(num_gp,1) = (i+1) * dx 1953 location(num_gp,2) = j * dy + 0.5 * dy 1954 location(num_gp,3) = k * dz + 0.5 * dz 1955 ei(num_gp) = e(k+1,j,i+1) 1956 dissi(num_gp) = diss(k+1,j,i+1) 1957 de_dxi(num_gp) = de_dx(k+1,j,i+1) 1958 de_dyi(num_gp) = 0.0 1959 de_dzi(num_gp) = de_dz(k+1,j,i+1) 1960 ENDIF 1961 1962 IF ( gp_outside_of_building(8) == 1 .AND. & 1963 gp_outside_of_building(7) == 0 ) THEN 1964 num_gp = num_gp + 1 1965 location(num_gp,1) = (i+1) * dx 1966 location(num_gp,2) = j * dy + 0.5 * dy 1967 location(num_gp,3) = k * dz + 0.5 * dz 1968 ei(num_gp) = e(k+1,j+1,i+1) 1969 dissi(num_gp) = diss(k+1,j+1,i+1) 1970 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 1971 de_dyi(num_gp) = 0.0 1972 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 1973 ENDIF 1974 1975 ! 1976 !-- If wall between gridpoint 4 and gridpoint 8, then 1977 !-- Neumann boundary condition has to be applied 1978 IF ( gp_outside_of_building(4) == 1 .AND. & 1979 gp_outside_of_building(8) == 0 ) THEN 1980 num_gp = num_gp + 1 1981 location(num_gp,1) = i * dx + 0.5 * dx 1982 location(num_gp,2) = (j+1) * dy 1983 location(num_gp,3) = k * dz + 0.5 * dz 1984 ei(num_gp) = e(k+1,j+1,i) 1985 dissi(num_gp) = diss(k+1,j+1,i) 1986 de_dxi(num_gp) = 0.0 1987 de_dyi(num_gp) = de_dy(k+1,j+1,i) 1988 de_dzi(num_gp) = de_dz(k+1,j+1,i) 1989 ENDIF 1990 1991 IF ( gp_outside_of_building(8) == 1 .AND. & 1992 gp_outside_of_building(4) == 0 ) THEN 1993 num_gp = num_gp + 1 1994 location(num_gp,1) = i * dx + 0.5 * dx 1995 location(num_gp,2) = (j+1) * dy 1996 location(num_gp,3) = k * dz + 0.5 * dz 1997 ei(num_gp) = e(k+1,j+1,i+1) 1998 dissi(num_gp) = diss(k+1,j+1,i+1) 1999 de_dxi(num_gp) = 0.0 2000 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 2001 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 2002 ENDIF 2003 2004 ! 2005 !-- If wall between gridpoint 3 and gridpoint 4, then 2006 !-- Neumann boundary condition has to be applied 2007 IF ( gp_outside_of_building(3) == 1 .AND. & 2008 gp_outside_of_building(4) == 0 ) THEN 2009 num_gp = num_gp + 1 2010 location(num_gp,1) = i * dx 2011 location(num_gp,2) = j * dy + 0.5 * dy 2012 location(num_gp,3) = k * dz + 0.5 * dz 2013 ei(num_gp) = e(k+1,j,i) 2014 dissi(num_gp) = diss(k+1,j,i) 2015 de_dxi(num_gp) = de_dx(k+1,j,i) 2016 de_dyi(num_gp) = 0.0 2017 de_dzi(num_gp) = de_dz(k+1,j,i) 2018 ENDIF 2019 2020 IF ( gp_outside_of_building(4) == 1 .AND. & 2021 gp_outside_of_building(3) == 0 ) THEN 2022 num_gp = num_gp + 1 2023 location(num_gp,1) = i * dx 2024 location(num_gp,2) = j * dy + 0.5 * dy 2025 location(num_gp,3) = k * dz + 0.5 * dz 2026 ei(num_gp) = e(k+1,j+1,i) 2027 dissi(num_gp) = diss(k+1,j+1,i) 2028 de_dxi(num_gp) = de_dx(k+1,j+1,i) 2029 de_dyi(num_gp) = 0.0 2030 de_dzi(num_gp) = de_dz(k+1,j+1,i) 2031 ENDIF 2032 2033 ! 2034 !-- If wall between gridpoint 1 and gridpoint 3, then 2035 !-- Neumann boundary condition has to be applied 2036 !-- (only one case as only building beneath is possible) 2037 IF ( gp_outside_of_building(1) == 0 .AND. & 2038 gp_outside_of_building(3) == 1 ) THEN 2039 num_gp = num_gp + 1 2040 location(num_gp,1) = i * dx 2041 location(num_gp,2) = j * dy 2042 location(num_gp,3) = k * dz 2043 ei(num_gp) = e(k+1,j,i) 2044 dissi(num_gp) = diss(k+1,j,i) 2045 de_dxi(num_gp) = de_dx(k+1,j,i) 2046 de_dyi(num_gp) = de_dy(k+1,j,i) 2047 de_dzi(num_gp) = 0.0 2048 ENDIF 2049 2050 ! 2051 !-- If wall between gridpoint 5 and gridpoint 7, then 2052 !-- Neumann boundary condition has to be applied 2053 !-- (only one case as only building beneath is possible) 2054 IF ( gp_outside_of_building(5) == 0 .AND. & 2055 gp_outside_of_building(7) == 1 ) THEN 2056 num_gp = num_gp + 1 2057 location(num_gp,1) = (i+1) * dx 2058 location(num_gp,2) = j * dy 2059 location(num_gp,3) = k * dz 2060 ei(num_gp) = e(k+1,j,i+1) 2061 dissi(num_gp) = diss(k+1,j,i+1) 2062 de_dxi(num_gp) = de_dx(k+1,j,i+1) 2063 de_dyi(num_gp) = de_dy(k+1,j,i+1) 2064 de_dzi(num_gp) = 0.0 2065 ENDIF 2066 2067 ! 2068 !-- If wall between gridpoint 2 and gridpoint 4, then 2069 !-- Neumann boundary condition has to be applied 2070 !-- (only one case as only building beneath is possible) 2071 IF ( gp_outside_of_building(2) == 0 .AND. & 2072 gp_outside_of_building(4) == 1 ) THEN 2073 num_gp = num_gp + 1 2074 location(num_gp,1) = i * dx 2075 location(num_gp,2) = (j+1) * dy 2076 location(num_gp,3) = k * dz 2077 ei(num_gp) = e(k+1,j+1,i) 2078 dissi(num_gp) = diss(k+1,j+1,i) 2079 de_dxi(num_gp) = de_dx(k+1,j+1,i) 2080 de_dyi(num_gp) = de_dy(k+1,j+1,i) 2081 de_dzi(num_gp) = 0.0 2082 ENDIF 2083 2084 ! 2085 !-- If wall between gridpoint 6 and gridpoint 8, then 2086 !-- Neumann boundary condition has to be applied 2087 !-- (only one case as only building beneath is possible) 2088 IF ( gp_outside_of_building(6) == 0 .AND. & 2089 gp_outside_of_building(8) == 1 ) THEN 2090 num_gp = num_gp + 1 2091 location(num_gp,1) = (i+1) * dx 2092 location(num_gp,2) = (j+1) * dy 2093 location(num_gp,3) = k * dz 2094 ei(num_gp) = e(k+1,j+1,i+1) 2095 dissi(num_gp) = diss(k+1,j+1,i+1) 2096 de_dxi(num_gp) = de_dx(k+1,j+1,i+1) 2097 de_dyi(num_gp) = de_dy(k+1,j+1,i+1) 2098 de_dzi(num_gp) = 0.0 2099 ENDIF 2100 2101 ! 2102 !-- Carry out the interpolation 2103 IF ( num_gp == 1 ) THEN 2104 ! 2105 !-- If only one of the gridpoints is situated outside of the 2106 !-- building, it follows that the values at the particle 2107 !-- location are the same as the gridpoint values 2108 e_int = ei(num_gp) 2109 diss_int = dissi(num_gp) 2110 de_dx_int = de_dxi(num_gp) 2111 de_dy_int = de_dyi(num_gp) 2112 de_dz_int = de_dzi(num_gp) 2113 ELSE IF ( num_gp > 1 ) THEN 2114 2115 d_sum = 0.0 2116 ! 2117 !-- Evaluation of the distances between the gridpoints 2118 !-- contributing to the interpolated values, and the particle 2119 !-- location 2120 DO agp = 1, num_gp 2121 d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2 & 2122 + ( particles(n)%y-location(agp,2) )**2 & 2123 + ( particles(n)%z-location(agp,3) )**2 2124 d_sum = d_sum + d_gp_pl(agp) 2125 ENDDO 2126 2127 ! 2128 !-- Finally the interpolation can be carried out 2129 e_int = 0.0 2130 diss_int = 0.0 2131 de_dx_int = 0.0 2132 de_dy_int = 0.0 2133 de_dz_int = 0.0 2134 DO agp = 1, num_gp 2135 e_int = e_int + ( d_sum - d_gp_pl(agp) ) * & 2136 ei(agp) / ( (num_gp-1) * d_sum ) 2137 diss_int = diss_int + ( d_sum - d_gp_pl(agp) ) * & 2138 dissi(agp) / ( (num_gp-1) * d_sum ) 2139 de_dx_int = de_dx_int + ( d_sum - d_gp_pl(agp) ) * & 2140 de_dxi(agp) / ( (num_gp-1) * d_sum ) 2141 de_dy_int = de_dy_int + ( d_sum - d_gp_pl(agp) ) * & 2142 de_dyi(agp) / ( (num_gp-1) * d_sum ) 2143 de_dz_int = de_dz_int + ( d_sum - d_gp_pl(agp) ) * & 2144 de_dzi(agp) / ( (num_gp-1) * d_sum ) 2145 ENDDO 2146 2147 ENDIF 2148 2149 ENDIF 2150 2151 ENDIF 2152 2153 ! 2154 !-- Vertically interpolate the horizontally averaged SGS TKE and 2155 !-- resolved-scale velocity variances and use the interpolated values 2156 !-- to calculate the coefficient fs, which is a measure of the ratio 2157 !-- of the subgrid-scale turbulent kinetic energy to the total amount 2158 !-- of turbulent kinetic energy. 2159 IF ( k == 0 ) THEN 2160 e_mean_int = hom(0,1,8,0) 2161 ELSE 2162 e_mean_int = hom(k,1,8,0) + & 2163 ( hom(k+1,1,8,0) - hom(k,1,8,0) ) / & 2164 ( zu(k+1) - zu(k) ) * & 2165 ( particles(n)%z - zu(k) ) 2166 ENDIF 2167 2168 kw = particles(n)%z / dz 2169 2170 IF ( k == 0 ) THEN 2171 aa = hom(k+1,1,30,0) * ( particles(n)%z / & 2172 ( 0.5 * ( zu(k+1) - zu(k) ) ) ) 2173 bb = hom(k+1,1,31,0) * ( particles(n)%z / & 2174 ( 0.5 * ( zu(k+1) - zu(k) ) ) ) 2175 cc = hom(kw+1,1,32,0) * ( particles(n)%z / & 2176 ( 1.0 * ( zw(kw+1) - zw(kw) ) ) ) 2177 ELSE 2178 aa = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) * & 2179 ( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) ) 2180 bb = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) * & 2181 ( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) ) 2182 cc = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) *& 2183 ( ( particles(n)%z - zw(kw) ) / ( zw(kw+1)-zw(kw) ) ) 2184 ENDIF 2185 2186 vv_int = ( 1.0 / 3.0 ) * ( aa + bb + cc ) 2187 2188 fs_int = ( 2.0 / 3.0 ) * e_mean_int / & 2189 ( vv_int + ( 2.0 / 3.0 ) * e_mean_int ) 2190 2191 ! 2192 !-- Calculate the Lagrangian timescale according to the suggestion of 2193 !-- Weil et al. (2004). 2194 lagr_timescale = ( 4.0 * e_int ) / & 2195 ( 3.0 * fs_int * c_0 * diss_int ) 2196 2197 ! 2198 !-- Calculate the next particle timestep. dt_gap is the time needed to 2199 !-- complete the current LES timestep. 2200 dt_gap = dt_3d - particles(n)%dt_sum 2201 dt_particle = MIN( dt_3d, 0.025 * lagr_timescale, dt_gap ) 2202 2203 ! 2204 !-- The particle timestep should not be too small in order to prevent 2205 !-- the number of particle timesteps of getting too large 2206 IF ( dt_particle < dt_min_part .AND. dt_min_part < dt_gap ) & 2207 THEN 2208 dt_particle = dt_min_part 2209 ENDIF 2210 2211 ! 2212 !-- Calculate the SGS velocity components 2213 IF ( particles(n)%age == 0.0 ) THEN 2214 ! 2215 !-- For new particles the SGS components are derived from the SGS 2216 !-- TKE. Limit the Gaussian random number to the interval 2217 !-- [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities 2218 !-- from becoming unrealistically large. 2219 particles(n)%rvar1 = SQRT( 2.0 * sgs_wfu_part * e_int ) * & 2220 ( random_gauss( iran_part, 5.0 ) & 2221 - 1.0 ) 2222 particles(n)%rvar2 = SQRT( 2.0 * sgs_wfv_part * e_int ) * & 2223 ( random_gauss( iran_part, 5.0 ) & 2224 - 1.0 ) 2225 particles(n)%rvar3 = SQRT( 2.0 * sgs_wfw_part * e_int ) * & 2226 ( random_gauss( iran_part, 5.0 ) & 2227 - 1.0 ) 2228 2229 ELSE 2230 2231 ! 2232 !-- Restriction of the size of the new timestep: compared to the 2233 !-- previous timestep the increase must not exceed 200% 2234 2235 dt_particle_m = particles(n)%age - particles(n)%age_m 2236 IF ( dt_particle > 2.0 * dt_particle_m ) THEN 2237 dt_particle = 2.0 * dt_particle_m 2238 ENDIF 2239 2240 ! 2241 !-- For old particles the SGS components are correlated with the 2242 !-- values from the previous timestep. Random numbers have also to 2243 !-- be limited (see above). 2244 !-- As negative values for the subgrid TKE are not allowed, the 2245 !-- change of the subgrid TKE with time cannot be smaller than 2246 !-- -e_int/dt_particle. This value is used as a lower boundary 2247 !-- value for the change of TKE 2248 2249 de_dt_min = - e_int / dt_particle 2250 2251 de_dt = ( e_int - particles(n)%e_m ) / dt_particle_m 2252 2253 IF ( de_dt < de_dt_min ) THEN 2254 de_dt = de_dt_min 2255 ENDIF 2256 2257 particles(n)%rvar1 = particles(n)%rvar1 - & 2258 fs_int * c_0 * diss_int * particles(n)%rvar1 * & 2259 dt_particle / ( 4.0 * sgs_wfu_part * e_int ) + & 2260 ( 2.0 * sgs_wfu_part * de_dt * & 2261 particles(n)%rvar1 / & 2262 ( 2.0 * sgs_wfu_part * e_int ) + de_dx_int & 2263 ) * dt_particle / 2.0 + & 2264 SQRT( fs_int * c_0 * diss_int ) * & 2265 ( random_gauss( iran_part, 5.0 ) - 1.0 ) * & 2266 SQRT( dt_particle ) 2267 2268 particles(n)%rvar2 = particles(n)%rvar2 - & 2269 fs_int * c_0 * diss_int * particles(n)%rvar2 * & 2270 dt_particle / ( 4.0 * sgs_wfv_part * e_int ) + & 2271 ( 2.0 * sgs_wfv_part * de_dt * & 2272 particles(n)%rvar2 / & 2273 ( 2.0 * sgs_wfv_part * e_int ) + de_dy_int & 2274 ) * dt_particle / 2.0 + & 2275 SQRT( fs_int * c_0 * diss_int ) * & 2276 ( random_gauss( iran_part, 5.0 ) - 1.0 ) * & 2277 SQRT( dt_particle ) 2278 2279 particles(n)%rvar3 = particles(n)%rvar3 - & 2280 fs_int * c_0 * diss_int * particles(n)%rvar3 * & 2281 dt_particle / ( 4.0 * sgs_wfw_part * e_int ) + & 2282 ( 2.0 * sgs_wfw_part * de_dt * & 2283 particles(n)%rvar3 / & 2284 ( 2.0 * sgs_wfw_part * e_int ) + de_dz_int & 2285 ) * dt_particle / 2.0 + & 2286 SQRT( fs_int * c_0 * diss_int ) * & 2287 ( random_gauss( iran_part, 5.0 ) - 1.0 ) * & 2288 SQRT( dt_particle ) 2289 2290 ENDIF 2291 2292 u_int = u_int + particles(n)%rvar1 2293 v_int = v_int + particles(n)%rvar2 2294 w_int = w_int + particles(n)%rvar3 2295 2296 ! 2297 !-- Store the SGS TKE of the current timelevel which is needed for 2298 !-- for calculating the SGS particle velocities at the next timestep 2299 particles(n)%e_m = e_int 2300 2301 ELSE 2302 ! 2303 !-- If no SGS velocities are used, only the particle timestep has to 2304 !-- be set 2305 dt_particle = dt_3d 2306 2307 ENDIF 2308 2309 ! 2310 !-- Remember the old age of the particle ( needed to prevent that a 2311 !-- particle crosses several PEs during one timestep and for the 2312 !-- evaluation of the subgrid particle velocity fluctuations ) 2313 particles(n)%age_m = particles(n)%age 2314 2315 2316 ! 2317 !-- Particle advection 2318 IF ( particle_groups(particles(n)%group)%density_ratio == 0.0 ) THEN 2319 ! 2320 !-- Pure passive transport (without particle inertia) 2321 particles(n)%x = particles(n)%x + u_int * dt_particle 2322 particles(n)%y = particles(n)%y + v_int * dt_particle 2323 particles(n)%z = particles(n)%z + w_int * dt_particle 2324 2325 particles(n)%speed_x = u_int 2326 particles(n)%speed_y = v_int 2327 particles(n)%speed_z = w_int 2328 2329 ELSE 2330 ! 2331 !-- Transport of particles with inertia 2332 particles(n)%x = particles(n)%x + particles(n)%speed_x * & 2333 dt_particle 2334 particles(n)%y = particles(n)%y + particles(n)%speed_y * & 2335 dt_particle 2336 particles(n)%z = particles(n)%z + particles(n)%speed_z * & 2337 dt_particle 2338 2339 ! 2340 !-- Update of the particle velocity 2341 dens_ratio = particle_groups(particles(n)%group)%density_ratio 2342 IF ( cloud_droplets ) THEN 2343 exp_arg = 4.5 * dens_ratio * molecular_viscosity / & 2344 ( particles(n)%radius )**2 * & 2345 ( 1.0 + 0.15 * ( 2.0 * particles(n)%radius * & 2346 SQRT( ( u_int - particles(n)%speed_x )**2 + & 2347 ( v_int - particles(n)%speed_y )**2 + & 2348 ( w_int - particles(n)%speed_z )**2 ) / & 2349 molecular_viscosity )**0.687 & 2350 ) 2351 exp_term = EXP( -exp_arg * dt_particle ) 2352 ELSEIF ( use_sgs_for_particles ) THEN 2353 exp_arg = particle_groups(particles(n)%group)%exp_arg 2354 exp_term = EXP( -exp_arg * dt_particle ) 2355 ELSE 2356 exp_arg = particle_groups(particles(n)%group)%exp_arg 2357 exp_term = particle_groups(particles(n)%group)%exp_term 2358 ENDIF 2359 particles(n)%speed_x = particles(n)%speed_x * exp_term + & 2360 u_int * ( 1.0 - exp_term ) 2361 particles(n)%speed_y = particles(n)%speed_y * exp_term + & 2362 v_int * ( 1.0 - exp_term ) 2363 particles(n)%speed_z = particles(n)%speed_z * exp_term + & 2364 ( w_int - ( 1.0 - dens_ratio ) * g / exp_arg ) & 2365 * ( 1.0 - exp_term ) 2366 ENDIF 2367 2368 ! 2369 !-- Increment the particle age and the total time that the particle 2370 !-- has advanced within the particle timestep procedure 2371 particles(n)%age = particles(n)%age + dt_particle 2372 particles(n)%dt_sum = particles(n)%dt_sum + dt_particle 2373 2374 ! 2375 !-- Check whether there is still a particle that has not yet completed 2376 !-- the total LES timestep 2377 IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8 ) THEN 2378 dt_3d_reached_l = .FALSE. 2379 ENDIF 2380 2381 ENDDO ! advection loop 178 ! 179 !-- Particle advection 180 CALL lpm_advec 2382 181 2383 182 ! 2384 183 !-- Particle reflection from walls 2385 CALL particle_boundary_conds184 CALL lpm_boundary_conds( 'walls' ) 2386 185 2387 186 ! 2388 187 !-- User-defined actions after the calculation of the new particle position 2389 CALL user_ advec_particles188 CALL user_lpm_advec 2390 189 2391 190 ! … … 2400 199 #endif 2401 200 2402 CALL cpu_log( log_point_s(44), ' advec_part_advec', 'stop' )201 CALL cpu_log( log_point_s(44), 'lpm_advec', 'stop' ) 2403 202 2404 203 ! 2405 204 !-- Increment time since last release 2406 205 IF ( dt_3d_reached ) time_prel = time_prel + dt_3d 2407 2408 ! WRITE ( 9, * ) '*** advec_particles: ##0.4'2409 ! CALL local_flush( 9 )2410 ! nd = 02411 ! DO n = 1, number_of_particles2412 ! IF ( .NOT. particle_mask(n) ) nd = nd + 12413 ! ENDDO2414 ! IF ( nd /= deleted_particles ) THEN2415 ! WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles2416 ! CALL local_flush( 9 )2417 ! CALL MPI_ABORT( comm2d, 9999, ierr )2418 ! ENDIF2419 206 2420 207 ! … … 2423 210 dt_3d_reached ) THEN 2424 211 2425 ! 2426 !-- Check, if particle storage must be extended 2427 IF ( number_of_particles + number_of_initial_particles > & 2428 maximum_number_of_particles ) THEN 2429 IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN 2430 message_string = 'maximum_number_of_particles ' // & 2431 'needs to be increased ' // & 2432 '&but this is not allowed with ' // & 2433 'netcdf_data_format < 3' 2434 CALL message( 'advec_particles', 'PA0146', 2, 2, -1, 6, 1 ) 2435 ELSE 2436 ! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory dt_prel' 2437 ! CALL local_flush( 9 ) 2438 CALL allocate_prt_memory( number_of_initial_particles ) 2439 ! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory dt_prel' 2440 ! CALL local_flush( 9 ) 2441 ENDIF 2442 ENDIF 2443 2444 ! 2445 !-- Check, if tail storage must be extended 2446 IF ( use_particle_tails ) THEN 2447 IF ( number_of_tails + number_of_initial_tails > & 2448 maximum_number_of_tails ) THEN 2449 IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN 2450 message_string = 'maximum_number_of_tails ' // & 2451 'needs to be increased ' // & 2452 '&but this is not allowed wi' // & 2453 'th netcdf_data_format < 3' 2454 CALL message( 'advec_particles', 'PA0147', 2, 2, -1, 6, 1 ) 2455 ELSE 2456 ! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory dt_prel' 2457 ! CALL local_flush( 9 ) 2458 CALL allocate_tail_memory( number_of_initial_tails ) 2459 ! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory dt_prel' 2460 ! CALL local_flush( 9 ) 2461 ENDIF 2462 ENDIF 2463 ENDIF 212 CALL lpm_release_set 2464 213 2465 214 ! … … 2467 216 !-- restart runs. 2468 217 time_prel = MOD( time_prel, MAX( dt_prel, dt_3d ) ) 2469 IF ( number_of_initial_particles /= 0 ) THEN2470 is = number_of_particles+12471 ie = number_of_particles+number_of_initial_particles2472 particles(is:ie) = initial_particles(1:number_of_initial_particles)2473 !2474 !-- Add random fluctuation to particle positions. Particles should2475 !-- remain in the subdomain.2476 IF ( random_start_position ) THEN2477 DO n = is, ie2478 IF ( psl(particles(n)%group) /= psr(particles(n)%group) ) &2479 THEN2480 particles(n)%x = particles(n)%x + &2481 ( random_function( iran_part ) - 0.5 ) *&2482 pdx(particles(n)%group)2483 IF ( particles(n)%x <= ( nxl - 0.5 ) * dx ) THEN2484 particles(n)%x = ( nxl - 0.4999999999 ) * dx2485 ELSEIF ( particles(n)%x >= ( nxr + 0.5 ) * dx ) THEN2486 particles(n)%x = ( nxr + 0.4999999999 ) * dx2487 ENDIF2488 ENDIF2489 IF ( pss(particles(n)%group) /= psn(particles(n)%group) ) &2490 THEN2491 particles(n)%y = particles(n)%y + &2492 ( random_function( iran_part ) - 0.5 ) *&2493 pdy(particles(n)%group)2494 IF ( particles(n)%y <= ( nys - 0.5 ) * dy ) THEN2495 particles(n)%y = ( nys - 0.4999999999 ) * dy2496 ELSEIF ( particles(n)%y >= ( nyn + 0.5 ) * dy ) THEN2497 particles(n)%y = ( nyn + 0.4999999999 ) * dy2498 ENDIF2499 ENDIF2500 IF ( psb(particles(n)%group) /= pst(particles(n)%group) ) &2501 THEN2502 particles(n)%z = particles(n)%z + &2503 ( random_function( iran_part ) - 0.5 ) *&2504 pdz(particles(n)%group)2505 ENDIF2506 ENDDO2507 ENDIF2508 2509 !2510 !-- Set the beginning of the new particle tails and their age2511 IF ( use_particle_tails ) THEN2512 DO n = is, ie2513 !2514 !-- New particles which should have a tail, already have got a2515 !-- provisional tail id unequal zero (see init_particles)2516 IF ( particles(n)%tail_id /= 0 ) THEN2517 number_of_tails = number_of_tails + 12518 nn = number_of_tails2519 particles(n)%tail_id = nn ! set the final tail id2520 particle_tail_coordinates(1,1,nn) = particles(n)%x2521 particle_tail_coordinates(1,2,nn) = particles(n)%y2522 particle_tail_coordinates(1,3,nn) = particles(n)%z2523 particle_tail_coordinates(1,4,nn) = particles(n)%class2524 particles(n)%tailpoints = 12525 IF ( minimum_tailpoint_distance /= 0.0 ) THEN2526 particle_tail_coordinates(2,1,nn) = particles(n)%x2527 particle_tail_coordinates(2,2,nn) = particles(n)%y2528 particle_tail_coordinates(2,3,nn) = particles(n)%z2529 particle_tail_coordinates(2,4,nn) = particles(n)%class2530 particle_tail_coordinates(1:2,5,nn) = 0.02531 particles(n)%tailpoints = 22532 ENDIF2533 ENDIF2534 ENDDO2535 ENDIF2536 ! WRITE ( 9, * ) '*** advec_particles: after setting the beginning of new tails'2537 ! CALL local_flush( 9 )2538 2539 number_of_particles = number_of_particles + &2540 number_of_initial_particles2541 ENDIF2542 218 2543 219 ENDIF 2544 220 2545 ! WRITE ( 9, * ) '*** advec_particles: ##0.5' 2546 ! CALL local_flush( 9 ) 2547 ! nd = 0 2548 ! DO n = 1, number_of_particles 2549 ! IF ( .NOT. particle_mask(n) ) nd = nd + 1 2550 ! ENDDO 2551 ! IF ( nd /= deleted_particles ) THEN 2552 ! WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles 2553 ! CALL local_flush( 9 ) 2554 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 2555 ! ENDIF 2556 2557 ! IF ( number_of_particles /= number_of_tails ) THEN 2558 ! WRITE (9,*) '--- advec_particles: #2' 2559 ! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails 2560 ! CALL local_flush( 9 ) 2561 ! ENDIF 2562 ! DO n = 1, number_of_particles 2563 ! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) & 2564 ! THEN 2565 ! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')' 2566 ! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')' 2567 ! CALL local_flush( 9 ) 2568 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 2569 ! ENDIF 2570 ! ENDDO 2571 2572 #if defined( __parallel ) 2573 ! 2574 !-- As soon as one particle has moved beyond the boundary of the domain, it 2575 !-- is included in the relevant transfer arrays and marked for subsequent 2576 !-- deletion on this PE. 2577 !-- First run for crossings in x direction. Find out first the number of 2578 !-- particles to be transferred and allocate temporary arrays needed to store 2579 !-- them. 2580 !-- For a one-dimensional decomposition along y, no transfer is necessary, 2581 !-- because the particle remains on the PE. 2582 trlp_count = 0 2583 trlpt_count = 0 2584 trrp_count = 0 2585 trrpt_count = 0 2586 IF ( pdims(1) /= 1 ) THEN 2587 ! 2588 !-- First calculate the storage necessary for sending and receiving the 2589 !-- data 2590 DO n = 1, number_of_particles 2591 i = ( particles(n)%x + 0.5 * dx ) * ddx 2592 ! 2593 !-- Above calculation does not work for indices less than zero 2594 IF ( particles(n)%x < -0.5 * dx ) i = -1 2595 2596 IF ( i < nxl ) THEN 2597 trlp_count = trlp_count + 1 2598 IF ( particles(n)%tail_id /= 0 ) trlpt_count = trlpt_count + 1 2599 ELSEIF ( i > nxr ) THEN 2600 trrp_count = trrp_count + 1 2601 IF ( particles(n)%tail_id /= 0 ) trrpt_count = trrpt_count + 1 2602 ENDIF 2603 ENDDO 2604 IF ( trlp_count == 0 ) trlp_count = 1 2605 IF ( trlpt_count == 0 ) trlpt_count = 1 2606 IF ( trrp_count == 0 ) trrp_count = 1 2607 IF ( trrpt_count == 0 ) trrpt_count = 1 2608 2609 ALLOCATE( trlp(trlp_count), trrp(trrp_count) ) 2610 2611 trlp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 2612 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 2613 0.0, 0, 0, 0, 0 ) 2614 trrp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 2615 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 2616 0.0, 0, 0, 0, 0 ) 2617 2618 IF ( use_particle_tails ) THEN 2619 ALLOCATE( trlpt(maximum_number_of_tailpoints,5,trlpt_count), & 2620 trrpt(maximum_number_of_tailpoints,5,trrpt_count) ) 2621 tlength = maximum_number_of_tailpoints * 5 2622 ENDIF 2623 2624 trlp_count = 0 2625 trlpt_count = 0 2626 trrp_count = 0 2627 trrpt_count = 0 2628 2629 ENDIF 2630 2631 ! WRITE ( 9, * ) '*** advec_particles: ##1' 2632 ! CALL local_flush( 9 ) 2633 ! nd = 0 2634 ! DO n = 1, number_of_particles 2635 ! IF ( .NOT. particle_mask(n) ) nd = nd + 1 2636 ! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) & 2637 ! THEN 2638 ! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')' 2639 ! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')' 2640 ! CALL local_flush( 9 ) 2641 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 2642 ! ENDIF 2643 ! ENDDO 2644 ! IF ( nd /= deleted_particles ) THEN 2645 ! WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles 2646 ! CALL local_flush( 9 ) 2647 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 2648 ! ENDIF 2649 2650 DO n = 1, number_of_particles 2651 2652 nn = particles(n)%tail_id 2653 2654 i = ( particles(n)%x + 0.5 * dx ) * ddx 2655 ! 2656 !-- Above calculation does not work for indices less than zero 2657 IF ( particles(n)%x < - 0.5 * dx ) i = -1 2658 2659 IF ( i < nxl ) THEN 2660 IF ( i < 0 ) THEN 2661 ! 2662 !-- Apply boundary condition along x 2663 IF ( ibc_par_lr == 0 ) THEN 2664 ! 2665 !-- Cyclic condition 2666 IF ( pdims(1) == 1 ) THEN 2667 particles(n)%x = ( nx + 1 ) * dx + particles(n)%x 2668 particles(n)%origin_x = ( nx + 1 ) * dx + & 2669 particles(n)%origin_x 2670 IF ( use_particle_tails .AND. nn /= 0 ) THEN 2671 i = particles(n)%tailpoints 2672 particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx & 2673 + particle_tail_coordinates(1:i,1,nn) 2674 ENDIF 2675 ELSE 2676 trlp_count = trlp_count + 1 2677 trlp(trlp_count) = particles(n) 2678 trlp(trlp_count)%x = ( nx + 1 ) * dx + trlp(trlp_count)%x 2679 trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x + & 2680 ( nx + 1 ) * dx 2681 particle_mask(n) = .FALSE. 2682 deleted_particles = deleted_particles + 1 2683 2684 IF ( trlp(trlp_count)%x >= (nx + 0.5)* dx - 1.e-12 ) THEN 2685 trlp(trlp_count)%x = trlp(trlp_count)%x - 1.e-10 2686 trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x & 2687 - 1 2688 ENDIF 2689 2690 IF ( use_particle_tails .AND. nn /= 0 ) THEN 2691 trlpt_count = trlpt_count + 1 2692 trlpt(:,:,trlpt_count) = & 2693 particle_tail_coordinates(:,:,nn) 2694 trlpt(:,1,trlpt_count) = ( nx + 1 ) * dx + & 2695 trlpt(:,1,trlpt_count) 2696 tail_mask(nn) = .FALSE. 2697 deleted_tails = deleted_tails + 1 2698 ENDIF 2699 ENDIF 2700 2701 ELSEIF ( ibc_par_lr == 1 ) THEN 2702 ! 2703 !-- Particle absorption 2704 particle_mask(n) = .FALSE. 2705 deleted_particles = deleted_particles + 1 2706 IF ( use_particle_tails .AND. nn /= 0 ) THEN 2707 tail_mask(nn) = .FALSE. 2708 deleted_tails = deleted_tails + 1 2709 ENDIF 2710 2711 ELSEIF ( ibc_par_lr == 2 ) THEN 2712 ! 2713 !-- Particle reflection 2714 particles(n)%x = -particles(n)%x 2715 particles(n)%speed_x = -particles(n)%speed_x 2716 2717 ENDIF 2718 ELSE 2719 ! 2720 !-- Store particle data in the transfer array, which will be send 2721 !-- to the neighbouring PE 2722 trlp_count = trlp_count + 1 2723 trlp(trlp_count) = particles(n) 2724 particle_mask(n) = .FALSE. 2725 deleted_particles = deleted_particles + 1 2726 2727 IF ( use_particle_tails .AND. nn /= 0 ) THEN 2728 trlpt_count = trlpt_count + 1 2729 trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn) 2730 tail_mask(nn) = .FALSE. 2731 deleted_tails = deleted_tails + 1 2732 ENDIF 2733 ENDIF 2734 2735 ELSEIF ( i > nxr ) THEN 2736 IF ( i > nx ) THEN 2737 ! 2738 !-- Apply boundary condition along x 2739 IF ( ibc_par_lr == 0 ) THEN 2740 ! 2741 !-- Cyclic condition 2742 IF ( pdims(1) == 1 ) THEN 2743 particles(n)%x = particles(n)%x - ( nx + 1 ) * dx 2744 particles(n)%origin_x = particles(n)%origin_x - & 2745 ( nx + 1 ) * dx 2746 IF ( use_particle_tails .AND. nn /= 0 ) THEN 2747 i = particles(n)%tailpoints 2748 particle_tail_coordinates(1:i,1,nn) = - ( nx+1 ) * dx & 2749 + particle_tail_coordinates(1:i,1,nn) 2750 ENDIF 2751 ELSE 2752 trrp_count = trrp_count + 1 2753 trrp(trrp_count) = particles(n) 2754 trrp(trrp_count)%x = trrp(trrp_count)%x - ( nx + 1 ) * dx 2755 trrp(trrp_count)%origin_x = trrp(trrp_count)%origin_x - & 2756 ( nx + 1 ) * dx 2757 particle_mask(n) = .FALSE. 2758 deleted_particles = deleted_particles + 1 2759 2760 IF ( use_particle_tails .AND. nn /= 0 ) THEN 2761 trrpt_count = trrpt_count + 1 2762 trrpt(:,:,trrpt_count) = & 2763 particle_tail_coordinates(:,:,nn) 2764 trrpt(:,1,trrpt_count) = trrpt(:,1,trrpt_count) - & 2765 ( nx + 1 ) * dx 2766 tail_mask(nn) = .FALSE. 2767 deleted_tails = deleted_tails + 1 2768 ENDIF 2769 ENDIF 2770 2771 ELSEIF ( ibc_par_lr == 1 ) THEN 2772 ! 2773 !-- Particle absorption 2774 particle_mask(n) = .FALSE. 2775 deleted_particles = deleted_particles + 1 2776 IF ( use_particle_tails .AND. nn /= 0 ) THEN 2777 tail_mask(nn) = .FALSE. 2778 deleted_tails = deleted_tails + 1 2779 ENDIF 2780 2781 ELSEIF ( ibc_par_lr == 2 ) THEN 2782 ! 2783 !-- Particle reflection 2784 particles(n)%x = 2 * ( nx * dx ) - particles(n)%x 2785 particles(n)%speed_x = -particles(n)%speed_x 2786 2787 ENDIF 2788 ELSE 2789 ! 2790 !-- Store particle data in the transfer array, which will be send 2791 !-- to the neighbouring PE 2792 trrp_count = trrp_count + 1 2793 trrp(trrp_count) = particles(n) 2794 particle_mask(n) = .FALSE. 2795 deleted_particles = deleted_particles + 1 2796 2797 IF ( use_particle_tails .AND. nn /= 0 ) THEN 2798 trrpt_count = trrpt_count + 1 2799 trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn) 2800 tail_mask(nn) = .FALSE. 2801 deleted_tails = deleted_tails + 1 2802 ENDIF 2803 ENDIF 2804 2805 ENDIF 2806 ENDDO 2807 2808 ! WRITE ( 9, * ) '*** advec_particles: ##2' 2809 ! CALL local_flush( 9 ) 2810 ! nd = 0 2811 ! DO n = 1, number_of_particles 2812 ! IF ( .NOT. particle_mask(n) ) nd = nd + 1 2813 ! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) & 2814 ! THEN 2815 ! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')' 2816 ! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')' 2817 ! CALL local_flush( 9 ) 2818 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 2819 ! ENDIF 2820 ! ENDDO 2821 ! IF ( nd /= deleted_particles ) THEN 2822 ! WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles 2823 ! CALL local_flush( 9 ) 2824 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 2825 ! ENDIF 2826 2827 ! 2828 !-- Send left boundary, receive right boundary (but first exchange how many 2829 !-- and check, if particle storage must be extended) 2830 IF ( pdims(1) /= 1 ) THEN 2831 2832 CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'start' ) 2833 CALL MPI_SENDRECV( trlp_count, 1, MPI_INTEGER, pleft, 0, & 2834 trrp_count_recv, 1, MPI_INTEGER, pright, 0, & 2835 comm2d, status, ierr ) 2836 2837 IF ( number_of_particles + trrp_count_recv > & 2838 maximum_number_of_particles ) & 2839 THEN 2840 IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN 2841 message_string = 'maximum_number_of_particles ' // & 2842 'needs to be increased ' // & 2843 '&but this is not allowed with ' // & 2844 'netcdf-data_format < 3' 2845 CALL message( 'advec_particles', 'PA0146', 2, 2, -1, 6, 1 ) 2846 ELSE 2847 ! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trrp' 2848 ! CALL local_flush( 9 ) 2849 CALL allocate_prt_memory( trrp_count_recv ) 2850 ! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trrp' 2851 ! CALL local_flush( 9 ) 2852 ENDIF 2853 ENDIF 2854 2855 CALL MPI_SENDRECV( trlp(1)%age, trlp_count, mpi_particle_type, & 2856 pleft, 1, particles(number_of_particles+1)%age, & 2857 trrp_count_recv, mpi_particle_type, pright, 1, & 2858 comm2d, status, ierr ) 2859 2860 IF ( use_particle_tails ) THEN 2861 2862 CALL MPI_SENDRECV( trlpt_count, 1, MPI_INTEGER, pleft, 0, & 2863 trrpt_count_recv, 1, MPI_INTEGER, pright, 0, & 2864 comm2d, status, ierr ) 2865 2866 IF ( number_of_tails+trrpt_count_recv > maximum_number_of_tails ) & 2867 THEN 2868 IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN 2869 message_string = 'maximum_number_of_tails ' // & 2870 'needs to be increased ' // & 2871 '&but this is not allowed wi'// & 2872 'th netcdf_data_format < 3' 2873 CALL message( 'advec_particles', 'PA0147', 2, 2, -1, 6, 1 ) 2874 ELSE 2875 ! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trrpt' 2876 ! CALL local_flush( 9 ) 2877 CALL allocate_tail_memory( trrpt_count_recv ) 2878 ! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trrpt' 2879 ! CALL local_flush( 9 ) 2880 ENDIF 2881 ENDIF 2882 2883 CALL MPI_SENDRECV( trlpt(1,1,1), trlpt_count*tlength, MPI_REAL, & 2884 pleft, 1, & 2885 particle_tail_coordinates(1,1,number_of_tails+1), & 2886 trrpt_count_recv*tlength, MPI_REAL, pright, 1, & 2887 comm2d, status, ierr ) 2888 ! 2889 !-- Update the tail ids for the transferred particles 2890 nn = number_of_tails 2891 DO n = number_of_particles+1, number_of_particles+trrp_count_recv 2892 IF ( particles(n)%tail_id /= 0 ) THEN 2893 nn = nn + 1 2894 particles(n)%tail_id = nn 2895 ENDIF 2896 ENDDO 2897 2898 ENDIF 2899 2900 number_of_particles = number_of_particles + trrp_count_recv 2901 number_of_tails = number_of_tails + trrpt_count_recv 2902 ! IF ( number_of_particles /= number_of_tails ) THEN 2903 ! WRITE (9,*) '--- advec_particles: #3' 2904 ! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails 2905 ! CALL local_flush( 9 ) 2906 ! ENDIF 2907 2908 ! 2909 !-- Send right boundary, receive left boundary 2910 CALL MPI_SENDRECV( trrp_count, 1, MPI_INTEGER, pright, 0, & 2911 trlp_count_recv, 1, MPI_INTEGER, pleft, 0, & 2912 comm2d, status, ierr ) 2913 2914 IF ( number_of_particles + trlp_count_recv > & 2915 maximum_number_of_particles ) & 2916 THEN 2917 IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN 2918 message_string = 'maximum_number_of_particles ' // & 2919 'needs to be increased ' // & 2920 '&but this is not allowed with '// & 2921 'netcdf_data_format < 3' 2922 CALL message( 'advec_particles', 'PA0146', 2, 2, -1, 6, 1 ) 2923 ELSE 2924 ! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trlp' 2925 ! CALL local_flush( 9 ) 2926 CALL allocate_prt_memory( trlp_count_recv ) 2927 ! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trlp' 2928 ! CALL local_flush( 9 ) 2929 ENDIF 2930 ENDIF 2931 2932 CALL MPI_SENDRECV( trrp(1)%age, trrp_count, mpi_particle_type, & 2933 pright, 1, particles(number_of_particles+1)%age, & 2934 trlp_count_recv, mpi_particle_type, pleft, 1, & 2935 comm2d, status, ierr ) 2936 2937 IF ( use_particle_tails ) THEN 2938 2939 CALL MPI_SENDRECV( trrpt_count, 1, MPI_INTEGER, pright, 0, & 2940 trlpt_count_recv, 1, MPI_INTEGER, pleft, 0, & 2941 comm2d, status, ierr ) 2942 2943 IF ( number_of_tails+trlpt_count_recv > maximum_number_of_tails ) & 2944 THEN 2945 IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN 2946 message_string = 'maximum_number_of_tails ' // & 2947 'needs to be increased ' // & 2948 '&but this is not allowed wi'// & 2949 'th netcdf_data_format < 3' 2950 CALL message( 'advec_particles', 'PA0147', 2, 2, -1, 6, 1 ) 2951 ELSE 2952 ! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trlpt' 2953 ! CALL local_flush( 9 ) 2954 CALL allocate_tail_memory( trlpt_count_recv ) 2955 ! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trlpt' 2956 ! CALL local_flush( 9 ) 2957 ENDIF 2958 ENDIF 2959 2960 CALL MPI_SENDRECV( trrpt(1,1,1), trrpt_count*tlength, MPI_REAL, & 2961 pright, 1, & 2962 particle_tail_coordinates(1,1,number_of_tails+1), & 2963 trlpt_count_recv*tlength, MPI_REAL, pleft, 1, & 2964 comm2d, status, ierr ) 2965 ! 2966 !-- Update the tail ids for the transferred particles 2967 nn = number_of_tails 2968 DO n = number_of_particles+1, number_of_particles+trlp_count_recv 2969 IF ( particles(n)%tail_id /= 0 ) THEN 2970 nn = nn + 1 2971 particles(n)%tail_id = nn 2972 ENDIF 2973 ENDDO 2974 2975 ENDIF 2976 2977 number_of_particles = number_of_particles + trlp_count_recv 2978 number_of_tails = number_of_tails + trlpt_count_recv 2979 ! IF ( number_of_particles /= number_of_tails ) THEN 2980 ! WRITE (9,*) '--- advec_particles: #4' 2981 ! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails 2982 ! CALL local_flush( 9 ) 2983 ! ENDIF 2984 2985 IF ( use_particle_tails ) THEN 2986 DEALLOCATE( trlpt, trrpt ) 2987 ENDIF 2988 DEALLOCATE( trlp, trrp ) 2989 2990 CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'pause' ) 2991 2992 ENDIF 2993 2994 ! WRITE ( 9, * ) '*** advec_particles: ##3' 2995 ! CALL local_flush( 9 ) 2996 ! nd = 0 2997 ! DO n = 1, number_of_particles 2998 ! IF ( .NOT. particle_mask(n) ) nd = nd + 1 2999 ! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) & 3000 ! THEN 3001 ! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')' 3002 ! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')' 3003 ! CALL local_flush( 9 ) 3004 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 3005 ! ENDIF 3006 ! ENDDO 3007 ! IF ( nd /= deleted_particles ) THEN 3008 ! WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles 3009 ! CALL local_flush( 9 ) 3010 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 3011 ! ENDIF 3012 3013 ! 3014 !-- Check whether particles have crossed the boundaries in y direction. Note 3015 !-- that this case can also apply to particles that have just been received 3016 !-- from the adjacent right or left PE. 3017 !-- Find out first the number of particles to be transferred and allocate 3018 !-- temporary arrays needed to store them. 3019 !-- For a one-dimensional decomposition along x, no transfer is necessary, 3020 !-- because the particle remains on the PE. 3021 trsp_count = 0 3022 trspt_count = 0 3023 trnp_count = 0 3024 trnpt_count = 0 3025 IF ( pdims(2) /= 1 ) THEN 3026 ! 3027 !-- First calculate the storage necessary for sending and receiving the 3028 !-- data 3029 DO n = 1, number_of_particles 3030 IF ( particle_mask(n) ) THEN 3031 j = ( particles(n)%y + 0.5 * dy ) * ddy 3032 ! 3033 !-- Above calculation does not work for indices less than zero 3034 IF ( particles(n)%y < -0.5 * dy ) j = -1 3035 3036 IF ( j < nys ) THEN 3037 trsp_count = trsp_count + 1 3038 IF ( particles(n)%tail_id /= 0 ) trspt_count = trspt_count+1 3039 ELSEIF ( j > nyn ) THEN 3040 trnp_count = trnp_count + 1 3041 IF ( particles(n)%tail_id /= 0 ) trnpt_count = trnpt_count+1 3042 ENDIF 3043 ENDIF 3044 ENDDO 3045 IF ( trsp_count == 0 ) trsp_count = 1 3046 IF ( trspt_count == 0 ) trspt_count = 1 3047 IF ( trnp_count == 0 ) trnp_count = 1 3048 IF ( trnpt_count == 0 ) trnpt_count = 1 3049 3050 ALLOCATE( trsp(trsp_count), trnp(trnp_count) ) 3051 3052 trsp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 3053 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 3054 0.0, 0, 0, 0, 0 ) 3055 trnp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 3056 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 3057 0.0, 0, 0, 0, 0 ) 3058 3059 IF ( use_particle_tails ) THEN 3060 ALLOCATE( trspt(maximum_number_of_tailpoints,5,trspt_count), & 3061 trnpt(maximum_number_of_tailpoints,5,trnpt_count) ) 3062 tlength = maximum_number_of_tailpoints * 5 3063 ENDIF 3064 3065 trsp_count = 0 3066 trspt_count = 0 3067 trnp_count = 0 3068 trnpt_count = 0 3069 3070 ENDIF 3071 3072 ! WRITE ( 9, * ) '*** advec_particles: ##4' 3073 ! CALL local_flush( 9 ) 3074 ! nd = 0 3075 ! DO n = 1, number_of_particles 3076 ! IF ( .NOT. particle_mask(n) ) nd = nd + 1 3077 ! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) & 3078 ! THEN 3079 ! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')' 3080 ! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')' 3081 ! CALL local_flush( 9 ) 3082 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 3083 ! ENDIF 3084 ! ENDDO 3085 ! IF ( nd /= deleted_particles ) THEN 3086 ! WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles 3087 ! CALL local_flush( 9 ) 3088 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 3089 ! ENDIF 3090 3091 DO n = 1, number_of_particles 3092 3093 nn = particles(n)%tail_id 3094 ! 3095 !-- Only those particles that have not been marked as 'deleted' may be 3096 !-- moved. 3097 IF ( particle_mask(n) ) THEN 3098 j = ( particles(n)%y + 0.5 * dy ) * ddy 3099 ! 3100 !-- Above calculation does not work for indices less than zero 3101 IF ( particles(n)%y < -0.5 * dy ) j = -1 3102 3103 IF ( j < nys ) THEN 3104 IF ( j < 0 ) THEN 3105 ! 3106 !-- Apply boundary condition along y 3107 IF ( ibc_par_ns == 0 ) THEN 3108 ! 3109 !-- Cyclic condition 3110 IF ( pdims(2) == 1 ) THEN 3111 particles(n)%y = ( ny + 1 ) * dy + particles(n)%y 3112 particles(n)%origin_y = ( ny + 1 ) * dy + & 3113 particles(n)%origin_y 3114 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3115 i = particles(n)%tailpoints 3116 particle_tail_coordinates(1:i,2,nn) = ( ny+1 ) * dy& 3117 + particle_tail_coordinates(1:i,2,nn) 3118 ENDIF 3119 ELSE 3120 trsp_count = trsp_count + 1 3121 trsp(trsp_count) = particles(n) 3122 trsp(trsp_count)%y = ( ny + 1 ) * dy + & 3123 trsp(trsp_count)%y 3124 trsp(trsp_count)%origin_y = trsp(trsp_count)%origin_y & 3125 + ( ny + 1 ) * dy 3126 particle_mask(n) = .FALSE. 3127 deleted_particles = deleted_particles + 1 3128 3129 IF ( trsp(trsp_count)%y >= (ny+0.5)* dy - 1.e-12 ) THEN 3130 trsp(trsp_count)%y = trsp(trsp_count)%y - 1.e-10 3131 trsp(trsp_count)%origin_y = & 3132 trsp(trsp_count)%origin_y - 1 3133 ENDIF 3134 3135 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3136 trspt_count = trspt_count + 1 3137 trspt(:,:,trspt_count) = & 3138 particle_tail_coordinates(:,:,nn) 3139 trspt(:,2,trspt_count) = ( ny + 1 ) * dy + & 3140 trspt(:,2,trspt_count) 3141 tail_mask(nn) = .FALSE. 3142 deleted_tails = deleted_tails + 1 3143 ENDIF 3144 ENDIF 3145 3146 ELSEIF ( ibc_par_ns == 1 ) THEN 3147 ! 3148 !-- Particle absorption 3149 particle_mask(n) = .FALSE. 3150 deleted_particles = deleted_particles + 1 3151 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3152 tail_mask(nn) = .FALSE. 3153 deleted_tails = deleted_tails + 1 3154 ENDIF 3155 3156 ELSEIF ( ibc_par_ns == 2 ) THEN 3157 ! 3158 !-- Particle reflection 3159 particles(n)%y = -particles(n)%y 3160 particles(n)%speed_y = -particles(n)%speed_y 3161 3162 ENDIF 3163 ELSE 3164 ! 3165 !-- Store particle data in the transfer array, which will be send 3166 !-- to the neighbouring PE 3167 trsp_count = trsp_count + 1 3168 trsp(trsp_count) = particles(n) 3169 particle_mask(n) = .FALSE. 3170 deleted_particles = deleted_particles + 1 3171 3172 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3173 trspt_count = trspt_count + 1 3174 trspt(:,:,trspt_count) = particle_tail_coordinates(:,:,nn) 3175 tail_mask(nn) = .FALSE. 3176 deleted_tails = deleted_tails + 1 3177 ENDIF 3178 ENDIF 3179 3180 ELSEIF ( j > nyn ) THEN 3181 IF ( j > ny ) THEN 3182 ! 3183 !-- Apply boundary condition along x 3184 IF ( ibc_par_ns == 0 ) THEN 3185 ! 3186 !-- Cyclic condition 3187 IF ( pdims(2) == 1 ) THEN 3188 particles(n)%y = particles(n)%y - ( ny + 1 ) * dy 3189 particles(n)%origin_y = particles(n)%origin_y - & 3190 ( ny + 1 ) * dy 3191 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3192 i = particles(n)%tailpoints 3193 particle_tail_coordinates(1:i,2,nn) = - (ny+1) * dy& 3194 + particle_tail_coordinates(1:i,2,nn) 3195 ENDIF 3196 ELSE 3197 trnp_count = trnp_count + 1 3198 trnp(trnp_count) = particles(n) 3199 trnp(trnp_count)%y = trnp(trnp_count)%y - & 3200 ( ny + 1 ) * dy 3201 trnp(trnp_count)%origin_y = trnp(trnp_count)%origin_y & 3202 - ( ny + 1 ) * dy 3203 particle_mask(n) = .FALSE. 3204 deleted_particles = deleted_particles + 1 3205 3206 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3207 trnpt_count = trnpt_count + 1 3208 trnpt(:,:,trnpt_count) = & 3209 particle_tail_coordinates(:,:,nn) 3210 trnpt(:,2,trnpt_count) = trnpt(:,2,trnpt_count) - & 3211 ( ny + 1 ) * dy 3212 tail_mask(nn) = .FALSE. 3213 deleted_tails = deleted_tails + 1 3214 ENDIF 3215 ENDIF 3216 3217 ELSEIF ( ibc_par_ns == 1 ) THEN 3218 ! 3219 !-- Particle absorption 3220 particle_mask(n) = .FALSE. 3221 deleted_particles = deleted_particles + 1 3222 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3223 tail_mask(nn) = .FALSE. 3224 deleted_tails = deleted_tails + 1 3225 ENDIF 3226 3227 ELSEIF ( ibc_par_ns == 2 ) THEN 3228 ! 3229 !-- Particle reflection 3230 particles(n)%y = 2 * ( ny * dy ) - particles(n)%y 3231 particles(n)%speed_y = -particles(n)%speed_y 3232 3233 ENDIF 3234 ELSE 3235 ! 3236 !-- Store particle data in the transfer array, which will be send 3237 !-- to the neighbouring PE 3238 trnp_count = trnp_count + 1 3239 trnp(trnp_count) = particles(n) 3240 particle_mask(n) = .FALSE. 3241 deleted_particles = deleted_particles + 1 3242 3243 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3244 trnpt_count = trnpt_count + 1 3245 trnpt(:,:,trnpt_count) = particle_tail_coordinates(:,:,nn) 3246 tail_mask(nn) = .FALSE. 3247 deleted_tails = deleted_tails + 1 3248 ENDIF 3249 ENDIF 3250 3251 ENDIF 3252 ENDIF 3253 ENDDO 3254 3255 ! WRITE ( 9, * ) '*** advec_particles: ##5' 3256 ! CALL local_flush( 9 ) 3257 ! nd = 0 3258 ! DO n = 1, number_of_particles 3259 ! IF ( .NOT. particle_mask(n) ) nd = nd + 1 3260 ! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) & 3261 ! THEN 3262 ! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')' 3263 ! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')' 3264 ! CALL local_flush( 9 ) 3265 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 3266 ! ENDIF 3267 ! ENDDO 3268 ! IF ( nd /= deleted_particles ) THEN 3269 ! WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles 3270 ! CALL local_flush( 9 ) 3271 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 3272 ! ENDIF 3273 3274 ! 3275 !-- Send front boundary, receive back boundary (but first exchange how many 3276 !-- and check, if particle storage must be extended) 3277 IF ( pdims(2) /= 1 ) THEN 3278 3279 CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'continue' ) 3280 CALL MPI_SENDRECV( trsp_count, 1, MPI_INTEGER, psouth, 0, & 3281 trnp_count_recv, 1, MPI_INTEGER, pnorth, 0, & 3282 comm2d, status, ierr ) 3283 3284 IF ( number_of_particles + trnp_count_recv > & 3285 maximum_number_of_particles ) & 3286 THEN 3287 IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN 3288 message_string = 'maximum_number_of_particles ' // & 3289 'needs to be increased ' // & 3290 '&but this is not allowed with '// & 3291 'netcdf_data_format < 3' 3292 CALL message( 'advec_particles', 'PA0146', 2, 2, -1, 6, 1 ) 3293 ELSE 3294 ! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trnp' 3295 ! CALL local_flush( 9 ) 3296 CALL allocate_prt_memory( trnp_count_recv ) 3297 ! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trnp' 3298 ! CALL local_flush( 9 ) 3299 ENDIF 3300 ENDIF 3301 3302 CALL MPI_SENDRECV( trsp(1)%age, trsp_count, mpi_particle_type, & 3303 psouth, 1, particles(number_of_particles+1)%age, & 3304 trnp_count_recv, mpi_particle_type, pnorth, 1, & 3305 comm2d, status, ierr ) 3306 3307 IF ( use_particle_tails ) THEN 3308 3309 CALL MPI_SENDRECV( trspt_count, 1, MPI_INTEGER, psouth, 0, & 3310 trnpt_count_recv, 1, MPI_INTEGER, pnorth, 0, & 3311 comm2d, status, ierr ) 3312 3313 IF ( number_of_tails+trnpt_count_recv > maximum_number_of_tails ) & 3314 THEN 3315 IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN 3316 message_string = 'maximum_number_of_tails ' // & 3317 'needs to be increased ' // & 3318 '&but this is not allowed wi' // & 3319 'th netcdf_data_format < 3' 3320 CALL message( 'advec_particles', 'PA0147', 2, 2, -1, 6, 1 ) 3321 ELSE 3322 ! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trnpt' 3323 ! CALL local_flush( 9 ) 3324 CALL allocate_tail_memory( trnpt_count_recv ) 3325 ! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trnpt' 3326 ! CALL local_flush( 9 ) 3327 ENDIF 3328 ENDIF 3329 3330 CALL MPI_SENDRECV( trspt(1,1,1), trspt_count*tlength, MPI_REAL, & 3331 psouth, 1, & 3332 particle_tail_coordinates(1,1,number_of_tails+1), & 3333 trnpt_count_recv*tlength, MPI_REAL, pnorth, 1, & 3334 comm2d, status, ierr ) 3335 3336 ! 3337 !-- Update the tail ids for the transferred particles 3338 nn = number_of_tails 3339 DO n = number_of_particles+1, number_of_particles+trnp_count_recv 3340 IF ( particles(n)%tail_id /= 0 ) THEN 3341 nn = nn + 1 3342 particles(n)%tail_id = nn 3343 ENDIF 3344 ENDDO 3345 3346 ENDIF 3347 3348 number_of_particles = number_of_particles + trnp_count_recv 3349 number_of_tails = number_of_tails + trnpt_count_recv 3350 ! IF ( number_of_particles /= number_of_tails ) THEN 3351 ! WRITE (9,*) '--- advec_particles: #5' 3352 ! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails 3353 ! CALL local_flush( 9 ) 3354 ! ENDIF 3355 3356 ! 3357 !-- Send back boundary, receive front boundary 3358 CALL MPI_SENDRECV( trnp_count, 1, MPI_INTEGER, pnorth, 0, & 3359 trsp_count_recv, 1, MPI_INTEGER, psouth, 0, & 3360 comm2d, status, ierr ) 3361 3362 IF ( number_of_particles + trsp_count_recv > & 3363 maximum_number_of_particles ) & 3364 THEN 3365 IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN 3366 message_string = 'maximum_number_of_particles ' // & 3367 'needs to be increased ' // & 3368 '&but this is not allowed with ' // & 3369 'netcdf_data_format < 3' 3370 CALL message( 'advec_particles', 'PA0146', 2, 2, -1, 6, 1 ) 3371 ELSE 3372 ! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trsp' 3373 ! CALL local_flush( 9 ) 3374 CALL allocate_prt_memory( trsp_count_recv ) 3375 ! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trsp' 3376 ! CALL local_flush( 9 ) 3377 ENDIF 3378 ENDIF 3379 3380 CALL MPI_SENDRECV( trnp(1)%age, trnp_count, mpi_particle_type, & 3381 pnorth, 1, particles(number_of_particles+1)%age, & 3382 trsp_count_recv, mpi_particle_type, psouth, 1, & 3383 comm2d, status, ierr ) 3384 3385 IF ( use_particle_tails ) THEN 3386 3387 CALL MPI_SENDRECV( trnpt_count, 1, MPI_INTEGER, pnorth, 0, & 3388 trspt_count_recv, 1, MPI_INTEGER, psouth, 0, & 3389 comm2d, status, ierr ) 3390 3391 IF ( number_of_tails+trspt_count_recv > maximum_number_of_tails ) & 3392 THEN 3393 IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN 3394 message_string = 'maximum_number_of_tails ' // & 3395 'needs to be increased ' // & 3396 '&but this is not allowed wi'// & 3397 'th NetCDF output switched on' 3398 CALL message( 'advec_particles', 'PA0147', 2, 2, -1, 6, 1 ) 3399 ELSE 3400 ! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trspt' 3401 ! CALL local_flush( 9 ) 3402 CALL allocate_tail_memory( trspt_count_recv ) 3403 ! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trspt' 3404 ! CALL local_flush( 9 ) 3405 ENDIF 3406 ENDIF 3407 3408 CALL MPI_SENDRECV( trnpt(1,1,1), trnpt_count*tlength, MPI_REAL, & 3409 pnorth, 1, & 3410 particle_tail_coordinates(1,1,number_of_tails+1), & 3411 trspt_count_recv*tlength, MPI_REAL, psouth, 1, & 3412 comm2d, status, ierr ) 3413 ! 3414 !-- Update the tail ids for the transferred particles 3415 nn = number_of_tails 3416 DO n = number_of_particles+1, number_of_particles+trsp_count_recv 3417 IF ( particles(n)%tail_id /= 0 ) THEN 3418 nn = nn + 1 3419 particles(n)%tail_id = nn 3420 ENDIF 3421 ENDDO 3422 3423 ENDIF 3424 3425 number_of_particles = number_of_particles + trsp_count_recv 3426 number_of_tails = number_of_tails + trspt_count_recv 3427 ! IF ( number_of_particles /= number_of_tails ) THEN 3428 ! WRITE (9,*) '--- advec_particles: #6' 3429 ! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails 3430 ! CALL local_flush( 9 ) 3431 ! ENDIF 3432 3433 IF ( use_particle_tails ) THEN 3434 DEALLOCATE( trspt, trnpt ) 3435 ENDIF 3436 DEALLOCATE( trsp, trnp ) 3437 3438 CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'stop' ) 3439 3440 ENDIF 3441 3442 ! WRITE ( 9, * ) '*** advec_particles: ##6' 3443 ! CALL local_flush( 9 ) 3444 ! nd = 0 3445 ! DO n = 1, number_of_particles 3446 ! IF ( .NOT. particle_mask(n) ) nd = nd + 1 3447 ! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) & 3448 ! THEN 3449 ! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')' 3450 ! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')' 3451 ! CALL local_flush( 9 ) 3452 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 3453 ! ENDIF 3454 ! ENDDO 3455 ! IF ( nd /= deleted_particles ) THEN 3456 ! WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles 3457 ! CALL local_flush( 9 ) 3458 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 3459 ! ENDIF 3460 3461 #else 3462 3463 ! 3464 !-- Apply boundary conditions 3465 DO n = 1, number_of_particles 3466 3467 nn = particles(n)%tail_id 3468 3469 IF ( particles(n)%x < -0.5 * dx ) THEN 3470 3471 IF ( ibc_par_lr == 0 ) THEN 3472 ! 3473 !-- Cyclic boundary. Relevant coordinate has to be changed. 3474 particles(n)%x = ( nx + 1 ) * dx + particles(n)%x 3475 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3476 i = particles(n)%tailpoints 3477 particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx + & 3478 particle_tail_coordinates(1:i,1,nn) 3479 ENDIF 3480 ELSEIF ( ibc_par_lr == 1 ) THEN 3481 ! 3482 !-- Particle absorption 3483 particle_mask(n) = .FALSE. 3484 deleted_particles = deleted_particles + 1 3485 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3486 tail_mask(nn) = .FALSE. 3487 deleted_tails = deleted_tails + 1 3488 ENDIF 3489 ELSEIF ( ibc_par_lr == 2 ) THEN 3490 ! 3491 !-- Particle reflection 3492 particles(n)%x = -dx - particles(n)%x 3493 particles(n)%speed_x = -particles(n)%speed_x 3494 ENDIF 3495 3496 ELSEIF ( particles(n)%x >= ( nx + 0.5 ) * dx ) THEN 3497 3498 IF ( ibc_par_lr == 0 ) THEN 3499 ! 3500 !-- Cyclic boundary. Relevant coordinate has to be changed. 3501 particles(n)%x = particles(n)%x - ( nx + 1 ) * dx 3502 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3503 i = particles(n)%tailpoints 3504 particle_tail_coordinates(1:i,1,nn) = - ( nx + 1 ) * dx + & 3505 particle_tail_coordinates(1:i,1,nn) 3506 ENDIF 3507 ELSEIF ( ibc_par_lr == 1 ) THEN 3508 ! 3509 !-- Particle absorption 3510 particle_mask(n) = .FALSE. 3511 deleted_particles = deleted_particles + 1 3512 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3513 tail_mask(nn) = .FALSE. 3514 deleted_tails = deleted_tails + 1 3515 ENDIF 3516 ELSEIF ( ibc_par_lr == 2 ) THEN 3517 ! 3518 !-- Particle reflection 3519 particles(n)%x = ( nx + 1 ) * dx - particles(n)%x 3520 particles(n)%speed_x = -particles(n)%speed_x 3521 ENDIF 3522 3523 ENDIF 3524 3525 IF ( particles(n)%y < -0.5 * dy ) THEN 3526 3527 IF ( ibc_par_ns == 0 ) THEN 3528 ! 3529 !-- Cyclic boundary. Relevant coordinate has to be changed. 3530 particles(n)%y = ( ny + 1 ) * dy + particles(n)%y 3531 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3532 i = particles(n)%tailpoints 3533 particle_tail_coordinates(1:i,2,nn) = ( ny + 1 ) * dy + & 3534 particle_tail_coordinates(1:i,2,nn) 3535 ENDIF 3536 ELSEIF ( ibc_par_ns == 1 ) THEN 3537 ! 3538 !-- Particle absorption 3539 particle_mask(n) = .FALSE. 3540 deleted_particles = deleted_particles + 1 3541 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3542 tail_mask(nn) = .FALSE. 3543 deleted_tails = deleted_tails + 1 3544 ENDIF 3545 ELSEIF ( ibc_par_ns == 2 ) THEN 3546 ! 3547 !-- Particle reflection 3548 particles(n)%y = -dy - particles(n)%y 3549 particles(n)%speed_y = -particles(n)%speed_y 3550 ENDIF 3551 3552 ELSEIF ( particles(n)%y >= ( ny + 0.5 ) * dy ) THEN 3553 3554 IF ( ibc_par_ns == 0 ) THEN 3555 ! 3556 !-- Cyclic boundary. Relevant coordinate has to be changed. 3557 particles(n)%y = particles(n)%y - ( ny + 1 ) * dy 3558 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3559 i = particles(n)%tailpoints 3560 particle_tail_coordinates(1:i,2,nn) = - ( ny + 1 ) * dy + & 3561 particle_tail_coordinates(1:i,2,nn) 3562 ENDIF 3563 ELSEIF ( ibc_par_ns == 1 ) THEN 3564 ! 3565 !-- Particle absorption 3566 particle_mask(n) = .FALSE. 3567 deleted_particles = deleted_particles + 1 3568 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3569 tail_mask(nn) = .FALSE. 3570 deleted_tails = deleted_tails + 1 3571 ENDIF 3572 ELSEIF ( ibc_par_ns == 2 ) THEN 3573 ! 3574 !-- Particle reflection 3575 particles(n)%y = ( ny + 1 ) * dy - particles(n)%y 3576 particles(n)%speed_y = -particles(n)%speed_y 3577 ENDIF 3578 3579 ENDIF 3580 ENDDO 3581 3582 #endif 221 ! 222 !-- Horizontal boundary conditions including exchange between subdmains 223 CALL lpm_exchange_horiz 3583 224 3584 225 ! 3585 226 !-- Apply boundary conditions to those particles that have crossed the top or 3586 227 !-- bottom boundary and delete those particles, which are older than allowed 3587 DO n = 1, number_of_particles 3588 3589 nn = particles(n)%tail_id 3590 3591 ! 3592 !-- Stop if particles have moved further than the length of one 3593 !-- PE subdomain (newly released particles have age = age_m!) 3594 IF ( particles(n)%age /= particles(n)%age_m ) THEN 3595 IF ( ABS(particles(n)%speed_x) > & 3596 ((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m) .OR. & 3597 ABS(particles(n)%speed_y) > & 3598 ((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) ) THEN 3599 3600 WRITE( message_string, * ) 'particle too fast. n = ', n 3601 CALL message( 'advec_particles', 'PA0148', 2, 2, -1, 6, 1 ) 3602 ENDIF 3603 ENDIF 3604 3605 IF ( particles(n)%age > particle_maximum_age .AND. & 3606 particle_mask(n) ) & 3607 THEN 3608 particle_mask(n) = .FALSE. 3609 deleted_particles = deleted_particles + 1 3610 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3611 tail_mask(nn) = .FALSE. 3612 deleted_tails = deleted_tails + 1 3613 ENDIF 3614 ENDIF 3615 3616 IF ( particles(n)%z >= zu(nz) .AND. particle_mask(n) ) THEN 3617 IF ( ibc_par_t == 1 ) THEN 3618 ! 3619 !-- Particle absorption 3620 particle_mask(n) = .FALSE. 3621 deleted_particles = deleted_particles + 1 3622 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3623 tail_mask(nn) = .FALSE. 3624 deleted_tails = deleted_tails + 1 3625 ENDIF 3626 ELSEIF ( ibc_par_t == 2 ) THEN 3627 ! 3628 !-- Particle reflection 3629 particles(n)%z = 2.0 * zu(nz) - particles(n)%z 3630 particles(n)%speed_z = -particles(n)%speed_z 3631 IF ( use_sgs_for_particles .AND. & 3632 particles(n)%rvar3 > 0.0 ) THEN 3633 particles(n)%rvar3 = -particles(n)%rvar3 3634 ENDIF 3635 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3636 particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - & 3637 particle_tail_coordinates(1,3,nn) 3638 ENDIF 3639 ENDIF 3640 ENDIF 3641 IF ( particles(n)%z < zw(0) .AND. particle_mask(n) ) THEN 3642 IF ( ibc_par_b == 1 ) THEN 3643 ! 3644 !-- Particle absorption 3645 particle_mask(n) = .FALSE. 3646 deleted_particles = deleted_particles + 1 3647 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3648 tail_mask(nn) = .FALSE. 3649 deleted_tails = deleted_tails + 1 3650 ENDIF 3651 ELSEIF ( ibc_par_b == 2 ) THEN 3652 ! 3653 !-- Particle reflection 3654 particles(n)%z = 2.0 * zw(0) - particles(n)%z 3655 particles(n)%speed_z = -particles(n)%speed_z 3656 IF ( use_sgs_for_particles .AND. & 3657 particles(n)%rvar3 < 0.0 ) THEN 3658 particles(n)%rvar3 = -particles(n)%rvar3 3659 ENDIF 3660 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3661 particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - & 3662 particle_tail_coordinates(1,3,nn) 3663 ENDIF 3664 IF ( use_particle_tails .AND. nn /= 0 ) THEN 3665 particle_tail_coordinates(1,3,nn) = 2.0 * zw(0) - & 3666 particle_tail_coordinates(1,3,nn) 3667 ENDIF 3668 ENDIF 3669 ENDIF 3670 ENDDO 3671 3672 ! WRITE ( 9, * ) '*** advec_particles: ##7' 3673 ! CALL local_flush( 9 ) 3674 ! nd = 0 3675 ! DO n = 1, number_of_particles 3676 ! IF ( .NOT. particle_mask(n) ) nd = nd + 1 3677 ! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) & 3678 ! THEN 3679 ! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')' 3680 ! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')' 3681 ! CALL local_flush( 9 ) 3682 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 3683 ! ENDIF 3684 ! ENDDO 3685 ! IF ( nd /= deleted_particles ) THEN 3686 ! WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles 3687 ! CALL local_flush( 9 ) 3688 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 3689 ! ENDIF 228 CALL lpm_boundary_conds( 'bottom/top' ) 3690 229 3691 230 ! … … 3693 232 !-- determine new number of particles 3694 233 IF ( number_of_particles > 0 .AND. deleted_particles > 0 ) THEN 3695 nn = 0 3696 nd = 0 3697 DO n = 1, number_of_particles 3698 IF ( particle_mask(n) ) THEN 3699 nn = nn + 1 3700 particles(nn) = particles(n) 3701 ELSE 3702 nd = nd + 1 3703 ENDIF 3704 ENDDO 3705 ! IF ( nd /= deleted_particles ) THEN 3706 ! WRITE (9,*) '*** advec_part nd=',nd,' deleted_particles=',deleted_particles 3707 ! CALL local_flush( 9 ) 3708 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 3709 ! ENDIF 3710 3711 number_of_particles = number_of_particles - deleted_particles 3712 ! 3713 !-- Pack the tails, store the new tail ids and re-assign it to the 3714 !-- respective 3715 !-- particles 3716 IF ( use_particle_tails ) THEN 3717 nn = 0 3718 nd = 0 3719 DO n = 1, number_of_tails 3720 IF ( tail_mask(n) ) THEN 3721 nn = nn + 1 3722 particle_tail_coordinates(:,:,nn) = & 3723 particle_tail_coordinates(:,:,n) 3724 new_tail_id(n) = nn 3725 ELSE 3726 nd = nd + 1 3727 ! WRITE (9,*) '+++ n=',n,' (of ',number_of_tails,' #oftails)' 3728 ! WRITE (9,*) ' id=',new_tail_id(n) 3729 ! CALL local_flush( 9 ) 3730 ENDIF 3731 ENDDO 3732 ENDIF 3733 3734 ! IF ( nd /= deleted_tails .AND. use_particle_tails ) THEN 3735 ! WRITE (9,*) '*** advec_part nd=',nd,' deleted_tails=',deleted_tails 3736 ! CALL local_flush( 9 ) 3737 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 3738 ! ENDIF 3739 3740 number_of_tails = number_of_tails - deleted_tails 3741 3742 ! nn = 0 3743 DO n = 1, number_of_particles 3744 IF ( particles(n)%tail_id /= 0 ) THEN 3745 ! nn = nn + 1 3746 ! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id > number_of_tails ) THEN 3747 ! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')' 3748 ! WRITE (9,*) ' tail_id=',particles(n)%tail_id 3749 ! WRITE (9,*) ' new_tail_id=', new_tail_id(particles(n)%tail_id), & 3750 ! ' of (',number_of_tails,')' 3751 ! CALL local_flush( 9 ) 3752 ! ENDIF 3753 particles(n)%tail_id = new_tail_id(particles(n)%tail_id) 3754 ENDIF 3755 ENDDO 3756 3757 ! IF ( nn /= number_of_tails .AND. use_particle_tails ) THEN 3758 ! WRITE (9,*) '*** advec_part #of_tails=',number_of_tails,' nn=',nn 3759 ! CALL local_flush( 9 ) 3760 ! DO n = 1, number_of_particles 3761 ! WRITE (9,*) 'prt# ',n,' tail_id=',particles(n)%tail_id, & 3762 ! ' x=',particles(n)%x, ' y=',particles(n)%y, & 3763 ! ' z=',particles(n)%z 3764 ! ENDDO 3765 ! CALL MPI_ABORT( comm2d, 9999, ierr ) 3766 ! ENDIF 3767 234 CALL lpm_pack_arrays 3768 235 ENDIF 3769 236 3770 ! IF ( number_of_particles /= number_of_tails ) THEN3771 ! WRITE (9,*) '--- advec_particles: #7'3772 ! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails3773 ! CALL local_flush( 9 )3774 ! ENDIF3775 ! WRITE ( 9, * ) '*** advec_particles: ##8'3776 ! CALL local_flush( 9 )3777 ! DO n = 1, number_of_particles3778 ! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &3779 ! THEN3780 ! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'3781 ! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'3782 ! CALL local_flush( 9 )3783 ! CALL MPI_ABORT( comm2d, 9999, ierr )3784 ! ENDIF3785 ! ENDDO3786 3787 ! WRITE ( 9, * ) '*** advec_particles: ##9'3788 ! CALL local_flush( 9 )3789 ! DO n = 1, number_of_particles3790 ! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &3791 ! THEN3792 ! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'3793 ! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'3794 ! CALL local_flush( 9 )3795 ! CALL MPI_ABORT( comm2d, 9999, ierr )3796 ! ENDIF3797 ! ENDDO3798 3799 !3800 !-- Accumulate the number of particles transferred between the subdomains3801 #if defined( __parallel )3802 trlp_count_sum = trlp_count_sum + trlp_count3803 trlp_count_recv_sum = trlp_count_recv_sum + trlp_count_recv3804 trrp_count_sum = trrp_count_sum + trrp_count3805 trrp_count_recv_sum = trrp_count_recv_sum + trrp_count_recv3806 trsp_count_sum = trsp_count_sum + trsp_count3807 trsp_count_recv_sum = trsp_count_recv_sum + trsp_count_recv3808 trnp_count_sum = trnp_count_sum + trnp_count3809 trnp_count_recv_sum = trnp_count_recv_sum + trnp_count_recv3810 #endif3811 3812 237 IF ( dt_3d_reached ) EXIT 3813 238 3814 !3815 !-- Initialize variables for the next (sub-) timestep, i.e. for marking those3816 !-- particles to be deleted after the timestep3817 particle_mask = .TRUE.3818 deleted_particles = 03819 trlp_count_recv = 03820 trnp_count_recv = 03821 trrp_count_recv = 03822 trsp_count_recv = 03823 trlpt_count_recv = 03824 trnpt_count_recv = 03825 trrpt_count_recv = 03826 trspt_count_recv = 03827 IF ( use_particle_tails ) THEN3828 tail_mask = .TRUE.3829 ENDIF3830 deleted_tails = 03831 3832 239 ENDDO ! timestep loop 240 3833 241 3834 242 ! … … 3836 244 time_sort_particles = time_sort_particles + dt_3d 3837 245 IF ( time_sort_particles >= dt_sort_particles ) THEN 3838 CALL sort_particles246 CALL lpm_sort_arrays 3839 247 time_sort_particles = MOD( time_sort_particles, & 3840 248 MAX( dt_sort_particles, dt_3d ) ) 3841 249 ENDIF 3842 250 3843 IF ( cloud_droplets ) THEN 3844 3845 CALL cpu_log( log_point_s(45), 'advec_part_reeval_we', 'start' ) 3846 3847 ql = 0.0; ql_v = 0.0; ql_vp = 0.0 3848 3849 ! 3850 !-- Calculate the liquid water content 3851 DO i = nxl, nxr 3852 DO j = nys, nyn 3853 DO k = nzb, nzt+1 3854 3855 ! 3856 !-- Calculate the total volume in the boxes (ql_v, weighting factor 3857 !-- has to beincluded) 3858 psi = prt_start_index(k,j,i) 3859 DO n = psi, psi+prt_count(k,j,i)-1 3860 ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * & 3861 particles(n)%radius**3 3862 ENDDO 3863 3864 ! 3865 !-- Calculate the liquid water content 3866 IF ( ql_v(k,j,i) /= 0.0 ) THEN 3867 ql(k,j,i) = ql(k,j,i) + rho_l * 1.33333333 * pi * & 3868 ql_v(k,j,i) / & 3869 ( rho_surface * dx * dy * dz ) 3870 3871 IF ( ql(k,j,i) < 0.0 ) THEN 3872 WRITE( message_string, * ) 'LWC out of range: ' , & 3873 ql(k,j,i) 3874 CALL message( 'advec_particles', '', 2, 2, -1, 6, 1 ) 3875 ENDIF 3876 3877 ELSE 3878 ql(k,j,i) = 0.0 3879 ENDIF 3880 3881 ENDDO 3882 ENDDO 3883 ENDDO 3884 3885 CALL cpu_log( log_point_s(45), 'advec_part_reeval_we', 'stop' ) 3886 3887 ENDIF 251 252 ! 253 !-- Calculate the new liquid water content for each grid box 254 IF ( cloud_droplets ) CALL lpm_calc_liquid_water_content 255 3888 256 3889 257 ! … … 3892 260 !-- particle attribute (class) is then used for storing the particle radius 3893 261 !-- class. 3894 IF ( collision_kernel == 'none' ) CALL set_particle_attributes 262 IF ( collision_kernel == 'none' ) CALL lpm_set_attributes 263 3895 264 3896 265 ! 3897 266 !-- Set particle attributes defined by the user 3898 CALL user_particle_attributes 3899 3900 ! 3901 !-- If necessary, add the actual particle positions to the particle tails 3902 IF ( use_particle_tails ) THEN 3903 3904 distance = 0.0 3905 DO n = 1, number_of_particles 3906 3907 nn = particles(n)%tail_id 3908 3909 IF ( nn /= 0 ) THEN 3910 ! 3911 !-- Calculate the distance between the actual particle position and the 3912 !-- next tailpoint 3913 ! WRITE ( 9, * ) '*** advec_particles: ##10.1 nn=',nn 3914 ! CALL local_flush( 9 ) 3915 IF ( minimum_tailpoint_distance /= 0.0 ) THEN 3916 distance = ( particle_tail_coordinates(1,1,nn) - & 3917 particle_tail_coordinates(2,1,nn) )**2 + & 3918 ( particle_tail_coordinates(1,2,nn) - & 3919 particle_tail_coordinates(2,2,nn) )**2 + & 3920 ( particle_tail_coordinates(1,3,nn) - & 3921 particle_tail_coordinates(2,3,nn) )**2 3922 ENDIF 3923 ! WRITE ( 9, * ) '*** advec_particles: ##10.2' 3924 ! CALL local_flush( 9 ) 3925 ! 3926 !-- First, increase the index of all existings tailpoints by one 3927 IF ( distance >= minimum_tailpoint_distance ) THEN 3928 DO i = particles(n)%tailpoints, 1, -1 3929 particle_tail_coordinates(i+1,:,nn) = & 3930 particle_tail_coordinates(i,:,nn) 3931 ENDDO 3932 ! 3933 !-- Increase the counter which contains the number of tailpoints. 3934 !-- This must always be smaller than the given maximum number of 3935 !-- tailpoints because otherwise the index bounds of 3936 !-- particle_tail_coordinates would be exceeded 3937 IF ( particles(n)%tailpoints < maximum_number_of_tailpoints-1 )& 3938 THEN 3939 particles(n)%tailpoints = particles(n)%tailpoints + 1 3940 ENDIF 3941 ENDIF 3942 ! WRITE ( 9, * ) '*** advec_particles: ##10.3' 3943 ! CALL local_flush( 9 ) 3944 ! 3945 !-- In any case, store the new point at the beginning of the tail 3946 particle_tail_coordinates(1,1,nn) = particles(n)%x 3947 particle_tail_coordinates(1,2,nn) = particles(n)%y 3948 particle_tail_coordinates(1,3,nn) = particles(n)%z 3949 particle_tail_coordinates(1,4,nn) = particles(n)%class 3950 ! WRITE ( 9, * ) '*** advec_particles: ##10.4' 3951 ! CALL local_flush( 9 ) 3952 ! 3953 !-- Increase the age of the tailpoints 3954 IF ( minimum_tailpoint_distance /= 0.0 ) THEN 3955 particle_tail_coordinates(2:particles(n)%tailpoints,5,nn) = & 3956 particle_tail_coordinates(2:particles(n)%tailpoints,5,nn) + & 3957 dt_3d 3958 ! 3959 !-- Delete the last tailpoint, if it has exceeded its maximum age 3960 IF ( particle_tail_coordinates(particles(n)%tailpoints,5,nn) > & 3961 maximum_tailpoint_age ) THEN 3962 particles(n)%tailpoints = particles(n)%tailpoints - 1 3963 ENDIF 3964 ENDIF 3965 ! WRITE ( 9, * ) '*** advec_particles: ##10.5' 3966 ! CALL local_flush( 9 ) 3967 3968 ENDIF 3969 3970 ENDDO 3971 3972 ENDIF 3973 ! WRITE ( 9, * ) '*** advec_particles: ##11' 3974 ! CALL local_flush( 9 ) 3975 3976 ! 3977 !-- Write particle statistics on file 3978 IF ( write_particle_statistics ) THEN 3979 CALL check_open( 80 ) 3980 #if defined( __parallel ) 3981 WRITE ( 80, 8000 ) current_timestep_number+1, simulated_time+dt_3d, & 3982 number_of_particles, pleft, trlp_count_sum, & 3983 trlp_count_recv_sum, pright, trrp_count_sum, & 3984 trrp_count_recv_sum, psouth, trsp_count_sum, & 3985 trsp_count_recv_sum, pnorth, trnp_count_sum, & 3986 trnp_count_recv_sum, maximum_number_of_particles 3987 CALL close_file( 80 ) 3988 #else 3989 WRITE ( 80, 8000 ) current_timestep_number+1, simulated_time+dt_3d, & 3990 number_of_particles, maximum_number_of_particles 3991 #endif 3992 ENDIF 3993 3994 CALL cpu_log( log_point(25), 'advec_particles', 'stop' ) 3995 3996 ! 3997 !-- Formats 3998 8000 FORMAT (I6,1X,F7.2,4X,I6,5X,4(I3,1X,I4,'/',I4,2X),6X,I6) 3999 4000 END SUBROUTINE advec_particles 4001 4002 4003 SUBROUTINE allocate_prt_memory( number_of_new_particles ) 4004 4005 !------------------------------------------------------------------------------! 4006 ! Description: 4007 ! ------------ 4008 ! Extend particle memory 4009 !------------------------------------------------------------------------------! 4010 4011 USE particle_attributes 4012 4013 IMPLICIT NONE 4014 4015 INTEGER :: new_maximum_number, number_of_new_particles 4016 4017 LOGICAL, DIMENSION(:), ALLOCATABLE :: tmp_particle_mask 4018 4019 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles 4020 4021 4022 new_maximum_number = maximum_number_of_particles + & 4023 MAX( 5*number_of_new_particles, number_of_initial_particles ) 4024 4025 IF ( write_particle_statistics ) THEN 4026 CALL check_open( 80 ) 4027 WRITE ( 80, '(''*** Request: '', I7, '' new_maximum_number(prt)'')' ) & 4028 new_maximum_number 4029 CALL close_file( 80 ) 4030 ENDIF 4031 4032 ALLOCATE( tmp_particles(maximum_number_of_particles), & 4033 tmp_particle_mask(maximum_number_of_particles) ) 4034 tmp_particles = particles 4035 tmp_particle_mask = particle_mask 4036 4037 DEALLOCATE( particles, particle_mask ) 4038 ALLOCATE( particles(new_maximum_number), & 4039 particle_mask(new_maximum_number) ) 4040 maximum_number_of_particles = new_maximum_number 4041 4042 particles(1:number_of_particles) = tmp_particles(1:number_of_particles) 4043 particle_mask(1:number_of_particles) = & 4044 tmp_particle_mask(1:number_of_particles) 4045 particle_mask(number_of_particles+1:maximum_number_of_particles) = .TRUE. 4046 DEALLOCATE( tmp_particles, tmp_particle_mask ) 4047 4048 END SUBROUTINE allocate_prt_memory 4049 4050 4051 SUBROUTINE allocate_tail_memory( number_of_new_tails ) 4052 4053 !------------------------------------------------------------------------------! 4054 ! Description: 4055 ! ------------ 4056 ! Extend tail memory 4057 !------------------------------------------------------------------------------! 4058 4059 USE particle_attributes 4060 4061 IMPLICIT NONE 4062 4063 INTEGER :: new_maximum_number, number_of_new_tails 4064 4065 LOGICAL, DIMENSION(maximum_number_of_tails) :: tmp_tail_mask 4066 4067 REAL, DIMENSION(maximum_number_of_tailpoints,5,maximum_number_of_tails) :: & 4068 tmp_tail 4069 4070 4071 new_maximum_number = maximum_number_of_tails + & 4072 MAX( 5*number_of_new_tails, number_of_initial_tails ) 4073 4074 IF ( write_particle_statistics ) THEN 4075 CALL check_open( 80 ) 4076 WRITE ( 80, '(''*** Request: '', I5, '' new_maximum_number(tails)'')' ) & 4077 new_maximum_number 4078 CALL close_file( 80 ) 4079 ENDIF 4080 WRITE (9,*) '*** Request: ',new_maximum_number,' new_maximum_number(tails)' 4081 ! CALL local_flush( 9 ) 4082 4083 tmp_tail(:,:,1:number_of_tails) = & 4084 particle_tail_coordinates(:,:,1:number_of_tails) 4085 tmp_tail_mask(1:number_of_tails) = tail_mask(1:number_of_tails) 4086 4087 DEALLOCATE( new_tail_id, particle_tail_coordinates, tail_mask ) 4088 ALLOCATE( new_tail_id(new_maximum_number), & 4089 particle_tail_coordinates(maximum_number_of_tailpoints,5, & 4090 new_maximum_number), & 4091 tail_mask(new_maximum_number) ) 4092 maximum_number_of_tails = new_maximum_number 4093 4094 particle_tail_coordinates = 0.0 4095 particle_tail_coordinates(:,:,1:number_of_tails) = & 4096 tmp_tail(:,:,1:number_of_tails) 4097 tail_mask(1:number_of_tails) = tmp_tail_mask(1:number_of_tails) 4098 tail_mask(number_of_tails+1:maximum_number_of_tails) = .TRUE. 4099 4100 END SUBROUTINE allocate_tail_memory 4101 4102 4103 SUBROUTINE output_particles_netcdf 4104 #if defined( __netcdf ) 4105 4106 USE control_parameters 4107 USE netcdf_control 4108 USE particle_attributes 4109 4110 IMPLICIT NONE 4111 4112 4113 CALL check_open( 108 ) 4114 4115 ! 4116 !-- Update the NetCDF time axis 4117 prt_time_count = prt_time_count + 1 4118 4119 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_time_prt, (/ simulated_time /), & 4120 start = (/ prt_time_count /), count = (/ 1 /) ) 4121 CALL handle_netcdf_error( 'output_particles_netcdf', 1 ) 4122 4123 ! 4124 !-- Output the real number of particles used 4125 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_rnop_prt, & 4126 (/ number_of_particles /), & 4127 start = (/ prt_time_count /), count = (/ 1 /) ) 4128 CALL handle_netcdf_error( 'output_particles_netcdf', 2 ) 4129 4130 ! 4131 !-- Output all particle attributes 4132 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(1), particles%age, & 4133 start = (/ 1, prt_time_count /), & 4134 count = (/ maximum_number_of_particles /) ) 4135 CALL handle_netcdf_error( 'output_particles_netcdf', 3 ) 4136 4137 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(2), particles%dvrp_psize, & 4138 start = (/ 1, prt_time_count /), & 4139 count = (/ maximum_number_of_particles /) ) 4140 CALL handle_netcdf_error( 'output_particles_netcdf', 4 ) 4141 4142 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(3), particles%origin_x, & 4143 start = (/ 1, prt_time_count /), & 4144 count = (/ maximum_number_of_particles /) ) 4145 CALL handle_netcdf_error( 'output_particles_netcdf', 5 ) 4146 4147 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(4), particles%origin_y, & 4148 start = (/ 1, prt_time_count /), & 4149 count = (/ maximum_number_of_particles /) ) 4150 CALL handle_netcdf_error( 'output_particles_netcdf', 6 ) 4151 4152 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(5), particles%origin_z, & 4153 start = (/ 1, prt_time_count /), & 4154 count = (/ maximum_number_of_particles /) ) 4155 CALL handle_netcdf_error( 'output_particles_netcdf', 7 ) 4156 4157 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(6), particles%radius, & 4158 start = (/ 1, prt_time_count /), & 4159 count = (/ maximum_number_of_particles /) ) 4160 CALL handle_netcdf_error( 'output_particles_netcdf', 8 ) 4161 4162 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(7), particles%speed_x, & 4163 start = (/ 1, prt_time_count /), & 4164 count = (/ maximum_number_of_particles /) ) 4165 CALL handle_netcdf_error( 'output_particles_netcdf', 9 ) 4166 4167 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(8), particles%speed_y, & 4168 start = (/ 1, prt_time_count /), & 4169 count = (/ maximum_number_of_particles /) ) 4170 CALL handle_netcdf_error( 'output_particles_netcdf', 10 ) 4171 4172 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(9), particles%speed_z, & 4173 start = (/ 1, prt_time_count /), & 4174 count = (/ maximum_number_of_particles /) ) 4175 CALL handle_netcdf_error( 'output_particles_netcdf', 11 ) 4176 4177 nc_stat = NF90_PUT_VAR( id_set_prt,id_var_prt(10),particles%weight_factor,& 4178 start = (/ 1, prt_time_count /), & 4179 count = (/ maximum_number_of_particles /) ) 4180 CALL handle_netcdf_error( 'output_particles_netcdf', 12 ) 4181 4182 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(11), particles%x, & 4183 start = (/ 1, prt_time_count /), & 4184 count = (/ maximum_number_of_particles /) ) 4185 CALL handle_netcdf_error( 'output_particles_netcdf', 13 ) 4186 4187 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(12), particles%y, & 4188 start = (/ 1, prt_time_count /), & 4189 count = (/ maximum_number_of_particles /) ) 4190 CALL handle_netcdf_error( 'output_particles_netcdf', 14 ) 4191 4192 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(13), particles%z, & 4193 start = (/ 1, prt_time_count /), & 4194 count = (/ maximum_number_of_particles /) ) 4195 CALL handle_netcdf_error( 'output_particles_netcdf', 15 ) 4196 4197 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(14), particles%class, & 4198 start = (/ 1, prt_time_count /), & 4199 count = (/ maximum_number_of_particles /) ) 4200 CALL handle_netcdf_error( 'output_particles_netcdf', 16 ) 4201 4202 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(15), particles%group, & 4203 start = (/ 1, prt_time_count /), & 4204 count = (/ maximum_number_of_particles /) ) 4205 CALL handle_netcdf_error( 'output_particles_netcdf', 17 ) 4206 4207 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(16), particles%tailpoints, & 4208 start = (/ 1, prt_time_count /), & 4209 count = (/ maximum_number_of_particles /) ) 4210 CALL handle_netcdf_error( 'output_particles_netcdf', 18 ) 4211 4212 nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(17), particles%tail_id, & 4213 start = (/ 1, prt_time_count /), & 4214 count = (/ maximum_number_of_particles /) ) 4215 CALL handle_netcdf_error( 'output_particles_netcdf', 19 ) 4216 4217 #endif 4218 END SUBROUTINE output_particles_netcdf 4219 4220 4221 SUBROUTINE write_particles 4222 4223 !------------------------------------------------------------------------------! 4224 ! Description: 4225 ! ------------ 4226 ! Write particle data on restart file 4227 !------------------------------------------------------------------------------! 4228 4229 USE control_parameters 4230 USE particle_attributes 4231 USE pegrid 4232 4233 IMPLICIT NONE 4234 4235 CHARACTER (LEN=10) :: particle_binary_version 4236 INTEGER :: i 4237 4238 ! 4239 !-- First open the output unit. 4240 IF ( myid_char == '' ) THEN 4241 OPEN ( 90, FILE='PARTICLE_RESTART_DATA_OUT'//myid_char, & 4242 FORM='UNFORMATTED') 4243 ELSE 4244 IF ( myid == 0 ) CALL local_system( 'mkdir PARTICLE_RESTART_DATA_OUT' ) 4245 #if defined( __parallel ) 4246 ! 4247 !-- Set a barrier in order to allow that thereafter all other processors 4248 !-- in the directory created by PE0 can open their file 4249 CALL MPI_BARRIER( comm2d, ierr ) 4250 #endif 4251 OPEN ( 90, FILE='PARTICLE_RESTART_DATA_OUT/'//myid_char, & 4252 FORM='UNFORMATTED' ) 4253 ENDIF 4254 4255 DO i = 0, io_blocks-1 4256 4257 IF ( i == io_group ) THEN 4258 4259 ! 4260 !-- Write the version number of the binary format. 4261 !-- Attention: After changes to the following output commands the version 4262 !-- --------- number of the variable particle_binary_version must be 4263 !-- changed! Also, the version number and the list of arrays 4264 !-- to be read in init_particles must be adjusted accordingly. 4265 particle_binary_version = '3.0' 4266 WRITE ( 90 ) particle_binary_version 4267 4268 ! 4269 !-- Write some particle parameters, the size of the particle arrays as 4270 !-- well as other dvrp-plot variables. 4271 WRITE ( 90 ) bc_par_b, bc_par_lr, bc_par_ns, bc_par_t, & 4272 maximum_number_of_particles, & 4273 maximum_number_of_tailpoints, maximum_number_of_tails, & 4274 number_of_initial_particles, number_of_particles, & 4275 number_of_particle_groups, number_of_tails, & 4276 particle_groups, time_prel, time_write_particle_data, & 4277 uniform_particles 4278 4279 IF ( number_of_initial_particles /= 0 ) WRITE ( 90 ) initial_particles 4280 4281 WRITE ( 90 ) prt_count, prt_start_index 4282 WRITE ( 90 ) particles 4283 4284 IF ( use_particle_tails ) THEN 4285 WRITE ( 90 ) particle_tail_coordinates 4286 ENDIF 4287 4288 CLOSE ( 90 ) 4289 4290 ENDIF 4291 4292 #if defined( __parallel ) 4293 CALL MPI_BARRIER( comm2d, ierr ) 4294 #endif 4295 4296 ENDDO 4297 4298 END SUBROUTINE write_particles 4299 4300 4301 4302 SUBROUTINE collision_efficiency( mean_r, r, e) 4303 !------------------------------------------------------------------------------! 4304 ! Description: 4305 ! ------------ 4306 ! Interpolate collision efficiency from table 4307 !------------------------------------------------------------------------------! 4308 4309 IMPLICIT NONE 4310 4311 INTEGER :: i, j, k 4312 4313 LOGICAL, SAVE :: first = .TRUE. 4314 4315 REAL :: aa, bb, cc, dd, dx, dy, e, gg, mean_r, mean_rm, r, rm, & 4316 x, y 4317 4318 REAL, DIMENSION(1:9), SAVE :: collected_r = 0.0 4319 REAL, DIMENSION(1:19), SAVE :: collector_r = 0.0 4320 REAL, DIMENSION(1:9,1:19), SAVE :: ef = 0.0 4321 4322 mean_rm = mean_r * 1.0E06 4323 rm = r * 1.0E06 4324 4325 IF ( first ) THEN 4326 collected_r = (/ 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0, 25.0 /) 4327 collector_r = (/ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 80.0, 100.0, 150.0,& 4328 200.0, 300.0, 400.0, 500.0, 600.0, 1000.0, 1400.0, & 4329 1800.0, 2400.0, 3000.0 /) 4330 ef(:,1) = (/0.017, 0.027, 0.037, 0.052, 0.052, 0.052, 0.052, 0.0, 0.0 /) 4331 ef(:,2) = (/0.001, 0.016, 0.027, 0.060, 0.12, 0.17, 0.17, 0.17, 0.0 /) 4332 ef(:,3) = (/0.001, 0.001, 0.02, 0.13, 0.28, 0.37, 0.54, 0.55, 0.47/) 4333 ef(:,4) = (/0.001, 0.001, 0.02, 0.23, 0.4, 0.55, 0.7, 0.75, 0.75/) 4334 ef(:,5) = (/0.01, 0.01, 0.03, 0.3, 0.4, 0.58, 0.73, 0.75, 0.79/) 4335 ef(:,6) = (/0.01, 0.01, 0.13, 0.38, 0.57, 0.68, 0.80, 0.86, 0.91/) 4336 ef(:,7) = (/0.01, 0.085, 0.23, 0.52, 0.68, 0.76, 0.86, 0.92, 0.95/) 4337 ef(:,8) = (/0.01, 0.14, 0.32, 0.60, 0.73, 0.81, 0.90, 0.94, 0.96/) 4338 ef(:,9) = (/0.025, 0.25, 0.43, 0.66, 0.78, 0.83, 0.92, 0.95, 0.96/) 4339 ef(:,10)= (/0.039, 0.3, 0.46, 0.69, 0.81, 0.87, 0.93, 0.95, 0.96/) 4340 ef(:,11)= (/0.095, 0.33, 0.51, 0.72, 0.82, 0.87, 0.93, 0.96, 0.97/) 4341 ef(:,12)= (/0.098, 0.36, 0.51, 0.73, 0.83, 0.88, 0.93, 0.96, 0.97/) 4342 ef(:,13)= (/0.1, 0.36, 0.52, 0.74, 0.83, 0.88, 0.93, 0.96, 0.97/) 4343 ef(:,14)= (/0.17, 0.4, 0.54, 0.72, 0.83, 0.88, 0.94, 0.98, 1.0 /) 4344 ef(:,15)= (/0.15, 0.37, 0.52, 0.74, 0.82, 0.88, 0.94, 0.98, 1.0 /) 4345 ef(:,16)= (/0.11, 0.34, 0.49, 0.71, 0.83, 0.88, 0.94, 0.95, 1.0 /) 4346 ef(:,17)= (/0.08, 0.29, 0.45, 0.68, 0.8, 0.86, 0.96, 0.94, 1.0 /) 4347 ef(:,18)= (/0.04, 0.22, 0.39, 0.62, 0.75, 0.83, 0.92, 0.96, 1.0 /) 4348 ef(:,19)= (/0.02, 0.16, 0.33, 0.55, 0.71, 0.81, 0.90, 0.94, 1.0 /) 4349 ENDIF 4350 4351 DO k = 1, 8 4352 IF ( collected_r(k) <= mean_rm ) i = k 4353 ENDDO 4354 4355 DO k = 1, 18 4356 IF ( collector_r(k) <= rm ) j = k 4357 ENDDO 4358 4359 IF ( rm < 10.0 ) THEN 4360 e = 0.0 4361 ELSEIF ( mean_rm < 2.0 ) THEN 4362 e = 0.001 4363 ELSEIF ( mean_rm >= 25.0 ) THEN 4364 IF( j <= 2 ) e = 0.0 4365 IF( j == 3 ) e = 0.47 4366 IF( j == 4 ) e = 0.8 4367 IF( j == 5 ) e = 0.9 4368 IF( j >=6 ) e = 1.0 4369 ELSEIF ( rm >= 3000.0 ) THEN 4370 IF( i == 1 ) e = 0.02 4371 IF( i == 2 ) e = 0.16 4372 IF( i == 3 ) e = 0.33 4373 IF( i == 4 ) e = 0.55 4374 IF( i == 5 ) e = 0.71 4375 IF( i == 6 ) e = 0.81 4376 IF( i == 7 ) e = 0.90 4377 IF( i >= 8 ) e = 0.94 4378 ELSE 4379 x = mean_rm - collected_r(i) 4380 y = rm - collector_r(j) 4381 dx = collected_r(i+1) - collected_r(i) 4382 dy = collector_r(j+1) - collector_r(j) 4383 aa = x**2 + y**2 4384 bb = ( dx - x )**2 + y**2 4385 cc = x**2 + ( dy - y )**2 4386 dd = ( dx - x )**2 + ( dy - y )**2 4387 gg = aa + bb + cc + dd 4388 4389 e = ( (gg-aa)*ef(i,j) + (gg-bb)*ef(i+1,j) + (gg-cc)*ef(i,j+1) + & 4390 (gg-dd)*ef(i+1,j+1) ) / (3.0*gg) 4391 ENDIF 4392 4393 END SUBROUTINE collision_efficiency 4394 4395 4396 4397 SUBROUTINE sort_particles 4398 4399 !------------------------------------------------------------------------------! 4400 ! Description: 4401 ! ------------ 4402 ! Sort particles in the sequence the grid boxes are stored in memory 4403 !------------------------------------------------------------------------------! 4404 4405 USE arrays_3d 4406 USE control_parameters 4407 USE cpulog 4408 USE grid_variables 4409 USE indices 4410 USE interfaces 4411 USE particle_attributes 4412 4413 IMPLICIT NONE 4414 4415 INTEGER :: i, ilow, j, k, n 4416 4417 TYPE(particle_type), DIMENSION(:), POINTER :: particles_temp 4418 4419 4420 CALL cpu_log( log_point_s(47), 'sort_particles', 'start' ) 4421 4422 ! 4423 !-- Initialize counters and set pointer of the temporary array into which 4424 !-- particles are sorted to free memory 4425 prt_count = 0 4426 sort_count = sort_count +1 4427 4428 SELECT CASE ( MOD( sort_count, 2 ) ) 4429 4430 CASE ( 0 ) 4431 4432 particles_temp => part_1 4433 ! part_1 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 4434 ! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 4435 ! 0.0, 0, 0, 0, 0 ) 4436 4437 CASE ( 1 ) 4438 4439 particles_temp => part_2 4440 ! part_2 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 4441 ! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 4442 ! 0.0, 0, 0, 0, 0 ) 4443 4444 END SELECT 4445 4446 ! 4447 !-- Count the particles per gridbox 4448 DO n = 1, number_of_particles 4449 4450 i = ( particles(n)%x + 0.5 * dx ) * ddx 4451 j = ( particles(n)%y + 0.5 * dy ) * ddy 4452 k = particles(n)%z / dz + 1 + offset_ocean_nzt 4453 ! only exact if equidistant 4454 4455 prt_count(k,j,i) = prt_count(k,j,i) + 1 4456 4457 IF ( i < nxl .OR. i > nxr .OR. j < nys .OR. j > nyn .OR. k < nzb+1 .OR. & 4458 k > nzt ) THEN 4459 WRITE( message_string, * ) ' particle out of range: i=', i, ' j=', & 4460 j, ' k=', k, & 4461 ' nxl=', nxl, ' nxr=', nxr, & 4462 ' nys=', nys, ' nyn=', nyn, & 4463 ' nzb=', nzb, ' nzt=', nzt 4464 CALL message( 'sort_particles', 'PA0149', 1, 2, 0, 6, 0 ) 4465 ENDIF 4466 4467 ENDDO 4468 4469 ! 4470 !-- Calculate the lower indices of those ranges of the particles-array 4471 !-- containing particles which belong to the same gridpox i,j,k 4472 ilow = 1 4473 DO i = nxl, nxr 4474 DO j = nys, nyn 4475 DO k = nzb+1, nzt 4476 prt_start_index(k,j,i) = ilow 4477 ilow = ilow + prt_count(k,j,i) 4478 ENDDO 4479 ENDDO 4480 ENDDO 4481 4482 ! 4483 !-- Sorting the particles 4484 DO n = 1, number_of_particles 4485 4486 i = ( particles(n)%x + 0.5 * dx ) * ddx 4487 j = ( particles(n)%y + 0.5 * dy ) * ddy 4488 k = particles(n)%z / dz + 1 + offset_ocean_nzt 4489 ! only exact if equidistant 4490 4491 particles_temp(prt_start_index(k,j,i)) = particles(n) 4492 4493 prt_start_index(k,j,i) = prt_start_index(k,j,i) + 1 4494 4495 ENDDO 4496 4497 ! 4498 !-- Redirect the particles pointer to the sorted array 4499 SELECT CASE ( MOD( sort_count, 2 ) ) 4500 4501 CASE ( 0 ) 4502 4503 particles => part_1 4504 4505 CASE ( 1 ) 4506 4507 particles => part_2 4508 4509 END SELECT 4510 4511 ! 4512 !-- Reset the index array to the actual start position 4513 DO i = nxl, nxr 4514 DO j = nys, nyn 4515 DO k = nzb+1, nzt 4516 prt_start_index(k,j,i) = prt_start_index(k,j,i) - prt_count(k,j,i) 4517 ENDDO 4518 ENDDO 4519 ENDDO 4520 4521 CALL cpu_log( log_point_s(47), 'sort_particles', 'stop' ) 4522 4523 END SUBROUTINE sort_particles 267 CALL user_lpm_set_attributes 268 269 270 ! 271 !-- If required, add the current particle positions to the particle tails 272 IF ( use_particle_tails ) CALL lpm_extend_tails 273 274 275 ! 276 !-- Write particle statistics (in particular the number of particles 277 !-- exchanged between the subdomains) on file 278 IF ( write_particle_statistics ) CALL lpm_write_exchange_statistics 279 280 CALL cpu_log( log_point(25), 'lpm', 'stop' ) 281 282 283 END SUBROUTINE lpm
Note: See TracChangeset
for help on using the changeset viewer.
