source: palm/tags/release-3.2b/SOURCE/init_particles.f90 @ 366

Last change on this file since 366 was 83, checked in by raasch, 17 years ago

New:
---

Changed:


PALM can be generally installed on any kind of Linux-, IBM-AIX-, or NEC-SX-system by adding appropriate settings to the configuration file.

Scripts are also running under the public domain ksh.

All system relevant compile and link options as well as the host identifier (local_host) are specified in the configuration file.

Filetransfer by ftp removed (options -f removed from mrun and mbuild).

Call of (system-)FLUSH routine moved to new routine local_flush.

return_addres and return_username are read from ENVPAR-NAMELIST-file instead of using local_getenv.

Preprocessor strings for different linux clusters changed to "lc", some preprocessor directives renamed (new: intel_openmp_bug), preprocessor directives for old systems removed

advec_particles, check_open, cpu_log, cpu_statistics, data_output_dvrp, flow_statistics, header, init_dvrp, init_particles, init_1d_model, init_dvrp, init_pegrid, local_getenv, local_system, local_tremain, local_tremain_ini, modules, palm, parin, run_control

new:
local_flush

mbuild, mrun

Errors:


  • Property svn:keywords set to Id
File size: 20.4 KB
Line 
1 SUBROUTINE init_particles
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: init_particles.f90 83 2007-04-19 16:27:07Z raasch $
11!
12! 82 2007-04-16 15:40:52Z raasch
13! Preprocessor directives for old systems removed
14!
15! 70 2007-03-18 23:46:30Z raasch
16! displacements for mpi_particle_type changed, age_m initialized,
17! particles-package is now part of the default code
18!
19! 16 2007-02-15 13:16:47Z raasch
20! Bugfix: MPI_REAL in MPI_ALLREDUCE replaced by MPI_INTEGER
21!
22! r4 | raasch | 2007-02-13 12:33:16 +0100 (Tue, 13 Feb 2007)
23! RCS Log replace by Id keyword, revision history cleaned up
24!
25! Revision 1.24  2007/02/11 13:00:17  raasch
26! Bugfix: allocation of tail_mask and new_tail_id in case of restart-runs
27! Bugfix: __ was missing in a cpp-directive
28!
29! Revision 1.1  1999/11/25 16:22:38  raasch
30! Initial revision
31!
32!
33! Description:
34! ------------
35! This routine initializes a set of particles and their attributes (position,
36! radius, ..). Advection of these particles is carried out by advec_particles,
37! plotting is done in data_output_dvrp.
38!------------------------------------------------------------------------------!
39
40    USE arrays_3d
41    USE control_parameters
42    USE grid_variables
43    USE indices
44    USE particle_attributes
45    USE pegrid
46    USE random_function_mod
47
48
49    IMPLICIT NONE
50
51    CHARACTER (LEN=10) ::  particle_binary_version, version_on_file
52
53    INTEGER ::  i, j, n, nn
54#if defined( __parallel )
55    INTEGER, DIMENSION(3) ::  blocklengths, displacements, types
56#endif
57    LOGICAL ::  uniform_particles_l
58    REAL    ::  factor, pos_x, pos_y, pos_z, value
59
60
61#if defined( __parallel )
62!
63!-- Define MPI derived datatype for FORTRAN datatype particle_type (see module
64!-- particle_attributes). Integer length is 4 byte, Real is 8 byte
65    blocklengths(1)  = 19;  blocklengths(2)  =   4;  blocklengths(3)  =   1
66    displacements(1) =  0;  displacements(2) = 152;  displacements(3) = 168
67
68    types(1) = MPI_REAL
69    types(2) = MPI_INTEGER
70    types(3) = MPI_UB
71    CALL MPI_TYPE_STRUCT( 3, blocklengths, displacements, types, &
72                          mpi_particle_type, ierr )
73    CALL MPI_TYPE_COMMIT( mpi_particle_type, ierr )
74#endif
75
76!
77!-- Check the number of particle groups.
78    IF ( number_of_particle_groups > max_number_of_particle_groups )  THEN
79       PRINT*, '+++ WARNING: init_particles: ', &
80                    'max_number_of_particle_groups =', &
81               max_number_of_particle_groups
82       PRINT*, '+++          number_of_particle_groups reset to ', &
83               max_number_of_particle_groups
84       number_of_particle_groups = max_number_of_particle_groups
85    ENDIF
86
87!
88!-- Set default start positions, if necessary
89    IF ( psl(1) == 9999999.9 )  psl(1) = -0.5 * dx
90    IF ( psr(1) == 9999999.9 )  psr(1) = ( nx + 0.5 ) * dx
91    IF ( pss(1) == 9999999.9 )  pss(1) = -0.5 * dy
92    IF ( psn(1) == 9999999.9 )  psn(1) = ( ny + 0.5 ) * dy
93    IF ( psb(1) == 9999999.9 )  psb(1) = zu(nz/2)
94    IF ( pst(1) == 9999999.9 )  pst(1) = psb(1)
95
96    IF ( pdx(1) == 9999999.9  .OR.  pdx(1) == 0.0 )  pdx(1) = dx
97    IF ( pdy(1) == 9999999.9  .OR.  pdy(1) == 0.0 )  pdy(1) = dy
98    IF ( pdz(1) == 9999999.9  .OR.  pdz(1) == 0.0 )  pdz(1) = zu(2) - zu(1)
99
100    DO  j = 2, number_of_particle_groups
101       IF ( psl(j) == 9999999.9 )  psl(j) = psl(j-1)
102       IF ( psr(j) == 9999999.9 )  psr(j) = psr(j-1)
103       IF ( pss(j) == 9999999.9 )  pss(j) = pss(j-1)
104       IF ( psn(j) == 9999999.9 )  psn(j) = psn(j-1)
105       IF ( psb(j) == 9999999.9 )  psb(j) = psb(j-1)
106       IF ( pst(j) == 9999999.9 )  pst(j) = pst(j-1)
107       IF ( pdx(j) == 9999999.9  .OR.  pdx(j) == 0.0 )  pdx(j) = pdx(j-1)
108       IF ( pdy(j) == 9999999.9  .OR.  pdy(j) == 0.0 )  pdy(j) = pdy(j-1)
109       IF ( pdz(j) == 9999999.9  .OR.  pdz(j) == 0.0 )  pdz(j) = pdz(j-1)
110    ENDDO
111
112!
113!-- For the first model run of a possible job chain initialize the
114!-- particles, otherwise read the particle data from file.
115    IF ( TRIM( initializing_actions ) == 'read_restart_data'  &
116         .AND.  read_particles_from_restartfile )  THEN
117
118!
119!--    Read particle data from previous model run.
120!--    First open the input unit.
121       IF ( myid_char == '' )  THEN
122          OPEN ( 90, FILE='PARTICLE_RESTART_DATA_IN'//myid_char, &
123                     FORM='UNFORMATTED' )
124       ELSE
125          OPEN ( 90, FILE='PARTICLE_RESTART_DATA_IN/'//myid_char, &
126                     FORM='UNFORMATTED' )
127       ENDIF
128
129!
130!--    First compare the version numbers
131       READ ( 90 )  version_on_file
132       particle_binary_version = '3.0'
133       IF ( TRIM( version_on_file ) /= TRIM( particle_binary_version ) )  THEN
134          IF ( myid == 0 )  THEN
135             PRINT*, '+++ init_particles: version mismatch concerning data ', &
136                     'from prior run'
137             PRINT*, '        version on file    = "', TRIM( version_on_file ),&
138                     '"'
139             PRINT*, '        version in program = "', &
140                     TRIM( particle_binary_version ), '"'
141          ENDIF
142          CALL local_stop
143       ENDIF
144
145!
146!--    Read some particle parameters and the size of the particle arrays,
147!--    allocate them and read their contents.
148       READ ( 90 )  bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,                  &
149                    maximum_number_of_particles, maximum_number_of_tailpoints, &
150                    maximum_number_of_tails, number_of_initial_particles,      &
151                    number_of_particles, number_of_particle_groups,            &
152                    number_of_tails, particle_groups, time_prel,               &
153                    time_write_particle_data, uniform_particles
154
155       IF ( number_of_initial_particles /= 0 )  THEN
156          ALLOCATE( initial_particles(1:number_of_initial_particles) )
157          READ ( 90 )  initial_particles
158       ENDIF
159
160       ALLOCATE( prt_count(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),       &
161                 prt_start_index(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
162                 particle_mask(maximum_number_of_particles),         &
163                 particles(maximum_number_of_particles) )
164
165       READ ( 90 )  prt_count, prt_start_index
166       READ ( 90 )  particles
167
168       IF ( use_particle_tails )  THEN
169          ALLOCATE( particle_tail_coordinates(maximum_number_of_tailpoints,5, &
170                    maximum_number_of_tails),                                 &
171                    new_tail_id(maximum_number_of_tails),                     &
172                    tail_mask(maximum_number_of_tails) )
173          READ ( 90 )  particle_tail_coordinates
174       ENDIF
175
176       CLOSE ( 90 )
177
178    ELSE
179
180!
181!--    Allocate particle arrays and set attributes of the initial set of
182!--    particles, which can be also periodically released at later times.
183!--    Also allocate array for particle tail coordinates, if needed.
184       ALLOCATE( prt_count(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),       &
185                 prt_start_index(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
186                 particle_mask(maximum_number_of_particles),         &
187                 particles(maximum_number_of_particles) )
188
189!
190!--    Initialize all particles with dummy values (otherwise errors may
191!--    occur within restart runs). The reason for this is still not clear
192!--    and may be presumably caused by errors in the respective user-interface.
193       particles = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
194                                  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
195                                  0.0, 0, 0, 0, 0 )
196       particle_groups = particle_groups_type( 0.0, 0.0, 0.0, 0.0 )
197
198!
199!--    Set the default particle size used for dvrp plots
200       IF ( dvrp_psize == 9999999.9 )  dvrp_psize = 0.2 * dx
201
202!
203!--    Set values for the density ratio and radius for all particle
204!--    groups, if necessary
205       IF ( density_ratio(1) == 9999999.9 )  density_ratio(1) = 0.0
206       IF ( radius(1)        == 9999999.9 )  radius(1) = 0.0
207       DO  i = 2, number_of_particle_groups
208          IF ( density_ratio(i) == 9999999.9 )  THEN
209             density_ratio(i) = density_ratio(i-1)
210          ENDIF
211          IF ( radius(i) == 9999999.9 )  radius(i) = radius(i-1)
212       ENDDO
213
214       DO  i = 1, number_of_particle_groups
215          IF ( density_ratio(i) /= 0.0  .AND.  radius(i) == 0 )  THEN
216             IF ( myid == 0 )  THEN
217                PRINT*, '+++ init_particles: particle group #', i, 'has a', &
218                        'density ratio /= 0 but radius = 0'
219             ENDIF
220             CALL local_stop
221          ENDIF
222          particle_groups(i)%density_ratio = density_ratio(i)
223          particle_groups(i)%radius        = radius(i)
224       ENDDO
225
226!
227!--    Calculate particle positions and store particle attributes, if
228!--    particle is situated on this PE
229       n = 0
230
231       DO  i = 1, number_of_particle_groups
232
233          pos_z = psb(i)
234
235          DO WHILE ( pos_z <= pst(i) )
236
237             pos_y = pss(i)
238
239             DO WHILE ( pos_y <= psn(i) )
240
241                IF ( pos_y >= ( nys - 0.5 ) * dy  .AND.  &
242                     pos_y <  ( nyn + 0.5 ) * dy )  THEN
243
244                   pos_x = psl(i)
245
246                   DO WHILE ( pos_x <= psr(i) )
247
248                      IF ( pos_x >= ( nxl - 0.5 ) * dx  .AND.  &
249                           pos_x <  ( nxr + 0.5 ) * dx )  THEN
250
251                         DO  j = 1, particles_per_point
252
253                            n = n + 1
254                            IF ( n > maximum_number_of_particles )  THEN
255                               PRINT*,'+++ init_particles: number of initial', &
256                                      ' particles (', n, ') exceeds'
257                               PRINT*,'    maximum_number_of_particles (',     &
258                                      maximum_number_of_particles, ') on PE ', &
259                                      myid
260#if defined( __parallel )
261                               CALL MPI_ABORT( comm2d, 9999, ierr )
262#else
263                               CALL local_stop
264#endif
265                            ENDIF
266                            particles(n)%x             = pos_x
267                            particles(n)%y             = pos_y
268                            particles(n)%z             = pos_z
269                            particles(n)%age           = 0.0
270                            particles(n)%age_m         = 0.0
271                            particles(n)%dt_sum        = 0.0
272                            particles(n)%dvrp_psize    = dvrp_psize
273                            particles(n)%e_m           = 0.0
274                            particles(n)%speed_x       = 0.0
275                            particles(n)%speed_x_sgs   = 0.0
276                            particles(n)%speed_y       = 0.0
277                            particles(n)%speed_y_sgs   = 0.0
278                            particles(n)%speed_z       = 0.0
279                            particles(n)%speed_z_sgs   = 0.0
280                            particles(n)%origin_x      = pos_x
281                            particles(n)%origin_y      = pos_y
282                            particles(n)%origin_z      = pos_z
283                            particles(n)%radius      = particle_groups(i)%radius
284                            particles(n)%weight_factor =initial_weighting_factor
285                            particles(n)%color         = 1
286                            particles(n)%group         = i
287                            particles(n)%tailpoints    = 0
288                            IF ( use_particle_tails  .AND. &
289                                 MOD( n, skip_particles_for_tail ) == 0 )  THEN
290                               number_of_tails         = number_of_tails + 1
291!
292!--                            This is a temporary provisional setting (see
293!--                            further below!)
294                               particles(n)%tail_id    = number_of_tails
295                            ELSE
296                               particles(n)%tail_id    = 0
297                            ENDIF
298
299                         ENDDO
300
301                      ENDIF
302
303                      pos_x = pos_x + pdx(i)
304
305                   ENDDO
306
307                ENDIF
308
309                pos_y = pos_y + pdy(i)
310
311             ENDDO
312
313             pos_z = pos_z + pdz(i)
314
315          ENDDO
316
317       ENDDO
318
319       number_of_initial_particles = n
320       number_of_particles         = n
321
322!
323!--    Calculate the number of particles and tails of the total domain
324#if defined( __parallel )
325       CALL MPI_ALLREDUCE( number_of_particles, total_number_of_particles, 1, &
326                           MPI_INTEGER, MPI_SUM, comm2d, ierr )
327       CALL MPI_ALLREDUCE( number_of_tails, total_number_of_tails, 1, &
328                           MPI_INTEGER, MPI_SUM, comm2d, ierr )
329#else
330       total_number_of_particles = number_of_particles
331       total_number_of_tails     = number_of_tails
332#endif
333
334!
335!--    Set a seed value for the random number generator to be exclusively
336!--    used for the particle code. The generated random numbers should be
337!--    different on the different PEs.
338       iran_part = iran_part + myid
339
340!
341!--    User modification of initial particles
342       CALL user_init_particles
343
344!
345!--    Store the initial set of particles for release at later times
346       IF ( number_of_initial_particles /= 0 )  THEN
347          ALLOCATE( initial_particles(1:number_of_initial_particles) )
348          initial_particles(1:number_of_initial_particles) = &
349                                        particles(1:number_of_initial_particles)
350       ENDIF
351
352!
353!--    Add random fluctuation to particle positions
354       IF ( random_start_position )  THEN
355
356          iran = iran + myid    ! Random positions should be different on
357                                ! different PEs
358
359          DO  n = 1, number_of_initial_particles
360             IF ( psl(particles(n)%group) /= psr(particles(n)%group) )  THEN
361                particles(n)%x = particles(n)%x + &
362                                 ( random_function( iran ) - 0.5 ) * &
363                                 pdx(particles(n)%group)
364                IF ( particles(n)%x  <=  ( nxl - 0.5 ) * dx )  THEN
365                   particles(n)%x = ( nxl - 0.4999999999 ) * dx
366                ELSEIF ( particles(n)%x  >=  ( nxr + 0.5 ) * dx )  THEN
367                   particles(n)%x = ( nxr + 0.4999999999 ) * dx
368                ENDIF
369             ENDIF
370             IF ( pss(particles(n)%group) /= psn(particles(n)%group) )  THEN
371                particles(n)%y = particles(n)%y + &
372                                 ( random_function( iran ) - 0.5 ) * &
373                                 pdy(particles(n)%group)
374                IF ( particles(n)%y  <=  ( nys - 0.5 ) * dy )  THEN
375                   particles(n)%y = ( nys - 0.4999999999 ) * dy
376                ELSEIF ( particles(n)%y  >=  ( nyn + 0.5 ) * dy )  THEN
377                   particles(n)%y = ( nyn + 0.4999999999 ) * dy
378                ENDIF
379             ENDIF
380             IF ( psb(particles(n)%group) /= pst(particles(n)%group) )  THEN
381                particles(n)%z = particles(n)%z + &
382                                 ( random_function( iran ) - 0.5 ) * &
383                                 pdz(particles(n)%group)
384             ENDIF
385          ENDDO
386       ENDIF
387
388!
389!--    Sort particles in the sequence the gridboxes are stored in the memory
390       CALL sort_particles
391
392!
393!--    Open file for statistical informations about particle conditions
394       IF ( write_particle_statistics )  THEN
395          CALL check_open( 80 )
396          WRITE ( 80, 8000 )  current_timestep_number, simulated_time, &
397                              number_of_initial_particles,             &
398                              maximum_number_of_particles
399          CALL close_file( 80 )
400       ENDIF
401
402!
403!--    Check if particles are really uniform in color and radius (dvrp_size)
404!--    (uniform_particles is preset TRUE)
405       IF ( uniform_particles )  THEN
406          IF ( number_of_initial_particles == 0 )  THEN
407             uniform_particles_l = .TRUE.
408          ELSE
409             n = number_of_initial_particles
410             IF ( MINVAL( particles(1:n)%dvrp_psize  ) ==     &
411                  MAXVAL( particles(1:n)%dvrp_psize  )  .AND. &
412                  MINVAL( particles(1:n)%color ) ==     &
413                  MAXVAL( particles(1:n)%color ) )  THEN
414                uniform_particles_l = .TRUE.
415             ELSE
416                uniform_particles_l = .FALSE.
417             ENDIF
418          ENDIF
419
420#if defined( __parallel )
421          CALL MPI_ALLREDUCE( uniform_particles_l, uniform_particles, 1, &
422                              MPI_LOGICAL, MPI_LAND, comm2d, ierr )
423#else
424          uniform_particles = uniform_particles_l
425#endif
426
427       ENDIF
428
429!
430!--    Set the beginning of the particle tails and their age
431       IF ( use_particle_tails )  THEN
432!
433!--       Choose the maximum number of tails significantly larger than the
434!--       one initially required
435          factor = 10.0
436          value  = number_of_tails
437          DO WHILE ( value / 10.0 >= 1.0 )
438             factor = factor * 10.0
439             value  = value / 10.0
440          ENDDO
441          maximum_number_of_tails = factor * INT( value )
442
443          ALLOCATE( particle_tail_coordinates(maximum_number_of_tailpoints,5, &
444                    maximum_number_of_tails),                                 &
445                    new_tail_id(maximum_number_of_tails),                     &
446                    tail_mask(maximum_number_of_tails) )
447
448          particle_tail_coordinates  = 0.0
449          minimum_tailpoint_distance = minimum_tailpoint_distance**2
450          number_of_initial_tails    = number_of_tails
451
452          nn = 0
453          DO  n = 1, number_of_particles
454!
455!--          Only for those particles marked above with a provisional tail_id
456!--          tails will be created. Particles now get their final tail_id.
457             IF ( particles(n)%tail_id /= 0 )  THEN
458
459                nn = nn + 1
460                particles(n)%tail_id = nn
461
462                particle_tail_coordinates(1,1,nn) = particles(n)%x
463                particle_tail_coordinates(1,2,nn) = particles(n)%y
464                particle_tail_coordinates(1,3,nn) = particles(n)%z
465                particle_tail_coordinates(1,4,nn) = particles(n)%color
466                particles(n)%tailpoints = 1
467                IF ( minimum_tailpoint_distance /= 0.0 )  THEN
468                   particle_tail_coordinates(2,1,nn) = particles(n)%x
469                   particle_tail_coordinates(2,2,nn) = particles(n)%y
470                   particle_tail_coordinates(2,3,nn) = particles(n)%z
471                   particle_tail_coordinates(2,4,nn) = particles(n)%color
472                   particle_tail_coordinates(1:2,5,nn) = 0.0
473                   particles(n)%tailpoints = 2
474                ENDIF
475
476             ENDIF
477          ENDDO
478       ENDIF
479
480!
481!--    Plot initial positions of particles (only if particle advection is
482!--    switched on from the beginning of the simulation (t=0))
483       IF ( particle_advection_start == 0.0 )  CALL data_output_dvrp
484
485    ENDIF
486
487!
488!-- Check boundary condition and set internal variables
489    SELECT CASE ( bc_par_b )
490   
491       CASE ( 'absorb' )
492          ibc_par_b = 1
493
494       CASE ( 'reflect' )
495          ibc_par_b = 2
496         
497       CASE DEFAULT
498          IF ( myid == 0 )  THEN
499             PRINT*,'+++ init_particles: unknown boundary condition ',   &
500                         'bc_par_b = "', TRIM( bc_par_b ), '"'
501          ENDIF
502          CALL local_stop
503         
504    END SELECT
505    SELECT CASE ( bc_par_t )
506   
507       CASE ( 'absorb' )
508          ibc_par_t = 1
509
510       CASE ( 'reflect' )
511          ibc_par_t = 2
512         
513       CASE DEFAULT
514          IF ( myid == 0 )  THEN
515             PRINT*,'+++ init_particles: unknown boundary condition ',   &
516                         'bc_par_t = "', TRIM( bc_par_t ), '"'
517          ENDIF
518          CALL local_stop
519         
520    END SELECT
521    SELECT CASE ( bc_par_lr )
522
523       CASE ( 'cyclic' )
524          ibc_par_lr = 0
525
526       CASE ( 'absorb' )
527          ibc_par_lr = 1
528
529       CASE ( 'reflect' )
530          ibc_par_lr = 2
531         
532       CASE DEFAULT
533          IF ( myid == 0 )  THEN
534             PRINT*,'+++ init_particles: unknown boundary condition ',   &
535                         'bc_par_lr = "', TRIM( bc_par_lr ), '"'
536          ENDIF
537          CALL local_stop
538         
539    END SELECT
540    SELECT CASE ( bc_par_ns )
541
542       CASE ( 'cyclic' )
543          ibc_par_ns = 0
544
545       CASE ( 'absorb' )
546          ibc_par_ns = 1
547
548       CASE ( 'reflect' )
549          ibc_par_ns = 2
550         
551       CASE DEFAULT
552          IF ( myid == 0 )  THEN
553             PRINT*,'+++ init_particles: unknown boundary condition ',   &
554                         'bc_par_ns = "', TRIM( bc_par_ns ), '"'
555          ENDIF
556          CALL local_stop
557         
558    END SELECT
559!
560!-- Formats
5618000 FORMAT (I6,1X,F7.2,4X,I6,71X,I6)
562
563 END SUBROUTINE init_particles
Note: See TracBrowser for help on using the repository browser.