Home

Context Navigation

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

Last change on this file since 1354 was 1354, checked in by heinze, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 38.7 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	!
[1354]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: microphysics.f90 1354 2014-04-08 15:22:57Z heinze $
	27	!
[1354]	28	! 1353 2014-04-08 15:21:23Z heinze
	29	! REAL constants provided with KIND-attribute
	30	!
[1347]	31	! 1346 2014-03-27 13:18:20Z heinze
	32	! Bugfix: REAL constants provided with KIND-attribute especially in call of
	33	! intrinsic function like MAX, MIN, SIGN
	34	!
[1335]	35	! 1334 2014-03-25 12:21:40Z heinze
	36	! Bugfix: REAL constants provided with KIND-attribute
	37	!
[1323]	38	! 1322 2014-03-20 16:38:49Z raasch
	39	! REAL constants defined as wp-kind
	40	!
[1321]	41	! 1320 2014-03-20 08:40:49Z raasch
[1320]	42	! ONLY-attribute added to USE-statements,
	43	! kind-parameters added to all INTEGER and REAL declaration statements,
	44	! kinds are defined in new module kinds,
	45	! comment fields (!:) to be used for variable explanations added to
	46	! all variable declaration statements
[1000]	47	!
[1242]	48	! 1241 2013-10-30 11:36:58Z heinze
	49	! hyp and rho have to be calculated at each time step if data from external
	50	! file LSF_DATA are used
	51	!
[1116]	52	! 1115 2013-03-26 18:16:16Z hoffmann
	53	! microphyical tendencies are calculated in microphysics_control in an optimized
	54	! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
	55	!
[1107]	56	! 1106 2013-03-04 05:31:38Z raasch
	57	! small changes in code formatting
	58	!
[1093]	59	! 1092 2013-02-02 11:24:22Z raasch
	60	! unused variables removed
	61	! file put under GPL
	62	!
[1066]	63	! 1065 2012-11-22 17:42:36Z hoffmann
	64	! Sedimentation process implemented according to Stevens and Seifert (2008).
[1115]	65	! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
[1066]	66	! and Stevens, 2010).
	67	!
[1054]	68	! 1053 2012-11-13 17:11:03Z hoffmann
	69	! initial revision
[1000]	70	!
	71	! Description:
	72	! ------------
	73	! Calculate cloud microphysics according to the two moment bulk
	74	! scheme by Seifert and Beheng (2006).
	75	!------------------------------------------------------------------------------!
	76
	77	PRIVATE
[1115]	78	PUBLIC microphysics_control
[1000]	79
[1115]	80	INTERFACE microphysics_control
	81	MODULE PROCEDURE microphysics_control
	82	MODULE PROCEDURE microphysics_control_ij
	83	END INTERFACE microphysics_control
[1022]	84
[1115]	85	INTERFACE adjust_cloud
	86	MODULE PROCEDURE adjust_cloud
	87	MODULE PROCEDURE adjust_cloud_ij
	88	END INTERFACE adjust_cloud
	89
[1000]	90	INTERFACE autoconversion
	91	MODULE PROCEDURE autoconversion
	92	MODULE PROCEDURE autoconversion_ij
	93	END INTERFACE autoconversion
	94
	95	INTERFACE accretion
	96	MODULE PROCEDURE accretion
	97	MODULE PROCEDURE accretion_ij
	98	END INTERFACE accretion
[1005]	99
	100	INTERFACE selfcollection_breakup
	101	MODULE PROCEDURE selfcollection_breakup
	102	MODULE PROCEDURE selfcollection_breakup_ij
	103	END INTERFACE selfcollection_breakup
[1012]	104
	105	INTERFACE evaporation_rain
	106	MODULE PROCEDURE evaporation_rain
	107	MODULE PROCEDURE evaporation_rain_ij
	108	END INTERFACE evaporation_rain
	109
	110	INTERFACE sedimentation_cloud
	111	MODULE PROCEDURE sedimentation_cloud
	112	MODULE PROCEDURE sedimentation_cloud_ij
	113	END INTERFACE sedimentation_cloud
[1000]	114
[1012]	115	INTERFACE sedimentation_rain
	116	MODULE PROCEDURE sedimentation_rain
	117	MODULE PROCEDURE sedimentation_rain_ij
	118	END INTERFACE sedimentation_rain
	119
[1000]	120	CONTAINS
	121
	122
	123	!------------------------------------------------------------------------------!
	124	! Call for all grid points
	125	!------------------------------------------------------------------------------!
[1115]	126	SUBROUTINE microphysics_control
[1022]	127
	128	USE arrays_3d
[1241]	129	USE cloud_parameters
[1115]	130	USE control_parameters
[1241]	131	USE grid_variables
[1115]	132	USE indices
[1320]	133	USE kinds
[1115]	134	USE statistics
	135
	136	IMPLICIT NONE
	137
[1320]	138	INTEGER(iwp) :: i !:
	139	INTEGER(iwp) :: j !:
	140	INTEGER(iwp) :: k !:
[1115]	141
	142	DO i = nxl, nxr
	143	DO j = nys, nyn
	144	DO k = nzb_s_inner(j,i)+1, nzt
	145
	146	ENDDO
	147	ENDDO
	148	ENDDO
	149
	150	END SUBROUTINE microphysics_control
	151
	152	SUBROUTINE adjust_cloud
	153
	154	USE arrays_3d
[1022]	155	USE cloud_parameters
	156	USE indices
[1320]	157	USE kinds
[1022]	158
	159	IMPLICIT NONE
	160
[1320]	161	INTEGER(iwp) :: i !:
	162	INTEGER(iwp) :: j !:
	163	INTEGER(iwp) :: k !:
[1022]	164
	165
	166	DO i = nxl, nxr
	167	DO j = nys, nyn
[1115]	168	DO k = nzb_s_inner(j,i)+1, nzt
[1022]	169
	170	ENDDO
	171	ENDDO
	172	ENDDO
	173
[1115]	174	END SUBROUTINE adjust_cloud
[1022]	175
[1106]	176
[1000]	177	SUBROUTINE autoconversion
	178
	179	USE arrays_3d
	180	USE cloud_parameters
[1115]	181	USE control_parameters
	182	USE grid_variables
[1000]	183	USE indices
[1320]	184	USE kinds
[1000]	185
	186	IMPLICIT NONE
	187
[1320]	188	INTEGER(iwp) :: i !:
	189	INTEGER(iwp) :: j !:
	190	INTEGER(iwp) :: k !:
[1000]	191
	192	DO i = nxl, nxr
	193	DO j = nys, nyn
[1115]	194	DO k = nzb_s_inner(j,i)+1, nzt
[1000]	195
	196	ENDDO
	197	ENDDO
	198	ENDDO
	199
	200	END SUBROUTINE autoconversion
	201
[1106]	202
[1005]	203	SUBROUTINE accretion
[1000]	204
	205	USE arrays_3d
	206	USE cloud_parameters
[1115]	207	USE control_parameters
[1000]	208	USE indices
[1320]	209	USE kinds
[1005]	210
[1000]	211	IMPLICIT NONE
	212
[1320]	213	INTEGER(iwp) :: i !:
	214	INTEGER(iwp) :: j !:
	215	INTEGER(iwp) :: k !:
[1000]	216
[1005]	217	DO i = nxl, nxr
	218	DO j = nys, nyn
[1115]	219	DO k = nzb_s_inner(j,i)+1, nzt
[1000]	220
[1005]	221	ENDDO
	222	ENDDO
[1000]	223	ENDDO
	224
[1005]	225	END SUBROUTINE accretion
[1000]	226
[1106]	227
[1005]	228	SUBROUTINE selfcollection_breakup
[1000]	229
	230	USE arrays_3d
	231	USE cloud_parameters
[1115]	232	USE control_parameters
[1000]	233	USE indices
[1320]	234	USE kinds
[1000]	235
	236	IMPLICIT NONE
	237
[1320]	238	INTEGER(iwp) :: i !:
	239	INTEGER(iwp) :: j !:
	240	INTEGER(iwp) :: k !:
[1000]	241
	242
	243	DO i = nxl, nxr
	244	DO j = nys, nyn
[1115]	245	DO k = nzb_s_inner(j,i)+1, nzt
[1000]	246
	247	ENDDO
	248	ENDDO
	249	ENDDO
	250
[1005]	251	END SUBROUTINE selfcollection_breakup
[1000]	252
[1106]	253
[1012]	254	SUBROUTINE evaporation_rain
[1000]	255
[1012]	256	USE arrays_3d
	257	USE cloud_parameters
	258	USE constants
[1115]	259	USE control_parameters
[1012]	260	USE indices
[1320]	261	USE kinds
[1012]	262
	263	IMPLICIT NONE
	264
[1320]	265	INTEGER(iwp) :: i !:
	266	INTEGER(iwp) :: j !:
	267	INTEGER(iwp) :: k !:
[1012]	268
	269	DO i = nxl, nxr
	270	DO j = nys, nyn
[1115]	271	DO k = nzb_s_inner(j,i)+1, nzt
[1012]	272
	273	ENDDO
	274	ENDDO
	275	ENDDO
	276
	277	END SUBROUTINE evaporation_rain
	278
[1106]	279
[1012]	280	SUBROUTINE sedimentation_cloud
	281
	282	USE arrays_3d
	283	USE cloud_parameters
	284	USE constants
[1115]	285	USE control_parameters
[1012]	286	USE indices
[1320]	287	USE kinds
[1012]	288
	289	IMPLICIT NONE
	290
[1320]	291	INTEGER(iwp) :: i !:
	292	INTEGER(iwp) :: j !:
	293	INTEGER(iwp) :: k !:
[1012]	294
	295	DO i = nxl, nxr
	296	DO j = nys, nyn
[1115]	297	DO k = nzb_s_inner(j,i)+1, nzt
[1012]	298
	299	ENDDO
	300	ENDDO
	301	ENDDO
	302
	303	END SUBROUTINE sedimentation_cloud
	304
[1106]	305
[1012]	306	SUBROUTINE sedimentation_rain
	307
	308	USE arrays_3d
	309	USE cloud_parameters
	310	USE constants
[1115]	311	USE control_parameters
[1012]	312	USE indices
[1320]	313	USE kinds
[1115]	314	USE statistics
[1012]	315
	316	IMPLICIT NONE
	317
[1320]	318	INTEGER(iwp) :: i !:
	319	INTEGER(iwp) :: j !:
	320	INTEGER(iwp) :: k !:
[1012]	321
	322	DO i = nxl, nxr
	323	DO j = nys, nyn
[1115]	324	DO k = nzb_s_inner(j,i)+1, nzt
[1012]	325
	326	ENDDO
	327	ENDDO
	328	ENDDO
	329
	330	END SUBROUTINE sedimentation_rain
	331
	332
[1000]	333	!------------------------------------------------------------------------------!
	334	! Call for grid point i,j
	335	!------------------------------------------------------------------------------!
[1022]	336
[1115]	337	SUBROUTINE microphysics_control_ij( i, j )
	338
[1320]	339	USE arrays_3d, &
	340	ONLY: hyp, nc_1d, nr, nr_1d, pt, pt_init, pt_1d, q, q_1d, qc, &
	341	qc_1d, qr, qr_1d, tend_nr, tend_pt, tend_q, tend_qr, zu
[1115]	342
[1320]	343	USE cloud_parameters, &
	344	ONLY: cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt
	345
	346	USE control_parameters, &
	347	ONLY: drizzle, dt_3d, dt_micro, g, intermediate_timestep_count, &
	348	large_scale_forcing, lsf_surf, precipitation, pt_surface, &
	349	rho_surface,surface_pressure
	350
	351	USE indices, &
	352	ONLY: nzb, nzt
	353
	354	USE kinds
	355
	356	USE statistics, &
	357	ONLY: weight_pres
	358
[1022]	359	IMPLICIT NONE
	360
[1320]	361	INTEGER(iwp) :: i !:
	362	INTEGER(iwp) :: j !:
	363	INTEGER(iwp) :: k !:
[1115]	364
[1320]	365	REAL(wp) :: t_surface !:
	366
[1241]	367	IF ( large_scale_forcing .AND. lsf_surf ) THEN
	368	!
	369	!-- Calculate:
	370	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
	371	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
	372	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
[1353]	373	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
[1241]	374	DO k = nzb, nzt+1
[1353]	375	hyp(k) = surface_pressure * 100.0_wp * &
[1334]	376	( (t_surface - g/cp * zu(k)) / t_surface )**(1.0_wp/0.286_wp)
[1353]	377	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
	378	t_d_pt(k) = 1.0_wp / pt_d_t(k)
[1241]	379	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
	380	ENDDO
	381	!
	382	!-- Compute reference density
[1353]	383	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
[1241]	384	ENDIF
	385
	386
[1115]	387	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
	388	!
	389	!-- Adjust unrealistic values
	390	IF ( precipitation ) CALL adjust_cloud( i,j )
	391	!
	392	!-- Use 1-d arrays
	393	q_1d(:) = q(:,j,i)
	394	pt_1d(:) = pt(:,j,i)
	395	qc_1d(:) = qc(:,j,i)
	396	nc_1d(:) = nc_const
	397	IF ( precipitation ) THEN
	398	qr_1d(:) = qr(:,j,i)
	399	nr_1d(:) = nr(:,j,i)
	400	ENDIF
	401	!
	402	!-- Compute cloud physics
	403	IF ( precipitation ) THEN
	404	CALL autoconversion( i,j )
	405	CALL accretion( i,j )
	406	CALL selfcollection_breakup( i,j )
	407	CALL evaporation_rain( i,j )
	408	CALL sedimentation_rain( i,j )
	409	ENDIF
	410
	411	IF ( drizzle ) CALL sedimentation_cloud( i,j )
	412	!
	413	!-- Derive tendencies
	414	tend_q(:,j,i) = ( q_1d(:) - q(:,j,i) ) / dt_micro
	415	tend_pt(:,j,i) = ( pt_1d(:) - pt(:,j,i) ) / dt_micro
	416	IF ( precipitation ) THEN
	417	tend_qr(:,j,i) = ( qr_1d(:) - qr(:,j,i) ) / dt_micro
	418	tend_nr(:,j,i) = ( nr_1d(:) - nr(:,j,i) ) / dt_micro
	419	ENDIF
	420
	421	END SUBROUTINE microphysics_control_ij
	422
	423	SUBROUTINE adjust_cloud_ij( i, j )
	424
[1320]	425	USE arrays_3d, &
	426	ONLY: qr, nr
[1115]	427
[1320]	428	USE cloud_parameters, &
	429	ONLY: eps_sb, xrmin, xrmax, hyrho, k_cc, x0
	430
	431	USE indices, &
	432	ONLY: nzb, nzb_s_inner, nzt
	433
	434	USE kinds
	435
[1115]	436	IMPLICIT NONE
	437
[1320]	438	INTEGER(iwp) :: i !:
	439	INTEGER(iwp) :: j !:
	440	INTEGER(iwp) :: k !:
[1115]	441	!
	442	!-- Adjust number of raindrops to avoid nonlinear effects in
	443	!-- sedimentation and evaporation of rain drops due to too small or
	444	!-- too big weights of rain drops (Stevens and Seifert, 2008).
	445	!-- The same procedure is applied to cloud droplets if they are determined
	446	!-- prognostically.
	447	DO k = nzb_s_inner(j,i)+1, nzt
[1022]	448
[1065]	449	IF ( qr(k,j,i) <= eps_sb ) THEN
[1353]	450	qr(k,j,i) = 0.0_wp
	451	nr(k,j,i) = 0.0_wp
[1065]	452	ELSE
[1022]	453	!
[1048]	454	!-- Adjust number of raindrops to avoid nonlinear effects in
	455	!-- sedimentation and evaporation of rain drops due to too small or
[1065]	456	!-- too big weights of rain drops (Stevens and Seifert, 2008).
	457	IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN
	458	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
	459	ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN
	460	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
[1048]	461	ENDIF
[1115]	462
[1022]	463	ENDIF
[1115]	464
[1022]	465	ENDDO
	466
[1115]	467	END SUBROUTINE adjust_cloud_ij
[1022]	468
[1106]	469
[1005]	470	SUBROUTINE autoconversion_ij( i, j )
[1000]	471
[1320]	472	USE arrays_3d, &
	473	ONLY: diss, dzu, nc_1d, nr_1d, qc_1d, qr_1d
[1115]	474
[1320]	475	USE cloud_parameters, &
	476	ONLY: a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3, &
	477	c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air, x0
	478
	479	USE control_parameters, &
	480	ONLY: dt_micro, rho_surface, turbulence
	481
	482	USE grid_variables, &
	483	ONLY: dx, dy
	484
	485	USE indices, &
	486	ONLY: nzb, nzb_s_inner, nzt
	487
	488	USE kinds
	489
[1000]	490	IMPLICIT NONE
	491
[1320]	492	INTEGER(iwp) :: i !:
	493	INTEGER(iwp) :: j !:
	494	INTEGER(iwp) :: k !:
[1000]	495
[1320]	496	REAL(wp) :: alpha_cc !:
	497	REAL(wp) :: autocon !:
	498	REAL(wp) :: epsilon !:
	499	REAL(wp) :: k_au !:
	500	REAL(wp) :: l_mix !:
	501	REAL(wp) :: nu_c !:
	502	REAL(wp) :: phi_au !:
	503	REAL(wp) :: r_cc !:
	504	REAL(wp) :: rc !:
	505	REAL(wp) :: re_lambda !:
	506	REAL(wp) :: selfcoll !:
	507	REAL(wp) :: sigma_cc !:
	508	REAL(wp) :: tau_cloud !:
	509	REAL(wp) :: xc !:
[1106]	510
[1353]	511	k_au = k_cc / ( 20.0_wp * x0 )
[1005]	512
[1115]	513	DO k = nzb_s_inner(j,i)+1, nzt
[1000]	514
[1115]	515	IF ( qc_1d(k) > eps_sb ) THEN
[1012]	516	!
[1048]	517	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[1353]	518	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
	519	tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
[1012]	520	!
	521	!-- Universal function for autoconversion process
	522	!-- (Seifert and Beheng, 2006):
[1353]	523	phi_au = 600.0_wp * tau_cloud*0.68_wp ( 1.0_wp - tau_cloud0.68_wp )3
[1012]	524	!
	525	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
[1353]	526	!-- (Use constant nu_c = 1.0_wp instead?)
	527	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
[1012]	528	!
	529	!-- Mean weight of cloud droplets:
[1115]	530	xc = hyrho(k) * qc_1d(k) / nc_1d(k)
[1012]	531	!
[1065]	532	!-- Parameterized turbulence effects on autoconversion (Seifert,
	533	!-- Nuijens and Stevens, 2010)
	534	IF ( turbulence ) THEN
	535	!
	536	!-- Weight averaged radius of cloud droplets:
[1353]	537	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
[1065]	538
[1353]	539	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
	540	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
	541	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
[1065]	542	!
	543	!-- Mixing length (neglecting distance to ground and stratification)
[1334]	544	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
[1065]	545	!
	546	!-- Limit dissipation rate according to Seifert, Nuijens and
	547	!-- Stevens (2010)
[1334]	548	epsilon = MIN( 0.06_wp, diss(k,j,i) )
[1065]	549	!
	550	!-- Compute Taylor-microscale Reynolds number:
[1353]	551	re_lambda = 6.0_wp / 11.0_wp * ( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
	552	SQRT( 15.0_wp / kin_vis_air ) * epsilon**( 1.0_wp / 6.0_wp )
[1065]	553	!
	554	!-- The factor of 1.0E4 is needed to convert the dissipation rate
	555	!-- from m2 s-3 to cm2 s-3.
[1353]	556	k_au = k_au * ( 1.0_wp + &
	557	epsilon * 1.0E4_wp * ( re_lambda * 1.0E-3_wp )*0.25_wp &
	558	( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) / &
[1065]	559	sigma_cc )**2 ) + beta_cc ) )
	560	ENDIF
	561	!
[1012]	562	!-- Autoconversion rate (Seifert and Beheng, 2006):
[1353]	563	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
	564	( nu_c + 1.0_wp )*2.0_wp qc_1d(k)*2.0_wp xc*2.0_wp &
	565	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2.0_wp ) &
[1115]	566	rho_surface
	567	autocon = MIN( autocon, qc_1d(k) / dt_micro )
[1106]	568
[1115]	569	qr_1d(k) = qr_1d(k) + autocon * dt_micro
	570	qc_1d(k) = qc_1d(k) - autocon * dt_micro
	571	nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
	572
[1005]	573	ENDIF
[1000]	574
	575	ENDDO
	576
[1005]	577	END SUBROUTINE autoconversion_ij
	578
[1106]	579
[1005]	580	SUBROUTINE accretion_ij( i, j )
	581
[1320]	582	USE arrays_3d, &
	583	ONLY: diss, qc_1d, qr_1d
[1115]	584
[1320]	585	USE cloud_parameters, &
	586	ONLY: eps_sb, hyrho, k_cr0
	587
	588	USE control_parameters, &
	589	ONLY: dt_micro, rho_surface, turbulence
	590
	591	USE indices, &
	592	ONLY: nzb, nzb_s_inner, nzt
	593
	594	USE kinds
	595
[1005]	596	IMPLICIT NONE
	597
[1320]	598	INTEGER(iwp) :: i !:
	599	INTEGER(iwp) :: j !:
	600	INTEGER(iwp) :: k !:
[1005]	601
[1320]	602	REAL(wp) :: accr !:
	603	REAL(wp) :: k_cr !:
	604	REAL(wp) :: phi_ac !:
	605	REAL(wp) :: tau_cloud !:
	606	REAL(wp) :: xc !:
	607
[1115]	608	DO k = nzb_s_inner(j,i)+1, nzt
	609	IF ( ( qc_1d(k) > eps_sb ) .AND. ( qr_1d(k) > eps_sb ) ) THEN
[1012]	610	!
[1048]	611	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[1353]	612	tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
[1012]	613	!
	614	!-- Universal function for accretion process
[1048]	615	!-- (Seifert and Beheng, 2001):
[1353]	616	phi_ac = tau_cloud / ( tau_cloud + 5.0E-5_wp )
[1334]	617	phi_ac = ( phi_ac2.0_wp )2.0_wp
[1012]	618	!
[1065]	619	!-- Parameterized turbulence effects on autoconversion (Seifert,
	620	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
	621	!-- convert the dissipation (diss) from m2 s-3 to cm2 s-3.
	622	IF ( turbulence ) THEN
[1353]	623	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
	624	MIN( 600.0_wp, diss(k,j,i) * 1.0E4_wp )**0.25_wp )
[1065]	625	ELSE
	626	k_cr = k_cr0
	627	ENDIF
	628	!
[1012]	629	!-- Accretion rate (Seifert and Beheng, 2006):
[1115]	630	accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * &
[1065]	631	SQRT( rho_surface * hyrho(k) )
[1115]	632	accr = MIN( accr, qc_1d(k) / dt_micro )
[1106]	633
[1115]	634	qr_1d(k) = qr_1d(k) + accr * dt_micro
	635	qc_1d(k) = qc_1d(k) - accr * dt_micro
	636
[1005]	637	ENDIF
[1106]	638
[1005]	639	ENDDO
	640
[1000]	641	END SUBROUTINE accretion_ij
	642
[1005]	643
	644	SUBROUTINE selfcollection_breakup_ij( i, j )
	645
[1320]	646	USE arrays_3d, &
	647	ONLY: nr_1d, qr_1d
	648
	649	USE cloud_parameters, &
	650	ONLY: dpirho_l, eps_sb, hyrho, k_br, k_rr
	651
	652	USE control_parameters, &
	653	ONLY: dt_micro, rho_surface
	654
	655	USE indices, &
	656	ONLY: nzb, nzb_s_inner, nzt
	657
	658	USE kinds
[1005]	659
	660	IMPLICIT NONE
	661
[1320]	662	INTEGER(iwp) :: i !:
	663	INTEGER(iwp) :: j !:
	664	INTEGER(iwp) :: k !:
[1005]	665
[1320]	666	REAL(wp) :: breakup !:
	667	REAL(wp) :: dr !:
	668	REAL(wp) :: phi_br !:
	669	REAL(wp) :: selfcoll !:
	670
[1115]	671	DO k = nzb_s_inner(j,i)+1, nzt
	672	IF ( qr_1d(k) > eps_sb ) THEN
[1012]	673	!
[1115]	674	!-- Selfcollection rate (Seifert and Beheng, 2001):
	675	selfcoll = k_rr * nr_1d(k) * qr_1d(k) * &
[1005]	676	SQRT( hyrho(k) * rho_surface )
[1012]	677	!
[1115]	678	!-- Weight averaged diameter of rain drops:
[1334]	679	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	680	!
[1048]	681	!-- Collisional breakup rate (Seifert, 2008):
[1353]	682	IF ( dr >= 0.3E-3_wp ) THEN
	683	phi_br = k_br * ( dr - 1.1E-3_wp )
	684	breakup = selfcoll * ( phi_br + 1.0_wp )
[1005]	685	ELSE
[1353]	686	breakup = 0.0_wp
[1005]	687	ENDIF
[1048]	688
[1115]	689	selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
	690	nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
[1106]	691
[1005]	692	ENDIF
	693	ENDDO
	694
	695	END SUBROUTINE selfcollection_breakup_ij
	696
[1106]	697
[1012]	698	SUBROUTINE evaporation_rain_ij( i, j )
[1022]	699	!
	700	!-- Evaporation of precipitable water. Condensation is neglected for
	701	!-- precipitable water.
[1012]	702
[1320]	703	USE arrays_3d, &
	704	ONLY: hyp, nr_1d, pt_1d, q_1d, qc_1d, qr_1d
[1048]	705
[1320]	706	USE cloud_parameters, &
	707	ONLY: a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,&
	708	dpirho_l, eps_sb, hyrho, kin_vis_air, k_st, l_d_cp, l_d_r, &
	709	l_v, rho_l, r_v, schmidt_p_1d3, thermal_conductivity_l, &
	710	t_d_pt, ventilation_effect
	711
	712	USE constants, &
	713	ONLY: pi
	714
	715	USE control_parameters, &
	716	ONLY: dt_micro
	717
	718	USE indices, &
	719	ONLY: nzb, nzb_s_inner, nzt
	720
	721	USE kinds
	722
[1012]	723	IMPLICIT NONE
	724
[1320]	725	INTEGER(iwp) :: i !:
	726	INTEGER(iwp) :: j !:
	727	INTEGER(iwp) :: k !:
[1012]	728
[1320]	729	REAL(wp) :: alpha !:
	730	REAL(wp) :: dr !:
	731	REAL(wp) :: e_s !:
	732	REAL(wp) :: evap !:
	733	REAL(wp) :: evap_nr !:
	734	REAL(wp) :: f_vent !:
	735	REAL(wp) :: g_evap !:
	736	REAL(wp) :: lambda_r !:
	737	REAL(wp) :: mu_r !:
	738	REAL(wp) :: mu_r_2 !:
	739	REAL(wp) :: mu_r_5d2 !:
	740	REAL(wp) :: nr_0 !:
	741	REAL(wp) :: q_s !:
	742	REAL(wp) :: sat !:
	743	REAL(wp) :: t_l !:
	744	REAL(wp) :: temp !:
	745	REAL(wp) :: xr !:
	746
[1115]	747	DO k = nzb_s_inner(j,i)+1, nzt
	748	IF ( qr_1d(k) > eps_sb ) THEN
[1012]	749	!
	750	!-- Actual liquid water temperature:
[1115]	751	t_l = t_d_pt(k) * pt_1d(k)
[1012]	752	!
	753	!-- Saturation vapor pressure at t_l:
[1353]	754	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / ( t_l - 35.86_wp ) )
[1012]	755	!
	756	!-- Computation of saturation humidity:
[1353]	757	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
	758	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
	759	q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
[1012]	760	!
[1106]	761	!-- Supersaturation:
[1353]	762	sat = MIN( 0.0_wp, ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp )
[1012]	763	!
	764	!-- Actual temperature:
[1115]	765	temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
	766
[1353]	767	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v / &
[1115]	768	( thermal_conductivity_l * temp ) + r_v * temp / &
	769	( diff_coeff_l * e_s ) )
[1012]	770	!
[1115]	771	!-- Mean weight of rain drops
	772	xr = hyrho(k) * qr_1d(k) / nr_1d(k)
[1012]	773	!
[1115]	774	!-- Weight averaged diameter of rain drops:
[1334]	775	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	776	!
[1049]	777	!-- Compute ventilation factor and intercept parameter
	778	!-- (Seifert and Beheng, 2006; Seifert, 2008):
[1048]	779	IF ( ventilation_effect ) THEN
[1115]	780	!
	781	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	782	!-- Stevens and Seifert, 2008):
[1353]	783	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
[1115]	784	!
	785	!-- Slope parameter of gamma distribution (Seifert, 2008):
[1353]	786	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	787	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
[1115]	788
[1353]	789	mu_r_2 = mu_r + 2.0_wp
	790	mu_r_5d2 = mu_r + 2.5_wp
[1048]	791	f_vent = a_vent * gamm( mu_r_2 ) * &
[1115]	792	lambda_r**( -mu_r_2 ) + &
[1048]	793	b_vent * schmidt_p_1d3 * &
	794	SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * &
[1115]	795	lambda_r*( -mu_r_5d2 ) &
[1353]	796	( 1.0_wp - 0.5_wp * ( b_term / a_term ) * &
[1115]	797	( lambda_r / &
	798	( c_term + lambda_r ) )**mu_r_5d2 - &
[1353]	799	0.125_wp * ( b_term / a_term )*2.0_wp &
[1115]	800	( lambda_r / &
[1353]	801	( 2.0_wp * c_term + lambda_r ) )**mu_r_5d2 - &
	802	0.0625_wp * ( b_term / a_term )*3.0_wp &
[1115]	803	( lambda_r / &
[1353]	804	( 3.0_wp * c_term + lambda_r ) )**mu_r_5d2 - &
	805	0.0390625_wp * ( b_term / a_term )*4.0_wp &
[1115]	806	( lambda_r / &
[1353]	807	( 4.0_wp * c_term + lambda_r ) )**mu_r_5d2 )
	808	nr_0 = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) / &
	809	gamm( mu_r + 1.0_wp )
[1048]	810	ELSE
[1353]	811	f_vent = 1.0_wp
[1115]	812	nr_0 = nr_1d(k) * dr
[1048]	813	ENDIF
[1012]	814	!
[1048]	815	!-- Evaporation rate of rain water content (Seifert and Beheng, 2006):
[1353]	816	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / &
[1048]	817	hyrho(k)
[1106]	818
[1115]	819	evap = MAX( evap, -qr_1d(k) / dt_micro )
	820	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
	821	-nr_1d(k) / dt_micro )
	822
	823	qr_1d(k) = qr_1d(k) + evap * dt_micro
	824	nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
[1012]	825	ENDIF
[1106]	826
[1012]	827	ENDDO
	828
	829	END SUBROUTINE evaporation_rain_ij
	830
[1106]	831
[1012]	832	SUBROUTINE sedimentation_cloud_ij( i, j )
	833
[1320]	834	USE arrays_3d, &
	835	ONLY: ddzu, dzu, nc_1d, pt_1d, q_1d, qc_1d
	836
	837	USE cloud_parameters, &
	838	ONLY: eps_sb, hyrho, k_st, l_d_cp, prr, pt_d_t, rho_l, sigma_gc
	839
	840	USE constants, &
	841	ONLY: pi
	842
	843	USE control_parameters, &
	844	ONLY: dt_do2d_xy, dt_micro, intermediate_timestep_count
	845
	846	USE indices, &
	847	ONLY: nzb, nzb_s_inner, nzt
	848
	849	USE kinds
[1012]	850
	851	IMPLICIT NONE
	852
[1320]	853	INTEGER(iwp) :: i !:
	854	INTEGER(iwp) :: j !:
	855	INTEGER(iwp) :: k !:
[1106]	856
[1320]	857	REAL(wp) :: sed_qc_const !:
[1115]	858
[1320]	859
	860	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc
	861
[1012]	862	!
	863	!-- Sedimentation of cloud droplets (Heus et al., 2010):
[1353]	864	sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l ))*( 2.0_wp / 3.0_wp ) &
	865	EXP( 5.0_wp * LOG( sigma_gc )**2 )
[1012]	866
[1353]	867	sed_qc(nzt+1) = 0.0_wp
[1012]	868
[1115]	869	DO k = nzt, nzb_s_inner(j,i)+1, -1
	870	IF ( qc_1d(k) > eps_sb ) THEN
[1334]	871	sed_qc(k) = sed_qc_const * nc_1d(k)*( -2.0_wp / 3.0_wp ) &
	872	( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )
[1115]	873	ELSE
[1353]	874	sed_qc(k) = 0.0_wp
[1012]	875	ENDIF
[1115]	876
	877	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) / &
	878	dt_micro + sed_qc(k+1) )
	879
	880	q_1d(k) = q_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
	881	hyrho(k) * dt_micro
	882	qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
	883	hyrho(k) * dt_micro
	884	pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
	885	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
	886
[1012]	887	ENDDO
	888
	889	END SUBROUTINE sedimentation_cloud_ij
	890
[1106]	891
[1012]	892	SUBROUTINE sedimentation_rain_ij( i, j )
	893
[1320]	894	USE arrays_3d, &
	895	ONLY: ddzu, dzu, nr_1d, pt_1d, q_1d, qr_1d
	896
	897	USE cloud_parameters, &
	898	ONLY: a_term, b_term, c_term, cof, dpirho_l, eps_sb, hyrho, &
	899	limiter_sedimentation, l_d_cp, precipitation_amount, prr, &
	900	pt_d_t, stp
	901
	902	USE control_parameters, &
	903	ONLY: dt_do2d_xy, dt_micro, dt_3d, intermediate_timestep_count, &
	904	intermediate_timestep_count_max, &
	905	precipitation_amount_interval, time_do2d_xy
	906
	907	USE indices, &
	908	ONLY: nzb, nzb_s_inner, nzt
	909
	910	USE kinds
	911
	912	USE statistics, &
	913	ONLY: weight_substep
[1012]	914
	915	IMPLICIT NONE
	916
[1320]	917	INTEGER(iwp) :: i !:
	918	INTEGER(iwp) :: j !:
	919	INTEGER(iwp) :: k !:
	920	INTEGER(iwp) :: k_run !:
[1012]	921
[1320]	922	REAL(wp) :: c_run !:
	923	REAL(wp) :: d_max !:
	924	REAL(wp) :: d_mean !:
	925	REAL(wp) :: d_min !:
	926	REAL(wp) :: dr !:
	927	REAL(wp) :: dt_sedi !:
	928	REAL(wp) :: flux !:
	929	REAL(wp) :: lambda_r !:
	930	REAL(wp) :: mu_r !:
	931	REAL(wp) :: z_run !:
	932
	933	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !:
	934	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !:
	935	REAL(wp), DIMENSION(nzb:nzt+1) :: d_nr !:
	936	REAL(wp), DIMENSION(nzb:nzt+1) :: d_qr !:
	937	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !:
	938	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !:
	939	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !:
	940	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !:
	941	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !:
	942	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !:
	943
	944
[1065]	945	!
	946	!-- Computation of sedimentation flux. Implementation according to Stevens
	947	!-- and Seifert (2008).
[1353]	948	IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp
[1012]	949	!
[1065]	950	!-- Compute velocities
	951	DO k = nzb_s_inner(j,i)+1, nzt
[1115]	952	IF ( qr_1d(k) > eps_sb ) THEN
	953	!
	954	!-- Weight averaged diameter of rain drops:
[1334]	955	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	956	!
	957	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	958	!-- Stevens and Seifert, 2008):
[1353]	959	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
[1115]	960	!
	961	!-- Slope parameter of gamma distribution (Seifert, 2008):
[1353]	962	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	963	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
[1115]	964
[1353]	965	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, a_term - b_term * ( 1.0_wp + &
	966	c_term / lambda_r )*( -1.0_wp ( mu_r + 1.0_wp ) ) ) )
	967	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, a_term - b_term * ( 1.0_wp + &
	968	c_term / lambda_r )*( -1.0_wp ( mu_r + 4.0_wp ) ) ) )
[1065]	969	ELSE
[1353]	970	w_nr(k) = 0.0_wp
	971	w_qr(k) = 0.0_wp
[1065]	972	ENDIF
	973	ENDDO
[1048]	974	!
[1065]	975	!-- Adjust boundary values
[1115]	976	w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
	977	w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
[1353]	978	w_nr(nzt+1) = 0.0_wp
	979	w_qr(nzt+1) = 0.0_wp
[1065]	980	!
	981	!-- Compute Courant number
[1115]	982	DO k = nzb_s_inner(j,i)+1, nzt
[1353]	983	c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) * &
[1115]	984	dt_micro * ddzu(k)
[1353]	985	c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * &
[1115]	986	dt_micro * ddzu(k)
	987	ENDDO
[1065]	988	!
	989	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
	990	IF ( limiter_sedimentation ) THEN
	991
[1115]	992	DO k = nzb_s_inner(j,i)+1, nzt
[1353]	993	d_mean = 0.5_wp * ( qr_1d(k+1) + qr_1d(k-1) )
[1115]	994	d_min = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
	995	d_max = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
[1065]	996
[1353]	997	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, 2.0_wp * d_max, &
[1065]	998	ABS( d_mean ) )
	999
[1353]	1000	d_mean = 0.5_wp * ( nr_1d(k+1) + nr_1d(k-1) )
[1115]	1001	d_min = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
	1002	d_max = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
[1065]	1003
[1353]	1004	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, 2.0_wp * d_max, &
[1065]	1005	ABS( d_mean ) )
[1022]	1006	ENDDO
[1048]	1007
[1065]	1008	ELSE
[1106]	1009
[1353]	1010	nr_slope = 0.0_wp
	1011	qr_slope = 0.0_wp
[1106]	1012
[1065]	1013	ENDIF
[1115]	1014
[1353]	1015	sed_nr(nzt+1) = 0.0_wp
	1016	sed_qr(nzt+1) = 0.0_wp
[1065]	1017	!
	1018	!-- Compute sedimentation flux
[1115]	1019	DO k = nzt, nzb_s_inner(j,i)+1, -1
[1065]	1020	!
	1021	!-- Sum up all rain drop number densities which contribute to the flux
	1022	!-- through k-1/2
[1353]	1023	flux = 0.0_wp
	1024	z_run = 0.0_wp ! height above z(k)
[1065]	1025	k_run = k
[1346]	1026	c_run = MIN( 1.0_wp, c_nr(k) )
[1353]	1027	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
[1065]	1028	flux = flux + hyrho(k_run) * &
[1353]	1029	( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) * &
	1030	0.5_wp ) * c_run * dzu(k_run)
[1065]	1031	z_run = z_run + dzu(k_run)
	1032	k_run = k_run + 1
[1346]	1033	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
[1022]	1034	ENDDO
	1035	!
[1065]	1036	!-- It is not allowed to sediment more rain drop number density than
	1037	!-- available
	1038	flux = MIN( flux, &
[1115]	1039	hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
[1065]	1040
[1115]	1041	sed_nr(k) = flux / dt_micro
	1042	nr_1d(k) = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / &
	1043	hyrho(k) * dt_micro
[1065]	1044	!
	1045	!-- Sum up all rain water content which contributes to the flux
	1046	!-- through k-1/2
[1353]	1047	flux = 0.0_wp
	1048	z_run = 0.0_wp ! height above z(k)
[1065]	1049	k_run = k
[1346]	1050	c_run = MIN( 1.0_wp, c_qr(k) )
[1106]	1051
[1353]	1052	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt-1 )
[1106]	1053
[1065]	1054	flux = flux + hyrho(k_run) * &
[1353]	1055	( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) * &
	1056	0.5_wp ) * c_run * dzu(k_run)
[1065]	1057	z_run = z_run + dzu(k_run)
	1058	k_run = k_run + 1
[1346]	1059	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
[1106]	1060
[1065]	1061	ENDDO
	1062	!
	1063	!-- It is not allowed to sediment more rain water content than available
	1064	flux = MIN( flux, &
[1115]	1065	hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
[1065]	1066
[1115]	1067	sed_qr(k) = flux / dt_micro
	1068
	1069	qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
	1070	hyrho(k) * dt_micro
	1071	q_1d(k) = q_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
	1072	hyrho(k) * dt_micro
	1073	pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
	1074	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
[1065]	1075	!
	1076	!-- Compute the rain rate
	1077	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * &
[1115]	1078	weight_substep(intermediate_timestep_count)
[1065]	1079	ENDDO
[1115]	1080
[1065]	1081	!
[1048]	1082	!-- Precipitation amount
	1083	IF ( intermediate_timestep_count == intermediate_timestep_count_max &
	1084	.AND. ( dt_do2d_xy - time_do2d_xy ) < &
	1085	precipitation_amount_interval ) THEN
[1012]	1086
[1048]	1087	precipitation_amount(j,i) = precipitation_amount(j,i) + &
[1115]	1088	prr(nzb_s_inner(j,i)+1,j,i) * &
	1089	hyrho(nzb_s_inner(j,i)+1) * dt_3d
[1048]	1090	ENDIF
	1091
[1012]	1092	END SUBROUTINE sedimentation_rain_ij
	1093
[1106]	1094
[1012]	1095	!
	1096	!-- This function computes the gamma function (Press et al., 1992).
	1097	!-- The gamma function is needed for the calculation of the evaporation
	1098	!-- of rain drops.
	1099	FUNCTION gamm( xx )
[1048]	1100
[1320]	1101	USE cloud_parameters, &
	1102	ONLY: cof, stp
	1103
	1104	USE kinds
	1105
[1012]	1106	IMPLICIT NONE
[1106]	1107
[1320]	1108	INTEGER(iwp) :: j !:
	1109
	1110	REAL(wp) :: gamm !:
	1111	REAL(wp) :: ser !:
	1112	REAL(wp) :: tmp !:
	1113	REAL(wp) :: x_gamm !:
	1114	REAL(wp) :: xx !:
	1115	REAL(wp) :: y_gamm !:
	1116
[1012]	1117	x_gamm = xx
	1118	y_gamm = x_gamm
[1353]	1119	tmp = x_gamm + 5.5_wp
	1120	tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
[1334]	1121	ser = 1.000000000190015_wp
[1106]	1122
	1123	DO j = 1, 6
[1353]	1124	y_gamm = y_gamm + 1.0_wp
[1012]	1125	ser = ser + cof( j ) / y_gamm
[1106]	1126	ENDDO
	1127
[1012]	1128	!
	1129	!-- Until this point the algorithm computes the logarithm of the gamma
	1130	!-- function. Hence, the exponential function is used.
	1131	! gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
	1132	gamm = EXP( tmp ) * stp * ser / x_gamm
[1106]	1133
[1012]	1134	RETURN
	1135
	1136	END FUNCTION gamm
	1137
	1138	END MODULE microphysics_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |