Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |