Home

Context Navigation

bulk_cloud_model_mod.f90 @ 4370

Last change on this file since 4370 was 4370, checked in by raasch, 5 years ago
bugfixes for previous commit: unused variables removed, Temperton-fft usage on GPU, openacc porting of vector version of Obukhov length calculation, collective read switched off on NEC to avoid hanging; some vector directives added in prognostic equations to force vectorization on Intel19 compiler, configuration files for NEC Aurora added
Property svn:keywords set to `Id`
File size: 203.0 KB

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