Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |