Home

Context Navigation

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

Last change on this file since 1822 was 1822, checked in by hoffmann, 8 years ago
changes in LPM and bulk cloud microphysics
Property svn:keywords set to `Id`
File size: 81.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |