Home

Context Navigation

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

Last change on this file since 1560 was 1360, checked in by hoffmann, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 34.9 KB

Rev	Line
[1359]	1	MODULE lpm_init_mod
[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	!
[1360]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: lpm_init.f90 1360 2014-04-11 17:20:32Z keck $
	27	!
[1360]	28	! 1359 2014-04-11 17:15:14Z hoffmann
	29	! New particle structure integrated.
	30	! Kind definition added to all floating point numbers.
	31	! lpm_init changed form a subroutine to a module.
	32	!
[1329]	33	! 1327 2014-03-21 11:00:16Z raasch
	34	! -netcdf_output
	35	!
[1323]	36	! 1322 2014-03-20 16:38:49Z raasch
	37	! REAL functions provided with KIND-attribute
	38	!
[1321]	39	! 1320 2014-03-20 08:40:49Z raasch
[1320]	40	! ONLY-attribute added to USE-statements,
	41	! kind-parameters added to all INTEGER and REAL declaration statements,
	42	! kinds are defined in new module kinds,
	43	! revision history before 2012 removed,
	44	! comment fields (!:) to be used for variable explanations added to
	45	! all variable declaration statements
	46	! bugfix: #if defined( __parallel ) added
[850]	47	!
[1315]	48	! 1314 2014-03-14 18:25:17Z suehring
	49	! Vertical logarithmic interpolation of horizontal particle speed for particles
	50	! between roughness height and first vertical grid level.
	51	!
[1093]	52	! 1092 2013-02-02 11:24:22Z raasch
	53	! unused variables removed
	54	!
[1037]	55	! 1036 2012-10-22 13:43:42Z raasch
	56	! code put under GPL (PALM 3.9)
	57	!
[850]	58	! 849 2012-03-15 10:35:09Z raasch
[849]	59	! routine renamed: init_particles -> lpm_init
	60	! de_dx, de_dy, de_dz are allocated here (instead of automatic arrays in
	61	! advec_particles),
	62	! sort_particles renamed lpm_sort_arrays, user_init_particles renamed lpm_init
[392]	63	!
[829]	64	! 828 2012-02-21 12:00:36Z raasch
	65	! call of init_kernels, particle feature color renamed class
	66	!
[826]	67	! 824 2012-02-17 09:09:57Z raasch
	68	! particle attributes speed_x\|y\|z_sgs renamed rvar1\|2\|3,
	69	! array particles implemented as pointer
	70	!
[668]	71	! 667 2010-12-23 12:06:00Z suehring/gryschka
	72	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation
	73	! of arrays.
	74	!
[1]	75	! Revision 1.1 1999/11/25 16:22:38 raasch
	76	! Initial revision
	77	!
	78	!
	79	! Description:
	80	! ------------
	81	! This routine initializes a set of particles and their attributes (position,
[849]	82	! radius, ..) which are used by the Lagrangian particle model (see lpm).
[1]	83	!------------------------------------------------------------------------------!
	84
[1320]	85	USE arrays_3d, &
	86	ONLY: de_dx, de_dy, de_dz, zu, zw, z0
	87
	88	USE cloud_parameters, &
	89	ONLY: curvature_solution_effects
	90
	91	USE control_parameters, &
[1359]	92	ONLY: cloud_droplets, current_timestep_number, dz, &
	93	initializing_actions, message_string, netcdf_data_format, &
	94	ocean, prandtl_layer, simulated_time
[1320]	95
	96	USE dvrp_variables, &
	97	ONLY: particle_color, particle_dvrpsize
	98
	99	USE grid_variables, &
[1359]	100	ONLY: ddx, dx, ddy, dy
[1320]	101
	102	USE indices, &
	103	ONLY: nx, nxl, nxlg, nxrg, nxr, ny, nyn, nys, nyng, nysg, nz, nzb, nzt
	104
	105	USE kinds
	106
	107	USE lpm_collision_kernels_mod, &
	108	ONLY: init_kernels
	109
	110	USE particle_attributes, &
[1359]	111	ONLY: alloc_factor, bc_par_b, bc_par_lr, bc_par_ns, bc_par_t, &
	112	block_offset, block_offset_def, collision_kernel, &
	113	density_ratio, dvrp_psize, grid_particles, &
	114	initial_weighting_factor, ibc_par_b, ibc_par_lr, ibc_par_ns, &
	115	ibc_par_t, iran_part, log_z_z0, &
	116	max_number_of_particle_groups, maximum_number_of_particles, &
	117	maximum_number_of_tailpoints, maximum_number_of_tails, &
	118	minimum_tailpoint_distance, min_nr_particle, &
	119	mpi_particle_type, new_tail_id, &
[1320]	120	number_of_initial_tails, number_of_particles, &
	121	number_of_particle_groups, number_of_sublayers, &
[1359]	122	number_of_tails, offset_ocean_nzt, offset_ocean_nzt_m1, &
	123	particles, particle_advection_start, particle_groups, &
	124	particle_groups_type, particles_per_point, &
[1320]	125	particle_tail_coordinates, particle_type, pdx, pdy, pdz, &
[1359]	126	prt_count, psb, psl, psn, psr, pss, pst, &
[1320]	127	radius, random_start_position, read_particles_from_restartfile,&
[1359]	128	skip_particles_for_tail, sort_count, &
	129	tail_mask, total_number_of_particles, total_number_of_tails, &
[1320]	130	use_particle_tails, use_sgs_for_particles, &
[1359]	131	write_particle_statistics, uniform_particles, zero_particle, &
	132	z0_av_global
[1320]	133
[1]	134	USE pegrid
	135
[1320]	136	USE random_function_mod, &
	137	ONLY: random_function
[1]	138
[1359]	139	IMPLICIT NONE
[1320]	140
[1359]	141	PRIVATE
	142
	143	INTEGER(iwp), PARAMETER :: PHASE_INIT = 1 !:
	144	INTEGER(iwp), PARAMETER, PUBLIC :: PHASE_RELEASE = 2 !:
	145
	146	INTERFACE lpm_init
	147	MODULE PROCEDURE lpm_init
	148	END INTERFACE lpm_init
	149
	150	INTERFACE lpm_create_particle
	151	MODULE PROCEDURE lpm_create_particle
	152	END INTERFACE lpm_create_particle
	153
	154	PUBLIC lpm_init, lpm_create_particle
	155
	156	CONTAINS
	157
	158	SUBROUTINE lpm_init
	159
	160	USE lpm_collision_kernels_mod, &
	161	ONLY: init_kernels
	162
[1]	163	IMPLICIT NONE
	164
[1320]	165	INTEGER(iwp) :: i !:
[1359]	166	INTEGER(iwp) :: ip !:
[1320]	167	INTEGER(iwp) :: j !:
[1359]	168	INTEGER(iwp) :: jp !:
[1320]	169	INTEGER(iwp) :: k !:
[1359]	170	INTEGER(iwp) :: kp !:
[1320]	171	INTEGER(iwp) :: n !:
	172	INTEGER(iwp) :: nn !:
	173
[1]	174	#if defined( __parallel )
[1320]	175	INTEGER(iwp), DIMENSION(3) :: blocklengths !:
	176	INTEGER(iwp), DIMENSION(3) :: displacements !:
	177	INTEGER(iwp), DIMENSION(3) :: types !:
[1]	178	#endif
[1320]	179	LOGICAL :: uniform_particles_l !:
[1]	180
[1359]	181	REAL(wp) :: height_int !:
	182	REAL(wp) :: height_p !:
	183	REAL(wp) :: z_p !:
	184	REAL(wp) :: z0_av_local !:
[1]	185
	186	#if defined( __parallel )
	187	!
	188	!-- Define MPI derived datatype for FORTRAN datatype particle_type (see module
[82]	189	!-- particle_attributes). Integer length is 4 byte, Real is 8 byte
[1359]	190	#if defined( __twocachelines )
	191	blocklengths(1) = 7; blocklengths(2) = 18; blocklengths(3) = 1
	192	displacements(1) = 0; displacements(2) = 64; displacements(3) = 128
[82]	193
[1359]	194	types(1) = MPI_REAL ! 64 bit words
	195	types(2) = MPI_INTEGER ! 32 Bit words
	196	types(3) = MPI_UB
	197	#else
	198	blocklengths(1) = 19; blocklengths(2) = 6; blocklengths(3) = 1
	199	displacements(1) = 0; displacements(2) = 152; displacements(3) = 176
	200
[1]	201	types(1) = MPI_REAL
	202	types(2) = MPI_INTEGER
	203	types(3) = MPI_UB
[1359]	204	#endif
[1]	205	CALL MPI_TYPE_STRUCT( 3, blocklengths, displacements, types, &
	206	mpi_particle_type, ierr )
	207	CALL MPI_TYPE_COMMIT( mpi_particle_type, ierr )
	208	#endif
	209
	210	!
[150]	211	!-- In case of oceans runs, the vertical index calculations need an offset,
	212	!-- because otherwise the k indices will become negative
	213	IF ( ocean ) THEN
	214	offset_ocean_nzt = nzt
	215	offset_ocean_nzt_m1 = nzt - 1
	216	ENDIF
	217
[1359]	218	!
	219	!-- Define block offsets for dividing a gridcell in 8 sub cells
[150]	220
[1359]	221	block_offset(0) = block_offset_def (-1,-1,-1)
	222	block_offset(1) = block_offset_def (-1,-1, 0)
	223	block_offset(2) = block_offset_def (-1, 0,-1)
	224	block_offset(3) = block_offset_def (-1, 0, 0)
	225	block_offset(4) = block_offset_def ( 0,-1,-1)
	226	block_offset(5) = block_offset_def ( 0,-1, 0)
	227	block_offset(6) = block_offset_def ( 0, 0,-1)
	228	block_offset(7) = block_offset_def ( 0, 0, 0)
[150]	229	!
[1]	230	!-- Check the number of particle groups.
	231	IF ( number_of_particle_groups > max_number_of_particle_groups ) THEN
[274]	232	WRITE( message_string, * ) 'max_number_of_particle_groups =', &
	233	max_number_of_particle_groups , &
[254]	234	'&number_of_particle_groups reset to ', &
	235	max_number_of_particle_groups
[849]	236	CALL message( 'lpm_init', 'PA0213', 0, 1, 0, 6, 0 )
[1]	237	number_of_particle_groups = max_number_of_particle_groups
	238	ENDIF
	239
	240	!
	241	!-- Set default start positions, if necessary
[1359]	242	IF ( psl(1) == 9999999.9_wp ) psl(1) = -0.5_wp * dx
	243	IF ( psr(1) == 9999999.9_wp ) psr(1) = ( nx + 0.5_wp ) * dx
	244	IF ( pss(1) == 9999999.9_wp ) pss(1) = -0.5_wp * dy
	245	IF ( psn(1) == 9999999.9_wp ) psn(1) = ( ny + 0.5_wp ) * dy
	246	IF ( psb(1) == 9999999.9_wp ) psb(1) = zu(nz/2)
	247	IF ( pst(1) == 9999999.9_wp ) pst(1) = psb(1)
[1]	248
[1359]	249	IF ( pdx(1) == 9999999.9_wp .OR. pdx(1) == 0.0_wp ) pdx(1) = dx
	250	IF ( pdy(1) == 9999999.9_wp .OR. pdy(1) == 0.0_wp ) pdy(1) = dy
	251	IF ( pdz(1) == 9999999.9_wp .OR. pdz(1) == 0.0_wp ) pdz(1) = zu(2) - zu(1)
[1]	252
	253	DO j = 2, number_of_particle_groups
[1359]	254	IF ( psl(j) == 9999999.9_wp ) psl(j) = psl(j-1)
	255	IF ( psr(j) == 9999999.9_wp ) psr(j) = psr(j-1)
	256	IF ( pss(j) == 9999999.9_wp ) pss(j) = pss(j-1)
	257	IF ( psn(j) == 9999999.9_wp ) psn(j) = psn(j-1)
	258	IF ( psb(j) == 9999999.9_wp ) psb(j) = psb(j-1)
	259	IF ( pst(j) == 9999999.9_wp ) pst(j) = pst(j-1)
	260	IF ( pdx(j) == 9999999.9_wp .OR. pdx(j) == 0.0_wp ) pdx(j) = pdx(j-1)
	261	IF ( pdy(j) == 9999999.9_wp .OR. pdy(j) == 0.0_wp ) pdy(j) = pdy(j-1)
	262	IF ( pdz(j) == 9999999.9_wp .OR. pdz(j) == 0.0_wp ) pdz(j) = pdz(j-1)
[1]	263	ENDDO
	264
	265	!
[849]	266	!-- Allocate arrays required for calculating particle SGS velocities
	267	IF ( use_sgs_for_particles ) THEN
	268	ALLOCATE( de_dx(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	269	de_dy(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	270	de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	271	ENDIF
	272
	273	!
[1314]	274	!-- Allocate array required for logarithmic vertical interpolation of
	275	!-- horizontal particle velocities between the surface and the first vertical
	276	!-- grid level. In order to avoid repeated CPU cost-intensive CALLS of
	277	!-- intrinsic FORTRAN procedure LOG(z/z0), LOG(z/z0) is precalculated for
	278	!-- several heights. Splitting into 20 sublayers turned out to be sufficient.
	279	!-- To obtain exact height levels of particles, linear interpolation is applied
	280	!-- (see lpm_advec.f90).
	281	IF ( prandtl_layer ) THEN
	282
	283	ALLOCATE ( log_z_z0(0:number_of_sublayers) )
	284	z_p = zu(nzb+1) - zw(nzb)
	285
	286	!
	287	!-- Calculate horizontal mean value of z0 used for logartihmic
	288	!-- interpolation. Note: this is not exact for heterogeneous z0.
	289	!-- However, sensitivity studies showed that the effect is
	290	!-- negligible.
	291	z0_av_local = SUM( z0(nys:nyn,nxl:nxr) )
[1359]	292	z0_av_global = 0.0_wp
[1314]	293
[1320]	294	#if defined( __parallel )
[1314]	295	CALL MPI_ALLREDUCE(z0_av_local, z0_av_global, 1, MPI_REAL, MPI_SUM, &
	296	comm2d, ierr )
[1320]	297	#else
	298	z0_av_global = z0_av_local
	299	#endif
[1314]	300
	301	z0_av_global = z0_av_global / ( ( ny + 1 ) * ( nx + 1 ) )
	302	!
	303	!-- Horizontal wind speed is zero below and at z0
[1359]	304	log_z_z0(0) = 0.0_wp
[1314]	305	!
	306	!-- Calculate vertical depth of the sublayers
[1322]	307	height_int = ( z_p - z0_av_global ) / REAL( number_of_sublayers, KIND=wp )
[1314]	308	!
	309	!-- Precalculate LOG(z/z0)
[1359]	310	height_p = 0.0_wp
[1314]	311	DO k = 1, number_of_sublayers
	312
	313	height_p = height_p + height_int
	314	log_z_z0(k) = LOG( height_p / z0_av_global )
	315
	316	ENDDO
	317
	318
	319	ENDIF
	320
	321	!
[1359]	322	!-- Check boundary condition and set internal variables
	323	SELECT CASE ( bc_par_b )
	324
	325	CASE ( 'absorb' )
	326	ibc_par_b = 1
	327
	328	CASE ( 'reflect' )
	329	ibc_par_b = 2
	330
	331	CASE DEFAULT
	332	WRITE( message_string, * ) 'unknown boundary condition ', &
	333	'bc_par_b = "', TRIM( bc_par_b ), '"'
	334	CALL message( 'lpm_init', 'PA0217', 1, 2, 0, 6, 0 )
	335
	336	END SELECT
	337	SELECT CASE ( bc_par_t )
	338
	339	CASE ( 'absorb' )
	340	ibc_par_t = 1
	341
	342	CASE ( 'reflect' )
	343	ibc_par_t = 2
	344
	345	CASE DEFAULT
	346	WRITE( message_string, * ) 'unknown boundary condition ', &
	347	'bc_par_t = "', TRIM( bc_par_t ), '"'
	348	CALL message( 'lpm_init', 'PA0218', 1, 2, 0, 6, 0 )
	349
	350	END SELECT
	351	SELECT CASE ( bc_par_lr )
	352
	353	CASE ( 'cyclic' )
	354	ibc_par_lr = 0
	355
	356	CASE ( 'absorb' )
	357	ibc_par_lr = 1
	358
	359	CASE ( 'reflect' )
	360	ibc_par_lr = 2
	361
	362	CASE DEFAULT
	363	WRITE( message_string, * ) 'unknown boundary condition ', &
	364	'bc_par_lr = "', TRIM( bc_par_lr ), '"'
	365	CALL message( 'lpm_init', 'PA0219', 1, 2, 0, 6, 0 )
	366
	367	END SELECT
	368	SELECT CASE ( bc_par_ns )
	369
	370	CASE ( 'cyclic' )
	371	ibc_par_ns = 0
	372
	373	CASE ( 'absorb' )
	374	ibc_par_ns = 1
	375
	376	CASE ( 'reflect' )
	377	ibc_par_ns = 2
	378
	379	CASE DEFAULT
	380	WRITE( message_string, * ) 'unknown boundary condition ', &
	381	'bc_par_ns = "', TRIM( bc_par_ns ), '"'
	382	CALL message( 'lpm_init', 'PA0220', 1, 2, 0, 6, 0 )
	383
	384	END SELECT
	385
	386	!
[828]	387	!-- Initialize collision kernels
	388	IF ( collision_kernel /= 'none' ) CALL init_kernels
	389
	390	!
[1]	391	!-- For the first model run of a possible job chain initialize the
[849]	392	!-- particles, otherwise read the particle data from restart file.
[1]	393	IF ( TRIM( initializing_actions ) == 'read_restart_data' &
	394	.AND. read_particles_from_restartfile ) THEN
	395
[849]	396	CALL lpm_read_restart_file
[1]	397
	398	ELSE
	399
	400	!
	401	!-- Allocate particle arrays and set attributes of the initial set of
	402	!-- particles, which can be also periodically released at later times.
	403	!-- Also allocate array for particle tail coordinates, if needed.
[1359]	404	ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	405	grid_particles(nzb+1:nzt,nys:nyn,nxl:nxr) )
[1]	406
[1359]	407	maximum_number_of_particles = 0
	408	number_of_particles = 0
[792]	409
	410	sort_count = 0
[1359]	411	prt_count = 0
[792]	412
[1]	413	!
	414	!-- Initialize all particles with dummy values (otherwise errors may
	415	!-- occur within restart runs). The reason for this is still not clear
	416	!-- and may be presumably caused by errors in the respective user-interface.
[1359]	417	#if defined( __twocachelines )
	418	zero_particle = particle_type( 0.0_wp, 0.0_sp, 0.0_sp, 0.0_sp, 0.0_sp, &
	419	0.0_sp, 0.0_sp, 0.0_wp, 0.0_wp, 0.0_wp, &
	420	0, .FALSE., 0.0_wp, 0.0_wp, 0.0_wp, &
	421	0.0_sp, 0.0_sp, 0.0_sp, 0.0_sp, 0.0_sp, &
	422	0.0_sp, 0, 0, 0, -1)
	423	#else
	424	zero_particle = particle_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
	425	0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
	426	0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
	427	0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0, 0, 0, &
	428	0, .FALSE., -1)
	429	#endif
	430	particle_groups = particle_groups_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
[1]	431
	432	!
	433	!-- Set the default particle size used for dvrp plots
[1359]	434	IF ( dvrp_psize == 9999999.9_wp ) dvrp_psize = 0.2_wp * dx
[1]	435
	436	!
	437	!-- Set values for the density ratio and radius for all particle
	438	!-- groups, if necessary
[1359]	439	IF ( density_ratio(1) == 9999999.9_wp ) density_ratio(1) = 0.0_wp
	440	IF ( radius(1) == 9999999.9_wp ) radius(1) = 0.0_wp
[1]	441	DO i = 2, number_of_particle_groups
[1359]	442	IF ( density_ratio(i) == 9999999.9_wp ) THEN
[1]	443	density_ratio(i) = density_ratio(i-1)
	444	ENDIF
[1359]	445	IF ( radius(i) == 9999999.9_wp ) radius(i) = radius(i-1)
[1]	446	ENDDO
	447
	448	DO i = 1, number_of_particle_groups
[1359]	449	IF ( density_ratio(i) /= 0.0_wp .AND. radius(i) == 0 ) THEN
[254]	450	WRITE( message_string, * ) 'particle group #', i, 'has a', &
	451	'density ratio /= 0 but radius = 0'
[849]	452	CALL message( 'lpm_init', 'PA0215', 1, 2, 0, 6, 0 )
[1]	453	ENDIF
	454	particle_groups(i)%density_ratio = density_ratio(i)
	455	particle_groups(i)%radius = radius(i)
	456	ENDDO
	457
[1359]	458	CALL lpm_create_particle (PHASE_INIT)
[1]	459	!
[1359]	460	!-- Set a seed value for the random number generator to be exclusively
	461	!-- used for the particle code. The generated random numbers should be
	462	!-- different on the different PEs.
	463	iran_part = iran_part + myid
	464
	465	!
	466	!-- User modification of initial particles
	467	CALL user_lpm_init
	468
	469	!
	470	!-- Open file for statistical informations about particle conditions
	471	IF ( write_particle_statistics ) THEN
	472	CALL check_open( 80 )
	473	WRITE ( 80, 8000 ) current_timestep_number, simulated_time, &
	474	number_of_particles, &
	475	maximum_number_of_particles
	476	CALL close_file( 80 )
	477	ENDIF
	478
	479	!
	480	!-- Check if particles are really uniform in color and radius (dvrp_size)
	481	!-- (uniform_particles is preset TRUE)
	482	IF ( uniform_particles ) THEN
	483	DO ip = nxl, nxr
	484	DO jp = nys, nyn
	485	DO kp = nzb+1, nzt
	486
	487	n = prt_count(kp,jp,ip)
	488	IF ( MINVAL( grid_particles(kp,jp,ip)%particles(1:n)%dvrp_psize ) == &
	489	MAXVAL( grid_particles(kp,jp,ip)%particles(1:n)%dvrp_psize ) .AND. &
	490	MINVAL( grid_particles(kp,jp,ip)%particles(1:n)%class ) == &
	491	MAXVAL( grid_particles(kp,jp,ip)%particles(1:n)%class ) ) THEN
	492	uniform_particles_l = .TRUE.
	493	ELSE
	494	uniform_particles_l = .FALSE.
	495	ENDIF
	496
	497	ENDDO
	498	ENDDO
	499	ENDDO
	500
	501	#if defined( __parallel )
	502	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	503	CALL MPI_ALLREDUCE( uniform_particles_l, uniform_particles, 1, &
	504	MPI_LOGICAL, MPI_LAND, comm2d, ierr )
	505	#else
	506	uniform_particles = uniform_particles_l
	507	#endif
	508
	509	ENDIF
	510
	511	!
	512	!-- Particles will probably become none-uniform, if their size and color
	513	!-- will be determined by flow variables
	514	IF ( particle_color /= 'none' .OR. particle_dvrpsize /= 'none' ) THEN
	515	uniform_particles = .FALSE.
	516	ENDIF
	517
	518	! !kk Not implemented aft individual particle array fort every gridcell
	519	! !
	520	! !-- Set the beginning of the particle tails and their age
	521	! IF ( use_particle_tails ) THEN
	522	! !
	523	! !-- Choose the maximum number of tails with respect to the maximum number
	524	! !-- of particles and skip_particles_for_tail
	525	! maximum_number_of_tails = maximum_number_of_particles / &
	526	! skip_particles_for_tail
	527	!
	528	! !
	529	! !-- Create a minimum number of tails in case that there is no tail
	530	! !-- initially (otherwise, index errors will occur when adressing the
	531	! !-- arrays below)
	532	! IF ( maximum_number_of_tails == 0 ) maximum_number_of_tails = 10
	533	!
	534	! ALLOCATE( particle_tail_coordinates(maximum_number_of_tailpoints,5, &
	535	! maximum_number_of_tails), &
	536	! new_tail_id(maximum_number_of_tails), &
	537	! tail_mask(maximum_number_of_tails) )
	538	!
	539	! particle_tail_coordinates = 0.0_wp
	540	! minimum_tailpoint_distance = minimum_tailpoint_distance**2
	541	! number_of_initial_tails = number_of_tails
	542	!
	543	! nn = 0
	544	! DO n = 1, number_of_particles
	545	! !
	546	! !-- Only for those particles marked above with a provisional tail_id
	547	! !-- tails will be created. Particles now get their final tail_id.
	548	! IF ( particles(n)%tail_id /= 0 ) THEN
	549	!
	550	! nn = nn + 1
	551	! particles(n)%tail_id = nn
	552	!
	553	! particle_tail_coordinates(1,1,nn) = particles(n)%x
	554	! particle_tail_coordinates(1,2,nn) = particles(n)%y
	555	! particle_tail_coordinates(1,3,nn) = particles(n)%z
	556	! particle_tail_coordinates(1,4,nn) = particles(n)%class
	557	! particles(n)%tailpoints = 1
	558	! IF ( minimum_tailpoint_distance /= 0.0_wp ) THEN
	559	! particle_tail_coordinates(2,1,nn) = particles(n)%x
	560	! particle_tail_coordinates(2,2,nn) = particles(n)%y
	561	! particle_tail_coordinates(2,3,nn) = particles(n)%z
	562	! particle_tail_coordinates(2,4,nn) = particles(n)%class
	563	! particle_tail_coordinates(1:2,5,nn) = 0.0_wp
	564	! particles(n)%tailpoints = 2
	565	! ENDIF
	566	!
	567	! ENDIF
	568	! ENDDO
	569	! ENDIF
	570	!
	571	! !
	572	! !-- Plot initial positions of particles (only if particle advection is
	573	! !-- switched on from the beginning of the simulation (t=0))
	574	! IF ( particle_advection_start == 0.0_wp ) CALL data_output_dvrp
	575
	576	ENDIF
	577
	578	!
	579	!-- To avoid programm abort, assign particles array to the local version of
	580	!-- first grid cell
	581	number_of_particles = prt_count(nzb+1,nys,nxl)
	582	particles => grid_particles(nzb+1,nys,nxl)%particles(1:number_of_particles)
	583
	584	!
	585	!-- Formats
	586	8000 FORMAT (I6,1X,F7.2,4X,I10,71X,I10)
	587
	588	END SUBROUTINE lpm_init
	589
	590	SUBROUTINE lpm_create_particle (phase)
	591
	592	USE lpm_exchange_horiz_mod, &
	593	ONLY: lpm_exchange_horiz, lpm_move_particle, realloc_particles_array
	594
	595	USE lpm_pack_arrays_mod, &
	596	ONLY: lpm_pack_all_arrays
	597
	598	IMPLICIT NONE
	599
	600	INTEGER(iwp) :: alloc_size !:
	601	INTEGER(iwp) :: i !:
	602	INTEGER(iwp) :: ip !:
	603	INTEGER(iwp) :: j !:
	604	INTEGER(iwp) :: jp !:
	605	INTEGER(iwp) :: kp !:
	606	INTEGER(iwp) :: loop_stride !:
	607	INTEGER(iwp) :: n !:
	608	INTEGER(iwp) :: new_size !:
	609	INTEGER(iwp) :: nn !:
	610
	611	INTEGER(iwp), INTENT(IN) :: phase !:
	612
	613	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_count !:
	614	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_start !:
	615
	616	LOGICAL :: first_stride !:
	617
	618	REAL(wp) :: pos_x !:
	619	REAL(wp) :: pos_y !:
	620	REAL(wp) :: pos_z !:
	621
	622	TYPE(particle_type),TARGET :: tmp_particle !:
	623
	624	!
	625	!-- Calculate particle positions and store particle attributes, if
	626	!-- particle is situated on this PE
	627	DO loop_stride = 1, 2
	628	first_stride = (loop_stride == 1)
	629	IF ( first_stride ) THEN
	630	local_count = 0 ! count number of particles
	631	ELSE
	632	local_count = prt_count ! Start address of new particles
	633	ENDIF
	634
[1]	635	n = 0
	636	DO i = 1, number_of_particle_groups
	637
	638	pos_z = psb(i)
	639
	640	DO WHILE ( pos_z <= pst(i) )
	641
	642	pos_y = pss(i)
	643
	644	DO WHILE ( pos_y <= psn(i) )
	645
[1359]	646	IF ( pos_y >= ( nys - 0.5_wp ) * dy .AND. &
	647	pos_y < ( nyn + 0.5_wp ) * dy ) THEN
[1]	648
	649	pos_x = psl(i)
	650
	651	DO WHILE ( pos_x <= psr(i) )
	652
[1359]	653	IF ( pos_x >= ( nxl - 0.5_wp ) * dx .AND. &
	654	pos_x < ( nxr + 0.5_wp ) * dx ) THEN
[1]	655
	656	DO j = 1, particles_per_point
	657
	658	n = n + 1
[1359]	659	#if defined( __twocachelines )
	660	tmp_particle%x = pos_x
	661	tmp_particle%y = pos_y
	662	tmp_particle%z = pos_z
	663	tmp_particle%age = 0.0_sp
	664	tmp_particle%age_m = 0.0_sp
	665	tmp_particle%dt_sum = 0.0_wp
	666	tmp_particle%dvrp_psize = dvrp_psize
	667	tmp_particle%e_m = 0.0_sp
	668	IF ( curvature_solution_effects ) THEN
	669	!
	670	!-- Initial values (internal timesteps, derivative)
	671	!-- for Rosenbrock method
	672	tmp_particle%rvar1 = 1.0E-12_wp
	673	tmp_particle%rvar2 = 1.0E-3_wp
	674	tmp_particle%rvar3 = -9999999.9_wp
	675	ELSE
	676	!
	677	!-- Initial values for SGS velocities
	678	tmp_particle%rvar1 = 0.0_wp
	679	tmp_particle%rvar2 = 0.0_wp
	680	tmp_particle%rvar3 = 0.0_wp
[1]	681	ENDIF
[1359]	682	tmp_particle%speed_x = 0.0_sp
	683	tmp_particle%speed_y = 0.0_sp
	684	tmp_particle%speed_z = 0.0_sp
	685	tmp_particle%origin_x = pos_x
	686	tmp_particle%origin_y = pos_y
	687	tmp_particle%origin_z = pos_z
	688	tmp_particle%radius = particle_groups(i)%radius
	689	tmp_particle%weight_factor = initial_weighting_factor
	690	tmp_particle%class = 1
	691	tmp_particle%group = i
	692	tmp_particle%tailpoints = 0
	693	tmp_particle%particle_mask = .TRUE.
	694	#else
	695	tmp_particle%x = pos_x
	696	tmp_particle%y = pos_y
	697	tmp_particle%z = pos_z
	698	tmp_particle%age = 0.0_wp
	699	tmp_particle%age_m = 0.0_wp
	700	tmp_particle%dt_sum = 0.0_wp
	701	tmp_particle%dvrp_psize = dvrp_psize
	702	tmp_particle%e_m = 0.0_wp
[824]	703	IF ( curvature_solution_effects ) THEN
	704	!
	705	!-- Initial values (internal timesteps, derivative)
	706	!-- for Rosenbrock method
[1359]	707	tmp_particle%rvar1 = 1.0E-12_wp
	708	tmp_particle%rvar2 = 1.0E-3_wp
	709	tmp_particle%rvar3 = -9999999.9_wp
[824]	710	ELSE
	711	!
	712	!-- Initial values for SGS velocities
[1359]	713	tmp_particle%rvar1 = 0.0_wp
	714	tmp_particle%rvar2 = 0.0_wp
	715	tmp_particle%rvar3 = 0.0_wp
[824]	716	ENDIF
[1359]	717	tmp_particle%speed_x = 0.0_wp
	718	tmp_particle%speed_y = 0.0_wp
	719	tmp_particle%speed_z = 0.0_wp
	720	tmp_particle%origin_x = pos_x
	721	tmp_particle%origin_y = pos_y
	722	tmp_particle%origin_z = pos_z
	723	tmp_particle%radius = particle_groups(i)%radius
	724	tmp_particle%weight_factor = initial_weighting_factor
	725	tmp_particle%class = 1
	726	tmp_particle%group = i
	727	tmp_particle%tailpoints = 0
	728	tmp_particle%particle_mask = .TRUE.
	729	#endif
[1]	730	IF ( use_particle_tails .AND. &
	731	MOD( n, skip_particles_for_tail ) == 0 ) THEN
	732	number_of_tails = number_of_tails + 1
	733	!
	734	!-- This is a temporary provisional setting (see
	735	!-- further below!)
[1359]	736	tmp_particle%tail_id = number_of_tails
[1]	737	ELSE
[1359]	738	tmp_particle%tail_id = 0
[1]	739	ENDIF
[1359]	740	ip = ( tmp_particle%x + 0.5_wp * dx ) * ddx
	741	jp = ( tmp_particle%y + 0.5_wp * dy ) * ddy
	742	kp = tmp_particle%z / dz + 1 + offset_ocean_nzt_m1
[1]	743
[1359]	744	local_count(kp,jp,ip) = local_count(kp,jp,ip) + 1
	745	IF ( .NOT. first_stride ) THEN
	746	IF ( ip < nxl .OR. jp < nys .OR. kp < nzb+1 ) THEN
	747	write(6,*) 'xl ',ip,jp,kp,nxl,nys,nzb+1
	748	ENDIF
	749	IF ( ip > nxr .OR. jp > nyn .OR. kp > nzt ) THEN
	750	write(6,*) 'xu ',ip,jp,kp,nxr,nyn,nzt
	751	ENDIF
	752	grid_particles(kp,jp,ip)%particles(local_count(kp,jp,ip)) = tmp_particle
	753	ENDIF
[1]	754	ENDDO
	755
	756	ENDIF
	757
	758	pos_x = pos_x + pdx(i)
	759
	760	ENDDO
	761
	762	ENDIF
	763
	764	pos_y = pos_y + pdy(i)
	765
	766	ENDDO
	767
	768	pos_z = pos_z + pdz(i)
	769
	770	ENDDO
	771
	772	ENDDO
	773
[1359]	774	IF ( first_stride ) THEN
	775	DO ip = nxl, nxr
	776	DO jp = nys, nyn
	777	DO kp = nzb+1, nzt
	778	IF ( phase == PHASE_INIT ) THEN
	779	IF ( local_count(kp,jp,ip) > 0 ) THEN
	780	alloc_size = MAX( INT( local_count(kp,jp,ip) * &
	781	( 1.0_wp + alloc_factor / 100.0_wp ) ), &
	782	min_nr_particle )
	783	ELSE
	784	alloc_size = min_nr_particle
	785	ENDIF
	786	ALLOCATE(grid_particles(kp,jp,ip)%particles(1:alloc_size))
	787	DO n = 1, alloc_size
	788	grid_particles(kp,jp,ip)%particles(n) = zero_particle
	789	ENDDO
	790	ELSEIF ( phase == PHASE_RELEASE ) THEN
	791	IF ( local_count(kp,jp,ip) > 0 ) THEN
	792	new_size = local_count(kp,jp,ip) + prt_count(kp,jp,ip)
	793	alloc_size = MAX( INT( new_size * ( 1.0_wp + &
	794	alloc_factor / 100.0_wp ) ), min_nr_particle )
	795	IF( alloc_size > SIZE( grid_particles(kp,jp,ip)%particles) ) THEN
	796	CALL realloc_particles_array(ip,jp,kp,alloc_size)
	797	ENDIF
	798	ENDIF
	799	ENDIF
	800	ENDDO
	801	ENDDO
	802	ENDDO
	803	ENDIF
	804	ENDDO
[1]	805
[1359]	806	local_start = prt_count+1
	807	prt_count = local_count
[1]	808	!
[1359]	809	!-- Add random fluctuation to particle positions
	810	IF ( random_start_position ) THEN
	811	DO ip = nxl, nxr
	812	DO jp = nys, nyn
	813	DO kp = nzb+1, nzt
	814	number_of_particles = prt_count(kp,jp,ip)
	815	IF ( number_of_particles <= 0 ) CYCLE
	816	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
[1]	817
[1359]	818	DO n = local_start(kp,jp,ip), number_of_particles !Move only new particles
	819	IF ( psl(particles(n)%group) /= psr(particles(n)%group) ) THEN
	820	particles(n)%x = particles(n)%x + &
	821	( random_function( iran_part ) - 0.5_wp ) * &
	822	pdx(particles(n)%group)
	823	ENDIF
	824	IF ( pss(particles(n)%group) /= psn(particles(n)%group) ) THEN
	825	particles(n)%y = particles(n)%y + &
	826	( random_function( iran_part ) - 0.5_wp ) * &
	827	pdy(particles(n)%group)
	828	ENDIF
	829	IF ( psb(particles(n)%group) /= pst(particles(n)%group) ) THEN
	830	particles(n)%z = particles(n)%z + &
	831	( random_function( iran_part ) - 0.5_wp ) * &
	832	pdz(particles(n)%group)
	833	ENDIF
	834	ENDDO
[1]	835	!
[1359]	836	!-- Identify particles located outside the model domain
	837	CALL lpm_boundary_conds( 'bottom/top' )
	838	ENDDO
	839	ENDDO
	840	ENDDO
[1]	841	!
[1359]	842	!-- Exchange particles between grid cells and processors
	843	CALL lpm_move_particle
	844	CALL lpm_exchange_horiz
[1]	845
[1359]	846	ENDIF
[1]	847	!
[1359]	848	!-- In case of random_start_position, delete particles identified by
	849	!-- lpm_exchange_horiz and lpm_boundary_conds. Then sort particles into blocks,
	850	!-- which is needed for a fast interpolation of the LES fields on the particle
	851	!-- position.
	852	CALL lpm_pack_all_arrays
[1]	853
	854	!
[1359]	855	!-- Determine maximum number of particles (i.e., all possible particles that
	856	!-- have been allocated) and the current number of particles
	857	DO ip = nxl, nxr
	858	DO jp = nys, nyn
	859	DO kp = nzb+1, nzt
	860	maximum_number_of_particles = maximum_number_of_particles &
	861	+ SIZE(grid_particles(kp,jp,ip)%particles)
	862	number_of_particles = number_of_particles &
	863	+ prt_count(kp,jp,ip)
[1]	864	ENDDO
[1359]	865	ENDDO
	866	ENDDO
[1]	867	!
[1359]	868	!-- Calculate the number of particles and tails of the total domain
[1]	869	#if defined( __parallel )
[1359]	870	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	871	CALL MPI_ALLREDUCE( number_of_particles, total_number_of_particles, 1, &
	872	MPI_INTEGER, MPI_SUM, comm2d, ierr )
	873	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	874	CALL MPI_ALLREDUCE( number_of_tails, total_number_of_tails, 1, &
	875	MPI_INTEGER, MPI_SUM, comm2d, ierr )
[1]	876	#else
[1359]	877	total_number_of_particles = number_of_particles
	878	total_number_of_tails = number_of_tails
[1]	879	#endif
	880
[1359]	881	RETURN
[1]	882
[1359]	883	END SUBROUTINE lpm_create_particle
[336]	884
[1359]	885	END MODULE lpm_init_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |