Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |