Home

Context Navigation

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

Last change on this file since 2214 was 2156, checked in by hoffmann, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 86.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |