Home

Context Navigation

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

Last change on this file since 2675 was 2608, checked in by schwenkel, 7 years ago
Inital revision of diagnostic_quantities_mod allows unified calculation of magnus equation and saturion mixing ratio
Property svn:keywords set to `Id`
File size: 122.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |