Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |