Changeset 4542 for palm


Ignore:
Timestamp:
May 19, 2020 3:45:12 PM (7 months ago)
Author:
raasch
Message:

files re-formatted to follow the PALM coding standard, redundant if statement removed

Location:
palm/trunk/SOURCE
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/bulk_cloud_model_mod.f90

    r4535 r4542  
    11!> @file bulk_cloud_model_mod.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
    9 !
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    13 !
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     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/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
     
    2524! -----------------
    2625! $Id$
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4535 2020-05-15 12:07:23Z raasch
    2729! bugfix for restart data format query
    28 ! 
     30!
    2931! 4533 2020-05-14 14:46:46Z schwenkel
    3032! Reformat intrinsic function
    31 ! 
     33!
    3234! 4532 2020-05-14 13:41:35Z schwenkel
    3335! Calculate mean droplet radius assuming gamma distibution for condensation
    34 ! 
     36!
    3537! 4521 2020-05-06 11:39:49Z schwenkel
    3638! Rename variable
    37 ! 
     39!
    3840! 4517 2020-05-03 14:29:30Z raasch
    3941! added restart with MPI-IO for reading local arrays
    40 ! 
     42!
    4143! 4506 2020-04-21 10:57:45Z schwenkel
    4244! Use correct magnus formula for liquid water temperature
    43 ! 
     45!
    4446! 4502 2020-04-17 16:14:16Z schwenkel
    4547! Implementation of ice microphysics
    46 ! 
     48!
    4749! 4495 2020-04-13 20:11:20Z raasch
    4850! restart data handling with MPI-IO added
    49 ! 
     51!
    5052! 4457 2020-03-11 14:20:43Z raasch
    5153! use statement for exchange horiz added
    52 ! 
     54!
    5355! 4418 2020-02-21 09:41:13Z raasch
    5456! bugfix for raindrop number adjustment
    55 ! 
     57!
    5658! 4370 2020-01-10 14:00:44Z raasch
    5759! vector directives added to force vectorization on Intel19 compiler
    58 ! 
     60!
    5961! 4360 2020-01-07 11:25:50Z suehring
    6062! Introduction of wall_flags_total_0, which currently sets bits based on static
    6163! topography information used in wall_flags_static_0
    62 ! 
     64!
    6365! 4329 2019-12-10 15:46:36Z motisi
    6466! Renamed wall_flags_0 to wall_flags_static_0
    65 ! 
     67!
    6668! 4289 2019-11-05 14:33:41Z knoop
    6769! Removed parameters precipitation and precipitation_amount_interval from namelist
    68 ! 
     70!
    6971! 4268 2019-10-17 11:29:38Z schwenkel
    7072! Introducing bcm_boundary_conditions
    71 ! 
     73!
    7274! 4182 2019-08-22 15:20:23Z scharf
    7375! Corrected "Former revisions" section
    74 ! 
     76!
    7577! 4168 2019-08-16 13:50:17Z suehring
    7678! Replace function get_topography_top_index by topo_top_ind
    77 ! 
     79!
    7880! 4110 2019-07-22 17:05:21Z suehring
    79 ! Pass integer flag array as well as boundary flags to WS scalar advection 
     81! Pass integer flag array as well as boundary flags to WS scalar advection
    8082! routine
    81 ! 
     83!
    8284! 4109 2019-07-22 17:00:34Z suehring
    8385! Added microphyics scheme 'morrision_no_rain'
    84 ! 
     86!
    8587! 3931 2019-04-24 16:34:28Z schwenkel
    8688! Added bcm_exchange_horiz which is called after non_transport_physics
    87 ! 
     89!
    8890! 3885 2019-04-11 11:29:34Z kanani
    89 ! Changes related to global restructuring of location messages and introduction 
     91! Changes related to global restructuring of location messages and introduction
    9092! of additional debug messages
    91 ! 
     93!
    9294! 3874 2019-04-08 16:53:48Z knoop
    9395! Implemented non_transport_physics module interfaces
    94 ! 
     96!
    9597! 3870 2019-04-08 13:44:34Z knoop
    9698! Moving prognostic equations of bcm into bulk_cloud_model_mod
    97 ! 
     99!
    98100! 3869 2019-04-08 11:54:20Z knoop
    99101! moving the furniture around ;-)
    100 ! 
     102!
    101103! 3786 2019-03-06 16:58:03Z raasch
    102104! unsed variables removed
    103 ! 
     105!
    104106! 3767 2019-02-27 08:18:02Z raasch
    105107! unused variable for file index removed from rrd-subroutines parameter list
    106 ! 
     108!
    107109! 3724 2019-02-06 16:28:23Z kanani
    108110! Correct double-used log_point_s unit
    109 ! 
     111!
    110112! 3700 2019-01-26 17:03:42Z knoop
    111113! nopointer option removed
     
    116118! ------------
    117119!> Calculate bulk cloud microphysics.
    118 !------------------------------------------------------------------------------!
     120!--------------------------------------------------------------------------------------------------!
    119121 MODULE bulk_cloud_model_mod
    120122
    121123
    122     USE advec_s_bc_mod,                                                        &
     124    USE advec_s_bc_mod,                                                                            &
    123125        ONLY:  advec_s_bc
    124126
    125     USE advec_s_pw_mod,                                                        &
     127    USE advec_s_pw_mod,                                                                            &
    126128        ONLY:  advec_s_pw
    127129
    128     USE advec_s_up_mod,                                                        &
     130    USE advec_s_up_mod,                                                                            &
    129131        ONLY:  advec_s_up
    130132
    131     USE advec_ws,                                                              &
     133    USE advec_ws,                                                                                  &
    132134        ONLY:  advec_s_ws
    133135
    134     USE arrays_3d,                                                             &
    135         ONLY:  ddzu, diss, dzu, dzw, hyp, hyrho,                               &
    136                nc, nc_1, nc_2, nc_3, nc_p, nr, nr_1, nr_2, nr_3, nr_p,         &
    137                precipitation_amount, prr, pt, d_exner, pt_init, q, ql, ql_1,   &
    138                qc, qc_1, qc_2, qc_3, qc_p, qr, qr_1, qr_2, qr_3, qr_p,         &
    139                exner, zu, tnc_m, tnr_m, tqc_m, tqr_m, tend, rdf_sc,            &
    140                flux_l_qc, flux_l_qr, flux_l_nc, flux_l_nr,                     &
    141                flux_s_qc, flux_s_qr, flux_s_nc, flux_s_nr,                     &
    142                diss_l_qc, diss_l_qr, diss_l_nc, diss_l_nr,                     &
    143                diss_s_qc, diss_s_qr, diss_s_nc, diss_s_nr,                     &
    144                ni, ni_1, ni_2, ni_3, ni_p, tni_m,                              &
    145                qi, qi_1, qi_2, qi_3, qi_p, tqi_m,                              &
    146                flux_l_qi, flux_l_ni, flux_s_qi, flux_s_ni,                     &
    147                diss_l_qi, diss_l_ni, diss_s_qi, diss_s_ni
    148 
    149 
    150     USE averaging,                                                             &
    151         ONLY:  nc_av, nr_av, prr_av, qc_av, ql_av, qr_av, ni_av, qi_av
    152 
    153     USE basic_constants_and_equations_mod,                                     &
    154         ONLY:  c_p, g, lv_d_cp, lv_d_rd, l_v, magnus, magnus_ice, magnus_tl,   &
    155                molecular_weight_of_solute,                                     &
    156                molecular_weight_of_water, pi, rho_l, rho_s, r_d, r_v, vanthoff,&
    157                exner_function, exner_function_invers, ideal_gas_law_rho,       &
    158                ideal_gas_law_rho_pt, barometric_formula, rd_d_rv, l_s,         &
    159                ls_d_cp
    160 
    161     USE control_parameters,                                                    &
    162         ONLY:  bc_dirichlet_l,                                                 &
    163                bc_dirichlet_n,                                                 &
    164                bc_dirichlet_r,                                                 &
    165                bc_dirichlet_s,                                                 &
    166                bc_radiation_l,                                                 &
    167                bc_radiation_n,                                                 &
    168                bc_radiation_r,                                                 &
    169                bc_radiation_s,                                                 &
    170                debug_output,                                                   &
    171                dt_3d, dt_do2d_xy, intermediate_timestep_count,                 &
    172                intermediate_timestep_count_max, large_scale_forcing,           &
    173                lsf_surf, pt_surface, restart_data_format_output, rho_surface,  &
    174                surface_pressure,                                               &
    175                time_do2d_xy, message_string, initializing_actions,             &
    176                ws_scheme_sca, scalar_advec, timestep_scheme, tsc,              &
    177                loop_optimization, simulated_time
    178 
    179     USE cpulog,                                                                &
     136    USE arrays_3d,                                                                                 &
     137        ONLY:  diss_s_qc, diss_s_qr, diss_s_nc, diss_s_nr,                                         &
     138               diss_l_qc, diss_l_qr, diss_l_nc, diss_l_nr,                                         &
     139               diss_l_qi, diss_l_ni, diss_s_qi, diss_s_ni,                                         &
     140               ddzu, diss, dzu, dzw, hyp, hyrho,                                                   &
     141               exner, zu, tnc_m, tnr_m, tqc_m, tqr_m, tend, rdf_sc,                                &
     142               flux_l_qc, flux_l_qr, flux_l_nc, flux_l_nr,                                         &
     143               flux_l_qi, flux_l_ni, flux_s_qi, flux_s_ni,                                         &
     144               flux_s_qc, flux_s_qr, flux_s_nc, flux_s_nr,                                         &
     145               nc, nc_1, nc_2, nc_3, nc_p, nr, nr_1, nr_2, nr_3, nr_p,                             &
     146               ni, ni_1, ni_2, ni_3, ni_p, tni_m,                                                  &
     147               precipitation_amount, prr, pt, d_exner, pt_init, q, ql, ql_1,                       &
     148               qc, qc_1, qc_2, qc_3, qc_p, qr, qr_1, qr_2, qr_3, qr_p,                             &
     149               qi, qi_1, qi_2, qi_3, qi_p, tqi_m
     150
     151
     152    USE averaging,                                                                                 &
     153        ONLY:  nc_av,ni_av, nr_av, prr_av, qc_av, ql_av, qi_av, qr_av
     154
     155    USE basic_constants_and_equations_mod,                                                         &
     156        ONLY:  c_p, g, lv_d_cp, lv_d_rd, l_v, magnus, magnus_ice, magnus_tl,                       &
     157               exner_function, exner_function_invers, ideal_gas_law_rho,                           &
     158               ideal_gas_law_rho_pt, barometric_formula, rd_d_rv, l_s,                             &
     159               ls_d_cp,                                                                            &
     160               molecular_weight_of_solute,                                                         &
     161               molecular_weight_of_water, pi, rho_l, rho_s, r_d, r_v, vanthoff
     162
     163    USE control_parameters,                                                                        &
     164        ONLY:  bc_dirichlet_l,                                                                     &
     165               bc_dirichlet_n,                                                                     &
     166               bc_dirichlet_r,                                                                     &
     167               bc_dirichlet_s,                                                                     &
     168               bc_radiation_l,                                                                     &
     169               bc_radiation_n,                                                                     &
     170               bc_radiation_r,                                                                     &
     171               bc_radiation_s,                                                                     &
     172               debug_output,                                                                       &
     173               dt_3d, dt_do2d_xy, intermediate_timestep_count,                                     &
     174               intermediate_timestep_count_max, large_scale_forcing,                               &
     175               loop_optimization, simulated_time,                                                  &
     176               lsf_surf, pt_surface, restart_data_format_output, rho_surface,                      &
     177               surface_pressure,                                                                   &
     178               time_do2d_xy, message_string, initializing_actions,                                 &
     179               ws_scheme_sca, scalar_advec, timestep_scheme, tsc
     180
     181    USE cpulog,                                                                                    &
    180182        ONLY:  cpu_log, log_point, log_point_s
    181183
    182     USE diffusion_s_mod,                                                       &
     184    USE diffusion_s_mod,                                                                           &
    183185        ONLY:  diffusion_s
    184186
    185     USE grid_variables,                                                        &
     187    USE grid_variables,                                                                            &
    186188        ONLY:  dx, dy
    187189
    188     USE indices,                                                               &
    189         ONLY:  advc_flags_s,                                                   &
    190                nbgp, nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb, nzt,     &
    191                topo_top_ind,                                                   &
     190    USE indices,                                                                                   &
     191        ONLY:  advc_flags_s,                                                                       &
     192               nbgp, nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb, nzt,                         &
     193               topo_top_ind,                                                                       &
    192194               wall_flags_total_0
    193195
    194196    USE kinds
    195197
    196     USE pegrid,                                                                &
     198    USE pegrid,                                                                                    &
    197199        ONLY:  threads_per_task
    198200
    199     USE restart_data_mpi_io_mod,                                               &
     201    USE restart_data_mpi_io_mod,                                                                   &
    200202        ONLY:  rd_mpi_io_check_array, rrd_mpi_io, wrd_mpi_io
    201203
    202     USE statistics,                                                            &
    203         ONLY:  weight_pres, weight_substep, sums_wsncs_ws_l, sums_wsnrs_ws_l,  &
    204                sums_wsqcs_ws_l, sums_wsqrs_ws_l,                               &
    205                sums_wsqis_ws_l, sums_wsnis_ws_l
    206 
    207     USE surface_mod,                                                           &
    208         ONLY :  bc_h,                                                          &
    209                 surf_bulk_cloud_model,                                         &
    210                 surf_microphysics_morrison, surf_microphysics_seifert,         &
    211                 surf_microphysics_ice_phase,                               &
    212                 surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
     204    USE statistics,                                                                                &
     205        ONLY:  sums_wsncs_ws_l, sums_wsnis_ws_l, sums_wsnrs_ws_l, sums_wsqcs_ws_l, sums_wsqis_ws_l,&
     206               sums_wsqrs_ws_l, weight_pres, weight_substep
     207
     208    USE surface_mod,                                                                               &
     209        ONLY :  bc_h,                                                                              &
     210                surf_bulk_cloud_model,                                                             &
     211                surf_microphysics_morrison, surf_microphysics_seifert,                             &
     212                surf_microphysics_ice_phase,                                                       &
     213                surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,                        &
    213214                surf_usm_v
    214215
     
    218219    CHARACTER (LEN=20)   ::  cloud_scheme = 'saturation_adjust'           !< namelist parameter
    219220
    220     LOGICAL ::  aerosol_nacl =.TRUE.                             !< nacl aerosol for bulk scheme
    221     LOGICAL ::  aerosol_c3h4o4 =.FALSE.                          !< malonic acid aerosol for bulk scheme
    222     LOGICAL ::  aerosol_nh4no3 =.FALSE.                          !< malonic acid aerosol for bulk scheme
    223 
    224     LOGICAL ::  bulk_cloud_model = .FALSE.                       !< namelist parameter
    225 
    226     LOGICAL ::  cloud_water_sedimentation = .FALSE.       !< cloud water sedimentation
    227     LOGICAL ::  curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory
    228     LOGICAL ::  ice_crystal_sedimentation = .FALSE.        !< flag for ice crystal sedimentation
    229     LOGICAL ::  limiter_sedimentation = .TRUE.            !< sedimentation limiter
    230     LOGICAL ::  collision_turbulence = .FALSE.            !< turbulence effects
    231     LOGICAL ::  ventilation_effect = .TRUE.               !< ventilation effect
    232 
    233     LOGICAL ::  call_microphysics_at_all_substeps = .FALSE.      !< namelist parameter
    234     LOGICAL ::  microphysics_ice_phase = .FALSE.              !< use ice microphysics scheme
    235     LOGICAL ::  microphysics_sat_adjust = .FALSE.                !< use saturation adjust bulk scheme?
    236     LOGICAL ::  microphysics_kessler = .FALSE.                   !< use kessler bulk scheme?
    237     LOGICAL ::  microphysics_morrison = .FALSE.                  !< use 2-moment Morrison (add. prog. eq. for nc and qc)
    238     LOGICAL ::  microphysics_seifert = .FALSE.                   !< use 2-moment Seifert and Beheng scheme
    239     LOGICAL ::  microphysics_morrison_no_rain = .FALSE.          !< use 2-moment Morrison
    240     LOGICAL ::  precipitation = .FALSE.                          !< namelist parameter
    241 
    242     REAL(wp) ::  precipitation_amount_interval = 9999999.9_wp    !< namelist parameter
    243 
    244     REAL(wp) ::  a_1 = 8.69E-4_wp          !< coef. in turb. parametrization (cm-2 s3)
    245     REAL(wp) ::  a_2 = -7.38E-5_wp         !< coef. in turb. parametrization (cm-2 s3)
    246     REAL(wp) ::  a_3 = -1.40E-2_wp         !< coef. in turb. parametrization
    247     REAL(wp) ::  a_term = 9.65_wp          !< coef. for terminal velocity (m s-1)
    248     REAL(wp) ::  a_vent = 0.78_wp          !< coef. for ventilation effect
    249     REAL(wp) ::  b_1 = 11.45E-6_wp         !< coef. in turb. parametrization (m)
    250     REAL(wp) ::  b_2 = 9.68E-6_wp          !< coef. in turb. parametrization (m)
    251     REAL(wp) ::  b_3 = 0.62_wp             !< coef. in turb. parametrization
    252     REAL(wp) ::  b_term = 9.8_wp           !< coef. for terminal velocity (m s-1)
    253     REAL(wp) ::  b_vent = 0.308_wp         !< coef. for ventilation effect
    254     REAL(wp) ::  beta_cc = 3.09E-4_wp      !< coef. in turb. parametrization (cm-2 s3)
    255     REAL(wp) ::  c_1 = 4.82E-6_wp          !< coef. in turb. parametrization (m)
    256     REAL(wp) ::  c_2 = 4.8E-6_wp           !< coef. in turb. parametrization (m)
    257     REAL(wp) ::  c_3 = 0.76_wp             !< coef. in turb. parametrization
    258     REAL(wp) ::  c_const = 0.93_wp         !< const. in Taylor-microscale Reynolds number
    259     REAL(wp) ::  c_evap = 0.7_wp           !< constant in evaporation
    260     REAL(wp) ::  c_term = 600.0_wp         !< coef. for terminal velocity (m-1)
    261     REAL(wp) ::  diff_coeff_l = 0.23E-4_wp  !< diffusivity of water vapor (m2 s-1)
    262     REAL(wp) ::  eps_sb = 1.0E-10_wp       !< threshold in two-moments scheme
    263     REAL(wp) ::  eps_mr = 0.0_wp           !< threshold for morrison scheme
    264     REAL(wp) ::  k_cc = 9.44E09_wp         !< const. cloud-cloud kernel (m3 kg-2 s-1)
    265     REAL(wp) ::  k_cr0 = 4.33_wp           !< const. cloud-rain kernel (m3 kg-1 s-1)
    266     REAL(wp) ::  k_rr = 7.12_wp            !< const. rain-rain kernel (m3 kg-1 s-1)
    267     REAL(wp) ::  k_br = 1000.0_wp          !< const. in breakup parametrization (m-1)
    268     REAL(wp) ::  k_st = 1.2E8_wp           !< const. in drizzle parametrization (m-1 s-1)
    269     REAL(wp) ::  kin_vis_air = 1.4086E-5_wp  !< kin. viscosity of air (m2 s-1)
    270     REAL(wp) ::  prec_time_const = 0.001_wp  !< coef. in Kessler scheme (s-1)
    271     REAL(wp) ::  ql_crit = 0.0005_wp       !< coef. in Kessler scheme (kg kg-1)
    272     REAL(wp) ::  schmidt_p_1d3=0.8921121_wp  !< Schmidt number**0.33333, 0.71**0.33333
    273     REAL(wp) ::  sigma_gc = 1.3_wp         !< geometric standard deviation cloud droplets
    274     REAL(wp) ::  thermal_conductivity_l = 2.43E-2_wp  !< therm. cond. air (J m-1 s-1 K-1)
    275     REAL(wp) ::  w_precipitation = 9.65_wp  !< maximum terminal velocity (m s-1)
    276     REAL(wp) ::  x0 = 2.6E-10_wp           !< separating drop mass (kg)
    277     REAL(wp) ::  ximin = 4.42E-14_wp       !< minimum mass of ice crystal (kg) (D~10µm)
    278     REAL(wp) ::  xcmin = 4.18E-15_wp       !< minimum cloud drop size (kg) (~ 1µm)
    279     REAL(wp) ::  xrmin = 2.6E-10_wp        !< minimum rain drop size (kg)
    280     REAL(wp) ::  xrmax = 5.0E-6_wp         !< maximum rain drop site (kg)
    281 
     221    LOGICAL ::  aerosol_nacl =.TRUE.                         !< nacl aerosol for bulk scheme
     222    LOGICAL ::  aerosol_c3h4o4 =.FALSE.                      !< malonic acid aerosol for bulk scheme
     223    LOGICAL ::  aerosol_nh4no3 =.FALSE.                      !< malonic acid aerosol for bulk scheme
     224    LOGICAL ::  bulk_cloud_model = .FALSE.                   !< namelist parameter
     225    LOGICAL ::  call_microphysics_at_all_substeps = .FALSE.  !< namelist parameter
     226    LOGICAL ::  cloud_water_sedimentation = .FALSE.          !< cloud water sedimentation
     227    LOGICAL ::  collision_turbulence = .FALSE.               !< turbulence effects
     228    LOGICAL ::  curvature_solution_effects_bulk = .FALSE.    !< flag for considering koehler theory
     229    LOGICAL ::  ice_crystal_sedimentation = .FALSE.          !< flag for ice crystal sedimentation
     230    LOGICAL ::  limiter_sedimentation = .TRUE.               !< sedimentation limiter
     231    LOGICAL ::  microphysics_ice_phase = .FALSE.             !< use ice microphysics scheme
     232    LOGICAL ::  microphysics_kessler = .FALSE.               !< use kessler bulk scheme?
     233    LOGICAL ::  microphysics_morrison = .FALSE.              !< use 2-moment Morrison
     234                                                             !< (add. prog. eq. for nc and qc)
     235    LOGICAL ::  microphysics_morrison_no_rain = .FALSE.      !< use 2-moment Morrison
     236    LOGICAL ::  microphysics_sat_adjust = .FALSE.            !< use saturation adjust bulk scheme?
     237    LOGICAL ::  microphysics_seifert = .FALSE.               !< use 2-moment Seifert and Beheng
     238                                                             !< scheme
     239    LOGICAL ::  precipitation = .FALSE.                      !< namelist parameter
     240    LOGICAL ::  ventilation_effect = .TRUE.                  !< ventilation effect
     241
     242
     243    REAL(wp) ::  a_1 = 8.69E-4_wp                !< coef. in turb. parametrization    (cm-2 s3)
     244    REAL(wp) ::  a_2 = -7.38E-5_wp               !< coef. in turb. parametrization    (cm-2 s3)
     245    REAL(wp) ::  a_3 = -1.40E-2_wp               !< coef. in turb. parametrization
     246    REAL(wp) ::  a_term = 9.65_wp                !< coef. for terminal velocity       (m s-1)
     247    REAL(wp) ::  a_vent = 0.78_wp                !< coef. for ventilation effect
     248    REAL(wp) ::  b_1 = 11.45E-6_wp               !< coef. in turb. parametrization    (m)
     249    REAL(wp) ::  b_2 = 9.68E-6_wp                !< coef. in turb. parametrization    (m)
     250    REAL(wp) ::  b_3 = 0.62_wp                   !< coef. in turb. parametrization
     251    REAL(wp) ::  b_term = 9.8_wp                 !< coef. for terminal velocity       (m s-1)
     252    REAL(wp) ::  b_vent = 0.308_wp               !< coef. for ventilation effect
     253    REAL(wp) ::  beta_cc = 3.09E-4_wp            !< coef. in turb. parametrization    (cm-2 s3)
     254    REAL(wp) ::  c_1 = 4.82E-6_wp                !< coef. in turb. parametrization    (m)
     255    REAL(wp) ::  c_2 = 4.8E-6_wp                 !< coef. in turb. parametrization    (m)
     256    REAL(wp) ::  c_3 = 0.76_wp                   !< coef. in turb. parametrization
     257    REAL(wp) ::  c_const = 0.93_wp               !< const. in Taylor-microscale Reynolds number
     258    REAL(wp) ::  c_evap = 0.7_wp                 !< constant in evaporation
    282259    REAL(wp) ::  c_sedimentation = 2.0_wp        !< Courant number of sedimentation process
     260    REAL(wp) ::  c_term = 600.0_wp               !< coef. for terminal velocity       (m-1)
     261    REAL(wp) ::  diff_coeff_l = 0.23E-4_wp       !< diffusivity of water vapor        (m2 s-1)
    283262    REAL(wp) ::  dpirho_l                        !< 6.0 / ( pi * rho_l )
    284263    REAL(wp) ::  dry_aerosol_radius = 0.05E-6_wp !< dry aerosol radius
    285264    REAL(wp) ::  dt_micro                        !< microphysics time step
     265    REAL(wp) ::  eps_sb = 1.0E-10_wp             !< threshold in two-moments scheme
     266    REAL(wp) ::  dt_precipitation = 100.0_wp     !< timestep precipitation (s)
     267    REAL(wp) ::  e_s                             !< saturation water vapor pressure
     268    REAL(wp) ::  e_si                            !< saturation water vapor pressure over ice
     269    REAL(wp) ::  eps_mr = 0.0_wp                 !< threshold for morrison scheme
    286270    REAL(wp) ::  in_init = 1000.0_wp             !< initial number of potential ice nucleii
    287     REAL(wp) ::  sigma_bulk = 2.0_wp             !< width of aerosol spectrum
     271    REAL(wp) ::  k_cc = 9.44E09_wp               !< const. cloud-cloud kernel         (m3 kg-2 s-1)
     272    REAL(wp) ::  k_cr0 = 4.33_wp                 !< const. cloud-rain kernel          (m3 kg-1 s-1)
     273    REAL(wp) ::  k_rr = 7.12_wp                  !< const. rain-rain kernel           (m3 kg-1 s-1)
     274    REAL(wp) ::  k_br = 1000.0_wp                !< const. in breakup parametrization (m-1)
     275    REAL(wp) ::  k_st = 1.2E8_wp                 !< const. in drizzle parametrization (m-1 s-1)
     276    REAL(wp) ::  kin_vis_air = 1.4086E-5_wp      !< kin. viscosity of air             (m2 s-1)
    288277    REAL(wp) ::  na_init = 100.0E6_wp            !< Total particle/aerosol concentration (cm-3)
    289278    REAL(wp) ::  nc_const = 70.0E6_wp            !< cloud droplet concentration
    290     REAL(wp) ::  dt_precipitation = 100.0_wp     !< timestep precipitation (s)
     279    REAL(wp) ::  pirho_l                         !< pi * rho_l / 6.0
     280    REAL(wp) ::  prec_time_const = 0.001_wp      !< coef. in Kessler scheme           (s-1)
     281    REAL(wp) ::  precipitation_amount_interval = 9999999.9_wp     !< namelist parameter
     282    REAL(wp) ::  ql_crit = 0.0005_wp             !< coef. in Kessler scheme           (kg kg-1)
     283    REAL(wp) ::  q_s                             !< saturation mixing ratio
     284    REAL(wp) ::  q_si                            !< saturation mixing ratio over ice
     285    REAL(wp) ::  sat                             !< supersaturation
     286    REAL(wp) ::  sat_ice                         !< supersaturation over ice
     287    REAL(wp) ::  start_ice_microphysics = 0.0_wp !< time from which on ice microhysics are executed
    291288    REAL(wp) ::  sed_qc_const                    !< const. for sedimentation of cloud water
    292     REAL(wp) ::  pirho_l                         !< pi * rho_l / 6.0
    293     REAL(wp) ::  start_ice_microphysics = 0.0_wp !< time from which on ice microhysics are executed
    294 
    295     REAL(wp) ::  e_s     !< saturation water vapor pressure
    296     REAL(wp) ::  e_si    !< saturation water vapor pressure over ice
    297     REAL(wp) ::  q_s     !< saturation mixing ratio
    298     REAL(wp) ::  q_si    !< saturation mixing ratio over ice
    299     REAL(wp) ::  sat     !< supersaturation
    300     REAL(wp) ::  sat_ice !< supersaturation over ice
    301     REAL(wp) ::  t_l     !< liquid-(ice) water temperature
     289    REAL(wp) ::  schmidt_p_1d3=0.8921121_wp      !< Schmidt number**0.33333, 0.71**0.33333
     290    REAL(wp) ::  sigma_bulk = 2.0_wp             !< width of aerosol spectrum
     291    REAL(wp) ::  sigma_gc = 1.3_wp               !< geometric standard deviation cloud droplets
     292    REAL(wp) ::  t_l                             !< liquid-(ice) water temperature
     293    REAL(wp) ::  thermal_conductivity_l = 2.43E-2_wp  !< therm. cond. air         (J m-1 s-1 K-1)
     294    REAL(wp) ::  w_precipitation = 9.65_wp       !< maximum terminal velocity         (m s-1)
     295    REAL(wp) ::  x0 = 2.6E-10_wp                 !< separating drop mass              (kg)
     296    REAL(wp) ::  ximin = 4.42E-14_wp             !< minimum mass of ice crystal       (kg) (D~10µm)
     297    REAL(wp) ::  xcmin = 4.18E-15_wp             !< minimum cloud drop size           (kg) (~ 1µm)
     298    REAL(wp) ::  xrmin = 2.6E-10_wp              !< minimum rain drop size            (kg)
     299    REAL(wp) ::  xrmax = 5.0E-6_wp               !< maximum rain drop site            (kg)
     300
    302301
    303302    SAVE
     
    305304    PRIVATE
    306305
    307     PUBLIC bcm_parin, &
    308            bcm_check_parameters, &
    309            bcm_check_data_output, &
    310            bcm_check_data_output_pr, &
    311            bcm_init_arrays, &
    312            bcm_init, &
    313            bcm_header, &
    314            bcm_actions, &
    315            bcm_non_advective_processes, &
    316            bcm_exchange_horiz, &
    317            bcm_prognostic_equations, &
    318            bcm_boundary_conditions, &
    319            bcm_3d_data_averaging, &
    320            bcm_data_output_2d, &
    321            bcm_data_output_3d, &
    322            bcm_swap_timelevel, &
    323            bcm_rrd_global, &
    324            bcm_rrd_local, &
    325            bcm_wrd_global, &
    326            bcm_wrd_local, &
     306    PUBLIC bcm_parin,                                                                              &
     307           bcm_check_parameters,                                                                   &
     308           bcm_check_data_output,                                                                  &
     309           bcm_check_data_output_pr,                                                               &
     310           bcm_init_arrays,                                                                        &
     311           bcm_init,                                                                               &
     312           bcm_header,                                                                             &
     313           bcm_actions,                                                                            &
     314           bcm_non_advective_processes,                                                            &
     315           bcm_exchange_horiz,                                                                     &
     316           bcm_prognostic_equations,                                                               &
     317           bcm_boundary_conditions,                                                                &
     318           bcm_3d_data_averaging,                                                                  &
     319           bcm_data_output_2d,                                                                     &
     320           bcm_data_output_3d,                                                                     &
     321           bcm_swap_timelevel,                                                                     &
     322           bcm_rrd_global,                                                                         &
     323           bcm_rrd_local,                                                                          &
     324           bcm_wrd_global,                                                                         &
     325           bcm_wrd_local,                                                                          &
    327326           calc_liquid_water_content
    328327
    329     PUBLIC call_microphysics_at_all_substeps, &
    330            cloud_water_sedimentation, &
    331            bulk_cloud_model, &
    332            cloud_scheme, &
    333            collision_turbulence, &
    334            dt_precipitation, &
    335            microphysics_morrison, &
    336            microphysics_morrison_no_rain, &
    337            microphysics_sat_adjust, &
    338            microphysics_seifert, &
    339            microphysics_ice_phase, &
    340            na_init, &
    341            nc_const, &
    342            precipitation, &
    343            sigma_gc, &
    344            start_ice_microphysics, &
    345            ice_crystal_sedimentation, &
    346            in_init
     328    PUBLIC bulk_cloud_model,                                                                      &
     329           call_microphysics_at_all_substeps,                                                      &
     330           cloud_water_sedimentation,                                                              &
     331           cloud_scheme,                                                                           &
     332           collision_turbulence,                                                                   &
     333           dt_precipitation,                                                                       &
     334           ice_crystal_sedimentation,                                                              &
     335           in_init,                                                                                &
     336           microphysics_morrison,                                                                  &
     337           microphysics_morrison_no_rain,                                                          &
     338           microphysics_sat_adjust,                                                                &
     339           microphysics_seifert,                                                                  &
     340           microphysics_ice_phase,                                                                &
     341           na_init,                                                                                &
     342           nc_const,                                                                              &
     343           precipitation,                                                                          &
     344           sigma_gc,                                                                              &
     345           start_ice_microphysics
    347346
    348347    INTERFACE bcm_parin
     
    438437
    439438
    440 !------------------------------------------------------------------------------!
     439!--------------------------------------------------------------------------------------------------!
    441440! Description:
    442441! ------------
    443442!> Parin for &bulk_cloud_parameters for the bulk cloud module
    444 !------------------------------------------------------------------------------!
     443!--------------------------------------------------------------------------------------------------!
    445444    SUBROUTINE bcm_parin
    446445
     
    448447       IMPLICIT NONE
    449448
    450        CHARACTER (LEN=80)  ::  line  !< dummy string that contains the current line of the parameter file
    451 
    452        NAMELIST /bulk_cloud_parameters/  &
    453           aerosol_bulk, &
    454           c_sedimentation, &
    455           call_microphysics_at_all_substeps, &
    456           bulk_cloud_model, &
    457           cloud_scheme, &
    458           cloud_water_sedimentation, &
    459           collision_turbulence, &
    460           curvature_solution_effects_bulk, &
    461           dry_aerosol_radius, &
    462           limiter_sedimentation, &
    463           na_init, &
    464           nc_const, &
    465           sigma_bulk, &
    466           ventilation_effect, &
    467           ice_crystal_sedimentation, &
    468           microphysics_ice_phase, &
    469           start_ice_microphysics, &
    470           in_init
     449       CHARACTER (LEN=80)  ::  line  !< dummy string that contains the current line of the parameter
     450                                     !< file
     451
     452       NAMELIST /bulk_cloud_parameters/                                                            &
     453          aerosol_bulk,                                                                            &
     454          bulk_cloud_model,                                                                        &
     455          c_sedimentation,                                                                         &
     456          call_microphysics_at_all_substeps,                                                       &
     457          cloud_scheme,                                                                            &
     458          cloud_water_sedimentation,                                                               &
     459          collision_turbulence,                                                                    &
     460          curvature_solution_effects_bulk,                                                         &
     461          dry_aerosol_radius,                                                                      &
     462          ice_crystal_sedimentation,                                                               &
     463          in_init,                                                                                 &
     464          limiter_sedimentation,                                                                   &
     465          microphysics_ice_phase,                                                                  &
     466          na_init,                                                                                 &
     467          nc_const,                                                                                &
     468          sigma_bulk,                                                                              &
     469          start_ice_microphysics,                                                                  &
     470          ventilation_effect
    471471
    472472       line = ' '
     
    492492
    493493
    494 !------------------------------------------------------------------------------!
     494!--------------------------------------------------------------------------------------------------!
    495495! Description:
    496496! ------------
    497497!> Check parameters routine for bulk cloud module
    498 !------------------------------------------------------------------------------!
     498!--------------------------------------------------------------------------------------------------!
    499499    SUBROUTINE bcm_check_parameters
    500500
     
    503503!
    504504!--    Check cloud scheme
    505 !--    This scheme considers only saturation adjustment,
    506 !--    i.e. water vapor surplus is converted into liquid
    507 !--    water. No other microphysical processes are considered
     505!--    This scheme considers only saturation adjustment, i.e. water vapor surplus is converted into
     506!--    liquid water. No other microphysical processes are considered.
    508507       IF ( cloud_scheme == 'saturation_adjust' )  THEN
    509508          microphysics_sat_adjust = .TRUE.
     
    513512          microphysics_morrison_no_rain = .FALSE.
    514513!
    515 !--    This scheme includes all process of the seifert
    516 !--    beheng scheme (2001,2006). Especially rain processes are
    517 !--    considered with prognostic quantities of qr and nr.
    518 !--    Cloud droplet concentration is assumed to be constant and
    519 !--    qc is diagnostic.
    520 !--    Technical remark: The switch 'microphysics_seifert' allocates
    521 !--    fields of qr and nr and enables all rain processes.
     514!--    This scheme includes all process of the seifert beheng scheme (2001,2006). Especially rain
     515!--    processes are considered with prognostic quantities of qr and nr.
     516!--    Cloud droplet concentration is assumed to be constant and qc is diagnostic.
     517!--    Technical remark: The switch 'microphysics_seifert' allocates fields of qr and nr and enables
     518!--    all rain processes.
    522519       ELSEIF ( cloud_scheme == 'seifert_beheng' )  THEN
    523520          microphysics_sat_adjust = .FALSE.
     
    528525          microphysics_morrison_no_rain = .FALSE.
    529526!
    530 !--    The kessler scheme is a simplified scheme without any
    531 !--    prognostic quantities for microphyical variables but
    532 !--    considering autoconversion.
     527!--    The kessler scheme is a simplified scheme without any prognostic quantities for microphyical
     528!--    variables but considering autoconversion.
    533529       ELSEIF ( cloud_scheme == 'kessler' )  THEN
    534530          microphysics_sat_adjust = .FALSE.
     
    539535          microphysics_morrison_no_rain = .FALSE.
    540536!
    541 !--    The morrison scheme is an extension of the seifer beheng scheme
    542 !--    including also relevant processes for cloud droplet size particles
    543 !--    such as activation and an diagnostic mehtod for diffusional growth.
    544 !--    I.e. here all processes of Seifert and Beheng as well as of the
    545 !--    morrision scheme are used. Therefore, ztis includes prognostic
    546 !--    quantities for qc and nc.
    547 !--    Technical remark: The switch 'microphysics_morrison' allocates
    548 !--    fields of qc and nc and enables diagnostic diffusional growth and
    549 !--    activation.
     537!--    The morrison scheme is an extension of the seifer beheng scheme including also relevant
     538!--    processes for cloud droplet size particles such as activation and an diagnostic mehtod for
     539!--    diffusional growth.
     540!--    I.e. here all processes of Seifert and Beheng as well as of the morrision scheme are used.
     541!--    Therefore, ztis includes prognostic quantities for qc and nc.
     542!--    Technical remark: The switch 'microphysics_morrison' allocates fields of qc and nc and
     543!--    enables diagnostic diffusional growth and activation.
    550544       ELSEIF ( cloud_scheme == 'morrison' )  THEN
    551545          microphysics_sat_adjust = .FALSE.
     
    554548          microphysics_morrison   = .TRUE.
    555549          precipitation           = .TRUE.
    556           microphysics_morrison_no_rain = .FALSE.   
    557 !         
    558 !--    The 'morrision_no_rain' scheme includes only processes of morrision scheme
    559 !--    without the rain processes of seifert beheng. Therfore, the prog. quantities
    560 !--    of qr and nr remain unallocated. This might be appropiate for cloud in which
    561 !--    the size distribution is narrow, e.g. fog.
     550          microphysics_morrison_no_rain = .FALSE.
     551!
     552!--    The 'morrision_no_rain' scheme includes only processes of morrision scheme without the rain
     553!--    processes of seifert beheng. Therfore, the prog. quantities of qr and nr remain unallocated.
     554!--    This might be appropiate for cloud in which the size distribution is narrow, e.g. fog.
    562555       ELSEIF ( cloud_scheme == 'morrison_no_rain' )  THEN
    563556          microphysics_sat_adjust = .FALSE.
     
    568561          precipitation           = .FALSE.
    569562       ELSE
    570           message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
     563          message_string = 'unknown cloud microphysics scheme cloud_scheme ="' //                  &
    571564                           TRIM( cloud_scheme ) // '"'
    572565          CALL message( 'check_parameters', 'PA0357', 1, 2, 0, 6, 0 )
     
    582575          ELSE
    583576             IF ( precipitation_amount_interval > dt_do2d_xy )  THEN
    584                 WRITE( message_string, * )  'precipitation_amount_interval = ',   &
    585                     precipitation_amount_interval, ' must not be larger than ',   &
    586                     'dt_do2d_xy = ', dt_do2d_xy
     577                WRITE( message_string, * )  'precipitation_amount_interval = ',                    &
     578                   precipitation_amount_interval, ' must not be larger than ',                     &
     579                   'dt_do2d_xy = ', dt_do2d_xy
    587580                CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 )
    588581             ENDIF
     
    590583       ENDIF
    591584
    592        ! TODO: find better sollution for circular dependency problem
     585       ! TODO: find better solution for circular dependency problem
    593586       surf_bulk_cloud_model = bulk_cloud_model
    594587       surf_microphysics_morrison = microphysics_morrison
     
    617610    END SUBROUTINE bcm_check_parameters
    618611
    619 !------------------------------------------------------------------------------!
     612!--------------------------------------------------------------------------------------------------!
    620613! Description:
    621614! ------------
    622615!> Check data output for bulk cloud module
    623 !------------------------------------------------------------------------------!
     616!--------------------------------------------------------------------------------------------------!
    624617    SUBROUTINE bcm_check_data_output( var, unit )
    625618
     
    633626          CASE ( 'nc' )
    634627             IF ( .NOT.  microphysics_morrison )  THEN
    635                 message_string = 'output of "' // TRIM( var ) // '" ' //       &
    636                                  'requires ' //                                &
     628                message_string = 'output of "' // TRIM( var ) // '" ' // 'requires ' //            &
    637629                                 'cloud_scheme = "morrison"'
    638630                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
     
    642634          CASE ( 'ni' )
    643635             IF ( .NOT.  microphysics_ice_phase )  THEN
    644                 message_string = 'output of "' // TRIM( var ) // '" ' //       &
    645                                  'requires ' //                                &
     636                message_string = 'output of "' // TRIM( var ) // '" ' // 'requires ' //            &
    646637                                 'microphysics_ice_phase = ".TRUE."'
    647638                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
     
    651642          CASE ( 'nr' )
    652643             IF ( .NOT.  microphysics_seifert )  THEN
    653                 message_string = 'output of "' // TRIM( var ) // '" ' //       &
    654                                  'requires ' //                                &
     644                message_string = 'output of "' // TRIM( var ) // '" ' // 'requires ' //            &
    655645                                 'cloud_scheme = "seifert_beheng"'
    656646                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
     
    660650          CASE ( 'prr' )
    661651             IF ( microphysics_sat_adjust )  THEN
    662                 message_string = 'output of "' // TRIM( var ) // '" ' //       &
    663                                  'is not available for ' //                    &
     652                message_string = 'output of "' // TRIM( var ) // '" ' // 'is not available for ' //&
    664653                                 'cloud_scheme = "saturation_adjust"'
    665654                CALL message( 'check_parameters', 'PA0423', 1, 2, 0, 6, 0 )
     
    672661          CASE ( 'qi' )
    673662             IF ( .NOT.  microphysics_ice_phase ) THEN
    674                 message_string = 'output of "' // TRIM( var ) // '" ' //       &
    675                                  'requires ' //                                &
     663                message_string = 'output of "' // TRIM( var ) // '" ' // 'requires ' //            &
    676664                                 'microphysics_ice_phase = ".TRUE."'
    677665                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
     
    681669          CASE ( 'qr' )
    682670             IF ( .NOT.  microphysics_seifert ) THEN
    683                 message_string = 'output of "' // TRIM( var ) // '" ' //       &
    684                                  'requires ' //                                &
     671                message_string = 'output of "' // TRIM( var ) // '" ' // 'requires ' //            &
    685672                                 'cloud_scheme = "seifert_beheng"'
    686673                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
     
    689676
    690677          CASE ( 'pra*' )
    691              IF ( .NOT. microphysics_kessler  .AND.                            &
    692                   .NOT. microphysics_seifert )  THEN
    693                 message_string = 'output of "' // TRIM( var ) // '" ' //       &
    694                                  'requires ' //                                &
     678             IF ( .NOT. microphysics_kessler  .AND.  .NOT. microphysics_seifert )  THEN
     679                message_string = 'output of "' // TRIM( var ) // '" ' // 'requires ' //            &
    695680                                 'cloud_scheme = "kessler" or "seifert_beheng"'
    696681                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
    697682             ENDIF
    698683! TODO: find sollution (maybe connected to flow_statistics redesign?)
    699 !              IF ( j == 1 )  THEN
    700 !                 message_string = 'temporal averaging of precipitation ' //     &
    701 !                           'amount "' // TRIM( var ) // '" is not possible'
    702 !                 CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
    703 !              ENDIF
     684!             IF ( j == 1 )  THEN
     685!                message_string = 'temporal averaging of precipitation ' //                        &
     686!                                 'amount "' // TRIM( var ) // '" is not possible'
     687!                CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
     688!             ENDIF
    704689             unit = 'mm'
    705690
    706691          CASE ( 'prr*' )
    707              IF ( .NOT. microphysics_kessler  .AND.                            &
    708                   .NOT. microphysics_seifert )  THEN
    709                 message_string = 'output of "' // TRIM( var ) // '"' //        &
    710                          ' requires' //                                        &
    711                          ' cloud_scheme = "kessler" or "seifert_beheng"'
     692             IF ( .NOT. microphysics_kessler  .AND.  .NOT. microphysics_seifert )  THEN
     693                message_string = 'output of "' // TRIM( var ) // '"' // ' requires' //             &
     694                                 ' cloud_scheme = "kessler" or "seifert_beheng"'
    712695                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
    713696             ENDIF
     
    723706
    724707
    725 !------------------------------------------------------------------------------!
     708!--------------------------------------------------------------------------------------------------!
    726709! Description:
    727710! ------------
    728711!> Check data output of profiles for bulk cloud module
    729 !------------------------------------------------------------------------------!
     712!--------------------------------------------------------------------------------------------------!
    730713    SUBROUTINE bcm_check_data_output_pr( variable, var_count, unit, dopr_unit )
    731714
    732        USE arrays_3d,                                                          &
     715       USE arrays_3d,                                                                              &
    733716           ONLY: zu
    734717
    735        USE control_parameters,                                                 &
     718       USE control_parameters,                                                                     &
    736719           ONLY: data_output_pr
    737720
    738        USE profil_parameter,                                                   &
     721       USE profil_parameter,                                                                       &
    739722           ONLY: dopr_index
    740723
    741        USE statistics,                                                         &
     724       USE statistics,                                                                             &
    742725           ONLY: hom, statistic_regions
    743726
     
    757740          CASE ( 'nc' )
    758741             IF ( .NOT.  microphysics_morrison )  THEN
    759                 message_string = 'data_output_pr = ' //                        &
    760                                  TRIM( data_output_pr(var_count) ) //          &
    761                                  ' is not implemented for' //                  &
    762                                  ' cloud_scheme /= morrison'
     742                message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) //       &
     743                                 ' is not implemented for' // ' cloud_scheme /= morrison'
    763744                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
    764745             ENDIF
     
    771752          CASE ( 'ni' )
    772753             IF ( .NOT.  microphysics_ice_phase )  THEN
    773                 message_string = 'data_output_pr = ' //                        &
    774                                  TRIM( data_output_pr(var_count) ) //          &
    775                                  ' is not implemented for' //                  &
    776                                  ' microphysics_ice_phase = ".F."'
     754                message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) //       &
     755                                 ' is not implemented for' // ' microphysics_ice_phase = ".F."'
    777756                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
    778757             ENDIF
     
    785764          CASE ( 'nr' )
    786765             IF ( .NOT.  microphysics_seifert )  THEN
    787                 message_string = 'data_output_pr = ' //                        &
    788                                  TRIM( data_output_pr(var_count) ) //          &
    789                                  ' is not implemented for' //                  &
    790                                  ' cloud_scheme /= seifert_beheng'
     766                message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) //       &
     767                                 ' is not implemented for' // ' cloud_scheme /= seifert_beheng'
    791768                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
    792769             ENDIF
     
    799776          CASE ( 'prr' )
    800777             IF ( microphysics_sat_adjust )  THEN
    801                 message_string = 'data_output_pr = ' //                        &
    802                                  TRIM( data_output_pr(var_count) ) //          &
    803                                  ' is not available for' //                    &
    804                                  ' cloud_scheme = saturation_adjust'
     778                message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) //       &
     779                                 ' is not available for' // ' cloud_scheme = saturation_adjust'
    805780                CALL message( 'check_parameters', 'PA0422', 1, 2, 0, 6, 0 )
    806              ENDIF 
     781             ENDIF
    807782             pr_index = 76
    808783             dopr_index(var_count) = pr_index
     
    820795          CASE ( 'qi' )
    821796             IF ( .NOT.  microphysics_ice_phase )  THEN
    822                 message_string = 'data_output_pr = ' //                        &
    823                                  TRIM( data_output_pr(var_count) ) //          &
    824                                  ' is not implemented for' //                  &
    825                                  ' microphysics_ice_phase = ".F."'
     797                message_string = 'data_output_pr = ' //  TRIM( data_output_pr(var_count) ) //      &
     798                                 ' is not implemented for' // ' microphysics_ice_phase = ".F."'
    826799                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
    827800             ENDIF
     
    834807          CASE ( 'qr' )
    835808             IF ( .NOT.  microphysics_seifert )  THEN
    836                 message_string = 'data_output_pr = ' //                        &
    837                                  TRIM( data_output_pr(var_count) ) //          &
    838                                  ' is not implemented for' //                  &
    839                                  ' cloud_scheme /= seifert_beheng'
     809                message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) //       &
     810                                 ' is not implemented for' // ' cloud_scheme /= seifert_beheng'
    840811                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
    841812             ENDIF
     
    854825
    855826
    856 !------------------------------------------------------------------------------!
     827!--------------------------------------------------------------------------------------------------!
    857828! Description:
    858829! ------------
    859830!> Allocate bulk cloud module arrays and define pointers
    860 !------------------------------------------------------------------------------!
     831!--------------------------------------------------------------------------------------------------!
    861832    SUBROUTINE bcm_init_arrays
    862833
    863        USE indices,                                                            &
     834       USE indices,                                                                                &
    864835           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
    865836
     
    887858!
    888859!--       3D-cloud drop water content, cloud drop concentration arrays
    889           ALLOCATE( nc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
    890                     nc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
    891                     nc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
    892                     qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
    893                     qc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
     860          ALLOCATE( nc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                           &
     861                    nc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                           &
     862                    nc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                           &
     863                    qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                           &
     864                    qc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                           &
    894865                    qc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    895866       ENDIF
     
    898869!
    899870!--       3D-rain water content, rain drop concentration arrays
    900           ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
    901                     nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
    902                     nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
    903                     qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
    904                     qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
     871          ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                           &
     872                    nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                           &
     873                    nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                           &
     874                    qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                           &
     875                    qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                           &
    905876                    qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    906877       ENDIF
     
    909880!
    910881!--       3D-cloud drop water content, cloud drop concentration arrays
    911           ALLOCATE( ni_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
    912                     ni_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
    913                     ni_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
    914                     qi_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
    915                     qi_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
     882          ALLOCATE( ni_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                           &
     883                    ni_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                           &
     884                    ni_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                           &
     885                    qi_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                           &
     886                    qi_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                           &
    916887                    qi_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    917888       ENDIF
     
    940911!
    941912!--    Arrays needed for reasons of speed optimization for cache version.
    942 !--    For the vector version the buffer arrays are not necessary,
    943 !--    because the the fluxes can swapped directly inside the loops of the
    944 !--    advection routines.
     913!--    For the vector version the buffer arrays are not necessary, because the the fluxes can
     914!--    swapped directly inside the loops of the advection routines.
    945915       IF ( loop_optimization /= 'vector' )  THEN
    946916          IF ( ws_scheme_sca )  THEN
    947917             IF ( microphysics_morrison )  THEN
    948                 ALLOCATE( flux_s_qc(nzb+1:nzt,0:threads_per_task-1),           &
    949                           diss_s_qc(nzb+1:nzt,0:threads_per_task-1),           &
    950                           flux_s_nc(nzb+1:nzt,0:threads_per_task-1),           &
     918                ALLOCATE( flux_s_qc(nzb+1:nzt,0:threads_per_task-1),                               &
     919                          diss_s_qc(nzb+1:nzt,0:threads_per_task-1),                               &
     920                          flux_s_nc(nzb+1:nzt,0:threads_per_task-1),                               &
    951921                          diss_s_nc(nzb+1:nzt,0:threads_per_task-1) )
    952                 ALLOCATE( flux_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
    953                           diss_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
    954                           flux_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
     922                ALLOCATE( flux_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                       &
     923                          diss_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                       &
     924                          flux_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                       &
    955925                          diss_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    956926             ENDIF
    957927             IF ( microphysics_seifert )  THEN
    958                 ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1),           &
    959                           diss_s_qr(nzb+1:nzt,0:threads_per_task-1),           &
    960                           flux_s_nr(nzb+1:nzt,0:threads_per_task-1),           &
     928                ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1),                               &
     929                          diss_s_qr(nzb+1:nzt,0:threads_per_task-1),                               &
     930                          flux_s_nr(nzb+1:nzt,0:threads_per_task-1),                               &
    961931                          diss_s_nr(nzb+1:nzt,0:threads_per_task-1) )
    962                 ALLOCATE( flux_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
    963                           diss_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
    964                           flux_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
     932                ALLOCATE( flux_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                       &
     933                          diss_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                       &
     934                          flux_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                       &
    965935                          diss_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    966936             ENDIF
    967937             IF ( microphysics_ice_phase )  THEN
    968                 ALLOCATE( flux_s_qi(nzb+1:nzt,0:threads_per_task-1),           &
    969                           diss_s_qi(nzb+1:nzt,0:threads_per_task-1),           &
    970                           flux_s_ni(nzb+1:nzt,0:threads_per_task-1),           &
     938                ALLOCATE( flux_s_qi(nzb+1:nzt,0:threads_per_task-1),                               &
     939                          diss_s_qi(nzb+1:nzt,0:threads_per_task-1),                               &
     940                          flux_s_ni(nzb+1:nzt,0:threads_per_task-1),                               &
    971941                          diss_s_ni(nzb+1:nzt,0:threads_per_task-1) )
    972                 ALLOCATE( flux_l_qi(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
    973                           diss_l_qi(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
    974                           flux_l_ni(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
     942                ALLOCATE( flux_l_qi(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                       &
     943                          diss_l_qi(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                       &
     944                          flux_l_ni(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                       &
    975945                          diss_l_ni(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    976946             ENDIF
     
    1002972
    1003973
    1004 !------------------------------------------------------------------------------!
     974!--------------------------------------------------------------------------------------------------!
    1005975! Description:
    1006976! ------------
    1007977!> Initialization of the bulk cloud module
    1008 !------------------------------------------------------------------------------!
     978!--------------------------------------------------------------------------------------------------!
    1009979    SUBROUTINE bcm_init
    1010980
     
    10481018             ENDIF
    10491019!
    1050 !--          Liquid water content and precipitation amount
    1051 !--          are zero at beginning of the simulation
     1020!--          Liquid water content and precipitation amount are zero at beginning of the simulation.
    10521021             ql = 0.0_wp
    10531022             qc = 0.0_wp
     
    10761045          ENDIF ! Only if not read_restart_data
    10771046!
    1078 !--       constant for the sedimentation of cloud water (2-moment cloud physics)
    1079           sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l )                &
    1080                              )**( 2.0_wp / 3.0_wp ) *                             &
     1047!--       Constant for the sedimentation of cloud water (2-moment cloud physics)
     1048          sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l ) )**( 2.0_wp / 3.0_wp ) *        &
    10811049                         EXP( 5.0_wp * LOG( sigma_gc )**2 )
    10821050
     
    10841052!--       Calculate timestep according to precipitation
    10851053          IF ( microphysics_seifert )  THEN
    1086              dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) /      &
    1087                                 w_precipitation
     1054             dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) / w_precipitation
    10881055          ENDIF
    10891056
     
    11221089
    11231090
    1124 !------------------------------------------------------------------------------!
     1091!--------------------------------------------------------------------------------------------------!
    11251092! Description:
    11261093! ------------
    11271094!> Header output for bulk cloud module
    1128 !------------------------------------------------------------------------------!
     1095!--------------------------------------------------------------------------------------------------!
    11291096    SUBROUTINE bcm_header ( io )
    11301097
     
    11661133
    11671134
    1168 1   FORMAT ( //' Bulk cloud module information:'/ &
    1169                ' ------------------------------------------'/ )
    1170 2   FORMAT ( '--> Bulk scheme with liquid water potential temperature and'/ &
     11351   FORMAT ( //' Bulk cloud module information:'/ ' ------------------------------------------'/ )
     11362   FORMAT ( '--> Bulk scheme with liquid water potential temperature and'/                        &
    11711137             '    total water content is used.' )
    117211383   FORMAT ( '--> Condensation is parameterized via 0% - or 100% scheme.' )
     
    11911157
    11921158
    1193 !------------------------------------------------------------------------------!
     1159!--------------------------------------------------------------------------------------------------!
    11941160! Description:
    11951161! ------------
    11961162!> Control of microphysics for all grid points
    1197 !------------------------------------------------------------------------------!
     1163!--------------------------------------------------------------------------------------------------!
    11981164    SUBROUTINE bcm_actions( location )
    11991165
     
    12301196
    12311197
    1232 !------------------------------------------------------------------------------!
     1198!--------------------------------------------------------------------------------------------------!
    12331199! Description:
    12341200! ------------
    12351201!> Control of microphysics for grid points i,j
    1236 !------------------------------------------------------------------------------!
     1202!--------------------------------------------------------------------------------------------------!
    12371203    SUBROUTINE bcm_actions_ij( i, j, location )
    12381204
    12391205
    1240     INTEGER(iwp),      INTENT(IN) ::  i         !< grid index in x-direction
    1241     INTEGER(iwp),      INTENT(IN) ::  j         !< grid index in y-direction
    12421206    CHARACTER (LEN=*), INTENT(IN) ::  location  !< call location string
    1243     INTEGER(iwp)  ::  dummy  !< call location string
    1244 
    1245     IF ( bulk_cloud_model    )   dummy = i + j
     1207
     1208    INTEGER(iwp)  ::  dummy                !< call location string
     1209
     1210    INTEGER(iwp), INTENT(IN) ::  i         !< grid index in x-direction
     1211    INTEGER(iwp), INTENT(IN) ::  j         !< grid index in y-direction
     1212
     1213
     1214    IF ( bulk_cloud_model )   dummy = i + j
    12461215
    12471216    SELECT CASE ( location )
     
    12761245
    12771246
    1278 !------------------------------------------------------------------------------!
     1247!--------------------------------------------------------------------------------------------------!
    12791248! Description:
    12801249! ------------
    12811250!> Control of microphysics for all grid points
    1282 !------------------------------------------------------------------------------!
     1251!--------------------------------------------------------------------------------------------------!
    12831252    SUBROUTINE bcm_non_advective_processes
    12841253
     
    12861255       CALL cpu_log( log_point(51), 'microphysics', 'start' )
    12871256
    1288        IF ( .NOT. microphysics_sat_adjust  .AND.                               &
    1289             ( intermediate_timestep_count == 1  .OR.                           &
    1290               call_microphysics_at_all_substeps ) )                            &
     1257       IF ( .NOT. microphysics_sat_adjust  .AND.  ( intermediate_timestep_count == 1  .OR.         &
     1258            call_microphysics_at_all_substeps ) )                                                  &
    12911259       THEN
    12921260
     
    12941262!
    12951263!--          Calculate vertical profile of the hydrostatic pressure (hyp)
    1296              hyp    = barometric_formula(zu, pt_surface *                      &
    1297                       exner_function(surface_pressure * 100.0_wp),             &
    1298                       surface_pressure * 100.0_wp)
     1264             hyp    = barometric_formula(zu, pt_surface *                                          &
     1265                      exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
    12991266             d_exner = exner_function_invers(hyp)
    13001267             exner = 1.0_wp / exner_function_invers(hyp)
     
    13021269!
    13031270!--          Compute reference density
    1304              rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp,      &
    1305                            pt_surface *                                        &
    1306                            exner_function(surface_pressure * 100.0_wp))
     1271             rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp,                          &
     1272                           pt_surface * exner_function(surface_pressure * 100.0_wp))
    13071273          ENDIF
    13081274
     
    13201286
    13211287!
    1322 !--       Compute cloud physics
     1288!--       Compute cloud physics.
    13231289!--       Here the the simple kessler scheme is used.
    13241290          IF ( microphysics_kessler )  THEN
     
    13271293
    13281294!
    1329 !--       Here the seifert beheng scheme is used. Cloud concentration is assumed to
    1330 !--       a constant value an qc a diagnostic value.
     1295!--       Here the seifert beheng scheme is used. Cloud concentration is assumed to a constant value
     1296!--       an qc a diagnostic value.
    13311297          ELSEIF ( microphysics_seifert  .AND.  .NOT. microphysics_morrison )  THEN
    13321298             CALL adjust_cloud
    1333              IF ( microphysics_ice_phase  .AND.  simulated_time > start_ice_microphysics )     &
    1334              CALL ice_nucleation
    1335              IF ( microphysics_ice_phase  .AND.  simulated_time > start_ice_microphysics )     &
    1336              CALL ice_deposition
     1299             IF ( microphysics_ice_phase  .AND.  simulated_time > start_ice_microphysics )  THEN
     1300                CALL ice_nucleation
     1301                CALL ice_deposition
     1302             ENDIF
    13371303             CALL autoconversion
    13381304             CALL accretion
     
    13411307             CALL sedimentation_rain
    13421308             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
    1343              IF ( microphysics_ice_phase  .AND.  simulated_time > start_ice_microphysics )     &
    1344              CALL adjust_ice
    1345              IF ( ice_crystal_sedimentation .AND. microphysics_ice_phase &
    1346                   .AND. simulated_time > start_ice_microphysics ) CALL sedimentation_ice
    1347 
    1348 !
    1349 !--       Here the morrison scheme is used. No rain processes are considered and qr and nr
    1350 !--       are not allocated
     1309             IF ( microphysics_ice_phase  .AND.  simulated_time > start_ice_microphysics )  THEN
     1310                CALL adjust_ice
     1311                IF ( ice_crystal_sedimentation )  CALL sedimentation_ice
     1312             ENDIF
     1313
     1314!
     1315!--       Here the morrison scheme is used. No rain processes are considered and qr and nr are not
     1316!--       allocated.
    13511317          ELSEIF ( microphysics_morrison_no_rain  .AND.  .NOT. microphysics_seifert )  THEN
    13521318             CALL activation
    13531319             CALL condensation
    1354              IF ( microphysics_ice_phase  .AND.  simulated_time > start_ice_microphysics )     &
    1355              CALL adjust_ice
    1356              IF ( microphysics_ice_phase  .AND.  simulated_time > start_ice_microphysics )     &
    1357              CALL ice_nucleation
    1358              IF ( microphysics_ice_phase  .AND.  simulated_time > start_ice_microphysics )     &
    1359              CALL ice_deposition
     1320             IF ( microphysics_ice_phase  .AND.  simulated_time > start_ice_microphysics )  THEN
     1321                CALL adjust_ice
     1322                CALL ice_nucleation
     1323                CALL ice_deposition
     1324             ENDIF
    13601325             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
    13611326
    13621327!
    1363 !--       Here the full morrison scheme is used and all processes of Seifert and Beheng are
    1364 !--       included
     1328!--       Here the full morrison scheme is used and all processes of Seifert and Beheng are included
    13651329          ELSEIF ( microphysics_morrison  .AND.  microphysics_seifert )  THEN
    13661330             CALL adjust_cloud
    13671331             CALL activation
    13681332             CALL condensation
    1369              IF ( microphysics_ice_phase  .AND.  simulated_time > start_ice_microphysics )     &
    1370              CALL adjust_ice
    1371              IF ( microphysics_ice_phase  .AND.  simulated_time > start_ice_microphysics )     &
    1372              CALL ice_nucleation
    1373              IF ( microphysics_ice_phase  .AND.  simulated_time > start_ice_microphysics )     &
    1374              CALL ice_deposition
     1333             IF ( microphysics_ice_phase  .AND.  simulated_time > start_ice_microphysics )  THEN
     1334                CALL adjust_ice
     1335                CALL ice_nucleation
     1336                CALL ice_deposition
     1337             ENDIF
    13751338             CALL autoconversion
    13761339             CALL accretion
     
    13911354
    13921355
    1393 !------------------------------------------------------------------------------!
     1356!--------------------------------------------------------------------------------------------------!
    13941357! Description:
    13951358! ------------
    13961359!> Control of microphysics for grid points i,j
    1397 !------------------------------------------------------------------------------!
     1360!--------------------------------------------------------------------------------------------------!
    13981361    SUBROUTINE bcm_non_advective_processes_ij( i, j )
    13991362
     
    14021365       INTEGER(iwp) ::  j                 !<
    14031366
    1404        IF ( .NOT. microphysics_sat_adjust  .AND.                               &
    1405             ( intermediate_timestep_count == 1  .OR.                           &
    1406               call_microphysics_at_all_substeps ) )                            &
     1367       IF ( .NOT. microphysics_sat_adjust  .AND.  ( intermediate_timestep_count == 1  .OR.         &
     1368            call_microphysics_at_all_substeps ) )                                                  &
    14071369       THEN
    14081370
     
    14101372!
    14111373!--          Calculate vertical profile of the hydrostatic pressure (hyp)
    1412              hyp    = barometric_formula(zu, pt_surface *                      &
    1413                       exner_function(surface_pressure * 100.0_wp),             &
    1414                       surface_pressure * 100.0_wp)
     1374             hyp    = barometric_formula(zu, pt_surface *                                          &
     1375                      exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
    14151376             d_exner = exner_function_invers(hyp)
    14161377             exner = 1.0_wp / exner_function_invers(hyp)
     
    14181379!
    14191380!--          Compute reference density
    1420              rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp,      &
    1421                            pt_surface *                                        &
    1422                            exner_function(surface_pressure * 100.0_wp))
     1381             rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp,                          &
     1382                           pt_surface * exner_function(surface_pressure * 100.0_wp))
    14231383          ENDIF
    14241384
     
    14421402
    14431403!
    1444 !--       Here the seifert beheng scheme is used. Cloud concentration is assumed to
    1445 !--       a constant value an qc a diagnostic value.
     1404!--       Here the seifert beheng scheme is used. Cloud concentration is assumed to a constant value
     1405!--       an qc a diagnostic value.
    14461406          ELSEIF ( microphysics_seifert  .AND.  .NOT. microphysics_morrison )  THEN
    14471407             CALL adjust_cloud_ij( i,j )
    1448              IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics )       &
     1408             IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics )  THEN
    14491409                 CALL ice_nucleation_ij( i,j )
    1450              IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics )       &
    14511410                 CALL ice_deposition_ij( i,j )
     1411             ENDIF
    14521412             CALL autoconversion_ij( i,j )
    14531413             CALL accretion_ij( i,j )
     
    14561416             CALL sedimentation_rain_ij( i,j )
    14571417             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud_ij( i,j )
    1458              IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics )       &
     1418             IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics )  THEN
    14591419                 CALL adjust_ice_ij ( i,j )
    1460              IF ( ice_crystal_sedimentation .AND. microphysics_ice_phase &
    1461                   .AND. simulated_time > start_ice_microphysics )  CALL sedimentation_ice_ij ( i,j )
    1462 !
    1463 !--       Here the morrison scheme is used. No rain processes are considered and qr and nr
    1464 !--       are not allocated
     1420                 IF ( ice_crystal_sedimentation )  CALL sedimentation_ice_ij ( i,j )
     1421             ENDIF
     1422
     1423!
     1424!--       Here the morrison scheme is used. No rain processes are considered and qr and nr are not
     1425!--       allocated.
    14651426          ELSEIF ( microphysics_morrison_no_rain  .AND.  .NOT. microphysics_seifert )  THEN
    14661427             CALL activation_ij( i,j )
    14671428             CALL condensation_ij( i,j )
    1468              IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics )       &
     1429             IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics )  THEN
    14691430                 CALL adjust_ice_ij ( i,j )
    1470              IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics )       &
    14711431                 CALL ice_nucleation_ij( i,j )
    1472              IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics )       &
    14731432                 CALL ice_deposition_ij( i,j )
     1433             ENDIF
    14741434             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud_ij( i,j )
    14751435
    14761436!
    1477 !--       Here the full morrison scheme is used and all processes of Seifert and Beheng are
    1478 !--       included
     1437!--       Here the full morrison scheme is used and all processes of Seifert and Beheng are included
    14791438          ELSEIF ( microphysics_morrison  .AND.  microphysics_seifert )  THEN
    14801439             CALL adjust_cloud_ij( i,j )
    14811440             CALL activation_ij( i,j )
    14821441             CALL condensation_ij( i,j )
    1483              IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics )       &
     1442             IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics )  THEN
    14841443                 CALL adjust_ice_ij ( i,j )
    1485              IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics )       &
    14861444                 CALL ice_nucleation_ij( i,j )
    1487              IF ( microphysics_ice_phase .AND. simulated_time > start_ice_microphysics )       &
    14881445                 CALL ice_deposition_ij( i,j )
     1446             ENDIF
    14891447             CALL autoconversion_ij( i,j )
    14901448             CALL accretion_ij( i,j )
     
    15011459
    15021460    END SUBROUTINE bcm_non_advective_processes_ij
    1503    
    1504    
    1505 !------------------------------------------------------------------------------!
     1461
     1462
     1463!--------------------------------------------------------------------------------------------------!
    15061464! Description:
    15071465! ------------
    15081466!> Control of microphysics for all grid points
    1509 !------------------------------------------------------------------------------!
     1467!--------------------------------------------------------------------------------------------------!
    15101468    SUBROUTINE bcm_exchange_horiz
    15111469
    1512        USE exchange_horiz_mod,                                                  &
     1470       USE exchange_horiz_mod,                                                                     &
    15131471           ONLY:  exchange_horiz
    15141472
    15151473
    1516        IF ( .NOT. microphysics_sat_adjust  .AND.                                &
    1517             ( intermediate_timestep_count == 1  .OR.                            &
    1518               call_microphysics_at_all_substeps ) )                             &
     1474       IF ( .NOT. microphysics_sat_adjust  .AND.  ( intermediate_timestep_count == 1  .OR.         &
     1475            call_microphysics_at_all_substeps ) )                                                  &
    15191476       THEN
    15201477          IF ( microphysics_morrison )  THEN
     
    15391496
    15401497
    1541 !------------------------------------------------------------------------------!
     1498!--------------------------------------------------------------------------------------------------!
    15421499! Description:
    15431500! ------------
    15441501!> Control of microphysics for all grid points
    1545 !------------------------------------------------------------------------------!
     1502!--------------------------------------------------------------------------------------------------!
    15461503    SUBROUTINE bcm_prognostic_equations
    15471504
     
    15541511
    15551512!
    1556 !--    If required, calculate prognostic equations for cloud water content
    1557 !--    and cloud drop concentration
     1513!--    If required, calculate prognostic equations for cloud water content and cloud drop
     1514!--    concentration.
    15581515       IF ( microphysics_morrison )  THEN
    15591516
     
    15811538             IF ( timestep_scheme(1:5) == 'runge' )  THEN
    15821539                IF ( ws_scheme_sca )  THEN
    1583                    CALL advec_s_ws( advc_flags_s, qc, 'qc',                    &
    1584                                     bc_dirichlet_l  .OR.  bc_radiation_l,      &
    1585                                     bc_dirichlet_n  .OR.  bc_radiation_n,      &
    1586                                     bc_dirichlet_r  .OR.  bc_radiation_r,      &
     1540                   CALL advec_s_ws( advc_flags_s, qc, 'qc',                                        &
     1541                                    bc_dirichlet_l  .OR.  bc_radiation_l,                          &
     1542                                    bc_dirichlet_n  .OR.  bc_radiation_n,                          &
     1543                                    bc_dirichlet_r  .OR.  bc_radiation_r,                          &
    15871544                                    bc_dirichlet_s  .OR.  bc_radiation_s )
    15881545                ELSE
     
    15941551          ENDIF
    15951552
    1596           CALL diffusion_s( qc,                                                &
    1597                             surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,          &
    1598                             surf_def_h(2)%qcsws,                               &
    1599                             surf_lsm_h%qcsws,    surf_usm_h%qcsws,             &
    1600                             surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,          &
    1601                             surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,          &
    1602                             surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,          &
    1603                             surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,          &
    1604                             surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,          &
     1553          CALL diffusion_s( qc,                                                                    &
     1554                            surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,                              &
     1555                            surf_def_h(2)%qcsws,                                                   &
     1556                            surf_lsm_h%qcsws,    surf_usm_h%qcsws,                                 &
     1557                            surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,                              &
     1558                            surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,                              &
     1559                            surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,                              &
     1560                            surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,                              &
     1561                            surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,                              &
    16051562                            surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
    16061563
     
    16121569                !DIR$ IVDEP
    16131570                DO  k = nzb+1, nzt
    1614                    qc_p(k,j,i) = qc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
    1615                                                       tsc(3) * tqc_m(k,j,i) )  &
    1616                                                     - tsc(5) * rdf_sc(k) *     &
    1617                                                                qc(k,j,i)       &
    1618                                              )                                 &
    1619                                     * MERGE( 1.0_wp, 0.0_wp,                   &
    1620                                              BTEST( wall_flags_total_0(k,j,i), 0 )   &
    1621                                           )
     1571                   qc_p(k,j,i) = qc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +                       &
     1572                                                      tsc(3) * tqc_m(k,j,i) )                      &
     1573                                                    - tsc(5) * rdf_sc(k) *                         &
     1574                                                               qc(k,j,i)                           &
     1575                                             )                                                     &
     1576                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 )   &
     1577                                         )
    16221578                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
    16231579                ENDDO
     
    16411597                   DO  j = nys, nyn
    16421598                      DO  k = nzb+1, nzt
    1643                          tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
    1644                                          + 5.3125_wp * tqc_m(k,j,i)
     1599                         tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) + 5.3125_wp * tqc_m(k,j,i)
    16451600                      ENDDO
    16461601                   ENDDO
     
    16731628             IF ( timestep_scheme(1:5) == 'runge' )  THEN
    16741629                IF ( ws_scheme_sca )  THEN
    1675                    CALL advec_s_ws( advc_flags_s, nc, 'nc',                    &
    1676                                     bc_dirichlet_l  .OR.  bc_radiation_l,      &
    1677                                     bc_dirichlet_n  .OR.  bc_radiation_n,      &
    1678                                     bc_dirichlet_r  .OR.  bc_radiation_r,      &
     1630                   CALL advec_s_ws( advc_flags_s, nc, 'nc',                                        &
     1631                                    bc_dirichlet_l  .OR.  bc_radiation_l,                          &
     1632                                    bc_dirichlet_n  .OR.  bc_radiation_n,                          &
     1633                                    bc_dirichlet_r  .OR.  bc_radiation_r,                          &
    16791634                                    bc_dirichlet_s  .OR.  bc_radiation_s )
    16801635                ELSE
     
    16861641          ENDIF
    16871642
    1688           CALL diffusion_s( nc,                                                &
    1689                             surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,          &
    1690                             surf_def_h(2)%ncsws,                               &
    1691                             surf_lsm_h%ncsws,    surf_usm_h%ncsws,             &
    1692                             surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,          &
    1693                             surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,          &
    1694                             surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,          &
    1695                             surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,          &
    1696                             surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,          &
     1643          CALL diffusion_s( nc,                                                                    &
     1644                            surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,                              &
     1645                            surf_def_h(2)%ncsws,                                                   &
     1646                            surf_lsm_h%ncsws,    surf_usm_h%ncsws,                                 &
     1647                            surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,                              &
     1648                            surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,                              &
     1649                            surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,                              &
     1650                            surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,                              &
     1651                            surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,                              &
    16971652                            surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
    16981653
     
    17041659                !DIR$ IVDEP
    17051660                DO  k = nzb+1, nzt
    1706                    nc_p(k,j,i) = nc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
    1707                                                       tsc(3) * tnc_m(k,j,i) )  &
    1708                                                     - tsc(5) * rdf_sc(k) *     &
    1709                                                                nc(k,j,i)       &
    1710                                              )                                 &
    1711                                    * MERGE( 1.0_wp, 0.0_wp,                    &
    1712                                              BTEST( wall_flags_total_0(k,j,i), 0 )   &
     1661                   nc_p(k,j,i) = nc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +                       &
     1662                                                      tsc(3) * tnc_m(k,j,i) )                      &
     1663                                                    - tsc(5) * rdf_sc(k) *                         &
     1664                                                               nc(k,j,i)                           &
     1665                                             )                                                     &
     1666                                   * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 )  &
    17131667                                          )
    17141668                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
     
    17331687                   DO  j = nys, nyn
    17341688                      DO  k = nzb+1, nzt
    1735                          tnc_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
    1736                                          + 5.3125_wp * tnc_m(k,j,i)
     1689                         tnc_m(k,j,i) =  -9.5625_wp * tend(k,j,i) + 5.3125_wp * tnc_m(k,j,i)
    17371690                      ENDDO
    17381691                   ENDDO
     
    17461699
    17471700!
    1748 !--    If required, calculate prognostic equations for ice crystal content
    1749 !--    and ice crystal concentration
     1701!--    If required, calculate prognostic equations for ice crystal content and ice crystal
     1702!--    concentration
    17501703       IF ( microphysics_ice_phase )  THEN
    17511704
     
    17731726             IF ( timestep_scheme(1:5) == 'runge' )  THEN
    17741727                IF ( ws_scheme_sca )  THEN
    1775                    CALL advec_s_ws( advc_flags_s, qi, 'qi',                    &
    1776                                     bc_dirichlet_l  .OR.  bc_radiation_l,      &
    1777                                     bc_dirichlet_n  .OR.  bc_radiation_n,      &
    1778                                     bc_dirichlet_r  .OR.  bc_radiation_r,      &
     1728                   CALL advec_s_ws( advc_flags_s, qi, 'qi',                                        &
     1729                                    bc_dirichlet_l  .OR.  bc_radiation_l,                          &
     1730                                    bc_dirichlet_n  .OR.  bc_radiation_n,                          &
     1731                                    bc_dirichlet_r  .OR.  bc_radiation_r,                          &
    17791732                                    bc_dirichlet_s  .OR.  bc_radiation_s )
    17801733                ELSE
     
    17861739          ENDIF
    17871740
    1788           CALL diffusion_s( qi,                                                &
    1789                             surf_def_h(0)%qisws, surf_def_h(1)%qisws,          &
    1790                             surf_def_h(2)%qisws,                               &
    1791                             surf_lsm_h%qisws,    surf_usm_h%qisws,             &
    1792                             surf_def_v(0)%qisws, surf_def_v(1)%qisws,          &
    1793                             surf_def_v(2)%qisws, surf_def_v(3)%qisws,          &
    1794                             surf_lsm_v(0)%qisws, surf_lsm_v(1)%qisws,          &
    1795                             surf_lsm_v(2)%qisws, surf_lsm_v(3)%qisws,          &
    1796                             surf_usm_v(0)%qisws, surf_usm_v(1)%qisws,          &
     1741          CALL diffusion_s( qi,                                                                    &
     1742                            surf_def_h(0)%qisws, surf_def_h(1)%qisws,                              &
     1743                            surf_def_h(2)%qisws,                                                   &
     1744                            surf_lsm_h%qisws,    surf_usm_h%qisws,                                 &
     1745                            surf_def_v(0)%qisws, surf_def_v(1)%qisws,                              &
     1746                            surf_def_v(2)%qisws, surf_def_v(3)%qisws,                              &
     1747                            surf_lsm_v(0)%qisws, surf_lsm_v(1)%qisws,                              &
     1748                            surf_lsm_v(2)%qisws, surf_lsm_v(3)%qisws,                              &
     1749                            surf_usm_v(0)%qisws, surf_usm_v(1)%qisws,                              &
    17971750                            surf_usm_v(2)%qisws, surf_usm_v(3)%qisws )
    17981751
     
    18041757                !DIR$ IVDEP
    18051758                DO  k = nzb+1, nzt
    1806                    qi_p(k,j,i) = qi(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
    1807                                                       tsc(3) * tqi_m(k,j,i) )  &
    1808                                                     - tsc(5) * rdf_sc(k) *     &
    1809                                                                qi(k,j,i)       &
    1810                                              )                                 &
    1811                                     * MERGE( 1.0_wp, 0.0_wp,                   &
    1812                                              BTEST( wall_flags_total_0(k,j,i), 0 )   &
     1759                   qi_p(k,j,i) = qi(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +                       &
     1760                                                      tsc(3) * tqi_m(k,j,i) )                      &
     1761                                                    - tsc(5) * rdf_sc(k) *                         &
     1762                                                               qi(k,j,i)                           &
     1763                                             )                                                     &
     1764                                    * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) &
    18131765                                          )
    18141766                   IF ( qi_p(k,j,i) < 0.0_wp )  qi_p(k,j,i) = 0.0_wp
     
    18331785                   DO  j = nys, nyn
    18341786                      DO  k = nzb+1, nzt
    1835                          tqi_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
    1836                                          + 5.3125_wp * tqi_m(k,j,i)
     1787                         tqi_m(k,j,i) =   -9.5625_wp * tend(k,j,i) + 5.3125_wp * tqi_m(k,j,i)
    18371788                      ENDDO
    18381789                   ENDDO
     
    18651816             IF ( timestep_scheme(1:5) == 'runge' )  THEN
    18661817                IF ( ws_scheme_sca )  THEN
    1867                    CALL advec_s_ws( advc_flags_s, ni, 'ni',                    &
    1868                                     bc_dirichlet_l  .OR.  bc_radiation_l,      &
    1869                                     bc_dirichlet_n  .OR.  bc_radiation_n,      &
    1870                                     bc_dirichlet_r  .OR.  bc_radiation_r,      &
     1818                   CALL advec_s_ws( advc_flags_s, ni, 'ni',                                        &
     1819                                    bc_dirichlet_l  .OR.  bc_radiation_l,                          &
     1820                                    bc_dirichlet_n  .OR.  bc_radiation_n,                          &
     1821                                    bc_dirichlet_r  .OR.  bc_radiation_r,                          &
    18711822                                    bc_dirichlet_s  .OR.  bc_radiation_s )
    18721823                ELSE
     
    18781829          ENDIF
    18791830
    1880           CALL diffusion_s( ni,                                                &
    1881                             surf_def_h(0)%nisws, surf_def_h(1)%nisws,          &
    1882                             surf_def_h(2)%nisws,                               &
    1883                             surf_lsm_h%nisws,    surf_usm_h%nisws,             &
    1884                             surf_def_v(0)%nisws, surf_def_v(1)%nisws,          &
    1885                             surf_def_v(2)%nisws, surf_def_v(3)%nisws,          &
    1886                             surf_lsm_v(0)%nisws, surf_lsm_v(1)%nisws,          &
    1887                             surf_lsm_v(2)%nisws, surf_lsm_v(3)%nisws,          &
    1888                             surf_usm_v(0)%nisws, surf_usm_v(1)%nisws,          &
     1831          CALL diffusion_s( ni,                                                                    &
     1832                            surf_def_h(0)%nisws, surf_def_h(1)%nisws,                              &
     1833                            surf_def_h(2)%nisws,                                                   &
     1834                            surf_lsm_h%nisws,    surf_usm_h%nisws,                                 &
     1835                            surf_def_v(0)%nisws, surf_def_v(1)%nisws,                              &
     1836                            surf_def_v(2)%nisws, surf_def_v(3)%nisws,                              &
     1837                            surf_lsm_v(0)%nisws, surf_lsm_v(1)%nisws,                              &
     1838                            surf_lsm_v(2)%nisws, surf_lsm_v(3)%nisws,                              &
     1839                            surf_usm_v(0)%nisws, surf_usm_v(1)%nisws,                              &
    18891840                            surf_usm_v(2)%nisws, surf_usm_v(3)%nisws )
    18901841
     
    18961847                !DIR$ IVDEP
    18971848                DO  k = nzb+1, nzt
    1898                    ni_p(k,j,i) = ni(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
    1899                                                       tsc(3) * tni_m(k,j,i) )  &
    1900                                                     - tsc(5) * rdf_sc(k) *     &
    1901                                                                ni(k,j,i)       &
    1902                                              )                                 &
    1903                                    * MERGE( 1.0_wp, 0.0_wp,                    &
    1904                                              BTEST( wall_flags_total_0(k,j,i), 0 )   &
     1849                   ni_p(k,j,i) = ni(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +                       &
     1850                                                      tsc(3) * tni_m(k,j,i) )                      &
     1851                                                    - tsc(5) * rdf_sc(k) *                         &
     1852                                                               ni(k,j,i)                           &
     1853                                             )                                                     &
     1854                                   * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 )  &
    19051855                                          )
    19061856                   IF ( ni_p(k,j,i) < 0.0_wp )  ni_p(k,j,i) = 0.0_wp
     
    19251875                   DO  j = nys, nyn
    19261876                      DO  k = nzb+1, nzt
    1927                          tni_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
    1928                                          + 5.3125_wp * tni_m(k,j,i)
     1877                         tni_m(k,j,i) =  -9.5625_wp * tend(k,j,i) + 5.3125_wp * tni_m(k,j,i)
    19291878                      ENDDO
    19301879                   ENDDO
     
    19381887
    19391888!
    1940 !--    If required, calculate prognostic equations for rain water content
    1941 !--    and rain drop concentration
     1889!--    If required, calculate prognostic equations for rain water content and rain drop
     1890!--    concentration.
    19421891       IF ( microphysics_seifert )  THEN
    19431892
     
    19651914             IF ( timestep_scheme(1:5) == 'runge' )  THEN
    19661915                IF ( ws_scheme_sca )  THEN
    1967                    CALL advec_s_ws( advc_flags_s, qr, 'qr',                    &
    1968                                     bc_dirichlet_l  .OR.  bc_radiation_l,      &
    1969                                     bc_dirichlet_n  .OR.  bc_radiation_n,      &
    1970                                     bc_dirichlet_r  .OR.  bc_radiation_r,      &
     1916                   CALL advec_s_ws( advc_flags_s, qr, 'qr',                                        &
     1917                                    bc_dirichlet_l  .OR.  bc_radiation_l,                          &
     1918                                    bc_dirichlet_n  .OR.  bc_radiation_n,                          &
     1919                                    bc_dirichlet_r  .OR.  bc_radiation_r,                          &
    19711920                                    bc_dirichlet_s  .OR.  bc_radiation_s )
    19721921                ELSE
     
    19781927          ENDIF
    19791928
    1980           CALL diffusion_s( qr,                                                &
    1981                             surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,          &
    1982                             surf_def_h(2)%qrsws,                               &
    1983                             surf_lsm_h%qrsws,    surf_usm_h%qrsws,             &
    1984                             surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,          &
    1985                             surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,          &
    1986                             surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,          &
    1987                             surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,          &
    1988                             surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,          &
     1929          CALL diffusion_s( qr,                                                                    &
     1930                            surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,                              &
     1931                            surf_def_h(2)%qrsws,                                                   &
     1932                            surf_lsm_h%qrsws,    surf_usm_h%qrsws,                                 &
     1933                            surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,                              &
     1934                            surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,                              &
     1935                            surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,                              &
     1936                            surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,                              &
     1937                            surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,                              &
    19891938                            surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
    19901939
     
    19961945                !DIR$ IVDEP
    19971946                DO  k = nzb+1, nzt
    1998                    qr_p(k,j,i) = qr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
    1999                                                       tsc(3) * tqr_m(k,j,i) )  &
    2000                                                     - tsc(5) * rdf_sc(k) *     &
    2001                                                                qr(k,j,i)       &
    2002                                              )                                 &
    2003                                     * MERGE( 1.0_wp, 0.0_wp,                   &
    2004                                              BTEST( wall_flags_total_0(k,j,i), 0 )   &
     1947                   qr_p(k,j,i) = qr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +                       &
     1948                                                      tsc(3) * tqr_m(k,j,i) )                      &
     1949                                                    - tsc(5) * rdf_sc(k) *                         &
     1950                                                               qr(k,j,i)                           &
     1951                                             )                                                     &
     1952                                    * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) &
    20051953                                          )
    20061954                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
     
    20251973                   DO  j = nys, nyn
    20261974                      DO  k = nzb+1, nzt
    2027                          tqr_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
    2028                                          + 5.3125_wp * tqr_m(k,j,i)
     1975                         tqr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) + 5.3125_wp * tqr_m(k,j,i)
    20291976                      ENDDO
    20301977                   ENDDO
     
    20572004             IF ( timestep_scheme(1:5) == 'runge' )  THEN
    20582005                IF ( ws_scheme_sca )  THEN
    2059                    CALL advec_s_ws( advc_flags_s, nr, 'nr',                    &
    2060                                     bc_dirichlet_l  .OR.  bc_radiation_l,      &
    2061                                     bc_dirichlet_n  .OR.  bc_radiation_n,      &
    2062                                     bc_dirichlet_r  .OR.  bc_radiation_r,      &
     2006                   CALL advec_s_ws( advc_flags_s, nr, 'nr',                                        &
     2007                                    bc_dirichlet_l  .OR.  bc_radiation_l,                          &
     2008                                    bc_dirichlet_n  .OR.  bc_radiation_n,                          &
     2009                                    bc_dirichlet_r  .OR.  bc_radiation_r,                          &
    20632010                                    bc_dirichlet_s  .OR.  bc_radiation_s )
    20642011                ELSE
     
    20702017          ENDIF
    20712018
    2072           CALL diffusion_s( nr,                                                &
    2073                             surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,          &
    2074                             surf_def_h(2)%nrsws,                               &
    2075                             surf_lsm_h%nrsws,    surf_usm_h%nrsws,             &
    2076                             surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,          &
    2077                             surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,          &
    2078                             surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,          &
    2079                             surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,          &
    2080                             surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,          &
     2019          CALL diffusion_s( nr,                                                                    &
     2020                            surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,                              &
     2021                            surf_def_h(2)%nrsws,                                                   &
     2022                            surf_lsm_h%nrsws,    surf_usm_h%nrsws,                                 &
     2023                            surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,                              &
     2024                            surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,                              &
     2025                            surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,                              &
     2026                            surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,                              &
     2027                            surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,                              &
    20812028                            surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
    20822029
     
    20882035                !DIR$ IVDEP
    20892036                DO  k = nzb+1, nzt
    2090                    nr_p(k,j,i) = nr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
    2091                                                       tsc(3) * tnr_m(k,j,i) )  &
    2092                                                     - tsc(5) * rdf_sc(k) *     &
    2093                                                                nr(k,j,i)       &
    2094                                              )                                 &
    2095                                    * MERGE( 1.0_wp, 0.0_wp,                    &
    2096                                              BTEST( wall_flags_total_0(k,j,i), 0 )   &
     2037                   nr_p(k,j,i) = nr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +                       &
     2038                                                      tsc(3) * tnr_m(k,j,i) )                      &
     2039                                                    - tsc(5) * rdf_sc(k) *                         &
     2040                                                               nr(k,j,i)                           &
     2041                                             )                                                     &
     2042                                   * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 )  &
    20972043                                          )
    20982044                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
     
    21322078
    21332079
    2134 !------------------------------------------------------------------------------!
     2080!--------------------------------------------------------------------------------------------------!
    21352081! Description:
    21362082! ------------
    21372083!> Control of microphysics for grid points i,j
    2138 !------------------------------------------------------------------------------!
     2084!--------------------------------------------------------------------------------------------------!
    21392085    SUBROUTINE bcm_prognostic_equations_ij( i, j, i_omp_start, tn )
    21402086
    21412087
    21422088       INTEGER(iwp), INTENT(IN) ::  i            !< grid index in x-direction
     2089       INTEGER(iwp), INTENT(IN) ::  i_omp_start  !< first loop index of i-loop in
     2090                                                 !< prognostic_equations
    21432091       INTEGER(iwp), INTENT(IN) ::  j            !< grid index in y-direction
    21442092       INTEGER(iwp)             ::  k            !< grid index in z-direction
    2145        INTEGER(iwp), INTENT(IN) ::  i_omp_start  !< first loop index of i-loop in prognostic_equations
    21462093       INTEGER(iwp), INTENT(IN) ::  tn           !< task number of openmp task
    21472094
    21482095!
    2149 !--    If required, calculate prognostic equations for cloud water content
    2150 !--    and cloud drop concentration
     2096!--    If required, calculate prognostic equations for cloud water content and cloud drop
     2097!--    concentration.
    21512098       IF ( microphysics_morrison )  THEN
    21522099!
     
    21562103          THEN
    21572104             IF ( ws_scheme_sca )  THEN
    2158                 CALL advec_s_ws( advc_flags_s, i, j, qc, 'qc', flux_s_qc,      &
    2159                                  diss_s_qc, flux_l_qc, diss_l_qc,              &
    2160                                  i_omp_start, tn,                              &
    2161                                  bc_dirichlet_l  .OR.  bc_radiation_l,         &
    2162                                  bc_dirichlet_n  .OR.  bc_radiation_n,         &
    2163                                  bc_dirichlet_r  .OR.  bc_radiation_r,         &
     2105                CALL advec_s_ws( advc_flags_s, i, j, qc, 'qc', flux_s_qc,                          &
     2106                                 diss_s_qc, flux_l_qc, diss_l_qc,                                  &
     2107                                 i_omp_start, tn,                                                  &
     2108                                 bc_dirichlet_l  .OR.  bc_radiation_l,                             &
     2109                                 bc_dirichlet_n  .OR.  bc_radiation_n,                             &
     2110                                 bc_dirichlet_r  .OR.  bc_radiation_r,                             &
    21642111                                 bc_dirichlet_s  .OR.  bc_radiation_s )
    21652112             ELSE
     
    21692116             CALL advec_s_up( i, j, qc )
    21702117          ENDIF
    2171           CALL diffusion_s( i, j, qc,                                   &
    2172                             surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,   &
    2173                             surf_def_h(2)%qcsws,                        &
    2174                             surf_lsm_h%qcsws,    surf_usm_h%qcsws,      &
    2175                             surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,   &
    2176                             surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,   &
    2177                             surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,   &
    2178                             surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,   &
    2179                             surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,   &
     2118          CALL diffusion_s( i, j, qc,                                                              &
     2119                            surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,                              &
     2120                            surf_def_h(2)%qcsws,                                                   &
     2121                            surf_lsm_h%qcsws,    surf_usm_h%qcsws,                                 &
     2122                            surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,                              &
     2123                            surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,                              &
     2124                            surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,                              &
     2125                            surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,                              &
     2126                            surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,                              &
    21802127                            surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
    21812128
     
    21832130!--       Prognostic equation for cloud water content
    21842131          DO  k = nzb+1, nzt
    2185              qc_p(k,j,i) = qc(k,j,i) + ( dt_3d *                         &
    2186                                                 ( tsc(2) * tend(k,j,i) + &
    2187                                                   tsc(3) * tqc_m(k,j,i) )&
    2188                                                 - tsc(5) * rdf_sc(k)     &
    2189                                                          * qc(k,j,i)     &
    2190                                        )                                 &
    2191                                  * MERGE( 1.0_wp, 0.0_wp,                &
    2192                                           BTEST( wall_flags_total_0(k,j,i), 0 )&
     2132             qc_p(k,j,i) = qc(k,j,i) + ( dt_3d *                                                   &
     2133                                                ( tsc(2) * tend(k,j,i) +                           &
     2134                                                  tsc(3) * tqc_m(k,j,i) )                          &
     2135                                                - tsc(5) * rdf_sc(k)                               &
     2136                                                         * qc(k,j,i)                               &
     2137                                       )                                                           &
     2138                                 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 )    &
    21932139                                        )
    21942140             IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
     
    22042150                      intermediate_timestep_count_max )  THEN
    22052151                DO  k = nzb+1, nzt
    2206                    tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
    2207                                      5.3125_wp * tqc_m(k,j,i)
     2152                   tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) + 5.3125_wp * tqc_m(k,j,i)
    22082153                ENDDO
    22092154             ENDIF
     
    22152160          IF ( timestep_scheme(1:5) == 'runge' )  THEN
    22162161             IF ( ws_scheme_sca )  THEN
    2217                 CALL advec_s_ws( advc_flags_s, i, j, nc, 'nc', flux_s_nc,      &
    2218                                  diss_s_nc, flux_l_nc, diss_l_nc,              &
    2219                                  i_omp_start, tn,                              &
    2220                                  bc_dirichlet_l  .OR.  bc_radiation_l,         &
    2221                                  bc_dirichlet_n  .OR.  bc_radiation_n,         &
    2222                                  bc_dirichlet_r  .OR.  bc_radiation_r,         &
     2162                CALL advec_s_ws( advc_flags_s, i, j, nc, 'nc', flux_s_nc,                          &
     2163                                 diss_s_nc, flux_l_nc, diss_l_nc,                                  &
     2164                                 i_omp_start, tn,                                                  &
     2165                                 bc_dirichlet_l  .OR.  bc_radiation_l,                             &
     2166                                 bc_dirichlet_n  .OR.  bc_radiation_n,                             &
     2167                                 bc_dirichlet_r  .OR.  bc_radiation_r,                             &
    22232168                                 bc_dirichlet_s  .OR.  bc_radiation_s )
    22242169             ELSE
     
    22282173             CALL advec_s_up( i, j, nc )
    22292174          ENDIF
    2230           CALL diffusion_s( i, j, nc,                                    &
    2231                             surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,    &
    2232                             surf_def_h(2)%ncsws,                         &
    2233                             surf_lsm_h%ncsws,    surf_usm_h%ncsws,       &
    2234                             surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,    &
    2235                             surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,    &
    2236                             surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,    &
    2237                             surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,    &
    2238                             surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,    &
     2175          CALL diffusion_s( i, j, nc,                                                              &
     2176                            surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,                              &
     2177                            surf_def_h(2)%ncsws,                                                   &
     2178                            surf_lsm_h%ncsws,    surf_usm_h%ncsws,                                 &
     2179                            surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,                              &
     2180                            surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,                              &
     2181                            surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,                              &
     2182                            surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,                              &
     2183                            surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,                              &
    22392184                            surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
    22402185
     
    22422187!--       Prognostic equation for cloud drop concentration
    22432188          DO  k = nzb+1, nzt
    2244              nc_p(k,j,i) = nc(k,j,i) + ( dt_3d *                         &
    2245                                                 ( tsc(2) * tend(k,j,i) + &
    2246                                                   tsc(3) * tnc_m(k,j,i) )&
    2247                                                 - tsc(5) * rdf_sc(k)     &
    2248                                                          * nc(k,j,i)     &
    2249                                        )                                 &
    2250                                  * MERGE( 1.0_wp, 0.0_wp,                &
    2251                                           BTEST( wall_flags_total_0(k,j,i), 0 )&
     2189             nc_p(k,j,i) = nc(k,j,i) + ( dt_3d *                                                   &
     2190                                                ( tsc(2) * tend(k,j,i) +                           &
     2191                                                  tsc(3) * tnc_m(k,j,i) )                          &
     2192                                                - tsc(5) * rdf_sc(k)                               &
     2193                                                         * nc(k,j,i)                               &
     2194                                       )                                                           &
     2195                                 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 )    &
    22522196                                        )
    22532197             IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
     
    22632207                      intermediate_timestep_count_max )  THEN
    22642208                DO  k = nzb+1, nzt
    2265                    tnc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
    2266                                      5.3125_wp * tnc_m(k,j,i)
     2209                   tnc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) + 5.3125_wp * tnc_m(k,j,i)
    22672210                ENDDO
    22682211             ENDIF
     
    22722215
    22732216!
    2274 !--    If required, calculate prognostic equations for ice crystal mixing ratio
    2275 !--    and ice crystal concentration
     2217!--    If required, calculate prognostic equations for ice crystal mixing ratio and ice crystal
     2218!--    concentration.
    22762219       IF ( microphysics_ice_phase )  THEN
    22772220!
     
    22812224          THEN
    22822225             IF ( ws_scheme_sca )  THEN
    2283                 CALL advec_s_ws( advc_flags_s, i, j, qi, 'qi', flux_s_qi,      &
    2284                                  diss_s_qi, flux_l_qi, diss_l_qi,              &
    2285                                  i_omp_start, tn,                              &
    2286                                  bc_dirichlet_l  .OR.  bc_radiation_l,         &
    2287                                  bc_dirichlet_n  .OR.  bc_radiation_n,         &
    2288                                  bc_dirichlet_r  .OR.  bc_radiation_r,         &
     2226                CALL advec_s_ws( advc_flags_s, i, j, qi, 'qi', flux_s_qi,                          &
     2227                                 diss_s_qi, flux_l_qi, diss_l_qi,                                  &
     2228                                 i_omp_start, tn,                                                  &
     2229                                 bc_dirichlet_l  .OR.  bc_radiation_l,                             &
     2230                                 bc_dirichlet_n  .OR.  bc_radiation_n,                             &
     2231                                 bc_dirichlet_r  .OR.  bc_radiation_r,                             &
    22892232                                 bc_dirichlet_s  .OR.  bc_radiation_s )
    22902233             ELSE
     
    22942237             CALL advec_s_up( i, j, qi )
    22952238          ENDIF
    2296           CALL diffusion_s( i, j, qi,                                   &
    2297                             surf_def_h(0)%qisws, surf_def_h(1)%qisws,   &
    2298                             surf_def_h(2)%qisws,                        &
    2299                             surf_lsm_h%qisws,    surf_usm_h%qisws,      &
    2300                             surf_def_v(0)%qisws, surf_def_v(1)%qisws,   &
    2301                             surf_def_v(2)%qisws, surf_def_v(3)%qisws,   &
    2302                             surf_lsm_v(0)%qisws, surf_lsm_v(1)%qisws,   &
    2303                             surf_lsm_v(2)%qisws, surf_lsm_v(3)%qisws,   &
    2304                             surf_usm_v(0)%qisws, surf_usm_v(1)%qisws,   &
     2239          CALL diffusion_s( i, j, qi,                                                              &
     2240                            surf_def_h(0)%qisws, surf_def_h(1)%qisws,                              &
     2241                            surf_def_h(2)%qisws,                                                   &
     2242                            surf_lsm_h%qisws,    surf_usm_h%qisws,                                 &
     2243                            surf_def_v(0)%qisws, surf_def_v(1)%qisws,                              &
     2244                            surf_def_v(2)%qisws, surf_def_v(3)%qisws,                              &
     2245                            surf_lsm_v(0)%qisws, surf_lsm_v(1)%qisws,                              &
     2246                            surf_lsm_v(2)%qisws, surf_lsm_v(3)%qisws,                              &
     2247                            surf_usm_v(0)%qisws, surf_usm_v(1)%qisws,                              &
    23052248                            surf_usm_v(2)%qisws, surf_usm_v(3)%qisws )
    23062249
     
    23082251!--       Prognostic equation for ice crystal mixing ratio
    23092252          DO  k = nzb+1, nzt
    2310              qi_p(k,j,i) = qi(k,j,i) + ( dt_3d *                         &
    2311                                                 ( tsc(2) * tend(k,j,i) + &
    2312                                                   tsc(3) * tqi_m(k,j,i) )&
    2313                                                 - tsc(5) * rdf_sc(k)     &
    2314                                                          * qi(k,j,i)     &
    2315                                        )                                 &
    2316                                  * MERGE( 1.0_wp, 0.0_wp,                &
    2317                                           BTEST( wall_flags_total_0(k,j,i), 0 )&
     2253             qi_p(k,j,i) = qi(k,j,i) + ( dt_3d *                                                   &
     2254                                                ( tsc(2) * tend(k,j,i) +                           &
     2255                                                  tsc(3) * tqi_m(k,j,i) )                          &
     2256                                                - tsc(5) * rdf_sc(k)                               &
     2257                                                         * qi(k,j,i)                               &
     2258                                       )                                                           &
     2259                                 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 )    &
    23182260                                        )
    23192261             IF ( qi_p(k,j,i) < 0.0_wp )  qi_p(k,j,i) = 0.0_wp
     
    23292271                      intermediate_timestep_count_max )  THEN
    23302272                DO  k = nzb+1, nzt
    2331                    tqi_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
    2332                                      5.3125_wp * tqi_m(k,j,i)
     2273                   tqi_m(k,j,i) =   -9.5625_wp * tend(k,j,i) + 5.3125_wp * tqi_m(k,j,i)
    23332274                ENDDO
    23342275             ENDIF
     
    23402281          IF ( timestep_scheme(1:5) == 'runge' )  THEN
    23412282             IF ( ws_scheme_sca )  THEN
    2342                 CALL advec_s_ws( advc_flags_s, i, j, ni, 'ni', flux_s_ni,      &
    2343                                  diss_s_ni, flux_l_ni, diss_l_ni,              &
    2344                                  i_omp_start, tn,                              &
    2345                                  bc_dirichlet_l  .OR.  bc_radiation_l,         &
    2346                                  bc_dirichlet_n  .OR.  bc_radiation_n,         &
    2347                                  bc_dirichlet_r  .OR.  bc_radiation_r,         &
     2283                CALL advec_s_ws( advc_flags_s, i, j, ni, 'ni', flux_s_ni,                          &
     2284                                 diss_s_ni, flux_l_ni, diss_l_ni,                                  &
     2285                                 i_omp_start, tn,                                                  &
     2286                                 bc_dirichlet_l  .OR.  bc_radiation_l,                             &
     2287                                 bc_dirichlet_n  .OR.  bc_radiation_n,                             &
     2288                                 bc_dirichlet_r  .OR.  bc_radiation_r,                             &
    23482289                                 bc_dirichlet_s  .OR.  bc_radiation_s )
    23492290             ELSE
     
    23532294             CALL advec_s_up( i, j, ni )
    23542295          ENDIF
    2355           CALL diffusion_s( i, j, ni,                                    &
    2356                             surf_def_h(0)%nisws, surf_def_h(1)%nisws,    &
    2357                             surf_def_h(2)%nisws,                         &
    2358                             surf_lsm_h%nisws,    surf_usm_h%nisws,       &
    2359                             surf_def_v(0)%nisws, surf_def_v(1)%nisws,    &
    2360                             surf_def_v(2)%nisws, surf_def_v(3)%nisws,    &
    2361                             surf_lsm_v(0)%nisws, surf_lsm_v(1)%nisws,    &
    2362                             surf_lsm_v(2)%nisws, surf_lsm_v(3)%nisws,    &
    2363                             surf_usm_v(0)%nisws, surf_usm_v(1)%nisws,    &
     2296          CALL diffusion_s( i, j, ni,                                                              &
     2297                            surf_def_h(0)%nisws, surf_def_h(1)%nisws,                              &
     2298                            surf_def_h(2)%nisws,                                                   &
     2299                            surf_lsm_h%nisws,    surf_usm_h%nisws,                                 &
     2300                            surf_def_v(0)%nisws, surf_def_v(1)%nisws,                              &
     2301                            surf_def_v(2)%nisws, surf_def_v(3)%nisws,                              &
     2302                            surf_lsm_v(0)%nisws, surf_lsm_v(1)%nisws,                              &
     2303                            surf_lsm_v(2)%nisws, surf_lsm_v(3)%nisws,                              &
     2304                            surf_usm_v(0)%nisws, surf_usm_v(1)%nisws,                              &
    23642305                            surf_usm_v(2)%nisws, surf_usm_v(3)%nisws )
    23652306
     
    23672308!--       Prognostic equation for ice crystal concentration
    23682309          DO  k = nzb+1, nzt
    2369              ni_p(k,j,i) = ni(k,j,i) + ( dt_3d *                         &
    2370                                                 ( tsc(2) * tend(k,j,i) + &
    2371                                                   tsc(3) * tni_m(k,j,i) )&
    2372                                                 - tsc(5) * rdf_sc(k)     &
    2373                                                          * ni(k,j,i)     &
    2374                                        )                                 &
    2375                                  * MERGE( 1.0_wp, 0.0_wp,                &
    2376                                           BTEST( wall_flags_total_0(k,j,i), 0 )&
     2310             ni_p(k,j,i) = ni(k,j,i) + ( dt_3d *                                                   &
     2311                                                ( tsc(2) * tend(k,j,i) +                           &
     2312                                                  tsc(3) * tni_m(k,j,i) )                          &
     2313                                                - tsc(5) * rdf_sc(k)                               &
     2314                                                         * ni(k,j,i)                               &
     2315                                       )                                                           &
     2316                                 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 )    &
    23772317                                        )
    23782318             IF ( ni_p(k,j,i) < 0.0_wp )  ni_p(k,j,i) = 0.0_wp
     
    23882328                      intermediate_timestep_count_max )  THEN
    23892329                DO  k = nzb+1, nzt
    2390                    tni_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
    2391                                      5.3125_wp * tni_m(k,j,i)
     2330                   tni_m(k,j,i) =   -9.5625_wp * tend(k,j,i) + 5.3125_wp * tni_m(k,j,i)
    23922331                ENDDO
    23932332             ENDIF
     
    23972336
    23982337!
    2399 !--    If required, calculate prognostic equations for rain water content
    2400 !--    and rain drop concentration
     2338!--    If required, calculate prognostic equations for rain water content and rain drop
     2339!--    concentration
    24012340       IF ( microphysics_seifert )  THEN
    24022341!
    2403 !--             Calculate prognostic equation for rain water content
     2342!--       Calculate prognostic equation for rain water content
    24042343          tend(:,j,i) = 0.0_wp
    24052344          IF ( timestep_scheme(1:5) == 'runge' ) &
    24062345          THEN
    24072346             IF ( ws_scheme_sca )  THEN
    2408                 CALL advec_s_ws( advc_flags_s, i, j, qr, 'qr', flux_s_qr,      &
    2409                                  diss_s_qr, flux_l_qr, diss_l_qr,              &
    2410                                  i_omp_start, tn,                              &
    2411                                  bc_dirichlet_l  .OR.  bc_radiation_l,         &
    2412                                  bc_dirichlet_n  .OR.  bc_radiation_n,         &
    2413                                  bc_dirichlet_r  .OR.  bc_radiation_r,         &
     2347                CALL advec_s_ws( advc_flags_s, i, j, qr, 'qr', flux_s_qr,                          &
     2348                                 diss_s_qr, flux_l_qr, diss_l_qr,                                  &
     2349                                 i_omp_start, tn,                                                  &
     2350                                 bc_dirichlet_l  .OR.  bc_radiation_l,                             &
     2351                                 bc_dirichlet_n  .OR.  bc_radiation_n,                             &
     2352                                 bc_dirichlet_r  .OR.  bc_radiation_r,                             &
    24142353                                 bc_dirichlet_s  .OR.  bc_radiation_s )
    24152354             ELSE
     
    24192358             CALL advec_s_up( i, j, qr )
    24202359          ENDIF
    2421           CALL diffusion_s( i, j, qr,                                   &
    2422                             surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,   &
    2423                             surf_def_h(2)%qrsws,                        &
    2424                             surf_lsm_h%qrsws,    surf_usm_h%qrsws,      &
    2425                             surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,   &
    2426                             surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,   &
    2427                             surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,   &
    2428                             surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,   &
    2429                             surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,   &
     2360          CALL diffusion_s( i, j, qr,                                                              &
     2361                            surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,                              &
     2362                            surf_def_h(2)%qrsws,                                                   &
     2363                            surf_lsm_h%qrsws,    surf_usm_h%qrsws,                                 &
     2364                            surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,                              &
     2365                            surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,                              &
     2366                            surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,                              &
     2367                            surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,                              &
     2368                            surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,                              &
    24302369                            surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
    24312370
     
    24332372!--       Prognostic equation for rain water content
    24342373          DO  k = nzb+1, nzt
    2435              qr_p(k,j,i) = qr(k,j,i) + ( dt_3d *                         &
    2436                                                 ( tsc(2) * tend(k,j,i) + &
    2437                                                   tsc(3) * tqr_m(k,j,i) )&
    2438                                                 - tsc(5) * rdf_sc(k)     &
    2439                                                          * qr(k,j,i)     &
    2440                                        )                                 &
    2441                                  * MERGE( 1.0_wp, 0.0_wp,                &
    2442                                           BTEST( wall_flags_total_0(k,j,i), 0 )&
     2374             qr_p(k,j,i) = qr(k,j,i) + ( dt_3d *                                                   &
     2375                                                ( tsc(2) * tend(k,j,i) +                           &
     2376                                                  tsc(3) * tqr_m(k,j,i) )                          &
     2377                                                - tsc(5) * rdf_sc(k)                               &
     2378                                                         * qr(k,j,i)                               &
     2379                                       )                                                           &
     2380                                 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 )    &
    24432381                                        )
    24442382             IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
     
    24542392                      intermediate_timestep_count_max )  THEN
    24552393                DO  k = nzb+1, nzt
    2456                    tqr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
    2457                                      5.3125_wp * tqr_m(k,j,i)
     2394                   tqr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) + 5.3125_wp * tqr_m(k,j,i)
    24582395                ENDDO
    24592396             ENDIF
     
    24652402          IF ( timestep_scheme(1:5) == 'runge' )  THEN
    24662403             IF ( ws_scheme_sca )  THEN
    2467                 CALL advec_s_ws( advc_flags_s, i, j, nr, 'nr', flux_s_nr,      &
    2468                                  diss_s_nr, flux_l_nr, diss_l_nr,              &
    2469                                  i_omp_start, tn,                              &
    2470                                  bc_dirichlet_l  .OR.  bc_radiation_l,         &
    2471                                  bc_dirichlet_n  .OR.  bc_radiation_n,         &
    2472                                  bc_dirichlet_r  .OR.  bc_radiation_r,         &
     2404                CALL advec_s_ws( advc_flags_s, i, j, nr, 'nr', flux_s_nr,                          &
     2405                                 diss_s_nr, flux_l_nr, diss_l_nr,                                  &
     2406                                 i_omp_start, tn,                                                  &
     2407                                 bc_dirichlet_l  .OR.  bc_radiation_l,                             &
     2408                                 bc_dirichlet_n  .OR.  bc_radiation_n,                             &
     2409                                 bc_dirichlet_r  .OR.  bc_radiation_r,                             &
    24732410                                 bc_dirichlet_s  .OR.  bc_radiation_s )
    24742411             ELSE
     
    24782415             CALL advec_s_up( i, j, nr )
    24792416          ENDIF
    2480           CALL diffusion_s( i, j, nr,                                    &
    2481                             surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,    &
    2482                             surf_def_h(2)%nrsws,                         &
    2483                             surf_lsm_h%nrsws,    surf_usm_h%nrsws,       &
    2484                             surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,    &
    2485                             surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,    &
    2486                             surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,    &
    2487                             surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,    &
    2488                             surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,    &
     2417          CALL diffusion_s( i, j, nr,                                                              &
     2418                            surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,                              &
     2419                            surf_def_h(2)%nrsws,                                                   &
     2420                            surf_lsm_h%nrsws,    surf_usm_h%nrsws,                                 &
     2421                            surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,                              &
     2422                            surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,                              &
     2423                            surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,                              &
     2424                            surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,                              &
     2425                            surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,                              &
    24892426                            surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
    24902427
     
    24922429!--       Prognostic equation for rain drop concentration
    24932430          DO  k = nzb+1, nzt
    2494              nr_p(k,j,i) = nr(k,j,i) + ( dt_3d *                         &
    2495                                                 ( tsc(2) * tend(k,j,i) + &
    2496                                                   tsc(3) * tnr_m(k,j,i) )&
    2497                                                 - tsc(5) * rdf_sc(k)     &
    2498                                                          * nr(k,j,i)     &
    2499                                        )                                 &
    2500                                  * MERGE( 1.0_wp, 0.0_wp,                &
    2501                                           BTEST( wall_flags_total_0(k,j,i), 0 )&
     2431             nr_p(k,j,i) = nr(k,j,i) + ( dt_3d *                                                   &
     2432                                                ( tsc(2) * tend(k,j,i) +                           &
     2433                                                  tsc(3) * tnr_m(k,j,i) )                          &
     2434                                                - tsc(5) * rdf_sc(k)                               &
     2435                                                         * nr(k,j,i)                               &
     2436                                       )                                                           &
     2437                                 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 )    &
    25022438                                        )
    25032439             IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
     
    25132449                      intermediate_timestep_count_max )  THEN
    25142450                DO  k = nzb+1, nzt
    2515                    tnr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
    2516                                      5.3125_wp * tnr_m(k,j,i)
     2451                   tnr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) + 5.3125_wp * tnr_m(k,j,i)
    25172452                ENDDO
    25182453             ENDIF
     
    25242459
    25252460
    2526 !------------------------------------------------------------------------------!
     2461!--------------------------------------------------------------------------------------------------!
    25272462! Description:
    25282463! ------------
    25292464!> Swapping of timelevels
    2530 !------------------------------------------------------------------------------!
     2465!--------------------------------------------------------------------------------------------------!
    25312466    SUBROUTINE bcm_swap_timelevel ( mod_count )
    25322467
     
    25772512
    25782513
    2579 !------------------------------------------------------------------------------!
     2514!--------------------------------------------------------------------------------------------------!
    25802515! Description: Boundary conditions of the bulk cloud module variables
    2581 !------------------------------------------------------------------------------!
     2516!--------------------------------------------------------------------------------------------------!
    25822517    SUBROUTINE bcm_boundary_conditions
    25832518
     
    25922527       IF ( microphysics_morrison )  THEN
    25932528!
    2594 !--       Surface conditions cloud water (Dirichlet)
    2595 !--       Run loop over all non-natural and natural walls. Note, in wall-datatype
    2596 !--       the k coordinate belongs to the atmospheric grid point, therefore, set
    2597 !--       qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls
     2529!--       Surface conditions cloud water (Dirichlet).
     2530!--       Run loop over all non-natural and natural walls. Note, in wall-datatype the k coordinate
     2531!--       belongs to the atmospheric grid point, therefore, set qr_p and nr_p at upward (k-1) and
     2532!--       downward-facing (k+1) walls.
    25982533          DO  l = 0, 1
    25992534          !$OMP PARALLEL DO PRIVATE( i, j, k )
     
    26152550       IF ( microphysics_ice_phase )  THEN
    26162551!
    2617 !--       Surface conditions ice crysral (Dirichlet)
    2618 !--       Run loop over all non-natural and natural walls. Note, in wall-datatype
    2619 !--       the k coordinate belongs to the atmospheric grid point, therefore, set
    2620 !--       qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls
     2552!--       Surface conditions ice crysral (Dirichlet).
     2553!--       Run loop over all non-natural and natural walls. Note, in wall-datatype the k coordinate
     2554!--       belongs to the atmospheric grid point, therefore, set qr_p and nr_p at upward (k-1) and
     2555!--       downward-facing (k+1) walls
    26212556          DO  l = 0, 1
    26222557          !$OMP PARALLEL DO PRIVATE( i, j, k )
     
    26392574       IF ( microphysics_seifert )  THEN
    26402575!
    2641 !--       Surface conditions rain water (Dirichlet)
    2642 !--       Run loop over all non-natural and natural walls. Note, in wall-datatype
    2643 !--       the k coordinate belongs to the atmospheric grid point, therefore, set
    2644 !--       qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls
     2576!--       Surface conditions rain water (Dirichlet).
     2577!--       Run loop over all non-natural and natural walls. Note, in wall-datatype the k coordinate
     2578!--       belongs to the atmospheric grid point, therefore, set qr_p and nr_p at upward (k-1) and
     2579!--       downward-facing (k+1) walls
    26452580          DO  l = 0, 1
    26462581          !$OMP PARALLEL DO PRIVATE( i, j, k )
     
    26622597!
    26632598!--    Lateral boundary conditions for scalar quantities at the outflow.
    2664 !--    Lateral oundary conditions for TKE and dissipation are set
    2665 !--    in tcm_boundary_conds.
     2599!--    Lateral oundary conditions for TKE and dissipation are set in tcm_boundary_conds.
    26662600       IF ( bc_radiation_s )  THEN
    26672601          IF ( microphysics_morrison )  THEN
     
    27202654    END SUBROUTINE bcm_boundary_conditions
    27212655
    2722 !------------------------------------------------------------------------------!
     2656!--------------------------------------------------------------------------------------------------!
    27232657!
    27242658! Description:
    27252659! ------------
    27262660!> Subroutine for averaging 3D data
    2727 !------------------------------------------------------------------------------!
     2661!--------------------------------------------------------------------------------------------------!
    27282662    SUBROUTINE bcm_3d_data_averaging( mode, variable )
    27292663
    2730        USE control_parameters,                                                 &
     2664       USE control_parameters,                                                                     &
    27312665           ONLY:  average_count_3d
    27322666
     
    27902724
    27912725             CASE ( 'nc' )
    2792                 IF ( ALLOCATED( nc_av ) ) THEN
     2726                IF ( ALLOCATED( nc_av ) )  THEN
    27932727                   DO  i = nxlg, nxrg
    27942728                      DO  j = nysg, nyng
     
    28012735
    28022736             CASE ( 'nr' )
    2803                 IF ( ALLOCATED( nr_av ) ) THEN
     2737                IF ( ALLOCATED( nr_av ) )  THEN
    28042738                   DO  i = nxlg, nxrg
    28052739                      DO  j = nysg, nyng
     
    28122746
    28132747             CASE ( 'prr' )
    2814                 IF ( ALLOCATED( prr_av ) ) THEN
     2748                IF ( ALLOCATED( prr_av ) )  THEN
    28152749                   DO  i = nxlg, nxrg
    28162750                      DO  j = nysg, nyng
     
    28232757
    28242758             CASE ( 'qc' )
    2825                 IF ( ALLOCATED( qc_av ) ) THEN
     2759                IF ( ALLOCATED( qc_av ) )  THEN
    28262760                   DO  i = nxlg, nxrg
    28272761                      DO  j = nysg, nyng
     
    28342768
    28352769             CASE ( 'ql' )
    2836                 IF ( ALLOCATED( ql_av ) ) THEN
     2770                IF ( ALLOCATED( ql_av ) )  THEN
    28372771                   DO  i = nxlg, nxrg
    28382772                      DO  j = nysg, nyng
     
    28452779
    28462780             CASE ( 'qr' )
    2847                 IF ( ALLOCATED( qr_av ) ) THEN
     2781                IF ( ALLOCATED( qr_av ) )  THEN
    28482782                   DO  i = nxlg, nxrg
    28492783                      DO  j = nysg, nyng
     
    28652799
    28662800             CASE ( 'nc' )
    2867                 IF ( ALLOCATED( nc_av ) ) THEN
     2801                IF ( ALLOCATED( nc_av ) )  THEN
    28682802                   DO  i = nxlg, nxrg
    28692803                      DO  j = nysg, nyng
    28702804                         DO  k = nzb, nzt+1
    2871                             nc_av(k,j,i) = nc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2805                            nc_av(k,j,i) = nc_av(k,j,i) / REAL( average_count_3d, KIND = wp )
    28722806                         ENDDO
    28732807                      ENDDO
     
    28762810
    28772811             CASE ( 'nr' )
    2878                 IF ( ALLOCATED( nr_av ) ) THEN
     2812                IF ( ALLOCATED( nr_av ) )  THEN
    28792813                   DO  i = nxlg, nxrg
    28802814                      DO  j = nysg, nyng
    28812815                         DO  k = nzb, nzt+1
    2882                             nr_av(k,j,i) = nr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2816                            nr_av(k,j,i) = nr_av(k,j,i) / REAL( average_count_3d, KIND = wp )
    28832817                         ENDDO
    28842818                      ENDDO
     
    28872821
    28882822             CASE ( 'prr' )
    2889                 IF ( ALLOCATED( prr_av ) ) THEN
     2823                IF ( ALLOCATED( prr_av ) )  THEN
    28902824                   DO  i = nxlg, nxrg
    28912825                      DO  j = nysg, nyng
    28922826                         DO  k = nzb, nzt+1
    2893                             prr_av(k,j,i) = prr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2827                            prr_av(k,j,i) = prr_av(k,j,i) / REAL( average_count_3d, KIND = wp )
    28942828                         ENDDO
    28952829                      ENDDO
     
    28982832
    28992833             CASE ( 'qc' )
    2900                 IF ( ALLOCATED( qc_av ) ) THEN
     2834                IF ( ALLOCATED( qc_av ) )  THEN
    29012835                   DO  i = nxlg, nxrg
    29022836                      DO  j = nysg, nyng
    29032837                         DO  k = nzb, nzt+1
    2904                             qc_av(k,j,i) = qc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2838                            qc_av(k,j,i) = qc_av(k,j,i) / REAL( average_count_3d, KIND = wp )
    29052839                         ENDDO
    29062840                      ENDDO
     
    29092843
    29102844             CASE ( 'ql' )
    2911                 IF ( ALLOCATED( ql_av ) ) THEN
     2845                IF ( ALLOCATED( ql_av ) )  THEN
    29122846                   DO  i = nxlg, nxrg
    29132847                      DO  j = nysg, nyng
    29142848                         DO  k = nzb, nzt+1
    2915                             ql_av(k,j,i) = ql_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2849                            ql_av(k,j,i) = ql_av(k,j,i) / REAL( average_count_3d, KIND = wp )
    29162850                         ENDDO
    29172851                      ENDDO
     
    29202854
    29212855             CASE ( 'qr' )
    2922                 IF ( ALLOCATED( qr_av ) ) THEN
     2856                IF ( ALLOCATED( qr_av ) )  THEN
    29232857                   DO  i = nxlg, nxrg
    29242858                      DO  j = nysg, nyng
    29252859                         DO  k = nzb, nzt+1
    2926                             qr_av(k,j,i) = qr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2860                            qr_av(k,j,i) = qr_av(k,j,i) / REAL( average_count_3d, KIND = wp )
    29272861                         ENDDO
    29282862                      ENDDO
     
    29402874
    29412875
    2942 !------------------------------------------------------------------------------!
     2876!--------------------------------------------------------------------------------------------------!
    29432877! Description:
    29442878! ------------
    29452879!> Define 2D output variables.
    2946 !------------------------------------------------------------------------------!
    2947  SUBROUTINE bcm_data_output_2d( av, variable, found, grid, mode, local_pf,     &
    2948                                 two_d, nzb_do, nzt_do )
     2880!--------------------------------------------------------------------------------------------------!
     2881 SUBROUTINE bcm_data_output_2d( av, variable, found, grid, mode, local_pf, two_d, nzb_do, nzt_do )
    29492882
    29502883
     
    29522885
    29532886    CHARACTER (LEN=*), INTENT(INOUT) ::  grid       !< name of vertical grid
    2954     CHARACTER (LEN=*), INTENT(IN) ::  mode       !< either 'xy', 'xz' or 'yz'
    2955     CHARACTER (LEN=*), INTENT(IN) ::  variable   !< name of variable
     2887    CHARACTER (LEN=*), INTENT(IN)    ::  mode       !< either 'xy', 'xz' or 'yz'
     2888    CHARACTER (LEN=*), INTENT(IN)    ::  variable   !< name of variable
    29562889
    29572890    INTEGER(iwp), INTENT(IN) ::  av        !< flag for (non-)average output
     
    29602893
    29612894    INTEGER(iwp) ::  flag_nr   !< number of masking flag
    2962 
    29632895    INTEGER(iwp) ::  i         !< loop index along x-direction
    29642896    INTEGER(iwp) ::  j         !< loop index along y-direction
    29652897    INTEGER(iwp) ::  k         !< loop index along z-direction
    29662898
     2899    LOGICAL ::  resorted  !< flag if output is already resorted
     2900
    29672901    LOGICAL, INTENT(INOUT) ::  found   !< flag if output variable is found
    2968     LOGICAL, INTENT(INOUT) ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
    2969     LOGICAL ::  resorted  !< flag if output is already resorted
     2902    LOGICAL, INTENT(INOUT) ::  two_d   !< flag parameter that indicates 2D variables
     2903                                       !<  (horizontal cross sections)
     2904
    29702905
    29712906    REAL(wp), PARAMETER ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
    29722907
    29732908    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do), INTENT(INOUT) ::  local_pf !< local
    2974        !< array to which output data is resorted to
     2909                                                        !< array to which output data is resorted to
    29752910
    29762911    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
     
    29882923             to_be_resorted => nc
    29892924          ELSE
    2990              IF ( .NOT. ALLOCATED( nc_av ) ) THEN
     2925             IF ( .NOT. ALLOCATED( nc_av ) )  THEN
    29912926                ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    29922927                nc_av = REAL( fill_value, KIND = wp )
     
    30002935             to_be_resorted => ni
    30012936          ELSE
    3002              IF ( .NOT. ALLOCATED( ni_av ) ) THEN
     2937             IF ( .NOT. ALLOCATED( ni_av ) )  THEN
    30032938                ALLOCATE( ni_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    30042939                ni_av = REAL( fill_value, KIND = wp )
     
    30122947             to_be_resorted => nr
    30132948          ELSE
    3014              IF ( .NOT. ALLOCATED( nr_av ) ) THEN
     2949             IF ( .NOT. ALLOCATED( nr_av ) )  THEN
    30152950                ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    30162951                nr_av = REAL( fill_value, KIND = wp )
     
    30432978             ENDDO
    30442979          ELSE
    3045              IF ( .NOT. ALLOCATED( prr_av ) ) THEN
     2980             IF ( .NOT. ALLOCATED( prr_av ) )  THEN
    30462981                ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    30472982                prr_av = REAL( fill_value, KIND = wp )
     
    30632998             to_be_resorted => qc
    30642999          ELSE
    3065              IF ( .NOT. ALLOCATED( qc_av ) ) THEN
     3000             IF ( .NOT. ALLOCATED( qc_av ) )  THEN
    30663001                ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    30673002                qc_av = REAL( fill_value, KIND = wp )
     
    30753010             to_be_resorted => qi
    30763011          ELSE
    3077              IF ( .NOT. ALLOCATED( qi_av ) ) THEN
     3012             IF ( .NOT. ALLOCATED( qi_av ) )  THEN
    30783013                ALLOCATE( qi_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    30793014                qi_av = REAL( fill_value, KIND = wp )
     
    30873022             to_be_resorted => ql
    30883023          ELSE
    3089              IF ( .NOT. ALLOCATED( ql_av ) ) THEN
     3024             IF ( .NOT. ALLOCATED( ql_av ) )  THEN
    30903025                ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    30913026                ql_av = REAL( fill_value, KIND = wp )
     
    30993034             to_be_resorted => qr
    31003035          ELSE
    3101              IF ( .NOT. ALLOCATED( qr_av ) ) THEN
     3036             IF ( .NOT. ALLOCATED( qr_av ) )  THEN
    31023037                ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    31033038                qr_av = REAL( fill_value, KIND = wp )
     
    31173052          DO  j = nys, nyn
    31183053             DO  k = nzb_do, nzt_do
    3119                 local_pf(i,j,k) = MERGE( &
    3120                                      to_be_resorted(k,j,i),                    &
    3121                                      REAL( fill_value, KIND = wp ),            &
    3122                                      BTEST( wall_flags_total_0(k,j,i), flag_nr )     &
    3123                                   )
     3054                local_pf(i,j,k) = MERGE(                                                           &
     3055                                         to_be_resorted(k,j,i),                                    &
     3056                                         REAL( fill_value, KIND = wp ),                            &
     3057                                         BTEST( wall_flags_total_0(k,j,i), flag_nr )               &
     3058                                       )
    31243059             ENDDO
    31253060          ENDDO
     
    31303065
    31313066
    3132 !------------------------------------------------------------------------------!
     3067!--------------------------------------------------------------------------------------------------!
    31333068! Description:
    31343069! ------------
    31353070!> Define 3D output variables.
    3136 !------------------------------------------------------------------------------!
     3071!--------------------------------------------------------------------------------------------------!
    31373072 SUBROUTINE bcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
    31383073
     
    31423077    CHARACTER (LEN=*), INTENT(IN) ::  variable   !< name of variable
    31433078
    3144     INTEGER(iwp), INTENT(IN) ::  av        !< flag for (non-)average output
    3145     INTEGER(iwp), INTENT(IN) ::  nzb_do    !< lower limit of the data output (usually 0)
    3146     INTEGER(iwp), INTENT(IN) ::  nzt_do    !< vertical upper limit of the data output (usually nz_do3d)
     3079    INTEGER(iwp), INTENT(IN) ::  av       !< flag for (non-)average output
     3080    INTEGER(iwp), INTENT(IN) ::  nzb_do   !< lower limit of the data output (usually 0)
     3081    INTEGER(iwp), INTENT(IN) ::  nzt_do   !< vertical upper limit of the data output
     3082                                          !< (usually nz_do3d)
    31473083
    31483084    INTEGER(iwp) ::  flag_nr   !< number of masking flag
    3149 
    31503085    INTEGER(iwp) ::  i         !< loop index along x-direction
    31513086    INTEGER(iwp) ::  j         !< loop index along y-direction
    31523087    INTEGER(iwp) ::  k         !< loop index along z-direction
    31533088
     3089    LOGICAL      ::  resorted  !< flag if output is already resorted
     3090
    31543091    LOGICAL, INTENT(INOUT) ::  found     !< flag if output variable is found
    3155     LOGICAL ::  resorted  !< flag if output is already resorted
    31563092
    31573093    REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
    31583094
    31593095    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do), INTENT(INOUT) ::  local_pf   !< local
    3160        !< array to which output data is resorted to
     3096                                                        !< array to which output data is resorted to
    31613097
    31623098    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
     
    31853121             to_be_resorted => ni
    31863122          ELSE
    3187              IF ( .NOT. ALLOCATED( ni_av ) ) THEN
     3123             IF ( .NOT. ALLOCATED( ni_av ) )  THEN
    31883124                ALLOCATE( ni_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    31893125                ni_av = REAL( fill_value, KIND = wp )
     
    31963132             to_be_resorted => nr
    31973133          ELSE
    3198              IF ( .NOT. ALLOCATED( nr_av ) ) THEN
     3134             IF ( .NOT. ALLOCATED( nr_av ) )  THEN
    31993135                ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    32003136                nr_av = REAL( fill_value, KIND = wp )
     
    32133149             ENDDO
    32143150          ELSE
    3215              IF ( .NOT. ALLOCATED( prr_av ) ) THEN
     3151             IF ( .NOT. ALLOCATED( prr_av ) )  THEN
    32163152                ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    32173153                prr_av = REAL( fill_value, KIND = wp )
     
    32313167             to_be_resorted => qc
    32323168          ELSE
    3233              IF ( .NOT. ALLOCATED( qc_av ) ) THEN
     3169             IF ( .NOT. ALLOCATED( qc_av ) )  THEN
    32343170                ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    32353171                qc_av = REAL( fill_value, KIND = wp )
     
    32423178             to_be_resorted => qi
    32433179          ELSE
    3244              IF ( .NOT. ALLOCATED( qi_av ) ) THEN
     3180             IF ( .NOT. ALLOCATED( qi_av ) )  THEN
    32453181                ALLOCATE( qi_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    32463182                qi_av = REAL( fill_value, KIND = wp )
     
    32533189             to_be_resorted => ql
    32543190          ELSE
    3255              IF ( .NOT. ALLOCATED( ql_av ) ) THEN
     3191             IF ( .NOT. ALLOCATED( ql_av ) )  THEN
    32563192                ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    32573193                ql_av = REAL( fill_value, KIND = wp )
     
    32643200             to_be_resorted => qr
    32653201          ELSE
    3266              IF ( .NOT. ALLOCATED( qr_av ) ) THEN
     3202             IF ( .NOT. ALLOCATED( qr_av ) )  THEN
    32673203                ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    32683204                qr_av = REAL( fill_value, KIND = wp )
     
    32813217          DO  j = nys, nyn
    32823218             DO  k = nzb_do, nzt_do
    3283                 local_pf(i,j,k) = MERGE(                                       &
    3284                                      to_be_resorted(k,j,i),                    &
    3285                                      REAL( fill_value, KIND = wp ),            &
    3286                                      BTEST( wall_flags_total_0(k,j,i), flag_nr )     &
    3287                                   )
     3219                local_pf(i,j,k) = MERGE(                                                           &
     3220                                     to_be_resorted(k,j,i),                                        &
     3221                                     REAL( fill_value, KIND = wp ),                                &
     3222                                     BTEST( wall_flags_total_0(k,j,i), flag_nr )                   &
     3223                                       )
    32883224             ENDDO
    32893225          ENDDO
     
    32943230
    32953231
    3296 !------------------------------------------------------------------------------!
     3232!--------------------------------------------------------------------------------------------------!
    32973233! Description:
    32983234! ------------
    32993235!> Read module-specific global restart data (Fortran binary format).
    3300 !------------------------------------------------------------------------------!
     3236!--------------------------------------------------------------------------------------------------!
    33013237    SUBROUTINE bcm_rrd_global_ftn( found )
    33023238
    33033239
    3304        USE control_parameters,                                                 &
     3240       USE control_parameters,                                                                     &
    33053241           ONLY: length, restart_string
    33063242
     
    33793315
    33803316
    3381 !------------------------------------------------------------------------------!
     3317!--------------------------------------------------------------------------------------------------!
    33823318! Description:
    33833319! ------------
    33843320!> Read module-specific global restart data (MPI-IO).
    3385 !------------------------------------------------------------------------------!
     3321!--------------------------------------------------------------------------------------------------!
    33863322    SUBROUTINE bcm_rrd_global_mpi
    33873323
     
    34103346
    34113347
    3412 !------------------------------------------------------------------------------!
     3348!--------------------------------------------------------------------------------------------------!
    34133349! Description:
    34143350! ------------
    34153351!> Read module-specific local restart data arrays (Fortran binary format).
    3416 !------------------------------------------------------------------------------!
    3417     SUBROUTINE bcm_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,          &
    3418                                   nxr_on_file, nynf, nync, nyn_on_file, nysf,      &
    3419                                   nysc, nys_on_file, tmp_2d, tmp_3d, found )
     3352!--------------------------------------------------------------------------------------------------!
     3353    SUBROUTINE bcm_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync, &
     3354                                  nyn_on_file, nysf, nysc, nys_on_file, tmp_2d, tmp_3d, found )
    34203355
    34213356
     
    34463381
    34473382       REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d   !<
     3383
    34483384       REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
    34493385
     
    34553391          CASE ( 'nc' )
    34563392             IF ( k == 1 )  READ ( 13 )  tmp_3d
    3457              nc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
     3393             nc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
    34583394                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    34593395
     
    34633399             ENDIF
    34643400             IF ( k == 1 )  READ ( 13 )  tmp_3d
    3465              nc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
     3401             nc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                    &
    34663402                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    34673403
    34683404          CASE ( 'ni' )
    34693405             IF ( k == 1 )  READ ( 13 )  tmp_3d
    3470              ni(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
     3406             ni(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
    34713407                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    34723408
     
    34763412             ENDIF
    34773413             IF ( k == 1 )  READ ( 13 )  tmp_3d
    3478              ni_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
     3414             ni_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                    &
    34793415                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    34803416
    34813417          CASE ( 'nr' )
    34823418             IF ( k == 1 )  READ ( 13 )  tmp_3d
    3483              nr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
     3419             nr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
    34843420                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    34853421
     
    34893425             ENDIF
    34903426             IF ( k == 1 )  READ ( 13 )  tmp_3d
    3491              nr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
     3427             nr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                    &
    34923428                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    34933429
    34943430          CASE ( 'precipitation_amount' )
    34953431             IF ( k == 1 )  READ ( 13 )  tmp_2d
    3496              precipitation_amount(nysc-nbgp:nync+nbgp,                   &
    3497                                   nxlc-nbgp:nxrc+nbgp)  =                &
    3498                    tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3432             precipitation_amount(nysc-nbgp:nync+nbgp, nxlc-nbgp:nxrc+nbgp) =                      &
     3433                tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    34993434
    35003435          CASE ( 'prr' )
     
    35033438             ENDIF
    35043439             IF ( k == 1 )  READ ( 13 )  tmp_3d
    3505              prr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =            &
     3440             prr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                      &
    35063441                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    35073442
     
    35113446             ENDIF
    35123447             IF ( k == 1 )  READ ( 13 )  tmp_3d
    3513              prr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
     3448             prr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                   &
    35143449                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    35153450
    35163451          CASE ( 'qc' )
    35173452             IF ( k == 1 )  READ ( 13 )  tmp_3d
    3518              qc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
     3453             qc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
    35193454                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    35203455
    35213456          CASE ( 'qi' )
    35223457             IF ( k == 1 )  READ ( 13 )  tmp_3d
    3523              qi(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
     3458             qi(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
    35243459                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    35253460
     
    35293464             ENDIF
    35303465             IF ( k == 1 )  READ ( 13 )  tmp_3d
    3531              qi_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
     3466             qi_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                    &
    35323467                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    35333468
     
    35373472             ENDIF
    35383473             IF ( k == 1 )  READ ( 13 )  tmp_3d
    3539              qc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
     3474             qc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                    &
    35403475                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    35413476
    35423477          CASE ( 'ql' )
    35433478             IF ( k == 1 )  READ ( 13 )  tmp_3d
    3544              ql(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
     3479             ql(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
    35453480                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    35463481
     
    35503485             ENDIF
    35513486             IF ( k == 1 )  READ ( 13 )  tmp_3d
    3552              ql_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
     3487             ql_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                    &
    35533488                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    35543489
    35553490          CASE ( 'qr' )
    35563491             IF ( k == 1 )  READ ( 13 )  tmp_3d
    3557              qr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
     3492             qr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
    35583493                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    35593494
     
    35633498             ENDIF
    35643499             IF ( k == 1 )  READ ( 13 )  tmp_3d
    3565              qr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
     3500             qr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                    &
    35663501                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    35673502!
     
    35763511
    35773512
    3578 !------------------------------------------------------------------------------!
     3513!--------------------------------------------------------------------------------------------------!
    35793514! Description:
    35803515! ------------
    35813516!> Read module-specific local restart data arrays (MPI-IO).
    3582 !------------------------------------------------------------------------------!
     3517!--------------------------------------------------------------------------------------------------!
    35833518    SUBROUTINE bcm_rrd_local_mpi
    35843519
     
    36783613
    36793614
    3680 !------------------------------------------------------------------------------!
     3615!--------------------------------------------------------------------------------------------------!
    36813616! Description:
    36823617! ------------
    36833618!> This routine writes the respective restart data for the bulk cloud module.
    3684 !------------------------------------------------------------------------------!
     3619!--------------------------------------------------------------------------------------------------!
    36853620    SUBROUTINE bcm_wrd_global
    36863621
     
    37453680          WRITE ( 14 )  ice_crystal_sedimentation
    37463681
    3747        ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
     3682       ELSEIF ( TRIM( restart_data_format_output(1:3) ) == 'mpi' )  THEN
    37483683
    37493684          CALL wrd_mpi_io( 'c_sedimentation', c_sedimentation )
     
    37713706
    37723707
    3773 !------------------------------------------------------------------------------!
     3708!--------------------------------------------------------------------------------------------------!
    37743709! Description:
    37753710! ------------
    37763711!> This routine writes the respective restart data for the bulk cloud module.
    3777 !------------------------------------------------------------------------------!
     3712!--------------------------------------------------------------------------------------------------!
    37783713    SUBROUTINE bcm_wrd_local
    37793714
     
    38683803          ENDIF
    38693804
    3870        ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
     3805       ELSEIF ( TRIM( restart_data_format_output(1:3) ) == 'mpi' )  THEN
    38713806
    38723807          IF ( ALLOCATED( prr ) )  CALL wrd_mpi_io( 'prr', prr )
     
    39003835    END SUBROUTINE bcm_wrd_local
    39013836
    3902 !------------------------------------------------------------------------------!
     3837!--------------------------------------------------------------------------------------------------!
    39033838! Description:
    39043839! ------------
    3905 !> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
    3906 !> evaporation of rain drops due to too small or too big weights
    3907 !> of rain drops (Stevens and Seifert, 2008).
    3908 !------------------------------------------------------------------------------!
     3840!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and evaporation of rain
     3841!> drops due to too small or too big weights of rain drops (Stevens and Seifert, 2008).
     3842!--------------------------------------------------------------------------------------------------!
    39093843    SUBROUTINE adjust_cloud
    39103844
     
    39383872                ENDIF
    39393873
    3940                 IF ( microphysics_morrison ) THEN
     3874                IF ( microphysics_morrison )  THEN
    39413875                   IF ( qc(k,j,i) <= eps_sb )  THEN
    39423876                      qc(k,j,i) = 0.0_wp
     
    39573891    END SUBROUTINE adjust_cloud
    39583892
    3959 !------------------------------------------------------------------------------!
     3893!--------------------------------------------------------------------------------------------------!
    39603894! Description:
    39613895! ------------
    3962 !> Adjust number of raindrops to avoid nonlinear effects in
    3963 !> sedimentation and evaporation of rain drops due to too small or
    3964 !> too big weights of rain drops (Stevens and Seifert, 2008).
    3965 !> The same procedure is applied to cloud droplets if they are determined
    3966 !> prognostically. Call for grid point i,j
    3967 !------------------------------------------------------------------------------!
     3896!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and evaporation of rain
     3897!> drops due to too small or too big weights of rain drops (Stevens and Seifert, 2008).
     3898!> The same procedure is applied to cloud droplets if they are determined prognostically. Call for
     3899!> grid point i,j
     3900!--------------------------------------------------------------------------------------------------!
    39683901    SUBROUTINE adjust_cloud_ij( i, j )
    39693902
     
    39813914          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    39823915
    3983           IF ( .NOT.  microphysics_morrison_no_rain )  THEN
     3916          IF ( .NOT. microphysics_morrison_no_rain )  THEN
    39843917             IF ( qr(k,j,i) <= eps_sb )  THEN
    39853918                qr(k,j,i) = 0.0_wp
     
    39873920             ELSE
    39883921!
    3989 !--             Adjust number of raindrops to avoid nonlinear effects in
    3990 !--             sedimentation and evaporation of rain drops due to too small or
    3991 !--             too big weights of rain drops (Stevens and Seifert, 2008).
     3922!--             Adjust number of raindrops to avoid nonlinear effects in sedimentation and
     3923!--             evaporation of rain drops due to too small or too big weights of rain drops
     3924!--            (Stevens and Seifert, 2008).
    39923925                IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
    39933926                   nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
     
    39983931          ENDIF
    39993932
    4000           IF ( microphysics_morrison ) THEN
     3933          IF ( microphysics_morrison )  THEN
    40013934             IF ( qc(k,j,i) <= eps_sb )  THEN
    40023935                qc(k,j,i) = 0.0_wp
     
    40133946    END SUBROUTINE adjust_cloud_ij
    40143947
    4015 !------------------------------------------------------------------------------!
     3948!--------------------------------------------------------------------------------------------------!
    40163949! Description:
    40173950! ------------
    4018 !> Adjust number of ice crystal to avoid nonlinear effects in sedimentation and
    4019 !> evaporation of ice crytals due to too small or too big weights
    4020 !> of ice crytals (Stevens and Seifert, 2008).
    4021 !------------------------------------------------------------------------------!
     3951!> Adjust number of ice crystal to avoid nonlinear effects in sedimentation and evaporation of ice
     3952!> crytals due to too small or too big weights of ice crytals (Stevens and Seifert, 2008).
     3953!--------------------------------------------------------------------------------------------------!
    40223954    SUBROUTINE adjust_ice
    40233955
     
    40543986    END SUBROUTINE adjust_ice
    40553987
    4056 !------------------------------------------------------------------------------!
     3988!--------------------------------------------------------------------------------------------------!
    40573989! Description:
    40583990! ------------
    4059 !> Adjust number of of ice crystal to avoid nonlinear effects in
    4060 !> sedimentation and evaporation of ice crystals due to too small or
    4061 !> too big weights of ice crytals (Stevens and Seifert, 2008).
    4062 !> The same procedure is applied to cloud droplets if they are determined
    4063 !> prognostically. Call for grid point i,j
    4064 !------------------------------------------------------------------------------!
     3991!> Adjust number of of ice crystal to avoid nonlinear effects in sedimentation and evaporation of
     3992!> ice crystals due to too small or too big weights of ice crytals (Stevens and Seifert, 2008).
     3993!> The same procedure is applied to cloud droplets if they are determined prognostically. Call for
     3994!> grid point i,j
     3995!--------------------------------------------------------------------------------------------------!
    40653996    SUBROUTINE adjust_ice_ij( i, j )
    40663997
     
    40904021
    40914022
    4092 !------------------------------------------------------------------------------!
     4023!--------------------------------------------------------------------------------------------------!
    40934024! Description:
    40944025! ------------
    40954026!> Calculate number of activated condensation nucleii after simple activation
    40964027!> scheme of Twomey, 1959.
    4097 !------------------------------------------------------------------------------!
     4028!--------------------------------------------------------------------------------------------------!
    40984029    SUBROUTINE activation
    40994030
     
    41084039       REAL(wp)     ::  beta_act          !<
    41094040       REAL(wp)     ::  bfactor           !<
     4041       REAL(wp)     ::  flag              !< flag to indicate first grid level above
    41104042       REAL(wp)     ::  k_act             !<
    41114043       REAL(wp)     ::  n_act             !<
     
    41164048       REAL(wp)     ::  sigma_act         !<
    41174049
    4118        REAL(wp) ::  flag                  !< flag to indicate first grid level above
    41194050
    41204051       CALL cpu_log( log_point_s(65), 'activation', 'start' )
     
    41314062                CALL supersaturation ( i, j, k )
    41324063!
    4133 !--             Prescribe parameters for activation
    4134 !--             (see: Bott + Trautmann, 2002, Atm. Res., 64)
     4064!--             Prescribe parameters for activation (see: Bott + Trautmann, 2002, Atm. Res., 64)
    41354065                k_act  = 0.7_wp
    41364066                activ  = 0.0_wp
    41374067
    41384068
    4139                 IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
     4069                IF ( sat > 0.0  .AND.  .NOT. curvature_solution_effects_bulk ) THEN
    41404070!
    41414071!--                Compute the number of activated Aerosols
     
    41434073                   n_act     = na_init * sat**k_act
    41444074!
    4145 !--                Compute the number of cloud droplets
    4146 !--                (see: Morrison + Grabowski, 2007, JAS, 64)
     4075!--                Compute the number of cloud droplets (see: Morrison + Grabowski, 2007, JAS, 64)
    41474076!                  activ = MAX( n_act - nc(k,j,i), 0.0_wp) / dt_micro
    41484077
     
    41514080!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
    41524081                   sat_max = 1.0_wp / 100.0_wp
    4153                    activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN       &
    4154                       ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /    &
    4155                        dt_micro
    4156                 ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk ) THEN
    4157 !
    4158 !--                Curvature effect (afactor) with surface tension
    4159 !--                parameterization by Straka (2009)
     4082                   activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) *                               &
     4083                             MIN( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) / dt_micro
     4084                ELSEIF ( sat > 0.0  .AND.  curvature_solution_effects_bulk )  THEN
     4085!
     4086!--                Curvature effect (afactor) with surface tension parameterization by Straka (2009)
    41604087                   sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
    41614088                   afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
    41624089!
    41634090!--                Solute effect (bfactor)
    4164                    bfactor = vanthoff * molecular_weight_of_water *            &
    4165                        rho_s / ( molecular_weight_of_solute * rho_l )
    4166 
    4167 !
    4168 !--                Prescribe power index that describes the soluble fraction
    4169 !--                of an aerosol particle (beta)
    4170 !--                (see: Morrison + Grabowski, 2007, JAS, 64)
     4091                   bfactor = vanthoff * molecular_weight_of_water * rho_s /                        &
     4092                             ( molecular_weight_of_solute * rho_l )
     4093
     4094!
     4095!--                Prescribe power index that describes the soluble fraction of an aerosol particle
     4096!--                (beta) (see: Morrison + Grabowski, 2007, JAS, 64)
    41714097                   beta_act  = 0.5_wp
    41724098                   sigma_act = sigma_bulk**( 1.0_wp + beta_act )
    41734099!
    4174 !--                Calculate mean geometric supersaturation (s_0) with
    4175 !--                parameterization by Khvorostyanov and Curry (2006)
    4176                    s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *    &
    4177                        ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
    4178 
    4179 !
    4180 !--                Calculate number of activated CCN as a function of
    4181 !--                supersaturation and taking Koehler theory into account
    4182 !--                (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
    4183                    n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(              &
    4184                       LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
     4100!--                Calculate mean geometric supersaturation (s_0) with parameterization by
     4101!--                Khvorostyanov and Curry (2006)
     4102                   s_0 = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) * ( 4.0_wp * afactor**3 / &
     4103                         ( 27.0_wp * bfactor ) )**0.5_wp
     4104
     4105!
     4106!--                Calculate number of activated CCN as a function of supersaturation and taking
     4107!--                Koehler theory into account (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
     4108                   n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                                  &
     4109                           LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
    41854110                   activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp )
    41864111                ENDIF
     
    41964121    END SUBROUTINE activation
    41974122
    4198 !------------------------------------------------------------------------------!
     4123!--------------------------------------------------------------------------------------------------!
    41994124! Description:
    42004125! ------------
    4201 !> Calculate number of activated condensation nucleii after simple activation
    4202 !> scheme of Twomey, 1959.
    4203 !------------------------------------------------------------------------------!
     4126!> Calculate number of activated condensation nucleii after simple activation scheme of
     4127!> Twomey, 1959.
     4128!--------------------------------------------------------------------------------------------------!
    42044129    SUBROUTINE activation_ij( i, j )
    42054130
     
    42314156          CALL supersaturation ( i, j, k )
    42324157!
    4233 !--       Prescribe parameters for activation
    4234 !--       (see: Bott + Trautmann, 2002, Atm. Res., 64)
     4158!--       Prescribe parameters for activation (see: Bott + Trautmann, 2002, Atm. Res., 64)
    42354159          k_act  = 0.7_wp
    42364160          activ  = 0.0_wp
     
    42504174!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
    42514175             sat_max = 0.8_wp / 100.0_wp
    4252              activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN             &
    4253                  ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /         &
    4254                   dt_micro
     4176             activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) *                                     &
     4177                       MIN( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) / dt_micro
    42554178
    42564179             nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)
    4257           ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk )  THEN
    4258 !
    4259 !--          Curvature effect (afactor) with surface tension
    4260 !--          parameterization by Straka (2009)
     4180          ELSEIF ( sat > 0.0  .AND.  curvature_solution_effects_bulk )  THEN
     4181!
     4182!--          Curvature effect (afactor) with surface tension parameterization by Straka (2009)
    42614183             sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
    42624184             afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
    42634185!
    42644186!--          Solute effect (bfactor)
    4265              bfactor = vanthoff * molecular_weight_of_water *                  &
    4266                   rho_s / ( molecular_weight_of_solute * rho_l )
    4267 
    4268 !
    4269 !--          Prescribe power index that describes the soluble fraction
    4270 !--          of an aerosol particle (beta).
    4271 !--          (see: Morrison + Grabowski, 2007, JAS, 64)
     4187             bfactor = vanthoff * molecular_weight_of_water * rho_s /                              &
     4188                       ( molecular_weight_of_solute * rho_l )
     4189
     4190!
     4191!--          Prescribe power index that describes the soluble fraction of an aerosol particle
     4192!--          (beta). (see: Morrison + Grabowski, 2007, JAS, 64)
    42724193             beta_act  = 0.5_wp
    42734194             sigma_act = sigma_bulk**( 1.0_wp + beta_act )
    42744195!
    4275 !--          Calculate mean geometric supersaturation (s_0) with
    4276 !--          parameterization by Khvorostyanov and Curry (2006)
    4277              s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *         &
    4278                ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
    4279 
    4280 !
    4281 !--          Calculate number of activated CCN as a function of
    4282 !--          supersaturation and taking Koehler theory into account
    4283 !--          (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
    4284              n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                    &
    4285                 LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
     4196!--          Calculate mean geometric supersaturation (s_0) with parameterization by Khvorostyanov
     4197!--          and Curry (2006)
     4198             s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *                             &
     4199                     ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
     4200
     4201!
     4202!--          Calculate number of activated CCN as a function of supersaturation and taking Koehler
     4203!--          theory into account (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
     4204             n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                                        &
     4205                     LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
    42864206             activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
    42874207
     
    42934213    END SUBROUTINE activation_ij
    42944214
    4295 !------------------------------------------------------------------------------!
     4215!--------------------------------------------------------------------------------------------------!
    42964216! Description:
    42974217! ------------
    4298 !> Calculate ice nucleation by applying the deposition-condensation formula as
    4299 !> given by Meyers et al 1992 and as described in Seifert and Beheng 2006
    4300 !------------------------------------------------------------------------------!
     4218!> Calculate ice nucleation by applying the deposition-condensation formula as given by
     4219!> Meyers et al 1992 and as described in Seifert and Beheng 2006
     4220!--------------------------------------------------------------------------------------------------!
    43014221    SUBROUTINE ice_nucleation
    43024222
     
    43264246                   CALL supersaturation_ice ( i, j, k )
    43274247                   nucle = 0.0_wp
    4328                    IF ( sat_ice >= 0.05_wp  .OR.  ql(k,j,i) >= 0.001E-3_wp  ) THEN
     4248                   IF ( sat_ice >= 0.05_wp  .OR.  ql(k,j,i) >= 0.001E-3_wp  )  THEN
    43294249!
    43304250!--                   Calculate ice nucleation
     
    43644284
    43654285
    4366 !------------------------------------------------------------------------------!
     4286!--------------------------------------------------------------------------------------------------!
    43674287! Description:
    43684288! ------------
    4369 !> Calculate ice nucleation by applying the deposition-condensation formula as
    4370 !> given by Meyers et al 1992 and as described in Seifert and Beheng 2006
    4371 !------------------------------------------------------------------------------!
     4289!> Calculate ice nucleation by applying the deposition-condensation formula as given by
     4290!> Meyers et al 1992 and as described in Seifert and Beheng 2006.
     4291!--------------------------------------------------------------------------------------------------!
    43724292    SUBROUTINE ice_nucleation_ij( i, j )
    43734293
     
    43934313             CALL supersaturation_ice ( i, j, k )
    43944314             nucle = 0.0_wp
    4395              IF ( sat_ice >= 0.05_wp  .OR.  ql(k,j,i) >= 0.001E-3_wp ) THEN
     4315             IF ( sat_ice >= 0.05_wp  .OR.  ql(k,j,i) >= 0.001E-3_wp )  THEN
    43964316!
    43974317!--             Calculate ice nucleation
     
    44094329             CALL supersaturation_ice ( i, j, k )
    44104330             nucle = 0.0_wp
    4411              IF ( sat_ice > 0.0 ) THEN
     4331             IF ( sat_ice > 0.0 )  THEN
    44124332!
    44134333!--             Calculate ice nucleation
     
    44224342    END SUBROUTINE ice_nucleation_ij
    44234343
    4424 !------------------------------------------------------------------------------!
     4344!--------------------------------------------------------------------------------------------------!
    44254345! Description:
    44264346! ------------
    4427 !> Calculate condensation rate for cloud water content (after Khairoutdinov and
    4428 !> Kogan, 2000).
    4429 !------------------------------------------------------------------------------!
     4347!> Calculate condensation rate for cloud water content (after Khairoutdinov and Kogan, 2000).
     4348!--------------------------------------------------------------------------------------------------!
    44304349    SUBROUTINE condensation
    44314350
     
    44554374                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    44564375!
    4457 !--             Call calculation of supersaturation 
     4376!--             Call calculation of supersaturation
    44584377                CALL supersaturation ( i, j, k )
    44594378!
    4460 !--             Actual temperature, t_l is calculated directly before
    4461 !--             in supersaturation
     4379!--             Actual temperature, t_l is calculated directly before in supersaturation
    44624380                IF ( microphysics_ice_phase ) THEN
    44634381                   temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i)
     
    44664384                ENDIF
    44674385
    4468                 g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
    4469                                     l_v / ( thermal_conductivity_l * temp )    &
    4470                                     + r_v * temp / ( diff_coeff_l * e_s )      &
    4471                                     )
     4386                g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *                            &
     4387                                      l_v / ( thermal_conductivity_l  * temp )                     &
     4388                                    + r_v * temp / ( diff_coeff_l * e_s )                          &
     4389                                  )
    44724390!
    44734391!--             Mean weight of cloud drops
     
    44754393!
    44764394!--             Calculating mean radius of cloud droplets assuming gamma distribution with shape
    4477 !--             parameter nu=1 (Seifert and Beheng, 2006). Tuning factor alpha_rc (intorduced
    4478 !--             in Seifert and Stevens, 2010 ) is switched off. Minimum radius is set to 1µm
     4395!--             parameter nu=1 (Seifert and Beheng, 2006). Tuning factor alpha_rc (introduced
     4396!--             in Seifert and Stevens, 2010 ) is switched off. Minimum radius is set to 1µm
    44794397!--             following (Seifert and Beheng, 2006, Kogan and Khairoutdinov, 2000,
    44804398!--             Seifert and Stevens, 2010)
    4481                 rc = MAX( alpha_rc  * GAMMA(nu + 1.33_wp) / GAMMA(nu + 1.0_wp) *                   &
    4482                         ( 3.0_wp * qc(k,j,i) /                                                     &
     4399                rc = MAX( alpha_rc  * GAMMA( nu + 1.33_wp ) / GAMMA( nu + 1.0_wp ) *               &
     4400                          ( 3.0_wp * qc(k,j,i) /                                                   &
    44834401                                   ( 4.0_wp * pi * rho_l  * ( nu + 2.0_wp ) * nc(k,j,i) )          &
    4484                         )**0.33_wp, 1.0E-6_wp )
     4402                          )**0.33_wp, 1.0E-6_wp )
    44854403!
    44864404!--             Condensation needs only to be calculated in supersaturated regions
     
    44914409!--                2007, and Seifert and Stevens, 2010)
    44924410                   cond      = 4.0_wp * pi * nc(k,j,i) * g_fac * sat * rc / hyrho(k)
    4493                    IF ( microphysics_seifert ) THEN
     4411                   IF ( microphysics_seifert )  THEN
    44944412                      cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
    4495                    ELSEIF ( microphysics_morrison_no_rain ) THEN
    4496                       cond_max  = q(k,j,i) - q_s - qc(k,j,i) 
     4413                   ELSEIF ( microphysics_morrison_no_rain )  THEN
     4414                      cond_max  = q(k,j,i) - q_s - qc(k,j,i)
    44974415                   ENDIF
    44984416                   cond      = MIN( cond, cond_max / dt_micro )
    44994417
    45004418                   qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag
    4501                 ELSEIF ( sat < 0.0_wp ) THEN
     4419                ELSEIF ( sat < 0.0_wp )  THEN
    45024420                   evap      = 4.0_wp * pi * nc(k,j,i) * g_fac * sat * rc / hyrho(k)
    45034421                   evap      = MAX( evap, -qc(k,j,i) / dt_micro )
     
    45164434    END SUBROUTINE condensation
    45174435
    4518 !------------------------------------------------------------------------------!
     4436!--------------------------------------------------------------------------------------------------!
    45194437! Description:
    45204438! ------------
    4521 !> Calculate condensation rate for cloud water content (after Khairoutdinov and
    4522 !> Kogan, 2000).
    4523 !------------------------------------------------------------------------------!
     4439!> Calculate condensation rate for cloud water content (after Khairoutdinov and Kogan, 2000).
     4440!--------------------------------------------------------------------------------------------------!
    45244441    SUBROUTINE condensation_ij( i, j )
    45254442
     
    45484465          CALL supersaturation ( i, j, k )
    45494466!
    4550 !--       Actual temperature, t_l is calculated directly before
    4551 !--       in supersaturation
     4467!--       Actual temperature, t_l is calculated directly before in supersaturation
    45524468          IF ( microphysics_ice_phase ) THEN
    45534469             temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i)
     
    45564472          ENDIF
    45574473
    4558           g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *              &
    4559                               l_v / ( thermal_conductivity_l * temp )          &
    4560                               + r_v * temp / ( diff_coeff_l * e_s )            &
     4474          g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *                                  &
     4475                                l_v / ( thermal_conductivity_l  * temp )                           &
     4476                              + r_v * temp / ( diff_coeff_l * e_s )                                &
    45614477                            )
    45624478!
    45634479!--       Mean weight of cloud drops
    4564           IF ( nc(k,j,i) <= 0.0_wp) CYCLE
     4480          IF ( nc(k,j,i) <= 0.0_wp)  CYCLE
    45654481!
    45664482!--       Calculating mean radius of cloud droplets assuming gamma distribution with shape
    4567 !--       parameter nu=1 (Seifert and Beheng, 2006). Tuning factor alpha_rc (intorduced
    4568 !--       in Seifert and Stevens, 2010 ) is switched off. Minimum radius is set to 1µm following
     4483!--       parameter nu=1 (Seifert and Beheng, 2006). Tuning factor alpha_rc (introduced
     4484!--       in Seifert and Stevens, 2010 ) is switched off. Minimum radius is set to 1µm following
    45694485!--       (Seifert and Beheng, 2006, Kogan and Khairoutdinov, 2000, Seifert and Stevens, 2010)
    4570           rc = MAX( alpha_rc  * GAMMA(nu + 1.33_wp) / GAMMA(nu + 1.0_wp) *                         &
    4571                ( 3.0_wp * qc(k,j,i) / ( 4.0_wp * pi * rho_l  * ( nu + 2.0_wp ) * nc(k,j,i) )       &
    4572                )**0.33_wp, 1.0E-6_wp )
     4486          rc = MAX( alpha_rc  * GAMMA( nu + 1.33_wp ) / GAMMA( nu + 1.0_wp ) *                     &
     4487                    ( 3.0_wp * qc(k,j,i) / ( 4.0_wp * pi * rho_l  * ( nu + 2.0_wp ) * nc(k,j,i) )  &
     4488                  )**0.33_wp, 1.0E-6_wp )
    45734489!
    45744490!--       Condensation needs only to be calculated in supersaturated regions
     
    45794495!--          and Seifert and Stevens, 2010)
    45804496             cond      = 4.0_wp * pi * nc(k,j,i) * g_fac * sat * rc / hyrho(k)
    4581              IF ( microphysics_seifert ) THEN
     4497             IF ( microphysics_seifert )  THEN
    45824498                cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
    4583              ELSEIF ( microphysics_morrison_no_rain ) THEN
    4584                 cond_max  = q(k,j,i) - q_s - qc(k,j,i) 
     4499             ELSEIF ( microphysics_morrison_no_rain )  THEN
     4500                cond_max  = q(k,j,i) - q_s - qc(k,j,i)
    45854501             ENDIF
    45864502             cond      = MIN( cond, cond_max / dt_micro )
    45874503
    45884504             qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag
    4589           ELSEIF ( sat < 0.0_wp ) THEN
     4505          ELSEIF ( sat < 0.0_wp )  THEN
    45904506             evap      = 4.0_wp * pi * nc(k,j,i) * g_fac * sat * rc / hyrho(k)
    45914507             evap      = MAX( evap, -qc(k,j,i) / dt_micro )
     
    45974513    END SUBROUTINE condensation_ij
    45984514
    4599 !------------------------------------------------------------------------------!
     4515!--------------------------------------------------------------------------------------------------!
    46004516! Description:
    46014517! ------------
    4602 !> Calculate the growth of ice particles by water vapor deposition (after
    4603 !> Seifert and Beheng, 2006).
    4604 !------------------------------------------------------------------------------!
     4518!> Calculate the growth of ice particles by water vapor deposition (after Seifert and Beheng, 2006).
     4519!--------------------------------------------------------------------------------------------------!
    46054520    SUBROUTINE ice_deposition
    46064521
     
    46114526       REAL(wp) ::  ac = 0.09_wp        !< parameter for ice capacitance
    46124527       REAL(wp) ::  bc = 0.33_wp        !< parameter for ice capacitance
    4613        REAL(wp) ::  fac_gamma = 0.76_wp !< parameter to describe spectral
    4614                                         !< distribution, here following gamma
    4615                                         !< size distribution with µ =1/3 and nu=0
    46164528       REAL(wp) ::  deposition_rate     !< depositions rate
    46174529       REAL(wp) ::  deposition_rate_max !< maximum deposition rate
     4530       REAL(wp) ::  fac_gamma = 0.76_wp !< parameter to describe spectral distribution, here
     4531                                        !< following gamma size distribution with µ =1/3 and nu=0
     4532       REAL(wp) ::  flag                !< flag to indicate first grid level above
     4533       REAL(wp) ::  gfac_dep            !< factor
    46184534       REAL(wp) ::  sublimation_rate    !< sublimations rate
    4619        REAL(wp) ::  gfac_dep            !< factor
    46204535       REAL(wp) ::  temp                !< actual temperature
    46214536       REAL(wp) ::  xi                  !< mean mass of ice crystal
    4622        REAL(wp) ::  flag                !< flag to indicate first grid level above
    46234537
    46244538       CALL cpu_log( log_point_s(95), 'ice deposition', 'start' )
     
    46374551                temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i)
    46384552
    4639                 IF ( temp >= 273.15_wp ) CYCLE
    4640 !
    4641 !--             calculating gfac_dep ( 1/ (Fk + Fd) ) see e.g.
    4642 !--             Rogers and Yau, 1989
    4643                 gfac_dep  = 1.0_wp / ( ( l_s / ( r_v * temp ) - 1.0_wp ) *     &
    4644                                     l_s / ( thermal_conductivity_l * temp )    &
    4645                                     + r_v * temp / ( diff_coeff_l * e_si )     &
    4646                                     )
    4647 !
    4648 !--             If there is nothing nucleated, than there is also no
    4649 !--             deposition (above -38°C)
    4650                 IF ( ni(k,j,i) <= 0.0_wp ) CYCLE
     4553                IF ( temp >= 273.15_wp )  CYCLE
     4554!
     4555!--             calculating gfac_dep ( 1/ (Fk + Fd) ) see e.g. Rogers and Yau, 1989
     4556                gfac_dep  = 1.0_wp / ( ( l_s / ( r_v * temp ) - 1.0_wp ) *                         &
     4557                                         l_s / ( thermal_conductivity_l  * temp )                  &
     4558                                    + r_v * temp / ( diff_coeff_l * e_si )                         &
     4559                                     )
     4560!
     4561!--             If there is nothing nucleated, than there is also no deposition (above -38°C)
     4562                IF ( ni(k,j,i) <= 0.0_wp )  CYCLE
    46514563!
    46524564!--             calculate mean mass of ice crystal
     
    46544566                xi = MAX( ( qi(k,j,i) * hyrho(k) / ni(k,j,i)), ximin )
    46554567!
    4656 !--             Condensation needs only to be calculated in supersaturated
    4657 !--             regions (regarding ice)
     4568!--             Condensation needs only to be calculated in supersaturated regions (regarding ice)
    46584569                IF ( sat_ice > 0.0_wp )  THEN
    46594570!
    4660 !--                Calculate deposition rate assuming ice crystal shape as
    4661 !--                prescribed in Ovchinnikov et al., 2014 and a gamma size
    4662 !--                distribution according to Seifert and Beheng with to
    4663 !--                µ =1/3 and nu=0
    4664                    deposition_rate  = 4.0_wp * pi * sat_ice * gfac_dep *       &
    4665                                       fac_gamma * ac * xi**bc * ni(k,j,i)
    4666                    IF ( microphysics_seifert ) THEN
    4667                       deposition_rate_max  = q(k,j,i) -                        &
    4668                                              q_si - qr(k,j,i) - qi(k,j,i)
     4571!--                Calculate deposition rate assuming ice crystal shape as prescribed in
     4572!--                Ovchinnikov et al., 2014 and a gamma size distribution according to
     4573!--                Seifert and Beheng with µ =1/3 and nu=0
     4574                   deposition_rate  = 4.0_wp * pi * sat_ice * gfac_dep * fac_gamma * ac * xi**bc * &
     4575                                      ni(k,j,i)
     4576                   IF ( microphysics_seifert )  THEN
     4577                      deposition_rate_max  = q(k,j,i) - q_si - qr(k,j,i) - qi(k,j,i)
    46694578                   ELSEIF ( microphysics_morrison_no_rain ) THEN
    4670                       deposition_rate_max  = q(k,j,i) -                        &
    4671                                              q_si - qc(k,j,i) - qi(k,j,i)
     4579                      deposition_rate_max  = q(k,j,i) - q_si - qc(k,j,i) - qi(k,j,i)
    46724580                   ENDIF
    4673                    deposition_rate = MIN( deposition_rate,                     &
    4674                                           deposition_rate_max / dt_micro )
    4675 
     4581                   deposition_rate = MIN( deposition_rate, deposition_rate_max / dt_micro )
    46764582                   qi(k,j,i) = qi(k,j,i) + deposition_rate * dt_micro * flag
    4677                 ELSEIF ( sat_ice < 0.0_wp ) THEN
    4678                    sublimation_rate = 4.0_wp * pi * sat_ice * gfac_dep *       &
    4679                                       fac_gamma * ac * xi**bc * ni(k,j,i)
    4680                    sublimation_rate = MAX( sublimation_rate,                   &
    4681                                            -qi(k,j,i) / dt_micro )
     4583                ELSEIF ( sat_ice < 0.0_wp )  THEN
     4584                   sublimation_rate = 4.0_wp * pi * sat_ice * gfac_dep * fac_gamma * ac * xi**bc * &
     4585                                      ni(k,j,i)
     4586                   sublimation_rate = MAX( sublimation_rate, - qi(k,j,i) / dt_micro )
    46824587                   qi(k,j,i) = qi(k,j,i) + sublimation_rate * dt_micro * flag
    46834588                ENDIF
     
    46904595    END SUBROUTINE ice_deposition
    46914596
    4692 !------------------------------------------------------------------------------!
     4597!--------------------------------------------------------------------------------------------------!
    46934598! Description:
    46944599! ------------
    4695 !> Calculate condensation rate for cloud water content (after Khairoutdinov and
    4696 !> Kogan, 2000).
    4697 !------------------------------------------------------------------------------!
     4600!> Calculate condensation rate for cloud water content (after Khairoutdinov and Kogan, 2000).
     4601!--------------------------------------------------------------------------------------------------!
    46984602    SUBROUTINE ice_deposition_ij( i, j )
    46994603
     
    47044608       REAL(wp) ::  ac = 0.09_wp        !< parameter for ice capacitance
    47054609       REAL(wp) ::  bc = 0.33_wp        !< parameter for ice capacitance
    4706        REAL(wp) ::  fac_gamma = 0.76_wp !< parameter to describe spectral
    4707                                         !< distribution, here following gamma
    4708                                         !< size distribution with nu=1, v=0
    47094610       REAL(wp) ::  deposition_rate     !< depositions rate
    47104611       REAL(wp) ::  deposition_rate_max !< maximum deposition rate
     4612       REAL(wp) ::  fac_gamma = 0.76_wp !< parameter to describe spectral distribution, here
     4613                                        !< following gamma size distribution with nu=1, v=0
     4614       REAL(wp) ::  flag                !< flag to indicate first grid level above
     4615       REAL(wp) ::  gfac_dep            !< factor
    47114616       REAL(wp) ::  sublimation_rate    !< sublimations rate
    4712        REAL(wp) ::  gfac_dep            !< factor
    47134617       REAL(wp) ::  temp                !< actual temperature
    47144618       REAL(wp) ::  xi                  !< mean mass of ice crystal
    4715        REAL(wp) ::  flag                !< flag to indicate first grid level above
    47164619
    47174620       DO  k = nzb+1, nzt
     
    47214624!
    47224625!--       Call calculation of supersaturation over a plane ice surface
    4723           CALL supersaturation_ice ( i, j, k )
     4626          CALL supersaturation_ice (i,j,k)
    47244627!
    47254628!--       Actual temperature:
    47264629          temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i)
    47274630
    4728           IF ( temp >= 273.15_wp ) CYCLE
     4631          IF ( temp >= 273.15_wp )  CYCLE
    47294632!
    47304633!--       calculating gfac_dep ( 1/ (Fk + Fd) ) see e.g.
    47314634!--       Rogers and Yau, 1989
    4732           gfac_dep  = 1.0_wp / ( ( l_s / ( r_v * temp ) - 1.0_wp ) *           &
    4733                               l_s / ( thermal_conductivity_l * temp )          &
    4734                               + r_v * temp / ( diff_coeff_l * e_si )           &
    4735                               )
    4736 !
    4737 !--       If there is nothing nucleated, than there is also no
    4738 !--       deposition (above -38°C)
    4739           IF ( ni(k,j,i) <= 0.0_wp ) CYCLE
     4635          gfac_dep  = 1.0_wp / ( ( l_s / ( r_v * temp ) - 1.0_wp ) *                               &
     4636                                   l_s / ( thermal_conductivity_l  * temp )                        &
     4637                                 + r_v * temp / ( diff_coeff_l * e_si )                            &
     4638                               )
     4639!
     4640!--       If there is nothing nucleated, than there is also no deposition (above -38°C)
     4641          IF ( ni(k,j,i) <= 0.0_wp )  CYCLE
    47404642!
    47414643!--       calculate mean mass of ice crystal
     
    47434645          xi = MAX( ( qi(k,j,i) * hyrho(k) / ni(k,j,i)), ximin )
    47444646!
    4745 !--       Condensation needs only to be calculated in supersaturated
    4746 !--       regions (regarding ice)
     4647!--       Condensation needs only to be calculated in supersaturated regions (regarding ice)
    47474648          IF ( sat_ice > 0.0_wp )  THEN
    47484649!
    4749 !--          Calculate deposition rate assuming ice crystal shape as
    4750 !--          prescribed in Ovchinnikov et al., 2014 and a gamma size
    4751 !--          distribution according to Seifert and Beheng with to
    4752 !--          µ =1/3 and nu=0
    4753              deposition_rate  = 4.0_wp * pi * sat_ice * gfac_dep *             &
    4754                                 fac_gamma * ac * xi**bc * ni(k,j,i)
    4755              IF ( microphysics_seifert ) THEN
    4756                 deposition_rate_max  = q(k,j,i) -                              &
    4757                                        q_si - qr(k,j,i) - qi(k,j,i)
    4758              ELSEIF ( microphysics_morrison_no_rain ) THEN
    4759                 deposition_rate_max  = q(k,j,i) -                              &
    4760                                        q_si - qc(k,j,i) - qi(k,j,i)
    4761              ENDIF
    4762              deposition_rate = MIN( deposition_rate,                           &
    4763                                     deposition_rate_max / dt_micro )
     4650!--          Calculate deposition rate assuming ice crystal shape as prescribed in
     4651!--          Ovchinnikov et al., 2014 and a gamma size distribution according to
     4652!--          Seifert and Beheng with µ =1/3 and nu=0
     4653             deposition_rate  = 4.0_wp * pi * sat_ice * gfac_dep * fac_gamma * ac * xi**bc *       &
     4654                                ni(k,j,i)
     4655             IF ( microphysics_seifert )  THEN
     4656                deposition_rate_max  = q(k,j,i) - q_si - qr(k,j,i) - qi(k,j,i)
     4657             ELSEIF ( microphysics_morrison_no_rain )  THEN
     4658                deposition_rate_max  = q(k,j,i) - q_si - qc(k,j,i) - qi(k,j,i)
     4659             ENDIF
     4660             deposition_rate = MIN( deposition_rate, deposition_rate_max / dt_micro )
    47644661
    47654662             qi(k,j,i) = qi(k,j,i) + deposition_rate * dt_micro * flag
    4766           ELSEIF ( sat_ice < 0.0_wp ) THEN
    4767              sublimation_rate = 4.0_wp * pi * sat_ice * gfac_dep *             &
    4768                                 fac_gamma * ac * xi**bc * ni(k,j,i)
    4769              sublimation_rate = MAX( sublimation_rate,                         &
    4770                                      -qi(k,j,i) / dt_micro )
     4663          ELSEIF ( sat_ice < 0.0_wp )  THEN
     4664             sublimation_rate = 4.0_wp * pi * sat_ice * gfac_dep * fac_gamma * ac * xi**bc *       &
     4665                                ni(k,j,i)
     4666             sublimation_rate = MAX( sublimation_rate, - qi(k,j,i) / dt_micro )
    47714667             qi(k,j,i) = qi(k,j,i) + sublimation_rate * dt_micro * flag
    47724668          ENDIF
     
    47764672
    47774673
    4778 !------------------------------------------------------------------------------!
     4674!--------------------------------------------------------------------------------------------------!
    47794675! Description:
    47804676! ------------
    47814677!> Autoconversion rate (Seifert and Beheng, 2006).
    4782 !------------------------------------------------------------------------------!
     4678!--------------------------------------------------------------------------------------------------!
    47834679    SUBROUTINE autoconversion
    47844680
     
    48144710                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    48154711
    4816                 IF ( microphysics_morrison ) THEN
     4712                IF ( microphysics_morrison )  THEN
    48174713                   nc_auto = nc(k,j,i)
    48184714                ELSE
     
    48264722!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
    48274723!--                (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
    4828                    tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) +         &
    4829                                     qc(k,j,i) ), 0.0_wp )
     4724                   tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) ), 0.0_wp )
    48304725!
    48314726!--                Universal function for autoconversion process
    48324727!--                (Seifert and Beheng, 2006):
    4833                    phi_au = 600.0_wp * tau_cloud**0.68_wp *                    &
    4834                             ( 1.0_wp - tau_cloud**0.68_wp )**3
     4728                   phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
    48354729!
    48364730!--                Shape parameter of gamma distribution (Geoffroy et al., 2010):
     
    48394733!
    48404734!--                Mean weight of cloud droplets:
    4841                    xc = MAX( hyrho(k) * qc(k,j,i) / nc_auto, xcmin) 
    4842 !
    4843 !--                Parameterized turbulence effects on autoconversion (Seifert,
    4844 !--                Nuijens and Stevens, 2010)
     4735                   xc = MAX( hyrho(k) * qc(k,j,i) / nc_auto, xcmin)
     4736!
     4737!--                Parameterized turbulence effects on autoconversion
     4738!--                (Seifert, Nuijens and Stevens, 2010)
    48454739                   IF ( collision_turbulence )  THEN
    48464740!
     
    48524746                      sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
    48534747!
    4854 !--                   Mixing length (neglecting distance to ground and
    4855 !--                   stratification)
     4748!--                   Mixing length (neglecting distance to ground and stratification)
    48564749                      l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
    48574750!
    4858 !--                   Limit dissipation rate according to Seifert, Nuijens and
    4859 !--                   Stevens (2010)
     4751!--                   Limit dissipation rate according to Seifert, Nuijens and Stevens (2010)
    48604752                      dissipation = MIN( 0.06_wp, diss(k,j,i) )
    48614753!
    48624754!--                   Compute Taylor-microscale Reynolds number:
    4863                       re_lambda = 6.0_wp / 11.0_wp *                           &
    4864                                   ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *   &
    4865                                   SQRT( 15.0_wp / kin_vis_air ) *              &
    4866                                   dissipation**( 1.0_wp / 6.0_wp )
    4867 !
    4868 !--                   The factor of 1.0E4 is needed to convert the dissipation
    4869 !--                   rate from m2 s-3 to cm2 s-3.
    4870                       k_au = k_au * ( 1.0_wp +                                 &
    4871                                       dissipation * 1.0E4_wp *                 &
    4872                                       ( re_lambda * 1.0E-3_wp )**0.25_wp *     &
    4873                                       ( alpha_cc * EXP( -1.0_wp * ( ( rc -     &
    4874                                                                       r_cc ) / &
    4875                                                         sigma_cc )**2          &
    4876                                                       ) + beta_cc              &
    4877                                       )                                        &
     4755                      re_lambda = 6.0_wp / 11.0_wp * ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *    &
     4756                                  SQRT( 15.0_wp / kin_vis_air ) * dissipation**( 1.0_wp / 6.0_wp )
     4757!
     4758!--                   The factor of 1.0E4 is needed to convert the dissipation rate from m2 s-3 to
     4759!--                   cm2 s-3.
     4760                      k_au = k_au * ( 1.0_wp +                                                     &
     4761                                      dissipation * 1.0E4_wp *                                     &
     4762                                      ( re_lambda * 1.0E-3_wp )**0.25_wp *                         &
     4763                                      ( alpha_cc * EXP( -1.0_wp * ( ( rc -                         &
     4764                                                                      r_cc ) /                     &
     4765                                                        sigma_cc )**2                              &
     4766                                                      ) + beta_cc                                  &
     4767                                      )                                                            &
    48784768                                    )
    48794769                   ENDIF
    48804770!
    48814771!--                Autoconversion rate (Seifert and Beheng, 2006):
    4882                    autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /    &
    4883                              ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *     &
    4884                              ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * &
    4885                              rho_surface
     4772                   autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp )    /                     &
     4773                             ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2     *                     &
     4774                             ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * rho_surface
    48864775                   autocon = MIN( autocon, qc(k,j,i) / dt_micro )
    48874776
    48884777                   qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag
    48894778                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag
    4890                    nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro  &
    4891                                                               * flag
    4892                    IF ( microphysics_morrison ) THEN
    4893                       nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *         &
    4894                                     autocon / x0 * hyrho(k) * dt_micro * flag )
     4779                   nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro * flag
     4780                   IF ( microphysics_morrison )  THEN
     4781                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *                             &
     4782                                  autocon / x0 * hyrho(k) * dt_micro * flag )
    48954783                   ENDIF
    48964784
     
    49064794
    49074795
    4908 !------------------------------------------------------------------------------!
     4796!--------------------------------------------------------------------------------------------------!
    49094797! Description:
    49104798! ------------
    49114799!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
    4912 !------------------------------------------------------------------------------!
     4800!--------------------------------------------------------------------------------------------------!
    49134801    SUBROUTINE autoconversion_ij( i, j )
    49144802
     
    49514839!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
    49524840!--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
    4953              tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) ),     &
    4954                               0.0_wp )
    4955 !
    4956 !--          Universal function for autoconversion process
    4957 !--          (Seifert and Beheng, 2006):
     4841             tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) ), 0.0_wp )
     4842!
     4843!--          Universal function for autoconversion process (Seifert and Beheng, 2006):
    49584844             phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
    49594845!
     
    49654851             xc = hyrho(k) * qc(k,j,i) / nc_auto
    49664852!
    4967 !--          Parameterized turbulence effects on autoconversion (Seifert,
    4968 !--          Nuijens and Stevens, 2010)
     4853!--          Parameterized turbulence effects on autoconversion (Seifert, Nuijens and Stevens, 2010)
    49694854             IF ( collision_turbulence )  THEN
    49704855!
     
    49794864                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
    49804865!
    4981 !--             Limit dissipation rate according to Seifert, Nuijens and
    4982 !--             Stevens (2010)
     4866!--             Limit dissipation rate according to Seifert, Nuijens and Stevens (2010)
    49834867                dissipation = MIN( 0.06_wp, diss(k,j,i) )
    49844868!
    49854869!--             Compute Taylor-microscale Reynolds number:
    4986                 re_lambda = 6.0_wp / 11.0_wp *                                 &
    4987                             ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *         &
    4988                             SQRT( 15.0_wp / kin_vis_air ) *                    &
    4989                             dissipation**( 1.0_wp / 6.0_wp )
    4990 !
    4991 !--             The factor of 1.0E4 is needed to convert the dissipation rate
    4992 !--             from m2 s-3 to cm2 s-3.
    4993                 k_au = k_au * ( 1.0_wp +                                       &
    4994                                 dissipation * 1.0E4_wp *                       &
    4995                                 ( re_lambda * 1.0E-3_wp )**0.25_wp *           &
    4996                                 ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /  &
    4997                                                   sigma_cc )**2                &
    4998                                                 ) + beta_cc                    &
    4999                                 )                                              &
     4870                re_lambda = 6.0_wp / 11.0_wp * ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *          &
     4871                            SQRT( 15.0_wp / kin_vis_air ) * dissipation**( 1.0_wp / 6.0_wp )
     4872!
     4873!--             The factor of 1.0E4 is needed to convert the dissipation rate from m2 s-3 to
     4874!--             cm2 s-3.
     4875                k_au = k_au * ( 1.0_wp +                                                           &
     4876                                dissipation * 1.0E4_wp *                                           &
     4877                                ( re_lambda * 1.0E-3_wp )**0.25_wp *                               &
     4878                                ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /                      &
     4879                                                  sigma_cc )**2                                    &
     4880                                                ) + beta_cc                                        &
     4881                                )                                                                  &
    50004882                              )
    50014883             ENDIF
    50024884!
    50034885!--          Autoconversion rate (Seifert and Beheng, 2006):
    5004              autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
    5005                        ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *           &
    5006                        ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
    5007                        rho_surface
     4886             autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp )    /                           &
     4887                       ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2     *                           &
     4888                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * rho_surface
    50084889             autocon = MIN( autocon, qc(k,j,i) / dt_micro )
    5009 
    50104890             qr(k,j,i) = qr(k,j,i) + autocon * dt_micro                 * flag
    50114891             qc(k,j,i) = qc(k,j,i) - autocon * dt_micro                 * flag
    50124892             nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro * flag
    5013              IF ( microphysics_morrison ) THEN
    5014                 nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *                &
    5015                                     autocon / x0 * hyrho(k) * dt_micro * flag )
     4893             IF ( microphysics_morrison )  THEN
     4894                nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *                                   &
     4895                            autocon / x0 * hyrho(k) * dt_micro * flag )
    50164896             ENDIF
    50174897
     
    50234903
    50244904
    5025 !------------------------------------------------------------------------------!
     4905!--------------------------------------------------------------------------------------------------!
    50264906! Description:
    50274907! ------------
    50284908!> Autoconversion process (Kessler, 1969).
    5029 !------------------------------------------------------------------------------!
     4909!--------------------------------------------------------------------------------------------------!
    50304910    SUBROUTINE autoconversion_kessler
    50314911
     
    50594939                qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
    50604940                q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro * flag
    5061                 pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp *     &
    5062                                         d_exner(k)              * flag
     4941                pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp *                         &
     4942                                        d_exner(k)             * flag
    50634943
    50644944!
     
    50724952   END SUBROUTINE autoconversion_kessler
    50734953
    5074 !------------------------------------------------------------------------------!
     4954!--------------------------------------------------------------------------------------------------!
    50754955! Description:
    50764956! ------------
    50774957!> Autoconversion process (Kessler, 1969).
    5078 !------------------------------------------------------------------------------!
     4958!--------------------------------------------------------------------------------------------------!
    50794959    SUBROUTINE autoconversion_kessler_ij( i, j )
    50804960
     
    50874967       INTEGER(iwp) ::  k_wall !< topography top index
    50884968
    5089        REAL(wp)    ::  dqdt_precip !<
     4969       REAL(wp)    ::  dqdt_precip       !<
    50904970       REAL(wp)    ::  flag              !< flag to indicate first grid level above surface
    50914971
     
    51064986          qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
    51074987          q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro * flag
    5108           pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp * d_exner(k)  &
    5109                                                                  * flag
     4988          pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp * d_exner(k) * flag
    51104989
    51114990!
     
    51184997
    51194998
    5120 !------------------------------------------------------------------------------!
     4999!--------------------------------------------------------------------------------------------------!
    51215000! Description:
    51225001! ------------
    51235002!> Accretion rate (Seifert and Beheng, 2006).
    5124 !------------------------------------------------------------------------------!
     5003!--------------------------------------------------------------------------------------------------!
    51255004    SUBROUTINE accretion
    51265005
     
    51495028                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    51505029
    5151                 IF ( microphysics_morrison ) THEN
     5030                IF ( microphysics_morrison )  THEN
    51525031                   nc_accr = nc(k,j,i)
    51535032                ELSE
     
    51555034                ENDIF
    51565035
    5157                 IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )    &
    5158                                              .AND.  ( nc_accr > eps_mr ) ) THEN
     5036                IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )                         &
     5037                                             .AND.  ( nc_accr > eps_mr   ) ) THEN
    51595038!
    51605039!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
     
    51695048                   xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
    51705049!
    5171 !--                Parameterized turbulence effects on autoconversion (Seifert,
    5172 !--                Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
    5173 !--                convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
     5050!--                Parameterized turbulence effects on autoconversion
     5051!--                (Seifert, Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to convert
     5052!--                the dissipation rate (diss) from m2 s-3 to cm2 s-3.
    51745053                   IF ( collision_turbulence )  THEN
    5175                       k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                      &
    5176 <