Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |