Home

Context Navigation

bulk_cloud_model_mod.f90 @ 3885

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