Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |