Home

Context Navigation

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

Last change on this file since 1060 was 1054, checked in by hoffmann, 12 years ago
last commit documented
File size: 19.7 KB

Rev	Line
[1000]	1	MODULE microphysics_mod
	2
	3	!------------------------------------------------------------------------------!
	4	! Current revisions:
	5	! -----------------
	6	!
	7	! Former revisions:
	8	! -----------------
[1052]	9	! $Id$
[1054]	10	!
	11	! 1053 2012-11-13 17:11:03Z hoffmann
	12	! initial revision
[1000]	13	!
	14	! Description:
	15	! ------------
	16	! Calculate cloud microphysics according to the two moment bulk
	17	! scheme by Seifert and Beheng (2006).
	18	!------------------------------------------------------------------------------!
	19
	20	PRIVATE
[1022]	21	PUBLIC dsd_properties, autoconversion, accretion, selfcollection_breakup, &
[1012]	22	evaporation_rain, sedimentation_cloud, sedimentation_rain
[1000]	23
[1022]	24	INTERFACE dsd_properties
	25	MODULE PROCEDURE dsd_properties
	26	MODULE PROCEDURE dsd_properties_ij
	27	END INTERFACE dsd_properties
	28
[1000]	29	INTERFACE autoconversion
	30	MODULE PROCEDURE autoconversion
	31	MODULE PROCEDURE autoconversion_ij
	32	END INTERFACE autoconversion
	33
	34	INTERFACE accretion
	35	MODULE PROCEDURE accretion
	36	MODULE PROCEDURE accretion_ij
	37	END INTERFACE accretion
[1005]	38
	39	INTERFACE selfcollection_breakup
	40	MODULE PROCEDURE selfcollection_breakup
	41	MODULE PROCEDURE selfcollection_breakup_ij
	42	END INTERFACE selfcollection_breakup
[1012]	43
	44	INTERFACE evaporation_rain
	45	MODULE PROCEDURE evaporation_rain
	46	MODULE PROCEDURE evaporation_rain_ij
	47	END INTERFACE evaporation_rain
	48
	49	INTERFACE sedimentation_cloud
	50	MODULE PROCEDURE sedimentation_cloud
	51	MODULE PROCEDURE sedimentation_cloud_ij
	52	END INTERFACE sedimentation_cloud
[1000]	53
[1012]	54	INTERFACE sedimentation_rain
	55	MODULE PROCEDURE sedimentation_rain
	56	MODULE PROCEDURE sedimentation_rain_ij
	57	END INTERFACE sedimentation_rain
	58
[1000]	59	CONTAINS
	60
	61
	62	!------------------------------------------------------------------------------!
	63	! Call for all grid points
	64	!------------------------------------------------------------------------------!
[1022]	65	SUBROUTINE dsd_properties
	66
	67	USE arrays_3d
	68	USE cloud_parameters
	69	USE constants
	70	USE indices
	71
	72	IMPLICIT NONE
	73
	74	INTEGER :: i, j, k
	75	REAL :: dqdt_precip
	76
	77
	78	DO i = nxl, nxr
	79	DO j = nys, nyn
	80	DO k = nzb_2d(j,i)+1, nzt
	81
	82	ENDDO
	83	ENDDO
	84	ENDDO
	85
	86	END SUBROUTINE dsd_properties
	87
[1000]	88	SUBROUTINE autoconversion
	89
	90	USE arrays_3d
	91	USE cloud_parameters
	92	USE constants
	93	USE indices
	94
	95	IMPLICIT NONE
	96
	97	INTEGER :: i, j, k
	98	REAL :: dqdt_precip
	99
	100
	101	DO i = nxl, nxr
	102	DO j = nys, nyn
	103	DO k = nzb_2d(j,i)+1, nzt
	104
	105	ENDDO
	106	ENDDO
	107	ENDDO
	108
	109	END SUBROUTINE autoconversion
	110
[1005]	111	SUBROUTINE accretion
[1000]	112
	113	USE arrays_3d
	114	USE cloud_parameters
	115	USE constants
	116	USE indices
[1005]	117
[1000]	118	IMPLICIT NONE
	119
	120	INTEGER :: i, j, k
[1005]	121	REAL :: dqdt_precip
[1000]	122
[1005]	123
	124	DO i = nxl, nxr
	125	DO j = nys, nyn
	126	DO k = nzb_2d(j,i)+1, nzt
[1000]	127
[1005]	128	ENDDO
	129	ENDDO
[1000]	130	ENDDO
	131
[1005]	132	END SUBROUTINE accretion
[1000]	133
[1005]	134	SUBROUTINE selfcollection_breakup
[1000]	135
	136	USE arrays_3d
	137	USE cloud_parameters
	138	USE constants
	139	USE indices
	140
	141	IMPLICIT NONE
	142
	143	INTEGER :: i, j, k
	144	REAL :: dqdt_precip
	145
	146
	147	DO i = nxl, nxr
	148	DO j = nys, nyn
	149	DO k = nzb_2d(j,i)+1, nzt
	150
	151	ENDDO
	152	ENDDO
	153	ENDDO
	154
[1005]	155	END SUBROUTINE selfcollection_breakup
[1000]	156
[1012]	157	SUBROUTINE evaporation_rain
[1000]	158
[1012]	159	USE arrays_3d
	160	USE cloud_parameters
	161	USE constants
	162	USE indices
	163
	164	IMPLICIT NONE
	165
	166	INTEGER :: i, j, k
	167	REAL :: dqdt_precip
	168
	169
	170	DO i = nxl, nxr
	171	DO j = nys, nyn
	172	DO k = nzb_2d(j,i)+1, nzt
	173
	174	ENDDO
	175	ENDDO
	176	ENDDO
	177
	178	END SUBROUTINE evaporation_rain
	179
	180	SUBROUTINE sedimentation_cloud
	181
	182	USE arrays_3d
	183	USE cloud_parameters
	184	USE constants
	185	USE indices
	186
	187	IMPLICIT NONE
	188
	189	INTEGER :: i, j, k
	190	REAL :: dqdt_precip
	191
	192
	193	DO i = nxl, nxr
	194	DO j = nys, nyn
	195	DO k = nzb_2d(j,i)+1, nzt
	196
	197	ENDDO
	198	ENDDO
	199	ENDDO
	200
	201	END SUBROUTINE sedimentation_cloud
	202
	203	SUBROUTINE sedimentation_rain
	204
	205	USE arrays_3d
	206	USE cloud_parameters
	207	USE constants
	208	USE indices
	209
	210	IMPLICIT NONE
	211
	212	INTEGER :: i, j, k
	213	REAL :: dqdt_precip
	214
	215
	216	DO i = nxl, nxr
	217	DO j = nys, nyn
	218	DO k = nzb_2d(j,i)+1, nzt
	219
	220	ENDDO
	221	ENDDO
	222	ENDDO
	223
	224
	225	END SUBROUTINE sedimentation_rain
	226
	227
[1000]	228	!------------------------------------------------------------------------------!
	229	! Call for grid point i,j
	230	!------------------------------------------------------------------------------!
[1022]	231	SUBROUTINE dsd_properties_ij( i, j )
	232
	233	USE arrays_3d
	234	USE cloud_parameters
	235	USE constants
	236	USE indices
	237	USE control_parameters
[1048]	238	USE user
[1022]	239
	240	IMPLICIT NONE
	241
	242	INTEGER :: i, j, k
[1048]	243	REAL :: dr_min = 2.0E-6, dr_max = 1.0E-3
[1022]	244
	245	DO k = nzb_2d(j,i)+1, nzt
[1048]	246	IF ( ( qr(k,j,i) > eps_sb ) .AND. ( nr(k,j,i) > eps_sb ) ) THEN
[1022]	247	!
[1048]	248	!-- Weight averaged diameter of rain drops:
[1049]	249	dr(k) = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * &
	250	dpirho_l )**( 1.0 / 3.0 )
[1022]	251	!
[1048]	252	!-- Adjust number of raindrops to avoid nonlinear effects in
	253	!-- sedimentation and evaporation of rain drops due to too small or
	254	!-- big diameters of rain drops (Stevens and Seifert, 2008).
	255	IF ( dr(k) < dr_min ) THEN
	256	nr(k,j,i) = qr(k,j,i) * hyrho(k) / dr_min*3 dpirho_l
	257	dr(k) = dr_min
	258	ELSEIF ( dr(k) > dr_max ) THEN
	259	nr(k,j,i) = qr(k,j,i) * hyrho(k) / dr_max*3 dpirho_l
	260	dr(k) = dr_max
	261	ENDIF
[1022]	262	!
[1049]	263	!-- Mean weight of rain drops (Seifert and Beheng, 2006):
	264	xr(k) = MIN( MAX( hyrho(k) * qr(k,j,i) / nr(k,j,i), xrmin ), xrmax)
[1048]	265	!
[1049]	266	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	267	!-- Stevens and Seifert, 2008):
	268	IF ( .NOT. mu_constant ) THEN
	269	mu_r(k) = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr(k) - 1.4E-3 ) ) )
	270	ELSE
	271	mu_r(k) = mu_constant_value
	272	ENDIF
	273	!
[1022]	274	!-- Slope parameter of gamma distribution (Seifert, 2008):
[1048]	275	lambda_r(k) = ( ( mu_r(k) + 3.0 ) * ( mu_r(k) + 2.0 ) * &
	276	( mu_r(k) + 1.0 ) )**( 1.0 / 3.0 ) / dr(k)
[1022]	277	ENDIF
	278	ENDDO
	279
	280	END SUBROUTINE dsd_properties_ij
	281
[1005]	282	SUBROUTINE autoconversion_ij( i, j )
[1000]	283
	284	USE arrays_3d
	285	USE cloud_parameters
	286	USE constants
	287	USE indices
[1005]	288	USE control_parameters
[1048]	289	USE statistics
[1000]	290
	291	IMPLICIT NONE
	292
	293	INTEGER :: i, j, k
[1048]	294	REAL :: k_au, autocon, phi_au, tau_cloud, xc, nu_c
[1000]	295
[1005]	296	k_au = k_cc / ( 20.0 * x0 )
	297
[1000]	298	DO k = nzb_2d(j,i)+1, nzt
	299
[1048]	300	IF ( ql(k,j,i) > 0.0 ) THEN
[1012]	301	!
[1048]	302	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[1012]	303	!-- (1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) ))
[1048]	304	tau_cloud = 1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) + 1.0E-20 )
[1012]	305	!
	306	!-- Universal function for autoconversion process
	307	!-- (Seifert and Beheng, 2006):
[1048]	308	phi_au = 600.0 * tau_cloud*0.68 ( 1.0 - tau_cloud0.68 )3
[1012]	309	!
	310	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
	311	!-- (Use constant nu_c = 1.0 instead?)
	312	nu_c = 1580.0 * hyrho(k) * ql(k,j,i) - 0.28
	313	!
	314	!-- Mean weight of cloud droplets:
[1005]	315	xc = hyrho(k) * ql(k,j,i) / nc
[1012]	316	!
	317	!-- Autoconversion rate (Seifert and Beheng, 2006):
[1005]	318	autocon = k_au * ( nu_c + 2.0 ) * ( nu_c + 4.0 ) / &
[1048]	319	( nu_c + 1.0 )*2 ql(k,j,i)*2 xc*2 &
	320	( 1.0 + phi_au / ( 1.0 - tau_cloud )*2 ) &
	321	rho_surface
	322	autocon = MIN( autocon, ql(k,j,i) / ( dt_3d * &
	323	weight_substep(intermediate_timestep_count) ) )
[1012]	324	!
[1022]	325	!-- Tendencies for q, qr, nr, pt:
[1005]	326	tend_qr(k,j,i) = tend_qr(k,j,i) + autocon
	327	tend_q(k,j,i) = tend_q(k,j,i) - autocon
	328	tend_nr(k,j,i) = tend_nr(k,j,i) + autocon / x0 * hyrho(k)
[1048]	329	tend_pt(k,j,i) = tend_pt(k,j,i) + autocon * l_d_cp * pt_d_t(k)
[1005]	330	ENDIF
[1000]	331
	332	ENDDO
	333
[1005]	334	END SUBROUTINE autoconversion_ij
	335
	336	SUBROUTINE accretion_ij( i, j )
	337
	338	USE arrays_3d
	339	USE cloud_parameters
	340	USE constants
	341	USE indices
	342	USE control_parameters
[1048]	343	USE statistics
[1005]	344
	345	IMPLICIT NONE
	346
	347	INTEGER :: i, j, k
[1048]	348	REAL :: accr, phi_ac, tau_cloud
[1005]	349
	350	DO k = nzb_2d(j,i)+1, nzt
[1048]	351	IF ( ( ql(k,j,i) > 0.0 ) .AND. ( qr(k,j,i) > eps_sb ) ) THEN
[1012]	352	!
[1048]	353	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
	354	tau_cloud = 1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) + 1.0E-20)
[1012]	355	!
	356	!-- Universal function for accretion process
[1048]	357	!-- (Seifert and Beheng, 2001):
	358	phi_ac = tau_cloud / ( tau_cloud + 5.0E-5 )
	359	phi_ac = ( phi_ac2 )2
[1012]	360	!
	361	!-- Accretion rate (Seifert and Beheng, 2006):
[1005]	362	accr = k_cr * ql(k,j,i) * qr(k,j,i) * phi_ac * &
	363	SQRT( rho_surface * hyrho(k) )
[1048]	364	accr = MIN( accr, ql(k,j,i) / ( dt_3d * &
	365	weight_substep(intermediate_timestep_count) ) )
[1012]	366	!
[1022]	367	!-- Tendencies for q, qr, pt:
[1005]	368	tend_qr(k,j,i) = tend_qr(k,j,i) + accr
	369	tend_q(k,j,i) = tend_q(k,j,i) - accr
[1048]	370	tend_pt(k,j,i) = tend_pt(k,j,i) + accr * l_d_cp * pt_d_t(k)
[1005]	371	ENDIF
	372	ENDDO
	373
[1000]	374	END SUBROUTINE accretion_ij
	375
[1005]	376
	377	SUBROUTINE selfcollection_breakup_ij( i, j )
	378
	379	USE arrays_3d
	380	USE cloud_parameters
	381	USE constants
	382	USE indices
	383	USE control_parameters
[1048]	384	USE statistics
[1005]	385
	386	IMPLICIT NONE
	387
	388	INTEGER :: i, j, k
[1048]	389	REAL :: selfcoll, breakup, phi_br, phi_sc
[1005]	390
	391	DO k = nzb_2d(j,i)+1, nzt
[1048]	392	IF ( ( qr(k,j,i) > eps_sb ) .AND. ( nr(k,j,i) > eps_sb ) ) THEN
[1012]	393	!
	394	!-- Selfcollection rate (Seifert and Beheng, 2006):
[1048]	395	!-- pirho_l**( 1.0 / 3.0 ) is necessary to convert [lambda_r] = m-1 to
	396	!-- kg**( 1.0 / 3.0 ).
	397	phi_sc = ( 1.0 + kappa_rr / lambda_r(k) * &
	398	pirho_l( 1.0 / 3.0 ) )( -9 )
	399
	400	selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * phi_sc * &
[1005]	401	SQRT( hyrho(k) * rho_surface )
[1012]	402	!
[1048]	403	!-- Collisional breakup rate (Seifert, 2008):
	404	IF ( dr(k) >= 0.3E-3 ) THEN
	405	phi_br = k_br * ( dr(k) - 1.1E-3 )
[1005]	406	breakup = selfcoll * ( phi_br + 1.0 )
	407	ELSE
	408	breakup = 0.0
	409	ENDIF
[1048]	410
	411	selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / ( dt_3d * &
	412	weight_substep(intermediate_timestep_count) ) )
[1012]	413	!
	414	!-- Tendency for nr:
[1048]	415	tend_nr(k,j,i) = tend_nr(k,j,i) + selfcoll
[1005]	416	ENDIF
	417	ENDDO
	418
	419	END SUBROUTINE selfcollection_breakup_ij
	420
[1012]	421	SUBROUTINE evaporation_rain_ij( i, j )
[1022]	422	!
	423	!-- Evaporation of precipitable water. Condensation is neglected for
	424	!-- precipitable water.
[1012]	425
	426	USE arrays_3d
	427	USE cloud_parameters
	428	USE constants
	429	USE indices
	430	USE control_parameters
[1048]	431	USE statistics
	432
[1012]	433	IMPLICIT NONE
	434
	435	INTEGER :: i, j, k
	436	REAL :: evap, alpha, e_s, q_s, t_l, sat, temp, g_evap, f_vent, &
[1049]	437	mu_r_2, mu_r_5d2, nr_0
[1012]	438
	439	DO k = nzb_2d(j,i)+1, nzt
[1048]	440	IF ( ( qr(k,j,i) > eps_sb ) .AND. ( nr(k,j,i) > eps_sb ) ) THEN
[1012]	441	!
	442	!-- Actual liquid water temperature:
	443	t_l = t_d_pt(k) * pt(k,j,i)
	444	!
	445	!-- Saturation vapor pressure at t_l:
	446	e_s = 610.78 * EXP( 17.269 * ( t_l - 273.16 ) / ( t_l - 35.86 ) )
	447	!
	448	!-- Computation of saturation humidity:
	449	q_s = 0.622 * e_s / ( hyp(k) - 0.378 * e_s )
	450	alpha = 0.622 * l_d_r * l_d_cp / ( t_l * t_l )
	451	q_s = q_s * ( 1.0 + alpha * q(k,j,i) ) / ( 1.0 + alpha * q_s )
	452	!
	453	!-- Oversaturation:
	454	sat = MIN( 0.0, ( q(k,j,i) - ql(k,j,i) ) / q_s - 1.0 )
	455	!
	456	!-- Actual temperature:
[1048]	457	temp = t_l + l_d_cp * ql(k,j,i)
[1012]	458	!
	459	!--
	460	g_evap = ( l_v / ( r_v * temp ) - 1.0 ) * l_v / &
[1048]	461	( thermal_conductivity_l * temp ) + rho_l * r_v * temp /&
[1022]	462	( diff_coeff_l * e_s )
[1012]	463	g_evap = 1.0 / g_evap
	464	!
[1049]	465	!-- Compute ventilation factor and intercept parameter
	466	!-- (Seifert and Beheng, 2006; Seifert, 2008):
[1048]	467	IF ( ventilation_effect ) THEN
	468	mu_r_2 = mu_r(k) + 2.0
	469	mu_r_5d2 = mu_r(k) + 2.5
	470	f_vent = a_vent * gamm( mu_r_2 ) * &
	471	lambda_r(k)**( -mu_r_2 ) + &
	472	b_vent * schmidt_p_1d3 * &
	473	SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * &
	474	lambda_r(k)*( -mu_r_5d2 ) &
	475	( 1.0 - 0.5 * ( b_term / a_term ) * &
	476	( lambda_r(k) / &
	477	( c_term + lambda_r(k) ) )**mu_r_5d2 - &
	478	0.125 * ( b_term / a_term )*2 &
	479	( lambda_r(k) / &
	480	( 2.0 * c_term + lambda_r(k) ) )**mu_r_5d2 - &
	481	0.0625 * ( b_term / a_term )*3 &
	482	( lambda_r(k) / &
	483	( 3.0 * c_term + lambda_r(k) ) )**mu_r_5d2 - &
	484	0.0390625 * ( b_term / a_term )*4 &
	485	( lambda_r(k) / &
	486	( 4.0 * c_term + lambda_r(k) ) )**mu_r_5d2 )
[1049]	487	nr_0 = nr(k,j,i) * lambda_r(k)**( mu_r(k) + 1.0 ) / &
	488	gamm( mu_r(k) + 1.0 )
[1048]	489	ELSE
	490	f_vent = 1.0
[1049]	491	nr_0 = nr(k,j,i) * dr(k)
[1048]	492	ENDIF
[1012]	493	!
[1048]	494	!-- Evaporation rate of rain water content (Seifert and Beheng, 2006):
[1049]	495	evap = 2.0 * pi * nr_0 * g_evap * f_vent * sat / &
[1048]	496	hyrho(k)
	497	evap = MAX( evap, -qr(k,j,i) / ( dt_3d * &
	498	weight_substep(intermediate_timestep_count) ) )
[1012]	499	!
[1022]	500	!-- Tendencies for q, qr, nr, pt:
[1012]	501	tend_qr(k,j,i) = tend_qr(k,j,i) + evap
	502	tend_q(k,j,i) = tend_q(k,j,i) - evap
[1049]	503	tend_nr(k,j,i) = tend_nr(k,j,i) + c_evap * evap / xr(k) * hyrho(k)
[1048]	504	tend_pt(k,j,i) = tend_pt(k,j,i) + evap * l_d_cp * pt_d_t(k)
[1012]	505	ENDIF
	506	ENDDO
	507
	508	END SUBROUTINE evaporation_rain_ij
	509
	510	SUBROUTINE sedimentation_cloud_ij( i, j )
	511
	512	USE arrays_3d
	513	USE cloud_parameters
	514	USE constants
	515	USE indices
	516	USE control_parameters
	517
	518	IMPLICIT NONE
	519
	520	INTEGER :: i, j, k
[1022]	521	REAL :: sed_q_const, sigma_gc = 1.3, k_st = 1.2E8
[1012]	522	!
	523	!-- Sedimentation of cloud droplets (Heus et al., 2010):
[1022]	524	sed_q_const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))*( 2.0 / 3.0 ) &
[1048]	525	EXP( 5.0 * LOG( sigma_gc )**2 )
[1012]	526
[1022]	527	sed_q = 0.0
[1012]	528
	529	DO k = nzb_2d(j,i)+1, nzt
[1048]	530	IF ( ql(k,j,i) > 0.0 ) THEN
[1022]	531	sed_q(k) = sed_q_const * nc*( -2.0 / 3.0 ) &
[1012]	532	( ql(k,j,i) * hyrho(k) )**( 5.0 / 3.0 )
	533	ENDIF
	534	ENDDO
	535	!
[1022]	536	!-- Tendency for q, pt:
[1012]	537	DO k = nzb_2d(j,i)+1, nzt
[1048]	538	tend_q(k,j,i) = tend_q(k,j,i) + ( sed_q(k+1) - sed_q(k) ) * &
	539	ddzu(k+1) / hyrho(k)
[1022]	540	tend_pt(k,j,i) = tend_pt(k,j,i) - ( sed_q(k+1) - sed_q(k) ) * &
[1048]	541	ddzu(k+1) / hyrho(k) * l_d_cp * pt_d_t(k)
[1012]	542	ENDDO
	543
	544	END SUBROUTINE sedimentation_cloud_ij
	545
	546	SUBROUTINE sedimentation_rain_ij( i, j )
	547
	548	USE arrays_3d
	549	USE cloud_parameters
	550	USE constants
	551	USE indices
	552	USE control_parameters
[1048]	553	USE statistics
[1012]	554
	555	IMPLICIT NONE
	556
[1022]	557	INTEGER :: i, j, k, n, n_substep
[1048]	558	REAL :: sed_nr_tend, sed_qr_tend
[1012]	559
[1048]	560	IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0
[1012]	561
[1022]	562	sed_nr = 0.0
	563	sed_qr = 0.0
	564
	565	DO k = nzb_2d(j,i)+1, nzt
[1048]	566	IF ( ( qr(k,j,i) > eps_sb ) .AND. ( nr(k,j,i) > eps_sb ) ) THEN
[1012]	567	!
[1022]	568	!-- Sedimentation of rain water content and rain drop concentration
	569	!-- according to Stevens and Seifert (2008):
[1048]	570	sed_nr(k) = MIN( 20.0, MAX( 0.1, a_term - b_term * ( 1.0 + &
	571	c_term / lambda_r(k) )*( -1.0 &
	572	( mu_r(k) + 1.0 ) ) ) ) * nr(k,j,i)
	573	sed_qr(k) = MIN( 20.0, MAX( 0.1, a_term - b_term * ( 1.0 + &
	574	c_term / lambda_r(k) )*( -1.0 &
	575	( mu_r(k) + 4.0 ) ) ) ) * qr(k,j,i) * hyrho(k)
	576	!
	577	!-- Computation of rain rate
	578	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * &
	579	weight_substep(intermediate_timestep_count)
[1022]	580	ENDIF
	581	ENDDO
	582
	583	DO k = nzb_2d(j,i)+1, nzt
[1048]	584	sed_nr_tend = MAX( ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1), &
	585	-nr(k,j,i) / ( dt_3d * &
	586	weight_substep(intermediate_timestep_count) ) )
	587	sed_qr_tend = MAX( ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
	588	hyrho(k), &
	589	-qr(k,j,i) / ( dt_3d * &
	590	weight_substep(intermediate_timestep_count) ) )
	591
	592	tend_nr(k,j,i) = tend_nr(k,j,i) + sed_nr_tend
	593	tend_qr(k,j,i) = tend_qr(k,j,i) + sed_qr_tend
[1022]	594	ENDDO
	595	!
[1048]	596	!-- Precipitation amount
	597	IF ( intermediate_timestep_count == intermediate_timestep_count_max &
	598	.AND. ( dt_do2d_xy - time_do2d_xy ) < &
	599	precipitation_amount_interval ) THEN
[1012]	600
[1048]	601	precipitation_amount(j,i) = precipitation_amount(j,i) + &
	602	prr(nzb_2d(j,i)+1,j,i) * &
	603	hyrho(nzb_2d(j,i)+1) * dt_3d
	604	ENDIF
	605
[1012]	606	END SUBROUTINE sedimentation_rain_ij
	607
	608	!
	609	!-- This function computes the gamma function (Press et al., 1992).
	610	!-- The gamma function is needed for the calculation of the evaporation
	611	!-- of rain drops.
	612	FUNCTION gamm( xx )
[1048]	613
	614	USE cloud_parameters
[1012]	615
	616	IMPLICIT NONE
	617
	618	REAL :: gamm, xx, &
[1048]	619	ser, tmp, x_gamm, y_gamm
[1012]	620	INTEGER :: j
	621
	622	x_gamm = xx
	623	y_gamm = x_gamm
	624	tmp = x_gamm + 5.5
	625	tmp = ( x_gamm + 0.5 ) * LOG( tmp ) - tmp
	626	ser = 1.000000000190015
	627	do j = 1, 6
	628	y_gamm = y_gamm + 1.0
	629	ser = ser + cof( j ) / y_gamm
	630	enddo
	631	!
	632	!-- Until this point the algorithm computes the logarithm of the gamma
	633	!-- function. Hence, the exponential function is used.
	634	! gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
	635	gamm = EXP( tmp ) * stp * ser / x_gamm
	636	RETURN
	637
	638	END FUNCTION gamm
	639
	640	END MODULE microphysics_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |