Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |