Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |