Home

Context Navigation

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

Last change on this file since 1092 was 1092, checked in by raasch, 11 years ago
unused variables remove from several routines
Property svn:keywords set to `Id`
File size: 24.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |