Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |