Home

Context Navigation

bulk_cloud_model_mod.f90 @ 4418

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