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

Last change on this file since 4110 was 4110, checked in by suehring, 5 years ago

last changes documented

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