Index: palm/trunk/SOURCE/bulk_cloud_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/bulk_cloud_model_mod.f90	(revision 3786)
+++ palm/trunk/SOURCE/bulk_cloud_model_mod.f90	(revision 3869)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! moving the furniture around ;-)
+! 
+! 3786 2019-03-06 16:58:03Z raasch
 ! unsed variables removed
 ! 
@@ -325,7 +328,8 @@
            bcm_check_data_output, &
            bcm_check_data_output_pr, &
-           bcm_header, &
            bcm_init_arrays, &
            bcm_init, &
+           bcm_header, &
+           bcm_actions, &
            bcm_3d_data_averaging, &
            bcm_data_output_2d, &
@@ -336,5 +340,4 @@
            bcm_wrd_global, &
            bcm_wrd_local, &
-           bcm_actions, &
            calc_liquid_water_content
 
@@ -370,47 +373,15 @@
     END INTERFACE bcm_check_data_output_pr
 
+    INTERFACE bcm_init_arrays
+       MODULE PROCEDURE bcm_init_arrays
+    END INTERFACE bcm_init_arrays
+
+    INTERFACE bcm_init
+       MODULE PROCEDURE bcm_init
+    END INTERFACE bcm_init
+
     INTERFACE bcm_header
        MODULE PROCEDURE bcm_header
     END INTERFACE bcm_header
-
-    INTERFACE bcm_init_arrays
-       MODULE PROCEDURE bcm_init_arrays
-    END INTERFACE bcm_init_arrays
-
-    INTERFACE bcm_init
-       MODULE PROCEDURE bcm_init
-    END INTERFACE bcm_init
-
-    INTERFACE bcm_3d_data_averaging
-       MODULE PROCEDURE bcm_3d_data_averaging
-    END INTERFACE bcm_3d_data_averaging
-
-    INTERFACE bcm_data_output_2d
-       MODULE PROCEDURE bcm_data_output_2d
-    END INTERFACE bcm_data_output_2d
-
-    INTERFACE bcm_data_output_3d
-       MODULE PROCEDURE bcm_data_output_3d
-    END INTERFACE bcm_data_output_3d
-
-    INTERFACE bcm_swap_timelevel
-       MODULE PROCEDURE bcm_swap_timelevel
-    END INTERFACE bcm_swap_timelevel
-
-    INTERFACE bcm_rrd_global
-       MODULE PROCEDURE bcm_rrd_global
-    END INTERFACE bcm_rrd_global
-
-    INTERFACE bcm_rrd_local
-       MODULE PROCEDURE bcm_rrd_local
-    END INTERFACE bcm_rrd_local
-
-    INTERFACE bcm_wrd_global
-       MODULE PROCEDURE bcm_wrd_global
-    END INTERFACE bcm_wrd_global
-
-    INTERFACE bcm_wrd_local
-       MODULE PROCEDURE bcm_wrd_local
-    END INTERFACE bcm_wrd_local
 
     INTERFACE bcm_actions
@@ -419,62 +390,35 @@
     END INTERFACE bcm_actions
 
-    INTERFACE adjust_cloud
-       MODULE PROCEDURE adjust_cloud
-       MODULE PROCEDURE adjust_cloud_ij
-    END INTERFACE adjust_cloud
-
-    INTERFACE activation
-       MODULE PROCEDURE activation
-       MODULE PROCEDURE activation_ij
-    END INTERFACE activation
-
-    INTERFACE condensation
-       MODULE PROCEDURE condensation
-       MODULE PROCEDURE condensation_ij
-    END INTERFACE condensation
-
-    INTERFACE autoconversion
-       MODULE PROCEDURE autoconversion
-       MODULE PROCEDURE autoconversion_ij
-    END INTERFACE autoconversion
-
-    INTERFACE autoconversion_kessler
-       MODULE PROCEDURE autoconversion_kessler
-       MODULE PROCEDURE autoconversion_kessler_ij
-    END INTERFACE autoconversion_kessler
-
-    INTERFACE accretion
-       MODULE PROCEDURE accretion
-       MODULE PROCEDURE accretion_ij
-    END INTERFACE accretion
-
-    INTERFACE selfcollection_breakup
-       MODULE PROCEDURE selfcollection_breakup
-       MODULE PROCEDURE selfcollection_breakup_ij
-    END INTERFACE selfcollection_breakup
-
-    INTERFACE evaporation_rain
-       MODULE PROCEDURE evaporation_rain
-       MODULE PROCEDURE evaporation_rain_ij
-    END INTERFACE evaporation_rain
-
-    INTERFACE sedimentation_cloud
-       MODULE PROCEDURE sedimentation_cloud
-       MODULE PROCEDURE sedimentation_cloud_ij
-    END INTERFACE sedimentation_cloud
-
-    INTERFACE sedimentation_rain
-       MODULE PROCEDURE sedimentation_rain
-       MODULE PROCEDURE sedimentation_rain_ij
-    END INTERFACE sedimentation_rain
-
-    INTERFACE calc_precipitation_amount
-       MODULE PROCEDURE calc_precipitation_amount
-       MODULE PROCEDURE calc_precipitation_amount_ij
-    END INTERFACE calc_precipitation_amount
-
-    INTERFACE supersaturation
-       MODULE PROCEDURE supersaturation
-    END INTERFACE supersaturation
+    INTERFACE bcm_swap_timelevel
+       MODULE PROCEDURE bcm_swap_timelevel
+    END INTERFACE bcm_swap_timelevel
+
+    INTERFACE bcm_3d_data_averaging
+       MODULE PROCEDURE bcm_3d_data_averaging
+    END INTERFACE bcm_3d_data_averaging
+
+    INTERFACE bcm_data_output_2d
+       MODULE PROCEDURE bcm_data_output_2d
+    END INTERFACE bcm_data_output_2d
+
+    INTERFACE bcm_data_output_3d
+       MODULE PROCEDURE bcm_data_output_3d
+    END INTERFACE bcm_data_output_3d
+
+    INTERFACE bcm_rrd_global
+       MODULE PROCEDURE bcm_rrd_global
+    END INTERFACE bcm_rrd_global
+
+    INTERFACE bcm_rrd_local
+       MODULE PROCEDURE bcm_rrd_local
+    END INTERFACE bcm_rrd_local
+
+    INTERFACE bcm_wrd_global
+       MODULE PROCEDURE bcm_wrd_global
+    END INTERFACE bcm_wrd_global
+
+    INTERFACE bcm_wrd_local
+       MODULE PROCEDURE bcm_wrd_local
+    END INTERFACE bcm_wrd_local
 
     INTERFACE calc_liquid_water_content
@@ -989,4 +933,194 @@
 ! Description:
 ! ------------
+!> Header output for bulk cloud module
+!------------------------------------------------------------------------------!
+    SUBROUTINE bcm_header ( io )
+
+
+       IMPLICIT NONE
+
+       INTEGER(iwp), INTENT(IN) ::  io  !< Unit of the output file
+
+!
+!--    Write bulk cloud module header
+       WRITE ( io, 1 )
+
+       WRITE ( io, 2 )
+       WRITE ( io, 3 )
+
+       IF ( microphysics_kessler )  THEN
+          WRITE ( io, 4 ) 'Kessler-Scheme'
+       ENDIF
+
+       IF ( microphysics_seifert )  THEN
+          WRITE ( io, 4 ) 'Seifert-Beheng-Scheme'
+          IF ( cloud_water_sedimentation )  WRITE ( io, 5 )
+          IF ( collision_turbulence )  WRITE ( io, 6 )
+          IF ( ventilation_effect )  WRITE ( io, 7 )
+          IF ( limiter_sedimentation )  WRITE ( io, 8 )
+       ENDIF
+
+       WRITE ( io, 20 )
+       WRITE ( io, 21 ) surface_pressure
+       WRITE ( io, 22 ) r_d
+       WRITE ( io, 23 ) rho_surface
+       WRITE ( io, 24 ) c_p
+       WRITE ( io, 25 ) l_v
+
+       IF ( microphysics_seifert )  THEN
+          WRITE ( io, 26 ) 1.0E-6_wp * nc_const
+          WRITE ( io, 27 ) c_sedimentation
+       ENDIF
+
+
+1   FORMAT ( //' Bulk cloud module information:'/ &
+               ' ------------------------------------------'/ )
+2   FORMAT ( '--> Bulk scheme with liquid water potential temperature and'/ &
+             '    total water content is used.' )
+3   FORMAT ( '--> Condensation is parameterized via 0% - or 100% scheme.' )
+4   FORMAT ( '--> Precipitation parameterization via ', A )
+
+5   FORMAT ( '--> Cloud water sedimentation parameterization via Stokes law' )
+6   FORMAT ( '--> Turbulence effects on precipitation process' )
+7   FORMAT ( '--> Ventilation effects on evaporation of rain drops' )
+8   FORMAT ( '--> Slope limiter used for sedimentation process' )
+
+20  FORMAT ( '--> Essential parameters:' )
+21  FORMAT ( '       Surface pressure             :   p_0   = ', F7.2, ' hPa')
+22  FORMAT ( '       Gas constant                 :   R     = ', F5.1, ' J/(kg K)')
+23  FORMAT ( '       Density of air               :   rho_0 = ', F6.3, ' kg/m**3')
+24  FORMAT ( '       Specific heat cap.           :   c_p   = ', F6.1, ' J/(kg K)')
+25  FORMAT ( '       Vapourization heat           :   L_v   = ', E9.2, ' J/kg')
+26  FORMAT ( '       Droplet density              :   N_c   = ', F6.1, ' 1/cm**3' )
+27  FORMAT ( '       Sedimentation Courant number :   C_s   = ', F4.1 )
+
+
+    END SUBROUTINE bcm_header
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Control of microphysics for all grid points
+!------------------------------------------------------------------------------!
+    SUBROUTINE bcm_actions
+
+       IMPLICIT NONE
+
+       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
+!
+!--       Calculate vertical profile of the hydrostatic pressure (hyp)
+          hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
+          d_exner = exner_function_invers(hyp)
+          exner = 1.0_wp / exner_function_invers(hyp)
+          hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
+!
+!--       Compute reference density
+          rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
+       ENDIF
+
+!
+!--    Compute length of time step
+       IF ( call_microphysics_at_all_substeps )  THEN
+          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
+       ELSE
+          dt_micro = dt_3d
+       ENDIF
+
+!
+!--    Reset precipitation rate
+       IF ( intermediate_timestep_count == 1 )  prr = 0.0_wp
+
+!
+!--    Compute cloud physics
+       IF ( microphysics_kessler )  THEN
+
+          CALL autoconversion_kessler
+          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
+
+       ELSEIF ( microphysics_seifert )  THEN
+
+          CALL adjust_cloud
+          IF ( microphysics_morrison )  CALL activation
+          IF ( microphysics_morrison )  CALL condensation
+          CALL autoconversion
+          CALL accretion
+          CALL selfcollection_breakup
+          CALL evaporation_rain
+          CALL sedimentation_rain
+          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
+
+       ENDIF
+
+       CALL calc_precipitation_amount
+
+    END SUBROUTINE bcm_actions
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Control of microphysics for grid points i,j
+!------------------------------------------------------------------------------!
+
+    SUBROUTINE bcm_actions_ij( i, j )
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !<
+       INTEGER(iwp) ::  j                 !<
+
+       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
+!
+!--       Calculate vertical profile of the hydrostatic pressure (hyp)
+          hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
+          d_exner = exner_function_invers(hyp)
+          exner = 1.0_wp / exner_function_invers(hyp)
+          hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
+!
+!--       Compute reference density
+          rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
+       ENDIF
+
+!
+!--    Compute length of time step
+       IF ( call_microphysics_at_all_substeps )  THEN
+          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
+       ELSE
+          dt_micro = dt_3d
+       ENDIF
+!
+!--    Reset precipitation rate
+       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
+
+!
+!--    Compute cloud physics
+       IF( microphysics_kessler )  THEN
+
+          CALL autoconversion_kessler_ij( i,j )
+          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud_ij( i,j )
+
+       ELSEIF ( microphysics_seifert )  THEN
+
+          CALL adjust_cloud_ij( i,j )
+          IF ( microphysics_morrison ) CALL activation_ij( i,j )
+          IF ( microphysics_morrison ) CALL condensation_ij( i,j )
+          CALL autoconversion_ij( i,j )
+          CALL accretion_ij( i,j )
+          CALL selfcollection_breakup_ij( i,j )
+          CALL evaporation_rain_ij( i,j )
+          CALL sedimentation_rain_ij( i,j )
+          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud_ij( i,j )
+
+       ENDIF
+
+       CALL calc_precipitation_amount_ij( i,j )
+
+    END SUBROUTINE bcm_actions_ij
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Swapping of timelevels
 !------------------------------------------------------------------------------!
@@ -1028,73 +1162,4 @@
 
     END SUBROUTINE bcm_swap_timelevel
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Header output for bulk cloud module
-!------------------------------------------------------------------------------!
-    SUBROUTINE bcm_header ( io )
-
-
-       IMPLICIT NONE
-
-       INTEGER(iwp), INTENT(IN) ::  io  !< Unit of the output file
-
-!
-!--    Write bulk cloud module header
-       WRITE ( io, 1 )
-
-       WRITE ( io, 2 )
-       WRITE ( io, 3 )
-
-       IF ( microphysics_kessler )  THEN
-          WRITE ( io, 4 ) 'Kessler-Scheme'
-       ENDIF
-
-       IF ( microphysics_seifert )  THEN
-          WRITE ( io, 4 ) 'Seifert-Beheng-Scheme'
-          IF ( cloud_water_sedimentation )  WRITE ( io, 5 )
-          IF ( collision_turbulence )  WRITE ( io, 6 )
-          IF ( ventilation_effect )  WRITE ( io, 7 )
-          IF ( limiter_sedimentation )  WRITE ( io, 8 )
-       ENDIF
-
-       WRITE ( io, 20 )
-       WRITE ( io, 21 ) surface_pressure
-       WRITE ( io, 22 ) r_d
-       WRITE ( io, 23 ) rho_surface
-       WRITE ( io, 24 ) c_p
-       WRITE ( io, 25 ) l_v
-
-       IF ( microphysics_seifert )  THEN
-          WRITE ( io, 26 ) 1.0E-6_wp * nc_const
-          WRITE ( io, 27 ) c_sedimentation
-       ENDIF
-
-
-1   FORMAT ( //' Bulk cloud module information:'/ &
-               ' ------------------------------------------'/ )
-2   FORMAT ( '--> Bulk scheme with liquid water potential temperature and'/ &
-             '    total water content is used.' )
-3   FORMAT ( '--> Condensation is parameterized via 0% - or 100% scheme.' )
-4   FORMAT ( '--> Precipitation parameterization via ', A )
-
-5   FORMAT ( '--> Cloud water sedimentation parameterization via Stokes law' )
-6   FORMAT ( '--> Turbulence effects on precipitation process' )
-7   FORMAT ( '--> Ventilation effects on evaporation of rain drops' )
-8   FORMAT ( '--> Slope limiter used for sedimentation process' )
-
-20  FORMAT ( '--> Essential parameters:' )
-21  FORMAT ( '       Surface pressure             :   p_0   = ', F7.2, ' hPa')
-22  FORMAT ( '       Gas constant                 :   R     = ', F5.1, ' J/(kg K)')
-23  FORMAT ( '       Density of air               :   rho_0 = ', F6.3, ' kg/m**3')
-24  FORMAT ( '       Specific heat cap.           :   c_p   = ', F6.1, ' J/(kg K)')
-25  FORMAT ( '       Vapourization heat           :   L_v   = ', E9.2, ' J/kg')
-26  FORMAT ( '       Droplet density              :   N_c   = ', F6.1, ' 1/cm**3' )
-27  FORMAT ( '       Sedimentation Courant number :   C_s   = ', F4.1 )
-
-
-    END SUBROUTINE bcm_header
 
 
@@ -1993,63 +2058,4 @@
     END SUBROUTINE bcm_wrd_local
 
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Control of microphysics for all grid points
-!------------------------------------------------------------------------------!
-    SUBROUTINE bcm_actions
-
-       IMPLICIT NONE
-
-       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
-!
-!--       Calculate vertical profile of the hydrostatic pressure (hyp)
-          hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
-          d_exner = exner_function_invers(hyp)
-          exner = 1.0_wp / exner_function_invers(hyp)
-          hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
-!
-!--       Compute reference density
-          rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
-       ENDIF
-
-!
-!--    Compute length of time step
-       IF ( call_microphysics_at_all_substeps )  THEN
-          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
-       ELSE
-          dt_micro = dt_3d
-       ENDIF
-
-!
-!--    Reset precipitation rate
-       IF ( intermediate_timestep_count == 1 )  prr = 0.0_wp
-
-!
-!--    Compute cloud physics
-       IF ( microphysics_kessler )  THEN
-
-          CALL autoconversion_kessler
-          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
-
-       ELSEIF ( microphysics_seifert )  THEN
-
-          CALL adjust_cloud
-          IF ( microphysics_morrison )  CALL activation
-          IF ( microphysics_morrison )  CALL condensation
-          CALL autoconversion
-          CALL accretion
-          CALL selfcollection_breakup
-          CALL evaporation_rain
-          CALL sedimentation_rain
-          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
-
-       ENDIF
-
-       CALL calc_precipitation_amount
-
-    END SUBROUTINE bcm_actions
-
 !------------------------------------------------------------------------------!
 ! Description:
@@ -2107,4 +2113,59 @@
 
     END SUBROUTINE adjust_cloud
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Adjust number of raindrops to avoid nonlinear effects in
+!> sedimentation and evaporation of rain drops due to too small or
+!> too big weights of rain drops (Stevens and Seifert, 2008).
+!> The same procedure is applied to cloud droplets if they are determined
+!> prognostically. Call for grid point i,j
+!------------------------------------------------------------------------------!
+    SUBROUTINE adjust_cloud_ij( i, j )
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !<
+       INTEGER(iwp) ::  j                 !<
+       INTEGER(iwp) ::  k                 !<
+
+       REAL(wp) ::  flag                  !< flag to indicate first grid level above surface
+
+       DO  k = nzb+1, nzt
+!
+!--       Predetermine flag to mask topography
+          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+
+          IF ( qr(k,j,i) <= eps_sb )  THEN
+             qr(k,j,i) = 0.0_wp
+             nr(k,j,i) = 0.0_wp
+          ELSE
+!
+!--          Adjust number of raindrops to avoid nonlinear effects in
+!--          sedimentation and evaporation of rain drops due to too small or
+!--          too big weights of rain drops (Stevens and Seifert, 2008).
+             IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
+                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
+             ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
+                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
+             ENDIF
+
+          ENDIF
+
+          IF ( microphysics_morrison ) THEN
+             IF ( qc(k,j,i) <= eps_sb )  THEN
+                qc(k,j,i) = 0.0_wp
+                nc(k,j,i) = 0.0_wp
+             ELSE
+                IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
+                   nc(k,j,i) = qc(k,j,i) * hyrho(k) / xamin * flag
+                ENDIF
+             ENDIF
+          ENDIF
+
+       ENDDO
+
+    END SUBROUTINE adjust_cloud_ij
 
 !------------------------------------------------------------------------------!
@@ -2214,4 +2275,101 @@
     END SUBROUTINE activation
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate number of activated condensation nucleii after simple activation
+!> scheme of Twomey, 1959.
+!------------------------------------------------------------------------------!
+    SUBROUTINE activation_ij( i, j )
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !<
+       INTEGER(iwp) ::  j                 !<
+       INTEGER(iwp) ::  k                 !<
+
+       REAL(wp)     ::  activ             !<
+       REAL(wp)     ::  afactor           !<
+       REAL(wp)     ::  beta_act          !<
+       REAL(wp)     ::  bfactor           !<
+       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
+       REAL(wp)     ::  k_act             !<
+       REAL(wp)     ::  n_act             !<
+       REAL(wp)     ::  n_ccn             !<
+       REAL(wp)     ::  s_0               !<
+       REAL(wp)     ::  sat_max           !<
+       REAL(wp)     ::  sigma             !<
+       REAL(wp)     ::  sigma_act         !<
+
+       DO  k = nzb+1, nzt
+!
+!--       Predetermine flag to mask topography
+          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+!
+!--       Call calculation of supersaturation
+          CALL supersaturation ( i, j, k )
+!
+!--       Prescribe parameters for activation
+!--       (see: Bott + Trautmann, 2002, Atm. Res., 64)
+          k_act  = 0.7_wp
+          activ  = 0.0_wp
+
+          IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk )  THEN
+!
+!--          Compute the number of activated Aerosols
+!--          (see: Twomey, 1959, Pure and applied Geophysics, 43)
+             n_act     = na_init * sat**k_act
+!
+!--          Compute the number of cloud droplets
+!--          (see: Morrison + Grabowski, 2007, JAS, 64)
+!            activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro
+
+!
+!--          Compute activation rate after Khairoutdinov and Kogan
+!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
+             sat_max = 0.8_wp / 100.0_wp
+             activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN             &
+                 ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /         &
+                  dt_micro
+
+             nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)
+          ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk )  THEN
+!
+!--          Curvature effect (afactor) with surface tension
+!--          parameterization by Straka (2009)
+             sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
+             afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
+!
+!--          Solute effect (bfactor)
+             bfactor = vanthoff * molecular_weight_of_water *                  &
+                  rho_s / ( molecular_weight_of_solute * rho_l )
+
+!
+!--          Prescribe power index that describes the soluble fraction
+!--          of an aerosol particle (beta).
+!--          (see: Morrison + Grabowski, 2007, JAS, 64)
+             beta_act  = 0.5_wp
+             sigma_act = sigma_bulk**( 1.0_wp + beta_act )
+!
+!--          Calculate mean geometric supersaturation (s_0) with
+!--          parameterization by Khvorostyanov and Curry (2006)
+             s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *         &
+               ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
+
+!
+!--          Calculate number of activated CCN as a function of
+!--          supersaturation and taking Koehler theory into account
+!--          (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
+             n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                    &
+                LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
+             activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
+
+             nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro * flag), na_init)
+          ENDIF
+
+       ENDDO
+
+    END SUBROUTINE activation_ij
+
 
 !------------------------------------------------------------------------------!
@@ -2297,4 +2455,76 @@
 
     END SUBROUTINE condensation
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate condensation rate for cloud water content (after Khairoutdinov and
+!> Kogan, 2000).
+!------------------------------------------------------------------------------!
+    SUBROUTINE condensation_ij( i, j )
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !<
+       INTEGER(iwp) ::  j                 !<
+       INTEGER(iwp) ::  k                 !<
+
+       REAL(wp)     ::  cond              !<
+       REAL(wp)     ::  cond_max          !<
+       REAL(wp)     ::  dc                !<
+       REAL(wp)     ::  evap              !<
+       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
+       REAL(wp)     ::  g_fac             !<
+       REAL(wp)     ::  nc_0              !<
+       REAL(wp)     ::  temp              !<
+       REAL(wp)     ::  xc                !<
+
+
+       DO  k = nzb+1, nzt
+!
+!--       Predetermine flag to mask topography
+          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+!
+!--       Call calculation of supersaturation
+          CALL supersaturation ( i, j, k )
+!
+!--       Actual temperature:
+          temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
+
+          g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
+                              l_v / ( thermal_conductivity_l * temp )    &
+                              + r_v * temp / ( diff_coeff_l * e_s )      &
+                            )
+!
+!--       Mean weight of cloud drops
+          IF ( nc(k,j,i) <= 0.0_wp) CYCLE
+          xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
+!
+!--       Weight averaged diameter of cloud drops:
+          dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
+!
+!--       Integral diameter of cloud drops
+          nc_0 = nc(k,j,i) * dc
+!
+!--       Condensation needs only to be calculated in supersaturated regions
+          IF ( sat > 0.0_wp )  THEN
+!
+!--          Condensation rate of cloud water content
+!--          after KK scheme.
+!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
+             cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
+             cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
+             cond      = MIN( cond, cond_max / dt_micro )
+
+             qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag
+          ELSEIF ( sat < 0.0_wp ) THEN
+             evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
+             evap      = MAX( evap, -qc(k,j,i) / dt_micro )
+
+             qc(k,j,i) = qc(k,j,i) + evap * dt_micro * flag
+          ENDIF
+       ENDDO
+
+    END SUBROUTINE condensation_ij
 
 
@@ -2432,4 +2662,121 @@
 ! Description:
 ! ------------
+!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
+!------------------------------------------------------------------------------!
+    SUBROUTINE autoconversion_ij( i, j )
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !<
+       INTEGER(iwp) ::  j                 !<
+       INTEGER(iwp) ::  k                 !<
+
+       REAL(wp)     ::  alpha_cc          !<
+       REAL(wp)     ::  autocon           !<
+       REAL(wp)     ::  dissipation       !<
+       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
+       REAL(wp)     ::  k_au              !<
+       REAL(wp)     ::  l_mix             !<
+       REAL(wp)     ::  nc_auto           !<
+       REAL(wp)     ::  nu_c              !<
+       REAL(wp)     ::  phi_au            !<
+       REAL(wp)     ::  r_cc              !<
+       REAL(wp)     ::  rc                !<
+       REAL(wp)     ::  re_lambda         !<
+       REAL(wp)     ::  sigma_cc          !<
+       REAL(wp)     ::  tau_cloud         !<
+       REAL(wp)     ::  xc                !<
+
+       DO  k = nzb+1, nzt
+!
+!--       Predetermine flag to mask topography
+          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+          IF ( microphysics_morrison ) THEN
+             nc_auto = nc(k,j,i)
+          ELSE
+             nc_auto = nc_const
+          ENDIF
+
+          IF ( qc(k,j,i) > eps_sb  .AND.  nc_auto > eps_mr )  THEN
+
+             k_au = k_cc / ( 20.0_wp * x0 )
+!
+!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
+!--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
+             tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) ),     &
+                              0.0_wp )
+!
+!--          Universal function for autoconversion process
+!--          (Seifert and Beheng, 2006):
+             phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
+!
+!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
+!--          (Use constant nu_c = 1.0_wp instead?)
+             nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
+!
+!--          Mean weight of cloud droplets:
+             xc = hyrho(k) * qc(k,j,i) / nc_auto
+!
+!--          Parameterized turbulence effects on autoconversion (Seifert,
+!--          Nuijens and Stevens, 2010)
+             IF ( collision_turbulence )  THEN
+!
+!--             Weight averaged radius of cloud droplets:
+                rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
+
+                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
+                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
+                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
+!
+!--             Mixing length (neglecting distance to ground and stratification)
+                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
+!
+!--             Limit dissipation rate according to Seifert, Nuijens and
+!--             Stevens (2010)
+                dissipation = MIN( 0.06_wp, diss(k,j,i) )
+!
+!--             Compute Taylor-microscale Reynolds number:
+                re_lambda = 6.0_wp / 11.0_wp *                                 &
+                            ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *         &
+                            SQRT( 15.0_wp / kin_vis_air ) *                    &
+                            dissipation**( 1.0_wp / 6.0_wp )
+!
+!--             The factor of 1.0E4 is needed to convert the dissipation rate
+!--             from m2 s-3 to cm2 s-3.
+                k_au = k_au * ( 1.0_wp +                                       &
+                                dissipation * 1.0E4_wp *                       &
+                                ( re_lambda * 1.0E-3_wp )**0.25_wp *           &
+                                ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /  &
+                                                  sigma_cc )**2                &
+                                                ) + beta_cc                    &
+                                )                                              &
+                              )
+             ENDIF
+!
+!--          Autoconversion rate (Seifert and Beheng, 2006):
+             autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
+                       ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *            &
+                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
+                       rho_surface
+             autocon = MIN( autocon, qc(k,j,i) / dt_micro )
+
+             qr(k,j,i) = qr(k,j,i) + autocon * dt_micro                 * flag
+             qc(k,j,i) = qc(k,j,i) - autocon * dt_micro                 * flag
+             nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro * flag
+             IF ( microphysics_morrison ) THEN
+                nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *                &
+                                    autocon / x0 * hyrho(k) * dt_micro * flag )
+             ENDIF
+
+          ENDIF
+
+       ENDDO
+
+    END SUBROUTINE autoconversion_ij
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Autoconversion process (Kessler, 1969).
 !------------------------------------------------------------------------------!
@@ -2477,4 +2824,49 @@
 
    END SUBROUTINE autoconversion_kessler
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Autoconversion process (Kessler, 1969).
+!------------------------------------------------------------------------------!
+    SUBROUTINE autoconversion_kessler_ij( i, j )
+
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i      !<
+       INTEGER(iwp) ::  j      !<
+       INTEGER(iwp) ::  k      !<
+       INTEGER(iwp) ::  k_wall !< topography top index
+
+       REAL(wp)    ::  dqdt_precip !<
+       REAL(wp)    ::  flag              !< flag to indicate first grid level above surface
+
+!
+!--    Determine vertical index of topography top
+       k_wall = get_topography_top_index_ji( j, i, 's' )
+       DO  k = nzb+1, nzt
+!
+!--       Predetermine flag to mask topography
+          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+
+          IF ( qc(k,j,i) > ql_crit )  THEN
+             dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
+          ELSE
+             dqdt_precip = 0.0_wp
+          ENDIF
+
+          qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
+          q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro * flag
+          pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp * d_exner(k)  &
+                                                                 * flag
+
+!
+!--       Compute the rain rate (stored on surface grid point)
+          prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
+
+       ENDDO
+
+    END SUBROUTINE autoconversion_kessler_ij
 
 
@@ -2564,4 +2956,82 @@
     END SUBROUTINE accretion
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
+!------------------------------------------------------------------------------!
+    SUBROUTINE accretion_ij( i, j )
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !<
+       INTEGER(iwp) ::  j                 !<
+       INTEGER(iwp) ::  k                 !<
+
+       REAL(wp)     ::  accr              !<
+       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
+       REAL(wp)     ::  k_cr              !<
+       REAL(wp)     ::  nc_accr           !<
+       REAL(wp)     ::  phi_ac            !<
+       REAL(wp)     ::  tau_cloud         !<
+       REAL(wp)     ::  xc                !<
+
+
+       DO  k = nzb+1, nzt
+!
+!--       Predetermine flag to mask topography
+          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+          IF ( microphysics_morrison ) THEN
+             nc_accr = nc(k,j,i)
+          ELSE
+             nc_accr = nc_const
+          ENDIF
+
+          IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )  .AND.  &
+               ( nc_accr > eps_mr ) )  THEN
+!
+!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
+             tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
+!
+!--          Universal function for accretion process
+!--          (Seifert and Beheng, 2001):
+             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
+
+!
+!--          Mean weight of cloud drops
+             xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
+!
+!--          Parameterized turbulence effects on autoconversion (Seifert,
+!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
+!--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
+             IF ( collision_turbulence )  THEN
+                k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                            &
+                                 MIN( 600.0_wp,                                &
+                                      diss(k,j,i) * 1.0E4_wp )**0.25_wp        &
+                               )
+             ELSE
+                k_cr = k_cr0
+             ENDIF
+!
+!--          Accretion rate (Seifert and Beheng, 2006):
+             accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *                      &
+                    SQRT( rho_surface * hyrho(k) )
+             accr = MIN( accr, qc(k,j,i) / dt_micro )
+
+             qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
+             qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
+             IF ( microphysics_morrison )  THEN
+                nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc *               &
+                                             hyrho(k) * dt_micro * flag        &
+                                         )
+             ENDIF
+
+
+          ENDIF
+
+       ENDDO
+
+    END SUBROUTINE accretion_ij
+
 
 !------------------------------------------------------------------------------!
@@ -2622,4 +3092,53 @@
 
     END SUBROUTINE selfcollection_breakup
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
+!------------------------------------------------------------------------------!
+    SUBROUTINE selfcollection_breakup_ij( i, j )
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !<
+       INTEGER(iwp) ::  j                 !<
+       INTEGER(iwp) ::  k                 !<
+
+       REAL(wp)     ::  breakup           !<
+       REAL(wp)     ::  dr                !<
+       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
+       REAL(wp)     ::  phi_br            !<
+       REAL(wp)     ::  selfcoll          !<
+
+       DO  k = nzb+1, nzt
+!
+!--       Predetermine flag to mask topography
+          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+
+          IF ( qr(k,j,i) > eps_sb )  THEN
+!
+!--          Selfcollection rate (Seifert and Beheng, 2001):
+             selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * SQRT( hyrho(k) * rho_surface )
+!
+!--          Weight averaged diameter of rain drops:
+             dr = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
+!
+!--          Collisional breakup rate (Seifert, 2008):
+             IF ( dr >= 0.3E-3_wp )  THEN
+                phi_br  = k_br * ( dr - 1.1E-3_wp )
+                breakup = selfcoll * ( phi_br + 1.0_wp )
+             ELSE
+                breakup = 0.0_wp
+             ENDIF
+
+             selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
+             nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
+
+          ENDIF
+       ENDDO
+
+    END SUBROUTINE selfcollection_breakup_ij
 
 
@@ -2753,958 +3272,4 @@
 ! Description:
 ! ------------
-!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
-!------------------------------------------------------------------------------!
-    SUBROUTINE sedimentation_cloud
-
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i             !<
-       INTEGER(iwp) ::  j             !<
-       INTEGER(iwp) ::  k             !<
-
-       REAL(wp) ::  flag              !< flag to mask topography grid points
-       REAL(wp) ::  nc_sedi           !<
-
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc !<
-
-
-       CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
-
-       sed_qc(nzt+1) = 0.0_wp
-       sed_nc(nzt+1) = 0.0_wp
-
-       DO  i = nxlg, nxrg
-          DO  j = nysg, nyng
-             DO  k = nzt, nzb+1, -1
-!
-!--             Predetermine flag to mask topography
-                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-
-                IF ( microphysics_morrison ) THEN
-                   nc_sedi = nc(k,j,i)
-                ELSE
-                   nc_sedi = nc_const
-                ENDIF
-
-!
-!--             Sedimentation fluxes for number concentration are only calculated
-!--             for cloud_scheme = 'morrison'
-                IF ( microphysics_morrison ) THEN
-                   IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
-                      sed_nc(k) = sed_qc_const *                               &
-                              ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *  &
-                              ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
-                   ELSE
-                      sed_nc(k) = 0.0_wp
-                   ENDIF
-
-                   sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *           &
-                                    nc(k,j,i) / dt_micro + sed_nc(k+1)         &
-                                  ) * flag
-
-                   nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *       &
-                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
-                ENDIF
-
-                IF ( qc(k,j,i) > eps_sb .AND.  nc_sedi > eps_mr )  THEN
-                   sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *  &
-                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * &
-                                                                           flag
-                ELSE
-                   sed_qc(k) = 0.0_wp
-                ENDIF
-
-                sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
-                                            dt_micro + sed_qc(k+1)             &
-                               ) * flag
-
-                q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) *          &
-                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
-                qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) *          &
-                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
-                pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) *          &
-                                        ddzu(k+1) / hyrho(k) * lv_d_cp *       &
-                                        d_exner(k) * dt_micro            * flag
-
-!
-!--             Compute the precipitation rate due to cloud (fog) droplets
-                IF ( call_microphysics_at_all_substeps )  THEN
-                   prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)             &
-                                * weight_substep(intermediate_timestep_count)  &
-                                * flag
-                ELSE
-                   prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
-                ENDIF
-
-             ENDDO
-          ENDDO
-       ENDDO
-
-       CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
-
-    END SUBROUTINE sedimentation_cloud
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Computation of sedimentation flux. Implementation according to Stevens
-!> and Seifert (2008). Code is based on UCLA-LES.
-!------------------------------------------------------------------------------!
-    SUBROUTINE sedimentation_rain
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i             !< running index x direction
-       INTEGER(iwp) ::  j             !< running index y direction
-       INTEGER(iwp) ::  k             !< running index z direction
-       INTEGER(iwp) ::  k_run         !<
-       INTEGER(iwp) ::  m             !< running index surface elements
-       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
-       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
-
-       REAL(wp)     ::  c_run                      !<
-       REAL(wp)     ::  d_max                      !<
-       REAL(wp)     ::  d_mean                     !<
-       REAL(wp)     ::  d_min                      !<
-       REAL(wp)     ::  dr                         !<
-       REAL(wp)     ::  flux                       !<
-       REAL(wp)     ::  flag                       !< flag to mask topography grid points
-       REAL(wp)     ::  lambda_r                   !<
-       REAL(wp)     ::  mu_r                       !<
-       REAL(wp)     ::  z_run                      !<
-
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
-       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
-
-       CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
-
-!
-!--    Compute velocities
-       DO  i = nxlg, nxrg
-          DO  j = nysg, nyng
-             DO  k = nzb+1, nzt
-!
-!--             Predetermine flag to mask topography
-                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-
-                IF ( qr(k,j,i) > eps_sb )  THEN
-!
-!--                Weight averaged diameter of rain drops:
-                   dr = ( hyrho(k) * qr(k,j,i) /                               &
-                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
-!
-!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
-!--                Stevens and Seifert, 2008):
-                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *                &
-                                                     ( dr - 1.4E-3_wp ) ) )
-!
-!--                Slope parameter of gamma distribution (Seifert, 2008):
-                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
-                                ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
-
-                   w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
-                                               a_term - b_term * ( 1.0_wp +    &
-                                                  c_term /                     &
-                                                  lambda_r )**( -1.0_wp *      &
-                                                  ( mu_r + 1.0_wp ) )          &
-                                              )                                &
-                                ) * flag
-
-                   w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
-                                               a_term - b_term * ( 1.0_wp +    &
-                                                  c_term /                     &
-                                                  lambda_r )**( -1.0_wp *      &
-                                                  ( mu_r + 4.0_wp ) )          &
-                                             )                                 &
-                                ) * flag
-                ELSE
-                   w_nr(k) = 0.0_wp
-                   w_qr(k) = 0.0_wp
-                ENDIF
-             ENDDO
-!
-!--          Adjust boundary values using surface data type.
-!--          Upward-facing
-             surf_s = bc_h(0)%start_index(j,i)
-             surf_e = bc_h(0)%end_index(j,i)
-             DO  m = surf_s, surf_e
-                k         = bc_h(0)%k(m)
-                w_nr(k-1) = w_nr(k)
-                w_qr(k-1) = w_qr(k)
-             ENDDO
-!
-!--          Downward-facing
-             surf_s = bc_h(1)%start_index(j,i)
-             surf_e = bc_h(1)%end_index(j,i)
-             DO  m = surf_s, surf_e
-                k         = bc_h(1)%k(m)
-                w_nr(k+1) = w_nr(k)
-                w_qr(k+1) = w_qr(k)
-             ENDDO
-!
-!--          Model top boundary value
-             w_nr(nzt+1) = 0.0_wp
-             w_qr(nzt+1) = 0.0_wp
-!
-!--          Compute Courant number
-             DO  k = nzb+1, nzt
-!
-!--             Predetermine flag to mask topography
-                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-
-                c_nr(k) = 0.25_wp * ( w_nr(k-1) +                              &
-                                      2.0_wp * w_nr(k) + w_nr(k+1) ) *         &
-                          dt_micro * ddzu(k) * flag
-                c_qr(k) = 0.25_wp * ( w_qr(k-1) +                              &
-                                      2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
-                          dt_micro * ddzu(k) * flag
-             ENDDO
-!
-!--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
-             IF ( limiter_sedimentation )  THEN
-
-                DO k = nzb+1, nzt
-!
-!--                Predetermine flag to mask topography
-                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-
-                   d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
-                   d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
-                   d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
-
-                   qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
-                                                              2.0_wp * d_max,  &
-                                                              ABS( d_mean ) )  &
-                                                      * flag
-
-                   d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
-                   d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
-                   d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
-
-                   nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
-                                                              2.0_wp * d_max,  &
-                                                              ABS( d_mean ) )
-                ENDDO
-
-             ELSE
-
-                nr_slope = 0.0_wp
-                qr_slope = 0.0_wp
-
-             ENDIF
-
-             sed_nr(nzt+1) = 0.0_wp
-             sed_qr(nzt+1) = 0.0_wp
-!
-!--          Compute sedimentation flux
-             DO  k = nzt, nzb+1, -1
-!
-!--             Predetermine flag to mask topography
-                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-!
-!--             Sum up all rain drop number densities which contribute to the flux
-!--             through k-1/2
-                flux  = 0.0_wp
-                z_run = 0.0_wp ! height above z(k)
-                k_run = k
-                c_run = MIN( 1.0_wp, c_nr(k) )
-                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
-                   flux  = flux + hyrho(k_run) *                               &
-                           ( nr(k_run,j,i) + nr_slope(k_run) *                 &
-                           ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)  &
-                                              * flag
-                   z_run = z_run + dzu(k_run) * flag
-                   k_run = k_run + 1          * flag
-                   c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )    &
-                                              * flag
-                ENDDO
-!
-!--             It is not allowed to sediment more rain drop number density than
-!--             available
-                flux = MIN( flux,                                              &
-                            hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) *    &
-                            dt_micro                                           &
-                          )
-
-                sed_nr(k) = flux / dt_micro * flag
-                nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *          &
-                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
-!
-!--             Sum up all rain water content which contributes to the flux
-!--             through k-1/2
-                flux  = 0.0_wp
-                z_run = 0.0_wp ! height above z(k)
-                k_run = k
-                c_run = MIN( 1.0_wp, c_qr(k) )
-
-                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
-
-                   flux  = flux + hyrho(k_run) * ( qr(k_run,j,i) +             &
-                                  qr_slope(k_run) * ( 1.0_wp - c_run ) *       &
-                                  0.5_wp ) * c_run * dzu(k_run) * flag
-                   z_run = z_run + dzu(k_run)                   * flag
-                   k_run = k_run + 1                            * flag
-                   c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )    &
-                                                                * flag
-
-                ENDDO
-!
-!--             It is not allowed to sediment more rain water content than
-!--             available
-                flux = MIN( flux,                                              &
-                            hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) *      &
-                            dt_micro                                           &
-                          )
-
-                sed_qr(k) = flux / dt_micro * flag
-
-                qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *          &
-                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
-                q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
-                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
-                pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
-                                        ddzu(k+1) / hyrho(k) * lv_d_cp *       &
-                                        d_exner(k) * dt_micro            * flag
-!
-!--             Compute the rain rate
-                IF ( call_microphysics_at_all_substeps )  THEN
-                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)             &
-                                * weight_substep(intermediate_timestep_count) &
-                                * flag
-                ELSE
-                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
-                ENDIF
-
-             ENDDO
-          ENDDO
-       ENDDO
-
-       CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
-
-    END SUBROUTINE sedimentation_rain
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Computation of the precipitation amount due to gravitational settling of
-!> rain and cloud (fog) droplets
-!------------------------------------------------------------------------------!
-    SUBROUTINE calc_precipitation_amount
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i             !< running index x direction
-       INTEGER(iwp) ::  j             !< running index y direction
-       INTEGER(iwp) ::  k             !< running index y direction
-       INTEGER(iwp) ::  m             !< running index surface elements
-
-       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
-            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
-            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
-       THEN
-!
-!--       Run over all upward-facing surface elements, i.e. non-natural,
-!--       natural and urban
-          DO  m = 1, bc_h(0)%ns
-             i = bc_h(0)%i(m)
-             j = bc_h(0)%j(m)
-             k = bc_h(0)%k(m)
-             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
-                                               prr(k,j,i) * hyrho(k) * dt_3d
-          ENDDO
-
-       ENDIF
-
-    END SUBROUTINE calc_precipitation_amount
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Control of microphysics for grid points i,j
-!------------------------------------------------------------------------------!
-
-    SUBROUTINE bcm_actions_ij( i, j )
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i                 !<
-       INTEGER(iwp) ::  j                 !<
-
-       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
-!
-!--       Calculate vertical profile of the hydrostatic pressure (hyp)
-          hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
-          d_exner = exner_function_invers(hyp)
-          exner = 1.0_wp / exner_function_invers(hyp)
-          hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
-!
-!--       Compute reference density
-          rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
-       ENDIF
-
-!
-!--    Compute length of time step
-       IF ( call_microphysics_at_all_substeps )  THEN
-          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
-       ELSE
-          dt_micro = dt_3d
-       ENDIF
-!
-!--    Reset precipitation rate
-       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
-
-!
-!--    Compute cloud physics
-       IF( microphysics_kessler )  THEN
-
-          CALL autoconversion_kessler( i,j )
-          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
-
-       ELSEIF ( microphysics_seifert )  THEN
-
-          CALL adjust_cloud( i,j )
-          IF ( microphysics_morrison ) CALL activation( i,j )
-          IF ( microphysics_morrison ) CALL condensation( i,j )
-          CALL autoconversion( i,j )
-          CALL accretion( i,j )
-          CALL selfcollection_breakup( i,j )
-          CALL evaporation_rain( i,j )
-          CALL sedimentation_rain( i,j )
-          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
-
-       ENDIF
-
-       CALL calc_precipitation_amount( i,j )
-
-    END SUBROUTINE bcm_actions_ij
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Adjust number of raindrops to avoid nonlinear effects in
-!> sedimentation and evaporation of rain drops due to too small or
-!> too big weights of rain drops (Stevens and Seifert, 2008).
-!> The same procedure is applied to cloud droplets if they are determined
-!> prognostically. Call for grid point i,j
-!------------------------------------------------------------------------------!
-    SUBROUTINE adjust_cloud_ij( i, j )
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i                 !<
-       INTEGER(iwp) ::  j                 !<
-       INTEGER(iwp) ::  k                 !<
-
-       REAL(wp) ::  flag                  !< flag to indicate first grid level above surface
-
-       DO  k = nzb+1, nzt
-!
-!--       Predetermine flag to mask topography
-          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-
-          IF ( qr(k,j,i) <= eps_sb )  THEN
-             qr(k,j,i) = 0.0_wp
-             nr(k,j,i) = 0.0_wp
-          ELSE
-!
-!--          Adjust number of raindrops to avoid nonlinear effects in
-!--          sedimentation and evaporation of rain drops due to too small or
-!--          too big weights of rain drops (Stevens and Seifert, 2008).
-             IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
-                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
-             ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
-                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
-             ENDIF
-
-          ENDIF
-
-          IF ( microphysics_morrison ) THEN
-             IF ( qc(k,j,i) <= eps_sb )  THEN
-                qc(k,j,i) = 0.0_wp
-                nc(k,j,i) = 0.0_wp
-             ELSE
-                IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
-                   nc(k,j,i) = qc(k,j,i) * hyrho(k) / xamin * flag
-                ENDIF
-             ENDIF
-          ENDIF
-
-       ENDDO
-
-    END SUBROUTINE adjust_cloud_ij
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Calculate number of activated condensation nucleii after simple activation
-!> scheme of Twomey, 1959.
-!------------------------------------------------------------------------------!
-    SUBROUTINE activation_ij( i, j )
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i                 !<
-       INTEGER(iwp) ::  j                 !<
-       INTEGER(iwp) ::  k                 !<
-
-       REAL(wp)     ::  activ             !<
-       REAL(wp)     ::  afactor           !<
-       REAL(wp)     ::  beta_act          !<
-       REAL(wp)     ::  bfactor           !<
-       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
-       REAL(wp)     ::  k_act             !<
-       REAL(wp)     ::  n_act             !<
-       REAL(wp)     ::  n_ccn             !<
-       REAL(wp)     ::  s_0               !<
-       REAL(wp)     ::  sat_max           !<
-       REAL(wp)     ::  sigma             !<
-       REAL(wp)     ::  sigma_act         !<
-
-       DO  k = nzb+1, nzt
-!
-!--       Predetermine flag to mask topography
-          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-!
-!--       Call calculation of supersaturation  
-          CALL supersaturation ( i, j, k )
-!
-!--       Prescribe parameters for activation
-!--       (see: Bott + Trautmann, 2002, Atm. Res., 64)
-          k_act  = 0.7_wp
-          activ  = 0.0_wp
-
-          IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk )  THEN
-!
-!--          Compute the number of activated Aerosols
-!--          (see: Twomey, 1959, Pure and applied Geophysics, 43)
-             n_act     = na_init * sat**k_act
-!
-!--          Compute the number of cloud droplets
-!--          (see: Morrison + Grabowski, 2007, JAS, 64)
-!            activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro
-
-!
-!--          Compute activation rate after Khairoutdinov and Kogan
-!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
-             sat_max = 0.8_wp / 100.0_wp
-             activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN             &
-                 ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /         &
-                  dt_micro
-
-             nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)
-          ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk )  THEN
-!
-!--          Curvature effect (afactor) with surface tension
-!--          parameterization by Straka (2009)
-             sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
-             afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
-!
-!--          Solute effect (bfactor)
-             bfactor = vanthoff * molecular_weight_of_water *                  &
-                  rho_s / ( molecular_weight_of_solute * rho_l )
-
-!
-!--          Prescribe power index that describes the soluble fraction
-!--          of an aerosol particle (beta).
-!--          (see: Morrison + Grabowski, 2007, JAS, 64)
-             beta_act  = 0.5_wp
-             sigma_act = sigma_bulk**( 1.0_wp + beta_act )
-!
-!--          Calculate mean geometric supersaturation (s_0) with
-!--          parameterization by Khvorostyanov and Curry (2006)
-             s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *         &
-               ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
-
-!
-!--          Calculate number of activated CCN as a function of
-!--          supersaturation and taking Koehler theory into account
-!--          (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
-             n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                    &
-                LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
-             activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
-
-             nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro * flag), na_init)
-          ENDIF
-
-       ENDDO
-
-    END SUBROUTINE activation_ij
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Calculate condensation rate for cloud water content (after Khairoutdinov and
-!> Kogan, 2000).
-!------------------------------------------------------------------------------!
-    SUBROUTINE condensation_ij( i, j )
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i                 !<
-       INTEGER(iwp) ::  j                 !<
-       INTEGER(iwp) ::  k                 !<
-
-       REAL(wp)     ::  cond              !<
-       REAL(wp)     ::  cond_max          !<
-       REAL(wp)     ::  dc                !<
-       REAL(wp)     ::  evap              !<
-       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
-       REAL(wp)     ::  g_fac             !<
-       REAL(wp)     ::  nc_0              !<
-       REAL(wp)     ::  temp              !<
-       REAL(wp)     ::  xc                !<
-
-
-       DO  k = nzb+1, nzt
-!
-!--       Predetermine flag to mask topography
-          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-!
-!--       Call calculation of supersaturation  
-          CALL supersaturation ( i, j, k )
-!
-!--       Actual temperature:
-          temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
-
-          g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
-                              l_v / ( thermal_conductivity_l * temp )    &
-                              + r_v * temp / ( diff_coeff_l * e_s )      &
-                            )
-!
-!--       Mean weight of cloud drops
-          IF ( nc(k,j,i) <= 0.0_wp) CYCLE
-          xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
-!
-!--       Weight averaged diameter of cloud drops:
-          dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
-!
-!--       Integral diameter of cloud drops
-          nc_0 = nc(k,j,i) * dc
-!
-!--       Condensation needs only to be calculated in supersaturated regions
-          IF ( sat > 0.0_wp )  THEN
-!
-!--          Condensation rate of cloud water content
-!--          after KK scheme.
-!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
-             cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
-             cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
-             cond      = MIN( cond, cond_max / dt_micro )
-
-             qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag
-          ELSEIF ( sat < 0.0_wp ) THEN
-             evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
-             evap      = MAX( evap, -qc(k,j,i) / dt_micro )
-
-             qc(k,j,i) = qc(k,j,i) + evap * dt_micro * flag
-          ENDIF
-       ENDDO
-
-    END SUBROUTINE condensation_ij
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
-!------------------------------------------------------------------------------!
-    SUBROUTINE autoconversion_ij( i, j )
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i                 !<
-       INTEGER(iwp) ::  j                 !<
-       INTEGER(iwp) ::  k                 !<
-
-       REAL(wp)     ::  alpha_cc          !<
-       REAL(wp)     ::  autocon           !<
-       REAL(wp)     ::  dissipation       !<
-       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
-       REAL(wp)     ::  k_au              !<
-       REAL(wp)     ::  l_mix             !<
-       REAL(wp)     ::  nc_auto           !<
-       REAL(wp)     ::  nu_c              !<
-       REAL(wp)     ::  phi_au            !<
-       REAL(wp)     ::  r_cc              !<
-       REAL(wp)     ::  rc                !<
-       REAL(wp)     ::  re_lambda         !<
-       REAL(wp)     ::  sigma_cc          !<
-       REAL(wp)     ::  tau_cloud         !<
-       REAL(wp)     ::  xc                !<
-
-       DO  k = nzb+1, nzt
-!
-!--       Predetermine flag to mask topography
-          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-          IF ( microphysics_morrison ) THEN
-             nc_auto = nc(k,j,i)
-          ELSE
-             nc_auto = nc_const
-          ENDIF
-
-          IF ( qc(k,j,i) > eps_sb  .AND.  nc_auto > eps_mr )  THEN
-
-             k_au = k_cc / ( 20.0_wp * x0 )
-!
-!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
-!--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
-             tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) ),     &
-                              0.0_wp )
-!
-!--          Universal function for autoconversion process
-!--          (Seifert and Beheng, 2006):
-             phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
-!
-!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
-!--          (Use constant nu_c = 1.0_wp instead?)
-             nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
-!
-!--          Mean weight of cloud droplets:
-             xc = hyrho(k) * qc(k,j,i) / nc_auto
-!
-!--          Parameterized turbulence effects on autoconversion (Seifert,
-!--          Nuijens and Stevens, 2010)
-             IF ( collision_turbulence )  THEN
-!
-!--             Weight averaged radius of cloud droplets:
-                rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
-
-                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
-                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
-                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
-!
-!--             Mixing length (neglecting distance to ground and stratification)
-                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
-!
-!--             Limit dissipation rate according to Seifert, Nuijens and
-!--             Stevens (2010)
-                dissipation = MIN( 0.06_wp, diss(k,j,i) )
-!
-!--             Compute Taylor-microscale Reynolds number:
-                re_lambda = 6.0_wp / 11.0_wp *                                 &
-                            ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *         &
-                            SQRT( 15.0_wp / kin_vis_air ) *                    &
-                            dissipation**( 1.0_wp / 6.0_wp )
-!
-!--             The factor of 1.0E4 is needed to convert the dissipation rate
-!--             from m2 s-3 to cm2 s-3.
-                k_au = k_au * ( 1.0_wp +                                       &
-                                dissipation * 1.0E4_wp *                       &
-                                ( re_lambda * 1.0E-3_wp )**0.25_wp *           &
-                                ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /  &
-                                                  sigma_cc )**2                &
-                                                ) + beta_cc                    &
-                                )                                              &
-                              )
-             ENDIF
-!
-!--          Autoconversion rate (Seifert and Beheng, 2006):
-             autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
-                       ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *            &
-                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
-                       rho_surface
-             autocon = MIN( autocon, qc(k,j,i) / dt_micro )
-
-             qr(k,j,i) = qr(k,j,i) + autocon * dt_micro                 * flag
-             qc(k,j,i) = qc(k,j,i) - autocon * dt_micro                 * flag
-             nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro * flag
-             IF ( microphysics_morrison ) THEN
-                nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *                &
-                                    autocon / x0 * hyrho(k) * dt_micro * flag )
-             ENDIF
-
-          ENDIF
-
-       ENDDO
-
-    END SUBROUTINE autoconversion_ij
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Autoconversion process (Kessler, 1969).
-!------------------------------------------------------------------------------!
-    SUBROUTINE autoconversion_kessler_ij( i, j )
-
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i      !<
-       INTEGER(iwp) ::  j      !<
-       INTEGER(iwp) ::  k      !<
-       INTEGER(iwp) ::  k_wall !< topography top index
-
-       REAL(wp)    ::  dqdt_precip !<
-       REAL(wp)    ::  flag              !< flag to indicate first grid level above surface
-
-!
-!--    Determine vertical index of topography top
-       k_wall = get_topography_top_index_ji( j, i, 's' )
-       DO  k = nzb+1, nzt
-!
-!--       Predetermine flag to mask topography
-          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-
-          IF ( qc(k,j,i) > ql_crit )  THEN
-             dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
-          ELSE
-             dqdt_precip = 0.0_wp
-          ENDIF
-
-          qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
-          q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro * flag
-          pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp * d_exner(k)  &
-                                                                 * flag
-
-!
-!--       Compute the rain rate (stored on surface grid point)
-          prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
-
-       ENDDO
-
-    END SUBROUTINE autoconversion_kessler_ij
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
-!------------------------------------------------------------------------------!
-    SUBROUTINE accretion_ij( i, j )
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i                 !<
-       INTEGER(iwp) ::  j                 !<
-       INTEGER(iwp) ::  k                 !<
-
-       REAL(wp)     ::  accr              !<
-       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
-       REAL(wp)     ::  k_cr              !<
-       REAL(wp)     ::  nc_accr           !<
-       REAL(wp)     ::  phi_ac            !<
-       REAL(wp)     ::  tau_cloud         !<
-       REAL(wp)     ::  xc                !<
-
-
-       DO  k = nzb+1, nzt
-!
-!--       Predetermine flag to mask topography
-          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-          IF ( microphysics_morrison ) THEN
-             nc_accr = nc(k,j,i)
-          ELSE
-             nc_accr = nc_const
-          ENDIF
-
-          IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )  .AND.  &
-               ( nc_accr > eps_mr ) )  THEN
-!
-!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
-             tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
-!
-!--          Universal function for accretion process
-!--          (Seifert and Beheng, 2001):
-             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
-
-!
-!--          Mean weight of cloud drops
-             xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
-!
-!--          Parameterized turbulence effects on autoconversion (Seifert,
-!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
-!--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
-             IF ( collision_turbulence )  THEN
-                k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                            &
-                                 MIN( 600.0_wp,                                &
-                                      diss(k,j,i) * 1.0E4_wp )**0.25_wp        &
-                               )
-             ELSE
-                k_cr = k_cr0
-             ENDIF
-!
-!--          Accretion rate (Seifert and Beheng, 2006):
-             accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *                      &
-                    SQRT( rho_surface * hyrho(k) )
-             accr = MIN( accr, qc(k,j,i) / dt_micro )
-
-             qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
-             qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
-             IF ( microphysics_morrison )  THEN
-                nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc *               &
-                                             hyrho(k) * dt_micro * flag        &
-                                         )
-             ENDIF
-
-
-          ENDIF
-
-       ENDDO
-
-    END SUBROUTINE accretion_ij
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
-!------------------------------------------------------------------------------!
-    SUBROUTINE selfcollection_breakup_ij( i, j )
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i                 !<
-       INTEGER(iwp) ::  j                 !<
-       INTEGER(iwp) ::  k                 !<
-
-       REAL(wp)     ::  breakup           !<
-       REAL(wp)     ::  dr                !<
-       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
-       REAL(wp)     ::  phi_br            !<
-       REAL(wp)     ::  selfcoll          !<
-
-       DO  k = nzb+1, nzt
-!
-!--       Predetermine flag to mask topography
-          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-
-          IF ( qr(k,j,i) > eps_sb )  THEN
-!
-!--          Selfcollection rate (Seifert and Beheng, 2001):
-             selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * SQRT( hyrho(k) * rho_surface )
-!
-!--          Weight averaged diameter of rain drops:
-             dr = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
-!
-!--          Collisional breakup rate (Seifert, 2008):
-             IF ( dr >= 0.3E-3_wp )  THEN
-                phi_br  = k_br * ( dr - 1.1E-3_wp )
-                breakup = selfcoll * ( phi_br + 1.0_wp )
-             ELSE
-                breakup = 0.0_wp
-             ENDIF
-
-             selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
-             nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
-
-          ENDIF
-       ENDDO
-
-    END SUBROUTINE selfcollection_breakup_ij
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
 !> Evaporation of precipitable water. Condensation is neglected for
 !> precipitable water. Call for grid point i,j
@@ -3739,5 +3304,5 @@
           IF ( qr(k,j,i) > eps_sb )  THEN
 !
-!--          Call calculation of supersaturation  
+!--          Call calculation of supersaturation
              CALL supersaturation ( i, j, k )
 !
@@ -3822,4 +3387,102 @@
 ! ------------
 !> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
+!------------------------------------------------------------------------------!
+    SUBROUTINE sedimentation_cloud
+
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i             !<
+       INTEGER(iwp) ::  j             !<
+       INTEGER(iwp) ::  k             !<
+
+       REAL(wp) ::  flag              !< flag to mask topography grid points
+       REAL(wp) ::  nc_sedi           !<
+
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc !<
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc !<
+
+
+       CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
+
+       sed_qc(nzt+1) = 0.0_wp
+       sed_nc(nzt+1) = 0.0_wp
+
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             DO  k = nzt, nzb+1, -1
+!
+!--             Predetermine flag to mask topography
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+
+                IF ( microphysics_morrison ) THEN
+                   nc_sedi = nc(k,j,i)
+                ELSE
+                   nc_sedi = nc_const
+                ENDIF
+
+!
+!--             Sedimentation fluxes for number concentration are only calculated
+!--             for cloud_scheme = 'morrison'
+                IF ( microphysics_morrison ) THEN
+                   IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
+                      sed_nc(k) = sed_qc_const *                               &
+                              ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *  &
+                              ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
+                   ELSE
+                      sed_nc(k) = 0.0_wp
+                   ENDIF
+
+                   sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *           &
+                                    nc(k,j,i) / dt_micro + sed_nc(k+1)         &
+                                  ) * flag
+
+                   nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *       &
+                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
+                ENDIF
+
+                IF ( qc(k,j,i) > eps_sb .AND.  nc_sedi > eps_mr )  THEN
+                   sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *  &
+                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * &
+                                                                           flag
+                ELSE
+                   sed_qc(k) = 0.0_wp
+                ENDIF
+
+                sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
+                                            dt_micro + sed_qc(k+1)             &
+                               ) * flag
+
+                q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) *          &
+                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
+                qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) *          &
+                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
+                pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) *          &
+                                        ddzu(k+1) / hyrho(k) * lv_d_cp *       &
+                                        d_exner(k) * dt_micro            * flag
+
+!
+!--             Compute the precipitation rate due to cloud (fog) droplets
+                IF ( call_microphysics_at_all_substeps )  THEN
+                   prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)             &
+                                * weight_substep(intermediate_timestep_count)  &
+                                * flag
+                ELSE
+                   prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
+                ENDIF
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+       CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
+
+    END SUBROUTINE sedimentation_cloud
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
 !> Call for grid point i,j
 !------------------------------------------------------------------------------!
@@ -3902,4 +3565,250 @@
 
     END SUBROUTINE sedimentation_cloud_ij
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Computation of sedimentation flux. Implementation according to Stevens
+!> and Seifert (2008). Code is based on UCLA-LES.
+!------------------------------------------------------------------------------!
+    SUBROUTINE sedimentation_rain
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i             !< running index x direction
+       INTEGER(iwp) ::  j             !< running index y direction
+       INTEGER(iwp) ::  k             !< running index z direction
+       INTEGER(iwp) ::  k_run         !<
+       INTEGER(iwp) ::  m             !< running index surface elements
+       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
+       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
+
+       REAL(wp)     ::  c_run                      !<
+       REAL(wp)     ::  d_max                      !<
+       REAL(wp)     ::  d_mean                     !<
+       REAL(wp)     ::  d_min                      !<
+       REAL(wp)     ::  dr                         !<
+       REAL(wp)     ::  flux                       !<
+       REAL(wp)     ::  flag                       !< flag to mask topography grid points
+       REAL(wp)     ::  lambda_r                   !<
+       REAL(wp)     ::  mu_r                       !<
+       REAL(wp)     ::  z_run                      !<
+
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
+
+       CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
+
+!
+!--    Compute velocities
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             DO  k = nzb+1, nzt
+!
+!--             Predetermine flag to mask topography
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+
+                IF ( qr(k,j,i) > eps_sb )  THEN
+!
+!--                Weight averaged diameter of rain drops:
+                   dr = ( hyrho(k) * qr(k,j,i) /                               &
+                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
+!
+!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
+!--                Stevens and Seifert, 2008):
+                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *                &
+                                                     ( dr - 1.4E-3_wp ) ) )
+!
+!--                Slope parameter of gamma distribution (Seifert, 2008):
+                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
+                                ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
+
+                   w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
+                                               a_term - b_term * ( 1.0_wp +    &
+                                                  c_term /                     &
+                                                  lambda_r )**( -1.0_wp *      &
+                                                  ( mu_r + 1.0_wp ) )          &
+                                              )                                &
+                                ) * flag
+
+                   w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
+                                               a_term - b_term * ( 1.0_wp +    &
+                                                  c_term /                     &
+                                                  lambda_r )**( -1.0_wp *      &
+                                                  ( mu_r + 4.0_wp ) )          &
+                                             )                                 &
+                                ) * flag
+                ELSE
+                   w_nr(k) = 0.0_wp
+                   w_qr(k) = 0.0_wp
+                ENDIF
+             ENDDO
+!
+!--          Adjust boundary values using surface data type.
+!--          Upward-facing
+             surf_s = bc_h(0)%start_index(j,i)
+             surf_e = bc_h(0)%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k         = bc_h(0)%k(m)
+                w_nr(k-1) = w_nr(k)
+                w_qr(k-1) = w_qr(k)
+             ENDDO
+!
+!--          Downward-facing
+             surf_s = bc_h(1)%start_index(j,i)
+             surf_e = bc_h(1)%end_index(j,i)
+             DO  m = surf_s, surf_e
+                k         = bc_h(1)%k(m)
+                w_nr(k+1) = w_nr(k)
+                w_qr(k+1) = w_qr(k)
+             ENDDO
+!
+!--          Model top boundary value
+             w_nr(nzt+1) = 0.0_wp
+             w_qr(nzt+1) = 0.0_wp
+!
+!--          Compute Courant number
+             DO  k = nzb+1, nzt
+!
+!--             Predetermine flag to mask topography
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+
+                c_nr(k) = 0.25_wp * ( w_nr(k-1) +                              &
+                                      2.0_wp * w_nr(k) + w_nr(k+1) ) *         &
+                          dt_micro * ddzu(k) * flag
+                c_qr(k) = 0.25_wp * ( w_qr(k-1) +                              &
+                                      2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
+                          dt_micro * ddzu(k) * flag
+             ENDDO
+!
+!--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
+             IF ( limiter_sedimentation )  THEN
+
+                DO k = nzb+1, nzt
+!
+!--                Predetermine flag to mask topography
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+
+                   d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
+                   d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
+                   d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
+
+                   qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
+                                                              2.0_wp * d_max,  &
+                                                              ABS( d_mean ) )  &
+                                                      * flag
+
+                   d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
+                   d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
+                   d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
+
+                   nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
+                                                              2.0_wp * d_max,  &
+                                                              ABS( d_mean ) )
+                ENDDO
+
+             ELSE
+
+                nr_slope = 0.0_wp
+                qr_slope = 0.0_wp
+
+             ENDIF
+
+             sed_nr(nzt+1) = 0.0_wp
+             sed_qr(nzt+1) = 0.0_wp
+!
+!--          Compute sedimentation flux
+             DO  k = nzt, nzb+1, -1
+!
+!--             Predetermine flag to mask topography
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+!
+!--             Sum up all rain drop number densities which contribute to the flux
+!--             through k-1/2
+                flux  = 0.0_wp
+                z_run = 0.0_wp ! height above z(k)
+                k_run = k
+                c_run = MIN( 1.0_wp, c_nr(k) )
+                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
+                   flux  = flux + hyrho(k_run) *                               &
+                           ( nr(k_run,j,i) + nr_slope(k_run) *                 &
+                           ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)  &
+                                              * flag
+                   z_run = z_run + dzu(k_run) * flag
+                   k_run = k_run + 1          * flag
+                   c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )    &
+                                              * flag
+                ENDDO
+!
+!--             It is not allowed to sediment more rain drop number density than
+!--             available
+                flux = MIN( flux,                                              &
+                            hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) *    &
+                            dt_micro                                           &
+                          )
+
+                sed_nr(k) = flux / dt_micro * flag
+                nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *          &
+                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
+!
+!--             Sum up all rain water content which contributes to the flux
+!--             through k-1/2
+                flux  = 0.0_wp
+                z_run = 0.0_wp ! height above z(k)
+                k_run = k
+                c_run = MIN( 1.0_wp, c_qr(k) )
+
+                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
+
+                   flux  = flux + hyrho(k_run) * ( qr(k_run,j,i) +             &
+                                  qr_slope(k_run) * ( 1.0_wp - c_run ) *       &
+                                  0.5_wp ) * c_run * dzu(k_run) * flag
+                   z_run = z_run + dzu(k_run)                   * flag
+                   k_run = k_run + 1                            * flag
+                   c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )    &
+                                                                * flag
+
+                ENDDO
+!
+!--             It is not allowed to sediment more rain water content than
+!--             available
+                flux = MIN( flux,                                              &
+                            hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) *      &
+                            dt_micro                                           &
+                          )
+
+                sed_qr(k) = flux / dt_micro * flag
+
+                qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *          &
+                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
+                q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
+                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
+                pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
+                                        ddzu(k+1) / hyrho(k) * lv_d_cp *       &
+                                        d_exner(k) * dt_micro            * flag
+!
+!--             Compute the rain rate
+                IF ( call_microphysics_at_all_substeps )  THEN
+                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)             &
+                                * weight_substep(intermediate_timestep_count) &
+                                * flag
+                ELSE
+                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
+                ENDIF
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+       CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
+
+    END SUBROUTINE sedimentation_rain
 
 
@@ -4123,4 +4032,39 @@
 
     END SUBROUTINE sedimentation_rain_ij
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Computation of the precipitation amount due to gravitational settling of
+!> rain and cloud (fog) droplets
+!------------------------------------------------------------------------------!
+    SUBROUTINE calc_precipitation_amount
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i             !< running index x direction
+       INTEGER(iwp) ::  j             !< running index y direction
+       INTEGER(iwp) ::  k             !< running index y direction
+       INTEGER(iwp) ::  m             !< running index surface elements
+
+       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
+            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
+            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
+       THEN
+!
+!--       Run over all upward-facing surface elements, i.e. non-natural,
+!--       natural and urban
+          DO  m = 1, bc_h(0)%ns
+             i = bc_h(0)%i(m)
+             j = bc_h(0)%j(m)
+             k = bc_h(0)%k(m)
+             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
+                                               prr(k,j,i) * hyrho(k) * dt_3d
+          ENDDO
+
+       ENDIF
+
+    END SUBROUTINE calc_precipitation_amount
 
 
