Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |