Home

Context Navigation

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

Last change on this file since 1074 was 1066, checked in by hoffmann, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 24.6 KB

Rev	Line
[1000]	1	MODULE microphysics_mod
	2
	3	!------------------------------------------------------------------------------!
	4	! Current revisions:
	5	! -----------------
	6	!
	7	! Former revisions:
	8	! -----------------
[1052]	9	! $Id: microphysics.f90 1066 2012-11-22 17:52:43Z hoffmann $
[1054]	10	!
[1066]	11	! 1065 2012-11-22 17:42:36Z hoffmann
	12	! Sedimentation process implemented according to Stevens and Seifert (2008).
	13	! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
	14	! and Stevens, 2010).
	15	!
[1054]	16	! 1053 2012-11-13 17:11:03Z hoffmann
	17	! initial revision
[1000]	18	!
	19	! Description:
	20	! ------------
	21	! Calculate cloud microphysics according to the two moment bulk
	22	! scheme by Seifert and Beheng (2006).
	23	!------------------------------------------------------------------------------!
	24
	25	PRIVATE
[1022]	26	PUBLIC dsd_properties, autoconversion, accretion, selfcollection_breakup, &
[1012]	27	evaporation_rain, sedimentation_cloud, sedimentation_rain
[1000]	28
[1022]	29	INTERFACE dsd_properties
	30	MODULE PROCEDURE dsd_properties
	31	MODULE PROCEDURE dsd_properties_ij
	32	END INTERFACE dsd_properties
	33
[1000]	34	INTERFACE autoconversion
	35	MODULE PROCEDURE autoconversion
	36	MODULE PROCEDURE autoconversion_ij
	37	END INTERFACE autoconversion
	38
	39	INTERFACE accretion
	40	MODULE PROCEDURE accretion
	41	MODULE PROCEDURE accretion_ij
	42	END INTERFACE accretion
[1005]	43
	44	INTERFACE selfcollection_breakup
	45	MODULE PROCEDURE selfcollection_breakup
	46	MODULE PROCEDURE selfcollection_breakup_ij
	47	END INTERFACE selfcollection_breakup
[1012]	48
	49	INTERFACE evaporation_rain
	50	MODULE PROCEDURE evaporation_rain
	51	MODULE PROCEDURE evaporation_rain_ij
	52	END INTERFACE evaporation_rain
	53
	54	INTERFACE sedimentation_cloud
	55	MODULE PROCEDURE sedimentation_cloud
	56	MODULE PROCEDURE sedimentation_cloud_ij
	57	END INTERFACE sedimentation_cloud
[1000]	58
[1012]	59	INTERFACE sedimentation_rain
	60	MODULE PROCEDURE sedimentation_rain
	61	MODULE PROCEDURE sedimentation_rain_ij
	62	END INTERFACE sedimentation_rain
	63
[1000]	64	CONTAINS
	65
	66
	67	!------------------------------------------------------------------------------!
	68	! Call for all grid points
	69	!------------------------------------------------------------------------------!
[1022]	70	SUBROUTINE dsd_properties
	71
	72	USE arrays_3d
	73	USE cloud_parameters
	74	USE constants
	75	USE indices
	76
	77	IMPLICIT NONE
	78
	79	INTEGER :: i, j, k
	80	REAL :: dqdt_precip
	81
	82
	83	DO i = nxl, nxr
	84	DO j = nys, nyn
	85	DO k = nzb_2d(j,i)+1, nzt
	86
	87	ENDDO
	88	ENDDO
	89	ENDDO
	90
	91	END SUBROUTINE dsd_properties
	92
[1000]	93	SUBROUTINE autoconversion
	94
	95	USE arrays_3d
	96	USE cloud_parameters
	97	USE constants
	98	USE indices
	99
	100	IMPLICIT NONE
	101
	102	INTEGER :: i, j, k
	103	REAL :: dqdt_precip
	104
	105
	106	DO i = nxl, nxr
	107	DO j = nys, nyn
	108	DO k = nzb_2d(j,i)+1, nzt
	109
	110	ENDDO
	111	ENDDO
	112	ENDDO
	113
	114	END SUBROUTINE autoconversion
	115
