Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 1064)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 1065)
@@ -20,4 +20,6 @@
 ! Current revisions:
 ! -----------------
+! Bugfix: It is not allowed to use cloud_scheme = seifert_beheng without 
+!         precipitation in order to save computational resources.
 !
 !
@@ -886,4 +888,11 @@
                         'loop_optimization = cache'
        CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 )
+    ENDIF 
+
+    IF ( cloud_physics  .AND. icloud_scheme == 0  .AND.                      &
+         .NOT. precipitation ) THEN
+       message_string = 'cloud_scheme = seifert_beheng requires ' // &
+                        'precipitation = .TRUE.'
+       CALL message( 'check_parameters', 'PA0363', 1, 2, 0, 6, 0 )
     ENDIF 
 
Index: palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_2d.f90	(revision 1064)
+++ palm/trunk/SOURCE/data_output_2d.f90	(revision 1065)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! -----------------
+! Bugfix: Output of cross sections of ql
 !
 ! Former revisions:
@@ -488,5 +489,5 @@
                       DO  i = nxlg, nxrg
                          DO  j = nysg, nyng
-                            DO  k = nzb, nz_do3d
+                            DO  k = nzb, nzt+1
                                local_pf(i,j,k) = ql(k,j,i) + qr(k,j,i)
                             ENDDO
@@ -501,5 +502,5 @@
                       DO  i = nxlg, nxrg
                          DO  j = nysg, nyng
-                            DO  k = nzb, nz_do3d
+                            DO  k = nzb, nzt+1
                                local_pf(i,j,k) = ql_av(k,j,i) + qr_av(k,j,i)
                             ENDDO
Index: palm/trunk/SOURCE/diffusion_e.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_e.f90	(revision 1064)
+++ palm/trunk/SOURCE/diffusion_e.f90	(revision 1065)
@@ -20,5 +20,7 @@
 ! Current revisions:
 ! -----------------
-!
+! Enabled the claculation of diss in case of turbulence = .TRUE. (parameterized 
+! effects of turbulence on autoconversion and accretion in two-moments cloud 
+! physics scheme). 
 !
 ! Former revisions:
@@ -183,5 +185,6 @@
 !--          Store dissipation if needed for calculating the sgs particle
 !--          velocities
-             IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
+             IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.               &
+                  turbulence )  THEN
                 DO  j = nys, nyn
                    DO  k = nzb_s_inner(j,i)+1, nzt
@@ -254,5 +257,6 @@
 !--          Store dissipation if needed for calculating the sgs particle
 !--          velocities
-             IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
+             IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.               &
+                  turbulence )  THEN
                 DO  j = nys, nyn
                    DO  k = nzb_s_inner(j,i)+1, nzt
@@ -268,5 +272,5 @@
 !
 !--    Boundary condition for dissipation
-       IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
+       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence )  THEN
           DO  i = nxl, nxr
              DO  j = nys, nyn
@@ -433,5 +437,6 @@
 !--                   Store dissipation if needed for calculating the sgs
 !--                   particle  velocities
-                      IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
+                      IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.      &
+                           turbulence )  THEN
                          diss(k,j,i) = dissipation
                       ENDIF
@@ -448,5 +453,5 @@
 !
 !--    Boundary condition for dissipation
-       IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
+       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence )  THEN
           !$acc kernels present( diss, nzb_s_inner )
           !$acc loop
@@ -539,5 +544,5 @@
 !
 !--    Store dissipation if needed for calculating the sgs particle velocities
-       IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
+       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence )  THEN
           DO  k = nzb_s_inner(j,i)+1, nzt
              diss(k,j,i) = dissipation(k)
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 1064)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 1065)
@@ -23,4 +23,5 @@
 ! Current revisions:
 ! ------------------
+! allocation of diss (dissipation rate) in case of turbulence = .TRUE. added
 !
 ! Former revisions:
@@ -467,5 +468,5 @@
 !-- 3D-array for storing the dissipation, needed for calculating the sgs
 !-- particle velocities
-    IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
+    IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence )  THEN
        ALLOCATE ( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     ELSE
Index: palm/trunk/SOURCE/init_cloud_physics.f90
===================================================================
--- palm/trunk/SOURCE/init_cloud_physics.f90	(revision 1064)
+++ palm/trunk/SOURCE/init_cloud_physics.f90	(revision 1065)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! -----------------
+! The Courant number of sedimentation can be controlled with c_sedimentation. 
 !
 ! Former revisions:
@@ -88,5 +89,6 @@
 !-- Calculate timestep according to precipitation
     IF ( icloud_scheme == 0  .AND.  precipitation )  THEN
-       dt_precipitation = MINVAL( dzu(nzb+2:nzt) ) / w_precipitation
+       dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) /        &
+                          w_precipitation
     ENDIF
 !
Index: palm/trunk/SOURCE/microphysics.f90
===================================================================
--- palm/trunk/SOURCE/microphysics.f90	(revision 1064)
+++ palm/trunk/SOURCE/microphysics.f90	(revision 1065)
@@ -4,4 +4,7 @@
 ! Current revisions:
 ! -----------------
+! Sedimentation process implemented according to Stevens and Seifert (2008).
+! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens 
+! and Stevens, 2010).
 !
 ! Former revisions:
@@ -241,8 +244,20 @@
 
        INTEGER ::  i, j, k
-       REAL    :: dr_min = 2.0E-6, dr_max = 1.0E-3
 
        DO  k = nzb_2d(j,i)+1, nzt
-          IF ( ( qr(k,j,i) > eps_sb ) .AND. ( nr(k,j,i) > eps_sb ) )  THEN
+
+          IF ( qr(k,j,i) <= eps_sb )  THEN
+             qr(k,j,i) = 0.0
+          ELSE
+!
+!--          Adjust number of raindrops to avoid nonlinear effects in 
+!--          sedimentation and evaporation of rain drops due to too small or
+!--          too big weights of rain drops (Stevens and Seifert, 2008).
+             IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
+                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
+             ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
+                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
+             ENDIF
+             xr(k) = hyrho(k) * qr(k,j,i) / nr(k,j,i)
 !
 !--          Weight averaged diameter of rain drops:
@@ -250,25 +265,7 @@
                        dpirho_l )**( 1.0 / 3.0 )
 !
-!--          Adjust number of raindrops to avoid nonlinear effects in 
-!--          sedimentation and evaporation of rain drops due to too small or
-!--          big diameters of rain drops (Stevens and Seifert, 2008).
-             IF ( dr(k) < dr_min )  THEN
-                nr(k,j,i) = qr(k,j,i) * hyrho(k) / dr_min**3 * dpirho_l
-                dr(k)     = dr_min
-             ELSEIF ( dr(k) > dr_max )  THEN
-                nr(k,j,i) = qr(k,j,i) * hyrho(k) / dr_max**3 * dpirho_l
-                dr(k)     = dr_max
-             ENDIF
-!
-!--          Mean weight of rain drops (Seifert and Beheng, 2006):
-             xr(k) = MIN( MAX( hyrho(k) * qr(k,j,i) / nr(k,j,i), xrmin ), xrmax)
-!
 !--          Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
 !--          Stevens and Seifert, 2008):
-             IF ( .NOT. mu_constant )  THEN
-                mu_r(k) = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr(k) - 1.4E-3 ) ) )
-             ELSE
-                mu_r(k) = mu_constant_value
-             ENDIF
+             mu_r(k) = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr(k) - 1.4E-3 ) ) )
 !
 !--          Slope parameter of gamma distribution (Seifert, 2008):
@@ -288,9 +285,11 @@
        USE control_parameters
        USE statistics
+       USE grid_variables
    
        IMPLICIT NONE
 
        INTEGER ::  i, j, k
-       REAL    ::  k_au, autocon, phi_au, tau_cloud, xc, nu_c
+       REAL    ::  k_au, autocon, phi_au, tau_cloud, xc, nu_c, rc,   &
+                   l_mix, re_lambda, alpha_cc, r_cc, sigma_cc, epsilon
 
        k_au = k_cc / ( 20.0 * x0 )
@@ -310,15 +309,45 @@
 !--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
 !--          (Use constant nu_c = 1.0 instead?)
-             nu_c      = 1580.0 * hyrho(k) * ql(k,j,i) - 0.28
+             nu_c      = 1.0 !MAX( 0.0, 1580.0 * hyrho(k) * ql(k,j,i) - 0.28 )
 !
 !--          Mean weight of cloud droplets:
              xc        = hyrho(k) * ql(k,j,i) / nc
 !
+!--          Parameterized turbulence effects on autoconversion (Seifert, 
+!--          Nuijens and Stevens, 2010)
+             IF ( turbulence )  THEN
+!
+!--             Weight averaged radius of cloud droplets:
+                rc = 0.5 * ( xc * dpirho_l )**( 1.0 / 3.0 )
+
+                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0 + a_3 * nu_c )
+                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0 + b_3 * nu_c )
+                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0 + c_3 * nu_c )
+!
+!--             Mixing length (neglecting distance to ground and stratification)
+                l_mix = ( dx * dy * dzu(k) )**( 1.0 / 3.0 )
+!
+!--             Limit dissipation rate according to Seifert, Nuijens and 
+!--             Stevens (2010)
+                epsilon = MIN( 0.06, diss(k,j,i) )
+!
+!--             Compute Taylor-microscale Reynolds number:
+                re_lambda = 6.0 / 11.0 * ( l_mix / c_const )**( 2.0 / 3.0 ) *  &
+                            SQRT( 15.0 / kin_vis_air ) * epsilon**( 1.0 / 6.0 )
+!
+!--             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 +                                          &
+                       epsilon * 1.0E4 * ( re_lambda * 1.0E-3 )**0.25 *        &
+                       ( alpha_cc * EXP( -1.0 * ( ( rc - r_cc ) /              &
+                       sigma_cc )**2 ) + beta_cc ) )
+             ENDIF
+!
 !--          Autoconversion rate (Seifert and Beheng, 2006):
-             autocon   = k_au * ( nu_c + 2.0 ) * ( nu_c + 4.0 ) /             &
-                         ( nu_c + 1.0 )**2 * ql(k,j,i)**2 * xc**2 *           &
-                         ( 1.0 + phi_au / ( 1.0 - tau_cloud )**2 ) *          &
+             autocon   = k_au * ( nu_c + 2.0 ) * ( nu_c + 4.0 ) /              &
+                         ( nu_c + 1.0 )**2 * ql(k,j,i)**2 * xc**2 *            &
+                         ( 1.0 + phi_au / ( 1.0 - tau_cloud )**2 ) *           &
                          rho_surface
-             autocon   = MIN( autocon, ql(k,j,i) / ( dt_3d *                  &
+             autocon   = MIN( autocon, ql(k,j,i) / ( dt_3d *                   &
                               weight_substep(intermediate_timestep_count) ) )
 !
@@ -346,5 +375,5 @@
 
        INTEGER ::  i, j, k
-       REAL    ::  accr, phi_ac, tau_cloud
+       REAL    ::  accr, phi_ac, tau_cloud, k_cr
 
        DO  k = nzb_2d(j,i)+1, nzt
@@ -356,12 +385,22 @@
 !--          Universal function for accretion process 
 !--          (Seifert and Beheng, 2001):
-             phi_ac    = tau_cloud / ( tau_cloud + 5.0E-5 ) 
-             phi_ac    = ( phi_ac**2 )**2
+             phi_ac = tau_cloud / ( tau_cloud + 5.0E-5 ) 
+             phi_ac = ( phi_ac**2 )**2
+!
+!--          Parameterized turbulence effects on autoconversion (Seifert, 
+!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to 
+!--          convert the dissipation (diss) from m2 s-3 to cm2 s-3.
+             IF ( turbulence )  THEN
+                k_cr = k_cr0 * ( 1.0 + 0.05 *                            &
+                                 MIN( 600.0, diss(k,j,i) * 1.0E4 )**0.25 )
+             ELSE
+                k_cr = k_cr0                       
+             ENDIF
 !
 !--          Accretion rate (Seifert and Beheng, 2006):
-             accr      = k_cr * ql(k,j,i) * qr(k,j,i) * phi_ac *              &
-                         SQRT( rho_surface * hyrho(k) )
-             accr      = MIN( accr, ql(k,j,i) / ( dt_3d *                     &
-                              weight_substep(intermediate_timestep_count) ) )
+             accr = k_cr * ql(k,j,i) * qr(k,j,i) * phi_ac *              &
+                    SQRT( rho_surface * hyrho(k) )
+             accr = MIN( accr, ql(k,j,i) / ( dt_3d *                     &
+                    weight_substep(intermediate_timestep_count) ) )
 !
 !--          Tendencies for q, qr, pt:
@@ -390,11 +429,11 @@
 
        DO  k = nzb_2d(j,i)+1, nzt
-          IF ( ( qr(k,j,i) > eps_sb )  .AND.  ( nr(k,j,i) > eps_sb ) )  THEN
+          IF ( qr(k,j,i) > eps_sb )  THEN
 !
 !--          Selfcollection rate (Seifert and Beheng, 2006):
 !--          pirho_l**( 1.0 / 3.0 ) is necessary to convert [lambda_r] = m-1 to 
 !--          kg**( 1.0 / 3.0 ).
-             phi_sc = ( 1.0 + kappa_rr / lambda_r(k) *                        &
-                      pirho_l**( 1.0 / 3.0 ) )**( -9 )
+             phi_sc = 1.0 !( 1.0 + kappa_rr / lambda_r(k) *                        &
+                      !pirho_l**( 1.0 / 3.0 ) )**( -9 )
 
              selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * phi_sc *               &
@@ -438,5 +477,5 @@
 
        DO  k = nzb_2d(j,i)+1, nzt
-          IF ( ( qr(k,j,i) > eps_sb )  .AND.  ( nr(k,j,i) > eps_sb ) )  THEN
+          IF ( qr(k,j,i) > eps_sb )  THEN
 !
 !--          Actual liquid water temperature:
@@ -555,42 +594,130 @@
        IMPLICIT NONE
 
-       INTEGER ::  i, j, k, n, n_substep
-       REAL    ::  sed_nr_tend, sed_qr_tend
+       INTEGER ::  i, j, k, k_run, n, n_substep
+       REAL    ::  c_run, d_max, d_mean, d_min, dt_sedi, flux, mean, z_run
+
+       REAL, DIMENSION(nzb:nzt) :: c_nr, c_qr, d_nr, d_qr, nr_slope, qr_slope, &
+                                   w_nr, w_qr
+!
+!--    Computation of sedimentation flux. Implementation according to Stevens
+!--    and Seifert (2008).
 
        IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0
 
-          sed_nr = 0.0
-          sed_qr = 0.0
-
-          DO  k = nzb_2d(j,i)+1, nzt
-             IF ( ( qr(k,j,i) > eps_sb )  .AND.  ( nr(k,j,i) > eps_sb ) )  THEN
-!
-!--             Sedimentation of rain water content and rain drop concentration
-!--             according to Stevens and Seifert (2008):
-                sed_nr(k) = MIN( 20.0, MAX( 0.1, a_term - b_term * ( 1.0 +    &
-                            c_term / lambda_r(k) )**( -1.0 *                  & 
-                            ( mu_r(k) + 1.0 ) ) ) ) * nr(k,j,i) 
-                sed_qr(k) = MIN( 20.0, MAX( 0.1, a_term - b_term * ( 1.0 +    &
-                            c_term / lambda_r(k) )**( -1.0 *                  &
-                            ( mu_r(k) + 4.0 ) ) ) ) * qr(k,j,i) * hyrho(k)
-!
-!--             Computation of rain rate
-                prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) *              &
-                             weight_substep(intermediate_timestep_count) 
-             ENDIF
-          ENDDO
-       
-          DO  k = nzb_2d(j,i)+1, nzt
-             sed_nr_tend = MAX( ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1),      &
-                                -nr(k,j,i) / ( dt_3d *                        &
-                                weight_substep(intermediate_timestep_count) ) )
-             sed_qr_tend = MAX( ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /     &
-                                hyrho(k),                                     &
-                                -qr(k,j,i) / ( dt_3d *                        &
-                                weight_substep(intermediate_timestep_count) ) )
-
-             tend_nr(k,j,i) = tend_nr(k,j,i) + sed_nr_tend
-             tend_qr(k,j,i) = tend_qr(k,j,i) + sed_qr_tend  
-          ENDDO
+       dt_sedi = dt_3d * weight_substep(intermediate_timestep_count)
+
+       w_nr = 0.0
+       w_qr = 0.0
+!
+!--    Compute velocities 
+       DO  k = nzb_s_inner(j,i)+1, nzt
+          IF ( qr(k,j,i) > eps_sb )  THEN
+             w_nr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 +          &
+                       c_term / lambda_r(k) )**( -1.0 * ( mu_r(k) + 1.0 ) ) ) )
+             w_qr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 +          &
+                       c_term / lambda_r(k) )**( -1.0 * ( mu_r(k) + 4.0 ) ) ) )
+          ELSE
+             w_nr(k) = 0.0
+             w_qr(k) = 0.0
+          ENDIF
+       ENDDO
+!
+!--    Adjust boundary values
+       w_nr(nzb_2d(j,i)) = w_nr(nzb_2d(j,i)+1)
+       w_qr(nzb_2d(j,i)) = w_qr(nzb_2d(j,i)+1)
+       w_nr(nzt) = w_nr(nzt-1)
+       w_qr(nzt) = w_qr(nzt-1)
+!
+!--    Compute Courant number
+       DO  k = nzb_s_inner(j,i)+1, nzt-1
+          c_nr(k) = 0.25 * ( w_nr(k-1) + 2.0 * w_nr(k) + w_nr(k+1) ) * &
+                    dt_sedi * ddzu(k)
+          c_qr(k) = 0.25 * ( w_qr(k-1) + 2.0 * w_qr(k) + w_qr(k+1) ) * &
+                    dt_sedi * ddzu(k)
+       ENDDO
+!
+!--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
+       IF ( limiter_sedimentation )  THEN
+
+          qr(nzb_s_inner(j,i),j,i) = 0.0
+          nr(nzb_s_inner(j,i),j,i) = 0.0
+          qr(nzt,j,i) = 0.0
+          nr(nzt,j,i) = 0.0
+
+          DO k = nzb_s_inner(j,i)+1, nzt-1
+             d_mean = 0.5 * ( qr(k+1,j,i) + qr(k-1,j,i) )
+             d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
+             d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
+
+             qr_slope(k) = SIGN(1.0, d_mean) * MIN ( 2.0 * d_min, 2.0 * d_max, &
+                                                     ABS( d_mean ) )
+
+             d_mean = 0.5 * ( nr(k+1,j,i) + nr(k-1,j,i) )
+             d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
+             d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
+
+             nr_slope(k) = SIGN(1.0, d_mean) * MIN ( 2.0 * d_min, 2.0 * d_max, &
+                                                     ABS( d_mean ) )
+          ENDDO
+
+       ELSE
+          nr_slope = 0.0
+          qr_slope = 0.0
+       ENDIF
+!
+!--    Compute sedimentation flux
+       DO  k = nzt-2, nzb_s_inner(j,i)+1, -1
+!
+!--       Sum up all rain drop number densities which contribute to the flux 
+!--       through k-1/2
+          flux  = 0.0
+          z_run = 0.0 ! height above z(k)
+          k_run = k
+          c_run = MIN( 1.0, c_nr(k) )
+          DO WHILE ( c_run > 0.0  .AND.  k_run <= nzt-1 )
+             flux  = flux + hyrho(k_run) *                                    &
+                     ( nr(k_run,j,i) + nr_slope(k_run) * ( 1.0 - c_run ) *    &
+                     0.5 ) * c_run * dzu(k_run)
+             z_run = z_run + dzu(k_run)
+             k_run = k_run + 1
+             c_run = MIN( 1.0, c_nr(k_run) - z_run * ddzu(k_run) )
+          ENDDO
+!
+!--       It is not allowed to sediment more rain drop number density than 
+!--       available
+          flux = MIN( flux,                                                   &
+                      hyrho(k) * dzu(k) * nr(k,j,i) + sed_nr(k+1) * dt_sedi )
+
+          sed_nr(k)      = flux / dt_sedi
+          tend_nr(k,j,i) = tend_nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *     &
+                                            ddzu(k+1) / hyrho(k)
+!
+!--       Sum up all rain water content which contributes to the flux 
+!--       through k-1/2
+          flux  = 0.0
+          z_run = 0.0 ! height above z(k)
+          k_run = k
+          c_run = MIN( 1.0, c_qr(k) )
+          DO WHILE ( c_run > 0.0  .AND.  k_run <= nzt-1 )
+             flux  = flux + hyrho(k_run) *                                    &
+                     ( qr(k_run,j,i) + qr_slope(k_run) * ( 1.0 - c_run ) *    &
+                     0.5 ) * c_run * dzu(k_run)
+             z_run = z_run + dzu(k_run)
+             k_run = k_run + 1
+             c_run = MIN( 1.0, c_qr(k_run) - z_run * ddzu(k_run) )
+          ENDDO
+!
+!--       It is not allowed to sediment more rain water content than available
+          flux = MIN( flux,                                                   &
+                      hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * dt_sedi )
+
+          sed_qr(k)      = flux / dt_sedi
+          tend_qr(k,j,i) = tend_qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *     &
+                                            ddzu(k+1) / hyrho(k)
+!
+!--       Compute the rain rate
+          prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) *                    &
+                       weight_substep(intermediate_timestep_count) 
+       ENDDO
 !
 !--    Precipitation amount
@@ -616,6 +743,5 @@
        IMPLICIT NONE 
  
-       REAL    ::  gamm, xx,                                        & 
-                   ser, tmp, x_gamm, y_gamm
+       REAL    ::  gamm, ser, tmp, x_gamm, xx, y_gamm
        INTEGER ::  j 
  
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 1064)
+++ palm/trunk/SOURCE/modules.f90	(revision 1065)
@@ -20,4 +20,8 @@
 ! Current revisions:
 ! -----------------
+! + c_sedimentation, limiter_sedimentation, turbulence, a_1, a_2, a_3, b_1, b_2, 
+! + b_3, c_1, c_2, c_3, beta_cc
+!
+! bottom boundary condition of qr, nr changed from Dirichlet to Neumann
 !
 ! Former revisions:
@@ -468,14 +472,27 @@
 
     LOGICAL ::  curvature_solution_effects = .FALSE., &
-                ventilation_effect = .FALSE., &
-                mu_constant = .FALSE.
-
-    REAL ::  a_vent = 0.78,      & ! coef. for ventilation effect 
-             a_term = 9.65,      & ! coef. for terminal velocity (m s-1)
-             b_vent = 0.308,     & ! coef. for ventilation effect
-             b_term = 9.8,       & ! coef. for terminal velocity (m s-1)
-             bfactor,            &
-             c_evap = 0.7,       & ! constant in evaporation
-             c_term = 600.0,     & ! coef. for terminal velocity (m-1)
+                limiter_sedimentation = .TRUE., &
+                ventilation_effect = .FALSE.
+                
+
+    REAL ::  a_1 = 8.69E-4,     & ! coef. in turb. parametrization (cm-2 s3)
+             a_2 = -7.38E-5,    & ! coef. in turb. parametrization (cm-2 s3)
+             a_3 = -1.40E-2,    & ! coef. in turb. parametrization
+             a_term = 9.65,     & ! coef. for terminal velocity (m s-1)
+             a_vent = 0.78,     & ! coef. for ventilation effect
+             b_1 = 11.45E-6,    & ! coef. in turb. parametrization (m)
+             b_2 = 9.68E-6,     & ! coef. in turb. parametrization (m)
+             b_3 = 0.62,        & ! coef. in turb. parametrization
+             b_term = 9.8,      & ! coef. for terminal velocity (m s-1)
+             b_vent = 0.308,    & ! coef. for ventilation effect
+             beta_cc = 3.09E-4, & ! coef. in turb. parametrization (cm-2 s3)
+             bfactor,           &
+             c_1 = 4.82E-6,     & ! coef. in turb. parametrization (m)
+             c_2 = 4.8E-6,      & ! coef. in turb. parametrization (m)
+             c_3 = 0.76,        & ! coef. in turb. parametrization
+             c_const = 0.93,    & ! const. in Taylor-microscale Reynolds number
+             c_evap = 0.7,      & ! constant in evaporation
+             c_sedimentation = 2.0, & ! Courant number of sedimentation process
+             c_term = 600.0,    & ! coef. for terminal velocity (m-1)
              cof(6) = (/ 76.18009172947146,      & ! coefficients in the 
                          -86.50532032941677,     & ! numerical
@@ -484,36 +501,35 @@
                          0.1208650973866179E-2,  &
                          -0.5395239384953E-5 /), &
-             cp = 1005.0,        & ! heat capacity of dry air (J kg-1 K-1) 
+             cp = 1005.0,       & ! heat capacity of dry air (J kg-1 K-1) 
              diff_coeff_l = 0.23E-4, & ! diffusivity of water vapor (m2 s-1)
              effective_coll_efficiency, &
-             eps_ros = 1.0E-4,   & ! accuracy of Rosenbrock method
-             eps_sb = 1.0E-20,   & ! threshold in two-moments scheme
-             k_cc = 4.44E09,     & ! const. rain-rain kernel (m3 kg-2 s-1) 
-             k_cr = 5.25,        & ! const. cloud-rain kernel (m3 kg-1 s-1) 
-             k_rr = 7.12,        & ! const. rain-rain kernel (m3 kg-1 s-1) 
-             k_br = 1000.,       & ! const. in breakup parametrization (m-1)
-             kappa_rr = 60.7,    & ! const. in collision kernel (kg-1/3)
+             eps_ros = 1.0E-4,  & ! accuracy of Rosenbrock method
+             eps_sb = 1.0E-20,  & ! threshold in two-moments scheme
+             k_cc = 9.44E09,    & ! const. rain-rain kernel (m3 kg-2 s-1) 
+             k_cr0 = 4.33,      & ! const. cloud-rain kernel (m3 kg-1 s-1) 
+             k_rr = 7.12,       & ! const. rain-rain kernel (m3 kg-1 s-1) 
+             k_br = 1000.,      & ! const. in breakup parametrization (m-1)
+             kappa_rr = 60.7,   & ! const. in collision kernel (kg-1/3)
              kin_vis_air = 1.4086E-5, & ! kin. viscosity of air (m2 s-1)
-             l_v = 2.5E+06,      & ! latent heat of vaporization (J kg-1)
+             l_v = 2.5E+06,     & ! latent heat of vaporization (J kg-1)
              l_d_cp, l_d_r, l_d_rv, & ! l_v / cp, l_v / r_d, l_v / r_v
              mass_of_solute = 1.0E-17, & ! soluted NaCl (kg)
              molecular_weight_of_solute = 0.05844, & ! mol. m. NaCl (kg mol-1)
              molecular_weight_of_water = 0.01801528, & ! mol. m. H2O (kg mol-1)
-             mu_constant_value = 0.0, & ! shape param. of gamma distribution
-             nc = 70.0E6,        & ! cloud droplet concentration 
+             nc = 70.0E6,       & ! cloud droplet concentration 
              prec_time_const = 0.001, & !coef. in Kessler scheme
-             pirho_l, dpirho_l,  & ! pi * rho_l / 6.0; 6.0 / ( pi * rho_l )
-             rho_l = 1.0E3,      & ! density of water (kg m-3)
-             ql_crit = 0.0005,   & ! coef. in Kessler scheme
-             r_d = 287.0,        & ! sp. gas const. dry air (J kg-1 K-1)
-             r_v = 461.51,       & ! sp. gas const. water vapor (J kg-1 K-1)
-             schmidt = 0.71,     & ! Schmidt number
-             schmidt_p_1d3,      & ! schmidt**( 1.0 / 3.0 )
+             pirho_l, dpirho_l, & ! pi * rho_l / 6.0; 6.0 / ( pi * rho_l )
+             rho_l = 1.0E3,     & ! density of water (kg m-3)
+             ql_crit = 0.0005,  & ! coef. in Kessler scheme
+             r_d = 287.0,       & ! sp. gas const. dry air (J kg-1 K-1)
+             r_v = 461.51,      & ! sp. gas const. water vapor (J kg-1 K-1)
+             schmidt = 0.71,    & ! Schmidt number
+             schmidt_p_1d3,     & ! schmidt**( 1.0 / 3.0 )
              stp = 2.5066282746310005, & ! parameter in gamma function
              thermal_conductivity_l = 2.43E-2, & ! therm. cond. air (J m-1 s-1 K-1)
-             vanthoff = 2.0,     & ! van't Hoff factor for NaCl
-             x0 = 2.6E-10,       & ! separating drop mass (kg)
-             xrmin = 2.6E-10,    & ! minimum rain drop size (kg)
-             xrmax = 5.0E-6,     & ! maximum rain drop site (kg)
+             vanthoff = 2.0,    & ! van't Hoff factor for NaCl
+             x0 = 2.6E-10,      & ! separating drop mass (kg)
+             xrmin = 2.6E-10,   & ! minimum rain drop size (kg)
+             xrmax = 5.0E-6,    & ! maximum rain drop site (kg)
              dt_precipitation = 100.0, & ! timestep precipitation (s)
              w_precipitation = 9.65      ! maximum terminal velocity (m s-1)
@@ -595,10 +611,10 @@
                              scalar_advec = 'ws-scheme'
     CHARACTER (LEN=20)   ::  bc_e_b = 'neumann', bc_lr = 'cyclic', &
-                             bc_nr_b = 'dirichlet', bc_nr_t = 'neumann', &
+                             bc_nr_b = 'neumann', bc_nr_t = 'neumann', &
                              bc_ns = 'cyclic', bc_p_b = 'neumann', &
                              bc_p_t = 'dirichlet', bc_pt_b = 'dirichlet', &
                              bc_pt_t = 'initial_gradient', &
                              bc_q_b = 'dirichlet', bc_q_t = 'neumann', &
-                             bc_qr_b = 'dirichlet', bc_qr_t = 'neumann',&
+                             bc_qr_b = 'neumann', bc_qr_t = 'neumann',&
                              bc_s_b = 'dirichlet', bc_s_t = 'neumann', &
                              bc_sa_t = 'neumann', &
@@ -710,5 +726,5 @@
                 do3d_compress = .FALSE., do_sum = .FALSE., &
                 dp_external = .FALSE., dp_smooth = .FALSE., &
-                drizzle = .TRUE., dt_fixed = .FALSE., &
+                drizzle = .FALSE., dt_fixed = .FALSE., &
                 dt_3d_reached, dt_3d_reached_l, exchange_mg = .FALSE., &
                 first_call_lpm = .TRUE., &
@@ -728,5 +744,6 @@
                 scalar_rayleigh_damping = .TRUE., sloping_surface = .FALSE., &
                 stop_dt = .FALSE., synchronous_exchange = .FALSE., &
-                terminate_run = .FALSE., turbulent_inflow = .FALSE., &
+                terminate_run = .FALSE., turbulence = .FALSE., &
+                turbulent_inflow = .FALSE., &
                 use_prescribed_profile_data = .FALSE., use_reference = .FALSE.,&
                 use_surface_fluxes = .FALSE., use_top_fluxes = .FALSE., &
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 1064)
+++ palm/trunk/SOURCE/parin.f90	(revision 1065)
@@ -20,4 +20,6 @@
 ! Current revisions:
 ! -----------------
+! +nc, c_sedimentation, limiter_sedimentation, turbulence
+! -mu_constant, mu_constant_value
 !
 ! Former revisions:
@@ -199,6 +201,7 @@
              call_psolver_at_all_substeps, canopy_mode, canyon_height, &
              canyon_width_x, canyon_width_y, canyon_wall_left, &
-             canyon_wall_south, cfl_factor, cloud_droplets, cloud_physics, &
-             cloud_scheme, collective_wait, conserve_volume_flow, &
+             canyon_wall_south, c_sedimentation, cfl_factor, cloud_droplets, &
+             cloud_physics, cloud_scheme, collective_wait, &
+             conserve_volume_flow, &
              conserve_volume_flow_mode, coupling_start_time, cthf, &
              curvature_solution_effects, cycle_mg, damp_level_1d, &
@@ -212,9 +215,8 @@
              initializing_actions, km_constant, lad_surface, &
              lad_vertical_gradient, lad_vertical_gradient_level, &
-             leaf_surface_concentration, &
+             leaf_surface_concentration, limiter_sedimentation, &
              loop_optimization, masking_method, mg_cycles, &
              mg_switch_to_pe0_level, mixing_length_1d, momentum_advec, &
-             mu_constant, mu_constant_value, &
-             netcdf_precision, neutral, ngsrb, nr_surface, &
+             nc, netcdf_precision, neutral, ngsrb, nr_surface, &
              nr_surface_initial_change, nr_vertical_gradient, &
              nr_vertical_gradient_level, nsor, &
@@ -241,5 +243,5 @@
              topography, topography_grid_convention, top_heatflux, &
              top_momentumflux_u, top_momentumflux_v, top_salinityflux, &
-             turbulent_inflow, ug_surface, ug_vertical_gradient, &
+             turbulence, turbulent_inflow, ug_surface, ug_vertical_gradient, &
              ug_vertical_gradient_level, use_surface_fluxes, &
              use_top_fluxes, use_ug_for_galilei_tr, use_upstream_for_tke, &
Index: palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- palm/trunk/SOURCE/read_var_list.f90	(revision 1064)
+++ palm/trunk/SOURCE/read_var_list.f90	(revision 1065)
@@ -20,4 +20,6 @@
 ! Current revisions:
 ! ------------------
+! +nc, c_sedimentation, limiter_sedimentation, turbulence
+! -mu_constant, mu_constant_value
 !
 ! Former revisions:
@@ -179,5 +181,5 @@
 !-- Make version number check first
     READ ( 13 )  version_on_file
-    binary_version = '3.6'
+    binary_version = '3.8'
     IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
        WRITE( message_string, * ) 'version mismatch concerning control ', &
@@ -342,4 +344,6 @@
           CASE ( 'canyon_wall_south' )
              READ ( 13 )  canyon_wall_south
+          CASE ( 'c_sedimentation' )
+             READ ( 13 )  c_sedimentation
           CASE ( 'cfl_factor' )
              READ ( 13 )  cfl_factor
@@ -447,4 +451,6 @@
           CASE ( 'leaf_surface_concentration' )
              READ ( 13 )  leaf_surface_concentration
+          CASE ( 'limiter_sedimentation' )
+             READ ( 13 )  limiter_sedimentation
           CASE ( 'loop_optimization' )
              READ ( 13 )  loop_optimization
@@ -464,8 +470,6 @@
           CASE ( 'momentum_advec' )
              READ ( 13 )  momentum_advec
-          CASE ( 'mu_constant' )
-             READ ( 13 )  mu_constant
-          CASE ( 'mu_constant_value' )
-             READ ( 13 )  mu_constant_value
+          CASE ( 'nc' )
+             READ ( 13 )  nc
           CASE ( 'netcdf_precision' )
              READ ( 13 )  netcdf_precision
@@ -678,4 +682,6 @@
           CASE ( 'tsc' )
              READ ( 13 )  tsc
+          CASE ( 'turbulence' )
+             READ ( 13 )  turbulence
           CASE ( 'turbulent_inflow' )
              READ ( 13 )  turbulent_inflow
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 1064)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 1065)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! -----------------
+! exchange of diss (dissipation rate) in case of turbulence = .TRUE. added
 !
 ! Former revisions:
@@ -255,5 +256,5 @@
              CALL exchange_horiz( ql_vp, nbgp )
           ENDIF
-          IF ( wang_kernel )  CALL exchange_horiz( diss, nbgp )
+          IF ( wang_kernel  .OR.  turbulence )  CALL exchange_horiz( diss, nbgp )
 
           CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
Index: palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- palm/trunk/SOURCE/write_var_list.f90	(revision 1064)
+++ palm/trunk/SOURCE/write_var_list.f90	(revision 1065)
@@ -20,4 +20,6 @@
 ! Current revisions:
 ! -----------------
+! +nc, c_sedimentation, turbulence, limiter_sedimentation
+! -mu_constant, mu_constant_value
 !
 ! Former revisions:
@@ -161,5 +163,5 @@
 
 
-    binary_version = '3.6'
+    binary_version = '3.8'
 
     WRITE ( 14 )  binary_version
@@ -264,4 +266,6 @@
     WRITE ( 14 )  'canyon_wall_south             '
     WRITE ( 14 )  canyon_wall_south
+    WRITE ( 14 )  'c_sedimentation               '
+    WRITE ( 14 )  c_sedimentation
     WRITE ( 14 )  'cfl_factor                    '
     WRITE ( 14 )  cfl_factor
@@ -368,4 +372,6 @@
     WRITE ( 14 )  'leaf_surface_concentration    '
     WRITE ( 14 )  leaf_surface_concentration
+    WRITE ( 14 )  'limiter_sedimentation         '
+    WRITE ( 14 )  limiter_sedimentation
     WRITE ( 14 )  'loop_optimization             '
     WRITE ( 14 )  loop_optimization
@@ -384,8 +390,6 @@
     WRITE ( 14 )  'momentum_advec                '
     WRITE ( 14 )  momentum_advec
-    WRITE ( 14 )  'mu_constant                   '
-    WRITE ( 14 )  mu_constant
-    WRITE ( 14 )  'mu_constant_value             '
-    WRITE ( 14 )  mu_constant_value
+    WRITE ( 14 )  'nc                            '
+    WRITE ( 14 )  nc
     WRITE ( 14 )  'netcdf_precision              '
     WRITE ( 14 )  netcdf_precision
@@ -596,4 +600,6 @@
     WRITE ( 14 )  'tsc                           '
     WRITE ( 14 )  tsc
+    WRITE ( 14 )  'turbulence                    '
+    WRITE ( 14 )  turbulence
     WRITE ( 14 )  'turbulent_inflow              '
     WRITE ( 14 )  turbulent_inflow
