Home

Context Navigation

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

Last change on this file since 1395 was 1362, checked in by hoffmann, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 72.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	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: microphysics.f90 1362 2014-04-16 15:19:12Z keck $
27	!
28	! 1361 2014-04-16 15:17:48Z hoffmann
29	! Bugfix in sedimentation_rain: Index corrected.
30	! Vectorized version of adjust_cloud added.
31	! Little reformatting of the code.
32	!
33	! 1353 2014-04-08 15:21:23Z heinze
34	! REAL constants provided with KIND-attribute
35	!
36	! 1346 2014-03-27 13:18:20Z heinze
37	! Bugfix: REAL constants provided with KIND-attribute especially in call of
38	! intrinsic function like MAX, MIN, SIGN
39	!
40	! 1334 2014-03-25 12:21:40Z heinze
41	! Bugfix: REAL constants provided with KIND-attribute
42	!
43	! 1322 2014-03-20 16:38:49Z raasch
44	! REAL constants defined as wp-kind
45	!
46	! 1320 2014-03-20 08:40:49Z raasch
47	! ONLY-attribute added to USE-statements,
48	! kind-parameters added to all INTEGER and REAL declaration statements,
49	! kinds are defined in new module kinds,
50	! comment fields (!:) to be used for variable explanations added to
51	! all variable declaration statements
52	!
53	! 1241 2013-10-30 11:36:58Z heinze
54	! hyp and rho have to be calculated at each time step if data from external
55	! file LSF_DATA are used
56	!
57	! 1115 2013-03-26 18:16:16Z hoffmann
58	! microphyical tendencies are calculated in microphysics_control in an optimized
59	! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
60	!
61	! 1106 2013-03-04 05:31:38Z raasch
62	! small changes in code formatting
63	!
64	! 1092 2013-02-02 11:24:22Z raasch
65	! unused variables removed
66	! file put under GPL
67	!
68	! 1065 2012-11-22 17:42:36Z hoffmann
69	! Sedimentation process implemented according to Stevens and Seifert (2008).
70	! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
71	! and Stevens, 2010).
72	!
73	! 1053 2012-11-13 17:11:03Z hoffmann
74	! initial revision
75	!
76	! Description:
77	! ------------
78	! Calculate cloud microphysics according to the two moment bulk
79	! scheme by Seifert and Beheng (2006).
80	!------------------------------------------------------------------------------!
81
82	PRIVATE
83	PUBLIC microphysics_control
84
85	INTERFACE microphysics_control
86	MODULE PROCEDURE microphysics_control
87	MODULE PROCEDURE microphysics_control_ij
88	END INTERFACE microphysics_control
89
90	INTERFACE adjust_cloud
91	MODULE PROCEDURE adjust_cloud
92	MODULE PROCEDURE adjust_cloud_ij
93	END INTERFACE adjust_cloud
94
95	INTERFACE autoconversion
96	MODULE PROCEDURE autoconversion
97	MODULE PROCEDURE autoconversion_ij
98	END INTERFACE autoconversion
99
100	INTERFACE accretion
101	MODULE PROCEDURE accretion
102	MODULE PROCEDURE accretion_ij
103	END INTERFACE accretion
104
105	INTERFACE selfcollection_breakup
106	MODULE PROCEDURE selfcollection_breakup
107	MODULE PROCEDURE selfcollection_breakup_ij
108	END INTERFACE selfcollection_breakup
109
110	INTERFACE evaporation_rain
111	MODULE PROCEDURE evaporation_rain
112	MODULE PROCEDURE evaporation_rain_ij
113	END INTERFACE evaporation_rain
114
115	INTERFACE sedimentation_cloud
116	MODULE PROCEDURE sedimentation_cloud
117	MODULE PROCEDURE sedimentation_cloud_ij
118	END INTERFACE sedimentation_cloud
119
120	INTERFACE sedimentation_rain
121	MODULE PROCEDURE sedimentation_rain
122	MODULE PROCEDURE sedimentation_rain_ij
123	END INTERFACE sedimentation_rain
124
125	CONTAINS
126
127
128	!------------------------------------------------------------------------------!
129	! Call for all grid points
130	!------------------------------------------------------------------------------!
131	SUBROUTINE microphysics_control
132
133	USE arrays_3d, &
134	ONLY: hyp, nr, pt, pt_init, q, qc, qr, zu
135
136	USE cloud_parameters, &
137	ONLY: cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt
138
139	USE control_parameters, &
140	ONLY: call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, &
141	g, intermediate_timestep_count, &
142	large_scale_forcing, lsf_surf, precipitation, pt_surface, &
143	rho_surface,surface_pressure
144
145	USE indices, &
146	ONLY: nzb, nzt
147
148	USE kinds
149
150	USE statistics, &
151	ONLY: weight_pres
152
153	IMPLICIT NONE
154
155	INTEGER(iwp) :: i !:
156	INTEGER(iwp) :: j !:
157	INTEGER(iwp) :: k !:
158
159	REAL(wp) :: t_surface !:
160
161	IF ( large_scale_forcing .AND. lsf_surf ) THEN
162	!
163	!-- Calculate:
164	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
165	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
166	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
167	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
168	DO k = nzb, nzt+1
169	hyp(k) = surface_pressure * 100.0_wp * &
170	( ( t_surface - g / cp * zu(k) ) / &
171	t_surface )**(1.0_wp / 0.286_wp)
172	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
173	t_d_pt(k) = 1.0_wp / pt_d_t(k)
174	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
175	ENDDO
176	!
177	!-- Compute reference density
178	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
179	ENDIF
180
181	!
182	!-- Compute length of time step
183	IF ( call_microphysics_at_all_substeps ) THEN
184	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
185	ELSE
186	dt_micro = dt_3d
187	ENDIF
188
189	!
190	!-- Compute cloud physics
191	IF ( precipitation ) THEN
192	CALL adjust_cloud
193	CALL autoconversion
194	CALL accretion
195	CALL selfcollection_breakup
196	CALL evaporation_rain
197	CALL sedimentation_rain
198	ENDIF
199
200	IF ( drizzle ) CALL sedimentation_cloud
201
202	END SUBROUTINE microphysics_control
203
204	SUBROUTINE adjust_cloud
205
206	USE arrays_3d, &
207	ONLY: qr, nr
208
209	USE cloud_parameters, &
210	ONLY: eps_sb, xrmin, xrmax, hyrho, k_cc, x0
211
212	USE cpulog, &
213	ONLY: cpu_log, log_point_s
214
215	USE indices, &
216	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
217
218	USE kinds
219
220	IMPLICIT NONE
221
222	INTEGER(iwp) :: i !:
223	INTEGER(iwp) :: j !:
224	INTEGER(iwp) :: k !:
225
226	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
227
228	DO i = nxl, nxr
229	DO j = nys, nyn
230	DO k = nzb_s_inner(j,i)+1, nzt
231	IF ( qr(k,j,i) <= eps_sb ) THEN
232	qr(k,j,i) = 0.0_wp
233	nr(k,j,i) = 0.0_wp
234	ELSE
235	!
236	!-- Adjust number of raindrops to avoid nonlinear effects in
237	!-- sedimentation and evaporation of rain drops due to too small
238	!-- or too big weights of rain drops (Stevens and Seifert, 2008).
239	IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN
240	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
241	ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN
242	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
243	ENDIF
244
245	ENDIF
246	ENDDO
247	ENDDO
248	ENDDO
249
250	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
251
252	END SUBROUTINE adjust_cloud
253
254
255	SUBROUTINE autoconversion
256
257	USE arrays_3d, &
258	ONLY: diss, dzu, nr, qc, qr
259
260	USE cloud_parameters, &
261	ONLY: a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3, &
262	c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air, &
263	nc_const, x0
264
265	USE control_parameters, &
266	ONLY: dt_micro, rho_surface, turbulence
267
268	USE cpulog, &
269	ONLY: cpu_log, log_point_s
270
271	USE grid_variables, &
272	ONLY: dx, dy
273
274	USE indices, &
275	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
276
277	USE kinds
278
279	IMPLICIT NONE
280
281	INTEGER(iwp) :: i !:
282	INTEGER(iwp) :: j !:
283	INTEGER(iwp) :: k !:
284
285	REAL(wp) :: alpha_cc !:
286	REAL(wp) :: autocon !:
287	REAL(wp) :: dissipation !:
288	REAL(wp) :: k_au !:
289	REAL(wp) :: l_mix !:
290	REAL(wp) :: nu_c !:
291	REAL(wp) :: phi_au !:
292	REAL(wp) :: r_cc !:
293	REAL(wp) :: rc !:
294	REAL(wp) :: re_lambda !:
295	REAL(wp) :: selfcoll !:
296	REAL(wp) :: sigma_cc !:
297	REAL(wp) :: tau_cloud !:
298	REAL(wp) :: xc !:
299
300	CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
301
302	DO i = nxl, nxr
303	DO j = nys, nyn
304	DO k = nzb_s_inner(j,i)+1, nzt
305
306	IF ( qc(k,j,i) > eps_sb ) THEN
307
308	k_au = k_cc / ( 20.0_wp * x0 )
309	!
310	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
311	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
312	tau_cloud = 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) )
313	!
314	!-- Universal function for autoconversion process
315	!-- (Seifert and Beheng, 2006):
316	phi_au = 600.0_wp * tau_cloud*0.68_wp &
317	( 1.0_wp - tau_cloud0.68_wp )3
318	!
319	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
320	!-- (Use constant nu_c = 1.0_wp instead?)
321	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
322	!
323	!-- Mean weight of cloud droplets:
324	xc = hyrho(k) * qc(k,j,i) / nc_const
325	!
326	!-- Parameterized turbulence effects on autoconversion (Seifert,
327	!-- Nuijens and Stevens, 2010)
328	IF ( turbulence ) THEN
329	!
330	!-- Weight averaged radius of cloud droplets:
331	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
332
333	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
334	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
335	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
336	!
337	!-- Mixing length (neglecting distance to ground and
338	!-- stratification)
339	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
340	!
341	!-- Limit dissipation rate according to Seifert, Nuijens and
342	!-- Stevens (2010)
343	dissipation = MIN( 0.06_wp, diss(k,j,i) )
344	!
345	!-- Compute Taylor-microscale Reynolds number:
346	re_lambda = 6.0_wp / 11.0_wp * &
347	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
348	SQRT( 15.0_wp / kin_vis_air ) * &
349	dissipation**( 1.0_wp / 6.0_wp )
350	!
351	!-- The factor of 1.0E4 is needed to convert the dissipation
352	!-- rate from m2 s-3 to cm2 s-3.
353	k_au = k_au * ( 1.0_wp + &
354	dissipation * 1.0E4_wp * &
355	( re_lambda * 1.0E-3_wp )*0.25_wp &
356	( alpha_cc * EXP( -1.0_wp * ( ( rc - &
357	r_cc ) / &
358	sigma_cc )**2 &
359	) + beta_cc &
360	) &
361	)
362	ENDIF
363	!
364	!-- Autoconversion rate (Seifert and Beheng, 2006):
365	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
366	( nu_c + 1.0_wp )*2 qc(k,j,i)*2 xc*2 &
367	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
368	rho_surface
369	autocon = MIN( autocon, qc(k,j,i) / dt_micro )
370
371	qr(k,j,i) = qr(k,j,i) + autocon * dt_micro
372	qc(k,j,i) = qc(k,j,i) - autocon * dt_micro
373	nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro
374
375	ENDIF
376
377	ENDDO
378	ENDDO
379	ENDDO
380
381	CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
382
383	END SUBROUTINE autoconversion
384
385
386	SUBROUTINE accretion
387
388	USE arrays_3d, &
389	ONLY: diss, qc, qr
390
391	USE cloud_parameters, &
392	ONLY: eps_sb, hyrho, k_cr0
393
394	USE control_parameters, &
395	ONLY: dt_micro, rho_surface, turbulence
396
397	USE cpulog, &
398	ONLY: cpu_log, log_point_s
399
400	USE indices, &
401	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
402
403	USE kinds
404
405	IMPLICIT NONE
406
407	INTEGER(iwp) :: i !:
408	INTEGER(iwp) :: j !:
409	INTEGER(iwp) :: k !:
410
411	REAL(wp) :: accr !:
412	REAL(wp) :: k_cr !:
413	REAL(wp) :: phi_ac !:
414	REAL(wp) :: tau_cloud !:
415	REAL(wp) :: xc !:
416
417	CALL cpu_log( log_point_s(56), 'accretion', 'start' )
418
419	DO i = nxl, nxr
420	DO j = nys, nyn
421	DO k = nzb_s_inner(j,i)+1, nzt
422
423	IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) ) THEN
424	!
425	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
426	tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
427	!
428	!-- Universal function for accretion process (Seifert and
429	!-- Beheng, 2001):
430	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
431	!
432	!-- Parameterized turbulence effects on autoconversion (Seifert,
433	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
434	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
435	IF ( turbulence ) THEN
436	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
437	MIN( 600.0_wp, &
438	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
439	)
440	ELSE
441	k_cr = k_cr0
442	ENDIF
443	!
444	!-- Accretion rate (Seifert and Beheng, 2006):
445	accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * &
446	SQRT( rho_surface * hyrho(k) )
447	accr = MIN( accr, qc(k,j,i) / dt_micro )
448
449	qr(k,j,i) = qr(k,j,i) + accr * dt_micro
450	qc(k,j,i) = qc(k,j,i) - accr * dt_micro
451
452	ENDIF
453
454	ENDDO
455	ENDDO
456	ENDDO
457
458	CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
459
460	END SUBROUTINE accretion
461
462
463	SUBROUTINE selfcollection_breakup
464
465	USE arrays_3d, &
466	ONLY: nr, qr
467
468	USE cloud_parameters, &
469	ONLY: dpirho_l, eps_sb, hyrho, k_br, k_rr
470
471	USE control_parameters, &
472	ONLY: dt_micro, rho_surface
473
474	USE cpulog, &
475	ONLY: cpu_log, log_point_s
476
477	USE indices, &
478	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
479
480	USE kinds
481
482	IMPLICIT NONE
483
484	INTEGER(iwp) :: i !:
485	INTEGER(iwp) :: j !:
486	INTEGER(iwp) :: k !:
487
488	REAL(wp) :: breakup !:
489	REAL(wp) :: dr !:
490	REAL(wp) :: phi_br !:
491	REAL(wp) :: selfcoll !:
492
493	CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
494
495	DO i = nxl, nxr
496	DO j = nys, nyn
497	DO k = nzb_s_inner(j,i)+1, nzt
498	IF ( qr(k,j,i) > eps_sb ) THEN
499	!
500	!-- Selfcollection rate (Seifert and Beheng, 2001):
501	selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * &
502	SQRT( hyrho(k) * rho_surface )
503	!
504	!-- Weight averaged diameter of rain drops:
505	dr = ( hyrho(k) * qr(k,j,i) / &
506	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
507	!
508	!-- Collisional breakup rate (Seifert, 2008):
509	IF ( dr >= 0.3E-3_wp ) THEN
510	phi_br = k_br * ( dr - 1.1E-3_wp )
511	breakup = selfcoll * ( phi_br + 1.0_wp )
512	ELSE
513	breakup = 0.0_wp
514	ENDIF
515
516	selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
517	nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro
518
519	ENDIF
520	ENDDO
521	ENDDO
522	ENDDO
523
524	CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
525
526	END SUBROUTINE selfcollection_breakup
527
528
529	SUBROUTINE evaporation_rain
530
531	!
532	!-- Evaporation of precipitable water. Condensation is neglected for
533	!-- precipitable water.
534
535	USE arrays_3d, &
536	ONLY: hyp, nr, pt, q, qc, qr
537
538	USE cloud_parameters, &
539	ONLY: a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,&
540	dpirho_l, eps_sb, hyrho, kin_vis_air, k_st, l_d_cp, l_d_r, &
541	l_v, rho_l, r_v, schmidt_p_1d3, thermal_conductivity_l, &
542	t_d_pt, ventilation_effect
543
544	USE constants, &
545	ONLY: pi
546
547	USE control_parameters, &
548	ONLY: dt_micro
549
550	USE cpulog, &
551	ONLY: cpu_log, log_point_s
552
553	USE indices, &
554	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
555
556	USE kinds
557
558	IMPLICIT NONE
559
560	INTEGER(iwp) :: i !:
561	INTEGER(iwp) :: j !:
562	INTEGER(iwp) :: k !:
563
564	REAL(wp) :: alpha !:
565	REAL(wp) :: dr !:
566	REAL(wp) :: e_s !:
567	REAL(wp) :: evap !:
568	REAL(wp) :: evap_nr !:
569	REAL(wp) :: f_vent !:
570	REAL(wp) :: g_evap !:
571	REAL(wp) :: lambda_r !:
572	REAL(wp) :: mu_r !:
573	REAL(wp) :: mu_r_2 !:
574	REAL(wp) :: mu_r_5d2 !:
575	REAL(wp) :: nr_0 !:
576	REAL(wp) :: q_s !:
577	REAL(wp) :: sat !:
578	REAL(wp) :: t_l !:
579	REAL(wp) :: temp !:
580	REAL(wp) :: xr !:
581
582	CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
583
584	DO i = nxl, nxr
585	DO j = nys, nyn
586	DO k = nzb_s_inner(j,i)+1, nzt
587	IF ( qr(k,j,i) > eps_sb ) THEN
588	!
589	!-- Actual liquid water temperature:
590	t_l = t_d_pt(k) * pt(k,j,i)
591	!
592	!-- Saturation vapor pressure at t_l:
593	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
594	( t_l - 35.86_wp ) &
595	)
596	!
597	!-- Computation of saturation humidity:
598	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
599	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
600	q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / &
601	( 1.0_wp + alpha * q_s )
602	!
603	!-- Supersaturation:
604	sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
605	!
606	!-- Evaporation needs only to be calculated in subsaturated regions
607	IF ( sat < 0.0_wp ) THEN
608	!
609	!-- Actual temperature:
610	temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
611
612	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
613	l_v / ( thermal_conductivity_l * temp ) &
614	+ r_v * temp / ( diff_coeff_l * e_s ) &
615	)
616	!
617	!-- Mean weight of rain drops
618	xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
619	!
620	!-- Weight averaged diameter of rain drops:
621	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
622	!
623	!-- Compute ventilation factor and intercept parameter
624	!-- (Seifert and Beheng, 2006; Seifert, 2008):
625	IF ( ventilation_effect ) THEN
626	!
627	!-- Shape parameter of gamma distribution (Milbrandt and Yau,
628	!-- 2005; Stevens and Seifert, 2008):
629	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
630	( dr - 1.4E-3_wp ) ) )
631	!
632	!-- Slope parameter of gamma distribution (Seifert, 2008):
633	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
634	( mu_r + 1.0_wp ) &
635	)**( 1.0_wp / 3.0_wp ) / dr
636
637	mu_r_2 = mu_r + 2.0_wp
638	mu_r_5d2 = mu_r + 2.5_wp
639
640	f_vent = a_vent * gamm( mu_r_2 ) * &
641	lambda_r*( -mu_r_2 ) + b_vent &
642	schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
643	gamm( mu_r_5d2 ) * lambda_r*( -mu_r_5d2 ) &
644	( 1.0_wp - &
645	0.5_wp * ( b_term / a_term ) * &
646	( lambda_r / ( c_term + lambda_r ) &
647	)**mu_r_5d2 - &
648	0.125_wp * ( b_term / a_term )*2 &
649	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
650	)**mu_r_5d2 - &
651	0.0625_wp * ( b_term / a_term )*3 &
652	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
653	)**mu_r_5d2 - &
654	0.0390625_wp * ( b_term / a_term )*4 &
655	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
656	)**mu_r_5d2 &
657	)
658
659	nr_0 = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) / &
660	gamm( mu_r + 1.0_wp )
661	ELSE
662	f_vent = 1.0_wp
663	nr_0 = nr(k,j,i) * dr
664	ENDIF
665	!
666	!-- Evaporation rate of rain water content (Seifert and
667	!-- Beheng, 2006):
668	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / &
669	hyrho(k)
670	evap = MAX( evap, -qr(k,j,i) / dt_micro )
671	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
672	-nr(k,j,i) / dt_micro )
673
674	qr(k,j,i) = qr(k,j,i) + evap * dt_micro
675	nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro
676
677	ENDIF
678	ENDIF
679
680	ENDDO
681	ENDDO
682	ENDDO
683
684	CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
685
686	END SUBROUTINE evaporation_rain
687
688
689	SUBROUTINE sedimentation_cloud
690
691	USE arrays_3d, &
692	ONLY: ddzu, dzu, pt, q, qc
693
694	USE cloud_parameters, &
695	ONLY: eps_sb, hyrho, l_d_cp, nc_const, pt_d_t, sed_qc_const
696
697	USE constants, &
698	ONLY: pi
699
700	USE control_parameters, &
701	ONLY: dt_do2d_xy, dt_micro, intermediate_timestep_count
702
703	USE cpulog, &
704	ONLY: cpu_log, log_point_s
705
706	USE indices, &
707	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
708
709	USE kinds
710
711	IMPLICIT NONE
712
713	INTEGER(iwp) :: i !:
714	INTEGER(iwp) :: j !:
715	INTEGER(iwp) :: k !:
716
717	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !:
718
719	CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
720
721	!
722	!-- Sedimentation of cloud droplets (Ackermann et al., 2009, MWR):
723	sed_qc(nzt+1) = 0.0_wp
724
725	DO i = nxl, nxr
726	DO j = nys, nyn
727	DO k = nzt, nzb_s_inner(j,i)+1, -1
728
729	IF ( qc(k,j,i) > eps_sb ) THEN
730	sed_qc(k) = sed_qc_const * nc_const*( -2.0_wp / 3.0_wp ) &
731	( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp )
732	ELSE
733	sed_qc(k) = 0.0_wp
734	ENDIF
735
736	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / &
737	dt_micro + sed_qc(k+1) &
738	)
739
740	q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
741	ddzu(k+1) / hyrho(k) * dt_micro
742	qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
743	ddzu(k+1) / hyrho(k) * dt_micro
744	pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * &
745	ddzu(k+1) / hyrho(k) * l_d_cp * &
746	pt_d_t(k) * dt_micro
747
748	ENDDO
749	ENDDO
750	ENDDO
751
752	CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
753
754	END SUBROUTINE sedimentation_cloud
755
756
757	SUBROUTINE sedimentation_rain
758
759	USE arrays_3d, &
760	ONLY: ddzu, dzu, nr, pt, q, qr
761
762	USE cloud_parameters, &
763	ONLY: a_term, b_term, c_term, cof, dpirho_l, eps_sb, hyrho, &
764	limiter_sedimentation, l_d_cp, precipitation_amount, prr, &
765	pt_d_t, stp
766
767	USE control_parameters, &
768	ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro, &
769	dt_3d, intermediate_timestep_count, &
770	intermediate_timestep_count_max, &
771	precipitation_amount_interval, time_do2d_xy
772
773	USE cpulog, &
774	ONLY: cpu_log, log_point_s
775
776	USE indices, &
777	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
778
779	USE kinds
780
781	USE statistics, &
782	ONLY: weight_substep
783
784	IMPLICIT NONE
785
786	INTEGER(iwp) :: i !:
787	INTEGER(iwp) :: j !:
788	INTEGER(iwp) :: k !:
789	INTEGER(iwp) :: k_run !:
790
791	REAL(wp) :: c_run !:
792	REAL(wp) :: d_max !:
793	REAL(wp) :: d_mean !:
794	REAL(wp) :: d_min !:
795	REAL(wp) :: dr !:
796	REAL(wp) :: dt_sedi !:
797	REAL(wp) :: flux !:
798	REAL(wp) :: lambda_r !:
799	REAL(wp) :: mu_r !:
800	REAL(wp) :: z_run !:
801
802	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !:
803	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !:
804	REAL(wp), DIMENSION(nzb:nzt+1) :: d_nr !:
805	REAL(wp), DIMENSION(nzb:nzt+1) :: d_qr !:
806	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !:
807	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !:
808	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !:
809	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !:
810	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !:
811	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !:
812
813	CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
814	!
815	!-- Computation of sedimentation flux. Implementation according to Stevens
816	!-- and Seifert (2008). Code is based on UCLA-LES.
817	IF ( intermediate_timestep_count == 1 ) prr(:,:,:) = 0.0_wp
818	!
819	!-- Compute velocities
820	DO i = nxl, nxr
821	DO j = nys, nyn
822	DO k = nzb_s_inner(j,i)+1, nzt
823	IF ( qr(k,j,i) > eps_sb ) THEN
824	!
825	!-- Weight averaged diameter of rain drops:
826	dr = ( hyrho(k) * qr(k,j,i) / &
827	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
828	!
829	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
830	!-- Stevens and Seifert, 2008):
831	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
832	( dr - 1.4E-3_wp ) ) )
833	!
834	!-- Slope parameter of gamma distribution (Seifert, 2008):
835	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
836	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
837
838	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
839	a_term - b_term * ( 1.0_wp + &
840	c_term / &
841	lambda_r )*( -1.0_wp &
842	( mu_r + 1.0_wp ) ) &
843	) &
844	)
845
846	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
847	a_term - b_term * ( 1.0_wp + &
848	c_term / &
849	lambda_r )*( -1.0_wp &
850	( mu_r + 4.0_wp ) ) &
851	) &
852	)
853	ELSE
854	w_nr(k) = 0.0_wp
855	w_qr(k) = 0.0_wp
856	ENDIF
857	ENDDO
858	!
859	!-- Adjust boundary values
860	w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
861	w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
862	w_nr(nzt+1) = 0.0_wp
863	w_qr(nzt+1) = 0.0_wp
864	!
865	!-- Compute Courant number
866	DO k = nzb_s_inner(j,i)+1, nzt
867	c_nr(k) = 0.25_wp * ( w_nr(k-1) + &
868	2.0_wp * w_nr(k) + w_nr(k+1) ) * &
869	dt_micro * ddzu(k)
870	c_qr(k) = 0.25_wp * ( w_qr(k-1) + &
871	2.0_wp * w_qr(k) + w_qr(k+1) ) * &
872	dt_micro * ddzu(k)
873	ENDDO
874	!
875	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
876	IF ( limiter_sedimentation ) THEN
877
878	DO k = nzb_s_inner(j,i)+1, nzt
879	d_mean = 0.5_wp * ( qr(k+1,j,i) + qr(k-1,j,i) )
880	d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
881	d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
882
883	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
884	2.0_wp * d_max, &
885	ABS( d_mean ) )
886
887	d_mean = 0.5_wp * ( nr(k+1,j,i) + nr(k-1,j,i) )
888	d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
889	d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
890
891	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
892	2.0_wp * d_max, &
893	ABS( d_mean ) )
894	ENDDO
895
896	ELSE
897
898	nr_slope = 0.0_wp
899	qr_slope = 0.0_wp
900
901	ENDIF
902
903	sed_nr(nzt+1) = 0.0_wp
904	sed_qr(nzt+1) = 0.0_wp
905	!
906	!-- Compute sedimentation flux
907	DO k = nzt, nzb_s_inner(j,i)+1, -1
908	!
909	!-- Sum up all rain drop number densities which contribute to the flux
910	!-- through k-1/2
911	flux = 0.0_wp
912	z_run = 0.0_wp ! height above z(k)
913	k_run = k
914	c_run = MIN( 1.0_wp, c_nr(k) )
915	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
916	flux = flux + hyrho(k_run) * &
917	( nr(k_run,j,i) + nr_slope(k_run) * &
918	( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)
919	z_run = z_run + dzu(k_run)
920	k_run = k_run + 1
921	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
922	ENDDO
923	!
924	!-- It is not allowed to sediment more rain drop number density than
925	!-- available
926	flux = MIN( flux, &
927	hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * &
928	dt_micro &
929	)
930
931	sed_nr(k) = flux / dt_micro
932	nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * &
933	ddzu(k+1) / hyrho(k) * dt_micro
934	!
935	!-- Sum up all rain water content which contributes to the flux
936	!-- through k-1/2
937	flux = 0.0_wp
938	z_run = 0.0_wp ! height above z(k)
939	k_run = k
940	c_run = MIN( 1.0_wp, c_qr(k) )
941
942	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
943
944	flux = flux + hyrho(k_run) * ( qr(k_run,j,i) + &
945	qr_slope(k_run) * ( 1.0_wp - c_run ) * &
946	0.5_wp ) * c_run * dzu(k_run)
947	z_run = z_run + dzu(k_run)
948	k_run = k_run + 1
949	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
950
951	ENDDO
952	!
953	!-- It is not allowed to sediment more rain water content than
954	!-- available
955	flux = MIN( flux, &
956	hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * &
957	dt_micro &
958	)
959
960	sed_qr(k) = flux / dt_micro
961
962	qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
963	ddzu(k+1) / hyrho(k) * dt_micro
964	q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
965	ddzu(k+1) / hyrho(k) * dt_micro
966	pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * &
967	ddzu(k+1) / hyrho(k) * l_d_cp * &
968	pt_d_t(k) * dt_micro
969	!
970	!-- Compute the rain rate
971	IF ( call_microphysics_at_all_substeps ) THEN
972	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * &
973	weight_substep(intermediate_timestep_count)
974	ELSE
975	prr(k,j,i) = sed_qr(k) / hyrho(k)
976	ENDIF
977
978	ENDDO
979	ENDDO
980	ENDDO
981
982	!
983	!-- Precipitation amount
984	IF ( intermediate_timestep_count == intermediate_timestep_count_max &
985	.AND. ( dt_do2d_xy - time_do2d_xy ) < &
986	precipitation_amount_interval ) THEN
987	DO i = nxl, nxr
988	DO j = nys, nyn
989	precipitation_amount(j,i) = precipitation_amount(j,i) + &
990	prr(nzb_s_inner(j,i)+1,j,i) * &
991	hyrho(nzb_s_inner(j,i)+1) * dt_3d
992	ENDDO
993	ENDDO
994	ENDIF
995
996	CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
997
998	END SUBROUTINE sedimentation_rain
999
1000
1001	!------------------------------------------------------------------------------!
1002	! Call for grid point i,j
1003	!------------------------------------------------------------------------------!
1004
1005	SUBROUTINE microphysics_control_ij( i, j )
1006
1007	USE arrays_3d, &
1008	ONLY: hyp, nc_1d, nr, nr_1d, pt, pt_init, pt_1d, q, q_1d, qc, &
1009	qc_1d, qr, qr_1d, zu
1010
1011	USE cloud_parameters, &
1012	ONLY: cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt
1013
1014	USE control_parameters, &
1015	ONLY: call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, &
1016	g, intermediate_timestep_count, large_scale_forcing, &
1017	lsf_surf, precipitation, pt_surface, &
1018	rho_surface,surface_pressure
1019
1020	USE indices, &
1021	ONLY: nzb, nzt
1022
1023	USE kinds
1024
1025	USE statistics, &
1026	ONLY: weight_pres
1027
1028	IMPLICIT NONE
1029
1030	INTEGER(iwp) :: i !:
1031	INTEGER(iwp) :: j !:
1032	INTEGER(iwp) :: k !:
1033
1034	REAL(wp) :: t_surface !:
1035
1036	IF ( large_scale_forcing .AND. lsf_surf ) THEN
1037	!
1038	!-- Calculate:
1039	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
1040	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
1041	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
1042	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
1043	DO k = nzb, nzt+1
1044	hyp(k) = surface_pressure * 100.0_wp * &
1045	( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
1046	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
1047	t_d_pt(k) = 1.0_wp / pt_d_t(k)
1048	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
1049	ENDDO
1050	!
1051	!-- Compute reference density
1052	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
1053	ENDIF
1054
1055	!
1056	!-- Compute length of time step
1057	IF ( call_microphysics_at_all_substeps ) THEN
1058	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
1059	ELSE
1060	dt_micro = dt_3d
1061	ENDIF
1062
1063	!
1064	!-- Use 1d arrays
1065	q_1d(:) = q(:,j,i)
1066	pt_1d(:) = pt(:,j,i)
1067	qc_1d(:) = qc(:,j,i)
1068	nc_1d(:) = nc_const
1069	IF ( precipitation ) THEN
1070	qr_1d(:) = qr(:,j,i)
1071	nr_1d(:) = nr(:,j,i)
1072	ENDIF
1073
1074	!
1075	!-- Compute cloud physics
1076	IF ( precipitation ) THEN
1077	CALL adjust_cloud( i,j )
1078	CALL autoconversion( i,j )
1079	CALL accretion( i,j )
1080	CALL selfcollection_breakup( i,j )
1081	CALL evaporation_rain( i,j )
1082	CALL sedimentation_rain( i,j )
1083	ENDIF
1084
1085	IF ( drizzle ) CALL sedimentation_cloud( i,j )
1086
1087	!
1088	!-- Store results on the 3d arrays
1089	q(:,j,i) = q_1d(:)
1090	pt(:,j,i) = pt_1d(:)
1091	IF ( precipitation ) THEN
1092	qr(:,j,i) = qr_1d(:)
1093	nr(:,j,i) = nr_1d(:)
1094	ENDIF
1095
1096	END SUBROUTINE microphysics_control_ij
1097
1098	SUBROUTINE adjust_cloud_ij( i, j )
1099
1100	USE arrays_3d, &
1101	ONLY: qr_1d, nr_1d
1102
1103	USE cloud_parameters, &
1104	ONLY: eps_sb, xrmin, xrmax, hyrho, k_cc, x0
1105
1106	USE indices, &
1107	ONLY: nzb, nzb_s_inner, nzt
1108
1109	USE kinds
1110
1111	IMPLICIT NONE
1112
1113	INTEGER(iwp) :: i !:
1114	INTEGER(iwp) :: j !:
1115	INTEGER(iwp) :: k !:
1116	!
1117	!-- Adjust number of raindrops to avoid nonlinear effects in
1118	!-- sedimentation and evaporation of rain drops due to too small or
1119	!-- too big weights of rain drops (Stevens and Seifert, 2008).
1120	!-- The same procedure is applied to cloud droplets if they are determined
1121	!-- prognostically.
1122	DO k = nzb_s_inner(j,i)+1, nzt
1123
1124	IF ( qr_1d(k) <= eps_sb ) THEN
1125	qr_1d(k) = 0.0_wp
1126	nr_1d(k) = 0.0_wp
1127	ELSE
1128	!
1129	!-- Adjust number of raindrops to avoid nonlinear effects in
1130	!-- sedimentation and evaporation of rain drops due to too small or
1131	!-- too big weights of rain drops (Stevens and Seifert, 2008).
1132	IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) ) THEN
1133	nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin
1134	ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) ) THEN
1135	nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax
1136	ENDIF
1137
1138	ENDIF
1139
1140	ENDDO
1141
1142	END SUBROUTINE adjust_cloud_ij
1143
1144
1145	SUBROUTINE autoconversion_ij( i, j )
1146
1147	USE arrays_3d, &
1148	ONLY: diss, dzu, nc_1d, nr_1d, qc_1d, qr_1d
1149
1150	USE cloud_parameters, &
1151	ONLY: a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3, &
1152	c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air, x0
1153
1154	USE control_parameters, &
1155	ONLY: dt_micro, rho_surface, turbulence
1156
1157	USE grid_variables, &
1158	ONLY: dx, dy
1159
1160	USE indices, &
1161	ONLY: nzb, nzb_s_inner, nzt
1162
1163	USE kinds
1164
1165	IMPLICIT NONE
1166
1167	INTEGER(iwp) :: i !:
1168	INTEGER(iwp) :: j !:
1169	INTEGER(iwp) :: k !:
1170
1171	REAL(wp) :: alpha_cc !:
1172	REAL(wp) :: autocon !:
1173	REAL(wp) :: dissipation !:
1174	REAL(wp) :: k_au !:
1175	REAL(wp) :: l_mix !:
1176	REAL(wp) :: nu_c !:
1177	REAL(wp) :: phi_au !:
1178	REAL(wp) :: r_cc !:
1179	REAL(wp) :: rc !:
1180	REAL(wp) :: re_lambda !:
1181	REAL(wp) :: selfcoll !:
1182	REAL(wp) :: sigma_cc !:
1183	REAL(wp) :: tau_cloud !:
1184	REAL(wp) :: xc !:
1185
1186	DO k = nzb_s_inner(j,i)+1, nzt
1187
1188	IF ( qc_1d(k) > eps_sb ) THEN
1189
1190	k_au = k_cc / ( 20.0_wp * x0 )
1191	!
1192	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
1193	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
1194	tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
1195	!
1196	!-- Universal function for autoconversion process
1197	!-- (Seifert and Beheng, 2006):
1198	phi_au = 600.0_wp * tau_cloud*0.68_wp ( 1.0_wp - tau_cloud0.68_wp )3
1199	!
1200	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
1201	!-- (Use constant nu_c = 1.0_wp instead?)
1202	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp )
1203	!
1204	!-- Mean weight of cloud droplets:
1205	xc = hyrho(k) * qc_1d(k) / nc_1d(k)
1206	!
1207	!-- Parameterized turbulence effects on autoconversion (Seifert,
1208	!-- Nuijens and Stevens, 2010)
1209	IF ( turbulence ) THEN
1210	!
1211	!-- Weight averaged radius of cloud droplets:
1212	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
1213
1214	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
1215	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
1216	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
1217	!
1218	!-- Mixing length (neglecting distance to ground and stratification)
1219	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
1220	!
1221	!-- Limit dissipation rate according to Seifert, Nuijens and
1222	!-- Stevens (2010)
1223	dissipation = MIN( 0.06_wp, diss(k,j,i) )
1224	!
1225	!-- Compute Taylor-microscale Reynolds number:
1226	re_lambda = 6.0_wp / 11.0_wp * &
1227	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
1228	SQRT( 15.0_wp / kin_vis_air ) * &
1229	dissipation**( 1.0_wp / 6.0_wp )
1230	!
1231	!-- The factor of 1.0E4 is needed to convert the dissipation rate
1232	!-- from m2 s-3 to cm2 s-3.
1233	k_au = k_au * ( 1.0_wp + &
1234	dissipation * 1.0E4_wp * &
1235	( re_lambda * 1.0E-3_wp )*0.25_wp &
1236	( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) / &
1237	sigma_cc )**2 &
1238	) + beta_cc &
1239	) &
1240	)
1241	ENDIF
1242	!
1243	!-- Autoconversion rate (Seifert and Beheng, 2006):
1244	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
1245	( nu_c + 1.0_wp )*2 qc_1d(k)*2 xc*2 &
1246	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
1247	rho_surface
1248	autocon = MIN( autocon, qc_1d(k) / dt_micro )
1249
1250	qr_1d(k) = qr_1d(k) + autocon * dt_micro
1251	qc_1d(k) = qc_1d(k) - autocon * dt_micro
1252	nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
1253
1254	ENDIF
1255
1256	ENDDO
1257
1258	END SUBROUTINE autoconversion_ij
1259
1260
1261	SUBROUTINE accretion_ij( i, j )
1262
1263	USE arrays_3d, &
1264	ONLY: diss, qc_1d, qr_1d
1265
1266	USE cloud_parameters, &
1267	ONLY: eps_sb, hyrho, k_cr0
1268
1269	USE control_parameters, &
1270	ONLY: dt_micro, rho_surface, turbulence
1271
1272	USE indices, &
1273	ONLY: nzb, nzb_s_inner, nzt
1274
1275	USE kinds
1276
1277	IMPLICIT NONE
1278
1279	INTEGER(iwp) :: i !:
1280	INTEGER(iwp) :: j !:
1281	INTEGER(iwp) :: k !:
1282
1283	REAL(wp) :: accr !:
1284	REAL(wp) :: k_cr !:
1285	REAL(wp) :: phi_ac !:
1286	REAL(wp) :: tau_cloud !:
1287	REAL(wp) :: xc !:
1288
1289	DO k = nzb_s_inner(j,i)+1, nzt
1290	IF ( ( qc_1d(k) > eps_sb ) .AND. ( qr_1d(k) > eps_sb ) ) THEN
1291	!
1292	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
1293	tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
1294	!
1295	!-- Universal function for accretion process
1296	!-- (Seifert and Beheng, 2001):
1297	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
1298	!
1299	!-- Parameterized turbulence effects on autoconversion (Seifert,
1300	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
1301	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
1302	IF ( turbulence ) THEN
1303	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
1304	MIN( 600.0_wp, &
1305	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
1306	)
1307	ELSE
1308	k_cr = k_cr0
1309	ENDIF
1310	!
1311	!-- Accretion rate (Seifert and Beheng, 2006):
1312	accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * SQRT( rho_surface * hyrho(k) )
1313	accr = MIN( accr, qc_1d(k) / dt_micro )
1314
1315	qr_1d(k) = qr_1d(k) + accr * dt_micro
1316	qc_1d(k) = qc_1d(k) - accr * dt_micro
1317
1318	ENDIF
1319
1320	ENDDO
1321
1322	END SUBROUTINE accretion_ij
1323
1324
1325	SUBROUTINE selfcollection_breakup_ij( i, j )
1326
1327	USE arrays_3d, &
1328	ONLY: nr_1d, qr_1d
1329
1330	USE cloud_parameters, &
1331	ONLY: dpirho_l, eps_sb, hyrho, k_br, k_rr
1332
1333	USE control_parameters, &
1334	ONLY: dt_micro, rho_surface
1335
1336	USE indices, &
1337	ONLY: nzb, nzb_s_inner, nzt
1338
1339	USE kinds
1340
1341	IMPLICIT NONE
1342
1343	INTEGER(iwp) :: i !:
1344	INTEGER(iwp) :: j !:
1345	INTEGER(iwp) :: k !:
1346
1347	REAL(wp) :: breakup !:
1348	REAL(wp) :: dr !:
1349	REAL(wp) :: phi_br !:
1350	REAL(wp) :: selfcoll !:
1351
1352	DO k = nzb_s_inner(j,i)+1, nzt
1353	IF ( qr_1d(k) > eps_sb ) THEN
1354	!
1355	!-- Selfcollection rate (Seifert and Beheng, 2001):
1356	selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
1357	!
1358	!-- Weight averaged diameter of rain drops:
1359	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
1360	!
1361	!-- Collisional breakup rate (Seifert, 2008):
1362	IF ( dr >= 0.3E-3_wp ) THEN
1363	phi_br = k_br * ( dr - 1.1E-3_wp )
1364	breakup = selfcoll * ( phi_br + 1.0_wp )
1365	ELSE
1366	breakup = 0.0_wp
1367	ENDIF
1368
1369	selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
1370	nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
1371
1372	ENDIF
1373	ENDDO
1374
1375	END SUBROUTINE selfcollection_breakup_ij
1376
1377
1378	SUBROUTINE evaporation_rain_ij( i, j )
1379	!
1380	!-- Evaporation of precipitable water. Condensation is neglected for
1381	!-- precipitable water.
1382
1383	USE arrays_3d, &
1384	ONLY: hyp, nr_1d, pt_1d, q_1d, qc_1d, qr_1d
1385
1386	USE cloud_parameters, &
1387	ONLY: a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,&
1388	dpirho_l, eps_sb, hyrho, kin_vis_air, k_st, l_d_cp, l_d_r, &
1389	l_v, rho_l, r_v, schmidt_p_1d3, thermal_conductivity_l, &
1390	t_d_pt, ventilation_effect
1391
1392	USE constants, &
1393	ONLY: pi
1394
1395	USE control_parameters, &
1396	ONLY: dt_micro
1397
1398	USE indices, &
1399	ONLY: nzb, nzb_s_inner, nzt
1400
1401	USE kinds
1402
1403	IMPLICIT NONE
1404
1405	INTEGER(iwp) :: i !:
1406	INTEGER(iwp) :: j !:
1407	INTEGER(iwp) :: k !:
1408
1409	REAL(wp) :: alpha !:
1410	REAL(wp) :: dr !:
1411	REAL(wp) :: e_s !:
1412	REAL(wp) :: evap !:
1413	REAL(wp) :: evap_nr !:
1414	REAL(wp) :: f_vent !:
1415	REAL(wp) :: g_evap !:
1416	REAL(wp) :: lambda_r !:
1417	REAL(wp) :: mu_r !:
1418	REAL(wp) :: mu_r_2 !:
1419	REAL(wp) :: mu_r_5d2 !:
1420	REAL(wp) :: nr_0 !:
1421	REAL(wp) :: q_s !:
1422	REAL(wp) :: sat !:
1423	REAL(wp) :: t_l !:
1424	REAL(wp) :: temp !:
1425	REAL(wp) :: xr !:
1426
1427	DO k = nzb_s_inner(j,i)+1, nzt
1428	IF ( qr_1d(k) > eps_sb ) THEN
1429	!
1430	!-- Actual liquid water temperature:
1431	t_l = t_d_pt(k) * pt_1d(k)
1432	!
1433	!-- Saturation vapor pressure at t_l:
1434	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
1435	( t_l - 35.86_wp ) &
1436	)
1437	!
1438	!-- Computation of saturation humidity:
1439	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
1440	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
1441	q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
1442	!
1443	!-- Supersaturation:
1444	sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
1445	!
1446	!-- Evaporation needs only to be calculated in subsaturated regions
1447	IF ( sat < 0.0_wp ) THEN
1448	!
1449	!-- Actual temperature:
1450	temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
1451
1452	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v / &
1453	( thermal_conductivity_l * temp ) + &
1454	r_v * temp / ( diff_coeff_l * e_s ) &
1455	)
1456	!
1457	!-- Mean weight of rain drops
1458	xr = hyrho(k) * qr_1d(k) / nr_1d(k)
1459	!
1460	!-- Weight averaged diameter of rain drops:
1461	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
1462	!
1463	!-- Compute ventilation factor and intercept parameter
1464	!-- (Seifert and Beheng, 2006; Seifert, 2008):
1465	IF ( ventilation_effect ) THEN
1466	!
1467	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1468	!-- Stevens and Seifert, 2008):
1469	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
1470	!
1471	!-- Slope parameter of gamma distribution (Seifert, 2008):
1472	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
1473	( mu_r + 1.0_wp ) &
1474	)**( 1.0_wp / 3.0_wp ) / dr
1475
1476	mu_r_2 = mu_r + 2.0_wp
1477	mu_r_5d2 = mu_r + 2.5_wp
1478
1479	f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) + &
1480	b_vent * schmidt_p_1d3 * &
1481	SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * &
1482	lambda_r*( -mu_r_5d2 ) &
1483	( 1.0_wp - &
1484	0.5_wp * ( b_term / a_term ) * &
1485	( lambda_r / ( c_term + lambda_r ) &
1486	)**mu_r_5d2 - &
1487	0.125_wp * ( b_term / a_term )*2 &
1488	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
1489	)**mu_r_5d2 - &
1490	0.0625_wp * ( b_term / a_term )*3 &
1491	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
1492	)**mu_r_5d2 - &
1493	0.0390625_wp * ( b_term / a_term )*4 &
1494	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
1495	)**mu_r_5d2 &
1496	)
1497
1498	nr_0 = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) / &
1499	gamm( mu_r + 1.0_wp )
1500	ELSE
1501	f_vent = 1.0_wp
1502	nr_0 = nr_1d(k) * dr
1503	ENDIF
1504	!
1505	!-- Evaporation rate of rain water content (Seifert and Beheng, 2006):
1506	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
1507	evap = MAX( evap, -qr_1d(k) / dt_micro )
1508	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
1509	-nr_1d(k) / dt_micro )
1510
1511	qr_1d(k) = qr_1d(k) + evap * dt_micro
1512	nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
1513
1514	ENDIF
1515	ENDIF
1516
1517	ENDDO
1518
1519	END SUBROUTINE evaporation_rain_ij
1520
1521
1522	SUBROUTINE sedimentation_cloud_ij( i, j )
1523
1524	USE arrays_3d, &
1525	ONLY: ddzu, dzu, nc_1d, pt_1d, q_1d, qc_1d
1526
1527	USE cloud_parameters, &
1528	ONLY: eps_sb, hyrho, l_d_cp, pt_d_t, sed_qc_const
1529
1530	USE constants, &
1531	ONLY: pi
1532
1533	USE control_parameters, &
1534	ONLY: dt_do2d_xy, dt_micro, intermediate_timestep_count
1535
1536	USE indices, &
1537	ONLY: nzb, nzb_s_inner, nzt
1538
1539	USE kinds
1540
1541	IMPLICIT NONE
1542
1543	INTEGER(iwp) :: i !:
1544	INTEGER(iwp) :: j !:
1545	INTEGER(iwp) :: k !:
1546
1547	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !:
1548
1549	!
1550	!-- Sedimentation of cloud droplets (Ackermann et al., 2009, MWR):
1551	sed_qc(nzt+1) = 0.0_wp
1552
1553	DO k = nzt, nzb_s_inner(j,i)+1, -1
1554	IF ( qc_1d(k) > eps_sb ) THEN
1555	sed_qc(k) = sed_qc_const * nc_1d(k)*( -2.0_wp / 3.0_wp ) &
1556	( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )
1557	ELSE
1558	sed_qc(k) = 0.0_wp
1559	ENDIF
1560
1561	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) / &
1562	dt_micro + sed_qc(k+1) &
1563	)
1564
1565	q_1d(k) = q_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
1566	hyrho(k) * dt_micro
1567	qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
1568	hyrho(k) * dt_micro
1569	pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
1570	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
1571
1572	ENDDO
1573
1574	END SUBROUTINE sedimentation_cloud_ij
1575
1576
1577	SUBROUTINE sedimentation_rain_ij( i, j )
1578
1579	USE arrays_3d, &
1580	ONLY: ddzu, dzu, nr_1d, pt_1d, q_1d, qr_1d
1581
1582	USE cloud_parameters, &
1583	ONLY: a_term, b_term, c_term, cof, dpirho_l, eps_sb, hyrho, &
1584	limiter_sedimentation, l_d_cp, precipitation_amount, prr, &
1585	pt_d_t, stp
1586
1587	USE control_parameters, &
1588	ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro, &
1589	dt_3d, intermediate_timestep_count, &
1590	intermediate_timestep_count_max, &
1591	precipitation_amount_interval, time_do2d_xy
1592
1593	USE indices, &
1594	ONLY: nzb, nzb_s_inner, nzt
1595
1596	USE kinds
1597
1598	USE statistics, &
1599	ONLY: weight_substep
1600
1601	IMPLICIT NONE
1602
1603	INTEGER(iwp) :: i !:
1604	INTEGER(iwp) :: j !:
1605	INTEGER(iwp) :: k !:
1606	INTEGER(iwp) :: k_run !:
1607
1608	REAL(wp) :: c_run !:
1609	REAL(wp) :: d_max !:
1610	REAL(wp) :: d_mean !:
1611	REAL(wp) :: d_min !:
1612	REAL(wp) :: dr !:
1613	REAL(wp) :: dt_sedi !:
1614	REAL(wp) :: flux !:
1615	REAL(wp) :: lambda_r !:
1616	REAL(wp) :: mu_r !:
1617	REAL(wp) :: z_run !:
1618
1619	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !:
1620	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !:
1621	REAL(wp), DIMENSION(nzb:nzt+1) :: d_nr !:
1622	REAL(wp), DIMENSION(nzb:nzt+1) :: d_qr !:
1623	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !:
1624	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !:
1625	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !:
1626	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !:
1627	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !:
1628	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !:
1629
1630
1631	!
1632	!-- Computation of sedimentation flux. Implementation according to Stevens
1633	!-- and Seifert (2008). Code is based on UCLA-LES.
1634	IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp
1635	!
1636	!-- Compute velocities
1637	DO k = nzb_s_inner(j,i)+1, nzt
1638	IF ( qr_1d(k) > eps_sb ) THEN
1639	!
1640	!-- Weight averaged diameter of rain drops:
1641	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
1642	!
1643	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1644	!-- Stevens and Seifert, 2008):
1645	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
1646	!
1647	!-- Slope parameter of gamma distribution (Seifert, 2008):
1648	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
1649	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
1650
1651	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
1652	a_term - b_term * ( 1.0_wp + &
1653	c_term / lambda_r )*( -1.0_wp &
1654	( mu_r + 1.0_wp ) ) &
1655	) &
1656	)
1657	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
1658	a_term - b_term * ( 1.0_wp + &
1659	c_term / lambda_r )*( -1.0_wp &
1660	( mu_r + 4.0_wp ) ) &
1661	) &
1662	)
1663	ELSE
1664	w_nr(k) = 0.0_wp
1665	w_qr(k) = 0.0_wp
1666	ENDIF
1667	ENDDO
1668	!
1669	!-- Adjust boundary values
1670	w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
1671	w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
1672	w_nr(nzt+1) = 0.0_wp
1673	w_qr(nzt+1) = 0.0_wp
1674	!
1675	!-- Compute Courant number
1676	DO k = nzb_s_inner(j,i)+1, nzt
1677	c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) * &
1678	dt_micro * ddzu(k)
1679	c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * &
1680	dt_micro * ddzu(k)
1681	ENDDO
1682	!
1683	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
1684	IF ( limiter_sedimentation ) THEN
1685
1686	DO k = nzb_s_inner(j,i)+1, nzt
1687	d_mean = 0.5_wp * ( qr_1d(k+1) + qr_1d(k-1) )
1688	d_min = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
1689	d_max = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
1690
1691	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
1692	2.0_wp * d_max, &
1693	ABS( d_mean ) )
1694
1695	d_mean = 0.5_wp * ( nr_1d(k+1) + nr_1d(k-1) )
1696	d_min = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
1697	d_max = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
1698
1699	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
1700	2.0_wp * d_max, &
1701	ABS( d_mean ) )
1702	ENDDO
1703
1704	ELSE
1705
1706	nr_slope = 0.0_wp
1707	qr_slope = 0.0_wp
1708
1709	ENDIF
1710
1711	sed_nr(nzt+1) = 0.0_wp
1712	sed_qr(nzt+1) = 0.0_wp
1713	!
1714	!-- Compute sedimentation flux
1715	DO k = nzt, nzb_s_inner(j,i)+1, -1
1716	!
1717	!-- Sum up all rain drop number densities which contribute to the flux
1718	!-- through k-1/2
1719	flux = 0.0_wp
1720	z_run = 0.0_wp ! height above z(k)
1721	k_run = k
1722	c_run = MIN( 1.0_wp, c_nr(k) )
1723	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
1724	flux = flux + hyrho(k_run) * &
1725	( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) * &
1726	0.5_wp ) * c_run * dzu(k_run)
1727	z_run = z_run + dzu(k_run)
1728	k_run = k_run + 1
1729	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
1730	ENDDO
1731	!
1732	!-- It is not allowed to sediment more rain drop number density than
1733	!-- available
1734	flux = MIN( flux, &
1735	hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
1736
1737	sed_nr(k) = flux / dt_micro
1738	nr_1d(k) = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / &
1739	hyrho(k) * dt_micro
1740	!
1741	!-- Sum up all rain water content which contributes to the flux
1742	!-- through k-1/2
1743	flux = 0.0_wp
1744	z_run = 0.0_wp ! height above z(k)
1745	k_run = k
1746	c_run = MIN( 1.0_wp, c_qr(k) )
1747
1748	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
1749
1750	flux = flux + hyrho(k_run) * &
1751	( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) * &
1752	0.5_wp ) * c_run * dzu(k_run)
1753	z_run = z_run + dzu(k_run)
1754	k_run = k_run + 1
1755	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
1756
1757	ENDDO
1758	!
1759	!-- It is not allowed to sediment more rain water content than available
1760	flux = MIN( flux, &
1761	hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
1762
1763	sed_qr(k) = flux / dt_micro
1764
1765	qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
1766	hyrho(k) * dt_micro
1767	q_1d(k) = q_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
1768	hyrho(k) * dt_micro
1769	pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
1770	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
1771	!
1772	!-- Compute the rain rate
1773	IF ( call_microphysics_at_all_substeps ) THEN
1774	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * &
1775	weight_substep(intermediate_timestep_count)
1776	ELSE
1777	prr(k,j,i) = sed_qr(k) / hyrho(k)
1778	ENDIF
1779
1780	ENDDO
1781
1782	!
1783	!-- Precipitation amount
1784	IF ( intermediate_timestep_count == intermediate_timestep_count_max &
1785	.AND. ( dt_do2d_xy - time_do2d_xy ) < &
1786	precipitation_amount_interval ) THEN
1787
1788	precipitation_amount(j,i) = precipitation_amount(j,i) + &
1789	prr(nzb_s_inner(j,i)+1,j,i) * &
1790	hyrho(nzb_s_inner(j,i)+1) * dt_3d
1791	ENDIF
1792
1793	END SUBROUTINE sedimentation_rain_ij
1794
1795	!------------------------------------------------------------------------------!
1796	! Call for all optimizations
1797	!------------------------------------------------------------------------------!
1798	!
1799	!-- This function computes the gamma function (Press et al., 1992).
1800	!-- The gamma function is needed for the calculation of the evaporation
1801	!-- of rain drops.
1802	FUNCTION gamm( xx )
1803
1804	USE cloud_parameters, &
1805	ONLY: cof, stp
1806
1807	USE kinds
1808
1809	IMPLICIT NONE
1810
1811	INTEGER(iwp) :: j !:
1812
1813	REAL(wp) :: gamm !:
1814	REAL(wp) :: ser !:
1815	REAL(wp) :: tmp !:
1816	REAL(wp) :: x_gamm !:
1817	REAL(wp) :: xx !:
1818	REAL(wp) :: y_gamm !:
1819
1820	x_gamm = xx
1821	y_gamm = x_gamm
1822	tmp = x_gamm + 5.5_wp
1823	tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
1824	ser = 1.000000000190015_wp
1825
1826	DO j = 1, 6
1827	y_gamm = y_gamm + 1.0_wp
1828	ser = ser + cof( j ) / y_gamm
1829	ENDDO
1830
1831	!
1832	!-- Until this point the algorithm computes the logarithm of the gamma
1833	!-- function. Hence, the exponential function is used.
1834	! gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
1835	gamm = EXP( tmp ) * stp * ser / x_gamm
1836
1837	RETURN
1838
1839	END FUNCTION gamm
1840
1841	END MODULE microphysics_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |