Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |