Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |