Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |