Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |