Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 3869)
+++ palm/trunk/SOURCE/Makefile	(revision 3870)
@@ -776,5 +776,4 @@
 	modules.o
 advec_ws.o: \
-	bulk_cloud_model_mod.o \
 	mod_kinds.o \
 	modules.o \
@@ -815,4 +814,9 @@
 	basic_constants_and_equations_mod.o \
 	cpulog_mod.o \
+	advec_ws.o \
+	advec_s_bc.o \
+	advec_s_pw.o \
+	advec_s_up.o \
+	diffusion_s.o \
 	mod_kinds.o \
 	modules.o \
Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 3869)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 3870)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Moving initialization of bcm specific flux arrays into bulk_cloud_model_mod
+! 
+! 3864 2019-04-05 09:01:56Z monakurppa
 ! Remove tailing white spaces
 ! 
@@ -325,14 +328,14 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  diss_l_diss, diss_l_e, diss_l_nc, diss_l_nr, diss_l_pt,      &
-                  diss_l_q, diss_l_qc, diss_l_qr, diss_l_s, diss_l_sa,         &
+           ONLY:  diss_l_diss, diss_l_e, diss_l_pt,      &
+                  diss_l_q, diss_l_s, diss_l_sa,         &
                   diss_l_u, diss_l_v, diss_l_w, flux_l_diss, flux_l_e,         & 
-                  flux_l_nc, flux_l_nr, flux_l_pt, flux_l_q, flux_l_qc,        &
-                  flux_l_qr, flux_l_s, flux_l_sa, flux_l_u, flux_l_v,          &
-                  flux_l_w, diss_s_diss, diss_s_e, diss_s_nc,  diss_s_nr,      &
-                  diss_s_pt, diss_s_q, diss_s_qc, diss_s_qr, diss_s_s,         &
+                  flux_l_pt, flux_l_q,      &
+                  flux_l_s, flux_l_sa, flux_l_u, flux_l_v,          &
+                  flux_l_w, diss_s_diss, diss_s_e,      &
+                  diss_s_pt, diss_s_q, diss_s_s,         &
                   diss_s_sa, diss_s_u, diss_s_v,  diss_s_w, flux_s_diss,       &
-                  flux_s_e, flux_s_nc, flux_s_nr, flux_s_pt, flux_s_q,         &
-                  flux_s_qc, flux_s_qr, flux_s_s, flux_s_sa, flux_s_u,         &
+                  flux_s_e, flux_s_pt, flux_s_q,         &
+                  flux_s_s, flux_s_sa, flux_s_u,         &
                   flux_s_v, flux_s_w
 
@@ -346,7 +349,4 @@
        USE kinds
 
-       USE bulk_cloud_model_mod,                                               &
-           ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
-
        USE pegrid
 
@@ -355,7 +355,7 @@
 
        USE statistics,                                                         &
-           ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsncs_ws_l,& 
-                  sums_wsnrs_ws_l,sums_wspts_ws_l, sums_wsqcs_ws_l,            &
-                  sums_wsqrs_ws_l, sums_wsqs_ws_l, sums_wsss_ws_l,             &
+           ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l,&
+                  sums_wspts_ws_l,            &
+                  sums_wsqs_ws_l, sums_wsss_ws_l,             &
                   sums_wssas_ws_l,  sums_wsss_ws_l, sums_wsus_ws_l,            &
                   sums_wsvs_ws_l
@@ -397,22 +397,8 @@
              sums_wsqs_ws_l = 0.0_wp
           ENDIF
-          
+
           IF ( passive_scalar )  THEN
              ALLOCATE( sums_wsss_ws_l(nzb:nzt+1,0:threads_per_task-1) )
              sums_wsss_ws_l = 0.0_wp
-          ENDIF
-
-          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
-             ALLOCATE( sums_wsqcs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
-             ALLOCATE( sums_wsncs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
-             sums_wsqcs_ws_l = 0.0_wp
-             sums_wsncs_ws_l = 0.0_wp
-          ENDIF
-
-          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
-             ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
-             ALLOCATE( sums_wsnrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
-             sums_wsqrs_ws_l = 0.0_wp
-             sums_wsnrs_ws_l = 0.0_wp
           ENDIF
 
@@ -477,5 +463,5 @@
                           diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
              ENDIF
-             
+
              IF ( passive_scalar )  THEN
                 ALLOCATE( flux_s_s(nzb+1:nzt,0:threads_per_task-1),            &
@@ -483,26 +469,4 @@
                 ALLOCATE( flux_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
                           diss_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
-             ENDIF
-
-             IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
-                ALLOCATE( flux_s_qc(nzb+1:nzt,0:threads_per_task-1),           &
-                          diss_s_qc(nzb+1:nzt,0:threads_per_task-1),           &
-                          flux_s_nc(nzb+1:nzt,0:threads_per_task-1),           &
-                          diss_s_nc(nzb+1:nzt,0:threads_per_task-1) )
-                ALLOCATE( flux_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
-                          diss_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
-                          flux_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
-                          diss_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) 
-             ENDIF                  
-
-             IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
-                ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1),           &
-                          diss_s_qr(nzb+1:nzt,0:threads_per_task-1),           &
-                          flux_s_nr(nzb+1:nzt,0:threads_per_task-1),           &
-                          diss_s_nr(nzb+1:nzt,0:threads_per_task-1) )
-                ALLOCATE( flux_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
-                          diss_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
-                          flux_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
-                          diss_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) 
              ENDIF
 
@@ -1140,8 +1104,5 @@
                   ws_scheme_sca, salsa
 
-       USE kinds  
-
-       USE bulk_cloud_model_mod,                                               &
-           ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
+       USE kinds
 
        USE salsa_util_mod,                                                     &
@@ -1149,7 +1110,7 @@
 
        USE statistics,                                                         &
-           ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsncs_ws_l,&
-                  sums_wsnrs_ws_l, sums_wspts_ws_l, sums_wsqcs_ws_l,           &
-                  sums_wsqrs_ws_l, sums_wsqs_ws_l, sums_wsss_ws_l,             &
+           ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, &
+                  sums_wspts_ws_l,           &
+                  sums_wsqs_ws_l, sums_wsss_ws_l,             &
                   sums_wssas_ws_l, sums_wsus_ws_l, sums_wsvs_ws_l
                    
@@ -1177,12 +1138,4 @@
           IF ( humidity       )  sums_wsqs_ws_l = 0.0_wp
           IF ( passive_scalar )  sums_wsss_ws_l = 0.0_wp
-          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
-             sums_wsqcs_ws_l = 0.0_wp
-             sums_wsncs_ws_l = 0.0_wp
-          ENDIF
-          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
-             sums_wsqrs_ws_l = 0.0_wp
-             sums_wsnrs_ws_l = 0.0_wp
-          ENDIF
           IF ( ocean_mode )  sums_wssas_ws_l = 0.0_wp
           IF ( salsa )       sums_salsa_ws_l = 0.0_wp
Index: palm/trunk/SOURCE/bulk_cloud_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/bulk_cloud_model_mod.f90	(revision 3869)
+++ palm/trunk/SOURCE/bulk_cloud_model_mod.f90	(revision 3870)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Moving prognostic equations of bcm into bulk_cloud_model_mod
+! 
+! 3869 2019-04-08 11:54:20Z knoop
 ! moving the furniture around ;-)
 ! 
@@ -199,4 +202,17 @@
  MODULE bulk_cloud_model_mod
 
+
+    USE advec_s_bc_mod,                                                        &
+        ONLY:  advec_s_bc
+
+    USE advec_s_pw_mod,                                                        &
+        ONLY:  advec_s_pw
+
+    USE advec_s_up_mod,                                                        &
+        ONLY:  advec_s_up
+
+    USE advec_ws,                                                              &
+        ONLY:  advec_s_ws
+
     USE arrays_3d,                                                             &
         ONLY:  ddzu, diss, dzu, dzw, hyp, hyrho,                               &
@@ -204,5 +220,9 @@
                precipitation_amount, prr, pt, d_exner, pt_init, q, ql, ql_1,   &
                qc, qc_1, qc_2, qc_3, qc_p, qr, qr_1, qr_2, qr_3, qr_p,         &
-               exner, zu, tnc_m, tnr_m, tqc_m, tqr_m
+               exner, zu, tnc_m, tnr_m, tqc_m, tqr_m, tend, rdf_sc, &
+               flux_l_qc, flux_l_qr, flux_l_nc, flux_l_nr, &
+               flux_s_qc, flux_s_qr, flux_s_nc, flux_s_nr, &
+               diss_l_qc, diss_l_qr, diss_l_nc, diss_l_nr, &
+               diss_s_qc, diss_s_qr, diss_s_nc, diss_s_nr
 
     USE averaging,                                                             &
@@ -219,8 +239,12 @@
                intermediate_timestep_count_max, large_scale_forcing,           &
                lsf_surf, pt_surface, rho_surface, surface_pressure,            &
-               time_do2d_xy, message_string, initializing_actions
+               time_do2d_xy, message_string, initializing_actions,             &
+               ws_scheme_sca, scalar_advec, timestep_scheme, tsc, loop_optimization
 
     USE cpulog,                                                                &
-        ONLY:  cpu_log, log_point_s
+        ONLY:  cpu_log, log_point, log_point_s
+
+    USE diffusion_s_mod,                                                       &
+        ONLY:  diffusion_s
 
     USE grid_variables,                                                        &
@@ -231,12 +255,16 @@
                wall_flags_0
 
-    USE  kinds
+    USE kinds
+
+    USE pegrid,                                                                &
+        ONLY:  threads_per_task
 
     USE statistics,                                                            &
-        ONLY:  weight_pres, weight_substep
+        ONLY:  weight_pres, weight_substep, sums_wsncs_ws_l, sums_wsnrs_ws_l, sums_wsqcs_ws_l, sums_wsqrs_ws_l
 
     USE surface_mod,                                                           &
         ONLY :  bc_h, get_topography_top_index_ji, surf_bulk_cloud_model,      &
-                surf_microphysics_morrison, surf_microphysics_seifert
+                surf_microphysics_morrison, surf_microphysics_seifert, &
+                surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
 
     IMPLICIT NONE
@@ -332,4 +360,6 @@
            bcm_header, &
            bcm_actions, &
+           bcm_actions_micro, &
+           bcm_prognostic_equations, &
            bcm_3d_data_averaging, &
            bcm_data_output_2d, &
@@ -389,4 +419,14 @@
        MODULE PROCEDURE bcm_actions_ij
     END INTERFACE bcm_actions
+
+    INTERFACE bcm_actions_micro
+       MODULE PROCEDURE bcm_actions_micro
+       MODULE PROCEDURE bcm_actions_micro_ij
+    END INTERFACE bcm_actions_micro
+
+    INTERFACE bcm_prognostic_equations
+       MODULE PROCEDURE bcm_prognostic_equations
+       MODULE PROCEDURE bcm_prognostic_equations_ij
+    END INTERFACE bcm_prognostic_equations
 
     INTERFACE bcm_swap_timelevel
@@ -807,4 +847,57 @@
        ENDIF
 
+       IF ( ws_scheme_sca )  THEN
+
+          IF ( microphysics_morrison )  THEN
+             ALLOCATE( sums_wsqcs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
+             ALLOCATE( sums_wsncs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
+             sums_wsqcs_ws_l = 0.0_wp
+             sums_wsncs_ws_l = 0.0_wp
+          ENDIF
+
+          IF ( microphysics_seifert )  THEN
+             ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
+             ALLOCATE( sums_wsnrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
+             sums_wsqrs_ws_l = 0.0_wp
+             sums_wsnrs_ws_l = 0.0_wp
+          ENDIF
+
+       ENDIF
+
+!
+!--    Arrays needed for reasons of speed optimization for cache version.
+!--    For the vector version the buffer arrays are not necessary,
+!--    because the the fluxes can swapped directly inside the loops of the
+!--    advection routines.
+       IF ( loop_optimization /= 'vector' )  THEN
+
+          IF ( ws_scheme_sca )  THEN
+
+             IF ( microphysics_morrison )  THEN
+                ALLOCATE( flux_s_qc(nzb+1:nzt,0:threads_per_task-1),           &
+                          diss_s_qc(nzb+1:nzt,0:threads_per_task-1),           &
+                          flux_s_nc(nzb+1:nzt,0:threads_per_task-1),           &
+                          diss_s_nc(nzb+1:nzt,0:threads_per_task-1) )
+                ALLOCATE( flux_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
+                          diss_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
+                          flux_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
+                          diss_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
+             ENDIF
+
+             IF ( microphysics_seifert )  THEN
+                ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1),           &
+                          diss_s_qr(nzb+1:nzt,0:threads_per_task-1),           &
+                          flux_s_nr(nzb+1:nzt,0:threads_per_task-1),           &
+                          diss_s_nr(nzb+1:nzt,0:threads_per_task-1) )
+                ALLOCATE( flux_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
+                          diss_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
+                          flux_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
+                          diss_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
+             ENDIF
+
+          ENDIF
+
+       ENDIF
+
 !
 !--    Initial assignment of the pointers
@@ -1004,5 +1097,82 @@
 !> Control of microphysics for all grid points
 !------------------------------------------------------------------------------!
-    SUBROUTINE bcm_actions
+    SUBROUTINE bcm_actions( location )
+
+
+    CHARACTER (LEN=*), INTENT(IN) ::  location !< call location string
+
+    SELECT CASE ( location )
+
+       CASE ( 'before_timestep' )
+
+          IF ( ws_scheme_sca )  THEN
+
+             IF ( microphysics_morrison )  THEN
+                sums_wsqcs_ws_l = 0.0_wp
+                sums_wsncs_ws_l = 0.0_wp
+             ENDIF
+             IF ( microphysics_seifert )  THEN
+                sums_wsqrs_ws_l = 0.0_wp
+                sums_wsnrs_ws_l = 0.0_wp
+             ENDIF
+
+          ENDIF
+
+       CASE DEFAULT
+          CONTINUE
+
+    END SELECT
+
+    END SUBROUTINE bcm_actions
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Control of microphysics for grid points i,j
+!------------------------------------------------------------------------------!
+
+    SUBROUTINE bcm_actions_ij( i, j, location )
+
+
+    INTEGER(iwp),      INTENT(IN) ::  i         !< grid index in x-direction
+    INTEGER(iwp),      INTENT(IN) ::  j         !< grid index in y-direction
+    CHARACTER (LEN=*), INTENT(IN) ::  location  !< call location string
+    INTEGER(iwp)  ::  dummy  !< call location string
+
+    IF ( bulk_cloud_model    )   dummy = i + j
+
+    SELECT CASE ( location )
+
+       CASE ( 'before_timestep' )
+
+          IF ( ws_scheme_sca )  THEN
+
+             IF ( microphysics_morrison )  THEN
+                sums_wsqcs_ws_l = 0.0_wp
+                sums_wsncs_ws_l = 0.0_wp
+             ENDIF
+             IF ( microphysics_seifert )  THEN
+                sums_wsqrs_ws_l = 0.0_wp
+                sums_wsnrs_ws_l = 0.0_wp
+             ENDIF
+
+          ENDIF
+
+       CASE DEFAULT
+          CONTINUE
+
+    END SELECT
+
+
+    END SUBROUTINE bcm_actions_ij
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Control of microphysics for all grid points
+!------------------------------------------------------------------------------!
+    SUBROUTINE bcm_actions_micro
 
        IMPLICIT NONE
@@ -1055,5 +1225,5 @@
        CALL calc_precipitation_amount
 
-    END SUBROUTINE bcm_actions
+    END SUBROUTINE bcm_actions_micro
 
 
@@ -1064,5 +1234,5 @@
 !------------------------------------------------------------------------------!
 
-    SUBROUTINE bcm_actions_ij( i, j )
+    SUBROUTINE bcm_actions_micro_ij( i, j )
 
        IMPLICIT NONE
@@ -1117,5 +1287,632 @@
        CALL calc_precipitation_amount_ij( i,j )
 
-    END SUBROUTINE bcm_actions_ij
+    END SUBROUTINE bcm_actions_micro_ij
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Control of microphysics for all grid points
+!------------------------------------------------------------------------------!
+    SUBROUTINE bcm_prognostic_equations
+
+
+       INTEGER(iwp) ::  i         !< grid index in x-direction
+       INTEGER(iwp) ::  j         !< grid index in y-direction
+       INTEGER(iwp) ::  k         !< grid index in z-direction
+
+       REAL(wp)     ::  sbt  !<
+
+!
+!--    If required, calculate prognostic equations for cloud water content
+!--    and cloud drop concentration
+       IF ( microphysics_morrison )  THEN
+
+          CALL cpu_log( log_point(67), 'qc-equation', 'start' )
+
+!
+!--       Calculate prognostic equation for cloud 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( qc, 'qc' )
+
+          ENDIF
+
+!
+!--       qc-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( qc, 'qc' )
+                ELSE
+                   CALL advec_s_pw( qc )
+                ENDIF
+             ELSE
+                CALL advec_s_up( qc )
+             ENDIF
+          ENDIF
+
+          CALL diffusion_s( qc,                                                &
+                            surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,          &
+                            surf_def_h(2)%qcsws,                               &
+                            surf_lsm_h%qcsws,    surf_usm_h%qcsws,             &
+                            surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,          &
+                            surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,          &
+                            surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,          &
+                            surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,          &
+                            surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,          &
+                            surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
+
+!
+!--       Prognostic equation for cloud water content
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+                   qc_p(k,j,i) = qc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
+                                                      tsc(3) * tqc_m(k,j,i) )  &
+                                                    - tsc(5) * rdf_sc(k) *     &
+                                                               qc(k,j,i)       &
+                                             )                                 &
+                                    * MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 0 )   &
+                                          )
+                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_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+1, nzt
+                         tqc_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+1, nzt
+                         tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
+                                         + 5.3125_wp * tqc_m(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+          ENDIF
+
+          CALL cpu_log( log_point(67), 'qc-equation', 'stop' )
+
+          CALL cpu_log( log_point(68), 'nc-equation', 'start' )
+!
+!--       Calculate prognostic equation for cloud 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( nc, 'nc' )
+
+          ENDIF
+
+!
+!--       nc-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( nc, 'nc' )
+                ELSE
+                   CALL advec_s_pw( nc )
+                ENDIF
+             ELSE
+                CALL advec_s_up( nc )
+             ENDIF
+          ENDIF
+
+          CALL diffusion_s( nc,                                                &
+                            surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,          &
+                            surf_def_h(2)%ncsws,                               &
+                            surf_lsm_h%ncsws,    surf_usm_h%ncsws,             &
+                            surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,          &
+                            surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,          &
+                            surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,          &
+                            surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,          &
+                            surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,          &
+                            surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
+
+!
+!--       Prognostic equation for cloud drop concentration
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+                   nc_p(k,j,i) = nc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
+                                                      tsc(3) * tnc_m(k,j,i) )  &
+                                                    - tsc(5) * rdf_sc(k) *     &
+                                                               nc(k,j,i)       &
+                                             )                                 &
+                                   * MERGE( 1.0_wp, 0.0_wp,                    &
+                                             BTEST( wall_flags_0(k,j,i), 0 )   &
+                                          )
+                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_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+1, nzt
+                         tnc_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+1, nzt
+                         tnc_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
+                                         + 5.3125_wp * tnc_m(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+          ENDIF
+
+          CALL cpu_log( log_point(68), 'nc-equation', 'stop' )
+
+       ENDIF
+!
+!--    If required, calculate prognostic equations for rain water content
+!--    and rain drop concentration
+       IF ( microphysics_seifert )  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,                                                &
+                            surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,          &
+                            surf_def_h(2)%qrsws,                               &
+                            surf_lsm_h%qrsws,    surf_usm_h%qrsws,             &
+                            surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,          &
+                            surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,          &
+                            surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,          &
+                            surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,          &
+                            surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,          &
+                            surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
+
+!
+!--       Prognostic equation for rain water content
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb+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)       &
+                                             )                                 &
+                                    * MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 0 )   &
+                                          )
+                   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+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+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,                                                &
+                            surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,          &
+                            surf_def_h(2)%nrsws,                               &
+                            surf_lsm_h%nrsws,    surf_usm_h%nrsws,             &
+                            surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,          &
+                            surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,          &
+                            surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,          &
+                            surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,          &
+                            surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,          &
+                            surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
+
+!
+!--       Prognostic equation for rain drop concentration
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb+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)       &
+                                             )                                 &
+                                   * MERGE( 1.0_wp, 0.0_wp,                    &
+                                             BTEST( wall_flags_0(k,j,i), 0 )   &
+                                          )
+                   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+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+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
+
+    END SUBROUTINE bcm_prognostic_equations
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Control of microphysics for grid points i,j
+!------------------------------------------------------------------------------!
+
+    SUBROUTINE bcm_prognostic_equations_ij( i, j, i_omp_start, tn )
+
+
+       INTEGER(iwp), INTENT(IN) ::  i            !< grid index in x-direction
+       INTEGER(iwp), INTENT(IN) ::  j            !< grid index in y-direction
+       INTEGER(iwp)             ::  k            !< grid index in z-direction
+       INTEGER(iwp), INTENT(IN) ::  i_omp_start  !< first loop index of i-loop in prognostic_equations
+       INTEGER(iwp), INTENT(IN) ::  tn           !< task number of openmp task
+
+!
+!--    If required, calculate prognostic equations for cloud water content
+!--    and cloud drop concentration
+       IF ( microphysics_morrison )  THEN
+!
+!--       Calculate prognostic equation for cloud water content
+          tend(:,j,i) = 0.0_wp
+          IF ( timestep_scheme(1:5) == 'runge' ) &
+          THEN
+             IF ( ws_scheme_sca )  THEN
+                CALL advec_s_ws( i, j, qc, 'qc', flux_s_qc,       &
+                                 diss_s_qc, flux_l_qc, diss_l_qc, &
+                                 i_omp_start, tn )
+             ELSE
+                CALL advec_s_pw( i, j, qc )
+             ENDIF
+          ELSE
+             CALL advec_s_up( i, j, qc )
+          ENDIF
+          CALL diffusion_s( i, j, qc,                                   &
+                            surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,   &
+                            surf_def_h(2)%qcsws,                        &
+                            surf_lsm_h%qcsws,    surf_usm_h%qcsws,      &
+                            surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,   &
+                            surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,   &
+                            surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,   &
+                            surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,   &
+                            surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,   &
+                            surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
+
+!
+!--       Prognostic equation for cloud water content
+          DO  k = nzb+1, nzt
+             qc_p(k,j,i) = qc(k,j,i) + ( dt_3d *                         &
+                                                ( tsc(2) * tend(k,j,i) + &
+                                                  tsc(3) * tqc_m(k,j,i) )&
+                                                - tsc(5) * rdf_sc(k)     &
+                                                         * qc(k,j,i)     &
+                                       )                                 &
+                                 * MERGE( 1.0_wp, 0.0_wp,                &
+                                          BTEST( wall_flags_0(k,j,i), 0 )&
+                                        )
+             IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
+          ENDDO
+!
+!--       Calculate tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                DO  k = nzb+1, nzt
+                   tqc_m(k,j,i) = tend(k,j,i)
+                ENDDO
+             ELSEIF ( intermediate_timestep_count < &
+                      intermediate_timestep_count_max )  THEN
+                DO  k = nzb+1, nzt
+                   tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
+                                     5.3125_wp * tqc_m(k,j,i)
+                ENDDO
+             ENDIF
+          ENDIF
+
+!
+!--       Calculate prognostic equation for cloud drop concentration.
+          tend(:,j,i) = 0.0_wp
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( ws_scheme_sca )  THEN
+                CALL advec_s_ws( i, j, nc, 'nc', flux_s_nc,    &
+                              diss_s_nc, flux_l_nc, diss_l_nc, &
+                              i_omp_start, tn )
+             ELSE
+                CALL advec_s_pw( i, j, nc )
+             ENDIF
+          ELSE
+             CALL advec_s_up( i, j, nc )
+          ENDIF
+          CALL diffusion_s( i, j, nc,                                    &
+                            surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,    &
+                            surf_def_h(2)%ncsws,                         &
+                            surf_lsm_h%ncsws,    surf_usm_h%ncsws,       &
+                            surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,    &
+                            surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,    &
+                            surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,    &
+                            surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,    &
+                            surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,    &
+                            surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
+
+!
+!--       Prognostic equation for cloud drop concentration
+          DO  k = nzb+1, nzt
+             nc_p(k,j,i) = nc(k,j,i) + ( dt_3d *                         &
+                                                ( tsc(2) * tend(k,j,i) + &
+                                                  tsc(3) * tnc_m(k,j,i) )&
+                                                - tsc(5) * rdf_sc(k)     &
+                                                         * nc(k,j,i)     &
+                                       )                                 &
+                                 * MERGE( 1.0_wp, 0.0_wp,                &
+                                          BTEST( wall_flags_0(k,j,i), 0 )&
+                                        )
+             IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
+          ENDDO
+!
+!--       Calculate tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                DO  k = nzb+1, nzt
+                   tnc_m(k,j,i) = tend(k,j,i)
+                ENDDO
+             ELSEIF ( intermediate_timestep_count < &
+                      intermediate_timestep_count_max )  THEN
+                DO  k = nzb+1, nzt
+                   tnc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
+                                     5.3125_wp * tnc_m(k,j,i)
+                ENDDO
+             ENDIF
+          ENDIF
+
+       ENDIF
+!
+!--    If required, calculate prognostic equations for rain water content
+!--    and rain drop concentration
+       IF ( microphysics_seifert )  THEN
+!
+!--             Calculate prognostic equation for rain water content
+          tend(:,j,i) = 0.0_wp
+          IF ( timestep_scheme(1:5) == 'runge' ) &
+          THEN
+             IF ( ws_scheme_sca )  THEN
+                CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr,       &
+                                 diss_s_qr, flux_l_qr, diss_l_qr, &
+                                 i_omp_start, tn )
+             ELSE
+                CALL advec_s_pw( i, j, qr )
+             ENDIF
+          ELSE
+             CALL advec_s_up( i, j, qr )
+          ENDIF
+          CALL diffusion_s( i, j, qr,                                   &
+                            surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,   &
+                            surf_def_h(2)%qrsws,                        &
+                            surf_lsm_h%qrsws,    surf_usm_h%qrsws,      &
+                            surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,   &
+                            surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,   &
+                            surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,   &
+                            surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,   &
+                            surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,   &
+                            surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
+
+!
+!--       Prognostic equation for rain water content
+          DO  k = nzb+1, nzt
+             qr_p(k,j,i) = qr(k,j,i) + ( dt_3d *                         &
+                                                ( tsc(2) * tend(k,j,i) + &
+                                                  tsc(3) * tqr_m(k,j,i) )&
+                                                - tsc(5) * rdf_sc(k)     &
+                                                         * qr(k,j,i)     &
+                                       )                                 &
+                                 * MERGE( 1.0_wp, 0.0_wp,                &
+                                          BTEST( wall_flags_0(k,j,i), 0 )&
+                                        )
+             IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
+          ENDDO
+!
+!--       Calculate tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                DO  k = nzb+1, nzt
+                   tqr_m(k,j,i) = tend(k,j,i)
+                ENDDO
+             ELSEIF ( intermediate_timestep_count < &
+                      intermediate_timestep_count_max )  THEN
+                DO  k = nzb+1, nzt
+                   tqr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
+                                     5.3125_wp * tqr_m(k,j,i)
+                ENDDO
+             ENDIF
+          ENDIF
+
+!
+!--       Calculate prognostic equation for rain drop concentration.
+          tend(:,j,i) = 0.0_wp
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( ws_scheme_sca )  THEN
+                CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr,    &
+                              diss_s_nr, flux_l_nr, diss_l_nr, &
+                              i_omp_start, tn )
+             ELSE
+                CALL advec_s_pw( i, j, nr )
+             ENDIF
+          ELSE
+             CALL advec_s_up( i, j, nr )
+          ENDIF
+          CALL diffusion_s( i, j, nr,                                    &
+                            surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,    &
+                            surf_def_h(2)%nrsws,                         &
+                            surf_lsm_h%nrsws,    surf_usm_h%nrsws,       &
+                            surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,    &
+                            surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,    &
+                            surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,    &
+                            surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,    &
+                            surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,    &
+                            surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
+
+!
+!--       Prognostic equation for rain drop concentration
+          DO  k = nzb+1, nzt
+             nr_p(k,j,i) = nr(k,j,i) + ( dt_3d *                         &
+                                                ( tsc(2) * tend(k,j,i) + &
+                                                  tsc(3) * tnr_m(k,j,i) )&
+                                                - tsc(5) * rdf_sc(k)     &
+                                                         * nr(k,j,i)     &
+                                       )                                 &
+                                 * MERGE( 1.0_wp, 0.0_wp,                &
+                                          BTEST( wall_flags_0(k,j,i), 0 )&
+                                        )
+             IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
+          ENDDO
+!
+!--       Calculate tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                DO  k = nzb+1, nzt
+                   tnr_m(k,j,i) = tend(k,j,i)
+                ENDDO
+             ELSEIF ( intermediate_timestep_count < &
+                      intermediate_timestep_count_max )  THEN
+                DO  k = nzb+1, nzt
+                   tnr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
+                                     5.3125_wp * tnr_m(k,j,i)
+                ENDDO
+             ENDIF
+          ENDIF
+
+       ENDIF
+
+    END SUBROUTINE bcm_prognostic_equations_ij
 
 
Index: palm/trunk/SOURCE/module_interface.f90
===================================================================
--- palm/trunk/SOURCE/module_interface.f90	(revision 3869)
+++ palm/trunk/SOURCE/module_interface.f90	(revision 3870)
@@ -146,4 +146,6 @@
                bcm_init,                                                       &
                bcm_header,                                                     &
+               bcm_actions,                                                    &
+               bcm_prognostic_equations,                                       &
                bcm_swap_timelevel,                                             &
                bcm_3d_data_averaging,                                          &
@@ -876,4 +878,5 @@
 
 
+    IF ( bulk_cloud_model    )  CALL bcm_actions( location )
     IF ( gust_module_enabled )  CALL gust_actions( location )
     IF ( user_module_enabled )  CALL user_actions( location )
@@ -896,4 +899,5 @@
 
 
+    IF ( bulk_cloud_model    )  CALL bcm_actions( i, j, location )
     IF ( gust_module_enabled )  CALL gust_actions( i, j, location )
     IF ( user_module_enabled )  CALL user_actions( i, j, location )
@@ -911,4 +915,5 @@
 
 
+    IF ( bulk_cloud_model    )  CALL bcm_prognostic_equations()
     IF ( gust_module_enabled )  CALL gust_prognostic_equations()
     IF ( ocean_mode          )  CALL ocean_prognostic_equations()
@@ -933,4 +938,5 @@
 
 
+    IF ( bulk_cloud_model    )  CALL bcm_prognostic_equations( i, j, i_omp_start, tn )
     IF ( gust_module_enabled )  CALL gust_prognostic_equations( i, j, i_omp_start, tn )
     IF ( ocean_mode          )  CALL ocean_prognostic_equations( i, j, i_omp_start, tn )
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 3869)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 3870)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Moving prognostic equations of bcm into bulk_cloud_model_mod
+! 
+! 3864 2019-04-05 09:01:56Z monakurppa
 ! Modifications made for salsa:
 ! - salsa_prognostic_equations moved to salsa_mod (and the call to 
@@ -332,5 +335,4 @@
  MODULE prognostic_equations_mod
 
-
     USE advec_s_bc_mod,                                                        &
         ONLY:  advec_s_bc
@@ -364,20 +366,19 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  diss_l_e, diss_l_nc, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qc, &
-               diss_l_qr, diss_l_s, diss_l_sa, diss_s_e, diss_s_nc, diss_s_nr, &
-               diss_s_pt, diss_s_q, diss_s_qc, diss_s_qr, diss_s_s, diss_s_sa, &
-               e, e_p, flux_s_e, flux_s_nc, flux_s_nr, flux_s_pt, flux_s_q,    &
-               flux_s_qc, flux_s_qr, flux_s_s, flux_s_sa, flux_l_e, flux_l_nc, &
-               flux_l_nr, flux_l_pt, flux_l_q, flux_l_qc, flux_l_qr, flux_l_s, &
-               flux_l_sa, nc, nc_p, nr, nr_p, pt, ptdf_x, ptdf_y, pt_init,     &
-               pt_p, prho, q, q_init, q_p, qc, qc_p, qr, qr_p, rdf, rdf_sc,    &
-               ref_state, rho_ocean, s, s_init, s_p, tend, te_m, tnc_m,        &
-               tnr_m, tpt_m, tq_m, tqc_m, tqr_m, ts_m, tu_m, tv_m, tw_m, u,    &
+        ONLY:  diss_l_e, diss_l_pt, diss_l_q,                                  &
+               diss_l_s, diss_l_sa, diss_s_e,                                  &
+               diss_s_pt, diss_s_q, diss_s_s, diss_s_sa,                       &
+               e, e_p, flux_s_e, flux_s_pt, flux_s_q,                          &
+               flux_s_s, flux_s_sa, flux_l_e,                                  &
+               flux_l_pt, flux_l_q, flux_l_s,                                  &
+               flux_l_sa, pt, ptdf_x, ptdf_y, pt_init,                         &
+               pt_p, prho, q, q_init, q_p, rdf, rdf_sc,                        &
+               ref_state, rho_ocean, s, s_init, s_p, tend, te_m,               &
+               tpt_m, tq_m, ts_m, tu_m, tv_m, tw_m, u,                         &
                ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
 
     USE bulk_cloud_model_mod,                                                  &
         ONLY:  call_microphysics_at_all_substeps, bulk_cloud_model,            &
-               bcm_actions, microphysics_sat_adjust,                           &
-               microphysics_morrison, microphysics_seifert
+               bcm_actions_micro, microphysics_sat_adjust
 
     USE buoyancy_mod,                                                          &
@@ -663,5 +664,5 @@
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
-             CALL bcm_actions( i, j )
+             CALL bcm_actions_micro( i, j )
            ENDDO
        ENDDO
@@ -1104,237 +1105,4 @@
              ENDIF
 
-!
-!--          If required, calculate prognostic equations for cloud water content
-!--          and cloud drop concentration
-             IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
-!
-!--             Calculate prognostic equation for cloud water content
-                tend(:,j,i) = 0.0_wp
-                IF ( timestep_scheme(1:5) == 'runge' ) &
-                THEN
-                   IF ( ws_scheme_sca )  THEN
-                      CALL advec_s_ws( i, j, qc, 'qc', flux_s_qc,       &
-                                       diss_s_qc, flux_l_qc, diss_l_qc, &
-                                       i_omp_start, tn )
-                   ELSE
-                      CALL advec_s_pw( i, j, qc )
-                   ENDIF
-                ELSE
-                   CALL advec_s_up( i, j, qc )
-                ENDIF
-                CALL diffusion_s( i, j, qc,                                   &
-                                  surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,   &
-                                  surf_def_h(2)%qcsws,                        &
-                                  surf_lsm_h%qcsws,    surf_usm_h%qcsws,      &  
-                                  surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,   &
-                                  surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,   &
-                                  surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,   &
-                                  surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,   &
-                                  surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,   &
-                                  surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
-
-!
-!--             Prognostic equation for cloud water content
-                DO  k = nzb+1, nzt
-                   qc_p(k,j,i) = qc(k,j,i) + ( dt_3d *                         &
-                                                      ( tsc(2) * tend(k,j,i) + &
-                                                        tsc(3) * tqc_m(k,j,i) )&
-                                                      - tsc(5) * rdf_sc(k)     &
-                                                               * qc(k,j,i)     &
-                                             )                                 &
-                                       * MERGE( 1.0_wp, 0.0_wp,                &
-                                                BTEST( wall_flags_0(k,j,i), 0 )&
-                                              ) 
-                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
-                ENDDO
-!
-!--             Calculate tendencies for the next Runge-Kutta step
-                IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                   IF ( intermediate_timestep_count == 1 )  THEN
-                      DO  k = nzb+1, nzt
-                         tqc_m(k,j,i) = tend(k,j,i)
-                      ENDDO
-                   ELSEIF ( intermediate_timestep_count < &
-                            intermediate_timestep_count_max )  THEN
-                      DO  k = nzb+1, nzt
-                         tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
-                                           5.3125_wp * tqc_m(k,j,i)
-                      ENDDO
-                   ENDIF
-                ENDIF
-
-!
-!--             Calculate prognostic equation for cloud drop concentration.
-                tend(:,j,i) = 0.0_wp
-                IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                   IF ( ws_scheme_sca )  THEN
-                      CALL advec_s_ws( i, j, nc, 'nc', flux_s_nc,    &
-                                    diss_s_nc, flux_l_nc, diss_l_nc, &
-                                    i_omp_start, tn )
-                   ELSE
-                      CALL advec_s_pw( i, j, nc )
-                   ENDIF
-                ELSE
-                   CALL advec_s_up( i, j, nc )
-                ENDIF
-                CALL diffusion_s( i, j, nc,                                    &
-                                  surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,    &
-                                  surf_def_h(2)%ncsws,                         &
-                                  surf_lsm_h%ncsws,    surf_usm_h%ncsws,       &
-                                  surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,    &
-                                  surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,    &
-                                  surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,    &
-                                  surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,    &
-                                  surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,    &
-                                  surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
-
-!
-!--             Prognostic equation for cloud drop concentration
-                DO  k = nzb+1, nzt
-                   nc_p(k,j,i) = nc(k,j,i) + ( dt_3d *                         &
-                                                      ( tsc(2) * tend(k,j,i) + &
-                                                        tsc(3) * tnc_m(k,j,i) )&
-                                                      - tsc(5) * rdf_sc(k)     &
-                                                               * nc(k,j,i)     &
-                                             )                                 &
-                                       * MERGE( 1.0_wp, 0.0_wp,                &
-                                                BTEST( wall_flags_0(k,j,i), 0 )&
-                                              )
-                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
-                ENDDO
-!
-!--             Calculate tendencies for the next Runge-Kutta step
-                IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                   IF ( intermediate_timestep_count == 1 )  THEN
-                      DO  k = nzb+1, nzt
-                         tnc_m(k,j,i) = tend(k,j,i)
-                      ENDDO
-                   ELSEIF ( intermediate_timestep_count < &
-                            intermediate_timestep_count_max )  THEN
-                      DO  k = nzb+1, nzt
-                         tnc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
-                                           5.3125_wp * tnc_m(k,j,i)
-                      ENDDO
-                   ENDIF
-                ENDIF
-
-             ENDIF
-!
-!--          If required, calculate prognostic equations for rain water content
-!--          and rain drop concentration
-             IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
-!
-!--             Calculate prognostic equation for rain water content
-                tend(:,j,i) = 0.0_wp
-                IF ( timestep_scheme(1:5) == 'runge' ) &
-                THEN
-                   IF ( ws_scheme_sca )  THEN
-                      CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr,       &
-                                       diss_s_qr, flux_l_qr, diss_l_qr, &
-                                       i_omp_start, tn )
-                   ELSE
-                      CALL advec_s_pw( i, j, qr )
-                   ENDIF
-                ELSE
-                   CALL advec_s_up( i, j, qr )
-                ENDIF
-                CALL diffusion_s( i, j, qr,                                   &
-                                  surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,   &
-                                  surf_def_h(2)%qrsws,                        &
-                                  surf_lsm_h%qrsws,    surf_usm_h%qrsws,      &  
-                                  surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,   &
-                                  surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,   &
-                                  surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,   &
-                                  surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,   &
-                                  surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,   &
-                                  surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
-
-!
-!--             Prognostic equation for rain water content
-                DO  k = nzb+1, nzt
-                   qr_p(k,j,i) = qr(k,j,i) + ( dt_3d *                         &
-                                                      ( tsc(2) * tend(k,j,i) + &
-                                                        tsc(3) * tqr_m(k,j,i) )&
-                                                      - tsc(5) * rdf_sc(k)     &
-                                                               * qr(k,j,i)     &
-                                             )                                 &
-                                       * MERGE( 1.0_wp, 0.0_wp,                &
-                                                BTEST( wall_flags_0(k,j,i), 0 )&
-                                              ) 
-                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
-                ENDDO
-!
-!--             Calculate tendencies for the next Runge-Kutta step
-                IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                   IF ( intermediate_timestep_count == 1 )  THEN
-                      DO  k = nzb+1, nzt
-                         tqr_m(k,j,i) = tend(k,j,i)
-                      ENDDO
-                   ELSEIF ( intermediate_timestep_count < &
-                            intermediate_timestep_count_max )  THEN
-                      DO  k = nzb+1, nzt
-                         tqr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
-                                           5.3125_wp * tqr_m(k,j,i)
-                      ENDDO
-                   ENDIF
-                ENDIF
-
-!
-!--             Calculate prognostic equation for rain drop concentration.
-                tend(:,j,i) = 0.0_wp
-                IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                   IF ( ws_scheme_sca )  THEN
-                      CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr,    &
-                                    diss_s_nr, flux_l_nr, diss_l_nr, &
-                                    i_omp_start, tn )
-                   ELSE
-                      CALL advec_s_pw( i, j, nr )
-                   ENDIF
-                ELSE
-                   CALL advec_s_up( i, j, nr )
-                ENDIF
-                CALL diffusion_s( i, j, nr,                                    &
-                                  surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,    &
-                                  surf_def_h(2)%nrsws,                         &
-                                  surf_lsm_h%nrsws,    surf_usm_h%nrsws,       &
-                                  surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,    &
-                                  surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,    &
-                                  surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,    &
-                                  surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,    &
-                                  surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,    &
-                                  surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
-
-!
-!--             Prognostic equation for rain drop concentration
-                DO  k = nzb+1, nzt
-                   nr_p(k,j,i) = nr(k,j,i) + ( dt_3d *                         &
-                                                      ( tsc(2) * tend(k,j,i) + &
-                                                        tsc(3) * tnr_m(k,j,i) )&
-                                                      - tsc(5) * rdf_sc(k)     &
-                                                               * nr(k,j,i)     &
-                                             )                                 &
-                                       * MERGE( 1.0_wp, 0.0_wp,                &
-                                                BTEST( wall_flags_0(k,j,i), 0 )&
-                                              )
-                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
-                ENDDO
-!
-!--             Calculate tendencies for the next Runge-Kutta step
-                IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                   IF ( intermediate_timestep_count == 1 )  THEN
-                      DO  k = nzb+1, nzt
-                         tnr_m(k,j,i) = tend(k,j,i)
-                      ENDDO
-                   ELSEIF ( intermediate_timestep_count < &
-                            intermediate_timestep_count_max )  THEN
-                      DO  k = nzb+1, nzt
-                         tnr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
-                                           5.3125_wp * tnr_m(k,j,i)
-                      ENDDO
-                   ENDIF
-                ENDIF
-
-             ENDIF
-
           ENDIF
 
@@ -1543,5 +1311,5 @@
        )  THEN
        CALL cpu_log( log_point(51), 'microphysics', 'start' )
-       CALL bcm_actions
+       CALL bcm_actions_micro
        CALL cpu_log( log_point(51), 'microphysics', 'stop' )
     ENDIF
@@ -2110,363 +1878,4 @@
        CALL cpu_log( log_point(29), 'q-equation', 'stop' )
 
-!
-!--    If required, calculate prognostic equations for cloud water content
-!--    and cloud drop concentration
-       IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
-
-          CALL cpu_log( log_point(67), 'qc-equation', 'start' )
-
-!
-!--       Calculate prognostic equation for cloud 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( qc, 'qc' )
-
-          ENDIF
-
-!
-!--       qc-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( qc, 'qc' )
-                ELSE
-                   CALL advec_s_pw( qc )
-                ENDIF
-             ELSE
-                CALL advec_s_up( qc )
-             ENDIF
-          ENDIF
-
-          CALL diffusion_s( qc,                                                &
-                            surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,          &
-                            surf_def_h(2)%qcsws,                               &
-                            surf_lsm_h%qcsws,    surf_usm_h%qcsws,             &
-                            surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,          &
-                            surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,          &
-                            surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,          &
-                            surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,          &
-                            surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,          &
-                            surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
-
-!
-!--       Prognostic equation for cloud water content
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb+1, nzt
-                   qc_p(k,j,i) = qc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
-                                                      tsc(3) * tqc_m(k,j,i) )  &
-                                                    - tsc(5) * rdf_sc(k) *     &
-                                                               qc(k,j,i)       &
-                                             )                                 &
-                                    * MERGE( 1.0_wp, 0.0_wp,                   &
-                                             BTEST( wall_flags_0(k,j,i), 0 )   &
-                                          )
-                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_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+1, nzt
-                         tqc_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+1, nzt
-                         tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
-                                         + 5.3125_wp * tqc_m(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-          ENDIF
-
-          CALL cpu_log( log_point(67), 'qc-equation', 'stop' )
-
-          CALL cpu_log( log_point(68), 'nc-equation', 'start' )
-!
-!--       Calculate prognostic equation for cloud 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( nc, 'nc' )
-
-          ENDIF
-
-!
-!--       nc-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( nc, 'nc' )
-                ELSE
-                   CALL advec_s_pw( nc )
-                ENDIF
-             ELSE
-                CALL advec_s_up( nc )
-             ENDIF
-          ENDIF
-
-          CALL diffusion_s( nc,                                                &
-                            surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,          &
-                            surf_def_h(2)%ncsws,                               &
-                            surf_lsm_h%ncsws,    surf_usm_h%ncsws,             & 
-                            surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,          &
-                            surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,          &
-                            surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,          &
-                            surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,          &
-                            surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,          &
-                            surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
-
-!
-!--       Prognostic equation for cloud drop concentration
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb+1, nzt
-                   nc_p(k,j,i) = nc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
-                                                      tsc(3) * tnc_m(k,j,i) )  &
-                                                    - tsc(5) * rdf_sc(k) *     &
-                                                               nc(k,j,i)       &
-                                             )                                 &
-                                   * MERGE( 1.0_wp, 0.0_wp,                    &
-                                             BTEST( wall_flags_0(k,j,i), 0 )   &
-                                          )
-                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_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+1, nzt
-                         tnc_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+1, nzt
-                         tnc_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
-                                         + 5.3125_wp * tnc_m(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-          ENDIF
-
-          CALL cpu_log( log_point(68), 'nc-equation', 'stop' )
-
-       ENDIF
-!
-!--    If required, calculate prognostic equations for rain water content
-!--    and rain drop concentration
-       IF ( bulk_cloud_model  .AND.  microphysics_seifert )  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,                                                &
-                            surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,          &
-                            surf_def_h(2)%qrsws,                               &
-                            surf_lsm_h%qrsws,    surf_usm_h%qrsws,             &
-                            surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,          &
-                            surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,          &
-                            surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,          &
-                            surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,          &
-                            surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,          &
-                            surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
-
-!
-!--       Prognostic equation for rain water content
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb+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)       &
-                                             )                                 &
-                                    * MERGE( 1.0_wp, 0.0_wp,                   &
-                                             BTEST( wall_flags_0(k,j,i), 0 )   &
-                                          )
-                   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+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+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,                                                &
-                            surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,          &
-                            surf_def_h(2)%nrsws,                               &
-                            surf_lsm_h%nrsws,    surf_usm_h%nrsws,             & 
-                            surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,          &
-                            surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,          &
-                            surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,          &
-                            surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,          &
-                            surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,          &
-                            surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
-
-!
-!--       Prognostic equation for rain drop concentration
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb+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)       &
-                                             )                                 &
-                                   * MERGE( 1.0_wp, 0.0_wp,                    &
-                                             BTEST( wall_flags_0(k,j,i), 0 )   &
-                                          )
-                   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+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+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
 !
