Home

Context Navigation

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

Last change on this file since 1095 was 1093, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 25.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |