Index: palm/trunk/SOURCE/advec_s_bc.f90
===================================================================
--- palm/trunk/SOURCE/advec_s_bc.f90	(revision 1347)
+++ palm/trunk/SOURCE/advec_s_bc.f90	(revision 1353)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! REAL constants provided with KIND-attribute
 !
 ! Former revisions:
@@ -170,15 +170,15 @@
 !-- Array sk_p requires 2 extra elements for each dimension
     ALLOCATE( sk_p(nzb-2:nzt+3,nys-3:nyn+3,nxl-3:nxr+3) )
-    sk_p = 0.0
+    sk_p = 0.0_wp
 
 !
 !-- Assign reciprocal values in order to avoid divisions later
-    f2    = 0.5
-    f4    = 0.25
-    f8    = 0.125
-    f12   = 0.8333333333333333E-01
-    f24   = 0.4166666666666666E-01
-    f48   = 0.2083333333333333E-01
-    f1920 = 0.5208333333333333E-03
+    f2    = 0.5_wp
+    f4    = 0.25_wp
+    f8    = 0.125_wp
+    f12   = 0.8333333333333333E-01_wp
+    f24   = 0.4166666666666666E-01_wp
+    f48   = 0.2083333333333333E-01_wp
+    f1920 = 0.5208333333333333E-03_wp
 
 !
@@ -236,13 +236,13 @@
 !
 !-- Initialise control density
-    d = 0.0
+    d = 0.0_wp
 
 !
 !-- Determine maxima of the first and second derivative in x-direction
-    fmax_l = 0.0
+    fmax_l = 0.0_wp
     DO  i = nxl, nxr
        DO  j = nys, nyn
           DO  k = nzb+1, nzt
-             numera = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) )
+             numera = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) )
              denomi  = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
              fmax_l(1) = MAX( fmax_l(1) , numera )
@@ -258,5 +258,5 @@
 #endif
 
-    fmax = 0.04 * fmax
+    fmax = 0.04_wp * fmax
 
 !
@@ -271,5 +271,5 @@
               sw(nzb+1:nzt,nxl-1:nxr+1)                                        &
             )
-    imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
+    imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
 
 !
@@ -287,16 +287,16 @@
        DO  i = nxl-1, nxr+1
           DO  k = nzb+1, nzt
-             a12(k,i) = 0.5 * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
-             a22(k,i) = 0.5 * ( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i)              &
-                                              + sk_p(k,j,i-1) )
-             a0(k,i) = ( 9.0 * sk_p(k,j,i+2) - 116.0 * sk_p(k,j,i+1)           &
-                         + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k,j,i-1)        &
-                         + 9.0 * sk_p(k,j,i-2) ) * f1920
-             a1(k,i) = ( -5.0 * sk_p(k,j,i+2) + 34.0 * sk_p(k,j,i+1)           &
-                         - 34.0 * sk_p(k,j,i-1) + 5.0 * sk_p(k,j,i-2)          &
+             a12(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
+             a22(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i)        &
+                                                 + sk_p(k,j,i-1) )
+             a0(k,i) = ( 9.0_wp * sk_p(k,j,i+2)    - 116.0_wp * sk_p(k,j,i+1)  &
+                         + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j,i-1)  &
+                         + 9.0_wp * sk_p(k,j,i-2) ) * f1920
+             a1(k,i) = ( -5.0_wp * sk_p(k,j,i+2)   + 34.0_wp * sk_p(k,j,i+1)   &
+                         - 34.0_wp * sk_p(k,j,i-1) + 5.0_wp * sk_p(k,j,i-2)    &
                        ) * f48
-             a2(k,i) = ( -3.0 * sk_p(k,j,i+2) + 36.0 * sk_p(k,j,i+1)           &
-                         - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k,j,i-1)           &
-                         - 3.0 * sk_p(k,j,i-2) ) * f48
+             a2(k,i) = ( -3.0_wp * sk_p(k,j,i+2) + 36.0_wp * sk_p(k,j,i+1)     &
+                         - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j,i-1)     &
+                         - 3.0_wp * sk_p(k,j,i-2) ) * f48
           ENDDO
        ENDDO
@@ -309,31 +309,31 @@
              cip  =  MAX( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
              cim  = -MIN( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
-             cipf = 1.0 - 2.0 * cip
-             cimf = 1.0 - 2.0 * cim
-             ip   =   a0(k,i)   * f2  * ( 1.0 - cipf )                         &
-                    + a1(k,i)   * f8  * ( 1.0 - cipf*cipf )                    &
-                    + a2(k,i)   * f24 * ( 1.0 - cipf*cipf*cipf )
-             im   =   a0(k,i+1) * f2  * ( 1.0 - cimf )                         &
-                    - a1(k,i+1) * f8  * ( 1.0 - cimf*cimf )                    &
-                    + a2(k,i+1) * f24 * ( 1.0 - cimf*cimf*cimf )
+             cipf = 1.0_wp - 2.0_wp * cip
+             cimf = 1.0_wp - 2.0_wp * cim
+             ip   =   a0(k,i)   * f2  * ( 1.0_wp - cipf )                      &
+                    + a1(k,i)   * f8  * ( 1.0_wp - cipf*cipf )                 &
+                    + a2(k,i)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
+             im   =   a0(k,i+1) * f2  * ( 1.0_wp - cimf )                      &
+                    - a1(k,i+1) * f8  * ( 1.0_wp - cimf*cimf )                 &
+                    + a2(k,i+1) * f24 * ( 1.0_wp - cimf*cimf*cimf )
              ip   = MAX( ip, 0.0_wp )
              im   = MAX( im, 0.0_wp )
-             ippb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15) )
-             impb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i+1) / (ip+im+1E-15) )
+             ippb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
+             impb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i+1) / (ip+im+1E-15_wp) )
 
              cip  =  MAX( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
              cim  = -MIN( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
-             cipf = 1.0 - 2.0 * cip
-             cimf = 1.0 - 2.0 * cim
-             ip   =   a0(k,i-1) * f2  * ( 1.0 - cipf )                         &
-                    + a1(k,i-1) * f8  * ( 1.0 - cipf*cipf )                    &
-                    + a2(k,i-1) * f24 * ( 1.0 - cipf*cipf*cipf )
-             im   =   a0(k,i)   * f2  * ( 1.0 - cimf )                         &
-                    - a1(k,i)   * f8  * ( 1.0 - cimf*cimf )                    &
-                    + a2(k,i)   * f24 * ( 1.0 - cimf*cimf*cimf )
+             cipf = 1.0_wp - 2.0_wp * cip
+             cimf = 1.0_wp - 2.0_wp * cim
+             ip   =   a0(k,i-1) * f2  * ( 1.0_wp - cipf )                      &
+                    + a1(k,i-1) * f8  * ( 1.0_wp - cipf*cipf )                 &
+                    + a2(k,i-1) * f24 * ( 1.0_wp - cipf*cipf*cipf )
+             im   =   a0(k,i)   * f2  * ( 1.0_wp - cimf )                      &
+                    - a1(k,i)   * f8  * ( 1.0_wp - cimf*cimf )                 &
+                    + a2(k,i)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
              ip   = MAX( ip, 0.0_wp )
              im   = MAX( im, 0.0_wp )
-             ipmb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i-1) / (ip+im+1E-15) )
-             immb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15) )
+             ipmb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i-1) / (ip+im+1E-15_wp) )
+             immb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
           ENDDO
        ENDDO
@@ -343,13 +343,13 @@
        DO  i = nxl-2, nxr+2
           DO  k = nzb+1, nzt
-             m1z = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) )
+             m1z = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) )
              m1n = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
-             IF ( m1n /= 0.0  .AND.  m1n >= m1z )  THEN
+             IF ( m1n /= 0.0_wp  .AND.  m1n >= m1z )  THEN
                 m1(k,i) = m1z / m1n
-                IF ( m1(k,i) /= 2.0  .AND.  m1n < fmax(2) )  m1(k,i) = 0.0
+                IF ( m1(k,i) /= 2.0_wp  .AND.  m1n < fmax(2) )  m1(k,i) = 0.0_wp
              ELSEIF ( m1n < m1z )  THEN
-                m1(k,i) = -1.0
+                m1(k,i) = -1.0_wp
              ELSE
-                m1(k,i) = 0.0
+                m1(k,i) = 0.0_wp
              ENDIF
           ENDDO
@@ -358,25 +358,25 @@
 !
 !--    Compute switch sw
-       sw = 0.0
+       sw = 0.0_wp
        DO  i = nxl-1, nxr+1
           DO  k = nzb+1, nzt
-             m2 = 2.0 * ABS( a1(k,i) - a12(k,i) ) /                            &
+             m2 = 2.0_wp * ABS( a1(k,i) - a12(k,i) ) /                         &
                   MAX( ABS( a1(k,i) + a12(k,i) ), 1E-35_wp )
-             IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) )  m2 = 0.0
-
-             m3 = 2.0 * ABS( a2(k,i) - a22(k,i) ) /                            &
+             IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) )  m2 = 0.0_wp
+
+             m3 = 2.0_wp * ABS( a2(k,i) - a22(k,i) ) /                         &
                   MAX( ABS( a2(k,i) + a22(k,i) ), 1E-35_wp )
-             IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) )  m3 = 0.0
-
-             t1 = 0.35
-             t2 = 0.35
-             IF ( m1(k,i) == -1.0 )  t2 = 0.12
+             IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) )  m3 = 0.0_wp
+
+             t1 = 0.35_wp
+             t2 = 0.35_wp
+             IF ( m1(k,i) == -1.0_wp )  t2 = 0.12_wp
 
 !--          *VOCL STMT,IF(10)
-             IF ( m1(k,i-1) == 1.0 .OR. m1(k,i) == 1.0 .OR. m1(k,i+1) == 1.0   &
-                  .OR.  m2 > t2  .OR.  m3 > t2  .OR.                           &
-                  ( m1(k,i) > t1  .AND.  m1(k,i-1) /= -1.0  .AND.              &
-                    m1(k,i) /= -1.0  .AND.  m1(k,i+1) /= -1.0 )                &
-                )  sw(k,i) = 1.0
+             IF ( m1(k,i-1) == 1.0_wp .OR. m1(k,i) == 1.0_wp                   &
+                  .OR. m1(k,i+1) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
+                  ( m1(k,i) > t1  .AND.  m1(k,i-1) /= -1.0_wp  .AND.           &
+                    m1(k,i) /= -1.0_wp  .AND.  m1(k,i+1) /= -1.0_wp )          &
+                )  sw(k,i) = 1.0_wp
           ENDDO
        ENDDO
@@ -389,7 +389,7 @@
 
 !--          *VOCL STMT,IF(10)
-             IF ( sw(k,i) == 1.0 )  THEN
+             IF ( sw(k,i) == 1.0_wp )  THEN
                 snenn = sk_p(k,j,i+1) - sk_p(k,j,i-1)
-                IF ( ABS( snenn ) < 1E-9 )  snenn = 1E-9
+                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
                 sterm = ( sk_p(k,j,i) - sk_p(k,j,i-1) ) / snenn
                 sterm = MIN( sterm, 0.9999_wp )
@@ -402,12 +402,13 @@
                 ippe(k,i) = sk_p(k,j,i-1) * cip + snenn * (                    &
                             aex(ix) * cip + bex(ix) / dex(ix) * (              &
-                            eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) &
+                            eex(ix) -                                          &
+                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
                                                                 )              &
                                                           )
-                IF ( sterm == 0.0001 )  ippe(k,i) = sk_p(k,j,i) * cip
-                IF ( sterm == 0.9999 )  ippe(k,i) = sk_p(k,j,i) * cip
+                IF ( sterm == 0.0001_wp )  ippe(k,i) = sk_p(k,j,i) * cip
+                IF ( sterm == 0.9999_wp )  ippe(k,i) = sk_p(k,j,i) * cip
 
                 snenn = sk_p(k,j,i-1) - sk_p(k,j,i+1)
-                IF ( ABS( snenn ) < 1E-9 )  snenn = 1E-9
+                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
                 sterm = ( sk_p(k,j,i) - sk_p(k,j,i+1) ) / snenn
                 sterm = MIN( sterm, 0.9999_wp )
@@ -420,15 +421,16 @@
                 imme(k,i) = sk_p(k,j,i+1) * cim + snenn * (                    &
                             aex(ix) * cim + bex(ix) / dex(ix) * (              &
-                            eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) &
+                            eex(ix) -                                          &
+                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
                                                                 )              &
                                                           )
-                IF ( sterm == 0.0001 )  imme(k,i) = sk_p(k,j,i) * cim
-                IF ( sterm == 0.9999 )  imme(k,i) = sk_p(k,j,i) * cim
+                IF ( sterm == 0.0001_wp )  imme(k,i) = sk_p(k,j,i) * cim
+                IF ( sterm == 0.9999_wp )  imme(k,i) = sk_p(k,j,i) * cim
              ENDIF
 
 !--          *VOCL STMT,IF(10)
-             IF ( sw(k,i+1) == 1.0 )  THEN
+             IF ( sw(k,i+1) == 1.0_wp )  THEN
                 snenn = sk_p(k,j,i) - sk_p(k,j,i+2)
-                IF ( ABS( snenn ) .LT. 1E-9 )  snenn = 1E-9
+                IF ( ABS( snenn ) .LT. 1E-9_wp )  snenn = 1E-9_wp
                 sterm = ( sk_p(k,j,i+1) - sk_p(k,j,i+2) ) / snenn
                 sterm = MIN( sterm, 0.9999_wp )
@@ -441,5 +443,6 @@
                 impe(k,i) = sk_p(k,j,i+2) * cim + snenn * (                    &
                             aex(ix) * cim + bex(ix) / dex(ix) * (              &
-                            eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) &
+                            eex(ix) -                                          &
+                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
                                                                 )              &
                                                           )
@@ -449,7 +452,7 @@
 
 !--          *VOCL STMT,IF(10)
-             IF ( sw(k,i-1) == 1.0 )  THEN
+             IF ( sw(k,i-1) == 1.0_wp )  THEN
                 snenn = sk_p(k,j,i) - sk_p(k,j,i-2)
-                IF ( ABS( snenn ) < 1E-9 )  snenn = 1E-9
+                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
                 sterm = ( sk_p(k,j,i-1) - sk_p(k,j,i-2) ) / snenn
                 sterm = MIN( sterm, 0.9999_wp )
@@ -462,5 +465,6 @@
                 ipme(k,i) = sk_p(k,j,i-2) * cip + snenn * (                    &
                             aex(ix) * cip + bex(ix) / dex(ix) * (              &
-                            eex(ix) - EXP( dex(ix)*0.5 * ( 1.0_wp - 2.0 * cip ) ) &
+                            eex(ix) -                                          &
+                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
                                                                 )              &
                                                           )
@@ -538,9 +542,9 @@
 !
 !-- Determine the maxima of the first and second derivative in y-direction
-    fmax_l = 0.0
+    fmax_l = 0.0_wp
     DO  i = nxl, nxr
        DO  j = nys, nyn
           DO  k = nzb+1, nzt
-             numera = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) )
+             numera = ABS( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) )
              denomi  = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
              fmax_l(1) = MAX( fmax_l(1) , numera )
@@ -556,5 +560,5 @@
 #endif
 
-    fmax = 0.04 * fmax
+    fmax = 0.04_wp * fmax
 
 !
@@ -569,5 +573,5 @@
               sw(nzb+1:nzt,nys-1:nyn+1)                                        &
             )
-    imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
+    imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
 
 !
@@ -579,16 +583,16 @@
        DO  j = nys-1, nyn+1
           DO  k = nzb+1, nzt
-             a12(k,j) = 0.5 * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
-             a22(k,j) = 0.5 * ( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i)              &
-                                              + sk_p(k,j-1,i) )
-             a0(k,j) = ( 9.0 * sk_p(k,j+2,i) - 116.0 * sk_p(k,j+1,i)           &
-                         + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k,j-1,i)        &
-                         + 9.0 * sk_p(k,j-2,i) ) * f1920
-             a1(k,j) = ( -5.0 * sk_p(k,j+2,i) + 34.0 * sk_p(k,j+1,i)           &
-                         - 34.0 * sk_p(k,j-1,i) + 5.0 * sk_p(k,j-2,i)          &
+             a12(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
+             a22(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i)        &
+                                                 + sk_p(k,j-1,i) )
+             a0(k,j) = ( 9.0_wp * sk_p(k,j+2,i)    - 116.0_wp * sk_p(k,j+1,i)  &
+                         + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j-1,i)  &
+                         + 9.0_wp * sk_p(k,j-2,i) ) * f1920
+             a1(k,j) = ( -5.0_wp   * sk_p(k,j+2,i) + 34.0_wp * sk_p(k,j+1,i)   &
+                         - 34.0_wp * sk_p(k,j-1,i) + 5.0_wp  * sk_p(k,j-2,i)   &
                        ) * f48
-             a2(k,j) = ( -3.0 * sk_p(k,j+2,i) + 36.0 * sk_p(k,j+1,i)           &
-                         - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k,j-1,i)           &
-                         - 3.0 * sk_p(k,j-2,i) ) * f48
+             a2(k,j) = ( -3.0_wp * sk_p(k,j+2,i) + 36.0_wp * sk_p(k,j+1,i)     &
+                         - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j-1,i)     &
+                         - 3.0_wp * sk_p(k,j-2,i) ) * f48
           ENDDO
        ENDDO
@@ -601,31 +605,31 @@
              cip  =  MAX( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
              cim  = -MIN( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
-             cipf = 1.0 - 2.0 * cip
-             cimf = 1.0 - 2.0 * cim
-             ip   =   a0(k,j)   * f2  * ( 1.0 - cipf )                         &
-                    + a1(k,j)   * f8  * ( 1.0 - cipf*cipf )                    &
-                    + a2(k,j)   * f24 * ( 1.0 - cipf*cipf*cipf )
-             im   =   a0(k,j+1) * f2  * ( 1.0 - cimf )                         &
-                    - a1(k,j+1) * f8  * ( 1.0 - cimf*cimf )                    &
-                    + a2(k,j+1) * f24 * ( 1.0 - cimf*cimf*cimf )
+             cipf = 1.0_wp - 2.0_wp * cip
+             cimf = 1.0_wp - 2.0_wp * cim
+             ip   =   a0(k,j)   * f2  * ( 1.0_wp - cipf )                      &
+                    + a1(k,j)   * f8  * ( 1.0_wp - cipf*cipf )                 &
+                    + a2(k,j)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
+             im   =   a0(k,j+1) * f2  * ( 1.0_wp - cimf )                      &
+                    - a1(k,j+1) * f8  * ( 1.0_wp - cimf*cimf )                 &
+                    + a2(k,j+1) * f24 * ( 1.0_wp - cimf*cimf*cimf )
              ip   = MAX( ip, 0.0_wp )
              im   = MAX( im, 0.0_wp )
-             ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15) )
-             impb(k,j) = im * MIN( 1.0_wp, sk_p(k,j+1,i) / (ip+im+1E-15) )
+             ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
+             impb(k,j) = im * MIN( 1.0_wp, sk_p(k,j+1,i) / (ip+im+1E-15_wp) )
 
              cip  =  MAX( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
              cim  = -MIN( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
-             cipf = 1.0 - 2.0 * cip
-             cimf = 1.0 - 2.0 * cim
-             ip   =   a0(k,j-1) * f2  * ( 1.0 - cipf )                         &
-                    + a1(k,j-1) * f8  * ( 1.0 - cipf*cipf )                    &
-                    + a2(k,j-1) * f24 * ( 1.0 - cipf*cipf*cipf )
-             im   =   a0(k,j)   * f2  * ( 1.0 - cimf )                         &
-                    - a1(k,j)   * f8  * ( 1.0 - cimf*cimf )                    &
-                    + a2(k,j)   * f24 * ( 1.0 - cimf*cimf*cimf )
+             cipf = 1.0_wp - 2.0_wp * cip
+             cimf = 1.0_wp - 2.0_wp * cim
+             ip   =   a0(k,j-1) * f2  * ( 1.0_wp - cipf )                      &
+                    + a1(k,j-1) * f8  * ( 1.0_wp - cipf*cipf )                 &
+                    + a2(k,j-1) * f24 * ( 1.0_wp - cipf*cipf*cipf )
+             im   =   a0(k,j)   * f2  * ( 1.0_wp - cimf )                      &
+                    - a1(k,j)   * f8  * ( 1.0_wp - cimf*cimf )                 &
+                    + a2(k,j)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
              ip   = MAX( ip, 0.0_wp )
              im   = MAX( im, 0.0_wp )
-             ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j-1,i) / (ip+im+1E-15) )
-             immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15) )
+             ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j-1,i) / (ip+im+1E-15_wp) )
+             immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
           ENDDO
        ENDDO
@@ -635,13 +639,13 @@
        DO  j = nys-2, nyn+2
           DO  k = nzb+1, nzt
-             m1z = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) )
+             m1z = ABS( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) )
              m1n = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
-             IF ( m1n /= 0.0  .AND.  m1n >= m1z )  THEN
+             IF ( m1n /= 0.0_wp  .AND.  m1n >= m1z )  THEN
                 m1(k,j) = m1z / m1n
-                IF ( m1(k,j) /= 2.0  .AND.  m1n < fmax(2) )  m1(k,j) = 0.0
+                IF ( m1(k,j) /= 2.0_wp  .AND.  m1n < fmax(2) )  m1(k,j) = 0.0_wp
              ELSEIF ( m1n < m1z )  THEN
-                m1(k,j) = -1.0
+                m1(k,j) = -1.0_wp
              ELSE
-                m1(k,j) = 0.0
+                m1(k,j) = 0.0_wp
              ENDIF
           ENDDO
@@ -650,25 +654,25 @@
 !
 !--    Compute switch sw
-       sw = 0.0
+       sw = 0.0_wp
        DO  j = nys-1, nyn+1
           DO  k = nzb+1, nzt
-             m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) /                            &
+             m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) /                            &
                   MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp )
-             IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0
-
-             m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) /                            &
+             IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0_wp
+
+             m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) /                            &
                   MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp )
-             IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0
-
-             t1 = 0.35
-             t2 = 0.35
-             IF ( m1(k,j) == -1.0 )  t2 = 0.12
+             IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0_wp
+
+             t1 = 0.35_wp
+             t2 = 0.35_wp
+             IF ( m1(k,j) == -1.0_wp )  t2 = 0.12_wp
 
 !--          *VOCL STMT,IF(10)
-             IF ( m1(k,j-1) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k,j+1) == 1.0   &
-                  .OR.  m2 > t2  .OR.  m3 > t2  .OR.                           &
-                  ( m1(k,j) > t1  .AND.  m1(k,j-1) /= -1.0  .AND.              &
-                    m1(k,j) /= -1.0  .AND.  m1(k,j+1) /= -1.0 )                &
-                )  sw(k,j) = 1.0
+             IF ( m1(k,j-1) == 1.0_wp .OR. m1(k,j) == 1.0_wp                   &
+                  .OR. m1(k,j+1) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
+                  ( m1(k,j) > t1  .AND.  m1(k,j-1) /= -1.0_wp  .AND.           &
+                    m1(k,j) /= -1.0_wp  .AND.  m1(k,j+1) /= -1.0_wp )          &
+                )  sw(k,j) = 1.0_wp
           ENDDO
        ENDDO
@@ -681,7 +685,7 @@
 
 !--          *VOCL STMT,IF(10)
-             IF ( sw(k,j) == 1.0 )  THEN
+             IF ( sw(k,j) == 1.0_wp )  THEN
                 snenn = sk_p(k,j+1,i) - sk_p(k,j-1,i)
-                IF ( ABS( snenn ) < 1E-9 )  snenn = 1E-9
+                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
                 sterm = ( sk_p(k,j,i) - sk_p(k,j-1,i) ) / snenn
                 sterm = MIN( sterm, 0.9999_wp )
@@ -694,5 +698,6 @@
                 ippe(k,j) = sk_p(k,j-1,i) * cip + snenn * (                    &
                             aex(ix) * cip + bex(ix) / dex(ix) * (              &
-                            eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) &
+                            eex(ix) -                                          &
+                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
                                                                 )              &
                                                           )
@@ -701,5 +706,5 @@
 
                 snenn = sk_p(k,j-1,i) - sk_p(k,j+1,i)
-                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9
+                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
                 sterm = ( sk_p(k,j,i) - sk_p(k,j+1,i) ) / snenn
                 sterm = MIN( sterm, 0.9999_wp )
@@ -712,5 +717,6 @@
                 imme(k,j) = sk_p(k,j+1,i) * cim + snenn * (                    &
                             aex(ix) * cim + bex(ix) / dex(ix) * (              &
-                            eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) &
+                            eex(ix) -                                          &
+                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
                                                                 )              &
                                                           )
@@ -720,7 +726,7 @@
 
 !--          *VOCL STMT,IF(10)
-             IF ( sw(k,j+1) == 1.0 )  THEN
+             IF ( sw(k,j+1) == 1.0_wp )  THEN
                 snenn = sk_p(k,j,i) - sk_p(k,j+2,i)
-                IF ( ABS( snenn ) .LT. 1E-9 )  snenn = 1E-9
+                IF ( ABS( snenn ) .LT. 1E-9_wp )  snenn = 1E-9_wp
                 sterm = ( sk_p(k,j+1,i) - sk_p(k,j+2,i) ) / snenn
                 sterm = MIN( sterm, 0.9999_wp )
@@ -733,5 +739,6 @@
                 impe(k,j) = sk_p(k,j+2,i) * cim + snenn * (                    &
                             aex(ix) * cim + bex(ix) / dex(ix) * (              &
-                            eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) &
+                            eex(ix) -                                          &
+                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
                                                                 )              &
                                                           )
@@ -741,7 +748,7 @@
 
 !--          *VOCL STMT,IF(10)
-             IF ( sw(k,j-1) == 1.0 )  THEN
+             IF ( sw(k,j-1) == 1.0_wp )  THEN
                 snenn = sk_p(k,j,i) - sk_p(k,j-2,i)
-                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9
+                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
                 sterm = ( sk_p(k,j-1,i) - sk_p(k,j-2,i) ) / snenn
                 sterm = MIN( sterm, 0.9999_wp )
@@ -754,5 +761,6 @@
                 ipme(k,j) = sk_p(k,j-2,i) * cip + snenn * (                    &
                             aex(ix) * cip + bex(ix) / dex(ix) * (              &
-                            eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) &
+                            eex(ix) -                                          &
+                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
                                                                 )              &
                                                           )
@@ -769,8 +777,8 @@
        DO  j = nys, nyn
           DO  k = nzb+1, nzt
-             fplus  = ( 1.0 - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j)  &
-                    - ( 1.0 - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j)
-             fminus = ( 1.0 - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j)  &
-                    - ( 1.0 - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
+             fplus  = ( 1.0_wp - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
+                    - ( 1.0_wp - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j)
+             fminus = ( 1.0_wp - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j) &
+                    - ( 1.0_wp - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
              tendcy = fplus - fminus
 !
@@ -781,6 +789,6 @@
 !--          Density correction because of possible remaining divergences
              d_new = d(k,j,i) - ( v(k,j+1,i) - v(k,j,i) ) * dt_3d * ddy
-             sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /    &
-                           ( 1.0 + d_new )
+             sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /  &
+                           ( 1.0_wp + d_new )
              d(k,j,i)  = d_new
           ENDDO
@@ -800,5 +808,5 @@
 !-- Initialise for the computation of heat fluxes (see below; required in
 !-- UP flow_statistics)
-    IF ( sk_char == 'pt' )  sums_wsts_bc_l = 0.0
+    IF ( sk_char == 'pt' )  sums_wsts_bc_l = 0.0_wp
 
 !
@@ -939,9 +947,9 @@
 !
 !-- Determine the maxima of the first and second derivative in z-direction
-    fmax_l = 0.0
+    fmax_l = 0.0_wp
     DO  i = nxl, nxr
        DO  j = nys, nyn
           DO  k = nzb, nzt+1
-             numera = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )
+             numera = ABS( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) )
              denomi  = ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) )
              fmax_l(1) = MAX( fmax_l(1) , numera )
@@ -957,5 +965,5 @@
 #endif
 
-    fmax = 0.04 * fmax
+    fmax = 0.04_wp * fmax
 
 !
@@ -970,5 +978,5 @@
               sw(nzb:nzt+1,nys:nyn)                                            &
             )
-    imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
+    imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
 
 !
@@ -980,16 +988,16 @@
        DO  j = nys, nyn
           DO  k = nzb, nzt+1
-             a12(k,j) = 0.5 * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
-             a22(k,j) = 0.5 * ( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i)              &
-                                              + sk_p(k-1,j,i) )
-             a0(k,j) = ( 9.0 * sk_p(k+2,j,i) - 116.0 * sk_p(k+1,j,i)           &
-                         + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k-1,j,i)        &
-                         + 9.0 * sk_p(k-2,j,i) ) * f1920
-             a1(k,j) = ( -5.0 * sk_p(k+2,j,i) + 34.0 * sk_p(k+1,j,i)           &
-                         - 34.0 * sk_p(k-1,j,i) + 5.0 * sk_p(k-2,j,i)          &
+             a12(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
+             a22(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i)        &
+                                                 + sk_p(k-1,j,i) )
+             a0(k,j) = ( 9.0_wp * sk_p(k+2,j,i)    - 116.0_wp * sk_p(k+1,j,i)  &
+                         + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k-1,j,i)  &
+                         + 9.0_wp * sk_p(k-2,j,i) ) * f1920
+             a1(k,j) = ( -5.0_wp   * sk_p(k+2,j,i) + 34.0_wp * sk_p(k+1,j,i)   &
+                         - 34.0_wp * sk_p(k-1,j,i) + 5.0_wp  * sk_p(k-2,j,i)   &
                        ) * f48
-             a2(k,j) = ( -3.0 * sk_p(k+2,j,i) + 36.0 * sk_p(k+1,j,i)           &
-                         - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k-1,j,i)           &
-                         - 3.0 * sk_p(k-2,j,i) ) * f48
+             a2(k,j) = ( -3.0_wp * sk_p(k+2,j,i) + 36.0_wp * sk_p(k+1,j,i)     &
+                         - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k-1,j,i)     &
+                         - 3.0_wp * sk_p(k-2,j,i) ) * f48
           ENDDO
        ENDDO
@@ -1002,31 +1010,31 @@
              cip  =  MAX( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
              cim  = -MIN( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
-             cipf = 1.0 - 2.0 * cip
-             cimf = 1.0 - 2.0 * cim
-             ip   =   a0(k,j)   * f2  * ( 1.0 - cipf )                         &
-                    + a1(k,j)   * f8  * ( 1.0 - cipf*cipf )                    &
-                    + a2(k,j)   * f24 * ( 1.0 - cipf*cipf*cipf )
-             im   =   a0(k+1,j) * f2  * ( 1.0 - cimf )                         &
-                    - a1(k+1,j) * f8  * ( 1.0 - cimf*cimf )                    &
-                    + a2(k+1,j) * f24 * ( 1.0 - cimf*cimf*cimf )
+             cipf = 1.0_wp - 2.0_wp * cip
+             cimf = 1.0_wp - 2.0_wp * cim
+             ip   =   a0(k,j)   * f2  * ( 1.0_wp - cipf )                      &
+                    + a1(k,j)   * f8  * ( 1.0_wp - cipf*cipf )                 &
+                    + a2(k,j)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
+             im   =   a0(k+1,j) * f2  * ( 1.0_wp - cimf )                      &
+                    - a1(k+1,j) * f8  * ( 1.0_wp - cimf*cimf )                 &
+                    + a2(k+1,j) * f24 * ( 1.0_wp - cimf*cimf*cimf )
              ip   = MAX( ip, 0.0_wp )
              im   = MAX( im, 0.0_wp )
-             ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15) )
-             impb(k,j) = im * MIN( 1.0_wp, sk_p(k+1,j,i) / (ip+im+1E-15) )
+             ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
+             impb(k,j) = im * MIN( 1.0_wp, sk_p(k+1,j,i) / (ip+im+1E-15_wp) )
 
              cip  =  MAX( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
              cim  = -MIN( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
-             cipf = 1.0 - 2.0 * cip
-             cimf = 1.0 - 2.0 * cim
-             ip   =   a0(k-1,j) * f2  * ( 1.0 - cipf )                         &
-                    + a1(k-1,j) * f8  * ( 1.0 - cipf*cipf )                    &
-                    + a2(k-1,j) * f24 * ( 1.0 - cipf*cipf*cipf )
-             im   =   a0(k,j)   * f2  * ( 1.0 - cimf )                         &
-                    - a1(k,j)   * f8  * ( 1.0 - cimf*cimf )                    &
-                    + a2(k,j)   * f24 * ( 1.0 - cimf*cimf*cimf )
+             cipf = 1.0_wp - 2.0_wp * cip
+             cimf = 1.0_wp - 2.0_wp * cim
+             ip   =   a0(k-1,j) * f2  * ( 1.0_wp - cipf )                      &
+                    + a1(k-1,j) * f8  * ( 1.0_wp - cipf*cipf )                 &
+                    + a2(k-1,j) * f24 * ( 1.0_wp - cipf*cipf*cipf )
+             im   =   a0(k,j)   * f2  * ( 1.0_wp - cimf )                      &
+                    - a1(k,j)   * f8  * ( 1.0_wp - cimf*cimf )                 &
+                    + a2(k,j)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
              ip   = MAX( ip, 0.0_wp )
              im   = MAX( im, 0.0_wp )
-             ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k-1,j,i) / (ip+im+1E-15) )
-             immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15) )
+             ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k-1,j,i) / (ip+im+1E-15_wp) )
+             immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
           ENDDO
        ENDDO
@@ -1036,13 +1044,13 @@
        DO  j = nys, nyn
           DO  k = nzb-1, nzt+2
-             m1z = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )
+             m1z = ABS( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) )
              m1n = ABS( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
-             IF ( m1n /= 0.0  .AND.  m1n >= m1z )  THEN
+             IF ( m1n /= 0.0_wp  .AND.  m1n >= m1z )  THEN
                 m1(k,j) = m1z / m1n
-                IF ( m1(k,j) /= 2.0  .AND.  m1n < fmax(2) )  m1(k,j) = 0.0
+                IF ( m1(k,j) /= 2.0_wp  .AND.  m1n < fmax(2) )  m1(k,j) = 0.0_wp
              ELSEIF ( m1n < m1z )  THEN
-                m1(k,j) = -1.0
+                m1(k,j) = -1.0_wp
              ELSE
-                m1(k,j) = 0.0
+                m1(k,j) = 0.0_wp
              ENDIF
           ENDDO
@@ -1051,25 +1059,25 @@
 !
 !--    Compute switch sw
-       sw = 0.0
+       sw = 0.0_wp
        DO  j = nys, nyn
           DO  k = nzb, nzt+1
-             m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) /                            &
+             m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) /                         &
                   MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp )
-             IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0
-
-             m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) /                            &
+             IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0_wp
+
+             m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) /                         &
                   MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp )
-             IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0
-
-             t1 = 0.35
-             t2 = 0.35
-             IF ( m1(k,j) == -1.0 )  t2 = 0.12
+             IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0_wp
+
+             t1 = 0.35_wp
+             t2 = 0.35_wp
+             IF ( m1(k,j) == -1.0_wp )  t2 = 0.12_wp
 
 !--          *VOCL STMT,IF(10)
-             IF ( m1(k-1,j) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k+1,j) == 1.0   &
-                  .OR.  m2 > t2  .OR.  m3 > t2  .OR.                           &
-                  ( m1(k,j) > t1  .AND.  m1(k-1,j) /= -1.0  .AND.              &
-                    m1(k,j) /= -1.0  .AND.  m1(k+1,j) /= -1.0 )                &
-                )  sw(k,j) = 1.0
+             IF ( m1(k-1,j) == 1.0_wp .OR. m1(k,j) == 1.0_wp                   &
+                  .OR. m1(k+1,j) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
+                  ( m1(k,j) > t1  .AND.  m1(k-1,j) /= -1.0_wp  .AND.           &
+                    m1(k,j) /= -1.0_wp  .AND.  m1(k+1,j) /= -1.0_wp )          &
+                )  sw(k,j) = 1.0_wp
           ENDDO
        ENDDO
@@ -1082,7 +1090,7 @@
 
 !--          *VOCL STMT,IF(10)
-             IF ( sw(k,j) == 1.0 )  THEN
+             IF ( sw(k,j) == 1.0_wp )  THEN
                 snenn = sk_p(k+1,j,i) - sk_p(k-1,j,i)
-                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9
+                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
                 sterm = ( sk_p(k,j,i) - sk_p(k-1,j,i) ) / snenn
                 sterm = MIN( sterm, 0.9999_wp )
@@ -1095,12 +1103,13 @@
                 ippe(k,j) = sk_p(k-1,j,i) * cip + snenn * (                    &
                             aex(ix) * cip + bex(ix) / dex(ix) * (              &
-                            eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) &
+                            eex(ix) -                                          &
+                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
                                                                 )              &
                                                           )
-                IF ( sterm == 0.0001 )  ippe(k,j) = sk_p(k,j,i) * cip
-                IF ( sterm == 0.9999 )  ippe(k,j) = sk_p(k,j,i) * cip
+                IF ( sterm == 0.0001_wp )  ippe(k,j) = sk_p(k,j,i) * cip
+                IF ( sterm == 0.9999_wp )  ippe(k,j) = sk_p(k,j,i) * cip
 
                 snenn = sk_p(k-1,j,i) - sk_p(k+1,j,i)
-                IF ( ABS( snenn ) < 1E-9 )  snenn = 1E-9
+                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
                 sterm = ( sk_p(k,j,i) - sk_p(k+1,j,i) ) / snenn
                 sterm = MIN( sterm, 0.9999_wp )
@@ -1113,5 +1122,6 @@
                 imme(k,j) = sk_p(k+1,j,i) * cim + snenn * (                    &
                             aex(ix) * cim + bex(ix) / dex(ix) * (              &
-                            eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) &
+                            eex(ix) -                                          &
+                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) &
                                                                 )              &
                                                           )
@@ -1121,7 +1131,7 @@
 
 !--          *VOCL STMT,IF(10)
-             IF ( sw(k+1,j) == 1.0 )  THEN
+             IF ( sw(k+1,j) == 1.0_wp )  THEN
                 snenn = sk_p(k,j,i) - sk_p(k+2,j,i)
-                IF ( ABS( snenn ) .LT. 1E-9 )  snenn = 1E-9
+                IF ( ABS( snenn ) .LT. 1E-9_wp )  snenn = 1E-9_wp
                 sterm = ( sk_p(k+1,j,i) - sk_p(k+2,j,i) ) / snenn
                 sterm = MIN( sterm, 0.9999_wp )
@@ -1134,5 +1144,6 @@
                 impe(k,j) = sk_p(k+2,j,i) * cim + snenn * (                    &
                             aex(ix) * cim + bex(ix) / dex(ix) * (              &
-                            eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) &
+                            eex(ix) -                                          &
+                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
                                                                 )              &
                                                           )
@@ -1142,7 +1153,7 @@
 
 !--          *VOCL STMT,IF(10)
-             IF ( sw(k-1,j) == 1.0 )  THEN
+             IF ( sw(k-1,j) == 1.0_wp )  THEN
                 snenn = sk_p(k,j,i) - sk_p(k-2,j,i)
-                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9
+                IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
                 sterm = ( sk_p(k-1,j,i) - sk_p(k-2,j,i) ) / snenn
                 sterm = MIN( sterm, 0.9999_wp )
@@ -1155,5 +1166,6 @@
                 ipme(k,j) = sk_p(k-2,j,i) * cip + snenn * (                    &
                             aex(ix) * cip + bex(ix) / dex(ix) * (              &
-                            eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) &
+                            eex(ix) -                                          &
+                            EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
                                                                 )              &
                                                           )
@@ -1170,8 +1182,8 @@
        DO  j = nys, nyn
           DO  k = nzb+1, nzt
-             fplus  = ( 1.0 - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j)  &
-                    - ( 1.0 - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j)
-             fminus = ( 1.0 - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j)  &
-                    - ( 1.0 - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
+             fplus  = ( 1.0_wp - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
+                    - ( 1.0_wp - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j)
+             fminus = ( 1.0_wp - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) &
+                    - ( 1.0_wp - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
              tendcy = fplus - fminus
 !
@@ -1182,6 +1194,6 @@
 !--          Density correction because of possible remaining divergences
              d_new = d(k,j,i) - ( w(k,j,i) - w(k-1,j,i) ) * dt_3d * ddzw(k)
-             sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /    &
-                           ( 1.0 + d_new )
+             sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /  &
+                           ( 1.0_wp + d_new )
 !
 !--          Store heat flux for subsequent statistics output.
