Home

Context Navigation

init_particles.f90 @ 514

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

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:

lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:

Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list

Property svn:keywords set to Id

File size: 21.8 KB

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

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

Context Navigation

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

Download in other formats: