Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 1047)
+++ palm/trunk/SOURCE/Makefile	(revision 1053)
@@ -20,5 +20,5 @@
 # Current revisions:
 # ------------------
-# 
+# +microphysics for two-moment cloud physics
 #
 # Former revisions:
@@ -159,5 +159,6 @@
 	user_parin.f90 user_read_restart_data.f90 \
 	user_spectra.f90 user_statistics.f90 wall_fluxes.f90 \
-	write_3d_binary.f90 write_compressed.f90 write_var_list.f90
+	write_3d_binary.f90 write_compressed.f90 write_var_list.f90 \
+	microphysics.f90
 
 OBJS =  advec_s_bc.o advec_s_pw.o advec_s_up.o advec_u_pw.o advec_u_up.o \
@@ -205,5 +206,6 @@
 	user_lpm_init.o user_lpm_set_attributes.o user_module.o user_parin.o \
 	user_read_restart_data.o user_spectra.o user_statistics.o \
-	wall_fluxes.o write_3d_binary.o write_compressed.o write_var_list.o
+	wall_fluxes.o write_3d_binary.o write_compressed.o write_var_list.o \
+	microphysics.o
 
 CC = cc
@@ -346,5 +348,6 @@
 	eqn_state_seawater.o impact_of_latent_heat.o \
         plant_canopy_model.o production_e.o \
-	subsidence.o user_actions.o
+	subsidence.o user_actions.o \
+	microphysics.o
 random_gauss.o: random_function.o
 read_3d_binary.o: modules.o random_function.o
@@ -394,2 +397,3 @@
 write_compressed.o: modules.o
 write_var_list.o: modules.o
+microphysics.o: modules.o
Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 1047)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 1053)
@@ -20,4 +20,7 @@
 ! Current revisions:
 ! ------------------
+! necessary expansions according to the two new prognostic equations (nr, qr) 
+! of the two-moment cloud physics scheme:
+! +flux_l_*, flux_s_*, diss_l_*, diss_s_*, sums_ws*s_ws_l 
 !
 ! Former revisions:
@@ -210,4 +213,11 @@
           ENDIF
 
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+             ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
+             ALLOCATE( sums_wsnrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
+             sums_wsqrs_ws_l = 0.0
+             sums_wsnrs_ws_l = 0.0
+          ENDIF
+
           IF ( ocean )  THEN
              ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) )
@@ -253,9 +263,20 @@
 
              IF ( humidity .OR. passive_scalar )  THEN
-                ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1),         &
+                ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1),          &
                           diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
-                ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
+                ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1),  &
                           diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
              ENDIF
+
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1),         &
+                          diss_s_qr(nzb+1:nzt,0:threads_per_task-1),         &
+                          flux_s_nr(nzb+1:nzt,0:threads_per_task-1),         &
+                          diss_s_nr(nzb+1:nzt,0:threads_per_task-1) )
+                ALLOCATE( flux_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
+                          diss_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
+                          flux_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
+                          diss_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) 
+             END IF
 
              IF ( ocean )  THEN
@@ -297,4 +318,8 @@
           sums_wspts_ws_l = 0.0
           IF ( humidity .OR. passive_scalar )  sums_wsqs_ws_l = 0.0
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+             sums_wsqrs_ws_l = 0.0
+             sums_wsnrs_ws_l = 0.0
+          ENDIF
           IF ( ocean )  sums_wssas_ws_l = 0.0
 
@@ -307,6 +332,6 @@
 ! Scalar advection - Call for grid point i,j
 !------------------------------------------------------------------------------!
-    SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char,swap_flux_y_local,  &
-                              swap_diss_y_local, swap_flux_x_local, &
+    SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char, swap_flux_y_local,  &
+                              swap_diss_y_local, swap_flux_x_local,  &
                               swap_diss_x_local, i_omp, tn )
 
@@ -727,4 +752,21 @@
                                  * weight_substep(intermediate_timestep_count)
              ENDDO
+
+          CASE ( 'qr' )
+
+             DO  k = nzb, nzt
+                sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_ws_l(k,tn) +               &
+                                      ( flux_t(k) + diss_t(k) )                &
+                                 * weight_substep(intermediate_timestep_count)
+             ENDDO
+
+          CASE ( 'nr' )
+
+             DO  k = nzb, nzt
+                sums_wsnrs_ws_l(k,tn)  = sums_wsnrs_ws_l(k,tn) +               &
+                                      ( flux_t(k) + diss_t(k) )                &
+                                 * weight_substep(intermediate_timestep_count)
+             ENDDO
+
 
          END SELECT
Index: palm/trunk/SOURCE/average_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/average_3d_data.f90	(revision 1047)
+++ palm/trunk/SOURCE/average_3d_data.f90	(revision 1053)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! averaging of nr, qr added
 !
 ! Former revisions:
@@ -117,4 +117,13 @@
              ENDDO
 
+          CASE ( 'nr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      nr_av(k,j,i) = nr_av(k,j,i) / REAL( average_count_3d )
+                   ENDDO
+                ENDDO
+             ENDDO
+
           CASE ( 'p' )
              DO  i = nxlg, nxrg
@@ -207,4 +216,13 @@
              ENDDO
 
+          CASE ( 'qr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      qr_av(k,j,i) = qr_av(k,j,i) / REAL( average_count_3d )
+                   ENDDO
+                ENDDO
+             ENDDO
+
           CASE ( 'qv' )
              DO  i = nxlg, nxrg
Index: palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- palm/trunk/SOURCE/boundary_conds.f90	(revision 1047)
+++ palm/trunk/SOURCE/boundary_conds.f90	(revision 1053)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-!
+! boundary conditions for the two new prognostic equations (nr, qr) of the 
+! two-moment cloud scheme
 !  
 ! Former revisions:
@@ -216,4 +217,43 @@
 !--       Top boundary
           q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
+
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+!             
+!--          Surface conditions for constant_humidity_flux
+             IF ( ibc_qr_b == 0 ) THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      qr_p(nzb_s_inner(j,i),j,i) = qr(nzb_s_inner(j,i),j,i)
+                   ENDDO
+                ENDDO
+             ELSE
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      qr_p(nzb_s_inner(j,i),j,i) = qr_p(nzb_s_inner(j,i)+1,j,i)
+                   ENDDO
+                ENDDO
+             ENDIF
+!
+!--          Top boundary
+             qr_p(nzt+1,:,:) = qr_p(nzt,:,:) + bc_qr_t_val * dzu(nzt+1)
+!             
+!--          Surface conditions for constant_humidity_flux
+             IF ( ibc_nr_b == 0 ) THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      nr_p(nzb_s_inner(j,i),j,i) = nr(nzb_s_inner(j,i),j,i)
+                   ENDDO
+                ENDDO
+             ELSE
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      nr_p(nzb_s_inner(j,i),j,i) = nr_p(nzb_s_inner(j,i)+1,j,i)
+                   ENDDO
+                ENDDO
+             ENDIF
+!
+!--          Top boundary
+             nr_p(nzt+1,:,:) = nr_p(nzt,:,:) + bc_nr_t_val * dzu(nzt+1)
+          ENDIF
 
        ENDIF
@@ -242,17 +282,41 @@
           pt_p(:,nys-1,:)     = pt_p(:,nys,:)
           IF ( .NOT. constant_diffusion     )  e_p(:,nys-1,:) = e_p(:,nys,:)
-          IF ( humidity .OR. passive_scalar )  q_p(:,nys-1,:) = q_p(:,nys,:)
+          IF ( humidity .OR. passive_scalar )  THEN
+             q_p(:,nys-1,:) = q_p(:,nys,:)
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                qr_p(:,nys-1,:) = qr_p(:,nys,:)
+                nr_p(:,nys-1,:) = nr_p(:,nys,:)
+             ENDIF
+          ENDIF
        ELSEIF ( outflow_n )  THEN
           pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
           IF ( .NOT. constant_diffusion     )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
-          IF ( humidity .OR. passive_scalar )  q_p(:,nyn+1,:) = q_p(:,nyn,:)
+          IF ( humidity .OR. passive_scalar )  THEN
+             q_p(:,nyn+1,:) = q_p(:,nyn,:)
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
+                nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
+             ENDIF
+          ENDIF
        ELSEIF ( outflow_l )  THEN
           pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
           IF ( .NOT. constant_diffusion     )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
-          IF ( humidity .OR. passive_scalar )  q_p(:,:,nxl-1) = q_p(:,:,nxl)
+          IF ( humidity .OR. passive_scalar )  THEN
+             q_p(:,:,nxl-1) = q_p(:,:,nxl)
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
+                nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
+             ENDIF
+          ENDIF
        ELSEIF ( outflow_r )  THEN
           pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
           IF ( .NOT. constant_diffusion     )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
-          IF ( humidity .OR. passive_scalar )  q_p(:,:,nxr+1) = q_p(:,:,nxr)
+          IF ( humidity .OR. passive_scalar )  THEN
+             q_p(:,:,nxr+1) = q_p(:,:,nxr)
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
+                nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
+             ENDIF
+          ENDIF
        ENDIF
 
Index: palm/trunk/SOURCE/calc_liquid_water_content.f90
===================================================================
--- palm/trunk/SOURCE/calc_liquid_water_content.f90	(revision 1047)
+++ palm/trunk/SOURCE/calc_liquid_water_content.f90	(revision 1053)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! -----------------
+! description expanded to the two-moment cloud scheme
 !
 ! Former revisions:
@@ -51,5 +52,7 @@
 ! Description:
 ! ------------
-! Calculation of the liquid water content (0%-or-100%-scheme)
+! Calculation of the liquid water content (0%-or-100%-scheme). This scheme is 
+! used by the one and the two moment cloud physics scheme. Using the two moment 
+! scheme, this calculation results in the cloud water content.
 !------------------------------------------------------------------------------!
 
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 1047)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 1053)
@@ -20,5 +20,11 @@
 ! Current revisions:
 ! -----------------
-!
+! necessary changes for the new two-moment cloud physics scheme added:
+! - check cloud physics scheme (Kessler or Seifert and Beheng)
+! - plant_canopy is not allowed
+! - currently, only cache loop_optimization is allowed 
+! - initial profiles of nr, qr
+! - boundary condition of nr, qr
+! - check output quantities (qr, nr, prr)
 !
 ! Former revisions:
@@ -617,5 +623,15 @@
 
     ENDIF
-
+!
+!-- Check cloud scheme
+    IF ( cloud_scheme == 'seifert_beheng' )  THEN
+       icloud_scheme = 0
+    ELSEIF ( cloud_scheme == 'kessler' )  THEN
+       icloud_scheme = 1
+    ELSE
+       message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
+                        TRIM( cloud_scheme ) // '"'
+       CALL message( 'check_parameters', 'PA0357', 1, 2, 0, 6, 0 )
+    ENDIF
 !
 !-- Check whether there are any illegal values
@@ -853,4 +869,17 @@
     ENDIF 
 
+    IF ( plant_canopy  .AND.  cloud_physics  .AND.  icloud_scheme == 0 ) THEN
+       message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' //  &
+                        ' seifert_beheng'
+       CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
+    ENDIF 
+
+    IF ( loop_optimization /= 'cache' .AND.  cloud_physics  .AND.            &
+         icloud_scheme == 0 ) THEN
+       message_string = 'cloud_scheme = seifert_beheng requires ' // &
+                        'loop_optimization = cache'
+       CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 )
+    ENDIF 
+
 !
 !-- In case of no model continuation run, check initialising parameters and
@@ -861,5 +890,14 @@
 !--    Initial profiles for 1D and 3D model, respectively (u,v further below)
        pt_init = pt_surface
-       IF ( humidity )        q_init  = q_surface
+       IF ( humidity )  THEN
+          q_init  = q_surface
+!
+!--       It is not allowed to choose initial profiles of rain water content 
+!--       and rain drop concentration. They are set to 0.0.
+          IF ( cloud_physics .AND. icloud_scheme == 0 ) THEN
+             qr_init = 0.0
+             nr_init = 0.0
+          ENDIF
+       ENDIF
        IF ( ocean )           sa_init = sa_surface
        IF ( passive_scalar )  q_init  = s_surface
@@ -1167,8 +1205,84 @@
           ENDIF
 
-!
-!--       Store humidity gradient at the top boundary for possile Neumann
-!--       boundary condition
-          bc_q_t_val = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+
+             i = 1
+             gradient = 0.0
+             qr_vertical_gradient_level_ind(1) = 0
+             DO  k = 1, nzt+1
+                IF ( i < 11 ) THEN
+                   IF ( qr_vertical_gradient_level(i) < zu(k)  .AND. &
+                        qr_vertical_gradient_level(i) >= 0.0 )  THEN
+                      gradient = qr_vertical_gradient(i) / 100.0
+                      qr_vertical_gradient_level_ind(i) = k - 1
+                      i = i + 1
+                   ENDIF
+                ENDIF
+                IF ( gradient /= 0.0 )  THEN
+                   IF ( k /= 1 )  THEN
+                      qr_init(k) = qr_init(k-1) + dzu(k) * gradient
+                   ELSE
+                      qr_init(k) = qr_init(k-1) + 0.5 * dzu(k) * gradient
+                   ENDIF
+                ELSE
+                   qr_init(k) = qr_init(k-1)
+                ENDIF
+!
+!--             Avoid negative rain water content
+                IF ( qr_init(k) < 0.0 )  THEN
+                   qr_init(k) = 0.0
+                ENDIF
+             ENDDO
+!
+!--          In case of no given rain water content gradients, choose zero gradient
+!--          conditions
+             IF ( qr_vertical_gradient_level(1) == -1.0 )  THEN
+                qr_vertical_gradient_level(1) = 0.0
+             ENDIF
+
+             i = 1
+             gradient = 0.0
+             nr_vertical_gradient_level_ind(1) = 0
+             DO  k = 1, nzt+1
+                IF ( i < 11 ) THEN
+                   IF ( nr_vertical_gradient_level(i) < zu(k)  .AND. &
+                        nr_vertical_gradient_level(i) >= 0.0 )  THEN
+                      gradient = nr_vertical_gradient(i) / 100.0
+                      nr_vertical_gradient_level_ind(i) = k - 1
+                      i = i + 1
+                   ENDIF
+                ENDIF
+                IF ( gradient /= 0.0 )  THEN
+                   IF ( k /= 1 )  THEN
+                      nr_init(k) = nr_init(k-1) + dzu(k) * gradient
+                   ELSE
+                      nr_init(k) = nr_init(k-1) + 0.5 * dzu(k) * gradient
+                   ENDIF
+                ELSE
+                   nr_init(k) = nr_init(k-1)
+                ENDIF
+!
+!--             Avoid negative rain water content
+                IF ( nr_init(k) < 0.0 )  THEN
+                   nr_init(k) = 0.0
+                ENDIF
+             ENDDO
+!
+!--          In case of no given rain water content gradients, choose zero gradient
+!--          conditions
+             IF ( nr_vertical_gradient_level(1) == -1.0 )  THEN
+                nr_vertical_gradient_level(1) = 0.0
+             ENDIF
+
+          ENDIF
+!
+!--       Store humidity, rain water content and rain drop concentration 
+!--       gradient at the top boundary for possile Neumann boundary condition
+          bc_q_t_val  = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
+            
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+             bc_qr_t_val = ( qr_init(nzt+1) - qr_init(nzt) ) / dzu(nzt+1)
+             bc_nr_t_val = ( nr_init(nzt+1) - nr_init(nzt) ) / dzu(nzt+1)
+          ENDIF
 
        ENDIF
@@ -1344,16 +1458,20 @@
 !-- Set wind speed in the Galilei-transformed system
     IF ( galilei_transformation )  THEN
-       IF ( use_ug_for_galilei_tr .AND.                &
-            ug_vertical_gradient_level(1) == 0.0 .AND. & 
-            vg_vertical_gradient_level(1) == 0.0 )  THEN
+       IF ( use_ug_for_galilei_tr .AND.                  &
+            ug_vertical_gradient_level(1) == 0.0  .AND.  &
+            ug_vertical_gradient(1) == 0.0  .AND.        & 
+            vg_vertical_gradient_level(1) == 0.0  .AND.  &
+            vg_vertical_gradient(1) == 0.0 )  THEN
           u_gtrans = ug_surface * 0.6
           v_gtrans = vg_surface * 0.6
-       ELSEIF ( use_ug_for_galilei_tr .AND.                &
-                ug_vertical_gradient_level(1) /= 0.0 )  THEN
+       ELSEIF ( use_ug_for_galilei_tr  .AND.                  &
+                ( ug_vertical_gradient_level(1) /= 0.0  .OR.  &
+                ug_vertical_gradient(1) /= 0.0 ) )  THEN
           message_string = 'baroclinicity (ug) not allowed simultaneously' // &
                            ' with galilei transformation'
           CALL message( 'check_parameters', 'PA0046', 1, 2, 0, 6, 0 )
-       ELSEIF ( use_ug_for_galilei_tr .AND.                &
-                vg_vertical_gradient_level(1) /= 0.0 )  THEN
+       ELSEIF ( use_ug_for_galilei_tr  .AND.                  &
+                ( vg_vertical_gradient_level(1) /= 0.0  .OR.  &
+                vg_vertical_gradient(1) /= 0.0 ) )  THEN
           message_string = 'baroclinicity (vg) not allowed simultaneously' // &
                            ' with galilei transformation'
@@ -1650,5 +1768,44 @@
           CALL message( 'check_parameters', 'PA0074', 1, 2, 0, 6, 0 )
        ENDIF
-       
+
+       IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+          IF ( bc_qr_b == 'dirichlet' )  THEN
+             ibc_qr_b = 0
+          ELSEIF ( bc_qr_b == 'neumann' )  THEN
+             ibc_qr_b = 1
+          ELSE
+             message_string = 'unknown boundary condition: bc_qr_b ="' // &
+                              TRIM( bc_qr_b ) // '"'
+             CALL message( 'check_parameters', 'PA0352', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( bc_qr_t == 'dirichlet' )  THEN
+             ibc_qr_t = 0
+          ELSEIF ( bc_qr_t == 'neumann' )  THEN
+             ibc_qr_t = 1
+          ELSE
+             message_string = 'unknown boundary condition: bc_qr_t ="' // &
+                              TRIM( bc_qr_t ) // '"'
+             CALL message( 'check_parameters', 'PA0353', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( bc_nr_b == 'dirichlet' )  THEN
+             ibc_nr_b = 0
+          ELSEIF ( bc_nr_b == 'neumann' )  THEN
+             ibc_nr_b = 1
+          ELSE
+             message_string = 'unknown boundary condition: bc_nr_b ="' // &
+                              TRIM( bc_nr_b ) // '"'
+             CALL message( 'check_parameters', 'PA0355', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( bc_nr_t == 'dirichlet' )  THEN
+             ibc_nr_t = 0
+          ELSEIF ( bc_nr_t == 'neumann' )  THEN
+             ibc_nr_t = 1
+          ELSE
+             message_string = 'unknown boundary condition: bc_nr_t ="' // &
+                              TRIM( bc_nr_t ) // '"'
+             CALL message( 'check_parameters', 'PA0356', 1, 2, 0, 6, 0 )
+          ENDIF
+       ENDIF       
+
     ENDIF
 
@@ -2483,4 +2640,78 @@
              hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
 
+          CASE ( 'nr' )
+             IF ( .NOT. cloud_physics )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for cloud_physics = .FALSE.'
+                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
+             ELSEIF ( icloud_scheme /= 0 )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for cloud_scheme /= seifert_beheng'
+                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 73
+                dopr_unit(i)  = '1/m3'
+                hom(:,2,73,:)  = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'qr' )
+             IF ( .NOT. cloud_physics )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for cloud_physics = .FALSE.'
+                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
+             ELSEIF ( icloud_scheme /= 0 )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for cloud_scheme /= seifert_beheng'
+                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 74
+                dopr_unit(i)  = 'kg/kg'
+                hom(:,2,74,:)  = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'qc' )
+             IF ( .NOT. cloud_physics )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for cloud_physics = .FALSE.'
+                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
+             ELSEIF ( icloud_scheme /= 0 )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for cloud_scheme /= seifert_beheng'
+                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 75
+                dopr_unit(i)  = 'kg/kg'
+                hom(:,2,75,:)  = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'prr' )
+             IF ( .NOT. cloud_physics )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for cloud_physics = .FALSE.'
+                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
+             ELSEIF ( icloud_scheme /= 0 )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for cloud_scheme /= seifert_beheng'
+                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
+             ELSEIF ( .NOT. precipitation )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for precipitation = .FALSE.'
+                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
+
+             ELSE
+                dopr_index(i) = 76
+                dopr_unit(i)  = 'kg/kg m/s'
+                hom(:,2,76,:)  = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
           CASE DEFAULT
 
@@ -2572,4 +2803,16 @@
              unit = 'K'
 
+          CASE ( 'nr' )
+             IF ( .NOT. cloud_physics )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' // &
+                         'res cloud_physics = .TRUE.'
+                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
+             ELSEIF ( icloud_scheme /= 0 )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' // &
+                         'res cloud_scheme = seifert_beheng'
+                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = '1/m3'
+
           CASE ( 'pc', 'pr' )
              IF ( .NOT. particle_advection )  THEN
@@ -2581,4 +2824,20 @@
              IF ( TRIM( var ) == 'pr' )  unit = 'm'
 
+          CASE ( 'prr' )
+             IF ( .NOT. cloud_physics )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' // &
+                         'res cloud_physics = .TRUE.'
+                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
+             ELSEIF ( icloud_scheme /= 0 )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' // &
+                         'res cloud_scheme = seifert_beheng'
+                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
+             ELSEIF ( .NOT. precipitation )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' // &
+                                 'res precipitation = .TRUE.'
+                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = 'kg/kg m/s'
+
           CASE ( 'q', 'vpt' )
              IF ( .NOT. humidity )  THEN
@@ -2590,4 +2849,16 @@
              IF ( TRIM( var ) == 'vpt' )  unit = 'K'
 
+          CASE ( 'qc' )
+             IF ( .NOT. cloud_physics )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' // &
+                         'res cloud_physics = .TRUE.'
+                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
+             ELSEIF ( icloud_scheme /= 0 ) THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' // &
+                         'res cloud_scheme = seifert_beheng'
+                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = 'kg/kg'
+
           CASE ( 'ql' )
              IF ( .NOT. ( cloud_physics  .OR.  cloud_droplets ) )  THEN
@@ -2607,4 +2878,16 @@
              IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
              IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
+
+          CASE ( 'qr' )
+             IF ( .NOT. cloud_physics )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' // &
+                         'res cloud_physics = .TRUE.'
+                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
+             ELSEIF ( icloud_scheme /= 0 ) THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' // &
+                         'res cloud_scheme = seifert_beheng'
+                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = 'kg/kg'
 
           CASE ( 'qv' )
Index: palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_2d.f90	(revision 1047)
+++ palm/trunk/SOURCE/data_output_2d.f90	(revision 1053)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! +qr, nr, qc and cross sections
 !
 ! Former revisions:
@@ -290,4 +290,12 @@
                 two_d = .TRUE.
                 level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'nr_xy', 'nr_xz', 'nr_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => nr
+                ELSE
+                   to_be_resorted => nr_av
+                ENDIF
+                IF ( mode == 'xy' )  level_z = zu
 
              CASE ( 'p_xy', 'p_xz', 'p_yz' )
@@ -376,22 +384,63 @@
 
              CASE ( 'prr*_xy' )        ! 2d-array
-                IF ( av == 0 )  THEN
-                   CALL exchange_horiz_2d( precipitation_rate )
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) =  precipitation_rate(j,i)
-                      ENDDO
-                   ENDDO
-                ELSE
-                   CALL exchange_horiz_2d( precipitation_rate_av )
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) =  precipitation_rate_av(j,i)
-                      ENDDO
-                   ENDDO
+                IF ( icloud_scheme == 1 )  THEN
+                   IF ( av == 0 )  THEN
+                      CALL exchange_horiz_2d( precipitation_rate )
+                      DO  i = nxlg, nxrg
+                         DO  j = nysg, nyng
+                            local_pf(i,j,nzb+1) =  precipitation_rate(j,i)
+                         ENDDO
+                      ENDDO
+                   ELSE
+                      CALL exchange_horiz_2d( precipitation_rate_av )
+                      DO  i = nxlg, nxrg
+                         DO  j = nysg, nyng
+                            local_pf(i,j,nzb+1) =  precipitation_rate_av(j,i)
+                         ENDDO
+                      ENDDO
+                   ENDIF
+                ELSE
+                   IF ( av == 0 )  THEN
+                      CALL exchange_horiz_2d( prr(nzb+1,:,:) )
+                      DO  i = nxlg, nxrg
+                         DO  j = nysg, nyng
+                            local_pf(i,j,nzb+1) = prr(nzb+1,j,i) * hyrho(nzb+1)
+                         ENDDO
+                      ENDDO
+                   ELSE
+                      CALL exchange_horiz_2d( prr_av(nzb+1,:,:) )
+                      DO  i = nxlg, nxrg
+                         DO  j = nysg, nyng
+                            local_pf(i,j,nzb+1) = prr_av(nzb+1,j,i) * hyrho(nzb+1)
+                         ENDDO
+                      ENDDO
+                   ENDIF
                 ENDIF
                 resorted = .TRUE.
                 two_d = .TRUE.
                 level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'prr_xy', 'prr_xz', 'prr_yz' )
+                IF ( av == 0 )  THEN
+                   CALL exchange_horiz( prr, nbgp )
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            local_pf(i,j,k) = prr(k,j,i)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ELSE
+                   CALL exchange_horiz( prr_av, nbgp )
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            local_pf(i,j,k) = prr_av(k,j,i)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                IF ( mode == 'xy' )  level_z = zu
 
              CASE ( 'pt_xy', 'pt_xz', 'pt_yz' )
@@ -424,9 +473,39 @@
                 IF ( mode == 'xy' )  level_z = zu
 
+             CASE ( 'qc_xy', 'qc_xz', 'qc_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => ql
+                ELSE
+                   to_be_resorted => ql_av
+                ENDIF
+                IF ( mode == 'xy' )  level_z = zu
+
              CASE ( 'ql_xy', 'ql_xz', 'ql_yz' )
                 IF ( av == 0 )  THEN
-                   to_be_resorted => ql
-                ELSE
-                   to_be_resorted => ql_av
+                   IF ( icloud_scheme == 0 )  THEN
+                      DO  i = nxlg, nxrg
+                         DO  j = nysg, nyng
+                            DO  k = nzb, nz_do3d
+                               local_pf(i,j,k) = ql(k,j,i) + qr(k,j,i)
+                            ENDDO
+                         ENDDO
+                      ENDDO
+                      resorted = .TRUE.
+                   ELSE
+                      to_be_resorted => ql
+                   ENDIF
+                ELSE
+                   IF ( icloud_scheme == 0 )  THEN
+                      DO  i = nxlg, nxrg
+                         DO  j = nysg, nyng
+                            DO  k = nzb, nz_do3d
+                               local_pf(i,j,k) = ql_av(k,j,i) + qr_av(k,j,i)
+                            ENDDO
+                         ENDDO
+                      ENDDO
+                      resorted = .TRUE.
+                   ELSE
+                      to_be_resorted => ql_av
+                   ENDIF
                 ENDIF
                 IF ( mode == 'xy' )  level_z = zu
@@ -478,4 +557,12 @@
                    CALL exchange_horiz( ql_vp_av, nbgp )
                    to_be_resorted => ql_vp
+                ENDIF
+                IF ( mode == 'xy' )  level_z = zu
+
+             CASE ( 'qr_xy', 'qr_xz', 'qr_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => qr
+                ELSE
+                   to_be_resorted => qr_av
                 ENDIF
                 IF ( mode == 'xy' )  level_z = zu
Index: palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_3d.f90	(revision 1047)
+++ palm/trunk/SOURCE/data_output_3d.f90	(revision 1053)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! +nr, qr, prr, qc and averaged quantities
 !
 ! Former revisions:
@@ -207,4 +207,11 @@
              ENDIF
 
+          CASE ( 'nr' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => nr
+             ELSE
+                to_be_resorted => nr_av
+             ENDIF
+
           CASE ( 'p' )
              IF ( av == 0 )  THEN
@@ -270,4 +277,26 @@
              ENDIF
 
+          CASE ( 'prr' )
+             IF ( av == 0 )  THEN
+                CALL exchange_horiz( prr, nbgp )
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      DO  k = nzb, nzt+1
+                         local_pf(i,j,k) = prr(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ELSE
+                CALL exchange_horiz( prr_av, nbgp )
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      DO  k = nzb, nzt+1
+                         local_pf(i,j,k) = prr_av(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             resorted = .TRUE.
+
           CASE ( 'pt' )
              IF ( av == 0 )  THEN
@@ -297,9 +326,38 @@
              ENDIF
 
+          CASE ( 'qc' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => ql
+             ELSE
+                to_be_resorted => ql_av
+             ENDIF
+
           CASE ( 'ql' )
              IF ( av == 0 )  THEN
-                to_be_resorted => ql
-             ELSE
-                to_be_resorted => ql_av
+                IF ( icloud_scheme == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nz_do3d
+                            local_pf(i,j,k) = ql(k,j,i) + qr(k,j,i)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                   resorted = .TRUE.
+                ELSE
+                   to_be_resorted => ql
+                ENDIF
+             ELSE
+                IF ( icloud_scheme == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nz_do3d
+                            local_pf(i,j,k) = ql_av(k,j,i) + qr_av(k,j,i)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                   resorted = .TRUE.
+                ELSE
+                   to_be_resorted => ql_av
+                ENDIF
              ENDIF
 
@@ -344,4 +402,11 @@
                 CALL exchange_horiz( ql_vp_av, nbgp )
                 to_be_resorted => ql_vp_av
+             ENDIF
+
+          CASE ( 'qr' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => qr
+             ELSE
+                to_be_resorted => qr_av
              ENDIF
 
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 1047)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 1053)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-!
+! necessary additions for two-moment cloud physics scheme:
+! +nr, qr, qc, prr
 !
 ! Former revisions:
@@ -748,14 +749,40 @@
 
 !
-!--             Buoyancy flux, water flux, humidity flux and liquid water
-!--             content
+!--             Buoyancy flux, water flux, humidity flux, liquid water
+!--             content, rain drop concentration and rain water content
                 IF ( humidity )  THEN
                    IF ( cloud_physics .OR. cloud_droplets )  THEN
-                      pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
+                      pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +        &
                                     vpt(k+1,j,i) - hom(k+1,1,44,sr) )
                       sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
                                                           rmask(j,i,sr)
-                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
-                                                          rmask(j,i,sr)
+                      IF ( .NOT. cloud_droplets )  THEN
+                         pts = 0.5 *                                          &
+                              ( ( q(k,j,i) - ql(k,j,i) ) -                    &
+                              hom(k,1,42,sr) +                                &
+                              ( q(k+1,j,i) - ql(k+1,j,i) ) -                  &
+                              hom(k+1,1,42,sr) )
+                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * &
+                                                             rmask(j,i,sr)
+                         IF ( icloud_scheme == 0  )  THEN
+                            sums_l(k,54,tn) = sums_l(k,54,tn) + ( ql(k,j,i) + &
+                                                                qr(k,j,i) ) * &
+                                                                rmask(j,i,sr)
+                            sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *   &
+                                                                rmask(j,i,sr)
+                            sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) *   &
+                                                                rmask(j,i,sr)
+                            sums_l(k,75,tn) = sums_l(k,75,tn) + ql(k,j,i) *   &
+                                                                rmask(j,i,sr)
+                            sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *  &
+                                                                rmask(j,i,sr)
+                         ELSE
+                            sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *   &
+                                                                rmask(j,i,sr)
+                         ENDIF
+                      ELSE
+                         sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *      &
+                                                             rmask(j,i,sr)
+                      ENDIF
                    ELSE
                       IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
@@ -765,5 +792,5 @@
                                                              rmask(j,i,sr)
                       ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
-                         sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) * &
+                         sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) *  &
                                              sums_l(k,17,tn) +      &
                                           0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
@@ -1105,4 +1132,8 @@
        hom(:,1,71,sr) = sums(:,71)     ! prho
        hom(:,1,72,sr) = hyp * 1E-4     ! hyp in dbar
+       hom(:,1,73,sr) = sums(:,73)     ! nr
+       hom(:,1,74,sr) = sums(:,74)     ! qr
+       hom(:,1,75,sr) = sums(:,75)     ! qc
+       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
 
        hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 1047)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 1053)
@@ -23,5 +23,8 @@
 ! Current revisions:
 ! ------------------
-!
+! allocation and initialisation of necessary data arrays for the two-moment 
+! cloud physics scheme the two new prognostic equations (nr, qr):
+! +dr, lambda_r, mu_r, sed_*, xr, *s, *sws, *swst, *, *_p, t*_m, *_1, *_2, *_3, 
+! +tend_*, prr
 !
 ! Former revisions:
@@ -353,4 +356,5 @@
 
           IF ( cloud_physics ) THEN
+
 !
 !--          Liquid water content
@@ -364,4 +368,51 @@
              ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg), &
                        precipitation_rate(nysg:nyng,nxlg:nxrg) )
+
+             IF ( icloud_scheme == 0 )  THEN
+!
+!--             1D-rain sedimentation fluxes and rain drop size distribution
+!--             properties 
+                ALLOCATE ( dr(nzb:nzt+1), lambda_r(nzb:nzt+1),  &
+                           mu_r(nzb:nzt+1), sed_q(nzb:nzt+1),   &
+                           sed_qr(nzb:nzt+1), sed_nr(nzb:nzt+1),& 
+                           xr(nzb:nzt+1) ) 
+!
+!--             2D-rain water content and rain drop concentration arrays
+                ALLOCATE ( qrs(nysg:nyng,nxlg:nxrg),   &
+                   qrsws(nysg:nyng,nxlg:nxrg),         &
+                   qrswst(nysg:nyng,nxlg:nxrg),        &
+                   nrs(nysg:nyng,nxlg:nxrg),           &
+                   nrsws(nysg:nyng,nxlg:nxrg),         &
+                   nrswst(nysg:nyng,nxlg:nxrg) )
+!
+!--             3D-rain water content, rain drop concentration arrays
+#if defined( __nopointer )
+                ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),         &
+                          nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                          qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),         &
+                          qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                          tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),      &
+                          tqr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+#else
+                ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                          nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                          nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                          qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                          qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                          qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+#endif
+!
+!--             3D-tendency arrays
+                ALLOCATE( tend_nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
+                          tend_pt(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
+                          tend_q(nzb:nzt+1,nysg:nyng,nxlg:nxrg),     &
+                          tend_qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+!
+!--             3d-precipitation rate
+                IF ( precipitation )  THEN
+                   ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+
+             ENDIF
           ENDIF
 
@@ -508,6 +559,14 @@
     IF ( humidity  .OR.  passive_scalar )  THEN
        q => q_1;  q_p => q_2;  tq_m => q_3
-       IF ( humidity )        vpt  => vpt_1
-       IF ( cloud_physics )   ql   => ql_1
+       IF ( humidity )  THEN
+          vpt  => vpt_1    
+          IF ( cloud_physics )  THEN
+             ql => ql_1
+             IF ( icloud_scheme == 0 )  THEN
+                qr => qr_1;  qr_p  => qr_2;  tqr_m  => qr_3
+                nr => nr_1;  nr_p  => nr_2;  tnr_m  => nr_3
+             ENDIF
+          ENDIF
+       ENDIF
        IF ( cloud_droplets )  THEN
           ql   => ql_1
@@ -564,4 +623,12 @@
                 ENDDO
              ENDDO
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      qr(:,j,i) = qr_init
+                      nr(:,j,i) = nr_init
+                   ENDDO
+                ENDDO
+             ENDIF
           ENDIF
 
@@ -606,5 +673,11 @@
 !--       This could actually be computed more accurately in the 1D model.
 !--       Update when opportunity arises!
-          IF ( humidity  .OR.  passive_scalar )  qs = 0.0
+          IF ( humidity  .OR.  passive_scalar )  THEN
+             qs = 0.0
+             IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
+                qrs = 0.0
+                nrs = 0.0
+             ENDIF
+          ENDIF
 
 !
@@ -671,4 +744,12 @@
                 ENDDO
              ENDDO
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      qr(:,j,i) = qr_init
+                      nr(:,j,i) = nr_init
+                   ENDDO
+                ENDDO
+             ENDIF
           ENDIF
 
@@ -811,4 +892,12 @@
 !--       Determine the near-surface water flux
           IF ( humidity  .OR.  passive_scalar )  THEN
+             IF ( cloud_physics .AND. icloud_scheme )  THEN
+                IF ( constant_waterflux_qr )  THEN
+                   qrsws = surface_waterflux_qr
+                ENDIF
+                IF (constant_waterflux_nr )  THEN
+                   nrsws = surface_waterflux_nr
+                ENDIF
+             ENDIF
              IF ( constant_waterflux )  THEN
                 qsws   = surface_waterflux
@@ -842,5 +931,11 @@
              tswst = top_heatflux
 
-             IF ( humidity  .OR.  passive_scalar )  qswst = 0.0
+             IF ( humidity  .OR.  passive_scalar )  THEN
+                qswst = 0.0
+                IF ( cloud_physics  .AND.  icloud_scheme == 0 ) THEN
+                   nrswst = 0.0
+                   qrswst = 0.0
+                ENDIF
+             ENDIF
 
              IF ( ocean )  THEN
@@ -877,4 +972,12 @@
           IF ( humidity  .OR.  passive_scalar )  THEN
              IF ( .NOT. constant_waterflux )  qsws   = 0.0
+             IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
+                IF ( .NOT. constant_waterflux_qr )  THEN
+                   qrsws = 0.0
+                ENDIF
+                IF ( .NOT. constant_waterflux_nr )  THEN
+                   nrsws = 0.0
+                ENDIF
+             ENDIF
           ENDIF
 
@@ -925,4 +1028,11 @@
             q_surface_initial_change /= 0.0 )  THEN
           q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+             IF ( qr_surface_initial_change /= 0.0 )  THEN
+                qr(nzb,:,:) = qr(nzb,:,:) + qr_surface_initial_change
+             ELSEIF ( nr_surface_initial_change /= 0.0 ) THEN
+                nr(nzb,:,:) = nr(nzb,:,:) + nr_surface_initial_change
+             ENDIF
+          ENDIF
        ENDIF
 
@@ -939,4 +1049,10 @@
           tq_m = 0.0
           q_p = q
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+             tqr_m = 0.0
+             qr_p = qr
+             tnr_m = 0.0
+             nr_p = nr
+          ENDIF
        ENDIF
 
@@ -1109,5 +1225,11 @@
 !--    including ghost points)
        e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
-       IF ( humidity  .OR.  passive_scalar )  q_p = q
+       IF ( humidity  .OR.  passive_scalar )  THEN
+          q_p = q
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+             qr_p = qr
+             nr_p = nr
+          ENDIF
+       ENDIF
        IF ( ocean )  sa_p = sa
 
@@ -1117,5 +1239,11 @@
 !--    there before they are set. 
        te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
-       IF ( humidity  .OR.  passive_scalar )  tq_m = 0.0
+       IF ( humidity  .OR.  passive_scalar )  THEN
+          tq_m = 0.0
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+             tqr_m = 0.0
+             tnr_m = 0.0
+          ENDIF
+       ENDIF
        IF ( ocean )  tsa_m = 0.0
 
Index: palm/trunk/SOURCE/init_cloud_physics.f90
===================================================================
--- palm/trunk/SOURCE/init_cloud_physics.f90	(revision 1047)
+++ palm/trunk/SOURCE/init_cloud_physics.f90	(revision 1053)
@@ -20,5 +20,9 @@
 ! Current revisions:
 ! -----------------
-! 
+! calculation of the maximum timestep according to the terminal velocity of rain 
+! drops in the two moment cloud scheme
+!
+! calculation of frequently used constants (pirho_l, dpirho_l, schmidt_p_1d3, 
+! hyrho) 
 !
 ! Former revisions:
@@ -66,5 +70,6 @@
     REAL    ::  t_surface
 
-    ALLOCATE( hyp(nzb:nzt+1), pt_d_t(nzb:nzt+1), t_d_pt(nzb:nzt+1) )
+    ALLOCATE( hyp(nzb:nzt+1), pt_d_t(nzb:nzt+1), t_d_pt(nzb:nzt+1),  &
+              hyrho(nzb:nzt+1) )
 
 !
@@ -73,4 +78,14 @@
     l_d_r  = l_v / r_d
     l_d_rv = l_v / r_v
+
+    schmidt_p_1d3 = schmidt**( 1.0 / 3.0 )
+
+    pirho_l  = pi * rho_l / 6.0
+    dpirho_l = 1.0 / pirho_l 
+!
+!-- Calculate timestep according to precipitation
+    IF ( icloud_scheme == 0  .AND.  precipitation )  THEN
+       dt_precipitation = MINVAL( dzu(nzb+2:nzt) ) / w_precipitation
+    ENDIF
 !
 !-- Calculate factor used in equation for droplet growth by condensation
@@ -94,5 +109,6 @@
                    ( (t_surface - g/cp * zu(k)) / t_surface )**(1.0/0.286)
        pt_d_t(k) = ( 100000.0 / hyp(k) )**0.286
-       t_d_pt(k) = 1.0 / pt_d_t(k)       
+       t_d_pt(k) = 1.0 / pt_d_t(k)
+       hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
     ENDDO
 
Index: palm/trunk/SOURCE/microphysics.f90
===================================================================
--- palm/trunk/SOURCE/microphysics.f90	(revision 1053)
+++ palm/trunk/SOURCE/microphysics.f90	(revision 1053)
@@ -0,0 +1,638 @@
+ MODULE microphysics_mod
+
+!------------------------------------------------------------------------------!
+! Current revisions:
+! -----------------
+! initial revision
+!
+! Former revisions:
+! -----------------
+! $Id$
+! 
+! Description:
+! ------------
+! Calculate cloud microphysics according to the two moment bulk 
+! scheme by Seifert and Beheng (2006).
+!------------------------------------------------------------------------------!
+
+    PRIVATE
+    PUBLIC dsd_properties, autoconversion, accretion, selfcollection_breakup, &
+           evaporation_rain, sedimentation_cloud, sedimentation_rain
+
+    INTERFACE dsd_properties
+       MODULE PROCEDURE dsd_properties
+       MODULE PROCEDURE dsd_properties_ij
+    END INTERFACE dsd_properties
+
+    INTERFACE autoconversion
+       MODULE PROCEDURE autoconversion
+       MODULE PROCEDURE autoconversion_ij
+    END INTERFACE autoconversion
+
+    INTERFACE accretion
+       MODULE PROCEDURE accretion
+       MODULE PROCEDURE accretion_ij
+    END INTERFACE accretion
+
+    INTERFACE selfcollection_breakup
+       MODULE PROCEDURE selfcollection_breakup
+       MODULE PROCEDURE selfcollection_breakup_ij
+    END INTERFACE selfcollection_breakup
+
+    INTERFACE evaporation_rain
+       MODULE PROCEDURE evaporation_rain
+       MODULE PROCEDURE evaporation_rain_ij
+    END INTERFACE evaporation_rain
+
+    INTERFACE sedimentation_cloud
+       MODULE PROCEDURE sedimentation_cloud
+       MODULE PROCEDURE sedimentation_cloud_ij
+    END INTERFACE sedimentation_cloud
+ 
+    INTERFACE sedimentation_rain
+       MODULE PROCEDURE sedimentation_rain
+       MODULE PROCEDURE sedimentation_rain_ij
+    END INTERFACE sedimentation_rain
+
+ CONTAINS
+
+
+!------------------------------------------------------------------------------!
+! Call for all grid points
+!------------------------------------------------------------------------------!
+    SUBROUTINE dsd_properties
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE indices
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  dqdt_precip
+
+ 
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_2d(j,i)+1, nzt
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE dsd_properties
+
+    SUBROUTINE autoconversion
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE indices
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  dqdt_precip
+
+ 
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_2d(j,i)+1, nzt
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE autoconversion
+
+    SUBROUTINE accretion
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE indices
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  dqdt_precip
+
+ 
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_2d(j,i)+1, nzt
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE accretion
+
+    SUBROUTINE selfcollection_breakup
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE indices
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  dqdt_precip
+
+ 
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_2d(j,i)+1, nzt
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE selfcollection_breakup
+
+    SUBROUTINE evaporation_rain
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE indices
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  dqdt_precip
+
+ 
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_2d(j,i)+1, nzt
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE evaporation_rain
+
+    SUBROUTINE sedimentation_cloud
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE indices
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  dqdt_precip
+
+ 
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_2d(j,i)+1, nzt
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE sedimentation_cloud
+
+    SUBROUTINE sedimentation_rain
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE indices
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  dqdt_precip
+
+ 
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_2d(j,i)+1, nzt
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+
+    END SUBROUTINE sedimentation_rain
+
+
+!------------------------------------------------------------------------------!
+! Call for grid point i,j
+!------------------------------------------------------------------------------!
+    SUBROUTINE dsd_properties_ij( i, j )
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE indices
+       USE control_parameters
+       USE user
+   
+       IMPLICIT NONE
+
+       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
+!
+!--          Weight averaged diameter of rain drops:
+             dr(k) = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) *                     &
+                       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
+!
+!--          Slope parameter of gamma distribution (Seifert, 2008):
+             lambda_r(k) = ( ( mu_r(k) + 3.0 ) * ( mu_r(k) + 2.0 ) *          &
+                             ( mu_r(k) + 1.0 ) )**( 1.0 / 3.0 ) / dr(k)
+          ENDIF
+       ENDDO
+
+    END SUBROUTINE dsd_properties_ij
+
+    SUBROUTINE autoconversion_ij( i, j )
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE indices
+       USE control_parameters
+       USE statistics
+   
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  k_au, autocon, phi_au, tau_cloud, xc, nu_c
+
+       k_au = k_cc / ( 20.0 * x0 )
+
+       DO  k = nzb_2d(j,i)+1, nzt
+
+          IF ( ql(k,j,i) > 0.0 )  THEN
+!
+!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
+!--          (1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) ))
+             tau_cloud = 1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) + 1.0E-20 )
+!
+!--          Universal function for autoconversion process 
+!--          (Seifert and Beheng, 2006):
+             phi_au    = 600.0 * tau_cloud**0.68 * ( 1.0 - tau_cloud**0.68 )**3
+!
+!--          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
+!
+!--          Mean weight of cloud droplets:
+             xc        = hyrho(k) * ql(k,j,i) / nc
+!
+!--          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 ) *          &
+                         rho_surface
+             autocon   = MIN( autocon, ql(k,j,i) / ( dt_3d *                  &
+                              weight_substep(intermediate_timestep_count) ) )
+!
+!--          Tendencies for q, qr, nr, pt:
+             tend_qr(k,j,i) = tend_qr(k,j,i) + autocon
+             tend_q(k,j,i)  = tend_q(k,j,i) - autocon
+             tend_nr(k,j,i) = tend_nr(k,j,i) + autocon / x0 * hyrho(k)
+             tend_pt(k,j,i) = tend_pt(k,j,i) + autocon * l_d_cp * pt_d_t(k) 
+          ENDIF
+
+       ENDDO
+
+    END SUBROUTINE autoconversion_ij
+
+    SUBROUTINE accretion_ij( i, j )
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE indices
+       USE control_parameters
+       USE statistics
+       
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  accr, phi_ac, tau_cloud
+
+       DO  k = nzb_2d(j,i)+1, nzt
+          IF ( ( ql(k,j,i) > 0.0 )  .AND.  ( qr(k,j,i) > eps_sb ) )  THEN
+!
+!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
+             tau_cloud = 1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) + 1.0E-20) 
+!
+!--          Universal function for accretion process 
+!--          (Seifert and Beheng, 2001):
+             phi_ac    = tau_cloud / ( tau_cloud + 5.0E-5 ) 
+             phi_ac    = ( phi_ac**2 )**2
+!
+!--          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) ) )
+!
+!--          Tendencies for q, qr, pt:
+             tend_qr(k,j,i) = tend_qr(k,j,i) + accr
+             tend_q(k,j,i)  = tend_q(k,j,i) - accr
+             tend_pt(k,j,i) = tend_pt(k,j,i) + accr * l_d_cp * pt_d_t(k) 
+          ENDIF
+       ENDDO
+
+    END SUBROUTINE accretion_ij
+
+
+    SUBROUTINE selfcollection_breakup_ij( i, j )
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE indices
+       USE control_parameters
+       USE statistics
+   
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  selfcoll, breakup, phi_br, phi_sc
+
+       DO  k = nzb_2d(j,i)+1, nzt
+          IF ( ( qr(k,j,i) > eps_sb )  .AND.  ( nr(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 )
+
+             selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * phi_sc *               &
+                        SQRT( hyrho(k) * rho_surface )
+!
+!--          Collisional breakup rate (Seifert, 2008):
+             IF ( dr(k) >= 0.3E-3 )  THEN
+                phi_br  = k_br * ( dr(k) - 1.1E-3 )
+                breakup = selfcoll * ( phi_br + 1.0 )
+             ELSE
+                breakup = 0.0
+             ENDIF
+
+             selfcoll =  MAX( breakup - selfcoll, -nr(k,j,i) / ( dt_3d *      &
+                              weight_substep(intermediate_timestep_count) ) )
+!
+!--          Tendency for nr:
+             tend_nr(k,j,i) = tend_nr(k,j,i) + selfcoll
+          ENDIF          
+       ENDDO
+
+    END SUBROUTINE selfcollection_breakup_ij
+
+    SUBROUTINE evaporation_rain_ij( i, j )
+!
+!--    Evaporation of precipitable water. Condensation is neglected for 
+!--    precipitable water.
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE indices
+       USE control_parameters
+       USE statistics
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  evap, alpha, e_s, q_s, t_l, sat, temp, g_evap, f_vent, &
+                   mu_r_2, mu_r_5d2, nr_0
+
+       DO  k = nzb_2d(j,i)+1, nzt
+          IF ( ( qr(k,j,i) > eps_sb )  .AND.  ( nr(k,j,i) > eps_sb ) )  THEN
+!
+!--          Actual liquid water temperature:
+             t_l = t_d_pt(k) * pt(k,j,i)
+!
+!--          Saturation vapor pressure at t_l:
+             e_s = 610.78 * EXP( 17.269 * ( t_l - 273.16 ) / ( t_l - 35.86 ) )
+!
+!--          Computation of saturation humidity:
+             q_s = 0.622 * e_s / ( hyp(k) - 0.378 * e_s )
+             alpha = 0.622 * l_d_r * l_d_cp / ( t_l * t_l )
+             q_s = q_s * ( 1.0 + alpha * q(k,j,i) ) / ( 1.0 + alpha * q_s )
+!
+!--          Oversaturation:
+             sat = MIN( 0.0, ( q(k,j,i) - ql(k,j,i) ) / q_s - 1.0 )
+!
+!--          Actual temperature:
+             temp = t_l + l_d_cp * ql(k,j,i)
+!
+!--          
+             g_evap = ( l_v / ( r_v * temp ) - 1.0 ) * l_v /                  &
+                      ( thermal_conductivity_l * temp ) + rho_l * r_v * temp /&
+                      ( diff_coeff_l * e_s )
+             g_evap = 1.0 / g_evap
+!
+!--          Compute ventilation factor and intercept parameter 
+!--          (Seifert and Beheng, 2006; Seifert, 2008):
+             IF ( ventilation_effect )  THEN
+                mu_r_2   = mu_r(k) + 2.0
+                mu_r_5d2 = mu_r(k) + 2.5 
+                f_vent = a_vent * gamm( mu_r_2 ) *                            &
+                         lambda_r(k)**( -mu_r_2 ) +                           &
+                         b_vent * schmidt_p_1d3 *                             &
+                         SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *    &
+                         lambda_r(k)**( -mu_r_5d2 ) *                         &
+                         ( 1.0 - 0.5 * ( b_term / a_term ) *                  &
+                         ( lambda_r(k) /                                      &
+                         (       c_term + lambda_r(k) ) )**mu_r_5d2 -         &
+                                 0.125 * ( b_term / a_term )**2 *             &
+                         ( lambda_r(k) /                                      &
+                         ( 2.0 * c_term + lambda_r(k) ) )**mu_r_5d2 -         &
+                                 0.0625 * ( b_term / a_term )**3 *            &
+                         ( lambda_r(k) /                                      &
+                         ( 3.0 * c_term + lambda_r(k) ) )**mu_r_5d2 -         &
+                                 0.0390625 * ( b_term / a_term )**4 *         &
+                         ( lambda_r(k) /                                      &
+                         ( 4.0 * c_term + lambda_r(k) ) )**mu_r_5d2 )
+                nr_0   = nr(k,j,i) * lambda_r(k)**( mu_r(k) + 1.0 ) /         &
+                         gamm( mu_r(k) + 1.0 ) 
+             ELSE
+                f_vent = 1.0
+                nr_0   = nr(k,j,i) * dr(k)
+             ENDIF
+!
+!--          Evaporation rate of rain water content (Seifert and Beheng, 2006):
+             evap = 2.0 * pi * nr_0 * g_evap * f_vent * sat /    &
+                    hyrho(k)
+             evap = MAX( evap, -qr(k,j,i) / ( dt_3d *                         &
+                         weight_substep(intermediate_timestep_count) ) )
+!
+!--          Tendencies for q, qr, nr, pt:
+             tend_qr(k,j,i) = tend_qr(k,j,i) + evap
+             tend_q(k,j,i)  = tend_q(k,j,i) - evap
+             tend_nr(k,j,i) = tend_nr(k,j,i) + c_evap * evap / xr(k) * hyrho(k)
+             tend_pt(k,j,i) = tend_pt(k,j,i) + evap * l_d_cp * pt_d_t(k) 
+          ENDIF          
+       ENDDO
+
+    END SUBROUTINE evaporation_rain_ij
+
+    SUBROUTINE sedimentation_cloud_ij( i, j )
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE indices
+       USE control_parameters
+       
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  sed_q_const, sigma_gc = 1.3, k_st = 1.2E8
+!
+!--    Sedimentation of cloud droplets (Heus et al., 2010):
+       sed_q_const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))**( 2.0 / 3.0 ) *    &
+                     EXP( 5.0 * LOG( sigma_gc )**2 )
+
+       sed_q = 0.0
+
+       DO  k = nzb_2d(j,i)+1, nzt
+          IF ( ql(k,j,i) > 0.0 )  THEN
+             sed_q(k) = sed_q_const * nc**( -2.0 / 3.0 ) *                    &
+                        ( ql(k,j,i) * hyrho(k) )**( 5.0 / 3.0 )
+          ENDIF
+       ENDDO
+!
+!--   Tendency for q, pt:
+       DO  k = nzb_2d(j,i)+1, nzt
+          tend_q(k,j,i)  = tend_q(k,j,i) + ( sed_q(k+1) - sed_q(k) ) *        &
+                           ddzu(k+1) / hyrho(k) 
+          tend_pt(k,j,i) = tend_pt(k,j,i) - ( sed_q(k+1) - sed_q(k) ) *       &
+                           ddzu(k+1) / hyrho(k) * l_d_cp * pt_d_t(k) 
+       ENDDO
+
+    END SUBROUTINE sedimentation_cloud_ij
+
+    SUBROUTINE sedimentation_rain_ij( i, j )
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE indices
+       USE control_parameters
+       USE statistics
+       
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k, n, n_substep
+       REAL    ::  sed_nr_tend, sed_qr_tend
+
+       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
+!
+!--    Precipitation amount
+       IF ( intermediate_timestep_count == intermediate_timestep_count_max    &
+            .AND.  ( dt_do2d_xy - time_do2d_xy ) <                            &
+            precipitation_amount_interval )  THEN
+
+          precipitation_amount(j,i) = precipitation_amount(j,i) +   &
+                                      prr(nzb_2d(j,i)+1,j,i) *      &
+                                      hyrho(nzb_2d(j,i)+1) * dt_3d
+       ENDIF
+
+    END SUBROUTINE sedimentation_rain_ij
+
+!
+!-- This function computes the gamma function (Press et al., 1992). 
+!-- The gamma function is needed for the calculation of the evaporation 
+!-- of rain drops.
+    FUNCTION gamm( xx ) 
+       
+       USE cloud_parameters
+    
+       IMPLICIT NONE 
+ 
+       REAL    ::  gamm, xx,                                        & 
+                   ser, tmp, x_gamm, y_gamm
+       INTEGER ::  j 
+ 
+       x_gamm = xx 
+       y_gamm = x_gamm 
+       tmp = x_gamm + 5.5 
+       tmp = ( x_gamm + 0.5 ) * LOG( tmp ) - tmp 
+       ser = 1.000000000190015 
+       do j = 1, 6 
+          y_gamm = y_gamm + 1.0 
+          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 
+
+ END MODULE microphysics_mod
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 1047)
+++ palm/trunk/SOURCE/modules.f90	(revision 1053)
@@ -20,5 +20,19 @@
 ! Current revisions:
 ! -----------------
-!
+! necessary expansions according to the two new prognostic equations (nr, qr) 
+! of the two-moment cloud physics scheme:
+! +*_init, flux_l_*, diss_l_*, flux_s_*, diss_s_*, *sws, *swst, tend_*, *, *_p
+! +t*_m, *_1, *_2, *_3, *_av, bc_*_b, bc_*_t, ibc_*_b, ibc_*_t, bc_*_t_val, 
+! +*_vertical_gradient, *_surface_initial_change, *_vertical_gradient_level, 
+! +*_vertical_gradient_level_ind, *_surface, constant_waterflux_*,  
+! +cloud_scheme, icloud_scheme, surface_waterflux_*, sums_ws*s_ws_l, wall_*flux
+!
+! constants for the two-moment scheme:
+! +a_vent, a_term, b_vent, b_term, c_evap, c_term, cof, eps_sb, k_cc, k_cr, k_rr, 
+! +k_br, kappa_rr, kin_vis_air, mu_constant_value, nc, pirho_l, dpirho_l, rho_1,
+! +schmidt, schmidt_p_1d3, stp, x0, xmin, xmax, dt_precipitation, w_precipitation
+!
+! steering parameters for the two_moment scheme:
+! +mu_constant, ventilation_effect
 !
 ! Former revisions:
@@ -369,37 +383,43 @@
           c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l, ddzu, ddzu_pres,     &
           dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor, lad, l_grid,      &
-          ptdf_x, ptdf_y, pt_init, q_init, rdf, rdf_sc, sa_init, ug, u_init,   &
-          u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, w_subs, zu, zw 
+          nr_init, ptdf_x, ptdf_y, pt_init, q_init, qr_init, rdf, rdf_sc,      &
+          sa_init, ug, u_init, u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, &
+          w_subs, zu, zw
 
     REAL, DIMENSION(:,:), ALLOCATABLE ::                                       &
-          c_u, c_v, c_w, diss_s_e, diss_s_pt, diss_s_q, diss_s_sa, diss_s_u,   &
-          diss_s_v, diss_s_w, dzu_mg, dzw_mg, flux_s_e, flux_s_pt, flux_s_q,   &
-          flux_s_sa, flux_s_u, flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg,        &
-          mean_inflow_profiles, pt_slope_ref, qs, qsws, qswst, qswst_remote,   &
-          rif, saswsb, saswst, shf, total_2d_a, total_2d_o, ts, tswst, us,     &
-          usws, uswst, vsws, vswst, z0, z0h
+          c_u, c_v, c_w, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q, diss_s_qr,  &
+          diss_s_sa, diss_s_u, diss_s_v, diss_s_w, dzu_mg, dzw_mg, flux_s_e,   &
+          flux_s_nr, flux_s_pt, flux_s_q, flux_s_qr, flux_s_sa, flux_s_u,      &
+          flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg, mean_inflow_profiles, nrs,  &
+          nrsws, nrswst, pt_slope_ref, qs, qsws, qswst, qswst_remote, qrs,     &
+          qrsws, qrswst, rif, saswsb, saswst, shf, total_2d_a, total_2d_o, ts, &
+          tswst, us, usws, uswst, vsws, vswst, z0, z0h
+          
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE ::                                     &
           canopy_heat_flux, cdc, d, de_dx, de_dy, de_dz, diss, diss_l_e,       &
-          diss_l_pt, diss_l_q, diss_l_sa, diss_l_u, diss_l_v, diss_l_w,        &
-          flux_l_e, flux_l_pt, flux_l_q, flux_l_sa, flux_l_u, flux_l_v,        &
-          flux_l_w, kh, km, lad_s, lad_u, lad_v, lad_w, lai, l_wall, p_loc,    &
-          sec, sls, tend, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n, v_m_r,     &
-          v_m_s, w_m_l, w_m_n, w_m_r, w_m_s
+          diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr, diss_l_sa, diss_l_u,      &
+          diss_l_v, diss_l_w, flux_l_e, flux_l_nr, flux_l_pt, flux_l_q,        &
+          flux_l_qr, flux_l_sa, flux_l_u, flux_l_v, flux_l_w, kh, km, lad_s,   &
+          lad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls, tend, tend_pt,    &
+          tend_nr, tend_q, tend_qr, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n,  &
+          v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s
+            
 
 #if defined( __nopointer )
     REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
-          e, e_p, p, prho, pt, pt_p, q, q_p, ql, ql_c, ql_v, ql_vp, rho, sa,   &
-          sa_p, te_m, tpt_m, tq_m, tsa_m, tu_m, tv_m, tw_m, u, u_p, v, v_p,    &
-          vpt, w, w_p
+          e, e_p, nr, nr_p, p, prho, pt, pt_p, q, q_p, ql, ql_c, ql_v, ql_vp,  &
+          qr, qr_p, rho, sa, sa_p, te_m, tnr_m, tpt_m, tq_m, tqr_m, tsa_m,     &
+          tu_m, tv_m, tw_m, u, u_p, v, v_p, vpt, w, w_p
 #else
     REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
-          e_1, e_2, e_3, p, prho_1, pt_1, pt_2, pt_3, q_1, q_2, q_3, ql_v,     &
-          ql_vp, ql_1, ql_2, rho_1, sa_1, sa_2, sa_3, u_1, u_2, u_3,           &
-          v_1, v_2, v_3, vpt_1, w_1, w_2, w_3
+          e_1, e_2, e_3, p, prho_1, nr_1, nr_2, nr_3, pt_1, pt_2, pt_3, q_1,   &
+          q_2, q_3, ql_v, ql_vp, ql_1, ql_2, qr_1, qr_2, qr_3, rho_1, sa_1,    &
+          sa_2, sa_3, u_1, u_2, u_3, v_1, v_2, v_3, vpt_1, w_1, w_2, w_3
 
     REAL, DIMENSION(:,:,:), POINTER ::                                         &
-          e, e_p, prho, pt, pt_p, q, q_p, ql, ql_c, rho, sa, sa_p, te_m,       &
-          tpt_m, tq_m, tsa_m, tu_m, tv_m, tw_m, u, u_p, v, v_p, vpt, w, w_p
+          e, e_p, nr, nr_p, prho, pt, pt_p, q, q_p, ql, ql_c, qr, qr_p, rho,   &
+          sa, sa_p, te_m, tnr_m, tpt_m, tq_m, tqr_m, tsa_m, tu_m, tv_m, tw_m,  &
+          u, u_p, v, v_p, vpt, w, w_p
 #endif
 
@@ -429,7 +449,8 @@
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
-          e_av, lpt_av, p_av, pc_av, pr_av, pt_av, q_av, ql_av, ql_c_av,     &
-          ql_v_av, ql_vp_av, qv_av, rho_av, s_av, sa_av, u_av, v_av, vpt_av, &
-          w_av
+          e_av, lpt_av, nr_av, p_av, pc_av, pr_av, prr_av, pt_av, q_av, ql_av, &
+          ql_c_av, ql_v_av, ql_vp_av, qr_av, qv_av, rho_av, s_av, sa_av, u_av, &
+          v_av, vpt_av, w_av
+  
  END MODULE averaging
 
@@ -444,21 +465,64 @@
 !------------------------------------------------------------------------------!
 
-    LOGICAL ::  curvature_solution_effects = .FALSE.
-
-    REAL  ::  bfactor, cp = 1005.0, diff_coeff_l = 0.23E-4,                    &
-              effective_coll_efficiency,                                       &
-              eps_ros = 1.0E-4,             &  ! accuracy of Rosenbrock method
-              l_d_cp, l_d_r, l_d_rv, l_v = 2.5E+06,                            &
-              mass_of_solute = 1.0E-17,              & ! soluted NaCl in kg
-              molecular_weight_of_solute = 0.05844,  & ! mol. mass NaCl (kg/mol)
-              molecular_weight_of_water = 0.01801528,& ! mol. mass H2O (kg/mol)
-              prec_time_const = 0.001, ql_crit = 0.0005, rho_l = 1.0E3,        &
-              r_d = 287.0, r_v = 461.51, thermal_conductivity_l = 2.43E-2,     &
-              vanthoff = 2.0                        ! van't Hoff factor (NaCl)
-
-    REAL, DIMENSION(:), ALLOCATABLE   ::  pt_d_t, t_d_pt
-
-    REAL, DIMENSION(:,:), ALLOCATABLE ::  precipitation_amount, &
-                                          precipitation_rate
+    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)
+             cof(6) = (/ 76.18009172947146,      & ! coefficients in the 
+                         -86.50532032941677,     & ! numerical
+                         24.01409824083091,      & ! calculation of the
+                         -1.231739572450155,     & ! gamma function
+                         0.1208650973866179E-2,  &
+                         -0.5395239384953E-5 /), &
+             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)
+             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_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 
+             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 )
+             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)
+             dt_precipitation = 100.0, & ! timestep precipitation (s)
+             w_precipitation = 9.65      ! maximum terminal velocity (m s-1)
+
+    REAL, DIMENSION(:), ALLOCATABLE     ::  dr, hyrho, lambda_r, mu_r, pt_d_t, &
+                                            sed_nr, sed_q, sed_qr, t_d_pt, xr  
+
+    REAL, DIMENSION(:,:), ALLOCATABLE   ::  precipitation_amount, &
+                                            precipitation_rate
+!
+!-- 3D array of precipitation rate
+    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  prr
 
     SAVE
@@ -529,12 +593,15 @@
                              scalar_advec = 'ws-scheme'
     CHARACTER (LEN=20)   ::  bc_e_b = 'neumann', bc_lr = 'cyclic', &
+                             bc_nr_b = 'dirichlet', 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_s_b = 'dirichlet', bc_s_t = 'neumann', &
                              bc_sa_t = 'neumann', &
                              bc_uv_b = 'dirichlet', bc_uv_t = 'dirichlet', &
                              canopy_mode = 'block', &
+                             cloud_scheme = 'seifert_beheng', &
                              coupling_mode = 'uncoupled', &
                              coupling_mode_remote = 'uncoupled', &
@@ -544,5 +611,5 @@
                              random_generator = 'numerical-recipes', &
                              return_addres, return_username, &
-                             timestep_scheme = 'runge-kutta-3'
+                             timestep_scheme = 'runge-kutta-3'                             
     CHARACTER (LEN=40)   ::  avs_data_file, output_format_netcdf, &
                              topography = 'flat'
@@ -577,6 +644,7 @@
                 dp_level_ind_b = 0, dvrp_filecount = 0, &
                 dz_stretch_level_index, gamma_mg, gathered_size, &
-                grid_level, ibc_e_b, ibc_p_b, ibc_p_t, ibc_pt_b, ibc_pt_t, &
-                ibc_q_b, ibc_q_t, ibc_sa_t, ibc_uv_b, ibc_uv_t, &
+                grid_level, ibc_e_b, ibc_nr_b, ibc_nr_t, ibc_p_b, ibc_p_t, &
+                ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_qr_b, ibc_qr_t, &
+                ibc_sa_t, ibc_uv_b, ibc_uv_t, icloud_scheme, &
                 inflow_disturbance_begin = -1, inflow_disturbance_end = -1, &
                 intermediate_timestep_count, intermediate_timestep_count_max, &
@@ -600,6 +668,8 @@
                 mask_size(max_masks,3) = -1, mask_size_l(max_masks,3) = -1, &
                 mask_start_l(max_masks,3) = -1, &
+                nr_vertical_gradient_level_ind(10) = -9999, &
                 pt_vertical_gradient_level_ind(10) = -9999, &
                 q_vertical_gradient_level_ind(10) = -9999, &
+                qr_vertical_gradient_level_ind(10) = -9999, &
                 sa_vertical_gradient_level_ind(10) = -9999, &
                 section(100,3), section_xy(100) = -9999, &
@@ -631,10 +701,12 @@
                 constant_top_momentumflux = .FALSE., &
                 constant_top_salinityflux = .TRUE., &
-                constant_waterflux = .TRUE., create_disturbances = .TRUE., &
+                constant_waterflux = .TRUE., constant_waterflux_nr = .TRUE., &
+                constant_waterflux_qr = .TRUE., create_disturbances = .TRUE., &
                 data_output_2d_on_each_pe = .TRUE., & 
                 dissipation_control = .FALSE., disturbance_created = .FALSE., &
                 do2d_at_begin = .FALSE., do3d_at_begin = .FALSE., &
                 do3d_compress = .FALSE., do_sum = .FALSE., &
-                dp_external = .FALSE., dp_smooth = .FALSE., dt_fixed = .FALSE., &
+                dp_external = .FALSE., dp_smooth = .FALSE., &
+                drizzle = .TRUE., dt_fixed = .FALSE., &
                 dt_3d_reached, dt_3d_reached_l, exchange_mg = .FALSE., &
                 first_call_lpm = .TRUE., &
@@ -668,5 +740,5 @@
              averaging_interval = 0.0, averaging_interval_pr = 9999999.9, &
              averaging_interval_sp = 9999999.9, bc_pt_t_val, bc_q_t_val, &
-             bottom_salinityflux = 0.0, &
+             bc_qr_t_val, bc_nr_t_val, bottom_salinityflux = 0.0, &
              building_height = 50.0, building_length_x = 50.0, &
              building_length_y = 50.0, building_wall_left = 9999999.9, &
@@ -699,4 +771,5 @@
              maximum_cpu_time_allowed = 0.0,  &
              molecular_viscosity = 1.461E-5, &
+             nr_surface = 0.0, nr_surface_initial_change = 0.0, &
              old_dt = 1.0E-10, omega = 7.29212E-5, omega_sor = 1.8, &
              particle_maximum_age = 9999999.9, &
@@ -707,4 +780,5 @@
              pt_surface = 300.0, pt_surface_initial_change = 0.0, &
              q_surface = 0.0, q_surface_initial_change = 0.0, &
+             qr_surface = 0.0, qr_surface_initial_change = 0.0, &
              rayleigh_damping_factor = -1.0, rayleigh_damping_height = -1.0, &
              recycling_width = 9999999.9, residual_limit = 1.0E-4, &
@@ -719,4 +793,5 @@
              surface_heatflux = 9999999.9, surface_pressure = 1013.25, &
              surface_scalarflux = 9999999.9, surface_waterflux = 9999999.9, &
+             surface_waterflux_nr = 0.0, surface_waterflux_qr = 0.0, &
              s_surface = 0.0, s_surface_initial_change = 0.0, &
              termination_time_needed = -1.0, time_coupling = 0.0, &
@@ -740,8 +815,12 @@
              lad_vertical_gradient_level(10) = -9999999.9, &
              mask_scale(3), &
+             nr_vertical_gradient(10) = 0.0, &
+             nr_vertical_gradient_level(10) = -1.0, &
              pt_vertical_gradient(10) = 0.0, &
              pt_vertical_gradient_level(10) = -9999999.9, &
              q_vertical_gradient(10) = 0.0, &
              q_vertical_gradient_level(10) = -1.0, &
+             qr_vertical_gradient(10) = 0.0, &
+             qr_vertical_gradient_level(10) = -1.0, &
              s_vertical_gradient(10) = 0.0, &
              s_vertical_gradient_level(10) = -1.0, &
@@ -759,5 +838,6 @@
              volume_flow(1:2) = 0.0, volume_flow_area(1:2) = 0.0, &
              volume_flow_initial(1:2) = 0.0, wall_heatflux(0:4) = 0.0, &
-             wall_humidityflux(0:4) = 0.0, wall_qflux(0:4) = 0.0, &
+             wall_humidityflux(0:4) = 0.0, wall_nrflux(0:4) = 0.0, & 
+             wall_qflux(0:4) = 0.0, wall_qrflux(0:4) = 0.0, &
              wall_salinityflux(0:4) = 0.0, wall_scalarflux(0:4) = 0.0, &
              subs_vertical_gradient(10) = 0.0, &
@@ -1532,6 +1612,10 @@
                                               sums_us2_ws_l, sums_vs2_ws_l,   &
                                               sums_ws2_ws_l,                  &
+                                              sums_wsnrs_ws_l,                &
                                               sums_wspts_ws_l,                &
-                                              sums_wssas_ws_l,sums_wsqs_ws_l
+                                              sums_wssas_ws_l,                &
+                                              sums_wsqs_ws_l,                 &
+                                              sums_wsqrs_ws_l
+                                              
     REAL, DIMENSION(:,:,:), ALLOCATABLE   ::  hom_sum, rmask, spectrum_x, &
                                               spectrum_y, sums_l, sums_l_l, &
Index: palm/trunk/SOURCE/netcdf.f90
===================================================================
--- palm/trunk/SOURCE/netcdf.f90	(revision 1047)
+++ palm/trunk/SOURCE/netcdf.f90	(revision 1053)
@@ -23,5 +23,5 @@
 ! Current revisions:
 ! ------------------
-!
+! +qr, nr, prr
 !
 ! Former revisions:
@@ -466,6 +466,7 @@
 !
 !--             Most variables are defined on the scalar grid
-                CASE ( 'e', 'lpt', 'p', 'pc', 'pr', 'pt', 'q', 'ql', 'ql_c', &
-                       'ql_v', 'ql_vp', 'qv', 'rho', 's', 'sa', 'vpt' )
+                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',&
+                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',    &
+                       'rho', 's', 'sa', 'vpt' )
 
                    grid_x = 'x'
@@ -994,6 +995,7 @@
 !
 !--             Most variables are defined on the scalar grid
-                CASE ( 'e', 'lpt', 'p', 'pc', 'pr', 'pt', 'q', 'ql', 'ql_c', &
-                       'ql_v', 'ql_vp', 'qv', 'rho', 's', 'sa', 'vpt' )
+                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',   &
+                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho',&
+                       's', 'sa', 'vpt' )
 
                    grid_x = 'x'
@@ -1556,7 +1558,8 @@
 !
 !--                   Most variables are defined on the zu grid
-                      CASE ( 'e_xy', 'lpt_xy', 'p_xy', 'pc_xy', 'pr_xy', 'pt_xy',&
-                             'q_xy', 'ql_xy', 'ql_c_xy', 'ql_v_xy', 'ql_vp_xy',  &
-                             'qv_xy', 'rho_xy', 's_xy', 'sa_xy', 'vpt_xy' )
+                      CASE ( 'e_xy', 'lpt_xy', 'nr_xy', 'p_xy', 'pc_xy', 'pr_xy',&
+                             'prr_xy', 'pt_xy', 'q_xy', 'qc_xy', 'ql_xy',        &
+                             'ql_c_xy', 'ql_v_xy', 'ql_vp_xy', 'qr_xy', 'qv_xy', &
+                             'rho_xy', 's_xy', 'sa_xy', 'vpt_xy' )
 
                          grid_x = 'x'
@@ -2142,7 +2145,8 @@
 !
 !--                Most variables are defined on the zu grid
-                   CASE ( 'e_xz', 'lpt_xz', 'p_xz', 'pc_xz', 'pr_xz', 'pt_xz', &
-                          'q_xz', 'ql_xz', 'ql_c_xz', 'ql_v_xz', 'ql_vp_xz',   &
-                          'qv_xz', 'rho_xz', 's_xz', 'sa_xz', 'vpt_xz' )
+                   CASE ( 'e_xz', 'lpt_xz', 'nr_xz', 'p_xz', 'pc_xz', 'pr_xz', &
+                          'prr_xz', 'pt_xz', 'q_xz', 'qc_xz', 'ql_xz',         &
+                          'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz', 'qv_xz',  &
+                          'rho_xz', 's_xz', 'sa_xz', 'vpt_xz' )
 
                       grid_x = 'x'
@@ -2728,7 +2732,8 @@
 !
 !--                Most variables are defined on the zu grid
-                   CASE ( 'e_yz', 'lpt_yz', 'p_yz', 'pc_yz', 'pr_yz', 'pt_yz',&
-                          'q_yz', 'ql_yz', 'ql_c_yz', 'ql_v_yz', 'ql_vp_yz',  &
-                          'qv_yz', 'rho_yz', 's_yz', 'sa_yz', 'vpt_yz' )
+                   CASE ( 'e_yz', 'lpt_yz', 'nr_yz', 'p_yz', 'pc_yz', 'pr_yz',&
+                          'prr_yz', 'pt_yz', 'q_yz', 'qc_yz', 'ql_yz',        &
+                          'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qr_yz', 'qv_yz', &
+                          'rho_yz', 's_yz', 'sa_yz', 'vpt_yz' )
 
                       grid_x = 'x'
@@ -3915,6 +3920,7 @@
 !
 !--             Most variables are defined on the zu levels
-                CASE ( 'e', 'lpt', 'p', 'pc', 'pr', 'pt', 'q', 'ql', 'ql_c',  &
-                       'ql_v', 'ql_vp', 'qv', 'rho', 's', 'sa', 'u', 'v', 'vpt' )
+                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',   &
+                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho',&
+                       's', 'sa', 'u', 'v', 'vpt' )
 
                    grid_z = 'zu'
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 1047)
+++ palm/trunk/SOURCE/parin.f90	(revision 1053)
@@ -20,5 +20,9 @@
 ! Current revisions:
 ! -----------------
-!
+! necessary expansions according to the two new prognostic equations (nr, qr) 
+! of the two-moment cloud physics scheme and steering parameters:
+! +*_init, *_surface, *_surface_initial_change, *_vertical_gradient, 
+! +*_vertical_gradient_level, surface_waterflux_*, 
+! +cloud_scheme, drizzle, mu_constant, mu_constant_value, ventilation_effect
 !
 ! Former revisions:
@@ -194,9 +198,10 @@
              canyon_width_x, canyon_width_y, canyon_wall_left, &
              canyon_wall_south, cfl_factor, cloud_droplets, cloud_physics, &
-             collective_wait, conserve_volume_flow, conserve_volume_flow_mode, &
-             coupling_start_time, cthf, curvature_solution_effects, &
-             cycle_mg, damp_level_1d, dissipation_1d, & !dissipation_control, &
+             cloud_scheme, collective_wait, conserve_volume_flow, &
+             conserve_volume_flow_mode, coupling_start_time, cthf, &
+             curvature_solution_effects, cycle_mg, damp_level_1d, &
+             dissipation_1d, & !dissipation_control, &
              dp_external, dp_level_b, dp_smooth, dpdxy, drag_coefficient, &
-             dt, dt_pr_1d, dt_run_control_1d, dx, dy, dz, dz_max, & 
+             drizzle, dt, dt_pr_1d, dt_run_control_1d, dx, dy, dz, dz_max, & 
              dz_stretch_factor, dz_stretch_level, e_init, e_min, end_time_1d, &
              fft_method, galilei_transformation, humidity, &
@@ -208,5 +213,8 @@
              loop_optimization, masking_method, mg_cycles, &
              mg_switch_to_pe0_level, mixing_length_1d, momentum_advec, &
-             netcdf_precision, neutral, ngsrb, nsor, &
+             mu_constant, mu_constant_value, &
+             netcdf_precision, neutral, ngsrb, nr_surface, &
+             nr_surface_initial_change, nr_vertical_gradient, &
+             nr_vertical_gradient_level, nsor, &
              nsor_ini, nx, ny, nz, ocean, omega, omega_sor, &
              passive_scalar, pch_index, phi, plant_canopy, prandtl_layer, &
@@ -215,5 +223,7 @@
              pt_surface_initial_change, pt_vertical_gradient, &
              pt_vertical_gradient_level, q_surface, q_surface_initial_change, &
-             q_vertical_gradient, q_vertical_gradient_level, radiation, &
+             q_vertical_gradient, q_vertical_gradient_level, &
+             qr_surface, qr_surface_initial_change, &
+             qr_vertical_gradient, qr_vertical_gradient_level, radiation, &
              random_generator, random_heatflux, rayleigh_damping_factor, &
              rayleigh_damping_height, recycling_width, residual_limit, &
@@ -223,5 +233,6 @@
              statistic_regions, subs_vertical_gradient, &
              subs_vertical_gradient_level, surface_heatflux, surface_pressure, &
-             surface_scalarflux, surface_waterflux, s_surface, &
+             surface_scalarflux, surface_waterflux, surface_waterflux_nr, &
+             surface_waterflux_qr, s_surface, &
              s_surface_initial_change, s_vertical_gradient, &
              s_vertical_gradient_level, timestep_scheme, &
@@ -232,8 +243,8 @@
              use_top_fluxes, use_ug_for_galilei_tr, use_upstream_for_tke, &
              uv_heights, u_bulk, u_profile, vg_surface, vg_vertical_gradient, &
-             vg_vertical_gradient_level, v_bulk, v_profile, wall_adjustment, &
-             wall_heatflux, wall_humidityflux, wall_scalarflux, z0h_factor
-
-
+             vg_vertical_gradient_level, v_bulk, v_profile, ventilation_effect, &
+             wall_adjustment, wall_heatflux, wall_humidityflux, wall_scalarflux, &
+             z0h_factor
+      
     NAMELIST /d3par/  averaging_interval, averaging_interval_pr, &
              create_disturbances, &
@@ -391,4 +402,5 @@
 !--                  also check the allocate statement in routine read_var_list
              ALLOCATE( lad(0:nz+1),pt_init(0:nz+1), q_init(0:nz+1),           &
+                       qr_init(0:nz+1), nr_init(0:nz+1),                      &
                        sa_init(0:nz+1), ug(0:nz+1), u_init(0:nz+1),           &
                        v_init(0:nz+1), vg(0:nz+1),                            &
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 1047)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 1053)
@@ -20,5 +20,8 @@
 ! Current revisions:
 ! -----------------
-!
+! implementation of two new prognostic equations for rain drop concentration (nr)
+! and rain water content (qr)
+!
+! currently, only available for cache loop optimization
 !
 ! Former revisions:
@@ -154,4 +157,5 @@
     USE diffusion_w_mod
     USE impact_of_latent_heat_mod
+    USE microphysics_mod
     USE plant_canopy_model_mod
     USE production_e_mod
@@ -409,4 +413,31 @@
 
 !
+!--      If required, calculate tendencies for total water content, rain water 
+!--      content, rain drop concentration and liquid temperature
+
+         IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+
+            tend_q(:,j,i)  = 0.0
+            tend_qr(:,j,i) = 0.0
+            tend_nr(:,j,i) = 0.0
+            tend_pt(:,j,i) = 0.0
+!
+!--         Droplet size distribution (dsd) properties are needed for the 
+!--         computation of selfcollection, breakup, evaporation and 
+!--         sedimentation of rain.  
+            IF ( precipitation )  THEN
+               CALL dsd_properties( i,j ) 
+               CALL autoconversion( i,j )
+               CALL accretion( i,j )
+               CALL selfcollection_breakup( i,j )
+               CALL evaporation_rain( i,j )
+               CALL sedimentation_rain( i,j )
+            ENDIF
+
+            IF ( drizzle )  CALL sedimentation_cloud( i,j )
+
+         ENDIF
+
+!
 !--       If required, compute prognostic equation for potential temperature
           IF ( .NOT. neutral )  THEN
@@ -433,8 +464,11 @@
              ENDIF
 
-!
-!--          If required compute impact of latent heat due to precipitation
-             IF ( precipitation )  THEN
-                CALL impact_of_latent_heat( i, j )
+!--          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
              ENDIF
 
@@ -555,16 +589,19 @@
              ENDIF
              CALL diffusion_s( i, j, q, qsws, qswst, wall_qflux )
-       
-!
-!--          If required compute decrease of total water content due to
-!--          precipitation
-             IF ( precipitation )  THEN
-                CALL calc_precipitation( i, j )
-             ENDIF
-
+     
+!
+!--          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
+             ENDIF
 !
 !--          Sink or source of scalar concentration due to canopy elements
              IF ( plant_canopy ) CALL plant_canopy_model( i, j, 5 )
 
+!
 !--          If required compute influence of large-scale subsidence/ascent
              IF ( large_scale_subsidence )  THEN
@@ -598,4 +635,119 @@
                    ENDDO
                 ENDIF
+             ENDIF
+
+!
+!--          If required, calculate prognostic equations for rain water content 
+!--          and rain drop concentration
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+!
+!--             Calculate prognostic equation for rain water content
+                tend(:,j,i) = 0.0
+                IF ( timestep_scheme(1:5) == 'runge' ) &
+                THEN
+                   IF ( ws_scheme_sca )  THEN
+                      CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr,       & 
+                                       diss_s_qr, flux_l_qr, diss_l_qr, &
+                                       i_omp_start, tn )
+                   ELSE 
+                      CALL advec_s_pw( i, j, qr )
+                   ENDIF
+                ELSE
+                   CALL advec_s_up( i, j, qr )
+                ENDIF
+                CALL diffusion_s( i, j, qr, qrsws, qrswst, wall_qrflux )
+
+!
+!--             Using microphysical tendencies (autoconversion, accretion, 
+!--             evaporation; if required: sedimentation)
+                tend(:,j,i) = tend(:,j,i) + tend_qr(:,j,i)
+
+!
+!--             If required, compute influence of large-scale subsidence/ascent
+                IF ( large_scale_subsidence )  THEN
+                   CALL subsidence( i, j, tend, qr, qr_init )
+                ENDIF
+
+!              CALL user_actions( i, j, 'qr-tendency' )
+
+!
+!--             Prognostic equation for rain water content
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +    &
+                                                     tsc(3) * tqr_m(k,j,i) )     &
+                                         - tsc(5) * rdf_sc(k) *                  &
+                                        ( qr(k,j,i) - qr_init(k) )
+                   IF ( qr_p(k,j,i) < 0.0 )  qr_p(k,j,i) = 0.1 * qr(k,j,i)
+                ENDDO
+!
+!--             Calculate tendencies for the next Runge-Kutta step
+                IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                   IF ( intermediate_timestep_count == 1 )  THEN
+                      DO  k = nzb_s_inner(j,i)+1, nzt
+                         tqr_m(k,j,i) = tend(k,j,i)
+                      ENDDO
+                   ELSEIF ( intermediate_timestep_count < &
+                            intermediate_timestep_count_max )  THEN
+                      DO  k = nzb_s_inner(j,i)+1, nzt
+                         tqr_m(k,j,i) = -9.5625 * tend(k,j,i) + &
+                                        5.3125 * tqr_m(k,j,i)
+                      ENDDO
+                   ENDIF
+                ENDIF
+
+!
+!--             Calculate prognostic equation for rain drop concentration.
+                tend(:,j,i) = 0.0
+                IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                   IF ( ws_scheme_sca )  THEN
+                      CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr,    & 
+                                    diss_s_nr, flux_l_nr, diss_l_nr, &
+                                    i_omp_start, tn )
+                   ELSE 
+                      CALL advec_s_pw( i, j, nr )
+                   ENDIF
+                ELSE
+                   CALL advec_s_up( i, j, nr )
+                ENDIF
+                CALL diffusion_s( i, j, nr, 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)
+
+!
+!--             If required, compute influence of large-scale subsidence/ascent
+                IF ( large_scale_subsidence )  THEN
+                   CALL subsidence( i, j, tend, nr, nr_init )
+                ENDIF
+
+!                CALL user_actions( i, j, 'nr-tendency' )
+
+!
+!--             Prognostic equation for rain drop concentration
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
+                                                  tsc(3) * tnr_m(k,j,i) )     &
+                                      - tsc(5) * rdf_sc(k) *                  &
+                                        ( nr(k,j,i) - nr_init(k) )
+                   IF ( nr_p(k,j,i) < 0.0 )  nr_p(k,j,i) = 0.1 * nr(k,j,i)
+                ENDDO
+!
+!--             Calculate tendencies for the next Runge-Kutta step
+                IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                   IF ( intermediate_timestep_count == 1 )  THEN
+                      DO  k = nzb_s_inner(j,i)+1, nzt
+                         tnr_m(k,j,i) = tend(k,j,i)
+                      ENDDO
+                   ELSEIF ( intermediate_timestep_count < &
+                            intermediate_timestep_count_max )  THEN
+                      DO  k = nzb_s_inner(j,i)+1, nzt
+                         tnr_m(k,j,i) = -9.5625 * tend(k,j,i) + &
+                                        5.3125 * tnr_m(k,j,i)
+                      ENDDO
+                   ENDIF
+                ENDIF
+
              ENDIF
 
Index: palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/read_3d_binary.f90	(revision 1047)
+++ palm/trunk/SOURCE/read_3d_binary.f90	(revision 1053)
@@ -20,5 +20,7 @@
 ! Current revisions:
 ! -----------------
-!
+! necessary expansions according to the two new prognostic equations (nr, qr) 
+! of the two-moment cloud physics scheme:
+! +prr, prr_av, *, *_av, *s, *sws, *swst
 !
 ! Former revisions:
@@ -279,5 +281,5 @@
 !--    First compare the version numbers
        READ ( 13 )  version_on_file
-       binary_version = '3.2'
+       binary_version = '3.3'
        IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
           WRITE( message_string, * ) 'version mismatch concerning data ',      &
@@ -457,4 +459,20 @@
                                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
+                CASE ( 'prr' )
+                   IF ( .NOT. ALLOCATED( prr ) )  THEN
+                      ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   prr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'prr_av' )
+                   IF ( .NOT. ALLOCATED( prr_av ) )  THEN
+                      ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   prr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
                 CASE ( 'precipitation_amount' )
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
@@ -508,4 +526,60 @@
                    ql_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
                                        tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'qr' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   qr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'nr' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   nr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'qr_av' )
+                   IF ( .NOT. ALLOCATED( qr_av ) )  THEN
+                      ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   qr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'nr_av' )
+                   IF ( .NOT. ALLOCATED( nr_av ) )  THEN
+                      ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   nr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'qrs' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   qrs(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'qrsws' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   qrsws(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'qrswst' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   qrswst(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'nrs' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   nrs(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'nrsws' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   nrsws(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'nrswst' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   nrswst(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 'ql_c_av' )
Index: palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- palm/trunk/SOURCE/read_var_list.f90	(revision 1047)
+++ palm/trunk/SOURCE/read_var_list.f90	(revision 1053)
@@ -20,5 +20,13 @@
 ! Current revisions:
 ! ------------------
-!
+! necessary expansions according to the two new prognostic equations (nr, qr) 
+! of the two-moment cloud physics scheme:
+! +bc_*_b, +bc_*_t, +bc_*_t_val, *_init, *_surface, *_surface_initial_change,
+! +*_vertical_gradient, +*_vertical_gradient_level, *_vertical_gradient_level_ind,
+! +surface_waterflux_*
+!
+! in addition, steering parameters parameters of the two-moment cloud physics 
+! scheme:   
+! +cloud_scheme, +drizzle, +mu_constant, +mu_constant_value, +ventilation_effect
 !
 ! Former revisions:
@@ -169,5 +177,5 @@
 !-- Make version number check first
     READ ( 13 )  version_on_file
-    binary_version = '3.5'
+    binary_version = '3.6'
     IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
        WRITE( message_string, * ) 'version mismatch concerning control ', &
@@ -235,5 +243,5 @@
        ALLOCATE( lad(0:nz+1), ug(0:nz+1), u_init(0:nz+1), vg(0:nz+1),    &
                  v_init(0:nz+1), pt_init(0:nz+1), q_init(0:nz+1),        &
-                 sa_init(0:nz+1),                                        &
+                 qr_init(0:nz+1), nr_init(0:nz+1), sa_init(0:nz+1),      &
                  hom(0:nz+1,2,pr_palm+max_pr_user,0:statistic_regions),  &
                  hom_sum(0:nz+1,pr_palm+max_pr_user,0:statistic_regions) )
@@ -268,4 +276,10 @@
           CASE ( 'bc_ns' )
              READ ( 13 )  bc_ns
+          CASE ( 'bc_nr_b' )
+             READ ( 13 )  bc_nr_b
+          CASE ( 'bc_nr_t' )
+             READ ( 13 )  bc_nr_t
+          CASE ( 'bc_nr_t_val' )
+             READ ( 13 )  bc_nr_t_val
           CASE ( 'bc_p_b' )
              READ ( 13 )  bc_p_b
@@ -284,4 +298,10 @@
           CASE ( 'bc_q_t_val' )
              READ ( 13 )  bc_q_t_val
+          CASE ( 'bc_qr_b' )
+             READ ( 13 )  bc_qr_b
+          CASE ( 'bc_qr_t' )
+             READ ( 13 )  bc_qr_t
+          CASE ( 'bc_qr_t_val' )
+             READ ( 13 )  bc_qr_t_val
           CASE ( 'bc_s_b' )
              READ ( 13 )  bc_s_b
@@ -326,4 +346,6 @@
           CASE ( 'cloud_physics' )
              READ ( 13 )  cloud_physics
+          CASE ( 'cloud_scheme' )
+             READ ( 13 )  cloud_scheme
           CASE ( 'collective_wait' )
              READ ( 13 )  collective_wait
@@ -356,4 +378,6 @@
           CASE ( 'drag_coefficient' )
              READ ( 13 )  drag_coefficient
+          CASE ( 'drizzle' )
+             READ ( 13 )  drizzle
           CASE ( 'dt_pr_1d' )
              READ ( 13 )  dt_pr_1d
@@ -438,4 +462,8 @@
           CASE ( 'momentum_advec' )
              READ ( 13 )  momentum_advec
+          CASE ( 'mu_constant' )
+             READ ( 13 )  mu_constant
+          CASE ( 'mu_constant_value' )
+             READ ( 13 )  mu_constant_value
           CASE ( 'netcdf_precision' )
              READ ( 13 )  netcdf_precision
@@ -444,4 +472,16 @@
           CASE ( 'ngsrb' )
              READ ( 13 )  ngsrb
+          CASE ( 'nr_init' )
+             READ ( 13 )  nr_init
+          CASE ( 'nr_surface' )
+             READ ( 13 )  nr_surface
+          CASE ( 'nr_surface_initial_change' )
+             READ ( 13 )  nr_surface_initial_change
+          CASE ( 'nr_vertical_gradient' )
+             READ ( 13 )  nr_vertical_gradient
+          CASE ( 'nr_vertical_gradient_level' )
+             READ ( 13 )  nr_vertical_gradient_level
+          CASE ( 'nr_vertical_gradient_level_ind' )
+             READ ( 13 )  nr_vertical_gradient_level_ind
           CASE ( 'nsor' )
              READ ( 13 )  nsor
@@ -510,4 +550,16 @@
           CASE ( 'q_vertical_gradient_level_ind' )
              READ ( 13 )  q_vertical_gradient_level_ind
+          CASE ( 'qr_init' )
+             READ ( 13 )  qr_init
+          CASE ( 'qr_surface' )
+             READ ( 13 )  qr_surface
+          CASE ( 'qr_surface_initial_change' )
+             READ ( 13 )  qr_surface_initial_change
+          CASE ( 'qr_vertical_gradient' )
+             READ ( 13 )  qr_vertical_gradient
+          CASE ( 'qr_vertical_gradient_level' )
+             READ ( 13 )  qr_vertical_gradient_level
+          CASE ( 'qr_vertical_gradient_level_ind' )
+             READ ( 13 )  qr_vertical_gradient_level_ind
           CASE ( 'radiation' )
              READ ( 13 )  radiation
@@ -555,5 +607,9 @@
              READ ( 13 )  surface_scalarflux             
           CASE ( 'surface_waterflux' )
-             READ ( 13 )  surface_waterflux             
+             READ ( 13 )  surface_waterflux      
+          CASE ( 'surface_waterflux_nr' )
+             READ ( 13 )  surface_waterflux_nr    
+          CASE ( 'surface_waterflux_qr' )
+             READ ( 13 )  surface_waterflux_qr    
           CASE ( 's_surface' )
              READ ( 13 )  s_surface
@@ -656,4 +712,6 @@
           CASE ( 'v_max_ijk' )
              READ ( 13 )  v_max_ijk
+          CASE ( 'ventilation_effect' )
+             READ ( 13 )  ventilation_effect
           CASE ( 'vg' )
              READ ( 13 )  vg
Index: palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 1047)
+++ palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 1053)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! +nr, prr, qr
 !
 ! Former revisions:
@@ -115,4 +115,10 @@
                 lwp_av = 0.0
 
+             CASE ( 'nr' )
+                IF ( .NOT. ALLOCATED( nr_av ) )  THEN
+                   ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                nr_av = 0.0
+
              CASE ( 'p' )
                 IF ( .NOT. ALLOCATED( p_av ) )  THEN
@@ -133,4 +139,10 @@
                 pr_av = 0.0
 
+             CASE ( 'prr' )
+                IF ( .NOT. ALLOCATED( prr_av ) )  THEN
+                   ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                prr_av = 0.0
+
              CASE ( 'prr*' )
                 IF ( .NOT. ALLOCATED( precipitation_rate_av ) )  THEN
@@ -175,4 +187,10 @@
                 ql_vp_av = 0.0
 
+             CASE ( 'qr' )
+                IF ( .NOT. ALLOCATED( qr_av ) )  THEN
+                   ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                qr_av = 0.0
+
              CASE ( 'qsws*' )
                 IF ( .NOT. ALLOCATED( qsws_av ) )  THEN
@@ -301,4 +319,13 @@
                    lwp_av(j,i) = lwp_av(j,i) + SUM( ql(nzb:nzt,j,i) * &
                                                     dzw(1:nzt+1) )
+                ENDDO
+             ENDDO
+
+          CASE ( 'nr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      nr_av(k,j,i) = nr_av(k,j,i) + nr(k,j,i)
+                   ENDDO
                 ENDDO
              ENDDO
@@ -423,4 +450,13 @@
              ENDDO
 
+          CASE ( 'qr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      qr_av(k,j,i) = qr_av(k,j,i) + qr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
           CASE ( 'qsws*' )
              DO  i = nxlg, nxrg
Index: palm/trunk/SOURCE/swap_timelevel.f90
===================================================================
--- palm/trunk/SOURCE/swap_timelevel.f90	(revision 1047)
+++ palm/trunk/SOURCE/swap_timelevel.f90	(revision 1053)
@@ -20,6 +20,6 @@
 ! Current revisions:
 ! -----------------
+! swap of timelevels for nr, qr added
 ! 
-!
 ! Former revisions:
 ! -----------------
@@ -91,5 +91,9 @@
     ENDIF
     IF ( humidity  .OR.  passive_scalar )  THEN
-       q = q_p
+       q = q_p             
+       IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
+          qr = qr_p
+          nr = nr_p
+       ENDIF
     ENDIF
 
@@ -116,4 +120,8 @@
           IF ( humidity  .OR.  passive_scalar )  THEN
              q => q_1;    q_p => q_2
+             IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
+                qr => qr_1;    qr_p => qr_2
+                nr => nr_1;    nr_p => nr_2
+             ENDIF
           ENDIF
 
@@ -135,4 +143,8 @@
           IF ( humidity  .OR.  passive_scalar )  THEN
              q => q_2;    q_p => q_1
+             IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
+                qr => qr_2;    qr_p => qr_1
+                nr => nr_2;    nr_p => nr_1
+             ENDIF
           ENDIF
 
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 1047)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 1053)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! exchange of ghost points for nr, qr added
 !
 ! Former revisions:
@@ -240,5 +240,11 @@
              CALL exchange_horiz( prho, nbgp )
           ENDIF
-          IF (humidity  .OR.  passive_scalar)  CALL exchange_horiz( q_p, nbgp )
+          IF (humidity  .OR.  passive_scalar)  THEN
+             CALL exchange_horiz( q_p, nbgp )
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                CALL exchange_horiz( qr_p, nbgp )
+                CALL exchange_horiz( nr_p, nbgp )
+             ENDIF
+          ENDIF
           IF ( cloud_droplets )  THEN
              CALL exchange_horiz( ql, nbgp )
Index: palm/trunk/SOURCE/timestep.f90
===================================================================
--- palm/trunk/SOURCE/timestep.f90	(revision 1047)
+++ palm/trunk/SOURCE/timestep.f90	(revision 1053)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-!
+! timestep is reduced in two-moment cloud scheme according to the maximum 
+! terminal velocity of rain drops
 !
 ! Former revisions:
@@ -75,4 +76,5 @@
 
     USE arrays_3d
+    USE cloud_parameters
     USE control_parameters
     USE cpulog
@@ -301,5 +303,6 @@
 !--    step minus a reduction (cfl_factor) to be on the safe side.
 !--    The time step must not exceed the maximum allowed value.
-       dt_3d = cfl_factor * MIN( dt_diff, dt_plant_canopy, dt_u, dt_v, dt_w )
+       dt_3d = cfl_factor * MIN( dt_diff, dt_plant_canopy, dt_u, dt_v, dt_w,   &
+                                 dt_precipitation )
        dt_3d = MIN( dt_3d, dt_max )
 
Index: palm/trunk/SOURCE/user_actions.f90
===================================================================
--- palm/trunk/SOURCE/user_actions.f90	(revision 1047)
+++ palm/trunk/SOURCE/user_actions.f90	(revision 1053)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! -----------------
+! +qr-tendency, nr-tendency
 !
 ! Former revisions:
@@ -136,4 +137,10 @@
 
 
+          CASE ( 'qr-tendency' )
+
+
+          CASE ( 'nr-tendency' )
+
+
           CASE DEFAULT
              message_string = 'unknown location "' // location // '"'
@@ -188,4 +195,10 @@
 
           CASE ( 'q-tendency' )
+
+
+          CASE ( 'qr-tendency' )
+
+
+          CASE ( 'nr-tendency' )
 
 
Index: palm/trunk/SOURCE/write_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/write_3d_binary.f90	(revision 1047)
+++ palm/trunk/SOURCE/write_3d_binary.f90	(revision 1053)
@@ -20,5 +20,7 @@
 ! Current revisions:
 ! -----------------
-!
+! necessary expansions according to the two new prognostic equations (nr, qr) 
+! of the two-moment cloud physics scheme:
+! +prr, prr_av, *, *_av, *s, *sws, *swst
 !
 ! Former revisions:
@@ -110,5 +112,5 @@
 !
 !-- Write arrays.
-    binary_version = '3.2'
+    binary_version = '3.3'
 
     WRITE ( 14 )  binary_version
@@ -145,4 +147,10 @@
        WRITE ( 14 )  'pr_av               ';  WRITE ( 14 )  pr_av
     ENDIF
+    IF ( ALLOCATED( prr ) )  THEN
+       WRITE ( 14 )  'prr                 ';  WRITE ( 14 )  prr
+    ENDIF
+    IF ( ALLOCATED( prr_av ) )  THEN
+       WRITE ( 14 )  'prr_av              ';  WRITE ( 14 )  prr_av
+    ENDIF
     IF ( ALLOCATED( precipitation_amount ) )  THEN
        WRITE ( 14 )  'precipitation_amount';  WRITE ( 14 )  precipitation_amount
@@ -166,4 +174,20 @@
              WRITE ( 14 )  'ql_av               ';  WRITE ( 14 )  ql_av
           ENDIF
+          IF ( icloud_scheme == 0 )  THEN
+             WRITE ( 14 )  'qr                  ';  WRITE ( 14 ) qr
+             WRITE ( 14 )  'nr                  ';  WRITE ( 14 ) nr
+             IF ( ALLOCATED( qr_av ) )  THEN
+                WRITE ( 14 )  'qr_av               ';  WRITE ( 14 )  qr_av
+             ENDIF
+             IF ( ALLOCATED( nr_av ) )  THEN
+                WRITE ( 14 )  'nr_av               ';  WRITE ( 14 )  nr_av
+             ENDIF
+             WRITE ( 14 )  'qrs                 ';  WRITE ( 14 ) qrs
+             WRITE ( 14 )  'qrsws               ';  WRITE ( 14 ) qrsws
+             WRITE ( 14 )  'qrswst              ';  WRITE ( 14 ) qrswst
+             WRITE ( 14 )  'nrs                 ';  WRITE ( 14 ) nrs
+             WRITE ( 14 )  'nrsws               ';  WRITE ( 14 ) nrsws
+             WRITE ( 14 )  'nrswst              ';  WRITE ( 14 ) nrswst
+          ENDIF
        ENDIF
        WRITE ( 14 )  'qs                  ';  WRITE ( 14 ) qs
Index: palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- palm/trunk/SOURCE/write_var_list.f90	(revision 1047)
+++ palm/trunk/SOURCE/write_var_list.f90	(revision 1053)
@@ -20,5 +20,13 @@
 ! Current revisions:
 ! -----------------
-!
+! necessary expansions according to the two new prognostic equations (nr, qr) 
+! of the two-moment cloud physics scheme:
+! +bc_*_b, bc_*_t, bc_*_t_val, *_init, *_surface, *_surface_initial_change,
+! +*_vertical_gradient, *_vertical_gradient_level, *_vertical_gradient_level_ind,
+! +surface_waterflux_*
+!
+! in addition, steering parameters parameters of the two-moment cloud physics 
+! scheme:   
+! +cloud_scheme, +drizzle, +mu_constant, +mu_constant_value, +ventilation_effect
 !
 ! Former revisions:
@@ -151,5 +159,5 @@
 
 
-    binary_version = '3.5'
+    binary_version = '3.6'
 
     WRITE ( 14 )  binary_version
@@ -190,4 +198,10 @@
     WRITE ( 14 )  'bc_ns                         '
     WRITE ( 14 )  bc_ns
+    WRITE ( 14 )  'bc_nr_b                       '
+    WRITE ( 14 )  bc_nr_b
+    WRITE ( 14 )  'bc_nr_t                       '
+    WRITE ( 14 )  bc_nr_t
+    WRITE ( 14 )  'bc_nr_t_val                   '
+    WRITE ( 14 )  bc_nr_t_val
     WRITE ( 14 )  'bc_p_b                        '
     WRITE ( 14 )  bc_p_b
@@ -206,4 +220,10 @@
     WRITE ( 14 )  'bc_q_t_val                    '
     WRITE ( 14 )  bc_q_t_val
+    WRITE ( 14 )  'bc_qr_b                       '
+    WRITE ( 14 )  bc_qr_b
+    WRITE ( 14 )  'bc_qr_t                       '
+    WRITE ( 14 )  bc_qr_t
+    WRITE ( 14 )  'bc_qr_t_val                   '
+    WRITE ( 14 )  bc_qr_t_val
     WRITE ( 14 )  'bc_s_b                        '
     WRITE ( 14 )  bc_s_b
@@ -248,4 +268,6 @@
     WRITE ( 14 )  'cloud_physics                 '
     WRITE ( 14 )  cloud_physics
+    WRITE ( 14 )  'cloud_scheme                  '
+    WRITE ( 14 )  cloud_scheme
     WRITE ( 14 )  'collective_wait               '
     WRITE ( 14 )  collective_wait
@@ -278,4 +300,6 @@
     WRITE ( 14 )  'drag_coefficient              '
     WRITE ( 14 )  drag_coefficient
+    WRITE ( 14 )  'drizzle                       '
+    WRITE ( 14 )  drizzle
     WRITE ( 14 )  'dt_pr_1d                      '
     WRITE ( 14 )  dt_pr_1d
@@ -358,4 +382,8 @@
     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 )  'netcdf_precision              '
     WRITE ( 14 )  netcdf_precision
@@ -364,4 +392,16 @@
     WRITE ( 14 )  'ngsrb                         '
     WRITE ( 14 )  ngsrb
+    WRITE ( 14 )  'nr_init                       '
+    WRITE ( 14 )  nr_init
+    WRITE ( 14 )  'nr_surface                    '
+    WRITE ( 14 )  nr_surface
+    WRITE ( 14 )  'nr_surface_initial_change     '
+    WRITE ( 14 )  nr_surface_initial_change
+    WRITE ( 14 )  'nr_vertical_gradient          '
+    WRITE ( 14 )  nr_vertical_gradient
+    WRITE ( 14 )  'nr_vertical_gradient_level    '
+    WRITE ( 14 )  nr_vertical_gradient_level
+    WRITE ( 14 )  'nr_vertical_gradient_level_ind'
+    WRITE ( 14 )  nr_vertical_gradient_level_ind
     WRITE ( 14 )  'nsor                          '
     WRITE ( 14 )  nsor
@@ -428,4 +468,16 @@
     WRITE ( 14 )  'q_vertical_gradient_level_ind '
     WRITE ( 14 )  q_vertical_gradient_level_ind
+    WRITE ( 14 )  'qr_init                       '
+    WRITE ( 14 )  qr_init
+    WRITE ( 14 )  'qr_surface                    '
+    WRITE ( 14 )  qr_surface
+    WRITE ( 14 )  'qr_surface_initial_change     '
+    WRITE ( 14 )  qr_surface_initial_change
+    WRITE ( 14 )  'qr_vertical_gradient          '
+    WRITE ( 14 )  qr_vertical_gradient
+    WRITE ( 14 )  'qr_vertical_gradient_level    '
+    WRITE ( 14 )  qr_vertical_gradient_level
+    WRITE ( 14 )  'qr_vertical_gradient_level_ind'
+    WRITE ( 14 )  qr_vertical_gradient_level_ind
     WRITE ( 14 )  'radiation                     '
     WRITE ( 14 )  radiation
@@ -474,4 +526,8 @@
     WRITE ( 14 )  'surface_waterflux             '
     WRITE ( 14 )  surface_waterflux    
+    WRITE ( 14 )  'surface_waterflux_nr          '
+    WRITE ( 14 )  surface_waterflux_nr    
+    WRITE ( 14 )  'surface_waterflux_qr          '
+    WRITE ( 14 )  surface_waterflux_qr    
     WRITE ( 14 )  's_surface                     '
     WRITE ( 14 )  s_surface
@@ -574,4 +630,6 @@
     WRITE ( 14 )  'v_max_ijk                     '
     WRITE ( 14 )  v_max_ijk
+    WRITE ( 14 )  'ventilation_effect            '
+    WRITE ( 14 )  ventilation_effect
     WRITE ( 14 )  'vg                            '
     WRITE ( 14 )  vg
