Index: palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3245)
+++ palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3274)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Modularization of all bulk cloud physics code components
+! 
+! 3245 2018-09-13 14:08:16Z knoop
 ! unused variables removed, shortest_distance has wp now
 ! 
@@ -154,8 +157,11 @@
 #endif
 
+    USE basic_constants_and_equations_mod,                                     &
+        ONLY:  g, kappa, lv_d_cp, lv_d_rd
+
     USE control_parameters,                                                    &
         ONLY:  constant_diffusion, dt_3d, e_init, humidity,                    &
                initializing_actions, intermediate_timestep_count,              &
-               intermediate_timestep_count_max, kappa, km_constant,            &
+               intermediate_timestep_count_max, km_constant,                   &
                les_dynamic, les_mw, ocean, plant_canopy, prandtl_number,       &
                prho_reference, pt_reference, rans_mode, rans_tke_e, rans_tke_l,&
@@ -950,5 +956,5 @@
  SUBROUTINE tcm_init_arrays
 
-    USE microphysics_mod,                                                      &
+    USE bulk_cloud_model_mod,                                                  &
         ONLY:  collision_turbulence
 
@@ -2598,11 +2604,8 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  ddzw, dd2zu, drho_air_zw, q, ql
-
-    USE cloud_parameters,                                                      &
-        ONLY:  l_d_cp, l_d_r, pt_d_t, t_d_pt
+        ONLY:  ddzw, dd2zu, drho_air_zw, q, ql, d_exner, exner
 
     USE control_parameters,                                                    &
-        ONLY:  cloud_droplets, cloud_physics, constant_flux_layer, g, neutral, &
+        ONLY:  cloud_droplets, constant_flux_layer, neutral,                   &
                rho_reference, use_single_reference_value, use_surface_fluxes,  &
                use_top_fluxes
@@ -2610,4 +2613,7 @@
     USE grid_variables,                                                        &
         ONLY:  ddx, dx, ddy, dy
+
+    USE bulk_cloud_model_mod,                                                  &
+        ONLY:  bulk_cloud_model
 
     USE surface_mod,                                                           &
@@ -3067,5 +3073,5 @@
 !--                Flag 9 is used to mask top fluxes, flag 30 to mask
 !--                surface fluxes
-                   IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+                   IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
                       k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                       k2 = 0.61_wp * pt(k,j,i)
@@ -3083,17 +3089,17 @@
                                           BTEST( wall_flags_0(k,j,i), 9 )      &
                                         )
-                   ELSE IF ( cloud_physics )  THEN
+                   ELSE IF ( bulk_cloud_model )  THEN
                       IF ( ql(k,j,i) == 0.0_wp )  THEN
                          k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                          k2 = 0.61_wp * pt(k,j,i)
                       ELSE
-                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                         temp  = theta * t_d_pt(k)
+                         theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
+                         temp  = theta * exner(k)
                          k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                  &
                                        ( q(k,j,i) - ql(k,j,i) ) *              &
-                              ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /         &
-                              ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *           &
+                              ( 1.0_wp + 0.622_wp * lv_d_rd / temp ) ) /       &
+                              ( 1.0_wp + 0.622_wp * lv_d_rd * lv_d_cp *        &
                               ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                         k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                         k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
                       ENDIF
                       tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *                  &
@@ -3144,20 +3150,20 @@
                          k = surf_def_h(l)%k(m)
 
-                         IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+                         IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
                             k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                             k2 = 0.61_wp * pt(k,j,i)
-                         ELSE IF ( cloud_physics )  THEN
+                         ELSE IF ( bulk_cloud_model )  THEN
                             IF ( ql(k,j,i) == 0.0_wp )  THEN
                                k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                                k2 = 0.61_wp * pt(k,j,i)
                             ELSE
-                               theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                               temp  = theta * t_d_pt(k)
+                               theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
+                               temp  = theta * exner(k)
                                k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *            &
                                           ( q(k,j,i) - ql(k,j,i) ) *           &
-                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /      &
-                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *        &
+                                 ( 1.0_wp + 0.622_wp * lv_d_rd / temp ) ) /    &
+                                 ( 1.0_wp + 0.622_wp * lv_d_rd * lv_d_cp *     &
                                  ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                               k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                               k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
                             ENDIF
                          ELSE IF ( cloud_droplets )  THEN
@@ -3181,20 +3187,20 @@
                       k = surf_lsm_h%k(m)
 
-                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+                      IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
                          k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                          k2 = 0.61_wp * pt(k,j,i)
-                      ELSE IF ( cloud_physics )  THEN
+                      ELSE IF ( bulk_cloud_model )  THEN
                          IF ( ql(k,j,i) == 0.0_wp )  THEN
                             k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                             k2 = 0.61_wp * pt(k,j,i)
                          ELSE
-                            theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                            temp  = theta * t_d_pt(k)
+                            theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
+                            temp  = theta * exner(k)
                             k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
                                           ( q(k,j,i) - ql(k,j,i) ) *           &
-                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /      &
-                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *        &
+                                 ( 1.0_wp + 0.622_wp * lv_d_rd / temp ) ) /    &
+                                 ( 1.0_wp + 0.622_wp * lv_d_rd * lv_d_cp *     &
                                  ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                            k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
                          ENDIF
                       ELSE IF ( cloud_droplets )  THEN
@@ -3217,20 +3223,20 @@
                       k = surf_usm_h%k(m)
 
-                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+                      IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
                          k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                          k2 = 0.61_wp * pt(k,j,i)
-                      ELSE IF ( cloud_physics )  THEN
+                      ELSE IF ( bulk_cloud_model )  THEN
                          IF ( ql(k,j,i) == 0.0_wp )  THEN
                             k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                             k2 = 0.61_wp * pt(k,j,i)
                          ELSE
-                            theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                            temp  = theta * t_d_pt(k)
+                            theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
+                            temp  = theta * exner(k)
                             k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
                                           ( q(k,j,i) - ql(k,j,i) ) *           &
-                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /      &
-                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *        &
+                                 ( 1.0_wp + 0.622_wp * lv_d_rd / temp ) ) /    &
+                                 ( 1.0_wp + 0.622_wp * lv_d_rd * lv_d_cp *     &
                                  ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                            k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
                          ENDIF
                       ELSE IF ( cloud_droplets )  THEN
@@ -3260,20 +3266,20 @@
                       k = surf_def_h(2)%k(m)
 
-                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+                      IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
                          k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                          k2 = 0.61_wp * pt(k,j,i)
-                      ELSE IF ( cloud_physics )  THEN
+                      ELSE IF ( bulk_cloud_model )  THEN
                          IF ( ql(k,j,i) == 0.0_wp )  THEN
                             k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                             k2 = 0.61_wp * pt(k,j,i)
                          ELSE
-                            theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                            temp  = theta * t_d_pt(k)
+                            theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
+                            temp  = theta * exner(k)
                             k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
                                        ( q(k,j,i) - ql(k,j,i) ) *              &
-                              ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /         &
-                              ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *           &
+                              ( 1.0_wp + 0.622_wp * lv_d_rd / temp ) ) /       &
+                              ( 1.0_wp + 0.622_wp * lv_d_rd * lv_d_cp *        &
                               ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                            k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
                          ENDIF
                       ELSE IF ( cloud_droplets )  THEN
@@ -3316,11 +3322,8 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  ddzw, dd2zu, drho_air_zw, q, ql
-
-    USE cloud_parameters,                                                      &
-        ONLY:  l_d_cp, l_d_r, pt_d_t, t_d_pt
+        ONLY:  ddzw, dd2zu, drho_air_zw, q, ql, d_exner, exner
 
     USE control_parameters,                                                    &
-        ONLY:  cloud_droplets, cloud_physics, constant_flux_layer, g, neutral, &
+        ONLY:  cloud_droplets, constant_flux_layer, neutral,                &
                rho_reference, use_single_reference_value, use_surface_fluxes,  &
                use_top_fluxes
@@ -3328,4 +3331,7 @@
     USE grid_variables,                                                        &
         ONLY:  ddx, dx, ddy, dy
+
+    USE bulk_cloud_model_mod,                                                  &
+        ONLY:  bulk_cloud_model
 
     USE surface_mod,                                                           &
@@ -3864,5 +3870,5 @@
 !
 !--          Flag 9 is used to mask top fluxes, flag 30 to mask surface fluxes
-             IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
+             IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets )  THEN
                 k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                 k2 = 0.61_wp * pt(k,j,i)
@@ -3880,17 +3886,17 @@
                                         )
 
-             ELSE IF ( cloud_physics )  THEN
+             ELSE IF ( bulk_cloud_model )  THEN
                 IF ( ql(k,j,i) == 0.0_wp )  THEN
                    k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                    k2 = 0.61_wp * pt(k,j,i)
                 ELSE
-                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                   temp  = theta * t_d_pt(k)
+                   theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
+                   temp  = theta * exner(k)
                    k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                        &
                                  ( q(k,j,i) - ql(k,j,i) ) *                    &
-                        ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /               &
-                        ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *                 &
+                        ( 1.0_wp + 0.622_wp * lv_d_rd / temp ) ) /             &
+                        ( 1.0_wp + 0.622_wp * lv_d_rd * lv_d_cp *              &
                         ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                   k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                   k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
                 ENDIF
                 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g /                    &
@@ -3934,20 +3940,20 @@
                    k = surf_def_h(l)%k(m)
 
-                   IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+                   IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
                       k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                       k2 = 0.61_wp * pt(k,j,i)
-                   ELSE IF ( cloud_physics )  THEN
+                   ELSE IF ( bulk_cloud_model )  THEN
                       IF ( ql(k,j,i) == 0.0_wp )  THEN
                          k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                          k2 = 0.61_wp * pt(k,j,i)
                       ELSE
-                        theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                        temp  = theta * t_d_pt(k)
+                        theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
+                        temp  = theta * exner(k)
                         k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                   &
                                    ( q(k,j,i) - ql(k,j,i) ) *                  &
-                          ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /             &
-                          ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *               & 
+                          ( 1.0_wp + 0.622_wp * lv_d_rd / temp ) ) /           &
+                          ( 1.0_wp + 0.622_wp * lv_d_rd * lv_d_cp *            &
                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                        k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                        k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
                       ENDIF
                    ELSE IF ( cloud_droplets )  THEN
@@ -3971,20 +3977,20 @@
                 k = surf_lsm_h%k(m)
 
-                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+                IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
                     k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                     k2 = 0.61_wp * pt(k,j,i)
-                ELSE IF ( cloud_physics )  THEN
+                ELSE IF ( bulk_cloud_model )  THEN
                     IF ( ql(k,j,i) == 0.0_wp )  THEN
                        k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                        k2 = 0.61_wp * pt(k,j,i)
                     ELSE
-                       theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                       temp  = theta * t_d_pt(k)
+                       theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
+                       temp  = theta * exner(k)
                        k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                    &
                                      ( q(k,j,i) - ql(k,j,i) ) *                &
-                            ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /           &
-                            ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *             &
+                            ( 1.0_wp + 0.622_wp * lv_d_rd / temp ) ) /         &
+                            ( 1.0_wp + 0.622_wp * lv_d_rd * lv_d_cp *          &
                             ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                       k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                       k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
                    ENDIF
                 ELSE IF ( cloud_droplets )  THEN
@@ -4007,20 +4013,20 @@
                 k = surf_usm_h%k(m)
 
-                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+                IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
                     k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                     k2 = 0.61_wp * pt(k,j,i)
-                ELSE IF ( cloud_physics )  THEN
+                ELSE IF ( bulk_cloud_model )  THEN
                     IF ( ql(k,j,i) == 0.0_wp )  THEN
                        k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                        k2 = 0.61_wp * pt(k,j,i)
                     ELSE
-                       theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                       temp  = theta * t_d_pt(k)
+                       theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
+                       temp  = theta * exner(k)
                        k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                    &
                                      ( q(k,j,i) - ql(k,j,i) ) *                &
-                            ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /           &
-                            ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *             &
+                            ( 1.0_wp + 0.622_wp * lv_d_rd / temp ) ) /         &
+                            ( 1.0_wp + 0.622_wp * lv_d_rd * lv_d_cp *          &
                             ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                       k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                       k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
                    ENDIF
                 ELSE IF ( cloud_droplets )  THEN
@@ -4047,20 +4053,20 @@
 
 
-                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
+                IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets )  THEN
                    k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                    k2 = 0.61_wp * pt(k,j,i)
-                ELSE IF ( cloud_physics )  THEN
+                ELSE IF ( bulk_cloud_model )  THEN
                    IF ( ql(k,j,i) == 0.0_wp )  THEN
                       k1 = 1.0_wp + 0.61_wp * q(k,j,i)
                       k2 = 0.61_wp * pt(k,j,i)
                    ELSE
-                      theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                      temp  = theta * t_d_pt(k)
+                      theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
+                      temp  = theta * exner(k)
                       k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                     &
                                  ( q(k,j,i) - ql(k,j,i) ) *                    &
-                        ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /               &
-                        ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *                 &
+                        ( 1.0_wp + 0.622_wp * lv_d_rd / temp ) ) /             &
+                        ( 1.0_wp + 0.622_wp * lv_d_rd * lv_d_cp *              &
                         ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
+                      k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
                    ENDIF
                 ELSE IF ( cloud_droplets )  THEN
@@ -4100,5 +4106,5 @@
         ONLY:  ddx2, ddy2
 
-    USE microphysics_mod,                                                      &
+    USE bulk_cloud_model_mod,                                                  &
         ONLY:  collision_turbulence
 
@@ -4237,5 +4243,5 @@
         ONLY:  ddx2, ddy2
         
-    USE microphysics_mod,                                                      &
+    USE bulk_cloud_model_mod,                                                  &
         ONLY:  collision_turbulence
 
@@ -4476,5 +4482,5 @@
 
     USE control_parameters,                                                    &
-        ONLY:  atmos_ocean_sign, g, use_single_reference_value,                &
+        ONLY:  atmos_ocean_sign, use_single_reference_value,                   &
                wall_adjustment, wall_adjustment_factor
 
@@ -4533,5 +4539,5 @@
 
     USE control_parameters,                                                    &
-        ONLY:  atmos_ocean_sign, g, use_single_reference_value
+        ONLY:  atmos_ocean_sign, use_single_reference_value
 
     IMPLICIT NONE
