Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 1360)
+++ /palm/trunk/SOURCE/Makefile	(revision 1361)
@@ -20,5 +20,5 @@
 # Current revisions:
 # ------------------
-# 
+# cpulog added to microphysics
 # 
 # Former revisions:
@@ -341,5 +341,5 @@
 ls_forcing.o: modules.o cpulog.o mod_kinds.o
 message.o: modules.o mod_kinds.o
-microphysics.o: modules.o mod_kinds.o
+microphysics.o: modules.o cpulog.o mod_kinds.o
 modules.o: modules.f90 mod_kinds.o
 mod_kinds.o: mod_kinds.f90
Index: /palm/trunk/SOURCE/advec_s_bc.f90
===================================================================
--- /palm/trunk/SOURCE/advec_s_bc.f90	(revision 1360)
+++ /palm/trunk/SOURCE/advec_s_bc.f90	(revision 1361)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! nr and qr added
 ! 
 ! Former revisions:
@@ -926,4 +926,48 @@
        ENDIF
 
+    ELSEIF ( sk_char == 'qr' )  THEN
+
+!
+!--    Rain water content boundary condition at the bottom boundary:
+!--    Dirichlet (fixed surface rain water content).
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
+             sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
+             sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
+          ENDDO
+       ENDDO
+
+!
+!--    Rain water content boundary condition at the top boundary: Dirichlet
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
+             sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
+          ENDDO
+       ENDDO
+
+    ELSEIF ( sk_char == 'nr' )  THEN
+
+!
+!--    Rain drop concentration boundary condition at the bottom boundary:
+!--    Dirichlet (fixed surface rain drop concentration).
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
+             sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
+             sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
+          ENDDO
+       ENDDO
+
+!
+!--    Rain drop concentration boundary condition at the top boundary: Dirichlet
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
+             sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
+          ENDDO
+       ENDDO
+
     ELSEIF ( sk_char == 'e' )  THEN
 
Index: /palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- /palm/trunk/SOURCE/advec_ws.f90	(revision 1360)
+++ /palm/trunk/SOURCE/advec_ws.f90	(revision 1361)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! accelerator and vector version for qr and nr added
 ! 
 ! Former revisions:
@@ -2233,5 +2233,5 @@
        USE statistics,                                                        &
            ONLY:  sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,           &
-                  weight_substep
+                  sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep
 
        IMPLICIT NONE
@@ -2656,4 +2656,16 @@
                         *   weight_substep(intermediate_timestep_count)
                     ENDDO
+                 CASE ( 'qr' )
+                    DO  k = nzb, nzt
+                       sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn)         &
+                        + ( flux_t(k) + diss_t(k) )                          &
+                        *   weight_substep(intermediate_timestep_count)
+                    ENDDO
+                 CASE ( 'nr' )
+                    DO  k = nzb, nzt
+                       sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn)         &
+                        + ( flux_t(k) + diss_t(k) )                          &
+                        *   weight_substep(intermediate_timestep_count)
+                    ENDDO
 
               END SELECT
@@ -2690,5 +2702,5 @@
 !        USE statistics,                                                       &
 !            ONLY:  sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,          &
-!                   weight_substep
+!                   sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep
 
        IMPLICIT NONE
@@ -2975,5 +2987,13 @@
 !                   CASE ( 'q' )
 !                      sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)           &
-!                      + ( flux_t + diss_t )                                &
+!                      + ( flux_t + diss_t )                                 &
+!                      *   weight_substep(intermediate_timestep_count)
+!                   CASE ( 'qr' )
+!                      sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn)         &
+!                      + ( flux_t + diss_t )                                 &
+!                      *   weight_substep(intermediate_timestep_count)
+!                   CASE ( 'nr' )
+!                      sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn)         &
+!                      + ( flux_t + diss_t )                                 &
 !                      *   weight_substep(intermediate_timestep_count)
 !
Index: /palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- /palm/trunk/SOURCE/boundary_conds.f90	(revision 1360)
+++ /palm/trunk/SOURCE/boundary_conds.f90	(revision 1361)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bottom and top boundary conditions of rain water content (qr) and 
+! rain drop concentration (nr) changed to Dirichlet
 ! 
 ! Former revisions:
@@ -271,18 +272,17 @@
        q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
 
-       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                    &
-            precipitation )  THEN
+       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
 !             
-!--       Surface conditions rain water (Neumann)
+!--       Surface conditions rain water (Dirichlet)
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                qr_p(nzb_s_inner(j,i),j,i) = qr_p(nzb_s_inner(j,i)+1,j,i)
-                nr_p(nzb_s_inner(j,i),j,i) = nr_p(nzb_s_inner(j,i)+1,j,i)
-             ENDDO
-          ENDDO
-!
-!--       Top boundary condition for rain water (Neumann)
-          qr_p(nzt+1,:,:) = qr_p(nzt,:,:)
-          nr_p(nzt+1,:,:) = nr_p(nzt,:,:)
+                qr_p(nzb_s_inner(j,i),j,i) = 0.0_wp
+                nr_p(nzb_s_inner(j,i),j,i) = 0.0_wp
+             ENDDO
+          ENDDO
+!
+!--       Top boundary condition for rain water (Dirichlet)
+          qr_p(nzt+1,:,:) = 0.0_wp
+          nr_p(nzt+1,:,:) = 0.0_wp
            
        ENDIF
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 1360)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 1361)
@@ -20,6 +20,7 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+! PA0363 removed
+! PA0362 changed 
+!
 ! Former revisions:
 ! -----------------
@@ -862,17 +863,11 @@
     ENDIF 
 
-    IF ( loop_optimization /= 'cache' .AND.  cloud_physics  .AND.            &
-         icloud_scheme == 0 ) THEN
+    IF ( .NOT. ( loop_optimization == 'cache'  .OR.                            &
+                 loop_optimization == 'vector' )                               &
+         .AND.  cloud_physics  .AND.  icloud_scheme == 0 )  THEN
        message_string = 'cloud_scheme = seifert_beheng requires ' // &
-                        'loop_optimization = cache'
+                        'loop_optimization = "cache" or "vector"'
        CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 )
     ENDIF 
-
-!    IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
-!         .NOT. precipitation  .AND.  .NOT. drizzle ) THEN
-!       message_string = 'cloud_scheme = seifert_beheng requires ' // &
-!                        'precipitation = .TRUE. or drizzle = .TRUE.'
-!       CALL message( 'check_parameters', 'PA0363', 1, 2, 0, 6, 0 )
-!    ENDIF 
 
 !
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 1360)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 1361)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! tend_* removed
+! Bugfix: w_subs is not allocated anymore if it is already allocated
 ! 
 ! Former revisions:
@@ -366,8 +367,4 @@
                 ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 #endif
-!
-!--             3D-tendency arrays
-                ALLOCATE( tend_pt(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
-                          tend_q(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 
                 IF ( precipitation )  THEN
@@ -375,9 +372,5 @@
 !--                1D-arrays 
                    ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) ) 
-!
-!
-!--                3D-tendency arrays
-                   ALLOCATE( tend_nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
-                             tend_qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
 !
 !--                2D-rain water content and rain drop concentration arrays
@@ -475,7 +468,8 @@
 !
 !-- 1D-array for large scale subsidence velocity
-    ALLOCATE ( w_subs(nzb:nzt+1) )
-    w_subs = 0.0_wp
-
+    IF ( .NOT. ALLOCATED( w_subs ) )  THEN
+       ALLOCATE ( w_subs(nzb:nzt+1) )
+       w_subs = 0.0_wp
+    ENDIF
 
 !
Index: /palm/trunk/SOURCE/init_cloud_physics.f90
===================================================================
--- /palm/trunk/SOURCE/init_cloud_physics.f90	(revision 1360)
+++ /palm/trunk/SOURCE/init_cloud_physics.f90	(revision 1361)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! sed_qc_const is now calculated here (2-moment microphysics)
 ! 
 ! Former revisions:
@@ -73,8 +73,8 @@
     USE cloud_parameters,                                                      &
         ONLY:  bfactor, cp, c_sedimentation, dpirho_l, dt_precipitation,       &
-               hyrho, l_d_cp, l_d_r, l_d_rv, l_v, mass_of_solute,              &
+               hyrho, k_st, l_d_cp, l_d_r, l_d_rv, l_v, mass_of_solute,  &
                molecular_weight_of_solute, molecular_weight_of_water, pirho_l, &
-               pt_d_t, rho_l, r_d, r_v, schmidt, schmidt_p_1d3, t_d_pt,        &
-               vanthoff, w_precipitation
+               pt_d_t, rho_l, r_d, r_v, sed_qc_const, schmidt, schmidt_p_1d3,  &
+               sigma_gc, t_d_pt, vanthoff, w_precipitation
         
     USE constants,                                                             &
@@ -108,5 +108,10 @@
 
     pirho_l  = pi * rho_l / 6.0_wp
-    dpirho_l = 1.0_wp / pirho_l 
+    dpirho_l = 1.0_wp / pirho_l
+!
+!-- constant fot the sedimentation of cloud water (2-moment cloud physics)
+    sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l )                   &
+                          )**( 2.0_wp / 3.0_wp ) *                             &
+                   EXP( 5.0_wp * LOG( sigma_gc )**2 )
 !
 !-- Calculate timestep according to precipitation
Index: /palm/trunk/SOURCE/microphysics.f90
===================================================================
--- /palm/trunk/SOURCE/microphysics.f90	(revision 1360)
+++ /palm/trunk/SOURCE/microphysics.f90	(revision 1361)
@@ -20,5 +20,7 @@
 ! Current revisions:
 ! ------------------
-! 
+! Bugfix in sedimentation_rain: Index corrected.
+! Vectorized version of adjust_cloud added.
+! Little reformatting of the code.
 ! 
 ! Former revisions:
@@ -126,11 +128,23 @@
     SUBROUTINE microphysics_control
 
-       USE arrays_3d
-       USE cloud_parameters
-       USE control_parameters
-       USE grid_variables
-       USE indices
+       USE arrays_3d,                                                          &
+           ONLY:  hyp, nr, pt, pt_init, q, qc, qr, zu
+
+       USE cloud_parameters,                                                   &
+           ONLY:  cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt
+
+       USE control_parameters,                                                 &
+           ONLY:  call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, &
+                  g, intermediate_timestep_count,                              &
+                  large_scale_forcing, lsf_surf, precipitation, pt_surface,    &
+                  rho_surface,surface_pressure
+
+       USE indices,                                                            &
+           ONLY:  nzb, nzt
+
        USE kinds
