Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |