Home

Context Navigation

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

Last change on this file since 1691 was 1691, checked in by maronga, 9 years ago
various bugfixes and modifications of the atmosphere-land-surface-radiation interaction. Completely re-written routine to calculate surface fluxes (surface_layer_fluxes.f90) that replaces prandtl_fluxes. Minor formatting corrections and renamings
Property svn:keywords set to `Id`
File size: 79.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |