Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |