Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |