Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |