Home

Context Navigation

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

Last change on this file since 3254 was 3241, checked in by raasch, 6 years ago
various changes to avoid compiler warnings (mainly removal of unused variables)
Property svn:keywords set to `Id`
File size: 121.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |