Index: /palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2546)
+++ /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2547)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! extended by cloud_droplets option
+! 
+! 2532 2017-10-11 16:00:46Z scharf
 ! bugfixes in data_output_3d
 ! 
@@ -281,6 +284,6 @@
 
     USE control_parameters,                                                    &
-        ONLY:  cloud_physics, coupling_start_time, dt_3d, end_time, humidity,  &
-               intermediate_timestep_count,                                    &
+        ONLY:  cloud_droplets, cloud_physics, coupling_start_time, dt_3d,      &
+               end_time, humidity, intermediate_timestep_count,                &
                initializing_actions, intermediate_timestep_count_max,          &
                land_surface, max_masks, precipitation, pt_surface,             &
@@ -1531,4 +1534,7 @@
           pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
           qv1 = q(k,j,i) - ql(k,j,i)
+       ELSEIF ( cloud_droplets ) THEN
+          pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
+          qv1 = q(k,j,i) 
        ELSE
           pt1 = pt(k,j,i)
@@ -1950,5 +1956,5 @@
 !
 !--       Calculate specific humidity at saturation
-          q_s = 0.622_wp * e_s / surface_pressure
+          q_s = 0.622_wp * e_s / surface_pressure - e_s 
 
           resistance = surf%r_a(m) / ( surf%r_a(m) + surf%r_s(m) )
@@ -2491,5 +2497,5 @@
              ENDIF
              
-             IF ( cloud_physics )  THEN
+             IF ( cloud_physics  .OR.  cloud_droplets )  THEN
                 pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
              ELSE
@@ -2585,5 +2591,5 @@
                 k   = surf_lsm_v(l)%k(m)
 
-                IF ( cloud_physics )  THEN
+                IF ( cloud_physics  .OR.  cloud_droplets )  THEN
                    pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
                 ELSE
Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 2546)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 2547)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! extended by cloud_droplets option, minor bugfix and correct calculation of
+! cloud droplet number concentration
+! 
+! 2544 2017-10-13 18:09:32Z maronga
 ! Moved date and time quantitis to separate module date_and_time_mod 
 ! 
@@ -173,5 +177,5 @@
  
     USE arrays_3d,                                                             &
-        ONLY:  dzw, hyp, pt, q, ql, zu, zw
+        ONLY:  dzw, hyp, nc, pt, q, ql, zu, zw
 
     USE cloud_parameters,                                                      &
@@ -183,6 +187,7 @@
     USE control_parameters,                                                    &
         ONLY:  cloud_droplets, cloud_physics, g, initializing_actions,         &
-               large_scale_forcing, lsf_surf, phi, pt_surface, rho_surface,    &
-               surface_pressure, time_since_reference_point
+               large_scale_forcing, lsf_surf, microphysics_morrison, phi,      & 
+               pt_surface, rho_surface, surface_pressure,                      &
+               time_since_reference_point
 
     USE date_and_time_mod,                                                     &
@@ -1252,5 +1257,5 @@
              rad_lw_out(0,j,i) = emis(j,i) * sigma_sb * (pt(k,j,i) * exn)**4
 
-             IF ( cloud_physics )  THEN
+             IF ( cloud_physics  .OR.  cloud_droplets )  THEN
                 pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
                 rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
@@ -1303,5 +1308,5 @@
              exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
 
-             IF ( cloud_physics )  THEN
+             IF ( cloud_physics  .OR.  cloud_droplets )  THEN
                 pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
                 rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
@@ -1450,6 +1455,7 @@
        INTEGER(iwp) :: i, j, k, n !< loop indices
 
-       REAL(wp)     ::  s_r2, &   !< weighted sum over all droplets with r^2
-                        s_r3      !< weighted sum over all droplets with r^3
+       REAL(wp)     ::  nc_rad, &    !< number concentration of cloud droplets
+                        s_r2,   &    !< weighted sum over all droplets with r^2
+                        s_r3         !< weighted sum over all droplets with r^3
 
 !
@@ -1490,4 +1496,10 @@
                    rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
                 ENDDO
+             ELSEIF ( cloud_droplets ) THEN
+                DO k = nzb+1, nzt+1
+                   rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp      &
+                                    )**0.286_wp + l_d_cp * ql(k,j,i)
+                   rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q(k,j,i) 
+                ENDDO     
              ELSE
                 DO k = nzb+1, nzt+1
@@ -1531,5 +1543,5 @@
              rrtm_icld   = 0
 
-             IF ( cloud_physics )  THEN
+             IF ( cloud_physics  .OR.  cloud_droplets )  THEN
                 DO k = nzb+1, nzt+1
                    rrtm_cliqwp(0,k) =  ql(k,j,i) * 1000.0_wp *                 &
@@ -1544,7 +1556,8 @@
 !--                   Calculate cloud droplet effective radius
                       IF ( cloud_physics )  THEN
+                         nc_rad = MERGE( nc(k,j,i), nc_const, microphysics_morrison )
                          rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)     &
                                            * rho_surface                       &
-                                           / ( 4.0_wp * pi * nc_const * rho_l )&
+                                           / ( 4.0_wp * pi * nc_rad * rho_l )&
                                            )**0.33333333333333_wp              &
                                            * EXP( LOG( sigma_gc )**2 )
Index: /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 2546)
+++ /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 2547)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! extended by cloud_droplets option
+! 
+! 2321 2017-07-24 15:57:07Z schwenkel
 ! Bugfix: Correct index in lookup table for Obukhov length 
 ! 
@@ -201,7 +204,7 @@
 
     USE control_parameters,                                                    &
-        ONLY:  cloud_physics, constant_heatflux, constant_scalarflux,          &
-               constant_waterflux, coupling_mode, g, humidity, ibc_e_b,        &
-               ibc_pt_b, initializing_actions, kappa,                          &
+        ONLY:  cloud_droplets, cloud_physics, constant_heatflux,               &
+               constant_scalarflux, constant_waterflux, coupling_mode, g,      &
+               humidity, ibc_e_b, ibc_pt_b, initializing_actions, kappa,       &
                intermediate_timestep_count, intermediate_timestep_count_max,   &
                land_surface, large_scale_forcing, lsf_surf,                    &
@@ -291,5 +294,5 @@
 !--    temperature and specific humidity at first grid level from the fields pt 
 !--    and q
-       IF ( cloud_physics )  THEN
+       IF ( cloud_physics  .OR.  cloud_droplets )  THEN
 !
 !--       First call for horizontal default-type surfaces (l=0 - upward facing, 
@@ -1396,5 +1399,5 @@
 
              ENDDO
-          ELSEIF ( cloud_physics )  THEN
+          ELSEIF ( cloud_physics  .OR.  cloud_droplets )  THEN
              !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
              DO  m = 1, surf%ns
@@ -1529,6 +1532,11 @@
           k   = surf%k(m)
 
-          surf%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
-          surf%qv1(m) = q(k,j,i) - ql(k,j,i)
+          IF ( cloud_physics ) THEN
+             surf%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
+             surf%qv1(m) = q(k,j,i) - ql(k,j,i)
+          ELSEIF( cloud_droplets ) THEN
+             surf%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
+             surf%qv1(m) = q(k,j,i) 
+          ENDIF 
 
        ENDDO
@@ -1584,5 +1592,5 @@
           ENDIF
 
-          IF ( cloud_physics )  THEN
+          IF ( cloud_physics  .OR.  cloud_droplets )  THEN
              !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
              DO  m = 1, surf%ns   
@@ -1685,5 +1693,5 @@
              ENDIF
 
-             IF ( cloud_physics )  THEN
+             IF ( cloud_physics  .OR.  cloud_droplets )  THEN
                !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
                 DO  m = 1, surf%ns   
Index: /palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_mod.f90	(revision 2546)
+++ /palm/trunk/SOURCE/surface_mod.f90	(revision 2547)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! extended by cloud_droplets option
+! 
+! 2508 2017-10-02 08:57:09Z suehring
 ! Minor formatting adjustment
 ! 
@@ -139,6 +142,6 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  z0q       !< roughness length for humidity
 
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt1       !< Specific humidity at first grid level (required for cloud_physics = .T.)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qv1       !< Potential temperature at first grid level (required for cloud_physics = .T.)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt1       !< Specific humidity at first grid level (required for cloud_physics = .T. or cloud_droplets = .T.)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qv1       !< Potential temperature at first grid level (required for cloud_physics = .T. or cloud_droplets = .T.)
 !
 !--    Define arrays for surface fluxes
@@ -745,5 +748,5 @@
 !--    When cloud physics is used, arrays for storing potential temperature and
 !--    specific humidity at first grid level are required. 
-       IF ( cloud_physics )  THEN
+       IF ( cloud_physics  .OR.  cloud_droplets )  THEN
           ALLOCATE ( surfaces%pt1(1:surfaces%ns) )
           ALLOCATE ( surfaces%qv1(1:surfaces%ns) )
