Index: /palm/trunk/SOURCE/calc_spectra.f90
===================================================================
--- /palm/trunk/SOURCE/calc_spectra.f90	(revision 1341)
+++ /palm/trunk/SOURCE/calc_spectra.f90	(revision 1342)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! -----------------
+! REAL constants defined as wp-kinds
 !
 ! Former revisions:
@@ -315,5 +316,5 @@
 !
 !-- Exponent for geometric average
-    exponent = 1.0 / ( ny + 1.0 )
+    exponent = 1.0_wp / ( ny + 1.0_wp )
 
 !
@@ -346,5 +347,5 @@
           DO  i = 0, nx/2
 
-             sums_spectra_l(i) = 1.0
+             sums_spectra_l(i) = 1.0_wp
              DO  j = nys_x, nyn_x
                 sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent
@@ -355,5 +356,5 @@
        ELSE
 
-          sums_spectra_l = 1.0
+          sums_spectra_l = 1.0_wp
 
        ENDIF
@@ -361,5 +362,5 @@
 !
 !--    Global sum of spectra on PE0 (from where they are written on file) 
-       sums_spectra(:,n) = 0.0
+       sums_spectra(:,n) = 0.0_wp
 #if defined( __parallel )   
        CALL MPI_BARRIER( comm2d, ierr )  ! Necessary?
@@ -380,7 +381,7 @@
 !--    Temperton fft results need to be normalized
        IF ( fft_method == 'temperton-algorithm' )  THEN
-          fac = nx + 1.0
+          fac = nx + 1.0_wp
        ELSE
-          fac = 1.0
+          fac = 1.0_wp
        ENDIF
        DO  i = 1, nx/2
@@ -454,5 +455,5 @@
 !
 !-- Exponent for geometric average
-    exponent = 1.0 / ( nx + 1.0 )
+    exponent = 1.0_wp / ( nx + 1.0_wp )
 
 !
@@ -485,5 +486,5 @@
           DO  j = 0, ny/2
 
-             sums_spectra_l(j) = 1.0
+             sums_spectra_l(j) = 1.0_wp
              DO  i = nxl_yd, nxr_yd
                 sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent
@@ -494,5 +495,5 @@
        ELSE
 
-          sums_spectra_l = 1.0
+          sums_spectra_l = 1.0_wp
 
        ENDIF
@@ -500,5 +501,5 @@
 !
 !--    Global sum of spectra on PE0 (from where they are written on file) 
-       sums_spectra(:,n) = 0.0
+       sums_spectra(:,n) = 0.0_wp
 #if defined( __parallel )   
        CALL MPI_BARRIER( comm2d, ierr )  ! Necessary?
@@ -520,7 +521,7 @@
 !--    Temperton fft results need to be normalized
        IF ( fft_method == 'temperton-algorithm' )  THEN
-          fac = ny + 1.0
+          fac = ny + 1.0_wp
        ELSE
-          fac = 1.0
+          fac = 1.0_wp
        ENDIF
        DO  j = 1, ny/2
Index: /palm/trunk/SOURCE/fft_xy.f90
===================================================================
--- /palm/trunk/SOURCE/fft_xy.f90	(revision 1341)
+++ /palm/trunk/SOURCE/fft_xy.f90	(revision 1342)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -242,6 +242,6 @@
        IF ( fft_method == 'system-specific' )  THEN
 
-          dnx = 1.0 / ( nx + 1.0 )
-          dny = 1.0 / ( ny + 1.0 )
+          dnx = 1.0_wp / ( nx + 1.0_wp )
+          dny = 1.0_wp / ( ny + 1.0_wp )
           sqr_dnx = SQRT( dnx )
           sqr_dny = SQRT( dny )
@@ -267,6 +267,6 @@
                     trig_yb(2*(ny+1)), trig_yf(2*(ny+1)) )
 
-          work_x = 0.0
-          work_y = 0.0
+          work_x = 0.0_wp
+          work_y = 0.0_wp
           nz1  = nz + MOD( nz+1, 2 )  ! odd nz slows down fft significantly
                                       ! when using the NEC ffts
@@ -433,10 +433,10 @@
                 DO  j = nys_x, nyn_x
 
-                   cwork(0) = CMPLX( ar(0,j,k), 0.0 )
+                   cwork(0) = CMPLX( ar(0,j,k), 0.0_wp )
                    DO  i = 1, (nx+1)/2 - 1
                       cwork(i)      = CMPLX( ar(i,j,k), -ar(nx+1-i,j,k) )
                       cwork(nx+1-i) = CMPLX( ar(i,j,k),  ar(nx+1-i,j,k) )
                    ENDDO
-                   cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
+                   cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp )
 
                    ishape = SHAPE( cwork )
@@ -494,6 +494,6 @@
                       work(2*i+1) = ar(nx+1-i,j,k)
                    ENDDO
-                   work(1)    = 0.0
-                   work(nx+2) = 0.0
+                   work(1)    = 0.0_wp
+                   work(nx+2) = 0.0_wp
 
                    CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
@@ -551,17 +551,17 @@
                    IF ( PRESENT( ar_2d ) )  THEN
 
-                      x_out(0) = CMPLX( ar_2d(0,j), 0.0 )
+                      x_out(0) = CMPLX( ar_2d(0,j), 0.0_wp )
                       DO  i = 1, (nx+1)/2 - 1
                          x_out(i) = CMPLX( ar_2d(i,j), ar_2d(nx+1-i,j) )
                       ENDDO
-                      x_out((nx+1)/2) = CMPLX( ar_2d((nx+1)/2,j), 0.0 )
+                      x_out((nx+1)/2) = CMPLX( ar_2d((nx+1)/2,j), 0.0_wp )
 
                    ELSE
 
-                      x_out(0) = CMPLX( ar(0,j,k), 0.0 )
+                      x_out(0) = CMPLX( ar(0,j,k), 0.0_wp )
                       DO  i = 1, (nx+1)/2 - 1
                          x_out(i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
                       ENDDO
-                      x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
+                      x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp )
 
                    ENDIF
@@ -614,6 +614,6 @@
                       work(2*i+1) = ar(nx+1-i,j,k)
                    ENDDO
-                   work(1) = 0.0
-                   work(nx+2) = 0.0
+                   work(1) = 0.0_wp
+                   work(nx+2) = 0.0_wp
 
                    CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx,      & 
@@ -667,6 +667,6 @@
                       work(2*i+1) = ar(nx+1-i,j,k)
                    ENDDO
-                   work(1) = 0.0
-                   work(nx+2) = 0.0
+                   work(1) = 0.0_wp
+                   work(nx+2) = 0.0_wp
 
                    CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
@@ -711,10 +711,10 @@
                 DO  j = nys_x, nyn_x
 
-                   ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0 )
+                   ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0_wp )
 
                    DO  i = 1, (nx+1)/2 - 1
                       ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
                    ENDDO
-                   ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
+                   ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp )
 
                 ENDDO
@@ -805,10 +805,10 @@
           ELSE
 
-             cwork(0) = CMPLX( ar(0), 0.0 )
+             cwork(0) = CMPLX( ar(0), 0.0_wp )
              DO  i = 1, (nx+1)/2 - 1
                 cwork(i)      = CMPLX( ar(i), -ar(nx+1-i) )
                 cwork(nx+1-i) = CMPLX( ar(i),  ar(nx+1-i) )
              ENDDO
-             cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
+             cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0_wp )
 
              ishape = SHAPE( cwork )
@@ -848,6 +848,6 @@
                 work(2*i+1) = ar(nx+1-i)
              ENDDO
-             work(1)    = 0.0
-             work(nx+2) = 0.0
+             work(1)    = 0.0_wp
+             work(nx+2) = 0.0_wp
 
              CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
@@ -873,9 +873,9 @@
          ELSE
 
-             x_out(0) = CMPLX( ar(0), 0.0 )
+             x_out(0) = CMPLX( ar(0), 0.0_wp )
              DO  i = 1, (nx+1)/2 - 1
                 x_out(i) = CMPLX( ar(i), ar(nx+1-i) )
              ENDDO
-             x_out((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
+             x_out((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0_wp )
 
              CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in)
@@ -908,6 +908,6 @@
                 work(2*i+1) = ar(nx+1-i)
              ENDDO
-             work(1) = 0.0
-             work(nx+2) = 0.0
+             work(1) = 0.0_wp
+             work(nx+2) = 0.0_wp
 
              CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
@@ -941,6 +941,6 @@
                 work(2*i+1) = ar(nx+1-i)
              ENDDO
-             work(1) = 0.0
-             work(nx+2) = 0.0
+             work(1) = 0.0_wp
+             work(nx+2) = 0.0_wp
 
              CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
@@ -1076,10 +1076,10 @@
                 DO  i = nxl_y_l, nxr_y_l
 
-                   cwork(0) = CMPLX( ar_tr(0,i,k), 0.0 )
+                   cwork(0) = CMPLX( ar_tr(0,i,k), 0.0_wp )
                    DO  j = 1, (ny+1)/2 - 1
                       cwork(j)      = CMPLX( ar_tr(j,i,k), -ar_tr(ny+1-j,i,k) )
                       cwork(ny+1-j) = CMPLX( ar_tr(j,i,k),  ar_tr(ny+1-j,i,k) )
                    ENDDO
-                   cwork((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0 )
+                   cwork((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp )
 
                    jshape = SHAPE( cwork )
@@ -1137,6 +1137,6 @@
                       work(2*j+1) = ar_tr(ny+1-j,i,k)
                    ENDDO
-                   work(1)    = 0.0
-                   work(ny+2) = 0.0
+                   work(1)    = 0.0_wp
+                   work(ny+2) = 0.0_wp
 
                    CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
@@ -1180,9 +1180,9 @@
                 DO  i = nxl_y_l, nxr_y_l
 
-                   y_out(0) = CMPLX( ar_tr(0,i,k), 0.0 )
+                   y_out(0) = CMPLX( ar_tr(0,i,k), 0.0_wp )
                    DO  j = 1, (ny+1)/2 - 1
                       y_out(j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k) )
                    ENDDO
-                   y_out((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0 )
+                   y_out((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp )
 
                    CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
@@ -1233,6 +1233,6 @@
                       work(2*j+1) = ar_tr(ny+1-j,i,k)
                    ENDDO
-                   work(1)    = 0.0
-                   work(ny+2) = 0.0
+                   work(1)    = 0.0_wp
+                   work(ny+2) = 0.0_wp
 
                    CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny,      &
@@ -1284,6 +1284,6 @@
                       work(2*j+1) = ar_tr(ny+1-j,i,k)
                    ENDDO
-                   work(1) = 0.0
-                   work(ny+2) = 0.0
+                   work(1) = 0.0_wp
+                   work(ny+2) = 0.0_wp
 
                    CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
@@ -1327,10 +1327,10 @@
                 DO  i = nxl_y, nxr_y
 
-                   ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0 )
+                   ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0_wp )
 
                    DO  j = 1, (ny+1)/2 - 1
                       ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
                    ENDDO
-                   ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
+                   ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0_wp )
 
                 ENDDO
@@ -1423,10 +1423,10 @@
           ELSE
 
-             cwork(0) = CMPLX( ar(0), 0.0 )
+             cwork(0) = CMPLX( ar(0), 0.0_wp )
              DO  j = 1, (ny+1)/2 - 1
                 cwork(j)      = CMPLX( ar(j), -ar(ny+1-j) )
                 cwork(ny+1-j) = CMPLX( ar(j),  ar(ny+1-j) )
              ENDDO
-             cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )
+             cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0_wp )
 
              jshape = SHAPE( cwork )
@@ -1466,6 +1466,6 @@
                 work(2*j+1) = ar(ny+1-j)
              ENDDO
-             work(1)    = 0.0
-             work(ny+2) = 0.0
+             work(1)    = 0.0_wp
+             work(ny+2) = 0.0_wp
 
              CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
@@ -1491,9 +1491,9 @@
           ELSE
 
-             y_out(0) = CMPLX( ar(0), 0.0 )
+             y_out(0) = CMPLX( ar(0), 0.0_wp )
              DO  j = 1, (ny+1)/2 - 1
                 y_out(j) = CMPLX( ar(j), ar(ny+1-j) )
              ENDDO
-             y_out((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )
+             y_out((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0_wp )
 
              CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
@@ -1526,6 +1526,6 @@
                 work(2*j+1) = ar(ny+1-j)
              ENDDO
-             work(1)    = 0.0
-             work(ny+2) = 0.0
+             work(1)    = 0.0_wp
+             work(ny+2) = 0.0_wp
 
              CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3,      &
@@ -1559,6 +1559,6 @@
                 work(2*j+1) = ar(ny+1-j)
              ENDDO
-             work(1) = 0.0
-             work(ny+2) = 0.0
+             work(1) = 0.0_wp
+             work(ny+2) = 0.0_wp
 
              CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
@@ -1621,5 +1621,5 @@
 
              ai(0:nx,1:nz) = ar(0:nx,1:nz)
-             ai(nx+1:,:)   = 0.0
+             ai(nx+1:,:)   = 0.0_wp
 
              CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
@@ -1643,6 +1643,6 @@
                    ai(2*i+1,k) = ar(nx+1-i,k)
                 ENDDO
-                ai(1,k) = 0.0
-                ai(nx+2,k) = 0.0
+                ai(1,k) = 0.0_wp
+                ai(nx+2,k) = 0.0_wp
              ENDDO
 
@@ -1669,5 +1669,5 @@
              ai(0:nx,1:nz) = ar(0:nx,1:nz)
              IF ( nz1 > nz )  THEN
-                ai(:,nz1) = 0.0
+                ai(:,nz1) = 0.0_wp
              ENDIF
 
@@ -1693,12 +1693,12 @@
 
              IF ( nz1 > nz )  THEN
-                work(:,nz1) = 0.0
+                work(:,nz1) = 0.0_wp
              ENDIF
              DO  k = 1, nz
-                work(1,k) = CMPLX( ar(0,k), 0.0 )
+                work(1,k) = CMPLX( ar(0,k), 0.0_wp )
                 DO  i = 1, (nx+1)/2 - 1
                    work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) )
                 ENDDO
-                work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )
+                work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0_wp )
              ENDDO
 
@@ -1764,5 +1764,5 @@
 
              ai(0:ny,1:nz) = ar(0:ny,1:nz)
-             ai(ny+1:,:)   = 0.0
+             ai(ny+1:,:)   = 0.0_wp
 
              CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
@@ -1786,6 +1786,6 @@
                    ai(2*j+1,k) = ar(ny+1-j,k)
                 ENDDO
-                ai(1,k) = 0.0
-                ai(ny+2,k) = 0.0
+                ai(1,k) = 0.0_wp
+                ai(ny+2,k) = 0.0_wp
              ENDDO
 
@@ -1812,5 +1812,5 @@
              ai(0:ny,1:nz) = ar(0:ny,1:nz)
              IF ( nz1 > nz )  THEN
-                ai(:,nz1) = 0.0
+                ai(:,nz1) = 0.0_wp
              ENDIF
 
@@ -1836,12 +1836,12 @@
 
              IF ( nz1 > nz )  THEN
-                work(:,nz1) = 0.0
+                work(:,nz1) = 0.0_wp
              ENDIF
              DO  k = 1, nz
-                work(1,k) = CMPLX( ar(0,k), 0.0 )
+                work(1,k) = CMPLX( ar(0,k), 0.0_wp )
                 DO  j = 1, (ny+1)/2 - 1
                    work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) )
                 ENDDO
-                work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )
+                work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0_wp )
              ENDDO
 
Index: /palm/trunk/SOURCE/pres.f90
===================================================================
--- /palm/trunk/SOURCE/pres.f90	(revision 1341)
+++ /palm/trunk/SOURCE/pres.f90	(revision 1342)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -135,6 +135,6 @@
 
 
-    ddt_3d = 1.0 / dt_3d
-    d_weight_pres = 1.0 / weight_pres(intermediate_timestep_count)
+    ddt_3d = 1.0_wp / dt_3d
+    d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count)
 
 !
@@ -180,6 +180,6 @@
     IF ( conserve_volume_flow  .AND.  ( outflow_l .OR. outflow_r ) )  THEN
 
-       volume_flow(1)   = 0.0
-       volume_flow_l(1) = 0.0
+       volume_flow(1)   = 0.0_wp
+       volume_flow_l(1) = 0.0_wp
 
        IF ( outflow_l )  THEN
@@ -219,6 +219,6 @@
     IF ( conserve_volume_flow  .AND.  ( outflow_n .OR. outflow_s ) )  THEN
 
-       volume_flow(2)   = 0.0
-       volume_flow_l(2) = 0.0
+       volume_flow(2)   = 0.0_wp
+       volume_flow_l(2) = 0.0_wp
 
        IF ( outflow_s )  THEN
@@ -257,6 +257,6 @@
 !-- Remove mean vertical velocity
     IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
-       IF ( simulated_time > 0.0 )  THEN ! otherwise nzb_w_inner not yet known
-          w_l = 0.0;  w_l_l = 0.0
+       IF ( simulated_time > 0.0_wp )  THEN ! otherwise nzb_w_inner not yet known
+          w_l = 0.0_wp;  w_l_l = 0.0_wp
           DO  i = nxl, nxr
              DO  j = nys, nyn
@@ -295,5 +295,5 @@
           DO  j = nys-1, nyn+1
              DO  k = nzb, nzt+1
-                d(k,j,i) = 0.0
+                d(k,j,i) = 0.0_wp
              ENDDO
           ENDDO
@@ -305,5 +305,5 @@
           DO  j = nys, nyn
              DO  k = nzb+1, nzt
-                d(k,j,i) = 0.0
+                d(k,j,i) = 0.0_wp
              ENDDO
           ENDDO
@@ -312,6 +312,6 @@
     ENDIF
 
-    localsum  = 0.0
-    threadsum = 0.0
+    localsum  = 0.0_wp
+    threadsum = 0.0_wp
 
 #if defined( __ibm )
@@ -432,5 +432,5 @@
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                tend(nzb_s_inner(j,i),j,i) = 0.0
+                tend(nzb_s_inner(j,i),j,i) = 0.0_wp
              ENDDO
           ENDDO
@@ -460,5 +460,5 @@
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                tend(nzt+1,j,i) = 0.0
+                tend(nzt+1,j,i) = 0.0_wp
              ENDDO
           ENDDO
@@ -577,6 +577,6 @@
 !-- pressure just computed
     IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc .OR. bc_ns_cyc ) )  THEN
-       volume_flow_l(1) = 0.0
-       volume_flow_l(2) = 0.0
+       volume_flow_l(1) = 0.0_wp
+       volume_flow_l(2) = 0.0_wp
     ENDIF
 
@@ -704,13 +704,13 @@
 !-- a possible PE-sum is computed in flow_statistics 
     CALL cpu_log( log_point_s(1), 'divergence', 'start' )
-    sums_divnew_l = 0.0
+    sums_divnew_l = 0.0_wp
 
 !
 !-- d must be reset to zero because it can contain nonzero values below the
 !-- topography
-    IF ( topography /= 'flat' )  d = 0.0
-
-    localsum  = 0.0
-    threadsum = 0.0
+    IF ( topography /= 'flat' )  d = 0.0_wp
+
+    localsum  = 0.0_wp
+    threadsum = 0.0_wp
 
     !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
Index: /palm/trunk/SOURCE/production_e.f90
===================================================================
--- /palm/trunk/SOURCE/production_e.f90	(revision 1341)
+++ /palm/trunk/SOURCE/production_e.f90	(revision 1342)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -173,27 +173,27 @@
              DO  k = nzb_diff_s_outer(j,i), nzt
 
-                dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-                dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+                dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+                dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                       dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                IF ( def < 0.0 )  def = 0.0
+                      dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                IF ( def < 0.0_wp )  def = 0.0_wp
 
                 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -210,17 +210,17 @@
              DO  j = nys, nyn
 
-                IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
+                IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) &
                 THEN
 
                    k = nzb_diff_s_inner(j,i) - 1
                    dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
-                   dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                  u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
+                   dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                     u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
                    dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
-                   dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                  v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
+                   dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                     v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
                    dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
 
-                   IF ( wall_e_y(j,i) /= 0.0 )  THEN
+                   IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
 !
 !--                   Inconsistency removed: as the thermal stratification is 
@@ -236,21 +236,21 @@
                       CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                           wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
-                      km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25 * &
-                                   0.5 * dy 
-                      IF ( km_neutral > 0.0 )  THEN
+                      km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25_wp * &
+                                   0.5_wp * dy 
+                      IF ( km_neutral > 0.0_wp )  THEN
                          dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
                          dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
                       ELSE
-                         dudy = 0.0
-                         dwdy = 0.0
+                         dudy = 0.0_wp
+                         dwdy = 0.0_wp
                       ENDIF
                    ELSE
-                      dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                      u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                      dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                      w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                      dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                         u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                      dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                         w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
                    ENDIF
 
-                   IF ( wall_e_x(j,i) /= 0.0 )  THEN
+                   IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
 !
 !--                   Inconsistency removed: as the thermal stratification is 
@@ -266,25 +266,25 @@
                       CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                           wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
-                      km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25 * &
-                                   0.5 * dx
-                      IF ( km_neutral > 0.0 )  THEN
+                      km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25_wp * &
+                                   0.5_wp * dx
+                      IF ( km_neutral > 0.0_wp )  THEN
                          dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
                          dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
                       ELSE
-                         dvdx = 0.0
-                         dwdx = 0.0
+                         dvdx = 0.0_wp
+                         dwdx = 0.0_wp
                       ENDIF
                    ELSE
-                      dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                      v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                      dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                      w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                      dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                         v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                      dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                         w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
                    ENDIF
 
-                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                          dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                   IF ( def < 0.0 )  def = 0.0
+                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                   IF ( def < 0.0_wp )  def = 0.0_wp
 
                    tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -301,12 +301,12 @@
 
                       dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
-                      dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                     u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-                      dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                      dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+                      dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                        u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+                      dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                      dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                        v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
                       dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
 
-                      IF ( wall_e_y(j,i) /= 0.0 )  THEN
+                      IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
 !
 !--                      Inconsistency removed: as the thermal stratification 
@@ -319,20 +319,20 @@
 !--                            validation has been available
                          km_neutral = kappa * ( usvs(k)**2 + & 
-                                                wsvs(k)**2 )**0.25 * 0.5 * dy
-                         IF ( km_neutral > 0.0 )  THEN
+                                                wsvs(k)**2 )**0.25_wp * 0.5_wp * dy
+                         IF ( km_neutral > 0.0_wp )  THEN
                             dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
                             dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
                          ELSE
-                            dudy = 0.0
-                            dwdy = 0.0
+                            dudy = 0.0_wp
+                            dwdy = 0.0_wp
                          ENDIF
                       ELSE
-                         dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                         u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                         dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                         w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                         dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                            u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                         dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                            w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
                       ENDIF
 
-                      IF ( wall_e_x(j,i) /= 0.0 )  THEN
+                      IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
 !
 !--                      Inconsistency removed: as the thermal stratification 
@@ -345,24 +345,24 @@
 !--                            validation has been available
                          km_neutral = kappa * ( vsus(k)**2 + & 
-                                                wsus(k)**2 )**0.25 * 0.5 * dx
-                         IF ( km_neutral > 0.0 )  THEN
+                                                wsus(k)**2 )**0.25_wp * 0.5_wp * dx
+                         IF ( km_neutral > 0.0_wp )  THEN
                             dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
                             dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
                          ELSE
-                            dvdx = 0.0
-                            dwdx = 0.0
+                            dvdx = 0.0_wp
+                            dwdx = 0.0_wp
                          ENDIF
                       ELSE
-                         dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                         v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                         dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                         w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                         dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                            v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                         dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                            w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
                       ENDIF
 
-                      def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                      def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                            dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 +  &
-                           dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                      IF ( def < 0.0 )  def = 0.0
+                           dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                      IF ( def < 0.0_wp )  def = 0.0_wp
 
                       tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -379,32 +379,32 @@
              DO  j = nys, nyn
 
-                IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
+                IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) &
                 THEN
 
                    k = nzb_diff_s_outer(j,i)-1
 
-                   dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                   dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                   dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-                   dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                   dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                   dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                   dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                   dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                   dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                   dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                   dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                   dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                       u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+                   dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                   dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                   dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                       v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+                   dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                   dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                   dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                          dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                   IF ( def < 0.0 )  def = 0.0
+                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                   IF ( def < 0.0_wp )  def = 0.0_wp
 
                    tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -420,68 +420,68 @@
              DO  j = nys, nyn
 
-                IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
+                IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) &
                 THEN
 
                    k = nzb_diff_s_inner(j,i)-1
 
-                   dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                   dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                   dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                   dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                   dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                       u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
+
+                   dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                   dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                   dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                       v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
+
+                   dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                   dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                   dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
+                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                   IF ( def < 0.0_wp )  def = 0.0_wp
+
+                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
+
+                ENDIF
+
+             ENDDO
+
+          ELSEIF ( use_surface_fluxes )  THEN
+
+             DO  j = nys, nyn
+
+                k = nzb_diff_s_outer(j,i)-1
+
+                dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
                                     u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                   dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                    u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
-
-                   dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                   u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+                dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
                                     v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                   dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                   dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                    v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
-
-                   dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+                dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
                                     w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                   dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
                                     w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                   dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
-                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                   IF ( def < 0.0 )  def = 0.0
-
-                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-
-                ENDIF
-
-             ENDDO
-
-          ELSEIF ( use_surface_fluxes )  THEN
-
-             DO  j = nys, nyn
-
-                k = nzb_diff_s_outer(j,i)-1
-
-                dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-                dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                       dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                IF ( def < 0.0 )  def = 0.0
+                      dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                IF ( def < 0.0_wp )  def = 0.0_wp
 
                 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -586,6 +586,6 @@
 
                       IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
-                         k1 = 1.0 + 0.61 * q(k,j,i)
-                         k2 = 0.61 * pt(k,j,i)
+                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                         k2 = 0.61_wp * pt(k,j,i)
                          tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
                                          g / vpt(k,j,i) *                      &
@@ -594,16 +594,16 @@
                                          ) * dd2zu(k)
                       ELSE IF ( cloud_physics )  THEN
-                         IF ( ql(k,j,i) == 0.0 )  THEN
-                            k1 = 1.0 + 0.61 * q(k,j,i)
-                            k2 = 0.61 * pt(k,j,i)
+                         IF ( ql(k,j,i) == 0.0_wp )  THEN
+                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                            k2 = 0.61_wp * pt(k,j,i)
                          ELSE
                             theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
                             temp  = theta * t_d_pt(k)
-                            k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                                       ( q(k,j,i) - ql(k,j,i) ) *          &
-                                 ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                                 ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
+                                          ( q(k,j,i) - ql(k,j,i) ) *          &
+                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
+                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
                                  ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                            k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
                          ENDIF
                          tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
@@ -613,6 +613,6 @@
                                          ) * dd2zu(k)
                       ELSE IF ( cloud_droplets )  THEN
-                         k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-                         k2 = 0.61 * pt(k,j,i)
+                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                         k2 = 0.61_wp * pt(k,j,i)
                          tend(k,j,i) = tend(k,j,i) -                          & 
                                        kh(k,j,i) * g / vpt(k,j,i) *           &
@@ -634,23 +634,23 @@
 
                       IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
-                         k1 = 1.0 + 0.61 * q(k,j,i)
-                         k2 = 0.61 * pt(k,j,i)
+                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                         k2 = 0.61_wp * pt(k,j,i)
                       ELSE IF ( cloud_physics )  THEN
-                         IF ( ql(k,j,i) == 0.0 )  THEN
-                            k1 = 1.0 + 0.61 * q(k,j,i)
-                            k2 = 0.61 * pt(k,j,i)
+                         IF ( ql(k,j,i) == 0.0_wp )  THEN
+                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                            k2 = 0.61_wp * pt(k,j,i)
                          ELSE
                             theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
                             temp  = theta * t_d_pt(k)
-                            k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                                       ( q(k,j,i) - ql(k,j,i) ) *          &
-                                 ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                                 ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
+                                          ( q(k,j,i) - ql(k,j,i) ) *          &
+                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
+                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
                                  ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                            k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
                          ENDIF
                       ELSE IF ( cloud_droplets )  THEN
-                         k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-                         k2 = 0.61 * pt(k,j,i)
+                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                         k2 = 0.61_wp * pt(k,j,i)
                       ENDIF
 
@@ -668,23 +668,23 @@
 
                       IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
-                         k1 = 1.0 + 0.61 * q(k,j,i)
-                         k2 = 0.61 * pt(k,j,i)
+                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                         k2 = 0.61_wp * pt(k,j,i)
                       ELSE IF ( cloud_physics )  THEN
-                         IF ( ql(k,j,i) == 0.0 )  THEN
-                            k1 = 1.0 + 0.61 * q(k,j,i)
-                            k2 = 0.61 * pt(k,j,i)
+                         IF ( ql(k,j,i) == 0.0_wp )  THEN
+                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                            k2 = 0.61_wp * pt(k,j,i)
                          ELSE
                             theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
                             temp  = theta * t_d_pt(k)
-                            k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                                       ( q(k,j,i) - ql(k,j,i) ) *          &
-                                 ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                                 ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
+                                          ( q(k,j,i) - ql(k,j,i) ) *          &
+                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
+                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
                                  ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                            k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
                          ENDIF
                       ELSE IF ( cloud_droplets )  THEN
-                         k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-                         k2 = 0.61 * pt(k,j,i)
+                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                         k2 = 0.61_wp * pt(k,j,i)
                       ENDIF
 
@@ -783,27 +783,27 @@
                 IF ( k >= nzb_diff_s_outer(j,i) )  THEN
 
-                   dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                   dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                   dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-                   dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                   dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                   dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                   dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                   dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                   dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                   dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                   dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                   dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                       u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+                   dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                   dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                   dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                       v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+                   dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                   dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                   dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                          dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                   IF ( def < 0.0 )  def = 0.0
+                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                   IF ( def < 0.0_wp )  def = 0.0_wp
 
                    tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -825,17 +825,17 @@
                 DO  k = 1, nzt
 
-                   IF ( ( wall_e_x(j,i) /= 0.0 ).OR.( wall_e_y(j,i) /= 0.0 ) ) &
+                   IF ( ( wall_e_x(j,i) /= 0.0_wp ).OR.( wall_e_y(j,i) /= 0.0_wp ) ) &
                    THEN
 
                       IF ( k == nzb_diff_s_inner(j,i) - 1 )  THEN
                          dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
-                         dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                        u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
+                         dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                           u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
                          dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
-                         dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                        v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
+                         dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                           v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
                          dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
 
-                         IF ( wall_e_y(j,i) /= 0.0 )  THEN
+                         IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
 !
 !--                         Inconsistency removed: as the thermal stratification is
@@ -852,21 +852,21 @@
 !                                                wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
                             km_neutral = kappa *                                    &
-                                        ( usvs(k,j,i)**2 + wsvs(k,j,i)**2 )**0.25 * &
-                                         0.5 * dy
-                            IF ( km_neutral > 0.0 )  THEN
+                                        ( usvs(k,j,i)**2 + wsvs(k,j,i)**2 )**0.25_wp * &
+                                         0.5_wp * dy
+                            IF ( km_neutral > 0.0_wp )  THEN
                                dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral
                                dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral
                             ELSE
-                               dudy = 0.0
-                               dwdy = 0.0
+                               dudy = 0.0_wp
+                               dwdy = 0.0_wp
                             ENDIF
                          ELSE
-                            dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                            u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                            dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                            w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                            dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                               u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                            dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                               w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
                          ENDIF
 
-                         IF ( wall_e_x(j,i) /= 0.0 )  THEN
+                         IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
 !
 !--                         Inconsistency removed: as the thermal stratification is
@@ -883,25 +883,25 @@
 !                                                wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
                             km_neutral = kappa *                                     &
-                                         ( vsus(k,j,i)**2 + wsus(k,j,i)**2 )**0.25 * &
-                                         0.5 * dx
-                            IF ( km_neutral > 0.0 )  THEN
+                                         ( vsus(k,j,i)**2 + wsus(k,j,i)**2 )**0.25_wp * &
+                                         0.5_wp * dx
+                            IF ( km_neutral > 0.0_wp )  THEN
                                dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral
                                dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral
                             ELSE
-                               dvdx = 0.0
-                               dwdx = 0.0
+                               dvdx = 0.0_wp
+                               dwdx = 0.0_wp
                             ENDIF
                          ELSE
-                            dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                            v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                            dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                            w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                            dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                               v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                            dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                               w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
                          ENDIF
 
-                         def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                         def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                                dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                               dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                         IF ( def < 0.0 )  def = 0.0
+                               dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                         IF ( def < 0.0_wp )  def = 0.0_wp
 
                          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -919,12 +919,12 @@
 
                          dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
-                         dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                        u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-                         dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                         dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                        v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+                         dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+                         dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                         dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
                          dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
 
-                         IF ( wall_e_y(j,i) /= 0.0 )  THEN
+                         IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
 !
 !--                         Inconsistency removed: as the thermal stratification
@@ -937,20 +937,20 @@
 !--                               validation has been available
                             km_neutral = kappa * ( usvs(k,j,i)**2 + &
-                                                   wsvs(k,j,i)**2 )**0.25 * 0.5 * dy
-                            IF ( km_neutral > 0.0 )  THEN
+                                                   wsvs(k,j,i)**2 )**0.25_wp * 0.5_wp * dy
+                            IF ( km_neutral > 0.0_wp )  THEN
                                dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral
                                dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral
                             ELSE
-                               dudy = 0.0
-                               dwdy = 0.0
+                               dudy = 0.0_wp
+                               dwdy = 0.0_wp
                             ENDIF
                          ELSE
-                            dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                            u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                            dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                            w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                            dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                               u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                            dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                               w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
                          ENDIF
 
-                         IF ( wall_e_x(j,i) /= 0.0 )  THEN
+                         IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
 !
 !--                         Inconsistency removed: as the thermal stratification
@@ -963,24 +963,24 @@
 !--                               validation has been available
                             km_neutral = kappa * ( vsus(k,j,i)**2 + &
-                                                   wsus(k,j,i)**2 )**0.25 * 0.5 * dx
-                            IF ( km_neutral > 0.0 )  THEN
+                                                   wsus(k,j,i)**2 )**0.25_wp * 0.5_wp * dx
+                            IF ( km_neutral > 0.0_wp )  THEN
                                dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral
                                dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral
                             ELSE
-                               dvdx = 0.0
-                               dwdx = 0.0
+                               dvdx = 0.0_wp
+                               dwdx = 0.0_wp
                             ENDIF
                          ELSE
-                            dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                            v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                            dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                            w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                            dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                               v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                            dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                               w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
                          ENDIF
 
-                         def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                         def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                               dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 +  &
-                              dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                         IF ( def < 0.0 )  def = 0.0
+                              dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                         IF ( def < 0.0_wp )  def = 0.0_wp
 
                          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -993,27 +993,27 @@
                       IF ( k == nzb_diff_s_outer(j,i)-1 )  THEN
 
-                         dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                         dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                          u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                         dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                          u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-                         dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                          v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                         dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                         dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                          v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                         dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                          w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                         dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                          w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                         dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                         def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                         dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                         dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                             u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                         dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                             u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+                         dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                             v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                         dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                         dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                             v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+                         dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                             w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                         dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                             w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                         dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                         def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                                dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                               dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                         IF ( def < 0.0 )  def = 0.0
+                               dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                         IF ( def < 0.0_wp )  def = 0.0_wp
 
                          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -1035,32 +1035,32 @@
                 DO  k = 1, nzt
 
-                   IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
+                   IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) &
                    THEN
 
                       IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
 
-                         dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                         dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                          u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                         dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                          u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
-
-                         dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                          v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                         dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                         dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                          v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
-
-                         dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                          w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                         dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                          w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                         dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                         def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                         dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                         dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                             u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                         dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                             u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
+
+                         dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                             v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                         dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                         dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                             v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
+
+                         dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                             w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                         dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                             w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                         dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                         def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                                dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                               dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                         IF ( def < 0.0 )  def = 0.0
+                               dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                         IF ( def < 0.0_wp )  def = 0.0_wp
 
                          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -1082,27 +1082,27 @@
                    IF ( k == nzb_diff_s_outer(j,i)-1 )  THEN
 
-                      dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-                      dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                      dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                                       u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-                      dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                      dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                      dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                                       v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-                      dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-                      dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-                      dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-                      def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                      dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                      dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                          u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                      dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                          u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+                      dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                          v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                      dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                      dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                          v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+                      dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                          w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                      dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                          w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                      dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                      def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                             dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                            dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                      IF ( def < 0.0 )  def = 0.0
+                            dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                      IF ( def < 0.0_wp )  def = 0.0_wp
 
                       tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -1232,6 +1232,6 @@
 !
 !                         IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
-!                            k1 = 1.0 + 0.61 * q(k,j,i)
-!                            k2 = 0.61 * pt(k,j,i)
+!                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+!                            k2 = 0.61_wp * pt(k,j,i)
 !                            tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
 !                                            g / vpt(k,j,i) *                      &
@@ -1240,16 +1240,16 @@
 !                                            ) * dd2zu(k)
 !                         ELSE IF ( cloud_physics )  THEN
-!                            IF ( ql(k,j,i) == 0.0 )  THEN
-!                               k1 = 1.0 + 0.61 * q(k,j,i)
-!                               k2 = 0.61 * pt(k,j,i)
+!                            IF ( ql(k,j,i) == 0.0_wp )  THEN
+!                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+!                               k2 = 0.61_wp * pt(k,j,i)
 !                            ELSE
 !                               theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
 !                               temp  = theta * t_d_pt(k)
-!                               k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-!                                          ( q(k,j,i) - ql(k,j,i) ) *          &
-!                                    ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-!                                    ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+!                               k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
+!                                             ( q(k,j,i) - ql(k,j,i) ) *          &
+!                                    ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
+!                                    ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
 !                                    ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-!                               k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+!                               k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
 !                            ENDIF
 !                            tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
@@ -1259,6 +1259,6 @@
 !                                            ) * dd2zu(k)
 !                         ELSE IF ( cloud_droplets )  THEN
-!                            k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-!                            k2 = 0.61 * pt(k,j,i)
+!                            k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+!                            k2 = 0.61_wp * pt(k,j,i)
 !                            tend(k,j,i) = tend(k,j,i) -                          &
 !                                          kh(k,j,i) * g / vpt(k,j,i) *           &
@@ -1289,23 +1289,23 @@
 !
 !                            IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
-!                               k1 = 1.0 + 0.61 * q(k,j,i)
-!                               k2 = 0.61 * pt(k,j,i)
+!                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+!                               k2 = 0.61_wp * pt(k,j,i)
 !                            ELSE IF ( cloud_physics )  THEN
-!                               IF ( ql(k,j,i) == 0.0 )  THEN
-!                                  k1 = 1.0 + 0.61 * q(k,j,i)
-!                                  k2 = 0.61 * pt(k,j,i)
+!                               IF ( ql(k,j,i) == 0.0_wp )  THEN
+!                                  k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+!                                  k2 = 0.61_wp * pt(k,j,i)
 !                               ELSE
 !                                  theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
 !                                  temp  = theta * t_d_pt(k)
-!                                  k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-!                                             ( q(k,j,i) - ql(k,j,i) ) *          &
-!                                       ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-!                                       ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+!                                  k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *        &
+!                                                ( q(k,j,i) - ql(k,j,i) ) *    &
+!                                       ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /&
+!                                       ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
 !                                       ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-!                                  k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+!                                  k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
 !                               ENDIF
 !                            ELSE IF ( cloud_droplets )  THEN
-!                               k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-!                               k2 = 0.61 * pt(k,j,i)
+!                               k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+!                               k2 = 0.61_wp * pt(k,j,i)
 !                            ENDIF
 !
@@ -1330,23 +1330,23 @@
 !
 !                            IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
-!                               k1 = 1.0 + 0.61 * q(k,j,i)
-!                               k2 = 0.61 * pt(k,j,i)
+!                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+!                               k2 = 0.61_wp * pt(k,j,i)
 !                            ELSE IF ( cloud_physics )  THEN
-!                               IF ( ql(k,j,i) == 0.0 )  THEN
-!                                  k1 = 1.0 + 0.61 * q(k,j,i)
-!                                  k2 = 0.61 * pt(k,j,i)
+!                               IF ( ql(k,j,i) == 0.0_wp )  THEN
+!                                  k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+!                                  k2 = 0.61_wp * pt(k,j,i)
 !                               ELSE
 !                                  theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
 !                                  temp  = theta * t_d_pt(k)
-!                                  k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-!                                             ( q(k,j,i) - ql(k,j,i) ) *          &
-!                                       ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-!                                       ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+!                                  k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *        &
+!                                                ( q(k,j,i) - ql(k,j,i) ) *    &
+!                                       ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /&
+!                                       ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
 !                                       ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-!                                  k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+!                                  k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
 !                               ENDIF
 !                            ELSE IF ( cloud_droplets )  THEN
-!                               k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-!                               k2 = 0.61 * pt(k,j,i)
+!                               k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+!                               k2 = 0.61_wp * pt(k,j,i)
 !                            ENDIF
 !
@@ -1426,27 +1426,27 @@
        DO  k = nzb_diff_s_outer(j,i), nzt
 
-          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
-                + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
-                + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-          IF ( def < 0.0 )  def = 0.0
+          dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+          dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+          dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+          dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+          dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+          dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+          dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+          dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+          dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+          def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
+                + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2    &
+                + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+          IF ( def < 0.0_wp )  def = 0.0_wp
 
           tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -1456,5 +1456,5 @@
        IF ( prandtl_layer )  THEN
 
-          IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) )  THEN
+          IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) )  THEN
 
 !
@@ -1465,12 +1465,12 @@
 
              dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
-             dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                            u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
-             dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-             dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                            v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
+             dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                               u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
+             dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+             dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                               v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
              dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
 
-             IF ( wall_e_y(j,i) /= 0.0 )  THEN
+             IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
 !
 !--             Inconsistency removed: as the thermal stratification 
@@ -1486,21 +1486,21 @@
                 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                     wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
-                km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25 * &
-                             0.5 * dy
-                IF ( km_neutral > 0.0 )  THEN
+                km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25_wp * &
+                             0.5_wp * dy
+                IF ( km_neutral > 0.0_wp )  THEN
                    dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
                    dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
                 ELSE
-                   dudy = 0.0
-                   dwdy = 0.0
+                   dudy = 0.0_wp
+                   dwdy = 0.0_wp
                 ENDIF
              ELSE
-                dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                   u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                   w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
              ENDIF
 
-             IF ( wall_e_x(j,i) /= 0.0 )  THEN
+             IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
 !
 !--             Inconsistency removed: as the thermal stratification 
@@ -1516,25 +1516,25 @@
                 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
                                     wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
-                km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25 * & 
-                             0.5 * dx
-                IF ( km_neutral > 0.0 )  THEN
+                km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25_wp * & 
+                             0.5_wp * dx
+                IF ( km_neutral > 0.0_wp )  THEN
                    dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
                    dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
                 ELSE
-                   dvdx = 0.0
-                   dwdx = 0.0
+                   dvdx = 0.0_wp
+                   dwdx = 0.0_wp
                 ENDIF
              ELSE
-                dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                   v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                   w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
              ENDIF
 
-             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+             def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                    dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                   dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-             IF ( def < 0.0 )  def = 0.0
+                   dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+             IF ( def < 0.0_wp )  def = 0.0_wp
 
              tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -1549,12 +1549,12 @@
 
                 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
-                dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                               u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-                dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-                dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                               v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+                dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                  u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+                dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                  v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
                 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
 
-                IF ( wall_e_y(j,i) /= 0.0 )  THEN
+                IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
 !
 !--                Inconsistency removed: as the thermal stratification 
@@ -1567,20 +1567,20 @@
 !--                      validation has been available
                    km_neutral = kappa * ( usvs(k)**2 + & 
-                                          wsvs(k)**2 )**0.25 * 0.5 * dy
-                   IF ( km_neutral > 0.0 )  THEN
+                                          wsvs(k)**2 )**0.25_wp * 0.5_wp * dy
+                   IF ( km_neutral > 0.0_wp )  THEN
                       dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
                       dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
                    ELSE
-                      dudy = 0.0
-                      dwdy = 0.0
+                      dudy = 0.0_wp
+                      dwdy = 0.0_wp
                    ENDIF
                 ELSE
-                   dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                                   u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-                   dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                                   w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                   dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                      u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                   dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                      w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
                 ENDIF
 
-                IF ( wall_e_x(j,i) /= 0.0 )  THEN
+                IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
 !
 !--                Inconsistency removed: as the thermal stratification 
@@ -1593,24 +1593,24 @@
 !--                      validation has been available
                    km_neutral = kappa * ( vsus(k)**2 + & 
-                                          wsus(k)**2 )**0.25 * 0.5 * dx
-                   IF ( km_neutral > 0.0 )  THEN
+                                          wsus(k)**2 )**0.25_wp * 0.5_wp * dx
+                   IF ( km_neutral > 0.0_wp )  THEN
                       dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
                       dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
                    ELSE
-                      dvdx = 0.0
-                      dwdx = 0.0
+                      dvdx = 0.0_wp
+                      dwdx = 0.0_wp
                    ENDIF
                 ELSE
-                   dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                                   v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-                   dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                                   w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                   dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                      v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                   dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                      w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
                 ENDIF
 
-                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                       dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-                IF ( def < 0.0 )  def = 0.0
+                      dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                IF ( def < 0.0_wp )  def = 0.0_wp
 
                 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -1623,21 +1623,21 @@
              k = nzb_diff_s_outer(j,i)-1
 
-             dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-             dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-             dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-             dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-             dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-             dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-             dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-             dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-             dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+             dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+             dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+             dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+             dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+             dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+             dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+             dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+             dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+             dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
 
              def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
@@ -1645,5 +1645,5 @@
                    dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
 
-             IF ( def < 0.0 )  def = 0.0
+             IF ( def < 0.0_wp )  def = 0.0_wp
 
              tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -1657,59 +1657,59 @@
              k = nzb_diff_s_inner(j,i)-1
 
-             dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-             dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+             dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+             dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+             dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                 u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
+
+             dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+             dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+             dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                 v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
+
+             dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+             dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+             dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+             def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
+                   + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
+                   + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+             IF ( def < 0.0_wp )  def = 0.0_wp
+
+             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
+
+          ENDIF
+
+       ELSEIF ( use_surface_fluxes )  THEN
+
+          k = nzb_diff_s_outer(j,i)-1
+
+          dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+          dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
                               u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-             dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                              u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
-
-             dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+          dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+          dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
                               v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-             dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-             dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                              v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
-
-             dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+          dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+          dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+          dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
                               w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-             dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+          dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
                               w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-             dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
-                   + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
-                   + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-             IF ( def < 0.0 )  def = 0.0
-
-             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
-
-          ENDIF
-
-       ELSEIF ( use_surface_fluxes )  THEN
-
-          k = nzb_diff_s_outer(j,i)-1
-
-          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
-          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
-                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
-          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
-                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
-
-          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
-                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
-          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
-          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
-                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
-
-          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
-                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
-          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
-                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
-          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
-
-          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+          dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+          def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
                 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
-                dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
-
-          IF ( def < 0.0 )  def = 0.0
+                dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+          IF ( def < 0.0_wp )  def = 0.0_wp
 
           tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
@@ -1794,6 +1794,6 @@
 
                 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
-                   k1 = 1.0 + 0.61 * q(k,j,i)
-                   k2 = 0.61 * pt(k,j,i)
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
                    tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *   &
                                          ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
@@ -1801,16 +1801,16 @@
                                          ) * dd2zu(k)
                 ELSE IF ( cloud_physics )  THEN
-                   IF ( ql(k,j,i) == 0.0 )  THEN
-                      k1 = 1.0 + 0.61 * q(k,j,i)
-                      k2 = 0.61 * pt(k,j,i)
+                   IF ( ql(k,j,i) == 0.0_wp )  THEN
+                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                      k2 = 0.61_wp * pt(k,j,i)
                    ELSE
                       theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
                       temp  = theta * t_d_pt(k)
-                      k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                                 ( q(k,j,i) - ql(k,j,i) ) *          &
-                           ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                           ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
+                                    ( q(k,j,i) - ql(k,j,i) ) *          &
+                           ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
+                           ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
                            ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                      k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
                    ENDIF
                    tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *   &
@@ -1819,6 +1819,6 @@
                                          ) * dd2zu(k)
                 ELSE IF ( cloud_droplets )  THEN
-                   k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-                   k2 = 0.61 * pt(k,j,i)
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
                    tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *  &
                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +    &
@@ -1833,23 +1833,23 @@
 
                 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
-                   k1 = 1.0 + 0.61 * q(k,j,i)
-                   k2 = 0.61 * pt(k,j,i)
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
                 ELSE IF ( cloud_physics )  THEN
-                   IF ( ql(k,j,i) == 0.0 )  THEN
-                      k1 = 1.0 + 0.61 * q(k,j,i)
-                      k2 = 0.61 * pt(k,j,i)
+                   IF ( ql(k,j,i) == 0.0_wp )  THEN
+                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                      k2 = 0.61_wp * pt(k,j,i)
                    ELSE
                       theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
                       temp  = theta * t_d_pt(k)
-                      k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                                 ( q(k,j,i) - ql(k,j,i) ) *          &
-                           ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                           ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
+                                    ( q(k,j,i) - ql(k,j,i) ) *          &
+                           ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
+                           ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
                            ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                      k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
                    ENDIF
                 ELSE IF ( cloud_droplets )  THEN
-                   k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-                   k2 = 0.61 * pt(k,j,i)
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
                 ENDIF
 
@@ -1862,23 +1862,23 @@
 
                 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
-                   k1 = 1.0 + 0.61 * q(k,j,i)
-                   k2 = 0.61 * pt(k,j,i)
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
                 ELSE IF ( cloud_physics )  THEN
-                   IF ( ql(k,j,i) == 0.0 )  THEN
-                      k1 = 1.0 + 0.61 * q(k,j,i)
-                      k2 = 0.61 * pt(k,j,i)
+                   IF ( ql(k,j,i) == 0.0_wp )  THEN
+                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
+                      k2 = 0.61_wp * pt(k,j,i)
                    ELSE
                       theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
                       temp  = theta * t_d_pt(k)
-                      k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                                 ( q(k,j,i) - ql(k,j,i) ) *          &
-                           ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                           ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
+                                    ( q(k,j,i) - ql(k,j,i) ) *          &
+                           ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
+                           ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
                            ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                      k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
                    ENDIF
                 ELSE IF ( cloud_droplets )  THEN
-                   k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
-                   k2 = 0.61 * pt(k,j,i)
+                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
+                   k2 = 0.61_wp * pt(k,j,i)
                 ENDIF
 
@@ -1917,6 +1917,6 @@
           IF ( first_call )  THEN
              ALLOCATE( u_0(nysg:nyng,nxlg:nxrg), v_0(nysg:nyng,nxlg:nxrg) )
-             u_0 = 0.0   ! just to avoid access of uninitialized memory
-             v_0 = 0.0   ! within exchange_horiz_2d
+             u_0 = 0.0_wp   ! just to avoid access of uninitialized memory
+             v_0 = 0.0_wp   ! within exchange_horiz_2d
              first_call = .FALSE.
           ENDIF
@@ -1942,10 +1942,10 @@
 
                 u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
-                                 ( 0.5 * ( km(ku,j,i) + km(ku,j,i-1) ) +       &
-                                   1.0E-20 )
+                                 ( 0.5_wp * ( km(ku,j,i) + km(ku,j,i-1) ) +    &
+                                   1.0E-20_wp )
 !                                  ( us(j,i) * kappa * zu(1) )
                 v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
-                                 ( 0.5 * ( km(kv,j,i) + km(kv,j-1,i) ) +       &
-                                   1.0E-20 )
+                                 ( 0.5_wp * ( km(kv,j,i) + km(kv,j-1,i) ) +    &
+                                   1.0E-20_wp )
 !                                  ( us(j,i) * kappa * zu(1) )
 
Index: /palm/trunk/SOURCE/random_function.f90
===================================================================
--- /palm/trunk/SOURCE/random_function.f90	(revision 1341)
+++ /palm/trunk/SOURCE/random_function.f90	(revision 1342)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -104,6 +104,6 @@
        REAL(wp)     ::  rnmx             !:
 
-       PARAMETER ( ia=16807, im=2147483647, am=1.0/im, iq=127773, ir=2836, &
-                   ntab=32, ndiv=1+(im-1)/ntab, eps=1.2e-7, rnmx=1.0-eps )
+       PARAMETER ( ia=16807, im=2147483647, am=1.0_wp/im, iq=127773, ir=2836, &
+                   ntab=32, ndiv=1+(im-1)/ntab, eps=1.2e-7_wp, rnmx=1.0_wp-eps )
 
        IF ( idum .le. 0  .or.  random_iy .eq. 0 )  THEN
Index: /palm/trunk/SOURCE/random_gauss.f90
===================================================================
--- /palm/trunk/SOURCE/random_gauss.f90	(revision 1341)
+++ /palm/trunk/SOURCE/random_gauss.f90	(revision 1342)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! REAL constants defined as wp-kind
 ! 
 ! Former revisions:
@@ -78,20 +78,20 @@
 
        IF ( iset == 0 )  THEN
-          rsq = 0.0
-          DO  WHILE ( rsq >= 1.0  .OR.  rsq == 0.0 )
-             v1  = 2.0 * random_function( idum ) - 1.0
-             v2  = 2.0 * random_function( idum ) - 1.0
+          rsq = 0.0_wp
+          DO  WHILE ( rsq >= 1.0_wp  .OR.  rsq == 0.0_wp )
+             v1  = 2.0_wp * random_function( idum ) - 1.0_wp
+             v2  = 2.0_wp * random_function( idum ) - 1.0_wp
              rsq = v1**2 + v2**2
           ENDDO
-          fac          = SQRT( -2.0 * LOG( rsq ) / rsq )
+          fac          = SQRT( -2.0_wp * LOG( rsq ) / rsq )
           gset         = v1 * fac
-          random_gauss = v2 * fac + 1.0
+          random_gauss = v2 * fac + 1.0_wp
           iset         = 1
        ELSE
-          random_gauss = gset + 1.0
+          random_gauss = gset + 1.0_wp
           iset         = 0
        ENDIF
 
-       IF ( ABS( random_gauss - 1.0 ) <= upper_limit )  EXIT
+       IF ( ABS( random_gauss - 1.0_wp ) <= upper_limit )  EXIT
 
     ENDDO
Index: /palm/trunk/SOURCE/temperton_fft.f90
===================================================================
--- /palm/trunk/SOURCE/temperton_fft.f90	(revision 1341)
+++ /palm/trunk/SOURCE/temperton_fft.f90	(revision 1342)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -167,5 +167,5 @@
 !OCL NOVREC
        DO j = 1, nvex
-          a(i+inc) = 0.5*a(i)
+          a(i+inc) = 0.5_wp*a(i)
           i = i + jump
        END DO
@@ -173,5 +173,5 @@
        i = istart + n*inc
        DO j = 1, nvex
-          a(i) = 0.5*a(i)
+          a(i) = 0.5_wp*a(i)
           i = i + jump
        END DO
@@ -223,6 +223,6 @@
 !OCL NOVREC
        DO j = 1, nvex
-          a(ix) = 0.0
-          a(ix+inc) = 0.0
+          a(ix) = 0.0_wp
+          a(ix+inc) = 0.0_wp
           ix = ix + jump
        END DO
@@ -285,5 +285,5 @@
        DO j = 1, nvex
           a(ix) = a(ix+inc)
-          a(ix+inc) = 0.0
+          a(ix+inc) = 0.0_wp
           ix = ix + jump
        END DO
@@ -291,5 +291,5 @@
        iz = istart + (n+1)*inc
        DO j = 1, nvex
-          a(iz) = 0.0
+          a(iz) = 0.0_wp
           iz = iz + jump
        END DO
@@ -552,5 +552,5 @@
     GO TO 170
 30  CONTINUE
-    z = 1.0_wp/REAL(n)
+    z = 1.0_wp/REAL(n,KIND=wp)
     DO l = 1, la
        i = ibase
@@ -590,5 +590,5 @@
        DO ijk = 1, lot
           c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
-          c(jb+j) = a(ia+i) - 0.5*(a(ib+i)+a(ic+i))
+          c(jb+j) = a(ia+i) - 0.5_wp*(a(ib+i)+a(ic+i))
           d(jb+j) = sin60*(a(ic+i)-a(ib+i))
           i = i + inc3
@@ -622,6 +622,6 @@
              a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c2*a(ic+i)+s2*b(ic+i))
              b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c2*b(ic+i)-s2*a(ic+i))
-             a2 = a(ia+i) - 0.5*a1
-             b2 = b(ia+i) - 0.5*b1
+             a2 = a(ia+i) - 0.5_wp*a1
+             b2 = b(ia+i) - 0.5_wp*b1
              a3 = sin60*((c1*a(ib+i)+s1*b(ib+i))-(c2*a(ic+i)+s2*b(ic+i)))
              b3 = sin60*((c1*b(ib+i)-s1*a(ib+i))-(c2*b(ic+i)-s2*a(ic+i)))
@@ -653,5 +653,5 @@
 !OCL NOVREC
        DO ijk = 1, lot
-          c(ja+j) = a(ia+i) + 0.5*(a(ib+i)-a(ic+i))
+          c(ja+j) = a(ia+i) + 0.5_wp*(a(ib+i)-a(ic+i))
           d(ja+j) = -sin60*(a(ib+i)+a(ic+i))
           c(jb+j) = a(ia+i) - (a(ib+i)-a(ic+i))
@@ -665,5 +665,5 @@
     GO TO 170
 60  CONTINUE
-    z = 1.0_wp/REAL(n)
+    z = 1.0_wp/REAL(n,KIND=wp)
     zsin60 = z*sin60
     DO l = 1, la
@@ -675,5 +675,5 @@
        DO ijk = 1, lot
           c(ja+j) = z*(a(ia+i)+(a(ib+i)+a(ic+i)))
-          c(jb+j) = z*(a(ia+i)-0.5*(a(ib+i)+a(ic+i)))
+          c(jb+j) = z*(a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))
           d(jb+j) = zsin60*(a(ic+i)-a(ib+i))
           i = i + inc3
@@ -794,5 +794,5 @@
     GO TO 170
 90  CONTINUE
-    z = 1.0_wp/REAL(n)
+    z = 1.0_wp/REAL(n,KIND=wp)
     DO l = 1, la
        i = ibase
@@ -841,5 +841,5 @@
           a2 = a(ic+i) + a(id+i)
           a4 = a(ic+i) - a(id+i)
-          a5 = a(ia+i) - 0.25*(a1+a2)
+          a5 = a(ia+i) - 0.25_wp*(a1+a2)
           a6 = qrt5*(a1-a2)
           c(ja+j) = a(ia+i) + (a1+a2)
@@ -892,7 +892,7 @@
              b2 = (c2*b(ic+i)-s2*a(ic+i)) + (c3*b(id+i)-s3*a(id+i))
              b4 = (c2*b(ic+i)-s2*a(ic+i)) - (c3*b(id+i)-s3*a(id+i))
-             a5 = a(ia+i) - 0.25*(a1+a2)
+             a5 = a(ia+i) - 0.25_wp*(a1+a2)
              a6 = qrt5*(a1-a2)
-             b5 = b(ia+i) - 0.25*(b1+b2)
+             b5 = b(ia+i) - 0.25_wp*(b1+b2)
              b6 = qrt5*(b1-b2)
              a10 = a5 + a6
@@ -941,5 +941,5 @@
           a2 = a(ic+i) + a(id+i)
           a4 = a(ic+i) - a(id+i)
-          a5 = a(ia+i) + 0.25*(a3-a4)
+          a5 = a(ia+i) + 0.25_wp*(a3-a4)
           a6 = qrt5*(a3+a4)
           c(ja+j) = a5 + a6
@@ -957,5 +957,5 @@
     GO TO 170
 120 CONTINUE
-    z = 1.0_wp/REAL(n)
+    z = 1.0_wp/REAL(n,KIND=wp)
     zqrt5 = z*qrt5
     zsin36 = z*sin36
@@ -972,5 +972,5 @@
           a2 = a(ic+i) + a(id+i)
           a4 = a(ic+i) - a(id+i)
-          a5 = z*(a(ia+i)-0.25*(a1+a2))
+          a5 = z*(a(ia+i)-0.25_wp*(a1+a2))
           a6 = zqrt5*(a1-a2)
           c(ja+j) = z*(a(ia+i)+(a1+a2))
@@ -1014,8 +1014,8 @@
           a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
           c(ja+j) = (a(ia+i)+a(id+i)) + a11
-          c(jc+j) = (a(ia+i)+a(id+i)-0.5*a11)
+          c(jc+j) = (a(ia+i)+a(id+i)-0.5_wp*a11)
           d(jc+j) = sin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
           a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
-          c(jb+j) = (a(ia+i)-a(id+i)) - 0.5*a11
+          c(jb+j) = (a(ia+i)-a(id+i)) - 0.5_wp*a11
           d(jb+j) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
           c(jd+j) = (a(ia+i)-a(id+i)) + a11
@@ -1071,8 +1071,8 @@
              b5 = c5*b(if+i) - s5*a(if+i)
              a11 = (a2+a5) + (a1+a4)
-             a20 = (a(ia+i)+a3) - 0.5*a11
+             a20 = (a(ia+i)+a3) - 0.5_wp*a11
              a21 = sin60*((a2+a5)-(a1+a4))
              b11 = (b2+b5) + (b1+b4)
-             b20 = (b(ia+i)+b3) - 0.5*b11
+             b20 = (b(ia+i)+b3) - 0.5_wp*b11
              b21 = sin60*((b2+b5)-(b1+b4))
              c(ja+j) = (a(ia+i)+a3) + a11
@@ -1083,8 +1083,8 @@
              d(je+j) = a21 - b20
              a11 = (a2-a5) + (a4-a1)
-             a20 = (a(ia+i)-a3) - 0.5*a11
+             a20 = (a(ia+i)-a3) - 0.5_wp*a11
              a21 = sin60*((a4-a1)-(a2-a5))
              b11 = (b5-b2) - (b4-b1)
-             b20 = (b3-b(ia+i)) - 0.5*b11
+             b20 = (b3-b(ia+i)) - 0.5_wp*b11
              b21 = sin60*((b5-b2)+(b4-b1))
              c(jb+j) = a20 - b21
@@ -1118,10 +1118,10 @@
 !OCL NOVREC
        DO ijk = 1, lot
-          c(ja+j) = (a(ia+i)+0.5*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i))
-          d(ja+j) = -(a(id+i)+0.5*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i))
+          c(ja+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i))
+          d(ja+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i))
           c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i))
           d(jb+j) = a(id+i) - (a(ib+i)+a(if+i))
-          c(jc+j) = (a(ia+i)+0.5*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i))
-          d(jc+j) = -(a(id+i)+0.5*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i))
+          c(jc+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i))
+          d(jc+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i))
           i = i + inc3
           j = j + inc4
@@ -1133,5 +1133,5 @@
     GO TO 170
 150 CONTINUE
-    z = 1.0_wp/REAL(n)
+    z = 1.0_wp/REAL(n,KIND=wp)
     zsin60 = z*sin60
     DO l = 1, la
@@ -1144,8 +1144,8 @@
           a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
           c(ja+j) = z*((a(ia+i)+a(id+i))+a11)
-          c(jc+j) = z*((a(ia+i)+a(id+i))-0.5*a11)
+          c(jc+j) = z*((a(ia+i)+a(id+i))-0.5_wp*a11)
           d(jc+j) = zsin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
           a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
-          c(jb+j) = z*((a(ia+i)-a(id+i))-0.5*a11)
+          c(jb+j) = z*((a(ia+i)-a(id+i))-0.5_wp*a11)
           d(jb+j) = zsin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
           c(jd+j) = z*((a(ia+i)-a(id+i))+a11)
@@ -1176,5 +1176,5 @@
     jd = jc + 2*m*inc2
     je = jd + 2*m*inc2
-    z = 1.0_wp/REAL(n)
+    z = 1.0_wp/REAL(n,KIND=wp)
     zsin45 = z*SQRT(0.5_wp)
 
@@ -1419,6 +1419,6 @@
 !OCL NOVREC
        DO ijk = 1, lot
-          c(ja+j) = 2.0*(a(ia+i)+a(ib+i))
-          c(jb+j) = 2.0*(a(ia+i)-a(ib+i))
+          c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
+          c(jb+j) = 2.0_wp*(a(ia+i)-a(ib+i))
           i = i + inc3
           j = j + inc4
@@ -1449,6 +1449,6 @@
        DO ijk = 1, lot
           c(ja+j) = a(ia+i) + a(ib+i)
-          c(jb+j) = (a(ia+i)-0.5*a(ib+i)) - (sin60*(b(ib+i)))
-          c(jc+j) = (a(ia+i)-0.5*a(ib+i)) + (sin60*(b(ib+i)))
+          c(jb+j) = (a(ia+i)-0.5_wp*a(ib+i)) - (sin60*(b(ib+i)))
+          c(jc+j) = (a(ia+i)-0.5_wp*a(ib+i)) + (sin60*(b(ib+i)))
           i = i + inc3
           j = j + inc4
@@ -1481,15 +1481,19 @@
              c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
              d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i))
-             c(jb+j) = c1*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
-                  &            b(ic+i)))) - s1*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
+             c(jb+j) = c1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
+                  &            b(ic+i))))                                      &
+                  &    - s1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
                   &            a(ic+i))))
-             d(jb+j) = s1*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
-                  &            b(ic+i)))) + c1*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
+             d(jb+j) = s1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
+                  &            b(ic+i))))                                      &
+                  &    + c1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
                   &            a(ic+i))))
-             c(jc+j) = c2*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
-                  &            b(ic+i)))) - s2*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
+             c(jc+j) = c2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
+                  &            b(ic+i))))                                      &
+                  &    - s2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
                   &            a(ic+i))))
-             d(jc+j) = s2*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
-                  &            b(ic+i)))) + c2*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
+             d(jc+j) = s2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
+                  &            b(ic+i))))                                      &
+                  &    + c2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
                   &            a(ic+i))))
              i = i + inc3
@@ -1515,6 +1519,6 @@
        DO ijk = 1, lot
           c(ja+j) = a(ia+i) + a(ib+i)
-          c(jb+j) = (0.5*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
-          c(jc+j) = -(0.5*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
+          c(jb+j) = (0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
+          c(jc+j) = -(0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
           i = i + inc3
           j = j + inc4
@@ -1526,15 +1530,15 @@
     GO TO 170
 60  CONTINUE
-    ssin60 = 2.0*sin60
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = 2.0*(a(ia+i)+a(ib+i))
-          c(jb+j) = (2.0*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i))
-          c(jc+j) = (2.0*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i))
+    ssin60 = 2.0_wp*sin60
+    DO l = 1, la
+       i = ibase
+       j = jbase
+!DIR$ IVDEP
+!CDIR NODEP
+!OCL NOVREC
+       DO ijk = 1, lot
+          c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
+          c(jb+j) = (2.0_wp*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i))
+          c(jc+j) = (2.0_wp*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i))
           i = i + inc3
           j = j + inc4
@@ -1659,8 +1663,8 @@
 !OCL NOVREC
        DO ijk = 1, lot
-          c(ja+j) = 2.0*((a(ia+i)+a(ic+i))+a(ib+i))
-          c(jb+j) = 2.0*((a(ia+i)-a(ic+i))-b(ib+i))
-          c(jc+j) = 2.0*((a(ia+i)+a(ic+i))-a(ib+i))
-          c(jd+j) = 2.0*((a(ia+i)-a(ic+i))+b(ib+i))
+          c(ja+j) = 2.0_wp*((a(ia+i)+a(ic+i))+a(ib+i))
+          c(jb+j) = 2.0_wp*((a(ia+i)-a(ic+i))-b(ib+i))
+          c(jc+j) = 2.0_wp*((a(ia+i)+a(ic+i))-a(ib+i))
+          c(jd+j) = 2.0_wp*((a(ia+i)-a(ic+i))+b(ib+i))
           i = i + inc3
           j = j + inc4
@@ -1695,11 +1699,11 @@
        DO ijk = 1, lot
           c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
-          c(jb+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) - &
+          c(jb+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) - &
                &          (sin72*b(ib+i)+sin36*b(ic+i))
-          c(jc+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) - &
+          c(jc+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) - &
                &          (sin36*b(ib+i)-sin72*b(ic+i))
-          c(jd+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) + &
+          c(jd+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) + &
                &          (sin36*b(ib+i)-sin72*b(ic+i))
-          c(je+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) + &
+          c(je+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) + &
                &          (sin72*b(ib+i)+sin36*b(ic+i))
           i = i + inc3
@@ -1740,11 +1744,11 @@
           DO ijk = 1, lot
 
-             a10(ijk) = (a(ia+i)-0.25*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + &
+             a10(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + &
                   &            qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
-             a20(ijk) = (a(ia+i)-0.25*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - &
+             a20(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - &
                   &            qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
-             b10(ijk) = (b(ia+i)-0.25*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + &
+             b10(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + &
                   &            qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
-             b20(ijk) = (b(ia+i)-0.25*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - &
+             b20(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - &
                   &            qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
              a11(ijk) = sin72*(b(ib+i)+b(ie+i)) + sin36*(b(ic+i)+b(id+i))
@@ -1788,11 +1792,11 @@
        DO ijk = 1, lot
           c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i)
-          c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
+          c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
                &          (sin36*b(ia+i)+sin72*b(ib+i))
-          c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
+          c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
                &          (sin36*b(ia+i)+sin72*b(ib+i))
-          c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
+          c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
                &          (sin72*b(ia+i)-sin36*b(ib+i))
-          c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
+          c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
                &          (sin72*b(ia+i)-sin36*b(ib+i))
           i = i + inc3
@@ -1815,12 +1819,12 @@
 !OCL NOVREC
        DO ijk = 1, lot
-          c(ja+j) = 2.0*(a(ia+i)+(a(ib+i)+a(ic+i)))
-          c(jb+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
+          c(ja+j) = 2.0_wp*(a(ia+i)+(a(ib+i)+a(ic+i)))
+          c(jb+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
                &          i))) - (ssin72*b(ib+i)+ssin36*b(ic+i))
-          c(jc+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
+          c(jc+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
                &          i))) - (ssin36*b(ib+i)-ssin72*b(ic+i))
-          c(jd+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
+          c(jd+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
                &          i))) + (ssin36*b(ib+i)-ssin72*b(ic+i))
-          c(je+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
+          c(je+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
                &          i))) + (ssin72*b(ib+i)+ssin36*b(ic+i))
           i = i + inc3
@@ -1859,11 +1863,11 @@
           c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i))
           c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i))
-          c(jb+j) = ((a(ia+i)-a(id+i))+0.5*(a(ib+i)-a(ic+i))) - (sin60*(b(ib+ &
+          c(jb+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) - (sin60*(b(ib+ &
                &          i)+b(ic+i)))
-          c(jf+j) = ((a(ia+i)-a(id+i))+0.5*(a(ib+i)-a(ic+i))) + (sin60*(b(ib+ &
+          c(jf+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) + (sin60*(b(ib+ &
                &          i)+b(ic+i)))
-          c(jc+j) = ((a(ia+i)+a(id+i))-0.5*(a(ib+i)+a(ic+i))) - (sin60*(b(ib+ &
+          c(jc+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) - (sin60*(b(ib+ &
                &          i)-b(ic+i)))
-          c(je+j) = ((a(ia+i)+a(id+i))-0.5*(a(ib+i)+a(ic+i))) + (sin60*(b(ib+ &
+          c(je+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) + (sin60*(b(ib+ &
                &          i)-b(ic+i)))
           i = i + inc3
@@ -1909,8 +1913,8 @@
 
              a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i))
-             a20(ijk) = (a(ia+i)+a(id+i)) - 0.5*a11(ijk)
+             a20(ijk) = (a(ia+i)+a(id+i)) - 0.5_wp*a11(ijk)
              a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i)))
              b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i))
-             b20(ijk) = (b(ia+i)-b(id+i)) - 0.5*b11(ijk)
+             b20(ijk) = (b(ia+i)-b(id+i)) - 0.5_wp*b11(ijk)
              b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i)))
 
@@ -1924,7 +1928,7 @@
              a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i))
              b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i))
-             a20(ijk) = (a(ia+i)-a(id+i)) - 0.5*a11(ijk)
+             a20(ijk) = (a(ia+i)-a(id+i)) - 0.5_wp*a11(ijk)
              a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
-             b20(ijk) = (b(ia+i)+b(id+i)) + 0.5*b11(ijk)
+             b20(ijk) = (b(ia+i)+b(id+i)) + 0.5_wp*b11(ijk)
              b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i)))
 
@@ -1964,8 +1968,8 @@
           c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i))
           c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i))
-          c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5*(b(ia+i)+b(ic+i))+b(ib+i))
-          c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5*(b(ia+i)+b(ic+i))+b(ib+i))
-          c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5*(a(ia+i)+a(ic+i))-a(ib+i))
-          c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5*(a(ia+i)+a(ic+i))-a(ib+i))
+          c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
+          c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
+          c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
+          c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
           i = i + inc3
           j = j + inc4
@@ -1977,21 +1981,21 @@
     GO TO 170
 150 CONTINUE
-    ssin60 = 2.0*sin60
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = (2.0*(a(ia+i)+a(id+i))) + (2.0*(a(ib+i)+a(ic+i)))
-          c(jd+j) = (2.0*(a(ia+i)-a(id+i))) - (2.0*(a(ib+i)-a(ic+i)))
-          c(jb+j) = (2.0*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+ &
+    ssin60 = 2.0_wp*sin60
+    DO l = 1, la
+       i = ibase
+       j = jbase
+!DIR$ IVDEP
+!CDIR NODEP
+!OCL NOVREC
+       DO ijk = 1, lot
+          c(ja+j) = (2.0_wp*(a(ia+i)+a(id+i))) + (2.0_wp*(a(ib+i)+a(ic+i)))
+          c(jd+j) = (2.0_wp*(a(ia+i)-a(id+i))) - (2.0_wp*(a(ib+i)-a(ic+i)))
+          c(jb+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+ &
                &          i)+b(ic+i)))
-          c(jf+j) = (2.0*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+ &
+          c(jf+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+ &
                &          i)+b(ic+i)))
-          c(jc+j) = (2.0*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+ &
+          c(jc+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+ &
                &          i)-b(ic+i)))
-          c(je+j) = (2.0*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+ &
+          c(je+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+ &
                &          i)-b(ic+i)))
           i = i + inc3
@@ -2030,16 +2034,16 @@
 !OCL NOVREC
        DO ijk = 1, lot
-          c(ja+j) = 2.0*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))
-          c(je+j) = 2.0*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))
-          c(jc+j) = 2.0*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))
-          c(jg+j) = 2.0*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))
-          c(jb+j) = 2.0*((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
-               &          i))-(b(ib+i)+b(id+i)))
-          c(jf+j) = 2.0*((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
-               &          i))-(b(ib+i)+b(id+i)))
-          c(jd+j) = 2.0*((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
-               &          i))+(b(ib+i)+b(id+i)))
-          c(jh+j) = 2.0*((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
-               &          i))+(b(ib+i)+b(id+i)))
+          c(ja+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))
+          c(je+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))
+          c(jc+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))
+          c(jg+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))
+          c(jb+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
+               &              i))-(b(ib+i)+b(id+i)))
+          c(jf+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
+               &              i))-(b(ib+i)+b(id+i)))
+          c(jd+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
+               &              i))+(b(ib+i)+b(id+i)))
+          c(jh+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
+               &              i))+(b(ib+i)+b(id+i)))
           i = i + inc3
           j = j + inc4
@@ -2114,9 +2118,9 @@
     ixxx = 1
 
-    del = 4.0_wp*ASIN(1.0_wp)/REAL(n)
+    del = 4.0_wp*ASIN(1.0_wp)/REAL(n,KIND=wp)
     nil = 0
     nhl = (n/2) - 1
     DO k = nil, nhl
-       angle = REAL(k)*del
+       angle = REAL(k,KIND=wp)*del
        trigs(2*k+1) = COS(angle)
        trigs(2*k+2) = SIN(angle)
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 1341)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 1342)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! REAL constants defined as wp-kind
 ! 
 ! Former revisions:
@@ -205,5 +205,5 @@
 !-- At the beginning of a simulation determine the time step as well as
 !-- determine and print out the run control parameters
-    IF ( simulated_time == 0.0 )  CALL timestep
+    IF ( simulated_time == 0.0_wp )  CALL timestep
 
     CALL run_control
@@ -221,5 +221,5 @@
           time_coupling = time_coupling - dt_coupling
        ENDDO
-       IF (time_coupling == 0.0 .AND. time_since_reference_point < dt_coupling)&
+       IF (time_coupling == 0.0_wp .AND. time_since_reference_point < dt_coupling)&
        THEN
           time_coupling = time_since_reference_point
@@ -243,5 +243,5 @@
 !
 !--    Determine size of next time step
-       IF ( simulated_time /= 0.0 )  THEN
+       IF ( simulated_time /= 0.0_wp )  THEN
           CALL timestep
        ENDIF
@@ -717,5 +717,5 @@
 !
 !--    If required, call flow_statistics for averaging in time
-       IF ( averaging_interval_pr /= 0.0  .AND.  &
+       IF ( averaging_interval_pr /= 0.0_wp  .AND.  &
             ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.  &
             simulated_time >= skip_time_dopr )  THEN
@@ -731,5 +731,5 @@
 !
 !--    Sum-up 3d-arrays for later output of time-averaged 2d/3d/masked data
-       IF ( averaging_interval /= 0.0  .AND.                                &
+       IF ( averaging_interval /= 0.0_wp  .AND.                                &
             ( dt_data_output_av - time_do_av ) <= averaging_interval  .AND. &
             simulated_time >= skip_time_data_output_av )                    &
@@ -745,5 +745,5 @@
 !
 !--    Calculate spectra for time averaging
-       IF ( averaging_interval_sp /= 0.0  .AND.  &
+       IF ( averaging_interval_sp /= 0.0_wp  .AND.  &
             ( dt_dosp - time_dosp ) <= averaging_interval_sp  .AND.  &
             simulated_time >= skip_time_dosp )  THEN
@@ -782,5 +782,5 @@
           IF ( dopr_n /= 0 )  CALL data_output_profiles
           time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
-          time_dopr_av = 0.0    ! due to averaging (see above)
+          time_dopr_av = 0.0_wp    ! due to averaging (see above)
        ENDIF
 
Index: /palm/trunk/SOURCE/time_to_string.f90
===================================================================
--- /palm/trunk/SOURCE/time_to_string.f90	(revision 1341)
+++ /palm/trunk/SOURCE/time_to_string.f90	(revision 1342)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! REAL constants defined as wp-kind
 ! 
 ! Former revisions:
@@ -61,7 +61,7 @@
 !
 !-- Calculate the number of hours, minutes, and seconds
-    hours     = INT( time / 3600.0 )
-    rest_time = time - hours * 3600
-    minutes   = INT( rest_time / 60.0 )
+    hours     = INT( time / 3600.0_wp )
+    rest_time = time - hours * 3600_wp
+    minutes   = INT( rest_time / 60.0_wp )
     seconds   = rest_time - minutes * 60
 
Index: /palm/trunk/SOURCE/timestep.f90
===================================================================
--- /palm/trunk/SOURCE/timestep.f90	(revision 1341)
+++ /palm/trunk/SOURCE/timestep.f90	(revision 1342)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -153,6 +153,6 @@
 !--       Horizontal averages already existent, just need to average them 
 !--       vertically.
-          u_gtrans = 0.0
-          v_gtrans = 0.0
+          u_gtrans = 0.0_wp
+          v_gtrans = 0.0_wp
           DO  k = nzb+1, nzt
              u_gtrans = u_gtrans + hom(k,1,1,0)
@@ -164,6 +164,6 @@
 !
 !--       Averaging over the entire model domain.
-          u_gtrans_l = 0.0
-          v_gtrans_l = 0.0
+          u_gtrans_l = 0.0_wp
+          v_gtrans_l = 0.0_wp
           !$acc parallel present( u, v )
           DO  i = nxl, nxr
@@ -196,10 +196,10 @@
 #if defined( __openacc )
     IF ( dt_fixed )  THEN  ! otherwise do it further below for better cache usage
-       u_max_l = -999999.9
-       u_min_l =  999999.9
-       v_max_l = -999999.9
-       v_min_l =  999999.9
-       w_max_l = -999999.9
-       w_min_l =  999999.9
+       u_max_l = -999999.9_wp
+       u_min_l =  999999.9_wp
+       v_max_l = -999999.9_wp
+       v_min_l =  999999.9_wp
+       w_max_l = -999999.9_wp
+       w_min_l =  999999.9_wp
        !$acc parallel present( u, v, w )
        DO  i = nxl, nxr
@@ -266,20 +266,20 @@
 !--    Calculate the maximum time step according to the CFL-criterion,
 !--    individually for each velocity component
-       dt_u_l  =  999999.9
-       dt_v_l  =  999999.9
-       dt_w_l  =  999999.9
-       u_max_l = -999999.9
-       u_min_l =  999999.9
-       v_max_l = -999999.9
-       v_min_l =  999999.9
-       w_max_l = -999999.9
-       w_min_l =  999999.9
+       dt_u_l  =  999999.9_wp
+       dt_v_l  =  999999.9_wp
+       dt_w_l  =  999999.9_wp
+       u_max_l = -999999.9_wp
+       u_min_l =  999999.9_wp
+       v_max_l = -999999.9_wp
+       v_min_l =  999999.9_wp
+       w_max_l = -999999.9_wp
+       w_min_l =  999999.9_wp
        !$acc parallel loop collapse(3) present( u, v, w )
        DO  i = nxl, nxr
           DO  j = nys, nyn
              DO  k = nzb+1, nzt
-                dt_u_l  = MIN( dt_u_l, ( dx     / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10 ) ) )
-                dt_v_l  = MIN( dt_v_l, ( dy     / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10 ) ) )
-                dt_w_l  = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) )            + 1.0E-10 ) ) )
+                dt_u_l  = MIN( dt_u_l, ( dx     / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10_wp ) ) )
+                dt_v_l  = MIN( dt_v_l, ( dy     / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10_wp ) ) )
+                dt_w_l  = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) )            + 1.0E-10_wp ) ) )
                 u_max_l = MAX( u_max_l, u(k,j,i) )
                 u_min_l = MIN( u_min_l, u(k,j,i) )
@@ -346,13 +346,13 @@
 !--    Calculate the maximum time step according to the CFL-criterion,
 !--    individually for each velocity component
-       dt_u_l = 999999.9
-       dt_v_l = 999999.9
-       dt_w_l = 999999.9
+       dt_u_l = 999999.9_wp
+       dt_v_l = 999999.9_wp
+       dt_w_l = 999999.9_wp
        DO  i = nxl, nxr
           DO  j = nys, nyn
              DO  k = nzb+1, nzt
-                dt_u_l = MIN( dt_u_l, ( dx     / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10 ) ) )
-                dt_v_l = MIN( dt_v_l, ( dy     / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10 ) ) )
-                dt_w_l = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) )            + 1.0E-10 ) ) )
+                dt_u_l = MIN( dt_u_l, ( dx     / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10_wp ) ) )
+                dt_v_l = MIN( dt_v_l, ( dy     / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10_wp ) ) )
+                dt_w_l = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) )            + 1.0E-10_wp ) ) )
              ENDDO
           ENDDO
@@ -382,8 +382,8 @@
 !--          in the Prandtl-layer friction term only dz/2 is used.
 !--          Experience from the old model seems to justify this.
-       dt_diff_l = 999999.0
+       dt_diff_l = 999999.0_wp
 
        DO  k = nzb+1, nzt
-           dxyz2_min(k) = MIN( dx2, dy2, dzw(k)*dzw(k) ) * 0.125
+           dxyz2_min(k) = MIN( dx2, dy2, dzw(k)*dzw(k) ) * 0.125_wp
        ENDDO
 
@@ -395,5 +395,5 @@
              DO  k = nzb+1, nzt
                 dt_diff_l = MIN( dt_diff_l, dxyz2_min(k) / &
-                                       ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20 ) )
+                                       ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20_wp ) )
              ENDDO
           ENDDO
@@ -414,5 +414,5 @@
        IF ( plant_canopy ) THEN
 
-          dt_plant_canopy_l = 0.0
+          dt_plant_canopy_l = 0.0_wp
           DO  i = nxl, nxr
              DO  j = nys, nyn
@@ -424,10 +424,10 @@
                                                  v(k,j+1,i)      +  &
                                                  v(k,j+1,i-1) )     &
-                                               / 4.0 )**2        +  &
+                                               / 4.0_wp )**2     +  &
                                              ( ( w(k-1,j,i-1)    +  &
                                                  w(k-1,j,i)      +  &
                                                  w(k,j,i-1)      +  &
                                                  w(k,j,i) )         &
-                                                 / 4.0 )**2 ) 
+                                                 / 4.0_wp )**2 ) 
                    IF ( dt_plant_canopy_u > dt_plant_canopy_l ) THEN
                       dt_plant_canopy_l = dt_plant_canopy_u  
@@ -438,5 +438,5 @@
                                                  u(k,j,i)        +  &
                                                  u(k,j,i+1) )       &
-                                               / 4.0 )**2        +  &
+                                               / 4.0_wp )**2     +  &
                                                  v(k,j,i)**2     +  &
                                              ( ( w(k-1,j-1,i)    +  &
@@ -444,5 +444,5 @@
                                                  w(k,j-1,i)      +  &
                                                  w(k,j,i) )         &
-                                                 / 4.0 )**2 ) 
+                                                 / 4.0_wp )**2 ) 
                    IF ( dt_plant_canopy_v > dt_plant_canopy_l ) THEN
                       dt_plant_canopy_l = dt_plant_canopy_v
@@ -453,10 +453,10 @@
                                                  u(k+1,j,i)      +  &
                                                  u(k+1,j,i+1) )     &
-                                               / 4.0 )**2        +  &
+                                               / 4.0_wp )**2     +  &
                                              ( ( v(k,j,i)        +  &
                                                  v(k,j+1,i)      +  &
                                                  v(k+1,j,i)      +  &
                                                  v(k+1,j+1,i) )     &
-                                               / 4.0 )**2        +  &
+                                               / 4.0_wp )**2        +  &
                                                  w(k,j,i)**2 )     
                    IF ( dt_plant_canopy_w > dt_plant_canopy_l ) THEN
@@ -467,8 +467,8 @@
           ENDDO 
 
-          IF ( dt_plant_canopy_l > 0.0 ) THEN
+          IF ( dt_plant_canopy_l > 0.0_wp ) THEN
 !
 !--          Invert dt_plant_canopy_l and apply a security timestep factor 0.1
-             dt_plant_canopy_l = 0.1 / dt_plant_canopy_l
+             dt_plant_canopy_l = 0.1_wp / dt_plant_canopy_l
           ELSE
 !
@@ -515,5 +515,5 @@
 !
 !--    Set flag if the time step becomes too small.
-       IF ( dt_3d < ( 0.00001 * dt_max ) )  THEN
+       IF ( dt_3d < ( 0.00001_wp * dt_max ) )  THEN
           stop_dt = .TRUE.
 
@@ -553,9 +553,9 @@
 !
 !--    Ensure a smooth value (two significant digits) of the timestep.
-       div = 1000.0
+       div = 1000.0_wp
        DO  WHILE ( dt_3d < div )
-          div = div / 10.0
+          div = div / 10.0_wp
        ENDDO
-       dt_3d = NINT( dt_3d * 100.0 / div ) * div / 100.0
+       dt_3d = NINT( dt_3d * 100.0_wp / div ) * div / 100.0_wp
 
 !
Index: /palm/trunk/SOURCE/tridia_solver.f90
===================================================================
--- /palm/trunk/SOURCE/tridia_solver.f90	(revision 1341)
+++ /palm/trunk/SOURCE/tridia_solver.f90	(revision 1342)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! REAL constants defined as wp-kind
 !
 ! Former revisions:
@@ -110,5 +110,5 @@
           ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1)
           ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1)
-          ddzuw(k,3) = -1.0 * &
+          ddzuw(k,3) = -1.0_wp * &
                        ( ddzu_pres(k+2) * ddzw(k+1) + ddzu_pres(k+1) * ddzw(k+1) )
        ENDDO
@@ -172,25 +172,25 @@
                 IF ( j >= 0  .AND.  j <= nnyh )  THEN
                    IF ( i >= 0  .AND.  i <= nnxh )  THEN
-                      ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
-                                                  REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
-                                2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
-                                                  REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
+                      ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
+                                            REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
+                                2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
+                                            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, KIND=wp ) ) ) / ( dx * dx ) + &
-                                2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
-                                                  REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
+                      ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
+                                            REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
+                                2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
+                                            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, KIND=wp ) ) ) / ( dx * dx ) + &
-                                2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
-                                                  REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
+                      ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
+                                            REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
+                                2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / &
+                                            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, KIND=wp ) ) ) / ( dx * dx ) + &
-                                2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
-                                                  REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
+                      ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
+                                            REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
+                                2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / &
+                                            REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
                    ENDIF
                 ENDIF
@@ -283,5 +283,5 @@
 
                    IF ( k == nz-1 )  THEN
-                      ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20 )
+                      ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20_wp )
                    ELSE
                       ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
@@ -301,5 +301,5 @@
                 !$acc kernels loop present( ar )
                 DO  k = 1, nz
-                   ar(nxl_z,nys_z,k) = 0.0
+                   ar(nxl_z,nys_z,k) = 0.0_wp
                 ENDDO
                 !$acc end kernels loop
@@ -367,5 +367,5 @@
 
                    IF ( k == nz-1 )  THEN
-                      ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20 )
+                      ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20_wp )
                    ELSE
                       ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
@@ -385,5 +385,5 @@
                 !$acc kernels loop present( ar )
                 DO  k = 1, nz
-                   ar(nxl_z,nys_z,k) = 0.0
+                   ar(nxl_z,nys_z,k) = 0.0_wp
                 ENDDO
              ENDIF
@@ -560,13 +560,13 @@
           DO  i = 0, nx
              IF ( i >= 0 .AND. i <= nnxh ) THEN
-                l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
-                                         REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
-                       2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
-                                         REAL( ny+1, KIND=wp ) ) ) * ddy2
+                l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
+                                   REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
+                       2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
+                                   REAL( ny+1, KIND=wp ) ) ) * ddy2
              ELSE
-                l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
-                                         REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
-                       2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
-                                         REAL( ny+1, KIND=wp ) ) ) * ddy2
+                l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
+                                   REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
+                       2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
+                                   REAL( ny+1, KIND=wp ) ) ) * ddy2
              ENDIF
           ENDDO
@@ -574,6 +574,6 @@
           DO  k = 0, nz-1
              DO  i = 0, nx
-                a = -1.0 * ddzu_pres(k+2) * ddzw(k+1)
-                c = -1.0 * ddzu_pres(k+1) * ddzw(k+1)
+                a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1)
+                c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1)
                 tri_for_1d(1,i,k) = a + c - l(i)
              ENDDO
@@ -660,5 +660,5 @@
 !--       the model domain.
           DO  i = 0, nx
-             ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20 )
+             ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20_wp )
           ENDDO
           DO  k = nz-2, 0, -1
@@ -676,5 +676,5 @@
              IF ( j == 0 )  THEN
                 DO  k = 1, nz
-                   ar(0,k) = 0.0
+                   ar(0,k) = 0.0_wp
                 ENDDO
              ENDIF
