source: palm/trunk/SOURCE/init_particles.f90 @ 4

Last change on this file since 4 was 4, checked in by raasch, 15 years ago

Id keyword set as property for all *.f90 files

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