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

Line
1	MODULE microphysics_mod
2
3	!------------------------------------------------------------------------------!
4	! Current revisions:
5	! -----------------
6	!
7	! Former revisions:
8	! -----------------
9	! $Id$
10	!
11	! 1053 2012-11-13 17:11:03Z hoffmann
12	! initial revision
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
21	PUBLIC dsd_properties, autoconversion, accretion, selfcollection_breakup, &
22	evaporation_rain, sedimentation_cloud, sedimentation_rain
23
24	INTERFACE dsd_properties
25	MODULE PROCEDURE dsd_properties
26	MODULE PROCEDURE dsd_properties_ij
27	END INTERFACE dsd_properties
28
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
38
39	INTERFACE selfcollection_breakup
40	MODULE PROCEDURE selfcollection_breakup
41	MODULE PROCEDURE selfcollection_breakup_ij
42	END INTERFACE selfcollection_breakup
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
53
54	INTERFACE sedimentation_rain
55	MODULE PROCEDURE sedimentation_rain
56	MODULE PROCEDURE sedimentation_rain_ij
57	END INTERFACE sedimentation_rain
58
59	CONTAINS
60
61
62	!------------------------------------------------------------------------------!
63	! Call for all grid points
64	!------------------------------------------------------------------------------!
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
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
111	SUBROUTINE accretion
112
113	USE arrays_3d
114	USE cloud_parameters
115	USE constants
116	USE indices
117
118	IMPLICIT NONE
119
120	INTEGER :: i, j, k
121	REAL :: dqdt_precip
122
123
124	DO i = nxl, nxr
125	DO j = nys, nyn
126	DO k = nzb_2d(j,i)+1, nzt
127
128	ENDDO
129	ENDDO
130	ENDDO
131
132	END SUBROUTINE accretion
133
134	SUBROUTINE selfcollection_breakup
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
155	END SUBROUTINE selfcollection_breakup
156
157	SUBROUTINE evaporation_rain
158
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
228	!------------------------------------------------------------------------------!
229	! Call for grid point i,j
230	!------------------------------------------------------------------------------!
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
238	USE user
239
240	IMPLICIT NONE
241
242	INTEGER :: i, j, k
243	REAL :: dr_min = 2.0E-6, dr_max = 1.0E-3
244
245	DO k = nzb_2d(j,i)+1, nzt
246	IF ( ( qr(k,j,i) > eps_sb ) .AND. ( nr(k,j,i) > eps_sb ) ) THEN
247	!
248	!-- Weight averaged diameter of rain drops:
249	dr(k) = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * &
250	dpirho_l )**( 1.0 / 3.0 )
251	!
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
262	!
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)
265	!
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	!
274	!-- Slope parameter of gamma distribution (Seifert, 2008):
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)
277	ENDIF
278	ENDDO
279
280	END SUBROUTINE dsd_properties_ij
281
282	SUBROUTINE autoconversion_ij( i, j )
283
284	USE arrays_3d
285	USE cloud_parameters
286	USE constants
287	USE indices
288	USE control_parameters
289	USE statistics
290
291	IMPLICIT NONE
292
293	INTEGER :: i, j, k
294	REAL :: k_au, autocon, phi_au, tau_cloud, xc, nu_c
295
296	k_au = k_cc / ( 20.0 * x0 )
297
298	DO k = nzb_2d(j,i)+1, nzt
299
300	IF ( ql(k,j,i) > 0.0 ) THEN
301	!
302	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
303	!-- (1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) ))
304	tau_cloud = 1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) + 1.0E-20 )
305	!
306	!-- Universal function for autoconversion process
307	!-- (Seifert and Beheng, 2006):
308	phi_au = 600.0 * tau_cloud*0.68 ( 1.0 - tau_cloud0.68 )3
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:
315	xc = hyrho(k) * ql(k,j,i) / nc
316	!
317	!-- Autoconversion rate (Seifert and Beheng, 2006):
318	autocon = k_au * ( nu_c + 2.0 ) * ( nu_c + 4.0 ) / &
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) ) )
324	!
325	!-- Tendencies for q, qr, nr, pt:
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)
329	tend_pt(k,j,i) = tend_pt(k,j,i) + autocon * l_d_cp * pt_d_t(k)
330	ENDIF
331
332	ENDDO
333
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
343	USE statistics
344
345	IMPLICIT NONE
346
347	INTEGER :: i, j, k
348	REAL :: accr, phi_ac, tau_cloud
349
350	DO k = nzb_2d(j,i)+1, nzt
351	IF ( ( ql(k,j,i) > 0.0 ) .AND. ( qr(k,j,i) > eps_sb ) ) THEN
352	!
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)
355	!
356	!-- Universal function for accretion process
357	!-- (Seifert and Beheng, 2001):
358	phi_ac = tau_cloud / ( tau_cloud + 5.0E-5 )
359	phi_ac = ( phi_ac2 )2
360	!
361	!-- Accretion rate (Seifert and Beheng, 2006):
362	accr = k_cr * ql(k,j,i) * qr(k,j,i) * phi_ac * &
363	SQRT( rho_surface * hyrho(k) )
364	accr = MIN( accr, ql(k,j,i) / ( dt_3d * &
365	weight_substep(intermediate_timestep_count) ) )
366	!
367	!-- Tendencies for q, qr, pt:
368	tend_qr(k,j,i) = tend_qr(k,j,i) + accr
369	tend_q(k,j,i) = tend_q(k,j,i) - accr
370	tend_pt(k,j,i) = tend_pt(k,j,i) + accr * l_d_cp * pt_d_t(k)
371	ENDIF
372	ENDDO
373
374	END SUBROUTINE accretion_ij
375
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
384	USE statistics
385
386	IMPLICIT NONE
387
388	INTEGER :: i, j, k
389	REAL :: selfcoll, breakup, phi_br, phi_sc
390
391	DO k = nzb_2d(j,i)+1, nzt
392	IF ( ( qr(k,j,i) > eps_sb ) .AND. ( nr(k,j,i) > eps_sb ) ) THEN
393	!
394	!-- Selfcollection rate (Seifert and Beheng, 2006):
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 * &
401	SQRT( hyrho(k) * rho_surface )
402	!
403	!-- Collisional breakup rate (Seifert, 2008):
404	IF ( dr(k) >= 0.3E-3 ) THEN
405	phi_br = k_br * ( dr(k) - 1.1E-3 )
406	breakup = selfcoll * ( phi_br + 1.0 )
407	ELSE
408	breakup = 0.0
409	ENDIF
410
411	selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / ( dt_3d * &
412	weight_substep(intermediate_timestep_count) ) )
413	!
414	!-- Tendency for nr:
415	tend_nr(k,j,i) = tend_nr(k,j,i) + selfcoll
416	ENDIF
417	ENDDO
418
419	END SUBROUTINE selfcollection_breakup_ij
420
421	SUBROUTINE evaporation_rain_ij( i, j )
422	!
423	!-- Evaporation of precipitable water. Condensation is neglected for
424	!-- precipitable water.
425
426	USE arrays_3d
427	USE cloud_parameters
428	USE constants
429	USE indices
430	USE control_parameters
431	USE statistics
432
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, &
437	mu_r_2, mu_r_5d2, nr_0
438
439	DO k = nzb_2d(j,i)+1, nzt
440	IF ( ( qr(k,j,i) > eps_sb ) .AND. ( nr(k,j,i) > eps_sb ) ) THEN
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:
457	temp = t_l + l_d_cp * ql(k,j,i)
458	!
459	!--
460	g_evap = ( l_v / ( r_v * temp ) - 1.0 ) * l_v / &
461	( thermal_conductivity_l * temp ) + rho_l * r_v * temp /&
462	( diff_coeff_l * e_s )
463	g_evap = 1.0 / g_evap
464	!
465	!-- Compute ventilation factor and intercept parameter
466	!-- (Seifert and Beheng, 2006; Seifert, 2008):
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 )
487	nr_0 = nr(k,j,i) * lambda_r(k)**( mu_r(k) + 1.0 ) / &
488	gamm( mu_r(k) + 1.0 )
489	ELSE
490	f_vent = 1.0
491	nr_0 = nr(k,j,i) * dr(k)
492	ENDIF
493	!
494	!-- Evaporation rate of rain water content (Seifert and Beheng, 2006):
495	evap = 2.0 * pi * nr_0 * g_evap * f_vent * sat / &
496	hyrho(k)
497	evap = MAX( evap, -qr(k,j,i) / ( dt_3d * &
498	weight_substep(intermediate_timestep_count) ) )
499	!
500	!-- Tendencies for q, qr, nr, pt:
501	tend_qr(k,j,i) = tend_qr(k,j,i) + evap
502	tend_q(k,j,i) = tend_q(k,j,i) - evap
503	tend_nr(k,j,i) = tend_nr(k,j,i) + c_evap * evap / xr(k) * hyrho(k)
504	tend_pt(k,j,i) = tend_pt(k,j,i) + evap * l_d_cp * pt_d_t(k)
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
521	REAL :: sed_q_const, sigma_gc = 1.3, k_st = 1.2E8
522	!
523	!-- Sedimentation of cloud droplets (Heus et al., 2010):
524	sed_q_const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))*( 2.0 / 3.0 ) &
525	EXP( 5.0 * LOG( sigma_gc )**2 )
526
527	sed_q = 0.0
528
529	DO k = nzb_2d(j,i)+1, nzt
530	IF ( ql(k,j,i) > 0.0 ) THEN
531	sed_q(k) = sed_q_const * nc*( -2.0 / 3.0 ) &
532	( ql(k,j,i) * hyrho(k) )**( 5.0 / 3.0 )
533	ENDIF
534	ENDDO
535	!
536	!-- Tendency for q, pt:
537	DO k = nzb_2d(j,i)+1, nzt
538	tend_q(k,j,i) = tend_q(k,j,i) + ( sed_q(k+1) - sed_q(k) ) * &
539	ddzu(k+1) / hyrho(k)
540	tend_pt(k,j,i) = tend_pt(k,j,i) - ( sed_q(k+1) - sed_q(k) ) * &
541	ddzu(k+1) / hyrho(k) * l_d_cp * pt_d_t(k)
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
553	USE statistics
554
555	IMPLICIT NONE
556
557	INTEGER :: i, j, k, n, n_substep
558	REAL :: sed_nr_tend, sed_qr_tend
559
560	IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0
561
562	sed_nr = 0.0
563	sed_qr = 0.0
564
565	DO k = nzb_2d(j,i)+1, nzt
566	IF ( ( qr(k,j,i) > eps_sb ) .AND. ( nr(k,j,i) > eps_sb ) ) THEN
567	!
568	!-- Sedimentation of rain water content and rain drop concentration
569	!-- according to Stevens and Seifert (2008):
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)
580	ENDIF
581	ENDDO
582
583	DO k = nzb_2d(j,i)+1, nzt
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
594	ENDDO
595	!
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
600
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
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 )
613
614	USE cloud_parameters
615
616	IMPLICIT NONE
617
618	REAL :: gamm, xx, &
619	ser, tmp, x_gamm, y_gamm
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 |