Index: palm/trunk/SOURCE/microphysics_mod.f90
===================================================================
--- palm/trunk/SOURCE/microphysics_mod.f90	(revision 2153)
+++ palm/trunk/SOURCE/microphysics_mod.f90	(revision 2155)
@@ -20,6 +20,6 @@
 ! Current revisions:
 ! ------------------
-! 
-! 
+! Bugfix in the calculation of microphysical quantities on ghost points.
+!
 ! Former revisions:
 ! -----------------
@@ -28,8 +28,8 @@
 ! 2031 2016-10-21 15:11:58Z knoop
 ! renamed variable rho to rho_ocean
-! 
+!
 ! 2000 2016-08-20 18:09:15Z knoop
 ! Forced header and separation lines into 80 columns
-! 
+!
 ! 1850 2016-04-08 13:29:27Z maronga
 ! Module renamed
@@ -52,7 +52,7 @@
 ! Added new routine calc_precipitation_amount. The routine now allows to account
 ! for precipitation due to sedimenation of cloud (fog) droplets
-! 
+!
 ! 1682 2015-10-07 23:56:08Z knoop
-! Code annotations made doxygen readable 
+! Code annotations made doxygen readable
 !
 ! 1646 2015-09-02 16:00:10Z hoffmann
@@ -63,12 +63,12 @@
 ! Vectorized version of adjust_cloud added.
 ! Little reformatting of the code.
-! 
+!
 ! 1353 2014-04-08 15:21:23Z heinze
 ! REAL constants provided with KIND-attribute
-! 
+!
 ! 1346 2014-03-27 13:18:20Z heinze
-! Bugfix: REAL constants provided with KIND-attribute especially in call of 
+! Bugfix: REAL constants provided with KIND-attribute especially in call of
 ! intrinsic function like MAX, MIN, SIGN
-! 
+!
 ! 1334 2014-03-25 12:21:40Z heinze
 ! Bugfix: REAL constants provided with KIND-attribute
@@ -106,5 +106,5 @@
 ! 1053 2012-11-13 17:11:03Z hoffmann
 ! initial revision
-! 
+!
 ! Description:
 ! ------------
@@ -140,5 +140,5 @@
     REAL(wp) ::  c_term = 600.0_wp         !< coef. for terminal velocity (m-1)
     REAL(wp) ::  diff_coeff_l = 0.23E-4_wp  !< diffusivity of water vapor (m2 s-1)
-    REAL(wp) ::  eps_sb = 1.0E-20_wp       !< threshold in two-moments scheme
+    REAL(wp) ::  eps_sb = 1.0E-10_wp       !< threshold in two-moments scheme
     REAL(wp) ::  k_cc = 9.44E09_wp         !< const. cloud-cloud kernel (m3 kg-2 s-1)
     REAL(wp) ::  k_cr0 = 4.33_wp           !< const. cloud-rain kernel (m3 kg-1 s-1)
@@ -221,5 +221,5 @@
        MODULE PROCEDURE sedimentation_cloud_ij
     END INTERFACE sedimentation_cloud
- 
+
     INTERFACE sedimentation_rain
        MODULE PROCEDURE sedimentation_rain
@@ -337,5 +337,5 @@
              pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
              t_d_pt(k) = 1.0_wp / pt_d_t(k)
-             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
+             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
           ENDDO
 
@@ -346,5 +346,5 @@
 
 !
-!--    Compute length of time step 
+!--    Compute length of time step
        IF ( call_microphysics_at_all_substeps )  THEN
           dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
@@ -383,6 +383,6 @@
 ! Description:
 ! ------------
-!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and 
-!> evaporation of rain drops due to too small or too big weights 
+!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
+!> evaporation of rain drops due to too small or too big weights
 !> of rain drops (Stevens and Seifert, 2008).
 !------------------------------------------------------------------------------!
@@ -399,5 +399,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt
 
        USE kinds
@@ -411,6 +411,6 @@
        CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
 
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              DO  k = nzb_s_inner(j,i)+1, nzt
                 IF ( qr(k,j,i) <= eps_sb )  THEN
@@ -456,5 +456,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt
 
        USE kinds
@@ -466,5 +466,5 @@
        INTEGER(iwp) ::  k                 !<
 
-       REAL(wp)     ::  alpha_cc          !<                   
+       REAL(wp)     ::  alpha_cc          !<
        REAL(wp)     ::  autocon           !<
        REAL(wp)     ::  dissipation       !<
@@ -482,6 +482,6 @@
        CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
 
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              DO  k = nzb_s_inner(j,i)+1, nzt
 
@@ -494,5 +494,5 @@
                    tau_cloud = 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) )
 !
-!--                Universal function for autoconversion process 
+!--                Universal function for autoconversion process
 !--                (Seifert and Beheng, 2006):
                    phi_au = 600.0_wp * tau_cloud**0.68_wp *                    &
@@ -506,5 +506,5 @@
                    xc = hyrho(k) * qc(k,j,i) / nc_const
 !
-!--                Parameterized turbulence effects on autoconversion (Seifert, 
+!--                Parameterized turbulence effects on autoconversion (Seifert,
 !--                Nuijens and Stevens, 2010)
                    IF ( collision_turbulence )  THEN
@@ -517,9 +517,9 @@
                       sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
 !
-!--                   Mixing length (neglecting distance to ground and 
+!--                   Mixing length (neglecting distance to ground and
 !--                   stratification)
                       l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
 !
-!--                   Limit dissipation rate according to Seifert, Nuijens and 
+!--                   Limit dissipation rate according to Seifert, Nuijens and
 !--                   Stevens (2010)
                       dissipation = MIN( 0.06_wp, diss(k,j,i) )
@@ -531,10 +531,10 @@
                                   dissipation**( 1.0_wp / 6.0_wp )
 !
-!--                   The factor of 1.0E4 is needed to convert the dissipation 
+!--                   The factor of 1.0E4 is needed to convert the dissipation
 !--                   rate from m2 s-3 to cm2 s-3.
                       k_au = k_au * ( 1.0_wp +                                 &
                                       dissipation * 1.0E4_wp *                 &
                                       ( re_lambda * 1.0E-3_wp )**0.25_wp *     &
-                                      ( alpha_cc * EXP( -1.0_wp * ( ( rc -     &      
+                                      ( alpha_cc * EXP( -1.0_wp * ( ( rc -     &
                                                                       r_cc ) / &
                                                         sigma_cc )**2          &
@@ -552,5 +552,5 @@
 
                    qr(k,j,i) = qr(k,j,i) + autocon * dt_micro
-                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro 
+                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro
                    nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro
 
@@ -580,5 +580,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxl, nxr, nyn, nys, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb_s_inner, nzt
 
        USE kinds
@@ -593,6 +593,6 @@
        REAL(wp)    ::  dqdt_precip !<
 
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              DO  k = nzb_s_inner(j,i)+1, nzt
 
@@ -640,5 +640,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt
 
        USE kinds
@@ -657,6 +657,6 @@
        CALL cpu_log( log_point_s(56), 'accretion', 'start' )
 
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              DO  k = nzb_s_inner(j,i)+1, nzt
 
@@ -664,12 +664,12 @@
 !
 !--                Intern time scale of coagulation (Seifert and Beheng, 2006):
-                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ) 
-!
-!--                Universal function for accretion process (Seifert and 
+                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
+!
+!--                Universal function for accretion process (Seifert and
 !--                Beheng, 2001):
                    phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
 !
-!--                Parameterized turbulence effects on autoconversion (Seifert, 
-!--                Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to 
+!--                Parameterized turbulence effects on autoconversion (Seifert,
+!--                Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
 !--                convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
                    IF ( collision_turbulence )  THEN
@@ -679,5 +679,5 @@
                                      )
                    ELSE
-                      k_cr = k_cr0                       
+                      k_cr = k_cr0
                    ENDIF
 !
@@ -687,5 +687,5 @@
                    accr = MIN( accr, qc(k,j,i) / dt_micro )
 
-                   qr(k,j,i) = qr(k,j,i) + accr * dt_micro 
+                   qr(k,j,i) = qr(k,j,i) + accr * dt_micro
                    qc(k,j,i) = qc(k,j,i) - accr * dt_micro
 
@@ -721,8 +721,8 @@
 
        USE indices,                                                            &
-           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt
 
        USE kinds
-   
+
        IMPLICIT NONE
 
@@ -738,6 +738,6 @@
        CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
 
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              DO  k = nzb_s_inner(j,i)+1, nzt
                 IF ( qr(k,j,i) > eps_sb )  THEN
@@ -762,5 +762,5 @@
                    nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro
 
-                ENDIF          
+                ENDIF
              ENDDO
           ENDDO
@@ -775,5 +775,5 @@
 ! Description:
 ! ------------
-!> Evaporation of precipitable water. Condensation is neglected for 
+!> Evaporation of precipitable water. Condensation is neglected for
 !> precipitable water.
 !------------------------------------------------------------------------------!
@@ -793,5 +793,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt
 
        USE kinds
@@ -823,6 +823,6 @@
        CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
 
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              DO  k = nzb_s_inner(j,i)+1, nzt
                 IF ( qr(k,j,i) > eps_sb )  THEN
@@ -850,5 +850,5 @@
 !--                   Actual temperature:
                       temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
-             
+
                       g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *     &
                                           l_v / ( thermal_conductivity_l * temp ) &
@@ -862,9 +862,9 @@
                       dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
 !
-!--                   Compute ventilation factor and intercept parameter 
+!--                   Compute ventilation factor and intercept parameter
 !--                   (Seifert and Beheng, 2006; Seifert, 2008):
                       IF ( ventilation_effect )  THEN
 !
-!--                      Shape parameter of gamma distribution (Milbrandt and Yau, 
+!--                      Shape parameter of gamma distribution (Milbrandt and Yau,
 !--                      2005; Stevens and Seifert, 2008):
                          mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *          &
@@ -877,5 +877,5 @@
 
                          mu_r_2   = mu_r + 2.0_wp
-                         mu_r_5d2 = mu_r + 2.5_wp 
+                         mu_r_5d2 = mu_r + 2.5_wp
 
                          f_vent = a_vent * gamm( mu_r_2 ) *                     &
@@ -893,5 +893,5 @@
                                     ( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
                                     )**mu_r_5d2 -                               &
-                                    0.0390625_wp * ( b_term / a_term )**4 *     & 
+                                    0.0390625_wp * ( b_term / a_term )**4 *     &
                                     ( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
                                     )**mu_r_5d2                                 &
@@ -899,5 +899,5 @@
 
                          nr_0   = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) /    &
-                                  gamm( mu_r + 1.0_wp ) 
+                                  gamm( mu_r + 1.0_wp )
                       ELSE
                          f_vent = 1.0_wp
@@ -905,5 +905,5 @@
                       ENDIF
    !
-   !--                Evaporation rate of rain water content (Seifert and 
+   !--                Evaporation rate of rain water content (Seifert and
    !--                Beheng, 2006):
                       evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat /   &
@@ -917,5 +917,5 @@
 
                    ENDIF
-                ENDIF          
+                ENDIF
 
              ENDDO
@@ -948,5 +948,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzb_s_inner, nzt
 
        USE kinds
@@ -954,5 +954,5 @@
        USE statistics,                                                         &
            ONLY:  weight_substep
-    
+
 
        IMPLICIT NONE
@@ -968,6 +968,6 @@
        sed_qc(nzt+1) = 0.0_wp
 
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              DO  k = nzt, nzb_s_inner(j,i)+1, -1
 
@@ -1029,5 +1029,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzb_s_inner, nzt
 
        USE kinds
@@ -1035,5 +1035,5 @@
        USE statistics,                                                         &
            ONLY:  weight_substep
-       
+
        IMPLICIT NONE
 
@@ -1065,7 +1065,7 @@
 
 !
-!--    Compute velocities 
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
+!--    Compute velocities
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              DO  k = nzb_s_inner(j,i)+1, nzt
                 IF ( qr(k,j,i) > eps_sb )  THEN
@@ -1119,5 +1119,5 @@
                                       2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
                           dt_micro * ddzu(k)
-             ENDDO     
+             ENDDO
 !
 !--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
@@ -1155,5 +1155,5 @@
              DO  k = nzt, nzb_s_inner(j,i)+1, -1
 !
-!--             Sum up all rain drop number densities which contribute to the flux 
+!--             Sum up all rain drop number densities which contribute to the flux
 !--             through k-1/2
                 flux  = 0.0_wp
@@ -1170,5 +1170,5 @@
                 ENDDO
 !
-!--             It is not allowed to sediment more rain drop number density than 
+!--             It is not allowed to sediment more rain drop number density than
 !--             available
                 flux = MIN( flux,                                              &
@@ -1181,5 +1181,5 @@
                                         ddzu(k+1) / hyrho(k) * dt_micro
 !
-!--             Sum up all rain water content which contributes to the flux 
+!--             Sum up all rain water content which contributes to the flux
 !--             through k-1/2
                 flux  = 0.0_wp
@@ -1199,5 +1199,5 @@
                 ENDDO
 !
-!--             It is not allowed to sediment more rain water content than 
+!--             It is not allowed to sediment more rain water content than
 !--             available
                 flux = MIN( flux,                                              &
@@ -1211,5 +1211,5 @@
                                         ddzu(k+1) / hyrho(k) * dt_micro
                 q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
-                                        ddzu(k+1) / hyrho(k) * dt_micro 
+                                        ddzu(k+1) / hyrho(k) * dt_micro
                 pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
                                         ddzu(k+1) / hyrho(k) * l_d_cp *        &
@@ -1330,5 +1330,5 @@
              pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
              t_d_pt(k) = 1.0_wp / pt_d_t(k)
-             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
+             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
           ENDDO
 !
@@ -1338,5 +1338,5 @@
 
 !
-!--    Compute length of time step 
+!--    Compute length of time step
        IF ( call_microphysics_at_all_substeps )  THEN
           dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
@@ -1395,5 +1395,5 @@
 ! Description:
 ! ------------
-!> Adjust number of raindrops to avoid nonlinear effects in 
+!> Adjust number of raindrops to avoid nonlinear effects in
 !> sedimentation and evaporation of rain drops due to too small or
 !> too big weights of rain drops (Stevens and Seifert, 2008).
@@ -1424,5 +1424,5 @@
           ELSE
 !
-!--          Adjust number of raindrops to avoid nonlinear effects in 
+!--          Adjust number of raindrops to avoid nonlinear effects in
 !--          sedimentation and evaporation of rain drops due to too small or
 !--          too big weights of rain drops (Stevens and Seifert, 2008).
@@ -1470,5 +1470,5 @@
        INTEGER(iwp) ::  k                 !<
 
-       REAL(wp)     ::  alpha_cc          !<                   
+       REAL(wp)     ::  alpha_cc          !<
        REAL(wp)     ::  autocon           !<
        REAL(wp)     ::  dissipation       !<
@@ -1494,5 +1494,5 @@
              tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
 !
-!--          Universal function for autoconversion process 
+!--          Universal function for autoconversion process
 !--          (Seifert and Beheng, 2006):
              phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
@@ -1505,5 +1505,5 @@
              xc = hyrho(k) * qc_1d(k) / nc_1d(k)
 !
-!--          Parameterized turbulence effects on autoconversion (Seifert, 
+!--          Parameterized turbulence effects on autoconversion (Seifert,
 !--          Nuijens and Stevens, 2010)
              IF ( collision_turbulence )  THEN
@@ -1519,5 +1519,5 @@
                 l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
 !
-!--             Limit dissipation rate according to Seifert, Nuijens and 
+!--             Limit dissipation rate according to Seifert, Nuijens and
 !--             Stevens (2010)
                 dissipation = MIN( 0.06_wp, diss(k,j,i) )
@@ -1549,5 +1549,5 @@
 
              qr_1d(k) = qr_1d(k) + autocon * dt_micro
-             qc_1d(k) = qc_1d(k) - autocon * dt_micro 
+             qc_1d(k) = qc_1d(k) - autocon * dt_micro
              nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
 
@@ -1642,12 +1642,12 @@
 !
 !--          Intern time scale of coagulation (Seifert and Beheng, 2006):
-             tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 
-!
-!--          Universal function for accretion process 
+             tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
+!
+!--          Universal function for accretion process
 !--          (Seifert and Beheng, 2001):
              phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
 !
-!--          Parameterized turbulence effects on autoconversion (Seifert, 
-!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to 
+!--          Parameterized turbulence effects on autoconversion (Seifert,
+!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
 !--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
              IF ( collision_turbulence )  THEN
@@ -1657,5 +1657,5 @@
                                )
              ELSE
-                k_cr = k_cr0                       
+                k_cr = k_cr0
              ENDIF
 !
@@ -1664,5 +1664,5 @@
              accr = MIN( accr, qc_1d(k) / dt_micro )
 
-             qr_1d(k) = qr_1d(k) + accr * dt_micro 
+             qr_1d(k) = qr_1d(k) + accr * dt_micro
              qc_1d(k) = qc_1d(k) - accr * dt_micro
 
@@ -1691,5 +1691,5 @@
 
        USE kinds
-   
+
        IMPLICIT NONE
 
@@ -1723,5 +1723,5 @@
              nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
 
-          ENDIF          
+          ENDIF
        ENDDO
 
@@ -1732,5 +1732,5 @@
 ! Description:
 ! ------------
-!> Evaporation of precipitable water. Condensation is neglected for 
+!> Evaporation of precipitable water. Condensation is neglected for
 !> precipitable water. Call for grid point i,j
 !------------------------------------------------------------------------------!
@@ -1799,5 +1799,5 @@
 !--             Actual temperature:
                 temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
-       
+
                 g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v /  &
                                     ( thermal_conductivity_l * temp ) +        &
@@ -1811,5 +1811,5 @@
                 dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
 !
-!--             Compute ventilation factor and intercept parameter 
+!--             Compute ventilation factor and intercept parameter
 !--             (Seifert and Beheng, 2006; Seifert, 2008):
                 IF ( ventilation_effect )  THEN
@@ -1825,7 +1825,7 @@
 
                    mu_r_2   = mu_r + 2.0_wp
-                   mu_r_5d2 = mu_r + 2.5_wp 
-
-                   f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) +  & 
+                   mu_r_5d2 = mu_r + 2.5_wp
+
+                   f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) +  &
                             b_vent * schmidt_p_1d3 *                           &
                             SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *  &
@@ -1841,5 +1841,5 @@
                               ( lambda_r / ( 3.0_wp * c_term + lambda_r )      &
                               )**mu_r_5d2 -                                    &
-                              0.0390625_wp * ( b_term / a_term )**4 *          & 
+                              0.0390625_wp * ( b_term / a_term )**4 *          &
                               ( lambda_r / ( 4.0_wp * c_term + lambda_r )      &
                               )**mu_r_5d2                                      &
@@ -1847,5 +1847,5 @@
 
                    nr_0   = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) /           &
-                            gamm( mu_r + 1.0_wp ) 
+                            gamm( mu_r + 1.0_wp )
                 ELSE
                    f_vent = 1.0_wp
@@ -1863,5 +1863,5 @@
 
              ENDIF
-          ENDIF          
+          ENDIF
 
        ENDDO
@@ -1873,5 +1873,5 @@
 ! Description:
 ! ------------
-!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR). 
+!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
 !> Call for grid point i,j
 !------------------------------------------------------------------------------!
@@ -1891,5 +1891,5 @@
 
        USE kinds
-       
+
        USE statistics,                                                         &
            ONLY:  weight_substep
@@ -1919,5 +1919,5 @@
           q_1d(k)  = q_1d(k)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
                                 hyrho(k) * dt_micro
-          qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      & 
+          qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
                                 hyrho(k) * dt_micro
           pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
@@ -1962,5 +1962,5 @@
        USE statistics,                                                         &
            ONLY:  weight_substep
-       
+
        IMPLICIT NONE
 
@@ -1990,5 +1990,5 @@
 
 !
-!--    Compute velocities 
+!--    Compute velocities
        DO  k = nzb_s_inner(j,i)+1, nzt
           IF ( qr_1d(k) > eps_sb )  THEN
@@ -2035,5 +2035,5 @@
           c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) *   &
                     dt_micro * ddzu(k)
-       ENDDO     
+       ENDDO
 !
 !--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
@@ -2071,5 +2071,5 @@
        DO  k = nzt, nzb_s_inner(j,i)+1, -1
 !
-!--       Sum up all rain drop number densities which contribute to the flux 
+!--       Sum up all rain drop number densities which contribute to the flux
 !--       through k-1/2
           flux  = 0.0_wp
@@ -2086,5 +2086,5 @@
           ENDDO
 !
-!--       It is not allowed to sediment more rain drop number density than 
+!--       It is not allowed to sediment more rain drop number density than
 !--       available
           flux = MIN( flux,                                                    &
@@ -2095,5 +2095,5 @@
                                     hyrho(k) * dt_micro
 !
-!--       Sum up all rain water content which contributes to the flux 
+!--       Sum up all rain water content which contributes to the flux
 !--       through k-1/2
           flux  = 0.0_wp
@@ -2122,5 +2122,5 @@
                                 hyrho(k) * dt_micro
           q_1d(k)  = q_1d(k)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
-                                hyrho(k) * dt_micro 
+                                hyrho(k) * dt_micro
           pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
                                 hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
@@ -2185,12 +2185,12 @@
 ! ------------
 !> This function computes the gamma function (Press et al., 1992).
-!> The gamma function is needed for the calculation of the evaporation 
+!> The gamma function is needed for the calculation of the evaporation
 !> of rain drops.
 !------------------------------------------------------------------------------!
-    FUNCTION gamm( xx ) 
+    FUNCTION gamm( xx )
 
        USE kinds
 
-       IMPLICIT NONE 
+       IMPLICIT NONE
 
        INTEGER(iwp) ::  j            !<
@@ -2218,18 +2218,18 @@
        ser = 1.000000000190015_wp
 
-       DO  j = 1, 6 
-          y_gamm = y_gamm + 1.0_wp 
-          ser    = ser + cof( j ) / y_gamm 
+       DO  j = 1, 6
+          y_gamm = y_gamm + 1.0_wp
+          ser    = ser + cof( j ) / y_gamm
        ENDDO
 
-! 
-!--    Until this point the algorithm computes the logarithm of the gamma 
-!--    function. Hence, the exponential function is used.  
-!       gamm = EXP( tmp + LOG( stp * ser / x_gamm ) ) 
-       gamm = EXP( tmp ) * stp * ser / x_gamm 
-
-       RETURN 
- 
-    END FUNCTION gamm 
+!
+!--    Until this point the algorithm computes the logarithm of the gamma
+!--    function. Hence, the exponential function is used.
+!       gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
+       gamm = EXP( tmp ) * stp * ser / x_gamm
+
+       RETURN
+
+    END FUNCTION gamm
 
  END MODULE microphysics_mod
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 2153)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 2155)
@@ -20,6 +20,6 @@
 ! Current revisions:
 ! ------------------
-! 
-! 
+! Bugfix in the calculation of microphysical quantities on ghost points.
+!
 ! Former revisions:
 ! -----------------
@@ -28,24 +28,24 @@
 ! 2118 2017-01-17 16:38:49Z raasch
 ! OpenACC version of subroutine removed
-! 
+!
 ! 2031 2016-10-21 15:11:58Z knoop
 ! renamed variable rho to rho_ocean
-! 
+!
 ! 2011 2016-09-19 17:29:57Z kanani
 ! Flag urban_surface is now defined in module control_parameters.
-! 
+!
 ! 2007 2016-08-24 15:47:17Z kanani
 ! Added pt tendency calculation based on energy balance at urban surfaces
 ! (new urban surface model)
-! 
+!
 ! 2000 2016-08-20 18:09:15Z knoop
 ! Forced header and separation lines into 80 columns
-! 
+!
 ! 1976 2016-07-27 13:28:04Z maronga
 ! Simplied calls to radiation model
-! 
+!
 ! 1960 2016-07-12 16:34:24Z suehring
 ! Separate humidity and passive scalar
-! 
+!
 ! 1914 2016-05-26 14:44:07Z witha
 ! Added calls for wind turbine model
@@ -53,26 +53,26 @@
 ! 1873 2016-04-18 14:50:06Z maronga
 ! Module renamed (removed _mod)
-! 
+!
 ! 1850 2016-04-08 13:29:27Z maronga
 ! Module renamed
-! 
+!
 ! 1826 2016-04-07 12:01:39Z maronga
 ! Renamed canopy model calls.
-! 
+!
 ! 1822 2016-04-07 07:49:42Z hoffmann
 ! Kessler microphysics scheme moved to microphysics.
 !
 ! 1757 2016-02-22 15:49:32Z maronga
-! 
+!
 ! 1691 2015-10-26 16:17:44Z maronga
 ! Added optional model spin-up without radiation / land surface model calls.
 ! Formatting corrections.
-! 
+!
 ! 1682 2015-10-07 23:56:08Z knoop
-! Code annotations made doxygen readable 
-! 
+! Code annotations made doxygen readable
+!
 ! 1585 2015-04-30 07:05:52Z maronga
 ! Added call for temperature tendency calculation due to radiative flux divergence
-! 
+!
 ! 1517 2015-01-07 19:12:25Z hoffmann
 ! advec_s_bc_mod addded, since advec_s_bc is now a module
@@ -80,5 +80,5 @@
 ! 1496 2014-12-02 17:25:50Z maronga
 ! Renamed "radiation" -> "cloud_top_radiation"
-! 
+!
 ! 1484 2014-10-21 10:53:05Z kanani
 ! Changes due to new module structure of the plant canopy model:
@@ -86,28 +86,28 @@
 ! Removed double-listing of use_upstream_for_tke in ONLY-list of module
 ! control_parameters
-! 
+!
 ! 1409 2014-05-23 12:11:32Z suehring
-! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary. 
+! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary.
 ! This ensures that left-hand side fluxes are also calculated for nxl in that
-! case, even though the solution at nxl is overwritten in boundary_conds() 
-! 
+! case, even though the solution at nxl is overwritten in boundary_conds()
+!
 ! 1398 2014-05-07 11:15:00Z heinze
 ! Rayleigh-damping for horizontal velocity components changed: instead of damping
-! against ug and vg, damping against u_init and v_init is used to allow for a 
+! against ug and vg, damping against u_init and v_init is used to allow for a
 ! homogenized treatment in case of nudging
-! 
+!
 ! 1380 2014-04-28 12:40:45Z heinze
-! Change order of calls for scalar prognostic quantities: 
-! ls_advec -> nudging -> subsidence since initial profiles 
-! 
+! Change order of calls for scalar prognostic quantities:
+! ls_advec -> nudging -> subsidence since initial profiles
+!
 ! 1374 2014-04-25 12:55:07Z raasch
 ! missing variables added to ONLY lists
-! 
+!
 ! 1365 2014-04-22 15:03:56Z boeske
-! Calls of ls_advec for large scale advection added, 
+! Calls of ls_advec for large scale advection added,
 ! subroutine subsidence is only called if use_subsidence_tendencies = .F.,
 ! new argument ls_index added to the calls of subsidence
 ! +ls_index
-! 
+!
 ! 1361 2014-04-16 15:17:48Z hoffmann
 ! Two-moment microphysics moved to the start of prognostic equations. This makes
@@ -117,16 +117,16 @@
 !
 ! Two-moment cloud physics added for vector and accelerator optimization.
-! 
+!
 ! 1353 2014-04-08 15:21:23Z heinze
 ! REAL constants provided with KIND-attribute
-! 
+!
 ! 1337 2014-03-25 15:11:48Z heinze
 ! Bugfix: REAL constants provided with KIND-attribute
-! 
+!
 ! 1332 2014-03-25 11:59:43Z suehring
-! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke 
-! 
+! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke
+!
 ! 1330 2014-03-24 17:29:32Z suehring
-! In case of SGS-particle velocity advection of TKE is also allowed with 
+! In case of SGS-particle velocity advection of TKE is also allowed with
 ! dissipative 5th-order scheme.
 !
@@ -161,5 +161,5 @@
 !
 ! 1115 2013-03-26 18:16:16Z hoffmann
-! optimized cloud physics: calculation of microphysical tendencies transfered 
+! optimized cloud physics: calculation of microphysical tendencies transfered
 ! to microphysics.f90; qr and nr are only calculated if precipitation is required
 !
@@ -211,5 +211,5 @@
  MODULE prognostic_equations_mod
 
- 
+
 
     USE arrays_3d,                                                             &
@@ -225,5 +225,5 @@
                te_m, tnr_m, tpt_m, tq_m, tqr_m, ts_m, tsa_m, tswst, tu_m, tv_m,&
                tw_m, u, ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
-        
+
     USE control_parameters,                                                    &
         ONLY:  call_microphysics_at_all_substeps, cloud_physics,               &
@@ -249,6 +249,6 @@
 
     USE indices,                                                               &
-        ONLY:  nxl, nxlu, nxr, nyn, nys, nysv, nzb_s_inner, nzb_u_inner,       &
-               nzb_v_inner, nzb_w_inner, nzt
+        ONLY:  nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,         &
+               nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt
 
     USE advec_ws,                                                              &
@@ -287,5 +287,5 @@
     USE calc_radiation_mod,                                                    &
         ONLY:  calc_radiation
-  
+
     USE coriolis_mod,                                                          &
         ONLY:  coriolis
@@ -370,5 +370,5 @@
 !> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
 !------------------------------------------------------------------------------!
- 
+
  SUBROUTINE prognostic_equations_cache
 
@@ -382,5 +382,5 @@
     INTEGER(iwp) ::  omp_get_thread_num  !<
     INTEGER(iwp) ::  tn = 0              !<
-    
+
     LOGICAL      ::  loop_start          !<
 
@@ -394,7 +394,21 @@
 !$OMP PARALLEL private (i,i_omp_start,j,k,loop_start,tn)
 
-!$  tn = omp_get_thread_num() 
+!$  tn = omp_get_thread_num()
     loop_start = .TRUE.
 !$OMP DO
+
+!
+!-- If required, calculate cloud microphysics
+    IF ( cloud_physics  .AND.  .NOT. microphysics_sat_adjust  .AND.            &
+         ( intermediate_timestep_count == 1  .OR.                              &
+           call_microphysics_at_all_substeps )                                 &
+       )  THEN
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             CALL microphysics_control( i, j )
+           END DO
+       END DO
+    ENDIF
+
     DO  i = nxl, nxr
 
@@ -404,16 +418,8 @@
        IF ( loop_start )  THEN
           loop_start  = .FALSE.
-          i_omp_start = i 
+          i_omp_start = i
        ENDIF
 
        DO  j = nys, nyn
-!
-!--       If required, calculate cloud microphysics
-          IF ( cloud_physics  .AND.  .NOT. microphysics_sat_adjust  .AND.      &
-               ( intermediate_timestep_count == 1  .OR.                        &
-                 call_microphysics_at_all_substeps )                           &
-             )  THEN
-             CALL microphysics_control( i, j )
-          ENDIF
 !
 !--       Tendency terms for u-velocity component
@@ -424,7 +430,7 @@
                 IF ( ws_scheme_mom )  THEN
                    CALL advec_u_ws( i, j, i_omp_start, tn )
-                ELSE 
+                ELSE
                    CALL advec_u_pw( i, j )
-                ENDIF 
+                ENDIF
              ELSE
                 CALL advec_u_up( i, j )
@@ -490,5 +496,5 @@
                 IF ( ws_scheme_mom )  THEN
                     CALL advec_v_ws( i, j, i_omp_start, tn )
-                ELSE 
+                ELSE
                     CALL advec_v_pw( i, j )
                 ENDIF
@@ -501,5 +507,5 @@
 !
 !--          Drag by plant canopy
-             IF ( plant_canopy )  CALL pcm_tendency( i, j, 2 )        
+             IF ( plant_canopy )  CALL pcm_tendency( i, j, 2 )
 
 !
@@ -551,5 +557,5 @@
              IF ( ws_scheme_mom )  THEN
                 CALL advec_w_ws( i, j, i_omp_start, tn )
-             ELSE 
+             ELSE
                 CALL advec_w_pw( i, j )
              END IF
@@ -646,9 +652,9 @@
              IF ( large_scale_forcing )  THEN
                 CALL ls_advec( i, j, simulated_time, 'pt' )
-             ENDIF     
+             ENDIF
 
 !
 !--          Nudging
-             IF ( nudging )  CALL nudge( i, j, simulated_time, 'pt' ) 
+             IF ( nudging )  CALL nudge( i, j, simulated_time, 'pt' )
 
 !
@@ -708,5 +714,5 @@
                     CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa,  &
                                 diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn  )
-                ELSE 
+                ELSE
                     CALL advec_s_pw( i, j, sa )
                 ENDIF
@@ -760,7 +766,7 @@
              THEN
                 IF ( ws_scheme_sca )  THEN
-                   CALL advec_s_ws( i, j, q, 'q', flux_s_q, & 
+                   CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
                                 diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
-                ELSE 
+                ELSE
                    CALL advec_s_pw( i, j, q )
                 ENDIF
@@ -782,5 +788,5 @@
 !
 !--          Nudging
-             IF ( nudging )  CALL nudge( i, j, simulated_time, 'q' ) 
+             IF ( nudging )  CALL nudge( i, j, simulated_time, 'q' )
 
 !
@@ -820,5 +826,5 @@
 
 !
-!--          If required, calculate prognostic equations for rain water content 
+!--          If required, calculate prognostic equations for rain water content
 !--          and rain drop concentration
              IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
@@ -829,8 +835,8 @@
                 THEN
                    IF ( ws_scheme_sca )  THEN
-                      CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr,       & 
+                      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 
+                   ELSE
                       CALL advec_s_pw( i, j, qr )
                    ENDIF
@@ -869,8 +875,8 @@
                 IF ( timestep_scheme(1:5) == 'runge' )  THEN
                    IF ( ws_scheme_sca )  THEN
-                      CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr,    & 
+                      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 
+                   ELSE
                       CALL advec_s_pw( i, j, nr )
                    ENDIF
@@ -907,5 +913,5 @@
 
           ENDIF
-          
+
 !
 !--       If required, compute prognostic equation for scalar
@@ -917,7 +923,7 @@
              THEN
                 IF ( ws_scheme_sca )  THEN
-                   CALL advec_s_ws( i, j, s, 's', flux_s_s, & 
+                   CALL advec_s_ws( i, j, s, 's', flux_s_s, &
                                 diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn )
-                ELSE 
+                ELSE
                    CALL advec_s_pw( i, j, s )
                 ENDIF
@@ -939,10 +945,10 @@
 !
 !--          Nudging, still need to be extended for scalars
-!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' ) 
+!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' )
 
 !
 !--          If required compute influence of large-scale subsidence/ascent.
-!--          Note, the last argument is of no meaning in this case, as it is 
-!--          only used in conjunction with large_scale_forcing, which is to 
+!--          Note, the last argument is of no meaning in this case, as it is
+!--          only used in conjunction with large_scale_forcing, which is to
 !--          date not implemented for scalars.
              IF ( large_scale_subsidence  .AND.                                &
@@ -980,7 +986,7 @@
              ENDIF
 
-          ENDIF          
-!
-!--       If required, compute prognostic equation for turbulent kinetic 
+          ENDIF
+!
+!--       If required, compute prognostic equation for turbulent kinetic
 !--       energy (TKE)
           IF ( .NOT. constant_diffusion )  THEN
@@ -990,5 +996,5 @@
              tend(:,j,i) = 0.0_wp
              IF ( timestep_scheme(1:5) == 'runge'  &
-                 .AND.  .NOT. use_upstream_for_tke )  THEN 
+                 .AND.  .NOT. use_upstream_for_tke )  THEN
                  IF ( ws_scheme_sca )  THEN
                      CALL advec_s_ws( i, j, e, 'e', flux_s_e, diss_s_e, &
@@ -1013,5 +1019,5 @@
 !
 !--          Additional sink term for flows through plant canopies
-             IF ( plant_canopy )  CALL pcm_tendency( i, j, 6 ) 
+             IF ( plant_canopy )  CALL pcm_tendency( i, j, 6 )
 
              CALL user_actions( i, j, 'e-tendency' )
@@ -1061,5 +1067,5 @@
 !> Version for vector machines
 !------------------------------------------------------------------------------!
- 
+
  SUBROUTINE prognostic_equations_vector
 
@@ -1176,5 +1182,5 @@
        IF ( ws_scheme_mom )  THEN
           CALL advec_v_ws
-       ELSE 
+       ELSE
           CALL advec_v_pw
        END IF
@@ -1389,5 +1395,5 @@
 !
 !--    Nudging
-       IF ( nudging )  CALL nudge( simulated_time, 'pt' ) 
+       IF ( nudging )  CALL nudge( simulated_time, 'pt' )
 
 !
@@ -1485,5 +1491,5 @@
 
        CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux )
-       
+
        CALL user_actions( 'sa-tendency' )
 
@@ -1573,5 +1579,5 @@
 
        CALL diffusion_s( q, qsws, qswst, wall_qflux )
-       
+
 !
 !--    Sink or source of humidity due to canopy elements
@@ -1586,5 +1592,5 @@
 !
 !--    Nudging
-       IF ( nudging )  CALL nudge( simulated_time, 'q' ) 
+       IF ( nudging )  CALL nudge( simulated_time, 'q' )
 
 !
@@ -1637,5 +1643,5 @@
 
 !
-!--    If required, calculate prognostic equations for rain water content 
+!--    If required, calculate prognostic equations for rain water content
 !--    and rain drop concentration
        IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
@@ -1675,6 +1681,4 @@
           CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux )
 
-          CALL user_actions( 'qr-tendency' )
-
 !
 !--       Prognostic equation for rain water content
@@ -1828,5 +1832,5 @@
 
        CALL diffusion_s( s, ssws, sswst, wall_sflux )
-       
+
 !
 !--    Sink or source of humidity due to canopy elements
@@ -1841,5 +1845,5 @@
 !
 !--    Nudging. Not implemented for scalars so far.
-!        IF ( nudging )  CALL nudge( simulated_time, 'q' ) 
+!        IF ( nudging )  CALL nudge( simulated_time, 'q' )
 
 !
@@ -1895,5 +1899,5 @@
     ENDIF
 !
-!-- If required, compute prognostic equation for turbulent kinetic 
+!-- If required, compute prognostic equation for turbulent kinetic
 !-- energy (TKE)
     IF ( .NOT. constant_diffusion )  THEN
@@ -1927,5 +1931,5 @@
                 IF ( ws_scheme_sca )  THEN
                    CALL advec_s_ws( e, 'e' )
-                ELSE 
+                ELSE
                    CALL advec_s_pw( e )
                 ENDIF
