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

Line
1	MODULE microphysics_mod
2
3	!------------------------------------------------------------------------------!
4	! Current revisions:
5	! -----------------
6	!
7	! Former revisions:
8	! -----------------
9	! $Id: microphysics.f90 1066 2012-11-22 17:52:43Z hoffmann $
10	!
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	!
16	! 1053 2012-11-13 17:11:03Z hoffmann
17	! initial revision
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
26	PUBLIC dsd_properties, autoconversion, accretion, selfcollection_breakup, &
27	evaporation_rain, sedimentation_cloud, sedimentation_rain
28
29	INTERFACE dsd_properties
30	MODULE PROCEDURE dsd_properties
31	MODULE PROCEDURE dsd_properties_ij
32	END INTERFACE dsd_properties
33
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
43
44	INTERFACE selfcollection_breakup
45	MODULE PROCEDURE selfcollection_breakup
46	MODULE PROCEDURE selfcollection_breakup_ij
47	END INTERFACE selfcollection_breakup
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
58
59	INTERFACE sedimentation_rain
60	MODULE PROCEDURE sedimentation_rain
61	MODULE PROCEDURE sedimentation_rain_ij
62	END INTERFACE sedimentation_rain
63
64	CONTAINS
65
66
67	!------------------------------------------------------------------------------!
68	! Call for all grid points
69	!------------------------------------------------------------------------------!
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
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
116	SUBROUTINE accretion
117
118	USE arrays_3d
119	USE cloud_parameters
120	USE constants
121	USE indices
122
123	IMPLICIT NONE
124
125	INTEGER :: i, j, k
126	REAL :: dqdt_precip
127
128
129	DO i = nxl, nxr
130	DO j = nys, nyn
131	DO k = nzb_2d(j,i)+1, nzt
132
133	ENDDO
134	ENDDO
135	ENDDO
136
137	END SUBROUTINE accretion
138
139	SUBROUTINE selfcollection_breakup
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
160	END SUBROUTINE selfcollection_breakup
161
162	SUBROUTINE evaporation_rain
163
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
233	!------------------------------------------------------------------------------!
234	! Call for grid point i,j
235	!------------------------------------------------------------------------------!
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
243	USE user
244
245	IMPLICIT NONE
246
247	INTEGER :: i, j, k
248
249	DO k = nzb_2d(j,i)+1, nzt
250
251	IF ( qr(k,j,i) <= eps_sb ) THEN
252	qr(k,j,i) = 0.0
253	ELSE
254	!
255	!-- Adjust number of raindrops to avoid nonlinear effects in
256	!-- sedimentation and evaporation of rain drops due to too small or
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
262	ENDIF
263	xr(k) = hyrho(k) * qr(k,j,i) / nr(k,j,i)
264	!
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 )
268	!
269	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
270	!-- Stevens and Seifert, 2008):
271	mu_r(k) = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr(k) - 1.4E-3 ) ) )
272	!
273	!-- Slope parameter of gamma distribution (Seifert, 2008):
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)
276	ENDIF
277	ENDDO
278
279	END SUBROUTINE dsd_properties_ij
280
281	SUBROUTINE autoconversion_ij( i, j )
282
283	USE arrays_3d
284	USE cloud_parameters
285	USE constants
286	USE indices
287	USE control_parameters
288	USE statistics
289	USE grid_variables
290
291	IMPLICIT NONE
292
293	INTEGER :: i, j, k
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
296
297	k_au = k_cc / ( 20.0 * x0 )
298
299	DO k = nzb_2d(j,i)+1, nzt
300
301	IF ( ql(k,j,i) > 0.0 ) THEN
302	!
303	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
304	!-- (1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) ))
305	tau_cloud = 1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) + 1.0E-20 )
306	!
307	!-- Universal function for autoconversion process
308	!-- (Seifert and Beheng, 2006):
309	phi_au = 600.0 * tau_cloud*0.68 ( 1.0 - tau_cloud0.68 )3
310	!
311	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
312	!-- (Use constant nu_c = 1.0 instead?)
313	nu_c = 1.0 !MAX( 0.0, 1580.0 * hyrho(k) * ql(k,j,i) - 0.28 )
314	!
315	!-- Mean weight of cloud droplets:
316	xc = hyrho(k) * ql(k,j,i) / nc
317	!
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	!
348	!-- Autoconversion rate (Seifert and Beheng, 2006):
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 ) &
352	rho_surface
353	autocon = MIN( autocon, ql(k,j,i) / ( dt_3d * &
354	weight_substep(intermediate_timestep_count) ) )
355	!
356	!-- Tendencies for q, qr, nr, pt:
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)
360	tend_pt(k,j,i) = tend_pt(k,j,i) + autocon * l_d_cp * pt_d_t(k)
361	ENDIF
362
363	ENDDO
364
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
374	USE statistics
375
376	IMPLICIT NONE
377
378	INTEGER :: i, j, k
379	REAL :: accr, phi_ac, tau_cloud, k_cr
380
381	DO k = nzb_2d(j,i)+1, nzt
382	IF ( ( ql(k,j,i) > 0.0 ) .AND. ( qr(k,j,i) > eps_sb ) ) THEN
383	!
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)
386	!
387	!-- Universal function for accretion process
388	!-- (Seifert and Beheng, 2001):
389	phi_ac = tau_cloud / ( tau_cloud + 5.0E-5 )
390	phi_ac = ( phi_ac2 )2
391	!
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	!
402	!-- Accretion rate (Seifert and Beheng, 2006):
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) ) )
407	!
408	!-- Tendencies for q, qr, pt:
409	tend_qr(k,j,i) = tend_qr(k,j,i) + accr
410	tend_q(k,j,i) = tend_q(k,j,i) - accr
411	tend_pt(k,j,i) = tend_pt(k,j,i) + accr * l_d_cp * pt_d_t(k)
412	ENDIF
413	ENDDO
414
415	END SUBROUTINE accretion_ij
416
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
425	USE statistics
426
427	IMPLICIT NONE
428
429	INTEGER :: i, j, k
430	REAL :: selfcoll, breakup, phi_br, phi_sc
431
432	DO k = nzb_2d(j,i)+1, nzt
433	IF ( qr(k,j,i) > eps_sb ) THEN
434	!
435	!-- Selfcollection rate (Seifert and Beheng, 2006):
436	!-- pirho_l**( 1.0 / 3.0 ) is necessary to convert [lambda_r] = m-1 to
437	!-- kg**( 1.0 / 3.0 ).
438	phi_sc = 1.0 !( 1.0 + kappa_rr / lambda_r(k) * &
439	!pirho_l( 1.0 / 3.0 ) )( -9 )
440
441	selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * phi_sc * &
442	SQRT( hyrho(k) * rho_surface )
443	!
444	!-- Collisional breakup rate (Seifert, 2008):
445	IF ( dr(k) >= 0.3E-3 ) THEN
446	phi_br = k_br * ( dr(k) - 1.1E-3 )
447	breakup = selfcoll * ( phi_br + 1.0 )
448	ELSE
449	breakup = 0.0
450	ENDIF
451
452	selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / ( dt_3d * &
453	weight_substep(intermediate_timestep_count) ) )
454	!
455	!-- Tendency for nr:
456	tend_nr(k,j,i) = tend_nr(k,j,i) + selfcoll
457	ENDIF
458	ENDDO
459
460	END SUBROUTINE selfcollection_breakup_ij
461
462	SUBROUTINE evaporation_rain_ij( i, j )
463	!
464	!-- Evaporation of precipitable water. Condensation is neglected for
465	!-- precipitable water.
466
467	USE arrays_3d
468	USE cloud_parameters
469	USE constants
470	USE indices
471	USE control_parameters
472	USE statistics
473
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, &
478	mu_r_2, mu_r_5d2, nr_0
479
480	DO k = nzb_2d(j,i)+1, nzt
481	IF ( qr(k,j,i) > eps_sb ) THEN
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:
498	temp = t_l + l_d_cp * ql(k,j,i)
499	!
500	!--
501	g_evap = ( l_v / ( r_v * temp ) - 1.0 ) * l_v / &
502	( thermal_conductivity_l * temp ) + rho_l * r_v * temp /&
503	( diff_coeff_l * e_s )
504	g_evap = 1.0 / g_evap
505	!
506	!-- Compute ventilation factor and intercept parameter
507	!-- (Seifert and Beheng, 2006; Seifert, 2008):
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 )
528	nr_0 = nr(k,j,i) * lambda_r(k)**( mu_r(k) + 1.0 ) / &
529	gamm( mu_r(k) + 1.0 )
530	ELSE
531	f_vent = 1.0
532	nr_0 = nr(k,j,i) * dr(k)
533	ENDIF
534	!
535	!-- Evaporation rate of rain water content (Seifert and Beheng, 2006):
536	evap = 2.0 * pi * nr_0 * g_evap * f_vent * sat / &
537	hyrho(k)
538	evap = MAX( evap, -qr(k,j,i) / ( dt_3d * &
539	weight_substep(intermediate_timestep_count) ) )
540	!
541	!-- Tendencies for q, qr, nr, pt:
542	tend_qr(k,j,i) = tend_qr(k,j,i) + evap
543	tend_q(k,j,i) = tend_q(k,j,i) - evap
544	tend_nr(k,j,i) = tend_nr(k,j,i) + c_evap * evap / xr(k) * hyrho(k)
545	tend_pt(k,j,i) = tend_pt(k,j,i) + evap * l_d_cp * pt_d_t(k)
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
562	REAL :: sed_q_const, sigma_gc = 1.3, k_st = 1.2E8
563	!
564	!-- Sedimentation of cloud droplets (Heus et al., 2010):
565	sed_q_const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))*( 2.0 / 3.0 ) &
566	EXP( 5.0 * LOG( sigma_gc )**2 )
567
568	sed_q = 0.0
569
570	DO k = nzb_2d(j,i)+1, nzt
571	IF ( ql(k,j,i) > 0.0 ) THEN
572	sed_q(k) = sed_q_const * nc*( -2.0 / 3.0 ) &
573	( ql(k,j,i) * hyrho(k) )**( 5.0 / 3.0 )
574	ENDIF
575	ENDDO
576	!
577	!-- Tendency for q, pt:
578	DO k = nzb_2d(j,i)+1, nzt
579	tend_q(k,j,i) = tend_q(k,j,i) + ( sed_q(k+1) - sed_q(k) ) * &
580	ddzu(k+1) / hyrho(k)
581	tend_pt(k,j,i) = tend_pt(k,j,i) - ( sed_q(k+1) - sed_q(k) ) * &
582	ddzu(k+1) / hyrho(k) * l_d_cp * pt_d_t(k)
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
594	USE statistics
595
596	IMPLICIT NONE
597
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
600
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
607	IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0
608
609	dt_sedi = dt_3d * weight_substep(intermediate_timestep_count)
610
611	w_nr = 0.0
612	w_qr = 0.0
613	!
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
626	!
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 ) )
663	ENDDO
664
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) )
686	ENDDO
687	!
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	!
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
729
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
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 )
742
743	USE cloud_parameters
744
745	IMPLICIT NONE
746
747	REAL :: gamm, ser, tmp, x_gamm, xx, y_gamm
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 |