Home

Context Navigation

source: palm/trunk/SOURCE/microphysics.f90 @ 1828

Last change on this file since 1828 was 1823, checked in by hoffmann, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 81.6 KB

Rev	Line
[1682]	1	!> @file microphysics.f90
[1093]	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
[1093]	17	!--------------------------------------------------------------------------------!
	18	!
[1000]	19	! Current revisions:
[1092]	20	! ------------------
[1822]	21	!
[1823]	22	!
[1321]	23	! Former revisions:
	24	! -----------------
	25	! $Id: microphysics.f90 1823 2016-04-07 08:57:52Z maronga $
	26	!
[1823]	27	! 1822 2016-04-07 07:49:42Z hoffmann
	28	! Unused variables removed.
	29	! Kessler scheme integrated.
	30	!
[1692]	31	! 1691 2015-10-26 16:17:44Z maronga
	32	! Added new routine calc_precipitation_amount. The routine now allows to account
	33	! for precipitation due to sedimenation of cloud (fog) droplets
	34	!
[1683]	35	! 1682 2015-10-07 23:56:08Z knoop
	36	! Code annotations made doxygen readable
	37	!
[1647]	38	! 1646 2015-09-02 16:00:10Z hoffmann
	39	! Bugfix: Wrong computation of d_mean.
	40	!
[1362]	41	! 1361 2014-04-16 15:17:48Z hoffmann
	42	! Bugfix in sedimentation_rain: Index corrected.
	43	! Vectorized version of adjust_cloud added.
	44	! Little reformatting of the code.
	45	!
[1354]	46	! 1353 2014-04-08 15:21:23Z heinze
	47	! REAL constants provided with KIND-attribute
	48	!
[1347]	49	! 1346 2014-03-27 13:18:20Z heinze
	50	! Bugfix: REAL constants provided with KIND-attribute especially in call of
	51	! intrinsic function like MAX, MIN, SIGN
	52	!
[1335]	53	! 1334 2014-03-25 12:21:40Z heinze
	54	! Bugfix: REAL constants provided with KIND-attribute
	55	!
[1323]	56	! 1322 2014-03-20 16:38:49Z raasch
	57	! REAL constants defined as wp-kind
	58	!
[1321]	59	! 1320 2014-03-20 08:40:49Z raasch
[1320]	60	! ONLY-attribute added to USE-statements,
	61	! kind-parameters added to all INTEGER and REAL declaration statements,
	62	! kinds are defined in new module kinds,
	63	! comment fields (!:) to be used for variable explanations added to
	64	! all variable declaration statements
[1000]	65	!
[1242]	66	! 1241 2013-10-30 11:36:58Z heinze
	67	! hyp and rho have to be calculated at each time step if data from external
	68	! file LSF_DATA are used
	69	!
[1116]	70	! 1115 2013-03-26 18:16:16Z hoffmann
	71	! microphyical tendencies are calculated in microphysics_control in an optimized
	72	! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
	73	!
[1107]	74	! 1106 2013-03-04 05:31:38Z raasch
	75	! small changes in code formatting
	76	!
[1093]	77	! 1092 2013-02-02 11:24:22Z raasch
	78	! unused variables removed
	79	! file put under GPL
	80	!
[1066]	81	! 1065 2012-11-22 17:42:36Z hoffmann
	82	! Sedimentation process implemented according to Stevens and Seifert (2008).
[1115]	83	! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
[1066]	84	! and Stevens, 2010).
	85	!
[1054]	86	! 1053 2012-11-13 17:11:03Z hoffmann
	87	! initial revision
[1000]	88	!
	89	! Description:
	90	! ------------
[1682]	91	!> Calculate cloud microphysics according to the two moment bulk
	92	!> scheme by Seifert and Beheng (2006).
[1000]	93	!------------------------------------------------------------------------------!
[1682]	94	MODULE microphysics_mod
	95
[1000]	96
	97	PRIVATE
[1115]	98	PUBLIC microphysics_control
[1000]	99
[1115]	100	INTERFACE microphysics_control
	101	MODULE PROCEDURE microphysics_control
	102	MODULE PROCEDURE microphysics_control_ij
	103	END INTERFACE microphysics_control
[1022]	104
[1115]	105	INTERFACE adjust_cloud
	106	MODULE PROCEDURE adjust_cloud
	107	MODULE PROCEDURE adjust_cloud_ij
	108	END INTERFACE adjust_cloud
	109
[1000]	110	INTERFACE autoconversion
	111	MODULE PROCEDURE autoconversion
	112	MODULE PROCEDURE autoconversion_ij
	113	END INTERFACE autoconversion
	114
[1822]	115	INTERFACE autoconversion_kessler
	116	MODULE PROCEDURE autoconversion_kessler
	117	MODULE PROCEDURE autoconversion_kessler_ij
	118	END INTERFACE autoconversion_kessler
	119
[1000]	120	INTERFACE accretion
	121	MODULE PROCEDURE accretion
	122	MODULE PROCEDURE accretion_ij
	123	END INTERFACE accretion
[1005]	124
	125	INTERFACE selfcollection_breakup
	126	MODULE PROCEDURE selfcollection_breakup
	127	MODULE PROCEDURE selfcollection_breakup_ij
	128	END INTERFACE selfcollection_breakup
[1012]	129
	130	INTERFACE evaporation_rain
	131	MODULE PROCEDURE evaporation_rain
	132	MODULE PROCEDURE evaporation_rain_ij
	133	END INTERFACE evaporation_rain
	134
	135	INTERFACE sedimentation_cloud
	136	MODULE PROCEDURE sedimentation_cloud
	137	MODULE PROCEDURE sedimentation_cloud_ij
	138	END INTERFACE sedimentation_cloud
[1000]	139
[1012]	140	INTERFACE sedimentation_rain
	141	MODULE PROCEDURE sedimentation_rain
	142	MODULE PROCEDURE sedimentation_rain_ij
	143	END INTERFACE sedimentation_rain
	144
[1691]	145	INTERFACE calc_precipitation_amount
	146	MODULE PROCEDURE calc_precipitation_amount
	147	MODULE PROCEDURE calc_precipitation_amount_ij
	148	END INTERFACE calc_precipitation_amount
	149
	150
[1000]	151	CONTAINS
	152
	153
	154	!------------------------------------------------------------------------------!
[1682]	155	! Description:
	156	! ------------
	157	!> Call for all grid points
[1000]	158	!------------------------------------------------------------------------------!
[1115]	159	SUBROUTINE microphysics_control
[1022]	160
[1361]	161	USE arrays_3d, &
[1822]	162	ONLY: hyp, pt_init, zu
[1361]	163
	164	USE cloud_parameters, &
[1822]	165	ONLY: cp, hyrho, prr, pt_d_t, r_d, t_d_pt
[1361]	166
	167	USE control_parameters, &
	168	ONLY: call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, &
[1822]	169	g, intermediate_timestep_count, large_scale_forcing, &
	170	lsf_surf, microphysics_kessler, microphysics_seifert, &
	171	pt_surface, rho_surface,surface_pressure
[1361]	172
	173	USE indices, &
	174	ONLY: nzb, nzt
	175
[1320]	176	USE kinds
[1115]	177
[1361]	178	USE statistics, &
	179	ONLY: weight_pres
	180
[1115]	181	IMPLICIT NONE
	182
[1682]	183	INTEGER(iwp) :: k !<
[1115]	184
[1682]	185	REAL(wp) :: t_surface !<
[1361]	186
	187	IF ( large_scale_forcing .AND. lsf_surf ) THEN
	188	!
	189	!-- Calculate:
	190	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
	191	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
	192	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
	193	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
	194	DO k = nzb, nzt+1
	195	hyp(k) = surface_pressure * 100.0_wp * &
	196	( ( t_surface - g / cp * zu(k) ) / &
	197	t_surface )**(1.0_wp / 0.286_wp)
	198	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
	199	t_d_pt(k) = 1.0_wp / pt_d_t(k)
	200	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
[1115]	201	ENDDO
[1822]	202
[1361]	203	!
	204	!-- Compute reference density
	205	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
	206	ENDIF
[1115]	207
[1361]	208	!
	209	!-- Compute length of time step
	210	IF ( call_microphysics_at_all_substeps ) THEN
	211	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
	212	ELSE
	213	dt_micro = dt_3d
	214	ENDIF
	215
	216	!
[1822]	217	!-- Reset precipitation rate
	218	IF ( intermediate_timestep_count == 1 ) prr = 0.0_wp
	219
	220	!
[1361]	221	!-- Compute cloud physics
[1822]	222	IF ( microphysics_kessler ) THEN
	223
	224	CALL autoconversion_kessler
	225
	226	ELSEIF ( microphysics_seifert ) THEN
	227
[1361]	228	CALL adjust_cloud
	229	CALL autoconversion
	230	CALL accretion
	231	CALL selfcollection_breakup
	232	CALL evaporation_rain
	233	CALL sedimentation_rain
	234
[1822]	235	IF ( drizzle ) CALL sedimentation_cloud
[1361]	236
[1691]	237	ENDIF
	238
[1822]	239	CALL calc_precipitation_amount
	240
[1115]	241	END SUBROUTINE microphysics_control
	242
[1682]	243	!------------------------------------------------------------------------------!
	244	! Description:
	245	! ------------
	246	!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
	247	!> evaporation of rain drops due to too small or too big weights
	248	!> of rain drops (Stevens and Seifert, 2008).
	249	!------------------------------------------------------------------------------!
[1115]	250	SUBROUTINE adjust_cloud
	251
[1361]	252	USE arrays_3d, &
	253	ONLY: qr, nr
	254
	255	USE cloud_parameters, &
[1822]	256	ONLY: eps_sb, xrmin, xrmax, hyrho
[1361]	257
	258	USE cpulog, &
	259	ONLY: cpu_log, log_point_s
	260
	261	USE indices, &
[1822]	262	ONLY: nxl, nxr, nys, nyn, nzb_s_inner, nzt
[1361]	263
[1320]	264	USE kinds
[1022]	265
	266	IMPLICIT NONE
	267
[1682]	268	INTEGER(iwp) :: i !<
	269	INTEGER(iwp) :: j !<
	270	INTEGER(iwp) :: k !<
[1022]	271
[1361]	272	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
	273
[1022]	274	DO i = nxl, nxr
	275	DO j = nys, nyn
[1115]	276	DO k = nzb_s_inner(j,i)+1, nzt
[1361]	277	IF ( qr(k,j,i) <= eps_sb ) THEN
	278	qr(k,j,i) = 0.0_wp
	279	nr(k,j,i) = 0.0_wp
	280	ELSE
	281	IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN
	282	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
	283	ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN
	284	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
	285	ENDIF
	286	ENDIF
[1022]	287	ENDDO
	288	ENDDO
	289	ENDDO
	290
[1361]	291	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
	292
[1115]	293	END SUBROUTINE adjust_cloud
[1022]	294
[1106]	295
[1682]	296	!------------------------------------------------------------------------------!
	297	! Description:
	298	! ------------
	299	!> Autoconversion rate (Seifert and Beheng, 2006).
	300	!------------------------------------------------------------------------------!
[1000]	301	SUBROUTINE autoconversion
	302
[1361]	303	USE arrays_3d, &
	304	ONLY: diss, dzu, nr, qc, qr
	305
	306	USE cloud_parameters, &
	307	ONLY: a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3, &
	308	c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air, &
	309	nc_const, x0
	310
	311	USE control_parameters, &
	312	ONLY: dt_micro, rho_surface, turbulence
	313
	314	USE cpulog, &
	315	ONLY: cpu_log, log_point_s
	316
	317	USE grid_variables, &
	318	ONLY: dx, dy
	319
	320	USE indices, &
[1822]	321	ONLY: nxl, nxr, nys, nyn, nzb_s_inner, nzt
[1361]	322
[1320]	323	USE kinds
[1000]	324
	325	IMPLICIT NONE
	326
[1682]	327	INTEGER(iwp) :: i !<
	328	INTEGER(iwp) :: j !<
	329	INTEGER(iwp) :: k !<
[1000]	330
[1682]	331	REAL(wp) :: alpha_cc !<
	332	REAL(wp) :: autocon !<
	333	REAL(wp) :: dissipation !<
	334	REAL(wp) :: k_au !<
	335	REAL(wp) :: l_mix !<
	336	REAL(wp) :: nu_c !<
	337	REAL(wp) :: phi_au !<
	338	REAL(wp) :: r_cc !<
	339	REAL(wp) :: rc !<
	340	REAL(wp) :: re_lambda !<
	341	REAL(wp) :: sigma_cc !<
	342	REAL(wp) :: tau_cloud !<
	343	REAL(wp) :: xc !<
[1361]	344
	345	CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
	346
[1000]	347	DO i = nxl, nxr
	348	DO j = nys, nyn
[1115]	349	DO k = nzb_s_inner(j,i)+1, nzt
[1000]	350
[1361]	351	IF ( qc(k,j,i) > eps_sb ) THEN
	352
	353	k_au = k_cc / ( 20.0_wp * x0 )
	354	!
	355	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
	356	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
	357	tau_cloud = 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) )
	358	!
	359	!-- Universal function for autoconversion process
	360	!-- (Seifert and Beheng, 2006):
	361	phi_au = 600.0_wp * tau_cloud*0.68_wp &
	362	( 1.0_wp - tau_cloud0.68_wp )3
	363	!
	364	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
	365	!-- (Use constant nu_c = 1.0_wp instead?)
	366	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
	367	!
	368	!-- Mean weight of cloud droplets:
	369	xc = hyrho(k) * qc(k,j,i) / nc_const
	370	!
	371	!-- Parameterized turbulence effects on autoconversion (Seifert,
	372	!-- Nuijens and Stevens, 2010)
	373	IF ( turbulence ) THEN
	374	!
	375	!-- Weight averaged radius of cloud droplets:
	376	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
	377
	378	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
	379	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
	380	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
	381	!
	382	!-- Mixing length (neglecting distance to ground and
	383	!-- stratification)
	384	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
	385	!
	386	!-- Limit dissipation rate according to Seifert, Nuijens and
	387	!-- Stevens (2010)
	388	dissipation = MIN( 0.06_wp, diss(k,j,i) )
	389	!
	390	!-- Compute Taylor-microscale Reynolds number:
	391	re_lambda = 6.0_wp / 11.0_wp * &
	392	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
	393	SQRT( 15.0_wp / kin_vis_air ) * &
	394	dissipation**( 1.0_wp / 6.0_wp )
	395	!
	396	!-- The factor of 1.0E4 is needed to convert the dissipation
	397	!-- rate from m2 s-3 to cm2 s-3.
	398	k_au = k_au * ( 1.0_wp + &
	399	dissipation * 1.0E4_wp * &
	400	( re_lambda * 1.0E-3_wp )*0.25_wp &
	401	( alpha_cc * EXP( -1.0_wp * ( ( rc - &
	402	r_cc ) / &
	403	sigma_cc )**2 &
	404	) + beta_cc &
	405	) &
	406	)
	407	ENDIF
	408	!
	409	!-- Autoconversion rate (Seifert and Beheng, 2006):
	410	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
	411	( nu_c + 1.0_wp )*2 qc(k,j,i)*2 xc*2 &
	412	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
	413	rho_surface
	414	autocon = MIN( autocon, qc(k,j,i) / dt_micro )
	415
	416	qr(k,j,i) = qr(k,j,i) + autocon * dt_micro
	417	qc(k,j,i) = qc(k,j,i) - autocon * dt_micro
	418	nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro
	419
	420	ENDIF
	421
[1000]	422	ENDDO
	423	ENDDO
	424	ENDDO
	425
[1361]	426	CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
	427
[1000]	428	END SUBROUTINE autoconversion
	429
[1106]	430
[1682]	431	!------------------------------------------------------------------------------!
	432	! Description:
	433	! ------------
[1822]	434	!> Autoconversion process (Kessler, 1969).
	435	!------------------------------------------------------------------------------!
	436	SUBROUTINE autoconversion_kessler
	437
	438	USE arrays_3d, &
	439	ONLY: dzw, pt, q, qc
	440
	441	USE cloud_parameters, &
	442	ONLY: l_d_cp, pt_d_t, prec_time_const, prr, ql_crit
	443
	444	USE control_parameters, &
	445	ONLY: dt_micro
	446
	447	USE indices, &
	448	ONLY: nxl, nxr, nyn, nys, nzb_2d, nzt
	449
	450	USE kinds
	451
	452
	453	IMPLICIT NONE
	454
	455	INTEGER(iwp) :: i !<
	456	INTEGER(iwp) :: j !<
	457	INTEGER(iwp) :: k !<
	458
	459	REAL(wp) :: dqdt_precip !<
	460
	461	DO i = nxl, nxr
	462	DO j = nys, nyn
	463	DO k = nzb_2d(j,i)+1, nzt
	464
	465	IF ( qc(k,j,i) > ql_crit ) THEN
	466	dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
	467	ELSE
	468	dqdt_precip = 0.0_wp
	469	ENDIF
	470
	471	qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro
	472	q(k,j,i) = q(k,j,i) - dqdt_precip * dt_micro
	473	pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k)
	474
	475	!
	476	!-- Compute the rain rate
	477	prr(nzb_2d(j,i)+1,j,i) = prr(nzb_2d(j,i)+1,j,i) + dqdt_precip * dzw(k)
	478
	479	ENDDO
	480	ENDDO
	481	ENDDO
	482
	483	END SUBROUTINE autoconversion_kessler
	484
	485
	486	!------------------------------------------------------------------------------!
	487	! Description:
	488	! ------------
[1682]	489	!> Accretion rate (Seifert and Beheng, 2006).
	490	!------------------------------------------------------------------------------!
[1005]	491	SUBROUTINE accretion
[1000]	492
[1361]	493	USE arrays_3d, &
	494	ONLY: diss, qc, qr
	495
	496	USE cloud_parameters, &
	497	ONLY: eps_sb, hyrho, k_cr0
	498
	499	USE control_parameters, &
	500	ONLY: dt_micro, rho_surface, turbulence
	501
	502	USE cpulog, &
	503	ONLY: cpu_log, log_point_s
	504
	505	USE indices, &
[1822]	506	ONLY: nxl, nxr, nys, nyn, nzb_s_inner, nzt
[1361]	507
[1320]	508	USE kinds
[1005]	509
[1000]	510	IMPLICIT NONE
	511
[1682]	512	INTEGER(iwp) :: i !<
	513	INTEGER(iwp) :: j !<
	514	INTEGER(iwp) :: k !<
[1000]	515
[1682]	516	REAL(wp) :: accr !<
	517	REAL(wp) :: k_cr !<
	518	REAL(wp) :: phi_ac !<
	519	REAL(wp) :: tau_cloud !<
[1361]	520
	521	CALL cpu_log( log_point_s(56), 'accretion', 'start' )
	522
[1005]	523	DO i = nxl, nxr
	524	DO j = nys, nyn
[1115]	525	DO k = nzb_s_inner(j,i)+1, nzt
[1000]	526
[1361]	527	IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) ) THEN
	528	!
	529	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
	530	tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
	531	!
	532	!-- Universal function for accretion process (Seifert and
	533	!-- Beheng, 2001):
	534	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
	535	!
	536	!-- Parameterized turbulence effects on autoconversion (Seifert,
	537	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
	538	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
	539	IF ( turbulence ) THEN
	540	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
	541	MIN( 600.0_wp, &
	542	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
	543	)
	544	ELSE
	545	k_cr = k_cr0
	546	ENDIF
	547	!
	548	!-- Accretion rate (Seifert and Beheng, 2006):
	549	accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * &
	550	SQRT( rho_surface * hyrho(k) )
	551	accr = MIN( accr, qc(k,j,i) / dt_micro )
	552
	553	qr(k,j,i) = qr(k,j,i) + accr * dt_micro
	554	qc(k,j,i) = qc(k,j,i) - accr * dt_micro
	555
	556	ENDIF
	557
[1005]	558	ENDDO
	559	ENDDO
[1000]	560	ENDDO
	561
[1361]	562	CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
	563
[1005]	564	END SUBROUTINE accretion
[1000]	565
[1106]	566
[1682]	567	!------------------------------------------------------------------------------!
	568	! Description:
	569	! ------------
	570	!> Collisional breakup rate (Seifert, 2008).
	571	!------------------------------------------------------------------------------!
[1005]	572	SUBROUTINE selfcollection_breakup
[1000]	573
[1361]	574	USE arrays_3d, &
	575	ONLY: nr, qr
	576
	577	USE cloud_parameters, &
	578	ONLY: dpirho_l, eps_sb, hyrho, k_br, k_rr
	579
	580	USE control_parameters, &
	581	ONLY: dt_micro, rho_surface
	582
	583	USE cpulog, &
	584	ONLY: cpu_log, log_point_s
	585
	586	USE indices, &
[1822]	587	ONLY: nxl, nxr, nys, nyn, nzb_s_inner, nzt
[1361]	588
[1320]	589	USE kinds
[1361]	590
[1000]	591	IMPLICIT NONE
	592
[1682]	593	INTEGER(iwp) :: i !<
	594	INTEGER(iwp) :: j !<
	595	INTEGER(iwp) :: k !<
[1000]	596
[1682]	597	REAL(wp) :: breakup !<
	598	REAL(wp) :: dr !<
	599	REAL(wp) :: phi_br !<
	600	REAL(wp) :: selfcoll !<
[1361]	601
	602	CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
	603
[1000]	604	DO i = nxl, nxr
	605	DO j = nys, nyn
[1115]	606	DO k = nzb_s_inner(j,i)+1, nzt
[1361]	607	IF ( qr(k,j,i) > eps_sb ) THEN
	608	!
	609	!-- Selfcollection rate (Seifert and Beheng, 2001):
	610	selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * &
	611	SQRT( hyrho(k) * rho_surface )
	612	!
	613	!-- Weight averaged diameter of rain drops:
	614	dr = ( hyrho(k) * qr(k,j,i) / &
	615	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
	616	!
	617	!-- Collisional breakup rate (Seifert, 2008):
	618	IF ( dr >= 0.3E-3_wp ) THEN
	619	phi_br = k_br * ( dr - 1.1E-3_wp )
	620	breakup = selfcoll * ( phi_br + 1.0_wp )
	621	ELSE
	622	breakup = 0.0_wp
	623	ENDIF
[1000]	624
[1361]	625	selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
	626	nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro
	627
	628	ENDIF
[1000]	629	ENDDO
	630	ENDDO
	631	ENDDO
	632
[1361]	633	CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
	634
[1005]	635	END SUBROUTINE selfcollection_breakup
[1000]	636
[1106]	637
[1682]	638	!------------------------------------------------------------------------------!
	639	! Description:
	640	! ------------
	641	!> Evaporation of precipitable water. Condensation is neglected for
	642	!> precipitable water.
	643	!------------------------------------------------------------------------------!
[1012]	644	SUBROUTINE evaporation_rain
[1000]	645
[1361]	646	USE arrays_3d, &
	647	ONLY: hyp, nr, pt, q, qc, qr
	648
	649	USE cloud_parameters, &
[1822]	650	ONLY: a_term, a_vent, b_term, b_vent, c_evap, c_term, &
	651	diff_coeff_l, dpirho_l, eps_sb, hyrho, kin_vis_air, &
	652	l_d_cp, l_d_r, l_v, r_v, schmidt_p_1d3, &
	653	thermal_conductivity_l, t_d_pt, ventilation_effect
[1361]	654
	655	USE constants, &
	656	ONLY: pi
	657
	658	USE control_parameters, &
	659	ONLY: dt_micro
	660
	661	USE cpulog, &
	662	ONLY: cpu_log, log_point_s
	663
	664	USE indices, &
[1822]	665	ONLY: nxl, nxr, nys, nyn, nzb_s_inner, nzt
[1361]	666
[1320]	667	USE kinds
[1012]	668
	669	IMPLICIT NONE
	670
[1682]	671	INTEGER(iwp) :: i !<
	672	INTEGER(iwp) :: j !<
	673	INTEGER(iwp) :: k !<
[1361]	674
[1682]	675	REAL(wp) :: alpha !<
	676	REAL(wp) :: dr !<
	677	REAL(wp) :: e_s !<
	678	REAL(wp) :: evap !<
	679	REAL(wp) :: evap_nr !<
	680	REAL(wp) :: f_vent !<
	681	REAL(wp) :: g_evap !<
	682	REAL(wp) :: lambda_r !<
	683	REAL(wp) :: mu_r !<
	684	REAL(wp) :: mu_r_2 !<
	685	REAL(wp) :: mu_r_5d2 !<
	686	REAL(wp) :: nr_0 !<
	687	REAL(wp) :: q_s !<
	688	REAL(wp) :: sat !<
	689	REAL(wp) :: t_l !<
	690	REAL(wp) :: temp !<
	691	REAL(wp) :: xr !<
[1361]	692
	693	CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
	694
[1012]	695	DO i = nxl, nxr
	696	DO j = nys, nyn
[1115]	697	DO k = nzb_s_inner(j,i)+1, nzt
[1361]	698	IF ( qr(k,j,i) > eps_sb ) THEN
	699	!
	700	!-- Actual liquid water temperature:
	701	t_l = t_d_pt(k) * pt(k,j,i)
	702	!
	703	!-- Saturation vapor pressure at t_l:
	704	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
	705	( t_l - 35.86_wp ) &
	706	)
	707	!
	708	!-- Computation of saturation humidity:
	709	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
	710	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
	711	q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / &
	712	( 1.0_wp + alpha * q_s )
	713	!
	714	!-- Supersaturation:
	715	sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
	716	!
	717	!-- Evaporation needs only to be calculated in subsaturated regions
	718	IF ( sat < 0.0_wp ) THEN
	719	!
	720	!-- Actual temperature:
	721	temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
	722
	723	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
	724	l_v / ( thermal_conductivity_l * temp ) &
	725	+ r_v * temp / ( diff_coeff_l * e_s ) &
	726	)
	727	!
	728	!-- Mean weight of rain drops
	729	xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
	730	!
	731	!-- Weight averaged diameter of rain drops:
	732	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
	733	!
	734	!-- Compute ventilation factor and intercept parameter
	735	!-- (Seifert and Beheng, 2006; Seifert, 2008):
	736	IF ( ventilation_effect ) THEN
	737	!
	738	!-- Shape parameter of gamma distribution (Milbrandt and Yau,
	739	!-- 2005; Stevens and Seifert, 2008):
	740	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
	741	( dr - 1.4E-3_wp ) ) )
	742	!
	743	!-- Slope parameter of gamma distribution (Seifert, 2008):
	744	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	745	( mu_r + 1.0_wp ) &
	746	)**( 1.0_wp / 3.0_wp ) / dr
[1012]	747
[1361]	748	mu_r_2 = mu_r + 2.0_wp
	749	mu_r_5d2 = mu_r + 2.5_wp
	750
	751	f_vent = a_vent * gamm( mu_r_2 ) * &
	752	lambda_r*( -mu_r_2 ) + b_vent &
	753	schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
	754	gamm( mu_r_5d2 ) * lambda_r*( -mu_r_5d2 ) &
	755	( 1.0_wp - &
	756	0.5_wp * ( b_term / a_term ) * &
	757	( lambda_r / ( c_term + lambda_r ) &
	758	)**mu_r_5d2 - &
	759	0.125_wp * ( b_term / a_term )*2 &
	760	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
	761	)**mu_r_5d2 - &
	762	0.0625_wp * ( b_term / a_term )*3 &
	763	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
	764	)**mu_r_5d2 - &
	765	0.0390625_wp * ( b_term / a_term )*4 &
	766	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
	767	)**mu_r_5d2 &
	768	)
	769
	770	nr_0 = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) / &
	771	gamm( mu_r + 1.0_wp )
	772	ELSE
	773	f_vent = 1.0_wp
	774	nr_0 = nr(k,j,i) * dr
	775	ENDIF
	776	!
	777	!-- Evaporation rate of rain water content (Seifert and
	778	!-- Beheng, 2006):
	779	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / &
	780	hyrho(k)
	781	evap = MAX( evap, -qr(k,j,i) / dt_micro )
	782	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
	783	-nr(k,j,i) / dt_micro )
	784
	785	qr(k,j,i) = qr(k,j,i) + evap * dt_micro
	786	nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro
	787
	788	ENDIF
	789	ENDIF
	790
[1012]	791	ENDDO
	792	ENDDO
	793	ENDDO
	794
[1361]	795	CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
	796
[1012]	797	END SUBROUTINE evaporation_rain
	798
[1106]	799
[1682]	800	!------------------------------------------------------------------------------!
	801	! Description:
	802	! ------------
	803	!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
	804	!------------------------------------------------------------------------------!
[1012]	805	SUBROUTINE sedimentation_cloud
	806
[1361]	807	USE arrays_3d, &
	808	ONLY: ddzu, dzu, pt, q, qc
	809
	810	USE cloud_parameters, &
[1691]	811	ONLY: eps_sb, hyrho, l_d_cp, nc_const, prr, pt_d_t, sed_qc_const
[1361]	812
	813	USE control_parameters, &
[1691]	814	ONLY: call_microphysics_at_all_substeps, dt_micro, &
[1822]	815	intermediate_timestep_count
[1361]	816
	817	USE cpulog, &
	818	ONLY: cpu_log, log_point_s
	819
	820	USE indices, &
	821	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
	822
[1320]	823	USE kinds
[1691]	824
	825	USE statistics, &
	826	ONLY: weight_substep
	827
	828
[1012]	829	IMPLICIT NONE
	830
[1682]	831	INTEGER(iwp) :: i !<
	832	INTEGER(iwp) :: j !<
	833	INTEGER(iwp) :: k !<
[1361]	834
[1682]	835	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
[1361]	836
	837	CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
	838
	839	sed_qc(nzt+1) = 0.0_wp
	840
[1012]	841	DO i = nxl, nxr
	842	DO j = nys, nyn
[1361]	843	DO k = nzt, nzb_s_inner(j,i)+1, -1
[1012]	844
[1361]	845	IF ( qc(k,j,i) > eps_sb ) THEN
	846	sed_qc(k) = sed_qc_const * nc_const*( -2.0_wp / 3.0_wp ) &
	847	( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp )
	848	ELSE
	849	sed_qc(k) = 0.0_wp
	850	ENDIF
	851
	852	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / &
	853	dt_micro + sed_qc(k+1) &
	854	)
	855
	856	q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
	857	ddzu(k+1) / hyrho(k) * dt_micro
	858	qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
	859	ddzu(k+1) / hyrho(k) * dt_micro
	860	pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * &
	861	ddzu(k+1) / hyrho(k) * l_d_cp * &
	862	pt_d_t(k) * dt_micro
	863
[1691]	864	!
	865	!-- Compute the precipitation rate due to cloud (fog) droplets
[1822]	866	IF ( call_microphysics_at_all_substeps ) THEN
	867	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) &
	868	* weight_substep(intermediate_timestep_count)
	869	ELSE
	870	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k)
[1691]	871	ENDIF
	872
[1012]	873	ENDDO
	874	ENDDO
	875	ENDDO
	876
[1361]	877	CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
	878
[1012]	879	END SUBROUTINE sedimentation_cloud
	880
[1106]	881
[1682]	882	!------------------------------------------------------------------------------!
	883	! Description:
	884	! ------------
	885	!> Computation of sedimentation flux. Implementation according to Stevens
	886	!> and Seifert (2008). Code is based on UCLA-LES.
	887	!------------------------------------------------------------------------------!
[1012]	888	SUBROUTINE sedimentation_rain
	889
[1361]	890	USE arrays_3d, &
	891	ONLY: ddzu, dzu, nr, pt, q, qr
	892
	893	USE cloud_parameters, &
[1822]	894	ONLY: a_term, b_term, c_term, dpirho_l, eps_sb, hyrho, &
	895	limiter_sedimentation, l_d_cp, prr, pt_d_t
[1361]	896
	897	USE control_parameters, &
[1691]	898	ONLY: call_microphysics_at_all_substeps, dt_micro, &
	899	intermediate_timestep_count
[1361]	900	USE cpulog, &
	901	ONLY: cpu_log, log_point_s
	902
	903	USE indices, &
	904	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
	905
[1320]	906	USE kinds
[1012]	907
[1361]	908	USE statistics, &
	909	ONLY: weight_substep
	910
[1012]	911	IMPLICIT NONE
	912
[1682]	913	INTEGER(iwp) :: i !<
	914	INTEGER(iwp) :: j !<
	915	INTEGER(iwp) :: k !<
	916	INTEGER(iwp) :: k_run !<
[1361]	917
[1682]	918	REAL(wp) :: c_run !<
	919	REAL(wp) :: d_max !<
	920	REAL(wp) :: d_mean !<
	921	REAL(wp) :: d_min !<
	922	REAL(wp) :: dr !<
	923	REAL(wp) :: flux !<
	924	REAL(wp) :: lambda_r !<
	925	REAL(wp) :: mu_r !<
	926	REAL(wp) :: z_run !<
[1361]	927
[1682]	928	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<
	929	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<
	930	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<
	931	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<
	932	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<
	933	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<
	934	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<
	935	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<
[1361]	936
	937	CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
[1682]	938
[1361]	939	!
	940	!-- Compute velocities
[1012]	941	DO i = nxl, nxr
	942	DO j = nys, nyn
[1115]	943	DO k = nzb_s_inner(j,i)+1, nzt
[1361]	944	IF ( qr(k,j,i) > eps_sb ) THEN
	945	!
	946	!-- Weight averaged diameter of rain drops:
	947	dr = ( hyrho(k) * qr(k,j,i) / &
	948	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
	949	!
	950	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	951	!-- Stevens and Seifert, 2008):
	952	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
	953	( dr - 1.4E-3_wp ) ) )
	954	!
	955	!-- Slope parameter of gamma distribution (Seifert, 2008):
	956	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	957	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
[1012]	958
[1361]	959	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	960	a_term - b_term * ( 1.0_wp + &
	961	c_term / &
	962	lambda_r )*( -1.0_wp &
	963	( mu_r + 1.0_wp ) ) &
	964	) &
	965	)
	966
	967	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	968	a_term - b_term * ( 1.0_wp + &
	969	c_term / &
	970	lambda_r )*( -1.0_wp &
	971	( mu_r + 4.0_wp ) ) &
	972	) &
	973	)
	974	ELSE
	975	w_nr(k) = 0.0_wp
	976	w_qr(k) = 0.0_wp
	977	ENDIF
[1012]	978	ENDDO
[1361]	979	!
	980	!-- Adjust boundary values
	981	w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
	982	w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
	983	w_nr(nzt+1) = 0.0_wp
	984	w_qr(nzt+1) = 0.0_wp
	985	!
	986	!-- Compute Courant number
	987	DO k = nzb_s_inner(j,i)+1, nzt
	988	c_nr(k) = 0.25_wp * ( w_nr(k-1) + &
	989	2.0_wp * w_nr(k) + w_nr(k+1) ) * &
	990	dt_micro * ddzu(k)
	991	c_qr(k) = 0.25_wp * ( w_qr(k-1) + &
	992	2.0_wp * w_qr(k) + w_qr(k+1) ) * &
	993	dt_micro * ddzu(k)
	994	ENDDO
	995	!
	996	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
	997	IF ( limiter_sedimentation ) THEN
	998
	999	DO k = nzb_s_inner(j,i)+1, nzt
[1646]	1000	d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
[1361]	1001	d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
	1002	d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
	1003
	1004	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	1005	2.0_wp * d_max, &
	1006	ABS( d_mean ) )
	1007
[1646]	1008	d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
[1361]	1009	d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
	1010	d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
	1011
	1012	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	1013	2.0_wp * d_max, &
	1014	ABS( d_mean ) )
	1015	ENDDO
	1016
	1017	ELSE
	1018
	1019	nr_slope = 0.0_wp
	1020	qr_slope = 0.0_wp
	1021
	1022	ENDIF
	1023
	1024	sed_nr(nzt+1) = 0.0_wp
	1025	sed_qr(nzt+1) = 0.0_wp
	1026	!
	1027	!-- Compute sedimentation flux
	1028	DO k = nzt, nzb_s_inner(j,i)+1, -1
	1029	!
	1030	!-- Sum up all rain drop number densities which contribute to the flux
	1031	!-- through k-1/2
	1032	flux = 0.0_wp
	1033	z_run = 0.0_wp ! height above z(k)
	1034	k_run = k
	1035	c_run = MIN( 1.0_wp, c_nr(k) )
	1036	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
	1037	flux = flux + hyrho(k_run) * &
	1038	( nr(k_run,j,i) + nr_slope(k_run) * &
	1039	( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)
	1040	z_run = z_run + dzu(k_run)
	1041	k_run = k_run + 1
	1042	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
	1043	ENDDO
	1044	!
	1045	!-- It is not allowed to sediment more rain drop number density than
	1046	!-- available
	1047	flux = MIN( flux, &
	1048	hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * &
	1049	dt_micro &
	1050	)
	1051
	1052	sed_nr(k) = flux / dt_micro
	1053	nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * &
	1054	ddzu(k+1) / hyrho(k) * dt_micro
	1055	!
	1056	!-- Sum up all rain water content which contributes to the flux
	1057	!-- through k-1/2
	1058	flux = 0.0_wp
	1059	z_run = 0.0_wp ! height above z(k)
	1060	k_run = k
	1061	c_run = MIN( 1.0_wp, c_qr(k) )
	1062
	1063	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
	1064
	1065	flux = flux + hyrho(k_run) * ( qr(k_run,j,i) + &
	1066	qr_slope(k_run) * ( 1.0_wp - c_run ) * &
	1067	0.5_wp ) * c_run * dzu(k_run)
	1068	z_run = z_run + dzu(k_run)
	1069	k_run = k_run + 1
	1070	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
	1071
	1072	ENDDO
	1073	!
	1074	!-- It is not allowed to sediment more rain water content than
	1075	!-- available
	1076	flux = MIN( flux, &
	1077	hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * &
	1078	dt_micro &
	1079	)
	1080
	1081	sed_qr(k) = flux / dt_micro
	1082
	1083	qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
	1084	ddzu(k+1) / hyrho(k) * dt_micro
	1085	q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
	1086	ddzu(k+1) / hyrho(k) * dt_micro
	1087	pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * &
	1088	ddzu(k+1) / hyrho(k) * l_d_cp * &
	1089	pt_d_t(k) * dt_micro
	1090	!
	1091	!-- Compute the rain rate
	1092	IF ( call_microphysics_at_all_substeps ) THEN
[1691]	1093	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) &
	1094	* weight_substep(intermediate_timestep_count)
[1361]	1095	ELSE
[1691]	1096	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)
[1361]	1097	ENDIF
	1098
	1099	ENDDO
[1012]	1100	ENDDO
	1101	ENDDO
	1102
[1691]	1103	CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
	1104
	1105	END SUBROUTINE sedimentation_rain
	1106
	1107
	1108	!------------------------------------------------------------------------------!
	1109	! Description:
	1110	! ------------
	1111	!> Computation of the precipitation amount due to gravitational settling of
	1112	!> rain and cloud (fog) droplets
	1113	!------------------------------------------------------------------------------!
	1114	SUBROUTINE calc_precipitation_amount
	1115
	1116	USE cloud_parameters, &
	1117	ONLY: hyrho, precipitation_amount, prr
	1118
	1119	USE control_parameters, &
	1120	ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d, &
	1121	intermediate_timestep_count, intermediate_timestep_count_max,&
	1122	precipitation_amount_interval, time_do2d_xy
	1123
	1124	USE indices, &
	1125	ONLY: nxl, nxr, nys, nyn, nzb_s_inner
	1126
	1127	USE kinds
	1128
	1129	IMPLICIT NONE
	1130
	1131	INTEGER(iwp) :: i !:
	1132	INTEGER(iwp) :: j !:
	1133
	1134
	1135	IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
	1136	( .NOT. call_microphysics_at_all_substeps .OR. &
	1137	intermediate_timestep_count == intermediate_timestep_count_max ) ) &
	1138	THEN
	1139
[1361]	1140	DO i = nxl, nxr
	1141	DO j = nys, nyn
[1691]	1142
[1361]	1143	precipitation_amount(j,i) = precipitation_amount(j,i) + &
	1144	prr(nzb_s_inner(j,i)+1,j,i) * &
	1145	hyrho(nzb_s_inner(j,i)+1) * dt_3d
[1691]	1146
[1361]	1147	ENDDO
	1148	ENDDO
	1149	ENDIF
	1150
[1691]	1151	END SUBROUTINE calc_precipitation_amount
[1361]	1152
[1012]	1153
[1000]	1154	!------------------------------------------------------------------------------!
[1682]	1155	! Description:
	1156	! ------------
	1157	!> Call for grid point i,j
[1000]	1158	!------------------------------------------------------------------------------!
[1022]	1159
[1115]	1160	SUBROUTINE microphysics_control_ij( i, j )
	1161
[1320]	1162	USE arrays_3d, &
[1361]	1163	ONLY: hyp, nc_1d, nr, nr_1d, pt, pt_init, pt_1d, q, q_1d, qc, &
	1164	qc_1d, qr, qr_1d, zu
[1115]	1165
[1320]	1166	USE cloud_parameters, &
[1822]	1167	ONLY: cp, hyrho, nc_const, prr, pt_d_t, r_d, t_d_pt
[1320]	1168
	1169	USE control_parameters, &
[1361]	1170	ONLY: call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, &
	1171	g, intermediate_timestep_count, large_scale_forcing, &
[1822]	1172	lsf_surf, microphysics_seifert, microphysics_kessler, &
	1173	pt_surface, rho_surface, surface_pressure
[1320]	1174
	1175	USE indices, &
	1176	ONLY: nzb, nzt
	1177
	1178	USE kinds
	1179
	1180	USE statistics, &
	1181	ONLY: weight_pres
	1182
[1022]	1183	IMPLICIT NONE
	1184
[1682]	1185	INTEGER(iwp) :: i !<
	1186	INTEGER(iwp) :: j !<
	1187	INTEGER(iwp) :: k !<
[1115]	1188
[1682]	1189	REAL(wp) :: t_surface !<
[1320]	1190
[1361]	1191	IF ( large_scale_forcing .AND. lsf_surf ) THEN
[1241]	1192	!
	1193	!-- Calculate:
	1194	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
	1195	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
	1196	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
[1353]	1197	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
[1241]	1198	DO k = nzb, nzt+1
[1353]	1199	hyp(k) = surface_pressure * 100.0_wp * &
[1361]	1200	( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
[1353]	1201	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
	1202	t_d_pt(k) = 1.0_wp / pt_d_t(k)
[1241]	1203	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
	1204	ENDDO
	1205	!
	1206	!-- Compute reference density
[1353]	1207	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
[1241]	1208	ENDIF
	1209
[1361]	1210	!
	1211	!-- Compute length of time step
	1212	IF ( call_microphysics_at_all_substeps ) THEN
	1213	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
	1214	ELSE
	1215	dt_micro = dt_3d
	1216	ENDIF
[1241]	1217
[1115]	1218	!
[1361]	1219	!-- Use 1d arrays
[1115]	1220	q_1d(:) = q(:,j,i)
	1221	pt_1d(:) = pt(:,j,i)
	1222	qc_1d(:) = qc(:,j,i)
	1223	nc_1d(:) = nc_const
[1822]	1224	IF ( microphysics_seifert ) THEN
[1115]	1225	qr_1d(:) = qr(:,j,i)
	1226	nr_1d(:) = nr(:,j,i)
	1227	ENDIF
[1361]	1228
[1115]	1229	!
[1822]	1230	!-- Reset precipitation rate
	1231	IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp
	1232
	1233	!
[1115]	1234	!-- Compute cloud physics
[1822]	1235	IF( microphysics_kessler ) THEN
	1236
	1237	CALL autoconversion_kessler( i,j )
	1238
	1239	ELSEIF ( microphysics_seifert ) THEN
	1240
	1241	CALL adjust_cloud( i,j )
[1115]	1242	CALL autoconversion( i,j )
	1243	CALL accretion( i,j )
	1244	CALL selfcollection_breakup( i,j )
	1245	CALL evaporation_rain( i,j )
	1246	CALL sedimentation_rain( i,j )
	1247
[1822]	1248	IF ( drizzle ) CALL sedimentation_cloud( i,j )
[1361]	1249
[1691]	1250	ENDIF
	1251
[1822]	1252	CALL calc_precipitation_amount( i,j )
	1253
[1115]	1254	!
[1361]	1255	!-- Store results on the 3d arrays
	1256	q(:,j,i) = q_1d(:)
	1257	pt(:,j,i) = pt_1d(:)
[1822]	1258	IF ( microphysics_seifert ) THEN
[1361]	1259	qr(:,j,i) = qr_1d(:)
	1260	nr(:,j,i) = nr_1d(:)
[1115]	1261	ENDIF
	1262
	1263	END SUBROUTINE microphysics_control_ij
	1264
[1682]	1265	!------------------------------------------------------------------------------!
	1266	! Description:
	1267	! ------------
	1268	!> Adjust number of raindrops to avoid nonlinear effects in
	1269	!> sedimentation and evaporation of rain drops due to too small or
	1270	!> too big weights of rain drops (Stevens and Seifert, 2008).
	1271	!> The same procedure is applied to cloud droplets if they are determined
	1272	!> prognostically. Call for grid point i,j
	1273	!------------------------------------------------------------------------------!
[1115]	1274	SUBROUTINE adjust_cloud_ij( i, j )
	1275
[1320]	1276	USE arrays_3d, &
[1361]	1277	ONLY: qr_1d, nr_1d
[1115]	1278
[1320]	1279	USE cloud_parameters, &
[1822]	1280	ONLY: eps_sb, xrmin, xrmax, hyrho
[1320]	1281
	1282	USE indices, &
[1822]	1283	ONLY: nzb_s_inner, nzt
[1320]	1284
	1285	USE kinds
	1286
[1115]	1287	IMPLICIT NONE
	1288
[1682]	1289	INTEGER(iwp) :: i !<
	1290	INTEGER(iwp) :: j !<
	1291	INTEGER(iwp) :: k !<
	1292
[1115]	1293	DO k = nzb_s_inner(j,i)+1, nzt
[1022]	1294
[1361]	1295	IF ( qr_1d(k) <= eps_sb ) THEN
	1296	qr_1d(k) = 0.0_wp
	1297	nr_1d(k) = 0.0_wp
[1065]	1298	ELSE
[1022]	1299	!
[1048]	1300	!-- Adjust number of raindrops to avoid nonlinear effects in
	1301	!-- sedimentation and evaporation of rain drops due to too small or
[1065]	1302	!-- too big weights of rain drops (Stevens and Seifert, 2008).
[1361]	1303	IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) ) THEN
	1304	nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin
	1305	ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) ) THEN
	1306	nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax
[1048]	1307	ENDIF
[1115]	1308
[1022]	1309	ENDIF
[1115]	1310
[1022]	1311	ENDDO
	1312
[1115]	1313	END SUBROUTINE adjust_cloud_ij
[1022]	1314
[1106]	1315
[1682]	1316	!------------------------------------------------------------------------------!
	1317	! Description:
	1318	! ------------
	1319	!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
	1320	!------------------------------------------------------------------------------!
[1005]	1321	SUBROUTINE autoconversion_ij( i, j )
[1000]	1322
[1320]	1323	USE arrays_3d, &
	1324	ONLY: diss, dzu, nc_1d, nr_1d, qc_1d, qr_1d
[1115]	1325
[1320]	1326	USE cloud_parameters, &
	1327	ONLY: a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3, &
[1822]	1328	c_const, dpirho_l, eps_sb, hyrho, kin_vis_air, k_cc, x0
[1320]	1329
	1330	USE control_parameters, &
	1331	ONLY: dt_micro, rho_surface, turbulence
	1332
	1333	USE grid_variables, &
	1334	ONLY: dx, dy
	1335
	1336	USE indices, &
[1822]	1337	ONLY: nzb_s_inner, nzt
[1320]	1338
	1339	USE kinds
	1340
[1000]	1341	IMPLICIT NONE
	1342
[1682]	1343	INTEGER(iwp) :: i !<
	1344	INTEGER(iwp) :: j !<
	1345	INTEGER(iwp) :: k !<
[1000]	1346
[1682]	1347	REAL(wp) :: alpha_cc !<
	1348	REAL(wp) :: autocon !<
	1349	REAL(wp) :: dissipation !<
	1350	REAL(wp) :: k_au !<
	1351	REAL(wp) :: l_mix !<
	1352	REAL(wp) :: nu_c !<
	1353	REAL(wp) :: phi_au !<
	1354	REAL(wp) :: r_cc !<
	1355	REAL(wp) :: rc !<
	1356	REAL(wp) :: re_lambda !<
	1357	REAL(wp) :: sigma_cc !<
	1358	REAL(wp) :: tau_cloud !<
	1359	REAL(wp) :: xc !<
[1106]	1360
[1115]	1361	DO k = nzb_s_inner(j,i)+1, nzt
[1000]	1362
[1115]	1363	IF ( qc_1d(k) > eps_sb ) THEN
[1361]	1364
	1365	k_au = k_cc / ( 20.0_wp * x0 )
[1012]	1366	!
[1048]	1367	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[1353]	1368	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
	1369	tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
[1012]	1370	!
	1371	!-- Universal function for autoconversion process
	1372	!-- (Seifert and Beheng, 2006):
[1361]	1373	phi_au = 600.0_wp * tau_cloud*0.68_wp ( 1.0_wp - tau_cloud0.68_wp )3
[1012]	1374	!
	1375	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
[1353]	1376	!-- (Use constant nu_c = 1.0_wp instead?)
[1361]	1377	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp )
[1012]	1378	!
	1379	!-- Mean weight of cloud droplets:
[1115]	1380	xc = hyrho(k) * qc_1d(k) / nc_1d(k)
[1012]	1381	!
[1065]	1382	!-- Parameterized turbulence effects on autoconversion (Seifert,
	1383	!-- Nuijens and Stevens, 2010)
	1384	IF ( turbulence ) THEN
	1385	!
	1386	!-- Weight averaged radius of cloud droplets:
[1353]	1387	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
[1065]	1388
[1353]	1389	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
	1390	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
	1391	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
[1065]	1392	!
	1393	!-- Mixing length (neglecting distance to ground and stratification)
[1334]	1394	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
[1065]	1395	!
	1396	!-- Limit dissipation rate according to Seifert, Nuijens and
	1397	!-- Stevens (2010)
[1361]	1398	dissipation = MIN( 0.06_wp, diss(k,j,i) )
[1065]	1399	!
	1400	!-- Compute Taylor-microscale Reynolds number:
[1361]	1401	re_lambda = 6.0_wp / 11.0_wp * &
	1402	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
	1403	SQRT( 15.0_wp / kin_vis_air ) * &
	1404	dissipation**( 1.0_wp / 6.0_wp )
[1065]	1405	!
	1406	!-- The factor of 1.0E4 is needed to convert the dissipation rate
	1407	!-- from m2 s-3 to cm2 s-3.
[1361]	1408	k_au = k_au * ( 1.0_wp + &
	1409	dissipation * 1.0E4_wp * &
	1410	( re_lambda * 1.0E-3_wp )*0.25_wp &
	1411	( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) / &
	1412	sigma_cc )**2 &
	1413	) + beta_cc &
	1414	) &
	1415	)
[1065]	1416	ENDIF
	1417	!
[1012]	1418	!-- Autoconversion rate (Seifert and Beheng, 2006):
[1361]	1419	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
	1420	( nu_c + 1.0_wp )*2 qc_1d(k)*2 xc*2 &
	1421	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
[1115]	1422	rho_surface
	1423	autocon = MIN( autocon, qc_1d(k) / dt_micro )
[1106]	1424
[1115]	1425	qr_1d(k) = qr_1d(k) + autocon * dt_micro
	1426	qc_1d(k) = qc_1d(k) - autocon * dt_micro
	1427	nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
	1428
[1005]	1429	ENDIF
[1000]	1430
	1431	ENDDO
	1432
[1005]	1433	END SUBROUTINE autoconversion_ij
	1434
[1822]	1435	!------------------------------------------------------------------------------!
	1436	! Description:
	1437	! ------------
	1438	!> Autoconversion process (Kessler, 1969).
	1439	!------------------------------------------------------------------------------!
	1440	SUBROUTINE autoconversion_kessler_ij( i, j )
[1106]	1441
[1822]	1442	USE arrays_3d, &
	1443	ONLY: dzw, pt_1d, q_1d, qc_1d
	1444
	1445	USE cloud_parameters, &
	1446	ONLY: l_d_cp, pt_d_t, prec_time_const, prr, ql_crit
	1447
	1448	USE control_parameters, &
	1449	ONLY: dt_micro
	1450
	1451	USE indices, &
	1452	ONLY: nzb_2d, nzt
	1453
	1454	USE kinds
	1455
	1456
	1457	IMPLICIT NONE
	1458
	1459	INTEGER(iwp) :: i !<
	1460	INTEGER(iwp) :: j !<
	1461	INTEGER(iwp) :: k !<
	1462
	1463	REAL(wp) :: dqdt_precip !<
	1464
	1465	DO k = nzb_2d(j,i)+1, nzt
	1466
	1467	IF ( qc_1d(k) > ql_crit ) THEN
	1468	dqdt_precip = prec_time_const * ( qc_1d(k) - ql_crit )
	1469	ELSE
	1470	dqdt_precip = 0.0_wp
	1471	ENDIF
	1472
	1473	qc_1d(k) = qc_1d(k) - dqdt_precip * dt_micro
	1474	q_1d(k) = q_1d(k) - dqdt_precip * dt_micro
	1475	pt_1d(k) = pt_1d(k) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k)
	1476
	1477	!
	1478	!-- Compute the rain rate
	1479	prr(nzb_2d(j,i)+1,j,i) = prr(nzb_2d(j,i)+1,j,i) + &
	1480	dqdt_precip * dzw(k)
	1481
	1482	ENDDO
	1483
	1484	END SUBROUTINE autoconversion_kessler_ij
	1485
[1682]	1486	!------------------------------------------------------------------------------!
	1487	! Description:
	1488	! ------------
	1489	!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
	1490	!------------------------------------------------------------------------------!
[1005]	1491	SUBROUTINE accretion_ij( i, j )
	1492
[1320]	1493	USE arrays_3d, &
	1494	ONLY: diss, qc_1d, qr_1d
[1115]	1495
[1320]	1496	USE cloud_parameters, &
	1497	ONLY: eps_sb, hyrho, k_cr0
	1498
	1499	USE control_parameters, &
	1500	ONLY: dt_micro, rho_surface, turbulence
	1501
	1502	USE indices, &
[1822]	1503	ONLY: nzb_s_inner, nzt
[1320]	1504
	1505	USE kinds
	1506
[1005]	1507	IMPLICIT NONE
	1508
[1682]	1509	INTEGER(iwp) :: i !<
	1510	INTEGER(iwp) :: j !<
	1511	INTEGER(iwp) :: k !<
[1005]	1512
[1682]	1513	REAL(wp) :: accr !<
	1514	REAL(wp) :: k_cr !<
	1515	REAL(wp) :: phi_ac !<
	1516	REAL(wp) :: tau_cloud !<
[1320]	1517
[1115]	1518	DO k = nzb_s_inner(j,i)+1, nzt
	1519	IF ( ( qc_1d(k) > eps_sb ) .AND. ( qr_1d(k) > eps_sb ) ) THEN
[1012]	1520	!
[1048]	1521	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[1353]	1522	tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
[1012]	1523	!
	1524	!-- Universal function for accretion process
[1048]	1525	!-- (Seifert and Beheng, 2001):
[1361]	1526	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
[1012]	1527	!
[1065]	1528	!-- Parameterized turbulence effects on autoconversion (Seifert,
	1529	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
[1361]	1530	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
[1065]	1531	IF ( turbulence ) THEN
[1361]	1532	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
	1533	MIN( 600.0_wp, &
	1534	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
	1535	)
[1065]	1536	ELSE
	1537	k_cr = k_cr0
	1538	ENDIF
	1539	!
[1012]	1540	!-- Accretion rate (Seifert and Beheng, 2006):
[1361]	1541	accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * SQRT( rho_surface * hyrho(k) )
[1115]	1542	accr = MIN( accr, qc_1d(k) / dt_micro )
[1106]	1543
[1115]	1544	qr_1d(k) = qr_1d(k) + accr * dt_micro
	1545	qc_1d(k) = qc_1d(k) - accr * dt_micro
	1546
[1005]	1547	ENDIF
[1106]	1548
[1005]	1549	ENDDO
	1550
[1000]	1551	END SUBROUTINE accretion_ij
	1552
[1005]	1553
[1682]	1554	!------------------------------------------------------------------------------!
	1555	! Description:
	1556	! ------------
	1557	!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
	1558	!------------------------------------------------------------------------------!
[1005]	1559	SUBROUTINE selfcollection_breakup_ij( i, j )
	1560
[1320]	1561	USE arrays_3d, &
	1562	ONLY: nr_1d, qr_1d
	1563
	1564	USE cloud_parameters, &
	1565	ONLY: dpirho_l, eps_sb, hyrho, k_br, k_rr
	1566
	1567	USE control_parameters, &
	1568	ONLY: dt_micro, rho_surface
	1569
	1570	USE indices, &
[1822]	1571	ONLY: nzb_s_inner, nzt
[1320]	1572
	1573	USE kinds
[1005]	1574
	1575	IMPLICIT NONE
	1576
[1682]	1577	INTEGER(iwp) :: i !<
	1578	INTEGER(iwp) :: j !<
	1579	INTEGER(iwp) :: k !<
[1005]	1580
[1682]	1581	REAL(wp) :: breakup !<
	1582	REAL(wp) :: dr !<
	1583	REAL(wp) :: phi_br !<
	1584	REAL(wp) :: selfcoll !<
[1320]	1585
[1115]	1586	DO k = nzb_s_inner(j,i)+1, nzt
	1587	IF ( qr_1d(k) > eps_sb ) THEN
[1012]	1588	!
[1115]	1589	!-- Selfcollection rate (Seifert and Beheng, 2001):
[1361]	1590	selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
[1012]	1591	!
[1115]	1592	!-- Weight averaged diameter of rain drops:
[1334]	1593	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	1594	!
[1048]	1595	!-- Collisional breakup rate (Seifert, 2008):
[1353]	1596	IF ( dr >= 0.3E-3_wp ) THEN
	1597	phi_br = k_br * ( dr - 1.1E-3_wp )
	1598	breakup = selfcoll * ( phi_br + 1.0_wp )
[1005]	1599	ELSE
[1353]	1600	breakup = 0.0_wp
[1005]	1601	ENDIF
[1048]	1602
[1115]	1603	selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
	1604	nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
[1106]	1605
[1005]	1606	ENDIF
	1607	ENDDO
	1608
	1609	END SUBROUTINE selfcollection_breakup_ij
	1610
[1106]	1611
[1682]	1612	!------------------------------------------------------------------------------!
	1613	! Description:
	1614	! ------------
	1615	!> Evaporation of precipitable water. Condensation is neglected for
	1616	!> precipitable water. Call for grid point i,j
	1617	!------------------------------------------------------------------------------!
[1012]	1618	SUBROUTINE evaporation_rain_ij( i, j )
	1619
[1320]	1620	USE arrays_3d, &
	1621	ONLY: hyp, nr_1d, pt_1d, q_1d, qc_1d, qr_1d
[1048]	1622
[1320]	1623	USE cloud_parameters, &
	1624	ONLY: a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,&
[1822]	1625	dpirho_l, eps_sb, hyrho, kin_vis_air, l_d_cp, l_d_r, &
	1626	l_v, r_v, schmidt_p_1d3, thermal_conductivity_l, &
[1320]	1627	t_d_pt, ventilation_effect
	1628
	1629	USE constants, &
	1630	ONLY: pi
	1631
	1632	USE control_parameters, &
	1633	ONLY: dt_micro
	1634
	1635	USE indices, &
[1822]	1636	ONLY: nzb_s_inner, nzt
[1320]	1637
	1638	USE kinds
	1639
[1012]	1640	IMPLICIT NONE
	1641
[1682]	1642	INTEGER(iwp) :: i !<
	1643	INTEGER(iwp) :: j !<
	1644	INTEGER(iwp) :: k !<
[1012]	1645
[1682]	1646	REAL(wp) :: alpha !<
	1647	REAL(wp) :: dr !<
	1648	REAL(wp) :: e_s !<
	1649	REAL(wp) :: evap !<
	1650	REAL(wp) :: evap_nr !<
	1651	REAL(wp) :: f_vent !<
	1652	REAL(wp) :: g_evap !<
	1653	REAL(wp) :: lambda_r !<
	1654	REAL(wp) :: mu_r !<
	1655	REAL(wp) :: mu_r_2 !<
	1656	REAL(wp) :: mu_r_5d2 !<
	1657	REAL(wp) :: nr_0 !<
	1658	REAL(wp) :: q_s !<
	1659	REAL(wp) :: sat !<
	1660	REAL(wp) :: t_l !<
	1661	REAL(wp) :: temp !<
	1662	REAL(wp) :: xr !<
[1320]	1663
[1115]	1664	DO k = nzb_s_inner(j,i)+1, nzt
	1665	IF ( qr_1d(k) > eps_sb ) THEN
[1012]	1666	!
	1667	!-- Actual liquid water temperature:
[1115]	1668	t_l = t_d_pt(k) * pt_1d(k)
[1012]	1669	!
	1670	!-- Saturation vapor pressure at t_l:
[1361]	1671	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
	1672	( t_l - 35.86_wp ) &
	1673	)
[1012]	1674	!
	1675	!-- Computation of saturation humidity:
[1361]	1676	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
[1353]	1677	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
[1361]	1678	q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
[1012]	1679	!
[1106]	1680	!-- Supersaturation:
[1361]	1681	sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
[1012]	1682	!
[1361]	1683	!-- Evaporation needs only to be calculated in subsaturated regions
	1684	IF ( sat < 0.0_wp ) THEN
[1012]	1685	!
[1361]	1686	!-- Actual temperature:
	1687	temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
	1688
	1689	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v / &
	1690	( thermal_conductivity_l * temp ) + &
	1691	r_v * temp / ( diff_coeff_l * e_s ) &
	1692	)
[1012]	1693	!
[1361]	1694	!-- Mean weight of rain drops
	1695	xr = hyrho(k) * qr_1d(k) / nr_1d(k)
[1115]	1696	!
[1361]	1697	!-- Weight averaged diameter of rain drops:
	1698	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	1699	!
[1361]	1700	!-- Compute ventilation factor and intercept parameter
	1701	!-- (Seifert and Beheng, 2006; Seifert, 2008):
	1702	IF ( ventilation_effect ) THEN
[1115]	1703	!
[1361]	1704	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	1705	!-- Stevens and Seifert, 2008):
	1706	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
	1707	!
	1708	!-- Slope parameter of gamma distribution (Seifert, 2008):
	1709	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	1710	( mu_r + 1.0_wp ) &
	1711	)**( 1.0_wp / 3.0_wp ) / dr
[1115]	1712
[1361]	1713	mu_r_2 = mu_r + 2.0_wp
	1714	mu_r_5d2 = mu_r + 2.5_wp
	1715
	1716	f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) + &
	1717	b_vent * schmidt_p_1d3 * &
	1718	SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * &
	1719	lambda_r*( -mu_r_5d2 ) &
	1720	( 1.0_wp - &
	1721	0.5_wp * ( b_term / a_term ) * &
	1722	( lambda_r / ( c_term + lambda_r ) &
	1723	)**mu_r_5d2 - &
	1724	0.125_wp * ( b_term / a_term )*2 &
	1725	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
	1726	)**mu_r_5d2 - &
	1727	0.0625_wp * ( b_term / a_term )*3 &
	1728	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
	1729	)**mu_r_5d2 - &
	1730	0.0390625_wp * ( b_term / a_term )*4 &
	1731	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
	1732	)**mu_r_5d2 &
	1733	)
	1734
	1735	nr_0 = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) / &
	1736	gamm( mu_r + 1.0_wp )
	1737	ELSE
	1738	f_vent = 1.0_wp
	1739	nr_0 = nr_1d(k) * dr
	1740	ENDIF
[1012]	1741	!
[1361]	1742	!-- Evaporation rate of rain water content (Seifert and Beheng, 2006):
	1743	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
	1744	evap = MAX( evap, -qr_1d(k) / dt_micro )
	1745	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
	1746	-nr_1d(k) / dt_micro )
[1106]	1747
[1361]	1748	qr_1d(k) = qr_1d(k) + evap * dt_micro
	1749	nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
[1115]	1750
[1361]	1751	ENDIF
[1012]	1752	ENDIF
[1106]	1753
[1012]	1754	ENDDO
	1755
	1756	END SUBROUTINE evaporation_rain_ij
	1757
[1106]	1758
[1682]	1759	!------------------------------------------------------------------------------!
	1760	! Description:
	1761	! ------------
	1762	!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
	1763	!> Call for grid point i,j
	1764	!------------------------------------------------------------------------------!
[1012]	1765	SUBROUTINE sedimentation_cloud_ij( i, j )
	1766
[1320]	1767	USE arrays_3d, &
	1768	ONLY: ddzu, dzu, nc_1d, pt_1d, q_1d, qc_1d
	1769
	1770	USE cloud_parameters, &
[1691]	1771	ONLY: eps_sb, hyrho, l_d_cp, prr, pt_d_t, sed_qc_const
[1320]	1772
	1773	USE control_parameters, &
[1822]	1774	ONLY: call_microphysics_at_all_substeps, dt_micro, &
	1775	intermediate_timestep_count
[1320]	1776
	1777	USE indices, &
	1778	ONLY: nzb, nzb_s_inner, nzt
	1779
	1780	USE kinds
[1012]	1781
[1691]	1782	USE statistics, &
	1783	ONLY: weight_substep
	1784
[1012]	1785	IMPLICIT NONE
	1786
[1682]	1787	INTEGER(iwp) :: i !<
	1788	INTEGER(iwp) :: j !<
	1789	INTEGER(iwp) :: k !<
[1106]	1790
[1682]	1791	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
[1115]	1792
[1353]	1793	sed_qc(nzt+1) = 0.0_wp
[1012]	1794
[1115]	1795	DO k = nzt, nzb_s_inner(j,i)+1, -1
	1796	IF ( qc_1d(k) > eps_sb ) THEN
[1361]	1797	sed_qc(k) = sed_qc_const * nc_1d(k)*( -2.0_wp / 3.0_wp ) &
	1798	( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )
[1115]	1799	ELSE
[1353]	1800	sed_qc(k) = 0.0_wp
[1012]	1801	ENDIF
[1115]	1802
[1361]	1803	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) / &
	1804	dt_micro + sed_qc(k+1) &
	1805	)
[1115]	1806
[1361]	1807	q_1d(k) = q_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
[1115]	1808	hyrho(k) * dt_micro
[1361]	1809	qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
[1115]	1810	hyrho(k) * dt_micro
[1361]	1811	pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
[1115]	1812	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
	1813
[1691]	1814	!
	1815	!-- Compute the precipitation rate of cloud (fog) droplets
[1822]	1816	IF ( call_microphysics_at_all_substeps ) THEN
	1817	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * &
[1691]	1818	weight_substep(intermediate_timestep_count)
[1822]	1819	ELSE
	1820	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k)
[1691]	1821	ENDIF
	1822
[1012]	1823	ENDDO
	1824
	1825	END SUBROUTINE sedimentation_cloud_ij
	1826
[1106]	1827
[1682]	1828	!------------------------------------------------------------------------------!
	1829	! Description:
	1830	! ------------
	1831	!> Computation of sedimentation flux. Implementation according to Stevens
	1832	!> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j
	1833	!------------------------------------------------------------------------------!
[1012]	1834	SUBROUTINE sedimentation_rain_ij( i, j )
	1835
[1320]	1836	USE arrays_3d, &
	1837	ONLY: ddzu, dzu, nr_1d, pt_1d, q_1d, qr_1d
	1838
	1839	USE cloud_parameters, &
[1822]	1840	ONLY: a_term, b_term, c_term, dpirho_l, eps_sb, hyrho, &
	1841	limiter_sedimentation, l_d_cp, prr, pt_d_t
[1320]	1842
	1843	USE control_parameters, &
[1822]	1844	ONLY: call_microphysics_at_all_substeps, dt_micro, &
	1845	intermediate_timestep_count
[1320]	1846
	1847	USE indices, &
	1848	ONLY: nzb, nzb_s_inner, nzt
	1849
	1850	USE kinds
	1851
	1852	USE statistics, &
	1853	ONLY: weight_substep
[1012]	1854
	1855	IMPLICIT NONE
	1856
[1682]	1857	INTEGER(iwp) :: i !<
	1858	INTEGER(iwp) :: j !<
	1859	INTEGER(iwp) :: k !<
	1860	INTEGER(iwp) :: k_run !<
[1012]	1861
[1682]	1862	REAL(wp) :: c_run !<
	1863	REAL(wp) :: d_max !<
	1864	REAL(wp) :: d_mean !<
	1865	REAL(wp) :: d_min !<
	1866	REAL(wp) :: dr !<
	1867	REAL(wp) :: flux !<
	1868	REAL(wp) :: lambda_r !<
	1869	REAL(wp) :: mu_r !<
	1870	REAL(wp) :: z_run !<
[1320]	1871
[1682]	1872	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<
	1873	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<
	1874	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<
	1875	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<
	1876	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<
	1877	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<
	1878	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<
	1879	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<
[1320]	1880
[1012]	1881	!
[1065]	1882	!-- Compute velocities
	1883	DO k = nzb_s_inner(j,i)+1, nzt
[1115]	1884	IF ( qr_1d(k) > eps_sb ) THEN
	1885	!
	1886	!-- Weight averaged diameter of rain drops:
[1334]	1887	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	1888	!
	1889	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	1890	!-- Stevens and Seifert, 2008):
[1353]	1891	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
[1115]	1892	!
	1893	!-- Slope parameter of gamma distribution (Seifert, 2008):
[1361]	1894	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	1895	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
[1115]	1896
[1361]	1897	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	1898	a_term - b_term * ( 1.0_wp + &
	1899	c_term / lambda_r )*( -1.0_wp &
	1900	( mu_r + 1.0_wp ) ) &
	1901	) &
	1902	)
	1903	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	1904	a_term - b_term * ( 1.0_wp + &
	1905	c_term / lambda_r )*( -1.0_wp &
	1906	( mu_r + 4.0_wp ) ) &
	1907	) &
	1908	)
[1065]	1909	ELSE
[1353]	1910	w_nr(k) = 0.0_wp
	1911	w_qr(k) = 0.0_wp
[1065]	1912	ENDIF
	1913	ENDDO
[1048]	1914	!
[1065]	1915	!-- Adjust boundary values
[1115]	1916	w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
	1917	w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
[1353]	1918	w_nr(nzt+1) = 0.0_wp
	1919	w_qr(nzt+1) = 0.0_wp
[1065]	1920	!
	1921	!-- Compute Courant number
[1115]	1922	DO k = nzb_s_inner(j,i)+1, nzt
[1361]	1923	c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) * &
[1115]	1924	dt_micro * ddzu(k)
[1361]	1925	c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * &
[1115]	1926	dt_micro * ddzu(k)
	1927	ENDDO
[1065]	1928	!
	1929	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
	1930	IF ( limiter_sedimentation ) THEN
	1931
[1115]	1932	DO k = nzb_s_inner(j,i)+1, nzt
[1646]	1933	d_mean = 0.5_wp * ( qr_1d(k+1) - qr_1d(k-1) )
[1115]	1934	d_min = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
	1935	d_max = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
[1065]	1936
[1361]	1937	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	1938	2.0_wp * d_max, &
	1939	ABS( d_mean ) )
[1065]	1940
[1646]	1941	d_mean = 0.5_wp * ( nr_1d(k+1) - nr_1d(k-1) )
[1115]	1942	d_min = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
	1943	d_max = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
[1065]	1944
[1361]	1945	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	1946	2.0_wp * d_max, &
	1947	ABS( d_mean ) )
[1022]	1948	ENDDO
[1048]	1949
[1065]	1950	ELSE
[1106]	1951
[1353]	1952	nr_slope = 0.0_wp
	1953	qr_slope = 0.0_wp
[1106]	1954
[1065]	1955	ENDIF
[1115]	1956
[1353]	1957	sed_nr(nzt+1) = 0.0_wp
	1958	sed_qr(nzt+1) = 0.0_wp
[1065]	1959	!
	1960	!-- Compute sedimentation flux
[1115]	1961	DO k = nzt, nzb_s_inner(j,i)+1, -1
[1065]	1962	!
	1963	!-- Sum up all rain drop number densities which contribute to the flux
	1964	!-- through k-1/2
[1353]	1965	flux = 0.0_wp
	1966	z_run = 0.0_wp ! height above z(k)
[1065]	1967	k_run = k
[1346]	1968	c_run = MIN( 1.0_wp, c_nr(k) )
[1353]	1969	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
[1361]	1970	flux = flux + hyrho(k_run) * &
	1971	( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) * &
[1353]	1972	0.5_wp ) * c_run * dzu(k_run)
[1065]	1973	z_run = z_run + dzu(k_run)
	1974	k_run = k_run + 1
[1346]	1975	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
[1022]	1976	ENDDO
	1977	!
[1065]	1978	!-- It is not allowed to sediment more rain drop number density than
	1979	!-- available
[1361]	1980	flux = MIN( flux, &
[1115]	1981	hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
[1065]	1982
[1115]	1983	sed_nr(k) = flux / dt_micro
[1361]	1984	nr_1d(k) = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / &
	1985	hyrho(k) * dt_micro
[1065]	1986	!
	1987	!-- Sum up all rain water content which contributes to the flux
	1988	!-- through k-1/2
[1353]	1989	flux = 0.0_wp
	1990	z_run = 0.0_wp ! height above z(k)
[1065]	1991	k_run = k
[1346]	1992	c_run = MIN( 1.0_wp, c_qr(k) )
[1106]	1993
[1361]	1994	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
[1106]	1995
[1361]	1996	flux = flux + hyrho(k_run) * &
	1997	( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) * &
[1353]	1998	0.5_wp ) * c_run * dzu(k_run)
[1065]	1999	z_run = z_run + dzu(k_run)
	2000	k_run = k_run + 1
[1346]	2001	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
[1106]	2002
[1065]	2003	ENDDO
	2004	!
	2005	!-- It is not allowed to sediment more rain water content than available
[1361]	2006	flux = MIN( flux, &
[1115]	2007	hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
[1065]	2008
[1115]	2009	sed_qr(k) = flux / dt_micro
	2010
[1361]	2011	qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
[1115]	2012	hyrho(k) * dt_micro
[1361]	2013	q_1d(k) = q_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
[1115]	2014	hyrho(k) * dt_micro
[1361]	2015	pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
[1115]	2016	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
[1065]	2017	!
	2018	!-- Compute the rain rate
[1361]	2019	IF ( call_microphysics_at_all_substeps ) THEN
[1691]	2020	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) &
	2021	* weight_substep(intermediate_timestep_count)
[1361]	2022	ELSE
[1691]	2023	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)
[1361]	2024	ENDIF
	2025
[1065]	2026	ENDDO
[1115]	2027
[1691]	2028	END SUBROUTINE sedimentation_rain_ij
[1012]	2029
[1691]	2030
	2031	!------------------------------------------------------------------------------!
	2032	! Description:
	2033	! ------------
	2034	!> This subroutine computes the precipitation amount due to gravitational
	2035	!> settling of rain and cloud (fog) droplets
	2036	!------------------------------------------------------------------------------!
	2037	SUBROUTINE calc_precipitation_amount_ij( i, j )
	2038
	2039	USE cloud_parameters, &
	2040	ONLY: hyrho, precipitation_amount, prr
	2041
	2042	USE control_parameters, &
	2043	ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d, &
	2044	intermediate_timestep_count, intermediate_timestep_count_max,&
[1822]	2045	precipitation_amount_interval, time_do2d_xy
[1691]	2046
	2047	USE indices, &
	2048	ONLY: nzb_s_inner
	2049
	2050	USE kinds
	2051
	2052	IMPLICIT NONE
	2053
	2054	INTEGER(iwp) :: i !:
	2055	INTEGER(iwp) :: j !:
	2056
	2057
	2058	IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
	2059	( .NOT. call_microphysics_at_all_substeps .OR. &
	2060	intermediate_timestep_count == intermediate_timestep_count_max ) ) &
	2061	THEN
	2062
[1361]	2063	precipitation_amount(j,i) = precipitation_amount(j,i) + &
	2064	prr(nzb_s_inner(j,i)+1,j,i) * &
[1115]	2065	hyrho(nzb_s_inner(j,i)+1) * dt_3d
[1048]	2066	ENDIF
	2067
[1691]	2068	END SUBROUTINE calc_precipitation_amount_ij
[1012]	2069
[1361]	2070	!------------------------------------------------------------------------------!
[1682]	2071	! Description:
	2072	! ------------
	2073	!> This function computes the gamma function (Press et al., 1992).
	2074	!> The gamma function is needed for the calculation of the evaporation
	2075	!> of rain drops.
[1361]	2076	!------------------------------------------------------------------------------!
[1012]	2077	FUNCTION gamm( xx )
[1048]	2078
[1320]	2079	USE cloud_parameters, &
	2080	ONLY: cof, stp
	2081
	2082	USE kinds
	2083
[1012]	2084	IMPLICIT NONE
[1106]	2085
[1682]	2086	INTEGER(iwp) :: j !<
[1320]	2087
[1682]	2088	REAL(wp) :: gamm !<
	2089	REAL(wp) :: ser !<
	2090	REAL(wp) :: tmp !<
	2091	REAL(wp) :: x_gamm !<
	2092	REAL(wp) :: xx !<
	2093	REAL(wp) :: y_gamm !<
[1320]	2094
[1012]	2095	x_gamm = xx
	2096	y_gamm = x_gamm
[1353]	2097	tmp = x_gamm + 5.5_wp
	2098	tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
[1334]	2099	ser = 1.000000000190015_wp
[1106]	2100
	2101	DO j = 1, 6
[1353]	2102	y_gamm = y_gamm + 1.0_wp
[1012]	2103	ser = ser + cof( j ) / y_gamm
[1106]	2104	ENDDO
	2105
[1012]	2106	!
	2107	!-- Until this point the algorithm computes the logarithm of the gamma
	2108	!-- function. Hence, the exponential function is used.
	2109	! gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
	2110	gamm = EXP( tmp ) * stp * ser / x_gamm
[1106]	2111
[1012]	2112	RETURN
	2113
	2114	END FUNCTION gamm
	2115
	2116	END MODULE microphysics_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |