Home

Context Navigation

bulk_cloud_model_mod.f90 @ 4110

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