Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |