Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |