Home

Context Navigation

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

Last change on this file since 1859 was 1851, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 86.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |