Home

Context Navigation

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

Last change on this file since 1320 was 1320, checked in by raasch, 10 years ago
ONLY-attribute added to USE-statements, kind-parameters added to all INTEGER and REAL declaration statements, kinds are defined in new module kinds, old module precision_kind is removed, revision history before 2012 removed, comment fields (!:) to be used for variable explanations added to all variable declaration statements
Property svn:keywords set to `Id`
File size: 37.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |