Home

Context Navigation

source: palm/trunk/SOURCE/lpm_droplet_collision.f90 @ 1320

Last change on this file since 1320 was 1320, checked in by raasch, 10 years ago
ONLY-attribute added to USE-statements, kind-parameters added to all INTEGER and REAL declaration statements, kinds are defined in new module kinds, old module precision_kind is removed, revision history before 2012 removed, comment fields (!:) to be used for variable explanations added to all variable declaration statements
Property svn:keywords set to `Id`
File size: 30.1 KB

Rev	Line
[849]	1	SUBROUTINE lpm_droplet_collision
	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	!
[849]	20	! Current revisions:
	21	! ------------------
[1320]	22	! ONLY-attribute added to USE-statements,
	23	! kind-parameters added to all INTEGER and REAL declaration statements,
	24	! kinds are defined in new module kinds,
	25	! revision history before 2012 removed,
	26	! comment fields (!:) to be used for variable explanations added to
	27	! all variable declaration statements
[1072]	28	!
	29	! Former revisions:
	30	! -----------------
	31	! $Id: lpm_droplet_collision.f90 1320 2014-03-20 08:40:49Z raasch $
	32	!
[1093]	33	! 1092 2013-02-02 11:24:22Z raasch
	34	! unused variables removed
	35	!
[1072]	36	! 1071 2012-11-29 16:54:55Z franke
[1071]	37	! Calculation of Hall and Wang kernel now uses collision-coalescence formulation
	38	! proposed by Wang instead of the continuous collection equation (for more
	39	! information about new method see PALM documentation)
	40	! Bugfix: message identifiers added
[849]	41	!
[1037]	42	! 1036 2012-10-22 13:43:42Z raasch
	43	! code put under GPL (PALM 3.9)
	44	!
[850]	45	! 849 2012-03-15 10:35:09Z raasch
	46	! initial revision (former part of advec_particles)
[849]	47	!
[850]	48	!
[849]	49	! Description:
	50	! ------------
[1071]	51	! Calculates change in droplet radius by collision. Droplet collision is
[849]	52	! calculated for each grid box seperately. Collision is parameterized by
	53	! using collision kernels. Three different kernels are available:
	54	! PALM kernel: Kernel is approximated using a method from Rogers and
	55	! Yau (1989, A Short Course in Cloud Physics, Pergamon Press).
	56	! All droplets smaller than the treated one are represented by
	57	! one droplet with mean features. Collision efficiencies are taken
	58	! from the respective table in Rogers and Yau.
	59	! Hall kernel: Kernel from Hall (1980, J. Atmos. Sci., 2486-2507), which
	60	! considers collision due to pure gravitational effects.
	61	! Wang kernel: Beside gravitational effects (treated with the Hall-kernel) also
	62	! the effects of turbulence on the collision are considered using
	63	! parameterizations of Ayala et al. (2008, New J. Phys., 10,
	64	! 075015) and Wang and Grabowski (2009, Atmos. Sci. Lett., 10,
	65	! 1-8). This kernel includes three possible effects of turbulence:
	66	! the modification of the relative velocity between the droplets,
	67	! the effect of preferential concentration, and the enhancement of
	68	! collision efficiencies.
	69	!------------------------------------------------------------------------------!
	70
[1320]	71	USE arrays_3d, &
	72	ONLY: diss, ql, ql_v, ql_vp, u, v, w, zu, zw
[849]	73
[1320]	74	USE cloud_parameters, &
	75	ONLY: effective_coll_efficiency
	76
	77	USE constants, &
	78	ONLY: pi
	79
	80	USE control_parameters, &
	81	ONLY: dt_3d, message_string, u_gtrans, v_gtrans, dz
	82
	83	USE cpulog, &
	84	ONLY: cpu_log, log_point_s
	85
	86	USE grid_variables, &
	87	ONLY: ddx, dx, ddy, dy
	88
	89	USE indices, &
	90	ONLY: nxl, nxr, nyn, nys, nzb, nzt
	91
	92	USE kinds
	93
	94	USE lpm_collision_kernels_mod, &
	95	ONLY: ckernel, collision_efficiency_rogers, recalculate_kernel
	96
	97	USE particle_attributes, &
	98	ONLY: deleted_particles, dissipation_classes, hall_kernel, &
	99	palm_kernel, particles, particle_mask, particle_type, &
	100	prt_count, prt_start_index, use_kernel_tables, wang_kernel
	101
[849]	102	IMPLICIT NONE
	103
[1320]	104	INTEGER(iwp) :: eclass !:
	105	INTEGER(iwp) :: i !:
	106	INTEGER(iwp) :: ii !:
	107	INTEGER(iwp) :: inc !:
	108	INTEGER(iwp) :: is !:
	109	INTEGER(iwp) :: j !:
	110	INTEGER(iwp) :: jj !:
	111	INTEGER(iwp) :: js !:
	112	INTEGER(iwp) :: k !:
	113	INTEGER(iwp) :: kk !:
	114	INTEGER(iwp) :: n !:
	115	INTEGER(iwp) :: pse !:
	116	INTEGER(iwp) :: psi !:
	117	INTEGER(iwp) :: rclass_l !:
	118	INTEGER(iwp) :: rclass_s !:
[849]	119
[1320]	120	REAL(wp) :: aa !:
	121	REAL(wp) :: bb !:
	122	REAL(wp) :: cc !:
	123	REAL(wp) :: dd !:
	124	REAL(wp) :: ddV !:
	125	REAL(wp) :: delta_r !:
	126	REAL(wp) :: delta_v !:
	127	REAL(wp) :: epsilon !:
	128	REAL(wp) :: gg !:
	129	REAL(wp) :: mean_r !:
	130	REAL(wp) :: ql_int !:
	131	REAL(wp) :: ql_int_l !:
	132	REAL(wp) :: ql_int_u !:
	133	REAL(wp) :: r3 !:
	134	REAL(wp) :: sl_r3 !:
	135	REAL(wp) :: sl_r4 !:
	136	REAL(wp) :: sum1 !:
	137	REAL(wp) :: sum2 !:
	138	REAL(wp) :: sum3 !:
	139	REAL(wp) :: u_int !:
	140	REAL(wp) :: u_int_l !:
	141	REAL(wp) :: u_int_u !:
	142	REAL(wp) :: v_int !:
	143	REAL(wp) :: v_int_l !:
	144	REAL(wp) :: v_int_u !:
	145	REAL(wp) :: w_int !:
	146	REAL(wp) :: w_int_l !:
	147	REAL(wp) :: w_int_u !:
	148	REAL(wp) :: x !:
	149	REAL(wp) :: y !:
[849]	150
[1320]	151	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad !:
	152	REAL(wp), DIMENSION(:), ALLOCATABLE :: weight !:
[849]	153
	154
[1320]	155	TYPE(particle_type) :: tmp_particle !:
	156
	157
	158
[849]	159	CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'start' )
	160
	161	DO i = nxl, nxr
	162	DO j = nys, nyn
	163	DO k = nzb+1, nzt
	164	!
	165	!-- Collision requires at least two particles in the box
	166	IF ( prt_count(k,j,i) > 1 ) THEN
	167	!
	168	!-- First, sort particles within the gridbox by their size,
	169	!-- using Shell's method (see Numerical Recipes)
	170	!-- NOTE: In case of using particle tails, the re-sorting of
	171	!-- ---- tails would have to be included here!
	172	psi = prt_start_index(k,j,i) - 1
	173	inc = 1
	174	DO WHILE ( inc <= prt_count(k,j,i) )
	175	inc = 3 * inc + 1
	176	ENDDO
	177
	178	DO WHILE ( inc > 1 )
	179	inc = inc / 3
	180	DO is = inc+1, prt_count(k,j,i)
	181	tmp_particle = particles(psi+is)
	182	js = is
	183	DO WHILE ( particles(psi+js-inc)%radius > &
	184	tmp_particle%radius )
	185	particles(psi+js) = particles(psi+js-inc)
	186	js = js - inc
	187	IF ( js <= inc ) EXIT
	188	ENDDO
	189	particles(psi+js) = tmp_particle
	190	ENDDO
	191	ENDDO
	192
	193	psi = prt_start_index(k,j,i)
	194	pse = psi + prt_count(k,j,i)-1
	195
	196	!
	197	!-- Now apply the different kernels
	198	IF ( ( hall_kernel .OR. wang_kernel ) .AND. &
	199	use_kernel_tables ) THEN
	200	!
	201	!-- Fast method with pre-calculated efficiencies for
	202	!-- discrete radius- and dissipation-classes.
	203	!
	204	!-- Determine dissipation class index of this gridbox
	205	IF ( wang_kernel ) THEN
	206	eclass = INT( diss(k,j,i) * 1.0E4 / 1000.0 * &
	207	dissipation_classes ) + 1
	208	epsilon = diss(k,j,i)
	209	ELSE
	210	epsilon = 0.0
	211	ENDIF
	212	IF ( hall_kernel .OR. epsilon * 1.0E4 < 0.001 ) THEN
	213	eclass = 0 ! Hall kernel is used
	214	ELSE
	215	eclass = MIN( dissipation_classes, eclass )
	216	ENDIF
	217
[1071]	218	!
	219	!-- Droplet collision are calculated using collision-coalescence
	220	!-- formulation proposed by Wang (see PALM documentation)
	221	!-- Since new radii after collision are defined by radii of all
	222	!-- droplets before collision, temporary fields for new radii and
	223	!-- weighting factors are needed
	224	ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
[849]	225
[1071]	226	rad = 0.0
	227	weight = 0.0
	228
	229	DO n = psi, pse, 1
	230
	231	sum1 = 0.0
	232	sum2 = 0.0
	233	sum3 = 0.0
	234
[849]	235	rclass_l = particles(n)%class
	236	!
[1071]	237	!-- Mass added due to collisions with smaller droplets
[849]	238	DO is = psi, n-1
	239
	240	rclass_s = particles(is)%class
[1071]	241	sum1 = sum1 + ( particles(is)%radius*3 &
	242	ckernel(rclass_l,rclass_s,eclass) * &
	243	particles(is)%weight_factor )
[849]	244
	245	ENDDO
	246	!
[1071]	247	!-- Rate of collisions with larger droplets
	248	DO is = n+1, pse
[849]	249
[1071]	250	rclass_s = particles(is)%class
	251	sum2 = sum2 + ( ckernel(rclass_l,rclass_s,eclass) * &
	252	particles(is)%weight_factor )
[849]	253
[1071]	254	ENDDO
[849]	255
[1071]	256	r3 = particles(n)%radius**3
	257	ddV = ddx * ddy / dz
	258	is = prt_start_index(k,j,i)
[849]	259	!
[1071]	260	!-- Change of the current weighting factor
	261	sum3 = 1 - dt_3d * ddV * &
	262	ckernel(rclass_l,rclass_l,eclass) * &
	263	( particles(n)%weight_factor - 1 ) * 0.5 - &
	264	dt_3d * ddV * sum2
	265	weight(n-is+1) = particles(n)%weight_factor * sum3
	266	!
	267	!-- Change of the current droplet radius
	268	rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/&
	269	sum3 )**0.33333333333333
[849]	270
[1071]	271	IF ( weight(n-is+1) < 0.0 ) THEN
	272	WRITE( message_string, * ) 'negative weighting', &
	273	'factor: ', weight(n-is+1)
	274	CALL message( 'lpm_droplet_collision', 'PA0028', &
	275	2, 2, -1, 6, 1 )
	276	ENDIF
[849]	277
[1071]	278	ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) &
	279	* rad(n-is+1)**3
[849]	280
[1071]	281	ENDDO
[849]	282
[1071]	283	particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
	284	particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
[849]	285
[1071]	286	DEALLOCATE(rad, weight)
[849]	287
	288	ELSEIF ( ( hall_kernel .OR. wang_kernel ) .AND. &
	289	.NOT. use_kernel_tables ) THEN
	290	!
	291	!-- Collision efficiencies are calculated for every new
	292	!-- grid box. First, allocate memory for kernel table.
	293	!-- Third dimension is 1, because table is re-calculated for
	294	!-- every new dissipation value.
	295	ALLOCATE( ckernel(prt_start_index(k,j,i): &
	296	prt_start_index(k,j,i)+prt_count(k,j,i)-1, &
	297	prt_start_index(k,j,i): &
	298	prt_start_index(k,j,i)+prt_count(k,j,i)-1,1:1) )
	299	!
	300	!-- Now calculate collision efficiencies for this box
	301	CALL recalculate_kernel( i, j, k )
	302
[1071]	303	!
	304	!-- Droplet collision are calculated using collision-coalescence
	305	!-- formulation proposed by Wang (see PALM documentation)
	306	!-- Since new radii after collision are defined by radii of all
	307	!-- droplets before collision, temporary fields for new radii and
	308	!-- weighting factors are needed
	309	ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
[849]	310
[1071]	311	rad = 0.0
	312	weight = 0.0
	313
	314	DO n = psi, pse, 1
	315
	316	sum1 = 0.0
	317	sum2 = 0.0
	318	sum3 = 0.0
[849]	319	!
[1071]	320	!-- Mass added due to collisions with smaller droplets
[849]	321	DO is = psi, n-1
[1071]	322	sum1 = sum1 + ( particles(is)%radius*3 &
	323	ckernel(n,is,1) * &
	324	particles(is)%weight_factor )
	325	ENDDO
[849]	326	!
[1071]	327	!-- Rate of collisions with larger droplets
	328	DO is = n+1, pse
	329	sum2 = sum2 + ( ckernel(n,is,1) * &
	330	particles(is)%weight_factor )
[849]	331	ENDDO
	332
[1071]	333	r3 = particles(n)%radius**3
	334	ddV = ddx * ddy / dz
	335	is = prt_start_index(k,j,i)
[849]	336	!
[1071]	337	!-- Change of the current weighting factor
	338	sum3 = 1 - dt_3d * ddV * &
	339	ckernel(n,n,1) * &
	340	( particles(n)%weight_factor - 1 ) * 0.5 - &
	341	dt_3d * ddV * sum2
	342	weight(n-is+1) = particles(n)%weight_factor * sum3
[849]	343	!
[1071]	344	!-- Change of the current droplet radius
	345	rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/&
	346	sum3 )**0.33333333333333
[849]	347
[1071]	348	IF ( weight(n-is+1) < 0.0 ) THEN
	349	WRITE( message_string, * ) 'negative weighting', &
	350	'factor: ', weight(n-is+1)
	351	CALL message( 'lpm_droplet_collision', 'PA0037', &
	352	2, 2, -1, 6, 1 )
[849]	353	ENDIF
	354
[1071]	355	ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) &
	356	* rad(n-is+1)**3
[849]	357
	358	ENDDO
	359
[1071]	360	particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
	361	particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
[849]	362
[1071]	363	DEALLOCATE( rad, weight, ckernel )
	364
[849]	365	ELSEIF ( palm_kernel ) THEN
	366	!
	367	!-- PALM collision kernel
	368	!
	369	!-- Calculate the mean radius of all those particles which
	370	!-- are of smaller size than the current particle and
	371	!-- use this radius for calculating the collision efficiency
	372	DO n = psi+prt_count(k,j,i)-1, psi+1, -1
	373
	374	sl_r3 = 0.0
	375	sl_r4 = 0.0
	376
	377	DO is = n-1, psi, -1
	378	IF ( particles(is)%radius < particles(n)%radius ) &
	379	THEN
	380	sl_r3 = sl_r3 + particles(is)%weight_factor &
	381	* particles(is)%radius**3
	382	sl_r4 = sl_r4 + particles(is)%weight_factor &
	383	* particles(is)%radius**4
	384	ENDIF
	385	ENDDO
	386
	387	IF ( ( sl_r3 ) > 0.0 ) THEN
	388	mean_r = ( sl_r4 ) / ( sl_r3 )
	389
	390	CALL collision_efficiency_rogers( mean_r, &
	391	particles(n)%radius, &
	392	effective_coll_efficiency )
	393
	394	ELSE
	395	effective_coll_efficiency = 0.0
	396	ENDIF
	397
	398	IF ( effective_coll_efficiency > 1.0 .OR. &
	399	effective_coll_efficiency < 0.0 ) &
	400	THEN
	401	WRITE( message_string, * ) 'collision_efficien' , &
	402	'cy out of range:' ,effective_coll_efficiency
	403	CALL message( 'lpm_droplet_collision', 'PA0145', 2, &
	404	2, -1, 6, 1 )
	405	ENDIF
	406
	407	!
	408	!-- Interpolation of ...
	409	ii = particles(n)%x * ddx
	410	jj = particles(n)%y * ddy
	411	kk = ( particles(n)%z + 0.5 * dz ) / dz
	412
	413	x = particles(n)%x - ii * dx
	414	y = particles(n)%y - jj * dy
	415	aa = x2 + y2
	416	bb = ( dx - x )2 + y2
	417	cc = x2 + ( dy - y )2
	418	dd = ( dx - x )2 + ( dy - y )2
	419	gg = aa + bb + cc + dd
	420
	421	ql_int_l = ( (gg-aa) * ql(kk,jj,ii) + (gg-bb) * &
	422	ql(kk,jj,ii+1) &
	423	+ (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) * &
	424	ql(kk,jj+1,ii+1) &
	425	) / ( 3.0 * gg )
	426
	427	ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii) + (gg-bb) * &
	428	ql(kk+1,jj,ii+1) &
	429	+ (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) * &
	430	ql(kk+1,jj+1,ii+1) &
	431	) / ( 3.0 * gg )
	432
	433	ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *&
	434	( ql_int_u - ql_int_l )
	435
	436	!
	437	!-- Interpolate u velocity-component
	438	ii = ( particles(n)%x + 0.5 * dx ) * ddx
	439	jj = particles(n)%y * ddy
	440	kk = ( particles(n)%z + 0.5 * dz ) / dz ! only if eqist
	441
	442	IF ( ( particles(n)%z - zu(kk) ) > (0.5*dz) ) kk = kk+1
	443
	444	x = particles(n)%x + ( 0.5 - ii ) * dx
	445	y = particles(n)%y - jj * dy
	446	aa = x2 + y2
	447	bb = ( dx - x )2 + y2
	448	cc = x2 + ( dy - y )2
	449	dd = ( dx - x )2 + ( dy - y )2
	450	gg = aa + bb + cc + dd
	451
	452	u_int_l = ( (gg-aa) * u(kk,jj,ii) + (gg-bb) * &
	453	u(kk,jj,ii+1) &
	454	+ (gg-cc) * u(kk,jj+1,ii) + (gg-dd) * &
	455	u(kk,jj+1,ii+1) &
	456	) / ( 3.0 * gg ) - u_gtrans
	457	IF ( kk+1 == nzt+1 ) THEN
	458	u_int = u_int_l
	459	ELSE
	460	u_int_u = ( (gg-aa) * u(kk+1,jj,ii) + (gg-bb) * &
	461	u(kk+1,jj,ii+1) &
	462	+ (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) * &
	463	u(kk+1,jj+1,ii+1) &
	464	) / ( 3.0 * gg ) - u_gtrans
	465	u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz &
	466	* ( u_int_u - u_int_l )
	467	ENDIF
	468
	469	!
	470	!-- Same procedure for interpolation of the v velocity-com-
	471	!-- ponent (adopt index k from u velocity-component)
	472	ii = particles(n)%x * ddx
	473	jj = ( particles(n)%y + 0.5 * dy ) * ddy
	474
	475	x = particles(n)%x - ii * dx
	476	y = particles(n)%y + ( 0.5 - jj ) * dy
	477	aa = x2 + y2
	478	bb = ( dx - x )2 + y2
	479	cc = x2 + ( dy - y )2
	480	dd = ( dx - x )2 + ( dy - y )2
	481	gg = aa + bb + cc + dd
	482
	483	v_int_l = ( ( gg-aa ) * v(kk,jj,ii) + ( gg-bb ) * &
	484	v(kk,jj,ii+1) &
	485	+ ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) * &
	486	v(kk,jj+1,ii+1) &
	487	) / ( 3.0 * gg ) - v_gtrans
	488	IF ( kk+1 == nzt+1 ) THEN
	489	v_int = v_int_l
	490	ELSE
	491	v_int_u = ( (gg-aa) * v(kk+1,jj,ii) + (gg-bb) * &
	492	v(kk+1,jj,ii+1) &
	493	+ (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) * &
	494	v(kk+1,jj+1,ii+1) &
	495	) / ( 3.0 * gg ) - v_gtrans
	496	v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz &
	497	* ( v_int_u - v_int_l )
	498	ENDIF
	499
	500	!
	501	!-- Same procedure for interpolation of the w velocity-com-
	502	!-- ponent (adopt index i from v velocity-component)
	503	jj = particles(n)%y * ddy
	504	kk = particles(n)%z / dz
	505
	506	x = particles(n)%x - ii * dx
	507	y = particles(n)%y - jj * dy
	508	aa = x2 + y2
	509	bb = ( dx - x )2 + y2
	510	cc = x2 + ( dy - y )2
	511	dd = ( dx - x )2 + ( dy - y )2
	512	gg = aa + bb + cc + dd
	513
	514	w_int_l = ( ( gg-aa ) * w(kk,jj,ii) + ( gg-bb ) * &
	515	w(kk,jj,ii+1) &
	516	+ ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) * &
	517	w(kk,jj+1,ii+1) &
	518	) / ( 3.0 * gg )
	519	IF ( kk+1 == nzt+1 ) THEN
	520	w_int = w_int_l
	521	ELSE
	522	w_int_u = ( (gg-aa) * w(kk+1,jj,ii) + (gg-bb) * &
	523	w(kk+1,jj,ii+1) &
	524	+ (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) * &
	525	w(kk+1,jj+1,ii+1) &
	526	) / ( 3.0 * gg )
	527	w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz &
	528	* ( w_int_u - w_int_l )
	529	ENDIF
	530
	531	!
	532	!-- Change in radius due to collision
	533	delta_r = effective_coll_efficiency / 3.0 &
	534	* pi * sl_r3 * ddx * ddy / dz &
	535	* SQRT( ( u_int - particles(n)%speed_x )**2 &
	536	+ ( v_int - particles(n)%speed_y )**2 &
	537	+ ( w_int - particles(n)%speed_z )**2 &
	538	) * dt_3d
	539	!
	540	!-- Change in volume due to collision
	541	delta_v = particles(n)%weight_factor &
	542	* ( ( particles(n)%radius + delta_r )**3 &
	543	- particles(n)%radius**3 )
	544
	545	!
	546	!-- Check if collected particles provide enough LWC for
	547	!-- volume change of collector particle
	548	IF ( delta_v >= sl_r3 .AND. sl_r3 > 0.0 ) THEN
	549
	550	delta_r = ( ( sl_r3/particles(n)%weight_factor ) &
	551	+ particles(n)%radius3 )( 1./3. ) &
	552	- particles(n)%radius
	553
	554	DO is = n-1, psi, -1
	555	IF ( particles(is)%radius < &
	556	particles(n)%radius ) THEN
	557	particles(is)%weight_factor = 0.0
	558	particle_mask(is) = .FALSE.
	559	deleted_particles = deleted_particles + 1
	560	ENDIF
	561	ENDDO
	562
	563	ELSE IF ( delta_v < sl_r3 .AND. sl_r3 > 0.0 ) THEN
	564
	565	DO is = n-1, psi, -1
	566	IF ( particles(is)%radius < particles(n)%radius &
	567	.AND. sl_r3 > 0.0 ) THEN
	568	particles(is)%weight_factor = &
	569	( ( particles(is)%weight_factor &
	570	* ( particles(is)%radius**3 ) ) &
	571	- ( delta_v &
	572	* particles(is)%weight_factor &
	573	* ( particles(is)%radius**3 ) &
	574	/ sl_r3 ) ) &
	575	/ ( particles(is)%radius**3 )
	576
	577	IF ( particles(is)%weight_factor < 0.0 ) THEN
	578	WRITE( message_string, * ) 'negative ', &
	579	'weighting factor: ', &
	580	particles(is)%weight_factor
[1071]	581	CALL message( 'lpm_droplet_collision', &
	582	'PA0039', &
[849]	583	2, 2, -1, 6, 1 )
	584	ENDIF
	585	ENDIF
	586	ENDDO
	587
	588	ENDIF
	589
	590	particles(n)%radius = particles(n)%radius + delta_r
	591	ql_vp(k,j,i) = ql_vp(k,j,i) + &
	592	particles(n)%weight_factor * &
	593	( particles(n)%radius**3 )
	594	ENDDO
	595
	596	ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor &
	597	* particles(psi)%radius**3
	598
[1071]	599	ENDIF ! collision kernel
[849]	600
	601	ELSE IF ( prt_count(k,j,i) == 1 ) THEN
	602
	603	psi = prt_start_index(k,j,i)
[1071]	604
	605	!
	606	!-- Calculate change of weighting factor due to self collision
	607	IF ( ( hall_kernel .OR. wang_kernel ) .AND. &
	608	use_kernel_tables ) THEN
	609
	610	IF ( wang_kernel ) THEN
	611	eclass = INT( diss(k,j,i) * 1.0E4 / 1000.0 * &
	612	dissipation_classes ) + 1
	613	epsilon = diss(k,j,i)
	614	ELSE
	615	epsilon = 0.0
	616	ENDIF
	617	IF ( hall_kernel .OR. epsilon * 1.0E4 < 0.001 ) THEN
	618	eclass = 0 ! Hall kernel is used
	619	ELSE
	620	eclass = MIN( dissipation_classes, eclass )
	621	ENDIF
	622
	623	ddV = ddx * ddy / dz
	624	rclass_l = particles(psi)%class
	625	sum3 = 1 - dt_3d * ddV * &
	626	( ckernel(rclass_l,rclass_l,eclass) * &
	627	( particles(psi)%weight_factor-1 ) * 0.5 )
	628
	629	particles(psi)%radius = ( particles(psi)%radius**3 / &
	630	sum3 )**0.33333333333333
	631	particles(psi)%weight_factor = particles(psi)%weight_factor &
	632	* sum3
	633
	634	ELSE IF ( ( hall_kernel .OR. wang_kernel ) .AND. &
	635	.NOT. use_kernel_tables ) THEN
	636	!
	637	!-- Collision efficiencies are calculated for every new
	638	!-- grid box. First, allocate memory for kernel table.
	639	!-- Third dimension is 1, because table is re-calculated for
	640	!-- every new dissipation value.
	641	ALLOCATE( ckernel(psi:psi, psi:psi, 1:1) )
	642	!
	643	!-- Now calculate collision efficiencies for this box
	644	CALL recalculate_kernel( i, j, k )
	645
	646	ddV = ddx * ddy / dz
	647	sum3 = 1 - dt_3d * ddV * ( ckernel(psi,psi,1) * &
	648	( particles(psi)%weight_factor - 1 ) * 0.5 )
	649
	650	particles(psi)%radius = ( particles(psi)%radius**3 / &
	651	sum3 )**0.33333333333333
	652	particles(psi)%weight_factor = particles(psi)%weight_factor &
	653	* sum3
	654
	655	DEALLOCATE( ckernel )
	656	ENDIF
	657
	658	ql_vp(k,j,i) = particles(psi)%weight_factor * &
[849]	659	particles(psi)%radius**3
	660	ENDIF
	661
	662	!
	663	!-- Check if condensation of LWC was conserved during collision
	664	!-- process
	665	IF ( ql_v(k,j,i) /= 0.0 ) THEN
	666	IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001 .OR. &
	667	ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999 ) THEN
	668	WRITE( message_string, * ) 'LWC is not conserved during',&
	669	' collision! ', &
	670	'LWC after condensation: ', &
	671	ql_v(k,j,i), &
	672	' LWC after collision: ', &
	673	ql_vp(k,j,i)
[1071]	674	CALL message( 'lpm_droplet_collision', 'PA0040', &
	675	2, 2, -1, 6, 1 )
[849]	676	ENDIF
	677	ENDIF
	678
	679	ENDDO
	680	ENDDO
	681	ENDDO
	682
	683	CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'stop' )
	684
	685
	686	END SUBROUTINE lpm_droplet_collision

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |