Home

Context Navigation

bulk_cloud_model_mod.f90 @ 3622

Last change on this file since 3622 was 3622, checked in by schwenkel, 6 years ago
Important bugfix in case of restart runs
Property svn:keywords set to `Id`
File size: 162.5 KB

Line
1	!> @file bulk_cloud_model_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: bulk_cloud_model_mod.f90 3622 2018-12-12 09:52:53Z schwenkel $
27	! Important bugfix in case of restart runs.
28	!
29	! 3507 2018-11-08 14:28:46Z schwenkel
30	! Minor bugfixes for bcm cache version and initializing prr
31	!
32	! 3452 2018-10-30 13:13:34Z schwenkel
33	! Bugfix for profiles output
34	!
35	! 3445 2018-10-29 12:23:02Z schwenkel
36	! Minor bugfix and use of subroutine for supersaturation calculation in case
37	! of cache version
38	!
39	! 3383 2018-10-19 14:22:58Z knoop
40	! Modularization of all bulk cloud physics code components
41	!
42	! unused variables removed
43	!
44	! 3026 2018-05-22 10:30:53Z schwenkel
45	! Changed the name specific humidity to mixing ratio, since we are computing
46	! mixing ratios.
47	!
48	! 2718 2018-01-02 08:49:38Z maronga
49	! Corrected "Former revisions" section
50	!
51	! 2701 2017-12-15 15:40:50Z suehring
52	! Changes from last commit documented
53	!
54	! 2698 2017-12-14 18:46:24Z suehring
55	! Bugfix in get_topography_top_index
56	!
57	! 2696 2017-12-14 17:12:51Z kanani
58	! Change in file header (GPL part)
59	!
60	! 2608 2017-11-13 14:04:26Z schwenkel
61	! Calculation of supersaturation in external module (diagnostic_quantities_mod).
62	! Change: correct calculation of saturation specific humidity to saturation
63	! mixing ratio (the factor of 0.378 vanishes).
64	!
65	! 2522 2017-10-05 14:20:37Z schwenkel
66	! Minor bugfix
67	!
68	! 2375 2017-08-29 14:10:28Z schwenkel
69	! Improved aerosol initilization and some minor bugfixes
70	! for droplet sedimenation
71	!
72	! 2318 2017-07-20 17:27:44Z suehring
73	! Get topography top index via Function call
74	!
75	! 2317 2017-07-20 17:27:19Z suehring
76	! s1 changed to log_sigma
77	!
78	! 2292 2017-06-20 09:51:42Z schwenkel
79	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
80	! includes two more prognostic equations for cloud drop concentration (nc)
81	! and cloud water content (qc).
82	! - The process of activation is parameterized with a simple Twomey
83	! activion scheme or with considering solution and curvature
84	! effects (Khvorostyanov and Curry ,2006).
85	! - The saturation adjustment scheme is replaced by the parameterization
86	! of condensation rates (Khairoutdinov and Kogan, 2000, Mon. Wea. Rev.,128).
87	! - All other microphysical processes of Seifert and Beheng are used.
88	! Additionally, in those processes the reduction of cloud number concentration
89	! is considered.
90	!
91	! 2233 2017-05-30 18:08:54Z suehring
92	!
93	! 2232 2017-05-30 17:47:52Z suehring
94	! Adjustments to new topography and surface concept
95	!
96	! 2155 2017-02-21 09:57:40Z hoffmann
97	! Bugfix in the calculation of microphysical quantities on ghost points.
98	!
99	! 2031 2016-10-21 15:11:58Z knoop
100	! renamed variable rho to rho_ocean
101	!
102	! 2000 2016-08-20 18:09:15Z knoop
103	! Forced header and separation lines into 80 columns
104	!
105	! 1850 2016-04-08 13:29:27Z maronga
106	! Module renamed
107	! Adapted for modularization of microphysics.
108	!
109	! 1845 2016-04-08 08:29:13Z raasch
110	! nzb_2d replaced by nzb_s_inner, Kessler precipitation is stored at surface
111	! point (instead of one point above surface)
112	!
113	! 1831 2016-04-07 13:15:51Z hoffmann
114	! turbulence renamed collision_turbulence,
115	! drizzle renamed cloud_water_sedimentation. cloud_water_sedimentation also
116	! avaialble for microphysics_kessler.
117	!
118	! 1822 2016-04-07 07:49:42Z hoffmann
119	! Unused variables removed.
120	! Kessler scheme integrated.
121	!
122	! 1691 2015-10-26 16:17:44Z maronga
123	! Added new routine calc_precipitation_amount. The routine now allows to account
124	! for precipitation due to sedimenation of cloud (fog) droplets
125	!
126	! 1682 2015-10-07 23:56:08Z knoop
127	! Code annotations made doxygen readable
128	!
129	! 1646 2015-09-02 16:00:10Z hoffmann
130	! Bugfix: Wrong computation of d_mean.
131	!
132	! 1361 2014-04-16 15:17:48Z hoffmann
133	! Bugfix in sedimentation_rain: Index corrected.
134	! Vectorized version of adjust_cloud added.
135	! Little reformatting of the code.
136	!
137	! 1353 2014-04-08 15:21:23Z heinze
138	! REAL constants provided with KIND-attribute
139	!
140	! 1346 2014-03-27 13:18:20Z heinze
141	! Bugfix: REAL constants provided with KIND-attribute especially in call of
142	! intrinsic function like MAX, MIN, SIGN
143	!
144	! 1334 2014-03-25 12:21:40Z heinze
145	! Bugfix: REAL constants provided with KIND-attribute
146	!
147	! 1322 2014-03-20 16:38:49Z raasch
148	! REAL constants defined as wp-kind
149	!
150	! 1320 2014-03-20 08:40:49Z raasch
151	! ONLY-attribute added to USE-statements,
152	! kind-parameters added to all INTEGER and REAL declaration statements,
153	! kinds are defined in new module kinds,
154	! comment fields (!:) to be used for variable explanations added to
155	! all variable declaration statements
156	!
157	! 1241 2013-10-30 11:36:58Z heinze
158	! hyp and rho_ocean have to be calculated at each time step if data from external
159	! file LSF_DATA are used
160	!
161	! 1115 2013-03-26 18:16:16Z hoffmann
162	! microphyical tendencies are calculated in bcm_actions in an optimized
163	! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
164	!
165	! 1106 2013-03-04 05:31:38Z raasch
166	! small changes in code formatting
167	!
168	! 1092 2013-02-02 11:24:22Z raasch
169	! unused variables removed
170	! file put under GPL
171	!
172	! 1065 2012-11-22 17:42:36Z hoffmann
173	! Sedimentation process implemented according to Stevens and Seifert (2008).
174	! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
175	! and Stevens, 2010).
176	!
177	! 1053 2012-11-13 17:11:03Z hoffmann
178	! initial revision
179	!
180	! Description:
181	! ------------
182	!> Calculate bulk cloud microphysics.
183	!------------------------------------------------------------------------------!
184	MODULE bulk_cloud_model_mod
185
186	USE arrays_3d, &
187	#if defined (__nopointer)
188	ONLY: ddzu, diss, dzu, dzw, hyp, hyrho, &
189	nc, nc_p, nr, nr_p, &
190	precipitation_amount, prr, pt, d_exner, pt_init, q, ql, &
191	qc, qc_p, qr, qr_p, &
192	exner, zu, tnc_m, tnr_m, tqc_m, tqr_m
193	#else
194	ONLY: ddzu, diss, dzu, dzw, hyp, hyrho, &
195	nc, nc_1, nc_2, nc_3, nc_p, nr, nr_1, nr_2, nr_3, nr_p, &
196	precipitation_amount, prr, pt, d_exner, pt_init, q, ql, ql_1, &
197	qc, qc_1, qc_2, qc_3, qc_p, qr, qr_1, qr_2, qr_3, qr_p, &
198	exner, zu, tnc_m, tnr_m, tqc_m, tqr_m
199	#endif
200
201	USE averaging, &
202	ONLY: nc_av, nr_av, prr_av, qc_av, ql_av, qr_av
203
204	USE basic_constants_and_equations_mod, &
205	ONLY: c_p, g, lv_d_cp, lv_d_rd, l_v, magnus, molecular_weight_of_solute,&
206	molecular_weight_of_water, pi, rho_l, rho_s, r_d, r_v, vanthoff,&
207	exner_function, exner_function_invers, ideal_gas_law_rho, &
208	ideal_gas_law_rho_pt, barometric_formula, rd_d_rv
209
210	USE control_parameters, &
211	ONLY: dt_3d, dt_do2d_xy, intermediate_timestep_count, &
212	intermediate_timestep_count_max, large_scale_forcing, &
213	lsf_surf, pt_surface, rho_surface, surface_pressure, &
214	time_do2d_xy, message_string, initializing_actions
215
216	USE cpulog, &
217	ONLY: cpu_log, log_point_s
218
219	USE grid_variables, &
220	ONLY: dx, dy
221
222	USE indices, &
223	ONLY: nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb, nzt, &
224	wall_flags_0
225
226	USE kinds
227
228	USE statistics, &
229	ONLY: weight_pres, weight_substep
230
231	USE surface_mod, &
232	ONLY : bc_h, get_topography_top_index_ji, surf_bulk_cloud_model, &
233	surf_microphysics_morrison, surf_microphysics_seifert
234
235	IMPLICIT NONE
236
237	CHARACTER (LEN=20) :: aerosol_bulk = 'nacl' !< namelist parameter
238	CHARACTER (LEN=20) :: cloud_scheme = 'saturation_adjust' !< namelist parameter
239
240	LOGICAL :: aerosol_nacl =.TRUE. !< nacl aerosol for bulk scheme
241	LOGICAL :: aerosol_c3h4o4 =.FALSE. !< malonic acid aerosol for bulk scheme
242	LOGICAL :: aerosol_nh4no3 =.FALSE. !< malonic acid aerosol for bulk scheme
243
244	LOGICAL :: bulk_cloud_model = .FALSE. !< namelist parameter
245
246	LOGICAL :: cloud_water_sedimentation = .FALSE. !< cloud water sedimentation
247	LOGICAL :: curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory
248	LOGICAL :: limiter_sedimentation = .TRUE. !< sedimentation limiter
249	LOGICAL :: collision_turbulence = .FALSE. !< turbulence effects
250	LOGICAL :: ventilation_effect = .TRUE. !< ventilation effect
251
252	LOGICAL :: call_microphysics_at_all_substeps = .FALSE. !< namelist parameter
253	LOGICAL :: microphysics_sat_adjust = .FALSE. !< use saturation adjust bulk scheme?
254	LOGICAL :: microphysics_kessler = .FALSE. !< use kessler bulk scheme?
255	LOGICAL :: microphysics_morrison = .FALSE. !< use 2-moment Morrison (add. prog. eq. for nc and qc)
256	LOGICAL :: microphysics_seifert = .FALSE. !< use 2-moment Seifert and Beheng scheme
257	LOGICAL :: precipitation = .FALSE. !< namelist parameter
258
259	REAL(wp) :: precipitation_amount_interval = 9999999.9_wp !< namelist parameter
260
261	REAL(wp) :: a_1 = 8.69E-4_wp !< coef. in turb. parametrization (cm-2 s3)
262	REAL(wp) :: a_2 = -7.38E-5_wp !< coef. in turb. parametrization (cm-2 s3)
263	REAL(wp) :: a_3 = -1.40E-2_wp !< coef. in turb. parametrization
264	REAL(wp) :: a_term = 9.65_wp !< coef. for terminal velocity (m s-1)
265	REAL(wp) :: a_vent = 0.78_wp !< coef. for ventilation effect
266	REAL(wp) :: b_1 = 11.45E-6_wp !< coef. in turb. parametrization (m)
267	REAL(wp) :: b_2 = 9.68E-6_wp !< coef. in turb. parametrization (m)
268	REAL(wp) :: b_3 = 0.62_wp !< coef. in turb. parametrization
269	REAL(wp) :: b_term = 9.8_wp !< coef. for terminal velocity (m s-1)
270	REAL(wp) :: b_vent = 0.308_wp !< coef. for ventilation effect
271	REAL(wp) :: beta_cc = 3.09E-4_wp !< coef. in turb. parametrization (cm-2 s3)
272	REAL(wp) :: c_1 = 4.82E-6_wp !< coef. in turb. parametrization (m)
273	REAL(wp) :: c_2 = 4.8E-6_wp !< coef. in turb. parametrization (m)
274	REAL(wp) :: c_3 = 0.76_wp !< coef. in turb. parametrization
275	REAL(wp) :: c_const = 0.93_wp !< const. in Taylor-microscale Reynolds number
276	REAL(wp) :: c_evap = 0.7_wp !< constant in evaporation
277	REAL(wp) :: c_term = 600.0_wp !< coef. for terminal velocity (m-1)
278	REAL(wp) :: diff_coeff_l = 0.23E-4_wp !< diffusivity of water vapor (m2 s-1)
279	REAL(wp) :: eps_sb = 1.0E-10_wp !< threshold in two-moments scheme
280	REAL(wp) :: eps_mr = 0.0_wp !< threshold for morrison scheme
281	REAL(wp) :: k_cc = 9.44E09_wp !< const. cloud-cloud kernel (m3 kg-2 s-1)
282	REAL(wp) :: k_cr0 = 4.33_wp !< const. cloud-rain kernel (m3 kg-1 s-1)
283	REAL(wp) :: k_rr = 7.12_wp !< const. rain-rain kernel (m3 kg-1 s-1)
284	REAL(wp) :: k_br = 1000.0_wp !< const. in breakup parametrization (m-1)
285	REAL(wp) :: k_st = 1.2E8_wp !< const. in drizzle parametrization (m-1 s-1)
286	REAL(wp) :: kin_vis_air = 1.4086E-5_wp !< kin. viscosity of air (m2 s-1)
287	REAL(wp) :: prec_time_const = 0.001_wp !< coef. in Kessler scheme (s-1)
288	REAL(wp) :: ql_crit = 0.0005_wp !< coef. in Kessler scheme (kg kg-1)
289	REAL(wp) :: schmidt_p_1d3=0.8921121_wp !< Schmidt number0.33333, 0.710.33333
290	REAL(wp) :: sigma_gc = 1.3_wp !< geometric standard deviation cloud droplets
291	REAL(wp) :: thermal_conductivity_l = 2.43E-2_wp !< therm. cond. air (J m-1 s-1 K-1)
292	REAL(wp) :: w_precipitation = 9.65_wp !< maximum terminal velocity (m s-1)
293	REAL(wp) :: x0 = 2.6E-10_wp !< separating drop mass (kg)
294	REAL(wp) :: xamin = 5.24E-19_wp !< average aerosol mass (kg) (~ 0.05Âµm)
295	REAL(wp) :: xcmin = 4.18E-15_wp !< minimum cloud drop size (kg) (~ 1Âµm)
296	REAL(wp) :: xrmin = 2.6E-10_wp !< minimum rain drop size (kg)
297	REAL(wp) :: xrmax = 5.0E-6_wp !< maximum rain drop site (kg)
298
299	REAL(wp) :: c_sedimentation = 2.0_wp !< Courant number of sedimentation process
300	REAL(wp) :: dpirho_l !< 6.0 / ( pi * rho_l )
301	REAL(wp) :: dry_aerosol_radius = 0.05E-6_wp !< dry aerosol radius
302	REAL(wp) :: dt_micro !< microphysics time step
303	REAL(wp) :: sigma_bulk = 2.0_wp !< width of aerosol spectrum
304	REAL(wp) :: na_init = 100.0E6_wp !< Total particle/aerosol concentration (cm-3)
305	REAL(wp) :: nc_const = 70.0E6_wp !< cloud droplet concentration
306	REAL(wp) :: dt_precipitation = 100.0_wp !< timestep precipitation (s)
307	REAL(wp) :: sed_qc_const !< const. for sedimentation of cloud water
308	REAL(wp) :: pirho_l !< pi * rho_l / 6.0;
309
310	REAL(wp) :: e_s !< saturation water vapor pressure
311	REAL(wp) :: q_s !< saturation mixing ratio
312	REAL(wp) :: sat !< supersaturation
313	REAL(wp) :: t_l !< actual temperature
314
315	SAVE
316
317	PRIVATE
318
319	PUBLIC bcm_parin, &
320	bcm_check_parameters, &
321	bcm_check_data_output, &
322	bcm_check_data_output_pr, &
323	bcm_header, &
324	bcm_init_arrays, &
325	bcm_init, &
326	bcm_3d_data_averaging, &
327	bcm_data_output_2d, &
328	bcm_data_output_3d, &
329	bcm_swap_timelevel, &
330	bcm_rrd_global, &
331	bcm_rrd_local, &
332	bcm_wrd_global, &
333	bcm_wrd_local, &
334	bcm_actions, &
335	calc_liquid_water_content
336
337	PUBLIC call_microphysics_at_all_substeps, &
338	cloud_water_sedimentation, &
339	bulk_cloud_model, &
340	cloud_scheme, &
341	collision_turbulence, &
342	dt_precipitation, &
343	microphysics_morrison, &
344	microphysics_sat_adjust, &
345	microphysics_seifert, &
346	na_init, &
347	nc_const, &
348	precipitation, &
349	sigma_gc
350
351
352	INTERFACE bcm_parin
353	MODULE PROCEDURE bcm_parin
354	END INTERFACE bcm_parin
355
356	INTERFACE bcm_check_parameters
357	MODULE PROCEDURE bcm_check_parameters
358	END INTERFACE bcm_check_parameters
359
360	INTERFACE bcm_check_data_output
361	MODULE PROCEDURE bcm_check_data_output
362	END INTERFACE bcm_check_data_output
363
364	INTERFACE bcm_check_data_output_pr
365	MODULE PROCEDURE bcm_check_data_output_pr
366	END INTERFACE bcm_check_data_output_pr
367
368	INTERFACE bcm_header
369	MODULE PROCEDURE bcm_header
370	END INTERFACE bcm_header
371
372	INTERFACE bcm_init_arrays
373	MODULE PROCEDURE bcm_init_arrays
374	END INTERFACE bcm_init_arrays
375
376	INTERFACE bcm_init
377	MODULE PROCEDURE bcm_init
378	END INTERFACE bcm_init
379
380	INTERFACE bcm_3d_data_averaging
381	MODULE PROCEDURE bcm_3d_data_averaging
382	END INTERFACE bcm_3d_data_averaging
383
384	INTERFACE bcm_data_output_2d
385	MODULE PROCEDURE bcm_data_output_2d
386	END INTERFACE bcm_data_output_2d
387
388	INTERFACE bcm_data_output_3d
389	MODULE PROCEDURE bcm_data_output_3d
390	END INTERFACE bcm_data_output_3d
391
392	INTERFACE bcm_swap_timelevel
393	MODULE PROCEDURE bcm_swap_timelevel
394	END INTERFACE bcm_swap_timelevel
395
396	INTERFACE bcm_rrd_global
397	MODULE PROCEDURE bcm_rrd_global
398	END INTERFACE bcm_rrd_global
399
400	INTERFACE bcm_rrd_local
401	MODULE PROCEDURE bcm_rrd_local
402	END INTERFACE bcm_rrd_local
403
404	INTERFACE bcm_wrd_global
405	MODULE PROCEDURE bcm_wrd_global
406	END INTERFACE bcm_wrd_global
407
408	INTERFACE bcm_wrd_local
409	MODULE PROCEDURE bcm_wrd_local
410	END INTERFACE bcm_wrd_local
411
412	INTERFACE bcm_actions
413	MODULE PROCEDURE bcm_actions
414	MODULE PROCEDURE bcm_actions_ij
415	END INTERFACE bcm_actions
416
417	INTERFACE adjust_cloud
418	MODULE PROCEDURE adjust_cloud
419	MODULE PROCEDURE adjust_cloud_ij
420	END INTERFACE adjust_cloud
421
422	INTERFACE activation
423	MODULE PROCEDURE activation
424	MODULE PROCEDURE activation_ij
425	END INTERFACE activation
426
427	INTERFACE condensation
428	MODULE PROCEDURE condensation
429	MODULE PROCEDURE condensation_ij
430	END INTERFACE condensation
431
432	INTERFACE autoconversion
433	MODULE PROCEDURE autoconversion
434	MODULE PROCEDURE autoconversion_ij
435	END INTERFACE autoconversion
436
437	INTERFACE autoconversion_kessler
438	MODULE PROCEDURE autoconversion_kessler
439	MODULE PROCEDURE autoconversion_kessler_ij
440	END INTERFACE autoconversion_kessler
441
442	INTERFACE accretion
443	MODULE PROCEDURE accretion
444	MODULE PROCEDURE accretion_ij
445	END INTERFACE accretion
446
447	INTERFACE selfcollection_breakup
448	MODULE PROCEDURE selfcollection_breakup
449	MODULE PROCEDURE selfcollection_breakup_ij
450	END INTERFACE selfcollection_breakup
451
452	INTERFACE evaporation_rain
453	MODULE PROCEDURE evaporation_rain
454	MODULE PROCEDURE evaporation_rain_ij
455	END INTERFACE evaporation_rain
456
457	INTERFACE sedimentation_cloud
458	MODULE PROCEDURE sedimentation_cloud
459	MODULE PROCEDURE sedimentation_cloud_ij
460	END INTERFACE sedimentation_cloud
461
462	INTERFACE sedimentation_rain
463	MODULE PROCEDURE sedimentation_rain
464	MODULE PROCEDURE sedimentation_rain_ij
465	END INTERFACE sedimentation_rain
466
467	INTERFACE calc_precipitation_amount
468	MODULE PROCEDURE calc_precipitation_amount
469	MODULE PROCEDURE calc_precipitation_amount_ij
470	END INTERFACE calc_precipitation_amount
471
472	INTERFACE supersaturation
473	MODULE PROCEDURE supersaturation
474	END INTERFACE supersaturation
475
476	INTERFACE calc_liquid_water_content
477	MODULE PROCEDURE calc_liquid_water_content
478	END INTERFACE calc_liquid_water_content
479
480	CONTAINS
481
482
483	!------------------------------------------------------------------------------!
484	! Description:
485	! ------------
486	!> Parin for &bulk_cloud_parameters for the bulk cloud module
487	!------------------------------------------------------------------------------!
488	SUBROUTINE bcm_parin
489
490
491	IMPLICIT NONE
492
493	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
494
495	NAMELIST /bulk_cloud_parameters/ &
496	aerosol_bulk, &
497	c_sedimentation, &
498	call_microphysics_at_all_substeps, &
499	bulk_cloud_model, &
500	cloud_scheme, &
501	cloud_water_sedimentation, &
502	collision_turbulence, &
503	curvature_solution_effects_bulk, &
504	dry_aerosol_radius, &
505	limiter_sedimentation, &
506	na_init, &
507	nc_const, &
508	precipitation, &
509	precipitation_amount_interval, &
510	sigma_bulk, &
511	ventilation_effect
512
513	line = ' '
514	!
515	!-- Try to find bulk cloud module namelist
516	REWIND ( 11 )
517	line = ' '
518	DO WHILE ( INDEX( line, '&bulk_cloud_parameters' ) == 0 )
519	READ ( 11, '(A)', END=10 ) line
520	ENDDO
521	BACKSPACE ( 11 )
522	!
523	!-- Read user-defined namelist
524	READ ( 11, bulk_cloud_parameters )
525	!
526	!-- Set flag that indicates that the bulk cloud module is switched on
527	!bulk_cloud_model = .TRUE.
528
529	10 CONTINUE
530
531
532	END SUBROUTINE bcm_parin
533
534
535	!------------------------------------------------------------------------------!
536	! Description:
537	! ------------
538	!> Check parameters routine for bulk cloud module
539	!------------------------------------------------------------------------------!
540	SUBROUTINE bcm_check_parameters
541
542
543	IMPLICIT NONE
544	!
545	!-- Check cloud scheme
546	IF ( cloud_scheme == 'saturation_adjust' ) THEN
547	microphysics_sat_adjust = .TRUE.
548	microphysics_seifert = .FALSE.
549	microphysics_kessler = .FALSE.
550	precipitation = .FALSE.
551	ELSEIF ( cloud_scheme == 'seifert_beheng' ) THEN
552	microphysics_sat_adjust = .FALSE.
553	microphysics_seifert = .TRUE.
554	microphysics_kessler = .FALSE.
555	microphysics_morrison = .FALSE.
556	precipitation = .TRUE.
557	ELSEIF ( cloud_scheme == 'kessler' ) THEN
558	microphysics_sat_adjust = .FALSE.
559	microphysics_seifert = .FALSE.
560	microphysics_kessler = .TRUE.
561	microphysics_morrison = .FALSE.
562	precipitation = .TRUE.
563	ELSEIF ( cloud_scheme == 'morrison' ) THEN
564	microphysics_sat_adjust = .FALSE.
565	microphysics_seifert = .TRUE.
566	microphysics_kessler = .FALSE.
567	microphysics_morrison = .TRUE.
568	precipitation = .TRUE.
569	ELSE
570	message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
571	TRIM( cloud_scheme ) // '"'
572	CALL message( 'check_parameters', 'PA0357', 1, 2, 0, 6, 0 )
573	ENDIF
574
575
576
577	!
578	!-- Set the default value for the integration interval of precipitation amount
579	IF ( microphysics_seifert .OR. microphysics_kessler ) THEN
580	IF ( precipitation_amount_interval == 9999999.9_wp ) THEN
581	precipitation_amount_interval = dt_do2d_xy
582	ELSE
583	IF ( precipitation_amount_interval > dt_do2d_xy ) THEN
584	WRITE( message_string, * ) 'precipitation_amount_interval = ', &
585	precipitation_amount_interval, ' must not be larger than ', &
586	'dt_do2d_xy = ', dt_do2d_xy
587	CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 )
588	ENDIF
589	ENDIF
590	ENDIF
591
592	! TODO: find better sollution for circular dependency problem
593	surf_bulk_cloud_model = bulk_cloud_model
594	surf_microphysics_morrison = microphysics_morrison
595	surf_microphysics_seifert = microphysics_seifert
596
597	!
598	!-- Check aerosol
599	IF ( aerosol_bulk == 'nacl' ) THEN
600	aerosol_nacl = .TRUE.
601	aerosol_c3h4o4 = .FALSE.
602	aerosol_nh4no3 = .FALSE.
603	ELSEIF ( aerosol_bulk == 'c3h4o4' ) THEN
604	aerosol_nacl = .FALSE.
605	aerosol_c3h4o4 = .TRUE.
606	aerosol_nh4no3 = .FALSE.
607	ELSEIF ( aerosol_bulk == 'nh4no3' ) THEN
608	aerosol_nacl = .FALSE.
609	aerosol_c3h4o4 = .FALSE.
610	aerosol_nh4no3 = .TRUE.
611	ELSE
612	message_string = 'unknown aerosol = "' // TRIM( aerosol_bulk ) // '"'
613	CALL message( 'check_parameters', 'PA0469', 1, 2, 0, 6, 0 )
614	ENDIF
615
616
617	END SUBROUTINE bcm_check_parameters
618
619	!------------------------------------------------------------------------------!
620	! Description:
621	! ------------
622	!> Check data output for bulk cloud module
623	!------------------------------------------------------------------------------!
624	SUBROUTINE bcm_check_data_output( var, unit )
625
626	IMPLICIT NONE
627
628	CHARACTER (LEN=*) :: unit !<
629	CHARACTER (LEN=*) :: var !<
630
631	SELECT CASE ( TRIM( var ) )
632
633	CASE ( 'nc' )
634	IF ( .NOT. microphysics_morrison ) THEN
635	message_string = 'output of "' // TRIM( var ) // '" ' // &
636	'requires ' // &
637	'cloud_scheme = "morrison"'
638	CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
639	ENDIF
640	unit = '1/m3'
641
642	CASE ( 'nr' )
643	IF ( .NOT. microphysics_seifert ) THEN
644	message_string = 'output of "' // TRIM( var ) // '" ' // &
645	'requires ' // &
646	'cloud_scheme = "seifert_beheng"'
647	CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
648	ENDIF
649	unit = '1/m3'
650
651	CASE ( 'prr' )
652	IF ( microphysics_sat_adjust ) THEN
653	message_string = 'output of "' // TRIM( var ) // '" ' // &
654	'is not available for ' // &
655	'cloud_scheme = "saturation_adjust"'
656	CALL message( 'check_parameters', 'PA0423', 1, 2, 0, 6, 0 )
657	ENDIF
658	unit = 'kg/kg m/s'
659
660	CASE ( 'qc' )
661	unit = 'kg/kg'
662
663	CASE ( 'qr' )
664	IF ( .NOT. microphysics_seifert ) THEN
665	message_string = 'output of "' // TRIM( var ) // '" ' // &
666	'requires ' // &
667	'cloud_scheme = "seifert_beheng"'
668	CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
669	ENDIF
670	unit = 'kg/kg'
671
672	CASE ( 'pra*' )
673	IF ( .NOT. microphysics_kessler .AND. &
674	.NOT. microphysics_seifert ) THEN
675	message_string = 'output of "' // TRIM( var ) // '" ' // &
676	'requires ' // &
677	'cloud_scheme = "kessler" or "seifert_beheng"'
678	CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
679	ENDIF
680	! TODO: find sollution (maybe connected to flow_statistics redesign?)
681	! IF ( j == 1 ) THEN
682	! message_string = 'temporal averaging of precipitation ' // &
683	! 'amount "' // TRIM( var ) // '" is not possible'
684	! CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
685	! ENDIF
686	unit = 'mm'
687
688	CASE ( 'prr*' )
689	IF ( .NOT. microphysics_kessler .AND. &
690	.NOT. microphysics_seifert ) THEN
691	message_string = 'output of "' // TRIM( var ) // '"' // &
692	' requires' // &
693	' cloud_scheme = "kessler" or "seifert_beheng"'
694	CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
695	ENDIF
696	unit = 'mm/s'
697
698	CASE DEFAULT
699	unit = 'illegal'
700
701	END SELECT
702
703
704	END SUBROUTINE bcm_check_data_output
705
706
707	!------------------------------------------------------------------------------!
708	! Description:
709	! ------------
710	!> Check data output of profiles for bulk cloud module
711	!------------------------------------------------------------------------------!
712	SUBROUTINE bcm_check_data_output_pr( variable, var_count, unit, dopr_unit )
713
714	USE arrays_3d, &
715	ONLY: zu
716
717	USE control_parameters, &
718	ONLY: data_output_pr
719
720	USE profil_parameter, &
721	ONLY: dopr_index
722
723	USE statistics, &
724	ONLY: hom, statistic_regions, pr_palm
725
726	IMPLICIT NONE
727
728	CHARACTER (LEN=*) :: unit !<
729	CHARACTER (LEN=*) :: variable !<
730	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
731
732	INTEGER(iwp) :: var_count !<
733	INTEGER(iwp) :: pr_index !<
734
735	SELECT CASE ( TRIM( variable ) )
736
737	! TODO: make index generic: pr_index = pr_palm+1
738
739	CASE ( 'nc' )
740	IF ( .NOT. microphysics_morrison ) THEN
741	message_string = 'data_output_pr = ' // &
742	TRIM( data_output_pr(var_count) ) // &
743	' is not implemented for' // &
744	' cloud_scheme /= morrison'
745	CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
746	ENDIF
747	pr_index = 123
748	dopr_index(var_count) = pr_index
749	dopr_unit = '1/m3'
750	unit = dopr_unit
751	hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 )
752
753	CASE ( 'nr' )
754	IF ( .NOT. microphysics_seifert ) THEN
755	message_string = 'data_output_pr = ' // &
756	TRIM( data_output_pr(var_count) ) // &
757	' is not implemented for' // &
758	' cloud_scheme /= seifert_beheng'
759	CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
760	ENDIF
761	pr_index = 73
762	dopr_index(var_count) = pr_index
763	dopr_unit = '1/m3'
764	unit = dopr_unit
765	hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 )
766
767	CASE ( 'prr' )
768	IF ( microphysics_sat_adjust ) THEN
769	message_string = 'data_output_pr = ' // &
770	TRIM( data_output_pr(var_count) ) // &
771	' is not available for' // &
772	' cloud_scheme = saturation_adjust'
773	CALL message( 'check_parameters', 'PA0422', 1, 2, 0, 6, 0 )
774	ENDIF
775	pr_index = 76
776	dopr_index(var_count) = pr_index
777	dopr_unit = 'kg/kg m/s'
778	unit = dopr_unit
779	hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 )
780	CASE ( 'qc' )
781	pr_index = 75
782	dopr_index(var_count) = pr_index
783	dopr_unit = 'kg/kg'
784	unit = dopr_unit
785	hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 )
786
787	CASE ( 'qr' )
788	IF ( .NOT. microphysics_seifert ) THEN
789	message_string = 'data_output_pr = ' // &
790	TRIM( data_output_pr(var_count) ) // &
791	' is not implemented for' // &
792	' cloud_scheme /= seifert_beheng'
793	CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
794	ENDIF
795	pr_index = 74
796	dopr_index(var_count) = pr_index
797	dopr_unit = 'kg/kg'
798	unit = dopr_unit
799	hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 )
800
801	CASE DEFAULT
802	unit = 'illegal'
803
804	END SELECT
805
806	END SUBROUTINE bcm_check_data_output_pr
807
808
809	!------------------------------------------------------------------------------!
810	! Description:
811	! ------------
812	!> Allocate bulk cloud module arrays and define pointers
813	!------------------------------------------------------------------------------!
814	SUBROUTINE bcm_init_arrays
815
816	USE indices, &
817	ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt
818
819
820	IMPLICIT NONE
821
822	INTEGER(iwp) :: i !<
823	INTEGER(iwp) :: j !<
824	INTEGER(iwp) :: k !<
825	!
826	!-- Liquid water content
827	#if defined( __nopointer )
828	ALLOCATE ( ql(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
829	#else
830	ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
831	#endif
832
833	!
834	!-- 3D-cloud water content
835	IF ( .NOT. microphysics_morrison ) THEN
836	#if defined( __nopointer )
837	ALLOCATE( qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
838	#else
839	ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
840	#endif
841	ENDIF
842	!
843	!-- Precipitation amount and rate (only needed if output is switched)
844	ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg) )
845
846	!
847	!-- 3d-precipitation rate
848	ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
849
850	IF ( microphysics_morrison ) THEN
851	!
852	!-- 3D-cloud drop water content, cloud drop concentration arrays
853	#if defined( __nopointer )
854	ALLOCATE( nc(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
855	nc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
856	qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
857	qc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
858	tnc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
859	tqc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
860	#else
861	ALLOCATE( nc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
862	nc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
863	nc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
864	qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
865	qc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
866	qc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
867	#endif
868	ENDIF
869
870	IF ( microphysics_seifert ) THEN
871	!
872	!-- 3D-rain water content, rain drop concentration arrays
873	#if defined( __nopointer )
874	ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
875	nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
876	qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
877	qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
878	tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
879	tqr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
880	#else
881	ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
882	nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
883	nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
884	qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
885	qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
886	qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
887	#endif
888	ENDIF
889
890	#if ! defined( __nopointer )
891	!
892	!-- Initial assignment of the pointers
893	ql => ql_1
894	IF ( .NOT. microphysics_morrison ) THEN
895	qc => qc_1
896	ENDIF
897	IF ( microphysics_morrison ) THEN
898	qc => qc_1; qc_p => qc_2; tqc_m => qc_3
899	nc => nc_1; nc_p => nc_2; tnc_m => nc_3
900	ENDIF
901	IF ( microphysics_seifert ) THEN
902	qr => qr_1; qr_p => qr_2; tqr_m => qr_3
903	nr => nr_1; nr_p => nr_2; tnr_m => nr_3
904	ENDIF
905	#endif
906
907
908	END SUBROUTINE bcm_init_arrays
909
910
911	!------------------------------------------------------------------------------!
912	! Description:
913	! ------------
914	!> Initialization of the bulk cloud module
915	!------------------------------------------------------------------------------!
916	SUBROUTINE bcm_init !( dots_label, dots_unit, dots_num, dots_max )
917
918	IMPLICIT NONE
919
920	INTEGER(iwp) :: i !<
921	INTEGER(iwp) :: j !<
922	INTEGER(iwp) :: k !<
923
924	! INTEGER(iwp) :: dots_num
925	! INTEGER(iwp) :: dots_max
926	! CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_unit
927	! CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label
928
929	CALL location_message( 'initializing bulk cloud module', .FALSE. )
930
931	IF ( bulk_cloud_model ) THEN
932	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
933	!
934	!-- Initialize the remaining quantities
935	IF ( microphysics_morrison ) THEN
936	DO i = nxlg, nxrg
937	DO j = nysg, nyng
938	qc(:,j,i) = 0.0_wp
939	nc(:,j,i) = 0.0_wp
940	ENDDO
941	ENDDO
942	ENDIF
943
944	IF ( microphysics_seifert ) THEN
945	DO i = nxlg, nxrg
946	DO j = nysg, nyng
947	qr(:,j,i) = 0.0_wp
948	nr(:,j,i) = 0.0_wp
949	ENDDO
950	ENDDO
951	ENDIF
952	!
953	!-- Liquid water content and precipitation amount
954	!-- are zero at beginning of the simulation
955	ql = 0.0_wp
956	qc = 0.0_wp
957	precipitation_amount = 0.0_wp
958	prr = 0.0_wp
959	!
960	!-- Initialize old and new time levels.
961	IF ( microphysics_morrison ) THEN
962	tqc_m = 0.0_wp
963	tnc_m = 0.0_wp
964	qc_p = qc
965	nc_p = nc
966	ENDIF
967	IF ( microphysics_seifert ) THEN
968	tqr_m = 0.0_wp
969	tnr_m = 0.0_wp
970	qr_p = qr
971	nr_p = nr
972	ENDIF
973	ENDIF ! Only if not read_restart_data
974	!
975	!-- constant for the sedimentation of cloud water (2-moment cloud physics)
976	sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l ) &
977	)*( 2.0_wp / 3.0_wp ) &
978	EXP( 5.0_wp * LOG( sigma_gc )**2 )
979
980	!
981	!-- Calculate timestep according to precipitation
982	IF ( microphysics_seifert ) THEN
983	dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) / &
984	w_precipitation
985	ENDIF
986
987	!
988	!-- Set constants for certain aerosol type
989	IF ( microphysics_morrison ) THEN
990	IF ( aerosol_nacl ) THEN
991	molecular_weight_of_solute = 0.05844_wp
992	rho_s = 2165.0_wp
993	vanthoff = 2.0_wp
994	ELSEIF ( aerosol_c3h4o4 ) THEN
995	molecular_weight_of_solute = 0.10406_wp
996	rho_s = 1600.0_wp
997	vanthoff = 1.37_wp
998	ELSEIF ( aerosol_nh4no3 ) THEN
999	molecular_weight_of_solute = 0.08004_wp
1000	rho_s = 1720.0_wp
1001	vanthoff = 2.31_wp
1002	ENDIF
1003	ENDIF
1004
1005	!
1006	!-- Pre-calculate frequently calculated fractions of pi and rho_l
1007	pirho_l = pi * rho_l / 6.0_wp
1008	dpirho_l = 1.0_wp / pirho_l
1009
1010	CALL location_message( 'finished', .TRUE. )
1011
1012	ELSE
1013
1014	CALL location_message( 'skipped', .TRUE. )
1015
1016	ENDIF
1017
1018	END SUBROUTINE bcm_init
1019
1020
1021	!------------------------------------------------------------------------------!
1022	! Description:
1023	! ------------
1024	!> Swapping of timelevels
1025	!------------------------------------------------------------------------------!
1026	SUBROUTINE bcm_swap_timelevel ( mod_count )
1027
1028	IMPLICIT NONE
1029
1030	INTEGER, INTENT(IN) :: mod_count
1031
1032	IF ( bulk_cloud_model ) THEN
1033
1034	#if defined( __nopointer )
1035	IF ( microphysics_morrison ) THEN
1036	qc = qc_p
1037	nc = nc_p
1038	ENDIF
1039	IF ( microphysics_seifert ) THEN
1040	qr = qr_p
1041	nr = nr_p
1042	ENDIF
1043	#else
1044	SELECT CASE ( mod_count )
1045
1046	CASE ( 0 )
1047
1048	IF ( microphysics_morrison ) THEN
1049	qc => qc_1; qc_p => qc_2
1050	nc => nc_1; nc_p => nc_2
1051	ENDIF
1052	IF ( microphysics_seifert ) THEN
1053	qr => qr_1; qr_p => qr_2
1054	nr => nr_1; nr_p => nr_2
1055	ENDIF
1056
1057	CASE ( 1 )
1058
1059	IF ( microphysics_morrison ) THEN
1060	qc => qc_2; qc_p => qc_1
1061	nc => nc_2; nc_p => nc_1
1062	ENDIF
1063	IF ( microphysics_seifert ) THEN
1064	qr => qr_2; qr_p => qr_1
1065	nr => nr_2; nr_p => nr_1
1066	ENDIF
1067
1068	END SELECT
1069	#endif
1070
1071	ENDIF
1072
1073	END SUBROUTINE bcm_swap_timelevel
1074
1075
1076	!------------------------------------------------------------------------------!
1077	! Description:
1078	! ------------
1079	!> Header output for bulk cloud module
1080	!------------------------------------------------------------------------------!
1081	SUBROUTINE bcm_header ( io )
1082
1083
1084	IMPLICIT NONE
1085
1086	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
1087
1088	!
1089	!-- Write bulk cloud module header
1090	WRITE ( io, 1 )
1091
1092	WRITE ( io, 2 )
1093	WRITE ( io, 3 )
1094
1095	IF ( microphysics_kessler ) THEN
1096	WRITE ( io, 4 ) 'Kessler-Scheme'
1097	ENDIF
1098
1099	IF ( microphysics_seifert ) THEN
1100	WRITE ( io, 4 ) 'Seifert-Beheng-Scheme'
1101	IF ( cloud_water_sedimentation ) WRITE ( io, 5 )
1102	IF ( collision_turbulence ) WRITE ( io, 6 )
1103	IF ( ventilation_effect ) WRITE ( io, 7 )
1104	IF ( limiter_sedimentation ) WRITE ( io, 8 )
1105	ENDIF
1106
1107	WRITE ( io, 20 )
1108	WRITE ( io, 21 ) surface_pressure
1109	WRITE ( io, 22 ) r_d
1110	WRITE ( io, 23 ) rho_surface
1111	WRITE ( io, 24 ) c_p
1112	WRITE ( io, 25 ) l_v
1113
1114	IF ( microphysics_seifert ) THEN
1115	WRITE ( io, 26 ) 1.0E-6_wp * nc_const
1116	WRITE ( io, 27 ) c_sedimentation
1117	ENDIF
1118
1119
1120	1 FORMAT ( //' Bulk cloud module information:'/ &
1121	' ------------------------------------------'/ )
1122	2 FORMAT ( '--> Bulk scheme with liquid water potential temperature and'/ &
1123	' total water content is used.' )
1124	3 FORMAT ( '--> Condensation is parameterized via 0% - or 100% scheme.' )
1125	4 FORMAT ( '--> Precipitation parameterization via ', A )
1126
1127	5 FORMAT ( '--> Cloud water sedimentation parameterization via Stokes law' )
1128	6 FORMAT ( '--> Turbulence effects on precipitation process' )
1129	7 FORMAT ( '--> Ventilation effects on evaporation of rain drops' )
1130	8 FORMAT ( '--> Slope limiter used for sedimentation process' )
1131
1132	20 FORMAT ( '--> Essential parameters:' )
1133	21 FORMAT ( ' Surface pressure : p_0 = ', F7.2, ' hPa')
1134	22 FORMAT ( ' Gas constant : R = ', F5.1, ' J/(kg K)')
1135	23 FORMAT ( ' Density of air : rho_0 = ', F6.3, ' kg/m**3')
1136	24 FORMAT ( ' Specific heat cap. : c_p = ', F6.1, ' J/(kg K)')
1137	25 FORMAT ( ' Vapourization heat : L_v = ', E9.2, ' J/kg')
1138	26 FORMAT ( ' Droplet density : N_c = ', F6.1, ' 1/cm**3' )
1139	27 FORMAT ( ' Sedimentation Courant number : C_s = ', F4.1 )
1140
1141
1142	END SUBROUTINE bcm_header
1143
1144
1145	!------------------------------------------------------------------------------!
1146	!
1147	! Description:
1148	! ------------
1149	!> Subroutine for averaging 3D data
1150	!------------------------------------------------------------------------------!
1151	SUBROUTINE bcm_3d_data_averaging( mode, variable )
1152
1153	USE control_parameters, &
1154	ONLY: average_count_3d
1155
1156	IMPLICIT NONE
1157
1158	CHARACTER (LEN=*) :: mode !<
1159	CHARACTER (LEN=*) :: variable !<
1160
1161	INTEGER(iwp) :: i !< local index
1162	INTEGER(iwp) :: j !< local index
1163	INTEGER(iwp) :: k !< local index
1164
1165	IF ( mode == 'allocate' ) THEN
1166
1167	SELECT CASE ( TRIM( variable ) )
1168
1169	CASE ( 'nc' )
1170	IF ( .NOT. ALLOCATED( nc_av ) ) THEN
1171	ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1172	ENDIF
1173	nc_av = 0.0_wp
1174
1175	CASE ( 'nr' )
1176	IF ( .NOT. ALLOCATED( nr_av ) ) THEN
1177	ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1178	ENDIF
1179	nr_av = 0.0_wp
1180
1181	CASE ( 'prr' )
1182	IF ( .NOT. ALLOCATED( prr_av ) ) THEN
1183	ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1184	ENDIF
1185	prr_av = 0.0_wp
1186
1187	CASE ( 'qc' )
1188	IF ( .NOT. ALLOCATED( qc_av ) ) THEN
1189	ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1190	ENDIF
1191	qc_av = 0.0_wp
1192
1193	CASE ( 'ql' )
1194	IF ( .NOT. ALLOCATED( ql_av ) ) THEN
1195	ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1196	ENDIF
1197	ql_av = 0.0_wp
1198
1199	CASE ( 'qr' )
1200	IF ( .NOT. ALLOCATED( qr_av ) ) THEN
1201	ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1202	ENDIF
1203	qr_av = 0.0_wp
1204
1205	CASE DEFAULT
1206	CONTINUE
1207
1208	END SELECT
1209
1210	ELSEIF ( mode == 'sum' ) THEN
1211
1212	SELECT CASE ( TRIM( variable ) )
1213
1214	CASE ( 'nc' )
1215	IF ( ALLOCATED( nc_av ) ) THEN
1216	DO i = nxlg, nxrg
1217	DO j = nysg, nyng
1218	DO k = nzb, nzt+1
1219	nc_av(k,j,i) = nc_av(k,j,i) + nc(k,j,i)
1220	ENDDO
1221	ENDDO
1222	ENDDO
1223	ENDIF
1224
1225	CASE ( 'nr' )
1226	IF ( ALLOCATED( nr_av ) ) THEN
1227	DO i = nxlg, nxrg
1228	DO j = nysg, nyng
1229	DO k = nzb, nzt+1
1230	nr_av(k,j,i) = nr_av(k,j,i) + nr(k,j,i)
1231	ENDDO
1232	ENDDO
1233	ENDDO
1234	ENDIF
1235
1236	CASE ( 'prr' )
1237	IF ( ALLOCATED( prr_av ) ) THEN
1238	DO i = nxlg, nxrg
1239	DO j = nysg, nyng
1240	DO k = nzb, nzt+1
1241	prr_av(k,j,i) = prr_av(k,j,i) + prr(k,j,i)
1242	ENDDO
1243	ENDDO
1244	ENDDO
1245	ENDIF
1246
1247	CASE ( 'qc' )
1248	IF ( ALLOCATED( qc_av ) ) THEN
1249	DO i = nxlg, nxrg
1250	DO j = nysg, nyng
1251	DO k = nzb, nzt+1
1252	qc_av(k,j,i) = qc_av(k,j,i) + qc(k,j,i)
1253	ENDDO
1254	ENDDO
1255	ENDDO
1256	ENDIF
1257
1258	CASE ( 'ql' )
1259	IF ( ALLOCATED( ql_av ) ) THEN
1260	DO i = nxlg, nxrg
1261	DO j = nysg, nyng
1262	DO k = nzb, nzt+1
1263	ql_av(k,j,i) = ql_av(k,j,i) + ql(k,j,i)
1264	ENDDO
1265	ENDDO
1266	ENDDO
1267	ENDIF
1268
1269	CASE ( 'qr' )
1270	IF ( ALLOCATED( qr_av ) ) THEN
1271	DO i = nxlg, nxrg
1272	DO j = nysg, nyng
1273	DO k = nzb, nzt+1
1274	qr_av(k,j,i) = qr_av(k,j,i) + qr(k,j,i)
1275	ENDDO
1276	ENDDO
1277	ENDDO
1278	ENDIF
1279
1280	CASE DEFAULT
1281	CONTINUE
1282
1283	END SELECT
1284
1285	ELSEIF ( mode == 'average' ) THEN
1286
1287	SELECT CASE ( TRIM( variable ) )
1288
1289	CASE ( 'nc' )
1290	IF ( ALLOCATED( nc_av ) ) THEN
1291	DO i = nxlg, nxrg
1292	DO j = nysg, nyng
1293	DO k = nzb, nzt+1
1294	nc_av(k,j,i) = nc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1295	ENDDO
1296	ENDDO
1297	ENDDO
1298	ENDIF
1299
1300	CASE ( 'nr' )
1301	IF ( ALLOCATED( nr_av ) ) THEN
1302	DO i = nxlg, nxrg
1303	DO j = nysg, nyng
1304	DO k = nzb, nzt+1
1305	nr_av(k,j,i) = nr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1306	ENDDO
1307	ENDDO
1308	ENDDO
1309	ENDIF
1310
1311	CASE ( 'prr' )
1312	IF ( ALLOCATED( prr_av ) ) THEN
1313	DO i = nxlg, nxrg
1314	DO j = nysg, nyng
1315	DO k = nzb, nzt+1
1316	prr_av(k,j,i) = prr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1317	ENDDO
1318	ENDDO
1319	ENDDO
1320	ENDIF
1321
1322	CASE ( 'qc' )
1323	IF ( ALLOCATED( qc_av ) ) THEN
1324	DO i = nxlg, nxrg
1325	DO j = nysg, nyng
1326	DO k = nzb, nzt+1
1327	qc_av(k,j,i) = qc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1328	ENDDO
1329	ENDDO
1330	ENDDO
1331	ENDIF
1332
1333	CASE ( 'ql' )
1334	IF ( ALLOCATED( ql_av ) ) THEN
1335	DO i = nxlg, nxrg
1336	DO j = nysg, nyng
1337	DO k = nzb, nzt+1
1338	ql_av(k,j,i) = ql_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1339	ENDDO
1340	ENDDO
1341	ENDDO
1342	ENDIF
1343
1344	CASE ( 'qr' )
1345	IF ( ALLOCATED( qr_av ) ) THEN
1346	DO i = nxlg, nxrg
1347	DO j = nysg, nyng
1348	DO k = nzb, nzt+1
1349	qr_av(k,j,i) = qr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1350	ENDDO
1351	ENDDO
1352	ENDDO
1353	ENDIF
1354
1355	CASE DEFAULT
1356	CONTINUE
1357
1358	END SELECT
1359
1360	ENDIF
1361
1362	END SUBROUTINE bcm_3d_data_averaging
1363
1364
1365	!------------------------------------------------------------------------------!
1366	! Description:
1367	! ------------
1368	!> Define 2D output variables.
1369	!------------------------------------------------------------------------------!
1370	SUBROUTINE bcm_data_output_2d( av, variable, found, grid, mode, local_pf, &
1371	two_d, nzb_do, nzt_do )
1372
1373
1374	IMPLICIT NONE
1375
1376	CHARACTER (LEN=*), INTENT(INOUT) :: grid !< name of vertical grid
1377	CHARACTER (LEN=*), INTENT(IN) :: mode !< either 'xy', 'xz' or 'yz'
1378	CHARACTER (LEN=*), INTENT(IN) :: variable !< name of variable
1379
1380	INTEGER(iwp), INTENT(IN) :: av !< flag for (non-)average output
1381	INTEGER(iwp), INTENT(IN) :: nzb_do !< vertical output index (bottom)
1382	INTEGER(iwp), INTENT(IN) :: nzt_do !< vertical output index (top)
1383
1384	INTEGER(iwp) :: flag_nr !< number of masking flag
1385
1386	INTEGER(iwp) :: i !< loop index along x-direction
1387	INTEGER(iwp) :: j !< loop index along y-direction
1388	INTEGER(iwp) :: k !< loop index along z-direction
1389
1390	LOGICAL, INTENT(INOUT) :: found !< flag if output variable is found
1391	LOGICAL, INTENT(INOUT) :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
1392	LOGICAL :: resorted !< flag if output is already resorted
1393
1394	REAL(wp), PARAMETER :: fill_value = -999.0_wp !< value for the _FillValue attribute
1395
1396	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do), INTENT(INOUT) :: local_pf !< local
1397	!< array to which output data is resorted to
1398
1399	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable
1400
1401	found = .TRUE.
1402	resorted = .FALSE.
1403	!
1404	!-- Set masking flag for topography for not resorted arrays
1405	flag_nr = 0 ! 0 = scalar, 1 = u, 2 = v, 3 = w
1406
1407	SELECT CASE ( TRIM( variable ) )
1408
1409	CASE ( 'nc_xy', 'nc_xz', 'nc_yz' )
1410	IF ( av == 0 ) THEN
1411	to_be_resorted => nc
1412	ELSE
1413	IF ( .NOT. ALLOCATED( nc_av ) ) THEN
1414	ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1415	nc_av = REAL( fill_value, KIND = wp )
1416	ENDIF
1417	to_be_resorted => nc_av
1418	ENDIF
1419	IF ( mode == 'xy' ) grid = 'zu'
1420
1421	CASE ( 'nr_xy', 'nr_xz', 'nr_yz' )
1422	IF ( av == 0 ) THEN
1423	to_be_resorted => nr
1424	ELSE
1425	IF ( .NOT. ALLOCATED( nr_av ) ) THEN
1426	ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1427	nr_av = REAL( fill_value, KIND = wp )
1428	ENDIF
1429	to_be_resorted => nr_av
1430	ENDIF
1431	IF ( mode == 'xy' ) grid = 'zu'
1432
1433	CASE ( 'pra*_xy' ) ! 2d-array / integral quantity => no av
1434	! CALL exchange_horiz_2d( precipitation_amount )
1435	DO i = nxl, nxr
1436	DO j = nys, nyn
1437	local_pf(i,j,nzb+1) = precipitation_amount(j,i)
1438	ENDDO
1439	ENDDO
1440	precipitation_amount = 0.0_wp ! reset for next integ. interval
1441	resorted = .TRUE.
1442	two_d = .TRUE.
1443	IF ( mode == 'xy' ) grid = 'zu1'
1444
1445	CASE ( 'prr_xy', 'prr_xz', 'prr_yz' )
1446	IF ( av == 0 ) THEN
1447	! CALL exchange_horiz( prr, nbgp )
1448	DO i = nxl, nxr
1449	DO j = nys, nyn
1450	DO k = nzb, nzt+1
1451	local_pf(i,j,k) = prr(k,j,i) * hyrho(nzb+1)
1452	ENDDO
1453	ENDDO
1454	ENDDO
1455	ELSE
1456	IF ( .NOT. ALLOCATED( prr_av ) ) THEN
1457	ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1458	prr_av = REAL( fill_value, KIND = wp )
1459	ENDIF
1460	! CALL exchange_horiz( prr_av, nbgp )
1461	DO i = nxl, nxr
1462	DO j = nys, nyn
1463	DO k = nzb, nzt+1
1464	local_pf(i,j,k) = prr_av(k,j,i) * hyrho(nzb+1)
1465	ENDDO
1466	ENDDO
1467	ENDDO
1468	ENDIF
1469	resorted = .TRUE.
1470	IF ( mode == 'xy' ) grid = 'zu'
1471
1472	CASE ( 'qc_xy', 'qc_xz', 'qc_yz' )
1473	IF ( av == 0 ) THEN
1474	to_be_resorted => qc
1475	ELSE
1476	IF ( .NOT. ALLOCATED( qc_av ) ) THEN
1477	ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1478	qc_av = REAL( fill_value, KIND = wp )
1479	ENDIF
1480	to_be_resorted => qc_av
1481	ENDIF
1482	IF ( mode == 'xy' ) grid = 'zu'
1483
1484	CASE ( 'ql_xy', 'ql_xz', 'ql_yz' )
1485	IF ( av == 0 ) THEN
1486	to_be_resorted => ql
1487	ELSE
1488	IF ( .NOT. ALLOCATED( ql_av ) ) THEN
1489	ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1490	ql_av = REAL( fill_value, KIND = wp )
1491	ENDIF
1492	to_be_resorted => ql_av
1493	ENDIF
1494	IF ( mode == 'xy' ) grid = 'zu'
1495
1496	CASE ( 'qr_xy', 'qr_xz', 'qr_yz' )
1497	IF ( av == 0 ) THEN
1498	to_be_resorted => qr
1499	ELSE
1500	IF ( .NOT. ALLOCATED( qr_av ) ) THEN
1501	ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1502	qr_av = REAL( fill_value, KIND = wp )
1503	ENDIF
1504	to_be_resorted => qr_av
1505	ENDIF
1506	IF ( mode == 'xy' ) grid = 'zu'
1507
1508	CASE DEFAULT
1509	found = .FALSE.
1510	grid = 'none'
1511
1512	END SELECT
1513
1514	IF ( found .AND. .NOT. resorted ) THEN
1515	DO i = nxl, nxr
1516	DO j = nys, nyn
1517	DO k = nzb_do, nzt_do
1518	local_pf(i,j,k) = MERGE( &
1519	to_be_resorted(k,j,i), &
1520	REAL( fill_value, KIND = wp ), &
1521	BTEST( wall_flags_0(k,j,i), flag_nr ) &
1522	)
1523	ENDDO
1524	ENDDO
1525	ENDDO
1526	ENDIF
1527
1528	END SUBROUTINE bcm_data_output_2d
1529
1530
1531	!------------------------------------------------------------------------------!
1532	! Description:
1533	! ------------
1534	!> Define 3D output variables.
1535	!------------------------------------------------------------------------------!
1536	SUBROUTINE bcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
1537
1538
1539	IMPLICIT NONE
1540
1541	CHARACTER (LEN=*), INTENT(IN) :: variable !< name of variable
1542
1543	INTEGER(iwp), INTENT(IN) :: av !< flag for (non-)average output
1544	INTEGER(iwp), INTENT(IN) :: nzb_do !< lower limit of the data output (usually 0)
1545	INTEGER(iwp), INTENT(IN) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
1546
1547	INTEGER(iwp) :: flag_nr !< number of masking flag
1548
1549	INTEGER(iwp) :: i !< loop index along x-direction
1550	INTEGER(iwp) :: j !< loop index along y-direction
1551	INTEGER(iwp) :: k !< loop index along z-direction
1552
1553	LOGICAL, INTENT(INOUT) :: found !< flag if output variable is found
1554	LOGICAL :: resorted !< flag if output is already resorted
1555
1556	REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute
1557
1558	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do), INTENT(INOUT) :: local_pf !< local
1559	!< array to which output data is resorted to
1560
1561	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable
1562
1563	found = .TRUE.
1564	resorted = .FALSE.
1565	!
1566	!-- Set masking flag for topography for not resorted arrays
1567	flag_nr = 0 ! 0 = scalar, 1 = u, 2 = v, 3 = w
1568
1569	SELECT CASE ( TRIM( variable ) )
1570
1571	CASE ( 'nc' )
1572	IF ( av == 0 ) THEN
1573	to_be_resorted => nc
1574	ELSE
1575	IF ( .NOT. ALLOCATED( nc_av ) ) THEN
1576	ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1577	nc_av = REAL( fill_value, KIND = wp )
1578	ENDIF
1579	to_be_resorted => nc_av
1580	ENDIF
1581
1582	CASE ( 'nr' )
1583	IF ( av == 0 ) THEN
1584	to_be_resorted => nr
1585	ELSE
1586	IF ( .NOT. ALLOCATED( nr_av ) ) THEN
1587	ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1588	nr_av = REAL( fill_value, KIND = wp )
1589	ENDIF
1590	to_be_resorted => nr_av
1591	ENDIF
1592
1593	CASE ( 'prr' )
1594	IF ( av == 0 ) THEN
1595	DO i = nxl, nxr
1596	DO j = nys, nyn
1597	DO k = nzb_do, nzt_do
1598	local_pf(i,j,k) = prr(k,j,i)
1599	ENDDO
1600	ENDDO
1601	ENDDO
1602	ELSE
1603	IF ( .NOT. ALLOCATED( prr_av ) ) THEN
1604	ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1605	prr_av = REAL( fill_value, KIND = wp )
1606	ENDIF
1607	DO i = nxl, nxr
1608	DO j = nys, nyn
1609	DO k = nzb_do, nzt_do
1610	local_pf(i,j,k) = prr_av(k,j,i)
1611	ENDDO
1612	ENDDO
1613	ENDDO
1614	ENDIF
1615	resorted = .TRUE.
1616
1617	CASE ( 'qc' )
1618	IF ( av == 0 ) THEN
1619	to_be_resorted => qc
1620	ELSE
1621	IF ( .NOT. ALLOCATED( qc_av ) ) THEN
1622	ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1623	qc_av = REAL( fill_value, KIND = wp )
1624	ENDIF
1625	to_be_resorted => qc_av
1626	ENDIF
1627
1628	CASE ( 'ql' )
1629	IF ( av == 0 ) THEN
1630	to_be_resorted => ql
1631	ELSE
1632	IF ( .NOT. ALLOCATED( ql_av ) ) THEN
1633	ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1634	ql_av = REAL( fill_value, KIND = wp )
1635	ENDIF
1636	to_be_resorted => ql_av
1637	ENDIF
1638
1639	CASE ( 'qr' )
1640	IF ( av == 0 ) THEN
1641	to_be_resorted => qr
1642	ELSE
1643	IF ( .NOT. ALLOCATED( qr_av ) ) THEN
1644	ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1645	qr_av = REAL( fill_value, KIND = wp )
1646	ENDIF
1647	to_be_resorted => qr_av
1648	ENDIF
1649
1650	CASE DEFAULT
1651	found = .FALSE.
1652
1653	END SELECT
1654
1655
1656	IF ( found .AND. .NOT. resorted ) THEN
1657	DO i = nxl, nxr
1658	DO j = nys, nyn
1659	DO k = nzb_do, nzt_do
1660	local_pf(i,j,k) = MERGE( &
1661	to_be_resorted(k,j,i), &
1662	REAL( fill_value, KIND = wp ), &
1663	BTEST( wall_flags_0(k,j,i), flag_nr ) &
1664	)
1665	ENDDO
1666	ENDDO
1667	ENDDO
1668	ENDIF
1669
1670	END SUBROUTINE bcm_data_output_3d
1671
1672
1673	!------------------------------------------------------------------------------!
1674	! Description:
1675	! ------------
1676	!> This routine reads the respective restart data for the bulk cloud module.
1677	!------------------------------------------------------------------------------!
1678	SUBROUTINE bcm_rrd_global( found )
1679
1680
1681	USE control_parameters, &
1682	ONLY: length, restart_string
1683
1684
1685	IMPLICIT NONE
1686
1687	LOGICAL, INTENT(OUT) :: found
1688
1689
1690	found = .TRUE.
1691
1692	SELECT CASE ( restart_string(1:length) )
1693
1694	CASE ( 'c_sedimentation' )
1695	READ ( 13 ) c_sedimentation
1696
1697	CASE ( 'bulk_cloud_model' )
1698	READ ( 13 ) bulk_cloud_model
1699
1700	CASE ( 'cloud_scheme' )
1701	READ ( 13 ) cloud_scheme
1702
1703	CASE ( 'cloud_water_sedimentation' )
1704	READ ( 13 ) cloud_water_sedimentation
1705
1706	CASE ( 'collision_turbulence' )
1707	READ ( 13 ) collision_turbulence
1708
1709	CASE ( 'limiter_sedimentation' )
1710	READ ( 13 ) limiter_sedimentation
1711
1712	CASE ( 'nc_const' )
1713	READ ( 13 ) nc_const
1714
1715	CASE ( 'precipitation' )
1716	READ ( 13 ) precipitation
1717
1718	CASE ( 'ventilation_effect' )
1719	READ ( 13 ) ventilation_effect
1720
1721	CASE ( 'na_init' )
1722	READ ( 13 ) na_init
1723
1724	CASE ( 'dry_aerosol_radius' )
1725	READ ( 13 ) dry_aerosol_radius
1726
1727	CASE ( 'sigma_bulk' )
1728	READ ( 13 ) sigma_bulk
1729
1730	CASE ( 'aerosol_bulk' )
1731	READ ( 13 ) aerosol_bulk
1732
1733	CASE ( 'curvature_solution_effects_bulk' )
1734	READ ( 13 ) curvature_solution_effects_bulk
1735
1736
1737	! CASE ( 'global_paramter' )
1738	! READ ( 13 ) global_parameter
1739	! CASE ( 'global_array' )
1740	! IF ( .NOT. ALLOCATED( global_array ) ) ALLOCATE( global_array(1:10) )
1741	! READ ( 13 ) global_array
1742
1743	CASE DEFAULT
1744
1745	found = .FALSE.
1746
1747	END SELECT
1748
1749
1750	END SUBROUTINE bcm_rrd_global
1751
1752
1753	!------------------------------------------------------------------------------!
1754	! Description:
1755	! ------------
1756	!> This routine reads the respective restart data for the bulk cloud module.
1757	!------------------------------------------------------------------------------!
1758	SUBROUTINE bcm_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, &
1759	nxr_on_file, nynf, nync, nyn_on_file, nysf, &
1760	nysc, nys_on_file, tmp_2d, tmp_3d, found )
1761
1762
1763	USE control_parameters
1764
1765	USE indices
1766
1767	USE pegrid
1768
1769
1770	IMPLICIT NONE
1771
1772	INTEGER(iwp) :: i !<
1773	INTEGER(iwp) :: k !<
1774	INTEGER(iwp) :: nxlc !<
1775	INTEGER(iwp) :: nxlf !<
1776	INTEGER(iwp) :: nxl_on_file !<
1777	INTEGER(iwp) :: nxrc !<
1778	INTEGER(iwp) :: nxrf !<
1779	INTEGER(iwp) :: nxr_on_file !<
1780	INTEGER(iwp) :: nync !<
1781	INTEGER(iwp) :: nynf !<
1782	INTEGER(iwp) :: nyn_on_file !<
1783	INTEGER(iwp) :: nysc !<
1784	INTEGER(iwp) :: nysf !<
1785	INTEGER(iwp) :: nys_on_file !<
1786
1787	LOGICAL, INTENT(OUT) :: found
1788
1789	REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d !<
1790	REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !<
1791
1792	!
1793	!-- Here the reading of user-defined restart data follows:
1794	!-- Sample for user-defined output
1795
1796
1797	found = .TRUE.
1798
1799	SELECT CASE ( restart_string(1:length) )
1800
1801	CASE ( 'nc' )
1802	IF ( k == 1 ) READ ( 13 ) tmp_3d
1803	nc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1804	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1805
1806	CASE ( 'nc_av' )
1807	IF ( .NOT. ALLOCATED( nc_av ) ) THEN
1808	ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1809	ENDIF
1810	IF ( k == 1 ) READ ( 13 ) tmp_3d
1811	nc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1812	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1813
1814	CASE ( 'nr' )
1815	IF ( k == 1 ) READ ( 13 ) tmp_3d
1816	nr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1817	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1818
1819	CASE ( 'nr_av' )
1820	IF ( .NOT. ALLOCATED( nr_av ) ) THEN
1821	ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1822	ENDIF
1823	IF ( k == 1 ) READ ( 13 ) tmp_3d
1824	nr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1825	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1826
1827	CASE ( 'precipitation_amount' )
1828	IF ( k == 1 ) READ ( 13 ) tmp_2d
1829	precipitation_amount(nysc-nbgp:nync+nbgp, &
1830	nxlc-nbgp:nxrc+nbgp) = &
1831	tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1832
1833	CASE ( 'prr' )
1834	IF ( .NOT. ALLOCATED( prr ) ) THEN
1835	ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1836	ENDIF
1837	IF ( k == 1 ) READ ( 13 ) tmp_3d
1838	prr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1839	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1840
1841	CASE ( 'prr_av' )
1842	IF ( .NOT. ALLOCATED( prr_av ) ) THEN
1843	ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1844	ENDIF
1845	IF ( k == 1 ) READ ( 13 ) tmp_3d
1846	prr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1847	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1848
1849	CASE ( 'qc' )
1850	IF ( k == 1 ) READ ( 13 ) tmp_3d
1851	qc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1852	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1853
1854	CASE ( 'qc_av' )
1855	IF ( .NOT. ALLOCATED( qc_av ) ) THEN
1856	ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1857	ENDIF
1858	IF ( k == 1 ) READ ( 13 ) tmp_3d
1859	qc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1860	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1861
1862	CASE ( 'ql' )
1863	IF ( k == 1 ) READ ( 13 ) tmp_3d
1864	ql(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1865	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1866
1867	CASE ( 'ql_av' )
1868	IF ( .NOT. ALLOCATED( ql_av ) ) THEN
1869	ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1870	ENDIF
1871	IF ( k == 1 ) READ ( 13 ) tmp_3d
1872	ql_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1873	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1874
1875	CASE ( 'qr' )
1876	IF ( k == 1 ) READ ( 13 ) tmp_3d
1877	qr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1878	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1879
1880	CASE ( 'qr_av' )
1881	IF ( .NOT. ALLOCATED( qr_av ) ) THEN
1882	ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1883	ENDIF
1884	IF ( k == 1 ) READ ( 13 ) tmp_3d
1885	qr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1886	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1887	!
1888	CASE DEFAULT
1889
1890	found = .FALSE.
1891
1892	END SELECT
1893
1894
1895	END SUBROUTINE bcm_rrd_local
1896
1897
1898	!------------------------------------------------------------------------------!
1899	! Description:
1900	! ------------
1901	!> This routine writes the respective restart data for the bulk cloud module.
1902	!------------------------------------------------------------------------------!
1903	SUBROUTINE bcm_wrd_global
1904
1905
1906	IMPLICIT NONE
1907
1908	CALL wrd_write_string( 'c_sedimentation' )
1909	WRITE ( 14 ) c_sedimentation
1910
1911	CALL wrd_write_string( 'bulk_cloud_model' )
1912	WRITE ( 14 ) bulk_cloud_model
1913
1914	CALL wrd_write_string( 'cloud_scheme' )
1915	WRITE ( 14 ) cloud_scheme
1916
1917	CALL wrd_write_string( 'cloud_water_sedimentation' )
1918	WRITE ( 14 ) cloud_water_sedimentation
1919
1920	CALL wrd_write_string( 'collision_turbulence' )
1921	WRITE ( 14 ) collision_turbulence
1922
1923	CALL wrd_write_string( 'limiter_sedimentation' )
1924	WRITE ( 14 ) limiter_sedimentation
1925
1926	CALL wrd_write_string( 'nc_const' )
1927	WRITE ( 14 ) nc_const
1928
1929	CALL wrd_write_string( 'precipitation' )
1930	WRITE ( 14 ) precipitation
1931
1932	CALL wrd_write_string( 'ventilation_effect' )
1933	WRITE ( 14 ) ventilation_effect
1934
1935	CALL wrd_write_string( 'na_init' )
1936	WRITE ( 14 ) na_init
1937
1938	CALL wrd_write_string( 'dry_aerosol_radius' )
1939	WRITE ( 14 ) dry_aerosol_radius
1940
1941	CALL wrd_write_string( 'sigma_bulk' )
1942	WRITE ( 14 ) sigma_bulk
1943
1944	CALL wrd_write_string( 'aerosol_bulk' )
1945	WRITE ( 14 ) aerosol_bulk
1946
1947	CALL wrd_write_string( 'curvature_solution_effects_bulk' )
1948	WRITE ( 14 ) curvature_solution_effects_bulk
1949
1950
1951	! needs preceeding allocation if array
1952	! CALL wrd_write_string( 'global_parameter' )
1953	! WRITE ( 14 ) global_parameter
1954
1955	! IF ( ALLOCATED( inflow_damping_factor ) ) THEN
1956	! CALL wrd_write_string( 'inflow_damping_factor' )
1957	! WRITE ( 14 ) inflow_damping_factor
1958	! ENDIF
1959
1960
1961	END SUBROUTINE bcm_wrd_global
1962
1963
1964	!------------------------------------------------------------------------------!
1965	! Description:
1966	! ------------
1967	!> This routine writes the respective restart data for the bulk cloud module.
1968	!------------------------------------------------------------------------------!
1969	SUBROUTINE bcm_wrd_local
1970
1971
1972	IMPLICIT NONE
1973
1974	IF ( ALLOCATED( prr ) ) THEN
1975	CALL wrd_write_string( 'prr' )
1976	WRITE ( 14 ) prr
1977	ENDIF
1978
1979	IF ( ALLOCATED( prr_av ) ) THEN
1980	CALL wrd_write_string( 'prr_av' )
1981	WRITE ( 14 ) prr_av
1982	ENDIF
1983
1984	IF ( ALLOCATED( precipitation_amount ) ) THEN
1985	CALL wrd_write_string( 'precipitation_amount' )
1986	WRITE ( 14 ) precipitation_amount
1987	ENDIF
1988
1989	CALL wrd_write_string( 'ql' )
1990	WRITE ( 14 ) ql
1991
1992	IF ( ALLOCATED( ql_av ) ) THEN
1993	CALL wrd_write_string( 'ql_av' )
1994	WRITE ( 14 ) ql_av
1995	ENDIF
1996
1997	CALL wrd_write_string( 'qc' )
1998	WRITE ( 14 ) qc
1999
2000	IF ( ALLOCATED( qc_av ) ) THEN
2001	CALL wrd_write_string( 'qc_av' )
2002	WRITE ( 14 ) qc_av
2003	ENDIF
2004
2005	IF ( microphysics_morrison ) THEN
2006
2007	CALL wrd_write_string( 'nc' )
2008	WRITE ( 14 ) nc
2009
2010	IF ( ALLOCATED( nc_av ) ) THEN
2011	CALL wrd_write_string( 'nc_av' )
2012	WRITE ( 14 ) nc_av
2013	ENDIF
2014
2015	ENDIF
2016
2017	IF ( microphysics_seifert ) THEN
2018
2019	CALL wrd_write_string( 'nr' )
2020	WRITE ( 14 ) nr
2021
2022	IF ( ALLOCATED( nr_av ) ) THEN
2023	CALL wrd_write_string( 'nr_av' )
2024	WRITE ( 14 ) nr_av
2025	ENDIF
2026
2027	CALL wrd_write_string( 'qr' )
2028	WRITE ( 14 ) qr
2029
2030	IF ( ALLOCATED( qr_av ) ) THEN
2031	CALL wrd_write_string( 'qr_av' )
2032	WRITE ( 14 ) qr_av
2033	ENDIF
2034
2035	ENDIF
2036
2037
2038	END SUBROUTINE bcm_wrd_local
2039
2040
2041	!------------------------------------------------------------------------------!
2042	! Description:
2043	! ------------
2044	!> Control of microphysics for all grid points
2045	!------------------------------------------------------------------------------!
2046	SUBROUTINE bcm_actions
2047
2048	IMPLICIT NONE
2049
2050	IF ( large_scale_forcing .AND. lsf_surf ) THEN
2051	!
2052	!-- Calculate vertical profile of the hydrostatic pressure (hyp)
2053	hyp = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
2054	d_exner = exner_function_invers(hyp)
2055	exner = 1.0_wp / exner_function_invers(hyp)
2056	hyrho = ideal_gas_law_rho_pt(hyp, pt_init)
2057	!
2058	!-- Compute reference density
2059	rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
2060	ENDIF
2061
2062	!
2063	!-- Compute length of time step
2064	IF ( call_microphysics_at_all_substeps ) THEN
2065	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
2066	ELSE
2067	dt_micro = dt_3d
2068	ENDIF
2069
2070	!
2071	!-- Reset precipitation rate
2072	IF ( intermediate_timestep_count == 1 ) prr = 0.0_wp
2073
2074	!
2075	!-- Compute cloud physics
2076	IF ( microphysics_kessler ) THEN
2077
2078	CALL autoconversion_kessler
2079	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud
2080
2081	ELSEIF ( microphysics_seifert ) THEN
2082
2083	CALL adjust_cloud
2084	IF ( microphysics_morrison ) CALL activation
2085	IF ( microphysics_morrison ) CALL condensation
2086	CALL autoconversion
2087	CALL accretion
2088	CALL selfcollection_breakup
2089	CALL evaporation_rain
2090	CALL sedimentation_rain
2091	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud
2092
2093	ENDIF
2094
2095	CALL calc_precipitation_amount
2096
2097	END SUBROUTINE bcm_actions
2098
2099	!------------------------------------------------------------------------------!
2100	! Description:
2101	! ------------
2102	!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
2103	!> evaporation of rain drops due to too small or too big weights
2104	!> of rain drops (Stevens and Seifert, 2008).
2105	!------------------------------------------------------------------------------!
2106	SUBROUTINE adjust_cloud
2107
2108	IMPLICIT NONE
2109
2110	INTEGER(iwp) :: i !<
2111	INTEGER(iwp) :: j !<
2112	INTEGER(iwp) :: k !<
2113
2114	REAL(wp) :: flag !< flag to indicate first grid level above
2115
2116	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
2117
2118	DO i = nxlg, nxrg
2119	DO j = nysg, nyng
2120	DO k = nzb+1, nzt
2121	!
2122	!-- Predetermine flag to mask topography
2123	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2124
2125	IF ( qr(k,j,i) <= eps_sb ) THEN
2126	qr(k,j,i) = 0.0_wp
2127	nr(k,j,i) = 0.0_wp
2128	ELSE
2129	IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN
2130	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
2131	ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN
2132	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
2133	ENDIF
2134	ENDIF
2135
2136	IF ( microphysics_morrison ) THEN
2137	IF ( qc(k,j,i) <= eps_sb ) THEN
2138	qc(k,j,i) = 0.0_wp
2139	nc(k,j,i) = 0.0_wp
2140	ELSE
2141	IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN
2142	nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin * flag
2143	ENDIF
2144	ENDIF
2145	ENDIF
2146
2147	ENDDO
2148	ENDDO
2149	ENDDO
2150
2151	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
2152
2153	END SUBROUTINE adjust_cloud
2154
2155	!------------------------------------------------------------------------------!
2156	! Description:
2157	! ------------
2158	!> Calculate number of activated condensation nucleii after simple activation
2159	!> scheme of Twomey, 1959.
2160	!------------------------------------------------------------------------------!
2161	SUBROUTINE activation
2162
2163	IMPLICIT NONE
2164
2165	INTEGER(iwp) :: i !<
2166	INTEGER(iwp) :: j !<
2167	INTEGER(iwp) :: k !<
2168
2169	REAL(wp) :: activ !<
2170	REAL(wp) :: afactor !<
2171	REAL(wp) :: beta_act !<
2172	REAL(wp) :: bfactor !<
2173	REAL(wp) :: k_act !<
2174	REAL(wp) :: n_act !<
2175	REAL(wp) :: n_ccn !<
2176	REAL(wp) :: s_0 !<
2177	REAL(wp) :: sat_max !<
2178	REAL(wp) :: sigma !<
2179	REAL(wp) :: sigma_act !<
2180
2181	REAL(wp) :: flag !< flag to indicate first grid level above
2182
2183	CALL cpu_log( log_point_s(65), 'activation', 'start' )
2184
2185	DO i = nxlg, nxrg
2186	DO j = nysg, nyng
2187	DO k = nzb+1, nzt
2188	!
2189	!-- Predetermine flag to mask topography
2190	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2191
2192	!
2193	!-- Call calculation of supersaturation located in subroutine
2194	CALL supersaturation ( i, j, k )
2195	!
2196	!-- Prescribe parameters for activation
2197	!-- (see: Bott + Trautmann, 2002, Atm. Res., 64)
2198	k_act = 0.7_wp
2199	activ = 0.0_wp
2200
2201
2202	IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
2203	!
2204	!-- Compute the number of activated Aerosols
2205	!-- (see: Twomey, 1959, Pure and applied Geophysics, 43)
2206	n_act = na_init * sat**k_act
2207	!
2208	!-- Compute the number of cloud droplets
2209	!-- (see: Morrison + Grabowski, 2007, JAS, 64)
2210	! activ = MAX( n_act - nc(k,j,i), 0.0_wp) / dt_micro
2211
2212	!
2213	!-- Compute activation rate after Khairoutdinov and Kogan
2214	!-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
2215	sat_max = 1.0_wp / 100.0_wp
2216	activ = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN &
2217	( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) / &
2218	dt_micro
2219	ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk ) THEN
2220	!
2221	!-- Curvature effect (afactor) with surface tension
2222	!-- parameterization by Straka (2009)
2223	sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
2224	afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
2225	!
2226	!-- Solute effect (bfactor)
2227	bfactor = vanthoff * molecular_weight_of_water * &
2228	rho_s / ( molecular_weight_of_solute * rho_l )
2229
2230	!
2231	!-- Prescribe power index that describes the soluble fraction
2232	!-- of an aerosol particle (beta)
2233	!-- (see: Morrison + Grabowski, 2007, JAS, 64)
2234	beta_act = 0.5_wp
2235	sigma_act = sigma_bulk**( 1.0_wp + beta_act )
2236	!
2237	!-- Calculate mean geometric supersaturation (s_0) with
2238	!-- parameterization by Khvorostyanov and Curry (2006)
2239	s_0 = dry_aerosol_radius *(- ( 1.0_wp + beta_act ) ) &
2240	( 4.0_wp * afactor*3 / ( 27.0_wp bfactor ) )**0.5_wp
2241
2242	!
2243	!-- Calculate number of activated CCN as a function of
2244	!-- supersaturation and taking Koehler theory into account
2245	!-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
2246	n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( &
2247	LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
2248	activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp )
2249	ENDIF
2250
2251	nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro * flag), na_init)
2252
2253	ENDDO
2254	ENDDO
2255	ENDDO
2256
2257	CALL cpu_log( log_point_s(65), 'activation', 'stop' )
2258
2259	END SUBROUTINE activation
2260
2261
2262	!------------------------------------------------------------------------------!
2263	! Description:
2264	! ------------
2265	!> Calculate condensation rate for cloud water content (after Khairoutdinov and
2266	!> Kogan, 2000).
2267	!------------------------------------------------------------------------------!
2268	SUBROUTINE condensation
2269
2270	IMPLICIT NONE
2271
2272	INTEGER(iwp) :: i !<
2273	INTEGER(iwp) :: j !<
2274	INTEGER(iwp) :: k !<
2275
2276	REAL(wp) :: cond !<
2277	REAL(wp) :: cond_max !<
2278	REAL(wp) :: dc !<
2279	REAL(wp) :: evap !<
2280	REAL(wp) :: g_fac !<
2281	REAL(wp) :: nc_0 !<
2282	REAL(wp) :: temp !<
2283	REAL(wp) :: xc !<
2284
2285	REAL(wp) :: flag !< flag to indicate first grid level above
2286
2287	CALL cpu_log( log_point_s(66), 'condensation', 'start' )
2288
2289	DO i = nxlg, nxrg
2290	DO j = nysg, nyng
2291	DO k = nzb+1, nzt
2292	!
2293	!-- Predetermine flag to mask topography
2294	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2295	!
2296	!-- Call calculation of supersaturation
2297	CALL supersaturation ( i, j, k )
2298	!
2299	!-- Actual temperature:
2300	temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
2301
2302	g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
2303	l_v / ( thermal_conductivity_l * temp ) &
2304	+ r_v * temp / ( diff_coeff_l * e_s ) &
2305	)
2306	!
2307	!-- Mean weight of cloud drops
2308	IF ( nc(k,j,i) <= 0.0_wp) CYCLE
2309	xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
2310	!
2311	!-- Weight averaged diameter of cloud drops:
2312	dc = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
2313	!
2314	!-- Integral diameter of cloud drops
2315	nc_0 = nc(k,j,i) * dc
2316	!
2317	!-- Condensation needs only to be calculated in supersaturated regions
2318	IF ( sat > 0.0_wp ) THEN
2319	!
2320	!-- Condensation rate of cloud water content
2321	!-- after KK scheme.
2322	!-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
2323	cond = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
2324	cond_max = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
2325	cond = MIN( cond, cond_max / dt_micro )
2326
2327	qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag
2328	ELSEIF ( sat < 0.0_wp ) THEN
2329	evap = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
2330	evap = MAX( evap, -qc(k,j,i) / dt_micro )
2331
2332	qc(k,j,i) = qc(k,j,i) + evap * dt_micro * flag
2333	ENDIF
2334	IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN
2335	nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin
2336	ENDIF
2337	ENDDO
2338	ENDDO
2339	ENDDO
2340
2341	CALL cpu_log( log_point_s(66), 'condensation', 'stop' )
2342
2343	END SUBROUTINE condensation
2344
2345
2346	!------------------------------------------------------------------------------!
2347	! Description:
2348	! ------------
2349	!> Autoconversion rate (Seifert and Beheng, 2006).
2350	!------------------------------------------------------------------------------!
2351	SUBROUTINE autoconversion
2352
2353	IMPLICIT NONE
2354
2355	INTEGER(iwp) :: i !<
2356	INTEGER(iwp) :: j !<
2357	INTEGER(iwp) :: k !<
2358
2359	REAL(wp) :: alpha_cc !<
2360	REAL(wp) :: autocon !<
2361	REAL(wp) :: dissipation !<
2362	REAL(wp) :: flag !< flag to mask topography grid points
2363	REAL(wp) :: k_au !<
2364	REAL(wp) :: l_mix !<
2365	REAL(wp) :: nc_auto !<
2366	REAL(wp) :: nu_c !<
2367	REAL(wp) :: phi_au !<
2368	REAL(wp) :: r_cc !<
2369	REAL(wp) :: rc !<
2370	REAL(wp) :: re_lambda !<
2371	REAL(wp) :: sigma_cc !<
2372	REAL(wp) :: tau_cloud !<
2373	REAL(wp) :: xc !<
2374
2375	CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
2376
2377	DO i = nxlg, nxrg
2378	DO j = nysg, nyng
2379	DO k = nzb+1, nzt
2380	!
2381	!-- Predetermine flag to mask topography
2382	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2383
2384	IF ( microphysics_morrison ) THEN
2385	nc_auto = nc(k,j,i)
2386	ELSE
2387	nc_auto = nc_const
2388	ENDIF
2389
2390	IF ( qc(k,j,i) > eps_sb .AND. nc_auto > eps_mr ) THEN
2391
2392	k_au = k_cc / ( 20.0_wp * x0 )
2393	!
2394	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
2395	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
2396	tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + &
2397	qc(k,j,i) ), 0.0_wp )
2398	!
2399	!-- Universal function for autoconversion process
2400	!-- (Seifert and Beheng, 2006):
2401	phi_au = 600.0_wp * tau_cloud*0.68_wp &
2402	( 1.0_wp - tau_cloud0.68_wp )3
2403	!
2404	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
2405	!-- (Use constant nu_c = 1.0_wp instead?)
2406	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
2407	!
2408	!-- Mean weight of cloud droplets:
2409	xc = MAX( hyrho(k) * qc(k,j,i) / nc_auto, xcmin)
2410	!
2411	!-- Parameterized turbulence effects on autoconversion (Seifert,
2412	!-- Nuijens and Stevens, 2010)
2413	IF ( collision_turbulence ) THEN
2414	!
2415	!-- Weight averaged radius of cloud droplets:
2416	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
2417
2418	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
2419	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
2420	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
2421	!
2422	!-- Mixing length (neglecting distance to ground and
2423	!-- stratification)
2424	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
2425	!
2426	!-- Limit dissipation rate according to Seifert, Nuijens and
2427	!-- Stevens (2010)
2428	dissipation = MIN( 0.06_wp, diss(k,j,i) )
2429	!
2430	!-- Compute Taylor-microscale Reynolds number:
2431	re_lambda = 6.0_wp / 11.0_wp * &
2432	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
2433	SQRT( 15.0_wp / kin_vis_air ) * &
2434	dissipation**( 1.0_wp / 6.0_wp )
2435	!
2436	!-- The factor of 1.0E4 is needed to convert the dissipation
2437	!-- rate from m2 s-3 to cm2 s-3.
2438	k_au = k_au * ( 1.0_wp + &
2439	dissipation * 1.0E4_wp * &
2440	( re_lambda * 1.0E-3_wp )*0.25_wp &
2441	( alpha_cc * EXP( -1.0_wp * ( ( rc - &
2442	r_cc ) / &
2443	sigma_cc )**2 &
2444	) + beta_cc &
2445	) &
2446	)
2447	ENDIF
2448	!
2449	!-- Autoconversion rate (Seifert and Beheng, 2006):
2450	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
2451	( nu_c + 1.0_wp )*2 qc(k,j,i)*2 xc*2 &
2452	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
2453	rho_surface
2454	autocon = MIN( autocon, qc(k,j,i) / dt_micro )
2455
2456	qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag
2457	qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag
2458	nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro &
2459	* flag
2460	IF ( microphysics_morrison ) THEN
2461	nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp * &
2462	autocon / x0 * hyrho(k) * dt_micro * flag )
2463	ENDIF
2464
2465	ENDIF
2466
2467	ENDDO
2468	ENDDO
2469	ENDDO
2470
2471	CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
2472
2473	END SUBROUTINE autoconversion
2474
2475
2476	!------------------------------------------------------------------------------!
2477	! Description:
2478	! ------------
2479	!> Autoconversion process (Kessler, 1969).
2480	!------------------------------------------------------------------------------!
2481	SUBROUTINE autoconversion_kessler
2482
2483
2484	IMPLICIT NONE
2485
2486	INTEGER(iwp) :: i !<
2487	INTEGER(iwp) :: j !<
2488	INTEGER(iwp) :: k !<
2489	INTEGER(iwp) :: k_wall !< topgraphy top index
2490
2491	REAL(wp) :: dqdt_precip !<
2492	REAL(wp) :: flag !< flag to mask topography grid points
2493
2494	DO i = nxlg, nxrg
2495	DO j = nysg, nyng
2496	!
2497	!-- Determine vertical index of topography top
2498	k_wall = get_topography_top_index_ji( j, i, 's' )
2499	DO k = nzb+1, nzt
2500	!
2501	!-- Predetermine flag to mask topography
2502	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2503
2504	IF ( qc(k,j,i) > ql_crit ) THEN
2505	dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
2506	ELSE
2507	dqdt_precip = 0.0_wp
2508	ENDIF
2509
2510	qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
2511	q(k,j,i) = q(k,j,i) - dqdt_precip * dt_micro * flag
2512	pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp * &
2513	d_exner(k) * flag
2514
2515	!
2516	!-- Compute the rain rate (stored on surface grid point)
2517	prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
2518
2519	ENDDO
2520	ENDDO
2521	ENDDO
2522
2523	END SUBROUTINE autoconversion_kessler
2524
2525
2526	!------------------------------------------------------------------------------!
2527	! Description:
2528	! ------------
2529	!> Accretion rate (Seifert and Beheng, 2006).
2530	!------------------------------------------------------------------------------!
2531	SUBROUTINE accretion
2532
2533	IMPLICIT NONE
2534
2535	INTEGER(iwp) :: i !<
2536	INTEGER(iwp) :: j !<
2537	INTEGER(iwp) :: k !<
2538
2539	REAL(wp) :: accr !<
2540	REAL(wp) :: flag !< flag to mask topography grid points
2541	REAL(wp) :: k_cr !<
2542	REAL(wp) :: nc_accr !<
2543	REAL(wp) :: phi_ac !<
2544	REAL(wp) :: tau_cloud !<
2545	REAL(wp) :: xc !<
2546
2547
2548	CALL cpu_log( log_point_s(56), 'accretion', 'start' )
2549
2550	DO i = nxlg, nxrg
2551	DO j = nysg, nyng
2552	DO k = nzb+1, nzt
2553	!
2554	!-- Predetermine flag to mask topography
2555	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2556
2557	IF ( microphysics_morrison ) THEN
2558	nc_accr = nc(k,j,i)
2559	ELSE
2560	nc_accr = nc_const
2561	ENDIF
2562
2563	IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) &
2564	.AND. ( nc_accr > eps_mr ) ) THEN
2565	!
2566	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
2567	tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
2568	!
2569	!-- Universal function for accretion process (Seifert and
2570	!-- Beheng, 2001):
2571	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
2572
2573	!
2574	!-- Mean weight of cloud drops
2575	xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
2576	!
2577	!-- Parameterized turbulence effects on autoconversion (Seifert,
2578	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
2579	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
2580	IF ( collision_turbulence ) THEN
2581	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
2582	MIN( 600.0_wp, &
2583	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
2584	)
2585	ELSE
2586	k_cr = k_cr0
2587	ENDIF
2588	!
2589	!-- Accretion rate (Seifert and Beheng, 2006):
2590	accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * &
2591	SQRT( rho_surface * hyrho(k) )
2592	accr = MIN( accr, qc(k,j,i) / dt_micro )
2593
2594	qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
2595	qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
2596	IF ( microphysics_morrison ) THEN
2597	nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), &
2598	accr / xc * hyrho(k) * dt_micro * flag)
2599	ENDIF
2600
2601	ENDIF
2602
2603	ENDDO
2604	ENDDO
2605	ENDDO
2606
2607	CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
2608
2609	END SUBROUTINE accretion
2610
2611
2612	!------------------------------------------------------------------------------!
2613	! Description:
2614	! ------------
2615	!> Collisional breakup rate (Seifert, 2008).
2616	!------------------------------------------------------------------------------!
2617	SUBROUTINE selfcollection_breakup
2618
2619	IMPLICIT NONE
2620
2621	INTEGER(iwp) :: i !<
2622	INTEGER(iwp) :: j !<
2623	INTEGER(iwp) :: k !<
2624
2625	REAL(wp) :: breakup !<
2626	REAL(wp) :: dr !<
2627	REAL(wp) :: flag !< flag to mask topography grid points
2628	REAL(wp) :: phi_br !<
2629	REAL(wp) :: selfcoll !<
2630
2631	CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
2632
2633	DO i = nxlg, nxrg
2634	DO j = nysg, nyng
2635	DO k = nzb+1, nzt
2636	!
2637	!-- Predetermine flag to mask topography
2638	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2639
2640	IF ( qr(k,j,i) > eps_sb ) THEN
2641	!
2642	!-- Selfcollection rate (Seifert and Beheng, 2001):
2643	selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * &
2644	SQRT( hyrho(k) * rho_surface )
2645	!
2646	!-- Weight averaged diameter of rain drops:
2647	dr = ( hyrho(k) * qr(k,j,i) / &
2648	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
2649	!
2650	!-- Collisional breakup rate (Seifert, 2008):
2651	IF ( dr >= 0.3E-3_wp ) THEN
2652	phi_br = k_br * ( dr - 1.1E-3_wp )
2653	breakup = selfcoll * ( phi_br + 1.0_wp )
2654	ELSE
2655	breakup = 0.0_wp
2656	ENDIF
2657
2658	selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
2659	nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
2660
2661	ENDIF
2662	ENDDO
2663	ENDDO
2664	ENDDO
2665
2666	CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
2667
2668	END SUBROUTINE selfcollection_breakup
2669
2670
2671	!------------------------------------------------------------------------------!
2672	! Description:
2673	! ------------
2674	!> Evaporation of precipitable water. Condensation is neglected for
2675	!> precipitable water.
2676	!------------------------------------------------------------------------------!
2677	SUBROUTINE evaporation_rain
2678
2679	IMPLICIT NONE
2680
2681	INTEGER(iwp) :: i !<
2682	INTEGER(iwp) :: j !<
2683	INTEGER(iwp) :: k !<
2684
2685	REAL(wp) :: dr !<
2686	REAL(wp) :: evap !<
2687	REAL(wp) :: evap_nr !<
2688	REAL(wp) :: f_vent !<
2689	REAL(wp) :: flag !< flag to mask topography grid points
2690	REAL(wp) :: g_evap !<
2691	REAL(wp) :: lambda_r !<
2692	REAL(wp) :: mu_r !<
2693	REAL(wp) :: mu_r_2 !<
2694	REAL(wp) :: mu_r_5d2 !<
2695	REAL(wp) :: nr_0 !<
2696	REAL(wp) :: temp !<
2697	REAL(wp) :: xr !<
2698
2699	CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
2700
2701	DO i = nxlg, nxrg
2702	DO j = nysg, nyng
2703	DO k = nzb+1, nzt
2704	!
2705	!-- Predetermine flag to mask topography
2706	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2707
2708	IF ( qr(k,j,i) > eps_sb ) THEN
2709
2710	!
2711	!-- Call calculation of supersaturation
2712	CALL supersaturation ( i, j, k )
2713	!
2714	!-- Evaporation needs only to be calculated in subsaturated regions
2715	IF ( sat < 0.0_wp ) THEN
2716	!
2717	!-- Actual temperature:
2718	temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
2719
2720	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
2721	l_v / ( thermal_conductivity_l * temp ) &
2722	+ r_v * temp / ( diff_coeff_l * e_s ) &
2723	)
2724	!
2725	!-- Mean weight of rain drops
2726	xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
2727	!
2728	!-- Weight averaged diameter of rain drops:
2729	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
2730	!
2731	!-- Compute ventilation factor and intercept parameter
2732	!-- (Seifert and Beheng, 2006; Seifert, 2008):
2733	IF ( ventilation_effect ) THEN
2734	!
2735	!-- Shape parameter of gamma distribution (Milbrandt and Yau,
2736	!-- 2005; Stevens and Seifert, 2008):
2737	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
2738	( dr - 1.4E-3_wp ) ) )
2739	!
2740	!-- Slope parameter of gamma distribution (Seifert, 2008):
2741	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
2742	( mu_r + 1.0_wp ) &
2743	)**( 1.0_wp / 3.0_wp ) / dr
2744
2745	mu_r_2 = mu_r + 2.0_wp
2746	mu_r_5d2 = mu_r + 2.5_wp
2747
2748	f_vent = a_vent * gamm( mu_r_2 ) * &
2749	lambda_r*( -mu_r_2 ) + b_vent &
2750	schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
2751	gamm( mu_r_5d2 ) * lambda_r*( -mu_r_5d2 ) &
2752	( 1.0_wp - &
2753	0.5_wp * ( b_term / a_term ) * &
2754	( lambda_r / ( c_term + lambda_r ) &
2755	)**mu_r_5d2 - &
2756	0.125_wp * ( b_term / a_term )*2 &
2757	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
2758	)**mu_r_5d2 - &
2759	0.0625_wp * ( b_term / a_term )*3 &
2760	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
2761	)**mu_r_5d2 - &
2762	0.0390625_wp * ( b_term / a_term )*4 &
2763	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
2764	)**mu_r_5d2 &
2765	)
2766
2767	nr_0 = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) / &
2768	gamm( mu_r + 1.0_wp )
2769	ELSE
2770	f_vent = 1.0_wp
2771	nr_0 = nr(k,j,i) * dr
2772	ENDIF
2773	!
2774	!-- Evaporation rate of rain water content (Seifert and
2775	!-- Beheng, 2006):
2776	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / &
2777	hyrho(k)
2778	evap = MAX( evap, -qr(k,j,i) / dt_micro )
2779	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
2780	-nr(k,j,i) / dt_micro )
2781
2782	qr(k,j,i) = qr(k,j,i) + evap * dt_micro * flag
2783	nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag
2784
2785	ENDIF
2786	ENDIF
2787
2788	ENDDO
2789	ENDDO
2790	ENDDO
2791
2792	CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
2793
2794	END SUBROUTINE evaporation_rain
2795
2796
2797	!------------------------------------------------------------------------------!
2798	! Description:
2799	! ------------
2800	!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
2801	!------------------------------------------------------------------------------!
2802	SUBROUTINE sedimentation_cloud
2803
2804
2805	IMPLICIT NONE
2806
2807	INTEGER(iwp) :: i !<
2808	INTEGER(iwp) :: j !<
2809	INTEGER(iwp) :: k !<
2810
2811	REAL(wp) :: flag !< flag to mask topography grid points
2812	REAL(wp) :: nc_sedi !<
2813
2814	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
2815	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nc !<
2816
2817
2818	CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
2819
2820	sed_qc(nzt+1) = 0.0_wp
2821	sed_nc(nzt+1) = 0.0_wp
2822
2823	DO i = nxlg, nxrg
2824	DO j = nysg, nyng
2825	DO k = nzt, nzb+1, -1
2826	!
2827	!-- Predetermine flag to mask topography
2828	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2829
2830	IF ( microphysics_morrison ) THEN
2831	nc_sedi = nc(k,j,i)
2832	ELSE
2833	nc_sedi = nc_const
2834	ENDIF
2835
2836	!
2837	!-- Sedimentation fluxes for number concentration are only calculated
2838	!-- for cloud_scheme = 'morrison'
2839	IF ( microphysics_morrison ) THEN
2840	IF ( qc(k,j,i) > eps_sb .AND. nc(k,j,i) > eps_mr ) THEN
2841	sed_nc(k) = sed_qc_const * &
2842	( qc(k,j,i) * hyrho(k) )*( 2.0_wp / 3.0_wp ) &
2843	( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
2844	ELSE
2845	sed_nc(k) = 0.0_wp
2846	ENDIF
2847
2848	sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) * &
2849	nc(k,j,i) / dt_micro + sed_nc(k+1) &
2850	) * flag
2851
2852	nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) * &
2853	ddzu(k+1) / hyrho(k) * dt_micro * flag
2854	ENDIF
2855
2856	IF ( qc(k,j,i) > eps_sb .AND. nc_sedi > eps_mr ) THEN
2857	sed_qc(k) = sed_qc_const * nc_sedi*( -2.0_wp / 3.0_wp ) &
2858	( qc(k,j,i) * hyrho(k) )*( 5.0_wp / 3.0_wp ) &
2859	flag
2860	ELSE
2861	sed_qc(k) = 0.0_wp
2862	ENDIF
2863
2864	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / &
2865	dt_micro + sed_qc(k+1) &
2866	) * flag
2867
2868	q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
2869	ddzu(k+1) / hyrho(k) * dt_micro * flag
2870	qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
2871	ddzu(k+1) / hyrho(k) * dt_micro * flag
2872	pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * &
2873	ddzu(k+1) / hyrho(k) * lv_d_cp * &
2874	d_exner(k) * dt_micro * flag
2875
2876	!
2877	!-- Compute the precipitation rate due to cloud (fog) droplets
2878	IF ( call_microphysics_at_all_substeps ) THEN
2879	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) &
2880	* weight_substep(intermediate_timestep_count) &
2881	* flag
2882	ELSE
2883	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
2884	ENDIF
2885
2886	ENDDO
2887	ENDDO
2888	ENDDO
2889
2890	CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
2891
2892	END SUBROUTINE sedimentation_cloud
2893
2894
2895	!------------------------------------------------------------------------------!
2896	! Description:
2897	! ------------
2898	!> Computation of sedimentation flux. Implementation according to Stevens
2899	!> and Seifert (2008). Code is based on UCLA-LES.
2900	!------------------------------------------------------------------------------!
2901	SUBROUTINE sedimentation_rain
2902
2903	IMPLICIT NONE
2904
2905	INTEGER(iwp) :: i !< running index x direction
2906	INTEGER(iwp) :: j !< running index y direction
2907	INTEGER(iwp) :: k !< running index z direction
2908	INTEGER(iwp) :: k_run !<
2909	INTEGER(iwp) :: m !< running index surface elements
2910	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
2911	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
2912
2913	REAL(wp) :: c_run !<
2914	REAL(wp) :: d_max !<
2915	REAL(wp) :: d_mean !<
2916	REAL(wp) :: d_min !<
2917	REAL(wp) :: dr !<
2918	REAL(wp) :: flux !<
2919	REAL(wp) :: flag !< flag to mask topography grid points
2920	REAL(wp) :: lambda_r !<
2921	REAL(wp) :: mu_r !<
2922	REAL(wp) :: z_run !<
2923
2924	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<
2925	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<
2926	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<
2927	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<
2928	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<
2929	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<
2930	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<
2931	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<
2932
2933	CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
2934
2935	!
2936	!-- Compute velocities
2937	DO i = nxlg, nxrg
2938	DO j = nysg, nyng
2939	DO k = nzb+1, nzt
2940	!
2941	!-- Predetermine flag to mask topography
2942	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2943
2944	IF ( qr(k,j,i) > eps_sb ) THEN
2945	!
2946	!-- Weight averaged diameter of rain drops:
2947	dr = ( hyrho(k) * qr(k,j,i) / &
2948	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
2949	!
2950	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
2951	!-- Stevens and Seifert, 2008):
2952	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
2953	( dr - 1.4E-3_wp ) ) )
2954	!
2955	!-- Slope parameter of gamma distribution (Seifert, 2008):
2956	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
2957	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
2958
2959	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
2960	a_term - b_term * ( 1.0_wp + &
2961	c_term / &
2962	lambda_r )*( -1.0_wp &
2963	( mu_r + 1.0_wp ) ) &
2964	) &
2965	) * flag
2966
2967	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
2968	a_term - b_term * ( 1.0_wp + &
2969	c_term / &
2970	lambda_r )*( -1.0_wp &
2971	( mu_r + 4.0_wp ) ) &
2972	) &
2973	) * flag
2974	ELSE
2975	w_nr(k) = 0.0_wp
2976	w_qr(k) = 0.0_wp
2977	ENDIF
2978	ENDDO
2979	!
2980	!-- Adjust boundary values using surface data type.
2981	!-- Upward-facing
2982	surf_s = bc_h(0)%start_index(j,i)
2983	surf_e = bc_h(0)%end_index(j,i)
2984	DO m = surf_s, surf_e
2985	k = bc_h(0)%k(m)
2986	w_nr(k-1) = w_nr(k)
2987	w_qr(k-1) = w_qr(k)
2988	ENDDO
2989	!
2990	!-- Downward-facing
2991	surf_s = bc_h(1)%start_index(j,i)
2992	surf_e = bc_h(1)%end_index(j,i)
2993	DO m = surf_s, surf_e
2994	k = bc_h(1)%k(m)
2995	w_nr(k+1) = w_nr(k)
2996	w_qr(k+1) = w_qr(k)
2997	ENDDO
2998	!
2999	!-- Model top boundary value
3000	w_nr(nzt+1) = 0.0_wp
3001	w_qr(nzt+1) = 0.0_wp
3002	!
3003	!-- Compute Courant number
3004	DO k = nzb+1, nzt
3005	!
3006	!-- Predetermine flag to mask topography
3007	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3008
3009	c_nr(k) = 0.25_wp * ( w_nr(k-1) + &
3010	2.0_wp * w_nr(k) + w_nr(k+1) ) * &
3011	dt_micro * ddzu(k) * flag
3012	c_qr(k) = 0.25_wp * ( w_qr(k-1) + &
3013	2.0_wp * w_qr(k) + w_qr(k+1) ) * &
3014	dt_micro * ddzu(k) * flag
3015	ENDDO
3016	!
3017	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
3018	IF ( limiter_sedimentation ) THEN
3019
3020	DO k = nzb+1, nzt
3021	!
3022	!-- Predetermine flag to mask topography
3023	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3024
3025	d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
3026	d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
3027	d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
3028
3029	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
3030	2.0_wp * d_max, &
3031	ABS( d_mean ) ) &
3032	* flag
3033
3034	d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
3035	d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
3036	d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
3037
3038	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
3039	2.0_wp * d_max, &
3040	ABS( d_mean ) )
3041	ENDDO
3042
3043	ELSE
3044
3045	nr_slope = 0.0_wp
3046	qr_slope = 0.0_wp
3047
3048	ENDIF
3049
3050	sed_nr(nzt+1) = 0.0_wp
3051	sed_qr(nzt+1) = 0.0_wp
3052	!
3053	!-- Compute sedimentation flux
3054	DO k = nzt, nzb+1, -1
3055	!
3056	!-- Predetermine flag to mask topography
3057	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3058	!
3059	!-- Sum up all rain drop number densities which contribute to the flux
3060	!-- through k-1/2
3061	flux = 0.0_wp
3062	z_run = 0.0_wp ! height above z(k)
3063	k_run = k
3064	c_run = MIN( 1.0_wp, c_nr(k) )
3065	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
3066	flux = flux + hyrho(k_run) * &
3067	( nr(k_run,j,i) + nr_slope(k_run) * &
3068	( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run) &
3069	* flag
3070	z_run = z_run + dzu(k_run) * flag
3071	k_run = k_run + 1 * flag
3072	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) &
3073	* flag
3074	ENDDO
3075	!
3076	!-- It is not allowed to sediment more rain drop number density than
3077	!-- available
3078	flux = MIN( flux, &
3079	hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * &
3080	dt_micro &
3081	)
3082
3083	sed_nr(k) = flux / dt_micro * flag
3084	nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * &
3085	ddzu(k+1) / hyrho(k) * dt_micro * flag
3086	!
3087	!-- Sum up all rain water content which contributes to the flux
3088	!-- through k-1/2
3089	flux = 0.0_wp
3090	z_run = 0.0_wp ! height above z(k)
3091	k_run = k
3092	c_run = MIN( 1.0_wp, c_qr(k) )
3093
3094	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
3095
3096	flux = flux + hyrho(k_run) * ( qr(k_run,j,i) + &
3097	qr_slope(k_run) * ( 1.0_wp - c_run ) * &
3098	0.5_wp ) * c_run * dzu(k_run) * flag
3099	z_run = z_run + dzu(k_run) * flag
3100	k_run = k_run + 1 * flag
3101	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) &
3102	* flag
3103
3104	ENDDO
3105	!
3106	!-- It is not allowed to sediment more rain water content than
3107	!-- available
3108	flux = MIN( flux, &
3109	hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * &
3110	dt_micro &
3111	)
3112
3113	sed_qr(k) = flux / dt_micro * flag
3114
3115	qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
3116	ddzu(k+1) / hyrho(k) * dt_micro * flag
3117	q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
3118	ddzu(k+1) / hyrho(k) * dt_micro * flag
3119	pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * &
3120	ddzu(k+1) / hyrho(k) * lv_d_cp * &
3121	d_exner(k) * dt_micro * flag
3122	!
3123	!-- Compute the rain rate
3124	IF ( call_microphysics_at_all_substeps ) THEN
3125	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) &
3126	* weight_substep(intermediate_timestep_count) &
3127	* flag
3128	ELSE
3129	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
3130	ENDIF
3131
3132	ENDDO
3133	ENDDO
3134	ENDDO
3135
3136	CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
3137
3138	END SUBROUTINE sedimentation_rain
3139
3140
3141	!------------------------------------------------------------------------------!
3142	! Description:
3143	! ------------
3144	!> Computation of the precipitation amount due to gravitational settling of
3145	!> rain and cloud (fog) droplets
3146	!------------------------------------------------------------------------------!
3147	SUBROUTINE calc_precipitation_amount
3148
3149	IMPLICIT NONE
3150
3151	INTEGER(iwp) :: i !< running index x direction
3152	INTEGER(iwp) :: j !< running index y direction
3153	INTEGER(iwp) :: k !< running index y direction
3154	INTEGER(iwp) :: m !< running index surface elements
3155
3156	IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
3157	( .NOT. call_microphysics_at_all_substeps .OR. &
3158	intermediate_timestep_count == intermediate_timestep_count_max ) ) &
3159	THEN
3160	!
3161	!-- Run over all upward-facing surface elements, i.e. non-natural,
3162	!-- natural and urban
3163	DO m = 1, bc_h(0)%ns
3164	i = bc_h(0)%i(m)
3165	j = bc_h(0)%j(m)
3166	k = bc_h(0)%k(m)
3167	precipitation_amount(j,i) = precipitation_amount(j,i) + &
3168	prr(k,j,i) * hyrho(k) * dt_3d
3169	ENDDO
3170
3171	ENDIF
3172
3173	END SUBROUTINE calc_precipitation_amount
3174
3175
3176	!------------------------------------------------------------------------------!
3177	! Description:
3178	! ------------
3179	!> Control of microphysics for grid points i,j
3180	!------------------------------------------------------------------------------!
3181
3182	SUBROUTINE bcm_actions_ij( i, j )
3183
3184	IMPLICIT NONE
3185
3186	INTEGER(iwp) :: i !<
3187	INTEGER(iwp) :: j !<
3188
3189	IF ( large_scale_forcing .AND. lsf_surf ) THEN
3190	!
3191	!-- Calculate vertical profile of the hydrostatic pressure (hyp)
3192	hyp = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
3193	d_exner = exner_function_invers(hyp)
3194	exner = 1.0_wp / exner_function_invers(hyp)
3195	hyrho = ideal_gas_law_rho_pt(hyp, pt_init)
3196	!
3197	!-- Compute reference density
3198	rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
3199	ENDIF
3200
3201	!
3202	!-- Compute length of time step
3203	IF ( call_microphysics_at_all_substeps ) THEN
3204	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
3205	ELSE
3206	dt_micro = dt_3d
3207	ENDIF
3208	!
3209	!-- Reset precipitation rate
3210	IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp
3211
3212	!
3213	!-- Compute cloud physics
3214	IF( microphysics_kessler ) THEN
3215
3216	CALL autoconversion_kessler( i,j )
3217	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud( i,j )
3218
3219	ELSEIF ( microphysics_seifert ) THEN
3220
3221	CALL adjust_cloud( i,j )
3222	IF ( microphysics_morrison ) CALL activation( i,j )
3223	IF ( microphysics_morrison ) CALL condensation( i,j )
3224	CALL autoconversion( i,j )
3225	CALL accretion( i,j )
3226	CALL selfcollection_breakup( i,j )
3227	CALL evaporation_rain( i,j )
3228	CALL sedimentation_rain( i,j )
3229	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud( i,j )
3230
3231	ENDIF
3232
3233	CALL calc_precipitation_amount( i,j )
3234
3235	END SUBROUTINE bcm_actions_ij
3236
3237	!------------------------------------------------------------------------------!
3238	! Description:
3239	! ------------
3240	!> Adjust number of raindrops to avoid nonlinear effects in
3241	!> sedimentation and evaporation of rain drops due to too small or
3242	!> too big weights of rain drops (Stevens and Seifert, 2008).
3243	!> The same procedure is applied to cloud droplets if they are determined
3244	!> prognostically. Call for grid point i,j
3245	!------------------------------------------------------------------------------!
3246	SUBROUTINE adjust_cloud_ij( i, j )
3247
3248	IMPLICIT NONE
3249
3250	INTEGER(iwp) :: i !<
3251	INTEGER(iwp) :: j !<
3252	INTEGER(iwp) :: k !<
3253
3254	REAL(wp) :: flag !< flag to indicate first grid level above surface
3255
3256	DO k = nzb+1, nzt
3257	!
3258	!-- Predetermine flag to mask topography
3259	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3260
3261	IF ( qr(k,j,i) <= eps_sb ) THEN
3262	qr(k,j,i) = 0.0_wp
3263	nr(k,j,i) = 0.0_wp
3264	ELSE
3265	!
3266	!-- Adjust number of raindrops to avoid nonlinear effects in
3267	!-- sedimentation and evaporation of rain drops due to too small or
3268	!-- too big weights of rain drops (Stevens and Seifert, 2008).
3269	IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN
3270	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
3271	ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN
3272	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
3273	ENDIF
3274
3275	ENDIF
3276
3277	IF ( microphysics_morrison ) THEN
3278	IF ( qc(k,j,i) <= eps_sb ) THEN
3279	qc(k,j,i) = 0.0_wp
3280	nc(k,j,i) = 0.0_wp
3281	ELSE
3282	IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN
3283	nc(k,j,i) = qc(k,j,i) * hyrho(k) / xamin * flag
3284	ENDIF
3285	ENDIF
3286	ENDIF
3287
3288	ENDDO
3289
3290	END SUBROUTINE adjust_cloud_ij
3291
3292	!------------------------------------------------------------------------------!
3293	! Description:
3294	! ------------
3295	!> Calculate number of activated condensation nucleii after simple activation
3296	!> scheme of Twomey, 1959.
3297	!------------------------------------------------------------------------------!
3298	SUBROUTINE activation_ij( i, j )
3299
3300	IMPLICIT NONE
3301
3302	INTEGER(iwp) :: i !<
3303	INTEGER(iwp) :: j !<
3304	INTEGER(iwp) :: k !<
3305
3306	REAL(wp) :: activ !<
3307	REAL(wp) :: afactor !<
3308	REAL(wp) :: beta_act !<
3309	REAL(wp) :: bfactor !<
3310	REAL(wp) :: flag !< flag to indicate first grid level above surface
3311	REAL(wp) :: k_act !<
3312	REAL(wp) :: n_act !<
3313	REAL(wp) :: n_ccn !<
3314	REAL(wp) :: s_0 !<
3315	REAL(wp) :: sat_max !<
3316	REAL(wp) :: sigma !<
3317	REAL(wp) :: sigma_act !<
3318
3319	DO k = nzb+1, nzt
3320	!
3321	!-- Predetermine flag to mask topography
3322	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3323	!
3324	!-- Call calculation of supersaturation
3325	CALL supersaturation ( i, j, k )
3326	!
3327	!-- Prescribe parameters for activation
3328	!-- (see: Bott + Trautmann, 2002, Atm. Res., 64)
3329	k_act = 0.7_wp
3330	activ = 0.0_wp
3331
3332	IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
3333	!
3334	!-- Compute the number of activated Aerosols
3335	!-- (see: Twomey, 1959, Pure and applied Geophysics, 43)
3336	n_act = na_init * sat**k_act
3337	!
3338	!-- Compute the number of cloud droplets
3339	!-- (see: Morrison + Grabowski, 2007, JAS, 64)
3340	! activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro
3341
3342	!
3343	!-- Compute activation rate after Khairoutdinov and Kogan
3344	!-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
3345	sat_max = 0.8_wp / 100.0_wp
3346	activ = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN &
3347	( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) / &
3348	dt_micro
3349
3350	nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)
3351	ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk ) THEN
3352	!
3353	!-- Curvature effect (afactor) with surface tension
3354	!-- parameterization by Straka (2009)
3355	sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
3356	afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
3357	!
3358	!-- Solute effect (bfactor)
3359	bfactor = vanthoff * molecular_weight_of_water * &
3360	rho_s / ( molecular_weight_of_solute * rho_l )
3361
3362	!
3363	!-- Prescribe power index that describes the soluble fraction
3364	!-- of an aerosol particle (beta).
3365	!-- (see: Morrison + Grabowski, 2007, JAS, 64)
3366	beta_act = 0.5_wp
3367	sigma_act = sigma_bulk**( 1.0_wp + beta_act )
3368	!
3369	!-- Calculate mean geometric supersaturation (s_0) with
3370	!-- parameterization by Khvorostyanov and Curry (2006)
3371	s_0 = dry_aerosol_radius *(- ( 1.0_wp + beta_act ) ) &
3372	( 4.0_wp * afactor*3 / ( 27.0_wp bfactor ) )**0.5_wp
3373
3374	!
3375	!-- Calculate number of activated CCN as a function of
3376	!-- supersaturation and taking Koehler theory into account
3377	!-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
3378	n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( &
3379	LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
3380	activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
3381
3382	nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro * flag), na_init)
3383	ENDIF
3384
3385	ENDDO
3386
3387	END SUBROUTINE activation_ij
3388
3389	!------------------------------------------------------------------------------!
3390	! Description:
3391	! ------------
3392	!> Calculate condensation rate for cloud water content (after Khairoutdinov and
3393	!> Kogan, 2000).
3394	!------------------------------------------------------------------------------!
3395	SUBROUTINE condensation_ij( i, j )
3396
3397	IMPLICIT NONE
3398
3399	INTEGER(iwp) :: i !<
3400	INTEGER(iwp) :: j !<
3401	INTEGER(iwp) :: k !<
3402
3403	REAL(wp) :: cond !<
3404	REAL(wp) :: cond_max !<
3405	REAL(wp) :: dc !<
3406	REAL(wp) :: evap !<
3407	REAL(wp) :: flag !< flag to indicate first grid level above surface
3408	REAL(wp) :: g_fac !<
3409	REAL(wp) :: nc_0 !<
3410	REAL(wp) :: temp !<
3411	REAL(wp) :: xc !<
3412
3413
3414	DO k = nzb+1, nzt
3415	!
3416	!-- Predetermine flag to mask topography
3417	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3418	!
3419	!-- Call calculation of supersaturation
3420	CALL supersaturation ( i, j, k )
3421	!
3422	!-- Actual temperature:
3423	temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
3424
3425	g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
3426	l_v / ( thermal_conductivity_l * temp ) &
3427	+ r_v * temp / ( diff_coeff_l * e_s ) &
3428	)
3429	!
3430	!-- Mean weight of cloud drops
3431	IF ( nc(k,j,i) <= 0.0_wp) CYCLE
3432	xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
3433	!
3434	!-- Weight averaged diameter of cloud drops:
3435	dc = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
3436	!
3437	!-- Integral diameter of cloud drops
3438	nc_0 = nc(k,j,i) * dc
3439	!
3440	!-- Condensation needs only to be calculated in supersaturated regions
3441	IF ( sat > 0.0_wp ) THEN
3442	!
3443	!-- Condensation rate of cloud water content
3444	!-- after KK scheme.
3445	!-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
3446	cond = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
3447	cond_max = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
3448	cond = MIN( cond, cond_max / dt_micro )
3449
3450	qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag
3451	ELSEIF ( sat < 0.0_wp ) THEN
3452	evap = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
3453	evap = MAX( evap, -qc(k,j,i) / dt_micro )
3454
3455	qc(k,j,i) = qc(k,j,i) + evap * dt_micro * flag
3456	ENDIF
3457	ENDDO
3458
3459	END SUBROUTINE condensation_ij
3460
3461
3462	!------------------------------------------------------------------------------!
3463	! Description:
3464	! ------------
3465	!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
3466	!------------------------------------------------------------------------------!
3467	SUBROUTINE autoconversion_ij( i, j )
3468
3469	IMPLICIT NONE
3470
3471	INTEGER(iwp) :: i !<
3472	INTEGER(iwp) :: j !<
3473	INTEGER(iwp) :: k !<
3474
3475	REAL(wp) :: alpha_cc !<
3476	REAL(wp) :: autocon !<
3477	REAL(wp) :: dissipation !<
3478	REAL(wp) :: flag !< flag to indicate first grid level above surface
3479	REAL(wp) :: k_au !<
3480	REAL(wp) :: l_mix !<
3481	REAL(wp) :: nc_auto !<
3482	REAL(wp) :: nu_c !<
3483	REAL(wp) :: phi_au !<
3484	REAL(wp) :: r_cc !<
3485	REAL(wp) :: rc !<
3486	REAL(wp) :: re_lambda !<
3487	REAL(wp) :: sigma_cc !<
3488	REAL(wp) :: tau_cloud !<
3489	REAL(wp) :: xc !<
3490
3491	DO k = nzb+1, nzt
3492	!
3493	!-- Predetermine flag to mask topography
3494	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3495	IF ( microphysics_morrison ) THEN
3496	nc_auto = nc(k,j,i)
3497	ELSE
3498	nc_auto = nc_const
3499	ENDIF
3500
3501	IF ( qc(k,j,i) > eps_sb .AND. nc_auto > eps_mr ) THEN
3502
3503	k_au = k_cc / ( 20.0_wp * x0 )
3504	!
3505	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
3506	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
3507	tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) ), &
3508	0.0_wp )
3509	!
3510	!-- Universal function for autoconversion process
3511	!-- (Seifert and Beheng, 2006):
3512	phi_au = 600.0_wp * tau_cloud*0.68_wp ( 1.0_wp - tau_cloud0.68_wp )3
3513	!
3514	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
3515	!-- (Use constant nu_c = 1.0_wp instead?)
3516	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
3517	!
3518	!-- Mean weight of cloud droplets:
3519	xc = hyrho(k) * qc(k,j,i) / nc_auto
3520	!
3521	!-- Parameterized turbulence effects on autoconversion (Seifert,
3522	!-- Nuijens and Stevens, 2010)
3523	IF ( collision_turbulence ) THEN
3524	!
3525	!-- Weight averaged radius of cloud droplets:
3526	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
3527
3528	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
3529	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
3530	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
3531	!
3532	!-- Mixing length (neglecting distance to ground and stratification)
3533	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
3534	!
3535	!-- Limit dissipation rate according to Seifert, Nuijens and
3536	!-- Stevens (2010)
3537	dissipation = MIN( 0.06_wp, diss(k,j,i) )
3538	!
3539	!-- Compute Taylor-microscale Reynolds number:
3540	re_lambda = 6.0_wp / 11.0_wp * &
3541	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
3542	SQRT( 15.0_wp / kin_vis_air ) * &
3543	dissipation**( 1.0_wp / 6.0_wp )
3544	!
3545	!-- The factor of 1.0E4 is needed to convert the dissipation rate
3546	!-- from m2 s-3 to cm2 s-3.
3547	k_au = k_au * ( 1.0_wp + &
3548	dissipation * 1.0E4_wp * &
3549	( re_lambda * 1.0E-3_wp )*0.25_wp &
3550	( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) / &
3551	sigma_cc )**2 &
3552	) + beta_cc &
3553	) &
3554	)
3555	ENDIF
3556	!
3557	!-- Autoconversion rate (Seifert and Beheng, 2006):
3558	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
3559	( nu_c + 1.0_wp )*2 qc(k,j,i)*2 xc*2 &
3560	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
3561	rho_surface
3562	autocon = MIN( autocon, qc(k,j,i) / dt_micro )
3563
3564	qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag
3565	qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag
3566	nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro * flag
3567	IF ( microphysics_morrison ) THEN
3568	nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp * &
3569	autocon / x0 * hyrho(k) * dt_micro * flag )
3570	ENDIF
3571
3572	ENDIF
3573
3574	ENDDO
3575
3576	END SUBROUTINE autoconversion_ij
3577
3578	!------------------------------------------------------------------------------!
3579	! Description:
3580	! ------------
3581	!> Autoconversion process (Kessler, 1969).
3582	!------------------------------------------------------------------------------!
3583	SUBROUTINE autoconversion_kessler_ij( i, j )
3584
3585
3586	IMPLICIT NONE
3587
3588	INTEGER(iwp) :: i !<
3589	INTEGER(iwp) :: j !<
3590	INTEGER(iwp) :: k !<
3591	INTEGER(iwp) :: k_wall !< topography top index
3592
3593	REAL(wp) :: dqdt_precip !<
3594	REAL(wp) :: flag !< flag to indicate first grid level above surface
3595
3596	!
3597	!-- Determine vertical index of topography top
3598	k_wall = get_topography_top_index_ji( j, i, 's' )
3599	DO k = nzb+1, nzt
3600	!
3601	!-- Predetermine flag to mask topography
3602	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3603
3604	IF ( qc(k,j,i) > ql_crit ) THEN
3605	dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
3606	ELSE
3607	dqdt_precip = 0.0_wp
3608	ENDIF
3609
3610	qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
3611	q(k,j,i) = q(k,j,i) - dqdt_precip * dt_micro * flag
3612	pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp * d_exner(k) &
3613	* flag
3614
3615	!
3616	!-- Compute the rain rate (stored on surface grid point)
3617	prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
3618
3619	ENDDO
3620
3621	END SUBROUTINE autoconversion_kessler_ij
3622
3623	!------------------------------------------------------------------------------!
3624	! Description:
3625	! ------------
3626	!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
3627	!------------------------------------------------------------------------------!
3628	SUBROUTINE accretion_ij( i, j )
3629
3630	IMPLICIT NONE
3631
3632	INTEGER(iwp) :: i !<
3633	INTEGER(iwp) :: j !<
3634	INTEGER(iwp) :: k !<
3635
3636	REAL(wp) :: accr !<
3637	REAL(wp) :: flag !< flag to indicate first grid level above surface
3638	REAL(wp) :: k_cr !<
3639	REAL(wp) :: nc_accr !<
3640	REAL(wp) :: phi_ac !<
3641	REAL(wp) :: tau_cloud !<
3642	REAL(wp) :: xc !<
3643
3644
3645	DO k = nzb+1, nzt
3646	!
3647	!-- Predetermine flag to mask topography
3648	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3649	IF ( microphysics_morrison ) THEN
3650	nc_accr = nc(k,j,i)
3651	ELSE
3652	nc_accr = nc_const
3653	ENDIF
3654
3655	IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) .AND. &
3656	( nc_accr > eps_mr ) ) THEN
3657	!
3658	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
3659	tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
3660	!
3661	!-- Universal function for accretion process
3662	!-- (Seifert and Beheng, 2001):
3663	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
3664
3665	!
3666	!-- Mean weight of cloud drops
3667	xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
3668	!
3669	!-- Parameterized turbulence effects on autoconversion (Seifert,
3670	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
3671	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
3672	IF ( collision_turbulence ) THEN
3673	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
3674	MIN( 600.0_wp, &
3675	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
3676	)
3677	ELSE
3678	k_cr = k_cr0
3679	ENDIF
3680	!
3681	!-- Accretion rate (Seifert and Beheng, 2006):
3682	accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * &
3683	SQRT( rho_surface * hyrho(k) )
3684	accr = MIN( accr, qc(k,j,i) / dt_micro )
3685
3686	qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
3687	qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
3688	IF ( microphysics_morrison ) THEN
3689	nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc * &
3690	hyrho(k) * dt_micro * flag &
3691	)
3692	ENDIF
3693
3694
3695	ENDIF
3696
3697	ENDDO
3698
3699	END SUBROUTINE accretion_ij
3700
3701
3702	!------------------------------------------------------------------------------!
3703	! Description:
3704	! ------------
3705	!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
3706	!------------------------------------------------------------------------------!
3707	SUBROUTINE selfcollection_breakup_ij( i, j )
3708
3709	IMPLICIT NONE
3710
3711	INTEGER(iwp) :: i !<
3712	INTEGER(iwp) :: j !<
3713	INTEGER(iwp) :: k !<
3714
3715	REAL(wp) :: breakup !<
3716	REAL(wp) :: dr !<
3717	REAL(wp) :: flag !< flag to indicate first grid level above surface
3718	REAL(wp) :: phi_br !<
3719	REAL(wp) :: selfcoll !<
3720
3721	DO k = nzb+1, nzt
3722	!
3723	!-- Predetermine flag to mask topography
3724	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3725
3726	IF ( qr(k,j,i) > eps_sb ) THEN
3727	!
3728	!-- Selfcollection rate (Seifert and Beheng, 2001):
3729	selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * SQRT( hyrho(k) * rho_surface )
3730	!
3731	!-- Weight averaged diameter of rain drops:
3732	dr = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
3733	!
3734	!-- Collisional breakup rate (Seifert, 2008):
3735	IF ( dr >= 0.3E-3_wp ) THEN
3736	phi_br = k_br * ( dr - 1.1E-3_wp )
3737	breakup = selfcoll * ( phi_br + 1.0_wp )
3738	ELSE
3739	breakup = 0.0_wp
3740	ENDIF
3741
3742	selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
3743	nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
3744
3745	ENDIF
3746	ENDDO
3747
3748	END SUBROUTINE selfcollection_breakup_ij
3749
3750
3751	!------------------------------------------------------------------------------!
3752	! Description:
3753	! ------------
3754	!> Evaporation of precipitable water. Condensation is neglected for
3755	!> precipitable water. Call for grid point i,j
3756	!------------------------------------------------------------------------------!
3757	SUBROUTINE evaporation_rain_ij( i, j )
3758
3759	IMPLICIT NONE
3760
3761	INTEGER(iwp) :: i !<
3762	INTEGER(iwp) :: j !<
3763	INTEGER(iwp) :: k !<
3764
3765	REAL(wp) :: dr !<
3766	REAL(wp) :: evap !<
3767	REAL(wp) :: evap_nr !<
3768	REAL(wp) :: f_vent !<
3769	REAL(wp) :: flag !< flag to indicate first grid level above surface
3770	REAL(wp) :: g_evap !<
3771	REAL(wp) :: lambda_r !<
3772	REAL(wp) :: mu_r !<
3773	REAL(wp) :: mu_r_2 !<
3774	REAL(wp) :: mu_r_5d2 !<
3775	REAL(wp) :: nr_0 !<
3776	REAL(wp) :: temp !<
3777	REAL(wp) :: xr !<
3778
3779	DO k = nzb+1, nzt
3780	!
3781	!-- Predetermine flag to mask topography
3782	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3783
3784	IF ( qr(k,j,i) > eps_sb ) THEN
3785	!
3786	!-- Call calculation of supersaturation
3787	CALL supersaturation ( i, j, k )
3788	!
3789	!-- Evaporation needs only to be calculated in subsaturated regions
3790	IF ( sat < 0.0_wp ) THEN
3791	!
3792	!-- Actual temperature:
3793	temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
3794
3795	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v / &
3796	( thermal_conductivity_l * temp ) + &
3797	r_v * temp / ( diff_coeff_l * e_s ) &
3798	)
3799	!
3800	!-- Mean weight of rain drops
3801	xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
3802	!
3803	!-- Weight averaged diameter of rain drops:
3804	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
3805	!
3806	!-- Compute ventilation factor and intercept parameter
3807	!-- (Seifert and Beheng, 2006; Seifert, 2008):
3808	IF ( ventilation_effect ) THEN
3809	!
3810	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
3811	!-- Stevens and Seifert, 2008):
3812	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
3813	!
3814	!-- Slope parameter of gamma distribution (Seifert, 2008):
3815	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
3816	( mu_r + 1.0_wp ) &
3817	)**( 1.0_wp / 3.0_wp ) / dr
3818
3819	mu_r_2 = mu_r + 2.0_wp
3820	mu_r_5d2 = mu_r + 2.5_wp
3821
3822	f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) + &
3823	b_vent * schmidt_p_1d3 * &
3824	SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * &
3825	lambda_r*( -mu_r_5d2 ) &
3826	( 1.0_wp - &
3827	0.5_wp * ( b_term / a_term ) * &
3828	( lambda_r / ( c_term + lambda_r ) &
3829	)**mu_r_5d2 - &
3830	0.125_wp * ( b_term / a_term )*2 &
3831	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
3832	)**mu_r_5d2 - &
3833	0.0625_wp * ( b_term / a_term )*3 &
3834	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
3835	)**mu_r_5d2 - &
3836	0.0390625_wp * ( b_term / a_term )*4 &
3837	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
3838	)**mu_r_5d2 &
3839	)
3840
3841	nr_0 = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) / &
3842	gamm( mu_r + 1.0_wp )
3843	ELSE
3844	f_vent = 1.0_wp
3845	nr_0 = nr(k,j,i) * dr
3846	ENDIF
3847	!
3848	!-- Evaporation rate of rain water content (Seifert and Beheng, 2006):
3849	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
3850	evap = MAX( evap, -qr(k,j,i) / dt_micro )
3851	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
3852	-nr(k,j,i) / dt_micro )
3853
3854	qr(k,j,i) = qr(k,j,i) + evap * dt_micro * flag
3855	nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag
3856
3857	ENDIF
3858	ENDIF
3859
3860	ENDDO
3861
3862	END SUBROUTINE evaporation_rain_ij
3863
3864
3865	!------------------------------------------------------------------------------!
3866	! Description:
3867	! ------------
3868	!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
3869	!> Call for grid point i,j
3870	!------------------------------------------------------------------------------!
3871	SUBROUTINE sedimentation_cloud_ij( i, j )
3872
3873	IMPLICIT NONE
3874
3875	INTEGER(iwp) :: i !<
3876	INTEGER(iwp) :: j !<
3877	INTEGER(iwp) :: k !<
3878
3879	REAL(wp) :: flag !< flag to indicate first grid level above surface
3880	REAL(wp) :: nc_sedi !<
3881
3882	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nc !<
3883	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
3884
3885	sed_qc(nzt+1) = 0.0_wp
3886	sed_nc(nzt+1) = 0.0_wp
3887
3888
3889	DO k = nzt, nzb+1, -1
3890	!
3891	!-- Predetermine flag to mask topography
3892	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3893	IF ( microphysics_morrison ) THEN
3894	nc_sedi = nc(k,j,i)
3895	ELSE
3896	nc_sedi = nc_const
3897	ENDIF
3898	!
3899	!-- Sedimentation fluxes for number concentration are only calculated
3900	!-- for cloud_scheme = 'morrison'
3901	IF ( microphysics_morrison ) THEN
3902	IF ( qc(k,j,i) > eps_sb .AND. nc(k,j,i) > eps_mr ) THEN
3903	sed_nc(k) = sed_qc_const * &
3904	( qc(k,j,i) * hyrho(k) )*( 2.0_wp / 3.0_wp ) &
3905	( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
3906	ELSE
3907	sed_nc(k) = 0.0_wp
3908	ENDIF
3909
3910	sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) * &
3911	nc(k,j,i) / dt_micro + sed_nc(k+1) &
3912	) * flag
3913
3914	nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) * &
3915	ddzu(k+1) / hyrho(k) * dt_micro * flag
3916	ENDIF
3917
3918	IF ( qc(k,j,i) > eps_sb .AND. nc_sedi > eps_mr ) THEN
3919	sed_qc(k) = sed_qc_const * nc_sedi*( -2.0_wp / 3.0_wp ) &
3920	( qc(k,j,i) * hyrho(k) )*( 5.0_wp / 3.0_wp ) flag
3921	ELSE
3922	sed_qc(k) = 0.0_wp
3923	ENDIF
3924
3925	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / &
3926	dt_micro + sed_qc(k+1) &
3927	) * flag
3928
3929	q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
3930	hyrho(k) * dt_micro * flag
3931	qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
3932	hyrho(k) * dt_micro * flag
3933	pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
3934	hyrho(k) * lv_d_cp * d_exner(k) * dt_micro &
3935	* flag
3936
3937	!
3938	!-- Compute the precipitation rate of cloud (fog) droplets
3939	IF ( call_microphysics_at_all_substeps ) THEN
3940	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * &
3941	weight_substep(intermediate_timestep_count) * flag
3942	ELSE
3943	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
3944	ENDIF
3945
3946	ENDDO
3947
3948	END SUBROUTINE sedimentation_cloud_ij
3949
3950
3951	!------------------------------------------------------------------------------!
3952	! Description:
3953	! ------------
3954	!> Computation of sedimentation flux. Implementation according to Stevens
3955	!> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j
3956	!------------------------------------------------------------------------------!
3957	SUBROUTINE sedimentation_rain_ij( i, j )
3958
3959	IMPLICIT NONE
3960
3961	INTEGER(iwp) :: i !< running index x direction
3962	INTEGER(iwp) :: j !< running index y direction
3963	INTEGER(iwp) :: k !< running index z direction
3964	INTEGER(iwp) :: k_run !<
3965	INTEGER(iwp) :: m !< running index surface elements
3966	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
3967	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
3968
3969	REAL(wp) :: c_run !<
3970	REAL(wp) :: d_max !<
3971	REAL(wp) :: d_mean !<
3972	REAL(wp) :: d_min !<
3973	REAL(wp) :: dr !<
3974	REAL(wp) :: flux !<
3975	REAL(wp) :: flag !< flag to indicate first grid level above surface
3976	REAL(wp) :: lambda_r !<
3977	REAL(wp) :: mu_r !<
3978	REAL(wp) :: z_run !<
3979
3980	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<
3981	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<
3982	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<
3983	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<
3984	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<
3985	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<
3986	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<
3987	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<
3988
3989	!
3990	!-- Compute velocities
3991	DO k = nzb+1, nzt
3992	!
3993	!-- Predetermine flag to mask topography
3994	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3995
3996	IF ( qr(k,j,i) > eps_sb ) THEN
3997	!
3998	!-- Weight averaged diameter of rain drops:
3999	dr = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
4000	!
4001	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
4002	!-- Stevens and Seifert, 2008):
4003	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
4004	!
4005	!-- Slope parameter of gamma distribution (Seifert, 2008):
4006	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
4007	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
4008
4009	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
4010	a_term - b_term * ( 1.0_wp + &
4011	c_term / lambda_r )*( -1.0_wp &
4012	( mu_r + 1.0_wp ) ) &
4013	) &
4014	) * flag
4015	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
4016	a_term - b_term * ( 1.0_wp + &
4017	c_term / lambda_r )*( -1.0_wp &
4018	( mu_r + 4.0_wp ) ) &
4019	) &
4020	) * flag
4021	ELSE
4022	w_nr(k) = 0.0_wp
4023	w_qr(k) = 0.0_wp
4024	ENDIF
4025	ENDDO
4026	!
4027	!-- Adjust boundary values using surface data type.
4028	!-- Upward facing non-natural
4029	surf_s = bc_h(0)%start_index(j,i)
4030	surf_e = bc_h(0)%end_index(j,i)
4031	DO m = surf_s, surf_e
4032	k = bc_h(0)%k(m)
4033	w_nr(k-1) = w_nr(k)
4034	w_qr(k-1) = w_qr(k)
4035	ENDDO
4036	!
4037	!-- Downward facing non-natural
4038	surf_s = bc_h(1)%start_index(j,i)
4039	surf_e = bc_h(1)%end_index(j,i)
4040	DO m = surf_s, surf_e
4041	k = bc_h(1)%k(m)
4042	w_nr(k+1) = w_nr(k)
4043	w_qr(k+1) = w_qr(k)
4044	ENDDO
4045	!
4046	!-- Neumann boundary condition at model top
4047	w_nr(nzt+1) = 0.0_wp
4048	w_qr(nzt+1) = 0.0_wp
4049	!
4050	!-- Compute Courant number
4051	DO k = nzb+1, nzt
4052	!
4053	!-- Predetermine flag to mask topography
4054	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4055
4056	c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) * &
4057	dt_micro * ddzu(k) * flag
4058	c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * &
4059	dt_micro * ddzu(k) * flag
4060	ENDDO
4061	!
4062	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
4063	IF ( limiter_sedimentation ) THEN
4064
4065	DO k = nzb+1, nzt
4066	!
4067	!-- Predetermine flag to mask topography
4068	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4069
4070	d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
4071	d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
4072	d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
4073
4074	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
4075	2.0_wp * d_max, &
4076	ABS( d_mean ) ) * flag
4077
4078	d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
4079	d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
4080	d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
4081
4082	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
4083	2.0_wp * d_max, &
4084	ABS( d_mean ) ) * flag
4085	ENDDO
4086
4087	ELSE
4088
4089	nr_slope = 0.0_wp
4090	qr_slope = 0.0_wp
4091
4092	ENDIF
4093
4094	sed_nr(nzt+1) = 0.0_wp
4095	sed_qr(nzt+1) = 0.0_wp
4096	!
4097	!-- Compute sedimentation flux
4098	DO k = nzt, nzb+1, -1
4099	!
4100	!-- Predetermine flag to mask topography
4101	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4102	!
4103	!-- Sum up all rain drop number densities which contribute to the flux
4104	!-- through k-1/2
4105	flux = 0.0_wp
4106	z_run = 0.0_wp ! height above z(k)
4107	k_run = k
4108	c_run = MIN( 1.0_wp, c_nr(k) )
4109	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
4110	flux = flux + hyrho(k_run) * &
4111	( nr(k_run,j,i) + nr_slope(k_run) * ( 1.0_wp - c_run ) * &
4112	0.5_wp ) * c_run * dzu(k_run) * flag
4113	z_run = z_run + dzu(k_run) * flag
4114	k_run = k_run + 1 * flag
4115	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) * flag
4116	ENDDO
4117	!
4118	!-- It is not allowed to sediment more rain drop number density than
4119	!-- available
4120	flux = MIN( flux, &
4121	hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * dt_micro )
4122
4123	sed_nr(k) = flux / dt_micro * flag
4124	nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / &
4125	hyrho(k) * dt_micro * flag
4126	!
4127	!-- Sum up all rain water content which contributes to the flux
4128	!-- through k-1/2
4129	flux = 0.0_wp
4130	z_run = 0.0_wp ! height above z(k)
4131	k_run = k
4132	c_run = MIN( 1.0_wp, c_qr(k) )
4133
4134	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
4135
4136	flux = flux + hyrho(k_run) * &
4137	( qr(k_run,j,i) + qr_slope(k_run) * ( 1.0_wp - c_run ) * &
4138	0.5_wp ) * c_run * dzu(k_run) * flag
4139	z_run = z_run + dzu(k_run) * flag
4140	k_run = k_run + 1 * flag
4141	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) * flag
4142
4143	ENDDO
4144	!
4145	!-- It is not allowed to sediment more rain water content than available
4146	flux = MIN( flux, &
4147	hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * dt_micro )
4148
4149	sed_qr(k) = flux / dt_micro * flag
4150
4151	qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
4152	hyrho(k) * dt_micro * flag
4153	q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
4154	hyrho(k) * dt_micro * flag
4155	pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
4156	hyrho(k) * lv_d_cp * d_exner(k) * dt_micro &
4157	* flag
4158	!
4159	!-- Compute the rain rate
4160	IF ( call_microphysics_at_all_substeps ) THEN
4161	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) &
4162	* weight_substep(intermediate_timestep_count) * flag
4163	ELSE
4164	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
4165	ENDIF
4166
4167	ENDDO
4168
4169	END SUBROUTINE sedimentation_rain_ij
4170
4171
4172	!------------------------------------------------------------------------------!
4173	! Description:
4174	! ------------
4175	!> This subroutine computes the precipitation amount due to gravitational
4176	!> settling of rain and cloud (fog) droplets
4177	!------------------------------------------------------------------------------!
4178	SUBROUTINE calc_precipitation_amount_ij( i, j )
4179
4180	IMPLICIT NONE
4181
4182	INTEGER(iwp) :: i !< running index x direction
4183	INTEGER(iwp) :: j !< running index y direction
4184	INTEGER(iwp) :: k !< running index z direction
4185	INTEGER(iwp) :: m !< running index surface elements
4186	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
4187	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
4188
4189	IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
4190	( .NOT. call_microphysics_at_all_substeps .OR. &
4191	intermediate_timestep_count == intermediate_timestep_count_max ) ) &
4192	THEN
4193
4194	surf_s = bc_h(0)%start_index(j,i)
4195	surf_e = bc_h(0)%end_index(j,i)
4196	DO m = surf_s, surf_e
4197	k = bc_h(0)%k(m)
4198	precipitation_amount(j,i) = precipitation_amount(j,i) + &
4199	prr(k,j,i) * hyrho(k) * dt_3d
4200	ENDDO
4201
4202	ENDIF
4203
4204	END SUBROUTINE calc_precipitation_amount_ij
4205
4206
4207	!------------------------------------------------------------------------------!
4208	! Description:
4209	! ------------
4210	!> Computation of the diagnostic supersaturation sat, actual liquid water
4211	!< temperature t_l and saturation water vapor mixing ratio q_s
4212	!------------------------------------------------------------------------------!
4213	SUBROUTINE supersaturation ( i,j,k )
4214
4215	IMPLICIT NONE
4216
4217	INTEGER(iwp) :: i !< running index
4218	INTEGER(iwp) :: j !< running index
4219	INTEGER(iwp) :: k !< running index
4220
4221	REAL(wp) :: alpha !< correction factor
4222	!
4223	!-- Actual liquid water temperature:
4224	t_l = exner(k) * pt(k,j,i)
4225	!
4226	!-- Calculate water vapor saturation pressure
4227	e_s = magnus( t_l )
4228	!
4229	!-- Computation of saturation mixing ratio:
4230	q_s = rd_d_rv * e_s / ( hyp(k) - e_s )
4231	!
4232	!-- Correction factor
4233	alpha = rd_d_rv * lv_d_rd * lv_d_cp / ( t_l * t_l )
4234	!
4235	!-- Correction of the approximated value
4236	!-- (see: Cuijpers + Duynkerke, 1993, JAS, 23)
4237	q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / ( 1.0_wp + alpha * q_s )
4238	!
4239	!-- Supersaturation:
4240	!-- Not in case of microphysics_kessler or microphysics_sat_adjust
4241	!-- since qr is unallocated
4242	IF ( .NOT. microphysics_kessler .AND. &
4243	.NOT. microphysics_sat_adjust ) THEN
4244	sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
4245	ENDIF
4246
4247	END SUBROUTINE supersaturation
4248
4249
4250	!------------------------------------------------------------------------------!
4251	! Description:
4252	! ------------
4253	!> Calculation of the liquid water content (0%-or-100%-scheme). This scheme is
4254	!> used by the one and the two moment cloud physics scheme. Using the two moment
4255	!> scheme, this calculation results in the cloud water content.
4256	!------------------------------------------------------------------------------!
4257	SUBROUTINE calc_liquid_water_content
4258
4259
4260
4261	IMPLICIT NONE
4262
4263	INTEGER(iwp) :: i !<
4264	INTEGER(iwp) :: j !<
4265	INTEGER(iwp) :: k !<
4266
4267	REAL(wp) :: flag !< flag to indicate first grid level above surface
4268
4269
4270	DO i = nxlg, nxrg
4271	DO j = nysg, nyng
4272	DO k = nzb+1, nzt
4273	!
4274	!-- Predetermine flag to mask topography
4275	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4276
4277	!
4278	!-- Call calculation of supersaturation located
4279	CALL supersaturation( i, j, k )
4280
4281	!
4282	!-- Compute the liquid water content
4283	IF ( microphysics_seifert .AND. .NOT. microphysics_morrison ) &
4284	THEN
4285	IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp ) THEN
4286	qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) ) * flag
4287	ql(k,j,i) = ( qc(k,j,i) + qr(k,j,i) ) * flag
4288	ELSE
4289	IF ( q(k,j,i) < qr(k,j,i) ) q(k,j,i) = qr(k,j,i)
4290	qc(k,j,i) = 0.0_wp
4291	ql(k,j,i) = qr(k,j,i) * flag
4292	ENDIF
4293	ELSEIF ( microphysics_morrison ) THEN
4294	ql(k,j,i) = qc(k,j,i) + qr(k,j,i) * flag
4295	ELSE
4296	IF ( ( q(k,j,i) - q_s ) > 0.0_wp ) THEN
4297	qc(k,j,i) = ( q(k,j,i) - q_s ) * flag
4298	ql(k,j,i) = qc(k,j,i) * flag
4299	ELSE
4300	qc(k,j,i) = 0.0_wp
4301	ql(k,j,i) = 0.0_wp
4302	ENDIF
4303	ENDIF
4304	ENDDO
4305	ENDDO
4306	ENDDO
4307
4308	END SUBROUTINE calc_liquid_water_content
4309
4310	!------------------------------------------------------------------------------!
4311	! Description:
4312	! ------------
4313	!> This function computes the gamma function (Press et al., 1992).
4314	!> The gamma function is needed for the calculation of the evaporation
4315	!> of rain drops.
4316	!------------------------------------------------------------------------------!
4317	FUNCTION gamm( xx )
4318
4319	IMPLICIT NONE
4320
4321	INTEGER(iwp) :: j !<
4322
4323