[1005]	116	SUBROUTINE accretion
[1000]	117
	118	USE arrays_3d
	119	USE cloud_parameters
	120	USE constants
	121	USE indices
[1005]	122
[1000]	123	IMPLICIT NONE
	124
	125	INTEGER :: i, j, k
[1005]	126	REAL :: dqdt_precip
[1000]	127
[1005]	128
	129	DO i = nxl, nxr
	130	DO j = nys, nyn
	131	DO k = nzb_2d(j,i)+1, nzt
[1000]	132
[1005]	133	ENDDO
	134	ENDDO
[1000]	135	ENDDO
	136
[1005]	137	END SUBROUTINE accretion
[1000]	138
[1005]	139	SUBROUTINE selfcollection_breakup
[1000]	140
	141	USE arrays_3d
	142	USE cloud_parameters
	143	USE constants
	144	USE indices
	145
	146	IMPLICIT NONE
	147
	148	INTEGER :: i, j, k
	149	REAL :: dqdt_precip
	150
	151
	152	DO i = nxl, nxr
	153	DO j = nys, nyn
	154	DO k = nzb_2d(j,i)+1, nzt
	155
	156	ENDDO
	157	ENDDO
	158	ENDDO
	159
[1005]	160	END SUBROUTINE selfcollection_breakup
[1000]	161
[1012]	162	SUBROUTINE evaporation_rain
[1000]	163
[1012]	164	USE arrays_3d
	165	USE cloud_parameters
	166	USE constants
	167	USE indices
	168
	169	IMPLICIT NONE
	170
	171	INTEGER :: i, j, k
	172	REAL :: dqdt_precip
	173
	174
	175	DO i = nxl, nxr
	176	DO j = nys, nyn
	177	DO k = nzb_2d(j,i)+1, nzt
	178
	179	ENDDO
	180	ENDDO
	181	ENDDO
	182
	183	END SUBROUTINE evaporation_rain
	184
	185	SUBROUTINE sedimentation_cloud
	186
	187	USE arrays_3d
	188	USE cloud_parameters
	189	USE constants
	190	USE indices
	191
	192	IMPLICIT NONE
	193
	194	INTEGER :: i, j, k
	195	REAL :: dqdt_precip
	196
	197
	198	DO i = nxl, nxr
	199	DO j = nys, nyn
	200	DO k = nzb_2d(j,i)+1, nzt
	201
	202	ENDDO
	203	ENDDO
	204	ENDDO
	205
	206	END SUBROUTINE sedimentation_cloud
	207
	208	SUBROUTINE sedimentation_rain
	209
	210	USE arrays_3d
	211	USE cloud_parameters
	212	USE constants
	213	USE indices
	214
	215	IMPLICIT NONE
	216
	217	INTEGER :: i, j, k
	218	REAL :: dqdt_precip
	219
	220
	221	DO i = nxl, nxr
	222	DO j = nys, nyn
	223	DO k = nzb_2d(j,i)+1, nzt
	224
	225	ENDDO
	226	ENDDO
	227	ENDDO
	228
	229
	230	END SUBROUTINE sedimentation_rain
	231
	232
[1000]	233	!------------------------------------------------------------------------------!
	234	! Call for grid point i,j
	235	!------------------------------------------------------------------------------!
[1022]	236	SUBROUTINE dsd_properties_ij( i, j )
	237
	238	USE arrays_3d
	239	USE cloud_parameters
	240	USE constants
	241	USE indices
	242	USE control_parameters
[1048]	243	USE user
[1022]	244
	245	IMPLICIT NONE
	246
	247	INTEGER :: i, j, k
	248
	249	DO k = nzb_2d(j,i)+1, nzt
[1065]	250
	251	IF ( qr(k,j,i) <= eps_sb ) THEN
	252	qr(k,j,i) = 0.0
	253	ELSE
[1022]	254	!
[1048]	255	!-- Adjust number of raindrops to avoid nonlinear effects in
	256	!-- sedimentation and evaporation of rain drops due to too small or
[1065]	257	!-- too big weights of rain drops (Stevens and Seifert, 2008).
	258	IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN
	259	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
	260	ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN
	261	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
[1048]	262	ENDIF
[1065]	263	xr(k) = hyrho(k) * qr(k,j,i) / nr(k,j,i)
[1022]	264	!
[1065]	265	!-- Weight averaged diameter of rain drops:
	266	dr(k) = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * &
	267	dpirho_l )**( 1.0 / 3.0 )
[1048]	268	!
[1049]	269	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	270	!-- Stevens and Seifert, 2008):
[1065]	271	mu_r(k) = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr(k) - 1.4E-3 ) ) )
[1049]	272	!
[1022]	273	!-- Slope parameter of gamma distribution (Seifert, 2008):
[1048]	274	lambda_r(k) = ( ( mu_r(k) + 3.0 ) * ( mu_r(k) + 2.0 ) * &
	275	( mu_r(k) + 1.0 ) )**( 1.0 / 3.0 ) / dr(k)
[1022]	276	ENDIF
	277	ENDDO
	278
	279	END SUBROUTINE dsd_properties_ij
	280
[1005]	281	SUBROUTINE autoconversion_ij( i, j )
[1000]	282
	283	USE arrays_3d
	284	USE cloud_parameters
	285	USE constants
	286	USE indices
[1005]	287	USE control_parameters
[1048]	288	USE statistics
[1065]	289	USE grid_variables
[1000]	290
	291	IMPLICIT NONE
	292
	293	INTEGER :: i, j, k
[1065]	294	REAL :: k_au, autocon, phi_au, tau_cloud, xc, nu_c, rc, &
	295	l_mix, re_lambda, alpha_cc, r_cc, sigma_cc, epsilon
[1000]	296
[1005]	297	k_au = k_cc / ( 20.0 * x0 )
	298
[1000]	299	DO k = nzb_2d(j,i)+1, nzt
	300
[1048]	301	IF ( ql(k,j,i) > 0.0 ) THEN
[1012]	302	!
[1048]	303	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[1012]	304	!-- (1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) ))
[1048]	305	tau_cloud = 1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) + 1.0E-20 )
[1012]	306	!
	307	!-- Universal function for autoconversion process
	308	!-- (Seifert and Beheng, 2006):
[1048]	309	phi_au = 600.0 * tau_cloud*0.68 ( 1.0 - tau_cloud0.68 )3
[1012]	310	!
	311	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
	312	!-- (Use constant nu_c = 1.0 instead?)
[1065]	313	nu_c = 1.0 !MAX( 0.0, 1580.0 * hyrho(k) * ql(k,j,i) - 0.28 )
[1012]	314	!
	315	!-- Mean weight of cloud droplets:
[1005]	316	xc = hyrho(k) * ql(k,j,i) / nc
[1012]	317	!
[1065]	318	!-- Parameterized turbulence effects on autoconversion (Seifert,
	319	!-- Nuijens and Stevens, 2010)
	320	IF ( turbulence ) THEN
	321	!
	322	!-- Weight averaged radius of cloud droplets:
	323	rc = 0.5 * ( xc * dpirho_l )**( 1.0 / 3.0 )
	324
	325	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0 + a_3 * nu_c )
	326	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0 + b_3 * nu_c )
	327	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0 + c_3 * nu_c )
	328	!
	329	!-- Mixing length (neglecting distance to ground and stratification)
	330	l_mix = ( dx * dy * dzu(k) )**( 1.0 / 3.0 )
	331	!
	332	!-- Limit dissipation rate according to Seifert, Nuijens and
	333	!-- Stevens (2010)
	334	epsilon = MIN( 0.06, diss(k,j,i) )
	335	!
	336	!-- Compute Taylor-microscale Reynolds number:
	337	re_lambda = 6.0 / 11.0 * ( l_mix / c_const )*( 2.0 / 3.0 ) &
	338	SQRT( 15.0 / kin_vis_air ) * epsilon**( 1.0 / 6.0 )
	339	!
	340	!-- The factor of 1.0E4 is needed to convert the dissipation rate
	341	!-- from m2 s-3 to cm2 s-3.
	342	k_au = k_au * ( 1.0 + &
	343	epsilon * 1.0E4 * ( re_lambda * 1.0E-3 )*0.25 &
	344	( alpha_cc * EXP( -1.0 * ( ( rc - r_cc ) / &
	345	sigma_cc )**2 ) + beta_cc ) )
	346	ENDIF
	347	!
[1012]	348	!-- Autoconversion rate (Seifert and Beheng, 2006):
[1065]	349	autocon = k_au * ( nu_c + 2.0 ) * ( nu_c + 4.0 ) / &
	350	( nu_c + 1.0 )*2 ql(k,j,i)*2 xc*2 &
	351	( 1.0 + phi_au / ( 1.0 - tau_cloud )*2 ) &
[1048]	352	rho_surface
[1065]	353	autocon = MIN( autocon, ql(k,j,i) / ( dt_3d * &
[1048]	354	weight_substep(intermediate_timestep_count) ) )
[1012]	355	!
[1022]	356	!-- Tendencies for q, qr, nr, pt:
[1005]	357	tend_qr(k,j,i) = tend_qr(k,j,i) + autocon
	358	tend_q(k,j,i) = tend_q(k,j,i) - autocon
	359	tend_nr(k,j,i) = tend_nr(k,j,i) + autocon / x0 * hyrho(k)
[1048]	360	tend_pt(k,j,i) = tend_pt(k,j,i) + autocon * l_d_cp * pt_d_t(k)
[1005]	361	ENDIF
[1000]	362
	363	ENDDO
	364
[1005]	365	END SUBROUTINE autoconversion_ij
	366
	367	SUBROUTINE accretion_ij( i, j )
	368
	369	USE arrays_3d
	370	USE cloud_parameters
	371	USE constants
	372	USE indices
	373	USE control_parameters
[1048]	374	USE statistics
[1005]	375
	376	IMPLICIT NONE
	377
	378	INTEGER :: i, j, k
[1065]	379	REAL :: accr, phi_ac, tau_cloud, k_cr
[1005]	380
	381	DO k = nzb_2d(j,i)+1, nzt
[1048]	382	IF ( ( ql(k,j,i) > 0.0 ) .AND. ( qr(k,j,i) > eps_sb ) ) THEN
[1012]	383	!
[1048]	384	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
	385	tau_cloud = 1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) + 1.0E-20)
[1012]	386	!
	387	!-- Universal function for accretion process
[1048]	388	!-- (Seifert and Beheng, 2001):
[1065]	389	phi_ac = tau_cloud / ( tau_cloud + 5.0E-5 )
	390	phi_ac = ( phi_ac2 )2
[1012]	391	!
[1065]	392	!-- Parameterized turbulence effects on autoconversion (Seifert,
	393	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
	394	!-- convert the dissipation (diss) from m2 s-3 to cm2 s-3.
	395	IF ( turbulence ) THEN
	396	k_cr = k_cr0 * ( 1.0 + 0.05 * &
	397	MIN( 600.0, diss(k,j,i) * 1.0E4 )**0.25 )
	398	ELSE
	399	k_cr = k_cr0
	400	ENDIF
	401	!
[1012]	402	!-- Accretion rate (Seifert and Beheng, 2006):
[1065]	403	accr = k_cr * ql(k,j,i) * qr(k,j,i) * phi_ac * &
	404	SQRT( rho_surface * hyrho(k) )
	405	accr = MIN( accr, ql(k,j,i) / ( dt_3d * &
	406	weight_substep(intermediate_timestep_count) ) )
[1012]	407	!
[1022]	408	!-- Tendencies for q, qr, pt:
[1005]	409	tend_qr(k,j,i) = tend_qr(k,j,i) + accr
	410	tend_q(k,j,i) = tend_q(k,j,i) - accr
[1048]	411	tend_pt(k,j,i) = tend_pt(k,j,i) + accr * l_d_cp * pt_d_t(k)
[1005]	412	ENDIF
	413	ENDDO
	414
[1000]	415	END SUBROUTINE accretion_ij
	416
[1005]	417
	418	SUBROUTINE selfcollection_breakup_ij( i, j )
	419
	420	USE arrays_3d
	421	USE cloud_parameters
	422	USE constants
	423	USE indices
	424	USE control_parameters
[1048]	425	USE statistics
[1005]	426
	427	IMPLICIT NONE
	428
	429	INTEGER :: i, j, k
[1048]	430	REAL :: selfcoll, breakup, phi_br, phi_sc
[1005]	431
	432	DO k = nzb_2d(j,i)+1, nzt
[1065]	433	IF ( qr(k,j,i) > eps_sb ) THEN
[1012]	434	!
	435	!-- Selfcollection rate (Seifert and Beheng, 2006):
[1048]	436	!-- pirho_l**( 1.0 / 3.0 ) is necessary to convert [lambda_r] = m-1 to
	437	!-- kg**( 1.0 / 3.0 ).
[1065]	438	phi_sc = 1.0 !( 1.0 + kappa_rr / lambda_r(k) * &
	439	!pirho_l( 1.0 / 3.0 ) )( -9 )
[1048]	440
	441	selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * phi_sc * &
[1005]	442	SQRT( hyrho(k) * rho_surface )
[1012]	443	!
[1048]	444	!-- Collisional breakup rate (Seifert, 2008):
	445	IF ( dr(k) >= 0.3E-3 ) THEN
	446	phi_br = k_br * ( dr(k) - 1.1E-3 )
[1005]	447	breakup = selfcoll * ( phi_br + 1.0 )
	448	ELSE
	449	breakup = 0.0
	450	ENDIF
[1048]	451
	452	selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / ( dt_3d * &
	453	weight_substep(intermediate_timestep_count) ) )
[1012]	454	!
	455	!-- Tendency for nr:
[1048]	456	tend_nr(k,j,i) = tend_nr(k,j,i) + selfcoll
[1005]	457	ENDIF
	458	ENDDO
	459
	460	END SUBROUTINE selfcollection_breakup_ij
	461
[1012]	462	SUBROUTINE evaporation_rain_ij( i, j )
[1022]	463	!
	464	!-- Evaporation of precipitable water. Condensation is neglected for
	465	!-- precipitable water.
[1012]	466
	467	USE arrays_3d
	468	USE cloud_parameters
	469	USE constants
	470	USE indices
	471	USE control_parameters
[1048]	472	USE statistics
	473
[1012]	474	IMPLICIT NONE
	475
	476	INTEGER :: i, j, k
	477	REAL :: evap, alpha, e_s, q_s, t_l, sat, temp, g_evap, f_vent, &
[1049]	478	mu_r_2, mu_r_5d2, nr_0
[1012]	479
	480	DO k = nzb_2d(j,i)+1, nzt
[1065]	481	IF ( qr(k,j,i) > eps_sb ) THEN
[1012]	482	!
	483	!-- Actual liquid water temperature:
	484	t_l = t_d_pt(k) * pt(k,j,i)
	485	!
	486	!-- Saturation vapor pressure at t_l:
	487	e_s = 610.78 * EXP( 17.269 * ( t_l - 273.16 ) / ( t_l - 35.86 ) )
	488	!
	489	!-- Computation of saturation humidity:
	490	q_s = 0.622 * e_s / ( hyp(k) - 0.378 * e_s )
	491	alpha = 0.622 * l_d_r * l_d_cp / ( t_l * t_l )
	492	q_s = q_s * ( 1.0 + alpha * q(k,j,i) ) / ( 1.0 + alpha * q_s )
	493	!
	494	!-- Oversaturation:
	495	sat = MIN( 0.0, ( q(k,j,i) - ql(k,j,i) ) / q_s - 1.0 )
	496	!
	497	!-- Actual temperature:
[1048]	498	temp = t_l + l_d_cp * ql(k,j,i)
[1012]	499	!
	500	!--
	501	g_evap = ( l_v / ( r_v * temp ) - 1.0 ) * l_v / &
[1048]	502	( thermal_conductivity_l * temp ) + rho_l * r_v * temp /&
[1022]	503	( diff_coeff_l * e_s )
[1012]	504	g_evap = 1.0 / g_evap
	505	!
[1049]	506	!-- Compute ventilation factor and intercept parameter
	507	!-- (Seifert and Beheng, 2006; Seifert, 2008):
[1048]	508	IF ( ventilation_effect ) THEN
	509	mu_r_2 = mu_r(k) + 2.0
	510	mu_r_5d2 = mu_r(k) + 2.5
	511	f_vent = a_vent * gamm( mu_r_2 ) * &
	512	lambda_r(k)**( -mu_r_2 ) + &
	513	b_vent * schmidt_p_1d3 * &
	514	SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * &
	515	lambda_r(k)*( -mu_r_5d2 ) &
	516	( 1.0 - 0.5 * ( b_term / a_term ) * &
	517	( lambda_r(k) / &
	518	( c_term + lambda_r(k) ) )**mu_r_5d2 - &
	519	0.125 * ( b_term / a_term )*2 &
	520	( lambda_r(k) / &
	521	( 2.0 * c_term + lambda_r(k) ) )**mu_r_5d2 - &
	522	0.0625 * ( b_term / a_term )*3 &
	523	( lambda_r(k) / &
	524	( 3.0 * c_term + lambda_r(k) ) )**mu_r_5d2 - &
	525	0.0390625 * ( b_term / a_term )*4 &
	526	( lambda_r(k) / &
	527	( 4.0 * c_term + lambda_r(k) ) )**mu_r_5d2 )
[1049]	528	nr_0 = nr(k,j,i) * lambda_r(k)**( mu_r(k) + 1.0 ) / &
	529	gamm( mu_r(k) + 1.0 )
[1048]	530	ELSE
	531	f_vent = 1.0
[1049]	532	nr_0 = nr(k,j,i) * dr(k)
[1048]	533	ENDIF
[1012]	534	!
[1048]	535	!-- Evaporation rate of rain water content (Seifert and Beheng, 2006):
[1049]	536	evap = 2.0 * pi * nr_0 * g_evap * f_vent * sat / &
[1048]	537	hyrho(k)
	538	evap = MAX( evap, -qr(k,j,i) / ( dt_3d * &
	539	weight_substep(intermediate_timestep_count) ) )
[1012]	540	!
[1022]	541	!-- Tendencies for q, qr, nr, pt:
[1012]	542	tend_qr(k,j,i) = tend_qr(k,j,i) + evap
	543	tend_q(k,j,i) = tend_q(k,j,i) - evap
[1049]	544	tend_nr(k,j,i) = tend_nr(k,j,i) + c_evap * evap / xr(k) * hyrho(k)
[1048]	545	tend_pt(k,j,i) = tend_pt(k,j,i) + evap * l_d_cp * pt_d_t(k)
[1012]	546	ENDIF
	547	ENDDO
	548
	549	END SUBROUTINE evaporation_rain_ij
	550
	551	SUBROUTINE sedimentation_cloud_ij( i, j )
	552
	553	USE arrays_3d
	554	USE cloud_parameters
	555	USE constants
	556	USE indices
	557	USE control_parameters
	558
	559	IMPLICIT NONE
	560
	561	INTEGER :: i, j, k
[1022]	562	REAL :: sed_q_const, sigma_gc = 1.3, k_st = 1.2E8
[1012]	563	!
	564	!-- Sedimentation of cloud droplets (Heus et al., 2010):
[1022]	565	sed_q_const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))*( 2.0 / 3.0 ) &
[1048]	566	EXP( 5.0 * LOG( sigma_gc )**2 )
[1012]	567
[1022]	568	sed_q = 0.0
[1012]	569
	570	DO k = nzb_2d(j,i)+1, nzt
[1048]	571	IF ( ql(k,j,i) > 0.0 ) THEN
[1022]	572	sed_q(k) = sed_q_const * nc*( -2.0 / 3.0 ) &
[1012]	573	( ql(k,j,i) * hyrho(k) )**( 5.0 / 3.0 )
	574	ENDIF
	575	ENDDO
	576	!
[1022]	577	!-- Tendency for q, pt:
[1012]	578	DO k = nzb_2d(j,i)+1, nzt
[1048]	579	tend_q(k,j,i) = tend_q(k,j,i) + ( sed_q(k+1) - sed_q(k) ) * &
	580	ddzu(k+1) / hyrho(k)
[1022]	581	tend_pt(k,j,i) = tend_pt(k,j,i) - ( sed_q(k+1) - sed_q(k) ) * &
[1048]	582	ddzu(k+1) / hyrho(k) * l_d_cp * pt_d_t(k)
[1012]	583	ENDDO
	584
	585	END SUBROUTINE sedimentation_cloud_ij
	586
	587	SUBROUTINE sedimentation_rain_ij( i, j )
	588
	589	USE arrays_3d
	590	USE cloud_parameters
	591	USE constants
	592	USE indices
	593	USE control_parameters
[1048]	594	USE statistics
[1012]	595
	596	IMPLICIT NONE
	597
[1065]	598	INTEGER :: i, j, k, k_run, n, n_substep
	599	REAL :: c_run, d_max, d_mean, d_min, dt_sedi, flux, mean, z_run
[1012]	600
[1065]	601	REAL, DIMENSION(nzb:nzt) :: c_nr, c_qr, d_nr, d_qr, nr_slope, qr_slope, &
	602	w_nr, w_qr
	603	!
	604	!-- Computation of sedimentation flux. Implementation according to Stevens
	605	!-- and Seifert (2008).
	606
[1048]	607	IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0
[1012]	608
[1065]	609	dt_sedi = dt_3d * weight_substep(intermediate_timestep_count)
[1022]	610
[1065]	611	w_nr = 0.0
	612	w_qr = 0.0
[1012]	613	!
[1065]	614	!-- Compute velocities
	615	DO k = nzb_s_inner(j,i)+1, nzt
	616	IF ( qr(k,j,i) > eps_sb ) THEN
	617	w_nr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 + &
	618	c_term / lambda_r(k) )*( -1.0 ( mu_r(k) + 1.0 ) ) ) )
	619	w_qr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 + &
	620	c_term / lambda_r(k) )*( -1.0 ( mu_r(k) + 4.0 ) ) ) )
	621	ELSE
	622	w_nr(k) = 0.0
	623	w_qr(k) = 0.0
	624	ENDIF
	625	ENDDO
[1048]	626	!
[1065]	627	!-- Adjust boundary values
	628	w_nr(nzb_2d(j,i)) = w_nr(nzb_2d(j,i)+1)
	629	w_qr(nzb_2d(j,i)) = w_qr(nzb_2d(j,i)+1)
	630	w_nr(nzt) = w_nr(nzt-1)
	631	w_qr(nzt) = w_qr(nzt-1)
	632	!
	633	!-- Compute Courant number
	634	DO k = nzb_s_inner(j,i)+1, nzt-1
	635	c_nr(k) = 0.25 * ( w_nr(k-1) + 2.0 * w_nr(k) + w_nr(k+1) ) * &
	636	dt_sedi * ddzu(k)
	637	c_qr(k) = 0.25 * ( w_qr(k-1) + 2.0 * w_qr(k) + w_qr(k+1) ) * &
	638	dt_sedi * ddzu(k)
	639	ENDDO
	640	!
	641	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
	642	IF ( limiter_sedimentation ) THEN
	643
	644	qr(nzb_s_inner(j,i),j,i) = 0.0
	645	nr(nzb_s_inner(j,i),j,i) = 0.0
	646	qr(nzt,j,i) = 0.0
	647	nr(nzt,j,i) = 0.0
	648
	649	DO k = nzb_s_inner(j,i)+1, nzt-1
	650	d_mean = 0.5 * ( qr(k+1,j,i) + qr(k-1,j,i) )
	651	d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
	652	d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
	653
	654	qr_slope(k) = SIGN(1.0, d_mean) * MIN ( 2.0 * d_min, 2.0 * d_max, &
	655	ABS( d_mean ) )
	656
	657	d_mean = 0.5 * ( nr(k+1,j,i) + nr(k-1,j,i) )
	658	d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
	659	d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
	660
	661	nr_slope(k) = SIGN(1.0, d_mean) * MIN ( 2.0 * d_min, 2.0 * d_max, &
	662	ABS( d_mean ) )
[1022]	663	ENDDO
[1048]	664
[1065]	665	ELSE
	666	nr_slope = 0.0
	667	qr_slope = 0.0
	668	ENDIF
	669	!
	670	!-- Compute sedimentation flux
	671	DO k = nzt-2, nzb_s_inner(j,i)+1, -1
	672	!
	673	!-- Sum up all rain drop number densities which contribute to the flux
	674	!-- through k-1/2
	675	flux = 0.0
	676	z_run = 0.0 ! height above z(k)
	677	k_run = k
	678	c_run = MIN( 1.0, c_nr(k) )
	679	DO WHILE ( c_run > 0.0 .AND. k_run <= nzt-1 )
	680	flux = flux + hyrho(k_run) * &
	681	( nr(k_run,j,i) + nr_slope(k_run) * ( 1.0 - c_run ) * &
	682	0.5 ) * c_run * dzu(k_run)
	683	z_run = z_run + dzu(k_run)
	684	k_run = k_run + 1
	685	c_run = MIN( 1.0, c_nr(k_run) - z_run * ddzu(k_run) )
[1022]	686	ENDDO
	687	!
[1065]	688	!-- It is not allowed to sediment more rain drop number density than
	689	!-- available
	690	flux = MIN( flux, &
	691	hyrho(k) * dzu(k) * nr(k,j,i) + sed_nr(k+1) * dt_sedi )
	692
	693	sed_nr(k) = flux / dt_sedi
	694	tend_nr(k,j,i) = tend_nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * &
	695	ddzu(k+1) / hyrho(k)
	696	!
	697	!-- Sum up all rain water content which contributes to the flux
	698	!-- through k-1/2
	699	flux = 0.0
	700	z_run = 0.0 ! height above z(k)
	701	k_run = k
	702	c_run = MIN( 1.0, c_qr(k) )
	703	DO WHILE ( c_run > 0.0 .AND. k_run <= nzt-1 )
	704	flux = flux + hyrho(k_run) * &
	705	( qr(k_run,j,i) + qr_slope(k_run) * ( 1.0 - c_run ) * &
	706	0.5 ) * c_run * dzu(k_run)
	707	z_run = z_run + dzu(k_run)
	708	k_run = k_run + 1
	709	c_run = MIN( 1.0, c_qr(k_run) - z_run * ddzu(k_run) )
	710	ENDDO
	711	!
	712	!-- It is not allowed to sediment more rain water content than available
	713	flux = MIN( flux, &
	714	hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * dt_sedi )
	715
	716	sed_qr(k) = flux / dt_sedi
	717	tend_qr(k,j,i) = tend_qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
	718	ddzu(k+1) / hyrho(k)
	719	!
	720	!-- Compute the rain rate
	721	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * &
	722	weight_substep(intermediate_timestep_count)
	723	ENDDO
	724	!
[1048]	725	!-- Precipitation amount
	726	IF ( intermediate_timestep_count == intermediate_timestep_count_max &
	727	.AND. ( dt_do2d_xy - time_do2d_xy ) < &
	728	precipitation_amount_interval ) THEN
[1012]	729
[1048]	730	precipitation_amount(j,i) = precipitation_amount(j,i) + &
	731	prr(nzb_2d(j,i)+1,j,i) * &
	732	hyrho(nzb_2d(j,i)+1) * dt_3d
	733	ENDIF
	734
[1012]	735	END SUBROUTINE sedimentation_rain_ij
	736
	737	!
	738	!-- This function computes the gamma function (Press et al., 1992).
	739	!-- The gamma function is needed for the calculation of the evaporation
	740	!-- of rain drops.
	741	FUNCTION gamm( xx )
[1048]	742
	743	USE cloud_parameters
[1012]	744
	745	IMPLICIT NONE
	746
[1065]	747	REAL :: gamm, ser, tmp, x_gamm, xx, y_gamm
[1012]	748	INTEGER :: j
	749
	750	x_gamm = xx
	751	y_gamm = x_gamm
	752	tmp = x_gamm + 5.5
	753	tmp = ( x_gamm + 0.5 ) * LOG( tmp ) - tmp
	754	ser = 1.000000000190015
	755	do j = 1, 6
	756	y_gamm = y_gamm + 1.0
	757	ser = ser + cof( j ) / y_gamm
	758	enddo
	759	!
	760	!-- Until this point the algorithm computes the logarithm of the gamma
	761	!-- function. Hence, the exponential function is used.
	762	! gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
	763	gamm = EXP( tmp ) * stp * ser / x_gamm
	764	RETURN
	765
	766	END FUNCTION gamm
	767
	768	END MODULE microphysics_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |