Home

Context Navigation

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

Last change on this file since 3137 was 3026, checked in by schwenkel, 7 years ago
Changed the name specific humidity to mixing ratio
Property svn:keywords set to `Id`
File size: 123.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |