Home

Context Navigation

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

Last change on this file since 2156 was 2123, checked in by hoffmann, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 43.8 KB

Rev	Line
[1873]	1	!> @file lpm_init.f90
[2000]	2	!------------------------------------------------------------------------------!
[1036]	3	! This file is part of PALM.
	4	!
[2000]	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
	8	! version.
[1036]	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	!
[2101]	17	! Copyright 1997-2017 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[254]	20	! Current revisions:
[1]	21	! -----------------
[2123]	22	!
[1930]	23	!
	24	! Former revisions:
	25	! -----------------
	26	! $Id: lpm_init.f90 2123 2017-01-18 12:34:59Z hoffmann $
	27	!
[2123]	28	! 2122 2017-01-18 12:22:54Z hoffmann
	29	! Improved initialization of equilibrium aerosol radii
	30	! Calculation of particle ID
	31	!
[2001]	32	! 2000 2016-08-20 18:09:15Z knoop
	33	! Forced header and separation lines into 80 columns
	34	!
[1930]	35	! 2016-06-09 16:25:25Z suehring
[1929]	36	! Bugfix in determining initial particle height and grid index in case of
	37	! seed_follows_topography.
	38	! Bugfix concerning random positions, ensure that particles do not move more
	39	! than one grid length.
	40	! Bugfix logarithmic interpolation.
	41	! Initial setting of sgs_wf_part.
[1321]	42	!
[1891]	43	! 1890 2016-04-22 08:52:11Z hoffmann
	44	! Initialization of aerosol equilibrium radius not possible in supersaturated
	45	! environments. Therefore, a maximum supersaturation of -1 % is assumed during
	46	! initialization.
	47	!
[1874]	48	! 1873 2016-04-18 14:50:06Z maronga
[1929]	49	! Module renamed (removed _mod
[1874]	50	!
[1872]	51	! 1871 2016-04-15 11:46:09Z hoffmann
	52	! Initialization of aerosols added.
	53	!
[1851]	54	! 1850 2016-04-08 13:29:27Z maronga
	55	! Module renamed
	56	!
[1832]	57	! 1831 2016-04-07 13:15:51Z hoffmann
	58	! curvature_solution_effects moved to particle_attributes
	59	!
[1823]	60	! 1822 2016-04-07 07:49:42Z hoffmann
	61	! Unused variables removed.
	62	!
[1784]	63	! 1783 2016-03-06 18:36:17Z raasch
	64	! netcdf module added
	65	!
[1726]	66	! 1725 2015-11-17 13:01:51Z hoffmann
	67	! Bugfix: Processor-dependent seed for random function is generated before it is
	68	! used.
	69	!
[1692]	70	! 1691 2015-10-26 16:17:44Z maronga
	71	! Renamed prandtl_layer to constant_flux_layer.
	72	!
[1686]	73	! 1685 2015-10-08 07:32:13Z raasch
	74	! bugfix concerning vertical index offset in case of ocean
	75	!
[1683]	76	! 1682 2015-10-07 23:56:08Z knoop
	77	! Code annotations made doxygen readable
	78	!
[1576]	79	! 1575 2015-03-27 09:56:27Z raasch
	80	! initial vertical particle position is allowed to follow the topography
	81	!
[1360]	82	! 1359 2014-04-11 17:15:14Z hoffmann
	83	! New particle structure integrated.
	84	! Kind definition added to all floating point numbers.
	85	! lpm_init changed form a subroutine to a module.
	86	!
[1329]	87	! 1327 2014-03-21 11:00:16Z raasch
	88	! -netcdf_output
	89	!
[1323]	90	! 1322 2014-03-20 16:38:49Z raasch
	91	! REAL functions provided with KIND-attribute
	92	!
[1321]	93	! 1320 2014-03-20 08:40:49Z raasch
[1320]	94	! ONLY-attribute added to USE-statements,
	95	! kind-parameters added to all INTEGER and REAL declaration statements,
	96	! kinds are defined in new module kinds,
	97	! revision history before 2012 removed,
	98	! comment fields (!:) to be used for variable explanations added to
	99	! all variable declaration statements
	100	! bugfix: #if defined( __parallel ) added
[850]	101	!
[1315]	102	! 1314 2014-03-14 18:25:17Z suehring
	103	! Vertical logarithmic interpolation of horizontal particle speed for particles
	104	! between roughness height and first vertical grid level.
	105	!
[1093]	106	! 1092 2013-02-02 11:24:22Z raasch
	107	! unused variables removed
	108	!
[1037]	109	! 1036 2012-10-22 13:43:42Z raasch
	110	! code put under GPL (PALM 3.9)
	111	!
[850]	112	! 849 2012-03-15 10:35:09Z raasch
[849]	113	! routine renamed: init_particles -> lpm_init
	114	! de_dx, de_dy, de_dz are allocated here (instead of automatic arrays in
	115	! advec_particles),
	116	! sort_particles renamed lpm_sort_arrays, user_init_particles renamed lpm_init
[392]	117	!
[829]	118	! 828 2012-02-21 12:00:36Z raasch
	119	! call of init_kernels, particle feature color renamed class
	120	!
[826]	121	! 824 2012-02-17 09:09:57Z raasch
	122	! particle attributes speed_x\|y\|z_sgs renamed rvar1\|2\|3,
	123	! array particles implemented as pointer
	124	!
[668]	125	! 667 2010-12-23 12:06:00Z suehring/gryschka
	126	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation
	127	! of arrays.
	128	!
[1]	129	! Revision 1.1 1999/11/25 16:22:38 raasch
	130	! Initial revision
	131	!
	132	!
	133	! Description:
	134	! ------------
[1682]	135	!> This routine initializes a set of particles and their attributes (position,
	136	!> radius, ..) which are used by the Lagrangian particle model (see lpm).
[1]	137	!------------------------------------------------------------------------------!
[1682]	138	MODULE lpm_init_mod
	139
[1]	140
[1320]	141	USE arrays_3d, &
	142	ONLY: de_dx, de_dy, de_dz, zu, zw, z0
	143
	144	USE control_parameters, &
[1691]	145	ONLY: cloud_droplets, constant_flux_layer, current_timestep_number, &
[1783]	146	dz, initializing_actions, message_string, ocean, simulated_time
[1320]	147
	148	USE grid_variables, &
[1359]	149	ONLY: ddx, dx, ddy, dy
[1320]	150
	151	USE indices, &
[1575]	152	ONLY: nx, nxl, nxlg, nxrg, nxr, ny, nyn, nys, nyng, nysg, nz, nzb, &
	153	nzb_w_inner, nzt
[1320]	154
	155	USE kinds
	156
	157	USE lpm_collision_kernels_mod, &
	158	ONLY: init_kernels
	159
[1783]	160	USE netcdf_interface, &
	161	ONLY: netcdf_data_format
	162
[1320]	163	USE particle_attributes, &
[1359]	164	ONLY: alloc_factor, bc_par_b, bc_par_lr, bc_par_ns, bc_par_t, &
	165	block_offset, block_offset_def, collision_kernel, &
[1831]	166	curvature_solution_effects, &
[1822]	167	density_ratio, grid_particles, &
[1359]	168	initial_weighting_factor, ibc_par_b, ibc_par_lr, ibc_par_ns, &
	169	ibc_par_t, iran_part, log_z_z0, &
	170	max_number_of_particle_groups, maximum_number_of_particles, &
[1822]	171	min_nr_particle, mpi_particle_type, &
	172	number_of_particles, &
[1320]	173	number_of_particle_groups, number_of_sublayers, &
[1822]	174	offset_ocean_nzt, offset_ocean_nzt_m1, &
[1359]	175	particles, particle_advection_start, particle_groups, &
	176	particle_groups_type, particles_per_point, &
[1822]	177	particle_type, pdx, pdy, pdz, &
[1359]	178	prt_count, psb, psl, psn, psr, pss, pst, &
[1320]	179	radius, random_start_position, read_particles_from_restartfile,&
[1929]	180	seed_follows_topography, sgs_wf_part, sort_count, &
[1822]	181	total_number_of_particles, &
	182	use_sgs_for_particles, &
[1359]	183	write_particle_statistics, uniform_particles, zero_particle, &
	184	z0_av_global
[1320]	185
[1]	186	USE pegrid
	187
[1320]	188	USE random_function_mod, &
	189	ONLY: random_function
[1]	190
[1359]	191	IMPLICIT NONE
[1320]	192
[1359]	193	PRIVATE
	194
[1682]	195	INTEGER(iwp), PARAMETER :: PHASE_INIT = 1 !<
	196	INTEGER(iwp), PARAMETER, PUBLIC :: PHASE_RELEASE = 2 !<
[1359]	197
	198	INTERFACE lpm_init
	199	MODULE PROCEDURE lpm_init
	200	END INTERFACE lpm_init
	201
	202	INTERFACE lpm_create_particle
	203	MODULE PROCEDURE lpm_create_particle
	204	END INTERFACE lpm_create_particle
	205
	206	PUBLIC lpm_init, lpm_create_particle
	207
[1929]	208	CONTAINS
[1359]	209
[1682]	210	!------------------------------------------------------------------------------!
	211	! Description:
	212	! ------------
	213	!> @todo Missing subroutine description.
	214	!------------------------------------------------------------------------------!
[1359]	215	SUBROUTINE lpm_init
	216
	217	USE lpm_collision_kernels_mod, &
	218	ONLY: init_kernels
	219
[1]	220	IMPLICIT NONE
	221
[1682]	222	INTEGER(iwp) :: i !<
	223	INTEGER(iwp) :: j !<
	224	INTEGER(iwp) :: k !<
[1320]	225
[1]	226	#if defined( __parallel )
[1682]	227	INTEGER(iwp), DIMENSION(3) :: blocklengths !<
	228	INTEGER(iwp), DIMENSION(3) :: displacements !<
	229	INTEGER(iwp), DIMENSION(3) :: types !<
[1]	230	#endif
	231
[1682]	232	REAL(wp) :: height_int !<
	233	REAL(wp) :: height_p !<
	234	REAL(wp) :: z_p !<
	235	REAL(wp) :: z0_av_local !<
[1]	236
	237	#if defined( __parallel )
	238	!
	239	!-- Define MPI derived datatype for FORTRAN datatype particle_type (see module
[82]	240	!-- particle_attributes). Integer length is 4 byte, Real is 8 byte
[1359]	241	blocklengths(1) = 19; blocklengths(2) = 6; blocklengths(3) = 1
	242	displacements(1) = 0; displacements(2) = 152; displacements(3) = 176
	243
[1]	244	types(1) = MPI_REAL
	245	types(2) = MPI_INTEGER
	246	types(3) = MPI_UB
	247	CALL MPI_TYPE_STRUCT( 3, blocklengths, displacements, types, &
	248	mpi_particle_type, ierr )
	249	CALL MPI_TYPE_COMMIT( mpi_particle_type, ierr )
	250	#endif
	251
	252	!
[150]	253	!-- In case of oceans runs, the vertical index calculations need an offset,
	254	!-- because otherwise the k indices will become negative
	255	IF ( ocean ) THEN
	256	offset_ocean_nzt = nzt
	257	offset_ocean_nzt_m1 = nzt - 1
	258	ENDIF
	259
[1359]	260	!
	261	!-- Define block offsets for dividing a gridcell in 8 sub cells
[150]	262
[1359]	263	block_offset(0) = block_offset_def (-1,-1,-1)
	264	block_offset(1) = block_offset_def (-1,-1, 0)
	265	block_offset(2) = block_offset_def (-1, 0,-1)
	266	block_offset(3) = block_offset_def (-1, 0, 0)
	267	block_offset(4) = block_offset_def ( 0,-1,-1)
	268	block_offset(5) = block_offset_def ( 0,-1, 0)
	269	block_offset(6) = block_offset_def ( 0, 0,-1)
	270	block_offset(7) = block_offset_def ( 0, 0, 0)
[150]	271	!
[1]	272	!-- Check the number of particle groups.
	273	IF ( number_of_particle_groups > max_number_of_particle_groups ) THEN
[274]	274	WRITE( message_string, * ) 'max_number_of_particle_groups =', &
	275	max_number_of_particle_groups , &
[254]	276	'&number_of_particle_groups reset to ', &
	277	max_number_of_particle_groups
[849]	278	CALL message( 'lpm_init', 'PA0213', 0, 1, 0, 6, 0 )
[1]	279	number_of_particle_groups = max_number_of_particle_groups
	280	ENDIF
	281
	282	!
	283	!-- Set default start positions, if necessary
[1359]	284	IF ( psl(1) == 9999999.9_wp ) psl(1) = -0.5_wp * dx
	285	IF ( psr(1) == 9999999.9_wp ) psr(1) = ( nx + 0.5_wp ) * dx
	286	IF ( pss(1) == 9999999.9_wp ) pss(1) = -0.5_wp * dy
	287	IF ( psn(1) == 9999999.9_wp ) psn(1) = ( ny + 0.5_wp ) * dy
	288	IF ( psb(1) == 9999999.9_wp ) psb(1) = zu(nz/2)
	289	IF ( pst(1) == 9999999.9_wp ) pst(1) = psb(1)
[1]	290
[1359]	291	IF ( pdx(1) == 9999999.9_wp .OR. pdx(1) == 0.0_wp ) pdx(1) = dx
	292	IF ( pdy(1) == 9999999.9_wp .OR. pdy(1) == 0.0_wp ) pdy(1) = dy
	293	IF ( pdz(1) == 9999999.9_wp .OR. pdz(1) == 0.0_wp ) pdz(1) = zu(2) - zu(1)
[1]	294
	295	DO j = 2, number_of_particle_groups
[1359]	296	IF ( psl(j) == 9999999.9_wp ) psl(j) = psl(j-1)
	297	IF ( psr(j) == 9999999.9_wp ) psr(j) = psr(j-1)
	298	IF ( pss(j) == 9999999.9_wp ) pss(j) = pss(j-1)
	299	IF ( psn(j) == 9999999.9_wp ) psn(j) = psn(j-1)
	300	IF ( psb(j) == 9999999.9_wp ) psb(j) = psb(j-1)
	301	IF ( pst(j) == 9999999.9_wp ) pst(j) = pst(j-1)
	302	IF ( pdx(j) == 9999999.9_wp .OR. pdx(j) == 0.0_wp ) pdx(j) = pdx(j-1)
	303	IF ( pdy(j) == 9999999.9_wp .OR. pdy(j) == 0.0_wp ) pdy(j) = pdy(j-1)
	304	IF ( pdz(j) == 9999999.9_wp .OR. pdz(j) == 0.0_wp ) pdz(j) = pdz(j-1)
[1]	305	ENDDO
	306
	307	!
[1929]	308	!-- Allocate arrays required for calculating particle SGS velocities.
	309	!-- Initialize prefactor required for stoachastic Weil equation.
[1822]	310	IF ( use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN
[849]	311	ALLOCATE( de_dx(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	312	de_dy(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	313	de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1929]	314
	315	sgs_wf_part = 1.0_wp / 3.0_wp
[849]	316	ENDIF
	317
	318	!
[1314]	319	!-- Allocate array required for logarithmic vertical interpolation of
	320	!-- horizontal particle velocities between the surface and the first vertical
	321	!-- grid level. In order to avoid repeated CPU cost-intensive CALLS of
	322	!-- intrinsic FORTRAN procedure LOG(z/z0), LOG(z/z0) is precalculated for
	323	!-- several heights. Splitting into 20 sublayers turned out to be sufficient.
	324	!-- To obtain exact height levels of particles, linear interpolation is applied
	325	!-- (see lpm_advec.f90).
[1691]	326	IF ( constant_flux_layer ) THEN
[1314]	327
	328	ALLOCATE ( log_z_z0(0:number_of_sublayers) )
	329	z_p = zu(nzb+1) - zw(nzb)
	330
	331	!
	332	!-- Calculate horizontal mean value of z0 used for logartihmic
	333	!-- interpolation. Note: this is not exact for heterogeneous z0.
	334	!-- However, sensitivity studies showed that the effect is
	335	!-- negligible.
	336	z0_av_local = SUM( z0(nys:nyn,nxl:nxr) )
[1359]	337	z0_av_global = 0.0_wp
[1314]	338
[1320]	339	#if defined( __parallel )
[1314]	340	CALL MPI_ALLREDUCE(z0_av_local, z0_av_global, 1, MPI_REAL, MPI_SUM, &
	341	comm2d, ierr )
[1320]	342	#else
	343	z0_av_global = z0_av_local
	344	#endif
[1314]	345
	346	z0_av_global = z0_av_global / ( ( ny + 1 ) * ( nx + 1 ) )
	347	!
	348	!-- Horizontal wind speed is zero below and at z0
[1359]	349	log_z_z0(0) = 0.0_wp
[1314]	350	!
	351	!-- Calculate vertical depth of the sublayers
[1322]	352	height_int = ( z_p - z0_av_global ) / REAL( number_of_sublayers, KIND=wp )
[1314]	353	!
	354	!-- Precalculate LOG(z/z0)
[1929]	355	height_p = z0_av_global
[1314]	356	DO k = 1, number_of_sublayers
	357
	358	height_p = height_p + height_int
	359	log_z_z0(k) = LOG( height_p / z0_av_global )
	360
	361	ENDDO
	362
	363	ENDIF
	364
	365	!
[1359]	366	!-- Check boundary condition and set internal variables
	367	SELECT CASE ( bc_par_b )
	368
	369	CASE ( 'absorb' )
	370	ibc_par_b = 1
	371
	372	CASE ( 'reflect' )
	373	ibc_par_b = 2
	374
	375	CASE DEFAULT
	376	WRITE( message_string, * ) 'unknown boundary condition ', &
	377	'bc_par_b = "', TRIM( bc_par_b ), '"'
	378	CALL message( 'lpm_init', 'PA0217', 1, 2, 0, 6, 0 )
	379
	380	END SELECT
	381	SELECT CASE ( bc_par_t )
	382
	383	CASE ( 'absorb' )
	384	ibc_par_t = 1
	385
	386	CASE ( 'reflect' )
	387	ibc_par_t = 2
	388
	389	CASE DEFAULT
	390	WRITE( message_string, * ) 'unknown boundary condition ', &
	391	'bc_par_t = "', TRIM( bc_par_t ), '"'
	392	CALL message( 'lpm_init', 'PA0218', 1, 2, 0, 6, 0 )
	393
	394	END SELECT
	395	SELECT CASE ( bc_par_lr )
	396
	397	CASE ( 'cyclic' )
	398	ibc_par_lr = 0
	399
	400	CASE ( 'absorb' )
	401	ibc_par_lr = 1
	402
	403	CASE ( 'reflect' )
	404	ibc_par_lr = 2
	405
	406	CASE DEFAULT
	407	WRITE( message_string, * ) 'unknown boundary condition ', &
	408	'bc_par_lr = "', TRIM( bc_par_lr ), '"'
	409	CALL message( 'lpm_init', 'PA0219', 1, 2, 0, 6, 0 )
	410
	411	END SELECT
	412	SELECT CASE ( bc_par_ns )
	413
	414	CASE ( 'cyclic' )
	415	ibc_par_ns = 0
	416
	417	CASE ( 'absorb' )
	418	ibc_par_ns = 1
	419
	420	CASE ( 'reflect' )
	421	ibc_par_ns = 2
	422
	423	CASE DEFAULT
	424	WRITE( message_string, * ) 'unknown boundary condition ', &
	425	'bc_par_ns = "', TRIM( bc_par_ns ), '"'
	426	CALL message( 'lpm_init', 'PA0220', 1, 2, 0, 6, 0 )
	427
	428	END SELECT
	429
	430	!
[828]	431	!-- Initialize collision kernels
	432	IF ( collision_kernel /= 'none' ) CALL init_kernels
	433
	434	!
[1]	435	!-- For the first model run of a possible job chain initialize the
[849]	436	!-- particles, otherwise read the particle data from restart file.
[1]	437	IF ( TRIM( initializing_actions ) == 'read_restart_data' &
	438	.AND. read_particles_from_restartfile ) THEN
	439
[849]	440	CALL lpm_read_restart_file
[1]	441
	442	ELSE
	443
	444	!
	445	!-- Allocate particle arrays and set attributes of the initial set of
	446	!-- particles, which can be also periodically released at later times.
[1359]	447	ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	448	grid_particles(nzb+1:nzt,nys:nyn,nxl:nxr) )
[1]	449
[1359]	450	maximum_number_of_particles = 0
	451	number_of_particles = 0
[792]	452
	453	sort_count = 0
[1359]	454	prt_count = 0
[792]	455
[1]	456	!
[2122]	457	!-- initialize counter for particle IDs
	458	grid_particles%id_counter = 0
	459
	460	!
[1]	461	!-- Initialize all particles with dummy values (otherwise errors may
	462	!-- occur within restart runs). The reason for this is still not clear
	463	!-- and may be presumably caused by errors in the respective user-interface.
[1359]	464	zero_particle = particle_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
	465	0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
	466	0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
[2122]	467	0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
	468	0, 0, 0, 0, .FALSE., -1 )
[1822]	469
[1359]	470	particle_groups = particle_groups_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
[1]	471
	472	!
	473	!-- Set values for the density ratio and radius for all particle
	474	!-- groups, if necessary
[1359]	475	IF ( density_ratio(1) == 9999999.9_wp ) density_ratio(1) = 0.0_wp
	476	IF ( radius(1) == 9999999.9_wp ) radius(1) = 0.0_wp
[1]	477	DO i = 2, number_of_particle_groups
[1359]	478	IF ( density_ratio(i) == 9999999.9_wp ) THEN
[1]	479	density_ratio(i) = density_ratio(i-1)
	480	ENDIF
[1359]	481	IF ( radius(i) == 9999999.9_wp ) radius(i) = radius(i-1)
[1]	482	ENDDO
	483
	484	DO i = 1, number_of_particle_groups
[1359]	485	IF ( density_ratio(i) /= 0.0_wp .AND. radius(i) == 0 ) THEN
[254]	486	WRITE( message_string, * ) 'particle group #', i, 'has a', &
	487	'density ratio /= 0 but radius = 0'
[849]	488	CALL message( 'lpm_init', 'PA0215', 1, 2, 0, 6, 0 )
[1]	489	ENDIF
	490	particle_groups(i)%density_ratio = density_ratio(i)
	491	particle_groups(i)%radius = radius(i)
	492	ENDDO
	493
	494	!
[1359]	495	!-- Set a seed value for the random number generator to be exclusively
	496	!-- used for the particle code. The generated random numbers should be
	497	!-- different on the different PEs.
	498	iran_part = iran_part + myid
	499
[1725]	500	CALL lpm_create_particle (PHASE_INIT)
[1359]	501	!
	502	!-- User modification of initial particles
	503	CALL user_lpm_init
	504
	505	!
	506	!-- Open file for statistical informations about particle conditions
	507	IF ( write_particle_statistics ) THEN
	508	CALL check_open( 80 )
	509	WRITE ( 80, 8000 ) current_timestep_number, simulated_time, &
	510	number_of_particles, &
	511	maximum_number_of_particles
	512	CALL close_file( 80 )
	513	ENDIF
	514
	515	ENDIF
	516
	517	!
	518	!-- To avoid programm abort, assign particles array to the local version of
	519	!-- first grid cell
	520	number_of_particles = prt_count(nzb+1,nys,nxl)
	521	particles => grid_particles(nzb+1,nys,nxl)%particles(1:number_of_particles)
	522	!
	523	!-- Formats
	524	8000 FORMAT (I6,1X,F7.2,4X,I10,71X,I10)
	525
	526	END SUBROUTINE lpm_init
	527
[1682]	528	!------------------------------------------------------------------------------!
	529	! Description:
	530	! ------------
	531	!> @todo Missing subroutine description.
	532	!------------------------------------------------------------------------------!
[1359]	533	SUBROUTINE lpm_create_particle (phase)
	534
	535	USE lpm_exchange_horiz_mod, &
	536	ONLY: lpm_exchange_horiz, lpm_move_particle, realloc_particles_array
	537
	538	USE lpm_pack_arrays_mod, &
	539	ONLY: lpm_pack_all_arrays
	540
[1871]	541	USE particle_attributes, &
[1929]	542	ONLY: deleted_particles, monodisperse_aerosols
[1871]	543
[1359]	544	IMPLICIT NONE
	545
[1929]	546	INTEGER(iwp) :: alloc_size !< relative increase of allocated memory for particles
	547	INTEGER(iwp) :: i !< loop variable ( particle groups )
	548	INTEGER(iwp) :: ip !< index variable along x
	549	INTEGER(iwp) :: j !< loop variable ( particles per point )
	550	INTEGER(iwp) :: jp !< index variable along y
	551	INTEGER(iwp) :: kp !< index variable along z
	552	INTEGER(iwp) :: loop_stride !< loop variable for initialization
	553	INTEGER(iwp) :: n !< loop variable ( number of particles )
	554	INTEGER(iwp) :: new_size !< new size of allocated memory for particles
[1359]	555
[1929]	556	INTEGER(iwp), INTENT(IN) :: phase !< mode of inititialization
[1359]	557
[1929]	558	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_count !< start address of new particle
	559	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_start !< start address of new particle
[1359]	560
[1929]	561	LOGICAL :: first_stride !< flag for initialization
[1359]	562
[1929]	563	REAL(wp) :: pos_x !< increment for particle position in x
	564	REAL(wp) :: pos_y !< increment for particle position in y
	565	REAL(wp) :: pos_z !< increment for particle position in z
	566	REAL(wp) :: rand_contr !< dummy argument for random position
[1359]	567
[1929]	568	TYPE(particle_type),TARGET :: tmp_particle !< temporary particle used for initialization
[1359]	569
	570	!
	571	!-- Calculate particle positions and store particle attributes, if
	572	!-- particle is situated on this PE
	573	DO loop_stride = 1, 2
	574	first_stride = (loop_stride == 1)
	575	IF ( first_stride ) THEN
	576	local_count = 0 ! count number of particles
	577	ELSE
	578	local_count = prt_count ! Start address of new particles
	579	ENDIF
	580
[1]	581	n = 0
	582	DO i = 1, number_of_particle_groups
	583
	584	pos_z = psb(i)
	585
	586	DO WHILE ( pos_z <= pst(i) )
	587
	588	pos_y = pss(i)
	589
	590	DO WHILE ( pos_y <= psn(i) )
	591
[1359]	592	IF ( pos_y >= ( nys - 0.5_wp ) * dy .AND. &
	593	pos_y < ( nyn + 0.5_wp ) * dy ) THEN
[1]	594
	595	pos_x = psl(i)
	596
[1575]	597	xloop: DO WHILE ( pos_x <= psr(i) )
[1]	598
[1359]	599	IF ( pos_x >= ( nxl - 0.5_wp ) * dx .AND. &
	600	pos_x < ( nxr + 0.5_wp ) * dx ) THEN
[1]	601
	602	DO j = 1, particles_per_point
	603
	604	n = n + 1
[1359]	605	tmp_particle%x = pos_x
	606	tmp_particle%y = pos_y
	607	tmp_particle%z = pos_z
	608	tmp_particle%age = 0.0_wp
	609	tmp_particle%age_m = 0.0_wp
	610	tmp_particle%dt_sum = 0.0_wp
[2122]	611	tmp_particle%user = 0.0_wp !unused, free for the user
[1359]	612	tmp_particle%e_m = 0.0_wp
[824]	613	IF ( curvature_solution_effects ) THEN
	614	!
	615	!-- Initial values (internal timesteps, derivative)
	616	!-- for Rosenbrock method
[1871]	617	tmp_particle%rvar1 = 1.0E-6_wp !last Rosenbrock timestep
	618	tmp_particle%rvar2 = 0.1E-6_wp !dry aerosol radius
[2122]	619	tmp_particle%rvar3 = -9999999.9_wp !unused in this configuration
[824]	620	ELSE
	621	!
	622	!-- Initial values for SGS velocities
[1359]	623	tmp_particle%rvar1 = 0.0_wp
	624	tmp_particle%rvar2 = 0.0_wp
	625	tmp_particle%rvar3 = 0.0_wp
[824]	626	ENDIF
[1359]	627	tmp_particle%speed_x = 0.0_wp
	628	tmp_particle%speed_y = 0.0_wp
	629	tmp_particle%speed_z = 0.0_wp
	630	tmp_particle%origin_x = pos_x
	631	tmp_particle%origin_y = pos_y
	632	tmp_particle%origin_z = pos_z
	633	tmp_particle%radius = particle_groups(i)%radius
	634	tmp_particle%weight_factor = initial_weighting_factor
	635	tmp_particle%class = 1
	636	tmp_particle%group = i
[2122]	637	tmp_particle%id1 = 0
	638	tmp_particle%id2 = 0
[1359]	639	tmp_particle%particle_mask = .TRUE.
[2122]	640	tmp_particle%block_nr = -1
[1]	641	!
[1575]	642	!-- Determine the grid indices of the particle position
[1359]	643	ip = ( tmp_particle%x + 0.5_wp * dx ) * ddx
	644	jp = ( tmp_particle%y + 0.5_wp * dy ) * ddy
[1685]	645	kp = tmp_particle%z / dz + 1 + offset_ocean_nzt
[1]	646
[1575]	647	IF ( seed_follows_topography ) THEN
	648	!
	649	!-- Particle height is given relative to topography
	650	kp = kp + nzb_w_inner(jp,ip)
[1929]	651	tmp_particle%z = tmp_particle%z + &
	652	zw(nzb_w_inner(jp,ip))
[1575]	653	IF ( kp > nzt ) THEN
	654	pos_x = pos_x + pdx(i)
	655	CYCLE xloop
	656	ENDIF
[1929]	657	ELSEIF ( .NOT. seed_follows_topography .AND. &
	658	tmp_particle%z <= zw(nzb_w_inner(jp,ip)) ) THEN
	659	pos_x = pos_x + pdx(i)
	660	CYCLE xloop
[1575]	661	ENDIF
	662
[1359]	663	local_count(kp,jp,ip) = local_count(kp,jp,ip) + 1
	664	IF ( .NOT. first_stride ) THEN
	665	IF ( ip < nxl .OR. jp < nys .OR. kp < nzb+1 ) THEN
	666	write(6,*) 'xl ',ip,jp,kp,nxl,nys,nzb+1
	667	ENDIF
	668	IF ( ip > nxr .OR. jp > nyn .OR. kp > nzt ) THEN
	669	write(6,*) 'xu ',ip,jp,kp,nxr,nyn,nzt
	670	ENDIF
	671	grid_particles(kp,jp,ip)%particles(local_count(kp,jp,ip)) = tmp_particle
[1929]	672
[1359]	673	ENDIF
[1]	674	ENDDO
	675
	676	ENDIF
	677
	678	pos_x = pos_x + pdx(i)
	679
[1575]	680	ENDDO xloop
[1]	681
	682	ENDIF
	683
	684	pos_y = pos_y + pdy(i)
	685
	686	ENDDO
	687
	688	pos_z = pos_z + pdz(i)
	689
	690	ENDDO
	691
	692	ENDDO
	693
[1359]	694	IF ( first_stride ) THEN
	695	DO ip = nxl, nxr
	696	DO jp = nys, nyn
	697	DO kp = nzb+1, nzt
	698	IF ( phase == PHASE_INIT ) THEN
	699	IF ( local_count(kp,jp,ip) > 0 ) THEN
	700	alloc_size = MAX( INT( local_count(kp,jp,ip) * &
	701	( 1.0_wp + alloc_factor / 100.0_wp ) ), &
	702	min_nr_particle )
	703	ELSE
	704	alloc_size = min_nr_particle
	705	ENDIF
	706	ALLOCATE(grid_particles(kp,jp,ip)%particles(1:alloc_size))
	707	DO n = 1, alloc_size
	708	grid_particles(kp,jp,ip)%particles(n) = zero_particle
	709	ENDDO
	710	ELSEIF ( phase == PHASE_RELEASE ) THEN
	711	IF ( local_count(kp,jp,ip) > 0 ) THEN
	712	new_size = local_count(kp,jp,ip) + prt_count(kp,jp,ip)
	713	alloc_size = MAX( INT( new_size * ( 1.0_wp + &
	714	alloc_factor / 100.0_wp ) ), min_nr_particle )
	715	IF( alloc_size > SIZE( grid_particles(kp,jp,ip)%particles) ) THEN
	716	CALL realloc_particles_array(ip,jp,kp,alloc_size)
	717	ENDIF
	718	ENDIF
	719	ENDIF
[1929]	720
[1359]	721	ENDDO
	722	ENDDO
	723	ENDDO
	724	ENDIF
[1929]	725
[1359]	726	ENDDO
[1]	727
[1359]	728	local_start = prt_count+1
	729	prt_count = local_count
[1871]	730
[1]	731	!
[2122]	732	!-- Calculate particle IDs
	733	DO ip = nxl, nxr
	734	DO jp = nys, nyn
	735	DO kp = nzb+1, nzt
	736	number_of_particles = prt_count(kp,jp,ip)
	737	IF ( number_of_particles <= 0 ) CYCLE
	738	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
	739
	740	DO n = local_start(kp,jp,ip), number_of_particles !only new particles
	741
	742	particles(n)%id1 = 10000 * grid_particles(kp,jp,ip)%id_counter + kp
	743	particles(n)%id2 = 10000 * jp + ip
	744
	745	grid_particles(kp,jp,ip)%id_counter = &
	746	grid_particles(kp,jp,ip)%id_counter + 1
	747
	748	ENDDO
	749
	750	ENDDO
	751	ENDDO
	752	ENDDO
	753
	754	!
[1871]	755	!-- Initialize aerosol background spectrum
	756	IF ( curvature_solution_effects .AND. .NOT. monodisperse_aerosols ) THEN
	757	CALL lpm_init_aerosols(local_start)
	758	ENDIF
	759
	760	!
[1929]	761	!-- Add random fluctuation to particle positions.
[1359]	762	IF ( random_start_position ) THEN
	763	DO ip = nxl, nxr
	764	DO jp = nys, nyn
	765	DO kp = nzb+1, nzt
	766	number_of_particles = prt_count(kp,jp,ip)
	767	IF ( number_of_particles <= 0 ) CYCLE
	768	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
[1929]	769	!
	770	!-- Move only new particles. Moreover, limit random fluctuation
	771	!-- in order to prevent that particles move more than one grid box,
	772	!-- which would lead to problems concerning particle exchange
	773	!-- between processors in case pdx/pdy are larger than dx/dy,
	774	!-- respectively.
	775	DO n = local_start(kp,jp,ip), number_of_particles
[1359]	776	IF ( psl(particles(n)%group) /= psr(particles(n)%group) ) THEN
[1929]	777	rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
	778	pdx(particles(n)%group)
[1359]	779	particles(n)%x = particles(n)%x + &
[1929]	780	MERGE( rand_contr, SIGN( dx, rand_contr ), &
	781	ABS( rand_contr ) < dx &
	782	)
[1359]	783	ENDIF
	784	IF ( pss(particles(n)%group) /= psn(particles(n)%group) ) THEN
[1929]	785	rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
	786	pdy(particles(n)%group)
[1359]	787	particles(n)%y = particles(n)%y + &
[1929]	788	MERGE( rand_contr, SIGN( dy, rand_contr ), &
	789	ABS( rand_contr ) < dy &
	790	)
[1359]	791	ENDIF
	792	IF ( psb(particles(n)%group) /= pst(particles(n)%group) ) THEN
[1929]	793	rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
	794	pdz(particles(n)%group)
[1359]	795	particles(n)%z = particles(n)%z + &
[1929]	796	MERGE( rand_contr, SIGN( dz, rand_contr ), &
	797	ABS( rand_contr ) < dz &
	798	)
[1359]	799	ENDIF
	800	ENDDO
[1]	801	!
[1929]	802	!-- Identify particles located outside the model domain and reflect
	803	!-- or absorb them if necessary.
[1359]	804	CALL lpm_boundary_conds( 'bottom/top' )
[1929]	805	!
	806	!-- Furthermore, remove particles located in topography. Note, as
	807	!-- the particle speed is still zero at this point, wall
	808	!-- reflection boundary conditions will not work in this case.
	809	particles => &
	810	grid_particles(kp,jp,ip)%particles(1:number_of_particles)
	811	DO n = local_start(kp,jp,ip), number_of_particles
	812	i = ( particles(n)%x + 0.5_wp * dx ) * ddx
	813	j = ( particles(n)%y + 0.5_wp * dy ) * ddy
	814	IF ( particles(n)%z <= zw(nzb_w_inner(j,i)) ) THEN
	815	particles(n)%particle_mask = .FALSE.
	816	deleted_particles = deleted_particles + 1
	817	ENDIF
	818	ENDDO
[1359]	819	ENDDO
	820	ENDDO
	821	ENDDO
[1]	822	!
[1359]	823	!-- Exchange particles between grid cells and processors
	824	CALL lpm_move_particle
	825	CALL lpm_exchange_horiz
[1]	826
[1359]	827	ENDIF
[1]	828	!
[1359]	829	!-- In case of random_start_position, delete particles identified by
	830	!-- lpm_exchange_horiz and lpm_boundary_conds. Then sort particles into blocks,
	831	!-- which is needed for a fast interpolation of the LES fields on the particle
	832	!-- position.
	833	CALL lpm_pack_all_arrays
[1]	834
	835	!
[1359]	836	!-- Determine maximum number of particles (i.e., all possible particles that
	837	!-- have been allocated) and the current number of particles
	838	DO ip = nxl, nxr
	839	DO jp = nys, nyn
	840	DO kp = nzb+1, nzt
	841	maximum_number_of_particles = maximum_number_of_particles &
	842	+ SIZE(grid_particles(kp,jp,ip)%particles)
	843	number_of_particles = number_of_particles &
	844	+ prt_count(kp,jp,ip)
[1]	845	ENDDO
[1359]	846	ENDDO
	847	ENDDO
[1]	848	!
[1822]	849	!-- Calculate the number of particles of the total domain
[1]	850	#if defined( __parallel )
[1359]	851	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	852	CALL MPI_ALLREDUCE( number_of_particles, total_number_of_particles, 1, &
	853	MPI_INTEGER, MPI_SUM, comm2d, ierr )
[1]	854	#else
[1359]	855	total_number_of_particles = number_of_particles
[1]	856	#endif
	857
[1359]	858	RETURN
[1]	859
[1359]	860	END SUBROUTINE lpm_create_particle
[336]	861
[1871]	862	SUBROUTINE lpm_init_aerosols(local_start)
	863
	864	USE arrays_3d, &
	865	ONLY: hyp, pt, q
	866
	867	USE cloud_parameters, &
[2122]	868	ONLY: l_d_rv, rho_l, r_v
[1871]	869
	870	USE constants, &
	871	ONLY: pi
	872
	873	USE kinds
	874
	875	USE particle_attributes, &
	876	ONLY: init_aerosol_probabilistic, molecular_weight_of_solute, &
	877	molecular_weight_of_water, n1, n2, n3, rho_s, rm1, rm2, rm3, &
	878	s1, s2, s3, vanthoff
	879
	880	IMPLICIT NONE
	881
	882	REAL(wp), DIMENSION(:), ALLOCATABLE :: cdf !< CDF of aerosol spectrum
	883	REAL(wp), DIMENSION(:), ALLOCATABLE :: r_temp !< dry aerosol radius spectrum
	884
[2122]	885	REAL(wp) :: afactor !< curvature effects
[1871]	886	REAL(wp) :: bfactor !< solute effects
	887	REAL(wp) :: dr !< width of radius bin
	888	REAL(wp) :: e_a !< vapor pressure
	889	REAL(wp) :: e_s !< saturation vapor pressure
	890	REAL(wp) :: n_init !< sum of all aerosol concentrations
	891	REAL(wp) :: pdf !< PDF of aerosol spectrum
	892	REAL(wp) :: rmin = 1.0e-8_wp !< minimum aerosol radius
	893	REAL(wp) :: rmax = 1.0e-6_wp !< maximum aerosol radius
	894	REAL(wp) :: rs_rand !< random number
	895	REAL(wp) :: r_mid !< mean radius
[2122]	896	REAL(wp) :: sigma !< surface tension
[1871]	897	REAL(wp) :: t_int !< temperature
	898	REAL(wp) :: weight_sum !< sum of all weighting factors
	899
[1890]	900	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) :: local_start !<
[1871]	901
	902	INTEGER(iwp) :: n !<
	903	INTEGER(iwp) :: nn !<
	904	INTEGER(iwp) :: no_bins = 999 !< number of bins
	905	INTEGER(iwp) :: ip !<
	906	INTEGER(iwp) :: jp !<
	907	INTEGER(iwp) :: kp !<
	908
	909	LOGICAL :: new_pdf = .FALSE. !< check if aerosol PDF has to be recalculated
	910
	911	!
	912	!-- Compute aerosol background distribution
	913	IF ( init_aerosol_probabilistic ) THEN
	914	ALLOCATE( cdf(0:no_bins), r_temp(0:no_bins) )
	915	DO n = 0, no_bins
	916	r_temp(n) = EXP( LOG(rmin) + ( LOG(rmax) - LOG(rmin ) ) / &
	917	REAL(no_bins, KIND=wp) * REAL(n, KIND=wp) )
	918
	919	cdf(n) = 0.0_wp
	920	n_init = n1 + n2 + n3
	921	IF ( n1 > 0.0_wp ) THEN
	922	cdf(n) = cdf(n) + n1 / n_init * ( 0.5_wp + 0.5_wp * &
	923	ERF( LOG( r_temp(n) / rm1 ) / &
	924	( SQRT(2.0_wp) * LOG(s1) ) &
	925	) )
	926	ENDIF
	927	IF ( n2 > 0.0_wp ) THEN
	928	cdf(n) = cdf(n) + n2 / n_init * ( 0.5_wp + 0.5_wp * &
	929	ERF( LOG( r_temp(n) / rm2 ) / &
	930	( SQRT(2.0_wp) * LOG(s2) ) &
	931	) )
	932	ENDIF
	933	IF ( n3 > 0.0_wp ) THEN
	934	cdf(n) = cdf(n) + n3 / n_init * ( 0.5_wp + 0.5_wp * &
	935	ERF( LOG( r_temp(n) / rm3 ) / &
	936	( SQRT(2.0_wp) * LOG(s3) ) &
	937	) )
	938	ENDIF
	939
	940	ENDDO
	941	ENDIF
	942
	943	DO ip = nxl, nxr
	944	DO jp = nys, nyn
	945	DO kp = nzb+1, nzt
	946
	947	number_of_particles = prt_count(kp,jp,ip)
	948	IF ( number_of_particles <= 0 ) CYCLE
	949	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
	950	!
	951	!-- Initialize the aerosols with a predefined spectral distribution
	952	!-- of the dry radius (logarithmically increasing bins) and a varying
	953	!-- weighting factor
	954	IF ( .NOT. init_aerosol_probabilistic ) THEN
	955
	956	new_pdf = .FALSE.
	957	IF ( .NOT. ALLOCATED( r_temp ) ) THEN
	958	new_pdf = .TRUE.
	959	ELSE
	960	IF ( SIZE( r_temp ) .NE. &
	961	number_of_particles - local_start(kp,jp,ip) + 2 ) THEN
	962	new_pdf = .TRUE.
	963	DEALLOCATE( r_temp )
	964	ENDIF
	965	ENDIF
	966
	967	IF ( new_pdf ) THEN
	968
	969	no_bins = number_of_particles + 1 - local_start(kp,jp,ip)
	970	ALLOCATE( r_temp(0:no_bins) )
	971
	972	DO n = 0, no_bins
	973	r_temp(n) = EXP( LOG(rmin) + ( LOG(rmax) - LOG(rmin ) ) / &
	974	REAL(no_bins, KIND=wp) * &
	975	REAL(n, KIND=wp) )
	976	ENDDO
	977
	978	ENDIF
	979
	980	!
	981	!-- Calculate radius and concentration of each aerosol
	982	DO n = local_start(kp,jp,ip), number_of_particles
	983
	984	nn = n - local_start(kp,jp,ip)
	985
	986	r_mid = SQRT( r_temp(nn) * r_temp(nn+1) )
	987	dr = r_temp(nn+1) - r_temp(nn)
	988
	989	pdf = 0.0_wp
	990	n_init = n1 + n2 + n3
	991	IF ( n1 > 0.0_wp ) THEN
	992	pdf = pdf + n1 / n_init * ( 1.0_wp / ( r_mid * LOG(s1) * &
	993	SQRT( 2.0_wp * pi ) &
	994	) * &
	995	EXP( -( LOG( r_mid / rm1 ) )**2 / &
	996	( 2.0_wp * LOG(s1)**2 ) &
	997	) &
	998	)
	999	ENDIF
	1000	IF ( n2 > 0.0_wp ) THEN
	1001	pdf = pdf + n2 / n_init * ( 1.0_wp / ( r_mid * LOG(s2) * &
	1002	SQRT( 2.0_wp * pi ) &
	1003	) * &
	1004	EXP( -( LOG( r_mid / rm2 ) )**2 / &
	1005	( 2.0_wp * LOG(s2)**2 ) &
	1006	) &
	1007	)
	1008	ENDIF
	1009	IF ( n3 > 0.0_wp ) THEN
	1010	pdf = pdf + n3 / n_init * ( 1.0_wp / ( r_mid * LOG(s3) * &
	1011	SQRT( 2.0_wp * pi ) &
	1012	) * &
	1013	EXP( -( LOG( r_mid / rm3 ) )**2 / &
	1014	( 2.0_wp * LOG(s3)**2 ) &
	1015	) &
	1016	)
	1017	ENDIF
	1018
	1019	particles(n)%rvar2 = r_mid
	1020	particles(n)%weight_factor = pdf * dr
	1021
	1022	END DO
	1023	!
	1024	!-- Adjust weighting factors to initialize the same number of aerosols
	1025	!-- in every grid box
	1026	weight_sum = SUM(particles(local_start(kp,jp,ip):number_of_particles)%weight_factor)
	1027
	1028	particles(local_start(kp,jp,ip):number_of_particles)%weight_factor = &
	1029	particles(local_start(kp,jp,ip):number_of_particles)%weight_factor / &
	1030	weight_sum * initial_weighting_factor * ( no_bins + 1 )
	1031
	1032	ENDIF
	1033	!
	1034	!-- Initialize the aerosols with a predefined weighting factor but
	1035	!-- a randomly choosen dry radius
	1036	IF ( init_aerosol_probabilistic ) THEN
	1037
	1038	DO n = local_start(kp,jp,ip), number_of_particles !only new particles
	1039
	1040	rs_rand = -1.0_wp
	1041	DO WHILE ( rs_rand .LT. cdf(0) .OR. rs_rand .GE. cdf(no_bins) )
	1042	rs_rand = random_function( iran_part )
	1043	ENDDO
	1044	!
	1045	!-- Determine aerosol dry radius by a random number generator
	1046	DO nn = 0, no_bins-1
	1047	IF ( cdf(nn) .LE. rs_rand .AND. cdf(nn+1) .GT. rs_rand ) THEN
	1048	particles(n)%rvar2 = r_temp(nn) + ( r_temp(nn+1) - r_temp(nn) ) / &
	1049	( cdf(nn+1) - cdf(nn) ) * ( rs_rand - cdf(nn) )
	1050	EXIT
	1051	ENDIF
	1052	ENDDO
	1053
	1054	ENDDO
	1055
	1056	ENDIF
	1057
	1058	!
	1059	!-- Set particle radius to equilibrium radius based on the environmental
	1060	!-- supersaturation (Khvorostyanov and Curry, 2007, JGR). This avoids
	1061	!-- the sometimes lengthy growth toward their equilibrium radius within
	1062	!-- the simulation.
	1063	t_int = pt(kp,jp,ip) * ( hyp(kp) / 100000.0_wp )**0.286_wp
	1064
	1065	e_s = 611.0_wp * EXP( l_d_rv * ( 3.6609E-3_wp - 1.0_wp / t_int ) )
	1066	e_a = q(kp,jp,ip) * hyp(kp) / ( 0.378_wp * q(kp,jp,ip) + 0.622_wp )
	1067
[2122]	1068	sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
	1069	afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
	1070
	1071	bfactor = vanthoff * molecular_weight_of_water * &
	1072	rho_s / ( molecular_weight_of_solute * rho_l )
[1871]	1073	!
[1890]	1074	!-- The formula is only valid for subsaturated environments. For
[2122]	1075	!-- supersaturations higher than -5 %, the supersaturation is set to -5%.
	1076	IF ( e_a / e_s >= 0.95_wp ) e_a = 0.95_wp * e_s
[1871]	1077
[2122]	1078	DO n = local_start(kp,jp,ip), number_of_particles !only new particles
	1079	!
	1080	!-- For details on this equation, see Eq. (14) of Khvorostyanov and
	1081	!-- Curry (2007, JGR)
	1082	particles(n)%radius = bfactor*0.3333333_wp &
	1083	particles(n)%rvar2 / ( 1.0_wp - e_a / e_s )**0.3333333_wp / &
	1084	( 1.0_wp + ( afactor / ( 3.0_wp * bfactor*0.3333333_wp &
	1085	particles(n)%rvar2 ) ) / &
	1086	( 1.0_wp - e_a / e_s )**0.6666666_wp &
	1087	)
[1871]	1088
[2122]	1089	ENDDO
[1871]	1090
	1091	ENDDO
	1092	ENDDO
	1093	ENDDO
	1094	!
	1095	!-- Deallocate used arrays
	1096	IF ( ALLOCATED(r_temp) ) DEALLOCATE( r_temp )
	1097	IF ( ALLOCATED(cdf) ) DEALLOCATE( cdf )
	1098
	1099	END SUBROUTINE lpm_init_aerosols
	1100
[1359]	1101	END MODULE lpm_init_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |