source: palm/trunk/SOURCE/bulk_cloud_model_mod.f90

Last change on this file was 4872, checked in by raasch, 7 months ago

internal switches removed from namelists

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