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

Last change on this file since 4506 was 4506, checked in by schwenkel, 14 months ago

Use magnus formula for liquid water temperature

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