Home

Context Navigation

source: palm/trunk/SOURCE/lpm_init.f90 @ 1344

Last change on this file since 1344 was 1329, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 26.1 KB

Rev	Line
[849]	1	SUBROUTINE lpm_init
[1]	2
[1036]	3	!--------------------------------------------------------------------------------!
	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[1310]	17	! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]	18	!--------------------------------------------------------------------------------!
	19	!
[254]	20	! Current revisions:
[1]	21	! -----------------
[1329]	22	!
	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: lpm_init.f90 1329 2014-03-21 11:09:15Z kanani $
	27	!
[1329]	28	! 1327 2014-03-21 11:00:16Z raasch
	29	! -netcdf_output
	30	!
[1323]	31	! 1322 2014-03-20 16:38:49Z raasch
	32	! REAL functions provided with KIND-attribute
	33	!
[1321]	34	! 1320 2014-03-20 08:40:49Z raasch
[1320]	35	! ONLY-attribute added to USE-statements,
	36	! kind-parameters added to all INTEGER and REAL declaration statements,
	37	! kinds are defined in new module kinds,
	38	! revision history before 2012 removed,
	39	! comment fields (!:) to be used for variable explanations added to
	40	! all variable declaration statements
	41	! bugfix: #if defined( __parallel ) added
[850]	42	!
[1315]	43	! 1314 2014-03-14 18:25:17Z suehring
	44	! Vertical logarithmic interpolation of horizontal particle speed for particles
	45	! between roughness height and first vertical grid level.
	46	!
[1093]	47	! 1092 2013-02-02 11:24:22Z raasch
	48	! unused variables removed
	49	!
[1037]	50	! 1036 2012-10-22 13:43:42Z raasch
	51	! code put under GPL (PALM 3.9)
	52	!
[850]	53	! 849 2012-03-15 10:35:09Z raasch
[849]	54	! routine renamed: init_particles -> lpm_init
	55	! de_dx, de_dy, de_dz are allocated here (instead of automatic arrays in
	56	! advec_particles),
	57	! sort_particles renamed lpm_sort_arrays, user_init_particles renamed lpm_init
[392]	58	!
[829]	59	! 828 2012-02-21 12:00:36Z raasch
	60	! call of init_kernels, particle feature color renamed class
	61	!
[826]	62	! 824 2012-02-17 09:09:57Z raasch
	63	! particle attributes speed_x\|y\|z_sgs renamed rvar1\|2\|3,
	64	! array particles implemented as pointer
	65	!
[668]	66	! 667 2010-12-23 12:06:00Z suehring/gryschka
	67	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation
	68	! of arrays.
	69	!
[1]	70	! Revision 1.1 1999/11/25 16:22:38 raasch
	71	! Initial revision
	72	!
	73	!
	74	! Description:
	75	! ------------
	76	! This routine initializes a set of particles and their attributes (position,
[849]	77	! radius, ..) which are used by the Lagrangian particle model (see lpm).
[1]	78	!------------------------------------------------------------------------------!
	79
[1320]	80	USE arrays_3d, &
	81	ONLY: de_dx, de_dy, de_dz, zu, zw, z0
	82
	83	USE cloud_parameters, &
	84	ONLY: curvature_solution_effects
	85
	86	USE control_parameters, &
	87	ONLY: cloud_droplets, current_timestep_number, initializing_actions, &
[1327]	88	message_string, netcdf_data_format, ocean, &
[1320]	89	prandtl_layer, simulated_time
	90
	91	USE dvrp_variables, &
	92	ONLY: particle_color, particle_dvrpsize
	93
	94	USE grid_variables, &
	95	ONLY: dx, dy
	96
	97	USE indices, &
	98	ONLY: nx, nxl, nxlg, nxrg, nxr, ny, nyn, nys, nyng, nysg, nz, nzb, nzt
	99
	100	USE kinds
	101
	102	USE lpm_collision_kernels_mod, &
	103	ONLY: init_kernels
	104
	105	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, &
	113	number_of_initial_tails, number_of_particles, &
	114	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, &
	118	particle_tail_coordinates, particle_type, pdx, pdy, pdz, &
	119	prt_count, prt_start_index, psb, psl, psn, psr, pss, pst, &
	120	radius, random_start_position, read_particles_from_restartfile,&
	121	skip_particles_for_tail, sort_count, tail_mask, &
	122	total_number_of_particles, total_number_of_tails, &
	123	use_particle_tails, use_sgs_for_particles, &
	124	write_particle_statistics, uniform_particles, z0_av_global
	125
[1]	126	USE pegrid
	127
[1320]	128	USE random_function_mod, &
	129	ONLY: random_function
[1]	130
[1320]	131
[1]	132	IMPLICIT NONE
	133
[1320]	134	INTEGER(iwp) :: i !:
	135	INTEGER(iwp) :: j !:
	136	INTEGER(iwp) :: k !:
	137	INTEGER(iwp) :: n !:
	138	INTEGER(iwp) :: nn !:
	139
[1]	140	#if defined( __parallel )
[1320]	141	INTEGER(iwp), DIMENSION(3) :: blocklengths !:
	142	INTEGER(iwp), DIMENSION(3) :: displacements !:
	143	INTEGER(iwp), DIMENSION(3) :: types !:
[1]	144	#endif
[1320]	145	LOGICAL :: uniform_particles_l !:
[1]	146
[1320]	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 !:
[1]	154
[1320]	155
	156
	157
[1]	158	#if defined( __parallel )
	159	!
	160	!-- Define MPI derived datatype for FORTRAN datatype particle_type (see module
[82]	161	!-- 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
	164
[1]	165	types(1) = MPI_REAL
	166	types(2) = MPI_INTEGER
	167	types(3) = MPI_UB
	168	CALL MPI_TYPE_STRUCT( 3, blocklengths, displacements, types, &
	169	mpi_particle_type, ierr )
	170	CALL MPI_TYPE_COMMIT( mpi_particle_type, ierr )
	171	#endif
	172
	173	!
[150]	174	!-- In case of oceans runs, the vertical index calculations need an offset,
	175	!-- because otherwise the k indices will become negative
	176	IF ( ocean ) THEN
	177	offset_ocean_nzt = nzt
	178	offset_ocean_nzt_m1 = nzt - 1
	179	ENDIF
	180
	181
	182	!
[1]	183	!-- Check the number of particle groups.
	184	IF ( number_of_particle_groups > max_number_of_particle_groups ) THEN
[274]	185	WRITE( message_string, * ) 'max_number_of_particle_groups =', &
	186	max_number_of_particle_groups , &
[254]	187	'&number_of_particle_groups reset to ', &
	188	max_number_of_particle_groups
[849]	189	CALL message( 'lpm_init', 'PA0213', 0, 1, 0, 6, 0 )
[1]	190	number_of_particle_groups = max_number_of_particle_groups
	191	ENDIF
	192
	193	!
	194	!-- Set default start positions, if necessary
	195	IF ( psl(1) == 9999999.9 ) psl(1) = -0.5 * dx
	196	IF ( psr(1) == 9999999.9 ) psr(1) = ( nx + 0.5 ) * dx
	197	IF ( pss(1) == 9999999.9 ) pss(1) = -0.5 * dy
	198	IF ( psn(1) == 9999999.9 ) psn(1) = ( ny + 0.5 ) * dy
	199	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) = dx
	203	IF ( pdy(1) == 9999999.9 .OR. pdy(1) == 0.0 ) pdy(1) = dy
	204	IF ( pdz(1) == 9999999.9 .OR. pdz(1) == 0.0 ) pdz(1) = zu(2) - zu(1)
	205
	206	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)
	216	ENDDO
	217
	218	!
[849]	219	!-- Allocate arrays required for calculating particle SGS velocities
	220	IF ( use_sgs_for_particles ) THEN
	221	ALLOCATE( de_dx(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	222	de_dy(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	223	de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	224	ENDIF
	225
	226	!
[1314]	227	!-- Allocate array required for logarithmic vertical interpolation of
	228	!-- horizontal particle velocities between the surface and the first vertical
	229	!-- grid level. In order to avoid repeated CPU cost-intensive CALLS of
	230	!-- intrinsic FORTRAN procedure LOG(z/z0), LOG(z/z0) is precalculated for
	231	!-- several heights. Splitting into 20 sublayers turned out to be sufficient.
	232	!-- To obtain exact height levels of particles, linear interpolation is applied
	233	!-- (see lpm_advec.f90).
	234	IF ( prandtl_layer ) THEN
	235
	236	ALLOCATE ( log_z_z0(0:number_of_sublayers) )
	237	z_p = zu(nzb+1) - zw(nzb)
	238
	239	!
	240	!-- Calculate horizontal mean value of z0 used for logartihmic
	241	!-- interpolation. Note: this is not exact for heterogeneous z0.
	242	!-- However, sensitivity studies showed that the effect is
	243	!-- negligible.
	244	z0_av_local = SUM( z0(nys:nyn,nxl:nxr) )
	245	z0_av_global = 0.0
	246
[1320]	247	#if defined( __parallel )
[1314]	248	CALL MPI_ALLREDUCE(z0_av_local, z0_av_global, 1, MPI_REAL, MPI_SUM, &
	249	comm2d, ierr )
[1320]	250	#else
	251	z0_av_global = z0_av_local
	252	#endif
[1314]	253
	254	z0_av_global = z0_av_global / ( ( ny + 1 ) * ( nx + 1 ) )
	255	!
	256	!-- Horizontal wind speed is zero below and at z0
	257	log_z_z0(0) = 0.0
	258	!
	259	!-- Calculate vertical depth of the sublayers
[1322]	260	height_int = ( z_p - z0_av_global ) / REAL( number_of_sublayers, KIND=wp )
[1314]	261	!
	262	!-- Precalculate LOG(z/z0)
	263	height_p = 0.0
	264	DO k = 1, number_of_sublayers
	265
	266	height_p = height_p + height_int
	267	log_z_z0(k) = LOG( height_p / z0_av_global )
	268
	269	ENDDO
	270
	271
	272	ENDIF
	273
	274	!
[828]	275	!-- Initialize collision kernels
	276	IF ( collision_kernel /= 'none' ) CALL init_kernels
	277
	278	!
[1]	279	!-- For the first model run of a possible job chain initialize the
[849]	280	!-- particles, otherwise read the particle data from restart file.
[1]	281	IF ( TRIM( initializing_actions ) == 'read_restart_data' &
	282	.AND. read_particles_from_restartfile ) THEN
	283
[849]	284	CALL lpm_read_restart_file
[1]	285
	286	ELSE
	287
	288	!
	289	!-- Allocate particle arrays and set attributes of the initial set of
	290	!-- particles, which can be also periodically released at later times.
	291	!-- Also allocate array for particle tail coordinates, if needed.
[667]	292	ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	293	prt_start_index(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
[792]	294	particle_mask(maximum_number_of_particles), &
	295	part_1(maximum_number_of_particles), &
	296	part_2(maximum_number_of_particles) )
[1]	297
[792]	298	particles => part_1
	299
	300	sort_count = 0
	301
[1]	302	!
	303	!-- Initialize all particles with dummy values (otherwise errors may
	304	!-- occur within restart runs). The reason for this is still not clear
	305	!-- 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, &
[57]	308	0.0, 0, 0, 0, 0 )
[1]	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 plots
	313	IF ( dvrp_psize == 9999999.9 ) dvrp_psize = 0.2 * dx
	314
	315	!
	316	!-- Set values for the density ratio and radius for all particle
	317	!-- groups, if necessary
	318	IF ( density_ratio(1) == 9999999.9 ) density_ratio(1) = 0.0
	319	IF ( radius(1) == 9999999.9 ) radius(1) = 0.0
	320	DO i = 2, number_of_particle_groups
	321	IF ( density_ratio(i) == 9999999.9 ) THEN
	322	density_ratio(i) = density_ratio(i-1)
	323	ENDIF
	324	IF ( radius(i) == 9999999.9 ) radius(i) = radius(i-1)
	325	ENDDO
	326
	327	DO i = 1, number_of_particle_groups
	328	IF ( density_ratio(i) /= 0.0 .AND. radius(i) == 0 ) THEN
[254]	329	WRITE( message_string, * ) 'particle group #', i, 'has a', &
	330	'density ratio /= 0 but radius = 0'
[849]	331	CALL message( 'lpm_init', 'PA0215', 1, 2, 0, 6, 0 )
[1]	332	ENDIF
	333	particle_groups(i)%density_ratio = density_ratio(i)
	334	particle_groups(i)%radius = radius(i)
	335	ENDDO
	336
	337	!
	338	!-- Calculate particle positions and store particle attributes, if
	339	!-- particle is situated on this PE
	340	n = 0
	341
	342	DO i = 1, number_of_particle_groups
	343
	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 ) THEN
	354
	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 ) THEN
	361
	362	DO j = 1, particles_per_point
	363
	364	n = n + 1
	365	IF ( n > maximum_number_of_particles ) THEN
[254]	366	WRITE( message_string, * ) 'number of initial', &
[274]	367	'particles (', n, ') exceeds', &
	368	'&maximum_number_of_particles (', &
	369	maximum_number_of_particles, ') on PE ', &
[254]	370	myid
[849]	371	CALL message( 'lpm_init', 'PA0216', 2, 2, -1, 6,&
	372	1 )
[1]	373	ENDIF
	374	particles(n)%x = pos_x
	375	particles(n)%y = pos_y
	376	particles(n)%z = pos_z
	377	particles(n)%age = 0.0
[57]	378	particles(n)%age_m = 0.0
[1]	379	particles(n)%dt_sum = 0.0
	380	particles(n)%dvrp_psize = dvrp_psize
	381	particles(n)%e_m = 0.0
[824]	382	IF ( curvature_solution_effects ) THEN
	383	!
	384	!-- Initial values (internal timesteps, derivative)
	385	!-- for Rosenbrock method
	386	particles(n)%rvar1 = 1.0E-12
	387	particles(n)%rvar2 = 1.0E-3
	388	particles(n)%rvar3 = -9999999.9
	389	ELSE
	390	!
	391	!-- Initial values for SGS velocities
	392	particles(n)%rvar1 = 0.0
	393	particles(n)%rvar2 = 0.0
	394	particles(n)%rvar3 = 0.0
	395	ENDIF
[1]	396	particles(n)%speed_x = 0.0
	397	particles(n)%speed_y = 0.0
	398	particles(n)%speed_z = 0.0
	399	particles(n)%origin_x = pos_x
	400	particles(n)%origin_y = pos_y
	401	particles(n)%origin_z = pos_z
	402	particles(n)%radius = particle_groups(i)%radius
	403	particles(n)%weight_factor =initial_weighting_factor
[828]	404	particles(n)%class = 1
[1]	405	particles(n)%group = i
	406	particles(n)%tailpoints = 0
	407	IF ( use_particle_tails .AND. &
	408	MOD( n, skip_particles_for_tail ) == 0 ) THEN
	409	number_of_tails = number_of_tails + 1
	410	!
	411	!-- This is a temporary provisional setting (see
	412	!-- further below!)
	413	particles(n)%tail_id = number_of_tails
	414	ELSE
	415	particles(n)%tail_id = 0
	416	ENDIF
	417
	418	ENDDO
	419
	420	ENDIF
	421
	422	pos_x = pos_x + pdx(i)
	423
	424	ENDDO
	425
	426	ENDIF
	427
	428	pos_y = pos_y + pdy(i)
	429
	430	ENDDO
	431
	432	pos_z = pos_z + pdz(i)
	433
	434	ENDDO
	435
	436	ENDDO
	437
	438	number_of_initial_particles = n
	439	number_of_particles = n
	440
	441	!
	442	!-- Calculate the number of particles and tails of the total domain
	443	#if defined( __parallel )
[622]	444	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	445	CALL MPI_ALLREDUCE( number_of_particles, total_number_of_particles, 1, &
[16]	446	MPI_INTEGER, MPI_SUM, comm2d, ierr )
[622]	447	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	448	CALL MPI_ALLREDUCE( number_of_tails, total_number_of_tails, 1, &
[16]	449	MPI_INTEGER, MPI_SUM, comm2d, ierr )
[1]	450	#else
	451	total_number_of_particles = number_of_particles
	452	total_number_of_tails = number_of_tails
	453	#endif
	454
	455	!
	456	!-- Set a seed value for the random number generator to be exclusively
	457	!-- used for the particle code. The generated random numbers should be
	458	!-- different on the different PEs.
	459	iran_part = iran_part + myid
	460
	461	!
	462	!-- User modification of initial particles
[849]	463	CALL user_lpm_init
[1]	464
	465	!
	466	!-- Store the initial set of particles for release at later times
	467	IF ( number_of_initial_particles /= 0 ) THEN
	468	ALLOCATE( initial_particles(1:number_of_initial_particles) )
	469	initial_particles(1:number_of_initial_particles) = &
	470	particles(1:number_of_initial_particles)
	471	ENDIF
	472
	473	!
	474	!-- Add random fluctuation to particle positions
	475	IF ( random_start_position ) THEN
	476
	477	DO n = 1, number_of_initial_particles
	478	IF ( psl(particles(n)%group) /= psr(particles(n)%group) ) THEN
	479	particles(n)%x = particles(n)%x + &
[106]	480	( random_function( iran_part ) - 0.5 ) * &
[1]	481	pdx(particles(n)%group)
	482	IF ( particles(n)%x <= ( nxl - 0.5 ) * dx ) THEN
	483	particles(n)%x = ( nxl - 0.4999999999 ) * dx
	484	ELSEIF ( particles(n)%x >= ( nxr + 0.5 ) * dx ) THEN
	485	particles(n)%x = ( nxr + 0.4999999999 ) * dx
	486	ENDIF
	487	ENDIF
	488	IF ( pss(particles(n)%group) /= psn(particles(n)%group) ) THEN
	489	particles(n)%y = particles(n)%y + &
[106]	490	( random_function( iran_part ) - 0.5 ) * &
[1]	491	pdy(particles(n)%group)
	492	IF ( particles(n)%y <= ( nys - 0.5 ) * dy ) THEN
	493	particles(n)%y = ( nys - 0.4999999999 ) * dy
	494	ELSEIF ( particles(n)%y >= ( nyn + 0.5 ) * dy ) THEN
	495	particles(n)%y = ( nyn + 0.4999999999 ) * dy
	496	ENDIF
	497	ENDIF
	498	IF ( psb(particles(n)%group) /= pst(particles(n)%group) ) THEN
	499	particles(n)%z = particles(n)%z + &
[106]	500	( random_function( iran_part ) - 0.5 ) * &
[1]	501	pdz(particles(n)%group)
	502	ENDIF
	503	ENDDO
	504	ENDIF
	505
	506	!
[117]	507	!-- Sort particles in the sequence the gridboxes are stored in the memory.
	508	!-- Only required if cloud droplets are used.
[849]	509	IF ( cloud_droplets ) CALL lpm_sort_arrays
[1]	510
	511	!
	512	!-- Open file for statistical informations about particle conditions
	513	IF ( write_particle_statistics ) THEN
	514	CALL check_open( 80 )
	515	WRITE ( 80, 8000 ) current_timestep_number, simulated_time, &
	516	number_of_initial_particles, &
	517	maximum_number_of_particles
	518	CALL close_file( 80 )
	519	ENDIF
	520
	521	!
	522	!-- Check if particles are really uniform in color and radius (dvrp_size)
	523	!-- (uniform_particles is preset TRUE)
	524	IF ( uniform_particles ) THEN
	525	IF ( number_of_initial_particles == 0 ) THEN
	526	uniform_particles_l = .TRUE.
	527	ELSE
	528	n = number_of_initial_particles
	529	IF ( MINVAL( particles(1:n)%dvrp_psize ) == &
	530	MAXVAL( particles(1:n)%dvrp_psize ) .AND. &
[828]	531	MINVAL( particles(1:n)%class ) == &
	532	MAXVAL( particles(1:n)%class ) ) THEN
[1]	533	uniform_particles_l = .TRUE.
	534	ELSE
	535	uniform_particles_l = .FALSE.
	536	ENDIF
	537	ENDIF
	538
	539	#if defined( __parallel )
[622]	540	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	541	CALL MPI_ALLREDUCE( uniform_particles_l, uniform_particles, 1, &
	542	MPI_LOGICAL, MPI_LAND, comm2d, ierr )
	543	#else
	544	uniform_particles = uniform_particles_l
	545	#endif
	546
	547	ENDIF
	548
	549	!
[336]	550	!-- Particles will probably become none-uniform, if their size and color
	551	!-- will be determined by flow variables
	552	IF ( particle_color /= 'none' .OR. particle_dvrpsize /= 'none' ) THEN
	553	uniform_particles = .FALSE.
	554	ENDIF
	555
	556	!
[1]	557	!-- Set the beginning of the particle tails and their age
	558	IF ( use_particle_tails ) THEN
	559	!
[276]	560	!-- Choose the maximum number of tails with respect to the maximum number
	561	!-- of particles and skip_particles_for_tail
	562	maximum_number_of_tails = maximum_number_of_particles / &
	563	skip_particles_for_tail
	564
[229]	565	!
	566	!-- Create a minimum number of tails in case that there is no tail
	567	!-- initially (otherwise, index errors will occur when adressing the
	568	!-- arrays below)
	569	IF ( maximum_number_of_tails == 0 ) maximum_number_of_tails = 10
[1]	570
	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.0
	577	minimum_tailpoint_distance = minimum_tailpoint_distance**2
	578	number_of_initial_tails = number_of_tails
	579
	580	nn = 0
	581	DO n = 1, number_of_particles
	582	!
	583	!-- Only for those particles marked above with a provisional tail_id
	584	!-- tails will be created. Particles now get their final tail_id.
	585	IF ( particles(n)%tail_id /= 0 ) THEN
	586
	587	nn = nn + 1
	588	particles(n)%tail_id = nn
	589
	590	particle_tail_coordinates(1,1,nn) = particles(n)%x
	591	particle_tail_coordinates(1,2,nn) = particles(n)%y
	592	particle_tail_coordinates(1,3,nn) = particles(n)%z
[828]	593	particle_tail_coordinates(1,4,nn) = particles(n)%class
[1]	594	particles(n)%tailpoints = 1
	595	IF ( minimum_tailpoint_distance /= 0.0 ) THEN
	596	particle_tail_coordinates(2,1,nn) = particles(n)%x
	597	particle_tail_coordinates(2,2,nn) = particles(n)%y
	598	particle_tail_coordinates(2,3,nn) = particles(n)%z
[828]	599	particle_tail_coordinates(2,4,nn) = particles(n)%class
[1]	600	particle_tail_coordinates(1:2,5,nn) = 0.0
	601	particles(n)%tailpoints = 2
	602	ENDIF
	603
	604	ENDIF
	605	ENDDO
	606	ENDIF
	607
	608	!
	609	!-- Plot initial positions of particles (only if particle advection is
	610	!-- switched on from the beginning of the simulation (t=0))
	611	IF ( particle_advection_start == 0.0 ) CALL data_output_dvrp
	612
	613	ENDIF
	614
	615	!
	616	!-- Check boundary condition and set internal variables
	617	SELECT CASE ( bc_par_b )
	618
	619	CASE ( 'absorb' )
	620	ibc_par_b = 1
	621
	622	CASE ( 'reflect' )
	623	ibc_par_b = 2
	624
	625	CASE DEFAULT
[254]	626	WRITE( message_string, * ) 'unknown boundary condition ', &
	627	'bc_par_b = "', TRIM( bc_par_b ), '"'
[849]	628	CALL message( 'lpm_init', 'PA0217', 1, 2, 0, 6, 0 )
[1]	629
	630	END SELECT
	631	SELECT CASE ( bc_par_t )
	632
	633	CASE ( 'absorb' )
	634	ibc_par_t = 1
	635
	636	CASE ( 'reflect' )
	637	ibc_par_t = 2
	638
	639	CASE DEFAULT
[254]	640	WRITE( message_string, * ) 'unknown boundary condition ', &
	641	'bc_par_t = "', TRIM( bc_par_t ), '"'
[849]	642	CALL message( 'lpm_init', 'PA0218', 1, 2, 0, 6, 0 )
[1]	643
	644	END SELECT
	645	SELECT CASE ( bc_par_lr )
	646
	647	CASE ( 'cyclic' )
	648	ibc_par_lr = 0
	649
	650	CASE ( 'absorb' )
	651	ibc_par_lr = 1
	652
	653	CASE ( 'reflect' )
	654	ibc_par_lr = 2
	655
	656	CASE DEFAULT
[254]	657	WRITE( message_string, * ) 'unknown boundary condition ', &
	658	'bc_par_lr = "', TRIM( bc_par_lr ), '"'
[849]	659	CALL message( 'lpm_init', 'PA0219', 1, 2, 0, 6, 0 )
[1]	660
	661	END SELECT
	662	SELECT CASE ( bc_par_ns )
	663
	664	CASE ( 'cyclic' )
	665	ibc_par_ns = 0
	666
	667	CASE ( 'absorb' )
	668	ibc_par_ns = 1
	669
	670	CASE ( 'reflect' )
	671	ibc_par_ns = 2
	672
	673	CASE DEFAULT
[254]	674	WRITE( message_string, * ) 'unknown boundary condition ', &
	675	'bc_par_ns = "', TRIM( bc_par_ns ), '"'
[849]	676	CALL message( 'lpm_init', 'PA0220', 1, 2, 0, 6, 0 )
[1]	677
	678	END SELECT
	679	!
	680	!-- Formats
	681	8000 FORMAT (I6,1X,F7.2,4X,I6,71X,I6)
	682
[849]	683	END SUBROUTINE lpm_init

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |