Index: /palm/trunk/SOURCE/diffusion_e.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_e.f90	(revision 1339)
+++ /palm/trunk/SOURCE/diffusion_e.f90	(revision 1340)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -148,7 +148,7 @@
                    dvar_dz = atmos_ocean_sign * &
                              ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
-                   IF ( dvar_dz > 0.0 ) THEN
-                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
-                                 SQRT( g / var_reference * dvar_dz ) + 1E-5
+                   IF ( dvar_dz > 0.0_wp ) THEN
+                      l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
+                                    SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
                    ELSE
                       l_stable = l_grid(k)
@@ -176,5 +176,5 @@
                 DO  k = nzb_s_inner(j,i)+1, nzt
 
-                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
+                    dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
                                        e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
 
@@ -220,7 +220,7 @@
                    dvar_dz = atmos_ocean_sign * &
                              ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
-                   IF ( dvar_dz > 0.0 ) THEN
-                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
-                                        SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
+                   IF ( dvar_dz > 0.0_wp ) THEN
+                      l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
+                                           SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
                    ELSE
                       l_stable = l_grid(k)
@@ -248,6 +248,6 @@
                 DO  k = nzb_s_inner(j,i)+1, nzt
 
-                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
-                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
+                    dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
+                                             e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
 
                     tend(k,j,i) = tend(k,j,i)                                  &
@@ -356,7 +356,7 @@
                       dvar_dz = atmos_ocean_sign * &
                                 ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
-                      IF ( dvar_dz > 0.0 ) THEN
-                         l_stable = 0.76 * SQRT( e(k,j,i) ) / &
-                                    SQRT( g / var_reference * dvar_dz ) + 1E-5
+                      IF ( dvar_dz > 0.0_wp ) THEN
+                         l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
+                                       SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
                       ELSE
                          l_stable = l_grid(k)
@@ -377,6 +377,6 @@
 !
 !--                   Calculate the tendency terms
-                      dissipation = ( 0.19 + 0.74 * l / ll ) * &
-                                    e(k,j,i) * SQRT( e(k,j,i) ) / l
+                      dissipation = ( 0.19_wp + 0.74_wp * l / ll ) * &
+                                          e(k,j,i) * SQRT( e(k,j,i) ) / l
 
                       tend(k,j,i) = tend(k,j,i)                                  &
@@ -423,7 +423,7 @@
                       dvar_dz = atmos_ocean_sign * &
                                 ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
-                      IF ( dvar_dz > 0.0 ) THEN
-                         l_stable = 0.76 * SQRT( e(k,j,i) ) / &
-                                           SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
+                      IF ( dvar_dz > 0.0_wp ) THEN
+                         l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
+                                              SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
                       ELSE
                          l_stable = l_grid(k)
@@ -444,6 +444,6 @@
 !
 !--                   Calculate the tendency terms
-                      dissipation = ( 0.19 + 0.74 * l / ll ) * &
-                                    e(k,j,i) * SQRT( e(k,j,i) ) / l
+                      dissipation = ( 0.19_wp + 0.74_wp * l / ll ) * &
+                                          e(k,j,i) * SQRT( e(k,j,i) ) / l
 
                       tend(k,j,i) = tend(k,j,i)                                &
@@ -541,11 +541,11 @@
           dvar_dz = atmos_ocean_sign * &
                     ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
-          IF ( dvar_dz > 0.0 ) THEN
+          IF ( dvar_dz > 0.0_wp ) THEN
              IF ( use_single_reference_value )  THEN
-                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
-                                  SQRT( g / var_reference * dvar_dz ) + 1E-5
+                l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
+                                     SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
              ELSE
-                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
-                                  SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
+                l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
+                                     SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
              ENDIF
           ELSE
@@ -566,6 +566,6 @@
 !
 !--       Calculate the tendency term
-          dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
-                           SQRT( e(k,j,i) ) / l(k)
+          dissipation(k) = ( 0.19_wp + 0.74_wp * l(k) / ll(k) ) * e(k,j,i) * &
+                                 SQRT( e(k,j,i) ) / l(k)
 
           tend(k,j,i) = tend(k,j,i)                                           &
Index: /palm/trunk/SOURCE/diffusion_s.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_s.f90	(revision 1339)
+++ /palm/trunk/SOURCE/diffusion_s.f90	(revision 1340)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -120,34 +120,34 @@
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                                          + 0.5 * (                            &
+                                          + 0.5_wp * (                         &
                         ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
                       - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
-                                                  ) * ddx2                     &
-                                          + 0.5 * (                            &
+                                                     ) * ddx2                  &
+                                          + 0.5_wp * (                         &
                         ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
                       - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
-                                                  ) * ddy2
+                                                     ) * ddy2
              ENDDO
 
 !
 !--          Apply prescribed horizontal wall heatflux where necessary
-             IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0 ) ) &
+             IF ( ( wall_w_x(j,i) .NE. 0.0_wp ) .OR. ( wall_w_y(j,i) .NE. 0.0_wp ) ) &
              THEN
                 DO  k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
 
                    tend(k,j,i) = tend(k,j,i)                                   &
-                                                + ( fwxp(j,i) * 0.5 *          &
-                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
-                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                 &
-                                                   -fwxm(j,i) * 0.5 *          &
+                                                + ( fwxp(j,i) * 0.5_wp *       &
+                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
+                        + ( 1.0_wp - fwxp(j,i) ) * wall_s_flux(1)              &
+                                                   -fwxm(j,i) * 0.5_wp *       &
                         ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
-                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                 &
+                        + ( 1.0_wp - fwxm(j,i) ) * wall_s_flux(2)              &
                                                   ) * ddx2                     &
-                                                + ( fwyp(j,i) * 0.5 *          &
-                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
-                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                 &
-                                                   -fwym(j,i) * 0.5 *          &
+                                                + ( fwyp(j,i) * 0.5_wp *       &
+                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
+                        + ( 1.0_wp - fwyp(j,i) ) * wall_s_flux(3)              &
+                                                   -fwym(j,i) * 0.5_wp *       &
                         ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
-                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                 &
+                        + ( 1.0_wp - fwym(j,i) ) * wall_s_flux(4)              &
                                                   ) * ddy2
                 ENDDO
@@ -161,8 +161,8 @@
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                                       + 0.5 * (                               &
+                                       + 0.5_wp * (                            &
             ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
           - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
-                                               ) * ddzw(k)
+                                                  ) * ddzw(k)
              ENDDO
 
@@ -175,7 +175,7 @@
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                                       + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )     &
-                                               * ( s(k+1,j,i)-s(k,j,i) )       &
-                                               * ddzu(k+1)                     &
+                                       + ( 0.5_wp * ( kh(k,j,i)+kh(k+1,j,i) )  &
+                                                  * ( s(k+1,j,i)-s(k,j,i) )    &
+                                                  * ddzu(k+1)                  &
                                            + s_flux_b(j,i)                     &
                                          ) * ddzw(k)
@@ -192,7 +192,7 @@
                 tend(k,j,i) = tend(k,j,i)                                      &
                                        + ( - s_flux_t(j,i)                     &
-                                           - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )   &
-                                                 * ( s(k,j,i)-s(k-1,j,i) )     &
-                                                 * ddzu(k)                     &
+                                           - 0.5_wp * ( kh(k-1,j,i)+kh(k,j,i) )&
+                                                    * ( s(k,j,i)-s(k-1,j,i) )  &
+                                                    * ddzu(k)                  &
                                          ) * ddzw(k)
 
@@ -251,12 +251,12 @@
 
                    tend(k,j,i) = tend(k,j,i)                                   &
-                                          + 0.5 * (                            &
+                                          + 0.5_wp * (                         &
                         ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
                       - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
-                                                  ) * ddx2                     &
-                                          + 0.5 * (                            &
+                                                     ) * ddx2                  &
+                                          + 0.5_wp * (                         &
                         ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
                       - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
-                                                  ) * ddy2
+                                                     ) * ddy2
                 ENDIF
              ENDDO
@@ -266,20 +266,20 @@
              DO  k = 1, nzt
                 IF ( k > nzb_s_inner(j,i)  .AND.  k <= nzb_s_outer(j,i)  .AND. &
-                     ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 ) )    &
+                     ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp ) )    &
                 THEN
                    tend(k,j,i) = tend(k,j,i)                                   &
-                                                + ( fwxp(j,i) * 0.5 *          &
-                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
-                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                 &
-                                                   -fwxm(j,i) * 0.5 *          &
+                                                + ( fwxp(j,i) * 0.5_wp *       &
+                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
+                        + ( 1.0_wp - fwxp(j,i) ) * wall_s_flux(1)              &
+                                                   -fwxm(j,i) * 0.5_wp *       &
                         ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
-                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                 &
+                        + ( 1.0_wp - fwxm(j,i) ) * wall_s_flux(2)              &
                                                   ) * ddx2                     &
-                                                + ( fwyp(j,i) * 0.5 *          &
-                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
-                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                 &
-                                                   -fwym(j,i) * 0.5 *          &
+                                                + ( fwyp(j,i) * 0.5_wp *       &
+                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
+                        + ( 1.0_wp - fwyp(j,i) ) * wall_s_flux(3)              &
+                                                   -fwym(j,i) * 0.5_wp *       &
                         ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
-                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                 &
+                        + ( 1.0_wp - fwym(j,i) ) * wall_s_flux(4)              &
                                                   ) * ddy2
                 ENDIF
@@ -293,8 +293,8 @@
                 IF ( k >= nzb_diff_s_inner(j,i) )  THEN
                    tend(k,j,i) = tend(k,j,i)                                   &
-                                       + 0.5 * (                               &
+                                       + 0.5_wp * (                            &
             ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
           - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
-                                               ) * ddzw(k)
+                                                  ) * ddzw(k)
                 ENDIF
              ENDDO
@@ -306,7 +306,7 @@
                 IF ( use_surface_fluxes  .AND.  k == nzb_s_inner(j,i)+1 )  THEN
                    tend(k,j,i) = tend(k,j,i)                                   &
-                                          + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )  &
-                                                  * ( s(k+1,j,i)-s(k,j,i) )    &
-                                                  * ddzu(k+1)                  &
+                                          + ( 0.5_wp * ( kh(k,j,i)+kh(k+1,j,i) )&
+                                                     * ( s(k+1,j,i)-s(k,j,i) ) &
+                                                     * ddzu(k+1)               &
                                               + s_flux_b(j,i)                  &
                                             ) * ddzw(k)
@@ -319,7 +319,7 @@
                    tend(k,j,i) = tend(k,j,i)                                   &
                                           + ( - s_flux_t(j,i)                  &
-                                              - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )&
-                                                    * ( s(k,j,i)-s(k-1,j,i) )  &
-                                                    * ddzu(k)                  &
+                                              - 0.5_wp * ( kh(k-1,j,i)+kh(k,j,i) )&
+                                                       * ( s(k,j,i)-s(k-1,j,i) )  &
+                                                       * ddzu(k)                  &
                                             ) * ddzw(k)
                 ENDIF
@@ -372,34 +372,34 @@
 
           tend(k,j,i) = tend(k,j,i)                                            &
-                                          + 0.5 * (                            &
+                                          + 0.5_wp * (                         &
                         ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
                       - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
-                                                  ) * ddx2                     &
-                                          + 0.5 * (                            &
+                                                     ) * ddx2                  &
+                                          + 0.5_wp * (                         &
                         ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
                       - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
-                                                  ) * ddy2
+                                                     ) * ddy2
        ENDDO
 
 !
 !--    Apply prescribed horizontal wall heatflux where necessary
-       IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0 ) )       &
+       IF ( ( wall_w_x(j,i) .NE. 0.0_wp ) .OR. ( wall_w_y(j,i) .NE. 0.0_wp ) )       &
        THEN
           DO  k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
 
              tend(k,j,i) = tend(k,j,i)                                         &
-                                                + ( fwxp(j,i) * 0.5 *          &
-                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
-                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                 &
-                                                   -fwxm(j,i) * 0.5 *          &
+                                                + ( fwxp(j,i) * 0.5_wp *       &
+                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) )  &
+                        + ( 1.0_wp - fwxp(j,i) ) * wall_s_flux(1)              &
+                                                   -fwxm(j,i) * 0.5_wp *       &
                         ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) )  &
-                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                 &
+                        + ( 1.0_wp - fwxm(j,i) ) * wall_s_flux(2)              &
                                                   ) * ddx2                     &
-                                                + ( fwyp(j,i) * 0.5 *          &
-                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
-                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                 &
-                                                   -fwym(j,i) * 0.5 *          &
+                                                + ( fwyp(j,i) * 0.5_wp *       &
+                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) )  &
+                        + ( 1.0_wp - fwyp(j,i) ) * wall_s_flux(3)              &
+                                                   -fwym(j,i) * 0.5_wp *       &
                         ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) )  &
-                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                 &
+                        + ( 1.0_wp - fwym(j,i) ) * wall_s_flux(4)              &
                                                   ) * ddy2
           ENDDO
@@ -413,8 +413,8 @@
 
           tend(k,j,i) = tend(k,j,i)                                            &
-                                       + 0.5 * (                               &
+                                       + 0.5_wp * (                            &
             ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
           - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
-                                               ) * ddzw(k)
+                                                  ) * ddzw(k)
        ENDDO
 
@@ -425,7 +425,7 @@
           k = nzb_s_inner(j,i)+1
 
-          tend(k,j,i) = tend(k,j,i) + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) )        &
-                                            * ( s(k+1,j,i)-s(k,j,i) )          &
-                                            * ddzu(k+1)                        &
+          tend(k,j,i) = tend(k,j,i) + ( 0.5_wp * ( kh(k,j,i)+kh(k+1,j,i) )     &
+                                               * ( s(k+1,j,i)-s(k,j,i) )       &
+                                               * ddzu(k+1)                     &
                                         + s_flux_b(j,i)                        &
                                       ) * ddzw(k)
@@ -440,7 +440,7 @@
 
           tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(j,i)                        &
-                                      - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )        &
-                                            * ( s(k,j,i)-s(k-1,j,i) )          &
-                                            * ddzu(k)                          &
+                                      - 0.5_wp * ( kh(k-1,j,i)+kh(k,j,i) )     &
+                                               * ( s(k,j,i)-s(k-1,j,i) )       &
+                                               * ddzu(k)                       &
                                       ) * ddzw(k)
 
Index: /palm/trunk/SOURCE/diffusion_u.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_u.f90	(revision 1339)
+++ /palm/trunk/SOURCE/diffusion_u.f90	(revision 1340)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -131,14 +131,14 @@
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmyp = 0.25 *                                                  &
+                kmyp = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
-                kmym = 0.25 *                                                  &
+                kmym = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                      & + 2.0 * (                                              &
+                      & + 2.0_wp * (                                           &
                       &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )    &
                       &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )    &
-                      &         ) * ddx2                                       &
+                      &            ) * ddx2                                    &
                       & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy         &
                       &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx         &
@@ -150,17 +150,17 @@
 !
 !--          Wall functions at the north and south walls, respectively
-             IF ( wall_u(j,i) /= 0.0 )  THEN
+             IF ( wall_u(j,i) /= 0.0_wp )  THEN
 
                 DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
-                   kmyp = 0.25 *                                               &
+                   kmyp = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
-                   kmym = 0.25 *                                               &
+                   kmym = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
 
                    tend(k,j,i) = tend(k,j,i)                                   &
-                                 + 2.0 * (                                     &
+                                 + 2.0_wp * (                                  &
                                        km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
                                      - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
-                                         ) * ddx2                              &
+                                            ) * ddx2                           &
                                  + (   fyp(j,i) * (                            &
                                   kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
@@ -182,7 +182,7 @@
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmzp = 0.25 *                                                  &
+                kmzp = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
-                kmzm = 0.25 *                                                  &
+                kmzm = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
 
@@ -211,7 +211,7 @@
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmzp = 0.25 *                                                  &
+                kmzp = 0.25_wp *                                               &
                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
-                kmzm = 0.25 *                                                  &
+                kmzm = 0.25_wp *                                               &
                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
 
@@ -231,7 +231,7 @@
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmzp = 0.25 *                                                  &
+                kmzp = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
-                kmzm = 0.25 *                                                  &
+                kmzm = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
 
@@ -303,14 +303,14 @@
 !
 !--                Interpolate eddy diffusivities on staggered gridpoints
-                   kmyp = 0.25 *                                               &
+                   kmyp = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
-                   kmym = 0.25 *                                               &
+                   kmym = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
 
                    tend(k,j,i) = tend(k,j,i)                                   &
-                         & + 2.0 * (                                           &
+                         & + 2.0_wp * (                                        &
                          &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   ) &
                          &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) ) &
-                         &         ) * ddx2                                    &
+                         &            ) * ddx2                                 &
                          & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy      &
                          &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx      &
@@ -325,16 +325,16 @@
              DO  k = 1, nzt
                 IF( k > nzb_u_inner(j,i)  .AND.  k <= nzb_u_outer(j,i)  .AND.  &
-                    wall_u(j,i) /= 0.0 )  THEN
-
-                   kmyp = 0.25 *                                               &
+                    wall_u(j,i) /= 0.0_wp )  THEN
+
+                   kmyp = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
-                   kmym = 0.25 *                                               &
+                   kmym = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
 
                    tend(k,j,i) = tend(k,j,i)                                   &
-                                 + 2.0 * (                                     &
+                                 + 2.0_wp * (                                  &
                                        km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
                                      - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
-                                         ) * ddx2                              &
+                                            ) * ddx2                           &
                                  + (   fyp(j,i) * (                            &
                                   kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
@@ -357,7 +357,7 @@
 !
 !--                Interpolate eddy diffusivities on staggered gridpoints
-                   kmzp = 0.25 *                                               &
+                   kmzp = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
-                   kmzm = 0.25 *                                               &
+                   kmzm = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
 
@@ -394,7 +394,7 @@
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmzp = 0.25 *                                                  &
+                kmzp = 0.25_wp *                                               &
                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
-                kmzm = 0.25 *                                                  &
+                kmzm = 0.25_wp *                                               &
                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
 
@@ -422,7 +422,7 @@
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmzp = 0.25 *                                                  &
+                kmzp = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
-                kmzm = 0.25 *                                                  &
+                kmzm = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
 
@@ -478,12 +478,12 @@
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
-          kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
-          kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
+          kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
+          kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
 
           tend(k,j,i) = tend(k,j,i)                                            &
-                      & + 2.0 * (                                              &
+                      & + 2.0_wp * (                                           &
                       &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )    &
                       &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )    &
-                      &         ) * ddx2                                       &
+                      &            ) * ddx2                                    &
                       & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy         &
                       &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx         &
@@ -495,5 +495,5 @@
 !
 !--    Wall functions at the north and south walls, respectively
-       IF ( wall_u(j,i) .NE. 0.0 )  THEN
+       IF ( wall_u(j,i) .NE. 0.0_wp )  THEN
 
 !
@@ -503,12 +503,12 @@
 
           DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
-             kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
-             kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
+             kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
+             kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
 
              tend(k,j,i) = tend(k,j,i)                                         &
-                                 + 2.0 * (                                     &
+                                 + 2.0_wp * (                                  &
                                        km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
                                      - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
-                                         ) * ddx2                              &
+                                            ) * ddx2                           &
                                  + (   fyp(j,i) * (                            &
                                   kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
@@ -530,6 +530,6 @@
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
-          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
-          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
+          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
+          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
 
           tend(k,j,i) = tend(k,j,i)                                            &
@@ -557,6 +557,6 @@
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
-          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
-          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
+          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
+          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
 
           tend(k,j,i) = tend(k,j,i)                                            &
@@ -575,7 +575,7 @@
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
-          kmzp = 0.25 *                                                        &
+          kmzp = 0.25_wp *                                                     &
                  ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
-          kmzm = 0.25 *                                                        &
+          kmzm = 0.25_wp *                                                     &
                  ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
 
Index: /palm/trunk/SOURCE/diffusion_v.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_v.f90	(revision 1339)
+++ /palm/trunk/SOURCE/diffusion_v.f90	(revision 1340)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -129,7 +129,7 @@
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmxp = 0.25 * &
+                kmxp = 0.25_wp * &
                        ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
-                kmxm = 0.25 * &
+                kmxm = 0.25_wp * &
                        ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
 
@@ -140,25 +140,25 @@
                       &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
                       &   ) * ddx                                              &
-                      & + 2.0 * (                                              &
+                      & + 2.0_wp * (                                           &
                       &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )      &
                       &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )      &
-                      &         ) * ddy2
+                      &            ) * ddy2
              ENDDO
 
 !
 !--          Wall functions at the left and right walls, respectively
-             IF ( wall_v(j,i) /= 0.0 )  THEN
+             IF ( wall_v(j,i) /= 0.0_wp )  THEN
 
                 DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
-                   kmxp = 0.25 *                                               &
+                   kmxp = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
-                   kmxm = 0.25 *                                               &
+                   kmxm = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
                    
                    tend(k,j,i) = tend(k,j,i)                                   &
-                                 + 2.0 * (                                     &
+                                 + 2.0_wp * (                                  &
                                        km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
                                      - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
-                                         ) * ddy2                              &
+                                            ) * ddy2                           &
                                  + (   fxp(j,i) * (                            &
                                   kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
@@ -180,7 +180,7 @@
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmzp = 0.25 * &
+                kmzp = 0.25_wp * &
                        ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
-                kmzm = 0.25 * &
+                kmzm = 0.25_wp * &
                        ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
 
@@ -209,7 +209,7 @@
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmzp = 0.25 *                                                  &
+                kmzp = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
-                kmzm = 0.25 *                                                  &
+                kmzm = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
 
@@ -229,7 +229,7 @@
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmzp = 0.25 *                                                  &
+                kmzp = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
-                kmzm = 0.25 *                                                  &
+                kmzm = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
 
@@ -301,7 +301,7 @@
 !
 !--                Interpolate eddy diffusivities on staggered gridpoints
-                   kmxp = 0.25 *                                               &
+                   kmxp = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
-                   kmxm = 0.25 *                                               &
+                   kmxm = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
 
@@ -312,8 +312,8 @@
                          &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy          &
                          &   ) * ddx                                           &
-                         & + 2.0 * (                                           &
+                         & + 2.0_wp * (                                        &
                          &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )   &
                          &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )   &
-                         &         ) * ddy2
+                         &            ) * ddy2
                 ENDIF
              ENDDO
@@ -323,16 +323,16 @@
              DO  k = 1, nzt
                 IF( k > nzb_v_inner(j,i)  .AND.  k <= nzb_v_outer(j,i)  .AND.  &
-                    wall_v(j,i) /= 0.0 )  THEN
-
-                   kmxp = 0.25 *                                               &
+                    wall_v(j,i) /= 0.0_wp )  THEN
+
+                   kmxp = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
-                   kmxm = 0.25 *                                               &
+                   kmxm = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
                    
                    tend(k,j,i) = tend(k,j,i)                                   &
-                                 + 2.0 * (                                     &
+                                 + 2.0_wp * (                                  &
                                        km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
                                      - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
-                                         ) * ddy2                              &
+                                            ) * ddy2                           &
                                  + (   fxp(j,i) * (                            &
                                   kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
@@ -355,7 +355,7 @@
 !
 !--                Interpolate eddy diffusivities on staggered gridpoints
-                   kmzp = 0.25 *                                               &
+                   kmzp = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
-                   kmzm = 0.25 *                                               &
+                   kmzm = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
 
@@ -392,7 +392,7 @@
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmzp = 0.25 *                                                  &
+                kmzp = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
-                kmzm = 0.25 *                                                  &
+                kmzm = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
 
@@ -420,7 +420,7 @@
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmzp = 0.25 *                                                  &
+                kmzp = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
-                kmzm = 0.25 *                                                  &
+                kmzm = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
 
@@ -476,6 +476,6 @@
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
-          kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
-          kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
+          kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
+          kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
 
           tend(k,j,i) = tend(k,j,i)                                            &
@@ -485,13 +485,13 @@
                       &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
                       &   ) * ddx                                              &
-                      & + 2.0 * (                                              &
+                      & + 2.0_wp * (                                           &
                       &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )      &
                       &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )      &
-                      &         ) * ddy2
+                      &            ) * ddy2
        ENDDO
 
 !
 !--    Wall functions at the left and right walls, respectively
-       IF ( wall_v(j,i) /= 0.0 )  THEN
+       IF ( wall_v(j,i) /= 0.0_wp )  THEN
 
 !
@@ -501,14 +501,14 @@
 
           DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
-             kmxp = 0.25 *                                                     &
+             kmxp = 0.25_wp *                                                  &
                     ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
-             kmxm = 0.25 *                                                     &
+             kmxm = 0.25_wp *                                                  &
                     ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
 
              tend(k,j,i) = tend(k,j,i)                                         &
-                                 + 2.0 * (                                     &
+                                 + 2.0_wp * (                                  &
                                        km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
                                      - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
-                                         ) * ddy2                              &
+                                            ) * ddy2                           &
                                  + (   fxp(j,i) * (                            &
                                   kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
@@ -530,6 +530,6 @@
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
-          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
-          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
+          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
 
           tend(k,j,i) = tend(k,j,i)                                            &
@@ -557,6 +557,6 @@
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
-          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
-          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
+          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
 
           tend(k,j,i) = tend(k,j,i)                                            &
@@ -575,7 +575,7 @@
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
-          kmzp = 0.25 * &
+          kmzp = 0.25_wp * &
                  ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
-          kmzm = 0.25 * &
+          kmzm = 0.25_wp * &
                  ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
 
Index: /palm/trunk/SOURCE/diffusion_w.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_w.f90	(revision 1339)
+++ /palm/trunk/SOURCE/diffusion_w.f90	(revision 1340)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -158,10 +158,10 @@
                       &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
                       & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
-                      &         ) * ddzu(k+1)
+                      &            ) * ddzu(k+1)
              ENDDO
 
 !
 !--          Wall functions at all vertical walls, where necessary
-             IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
+             IF ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp )  THEN
 
                 DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
@@ -198,8 +198,8 @@
                                      + wall_w_y(j,i) * wsvs(k,j,i)             &
                                    ) * ddy                                     &
-                                 + 2.0 * (                                     &
+                                 + 2.0_wp * (                                  &
                            km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
                          - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
-                                         ) * ddzu(k+1)
+                                            ) * ddzu(k+1)
                 ENDDO
              ENDIF
@@ -288,5 +288,5 @@
                          &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
                          & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
-                         &         ) * ddzu(k+1)
+                         &            ) * ddzu(k+1)
                 ENDIF
              ENDDO
@@ -297,5 +297,5 @@
 
                 IF ( k > nzb_w_inner(j,i)  .AND.  k <= nzb_w_outer(j,i)  .AND. &
-                     wall_w_x(j,i) /= 0.0  .AND.  wall_w_y(j,i) /= 0.0 )  THEN
+                     wall_w_x(j,i) /= 0.0_wp  .AND.  wall_w_y(j,i) /= 0.0_wp )  THEN
 !
 !--                Interpolate eddy diffusivities on staggered gridpoints
@@ -330,8 +330,8 @@
                                      + wall_w_y(j,i) * wsvs(k,j,i)             &
                                    ) * ddy                                     &
-                                 + 2.0 * (                                     &
+                                 + 2.0_wp * (                                  &
                            km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
                          - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
-                                         ) * ddzu(k+1)
+                                            ) * ddzu(k+1)
                 ENDIF
              ENDDO
@@ -400,25 +400,25 @@
                       &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
                       & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
-                      &         ) * ddzu(k+1)
+                      &            ) * ddzu(k+1)
        ENDDO
 
 !
 !--    Wall functions at all vertical walls, where necessary
-       IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
+       IF ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp )  THEN
 
 !
 !--       Calculate the horizontal momentum fluxes w'u' and/or w'v'
-          IF ( wall_w_x(j,i) /= 0.0 )  THEN
+          IF ( wall_w_x(j,i) /= 0.0_wp )  THEN
              CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i),     &
                                wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
           ELSE
-             wsus = 0.0
+             wsus = 0.0_wp
           ENDIF
 
-          IF ( wall_w_y(j,i) /= 0.0 )  THEN
+          IF ( wall_w_y(j,i) /= 0.0_wp )  THEN
              CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i),     &
                                wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
           ELSE
-             wsvs = 0.0
+             wsvs = 0.0_wp
           ENDIF
 
@@ -452,8 +452,8 @@
                                      + wall_w_y(j,i) * wsvs(k)                 &
                                    ) * ddy                                     &
-                                 + 2.0 * (                                     &
+                                 + 2.0_wp * (                                  &
                            km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
                          - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
-                                         ) * ddzu(k+1)
+                                            ) * ddzu(k+1)
           ENDDO
        ENDIF
Index: /palm/trunk/SOURCE/diffusivities.f90
===================================================================
--- /palm/trunk/SOURCE/diffusivities.f90	(revision 1339)
+++ /palm/trunk/SOURCE/diffusivities.f90	(revision 1340)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -99,5 +99,5 @@
 !
 !-- Initialization for calculation of the mixing length profile
-    sums_l_l = 0.0
+    sums_l_l = 0.0_wp
 
 !
@@ -113,5 +113,5 @@
 !
 !-- Introduce an optional minimum tke
-    IF ( e_min > 0.0 )  THEN
+    IF ( e_min > 0.0_wp )  THEN
        !$OMP DO
        !$acc loop
@@ -142,11 +142,11 @@
                 dvar_dz = atmos_ocean_sign * &  ! inverse effect of pt/rho gradient
                           ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
-                IF ( dvar_dz > 0.0 ) THEN
+                IF ( dvar_dz > 0.0_wp ) THEN
                    IF ( use_single_reference_value )  THEN
                       l_stable = 0.76_wp * sqrt_e /                            &
-                                 SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
+                                    SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
                    ELSE
                       l_stable = 0.76_wp * sqrt_e /                            &
-                                 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
+                                    SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
                    ENDIF
                 ELSE
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 1339)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 1340)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -441,5 +441,5 @@
        IF ( humidity_remote )  THEN
           ALLOCATE( qswst_remote(nysg:nyng,nxlg:nxrg))
-          qswst_remote = 0.0
+          qswst_remote = 0.0_wp
        ENDIF
     ENDIF
@@ -453,9 +453,9 @@
     ENDIF
 
-    IF ( dt_dosp /= 9999999.9 )  THEN
+    IF ( dt_dosp /= 9999999.9_wp )  THEN
        ALLOCATE( spectrum_x( 1:nx/2, 1:10, 1:10 ), &
                  spectrum_y( 1:ny/2, 1:10, 1:10 ) )
-       spectrum_x = 0.0
-       spectrum_y = 0.0
+       spectrum_x = 0.0_wp
+       spectrum_y = 0.0_wp
     ENDIF
 
@@ -463,5 +463,5 @@
 !-- 1D-array for large scale subsidence velocity
     ALLOCATE ( w_subs(nzb:nzt+1) )
-    w_subs = 0.0
+    w_subs = 0.0_wp
 
 
@@ -480,5 +480,5 @@
        ENDIF
 
-       IF ( cthf /= 0.0 ) THEN
+       IF ( cthf /= 0.0_wp ) THEN
           ALLOCATE ( lai(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
                      canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
@@ -491,5 +491,5 @@
     IF ( topography /= 'flat' )  THEN
        ALLOCATE( rif_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg,1:4) )
-       rif_wall = 0.0
+       rif_wall = 0.0_wp
     ENDIF
 
@@ -578,6 +578,6 @@
     ALLOCATE( weight_substep(1:intermediate_timestep_count_max), &
               weight_pres(1:intermediate_timestep_count_max) )
-    weight_substep = 1.0
-    weight_pres    = 1.0
+    weight_substep = 1.0_wp
+    weight_pres    = 1.0_wp
     intermediate_timestep_count = 1  ! needed when simulated_time = 0.0
        
@@ -617,6 +617,6 @@
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
-                      qr(:,j,i) = 0.0
-                      nr(:,j,i) = 0.0
+                      qr(:,j,i) = 0.0_wp
+                      nr(:,j,i) = 0.0_wp
                    ENDDO
                 ENDDO
@@ -640,24 +640,24 @@
              IF ( prandtl_layer )  THEN
                 rif  = rif1d(nzb+1)
-                ts   = 0.0  ! could actually be computed more accurately in the
-                            ! 1D model. Update when opportunity arises.
+                ts   = 0.0_wp  ! could actually be computed more accurately in the
+                               ! 1D model. Update when opportunity arises.
                 us   = us1d
                 usws = usws1d
                 vsws = vsws1d
              ELSE
-                ts   = 0.0  ! must be set, because used in
-                rif  = 0.0  ! flowste
-                us   = 0.0
-                usws = 0.0
-                vsws = 0.0
+                ts   = 0.0_wp  ! must be set, because used in
+                rif  = 0.0_wp  ! flowste
+                us   = 0.0_wp
+                usws = 0.0_wp
+                vsws = 0.0_wp
              ENDIF
 
           ELSE
-             e    = 0.0  ! must be set, because used in
-             rif  = 0.0  ! flowste
-             ts   = 0.0
-             us   = 0.0
-             usws = 0.0
-             vsws = 0.0
+             e    = 0.0_wp  ! must be set, because used in
+             rif  = 0.0_wp  ! flowste
+             ts   = 0.0_wp
+             us   = 0.0_wp
+             usws = 0.0_wp
+             vsws = 0.0_wp
           ENDIF
           uswst = top_momentumflux_u
@@ -669,9 +669,9 @@
 !--       Update when opportunity arises!
           IF ( humidity  .OR.  passive_scalar )  THEN
-             qs = 0.0
+             qs = 0.0_wp
              IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
                   precipitation )  THEN
-                qrs = 0.0
-                nrs = 0.0
+                qrs = 0.0_wp
+                nrs = 0.0_wp
              ENDIF
           ENDIF
@@ -682,6 +682,6 @@
              DO  i = nxl-1, nxr+1
                 DO  j = nys-1, nyn+1
-                   u(nzb:nzb_u_inner(j,i),j,i) = 0.0
-                   v(nzb:nzb_v_inner(j,i),j,i) = 0.0
+                   u(nzb:nzb_u_inner(j,i),j,i) = 0.0_wp
+                   v(nzb:nzb_v_inner(j,i),j,i) = 0.0_wp
                 ENDDO
              ENDDO
@@ -745,6 +745,6 @@
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0
-                v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0
+                u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0_wp
+                v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0_wp
              ENDDO
           ENDDO
@@ -764,6 +764,6 @@
                    DO  i = nxlg, nxrg
                       DO  j = nysg, nyng
-                         qr(:,j,i) = 0.0
-                         nr(:,j,i) = 0.0
+                         qr(:,j,i) = 0.0_wp
+                         nr(:,j,i) = 0.0_wp
                       ENDDO
                    ENDDO
@@ -784,8 +784,8 @@
              km   = km_constant
              kh   = km / prandtl_number
-             e    = 0.0
-          ELSEIF ( e_init > 0.0 )  THEN
+             e    = 0.0_wp
+          ELSEIF ( e_init > 0.0_wp )  THEN
              DO  k = nzb+1, nzt
-                km(k,:,:) = 0.1 * l_grid(k) * SQRT( e_init )
+                km(k,:,:) = 0.1_wp * l_grid(k) * SQRT( e_init )
              ENDDO
              km(nzb,:,:)   = km(nzb+1,:,:)
@@ -795,22 +795,22 @@
           ELSE
              IF ( .NOT. ocean )  THEN
-                kh   = 0.01   ! there must exist an initial diffusion, because 
-                km   = 0.01   ! otherwise no TKE would be produced by the
+                kh   = 0.01_wp   ! there must exist an initial diffusion, because 
+                km   = 0.01_wp   ! otherwise no TKE would be produced by the
                               ! production terms, as long as not yet
                               ! e = (u*/cm)**2 at k=nzb+1
              ELSE
-                kh   = 0.00001
-                km   = 0.00001
-             ENDIF
-             e    = 0.0
-          ENDIF
-          rif   = 0.0
-          ts    = 0.0
-          us    = 0.0
-          usws  = 0.0
+                kh   = 0.00001_wp
+                km   = 0.00001_wp
+             ENDIF
+             e    = 0.0_wp
+          ENDIF
+          rif   = 0.0_wp
+          ts    = 0.0_wp
+          us    = 0.0_wp
+          usws  = 0.0_wp
           uswst = top_momentumflux_u
-          vsws  = 0.0
+          vsws  = 0.0_wp
           vswst = top_momentumflux_v
-          IF ( humidity  .OR.  passive_scalar ) qs = 0.0
+          IF ( humidity  .OR.  passive_scalar ) qs = 0.0_wp
 
 !
@@ -829,6 +829,6 @@
 !--    Bottom boundary
        IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2  )  THEN
-          u(nzb,:,:) = 0.0
-          v(nzb,:,:) = 0.0
+          u(nzb,:,:) = 0.0_wp
+          v(nzb,:,:) = 0.0_wp
        ENDIF
 
@@ -836,11 +836,11 @@
 !--    Apply channel flow boundary condition
        IF ( TRIM( bc_uv_t ) == 'dirichlet_0' )  THEN
-          u(nzt+1,:,:) = 0.0
-          v(nzt+1,:,:) = 0.0
+          u(nzt+1,:,:) = 0.0_wp
+          v(nzt+1,:,:) = 0.0_wp
        ENDIF
 
 !
 !--    Calculate virtual potential temperature
-       IF ( humidity ) vpt = pt * ( 1.0 + 0.61 * q )
+       IF ( humidity ) vpt = pt * ( 1.0_wp + 0.61_wp * q )
 
 !
@@ -849,6 +849,6 @@
        hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
        IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2)  THEN
-          hom(nzb,1,5,:) = 0.0   
-          hom(nzb,1,6,:) = 0.0  
+          hom(nzb,1,5,:) = 0.0_wp
+          hom(nzb,1,6,:) = 0.0_wp
        ENDIF
        hom(:,1,7,:)  = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
@@ -919,6 +919,6 @@
              IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
                   precipitation )  THEN
-                qrsws = 0.0
-                nrsws = 0.0
+                qrsws = 0.0_wp
+                nrsws = 0.0_wp
              ENDIF
              IF ( constant_waterflux )  THEN
@@ -954,9 +954,9 @@
 
              IF ( humidity  .OR.  passive_scalar )  THEN
-                qswst = 0.0
+                qswst = 0.0_wp
                 IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
                      precipitation ) THEN
-                   nrswst = 0.0
-                   qrswst = 0.0
+                   nrswst = 0.0_wp
+                   qrswst = 0.0_wp
                 ENDIF
              ENDIF
@@ -971,5 +971,5 @@
 !--       Initialization in case of a coupled model run
           IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
-             tswst = 0.0
+             tswst = 0.0_wp
           ENDIF
 
@@ -990,13 +990,13 @@
 !--          to be zero before the first time step. It approaches its correct
 !--          value in the course of the first few time steps.
-             shf   = 0.0
+             shf   = 0.0_wp
           ENDIF
 
           IF ( humidity  .OR.  passive_scalar )  THEN
-             IF ( .NOT. constant_waterflux )  qsws   = 0.0
+             IF ( .NOT. constant_waterflux )  qsws   = 0.0_wp
              IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
                   precipitation )  THEN
-                qrsws = 0.0
-                nrsws = 0.0
+                qrsws = 0.0_wp
+                nrsws = 0.0_wp
              ENDIF
           ENDIF
@@ -1023,13 +1023,13 @@
 !
 !--    For the moment, vertical velocity is zero
-       w = 0.0
+       w = 0.0_wp
 
 !
 !--    Initialize array sums (must be defined in first call of pres)
-       sums = 0.0
+       sums = 0.0_wp
 
 !
 !--    In case of iterative solvers, p must get an initial value
-       IF ( psolver == 'multigrid'  .OR.  psolver == 'sor' )  p = 0.0
+       IF ( psolver == 'multigrid'  .OR.  psolver == 'sor' )  p = 0.0_wp
 
 !
@@ -1037,8 +1037,8 @@
 !--    are zero at beginning of the simulation
        IF ( cloud_physics )  THEN
-          ql = 0.0
-          IF ( precipitation )  precipitation_amount = 0.0
+          ql = 0.0_wp
+          IF ( precipitation )  precipitation_amount = 0.0_wp
           IF ( icloud_scheme == 0 )  THEN
-             qc = 0.0
+             qc = 0.0_wp
              nc_1d = nc_const
           ENDIF 
@@ -1058,5 +1058,5 @@
 !
 !--    If required, change the surface temperature at the start of the 3D run
-       IF ( pt_surface_initial_change /= 0.0 )  THEN
+       IF ( pt_surface_initial_change /= 0.0_wp )  THEN
           pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
        ENDIF
@@ -1066,5 +1066,5 @@
 !--    run
        IF ( ( humidity .OR. passive_scalar ) .AND. &
-            q_surface_initial_change /= 0.0 )  THEN
+            q_surface_initial_change /= 0.0_wp )  THEN
           q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
        ENDIF
@@ -1075,15 +1075,15 @@
 !
 !--    Initialize old and new time levels.
-       te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
+       te_m = 0.0_wp; tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
        e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
 
        IF ( humidity  .OR.  passive_scalar )  THEN
-          tq_m = 0.0
+          tq_m = 0.0_wp
           q_p = q
           IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
                precipitation )  THEN
-             tqr_m = 0.0
+             tqr_m = 0.0_wp
              qr_p = qr
-             tnr_m = 0.0
+             tnr_m = 0.0_wp
              nr_p = nr
           ENDIF
@@ -1091,5 +1091,5 @@
 
        IF ( ocean )  THEN
-          tsa_m = 0.0
+          tsa_m = 0.0_wp
           sa_p  = sa
        ENDIF
@@ -1171,5 +1171,5 @@
                    u(k,j,nxlg:-1)  = mean_inflow_profiles(k,1)
                    v(k,j,nxlg:-1)  = mean_inflow_profiles(k,2)
-                   w(k,j,nxlg:-1)  = 0.0
+                   w(k,j,nxlg:-1)  = 0.0_wp
                    pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
                    e(k,j,nxlg:-1)  = mean_inflow_profiles(k,5)
@@ -1183,10 +1183,10 @@
 !--       vertically because otherwise the turbulent inflow layer will grow
 !--       in time.
-          IF ( inflow_damping_height == 9999999.9 )  THEN
+          IF ( inflow_damping_height == 9999999.9_wp )  THEN
 !
 !--          Default: use the inversion height calculated by the prerun; if
 !--          this is zero, inflow_damping_height must be explicitly
 !--          specified.
-             IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0 )  THEN
+             IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0_wp )  THEN
                 inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
              ELSE
@@ -1199,9 +1199,9 @@
           ENDIF
 
-          IF ( inflow_damping_width == 9999999.9 )  THEN
+          IF ( inflow_damping_width == 9999999.9_wp )  THEN
 !
 !--          Default for the transition range: one tenth of the undamped
 !--          layer
-             inflow_damping_width = 0.1 * inflow_damping_height
+             inflow_damping_width = 0.1_wp * inflow_damping_height
 
           ENDIF
@@ -1212,11 +1212,11 @@
 
              IF ( zu(k) <= inflow_damping_height )  THEN
-                inflow_damping_factor(k) = 1.0
+                inflow_damping_factor(k) = 1.0_wp
              ELSEIF ( zu(k) <= ( inflow_damping_height + inflow_damping_width ) )  THEN
-                inflow_damping_factor(k) = 1.0 -                               &
+                inflow_damping_factor(k) = 1.0_wp -                            &
                                            ( zu(k) - inflow_damping_height ) / &
                                            inflow_damping_width
              ELSE
-                inflow_damping_factor(k) = 0.0
+                inflow_damping_factor(k) = 0.0_wp
              ENDIF
 
@@ -1235,13 +1235,13 @@
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                u  (nzb:nzb_u_inner(j,i),j,i)   = 0.0
-                v  (nzb:nzb_v_inner(j,i),j,i)   = 0.0
-                w  (nzb:nzb_w_inner(j,i),j,i)   = 0.0
-                e  (nzb:nzb_w_inner(j,i),j,i)   = 0.0
-                tu_m(nzb:nzb_u_inner(j,i),j,i)  = 0.0
-                tv_m(nzb:nzb_v_inner(j,i),j,i)  = 0.0
-                tw_m(nzb:nzb_w_inner(j,i),j,i)  = 0.0
-                te_m(nzb:nzb_w_inner(j,i),j,i)  = 0.0
-                tpt_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
+                u  (nzb:nzb_u_inner(j,i),j,i)   = 0.0_wp
+                v  (nzb:nzb_v_inner(j,i),j,i)   = 0.0_wp
+                w  (nzb:nzb_w_inner(j,i),j,i)   = 0.0_wp
+                e  (nzb:nzb_w_inner(j,i),j,i)   = 0.0_wp
+                tu_m(nzb:nzb_u_inner(j,i),j,i)  = 0.0_wp
+                tv_m(nzb:nzb_v_inner(j,i),j,i)  = 0.0_wp
+                tw_m(nzb:nzb_w_inner(j,i),j,i)  = 0.0_wp
+                te_m(nzb:nzb_w_inner(j,i),j,i)  = 0.0_wp
+                tpt_m(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp
              ENDDO
           ENDDO
@@ -1272,14 +1272,14 @@
 !--    have to be predefined here because they are used (but multiplied with 0)
 !--    there before they are set. 
-       te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
+       te_m = 0.0_wp; tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
        IF ( humidity  .OR.  passive_scalar )  THEN
-          tq_m = 0.0
+          tq_m = 0.0_wp
           IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
                precipitation )  THEN
-             tqr_m = 0.0
-             tnr_m = 0.0
-          ENDIF
-       ENDIF
-       IF ( ocean )  tsa_m = 0.0
+             tqr_m = 0.0_wp
+             tnr_m = 0.0_wp
+          ENDIF
+       ENDIF
+       IF ( ocean )  tsa_m = 0.0_wp
 
     ELSE
@@ -1323,6 +1323,6 @@
        IF ( use_prescribed_profile_data )  THEN
 
-          volume_flow_initial_l = 0.0
-          volume_flow_area_l    = 0.0
+          volume_flow_initial_l = 0.0_wp
+          volume_flow_area_l    = 0.0_wp
 
           IF ( nxr == nx )  THEN
@@ -1359,6 +1359,6 @@
        ELSEIF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
 
-          volume_flow_initial_l = 0.0
-          volume_flow_area_l    = 0.0
+          volume_flow_initial_l = 0.0_wp
+          volume_flow_area_l    = 0.0_wp
 
           IF ( nxr == nx )  THEN
@@ -1395,6 +1395,6 @@
        ELSEIF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
 
-          volume_flow_initial_l = 0.0
-          volume_flow_area_l    = 0.0
+          volume_flow_initial_l = 0.0_wp
+          volume_flow_area_l    = 0.0_wp
 
           IF ( nxr == nx )  THEN
@@ -1513,5 +1513,5 @@
           DO  j = nys, nyn
              DO  k = nzb, nzt+1 
-                IF ( lad_s(k,j,i) > 0.0 )  THEN
+                IF ( lad_s(k,j,i) > 0.0_wp )  THEN
                    lad_u(k,j,i)   = lad_s(k,j,i) 
                    lad_u(k,j,i+1) = lad_s(k,j,i)
@@ -1521,10 +1521,10 @@
              ENDDO
              DO  k = nzb, nzt
-                lad_w(k,j,i) = 0.5 * ( lad_s(k+1,j,i) + lad_s(k,j,i) )
+                lad_w(k,j,i) = 0.5_wp * ( lad_s(k+1,j,i) + lad_s(k,j,i) )
              ENDDO
           ENDDO
        ENDDO
 
-       lad_w(pch_index,:,:) = 0.0
+       lad_w(pch_index,:,:) = 0.0_wp
        lad_w(nzt+1,:,:)     = lad_w(nzt,:,:)
 
@@ -1535,17 +1535,17 @@
 !
 !--    Initialisation of the canopy heat source distribution
-       IF ( cthf /= 0.0 )  THEN
+       IF ( cthf /= 0.0_wp )  THEN
 !
 !--       Piecewise evaluation of the leaf area index by 
 !--       integration of the leaf area density
-          lai(:,:,:) = 0.0
+          lai(:,:,:) = 0.0_wp
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
                 DO  k = pch_index-1, 0, -1
                    lai(k,j,i) = lai(k+1,j,i) +                   &
-                                ( 0.5 * ( lad_w(k+1,j,i) +       &
+                                ( 0.5_wp * ( lad_w(k+1,j,i) +    &
                                           lad_s(k+1,j,i) ) *     &
                                   ( zw(k+1) - zu(k+1) ) )  +     &
-                                ( 0.5 * ( lad_w(k,j,i)   +       &
+                                ( 0.5_wp * ( lad_w(k,j,i)   +    &
                                           lad_s(k+1,j,i) ) *     &
                                   ( zu(k+1) - zw(k) ) )
@@ -1561,5 +1561,5 @@
                 DO  k = 0, pch_index
                    canopy_heat_flux(k,j,i) = cthf *                    &
-                                             exp( -0.6 * lai(k,j,i) )
+                                             exp( -0.6_wp * lai(k,j,i) )
                 ENDDO
              ENDDO
@@ -1577,5 +1577,5 @@
 !
 !-- If required, initialize dvrp-software
-    IF ( dt_dvrp /= 9999999.9 )  CALL init_dvrp
+    IF ( dt_dvrp /= 9999999.9_wp )  CALL init_dvrp
 
     IF ( ocean )  THEN
@@ -1624,6 +1624,6 @@
     ELSE                                     ! for Euler-method
 
-       weight_substep(1) = 1.0      
-       weight_pres(1)    = 1.0                   
+       weight_substep(1) = 1.0_wp      
+       weight_pres(1)    = 1.0_wp                   
 
     ENDIF
@@ -1631,13 +1631,13 @@
 !
 !-- Initialize Rayleigh damping factors
-    rdf    = 0.0
-    rdf_sc = 0.0
-    IF ( rayleigh_damping_factor /= 0.0 )  THEN
+    rdf    = 0.0_wp
+    rdf_sc = 0.0_wp
+    IF ( rayleigh_damping_factor /= 0.0_wp )  THEN
        IF ( .NOT. ocean )  THEN
           DO  k = nzb+1, nzt
              IF ( zu(k) >= rayleigh_damping_height )  THEN
                 rdf(k) = rayleigh_damping_factor * &
-                      ( SIN( pi * 0.5 * ( zu(k) - rayleigh_damping_height )    &
-                                      / ( zu(nzt) - rayleigh_damping_height ) )&
+                      ( SIN( pi * 0.5_wp * ( zu(k) - rayleigh_damping_height ) &
+                                         / ( zu(nzt) - rayleigh_damping_height ) )&
                       )**2
              ENDIF
@@ -1647,6 +1647,6 @@
              IF ( zu(k) <= rayleigh_damping_height )  THEN
                 rdf(k) = rayleigh_damping_factor * &
-                      ( SIN( pi * 0.5 * ( rayleigh_damping_height - zu(k) )    &
-                                      / ( rayleigh_damping_height - zu(nzb+1)))&
+                      ( SIN( pi * 0.5_wp * ( rayleigh_damping_height - zu(k) ) &
+                                         / ( rayleigh_damping_height - zu(nzb+1)))&
                       )**2
              ENDIF
@@ -1659,5 +1659,5 @@
 !-- Initialize the starting level and the vertical smoothing factor used for 
 !-- the external pressure gradient
-    dp_smooth_factor = 1.0
+    dp_smooth_factor = 1.0_wp
     IF ( dp_external )  THEN
 !
@@ -1670,9 +1670,9 @@
        ENDIF
        IF ( dp_smooth )  THEN
-          dp_smooth_factor(:dp_level_ind_b) = 0.0
+          dp_smooth_factor(:dp_level_ind_b) = 0.0_wp
           DO  k = dp_level_ind_b+1, nzt
-             dp_smooth_factor(k) = 0.5 * ( 1.0 + SIN( pi * &
-                  ( REAL( k - dp_level_ind_b, KIND=wp ) /  &
-                    REAL( nzt - dp_level_ind_b, KIND=wp ) - 0.5 ) ) )
+             dp_smooth_factor(k) = 0.5_wp * ( 1.0_wp + SIN( pi *               &
+                        ( REAL( k - dp_level_ind_b, KIND=wp ) /                &
+                          REAL( nzt - dp_level_ind_b, KIND=wp ) - 0.5_wp ) ) )
           ENDDO
        ENDIF
@@ -1683,12 +1683,12 @@
 !-- non-cyclic lateral boundaries. The damping zone has the maximum value 
 !-- at the inflow boundary and decreases to zero at pt_damping_width.
-    ptdf_x = 0.0
-    ptdf_y = 0.0
+    ptdf_x = 0.0_wp
+    ptdf_y = 0.0_wp
     IF ( bc_lr_dirrad )  THEN
        DO  i = nxl, nxr
           IF ( ( i * dx ) < pt_damping_width )  THEN
-             ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5 *              &
-                         REAL( pt_damping_width - i * dx, KIND=wp ) / (                 &
-                         REAL( pt_damping_width, KIND=wp )            ) ) )**2 
+             ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5_wp *              &
+                            REAL( pt_damping_width - i * dx, KIND=wp ) / (     &
+                            REAL( pt_damping_width, KIND=wp )            ) ) )**2 
           ENDIF
        ENDDO
@@ -1697,7 +1697,7 @@
           IF ( ( i * dx ) > ( nx * dx - pt_damping_width ) )  THEN
              ptdf_x(i) = pt_damping_factor *                                   &
-                         SIN( pi * 0.5 *                                       &
-                              ( ( i - nx ) * dx + pt_damping_width ) /         &
-                              REAL( pt_damping_width, KIND=wp ) )**2
+                         SIN( pi * 0.5_wp *                                    &
+                                 ( ( i - nx ) * dx + pt_damping_width ) /      &
+                                 REAL( pt_damping_width, KIND=wp ) )**2
           ENDIF
        ENDDO 
@@ -1706,7 +1706,7 @@
           IF ( ( j * dy ) > ( ny * dy - pt_damping_width ) )  THEN
              ptdf_y(j) = pt_damping_factor *                                   &
-                         SIN( pi * 0.5 *                                       &
-                              ( ( j - ny ) * dy + pt_damping_width ) /         &
-                              REAL( pt_damping_width, KIND=wp ) )**2
+                         SIN( pi * 0.5_wp *                                    &
+                                 ( ( j - ny ) * dy + pt_damping_width ) /      &
+                                 REAL( pt_damping_width, KIND=wp ) )**2
           ENDIF
        ENDDO 
@@ -1715,7 +1715,7 @@
           IF ( ( j * dy ) < pt_damping_width )  THEN
              ptdf_y(j) = pt_damping_factor *                                   &
-                         SIN( pi * 0.5 *                                       &
-                             ( pt_damping_width - j * dy ) /                   &
-                             REAL( pt_damping_width, KIND=wp ) )**2
+                         SIN( pi * 0.5_wp *                                    &
+                                ( pt_damping_width - j * dy ) /                &
+                                REAL( pt_damping_width, KIND=wp ) )**2
           ENDIF
        ENDDO
@@ -1727,9 +1727,9 @@
 !-- are called from flow_statistics (or - depending on the chosen model run -
 !-- are never initialized)
-    sums_divnew_l      = 0.0
-    sums_divold_l      = 0.0
-    sums_l_l           = 0.0
-    sums_up_fraction_l = 0.0
-    sums_wsts_bc_l     = 0.0
+    sums_divnew_l      = 0.0_wp
+    sums_divold_l      = 0.0_wp
+    sums_l_l           = 0.0_wp
+    sums_up_fraction_l = 0.0_wp
+    sums_wsts_bc_l     = 0.0_wp
 
 !
@@ -1737,7 +1737,7 @@
 !-- Ghost points are excluded because counting values at the ghost boundaries
 !-- would bias the statistics
-    rmask = 1.0
-    rmask(:,nxlg:nxl-1,:) = 0.0;  rmask(:,nxr+1:nxrg,:) = 0.0
-    rmask(nysg:nys-1,:,:) = 0.0;  rmask(nyn+1:nyng,:,:) = 0.0
+    rmask = 1.0_wp
+    rmask(:,nxlg:nxl-1,:) = 0.0_wp;  rmask(:,nxr+1:nxrg,:) = 0.0_wp
+    rmask(nysg:nys-1,:,:) = 0.0_wp;  rmask(nyn+1:nyng,:,:) = 0.0_wp
 
 !
@@ -1769,5 +1769,5 @@
     ngp_2dh_l         = 0
     ngp_2dh           = 0
-    ngp_3d_inner_l    = 0.0
+    ngp_3d_inner_l    = 0.0_wp
     ngp_3d_inner      = 0
     ngp_3d            = 0
@@ -1777,5 +1777,5 @@
        DO  i = nxl, nxr
           DO  j = nys, nyn
-             IF ( rmask(j,i,sr) == 1.0 )  THEN
+             IF ( rmask(j,i,sr) == 1.0_wp )  THEN
 !
 !--             All xy-grid points
Index: /palm/trunk/SOURCE/package_parin.f90
===================================================================
--- /palm/trunk/SOURCE/package_parin.f90	(revision 1339)
+++ /palm/trunk/SOURCE/package_parin.f90	(revision 1340)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! -----------------
+! REAL constants defined as wp-kinds
 !
 ! Former revisions:
@@ -196,5 +197,5 @@
 !-- Default setting of dt_dosp here (instead of check_parameters), because its
 !-- current value is needed in init_pegrid
-    IF ( dt_dosp == 9999999.9 )  dt_dosp = dt_data_output
+    IF ( dt_dosp == 9999999.9_wp )  dt_dosp = dt_data_output
 
  30 CONTINUE
Index: /palm/trunk/SOURCE/plant_canopy_model.f90
===================================================================
--- /palm/trunk/SOURCE/plant_canopy_model.f90	(revision 1339)
+++ /palm/trunk/SOURCE/plant_canopy_model.f90	(revision 1340)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -100,10 +100,10 @@
                                               v(k,j+1,i)      +  &
                                               v(k,j+1,i-1) )     &
-                                            / 4.0 )**2        +  &
+                                            / 4.0_wp )**2     +  &
                                           ( ( w(k-1,j,i-1)    +  &
                                               w(k-1,j,i)      +  &
                                               w(k,j,i-1)      +  &
                                               w(k,j,i) )         &
-                                            / 4.0 )**2 )      *  &
+                                            / 4.0_wp )**2 )   *  &
                                     u(k,j,i)
                    ENDDO
@@ -123,5 +123,5 @@
                                               u(k,j,i)        +  &
                                               u(k,j,i+1) )       &
-                                            / 4.0 )**2        +  &
+                                            / 4.0_wp )**2     +  &
                                               v(k,j,i)**2     +  &
                                           ( ( w(k-1,j-1,i)    +  &
@@ -129,5 +129,5 @@
                                               w(k,j-1,i)      +  &
                                               w(k,j,i) )         &
-                                            / 4.0 )**2 )      *  &
+                                            / 4.0_wp )**2 )   *  &
                                     v(k,j,i)  
                    ENDDO
@@ -147,10 +147,10 @@
                                               u(k+1,j,i)      +  &
                                               u(k+1,j,i+1) )     &
-                                            / 4.0 )**2        +  &
+                                            / 4.0_wp )**2     +  &
                                           ( ( v(k,j,i)        +  &
                                               v(k,j+1,i)      +  &
                                               v(k+1,j,i)      +  &
                                               v(k+1,j+1,i) )     &
-                                            / 4.0 )**2        +  &
+                                            / 4.0_wp )**2     +  &
                                               w(k,j,i)**2 )   *  &
                                     w(k,j,i)
@@ -183,11 +183,11 @@
                                     SQRT( ( ( u(k,j,i)        +       &
                                               u(k,j,i+1) )            &
-                                            / 2.0 )**2        +       &
+                                            / 2.0_wp )**2     +       &
                                           ( ( v(k,j,i)        +       &
                                               v(k,j+1,i) )            &
-                                            / 2.0 )**2        +       &
+                                            / 2.0_wp )**2     +       &
                                           ( ( w(k-1,j,i)      +       & 
                                               w(k,j,i) )              &
-                                            / 2.0 )**2 )      *       &
+                                            / 2.0_wp )**2 )   *       &
                                     ( q(k,j,i) - sls(k,j,i) )
                    ENDDO
@@ -202,14 +202,14 @@
                    DO  k = nzb_s_inner(j,i)+1, pch_index
                       tend(k,j,i) = tend(k,j,i) -                     &
-                                    2.0 * cdc(k,j,i) * lad_s(k,j,i) * &
+                                    2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * &
                                     SQRT( ( ( u(k,j,i)              + &
                                               u(k,j,i+1) )            &
-                                            / 2.0 )**2              + &
+                                            / 2.0_wp )**2           + &
                                           ( ( v(k,j,i)              + &
                                               v(k,j+1,i) )            &
-                                            / 2.0 )**2              + &
+                                            / 2.0_wp )**2           + &
                                           ( ( w(k,j,i)              + &
                                               w(k+1,j,i) )            &
-                                            / 2.0 )**2 )            * &
+                                            / 2.0_wp )**2 )         * &
                                     e(k,j,i)
                    ENDDO
@@ -267,10 +267,10 @@
                                         v(k,j+1,i)  +        &
                                         v(k,j+1,i-1) )       &
-                                      / 4.0 )**2    +        &
+                                      / 4.0_wp )**2 +        &
                                     ( ( w(k-1,j,i-1) +       &
                                         w(k-1,j,i)   +       &
                                         w(k,j,i-1)   +       &
                                         w(k,j,i) )           &
-                                      / 4.0 )**2 ) *         &
+                                      / 4.0_wp )**2 ) *      &
                               u(k,j,i)
           ENDDO
@@ -286,5 +286,5 @@
                                         u(k,j,i)     +       &
                                         u(k,j,i+1) )         &
-                                      / 4.0 )**2     +       &
+                                      / 4.0_wp )**2  +       &
                                         v(k,j,i)**2  +       &
                                     ( ( w(k-1,j-1,i) +       &
@@ -292,5 +292,5 @@
                                         w(k,j-1,i)   +       &
                                         w(k,j,i) )           &
-                                      / 4.0 )**2 ) *         &
+                                      / 4.0_wp )**2 ) *      &
                               v(k,j,i)
           ENDDO
@@ -306,10 +306,10 @@
                                         u(k+1,j,i)  +        &
                                         u(k+1,j,i+1) )       &
-                                      / 4.0 )**2    +        &
+                                      / 4.0_wp )**2 +        &
                                     ( ( v(k,j,i)    +        &
                                         v(k,j+1,i)  +        &
                                         v(k+1,j,i)  +        &
                                         v(k+1,j+1,i) )       &
-                                      / 4.0 )**2    +        &
+                                      / 4.0_wp )**2 +        &
                                         w(k,j,i)**2 ) *      &
                               w(k,j,i)
@@ -336,11 +336,11 @@
                               SQRT( ( ( u(k,j,i)        +       &
                                         u(k,j,i+1) )            &
-                                      / 2.0 )**2        +       &
+                                      / 2.0_wp )**2     +       &
                                     ( ( v(k,j,i)        +       &
                                         v(k,j+1,i) )            &
-                                      / 2.0 )**2        +       &
+                                      / 2.0_wp )**2     +       &
                                     ( ( w(k-1,j,i)      +       &
                                         w(k,j,i) )              &
-                                      / 2.0 )**2 )      *       &
+                                      / 2.0_wp )**2 )   *       &
                               ( q(k,j,i) - sls(k,j,i) )
              ENDDO   
@@ -351,14 +351,14 @@
           DO  k = nzb_s_inner(j,i)+1, pch_index   
              tend(k,j,i) = tend(k,j,i) -                        &
-                              2.0 * cdc(k,j,i) * lad_s(k,j,i) * &
+                              2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * &
                               SQRT( ( ( u(k,j,i)           +    &
                                         u(k,j,i+1) )            &
-                                      / 2.0 )**2           +    &  
+                                      / 2.0_wp )**2        +    &  
                                     ( ( v(k,j,i)           +    &
                                         v(k,j+1,i) )            &
-                                      / 2.0 )**2           +    &
+                                      / 2.0_wp )**2        +    &
                                     ( ( w(k,j,i)           +    &
                                         w(k+1,j,i) )            &
-                                      / 2.0 )**2 )         *    &
+                                      / 2.0_wp )**2 )      *    &
                               e(k,j,i)
           ENDDO
Index: /palm/trunk/SOURCE/prandtl_fluxes.f90
===================================================================
--- /palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 1339)
+++ /palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 1340)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -105,10 +105,10 @@
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
-             ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30 )
+             ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30_wp )
 !
 !--          ts must be limited, because otherwise overflow may occur in case of
 !--          us=0 when computing rif further below
-             IF ( ts(j,i) < -1.05E5 )  ts(j,i) = -1.0E5
-             IF ( ts(j,i) >   1.0E5 )  ts(j,i) =  1.0E5
+             IF ( ts(j,i) < -1.05E5_wp )  ts(j,i) = -1.0E5_wp
+             IF ( ts(j,i) >   1.0E5_wp )  ts(j,i) =  1.0E5_wp
           ENDDO
        ENDDO
@@ -131,20 +131,20 @@
              z_p = zu(k+1) - zw(k)
 
-             IF ( rif(j,i) >= 0.0 )  THEN
+             IF ( rif(j,i) >= 0.0_wp )  THEN
 !
 !--             Stable stratification
                 ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (           &
                                   LOG( z_p / z0h(j,i) ) +                   &
-                                  5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
+                                  5.0_wp * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
                                                                 )
              ELSE
 !
 !--             Unstable stratification
-                a = SQRT( 1.0 - 16.0 * rif(j,i) )
-                b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p )
+                a = SQRT( 1.0_wp - 16.0_wp * rif(j,i) )
+                b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p )
 
                 ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) /  (          &
                           LOG( z_p / z0h(j,i) ) -                           &
-                          2.0 * LOG( ( 1.0 + a ) / ( 1.0 + b ) ) )
+                          2.0_wp * LOG( ( 1.0_wp + a ) / ( 1.0_wp + b ) ) )
              ENDIF
 
@@ -163,5 +163,5 @@
              z_p = zu(k+1) - zw(k)
              rif(j,i) = z_p * kappa * g * ts(j,i) / &
-                        ( pt(k+1,j,i) * ( us(j,i)**2 + 1E-30 ) )
+                        ( pt(k+1,j,i) * ( us(j,i)**2 + 1E-30_wp ) )
 !
 !--          Limit the value range of the Richardson numbers.
@@ -182,6 +182,6 @@
              z_p = zu(k+1) - zw(k)
              rif(j,i) = z_p * kappa * g *                            &
-                        ( ts(j,i) + 0.61 * pt(k+1,j,i) * qs(j,i) ) / &
-                        ( vpt(k+1,j,i) * ( us(j,i)**2 + 1E-30 ) )
+                        ( ts(j,i) + 0.61_wp * pt(k+1,j,i) * qs(j,i) ) / &
+                        ( vpt(k+1,j,i) * ( us(j,i)**2 + 1E-30_wp ) )
 !
 !--          Limit the value range of the Richardson numbers.
@@ -209,28 +209,28 @@
 !--       Compute the absolute value of the horizontal velocity 
 !--       (relative to the surface)
-          uv_total = SQRT( ( 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1)        &
-                                   - u(k,j,i)   - u(k,j,i+1) ) )**2 + &
-                           ( 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i)        &
-                                   - v(k,j,i)   - v(k,j+1,i) ) )**2 )    
-
-
-          IF ( rif(j,i) >= 0.0 )  THEN
+          uv_total = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)        &
+                                      - u(k,j,i)   - u(k,j,i+1) ) )**2 + &
+                           ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)        &
+                                      - v(k,j,i)   - v(k,j+1,i) ) )**2 )    
+
+
+          IF ( rif(j,i) >= 0.0_wp )  THEN
 !
 !--          Stable stratification
              us(j,i) = kappa * uv_total / (                                &
                                   LOG( z_p / z0(j,i) ) +                   &
-                                  5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
+                                  5.0_wp * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
                                           )
           ELSE
 !
 !--          Unstable stratification
-             a = SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
-             b = SQRT( SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) )
+             a = SQRT( SQRT( 1.0_wp - 16.0_wp * rif(j,i) ) )
+             b = SQRT( SQRT( 1.0_wp - 16.0_wp * rif(j,i) / z_p * z0(j,i) ) )
 
              us(j,i) = kappa * uv_total / (                                  &
                        LOG( z_p / z0(j,i) ) -                                &
-                       LOG( ( 1.0 + a )**2 * ( 1.0 + a**2 ) / (              &
-                            ( 1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +          &
-                            2.0 * ( ATAN( a ) - ATAN( b ) )                  &
+                       LOG( ( 1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (              &
+                            ( 1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +          &
+                               2.0_wp * ( ATAN( a ) - ATAN( b ) )                  &
                                            )
           ENDIF
@@ -258,26 +258,26 @@
 !
 !--       Compute Richardson-flux number for this point
-          rifm = 0.5 * ( rif(j,i-1) + rif(j,i) )
-          IF ( rifm >= 0.0 )  THEN
+          rifm = 0.5_wp * ( rif(j,i-1) + rif(j,i) )
+          IF ( rifm >= 0.0_wp )  THEN
 !
 !--          Stable stratification
              usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) )/ (              &
                                      LOG( z_p / z0(j,i) ) +               &
-                                     5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
-                                              )
+                                     5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p &
+                                                            )
           ELSE
 !
 !--          Unstable stratification
-             a = SQRT( SQRT( 1.0 - 16.0 * rifm ) )
-             b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
+             a = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm ) )
+             b = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm / z_p * z0(j,i) ) )
 
              usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) ) / (            &
                          LOG( z_p / z0(j,i) ) -                           &
-                         LOG( (1.0 + a )**2 * ( 1.0 + a**2 ) / (          &
-                              (1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +      &
-                              2.0 * ( ATAN( a ) - ATAN( b ) )             &
+                         LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (          &
+                              (1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +      &
+                                 2.0_wp * ( ATAN( a ) - ATAN( b ) )             &
                                                  )
           ENDIF
-          usws(j,i) = -usws(j,i) * 0.5 * ( us(j,i-1) + us(j,i) )
+          usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
        ENDDO
     ENDDO
@@ -296,26 +296,26 @@
 !
 !--       Compute Richardson-flux number for this point
-          rifm = 0.5 * ( rif(j-1,i) + rif(j,i) )
-          IF ( rifm >= 0.0 )  THEN
+          rifm = 0.5_wp * ( rif(j-1,i) + rif(j,i) )
+          IF ( rifm >= 0.0_wp )  THEN
 !
 !--          Stable stratification
              vsws(j,i) = kappa * ( v(k+1,j,i) -  v(k,j,i) ) / (           &
                                      LOG( z_p / z0(j,i) ) +               &
-                                     5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
-                                              )
+                                     5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p &
+                                                              )
           ELSE
 !
 !--          Unstable stratification
-             a = SQRT( SQRT( 1.0 - 16.0 * rifm ) )
-             b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
+             a = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm ) )
+             b = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm / z_p * z0(j,i) ) )
 
              vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) ) / (            &
                          LOG( z_p / z0(j,i) ) -                           &
-                         LOG( (1.0 + a )**2 * ( 1.0 + a**2 ) / (          &
-                              (1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +      &
-                              2.0 * ( ATAN( a ) - ATAN( b ) )             &
+                         LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (          &
+                              (1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +      &
+                                 2.0_wp * ( ATAN( a ) - ATAN( b ) )             &
                                                  )
           ENDIF
-          vsws(j,i) = -vsws(j,i) * 0.5 * ( us(j-1,i) + us(j,i) )
+          vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j-1,i) + us(j,i) )
        ENDDO
     ENDDO
@@ -331,5 +331,5 @@
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30 )
+                qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp )
              ENDDO
           ENDDO
@@ -355,24 +355,24 @@
 !--             in case of precursor runs)
                 IF ( coupled_run )  THEN
-                   e_q = 6.1 * &
-                        EXP( 0.07 * ( MIN(pt(0,j,i),pt(1,j,i)) - 273.15 ) )
-                   q(k,j,i) = 0.622 * e_q / ( surface_pressure - e_q )
+                   e_q = 6.1_wp * &
+                           EXP( 0.07_wp * ( MIN(pt(0,j,i),pt(1,j,i)) - 273.15_wp ) )
+                   q(k,j,i) = 0.622_wp * e_q / ( surface_pressure - e_q )
                 ENDIF
-                IF ( rif(j,i) >= 0.0 )  THEN
+                IF ( rif(j,i) >= 0.0_wp )  THEN
 !
 !--                Stable stratification
                    qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (          &
                                   LOG( z_p / z0h(j,i) ) +                   &
-                                  5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
+                                  5.0_wp * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
                                                                  )
                 ELSE
 !
 !--                Unstable stratification
-                   a = SQRT( 1.0 - 16.0 * rif(j,i) ) 
-                   b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p ) 
+                   a = SQRT( 1.0_wp - 16.0_wp * rif(j,i) ) 
+                   b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p ) 
  
                    qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) /   (        &
                              LOG( z_p / z0h(j,i) ) -                        &
-                              2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
+                              2.0_wp * LOG( (1.0_wp + a ) / ( 1.0_wp + b ) ) )
                 ENDIF
 
@@ -429,8 +429,8 @@
           !$acc loop independent
           DO  j = nysg, nyng
-             e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1 )**2
+             e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1_wp )**2
 !
 !--          As a test: cm = 0.4
-!            e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4 )**2
+!            e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4_wp )**2
              e(nzb_s_inner(j,i),j,i)   = e(nzb_s_inner(j,i)+1,j,i)
           ENDDO
