Home

Context Navigation

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

Last change on this file since 2610 was 2608, checked in by schwenkel, 7 years ago
Inital revision of diagnostic_quantities_mod allows unified calculation of magnus equation and saturion mixing ratio
Property svn:keywords set to `Id`
File size: 122.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |