Home

Context Navigation

bulk_cloud_model_mod.f90 @ 4457

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