Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |