Home

Context Navigation

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

Last change on this file since 1745 was 1683, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 27.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |