Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |