Home

Context Navigation

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

Last change on this file since 1346 was 1346, checked in by heinze, 10 years ago
Bugfix: REAL constants provided with KIND-attribute especially in call of intrinsic function like MAX, MIN, SIGN
Property svn:keywords set to `Id`
File size: 38.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |