Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |