source: palm/trunk/SOURCE/bulk_cloud_model_mod.f90 @ 4109

Last change on this file since 4109 was 4109, checked in by suehring, 2 years ago

Control discretization of advection term: separate initialization of WS advection flags for momentum and scalars. In this context, resort the bits and do some minor formatting; Make initialization of scalar-advection flags more flexible, i.e. introduce an arguemnt list to indicate non-cyclic boundaries (required for decycled scalars such as chemical species or aerosols); Introduce extended 'degradation zones', where horizontal advection of passive scalars is discretized by first-order scheme at all grid points that in the vicinity of buildings (<= 3 grid points). Even though no building is within the numerical stencil, first-order scheme is used. At fourth and fifth grid point the order of the horizontal advection scheme is successively upgraded. These extended degradation zones are used to avoid stationary numerical oscillations, which are responsible for high concentration maxima that may appear under shear-free stable conditions. Therefore, an additional 3D interger array used to store control flags is introduced; Change interface for scalar advection routine; Bugfix, avoid uninitialized value sk_num in vector version of WS scalar advection; Chemistry: Decycling boundary conditions are only set at the ghost points not on the prognostic grid points; Land-surface model: Relax checks for non-consistent initialization in case static or dynamic input is provided. For example, soil_temperature or deep_soil_temperature is not mandatory any more if dynamic input is available. Also, improper settings of x_type in namelist are only checked if no static file is available.

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