Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |