Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |