Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |