Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |