Home

Context Navigation

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

Last change on this file since 3159 was 3026, checked in by schwenkel, 6 years ago
Changed the name specific humidity to mixing ratio
Property svn:keywords set to `Id`
File size: 123.1 KB

Rev	Line
[1850]	1	!> @file microphysics_mod.f90
[2000]	2	!------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1093]	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	!
[2718]	17	! Copyright 1997-2018 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1093]	19	!
[1000]	20	! Current revisions:
[1092]	21	! ------------------
[2318]	22	!
[2701]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: microphysics_mod.f90 3026 2018-05-22 10:30:53Z sward $
[3026]	27	! Changed the name specific humidity to mixing ratio, since we are computing
	28	! mixing ratios.
	29	!
	30	! 2718 2018-01-02 08:49:38Z maronga
[2716]	31	! Corrected "Former revisions" section
[2701]	32	!
[2716]	33	! 2701 2017-12-15 15:40:50Z suehring
	34	! Changes from last commit documented
	35	!
[2701]	36	! 2698 2017-12-14 18:46:24Z suehring
[2716]	37	! Bugfix in get_topography_top_index
	38	!
	39	! 2696 2017-12-14 17:12:51Z kanani
	40	! Change in file header (GPL part)
	41	!
	42	! 2608 2017-11-13 14:04:26Z schwenkel
[2608]	43	! Calculation of supersaturation in external module (diagnostic_quantities_mod).
	44	! Change: correct calculation of saturation specific humidity to saturation
	45	! mixing ratio (the factor of 0.378 vanishes).
	46	!
	47	! 2522 2017-10-05 14:20:37Z schwenkel
[2522]	48	! Minor bugfix
	49	!
	50	! 2375 2017-08-29 14:10:28Z schwenkel
[2375]	51	! Improved aerosol initilization and some minor bugfixes
	52	! for droplet sedimenation
	53	!
	54	! 2318 2017-07-20 17:27:44Z suehring
[2318]	55	! Get topography top index via Function call
	56	!
	57	! 2317 2017-07-20 17:27:19Z suehring
[2312]	58	! s1 changed to log_sigma
	59	!
	60	! 2292 2017-06-20 09:51:42Z schwenkel
	61	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
	62	! includes two more prognostic equations for cloud drop concentration (nc)
	63	! and cloud water content (qc).
	64	! - The process of activation is parameterized with a simple Twomey
	65	! activion scheme or with considering solution and curvature
[2292]	66	! effects (Khvorostyanov and Curry ,2006).
	67	! - The saturation adjustment scheme is replaced by the parameterization
	68	! of condensation rates (Khairoutdinov and Kogan, 2000, Mon. Wea. Rev.,128).
	69	! - All other microphysical processes of Seifert and Beheng are used.
[2312]	70	! Additionally, in those processes the reduction of cloud number concentration
	71	! is considered.
	72	!
[2292]	73	! 2233 2017-05-30 18:08:54Z suehring
[1321]	74	!
[2233]	75	! 2232 2017-05-30 17:47:52Z suehring
	76	! Adjustments to new topography and surface concept
[2312]	77	!
[2156]	78	! 2155 2017-02-21 09:57:40Z hoffmann
	79	! Bugfix in the calculation of microphysical quantities on ghost points.
	80	!
[2032]	81	! 2031 2016-10-21 15:11:58Z knoop
	82	! renamed variable rho to rho_ocean
[2155]	83	!
[2001]	84	! 2000 2016-08-20 18:09:15Z knoop
	85	! Forced header and separation lines into 80 columns
[2155]	86	!
[1851]	87	! 1850 2016-04-08 13:29:27Z maronga
	88	! Module renamed
	89	! Adapted for modularization of microphysics.
	90	!
[1846]	91	! 1845 2016-04-08 08:29:13Z raasch
	92	! nzb_2d replaced by nzb_s_inner, Kessler precipitation is stored at surface
	93	! point (instead of one point above surface)
	94	!
[1832]	95	! 1831 2016-04-07 13:15:51Z hoffmann
	96	! turbulence renamed collision_turbulence,
	97	! drizzle renamed cloud_water_sedimentation. cloud_water_sedimentation also
	98	! avaialble for microphysics_kessler.
	99	!
[1823]	100	! 1822 2016-04-07 07:49:42Z hoffmann
	101	! Unused variables removed.
	102	! Kessler scheme integrated.
	103	!
[1692]	104	! 1691 2015-10-26 16:17:44Z maronga
	105	! Added new routine calc_precipitation_amount. The routine now allows to account
	106	! for precipitation due to sedimenation of cloud (fog) droplets
[2155]	107	!
[1683]	108	! 1682 2015-10-07 23:56:08Z knoop
[2155]	109	! Code annotations made doxygen readable
[1683]	110	!
[1647]	111	! 1646 2015-09-02 16:00:10Z hoffmann
	112	! Bugfix: Wrong computation of d_mean.
	113	!
[1362]	114	! 1361 2014-04-16 15:17:48Z hoffmann
	115	! Bugfix in sedimentation_rain: Index corrected.
	116	! Vectorized version of adjust_cloud added.
	117	! Little reformatting of the code.
[2155]	118	!
[1354]	119	! 1353 2014-04-08 15:21:23Z heinze
	120	! REAL constants provided with KIND-attribute
[2155]	121	!
[1347]	122	! 1346 2014-03-27 13:18:20Z heinze
[2155]	123	! Bugfix: REAL constants provided with KIND-attribute especially in call of
[1347]	124	! intrinsic function like MAX, MIN, SIGN
[2155]	125	!
[1335]	126	! 1334 2014-03-25 12:21:40Z heinze
	127	! Bugfix: REAL constants provided with KIND-attribute
	128	!
[1323]	129	! 1322 2014-03-20 16:38:49Z raasch
	130	! REAL constants defined as wp-kind
	131	!
[1321]	132	! 1320 2014-03-20 08:40:49Z raasch
[1320]	133	! ONLY-attribute added to USE-statements,
	134	! kind-parameters added to all INTEGER and REAL declaration statements,
	135	! kinds are defined in new module kinds,
	136	! comment fields (!:) to be used for variable explanations added to
	137	! all variable declaration statements
[1000]	138	!
[1242]	139	! 1241 2013-10-30 11:36:58Z heinze
[2031]	140	! hyp and rho_ocean have to be calculated at each time step if data from external
[1242]	141	! file LSF_DATA are used
	142	!
[1116]	143	! 1115 2013-03-26 18:16:16Z hoffmann
	144	! microphyical tendencies are calculated in microphysics_control in an optimized
	145	! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
	146	!
[1107]	147	! 1106 2013-03-04 05:31:38Z raasch
	148	! small changes in code formatting
	149	!
[1093]	150	! 1092 2013-02-02 11:24:22Z raasch
	151	! unused variables removed
	152	! file put under GPL
	153	!
[1066]	154	! 1065 2012-11-22 17:42:36Z hoffmann
	155	! Sedimentation process implemented according to Stevens and Seifert (2008).
[1115]	156	! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
[1066]	157	! and Stevens, 2010).
	158	!
[1054]	159	! 1053 2012-11-13 17:11:03Z hoffmann
	160	! initial revision
[2155]	161	!
[1000]	162	! Description:
	163	! ------------
[1849]	164	!> Calculate bilk cloud microphysics.
[1000]	165	!------------------------------------------------------------------------------!
[1682]	166	MODULE microphysics_mod
[1000]	167
[1849]	168	USE kinds
	169
	170	IMPLICIT NONE
	171
[2292]	172	LOGICAL :: cloud_water_sedimentation = .FALSE. !< cloud water sedimentation
	173	LOGICAL :: curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory
	174	LOGICAL :: limiter_sedimentation = .TRUE. !< sedimentation limiter
	175	LOGICAL :: collision_turbulence = .FALSE. !< turbulence effects
	176	LOGICAL :: ventilation_effect = .TRUE. !< ventilation effect
[1849]	177
	178	REAL(wp) :: a_1 = 8.69E-4_wp !< coef. in turb. parametrization (cm-2 s3)
	179	REAL(wp) :: a_2 = -7.38E-5_wp !< coef. in turb. parametrization (cm-2 s3)
	180	REAL(wp) :: a_3 = -1.40E-2_wp !< coef. in turb. parametrization
	181	REAL(wp) :: a_term = 9.65_wp !< coef. for terminal velocity (m s-1)
	182	REAL(wp) :: a_vent = 0.78_wp !< coef. for ventilation effect
	183	REAL(wp) :: b_1 = 11.45E-6_wp !< coef. in turb. parametrization (m)
	184	REAL(wp) :: b_2 = 9.68E-6_wp !< coef. in turb. parametrization (m)
	185	REAL(wp) :: b_3 = 0.62_wp !< coef. in turb. parametrization
	186	REAL(wp) :: b_term = 9.8_wp !< coef. for terminal velocity (m s-1)
	187	REAL(wp) :: b_vent = 0.308_wp !< coef. for ventilation effect
	188	REAL(wp) :: beta_cc = 3.09E-4_wp !< coef. in turb. parametrization (cm-2 s3)
	189	REAL(wp) :: c_1 = 4.82E-6_wp !< coef. in turb. parametrization (m)
	190	REAL(wp) :: c_2 = 4.8E-6_wp !< coef. in turb. parametrization (m)
	191	REAL(wp) :: c_3 = 0.76_wp !< coef. in turb. parametrization
	192	REAL(wp) :: c_const = 0.93_wp !< const. in Taylor-microscale Reynolds number
	193	REAL(wp) :: c_evap = 0.7_wp !< constant in evaporation
	194	REAL(wp) :: c_term = 600.0_wp !< coef. for terminal velocity (m-1)
	195	REAL(wp) :: diff_coeff_l = 0.23E-4_wp !< diffusivity of water vapor (m2 s-1)
[2155]	196	REAL(wp) :: eps_sb = 1.0E-10_wp !< threshold in two-moments scheme
[2292]	197	REAL(wp) :: eps_mr = 0.0_wp !< threshold for morrison scheme
[1849]	198	REAL(wp) :: k_cc = 9.44E09_wp !< const. cloud-cloud kernel (m3 kg-2 s-1)
	199	REAL(wp) :: k_cr0 = 4.33_wp !< const. cloud-rain kernel (m3 kg-1 s-1)
	200	REAL(wp) :: k_rr = 7.12_wp !< const. rain-rain kernel (m3 kg-1 s-1)
	201	REAL(wp) :: k_br = 1000.0_wp !< const. in breakup parametrization (m-1)
	202	REAL(wp) :: k_st = 1.2E8_wp !< const. in drizzle parametrization (m-1 s-1)
	203	REAL(wp) :: kappa_rr = 60.7_wp !< const. in collision kernel (kg-1/3)
	204	REAL(wp) :: kin_vis_air = 1.4086E-5_wp !< kin. viscosity of air (m2 s-1)
	205	REAL(wp) :: prec_time_const = 0.001_wp !< coef. in Kessler scheme (s-1)
	206	REAL(wp) :: ql_crit = 0.0005_wp !< coef. in Kessler scheme (kg kg-1)
	207	REAL(wp) :: schmidt_p_1d3=0.8921121_wp !< Schmidt number0.33333, 0.710.33333
	208	REAL(wp) :: sigma_gc = 1.3_wp !< geometric standard deviation cloud droplets
	209	REAL(wp) :: thermal_conductivity_l = 2.43E-2_wp !< therm. cond. air (J m-1 s-1 K-1)
	210	REAL(wp) :: w_precipitation = 9.65_wp !< maximum terminal velocity (m s-1)
	211	REAL(wp) :: x0 = 2.6E-10_wp !< separating drop mass (kg)
[2292]	212	REAL(wp) :: xamin = 5.24E-19_wp !< average aerosol mass (kg) (~ 0.05Âµm)
	213	REAL(wp) :: xcmin = 4.18E-15_wp !< minimum cloud drop size (kg) (~ 1Âµm)
	214	REAL(wp) :: xcmax = 2.6E-10_wp !< maximum cloud drop size (kg) (~ 40Âµm)
[1849]	215	REAL(wp) :: xrmin = 2.6E-10_wp !< minimum rain drop size (kg)
	216	REAL(wp) :: xrmax = 5.0E-6_wp !< maximum rain drop site (kg)
	217
[2375]	218	REAL(wp) :: c_sedimentation = 2.0_wp !< Courant number of sedimentation process
	219	REAL(wp) :: dpirho_l !< 6.0 / ( pi * rho_l )
	220	REAL(wp) :: dry_aerosol_radius = 0.05E-6_wp !< dry aerosol radius
	221	REAL(wp) :: dt_micro !< microphysics time step
	222	REAL(wp) :: sigma_bulk = 2.0_wp !< width of aerosol spectrum
	223	REAL(wp) :: na_init = 100.0E6_wp !< Total particle/aerosol concentration (cm-3)
	224	REAL(wp) :: nc_const = 70.0E6_wp !< cloud droplet concentration
	225	REAL(wp) :: dt_precipitation = 100.0_wp !< timestep precipitation (s)
	226	REAL(wp) :: sed_qc_const !< const. for sedimentation of cloud water
	227	REAL(wp) :: pirho_l !< pi * rho_l / 6.0;
[1849]	228
	229	REAL(wp), DIMENSION(:), ALLOCATABLE :: nc_1d !<
	230	REAL(wp), DIMENSION(:), ALLOCATABLE :: nr_1d !<
	231	REAL(wp), DIMENSION(:), ALLOCATABLE :: pt_1d !<
	232	REAL(wp), DIMENSION(:), ALLOCATABLE :: qc_1d !<
	233	REAL(wp), DIMENSION(:), ALLOCATABLE :: qr_1d !<
	234	REAL(wp), DIMENSION(:), ALLOCATABLE :: q_1d !<
	235
	236	SAVE
	237
[1000]	238	PRIVATE
[1849]	239	PUBLIC microphysics_control, microphysics_init
[1000]	240
[2292]	241	PUBLIC cloud_water_sedimentation, collision_turbulence, &
[2375]	242	curvature_solution_effects_bulk, c_sedimentation, &
	243	dry_aerosol_radius, dt_precipitation, &
	244	limiter_sedimentation, na_init, nc_const, sigma_gc, sigma_bulk, &
[1849]	245	ventilation_effect
	246
[2312]	247
[1115]	248	INTERFACE microphysics_control
	249	MODULE PROCEDURE microphysics_control
	250	MODULE PROCEDURE microphysics_control_ij
	251	END INTERFACE microphysics_control
[1022]	252
[1115]	253	INTERFACE adjust_cloud
	254	MODULE PROCEDURE adjust_cloud
	255	MODULE PROCEDURE adjust_cloud_ij
	256	END INTERFACE adjust_cloud
	257
[2292]	258	INTERFACE activation
	259	MODULE PROCEDURE activation
	260	MODULE PROCEDURE activation_ij
	261	END INTERFACE activation
	262
	263	INTERFACE condensation
	264	MODULE PROCEDURE condensation
	265	MODULE PROCEDURE condensation_ij
	266	END INTERFACE condensation
	267
[1000]	268	INTERFACE autoconversion
	269	MODULE PROCEDURE autoconversion
	270	MODULE PROCEDURE autoconversion_ij
	271	END INTERFACE autoconversion
	272
[1822]	273	INTERFACE autoconversion_kessler
	274	MODULE PROCEDURE autoconversion_kessler
	275	MODULE PROCEDURE autoconversion_kessler_ij
	276	END INTERFACE autoconversion_kessler
	277
[1000]	278	INTERFACE accretion
	279	MODULE PROCEDURE accretion
	280	MODULE PROCEDURE accretion_ij
	281	END INTERFACE accretion
[1005]	282
	283	INTERFACE selfcollection_breakup
	284	MODULE PROCEDURE selfcollection_breakup
	285	MODULE PROCEDURE selfcollection_breakup_ij
	286	END INTERFACE selfcollection_breakup
[1012]	287
	288	INTERFACE evaporation_rain
	289	MODULE PROCEDURE evaporation_rain
	290	MODULE PROCEDURE evaporation_rain_ij
	291	END INTERFACE evaporation_rain
	292
	293	INTERFACE sedimentation_cloud
	294	MODULE PROCEDURE sedimentation_cloud
	295	MODULE PROCEDURE sedimentation_cloud_ij
	296	END INTERFACE sedimentation_cloud
[2155]	297
[1012]	298	INTERFACE sedimentation_rain
	299	MODULE PROCEDURE sedimentation_rain
	300	MODULE PROCEDURE sedimentation_rain_ij
	301	END INTERFACE sedimentation_rain
	302
[1691]	303	INTERFACE calc_precipitation_amount
	304	MODULE PROCEDURE calc_precipitation_amount
	305	MODULE PROCEDURE calc_precipitation_amount_ij
	306	END INTERFACE calc_precipitation_amount
	307
[1000]	308	CONTAINS
[1849]	309	!------------------------------------------------------------------------------!
	310	! Description:
	311	! ------------
	312	!> Initialization of bulk microphysics
	313	!------------------------------------------------------------------------------!
	314	SUBROUTINE microphysics_init
[1000]	315
[1849]	316	USE arrays_3d, &
	317	ONLY: dzu
[1000]	318
[1849]	319	USE constants, &
	320	ONLY: pi
	321
	322	USE cloud_parameters, &
[2375]	323	ONLY: molecular_weight_of_solute, rho_l, rho_s, vanthoff
[1849]	324
	325	USE control_parameters, &
[2375]	326	ONLY: aerosol_nacl, aerosol_c3h4o4 , aerosol_nh4no3, &
	327	microphysics_morrison, microphysics_seifert
[1849]	328
	329	USE indices, &
	330	ONLY: nzb, nzt
	331
	332	IMPLICIT NONE
	333
	334	!
	335	!-- constant for the sedimentation of cloud water (2-moment cloud physics)
	336	sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l ) &
	337	)*( 2.0_wp / 3.0_wp ) &
	338	EXP( 5.0_wp * LOG( sigma_gc )**2 )
	339
	340	!
	341	!-- Calculate timestep according to precipitation
	342	IF ( microphysics_seifert ) THEN
	343	dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) / &
	344	w_precipitation
	345	ENDIF
	346
	347	!
[2375]	348	!-- Set constants for certain aerosol type
	349	IF ( microphysics_morrison ) THEN
	350	IF ( aerosol_nacl ) THEN
	351	molecular_weight_of_solute = 0.05844_wp
	352	rho_s = 2165.0_wp
	353	vanthoff = 2.0_wp
	354	ELSEIF ( aerosol_c3h4o4 ) THEN
	355	molecular_weight_of_solute = 0.10406_wp
	356	rho_s = 1600.0_wp
	357	vanthoff = 1.37_wp
	358	ELSEIF ( aerosol_nh4no3 ) THEN
	359	molecular_weight_of_solute = 0.08004_wp
	360	rho_s = 1720.0_wp
	361	vanthoff = 2.31_wp
	362	ENDIF
	363	ENDIF
	364
	365	!
[1849]	366	!-- Pre-calculate frequently calculated fractions of pi and rho_l
	367	pirho_l = pi * rho_l / 6.0_wp
	368	dpirho_l = 1.0_wp / pirho_l
	369
	370	!
	371	!-- Allocate 1D microphysics arrays
[2292]	372	ALLOCATE ( pt_1d(nzb:nzt+1), q_1d(nzb:nzt+1), &
[1849]	373	qc_1d(nzb:nzt+1) )
	374
[2292]	375	IF ( microphysics_morrison .OR. microphysics_seifert ) THEN
	376	ALLOCATE ( nc_1d(nzb:nzt+1) )
	377	ENDIF
	378
[1849]	379	IF ( microphysics_seifert ) THEN
	380	ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) )
	381	ENDIF
	382
	383	END SUBROUTINE microphysics_init
	384
	385
[1000]	386	!------------------------------------------------------------------------------!
[1682]	387	! Description:
	388	! ------------
[1849]	389	!> Control of microphysics for all grid points
[1000]	390	!------------------------------------------------------------------------------!
[1115]	391	SUBROUTINE microphysics_control
[1022]	392
[1361]	393	USE arrays_3d, &
[1849]	394	ONLY: hyp, pt_init, prr, zu
[1361]	395
	396	USE cloud_parameters, &
[1849]	397	ONLY: cp, hyrho, pt_d_t, r_d, t_d_pt
[1361]	398
	399	USE control_parameters, &
[1849]	400	ONLY: call_microphysics_at_all_substeps, dt_3d, g, &
	401	intermediate_timestep_count, large_scale_forcing, &
[2292]	402	lsf_surf, microphysics_kessler, microphysics_morrison, &
	403	microphysics_seifert, pt_surface, &
	404	rho_surface,surface_pressure
[1361]	405
	406	USE indices, &
	407	ONLY: nzb, nzt
	408
[1320]	409	USE kinds
[1115]	410
[1361]	411	USE statistics, &
	412	ONLY: weight_pres
	413
[1115]	414	IMPLICIT NONE
	415
[1682]	416	INTEGER(iwp) :: k !<
[1115]	417
[1682]	418	REAL(wp) :: t_surface !<
[1361]	419
	420	IF ( large_scale_forcing .AND. lsf_surf ) THEN
	421	!
	422	!-- Calculate:
	423	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
	424	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
	425	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
	426	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
	427	DO k = nzb, nzt+1
	428	hyp(k) = surface_pressure * 100.0_wp * &
	429	( ( t_surface - g / cp * zu(k) ) / &
	430	t_surface )**(1.0_wp / 0.286_wp)
	431	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
	432	t_d_pt(k) = 1.0_wp / pt_d_t(k)
[2155]	433	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
[1115]	434	ENDDO
[1822]	435
[1361]	436	!
	437	!-- Compute reference density
	438	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
	439	ENDIF
[1115]	440
[1361]	441	!
[2155]	442	!-- Compute length of time step
[1361]	443	IF ( call_microphysics_at_all_substeps ) THEN
	444	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
	445	ELSE
	446	dt_micro = dt_3d
	447	ENDIF
	448
	449	!
[1822]	450	!-- Reset precipitation rate
	451	IF ( intermediate_timestep_count == 1 ) prr = 0.0_wp
	452
	453	!
[1361]	454	!-- Compute cloud physics
[1822]	455	IF ( microphysics_kessler ) THEN
	456
	457	CALL autoconversion_kessler
[1831]	458	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud
[1822]	459
	460	ELSEIF ( microphysics_seifert ) THEN
	461
[1361]	462	CALL adjust_cloud
[2292]	463	IF ( microphysics_morrison ) CALL activation
	464	IF ( microphysics_morrison ) CALL condensation
[1361]	465	CALL autoconversion
	466	CALL accretion
	467	CALL selfcollection_breakup
	468	CALL evaporation_rain
	469	CALL sedimentation_rain
[1831]	470	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud
[1361]	471
[1691]	472	ENDIF
	473
[1822]	474	CALL calc_precipitation_amount
	475
[1115]	476	END SUBROUTINE microphysics_control
	477
[1682]	478	!------------------------------------------------------------------------------!
	479	! Description:
	480	! ------------
[2155]	481	!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
	482	!> evaporation of rain drops due to too small or too big weights
[1682]	483	!> of rain drops (Stevens and Seifert, 2008).
	484	!------------------------------------------------------------------------------!
[1115]	485	SUBROUTINE adjust_cloud
	486
[1361]	487	USE arrays_3d, &
[2292]	488	ONLY: qc, nc, qr, nr
[1361]	489
	490	USE cloud_parameters, &
[1849]	491	ONLY: hyrho
[1361]	492
[2292]	493	USE control_parameters, &
	494	ONLY: microphysics_morrison
	495
[1361]	496	USE cpulog, &
	497	ONLY: cpu_log, log_point_s
	498
	499	USE indices, &
[2317]	500	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]	501
[1320]	502	USE kinds
[1022]	503
	504	IMPLICIT NONE
	505
[1682]	506	INTEGER(iwp) :: i !<
	507	INTEGER(iwp) :: j !<
	508	INTEGER(iwp) :: k !<
[1022]	509
[1361]	510	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
	511
[2155]	512	DO i = nxlg, nxrg
	513	DO j = nysg, nyng
[2232]	514	DO k = nzb+1, nzt
[1361]	515	IF ( qr(k,j,i) <= eps_sb ) THEN
	516	qr(k,j,i) = 0.0_wp
	517	nr(k,j,i) = 0.0_wp
	518	ELSE
	519	IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN
[2232]	520	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * &
[2312]	521	MERGE( 1.0_wp, 0.0_wp, &
[2232]	522	BTEST( wall_flags_0(k,j,i), 0 ) )
[1361]	523	ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN
[2232]	524	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * &
[2312]	525	MERGE( 1.0_wp, 0.0_wp, &
[2232]	526	BTEST( wall_flags_0(k,j,i), 0 ) )
[1361]	527	ENDIF
	528	ENDIF
[2292]	529	IF ( microphysics_morrison ) THEN
	530	IF ( qc(k,j,i) <= eps_sb ) THEN
	531	qc(k,j,i) = 0.0_wp
	532	nc(k,j,i) = 0.0_wp
	533	ELSE
	534	IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN
	535	nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin * &
[2312]	536	MERGE( 1.0_wp, 0.0_wp, &
[2292]	537	BTEST( wall_flags_0(k,j,i), 0 ) )
	538	ENDIF
	539	ENDIF
	540	ENDIF
[1022]	541	ENDDO
	542	ENDDO
	543	ENDDO
	544
[1361]	545	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
	546
[1115]	547	END SUBROUTINE adjust_cloud
[1022]	548
[2292]	549	!------------------------------------------------------------------------------!
	550	! Description:
	551	! ------------
	552	!> Calculate number of activated condensation nucleii after simple activation
	553	!> scheme of Twomey, 1959.
	554	!------------------------------------------------------------------------------!
	555	SUBROUTINE activation
[1106]	556
[2292]	557	USE arrays_3d, &
	558	ONLY: hyp, nc, nr, pt, q, qc, qr
	559
	560	USE cloud_parameters, &
[2375]	561	ONLY: hyrho, l_d_cp, l_d_r, l_v, molecular_weight_of_solute, &
	562	molecular_weight_of_water, rho_l, rho_s, r_v, t_d_pt, &
	563	vanthoff
[2292]	564
	565	USE constants, &
	566	ONLY: pi
	567
	568	USE cpulog, &
	569	ONLY: cpu_log, log_point_s
	570
[2608]	571	USE diagnostic_quantities_mod, &
	572	ONLY: e_s, magnus, q_s, sat, supersaturation, t_l
	573
[2292]	574	USE indices, &
	575	ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt
	576
	577	USE kinds
	578
	579	USE control_parameters, &
	580	ONLY: simulated_time
	581
	582	IMPLICIT NONE
	583
	584	INTEGER(iwp) :: i !<
	585	INTEGER(iwp) :: j !<
	586	INTEGER(iwp) :: k !<
	587
	588	REAL(wp) :: activ !<
	589	REAL(wp) :: afactor !<
	590	REAL(wp) :: beta_act !<
	591	REAL(wp) :: bfactor !<
	592	REAL(wp) :: k_act !<
	593	REAL(wp) :: n_act !<
	594	REAL(wp) :: n_ccn !<
	595	REAL(wp) :: s_0 !<
	596	REAL(wp) :: sat_max !<
	597	REAL(wp) :: sigma !<
[2375]	598	REAL(wp) :: sigma_act !<
[2292]	599
	600	CALL cpu_log( log_point_s(65), 'activation', 'start' )
	601
	602	DO i = nxlg, nxrg
	603	DO j = nysg, nyng
	604	DO k = nzb+1, nzt
	605
	606	!
[2608]	607	!-- Call calculation of supersaturation located
	608	!-- in diagnostic_quantities_mod
	609	CALL supersaturation ( i, j, k )
[2292]	610	!
	611	!-- Prescribe parameters for activation
	612	!-- (see: Bott + Trautmann, 2002, Atm. Res., 64)
[2312]	613	k_act = 0.7_wp
[2522]	614	activ = 0.0_wp
[2292]	615
[2522]	616
	617	IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
[2292]	618	!
[2312]	619	!-- Compute the number of activated Aerosols
[2292]	620	!-- (see: Twomey, 1959, Pure and applied Geophysics, 43)
	621	n_act = na_init * sat**k_act
	622	!
	623	!-- Compute the number of cloud droplets
	624	!-- (see: Morrison + Grabowski, 2007, JAS, 64)
	625	! activ = MAX( n_act - nc(k,j,i), 0.0_wp) / dt_micro
	626
	627	!
	628	!-- Compute activation rate after Khairoutdinov and Kogan
[2312]	629	!-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
	630	sat_max = 1.0_wp / 100.0_wp
	631	activ = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN &
	632	( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) / &
[2292]	633	dt_micro
[2522]	634	ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk ) THEN
[2292]	635	!
[2312]	636	!-- Curvature effect (afactor) with surface tension
[2292]	637	!-- parameterization by Straka (2009)
[2608]	638	sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
	639	afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
[2292]	640	!
	641	!-- Solute effect (bfactor)
	642	bfactor = vanthoff * molecular_weight_of_water * &
	643	rho_s / ( molecular_weight_of_solute * rho_l )
	644
	645	!
[2312]	646	!-- Prescribe power index that describes the soluble fraction
[2375]	647	!-- of an aerosol particle (beta)
[2292]	648	!-- (see: Morrison + Grabowski, 2007, JAS, 64)
[2375]	649	beta_act = 0.5_wp
	650	sigma_act = sigma_bulk**( 1.0_wp + beta_act )
[2292]	651	!
[2312]	652	!-- Calculate mean geometric supersaturation (s_0) with
[2292]	653	!-- parameterization by Khvorostyanov and Curry (2006)
[2375]	654	s_0 = dry_aerosol_radius *(- ( 1.0_wp + beta_act ) ) &
[2292]	655	( 4.0_wp * afactor*3 / ( 27.0_wp bfactor ) )**0.5_wp
	656
	657	!
[2312]	658	!-- Calculate number of activated CCN as a function of
[2292]	659	!-- supersaturation and taking Koehler theory into account
[2312]	660	!-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
	661	n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( &
[2375]	662	LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
[2292]	663	activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp )
	664	ENDIF
	665
	666	nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)
	667
	668	ENDDO
	669	ENDDO
	670	ENDDO
	671
	672	CALL cpu_log( log_point_s(65), 'activation', 'stop' )
	673
	674	END SUBROUTINE activation
	675
	676
[1682]	677	!------------------------------------------------------------------------------!
	678	! Description:
	679	! ------------
[2312]	680	!> Calculate condensation rate for cloud water content (after Khairoutdinov and
[2292]	681	!> Kogan, 2000).
	682	!------------------------------------------------------------------------------!
	683	SUBROUTINE condensation
	684
	685	USE arrays_3d, &
	686	ONLY: hyp, nr, pt, q, qc, qr, nc
	687
	688	USE cloud_parameters, &
	689	ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
	690
	691	USE constants, &
	692	ONLY: pi
	693
	694	USE cpulog, &
	695	ONLY: cpu_log, log_point_s
	696
[2608]	697	USE diagnostic_quantities_mod, &
	698	ONLY: e_s, magnus, q_s, sat, supersaturation, t_l
	699
[2292]	700	USE indices, &
	701	ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt
	702
	703	USE kinds
	704
	705	USE control_parameters, &
	706	ONLY: simulated_time
	707
	708	IMPLICIT NONE
	709
	710	INTEGER(iwp) :: i !<
	711	INTEGER(iwp) :: j !<
	712	INTEGER(iwp) :: k !<
	713
	714	REAL(wp) :: cond !<
	715	REAL(wp) :: cond_max !<
	716	REAL(wp) :: dc !<
	717	REAL(wp) :: evap !<
	718	REAL(wp) :: evap_nc !<
	719	REAL(wp) :: g_fac !<
	720	REAL(wp) :: nc_0 !<
	721	REAL(wp) :: temp !<
	722	REAL(wp) :: xc !<
	723
	724	CALL cpu_log( log_point_s(66), 'condensation', 'start' )
	725
	726	DO i = nxlg, nxrg
	727	DO j = nysg, nyng
	728	DO k = nzb+1, nzt
	729	!
[2608]	730	!-- Call calculation of supersaturation located
	731	!-- in diagnostic_quantities_mod
	732	CALL supersaturation ( i, j, k )
[2292]	733	!
	734	!-- Actual temperature:
	735	temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
	736
	737	g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
	738	l_v / ( thermal_conductivity_l * temp ) &
	739	+ r_v * temp / ( diff_coeff_l * e_s ) &
	740	)
	741	!
	742	!-- Mean weight of cloud drops
	743	IF ( nc(k,j,i) <= 0.0_wp) CYCLE
[2312]	744	xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
[2292]	745	!
	746	!-- Weight averaged diameter of cloud drops:
	747	dc = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
	748	!
	749	!-- Integral diameter of cloud drops
[2312]	750	nc_0 = nc(k,j,i) * dc
[2292]	751	!
	752	!-- Condensation needs only to be calculated in supersaturated regions
	753	IF ( sat > 0.0_wp ) THEN
	754	!
	755	!-- Condensation rate of cloud water content
	756	!-- after KK scheme.
	757	!-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
	758	cond = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
	759	cond_max = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
	760	cond = MIN( cond, cond_max / dt_micro )
[2312]	761
[2292]	762	qc(k,j,i) = qc(k,j,i) + cond * dt_micro
	763	ELSEIF ( sat < 0.0_wp ) THEN
	764	evap = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
	765	evap = MAX( evap, -qc(k,j,i) / dt_micro )
	766
	767	qc(k,j,i) = qc(k,j,i) + evap * dt_micro
	768	ENDIF
	769	IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN
	770	nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin
	771	ENDIF
	772	ENDDO
	773	ENDDO
	774	ENDDO
	775
	776	CALL cpu_log( log_point_s(66), 'condensation', 'stop' )
	777
	778	END SUBROUTINE condensation
	779
	780
	781	!------------------------------------------------------------------------------!
	782	! Description:
	783	! ------------
[1682]	784	!> Autoconversion rate (Seifert and Beheng, 2006).
	785	!------------------------------------------------------------------------------!
[1000]	786	SUBROUTINE autoconversion
	787
[1361]	788	USE arrays_3d, &
[2292]	789	ONLY: diss, dzu, nc, nr, qc, qr
[1361]	790
	791	USE cloud_parameters, &
[1849]	792	ONLY: hyrho
[1361]	793
	794	USE control_parameters, &
[2292]	795	ONLY: microphysics_morrison, rho_surface
[1361]	796
	797	USE cpulog, &
	798	ONLY: cpu_log, log_point_s
	799
	800	USE grid_variables, &
	801	ONLY: dx, dy
	802
	803	USE indices, &
[2317]	804	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]	805
[1320]	806	USE kinds
[1000]	807
	808	IMPLICIT NONE
	809
[1682]	810	INTEGER(iwp) :: i !<
	811	INTEGER(iwp) :: j !<
	812	INTEGER(iwp) :: k !<
[1000]	813
[2155]	814	REAL(wp) :: alpha_cc !<
[1682]	815	REAL(wp) :: autocon !<
	816	REAL(wp) :: dissipation !<
[2232]	817	REAL(wp) :: flag !< flag to mask topography grid points
[1682]	818	REAL(wp) :: k_au !<
	819	REAL(wp) :: l_mix !<
[2292]	820	REAL(wp) :: nc_auto !<
[1682]	821	REAL(wp) :: nu_c !<
	822	REAL(wp) :: phi_au !<
	823	REAL(wp) :: r_cc !<
	824	REAL(wp) :: rc !<
	825	REAL(wp) :: re_lambda !<
	826	REAL(wp) :: sigma_cc !<
	827	REAL(wp) :: tau_cloud !<
	828	REAL(wp) :: xc !<
[1361]	829
	830	CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
	831
[2155]	832	DO i = nxlg, nxrg
	833	DO j = nysg, nyng
[2232]	834	DO k = nzb+1, nzt
	835	!
	836	!-- Predetermine flag to mask topography
	837	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1000]	838
[2292]	839	nc_auto = MERGE( nc(k,j,i), nc_const, microphysics_morrison )
	840
[2522]	841	IF ( qc(k,j,i) > eps_sb .AND. nc_auto > eps_mr ) THEN
[1361]	842
	843	k_au = k_cc / ( 20.0_wp * x0 )
	844	!
	845	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
	846	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
[2522]	847	tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + &
	848	qc(k,j,i) ), 0.0_wp )
[1361]	849	!
[2155]	850	!-- Universal function for autoconversion process
[1361]	851	!-- (Seifert and Beheng, 2006):
	852	phi_au = 600.0_wp * tau_cloud*0.68_wp &
	853	( 1.0_wp - tau_cloud0.68_wp )3
	854	!
	855	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
	856	!-- (Use constant nu_c = 1.0_wp instead?)
	857	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
	858	!
	859	!-- Mean weight of cloud droplets:
[2522]	860	xc = MAX( hyrho(k) * qc(k,j,i) / nc_auto, xcmin)
[1361]	861	!
[2155]	862	!-- Parameterized turbulence effects on autoconversion (Seifert,
[1361]	863	!-- Nuijens and Stevens, 2010)
[1831]	864	IF ( collision_turbulence ) THEN
[1361]	865	!
	866	!-- Weight averaged radius of cloud droplets:
	867	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
	868
	869	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
	870	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
	871	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
	872	!
[2155]	873	!-- Mixing length (neglecting distance to ground and
[1361]	874	!-- stratification)
	875	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
	876	!
[2155]	877	!-- Limit dissipation rate according to Seifert, Nuijens and
[1361]	878	!-- Stevens (2010)
	879	dissipation = MIN( 0.06_wp, diss(k,j,i) )
	880	!
	881	!-- Compute Taylor-microscale Reynolds number:
	882	re_lambda = 6.0_wp / 11.0_wp * &
	883	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
	884	SQRT( 15.0_wp / kin_vis_air ) * &
	885	dissipation**( 1.0_wp / 6.0_wp )
	886	!
[2155]	887	!-- The factor of 1.0E4 is needed to convert the dissipation
[1361]	888	!-- rate from m2 s-3 to cm2 s-3.
	889	k_au = k_au * ( 1.0_wp + &
	890	dissipation * 1.0E4_wp * &
	891	( re_lambda * 1.0E-3_wp )*0.25_wp &
[2155]	892	( alpha_cc * EXP( -1.0_wp * ( ( rc - &
[1361]	893	r_cc ) / &
	894	sigma_cc )**2 &
	895	) + beta_cc &
	896	) &
	897	)
	898	ENDIF
	899	!
	900	!-- Autoconversion rate (Seifert and Beheng, 2006):
	901	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
	902	( nu_c + 1.0_wp )*2 qc(k,j,i)*2 xc*2 &
	903	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
	904	rho_surface
	905	autocon = MIN( autocon, qc(k,j,i) / dt_micro )
	906
[2232]	907	qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag
	908	qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag
	909	nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro &
	910	* flag
[2292]	911	IF ( microphysics_morrison ) THEN
[2312]	912	nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp * &
[2292]	913	autocon / x0 * hyrho(k) * dt_micro * flag )
	914	ENDIF
[1361]	915
	916	ENDIF
	917
[1000]	918	ENDDO
	919	ENDDO
	920	ENDDO
	921
[1361]	922	CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
	923
[1000]	924	END SUBROUTINE autoconversion
	925
[1106]	926
[1682]	927	!------------------------------------------------------------------------------!
	928	! Description:
	929	! ------------
[1822]	930	!> Autoconversion process (Kessler, 1969).
	931	!------------------------------------------------------------------------------!
	932	SUBROUTINE autoconversion_kessler
	933
	934	USE arrays_3d, &
[1849]	935	ONLY: dzw, pt, prr, q, qc
[1822]	936
	937	USE cloud_parameters, &
[1849]	938	ONLY: l_d_cp, pt_d_t
[1822]	939
	940	USE indices, &
[2317]	941	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1822]	942
	943	USE kinds
	944
[2317]	945	USE surface_mod, &
[2698]	946	ONLY: get_topography_top_index_ji
[1822]	947
[2317]	948
[1822]	949	IMPLICIT NONE
	950
[2232]	951	INTEGER(iwp) :: i !<
	952	INTEGER(iwp) :: j !<
	953	INTEGER(iwp) :: k !<
	954	INTEGER(iwp) :: k_wall !< topgraphy top index
[1822]	955
	956	REAL(wp) :: dqdt_precip !<
[2232]	957	REAL(wp) :: flag !< flag to mask topography grid points
[1822]	958
[2155]	959	DO i = nxlg, nxrg
	960	DO j = nysg, nyng
[2232]	961	!
	962	!-- Determine vertical index of topography top
[2698]	963	k_wall = get_topography_top_index_ji( j, i, 's' )
[2232]	964	DO k = nzb+1, nzt
	965	!
	966	!-- Predetermine flag to mask topography
	967	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1822]	968
	969	IF ( qc(k,j,i) > ql_crit ) THEN
	970	dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
	971	ELSE
	972	dqdt_precip = 0.0_wp
	973	ENDIF
	974
[2232]	975	qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
	976	q(k,j,i) = q(k,j,i) - dqdt_precip * dt_micro * flag
[1845]	977	pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * l_d_cp * &
[2232]	978	pt_d_t(k) * flag
[1822]	979
	980	!
[1845]	981	!-- Compute the rain rate (stored on surface grid point)
[2232]	982	prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
[1822]	983
	984	ENDDO
	985	ENDDO
	986	ENDDO
	987
	988	END SUBROUTINE autoconversion_kessler
	989
	990
	991	!------------------------------------------------------------------------------!
	992	! Description:
	993	! ------------
[1682]	994	!> Accretion rate (Seifert and Beheng, 2006).
	995	!------------------------------------------------------------------------------!
[1005]	996	SUBROUTINE accretion
[1000]	997
[1361]	998	USE arrays_3d, &
[2292]	999	ONLY: diss, qc, qr, nc
[1361]	1000
	1001	USE cloud_parameters, &
[1849]	1002	ONLY: hyrho
[1361]	1003
	1004	USE control_parameters, &
[2292]	1005	ONLY: microphysics_morrison, rho_surface
[1361]	1006
	1007	USE cpulog, &
	1008	ONLY: cpu_log, log_point_s
	1009
	1010	USE indices, &
[2317]	1011	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]	1012
[1320]	1013	USE kinds
[1005]	1014
[1000]	1015	IMPLICIT NONE
	1016
[1682]	1017	INTEGER(iwp) :: i !<
	1018	INTEGER(iwp) :: j !<
	1019	INTEGER(iwp) :: k !<
[1000]	1020
[1682]	1021	REAL(wp) :: accr !<
[2232]	1022	REAL(wp) :: flag !< flag to mask topography grid points
[1682]	1023	REAL(wp) :: k_cr !<
[2292]	1024	REAL(wp) :: nc_accr !<
[1682]	1025	REAL(wp) :: phi_ac !<
	1026	REAL(wp) :: tau_cloud !<
[2292]	1027	REAL(wp) :: xc !<
[1361]	1028
[2292]	1029
[1361]	1030	CALL cpu_log( log_point_s(56), 'accretion', 'start' )
	1031
[2155]	1032	DO i = nxlg, nxrg
	1033	DO j = nysg, nyng
[2232]	1034	DO k = nzb+1, nzt
	1035	!
	1036	!-- Predetermine flag to mask topography
	1037	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1000]	1038
[2292]	1039	nc_accr = MERGE( nc(k,j,i), nc_const, microphysics_morrison )
	1040
	1041	IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) &
	1042	.AND. ( nc_accr > eps_mr ) ) THEN
[1361]	1043	!
	1044	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[2155]	1045	tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
[1361]	1046	!
[2155]	1047	!-- Universal function for accretion process (Seifert and
[1361]	1048	!-- Beheng, 2001):
	1049	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
[2292]	1050
[1361]	1051	!
[2292]	1052	!-- Mean weight of cloud drops
[2312]	1053	xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
[2292]	1054	!
[2155]	1055	!-- Parameterized turbulence effects on autoconversion (Seifert,
	1056	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
[1361]	1057	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
[1831]	1058	IF ( collision_turbulence ) THEN
[1361]	1059	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
	1060	MIN( 600.0_wp, &
	1061	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
	1062	)
	1063	ELSE
[2155]	1064	k_cr = k_cr0
[1361]	1065	ENDIF
	1066	!
	1067	!-- Accretion rate (Seifert and Beheng, 2006):
	1068	accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * &
	1069	SQRT( rho_surface * hyrho(k) )
	1070	accr = MIN( accr, qc(k,j,i) / dt_micro )
	1071
[2232]	1072	qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
	1073	qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
[2292]	1074	IF ( microphysics_morrison ) THEN
	1075	nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), &
	1076	accr / xc * hyrho(k) * dt_micro * flag)
	1077	ENDIF
[1361]	1078
	1079	ENDIF
	1080
[1005]	1081	ENDDO
	1082	ENDDO
[1000]	1083	ENDDO
	1084
[1361]	1085	CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
	1086
[1005]	1087	END SUBROUTINE accretion
[1000]	1088
[1106]	1089
[1682]	1090	!------------------------------------------------------------------------------!
	1091	! Description:
	1092	! ------------
	1093	!> Collisional breakup rate (Seifert, 2008).
	1094	!------------------------------------------------------------------------------!
[1005]	1095	SUBROUTINE selfcollection_breakup
[1000]	1096
[1361]	1097	USE arrays_3d, &
	1098	ONLY: nr, qr
	1099
	1100	USE cloud_parameters, &
[1849]	1101	ONLY: hyrho
[1361]	1102
	1103	USE control_parameters, &
[1849]	1104	ONLY: rho_surface
[1361]	1105
	1106	USE cpulog, &
	1107	ONLY: cpu_log, log_point_s
	1108
	1109	USE indices, &
[2317]	1110	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]	1111
[1320]	1112	USE kinds
[2155]	1113
[1000]	1114	IMPLICIT NONE
	1115
[1682]	1116	INTEGER(iwp) :: i !<
	1117	INTEGER(iwp) :: j !<
	1118	INTEGER(iwp) :: k !<
[1000]	1119
[1682]	1120	REAL(wp) :: breakup !<
	1121	REAL(wp) :: dr !<
[2232]	1122	REAL(wp) :: flag !< flag to mask topography grid points
[1682]	1123	REAL(wp) :: phi_br !<
	1124	REAL(wp) :: selfcoll !<
[1361]	1125
	1126	CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
	1127
[2155]	1128	DO i = nxlg, nxrg
	1129	DO j = nysg, nyng
[2232]	1130	DO k = nzb+1, nzt
	1131	!
	1132	!-- Predetermine flag to mask topography
	1133	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1134
[1361]	1135	IF ( qr(k,j,i) > eps_sb ) THEN
	1136	!
	1137	!-- Selfcollection rate (Seifert and Beheng, 2001):
	1138	selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * &
	1139	SQRT( hyrho(k) * rho_surface )
	1140	!
	1141	!-- Weight averaged diameter of rain drops:
	1142	dr = ( hyrho(k) * qr(k,j,i) / &
	1143	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
	1144	!
	1145	!-- Collisional breakup rate (Seifert, 2008):
	1146	IF ( dr >= 0.3E-3_wp ) THEN
	1147	phi_br = k_br * ( dr - 1.1E-3_wp )
	1148	breakup = selfcoll * ( phi_br + 1.0_wp )
	1149	ELSE
	1150	breakup = 0.0_wp
	1151	ENDIF
[1000]	1152
[1361]	1153	selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
[2232]	1154	nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
[1361]	1155
[2155]	1156	ENDIF
[1000]	1157	ENDDO
	1158	ENDDO
	1159	ENDDO
	1160
[1361]	1161	CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
	1162
[1005]	1163	END SUBROUTINE selfcollection_breakup
[1000]	1164
[1106]	1165
[1682]	1166	!------------------------------------------------------------------------------!
	1167	! Description:
	1168	! ------------
[2155]	1169	!> Evaporation of precipitable water. Condensation is neglected for
[1682]	1170	!> precipitable water.
	1171	!------------------------------------------------------------------------------!
[1012]	1172	SUBROUTINE evaporation_rain
[1000]	1173
[1361]	1174	USE arrays_3d, &
	1175	ONLY: hyp, nr, pt, q, qc, qr
	1176
	1177	USE cloud_parameters, &
[1849]	1178	ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
[1361]	1179
	1180	USE constants, &
	1181	ONLY: pi
	1182
	1183	USE cpulog, &
	1184	ONLY: cpu_log, log_point_s
	1185
[2608]	1186	USE diagnostic_quantities_mod, &
	1187	ONLY: e_s, magnus, q_s, sat, supersaturation, t_l
	1188
[1361]	1189	USE indices, &
[2317]	1190	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]	1191
[1320]	1192	USE kinds
[1012]	1193
	1194	IMPLICIT NONE
	1195
[1682]	1196	INTEGER(iwp) :: i !<
	1197	INTEGER(iwp) :: j !<
	1198	INTEGER(iwp) :: k !<
[1361]	1199
[1682]	1200	REAL(wp) :: dr !<
	1201	REAL(wp) :: evap !<
	1202	REAL(wp) :: evap_nr !<
	1203	REAL(wp) :: f_vent !<
[2232]	1204	REAL(wp) :: flag !< flag to mask topography grid points
[1682]	1205	REAL(wp) :: g_evap !<
	1206	REAL(wp) :: lambda_r !<
	1207	REAL(wp) :: mu_r !<
	1208	REAL(wp) :: mu_r_2 !<
	1209	REAL(wp) :: mu_r_5d2 !<
	1210	REAL(wp) :: nr_0 !<
	1211	REAL(wp) :: temp !<
	1212	REAL(wp) :: xr !<
[1361]	1213
	1214	CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
	1215
[2155]	1216	DO i = nxlg, nxrg
	1217	DO j = nysg, nyng
[2232]	1218	DO k = nzb+1, nzt
	1219	!
	1220	!-- Predetermine flag to mask topography
	1221	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1222
[1361]	1223	IF ( qr(k,j,i) > eps_sb ) THEN
[2608]	1224
[1361]	1225	!
[2608]	1226	!-- Call calculation of supersaturation located
	1227	!-- in diagnostic_quantities_mod
	1228	CALL supersaturation ( i, j, k )
[1361]	1229	!
	1230	!-- Evaporation needs only to be calculated in subsaturated regions
	1231	IF ( sat < 0.0_wp ) THEN
	1232	!
	1233	!-- Actual temperature:
	1234	temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
[2155]	1235
[1361]	1236	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
	1237	l_v / ( thermal_conductivity_l * temp ) &
	1238	+ r_v * temp / ( diff_coeff_l * e_s ) &
	1239	)
	1240	!
	1241	!-- Mean weight of rain drops
	1242	xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
	1243	!
	1244	!-- Weight averaged diameter of rain drops:
	1245	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
	1246	!
[2155]	1247	!-- Compute ventilation factor and intercept parameter
[1361]	1248	!-- (Seifert and Beheng, 2006; Seifert, 2008):
	1249	IF ( ventilation_effect ) THEN
	1250	!
[2155]	1251	!-- Shape parameter of gamma distribution (Milbrandt and Yau,
[1361]	1252	!-- 2005; Stevens and Seifert, 2008):
	1253	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
	1254	( dr - 1.4E-3_wp ) ) )
	1255	!
	1256	!-- Slope parameter of gamma distribution (Seifert, 2008):
	1257	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	1258	( mu_r + 1.0_wp ) &
	1259	)**( 1.0_wp / 3.0_wp ) / dr
[1012]	1260
[1361]	1261	mu_r_2 = mu_r + 2.0_wp
[2155]	1262	mu_r_5d2 = mu_r + 2.5_wp
[1361]	1263
	1264	f_vent = a_vent * gamm( mu_r_2 ) * &
	1265	lambda_r*( -mu_r_2 ) + b_vent &
	1266	schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
	1267	gamm( mu_r_5d2 ) * lambda_r*( -mu_r_5d2 ) &
	1268	( 1.0_wp - &
	1269	0.5_wp * ( b_term / a_term ) * &
	1270	( lambda_r / ( c_term + lambda_r ) &
	1271	)**mu_r_5d2 - &
	1272	0.125_wp * ( b_term / a_term )*2 &
	1273	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
	1274	)**mu_r_5d2 - &
	1275	0.0625_wp * ( b_term / a_term )*3 &
	1276	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
	1277	)**mu_r_5d2 - &
[2155]	1278	0.0390625_wp * ( b_term / a_term )*4 &
[1361]	1279	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
	1280	)**mu_r_5d2 &
	1281	)
	1282
	1283	nr_0 = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) / &
[2155]	1284	gamm( mu_r + 1.0_wp )
[1361]	1285	ELSE
	1286	f_vent = 1.0_wp
	1287	nr_0 = nr(k,j,i) * dr
	1288	ENDIF
	1289	!
[2155]	1290	!-- Evaporation rate of rain water content (Seifert and
[1361]	1291	!-- Beheng, 2006):
	1292	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / &
	1293	hyrho(k)
	1294	evap = MAX( evap, -qr(k,j,i) / dt_micro )
	1295	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
	1296	-nr(k,j,i) / dt_micro )
	1297
[2232]	1298	qr(k,j,i) = qr(k,j,i) + evap * dt_micro * flag
	1299	nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag
[1361]	1300
	1301	ENDIF
[2155]	1302	ENDIF
[1361]	1303
[1012]	1304	ENDDO
	1305	ENDDO
	1306	ENDDO
	1307
[1361]	1308	CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
	1309
[1012]	1310	END SUBROUTINE evaporation_rain
	1311
[1106]	1312
[1682]	1313	!------------------------------------------------------------------------------!
	1314	! Description:
	1315	! ------------
	1316	!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
	1317	!------------------------------------------------------------------------------!
[1012]	1318	SUBROUTINE sedimentation_cloud
	1319
[1361]	1320	USE arrays_3d, &
[2292]	1321	ONLY: ddzu, dzu, nc, pt, prr, q, qc
[1361]	1322
	1323	USE cloud_parameters, &
[1849]	1324	ONLY: hyrho, l_d_cp, pt_d_t
[1361]	1325
	1326	USE control_parameters, &
[2312]	1327	ONLY: call_microphysics_at_all_substeps, &
[2292]	1328	intermediate_timestep_count, microphysics_morrison
[1361]	1329
	1330	USE cpulog, &
	1331	ONLY: cpu_log, log_point_s
	1332
	1333	USE indices, &
[2317]	1334	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]	1335
[1320]	1336	USE kinds
[1691]	1337
	1338	USE statistics, &
	1339	ONLY: weight_substep
	1340
[2155]	1341
[1012]	1342	IMPLICIT NONE
	1343
[1849]	1344	INTEGER(iwp) :: i !<
	1345	INTEGER(iwp) :: j !<
	1346	INTEGER(iwp) :: k !<
[1361]	1347
[2292]	1348	REAL(wp) :: flag !< flag to mask topography grid points
	1349	REAL(wp) :: nc_sedi !<
	1350
[2232]	1351	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
[2292]	1352	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nc !<
[1361]	1353
[2292]	1354
[1361]	1355	CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
	1356
	1357	sed_qc(nzt+1) = 0.0_wp
[2375]	1358	sed_nc(nzt+1) = 0.0_wp
[1361]	1359
[2155]	1360	DO i = nxlg, nxrg
	1361	DO j = nysg, nyng
[2232]	1362	DO k = nzt, nzb+1, -1
	1363	!
	1364	!-- Predetermine flag to mask topography
	1365	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[2292]	1366	nc_sedi = MERGE ( nc(k,j,i), nc_const, microphysics_morrison )
[1012]	1367
[2292]	1368	!
	1369	!-- Sedimentation fluxes for number concentration are only calculated
	1370	!-- for cloud_scheme = 'morrison'
	1371	IF ( microphysics_morrison ) THEN
	1372	IF ( qc(k,j,i) > eps_sb .AND. nc(k,j,i) > eps_mr ) THEN
[2312]	1373	sed_nc(k) = sed_qc_const * &
[2292]	1374	( qc(k,j,i) * hyrho(k) )*( 2.0_wp / 3.0_wp ) &
	1375	( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
	1376	ELSE
	1377	sed_nc(k) = 0.0_wp
	1378	ENDIF
	1379
	1380	sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) * &
	1381	nc(k,j,i) / dt_micro + sed_nc(k+1) &
	1382	) * flag
	1383
	1384	nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) * &
	1385	ddzu(k+1) / hyrho(k) * dt_micro * flag
	1386	ENDIF
	1387
[2375]	1388	IF ( qc(k,j,i) > eps_sb .AND. nc_sedi > eps_mr ) THEN
[2292]	1389	sed_qc(k) = sed_qc_const * nc_sedi*( -2.0_wp / 3.0_wp ) &
[2232]	1390	( qc(k,j,i) * hyrho(k) )*( 5.0_wp / 3.0_wp ) &
	1391	flag
[1361]	1392	ELSE
	1393	sed_qc(k) = 0.0_wp
	1394	ENDIF
	1395
	1396	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / &
	1397	dt_micro + sed_qc(k+1) &
[2232]	1398	) * flag
[1361]	1399
	1400	q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
[2232]	1401	ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]	1402	qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
[2232]	1403	ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]	1404	pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * &
	1405	ddzu(k+1) / hyrho(k) * l_d_cp * &
[2232]	1406	pt_d_t(k) * dt_micro * flag
[1361]	1407
[1691]	1408	!
	1409	!-- Compute the precipitation rate due to cloud (fog) droplets
[1822]	1410	IF ( call_microphysics_at_all_substeps ) THEN
	1411	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) &
[2232]	1412	* weight_substep(intermediate_timestep_count) &
	1413	* flag
[1822]	1414	ELSE
[2232]	1415	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
[1691]	1416	ENDIF
	1417
[1012]	1418	ENDDO
	1419	ENDDO
	1420	ENDDO
	1421
[1361]	1422	CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
	1423
[1012]	1424	END SUBROUTINE sedimentation_cloud
	1425
[1106]	1426
[1682]	1427	!------------------------------------------------------------------------------!
	1428	! Description:
	1429	! ------------
	1430	!> Computation of sedimentation flux. Implementation according to Stevens
	1431	!> and Seifert (2008). Code is based on UCLA-LES.
	1432	!------------------------------------------------------------------------------!
[1012]	1433	SUBROUTINE sedimentation_rain
	1434
[1361]	1435	USE arrays_3d, &
[1849]	1436	ONLY: ddzu, dzu, nr, pt, prr, q, qr
[1361]	1437
	1438	USE cloud_parameters, &
[1849]	1439	ONLY: hyrho, l_d_cp, pt_d_t
[1361]	1440
	1441	USE control_parameters, &
[1849]	1442	ONLY: call_microphysics_at_all_substeps, intermediate_timestep_count
[1361]	1443	USE cpulog, &
	1444	ONLY: cpu_log, log_point_s
	1445
	1446	USE indices, &
[2317]	1447	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]	1448
[1320]	1449	USE kinds
[1012]	1450
[1361]	1451	USE statistics, &
	1452	ONLY: weight_substep
[2155]	1453
[2232]	1454	USE surface_mod, &
	1455	ONLY : bc_h
[2312]	1456
[1012]	1457	IMPLICIT NONE
	1458
[2232]	1459	INTEGER(iwp) :: i !< running index x direction
	1460	INTEGER(iwp) :: j !< running index y direction
	1461	INTEGER(iwp) :: k !< running index z direction
[2312]	1462	INTEGER(iwp) :: k_run !<
[2232]	1463	INTEGER(iwp) :: l !< running index of surface type
	1464	INTEGER(iwp) :: m !< running index surface elements
	1465	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	1466	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1361]	1467
[1682]	1468	REAL(wp) :: c_run !<
	1469	REAL(wp) :: d_max !<
	1470	REAL(wp) :: d_mean !<
	1471	REAL(wp) :: d_min !<
	1472	REAL(wp) :: dr !<
	1473	REAL(wp) :: flux !<
[2232]	1474	REAL(wp) :: flag !< flag to mask topography grid points
[1682]	1475	REAL(wp) :: lambda_r !<
	1476	REAL(wp) :: mu_r !<
	1477	REAL(wp) :: z_run !<
[1361]	1478
[1682]	1479	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<
	1480	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<
	1481	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<
	1482	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<
	1483	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<
	1484	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<
	1485	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<
	1486	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<
[1361]	1487
	1488	CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
[1682]	1489
[1361]	1490	!
[2155]	1491	!-- Compute velocities
	1492	DO i = nxlg, nxrg
	1493	DO j = nysg, nyng
[2232]	1494	DO k = nzb+1, nzt
	1495	!
	1496	!-- Predetermine flag to mask topography
	1497	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1498
[1361]	1499	IF ( qr(k,j,i) > eps_sb ) THEN
	1500	!
	1501	!-- Weight averaged diameter of rain drops:
	1502	dr = ( hyrho(k) * qr(k,j,i) / &
	1503	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
	1504	!
	1505	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	1506	!-- Stevens and Seifert, 2008):
	1507	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
	1508	( dr - 1.4E-3_wp ) ) )
	1509	!
	1510	!-- Slope parameter of gamma distribution (Seifert, 2008):
	1511	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	1512	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
[1012]	1513
[1361]	1514	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	1515	a_term - b_term * ( 1.0_wp + &
	1516	c_term / &
	1517	lambda_r )*( -1.0_wp &
	1518	( mu_r + 1.0_wp ) ) &
	1519	) &
[2232]	1520	) * flag
[1361]	1521
	1522	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	1523	a_term - b_term * ( 1.0_wp + &
	1524	c_term / &
	1525	lambda_r )*( -1.0_wp &
	1526	( mu_r + 4.0_wp ) ) &
	1527	) &
[2232]	1528	) * flag
[1361]	1529	ELSE
	1530	w_nr(k) = 0.0_wp
	1531	w_qr(k) = 0.0_wp
	1532	ENDIF
[1012]	1533	ENDDO
[1361]	1534	!
[2312]	1535	!-- Adjust boundary values using surface data type.
[2232]	1536	!-- Upward-facing
[2312]	1537	surf_s = bc_h(0)%start_index(j,i)
[2232]	1538	surf_e = bc_h(0)%end_index(j,i)
	1539	DO m = surf_s, surf_e
	1540	k = bc_h(0)%k(m)
	1541	w_nr(k-1) = w_nr(k)
	1542	w_qr(k-1) = w_qr(k)
	1543	ENDDO
	1544	!
	1545	!-- Downward-facing
[2312]	1546	surf_s = bc_h(1)%start_index(j,i)
[2232]	1547	surf_e = bc_h(1)%end_index(j,i)
	1548	DO m = surf_s, surf_e
	1549	k = bc_h(1)%k(m)
	1550	w_nr(k+1) = w_nr(k)
	1551	w_qr(k+1) = w_qr(k)
	1552	ENDDO
	1553	!
	1554	!-- Model top boundary value
[1361]	1555	w_nr(nzt+1) = 0.0_wp
	1556	w_qr(nzt+1) = 0.0_wp
	1557	!
	1558	!-- Compute Courant number
[2232]	1559	DO k = nzb+1, nzt
	1560	!
	1561	!-- Predetermine flag to mask topography
	1562	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1563
[1361]	1564	c_nr(k) = 0.25_wp * ( w_nr(k-1) + &
	1565	2.0_wp * w_nr(k) + w_nr(k+1) ) * &
[2232]	1566	dt_micro * ddzu(k) * flag
[1361]	1567	c_qr(k) = 0.25_wp * ( w_qr(k-1) + &
	1568	2.0_wp * w_qr(k) + w_qr(k+1) ) * &
[2232]	1569	dt_micro * ddzu(k) * flag
[2155]	1570	ENDDO
[1361]	1571	!
	1572	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
	1573	IF ( limiter_sedimentation ) THEN
	1574
[2232]	1575	DO k = nzb+1, nzt
	1576	!
	1577	!-- Predetermine flag to mask topography
	1578	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1579
[1646]	1580	d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
[1361]	1581	d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
	1582	d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
	1583
	1584	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	1585	2.0_wp * d_max, &
[2232]	1586	ABS( d_mean ) ) &
	1587	* flag
[1361]	1588
[1646]	1589	d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
[1361]	1590	d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
	1591	d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
	1592
	1593	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	1594	2.0_wp * d_max, &
	1595	ABS( d_mean ) )
	1596	ENDDO
	1597
	1598	ELSE
	1599
	1600	nr_slope = 0.0_wp
	1601	qr_slope = 0.0_wp
	1602
	1603	ENDIF
	1604
	1605	sed_nr(nzt+1) = 0.0_wp
	1606	sed_qr(nzt+1) = 0.0_wp
	1607	!
	1608	!-- Compute sedimentation flux
[2232]	1609	DO k = nzt, nzb+1, -1
[1361]	1610	!
[2232]	1611	!-- Predetermine flag to mask topography
	1612	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1613	!
[2312]	1614	!-- Sum up all rain drop number densities which contribute to the flux
[1361]	1615	!-- through k-1/2
	1616	flux = 0.0_wp
	1617	z_run = 0.0_wp ! height above z(k)
	1618	k_run = k
	1619	c_run = MIN( 1.0_wp, c_nr(k) )
	1620	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
	1621	flux = flux + hyrho(k_run) * &
	1622	( nr(k_run,j,i) + nr_slope(k_run) * &
[2232]	1623	( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run) &
	1624	* flag
	1625	z_run = z_run + dzu(k_run) * flag
	1626	k_run = k_run + 1 * flag
	1627	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) &
	1628	* flag
[1361]	1629	ENDDO
	1630	!
[2155]	1631	!-- It is not allowed to sediment more rain drop number density than
[1361]	1632	!-- available
	1633	flux = MIN( flux, &
	1634	hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * &
	1635	dt_micro &
	1636	)
	1637
[2232]	1638	sed_nr(k) = flux / dt_micro * flag
[1361]	1639	nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * &
[2232]	1640	ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]	1641	!
[2155]	1642	!-- Sum up all rain water content which contributes to the flux
[1361]	1643	!-- through k-1/2
	1644	flux = 0.0_wp
	1645	z_run = 0.0_wp ! height above z(k)
	1646	k_run = k
	1647	c_run = MIN( 1.0_wp, c_qr(k) )
	1648
	1649	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
	1650
	1651	flux = flux + hyrho(k_run) * ( qr(k_run,j,i) + &
	1652	qr_slope(k_run) * ( 1.0_wp - c_run ) * &
[2232]	1653	0.5_wp ) * c_run * dzu(k_run) * flag
	1654	z_run = z_run + dzu(k_run) * flag
	1655	k_run = k_run + 1 * flag
	1656	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) &
	1657	* flag
[1361]	1658
	1659	ENDDO
	1660	!
[2155]	1661	!-- It is not allowed to sediment more rain water content than
[1361]	1662	!-- available
	1663	flux = MIN( flux, &
	1664	hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * &
	1665	dt_micro &
	1666	)
	1667
[2232]	1668	sed_qr(k) = flux / dt_micro * flag
[1361]	1669
	1670	qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
[2232]	1671	ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]	1672	q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
[2232]	1673	ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]	1674	pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * &
	1675	ddzu(k+1) / hyrho(k) * l_d_cp * &
[2232]	1676	pt_d_t(k) * dt_micro * flag
[1361]	1677	!
	1678	!-- Compute the rain rate
	1679	IF ( call_microphysics_at_all_substeps ) THEN
[1691]	1680	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) &
[2232]	1681	* weight_substep(intermediate_timestep_count) &
	1682	* flag
[1361]	1683	ELSE
[2232]	1684	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
[1361]	1685	ENDIF
	1686
	1687	ENDDO
[1012]	1688	ENDDO
	1689	ENDDO
	1690
[1691]	1691	CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
	1692
	1693	END SUBROUTINE sedimentation_rain
	1694
	1695
	1696	!------------------------------------------------------------------------------!
	1697	! Description:
	1698	! ------------
	1699	!> Computation of the precipitation amount due to gravitational settling of
	1700	!> rain and cloud (fog) droplets
	1701	!------------------------------------------------------------------------------!
	1702	SUBROUTINE calc_precipitation_amount
	1703
[1849]	1704	USE arrays_3d, &
	1705	ONLY: precipitation_amount, prr
	1706
[1691]	1707	USE cloud_parameters, &
[1849]	1708	ONLY: hyrho
[1691]	1709
	1710	USE control_parameters, &
	1711	ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d, &
	1712	intermediate_timestep_count, intermediate_timestep_count_max,&
	1713	precipitation_amount_interval, time_do2d_xy
	1714
	1715	USE indices, &
[2232]	1716	ONLY: nxl, nxr, nys, nyn, nzb, nzt, wall_flags_0
[1691]	1717
	1718	USE kinds
	1719
[2232]	1720	USE surface_mod, &
	1721	ONLY : bc_h
	1722
[1691]	1723	IMPLICIT NONE
	1724
[2232]	1725	INTEGER(iwp) :: i !< running index x direction
	1726	INTEGER(iwp) :: j !< running index y direction
	1727	INTEGER(iwp) :: k !< running index y direction
	1728	INTEGER(iwp) :: m !< running index surface elements
	1729	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	1730	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1691]	1731
	1732	IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
	1733	( .NOT. call_microphysics_at_all_substeps .OR. &
	1734	intermediate_timestep_count == intermediate_timestep_count_max ) ) &
	1735	THEN
[2232]	1736	!
[2312]	1737	!-- Run over all upward-facing surface elements, i.e. non-natural,
[2232]	1738	!-- natural and urban
	1739	DO m = 1, bc_h(0)%ns
[2312]	1740	i = bc_h(0)%i(m)
[2232]	1741	j = bc_h(0)%j(m)
	1742	k = bc_h(0)%k(m)
	1743	precipitation_amount(j,i) = precipitation_amount(j,i) + &
	1744	prr(k,j,i) * hyrho(k) * dt_3d
	1745	ENDDO
[1691]	1746
[1361]	1747	ENDIF
	1748
[1691]	1749	END SUBROUTINE calc_precipitation_amount
[1361]	1750
[1012]	1751
[1000]	1752	!------------------------------------------------------------------------------!
[1682]	1753	! Description:
	1754	! ------------
[1849]	1755	!> Control of microphysics for grid points i,j
[1000]	1756	!------------------------------------------------------------------------------!
[1022]	1757
[1115]	1758	SUBROUTINE microphysics_control_ij( i, j )
	1759
[1320]	1760	USE arrays_3d, &
[2292]	1761	ONLY: hyp, nc, nr, pt, pt_init, prr, q, qc, qr, zu
[1115]	1762
[1320]	1763	USE cloud_parameters, &
[1849]	1764	ONLY: cp, hyrho, pt_d_t, r_d, t_d_pt
[1320]	1765
	1766	USE control_parameters, &
[1849]	1767	ONLY: call_microphysics_at_all_substeps, dt_3d, g, &
	1768	intermediate_timestep_count, large_scale_forcing, &
[2292]	1769	lsf_surf, microphysics_morrison, microphysics_seifert, &
	1770	microphysics_kessler, pt_surface, rho_surface, &
	1771	surface_pressure
[1320]	1772
	1773	USE indices, &
	1774	ONLY: nzb, nzt
	1775
	1776	USE kinds
	1777
	1778	USE statistics, &
	1779	ONLY: weight_pres
	1780
[1022]	1781	IMPLICIT NONE
	1782
[1682]	1783	INTEGER(iwp) :: i !<
	1784	INTEGER(iwp) :: j !<
	1785	INTEGER(iwp) :: k !<
[1115]	1786
[1682]	1787	REAL(wp) :: t_surface !<
[1320]	1788
[1361]	1789	IF ( large_scale_forcing .AND. lsf_surf ) THEN
[1241]	1790	!
	1791	!-- Calculate:
	1792	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
	1793	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
	1794	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
[1353]	1795	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
[1241]	1796	DO k = nzb, nzt+1
[1353]	1797	hyp(k) = surface_pressure * 100.0_wp * &
[1361]	1798	( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
[1353]	1799	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
	1800	t_d_pt(k) = 1.0_wp / pt_d_t(k)
[2155]	1801	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
[1241]	1802	ENDDO
	1803	!
	1804	!-- Compute reference density
[1353]	1805	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
[1241]	1806	ENDIF
	1807
[1361]	1808	!
[2155]	1809	!-- Compute length of time step
[1361]	1810	IF ( call_microphysics_at_all_substeps ) THEN
	1811	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
	1812	ELSE
	1813	dt_micro = dt_3d
	1814	ENDIF
[1241]	1815
[1115]	1816	!
[1361]	1817	!-- Use 1d arrays
[1115]	1818	q_1d(:) = q(:,j,i)
	1819	pt_1d(:) = pt(:,j,i)
	1820	qc_1d(:) = qc(:,j,i)
[1822]	1821	IF ( microphysics_seifert ) THEN
[1115]	1822	qr_1d(:) = qr(:,j,i)
	1823	nr_1d(:) = nr(:,j,i)
	1824	ENDIF
[2292]	1825	IF ( microphysics_morrison ) THEN
	1826	nc_1d(:) = nc(:,j,i)
	1827	ENDIF
[1361]	1828
[2292]	1829
[1115]	1830	!
[1822]	1831	!-- Reset precipitation rate
	1832	IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp
	1833
	1834	!
[1115]	1835	!-- Compute cloud physics
[1822]	1836	IF( microphysics_kessler ) THEN
	1837
	1838	CALL autoconversion_kessler( i,j )
[1831]	1839	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud( i,j )
[1822]	1840
	1841	ELSEIF ( microphysics_seifert ) THEN
	1842
	1843	CALL adjust_cloud( i,j )
[2292]	1844	IF ( microphysics_morrison ) CALL activation( i,j )
	1845	IF ( microphysics_morrison ) CALL condensation( i,j )
[1115]	1846	CALL autoconversion( i,j )
	1847	CALL accretion( i,j )
	1848	CALL selfcollection_breakup( i,j )
	1849	CALL evaporation_rain( i,j )
	1850	CALL sedimentation_rain( i,j )
[1831]	1851	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud( i,j )
[1115]	1852
[1691]	1853	ENDIF
	1854
[1822]	1855	CALL calc_precipitation_amount( i,j )
	1856
[1115]	1857	!
[1361]	1858	!-- Store results on the 3d arrays
	1859	q(:,j,i) = q_1d(:)
	1860	pt(:,j,i) = pt_1d(:)
[2292]	1861	IF ( microphysics_morrison ) THEN
	1862	qc(:,j,i) = qc_1d(:)
	1863	nc(:,j,i) = nc_1d(:)
	1864	ENDIF
[1822]	1865	IF ( microphysics_seifert ) THEN
[1361]	1866	qr(:,j,i) = qr_1d(:)
	1867	nr(:,j,i) = nr_1d(:)
[1115]	1868	ENDIF
	1869
	1870	END SUBROUTINE microphysics_control_ij
	1871
[1682]	1872	!------------------------------------------------------------------------------!
	1873	! Description:
	1874	! ------------
[2155]	1875	!> Adjust number of raindrops to avoid nonlinear effects in
[1682]	1876	!> sedimentation and evaporation of rain drops due to too small or
	1877	!> too big weights of rain drops (Stevens and Seifert, 2008).
	1878	!> The same procedure is applied to cloud droplets if they are determined
	1879	!> prognostically. Call for grid point i,j
	1880	!------------------------------------------------------------------------------!
[1115]	1881	SUBROUTINE adjust_cloud_ij( i, j )
	1882
[1320]	1883	USE cloud_parameters, &
[1849]	1884	ONLY: hyrho
[1320]	1885
[2292]	1886	USE control_parameters, &
	1887	ONLY: microphysics_morrison
	1888
[1320]	1889	USE indices, &
[2232]	1890	ONLY: nzb, nzt, wall_flags_0
[1320]	1891
	1892	USE kinds
	1893
[1115]	1894	IMPLICIT NONE
	1895
[1682]	1896	INTEGER(iwp) :: i !<
	1897	INTEGER(iwp) :: j !<
	1898	INTEGER(iwp) :: k !<
	1899
[2232]	1900	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1022]	1901
[2232]	1902	DO k = nzb+1, nzt
	1903	!
	1904	!-- Predetermine flag to mask topography
	1905	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1906
[1361]	1907	IF ( qr_1d(k) <= eps_sb ) THEN
	1908	qr_1d(k) = 0.0_wp
	1909	nr_1d(k) = 0.0_wp
[1065]	1910	ELSE
[1022]	1911	!
[2155]	1912	!-- Adjust number of raindrops to avoid nonlinear effects in
[1048]	1913	!-- sedimentation and evaporation of rain drops due to too small or
[1065]	1914	!-- too big weights of rain drops (Stevens and Seifert, 2008).
[1361]	1915	IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) ) THEN
[2232]	1916	nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin * flag
[1361]	1917	ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) ) THEN
[2232]	1918	nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax * flag
[1048]	1919	ENDIF
[1115]	1920
[1022]	1921	ENDIF
[1115]	1922
[2292]	1923	IF ( microphysics_morrison ) THEN
	1924	IF ( qc_1d(k) <= eps_sb ) THEN
	1925	qc_1d(k) = 0.0_wp
	1926	nc_1d(k) = 0.0_wp
	1927	ELSE
	1928	IF ( nc_1d(k) * xcmin > qc_1d(k) * hyrho(k) ) THEN
	1929	nc_1d(k) = qc_1d(k) * hyrho(k) / xamin * flag
	1930	ENDIF
	1931	ENDIF
	1932	ENDIF
	1933
[1022]	1934	ENDDO
	1935
[1115]	1936	END SUBROUTINE adjust_cloud_ij
[1022]	1937
[2292]	1938	!------------------------------------------------------------------------------!
	1939	! Description:
	1940	! ------------
	1941	!> Calculate number of activated condensation nucleii after simple activation
	1942	!> scheme of Twomey, 1959.
	1943	!------------------------------------------------------------------------------!
	1944	SUBROUTINE activation_ij( i, j )
[1106]	1945
[2292]	1946	USE arrays_3d, &
	1947	ONLY: hyp, nr, pt, q, qc, qr, nc
	1948
	1949	USE cloud_parameters, &
[2375]	1950	ONLY: hyrho, l_d_cp, l_d_r, l_v, molecular_weight_of_solute, &
	1951	molecular_weight_of_water, rho_l, rho_s, r_v, t_d_pt, &
	1952	vanthoff
[2292]	1953
	1954	USE constants, &
	1955	ONLY: pi
	1956
	1957	USE cpulog, &
	1958	ONLY: cpu_log, log_point_s
	1959
	1960	USE indices, &
	1961	ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt
	1962
	1963	USE kinds
	1964
	1965	USE control_parameters, &
	1966	ONLY: simulated_time
	1967
	1968	IMPLICIT NONE
	1969
	1970	INTEGER(iwp) :: i !<
	1971	INTEGER(iwp) :: j !<
	1972	INTEGER(iwp) :: k !<
	1973
	1974	REAL(wp) :: activ !<
	1975	REAL(wp) :: afactor !<
	1976	REAL(wp) :: alpha !<
	1977	REAL(wp) :: beta_act !<
	1978	REAL(wp) :: bfactor !<
	1979	REAL(wp) :: e_s !<
	1980	REAL(wp) :: k_act !<
	1981	REAL(wp) :: n_act !<
	1982	REAL(wp) :: n_ccn !<
	1983	REAL(wp) :: q_s !<
	1984	REAL(wp) :: s_0 !<
	1985	REAL(wp) :: sat !<
	1986	REAL(wp) :: sat_max !<
	1987	REAL(wp) :: sigma !<
[2375]	1988	REAL(wp) :: sigma_act !<
[2292]	1989	REAL(wp) :: t_int !<
	1990	REAL(wp) :: t_l !<
	1991
	1992	DO k = nzb+1, nzt
	1993	!
	1994	!-- Actual liquid water temperature:
	1995	t_l = t_d_pt(k) * pt_1d(k)
	1996
	1997	!
	1998	!-- Calculate actual temperature
	1999	t_int = pt_1d(k) * ( hyp(k) / 100000.0_wp )**0.286_wp
	2000	!
	2001	!-- Saturation vapor pressure at t_l:
	2002	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
	2003	( t_l - 35.86_wp ) &
	2004	)
	2005	!
[3026]	2006	!-- Computation of saturation mixing ratio:
	2007	q_s = 0.622_wp * e_s / ( hyp(k) - e_s )
[2292]	2008	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
	2009	q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / &
	2010	( 1.0_wp + alpha * q_s )
	2011
	2012	!-- Supersaturation:
	2013	sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
	2014
	2015	!
	2016	!-- Prescribe parameters for activation
	2017	!-- (see: Bott + Trautmann, 2002, Atm. Res., 64)
[2312]	2018	k_act = 0.7_wp
[2522]	2019	activ = 0.0_wp
[2292]	2020
	2021	IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
	2022	!
[2312]	2023	!-- Compute the number of activated Aerosols
[2292]	2024	!-- (see: Twomey, 1959, Pure and applied Geophysics, 43)
	2025	n_act = na_init * sat**k_act
	2026	!
	2027	!-- Compute the number of cloud droplets
	2028	!-- (see: Morrison + Grabowski, 2007, JAS, 64)
	2029	! activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro
	2030
	2031	!
	2032	!-- Compute activation rate after Khairoutdinov and Kogan
[2312]	2033	!-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
	2034	sat_max = 0.8_wp / 100.0_wp
	2035	activ = MAX( 0.0_wp, ( (na_init + nc_1d(k) ) * MIN &
	2036	( 1.0_wp, ( sat / sat_max )**k_act) - nc_1d(k) ) ) / &
[2292]	2037	dt_micro
	2038
	2039	nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)
	2040	ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk ) THEN
	2041	!
[2312]	2042	!-- Curvature effect (afactor) with surface tension
[2292]	2043	!-- parameterization by Straka (2009)
	2044	sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
	2045	afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
	2046	!
	2047	!-- Solute effect (bfactor)
	2048	bfactor = vanthoff * molecular_weight_of_water * &
	2049	rho_s / ( molecular_weight_of_solute * rho_l )
	2050
	2051	!
[2312]	2052	!-- Prescribe power index that describes the soluble fraction
[2375]	2053	!-- of an aerosol particle (beta).
[2292]	2054	!-- (see: Morrison + Grabowski, 2007, JAS, 64)
[2375]	2055	beta_act = 0.5_wp
	2056	sigma_act = sigma_bulk**( 1.0_wp + beta_act )
[2292]	2057	!
[2312]	2058	!-- Calculate mean geometric supersaturation (s_0) with
[2292]	2059	!-- parameterization by Khvorostyanov and Curry (2006)
[2375]	2060	s_0 = dry_aerosol_radius *(- ( 1.0_wp + beta_act ) ) &
[2292]	2061	( 4.0_wp * afactor*3 / ( 27.0_wp bfactor ) )**0.5_wp
	2062
	2063	!
[2312]	2064	!-- Calculate number of activated CCN as a function of
[2292]	2065	!-- supersaturation and taking Koehler theory into account
	2066	!-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
[2375]	2067	n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( &
	2068	LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
[2292]	2069	activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
[2312]	2070
	2071	nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)
[2292]	2072	ENDIF
	2073
	2074	ENDDO
	2075
	2076	END SUBROUTINE activation_ij
	2077
[1682]	2078	!------------------------------------------------------------------------------!
	2079	! Description:
	2080	! ------------
[2312]	2081	!> Calculate condensation rate for cloud water content (after Khairoutdinov and
[2292]	2082	!> Kogan, 2000).
	2083	!------------------------------------------------------------------------------!
	2084	SUBROUTINE condensation_ij( i, j )
	2085
	2086	USE arrays_3d, &
	2087	ONLY: hyp, nr, pt, q, qc, qr, nc
	2088
	2089	USE cloud_parameters, &
	2090	ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
	2091
	2092	USE constants, &
	2093	ONLY: pi
	2094
	2095	USE cpulog, &
	2096	ONLY: cpu_log, log_point_s
	2097
	2098	USE indices, &
	2099	ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt
	2100
	2101	USE kinds
	2102
	2103	USE control_parameters, &
	2104	ONLY: simulated_time
	2105
	2106	IMPLICIT NONE
	2107
	2108	INTEGER(iwp) :: i !<
	2109	INTEGER(iwp) :: j !<
	2110	INTEGER(iwp) :: k !<
	2111
	2112	REAL(wp) :: alpha !<
	2113	REAL(wp) :: cond !<
	2114	REAL(wp) :: cond_max !<
	2115	REAL(wp) :: dc !<
	2116	REAL(wp) :: e_s !<
	2117	REAL(wp) :: evap !<
	2118	REAL(wp) :: evap_nc !<
	2119	REAL(wp) :: g_fac !<
	2120	REAL(wp) :: nc_0 !<
	2121	REAL(wp) :: q_s !<
	2122	REAL(wp) :: sat !<
	2123	REAL(wp) :: t_l !<
	2124	REAL(wp) :: temp !<
	2125	REAL(wp) :: xc !<
	2126
	2127
	2128	DO k = nzb+1, nzt
	2129	!
	2130	!-- Actual liquid water temperature:
	2131	t_l = t_d_pt(k) * pt_1d(k)
	2132	!
	2133	!-- Saturation vapor pressure at t_l:
	2134	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
	2135	( t_l - 35.86_wp ) &
	2136	)
	2137	!
[3026]	2138	!-- Computation of saturation mixing ratio:
	2139	q_s = 0.622_wp * e_s / ( hyp(k) - e_s )
[2292]	2140	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
	2141	q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / &
	2142	( 1.0_wp + alpha * q_s )
	2143
	2144	!-- Supersaturation:
	2145	sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
	2146
	2147
	2148	!
	2149	!-- Actual temperature:
	2150	temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
	2151
	2152	g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
	2153	l_v / ( thermal_conductivity_l * temp ) &
	2154	+ r_v * temp / ( diff_coeff_l * e_s ) &
	2155	)
	2156	!
	2157	!-- Mean weight of cloud drops
	2158	IF ( nc_1d(k) <= 0.0_wp) CYCLE
[2312]	2159	xc = MAX( (hyrho(k) * qc_1d(k) / nc_1d(k)), xcmin)
[2292]	2160	!
	2161	!-- Weight averaged diameter of cloud drops:
	2162	dc = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
	2163	!
	2164	!-- Integral diameter of cloud drops
[2312]	2165	nc_0 = nc_1d(k) * dc
[2292]	2166	!
	2167	!-- Condensation needs only to be calculated in supersaturated regions
	2168	IF ( sat > 0.0_wp ) THEN
	2169	!
	2170	!-- Condensation rate of cloud water content
	2171	!-- after KK scheme.
	2172	!-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
	2173	cond = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
	2174	cond_max = q_1d(k) - q_s - qc_1d(k) - qr_1d(k)
	2175	cond = MIN( cond, cond_max / dt_micro )
[2312]	2176
[2292]	2177	qc_1d(k) = qc_1d(k) + cond * dt_micro
	2178	ELSEIF ( sat < 0.0_wp ) THEN
	2179	evap = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
	2180	evap = MAX( evap, -qc_1d(k) / dt_micro )
	2181
	2182	qc_1d(k) = qc_1d(k) + evap * dt_micro
	2183	ENDIF
	2184	ENDDO
	2185
	2186	END SUBROUTINE condensation_ij
	2187
	2188
	2189	!------------------------------------------------------------------------------!
	2190	! Description:
	2191	! ------------
[1682]	2192	!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
	2193	!------------------------------------------------------------------------------!
[1005]	2194	SUBROUTINE autoconversion_ij( i, j )
[1000]	2195
[1320]	2196	USE arrays_3d, &
[1849]	2197	ONLY: diss, dzu
[1115]	2198
[1320]	2199	USE cloud_parameters, &
[1849]	2200	ONLY: hyrho
[1320]	2201
	2202	USE control_parameters, &
[2292]	2203	ONLY: microphysics_morrison, rho_surface
[1320]	2204
	2205	USE grid_variables, &
	2206	ONLY: dx, dy
	2207
	2208	USE indices, &
[2232]	2209	ONLY: nzb, nzt, wall_flags_0
[1320]	2210
	2211	USE kinds
	2212
[1000]	2213	IMPLICIT NONE
	2214
[1682]	2215	INTEGER(iwp) :: i !<
	2216	INTEGER(iwp) :: j !<
	2217	INTEGER(iwp) :: k !<
[1000]	2218
[2155]	2219	REAL(wp) :: alpha_cc !<
[1682]	2220	REAL(wp) :: autocon !<
	2221	REAL(wp) :: dissipation !<
[2232]	2222	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1682]	2223	REAL(wp) :: k_au !<
	2224	REAL(wp) :: l_mix !<
[2292]	2225	REAL(wp) :: nc_auto !<
[1682]	2226	REAL(wp) :: nu_c !<
	2227	REAL(wp) :: phi_au !<
	2228	REAL(wp) :: r_cc !<
	2229	REAL(wp) :: rc !<
	2230	REAL(wp) :: re_lambda !<
	2231	REAL(wp) :: sigma_cc !<
	2232	REAL(wp) :: tau_cloud !<
	2233	REAL(wp) :: xc !<
[1106]	2234
[2232]	2235	DO k = nzb+1, nzt
	2236	!
	2237	!-- Predetermine flag to mask topography
	2238	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[2292]	2239	nc_auto = MERGE ( nc_1d(k), nc_const, microphysics_morrison )
[1000]	2240
[2522]	2241	IF ( qc_1d(k) > eps_sb .AND. nc_auto > eps_mr ) THEN
[1361]	2242
	2243	k_au = k_cc / ( 20.0_wp * x0 )
[1012]	2244	!
[1048]	2245	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[1353]	2246	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
[2522]	2247	tau_cloud = MAX( 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) ), &
	2248	0.0_wp )
[1012]	2249	!
[2155]	2250	!-- Universal function for autoconversion process
[1012]	2251	!-- (Seifert and Beheng, 2006):
[1361]	2252	phi_au = 600.0_wp * tau_cloud*0.68_wp ( 1.0_wp - tau_cloud0.68_wp )3
[1012]	2253	!
	2254	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
[1353]	2255	!-- (Use constant nu_c = 1.0_wp instead?)
[1361]	2256	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp )
[1012]	2257	!
	2258	!-- Mean weight of cloud droplets:
[2292]	2259	xc = hyrho(k) * qc_1d(k) / nc_auto
[1012]	2260	!
[2155]	2261	!-- Parameterized turbulence effects on autoconversion (Seifert,
[1065]	2262	!-- Nuijens and Stevens, 2010)
[1831]	2263	IF ( collision_turbulence ) THEN
[1065]	2264	!
	2265	!-- Weight averaged radius of cloud droplets:
[1353]	2266	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
[1065]	2267
[1353]	2268	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
	2269	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
	2270	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
[1065]	2271	!
	2272	!-- Mixing length (neglecting distance to ground and stratification)
[1334]	2273	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
[1065]	2274	!
[2155]	2275	!-- Limit dissipation rate according to Seifert, Nuijens and
[1065]	2276	!-- Stevens (2010)
[1361]	2277	dissipation = MIN( 0.06_wp, diss(k,j,i) )
[1065]	2278	!
	2279	!-- Compute Taylor-microscale Reynolds number:
[1361]	2280	re_lambda = 6.0_wp / 11.0_wp * &
	2281	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
	2282	SQRT( 15.0_wp / kin_vis_air ) * &
	2283	dissipation**( 1.0_wp / 6.0_wp )
[1065]	2284	!
	2285	!-- The factor of 1.0E4 is needed to convert the dissipation rate
	2286	!-- from m2 s-3 to cm2 s-3.
[1361]	2287	k_au = k_au * ( 1.0_wp + &
	2288	dissipation * 1.0E4_wp * &
	2289	( re_lambda * 1.0E-3_wp )*0.25_wp &
	2290	( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) / &
	2291	sigma_cc )**2 &
	2292	) + beta_cc &
	2293	) &
	2294	)
[1065]	2295	ENDIF
	2296	!
[1012]	2297	!-- Autoconversion rate (Seifert and Beheng, 2006):
[1361]	2298	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
	2299	( nu_c + 1.0_wp )*2 qc_1d(k)*2 xc*2 &
	2300	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
[1115]	2301	rho_surface
	2302	autocon = MIN( autocon, qc_1d(k) / dt_micro )
[1106]	2303
[2232]	2304	qr_1d(k) = qr_1d(k) + autocon * dt_micro * flag
	2305	qc_1d(k) = qc_1d(k) - autocon * dt_micro * flag
	2306	nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro * flag
[2292]	2307	IF ( microphysics_morrison ) THEN
[2312]	2308	nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), 2.0_wp * &
[2292]	2309	autocon / x0 * hyrho(k) * dt_micro * flag )
	2310	ENDIF
[1115]	2311
[1005]	2312	ENDIF
[1000]	2313
	2314	ENDDO
	2315
[1005]	2316	END SUBROUTINE autoconversion_ij
	2317
[1822]	2318	!------------------------------------------------------------------------------!
	2319	! Description:
	2320	! ------------
	2321	!> Autoconversion process (Kessler, 1969).
	2322	!------------------------------------------------------------------------------!
	2323	SUBROUTINE autoconversion_kessler_ij( i, j )
[1106]	2324
[1822]	2325	USE arrays_3d, &
[1849]	2326	ONLY: dzw, prr
[1822]	2327
	2328	USE cloud_parameters, &
[1849]	2329	ONLY: l_d_cp, pt_d_t
[1822]	2330
	2331	USE indices, &
[2317]	2332	ONLY: nzb, nzt, wall_flags_0
[1822]	2333
	2334	USE kinds
	2335
[2317]	2336	USE surface_mod, &
[2698]	2337	ONLY: get_topography_top_index_ji
[1822]	2338
[2317]	2339
[1822]	2340	IMPLICIT NONE
	2341
[2232]	2342	INTEGER(iwp) :: i !<
	2343	INTEGER(iwp) :: j !<
	2344	INTEGER(iwp) :: k !<
	2345	INTEGER(iwp) :: k_wall !< topography top index
[1822]	2346
	2347	REAL(wp) :: dqdt_precip !<
[2232]	2348	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1822]	2349
[2232]	2350	!
	2351	!-- Determine vertical index of topography top
[2698]	2352	k_wall = get_topography_top_index_ji( j, i, 's' )
[2232]	2353	DO k = nzb+1, nzt
	2354	!
	2355	!-- Predetermine flag to mask topography
	2356	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1822]	2357
	2358	IF ( qc_1d(k) > ql_crit ) THEN
	2359	dqdt_precip = prec_time_const * ( qc_1d(k) - ql_crit )
	2360	ELSE
	2361	dqdt_precip = 0.0_wp
	2362	ENDIF
	2363
[2232]	2364	qc_1d(k) = qc_1d(k) - dqdt_precip * dt_micro * flag
	2365	q_1d(k) = q_1d(k) - dqdt_precip * dt_micro * flag
	2366	pt_1d(k) = pt_1d(k) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k) * flag
[1822]	2367
	2368	!
[1845]	2369	!-- Compute the rain rate (stored on surface grid point)
[2232]	2370	prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
[1822]	2371
	2372	ENDDO
	2373
	2374	END SUBROUTINE autoconversion_kessler_ij
	2375
[1682]	2376	!------------------------------------------------------------------------------!
	2377	! Description:
	2378	! ------------
	2379	!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
	2380	!------------------------------------------------------------------------------!
[1005]	2381	SUBROUTINE accretion_ij( i, j )
	2382
[1320]	2383	USE arrays_3d, &
[1849]	2384	ONLY: diss
[1115]	2385
[1320]	2386	USE cloud_parameters, &
[1849]	2387	ONLY: hyrho
[1320]	2388
	2389	USE control_parameters, &
[2292]	2390	ONLY: microphysics_morrison, rho_surface
[1320]	2391
	2392	USE indices, &
[2232]	2393	ONLY: nzb, nzt, wall_flags_0
[1320]	2394
	2395	USE kinds
	2396
[1005]	2397	IMPLICIT NONE
	2398
[1682]	2399	INTEGER(iwp) :: i !<
	2400	INTEGER(iwp) :: j !<
	2401	INTEGER(iwp) :: k !<
[1005]	2402
[1682]	2403	REAL(wp) :: accr !<
[2232]	2404	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1682]	2405	REAL(wp) :: k_cr !<
[2292]	2406	REAL(wp) :: nc_accr !<
[1682]	2407	REAL(wp) :: phi_ac !<
	2408	REAL(wp) :: tau_cloud !<
[2292]	2409	REAL(wp) :: xc !<
[1320]	2410
[2292]	2411
[2232]	2412	DO k = nzb+1, nzt
	2413	!
	2414	!-- Predetermine flag to mask topography
	2415	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[2292]	2416	nc_accr = MERGE ( nc_1d(k), nc_const, microphysics_morrison )
[2232]	2417
[2292]	2418	IF ( ( qc_1d(k) > eps_sb ) .AND. ( qr_1d(k) > eps_sb ) .AND. &
	2419	( nc_accr > eps_mr ) ) THEN
[1012]	2420	!
[1048]	2421	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[2155]	2422	tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
[1012]	2423	!
[2155]	2424	!-- Universal function for accretion process
[1048]	2425	!-- (Seifert and Beheng, 2001):
[1361]	2426	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
[2292]	2427
[1012]	2428	!
[2292]	2429	!-- Mean weight of cloud drops
[2312]	2430	xc = MAX( (hyrho(k) * qc_1d(k) / nc_accr), xcmin)
[2292]	2431	!
[2155]	2432	!-- Parameterized turbulence effects on autoconversion (Seifert,
	2433	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
[1361]	2434	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
[1831]	2435	IF ( collision_turbulence ) THEN
[1361]	2436	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
	2437	MIN( 600.0_wp, &
	2438	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
	2439	)
[1065]	2440	ELSE
[2155]	2441	k_cr = k_cr0
[1065]	2442	ENDIF
	2443	!
[1012]	2444	!-- Accretion rate (Seifert and Beheng, 2006):
[2292]	2445	accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * &
	2446	SQRT( rho_surface * hyrho(k) )
[1115]	2447	accr = MIN( accr, qc_1d(k) / dt_micro )
[1106]	2448
[2232]	2449	qr_1d(k) = qr_1d(k) + accr * dt_micro * flag
	2450	qc_1d(k) = qc_1d(k) - accr * dt_micro * flag
[2292]	2451	IF ( microphysics_morrison ) THEN
	2452	nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), accr / xc * &
	2453	hyrho(k) * dt_micro * flag &
	2454	)
	2455	ENDIF
[1115]	2456
[2292]	2457
[1005]	2458	ENDIF
[1106]	2459
[1005]	2460	ENDDO
	2461
[1000]	2462	END SUBROUTINE accretion_ij
	2463
[1005]	2464
[1682]	2465	!------------------------------------------------------------------------------!
	2466	! Description:
	2467	! ------------
	2468	!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
	2469	!------------------------------------------------------------------------------!
[1005]	2470	SUBROUTINE selfcollection_breakup_ij( i, j )
	2471
[1320]	2472	USE cloud_parameters, &
[1849]	2473	ONLY: hyrho
[1320]	2474
	2475	USE control_parameters, &
[1849]	2476	ONLY: rho_surface
[1320]	2477
	2478	USE indices, &
[2232]	2479	ONLY: nzb, nzt, wall_flags_0
[1320]	2480
	2481	USE kinds
[2155]	2482
[1005]	2483	IMPLICIT NONE
	2484
[1682]	2485	INTEGER(iwp) :: i !<
	2486	INTEGER(iwp) :: j !<
	2487	INTEGER(iwp) :: k !<
[1005]	2488
[1682]	2489	REAL(wp) :: breakup !<
	2490	REAL(wp) :: dr !<
[2232]	2491	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1682]	2492	REAL(wp) :: phi_br !<
	2493	REAL(wp) :: selfcoll !<
[1320]	2494
[2232]	2495	DO k = nzb+1, nzt
	2496	!
	2497	!-- Predetermine flag to mask topography
	2498	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2499
[1115]	2500	IF ( qr_1d(k) > eps_sb ) THEN
[1012]	2501	!
[1115]	2502	!-- Selfcollection rate (Seifert and Beheng, 2001):
[1361]	2503	selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
[1012]	2504	!
[1115]	2505	!-- Weight averaged diameter of rain drops:
[1334]	2506	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	2507	!
[1048]	2508	!-- Collisional breakup rate (Seifert, 2008):
[1353]	2509	IF ( dr >= 0.3E-3_wp ) THEN
	2510	phi_br = k_br * ( dr - 1.1E-3_wp )
	2511	breakup = selfcoll * ( phi_br + 1.0_wp )
[1005]	2512	ELSE
[1353]	2513	breakup = 0.0_wp
[1005]	2514	ENDIF
[1048]	2515
[1115]	2516	selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
[2232]	2517	nr_1d(k) = nr_1d(k) + selfcoll * dt_micro * flag
[1106]	2518
[2155]	2519	ENDIF
[1005]	2520	ENDDO
	2521
	2522	END SUBROUTINE selfcollection_breakup_ij
	2523
[1106]	2524
[1682]	2525	!------------------------------------------------------------------------------!
	2526	! Description:
	2527	! ------------
[2155]	2528	!> Evaporation of precipitable water. Condensation is neglected for
[1682]	2529	!> precipitable water. Call for grid point i,j
	2530	!------------------------------------------------------------------------------!
[1012]	2531	SUBROUTINE evaporation_rain_ij( i, j )
	2532
[1320]	2533	USE arrays_3d, &
[1849]	2534	ONLY: hyp
[1048]	2535
[1320]	2536	USE cloud_parameters, &
[1849]	2537	ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
[1320]	2538
	2539	USE constants, &
	2540	ONLY: pi
	2541
	2542	USE indices, &
[2232]	2543	ONLY: nzb, nzt, wall_flags_0
[1320]	2544
	2545	USE kinds
	2546
[1012]	2547	IMPLICIT NONE
	2548
[1682]	2549	INTEGER(iwp) :: i !<
	2550	INTEGER(iwp) :: j !<
	2551	INTEGER(iwp) :: k !<
[1012]	2552
[1682]	2553	REAL(wp) :: alpha !<
	2554	REAL(wp) :: dr !<
	2555	REAL(wp) :: e_s !<
	2556	REAL(wp) :: evap !<
	2557	REAL(wp) :: evap_nr !<
	2558	REAL(wp) :: f_vent !<
[2232]	2559	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1682]	2560	REAL(wp) :: g_evap !<
	2561	REAL(wp) :: lambda_r !<
	2562	REAL(wp) :: mu_r !<
	2563	REAL(wp) :: mu_r_2 !<
	2564	REAL(wp) :: mu_r_5d2 !<
	2565	REAL(wp) :: nr_0 !<
	2566	REAL(wp) :: q_s !<
	2567	REAL(wp) :: sat !<
	2568	REAL(wp) :: t_l !<
	2569	REAL(wp) :: temp !<
	2570	REAL(wp) :: xr !<
[1320]	2571
[2232]	2572	DO k = nzb+1, nzt
	2573	!
	2574	!-- Predetermine flag to mask topography
	2575	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2576
[1115]	2577	IF ( qr_1d(k) > eps_sb ) THEN
[1012]	2578	!
	2579	!-- Actual liquid water temperature:
[1115]	2580	t_l = t_d_pt(k) * pt_1d(k)
[1012]	2581	!
	2582	!-- Saturation vapor pressure at t_l:
[1361]	2583	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
	2584	( t_l - 35.86_wp ) &
	2585	)
[1012]	2586	!
[3026]	2587	!-- Computation of saturation mixing ratio:
	2588	q_s = 0.622_wp * e_s / ( hyp(k) - e_s )
[1353]	2589	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
[1361]	2590	q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
[1012]	2591	!
[1106]	2592	!-- Supersaturation:
[1361]	2593	sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
[1012]	2594	!
[1361]	2595	!-- Evaporation needs only to be calculated in subsaturated regions
	2596	IF ( sat < 0.0_wp ) THEN
[1012]	2597	!
[1361]	2598	!-- Actual temperature:
	2599	temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
[2155]	2600
[1361]	2601	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v / &
	2602	( thermal_conductivity_l * temp ) + &
	2603	r_v * temp / ( diff_coeff_l * e_s ) &
	2604	)
[1012]	2605	!
[1361]	2606	!-- Mean weight of rain drops
	2607	xr = hyrho(k) * qr_1d(k) / nr_1d(k)
[1115]	2608	!
[1361]	2609	!-- Weight averaged diameter of rain drops:
	2610	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	2611	!
[2155]	2612	!-- Compute ventilation factor and intercept parameter
[1361]	2613	!-- (Seifert and Beheng, 2006; Seifert, 2008):
	2614	IF ( ventilation_effect ) THEN
[1115]	2615	!
[1361]	2616	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	2617	!-- Stevens and Seifert, 2008):
	2618	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
	2619	!
	2620	!-- Slope parameter of gamma distribution (Seifert, 2008):
	2621	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	2622	( mu_r + 1.0_wp ) &
	2623	)**( 1.0_wp / 3.0_wp ) / dr
[1115]	2624
[1361]	2625	mu_r_2 = mu_r + 2.0_wp
[2155]	2626	mu_r_5d2 = mu_r + 2.5_wp
[1361]	2627
[2155]	2628	f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) + &
[1361]	2629	b_vent * schmidt_p_1d3 * &
	2630	SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * &
	2631	lambda_r*( -mu_r_5d2 ) &
	2632	( 1.0_wp - &
	2633	0.5_wp * ( b_term / a_term ) * &
	2634	( lambda_r / ( c_term + lambda_r ) &
	2635	)**mu_r_5d2 - &
	2636	0.125_wp * ( b_term / a_term )*2 &
	2637	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
	2638	)**mu_r_5d2 - &
	2639	0.0625_wp * ( b_term / a_term )*3 &
	2640	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
	2641	)**mu_r_5d2 - &
[2155]	2642	0.0390625_wp * ( b_term / a_term )*4 &
[1361]	2643	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
	2644	)**mu_r_5d2 &
	2645	)
	2646
	2647	nr_0 = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) / &
[2155]	2648	gamm( mu_r + 1.0_wp )
[1361]	2649	ELSE
	2650	f_vent = 1.0_wp
	2651	nr_0 = nr_1d(k) * dr
	2652	ENDIF
[1012]	2653	!
[1361]	2654	!-- Evaporation rate of rain water content (Seifert and Beheng, 2006):
	2655	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
	2656	evap = MAX( evap, -qr_1d(k) / dt_micro )
	2657	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
	2658	-nr_1d(k) / dt_micro )
[1106]	2659
[2232]	2660	qr_1d(k) = qr_1d(k) + evap * dt_micro * flag
	2661	nr_1d(k) = nr_1d(k) + evap_nr * dt_micro * flag
[1115]	2662
[1361]	2663	ENDIF
[2155]	2664	ENDIF
[1106]	2665
[1012]	2666	ENDDO
	2667
	2668	END SUBROUTINE evaporation_rain_ij
	2669
[1106]	2670
[1682]	2671	!------------------------------------------------------------------------------!
	2672	! Description:
	2673	! ------------
[2155]	2674	!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
[1682]	2675	!> Call for grid point i,j
	2676	!------------------------------------------------------------------------------!
[1012]	2677	SUBROUTINE sedimentation_cloud_ij( i, j )
	2678
[1320]	2679	USE arrays_3d, &
[1849]	2680	ONLY: ddzu, dzu, prr
[1320]	2681
	2682	USE cloud_parameters, &
[1849]	2683	ONLY: hyrho, l_d_cp, pt_d_t
[1320]	2684
	2685	USE control_parameters, &
[2292]	2686	ONLY: call_microphysics_at_all_substeps, &
	2687	intermediate_timestep_count, microphysics_morrison
[1320]	2688
	2689	USE indices, &
[2232]	2690	ONLY: nzb, nzb, nzt, wall_flags_0
[1320]	2691
	2692	USE kinds
[2155]	2693
[1691]	2694	USE statistics, &
	2695	ONLY: weight_substep
	2696
[1012]	2697	IMPLICIT NONE
	2698
[1849]	2699	INTEGER(iwp) :: i !<
	2700	INTEGER(iwp) :: j !<
	2701	INTEGER(iwp) :: k !<
[1106]	2702
[2292]	2703	REAL(wp) :: flag !< flag to indicate first grid level above surface
	2704	REAL(wp) :: nc_sedi !<
	2705
	2706	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nc !<
[2232]	2707	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
[1115]	2708
[1353]	2709	sed_qc(nzt+1) = 0.0_wp
[2375]	2710	sed_nc(nzt+1) = 0.0_wp
[1012]	2711
[2375]	2712
[2232]	2713	DO k = nzt, nzb+1, -1
	2714	!
	2715	!-- Predetermine flag to mask topography
	2716	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[2292]	2717	nc_sedi = MERGE( nc_1d(k), nc_const, microphysics_morrison )
	2718	!
	2719	!-- Sedimentation fluxes for number concentration are only calculated
	2720	!-- for cloud_scheme = 'morrison'
	2721	IF ( microphysics_morrison ) THEN
	2722	IF ( qc_1d(k) > eps_sb .AND. nc_1d(k) > eps_mr ) THEN
[2312]	2723	sed_nc(k) = sed_qc_const * &
[2292]	2724	( qc_1d(k) * hyrho(k) )*( 2.0_wp / 3.0_wp ) &
	2725	( nc_1d(k) )**( 1.0_wp / 3.0_wp )
	2726	ELSE
	2727	sed_nc(k) = 0.0_wp
	2728	ENDIF
[2232]	2729
[2292]	2730	sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) * &
	2731	nc_1d(k) / dt_micro + sed_nc(k+1) &
	2732	) * flag
	2733
	2734	nc_1d(k) = nc_1d(k) + ( sed_nc(k+1) - sed_nc(k) ) * &
	2735	ddzu(k+1) / hyrho(k) * dt_micro * flag
	2736	ENDIF
	2737
[2375]	2738	IF ( qc_1d(k) > eps_sb .AND. nc_sedi > eps_mr ) THEN
[2292]	2739	sed_qc(k) = sed_qc_const * nc_sedi*( -2.0_wp / 3.0_wp ) &
[2232]	2740	( qc_1d(k) * hyrho(k) )*( 5.0_wp / 3.0_wp ) flag
[1115]	2741	ELSE
[1353]	2742	sed_qc(k) = 0.0_wp
[1012]	2743	ENDIF
[1115]	2744
[1361]	2745	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) / &
	2746	dt_micro + sed_qc(k+1) &
[2232]	2747	) * flag
[1115]	2748
[1361]	2749	q_1d(k) = q_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
[2232]	2750	hyrho(k) * dt_micro * flag
[2155]	2751	qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
[2232]	2752	hyrho(k) * dt_micro * flag
[1361]	2753	pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
[2232]	2754	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro * flag
[1115]	2755
[1691]	2756	!
	2757	!-- Compute the precipitation rate of cloud (fog) droplets
[1822]	2758	IF ( call_microphysics_at_all_substeps ) THEN
[2232]	2759	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * &
	2760	weight_substep(intermediate_timestep_count) * flag
[1822]	2761	ELSE
[2232]	2762	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
[1691]	2763	ENDIF
	2764
[1012]	2765	ENDDO
	2766
	2767	END SUBROUTINE sedimentation_cloud_ij
	2768
[1106]	2769
[1682]	2770	!------------------------------------------------------------------------------!
	2771	! Description:
	2772	! ------------
	2773	!> Computation of sedimentation flux. Implementation according to Stevens
	2774	!> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j
	2775	!------------------------------------------------------------------------------!
[1012]	2776	SUBROUTINE sedimentation_rain_ij( i, j )
	2777
[1320]	2778	USE arrays_3d, &
[1849]	2779	ONLY: ddzu, dzu, prr
[1320]	2780
	2781	USE cloud_parameters, &
[1849]	2782	ONLY: hyrho, l_d_cp, pt_d_t
[1320]	2783
	2784	USE control_parameters, &
[1849]	2785	ONLY: call_microphysics_at_all_substeps, intermediate_timestep_count
[1320]	2786
	2787	USE indices, &
[2232]	2788	ONLY: nzb, nzb, nzt, wall_flags_0
[1320]	2789
	2790	USE kinds
	2791
	2792	USE statistics, &
	2793	ONLY: weight_substep
[2155]	2794
[2232]	2795	USE surface_mod, &
	2796	ONLY : bc_h
[2312]	2797
[1012]	2798	IMPLICIT NONE
	2799
[2232]	2800	INTEGER(iwp) :: i !< running index x direction
	2801	INTEGER(iwp) :: j !< running index y direction
	2802	INTEGER(iwp) :: k !< running index z direction
	2803	INTEGER(iwp) :: k_run !<
	2804	INTEGER(iwp) :: m !< running index surface elements
	2805	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	2806	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1012]	2807
[1682]	2808	REAL(wp) :: c_run !<
	2809	REAL(wp) :: d_max !<
	2810	REAL(wp) :: d_mean !<
	2811	REAL(wp) :: d_min !<
	2812	REAL(wp) :: dr !<
	2813	REAL(wp) :: flux !<
[2232]	2814	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1682]	2815	REAL(wp) :: lambda_r !<
	2816	REAL(wp) :: mu_r !<
	2817	REAL(wp) :: z_run !<
[1320]	2818
[1682]	2819	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<
	2820	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<
	2821	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<
	2822	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<
	2823	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<
	2824	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<
	2825	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<
	2826	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<
[1320]	2827
[1012]	2828	!
[2155]	2829	!-- Compute velocities
[2232]	2830	DO k = nzb+1, nzt
	2831	!
	2832	!-- Predetermine flag to mask topography
	2833	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2834
[1115]	2835	IF ( qr_1d(k) > eps_sb ) THEN
	2836	!
	2837	!-- Weight averaged diameter of rain drops:
[1334]	2838	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	2839	!
	2840	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	2841	!-- Stevens and Seifert, 2008):
[1353]	2842	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
[1115]	2843	!
	2844	!-- Slope parameter of gamma distribution (Seifert, 2008):
[1361]	2845	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	2846	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
[1115]	2847
[1361]	2848	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	2849	a_term - b_term * ( 1.0_wp + &
	2850	c_term / lambda_r )*( -1.0_wp &
	2851	( mu_r + 1.0_wp ) ) &
	2852	) &
[2232]	2853	) * flag
[1361]	2854	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	2855	a_term - b_term * ( 1.0_wp + &
	2856	c_term / lambda_r )*( -1.0_wp &
	2857	( mu_r + 4.0_wp ) ) &
	2858	) &
[2232]	2859	) * flag
[1065]	2860	ELSE
[1353]	2861	w_nr(k) = 0.0_wp
	2862	w_qr(k) = 0.0_wp
[1065]	2863	ENDIF
	2864	ENDDO
[1048]	2865	!
[2312]	2866	!-- Adjust boundary values using surface data type.
[2232]	2867	!-- Upward facing non-natural
[2312]	2868	surf_s = bc_h(0)%start_index(j,i)
[2232]	2869	surf_e = bc_h(0)%end_index(j,i)
	2870	DO m = surf_s, surf_e
	2871	k = bc_h(0)%k(m)
	2872	w_nr(k-1) = w_nr(k)
	2873	w_qr(k-1) = w_qr(k)
	2874	ENDDO
	2875	!
	2876	!-- Downward facing non-natural
[2312]	2877	surf_s = bc_h(1)%start_index(j,i)
[2232]	2878	surf_e = bc_h(1)%end_index(j,i)
	2879	DO m = surf_s, surf_e
	2880	k = bc_h(1)%k(m)
	2881	w_nr(k+1) = w_nr(k)
	2882	w_qr(k+1) = w_qr(k)
	2883	ENDDO
	2884	!
	2885	!-- Neumann boundary condition at model top
[1353]	2886	w_nr(nzt+1) = 0.0_wp
	2887	w_qr(nzt+1) = 0.0_wp
[1065]	2888	!
	2889	!-- Compute Courant number
[2232]	2890	DO k = nzb+1, nzt
	2891	!
	2892	!-- Predetermine flag to mask topography
	2893	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2894
[1361]	2895	c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) * &
[2232]	2896	dt_micro * ddzu(k) * flag
[1361]	2897	c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * &
[2232]	2898	dt_micro * ddzu(k) * flag
[2155]	2899	ENDDO
[1065]	2900	!
	2901	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
	2902	IF ( limiter_sedimentation ) THEN
	2903
[2232]	2904	DO k = nzb+1, nzt
	2905	!
	2906	!-- Predetermine flag to mask topography
	2907	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2908
[1646]	2909	d_mean = 0.5_wp * ( qr_1d(k+1) - qr_1d(k-1) )
[1115]	2910	d_min = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
	2911	d_max = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
[1065]	2912
[1361]	2913	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	2914	2.0_wp * d_max, &
[2232]	2915	ABS( d_mean ) ) * flag
[1065]	2916
[1646]	2917	d_mean = 0.5_wp * ( nr_1d(k+1) - nr_1d(k-1) )
[1115]	2918	d_min = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
	2919	d_max = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
[1065]	2920
[1361]	2921	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	2922	2.0_wp * d_max, &
[2232]	2923	ABS( d_mean ) ) * flag
[1022]	2924	ENDDO
[1048]	2925
[1065]	2926	ELSE
[1106]	2927
[1353]	2928	nr_slope = 0.0_wp
	2929	qr_slope = 0.0_wp
[1106]	2930
[1065]	2931	ENDIF
[1115]	2932
[1353]	2933	sed_nr(nzt+1) = 0.0_wp
	2934	sed_qr(nzt+1) = 0.0_wp
[1065]	2935	!
	2936	!-- Compute sedimentation flux
[2232]	2937	DO k = nzt, nzb+1, -1
[1065]	2938	!
[2232]	2939	!-- Predetermine flag to mask topography
	2940	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2941	!
[2312]	2942	!-- Sum up all rain drop number densities which contribute to the flux
[1065]	2943	!-- through k-1/2
[1353]	2944	flux = 0.0_wp
	2945	z_run = 0.0_wp ! height above z(k)
[1065]	2946	k_run = k
[1346]	2947	c_run = MIN( 1.0_wp, c_nr(k) )
[1353]	2948	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
[1361]	2949	flux = flux + hyrho(k_run) * &
	2950	( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) * &
[2232]	2951	0.5_wp ) * c_run * dzu(k_run) * flag
	2952	z_run = z_run + dzu(k_run) * flag
	2953	k_run = k_run + 1 * flag
	2954	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) * flag
[1022]	2955	ENDDO
	2956	!
[2155]	2957	!-- It is not allowed to sediment more rain drop number density than
[1065]	2958	!-- available
[1361]	2959	flux = MIN( flux, &
[1115]	2960	hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
[1065]	2961
[2232]	2962	sed_nr(k) = flux / dt_micro * flag
[1361]	2963	nr_1d(k) = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / &
[2232]	2964	hyrho(k) * dt_micro * flag
[1065]	2965	!
[2155]	2966	!-- Sum up all rain water content which contributes to the flux
[1065]	2967	!-- through k-1/2
[1353]	2968	flux = 0.0_wp
	2969	z_run = 0.0_wp ! height above z(k)
[1065]	2970	k_run = k
[1346]	2971	c_run = MIN( 1.0_wp, c_qr(k) )
[1106]	2972
[1361]	2973	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
[1106]	2974
[1361]	2975	flux = flux + hyrho(k_run) * &
	2976	( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) * &
[2232]	2977	0.5_wp ) * c_run * dzu(k_run) * flag
	2978	z_run = z_run + dzu(k_run) * flag
	2979	k_run = k_run + 1 * flag
	2980	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) * flag
[1106]	2981
[1065]	2982	ENDDO
	2983	!
	2984	!-- It is not allowed to sediment more rain water content than available
[1361]	2985	flux = MIN( flux, &
[1115]	2986	hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
[1065]	2987
[2232]	2988	sed_qr(k) = flux / dt_micro * flag
[1115]	2989
[1361]	2990	qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
[2232]	2991	hyrho(k) * dt_micro * flag
[1361]	2992	q_1d(k) = q_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
[2232]	2993	hyrho(k) * dt_micro * flag
[1361]	2994	pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
[2232]	2995	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro * flag
[1065]	2996	!
	2997	!-- Compute the rain rate
[1361]	2998	IF ( call_microphysics_at_all_substeps ) THEN
[1691]	2999	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) &
[2232]	3000	* weight_substep(intermediate_timestep_count) * flag
[1361]	3001	ELSE
[2232]	3002	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
[1361]	3003	ENDIF
	3004
[1065]	3005	ENDDO
[1115]	3006
[1691]	3007	END SUBROUTINE sedimentation_rain_ij
[1012]	3008
[1691]	3009
	3010	!------------------------------------------------------------------------------!
	3011	! Description:
	3012	! ------------
	3013	!> This subroutine computes the precipitation amount due to gravitational
	3014	!> settling of rain and cloud (fog) droplets
	3015	!------------------------------------------------------------------------------!
	3016	SUBROUTINE calc_precipitation_amount_ij( i, j )
	3017
[1849]	3018	USE arrays_3d, &
	3019	ONLY: precipitation_amount, prr
	3020
[1691]	3021	USE cloud_parameters, &
[1849]	3022	ONLY: hyrho
[1691]	3023
	3024	USE control_parameters, &
	3025	ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d, &
	3026	intermediate_timestep_count, intermediate_timestep_count_max,&
[1822]	3027	precipitation_amount_interval, time_do2d_xy
[1691]	3028
	3029	USE indices, &
[2232]	3030	ONLY: nzb, nzt, wall_flags_0
[1691]	3031
	3032	USE kinds
	3033
[2232]	3034	USE surface_mod, &
	3035	ONLY : bc_h
	3036
[1691]	3037	IMPLICIT NONE
	3038
[2232]	3039	INTEGER(iwp) :: i !< running index x direction
	3040	INTEGER(iwp) :: j !< running index y direction
	3041	INTEGER(iwp) :: k !< running index z direction
	3042	INTEGER(iwp) :: m !< running index surface elements
	3043	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	3044	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1691]	3045
	3046	IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
	3047	( .NOT. call_microphysics_at_all_substeps .OR. &
	3048	intermediate_timestep_count == intermediate_timestep_count_max ) ) &
	3049	THEN
	3050
[2312]	3051	surf_s = bc_h(0)%start_index(j,i)
[2232]	3052	surf_e = bc_h(0)%end_index(j,i)
	3053	DO m = surf_s, surf_e
	3054	k = bc_h(0)%k(m)
	3055	precipitation_amount(j,i) = precipitation_amount(j,i) + &
	3056	prr(k,j,i) * hyrho(k) * dt_3d
	3057	ENDDO
	3058
[1048]	3059	ENDIF
	3060
[1691]	3061	END SUBROUTINE calc_precipitation_amount_ij
[1012]	3062
[1361]	3063	!------------------------------------------------------------------------------!
[1682]	3064	! Description:
	3065	! ------------
	3066	!> This function computes the gamma function (Press et al., 1992).
[2155]	3067	!> The gamma function is needed for the calculation of the evaporation
[1682]	3068	!> of rain drops.
[1361]	3069	!------------------------------------------------------------------------------!
[2155]	3070	FUNCTION gamm( xx )
[1320]	3071
	3072	USE kinds
	3073
[2155]	3074	IMPLICIT NONE
[1106]	3075
[1682]	3076	INTEGER(iwp) :: j !<
[1320]	3077
[1682]	3078	REAL(wp) :: gamm !<
	3079	REAL(wp) :: ser !<
	3080	REAL(wp) :: tmp !<
	3081	REAL(wp) :: x_gamm !<
	3082	REAL(wp) :: xx !<
	3083	REAL(wp) :: y_gamm !<
[1320]	3084
[1849]	3085
	3086	REAL(wp), PARAMETER :: stp = 2.5066282746310005_wp !<
	3087	REAL(wp), PARAMETER :: cof(6) = (/ 76.18009172947146_wp, &
	3088	-86.50532032941677_wp, &
	3089	24.01409824083091_wp, &
	3090	-1.231739572450155_wp, &
	3091	0.1208650973866179E-2_wp, &
	3092	-0.5395239384953E-5_wp /) !<
	3093
	3094	x_gamm = xx
	3095	y_gamm = x_gamm
[1353]	3096	tmp = x_gamm + 5.5_wp
[1849]	3097	tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
[1334]	3098	ser = 1.000000000190015_wp
[1106]	3099
[2155]	3100	DO j = 1, 6
	3101	y_gamm = y_gamm + 1.0_wp
	3102	ser = ser + cof( j ) / y_gamm
[1106]	3103	ENDDO
	3104
[2155]	3105	!
	3106	!-- Until this point the algorithm computes the logarithm of the gamma
	3107	!-- function. Hence, the exponential function is used.
	3108	! gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
	3109	gamm = EXP( tmp ) * stp * ser / x_gamm
[1106]	3110
[2155]	3111	RETURN
[1012]	3112
[2155]	3113	END FUNCTION gamm
	3114
[1012]	3115	END MODULE microphysics_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |