Home

Context Navigation

bulk_cloud_model_mod.f90 @ 4027

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