Home

Context Navigation

source: palm/trunk/SOURCE/microphysics_mod.f90 @ 2106

Last change on this file since 2106 was 2101, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 86.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |