Home

Context Navigation

init_particles.f90 @ 92

Last change on this file since 92 was 83, checked in by raasch, 18 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
561	8000 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.

Context Navigation

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

Download in other formats: