Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |