Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 1321)
+++ /palm/trunk/SOURCE/Makefile	(revision 1322)
@@ -269,5 +269,5 @@
 diffusivities.o: modules.o mod_kinds.o
 disturb_field.o: modules.o cpulog.o mod_kinds.o random_function.o
-disturb_heatflux.o: modules.o mod_kinds.o
+disturb_heatflux.o: modules.o cpulog.o mod_kinds.o
 eqn_state_seawater.o: modules.o mod_kinds.o
 exchange_horiz.o: modules.o cpulog.o mod_kinds.o
Index: /palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- /palm/trunk/SOURCE/advec_ws.f90	(revision 1321)
+++ /palm/trunk/SOURCE/advec_ws.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -195,10 +195,10 @@
 !
 !--    Set the appropriate factors for scalar and momentum advection.
-       adv_sca_5 = 1./60.
-       adv_sca_3 = 1./12.
-       adv_sca_1 = 1./2.
-       adv_mom_5 = 1./120.
-       adv_mom_3 = 1./24.
-       adv_mom_1 = 1./4.
+       adv_sca_5 = 1.0_wp /  60.0_wp
+       adv_sca_3 = 1.0_wp /  12.0_wp
+       adv_sca_1 = 1.0_wp /   2.0_wp
+       adv_mom_5 = 1.0_wp / 120.0_wp
+       adv_mom_3 = 1.0_wp /  24.0_wp
+       adv_mom_1 = 1.0_wp /   4.0_wp
 !         
 !--    Arrays needed for statical evaluation of fluxes.
@@ -1174,8 +1174,8 @@
                               + ( flux_r(k) *                                 &
                                 ( u_comp(k) - 2.0 * hom(k,1,1,0) )            &
-                              / ( u_comp(k) - gu + 1.0E-20      )             &
+                              / ( u_comp(k) - gu + 1.0E-20_wp   )             &
                               +   diss_r(k) *                                 &
                                   ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )       &
-                              / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )         &
+                              / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) )      &
                               *   weight_substep(intermediate_timestep_count)
 !
@@ -1278,8 +1278,8 @@
                               + ( flux_r(k) *                                 &
                                 ( u_comp(k) - 2.0 * hom(k,1,1,0) )            &
-                              / ( u_comp(k) - gu + 1.0E-20      )             &
+                              / ( u_comp(k) - gu + 1.0E-20_wp   )             &
                               +   diss_r(k) *                                 &
                                   ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )       &
-                              / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )         &
+                              / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) )      &
                               *   weight_substep(intermediate_timestep_count)
 !
@@ -1631,8 +1631,8 @@
              + ( flux_n(k)                                                    &
              * ( v_comp(k) - 2.0 * hom(k,1,2,0) )                             &
-             / ( v_comp(k) - gv + 1.0E-20 )                                   &
+             / ( v_comp(k) - gv + 1.0E-20_wp )                                &
              +   diss_n(k)                                                    &
              *   ABS( v_comp(k) - 2.0 * hom(k,1,2,0) )                        &
-             / ( ABS( v_comp(k) - gv ) +1.0E-20 ) )                           &
+             / ( ABS( v_comp(k) - gv ) +1.0E-20_wp ) )                        &
              *   weight_substep(intermediate_timestep_count)
 !
@@ -1740,8 +1740,8 @@
              + ( flux_n(k)                                                    &
              * ( v_comp(k) - 2.0 * hom(k,1,2,0) )                             &
-             / ( v_comp(k) - gv + 1.0E-20 )                                   &
+             / ( v_comp(k) - gv + 1.0E-20_wp )                                &
              +   diss_n(k)                                                    &
              *   ABS( v_comp(k) - 2.0 * hom(k,1,2,0) )                        &
-             / ( ABS( v_comp(k) - gv ) +1.0E-20 ) )                           &
+             / ( ABS( v_comp(k) - gv ) +1.0E-20_wp ) )                        &
              *   weight_substep(intermediate_timestep_count)
 !
@@ -3312,8 +3312,8 @@
                               + ( flux_r(k) *                                 &
                                 ( u_comp(k) - 2.0 * hom(k,1,1,0) )            &
-                              / ( u_comp(k) - gu + 1.0E-20      )             &
+                              / ( u_comp(k) - gu + 1.0E-20_wp    )            &
                               +   diss_r(k) *                                 &
                                   ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )       &
-                              / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )         &
+                              / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) )      &
                               *   weight_substep(intermediate_timestep_count)
 !
@@ -3418,8 +3418,8 @@
                               + ( flux_r(k) *                                 &
                                 ( u_comp(k) - 2.0 * hom(k,1,1,0) )            &
-                              / ( u_comp(k) - gu + 1.0E-20      )             &
+                              / ( u_comp(k) - gu + 1.0E-20_wp   )             &
                               +   diss_r(k) *                                 &
                                   ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )       &
-                              / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )         &
+                              / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) )      &
                               *   weight_substep(intermediate_timestep_count)
 !
@@ -3749,8 +3749,8 @@
 !                              + ( flux_r    *                                 &
 !                                ( u_comp    - 2.0 * hom(k,1,1,0) )            &
-!                              / ( u_comp    - gu + 1.0E-20      )             &
+!                              / ( u_comp    - gu + 1.0E-20_wp   )             &
 !                              +   diss_r    *                                 &
 !                                  ABS( u_comp    - 2.0 * hom(k,1,1,0) )       &
-!                              / ( ABS( u_comp    - gu ) + 1.0E-20 ) )         &
+!                              / ( ABS( u_comp    - gu ) + 1.0E-20_wp ) )      &
 !                              *   weight_substep(intermediate_timestep_count)
 !
@@ -4108,8 +4108,8 @@
                       + ( flux_n(k)                                           &
                       * ( v_comp(k) - 2.0 * hom(k,1,2,0) )                    &
-                      / ( v_comp(k) - gv + 1.0E-20 )                          &
+                      / ( v_comp(k) - gv + 1.0E-20_wp )                       &
                       +   diss_n(k)                                           &
                       *   ABS( v_comp(k) - 2.0 * hom(k,1,2,0) )               &
-                      / ( ABS( v_comp(k) - gv ) +1.0E-20 ) )                  &
+                      / ( ABS( v_comp(k) - gv ) +1.0E-20_wp ) )               &
                       *   weight_substep(intermediate_timestep_count)
 !
@@ -4221,8 +4221,8 @@
                          + ( flux_n(k)                                        &
                          * ( v_comp(k) - 2.0 * hom(k,1,2,0) )                 &
-                         / ( v_comp(k) - gv + 1.0E-20 )                       &
+                         / ( v_comp(k) - gv + 1.0E-20_wp )                    &
                          +   diss_n(k)                                        &
                          *   ABS( v_comp(k) - 2.0 * hom(k,1,2,0) )            &
-                         / ( ABS( v_comp(k) - gv ) +1.0E-20 ) )               &
+                         / ( ABS( v_comp(k) - gv ) +1.0E-20_wp ) )            &
                          *   weight_substep(intermediate_timestep_count)
 !
@@ -4554,8 +4554,8 @@
 !                      + ( flux_n                                           &
 !                      * ( v_comp - 2.0 * hom(k,1,2,0) )                    &
-!                      / ( v_comp - gv + 1.0E-20 )                          &
+!                      / ( v_comp - gv + 1.0E-20_wp )                       &
 !                      +   diss_n                                           &
 !                      *   ABS( v_comp - 2.0 * hom(k,1,2,0) )               &
-!                      / ( ABS( v_comp - gv ) +1.0E-20 ) )                  &
+!                      / ( ABS( v_comp - gv ) +1.0E-20_wp ) )               &
 !                      *   weight_substep(intermediate_timestep_count)
 !
Index: /palm/trunk/SOURCE/average_3d_data.f90
===================================================================
--- /palm/trunk/SOURCE/average_3d_data.f90	(revision 1321)
+++ /palm/trunk/SOURCE/average_3d_data.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! REAL functions provided with KIND-attribute
 !
 ! Former revisions:
@@ -99,5 +99,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      e_av(k,j,i) = e_av(k,j,i) / REAL( average_count_3d )
+                      e_av(k,j,i) = e_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -107,5 +107,5 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   qsws_av(j,i) = qsws_av(j,i) / REAL( average_count_3d )
+                   qsws_av(j,i) = qsws_av(j,i) / REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
@@ -115,5 +115,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      lpt_av(k,j,i) = lpt_av(k,j,i) / REAL( average_count_3d )
+                      lpt_av(k,j,i) = lpt_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -123,5 +123,5 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   lwp_av(j,i) = lwp_av(j,i) / REAL( average_count_3d )
+                   lwp_av(j,i) = lwp_av(j,i) / REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
@@ -131,5 +131,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      nr_av(k,j,i) = nr_av(k,j,i) / REAL( average_count_3d )
+                      nr_av(k,j,i) = nr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -140,5 +140,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      p_av(k,j,i) = p_av(k,j,i) / REAL( average_count_3d )
+                      p_av(k,j,i) = p_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -149,5 +149,5 @@
                 DO  j = nys, nyn
                    DO  k = nzb, nzt+1
-                      pc_av(k,j,i) = pc_av(k,j,i) / REAL( average_count_3d )
+                      pc_av(k,j,i) = pc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -158,5 +158,5 @@
                 DO  j = nys, nyn
                    DO  k = nzb, nzt+1
-                      pr_av(k,j,i) = pr_av(k,j,i) / REAL( average_count_3d )
+                      pr_av(k,j,i) = pr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -167,5 +167,5 @@
                 DO  j = nysg, nyng
                    precipitation_rate_av(j,i) = precipitation_rate_av(j,i) /   &
-                                                REAL( average_count_3d )
+                                                REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
@@ -175,5 +175,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      pt_av(k,j,i) = pt_av(k,j,i) / REAL( average_count_3d )
+                      pt_av(k,j,i) = pt_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -184,5 +184,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      q_av(k,j,i) = q_av(k,j,i) / REAL( average_count_3d )
+                      q_av(k,j,i) = q_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -193,5 +193,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      qc_av(k,j,i) = qc_av(k,j,i) / REAL( average_count_3d )
+                      qc_av(k,j,i) = qc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -202,5 +202,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      ql_av(k,j,i) = ql_av(k,j,i) / REAL( average_count_3d )
+                      ql_av(k,j,i) = ql_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -211,5 +211,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      ql_c_av(k,j,i) = ql_c_av(k,j,i) / REAL( average_count_3d )
+                      ql_c_av(k,j,i) = ql_c_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -220,5 +220,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      ql_v_av(k,j,i) = ql_v_av(k,j,i) / REAL( average_count_3d )
+                      ql_v_av(k,j,i) = ql_v_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -230,5 +230,5 @@
                    DO  k = nzb, nzt+1
                       ql_vp_av(k,j,i) = ql_vp_av(k,j,i) /                      &
-                                        REAL( average_count_3d )
+                                        REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -239,5 +239,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      qr_av(k,j,i) = qr_av(k,j,i) / REAL( average_count_3d )
+                      qr_av(k,j,i) = qr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -248,5 +248,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      qv_av(k,j,i) = qv_av(k,j,i) / REAL( average_count_3d )
+                      qv_av(k,j,i) = qv_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -257,5 +257,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      rho_av(k,j,i) = rho_av(k,j,i) / REAL( average_count_3d )
+                      rho_av(k,j,i) = rho_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -266,5 +266,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      s_av(k,j,i) = s_av(k,j,i) / REAL( average_count_3d )
+                      s_av(k,j,i) = s_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -275,5 +275,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      sa_av(k,j,i) = sa_av(k,j,i) / REAL( average_count_3d )
+                      sa_av(k,j,i) = sa_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -283,5 +283,5 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   shf_av(j,i) = shf_av(j,i) / REAL( average_count_3d )
+                   shf_av(j,i) = shf_av(j,i) / REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
@@ -290,5 +290,5 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   ts_av(j,i) = ts_av(j,i) / REAL( average_count_3d )
+                   ts_av(j,i) = ts_av(j,i) / REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
@@ -298,5 +298,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      u_av(k,j,i) = u_av(k,j,i) / REAL( average_count_3d )
+                      u_av(k,j,i) = u_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -306,5 +306,5 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   us_av(j,i) = us_av(j,i) / REAL( average_count_3d )
+                   us_av(j,i) = us_av(j,i) / REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
@@ -314,5 +314,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      v_av(k,j,i) = v_av(k,j,i) / REAL( average_count_3d )
+                      v_av(k,j,i) = v_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -323,5 +323,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      vpt_av(k,j,i) = vpt_av(k,j,i) / REAL( average_count_3d )
+                      vpt_av(k,j,i) = vpt_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -332,5 +332,5 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      w_av(k,j,i) = w_av(k,j,i) / REAL( average_count_3d )
+                      w_av(k,j,i) = w_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -340,5 +340,5 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   z0_av(j,i) = z0_av(j,i) / REAL( average_count_3d )
+                   z0_av(j,i) = z0_av(j,i) / REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
@@ -347,5 +347,5 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   z0h_av(j,i) = z0h_av(j,i) / REAL( average_count_3d )
+                   z0h_av(j,i) = z0h_av(j,i) / REAL( average_count_3d, KIND=wp )
                 ENDDO
              ENDDO
Index: /palm/trunk/SOURCE/calc_radiation.f90
===================================================================
--- /palm/trunk/SOURCE/calc_radiation.f90	(revision 1321)
+++ /palm/trunk/SOURCE/calc_radiation.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! exponent 4.0 changed to integer
 !
 ! Former revisions:
@@ -127,5 +127,5 @@
 
              temperature     = pt(nzb,j,i) * t_d_pt(nzb) + l_d_cp * ql(nzb,j,i)
-             blackbody_emission(nzb) = sigma * temperature**4.0
+             blackbody_emission(nzb) = sigma * temperature**4
 
              DO  k = nzb_2d(j,i)+1, nzt
@@ -139,5 +139,5 @@
 
                 temperature     = pt(k,j,i) * t_d_pt(k) + l_d_cp * ql(k,j,i)
-                blackbody_emission(k) = sigma * temperature**4.0
+                blackbody_emission(k) = sigma * temperature**4
 
              ENDDO
@@ -149,5 +149,5 @@
              temperature       = pt(nzt+1,j,i) * t_d_pt(nzt+1) + l_d_cp *      &
                                  ql(nzt+1,j,i)
-             blackbody_emission(nzt+1) = sigma * temperature**4.0
+             blackbody_emission(nzt+1) = sigma * temperature**4
 
 !
@@ -267,5 +267,5 @@
 
        temperature     = pt(nzb,j,i) * t_d_pt(nzb) + l_d_cp * ql(nzb,j,i)
-       blackbody_emission(nzb) = sigma * temperature**4.0
+       blackbody_emission(nzb) = sigma * temperature**4
 
        DO  k = nzb_2d(j,i)+1, nzt
@@ -277,5 +277,5 @@
 
           temperature     = pt(k,j,i) * t_d_pt(k) + l_d_cp * ql(k,j,i)
-          blackbody_emission(k) = sigma * temperature**4.0
+          blackbody_emission(k) = sigma * temperature**4
 
        ENDDO
@@ -286,5 +286,5 @@
        temperature       = pt(nzt+1,j,i) * t_d_pt(nzt+1) + l_d_cp *            &
                            ql(nzt+1,j,i)
-       blackbody_emission(nzt+1) = sigma * temperature**4.0
+       blackbody_emission(nzt+1) = sigma * temperature**4
 
 !
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 1321)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! some REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -259,5 +259,5 @@
     IF ( coupling_mode /= 'uncoupled')  THEN
 
-       IF ( dt_coupling == 9999999.9 )  THEN
+       IF ( dt_coupling == 9999999.9_wp )  THEN
           message_string = 'dt_coupling is not set but required for coup' // &
                            'ling mode "' //  TRIM( coupling_mode ) // '"'
@@ -920,5 +920,5 @@
 !
 !--    In case of no given gradients for ug, choose a zero gradient
-       IF ( ug_vertical_gradient_level(1) == -9999999.9 )  THEN
+       IF ( ug_vertical_gradient_level(1) == -9999999.9_wp )  THEN
           ug_vertical_gradient_level(1) = 0.0
        ENDIF  
@@ -984,5 +984,5 @@
 !
 !--    In case of no given gradients for vg, choose a zero gradient
-       IF ( vg_vertical_gradient_level(1) == -9999999.9 )  THEN
+       IF ( vg_vertical_gradient_level(1) == -9999999.9_wp )  THEN
           vg_vertical_gradient_level(1) = 0.0
        ENDIF
@@ -991,5 +991,5 @@
 !--    Let the initial wind profiles be the calculated ug/vg profiles or
 !--    interpolate them from wind profile data (if given)
-       IF ( u_profile(1) == 9999999.9  .AND.  v_profile(1) == 9999999.9 )  THEN
+       IF ( u_profile(1) == 9999999.9_wp  .AND.  v_profile(1) == 9999999.9_wp )  THEN
 
           u_init = ug
@@ -1018,5 +1018,5 @@
              ENDIF
 
-             IF ( kk < 100 .AND. uv_heights(kk+1) /= 9999999.9 )  THEN
+             IF ( kk < 100 .AND. uv_heights(kk+1) /= 9999999.9_wp )  THEN
                 u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
                                        ( uv_heights(kk+1) - uv_heights(kk) ) * &
@@ -1100,5 +1100,5 @@
 !--    In case of no given temperature gradients, choose gradient of neutral
 !--    stratification
-       IF ( pt_vertical_gradient_level(1) == -9999999.9 )  THEN
+       IF ( pt_vertical_gradient_level(1) == -9999999.9_wp )  THEN
           pt_vertical_gradient_level(1) = 0.0
        ENDIF
@@ -1236,5 +1236,5 @@
 !--       In case of no given leaf area density gradients, choose a vanishing
 !--       gradient
-          IF ( lad_vertical_gradient_level(1) == -9999999.9 ) THEN
+          IF ( lad_vertical_gradient_level(1) == -9999999.9_wp ) THEN
              lad_vertical_gradient_level(1) = 0.0
           ENDIF
@@ -1247,5 +1247,5 @@
 !-- Initialize large scale subsidence if required
     If ( large_scale_subsidence )  THEN
-       IF ( subs_vertical_gradient_level(1) /= -9999999.9 .AND. &
+       IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp .AND. &
                                      .NOT. large_scale_forcing )  THEN
           CALL init_w_subsidence
@@ -1255,5 +1255,5 @@
 !--    are read in from file LSF_DATA
 
-       IF ( subs_vertical_gradient_level(1) == -9999999.9 .AND. &
+       IF ( subs_vertical_gradient_level(1) == -9999999.9_wp .AND. &
                                      .NOT. large_scale_forcing )  THEN
           message_string = 'There is no default large scale vertical ' // &
@@ -1264,5 +1264,5 @@
        ENDIF
     ELSE
-        IF ( subs_vertical_gradient_level(1) /= -9999999.9 )  THEN
+        IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp )  THEN
            message_string = 'Enable usage of large scale subsidence by ' // &
                             'setting large_scale_subsidence = .T..'
@@ -1284,5 +1284,5 @@
     ELSEIF ( reference_state == 'single_value' )  THEN
        use_single_reference_value = .TRUE.
-       IF ( pt_reference == 9999999.9 )  pt_reference = pt_surface
+       IF ( pt_reference == 9999999.9_wp )  pt_reference = pt_surface
        vpt_reference = pt_reference * ( 1.0 + 0.61 * q_surface )
     ELSE
@@ -1528,5 +1528,5 @@
     ENDIF
 
-    IF ( surface_heatflux == 9999999.9  )  THEN
+    IF ( surface_heatflux == 9999999.9_wp  )  THEN
        constant_heatflux = .FALSE.
        IF ( large_scale_forcing )  THEN
@@ -1548,9 +1548,9 @@
     ENDIF
 
-    IF ( top_heatflux     == 9999999.9 )  constant_top_heatflux = .FALSE.
+    IF ( top_heatflux     == 9999999.9_wp )  constant_top_heatflux = .FALSE.
 
     IF ( neutral )  THEN
 
-       IF ( surface_heatflux /= 0.0  .AND.  surface_heatflux /= 9999999.9 ) &
+       IF ( surface_heatflux /= 0.0  .AND.  surface_heatflux /= 9999999.9_wp ) &
        THEN
           message_string = 'heatflux must not be set for pure neutral flow'
@@ -1558,5 +1558,5 @@
        ENDIF
 
-       IF ( top_heatflux /= 0.0  .AND.  top_heatflux /= 9999999.9 ) &
+       IF ( top_heatflux /= 0.0  .AND.  top_heatflux /= 9999999.9_wp ) &
        THEN
           message_string = 'heatflux must not be set for pure neutral flow'
@@ -1566,9 +1566,9 @@
     ENDIF
 
-    IF ( top_momentumflux_u /= 9999999.9  .AND.  &
-         top_momentumflux_v /= 9999999.9 )  THEN
+    IF ( top_momentumflux_u /= 9999999.9_wp  .AND.  &
+         top_momentumflux_v /= 9999999.9_wp )  THEN
        constant_top_momentumflux = .TRUE.
-    ELSEIF (  .NOT. ( top_momentumflux_u == 9999999.9  .AND.  &
-           top_momentumflux_v == 9999999.9 ) )  THEN
+    ELSEIF (  .NOT. ( top_momentumflux_u == 9999999.9_wp  .AND.  &
+           top_momentumflux_v == 9999999.9_wp ) )  THEN
        message_string = 'both, top_momentumflux_u AND top_momentumflux_v ' // &
                         'must be set'
@@ -1617,6 +1617,6 @@
        ENDIF
 
-       IF ( top_salinityflux == 9999999.9 )  constant_top_salinityflux = .FALSE.
-       IF ( ibc_sa_t == 1  .AND.   top_salinityflux == 9999999.9 )  THEN
+       IF ( top_salinityflux == 9999999.9_wp )  constant_top_salinityflux = .FALSE.
+       IF ( ibc_sa_t == 1  .AND.   top_salinityflux == 9999999.9_wp )  THEN
           message_string = 'boundary condition: bc_sa_t = "' // &
                            TRIM( bc_sa_t ) // '" requires to set ' // &
@@ -1667,5 +1667,5 @@
        ENDIF
 
-       IF ( surface_waterflux == 9999999.9  )  THEN
+       IF ( surface_waterflux == 9999999.9_wp  )  THEN
           constant_waterflux = .FALSE.
           IF ( large_scale_forcing )  THEN
@@ -1817,14 +1817,14 @@
 !-- Set the default intervals for data output, if necessary
 !-- NOTE: dt_dosp has already been set in package_parin
-    IF ( dt_data_output /= 9999999.9 )  THEN
-       IF ( dt_dopr           == 9999999.9 )  dt_dopr           = dt_data_output
-       IF ( dt_dopts          == 9999999.9 )  dt_dopts          = dt_data_output
-       IF ( dt_do2d_xy        == 9999999.9 )  dt_do2d_xy        = dt_data_output
-       IF ( dt_do2d_xz        == 9999999.9 )  dt_do2d_xz        = dt_data_output
-       IF ( dt_do2d_yz        == 9999999.9 )  dt_do2d_yz        = dt_data_output
-       IF ( dt_do3d           == 9999999.9 )  dt_do3d           = dt_data_output
-       IF ( dt_data_output_av == 9999999.9 )  dt_data_output_av = dt_data_output
+    IF ( dt_data_output /= 9999999.9_wp )  THEN
+       IF ( dt_dopr           == 9999999.9_wp )  dt_dopr           = dt_data_output
+       IF ( dt_dopts          == 9999999.9_wp )  dt_dopts          = dt_data_output
+       IF ( dt_do2d_xy        == 9999999.9_wp )  dt_do2d_xy        = dt_data_output
+       IF ( dt_do2d_xz        == 9999999.9_wp )  dt_do2d_xz        = dt_data_output
+       IF ( dt_do2d_yz        == 9999999.9_wp )  dt_do2d_yz        = dt_data_output
+       IF ( dt_do3d           == 9999999.9_wp )  dt_do3d           = dt_data_output
+       IF ( dt_data_output_av == 9999999.9_wp )  dt_data_output_av = dt_data_output
        DO  mid = 1, max_masks
-          IF ( dt_domask(mid) == 9999999.9 )  dt_domask(mid)    = dt_data_output
+          IF ( dt_domask(mid) == 9999999.9_wp )  dt_domask(mid)    = dt_data_output
        ENDDO
     ENDIF
@@ -1832,20 +1832,20 @@
 !
 !-- Set the default skip time intervals for data output, if necessary
-    IF ( skip_time_dopr    == 9999999.9 ) &
+    IF ( skip_time_dopr    == 9999999.9_wp ) &
                                        skip_time_dopr    = skip_time_data_output
-    IF ( skip_time_dosp    == 9999999.9 ) &
+    IF ( skip_time_dosp    == 9999999.9_wp ) &
                                        skip_time_dosp    = skip_time_data_output
-    IF ( skip_time_do2d_xy == 9999999.9 ) &
+    IF ( skip_time_do2d_xy == 9999999.9_wp ) &
                                        skip_time_do2d_xy = skip_time_data_output
-    IF ( skip_time_do2d_xz == 9999999.9 ) &
+    IF ( skip_time_do2d_xz == 9999999.9_wp ) &
                                        skip_time_do2d_xz = skip_time_data_output
-    IF ( skip_time_do2d_yz == 9999999.9 ) &
+    IF ( skip_time_do2d_yz == 9999999.9_wp ) &
                                        skip_time_do2d_yz = skip_time_data_output
-    IF ( skip_time_do3d    == 9999999.9 ) &
+    IF ( skip_time_do3d    == 9999999.9_wp ) &
                                        skip_time_do3d    = skip_time_data_output
-    IF ( skip_time_data_output_av == 9999999.9 ) &
+    IF ( skip_time_data_output_av == 9999999.9_wp ) &
                                 skip_time_data_output_av = skip_time_data_output
     DO  mid = 1, max_masks
-       IF ( skip_time_domask(mid) == 9999999.9 ) &
+       IF ( skip_time_domask(mid) == 9999999.9_wp ) &
                                 skip_time_domask(mid)    = skip_time_data_output
     ENDDO
@@ -1860,5 +1860,5 @@
     ENDIF
 
-    IF ( averaging_interval_pr == 9999999.9 )  THEN
+    IF ( averaging_interval_pr == 9999999.9_wp )  THEN
        averaging_interval_pr = averaging_interval
     ENDIF
@@ -1870,5 +1870,5 @@
     ENDIF
 
-    IF ( averaging_interval_sp == 9999999.9 )  THEN
+    IF ( averaging_interval_sp == 9999999.9_wp )  THEN
        averaging_interval_sp = averaging_interval
     ENDIF
@@ -1882,5 +1882,5 @@
 !
 !-- Set the default interval for profiles entering the temporal average
-    IF ( dt_averaging_input_pr == 9999999.9 )  THEN
+    IF ( dt_averaging_input_pr == 9999999.9_wp )  THEN
        dt_averaging_input_pr = dt_averaging_input
     ENDIF
@@ -1889,5 +1889,5 @@
 !-- Set the default interval for the output of timeseries to a reasonable
 !-- value (tries to minimize the number of calls of flow_statistics)
-    IF ( dt_dots == 9999999.9 )  THEN
+    IF ( dt_dots == 9999999.9_wp )  THEN
        IF ( averaging_interval_pr == 0.0 )  THEN
           dt_dots = MIN( dt_run_control, dt_dopr )
@@ -1916,5 +1916,5 @@
 !-- Set the default value for the integration interval of precipitation amount
     IF ( precipitation )  THEN
-       IF ( precipitation_amount_interval == 9999999.9 )  THEN
+       IF ( precipitation_amount_interval == 9999999.9_wp )  THEN
           precipitation_amount_interval = dt_do2d_xy
        ELSE
@@ -3277,5 +3277,5 @@
 !
 !-- Determine upper and lower hight level indices for random perturbations
-    IF ( disturbance_level_b == -9999999.9 )  THEN
+    IF ( disturbance_level_b == -9999999.9_wp )  THEN
        IF ( ocean ) THEN
           disturbance_level_b     = zu((nzt*2)/3)
@@ -3302,5 +3302,5 @@
     ENDIF
 
-    IF ( disturbance_level_t == -9999999.9 )  THEN
+    IF ( disturbance_level_t == -9999999.9_wp )  THEN
        IF ( ocean )  THEN
           disturbance_level_t     = zu(nzt-3)
@@ -3419,5 +3419,5 @@
 !-- In case of turbulent inflow calculate the index of the recycling plane
     IF ( turbulent_inflow )  THEN
-       IF ( recycling_width == 9999999.9 )  THEN
+       IF ( recycling_width == 9999999.9_wp )  THEN
 !
 !--       Set the default value for the width of the recycling domain
@@ -3482,5 +3482,5 @@
 !-- Set time for the next user defined restart (time_restart is the
 !-- internal parameter for steering restart events)
-    IF ( restart_time /= 9999999.9 )  THEN
+    IF ( restart_time /= 9999999.9_wp )  THEN
        IF ( restart_time > time_since_reference_point )  THEN
           time_restart = restart_time
@@ -3490,5 +3490,5 @@
 !--    In case of a restart run, set internal parameter to default (no restart)
 !--    if the NAMELIST-parameter restart_time is omitted
-       time_restart = 9999999.9
+       time_restart = 9999999.9_wp
     ENDIF
 
Index: /palm/trunk/SOURCE/cpulog.f90
===================================================================
--- /palm/trunk/SOURCE/cpulog.f90	(revision 1321)
+++ /palm/trunk/SOURCE/cpulog.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL functions provided with KIND-attribute
 !
 ! Former revisions:
@@ -165,7 +165,7 @@
 #if defined( __lc ) || defined( __decalpha ) || defined( __nec )
        CALL SYSTEM_CLOCK( count, count_rate )
-       mtime = REAL( count ) / REAL( count_rate )
+       mtime = REAL( count, KIND=wp ) / REAL( count_rate, KIND=wp )
 #elif defined( __ibm )
-       mtime = IRTC( ) * 1E-9
+       mtime = IRTC( ) * 1E-9_wp
 #else
        message_string = 'no time measurement defined on this host'
@@ -322,6 +322,6 @@
 !--       Get total time in order to calculate CPU-time per gridpoint and timestep
           IF ( nr_timesteps_this_run /= 0 )  THEN
-             average_cputime = log_point(1)%sum / REAL( (nx+1) * (ny+1) * nz ) / &
-                               REAL( nr_timesteps_this_run ) * 1E6  ! in micro-sec
+             average_cputime = log_point(1)%sum / REAL( (nx+1) * (ny+1) * nz, KIND=wp ) / &
+                               REAL( nr_timesteps_this_run, KIND=wp ) * 1E6_wp  ! in micro-sec
           ELSE
              average_cputime = -1.0
Index: /palm/trunk/SOURCE/data_output_profiles.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_profiles.f90	(revision 1321)
+++ /palm/trunk/SOURCE/data_output_profiles.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL functions provided with KIND-attribute
 !
 ! Former revisions:
@@ -102,5 +102,5 @@
     ELSE
        IF ( average_count_pr > 0 )  THEN
-          hom_sum = hom_sum / REAL( average_count_pr )
+          hom_sum = hom_sum / REAL( average_count_pr, KIND=wp )
        ELSE
 !
Index: /palm/trunk/SOURCE/data_output_spectra.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_spectra.f90	(revision 1321)
+++ /palm/trunk/SOURCE/data_output_spectra.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL functions provided with KIND-attribute
 !
 ! Former revisions:
@@ -113,6 +113,6 @@
        ENDIF
        IF ( average_count_sp /= 1 )  THEN
-          spectrum_x = spectrum_x / REAL( average_count_sp )
-          spectrum_y = spectrum_y / REAL( average_count_sp )
+          spectrum_x = spectrum_x / REAL( average_count_sp, KIND=wp )
+          spectrum_y = spectrum_y / REAL( average_count_sp, KIND=wp )
           average_count_sp = 0
        ENDIF
Index: /palm/trunk/SOURCE/diffusion_w.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_w.f90	(revision 1321)
+++ /palm/trunk/SOURCE/diffusion_w.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -132,11 +132,11 @@
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmxp = 0.25 * &
+                kmxp = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
-                kmxm = 0.25 * &
+                kmxm = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
-                kmyp = 0.25 * &
+                kmyp = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
-                kmym = 0.25 * &
+                kmym = 0.25_wp *                                               &
                        ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
 
@@ -152,5 +152,5 @@
                       &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
                       &   ) * 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)  &
@@ -165,11 +165,11 @@
 !
 !--                Interpolate eddy diffusivities on staggered gridpoints
-                   kmxp = 0.25 * &
+                   kmxp = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
-                   kmxm = 0.25 * &
+                   kmxm = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
-                   kmyp = 0.25 * &
+                   kmyp = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
-                   kmym = 0.25 * &
+                   kmym = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
 
@@ -262,11 +262,11 @@
 !
 !--                Interpolate eddy diffusivities on staggered gridpoints
-                   kmxp = 0.25 * &
+                   kmxp = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
-                   kmxm = 0.25 * &
+                   kmxm = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
-                   kmyp = 0.25 * &
+                   kmyp = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
-                   kmym = 0.25 * &
+                   kmym = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
 
@@ -282,5 +282,5 @@
                          &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)    &
                          &   ) * 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) &
@@ -297,11 +297,11 @@
 !
 !--                Interpolate eddy diffusivities on staggered gridpoints
-                   kmxp = 0.25 * &
+                   kmxp = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
-                   kmxm = 0.25 * &
+                   kmxm = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
-                   kmyp = 0.25 * &
+                   kmyp = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
-                   kmym = 0.25 * &
+                   kmym = 0.25_wp *                                            &
                           ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
 
@@ -378,8 +378,8 @@
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
-          kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
-          kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
-          kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
-          kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+          kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
+          kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
+          kmyp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
+          kmym = 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)                                            &
@@ -394,5 +394,5 @@
                       &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
                       &   ) * 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)  &
@@ -423,8 +423,8 @@
 !
 !--          Interpolate eddy diffusivities on staggered gridpoints
-             kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
-             kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
-             kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
-             kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+             kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
+             kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
+             kmyp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
+             kmym = 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)                                         &
Index: /palm/trunk/SOURCE/diffusivities.f90
===================================================================
--- /palm/trunk/SOURCE/diffusivities.f90	(revision 1321)
+++ /palm/trunk/SOURCE/diffusivities.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL constants defined as wp-kind 
 !
 ! Former revisions:
@@ -141,9 +141,9 @@
                 IF ( dvar_dz > 0.0 ) THEN
                    IF ( use_single_reference_value )  THEN
-                      l_stable = 0.76 * sqrt_e /                               &
-                                 SQRT( g / var_reference * dvar_dz ) + 1E-5
+                      l_stable = 0.76_wp * sqrt_e /                            &
+                                 SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
                    ELSE
-                      l_stable = 0.76 * sqrt_e /                               &
-                                 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
+                      l_stable = 0.76_wp * sqrt_e /                            &
+                                 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
                    ENDIF
                 ELSE
@@ -162,6 +162,6 @@
       !
       !--       Compute diffusion coefficients for momentum and heat
-                km(k,j,i) = 0.1 * l * sqrt_e
-                kh(k,j,i) = ( 1.0 + 2.0 * l / ll ) * km(k,j,i)
+                km(k,j,i) = 0.1_wp * l * sqrt_e
+                kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / ll ) * km(k,j,i)
 
              ENDIF
Index: /palm/trunk/SOURCE/disturb_field.f90
===================================================================
--- /palm/trunk/SOURCE/disturb_field.f90	(revision 1321)
+++ /palm/trunk/SOURCE/disturb_field.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -98,6 +98,6 @@
           DO  j = dist_nys(dist_range), dist_nyn(dist_range)
              DO  k = disturbance_level_ind_b, disturbance_level_ind_t
-                randomnumber = 3.0 * disturbance_amplitude *                   &
-                               ( random_function( iran ) - 0.5 )
+                randomnumber = 3.0_wp * disturbance_amplitude *                &
+                               ( random_function( iran ) - 0.5_wp )
                 IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.         &
                      nyn >= j )                                                &
@@ -113,10 +113,10 @@
              DO  k = disturbance_level_ind_b, disturbance_level_ind_t
 #if defined( __nec )
-                randomnumber = 3.0 * disturbance_amplitude *                   &
-                               ( RANDOM( 0 ) - 0.5 )
+                randomnumber = 3.0_wp * disturbance_amplitude *                &
+                               ( RANDOM( 0 ) - 0.5_wp )
 #else
                 CALL RANDOM_NUMBER( randomnumber )
-                randomnumber = 3.0 * disturbance_amplitude *                   &
-                                ( randomnumber - 0.5 )
+                randomnumber = 3.0_wp * disturbance_amplitude *                &
+                                ( randomnumber - 0.5_wp )
 #endif
                 IF ( nxl <= i .AND. nxr >= i .AND. nys <= j .AND. nyn >= j )   &
@@ -145,10 +145,10 @@
              dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1)                  &
                             + dist1(k,j+1,i) + dist1(k+1,j,i)                  &
-                            ) / 12.0
+                            ) / 12.0_wp
           ENDDO
           DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
               dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i)  &
-                            + 6.0 * dist1(k,j,i)                               &
-                            ) / 12.0
+                            + 6.0_wp * dist1(k,j,i)                            &
+                            ) / 12.0_wp
           ENDDO
         ENDDO
@@ -165,6 +165,6 @@
              dist1(k,j,i) = ( dist2(k,j,i-1) + dist2(k,j,i+1) + dist2(k,j-1,i) &
                             + dist2(k,j+1,i) + dist2(k+1,j,i) + dist2(k-1,j,i) &
-                            + 6.0 * dist2(k,j,i)                               &
-                            ) / 12.0
+                            + 6.0_wp * dist2(k,j,i)                            &
+                            ) / 12.0_wp
           ENDDO
        ENDDO
Index: /palm/trunk/SOURCE/disturb_heatflux.f90
===================================================================
--- /palm/trunk/SOURCE/disturb_heatflux.f90	(revision 1321)
+++ /palm/trunk/SOURCE/disturb_heatflux.f90	(revision 1322)
@@ -19,6 +19,6 @@
 !
 ! Current revisions:
-! -----------------
-! 
+! ------------------
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -79,5 +79,5 @@
     DO  i = 0, nx
        DO  j = 0, ny
-          randomnumber = random_gauss( iran, 5.0 )
+          randomnumber = random_gauss( iran, 5.0_wp )
           IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.  nyn >= j )   &
           THEN
Index: /palm/trunk/SOURCE/fft_xy.f90
===================================================================
--- /palm/trunk/SOURCE/fft_xy.f90	(revision 1321)
+++ /palm/trunk/SOURCE/fft_xy.f90	(revision 1322)
@@ -16,9 +16,9 @@
 !
 ! Copyright 1997-2014 Leibniz Universitaet Hannover
-!--------------------------------------------------------------------------------!
+!------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! REAL functions provided with KIND-attribute
 !
 ! Former revisions:
@@ -413,5 +413,5 @@
 
                    DO  i = 0, (nx+1)/2
-                      ar(i,j,k) = REAL( cwork(i) )
+                      ar(i,j,k) = REAL( cwork(i), KIND=wp )
                    ENDDO
                    DO  i = 1, (nx+1)/2 - 1
@@ -441,5 +441,5 @@
 
                    DO  i = 0, nx
-                      ar(i,j,k) = REAL( cwork(i) )
+                      ar(i,j,k) = REAL( cwork(i), KIND=wp )
                    ENDDO
 
@@ -519,5 +519,5 @@
 
                       DO  i = 0, (nx+1)/2
-                         ar_2d(i,j) = REAL( x_out(i) ) / ( nx+1 )
+                         ar_2d(i,j) = REAL( x_out(i), KIND=wp ) / ( nx+1 )
                       ENDDO
                       DO  i = 1, (nx+1)/2 - 1
@@ -528,5 +528,5 @@
 
                       DO  i = 0, (nx+1)/2
-                         ar(i,j,k) = REAL( x_out(i) ) / ( nx+1 )
+                         ar(i,j,k) = REAL( x_out(i), KIND=wp ) / ( nx+1 )
                       ENDDO
                       DO  i = 1, (nx+1)/2 - 1
@@ -689,5 +689,5 @@
 
                    DO  i = 0, (nx+1)/2
-                      ar(i,j,k)      = REAL( ar_tmp(i,j,k) )  * dnx
+                      ar(i,j,k)      = REAL( ar_tmp(i,j,k), KIND=wp )  * dnx
                    ENDDO
 
@@ -794,5 +794,5 @@
              CALL FFTN( cwork, ishape )
              DO  i = 0, (nx+1)/2
-                ar(i) = REAL( cwork(i) )
+                ar(i) = REAL( cwork(i), KIND=wp )
              ENDDO
              DO  i = 1, (nx+1)/2 - 1
@@ -813,5 +813,5 @@
 
              DO  i = 0, nx
-                ar(i) = REAL( cwork(i) )
+                ar(i) = REAL( cwork(i), KIND=wp )
              ENDDO
 
@@ -862,5 +862,5 @@
 
              DO  i = 0, (nx+1)/2
-                ar(i) = REAL( x_out(i) ) / ( nx+1 )
+                ar(i) = REAL( x_out(i), KIND=wp ) / ( nx+1 )
              ENDDO
              DO  i = 1, (nx+1)/2 - 1
@@ -1056,5 +1056,5 @@
 
                    DO  j = 0, (ny+1)/2
-                      ar_tr(j,i,k) = REAL( cwork(j) )
+                      ar_tr(j,i,k) = REAL( cwork(j), KIND=wp )
                    ENDDO
                    DO  j = 1, (ny+1)/2 - 1
@@ -1084,5 +1084,5 @@
 
                    DO  j = 0, ny
-                      ar(j,i,k) = REAL( cwork(j) )
+                      ar(j,i,k) = REAL( cwork(j), KIND=wp )
                    ENDDO
 
@@ -1160,5 +1160,5 @@
 
                    DO  j = 0, (ny+1)/2
-                      ar_tr(j,i,k) = REAL( y_out(j) ) / (ny+1)
+                      ar_tr(j,i,k) = REAL( y_out(j), KIND=wp ) / (ny+1)
                    ENDDO
                    DO  j = 1, (ny+1)/2 - 1
@@ -1305,5 +1305,5 @@
 
                    DO  j = 0, (ny+1)/2
-                      ar(j,i,k)      = REAL( ar_tmp(j,i,k) )  * dny
+                      ar(j,i,k)      = REAL( ar_tmp(j,i,k), KIND=wp )  * dny
                    ENDDO
 
@@ -1412,5 +1412,5 @@
 
              DO  j = 0, (ny+1)/2
-                ar(j) = REAL( cwork(j) )
+                ar(j) = REAL( cwork(j), KIND=wp )
              ENDDO
              DO  j = 1, (ny+1)/2 - 1
@@ -1431,5 +1431,5 @@
 
              DO  j = 0, ny
-                ar(j) = REAL( cwork(j) )
+                ar(j) = REAL( cwork(j), KIND=wp )
              ENDDO
 
@@ -1480,5 +1480,5 @@
 
              DO  j = 0, (ny+1)/2
-                ar(j) = REAL( y_out(j) ) / (ny+1)
+                ar(j) = REAL( y_out(j), KIND=wp ) / (ny+1)
              ENDDO
              DO  j = 1, (ny+1)/2 - 1
@@ -1674,5 +1674,5 @@
              DO  k = 1, nz
                 DO  i = 0, (nx+1)/2
-                   ar(i,k) = REAL( work(i+1,k) )
+                   ar(i,k) = REAL( work(i+1,k), KIND=wp )
                 ENDDO
                 DO  i = 1, (nx+1)/2 - 1
@@ -1817,5 +1817,5 @@
              DO  k = 1, nz
                 DO  j = 0, (ny+1)/2
-                   ar(j,k) = REAL( work(j+1,k) )
+                   ar(j,k) = REAL( work(j+1,k), KIND=wp )
                 ENDDO
                 DO  j = 1, (ny+1)/2 - 1
Index: /palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- /palm/trunk/SOURCE/flow_statistics.f90	(revision 1321)
+++ /palm/trunk/SOURCE/flow_statistics.f90	(revision 1322)
@@ -21,5 +21,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -1193,5 +1193,5 @@
           dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
        ENDDO
-       dptdz_threshold = 0.2 / 100.0
+       dptdz_threshold = 0.2_wp / 100.0_wp
 
        IF ( ocean )  THEN
@@ -2859,5 +2859,5 @@
           dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
        ENDDO
-       dptdz_threshold = 0.2 / 100.0
+       dptdz_threshold = 0.2_wp / 100.0_wp
 
        IF ( ocean )  THEN
Index: /palm/trunk/SOURCE/header.f90
===================================================================
--- /palm/trunk/SOURCE/header.f90	(revision 1321)
+++ /palm/trunk/SOURCE/header.f90	(revision 1322)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL functions provided with KIND-attribute,
+! some REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -454,7 +455,7 @@
     WRITE ( io, 203 )  simulated_time_at_begin, end_time
 
-    IF ( time_restart /= 9999999.9  .AND. &
+    IF ( time_restart /= 9999999.9_wp  .AND. &
          simulated_time_at_begin == simulated_time )  THEN
-       IF ( dt_restart == 9999999.9 )  THEN
+       IF ( dt_restart == 9999999.9_wp )  THEN
           WRITE ( io, 204 )  ' Restart at:       ',time_restart
        ELSE
@@ -472,9 +473,9 @@
                                               simulated_time_at_begin )
        ENDIF
-       WRITE ( io, 206 )  simulated_time, log_point_s(10)%sum, &
-                          log_point_s(10)%sum / REAL( i ),     &
+       WRITE ( io, 206 )  simulated_time, log_point_s(10)%sum,      &
+                          log_point_s(10)%sum / REAL( i, KIND=wp ), &
                           cpuseconds_per_simulated_second
-       IF ( time_restart /= 9999999.9  .AND.  time_restart < end_time )  THEN
-          IF ( dt_restart == 9999999.9 )  THEN
+       IF ( time_restart /= 9999999.9_wp  .AND.  time_restart < end_time )  THEN
+          IF ( dt_restart == 9999999.9_wp )  THEN
              WRITE ( io, 204 )  ' Next restart at:     ',time_restart
           ELSE
@@ -524,5 +525,5 @@
           bh  = INT( building_height / dz )
 
-          IF ( building_wall_left == 9999999.9 )  THEN
+          IF ( building_wall_left == 9999999.9_wp )  THEN
              building_wall_left = ( nx + 1 - blx ) / 2 * dx
           ENDIF
@@ -530,5 +531,5 @@
           bxr = bxl + blx
 
-          IF ( building_wall_south == 9999999.9 )  THEN
+          IF ( building_wall_south == 9999999.9_wp )  THEN
              building_wall_south = ( ny + 1 - bly ) / 2 * dy
           ENDIF
@@ -541,9 +542,9 @@
        CASE ( 'single_street_canyon' )
           ch  = NINT( canyon_height / dz )
-          IF ( canyon_width_x /= 9999999.9 )  THEN
+          IF ( canyon_width_x /= 9999999.9_wp )  THEN
 !
 !--          Street canyon in y direction
              cwx = NINT( canyon_width_x / dx )
-             IF ( canyon_wall_left == 9999999.9 )  THEN
+             IF ( canyon_wall_left == 9999999.9_wp )  THEN
                 canyon_wall_left = ( nx + 1 - cwx ) / 2 * dx
              ENDIF
@@ -552,9 +553,9 @@
              WRITE ( io, 272 )  'y', canyon_height, ch, 'u', cxl, cxr
 
-          ELSEIF ( canyon_width_y /= 9999999.9 )  THEN
+          ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
 !
 !--          Street canyon in x direction
              cwy = NINT( canyon_width_y / dy )
-             IF ( canyon_wall_south == 9999999.9 )  THEN
+             IF ( canyon_wall_south == 9999999.9_wp )  THEN
                 canyon_wall_south = ( ny + 1 - cwy ) / 2 * dy
              ENDIF
@@ -1188,5 +1189,5 @@
 !
 !-- Timeseries
-    IF ( dt_dots /= 9999999.9 )  THEN
+    IF ( dt_dots /= 9999999.9_wp )  THEN
        WRITE ( io, 340 )
 
@@ -1202,5 +1203,5 @@
 !
 !-- Dvrp-output
-    IF ( dt_dvrp /= 9999999.9 )  THEN
+    IF ( dt_dvrp /= 9999999.9_wp )  THEN
        WRITE ( io, 360 )  dt_dvrp, TRIM( dvrp_output ), TRIM( dvrp_host ), &
                           TRIM( dvrp_username ), TRIM( dvrp_directory )
@@ -1254,5 +1255,5 @@
 !
 !-- Spectra output
-    IF ( dt_dosp /= 9999999.9 )  THEN
+    IF ( dt_dosp /= 9999999.9_wp )  THEN
        WRITE ( io, 370 )
 
@@ -1377,5 +1378,5 @@
 !
 !-- Initial wind profiles
-    IF ( u_profile(1) /= 9999999.9 )  WRITE ( io, 427 )
+    IF ( u_profile(1) /= 9999999.9_wp )  WRITE ( io, 427 )
 
 !
@@ -1627,5 +1628,5 @@
           ENDIF
        ENDIF
-       IF ( dt_write_particle_data /= 9999999.9 )  THEN
+       IF ( dt_write_particle_data /= 9999999.9_wp )  THEN
           WRITE ( io, 485 )  dt_write_particle_data
           output_format = ' '
@@ -1641,5 +1642,5 @@
           WRITE ( io, 344 )  output_format
        ENDIF
-       IF ( dt_dopts /= 9999999.9 )  WRITE ( io, 494 )  dt_dopts
+       IF ( dt_dopts /= 9999999.9_wp )  WRITE ( io, 494 )  dt_dopts
        IF ( write_particle_statistics )  WRITE ( io, 486 )
 
@@ -1647,5 +1648,5 @@
 
        DO  i = 1, number_of_particle_groups
-          IF ( i == 1  .AND.  density_ratio(i) == 9999999.9 )  THEN
+          IF ( i == 1  .AND.  density_ratio(i) == 9999999.9_wp )  THEN
              WRITE ( io, 490 )  i, 0.0
              WRITE ( io, 492 )
Index: /palm/trunk/SOURCE/init_1d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_1d_model.f90	(revision 1321)
+++ /palm/trunk/SOURCE/init_1d_model.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL functions provided with KIND-attribute
 !
 ! Former revisions:
@@ -777,5 +777,5 @@
           energy = energy + 0.5 * ( u1d(k)**2 + v1d(k)**2 )
        ENDDO
-       energy = energy / REAL( nzt - nzb + 1 )
+       energy = energy / REAL( nzt - nzb + 1, KIND=wp )
 
        uv_total = SQRT( u1d(nzb+1)**2 + v1d(nzb+1)**2 )
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 1321)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 1322)
@@ -20,6 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
-!
+! REAL constants defined as wp-kind
 ! module interfaces removed
 !
@@ -1604,19 +1603,19 @@
     IF ( TRIM(timestep_scheme) == 'runge-kutta-3' )  THEN      ! for RK3-method
 
-       weight_substep(1) = 1./6.
-       weight_substep(2) = 3./10.
-       weight_substep(3) = 8./15.
-
-       weight_pres(1)    = 1./3.
-       weight_pres(2)    = 5./12.
-       weight_pres(3)    = 1./4.
+       weight_substep(1) = 1._wp/6._wp
+       weight_substep(2) = 3._wp/10._wp
+       weight_substep(3) = 8._wp/15._wp
+
+       weight_pres(1)    = 1._wp/3._wp
+       weight_pres(2)    = 5._wp/12._wp
+       weight_pres(3)    = 1._wp/4._wp
 
     ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' )  THEN  ! for RK2-method
 
-       weight_substep(1) = 1./2.
-       weight_substep(2) = 1./2.
+       weight_substep(1) = 1._wp/2._wp
+       weight_substep(2) = 1._wp/2._wp
           
-       weight_pres(1)    = 1./2.
-       weight_pres(2)    = 1./2.        
+       weight_pres(1)    = 1._wp/2._wp
+       weight_pres(2)    = 1._wp/2._wp        
 
     ELSE                                     ! for Euler-method
@@ -1671,6 +1670,6 @@
           DO  k = dp_level_ind_b+1, nzt
              dp_smooth_factor(k) = 0.5 * ( 1.0 + SIN( pi * &
-                  ( REAL( k - dp_level_ind_b ) /  &
-                    REAL( nzt - dp_level_ind_b ) - 0.5 ) ) )
+                  ( REAL( k - dp_level_ind_b, KIND=wp ) /  &
+                    REAL( nzt - dp_level_ind_b, KIND=wp ) - 0.5 ) ) )
           ENDDO
        ENDIF
@@ -1686,7 +1685,7 @@
        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 ) / (        &
-                         REAL( pt_damping_width )            ) ) )**2 
+             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 
           ENDIF
        ENDDO
@@ -1694,7 +1693,8 @@
        DO  i = nxl, nxr
           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 ) )**2
+             ptdf_x(i) = pt_damping_factor *                                   &
+                         SIN( pi * 0.5 *                                       &
+                              ( ( i - nx ) * dx + pt_damping_width ) /         &
+                              REAL( pt_damping_width, KIND=wp ) )**2
           ENDIF
        ENDDO 
@@ -1702,7 +1702,8 @@
        DO  j = nys, nyn
           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 ) )**2
+             ptdf_y(j) = pt_damping_factor *                                   &
+                         SIN( pi * 0.5 *                                       &
+                              ( ( j - ny ) * dy + pt_damping_width ) /         &
+                              REAL( pt_damping_width, KIND=wp ) )**2
           ENDIF
        ENDDO 
@@ -1710,7 +1711,8 @@
        DO  j = nys, nyn
           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 ) )**2
+             ptdf_y(j) = pt_damping_factor *                                   &
+                         SIN( pi * 0.5 *                                       &
+                             ( pt_damping_width - j * dy ) /                   &
+                             REAL( pt_damping_width, KIND=wp ) )**2
           ENDIF
        ENDDO
Index: /palm/trunk/SOURCE/init_advec.f90
===================================================================
--- /palm/trunk/SOURCE/init_advec.f90	(revision 1321)
+++ /palm/trunk/SOURCE/init_advec.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -85,5 +85,5 @@
        ALLOCATE( aex(intervals), bex(intervals), dex(intervals), eex(intervals) )
 
-       delt  = 1.0 / REAL( intervals )
+       delt  = 1.0_wp / REAL( intervals, KIND=wp )
        sterm = delt * 0.5
 
Index: /palm/trunk/SOURCE/init_cloud_physics.f90
===================================================================
--- /palm/trunk/SOURCE/init_cloud_physics.f90	(revision 1321)
+++ /palm/trunk/SOURCE/init_cloud_physics.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -96,5 +96,5 @@
     l_d_rv = l_v / r_v
 
-    schmidt_p_1d3 = schmidt**( 1.0 / 3.0 )
+    schmidt_p_1d3 = schmidt**( 1.0_wp / 3.0_wp )
 
     pirho_l  = pi * rho_l / 6.0
@@ -125,5 +125,5 @@
        ENDIF
        hyp(k)    = surface_pressure * 100.0 * &
-                   ( (t_surface - g/cp * zu(k)) / t_surface )**(1.0/0.286)
+                   ( (t_surface - g/cp * zu(k)) / t_surface )**(1.0_wp/0.286_wp)
        pt_d_t(k) = ( 100000.0 / hyp(k) )**0.286
        t_d_pt(k) = 1.0 / pt_d_t(k)
Index: /palm/trunk/SOURCE/init_dvrp.f90
===================================================================
--- /palm/trunk/SOURCE/init_dvrp.f90	(revision 1321)
+++ /palm/trunk/SOURCE/init_dvrp.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -710,13 +710,13 @@
 !--       Define a default colourtable for particles
           DO  i = 1, 11
-             interval_values_dvrp_prt(1,i) = i - 1.0
-             interval_values_dvrp_prt(2,i) = REAL( i )
-             interval_h_dvrp_prt(:,i) = 270.0 - ( i - 1.0 ) * 9.0
+             interval_values_dvrp_prt(1,i) = i - 1.0_wp
+             interval_values_dvrp_prt(2,i) = REAL( i, KIND=wp )
+             interval_h_dvrp_prt(:,i) = 270.0_wp - ( i - 1.0_wp ) * 9.0_wp
           ENDDO
 
           DO  i = 12, 22
-             interval_values_dvrp_prt(1,i) = i - 1.0
-             interval_values_dvrp_prt(2,i) = REAL( i )
-             interval_h_dvrp_prt(:,i) = 70.0 - ( i - 12.0 ) * 9.5
+             interval_values_dvrp_prt(1,i) = i - 1.0_wp
+             interval_values_dvrp_prt(2,i) = REAL( i, KIND=wp )
+             interval_h_dvrp_prt(:,i) = 70.0_wp - ( i - 12.0_wp ) * 9.5_wp
           ENDDO
 
Index: /palm/trunk/SOURCE/init_grid.f90
===================================================================
--- /palm/trunk/SOURCE/init_grid.f90	(revision 1321)
+++ /palm/trunk/SOURCE/init_grid.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -377,5 +377,5 @@
 !-- Compute the grid-dependent mixing length.
     DO  k = 1, nzt
-       l_grid(k)  = ( dx * dy * dzw(k) )**0.33333333333333
+       l_grid(k)  = ( dx * dy * dzw(k) )**0.33333333333333_wp
     ENDDO
 
@@ -516,5 +516,5 @@
           bh  = NINT( building_height / dz )
 
-          IF ( building_wall_left == 9999999.9 )  THEN
+          IF ( building_wall_left == 9999999.9_wp )  THEN
              building_wall_left = ( nx + 1 - blx ) / 2 * dx
           ENDIF
@@ -522,5 +522,5 @@
           bxr = bxl + blx
 
-          IF ( building_wall_south == 9999999.9 )  THEN
+          IF ( building_wall_south == 9999999.9_wp )  THEN
              building_wall_south = ( ny + 1 - bly ) / 2 * dy
           ENDIF
@@ -547,9 +547,9 @@
 !--       Single quasi-2D street canyon of infinite length in x or y direction.
 !--       The canyon is centered in the other direction by default.
-          IF ( canyon_width_x /= 9999999.9 )  THEN
+          IF ( canyon_width_x /= 9999999.9_wp )  THEN
 !
 !--          Street canyon in y direction
              cwx = NINT( canyon_width_x / dx )
-             IF ( canyon_wall_left == 9999999.9 )  THEN
+             IF ( canyon_wall_left == 9999999.9_wp )  THEN
                 canyon_wall_left = ( nx + 1 - cwx ) / 2 * dx
              ENDIF
@@ -557,9 +557,9 @@
              cxr = cxl + cwx
 
-          ELSEIF ( canyon_width_y /= 9999999.9 )  THEN
+          ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
 !
 !--          Street canyon in x direction
              cwy = NINT( canyon_width_y / dy )
-             IF ( canyon_wall_south == 9999999.9 )  THEN
+             IF ( canyon_wall_south == 9999999.9_wp )  THEN
                 canyon_wall_south = ( ny + 1 - cwy ) / 2 * dy
              ENDIF
@@ -578,5 +578,5 @@
 !
 !--       Street canyon size has to meet some requirements
-          IF ( canyon_width_x /= 9999999.9 )  THEN
+          IF ( canyon_width_x /= 9999999.9_wp )  THEN
              IF ( ( cxl < 1 ) .OR. ( cxr > nx-1 ) .OR. ( cwx < 3 ) .OR.  &
                ( ch < 3 ) )  THEN
@@ -587,5 +587,5 @@
                 CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 ) 
              ENDIF
-          ELSEIF ( canyon_width_y /= 9999999.9 )  THEN
+          ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
              IF ( ( cys < 1 ) .OR. ( cyn > ny-1 ) .OR. ( cwy < 3 ) .OR.  &
                ( ch < 3 ) )  THEN
@@ -597,5 +597,5 @@
              ENDIF
           ENDIF
-          IF ( canyon_width_x /= 9999999.9 .AND. canyon_width_y /= 9999999.9 ) &
+          IF ( canyon_width_x /= 9999999.9_wp .AND. canyon_width_y /= 9999999.9_wp ) &
                THEN
              message_string = 'inconsistent canyon parameters:' //     &  
@@ -606,7 +606,7 @@
 
           nzb_local = ch
-          IF ( canyon_width_x /= 9999999.9 )  THEN
+          IF ( canyon_width_x /= 9999999.9_wp )  THEN
              nzb_local(:,cxl+1:cxr-1) = 0
-          ELSEIF ( canyon_width_y /= 9999999.9 )  THEN
+          ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
              nzb_local(cys+1:cyn-1,:) = 0
           ENDIF
Index: /palm/trunk/SOURCE/init_ocean.f90
===================================================================
--- /palm/trunk/SOURCE/init_ocean.f90	(revision 1321)
+++ /palm/trunk/SOURCE/init_ocean.f90	(revision 1322)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! -----------------
+! REAL constants defined as wp_kind
 !
 ! Former revisions:
@@ -87,5 +88,5 @@
 !-- and the reference density (used later in buoyancy term)
 !-- First step: Calculate pressure using reference density
-    hyp(nzt+1) = surface_pressure * 100.0
+    hyp(nzt+1) = surface_pressure * 100.0_wp
 
     hyp(nzt)      = hyp(nzt+1) + rho_surface * g * 0.5 * dzu(nzt+1)
Index: /palm/trunk/SOURCE/init_pegrid.f90
===================================================================
--- /palm/trunk/SOURCE/init_pegrid.f90	(revision 1321)
+++ /palm/trunk/SOURCE/init_pegrid.f90	(revision 1322)
@@ -19,5 +19,6 @@
 !
 ! Current revisions:
-! -----------------
+! ------------------
+! REAL functions provided with KIND-attribute
 !
 ! Former revisions:
@@ -174,5 +175,5 @@
 !
 !--    Automatic determination of the topology
-       numproc_sqr = SQRT( REAL( numprocs ) )
+       numproc_sqr = SQRT( REAL( numprocs, KIND=wp ) )
        pdims(1)    = MAX( numproc_sqr , 1 )
        DO  WHILE ( MOD( numprocs , pdims(1) ) /= 0 )
Index: /palm/trunk/SOURCE/init_pt_anomaly.f90
===================================================================
--- /palm/trunk/SOURCE/init_pt_anomaly.f90	(revision 1321)
+++ /palm/trunk/SOURCE/init_pt_anomaly.f90	(revision 1322)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! -----------------
+! REAL constants defined as wp_kind
 !
 ! Former revisions:
@@ -93,5 +94,5 @@
              radius = SQRT( x**2 + y**2 + z**2 )
              IF ( radius <= rc )  THEN
-                betrag = 5.0 * EXP( -( radius * 0.001 / 2.0 )**2 )
+                betrag = 5.0 * EXP( -( radius * 0.001_wp / 2.0_wp )**2 )
              ELSE
                 betrag = 0.0
Index: /palm/trunk/SOURCE/init_rankine.f90
===================================================================
--- /palm/trunk/SOURCE/init_rankine.f90	(revision 1321)
+++ /palm/trunk/SOURCE/init_rankine.f90	(revision 1322)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! -----------------
+! REAL constants defined as wp_kind
 !
 ! Former revisions:
@@ -111,7 +112,7 @@
           radius = SQRT( x**2 + y**2 )
           IF ( radius <= 2.0 * rc )  THEN
-             betrag = radius / ( 2.0 * rc ) * 0.08
+             betrag = radius / ( 2.0 * rc ) * 0.08_wp
           ELSEIF ( radius > 2.0 * rc  .AND.  radius < 8.0 * rc )  THEN
-             betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 )
+             betrag = 0.08 * EXP( -( radius - 2.0_wp * rc ) / 2.0_wp )
           ELSE
              betrag = 0.0
@@ -120,7 +121,7 @@
           IF ( x == 0.0 )  THEN
              IF ( y > 0.0 )  THEN
-                alpha = pi / 2.0
+                alpha = pi / 2.0_wp
              ELSEIF ( y < 0.0 )  THEN
-                alpha = 3.0 * pi / 2.0
+                alpha = 3.0 * pi / 2.0_wp
              ENDIF
           ELSE
@@ -129,5 +130,5 @@
              ELSE
                 IF ( y < 0.0 )  THEN
-                   alpha = ATAN( y / x ) + 2.0 * pi
+                   alpha = ATAN( y / x ) + 2.0_wp * pi
                 ELSE
                    alpha = ATAN( y / x )
@@ -152,7 +153,7 @@
           radius = SQRT( x**2 + y**2 )
           IF ( radius <= 2.0 * rc )  THEN
-             betrag = radius / ( 2.0 * rc ) * 0.08
+             betrag = radius / ( 2.0_wp * rc ) * 0.08_wp
           ELSEIF ( radius > 2.0 * rc  .AND.  radius < 8.0 * rc )  THEN
-             betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 )
+             betrag = 0.08 * EXP( -( radius - 2.0_wp * rc ) / 2.0_wp )
           ELSE
              betrag = 0.0
@@ -161,7 +162,7 @@
           IF ( x == 0.0 )  THEN
              IF ( y > 0.0 )  THEN
-                alpha = pi / 2.0
+                alpha = pi / 2.0_wp
              ELSEIF ( y < 0.0 )  THEN
-                alpha = 3.0 * pi / 2.0
+                alpha = 3.0 * pi / 2.0_wp
              ENDIF
           ELSE
@@ -170,5 +171,5 @@
              ELSE
                 IF ( y < 0.0 )  THEN
-                   alpha = ATAN( y / x ) + 2.0 * pi
+                   alpha = ATAN( y / x ) + 2.0_wp * pi
                 ELSE
                    alpha = ATAN( y / x )
Index: /palm/trunk/SOURCE/init_slope.f90
===================================================================
--- /palm/trunk/SOURCE/init_slope.f90	(revision 1321)
+++ /palm/trunk/SOURCE/init_slope.f90	(revision 1322)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! -----------------
+! REAL constants defined as wp_kind
 !
 ! Former revisions:
@@ -109,7 +110,7 @@
 !
 !--       Compute temperatures in the rotated coordinate system
-          alpha    = alpha + alpha_surface / 180.0 * pi
+          alpha    = alpha + alpha_surface / 180.0_wp * pi
           pt_value = pt_surface + radius * SIN( alpha ) * &
-                                  pt_vertical_gradient(1) / 100.0
+                                  pt_vertical_gradient(1) / 100.0_wp
           pt_slope_ref(k,i) = pt_value
        ENDDO                
@@ -120,5 +121,5 @@
 !-- used for the cyclic boundary in x-direction
     pt_slope_offset = (nx+1) * dx * sin_alpha_surface * &
-                      pt_vertical_gradient(1) / 100.0
+                      pt_vertical_gradient(1) / 100.0_wp
 
 
Index: /palm/trunk/SOURCE/local_tremain.f90
===================================================================
--- /palm/trunk/SOURCE/local_tremain.f90	(revision 1321)
+++ /palm/trunk/SOURCE/local_tremain.f90	(revision 1322)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! -----------------
+! REAL constants defined as wp_kind
 !
 ! Former revisions:
@@ -69,5 +70,5 @@
 #if defined( __ibm )
 
-    actual_wallclock_time = IRTC( ) * 1E-9
+    actual_wallclock_time = IRTC( ) * 1E-9_wp
     remaining_time = maximum_cpu_time_allowed - &
                      ( actual_wallclock_time - initial_wallclock_time )
@@ -76,5 +77,5 @@
 
     CALL SYSTEM_CLOCK( count, count_rate )
-    actual_wallclock_time = REAL( count ) / REAL( count_rate )
+    actual_wallclock_time = REAL( count, KIND=wp ) / REAL( count_rate, KIND=wp )
     remaining_time = maximum_cpu_time_allowed - &
                      ( actual_wallclock_time - initial_wallclock_time )
Index: /palm/trunk/SOURCE/local_tremain_ini.f90
===================================================================
--- /palm/trunk/SOURCE/local_tremain_ini.f90	(revision 1321)
+++ /palm/trunk/SOURCE/local_tremain_ini.f90	(revision 1322)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! -----------------
+! REAL constants defined as wp_kind
 !
 ! Former revisions:
@@ -67,10 +68,10 @@
 #if defined( __ibm )
 
-    initial_wallclock_time = IRTC( ) * 1E-9
+    initial_wallclock_time = IRTC( ) * 1E-9_wp
 
 #elif defined( __lc )
 
     CALL SYSTEM_CLOCK( count, count_rate )
-    initial_wallclock_time = REAL( count ) / REAL( count_rate )
+    initial_wallclock_time = REAL( count, KIND=wp ) / REAL( count_rate, KIND=wp )
 
 #else
Index: /palm/trunk/SOURCE/lpm_advec.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_advec.f90	(revision 1321)
+++ /palm/trunk/SOURCE/lpm_advec.f90	(revision 1322)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! ------------------
+! REAL constants defined as wp_kind
 !
 ! Former revisions:
@@ -186,6 +187,6 @@
 !
 !--          Determine the sublayer. Further used as index.
-             height_p     = ( particles(n)%z - z0_av_global )           &
-                                  * REAL( number_of_sublayers )         &
+             height_p     = ( particles(n)%z - z0_av_global )            &
+                                  * REAL( number_of_sublayers, KIND=wp ) &
                                   * d_z_p_z0 
 !
@@ -210,5 +211,5 @@
              us_int = 0.5 * ( us(j,i) + us(j,i-1) )  
 
-             u_int  = -usws(j,i) / ( us_int * kappa + 1E-10 )           & 
+             u_int  = -usws(j,i) / ( us_int * kappa + 1E-10_wp )           & 
                       * log_z_z0_int 
 
@@ -277,5 +278,5 @@
               us_int = 0.5 * ( us(j,i) + us(j-1,i) )    
 
-              v_int  = -vsws(j,i) / ( us_int * kappa + 1E-10 )             &
+              v_int  = -vsws(j,i) / ( us_int * kappa + 1E-10_wp )             &
                        * log_z_z0_int 
 
@@ -1093,8 +1094,8 @@
           ENDIF
 
-          vv_int = ( 1.0 / 3.0 ) * ( aa + bb + cc )
-
-          fs_int = ( 2.0 / 3.0 ) * e_mean_int / &
-                      ( vv_int + ( 2.0 / 3.0 ) * e_mean_int )
+          vv_int = ( 1.0_wp / 3.0_wp ) * ( aa + bb + cc )
+
+          fs_int = ( 2.0_wp / 3.0_wp ) * e_mean_int / &
+                      ( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int )
 
 !
@@ -1126,9 +1127,9 @@
 !--          from becoming unrealistically large.
              particles(n)%rvar1 = SQRT( 2.0 * sgs_wfu_part * e_int ) * &
-                                        ( random_gauss( iran_part, 5.0 ) - 1.0 )
+                                        ( random_gauss( iran_part, 5.0_wp ) - 1.0_wp )
              particles(n)%rvar2 = SQRT( 2.0 * sgs_wfv_part * e_int ) * &
-                                        ( random_gauss( iran_part, 5.0 ) - 1.0 )
+                                        ( random_gauss( iran_part, 5.0_wp ) - 1.0_wp )
              particles(n)%rvar3 = SQRT( 2.0 * sgs_wfw_part * e_int ) * &
-                                        ( random_gauss( iran_part, 5.0 ) - 1.0 )
+                                        ( random_gauss( iran_part, 5.0_wp ) - 1.0_wp )
 
           ELSE
@@ -1168,5 +1169,5 @@
                                   ) * dt_particle / 2.0          +             &
                                   SQRT( fs_int * c_0 * diss_int ) *            &
-                                  ( random_gauss( iran_part, 5.0 ) - 1.0 ) *   &
+                                  ( random_gauss( iran_part, 5.0_wp ) - 1.0_wp ) *   &
                                   SQRT( dt_particle )
 
@@ -1177,7 +1178,7 @@
                                     particles(n)%rvar2 /                       &
                                     ( 2.0 * sgs_wfv_part * e_int ) + de_dy_int &
-                                  ) * dt_particle / 2.0          +             &
+                                  ) * dt_particle / 2.0_wp          +             &
                                   SQRT( fs_int * c_0 * diss_int ) *            &
-                                  ( random_gauss( iran_part, 5.0 ) - 1.0 ) *   &
+                                  ( random_gauss( iran_part, 5.0_wp ) - 1.0_wp ) *   &
                                   SQRT( dt_particle )
 
@@ -1188,7 +1189,7 @@
                                     particles(n)%rvar3 /                       &
                                     ( 2.0 * sgs_wfw_part * e_int ) + de_dz_int &
-                                  ) * dt_particle / 2.0          +             &
+                                  ) * dt_particle / 2.0_wp         +             &
                                   SQRT( fs_int * c_0 * diss_int ) *            &
-                                  ( random_gauss( iran_part, 5.0 ) - 1.0 ) *   &
+                                  ( random_gauss( iran_part, 5.0_wp ) - 1.0_wp ) *   &
                                   SQRT( dt_particle )
 
@@ -1252,5 +1253,5 @@
                                 ( v_int - particles(n)%speed_y )**2 +   &
                                 ( w_int - particles(n)%speed_z )**2 ) / &
-                                         molecular_viscosity )**0.687   &
+                                         molecular_viscosity )**0.687_wp   &
                         )
              exp_term = EXP( -exp_arg * dt_particle )
Index: /palm/trunk/SOURCE/lpm_collision_kernels.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_collision_kernels.f90	(revision 1321)
+++ /palm/trunk/SOURCE/lpm_collision_kernels.f90	(revision 1322)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! -----------------
+! REAL constants defined as wp_kind
 !
 ! Former revisions:
@@ -167,6 +168,6 @@
 !--       Calculate the radius class bounds with logarithmic distances
 !--       in the interval [1.0E-6, 2.0E-4] m
-          rclass_lbound = LOG( 1.0E-6 )
-          rclass_ubound = LOG( 2.0E-4 )
+          rclass_lbound = LOG( 1.0E-6_wp )
+          rclass_ubound = LOG( 2.0E-4_wp )
           radclass(1)   = 1.0E-6
           DO  i = 2, radius_classes
@@ -182,5 +183,5 @@
 !--       Set the class bounds for dissipation in interval [0.0, 0.1] m**2/s**3
           DO  i = 1, dissipation_classes
-             epsclass(i) = 0.1 * REAL( i ) / dissipation_classes
+             epsclass(i) = 0.1 * REAL( i, KIND=wp ) / dissipation_classes
 !             IF ( myid == 0 )  THEN
 !                PRINT*, 'i=', i, ' eps = ', epsclass(i)
@@ -197,5 +198,5 @@
 
              epsilon = epsclass(k)
-             urms    = 2.02 * ( epsilon / 0.04 )**( 1.0 / 3.0 )
+             urms    = 2.02 * ( epsilon / 0.04_wp )**( 1.0_wp / 3.0_wp )
 
              CALL turbsd
@@ -313,5 +314,5 @@
           epsilon = 0.0
        ENDIF
-       urms    = 2.02 * ( epsilon / 0.04 )**( 0.33333333333 )
+       urms    = 2.02 * ( epsilon / 0.04_wp )**( 0.33333333333_wp )
 
        IF ( wang_kernel  .AND.  epsilon > 1.0E-7 )  THEN
@@ -431,14 +432,14 @@
        ENDIF
 
-       lambda    = urms * SQRT( 15.0 * molecular_viscosity / epsilon )    ! in m
-       lambda_re = urms**2 * SQRT( 15.0 / epsilon / molecular_viscosity )
+       lambda    = urms * SQRT( 15.0_wp * molecular_viscosity / epsilon ) ! in m
+       lambda_re = urms**2 * SQRT( 15.0_wp / epsilon / molecular_viscosity )
        tl        = urms**2 / epsilon                       ! in s
        lf        = 0.5 * urms**3 / epsilon                 ! in m
        tauk      = SQRT( molecular_viscosity / epsilon )                  ! in s
-       eta       = ( molecular_viscosity**3 / epsilon )**0.25             ! in m
+       eta       = ( molecular_viscosity**3 / epsilon )**0.25_wp          ! in m
        vk        = eta / tauk
 
        ao = ( 11.0 + 7.0 * lambda_re ) / ( 205.0 + lambda_re )
-       tt = SQRT( 2.0 * lambda_re / ( SQRT( 15.0 ) * ao ) ) * tauk   ! in s
+       tt = SQRT( 2.0 * lambda_re / ( SQRT( 15.0_wp ) * ao ) ) * tauk   ! in s
 
        CALL fallg    ! gives winf in m/s
@@ -452,5 +453,5 @@
 !--    Calculate wr (from Aayala 2008b, page 38f)
        z   = tt / tl
-       be  = SQRT( 2.0 ) * lambda / lf
+       be  = SQRT( 2.0_wp ) * lambda / lf
        bbb = SQRT( 1.0 - 2.0 * be**2 )
        d1  = ( 1.0 + bbb ) / ( 2.0 * bbb )
@@ -504,6 +505,6 @@
              wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0 * v1v2xy   ! in m**2/s**2
              IF ( wrtur2xy < 0.0 )  wrtur2xy = 0.0
-             wrgrav2  = pi / 8.0 * ( winf(j) - winf(i) )**2
-             wrfin    = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) )   ! in m/s
+             wrgrav2  = pi / 8.0_wp * ( winf(j) - winf(i) )**2
+             wrfin    = SQRT( ( 2.0_wp / pi ) * ( wrtur2xy + wrgrav2) ) ! in m/s
 
 !
@@ -518,9 +519,9 @@
                    5.3406 * sst
              IF ( xx < 0.0 )  xx = 0.0
-             yy = 0.1886 * EXP( 20.306 / lambda_re )
+             yy = 0.1886 * EXP( 20.306_wp / lambda_re )
 
              c1_gr  =  xx / ( g / vk * tauk )**yy
 
-             ao_gr  = ao + ( pi / 8.0) * ( g / vk * tauk )**2
+             ao_gr  = ao + ( pi / 8.0_wp) * ( g / vk * tauk )**2
              fao_gr = 20.115 * SQRT( ao_gr / lambda_re )
              rc     = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta   ! in cm
@@ -652,5 +653,5 @@
           stb   = 32.0 * rho_a * ( rho_l - rho_a) * g / (3.0 * eta * eta)
           phy   = sigma**3 * rho_a**2 / ( eta**4 * g * ( rho_l - rho_a ) )
-          py    = phy**( 1.0 / 6.0 )
+          py    = phy**( 1.0_wp / 6.0_wp )
 
        ENDIF
@@ -685,5 +686,5 @@
              ENDIF
 
-             x = LOG( 16.0 * bond * py / 3.0 )
+             x = LOG( 16.0 * bond * py / 3.0_wp )
              y = 0.0
 
@@ -695,5 +696,5 @@
 
              IF ( radclass(j) > 0.0035 )  THEN
-                winf(j) = xrey * eta / ( 2.0 * rho_a * 0.0035 )
+                winf(j) = xrey * eta / ( 2.0 * rho_a * 0.0035_wp )
              ELSE
                 winf(j) = xrey * eta / ( 2.0 * rho_a * radclass(j) )
@@ -799,5 +800,5 @@
        ALLOCATE( ira(1:radius_classes) )
        DO  j = 1, radius_classes
-          particle_radius = radclass(j) * 1.0E6
+          particle_radius = radclass(j) * 1.0E6_wp
           DO  k = 1, 15
              IF ( particle_radius < r0(k) )  THEN
@@ -822,5 +823,5 @@
              IF ( ir < 16 )  THEN
                 IF ( ir >= 2 )  THEN
-                   pp = ( ( radclass(j) * 1.0E06 ) - r0(ir-1) ) / &
+                   pp = ( ( radclass(j) * 1.0E06_wp ) - r0(ir-1) ) / &
                         ( r0(ir) - r0(ir-1) )
                    qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
@@ -930,5 +931,5 @@
 
        DO  j = 1, radius_classes
-          particle_radius = radclass(j) * 1.0E6
+          particle_radius = radclass(j) * 1.0E6_wp
           DO  k = 1, 7
              IF ( particle_radius < r0(k) )  THEN
@@ -1039,6 +1040,6 @@
        REAL(wp), DIMENSION(1:9,1:19), SAVE ::  ef = 0.0          !:
 
-       mean_rm = mean_r * 1.0E06
-       rm      = r      * 1.0E06
+       mean_rm = mean_r * 1.0E06_wp
+       rm      = r      * 1.0E06_wp
 
        IF ( first )  THEN
Index: /palm/trunk/SOURCE/lpm_droplet_collision.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_droplet_collision.f90	(revision 1321)
+++ /palm/trunk/SOURCE/lpm_droplet_collision.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! REAL constants defined as wp_kind
 !
 ! Former revisions:
@@ -207,5 +207,5 @@
 !--                Determine dissipation class index of this gridbox
                    IF ( wang_kernel )  THEN
-                      eclass = INT( diss(k,j,i) * 1.0E4 / 1000.0 * &
+                      eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * &
                                     dissipation_classes ) + 1
                       epsilon = diss(k,j,i)
@@ -213,5 +213,5 @@
                       epsilon = 0.0
                    ENDIF
-                   IF ( hall_kernel .OR. epsilon * 1.0E4 < 0.001 )  THEN
+                   IF ( hall_kernel .OR. epsilon * 1.0E4_wp < 0.001 )  THEN
                       eclass = 0   ! Hall kernel is used
                    ELSE
@@ -270,5 +270,5 @@
 !--                   Change of the current droplet radius
                       rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/&
-                                       sum3 )**0.33333333333333
+                                       sum3 )**0.33333333333333_wp
 
                       IF ( weight(n-is+1) < 0.0 )  THEN
@@ -347,5 +347,5 @@
 !--                   Change of the current droplet radius
                       rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/&
-                                       sum3 )**0.33333333333333
+                                       sum3 )**0.33333333333333_wp
 
                       IF ( weight(n-is+1) < 0.0 )  THEN
@@ -534,5 +534,5 @@
 !
 !--                   Change in radius due to collision
-                      delta_r = effective_coll_efficiency / 3.0             &
+                      delta_r = effective_coll_efficiency / 3.0_wp          &
                                 * pi * sl_r3 * ddx * ddy / dz               &
                                 * SQRT( ( u_int - particles(n)%speed_x )**2 &
@@ -551,6 +551,6 @@
                       IF ( delta_v >= sl_r3  .AND.  sl_r3 > 0.0 )  THEN
 
-                         delta_r = ( ( sl_r3/particles(n)%weight_factor )   &
-                                     + particles(n)%radius**3 )**( 1./3. )  &
+                         delta_r = ( ( sl_r3/particles(n)%weight_factor )           &
+                                     + particles(n)%radius**3 )**( 1.0_wp/3.0_wp )  &
                                      - particles(n)%radius
 
@@ -612,5 +612,5 @@
 
                    IF ( wang_kernel )  THEN
-                      eclass = INT( diss(k,j,i) * 1.0E4 / 1000.0 * &
+                      eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * &
                                     dissipation_classes ) + 1
                       epsilon = diss(k,j,i)
@@ -618,5 +618,5 @@
                       epsilon = 0.0
                    ENDIF
-                   IF ( hall_kernel .OR. epsilon * 1.0E4 < 0.001 )  THEN
+                   IF ( hall_kernel .OR. epsilon * 1.0E4_wp < 0.001 )  THEN
                       eclass = 0   ! Hall kernel is used
                    ELSE
@@ -631,5 +631,5 @@
 
                    particles(psi)%radius = ( particles(psi)%radius**3 / &
-                                             sum3 )**0.33333333333333
+                                             sum3 )**0.33333333333333_wp
                    particles(psi)%weight_factor = particles(psi)%weight_factor &
                                                   * sum3
@@ -652,5 +652,5 @@
 
                    particles(psi)%radius = ( particles(psi)%radius**3 / &
-                                             sum3 )**0.33333333333333
+                                             sum3 )**0.33333333333333_wp
                    particles(psi)%weight_factor = particles(psi)%weight_factor &
                                                   * sum3
Index: /palm/trunk/SOURCE/lpm_droplet_condensation.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_droplet_condensation.f90	(revision 1321)
+++ /palm/trunk/SOURCE/lpm_droplet_condensation.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -153,27 +153,27 @@
  
 !-- Parameters for Rosenbrock method
-    REAL(wp), PARAMETER ::  a21 = 2.0          !:
-    REAL(wp), PARAMETER ::  a31 = 48.0/25.0    !:
-    REAL(wp), PARAMETER ::  a32 = 6.0/25.0     !:
-    REAL(wp), PARAMETER ::  b1 = 19.0/9.0      !:
-    REAL(wp), PARAMETER ::  b2 = 0.5           !:
-    REAL(wp), PARAMETER ::  b3 = 25.0/108.0    !:
-    REAL(wp), PARAMETER ::  b4 = 125.0/108.0   !:
-    REAL(wp), PARAMETER ::  c21 = -8.0         !:
-    REAL(wp), PARAMETER ::  c31 = 372.0/25.0   !:
-    REAL(wp), PARAMETER ::  c32 = 12.0/5.0     !:
-    REAL(wp), PARAMETER ::  c41 = -112.0/125.0 !:
-    REAL(wp), PARAMETER ::  c42 = -54.0/125.0  !:
-    REAL(wp), PARAMETER ::  c43 = -2.0/5.0     !:
-    REAL(wp), PARAMETER ::  errcon = 0.1296    !:
-    REAL(wp), PARAMETER ::  e1 = 17.0/54.0     !:
-    REAL(wp), PARAMETER ::  e2 = 7.0/36.0      !:
-    REAL(wp), PARAMETER ::  e3 = 0.0           !:
-    REAL(wp), PARAMETER ::  e4 = 125.0/108.0   !:
-    REAL(wp), PARAMETER ::  gam = 0.5          !:
-    REAL(wp), PARAMETER ::  grow = 1.5         !:
-    REAL(wp), PARAMETER ::  pgrow = -0.25      !:
-    REAL(wp), PARAMETER ::  pshrnk = -1.0/3.0  !:
-    REAL(wp), PARAMETER ::  shrnk = 0.5        !:
+    REAL(wp), PARAMETER ::  a21 = 2.0             !:
+    REAL(wp), PARAMETER ::  a31 = 48.0/25.0_wp    !:
+    REAL(wp), PARAMETER ::  a32 = 6.0/25.0_wp     !:
+    REAL(wp), PARAMETER ::  b1 = 19.0/9.0_wp      !:
+    REAL(wp), PARAMETER ::  b2 = 0.5              !:
+    REAL(wp), PARAMETER ::  b3 = 25.0/108.0_wp    !:
+    REAL(wp), PARAMETER ::  b4 = 125.0/108.0_wp   !:
+    REAL(wp), PARAMETER ::  c21 = -8.0            !:
+    REAL(wp), PARAMETER ::  c31 = 372.0/25.0_wp   !:
+    REAL(wp), PARAMETER ::  c32 = 12.0/5.0_wp     !:
+    REAL(wp), PARAMETER ::  c41 = -112.0/125.0_wp !:
+    REAL(wp), PARAMETER ::  c42 = -54.0/125.0_wp  !:
+    REAL(wp), PARAMETER ::  c43 = -2.0/5.0_wp     !:
+    REAL(wp), PARAMETER ::  errcon = 0.1296       !:
+    REAL(wp), PARAMETER ::  e1 = 17.0/54.0_wp     !:
+    REAL(wp), PARAMETER ::  e2 = 7.0/36.0_wp      !:
+    REAL(wp), PARAMETER ::  e3 = 0.0              !:
+    REAL(wp), PARAMETER ::  e4 = 125.0/108.0_wp   !:
+    REAL(wp), PARAMETER ::  gam = 0.5             !:
+    REAL(wp), PARAMETER ::  grow = 1.5            !:
+    REAL(wp), PARAMETER ::  pgrow = -0.25         !:
+    REAL(wp), PARAMETER ::  pshrnk = -1.0/3.0_wp  !:
+    REAL(wp), PARAMETER ::  shrnk = 0.5           !:
 
 
Index: /palm/trunk/SOURCE/lpm_init.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_init.f90	(revision 1321)
+++ /palm/trunk/SOURCE/lpm_init.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! REAL functions provided with KIND-attribute
 !
 ! Former revisions:
@@ -252,5 +252,5 @@
 !
 !--    Calculate vertical depth of the sublayers
-       height_int  = ( z_p - z0_av_global ) / REAL( number_of_sublayers ) 
+       height_int  = ( z_p - z0_av_global ) / REAL( number_of_sublayers, KIND=wp )
 !
 !--    Precalculate LOG(z/z0)
Index: /palm/trunk/SOURCE/microphysics.f90
===================================================================
--- /palm/trunk/SOURCE/microphysics.f90	(revision 1321)
+++ /palm/trunk/SOURCE/microphysics.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -522,5 +522,5 @@
 !
 !--             Weight averaged radius of cloud droplets:
-                rc = 0.5 * ( xc * dpirho_l )**( 1.0 / 3.0 )
+                rc = 0.5 * ( xc * dpirho_l )**( 1.0 / 3.0_wp )
 
                 alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0 + a_3 * nu_c )
@@ -529,5 +529,5 @@
 !
 !--             Mixing length (neglecting distance to ground and stratification)
-                l_mix = ( dx * dy * dzu(k) )**( 1.0 / 3.0 )
+                l_mix = ( dx * dy * dzu(k) )**( 1.0 / 3.0_wp )
 !
 !--             Limit dissipation rate according to Seifert, Nuijens and 
@@ -536,6 +536,6 @@
 !
 !--             Compute Taylor-microscale Reynolds number:
-                re_lambda = 6.0 / 11.0 * ( l_mix / c_const )**( 2.0 / 3.0 ) *  &
-                            SQRT( 15.0 / kin_vis_air ) * epsilon**( 1.0 / 6.0 )
+                re_lambda = 6.0 / 11.0 * ( l_mix / c_const )**( 2.0 / 3.0_wp ) *  &
+                            SQRT( 15.0 / kin_vis_air ) * epsilon**( 1.0 / 6.0_wp )
 !
 !--             The factor of 1.0E4 is needed to convert the dissipation rate
@@ -664,5 +664,5 @@
 !
 !--          Weight averaged diameter of rain drops:
-             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0 / 3.0 )
+             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0 / 3.0_wp )
 !
 !--          Collisional breakup rate (Seifert, 2008):
@@ -760,5 +760,5 @@
 !
 !--          Weight averaged diameter of rain drops:
-             dr = ( xr * dpirho_l )**( 1.0 / 3.0 )
+             dr = ( xr * dpirho_l )**( 1.0 / 3.0_wp )
 !
 !--          Compute ventilation factor and intercept parameter 
@@ -772,5 +772,5 @@
 !--             Slope parameter of gamma distribution (Seifert, 2008):
                 lambda_r = ( ( mu_r + 3.0 ) * ( mu_r + 2.0 ) *       &
-                             ( mu_r + 1.0 ) )**( 1.0 / 3.0 ) / dr
+                             ( mu_r + 1.0 ) )**( 1.0 / 3.0_wp ) / dr
 
                 mu_r_2   = mu_r + 2.0
@@ -849,5 +849,5 @@
 !
 !--    Sedimentation of cloud droplets (Heus et al., 2010):
-       sed_qc_const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))**( 2.0 / 3.0 ) *   &
+       sed_qc_const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))**( 2.0 / 3.0_wp ) *   &
                      EXP( 5.0 * LOG( sigma_gc )**2 )
 
@@ -856,6 +856,6 @@
        DO  k = nzt, nzb_s_inner(j,i)+1, -1
           IF ( qc_1d(k) > eps_sb )  THEN
-             sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0 / 3.0 ) * &
-                        ( qc_1d(k) * hyrho(k) )**( 5.0 / 3.0 )
+             sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0 / 3.0_wp ) * &
+                        ( qc_1d(k) * hyrho(k) )**( 5.0 / 3.0_wp )
           ELSE
              sed_qc(k) = 0.0
@@ -940,5 +940,5 @@
 !
 !--          Weight averaged diameter of rain drops:
-             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0 / 3.0 )
+             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0 / 3.0_wp )
 !
 !--          Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
@@ -948,5 +948,5 @@
 !--          Slope parameter of gamma distribution (Seifert, 2008):
              lambda_r = ( ( mu_r + 3.0 ) * ( mu_r + 2.0 ) *          &
-                        ( mu_r + 1.0 ) )**( 1.0 / 3.0 ) / dr
+                        ( mu_r + 1.0 ) )**( 1.0 / 3.0_wp ) / dr
 
              w_nr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 +          &
Index: /palm/trunk/SOURCE/netcdf.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf.f90	(revision 1321)
+++ /palm/trunk/SOURCE/netcdf.f90	(revision 1322)
@@ -23,5 +23,5 @@
 ! Current revisions:
 ! ------------------
-!
+! Forgotten ONLY-attribute added to USE-statements
 !
 ! Former revisions:
@@ -4924,7 +4924,11 @@
 !------------------------------------------------------------------------------!
 
-    USE control_parameters
+    USE control_parameters,                                                    &
+        ONLY:  netcdf_data_format
+
     USE netcdf
+
     USE netcdf_control
+
     USE pegrid
 
@@ -4962,7 +4966,11 @@
 !------------------------------------------------------------------------------!
 
-    USE control_parameters
+    USE control_parameters,                                                    &
+        ONLY:  message_string
+
     USE netcdf
+
     USE netcdf_control
+
     USE pegrid
 
Index: /palm/trunk/SOURCE/poismg.f90
===================================================================
--- /palm/trunk/SOURCE/poismg.f90	(revision 1321)
+++ /palm/trunk/SOURCE/poismg.f90	(revision 1322)
@@ -23,5 +23,5 @@
 ! Current revisions:
 ! -----------------
-!
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -505,5 +505,5 @@
                                         ( r(k,j,i) - r(k-1,j+1,i+1) )
 
-             f_mg(kc,jc,ic) = 1.0 / 64.0 * (                            &
+             f_mg(kc,jc,ic) = 1.0 / 64.0_wp * (                         &
                               8.0 * r(k,j,i)                            &
                             + 4.0 * ( rkjim   + rkjip   +               &
@@ -521,7 +521,7 @@
                             +       ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
                                       r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
-                                           ) 
-
-!             f_mg(kc,jc,ic) = 1.0 / 64.0 * (                            &
+                                              )
+
+!             f_mg(kc,jc,ic) = 1.0 / 64.0_wp * (                         &
 !                              8.0 * r(k,j,i)                            &
 !                            + 4.0 * ( r(k,j,i-1)     + r(k,j,i+1)     + &
@@ -539,5 +539,5 @@
 !                            +       ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
 !                                      r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
-!                                           ) 
+!                                              )
           ENDDO
        ENDDO
Index: /palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- /palm/trunk/SOURCE/read_var_list.f90	(revision 1321)
+++ /palm/trunk/SOURCE/read_var_list.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! REAL functions provided with KIND-attribute
 !
 ! Former revisions:
@@ -906,5 +906,5 @@
 !-- Calculate the temporal average of vertical profiles, if neccessary
     IF ( average_count_pr /= 0 )  THEN
-       hom_sum = hom_sum / REAL( average_count_pr )
+       hom_sum = hom_sum / REAL( average_count_pr, KIND=wp )
     ENDIF
 
Index: /palm/trunk/SOURCE/surface_coupler.f90
===================================================================
--- /palm/trunk/SOURCE/surface_coupler.f90	(revision 1321)
+++ /palm/trunk/SOURCE/surface_coupler.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! REAL constants defined as wp-kind
 ! 
 ! Former revisions:
@@ -386,5 +386,5 @@
 !--       * latent heat of vaporization in m2/s2, or 540 cal/g, or 40.65 kJ/mol
 !--       /(rho_atm(=1.0)*c_p)
-          tswst = tswst + qswst_remote * 2.2626108E6 / 1005.0
+          tswst = tswst + qswst_remote * 2.2626108E6_wp / 1005.0_wp
 !
 !--        ...and convert it to a salinity flux at the sea surface (top)
@@ -399,5 +399,5 @@
 !--    (constants are the specific heat capacities for air and water)
 !--    now tswst is the ocean top heat flux
-       tswst = tswst / rho(nzt,:,:) * 1005.0 / 4218.0
+       tswst = tswst / rho(nzt,:,:) * 1005.0_wp / 4218.0_wp
 
 !
Index: /palm/trunk/SOURCE/temperton_fft.f90
===================================================================
--- /palm/trunk/SOURCE/temperton_fft.f90	(revision 1321)
+++ /palm/trunk/SOURCE/temperton_fft.f90	(revision 1322)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL constants defined as wp-kind
 ! 
 ! Former revisions:
@@ -1174,5 +1174,5 @@
     je = jd + 2*m*inc2
     z = 1.0_wp/REAL(n)
-    zsin45 = z*SQRT(0.5)
+    zsin45 = z*SQRT(0.5_wp)
 
     DO l = 1, la
@@ -1628,5 +1628,5 @@
 80  CONTINUE
     ibase = 0
-    sin45 = SQRT(0.5)
+    sin45 = SQRT(0.5_wp)
     DO l = 1, la
        i = ibase
Index: /palm/trunk/SOURCE/timestep.f90
===================================================================
--- /palm/trunk/SOURCE/timestep.f90	(revision 1321)
+++ /palm/trunk/SOURCE/timestep.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! REAL functions provided with KIND-attribute
 ! 
 ! Former revisions:
@@ -156,6 +156,6 @@
              v_gtrans = v_gtrans + hom(k,1,2,0)
           ENDDO
-          u_gtrans = u_gtrans / REAL( nzt - nzb )
-          v_gtrans = v_gtrans / REAL( nzt - nzb )
+          u_gtrans = u_gtrans / REAL( nzt - nzb, KIND=wp )
+          v_gtrans = v_gtrans / REAL( nzt - nzb, KIND=wp )
        ELSE
 !
@@ -173,12 +173,12 @@
           ENDDO
           !$acc end parallel
-          uv_gtrans_l(1) = u_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb) )
-          uv_gtrans_l(2) = v_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb) )
+          uv_gtrans_l(1) = u_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb), KIND=wp )
+          uv_gtrans_l(2) = v_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb), KIND=wp )
 #if defined( __parallel )
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
           CALL MPI_ALLREDUCE( uv_gtrans_l, uv_gtrans, 2, MPI_REAL, MPI_SUM, &
                               comm2d, ierr )
-          u_gtrans = uv_gtrans(1) / REAL( numprocs )
-          v_gtrans = uv_gtrans(2) / REAL( numprocs )
+          u_gtrans = uv_gtrans(1) / REAL( numprocs, KIND=wp )
+          v_gtrans = uv_gtrans(2) / REAL( numprocs, KIND=wp )
 #else
           u_gtrans = uv_gtrans_l(1)
Index: /palm/trunk/SOURCE/timestep_scheme_steering.f90
===================================================================
--- /palm/trunk/SOURCE/timestep_scheme_steering.f90	(revision 1321)
+++ /palm/trunk/SOURCE/timestep_scheme_steering.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! REAL constants defined as wp-kind
 ! 
 ! Former revisions:
@@ -48,4 +48,6 @@
         ONLY:  intermediate_timestep_count, timestep_scheme, tsc
 
+    USE kinds
+
     IMPLICIT NONE
 
@@ -57,15 +59,15 @@
        IF ( timestep_scheme == 'runge-kutta-2' )  THEN
           IF ( intermediate_timestep_count == 1 )  THEN
-             tsc(1:5) = (/ 1.0, 1.0,  0.0, 0.0, 0.0 /)
+             tsc(1:5) = (/ 1.0_wp, 1.0_wp,  0.0_wp, 0.0_wp, 0.0_wp /)
           ELSE
-             tsc(1:5) = (/ 1.0, 0.5, -0.5, 0.0, 1.0 /)
+             tsc(1:5) = (/ 1.0_wp, 0.5_wp, -0.5_wp, 0.0_wp, 1.0_wp /)
           ENDIF
        ELSE
           IF ( intermediate_timestep_count == 1 )  THEN
-             tsc(1:5) = (/ 1.0,   1.0/3.0,        0.0, 0.0, 0.0 /)
+             tsc(1:5) = (/ 1.0_wp,  1.0_wp /  3.0_wp,           0.0_wp, 0.0_wp, 0.0_wp /)
           ELSEIF ( intermediate_timestep_count == 2 )  THEN
-             tsc(1:5) = (/ 1.0, 15.0/16.0, -25.0/48.0, 0.0, 0.0 /)
+             tsc(1:5) = (/ 1.0_wp, 15.0_wp / 16.0_wp, -25.0_wp/48.0_wp, 0.0_wp, 0.0_wp /)
           ELSE
-             tsc(1:5) = (/ 1.0,  8.0/15.0,   1.0/15.0, 0.0, 1.0 /)
+             tsc(1:5) = (/ 1.0_wp,  8.0_wp / 15.0_wp,   1.0_wp/15.0_wp, 0.0_wp, 1.0_wp /)
           ENDIF          
        ENDIF
@@ -74,5 +76,5 @@
 !
 !--    Euler scheme
-       tsc(1:5) = (/ 1.0, 1.0, 0.0, 0.0, 1.0 /)
+       tsc(1:5) = (/ 1.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 1.0_wp /)
 
     ENDIF
Index: /palm/trunk/SOURCE/tridia_solver.f90
===================================================================
--- /palm/trunk/SOURCE/tridia_solver.f90	(revision 1321)
+++ /palm/trunk/SOURCE/tridia_solver.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! REAL functions provided with KIND-attribute
 ! 
 ! Former revisions:
@@ -170,24 +170,24 @@
                    IF ( i >= 0  .AND.  i <= nnxh )  THEN
                       ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
-                                                  REAL( nx+1 ) ) ) / ( dx * dx ) + &
+                                                  REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
                                 2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
-                                                  REAL( ny+1 ) ) ) / ( dy * dy )
+                                                  REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
                    ELSE
                       ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
-                                                  REAL( nx+1 ) ) ) / ( dx * dx ) + &
+                                                  REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
                                 2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
-                                                  REAL( ny+1 ) ) ) / ( dy * dy )
+                                                  REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
                    ENDIF
                 ELSE
                    IF ( i >= 0  .AND.  i <= nnxh )  THEN
                       ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
-                                                  REAL( nx+1 ) ) ) / ( dx * dx ) + &
+                                                  REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
                                 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
-                                                  REAL( ny+1 ) ) ) / ( dy * dy )
+                                                  REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
                    ELSE
                       ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
-                                                  REAL( nx+1 ) ) ) / ( dx * dx ) + &
+                                                  REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
                                 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
-                                                  REAL( ny+1 ) ) ) / ( dy * dy )
+                                                  REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
                    ENDIF
                 ENDIF
@@ -558,12 +558,12 @@
              IF ( i >= 0 .AND. i <= nnxh ) THEN
                 l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
-                                         REAL( nx+1 ) ) ) * ddx2 + &
+                                         REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
                        2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
-                                         REAL( ny+1 ) ) ) * ddy2
+                                         REAL( ny+1, KIND=wp ) ) ) * ddy2
              ELSE
                 l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
-                                         REAL( nx+1 ) ) ) * ddx2 + &
+                                         REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
                        2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
-                                         REAL( ny+1 ) ) ) * ddy2
+                                         REAL( ny+1, KIND=wp ) ) ) * ddy2
              ENDIF
           ENDDO
Index: /palm/trunk/SOURCE/user_3d_data_averaging.f90
===================================================================
--- /palm/trunk/SOURCE/user_3d_data_averaging.f90	(revision 1321)
+++ /palm/trunk/SOURCE/user_3d_data_averaging.f90	(revision 1322)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL functions provided with KIND-attribute
 ! 
 ! Former revisions:
@@ -118,5 +118,5 @@
 !                DO  j = nysg, nyng
 !                   DO  k = nzb, nzt+1
-!                      u2_av(k,j,i) = u2_av(k,j,i) / REAL( average_count_3d )
+!                      u2_av(k,j,i) = u2_av(k,j,i) / REAL( average_count_3d, KIND=wp )
 !                   ENDDO
 !                ENDDO
