Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |