Home

Context Navigation

source: palm/trunk/SOURCE/microphysics_mod.f90 @ 2124

Last change on this file since 2124 was 2101, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 86.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |