Home

Context Navigation

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

Last change on this file since 1065 was 1065, checked in by hoffmann, 11 years ago
cloud physics: rain sedimentation and turbulence effects
Property svn:keywords set to `Id`
File size: 24.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |