Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |