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

Last change on this file since 4517 was 4517, checked in by raasch, 14 months ago

added restart with MPI-IO for reading local arrays

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