Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |