Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 1323)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 1353)
@@ -21,5 +21,5 @@
 ! Current revisions:
 ! -----------------
-!
+! REAL constants provided with KIND-attribute 
 !
 ! Former revisions:
@@ -170,5 +170,5 @@
 !
 !--    Initialize (local) summation array
-       sums_l = 0.0
+       sums_l = 0.0_wp
 
 !
@@ -202,9 +202,9 @@
              sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2 
              sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2 
-             sums_l(:,34,i) = sums_l(:,34,i) + 0.5 *                           & 
+             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        & 
                               ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
                                 sums_ws2_ws_l(:,i) )    ! e*
              DO  k = nzb, nzt
-                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5 * (    &
+                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5_wp * ( &
                                                       sums_us2_ws_l(k,i) +     &
                                                       sums_vs2_ws_l(k,i) +     &
@@ -446,5 +446,5 @@
        DO  i = nxl, nxr
           DO  j =  nys, nyn
-             sums_l_etot = 0.0
+             sums_l_etot = 0.0_wp
              DO  k = nzb_s_inner(j,i), nzt+1
 !
@@ -470,5 +470,5 @@
 
                 sums_l_etot  = sums_l_etot + &
-                                        0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 +    &
+                                        0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 + &
                                         w(k,j,i)**2 ) * rmask(j,i,sr)
              ENDDO
@@ -503,5 +503,5 @@
           DO  i = nxl, nxr
              DO  j =  nys, nyn
-                sums_l_eper = 0.0
+                sums_l_eper = 0.0_wp
                 DO  k = nzb_s_inner(j,i), nzt+1
                    u2   = u(k,j,i)**2
@@ -517,11 +517,11 @@
 !--             Perturbation energy
 
-                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 *       &
+                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5_wp *                &
                                   ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
-                   sums_l_eper  = sums_l_eper + &
-                                        0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
-
-                ENDDO
-                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)   &
+                   sums_l_eper     = sums_l_eper +                             &
+                                     0.5_wp * ( ust2+vst2+w2 ) * rmask(j,i,sr)
+
+                ENDDO
+                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)            &
                      + sums_l_eper
              ENDDO
@@ -547,5 +547,5 @@
 !
 !--             Momentum flux w"u"
-                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25 * (                   &
+                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25_wp * (                &
                                km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
                                                            ) * (               &
@@ -555,5 +555,5 @@
 !
 !--             Momentum flux w"v"
-                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25 * (                   &
+                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25_wp * (                &
                                km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
                                                            ) * (               &
@@ -564,5 +564,5 @@
 !--             Heat flux w"pt"
                 sums_l(k,16,tn) = sums_l(k,16,tn)                              &
-                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
+                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
                                                * ( pt(k+1,j,i) - pt(k,j,i) )   &
                                                * ddzu(k+1) * rmask(j,i,sr)
@@ -573,5 +573,5 @@
                 IF ( ocean )  THEN
                    sums_l(k,65,tn) = sums_l(k,65,tn)                           &
-                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
+                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
                                                * ( sa(k+1,j,i) - sa(k,j,i) )   &
                                                * ddzu(k+1) * rmask(j,i,sr)
@@ -582,9 +582,9 @@
                 IF ( humidity ) THEN
                    sums_l(k,45,tn) = sums_l(k,45,tn)                           &
-                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
+                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
                                                * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
                                                * ddzu(k+1) * rmask(j,i,sr)
                    sums_l(k,48,tn) = sums_l(k,48,tn)                           &
-                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
+                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
                                                * ( q(k+1,j,i) - q(k,j,i) )     &
                                                * ddzu(k+1) * rmask(j,i,sr)
@@ -592,5 +592,5 @@
                    IF ( cloud_physics ) THEN
                       sums_l(k,51,tn) = sums_l(k,51,tn)                        &
-                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
+                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
                                                * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
                                                 - ( q(k,j,i) - ql(k,j,i) ) )   &
@@ -603,5 +603,5 @@
                 IF ( passive_scalar )  THEN
                    sums_l(k,48,tn) = sums_l(k,48,tn)                           &
-                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
+                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
                                                * ( q(k+1,j,i) - q(k,j,i) )     &
                                                * ddzu(k+1) * rmask(j,i,sr)
@@ -620,7 +620,7 @@
                                     shf(j,i)  * rmask(j,i,sr)     ! w"pt"
                 sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
-                                    0.0 * rmask(j,i,sr)           ! u"pt"
+                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
                 sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
-                                    0.0 * rmask(j,i,sr)           ! v"pt"
+                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
                 IF ( ocean )  THEN
                    sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
@@ -628,15 +628,15 @@
                 ENDIF
                 IF ( humidity )  THEN
-                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
+                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
                                        qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
-                   sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (           &
-                                       ( 1.0 + 0.61 * q(nzb,j,i) ) *   &
-                                       shf(j,i) + 0.61 * pt(nzb,j,i) * &
+                   sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
+                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
+                                       shf(j,i) + 0.61_wp * pt(nzb,j,i) *      &
                                                   qsws(j,i) )
                    IF ( cloud_droplets )  THEN
-                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (        &
-                                         ( 1.0 + 0.61 * q(nzb,j,i) -   &
-                                           ql(nzb,j,i) ) * shf(j,i) +  &
-                                           0.61 * pt(nzb,j,i) * qsws(j,i) )
+                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
+                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
+                                           ql(nzb,j,i) ) * shf(j,i) +          &
+                                           0.61_wp * pt(nzb,j,i) * qsws(j,i) )
                    ENDIF
                    IF ( cloud_physics )  THEN
@@ -663,7 +663,7 @@
                                     tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
                 sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
-                                    0.0 * rmask(j,i,sr)           ! u"pt"
+                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
                 sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
-                                    0.0 * rmask(j,i,sr)           ! v"pt" 
+                                    0.0_wp * rmask(j,i,sr)        ! v"pt" 
 
                 IF ( ocean )  THEN
@@ -672,15 +672,15 @@
                 ENDIF
                 IF ( humidity )  THEN
-                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
+                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
                                        qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
-                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (             &
-                                       ( 1.0 + 0.61 * q(nzt,j,i) ) *     &
-                                       tswst(j,i) + 0.61 * pt(nzt,j,i) * &
-                                                  qswst(j,i) )
+                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                   &
+                                       ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *     &
+                                       tswst(j,i) + 0.61_wp * pt(nzt,j,i) *    &
+                                                             qswst(j,i) )
                    IF ( cloud_droplets )  THEN
-                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
-                                          ( 1.0 + 0.61 * q(nzt,j,i) -     &
-                                            ql(nzt,j,i) ) * tswst(j,i) +  &
-                                            0.61 * pt(nzt,j,i) * qswst(j,i) )
+                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                &
+                                          ( 1.0_wp + 0.61_wp * q(nzt,j,i) -    &
+                                            ql(nzt,j,i) ) * tswst(j,i) +       &
+                                            0.61_wp * pt(nzt,j,i) * qswst(j,i) )
                    ENDIF
                    IF ( cloud_physics )  THEN
@@ -703,15 +703,15 @@
 !--                rearranged according to the staggered grid.
              DO  k = nzb_s_inner(j,i), nzt
-                ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
-                              u(k+1,j,i) - hom(k+1,1,1,sr) )
-                vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
-                              v(k+1,j,i) - hom(k+1,1,2,sr) )
-                pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
-                              pt(k+1,j,i) - hom(k+1,1,4,sr) )
+                ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                  &
+                                 u(k+1,j,i) - hom(k+1,1,1,sr) )
+                vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                  &
+                                 v(k+1,j,i) - hom(k+1,1,2,sr) )
+                pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +                 &
+                                 pt(k+1,j,i) - hom(k+1,1,4,sr) )
 
 !--             Higher moments
-                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 * &
+                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 *        &
                                                     rmask(j,i,sr)
-                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) * &
+                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) *        &
                                                     rmask(j,i,sr)
 
@@ -721,12 +721,12 @@
                 IF ( ocean )  THEN
                    IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
-                      pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
-                                 sa(k+1,j,i) - hom(k+1,1,23,sr) )
-                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
+                      pts = 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +          &
+                                       sa(k+1,j,i) - hom(k+1,1,23,sr) )
+                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) *     &
+                                        rmask(j,i,sr)
+                   ENDIF
+                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) *            &
                                                        rmask(j,i,sr)
-                   ENDIF
-                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
-                                                       rmask(j,i,sr)
-                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) * &
+                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) *           &
                                                        rmask(j,i,sr)
                 ENDIF
@@ -737,10 +737,10 @@
                 IF ( humidity )  THEN
                    IF ( cloud_physics .OR. cloud_droplets )  THEN
-                      pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +        &
+                      pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +         &
                                     vpt(k+1,j,i) - hom(k+1,1,44,sr) )
-                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
+                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
                                                           rmask(j,i,sr)
                       IF ( .NOT. cloud_droplets )  THEN
-                         pts = 0.5 *                                           &
+                         pts = 0.5_wp *                                        &
                               ( ( q(k,j,i) - ql(k,j,i) ) -                     &
                               hom(k,1,42,sr) +                                 &
@@ -772,12 +772,14 @@
                    ELSE
                       IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
-                         pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
-                                       vpt(k+1,j,i) - hom(k+1,1,44,sr) )
-                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
+                         pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +      &
+                                          vpt(k+1,j,i) - hom(k+1,1,44,sr) )
+                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
                                                              rmask(j,i,sr)
                       ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
-                         sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) *  &
-                                             sums_l(k,17,tn) +      &
-                                          0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
+                         sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp *                & 
+                                             hom(k,1,41,sr) ) *                &
+                                           sums_l(k,17,tn) +                   &
+                                           0.61_wp * hom(k,1,4,sr) *           &
+                                           sums_l(k,49,tn)
                       END IF
                    END IF
@@ -785,9 +787,9 @@
 !
 !--             Passive scalar flux
-                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca   &
+                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca                &
                      .OR. sr /= 0 ) )  THEN
-                   pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
-                                 q(k+1,j,i) - hom(k+1,1,41,sr) )
-                   sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
+                   pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +              &
+                                    q(k+1,j,i) - hom(k+1,1,41,sr) )
+                   sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *        &
                                                        rmask(j,i,sr)
                 ENDIF
@@ -796,6 +798,6 @@
 !--             Energy flux w*e*
 !--             has to be adjusted
-                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *          &
-                                             ( ust**2 + vst**2 + w(k,j,i)**2 )&
+                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp *        &
+                                             ( ust**2 + vst**2 + w(k,j,i)**2 ) &
                                              * rmask(j,i,sr)
              ENDDO
@@ -811,17 +813,17 @@
             DO  j = nys, nyn
                DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
-                  ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
-                              u(k+1,j,i) - hom(k+1,1,1,sr) )
-                  vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
-                              v(k+1,j,i) - hom(k+1,1,2,sr) )
+                  ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                &
+                                   u(k+1,j,i) - hom(k+1,1,1,sr) )
+                  vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                &
+                                   v(k+1,j,i) - hom(k+1,1,2,sr) )
 !
 !--               Momentum flux w*u*
-                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                   &
-                                                    ( w(k,j,i-1) + w(k,j,i) ) &
+                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                 &
+                                                    ( w(k,j,i-1) + w(k,j,i) )  &
                                                     * ust * rmask(j,i,sr)
 !
 !--               Momentum flux w*v*
-                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 *                   &
-                                                    ( w(k,j-1,i) + w(k,j,i) ) &
+                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                 &
+                                                    ( w(k,j-1,i) + w(k,j,i) )  &
                                                     * vst * rmask(j,i,sr)
                ENDDO
@@ -837,13 +839,13 @@
 !
 !--               Vertical heat flux
-                  sums_l(k,17,tn) = sums_l(k,17,tn) +  0.5 * &
-                           ( pt(k,j,i)   - hom(k,1,4,sr) + &
-                           pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
+                  sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp *                 &
+                           ( pt(k,j,i)   - hom(k,1,4,sr) +                     &
+                             pt(k+1,j,i) - hom(k+1,1,4,sr) )                   &
                            * w(k,j,i) * rmask(j,i,sr)
                   IF ( humidity )  THEN
-                     pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
-                                q(k+1,j,i) - hom(k+1,1,41,sr) )
-                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
-                                                      rmask(j,i,sr)
+                     pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +            &
+                                      q(k+1,j,i) - hom(k+1,1,41,sr) )
+                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *      &
+                                       rmask(j,i,sr)
                   ENDIF
                ENDDO
@@ -865,7 +867,7 @@
 !--    First calculate the products, then the divergence.
 !--    Calculation is time consuming. Do it only, if profiles shall be plotted.
-       IF ( hom(nzb+1,2,55,0) /= 0.0  .OR.  hom(nzb+1,2,68,0) /= 0.0 )  THEN
-
-          sums_ll = 0.0  ! local array
+       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )  THEN
+
+          sums_ll = 0.0_wp  ! local array
 
           !$OMP DO
@@ -874,14 +876,14 @@
                 DO  k = nzb_s_inner(j,i)+1, nzt
 
-                   sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * (        &
-                  ( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)   &
-                           - 0.5 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
-                           ) )**2                                          &
-                + ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)   &
-                           - 0.5 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
-                           ) )**2                                          &
-                   + w(k,j,i)**2                                  )
-
-                   sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
+                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
+                  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)    &
+                            - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
+                              ) )**2                                           &
+                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)    &
+                            - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
+                              ) )**2                                           &
+                + w(k,j,i)**2                                        )
+
+                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
                                                * ( p(k,j,i) + p(k+1,j,i) )
 
@@ -889,8 +891,8 @@
              ENDDO
           ENDDO
-          sums_ll(0,1)     = 0.0    ! because w is zero at the bottom
-          sums_ll(nzt+1,1) = 0.0
-          sums_ll(0,2)     = 0.0
-          sums_ll(nzt+1,2) = 0.0
+          sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
+          sums_ll(nzt+1,1) = 0.0_wp
+          sums_ll(0,2)     = 0.0_wp
+          sums_ll(nzt+1,2) = 0.0_wp
 
           DO  k = nzb+1, nzt
@@ -901,5 +903,5 @@
           sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
           sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
-          sums_l(nzb,68,tn) = 0.0    ! because w* = 0 at nzb
+          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
 
        ENDIF
@@ -907,5 +909,5 @@
 !
 !--    Divergence of vertical flux of SGS TKE and the flux itself (69)
-       IF ( hom(nzb+1,2,57,0) /= 0.0  .OR.  hom(nzb+1,2,69,0) /= 0.0 )  THEN
+       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )  THEN
 
           !$OMP DO
@@ -914,12 +916,12 @@
                 DO  k = nzb_s_inner(j,i)+1, nzt
 
-                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * (                 &
+                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
                    (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
                  - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
-                                                             ) * ddzw(k)
-
-                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * (                 &
+                                                                ) * ddzw(k)
+
+                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
                    (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
-                                                              )
+                                                                )
 
                 ENDDO
@@ -934,5 +936,5 @@
 !--    Horizontal heat fluxes (subgrid, resolved, total). 
 !--    Do it only, if profiles shall be plotted. 
-       IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
+       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
 
           !$OMP DO
@@ -942,9 +944,9 @@
 !
 !--                Subgrid horizontal heat fluxes u"pt", v"pt"
-                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 *                   &
+                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
                                                    ( kh(k,j,i) + kh(k,j,i-1) ) &
                                                  * ( pt(k,j,i-1) - pt(k,j,i) ) &
                                                  * ddx * rmask(j,i,sr)
-                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 *                   &
+                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
                                                    ( kh(k,j,i) + kh(k,j-1,i) ) &
                                                  * ( pt(k,j-1,i) - pt(k,j,i) ) &
@@ -954,11 +956,11 @@
                    sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
                                                   ( u(k,j,i) - hom(k,1,1,sr) ) &
-                                       * 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
+                                    * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
                                                  pt(k,j,i)   - hom(k,1,4,sr) )
-                   pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
-                                 pt(k,j,i)   - hom(k,1,4,sr) )
+                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
+                                    pt(k,j,i)   - hom(k,1,4,sr) )
                    sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
                                                   ( v(k,j,i) - hom(k,1,2,sr) ) &
-                                       * 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
+                                    * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
                                                  pt(k,j,i)   - hom(k,1,4,sr) )
                 ENDDO
@@ -967,10 +969,10 @@
 !
 !--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
-          sums_l(nzb,58,tn) = 0.0
-          sums_l(nzb,59,tn) = 0.0
-          sums_l(nzb,60,tn) = 0.0
-          sums_l(nzb,61,tn) = 0.0
-          sums_l(nzb,62,tn) = 0.0
-          sums_l(nzb,63,tn) = 0.0
+          sums_l(nzb,58,tn) = 0.0_wp
+          sums_l(nzb,59,tn) = 0.0_wp
+          sums_l(nzb,60,tn) = 0.0_wp
+          sums_l(nzb,61,tn) = 0.0_wp
+          sums_l(nzb,62,tn) = 0.0_wp
+          sums_l(nzb,63,tn) = 0.0_wp
 
        ENDIF
@@ -1086,5 +1088,5 @@
        hom(:,1,37,sr) = sums(:,37)     ! w*e*
        hom(:,1,38,sr) = sums(:,38)     ! w*3
-       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20 )**1.5   ! Sw
+       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
        hom(:,1,40,sr) = sums(:,40)     ! p
        hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
@@ -1115,5 +1117,5 @@
        hom(:,1,70,sr) = sums(:,70)     ! q*2
        hom(:,1,71,sr) = sums(:,71)     ! prho
-       hom(:,1,72,sr) = hyp * 1E-4     ! hyp in dbar
+       hom(:,1,72,sr) = hyp * 1E-4_wp  ! hyp in dbar
        hom(:,1,73,sr) = sums(:,73)     ! nr
        hom(:,1,74,sr) = sums(:,74)     ! qr
@@ -1143,18 +1145,18 @@
 !--    is less than 1.5 times the height where the heat flux becomes negative
 !--    (positive) for the first time.
-       z_i(1) = 0.0
+       z_i(1) = 0.0_wp
        first = .TRUE.
 
        IF ( ocean )  THEN
           DO  k = nzt, nzb+1, -1
-             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
-                .AND. abs(hom(k,1,18,sr)) > 1.0E-8)  THEN
+             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
+                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp)  THEN
                 first = .FALSE.
                 height = zw(k)
              ENDIF
-             IF ( hom(k,1,18,sr) < 0.0  .AND. &
-                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
+             IF ( hom(k,1,18,sr) < 0.0_wp  .AND.                               &
+                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND.                        &
                   hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
-                IF ( zw(k) < 1.5 * height )  THEN
+                IF ( zw(k) < 1.5_wp * height )  THEN
                    z_i(1) = zw(k)
                 ELSE
@@ -1166,13 +1168,13 @@
        ELSE
           DO  k = nzb, nzt-1
-             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
-                .AND. abs(hom(k,1,18,sr)) > 1.0E-8 )  THEN
+             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
+                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
                 first = .FALSE.
                 height = zw(k)
              ENDIF
-             IF ( hom(k,1,18,sr) < 0.0  .AND. & 
-                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
+             IF ( hom(k,1,18,sr) < 0.0_wp  .AND.                               & 
+                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND.                        &
                   hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
-                IF ( zw(k) < 1.5 * height )  THEN
+                IF ( zw(k) < 1.5_wp * height )  THEN
                    z_i(1) = zw(k)
                 ELSE
@@ -1192,5 +1194,5 @@
 !--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
 !--             ocean case! 
-       z_i(2) = 0.0
+       z_i(2) = 0.0_wp
        DO  k = nzb+1, nzt+1
           dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
@@ -1225,14 +1227,14 @@
 !--    the characteristic convective boundary layer temperature.
 !--    The horizontal average at nzb+1 is input for the average temperature.
-       IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
-           .AND.  z_i(1) /= 0.0 )  THEN
-          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
-                                       hom(nzb,1,18,sr) *      &
-                                       ABS( z_i(1) ) )**0.333333333
+       IF ( hom(nzb,1,18,sr) > 0.0_wp .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8_wp  &
+           .AND.  z_i(1) /= 0.0_wp )  THEN
+          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) *                 &
+                                       hom(nzb,1,18,sr) *                      &
+                                       ABS( z_i(1) ) )**0.333333333_wp
 !--       so far this only works if Prandtl layer is used
           hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
        ELSE
-          hom(nzb+8,1,pr_palm,sr)  = 0.0
-          hom(nzb+11,1,pr_palm,sr) = 0.0
+          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
+          hom(nzb+11,1,pr_palm,sr) = 0.0_wp
        ENDIF
 
@@ -1261,9 +1263,9 @@
        ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
 
-       IF ( ts_value(5,sr) /= 0.0 )  THEN
+       IF ( ts_value(5,sr) /= 0.0_wp )  THEN
           ts_value(22,sr) = ts_value(4,sr)**2 / &
                             ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
        ELSE
-          ts_value(22,sr) = 10000.0
+          ts_value(22,sr) = 10000.0_wp
        ENDIF
 
@@ -1278,5 +1280,5 @@
 !-- If required, sum up horizontal averages for subsequent time averaging
     IF ( do_sum )  THEN
-       IF ( average_count_pr == 0 )  hom_sum = 0.0
+       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
        hom_sum = hom_sum + hom(:,1,:,:)
        average_count_pr = average_count_pr + 1
@@ -1389,5 +1391,5 @@
 !
 !--    Initialize (local) summation array
-       sums_l = 0.0
+       sums_l = 0.0_wp
 
 !
@@ -1421,12 +1423,12 @@
              sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
              sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
-             sums_l(:,34,i) = sums_l(:,34,i) + 0.5 *                        &
-                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +   &
+             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        &
+                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
                                 sums_ws2_ws_l(:,i) )    ! e*
              DO  k = nzb, nzt
-                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5 * ( &
-                                                      sums_us2_ws_l(k,i) +  &
-                                                      sums_vs2_ws_l(k,i) +  &
-                                                      sums_ws2_ws_l(k,i) )
+                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5_wp * ( &
+                                                      sums_us2_ws_l(k,i) +     &
+                                                      sums_vs2_ws_l(k,i) +     &
+                                                      sums_ws2_ws_l(k,i)     )
              ENDDO
           ENDDO
@@ -1593,30 +1595,30 @@
        !$acc update host( sums_l )
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
+       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL,  &
                            MPI_SUM, comm2d, ierr )
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
+       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL,  &
                            MPI_SUM, comm2d, ierr )
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
+       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL,  &
                            MPI_SUM, comm2d, ierr )
        IF ( ocean )  THEN
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, &
+          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb,       &
                               MPI_REAL, MPI_SUM, comm2d, ierr )
        ENDIF
        IF ( humidity ) THEN
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
+          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb,       &
                               MPI_REAL, MPI_SUM, comm2d, ierr )
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
+          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
                               MPI_REAL, MPI_SUM, comm2d, ierr )
           IF ( cloud_physics ) THEN
              IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
+             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb,    &
                                  MPI_REAL, MPI_SUM, comm2d, ierr )
              IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
+             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb,    &
                                  MPI_REAL, MPI_SUM, comm2d, ierr )
           ENDIF
@@ -1625,5 +1627,5 @@
        IF ( passive_scalar )  THEN
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
+          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
                               MPI_REAL, MPI_SUM, comm2d, ierr )
        ENDIF
@@ -1781,5 +1783,6 @@
           DO  j =  nys, nyn
              DO  k = nzb, nzt+1
-                s1 = s1 + 0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 + w(k,j,i)**2 ) * &
+                s1 = s1 + 0.5_wp *                                             & 
+                          ( u(k,j,i)**2 + v(k,j,i)**2 + w(k,j,i)**2 ) *        &
                           rmask(j,i,sr) * rflags_invers(j,i,k+1)
              ENDDO
@@ -1842,5 +1845,5 @@
 !
 !--                Perturbation energy
-                   s4 = s4 + 0.5 * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) * &
+                   s4 = s4 + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) *   &
                              rflags_invers(j,i,k+1)
                 ENDDO
@@ -1889,5 +1892,5 @@
 !
 !--             Momentum flux w"u"
-                s1 = s1 - 0.25 * (                   &
+                s1 = s1 - 0.25_wp * (                                          &
                                km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
                                                            ) * (               &
@@ -1898,5 +1901,5 @@
 !
 !--             Momentum flux w"v"
-                s2 = s2 - 0.25 * (                   &
+                s2 = s2 - 0.25_wp * (                                          &
                                km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
                                                            ) * (               &
@@ -1907,7 +1910,8 @@
 !
 !--             Heat flux w"pt"
-                s3 = s3 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
-                              * ( pt(k+1,j,i) - pt(k,j,i) )   &
-                              * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                s3 = s3 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                 &
+                                 * ( pt(k+1,j,i) - pt(k,j,i) )                 &
+                                 * ddzu(k+1) * rmask(j,i,sr)                   &
+                                 * rflags_invers(j,i,k+1)
              ENDDO
           ENDDO
@@ -1926,7 +1930,8 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
-                                 * ( sa(k+1,j,i) - sa(k,j,i) )   &
-                                 * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
+                                    * ( sa(k+1,j,i) - sa(k,j,i) )              &
+                                    * ddzu(k+1) * rmask(j,i,sr)                & 
+                                    * rflags_invers(j,i,k+1)
                 ENDDO
              ENDDO
@@ -1945,10 +1950,12 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
-                                 * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
-                                 * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
-                   s2 = s2 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
-                                 * ( q(k+1,j,i) - q(k,j,i) )     &
-                                 * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
+                                    * ( vpt(k+1,j,i) - vpt(k,j,i) )            &
+                                    * ddzu(k+1) * rmask(j,i,sr) 
+                                    * rflags_invers(j,i,k+1)
+                   s2 = s2 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
+                                    * ( q(k+1,j,i) - q(k,j,i) )                &
+                                    * ddzu(k+1) * rmask(j,i,sr)                &
+                                    * rflags_invers(j,i,k+1)
                 ENDDO
              ENDDO
@@ -1965,8 +1972,9 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )    &
-                                    * ( ( q(k+1,j,i) - ql(k+1,j,i) ) &
-                                      - ( q(k,j,i) - ql(k,j,i) ) )   &
-                                    * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                      s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )           &
+                                       *  ( ( q(k+1,j,i) - ql(k+1,j,i) )       &
+                                          - ( q(k,j,i) - ql(k,j,i) ) )         &
+                                       * ddzu(k+1) * rmask(j,i,sr)             & 
+                                       * rflags_invers(j,i,k+1)
                    ENDDO
                 ENDDO
@@ -1987,7 +1995,8 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
-                                 * ( q(k+1,j,i) - q(k,j,i) )     &
-                                 * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
+                                    * ( q(k+1,j,i) - q(k,j,i) )                &
+                                    * ddzu(k+1) * rmask(j,i,sr)                &
+                                    * rflags_invers(j,i,k+1)
                 ENDDO
              ENDDO
@@ -2010,6 +2019,6 @@
                 s2 = s2 + vsws(j,i) * rmask(j,i,sr)     ! w"v"
                 s3 = s3 + shf(j,i)  * rmask(j,i,sr)     ! w"pt"
-                s4 = s4 + 0.0 * rmask(j,i,sr)           ! u"pt"
-                s5 = s5 + 0.0 * rmask(j,i,sr)           ! v"pt"
+                s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
+                s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
              ENDDO
           ENDDO
@@ -2044,6 +2053,6 @@
                 DO  j =  nys, nyn
                    s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
-                   s2 = s2 + ( ( 1.0 + 0.61 * q(nzb,j,i) ) * shf(j,i) &
-                               + 0.61 * pt(nzb,j,i) * qsws(j,i) )
+                   s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * shf(j,i)    &
+                               + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
                 ENDDO
              ENDDO
@@ -2059,6 +2068,7 @@
                 DO  i = nxl, nxr
                    DO  j =  nys, nyn
-                      s1 = s1 + ( ( 1.0 + 0.61 * q(nzb,j,i) - ql(nzb,j,i) ) * &
-                                  shf(j,i) + 0.61 * pt(nzb,j,i) * qsws(j,i) )
+                      s1 = s1 + ( ( 1.0_wp +                                   &
+                                    0.61_wp * q(nzb,j,i) - ql(nzb,j,i) ) *     &
+                                  shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
                    ENDDO
                 ENDDO
@@ -2116,6 +2126,6 @@
                 s2 = s2 + vswst(j,i) * rmask(j,i,sr)    ! w"v"
                 s3 = s3 + tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
-                s4 = s4 + 0.0 * rmask(j,i,sr)           ! u"pt"
-                s5 = s5 + 0.0 * rmask(j,i,sr)           ! v"pt"
+                s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
+                s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
              ENDDO
           ENDDO
@@ -2150,6 +2160,6 @@
                 DO  j =  nys, nyn
                    s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
-                   s2 = s2 + ( ( 1.0 + 0.61 * q(nzt,j,i) ) * tswst(j,i) + &
-                               0.61 * pt(nzt,j,i) * qswst(j,i) )
+                   s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * tswst(j,i) +&
+                                 0.61_wp * pt(nzt,j,i) * qswst(j,i) )
                 ENDDO
              ENDDO
@@ -2165,6 +2175,8 @@
                 DO  i = nxl, nxr
                    DO  j =  nys, nyn
-                      s1 = s1 + ( ( 1.0 + 0.61 * q(nzt,j,i) - ql(nzt,j,i) ) * &
-                                  tswst(j,i) + 0.61 * pt(nzt,j,i) * qswst(j,i) )
+                      s1 = s1 + ( ( 1.0_wp +                                   &
+                                    0.61_wp * q(nzt,j,i) - ql(nzt,j,i) ) *     &
+                                  tswst(j,i) +                                 &
+                                  0.61_wp * pt(nzt,j,i) * qswst(j,i) )
                    ENDDO
                 ENDDO
@@ -2220,10 +2232,10 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
-                              u(k+1,j,i) - hom(k+1,1,1,sr) )
-                vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
-                              v(k+1,j,i) - hom(k+1,1,2,sr) )
-                pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
-                              pt(k+1,j,i) - hom(k+1,1,4,sr) )
+                ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) + &
+                                 u(k+1,j,i) - hom(k+1,1,1,sr) )
+                vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) + &
+                                 v(k+1,j,i) - hom(k+1,1,2,sr) )
+                pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) + &
+                                 pt(k+1,j,i) - hom(k+1,1,4,sr) )
 !
 !--             Higher moments
@@ -2232,5 +2244,5 @@
 !
 !--             Energy flux w*e* (has to be adjusted?)
-                s3 = s3 + w(k,j,i) * 0.5 * ( ust**2 + vst**2 + w(k,j,i)**2 ) &
+                s3 = s3 + w(k,j,i) * 0.5_wp * ( ust**2 + vst**2 + w(k,j,i)**2 )&
                                    * rmask(j,i,sr) * rflags_invers(j,i,k+1)
              ENDDO
@@ -2254,7 +2266,8 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      s1 = s1 + 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) +   &
-                                        sa(k+1,j,i) - hom(k+1,1,23,sr) ) &
-                                    * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                      s1 = s1 + 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +      &
+                                           sa(k+1,j,i) - hom(k+1,1,23,sr) )    &
+                                       * w(k,j,i) * rmask(j,i,sr)              & 
+                                       * rflags_invers(j,i,k+1)
                    ENDDO
                 ENDDO
@@ -2293,7 +2306,7 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      s1 = s1 + 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +     &
-                                        vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
-                                      w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                      s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +     &
+                                           vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
+                                         w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
                    ENDDO
                 ENDDO
@@ -2309,7 +2322,7 @@
                    DO  i = nxl, nxr
                       DO  j = nys, nyn
-                         s1 = s1 + 0.5 * ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) +   &
-                                           ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
-                                       * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                         s1 = s1 + 0.5_wp * ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) +   &
+                                              ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
+                                          * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
                       ENDDO
                    ENDDO
@@ -2395,7 +2408,7 @@
                    DO  i = nxl, nxr
                       DO  j = nys, nyn
-                         s1 = s1 + 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +   &
-                                           vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
-                                       * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                         s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +   &
+                                              vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
+                                          * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
                       ENDDO
                    ENDDO
@@ -2408,6 +2421,6 @@
                 !$acc parallel loop present( hom, sums_l )
                 DO  k = nzb, nzt_diff
-                   sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
-                                             0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
+                   sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
+                                                0.61_wp * hom(k,1,4,sr) * sums_l(k,49,tn)
                 ENDDO
                 !$acc end parallel loop
@@ -2427,7 +2440,8 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   s1 = s1 + 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) +   &
-                                     q(k+1,j,i) - hom(k+1,1,41,sr) ) &
-                                 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +          &
+                                        q(k+1,j,i) - hom(k+1,1,41,sr) )        &
+                                    * w(k,j,i) * rmask(j,i,sr)                 &
+                                    * rflags_invers(j,i,k+1)
                 ENDDO
              ENDDO
@@ -2450,16 +2464,18 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
-                               u(k+1,j,i) - hom(k+1,1,1,sr) )
-                   vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
-                               v(k+1,j,i) - hom(k+1,1,2,sr) )
+                   ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +               &
+                                    u(k+1,j,i) - hom(k+1,1,1,sr) )
+                   vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +               &
+                                    v(k+1,j,i) - hom(k+1,1,2,sr) )
 !
 !--                Momentum flux w*u*
-                   s1 = s1 + 0.5 * ( w(k,j,i-1) + w(k,j,i) ) &
-                                 * ust * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   s1 = s1 + 0.5_wp * ( w(k,j,i-1) + w(k,j,i) )                &
+                                    * ust * rmask(j,i,sr)                      &
+                                    * rflags_invers(j,i,k+1)
 !
 !--                Momentum flux w*v*
-                   s2 = s2 + 0.5 * ( w(k,j-1,i) + w(k,j,i) ) &
-                                 * vst * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   s2 = s2 + 0.5_wp * ( w(k,j-1,i) + w(k,j,i) )                &
+                                    * vst * rmask(j,i,sr)                      & 
+                                    * rflags_invers(j,i,k+1)
                 ENDDO
              ENDDO
@@ -2481,7 +2497,8 @@
 !
 !--                Vertical heat flux
-                   s1 = s1 + 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
-                                     pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
-                                 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   s1 = s1 + 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +          &
+                                        pt(k+1,j,i) - hom(k+1,1,4,sr) )        &
+                                    * w(k,j,i) * rmask(j,i,sr)                 & 
+                                    * rflags_invers(j,i,k+1)
                 ENDDO
              ENDDO
@@ -2497,7 +2514,8 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      s1 = s1 + 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) +   &
-                                        q(k+1,j,i) - hom(k+1,1,41,sr) ) &
-                                    * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                      s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +       &
+                                           q(k+1,j,i) - hom(k+1,1,41,sr) )     &
+                                       * w(k,j,i) * rmask(j,i,sr)              &
+                                       * rflags_invers(j,i,k+1)
                    ENDDO
                 ENDDO
@@ -2525,8 +2543,8 @@
 !--    First calculate the products, then the divergence.
 !--    Calculation is time consuming. Do it only, if profiles shall be plotted.
-       IF ( hom(nzb+1,2,55,0) /= 0.0  .OR.  hom(nzb+1,2,68,0) /= 0.0 )  THEN
+       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )  THEN
 
           STOP '+++ openACC porting for vertical flux div of resolved scale TKE in flow_statistics is still missing'
-          sums_ll = 0.0  ! local array
+          sums_ll = 0.0_wp  ! local array
 
           !$OMP DO
@@ -2535,14 +2553,14 @@
                 DO  k = nzb_s_inner(j,i)+1, nzt
 
-                   sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * (        &
-                  ( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)   &
-                           - 0.5 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
-                           ) )**2                                          &
-                + ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)   &
-                           - 0.5 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
-                           ) )**2                                          &
-                   + w(k,j,i)**2                                  )
-
-                   sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
+                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
+                  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)    &
+                              - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )   &
+                              ) )**2                                           &
+                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)    &
+                              - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )   &
+                              ) )**2                                           &
+                + w(k,j,i)**2                                        )
+
+                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
                                                * ( p(k,j,i) + p(k+1,j,i) )
 
@@ -2550,8 +2568,8 @@
              ENDDO
           ENDDO
-          sums_ll(0,1)     = 0.0    ! because w is zero at the bottom
-          sums_ll(nzt+1,1) = 0.0
-          sums_ll(0,2)     = 0.0
-          sums_ll(nzt+1,2) = 0.0
+          sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
+          sums_ll(nzt+1,1) = 0.0_wp
+          sums_ll(0,2)     = 0.0_wp
+          sums_ll(nzt+1,2) = 0.0_wp
 
           DO  k = nzb+1, nzt
@@ -2562,5 +2580,5 @@
           sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
           sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
-          sums_l(nzb,68,tn) = 0.0    ! because w* = 0 at nzb
+          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
 
        ENDIF
@@ -2568,5 +2586,5 @@
 !
 !--    Divergence of vertical flux of SGS TKE and the flux itself (69)
-       IF ( hom(nzb+1,2,57,0) /= 0.0  .OR.  hom(nzb+1,2,69,0) /= 0.0 )  THEN
+       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )  THEN
 
           STOP '+++ openACC porting for vertical flux div of SGS TKE in flow_statistics is still missing'
@@ -2576,12 +2594,12 @@
                 DO  k = nzb_s_inner(j,i)+1, nzt
 
-                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * (                 &
+                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
                    (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
                  - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
-                                                             ) * ddzw(k)
-
-                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * (                 &
+                                                                ) * ddzw(k)
+
+                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
                    (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
-                                                              )
+                                                                )
 
                 ENDDO
@@ -2596,5 +2614,5 @@
 !--    Horizontal heat fluxes (subgrid, resolved, total).
 !--    Do it only, if profiles shall be plotted.
-       IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
+       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
 
           STOP '+++ openACC porting for horizontal flux calculation in flow_statistics is still missing'
@@ -2605,9 +2623,9 @@
 !
 !--                Subgrid horizontal heat fluxes u"pt", v"pt"
-                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 *                   &
+                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
                                                    ( kh(k,j,i) + kh(k,j,i-1) ) &
                                                  * ( pt(k,j,i-1) - pt(k,j,i) ) &
                                                  * ddx * rmask(j,i,sr)
-                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 *                   &
+                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
                                                    ( kh(k,j,i) + kh(k,j-1,i) ) &
                                                  * ( pt(k,j-1,i) - pt(k,j,i) ) &
@@ -2616,13 +2634,13 @@
 !--                Resolved horizontal heat fluxes u*pt*, v*pt*
                    sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
-                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
-                                       * 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
-                                                 pt(k,j,i)   - hom(k,1,4,sr) )
-                   pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
-                                 pt(k,j,i)   - hom(k,1,4,sr) )
+                                     ( u(k,j,i) - hom(k,1,1,sr) ) * 0.5_wp *   &
+                                     ( pt(k,j,i-1) - hom(k,1,4,sr) +           &
+                                       pt(k,j,i)   - hom(k,1,4,sr) )
+                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
+                                    pt(k,j,i)   - hom(k,1,4,sr) )
                    sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
-                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
-                                       * 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
-                                                 pt(k,j,i)   - hom(k,1,4,sr) )
+                                     ( v(k,j,i) - hom(k,1,2,sr) ) * 0.5_wp *   &
+                                     ( pt(k,j-1,i) - hom(k,1,4,sr) +           &
+                                       pt(k,j,i)   - hom(k,1,4,sr) )
                 ENDDO
              ENDDO
@@ -2630,10 +2648,10 @@
 !
 !--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
-          sums_l(nzb,58,tn) = 0.0
-          sums_l(nzb,59,tn) = 0.0
-          sums_l(nzb,60,tn) = 0.0
-          sums_l(nzb,61,tn) = 0.0
-          sums_l(nzb,62,tn) = 0.0
-          sums_l(nzb,63,tn) = 0.0
+          sums_l(nzb,58,tn) = 0.0_wp
+          sums_l(nzb,59,tn) = 0.0_wp
+          sums_l(nzb,60,tn) = 0.0_wp
+          sums_l(nzb,61,tn) = 0.0_wp
+          sums_l(nzb,62,tn) = 0.0_wp
+          sums_l(nzb,63,tn) = 0.0_wp
 
        ENDIF
@@ -2752,5 +2770,5 @@
        hom(:,1,37,sr) = sums(:,37)     ! w*e*
        hom(:,1,38,sr) = sums(:,38)     ! w*3
-       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20 )**1.5   ! Sw
+       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
        hom(:,1,40,sr) = sums(:,40)     ! p
        hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
@@ -2781,5 +2799,5 @@
        hom(:,1,70,sr) = sums(:,70)     ! q*2
        hom(:,1,71,sr) = sums(:,71)     ! prho
-       hom(:,1,72,sr) = hyp * 1E-4     ! hyp in dbar
+       hom(:,1,72,sr) = hyp * 1E-4_wp     ! hyp in dbar
        hom(:,1,73,sr) = sums(:,73)     ! nr
        hom(:,1,74,sr) = sums(:,74)     ! qr
@@ -2809,18 +2827,18 @@
 !--    is less than 1.5 times the height where the heat flux becomes negative
 !--    (positive) for the first time.
-       z_i(1) = 0.0
+       z_i(1) = 0.0_wp
        first = .TRUE.
 
        IF ( ocean )  THEN
           DO  k = nzt, nzb+1, -1
-             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
-                .AND. abs(hom(k,1,18,sr)) > 1.0E-8)  THEN
+             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
+                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
                 first = .FALSE.
                 height = zw(k)
              ENDIF
-             IF ( hom(k,1,18,sr) < 0.0  .AND. &
-                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
+             IF ( hom(k,1,18,sr) < 0.0_wp  .AND.                               &
+                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND.                        &
                   hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
-                IF ( zw(k) < 1.5 * height )  THEN
+                IF ( zw(k) < 1.5_wp * height )  THEN
                    z_i(1) = zw(k)
                 ELSE
@@ -2832,13 +2850,13 @@
        ELSE
           DO  k = nzb, nzt-1
-             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
-                .AND. abs(hom(k,1,18,sr)) > 1.0E-8 )  THEN
+             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
+                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
                 first = .FALSE.
                 height = zw(k)
              ENDIF
              IF ( hom(k,1,18,sr) < 0.0  .AND. &
-                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
+                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND. &
                   hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
-                IF ( zw(k) < 1.5 * height )  THEN
+                IF ( zw(k) < 1.5_wp * height )  THEN
                    z_i(1) = zw(k)
                 ELSE
@@ -2858,5 +2876,5 @@
 !--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
 !--             ocean case!
-       z_i(2) = 0.0
+       z_i(2) = 0.0_wp
        DO  k = nzb+1, nzt+1
           dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
@@ -2891,14 +2909,14 @@
 !--    the characteristic convective boundary layer temperature.
 !--    The horizontal average at nzb+1 is input for the average temperature.
-       IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
-           .AND.  z_i(1) /= 0.0 )  THEN
-          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
-                                       hom(nzb,1,18,sr) *      &
-                                       ABS( z_i(1) ) )**0.333333333
+       IF ( hom(nzb,1,18,sr) > 0.0_wp .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8_wp  &
+           .AND.  z_i(1) /= 0.0_wp )  THEN
+          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) *                 &
+                                       hom(nzb,1,18,sr) *                      &
+                                       ABS( z_i(1) ) )**0.333333333_wp
 !--       so far this only works if Prandtl layer is used
           hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
        ELSE
-          hom(nzb+8,1,pr_palm,sr)  = 0.0
-          hom(nzb+11,1,pr_palm,sr) = 0.0
+          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
+          hom(nzb+11,1,pr_palm,sr) = 0.0_wp
        ENDIF
 
@@ -2927,9 +2945,9 @@
        ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
 
-       IF ( ts_value(5,sr) /= 0.0 )  THEN
-          ts_value(22,sr) = ts_value(4,sr)**2 / &
+       IF ( ts_value(5,sr) /= 0.0_wp )  THEN
+          ts_value(22,sr) = ts_value(4,sr)**2_wp / &
                             ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
        ELSE
-          ts_value(22,sr) = 10000.0
+          ts_value(22,sr) = 10000.0_wp
        ENDIF
 
@@ -2946,5 +2964,5 @@
 !-- If required, sum up horizontal averages for subsequent time averaging
     IF ( do_sum )  THEN
-       IF ( average_count_pr == 0 )  hom_sum = 0.0
+       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
        hom_sum = hom_sum + hom(:,1,:,:)
        average_count_pr = average_count_pr + 1
