Home

Context Navigation

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

Last change on this file since 1683 was 1683, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 75.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |