Index: palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 3271)
+++ palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 3274)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Modularization of all bulk cloud physics code components
+! 
+! 3271 2018-09-24 08:20:34Z suehring
 ! Several bugfixes:
 ! - Initialization of pt_surface array with soil temperature in the uppermost
@@ -474,17 +477,17 @@
  
     USE arrays_3d,                                                             &
-        ONLY:  hyp, pt, prr, q, q_p, ql, vpt, u, v, w
+        ONLY:  hyp, pt, prr, q, q_p, ql, vpt, u, v, w, hyrho, exner, d_exner
+
+    USE basic_constants_and_equations_mod,                                     &
+        ONLY:  c_p, g, lv_d_cp, l_v, magnus, rho_l, r_d, r_v
 
     USE calc_mean_profile_mod,                                                 &
         ONLY:  calc_mean_profile
 
-    USE cloud_parameters,                                                      &
-        ONLY:  cp, hyrho, l_d_cp, l_d_r, l_v, pt_d_t, rho_l, r_d, r_v
-
     USE control_parameters,                                                    &
-        ONLY:  cloud_droplets, cloud_physics, coupling_start_time, dt_3d,      &
+        ONLY:  cloud_droplets, coupling_start_time, dt_3d,      &
                end_time, humidity, intermediate_timestep_count,                &
                initializing_actions, intermediate_timestep_count_max,          &
-               land_surface, max_masks, precipitation, pt_surface,             &
+               land_surface, max_masks, pt_surface,             &
                rho_surface, spinup, spinup_pt_mean, spinup_time,               &
                surface_pressure, timestep_scheme, tsc,                         &
@@ -493,4 +496,7 @@
     USE indices,                                                               &
         ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb
+
+    USE bulk_cloud_model_mod,                                                  &
+        ONLY: bulk_cloud_model, precipitation
 
     USE netcdf_data_input_mod,                                                 &
@@ -584,5 +590,4 @@
                 deep_soil_temperature =  9999999.9_wp,  & !< Deep soil temperature (bottom boundary condition)
                 drho_l_lv,                              & !< (rho_l * l_v)**-1
-                exn,                                    & !< value of the Exner function
                 field_capacity = 9999999.9_wp,          & !< NAMELIST m_fc
                 f_shortwave_incoming = 9999999.9_wp,    & !< NAMELIST f_sw_in
@@ -604,5 +609,4 @@
                 rho_cp,                                 & !< rho_surface * cp
                 rho_lv,                                 & !< rho_ocean * l_v
-                rd_d_rv,                                & !< r_d / r_v
                 saturation_moisture = 9999999.9_wp,     & !< NAMELIST m_sat
                 skip_time_do_lsm = 0.0_wp,              & !< LSM is not called before this time
@@ -1734,7 +1738,4 @@
  SUBROUTINE lsm_energy_balance( horizontal, l )
 
-    USE diagnostic_quantities_mod,                                             &
-        ONLY:  magnus 
-
     USE pegrid
 
@@ -1816,8 +1817,4 @@
     j_off = surf%joff
     i_off = surf%ioff
-    
-!
-!-- Calculate the exner function for the current time step
-    exn = ( surface_pressure / 1000.0_wp )**0.286_wp
 
     DO  m = 1, surf%ns
@@ -1889,9 +1886,9 @@
 !--    time step is used here. Note that this formulation is the
 !--    equivalent to the ECMWF formulation using drag coefficients
-!        IF ( cloud_physics )  THEN
-!           pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
+!        IF ( bulk_cloud_model )  THEN
+!           pt1 = pt(k,j,i) + lv_d_cp * d_exner(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)
+!           pt1 = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
 !           qv1 = q(k,j,i) 
 !        ELSE
@@ -2097,5 +2094,5 @@
 !--       Denominator of the prognostic equation
           coef_2 = surf%rad_lw_out_change_0(m) + f_qsws * dq_s_dt              &
-                   + lambda_surface + f_shf / exn
+                   + lambda_surface + f_shf / exner(nzb)
        ELSE
 !
@@ -2107,5 +2104,5 @@
 !
 !--       Denominator of the prognostic equation
-          coef_2 = surf%rad_lw_out_change_0(m) + lambda_surface + f_shf / exn
+          coef_2 = surf%rad_lw_out_change_0(m) + lambda_surface + f_shf / exner(nzb)
 
        ENDIF
@@ -2160,6 +2157,6 @@
 
 
-!        pt(k+k_off,j+j_off,i+i_off) = surf_t_surface_p%var_1d(m) / exn  !is actually no air temperature
-       surf%pt_surface(m)          = surf_t_surface_p%var_1d(m) / exn
+!        pt(k+k_off,j+j_off,i+i_off) = surf_t_surface_p%var_1d(m) / exner(nzb)  !is actually no air temperature
+       surf%pt_surface(m)          = surf_t_surface_p%var_1d(m) / exner(nzb)
 
 !
@@ -2177,5 +2174,5 @@
                                              - surf_t_soil%var_2d(nzb_soil,m) )
 
-       surf%shf(m) = - f_shf * ( surf%pt1(m) - surf%pt_surface(m) ) / cp
+       surf%shf(m) = - f_shf * ( surf%pt1(m) - surf%pt_surface(m) ) / c_p
 
        IF ( humidity )  THEN
@@ -2330,9 +2327,9 @@
     SUBROUTINE calc_q_surface
 
-       USE diagnostic_quantities_mod
-
        IMPLICIT NONE
 
-       REAL(wp) :: resistance    !< aerodynamic and soil resistance term
+       REAL(wp) ::  e_s           !< saturation water vapor pressure
+       REAL(wp) ::  q_s           !< saturation mixing ratio
+       REAL(wp) ::  resistance    !< aerodynamic and soil resistance term
 
        DO  m = 1, surf%ns
@@ -2353,5 +2350,5 @@
 !
 !--       Calculate mixing ratio at surface
-          IF ( cloud_physics )  THEN
+          IF ( bulk_cloud_model )  THEN
              q(k+k_off,j+j_off,i+i_off) = resistance * q_s +                   &
                                         ( 1.0_wp - resistance ) *              &
@@ -2471,5 +2468,5 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE lsm_init
-    
+
        USE control_parameters,                                                 &
            ONLY:  message_string
@@ -2501,19 +2498,14 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  bound, bound_root_fr  !< temporary arrays for storing index bounds
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  pr_soil_init !< temporary array used for averaging soil profiles
-
-!
-!--    Calculate Exner function
-       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
 !
 !--    If no cloud physics is used, rho_surface has not been calculated before
-       IF (  .NOT.  cloud_physics  .AND.  .NOT.  cloud_droplets )  THEN
+       IF (  .NOT.  bulk_cloud_model  .AND.  .NOT.  cloud_droplets )  THEN
           CALL calc_mean_profile( pt, 4 )
-          rho_surface = surface_pressure * 100.0_wp / ( r_d * hom(topo_min_level+1,1,4,0) * exn )
+          rho_surface = hyp(nzb) / ( r_d * hom(topo_min_level+1,1,4,0) * exner(nzb) )
        ENDIF
 
 !
 !--    Calculate frequently used parameters
-       rho_cp    = cp * rho_surface
-       rd_d_rv   = r_d / r_v
+       rho_cp    = c_p * rho_surface
        rho_lv    = rho_surface * l_v
        drho_l_lv = 1.0_wp / (rho_l * l_v)
@@ -4451,8 +4443,8 @@
 !--          uppermost layer
              t_surface_h%var_1d(m)    = t_soil_h%var_2d(nzb_soil,m)
-             surf_lsm_h%pt_surface(m) = t_soil_h%var_2d(nzb_soil,m) / exn
+             surf_lsm_h%pt_surface(m) = t_soil_h%var_2d(nzb_soil,m) / exner(nzb)
              
-             IF ( cloud_physics  .OR. cloud_droplets ) THEN
-                surf_lsm_h%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
+             IF ( bulk_cloud_model  .OR. cloud_droplets ) THEN
+                surf_lsm_h%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
              ELSE
                 surf_lsm_h%pt1(m) = pt(k,j,i)
@@ -4480,8 +4472,8 @@
 !--             uppermost layer
                 t_surface_v(l)%var_1d(m)      = t_soil_v(l)%var_2d(nzb_soil,m)
-                surf_lsm_v(l)%pt_surface(m)   = t_soil_v(l)%var_2d(nzb_soil,m) / exn
-
-                IF ( cloud_physics  .OR. cloud_droplets ) THEN
-                   surf_lsm_v(l)%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
+                surf_lsm_v(l)%pt_surface(m)   = t_soil_v(l)%var_2d(nzb_soil,m) / exner(nzb)
+
+                IF ( bulk_cloud_model  .OR. cloud_droplets ) THEN
+                   surf_lsm_v(l)%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
                 ELSE
                    surf_lsm_v(l)%pt1(m) = pt(k,j,i)
@@ -7078,5 +7070,5 @@
 
        USE control_parameters,                                                 &
-           ONLY: g, message_string, molecular_viscosity
+           ONLY: message_string, molecular_viscosity
 
        IMPLICIT NONE
