Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 3864)
+++ 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
 !
