Home

Context Navigation

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

Last change on this file since 2600 was 2522, checked in by schwenkel, 7 years ago
minor bugfix
Property svn:keywords set to `Id`
File size: 124.6 KB

Line
1	!> @file microphysics_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
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 2522 2017-10-05 14:20:37Z raasch $
27	! Minor bugfix
28	!
29	! 2375 2017-08-29 14:10:28Z schwenkel
30	! Improved aerosol initilization and some minor bugfixes
31	! for droplet sedimenation
32	!
33	! 2318 2017-07-20 17:27:44Z suehring
34	! Get topography top index via Function call
35	!
36	! 2317 2017-07-20 17:27:19Z suehring
37	! s1 changed to log_sigma
38	!
39	! 2292 2017-06-20 09:51:42Z schwenkel
40	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
41	! includes two more prognostic equations for cloud drop concentration (nc)
42	! and cloud water content (qc).
43	! - The process of activation is parameterized with a simple Twomey
44	! activion scheme or with considering solution and curvature
45	! effects (Khvorostyanov and Curry ,2006).
46	! - The saturation adjustment scheme is replaced by the parameterization
47	! of condensation rates (Khairoutdinov and Kogan, 2000, Mon. Wea. Rev.,128).
48	! - All other microphysical processes of Seifert and Beheng are used.
49	! Additionally, in those processes the reduction of cloud number concentration
50	! is considered.
51	!
52	! 2233 2017-05-30 18:08:54Z suehring
53	!
54	! 2232 2017-05-30 17:47:52Z suehring
55	! Adjustments to new topography and surface concept
56	!
57	! 2155 2017-02-21 09:57:40Z hoffmann
58	! Bugfix in the calculation of microphysical quantities on ghost points.
59	!
60	! 2031 2016-10-21 15:11:58Z knoop
61	! renamed variable rho to rho_ocean
62	!
63	! 2000 2016-08-20 18:09:15Z knoop
64	! Forced header and separation lines into 80 columns
65	!
66	! 1850 2016-04-08 13:29:27Z maronga
67	! Module renamed
68	! Adapted for modularization of microphysics.
69	!
70	! 1845 2016-04-08 08:29:13Z raasch
71	! nzb_2d replaced by nzb_s_inner, Kessler precipitation is stored at surface
72	! point (instead of one point above surface)
73	!
74	! 1831 2016-04-07 13:15:51Z hoffmann
75	! turbulence renamed collision_turbulence,
76	! drizzle renamed cloud_water_sedimentation. cloud_water_sedimentation also
77	! avaialble for microphysics_kessler.
78	!
79	! 1822 2016-04-07 07:49:42Z hoffmann
80	! Unused variables removed.
81	! Kessler scheme integrated.
82	!
83	! 1691 2015-10-26 16:17:44Z maronga
84	! Added new routine calc_precipitation_amount. The routine now allows to account
85	! for precipitation due to sedimenation of cloud (fog) droplets
86	!
87	! 1682 2015-10-07 23:56:08Z knoop
88	! Code annotations made doxygen readable
89	!
90	! 1646 2015-09-02 16:00:10Z hoffmann
91	! Bugfix: Wrong computation of d_mean.
92	!
93	! 1361 2014-04-16 15:17:48Z hoffmann
94	! Bugfix in sedimentation_rain: Index corrected.
95	! Vectorized version of adjust_cloud added.
96	! Little reformatting of the code.
97	!
98	! 1353 2014-04-08 15:21:23Z heinze
99	! REAL constants provided with KIND-attribute
100	!
101	! 1346 2014-03-27 13:18:20Z heinze
102	! Bugfix: REAL constants provided with KIND-attribute especially in call of
103	! intrinsic function like MAX, MIN, SIGN
104	!
105	! 1334 2014-03-25 12:21:40Z heinze
106	! Bugfix: REAL constants provided with KIND-attribute
107	!
108	! 1322 2014-03-20 16:38:49Z raasch
109	! REAL constants defined as wp-kind
110	!
111	! 1320 2014-03-20 08:40:49Z raasch
112	! ONLY-attribute added to USE-statements,
113	! kind-parameters added to all INTEGER and REAL declaration statements,
114	! kinds are defined in new module kinds,
115	! comment fields (!:) to be used for variable explanations added to
116	! all variable declaration statements
117	!
118	! 1241 2013-10-30 11:36:58Z heinze
119	! hyp and rho_ocean have to be calculated at each time step if data from external
120	! file LSF_DATA are used
121	!
122	! 1115 2013-03-26 18:16:16Z hoffmann
123	! microphyical tendencies are calculated in microphysics_control in an optimized
124	! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
125	!
126	! 1106 2013-03-04 05:31:38Z raasch
127	! small changes in code formatting
128	!
129	! 1092 2013-02-02 11:24:22Z raasch
130	! unused variables removed
131	! file put under GPL
132	!
133	! 1065 2012-11-22 17:42:36Z hoffmann
134	! Sedimentation process implemented according to Stevens and Seifert (2008).
135	! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
136	! and Stevens, 2010).
137	!
138	! 1053 2012-11-13 17:11:03Z hoffmann
139	! initial revision
140	!
141	! Description:
142	! ------------
143	!> Calculate bilk cloud microphysics.
144	!------------------------------------------------------------------------------!
145	MODULE microphysics_mod
146
147	USE kinds
148
149	IMPLICIT NONE
150
151	LOGICAL :: cloud_water_sedimentation = .FALSE. !< cloud water sedimentation
152	LOGICAL :: curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory
153	LOGICAL :: limiter_sedimentation = .TRUE. !< sedimentation limiter
154	LOGICAL :: collision_turbulence = .FALSE. !< turbulence effects
155	LOGICAL :: ventilation_effect = .TRUE. !< ventilation effect
156
157	REAL(wp) :: a_1 = 8.69E-4_wp !< coef. in turb. parametrization (cm-2 s3)
158	REAL(wp) :: a_2 = -7.38E-5_wp !< coef. in turb. parametrization (cm-2 s3)
159	REAL(wp) :: a_3 = -1.40E-2_wp !< coef. in turb. parametrization
160	REAL(wp) :: a_term = 9.65_wp !< coef. for terminal velocity (m s-1)
161	REAL(wp) :: a_vent = 0.78_wp !< coef. for ventilation effect
162	REAL(wp) :: b_1 = 11.45E-6_wp !< coef. in turb. parametrization (m)
163	REAL(wp) :: b_2 = 9.68E-6_wp !< coef. in turb. parametrization (m)
164	REAL(wp) :: b_3 = 0.62_wp !< coef. in turb. parametrization
165	REAL(wp) :: b_term = 9.8_wp !< coef. for terminal velocity (m s-1)
166	REAL(wp) :: b_vent = 0.308_wp !< coef. for ventilation effect
167	REAL(wp) :: beta_cc = 3.09E-4_wp !< coef. in turb. parametrization (cm-2 s3)
168	REAL(wp) :: c_1 = 4.82E-6_wp !< coef. in turb. parametrization (m)
169	REAL(wp) :: c_2 = 4.8E-6_wp !< coef. in turb. parametrization (m)
170	REAL(wp) :: c_3 = 0.76_wp !< coef. in turb. parametrization
171	REAL(wp) :: c_const = 0.93_wp !< const. in Taylor-microscale Reynolds number
172	REAL(wp) :: c_evap = 0.7_wp !< constant in evaporation
173	REAL(wp) :: c_term = 600.0_wp !< coef. for terminal velocity (m-1)
174	REAL(wp) :: diff_coeff_l = 0.23E-4_wp !< diffusivity of water vapor (m2 s-1)
175	REAL(wp) :: eps_sb = 1.0E-10_wp !< threshold in two-moments scheme
176	REAL(wp) :: eps_mr = 0.0_wp !< threshold for morrison scheme
177	REAL(wp) :: k_cc = 9.44E09_wp !< const. cloud-cloud kernel (m3 kg-2 s-1)
178	REAL(wp) :: k_cr0 = 4.33_wp !< const. cloud-rain kernel (m3 kg-1 s-1)
179	REAL(wp) :: k_rr = 7.12_wp !< const. rain-rain kernel (m3 kg-1 s-1)
180	REAL(wp) :: k_br = 1000.0_wp !< const. in breakup parametrization (m-1)
181	REAL(wp) :: k_st = 1.2E8_wp !< const. in drizzle parametrization (m-1 s-1)
182	REAL(wp) :: kappa_rr = 60.7_wp !< const. in collision kernel (kg-1/3)
183	REAL(wp) :: kin_vis_air = 1.4086E-5_wp !< kin. viscosity of air (m2 s-1)
184	REAL(wp) :: prec_time_const = 0.001_wp !< coef. in Kessler scheme (s-1)
185	REAL(wp) :: ql_crit = 0.0005_wp !< coef. in Kessler scheme (kg kg-1)
186	REAL(wp) :: schmidt_p_1d3=0.8921121_wp !< Schmidt number0.33333, 0.710.33333
187	REAL(wp) :: sigma_gc = 1.3_wp !< geometric standard deviation cloud droplets
188	REAL(wp) :: thermal_conductivity_l = 2.43E-2_wp !< therm. cond. air (J m-1 s-1 K-1)
189	REAL(wp) :: w_precipitation = 9.65_wp !< maximum terminal velocity (m s-1)
190	REAL(wp) :: x0 = 2.6E-10_wp !< separating drop mass (kg)
191	REAL(wp) :: xamin = 5.24E-19_wp !< average aerosol mass (kg) (~ 0.05Âµm)
192	REAL(wp) :: xcmin = 4.18E-15_wp !< minimum cloud drop size (kg) (~ 1Âµm)
193	REAL(wp) :: xcmax = 2.6E-10_wp !< maximum cloud drop size (kg) (~ 40Âµm)
194	REAL(wp) :: xrmin = 2.6E-10_wp !< minimum rain drop size (kg)
195	REAL(wp) :: xrmax = 5.0E-6_wp !< maximum rain drop site (kg)
196
197	REAL(wp) :: c_sedimentation = 2.0_wp !< Courant number of sedimentation process
198	REAL(wp) :: dpirho_l !< 6.0 / ( pi * rho_l )
199	REAL(wp) :: dry_aerosol_radius = 0.05E-6_wp !< dry aerosol radius
200	REAL(wp) :: dt_micro !< microphysics time step
201	REAL(wp) :: sigma_bulk = 2.0_wp !< width of aerosol spectrum
202	REAL(wp) :: na_init = 100.0E6_wp !< Total particle/aerosol concentration (cm-3)
203	REAL(wp) :: nc_const = 70.0E6_wp !< cloud droplet concentration
204	REAL(wp) :: dt_precipitation = 100.0_wp !< timestep precipitation (s)
205	REAL(wp) :: sed_qc_const !< const. for sedimentation of cloud water
206	REAL(wp) :: pirho_l !< pi * rho_l / 6.0;
207
208	REAL(wp), DIMENSION(:), ALLOCATABLE :: nc_1d !<
209	REAL(wp), DIMENSION(:), ALLOCATABLE :: nr_1d !<
210	REAL(wp), DIMENSION(:), ALLOCATABLE :: pt_1d !<
211	REAL(wp), DIMENSION(:), ALLOCATABLE :: qc_1d !<
212	REAL(wp), DIMENSION(:), ALLOCATABLE :: qr_1d !<
213	REAL(wp), DIMENSION(:), ALLOCATABLE :: q_1d !<
214
215	SAVE
216
217	PRIVATE
218	PUBLIC microphysics_control, microphysics_init
219
220	PUBLIC cloud_water_sedimentation, collision_turbulence, &
221	curvature_solution_effects_bulk, c_sedimentation, &
222	dry_aerosol_radius, dt_precipitation, &
223	limiter_sedimentation, na_init, nc_const, sigma_gc, sigma_bulk, &
224	ventilation_effect
225
226
227	INTERFACE microphysics_control
228	MODULE PROCEDURE microphysics_control
229	MODULE PROCEDURE microphysics_control_ij
230	END INTERFACE microphysics_control
231
232	INTERFACE adjust_cloud
233	MODULE PROCEDURE adjust_cloud
234	MODULE PROCEDURE adjust_cloud_ij
235	END INTERFACE adjust_cloud
236
237	INTERFACE activation
238	MODULE PROCEDURE activation
239	MODULE PROCEDURE activation_ij
240	END INTERFACE activation
241
242	INTERFACE condensation
243	MODULE PROCEDURE condensation
244	MODULE PROCEDURE condensation_ij
245	END INTERFACE condensation
246
247	INTERFACE autoconversion
248	MODULE PROCEDURE autoconversion
249	MODULE PROCEDURE autoconversion_ij
250	END INTERFACE autoconversion
251
252	INTERFACE autoconversion_kessler
253	MODULE PROCEDURE autoconversion_kessler
254	MODULE PROCEDURE autoconversion_kessler_ij
255	END INTERFACE autoconversion_kessler
256
257	INTERFACE accretion
258	MODULE PROCEDURE accretion
259	MODULE PROCEDURE accretion_ij
260	END INTERFACE accretion
261
262	INTERFACE selfcollection_breakup
263	MODULE PROCEDURE selfcollection_breakup
264	MODULE PROCEDURE selfcollection_breakup_ij
265	END INTERFACE selfcollection_breakup
266
267	INTERFACE evaporation_rain
268	MODULE PROCEDURE evaporation_rain
269	MODULE PROCEDURE evaporation_rain_ij
270	END INTERFACE evaporation_rain
271
272	INTERFACE sedimentation_cloud
273	MODULE PROCEDURE sedimentation_cloud
274	MODULE PROCEDURE sedimentation_cloud_ij
275	END INTERFACE sedimentation_cloud
276
277	INTERFACE sedimentation_rain
278	MODULE PROCEDURE sedimentation_rain
279	MODULE PROCEDURE sedimentation_rain_ij
280	END INTERFACE sedimentation_rain
281
282	INTERFACE calc_precipitation_amount
283	MODULE PROCEDURE calc_precipitation_amount
284	MODULE PROCEDURE calc_precipitation_amount_ij
285	END INTERFACE calc_precipitation_amount
286
287	CONTAINS
288	!------------------------------------------------------------------------------!
289	! Description:
290	! ------------
291	!> Initialization of bulk microphysics
292	!------------------------------------------------------------------------------!
293	SUBROUTINE microphysics_init
294
295	USE arrays_3d, &
296	ONLY: dzu
297
298	USE constants, &
299	ONLY: pi
300
301	USE cloud_parameters, &
302	ONLY: molecular_weight_of_solute, rho_l, rho_s, vanthoff
303
304	USE control_parameters, &
305	ONLY: aerosol_nacl, aerosol_c3h4o4 , aerosol_nh4no3, &
306	microphysics_morrison, microphysics_seifert
307
308	USE indices, &
309	ONLY: nzb, nzt
310
311	IMPLICIT NONE
312
313	!
314	!-- constant for the sedimentation of cloud water (2-moment cloud physics)
315	sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l ) &
316	)*( 2.0_wp / 3.0_wp ) &
317	EXP( 5.0_wp * LOG( sigma_gc )**2 )
318
319	!
320	!-- Calculate timestep according to precipitation
321	IF ( microphysics_seifert ) THEN
322	dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) / &
323	w_precipitation
324	ENDIF
325
326	!
327	!-- Set constants for certain aerosol type
328	IF ( microphysics_morrison ) THEN
329	IF ( aerosol_nacl ) THEN
330	molecular_weight_of_solute = 0.05844_wp
331	rho_s = 2165.0_wp
332	vanthoff = 2.0_wp
333	ELSEIF ( aerosol_c3h4o4 ) THEN
334	molecular_weight_of_solute = 0.10406_wp
335	rho_s = 1600.0_wp
336	vanthoff = 1.37_wp
337	ELSEIF ( aerosol_nh4no3 ) THEN
338	molecular_weight_of_solute = 0.08004_wp
339	rho_s = 1720.0_wp
340	vanthoff = 2.31_wp
341	ENDIF
342	ENDIF
343
344	!
345	!-- Pre-calculate frequently calculated fractions of pi and rho_l
346	pirho_l = pi * rho_l / 6.0_wp
347	dpirho_l = 1.0_wp / pirho_l
348
349	!
350	!-- Allocate 1D microphysics arrays
351	ALLOCATE ( pt_1d(nzb:nzt+1), q_1d(nzb:nzt+1), &
352	qc_1d(nzb:nzt+1) )
353
354	IF ( microphysics_morrison .OR. microphysics_seifert ) THEN
355	ALLOCATE ( nc_1d(nzb:nzt+1) )
356	ENDIF
357
358	IF ( microphysics_seifert ) THEN
359	ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) )
360	ENDIF
361
362	END SUBROUTINE microphysics_init
363
364
365	!------------------------------------------------------------------------------!
366	! Description:
367	! ------------
368	!> Control of microphysics for all grid points
369	!------------------------------------------------------------------------------!
370	SUBROUTINE microphysics_control
371
372	USE arrays_3d, &
373	ONLY: hyp, pt_init, prr, zu
374
375	USE cloud_parameters, &
376	ONLY: cp, hyrho, pt_d_t, r_d, t_d_pt
377
378	USE control_parameters, &
379	ONLY: call_microphysics_at_all_substeps, dt_3d, g, &
380	intermediate_timestep_count, large_scale_forcing, &
381	lsf_surf, microphysics_kessler, microphysics_morrison, &
382	microphysics_seifert, pt_surface, &
383	rho_surface,surface_pressure
384
385	USE indices, &
386	ONLY: nzb, nzt
387
388	USE kinds
389
390	USE statistics, &
391	ONLY: weight_pres
392
393	IMPLICIT NONE
394
395	INTEGER(iwp) :: k !<
396
397	REAL(wp) :: t_surface !<
398
399	IF ( large_scale_forcing .AND. lsf_surf ) THEN
400	!
401	!-- Calculate:
402	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
403	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
404	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
405	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
406	DO k = nzb, nzt+1
407	hyp(k) = surface_pressure * 100.0_wp * &
408	( ( t_surface - g / cp * zu(k) ) / &
409	t_surface )**(1.0_wp / 0.286_wp)
410	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
411	t_d_pt(k) = 1.0_wp / pt_d_t(k)
412	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
413	ENDDO
414
415	!
416	!-- Compute reference density
417	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
418	ENDIF
419
420	!
421	!-- Compute length of time step
422	IF ( call_microphysics_at_all_substeps ) THEN
423	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
424	ELSE
425	dt_micro = dt_3d
426	ENDIF
427
428	!
429	!-- Reset precipitation rate
430	IF ( intermediate_timestep_count == 1 ) prr = 0.0_wp
431
432	!
433	!-- Compute cloud physics
434	IF ( microphysics_kessler ) THEN
435
436	CALL autoconversion_kessler
437	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud
438
439	ELSEIF ( microphysics_seifert ) THEN
440
441	CALL adjust_cloud
442	IF ( microphysics_morrison ) CALL activation
443	IF ( microphysics_morrison ) CALL condensation
444	CALL autoconversion
445	CALL accretion
446	CALL selfcollection_breakup
447	CALL evaporation_rain
448	CALL sedimentation_rain
449	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud
450
451	ENDIF
452
453	CALL calc_precipitation_amount
454
455	END SUBROUTINE microphysics_control
456
457	!------------------------------------------------------------------------------!
458	! Description:
459	! ------------
460	!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
461	!> evaporation of rain drops due to too small or too big weights
462	!> of rain drops (Stevens and Seifert, 2008).
463	!------------------------------------------------------------------------------!
464	SUBROUTINE adjust_cloud
465
466	USE arrays_3d, &
467	ONLY: qc, nc, qr, nr
468
469	USE cloud_parameters, &
470	ONLY: hyrho
471
472	USE control_parameters, &
473	ONLY: microphysics_morrison
474
475	USE cpulog, &
476	ONLY: cpu_log, log_point_s
477
478	USE indices, &
479	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
480
481	USE kinds
482
483	IMPLICIT NONE
484
485	INTEGER(iwp) :: i !<
486	INTEGER(iwp) :: j !<
487	INTEGER(iwp) :: k !<
488
489	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
490
491	DO i = nxlg, nxrg
492	DO j = nysg, nyng
493	DO k = nzb+1, nzt
494	IF ( qr(k,j,i) <= eps_sb ) THEN
495	qr(k,j,i) = 0.0_wp
496	nr(k,j,i) = 0.0_wp
497	ELSE
498	IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN
499	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * &
500	MERGE( 1.0_wp, 0.0_wp, &
501	BTEST( wall_flags_0(k,j,i), 0 ) )
502	ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN
503	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * &
504	MERGE( 1.0_wp, 0.0_wp, &
505	BTEST( wall_flags_0(k,j,i), 0 ) )
506	ENDIF
507	ENDIF
508	IF ( microphysics_morrison ) THEN
509	IF ( qc(k,j,i) <= eps_sb ) THEN
510	qc(k,j,i) = 0.0_wp
511	nc(k,j,i) = 0.0_wp
512	ELSE
513	IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN
514	nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin * &
515	MERGE( 1.0_wp, 0.0_wp, &
516	BTEST( wall_flags_0(k,j,i), 0 ) )
517	ENDIF
518	ENDIF
519	ENDIF
520	ENDDO
521	ENDDO
522	ENDDO
523
524	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
525
526	END SUBROUTINE adjust_cloud
527
528	!------------------------------------------------------------------------------!
529	! Description:
530	! ------------
531	!> Calculate number of activated condensation nucleii after simple activation
532	!> scheme of Twomey, 1959.
533	!------------------------------------------------------------------------------!
534	SUBROUTINE activation
535
536	USE arrays_3d, &
537	ONLY: hyp, nc, nr, pt, q, qc, qr
538
539	USE cloud_parameters, &
540	ONLY: hyrho, l_d_cp, l_d_r, l_v, molecular_weight_of_solute, &
541	molecular_weight_of_water, rho_l, rho_s, r_v, t_d_pt, &
542	vanthoff
543
544	USE constants, &
545	ONLY: pi
546
547	USE cpulog, &
548	ONLY: cpu_log, log_point_s
549
550	USE indices, &
551	ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt
552
553	USE kinds
554
555	USE control_parameters, &
556	ONLY: simulated_time
557
558	IMPLICIT NONE
559
560	INTEGER(iwp) :: i !<
561	INTEGER(iwp) :: j !<
562	INTEGER(iwp) :: k !<
563
564	REAL(wp) :: activ !<
565	REAL(wp) :: afactor !<
566	REAL(wp) :: alpha !<
567	REAL(wp) :: beta_act !<
568	REAL(wp) :: bfactor !<
569	REAL(wp) :: e_s !<
570	REAL(wp) :: k_act !<
571	REAL(wp) :: n_act !<
572	REAL(wp) :: n_ccn !<
573	REAL(wp) :: q_s !<
574	REAL(wp) :: s_0 !<
575	REAL(wp) :: sat !<
576	REAL(wp) :: sat_max !<
577	REAL(wp) :: sigma !<
578	REAL(wp) :: sigma_act !<
579	REAL(wp) :: t_int !<
580	REAL(wp) :: t_l !<
581
582	CALL cpu_log( log_point_s(65), 'activation', 'start' )
583
584	DO i = nxlg, nxrg
585	DO j = nysg, nyng
586	DO k = nzb+1, nzt
587
588	!
589	!-- Actual liquid water temperature:
590	t_l = t_d_pt(k) * pt(k,j,i)
591
592	!
593	!-- Calculate actual temperature
594	t_int = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
595	!
596	!-- Saturation vapor pressure at t_l:
597	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
598	( t_l - 35.86_wp ) &
599	)
600	!
601	!-- Computation of saturation humidity:
602	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
603	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
604	q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / &
605	( 1.0_wp + alpha * q_s )
606
607	!-- Supersaturation:
608	sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
609
610	!
611	!-- Prescribe parameters for activation
612	!-- (see: Bott + Trautmann, 2002, Atm. Res., 64)
613	k_act = 0.7_wp
614	activ = 0.0_wp
615
616
617	IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
618	!
619	!-- Compute the number of activated Aerosols
620	!-- (see: Twomey, 1959, Pure and applied Geophysics, 43)
621	n_act = na_init * sat**k_act
622	!
623	!-- Compute the number of cloud droplets
624	!-- (see: Morrison + Grabowski, 2007, JAS, 64)
625	! activ = MAX( n_act - nc(k,j,i), 0.0_wp) / dt_micro
626
627	!
628	!-- Compute activation rate after Khairoutdinov and Kogan
629	!-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
630	sat_max = 1.0_wp / 100.0_wp
631	activ = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN &
632	( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) / &
633	dt_micro
634	ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk ) THEN
635	!
636	!-- Curvature effect (afactor) with surface tension
637	!-- parameterization by Straka (2009)
638	sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
639	afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
640	!
641	!-- Solute effect (bfactor)
642	bfactor = vanthoff * molecular_weight_of_water * &
643	rho_s / ( molecular_weight_of_solute * rho_l )
644
645	!
646	!-- Prescribe power index that describes the soluble fraction
647	!-- of an aerosol particle (beta)
648	!-- (see: Morrison + Grabowski, 2007, JAS, 64)
649	beta_act = 0.5_wp
650	sigma_act = sigma_bulk**( 1.0_wp + beta_act )
651	!
652	!-- Calculate mean geometric supersaturation (s_0) with
653	!-- parameterization by Khvorostyanov and Curry (2006)
654	s_0 = dry_aerosol_radius *(- ( 1.0_wp + beta_act ) ) &
655	( 4.0_wp * afactor*3 / ( 27.0_wp bfactor ) )**0.5_wp
656
657	!
658	!-- Calculate number of activated CCN as a function of
659	!-- supersaturation and taking Koehler theory into account
660	!-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
661	n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( &
662	LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
663	activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp )
664	ENDIF
665
666	nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)
667
668	ENDDO
669	ENDDO
670	ENDDO
671
672	CALL cpu_log( log_point_s(65), 'activation', 'stop' )
673
674	END SUBROUTINE activation
675
676
677	!------------------------------------------------------------------------------!
678	! Description:
679	! ------------
680	!> Calculate condensation rate for cloud water content (after Khairoutdinov and
681	!> Kogan, 2000).
682	!------------------------------------------------------------------------------!
683	SUBROUTINE condensation
684
685	USE arrays_3d, &
686	ONLY: hyp, nr, pt, q, qc, qr, nc
687
688	USE cloud_parameters, &
689	ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
690
691	USE constants, &
692	ONLY: pi
693
694	USE cpulog, &
695	ONLY: cpu_log, log_point_s
696
697	USE indices, &
698	ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt
699
700	USE kinds
701
702	USE control_parameters, &
703	ONLY: simulated_time
704
705	IMPLICIT NONE
706
707	INTEGER(iwp) :: i !<
708	INTEGER(iwp) :: j !<
709	INTEGER(iwp) :: k !<
710
711	REAL(wp) :: alpha !<
712	REAL(wp) :: cond !<
713	REAL(wp) :: cond_max !<
714	REAL(wp) :: dc !<
715	REAL(wp) :: e_s !<
716	REAL(wp) :: evap !<
717	REAL(wp) :: evap_nc !<
718	REAL(wp) :: g_fac !<
719	REAL(wp) :: nc_0 !<
720	REAL(wp) :: q_s !<
721	REAL(wp) :: sat !<
722	REAL(wp) :: t_l !<
723	REAL(wp) :: temp !<
724	REAL(wp) :: xc !<
725
726	CALL cpu_log( log_point_s(66), 'condensation', 'start' )
727
728	DO i = nxlg, nxrg
729	DO j = nysg, nyng
730	DO k = nzb+1, nzt
731	!
732	!-- Actual liquid water temperature:
733	t_l = t_d_pt(k) * pt(k,j,i)
734	!
735	!-- Saturation vapor pressure at t_l:
736	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
737	( t_l - 35.86_wp ) &
738	)
739	!
740	!-- Computation of saturation humidity:
741	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
742	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
743	q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / &
744	( 1.0_wp + alpha * q_s )
745
746	!-- Supersaturation:
747	sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
748
749	!
750	!-- Actual temperature:
751	temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
752
753	g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
754	l_v / ( thermal_conductivity_l * temp ) &
755	+ r_v * temp / ( diff_coeff_l * e_s ) &
756	)
757	!
758	!-- Mean weight of cloud drops
759	IF ( nc(k,j,i) <= 0.0_wp) CYCLE
760	xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
761	!
762	!-- Weight averaged diameter of cloud drops:
763	dc = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
764	!
765	!-- Integral diameter of cloud drops
766	nc_0 = nc(k,j,i) * dc
767	!
768	!-- Condensation needs only to be calculated in supersaturated regions
769	IF ( sat > 0.0_wp ) THEN
770	!
771	!-- Condensation rate of cloud water content
772	!-- after KK scheme.
773	!-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
774	cond = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
775	cond_max = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
776	cond = MIN( cond, cond_max / dt_micro )
777
778	qc(k,j,i) = qc(k,j,i) + cond * dt_micro
779	ELSEIF ( sat < 0.0_wp ) THEN
780	evap = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
781	evap = MAX( evap, -qc(k,j,i) / dt_micro )
782
783	qc(k,j,i) = qc(k,j,i) + evap * dt_micro
784	ENDIF
785	IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN
786	nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin
787	ENDIF
788	ENDDO
789	ENDDO
790	ENDDO
791
792	CALL cpu_log( log_point_s(66), 'condensation', 'stop' )
793
794	END SUBROUTINE condensation
795
796
797	!------------------------------------------------------------------------------!
798	! Description:
799	! ------------
800	!> Autoconversion rate (Seifert and Beheng, 2006).
801	!------------------------------------------------------------------------------!
802	SUBROUTINE autoconversion
803
804	USE arrays_3d, &
805	ONLY: diss, dzu, nc, nr, qc, qr
806
807	USE cloud_parameters, &
808	ONLY: hyrho
809
810	USE control_parameters, &
811	ONLY: microphysics_morrison, rho_surface
812
813	USE cpulog, &
814	ONLY: cpu_log, log_point_s
815
816	USE grid_variables, &
817	ONLY: dx, dy
818
819	USE indices, &
820	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
821
822	USE kinds
823
824	IMPLICIT NONE
825
826	INTEGER(iwp) :: i !<
827	INTEGER(iwp) :: j !<
828	INTEGER(iwp) :: k !<
829
830	REAL(wp) :: alpha_cc !<
831	REAL(wp) :: autocon !<
832	REAL(wp) :: dissipation !<
833	REAL(wp) :: flag !< flag to mask topography grid points
834	REAL(wp) :: k_au !<
835	REAL(wp) :: l_mix !<
836	REAL(wp) :: nc_auto !<
837	REAL(wp) :: nu_c !<
838	REAL(wp) :: phi_au !<
839	REAL(wp) :: r_cc !<
840	REAL(wp) :: rc !<
841	REAL(wp) :: re_lambda !<
842	REAL(wp) :: sigma_cc !<
843	REAL(wp) :: tau_cloud !<
844	REAL(wp) :: xc !<
845
846	CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
847
848	DO i = nxlg, nxrg
849	DO j = nysg, nyng
850	DO k = nzb+1, nzt
851	!
852	!-- Predetermine flag to mask topography
853	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
854
855	nc_auto = MERGE( nc(k,j,i), nc_const, microphysics_morrison )
856
857	IF ( qc(k,j,i) > eps_sb .AND. nc_auto > eps_mr ) THEN
858
859	k_au = k_cc / ( 20.0_wp * x0 )
860	!
861	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
862	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
863	tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + &
864	qc(k,j,i) ), 0.0_wp )
865	!
866	!-- Universal function for autoconversion process
867	!-- (Seifert and Beheng, 2006):
868	phi_au = 600.0_wp * tau_cloud*0.68_wp &
869	( 1.0_wp - tau_cloud0.68_wp )3
870	!
871	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
872	!-- (Use constant nu_c = 1.0_wp instead?)
873	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
874	!
875	!-- Mean weight of cloud droplets:
876	xc = MAX( hyrho(k) * qc(k,j,i) / nc_auto, xcmin)
877	!
878	!-- Parameterized turbulence effects on autoconversion (Seifert,
879	!-- Nuijens and Stevens, 2010)
880	IF ( collision_turbulence ) THEN
881	!
882	!-- Weight averaged radius of cloud droplets:
883	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
884
885	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
886	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
887	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
888	!
889	!-- Mixing length (neglecting distance to ground and
890	!-- stratification)
891	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
892	!
893	!-- Limit dissipation rate according to Seifert, Nuijens and
894	!-- Stevens (2010)
895	dissipation = MIN( 0.06_wp, diss(k,j,i) )
896	!
897	!-- Compute Taylor-microscale Reynolds number:
898	re_lambda = 6.0_wp / 11.0_wp * &
899	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
900	SQRT( 15.0_wp / kin_vis_air ) * &
901	dissipation**( 1.0_wp / 6.0_wp )
902	!
903	!-- The factor of 1.0E4 is needed to convert the dissipation
904	!-- rate from m2 s-3 to cm2 s-3.
905	k_au = k_au * ( 1.0_wp + &
906	dissipation * 1.0E4_wp * &
907	( re_lambda * 1.0E-3_wp )*0.25_wp &
908	( alpha_cc * EXP( -1.0_wp * ( ( rc - &
909	r_cc ) / &
910	sigma_cc )**2 &
911	) + beta_cc &
912	) &
913	)
914	ENDIF
915	!
916	!-- Autoconversion rate (Seifert and Beheng, 2006):
917	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
918	( nu_c + 1.0_wp )*2 qc(k,j,i)*2 xc*2 &
919	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
920	rho_surface
921	autocon = MIN( autocon, qc(k,j,i) / dt_micro )
922
923	qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag
924	qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag
925	nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro &
926	* flag
927	IF ( microphysics_morrison ) THEN
928	nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp * &
929	autocon / x0 * hyrho(k) * dt_micro * flag )
930	ENDIF
931
932	ENDIF
933
934	ENDDO
935	ENDDO
936	ENDDO
937
938	CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
939
940	END SUBROUTINE autoconversion
941
942
943	!------------------------------------------------------------------------------!
944	! Description:
945	! ------------
946	!> Autoconversion process (Kessler, 1969).
947	!------------------------------------------------------------------------------!
948	SUBROUTINE autoconversion_kessler
949
950	USE arrays_3d, &
951	ONLY: dzw, pt, prr, q, qc
952
953	USE cloud_parameters, &
954	ONLY: l_d_cp, pt_d_t
955
956	USE indices, &
957	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
958
959	USE kinds
960
961	USE surface_mod, &
962	ONLY: get_topography_top_index
963
964
965	IMPLICIT NONE
966
967	INTEGER(iwp) :: i !<
968	INTEGER(iwp) :: j !<
969	INTEGER(iwp) :: k !<
970	INTEGER(iwp) :: k_wall !< topgraphy top index
971
972	REAL(wp) :: dqdt_precip !<
973	REAL(wp) :: flag !< flag to mask topography grid points
974
975	DO i = nxlg, nxrg
976	DO j = nysg, nyng
977	!
978	!-- Determine vertical index of topography top
979	k_wall = get_topography_top_index( j, i, 's' )
980	DO k = nzb+1, nzt
981	!
982	!-- Predetermine flag to mask topography
983	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
984
985	IF ( qc(k,j,i) > ql_crit ) THEN
986	dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
987	ELSE
988	dqdt_precip = 0.0_wp
989	ENDIF
990
991	qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
992	q(k,j,i) = q(k,j,i) - dqdt_precip * dt_micro * flag
993	pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * l_d_cp * &
994	pt_d_t(k) * flag
995
996	!
997	!-- Compute the rain rate (stored on surface grid point)
998	prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
999
1000	ENDDO
1001	ENDDO
1002	ENDDO
1003
1004	END SUBROUTINE autoconversion_kessler
1005
1006
1007	!------------------------------------------------------------------------------!
1008	! Description:
1009	! ------------
1010	!> Accretion rate (Seifert and Beheng, 2006).
1011	!------------------------------------------------------------------------------!
1012	SUBROUTINE accretion
1013
1014	USE arrays_3d, &
1015	ONLY: diss, qc, qr, nc
1016
1017	USE cloud_parameters, &
1018	ONLY: hyrho
1019
1020	USE control_parameters, &
1021	ONLY: microphysics_morrison, rho_surface
1022
1023	USE cpulog, &
1024	ONLY: cpu_log, log_point_s
1025
1026	USE indices, &
1027	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
1028
1029	USE kinds
1030
1031	IMPLICIT NONE
1032
1033	INTEGER(iwp) :: i !<
1034	INTEGER(iwp) :: j !<
1035	INTEGER(iwp) :: k !<
1036
1037	REAL(wp) :: accr !<
1038	REAL(wp) :: flag !< flag to mask topography grid points
1039	REAL(wp) :: k_cr !<
1040	REAL(wp) :: nc_accr !<
1041	REAL(wp) :: phi_ac !<
1042	REAL(wp) :: tau_cloud !<
1043	REAL(wp) :: xc !<
1044
1045
1046	CALL cpu_log( log_point_s(56), 'accretion', 'start' )
1047
1048	DO i = nxlg, nxrg
1049	DO j = nysg, nyng
1050	DO k = nzb+1, nzt
1051	!
1052	!-- Predetermine flag to mask topography
1053	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1054
1055	nc_accr = MERGE( nc(k,j,i), nc_const, microphysics_morrison )
1056
1057	IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) &
1058	.AND. ( nc_accr > eps_mr ) ) THEN
1059	!
1060	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
1061	tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
1062	!
1063	!-- Universal function for accretion process (Seifert and
1064	!-- Beheng, 2001):
1065	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
1066
1067	!
1068	!-- Mean weight of cloud drops
1069	xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
1070	!
1071	!-- Parameterized turbulence effects on autoconversion (Seifert,
1072	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
1073	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
1074	IF ( collision_turbulence ) THEN
1075	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
1076	MIN( 600.0_wp, &
1077	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
1078	)
1079	ELSE
1080	k_cr = k_cr0
1081	ENDIF
1082	!
1083	!-- Accretion rate (Seifert and Beheng, 2006):
1084	accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * &
1085	SQRT( rho_surface * hyrho(k) )
1086	accr = MIN( accr, qc(k,j,i) / dt_micro )
1087
1088	qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
1089	qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
1090	IF ( microphysics_morrison ) THEN
1091	nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), &
1092	accr / xc * hyrho(k) * dt_micro * flag)
1093	ENDIF
1094
1095	ENDIF
1096
1097	ENDDO
1098	ENDDO
1099	ENDDO
1100
1101	CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
1102
1103	END SUBROUTINE accretion
1104
1105
1106	!------------------------------------------------------------------------------!
1107	! Description:
1108	! ------------
1109	!> Collisional breakup rate (Seifert, 2008).
1110	!------------------------------------------------------------------------------!
1111	SUBROUTINE selfcollection_breakup
1112
1113	USE arrays_3d, &
1114	ONLY: nr, qr
1115
1116	USE cloud_parameters, &
1117	ONLY: hyrho
1118
1119	USE control_parameters, &
1120	ONLY: rho_surface
1121
1122	USE cpulog, &
1123	ONLY: cpu_log, log_point_s
1124
1125	USE indices, &
1126	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
1127
1128	USE kinds
1129
1130	IMPLICIT NONE
1131
1132	INTEGER(iwp) :: i !<
1133	INTEGER(iwp) :: j !<
1134	INTEGER(iwp) :: k !<
1135
1136	REAL(wp) :: breakup !<
1137	REAL(wp) :: dr !<
1138	REAL(wp) :: flag !< flag to mask topography grid points
1139	REAL(wp) :: phi_br !<
1140	REAL(wp) :: selfcoll !<
1141
1142	CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
1143
1144	DO i = nxlg, nxrg
1145	DO j = nysg, nyng
1146	DO k = nzb+1, nzt
1147	!
1148	!-- Predetermine flag to mask topography
1149	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1150
1151	IF ( qr(k,j,i) > eps_sb ) THEN
1152	!
1153	!-- Selfcollection rate (Seifert and Beheng, 2001):
1154	selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * &
1155	SQRT( hyrho(k) * rho_surface )
1156	!
1157	!-- Weight averaged diameter of rain drops:
1158	dr = ( hyrho(k) * qr(k,j,i) / &
1159	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
1160	!
1161	!-- Collisional breakup rate (Seifert, 2008):
1162	IF ( dr >= 0.3E-3_wp ) THEN
1163	phi_br = k_br * ( dr - 1.1E-3_wp )
1164	breakup = selfcoll * ( phi_br + 1.0_wp )
1165	ELSE
1166	breakup = 0.0_wp
1167	ENDIF
1168
1169	selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
1170	nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
1171
1172	ENDIF
1173	ENDDO
1174	ENDDO
1175	ENDDO
1176
1177	CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
1178
1179	END SUBROUTINE selfcollection_breakup
1180
1181
1182	!------------------------------------------------------------------------------!
1183	! Description:
1184	! ------------
1185	!> Evaporation of precipitable water. Condensation is neglected for
1186	!> precipitable water.
1187	!------------------------------------------------------------------------------!
1188	SUBROUTINE evaporation_rain
1189
1190	USE arrays_3d, &
1191	ONLY: hyp, nr, pt, q, qc, qr
1192
1193	USE cloud_parameters, &
1194	ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
1195
1196	USE constants, &
1197	ONLY: pi
1198
1199	USE cpulog, &
1200	ONLY: cpu_log, log_point_s
1201
1202	USE indices, &
1203	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
1204
1205	USE kinds
1206
1207	IMPLICIT NONE
1208
1209	INTEGER(iwp) :: i !<
1210	INTEGER(iwp) :: j !<
1211	INTEGER(iwp) :: k !<
1212
1213	REAL(wp) :: alpha !<
1214	REAL(wp) :: dr !<
1215	REAL(wp) :: e_s !<
1216	REAL(wp) :: evap !<
1217	REAL(wp) :: evap_nr !<
1218	REAL(wp) :: f_vent !<
1219	REAL(wp) :: flag !< flag to mask topography grid points
1220	REAL(wp) :: g_evap !<
1221	REAL(wp) :: lambda_r !<
1222	REAL(wp) :: mu_r !<
1223	REAL(wp) :: mu_r_2 !<
1224	REAL(wp) :: mu_r_5d2 !<
1225	REAL(wp) :: nr_0 !<
1226	REAL(wp) :: q_s !<
1227	REAL(wp) :: sat !<
1228	REAL(wp) :: t_l !<
1229	REAL(wp) :: temp !<
1230	REAL(wp) :: xr !<
1231
1232	CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
1233
1234	DO i = nxlg, nxrg
1235	DO j = nysg, nyng
1236	DO k = nzb+1, nzt
1237	!
1238	!-- Predetermine flag to mask topography
1239	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1240
1241	IF ( qr(k,j,i) > eps_sb ) THEN
1242	!
1243	!-- Actual liquid water temperature:
1244	t_l = t_d_pt(k) * pt(k,j,i)
1245	!
1246	!-- Saturation vapor pressure at t_l:
1247	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
1248	( t_l - 35.86_wp ) &
1249	)
1250	!
1251	!-- Computation of saturation humidity:
1252	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
1253	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
1254	q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / &
1255	( 1.0_wp + alpha * q_s )
1256	!
1257	!-- Supersaturation:
1258	sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
1259	!
1260	!-- Evaporation needs only to be calculated in subsaturated regions
1261	IF ( sat < 0.0_wp ) THEN
1262	!
1263	!-- Actual temperature:
1264	temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
1265
1266	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
1267	l_v / ( thermal_conductivity_l * temp ) &
1268	+ r_v * temp / ( diff_coeff_l * e_s ) &
1269	)
1270	!
1271	!-- Mean weight of rain drops
1272	xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
1273	!
1274	!-- Weight averaged diameter of rain drops:
1275	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
1276	!
1277	!-- Compute ventilation factor and intercept parameter
1278	!-- (Seifert and Beheng, 2006; Seifert, 2008):
1279	IF ( ventilation_effect ) THEN
1280	!
1281	!-- Shape parameter of gamma distribution (Milbrandt and Yau,
1282	!-- 2005; Stevens and Seifert, 2008):
1283	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
1284	( dr - 1.4E-3_wp ) ) )
1285	!
1286	!-- Slope parameter of gamma distribution (Seifert, 2008):
1287	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
1288	( mu_r + 1.0_wp ) &
1289	)**( 1.0_wp / 3.0_wp ) / dr
1290
1291	mu_r_2 = mu_r + 2.0_wp
1292	mu_r_5d2 = mu_r + 2.5_wp
1293
1294	f_vent = a_vent * gamm( mu_r_2 ) * &
1295	lambda_r*( -mu_r_2 ) + b_vent &
1296	schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
1297	gamm( mu_r_5d2 ) * lambda_r*( -mu_r_5d2 ) &
1298	( 1.0_wp - &
1299	0.5_wp * ( b_term / a_term ) * &
1300	( lambda_r / ( c_term + lambda_r ) &
1301	)**mu_r_5d2 - &
1302	0.125_wp * ( b_term / a_term )*2 &
1303	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
1304	)**mu_r_5d2 - &
1305	0.0625_wp * ( b_term / a_term )*3 &
1306	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
1307	)**mu_r_5d2 - &
1308	0.0390625_wp * ( b_term / a_term )*4 &
1309	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
1310	)**mu_r_5d2 &
1311	)
1312
1313	nr_0 = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) / &
1314	gamm( mu_r + 1.0_wp )
1315	ELSE
1316	f_vent = 1.0_wp
1317	nr_0 = nr(k,j,i) * dr
1318	ENDIF
1319	!
1320	!-- Evaporation rate of rain water content (Seifert and
1321	!-- Beheng, 2006):
1322	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / &
1323	hyrho(k)
1324	evap = MAX( evap, -qr(k,j,i) / dt_micro )
1325	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
1326	-nr(k,j,i) / dt_micro )
1327
1328	qr(k,j,i) = qr(k,j,i) + evap * dt_micro * flag
1329	nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag
1330
1331	ENDIF
1332	ENDIF
1333
1334	ENDDO
1335	ENDDO
1336	ENDDO
1337
1338	CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
1339
1340	END SUBROUTINE evaporation_rain
1341
1342
1343	!------------------------------------------------------------------------------!
1344	! Description:
1345	! ------------
1346	!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
1347	!------------------------------------------------------------------------------!
1348	SUBROUTINE sedimentation_cloud
1349
1350	USE arrays_3d, &
1351	ONLY: ddzu, dzu, nc, pt, prr, q, qc
1352
1353	USE cloud_parameters, &
1354	ONLY: hyrho, l_d_cp, pt_d_t
1355
1356	USE control_parameters, &
1357	ONLY: call_microphysics_at_all_substeps, &
1358	intermediate_timestep_count, microphysics_morrison
1359
1360	USE cpulog, &
1361	ONLY: cpu_log, log_point_s
1362
1363	USE indices, &
1364	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
1365
1366	USE kinds
1367
1368	USE statistics, &
1369	ONLY: weight_substep
1370
1371
1372	IMPLICIT NONE
1373
1374	INTEGER(iwp) :: i !<
1375	INTEGER(iwp) :: j !<
1376	INTEGER(iwp) :: k !<
1377
1378	REAL(wp) :: flag !< flag to mask topography grid points
1379	REAL(wp) :: nc_sedi !<
1380
1381	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
1382	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nc !<
1383
1384
1385	CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
1386
1387	sed_qc(nzt+1) = 0.0_wp
1388	sed_nc(nzt+1) = 0.0_wp
1389
1390	DO i = nxlg, nxrg
1391	DO j = nysg, nyng
1392	DO k = nzt, nzb+1, -1
1393	!
1394	!-- Predetermine flag to mask topography
1395	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1396	nc_sedi = MERGE ( nc(k,j,i), nc_const, microphysics_morrison )
1397
1398	!
1399	!-- Sedimentation fluxes for number concentration are only calculated
1400	!-- for cloud_scheme = 'morrison'
1401	IF ( microphysics_morrison ) THEN
1402	IF ( qc(k,j,i) > eps_sb .AND. nc(k,j,i) > eps_mr ) THEN
1403	sed_nc(k) = sed_qc_const * &
1404	( qc(k,j,i) * hyrho(k) )*( 2.0_wp / 3.0_wp ) &
1405	( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
1406	ELSE
1407	sed_nc(k) = 0.0_wp
1408	ENDIF
1409
1410	sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) * &
1411	nc(k,j,i) / dt_micro + sed_nc(k+1) &
1412	) * flag
1413
1414	nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) * &
1415	ddzu(k+1) / hyrho(k) * dt_micro * flag
1416	ENDIF
1417
1418	IF ( qc(k,j,i) > eps_sb .AND. nc_sedi > eps_mr ) THEN
1419	sed_qc(k) = sed_qc_const * nc_sedi*( -2.0_wp / 3.0_wp ) &
1420	( qc(k,j,i) * hyrho(k) )*( 5.0_wp / 3.0_wp ) &
1421	flag
1422	ELSE
1423	sed_qc(k) = 0.0_wp
1424	ENDIF
1425
1426	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / &
1427	dt_micro + sed_qc(k+1) &
1428	) * flag
1429
1430	q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
1431	ddzu(k+1) / hyrho(k) * dt_micro * flag
1432	qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
1433	ddzu(k+1) / hyrho(k) * dt_micro * flag
1434	pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * &
1435	ddzu(k+1) / hyrho(k) * l_d_cp * &
1436	pt_d_t(k) * dt_micro * flag
1437
1438	!
1439	!-- Compute the precipitation rate due to cloud (fog) droplets
1440	IF ( call_microphysics_at_all_substeps ) THEN
1441	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) &
1442	* weight_substep(intermediate_timestep_count) &
1443	* flag
1444	ELSE
1445	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
1446	ENDIF
1447
1448	ENDDO
1449	ENDDO
1450	ENDDO
1451
1452	CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
1453
1454	END SUBROUTINE sedimentation_cloud
1455
1456
1457	!------------------------------------------------------------------------------!
1458	! Description:
1459	! ------------
1460	!> Computation of sedimentation flux. Implementation according to Stevens
1461	!> and Seifert (2008). Code is based on UCLA-LES.
1462	!------------------------------------------------------------------------------!
1463	SUBROUTINE sedimentation_rain
1464
1465	USE arrays_3d, &
1466	ONLY: ddzu, dzu, nr, pt, prr, q, qr
1467
1468	USE cloud_parameters, &
1469	ONLY: hyrho, l_d_cp, pt_d_t
1470
1471	USE control_parameters, &
1472	ONLY: call_microphysics_at_all_substeps, intermediate_timestep_count
1473	USE cpulog, &
1474	ONLY: cpu_log, log_point_s
1475
1476	USE indices, &
1477	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
1478
1479	USE kinds
1480
1481	USE statistics, &
1482	ONLY: weight_substep
1483
1484	USE surface_mod, &
1485	ONLY : bc_h
1486
1487	IMPLICIT NONE
1488
1489	INTEGER(iwp) :: i !< running index x direction
1490	INTEGER(iwp) :: j !< running index y direction
1491	INTEGER(iwp) :: k !< running index z direction
1492	INTEGER(iwp) :: k_run !<
1493	INTEGER(iwp) :: l !< running index of surface type
1494	INTEGER(iwp) :: m !< running index surface elements
1495	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
1496	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
1497
1498	REAL(wp) :: c_run !<
1499	REAL(wp) :: d_max !<
1500	REAL(wp) :: d_mean !<
1501	REAL(wp) :: d_min !<
1502	REAL(wp) :: dr !<
1503	REAL(wp) :: flux !<
1504	REAL(wp) :: flag !< flag to mask topography grid points
1505	REAL(wp) :: lambda_r !<
1506	REAL(wp) :: mu_r !<
1507	REAL(wp) :: z_run !<
1508
1509	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<
1510	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<
1511	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<
1512	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<
1513	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<
1514	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<
1515	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<
1516	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<
1517
1518	CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
1519
1520	!
1521	!-- Compute velocities
1522	DO i = nxlg, nxrg
1523	DO j = nysg, nyng
1524	DO k = nzb+1, nzt
1525	!
1526	!-- Predetermine flag to mask topography
1527	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1528
1529	IF ( qr(k,j,i) > eps_sb ) THEN
1530	!
1531	!-- Weight averaged diameter of rain drops:
1532	dr = ( hyrho(k) * qr(k,j,i) / &
1533	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
1534	!
1535	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1536	!-- Stevens and Seifert, 2008):
1537	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
1538	( dr - 1.4E-3_wp ) ) )
1539	!
1540	!-- Slope parameter of gamma distribution (Seifert, 2008):
1541	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
1542	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
1543
1544	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
1545	a_term - b_term * ( 1.0_wp + &
1546	c_term / &
1547	lambda_r )*( -1.0_wp &
1548	( mu_r + 1.0_wp ) ) &
1549	) &
1550	) * flag
1551
1552	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
1553	a_term - b_term * ( 1.0_wp + &
1554	c_term / &
1555	lambda_r )*( -1.0_wp &
1556	( mu_r + 4.0_wp ) ) &
1557	) &
1558	) * flag
1559	ELSE
1560	w_nr(k) = 0.0_wp
1561	w_qr(k) = 0.0_wp
1562	ENDIF
1563	ENDDO
1564	!
1565	!-- Adjust boundary values using surface data type.
1566	!-- Upward-facing
1567	surf_s = bc_h(0)%start_index(j,i)
1568	surf_e = bc_h(0)%end_index(j,i)
1569	DO m = surf_s, surf_e
1570	k = bc_h(0)%k(m)
1571	w_nr(k-1) = w_nr(k)
1572	w_qr(k-1) = w_qr(k)
1573	ENDDO
1574	!
1575	!-- Downward-facing
1576	surf_s = bc_h(1)%start_index(j,i)
1577	surf_e = bc_h(1)%end_index(j,i)
1578	DO m = surf_s, surf_e
1579	k = bc_h(1)%k(m)
1580	w_nr(k+1) = w_nr(k)
1581	w_qr(k+1) = w_qr(k)
1582	ENDDO
1583	!
1584	!-- Model top boundary value
1585	w_nr(nzt+1) = 0.0_wp
1586	w_qr(nzt+1) = 0.0_wp
1587	!
1588	!-- Compute Courant number
1589	DO k = nzb+1, nzt
1590	!
1591	!-- Predetermine flag to mask topography
1592	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1593
1594	c_nr(k) = 0.25_wp * ( w_nr(k-1) + &
1595	2.0_wp * w_nr(k) + w_nr(k+1) ) * &
1596	dt_micro * ddzu(k) * flag
1597	c_qr(k) = 0.25_wp * ( w_qr(k-1) + &
1598	2.0_wp * w_qr(k) + w_qr(k+1) ) * &
1599	dt_micro * ddzu(k) * flag
1600	ENDDO
1601	!
1602	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
1603	IF ( limiter_sedimentation ) THEN
1604
1605	DO k = nzb+1, nzt
1606	!
1607	!-- Predetermine flag to mask topography
1608	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1609
1610	d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
1611	d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
1612	d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
1613
1614	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
1615	2.0_wp * d_max, &
1616	ABS( d_mean ) ) &
1617	* flag
1618
1619	d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
1620	d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
1621	d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
1622
1623	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
1624	2.0_wp * d_max, &
1625	ABS( d_mean ) )
1626	ENDDO
1627
1628	ELSE
1629
1630	nr_slope = 0.0_wp
1631	qr_slope = 0.0_wp
1632
1633	ENDIF
1634
1635	sed_nr(nzt+1) = 0.0_wp
1636	sed_qr(nzt+1) = 0.0_wp
1637	!
1638	!-- Compute sedimentation flux
1639	DO k = nzt, nzb+1, -1
1640	!
1641	!-- Predetermine flag to mask topography
1642	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1643	!
1644	!-- Sum up all rain drop number densities which contribute to the flux
1645	!-- through k-1/2
1646	flux = 0.0_wp
1647	z_run = 0.0_wp ! height above z(k)
1648	k_run = k
1649	c_run = MIN( 1.0_wp, c_nr(k) )
1650	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
1651	flux = flux + hyrho(k_run) * &
1652	( nr(k_run,j,i) + nr_slope(k_run) * &
1653	( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run) &
1654	* flag
1655	z_run = z_run + dzu(k_run) * flag
1656	k_run = k_run + 1 * flag
1657	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) &
1658	* flag
1659	ENDDO
1660	!
1661	!-- It is not allowed to sediment more rain drop number density than
1662	!-- available
1663	flux = MIN( flux, &
1664	hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * &
1665	dt_micro &
1666	)
1667
1668	sed_nr(k) = flux / dt_micro * flag
1669	nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * &
1670	ddzu(k+1) / hyrho(k) * dt_micro * flag
1671	!
1672	!-- Sum up all rain water content which contributes to the flux
1673	!-- through k-1/2
1674	flux = 0.0_wp
1675	z_run = 0.0_wp ! height above z(k)
1676	k_run = k
1677	c_run = MIN( 1.0_wp, c_qr(k) )
1678
1679	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
1680
1681	flux = flux + hyrho(k_run) * ( qr(k_run,j,i) + &
1682	qr_slope(k_run) * ( 1.0_wp - c_run ) * &
1683	0.5_wp ) * c_run * dzu(k_run) * flag
1684	z_run = z_run + dzu(k_run) * flag
1685	k_run = k_run + 1 * flag
1686	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) &
1687	* flag
1688
1689	ENDDO
1690	!
1691	!-- It is not allowed to sediment more rain water content than
1692	!-- available
1693	flux = MIN( flux, &
1694	hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * &
1695	dt_micro &
1696	)
1697
1698	sed_qr(k) = flux / dt_micro * flag
1699
1700	qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
1701	ddzu(k+1) / hyrho(k) * dt_micro * flag
1702	q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
1703	ddzu(k+1) / hyrho(k) * dt_micro * flag
1704	pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * &
1705	ddzu(k+1) / hyrho(k) * l_d_cp * &
1706	pt_d_t(k) * dt_micro * flag
1707	!
1708	!-- Compute the rain rate
1709	IF ( call_microphysics_at_all_substeps ) THEN
1710	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) &
1711	* weight_substep(intermediate_timestep_count) &
1712	* flag
1713	ELSE
1714	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
1715	ENDIF
1716
1717	ENDDO
1718	ENDDO
1719	ENDDO
1720
1721	CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
1722
1723	END SUBROUTINE sedimentation_rain
1724
1725
1726	!------------------------------------------------------------------------------!
1727	! Description:
1728	! ------------
1729	!> Computation of the precipitation amount due to gravitational settling of
1730	!> rain and cloud (fog) droplets
1731	!------------------------------------------------------------------------------!
1732	SUBROUTINE calc_precipitation_amount
1733
1734	USE arrays_3d, &
1735	ONLY: precipitation_amount, prr
1736
1737	USE cloud_parameters, &
1738	ONLY: hyrho
1739
1740	USE control_parameters, &
1741	ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d, &
1742	intermediate_timestep_count, intermediate_timestep_count_max,&
1743	precipitation_amount_interval, time_do2d_xy
1744
1745	USE indices, &
1746	ONLY: nxl, nxr, nys, nyn, nzb, nzt, wall_flags_0
1747
1748	USE kinds
1749
1750	USE surface_mod, &
1751	ONLY : bc_h
1752
1753	IMPLICIT NONE
1754
1755	INTEGER(iwp) :: i !< running index x direction
1756	INTEGER(iwp) :: j !< running index y direction
1757	INTEGER(iwp) :: k !< running index y direction
1758	INTEGER(iwp) :: m !< running index surface elements
1759	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
1760	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
1761
1762	IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
1763	( .NOT. call_microphysics_at_all_substeps .OR. &
1764	intermediate_timestep_count == intermediate_timestep_count_max ) ) &
1765	THEN
1766	!
1767	!-- Run over all upward-facing surface elements, i.e. non-natural,
1768	!-- natural and urban
1769	DO m = 1, bc_h(0)%ns
1770	i = bc_h(0)%i(m)
1771	j = bc_h(0)%j(m)
1772	k = bc_h(0)%k(m)
1773	precipitation_amount(j,i) = precipitation_amount(j,i) + &
1774	prr(k,j,i) * hyrho(k) * dt_3d
1775	ENDDO
1776
1777	ENDIF
1778
1779	END SUBROUTINE calc_precipitation_amount
1780
1781
1782	!------------------------------------------------------------------------------!
1783	! Description:
1784	! ------------
1785	!> Control of microphysics for grid points i,j
1786	!------------------------------------------------------------------------------!
1787
1788	SUBROUTINE microphysics_control_ij( i, j )
1789
1790	USE arrays_3d, &
1791	ONLY: hyp, nc, nr, pt, pt_init, prr, q, qc, qr, zu
1792
1793	USE cloud_parameters, &
1794	ONLY: cp, hyrho, pt_d_t, r_d, t_d_pt
1795
1796	USE control_parameters, &
1797	ONLY: call_microphysics_at_all_substeps, dt_3d, g, &
1798	intermediate_timestep_count, large_scale_forcing, &
1799	lsf_surf, microphysics_morrison, microphysics_seifert, &
1800	microphysics_kessler, pt_surface, rho_surface, &
1801	surface_pressure
1802
1803	USE indices, &
1804	ONLY: nzb, nzt
1805
1806	USE kinds
1807
1808	USE statistics, &
1809	ONLY: weight_pres
1810
1811	IMPLICIT NONE
1812
1813	INTEGER(iwp) :: i !<
1814	INTEGER(iwp) :: j !<
1815	INTEGER(iwp) :: k !<
1816
1817	REAL(wp) :: t_surface !<
1818
1819	IF ( large_scale_forcing .AND. lsf_surf ) THEN
1820	!
1821	!-- Calculate:
1822	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
1823	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
1824	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
1825	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
1826	DO k = nzb, nzt+1
1827	hyp(k) = surface_pressure * 100.0_wp * &
1828	( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
1829	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
1830	t_d_pt(k) = 1.0_wp / pt_d_t(k)
1831	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
1832	ENDDO
1833	!
1834	!-- Compute reference density
1835	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
1836	ENDIF
1837
1838	!
1839	!-- Compute length of time step
1840	IF ( call_microphysics_at_all_substeps ) THEN
1841	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
1842	ELSE
1843	dt_micro = dt_3d
1844	ENDIF
1845
1846	!
1847	!-- Use 1d arrays
1848	q_1d(:) = q(:,j,i)
1849	pt_1d(:) = pt(:,j,i)
1850	qc_1d(:) = qc(:,j,i)
1851	IF ( microphysics_seifert ) THEN
1852	qr_1d(:) = qr(:,j,i)
1853	nr_1d(:) = nr(:,j,i)
1854	ENDIF
1855	IF ( microphysics_morrison ) THEN
1856	nc_1d(:) = nc(:,j,i)
1857	ENDIF
1858
1859
1860	!
1861	!-- Reset precipitation rate
1862	IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp
1863
1864	!
1865	!-- Compute cloud physics
1866	IF( microphysics_kessler ) THEN
1867
1868	CALL autoconversion_kessler( i,j )
1869	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud( i,j )
1870
1871	ELSEIF ( microphysics_seifert ) THEN
1872
1873	CALL adjust_cloud( i,j )
1874	IF ( microphysics_morrison ) CALL activation( i,j )
1875	IF ( microphysics_morrison ) CALL condensation( i,j )
1876	CALL autoconversion( i,j )
1877	CALL accretion( i,j )
1878	CALL selfcollection_breakup( i,j )
1879	CALL evaporation_rain( i,j )
1880	CALL sedimentation_rain( i,j )
1881	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud( i,j )
1882
1883	ENDIF
1884
1885	CALL calc_precipitation_amount( i,j )
1886
1887	!
1888	!-- Store results on the 3d arrays
1889	q(:,j,i) = q_1d(:)
1890	pt(:,j,i) = pt_1d(:)
1891	IF ( microphysics_morrison ) THEN
1892	qc(:,j,i) = qc_1d(:)
1893	nc(:,j,i) = nc_1d(:)
1894	ENDIF
1895	IF ( microphysics_seifert ) THEN
1896	qr(:,j,i) = qr_1d(:)
1897	nr(:,j,i) = nr_1d(:)
1898	ENDIF
1899
1900	END SUBROUTINE microphysics_control_ij
1901
1902	!------------------------------------------------------------------------------!
1903	! Description:
1904	! ------------
1905	!> Adjust number of raindrops to avoid nonlinear effects in
1906	!> sedimentation and evaporation of rain drops due to too small or
1907	!> too big weights of rain drops (Stevens and Seifert, 2008).
1908	!> The same procedure is applied to cloud droplets if they are determined
1909	!> prognostically. Call for grid point i,j
1910	!------------------------------------------------------------------------------!
1911	SUBROUTINE adjust_cloud_ij( i, j )
1912
1913	USE cloud_parameters, &
1914	ONLY: hyrho
1915
1916	USE control_parameters, &
1917	ONLY: microphysics_morrison
1918
1919	USE indices, &
1920	ONLY: nzb, nzt, wall_flags_0
1921
1922	USE kinds
1923
1924	IMPLICIT NONE
1925
1926	INTEGER(iwp) :: i !<
1927	INTEGER(iwp) :: j !<
1928	INTEGER(iwp) :: k !<
1929
1930	REAL(wp) :: flag !< flag to indicate first grid level above surface
1931
1932	DO k = nzb+1, nzt
1933	!
1934	!-- Predetermine flag to mask topography
1935	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1936
1937	IF ( qr_1d(k) <= eps_sb ) THEN
1938	qr_1d(k) = 0.0_wp
1939	nr_1d(k) = 0.0_wp
1940	ELSE
1941	!
1942	!-- Adjust number of raindrops to avoid nonlinear effects in
1943	!-- sedimentation and evaporation of rain drops due to too small or
1944	!-- too big weights of rain drops (Stevens and Seifert, 2008).
1945	IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) ) THEN
1946	nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin * flag
1947	ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) ) THEN
1948	nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax * flag
1949	ENDIF
1950
1951	ENDIF
1952
1953	IF ( microphysics_morrison ) THEN
1954	IF ( qc_1d(k) <= eps_sb ) THEN
1955	qc_1d(k) = 0.0_wp
1956	nc_1d(k) = 0.0_wp
1957	ELSE
1958	IF ( nc_1d(k) * xcmin > qc_1d(k) * hyrho(k) ) THEN
1959	nc_1d(k) = qc_1d(k) * hyrho(k) / xamin * flag
1960	ENDIF
1961	ENDIF
1962	ENDIF
1963
1964	ENDDO
1965
1966	END SUBROUTINE adjust_cloud_ij
1967
1968	!------------------------------------------------------------------------------!
1969	! Description:
1970	! ------------
1971	!> Calculate number of activated condensation nucleii after simple activation
1972	!> scheme of Twomey, 1959.
1973	!------------------------------------------------------------------------------!
1974	SUBROUTINE activation_ij( i, j )
1975
1976	USE arrays_3d, &
1977	ONLY: hyp, nr, pt, q, qc, qr, nc
1978
1979	USE cloud_parameters, &
1980	ONLY: hyrho, l_d_cp, l_d_r, l_v, molecular_weight_of_solute, &
1981	molecular_weight_of_water, rho_l, rho_s, r_v, t_d_pt, &
1982	vanthoff
1983
1984	USE constants, &
1985	ONLY: pi
1986
1987	USE cpulog, &
1988	ONLY: cpu_log, log_point_s
1989
1990	USE indices, &
1991	ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt
1992
1993	USE kinds
1994
1995	USE control_parameters, &
1996	ONLY: simulated_time
1997
1998	IMPLICIT NONE
1999
2000	INTEGER(iwp) :: i !<
2001	INTEGER(iwp) :: j !<
2002	INTEGER(iwp) :: k !<
2003
2004	REAL(wp) :: activ !<
2005	REAL(wp) :: afactor !<
2006	REAL(wp) :: alpha !<
2007	REAL(wp) :: beta_act !<
2008	REAL(wp) :: bfactor !<
2009	REAL(wp) :: e_s !<
2010	REAL(wp) :: k_act !<
2011	REAL(wp) :: n_act !<
2012	REAL(wp) :: n_ccn !<
2013	REAL(wp) :: q_s !<
2014	REAL(wp) :: s_0 !<
2015	REAL(wp) :: sat !<
2016	REAL(wp) :: sat_max !<
2017	REAL(wp) :: sigma !<
2018	REAL(wp) :: sigma_act !<
2019	REAL(wp) :: t_int !<
2020	REAL(wp) :: t_l !<
2021
2022	DO k = nzb+1, nzt
2023	!
2024	!-- Actual liquid water temperature:
2025	t_l = t_d_pt(k) * pt_1d(k)
2026
2027	!
2028	!-- Calculate actual temperature
2029	t_int = pt_1d(k) * ( hyp(k) / 100000.0_wp )**0.286_wp
2030	!
2031	!-- Saturation vapor pressure at t_l:
2032	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
2033	( t_l - 35.86_wp ) &
2034	)
2035	!
2036	!-- Computation of saturation humidity:
2037	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
2038	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
2039	q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / &
2040	( 1.0_wp + alpha * q_s )
2041
2042	!-- Supersaturation:
2043	sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
2044
2045	!
2046	!-- Prescribe parameters for activation
2047	!-- (see: Bott + Trautmann, 2002, Atm. Res., 64)
2048	k_act = 0.7_wp
2049	activ = 0.0_wp
2050
2051	IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
2052	!
2053	!-- Compute the number of activated Aerosols
2054	!-- (see: Twomey, 1959, Pure and applied Geophysics, 43)
2055	n_act = na_init * sat**k_act
2056	!
2057	!-- Compute the number of cloud droplets
2058	!-- (see: Morrison + Grabowski, 2007, JAS, 64)
2059	! activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro
2060
2061	!
2062	!-- Compute activation rate after Khairoutdinov and Kogan
2063	!-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
2064	sat_max = 0.8_wp / 100.0_wp
2065	activ = MAX( 0.0_wp, ( (na_init + nc_1d(k) ) * MIN &
2066	( 1.0_wp, ( sat / sat_max )**k_act) - nc_1d(k) ) ) / &
2067	dt_micro
2068
2069	nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)
2070	ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk ) THEN
2071	!
2072	!-- Curvature effect (afactor) with surface tension
2073	!-- parameterization by Straka (2009)
2074	sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
2075	afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
2076	!
2077	!-- Solute effect (bfactor)
2078	bfactor = vanthoff * molecular_weight_of_water * &
2079	rho_s / ( molecular_weight_of_solute * rho_l )
2080
2081	!
2082	!-- Prescribe power index that describes the soluble fraction
2083	!-- of an aerosol particle (beta).
2084	!-- (see: Morrison + Grabowski, 2007, JAS, 64)
2085	beta_act = 0.5_wp
2086	sigma_act = sigma_bulk**( 1.0_wp + beta_act )
2087	!
2088	!-- Calculate mean geometric supersaturation (s_0) with
2089	!-- parameterization by Khvorostyanov and Curry (2006)
2090	s_0 = dry_aerosol_radius *(- ( 1.0_wp + beta_act ) ) &
2091	( 4.0_wp * afactor*3 / ( 27.0_wp bfactor ) )**0.5_wp
2092
2093	!
2094	!-- Calculate number of activated CCN as a function of
2095	!-- supersaturation and taking Koehler theory into account
2096	!-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
2097	n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( &
2098	LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
2099	activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
2100
2101	nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)
2102	ENDIF
2103
2104	ENDDO
2105
2106	END SUBROUTINE activation_ij
2107
2108	!------------------------------------------------------------------------------!
2109	! Description:
2110	! ------------
2111	!> Calculate condensation rate for cloud water content (after Khairoutdinov and
2112	!> Kogan, 2000).
2113	!------------------------------------------------------------------------------!
2114	SUBROUTINE condensation_ij( i, j )
2115
2116	USE arrays_3d, &
2117	ONLY: hyp, nr, pt, q, qc, qr, nc
2118
2119	USE cloud_parameters, &
2120	ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
2121
2122	USE constants, &
2123	ONLY: pi
2124
2125	USE cpulog, &
2126	ONLY: cpu_log, log_point_s
2127
2128	USE indices, &
2129	ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt
2130
2131	USE kinds
2132
2133	USE control_parameters, &
2134	ONLY: simulated_time
2135
2136	IMPLICIT NONE
2137
2138	INTEGER(iwp) :: i !<
2139	INTEGER(iwp) :: j !<
2140	INTEGER(iwp) :: k !<
2141
2142	REAL(wp) :: alpha !<
2143	REAL(wp) :: cond !<
2144	REAL(wp) :: cond_max !<
2145	REAL(wp) :: dc !<
2146	REAL(wp) :: e_s !<
2147	REAL(wp) :: evap !<
2148	REAL(wp) :: evap_nc !<
2149	REAL(wp) :: g_fac !<
2150	REAL(wp) :: nc_0 !<
2151	REAL(wp) :: q_s !<
2152	REAL(wp) :: sat !<
2153	REAL(wp) :: t_l !<
2154	REAL(wp) :: temp !<
2155	REAL(wp) :: xc !<
2156
2157
2158	DO k = nzb+1, nzt
2159	!
2160	!-- Actual liquid water temperature:
2161	t_l = t_d_pt(k) * pt_1d(k)
2162	!
2163	!-- Saturation vapor pressure at t_l:
2164	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
2165	( t_l - 35.86_wp ) &
2166	)
2167	!
2168	!-- Computation of saturation humidity:
2169	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
2170	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
2171	q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / &
2172	( 1.0_wp + alpha * q_s )
2173
2174	!-- Supersaturation:
2175	sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
2176
2177
2178	!
2179	!-- Actual temperature:
2180	temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
2181
2182	g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
2183	l_v / ( thermal_conductivity_l * temp ) &
2184	+ r_v * temp / ( diff_coeff_l * e_s ) &
2185	)
2186	!
2187	!-- Mean weight of cloud drops
2188	IF ( nc_1d(k) <= 0.0_wp) CYCLE
2189	xc = MAX( (hyrho(k) * qc_1d(k) / nc_1d(k)), xcmin)
2190	!
2191	!-- Weight averaged diameter of cloud drops:
2192	dc = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
2193	!
2194	!-- Integral diameter of cloud drops
2195	nc_0 = nc_1d(k) * dc
2196	!
2197	!-- Condensation needs only to be calculated in supersaturated regions
2198	IF ( sat > 0.0_wp ) THEN
2199	!
2200	!-- Condensation rate of cloud water content
2201	!-- after KK scheme.
2202	!-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
2203	cond = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
2204	cond_max = q_1d(k) - q_s - qc_1d(k) - qr_1d(k)
2205	cond = MIN( cond, cond_max / dt_micro )
2206
2207	qc_1d(k) = qc_1d(k) + cond * dt_micro
2208	ELSEIF ( sat < 0.0_wp ) THEN
2209	evap = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
2210	evap = MAX( evap, -qc_1d(k) / dt_micro )
2211
2212	qc_1d(k) = qc_1d(k) + evap * dt_micro
2213	ENDIF
2214	ENDDO
2215
2216	END SUBROUTINE condensation_ij
2217
2218
2219	!------------------------------------------------------------------------------!
2220	! Description:
2221	! ------------
2222	!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
2223	!------------------------------------------------------------------------------!
2224	SUBROUTINE autoconversion_ij( i, j )
2225
2226	USE arrays_3d, &
2227	ONLY: diss, dzu
2228
2229	USE cloud_parameters, &
2230	ONLY: hyrho
2231
2232	USE control_parameters, &
2233	ONLY: microphysics_morrison, rho_surface
2234
2235	USE grid_variables, &
2236	ONLY: dx, dy
2237
2238	USE indices, &
2239	ONLY: nzb, nzt, wall_flags_0
2240
2241	USE kinds
2242
2243	IMPLICIT NONE
2244
2245	INTEGER(iwp) :: i !<
2246	INTEGER(iwp) :: j !<
2247	INTEGER(iwp) :: k !<
2248
2249	REAL(wp) :: alpha_cc !<
2250	REAL(wp) :: autocon !<
2251	REAL(wp) :: dissipation !<
2252	REAL(wp) :: flag !< flag to indicate first grid level above surface
2253	REAL(wp) :: k_au !<
2254	REAL(wp) :: l_mix !<
2255	REAL(wp) :: nc_auto !<
2256	REAL(wp) :: nu_c !<
2257	REAL(wp) :: phi_au !<
2258	REAL(wp) :: r_cc !<
2259	REAL(wp) :: rc !<
2260	REAL(wp) :: re_lambda !<
2261	REAL(wp) :: sigma_cc !<
2262	REAL(wp) :: tau_cloud !<
2263	REAL(wp) :: xc !<
2264
2265	DO k = nzb+1, nzt
2266	!
2267	!-- Predetermine flag to mask topography
2268	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2269	nc_auto = MERGE ( nc_1d(k), nc_const, microphysics_morrison )
2270
2271	IF ( qc_1d(k) > eps_sb .AND. nc_auto > eps_mr ) THEN
2272
2273	k_au = k_cc / ( 20.0_wp * x0 )
2274	!
2275	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
2276	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
2277	tau_cloud = MAX( 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) ), &
2278	0.0_wp )
2279	!
2280	!-- Universal function for autoconversion process
2281	!-- (Seifert and Beheng, 2006):
2282	phi_au = 600.0_wp * tau_cloud*0.68_wp ( 1.0_wp - tau_cloud0.68_wp )3
2283	!
2284	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
2285	!-- (Use constant nu_c = 1.0_wp instead?)
2286	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp )
2287	!
2288	!-- Mean weight of cloud droplets:
2289	xc = hyrho(k) * qc_1d(k) / nc_auto
2290	!
2291	!-- Parameterized turbulence effects on autoconversion (Seifert,
2292	!-- Nuijens and Stevens, 2010)
2293	IF ( collision_turbulence ) THEN
2294	!
2295	!-- Weight averaged radius of cloud droplets:
2296	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
2297
2298	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
2299	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
2300	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
2301	!
2302	!-- Mixing length (neglecting distance to ground and stratification)
2303	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
2304	!
2305	!-- Limit dissipation rate according to Seifert, Nuijens and
2306	!-- Stevens (2010)
2307	dissipation = MIN( 0.06_wp, diss(k,j,i) )
2308	!
2309	!-- Compute Taylor-microscale Reynolds number:
2310	re_lambda = 6.0_wp / 11.0_wp * &
2311	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
2312	SQRT( 15.0_wp / kin_vis_air ) * &
2313	dissipation**( 1.0_wp / 6.0_wp )
2314	!
2315	!-- The factor of 1.0E4 is needed to convert the dissipation rate
2316	!-- from m2 s-3 to cm2 s-3.
2317	k_au = k_au * ( 1.0_wp + &
2318	dissipation * 1.0E4_wp * &
2319	( re_lambda * 1.0E-3_wp )*0.25_wp &
2320	( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) / &
2321	sigma_cc )**2 &
2322	) + beta_cc &
2323	) &
2324	)
2325	ENDIF
2326	!
2327	!-- Autoconversion rate (Seifert and Beheng, 2006):
2328	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
2329	( nu_c + 1.0_wp )*2 qc_1d(k)*2 xc*2 &
2330	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
2331	rho_surface
2332	autocon = MIN( autocon, qc_1d(k) / dt_micro )
2333
2334	qr_1d(k) = qr_1d(k) + autocon * dt_micro * flag
2335	qc_1d(k) = qc_1d(k) - autocon * dt_micro * flag
2336	nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro * flag
2337	IF ( microphysics_morrison ) THEN
2338	nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), 2.0_wp * &
2339	autocon / x0 * hyrho(k) * dt_micro * flag )
2340	ENDIF
2341
2342	ENDIF
2343
2344	ENDDO
2345
2346	END SUBROUTINE autoconversion_ij
2347
2348	!------------------------------------------------------------------------------!
2349	! Description:
2350	! ------------
2351	!> Autoconversion process (Kessler, 1969).
2352	!------------------------------------------------------------------------------!
2353	SUBROUTINE autoconversion_kessler_ij( i, j )
2354
2355	USE arrays_3d, &
2356	ONLY: dzw, prr
2357
2358	USE cloud_parameters, &
2359	ONLY: l_d_cp, pt_d_t
2360
2361	USE indices, &
2362	ONLY: nzb, nzt, wall_flags_0
2363
2364	USE kinds
2365
2366	USE surface_mod, &
2367	ONLY: get_topography_top_index
2368
2369
2370	IMPLICIT NONE
2371
2372	INTEGER(iwp) :: i !<
2373	INTEGER(iwp) :: j !<
2374	INTEGER(iwp) :: k !<
2375	INTEGER(iwp) :: k_wall !< topography top index
2376
2377	REAL(wp) :: dqdt_precip !<
2378	REAL(wp) :: flag !< flag to indicate first grid level above surface
2379
2380	!
2381	!-- Determine vertical index of topography top
2382	k_wall = get_topography_top_index( j, i, 's' )
2383	DO k = nzb+1, nzt
2384	!
2385	!-- Predetermine flag to mask topography
2386	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2387
2388	IF ( qc_1d(k) > ql_crit ) THEN
2389	dqdt_precip = prec_time_const * ( qc_1d(k) - ql_crit )
2390	ELSE
2391	dqdt_precip = 0.0_wp
2392	ENDIF
2393
2394	qc_1d(k) = qc_1d(k) - dqdt_precip * dt_micro * flag
2395	q_1d(k) = q_1d(k) - dqdt_precip * dt_micro * flag
2396	pt_1d(k) = pt_1d(k) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k) * flag
2397
2398	!
2399	!-- Compute the rain rate (stored on surface grid point)
2400	prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
2401
2402	ENDDO
2403
2404	END SUBROUTINE autoconversion_kessler_ij
2405
2406	!------------------------------------------------------------------------------!
2407	! Description:
2408	! ------------
2409	!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
2410	!------------------------------------------------------------------------------!
2411	SUBROUTINE accretion_ij( i, j )
2412
2413	USE arrays_3d, &
2414	ONLY: diss
2415
2416	USE cloud_parameters, &
2417	ONLY: hyrho
2418
2419	USE control_parameters, &
2420	ONLY: microphysics_morrison, rho_surface
2421
2422	USE indices, &
2423	ONLY: nzb, nzt, wall_flags_0
2424
2425	USE kinds
2426
2427	IMPLICIT NONE
2428
2429	INTEGER(iwp) :: i !<
2430	INTEGER(iwp) :: j !<
2431	INTEGER(iwp) :: k !<
2432
2433	REAL(wp) :: accr !<
2434	REAL(wp) :: flag !< flag to indicate first grid level above surface
2435	REAL(wp) :: k_cr !<
2436	REAL(wp) :: nc_accr !<
2437	REAL(wp) :: phi_ac !<
2438	REAL(wp) :: tau_cloud !<
2439	REAL(wp) :: xc !<
2440
2441
2442	DO k = nzb+1, nzt
2443	!
2444	!-- Predetermine flag to mask topography
2445	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2446	nc_accr = MERGE ( nc_1d(k), nc_const, microphysics_morrison )
2447
2448	IF ( ( qc_1d(k) > eps_sb ) .AND. ( qr_1d(k) > eps_sb ) .AND. &
2449	( nc_accr > eps_mr ) ) THEN
2450	!
2451	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
2452	tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
2453	!
2454	!-- Universal function for accretion process
2455	!-- (Seifert and Beheng, 2001):
2456	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
2457
2458	!
2459	!-- Mean weight of cloud drops
2460	xc = MAX( (hyrho(k) * qc_1d(k) / nc_accr), xcmin)
2461	!
2462	!-- Parameterized turbulence effects on autoconversion (Seifert,
2463	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
2464	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
2465	IF ( collision_turbulence ) THEN
2466	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
2467	MIN( 600.0_wp, &
2468	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
2469	)
2470	ELSE
2471	k_cr = k_cr0
2472	ENDIF
2473	!
2474	!-- Accretion rate (Seifert and Beheng, 2006):
2475	accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * &
2476	SQRT( rho_surface * hyrho(k) )
2477	accr = MIN( accr, qc_1d(k) / dt_micro )
2478
2479	qr_1d(k) = qr_1d(k) + accr * dt_micro * flag
2480	qc_1d(k) = qc_1d(k) - accr * dt_micro * flag
2481	IF ( microphysics_morrison ) THEN
2482	nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), accr / xc * &
2483	hyrho(k) * dt_micro * flag &
2484	)
2485	ENDIF
2486
2487
2488	ENDIF
2489
2490	ENDDO
2491
2492	END SUBROUTINE accretion_ij
2493
2494
2495	!------------------------------------------------------------------------------!
2496	! Description:
2497	! ------------
2498	!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
2499	!------------------------------------------------------------------------------!
2500	SUBROUTINE selfcollection_breakup_ij( i, j )
2501
2502	USE cloud_parameters, &
2503	ONLY: hyrho
2504
2505	USE control_parameters, &
2506	ONLY: rho_surface
2507
2508	USE indices, &
2509	ONLY: nzb, nzt, wall_flags_0
2510
2511	USE kinds
2512
2513	IMPLICIT NONE
2514
2515	INTEGER(iwp) :: i !<
2516	INTEGER(iwp) :: j !<
2517	INTEGER(iwp) :: k !<
2518
2519	REAL(wp) :: breakup !<
2520	REAL(wp) :: dr !<
2521	REAL(wp) :: flag !< flag to indicate first grid level above surface
2522	REAL(wp) :: phi_br !<
2523	REAL(wp) :: selfcoll !<
2524
2525	DO k = nzb+1, nzt
2526	!
2527	!-- Predetermine flag to mask topography
2528	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2529
2530	IF ( qr_1d(k) > eps_sb ) THEN
2531	!
2532	!-- Selfcollection rate (Seifert and Beheng, 2001):
2533	selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
2534	!
2535	!-- Weight averaged diameter of rain drops:
2536	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
2537	!
2538	!-- Collisional breakup rate (Seifert, 2008):
2539	IF ( dr >= 0.3E-3_wp ) THEN
2540	phi_br = k_br * ( dr - 1.1E-3_wp )
2541	breakup = selfcoll * ( phi_br + 1.0_wp )
2542	ELSE
2543	breakup = 0.0_wp
2544	ENDIF
2545
2546	selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
2547	nr_1d(k) = nr_1d(k) + selfcoll * dt_micro * flag
2548
2549	ENDIF
2550	ENDDO
2551
2552	END SUBROUTINE selfcollection_breakup_ij
2553
2554
2555	!------------------------------------------------------------------------------!
2556	! Description:
2557	! ------------
2558	!> Evaporation of precipitable water. Condensation is neglected for
2559	!> precipitable water. Call for grid point i,j
2560	!------------------------------------------------------------------------------!
2561	SUBROUTINE evaporation_rain_ij( i, j )
2562
2563	USE arrays_3d, &
2564	ONLY: hyp
2565
2566	USE cloud_parameters, &
2567	ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
2568
2569	USE constants, &
2570	ONLY: pi
2571
2572	USE indices, &
2573	ONLY: nzb, nzt, wall_flags_0
2574
2575	USE kinds
2576
2577	IMPLICIT NONE
2578
2579	INTEGER(iwp) :: i !<
2580	INTEGER(iwp) :: j !<
2581	INTEGER(iwp) :: k !<
2582
2583	REAL(wp) :: alpha !<
2584	REAL(wp) :: dr !<
2585	REAL(wp) :: e_s !<
2586	REAL(wp) :: evap !<
2587	REAL(wp) :: evap_nr !<
2588	REAL(wp) :: f_vent !<
2589	REAL(wp) :: flag !< flag to indicate first grid level above surface
2590	REAL(wp) :: g_evap !<
2591	REAL(wp) :: lambda_r !<
2592	REAL(wp) :: mu_r !<
2593	REAL(wp) :: mu_r_2 !<
2594	REAL(wp) :: mu_r_5d2 !<
2595	REAL(wp) :: nr_0 !<
2596	REAL(wp) :: q_s !<
2597	REAL(wp) :: sat !<
2598	REAL(wp) :: t_l !<
2599	REAL(wp) :: temp !<
2600	REAL(wp) :: xr !<
2601
2602	DO k = nzb+1, nzt
2603	!
2604	!-- Predetermine flag to mask topography
2605	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2606
2607	IF ( qr_1d(k) > eps_sb ) THEN
2608	!
2609	!-- Actual liquid water temperature:
2610	t_l = t_d_pt(k) * pt_1d(k)
2611	!
2612	!-- Saturation vapor pressure at t_l:
2613	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
2614	( t_l - 35.86_wp ) &
2615	)
2616	!
2617	!-- Computation of saturation humidity:
2618	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
2619	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
2620	q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
2621	!
2622	!-- Supersaturation:
2623	sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
2624	!
2625	!-- Evaporation needs only to be calculated in subsaturated regions
2626	IF ( sat < 0.0_wp ) THEN
2627	!
2628	!-- Actual temperature:
2629	temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
2630
2631	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v / &
2632	( thermal_conductivity_l * temp ) + &
2633	r_v * temp / ( diff_coeff_l * e_s ) &
2634	)
2635	!
2636	!-- Mean weight of rain drops
2637	xr = hyrho(k) * qr_1d(k) / nr_1d(k)
2638	!
2639	!-- Weight averaged diameter of rain drops:
2640	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
2641	!
2642	!-- Compute ventilation factor and intercept parameter
2643	!-- (Seifert and Beheng, 2006; Seifert, 2008):
2644	IF ( ventilation_effect ) THEN
2645	!
2646	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
2647	!-- Stevens and Seifert, 2008):
2648	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
2649	!
2650	!-- Slope parameter of gamma distribution (Seifert, 2008):
2651	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
2652	( mu_r + 1.0_wp ) &
2653	)**( 1.0_wp / 3.0_wp ) / dr
2654
2655	mu_r_2 = mu_r + 2.0_wp
2656	mu_r_5d2 = mu_r + 2.5_wp
2657
2658	f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) + &
2659	b_vent * schmidt_p_1d3 * &
2660	SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * &
2661	lambda_r*( -mu_r_5d2 ) &
2662	( 1.0_wp - &
2663	0.5_wp * ( b_term / a_term ) * &
2664	( lambda_r / ( c_term + lambda_r ) &
2665	)**mu_r_5d2 - &
2666	0.125_wp * ( b_term / a_term )*2 &
2667	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
2668	)**mu_r_5d2 - &
2669	0.0625_wp * ( b_term / a_term )*3 &
2670	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
2671	)**mu_r_5d2 - &
2672	0.0390625_wp * ( b_term / a_term )*4 &
2673	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
2674	)**mu_r_5d2 &
2675	)
2676
2677	nr_0 = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) / &
2678	gamm( mu_r + 1.0_wp )
2679	ELSE
2680	f_vent = 1.0_wp
2681	nr_0 = nr_1d(k) * dr
2682	ENDIF
2683	!
2684	!-- Evaporation rate of rain water content (Seifert and Beheng, 2006):
2685	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
2686	evap = MAX( evap, -qr_1d(k) / dt_micro )
2687	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
2688	-nr_1d(k) / dt_micro )
2689
2690	qr_1d(k) = qr_1d(k) + evap * dt_micro * flag
2691	nr_1d(k) = nr_1d(k) + evap_nr * dt_micro * flag
2692
2693	ENDIF
2694	ENDIF
2695
2696	ENDDO
2697
2698	END SUBROUTINE evaporation_rain_ij
2699
2700
2701	!------------------------------------------------------------------------------!
2702	! Description:
2703	! ------------
2704	!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
2705	!> Call for grid point i,j
2706	!------------------------------------------------------------------------------!
2707	SUBROUTINE sedimentation_cloud_ij( i, j )
2708
2709	USE arrays_3d, &
2710	ONLY: ddzu, dzu, prr
2711
2712	USE cloud_parameters, &
2713	ONLY: hyrho, l_d_cp, pt_d_t
2714
2715	USE control_parameters, &
2716	ONLY: call_microphysics_at_all_substeps, &
2717	intermediate_timestep_count, microphysics_morrison
2718
2719	USE indices, &
2720	ONLY: nzb, nzb, nzt, wall_flags_0
2721
2722	USE kinds
2723
2724	USE statistics, &
2725	ONLY: weight_substep
2726
2727	IMPLICIT NONE
2728
2729	INTEGER(iwp) :: i !<
2730	INTEGER(iwp) :: j !<
2731	INTEGER(iwp) :: k !<
2732
2733	REAL(wp) :: flag !< flag to indicate first grid level above surface
2734	REAL(wp) :: nc_sedi !<
2735
2736	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nc !<
2737	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
2738
2739	sed_qc(nzt+1) = 0.0_wp
2740	sed_nc(nzt+1) = 0.0_wp
2741
2742
2743	DO k = nzt, nzb+1, -1
2744	!
2745	!-- Predetermine flag to mask topography
2746	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2747	nc_sedi = MERGE( nc_1d(k), nc_const, microphysics_morrison )
2748	!
2749	!-- Sedimentation fluxes for number concentration are only calculated
2750	!-- for cloud_scheme = 'morrison'
2751	IF ( microphysics_morrison ) THEN
2752	IF ( qc_1d(k) > eps_sb .AND. nc_1d(k) > eps_mr ) THEN
2753	sed_nc(k) = sed_qc_const * &
2754	( qc_1d(k) * hyrho(k) )*( 2.0_wp / 3.0_wp ) &
2755	( nc_1d(k) )**( 1.0_wp / 3.0_wp )
2756	ELSE
2757	sed_nc(k) = 0.0_wp
2758	ENDIF
2759
2760	sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) * &
2761	nc_1d(k) / dt_micro + sed_nc(k+1) &
2762	) * flag
2763
2764	nc_1d(k) = nc_1d(k) + ( sed_nc(k+1) - sed_nc(k) ) * &
2765	ddzu(k+1) / hyrho(k) * dt_micro * flag
2766	ENDIF
2767
2768	IF ( qc_1d(k) > eps_sb .AND. nc_sedi > eps_mr ) THEN
2769	sed_qc(k) = sed_qc_const * nc_sedi*( -2.0_wp / 3.0_wp ) &
2770	( qc_1d(k) * hyrho(k) )*( 5.0_wp / 3.0_wp ) flag
2771	ELSE
2772	sed_qc(k) = 0.0_wp
2773	ENDIF
2774
2775	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) / &
2776	dt_micro + sed_qc(k+1) &
2777	) * flag
2778
2779	q_1d(k) = q_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
2780	hyrho(k) * dt_micro * flag
2781	qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
2782	hyrho(k) * dt_micro * flag
2783	pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
2784	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro * flag
2785
2786	!
2787	!-- Compute the precipitation rate of cloud (fog) droplets
2788	IF ( call_microphysics_at_all_substeps ) THEN
2789	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * &
2790	weight_substep(intermediate_timestep_count) * flag
2791	ELSE
2792	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
2793	ENDIF
2794
2795	ENDDO
2796
2797	END SUBROUTINE sedimentation_cloud_ij
2798
2799
2800	!------------------------------------------------------------------------------!
2801	! Description:
2802	! ------------
2803	!> Computation of sedimentation flux. Implementation according to Stevens
2804	!> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j
2805	!------------------------------------------------------------------------------!
2806	SUBROUTINE sedimentation_rain_ij( i, j )
2807
2808	USE arrays_3d, &
2809	ONLY: ddzu, dzu, prr
2810
2811	USE cloud_parameters, &
2812	ONLY: hyrho, l_d_cp, pt_d_t
2813
2814	USE control_parameters, &
2815	ONLY: call_microphysics_at_all_substeps, intermediate_timestep_count
2816
2817	USE indices, &
2818	ONLY: nzb, nzb, nzt, wall_flags_0
2819
2820	USE kinds
2821
2822	USE statistics, &
2823	ONLY: weight_substep
2824
2825	USE surface_mod, &
2826	ONLY : bc_h
2827
2828	IMPLICIT NONE
2829
2830	INTEGER(iwp) :: i !< running index x direction
2831	INTEGER(iwp) :: j !< running index y direction
2832	INTEGER(iwp) :: k !< running index z direction
2833	INTEGER(iwp) :: k_run !<
2834	INTEGER(iwp) :: m !< running index surface elements
2835	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
2836	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
2837
2838	REAL(wp) :: c_run !<
2839	REAL(wp) :: d_max !<
2840	REAL(wp) :: d_mean !<
2841	REAL(wp) :: d_min !<
2842	REAL(wp) :: dr !<
2843	REAL(wp) :: flux !<
2844	REAL(wp) :: flag !< flag to indicate first grid level above surface
2845	REAL(wp) :: lambda_r !<
2846	REAL(wp) :: mu_r !<
2847	REAL(wp) :: z_run !<
2848
2849	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<
2850	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<
2851	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<
2852	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<
2853	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<
2854	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<
2855	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<
2856	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<
2857
2858	!
2859	!-- Compute velocities
2860	DO k = nzb+1, nzt
2861	!
2862	!-- Predetermine flag to mask topography
2863	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2864
2865	IF ( qr_1d(k) > eps_sb ) THEN
2866	!
2867	!-- Weight averaged diameter of rain drops:
2868	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
2869	!
2870	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
2871	!-- Stevens and Seifert, 2008):
2872	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
2873	!
2874	!-- Slope parameter of gamma distribution (Seifert, 2008):
2875	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
2876	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
2877
2878	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
2879	a_term - b_term * ( 1.0_wp + &
2880	c_term / lambda_r )*( -1.0_wp &
2881	( mu_r + 1.0_wp ) ) &
2882	) &
2883	) * flag
2884	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
2885	a_term - b_term * ( 1.0_wp + &
2886	c_term / lambda_r )*( -1.0_wp &
2887	( mu_r + 4.0_wp ) ) &
2888	) &
2889	) * flag
2890	ELSE
2891	w_nr(k) = 0.0_wp
2892	w_qr(k) = 0.0_wp
2893	ENDIF
2894	ENDDO
2895	!
2896	!-- Adjust boundary values using surface data type.
2897	!-- Upward facing non-natural
2898	surf_s = bc_h(0)%start_index(j,i)
2899	surf_e = bc_h(0)%end_index(j,i)
2900	DO m = surf_s, surf_e
2901	k = bc_h(0)%k(m)
2902	w_nr(k-1) = w_nr(k)
2903	w_qr(k-1) = w_qr(k)
2904	ENDDO
2905	!
2906	!-- Downward facing non-natural
2907	surf_s = bc_h(1)%start_index(j,i)
2908	surf_e = bc_h(1)%end_index(j,i)
2909	DO m = surf_s, surf_e
2910	k = bc_h(1)%k(m)
2911	w_nr(k+1) = w_nr(k)
2912	w_qr(k+1) = w_qr(k)
2913	ENDDO
2914	!
2915	!-- Neumann boundary condition at model top
2916	w_nr(nzt+1) = 0.0_wp
2917	w_qr(nzt+1) = 0.0_wp
2918	!
2919	!-- Compute Courant number
2920	DO k = nzb+1, nzt
2921	!
2922	!-- Predetermine flag to mask topography
2923	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2924
2925	c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) * &
2926	dt_micro * ddzu(k) * flag
2927	c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * &
2928	dt_micro * ddzu(k) * flag
2929	ENDDO
2930	!
2931	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
2932	IF ( limiter_sedimentation ) THEN
2933
2934	DO k = nzb+1, nzt
2935	!
2936	!-- Predetermine flag to mask topography
2937	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2938
2939	d_mean = 0.5_wp * ( qr_1d(k+1) - qr_1d(k-1) )
2940	d_min = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
2941	d_max = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
2942
2943	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
2944	2.0_wp * d_max, &
2945	ABS( d_mean ) ) * flag
2946
2947	d_mean = 0.5_wp * ( nr_1d(k+1) - nr_1d(k-1) )
2948	d_min = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
2949	d_max = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
2950
2951	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
2952	2.0_wp * d_max, &
2953	ABS( d_mean ) ) * flag
2954	ENDDO
2955
2956	ELSE
2957
2958	nr_slope = 0.0_wp
2959	qr_slope = 0.0_wp
2960
2961	ENDIF
2962
2963	sed_nr(nzt+1) = 0.0_wp
2964	sed_qr(nzt+1) = 0.0_wp
2965	!
2966	!-- Compute sedimentation flux
2967	DO k = nzt, nzb+1, -1
2968	!
2969	!-- Predetermine flag to mask topography
2970	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2971	!
2972	!-- Sum up all rain drop number densities which contribute to the flux
2973	!-- through k-1/2
2974	flux = 0.0_wp
2975	z_run = 0.0_wp ! height above z(k)
2976	k_run = k
2977	c_run = MIN( 1.0_wp, c_nr(k) )
2978	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
2979	flux = flux + hyrho(k_run) * &
2980	( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) * &
2981	0.5_wp ) * c_run * dzu(k_run) * flag
2982	z_run = z_run + dzu(k_run) * flag
2983	k_run = k_run + 1 * flag
2984	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) * flag
2985	ENDDO
2986	!
2987	!-- It is not allowed to sediment more rain drop number density than
2988	!-- available
2989	flux = MIN( flux, &
2990	hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
2991
2992	sed_nr(k) = flux / dt_micro * flag
2993	nr_1d(k) = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / &
2994	hyrho(k) * dt_micro * flag
2995	!
2996	!-- Sum up all rain water content which contributes to the flux
2997	!-- through k-1/2
2998	flux = 0.0_wp
2999	z_run = 0.0_wp ! height above z(k)
3000	k_run = k
3001	c_run = MIN( 1.0_wp, c_qr(k) )
3002
3003	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
3004
3005	flux = flux + hyrho(k_run) * &
3006	( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) * &
3007	0.5_wp ) * c_run * dzu(k_run) * flag
3008	z_run = z_run + dzu(k_run) * flag
3009	k_run = k_run + 1 * flag
3010	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) * flag
3011
3012	ENDDO
3013	!
3014	!-- It is not allowed to sediment more rain water content than available
3015	flux = MIN( flux, &
3016	hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
3017
3018	sed_qr(k) = flux / dt_micro * flag
3019
3020	qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
3021	hyrho(k) * dt_micro * flag
3022	q_1d(k) = q_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
3023	hyrho(k) * dt_micro * flag
3024	pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
3025	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro * flag
3026	!
3027	!-- Compute the rain rate
3028	IF ( call_microphysics_at_all_substeps ) THEN
3029	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) &
3030	* weight_substep(intermediate_timestep_count) * flag
3031	ELSE
3032	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
3033	ENDIF
3034
3035	ENDDO
3036
3037	END SUBROUTINE sedimentation_rain_ij
3038
3039
3040	!------------------------------------------------------------------------------!
3041	! Description:
3042	! ------------
3043	!> This subroutine computes the precipitation amount due to gravitational
3044	!> settling of rain and cloud (fog) droplets
3045	!------------------------------------------------------------------------------!
3046	SUBROUTINE calc_precipitation_amount_ij( i, j )
3047
3048	USE arrays_3d, &
3049	ONLY: precipitation_amount, prr
3050
3051	USE cloud_parameters, &
3052	ONLY: hyrho
3053
3054	USE control_parameters, &
3055	ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d, &
3056	intermediate_timestep_count, intermediate_timestep_count_max,&
3057	precipitation_amount_interval, time_do2d_xy
3058
3059	USE indices, &
3060	ONLY: nzb, nzt, wall_flags_0
3061
3062	USE kinds
3063
3064	USE surface_mod, &
3065	ONLY : bc_h
3066
3067	IMPLICIT NONE
3068
3069	INTEGER(iwp) :: i !< running index x direction
3070	INTEGER(iwp) :: j !< running index y direction
3071	INTEGER(iwp) :: k !< running index z direction
3072	INTEGER(iwp) :: m !< running index surface elements
3073	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
3074	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
3075
3076	IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
3077	( .NOT. call_microphysics_at_all_substeps .OR. &
3078	intermediate_timestep_count == intermediate_timestep_count_max ) ) &
3079	THEN
3080
3081	surf_s = bc_h(0)%start_index(j,i)
3082	surf_e = bc_h(0)%end_index(j,i)
3083	DO m = surf_s, surf_e
3084	k = bc_h(0)%k(m)
3085	precipitation_amount(j,i) = precipitation_amount(j,i) + &
3086	prr(k,j,i) * hyrho(k) * dt_3d
3087	ENDDO
3088
3089	ENDIF
3090
3091	END SUBROUTINE calc_precipitation_amount_ij
3092
3093	!------------------------------------------------------------------------------!
3094	! Description:
3095	! ------------
3096	!> This function computes the gamma function (Press et al., 1992).
3097	!> The gamma function is needed for the calculation of the evaporation
3098	!> of rain drops.
3099	!------------------------------------------------------------------------------!
3100	FUNCTION gamm( xx )
3101
3102	USE kinds
3103
3104	IMPLICIT NONE
3105
3106	INTEGER(iwp) :: j !<
3107
3108	REAL(wp) :: gamm !<
3109	REAL(wp) :: ser !<
3110	REAL(wp) :: tmp !<
3111	REAL(wp) :: x_gamm !<
3112	REAL(wp) :: xx !<
3113	REAL(wp) :: y_gamm !<
3114
3115
3116	REAL(wp), PARAMETER :: stp = 2.5066282746310005_wp !<
3117	REAL(wp), PARAMETER :: cof(6) = (/ 76.18009172947146_wp, &
3118	-86.50532032941677_wp, &
3119	24.01409824083091_wp, &
3120	-1.231739572450155_wp, &
3121	0.1208650973866179E-2_wp, &
3122	-0.5395239384953E-5_wp /) !<
3123
3124	x_gamm = xx
3125	y_gamm = x_gamm
3126	tmp = x_gamm + 5.5_wp
3127	tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
3128	ser = 1.000000000190015_wp
3129
3130	DO j = 1, 6
3131	y_gamm = y_gamm + 1.0_wp
3132	ser = ser + cof( j ) / y_gamm
3133	ENDDO
3134
3135	!
3136	!-- Until this point the algorithm computes the logarithm of the gamma
3137	!-- function. Hence, the exponential function is used.
3138	! gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
3139	gamm = EXP( tmp ) * stp * ser / x_gamm
3140
3141	RETURN
3142
3143	END FUNCTION gamm
3144
3145	END MODULE microphysics_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |