Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |