Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 1321)
+++ palm/trunk/SOURCE/production_e.f90	(revision 1342)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -173,27 +173,27 @@
              DO  k = nzb_diff_s_outer(j,i), nzt
 
-                dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-                dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+                dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+                dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                       dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                IF ( def < 0.0 )  def = 0.0
+                      dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                IF ( def < 0.0_wp )  def = 0.0_wp
 
                 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -210,17 +210,17 @@
              DO  j = nys, nyn
 
-                IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
+                IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) &
                 THEN
 
                    k = nzb_diff_s_inner(j,i) - 1
                    dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
-                   dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                  u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
+                   dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                     u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
                    dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
-                   dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                  v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
+                   dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                     v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
                    dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
 
-                   IF ( wall_e_y(j,i) /= 0.0 )  THEN
+                   IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
 !
 !--                   Inconsistency removed: as the thermal stratification is 
@@ -236,21 +236,21 @@
                       CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                           wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
-                      km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25 * &
-                                   0.5 * dy 
-                      IF ( km_neutral > 0.0 )  THEN
+                      km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25_wp * &
+                                   0.5_wp * dy 
+                      IF ( km_neutral > 0.0_wp )  THEN
                          dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
                          dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
                       ELSE
-                         dudy = 0.0
-                         dwdy = 0.0
+                         dudy = 0.0_wp
+                         dwdy = 0.0_wp
                       ENDIF
                    ELSE
-                      dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                      u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                      dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                      w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                      dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                         u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                      dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                         w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
                    ENDIF
 
-                   IF ( wall_e_x(j,i) /= 0.0 )  THEN
+                   IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
 !
 !--                   Inconsistency removed: as the thermal stratification is 
@@ -266,25 +266,25 @@
                       CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                           wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
-                      km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25 * &
-                                   0.5 * dx
-                      IF ( km_neutral > 0.0 )  THEN
+                      km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25_wp * &
+                                   0.5_wp * dx
+                      IF ( km_neutral > 0.0_wp )  THEN
                          dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
                          dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
                       ELSE
-                         dvdx = 0.0
-                         dwdx = 0.0
+                         dvdx = 0.0_wp
+                         dwdx = 0.0_wp
                       ENDIF
                    ELSE
-                      dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                      v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                      dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                      w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                      dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                         v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                      dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                         w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
                    ENDIF
 
-                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                          dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                   IF ( def < 0.0 )  def = 0.0
+                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                   IF ( def < 0.0_wp )  def = 0.0_wp
 
                    tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -301,12 +301,12 @@
 
                       dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
-                      dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                     u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-                      dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                      dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+                      dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                        u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+                      dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                      dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                        v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
                       dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
 
-                      IF ( wall_e_y(j,i) /= 0.0 )  THEN
+                      IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
 !
 !--                      Inconsistency removed: as the thermal stratification 
@@ -319,20 +319,20 @@
 !--                            validation has been available
                          km_neutral = kappa * ( usvs(k)**2 + & 
-                                                wsvs(k)**2 )**0.25 * 0.5 * dy
-                         IF ( km_neutral > 0.0 )  THEN
+                                                wsvs(k)**2 )**0.25_wp * 0.5_wp * dy
+                         IF ( km_neutral > 0.0_wp )  THEN
                             dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
                             dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
                          ELSE
-                            dudy = 0.0
-                            dwdy = 0.0
+                            dudy = 0.0_wp
+                            dwdy = 0.0_wp
                          ENDIF
                       ELSE
-                         dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                         u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                         dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                         w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                         dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                            u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                         dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                            w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
                       ENDIF
 
-                      IF ( wall_e_x(j,i) /= 0.0 )  THEN
+                      IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
 !
 !--                      Inconsistency removed: as the thermal stratification 
@@ -345,24 +345,24 @@
 !--                            validation has been available
                          km_neutral = kappa * ( vsus(k)**2 + & 
-                                                wsus(k)**2 )**0.25 * 0.5 * dx
-                         IF ( km_neutral > 0.0 )  THEN
+                                                wsus(k)**2 )**0.25_wp * 0.5_wp * dx
+                         IF ( km_neutral > 0.0_wp )  THEN
                             dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
                             dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
                          ELSE
-                            dvdx = 0.0
-                            dwdx = 0.0
+                            dvdx = 0.0_wp
+                            dwdx = 0.0_wp
                          ENDIF
                       ELSE
-                         dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                         v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                         dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                         w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                         dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                            v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                         dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                            w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
                       ENDIF
 
-                      def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                      def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                            dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 +  &
-                           dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                      IF ( def < 0.0 )  def = 0.0
+                           dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                      IF ( def < 0.0_wp )  def = 0.0_wp
 
                       tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -379,32 +379,32 @@
              DO  j = nys, nyn
 
-                IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
+                IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) &
                 THEN
 
                    k = nzb_diff_s_outer(j,i)-1
 
-                   dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                   dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                   dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-                   dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                   dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                   dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                   dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                   dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                   dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                   dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                   dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                   dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                       u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+                   dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                   dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                   dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                       v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+                   dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                   dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                   dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                          dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                   IF ( def < 0.0 )  def = 0.0
+                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                   IF ( def < 0.0_wp )  def = 0.0_wp
 
                    tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -420,68 +420,68 @@
              DO  j = nys, nyn
 
-                IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
+                IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) &
                 THEN
 
                    k = nzb_diff_s_inner(j,i)-1
 
-                   dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                   dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                   dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                   dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                   dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                       u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
+
+                   dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                   dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                   dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                       v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
+
+                   dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                   dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                   dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
+                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                   IF ( def < 0.0_wp )  def = 0.0_wp
+
+                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
+
+                ENDIF
+
+             ENDDO
+
+          ELSEIF ( use_surface_fluxes )  THEN
+
+             DO  j = nys, nyn
+
+                k = nzb_diff_s_outer(j,i)-1
+
+                dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
                                     u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                   dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                    u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
-
-                   dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                   u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+                dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
                                     v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                   dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                   dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                    v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
-
-                   dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+                dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
                                     w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                   dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
                                     w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                   dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
-                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                   IF ( def < 0.0 )  def = 0.0
-
-                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-
-                ENDIF
-
-             ENDDO
-
-          ELSEIF ( use_surface_fluxes )  THEN
-
-             DO  j = nys, nyn
-
-                k = nzb_diff_s_outer(j,i)-1
-
-                dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-                dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                       dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                IF ( def < 0.0 )  def = 0.0
+                      dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                IF ( def < 0.0_wp )  def = 0.0_wp
 
                 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -586,6 +586,6 @@
 
                       IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
-                         k1 = 1.0 + 0.61 * q(k,j,i)
-                         k2 = 0.61 * pt(k,j,i)
+                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                         k2 = 0.61_wp * pt(k,j,i)
                          tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
                                          g / vpt(k,j,i) *                      &
@@ -594,16 +594,16 @@
                                          ) * dd2zu(k)
                       ELSE IF ( cloud_physics )  THEN
-                         IF ( ql(k,j,i) == 0.0 )  THEN
-                            k1 = 1.0 + 0.61 * q(k,j,i)
-                            k2 = 0.61 * pt(k,j,i)
+                         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)
-                            k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                                       ( q(k,j,i) - ql(k,j,i) ) *          &
-                                 ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                                 ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                            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 *          &
                                  ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                            k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
                          ENDIF
                          tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
@@ -613,6 +613,6 @@
                                          ) * dd2zu(k)
                       ELSE IF ( cloud_droplets )  THEN
-                         k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-                         k2 = 0.61 * pt(k,j,i)
+                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                         k2 = 0.61_wp * pt(k,j,i)
                          tend(k,j,i) = tend(k,j,i) -                          & 
                                        kh(k,j,i) * g / vpt(k,j,i) *           &
@@ -634,23 +634,23 @@
 
                       IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
-                         k1 = 1.0 + 0.61 * q(k,j,i)
-                         k2 = 0.61 * pt(k,j,i)
+                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                         k2 = 0.61_wp * pt(k,j,i)
                       ELSE IF ( cloud_physics )  THEN
-                         IF ( ql(k,j,i) == 0.0 )  THEN
-                            k1 = 1.0 + 0.61 * q(k,j,i)
-                            k2 = 0.61 * pt(k,j,i)
+                         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)
-                            k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                                       ( q(k,j,i) - ql(k,j,i) ) *          &
-                                 ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                                 ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                            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 *          &
                                  ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                            k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
                          ENDIF
                       ELSE IF ( cloud_droplets )  THEN
-                         k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-                         k2 = 0.61 * pt(k,j,i)
+                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                         k2 = 0.61_wp * pt(k,j,i)
                       ENDIF
 
@@ -668,23 +668,23 @@
 
                       IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
-                         k1 = 1.0 + 0.61 * q(k,j,i)
-                         k2 = 0.61 * pt(k,j,i)
+                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                         k2 = 0.61_wp * pt(k,j,i)
                       ELSE IF ( cloud_physics )  THEN
-                         IF ( ql(k,j,i) == 0.0 )  THEN
-                            k1 = 1.0 + 0.61 * q(k,j,i)
-                            k2 = 0.61 * pt(k,j,i)
+                         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)
-                            k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                                       ( q(k,j,i) - ql(k,j,i) ) *          &
-                                 ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                                 ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                            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 *          &
                                  ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                            k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
                          ENDIF
                       ELSE IF ( cloud_droplets )  THEN
-                         k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-                         k2 = 0.61 * pt(k,j,i)
+                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                         k2 = 0.61_wp * pt(k,j,i)
                       ENDIF
 
@@ -783,27 +783,27 @@
                 IF ( k >= nzb_diff_s_outer(j,i) )  THEN
 
-                   dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                   dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                   dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-                   dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                   dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                   dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                   dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                   dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                   dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                   dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                   dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                   dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                       u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+                   dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                   dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                   dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                       v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+                   dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                   dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                   dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                          dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                   IF ( def < 0.0 )  def = 0.0
+                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                   IF ( def < 0.0_wp )  def = 0.0_wp
 
                    tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -825,17 +825,17 @@
                 DO  k = 1, nzt
 
-                   IF ( ( wall_e_x(j,i) /= 0.0 ).OR.( wall_e_y(j,i) /= 0.0 ) ) &
+                   IF ( ( wall_e_x(j,i) /= 0.0_wp ).OR.( wall_e_y(j,i) /= 0.0_wp ) ) &
                    THEN
 
                       IF ( k == nzb_diff_s_inner(j,i) - 1 )  THEN
                          dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
-                         dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                        u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
+                         dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                           u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
                          dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
-                         dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                        v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
+                         dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                           v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
                          dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
 
-                         IF ( wall_e_y(j,i) /= 0.0 )  THEN
+                         IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
 !
 !--                         Inconsistency removed: as the thermal stratification is
@@ -852,21 +852,21 @@
 !                                                wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
                             km_neutral = kappa *                                    &
-                                        ( usvs(k,j,i)**2 + wsvs(k,j,i)**2 )**0.25 * &
-                                         0.5 * dy
-                            IF ( km_neutral > 0.0 )  THEN
+                                        ( usvs(k,j,i)**2 + wsvs(k,j,i)**2 )**0.25_wp * &
+                                         0.5_wp * dy
+                            IF ( km_neutral > 0.0_wp )  THEN
                                dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral
                                dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral
                             ELSE
-                               dudy = 0.0
-                               dwdy = 0.0
+                               dudy = 0.0_wp
+                               dwdy = 0.0_wp
                             ENDIF
                          ELSE
-                            dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                            u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                            dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                            w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                            dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                               u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                            dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                               w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
                          ENDIF
 
-                         IF ( wall_e_x(j,i) /= 0.0 )  THEN
+                         IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
 !
 !--                         Inconsistency removed: as the thermal stratification is
@@ -883,25 +883,25 @@
 !                                                wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
                             km_neutral = kappa *                                     &
-                                         ( vsus(k,j,i)**2 + wsus(k,j,i)**2 )**0.25 * &
-                                         0.5 * dx
-                            IF ( km_neutral > 0.0 )  THEN
+                                         ( vsus(k,j,i)**2 + wsus(k,j,i)**2 )**0.25_wp * &
+                                         0.5_wp * dx
+                            IF ( km_neutral > 0.0_wp )  THEN
                                dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral
                                dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral
                             ELSE
-                               dvdx = 0.0
-                               dwdx = 0.0
+                               dvdx = 0.0_wp
+                               dwdx = 0.0_wp
                             ENDIF
                          ELSE
-                            dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                            v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                            dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                            w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                            dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                               v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                            dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                               w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
                          ENDIF
 
-                         def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                         def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                                dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                               dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                         IF ( def < 0.0 )  def = 0.0
+                               dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                         IF ( def < 0.0_wp )  def = 0.0_wp
 
                          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -919,12 +919,12 @@
 
                          dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
-                         dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                        u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-                         dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                         dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                        v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+                         dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+                         dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                         dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
                          dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
 
-                         IF ( wall_e_y(j,i) /= 0.0 )  THEN
+                         IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
 !
 !--                         Inconsistency removed: as the thermal stratification
@@ -937,20 +937,20 @@
 !--                               validation has been available
                             km_neutral = kappa * ( usvs(k,j,i)**2 + &
-                                                   wsvs(k,j,i)**2 )**0.25 * 0.5 * dy
-                            IF ( km_neutral > 0.0 )  THEN
+                                                   wsvs(k,j,i)**2 )**0.25_wp * 0.5_wp * dy
+                            IF ( km_neutral > 0.0_wp )  THEN
                                dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral
                                dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral
                             ELSE
-                               dudy = 0.0
-                               dwdy = 0.0
+                               dudy = 0.0_wp
+                               dwdy = 0.0_wp
                             ENDIF
                          ELSE
-                            dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                            u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                            dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                            w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                            dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                               u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                            dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                               w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
                          ENDIF
 
-                         IF ( wall_e_x(j,i) /= 0.0 )  THEN
+                         IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
 !
 !--                         Inconsistency removed: as the thermal stratification
@@ -963,24 +963,24 @@
 !--                               validation has been available
                             km_neutral = kappa * ( vsus(k,j,i)**2 + &
-                                                   wsus(k,j,i)**2 )**0.25 * 0.5 * dx
-                            IF ( km_neutral > 0.0 )  THEN
+                                                   wsus(k,j,i)**2 )**0.25_wp * 0.5_wp * dx
+                            IF ( km_neutral > 0.0_wp )  THEN
                                dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral
                                dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral
                             ELSE
-                               dvdx = 0.0
-                               dwdx = 0.0
+                               dvdx = 0.0_wp
+                               dwdx = 0.0_wp
                             ENDIF
                          ELSE
-                            dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                            v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                            dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                            w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                            dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                               v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                            dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                               w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
                          ENDIF
 
-                         def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                         def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                               dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 +  &
-                              dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                         IF ( def < 0.0 )  def = 0.0
+                              dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                         IF ( def < 0.0_wp )  def = 0.0_wp
 
                          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -993,27 +993,27 @@
                       IF ( k == nzb_diff_s_outer(j,i)-1 )  THEN
 
-                         dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                         dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                          u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                         dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                          u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-                         dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                          v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                         dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                         dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                          v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                         dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                          w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                         dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                          w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                         dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                         def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                         dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                         dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                             u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                         dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                             u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+                         dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                             v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                         dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                         dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                             v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+                         dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                             w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                         dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                             w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                         dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                         def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                                dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                               dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                         IF ( def < 0.0 )  def = 0.0
+                               dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                         IF ( def < 0.0_wp )  def = 0.0_wp
 
                          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -1035,32 +1035,32 @@
                 DO  k = 1, nzt
 
-                   IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
+                   IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) &
                    THEN
 
                       IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
 
-                         dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                         dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                          u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                         dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                          u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
-
-                         dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                          v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                         dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                         dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                          v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
-
-                         dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                          w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                         dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                          w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                         dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                         def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                         dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                         dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                             u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                         dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                             u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
+
+                         dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                             v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                         dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                         dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                             v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
+
+                         dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                             w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                         dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                             w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                         dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                         def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                                dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                               dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                         IF ( def < 0.0 )  def = 0.0
+                               dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                         IF ( def < 0.0_wp )  def = 0.0_wp
 
                          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -1082,27 +1082,27 @@
                    IF ( k == nzb_diff_s_outer(j,i)-1 )  THEN
 
-                      dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                      dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                      dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                       u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-                      dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                      dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                      dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                       v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                      dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                      dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                      dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                      def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                      dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                      dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                          u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                      dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                          u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+                      dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                          v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                      dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                      dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                          v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+                      dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                          w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                      dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                          w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                      dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                      def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                             dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                            dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                      IF ( def < 0.0 )  def = 0.0
+                            dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                      IF ( def < 0.0_wp )  def = 0.0_wp
 
                       tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -1232,6 +1232,6 @@
 !
 !                         IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
-!                            k1 = 1.0 + 0.61 * q(k,j,i)
-!                            k2 = 0.61 * pt(k,j,i)
+!                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+!                            k2 = 0.61_wp * pt(k,j,i)
 !                            tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
 !                                            g / vpt(k,j,i) *                      &
@@ -1240,16 +1240,16 @@
 !                                            ) * dd2zu(k)
 !                         ELSE IF ( cloud_physics )  THEN
-!                            IF ( ql(k,j,i) == 0.0 )  THEN
-!                               k1 = 1.0 + 0.61 * q(k,j,i)
-!                               k2 = 0.61 * pt(k,j,i)
+!                            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)
-!                               k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-!                                          ( q(k,j,i) - ql(k,j,i) ) *          &
-!                                    ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-!                                    ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+!                               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 *          &
 !                                    ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-!                               k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+!                               k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
 !                            ENDIF
 !                            tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
@@ -1259,6 +1259,6 @@
 !                                            ) * dd2zu(k)
 !                         ELSE IF ( cloud_droplets )  THEN
-!                            k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-!                            k2 = 0.61 * pt(k,j,i)
+!                            k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+!                            k2 = 0.61_wp * pt(k,j,i)
 !                            tend(k,j,i) = tend(k,j,i) -                          &
 !                                          kh(k,j,i) * g / vpt(k,j,i) *           &
@@ -1289,23 +1289,23 @@
 !
 !                            IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
-!                               k1 = 1.0 + 0.61 * q(k,j,i)
-!                               k2 = 0.61 * pt(k,j,i)
+!                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+!                               k2 = 0.61_wp * pt(k,j,i)
 !                            ELSE IF ( cloud_physics )  THEN
-!                               IF ( ql(k,j,i) == 0.0 )  THEN
-!                                  k1 = 1.0 + 0.61 * q(k,j,i)
-!                                  k2 = 0.61 * pt(k,j,i)
+!                               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)
-!                                  k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-!                                             ( q(k,j,i) - ql(k,j,i) ) *          &
-!                                       ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-!                                       ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+!                                  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 * &
 !                                       ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-!                                  k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+!                                  k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
 !                               ENDIF
 !                            ELSE IF ( cloud_droplets )  THEN
-!                               k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-!                               k2 = 0.61 * pt(k,j,i)
+!                               k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+!                               k2 = 0.61_wp * pt(k,j,i)
 !                            ENDIF
 !
@@ -1330,23 +1330,23 @@
 !
 !                            IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
-!                               k1 = 1.0 + 0.61 * q(k,j,i)
-!                               k2 = 0.61 * pt(k,j,i)
+!                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+!                               k2 = 0.61_wp * pt(k,j,i)
 !                            ELSE IF ( cloud_physics )  THEN
-!                               IF ( ql(k,j,i) == 0.0 )  THEN
-!                                  k1 = 1.0 + 0.61 * q(k,j,i)
-!                                  k2 = 0.61 * pt(k,j,i)
+!                               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)
-!                                  k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-!                                             ( q(k,j,i) - ql(k,j,i) ) *          &
-!                                       ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-!                                       ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+!                                  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 * &
 !                                       ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-!                                  k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+!                                  k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
 !                               ENDIF
 !                            ELSE IF ( cloud_droplets )  THEN
-!                               k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-!                               k2 = 0.61 * pt(k,j,i)
+!                               k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+!                               k2 = 0.61_wp * pt(k,j,i)
 !                            ENDIF
 !
@@ -1426,27 +1426,27 @@
        DO  k = nzb_diff_s_outer(j,i), nzt
 
-          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
-                + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
-                + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-          IF ( def < 0.0 )  def = 0.0
+          dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+          dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+          dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+          dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+          dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+          dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+          dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+          dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+          dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+          def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
+                + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2    &
+                + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+          IF ( def < 0.0_wp )  def = 0.0_wp
 
           tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -1456,5 +1456,5 @@
        IF ( prandtl_layer )  THEN
 
-          IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) )  THEN
+          IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) )  THEN
 
 !
@@ -1465,12 +1465,12 @@
 
              dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
-             dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                            u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
-             dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-             dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                            v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
+             dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                               u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
+             dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+             dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                               v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
              dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
 
-             IF ( wall_e_y(j,i) /= 0.0 )  THEN
+             IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
 !
 !--             Inconsistency removed: as the thermal stratification 
@@ -1486,21 +1486,21 @@
                 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                     wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
-                km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25 * &
-                             0.5 * dy
-                IF ( km_neutral > 0.0 )  THEN
+                km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25_wp * &
+                             0.5_wp * dy
+                IF ( km_neutral > 0.0_wp )  THEN
                    dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
                    dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
                 ELSE
-                   dudy = 0.0
-                   dwdy = 0.0
+                   dudy = 0.0_wp
+                   dwdy = 0.0_wp
                 ENDIF
              ELSE
-                dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                   u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                   w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
              ENDIF
 
-             IF ( wall_e_x(j,i) /= 0.0 )  THEN
+             IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
 !
 !--             Inconsistency removed: as the thermal stratification 
@@ -1516,25 +1516,25 @@
                 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                     wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
-                km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25 * & 
-                             0.5 * dx
-                IF ( km_neutral > 0.0 )  THEN
+                km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25_wp * & 
+                             0.5_wp * dx
+                IF ( km_neutral > 0.0_wp )  THEN
                    dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
                    dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
                 ELSE
-                   dvdx = 0.0
-                   dwdx = 0.0
+                   dvdx = 0.0_wp
+                   dwdx = 0.0_wp
                 ENDIF
              ELSE
-                dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                   v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                   w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
              ENDIF
 
-             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+             def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                    dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                   dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-             IF ( def < 0.0 )  def = 0.0
+                   dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+             IF ( def < 0.0_wp )  def = 0.0_wp
 
              tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -1549,12 +1549,12 @@
 
                 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
-                dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                               u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-                dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                               v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+                dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                  u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+                dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                  v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
                 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
 
-                IF ( wall_e_y(j,i) /= 0.0 )  THEN
+                IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
 !
 !--                Inconsistency removed: as the thermal stratification 
@@ -1567,20 +1567,20 @@
 !--                      validation has been available
                    km_neutral = kappa * ( usvs(k)**2 + & 
-                                          wsvs(k)**2 )**0.25 * 0.5 * dy
-                   IF ( km_neutral > 0.0 )  THEN
+                                          wsvs(k)**2 )**0.25_wp * 0.5_wp * dy
+                   IF ( km_neutral > 0.0_wp )  THEN
                       dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
                       dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
                    ELSE
-                      dudy = 0.0
-                      dwdy = 0.0
+                      dudy = 0.0_wp
+                      dwdy = 0.0_wp
                    ENDIF
                 ELSE
-                   dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                   u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                   dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                   w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                   dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                      u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                   dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                      w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
                 ENDIF
 
-                IF ( wall_e_x(j,i) /= 0.0 )  THEN
+                IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
 !
 !--                Inconsistency removed: as the thermal stratification 
@@ -1593,24 +1593,24 @@
 !--                      validation has been available
                    km_neutral = kappa * ( vsus(k)**2 + & 
-                                          wsus(k)**2 )**0.25 * 0.5 * dx
-                   IF ( km_neutral > 0.0 )  THEN
+                                          wsus(k)**2 )**0.25_wp * 0.5_wp * dx
+                   IF ( km_neutral > 0.0_wp )  THEN
                       dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
                       dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
                    ELSE
-                      dvdx = 0.0
-                      dwdx = 0.0
+                      dvdx = 0.0_wp
+                      dwdx = 0.0_wp
                    ENDIF
                 ELSE
-                   dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                   v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                   dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                   w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                   dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                      v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                   dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                      w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
                 ENDIF
 
-                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                       dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                IF ( def < 0.0 )  def = 0.0
+                      dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                IF ( def < 0.0_wp )  def = 0.0_wp
 
                 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -1623,21 +1623,21 @@
              k = nzb_diff_s_outer(j,i)-1
 
-             dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-             dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-             dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-             dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-             dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-             dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-             dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-             dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-             dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+             dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+             dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+             dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+             dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+             dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+             dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+             dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+             dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+             dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
 
              def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
@@ -1645,5 +1645,5 @@
                    dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
 
-             IF ( def < 0.0 )  def = 0.0
+             IF ( def < 0.0_wp )  def = 0.0_wp
 
              tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -1657,59 +1657,59 @@
              k = nzb_diff_s_inner(j,i)-1
 
-             dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-             dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+             dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+             dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+             dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                 u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
+
+             dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+             dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+             dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                 v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
+
+             dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+             dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+             dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+             def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
+                   + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
+                   + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+             IF ( def < 0.0_wp )  def = 0.0_wp
+
+             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
+
+          ENDIF
+
+       ELSEIF ( use_surface_fluxes )  THEN
+
+          k = nzb_diff_s_outer(j,i)-1
+
+          dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+          dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
                               u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-             dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                              u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
-
-             dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+          dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+          dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
                               v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-             dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-             dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                              v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
-
-             dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+          dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+          dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+          dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
                               w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-             dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+          dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
                               w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-             dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
-                   + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
-                   + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-             IF ( def < 0.0 )  def = 0.0
-
-             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-
-          ENDIF
-
-       ELSEIF ( use_surface_fluxes )  THEN
-
-          k = nzb_diff_s_outer(j,i)-1
-
-          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+          dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+          def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-          IF ( def < 0.0 )  def = 0.0
+                dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+          IF ( def < 0.0_wp )  def = 0.0_wp
 
           tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -1794,6 +1794,6 @@
 
                 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
-                   k1 = 1.0 + 0.61 * q(k,j,i)
-                   k2 = 0.61 * pt(k,j,i)
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
                    tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *   &
                                          ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
@@ -1801,16 +1801,16 @@
                                          ) * dd2zu(k)
                 ELSE IF ( cloud_physics )  THEN
-                   IF ( ql(k,j,i) == 0.0 )  THEN
-                      k1 = 1.0 + 0.61 * q(k,j,i)
-                      k2 = 0.61 * pt(k,j,i)
+                   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)
-                      k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                                 ( q(k,j,i) - ql(k,j,i) ) *          &
-                           ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                           ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                      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 *          &
                            ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                      k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
                    ENDIF
                    tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *   &
@@ -1819,6 +1819,6 @@
                                          ) * dd2zu(k)
                 ELSE IF ( cloud_droplets )  THEN
-                   k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-                   k2 = 0.61 * pt(k,j,i)
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
                    tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *  &
                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +    &
@@ -1833,23 +1833,23 @@
 
                 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
-                   k1 = 1.0 + 0.61 * q(k,j,i)
-                   k2 = 0.61 * pt(k,j,i)
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
                 ELSE IF ( cloud_physics )  THEN
-                   IF ( ql(k,j,i) == 0.0 )  THEN
-                      k1 = 1.0 + 0.61 * q(k,j,i)
-                      k2 = 0.61 * pt(k,j,i)
+                   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)
-                      k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                                 ( q(k,j,i) - ql(k,j,i) ) *          &
-                           ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                           ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                      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 *          &
                            ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                      k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
                    ENDIF
                 ELSE IF ( cloud_droplets )  THEN
-                   k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-                   k2 = 0.61 * pt(k,j,i)
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
                 ENDIF
 
@@ -1862,23 +1862,23 @@
 
                 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
-                   k1 = 1.0 + 0.61 * q(k,j,i)
-                   k2 = 0.61 * pt(k,j,i)
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
                 ELSE IF ( cloud_physics )  THEN
-                   IF ( ql(k,j,i) == 0.0 )  THEN
-                      k1 = 1.0 + 0.61 * q(k,j,i)
-                      k2 = 0.61 * pt(k,j,i)
+                   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)
-                      k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                                 ( q(k,j,i) - ql(k,j,i) ) *          &
-                           ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                           ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                      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 *          &
                            ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                      k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
                    ENDIF
                 ELSE IF ( cloud_droplets )  THEN
-                   k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-                   k2 = 0.61 * pt(k,j,i)
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
                 ENDIF
 
@@ -1917,6 +1917,6 @@
           IF ( first_call )  THEN
              ALLOCATE( u_0(nysg:nyng,nxlg:nxrg), v_0(nysg:nyng,nxlg:nxrg) )
-             u_0 = 0.0   ! just to avoid access of uninitialized memory
-             v_0 = 0.0   ! within exchange_horiz_2d
+             u_0 = 0.0_wp   ! just to avoid access of uninitialized memory
+             v_0 = 0.0_wp   ! within exchange_horiz_2d
              first_call = .FALSE.
           ENDIF
@@ -1942,10 +1942,10 @@
 
                 u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
-                                 ( 0.5 * ( km(ku,j,i) + km(ku,j,i-1) ) +       &
-                                   1.0E-20 )
+                                 ( 0.5_wp * ( km(ku,j,i) + km(ku,j,i-1) ) +    &
+                                   1.0E-20_wp )
 !                                  ( us(j,i) * kappa * zu(1) )
                 v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
-                                 ( 0.5 * ( km(kv,j,i) + km(kv,j-1,i) ) +       &
-                                   1.0E-20 )
+                                 ( 0.5_wp * ( km(kv,j,i) + km(kv,j-1,i) ) +    &
+                                   1.0E-20_wp )
 !                                  ( us(j,i) * kappa * zu(1) )
 
