Changeset 4143 for palm/trunk
- Timestamp:
- Aug 5, 2019 3:14:53 PM (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lagrangian_particle_model_mod.f90
r4122 r4143 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Rename variable and change select case to if statement 28 ! 29 ! 4122 2019-07-26 13:11:56Z schwenkel 27 30 ! Implement reset method as bottom boundary condition 28 31 ! … … 278 281 CHARACTER(LEN=5) :: splitting_mode = 'const' !< splitting mode 279 282 280 CHARACTER(LEN=25) :: particle_ interpolation = 'trilinear'!< interpolation method for calculatin the particle283 CHARACTER(LEN=25) :: particle_advection_interpolation = 'trilinear' !< interpolation method for calculatin the particle 281 284 282 285 INTEGER(iwp) :: deleted_particles = 0 !< number of deleted particles per time step … … 315 318 LOGICAL :: use_kernel_tables = .FALSE. !< parameter, which turns on the use of precalculated collision kernels 316 319 LOGICAL :: write_particle_statistics = .FALSE. !< namelist parameter (see documentation) 320 LOGICAL :: interpolation_simple_predictor = .FALSE. !< flag for simple particle advection interpolation with predictor step 321 LOGICAL :: interpolation_simple_corrector = .FALSE. !< flag for simple particle advection interpolation with corrector step 322 LOGICAL :: interpolation_trilinear = .FALSE. !< flag for trilinear particle advection interpolation 317 323 318 324 LOGICAL, DIMENSION(max_number_of_particle_groups) :: vertical_particle_advection = .TRUE. !< Switch for vertical particle transport … … 555 561 particles_per_point, & 556 562 particle_advection_start, & 557 particle_ interpolation, &563 particle_advection_interpolation, & 558 564 particle_maximum_age, & 559 565 pdx, & … … 615 621 particles_per_point, & 616 622 particle_advection_start, & 617 particle_ interpolation, &623 particle_advection_interpolation, & 618 624 particle_maximum_age, & 619 625 pdx, & … … 862 868 !-- it can be combined as the sgs-velocites for active particles are 863 869 !-- calculated differently, i.e. no subboxes are needed. 864 IF ( .NOT. TRIM(particle_ interpolation) == 'trilinear' .AND. &870 IF ( .NOT. TRIM(particle_advection_interpolation) == 'trilinear' .AND. & 865 871 use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN 866 872 message_string = 'subrgrid scale velocities in combination with ' // & … … 1012 1018 de_dx = 0.0_wp 1013 1019 de_dy = 0.0_wp 1014 de_dz = 0.0_wp 1015 1020 de_dz = 0.0_wp 1021 1016 1022 sgs_wf_part = 1.0_wp / 3.0_wp 1017 1023 ENDIF … … 1063 1069 ENDDO 1064 1070 1071 ENDIF 1072 1073 ! 1074 !-- Check which particle interpolation method should be used 1075 IF ( TRIM(particle_advection_interpolation) == 'trilinear' ) THEN 1076 interpolation_simple_corrector = .FALSE. 1077 interpolation_simple_predictor = .FALSE. 1078 interpolation_trilinear = .TRUE. 1079 ELSEIF ( TRIM(particle_advection_interpolation) == 'simple_corrector' ) THEN 1080 interpolation_simple_corrector = .TRUE. 1081 interpolation_simple_predictor = .FALSE. 1082 interpolation_trilinear = .FALSE. 1083 ELSEIF ( TRIM(particle_advection_interpolation) == 'simple_predictor' ) THEN 1084 interpolation_simple_corrector = .FALSE. 1085 interpolation_simple_predictor = .TRUE. 1086 interpolation_trilinear = .FALSE. 1065 1087 ENDIF 1066 1088 … … 1309 1331 ! 1310 1332 !-- Calculate initial_weighting_factor diagnostically 1311 IF ( number_concentration /= -1.0_wp .AND. number_concentration > 0.0_wp ) THEN1312 initial_weighting_factor = number_concentration * &1333 IF ( number_concentration /= -1.0_wp .AND. number_concentration > 0.0_wp ) THEN 1334 initial_weighting_factor = number_concentration * & 1313 1335 pdx(1) * pdy(1) * pdz(1) 1314 1336 END IF … … 1322 1344 DO WHILE ( pos_y <= psn(i) ) 1323 1345 IF ( pos_y >= nys * dy .AND. & 1324 pos_y < ( nyn + 1 ) * dy ) THEN1346 pos_y < ( nyn + 1 ) * dy ) THEN 1325 1347 pos_x = psl(i) 1326 1348 xloop: DO WHILE ( pos_x <= psr(i) ) 1327 1349 IF ( pos_x >= nxl * dx .AND. & 1328 pos_x < ( nxr + 1) * dx ) THEN1350 pos_x < ( nxr + 1) * dx ) THEN 1329 1351 DO j = 1, particles_per_point 1330 1352 n = n + 1 … … 1366 1388 !-- In case of stretching the actual k index is found iteratively 1367 1389 IF ( dz_stretch_level .NE. -9999999.9_wp .OR. & 1368 dz_stretch_level_start(1) .NE. -9999999.9_wp ) THEN1390 dz_stretch_level_start(1) .NE. -9999999.9_wp ) THEN 1369 1391 kp = MINLOC( ABS( tmp_particle%z - zu ), DIM = 1 ) - 1 1370 1392 ELSE … … 1626 1648 ! 1627 1649 !-- Set constants for different aerosol species 1628 IF ( TRIM(aero_species) .EQ. 'nacl' ) THEN1650 IF ( TRIM(aero_species) .EQ. 'nacl' ) THEN 1629 1651 molecular_weight_of_solute = 0.05844_wp 1630 1652 rho_s = 2165.0_wp 1631 1653 vanthoff = 2.0_wp 1632 ELSEIF ( TRIM(aero_species) .EQ. 'c3h4o4' ) THEN1654 ELSEIF ( TRIM(aero_species) .EQ. 'c3h4o4' ) THEN 1633 1655 molecular_weight_of_solute = 0.10406_wp 1634 1656 rho_s = 1600.0_wp 1635 1657 vanthoff = 1.37_wp 1636 ELSEIF ( TRIM(aero_species) .EQ. 'nh4o3' ) THEN1658 ELSEIF ( TRIM(aero_species) .EQ. 'nh4o3' ) THEN 1637 1659 molecular_weight_of_solute = 0.08004_wp 1638 1660 rho_s = 1720.0_wp … … 2153 2175 2154 2176 IF ( .NOT. first_loop_stride .AND. & 2155 grid_particles(k,j,i)%time_loop_done ) CYCLE2177 grid_particles(k,j,i)%time_loop_done ) CYCLE 2156 2178 2157 2179 particles => grid_particles(k,j,i)%particles(1:number_of_particles) … … 2187 2209 ! 2188 2210 !-- Particle advection 2189 CALL lpm_advec( TRIM(particle_interpolation),i,j,k)2211 CALL lpm_advec(i,j,k) 2190 2212 ! 2191 2213 !-- Particle reflection from walls. Only applied if the particles … … 2243 2265 !-- Apply splitting and merging algorithm 2244 2266 IF ( cloud_droplets ) THEN 2245 IF ( splitting ) THEN2267 IF ( splitting ) THEN 2246 2268 CALL lpm_splitting 2247 2269 ENDIF 2248 IF ( merging ) THEN2270 IF ( merging ) THEN 2249 2271 CALL lpm_merging 2250 2272 ENDIF … … 2392 2414 CALL close_file( 80 ) 2393 2415 2394 IF ( number_of_particles > 0 ) THEN2416 IF ( number_of_particles > 0 ) THEN 2395 2417 WRITE(9,*) 'number_of_particles ', number_of_particles, & 2396 2418 current_timestep_number + 1, simulated_time + dt_3d … … 3131 3153 WRITE ( 14 ) curvature_solution_effects 3132 3154 3155 CALL wrd_write_string( 'interpolation_simple_corrector' ) 3156 WRITE ( 14 ) interpolation_simple_corrector 3157 3158 CALL wrd_write_string( 'interpolation_simple_predictor' ) 3159 WRITE ( 14 ) interpolation_simple_predictor 3160 3161 CALL wrd_write_string( 'interpolation_trilinear' ) 3162 WRITE ( 14 ) interpolation_trilinear 3163 3133 3164 END SUBROUTINE lpm_wrd_global 3134 3165 … … 3152 3183 CASE ( 'curvature_solution_effects' ) 3153 3184 READ ( 13 ) curvature_solution_effects 3154 3185 3186 CASE ( 'interpolation_simple_corrector' ) 3187 READ ( 13 ) interpolation_simple_corrector 3188 3189 CASE ( 'interpolation_simple_predictor' ) 3190 READ ( 13 ) interpolation_simple_predictor 3191 3192 CASE ( 'interpolation_trilinear' ) 3193 READ ( 13 ) interpolation_trilinear 3194 3155 3195 ! CASE ( 'global_paramter' ) 3156 3196 ! READ ( 13 ) global_parameter … … 3176 3216 !> this submodule should be excluded as an own file. 3177 3217 !------------------------------------------------------------------------------! 3178 SUBROUTINE lpm_advec (interpolation_method,ip,jp,kp) 3179 3180 CHARACTER (LEN=*), INTENT(IN) :: interpolation_method !< 3218 SUBROUTINE lpm_advec (ip,jp,kp) 3219 3181 3220 LOGICAL :: subbox_at_wall !< flag to see if the current subgridbox is adjacent to a wall 3182 3221 … … 3307 3346 dt_particle = dt_3d 3308 3347 3309 3310 SELECT CASE ( interpolation_method ) 3311 3312 ! 3313 !-- This case uses a simple interpolation method for the particle velocites, 3314 !-- and applying a predictor-corrector method. @attention: for the corrector 3315 !-- step the velocities of t(n+1) are required. However, at this moment of 3316 !-- the time integration they are not free of divergence. This interpolation 3317 !-- method is described in more detail in Grabowski et al., 2018 (GMD). 3318 CASE ( 'simple_corrector' ) 3319 ! 3320 !-- Predictor step 3321 kkw = kp - 1 3322 DO n = 1, number_of_particles 3323 3324 alpha = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp ) 3325 u_int(n) = u(kp,jp,ip) * ( 1.0_wp - alpha ) + u(kp,jp,ip+1) * alpha 3326 3327 beta = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp ) 3328 v_int(n) = v(kp,jp,ip) * ( 1.0_wp - beta ) + v(kp,jp+1,ip) * beta 3329 3330 gamma = MAX( MIN( ( particles(n)%z - zw(kkw) ) / & 3348 ! 3349 !-- This case uses a simple interpolation method for the particle velocites, 3350 !-- and applying a predictor-corrector method. @attention: for the corrector 3351 !-- step the velocities of t(n+1) are required. However, at this moment of 3352 !-- the time integration they are not free of divergence. This interpolation 3353 !-- method is described in more detail in Grabowski et al., 2018 (GMD). 3354 IF ( interpolation_simple_corrector ) THEN 3355 ! 3356 !-- Predictor step 3357 kkw = kp - 1 3358 DO n = 1, number_of_particles 3359 3360 alpha = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp ) 3361 u_int(n) = u(kp,jp,ip) * ( 1.0_wp - alpha ) + u(kp,jp,ip+1) * alpha 3362 3363 beta = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp ) 3364 v_int(n) = v(kp,jp,ip) * ( 1.0_wp - beta ) + v(kp,jp+1,ip) * beta 3365 3366 gamma = MAX( MIN( ( particles(n)%z - zw(kkw) ) / & 3367 ( zw(kkw+1) - zw(kkw) ), 1.0_wp ), 0.0_wp ) 3368 w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma 3369 3370 ENDDO 3371 ! 3372 !-- Corrector step 3373 DO n = 1, number_of_particles 3374 3375 IF ( .NOT. particles(n)%particle_mask ) CYCLE 3376 3377 DO nn = 1, iteration_steps 3378 3379 ! 3380 !-- Guess new position 3381 xp = particles(n)%x + u_int(n) * dt_particle(n) 3382 yp = particles(n)%y + v_int(n) * dt_particle(n) 3383 zp = particles(n)%z + w_int(n) * dt_particle(n) 3384 ! 3385 !-- x direction 3386 i_next = FLOOR( xp * ddx , KIND=iwp) 3387 alpha = MAX( MIN( ( xp - i_next * dx ) * ddx, 1.0_wp ), 0.0_wp ) 3388 ! 3389 !-- y direction 3390 j_next = FLOOR( yp * ddy ) 3391 beta = MAX( MIN( ( yp - j_next * dy ) * ddy, 1.0_wp ), 0.0_wp ) 3392 ! 3393 !-- z_direction 3394 k_next = MAX( MIN( FLOOR( zp / (zw(kkw+1)-zw(kkw)) ), nzt ), 0) 3395 gamma = MAX( MIN( ( zp - zw(k_next) ) / & 3396 ( zw(k_next+1) - zw(k_next) ), 1.0_wp ), 0.0_wp ) 3397 ! 3398 !-- Calculate part of the corrector step 3399 unext = u_p(k_next+1, j_next, i_next) * ( 1.0_wp - alpha ) + & 3400 u_p(k_next+1, j_next, i_next+1) * alpha 3401 3402 vnext = v_p(k_next+1, j_next, i_next) * ( 1.0_wp - beta ) + & 3403 v_p(k_next+1, j_next+1, i_next ) * beta 3404 3405 wnext = w_p(k_next, j_next, i_next) * ( 1.0_wp - gamma ) + & 3406 w_p(k_next+1, j_next, i_next ) * gamma 3407 3408 ! 3409 !-- Calculate interpolated particle velocity with predictor 3410 !-- corrector step. u_int, v_int and w_int describes the part of 3411 !-- the predictor step. unext, vnext and wnext is the part of the 3412 !-- corrector step. The resulting new position is set below. The 3413 !-- implementation is based on Grabowski et al., 2018 (GMD). 3414 u_int(n) = 0.5_wp * ( u_int(n) + unext ) 3415 v_int(n) = 0.5_wp * ( v_int(n) + vnext ) 3416 w_int(n) = 0.5_wp * ( w_int(n) + wnext ) 3417 3418 ENDDO 3419 ENDDO 3420 ! 3421 !-- This case uses a simple interpolation method for the particle velocites, 3422 !-- and applying a predictor. 3423 ELSEIF ( interpolation_simple_predictor ) THEN 3424 ! 3425 !-- The particle position for the w velociy is based on the value of kp and kp-1 3426 kkw = kp - 1 3427 DO n = 1, number_of_particles 3428 IF ( .NOT. particles(n)%particle_mask ) CYCLE 3429 3430 alpha = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp ) 3431 u_int(n) = u(kp,jp,ip) * ( 1.0_wp - alpha ) + u(kp,jp,ip+1) * alpha 3432 3433 beta = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp ) 3434 v_int(n) = v(kp,jp,ip) * ( 1.0_wp - beta ) + v(kp,jp+1,ip) * beta 3435 3436 gamma = MAX( MIN( ( particles(n)%z - zw(kkw) ) / & 3331 3437 ( zw(kkw+1) - zw(kkw) ), 1.0_wp ), 0.0_wp ) 3332 w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma 3333 3438 w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma 3439 ENDDO 3440 ! 3441 !-- The trilinear interpolation. 3442 ELSEIF ( interpolation_trilinear ) THEN 3443 3444 start_index = grid_particles(kp,jp,ip)%start_index 3445 end_index = grid_particles(kp,jp,ip)%end_index 3446 3447 DO nb = 0, 7 3448 ! 3449 !-- Interpolate u velocity-component 3450 i = ip 3451 j = jp + block_offset(nb)%j_off 3452 k = kp + block_offset(nb)%k_off 3453 3454 DO n = start_index(nb), end_index(nb) 3455 ! 3456 !-- Interpolation of the u velocity component onto particle position. 3457 !-- Particles are interpolation bi-linearly in the horizontal and a 3458 !-- linearly in the vertical. An exception is made for particles below 3459 !-- the first vertical grid level in case of a prandtl layer. In this 3460 !-- case the horizontal particle velocity components are determined using 3461 !-- Monin-Obukhov relations (if branch). 3462 !-- First, check if particle is located below first vertical grid level 3463 !-- above topography (Prandtl-layer height) 3464 !-- Determine vertical index of topography top 3465 k_wall = get_topography_top_index_ji( jp, ip, 's' ) 3466 3467 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN 3468 ! 3469 !-- Resolved-scale horizontal particle velocity is zero below z0. 3470 IF ( zv(n) - zw(k_wall) < z0_av_global ) THEN 3471 u_int(n) = 0.0_wp 3472 ELSE 3473 ! 3474 !-- Determine the sublayer. Further used as index. 3475 height_p = ( zv(n) - zw(k_wall) - z0_av_global ) & 3476 * REAL( number_of_sublayers, KIND=wp ) & 3477 * d_z_p_z0 3478 ! 3479 !-- Calculate LOG(z/z0) for exact particle height. Therefore, 3480 !-- interpolate linearly between precalculated logarithm. 3481 log_z_z0_int = log_z_z0(INT(height_p)) & 3482 + ( height_p - INT(height_p) ) & 3483 * ( log_z_z0(INT(height_p)+1) & 3484 - log_z_z0(INT(height_p)) & 3485 ) 3486 ! 3487 !-- Get friction velocity and momentum flux from new surface data 3488 !-- types. 3489 IF ( surf_def_h(0)%start_index(jp,ip) <= & 3490 surf_def_h(0)%end_index(jp,ip) ) THEN 3491 surf_start = surf_def_h(0)%start_index(jp,ip) 3492 !-- Limit friction velocity. In narrow canyons or holes the 3493 !-- friction velocity can become very small, resulting in a too 3494 !-- large particle speed. 3495 us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 3496 usws_int = surf_def_h(0)%usws(surf_start) 3497 ELSEIF ( surf_lsm_h%start_index(jp,ip) <= & 3498 surf_lsm_h%end_index(jp,ip) ) THEN 3499 surf_start = surf_lsm_h%start_index(jp,ip) 3500 us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 3501 usws_int = surf_lsm_h%usws(surf_start) 3502 ELSEIF ( surf_usm_h%start_index(jp,ip) <= & 3503 surf_usm_h%end_index(jp,ip) ) THEN 3504 surf_start = surf_usm_h%start_index(jp,ip) 3505 us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 3506 usws_int = surf_usm_h%usws(surf_start) 3507 ENDIF 3508 ! 3509 !-- Neutral solution is applied for all situations, e.g. also for 3510 !-- unstable and stable situations. Even though this is not exact 3511 !-- this saves a lot of CPU time since several calls of intrinsic 3512 !-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified 3513 !-- as sensitivity studies revealed no significant effect of 3514 !-- using the neutral solution also for un/stable situations. 3515 u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp ) & 3516 * log_z_z0_int - u_gtrans 3517 ENDIF 3518 ! 3519 !-- Particle above the first grid level. Bi-linear interpolation in the 3520 !-- horizontal and linear interpolation in the vertical direction. 3521 ELSE 3522 x = xv(n) - i * dx 3523 y = yv(n) + ( 0.5_wp - j ) * dy 3524 aa = x**2 + y**2 3525 bb = ( dx - x )**2 + y**2 3526 cc = x**2 + ( dy - y )**2 3527 dd = ( dx - x )**2 + ( dy - y )**2 3528 gg = aa + bb + cc + dd 3529 3530 u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) & 3531 + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * & 3532 u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans 3533 3534 IF ( k == nzt ) THEN 3535 u_int(n) = u_int_l 3536 ELSE 3537 u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) & 3538 + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * & 3539 u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans 3540 u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dzw(k+1) * & 3541 ( u_int_u - u_int_l ) 3542 ENDIF 3543 ENDIF 3334 3544 ENDDO 3335 3336 ! 3337 !-- Corrector step 3338 DO n = 1, number_of_particles 3339 3340 IF ( .NOT. particles(n)%particle_mask ) CYCLE 3341 3342 DO nn = 1, iteration_steps 3343 3344 ! 3345 !-- Guess new position 3346 xp = particles(n)%x + u_int(n) * dt_particle(n) 3347 yp = particles(n)%y + v_int(n) * dt_particle(n) 3348 zp = particles(n)%z + w_int(n) * dt_particle(n) 3349 ! 3350 !-- x direction 3351 i_next = FLOOR( xp * ddx , KIND=iwp) 3352 alpha = MAX( MIN( ( xp - i_next * dx ) * ddx, 1.0_wp ), 0.0_wp ) 3353 ! 3354 !-- y direction 3355 j_next = FLOOR( yp * ddy ) 3356 beta = MAX( MIN( ( yp - j_next * dy ) * ddy, 1.0_wp ), 0.0_wp ) 3357 ! 3358 !-- z_direction 3359 k_next = MAX( MIN( FLOOR( zp / (zw(kkw+1)-zw(kkw)) ), nzt ), 0) 3360 gamma = MAX( MIN( ( zp - zw(k_next) ) / & 3361 ( zw(k_next+1) - zw(k_next) ), 1.0_wp ), 0.0_wp ) 3362 ! 3363 !-- Calculate part of the corrector step 3364 unext = u_p(k_next+1, j_next, i_next) * ( 1.0_wp - alpha ) + & 3365 u_p(k_next+1, j_next, i_next+1) * alpha 3366 3367 vnext = v_p(k_next+1, j_next, i_next) * ( 1.0_wp - beta ) + & 3368 v_p(k_next+1, j_next+1, i_next ) * beta 3369 3370 wnext = w_p(k_next, j_next, i_next) * ( 1.0_wp - gamma ) + & 3371 w_p(k_next+1, j_next, i_next ) * gamma 3372 3373 ! 3374 !-- Calculate interpolated particle velocity with predictor 3375 !-- corrector step. u_int, v_int and w_int describes the part of 3376 !-- the predictor step. unext, vnext and wnext is the part of the 3377 !-- corrector step. The resulting new position is set below. The 3378 !-- implementation is based on Grabowski et al., 2018 (GMD). 3379 u_int(n) = 0.5_wp * ( u_int(n) + unext ) 3380 v_int(n) = 0.5_wp * ( v_int(n) + vnext ) 3381 w_int(n) = 0.5_wp * ( w_int(n) + wnext ) 3382 3383 ENDDO 3545 ! 3546 !-- Same procedure for interpolation of the v velocity-component 3547 i = ip + block_offset(nb)%i_off 3548 j = jp 3549 k = kp + block_offset(nb)%k_off 3550 3551 DO n = start_index(nb), end_index(nb) 3552 ! 3553 !-- Determine vertical index of topography top 3554 k_wall = get_topography_top_index_ji( jp,ip, 's' ) 3555 3556 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN 3557 IF ( zv(n) - zw(k_wall) < z0_av_global ) THEN 3558 ! 3559 !-- Resolved-scale horizontal particle velocity is zero below z0. 3560 v_int(n) = 0.0_wp 3561 ELSE 3562 ! 3563 !-- Determine the sublayer. Further used as index. Please note, 3564 !-- logarithmus can not be reused from above, as in in case of 3565 !-- topography particle on u-grid can be above surface-layer height, 3566 !-- whereas it can be below on v-grid. 3567 height_p = ( zv(n) - zw(k_wall) - z0_av_global ) & 3568 * REAL( number_of_sublayers, KIND=wp ) & 3569 * d_z_p_z0 3570 ! 3571 !-- Calculate LOG(z/z0) for exact particle height. Therefore, 3572 !-- interpolate linearly between precalculated logarithm. 3573 log_z_z0_int = log_z_z0(INT(height_p)) & 3574 + ( height_p - INT(height_p) ) & 3575 * ( log_z_z0(INT(height_p)+1) & 3576 - log_z_z0(INT(height_p)) & 3577 ) 3578 ! 3579 !-- Get friction velocity and momentum flux from new surface data 3580 !-- types. 3581 IF ( surf_def_h(0)%start_index(jp,ip) <= & 3582 surf_def_h(0)%end_index(jp,ip) ) THEN 3583 surf_start = surf_def_h(0)%start_index(jp,ip) 3584 !-- Limit friction velocity. In narrow canyons or holes the 3585 !-- friction velocity can become very small, resulting in a too 3586 !-- large particle speed. 3587 us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 3588 vsws_int = surf_def_h(0)%vsws(surf_start) 3589 ELSEIF ( surf_lsm_h%start_index(jp,ip) <= & 3590 surf_lsm_h%end_index(jp,ip) ) THEN 3591 surf_start = surf_lsm_h%start_index(jp,ip) 3592 us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 3593 vsws_int = surf_lsm_h%vsws(surf_start) 3594 ELSEIF ( surf_usm_h%start_index(jp,ip) <= & 3595 surf_usm_h%end_index(jp,ip) ) THEN 3596 surf_start = surf_usm_h%start_index(jp,ip) 3597 us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 3598 vsws_int = surf_usm_h%vsws(surf_start) 3599 ENDIF 3600 ! 3601 !-- Neutral solution is applied for all situations, e.g. also for 3602 !-- unstable and stable situations. Even though this is not exact 3603 !-- this saves a lot of CPU time since several calls of intrinsic 3604 !-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified 3605 !-- as sensitivity studies revealed no significant effect of 3606 !-- using the neutral solution also for un/stable situations. 3607 v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp ) & 3608 * log_z_z0_int - v_gtrans 3609 3610 ENDIF 3611 ELSE 3612 x = xv(n) + ( 0.5_wp - i ) * dx 3613 y = yv(n) - j * dy 3614 aa = x**2 + y**2 3615 bb = ( dx - x )**2 + y**2 3616 cc = x**2 + ( dy - y )**2 3617 dd = ( dx - x )**2 + ( dy - y )**2 3618 gg = aa + bb + cc + dd 3619 3620 v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) & 3621 + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) & 3622 ) / ( 3.0_wp * gg ) - v_gtrans 3623 3624 IF ( k == nzt ) THEN 3625 v_int(n) = v_int_l 3626 ELSE 3627 v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) & 3628 + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) & 3629 ) / ( 3.0_wp * gg ) - v_gtrans 3630 v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dzw(k+1) * & 3631 ( v_int_u - v_int_l ) 3632 ENDIF 3633 ENDIF 3384 3634 ENDDO 3385 3386 3387 ! 3388 !-- This case uses a simple interpolation method for the particle velocites, 3389 !-- and applying a predictor. 3390 CASE ( 'simple_predictor' ) 3391 ! 3392 !-- The particle position for the w velociy is based on the value of kp and kp-1 3393 kkw = kp - 1 3394 DO n = 1, number_of_particles 3395 IF ( .NOT. particles(n)%particle_mask ) CYCLE 3396 3397 alpha = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp ) 3398 u_int(n) = u(kp,jp,ip) * ( 1.0_wp - alpha ) + u(kp,jp,ip+1) * alpha 3399 3400 beta = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp ) 3401 v_int(n) = v(kp,jp,ip) * ( 1.0_wp - beta ) + v(kp,jp+1,ip) * beta 3402 3403 gamma = MAX( MIN( ( particles(n)%z - zw(kkw) ) / & 3404 ( zw(kkw+1) - zw(kkw) ), 1.0_wp ), 0.0_wp ) 3405 w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma 3635 ! 3636 !-- Same procedure for interpolation of the w velocity-component 3637 i = ip + block_offset(nb)%i_off 3638 j = jp + block_offset(nb)%j_off 3639 k = kp - 1 3640 3641 DO n = start_index(nb), end_index(nb) 3642 IF ( vertical_particle_advection(particles(n)%group) ) THEN 3643 x = xv(n) + ( 0.5_wp - i ) * dx 3644 y = yv(n) + ( 0.5_wp - j ) * dy 3645 aa = x**2 + y**2 3646 bb = ( dx - x )**2 + y**2 3647 cc = x**2 + ( dy - y )**2 3648 dd = ( dx - x )**2 + ( dy - y )**2 3649 gg = aa + bb + cc + dd 3650 3651 w_int_l = ( ( gg - aa ) * w(k,j,i) + ( gg - bb ) * w(k,j,i+1) & 3652 + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) & 3653 ) / ( 3.0_wp * gg ) 3654 3655 IF ( k == nzt ) THEN 3656 w_int(n) = w_int_l 3657 ELSE 3658 w_int_u = ( ( gg-aa ) * w(k+1,j,i) + & 3659 ( gg-bb ) * w(k+1,j,i+1) + & 3660 ( gg-cc ) * w(k+1,j+1,i) + & 3661 ( gg-dd ) * w(k+1,j+1,i+1) & 3662 ) / ( 3.0_wp * gg ) 3663 w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dzw(k+1) * & 3664 ( w_int_u - w_int_l ) 3665 ENDIF 3666 ELSE 3667 w_int(n) = 0.0_wp 3668 ENDIF 3406 3669 ENDDO 3407 ! 3408 !-- The trilinear interpolation. 3409 CASE ( 'trilinear' ) 3410 3411 start_index = grid_particles(kp,jp,ip)%start_index 3412 end_index = grid_particles(kp,jp,ip)%end_index 3413 3414 DO nb = 0, 7 3415 ! 3416 !-- Interpolate u velocity-component 3417 i = ip 3418 j = jp + block_offset(nb)%j_off 3419 k = kp + block_offset(nb)%k_off 3420 3421 DO n = start_index(nb), end_index(nb) 3422 ! 3423 !-- Interpolation of the u velocity component onto particle position. 3424 !-- Particles are interpolation bi-linearly in the horizontal and a 3425 !-- linearly in the vertical. An exception is made for particles below 3426 !-- the first vertical grid level in case of a prandtl layer. In this 3427 !-- case the horizontal particle velocity components are determined using 3428 !-- Monin-Obukhov relations (if branch). 3429 !-- First, check if particle is located below first vertical grid level 3430 !-- above topography (Prandtl-layer height) 3431 !-- Determine vertical index of topography top 3432 k_wall = get_topography_top_index_ji( jp, ip, 's' ) 3433 3434 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN 3435 ! 3436 !-- Resolved-scale horizontal particle velocity is zero below z0. 3437 IF ( zv(n) - zw(k_wall) < z0_av_global ) THEN 3438 u_int(n) = 0.0_wp 3439 ELSE 3440 ! 3441 !-- Determine the sublayer. Further used as index. 3442 height_p = ( zv(n) - zw(k_wall) - z0_av_global ) & 3443 * REAL( number_of_sublayers, KIND=wp ) & 3444 * d_z_p_z0 3445 ! 3446 !-- Calculate LOG(z/z0) for exact particle height. Therefore, 3447 !-- interpolate linearly between precalculated logarithm. 3448 log_z_z0_int = log_z_z0(INT(height_p)) & 3449 + ( height_p - INT(height_p) ) & 3450 * ( log_z_z0(INT(height_p)+1) & 3451 - log_z_z0(INT(height_p)) & 3452 ) 3453 ! 3454 !-- Get friction velocity and momentum flux from new surface data 3455 !-- types. 3456 IF ( surf_def_h(0)%start_index(jp,ip) <= & 3457 surf_def_h(0)%end_index(jp,ip) ) THEN 3458 surf_start = surf_def_h(0)%start_index(jp,ip) 3459 !-- Limit friction velocity. In narrow canyons or holes the 3460 !-- friction velocity can become very small, resulting in a too 3461 !-- large particle speed. 3462 us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 3463 usws_int = surf_def_h(0)%usws(surf_start) 3464 ELSEIF ( surf_lsm_h%start_index(jp,ip) <= & 3465 surf_lsm_h%end_index(jp,ip) ) THEN 3466 surf_start = surf_lsm_h%start_index(jp,ip) 3467 us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 3468 usws_int = surf_lsm_h%usws(surf_start) 3469 ELSEIF ( surf_usm_h%start_index(jp,ip) <= & 3470 surf_usm_h%end_index(jp,ip) ) THEN 3471 surf_start = surf_usm_h%start_index(jp,ip) 3472 us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 3473 usws_int = surf_usm_h%usws(surf_start) 3474 ENDIF 3475 ! 3476 !-- Neutral solution is applied for all situations, e.g. also for 3477 !-- unstable and stable situations. Even though this is not exact 3478 !-- this saves a lot of CPU time since several calls of intrinsic 3479 !-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified 3480 !-- as sensitivity studies revealed no significant effect of 3481 !-- using the neutral solution also for un/stable situations. 3482 u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp ) & 3483 * log_z_z0_int - u_gtrans 3484 ENDIF 3485 ! 3486 !-- Particle above the first grid level. Bi-linear interpolation in the 3487 !-- horizontal and linear interpolation in the vertical direction. 3488 ELSE 3489 x = xv(n) - i * dx 3490 y = yv(n) + ( 0.5_wp - j ) * dy 3491 aa = x**2 + y**2 3492 bb = ( dx - x )**2 + y**2 3493 cc = x**2 + ( dy - y )**2 3494 dd = ( dx - x )**2 + ( dy - y )**2 3495 gg = aa + bb + cc + dd 3496 3497 u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) & 3498 + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * & 3499 u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans 3500 3501 IF ( k == nzt ) THEN 3502 u_int(n) = u_int_l 3503 ELSE 3504 u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) & 3505 + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * & 3506 u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans 3507 u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dzw(k+1) * & 3508 ( u_int_u - u_int_l ) 3509 ENDIF 3510 ENDIF 3511 ENDDO 3512 ! 3513 !-- Same procedure for interpolation of the v velocity-component 3514 i = ip + block_offset(nb)%i_off 3515 j = jp 3516 k = kp + block_offset(nb)%k_off 3517 3518 DO n = start_index(nb), end_index(nb) 3519 ! 3520 !-- Determine vertical index of topography top 3521 k_wall = get_topography_top_index_ji( jp,ip, 's' ) 3522 3523 IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN 3524 IF ( zv(n) - zw(k_wall) < z0_av_global ) THEN 3525 ! 3526 !-- Resolved-scale horizontal particle velocity is zero below z0. 3527 v_int(n) = 0.0_wp 3528 ELSE 3529 ! 3530 !-- Determine the sublayer. Further used as index. Please note, 3531 !-- logarithmus can not be reused from above, as in in case of 3532 !-- topography particle on u-grid can be above surface-layer height, 3533 !-- whereas it can be below on v-grid. 3534 height_p = ( zv(n) - zw(k_wall) - z0_av_global ) & 3535 * REAL( number_of_sublayers, KIND=wp ) & 3536 * d_z_p_z0 3537 ! 3538 !-- Calculate LOG(z/z0) for exact particle height. Therefore, 3539 !-- interpolate linearly between precalculated logarithm. 3540 log_z_z0_int = log_z_z0(INT(height_p)) & 3541 + ( height_p - INT(height_p) ) & 3542 * ( log_z_z0(INT(height_p)+1) & 3543 - log_z_z0(INT(height_p)) & 3544 ) 3545 ! 3546 !-- Get friction velocity and momentum flux from new surface data 3547 !-- types. 3548 IF ( surf_def_h(0)%start_index(jp,ip) <= & 3549 surf_def_h(0)%end_index(jp,ip) ) THEN 3550 surf_start = surf_def_h(0)%start_index(jp,ip) 3551 !-- Limit friction velocity. In narrow canyons or holes the 3552 !-- friction velocity can become very small, resulting in a too 3553 !-- large particle speed. 3554 us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 3555 vsws_int = surf_def_h(0)%vsws(surf_start) 3556 ELSEIF ( surf_lsm_h%start_index(jp,ip) <= & 3557 surf_lsm_h%end_index(jp,ip) ) THEN 3558 surf_start = surf_lsm_h%start_index(jp,ip) 3559 us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 3560 vsws_int = surf_lsm_h%vsws(surf_start) 3561 ELSEIF ( surf_usm_h%start_index(jp,ip) <= & 3562 surf_usm_h%end_index(jp,ip) ) THEN 3563 surf_start = surf_usm_h%start_index(jp,ip) 3564 us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 3565 vsws_int = surf_usm_h%vsws(surf_start) 3566 ENDIF 3567 ! 3568 !-- Neutral solution is applied for all situations, e.g. also for 3569 !-- unstable and stable situations. Even though this is not exact 3570 !-- this saves a lot of CPU time since several calls of intrinsic 3571 !-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified 3572 !-- as sensitivity studies revealed no significant effect of 3573 !-- using the neutral solution also for un/stable situations. 3574 v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp ) & 3575 * log_z_z0_int - v_gtrans 3576 3577 ENDIF 3578 ELSE 3579 x = xv(n) + ( 0.5_wp - i ) * dx 3580 y = yv(n) - j * dy 3581 aa = x**2 + y**2 3582 bb = ( dx - x )**2 + y**2 3583 cc = x**2 + ( dy - y )**2 3584 dd = ( dx - x )**2 + ( dy - y )**2 3585 gg = aa + bb + cc + dd 3586 3587 v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) & 3588 + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) & 3589 ) / ( 3.0_wp * gg ) - v_gtrans 3590 3591 IF ( k == nzt ) THEN 3592 v_int(n) = v_int_l 3593 ELSE 3594 v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) & 3595 + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) & 3596 ) / ( 3.0_wp * gg ) - v_gtrans 3597 v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dzw(k+1) * & 3598 ( v_int_u - v_int_l ) 3599 ENDIF 3600 ENDIF 3601 ENDDO 3602 ! 3603 !-- Same procedure for interpolation of the w velocity-component 3604 i = ip + block_offset(nb)%i_off 3605 j = jp + block_offset(nb)%j_off 3606 k = kp - 1 3607 3608 DO n = start_index(nb), end_index(nb) 3609 IF ( vertical_particle_advection(particles(n)%group) ) THEN 3610 x = xv(n) + ( 0.5_wp - i ) * dx 3611 y = yv(n) + ( 0.5_wp - j ) * dy 3612 aa = x**2 + y**2 3613 bb = ( dx - x )**2 + y**2 3614 cc = x**2 + ( dy - y )**2 3615 dd = ( dx - x )**2 + ( dy - y )**2 3616 gg = aa + bb + cc + dd 3617 3618 w_int_l = ( ( gg - aa ) * w(k,j,i) + ( gg - bb ) * w(k,j,i+1) & 3619 + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) & 3620 ) / ( 3.0_wp * gg ) 3621 3622 IF ( k == nzt ) THEN 3623 w_int(n) = w_int_l 3624 ELSE 3625 w_int_u = ( ( gg-aa ) * w(k+1,j,i) + & 3626 ( gg-bb ) * w(k+1,j,i+1) + & 3627 ( gg-cc ) * w(k+1,j+1,i) + & 3628 ( gg-dd ) * w(k+1,j+1,i+1) & 3629 ) / ( 3.0_wp * gg ) 3630 w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dzw(k+1) * & 3631 ( w_int_u - w_int_l ) 3632 ENDIF 3633 ELSE 3634 w_int(n) = 0.0_wp 3635 ENDIF 3636 ENDDO 3637 ENDDO 3638 3639 CASE DEFAULT 3640 WRITE( message_string, * ) 'unknown particle velocity interpolation method = "', & 3641 TRIM( interpolation_method ), '"' 3642 CALL message( 'lpm_advec', 'PA0660', 1, 2, 0, 6, 0 ) 3643 3644 END SELECT 3670 ENDDO 3671 ENDIF 3645 3672 3646 3673 !-- Interpolate and calculate quantities needed for calculating the SGS … … 3653 3680 ! 3654 3681 !-- In case of topography check if subbox is adjacent to a wall 3655 IF ( .NOT. topography == 'flat' ) THEN3682 IF ( .NOT. topography == 'flat' ) THEN 3656 3683 i = ip + MERGE( -1_iwp , 1_iwp, BTEST( nb, 2 ) ) 3657 3684 j = jp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 1 ) ) … … 3664 3691 ENDIF 3665 3692 ENDIF 3666 IF ( subbox_at_wall ) THEN3693 IF ( subbox_at_wall ) THEN 3667 3694 e_int(start_index(nb):end_index(nb)) = e(kp,jp,ip) 3668 3695 diss_int(start_index(nb):end_index(nb)) = diss(kp,jp,ip) … … 3960 3987 (ABS( u_int(n) + rvar1_temp(n) ) > (dx/dt_particle(n)) .OR. & 3961 3988 ABS( v_int(n) + rvar2_temp(n) ) > (dy/dt_particle(n)) .OR. & 3962 ABS( w_int(n) + rvar3_temp(n) ) > (dz_temp/dt_particle(n)))) THEN3989 ABS( w_int(n) + rvar3_temp(n) ) > (dz_temp/dt_particle(n)))) THEN 3963 3990 3964 3991 dt_particle(n) = 0.9_wp * MIN( & … … 4027 4054 !-- number_of_particles. @todo find a more generic way to write this loop or 4028 4055 !-- delete trilinear interpolation 4029 IF ( TRIM(particle_interpolation) == 'trilinear') THEN4056 IF ( interpolation_trilinear ) THEN 4030 4057 subbox_start = 0 4031 4058 subbox_end = 7 … … 4039 4066 !-- from 1 to 1. 4040 4067 DO nb = subbox_start, subbox_end 4041 IF ( TRIM(particle_interpolation) == 'trilinear') THEN4068 IF ( interpolation_trilinear ) THEN 4042 4069 particle_start = start_index(nb) 4043 4070 particle_end = end_index(nb) … … 4141 4168 !-- Decide whether the particle loop runs over the subboxes or only over 1, 4142 4169 !-- number_of_particles. This depends on the selected interpolation method. 4143 IF ( TRIM(particle_interpolation) == 'trilinear') THEN4170 IF ( interpolation_trilinear ) THEN 4144 4171 subbox_start = 0 4145 4172 subbox_end = 7 … … 4152 4179 !-- from 1 to 1. 4153 4180 DO nb = subbox_start, subbox_end 4154 IF ( TRIM(particle_interpolation) == 'trilinear') THEN4181 IF ( interpolation_trilinear ) THEN 4155 4182 particle_start = start_index(nb) 4156 4183 particle_end = end_index(nb) … … 4245 4272 !-- Decide whether the particle loop runs over the subboxes or only over 1, 4246 4273 !-- number_of_particles. This depends on the selected interpolation method. 4247 IF ( TRIM(particle_interpolation) == 'trilinear') THEN4274 IF ( interpolation_trilinear ) THEN 4248 4275 subbox_start = 0 4249 4276 subbox_end = 7 … … 4253 4280 ENDIF 4254 4281 DO nb = subbox_start, subbox_end 4255 IF ( TRIM(particle_interpolation) == 'trilinear') THEN4282 IF ( interpolation_trilinear ) THEN 4256 4283 particle_start = start_index(nb) 4257 4284 particle_end = end_index(nb) … … 4505 4532 particles(n)%z = reset_top * ( 1.0 + ( ran_val / 10.0_wp) ) 4506 4533 particles(n)%speed_z = 0.0_wp 4507 IF ( curvature_solution_effects ) THEN4534 IF ( curvature_solution_effects ) THEN 4508 4535 particles(n)%radius = particles(n)%aux1 4509 4536 ELSE … … 4566 4593 ENDIF 4567 4594 4568 IF ( prt_z > pos_z_old ) THEN4595 IF ( prt_z > pos_z_old ) THEN 4569 4596 zwall = zw(k) 4570 4597 ELSE … … 4648 4675 reach_y(t_index) = .FALSE. 4649 4676 reach_z(t_index) = .TRUE. 4650 IF( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp) THEN4677 IF( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp) THEN 4651 4678 t_index = t_index + 1 4652 4679 cross_wall_z = .TRUE. … … 5001 5028 ENDDO 5002 5029 5003 IF( .NOT. curvature_solution_effects ) then5030 IF( .NOT. curvature_solution_effects ) THEN 5004 5031 ! 5005 5032 !-- Use analytic model for diffusional growth including gas-kinetic … … 6136 6163 6137 6164 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ira !< 6138 6165 6139 6166 LOGICAL, SAVE :: first = .TRUE. !< 6140 6167 … … 6149 6176 REAL(wp), DIMENSION(1:11), SAVE :: rat !< 6150 6177 REAL(wp), DIMENSION(1:7), SAVE :: r0 !< 6151 6178 6152 6179 REAL(wp), DIMENSION(1:7,1:11), SAVE :: ecoll_100 !< 6153 6180 REAL(wp), DIMENSION(1:7,1:11), SAVE :: ecoll_400 !< … … 7010 7037 DO kp = nzb+1, nzt 7011 7038 number_of_particles = prt_count(kp,jp,ip) 7012 IF ( number_of_particles <= 0 ) CYCLE7039 IF ( number_of_particles <= 0 ) CYCLE 7013 7040 particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles) 7014 7041 DO n = 1, number_of_particles … … 7589 7616 !-- In case of stretching the actual k index must be found 7590 7617 IF ( dz_stretch_level .NE. -9999999.9_wp .OR. & 7591 dz_stretch_level_start(1) .NE. -9999999.9_wp ) THEN7618 dz_stretch_level_start(1) .NE. -9999999.9_wp ) THEN 7592 7619 kp = MINLOC( ABS( particle_array(n)%z - zu ), DIM = 1 ) - 1 7593 7620 ELSE … … 7992 8019 ! ----------- 7993 8020 !> Routine for the whole processor 7994 !> Sort all particles into the 8 respective subgrid boxes and free space of 7995 !> particles which has been marked for deletion 8021 !> Sort all particles into the 8 respective subgrid boxes (in case of trilinear 8022 !> interpolation method) and free space of particles which has been marked for 8023 !> deletion. 7996 8024 !------------------------------------------------------------------------------! 7997 8025 SUBROUTINE lpm_sort_and_delete … … 8013 8041 8014 8042 CALL cpu_log( log_point_s(51), 'lpm_sort_and_delete', 'start' ) 8015 IF ( TRIM(particle_interpolation) == 'trilinear' )THEN8043 IF ( interpolation_trilinear ) THEN 8016 8044 DO ip = nxl, nxr 8017 8045 DO jp = nys, nyn … … 8047 8075 ENDIF 8048 8076 ENDDO 8049 8077 ! 8078 !-- Delete and resort particles by overwritting and set 8079 !-- the number_of_particles to the actual value. 8050 8080 nn = 0 8051 8081 DO is = 0,7 … … 8066 8096 8067 8097 !-- In case of the simple interpolation method the particles must not 8068 !-- be sorted in subboxes but particles marked for deletion must be 8069 !-- deleted and number of particles must be recalculated 8098 !-- be sorted in subboxes. Particles marked for deletion however, must be 8099 !-- deleted and number of particles must be recalculated as it is also 8100 !-- done for the trilinear particle advection interpolation method. 8070 8101 ELSE 8071 8102 … … 8078 8109 particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles) 8079 8110 ! 8080 !-- repack particles array, i.e. delete particles and recalculate8111 !-- Repack particles array, i.e. delete particles and recalculate 8081 8112 !-- number of particles 8082 8113 CALL lpm_pack
Note: See TracChangeset
for help on using the changeset viewer.
