Home

Context Navigation

source: palm/trunk/SOURCE/microphysics.f90 @ 1532

Last change on this file since 1532 was 1362, checked in by hoffmann, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 72.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |