Changeset 1929 for palm/trunk
- Timestamp:
- Jun 9, 2016 4:25:25 PM (8 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 11 edited
-
check_parameters.f90 (modified) (2 diffs)
-
lpm.f90 (modified) (5 diffs)
-
lpm_advec.f90 (modified) (38 diffs)
-
lpm_boundary_conds.f90 (modified) (3 diffs)
-
lpm_exchange_horiz.f90 (modified) (35 diffs)
-
lpm_init.f90 (modified) (17 diffs)
-
lpm_init_sgs_tke.f90 (modified) (8 diffs)
-
mod_particle_attributes.f90 (modified) (2 diffs)
-
pres.f90 (modified) (7 diffs)
-
surface_layer_fluxes_mod.f90 (modified) (2 diffs)
-
wind_turbine_model_mod.f90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/check_parameters.f90
r1919 r1929 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Bugfix in check for use_upstream_for_tke 22 22 ! 23 23 ! Former revisions: … … 894 894 IF ( use_sgs_for_particles .AND. .NOT. cloud_droplets .AND. & 895 895 .NOT. use_upstream_for_tke .AND. & 896 ( scalar_advec /= 'ws-scheme' . OR. scalar_advec /= 'ws-scheme-mono' )&896 ( scalar_advec /= 'ws-scheme' .AND. scalar_advec /= 'ws-scheme-mono' )& 897 897 ) THEN 898 898 use_upstream_for_tke = .TRUE. -
palm/trunk/SOURCE/lpm.f90
r1823 r1929 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! Call wall boundary conditions only if particles are in the vertical range of 22 ! topography. 22 23 ! 23 24 ! Former revisions: … … 111 112 112 113 USE indices, & 113 ONLY: nxl, nxr, nys, nyn, nzb, nz t114 ONLY: nxl, nxr, nys, nyn, nzb, nzb_max, nzt, nzb_w_inner 114 115 115 116 USE kinds … … 300 301 CALL lpm_advec(i,j,k) 301 302 ! 302 !-- Particle reflection from walls 303 IF ( topography /= 'flat' ) THEN 303 !-- Particle reflection from walls. Only applied if the particles 304 !-- are in the vertical range of the topography. (Here, some 305 !-- optimization is still possible.) 306 IF ( topography /= 'flat' .AND. k < nzb_max + 2 ) THEN 304 307 CALL lpm_boundary_conds( 'walls' ) 305 308 ENDIF … … 334 337 steps = steps + 1 335 338 dt_3d_reached_l = ALL(grid_particles(:,:,:)%time_loop_done) 336 337 339 ! 338 340 !-- Find out, if all particles on every PE have completed the LES timestep … … 359 361 !-- Horizontal boundary conditions including exchange between subdmains 360 362 CALL lpm_exchange_horiz 361 362 363 ! 363 364 !-- Pack particles (eliminate those marked for deletion), -
palm/trunk/SOURCE/lpm_advec.f90
r1889 r1929 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! Put stochastic equation in an extra subroutine. 22 ! Set flag for stochastic equation to communicate whether a particle is near 23 ! topography. This case, memory and drift term are disabled in the Weil equation. 24 ! 25 ! Enable vertical logarithmic interpolation also above topography. This case, 26 ! set a lower limit for the friction velocity, as it can become very small 27 ! in narrow street canyons, leading to too large particle speeds. 22 28 ! 23 29 ! Former revisions: … … 110 116 ONLY: block_offset, c_0, dt_min_part, grid_particles, & 111 117 iran_part, log_z_z0, number_of_particles, number_of_sublayers, & 112 particles, particle_groups, offset_ocean_nzt, sgs_wfu_part, & 113 sgs_wfv_part, sgs_wfw_part, use_sgs_for_particles, & 114 vertical_particle_advection, z0_av_global 118 particles, particle_groups, offset_ocean_nzt, sgs_wf_part, & 119 use_sgs_for_particles, vertical_particle_advection, z0_av_global 115 120 116 121 USE statistics, & … … 119 124 IMPLICIT NONE 120 125 121 INTEGER(iwp) :: agp !< 122 INTEGER(iwp) :: gp_outside_of_building(1:8) !< 123 INTEGER(iwp) :: i !< 124 INTEGER(iwp) :: ip !< 125 INTEGER(iwp) :: j !< 126 INTEGER(iwp) :: jp !< 127 INTEGER(iwp) :: k !< 128 INTEGER(iwp) :: kp !< 129 INTEGER(iwp) :: kw !< 130 INTEGER(iwp) :: n !< 131 INTEGER(iwp) :: nb !< 132 INTEGER(iwp) :: num_gp !< 133 134 INTEGER(iwp), DIMENSION(0:7) :: start_index !< 135 INTEGER(iwp), DIMENSION(0:7) :: end_index !< 136 137 REAL(wp) :: aa !< 138 REAL(wp) :: bb !< 139 REAL(wp) :: cc !< 140 REAL(wp) :: d_sum !< 141 REAL(wp) :: d_z_p_z0 !< 142 REAL(wp) :: dd !< 143 REAL(wp) :: de_dx_int_l !< 144 REAL(wp) :: de_dx_int_u !< 145 REAL(wp) :: de_dy_int_l !< 146 REAL(wp) :: de_dy_int_u !< 147 REAL(wp) :: de_dt !< 148 REAL(wp) :: de_dt_min !< 149 REAL(wp) :: de_dz_int_l !< 150 REAL(wp) :: de_dz_int_u !< 126 INTEGER(iwp) :: agp !< loop variable 127 INTEGER(iwp) :: gp_outside_of_building(1:8) !< number of grid points used for particle interpolation in case of topography 128 INTEGER(iwp) :: i !< index variable along x 129 INTEGER(iwp) :: ip !< index variable along x 130 INTEGER(iwp) :: ilog !< index variable along x 131 INTEGER(iwp) :: j !< index variable along y 132 INTEGER(iwp) :: jp !< index variable along y 133 INTEGER(iwp) :: jlog !< index variable along y 134 INTEGER(iwp) :: k !< index variable along z 135 INTEGER(iwp) :: kp !< index variable along z 136 INTEGER(iwp) :: kw !< index variable along z 137 INTEGER(iwp) :: n !< loop variable over all particles in a grid box 138 INTEGER(iwp) :: nb !< block number particles are sorted in 139 INTEGER(iwp) :: num_gp !< number of adjacent grid points inside topography 140 141 INTEGER(iwp), DIMENSION(0:7) :: start_index !< start particle index for current block 142 INTEGER(iwp), DIMENSION(0:7) :: end_index !< start particle index for current block 143 144 REAL(wp) :: aa !< dummy argument for horizontal particle interpolation 145 REAL(wp) :: bb !< dummy argument for horizontal particle interpolation 146 REAL(wp) :: cc !< dummy argument for horizontal particle interpolation 147 REAL(wp) :: d_sum !< dummy argument for horizontal particle interpolation in case of topography 148 REAL(wp) :: d_z_p_z0 !< inverse of interpolation length for logarithmic interpolation 149 REAL(wp) :: dd !< dummy argument for horizontal particle interpolation 150 REAL(wp) :: de_dx_int_l !< x/y-interpolated TKE gradient (x) at particle position at lower vertical level 151 REAL(wp) :: de_dx_int_u !< x/y-interpolated TKE gradient (x) at particle position at upper vertical level 152 REAL(wp) :: de_dy_int_l !< x/y-interpolated TKE gradient (y) at particle position at lower vertical level 153 REAL(wp) :: de_dy_int_u !< x/y-interpolated TKE gradient (y) at particle position at upper vertical level 154 REAL(wp) :: de_dt !< temporal derivative of TKE experienced by the particle 155 REAL(wp) :: de_dt_min !< lower level for temporal TKE derivative 156 REAL(wp) :: de_dz_int_l !< x/y-interpolated TKE gradient (z) at particle position at lower vertical level 157 REAL(wp) :: de_dz_int_u !< x/y-interpolated TKE gradient (z) at particle position at upper vertical level 151 158 REAL(wp) :: diameter !< diamter of droplet 152 REAL(wp) :: diss_int_l !< 153 REAL(wp) :: diss_int_u !< 154 REAL(wp) :: dt_gap !< 155 REAL(wp) :: dt_particle_m !< 156 REAL(wp) :: e_int_l !< 157 REAL(wp) :: e_int_u !< 158 REAL(wp) :: e_mean_int !< 159 REAL(wp) :: diss_int_l !< x/y-interpolated dissipation at particle position at lower vertical level 160 REAL(wp) :: diss_int_u !< x/y-interpolated dissipation at particle position at upper vertical level 161 REAL(wp) :: dt_gap !< remaining time until particle time integration reaches LES time 162 REAL(wp) :: dt_particle_m !< previous particle time step 163 REAL(wp) :: e_int_l !< x/y-interpolated TKE at particle position at lower vertical level 164 REAL(wp) :: e_int_u !< x/y-interpolated TKE at particle position at upper vertical level 165 REAL(wp) :: e_mean_int !< horizontal mean TKE at particle height 159 166 REAL(wp) :: exp_arg !< 160 167 REAL(wp) :: exp_term !< 161 REAL(wp) :: gg !< 162 REAL(wp) :: height_p !< 168 REAL(wp) :: gg !< dummy argument for horizontal particle interpolation 169 REAL(wp) :: height_p !< dummy argument for logarithmic interpolation 163 170 REAL(wp) :: lagr_timescale !< Lagrangian timescale 164 REAL(wp) :: location(1:30,1:3) !< 171 REAL(wp) :: location(1:30,1:3) !< wall locations 172 REAL(wp) :: log_z_z0_int !< logarithmus used for surface_layer interpolation 165 173 REAL(wp) :: random_gauss !< 166 174 REAL(wp) :: RL !< Lagrangian autocorrelation coefficient … … 169 177 REAL(wp) :: rg3 !< Gaussian distributed random number 170 178 REAL(wp) :: sigma !< velocity standard deviation 171 REAL(wp) :: u_int_l !< 172 REAL(wp) :: u_int_u !< 173 REAL(wp) :: us_int !< 174 REAL(wp) :: v_int_l !< 175 REAL(wp) :: v_int_u !< 179 REAL(wp) :: u_int_l !< x/y-interpolated u-component at particle position at lower vertical level 180 REAL(wp) :: u_int_u !< x/y-interpolated u-component at particle position at upper vertical level 181 REAL(wp) :: us_int !< friction velocity at particle grid box 182 REAL(wp) :: v_int_l !< x/y-interpolated v-component at particle position at lower vertical level 183 REAL(wp) :: v_int_u !< x/y-interpolated v-component at particle position at upper vertical level 176 184 REAL(wp) :: vv_int !< 177 REAL(wp) :: w_int_l !< 178 REAL(wp) :: w_int_u !< 185 REAL(wp) :: w_int_l !< x/y-interpolated w-component at particle position at lower vertical level 186 REAL(wp) :: w_int_u !< x/y-interpolated w-component at particle position at upper vertical level 179 187 REAL(wp) :: w_s !< terminal velocity of droplets 180 REAL(wp) :: x !< 181 REAL(wp) :: y !< 182 REAL(wp) :: z_p !< 188 REAL(wp) :: x !< dummy argument for horizontal particle interpolation 189 REAL(wp) :: y !< dummy argument for horizontal particle interpolation 190 REAL(wp) :: z_p !< surface layer height (0.5 dz) 183 191 184 192 REAL(wp), PARAMETER :: a_rog = 9.65_wp !< parameter for fall velocity … … 189 197 REAL(wp), PARAMETER :: d0_rog = 0.745_wp !< separation diameter 190 198 191 REAL(wp), DIMENSION(1:30) :: d_gp_pl !< 192 REAL(wp), DIMENSION(1:30) :: de_dxi !< 193 REAL(wp), DIMENSION(1:30) :: de_dyi !< 194 REAL(wp), DIMENSION(1:30) :: de_dzi !< 195 REAL(wp), DIMENSION(1:30) :: dissi !< 196 REAL(wp), DIMENSION(1:30) :: ei !< 197 199 REAL(wp), DIMENSION(1:30) :: d_gp_pl !< dummy argument for particle interpolation scheme in case of topography 200 REAL(wp), DIMENSION(1:30) :: de_dxi !< horizontal TKE gradient along x at adjacent wall 201 REAL(wp), DIMENSION(1:30) :: de_dyi !< horizontal TKE gradient along y at adjacent wall 202 REAL(wp), DIMENSION(1:30) :: de_dzi !< horizontal TKE gradient along z at adjacent wall 203 REAL(wp), DIMENSION(1:30) :: dissi !< dissipation at adjacent wall 204 REAL(wp), DIMENSION(1:30) :: ei !< TKE at adjacent wall 205 206 REAL(wp), DIMENSION(number_of_particles) :: term_1_2 !< flag to communicate whether a particle is near topography or not 198 207 REAL(wp), DIMENSION(number_of_particles) :: dens_ratio !< 199 REAL(wp), DIMENSION(number_of_particles) :: de_dx_int !< 200 REAL(wp), DIMENSION(number_of_particles) :: de_dy_int !< 201 REAL(wp), DIMENSION(number_of_particles) :: de_dz_int !< 202 REAL(wp), DIMENSION(number_of_particles) :: diss_int !< 203 REAL(wp), DIMENSION(number_of_particles) :: dt_particle !< 204 REAL(wp), DIMENSION(number_of_particles) :: e_int !< 205 REAL(wp), DIMENSION(number_of_particles) :: fs_int !< 206 REAL(wp), DIMENSION(number_of_particles) :: log_z_z0_int !< 207 REAL(wp), DIMENSION(number_of_particles) :: u_int !< 208 REAL(wp), DIMENSION(number_of_particles) :: v_int !< 209 REAL(wp), DIMENSION(number_of_particles) :: w_int !< 210 REAL(wp), DIMENSION(number_of_particles) :: xv !< 211 REAL(wp), DIMENSION(number_of_particles) :: yv !< 212 REAL(wp), DIMENSION(number_of_particles) :: zv !< 213 214 REAL(wp), DIMENSION(number_of_particles, 3) :: rg !< 208 REAL(wp), DIMENSION(number_of_particles) :: de_dx_int !< horizontal TKE gradient along x at particle position 209 REAL(wp), DIMENSION(number_of_particles) :: de_dy_int !< horizontal TKE gradient along y at particle position 210 REAL(wp), DIMENSION(number_of_particles) :: de_dz_int !< horizontal TKE gradient along z at particle position 211 REAL(wp), DIMENSION(number_of_particles) :: diss_int !< dissipation at particle position 212 REAL(wp), DIMENSION(number_of_particles) :: dt_particle !< particle time step 213 REAL(wp), DIMENSION(number_of_particles) :: e_int !< TKE at particle position 214 REAL(wp), DIMENSION(number_of_particles) :: fs_int !< weighting factor for subgrid-scale particle speed 215 REAL(wp), DIMENSION(number_of_particles) :: u_int !< u-component of particle speed 216 REAL(wp), DIMENSION(number_of_particles) :: v_int !< v-component of particle speed 217 REAL(wp), DIMENSION(number_of_particles) :: w_int !< w-component of particle speed 218 REAL(wp), DIMENSION(number_of_particles) :: xv !< x-position 219 REAL(wp), DIMENSION(number_of_particles) :: yv !< y-position 220 REAL(wp), DIMENSION(number_of_particles) :: zv !< z-position 221 222 REAL(wp), DIMENSION(number_of_particles, 3) :: rg !< vector of Gaussian distributed random numbers 215 223 216 224 CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' ) … … 236 244 j = jp + block_offset(nb)%j_off 237 245 k = kp + block_offset(nb)%k_off 246 238 247 239 248 ! … … 249 258 !-- First, check if particle is located below first vertical grid level 250 259 !-- (Prandtl-layer height) 251 IF ( constant_flux_layer .AND. particles(n)%z < z_p ) THEN 260 ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx 261 jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy 262 263 IF ( constant_flux_layer .AND. zv(n) - zw(nzb_s_inner(jlog,ilog)) < z_p ) THEN 252 264 ! 253 265 !-- Resolved-scale horizontal particle velocity is zero below z0. 254 IF ( particles(n)%z< z0_av_global ) THEN266 IF ( zv(n) - zw(nzb_s_inner(jlog,ilog)) < z0_av_global ) THEN 255 267 u_int(n) = 0.0_wp 256 268 ELSE 257 269 ! 258 270 !-- Determine the sublayer. Further used as index. 259 height_p = ( particles(n)%z- z0_av_global ) &271 height_p = ( zv(n) - zw(nzb_s_inner(jlog,ilog)) - z0_av_global ) & 260 272 * REAL( number_of_sublayers, KIND=wp ) & 261 273 * d_z_p_z0 … … 263 275 !-- Calculate LOG(z/z0) for exact particle height. Therefore, 264 276 !-- interpolate linearly between precalculated logarithm. 265 log_z_z0_int (n) = log_z_z0(INT(height_p))&277 log_z_z0_int = log_z_z0(INT(height_p)) & 266 278 + ( height_p - INT(height_p) ) & 267 279 * ( log_z_z0(INT(height_p)+1) & 268 280 - log_z_z0(INT(height_p)) & 269 281 ) 282 ! 283 !-- Limit friction velocity. In narrow canyons or holes the 284 !-- friction velocity can become very small, resulting in a too 285 !-- large particle speed. 286 us_int = MAX( 0.5_wp * ( us(jlog,ilog) + us(jlog,ilog-1) ), & 287 0.01_wp ) 270 288 ! 271 289 !-- Neutral solution is applied for all situations, e.g. also for … … 275 293 !-- as sensitivity studies revealed no significant effect of 276 294 !-- using the neutral solution also for un/stable situations. 277 !-- Calculated left and bottom index on u grid. 278 us_int = 0.5_wp * ( us(j,i) + us(j,i-1) ) 279 280 u_int(n) = -usws(j,i) / ( us_int * kappa + 1E-10_wp ) & 281 * log_z_z0_int(n) - u_gtrans 282 295 u_int(n) = -usws(jlog,ilog) / ( us_int * kappa + 1E-10_wp ) & 296 * log_z_z0_int - u_gtrans 297 283 298 ENDIF 284 299 ! … … 308 323 ( u_int_u - u_int_l ) 309 324 ENDIF 325 310 326 ENDIF 311 327 … … 319 335 DO n = start_index(nb), end_index(nb) 320 336 321 IF ( constant_flux_layer .AND. particles(n)%z < z_p ) THEN 322 323 IF ( particles(n)%z < z0_av_global ) THEN 337 ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx 338 jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy 339 IF ( constant_flux_layer .AND. zv(n) - zw(nzb_s_inner(jlog,ilog)) < z_p ) THEN 340 341 IF ( zv(n) - zw(nzb_s_inner(jlog,ilog)) < z0_av_global ) THEN 324 342 ! 325 343 !-- Resolved-scale horizontal particle velocity is zero below z0. 326 344 v_int(n) = 0.0_wp 327 345 ELSE 346 ! 347 !-- Determine the sublayer. Further used as index. Please note, 348 !-- logarithmus can not be reused from above, as in in case of 349 !-- topography particle on u-grid can be above surface-layer height, 350 !-- whereas it can be below on v-grid. 351 height_p = ( zv(n) - zw(nzb_s_inner(jlog,ilog)) - z0_av_global ) & 352 * REAL( number_of_sublayers, KIND=wp ) & 353 * d_z_p_z0 354 ! 355 !-- Calculate LOG(z/z0) for exact particle height. Therefore, 356 !-- interpolate linearly between precalculated logarithm. 357 log_z_z0_int = log_z_z0(INT(height_p)) & 358 + ( height_p - INT(height_p) ) & 359 * ( log_z_z0(INT(height_p)+1) & 360 - log_z_z0(INT(height_p)) & 361 ) 362 ! 363 !-- Limit friction velocity. In narrow canyons or holes the 364 !-- friction velocity can become very small, resulting in a too 365 !-- large particle speed. 366 us_int = MAX( 0.5_wp * ( us(jlog,ilog) + us(jlog-1,ilog) ), & 367 0.01_wp ) 328 368 ! 329 369 !-- Neutral solution is applied for all situations, e.g. also for … … 333 373 !-- as sensitivity studies revealed no significant effect of 334 374 !-- using the neutral solution also for un/stable situations. 335 !-- Calculated left and bottom index on v grid. 336 us_int = 0.5_wp * ( us(j,i) + us(j-1,i) ) 337 338 v_int(n) = -vsws(j,i) / ( us_int * kappa + 1E-10_wp ) & 339 * log_z_z0_int(n) - v_gtrans 340 ENDIF 375 v_int(n) = -vsws(jlog,ilog) / ( us_int * kappa + 1E-10_wp ) & 376 * log_z_z0_int - v_gtrans 377 378 ENDIF 379 341 380 ELSE 342 381 x = xv(n) - i * dx … … 361 400 ( v_int_u - v_int_l ) 362 401 ENDIF 402 363 403 ENDIF 364 404 … … 367 407 i = ip + block_offset(nb)%i_off 368 408 j = jp + block_offset(nb)%j_off 369 k = kp -1409 k = kp - 1 370 410 ! 371 411 !-- Same procedure for interpolation of the w velocity-component … … 535 575 ENDIF 536 576 577 ! 578 !-- Set flag for stochastic equation. 579 term_1_2(n) = 1.0_wp 580 537 581 ENDDO 538 582 ENDDO … … 541 585 542 586 DO n = 1, number_of_particles 543 544 587 i = particles(n)%x * ddx 545 588 j = particles(n)%y * ddy … … 575 618 de_dzi(num_gp) = de_dz(k,j,i) 576 619 ENDIF 577 578 IF ( k > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) & 579 THEN 620 IF ( k > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) THEN 580 621 num_gp = num_gp + 1 581 622 gp_outside_of_building(2) = 1 … … 603 644 ENDIF 604 645 605 IF ( k+1 > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) & 606 THEN 646 IF ( k+1 > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) THEN 607 647 num_gp = num_gp + 1 608 648 gp_outside_of_building(4) = 1 … … 617 657 ENDIF 618 658 619 IF ( k > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) & 620 THEN 659 IF ( k > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) THEN 621 660 num_gp = num_gp + 1 622 661 gp_outside_of_building(5) = 1 … … 631 670 ENDIF 632 671 633 IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) & 634 THEN 672 IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) THEN 635 673 num_gp = num_gp + 1 636 674 gp_outside_of_building(6) = 1 … … 645 683 ENDIF 646 684 647 IF ( k+1 > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) & 648 THEN 685 IF ( k+1 > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) THEN 649 686 num_gp = num_gp + 1 650 687 gp_outside_of_building(7) = 1 … … 659 696 ENDIF 660 697 661 IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0)& 662 THEN 698 IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0) THEN 663 699 num_gp = num_gp + 1 664 700 gp_outside_of_building(8) = 1 … … 672 708 de_dzi(num_gp) = de_dz(k+1,j+1,i+1) 673 709 ENDIF 674 675 710 ! 676 711 !-- If all gridpoints are situated outside of a building, then the … … 685 720 dd = ( dx - x )**2 + ( dy - y )**2 686 721 gg = aa + bb + cc + dd 687 722 688 723 e_int_l = ( ( gg - aa ) * e(k,j,i) + ( gg - bb ) * e(k,j,i+1) & 689 724 + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) & 690 725 ) / ( 3.0_wp * gg ) 691 726 692 727 IF ( k == nzt ) THEN 693 728 e_int(n) = e_int_l … … 699 734 ) / ( 3.0_wp * gg ) 700 735 e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * & 701 ( e_int_u - e_int_l )702 ENDIF 703 ! 736 ( e_int_u - e_int_l ) 737 ENDIF 738 ! 704 739 !-- Needed to avoid NaN particle velocities (this might not be 705 740 !-- required any more) … … 757 792 ( gg - dd ) * de_dz(k,j+1,i+1) & 758 793 ) / ( 3.0_wp * gg ) 759 794 760 795 IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN 761 796 de_dz_int(n) = de_dz_int_l … … 778 813 ( gg - dd ) * diss(k,j+1,i+1) & 779 814 ) / ( 3.0_wp * gg ) 780 815 781 816 IF ( k == nzt ) THEN 782 817 diss_int(n) = diss_int_l … … 790 825 ( diss_int_u - diss_int_l ) 791 826 ENDIF 792 827 ! 828 !-- Set flag for stochastic equation. 829 term_1_2(n) = 1.0_wp 830 793 831 ELSE 794 832 795 833 ! 796 834 !-- If wall between gridpoint 1 and gridpoint 5, then … … 810 848 811 849 IF ( gp_outside_of_building(5) == 1 .AND. & 812 gp_outside_of_building(1) == 0 ) THEN850 gp_outside_of_building(1) == 0 ) THEN 813 851 num_gp = num_gp + 1 814 852 location(num_gp,1) = i * dx + 0.5_wp * dx … … 893 931 de_dxi(num_gp) = de_dx(k,j,i) 894 932 de_dyi(num_gp) = 0.0_wp 895 de_dzi(num_gp) = de_dz(k,j,i) 933 de_dzi(num_gp) = de_dz(k,j,i) 896 934 ENDIF 897 935 … … 1076 1114 ENDIF 1077 1115 1078 ! 1116 ! 1079 1117 !-- If wall between gridpoint 6 and gridpoint 8, then 1080 1118 !-- Neumann boundary condition has to be applied … … 1092 1130 de_dzi(num_gp) = 0.0_wp 1093 1131 ENDIF 1094 1132 1095 1133 ! 1096 1134 !-- Carry out the interpolation 1097 1135 IF ( num_gp == 1 ) THEN 1098 ! 1136 ! 1099 1137 !-- If only one of the gridpoints is situated outside of the 1100 1138 !-- building, it follows that the values at the particle 1101 1139 !-- location are the same as the gridpoint values 1102 e_int(n) = ei(num_gp) 1103 diss_int(n) = dissi(num_gp) 1140 e_int(n) = ei(num_gp) 1141 diss_int(n) = dissi(num_gp) 1104 1142 de_dx_int(n) = de_dxi(num_gp) 1105 1143 de_dy_int(n) = de_dyi(num_gp) 1106 1144 de_dz_int(n) = de_dzi(num_gp) 1145 ! 1146 !-- Set flag for stochastic equation which disables calculation 1147 !-- of drift and memory term near topography. 1148 term_1_2(n) = 0.0_wp 1107 1149 ELSE IF ( num_gp > 1 ) THEN 1108 1150 1109 1151 d_sum = 0.0_wp 1110 ! 1152 ! 1111 1153 !-- Evaluation of the distances between the gridpoints 1112 1154 !-- contributing to the interpolated values, and the particle … … 1118 1160 d_sum = d_sum + d_gp_pl(agp) 1119 1161 ENDDO 1120 1162 1121 1163 ! 1122 1164 !-- Finally the interpolation can be carried out 1123 1165 e_int(n) = 0.0_wp 1124 1166 diss_int(n) = 0.0_wp 1125 de_dx_int(n) = 0.0_wp 1126 de_dy_int(n) = 0.0_wp 1127 de_dz_int(n) = 0.0_wp 1167 de_dx_int(n) = 0.0_wp 1168 de_dy_int(n) = 0.0_wp 1169 de_dz_int(n) = 0.0_wp 1128 1170 DO agp = 1, num_gp 1129 1171 e_int(n) = e_int(n) + ( d_sum - d_gp_pl(agp) ) * & … … 1138 1180 de_dzi(agp) / ( (num_gp-1) * d_sum ) 1139 1181 ENDDO 1140 1141 ENDIF 1142 1182 1183 ENDIF 1184 e_int(n) = MAX( 1E-20_wp, e_int(n) ) 1185 diss_int(n) = MAX( 1E-20_wp, diss_int(n) ) 1186 de_dx_int(n) = MAX( 1E-20_wp, de_dx_int(n) ) 1187 de_dy_int(n) = MAX( 1E-20_wp, de_dy_int(n) ) 1188 de_dz_int(n) = MAX( 1E-20_wp, de_dz_int(n) ) 1189 ! 1190 !-- Set flag for stochastic equation which disables calculation 1191 !-- of drift and memory term near topography. 1192 term_1_2(n) = 0.0_wp 1143 1193 ENDIF 1144 1194 ENDDO … … 1209 1259 ! 1210 1260 !-- Calculate the Lagrangian timescale according to Weil et al. (2004). 1211 lagr_timescale = ( 4.0_wp * e_int(n) ) / &1212 ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) )1261 lagr_timescale = ( 4.0_wp * e_int(n) + 1E-20_wp ) / & 1262 ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp ) 1213 1263 1214 1264 ! … … 1233 1283 !-- [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities 1234 1284 !-- from becoming unrealistically large. 1235 particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf u_part * e_int(n)) * &1285 particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) * & 1236 1286 ( rg(n,1) - 1.0_wp ) 1237 particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf v_part * e_int(n)) * &1287 particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) * & 1238 1288 ( rg(n,2) - 1.0_wp ) 1239 particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf w_part * e_int(n)) * &1289 particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) * & 1240 1290 ( rg(n,3) - 1.0_wp ) 1241 1291 … … 1267 1317 ENDIF 1268 1318 1269 particles(n)%rvar1 = particles(n)%rvar1 - fs_int(n) * c_0 * & 1270 diss_int(n) * particles(n)%rvar1 * dt_particle(n) / & 1271 ( 4.0_wp * sgs_wfu_part * e_int(n) ) + & 1272 ( 2.0_wp * sgs_wfu_part * de_dt * & 1273 particles(n)%rvar1 / & 1274 ( 2.0_wp * sgs_wfu_part * e_int(n) ) + & 1275 de_dx_int(n) & 1276 ) * dt_particle(n) / 2.0_wp + & 1277 SQRT( fs_int(n) * c_0 * diss_int(n) ) * & 1278 ( rg(n,1) - 1.0_wp ) * & 1279 SQRT( dt_particle(n) ) 1280 1281 particles(n)%rvar2 = particles(n)%rvar2 - fs_int(n) * c_0 * & 1282 diss_int(n) * particles(n)%rvar2 * dt_particle(n) / & 1283 ( 4.0_wp * sgs_wfv_part * e_int(n) ) + & 1284 ( 2.0_wp * sgs_wfv_part * de_dt * & 1285 particles(n)%rvar2 / & 1286 ( 2.0_wp * sgs_wfv_part * e_int(n) ) + & 1287 de_dy_int(n) & 1288 ) * dt_particle(n) / 2.0_wp + & 1289 SQRT( fs_int(n) * c_0 * diss_int(n) ) * & 1290 ( rg(n,2) - 1.0_wp ) * & 1291 SQRT( dt_particle(n) ) 1292 1293 particles(n)%rvar3 = particles(n)%rvar3 - fs_int(n) * c_0 * & 1294 diss_int(n) * particles(n)%rvar3 * dt_particle(n) / & 1295 ( 4.0_wp * sgs_wfw_part * e_int(n) ) + & 1296 ( 2.0_wp * sgs_wfw_part * de_dt * & 1297 particles(n)%rvar3 / & 1298 ( 2.0_wp * sgs_wfw_part * e_int(n) ) + & 1299 de_dz_int(n) & 1300 ) * dt_particle(n) / 2.0_wp + & 1301 SQRT( fs_int(n) * c_0 * diss_int(n) ) * & 1302 ( rg(n,3) - 1.0_wp ) * & 1303 SQRT( dt_particle(n) ) 1319 CALL weil_stochastic_eq(particles(n)%rvar1, fs_int(n), e_int(n), & 1320 de_dx_int(n), de_dt, diss_int(n), & 1321 dt_particle(n), rg(n,1), term_1_2(n) ) 1322 1323 CALL weil_stochastic_eq(particles(n)%rvar2, fs_int(n), e_int(n), & 1324 de_dy_int(n), de_dt, diss_int(n), & 1325 dt_particle(n), rg(n,2), term_1_2(n) ) 1326 1327 CALL weil_stochastic_eq(particles(n)%rvar3, fs_int(n), e_int(n), & 1328 de_dz_int(n), de_dt, diss_int(n), & 1329 dt_particle(n), rg(n,3), term_1_2(n) ) 1304 1330 1305 1331 ENDIF 1332 1306 1333 u_int(n) = u_int(n) + particles(n)%rvar1 1307 1334 v_int(n) = v_int(n) + particles(n)%rvar2 1308 1335 w_int(n) = w_int(n) + particles(n)%rvar3 1309 1310 1336 ! 1311 1337 !-- Store the SGS TKE of the current timelevel which is needed for … … 1507 1533 CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' ) 1508 1534 1535 1509 1536 END SUBROUTINE lpm_advec 1537 1538 ! Description: 1539 ! ------------ 1540 !> Calculation of subgrid-scale particle speed using the stochastic model 1541 !> of Weil et al. (2004, JAS, 61, 2877-2887). 1542 !------------------------------------------------------------------------------! 1543 SUBROUTINE weil_stochastic_eq( v_sgs, fs_n, e_n, dedxi_n, dedt_n, diss_n, & 1544 dt_n, rg_n, fac ) 1545 1546 USE kinds 1547 1548 USE particle_attributes, & 1549 ONLY: c_0, sgs_wf_part 1550 1551 IMPLICIT NONE 1552 1553 REAL(wp) :: a1 !< dummy argument 1554 REAL(wp) :: dedt_n !< time derivative of TKE at particle position 1555 REAL(wp) :: dedxi_n !< horizontal derivative of TKE at particle position 1556 REAL(wp) :: diss_n !< dissipation at particle position 1557 REAL(wp) :: dt_n !< particle timestep 1558 REAL(wp) :: e_n !< TKE at particle position 1559 REAL(wp) :: fac !< flag to identify adjacent topography 1560 REAL(wp) :: fs_n !< weighting factor to prevent that subgrid-scale particle speed becomes too large 1561 REAL(wp) :: sgs_w !< constant (1/3) 1562 REAL(wp) :: rg_n !< random number 1563 REAL(wp) :: term1 !< memory term 1564 REAL(wp) :: term2 !< drift correction term 1565 REAL(wp) :: term3 !< random term 1566 REAL(wp) :: v_sgs !< subgrid-scale velocity component 1567 1568 ! 1569 !-- Please note, terms 1 and 2 (drift and memory term, respectively) are 1570 !-- multiplied by a flag to switch of both terms near topography. 1571 !-- This is necessary, as both terms may cause a subgrid-scale velocity build up 1572 !-- if particles are trapped in regions with very small TKE, e.g. in narrow street 1573 !-- canyons resolved by only a few grid points. Hence, term 1 and term 2 are 1574 !-- disabled if one of the adjacent grid points belongs to topography. 1575 !-- Moreover, in this case, the previous subgrid-scale component is also set 1576 !-- to zero. 1577 1578 a1 = fs_n * c_0 * diss_n 1579 ! 1580 !-- Memory term 1581 term1 = - a1 * v_sgs * dt_n / ( 4.0_wp * sgs_wf_part * e_n + 1E-20_wp ) & 1582 * fac 1583 ! 1584 !-- Drift correction term 1585 term2 = ( ( dedt_n * v_sgs / e_n ) + dedxi_n ) * 0.5_wp * dt_n & 1586 * fac 1587 ! 1588 !-- Random term 1589 term3 = SQRT( MAX( a1, 1E-20 ) ) * ( rg_n - 1.0_wp ) * SQRT( dt_n ) 1590 ! 1591 !-- In cese one of the adjacent grid-boxes belongs to topograhy, the previous 1592 !-- subgrid-scale velocity component is set to zero, in order to prevent a 1593 !-- velocity build-up. 1594 1595 !-- This case, set also previous subgrid-scale component to zero. 1596 v_sgs = v_sgs * fac + term1 + term2 + term3 1597 1598 END SUBROUTINE weil_stochastic_eq -
palm/trunk/SOURCE/lpm_boundary_conds.f90
r1823 r1929 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Rewritten wall reflection 22 22 ! 23 23 ! Former revisions: … … 102 102 CHARACTER (LEN=*) :: range !< 103 103 104 INTEGER(iwp) :: i !< 105 INTEGER(iwp) :: inc !< 106 INTEGER(iwp) :: ir !< 107 INTEGER(iwp) :: i1 !< 108 INTEGER(iwp) :: i2 !< 109 INTEGER(iwp) :: i3 !< 110 INTEGER(iwp) :: i5 !< 111 INTEGER(iwp) :: j !< 112 INTEGER(iwp) :: jr !< 113 INTEGER(iwp) :: j1 !< 114 INTEGER(iwp) :: j2 !< 115 INTEGER(iwp) :: j3 !< 116 INTEGER(iwp) :: j5 !< 117 INTEGER(iwp) :: k !< 118 INTEGER(iwp) :: k1 !< 119 INTEGER(iwp) :: k2 !< 120 INTEGER(iwp) :: k3 !< 121 INTEGER(iwp) :: k5 !< 122 INTEGER(iwp) :: n !< 123 INTEGER(iwp) :: t_index !< 124 INTEGER(iwp) :: t_index_number !< 104 INTEGER(iwp) :: inc !< dummy for sorting algorithmus 105 INTEGER(iwp) :: ir !< dummy for sorting algorithmus 106 INTEGER(iwp) :: i1 !< grid index (x) of old particle position 107 INTEGER(iwp) :: i2 !< grid index (x) of current particle position 108 INTEGER(iwp) :: i3 !< grid index (x) of intermediate particle position 109 INTEGER(iwp) :: jr !< dummy for sorting algorithmus 110 INTEGER(iwp) :: j1 !< grid index (y) of old particle position 111 INTEGER(iwp) :: j2 !< grid index (x) of current particle position 112 INTEGER(iwp) :: j3 !< grid index (x) of intermediate particle position 113 INTEGER(iwp) :: n !< particle number 114 INTEGER(iwp) :: t_index !< running index for intermediate particle timesteps in reflection algorithmus 115 INTEGER(iwp) :: t_index_number !< number of intermediate particle timesteps in reflection algorithmus 116 INTEGER(iwp) :: tmp_x !< dummy for sorting algorithmus 117 INTEGER(iwp) :: tmp_y !< dummy for sorting algorithmus 118 119 INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions 120 INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (x) of intermediate particle positions 125 121 126 LOGICAL :: reflect_x !< 127 LOGICAL :: reflect_y !< 128 LOGICAL :: reflect_z !< 129 130 REAL(wp) :: dt_particle !< 131 REAL(wp) :: pos_x !< 132 REAL(wp) :: pos_x_old !< 133 REAL(wp) :: pos_y !< 134 REAL(wp) :: pos_y_old !< 135 REAL(wp) :: pos_z !< 136 REAL(wp) :: pos_z_old !< 137 REAL(wp) :: prt_x !< 138 REAL(wp) :: prt_y !< 139 REAL(wp) :: prt_z !< 140 REAL(wp) :: t(1:200) !< 141 REAL(wp) :: tmp_t !< 142 REAL(wp) :: xline !< 143 REAL(wp) :: yline !< 144 REAL(wp) :: zline !< 122 LOGICAL :: cross_wall_x !< flag to check if particle reflection along x is necessary 123 LOGICAL :: cross_wall_y !< flag to check if particle reflection along y is necessary 124 LOGICAL :: downwards !< flag to check if particle reflection along z is necessary (only if particle move downwards) 125 LOGICAL :: reflect_x !< flag to check if particle is already reflected along x 126 LOGICAL :: reflect_y !< flag to check if particle is already reflected along y 127 LOGICAL :: reflect_z !< flag to check if particle is already reflected along z 128 LOGICAL :: tmp_reach_x !< dummy for sorting algorithmus 129 LOGICAL :: tmp_reach_y !< dummy for sorting algorithmus 130 LOGICAL :: tmp_reach_z !< dummy for sorting algorithmus 131 LOGICAL :: x_wall_reached !< flag to check if particle has already reached wall 132 LOGICAL :: y_wall_reached !< flag to check if particle has already reached wall 133 134 LOGICAL, DIMENSION(0:10) :: reach_x !< flag to check if particle is at a yz-wall 135 LOGICAL, DIMENSION(0:10) :: reach_y !< flag to check if particle is at a xz-wall 136 LOGICAL, DIMENSION(0:10) :: reach_z !< flag to check if particle is at a xy-wall 137 138 REAL(wp) :: dt_particle !< particle timestep 139 REAL(wp) :: dum !< dummy argument 140 REAL(wp) :: eps = 1E-10_wp !< security number to check if particle has reached a wall 141 REAL(wp) :: pos_x !< intermediate particle position (x) 142 REAL(wp) :: pos_x_old !< particle position (x) at previous particle timestep 143 REAL(wp) :: pos_y !< intermediate particle position (y) 144 REAL(wp) :: pos_y_old !< particle position (y) at previous particle timestep 145 REAL(wp) :: pos_z !< intermediate particle position (z) 146 REAL(wp) :: pos_z_old !< particle position (z) at previous particle timestep 147 REAL(wp) :: prt_x !< current particle position (x) 148 REAL(wp) :: prt_y !< current particle position (y) 149 REAL(wp) :: prt_z !< current particle position (z) 150 REAL(wp) :: t_old !< previous reflection time 151 REAL(wp) :: tmp_t !< dummy for sorting algorithmus 152 REAL(wp) :: xwall !< location of wall in x 153 REAL(wp) :: ywall !< location of wall in y 154 REAL(wp) :: zwall1 !< location of wall in z (old grid box) 155 REAL(wp) :: zwall2 !< location of wall in z (current grid box) 156 REAL(wp) :: zwall3 !< location of wall in z (old y, current x) 157 REAL(wp) :: zwall4 !< location of wall in z (current y, old x) 158 159 REAL(wp), DIMENSION(0:10) :: t !< reflection time 160 145 161 146 162 IF ( range == 'bottom/top' ) THEN … … 211 227 ELSEIF ( range == 'walls' ) THEN 212 228 229 213 230 CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'start' ) 214 231 215 reflect_x = .FALSE.216 reflect_y = .FALSE.217 reflect_z = .FALSE.218 219 232 DO n = 1, number_of_particles 220 233 ! 234 !-- Recalculate particle timestep 221 235 dt_particle = particles(n)%age - particles(n)%age_m 222 236 ! 237 !-- Obtain x/y indices for current particle position 223 238 i2 = ( particles(n)%x + 0.5_wp * dx ) * ddx 224 239 j2 = ( particles(n)%y + 0.5_wp * dy ) * ddy 225 k2 = particles(n)%z / dz + 1 + offset_ocean_nzt_m1 226 240 ! 241 !-- Save current particle positions 227 242 prt_x = particles(n)%x 228 243 prt_y = particles(n)%y 229 244 prt_z = particles(n)%z 230 231 ! 232 !-- If the particle position is below the surface, it has to be reflected 233 IF ( k2 <= nzb_s_inner(j2,i2) .AND. nzb_s_inner(j2,i2) /=0 ) THEN 234 235 pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle 236 pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle 237 pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle 238 i1 = ( pos_x_old + 0.5_wp * dx ) * ddx 239 j1 = ( pos_y_old + 0.5_wp * dy ) * ddy 240 k1 = pos_z_old / dz + offset_ocean_nzt_m1 241 242 ! 243 !-- Case 1 244 IF ( particles(n)%x > pos_x_old .AND. particles(n)%y > pos_y_old )& 245 THEN 246 t_index = 1 247 248 DO i = i1, i2 249 xline = i * dx + 0.5_wp * dx 250 t(t_index) = ( xline - pos_x_old ) / & 251 ( particles(n)%x - pos_x_old ) 252 t_index = t_index + 1 245 ! 246 !-- Recalculate old particle positions 247 pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle 248 pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle 249 pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle 250 ! 251 !-- Obtain x/y indices for old particle positions 252 i1 = ( pos_x_old + 0.5_wp * dx ) * ddx 253 j1 = ( pos_y_old + 0.5_wp * dy ) * ddy 254 ! 255 !-- Determine horizontal as well as vertical walls at which particle can 256 !-- be potentially reflected. 257 !-- Start with walls aligned in yz layer. 258 !-- Wall to the right 259 IF ( prt_x > pos_x_old ) THEN 260 xwall = ( i1 + 0.5_wp ) * dx 261 ! 262 !-- Wall to the left 263 ELSE 264 xwall = ( i1 - 0.5_wp ) * dx 265 ENDIF 266 ! 267 !-- Walls aligned in xz layer 268 !-- Wall to the north 269 IF ( prt_y > pos_y_old ) THEN 270 ywall = ( j1 + 0.5_wp ) * dy 271 !-- Wall to the south 272 ELSE 273 ywall = ( j1 - 0.5_wp ) * dy 274 ENDIF 275 ! 276 !-- Walls aligned in xy layer at which particle can be possiblly reflected 277 zwall1 = zw(nzb_s_inner(j2,i2)) 278 zwall2 = zw(nzb_s_inner(j1,i1)) 279 zwall3 = zw(nzb_s_inner(j1,i2)) 280 zwall4 = zw(nzb_s_inner(j2,i1)) 281 ! 282 !-- Initialize flags to check if particle reflection is necessary 283 downwards = .FALSE. 284 cross_wall_x = .FALSE. 285 cross_wall_y = .FALSE. 286 ! 287 !-- Initialize flags to check if a wall is reached 288 reach_x = .FALSE. 289 reach_y = .FALSE. 290 reach_z = .FALSE. 291 ! 292 !-- Initialize flags to check if a particle was already reflected 293 reflect_x = .FALSE. 294 reflect_y = .FALSE. 295 reflect_z = .FALSE. 296 ! 297 !-- Initialize flags to check if a vertical wall is already crossed. 298 !-- ( Required to obtain correct indices. ) 299 x_wall_reached = .FALSE. 300 y_wall_reached = .FALSE. 301 ! 302 !-- Initialize time array 303 t = 0.0_wp 304 ! 305 !-- Check if particle can reach any wall. This case, calculate the 306 !-- fractional time needed to reach this wall. Store this fractional 307 !-- timestep in array t. Moreover, store indices for these grid 308 !-- boxes where the respective wall belongs to. 309 !-- Start with x-direction. 310 t_index = 1 311 t(t_index) = ( xwall - pos_x_old ) & 312 / MERGE( MAX( prt_x - pos_x_old, 1E-30_wp ), & 313 MIN( prt_x - pos_x_old, -1E-30_wp ), & 314 prt_x > pos_x_old ) 315 x_ind(t_index) = i2 316 y_ind(t_index) = j1 317 reach_x(t_index) = .TRUE. 318 reach_y(t_index) = .FALSE. 319 reach_z(t_index) = .FALSE. 320 ! 321 !-- Store these values only if particle really reaches any wall. t must 322 !-- be in a interval between [0:1]. 323 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) THEN 324 t_index = t_index + 1 325 cross_wall_x = .TRUE. 326 ENDIF 327 ! 328 !-- y-direction 329 t(t_index) = ( ywall - pos_y_old ) & 330 / MERGE( MAX( prt_y - pos_y_old, 1E-30_wp ), & 331 MIN( prt_y - pos_y_old, -1E-30_wp ), & 332 prt_y > pos_y_old ) 333 x_ind(t_index) = i1 334 y_ind(t_index) = j2 335 reach_x(t_index) = .FALSE. 336 reach_y(t_index) = .TRUE. 337 reach_z(t_index) = .FALSE. 338 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) THEN 339 t_index = t_index + 1 340 cross_wall_y = .TRUE. 341 ENDIF 342 ! 343 !-- z-direction 344 !-- At first, check if particle moves downwards. Only in this case a 345 !-- particle can be reflected vertically. 346 IF ( prt_z < pos_z_old ) THEN 347 348 downwards = .TRUE. 349 dum = 1.0_wp / MERGE( MAX( prt_z - pos_z_old, 1E-30_wp ), & 350 MIN( prt_z - pos_z_old, -1E-30_wp ), & 351 prt_z > pos_z_old ) 352 353 t(t_index) = ( zwall1 - pos_z_old ) * dum 354 x_ind(t_index) = i2 355 y_ind(t_index) = j2 356 reach_x(t_index) = .FALSE. 357 reach_y(t_index) = .FALSE. 358 reach_z(t_index) = .TRUE. 359 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) & 360 t_index = t_index + 1 361 362 reach_x(t_index) = .FALSE. 363 reach_y(t_index) = .FALSE. 364 reach_z(t_index) = .TRUE. 365 t(t_index) = ( zwall2 - pos_z_old ) * dum 366 x_ind(t_index) = i1 367 y_ind(t_index) = j1 368 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) & 369 t_index = t_index + 1 370 371 reach_x(t_index) = .FALSE. 372 reach_y(t_index) = .FALSE. 373 reach_z(t_index) = .TRUE. 374 t(t_index) = ( zwall3 - pos_z_old ) * dum 375 x_ind(t_index) = i2 376 y_ind(t_index) = j1 377 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) & 378 t_index = t_index + 1 379 380 reach_x(t_index) = .FALSE. 381 reach_y(t_index) = .FALSE. 382 reach_z(t_index) = .TRUE. 383 t(t_index) = ( zwall4 - pos_z_old ) * dum 384 x_ind(t_index) = i1 385 y_ind(t_index) = j2 386 IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp ) & 387 t_index = t_index + 1 388 389 END IF 390 t_index_number = t_index - 1 391 ! 392 !-- Carry out reflection only if particle reaches any wall 393 IF ( cross_wall_x .OR. cross_wall_y .OR. downwards ) THEN 394 ! 395 !-- Sort fractional timesteps in ascending order. Also sort the 396 !-- corresponding indices and flag according to the time interval a 397 !-- particle reaches the respective wall. 398 inc = 1 399 jr = 1 400 DO WHILE ( inc <= t_index_number ) 401 inc = 3 * inc + 1 402 ENDDO 403 404 DO WHILE ( inc > 1 ) 405 inc = inc / 3 406 DO ir = inc+1, t_index_number 407 tmp_t = t(ir) 408 tmp_x = x_ind(ir) 409 tmp_y = y_ind(ir) 410 tmp_reach_x = reach_x(ir) 411 tmp_reach_y = reach_y(ir) 412 tmp_reach_z = reach_z(ir) 413 jr = ir 414 DO WHILE ( t(jr-inc) > tmp_t ) 415 t(jr) = t(jr-inc) 416 x_ind(jr) = x_ind(jr-inc) 417 y_ind(jr) = y_ind(jr-inc) 418 reach_x(jr) = reach_x(jr-inc) 419 reach_y(jr) = reach_y(jr-inc) 420 reach_z(jr) = reach_z(jr-inc) 421 jr = jr - inc 422 IF ( jr <= inc ) EXIT 423 ENDDO 424 t(jr) = tmp_t 425 x_ind(jr) = tmp_x 426 y_ind(jr) = tmp_y 427 reach_x(jr) = tmp_reach_x 428 reach_y(jr) = tmp_reach_y 429 reach_z(jr) = tmp_reach_z 253 430 ENDDO 254 255 DO j = j1, j2 256 yline = j * dy + 0.5_wp * dy 257 t(t_index) = ( yline - pos_y_old ) / & 258 ( particles(n)%y - pos_y_old ) 259 t_index = t_index + 1 260 ENDDO 261 262 IF ( particles(n)%z < pos_z_old ) THEN 263 DO k = k1, k2 , -1 264 zline = k * dz 265 t(t_index) = ( pos_z_old - zline ) / & 266 ( pos_z_old - particles(n)%z ) 267 t_index = t_index + 1 268 ENDDO 431 ENDDO 432 ! 433 !-- Initialize temporary particle positions 434 pos_x = pos_x_old 435 pos_y = pos_y_old 436 pos_z = pos_z_old 437 ! 438 !-- Loop over all times a particle possibly moves into a new grid box 439 t_old = 0.0_wp 440 DO t_index = 1, t_index_number 441 ! 442 !-- Calculate intermediate particle position according to the 443 !-- timesteps a particle reaches any wall. 444 pos_x = pos_x + ( t(t_index) - t_old ) * dt_particle & 445 * particles(n)%speed_x 446 pos_y = pos_y + ( t(t_index) - t_old ) * dt_particle & 447 * particles(n)%speed_y 448 pos_z = pos_z + ( t(t_index) - t_old ) * dt_particle & 449 * particles(n)%speed_z 450 ! 451 !-- Obtain x/y grid indices for intermediate particle position from 452 !-- sorted index array 453 i3 = x_ind(t_index) 454 j3 = y_ind(t_index) 455 ! 456 !-- Check which wall is already reached 457 IF ( .NOT. x_wall_reached ) x_wall_reached = reach_x(t_index) 458 IF ( .NOT. y_wall_reached ) y_wall_reached = reach_y(t_index) 459 ! 460 !-- Check if a particle needs to be reflected at any yz-wall. If 461 !-- necessary, carry out reflection. Please note, a security 462 !-- constant is required, as the particle position do not 463 !-- necessarily exactly match the wall location due to rounding 464 !-- errors. 465 IF ( ABS( pos_x - xwall ) < eps .AND. & 466 pos_z <= zw(nzb_s_inner(j3,i3)) .AND. & 467 reach_x(t_index) .AND. & 468 .NOT. reflect_x ) THEN 469 ! 470 !-- Reflection in x-direction. 471 !-- Ensure correct reflection by MIN/MAX functions, depending on 472 !-- direction of particle transport. 473 !-- Due to rounding errors pos_x do not exactly matches the wall 474 !-- location, leading to erroneous reflection. 475 pos_x = MERGE( MIN( 2.0_wp * xwall - pos_x, xwall ), & 476 MAX( 2.0_wp * xwall - pos_x, xwall ), & 477 particles(n)%x > xwall ) 478 ! 479 !-- Change sign of particle speed 480 particles(n)%speed_x = - particles(n)%speed_x 481 ! 482 !-- Change also sign of subgrid-scale particle speed 483 particles(n)%rvar1 = - particles(n)%rvar1 484 ! 485 !-- Set flag that reflection along x is already done 486 reflect_x = .TRUE. 487 ! 488 !-- As particle do not crosses any further yz-wall during 489 !-- this timestep, set further x-indices to the current one. 490 x_ind(t_index:t_index_number) = i1 491 ! 492 !-- If particle already reached the wall but was not reflected, 493 !-- set further x-indices to the new one. 494 ELSEIF ( x_wall_reached .AND. .NOT. reflect_x ) THEN 495 x_ind(t_index:t_index_number) = i2 269 496 ENDIF 270 271 t_index_number = t_index - 1 272 273 ! 274 !-- Sorting t 275 inc = 1 276 jr = 1 277 DO WHILE ( inc <= t_index_number ) 278 inc = 3 * inc + 1 279 ENDDO 280 281 DO WHILE ( inc > 1 ) 282 inc = inc / 3 283 DO ir = inc+1, t_index_number 284 tmp_t = t(ir) 285 jr = ir 286 DO WHILE ( t(jr-inc) > tmp_t ) 287 t(jr) = t(jr-inc) 288 jr = jr - inc 289 IF ( jr <= inc ) EXIT 290 ENDDO 291 t(jr) = tmp_t 292 ENDDO 293 ENDDO 294 295 case1: DO t_index = 1, t_index_number 296 297 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old ) 298 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old ) 299 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old ) 300 301 i3 = ( pos_x + 0.5_wp * dx ) * ddx 302 j3 = ( pos_y + 0.5_wp * dy ) * ddy 303 k3 = pos_z / dz + offset_ocean_nzt_m1 304 305 i5 = pos_x * ddx 306 j5 = pos_y * ddy 307 k5 = pos_z / dz + offset_ocean_nzt_m1 308 309 IF ( k5 <= nzb_s_inner(j5,i3) .AND. & 310 nzb_s_inner(j5,i3) /= 0 ) THEN 311 312 IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. & 313 k3 == nzb_s_inner(j5,i3) ) THEN 314 reflect_z = .TRUE. 315 EXIT case1 316 ENDIF 497 ! 498 !-- Check if a particle needs to be reflected at any xz-wall. If 499 !-- necessary, carry out reflection. 500 IF ( ABS( pos_y - ywall ) < eps .AND. & 501 pos_z <= zw(nzb_s_inner(j3,i3)) .AND. & 502 reach_y(t_index) .AND. & 503 .NOT. reflect_y ) THEN 504 505 pos_y = MERGE( MIN( 2.0_wp * ywall - pos_y, ywall ), & 506 MAX( 2.0_wp * ywall - pos_y, ywall ), & 507 particles(n)%y > ywall ) 508 509 particles(n)%speed_y = - particles(n)%speed_y 510 particles(n)%rvar2 = - particles(n)%rvar2 511 512 reflect_y = .TRUE. 513 y_ind(t_index:t_index_number) = j1 514 515 ELSEIF ( y_wall_reached .AND. .NOT. reflect_y ) THEN 516 y_ind(t_index:t_index_number) = j2 517 ENDIF 518 ! 519 !-- Check if a particle needs to be reflected at any xy-wall. If 520 !-- necessary, carry out reflection. 521 IF ( downwards .AND. reach_z(t_index) .AND. & 522 .NOT. reflect_z ) THEN 523 IF ( pos_z - zw(nzb_s_inner(j3,i3)) < eps ) THEN 524 525 pos_z = MAX( 2.0_wp * zw(nzb_s_inner(j3,i3)) - pos_z, & 526 zw(nzb_s_inner(j3,i3)) ) 527 528 particles(n)%speed_z = - particles(n)%speed_z 529 particles(n)%rvar3 = - particles(n)%rvar3 530 531 reflect_z = .TRUE. 317 532 318 533 ENDIF 319 534 320 IF ( k5 <= nzb_s_inner(j3,i5) .AND. &321 nzb_s_inner(j3,i5) /= 0 ) THEN322 323 IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. &324 k3 == nzb_s_inner(j3,i5) ) THEN325 reflect_z = .TRUE.326 EXIT case1327 ENDIF328 329 ENDIF330 331 IF ( k3 <= nzb_s_inner(j3,i3) .AND. &332 nzb_s_inner(j3,i3) /= 0 ) THEN333 334 IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. &335 k3 == nzb_s_inner(j3,i3) ) THEN336 reflect_z = .TRUE.337 EXIT case1338 ENDIF339 340 IF ( pos_y == ( j3 * dy - 0.5_wp * dy ) .AND. &341 pos_z < nzb_s_inner(j3,i3) * dz ) THEN342 reflect_y = .TRUE.343 EXIT case1344 ENDIF345 346 IF ( pos_x == ( i3 * dx - 0.5_wp * dx ) .AND. &347 pos_z < nzb_s_inner(j3,i3) * dz ) THEN348 reflect_x = .TRUE.349 EXIT case1350 ENDIF351 352 ENDIF353 354 ENDDO case1355 356 !357 !-- Case 2358 ELSEIF ( particles(n)%x > pos_x_old .AND. &359 particles(n)%y < pos_y_old ) THEN360 361 t_index = 1362 363 DO i = i1, i2364 xline = i * dx + 0.5_wp * dx365 t(t_index) = ( xline - pos_x_old ) / &366 ( particles(n)%x - pos_x_old )367 t_index = t_index + 1368 ENDDO369 370 DO j = j1, j2, -1371 yline = j * dy - 0.5_wp * dy372 t(t_index) = ( pos_y_old - yline ) / &373 ( pos_y_old - particles(n)%y )374 t_index = t_index + 1375 ENDDO376 377 IF ( particles(n)%z < pos_z_old ) THEN378 DO k = k1, k2 , -1379 zline = k * dz380 t(t_index) = ( pos_z_old - zline ) / &381 ( pos_z_old - particles(n)%z )382 t_index = t_index + 1383 ENDDO384 535 ENDIF 385 t_index_number = t_index-1 386 387 ! 388 !-- Sorting t 389 inc = 1 390 jr = 1 391 DO WHILE ( inc <= t_index_number ) 392 inc = 3 * inc + 1 393 ENDDO 394 395 DO WHILE ( inc > 1 ) 396 inc = inc / 3 397 DO ir = inc+1, t_index_number 398 tmp_t = t(ir) 399 jr = ir 400 DO WHILE ( t(jr-inc) > tmp_t ) 401 t(jr) = t(jr-inc) 402 jr = jr - inc 403 IF ( jr <= inc ) EXIT 404 ENDDO 405 t(jr) = tmp_t 406 ENDDO 407 ENDDO 408 409 case2: DO t_index = 1, t_index_number 410 411 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old ) 412 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old ) 413 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old ) 414 415 i3 = ( pos_x + 0.5_wp * dx ) * ddx 416 j3 = ( pos_y + 0.5_wp * dy ) * ddy 417 k3 = pos_z / dz + offset_ocean_nzt_m1 418 419 i5 = pos_x * ddx 420 j5 = pos_y * ddy 421 k5 = pos_z / dz + offset_ocean_nzt_m1 422 423 IF ( k5 <= nzb_s_inner(j3,i5) .AND. & 424 nzb_s_inner(j3,i5) /= 0 ) THEN 425 426 IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. & 427 k3 == nzb_s_inner(j3,i5) ) THEN 428 reflect_z = .TRUE. 429 EXIT case2 430 ENDIF 431 432 ENDIF 433 434 IF ( k3 <= nzb_s_inner(j3,i3) .AND. & 435 nzb_s_inner(j3,i3) /= 0 ) THEN 436 437 IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. & 438 k3 == nzb_s_inner(j3,i3) ) THEN 439 reflect_z = .TRUE. 440 EXIT case2 441 ENDIF 442 443 IF ( pos_x == ( i3 * dx - 0.5_wp * dx ) .AND. & 444 pos_z < nzb_s_inner(j3,i3) * dz ) THEN 445 reflect_x = .TRUE. 446 EXIT case2 447 ENDIF 448 449 ENDIF 450 451 IF ( k5 <= nzb_s_inner(j5,i3) .AND. & 452 nzb_s_inner(j5,i3) /= 0 ) THEN 453 454 IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. & 455 k3 == nzb_s_inner(j5,i3) ) THEN 456 reflect_z = .TRUE. 457 EXIT case2 458 ENDIF 459 460 IF ( pos_y == ( j5 * dy + 0.5_wp * dy ) .AND. & 461 pos_z < nzb_s_inner(j5,i3) * dz ) THEN 462 reflect_y = .TRUE. 463 EXIT case2 464 ENDIF 465 466 ENDIF 467 468 ENDDO case2 469 470 ! 471 !-- Case 3 472 ELSEIF ( particles(n)%x < pos_x_old .AND. & 473 particles(n)%y > pos_y_old ) THEN 474 475 t_index = 1 476 477 DO i = i1, i2, -1 478 xline = i * dx - 0.5_wp * dx 479 t(t_index) = ( pos_x_old - xline ) / & 480 ( pos_x_old - particles(n)%x ) 481 t_index = t_index + 1 482 ENDDO 483 484 DO j = j1, j2 485 yline = j * dy + 0.5_wp * dy 486 t(t_index) = ( yline - pos_y_old ) / & 487 ( particles(n)%y - pos_y_old ) 488 t_index = t_index + 1 489 ENDDO 490 491 IF ( particles(n)%z < pos_z_old ) THEN 492 DO k = k1, k2 , -1 493 zline = k * dz 494 t(t_index) = ( pos_z_old - zline ) / & 495 ( pos_z_old - particles(n)%z ) 496 t_index = t_index + 1 497 ENDDO 498 ENDIF 499 t_index_number = t_index - 1 500 501 ! 502 !-- Sorting t 503 inc = 1 504 jr = 1 505 506 DO WHILE ( inc <= t_index_number ) 507 inc = 3 * inc + 1 508 ENDDO 509 510 DO WHILE ( inc > 1 ) 511 inc = inc / 3 512 DO ir = inc+1, t_index_number 513 tmp_t = t(ir) 514 jr = ir 515 DO WHILE ( t(jr-inc) > tmp_t ) 516 t(jr) = t(jr-inc) 517 jr = jr - inc 518 IF ( jr <= inc ) EXIT 519 ENDDO 520 t(jr) = tmp_t 521 ENDDO 522 ENDDO 523 524 case3: DO t_index = 1, t_index_number 525 526 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old ) 527 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old ) 528 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old ) 529 530 i3 = ( pos_x + 0.5_wp * dx ) * ddx 531 j3 = ( pos_y + 0.5_wp * dy ) * ddy 532 k3 = pos_z / dz + offset_ocean_nzt_m1 533 534 i5 = pos_x * ddx 535 j5 = pos_y * ddy 536 k5 = pos_z / dz + offset_ocean_nzt_m1 537 538 IF ( k5 <= nzb_s_inner(j5,i3) .AND. & 539 nzb_s_inner(j5,i3) /= 0 ) THEN 540 541 IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. & 542 k3 == nzb_s_inner(j5,i3) ) THEN 543 reflect_z = .TRUE. 544 EXIT case3 545 ENDIF 546 547 ENDIF 548 549 IF ( k3 <= nzb_s_inner(j3,i3) .AND. & 550 nzb_s_inner(j3,i3) /= 0 ) THEN 551 552 IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. & 553 k3 == nzb_s_inner(j3,i3) ) THEN 554 reflect_z = .TRUE. 555 EXIT case3 556 ENDIF 557 558 IF ( pos_y == ( j3 * dy - 0.5_wp * dy ) .AND. & 559 pos_z < nzb_s_inner(j3,i3) * dz ) THEN 560 reflect_y = .TRUE. 561 EXIT case3 562 ENDIF 563 564 ENDIF 565 566 IF ( k5 <= nzb_s_inner(j3,i5) .AND. & 567 nzb_s_inner(j3,i5) /= 0 ) THEN 568 569 IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. & 570 k3 == nzb_s_inner(j3,i5) ) THEN 571 reflect_z = .TRUE. 572 EXIT case3 573 ENDIF 574 575 IF ( pos_x == ( i5 * dx + 0.5_wp * dx ) .AND. & 576 pos_z < nzb_s_inner(j3,i5) * dz ) THEN 577 reflect_x = .TRUE. 578 EXIT case3 579 ENDIF 580 581 ENDIF 582 583 ENDDO case3 584 585 ! 586 !-- Case 4 587 ELSEIF ( particles(n)%x < pos_x_old .AND. & 588 particles(n)%y < pos_y_old ) THEN 589 590 t_index = 1 591 592 DO i = i1, i2, -1 593 xline = i * dx - 0.5_wp * dx 594 t(t_index) = ( pos_x_old - xline ) / & 595 ( pos_x_old - particles(n)%x ) 596 t_index = t_index + 1 597 ENDDO 598 599 DO j = j1, j2, -1 600 yline = j * dy - 0.5_wp * dy 601 t(t_index) = ( pos_y_old - yline ) / & 602 ( pos_y_old - particles(n)%y ) 603 t_index = t_index + 1 604 ENDDO 605 606 IF ( particles(n)%z < pos_z_old ) THEN 607 DO k = k1, k2 , -1 608 zline = k * dz 609 t(t_index) = ( pos_z_old - zline ) / & 610 ( pos_z_old-particles(n)%z ) 611 t_index = t_index + 1 612 ENDDO 613 ENDIF 614 t_index_number = t_index-1 615 616 ! 617 !-- Sorting t 618 inc = 1 619 jr = 1 620 621 DO WHILE ( inc <= t_index_number ) 622 inc = 3 * inc + 1 623 ENDDO 624 625 DO WHILE ( inc > 1 ) 626 inc = inc / 3 627 DO ir = inc+1, t_index_number 628 tmp_t = t(ir) 629 jr = ir 630 DO WHILE ( t(jr-inc) > tmp_t ) 631 t(jr) = t(jr-inc) 632 jr = jr - inc 633 IF ( jr <= inc ) EXIT 634 ENDDO 635 t(jr) = tmp_t 636 ENDDO 637 ENDDO 638 639 case4: DO t_index = 1, t_index_number 640 641 pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old ) 642 pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old ) 643 pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old ) 644 645 i3 = ( pos_x + 0.5_wp * dx ) * ddx 646 j3 = ( pos_y + 0.5_wp * dy ) * ddy 647 k3 = pos_z / dz + offset_ocean_nzt_m1 648 649 i5 = pos_x * ddx 650 j5 = pos_y * ddy 651 k5 = pos_z / dz + offset_ocean_nzt_m1 652 653 IF ( k3 <= nzb_s_inner(j3,i3) .AND. & 654 nzb_s_inner(j3,i3) /= 0 ) THEN 655 656 IF ( pos_z == nzb_s_inner(j3,i3) * dz .AND. & 657 k3 == nzb_s_inner(j3,i3) ) THEN 658 reflect_z = .TRUE. 659 EXIT case4 660 ENDIF 661 662 ENDIF 663 664 IF ( k5 <= nzb_s_inner(j3,i5) .AND. & 665 nzb_s_inner(j3,i5) /= 0 ) THEN 666 667 IF ( pos_z == nzb_s_inner(j3,i5) * dz .AND. & 668 k3 == nzb_s_inner(j3,i5) ) THEN 669 reflect_z = .TRUE. 670 EXIT case4 671 ENDIF 672 673 IF ( pos_x == ( i5 * dx + 0.5_wp * dx ) .AND. & 674 nzb_s_inner(j3,i5) /=0 .AND. & 675 pos_z < nzb_s_inner(j3,i5) * dz ) THEN 676 reflect_x = .TRUE. 677 EXIT case4 678 ENDIF 679 680 ENDIF 681 682 IF ( k5 <= nzb_s_inner(j5,i3) .AND. & 683 nzb_s_inner(j5,i3) /= 0 ) THEN 684 685 IF ( pos_z == nzb_s_inner(j5,i3) * dz .AND. & 686 k5 == nzb_s_inner(j5,i3) ) THEN 687 reflect_z = .TRUE. 688 EXIT case4 689 ENDIF 690 691 IF ( pos_y == ( j5 * dy + 0.5_wp * dy ) .AND. & 692 nzb_s_inner(j5,i3) /= 0 .AND. & 693 pos_z < nzb_s_inner(j5,i3) * dz ) THEN 694 reflect_y = .TRUE. 695 EXIT case4 696 ENDIF 697 698 ENDIF 699 700 ENDDO case4 536 ! 537 !-- Swap time 538 t_old = t(t_index) 539 540 ENDDO 541 ! 542 !-- If a particle was reflected, calculate final position from last 543 !-- intermediate position. 544 IF ( reflect_x .OR. reflect_y .OR. reflect_z ) THEN 545 546 particles(n)%x = pos_x + ( 1.0_wp - t_old ) * dt_particle & 547 * particles(n)%speed_x 548 particles(n)%y = pos_y + ( 1.0_wp - t_old ) * dt_particle & 549 * particles(n)%speed_y 550 particles(n)%z = pos_z + ( 1.0_wp - t_old ) * dt_particle & 551 * particles(n)%speed_z 701 552 702 553 ENDIF 703 554 704 ENDIF ! Check, if particle position is below the surface 705 706 ! 707 !-- Do the respective reflection, in case that one of the above conditions 708 !-- is found to be true 709 IF ( reflect_z ) THEN 710 711 particles(n)%z = 2.0_wp * pos_z - prt_z 712 particles(n)%speed_z = - particles(n)%speed_z 713 714 IF ( use_sgs_for_particles ) THEN 715 particles(n)%rvar3 = - particles(n)%rvar3 716 ENDIF 717 reflect_z = .FALSE. 718 719 ELSEIF ( reflect_y ) THEN 720 721 particles(n)%y = 2.0_wp * pos_y - prt_y 722 particles(n)%speed_y = - particles(n)%speed_y 723 724 IF ( use_sgs_for_particles ) THEN 725 particles(n)%rvar2 = - particles(n)%rvar2 726 ENDIF 727 reflect_y = .FALSE. 728 729 ELSEIF ( reflect_x ) THEN 730 731 particles(n)%x = 2.0_wp * pos_x - prt_x 732 particles(n)%speed_x = - particles(n)%speed_x 733 734 IF ( use_sgs_for_particles ) THEN 735 particles(n)%rvar1 = - particles(n)%rvar1 736 ENDIF 737 reflect_x = .FALSE. 738 739 ENDIF 740 555 ENDIF 556 741 557 ENDDO 742 558 -
palm/trunk/SOURCE/lpm_exchange_horiz.f90
r1874 r1929 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! Bugfixes: 22 ! - reallocation of new particles 23 ! ( did not work for small number of min_nr_particle ) 24 ! - dynamical reallocation of north-south exchange arrays ( particles got lost ) 25 ! - north-south exchange ( nr_move_north, nr_move_south were overwritten by zero ) 26 ! - horizontal particle boundary conditions in serial mode 27 ! 28 ! Remove unused variables 29 ! Descriptions in variable declaration blocks added 22 30 ! 23 31 ! Former revisions: … … 121 129 INTEGER(iwp) :: nr_move_south !< 122 130 123 TYPE(particle_type), DIMENSION( NR_2_direction_move):: move_also_north124 TYPE(particle_type), DIMENSION( NR_2_direction_move):: move_also_south131 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: move_also_north 132 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: move_also_south 125 133 126 134 SAVE … … 161 169 IMPLICIT NONE 162 170 163 INTEGER(iwp) :: i !< 164 INTEGER(iwp) :: ip !< 165 INTEGER(iwp) :: j !< 166 INTEGER(iwp) :: jp !< 167 INTEGER(iwp) :: kp !< 168 INTEGER(iwp) :: n !< 169 INTEGER(iwp) :: trlp_count !< 170 INTEGER(iwp) :: trlp_count_recv !< 171 INTEGER(iwp) :: trlpt_count !< 172 INTEGER(iwp) :: trlpt_count_recv !< 173 INTEGER(iwp) :: trnp_count !< 174 INTEGER(iwp) :: trnp_count_recv !< 175 INTEGER(iwp) :: trnpt_count !< 176 INTEGER(iwp) :: trnpt_count_recv !< 177 INTEGER(iwp) :: trrp_count !< 178 INTEGER(iwp) :: trrp_count_recv !< 179 INTEGER(iwp) :: trrpt_count !< 180 INTEGER(iwp) :: trrpt_count_recv !< 181 INTEGER(iwp) :: trsp_count !< 182 INTEGER(iwp) :: trsp_count_recv !< 183 INTEGER(iwp) :: trspt_count !< 184 INTEGER(iwp) :: trspt_count_recv !< 185 186 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: rvlp !< 187 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: rvnp !< 188 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: rvrp !< 189 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: rvsp !< 190 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: trlp !< 191 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: trnp !< 192 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: trrp !< 193 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: trsp !< 171 INTEGER(iwp) :: i !< grid index (x) of particle positition 172 INTEGER(iwp) :: ip !< index variable along x 173 INTEGER(iwp) :: j !< grid index (y) of particle positition 174 INTEGER(iwp) :: jp !< index variable along y 175 INTEGER(iwp) :: kp !< index variable along z 176 INTEGER(iwp) :: n !< particle index variable 177 INTEGER(iwp) :: trlp_count !< number of particles send to left PE 178 INTEGER(iwp) :: trlp_count_recv !< number of particles receive from right PE 179 INTEGER(iwp) :: trnp_count !< number of particles send to north PE 180 INTEGER(iwp) :: trnp_count_recv !< number of particles receive from south PE 181 INTEGER(iwp) :: trrp_count !< number of particles send to right PE 182 INTEGER(iwp) :: trrp_count_recv !< number of particles receive from left PE 183 INTEGER(iwp) :: trsp_count !< number of particles send to south PE 184 INTEGER(iwp) :: trsp_count_recv !< number of particles receive from north PE 185 186 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: rvlp !< particles received from right PE 187 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: rvnp !< particles received from south PE 188 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: rvrp !< particles received from left PE 189 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: rvsp !< particles received from north PE 190 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: trlp !< particles send to left PE 191 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: trnp !< particles send to north PE 192 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: trrp !< particles send to right PE 193 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: trsp !< particles send to south PE 194 194 195 195 CALL cpu_log( log_point_s(23), 'lpm_exchange_horiz', 'start' ) … … 209 209 !-- be adjusted. 210 210 trlp_count = 0 211 trlpt_count = 0212 211 trrp_count = 0 213 trrpt_count = 0214 212 215 213 trlp_count_recv = 0 216 trlpt_count_recv = 0217 214 trrp_count_recv = 0 218 trrpt_count_recv = 0219 215 220 216 IF ( pdims(1) /= 1 ) THEN … … 232 228 IF ( particles(n)%particle_mask ) THEN 233 229 i = ( particles(n)%x + 0.5_wp * dx ) * ddx 234 235 !--Above calculation does not work for indices less than zero230 ! 231 !-- Above calculation does not work for indices less than zero 236 232 IF ( particles(n)%x < -0.5_wp * dx ) i = -1 237 233 … … 249 245 250 246 IF ( trlp_count == 0 ) trlp_count = 1 251 IF ( trlpt_count == 0 ) trlpt_count = 1252 247 IF ( trrp_count == 0 ) trrp_count = 1 253 IF ( trrpt_count == 0 ) trrpt_count = 1254 248 255 249 ALLOCATE( trlp(trlp_count), trrp(trrp_count) ) … … 259 253 260 254 trlp_count = 0 261 trlpt_count = 0262 255 trrp_count = 0 263 trrpt_count = 0264 256 265 257 ENDIF 266 258 ! 267 259 !-- Compute only first (nxl) and last (nxr) loop iterration 268 DO ip = nxl, nxr, nxr-nxl260 DO ip = nxl, nxr, nxr-nxl 269 261 DO jp = nys, nyn 270 262 DO kp = nzb+1, nzt … … 279 271 280 272 i = ( particles(n)%x + 0.5_wp * dx ) * ddx 281 282 !--Above calculation does not work for indices less than zero273 ! 274 !-- Above calculation does not work for indices less than zero 283 275 IF ( particles(n)%x < - 0.5_wp * dx ) i = -1 284 276 285 277 IF ( i < nxl ) THEN 286 278 IF ( i < 0 ) THEN 287 288 279 ! 280 !-- Apply boundary condition along x 289 281 IF ( ibc_par_lr == 0 ) THEN 290 291 !--Cyclic condition282 ! 283 !-- Cyclic condition 292 284 IF ( pdims(1) == 1 ) THEN 293 285 particles(n)%x = ( nx + 1 ) * dx + particles(n)%x … … 312 304 313 305 ELSEIF ( ibc_par_lr == 1 ) THEN 314 315 !--Particle absorption306 ! 307 !-- Particle absorption 316 308 particles(n)%particle_mask = .FALSE. 317 309 deleted_particles = deleted_particles + 1 318 310 319 311 ELSEIF ( ibc_par_lr == 2 ) THEN 320 321 !--Particle reflection312 ! 313 !-- Particle reflection 322 314 particles(n)%x = -particles(n)%x 323 315 particles(n)%speed_x = -particles(n)%speed_x … … 325 317 ENDIF 326 318 ELSE 327 328 !--Store particle data in the transfer array, which will be329 !--send to the neighbouring PE319 ! 320 !-- Store particle data in the transfer array, which will be 321 !-- send to the neighbouring PE 330 322 trlp_count = trlp_count + 1 331 323 trlp(trlp_count) = particles(n) … … 337 329 ELSEIF ( i > nxr ) THEN 338 330 IF ( i > nx ) THEN 339 340 !--Apply boundary condition along x331 ! 332 !-- Apply boundary condition along x 341 333 IF ( ibc_par_lr == 0 ) THEN 342 343 !--Cyclic condition334 ! 335 !-- Cyclic condition 344 336 IF ( pdims(1) == 1 ) THEN 345 337 particles(n)%x = particles(n)%x - ( nx + 1 ) * dx … … 358 350 359 351 ELSEIF ( ibc_par_lr == 1 ) THEN 360 361 !--Particle absorption352 ! 353 !-- Particle absorption 362 354 particles(n)%particle_mask = .FALSE. 363 355 deleted_particles = deleted_particles + 1 364 356 365 357 ELSEIF ( ibc_par_lr == 2 ) THEN 366 367 !--Particle reflection358 ! 359 !-- Particle reflection 368 360 particles(n)%x = 2 * ( nx * dx ) - particles(n)%x 369 361 particles(n)%speed_x = -particles(n)%speed_x … … 371 363 ENDIF 372 364 ELSE 373 374 !--Store particle data in the transfer array, which will be send375 !--to the neighbouring PE365 ! 366 !-- Store particle data in the transfer array, which will be send 367 !-- to the neighbouring PE 376 368 trrp_count = trrp_count + 1 377 369 trrp(trrp_count) = particles(n) … … 383 375 ENDIF 384 376 ENDIF 377 385 378 ENDDO 386 379 ENDDO … … 389 382 390 383 ! 384 !-- Allocate arrays required for north-south exchange, as these 385 !-- are used directly after particles are exchange along x-direction. 386 ALLOCATE( move_also_north(1:NR_2_direction_move) ) 387 ALLOCATE( move_also_south(1:NR_2_direction_move) ) 388 389 nr_move_north = 0 390 nr_move_south = 0 391 ! 391 392 !-- Send left boundary, receive right boundary (but first exchange how many 392 393 !-- and check, if particle storage must be extended) … … 427 428 428 429 ENDIF 429 430 430 431 431 ! … … 435 435 !-- Find out first the number of particles to be transferred and allocate 436 436 !-- temporary arrays needed to store them. 437 !-- For a one-dimensional decomposition along x, no transfer is necessary,437 !-- For a one-dimensional decomposition along y, no transfer is necessary, 438 438 !-- because the particle remains on the PE. 439 439 trsp_count = nr_move_south 440 trspt_count = 0441 440 trnp_count = nr_move_north 442 trnpt_count = 0443 441 444 442 trsp_count_recv = 0 445 trspt_count_recv = 0446 443 trnp_count_recv = 0 447 trnpt_count_recv = 0448 444 449 445 IF ( pdims(2) /= 1 ) THEN … … 452 448 !-- data 453 449 DO ip = nxl, nxr 454 DO jp = nys, nyn, nyn-nys 450 DO jp = nys, nyn, nyn-nys !compute only first (nys) and last (nyn) loop iterration 455 451 DO kp = nzb+1, nzt 456 452 number_of_particles = prt_count(kp,jp,ip) … … 476 472 477 473 IF ( trsp_count == 0 ) trsp_count = 1 478 IF ( trspt_count == 0 ) trspt_count = 1479 474 IF ( trnp_count == 0 ) trnp_count = 1 480 IF ( trnpt_count == 0 ) trnpt_count = 1481 475 482 476 ALLOCATE( trsp(trsp_count), trnp(trnp_count) ) … … 486 480 487 481 trsp_count = nr_move_south 488 trspt_count = 0489 482 trnp_count = nr_move_north 490 trnpt_count = 0 491 483 492 484 trsp(1:nr_move_south) = move_also_south(1:nr_move_south) 493 485 trnp(1:nr_move_north) = move_also_north(1:nr_move_north) … … 506 498 !-- be moved. 507 499 IF ( particles(n)%particle_mask ) THEN 500 508 501 j = ( particles(n)%y + 0.5_wp * dy ) * ddy 509 502 ! … … 521 514 particles(n)%y = ( ny + 1 ) * dy + particles(n)%y 522 515 particles(n)%origin_y = ( ny + 1 ) * dy + & 523 particles(n)%origin_y516 particles(n)%origin_y 524 517 ELSE 525 trsp_count = trsp_count + 1526 trsp(trsp_count) = particles(n)518 trsp_count = trsp_count + 1 519 trsp(trsp_count) = particles(n) 527 520 trsp(trsp_count)%y = ( ny + 1 ) * dy + & 528 trsp(trsp_count)%y521 trsp(trsp_count)%y 529 522 trsp(trsp_count)%origin_y = trsp(trsp_count)%origin_y & 530 + ( ny + 1 ) * dy523 + ( ny + 1 ) * dy 531 524 particles(n)%particle_mask = .FALSE. 532 525 deleted_particles = deleted_particles + 1 … … 536 529 !++ why is 1 subtracted in next statement??? 537 530 trsp(trsp_count)%origin_y = & 538 trsp(trsp_count)%origin_y - 1531 trsp(trsp_count)%origin_y - 1 539 532 ENDIF 540 533 … … 545 538 !-- Particle absorption 546 539 particles(n)%particle_mask = .FALSE. 547 deleted_particles = deleted_particles + 1540 deleted_particles = deleted_particles + 1 548 541 549 542 ELSEIF ( ibc_par_ns == 2 ) THEN … … 568 561 IF ( j > ny ) THEN 569 562 ! 570 !-- Apply boundary condition along x563 !-- Apply boundary condition along y 571 564 IF ( ibc_par_ns == 0 ) THEN 572 565 ! 573 566 !-- Cyclic condition 574 567 IF ( pdims(2) == 1 ) THEN 575 particles(n)%y = particles(n)%y - ( ny + 1 ) * dy568 particles(n)%y = particles(n)%y - ( ny + 1 ) * dy 576 569 particles(n)%origin_y = & 577 particles(n)%origin_y - ( ny + 1 ) * dy570 particles(n)%origin_y - ( ny + 1 ) * dy 578 571 ELSE 579 trnp_count = trnp_count + 1580 trnp(trnp_count) = particles(n)572 trnp_count = trnp_count + 1 573 trnp(trnp_count) = particles(n) 581 574 trnp(trnp_count)%y = & 582 trnp(trnp_count)%y - ( ny + 1 ) * dy575 trnp(trnp_count)%y - ( ny + 1 ) * dy 583 576 trnp(trnp_count)%origin_y = & 584 trnp(trnp_count)%origin_y - ( ny + 1 ) * dy577 trnp(trnp_count)%origin_y - ( ny + 1 ) * dy 585 578 particles(n)%particle_mask = .FALSE. 586 deleted_particles = deleted_particles + 1579 deleted_particles = deleted_particles + 1 587 580 ENDIF 588 581 … … 628 621 629 622 ALLOCATE(rvnp(MAX(1,trnp_count_recv))) 630 623 631 624 CALL MPI_SENDRECV( trsp(1)%radius, trsp_count, mpi_particle_type, & 632 625 psouth, 1, rvnp(1)%radius, & … … 661 654 ENDIF 662 655 656 DEALLOCATE( move_also_north ) 657 DEALLOCATE( move_also_south ) 658 663 659 #else 664 660 665 ! 666 !-- Apply boundary conditions 667 DO n = 1, number_of_particles 668 669 IF ( particles(n)%x < -0.5_wp * dx ) THEN 670 671 IF ( ibc_par_lr == 0 ) THEN 672 ! 673 !-- Cyclic boundary. Relevant coordinate has to be changed. 674 particles(n)%x = ( nx + 1 ) * dx + particles(n)%x 675 676 ELSEIF ( ibc_par_lr == 1 ) THEN 677 ! 678 !-- Particle absorption 679 particles(n)%particle_mask = .FALSE. 680 deleted_particles = deleted_particles + 1 681 682 ELSEIF ( ibc_par_lr == 2 ) THEN 683 ! 684 !-- Particle reflection 685 particles(n)%x = -dx - particles(n)%x 686 particles(n)%speed_x = -particles(n)%speed_x 687 ENDIF 688 689 ELSEIF ( particles(n)%x >= ( nx + 0.5_wp ) * dx ) THEN 690 691 IF ( ibc_par_lr == 0 ) THEN 692 ! 693 !-- Cyclic boundary. Relevant coordinate has to be changed. 694 particles(n)%x = particles(n)%x - ( nx + 1 ) * dx 695 696 ELSEIF ( ibc_par_lr == 1 ) THEN 697 ! 698 !-- Particle absorption 699 particles(n)%particle_mask = .FALSE. 700 deleted_particles = deleted_particles + 1 701 702 ELSEIF ( ibc_par_lr == 2 ) THEN 703 ! 704 !-- Particle reflection 705 particles(n)%x = ( nx + 1 ) * dx - particles(n)%x 706 particles(n)%speed_x = -particles(n)%speed_x 707 ENDIF 708 709 ENDIF 710 711 IF ( particles(n)%y < -0.5_wp * dy ) THEN 712 713 IF ( ibc_par_ns == 0 ) THEN 714 ! 715 !-- Cyclic boundary. Relevant coordinate has to be changed. 716 particles(n)%y = ( ny + 1 ) * dy + particles(n)%y 717 718 ELSEIF ( ibc_par_ns == 1 ) THEN 719 ! 720 !-- Particle absorption 721 particles(n)%particle_mask = .FALSE. 722 deleted_particles = deleted_particles + 1 723 724 ELSEIF ( ibc_par_ns == 2 ) THEN 725 ! 726 !-- Particle reflection 727 particles(n)%y = -dy - particles(n)%y 728 particles(n)%speed_y = -particles(n)%speed_y 729 ENDIF 730 731 ELSEIF ( particles(n)%y >= ( ny + 0.5_wp ) * dy ) THEN 732 733 IF ( ibc_par_ns == 0 ) THEN 734 ! 735 !-- Cyclic boundary. Relevant coordinate has to be changed. 736 particles(n)%y = particles(n)%y - ( ny + 1 ) * dy 737 738 ELSEIF ( ibc_par_ns == 1 ) THEN 739 ! 740 !-- Particle absorption 741 particles(n)%particle_mask = .FALSE. 742 deleted_particles = deleted_particles + 1 743 744 ELSEIF ( ibc_par_ns == 2 ) THEN 745 ! 746 !-- Particle reflection 747 particles(n)%y = ( ny + 1 ) * dy - particles(n)%y 748 particles(n)%speed_y = -particles(n)%speed_y 749 ENDIF 750 751 ENDIF 661 DO ip = nxl, nxr, nxr-nxl 662 DO jp = nys, nyn 663 DO kp = nzb+1, nzt 664 number_of_particles = prt_count(kp,jp,ip) 665 IF ( number_of_particles <= 0 ) CYCLE 666 particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles) 667 DO n = 1, number_of_particles 668 ! 669 !-- Apply boundary conditions 670 671 IF ( particles(n)%x < -0.5_wp * dx ) THEN 672 673 IF ( ibc_par_lr == 0 ) THEN 674 ! 675 !-- Cyclic boundary. Relevant coordinate has to be changed. 676 particles(n)%x = ( nx + 1 ) * dx + particles(n)%x 677 678 ELSEIF ( ibc_par_lr == 1 ) THEN 679 ! 680 !-- Particle absorption 681 particles(n)%particle_mask = .FALSE. 682 deleted_particles = deleted_particles + 1 683 684 ELSEIF ( ibc_par_lr == 2 ) THEN 685 ! 686 !-- Particle reflection 687 particles(n)%x = -dx - particles(n)%x 688 particles(n)%speed_x = -particles(n)%speed_x 689 ENDIF 690 691 ELSEIF ( particles(n)%x >= ( nx + 0.5_wp ) * dx ) THEN 692 693 IF ( ibc_par_lr == 0 ) THEN 694 ! 695 !-- Cyclic boundary. Relevant coordinate has to be changed. 696 particles(n)%x = particles(n)%x - ( nx + 1 ) * dx 697 698 ELSEIF ( ibc_par_lr == 1 ) THEN 699 ! 700 !-- Particle absorption 701 particles(n)%particle_mask = .FALSE. 702 deleted_particles = deleted_particles + 1 703 704 ELSEIF ( ibc_par_lr == 2 ) THEN 705 ! 706 !-- Particle reflection 707 particles(n)%x = ( nx + 1 ) * dx - particles(n)%x 708 particles(n)%speed_x = -particles(n)%speed_x 709 ENDIF 710 711 ENDIF 712 ENDDO 713 ENDDO 714 ENDDO 752 715 ENDDO 753 716 717 DO ip = nxl, nxr 718 DO jp = nys, nyn, nyn-nys 719 DO kp = nzb+1, nzt 720 number_of_particles = prt_count(kp,jp,ip) 721 IF ( number_of_particles <= 0 ) CYCLE 722 particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles) 723 DO n = 1, number_of_particles 724 725 IF ( particles(n)%y < -0.5_wp * dy ) THEN 726 727 IF ( ibc_par_ns == 0 ) THEN 728 ! 729 !-- Cyclic boundary. Relevant coordinate has to be changed. 730 particles(n)%y = ( ny + 1 ) * dy + particles(n)%y 731 732 ELSEIF ( ibc_par_ns == 1 ) THEN 733 ! 734 !-- Particle absorption 735 particles(n)%particle_mask = .FALSE. 736 deleted_particles = deleted_particles + 1 737 738 ELSEIF ( ibc_par_ns == 2 ) THEN 739 ! 740 !-- Particle reflection 741 particles(n)%y = -dy - particles(n)%y 742 particles(n)%speed_y = -particles(n)%speed_y 743 ENDIF 744 745 ELSEIF ( particles(n)%y >= ( ny + 0.5_wp ) * dy ) THEN 746 747 IF ( ibc_par_ns == 0 ) THEN 748 ! 749 !-- Cyclic boundary. Relevant coordinate has to be changed. 750 particles(n)%y = particles(n)%y - ( ny + 1 ) * dy 751 752 ELSEIF ( ibc_par_ns == 1 ) THEN 753 ! 754 !-- Particle absorption 755 particles(n)%particle_mask = .FALSE. 756 deleted_particles = deleted_particles + 1 757 758 ELSEIF ( ibc_par_ns == 2 ) THEN 759 ! 760 !-- Particle reflection 761 particles(n)%y = ( ny + 1 ) * dy - particles(n)%y 762 particles(n)%speed_y = -particles(n)%speed_y 763 ENDIF 764 765 ENDIF 766 767 ENDDO 768 ENDDO 769 ENDDO 770 ENDDO 754 771 #endif 755 772 … … 782 799 IMPLICIT NONE 783 800 784 INTEGER(iwp) :: ip !< 785 INTEGER(iwp) :: jp !< 786 INTEGER(iwp) :: kp !< 787 INTEGER(iwp) :: n !< 788 INTEGER(iwp) :: pindex !< 801 INTEGER(iwp) :: ip !< grid index (x) of particle 802 INTEGER(iwp) :: jp !< grid index (x) of particle 803 INTEGER(iwp) :: kp !< grid index (x) of particle 804 INTEGER(iwp) :: n !< index variable of particle 805 INTEGER(iwp) :: pindex !< dummy argument for new number of particles per grid box 789 806 790 807 LOGICAL :: pack_done !< 791 808 792 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_array 809 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_array !< new particles in a grid box 810 TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: temp_ns !< temporary particle array for reallocation 793 811 794 812 pack_done = .FALSE. 795 813 796 nr_move_north = 0797 nr_move_south = 0798 799 814 DO n = 1, SIZE(particle_array) 815 816 IF ( .NOT. particle_array(n)%particle_mask ) CYCLE 817 800 818 ip = ( particle_array(n)%x + 0.5_wp * dx ) * ddx 801 819 jp = ( particle_array(n)%y + 0.5_wp * dy ) * ddy 802 kp = particle_array(n)%z / dz + 1 + offset_ocean_nzt820 kp = particle_array(n)%z / dz + 1 + offset_ocean_nzt 803 821 804 822 IF ( ip >= nxl .AND. ip <= nxr .AND. jp >= nys .AND. jp <= nyn & … … 824 842 prt_count(kp,jp,ip) = pindex 825 843 ELSE 826 IF ( jp == nys - 1 ) THEN844 IF ( jp <= nys - 1 ) THEN 827 845 nr_move_south = nr_move_south+1 846 ! 847 !-- Before particle information is swapped to exchange-array, check 848 !-- if enough memory is allocated. If required, reallocate exchange 849 !-- array. 850 IF ( nr_move_south > SIZE(move_also_south) ) THEN 851 ! 852 !-- At first, allocate further temporary array to swap particle 853 !-- information. 854 ALLOCATE( temp_ns(SIZE(move_also_south)+NR_2_direction_move) ) 855 temp_ns(1:nr_move_south-1) = move_also_south(1:nr_move_south-1) 856 DEALLOCATE( move_also_south ) 857 ALLOCATE( move_also_south(SIZE(temp_ns)) ) 858 move_also_south(1:nr_move_south-1) = temp_ns(1:nr_move_south-1) 859 DEALLOCATE( temp_ns ) 860 861 ENDIF 862 828 863 move_also_south(nr_move_south) = particle_array(n) 864 829 865 IF ( jp == -1 ) THEN 830 866 move_also_south(nr_move_south)%y = & … … 833 869 move_also_south(nr_move_south)%origin_y + ( ny + 1 ) * dy 834 870 ENDIF 835 ELSEIF ( jp == nyn+1 ) THEN871 ELSEIF ( jp >= nyn+1 ) THEN 836 872 nr_move_north = nr_move_north+1 873 ! 874 !-- Before particle information is swapped to exchange-array, check 875 !-- if enough memory is allocated. If required, reallocate exchange 876 !-- array. 877 IF ( nr_move_north > SIZE(move_also_north) ) THEN 878 ! 879 !-- At first, allocate further temporary array to swap particle 880 !-- information. 881 ALLOCATE( temp_ns(SIZE(move_also_north)+NR_2_direction_move) ) 882 temp_ns(1:nr_move_north-1) = move_also_south(1:nr_move_north-1) 883 DEALLOCATE( move_also_north ) 884 ALLOCATE( move_also_north(SIZE(temp_ns)) ) 885 move_also_north(1:nr_move_north-1) = temp_ns(1:nr_move_north-1) 886 DEALLOCATE( temp_ns ) 887 888 ENDIF 889 837 890 move_also_north(nr_move_north) = particle_array(n) 838 891 IF ( jp == ny+1 ) THEN … … 867 920 IMPLICIT NONE 868 921 869 INTEGER(iwp) :: i !< 870 INTEGER(iwp) :: ip !< 871 INTEGER(iwp) :: j !< 872 INTEGER(iwp) :: jp !< 873 INTEGER(iwp) :: k !< 874 INTEGER(iwp) :: kp !< 875 INTEGER(iwp) :: n !< 876 INTEGER(iwp) :: np_old_cell !< 877 INTEGER(iwp) :: n_start !< 878 INTEGER(iwp) :: pindex !< 922 INTEGER(iwp) :: i !< grid index (x) of particle position 923 INTEGER(iwp) :: ip !< index variable along x 924 INTEGER(iwp) :: j !< grid index (y) of particle position 925 INTEGER(iwp) :: jp !< index variable along y 926 INTEGER(iwp) :: k !< grid index (z) of particle position 927 INTEGER(iwp) :: kp !< index variable along z 928 INTEGER(iwp) :: n !< index variable for particle array 929 INTEGER(iwp) :: np_old_cell !< number of particles per grid box before moving 930 INTEGER(iwp) :: n_start !< start index 931 INTEGER(iwp) :: pindex !< dummy argument for number of new particle per grid box 879 932 880 933 LOGICAL :: pack_done !< 881 934 882 TYPE(particle_type), DIMENSION(:), POINTER :: particles_old_cell !< 935 TYPE(particle_type), DIMENSION(:), POINTER :: particles_old_cell !< particles before moving 883 936 884 937 CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'start' ) … … 947 1000 IF ( pack_done ) THEN 948 1001 CALL realloc_particles_array(i,j,k) 949 pindex = prt_count(k,j,i)+1950 1002 ELSE 951 1003 CALL lpm_pack_arrays 952 1004 prt_count(k,j,i) = number_of_particles 953 954 pindex = prt_count(k,j,i)+1955 1005 ! 956 1006 !-- If number of particles in the new grid box … … 959 1009 IF ( pindex > SIZE(grid_particles(k,j,i)%particles) ) THEN 960 1010 CALL realloc_particles_array(i,j,k) 961 pindex = prt_count(k,j,i)+1962 1011 ENDIF 963 1012 … … 987 1036 ! Description: 988 1037 ! ------------ 989 !> @todo Missing subroutine description. 1038 !> If the allocated memory for the particle array do not suffice to add arriving 1039 !> particles from neighbour grid cells, this subrouting reallocates the 1040 !> particle array to assure enough memory is available. 990 1041 !------------------------------------------------------------------------------! 991 1042 SUBROUTINE realloc_particles_array (i,j,k,size_in) … … 1012 1063 ENDIF 1013 1064 1014 new_size = MAX( new_size, min_nr_particle )1065 new_size = MAX( new_size, min_nr_particle, old_size + 1 ) 1015 1066 1016 1067 IF ( old_size <= 500 ) THEN -
palm/trunk/SOURCE/lpm_init.f90
r1891 r1929 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Bugfix in determining initial particle height and grid index in case of 22 ! seed_follows_topography. 23 ! Bugfix concerning random positions, ensure that particles do not move more 24 ! than one grid length. 25 ! Bugfix logarithmic interpolation. 26 ! Initial setting of sgs_wf_part. 22 27 ! 23 28 ! Former revisions: … … 31 36 ! 32 37 ! 1873 2016-04-18 14:50:06Z maronga 33 ! Module renamed (removed _mod) 38 ! Module renamed (removed _mod 39 34 40 ! 35 41 ! 1871 2016-04-15 11:46:09Z hoffmann … … 162 168 prt_count, psb, psl, psn, psr, pss, pst, & 163 169 radius, random_start_position, read_particles_from_restartfile,& 164 seed_follows_topography, s ort_count,&170 seed_follows_topography, sgs_wf_part, sort_count, & 165 171 total_number_of_particles, & 166 172 use_sgs_for_particles, & … … 190 196 PUBLIC lpm_init, lpm_create_particle 191 197 192 CONTAINS198 CONTAINS 193 199 194 200 !------------------------------------------------------------------------------! … … 290 296 291 297 ! 292 !-- Allocate arrays required for calculating particle SGS velocities 298 !-- Allocate arrays required for calculating particle SGS velocities. 299 !-- Initialize prefactor required for stoachastic Weil equation. 293 300 IF ( use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN 294 301 ALLOCATE( de_dx(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 295 302 de_dy(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 296 303 de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 304 305 sgs_wf_part = 1.0_wp / 3.0_wp 297 306 ENDIF 298 307 … … 334 343 ! 335 344 !-- Precalculate LOG(z/z0) 336 height_p = 0.0_wp345 height_p = z0_av_global 337 346 DO k = 1, number_of_sublayers 338 347 … … 341 350 342 351 ENDDO 343 344 352 345 353 ENDIF … … 444 452 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & 445 453 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0, 0, 0, & 446 0, .FALSE., -1 )454 0, .FALSE., -1 ) 447 455 448 456 particle_groups = particle_groups_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp ) … … 477 485 478 486 CALL lpm_create_particle (PHASE_INIT) 479 480 487 ! 481 488 !-- User modification of initial particles … … 499 506 number_of_particles = prt_count(nzb+1,nys,nxl) 500 507 particles => grid_particles(nzb+1,nys,nxl)%particles(1:number_of_particles) 501 502 508 ! 503 509 !-- Formats … … 520 526 521 527 USE particle_attributes, & 522 ONLY: monodisperse_aerosols528 ONLY: deleted_particles, monodisperse_aerosols 523 529 524 530 IMPLICIT NONE 525 531 526 INTEGER(iwp) :: alloc_size !< 527 INTEGER(iwp) :: i !< 528 INTEGER(iwp) :: ip !< 529 INTEGER(iwp) :: j !< 530 INTEGER(iwp) :: jp !< 531 INTEGER(iwp) :: kp !< 532 INTEGER(iwp) :: loop_stride !< 533 INTEGER(iwp) :: n !< 534 INTEGER(iwp) :: new_size !< 535 536 INTEGER(iwp), INTENT(IN) :: phase !< 537 538 INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_count !< 539 INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_start !< 540 541 LOGICAL :: first_stride !< 542 543 REAL(wp) :: pos_x !< 544 REAL(wp) :: pos_y !< 545 REAL(wp) :: pos_z !< 546 547 TYPE(particle_type),TARGET :: tmp_particle !< 532 INTEGER(iwp) :: alloc_size !< relative increase of allocated memory for particles 533 INTEGER(iwp) :: i !< loop variable ( particle groups ) 534 INTEGER(iwp) :: ip !< index variable along x 535 INTEGER(iwp) :: j !< loop variable ( particles per point ) 536 INTEGER(iwp) :: jp !< index variable along y 537 INTEGER(iwp) :: kp !< index variable along z 538 INTEGER(iwp) :: loop_stride !< loop variable for initialization 539 INTEGER(iwp) :: n !< loop variable ( number of particles ) 540 INTEGER(iwp) :: new_size !< new size of allocated memory for particles 541 542 INTEGER(iwp), INTENT(IN) :: phase !< mode of inititialization 543 544 INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_count !< start address of new particle 545 INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_start !< start address of new particle 546 547 LOGICAL :: first_stride !< flag for initialization 548 549 REAL(wp) :: pos_x !< increment for particle position in x 550 REAL(wp) :: pos_y !< increment for particle position in y 551 REAL(wp) :: pos_z !< increment for particle position in z 552 REAL(wp) :: rand_contr !< dummy argument for random position 553 554 TYPE(particle_type),TARGET :: tmp_particle !< temporary particle used for initialization 548 555 549 556 ! … … 618 625 tmp_particle%tail_id = 0 !unused 619 626 627 620 628 ! 621 629 !-- Determine the grid indices of the particle position … … 628 636 !-- Particle height is given relative to topography 629 637 kp = kp + nzb_w_inner(jp,ip) 630 tmp_particle%z = tmp_particle%z + zw(kp) 638 tmp_particle%z = tmp_particle%z + & 639 zw(nzb_w_inner(jp,ip)) 631 640 IF ( kp > nzt ) THEN 632 641 pos_x = pos_x + pdx(i) 633 642 CYCLE xloop 634 643 ENDIF 644 ELSEIF ( .NOT. seed_follows_topography .AND. & 645 tmp_particle%z <= zw(nzb_w_inner(jp,ip)) ) THEN 646 pos_x = pos_x + pdx(i) 647 CYCLE xloop 635 648 ENDIF 636 649 … … 644 657 ENDIF 645 658 grid_particles(kp,jp,ip)%particles(local_count(kp,jp,ip)) = tmp_particle 659 646 660 ENDIF 647 661 ENDDO … … 691 705 ENDIF 692 706 ENDIF 707 693 708 ENDDO 694 709 ENDDO 695 710 ENDDO 696 711 ENDIF 712 697 713 ENDDO 698 714 … … 707 723 708 724 ! 709 !-- Add random fluctuation to particle positions 725 !-- Add random fluctuation to particle positions. 710 726 IF ( random_start_position ) THEN 711 727 DO ip = nxl, nxr … … 715 731 IF ( number_of_particles <= 0 ) CYCLE 716 732 particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles) 717 718 DO n = local_start(kp,jp,ip), number_of_particles !Move only new particles 733 ! 734 !-- Move only new particles. Moreover, limit random fluctuation 735 !-- in order to prevent that particles move more than one grid box, 736 !-- which would lead to problems concerning particle exchange 737 !-- between processors in case pdx/pdy are larger than dx/dy, 738 !-- respectively. 739 DO n = local_start(kp,jp,ip), number_of_particles 719 740 IF ( psl(particles(n)%group) /= psr(particles(n)%group) ) THEN 741 rand_contr = ( random_function( iran_part ) - 0.5_wp ) * & 742 pdx(particles(n)%group) 720 743 particles(n)%x = particles(n)%x + & 721 ( random_function( iran_part ) - 0.5_wp ) * & 722 pdx(particles(n)%group) 744 MERGE( rand_contr, SIGN( dx, rand_contr ), & 745 ABS( rand_contr ) < dx & 746 ) 723 747 ENDIF 724 748 IF ( pss(particles(n)%group) /= psn(particles(n)%group) ) THEN 749 rand_contr = ( random_function( iran_part ) - 0.5_wp ) * & 750 pdy(particles(n)%group) 725 751 particles(n)%y = particles(n)%y + & 726 ( random_function( iran_part ) - 0.5_wp ) * & 727 pdy(particles(n)%group) 752 MERGE( rand_contr, SIGN( dy, rand_contr ), & 753 ABS( rand_contr ) < dy & 754 ) 728 755 ENDIF 729 756 IF ( psb(particles(n)%group) /= pst(particles(n)%group) ) THEN 757 rand_contr = ( random_function( iran_part ) - 0.5_wp ) * & 758 pdz(particles(n)%group) 730 759 particles(n)%z = particles(n)%z + & 731 ( random_function( iran_part ) - 0.5_wp ) * & 732 pdz(particles(n)%group) 760 MERGE( rand_contr, SIGN( dz, rand_contr ), & 761 ABS( rand_contr ) < dz & 762 ) 733 763 ENDIF 734 764 ENDDO 735 765 ! 736 !-- Identify particles located outside the model domain 766 !-- Identify particles located outside the model domain and reflect 767 !-- or absorb them if necessary. 737 768 CALL lpm_boundary_conds( 'bottom/top' ) 769 ! 770 !-- Furthermore, remove particles located in topography. Note, as 771 !-- the particle speed is still zero at this point, wall 772 !-- reflection boundary conditions will not work in this case. 773 particles => & 774 grid_particles(kp,jp,ip)%particles(1:number_of_particles) 775 DO n = local_start(kp,jp,ip), number_of_particles 776 i = ( particles(n)%x + 0.5_wp * dx ) * ddx 777 j = ( particles(n)%y + 0.5_wp * dy ) * ddy 778 IF ( particles(n)%z <= zw(nzb_w_inner(j,i)) ) THEN 779 particles(n)%particle_mask = .FALSE. 780 deleted_particles = deleted_particles + 1 781 ENDIF 782 ENDDO 738 783 ENDDO 739 784 ENDDO -
palm/trunk/SOURCE/lpm_init_sgs_tke.f90
r1823 r1929 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! sgs_wfu_par, sgs_wfv_par and sgs_wfw_par are replaced by sgs_wf_par 22 22 ! 23 23 ! Former revisions: … … 72 72 73 73 USE particle_attributes, & 74 ONLY: sgs_wf u_part, sgs_wfv_part, sgs_wfw_part74 ONLY: sgs_wf_part 75 75 76 76 USE pegrid … … 95 95 k > nzb_s_inner(j,i+1) ) & 96 96 THEN 97 de_dx(k,j,i) = 2.0_wp * sgs_wf u_part *&97 de_dx(k,j,i) = 2.0_wp * sgs_wf_part * & 98 98 ( e(k,j,i+1) - e(k,j,i) ) * ddx 99 99 ELSEIF ( k > nzb_s_inner(j,i-1) .AND. k > nzb_s_inner(j,i) & 100 100 .AND. k <= nzb_s_inner(j,i+1) ) & 101 101 THEN 102 de_dx(k,j,i) = 2.0_wp * sgs_wf u_part *&102 de_dx(k,j,i) = 2.0_wp * sgs_wf_part * & 103 103 ( e(k,j,i) - e(k,j,i-1) ) * ddx 104 104 ELSEIF ( k < nzb_s_inner(j,i) .AND. k < nzb_s_inner(j,i+1) ) & … … 109 109 de_dx(k,j,i) = 0.0_wp 110 110 ELSE 111 de_dx(k,j,i) = sgs_wf u_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx111 de_dx(k,j,i) = sgs_wf_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx 112 112 ENDIF 113 113 … … 115 115 k > nzb_s_inner(j+1,i) ) & 116 116 THEN 117 de_dy(k,j,i) = 2.0_wp * sgs_wf v_part *&117 de_dy(k,j,i) = 2.0_wp * sgs_wf_part * & 118 118 ( e(k,j+1,i) - e(k,j,i) ) * ddy 119 119 ELSEIF ( k > nzb_s_inner(j-1,i) .AND. k > nzb_s_inner(j,i) & 120 120 .AND. k <= nzb_s_inner(j+1,i) ) & 121 121 THEN 122 de_dy(k,j,i) = 2.0_wp * sgs_wf v_part *&122 de_dy(k,j,i) = 2.0_wp * sgs_wf_part * & 123 123 ( e(k,j,i) - e(k,j-1,i) ) * ddy 124 124 ELSEIF ( k < nzb_s_inner(j,i) .AND. k < nzb_s_inner(j+1,i) ) & … … 129 129 de_dy(k,j,i) = 0.0_wp 130 130 ELSE 131 de_dy(k,j,i) = sgs_wf v_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy131 de_dy(k,j,i) = sgs_wf_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy 132 132 ENDIF 133 133 … … 142 142 143 143 DO k = nzb_s_inner(j,i)+2, nzt-1 144 de_dz(k,j,i) = 2.0_wp * sgs_wf w_part *&144 de_dz(k,j,i) = 2.0_wp * sgs_wf_part * & 145 145 ( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1)-zu(k-1) ) 146 146 ENDDO … … 148 148 k = nzb_s_inner(j,i) 149 149 de_dz(nzb:k,j,i) = 0.0_wp 150 de_dz(k+1,j,i) = 2.0_wp * sgs_wf w_part *&150 de_dz(k+1,j,i) = 2.0_wp * sgs_wf_part * & 151 151 ( e(k+2,j,i) - e(k+1,j,i) ) / ( zu(k+2) - zu(k+1) ) 152 152 de_dz(nzt,j,i) = 0.0_wp -
palm/trunk/SOURCE/mod_particle_attributes.f90
r1872 r1929 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! -sgs_wfu_par, sgs_wfv_par, sgs_wfw_par 22 ! + sgs_wf_par 22 23 ! 23 24 ! Former revisions: … … 109 110 initial_weighting_factor = 1.0_wp, & 110 111 particle_advection_start = 0.0_wp, & 111 sgs_wfu_part = 0.3333333_wp, sgs_wfv_part = 0.3333333_wp, & 112 sgs_wfw_part = 0.3333333_wp, time_prel = 0.0_wp, & 113 time_sort_particles = 0.0_wp, & 112 sgs_wf_part, time_prel = 0.0_wp, time_sort_particles = 0.0_wp,& 114 113 time_write_particle_data = 0.0_wp, z0_av_global 115 114 -
palm/trunk/SOURCE/pres.f90
r1919 r1929 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! Bugfix: weight_substep for initial call, replace by local variable 22 22 ! 23 23 ! Former revisions: … … 162 162 REAL(wp) :: threadsum !< 163 163 REAL(wp) :: weight_pres_l !< 164 REAL(wp) :: weight_substep_l !< 164 165 165 166 REAL(wp), DIMENSION(1:3) :: volume_flow_l !< … … 178 179 ! 179 180 !-- If pres is called before initial time step 180 weight_pres_l = 1.0_wp 181 d_weight_pres = 1.0_wp 181 weight_pres_l = 1.0_wp 182 d_weight_pres = 1.0_wp 183 weight_substep_l = 1.0_wp 182 184 ELSE 183 weight_pres_l = weight_pres(intermediate_timestep_count) 184 d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count) 185 weight_pres_l = weight_pres(intermediate_timestep_count) 186 d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count) 187 weight_substep_l = weight_substep(intermediate_timestep_count) 185 188 ENDIF 186 189 … … 591 594 !$OMP PARALLEL PRIVATE (i,j,k) 592 595 !$OMP DO 593 !$acc kernels present( p, tend, weight_substep )596 !$acc kernels present( p, tend, weight_substep_l ) 594 597 !$acc loop independent 595 598 DO i = nxl-1, nxr+1 … … 599 602 DO k = nzb, nzt+1 600 603 p(k,j,i) = tend(k,j,i) * & 601 weight_substep (intermediate_timestep_count)604 weight_substep_l 602 605 ENDDO 603 606 ENDDO … … 609 612 !$OMP PARALLEL PRIVATE (i,j,k) 610 613 !$OMP DO 611 !$acc kernels present( p, tend, weight_substep )614 !$acc kernels present( p, tend, weight_substep_l ) 612 615 !$acc loop independent 613 616 DO i = nxl-1, nxr+1 … … 617 620 DO k = nzb, nzt+1 618 621 p(k,j,i) = p(k,j,i) + tend(k,j,i) * & 619 weight_substep (intermediate_timestep_count)622 weight_substep_l 620 623 ENDDO 621 624 ENDDO -
palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
r1921 r1929 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! Bugfix: avoid segmentation fault in case one grid point is horizontally 22 ! completely surrounded by topography 22 23 ! 23 24 ! Former revisions: … … 726 727 !-- consequence would result in very large shear stresses and very 727 728 !-- small momentum fluxes (both are generally unrealistic). 728 IF ( ( z_mo / ol(j,i) ) < zeta_min ) ol(j,i) = z_mo / zeta_min 729 IF ( ( z_mo / ol(j,i) ) > zeta_max ) ol(j,i) = z_mo / zeta_max 729 IF ( ( z_mo / ( ol(j,i) + 1E-30_wp ) ) < zeta_min ) & 730 ol(j,i) = z_mo / zeta_min 731 IF ( ( z_mo / ( ol(j,i) + 1E-30_wp ) ) > zeta_max ) & 732 ol(j,i) = z_mo / zeta_max 730 733 ENDDO 731 734 ENDDO -
palm/trunk/SOURCE/wind_turbine_model_mod.f90
r1917 r1929 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Bugfix: added preprocessor directives for parallel and serial mode 22 22 ! 23 23 ! Former revisions: … … 2058 2058 2059 2059 ! 2060 !-- Transfer of information to the other cpus 2060 !-- Transfer of information to the other cpus 2061 #if defined( __parallel ) 2061 2062 CALL MPI_ALLREDUCE( u_inflow_l, u_inflow, nturbines, MPI_REAL, & 2062 2063 MPI_SUM, comm2d, ierr ) … … 2065 2066 CALL MPI_ALLREDUCE( phi_yaw_l, phi_yaw, nturbines, MPI_REAL, & 2066 2067 MPI_SUM, comm2d, ierr ) 2067 2068 #else 2069 u_inflow = u_inflow_l 2070 wdir = wdir_l 2071 phi_yaw = phi_yaw_l 2072 #endif 2068 2073 DO inot = 1, nturbines 2069 2074 ! … … 2080 2085 ! CALL MPI_ALLREDUCE( omega_gen, omega_gen_old, nturbines, & 2081 2086 ! MPI_REAL,MPI_SUM, comm2d, ierr ) 2087 #if defined( __parallel ) 2082 2088 CALL MPI_ALLREDUCE( torque_gen, torque_gen_old, nturbines, & 2083 2089 MPI_REAL, MPI_SUM, comm2d, ierr ) … … 2086 2092 CALL MPI_ALLREDUCE( omega_gen_f, omega_gen_f_old, nturbines, & 2087 2093 MPI_REAL, MPI_SUM, comm2d, ierr ) 2094 #else 2095 torque_gen_old = torque_gen 2096 omega_rot = omega_rot_l 2097 omega_gen_f_old = omega_gen_f 2098 #endif 2088 2099 2089 2100 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