-       USE statistics
+
+       USE statistics,                                                         &
+           ONLY:  weight_pres
 
        IMPLICIT NONE
@@ -139,231 +153,8 @@
        INTEGER(iwp) ::  j                 !:
        INTEGER(iwp) ::  k                 !:
- 
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb_s_inner(j,i)+1, nzt
-
-             ENDDO
-          ENDDO
-       ENDDO
-
-    END SUBROUTINE microphysics_control
-
-    SUBROUTINE adjust_cloud
-
-       USE arrays_3d
-       USE cloud_parameters
-       USE indices
-       USE kinds
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i                 !:
-       INTEGER(iwp) ::  j                 !:
-       INTEGER(iwp) ::  k                 !:
-
- 
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb_s_inner(j,i)+1, nzt
-
-             ENDDO
-          ENDDO
-       ENDDO
-
-    END SUBROUTINE adjust_cloud
-
-
-    SUBROUTINE autoconversion
-
-       USE arrays_3d
-       USE cloud_parameters
-       USE control_parameters
-       USE grid_variables
-       USE indices
-       USE kinds
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i                 !:
-       INTEGER(iwp) ::  j                 !:
-       INTEGER(iwp) ::  k                 !:
-
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb_s_inner(j,i)+1, nzt
-
-             ENDDO
-          ENDDO
-       ENDDO
-
-    END SUBROUTINE autoconversion
-
-
-    SUBROUTINE accretion
-
-       USE arrays_3d
-       USE cloud_parameters
-       USE control_parameters
-       USE indices
-       USE kinds
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i                 !:
-       INTEGER(iwp) ::  j                 !:
-       INTEGER(iwp) ::  k                 !:
-
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb_s_inner(j,i)+1, nzt
-
-             ENDDO
-          ENDDO
-       ENDDO
-
-    END SUBROUTINE accretion
-
-
-    SUBROUTINE selfcollection_breakup
-
-       USE arrays_3d
-       USE cloud_parameters
-       USE control_parameters
-       USE indices
-       USE kinds
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i                 !:
-       INTEGER(iwp) ::  j                 !:
-       INTEGER(iwp) ::  k                 !:
-
- 
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb_s_inner(j,i)+1, nzt
-
-             ENDDO
-          ENDDO
-       ENDDO
-
-    END SUBROUTINE selfcollection_breakup
-
-
-    SUBROUTINE evaporation_rain
-
-       USE arrays_3d
-       USE cloud_parameters
-       USE constants
-       USE control_parameters
-       USE indices
-       USE kinds
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i                 !:
-       INTEGER(iwp) ::  j                 !:
-       INTEGER(iwp) ::  k                 !:
- 
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb_s_inner(j,i)+1, nzt
-
-             ENDDO
-          ENDDO
-       ENDDO
-
-    END SUBROUTINE evaporation_rain
-
-
-    SUBROUTINE sedimentation_cloud
-
-       USE arrays_3d
-       USE cloud_parameters
-       USE constants
-       USE control_parameters
-       USE indices
-       USE kinds
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i                 !:
-       INTEGER(iwp) ::  j                 !:
-       INTEGER(iwp) ::  k                 !:
- 
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb_s_inner(j,i)+1, nzt
-
-             ENDDO
-          ENDDO
-       ENDDO
-
-    END SUBROUTINE sedimentation_cloud
-
-
-    SUBROUTINE sedimentation_rain
-
-       USE arrays_3d
-       USE cloud_parameters
-       USE constants
-       USE control_parameters
-       USE indices
-       USE kinds
-       USE statistics
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i                 !:
-       INTEGER(iwp) ::  j                 !:
-       INTEGER(iwp) ::  k                 !:
- 
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb_s_inner(j,i)+1, nzt
-
-             ENDDO
-          ENDDO
-       ENDDO
-
-    END SUBROUTINE sedimentation_rain
-
-
-!------------------------------------------------------------------------------!
-! Call for grid point i,j
-!------------------------------------------------------------------------------!
-
-    SUBROUTINE microphysics_control_ij( i, j )
-
-       USE arrays_3d,                                                          &
-           ONLY:  hyp, nc_1d,  nr, nr_1d, pt, pt_init, pt_1d, q, q_1d, qc,     &
-                  qc_1d, qr, qr_1d, tend_nr, tend_pt, tend_q, tend_qr, zu
-
-       USE cloud_parameters,                                                   &
-           ONLY:  cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt
-
-       USE control_parameters,                                                 &
-           ONLY:  drizzle, dt_3d, dt_micro, g, intermediate_timestep_count,    &
-                  large_scale_forcing, lsf_surf, precipitation, pt_surface,    &
-                  rho_surface,surface_pressure
-
-       USE indices,                                                            &
-           ONLY:  nzb, nzt
-
-       USE kinds
-
-       USE statistics,                                                         &
-           ONLY:  weight_pres
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i                 !:
-       INTEGER(iwp) ::  j                 !:
-       INTEGER(iwp) ::  k                 !:
 
        REAL(wp)     ::  t_surface         !:
 
-       IF ( large_scale_forcing .AND. lsf_surf ) THEN
+       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
 !
 !--       Calculate:
@@ -374,5 +165,6 @@
           DO  k = nzb, nzt+1
              hyp(k)    = surface_pressure * 100.0_wp * &
-                         ( (t_surface - g/cp * zu(k)) / t_surface )**(1.0_wp/0.286_wp)
+                         ( ( t_surface - g / cp * zu(k) ) /                    &
+                         t_surface )**(1.0_wp / 0.286_wp)
              pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
              t_d_pt(k) = 1.0_wp / pt_d_t(k)
@@ -384,42 +176,28 @@
        ENDIF
 
-
-       dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
-!
-!--    Adjust unrealistic values
-       IF ( precipitation )  CALL adjust_cloud( i,j ) 
-!
-!--    Use 1-d arrays
-       q_1d(:)  = q(:,j,i)
-       pt_1d(:) = pt(:,j,i)
-       qc_1d(:) = qc(:,j,i)
-       nc_1d(:) = nc_const
-       IF ( precipitation )  THEN
-          qr_1d(:) = qr(:,j,i)
-          nr_1d(:) = nr(:,j,i)
+!
+!--    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
+
 !
 !--    Compute cloud physics
        IF ( precipitation )  THEN
-          CALL autoconversion( i,j )
-          CALL accretion( i,j )
-          CALL selfcollection_breakup( i,j )
-          CALL evaporation_rain( i,j )
-          CALL sedimentation_rain( i,j )
+          CALL adjust_cloud
+          CALL autoconversion
+          CALL accretion
+          CALL selfcollection_breakup
+          CALL evaporation_rain
+          CALL sedimentation_rain
        ENDIF
 
-       IF ( drizzle )  CALL sedimentation_cloud( i,j )
-!
-!--    Derive tendencies
-       tend_q(:,j,i)  = ( q_1d(:) - q(:,j,i) ) / dt_micro
-       tend_pt(:,j,i) = ( pt_1d(:) - pt(:,j,i) ) / dt_micro
-       IF ( precipitation )  THEN
-          tend_qr(:,j,i) = ( qr_1d(:) - qr(:,j,i) ) / dt_micro
-          tend_nr(:,j,i) = ( nr_1d(:) - nr(:,j,i) ) / dt_micro
-       ENDIF
-
-    END SUBROUTINE microphysics_control_ij
-
-    SUBROUTINE adjust_cloud_ij( i, j )
+       IF ( drizzle )  CALL sedimentation_cloud
+
+    END SUBROUTINE microphysics_control
+
+    SUBROUTINE adjust_cloud
 
        USE arrays_3d,                                                          &
@@ -429,6 +207,9 @@
            ONLY:  eps_sb, xrmin, xrmax, hyrho, k_cc, x0
 
+       USE cpulog,                                                             &
+           ONLY:  cpu_log, log_point_s
+
        USE indices,                                                            &
-           ONLY:  nzb, nzb_s_inner, nzt
+           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
 
        USE kinds
@@ -439,50 +220,55 @@
        INTEGER(iwp) ::  j                 !:
        INTEGER(iwp) ::  k                 !:
-!
-!--    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.  
-       DO  k = nzb_s_inner(j,i)+1, nzt
-
-          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
-             ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
-                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
-             ENDIF
-
-          ENDIF
-
+
+       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                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
+                   ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
+                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
+                   ENDIF
+
+                ENDIF
+             ENDDO
+          ENDDO
        ENDDO
 
-    END SUBROUTINE adjust_cloud_ij
-
-
-    SUBROUTINE autoconversion_ij( i, j )
+       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
+
+    END SUBROUTINE adjust_cloud
+
+
+    SUBROUTINE autoconversion
 
        USE arrays_3d,                                                          &
-           ONLY:  diss, dzu, nc_1d, nr_1d, qc_1d, qr_1d
+           ONLY:  diss, dzu, nr, qc, qr
 
        USE cloud_parameters,                                                   &
            ONLY:  a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3,        &
-                  c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air, x0
+                  c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air,         &
+                  nc_const, x0
 
        USE control_parameters,                                                 &
            ONLY:  dt_micro, rho_surface, turbulence
 
+       USE cpulog,                                                             &
+           ONLY:  cpu_log, log_point_s
+
        USE grid_variables,                                                     &
            ONLY:  dx, dy
 
        USE indices,                                                            &
-           ONLY:  nzb, nzb_s_inner, nzt
+           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
 
        USE kinds
@@ -496,5 +282,5 @@
        REAL(wp)     ::  alpha_cc          !:                   
        REAL(wp)     ::  autocon           !:
-       REAL(wp)     ::  epsilon           !:
+       REAL(wp)     ::  dissipation       !:
        REAL(wp)     ::  k_au              !:
        REAL(wp)     ::  l_mix             !:
@@ -509,77 +295,94 @@
        REAL(wp)     ::  xc                !:
 
-       k_au = k_cc / ( 20.0_wp * x0 )
-
-       DO  k = nzb_s_inner(j,i)+1, nzt
-
-          IF ( qc_1d(k) > eps_sb )  THEN
-!
-!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
-!--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
-             tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
-!
-!--          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_1d(k) / nc_1d(k)
-!
-!--          Parameterized turbulence effects on autoconversion (Seifert, 
-!--          Nuijens and Stevens, 2010)
-             IF ( 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)
-                epsilon = 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 ) * epsilon**( 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 +                                             &
-                       epsilon * 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.0_wp * qc_1d(k)**2.0_wp * xc**2.0_wp *   &
-                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2.0_wp ) *        &
-                       rho_surface
-             autocon = MIN( autocon, qc_1d(k) / dt_micro )
-
-             qr_1d(k) = qr_1d(k) + autocon * dt_micro
-             qc_1d(k) = qc_1d(k) - autocon * dt_micro 
-             nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
-
-          ENDIF
-
+       CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_s_inner(j,i)+1, nzt
+
+                IF ( qc(k,j,i) > eps_sb )  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 = 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) )
+!
+!--                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_const
+!
+!--                Parameterized turbulence effects on autoconversion (Seifert, 
+!--                Nuijens and Stevens, 2010)
+                   IF ( 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
+                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro 
+                   nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro
+
+                ENDIF
+
+             ENDDO
+          ENDDO
        ENDDO
 
-    END SUBROUTINE autoconversion_ij
-
-
-    SUBROUTINE accretion_ij( i, j )
+       CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
+
+    END SUBROUTINE autoconversion
+
+
+    SUBROUTINE accretion
 
        USE arrays_3d,                                                          &
-           ONLY:  diss, qc_1d, qr_1d
+           ONLY:  diss, qc, qr
 
        USE cloud_parameters,                                                   &
@@ -589,6 +392,9 @@
            ONLY:  dt_micro, rho_surface, turbulence
 
+       USE cpulog,                                                             &
+           ONLY:  cpu_log, log_point_s
+
        USE indices,                                                            &
-           ONLY:  nzb, nzb_s_inner, nzt
+           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
 
        USE kinds
@@ -606,44 +412,54 @@
        REAL(wp)     ::  xc                !:
 
-       DO  k = nzb_s_inner(j,i)+1, nzt
-          IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb ) )  THEN
-!
-!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
-             tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 
-!
-!--          Universal function for accretion process 
-!--          (Seifert and Beheng, 2001):
-             phi_ac = tau_cloud / ( tau_cloud + 5.0E-5_wp ) 
-             phi_ac = ( phi_ac**2.0_wp )**2.0_wp
-!
-!--          Parameterized turbulence effects on autoconversion (Seifert, 
-!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to 
-!--          convert the dissipation (diss) from m2 s-3 to cm2 s-3.
-             IF ( 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_1d(k) * qr_1d(k) * phi_ac *                 &
-                    SQRT( rho_surface * hyrho(k) )
-             accr = MIN( accr, qc_1d(k) / dt_micro )
-
-             qr_1d(k) = qr_1d(k) + accr * dt_micro 
-             qc_1d(k) = qc_1d(k) - accr * dt_micro
-
-          ENDIF
-
+       CALL cpu_log( log_point_s(56), 'accretion', 'start' )
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_s_inner(j,i)+1, nzt
+
+                IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb ) )  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
+!
+!--                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 ( 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 
+                   qc(k,j,i) = qc(k,j,i) - accr * dt_micro
+
+                ENDIF
+
+             ENDDO
+          ENDDO
        ENDDO
 
-    END SUBROUTINE accretion_ij
-
-
-    SUBROUTINE selfcollection_breakup_ij( i, j )
+       CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
+
+    END SUBROUTINE accretion
+
+
+    SUBROUTINE selfcollection_breakup
 
        USE arrays_3d,                                                          &
-           ONLY:  nr_1d, qr_1d
+           ONLY:  nr, qr
 
        USE cloud_parameters,                                                   &
@@ -653,6 +469,9 @@
            ONLY:  dt_micro, rho_surface
 
+       USE cpulog,                                                             &
+           ONLY:  cpu_log, log_point_s
+
        USE indices,                                                            &
-           ONLY:  nzb, nzb_s_inner, nzt
+           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
 
        USE kinds
@@ -669,32 +488,42 @@
        REAL(wp)     ::  selfcoll          !:
 
-       DO  k = nzb_s_inner(j,i)+1, nzt
-          IF ( qr_1d(k) > eps_sb )  THEN
-!
-!--          Selfcollection rate (Seifert and Beheng, 2001):
-             selfcoll = k_rr * nr_1d(k) * qr_1d(k) *         &
-                        SQRT( hyrho(k) * rho_surface )
-!
-!--          Weight averaged diameter of rain drops:
-             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * 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_1d(k) / dt_micro )
-             nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
-
-          ENDIF          
+       CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                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
+
+                ENDIF          
+             ENDDO
+          ENDDO
        ENDDO
 
-    END SUBROUTINE selfcollection_breakup_ij
-
-
-    SUBROUTINE evaporation_rain_ij( i, j )
+       CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
+
+    END SUBROUTINE selfcollection_breakup
+
+
+    SUBROUTINE evaporation_rain
+
 !
 !--    Evaporation of precipitable water. Condensation is neglected for 
@@ -702,5 +531,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  hyp, nr_1d, pt_1d, q_1d,  qc_1d, qr_1d
+           ONLY:  hyp, nr, pt, q,  qc, qr
 
        USE cloud_parameters,                                                   &
@@ -716,6 +545,9 @@
            ONLY:  dt_micro
 
+       USE cpulog,                                                             &
+           ONLY:  cpu_log, log_point_s
+
        USE indices,                                                            &
-           ONLY:  nzb, nzb_s_inner, nzt
+           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
 
        USE kinds
@@ -745,96 +577,118 @@
        REAL(wp)     ::  xr                !:
 
-       DO  k = nzb_s_inner(j,i)+1, nzt
-          IF ( qr_1d(k) > eps_sb )  THEN
-!
-!--          Actual liquid water temperature:
-             t_l = t_d_pt(k) * pt_1d(k)
-!
-!--          Saturation vapor pressure at t_l:
-             e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / ( t_l - 35.86_wp ) )
-!
-!--          Computation of saturation humidity:
-             q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
-             alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
-             q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
-!
-!--          Supersaturation:
-             sat = MIN( 0.0_wp, ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp )
-!
-!--          Actual temperature:
-             temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
-    
-             g_evap = 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 rain drops
-             xr = hyrho(k) * qr_1d(k) / nr_1d(k)
-!
-!--          Weight averaged diameter of rain drops:
-             dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
-!
-!--          Compute ventilation factor and intercept parameter 
-!--          (Seifert and Beheng, 2006; Seifert, 2008):
-             IF ( ventilation_effect )  THEN
-!
-!--             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
-
-                mu_r_2   = mu_r + 2.0_wp
-                mu_r_5d2 = mu_r + 2.5_wp 
-                f_vent = a_vent * gamm( mu_r_2 ) *                            &
-                         lambda_r**( -mu_r_2 ) +                              &
-                         b_vent * schmidt_p_1d3 *                             &
-                         SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *    &
-                         lambda_r**( -mu_r_5d2 ) *                            &
-                         ( 1.0_wp - 0.5_wp * ( b_term / a_term ) *            &
-                         ( lambda_r /                                         &
-                         (       c_term + lambda_r ) )**mu_r_5d2 -            &
-                                 0.125_wp * ( b_term / a_term )**2.0_wp *     &
-                         ( lambda_r /                                         &
-                         ( 2.0_wp * c_term + lambda_r ) )**mu_r_5d2 -         &
-                                 0.0625_wp * ( b_term / a_term )**3.0_wp *    &
-                         ( lambda_r /                                         &
-                         ( 3.0_wp * c_term + lambda_r ) )**mu_r_5d2 -         &
-                                 0.0390625_wp * ( b_term / a_term )**4.0_wp * &
-                         ( lambda_r /                                         &
-                         ( 4.0_wp * c_term + lambda_r ) )**mu_r_5d2 )
-                nr_0   = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) /             &
-                         gamm( mu_r + 1.0_wp ) 
-             ELSE
-                f_vent = 1.0_wp
-                nr_0   = nr_1d(k) * dr
-             ENDIF
-!
-!--          Evaporation rate of rain water content (Seifert and Beheng, 2006):
-             evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat /    &
-                    hyrho(k)
-
-             evap    = MAX( evap, -qr_1d(k) / dt_micro )
-             evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
-                            -nr_1d(k) / dt_micro )
-
-             qr_1d(k) = qr_1d(k) + evap * dt_micro
-             nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
-          ENDIF          
-
+       CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                IF ( qr(k,j,i) > eps_sb )  THEN
+!
+!--                Actual liquid water temperature:
+                   t_l = t_d_pt(k) * pt(k,j,i)
+!
+!--                Saturation vapor pressure at t_l:
+                   e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /    &
+                                          ( t_l - 35.86_wp )                   &
+                                        )
+!
+!--                Computation of saturation humidity:
+                   q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
+                   alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
+                   q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) /               &
+                           ( 1.0_wp + alpha * q_s )
+!
+!--                Supersaturation:
+                   sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
+!
+!--                Evaporation needs only to be calculated in subsaturated regions
+                   IF ( sat < 0.0_wp )  THEN
+!
+!--                   Actual temperature:
+                      temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
+             
+                      g_evap = 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 rain drops
+                      xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
+!
+!--                   Weight averaged diameter of rain drops:
+                      dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
+!
+!--                   Compute ventilation factor and intercept parameter 
+!--                   (Seifert and Beheng, 2006; Seifert, 2008):
+                      IF ( ventilation_effect )  THEN
+!
+!--                      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
+
+                         mu_r_2   = mu_r + 2.0_wp
+                         mu_r_5d2 = mu_r + 2.5_wp 
+
+                         f_vent = a_vent * gamm( mu_r_2 ) *                     &
+                                  lambda_r**( -mu_r_2 ) + b_vent *              &
+                                  schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
+                                  gamm( mu_r_5d2 ) * lambda_r**( -mu_r_5d2 ) *  &
+                                  ( 1.0_wp -                                    &
+                                    0.5_wp * ( b_term / a_term ) *              &
+                                    ( lambda_r / ( c_term + lambda_r )          &
+                                    )**mu_r_5d2 -                               &
+                                    0.125_wp * ( b_term / a_term )**2 *         &
+                                    ( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
+                                    )**mu_r_5d2 -                               &
+                                    0.0625_wp * ( b_term / a_term )**3 *        &
+                                    ( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
+                                    )**mu_r_5d2 -                               &
+                                    0.0390625_wp * ( b_term / a_term )**4 *     & 
+                                    ( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
+                                    )**mu_r_5d2                                 &
+                                  )
+
+                         nr_0   = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) /    &
+                                  gamm( mu_r + 1.0_wp ) 
+                      ELSE
+                         f_vent = 1.0_wp
+                         nr_0   = nr(k,j,i) * dr
+                      ENDIF
+   !
+   !--                Evaporation rate of rain water content (Seifert and 
+   !--                Beheng, 2006):
+                      evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat /   &
+                                hyrho(k)
+                      evap    = MAX( evap, -qr(k,j,i) / dt_micro )
+                      evap_nr = MAX( c_evap * evap / xr * hyrho(k),            &
+                                     -nr(k,j,i) / dt_micro )
+
+                      qr(k,j,i) = qr(k,j,i) + evap    * dt_micro
+                      nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro
+
+                   ENDIF
+                ENDIF          
+
+             ENDDO
+          ENDDO
        ENDDO
 
-    END SUBROUTINE evaporation_rain_ij
-
-
-    SUBROUTINE sedimentation_cloud_ij( i, j )
+       CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
+
+    END SUBROUTINE evaporation_rain
+
+
+    SUBROUTINE sedimentation_cloud
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, dzu, nc_1d, pt_1d, q_1d, qc_1d
+           ONLY:  ddzu, dzu, pt, q, qc
 
        USE cloud_parameters,                                                   &
-           ONLY:  eps_sb, hyrho, k_st, l_d_cp, prr, pt_d_t, rho_l, sigma_gc
+           ONLY:  eps_sb, hyrho, l_d_cp, nc_const, pt_d_t, sed_qc_const
 
        USE constants,                                                          &
@@ -844,6 +698,9 @@
            ONLY:  dt_do2d_xy, dt_micro, intermediate_timestep_count
 
+       USE cpulog,                                                             &
+           ONLY:  cpu_log, log_point_s
+
        USE indices,                                                            &
-           ONLY:  nzb, nzb_s_inner, nzt
+           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
 
        USE kinds
@@ -855,43 +712,48 @@
        INTEGER(iwp) ::  k                 !:
 
-       REAL(wp)     ::  sed_qc_const      !:
-
-
-       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc
-
-!
-!--    Sedimentation of cloud droplets (Heus et al., 2010):
-       sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l ))**( 2.0_wp / 3.0_wp ) *   &
-                     EXP( 5.0_wp * LOG( sigma_gc )**2 )
-
+       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !:
+
+       CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
+
+!
+!--    Sedimentation of cloud droplets (Ackermann et al., 2009, MWR):
        sed_qc(nzt+1) = 0.0_wp
 
-       DO  k = nzt, nzb_s_inner(j,i)+1, -1
-          IF ( qc_1d(k) > eps_sb )  THEN
-             sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0_wp / 3.0_wp ) * &
-                        ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )
-          ELSE
-             sed_qc(k) = 0.0_wp
-          ENDIF
-
-          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) /     &
-                                      dt_micro + sed_qc(k+1) )
-
-          q_1d(k)  = q_1d(k)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /  &
-                                hyrho(k) * dt_micro
-          qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 
-                                hyrho(k) * dt_micro
-          pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
-                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
-
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzt, nzb_s_inner(j,i)+1, -1
+
+                IF ( qc(k,j,i) > eps_sb )  THEN
+                   sed_qc(k) = sed_qc_const * nc_const**( -2.0_wp / 3.0_wp ) * &
+                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp )
+                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)             &
+                               )
+
+                q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) *          &
+                                        ddzu(k+1) / hyrho(k) * dt_micro
+                qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) *          &
+                                        ddzu(k+1) / hyrho(k) * dt_micro
+                pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) *          &
+                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
+                                        pt_d_t(k) * dt_micro
+
+             ENDDO
+          ENDDO
        ENDDO
 
-    END SUBROUTINE sedimentation_cloud_ij
-
-
-    SUBROUTINE sedimentation_rain_ij( i, j )
+       CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
+
+    END SUBROUTINE sedimentation_cloud
+
+
+    SUBROUTINE sedimentation_rain
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, dzu, nr_1d, pt_1d, q_1d, qr_1d
+           ONLY:  ddzu, dzu, nr, pt, q, qr
 
        USE cloud_parameters,                                                   &
@@ -901,10 +763,14 @@
 
        USE control_parameters,                                                 &
-           ONLY:  dt_do2d_xy, dt_micro, dt_3d, intermediate_timestep_count,    &
+           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro,     &
+                  dt_3d, intermediate_timestep_count,                          &
                   intermediate_timestep_count_max,                             &
                   precipitation_amount_interval, time_do2d_xy
 
+       USE cpulog,                                                             &
+           ONLY:  cpu_log, log_point_s
+
        USE indices,                                                            &
-           ONLY:  nzb, nzb_s_inner, nzt
+           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
 
        USE kinds
@@ -931,19 +797,836 @@
        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) ::  d_nr      !:
-      REAL(wp), DIMENSION(nzb:nzt+1) ::  d_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      !:
-
-
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !:
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !:
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_nr     !:
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_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' )
 !
 !--    Computation of sedimentation flux. Implementation according to Stevens
-!--    and Seifert (2008).
+!--    and Seifert (2008). Code is based on UCLA-LES.
+       IF ( intermediate_timestep_count == 1 )  prr(:,:,:) = 0.0_wp
+!
+!--    Compute velocities 
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                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 ) )          &
+                                              )                                &
+                                )
+
+                   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 ) )          &
+                                             )                                 &
+                                )
+                ELSE
+                   w_nr(k) = 0.0_wp
+                   w_qr(k) = 0.0_wp
+                ENDIF
+             ENDDO
+!
+!--          Adjust boundary values
+             w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
+             w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
+             w_nr(nzt+1) = 0.0_wp
+             w_qr(nzt+1) = 0.0_wp
+!
+!--          Compute Courant number
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                c_nr(k) = 0.25_wp * ( w_nr(k-1) +                              &
+                                      2.0_wp * w_nr(k) + w_nr(k+1) ) *         &
+                          dt_micro * ddzu(k)
+                c_qr(k) = 0.25_wp * ( w_qr(k-1) +                              &
+                                      2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
+                          dt_micro * ddzu(k)
+             ENDDO     
+!
+!--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
+             IF ( limiter_sedimentation )  THEN
+
+                DO k = nzb_s_inner(j,i)+1, nzt
+                   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 ) )
+
+                   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_s_inner(j,i)+1, -1
+!
+!--             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)
+                   z_run = z_run + dzu(k_run)
+                   k_run = k_run + 1
+                   c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
+                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
+                nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *          &
+                                        ddzu(k+1) / hyrho(k) * dt_micro
+!
+!--             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)
+                   z_run = z_run + dzu(k_run)
+                   k_run = k_run + 1
+                   c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
+
+                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
+
+                qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *          &
+                                        ddzu(k+1) / hyrho(k) * dt_micro
+                q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
+                                        ddzu(k+1) / hyrho(k) * dt_micro 
+                pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
+                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
+                                        pt_d_t(k) * dt_micro
+!
+!--             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)
+                ELSE
+                   prr(k,j,i) = sed_qr(k) / hyrho(k)
+                ENDIF
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+!
+!--    Precipitation amount
+       IF ( intermediate_timestep_count == intermediate_timestep_count_max     &
+            .AND.  ( dt_do2d_xy - time_do2d_xy ) <                             &
+                   precipitation_amount_interval )  THEN
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                precipitation_amount(j,i) = precipitation_amount(j,i) +        &
+                                            prr(nzb_s_inner(j,i)+1,j,i) *      &
+                                            hyrho(nzb_s_inner(j,i)+1) * dt_3d
+             ENDDO
+          ENDDO
+       ENDIF
+
+       CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
+
+    END SUBROUTINE sedimentation_rain
+
+
+!------------------------------------------------------------------------------!
+! Call for grid point i,j
+!------------------------------------------------------------------------------!
+
+    SUBROUTINE microphysics_control_ij( i, j )
+
+       USE arrays_3d,                                                          &
+           ONLY:  hyp, nc_1d, nr, nr_1d, pt, pt_init, pt_1d, q, q_1d, qc,      &
+                  qc_1d, qr, qr_1d, zu
+
+       USE cloud_parameters,                                                   &
+           ONLY:  cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt
+
+       USE control_parameters,                                                 &
+           ONLY:  call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, &
+                  g, intermediate_timestep_count, large_scale_forcing,         &
+                  lsf_surf, precipitation, pt_surface,                         &
+                  rho_surface,surface_pressure
+
+       USE indices,                                                            &
+           ONLY:  nzb, nzt
+
+       USE kinds
+
+       USE statistics,                                                         &
+           ONLY:  weight_pres
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !:
+       INTEGER(iwp) ::  j                 !:
+       INTEGER(iwp) ::  k                 !:
+
+       REAL(wp)     ::  t_surface         !:
+
+       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
+!
+!--       Calculate:
+!--       pt / t : ratio of potential and actual temperature (pt_d_t)
+!--       t / pt : ratio of actual and potential temperature (t_d_pt)
+!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
+          t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
+          DO  k = nzb, nzt+1
+             hyp(k)    = surface_pressure * 100.0_wp * &
+                         ( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
+             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
+             t_d_pt(k) = 1.0_wp / pt_d_t(k)
+             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
+          ENDDO
+!
+!--       Compute reference density
+          rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
+       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
+
+!
+!--    Use 1d arrays
+       q_1d(:)  = q(:,j,i)
+       pt_1d(:) = pt(:,j,i)
+       qc_1d(:) = qc(:,j,i)
+       nc_1d(:) = nc_const
+       IF ( precipitation )  THEN
+          qr_1d(:) = qr(:,j,i)
+          nr_1d(:) = nr(:,j,i)
+       ENDIF
+
+!
+!--    Compute cloud physics
+       IF ( precipitation )  THEN
+          CALL adjust_cloud( 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 )
+       ENDIF
+
+       IF ( drizzle )  CALL sedimentation_cloud( i,j )
+
+!
+!--    Store results on the 3d arrays
+       q(:,j,i)  = q_1d(:)
+       pt(:,j,i) = pt_1d(:)
+       IF ( precipitation )  THEN
+          qr(:,j,i) = qr_1d(:)
+          nr(:,j,i) = nr_1d(:)
+       ENDIF
+
+    END SUBROUTINE microphysics_control_ij
+
+    SUBROUTINE adjust_cloud_ij( i, j )
+
+       USE arrays_3d,                                                          &
+           ONLY:  qr_1d, nr_1d
+
+       USE cloud_parameters,                                                   &
+           ONLY:  eps_sb, xrmin, xrmax, hyrho, k_cc, x0
+
+       USE indices,                                                            &
+           ONLY:  nzb, nzb_s_inner, nzt
+
+       USE kinds
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !:
+       INTEGER(iwp) ::  j                 !:
+       INTEGER(iwp) ::  k                 !:
+!
+!--    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.  
+       DO  k = nzb_s_inner(j,i)+1, nzt
+
+          IF ( qr_1d(k) <= eps_sb )  THEN
+             qr_1d(k) = 0.0_wp
+             nr_1d(k) = 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_1d(k) * xrmin > qr_1d(k) * hyrho(k) )  THEN
+                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin
+             ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) )  THEN
+                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax
+             ENDIF
+
+          ENDIF
+
+       ENDDO
+
+    END SUBROUTINE adjust_cloud_ij
+
+
+    SUBROUTINE autoconversion_ij( i, j )
+
+       USE arrays_3d,                                                          &
+           ONLY:  diss, dzu, nc_1d, nr_1d, qc_1d, qr_1d
+
+       USE cloud_parameters,                                                   &
+           ONLY:  a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3,        &
+                  c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air, x0
+
+       USE control_parameters,                                                 &
+           ONLY:  dt_micro, rho_surface, turbulence
+
+       USE grid_variables,                                                     &
+           ONLY:  dx, dy
+
+       USE indices,                                                            &
+           ONLY:  nzb, nzb_s_inner, nzt
+
+       USE kinds
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !:
+       INTEGER(iwp) ::  j                 !:
+       INTEGER(iwp) ::  k                 !:
+
+       REAL(wp)     ::  alpha_cc          !:                   
+       REAL(wp)     ::  autocon           !:
+       REAL(wp)     ::  dissipation       !:
+       REAL(wp)     ::  k_au              !:
+       REAL(wp)     ::  l_mix             !:
+       REAL(wp)     ::  nu_c              !:
+       REAL(wp)     ::  phi_au            !:
+       REAL(wp)     ::  r_cc              !:
+       REAL(wp)     ::  rc                !:
+       REAL(wp)     ::  re_lambda         !:
+       REAL(wp)     ::  selfcoll          !:
+       REAL(wp)     ::  sigma_cc          !:
+       REAL(wp)     ::  tau_cloud         !:
+       REAL(wp)     ::  xc                !:
+
+       DO  k = nzb_s_inner(j,i)+1, nzt
+
+          IF ( qc_1d(k) > eps_sb )  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_1d(k) ))
+             tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
+!
+!--          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_1d(k) - 0.28_wp )
+!
+!--          Mean weight of cloud droplets:
+             xc = hyrho(k) * qc_1d(k) / nc_1d(k)
+!
+!--          Parameterized turbulence effects on autoconversion (Seifert, 
+!--          Nuijens and Stevens, 2010)
+             IF ( 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_1d(k)**2 * xc**2 *            &
+                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
+                       rho_surface
+             autocon = MIN( autocon, qc_1d(k) / dt_micro )
+
+             qr_1d(k) = qr_1d(k) + autocon * dt_micro
+             qc_1d(k) = qc_1d(k) - autocon * dt_micro 
+             nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
+
+          ENDIF
+
+       ENDDO
+
+    END SUBROUTINE autoconversion_ij
+
+
+    SUBROUTINE accretion_ij( i, j )
+
+       USE arrays_3d,                                                          &
+           ONLY:  diss, qc_1d, qr_1d
+
+       USE cloud_parameters,                                                   &
+           ONLY:  eps_sb, hyrho, k_cr0
+
+       USE control_parameters,                                                 &
+           ONLY:  dt_micro, rho_surface, turbulence
+
+       USE indices,                                                            &
+           ONLY:  nzb, nzb_s_inner, nzt
+
+       USE kinds
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !:
+       INTEGER(iwp) ::  j                 !:
+       INTEGER(iwp) ::  k                 !:
+
+       REAL(wp)     ::  accr              !:
+       REAL(wp)     ::  k_cr              !:
+       REAL(wp)     ::  phi_ac            !:
+       REAL(wp)     ::  tau_cloud         !:
+       REAL(wp)     ::  xc                !:
+
+       DO  k = nzb_s_inner(j,i)+1, nzt
+          IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb ) )  THEN
+!
+!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
+             tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 
+!
+!--          Universal function for accretion process 
+!--          (Seifert and Beheng, 2001):
+             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
+!
+!--          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 ( 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_1d(k) * qr_1d(k) * phi_ac * SQRT( rho_surface * hyrho(k) )
+             accr = MIN( accr, qc_1d(k) / dt_micro )
+
+             qr_1d(k) = qr_1d(k) + accr * dt_micro 
+             qc_1d(k) = qc_1d(k) - accr * dt_micro
+
+          ENDIF
+
+       ENDDO
+
+    END SUBROUTINE accretion_ij
+
+
+    SUBROUTINE selfcollection_breakup_ij( i, j )
+
+       USE arrays_3d,                                                          &
+           ONLY:  nr_1d, qr_1d
+
+       USE cloud_parameters,                                                   &
+           ONLY:  dpirho_l, eps_sb, hyrho, k_br, k_rr
+
+       USE control_parameters,                                                 &
+           ONLY:  dt_micro, rho_surface
+
+       USE indices,                                                            &
+           ONLY:  nzb, nzb_s_inner, nzt
+
+       USE kinds
+   
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !:
+       INTEGER(iwp) ::  j                 !:
+       INTEGER(iwp) ::  k                 !:
+
+       REAL(wp)     ::  breakup           !:
+       REAL(wp)     ::  dr                !:
+       REAL(wp)     ::  phi_br            !:
+       REAL(wp)     ::  selfcoll          !:
+
+       DO  k = nzb_s_inner(j,i)+1, nzt
+          IF ( qr_1d(k) > eps_sb )  THEN
+!
+!--          Selfcollection rate (Seifert and Beheng, 2001):
+             selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
+!
+!--          Weight averaged diameter of rain drops:
+             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * 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_1d(k) / dt_micro )
+             nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
+
+          ENDIF          
+       ENDDO
+
+    END SUBROUTINE selfcollection_breakup_ij
+
+
+    SUBROUTINE evaporation_rain_ij( i, j )
+!
+!--    Evaporation of precipitable water. Condensation is neglected for 
+!--    precipitable water.
+
+       USE arrays_3d,                                                          &
+           ONLY:  hyp, nr_1d, pt_1d, q_1d,  qc_1d, qr_1d
+
+       USE cloud_parameters,                                                   &
+           ONLY:  a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,&
+                  dpirho_l, eps_sb, hyrho, kin_vis_air, k_st, l_d_cp, l_d_r,   &
+                  l_v, rho_l, r_v, schmidt_p_1d3, thermal_conductivity_l,      &
+                  t_d_pt, ventilation_effect
+
+       USE constants,                                                          &
+           ONLY:  pi
+
+       USE control_parameters,                                                 &
+           ONLY:  dt_micro
+
+       USE indices,                                                            &
+           ONLY:  nzb, nzb_s_inner, nzt
+
+       USE kinds
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !:
+       INTEGER(iwp) ::  j                 !:
+       INTEGER(iwp) ::  k                 !:
+
+       REAL(wp)     ::  alpha             !:
+       REAL(wp)     ::  dr                !:
+       REAL(wp)     ::  e_s               !:
+       REAL(wp)     ::  evap              !:
+       REAL(wp)     ::  evap_nr           !:
+       REAL(wp)     ::  f_vent            !:
+       REAL(wp)     ::  g_evap            !:
+       REAL(wp)     ::  lambda_r          !:
+       REAL(wp)     ::  mu_r              !:
+       REAL(wp)     ::  mu_r_2            !:
+       REAL(wp)     ::  mu_r_5d2          !:
+       REAL(wp)     ::  nr_0              !:
+       REAL(wp)     ::  q_s               !:
+       REAL(wp)     ::  sat               !:
+       REAL(wp)     ::  t_l               !:
+       REAL(wp)     ::  temp              !:
+       REAL(wp)     ::  xr                !:
+
+       DO  k = nzb_s_inner(j,i)+1, nzt
+          IF ( qr_1d(k) > eps_sb )  THEN
+!
+!--          Actual liquid water temperature:
+             t_l = t_d_pt(k) * pt_1d(k)
+!
+!--          Saturation vapor pressure at t_l:
+             e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /          &
+                                    ( t_l - 35.86_wp )                         &
+                                  )
+!
+!--          Computation of saturation humidity:
+             q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
+             alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
+             q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
+!
+!--          Supersaturation:
+             sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
+!
+!--          Evaporation needs only to be calculated in subsaturated regions
+             IF ( sat < 0.0_wp )  THEN
+!
+!--             Actual temperature:
+                temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
+       
+                g_evap = 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 rain drops
+                xr = hyrho(k) * qr_1d(k) / nr_1d(k)
+!
+!--             Weight averaged diameter of rain drops:
+                dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
+!
+!--             Compute ventilation factor and intercept parameter 
+!--             (Seifert and Beheng, 2006; Seifert, 2008):
+                IF ( ventilation_effect )  THEN
+!
+!--                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
+
+                   mu_r_2   = mu_r + 2.0_wp
+                   mu_r_5d2 = mu_r + 2.5_wp 
+
+                   f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) +  & 
+                            b_vent * schmidt_p_1d3 *                           &
+                            SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *  &
+                            lambda_r**( -mu_r_5d2 ) *                          &
+                            ( 1.0_wp -                                         &
+                              0.5_wp * ( b_term / a_term ) *                   &
+                              ( lambda_r / ( c_term + lambda_r )               &
+                              )**mu_r_5d2 -                                    &
+                              0.125_wp * ( b_term / a_term )**2 *              &
+                              ( lambda_r / ( 2.0_wp * c_term + lambda_r )      &
+                              )**mu_r_5d2 -                                    &
+                              0.0625_wp * ( b_term / a_term )**3 *             &
+                              ( lambda_r / ( 3.0_wp * c_term + lambda_r )      &
+                              )**mu_r_5d2 -                                    &
+                              0.0390625_wp * ( b_term / a_term )**4 *          & 
+                              ( lambda_r / ( 4.0_wp * c_term + lambda_r )      &
+                              )**mu_r_5d2                                      &
+                            )
+
+                   nr_0   = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) /           &
+                            gamm( mu_r + 1.0_wp ) 
+                ELSE
+                   f_vent = 1.0_wp
+                   nr_0   = nr_1d(k) * dr
+                ENDIF
+!
+!--             Evaporation rate of rain water content (Seifert and Beheng, 2006):
+                evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
+                evap    = MAX( evap, -qr_1d(k) / dt_micro )
+                evap_nr = MAX( c_evap * evap / xr * hyrho(k),                  &
+                               -nr_1d(k) / dt_micro )
+
+                qr_1d(k) = qr_1d(k) + evap    * dt_micro
+                nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
+
+             ENDIF
+          ENDIF          
+
+       ENDDO
+
+    END SUBROUTINE evaporation_rain_ij
+
+
+    SUBROUTINE sedimentation_cloud_ij( i, j )
+
+       USE arrays_3d,                                                          &
+           ONLY:  ddzu, dzu, nc_1d, pt_1d, q_1d, qc_1d
+
+       USE cloud_parameters,                                                   &
+           ONLY:  eps_sb, hyrho, l_d_cp, pt_d_t, sed_qc_const
+
+       USE constants,                                                          &
+           ONLY:  pi
+
+       USE control_parameters,                                                 &
+           ONLY:  dt_do2d_xy, dt_micro, intermediate_timestep_count
+
+       USE indices,                                                            &
+           ONLY:  nzb, nzb_s_inner, nzt
+
+       USE kinds
+       
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !:
+       INTEGER(iwp) ::  j                 !:
+       INTEGER(iwp) ::  k                 !:
+
+       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc  !:
+
+!
+!--    Sedimentation of cloud droplets (Ackermann et al., 2009, MWR):
+       sed_qc(nzt+1) = 0.0_wp
+
+       DO  k = nzt, nzb_s_inner(j,i)+1, -1
+          IF ( qc_1d(k) > eps_sb )  THEN
+             sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0_wp / 3.0_wp ) *       &
+                         ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )
+          ELSE
+             sed_qc(k) = 0.0_wp
+          ENDIF
+
+          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) /          &
+                                      dt_micro + sed_qc(k+1)                   &
+                         )
+
+          q_1d(k)  = q_1d(k)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
+                                hyrho(k) * dt_micro
+          qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      & 
+                                hyrho(k) * dt_micro
+          pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
+                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
+
+       ENDDO
+
+    END SUBROUTINE sedimentation_cloud_ij
+
+
+    SUBROUTINE sedimentation_rain_ij( i, j )
+
+       USE arrays_3d,                                                          &
+           ONLY:  ddzu, dzu, nr_1d, pt_1d, q_1d, qr_1d
+
+       USE cloud_parameters,                                                   &
+           ONLY:  a_term, b_term, c_term, cof, dpirho_l, eps_sb, hyrho,        &
+                  limiter_sedimentation, l_d_cp, precipitation_amount, prr,    &
+                  pt_d_t, stp
+
+       USE control_parameters,                                                 &
+           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro,     &
+                  dt_3d, intermediate_timestep_count,                          &
+                  intermediate_timestep_count_max,                             &
+                  precipitation_amount_interval, time_do2d_xy
+
+       USE indices,                                                            &
+           ONLY:  nzb, nzb_s_inner, nzt
+
+       USE kinds
+
+       USE statistics,                                                         &
+           ONLY:  weight_substep
+       
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                          !:
+       INTEGER(iwp) ::  j                          !:
+       INTEGER(iwp) ::  k                          !:
+       INTEGER(iwp) ::  k_run                      !:
+
+       REAL(wp)     ::  c_run                      !:
+       REAL(wp)     ::  d_max                      !:
+       REAL(wp)     ::  d_mean                     !:
+       REAL(wp)     ::  d_min                      !:
+       REAL(wp)     ::  dr                         !:
+       REAL(wp)     ::  dt_sedi                    !:
+       REAL(wp)     ::  flux                       !:
+       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) ::  d_nr     !:
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  d_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     !:
+
+
+!
+!--    Computation of sedimentation flux. Implementation according to Stevens
+!--    and Seifert (2008). Code is based on UCLA-LES.
        IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
 !
@@ -960,11 +1643,19 @@
 !
 !--          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 ) ) ) )
-             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 ) ) ) )
+             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 ) )                &
+                                        )                                      &
+                          )
+             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 ) )                &
+                                       )                                       &
+                          )
           ELSE
              w_nr(k) = 0.0_wp
@@ -981,7 +1672,7 @@
 !--    Compute Courant number
        DO  k = nzb_s_inner(j,i)+1, nzt
-          c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) * &
+          c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) *   &
                     dt_micro * ddzu(k)
-          c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * &
+          c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) *   &
                     dt_micro * ddzu(k)
        ENDDO     
@@ -995,6 +1686,7 @@
              d_max  = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
 
-             qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, 2.0_wp * d_max, &
-                                                     ABS( d_mean ) )
+             qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
+                                                        2.0_wp * d_max,        &
+                                                        ABS( d_mean ) )
 
              d_mean = 0.5_wp * ( nr_1d(k+1) + nr_1d(k-1) )
@@ -1002,6 +1694,7 @@
              d_max  = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
 
-             nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, 2.0_wp * d_max, &
-                                                     ABS( d_mean ) )
+             nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
+                                                        2.0_wp * d_max,        &
+                                                        ABS( d_mean ) )
           ENDDO
 
@@ -1026,6 +1719,6 @@
           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_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) *     &
+             flux  = flux + hyrho(k_run) *                                     &
+                     ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) *   &
                      0.5_wp ) * c_run * dzu(k_run)
              z_run = z_run + dzu(k_run)
@@ -1036,10 +1729,10 @@
 !--       It is not allowed to sediment more rain drop number density than 
 !--       available
-          flux = MIN( flux,                                                   &
+          flux = MIN( flux,                                                    &
                       hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
 
           sed_nr(k) = flux / dt_micro
-          nr_1d(k)  = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /    &
-                                 hyrho(k) * dt_micro
+          nr_1d(k)  = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /     &
+                                    hyrho(k) * dt_micro
 !
 !--       Sum up all rain water content which contributes to the flux 
@@ -1050,8 +1743,8 @@
           c_run = MIN( 1.0_wp, c_qr(k) )
 
-          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt-1 )
-
-             flux  = flux + hyrho(k_run) *                                    &
-                     ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) *    &
+          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
+
+             flux  = flux + hyrho(k_run) *                                     &
+                     ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) *   &
                      0.5_wp ) * c_run * dzu(k_run)
              z_run = z_run + dzu(k_run)
@@ -1062,29 +1755,34 @@
 !
 !--       It is not allowed to sediment more rain water content than available
-          flux = MIN( flux,                                                   &
+          flux = MIN( flux,                                                    &
                       hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
 
           sed_qr(k) = flux / dt_micro
 
-          qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
+          qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
                                 hyrho(k) * dt_micro
-          q_1d(k)  = q_1d(k)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
+          q_1d(k)  = q_1d(k)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
                                 hyrho(k) * dt_micro 
-          pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
+          pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
                                 hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
 !
 !--       Compute the rain rate
-          prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) *                    &
-                       weight_substep(intermediate_timestep_count)
+          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)
+          ELSE
+             prr(k,j,i) = sed_qr(k) / hyrho(k)
+          ENDIF
+
        ENDDO
 
 !
 !--    Precipitation amount
-       IF ( intermediate_timestep_count == intermediate_timestep_count_max    &
-            .AND.  ( dt_do2d_xy - time_do2d_xy ) <                            &
-            precipitation_amount_interval )  THEN
-
-          precipitation_amount(j,i) = precipitation_amount(j,i) +   &
-                                      prr(nzb_s_inner(j,i)+1,j,i) *      &
+       IF ( intermediate_timestep_count == intermediate_timestep_count_max     &
+            .AND.  ( dt_do2d_xy - time_do2d_xy ) <                             &
+                   precipitation_amount_interval )  THEN
+
+          precipitation_amount(j,i) = precipitation_amount(j,i) +              &
+                                      prr(nzb_s_inner(j,i)+1,j,i) *            &
                                       hyrho(nzb_s_inner(j,i)+1) * dt_3d
        ENDIF
@@ -1092,5 +1790,7 @@
     END SUBROUTINE sedimentation_rain_ij
 
-
+!------------------------------------------------------------------------------!
+! Call for all optimizations
+!------------------------------------------------------------------------------!
 !
 !-- This function computes the gamma function (Press et al., 1992). 
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 1360)
+++ /palm/trunk/SOURCE/modules.f90	(revision 1361)
@@ -20,5 +20,7 @@
 ! Current revisions:
 ! ------------------
-! 
+! tend_* removed
+! call_microphysics_at_all_substeps added
+! default of drizzle set to true
 ! 
 ! Former revisions:
@@ -287,7 +289,7 @@
           diss_l_v, diss_l_w, flux_l_e, flux_l_nr, flux_l_pt, flux_l_q,        &
           flux_l_qr, flux_l_sa, flux_l_u, flux_l_v, flux_l_w, kh, km, lad_s,   &
-          lad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls, tend, tend_pt,    &
-          tend_nr, tend_q, tend_qr, tric, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l,   &
-          v_m_n, v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s
+          lad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls, tend, tric,       &
+          u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n, v_m_r, v_m_s, w_m_l,       &
+          w_m_n, w_m_r, w_m_s
 
 #if defined( __nopointer )
@@ -368,5 +370,5 @@
                  b_vent = 0.308_wp,    & !: coef. for ventilation effect
                  beta_cc = 3.09E-4_wp, & !: coef. in turb. parametrization (cm-2 s3)
-                 bfactor,           &
+                 bfactor,              & !: solution effect on diffusional growth
                  c_1 = 4.82E-6_wp,     & !: coef. in turb. parametrization (m)
                  c_2 = 4.8E-6_wp,      & !: coef. in turb. parametrization (m)
@@ -384,5 +386,5 @@
                 cp = 1005.0_wp,       & !: heat capacity of dry air (J kg-1 K-1) 
                 diff_coeff_l = 0.23E-4_wp, & !: diffusivity of water vapor (m2 s-1)
-                effective_coll_efficiency, & !:
+                effective_coll_efficiency, & !: effective collision efficiency
                 eps_ros = 1.0E-4_wp,  & !: accuracy of Rosenbrock method
                 eps_sb = 1.0E-20_wp,  & !: threshold in two-moments scheme
@@ -408,4 +410,5 @@
                 schmidt = 0.71_wp,    & !: Schmidt number
                 schmidt_p_1d3,        & !: schmidt**( 1.0 / 3.0 )
+                sed_qc_const,         & !: const. for sedimentation of cloud water
                 sigma_gc = 1.3_wp,    & !: log-normal geometric standard deviation
                 stp = 2.5066282746310005_wp, & !: parameter in gamma function
@@ -592,4 +595,5 @@
                 bc_lr_raddir = .FALSE., bc_ns_cyc = .TRUE., &
                 bc_ns_dirrad = .FALSE., bc_ns_raddir = .FALSE.,&
+                call_microphysics_at_all_substeps = .FALSE., &
                 call_psolver_at_all_substeps = .TRUE., &
                 cloud_droplets = .FALSE., cloud_physics = .FALSE., &
@@ -604,5 +608,5 @@
                 do_sum = .FALSE., &
                 dp_external = .FALSE., dp_smooth = .FALSE., &
-                drizzle = .FALSE., dt_fixed = .FALSE., &
+                drizzle = .TRUE., dt_fixed = .FALSE., &
                 dt_3d_reached, dt_3d_reached_l, exchange_mg = .FALSE., &
                 first_call_lpm = .TRUE., &
Index: /palm/trunk/SOURCE/parin.f90
===================================================================
--- /palm/trunk/SOURCE/parin.f90	(revision 1360)
+++ /palm/trunk/SOURCE/parin.f90	(revision 1361)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! +call_microphysics_at_all_substeps
 ! 
 ! Former revisions:
@@ -148,5 +148,6 @@
                bottom_salinityflux, building_height, building_length_x,        &
                building_length_y, building_wall_left, building_wall_south,     &
-               call_psolver_at_all_substeps, canopy_mode, canyon_height,       &
+               call_microphysics_at_all_substeps, call_psolver_at_all_substeps,&
+               canopy_mode, canyon_height,                                     &
                canyon_width_x, canyon_width_y, canyon_wall_left,               &
                canyon_wall_south, cfl_factor,                                  &
@@ -248,5 +249,6 @@
              bottom_salinityflux, building_height, building_length_x, &
              building_length_y, building_wall_left, building_wall_south, &
-             call_psolver_at_all_substeps, canopy_mode, canyon_height, &
+             call_psolver_at_all_substeps, call_microphysics_at_all_substeps, &
+             canopy_mode, canyon_height, &
              canyon_width_x, canyon_width_y, canyon_wall_left, &
              canyon_wall_south, c_sedimentation, cfl_factor, cloud_droplets, &
Index: /palm/trunk/SOURCE/prandtl_fluxes.f90
===================================================================
--- /palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 1360)
+++ /palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 1361)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix: calculation of turbulent fluxes of rain water content (qrsws) and rain 
+! drop concentration (nrsws) added
 ! 
 ! Former revisions:
@@ -64,12 +65,13 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  e, pt, q, qs, qsws, rif, shf, ts, u, us, usws, v, vpt, vsws,    &
-               zu, zw, z0, z0h
+        ONLY:  e, nr, nrs, nrsws, pt, q, qr, qrs, qrsws, qs, qsws, rif, shf,   &
+               ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h
 
     USE control_parameters,                                                    &
-        ONLY:  constant_heatflux, constant_waterflux, coupling_mode, g,        &
-               humidity, ibc_e_b, kappa, large_scale_forcing, lsf_surf,        &
-               passive_scalar, pt_surface, q_surface, rif_max, rif_min,        &
-               run_coupled, surface_pressure
+        ONLY:  cloud_physics, constant_heatflux, constant_waterflux,           &
+               coupling_mode, g, humidity, ibc_e_b, icloud_scheme, kappa,      &
+               large_scale_forcing, lsf_surf, passive_scalar, precipitation,   &
+               pt_surface, q_surface, rif_max, rif_min, run_coupled,           &
+               surface_pressure
 
     USE indices,                                                               &
@@ -96,6 +98,7 @@
 !
 !-- Data information for accelerators
-    !$acc data present( e, nzb_u_inner, nzb_v_inner, nzb_s_inner, pt, q, qs ) &
-    !$acc      present( qsws, rif, shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h )
+    !$acc data present( e, nrsws, nzb_u_inner, nzb_v_inner, nzb_s_inner, pt )  &
+    !$acc      present( q, qs, qsws, qrsws, rif, shf, ts, u, us, usws, v )     &
+    !$acc      present( vpt, vsws, zu, zw, z0, z0h )
 ! 
 !-- Compute theta*
@@ -383,4 +386,47 @@
           ENDDO
        ENDIF
+
+       IF ( cloud_physics  .AND.  icloud_scheme == 0  &
+            .AND.  precipitation )  THEN
+
+          !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
+          !$acc kernels loop independent
+          DO  i = nxlg, nxrg
+             !$acc loop independent
+             DO  j = nysg, nyng
+
+                k   = nzb_s_inner(j,i)
+                z_p = zu(k+1) - zw(k)
+
+                IF ( rif(j,i) >= 0.0 )  THEN
+!
+!--                Stable stratification
+                   qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) / (          &
+                                  LOG( z_p / z0h(j,i) ) +                      &
+                                  5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p )
+                   nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) / (          &
+                                  LOG( z_p / z0h(j,i) ) +                      &
+                                  5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p )
+
+                ELSE
+!
+!--                Unstable stratification
+                   a = SQRT( 1.0 - 16.0 * rif(j,i) ) 
+                   b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p ) 
+ 
+                   qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) /   (        &
+                             LOG( z_p / z0h(j,i) ) -                           &
+                             2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
+                   nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) /   (        &
+                             LOG( z_p / z0h(j,i) ) -                           &
+                             2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+       ENDIF
+
     ENDIF
 
@@ -396,4 +442,11 @@
        CALL exchange_horiz_2d( qsws )
        !$acc update device( qsws )
+       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
+            precipitation )  THEN
+          !$acc update host( qrsws, nrsws )
+          CALL exchange_horiz_2d( qrsws )
+          CALL exchange_horiz_2d( nrsws )
+          !$acc update device( qrsws, nrsws )
+       ENDIF
     ENDIF
 
@@ -425,4 +478,18 @@
 
 !
+!-- Compute (turbulent) fluxes of rain water content and rain drop concentartion
+    IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
+       !$OMP PARALLEL DO
+       !$acc kernels loop independent
+       DO  i = nxlg, nxrg
+          !$acc loop independent
+          DO  j = nysg, nyng
+             qrsws(j,i) = -qrs(j,i) * us(j,i)
+             nrsws(j,i) = -nrs(j,i) * us(j,i)
+          ENDDO
+       ENDDO
+    ENDIF
+
+!
 !-- Bottom boundary condition for the TKE
     IF ( ibc_e_b == 2 )  THEN
Index: /palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- /palm/trunk/SOURCE/prognostic_equations.f90	(revision 1360)
+++ /palm/trunk/SOURCE/prognostic_equations.f90	(revision 1361)
@@ -20,5 +20,10 @@
 ! Current revisions:
 ! ------------------
-! 
+! Two-moment microphysics moved to the start of prognostic equations. This makes
+! the 3d arrays for tend_q, tend_qr, tend_pt and tend_pt redundant.
+! Additionally, it is allowed to call the microphysics just once during the time
+! step (not at each sub-time step).
+!
+! Two-moment cloud physics added for vector and accelerator optimization.
 ! 
 ! Former revisions:
@@ -126,11 +131,11 @@
                nrswst, pt, ptdf_x, ptdf_y, pt_init, pt_p, prho, q, q_init,     &
                q_p, qsws, qswst, qr, qr_p, qrsws, qrswst, rdf, rdf_sc, rho,    &
-               sa, sa_init, sa_p, saswsb, saswst, shf, tend, tend_nr,          &
-               tend_pt, tend_q, tend_qr, te_m, tnr_m, tpt_m, tq_m, tqr_m,      &
-               tsa_m, tswst, tu_m, tv_m, tw_m, u, ug, u_p, v, vg, vpt, v_p,    &
-               w, w_p
+               sa, sa_init, sa_p, saswsb, saswst, shf, tend, te_m, tnr_m,      &
+               tpt_m, tq_m, tqr_m, tsa_m, tswst, tu_m, tv_m, tw_m, u, ug, u_p, &
+               v, vg, vpt, v_p, w, w_p
         
     USE control_parameters,                                                    &
-        ONLY:  cloud_physics, constant_diffusion, cthf, dp_external,           &
+        ONLY:  call_microphysics_at_all_substeps, cloud_physics,               &
+               constant_diffusion, cthf, dp_external,                          &
                dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, humidity,       &
                icloud_scheme, inflow_l, intermediate_timestep_count,           &
@@ -301,4 +306,12 @@
  
        DO  j = nys, nyn
+!
+!--       If required, calculate cloud microphysical impacts (two-moment scheme)
+          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
+               ( intermediate_timestep_count == 1  .OR.                        &
+                 call_microphysics_at_all_substeps )                           &
+             )  THEN
+             CALL microphysics_control( i, j )
+          ENDIF
 !
 !--       Tendency terms for u-velocity component
@@ -481,8 +494,5 @@
              ENDIF
           ENDIF
-!
-!--       If required, calculate tendencies for total water content, liquid water 
-!--       potential temperature, rain water content and rain drop concentration
-          IF ( cloud_physics  .AND.  icloud_scheme == 0 )  CALL microphysics_control( i, j )
+
 !
 !--       If required, compute prognostic equation for potential temperature
@@ -511,11 +521,8 @@
 
 !
-!--          Using microphysical tendencies (latent heat)
-             IF ( cloud_physics )  THEN
-                IF ( icloud_scheme == 0 )  THEN
-                   tend(:,j,i) = tend(:,j,i) + tend_pt(:,j,i)
-                ELSEIF ( icloud_scheme == 1  .AND.  precipitation )  THEN
-                   CALL impact_of_latent_heat( i, j ) 
-                ENDIF
+!--          If required compute impact of latent heat due to precipitation
+             IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.              &
+                  precipitation )  THEN
+                CALL impact_of_latent_heat( i, j ) 
              ENDIF
 
@@ -641,11 +648,9 @@
      
 !
-!--          Using microphysical tendencies
-             IF ( cloud_physics )  THEN
-                IF ( icloud_scheme == 0 )  THEN
-                   tend(:,j,i) = tend(:,j,i) + tend_q(:,j,i)
-                ELSEIF ( icloud_scheme == 1  .AND.  precipitation )  THEN
-                   CALL calc_precipitation( i, j )
-                ENDIF
+!--          If required compute decrease of total water content due to
+!--          precipitation
+             IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.              &
+                  precipitation )  THEN
+                CALL calc_precipitation( i, j )             
              ENDIF
 !
@@ -712,8 +717,4 @@
                 ENDIF
                 CALL diffusion_s( i, j, qr, qrsws, qrswst, wall_qrflux )
-!
-!--             Using microphysical tendencies (autoconversion, accretion, 
-!--             evaporation; if required: sedimentation)
-                tend(:,j,i) = tend(:,j,i) + tend_qr(:,j,i)
 
                 CALL user_actions( i, j, 'qr-tendency' )
@@ -757,9 +758,4 @@
                 ENDIF
                 CALL diffusion_s( i, j, nr, nrsws, nrswst, wall_nrflux )
-!
-!--             Using microphysical tendencies (autoconversion, accretion, 
-!--             selfcollection, breakup, evaporation; 
-!--             if required: sedimentation)
-                tend(:,j,i) = tend(:,j,i) + tend_nr(:,j,i)
 
                 CALL user_actions( i, j, 'nr-tendency' )
@@ -883,4 +879,15 @@
 
 !
+!-- If required, calculate cloud microphysical impacts (two-moment scheme)
+    IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                       &
+         ( intermediate_timestep_count == 1  .OR.                              &
+           call_microphysics_at_all_substeps )                                 &
+       )  THEN
+       CALL cpu_log( log_point(51), 'microphysics', 'start' )
+       CALL microphysics_control
+       CALL cpu_log( log_point(51), 'microphysics', 'stop' )
+    ENDIF
+
+!
 !-- u-velocity component
     CALL cpu_log( log_point(5), 'u-equation', 'start' )
@@ -1156,5 +1163,5 @@
 !
 !--    If required compute impact of latent heat due to precipitation
-       IF ( precipitation )  THEN
+       IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.  precipitation )  THEN
           CALL impact_of_latent_heat
        ENDIF
@@ -1348,5 +1355,5 @@
 !--    If required compute decrease of total water content due to
 !--    precipitation
-       IF ( precipitation )  THEN
+       IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.  precipitation )  THEN
           CALL calc_precipitation
        ENDIF
@@ -1406,4 +1413,160 @@
 
        CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
+
+!
+!--    If required, calculate prognostic equations for rain water content 
+!--    and rain drop concentration
+       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
+
+          CALL cpu_log( log_point(52), 'qr-equation', 'start' )
+
+!
+!--       Calculate prognostic equation for rain water content
+          sbt = tsc(2)
+          IF ( scalar_advec == 'bc-scheme' )  THEN
+
+             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--             Bott-Chlond scheme always uses Euler time step. Thus:
+                sbt = 1.0_wp
+             ENDIF
+             tend = 0.0_wp
+             CALL advec_s_bc( qr, 'qr' )
+
+          ENDIF
+
+!
+!--       qr-tendency terms with no communication
+          IF ( scalar_advec /= 'bc-scheme' )  THEN
+             tend = 0.0_wp
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( ws_scheme_sca )  THEN
+                   CALL advec_s_ws( qr, 'qr' )
+                ELSE
+                   CALL advec_s_pw( qr )
+                ENDIF
+             ELSE
+                CALL advec_s_up( qr )
+             ENDIF
+          ENDIF
+
+          CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux )
+
+          CALL user_actions( 'qr-tendency' )
+
+!
+!--       Prognostic equation for rain water content
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +     &
+                                                     tsc(3) * tqr_m(k,j,i) )   &
+                                           - tsc(5) * rdf_sc(k) * qr(k,j,i)
+                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
+                ENDDO
+             ENDDO
+          ENDDO
+
+!
+!--       Calculate tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb_s_inner(j,i)+1, nzt
+                         tqr_m(k,j,i) = tend(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ELSEIF ( intermediate_timestep_count < &
+                      intermediate_timestep_count_max )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb_s_inner(j,i)+1, nzt
+                         tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * &
+                                                                   tqr_m(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+          ENDIF
+
+          CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
+          CALL cpu_log( log_point(53), 'nr-equation', 'start' )
+
+!
+!--       Calculate prognostic equation for rain drop concentration
+          sbt = tsc(2)
+          IF ( scalar_advec == 'bc-scheme' )  THEN
+
+             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--             Bott-Chlond scheme always uses Euler time step. Thus:
+                sbt = 1.0_wp
+             ENDIF
+             tend = 0.0_wp
+             CALL advec_s_bc( nr, 'nr' )
+
+          ENDIF
+
+!
+!--       nr-tendency terms with no communication
+          IF ( scalar_advec /= 'bc-scheme' )  THEN
+             tend = 0.0_wp
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( ws_scheme_sca )  THEN
+                   CALL advec_s_ws( nr, 'nr' )
+                ELSE
+                   CALL advec_s_pw( nr )
+                ENDIF
+             ELSE
+                CALL advec_s_up( nr )
+             ENDIF
+          ENDIF
+
+          CALL diffusion_s( nr, nrsws, nrswst, wall_nrflux )
+
+          CALL user_actions( 'nr-tendency' )
+
+!
+!--       Prognostic equation for rain drop concentration
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +     &
+                                                     tsc(3) * tnr_m(k,j,i) )   &
+                                           - tsc(5) * rdf_sc(k) * nr(k,j,i)
+                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
+                ENDDO
+             ENDDO
+          ENDDO
+
+!
+!--       Calculate tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb_s_inner(j,i)+1, nzt
+                         tnr_m(k,j,i) = tend(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ELSEIF ( intermediate_timestep_count < &
+                      intermediate_timestep_count_max )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb_s_inner(j,i)+1, nzt
+                         tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * &
+                                                                   tnr_m(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+          ENDIF
+
+          CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
+
+       ENDIF
 
     ENDIF
@@ -1510,5 +1673,4 @@
     ENDIF
 
-
  END SUBROUTINE prognostic_equations_vector
 
@@ -1539,4 +1701,15 @@
           runge_step = 2
        ENDIF
+    ENDIF
+
+!
+!-- If required, calculate cloud microphysical impacts (two-moment scheme)
+    IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                       &
+         ( intermediate_timestep_count == 1  .OR.                              &
+           call_microphysics_at_all_substeps )                                 &
+       )  THEN
+       CALL cpu_log( log_point(51), 'microphysics', 'start' )
+       CALL microphysics_control
+       CALL cpu_log( log_point(51), 'microphysics', 'stop' )
     ENDIF
 
@@ -1791,5 +1964,5 @@
 !
 !--    If required compute impact of latent heat due to precipitation
-       IF ( precipitation )  THEN
+       IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.  precipitation )  THEN
           CALL impact_of_latent_heat
        ENDIF
@@ -1957,5 +2130,5 @@
 !--    If required compute decrease of total water content due to
 !--    precipitation
-       IF ( precipitation )  THEN
+       IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.  precipitation )  THEN
           CALL calc_precipitation
        ENDIF
@@ -1999,4 +2172,127 @@
 
        CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
+
+!
+!--    If required, calculate prognostic equations for rain water content 
+!--    and rain drop concentration
+       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
+
+          CALL cpu_log( log_point(52), 'qr-equation', 'start' )
+!
+!--       qr-tendency terms with communication
+          sbt = tsc(2)
+          IF ( scalar_advec == 'bc-scheme' )  THEN
+
+             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--             Bott-Chlond scheme always uses Euler time step. Thus:
+                sbt = 1.0_wp
+             ENDIF
+             tend = 0.0_wp
+             CALL advec_s_bc( qr, 'qr' )
+
+          ENDIF
+
+!
+!--       qr-tendency terms with no communication
+          IF ( scalar_advec /= 'bc-scheme' )  THEN
+             tend = 0.0_wp
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( ws_scheme_sca )  THEN
+                   CALL advec_s_ws( qr, 'qr' )
+                ELSE
+                   CALL advec_s_pw( qr )
+                ENDIF
+             ELSE
+                CALL advec_s_up( qr )
+             ENDIF
+          ENDIF
+
+          CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux )
+
+          CALL user_actions( 'qr-tendency' )
+
+!
+!--       Prognostic equation for rain water content
+          DO  i = i_left, i_right
+             DO  j = j_south, j_north
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +     &
+                                                       tsc(3) * tqr_m(k,j,i) ) &
+                                           - tsc(5) * rdf_sc(k) * qr(k,j,i)
+                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
+!
+!--                Tendencies for the next Runge-Kutta step
+                   IF ( runge_step == 1 )  THEN
+                      tqr_m(k,j,i) = tend(k,j,i)
+                   ELSEIF ( runge_step == 2 )  THEN
+                      tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp *    &
+                                                                tqr_m(k,j,i)
+                   ENDIF
+                ENDDO
+             ENDDO
+          ENDDO
+
+          CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
+          CALL cpu_log( log_point(53), 'nr-equation', 'start' )
+
+!
+!--       nr-tendency terms with communication
+          sbt = tsc(2)
+          IF ( scalar_advec == 'bc-scheme' )  THEN
+
+             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--             Bott-Chlond scheme always uses Euler time step. Thus:
+                sbt = 1.0_wp
+             ENDIF
+             tend = 0.0_wp
+             CALL advec_s_bc( nr, 'nr' )
+
+          ENDIF
+
+!
+!--       nr-tendency terms with no communication
+          IF ( scalar_advec /= 'bc-scheme' )  THEN
+             tend = 0.0_wp
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( ws_scheme_sca )  THEN
+                   CALL advec_s_ws( nr, 'nr' )
+                ELSE
+                   CALL advec_s_pw( nr )
+                ENDIF
+             ELSE
+                CALL advec_s_up( nr )
+             ENDIF
+          ENDIF
+
+          CALL diffusion_s( nr, nrsws, nrswst, wall_nrflux )
+
+          CALL user_actions( 'nr-tendency' )
+
+!
+!--       Prognostic equation for rain drop concentration
+          DO  i = i_left, i_right
+             DO  j = j_south, j_north
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +     &
+                                                       tsc(3) * tnr_m(k,j,i) ) &
+                                           - tsc(5) * rdf_sc(k) * nr(k,j,i)
+                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
+!
+!--                Tendencies for the next Runge-Kutta step
+                   IF ( runge_step == 1 )  THEN
+                      tnr_m(k,j,i) = tend(k,j,i)
+                   ELSEIF ( runge_step == 2 )  THEN
+                      tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp *    &
+                                                                tnr_m(k,j,i)
+                   ENDIF
+                ENDDO
+             ENDDO
+          ENDDO
+
+          CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
+
+       ENDIF
 
     ENDIF
@@ -2095,5 +2391,4 @@
     ENDIF
 
-
  END SUBROUTINE prognostic_equations_acc
 
