Home

Context Navigation

bulk_cloud_model_mod.f90 @ 4521

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