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

Last change on this file since 4598 was 4577, checked in by raasch, 10 months ago

further re-formatting to follow the PALM coding standard

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