Home

Context Navigation

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

Last change on this file since 2312 was 2312, checked in by hoffmann, 8 years ago
various improvements of the LCM
Property svn:keywords set to `Id`
File size: 123.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |