Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |