Home

Context Navigation

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

Last change on this file since 2230 was 2156, checked in by hoffmann, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 86.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |