Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 1354)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 1361)
@@ -20,5 +20,10 @@
 ! Current revisions:
 ! ------------------
-! 
+! Two-moment microphysics moved to the start of prognostic equations. This makes
+! the 3d arrays for tend_q, tend_qr, tend_pt and tend_pt redundant.
+! Additionally, it is allowed to call the microphysics just once during the time
+! step (not at each sub-time step).
+!
+! Two-moment cloud physics added for vector and accelerator optimization.
 ! 
 ! Former revisions:
@@ -126,11 +131,11 @@
                nrswst, pt, ptdf_x, ptdf_y, pt_init, pt_p, prho, q, q_init,     &
                q_p, qsws, qswst, qr, qr_p, qrsws, qrswst, rdf, rdf_sc, rho,    &
-               sa, sa_init, sa_p, saswsb, saswst, shf, tend, tend_nr,          &
-               tend_pt, tend_q, tend_qr, te_m, tnr_m, tpt_m, tq_m, tqr_m,      &
-               tsa_m, tswst, tu_m, tv_m, tw_m, u, ug, u_p, v, vg, vpt, v_p,    &
-               w, w_p
+               sa, sa_init, sa_p, saswsb, saswst, shf, tend, te_m, tnr_m,      &
+               tpt_m, tq_m, tqr_m, tsa_m, tswst, tu_m, tv_m, tw_m, u, ug, u_p, &
+               v, vg, vpt, v_p, w, w_p
         
     USE control_parameters,                                                    &
-        ONLY:  cloud_physics, constant_diffusion, cthf, dp_external,           &
+        ONLY:  call_microphysics_at_all_substeps, cloud_physics,               &
+               constant_diffusion, cthf, dp_external,                          &
                dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, humidity,       &
                icloud_scheme, inflow_l, intermediate_timestep_count,           &
@@ -301,4 +306,12 @@
  
        DO  j = nys, nyn
+!
+!--       If required, calculate cloud microphysical impacts (two-moment scheme)
+          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
+               ( intermediate_timestep_count == 1  .OR.                        &
+                 call_microphysics_at_all_substeps )                           &
+             )  THEN
+             CALL microphysics_control( i, j )
+          ENDIF
 !
 !--       Tendency terms for u-velocity component
@@ -481,8 +494,5 @@
              ENDIF
           ENDIF
-!
-!--       If required, calculate tendencies for total water content, liquid water 
-!--       potential temperature, rain water content and rain drop concentration
-          IF ( cloud_physics  .AND.  icloud_scheme == 0 )  CALL microphysics_control( i, j )
+
 !
 !--       If required, compute prognostic equation for potential temperature
@@ -511,11 +521,8 @@
 
 !
-!--          Using microphysical tendencies (latent heat)
-             IF ( cloud_physics )  THEN
-                IF ( icloud_scheme == 0 )  THEN
-                   tend(:,j,i) = tend(:,j,i) + tend_pt(:,j,i)
-                ELSEIF ( icloud_scheme == 1  .AND.  precipitation )  THEN
-                   CALL impact_of_latent_heat( i, j ) 
-                ENDIF
+!--          If required compute impact of latent heat due to precipitation
+             IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.              &
+                  precipitation )  THEN
+                CALL impact_of_latent_heat( i, j ) 
              ENDIF
 
@@ -641,11 +648,9 @@
      
 !
-!--          Using microphysical tendencies
-             IF ( cloud_physics )  THEN
-                IF ( icloud_scheme == 0 )  THEN
-                   tend(:,j,i) = tend(:,j,i) + tend_q(:,j,i)
-                ELSEIF ( icloud_scheme == 1  .AND.  precipitation )  THEN
-                   CALL calc_precipitation( i, j )
-                ENDIF
+!--          If required compute decrease of total water content due to
+!--          precipitation
+             IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.              &
+                  precipitation )  THEN
+                CALL calc_precipitation( i, j )             
              ENDIF
 !
@@ -712,8 +717,4 @@
                 ENDIF
                 CALL diffusion_s( i, j, qr, qrsws, qrswst, wall_qrflux )
-!
-!--             Using microphysical tendencies (autoconversion, accretion, 
-!--             evaporation; if required: sedimentation)
-                tend(:,j,i) = tend(:,j,i) + tend_qr(:,j,i)
 
                 CALL user_actions( i, j, 'qr-tendency' )
@@ -757,9 +758,4 @@
                 ENDIF
                 CALL diffusion_s( i, j, nr, nrsws, nrswst, wall_nrflux )
-!
-!--             Using microphysical tendencies (autoconversion, accretion, 
-!--             selfcollection, breakup, evaporation; 
-!--             if required: sedimentation)
-                tend(:,j,i) = tend(:,j,i) + tend_nr(:,j,i)
 
                 CALL user_actions( i, j, 'nr-tendency' )
@@ -883,4 +879,15 @@
 
 !
+!-- If required, calculate cloud microphysical impacts (two-moment scheme)
+    IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                       &
+         ( intermediate_timestep_count == 1  .OR.                              &
+           call_microphysics_at_all_substeps )                                 &
+       )  THEN
+       CALL cpu_log( log_point(51), 'microphysics', 'start' )
+       CALL microphysics_control
+       CALL cpu_log( log_point(51), 'microphysics', 'stop' )
+    ENDIF
+
+!
 !-- u-velocity component
     CALL cpu_log( log_point(5), 'u-equation', 'start' )
@@ -1156,5 +1163,5 @@
 !
 !--    If required compute impact of latent heat due to precipitation
-       IF ( precipitation )  THEN
+       IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.  precipitation )  THEN
           CALL impact_of_latent_heat
        ENDIF
@@ -1348,5 +1355,5 @@
 !--    If required compute decrease of total water content due to
 !--    precipitation
-       IF ( precipitation )  THEN
+       IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.  precipitation )  THEN
           CALL calc_precipitation
        ENDIF
@@ -1406,4 +1413,160 @@
 
        CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
+
+!
+!--    If required, calculate prognostic equations for rain water content 
+!--    and rain drop concentration
+       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
+
+          CALL cpu_log( log_point(52), 'qr-equation', 'start' )
+
+!
+!--       Calculate prognostic equation for rain water content
+          sbt = tsc(2)
+          IF ( scalar_advec == 'bc-scheme' )  THEN
+
+             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--             Bott-Chlond scheme always uses Euler time step. Thus:
+                sbt = 1.0_wp
+             ENDIF
+             tend = 0.0_wp
+             CALL advec_s_bc( qr, 'qr' )
+
+          ENDIF
+
+!
+!--       qr-tendency terms with no communication
+          IF ( scalar_advec /= 'bc-scheme' )  THEN
+             tend = 0.0_wp
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( ws_scheme_sca )  THEN
+                   CALL advec_s_ws( qr, 'qr' )
+                ELSE
+                   CALL advec_s_pw( qr )
+                ENDIF
+             ELSE
+                CALL advec_s_up( qr )
+             ENDIF
+          ENDIF
+
+          CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux )
+
+          CALL user_actions( 'qr-tendency' )
+
+!
+!--       Prognostic equation for rain water content
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +     &
+                                                     tsc(3) * tqr_m(k,j,i) )   &
+                                           - tsc(5) * rdf_sc(k) * qr(k,j,i)
+                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
+                ENDDO
+             ENDDO
+          ENDDO
+
+!
+!--       Calculate tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb_s_inner(j,i)+1, nzt
+                         tqr_m(k,j,i) = tend(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ELSEIF ( intermediate_timestep_count < &
+                      intermediate_timestep_count_max )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb_s_inner(j,i)+1, nzt
+                         tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * &
+                                                                   tqr_m(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+          ENDIF
+
+          CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
+          CALL cpu_log( log_point(53), 'nr-equation', 'start' )
+
+!
+!--       Calculate prognostic equation for rain drop concentration
+          sbt = tsc(2)
+          IF ( scalar_advec == 'bc-scheme' )  THEN
+
+             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--             Bott-Chlond scheme always uses Euler time step. Thus:
+                sbt = 1.0_wp
+             ENDIF
+             tend = 0.0_wp
+             CALL advec_s_bc( nr, 'nr' )
+
+          ENDIF
+
+!
+!--       nr-tendency terms with no communication
+          IF ( scalar_advec /= 'bc-scheme' )  THEN
+             tend = 0.0_wp
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( ws_scheme_sca )  THEN
+                   CALL advec_s_ws( nr, 'nr' )
+                ELSE
+                   CALL advec_s_pw( nr )
+                ENDIF
+             ELSE
+                CALL advec_s_up( nr )
+             ENDIF
+          ENDIF
+
+          CALL diffusion_s( nr, nrsws, nrswst, wall_nrflux )
+
+          CALL user_actions( 'nr-tendency' )
+
+!
+!--       Prognostic equation for rain drop concentration
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +     &
+                                                     tsc(3) * tnr_m(k,j,i) )   &
+                                           - tsc(5) * rdf_sc(k) * nr(k,j,i)
+                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
+                ENDDO
+             ENDDO
+          ENDDO
+
+!
+!--       Calculate tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb_s_inner(j,i)+1, nzt
+                         tnr_m(k,j,i) = tend(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ELSEIF ( intermediate_timestep_count < &
+                      intermediate_timestep_count_max )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb_s_inner(j,i)+1, nzt
+                         tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * &
+                                                                   tnr_m(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+          ENDIF
+
+          CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
+
+       ENDIF
 
     ENDIF
@@ -1510,5 +1673,4 @@
     ENDIF
 
-
  END SUBROUTINE prognostic_equations_vector
 
@@ -1539,4 +1701,15 @@
           runge_step = 2
        ENDIF
+    ENDIF
+
+!
+!-- If required, calculate cloud microphysical impacts (two-moment scheme)
+    IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                       &
+         ( intermediate_timestep_count == 1  .OR.                              &
+           call_microphysics_at_all_substeps )                                 &
+       )  THEN
+       CALL cpu_log( log_point(51), 'microphysics', 'start' )
+       CALL microphysics_control
+       CALL cpu_log( log_point(51), 'microphysics', 'stop' )
     ENDIF
 
@@ -1791,5 +1964,5 @@
 !
 !--    If required compute impact of latent heat due to precipitation
-       IF ( precipitation )  THEN
+       IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.  precipitation )  THEN
           CALL impact_of_latent_heat
        ENDIF
@@ -1957,5 +2130,5 @@
 !--    If required compute decrease of total water content due to
 !--    precipitation
-       IF ( precipitation )  THEN
+       IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.  precipitation )  THEN
           CALL calc_precipitation
        ENDIF
@@ -1999,4 +2172,127 @@
 
        CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
+
+!
+!--    If required, calculate prognostic equations for rain water content 
+!--    and rain drop concentration
+       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
+
+          CALL cpu_log( log_point(52), 'qr-equation', 'start' )
+!
+!--       qr-tendency terms with communication
+          sbt = tsc(2)
+          IF ( scalar_advec == 'bc-scheme' )  THEN
+
+             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--             Bott-Chlond scheme always uses Euler time step. Thus:
+                sbt = 1.0_wp
+             ENDIF
+             tend = 0.0_wp
+             CALL advec_s_bc( qr, 'qr' )
+
+          ENDIF
+
+!
+!--       qr-tendency terms with no communication
+          IF ( scalar_advec /= 'bc-scheme' )  THEN
+             tend = 0.0_wp
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( ws_scheme_sca )  THEN
+                   CALL advec_s_ws( qr, 'qr' )
+                ELSE
+                   CALL advec_s_pw( qr )
+                ENDIF
+             ELSE
+                CALL advec_s_up( qr )
+             ENDIF
+          ENDIF
+
+          CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux )
+
+          CALL user_actions( 'qr-tendency' )
+
+!
+!--       Prognostic equation for rain water content
+          DO  i = i_left, i_right
+             DO  j = j_south, j_north
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +     &
+                                                       tsc(3) * tqr_m(k,j,i) ) &
+                                           - tsc(5) * rdf_sc(k) * qr(k,j,i)
+                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
+!
+!--                Tendencies for the next Runge-Kutta step
+                   IF ( runge_step == 1 )  THEN
+                      tqr_m(k,j,i) = tend(k,j,i)
+                   ELSEIF ( runge_step == 2 )  THEN
+                      tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp *    &
+                                                                tqr_m(k,j,i)
+                   ENDIF
+                ENDDO
+             ENDDO
+          ENDDO
+
+          CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
+          CALL cpu_log( log_point(53), 'nr-equation', 'start' )
+
+!
+!--       nr-tendency terms with communication
+          sbt = tsc(2)
+          IF ( scalar_advec == 'bc-scheme' )  THEN
+
+             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--             Bott-Chlond scheme always uses Euler time step. Thus:
+                sbt = 1.0_wp
+             ENDIF
+             tend = 0.0_wp
+             CALL advec_s_bc( nr, 'nr' )
+
+          ENDIF
+
+!
+!--       nr-tendency terms with no communication
+          IF ( scalar_advec /= 'bc-scheme' )  THEN
+             tend = 0.0_wp
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( ws_scheme_sca )  THEN
+                   CALL advec_s_ws( nr, 'nr' )
+                ELSE
+                   CALL advec_s_pw( nr )
+                ENDIF
+             ELSE
+                CALL advec_s_up( nr )
+             ENDIF
+          ENDIF
+
+          CALL diffusion_s( nr, nrsws, nrswst, wall_nrflux )
+
+          CALL user_actions( 'nr-tendency' )
+
+!
+!--       Prognostic equation for rain drop concentration
+          DO  i = i_left, i_right
+             DO  j = j_south, j_north
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +     &
+                                                       tsc(3) * tnr_m(k,j,i) ) &
+                                           - tsc(5) * rdf_sc(k) * nr(k,j,i)
+                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
+!
+!--                Tendencies for the next Runge-Kutta step
+                   IF ( runge_step == 1 )  THEN
+                      tnr_m(k,j,i) = tend(k,j,i)
+                   ELSEIF ( runge_step == 2 )  THEN
+                      tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp *    &
+                                                                tnr_m(k,j,i)
+                   ENDIF
+                ENDDO
+             ENDDO
+          ENDDO
+
+          CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
+
+       ENDIF
 
     ENDIF
@@ -2095,5 +2391,4 @@
     ENDIF
 
-
  END SUBROUTINE prognostic_equations_acc
 
