Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |