Changeset 1359 for palm/trunk/SOURCE/lpm_init.f90
- Timestamp:
- Apr 11, 2014 5:15:14 PM (10 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/lpm_init.f90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_init.f90
r1329 r1359 1 SUBROUTINE lpm_init1 MODULE lpm_init_mod 2 2 3 3 !--------------------------------------------------------------------------------! … … 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! New particle structure integrated. 23 ! Kind definition added to all floating point numbers. 24 ! lpm_init changed form a subroutine to a module. 23 25 ! 24 26 ! Former revisions: … … 85 87 86 88 USE control_parameters, & 87 ONLY: cloud_droplets, current_timestep_number, initializing_actions,&88 message_string, netcdf_data_format, ocean,&89 prandtl_layer, simulated_time89 ONLY: cloud_droplets, current_timestep_number, dz, & 90 initializing_actions, message_string, netcdf_data_format, & 91 ocean, prandtl_layer, simulated_time 90 92 91 93 USE dvrp_variables, & … … 93 95 94 96 USE grid_variables, & 95 ONLY: d x, dy97 ONLY: ddx, dx, ddy, dy 96 98 97 99 USE indices, & … … 104 106 105 107 USE particle_attributes, & 106 ONLY: bc_par_b, bc_par_lr, bc_par_ns, bc_par_t, collision_kernel, & 107 density_ratio, dvrp_psize, initial_weighting_factor, ibc_par_b,& 108 ibc_par_lr, ibc_par_ns, ibc_par_t, initial_particles, & 109 iran_part, log_z_z0, max_number_of_particle_groups, & 110 maximum_number_of_particles, maximum_number_of_tailpoints, & 111 minimum_tailpoint_distance, maximum_number_of_tails, & 112 mpi_particle_type, new_tail_id, number_of_initial_particles, & 108 ONLY: alloc_factor, bc_par_b, bc_par_lr, bc_par_ns, bc_par_t, & 109 block_offset, block_offset_def, collision_kernel, & 110 density_ratio, dvrp_psize, grid_particles, & 111 initial_weighting_factor, ibc_par_b, ibc_par_lr, ibc_par_ns, & 112 ibc_par_t, iran_part, log_z_z0, & 113 max_number_of_particle_groups, maximum_number_of_particles, & 114 maximum_number_of_tailpoints, maximum_number_of_tails, & 115 minimum_tailpoint_distance, min_nr_particle, & 116 mpi_particle_type, new_tail_id, & 113 117 number_of_initial_tails, number_of_particles, & 114 118 number_of_particle_groups, number_of_sublayers, & 115 number_of_tails, offset_ocean_nzt, offset_ocean_nzt_m1, part_1,&116 part _2, particles, particle_advection_start, particle_groups,&117 particle_groups_type, particle _mask, particles_per_point,&119 number_of_tails, offset_ocean_nzt, offset_ocean_nzt_m1, & 120 particles, particle_advection_start, particle_groups, & 121 particle_groups_type, particles_per_point, & 118 122 particle_tail_coordinates, particle_type, pdx, pdy, pdz, & 119 prt_count, p rt_start_index, psb, psl, psn, psr, pss, pst,&123 prt_count, psb, psl, psn, psr, pss, pst, & 120 124 radius, random_start_position, read_particles_from_restartfile,& 121 skip_particles_for_tail, sort_count, tail_mask,&122 t otal_number_of_particles, total_number_of_tails,&125 skip_particles_for_tail, sort_count, & 126 tail_mask, total_number_of_particles, total_number_of_tails, & 123 127 use_particle_tails, use_sgs_for_particles, & 124 write_particle_statistics, uniform_particles, z0_av_global 128 write_particle_statistics, uniform_particles, zero_particle, & 129 z0_av_global 125 130 126 131 USE pegrid … … 129 134 ONLY: random_function 130 135 131 132 136 IMPLICIT NONE 133 137 138 PRIVATE 139 140 INTEGER(iwp), PARAMETER :: PHASE_INIT = 1 !: 141 INTEGER(iwp), PARAMETER, PUBLIC :: PHASE_RELEASE = 2 !: 142 143 INTERFACE lpm_init 144 MODULE PROCEDURE lpm_init 145 END INTERFACE lpm_init 146 147 INTERFACE lpm_create_particle 148 MODULE PROCEDURE lpm_create_particle 149 END INTERFACE lpm_create_particle 150 151 PUBLIC lpm_init, lpm_create_particle 152 153 CONTAINS 154 155 SUBROUTINE lpm_init 156 157 USE lpm_collision_kernels_mod, & 158 ONLY: init_kernels 159 160 IMPLICIT NONE 161 134 162 INTEGER(iwp) :: i !: 163 INTEGER(iwp) :: ip !: 135 164 INTEGER(iwp) :: j !: 165 INTEGER(iwp) :: jp !: 136 166 INTEGER(iwp) :: k !: 167 INTEGER(iwp) :: kp !: 137 168 INTEGER(iwp) :: n !: 138 169 INTEGER(iwp) :: nn !: … … 145 176 LOGICAL :: uniform_particles_l !: 146 177 147 REAL(wp) :: height_int !: 148 REAL(wp) :: height_p !: 149 REAL(wp) :: pos_x !: 150 REAL(wp) :: pos_y !: 151 REAL(wp) :: pos_z !: 152 REAL(wp) :: z_p !: 153 REAL(wp) :: z0_av_local = 0.0 !: 154 155 156 178 REAL(wp) :: height_int !: 179 REAL(wp) :: height_p !: 180 REAL(wp) :: z_p !: 181 REAL(wp) :: z0_av_local !: 157 182 158 183 #if defined( __parallel ) … … 160 185 !-- Define MPI derived datatype for FORTRAN datatype particle_type (see module 161 186 !-- particle_attributes). Integer length is 4 byte, Real is 8 byte 162 blocklengths(1) = 19; blocklengths(2) = 4; blocklengths(3) = 1 163 displacements(1) = 0; displacements(2) = 152; displacements(3) = 168 187 #if defined( __twocachelines ) 188 blocklengths(1) = 7; blocklengths(2) = 18; blocklengths(3) = 1 189 displacements(1) = 0; displacements(2) = 64; displacements(3) = 128 190 191 types(1) = MPI_REAL ! 64 bit words 192 types(2) = MPI_INTEGER ! 32 Bit words 193 types(3) = MPI_UB 194 #else 195 blocklengths(1) = 19; blocklengths(2) = 6; blocklengths(3) = 1 196 displacements(1) = 0; displacements(2) = 152; displacements(3) = 176 164 197 165 198 types(1) = MPI_REAL 166 199 types(2) = MPI_INTEGER 167 200 types(3) = MPI_UB 201 #endif 168 202 CALL MPI_TYPE_STRUCT( 3, blocklengths, displacements, types, & 169 203 mpi_particle_type, ierr ) … … 179 213 ENDIF 180 214 181 215 ! 216 !-- Define block offsets for dividing a gridcell in 8 sub cells 217 218 block_offset(0) = block_offset_def (-1,-1,-1) 219 block_offset(1) = block_offset_def (-1,-1, 0) 220 block_offset(2) = block_offset_def (-1, 0,-1) 221 block_offset(3) = block_offset_def (-1, 0, 0) 222 block_offset(4) = block_offset_def ( 0,-1,-1) 223 block_offset(5) = block_offset_def ( 0,-1, 0) 224 block_offset(6) = block_offset_def ( 0, 0,-1) 225 block_offset(7) = block_offset_def ( 0, 0, 0) 182 226 ! 183 227 !-- Check the number of particle groups. … … 193 237 ! 194 238 !-- Set default start positions, if necessary 195 IF ( psl(1) == 9999999.9 ) psl(1) = -0.5* dx196 IF ( psr(1) == 9999999.9 ) psr(1) = ( nx + 0.5) * dx197 IF ( pss(1) == 9999999.9 ) pss(1) = -0.5* dy198 IF ( psn(1) == 9999999.9 ) psn(1) = ( ny + 0.5) * dy199 IF ( psb(1) == 9999999.9 ) psb(1) = zu(nz/2)200 IF ( pst(1) == 9999999.9 ) pst(1) = psb(1)201 202 IF ( pdx(1) == 9999999.9 .OR. pdx(1) == 0.0) pdx(1) = dx203 IF ( pdy(1) == 9999999.9 .OR. pdy(1) == 0.0) pdy(1) = dy204 IF ( pdz(1) == 9999999.9 .OR. pdz(1) == 0.0) pdz(1) = zu(2) - zu(1)239 IF ( psl(1) == 9999999.9_wp ) psl(1) = -0.5_wp * dx 240 IF ( psr(1) == 9999999.9_wp ) psr(1) = ( nx + 0.5_wp ) * dx 241 IF ( pss(1) == 9999999.9_wp ) pss(1) = -0.5_wp * dy 242 IF ( psn(1) == 9999999.9_wp ) psn(1) = ( ny + 0.5_wp ) * dy 243 IF ( psb(1) == 9999999.9_wp ) psb(1) = zu(nz/2) 244 IF ( pst(1) == 9999999.9_wp ) pst(1) = psb(1) 245 246 IF ( pdx(1) == 9999999.9_wp .OR. pdx(1) == 0.0_wp ) pdx(1) = dx 247 IF ( pdy(1) == 9999999.9_wp .OR. pdy(1) == 0.0_wp ) pdy(1) = dy 248 IF ( pdz(1) == 9999999.9_wp .OR. pdz(1) == 0.0_wp ) pdz(1) = zu(2) - zu(1) 205 249 206 250 DO j = 2, number_of_particle_groups 207 IF ( psl(j) == 9999999.9 ) psl(j) = psl(j-1)208 IF ( psr(j) == 9999999.9 ) psr(j) = psr(j-1)209 IF ( pss(j) == 9999999.9 ) pss(j) = pss(j-1)210 IF ( psn(j) == 9999999.9 ) psn(j) = psn(j-1)211 IF ( psb(j) == 9999999.9 ) psb(j) = psb(j-1)212 IF ( pst(j) == 9999999.9 ) pst(j) = pst(j-1)213 IF ( pdx(j) == 9999999.9 .OR. pdx(j) == 0.0) pdx(j) = pdx(j-1)214 IF ( pdy(j) == 9999999.9 .OR. pdy(j) == 0.0) pdy(j) = pdy(j-1)215 IF ( pdz(j) == 9999999.9 .OR. pdz(j) == 0.0) pdz(j) = pdz(j-1)251 IF ( psl(j) == 9999999.9_wp ) psl(j) = psl(j-1) 252 IF ( psr(j) == 9999999.9_wp ) psr(j) = psr(j-1) 253 IF ( pss(j) == 9999999.9_wp ) pss(j) = pss(j-1) 254 IF ( psn(j) == 9999999.9_wp ) psn(j) = psn(j-1) 255 IF ( psb(j) == 9999999.9_wp ) psb(j) = psb(j-1) 256 IF ( pst(j) == 9999999.9_wp ) pst(j) = pst(j-1) 257 IF ( pdx(j) == 9999999.9_wp .OR. pdx(j) == 0.0_wp ) pdx(j) = pdx(j-1) 258 IF ( pdy(j) == 9999999.9_wp .OR. pdy(j) == 0.0_wp ) pdy(j) = pdy(j-1) 259 IF ( pdz(j) == 9999999.9_wp .OR. pdz(j) == 0.0_wp ) pdz(j) = pdz(j-1) 216 260 ENDDO 217 261 … … 243 287 !-- negligible. 244 288 z0_av_local = SUM( z0(nys:nyn,nxl:nxr) ) 245 z0_av_global = 0.0 289 z0_av_global = 0.0_wp 246 290 247 291 #if defined( __parallel ) … … 255 299 ! 256 300 !-- Horizontal wind speed is zero below and at z0 257 log_z_z0(0) = 0.0 301 log_z_z0(0) = 0.0_wp 258 302 ! 259 303 !-- Calculate vertical depth of the sublayers … … 261 305 ! 262 306 !-- Precalculate LOG(z/z0) 263 height_p = 0.0 307 height_p = 0.0_wp 264 308 DO k = 1, number_of_sublayers 265 309 … … 269 313 ENDDO 270 314 271 272 ENDIF273 274 !275 !-- Initialize collision kernels276 IF ( collision_kernel /= 'none' ) CALL init_kernels277 278 !279 !-- For the first model run of a possible job chain initialize the280 !-- particles, otherwise read the particle data from restart file.281 IF ( TRIM( initializing_actions ) == 'read_restart_data' &282 .AND. read_particles_from_restartfile ) THEN283 284 CALL lpm_read_restart_file285 286 ELSE287 288 !289 !-- Allocate particle arrays and set attributes of the initial set of290 !-- particles, which can be also periodically released at later times.291 !-- Also allocate array for particle tail coordinates, if needed.292 ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &293 prt_start_index(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &294 particle_mask(maximum_number_of_particles), &295 part_1(maximum_number_of_particles), &296 part_2(maximum_number_of_particles) )297 298 particles => part_1299 300 sort_count = 0301 302 !303 !-- Initialize all particles with dummy values (otherwise errors may304 !-- occur within restart runs). The reason for this is still not clear305 !-- and may be presumably caused by errors in the respective user-interface.306 particles = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &307 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &308 0.0, 0, 0, 0, 0 )309 particle_groups = particle_groups_type( 0.0, 0.0, 0.0, 0.0 )310 311 !312 !-- Set the default particle size used for dvrp plots313 IF ( dvrp_psize == 9999999.9 ) dvrp_psize = 0.2 * dx314 315 !316 !-- Set values for the density ratio and radius for all particle317 !-- groups, if necessary318 IF ( density_ratio(1) == 9999999.9 ) density_ratio(1) = 0.0319 IF ( radius(1) == 9999999.9 ) radius(1) = 0.0320 DO i = 2, number_of_particle_groups321 IF ( density_ratio(i) == 9999999.9 ) THEN322 density_ratio(i) = density_ratio(i-1)323 ENDIF324 IF ( radius(i) == 9999999.9 ) radius(i) = radius(i-1)325 ENDDO326 327 DO i = 1, number_of_particle_groups328 IF ( density_ratio(i) /= 0.0 .AND. radius(i) == 0 ) THEN329 WRITE( message_string, * ) 'particle group #', i, 'has a', &330 'density ratio /= 0 but radius = 0'331 CALL message( 'lpm_init', 'PA0215', 1, 2, 0, 6, 0 )332 ENDIF333 particle_groups(i)%density_ratio = density_ratio(i)334 particle_groups(i)%radius = radius(i)335 ENDDO336 337 !338 !-- Calculate particle positions and store particle attributes, if339 !-- particle is situated on this PE340 n = 0341 342 DO i = 1, number_of_particle_groups343 344 pos_z = psb(i)345 346 DO WHILE ( pos_z <= pst(i) )347 348 pos_y = pss(i)349 350 DO WHILE ( pos_y <= psn(i) )351 352 IF ( pos_y >= ( nys - 0.5 ) * dy .AND. &353 pos_y < ( nyn + 0.5 ) * dy ) THEN354 355 pos_x = psl(i)356 357 DO WHILE ( pos_x <= psr(i) )358 359 IF ( pos_x >= ( nxl - 0.5 ) * dx .AND. &360 pos_x < ( nxr + 0.5 ) * dx ) THEN361 362 DO j = 1, particles_per_point363 364 n = n + 1365 IF ( n > maximum_number_of_particles ) THEN366 WRITE( message_string, * ) 'number of initial', &367 'particles (', n, ') exceeds', &368 '&maximum_number_of_particles (', &369 maximum_number_of_particles, ') on PE ', &370 myid371 CALL message( 'lpm_init', 'PA0216', 2, 2, -1, 6,&372 1 )373 ENDIF374 particles(n)%x = pos_x375 particles(n)%y = pos_y376 particles(n)%z = pos_z377 particles(n)%age = 0.0378 particles(n)%age_m = 0.0379 particles(n)%dt_sum = 0.0380 particles(n)%dvrp_psize = dvrp_psize381 particles(n)%e_m = 0.0382 IF ( curvature_solution_effects ) THEN383 !384 !-- Initial values (internal timesteps, derivative)385 !-- for Rosenbrock method386 particles(n)%rvar1 = 1.0E-12387 particles(n)%rvar2 = 1.0E-3388 particles(n)%rvar3 = -9999999.9389 ELSE390 !391 !-- Initial values for SGS velocities392 particles(n)%rvar1 = 0.0393 particles(n)%rvar2 = 0.0394 particles(n)%rvar3 = 0.0395 ENDIF396 particles(n)%speed_x = 0.0397 particles(n)%speed_y = 0.0398 particles(n)%speed_z = 0.0399 particles(n)%origin_x = pos_x400 particles(n)%origin_y = pos_y401 particles(n)%origin_z = pos_z402 particles(n)%radius = particle_groups(i)%radius403 particles(n)%weight_factor =initial_weighting_factor404 particles(n)%class = 1405 particles(n)%group = i406 particles(n)%tailpoints = 0407 IF ( use_particle_tails .AND. &408 MOD( n, skip_particles_for_tail ) == 0 ) THEN409 number_of_tails = number_of_tails + 1410 !411 !-- This is a temporary provisional setting (see412 !-- further below!)413 particles(n)%tail_id = number_of_tails414 ELSE415 particles(n)%tail_id = 0416 ENDIF417 418 ENDDO419 420 ENDIF421 422 pos_x = pos_x + pdx(i)423 424 ENDDO425 426 ENDIF427 428 pos_y = pos_y + pdy(i)429 430 ENDDO431 432 pos_z = pos_z + pdz(i)433 434 ENDDO435 436 ENDDO437 438 number_of_initial_particles = n439 number_of_particles = n440 441 !442 !-- Calculate the number of particles and tails of the total domain443 #if defined( __parallel )444 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )445 CALL MPI_ALLREDUCE( number_of_particles, total_number_of_particles, 1, &446 MPI_INTEGER, MPI_SUM, comm2d, ierr )447 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )448 CALL MPI_ALLREDUCE( number_of_tails, total_number_of_tails, 1, &449 MPI_INTEGER, MPI_SUM, comm2d, ierr )450 #else451 total_number_of_particles = number_of_particles452 total_number_of_tails = number_of_tails453 #endif454 455 !456 !-- Set a seed value for the random number generator to be exclusively457 !-- used for the particle code. The generated random numbers should be458 !-- different on the different PEs.459 iran_part = iran_part + myid460 461 !462 !-- User modification of initial particles463 CALL user_lpm_init464 465 !466 !-- Store the initial set of particles for release at later times467 IF ( number_of_initial_particles /= 0 ) THEN468 ALLOCATE( initial_particles(1:number_of_initial_particles) )469 initial_particles(1:number_of_initial_particles) = &470 particles(1:number_of_initial_particles)471 ENDIF472 473 !474 !-- Add random fluctuation to particle positions475 IF ( random_start_position ) THEN476 477 DO n = 1, number_of_initial_particles478 IF ( psl(particles(n)%group) /= psr(particles(n)%group) ) THEN479 particles(n)%x = particles(n)%x + &480 ( random_function( iran_part ) - 0.5 ) * &481 pdx(particles(n)%group)482 IF ( particles(n)%x <= ( nxl - 0.5 ) * dx ) THEN483 particles(n)%x = ( nxl - 0.4999999999 ) * dx484 ELSEIF ( particles(n)%x >= ( nxr + 0.5 ) * dx ) THEN485 particles(n)%x = ( nxr + 0.4999999999 ) * dx486 ENDIF487 ENDIF488 IF ( pss(particles(n)%group) /= psn(particles(n)%group) ) THEN489 particles(n)%y = particles(n)%y + &490 ( random_function( iran_part ) - 0.5 ) * &491 pdy(particles(n)%group)492 IF ( particles(n)%y <= ( nys - 0.5 ) * dy ) THEN493 particles(n)%y = ( nys - 0.4999999999 ) * dy494 ELSEIF ( particles(n)%y >= ( nyn + 0.5 ) * dy ) THEN495 particles(n)%y = ( nyn + 0.4999999999 ) * dy496 ENDIF497 ENDIF498 IF ( psb(particles(n)%group) /= pst(particles(n)%group) ) THEN499 particles(n)%z = particles(n)%z + &500 ( random_function( iran_part ) - 0.5 ) * &501 pdz(particles(n)%group)502 ENDIF503 ENDDO504 ENDIF505 506 !507 !-- Sort particles in the sequence the gridboxes are stored in the memory.508 !-- Only required if cloud droplets are used.509 IF ( cloud_droplets ) CALL lpm_sort_arrays510 511 !512 !-- Open file for statistical informations about particle conditions513 IF ( write_particle_statistics ) THEN514 CALL check_open( 80 )515 WRITE ( 80, 8000 ) current_timestep_number, simulated_time, &516 number_of_initial_particles, &517 maximum_number_of_particles518 CALL close_file( 80 )519 ENDIF520 521 !522 !-- Check if particles are really uniform in color and radius (dvrp_size)523 !-- (uniform_particles is preset TRUE)524 IF ( uniform_particles ) THEN525 IF ( number_of_initial_particles == 0 ) THEN526 uniform_particles_l = .TRUE.527 ELSE528 n = number_of_initial_particles529 IF ( MINVAL( particles(1:n)%dvrp_psize ) == &530 MAXVAL( particles(1:n)%dvrp_psize ) .AND. &531 MINVAL( particles(1:n)%class ) == &532 MAXVAL( particles(1:n)%class ) ) THEN533 uniform_particles_l = .TRUE.534 ELSE535 uniform_particles_l = .FALSE.536 ENDIF537 ENDIF538 539 #if defined( __parallel )540 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )541 CALL MPI_ALLREDUCE( uniform_particles_l, uniform_particles, 1, &542 MPI_LOGICAL, MPI_LAND, comm2d, ierr )543 #else544 uniform_particles = uniform_particles_l545 #endif546 547 ENDIF548 549 !550 !-- Particles will probably become none-uniform, if their size and color551 !-- will be determined by flow variables552 IF ( particle_color /= 'none' .OR. particle_dvrpsize /= 'none' ) THEN553 uniform_particles = .FALSE.554 ENDIF555 556 !557 !-- Set the beginning of the particle tails and their age558 IF ( use_particle_tails ) THEN559 !560 !-- Choose the maximum number of tails with respect to the maximum number561 !-- of particles and skip_particles_for_tail562 maximum_number_of_tails = maximum_number_of_particles / &563 skip_particles_for_tail564 565 !566 !-- Create a minimum number of tails in case that there is no tail567 !-- initially (otherwise, index errors will occur when adressing the568 !-- arrays below)569 IF ( maximum_number_of_tails == 0 ) maximum_number_of_tails = 10570 571 ALLOCATE( particle_tail_coordinates(maximum_number_of_tailpoints,5, &572 maximum_number_of_tails), &573 new_tail_id(maximum_number_of_tails), &574 tail_mask(maximum_number_of_tails) )575 576 particle_tail_coordinates = 0.0577 minimum_tailpoint_distance = minimum_tailpoint_distance**2578 number_of_initial_tails = number_of_tails579 580 nn = 0581 DO n = 1, number_of_particles582 !583 !-- Only for those particles marked above with a provisional tail_id584 !-- tails will be created. Particles now get their final tail_id.585 IF ( particles(n)%tail_id /= 0 ) THEN586 587 nn = nn + 1588 particles(n)%tail_id = nn589 590 particle_tail_coordinates(1,1,nn) = particles(n)%x591 particle_tail_coordinates(1,2,nn) = particles(n)%y592 particle_tail_coordinates(1,3,nn) = particles(n)%z593 particle_tail_coordinates(1,4,nn) = particles(n)%class594 particles(n)%tailpoints = 1595 IF ( minimum_tailpoint_distance /= 0.0 ) THEN596 particle_tail_coordinates(2,1,nn) = particles(n)%x597 particle_tail_coordinates(2,2,nn) = particles(n)%y598 particle_tail_coordinates(2,3,nn) = particles(n)%z599 particle_tail_coordinates(2,4,nn) = particles(n)%class600 particle_tail_coordinates(1:2,5,nn) = 0.0601 particles(n)%tailpoints = 2602 ENDIF603 604 ENDIF605 ENDDO606 ENDIF607 608 !609 !-- Plot initial positions of particles (only if particle advection is610 !-- switched on from the beginning of the simulation (t=0))611 IF ( particle_advection_start == 0.0 ) CALL data_output_dvrp612 315 613 316 ENDIF … … 677 380 678 381 END SELECT 382 383 ! 384 !-- Initialize collision kernels 385 IF ( collision_kernel /= 'none' ) CALL init_kernels 386 387 ! 388 !-- For the first model run of a possible job chain initialize the 389 !-- particles, otherwise read the particle data from restart file. 390 IF ( TRIM( initializing_actions ) == 'read_restart_data' & 391 .AND. read_particles_from_restartfile ) THEN 392 393 CALL lpm_read_restart_file 394 395 ELSE 396 397 ! 398 !-- Allocate particle arrays and set attributes of the initial set of 399 !-- particles, which can be also periodically released at later times. 400 !-- Also allocate array for particle tail coordinates, if needed. 401 ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 402 grid_particles(nzb+1:nzt,nys:nyn,nxl:nxr) ) 403 404 maximum_number_of_particles = 0 405 number_of_particles = 0 406 407 sort_count = 0 408 prt_count = 0 409 410 ! 411 !-- Initialize all particles with dummy values (otherwise errors may 412 !-- occur within restart runs). The reason for this is still not clear 413 !-- and may be presumably caused by errors in the respective user-interface. 414 #if defined( __twocachelines ) 415 zero_particle = particle_type( 0.0_wp, 0.0_sp, 0.0_sp, 0.0_sp, 0.0_sp, & 416 0.0_sp, 0.0_sp, 0.0_wp, 0.0_wp, 0.0_wp, & 417 0, .FALSE., 0.0_wp, 0.0_wp, 0.0_wp, & 418 0.0_sp, 0.0_sp, 0.0_sp, 0.0_sp, 0.0_sp, & 419 0.0_sp, 0, 0, 0, -1) 420 #else 421 zero_particle = particle_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & 422 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & 423 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, & 424 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0, 0, 0, & 425 0, .FALSE., -1) 426 #endif 427 particle_groups = particle_groups_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp ) 428 429 ! 430 !-- Set the default particle size used for dvrp plots 431 IF ( dvrp_psize == 9999999.9_wp ) dvrp_psize = 0.2_wp * dx 432 433 ! 434 !-- Set values for the density ratio and radius for all particle 435 !-- groups, if necessary 436 IF ( density_ratio(1) == 9999999.9_wp ) density_ratio(1) = 0.0_wp 437 IF ( radius(1) == 9999999.9_wp ) radius(1) = 0.0_wp 438 DO i = 2, number_of_particle_groups 439 IF ( density_ratio(i) == 9999999.9_wp ) THEN 440 density_ratio(i) = density_ratio(i-1) 441 ENDIF 442 IF ( radius(i) == 9999999.9_wp ) radius(i) = radius(i-1) 443 ENDDO 444 445 DO i = 1, number_of_particle_groups 446 IF ( density_ratio(i) /= 0.0_wp .AND. radius(i) == 0 ) THEN 447 WRITE( message_string, * ) 'particle group #', i, 'has a', & 448 'density ratio /= 0 but radius = 0' 449 CALL message( 'lpm_init', 'PA0215', 1, 2, 0, 6, 0 ) 450 ENDIF 451 particle_groups(i)%density_ratio = density_ratio(i) 452 particle_groups(i)%radius = radius(i) 453 ENDDO 454 455 CALL lpm_create_particle (PHASE_INIT) 456 ! 457 !-- Set a seed value for the random number generator to be exclusively 458 !-- used for the particle code. The generated random numbers should be 459 !-- different on the different PEs. 460 iran_part = iran_part + myid 461 462 ! 463 !-- User modification of initial particles 464 CALL user_lpm_init 465 466 ! 467 !-- Open file for statistical informations about particle conditions 468 IF ( write_particle_statistics ) THEN 469 CALL check_open( 80 ) 470 WRITE ( 80, 8000 ) current_timestep_number, simulated_time, & 471 number_of_particles, & 472 maximum_number_of_particles 473 CALL close_file( 80 ) 474 ENDIF 475 476 ! 477 !-- Check if particles are really uniform in color and radius (dvrp_size) 478 !-- (uniform_particles is preset TRUE) 479 IF ( uniform_particles ) THEN 480 DO ip = nxl, nxr 481 DO jp = nys, nyn 482 DO kp = nzb+1, nzt 483 484 n = prt_count(kp,jp,ip) 485 IF ( MINVAL( grid_particles(kp,jp,ip)%particles(1:n)%dvrp_psize ) == & 486 MAXVAL( grid_particles(kp,jp,ip)%particles(1:n)%dvrp_psize ) .AND. & 487 MINVAL( grid_particles(kp,jp,ip)%particles(1:n)%class ) == & 488 MAXVAL( grid_particles(kp,jp,ip)%particles(1:n)%class ) ) THEN 489 uniform_particles_l = .TRUE. 490 ELSE 491 uniform_particles_l = .FALSE. 492 ENDIF 493 494 ENDDO 495 ENDDO 496 ENDDO 497 498 #if defined( __parallel ) 499 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 500 CALL MPI_ALLREDUCE( uniform_particles_l, uniform_particles, 1, & 501 MPI_LOGICAL, MPI_LAND, comm2d, ierr ) 502 #else 503 uniform_particles = uniform_particles_l 504 #endif 505 506 ENDIF 507 508 ! 509 !-- Particles will probably become none-uniform, if their size and color 510 !-- will be determined by flow variables 511 IF ( particle_color /= 'none' .OR. particle_dvrpsize /= 'none' ) THEN 512 uniform_particles = .FALSE. 513 ENDIF 514 515 ! !kk Not implemented aft individual particle array fort every gridcell 516 ! ! 517 ! !-- Set the beginning of the particle tails and their age 518 ! IF ( use_particle_tails ) THEN 519 ! ! 520 ! !-- Choose the maximum number of tails with respect to the maximum number 521 ! !-- of particles and skip_particles_for_tail 522 ! maximum_number_of_tails = maximum_number_of_particles / & 523 ! skip_particles_for_tail 524 ! 525 ! ! 526 ! !-- Create a minimum number of tails in case that there is no tail 527 ! !-- initially (otherwise, index errors will occur when adressing the 528 ! !-- arrays below) 529 ! IF ( maximum_number_of_tails == 0 ) maximum_number_of_tails = 10 530 ! 531 ! ALLOCATE( particle_tail_coordinates(maximum_number_of_tailpoints,5, & 532 ! maximum_number_of_tails), & 533 ! new_tail_id(maximum_number_of_tails), & 534 ! tail_mask(maximum_number_of_tails) ) 535 ! 536 ! particle_tail_coordinates = 0.0_wp 537 ! minimum_tailpoint_distance = minimum_tailpoint_distance**2 538 ! number_of_initial_tails = number_of_tails 539 ! 540 ! nn = 0 541 ! DO n = 1, number_of_particles 542 ! ! 543 ! !-- Only for those particles marked above with a provisional tail_id 544 ! !-- tails will be created. Particles now get their final tail_id. 545 ! IF ( particles(n)%tail_id /= 0 ) THEN 546 ! 547 ! nn = nn + 1 548 ! particles(n)%tail_id = nn 549 ! 550 ! particle_tail_coordinates(1,1,nn) = particles(n)%x 551 ! particle_tail_coordinates(1,2,nn) = particles(n)%y 552 ! particle_tail_coordinates(1,3,nn) = particles(n)%z 553 ! particle_tail_coordinates(1,4,nn) = particles(n)%class 554 ! particles(n)%tailpoints = 1 555 ! IF ( minimum_tailpoint_distance /= 0.0_wp ) THEN 556 ! particle_tail_coordinates(2,1,nn) = particles(n)%x 557 ! particle_tail_coordinates(2,2,nn) = particles(n)%y 558 ! particle_tail_coordinates(2,3,nn) = particles(n)%z 559 ! particle_tail_coordinates(2,4,nn) = particles(n)%class 560 ! particle_tail_coordinates(1:2,5,nn) = 0.0_wp 561 ! particles(n)%tailpoints = 2 562 ! ENDIF 563 ! 564 ! ENDIF 565 ! ENDDO 566 ! ENDIF 567 ! 568 ! ! 569 ! !-- Plot initial positions of particles (only if particle advection is 570 ! !-- switched on from the beginning of the simulation (t=0)) 571 ! IF ( particle_advection_start == 0.0_wp ) CALL data_output_dvrp 572 573 ENDIF 574 575 ! 576 !-- To avoid programm abort, assign particles array to the local version of 577 !-- first grid cell 578 number_of_particles = prt_count(nzb+1,nys,nxl) 579 particles => grid_particles(nzb+1,nys,nxl)%particles(1:number_of_particles) 580 679 581 ! 680 582 !-- Formats 681 8000 FORMAT (I6,1X,F7.2,4X,I 6,71X,I6)583 8000 FORMAT (I6,1X,F7.2,4X,I10,71X,I10) 682 584 683 585 END SUBROUTINE lpm_init 586 587 SUBROUTINE lpm_create_particle (phase) 588 589 USE lpm_exchange_horiz_mod, & 590 ONLY: lpm_exchange_horiz, lpm_move_particle, realloc_particles_array 591 592 USE lpm_pack_arrays_mod, & 593 ONLY: lpm_pack_all_arrays 594 595 IMPLICIT NONE 596 597 INTEGER(iwp) :: alloc_size !: 598 INTEGER(iwp) :: i !: 599 INTEGER(iwp) :: ip !: 600 INTEGER(iwp) :: j !: 601 INTEGER(iwp) :: jp !: 602 INTEGER(iwp) :: kp !: 603 INTEGER(iwp) :: loop_stride !: 604 INTEGER(iwp) :: n !: 605 INTEGER(iwp) :: new_size !: 606 INTEGER(iwp) :: nn !: 607 608 INTEGER(iwp), INTENT(IN) :: phase !: 609 610 INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_count !: 611 INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_start !: 612 613 LOGICAL :: first_stride !: 614 615 REAL(wp) :: pos_x !: 616 REAL(wp) :: pos_y !: 617 REAL(wp) :: pos_z !: 618 619 TYPE(particle_type),TARGET :: tmp_particle !: 620 621 ! 622 !-- Calculate particle positions and store particle attributes, if 623 !-- particle is situated on this PE 624 DO loop_stride = 1, 2 625 first_stride = (loop_stride == 1) 626 IF ( first_stride ) THEN 627 local_count = 0 ! count number of particles 628 ELSE 629 local_count = prt_count ! Start address of new particles 630 ENDIF 631 632 n = 0 633 DO i = 1, number_of_particle_groups 634 635 pos_z = psb(i) 636 637 DO WHILE ( pos_z <= pst(i) ) 638 639 pos_y = pss(i) 640 641 DO WHILE ( pos_y <= psn(i) ) 642 643 IF ( pos_y >= ( nys - 0.5_wp ) * dy .AND. & 644 pos_y < ( nyn + 0.5_wp ) * dy ) THEN 645 646 pos_x = psl(i) 647 648 DO WHILE ( pos_x <= psr(i) ) 649 650 IF ( pos_x >= ( nxl - 0.5_wp ) * dx .AND. & 651 pos_x < ( nxr + 0.5_wp ) * dx ) THEN 652 653 DO j = 1, particles_per_point 654 655 n = n + 1 656 #if defined( __twocachelines ) 657 tmp_particle%x = pos_x 658 tmp_particle%y = pos_y 659 tmp_particle%z = pos_z 660 tmp_particle%age = 0.0_sp 661 tmp_particle%age_m = 0.0_sp 662 tmp_particle%dt_sum = 0.0_wp 663 tmp_particle%dvrp_psize = dvrp_psize 664 tmp_particle%e_m = 0.0_sp 665 IF ( curvature_solution_effects ) THEN 666 ! 667 !-- Initial values (internal timesteps, derivative) 668 !-- for Rosenbrock method 669 tmp_particle%rvar1 = 1.0E-12_wp 670 tmp_particle%rvar2 = 1.0E-3_wp 671 tmp_particle%rvar3 = -9999999.9_wp 672 ELSE 673 ! 674 !-- Initial values for SGS velocities 675 tmp_particle%rvar1 = 0.0_wp 676 tmp_particle%rvar2 = 0.0_wp 677 tmp_particle%rvar3 = 0.0_wp 678 ENDIF 679 tmp_particle%speed_x = 0.0_sp 680 tmp_particle%speed_y = 0.0_sp 681 tmp_particle%speed_z = 0.0_sp 682 tmp_particle%origin_x = pos_x 683 tmp_particle%origin_y = pos_y 684 tmp_particle%origin_z = pos_z 685 tmp_particle%radius = particle_groups(i)%radius 686 tmp_particle%weight_factor = initial_weighting_factor 687 tmp_particle%class = 1 688 tmp_particle%group = i 689 tmp_particle%tailpoints = 0 690 tmp_particle%particle_mask = .TRUE. 691 #else 692 tmp_particle%x = pos_x 693 tmp_particle%y = pos_y 694 tmp_particle%z = pos_z 695 tmp_particle%age = 0.0_wp 696 tmp_particle%age_m = 0.0_wp 697 tmp_particle%dt_sum = 0.0_wp 698 tmp_particle%dvrp_psize = dvrp_psize 699 tmp_particle%e_m = 0.0_wp 700 IF ( curvature_solution_effects ) THEN 701 ! 702 !-- Initial values (internal timesteps, derivative) 703 !-- for Rosenbrock method 704 tmp_particle%rvar1 = 1.0E-12_wp 705 tmp_particle%rvar2 = 1.0E-3_wp 706 tmp_particle%rvar3 = -9999999.9_wp 707 ELSE 708 ! 709 !-- Initial values for SGS velocities 710 tmp_particle%rvar1 = 0.0_wp 711 tmp_particle%rvar2 = 0.0_wp 712 tmp_particle%rvar3 = 0.0_wp 713 ENDIF 714 tmp_particle%speed_x = 0.0_wp 715 tmp_particle%speed_y = 0.0_wp 716 tmp_particle%speed_z = 0.0_wp 717 tmp_particle%origin_x = pos_x 718 tmp_particle%origin_y = pos_y 719 tmp_particle%origin_z = pos_z 720 tmp_particle%radius = particle_groups(i)%radius 721 tmp_particle%weight_factor = initial_weighting_factor 722 tmp_particle%class = 1 723 tmp_particle%group = i 724 tmp_particle%tailpoints = 0 725 tmp_particle%particle_mask = .TRUE. 726 #endif 727 IF ( use_particle_tails .AND. & 728 MOD( n, skip_particles_for_tail ) == 0 ) THEN 729 number_of_tails = number_of_tails + 1 730 ! 731 !-- This is a temporary provisional setting (see 732 !-- further below!) 733 tmp_particle%tail_id = number_of_tails 734 ELSE 735 tmp_particle%tail_id = 0 736 ENDIF 737 ip = ( tmp_particle%x + 0.5_wp * dx ) * ddx 738 jp = ( tmp_particle%y + 0.5_wp * dy ) * ddy 739 kp = tmp_particle%z / dz + 1 + offset_ocean_nzt_m1 740 741 local_count(kp,jp,ip) = local_count(kp,jp,ip) + 1 742 IF ( .NOT. first_stride ) THEN 743 IF ( ip < nxl .OR. jp < nys .OR. kp < nzb+1 ) THEN 744 write(6,*) 'xl ',ip,jp,kp,nxl,nys,nzb+1 745 ENDIF 746 IF ( ip > nxr .OR. jp > nyn .OR. kp > nzt ) THEN 747 write(6,*) 'xu ',ip,jp,kp,nxr,nyn,nzt 748 ENDIF 749 grid_particles(kp,jp,ip)%particles(local_count(kp,jp,ip)) = tmp_particle 750 ENDIF 751 ENDDO 752 753 ENDIF 754 755 pos_x = pos_x + pdx(i) 756 757 ENDDO 758 759 ENDIF 760 761 pos_y = pos_y + pdy(i) 762 763 ENDDO 764 765 pos_z = pos_z + pdz(i) 766 767 ENDDO 768 769 ENDDO 770 771 IF ( first_stride ) THEN 772 DO ip = nxl, nxr 773 DO jp = nys, nyn 774 DO kp = nzb+1, nzt 775 IF ( phase == PHASE_INIT ) THEN 776 IF ( local_count(kp,jp,ip) > 0 ) THEN 777 alloc_size = MAX( INT( local_count(kp,jp,ip) * & 778 ( 1.0_wp + alloc_factor / 100.0_wp ) ), & 779 min_nr_particle ) 780 ELSE 781 alloc_size = min_nr_particle 782 ENDIF 783 ALLOCATE(grid_particles(kp,jp,ip)%particles(1:alloc_size)) 784 DO n = 1, alloc_size 785 grid_particles(kp,jp,ip)%particles(n) = zero_particle 786 ENDDO 787 ELSEIF ( phase == PHASE_RELEASE ) THEN 788 IF ( local_count(kp,jp,ip) > 0 ) THEN 789 new_size = local_count(kp,jp,ip) + prt_count(kp,jp,ip) 790 alloc_size = MAX( INT( new_size * ( 1.0_wp + & 791 alloc_factor / 100.0_wp ) ), min_nr_particle ) 792 IF( alloc_size > SIZE( grid_particles(kp,jp,ip)%particles) ) THEN 793 CALL realloc_particles_array(ip,jp,kp,alloc_size) 794 ENDIF 795 ENDIF 796 ENDIF 797 ENDDO 798 ENDDO 799 ENDDO 800 ENDIF 801 ENDDO 802 803 local_start = prt_count+1 804 prt_count = local_count 805 ! 806 !-- Add random fluctuation to particle positions 807 IF ( random_start_position ) THEN 808 DO ip = nxl, nxr 809 DO jp = nys, nyn 810 DO kp = nzb+1, nzt 811 number_of_particles = prt_count(kp,jp,ip) 812 IF ( number_of_particles <= 0 ) CYCLE 813 particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles) 814 815 DO n = local_start(kp,jp,ip), number_of_particles !Move only new particles 816 IF ( psl(particles(n)%group) /= psr(particles(n)%group) ) THEN 817 particles(n)%x = particles(n)%x + & 818 ( random_function( iran_part ) - 0.5_wp ) * & 819 pdx(particles(n)%group) 820 ENDIF 821 IF ( pss(particles(n)%group) /= psn(particles(n)%group) ) THEN 822 particles(n)%y = particles(n)%y + & 823 ( random_function( iran_part ) - 0.5_wp ) * & 824 pdy(particles(n)%group) 825 ENDIF 826 IF ( psb(particles(n)%group) /= pst(particles(n)%group) ) THEN 827 particles(n)%z = particles(n)%z + & 828 ( random_function( iran_part ) - 0.5_wp ) * & 829 pdz(particles(n)%group) 830 ENDIF 831 ENDDO 832 ! 833 !-- Identify particles located outside the model domain 834 CALL lpm_boundary_conds( 'bottom/top' ) 835 ENDDO 836 ENDDO 837 ENDDO 838 ! 839 !-- Exchange particles between grid cells and processors 840 CALL lpm_move_particle 841 CALL lpm_exchange_horiz 842 843 ENDIF 844 ! 845 !-- In case of random_start_position, delete particles identified by 846 !-- lpm_exchange_horiz and lpm_boundary_conds. Then sort particles into blocks, 847 !-- which is needed for a fast interpolation of the LES fields on the particle 848 !-- position. 849 CALL lpm_pack_all_arrays 850 851 ! 852 !-- Determine maximum number of particles (i.e., all possible particles that 853 !-- have been allocated) and the current number of particles 854 DO ip = nxl, nxr 855 DO jp = nys, nyn 856 DO kp = nzb+1, nzt 857 maximum_number_of_particles = maximum_number_of_particles & 858 + SIZE(grid_particles(kp,jp,ip)%particles) 859 number_of_particles = number_of_particles & 860 + prt_count(kp,jp,ip) 861 ENDDO 862 ENDDO 863 ENDDO 864 ! 865 !-- Calculate the number of particles and tails of the total domain 866 #if defined( __parallel ) 867 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 868 CALL MPI_ALLREDUCE( number_of_particles, total_number_of_particles, 1, & 869 MPI_INTEGER, MPI_SUM, comm2d, ierr ) 870 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 871 CALL MPI_ALLREDUCE( number_of_tails, total_number_of_tails, 1, & 872 MPI_INTEGER, MPI_SUM, comm2d, ierr ) 873 #else 874 total_number_of_particles = number_of_particles 875 total_number_of_tails = number_of_tails 876 #endif 877 878 RETURN 879 880 END SUBROUTINE lpm_create_particle 881 882 END MODULE lpm_init_mod
Note: See TracChangeset
for help on using the changeset viewer.
