Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |