Home

Context Navigation

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

Last change on this file since 1053 was 1053, checked in by hoffmann, 11 years ago
two-moment cloud physics implemented
File size: 19.7 KB

Line
1	MODULE microphysics_mod
2
3	!------------------------------------------------------------------------------!
4	! Current revisions:
5	! -----------------
6	! initial revision
7	!
8	! Former revisions:
9	! -----------------
10	! $Id$
11	!
12	! Description:
13	! ------------
14	! Calculate cloud microphysics according to the two moment bulk
15	! scheme by Seifert and Beheng (2006).
16	!------------------------------------------------------------------------------!
17
18	PRIVATE
19	PUBLIC dsd_properties, autoconversion, accretion, selfcollection_breakup, &
20	evaporation_rain, sedimentation_cloud, sedimentation_rain
21
22	INTERFACE dsd_properties
23	MODULE PROCEDURE dsd_properties
24	MODULE PROCEDURE dsd_properties_ij
25	END INTERFACE dsd_properties
26
27	INTERFACE autoconversion
28	MODULE PROCEDURE autoconversion
29	MODULE PROCEDURE autoconversion_ij
30	END INTERFACE autoconversion
31
32	INTERFACE accretion
33	MODULE PROCEDURE accretion
34	MODULE PROCEDURE accretion_ij
35	END INTERFACE accretion
36
37	INTERFACE selfcollection_breakup
38	MODULE PROCEDURE selfcollection_breakup
39	MODULE PROCEDURE selfcollection_breakup_ij
40	END INTERFACE selfcollection_breakup
41
42	INTERFACE evaporation_rain
43	MODULE PROCEDURE evaporation_rain
44	MODULE PROCEDURE evaporation_rain_ij
45	END INTERFACE evaporation_rain
46
47	INTERFACE sedimentation_cloud
48	MODULE PROCEDURE sedimentation_cloud
49	MODULE PROCEDURE sedimentation_cloud_ij
50	END INTERFACE sedimentation_cloud
51
52	INTERFACE sedimentation_rain
53	MODULE PROCEDURE sedimentation_rain
54	MODULE PROCEDURE sedimentation_rain_ij
55	END INTERFACE sedimentation_rain
56
57	CONTAINS
58
59
60	!------------------------------------------------------------------------------!
61	! Call for all grid points
62	!------------------------------------------------------------------------------!
63	SUBROUTINE dsd_properties
64
65	USE arrays_3d
66	USE cloud_parameters
67	USE constants
68	USE indices
69
70	IMPLICIT NONE
71
72	INTEGER :: i, j, k
73	REAL :: dqdt_precip
74
75
76	DO i = nxl, nxr
77	DO j = nys, nyn
78	DO k = nzb_2d(j,i)+1, nzt
79
80	ENDDO
81	ENDDO
82	ENDDO
83
84	END SUBROUTINE dsd_properties
85
86	SUBROUTINE autoconversion
87
88	USE arrays_3d
89	USE cloud_parameters
90	USE constants
91	USE indices
92
93	IMPLICIT NONE
94
95	INTEGER :: i, j, k
96	REAL :: dqdt_precip
97
98
99	DO i = nxl, nxr
100	DO j = nys, nyn
101	DO k = nzb_2d(j,i)+1, nzt
102
103	ENDDO
104	ENDDO
105	ENDDO
106
107	END SUBROUTINE autoconversion
108
109	SUBROUTINE accretion
110
111	USE arrays_3d
112	USE cloud_parameters
113	USE constants
114	USE indices
115
116	IMPLICIT NONE
117
118	INTEGER :: i, j, k
119	REAL :: dqdt_precip
120
121
122	DO i = nxl, nxr
123	DO j = nys, nyn
124	DO k = nzb_2d(j,i)+1, nzt
125
126	ENDDO
127	ENDDO
128	ENDDO
129
130	END SUBROUTINE accretion
131
132	SUBROUTINE selfcollection_breakup
133
134	USE arrays_3d
135	USE cloud_parameters
136	USE constants
137	USE indices
138
139	IMPLICIT NONE
140
141	INTEGER :: i, j, k
142	REAL :: dqdt_precip
143
144
145	DO i = nxl, nxr
146	DO j = nys, nyn
147	DO k = nzb_2d(j,i)+1, nzt
148
149	ENDDO
150	ENDDO
151	ENDDO
152
153	END SUBROUTINE selfcollection_breakup
154
155	SUBROUTINE evaporation_rain
156
157	USE arrays_3d
158	USE cloud_parameters
159	USE constants
160	USE indices
161
162	IMPLICIT NONE
163
164	INTEGER :: i, j, k
165	REAL :: dqdt_precip
166
167
168	DO i = nxl, nxr
169	DO j = nys, nyn
170	DO k = nzb_2d(j,i)+1, nzt
171
172	ENDDO
173	ENDDO
174	ENDDO
175
176	END SUBROUTINE evaporation_rain
177
178	SUBROUTINE sedimentation_cloud
179
180	USE arrays_3d
181	USE cloud_parameters
182	USE constants
183	USE indices
184
185	IMPLICIT NONE
186
187	INTEGER :: i, j, k
188	REAL :: dqdt_precip
189
190
191	DO i = nxl, nxr
192	DO j = nys, nyn
193	DO k = nzb_2d(j,i)+1, nzt
194
195	ENDDO
196	ENDDO
197	ENDDO
198
199	END SUBROUTINE sedimentation_cloud
200
201	SUBROUTINE sedimentation_rain
202
203	USE arrays_3d
204	USE cloud_parameters
205	USE constants
206	USE indices
207
208	IMPLICIT NONE
209
210	INTEGER :: i, j, k
211	REAL :: dqdt_precip
212
213
214	DO i = nxl, nxr
215	DO j = nys, nyn
216	DO k = nzb_2d(j,i)+1, nzt
217
218	ENDDO
219	ENDDO
220	ENDDO
221
222
223	END SUBROUTINE sedimentation_rain
224
225
226	!------------------------------------------------------------------------------!
227	! Call for grid point i,j
228	!------------------------------------------------------------------------------!
229	SUBROUTINE dsd_properties_ij( i, j )
230
231	USE arrays_3d
232	USE cloud_parameters
233	USE constants
234	USE indices
235	USE control_parameters
236	USE user
237
238	IMPLICIT NONE
239
240	INTEGER :: i, j, k
241	REAL :: dr_min = 2.0E-6, dr_max = 1.0E-3
242
243	DO k = nzb_2d(j,i)+1, nzt
244	IF ( ( qr(k,j,i) > eps_sb ) .AND. ( nr(k,j,i) > eps_sb ) ) THEN
245	!
246	!-- Weight averaged diameter of rain drops:
247	dr(k) = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * &
248	dpirho_l )**( 1.0 / 3.0 )
249	!
250	!-- Adjust number of raindrops to avoid nonlinear effects in
251	!-- sedimentation and evaporation of rain drops due to too small or
252	!-- big diameters of rain drops (Stevens and Seifert, 2008).
253	IF ( dr(k) < dr_min ) THEN
254	nr(k,j,i) = qr(k,j,i) * hyrho(k) / dr_min*3 dpirho_l
255	dr(k) = dr_min
256	ELSEIF ( dr(k) > dr_max ) THEN
257	nr(k,j,i) = qr(k,j,i) * hyrho(k) / dr_max*3 dpirho_l
258	dr(k) = dr_max
259	ENDIF
260	!
261	!-- Mean weight of rain drops (Seifert and Beheng, 2006):
262	xr(k) = MIN( MAX( hyrho(k) * qr(k,j,i) / nr(k,j,i), xrmin ), xrmax)
263	!
264	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
265	!-- Stevens and Seifert, 2008):
266	IF ( .NOT. mu_constant ) THEN
267	mu_r(k) = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr(k) - 1.4E-3 ) ) )
268	ELSE
269	mu_r(k) = mu_constant_value
270	ENDIF
271	!
272	!-- Slope parameter of gamma distribution (Seifert, 2008):
273	lambda_r(k) = ( ( mu_r(k) + 3.0 ) * ( mu_r(k) + 2.0 ) * &
274	( mu_r(k) + 1.0 ) )**( 1.0 / 3.0 ) / dr(k)
275	ENDIF
276	ENDDO
277
278	END SUBROUTINE dsd_properties_ij
279
280	SUBROUTINE autoconversion_ij( i, j )
281
282	USE arrays_3d
283	USE cloud_parameters
284	USE constants
285	USE indices
286	USE control_parameters
287	USE statistics
288
289	IMPLICIT NONE
290
291	INTEGER :: i, j, k
292	REAL :: k_au, autocon, phi_au, tau_cloud, xc, nu_c
293
294	k_au = k_cc / ( 20.0 * x0 )
295
296	DO k = nzb_2d(j,i)+1, nzt
297
298	IF ( ql(k,j,i) > 0.0 ) THEN
299	!
300	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
301	!-- (1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) ))
302	tau_cloud = 1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) + 1.0E-20 )
303	!
304	!-- Universal function for autoconversion process
305	!-- (Seifert and Beheng, 2006):
306	phi_au = 600.0 * tau_cloud*0.68 ( 1.0 - tau_cloud0.68 )3
307	!
308	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
309	!-- (Use constant nu_c = 1.0 instead?)
310	nu_c = 1580.0 * hyrho(k) * ql(k,j,i) - 0.28
311	!
312	!-- Mean weight of cloud droplets:
313	xc = hyrho(k) * ql(k,j,i) / nc
314	!
315	!-- Autoconversion rate (Seifert and Beheng, 2006):
316	autocon = k_au * ( nu_c + 2.0 ) * ( nu_c + 4.0 ) / &
317	( nu_c + 1.0 )*2 ql(k,j,i)*2 xc*2 &
318	( 1.0 + phi_au / ( 1.0 - tau_cloud )*2 ) &
319	rho_surface
320	autocon = MIN( autocon, ql(k,j,i) / ( dt_3d * &
321	weight_substep(intermediate_timestep_count) ) )
322	!
323	!-- Tendencies for q, qr, nr, pt:
324	tend_qr(k,j,i) = tend_qr(k,j,i) + autocon
325	tend_q(k,j,i) = tend_q(k,j,i) - autocon
326	tend_nr(k,j,i) = tend_nr(k,j,i) + autocon / x0 * hyrho(k)
327	tend_pt(k,j,i) = tend_pt(k,j,i) + autocon * l_d_cp * pt_d_t(k)
328	ENDIF
329
330	ENDDO
331
332	END SUBROUTINE autoconversion_ij
333
334	SUBROUTINE accretion_ij( i, j )
335
336	USE arrays_3d
337	USE cloud_parameters
338	USE constants
339	USE indices
340	USE control_parameters
341	USE statistics
342
343	IMPLICIT NONE
344
345	INTEGER :: i, j, k
346	REAL :: accr, phi_ac, tau_cloud
347
348	DO k = nzb_2d(j,i)+1, nzt
349	IF ( ( ql(k,j,i) > 0.0 ) .AND. ( qr(k,j,i) > eps_sb ) ) THEN
350	!
351	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
352	tau_cloud = 1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) + 1.0E-20)
353	!
354	!-- Universal function for accretion process
355	!-- (Seifert and Beheng, 2001):
356	phi_ac = tau_cloud / ( tau_cloud + 5.0E-5 )
357	phi_ac = ( phi_ac2 )2
358	!
359	!-- Accretion rate (Seifert and Beheng, 2006):
360	accr = k_cr * ql(k,j,i) * qr(k,j,i) * phi_ac * &
361	SQRT( rho_surface * hyrho(k) )
362	accr = MIN( accr, ql(k,j,i) / ( dt_3d * &
363	weight_substep(intermediate_timestep_count) ) )
364	!
365	!-- Tendencies for q, qr, pt:
366	tend_qr(k,j,i) = tend_qr(k,j,i) + accr
367	tend_q(k,j,i) = tend_q(k,j,i) - accr
368	tend_pt(k,j,i) = tend_pt(k,j,i) + accr * l_d_cp * pt_d_t(k)
369	ENDIF
370	ENDDO
371
372	END SUBROUTINE accretion_ij
373
374
375	SUBROUTINE selfcollection_breakup_ij( i, j )
376
377	USE arrays_3d
378	USE cloud_parameters
379	USE constants
380	USE indices
381	USE control_parameters
382	USE statistics
383
384	IMPLICIT NONE
385
386	INTEGER :: i, j, k
387	REAL :: selfcoll, breakup, phi_br, phi_sc
388
389	DO k = nzb_2d(j,i)+1, nzt
390	IF ( ( qr(k,j,i) > eps_sb ) .AND. ( nr(k,j,i) > eps_sb ) ) THEN
391	!
392	!-- Selfcollection rate (Seifert and Beheng, 2006):
393	!-- pirho_l**( 1.0 / 3.0 ) is necessary to convert [lambda_r] = m-1 to
394	!-- kg**( 1.0 / 3.0 ).
395	phi_sc = ( 1.0 + kappa_rr / lambda_r(k) * &
396	pirho_l( 1.0 / 3.0 ) )( -9 )
397
398	selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * phi_sc * &
399	SQRT( hyrho(k) * rho_surface )
400	!
401	!-- Collisional breakup rate (Seifert, 2008):
402	IF ( dr(k) >= 0.3E-3 ) THEN
403	phi_br = k_br * ( dr(k) - 1.1E-3 )
404	breakup = selfcoll * ( phi_br + 1.0 )
405	ELSE
406	breakup = 0.0
407	ENDIF
408
409	selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / ( dt_3d * &
410	weight_substep(intermediate_timestep_count) ) )
411	!
412	!-- Tendency for nr:
413	tend_nr(k,j,i) = tend_nr(k,j,i) + selfcoll
414	ENDIF
415	ENDDO
416
417	END SUBROUTINE selfcollection_breakup_ij
418
419	SUBROUTINE evaporation_rain_ij( i, j )
420	!
421	!-- Evaporation of precipitable water. Condensation is neglected for
422	!-- precipitable water.
423
424	USE arrays_3d
425	USE cloud_parameters
426	USE constants
427	USE indices
428	USE control_parameters
429	USE statistics
430
431	IMPLICIT NONE
432
433	INTEGER :: i, j, k
434	REAL :: evap, alpha, e_s, q_s, t_l, sat, temp, g_evap, f_vent, &
435	mu_r_2, mu_r_5d2, nr_0
436
437	DO k = nzb_2d(j,i)+1, nzt
438	IF ( ( qr(k,j,i) > eps_sb ) .AND. ( nr(k,j,i) > eps_sb ) ) THEN
439	!
440	!-- Actual liquid water temperature:
441	t_l = t_d_pt(k) * pt(k,j,i)
442	!
443	!-- Saturation vapor pressure at t_l:
444	e_s = 610.78 * EXP( 17.269 * ( t_l - 273.16 ) / ( t_l - 35.86 ) )
445	!
446	!-- Computation of saturation humidity:
447	q_s = 0.622 * e_s / ( hyp(k) - 0.378 * e_s )
448	alpha = 0.622 * l_d_r * l_d_cp / ( t_l * t_l )
449	q_s = q_s * ( 1.0 + alpha * q(k,j,i) ) / ( 1.0 + alpha * q_s )
450	!
451	!-- Oversaturation:
452	sat = MIN( 0.0, ( q(k,j,i) - ql(k,j,i) ) / q_s - 1.0 )
453	!
454	!-- Actual temperature:
455	temp = t_l + l_d_cp * ql(k,j,i)
456	!
457	!--
458	g_evap = ( l_v / ( r_v * temp ) - 1.0 ) * l_v / &
459	( thermal_conductivity_l * temp ) + rho_l * r_v * temp /&
460	( diff_coeff_l * e_s )
461	g_evap = 1.0 / g_evap
462	!
463	!-- Compute ventilation factor and intercept parameter
464	!-- (Seifert and Beheng, 2006; Seifert, 2008):
465	IF ( ventilation_effect ) THEN
466	mu_r_2 = mu_r(k) + 2.0
467	mu_r_5d2 = mu_r(k) + 2.5
468	f_vent = a_vent * gamm( mu_r_2 ) * &
469	lambda_r(k)**( -mu_r_2 ) + &
470	b_vent * schmidt_p_1d3 * &
471	SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * &
472	lambda_r(k)*( -mu_r_5d2 ) &
473	( 1.0 - 0.5 * ( b_term / a_term ) * &
474	( lambda_r(k) / &
475	( c_term + lambda_r(k) ) )**mu_r_5d2 - &
476	0.125 * ( b_term / a_term )*2 &
477	( lambda_r(k) / &
478	( 2.0 * c_term + lambda_r(k) ) )**mu_r_5d2 - &
479	0.0625 * ( b_term / a_term )*3 &
480	( lambda_r(k) / &
481	( 3.0 * c_term + lambda_r(k) ) )**mu_r_5d2 - &
482	0.0390625 * ( b_term / a_term )*4 &
483	( lambda_r(k) / &
484	( 4.0 * c_term + lambda_r(k) ) )**mu_r_5d2 )
485	nr_0 = nr(k,j,i) * lambda_r(k)**( mu_r(k) + 1.0 ) / &
486	gamm( mu_r(k) + 1.0 )
487	ELSE
488	f_vent = 1.0
489	nr_0 = nr(k,j,i) * dr(k)
490	ENDIF
491	!
492	!-- Evaporation rate of rain water content (Seifert and Beheng, 2006):
493	evap = 2.0 * pi * nr_0 * g_evap * f_vent * sat / &
494	hyrho(k)
495	evap = MAX( evap, -qr(k,j,i) / ( dt_3d * &
496	weight_substep(intermediate_timestep_count) ) )
497	!
498	!-- Tendencies for q, qr, nr, pt:
499	tend_qr(k,j,i) = tend_qr(k,j,i) + evap
500	tend_q(k,j,i) = tend_q(k,j,i) - evap
501	tend_nr(k,j,i) = tend_nr(k,j,i) + c_evap * evap / xr(k) * hyrho(k)
502	tend_pt(k,j,i) = tend_pt(k,j,i) + evap * l_d_cp * pt_d_t(k)
503	ENDIF
504	ENDDO
505
506	END SUBROUTINE evaporation_rain_ij
507
508	SUBROUTINE sedimentation_cloud_ij( i, j )
509
510	USE arrays_3d
511	USE cloud_parameters
512	USE constants
513	USE indices
514	USE control_parameters
515
516	IMPLICIT NONE
517
518	INTEGER :: i, j, k
519	REAL :: sed_q_const, sigma_gc = 1.3, k_st = 1.2E8
520	!
521	!-- Sedimentation of cloud droplets (Heus et al., 2010):
522	sed_q_const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))*( 2.0 / 3.0 ) &
523	EXP( 5.0 * LOG( sigma_gc )**2 )
524
525	sed_q = 0.0
526
527	DO k = nzb_2d(j,i)+1, nzt
528	IF ( ql(k,j,i) > 0.0 ) THEN
529	sed_q(k) = sed_q_const * nc*( -2.0 / 3.0 ) &
530	( ql(k,j,i) * hyrho(k) )**( 5.0 / 3.0 )
531	ENDIF
532	ENDDO
533	!
534	!-- Tendency for q, pt:
535	DO k = nzb_2d(j,i)+1, nzt
536	tend_q(k,j,i) = tend_q(k,j,i) + ( sed_q(k+1) - sed_q(k) ) * &
537	ddzu(k+1) / hyrho(k)
538	tend_pt(k,j,i) = tend_pt(k,j,i) - ( sed_q(k+1) - sed_q(k) ) * &
539	ddzu(k+1) / hyrho(k) * l_d_cp * pt_d_t(k)
540	ENDDO
541
542	END SUBROUTINE sedimentation_cloud_ij
543
544	SUBROUTINE sedimentation_rain_ij( i, j )
545
546	USE arrays_3d
547	USE cloud_parameters
548	USE constants
549	USE indices
550	USE control_parameters
551	USE statistics
552
553	IMPLICIT NONE
554
555	INTEGER :: i, j, k, n, n_substep
556	REAL :: sed_nr_tend, sed_qr_tend
557
558	IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0
559
560	sed_nr = 0.0
561	sed_qr = 0.0
562
563	DO k = nzb_2d(j,i)+1, nzt
564	IF ( ( qr(k,j,i) > eps_sb ) .AND. ( nr(k,j,i) > eps_sb ) ) THEN
565	!
566	!-- Sedimentation of rain water content and rain drop concentration
567	!-- according to Stevens and Seifert (2008):
568	sed_nr(k) = MIN( 20.0, MAX( 0.1, a_term - b_term * ( 1.0 + &
569	c_term / lambda_r(k) )*( -1.0 &
570	( mu_r(k) + 1.0 ) ) ) ) * nr(k,j,i)
571	sed_qr(k) = MIN( 20.0, MAX( 0.1, a_term - b_term * ( 1.0 + &
572	c_term / lambda_r(k) )*( -1.0 &
573	( mu_r(k) + 4.0 ) ) ) ) * qr(k,j,i) * hyrho(k)
574	!
575	!-- Computation of rain rate
576	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * &
577	weight_substep(intermediate_timestep_count)
578	ENDIF
579	ENDDO
580
581	DO k = nzb_2d(j,i)+1, nzt
582	sed_nr_tend = MAX( ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1), &
583	-nr(k,j,i) / ( dt_3d * &
584	weight_substep(intermediate_timestep_count) ) )
585	sed_qr_tend = MAX( ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
586	hyrho(k), &
587	-qr(k,j,i) / ( dt_3d * &
588	weight_substep(intermediate_timestep_count) ) )
589
590	tend_nr(k,j,i) = tend_nr(k,j,i) + sed_nr_tend
591	tend_qr(k,j,i) = tend_qr(k,j,i) + sed_qr_tend
592	ENDDO
593	!
594	!-- Precipitation amount
595	IF ( intermediate_timestep_count == intermediate_timestep_count_max &
596	.AND. ( dt_do2d_xy - time_do2d_xy ) < &
597	precipitation_amount_interval ) THEN
598
599	precipitation_amount(j,i) = precipitation_amount(j,i) + &
600	prr(nzb_2d(j,i)+1,j,i) * &
601	hyrho(nzb_2d(j,i)+1) * dt_3d
602	ENDIF
603
604	END SUBROUTINE sedimentation_rain_ij
605
606	!
607	!-- This function computes the gamma function (Press et al., 1992).
608	!-- The gamma function is needed for the calculation of the evaporation
609	!-- of rain drops.
610	FUNCTION gamm( xx )
611
612	USE cloud_parameters
613
614	IMPLICIT NONE
615
616	REAL :: gamm, xx, &
617	ser, tmp, x_gamm, y_gamm
618	INTEGER :: j
619
620	x_gamm = xx
621	y_gamm = x_gamm
622	tmp = x_gamm + 5.5
623	tmp = ( x_gamm + 0.5 ) * LOG( tmp ) - tmp
624	ser = 1.000000000190015
625	do j = 1, 6
626	y_gamm = y_gamm + 1.0
627	ser = ser + cof( j ) / y_gamm
628	enddo
629	!
630	!-- Until this point the algorithm computes the logarithm of the gamma
631	!-- function. Hence, the exponential function is used.
632	! gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
633	gamm = EXP( tmp ) * stp * ser / x_gamm
634	RETURN
635
636	END FUNCTION gamm
637
638	END MODULE microphysics_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |