Home

Context Navigation

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

Last change on this file since 2271 was 2233, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 95.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 2233 2017-05-30 18:08:54Z sward $
	27	!
[2233]	28	! 2232 2017-05-30 17:47:52Z suehring
	29	! Adjustments to new topography and surface concept
	30	!
[2156]	31	! 2155 2017-02-21 09:57:40Z hoffmann
	32	! Bugfix in the calculation of microphysical quantities on ghost points.
	33	!
[2032]	34	! 2031 2016-10-21 15:11:58Z knoop
	35	! renamed variable rho to rho_ocean
[2155]	36	!
[2001]	37	! 2000 2016-08-20 18:09:15Z knoop
	38	! Forced header and separation lines into 80 columns
[2155]	39	!
[1851]	40	! 1850 2016-04-08 13:29:27Z maronga
	41	! Module renamed
	42	! Adapted for modularization of microphysics.
	43	!
[1846]	44	! 1845 2016-04-08 08:29:13Z raasch
	45	! nzb_2d replaced by nzb_s_inner, Kessler precipitation is stored at surface
	46	! point (instead of one point above surface)
	47	!
[1832]	48	! 1831 2016-04-07 13:15:51Z hoffmann
	49	! turbulence renamed collision_turbulence,
	50	! drizzle renamed cloud_water_sedimentation. cloud_water_sedimentation also
	51	! avaialble for microphysics_kessler.
	52	!
[1823]	53	! 1822 2016-04-07 07:49:42Z hoffmann
	54	! Unused variables removed.
	55	! Kessler scheme integrated.
	56	!
[1692]	57	! 1691 2015-10-26 16:17:44Z maronga
	58	! Added new routine calc_precipitation_amount. The routine now allows to account
	59	! for precipitation due to sedimenation of cloud (fog) droplets
[2155]	60	!
[1683]	61	! 1682 2015-10-07 23:56:08Z knoop
[2155]	62	! Code annotations made doxygen readable
[1683]	63	!
[1647]	64	! 1646 2015-09-02 16:00:10Z hoffmann
	65	! Bugfix: Wrong computation of d_mean.
	66	!
[1362]	67	! 1361 2014-04-16 15:17:48Z hoffmann
	68	! Bugfix in sedimentation_rain: Index corrected.
	69	! Vectorized version of adjust_cloud added.
	70	! Little reformatting of the code.
[2155]	71	!
[1354]	72	! 1353 2014-04-08 15:21:23Z heinze
	73	! REAL constants provided with KIND-attribute
[2155]	74	!
[1347]	75	! 1346 2014-03-27 13:18:20Z heinze
[2155]	76	! Bugfix: REAL constants provided with KIND-attribute especially in call of
[1347]	77	! intrinsic function like MAX, MIN, SIGN
[2155]	78	!
[1335]	79	! 1334 2014-03-25 12:21:40Z heinze
	80	! Bugfix: REAL constants provided with KIND-attribute
	81	!
[1323]	82	! 1322 2014-03-20 16:38:49Z raasch
	83	! REAL constants defined as wp-kind
	84	!
[1321]	85	! 1320 2014-03-20 08:40:49Z raasch
[1320]	86	! ONLY-attribute added to USE-statements,
	87	! kind-parameters added to all INTEGER and REAL declaration statements,
	88	! kinds are defined in new module kinds,
	89	! comment fields (!:) to be used for variable explanations added to
	90	! all variable declaration statements
[1000]	91	!
[1242]	92	! 1241 2013-10-30 11:36:58Z heinze
[2031]	93	! hyp and rho_ocean have to be calculated at each time step if data from external
[1242]	94	! file LSF_DATA are used
	95	!
[1116]	96	! 1115 2013-03-26 18:16:16Z hoffmann
	97	! microphyical tendencies are calculated in microphysics_control in an optimized
	98	! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
	99	!
[1107]	100	! 1106 2013-03-04 05:31:38Z raasch
	101	! small changes in code formatting
	102	!
[1093]	103	! 1092 2013-02-02 11:24:22Z raasch
	104	! unused variables removed
	105	! file put under GPL
	106	!
[1066]	107	! 1065 2012-11-22 17:42:36Z hoffmann
	108	! Sedimentation process implemented according to Stevens and Seifert (2008).
[1115]	109	! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
[1066]	110	! and Stevens, 2010).
	111	!
[1054]	112	! 1053 2012-11-13 17:11:03Z hoffmann
	113	! initial revision
[2155]	114	!
[1000]	115	! Description:
	116	! ------------
[1849]	117	!> Calculate bilk cloud microphysics.
[1000]	118	!------------------------------------------------------------------------------!
[1682]	119	MODULE microphysics_mod
[1000]	120
[1849]	121	USE kinds
	122
	123	IMPLICIT NONE
	124
	125	LOGICAL :: cloud_water_sedimentation = .FALSE. !< cloud water sedimentation
	126	LOGICAL :: limiter_sedimentation = .TRUE. !< sedimentation limiter
	127	LOGICAL :: collision_turbulence = .FALSE. !< turbulence effects
	128	LOGICAL :: ventilation_effect = .TRUE. !< ventilation effect
	129
	130	REAL(wp) :: a_1 = 8.69E-4_wp !< coef. in turb. parametrization (cm-2 s3)
	131	REAL(wp) :: a_2 = -7.38E-5_wp !< coef. in turb. parametrization (cm-2 s3)
	132	REAL(wp) :: a_3 = -1.40E-2_wp !< coef. in turb. parametrization
	133	REAL(wp) :: a_term = 9.65_wp !< coef. for terminal velocity (m s-1)
	134	REAL(wp) :: a_vent = 0.78_wp !< coef. for ventilation effect
	135	REAL(wp) :: b_1 = 11.45E-6_wp !< coef. in turb. parametrization (m)
	136	REAL(wp) :: b_2 = 9.68E-6_wp !< coef. in turb. parametrization (m)
	137	REAL(wp) :: b_3 = 0.62_wp !< coef. in turb. parametrization
	138	REAL(wp) :: b_term = 9.8_wp !< coef. for terminal velocity (m s-1)
	139	REAL(wp) :: b_vent = 0.308_wp !< coef. for ventilation effect
	140	REAL(wp) :: beta_cc = 3.09E-4_wp !< coef. in turb. parametrization (cm-2 s3)
	141	REAL(wp) :: c_1 = 4.82E-6_wp !< coef. in turb. parametrization (m)
	142	REAL(wp) :: c_2 = 4.8E-6_wp !< coef. in turb. parametrization (m)
	143	REAL(wp) :: c_3 = 0.76_wp !< coef. in turb. parametrization
	144	REAL(wp) :: c_const = 0.93_wp !< const. in Taylor-microscale Reynolds number
	145	REAL(wp) :: c_evap = 0.7_wp !< constant in evaporation
	146	REAL(wp) :: c_term = 600.0_wp !< coef. for terminal velocity (m-1)
	147	REAL(wp) :: diff_coeff_l = 0.23E-4_wp !< diffusivity of water vapor (m2 s-1)
[2155]	148	REAL(wp) :: eps_sb = 1.0E-10_wp !< threshold in two-moments scheme
[1849]	149	REAL(wp) :: k_cc = 9.44E09_wp !< const. cloud-cloud kernel (m3 kg-2 s-1)
	150	REAL(wp) :: k_cr0 = 4.33_wp !< const. cloud-rain kernel (m3 kg-1 s-1)
	151	REAL(wp) :: k_rr = 7.12_wp !< const. rain-rain kernel (m3 kg-1 s-1)
	152	REAL(wp) :: k_br = 1000.0_wp !< const. in breakup parametrization (m-1)
	153	REAL(wp) :: k_st = 1.2E8_wp !< const. in drizzle parametrization (m-1 s-1)
	154	REAL(wp) :: kappa_rr = 60.7_wp !< const. in collision kernel (kg-1/3)
	155	REAL(wp) :: kin_vis_air = 1.4086E-5_wp !< kin. viscosity of air (m2 s-1)
	156	REAL(wp) :: prec_time_const = 0.001_wp !< coef. in Kessler scheme (s-1)
	157	REAL(wp) :: ql_crit = 0.0005_wp !< coef. in Kessler scheme (kg kg-1)
	158	REAL(wp) :: schmidt_p_1d3=0.8921121_wp !< Schmidt number0.33333, 0.710.33333
	159	REAL(wp) :: sigma_gc = 1.3_wp !< geometric standard deviation cloud droplets
	160	REAL(wp) :: thermal_conductivity_l = 2.43E-2_wp !< therm. cond. air (J m-1 s-1 K-1)
	161	REAL(wp) :: w_precipitation = 9.65_wp !< maximum terminal velocity (m s-1)
	162	REAL(wp) :: x0 = 2.6E-10_wp !< separating drop mass (kg)
	163	REAL(wp) :: xrmin = 2.6E-10_wp !< minimum rain drop size (kg)
	164	REAL(wp) :: xrmax = 5.0E-6_wp !< maximum rain drop site (kg)
	165
	166	REAL(wp) :: c_sedimentation = 2.0_wp !< Courant number of sedimentation process
	167	REAL(wp) :: dpirho_l !< 6.0 / ( pi * rho_l )
	168	REAL(wp) :: dt_micro !< microphysics time step
	169	REAL(wp) :: nc_const = 70.0E6_wp !< cloud droplet concentration
	170	REAL(wp) :: dt_precipitation = 100.0_wp !< timestep precipitation (s)
	171	REAL(wp) :: sed_qc_const !< const. for sedimentation of cloud water
	172	REAL(wp) :: pirho_l !< pi * rho_l / 6.0;
	173
	174	REAL(wp), DIMENSION(:), ALLOCATABLE :: nc_1d !<
	175	REAL(wp), DIMENSION(:), ALLOCATABLE :: nr_1d !<
	176	REAL(wp), DIMENSION(:), ALLOCATABLE :: pt_1d !<
	177	REAL(wp), DIMENSION(:), ALLOCATABLE :: qc_1d !<
	178	REAL(wp), DIMENSION(:), ALLOCATABLE :: qr_1d !<
	179	REAL(wp), DIMENSION(:), ALLOCATABLE :: q_1d !<
	180
	181	SAVE
	182
[1000]	183	PRIVATE
[1849]	184	PUBLIC microphysics_control, microphysics_init
[1000]	185
[1849]	186	PUBLIC cloud_water_sedimentation, collision_turbulence, c_sedimentation, &
	187	dt_precipitation, limiter_sedimentation, nc_const, sigma_gc, &
	188	ventilation_effect
	189
[1115]	190	INTERFACE microphysics_control
	191	MODULE PROCEDURE microphysics_control
	192	MODULE PROCEDURE microphysics_control_ij
	193	END INTERFACE microphysics_control
[1022]	194
[1115]	195	INTERFACE adjust_cloud
	196	MODULE PROCEDURE adjust_cloud
	197	MODULE PROCEDURE adjust_cloud_ij
	198	END INTERFACE adjust_cloud
	199
[1000]	200	INTERFACE autoconversion
	201	MODULE PROCEDURE autoconversion
	202	MODULE PROCEDURE autoconversion_ij
	203	END INTERFACE autoconversion
	204
[1822]	205	INTERFACE autoconversion_kessler
	206	MODULE PROCEDURE autoconversion_kessler
	207	MODULE PROCEDURE autoconversion_kessler_ij
	208	END INTERFACE autoconversion_kessler
	209
[1000]	210	INTERFACE accretion
	211	MODULE PROCEDURE accretion
	212	MODULE PROCEDURE accretion_ij
	213	END INTERFACE accretion
[1005]	214
	215	INTERFACE selfcollection_breakup
	216	MODULE PROCEDURE selfcollection_breakup
	217	MODULE PROCEDURE selfcollection_breakup_ij
	218	END INTERFACE selfcollection_breakup
[1012]	219
	220	INTERFACE evaporation_rain
	221	MODULE PROCEDURE evaporation_rain
	222	MODULE PROCEDURE evaporation_rain_ij
	223	END INTERFACE evaporation_rain
	224
	225	INTERFACE sedimentation_cloud
	226	MODULE PROCEDURE sedimentation_cloud
	227	MODULE PROCEDURE sedimentation_cloud_ij
	228	END INTERFACE sedimentation_cloud
[2155]	229
[1012]	230	INTERFACE sedimentation_rain
	231	MODULE PROCEDURE sedimentation_rain
	232	MODULE PROCEDURE sedimentation_rain_ij
	233	END INTERFACE sedimentation_rain
	234
[1691]	235	INTERFACE calc_precipitation_amount
	236	MODULE PROCEDURE calc_precipitation_amount
	237	MODULE PROCEDURE calc_precipitation_amount_ij
	238	END INTERFACE calc_precipitation_amount
	239
[1000]	240	CONTAINS
[1849]	241	!------------------------------------------------------------------------------!
	242	! Description:
	243	! ------------
	244	!> Initialization of bulk microphysics
	245	!------------------------------------------------------------------------------!
	246	SUBROUTINE microphysics_init
[1000]	247
[1849]	248	USE arrays_3d, &
	249	ONLY: dzu
[1000]	250
[1849]	251	USE constants, &
	252	ONLY: pi
	253
	254	USE cloud_parameters, &
	255	ONLY: rho_l
	256
	257	USE control_parameters, &
	258	ONLY: microphysics_seifert
	259
	260	USE indices, &
	261	ONLY: nzb, nzt
	262
	263	IMPLICIT NONE
	264
	265	!
	266	!-- constant for the sedimentation of cloud water (2-moment cloud physics)
	267	sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l ) &
	268	)*( 2.0_wp / 3.0_wp ) &
	269	EXP( 5.0_wp * LOG( sigma_gc )**2 )
	270
	271	!
	272	!-- Calculate timestep according to precipitation
	273	IF ( microphysics_seifert ) THEN
	274	dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) / &
	275	w_precipitation
	276	ENDIF
	277
	278	!
	279	!-- Pre-calculate frequently calculated fractions of pi and rho_l
	280	pirho_l = pi * rho_l / 6.0_wp
	281	dpirho_l = 1.0_wp / pirho_l
	282
	283	!
	284	!-- Allocate 1D microphysics arrays
	285	ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1), q_1d(nzb:nzt+1), &
	286	qc_1d(nzb:nzt+1) )
	287
	288	IF ( microphysics_seifert ) THEN
	289	ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) )
	290	ENDIF
	291
	292	!
	293	!-- Initialze nc_1d with nc_const
	294	nc_1d = nc_const
	295
	296	END SUBROUTINE microphysics_init
	297
	298
[1000]	299	!------------------------------------------------------------------------------!
[1682]	300	! Description:
	301	! ------------
[1849]	302	!> Control of microphysics for all grid points
[1000]	303	!------------------------------------------------------------------------------!
[1115]	304	SUBROUTINE microphysics_control
[1022]	305
[1361]	306	USE arrays_3d, &
[1849]	307	ONLY: hyp, pt_init, prr, zu
[1361]	308
	309	USE cloud_parameters, &
[1849]	310	ONLY: cp, hyrho, pt_d_t, r_d, t_d_pt
[1361]	311
	312	USE control_parameters, &
[1849]	313	ONLY: call_microphysics_at_all_substeps, dt_3d, g, &
	314	intermediate_timestep_count, large_scale_forcing, &
[1822]	315	lsf_surf, microphysics_kessler, microphysics_seifert, &
	316	pt_surface, rho_surface,surface_pressure
[1361]	317
	318	USE indices, &
	319	ONLY: nzb, nzt
	320
[1320]	321	USE kinds
[1115]	322
[1361]	323	USE statistics, &
	324	ONLY: weight_pres
	325
[1115]	326	IMPLICIT NONE
	327
[1682]	328	INTEGER(iwp) :: k !<
[1115]	329
[1682]	330	REAL(wp) :: t_surface !<
[1361]	331
	332	IF ( large_scale_forcing .AND. lsf_surf ) THEN
	333	!
	334	!-- Calculate:
	335	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
	336	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
	337	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
	338	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
	339	DO k = nzb, nzt+1
	340	hyp(k) = surface_pressure * 100.0_wp * &
	341	( ( t_surface - g / cp * zu(k) ) / &
	342	t_surface )**(1.0_wp / 0.286_wp)
	343	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
	344	t_d_pt(k) = 1.0_wp / pt_d_t(k)
[2155]	345	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
[1115]	346	ENDDO
[1822]	347
[1361]	348	!
	349	!-- Compute reference density
	350	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
	351	ENDIF
[1115]	352
[1361]	353	!
[2155]	354	!-- Compute length of time step
[1361]	355	IF ( call_microphysics_at_all_substeps ) THEN
	356	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
	357	ELSE
	358	dt_micro = dt_3d
	359	ENDIF
	360
	361	!
[1822]	362	!-- Reset precipitation rate
	363	IF ( intermediate_timestep_count == 1 ) prr = 0.0_wp
	364
	365	!
[1361]	366	!-- Compute cloud physics
[1822]	367	IF ( microphysics_kessler ) THEN
	368
	369	CALL autoconversion_kessler
[1831]	370	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud
[1822]	371
	372	ELSEIF ( microphysics_seifert ) THEN
	373
[1361]	374	CALL adjust_cloud
	375	CALL autoconversion
	376	CALL accretion
	377	CALL selfcollection_breakup
	378	CALL evaporation_rain
	379	CALL sedimentation_rain
[1831]	380	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud
[1361]	381
[1691]	382	ENDIF
	383
[1822]	384	CALL calc_precipitation_amount
	385
[1115]	386	END SUBROUTINE microphysics_control
	387
[1682]	388	!------------------------------------------------------------------------------!
	389	! Description:
	390	! ------------
[2155]	391	!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
	392	!> evaporation of rain drops due to too small or too big weights
[1682]	393	!> of rain drops (Stevens and Seifert, 2008).
	394	!------------------------------------------------------------------------------!
[1115]	395	SUBROUTINE adjust_cloud
	396
[1361]	397	USE arrays_3d, &
	398	ONLY: qr, nr
	399
	400	USE cloud_parameters, &
[1849]	401	ONLY: hyrho
[1361]	402
	403	USE cpulog, &
	404	ONLY: cpu_log, log_point_s
	405
	406	USE indices, &
[2232]	407	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
[1361]	408
[1320]	409	USE kinds
[1022]	410
	411	IMPLICIT NONE
	412
[1682]	413	INTEGER(iwp) :: i !<
	414	INTEGER(iwp) :: j !<
	415	INTEGER(iwp) :: k !<
[1022]	416
[1361]	417	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
	418
[2155]	419	DO i = nxlg, nxrg
	420	DO j = nysg, nyng
[2232]	421	DO k = nzb+1, nzt
[1361]	422	IF ( qr(k,j,i) <= eps_sb ) THEN
	423	qr(k,j,i) = 0.0_wp
	424	nr(k,j,i) = 0.0_wp
	425	ELSE
	426	IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN
[2232]	427	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * &
	428	MERGE( 1.0_wp, 0.0_wp, &
	429	BTEST( wall_flags_0(k,j,i), 0 ) )
[1361]	430	ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN
[2232]	431	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * &
	432	MERGE( 1.0_wp, 0.0_wp, &
	433	BTEST( wall_flags_0(k,j,i), 0 ) )
[1361]	434	ENDIF
	435	ENDIF
[1022]	436	ENDDO
	437	ENDDO
	438	ENDDO
	439
[1361]	440	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
	441
[1115]	442	END SUBROUTINE adjust_cloud
[1022]	443
[1106]	444
[1682]	445	!------------------------------------------------------------------------------!
	446	! Description:
	447	! ------------
	448	!> Autoconversion rate (Seifert and Beheng, 2006).
	449	!------------------------------------------------------------------------------!
[1000]	450	SUBROUTINE autoconversion
	451
[1361]	452	USE arrays_3d, &
	453	ONLY: diss, dzu, nr, qc, qr
	454
	455	USE cloud_parameters, &
[1849]	456	ONLY: hyrho
[1361]	457
	458	USE control_parameters, &
[1849]	459	ONLY: rho_surface
[1361]	460
	461	USE cpulog, &
	462	ONLY: cpu_log, log_point_s
	463
	464	USE grid_variables, &
	465	ONLY: dx, dy
	466
	467	USE indices, &
[2232]	468	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
[1361]	469
[1320]	470	USE kinds
[1000]	471
	472	IMPLICIT NONE
	473
[1682]	474	INTEGER(iwp) :: i !<
	475	INTEGER(iwp) :: j !<
	476	INTEGER(iwp) :: k !<
[1000]	477
[2155]	478	REAL(wp) :: alpha_cc !<
[1682]	479	REAL(wp) :: autocon !<
	480	REAL(wp) :: dissipation !<
[2232]	481	REAL(wp) :: flag !< flag to mask topography grid points
[1682]	482	REAL(wp) :: k_au !<
	483	REAL(wp) :: l_mix !<
	484	REAL(wp) :: nu_c !<
	485	REAL(wp) :: phi_au !<
	486	REAL(wp) :: r_cc !<
	487	REAL(wp) :: rc !<
	488	REAL(wp) :: re_lambda !<
	489	REAL(wp) :: sigma_cc !<
	490	REAL(wp) :: tau_cloud !<
	491	REAL(wp) :: xc !<
[1361]	492
	493	CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
	494
[2155]	495	DO i = nxlg, nxrg
	496	DO j = nysg, nyng
[2232]	497	DO k = nzb+1, nzt
	498	!
	499	!-- Predetermine flag to mask topography
	500	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1000]	501
[1361]	502	IF ( qc(k,j,i) > eps_sb ) THEN
	503
	504	k_au = k_cc / ( 20.0_wp * x0 )
	505	!
	506	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
	507	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
	508	tau_cloud = 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) )
	509	!
[2155]	510	!-- Universal function for autoconversion process
[1361]	511	!-- (Seifert and Beheng, 2006):
	512	phi_au = 600.0_wp * tau_cloud*0.68_wp &
	513	( 1.0_wp - tau_cloud0.68_wp )3
	514	!
	515	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
	516	!-- (Use constant nu_c = 1.0_wp instead?)
	517	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
	518	!
	519	!-- Mean weight of cloud droplets:
	520	xc = hyrho(k) * qc(k,j,i) / nc_const
	521	!
[2155]	522	!-- Parameterized turbulence effects on autoconversion (Seifert,
[1361]	523	!-- Nuijens and Stevens, 2010)
[1831]	524	IF ( collision_turbulence ) THEN
[1361]	525	!
	526	!-- Weight averaged radius of cloud droplets:
	527	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
	528
	529	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
	530	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
	531	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
	532	!
[2155]	533	!-- Mixing length (neglecting distance to ground and
[1361]	534	!-- stratification)
	535	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
	536	!
[2155]	537	!-- Limit dissipation rate according to Seifert, Nuijens and
[1361]	538	!-- Stevens (2010)
	539	dissipation = MIN( 0.06_wp, diss(k,j,i) )
	540	!
	541	!-- Compute Taylor-microscale Reynolds number:
	542	re_lambda = 6.0_wp / 11.0_wp * &
	543	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
	544	SQRT( 15.0_wp / kin_vis_air ) * &
	545	dissipation**( 1.0_wp / 6.0_wp )
	546	!
[2155]	547	!-- The factor of 1.0E4 is needed to convert the dissipation
[1361]	548	!-- rate from m2 s-3 to cm2 s-3.
	549	k_au = k_au * ( 1.0_wp + &
	550	dissipation * 1.0E4_wp * &
	551	( re_lambda * 1.0E-3_wp )*0.25_wp &
[2155]	552	( alpha_cc * EXP( -1.0_wp * ( ( rc - &
[1361]	553	r_cc ) / &
	554	sigma_cc )**2 &
	555	) + beta_cc &
	556	) &
	557	)
	558	ENDIF
	559	!
	560	!-- Autoconversion rate (Seifert and Beheng, 2006):
	561	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
	562	( nu_c + 1.0_wp )*2 qc(k,j,i)*2 xc*2 &
	563	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
	564	rho_surface
	565	autocon = MIN( autocon, qc(k,j,i) / dt_micro )
	566
[2232]	567	qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag
	568	qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag
	569	nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro &
	570	* flag
[1361]	571
	572	ENDIF
	573
[1000]	574	ENDDO
	575	ENDDO
	576	ENDDO
	577
[1361]	578	CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
	579
[1000]	580	END SUBROUTINE autoconversion
	581
[1106]	582
[1682]	583	!------------------------------------------------------------------------------!
	584	! Description:
	585	! ------------
[1822]	586	!> Autoconversion process (Kessler, 1969).
	587	!------------------------------------------------------------------------------!
	588	SUBROUTINE autoconversion_kessler
	589
	590	USE arrays_3d, &
[1849]	591	ONLY: dzw, pt, prr, q, qc
[1822]	592
	593	USE cloud_parameters, &
[1849]	594	ONLY: l_d_cp, pt_d_t
[1822]	595
	596	USE indices, &
[2232]	597	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
[1822]	598
	599	USE kinds
	600
	601
	602	IMPLICIT NONE
	603
[2232]	604	INTEGER(iwp) :: i !<
	605	INTEGER(iwp) :: j !<
	606	INTEGER(iwp) :: k !<
	607	INTEGER(iwp) :: k_wall !< topgraphy top index
[1822]	608
	609	REAL(wp) :: dqdt_precip !<
[2232]	610	REAL(wp) :: flag !< flag to mask topography grid points
[1822]	611
[2155]	612	DO i = nxlg, nxrg
	613	DO j = nysg, nyng
[2232]	614	!
	615	!-- Determine vertical index of topography top
	616	k_wall = MAXLOC( &
	617	MERGE( 1, 0, &
	618	BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 ) &
	619	), DIM = 1 &
	620	) - 1
	621	DO k = nzb+1, nzt
	622	!
	623	!-- Predetermine flag to mask topography
	624	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1822]	625
	626	IF ( qc(k,j,i) > ql_crit ) THEN
	627	dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
	628	ELSE
	629	dqdt_precip = 0.0_wp
	630	ENDIF
	631
[2232]	632	qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
	633	q(k,j,i) = q(k,j,i) - dqdt_precip * dt_micro * flag
[1845]	634	pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * l_d_cp * &
[2232]	635	pt_d_t(k) * flag
[1822]	636
	637	!
[1845]	638	!-- Compute the rain rate (stored on surface grid point)
[2232]	639	prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
[1822]	640
	641	ENDDO
	642	ENDDO
	643	ENDDO
	644
	645	END SUBROUTINE autoconversion_kessler
	646
	647
	648	!------------------------------------------------------------------------------!
	649	! Description:
	650	! ------------
[1682]	651	!> Accretion rate (Seifert and Beheng, 2006).
	652	!------------------------------------------------------------------------------!
[1005]	653	SUBROUTINE accretion
[1000]	654
[1361]	655	USE arrays_3d, &
	656	ONLY: diss, qc, qr
	657
	658	USE cloud_parameters, &
[1849]	659	ONLY: hyrho
[1361]	660
	661	USE control_parameters, &
[1849]	662	ONLY: rho_surface
[1361]	663
	664	USE cpulog, &
	665	ONLY: cpu_log, log_point_s
	666
	667	USE indices, &
[2232]	668	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
[1361]	669
[1320]	670	USE kinds
[1005]	671
[1000]	672	IMPLICIT NONE
	673
[1682]	674	INTEGER(iwp) :: i !<
	675	INTEGER(iwp) :: j !<
	676	INTEGER(iwp) :: k !<
[1000]	677
[1682]	678	REAL(wp) :: accr !<
[2232]	679	REAL(wp) :: flag !< flag to mask topography grid points
[1682]	680	REAL(wp) :: k_cr !<
	681	REAL(wp) :: phi_ac !<
	682	REAL(wp) :: tau_cloud !<
[1361]	683
	684	CALL cpu_log( log_point_s(56), 'accretion', 'start' )
	685
[2155]	686	DO i = nxlg, nxrg
	687	DO j = nysg, nyng
[2232]	688	DO k = nzb+1, nzt
	689	!
	690	!-- Predetermine flag to mask topography
	691	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1000]	692
[1361]	693	IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) ) THEN
	694	!
	695	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[2155]	696	tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
[1361]	697	!
[2155]	698	!-- Universal function for accretion process (Seifert and
[1361]	699	!-- Beheng, 2001):
	700	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
	701	!
[2155]	702	!-- Parameterized turbulence effects on autoconversion (Seifert,
	703	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
[1361]	704	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
[1831]	705	IF ( collision_turbulence ) THEN
[1361]	706	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
	707	MIN( 600.0_wp, &
	708	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
	709	)
	710	ELSE
[2155]	711	k_cr = k_cr0
[1361]	712	ENDIF
	713	!
	714	!-- Accretion rate (Seifert and Beheng, 2006):
	715	accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * &
	716	SQRT( rho_surface * hyrho(k) )
	717	accr = MIN( accr, qc(k,j,i) / dt_micro )
	718
[2232]	719	qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
	720	qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
[1361]	721
	722	ENDIF
	723
[1005]	724	ENDDO
	725	ENDDO
[1000]	726	ENDDO
	727
[1361]	728	CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
	729
[1005]	730	END SUBROUTINE accretion
[1000]	731
[1106]	732
[1682]	733	!------------------------------------------------------------------------------!
	734	! Description:
	735	! ------------
	736	!> Collisional breakup rate (Seifert, 2008).
	737	!------------------------------------------------------------------------------!
[1005]	738	SUBROUTINE selfcollection_breakup
[1000]	739
[1361]	740	USE arrays_3d, &
	741	ONLY: nr, qr
	742
	743	USE cloud_parameters, &
[1849]	744	ONLY: hyrho
[1361]	745
	746	USE control_parameters, &
[1849]	747	ONLY: rho_surface
[1361]	748
	749	USE cpulog, &
	750	ONLY: cpu_log, log_point_s
	751
	752	USE indices, &
[2232]	753	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
[1361]	754
[1320]	755	USE kinds
[2155]	756
[1000]	757	IMPLICIT NONE
	758
[1682]	759	INTEGER(iwp) :: i !<
	760	INTEGER(iwp) :: j !<
	761	INTEGER(iwp) :: k !<
[1000]	762
[1682]	763	REAL(wp) :: breakup !<
	764	REAL(wp) :: dr !<
[2232]	765	REAL(wp) :: flag !< flag to mask topography grid points
[1682]	766	REAL(wp) :: phi_br !<
	767	REAL(wp) :: selfcoll !<
[1361]	768
	769	CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
	770
[2155]	771	DO i = nxlg, nxrg
	772	DO j = nysg, nyng
[2232]	773	DO k = nzb+1, nzt
	774	!
	775	!-- Predetermine flag to mask topography
	776	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	777
[1361]	778	IF ( qr(k,j,i) > eps_sb ) THEN
	779	!
	780	!-- Selfcollection rate (Seifert and Beheng, 2001):
	781	selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * &
	782	SQRT( hyrho(k) * rho_surface )
	783	!
	784	!-- Weight averaged diameter of rain drops:
	785	dr = ( hyrho(k) * qr(k,j,i) / &
	786	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
	787	!
	788	!-- Collisional breakup rate (Seifert, 2008):
	789	IF ( dr >= 0.3E-3_wp ) THEN
	790	phi_br = k_br * ( dr - 1.1E-3_wp )
	791	breakup = selfcoll * ( phi_br + 1.0_wp )
	792	ELSE
	793	breakup = 0.0_wp
	794	ENDIF
[1000]	795
[1361]	796	selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
[2232]	797	nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
[1361]	798
[2155]	799	ENDIF
[1000]	800	ENDDO
	801	ENDDO
	802	ENDDO
	803
[1361]	804	CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
	805
[1005]	806	END SUBROUTINE selfcollection_breakup
[1000]	807
[1106]	808
[1682]	809	!------------------------------------------------------------------------------!
	810	! Description:
	811	! ------------
[2155]	812	!> Evaporation of precipitable water. Condensation is neglected for
[1682]	813	!> precipitable water.
	814	!------------------------------------------------------------------------------!
[1012]	815	SUBROUTINE evaporation_rain
[1000]	816
[1361]	817	USE arrays_3d, &
	818	ONLY: hyp, nr, pt, q, qc, qr
	819
	820	USE cloud_parameters, &
[1849]	821	ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
[1361]	822
	823	USE constants, &
	824	ONLY: pi
	825
	826	USE cpulog, &
	827	ONLY: cpu_log, log_point_s
	828
	829	USE indices, &
[2232]	830	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
[1361]	831
[1320]	832	USE kinds
[1012]	833
	834	IMPLICIT NONE
	835
[1682]	836	INTEGER(iwp) :: i !<
	837	INTEGER(iwp) :: j !<
	838	INTEGER(iwp) :: k !<
[1361]	839
[1682]	840	REAL(wp) :: alpha !<
	841	REAL(wp) :: dr !<
	842	REAL(wp) :: e_s !<
	843	REAL(wp) :: evap !<
	844	REAL(wp) :: evap_nr !<
	845	REAL(wp) :: f_vent !<
[2232]	846	REAL(wp) :: flag !< flag to mask topography grid points
[1682]	847	REAL(wp) :: g_evap !<
	848	REAL(wp) :: lambda_r !<
	849	REAL(wp) :: mu_r !<
	850	REAL(wp) :: mu_r_2 !<
	851	REAL(wp) :: mu_r_5d2 !<
	852	REAL(wp) :: nr_0 !<
	853	REAL(wp) :: q_s !<
	854	REAL(wp) :: sat !<
	855	REAL(wp) :: t_l !<
	856	REAL(wp) :: temp !<
	857	REAL(wp) :: xr !<
[1361]	858
	859	CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
	860
[2155]	861	DO i = nxlg, nxrg
	862	DO j = nysg, nyng
[2232]	863	DO k = nzb+1, nzt
	864	!
	865	!-- Predetermine flag to mask topography
	866	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	867
[1361]	868	IF ( qr(k,j,i) > eps_sb ) THEN
	869	!
	870	!-- Actual liquid water temperature:
	871	t_l = t_d_pt(k) * pt(k,j,i)
	872	!
	873	!-- Saturation vapor pressure at t_l:
	874	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
	875	( t_l - 35.86_wp ) &
	876	)
	877	!
	878	!-- Computation of saturation humidity:
	879	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
	880	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
	881	q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / &
	882	( 1.0_wp + alpha * q_s )
	883	!
	884	!-- Supersaturation:
	885	sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
	886	!
	887	!-- Evaporation needs only to be calculated in subsaturated regions
	888	IF ( sat < 0.0_wp ) THEN
	889	!
	890	!-- Actual temperature:
	891	temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
[2155]	892
[1361]	893	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
	894	l_v / ( thermal_conductivity_l * temp ) &
	895	+ r_v * temp / ( diff_coeff_l * e_s ) &
	896	)
	897	!
	898	!-- Mean weight of rain drops
	899	xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
	900	!
	901	!-- Weight averaged diameter of rain drops:
	902	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
	903	!
[2155]	904	!-- Compute ventilation factor and intercept parameter
[1361]	905	!-- (Seifert and Beheng, 2006; Seifert, 2008):
	906	IF ( ventilation_effect ) THEN
	907	!
[2155]	908	!-- Shape parameter of gamma distribution (Milbrandt and Yau,
[1361]	909	!-- 2005; Stevens and Seifert, 2008):
	910	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
	911	( dr - 1.4E-3_wp ) ) )
	912	!
	913	!-- Slope parameter of gamma distribution (Seifert, 2008):
	914	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	915	( mu_r + 1.0_wp ) &
	916	)**( 1.0_wp / 3.0_wp ) / dr
[1012]	917
[1361]	918	mu_r_2 = mu_r + 2.0_wp
[2155]	919	mu_r_5d2 = mu_r + 2.5_wp
[1361]	920
	921	f_vent = a_vent * gamm( mu_r_2 ) * &
	922	lambda_r*( -mu_r_2 ) + b_vent &
	923	schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
	924	gamm( mu_r_5d2 ) * lambda_r*( -mu_r_5d2 ) &
	925	( 1.0_wp - &
	926	0.5_wp * ( b_term / a_term ) * &
	927	( lambda_r / ( c_term + lambda_r ) &
	928	)**mu_r_5d2 - &
	929	0.125_wp * ( b_term / a_term )*2 &
	930	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
	931	)**mu_r_5d2 - &
	932	0.0625_wp * ( b_term / a_term )*3 &
	933	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
	934	)**mu_r_5d2 - &
[2155]	935	0.0390625_wp * ( b_term / a_term )*4 &
[1361]	936	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
	937	)**mu_r_5d2 &
	938	)
	939
	940	nr_0 = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) / &
[2155]	941	gamm( mu_r + 1.0_wp )
[1361]	942	ELSE
	943	f_vent = 1.0_wp
	944	nr_0 = nr(k,j,i) * dr
	945	ENDIF
	946	!
[2155]	947	!-- Evaporation rate of rain water content (Seifert and
[1361]	948	!-- Beheng, 2006):
	949	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / &
	950	hyrho(k)
	951	evap = MAX( evap, -qr(k,j,i) / dt_micro )
	952	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
	953	-nr(k,j,i) / dt_micro )
	954
[2232]	955	qr(k,j,i) = qr(k,j,i) + evap * dt_micro * flag
	956	nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag
[1361]	957
	958	ENDIF
[2155]	959	ENDIF
[1361]	960
[1012]	961	ENDDO
	962	ENDDO
	963	ENDDO
	964
[1361]	965	CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
	966
[1012]	967	END SUBROUTINE evaporation_rain
	968
[1106]	969
[1682]	970	!------------------------------------------------------------------------------!
	971	! Description:
	972	! ------------
	973	!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
	974	!------------------------------------------------------------------------------!
[1012]	975	SUBROUTINE sedimentation_cloud
	976
[1361]	977	USE arrays_3d, &
[1849]	978	ONLY: ddzu, dzu, pt, prr, q, qc
[1361]	979
	980	USE cloud_parameters, &
[1849]	981	ONLY: hyrho, l_d_cp, pt_d_t
[1361]	982
	983	USE control_parameters, &
[1849]	984	ONLY: call_microphysics_at_all_substeps, intermediate_timestep_count
[1361]	985
	986	USE cpulog, &
	987	ONLY: cpu_log, log_point_s
	988
	989	USE indices, &
[2232]	990	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
[1361]	991
[1320]	992	USE kinds
[1691]	993
	994	USE statistics, &
	995	ONLY: weight_substep
	996
[2155]	997
[1012]	998	IMPLICIT NONE
	999
[1849]	1000	INTEGER(iwp) :: i !<
	1001	INTEGER(iwp) :: j !<
	1002	INTEGER(iwp) :: k !<
[1361]	1003
[2232]	1004	REAL(wp) :: flag !< flag to mask topography grid points
	1005	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
[1361]	1006
	1007	CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
	1008
	1009	sed_qc(nzt+1) = 0.0_wp
	1010
[2155]	1011	DO i = nxlg, nxrg
	1012	DO j = nysg, nyng
[2232]	1013	DO k = nzt, nzb+1, -1
	1014	!
	1015	!-- Predetermine flag to mask topography
	1016	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1012]	1017
[1361]	1018	IF ( qc(k,j,i) > eps_sb ) THEN
	1019	sed_qc(k) = sed_qc_const * nc_const*( -2.0_wp / 3.0_wp ) &
[2232]	1020	( qc(k,j,i) * hyrho(k) )*( 5.0_wp / 3.0_wp ) &
	1021	flag
[1361]	1022	ELSE
	1023	sed_qc(k) = 0.0_wp
	1024	ENDIF
	1025
	1026	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / &
	1027	dt_micro + sed_qc(k+1) &
[2232]	1028	) * flag
[1361]	1029
	1030	q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
[2232]	1031	ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]	1032	qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
[2232]	1033	ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]	1034	pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * &
	1035	ddzu(k+1) / hyrho(k) * l_d_cp * &
[2232]	1036	pt_d_t(k) * dt_micro * flag
[1361]	1037
[1691]	1038	!
	1039	!-- Compute the precipitation rate due to cloud (fog) droplets
[1822]	1040	IF ( call_microphysics_at_all_substeps ) THEN
	1041	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) &
[2232]	1042	* weight_substep(intermediate_timestep_count) &
	1043	* flag
[1822]	1044	ELSE
[2232]	1045	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
[1691]	1046	ENDIF
	1047
[1012]	1048	ENDDO
	1049	ENDDO
	1050	ENDDO
	1051
[1361]	1052	CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
	1053
[1012]	1054	END SUBROUTINE sedimentation_cloud
	1055
[1106]	1056
[1682]	1057	!------------------------------------------------------------------------------!
	1058	! Description:
	1059	! ------------
	1060	!> Computation of sedimentation flux. Implementation according to Stevens
	1061	!> and Seifert (2008). Code is based on UCLA-LES.
	1062	!------------------------------------------------------------------------------!
[1012]	1063	SUBROUTINE sedimentation_rain
	1064
[1361]	1065	USE arrays_3d, &
[1849]	1066	ONLY: ddzu, dzu, nr, pt, prr, q, qr
[1361]	1067
	1068	USE cloud_parameters, &
[1849]	1069	ONLY: hyrho, l_d_cp, pt_d_t
[1361]	1070
	1071	USE control_parameters, &
[1849]	1072	ONLY: call_microphysics_at_all_substeps, intermediate_timestep_count
[1361]	1073	USE cpulog, &
	1074	ONLY: cpu_log, log_point_s
	1075
	1076	USE indices, &
[2232]	1077	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
[1361]	1078
[1320]	1079	USE kinds
[1012]	1080
[1361]	1081	USE statistics, &
	1082	ONLY: weight_substep
[2155]	1083
[2232]	1084	USE surface_mod, &
	1085	ONLY : bc_h
	1086
[1012]	1087	IMPLICIT NONE
	1088
[2232]	1089	INTEGER(iwp) :: i !< running index x direction
	1090	INTEGER(iwp) :: j !< running index y direction
	1091	INTEGER(iwp) :: k !< running index z direction
	1092	INTEGER(iwp) :: k_run !<
	1093	INTEGER(iwp) :: l !< running index of surface type
	1094	INTEGER(iwp) :: m !< running index surface elements
	1095	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	1096	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1361]	1097
[1682]	1098	REAL(wp) :: c_run !<
	1099	REAL(wp) :: d_max !<
	1100	REAL(wp) :: d_mean !<
	1101	REAL(wp) :: d_min !<
	1102	REAL(wp) :: dr !<
	1103	REAL(wp) :: flux !<
[2232]	1104	REAL(wp) :: flag !< flag to mask topography grid points
[1682]	1105	REAL(wp) :: lambda_r !<
	1106	REAL(wp) :: mu_r !<
	1107	REAL(wp) :: z_run !<
[1361]	1108
[1682]	1109	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<
	1110	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<
	1111	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<
	1112	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<
	1113	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<
	1114	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<
	1115	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<
	1116	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<
[1361]	1117
	1118	CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
[1682]	1119
[1361]	1120	!
[2155]	1121	!-- Compute velocities
	1122	DO i = nxlg, nxrg
	1123	DO j = nysg, nyng
[2232]	1124	DO k = nzb+1, nzt
	1125	!
	1126	!-- Predetermine flag to mask topography
	1127	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1128
[1361]	1129	IF ( qr(k,j,i) > eps_sb ) THEN
	1130	!
	1131	!-- Weight averaged diameter of rain drops:
	1132	dr = ( hyrho(k) * qr(k,j,i) / &
	1133	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
	1134	!
	1135	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	1136	!-- Stevens and Seifert, 2008):
	1137	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
	1138	( dr - 1.4E-3_wp ) ) )
	1139	!
	1140	!-- Slope parameter of gamma distribution (Seifert, 2008):
	1141	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	1142	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
[1012]	1143
[1361]	1144	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	1145	a_term - b_term * ( 1.0_wp + &
	1146	c_term / &
	1147	lambda_r )*( -1.0_wp &
	1148	( mu_r + 1.0_wp ) ) &
	1149	) &
[2232]	1150	) * flag
[1361]	1151
	1152	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	1153	a_term - b_term * ( 1.0_wp + &
	1154	c_term / &
	1155	lambda_r )*( -1.0_wp &
	1156	( mu_r + 4.0_wp ) ) &
	1157	) &
[2232]	1158	) * flag
[1361]	1159	ELSE
	1160	w_nr(k) = 0.0_wp
	1161	w_qr(k) = 0.0_wp
	1162	ENDIF
[1012]	1163	ENDDO
[1361]	1164	!
[2232]	1165	!-- Adjust boundary values using surface data type.
	1166	!-- Upward-facing
	1167	surf_s = bc_h(0)%start_index(j,i)
	1168	surf_e = bc_h(0)%end_index(j,i)
	1169	DO m = surf_s, surf_e
	1170	k = bc_h(0)%k(m)
	1171	w_nr(k-1) = w_nr(k)
	1172	w_qr(k-1) = w_qr(k)
	1173	ENDDO
	1174	!
	1175	!-- Downward-facing
	1176	surf_s = bc_h(1)%start_index(j,i)
	1177	surf_e = bc_h(1)%end_index(j,i)
	1178	DO m = surf_s, surf_e
	1179	k = bc_h(1)%k(m)
	1180	w_nr(k+1) = w_nr(k)
	1181	w_qr(k+1) = w_qr(k)
	1182	ENDDO
	1183	!
	1184	!-- Model top boundary value
[1361]	1185	w_nr(nzt+1) = 0.0_wp
	1186	w_qr(nzt+1) = 0.0_wp
	1187	!
	1188	!-- Compute Courant number
[2232]	1189	DO k = nzb+1, nzt
	1190	!
	1191	!-- Predetermine flag to mask topography
	1192	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1193
[1361]	1194	c_nr(k) = 0.25_wp * ( w_nr(k-1) + &
	1195	2.0_wp * w_nr(k) + w_nr(k+1) ) * &
[2232]	1196	dt_micro * ddzu(k) * flag
[1361]	1197	c_qr(k) = 0.25_wp * ( w_qr(k-1) + &
	1198	2.0_wp * w_qr(k) + w_qr(k+1) ) * &
[2232]	1199	dt_micro * ddzu(k) * flag
[2155]	1200	ENDDO
[1361]	1201	!
	1202	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
	1203	IF ( limiter_sedimentation ) THEN
	1204
[2232]	1205	DO k = nzb+1, nzt
	1206	!
	1207	!-- Predetermine flag to mask topography
	1208	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1209
[1646]	1210	d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
[1361]	1211	d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
	1212	d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
	1213
	1214	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	1215	2.0_wp * d_max, &
[2232]	1216	ABS( d_mean ) ) &
	1217	* flag
[1361]	1218
[1646]	1219	d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
[1361]	1220	d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
	1221	d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
	1222
	1223	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	1224	2.0_wp * d_max, &
	1225	ABS( d_mean ) )
	1226	ENDDO
	1227
	1228	ELSE
	1229
	1230	nr_slope = 0.0_wp
	1231	qr_slope = 0.0_wp
	1232
	1233	ENDIF
	1234
	1235	sed_nr(nzt+1) = 0.0_wp
	1236	sed_qr(nzt+1) = 0.0_wp
	1237	!
	1238	!-- Compute sedimentation flux
[2232]	1239	DO k = nzt, nzb+1, -1
[1361]	1240	!
[2232]	1241	!-- Predetermine flag to mask topography
	1242	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1243	!
	1244	!-- Sum up all rain drop number densities which contribute to the flux
[1361]	1245	!-- through k-1/2
	1246	flux = 0.0_wp
	1247	z_run = 0.0_wp ! height above z(k)
	1248	k_run = k
	1249	c_run = MIN( 1.0_wp, c_nr(k) )
	1250	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
	1251	flux = flux + hyrho(k_run) * &
	1252	( nr(k_run,j,i) + nr_slope(k_run) * &
[2232]	1253	( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run) &
	1254	* flag
	1255	z_run = z_run + dzu(k_run) * flag
	1256	k_run = k_run + 1 * flag
	1257	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) &
	1258	* flag
[1361]	1259	ENDDO
	1260	!
[2155]	1261	!-- It is not allowed to sediment more rain drop number density than
[1361]	1262	!-- available
	1263	flux = MIN( flux, &
	1264	hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * &
	1265	dt_micro &
	1266	)
	1267
[2232]	1268	sed_nr(k) = flux / dt_micro * flag
[1361]	1269	nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * &
[2232]	1270	ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]	1271	!
[2155]	1272	!-- Sum up all rain water content which contributes to the flux
[1361]	1273	!-- through k-1/2
	1274	flux = 0.0_wp
	1275	z_run = 0.0_wp ! height above z(k)
	1276	k_run = k
	1277	c_run = MIN( 1.0_wp, c_qr(k) )
	1278
	1279	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
	1280
	1281	flux = flux + hyrho(k_run) * ( qr(k_run,j,i) + &
	1282	qr_slope(k_run) * ( 1.0_wp - c_run ) * &
[2232]	1283	0.5_wp ) * c_run * dzu(k_run) * flag
	1284	z_run = z_run + dzu(k_run) * flag
	1285	k_run = k_run + 1 * flag
	1286	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) &
	1287	* flag
[1361]	1288
	1289	ENDDO
	1290	!
[2155]	1291	!-- It is not allowed to sediment more rain water content than
[1361]	1292	!-- available
	1293	flux = MIN( flux, &
	1294	hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * &
	1295	dt_micro &
	1296	)
	1297
[2232]	1298	sed_qr(k) = flux / dt_micro * flag
[1361]	1299
	1300	qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
[2232]	1301	ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]	1302	q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
[2232]	1303	ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]	1304	pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * &
	1305	ddzu(k+1) / hyrho(k) * l_d_cp * &
[2232]	1306	pt_d_t(k) * dt_micro * flag
[1361]	1307	!
	1308	!-- Compute the rain rate
	1309	IF ( call_microphysics_at_all_substeps ) THEN
[1691]	1310	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) &
[2232]	1311	* weight_substep(intermediate_timestep_count) &
	1312	* flag
[1361]	1313	ELSE
[2232]	1314	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
[1361]	1315	ENDIF
	1316
	1317	ENDDO
[1012]	1318	ENDDO
	1319	ENDDO
	1320
[1691]	1321	CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
	1322
	1323	END SUBROUTINE sedimentation_rain
	1324
	1325
	1326	!------------------------------------------------------------------------------!
	1327	! Description:
	1328	! ------------
	1329	!> Computation of the precipitation amount due to gravitational settling of
	1330	!> rain and cloud (fog) droplets
	1331	!------------------------------------------------------------------------------!
	1332	SUBROUTINE calc_precipitation_amount
	1333
[1849]	1334	USE arrays_3d, &
	1335	ONLY: precipitation_amount, prr
	1336
[1691]	1337	USE cloud_parameters, &
[1849]	1338	ONLY: hyrho
[1691]	1339
	1340	USE control_parameters, &
	1341	ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d, &
	1342	intermediate_timestep_count, intermediate_timestep_count_max,&
	1343	precipitation_amount_interval, time_do2d_xy
	1344
	1345	USE indices, &
[2232]	1346	ONLY: nxl, nxr, nys, nyn, nzb, nzt, wall_flags_0
[1691]	1347
	1348	USE kinds
	1349
[2232]	1350	USE surface_mod, &
	1351	ONLY : bc_h
	1352
[1691]	1353	IMPLICIT NONE
	1354
[2232]	1355	INTEGER(iwp) :: i !< running index x direction
	1356	INTEGER(iwp) :: j !< running index y direction
	1357	INTEGER(iwp) :: k !< running index y direction
	1358	INTEGER(iwp) :: m !< running index surface elements
	1359	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	1360	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1691]	1361
	1362	IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
	1363	( .NOT. call_microphysics_at_all_substeps .OR. &
	1364	intermediate_timestep_count == intermediate_timestep_count_max ) ) &
	1365	THEN
[2232]	1366	!
	1367	!-- Run over all upward-facing surface elements, i.e. non-natural,
	1368	!-- natural and urban
	1369	DO m = 1, bc_h(0)%ns
	1370	i = bc_h(0)%i(m)
	1371	j = bc_h(0)%j(m)
	1372	k = bc_h(0)%k(m)
	1373	precipitation_amount(j,i) = precipitation_amount(j,i) + &
	1374	prr(k,j,i) * hyrho(k) * dt_3d
	1375	ENDDO
[1691]	1376
[1361]	1377	ENDIF
	1378
[1691]	1379	END SUBROUTINE calc_precipitation_amount
[1361]	1380
[1012]	1381
[1000]	1382	!------------------------------------------------------------------------------!
[1682]	1383	! Description:
	1384	! ------------
[1849]	1385	!> Control of microphysics for grid points i,j
[1000]	1386	!------------------------------------------------------------------------------!
[1022]	1387
[1115]	1388	SUBROUTINE microphysics_control_ij( i, j )
	1389
[1320]	1390	USE arrays_3d, &
[1849]	1391	ONLY: hyp, nr, pt, pt_init, prr, q, qc, qr, zu
[1115]	1392
[1320]	1393	USE cloud_parameters, &
[1849]	1394	ONLY: cp, hyrho, pt_d_t, r_d, t_d_pt
[1320]	1395
	1396	USE control_parameters, &
[1849]	1397	ONLY: call_microphysics_at_all_substeps, dt_3d, g, &
	1398	intermediate_timestep_count, large_scale_forcing, &
[1822]	1399	lsf_surf, microphysics_seifert, microphysics_kessler, &
	1400	pt_surface, rho_surface, surface_pressure
[1320]	1401
	1402	USE indices, &
	1403	ONLY: nzb, nzt
	1404
	1405	USE kinds
	1406
	1407	USE statistics, &
	1408	ONLY: weight_pres
	1409
[1022]	1410	IMPLICIT NONE
	1411
[1682]	1412	INTEGER(iwp) :: i !<
	1413	INTEGER(iwp) :: j !<
	1414	INTEGER(iwp) :: k !<
[1115]	1415
[1682]	1416	REAL(wp) :: t_surface !<
[1320]	1417
[1361]	1418	IF ( large_scale_forcing .AND. lsf_surf ) THEN
[1241]	1419	!
	1420	!-- Calculate:
	1421	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
	1422	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
	1423	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
[1353]	1424	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
[1241]	1425	DO k = nzb, nzt+1
[1353]	1426	hyp(k) = surface_pressure * 100.0_wp * &
[1361]	1427	( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
[1353]	1428	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
	1429	t_d_pt(k) = 1.0_wp / pt_d_t(k)
[2155]	1430	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
[1241]	1431	ENDDO
	1432	!
	1433	!-- Compute reference density
[1353]	1434	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
[1241]	1435	ENDIF
	1436
[1361]	1437	!
[2155]	1438	!-- Compute length of time step
[1361]	1439	IF ( call_microphysics_at_all_substeps ) THEN
	1440	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
	1441	ELSE
	1442	dt_micro = dt_3d
	1443	ENDIF
[1241]	1444
[1115]	1445	!
[1361]	1446	!-- Use 1d arrays
[1115]	1447	q_1d(:) = q(:,j,i)
	1448	pt_1d(:) = pt(:,j,i)
	1449	qc_1d(:) = qc(:,j,i)
	1450	nc_1d(:) = nc_const
[1822]	1451	IF ( microphysics_seifert ) THEN
[1115]	1452	qr_1d(:) = qr(:,j,i)
	1453	nr_1d(:) = nr(:,j,i)
	1454	ENDIF
[1361]	1455
[1115]	1456	!
[1822]	1457	!-- Reset precipitation rate
	1458	IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp
	1459
	1460	!
[1115]	1461	!-- Compute cloud physics
[1822]	1462	IF( microphysics_kessler ) THEN
	1463
	1464	CALL autoconversion_kessler( i,j )
[1831]	1465	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud( i,j )
[1822]	1466
	1467	ELSEIF ( microphysics_seifert ) THEN
	1468
	1469	CALL adjust_cloud( i,j )
[1115]	1470	CALL autoconversion( i,j )
	1471	CALL accretion( i,j )
	1472	CALL selfcollection_breakup( i,j )
	1473	CALL evaporation_rain( i,j )
	1474	CALL sedimentation_rain( i,j )
[1831]	1475	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud( i,j )
[1115]	1476
[1691]	1477	ENDIF
	1478
[1822]	1479	CALL calc_precipitation_amount( i,j )
	1480
[1115]	1481	!
[1361]	1482	!-- Store results on the 3d arrays
	1483	q(:,j,i) = q_1d(:)
	1484	pt(:,j,i) = pt_1d(:)
[1822]	1485	IF ( microphysics_seifert ) THEN
[1361]	1486	qr(:,j,i) = qr_1d(:)
	1487	nr(:,j,i) = nr_1d(:)
[1115]	1488	ENDIF
	1489
	1490	END SUBROUTINE microphysics_control_ij
	1491
[1682]	1492	!------------------------------------------------------------------------------!
	1493	! Description:
	1494	! ------------
[2155]	1495	!> Adjust number of raindrops to avoid nonlinear effects in
[1682]	1496	!> sedimentation and evaporation of rain drops due to too small or
	1497	!> too big weights of rain drops (Stevens and Seifert, 2008).
	1498	!> The same procedure is applied to cloud droplets if they are determined
	1499	!> prognostically. Call for grid point i,j
	1500	!------------------------------------------------------------------------------!
[1115]	1501	SUBROUTINE adjust_cloud_ij( i, j )
	1502
[1320]	1503	USE cloud_parameters, &
[1849]	1504	ONLY: hyrho
[1320]	1505
	1506	USE indices, &
[2232]	1507	ONLY: nzb, nzt, wall_flags_0
[1320]	1508
	1509	USE kinds
	1510
[1115]	1511	IMPLICIT NONE
	1512
[1682]	1513	INTEGER(iwp) :: i !<
	1514	INTEGER(iwp) :: j !<
	1515	INTEGER(iwp) :: k !<
	1516
[2232]	1517	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1022]	1518
[2232]	1519	DO k = nzb+1, nzt
	1520	!
	1521	!-- Predetermine flag to mask topography
	1522	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1523
[1361]	1524	IF ( qr_1d(k) <= eps_sb ) THEN
	1525	qr_1d(k) = 0.0_wp
	1526	nr_1d(k) = 0.0_wp
[1065]	1527	ELSE
[1022]	1528	!
[2155]	1529	!-- Adjust number of raindrops to avoid nonlinear effects in
[1048]	1530	!-- sedimentation and evaporation of rain drops due to too small or
[1065]	1531	!-- too big weights of rain drops (Stevens and Seifert, 2008).
[1361]	1532	IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) ) THEN
[2232]	1533	nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin * flag
[1361]	1534	ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) ) THEN
[2232]	1535	nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax * flag
[1048]	1536	ENDIF
[1115]	1537
[1022]	1538	ENDIF
[1115]	1539
[1022]	1540	ENDDO
	1541
[1115]	1542	END SUBROUTINE adjust_cloud_ij
[1022]	1543
[1106]	1544
[1682]	1545	!------------------------------------------------------------------------------!
	1546	! Description:
	1547	! ------------
	1548	!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
	1549	!------------------------------------------------------------------------------!
[1005]	1550	SUBROUTINE autoconversion_ij( i, j )
[1000]	1551
[1320]	1552	USE arrays_3d, &
[1849]	1553	ONLY: diss, dzu
[1115]	1554
[1320]	1555	USE cloud_parameters, &
[1849]	1556	ONLY: hyrho
[1320]	1557
	1558	USE control_parameters, &
[1849]	1559	ONLY: rho_surface
[1320]	1560
	1561	USE grid_variables, &
	1562	ONLY: dx, dy
	1563
	1564	USE indices, &
[2232]	1565	ONLY: nzb, nzt, wall_flags_0
[1320]	1566
	1567	USE kinds
	1568
[1000]	1569	IMPLICIT NONE
	1570
[1682]	1571	INTEGER(iwp) :: i !<
	1572	INTEGER(iwp) :: j !<
	1573	INTEGER(iwp) :: k !<
[1000]	1574
[2155]	1575	REAL(wp) :: alpha_cc !<
[1682]	1576	REAL(wp) :: autocon !<
	1577	REAL(wp) :: dissipation !<
[2232]	1578	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1682]	1579	REAL(wp) :: k_au !<
	1580	REAL(wp) :: l_mix !<
	1581	REAL(wp) :: nu_c !<
	1582	REAL(wp) :: phi_au !<
	1583	REAL(wp) :: r_cc !<
	1584	REAL(wp) :: rc !<
	1585	REAL(wp) :: re_lambda !<
	1586	REAL(wp) :: sigma_cc !<
	1587	REAL(wp) :: tau_cloud !<
	1588	REAL(wp) :: xc !<
[1106]	1589
[2232]	1590	DO k = nzb+1, nzt
	1591	!
	1592	!-- Predetermine flag to mask topography
	1593	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1000]	1594
[1115]	1595	IF ( qc_1d(k) > eps_sb ) THEN
[1361]	1596
	1597	k_au = k_cc / ( 20.0_wp * x0 )
[1012]	1598	!
[1048]	1599	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[1353]	1600	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
	1601	tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
[1012]	1602	!
[2155]	1603	!-- Universal function for autoconversion process
[1012]	1604	!-- (Seifert and Beheng, 2006):
[1361]	1605	phi_au = 600.0_wp * tau_cloud*0.68_wp ( 1.0_wp - tau_cloud0.68_wp )3
[1012]	1606	!
	1607	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
[1353]	1608	!-- (Use constant nu_c = 1.0_wp instead?)
[1361]	1609	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp )
[1012]	1610	!
	1611	!-- Mean weight of cloud droplets:
[1115]	1612	xc = hyrho(k) * qc_1d(k) / nc_1d(k)
[1012]	1613	!
[2155]	1614	!-- Parameterized turbulence effects on autoconversion (Seifert,
[1065]	1615	!-- Nuijens and Stevens, 2010)
[1831]	1616	IF ( collision_turbulence ) THEN
[1065]	1617	!
	1618	!-- Weight averaged radius of cloud droplets:
[1353]	1619	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
[1065]	1620
[1353]	1621	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
	1622	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
	1623	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
[1065]	1624	!
	1625	!-- Mixing length (neglecting distance to ground and stratification)
[1334]	1626	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
[1065]	1627	!
[2155]	1628	!-- Limit dissipation rate according to Seifert, Nuijens and
[1065]	1629	!-- Stevens (2010)
[1361]	1630	dissipation = MIN( 0.06_wp, diss(k,j,i) )
[1065]	1631	!
	1632	!-- Compute Taylor-microscale Reynolds number:
[1361]	1633	re_lambda = 6.0_wp / 11.0_wp * &
	1634	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
	1635	SQRT( 15.0_wp / kin_vis_air ) * &
	1636	dissipation**( 1.0_wp / 6.0_wp )
[1065]	1637	!
	1638	!-- The factor of 1.0E4 is needed to convert the dissipation rate
	1639	!-- from m2 s-3 to cm2 s-3.
[1361]	1640	k_au = k_au * ( 1.0_wp + &
	1641	dissipation * 1.0E4_wp * &
	1642	( re_lambda * 1.0E-3_wp )*0.25_wp &
	1643	( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) / &
	1644	sigma_cc )**2 &
	1645	) + beta_cc &
	1646	) &
	1647	)
[1065]	1648	ENDIF
	1649	!
[1012]	1650	!-- Autoconversion rate (Seifert and Beheng, 2006):
[1361]	1651	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
	1652	( nu_c + 1.0_wp )*2 qc_1d(k)*2 xc*2 &
	1653	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
[1115]	1654	rho_surface
	1655	autocon = MIN( autocon, qc_1d(k) / dt_micro )
[1106]	1656
[2232]	1657	qr_1d(k) = qr_1d(k) + autocon * dt_micro * flag
	1658	qc_1d(k) = qc_1d(k) - autocon * dt_micro * flag
	1659	nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro * flag
[1115]	1660
[1005]	1661	ENDIF
[1000]	1662
	1663	ENDDO
	1664
[1005]	1665	END SUBROUTINE autoconversion_ij
	1666
[1822]	1667	!------------------------------------------------------------------------------!
	1668	! Description:
	1669	! ------------
	1670	!> Autoconversion process (Kessler, 1969).
	1671	!------------------------------------------------------------------------------!
	1672	SUBROUTINE autoconversion_kessler_ij( i, j )
[1106]	1673
[1822]	1674	USE arrays_3d, &
[1849]	1675	ONLY: dzw, prr
[1822]	1676
	1677	USE cloud_parameters, &
[1849]	1678	ONLY: l_d_cp, pt_d_t
[1822]	1679
	1680	USE indices, &
[2232]	1681	ONLY: nzb, nzb_max, nzt, wall_flags_0
[1822]	1682
	1683	USE kinds
	1684
	1685
	1686	IMPLICIT NONE
	1687
[2232]	1688	INTEGER(iwp) :: i !<
	1689	INTEGER(iwp) :: j !<
	1690	INTEGER(iwp) :: k !<
	1691	INTEGER(iwp) :: k_wall !< topography top index
[1822]	1692
	1693	REAL(wp) :: dqdt_precip !<
[2232]	1694	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1822]	1695
[2232]	1696	!
	1697	!-- Determine vertical index of topography top
	1698	k_wall = MAXLOC( &
	1699	MERGE( 1, 0, &
	1700	BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 ) &
	1701	), DIM = 1 &
	1702	) - 1
	1703	DO k = nzb+1, nzt
	1704	!
	1705	!-- Predetermine flag to mask topography
	1706	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1822]	1707
	1708	IF ( qc_1d(k) > ql_crit ) THEN
	1709	dqdt_precip = prec_time_const * ( qc_1d(k) - ql_crit )
	1710	ELSE
	1711	dqdt_precip = 0.0_wp
	1712	ENDIF
	1713
[2232]	1714	qc_1d(k) = qc_1d(k) - dqdt_precip * dt_micro * flag
	1715	q_1d(k) = q_1d(k) - dqdt_precip * dt_micro * flag
	1716	pt_1d(k) = pt_1d(k) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k) * flag
[1822]	1717
	1718	!
[1845]	1719	!-- Compute the rain rate (stored on surface grid point)
[2232]	1720	prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
[1822]	1721
	1722	ENDDO
	1723
	1724	END SUBROUTINE autoconversion_kessler_ij
	1725
[1682]	1726	!------------------------------------------------------------------------------!
	1727	! Description:
	1728	! ------------
	1729	!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
	1730	!------------------------------------------------------------------------------!
[1005]	1731	SUBROUTINE accretion_ij( i, j )
	1732
[1320]	1733	USE arrays_3d, &
[1849]	1734	ONLY: diss
[1115]	1735
[1320]	1736	USE cloud_parameters, &
[1849]	1737	ONLY: hyrho
[1320]	1738
	1739	USE control_parameters, &
[1849]	1740	ONLY: rho_surface
[1320]	1741
	1742	USE indices, &
[2232]	1743	ONLY: nzb, nzt, wall_flags_0
[1320]	1744
	1745	USE kinds
	1746
[1005]	1747	IMPLICIT NONE
	1748
[1682]	1749	INTEGER(iwp) :: i !<
	1750	INTEGER(iwp) :: j !<
	1751	INTEGER(iwp) :: k !<
[1005]	1752
[1682]	1753	REAL(wp) :: accr !<
[2232]	1754	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1682]	1755	REAL(wp) :: k_cr !<
	1756	REAL(wp) :: phi_ac !<
	1757	REAL(wp) :: tau_cloud !<
[1320]	1758
[2232]	1759	DO k = nzb+1, nzt
	1760	!
	1761	!-- Predetermine flag to mask topography
	1762	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1763
[1115]	1764	IF ( ( qc_1d(k) > eps_sb ) .AND. ( qr_1d(k) > eps_sb ) ) THEN
[1012]	1765	!
[1048]	1766	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[2155]	1767	tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
[1012]	1768	!
[2155]	1769	!-- Universal function for accretion process
[1048]	1770	!-- (Seifert and Beheng, 2001):
[1361]	1771	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
[1012]	1772	!
[2155]	1773	!-- Parameterized turbulence effects on autoconversion (Seifert,
	1774	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
[1361]	1775	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
[1831]	1776	IF ( collision_turbulence ) THEN
[1361]	1777	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
	1778	MIN( 600.0_wp, &
	1779	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
	1780	)
[1065]	1781	ELSE
[2155]	1782	k_cr = k_cr0
[1065]	1783	ENDIF
	1784	!
[1012]	1785	!-- Accretion rate (Seifert and Beheng, 2006):
[1361]	1786	accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * SQRT( rho_surface * hyrho(k) )
[1115]	1787	accr = MIN( accr, qc_1d(k) / dt_micro )
[1106]	1788
[2232]	1789	qr_1d(k) = qr_1d(k) + accr * dt_micro * flag
	1790	qc_1d(k) = qc_1d(k) - accr * dt_micro * flag
[1115]	1791
[1005]	1792	ENDIF
[1106]	1793
[1005]	1794	ENDDO
	1795
[1000]	1796	END SUBROUTINE accretion_ij
	1797
[1005]	1798
[1682]	1799	!------------------------------------------------------------------------------!
	1800	! Description:
	1801	! ------------
	1802	!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
	1803	!------------------------------------------------------------------------------!
[1005]	1804	SUBROUTINE selfcollection_breakup_ij( i, j )
	1805
[1320]	1806	USE cloud_parameters, &
[1849]	1807	ONLY: hyrho
[1320]	1808
	1809	USE control_parameters, &
[1849]	1810	ONLY: rho_surface
[1320]	1811
	1812	USE indices, &
[2232]	1813	ONLY: nzb, nzt, wall_flags_0
[1320]	1814
	1815	USE kinds
[2155]	1816
[1005]	1817	IMPLICIT NONE
	1818
[1682]	1819	INTEGER(iwp) :: i !<
	1820	INTEGER(iwp) :: j !<
	1821	INTEGER(iwp) :: k !<
[1005]	1822
[1682]	1823	REAL(wp) :: breakup !<
	1824	REAL(wp) :: dr !<
[2232]	1825	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1682]	1826	REAL(wp) :: phi_br !<
	1827	REAL(wp) :: selfcoll !<
[1320]	1828
[2232]	1829	DO k = nzb+1, nzt
	1830	!
	1831	!-- Predetermine flag to mask topography
	1832	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1833
[1115]	1834	IF ( qr_1d(k) > eps_sb ) THEN
[1012]	1835	!
[1115]	1836	!-- Selfcollection rate (Seifert and Beheng, 2001):
[1361]	1837	selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
[1012]	1838	!
[1115]	1839	!-- Weight averaged diameter of rain drops:
[1334]	1840	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	1841	!
[1048]	1842	!-- Collisional breakup rate (Seifert, 2008):
[1353]	1843	IF ( dr >= 0.3E-3_wp ) THEN
	1844	phi_br = k_br * ( dr - 1.1E-3_wp )
	1845	breakup = selfcoll * ( phi_br + 1.0_wp )
[1005]	1846	ELSE
[1353]	1847	breakup = 0.0_wp
[1005]	1848	ENDIF
[1048]	1849
[1115]	1850	selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
[2232]	1851	nr_1d(k) = nr_1d(k) + selfcoll * dt_micro * flag
[1106]	1852
[2155]	1853	ENDIF
[1005]	1854	ENDDO
	1855
	1856	END SUBROUTINE selfcollection_breakup_ij
	1857
[1106]	1858
[1682]	1859	!------------------------------------------------------------------------------!
	1860	! Description:
	1861	! ------------
[2155]	1862	!> Evaporation of precipitable water. Condensation is neglected for
[1682]	1863	!> precipitable water. Call for grid point i,j
	1864	!------------------------------------------------------------------------------!
[1012]	1865	SUBROUTINE evaporation_rain_ij( i, j )
	1866
[1320]	1867	USE arrays_3d, &
[1849]	1868	ONLY: hyp
[1048]	1869
[1320]	1870	USE cloud_parameters, &
[1849]	1871	ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
[1320]	1872
	1873	USE constants, &
	1874	ONLY: pi
	1875
	1876	USE indices, &
[2232]	1877	ONLY: nzb, nzt, wall_flags_0
[1320]	1878
	1879	USE kinds
	1880
[1012]	1881	IMPLICIT NONE
	1882
[1682]	1883	INTEGER(iwp) :: i !<
	1884	INTEGER(iwp) :: j !<
	1885	INTEGER(iwp) :: k !<
[1012]	1886
[1682]	1887	REAL(wp) :: alpha !<
	1888	REAL(wp) :: dr !<
	1889	REAL(wp) :: e_s !<
	1890	REAL(wp) :: evap !<
	1891	REAL(wp) :: evap_nr !<
	1892	REAL(wp) :: f_vent !<
[2232]	1893	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1682]	1894	REAL(wp) :: g_evap !<
	1895	REAL(wp) :: lambda_r !<
	1896	REAL(wp) :: mu_r !<
	1897	REAL(wp) :: mu_r_2 !<
	1898	REAL(wp) :: mu_r_5d2 !<
	1899	REAL(wp) :: nr_0 !<
	1900	REAL(wp) :: q_s !<
	1901	REAL(wp) :: sat !<
	1902	REAL(wp) :: t_l !<
	1903	REAL(wp) :: temp !<
	1904	REAL(wp) :: xr !<
[1320]	1905
[2232]	1906	DO k = nzb+1, nzt
	1907	!
	1908	!-- Predetermine flag to mask topography
	1909	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1910
[1115]	1911	IF ( qr_1d(k) > eps_sb ) THEN
[1012]	1912	!
	1913	!-- Actual liquid water temperature:
[1115]	1914	t_l = t_d_pt(k) * pt_1d(k)
[1012]	1915	!
	1916	!-- Saturation vapor pressure at t_l:
[1361]	1917	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
	1918	( t_l - 35.86_wp ) &
	1919	)
[1012]	1920	!
	1921	!-- Computation of saturation humidity:
[1361]	1922	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
[1353]	1923	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
[1361]	1924	q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
[1012]	1925	!
[1106]	1926	!-- Supersaturation:
[1361]	1927	sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
[1012]	1928	!
[1361]	1929	!-- Evaporation needs only to be calculated in subsaturated regions
	1930	IF ( sat < 0.0_wp ) THEN
[1012]	1931	!
[1361]	1932	!-- Actual temperature:
	1933	temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
[2155]	1934
[1361]	1935	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v / &
	1936	( thermal_conductivity_l * temp ) + &
	1937	r_v * temp / ( diff_coeff_l * e_s ) &
	1938	)
[1012]	1939	!
[1361]	1940	!-- Mean weight of rain drops
	1941	xr = hyrho(k) * qr_1d(k) / nr_1d(k)
[1115]	1942	!
[1361]	1943	!-- Weight averaged diameter of rain drops:
	1944	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	1945	!
[2155]	1946	!-- Compute ventilation factor and intercept parameter
[1361]	1947	!-- (Seifert and Beheng, 2006; Seifert, 2008):
	1948	IF ( ventilation_effect ) THEN
[1115]	1949	!
[1361]	1950	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	1951	!-- Stevens and Seifert, 2008):
	1952	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
	1953	!
	1954	!-- Slope parameter of gamma distribution (Seifert, 2008):
	1955	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	1956	( mu_r + 1.0_wp ) &
	1957	)**( 1.0_wp / 3.0_wp ) / dr
[1115]	1958
[1361]	1959	mu_r_2 = mu_r + 2.0_wp
[2155]	1960	mu_r_5d2 = mu_r + 2.5_wp
[1361]	1961
[2155]	1962	f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) + &
[1361]	1963	b_vent * schmidt_p_1d3 * &
	1964	SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * &
	1965	lambda_r*( -mu_r_5d2 ) &
	1966	( 1.0_wp - &
	1967	0.5_wp * ( b_term / a_term ) * &
	1968	( lambda_r / ( c_term + lambda_r ) &
	1969	)**mu_r_5d2 - &
	1970	0.125_wp * ( b_term / a_term )*2 &
	1971	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
	1972	)**mu_r_5d2 - &
	1973	0.0625_wp * ( b_term / a_term )*3 &
	1974	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
	1975	)**mu_r_5d2 - &
[2155]	1976	0.0390625_wp * ( b_term / a_term )*4 &
[1361]	1977	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
	1978	)**mu_r_5d2 &
	1979	)
	1980
	1981	nr_0 = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) / &
[2155]	1982	gamm( mu_r + 1.0_wp )
[1361]	1983	ELSE
	1984	f_vent = 1.0_wp
	1985	nr_0 = nr_1d(k) * dr
	1986	ENDIF
[1012]	1987	!
[1361]	1988	!-- Evaporation rate of rain water content (Seifert and Beheng, 2006):
	1989	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
	1990	evap = MAX( evap, -qr_1d(k) / dt_micro )
	1991	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
	1992	-nr_1d(k) / dt_micro )
[1106]	1993
[2232]	1994	qr_1d(k) = qr_1d(k) + evap * dt_micro * flag
	1995	nr_1d(k) = nr_1d(k) + evap_nr * dt_micro * flag
[1115]	1996
[1361]	1997	ENDIF
[2155]	1998	ENDIF
[1106]	1999
[1012]	2000	ENDDO
	2001
	2002	END SUBROUTINE evaporation_rain_ij
	2003
[1106]	2004
[1682]	2005	!------------------------------------------------------------------------------!
	2006	! Description:
	2007	! ------------
[2155]	2008	!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
[1682]	2009	!> Call for grid point i,j
	2010	!------------------------------------------------------------------------------!
[1012]	2011	SUBROUTINE sedimentation_cloud_ij( i, j )
	2012
[1320]	2013	USE arrays_3d, &
[1849]	2014	ONLY: ddzu, dzu, prr
[1320]	2015
	2016	USE cloud_parameters, &
[1849]	2017	ONLY: hyrho, l_d_cp, pt_d_t
[1320]	2018
	2019	USE control_parameters, &
[1849]	2020	ONLY: call_microphysics_at_all_substeps, intermediate_timestep_count
[1320]	2021
	2022	USE indices, &
[2232]	2023	ONLY: nzb, nzb, nzt, wall_flags_0
[1320]	2024
	2025	USE kinds
[2155]	2026
[1691]	2027	USE statistics, &
	2028	ONLY: weight_substep
	2029
[1012]	2030	IMPLICIT NONE
	2031
[1849]	2032	INTEGER(iwp) :: i !<
	2033	INTEGER(iwp) :: j !<
	2034	INTEGER(iwp) :: k !<
[1106]	2035
[2232]	2036	REAL(wp) :: flag !< flag to indicate first grid level above surface
	2037	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
[1115]	2038
[1353]	2039	sed_qc(nzt+1) = 0.0_wp
[1012]	2040
[2232]	2041	DO k = nzt, nzb+1, -1
	2042	!
	2043	!-- Predetermine flag to mask topography
	2044	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2045
[1115]	2046	IF ( qc_1d(k) > eps_sb ) THEN
[1361]	2047	sed_qc(k) = sed_qc_const * nc_1d(k)*( -2.0_wp / 3.0_wp ) &
[2232]	2048	( qc_1d(k) * hyrho(k) )*( 5.0_wp / 3.0_wp ) flag
[1115]	2049	ELSE
[1353]	2050	sed_qc(k) = 0.0_wp
[1012]	2051	ENDIF
[1115]	2052
[1361]	2053	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) / &
	2054	dt_micro + sed_qc(k+1) &
[2232]	2055	) * flag
[1115]	2056
[1361]	2057	q_1d(k) = q_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
[2232]	2058	hyrho(k) * dt_micro * flag
[2155]	2059	qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
[2232]	2060	hyrho(k) * dt_micro * flag
[1361]	2061	pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
[2232]	2062	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro * flag
[1115]	2063
[1691]	2064	!
	2065	!-- Compute the precipitation rate of cloud (fog) droplets
[1822]	2066	IF ( call_microphysics_at_all_substeps ) THEN
[2232]	2067	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * &
	2068	weight_substep(intermediate_timestep_count) * flag
[1822]	2069	ELSE
[2232]	2070	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
[1691]	2071	ENDIF
	2072
[1012]	2073	ENDDO
	2074
	2075	END SUBROUTINE sedimentation_cloud_ij
	2076
[1106]	2077
[1682]	2078	!------------------------------------------------------------------------------!
	2079	! Description:
	2080	! ------------
	2081	!> Computation of sedimentation flux. Implementation according to Stevens
	2082	!> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j
	2083	!------------------------------------------------------------------------------!
[1012]	2084	SUBROUTINE sedimentation_rain_ij( i, j )
	2085
[1320]	2086	USE arrays_3d, &
[1849]	2087	ONLY: ddzu, dzu, prr
[1320]	2088
	2089	USE cloud_parameters, &
[1849]	2090	ONLY: hyrho, l_d_cp, pt_d_t
[1320]	2091
	2092	USE control_parameters, &
[1849]	2093	ONLY: call_microphysics_at_all_substeps, intermediate_timestep_count
[1320]	2094
	2095	USE indices, &
[2232]	2096	ONLY: nzb, nzb, nzt, wall_flags_0
[1320]	2097
	2098	USE kinds
	2099
	2100	USE statistics, &
	2101	ONLY: weight_substep
[2155]	2102
[2232]	2103	USE surface_mod, &
	2104	ONLY : bc_h
	2105
[1012]	2106	IMPLICIT NONE
	2107
[2232]	2108	INTEGER(iwp) :: i !< running index x direction
	2109	INTEGER(iwp) :: j !< running index y direction
	2110	INTEGER(iwp) :: k !< running index z direction
	2111	INTEGER(iwp) :: k_run !<
	2112	INTEGER(iwp) :: m !< running index surface elements
	2113	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	2114	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1012]	2115
[1682]	2116	REAL(wp) :: c_run !<
	2117	REAL(wp) :: d_max !<
	2118	REAL(wp) :: d_mean !<
	2119	REAL(wp) :: d_min !<
	2120	REAL(wp) :: dr !<
	2121	REAL(wp) :: flux !<
[2232]	2122	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1682]	2123	REAL(wp) :: lambda_r !<
	2124	REAL(wp) :: mu_r !<
	2125	REAL(wp) :: z_run !<
[1320]	2126
[1682]	2127	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<
	2128	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<
	2129	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<
	2130	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<
	2131	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<
	2132	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<
	2133	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<
	2134	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<
[1320]	2135
[1012]	2136	!
[2155]	2137	!-- Compute velocities
[2232]	2138	DO k = nzb+1, nzt
	2139	!
	2140	!-- Predetermine flag to mask topography
	2141	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2142
[1115]	2143	IF ( qr_1d(k) > eps_sb ) THEN
	2144	!
	2145	!-- Weight averaged diameter of rain drops:
[1334]	2146	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	2147	!
	2148	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	2149	!-- Stevens and Seifert, 2008):
[1353]	2150	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
[1115]	2151	!
	2152	!-- Slope parameter of gamma distribution (Seifert, 2008):
[1361]	2153	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	2154	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
[1115]	2155
[1361]	2156	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	2157	a_term - b_term * ( 1.0_wp + &
	2158	c_term / lambda_r )*( -1.0_wp &
	2159	( mu_r + 1.0_wp ) ) &
	2160	) &
[2232]	2161	) * flag
[1361]	2162	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	2163	a_term - b_term * ( 1.0_wp + &
	2164	c_term / lambda_r )*( -1.0_wp &
	2165	( mu_r + 4.0_wp ) ) &
	2166	) &
[2232]	2167	) * flag
[1065]	2168	ELSE
[1353]	2169	w_nr(k) = 0.0_wp
	2170	w_qr(k) = 0.0_wp
[1065]	2171	ENDIF
	2172	ENDDO
[1048]	2173	!
[2232]	2174	!-- Adjust boundary values using surface data type.
	2175	!-- Upward facing non-natural
	2176	surf_s = bc_h(0)%start_index(j,i)
	2177	surf_e = bc_h(0)%end_index(j,i)
	2178	DO m = surf_s, surf_e
	2179	k = bc_h(0)%k(m)
	2180	w_nr(k-1) = w_nr(k)
	2181	w_qr(k-1) = w_qr(k)
	2182	ENDDO
	2183	!
	2184	!-- Downward facing non-natural
	2185	surf_s = bc_h(1)%start_index(j,i)
	2186	surf_e = bc_h(1)%end_index(j,i)
	2187	DO m = surf_s, surf_e
	2188	k = bc_h(1)%k(m)
	2189	w_nr(k+1) = w_nr(k)
	2190	w_qr(k+1) = w_qr(k)
	2191	ENDDO
	2192	!
	2193	!-- Neumann boundary condition at model top
[1353]	2194	w_nr(nzt+1) = 0.0_wp
	2195	w_qr(nzt+1) = 0.0_wp
[1065]	2196	!
	2197	!-- Compute Courant number
[2232]	2198	DO k = nzb+1, nzt
	2199	!
	2200	!-- Predetermine flag to mask topography
	2201	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2202
[1361]	2203	c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) * &
[2232]	2204	dt_micro * ddzu(k) * flag
[1361]	2205	c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * &
[2232]	2206	dt_micro * ddzu(k) * flag
[2155]	2207	ENDDO
[1065]	2208	!
	2209	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
	2210	IF ( limiter_sedimentation ) THEN
	2211
[2232]	2212	DO k = nzb+1, nzt
	2213	!
	2214	!-- Predetermine flag to mask topography
	2215	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2216
[1646]	2217	d_mean = 0.5_wp * ( qr_1d(k+1) - qr_1d(k-1) )
[1115]	2218	d_min = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
	2219	d_max = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
[1065]	2220
[1361]	2221	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	2222	2.0_wp * d_max, &
[2232]	2223	ABS( d_mean ) ) * flag
[1065]	2224
[1646]	2225	d_mean = 0.5_wp * ( nr_1d(k+1) - nr_1d(k-1) )
[1115]	2226	d_min = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
	2227	d_max = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
[1065]	2228
[1361]	2229	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	2230	2.0_wp * d_max, &
[2232]	2231	ABS( d_mean ) ) * flag
[1022]	2232	ENDDO
[1048]	2233
[1065]	2234	ELSE
[1106]	2235
[1353]	2236	nr_slope = 0.0_wp
	2237	qr_slope = 0.0_wp
[1106]	2238
[1065]	2239	ENDIF
[1115]	2240
[1353]	2241	sed_nr(nzt+1) = 0.0_wp
	2242	sed_qr(nzt+1) = 0.0_wp
[1065]	2243	!
	2244	!-- Compute sedimentation flux
[2232]	2245	DO k = nzt, nzb+1, -1
[1065]	2246	!
[2232]	2247	!-- Predetermine flag to mask topography
	2248	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2249	!
	2250	!-- Sum up all rain drop number densities which contribute to the flux
[1065]	2251	!-- through k-1/2
[1353]	2252	flux = 0.0_wp
	2253	z_run = 0.0_wp ! height above z(k)
[1065]	2254	k_run = k
[1346]	2255	c_run = MIN( 1.0_wp, c_nr(k) )
[1353]	2256	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
[1361]	2257	flux = flux + hyrho(k_run) * &
	2258	( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) * &
[2232]	2259	0.5_wp ) * c_run * dzu(k_run) * flag
	2260	z_run = z_run + dzu(k_run) * flag
	2261	k_run = k_run + 1 * flag
	2262	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) * flag
[1022]	2263	ENDDO
	2264	!
[2155]	2265	!-- It is not allowed to sediment more rain drop number density than
[1065]	2266	!-- available
[1361]	2267	flux = MIN( flux, &
[1115]	2268	hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
[1065]	2269
[2232]	2270	sed_nr(k) = flux / dt_micro * flag
[1361]	2271	nr_1d(k) = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / &
[2232]	2272	hyrho(k) * dt_micro * flag
[1065]	2273	!
[2155]	2274	!-- Sum up all rain water content which contributes to the flux
[1065]	2275	!-- through k-1/2
[1353]	2276	flux = 0.0_wp
	2277	z_run = 0.0_wp ! height above z(k)
[1065]	2278	k_run = k
[1346]	2279	c_run = MIN( 1.0_wp, c_qr(k) )
[1106]	2280
[1361]	2281	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
[1106]	2282
[1361]	2283	flux = flux + hyrho(k_run) * &
	2284	( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) * &
[2232]	2285	0.5_wp ) * c_run * dzu(k_run) * flag
	2286	z_run = z_run + dzu(k_run) * flag
	2287	k_run = k_run + 1 * flag
	2288	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) * flag
[1106]	2289
[1065]	2290	ENDDO
	2291	!
	2292	!-- It is not allowed to sediment more rain water content than available
[1361]	2293	flux = MIN( flux, &
[1115]	2294	hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
[1065]	2295
[2232]	2296	sed_qr(k) = flux / dt_micro * flag
[1115]	2297
[1361]	2298	qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
[2232]	2299	hyrho(k) * dt_micro * flag
[1361]	2300	q_1d(k) = q_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
[2232]	2301	hyrho(k) * dt_micro * flag
[1361]	2302	pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
[2232]	2303	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro * flag
[1065]	2304	!
	2305	!-- Compute the rain rate
[1361]	2306	IF ( call_microphysics_at_all_substeps ) THEN
[1691]	2307	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) &
[2232]	2308	* weight_substep(intermediate_timestep_count) * flag
[1361]	2309	ELSE
[2232]	2310	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
[1361]	2311	ENDIF
	2312
[1065]	2313	ENDDO
[1115]	2314
[1691]	2315	END SUBROUTINE sedimentation_rain_ij
[1012]	2316
[1691]	2317
	2318	!------------------------------------------------------------------------------!
	2319	! Description:
	2320	! ------------
	2321	!> This subroutine computes the precipitation amount due to gravitational
	2322	!> settling of rain and cloud (fog) droplets
	2323	!------------------------------------------------------------------------------!
	2324	SUBROUTINE calc_precipitation_amount_ij( i, j )
	2325
[1849]	2326	USE arrays_3d, &
	2327	ONLY: precipitation_amount, prr
	2328
[1691]	2329	USE cloud_parameters, &
[1849]	2330	ONLY: hyrho
[1691]	2331
	2332	USE control_parameters, &
	2333	ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d, &
	2334	intermediate_timestep_count, intermediate_timestep_count_max,&
[1822]	2335	precipitation_amount_interval, time_do2d_xy
[1691]	2336
	2337	USE indices, &
[2232]	2338	ONLY: nzb, nzt, wall_flags_0
[1691]	2339
	2340	USE kinds
	2341
[2232]	2342	USE surface_mod, &
	2343	ONLY : bc_h
	2344
[1691]	2345	IMPLICIT NONE
	2346
[2232]	2347	INTEGER(iwp) :: i !< running index x direction
	2348	INTEGER(iwp) :: j !< running index y direction
	2349	INTEGER(iwp) :: k !< running index z direction
	2350	INTEGER(iwp) :: m !< running index surface elements
	2351	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	2352	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1691]	2353
	2354	IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
	2355	( .NOT. call_microphysics_at_all_substeps .OR. &
	2356	intermediate_timestep_count == intermediate_timestep_count_max ) ) &
	2357	THEN
	2358
[2232]	2359	surf_s = bc_h(0)%start_index(j,i)
	2360	surf_e = bc_h(0)%end_index(j,i)
	2361	DO m = surf_s, surf_e
	2362	k = bc_h(0)%k(m)
	2363	precipitation_amount(j,i) = precipitation_amount(j,i) + &
	2364	prr(k,j,i) * hyrho(k) * dt_3d
	2365	ENDDO
	2366
[1048]	2367	ENDIF
	2368
[1691]	2369	END SUBROUTINE calc_precipitation_amount_ij
[1012]	2370
[1361]	2371	!------------------------------------------------------------------------------!
[1682]	2372	! Description:
	2373	! ------------
	2374	!> This function computes the gamma function (Press et al., 1992).
[2155]	2375	!> The gamma function is needed for the calculation of the evaporation
[1682]	2376	!> of rain drops.
[1361]	2377	!------------------------------------------------------------------------------!
[2155]	2378	FUNCTION gamm( xx )
[1320]	2379
	2380	USE kinds
	2381
[2155]	2382	IMPLICIT NONE
[1106]	2383
[1682]	2384	INTEGER(iwp) :: j !<
[1320]	2385
[1682]	2386	REAL(wp) :: gamm !<
	2387	REAL(wp) :: ser !<
	2388	REAL(wp) :: tmp !<
	2389	REAL(wp) :: x_gamm !<
	2390	REAL(wp) :: xx !<
	2391	REAL(wp) :: y_gamm !<
[1320]	2392
[1849]	2393
	2394	REAL(wp), PARAMETER :: stp = 2.5066282746310005_wp !<
	2395	REAL(wp), PARAMETER :: cof(6) = (/ 76.18009172947146_wp, &
	2396	-86.50532032941677_wp, &
	2397	24.01409824083091_wp, &
	2398	-1.231739572450155_wp, &
	2399	0.1208650973866179E-2_wp, &
	2400	-0.5395239384953E-5_wp /) !<
	2401
	2402	x_gamm = xx
	2403	y_gamm = x_gamm
[1353]	2404	tmp = x_gamm + 5.5_wp
[1849]	2405	tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
[1334]	2406	ser = 1.000000000190015_wp
[1106]	2407
[2155]	2408	DO j = 1, 6
	2409	y_gamm = y_gamm + 1.0_wp
	2410	ser = ser + cof( j ) / y_gamm
[1106]	2411	ENDDO
	2412
[2155]	2413	!
	2414	!-- Until this point the algorithm computes the logarithm of the gamma
	2415	!-- function. Hence, the exponential function is used.
	2416	! gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
	2417	gamm = EXP( tmp ) * stp * ser / x_gamm
[1106]	2418
[2155]	2419	RETURN
[1012]	2420
[2155]	2421	END FUNCTION gamm
	2422
[1012]	2423	END MODULE microphysics_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |