Home

Context Navigation

bulk_cloud_model_mod.f90 @ 4289

Last change on this file since 4289 was 4289, checked in by knoop, 4 years ago
Removed parameters precipitation and precipitation_amount_interval from bulk cloud model namelist
Property svn:keywords set to `Id`
File size: 202.0 KB

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