Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |