Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |