Home

Context Navigation

source: palm/trunk/SOURCE/lpm_droplet_condensation.f90 @ 1822

Last change on this file since 1822 was 1822, checked in by hoffmann, 8 years ago
changes in LPM and bulk cloud microphysics
Property svn:keywords set to `Id`
File size: 23.8 KB

Rev	Line
[1682]	1	!> @file lpm_droplet_condensation.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	!
[1818]	16	! Copyright 1997-2016 Leibniz Universitaet Hannover
[1036]	17	!--------------------------------------------------------------------------------!
	18	!
[849]	19	! Current revisions:
	20	! ------------------
[1822]	21	! Unused variables removed.
[1360]	22	!
[1321]	23	! Former revisions:
	24	! -----------------
	25	! $Id: lpm_droplet_condensation.f90 1822 2016-04-07 07:49:42Z hoffmann $
	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	!
[1347]	34	! 1346 2014-03-27 13:18:20Z heinze
	35	! Bugfix: REAL constants provided with KIND-attribute especially in call of
	36	! intrinsic function like MAX, MIN, SIGN
	37	!
[1323]	38	! 1322 2014-03-20 16:38:49Z raasch
	39	! REAL constants defined as wp-kind
	40	!
[1321]	41	! 1320 2014-03-20 08:40:49Z raasch
[1320]	42	! ONLY-attribute added to USE-statements,
	43	! kind-parameters added to all INTEGER and REAL declaration statements,
	44	! kinds are defined in new module kinds,
	45	! comment fields (!:) to be used for variable explanations added to
	46	! all variable declaration statements
[1072]	47	!
[1319]	48	! 1318 2014-03-17 13:35:16Z raasch
	49	! module interfaces removed
	50	!
[1093]	51	! 1092 2013-02-02 11:24:22Z raasch
	52	! unused variables removed
	53	!
[1072]	54	! 1071 2012-11-29 16:54:55Z franke
[1071]	55	! Ventilation effect for evaporation of large droplets included
	56	! Check for unreasonable results included in calculation of Rosenbrock method
	57	! since physically unlikely results were observed and for the same
	58	! reason the first internal time step in Rosenbrock method should be < 1.0E02 in
	59	! case of evaporation
	60	! Unnecessary calculation of ql_int removed
	61	! Unnecessary calculations in Rosenbrock method (d2rdt2, drdt_m, dt_ros_last)
	62	! removed
	63	! Bugfix: factor in calculation of surface tension changed from 0.00155 to
	64	! 0.000155
[849]	65	!
[1037]	66	! 1036 2012-10-22 13:43:42Z raasch
	67	! code put under GPL (PALM 3.9)
	68	!
[850]	69	! 849 2012-03-15 10:35:09Z raasch
	70	! initial revision (former part of advec_particles)
[849]	71	!
[850]	72	!
[849]	73	! Description:
	74	! ------------
[1682]	75	!> Calculates change in droplet radius by condensation/evaporation, using
	76	!> either an analytic formula or by numerically integrating the radius growth
	77	!> equation including curvature and solution effects using Rosenbrocks method
	78	!> (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
	79	!> The analytical formula and growth equation follow those given in
	80	!> Rogers and Yau (A short course in cloud physics, 3rd edition, p. 102/103).
[849]	81	!------------------------------------------------------------------------------!
[1682]	82	SUBROUTINE lpm_droplet_condensation (ip,jp,kp)
	83
[849]	84
[1320]	85	USE arrays_3d, &
	86	ONLY: hyp, pt, q, ql_c, ql_v, zu
[849]	87
[1320]	88	USE cloud_parameters, &
[1822]	89	ONLY: bfactor, curvature_solution_effects, eps_ros, l_d_rv, l_v, &
	90	rho_l, r_v
[849]	91
[1320]	92	USE constants, &
	93	ONLY: pi
[849]	94
[1320]	95	USE control_parameters, &
[1822]	96	ONLY: dt_3d, dz, message_string, molecular_viscosity, rho_surface
	97
[1320]	98	USE cpulog, &
	99	ONLY: cpu_log, log_point_s
[849]	100
[1320]	101	USE grid_variables, &
[1822]	102	ONLY: dx, dy
[1071]	103
[1320]	104	USE lpm_collision_kernels_mod, &
	105	ONLY: rclass_lbound, rclass_ubound
[849]	106
[1320]	107	USE kinds
	108
	109	USE particle_attributes, &
[1359]	110	ONLY: block_offset, grid_particles, hall_kernel, number_of_particles, &
[1822]	111	particles, radius_classes, use_kernel_tables, wang_kernel
[1320]	112
	113
	114	IMPLICIT NONE
	115
[1682]	116	INTEGER(iwp) :: i !<
	117	INTEGER(iwp) :: ip !<
	118	INTEGER(iwp) :: internal_timestep_count !<
	119	INTEGER(iwp) :: j !<
	120	INTEGER(iwp) :: jp !<
	121	INTEGER(iwp) :: jtry !<
	122	INTEGER(iwp) :: k !<
	123	INTEGER(iwp) :: kp !<
	124	INTEGER(iwp) :: n !<
	125	INTEGER(iwp) :: nb !<
	126	INTEGER(iwp) :: ros_count !<
[1320]	127
[1682]	128	INTEGER(iwp), PARAMETER :: maxtry = 40 !<
[1320]	129
[1682]	130	INTEGER(iwp), DIMENSION(0:7) :: end_index !<
	131	INTEGER(iwp), DIMENSION(0:7) :: start_index !<
[1320]	132
[1682]	133	LOGICAL :: repeat !<
[1359]	134
[1682]	135	LOGICAL, DIMENSION(number_of_particles) :: flag_1 !<
[1359]	136
[1682]	137	REAL(wp) :: aa !<
	138	REAL(wp) :: afactor !<
	139	REAL(wp) :: arg !<
	140	REAL(wp) :: bb !<
	141	REAL(wp) :: cc !<
	142	REAL(wp) :: dd !<
	143	REAL(wp) :: ddenom !<
	144	REAL(wp) :: delta_r !<
	145	REAL(wp) :: drdt !<
	146	REAL(wp) :: drdt_ini !<
	147	REAL(wp) :: dt_ros !<
	148	REAL(wp) :: dt_ros_next !<
	149	REAL(wp) :: dt_ros_sum !<
	150	REAL(wp) :: dt_ros_sum_ini !<
	151	REAL(wp) :: d2rdtdr !<
	152	REAL(wp) :: errmax !<
	153	REAL(wp) :: err_ros !<
	154	REAL(wp) :: g1 !<
	155	REAL(wp) :: g2 !<
	156	REAL(wp) :: g3 !<
	157	REAL(wp) :: g4 !<
	158	REAL(wp) :: gg !<
	159	REAL(wp) :: pt_int !<
	160	REAL(wp) :: pt_int_l !<
	161	REAL(wp) :: pt_int_u !<
	162	REAL(wp) :: q_int !<
	163	REAL(wp) :: q_int_l !<
	164	REAL(wp) :: q_int_u !<
	165	REAL(wp) :: r_ros !<
	166	REAL(wp) :: r_ros_ini !<
	167	REAL(wp) :: sigma !<
	168	REAL(wp) :: x !<
	169	REAL(wp) :: y !<
	170	REAL(wp) :: re_p !<
[1320]	171
[849]	172	!-- Parameters for Rosenbrock method
[1682]	173	REAL(wp), PARAMETER :: a21 = 2.0_wp !<
	174	REAL(wp), PARAMETER :: a31 = 48.0_wp / 25.0_wp !<
	175	REAL(wp), PARAMETER :: a32 = 6.0_wp / 25.0_wp !<
	176	REAL(wp), PARAMETER :: b1 = 19.0_wp / 9.0_wp !<
	177	REAL(wp), PARAMETER :: b2 = 0.5_wp !<
	178	REAL(wp), PARAMETER :: b3 = 25.0_wp / 108.0_wp !<
	179	REAL(wp), PARAMETER :: b4 = 125.0_wp / 108.0_wp !<
	180	REAL(wp), PARAMETER :: c21 = -8.0_wp !<
	181	REAL(wp), PARAMETER :: c31 = 372.0_wp / 25.0_wp !<
	182	REAL(wp), PARAMETER :: c32 = 12.0_wp / 5.0_wp !<
	183	REAL(wp), PARAMETER :: c41 = -112.0_wp / 125.0_wp !<
	184	REAL(wp), PARAMETER :: c42 = -54.0_wp / 125.0_wp !<
	185	REAL(wp), PARAMETER :: c43 = -2.0_wp / 5.0_wp !<
	186	REAL(wp), PARAMETER :: errcon = 0.1296_wp !<
	187	REAL(wp), PARAMETER :: e1 = 17.0_wp / 54.0_wp !<
	188	REAL(wp), PARAMETER :: e2 = 7.0_wp / 36.0_wp !<
	189	REAL(wp), PARAMETER :: e3 = 0.0_wp !<
	190	REAL(wp), PARAMETER :: e4 = 125.0_wp / 108.0_wp !<
	191	REAL(wp), PARAMETER :: gam = 0.5_wp !<
	192	REAL(wp), PARAMETER :: grow = 1.5_wp !<
	193	REAL(wp), PARAMETER :: pgrow = -0.25_wp !<
	194	REAL(wp), PARAMETER :: pshrnk = -1.0_wp /3.0_wp !<
	195	REAL(wp), PARAMETER :: shrnk = 0.5_wp !<
[849]	196
[1682]	197	REAL(wp), DIMENSION(number_of_particles) :: afactor_v !<
	198	REAL(wp), DIMENSION(number_of_particles) :: diff_coeff_v !<
	199	REAL(wp), DIMENSION(number_of_particles) :: e_s !<
	200	REAL(wp), DIMENSION(number_of_particles) :: e_a !<
	201	REAL(wp), DIMENSION(number_of_particles) :: new_r !<
	202	REAL(wp), DIMENSION(number_of_particles) :: p_int !<
	203	REAL(wp), DIMENSION(number_of_particles) :: thermal_conductivity_v !<
	204	REAL(wp), DIMENSION(number_of_particles) :: t_int !<
	205	REAL(wp), DIMENSION(number_of_particles) :: xv !<
	206	REAL(wp), DIMENSION(number_of_particles) :: yv !<
	207	REAL(wp), DIMENSION(number_of_particles) :: zv !<
[849]	208
[1320]	209
[849]	210	CALL cpu_log( log_point_s(42), 'lpm_droplet_condens', 'start' )
	211
[1359]	212	start_index = grid_particles(kp,jp,ip)%start_index
	213	end_index = grid_particles(kp,jp,ip)%end_index
	214
	215	xv = particles(1:number_of_particles)%x
	216	yv = particles(1:number_of_particles)%y
	217	zv = particles(1:number_of_particles)%z
	218
	219	DO nb = 0,7
	220
	221	i = ip + block_offset(nb)%i_off
	222	j = jp + block_offset(nb)%j_off
	223	k = kp + block_offset(nb)%k_off
	224
	225	DO n = start_index(nb), end_index(nb)
[849]	226	!
	227	!-- Interpolate temperature and humidity.
[1359]	228	x = xv(n) - i * dx
	229	y = yv(n) - j * dy
	230	aa = x2 + y2
	231	bb = ( dx - x )2 + y2
	232	cc = x2 + ( dy - y )2
	233	dd = ( dx - x )2 + ( dy - y )2
	234	gg = aa + bb + cc + dd
[849]	235
[1359]	236	pt_int_l = ( ( gg - aa ) * pt(k,j,i) + ( gg - bb ) * pt(k,j,i+1) &
	237	+ ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) &
	238	) / ( 3.0_wp * gg )
[849]	239
[1359]	240	pt_int_u = ( ( gg-aa ) * pt(k+1,j,i) + ( gg-bb ) * pt(k+1,j,i+1) &
	241	+ ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) &
	242	) / ( 3.0_wp * gg )
[849]	243
[1359]	244	pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * &
	245	( pt_int_u - pt_int_l )
[849]	246
[1359]	247	q_int_l = ( ( gg - aa ) * q(k,j,i) + ( gg - bb ) * q(k,j,i+1) &
	248	+ ( gg - cc ) * q(k,j+1,i) + ( gg - dd ) * q(k,j+1,i+1) &
	249	) / ( 3.0_wp * gg )
[849]	250
[1359]	251	q_int_u = ( ( gg-aa ) * q(k+1,j,i) + ( gg-bb ) * q(k+1,j,i+1) &
	252	+ ( gg-cc ) * q(k+1,j+1,i) + ( gg-dd ) * q(k+1,j+1,i+1) &
	253	) / ( 3.0_wp * gg )
[849]	254
[1359]	255	q_int = q_int_l + ( zv(n) - zu(k) ) / dz * &
	256	( q_int_u - q_int_l )
[849]	257
	258	!
[1359]	259	!-- Calculate real temperature and saturation vapor pressure
	260	p_int(n) = hyp(k) + ( particles(n)%z - zu(k) ) / dz * &
	261	( hyp(k+1)-hyp(k) )
	262	t_int(n) = pt_int * ( p_int(n) / 100000.0_wp )**0.286_wp
[849]	263
[1359]	264	e_s(n) = 611.0_wp * EXP( l_d_rv * ( 3.6609E-3_wp - 1.0_wp / &
	265	t_int(n) ) )
[849]	266
	267	!
[1359]	268	!-- Current vapor pressure
	269	e_a(n) = q_int * p_int(n) / ( 0.378_wp * q_int + 0.622_wp )
[849]	270
[1359]	271	ENDDO
	272	ENDDO
	273
	274	new_r = 0.0_wp
	275	flag_1 = .false.
	276
	277	DO n = 1, number_of_particles
[849]	278	!
	279	!-- Change in radius by condensation/evaporation
[1359]	280	IF ( particles(n)%radius >= 4.0E-5_wp .AND. &
	281	e_a(n)/e_s(n) < 1.0_wp ) THEN
[849]	282	!
[1359]	283	!-- Approximation for large radii, where curvature and solution effects
	284	!-- can be neglected but ventilation effect has to be included in case of
	285	!-- evaporation.
	286	!-- First calculate the droplet's Reynolds number
	287	re_p = 2.0_wp * particles(n)%radius * ABS( particles(n)%speed_z ) / &
	288	molecular_viscosity
[1071]	289	!
[1359]	290	!-- Ventilation coefficient (Rogers and Yau, 1989):
	291	IF ( re_p > 2.5_wp ) THEN
	292	afactor_v(n) = 0.78_wp + 0.28_wp * SQRT( re_p )
[1071]	293	ELSE
[1359]	294	afactor_v(n) = 1.0_wp + 0.09_wp * re_p
[1071]	295	ENDIF
[1359]	296	flag_1(n) = .TRUE.
	297	ELSEIF ( particles(n)%radius >= 1.0E-6_wp .OR. &
[1071]	298	.NOT. curvature_solution_effects ) THEN
	299	!
[1359]	300	!-- Approximation for larger radii in case that curvature and solution
	301	!-- effects are neglected and ventilation effects does not play a role
	302	afactor_v(n) = 1.0_wp
	303	flag_1(n) = .TRUE.
[849]	304	ENDIF
[1359]	305	ENDDO
[849]	306
[1359]	307	DO n = 1, number_of_particles
[849]	308	!
[1359]	309	!-- Thermal conductivity for water (from Rogers and Yau, Table 7.1),
	310	!-- diffusivity for water vapor (after Hall und Pruppacher, 1976)
	311	thermal_conductivity_v(n) = 7.94048E-05_wp * t_int(n) + 0.00227011_wp
	312	diff_coeff_v(n) = 0.211E-4_wp * &
	313	( t_int(n) / 273.15_wp )*1.94_wp ( 101325.0_wp / p_int(n))
	314
	315	IF(flag_1(n)) then
	316	arg = particles(n)%radius*2 + 2.0_wp dt_3d * afactor_v(n) * &
	317	( e_a(n) / e_s(n) - 1.0_wp ) / &
	318	( ( l_d_rv / t_int(n) - 1.0_wp ) * l_v * rho_l / t_int(n) / &
	319	thermal_conductivity_v(n) + &
	320	rho_l * r_v * t_int(n) / diff_coeff_v(n) / e_s(n) )
	321
	322	arg = MAX( arg, 1.0E-16_wp )
	323	new_r(n) = SQRT( arg )
	324	ENDIF
	325	ENDDO
	326
	327	DO n = 1, number_of_particles
	328	IF ( curvature_solution_effects .AND. &
	329	( ( particles(n)%radius < 1.0E-6_wp ) .OR. &
	330	( new_r(n) < 1.0E-6_wp ) ) ) THEN
	331	!
[849]	332	!-- Curvature and solutions effects are included in growth equation.
	333	!-- Change in Radius is calculated with 4th-order Rosenbrock method
	334	!-- for stiff o.d.e's with monitoring local truncation error to adjust
	335	!-- stepsize (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
	336	!-- For larger radii the simple analytic method (see ELSE) gives
	337	!-- almost the same results.
[1071]	338
	339	ros_count = 0
	340	repeat = .TRUE.
[849]	341	!
[1071]	342	!-- Carry out the Rosenbrock algorithm. In case of unreasonable results
	343	!-- the switch "repeat" will be set true and the algorithm will be carried
	344	!-- out again with the internal time step set to its initial (small) value.
[1359]	345	!-- Unreasonable results may occur if the external conditions, especially
	346	!-- the supersaturation, has significantly changed compared to the last
	347	!-- PALM timestep.
[1071]	348	DO WHILE ( repeat )
[849]	349
[1071]	350	repeat = .FALSE.
	351	!
[1359]	352	!-- Surface tension (Straka, 2009):
	353	sigma = 0.0761_wp - 0.000155_wp * ( t_int(n) - 273.15_wp )
[849]	354
[1071]	355	r_ros = particles(n)%radius
[1359]	356	dt_ros_sum = 0.0_wp ! internal integrated time (s)
[1071]	357	internal_timestep_count = 0
[849]	358
[1359]	359	ddenom = 1.0_wp / ( rho_l * r_v * t_int(n) / ( e_s(n) * &
	360	diff_coeff_v(n) ) + ( l_v / &
	361	( r_v * t_int(n) ) - 1.0_wp ) * &
	362	rho_l * l_v / ( thermal_conductivity_v(n) * &
	363	t_int(n) ) &
	364	)
[1071]	365
[1359]	366	afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int(n) )
[1071]	367
[849]	368	!
[1071]	369	!-- Take internal time step values from the end of last PALM time step
	370	dt_ros_next = particles(n)%rvar1
	371
[849]	372	!
[1071]	373	!-- Internal time step should not be > 1.0E-2 in case of evaporation
	374	!-- because larger values may lead to secondary solutions which are
	375	!-- physically unlikely
[1359]	376	IF ( dt_ros_next > 1.0E-2_wp .AND. e_a(n)/e_s(n) < 1.0_wp ) THEN
	377	dt_ros_next = 1.0E-3_wp
[1071]	378	ENDIF
[849]	379	!
[1071]	380	!-- If calculation of Rosenbrock method is repeated due to unreasonalble
	381	!-- results during previous try the initial internal time step has to be
	382	!-- reduced
	383	IF ( ros_count > 1 ) THEN
[1359]	384	dt_ros_next = dt_ros_next - ( 0.2_wp * dt_ros_next )
[1071]	385	ELSEIF ( ros_count > 5 ) THEN
[849]	386	!
[1071]	387	!-- Prevent creation of infinite loop
	388	message_string = 'ros_count > 5 in Rosenbrock method'
	389	CALL message( 'lpm_droplet_condensation', 'PA0018', 2, 2, &
	390	0, 6, 0 )
	391	ENDIF
	392
[849]	393	!
[1071]	394	!-- Internal time step must not be larger than PALM time step
	395	dt_ros = MIN( dt_ros_next, dt_3d )
	396	!
	397	!-- Integrate growth equation in time unless PALM time step is reached
	398	DO WHILE ( dt_ros_sum < dt_3d )
[849]	399
[1071]	400	internal_timestep_count = internal_timestep_count + 1
[849]	401
	402	!
[1071]	403	!-- Derivative at starting value
[1359]	404	drdt = ddenom / r_ros * ( e_a(n) / e_s(n) - 1.0_wp - afactor / &
	405	r_ros + bfactor / r_ros**3 )
[1071]	406	drdt_ini = drdt
	407	dt_ros_sum_ini = dt_ros_sum
	408	r_ros_ini = r_ros
[849]	409
	410	!
[1071]	411	!-- Calculate radial derivative of dr/dt
[1359]	412	d2rdtdr = ddenom * ( ( 1.0_wp - e_a(n)/e_s(n) ) / r_ros**2 + &
	413	2.0_wp * afactor / r_ros**3 - &
	414	4.0_wp * bfactor / r_ros**5 )
[849]	415	!
[1071]	416	!-- Adjust stepsize unless required accuracy is reached
	417	DO jtry = 1, maxtry+1
[849]	418
[1071]	419	IF ( jtry == maxtry+1 ) THEN
	420	message_string = 'maxtry > 40 in Rosenbrock method'
[1359]	421	CALL message( 'lpm_droplet_condensation', 'PA0347', 2, &
	422	2, 0, 6, 0 )
[1071]	423	ENDIF
[849]	424
[1359]	425	aa = 1.0_wp / ( gam * dt_ros ) - d2rdtdr
[1071]	426	g1 = drdt_ini / aa
	427	r_ros = r_ros_ini + a21 * g1
[1359]	428	drdt = ddenom / r_ros * ( e_a(n) / e_s(n) - 1.0_wp - &
	429	afactor / r_ros + &
[1071]	430	bfactor / r_ros**3 )
[849]	431
[1071]	432	g2 = ( drdt + c21 * g1 / dt_ros )&
	433	/ aa
	434	r_ros = r_ros_ini + a31 * g1 + a32 * g2
[1359]	435	drdt = ddenom / r_ros * ( e_a(n) / e_s(n) - 1.0_wp - &
	436	afactor / r_ros + &
[1071]	437	bfactor / r_ros**3 )
[849]	438
[1071]	439	g3 = ( drdt + &
	440	( c31 * g1 + c32 * g2 ) / dt_ros ) / aa
	441	g4 = ( drdt + &
	442	( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa
	443	r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4
[849]	444
[1071]	445	dt_ros_sum = dt_ros_sum_ini + dt_ros
[849]	446
[1071]	447	IF ( dt_ros_sum == dt_ros_sum_ini ) THEN
	448	message_string = 'zero stepsize in Rosenbrock method'
[1359]	449	CALL message( 'lpm_droplet_condensation', 'PA0348', 2, &
	450	2, 0, 6, 0 )
[1071]	451	ENDIF
[849]	452	!
[1071]	453	!-- Calculate error
[1359]	454	err_ros = e1 * g1 + e2 * g2 + e3 * g3 + e4 * g4
	455	errmax = 0.0_wp
[1071]	456	errmax = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros
[849]	457	!
[1071]	458	!-- Leave loop if accuracy is sufficient, otherwise try again
	459	!-- with a reduced stepsize
[1359]	460	IF ( errmax <= 1.0_wp ) THEN
[1071]	461	EXIT
	462	ELSE
[1359]	463	dt_ros = SIGN( MAX( ABS( 0.9_wp * dt_ros * &
	464	errmax**pshrnk ), &
	465	shrnk * ABS( dt_ros ) ), dt_ros )
[1071]	466	ENDIF
	467
	468	ENDDO ! loop for stepsize adjustment
	469
	470	!
	471	!-- Calculate next internal time step
	472	IF ( errmax > errcon ) THEN
[1359]	473	dt_ros_next = 0.9_wp * dt_ros * errmax**pgrow
[849]	474	ELSE
[1071]	475	dt_ros_next = grow * dt_ros
[849]	476	ENDIF
	477
[1071]	478	!
	479	!-- Estimated time step is reduced if the PALM time step is exceeded
	480	IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d ) THEN
	481	dt_ros = dt_3d - dt_ros_sum
	482	ELSE
	483	dt_ros = dt_ros_next
	484	ENDIF
[849]	485
[1071]	486	ENDDO
[849]	487	!
[1071]	488	!-- Store internal time step value for next PALM step
	489	particles(n)%rvar1 = dt_ros_next
[849]	490
[1359]	491	new_r(n) = r_ros
[849]	492	!
[1071]	493	!-- Radius should not fall below 1E-8 because Rosenbrock method may
	494	!-- lead to errors otherwise
[1359]	495	new_r(n) = MAX( new_r(n), 1.0E-8_wp )
[1071]	496	!
	497	!-- Check if calculated droplet radius change is reasonable since in
	498	!-- case of droplet evaporation the Rosenbrock method may lead to
	499	!-- secondary solutions which are physically unlikely.
	500	!-- Due to the solution effect the droplets may grow for relative
[1359]	501	!-- humidities below 100%, but change of radius should not be too
	502	!-- large. In case of unreasonable droplet growth the Rosenbrock
	503	!-- method is recalculated using a smaller initial time step.
[1071]	504	!-- Limiting values are tested for droplets down to 1.0E-7
[1359]	505	IF ( new_r(n) - particles(n)%radius >= 3.0E-7_wp .AND. &
	506	e_a(n)/e_s(n) < 0.97_wp ) THEN
[1071]	507	ros_count = ros_count + 1
	508	repeat = .TRUE.
[849]	509	ENDIF
	510
[1071]	511	ENDDO ! Rosenbrock method
[849]	512
	513	ENDIF
	514
[1359]	515	delta_r = new_r(n) - particles(n)%radius
[849]	516
	517	!
	518	!-- Sum up the change in volume of liquid water for the respective grid
	519	!-- volume (this is needed later in lpm_calc_liquid_water_content for
	520	!-- calculating the release of latent heat)
[1359]	521	i = ip
	522	j = jp
	523	k = kp
[849]	524	! only exact if equidistant
	525
[1359]	526	ql_c(k,j,i) = ql_c(k,j,i) + particles(n)%weight_factor * &
	527	rho_l * 1.33333333_wp * pi * &
	528	( new_r(n)3 - particles(n)%radius3 ) / &
[849]	529	( rho_surface * dx * dy * dz )
[1359]	530	IF ( ql_c(k,j,i) > 100.0_wp ) THEN
[849]	531	WRITE( message_string, * ) 'k=',k,' j=',j,' i=',i, &
	532	' ql_c=',ql_c(k,j,i), ' &part(',n,')%wf=', &
	533	particles(n)%weight_factor,' delta_r=',delta_r
	534	CALL message( 'lpm_droplet_condensation', 'PA0143', 2, 2, -1, 6, 1 )
	535	ENDIF
	536
	537	!
	538	!-- Change the droplet radius
[1359]	539	IF ( ( new_r(n) - particles(n)%radius ) < 0.0_wp .AND. &
	540	new_r(n) < 0.0_wp ) THEN
	541	WRITE( message_string, * ) '#1 k=',k,' j=',j,' i=',i, &
	542	' e_s=',e_s(n), ' e_a=',e_a(n),' t_int=',t_int(n), &
	543	' &delta_r=',delta_r, &
[849]	544	' particle_radius=',particles(n)%radius
	545	CALL message( 'lpm_droplet_condensation', 'PA0144', 2, 2, -1, 6, 1 )
	546	ENDIF
	547
	548	!
	549	!-- Sum up the total volume of liquid water (needed below for
	550	!-- re-calculating the weighting factors)
[1359]	551	ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * new_r(n)**3
[849]	552
[1359]	553	particles(n)%radius = new_r(n)
[849]	554
	555	!
	556	!-- Determine radius class of the particle needed for collision
[1359]	557	IF ( ( hall_kernel .OR. wang_kernel ) .AND. use_kernel_tables ) &
[849]	558	THEN
[1359]	559	particles(n)%class = ( LOG( new_r(n) ) - rclass_lbound ) / &
	560	( rclass_ubound - rclass_lbound ) * &
[849]	561	radius_classes
	562	particles(n)%class = MIN( particles(n)%class, radius_classes )
	563	particles(n)%class = MAX( particles(n)%class, 1 )
	564	ENDIF
	565
	566	ENDDO
	567
	568	CALL cpu_log( log_point_s(42), 'lpm_droplet_condens', 'stop' )
	569
	570
	571	END SUBROUTINE lpm_droplet_condensation

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |