Ignore:
Timestamp:
Apr 8, 2014 3:21:23 PM (7 years ago)
Author:
heinze
Message:

REAL constants provided with KIND-attribute

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/flow_statistics.f90

    r1323 r1353  
    2121! Current revisions:
    2222! -----------------
    23 !
     23! REAL constants provided with KIND-attribute
    2424!
    2525! Former revisions:
     
    170170!
    171171!--    Initialize (local) summation array
    172        sums_l = 0.0
     172       sums_l = 0.0_wp
    173173
    174174!
     
    202202             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
    203203             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
    204              sums_l(:,34,i) = sums_l(:,34,i) + 0.5 *                           &
     204             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        &
    205205                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
    206206                                sums_ws2_ws_l(:,i) )    ! e*
    207207             DO  k = nzb, nzt
    208                 sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5 * (    &
     208                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5_wp * ( &
    209209                                                      sums_us2_ws_l(k,i) +     &
    210210                                                      sums_vs2_ws_l(k,i) +     &
     
    446446       DO  i = nxl, nxr
    447447          DO  j =  nys, nyn
    448              sums_l_etot = 0.0
     448             sums_l_etot = 0.0_wp
    449449             DO  k = nzb_s_inner(j,i), nzt+1
    450450!
     
    470470
    471471                sums_l_etot  = sums_l_etot + &
    472                                         0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 +    &
     472                                        0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 + &
    473473                                        w(k,j,i)**2 ) * rmask(j,i,sr)
    474474             ENDDO
     
    503503          DO  i = nxl, nxr
    504504             DO  j =  nys, nyn
    505                 sums_l_eper = 0.0
     505                sums_l_eper = 0.0_wp
    506506                DO  k = nzb_s_inner(j,i), nzt+1
    507507                   u2   = u(k,j,i)**2
     
    517517!--             Perturbation energy
    518518
    519                    sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 *       &
     519                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5_wp *                &
    520520                                  ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
    521                    sums_l_eper  = sums_l_eper + &
    522                                         0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
    523 
    524                 ENDDO
    525                 sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)   &
     521                   sums_l_eper     = sums_l_eper +                            &
     522                                     0.5_wp * ( ust2+vst2+w2 ) * rmask(j,i,sr)
     523
     524                ENDDO
     525                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)            &
    526526                     + sums_l_eper
    527527             ENDDO
     
    547547!
    548548!--             Momentum flux w"u"
    549                 sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25 * (                   &
     549                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25_wp * (                &
    550550                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
    551551                                                           ) * (               &
     
    555555!
    556556!--             Momentum flux w"v"
    557                 sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25 * (                   &
     557                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25_wp * (                &
    558558                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
    559559                                                           ) * (               &
     
    564564!--             Heat flux w"pt"
    565565                sums_l(k,16,tn) = sums_l(k,16,tn)                              &
    566                                          - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
     566                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    567567                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
    568568                                               * ddzu(k+1) * rmask(j,i,sr)
     
    573573                IF ( ocean )  THEN
    574574                   sums_l(k,65,tn) = sums_l(k,65,tn)                           &
    575                                          - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
     575                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    576576                                               * ( sa(k+1,j,i) - sa(k,j,i) )   &
    577577                                               * ddzu(k+1) * rmask(j,i,sr)
     
    582582                IF ( humidity ) THEN
    583583                   sums_l(k,45,tn) = sums_l(k,45,tn)                           &
    584                                          - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
     584                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    585585                                               * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
    586586                                               * ddzu(k+1) * rmask(j,i,sr)
    587587                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
    588                                          - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
     588                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    589589                                               * ( q(k+1,j,i) - q(k,j,i) )     &
    590590                                               * ddzu(k+1) * rmask(j,i,sr)
     
    592592                   IF ( cloud_physics ) THEN
    593593                      sums_l(k,51,tn) = sums_l(k,51,tn)                        &
    594                                          - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
     594                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    595595                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
    596596                                                - ( q(k,j,i) - ql(k,j,i) ) )   &
     
    603603                IF ( passive_scalar )  THEN
    604604                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
    605                                          - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
     605                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    606606                                               * ( q(k+1,j,i) - q(k,j,i) )     &
    607607                                               * ddzu(k+1) * rmask(j,i,sr)
     
    620620                                    shf(j,i)  * rmask(j,i,sr)     ! w"pt"
    621621                sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
    622                                     0.0 * rmask(j,i,sr)           ! u"pt"
     622                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
    623623                sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
    624                                     0.0 * rmask(j,i,sr)           ! v"pt"
     624                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
    625625                IF ( ocean )  THEN
    626626                   sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
     
    628628                ENDIF
    629629                IF ( humidity )  THEN
    630                    sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
     630                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
    631631                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
    632                    sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (           &
    633                                        ( 1.0 + 0.61 * q(nzb,j,i) ) *   &
    634                                        shf(j,i) + 0.61 * pt(nzb,j,i) * &
     632                   sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
     633                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
     634                                       shf(j,i) + 0.61_wp * pt(nzb,j,i) *      &
    635635                                                  qsws(j,i) )
    636636                   IF ( cloud_droplets )  THEN
    637                       sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (        &
    638                                          ( 1.0 + 0.61 * q(nzb,j,i) -   &
    639                                            ql(nzb,j,i) ) * shf(j,i) +  &
    640                                            0.61 * pt(nzb,j,i) * qsws(j,i) )
     637                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
     638                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
     639                                           ql(nzb,j,i) ) * shf(j,i) +          &
     640                                           0.61_wp * pt(nzb,j,i) * qsws(j,i) )
    641641                   ENDIF
    642642                   IF ( cloud_physics )  THEN
     
    663663                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
    664664                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
    665                                     0.0 * rmask(j,i,sr)           ! u"pt"
     665                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
    666666                sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
    667                                     0.0 * rmask(j,i,sr)           ! v"pt"
     667                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
    668668
    669669                IF ( ocean )  THEN
     
    672672                ENDIF
    673673                IF ( humidity )  THEN
    674                    sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
     674                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
    675675                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
    676                    sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (             &
    677                                        ( 1.0 + 0.61 * q(nzt,j,i) ) *     &
    678                                        tswst(j,i) + 0.61 * pt(nzt,j,i) * &
    679                                                   qswst(j,i) )
     676                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                   &
     677                                       ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *     &
     678                                       tswst(j,i) + 0.61_wp * pt(nzt,j,i) *    &
     679                                                             qswst(j,i) )
    680680                   IF ( cloud_droplets )  THEN
    681                       sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
    682                                           ( 1.0 + 0.61 * q(nzt,j,i) -     &
    683                                             ql(nzt,j,i) ) * tswst(j,i) +  &
    684                                             0.61 * pt(nzt,j,i) * qswst(j,i) )
     681                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                &
     682                                          ( 1.0_wp + 0.61_wp * q(nzt,j,i) -    &
     683                                            ql(nzt,j,i) ) * tswst(j,i) +       &
     684                                            0.61_wp * pt(nzt,j,i) * qswst(j,i) )
    685685                   ENDIF
    686686                   IF ( cloud_physics )  THEN
     
    703703!--                rearranged according to the staggered grid.
    704704             DO  k = nzb_s_inner(j,i), nzt
    705                 ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
    706                               u(k+1,j,i) - hom(k+1,1,1,sr) )
    707                 vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
    708                               v(k+1,j,i) - hom(k+1,1,2,sr) )
    709                 pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
    710                               pt(k+1,j,i) - hom(k+1,1,4,sr) )
     705                ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                  &
     706                                 u(k+1,j,i) - hom(k+1,1,1,sr) )
     707                vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                  &
     708                                 v(k+1,j,i) - hom(k+1,1,2,sr) )
     709                pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +                &
     710                                 pt(k+1,j,i) - hom(k+1,1,4,sr) )
    711711
    712712!--             Higher moments
    713                 sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 * &
     713                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 *        &
    714714                                                    rmask(j,i,sr)
    715                 sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) * &
     715                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) *        &
    716716                                                    rmask(j,i,sr)
    717717
     
    721721                IF ( ocean )  THEN
    722722                   IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
    723                       pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
    724                                  sa(k+1,j,i) - hom(k+1,1,23,sr) )
    725                       sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
     723                      pts = 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +          &
     724                                       sa(k+1,j,i) - hom(k+1,1,23,sr) )
     725                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) *     &
     726                                        rmask(j,i,sr)
     727                   ENDIF
     728                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) *            &
    726729                                                       rmask(j,i,sr)
    727                    ENDIF
    728                    sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
    729                                                        rmask(j,i,sr)
    730                    sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) * &
     730                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) *           &
    731731                                                       rmask(j,i,sr)
    732732                ENDIF
     
    737737                IF ( humidity )  THEN
    738738                   IF ( cloud_physics .OR. cloud_droplets )  THEN
    739                       pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +        &
     739                      pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +         &
    740740                                    vpt(k+1,j,i) - hom(k+1,1,44,sr) )
    741                       sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
     741                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
    742742                                                          rmask(j,i,sr)
    743743                      IF ( .NOT. cloud_droplets )  THEN
    744                          pts = 0.5 *                                           &
     744                         pts = 0.5_wp *                                        &
    745745                              ( ( q(k,j,i) - ql(k,j,i) ) -                     &
    746746                              hom(k,1,42,sr) +                                 &
     
    772772                   ELSE
    773773                      IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
    774                          pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
    775                                        vpt(k+1,j,i) - hom(k+1,1,44,sr) )
    776                          sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
     774                         pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +      &
     775                                          vpt(k+1,j,i) - hom(k+1,1,44,sr) )
     776                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
    777777                                                             rmask(j,i,sr)
    778778                      ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
    779                          sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) *  &
    780                                              sums_l(k,17,tn) +      &
    781                                           0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
     779                         sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp *                &
     780                                             hom(k,1,41,sr) ) *                &
     781                                           sums_l(k,17,tn) +                   &
     782                                           0.61_wp * hom(k,1,4,sr) *           &
     783                                           sums_l(k,49,tn)
    782784                      END IF
    783785                   END IF
     
    785787!
    786788!--             Passive scalar flux
    787                 IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca   &
     789                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca                &
    788790                     .OR. sr /= 0 ) )  THEN
    789                    pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
    790                                  q(k+1,j,i) - hom(k+1,1,41,sr) )
    791                    sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
     791                   pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +              &
     792                                    q(k+1,j,i) - hom(k+1,1,41,sr) )
     793                   sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *        &
    792794                                                       rmask(j,i,sr)
    793795                ENDIF
     
    796798!--             Energy flux w*e*
    797799!--             has to be adjusted
    798                 sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *          &
    799                                              ( ust**2 + vst**2 + w(k,j,i)**2 )&
     800                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp *        &
     801                                             ( ust**2 + vst**2 + w(k,j,i)**2 ) &
    800802                                             * rmask(j,i,sr)
    801803             ENDDO
     
    811813            DO  j = nys, nyn
    812814               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
    813                   ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
    814                               u(k+1,j,i) - hom(k+1,1,1,sr) )
    815                   vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
    816                               v(k+1,j,i) - hom(k+1,1,2,sr) )
     815                  ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                &
     816                                   u(k+1,j,i) - hom(k+1,1,1,sr) )
     817                  vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                &
     818                                   v(k+1,j,i) - hom(k+1,1,2,sr) )
    817819!
    818820!--               Momentum flux w*u*
    819                   sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                   &
    820                                                     ( w(k,j,i-1) + w(k,j,i) ) &
     821                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                 &
     822                                                    ( w(k,j,i-1) + w(k,j,i) )  &
    821823                                                    * ust * rmask(j,i,sr)
    822824!
    823825!--               Momentum flux w*v*
    824                   sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 *                   &
    825                                                     ( w(k,j-1,i) + w(k,j,i) ) &
     826                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                 &
     827                                                    ( w(k,j-1,i) + w(k,j,i) )  &
    826828                                                    * vst * rmask(j,i,sr)
    827829               ENDDO
     
    837839!
    838840!--               Vertical heat flux
    839                   sums_l(k,17,tn) = sums_l(k,17,tn) +  0.5 * &
    840                            ( pt(k,j,i)   - hom(k,1,4,sr) + &
    841                            pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
     841                  sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp *                &
     842                           ( pt(k,j,i)   - hom(k,1,4,sr) +                     &
     843                             pt(k+1,j,i) - hom(k+1,1,4,sr) )                  &
    842844                           * w(k,j,i) * rmask(j,i,sr)
    843845                  IF ( humidity )  THEN
    844                      pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
    845                                 q(k+1,j,i) - hom(k+1,1,41,sr) )
    846                      sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
    847                                                       rmask(j,i,sr)
     846                     pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +            &
     847                                      q(k+1,j,i) - hom(k+1,1,41,sr) )
     848                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *      &
     849                                       rmask(j,i,sr)
    848850                  ENDIF
    849851               ENDDO
     
    865867!--    First calculate the products, then the divergence.
    866868!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
    867        IF ( hom(nzb+1,2,55,0) /= 0.0  .OR.  hom(nzb+1,2,68,0) /= 0.0 )  THEN
    868 
    869           sums_ll = 0.0  ! local array
     869       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )  THEN
     870
     871          sums_ll = 0.0_wp  ! local array
    870872
    871873          !$OMP DO
     
    874876                DO  k = nzb_s_inner(j,i)+1, nzt
    875877
    876                    sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * (        &
    877                   ( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)   &
    878                            - 0.5 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
    879                            ) )**2                                          &
    880                 + ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)   &
    881                            - 0.5 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
    882                            ) )**2                                          &
    883                    + w(k,j,i)**2                                  )
    884 
    885                    sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
     878                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
     879                  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)    &
     880                            - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
     881                              ) )**2                                           &
     882                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)    &
     883                            - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
     884                              ) )**2                                           &
     885                + w(k,j,i)**2                                        )
     886
     887                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)            &
    886888                                               * ( p(k,j,i) + p(k+1,j,i) )
    887889
     
    889891             ENDDO
    890892          ENDDO
    891           sums_ll(0,1)     = 0.0    ! because w is zero at the bottom
    892           sums_ll(nzt+1,1) = 0.0
    893           sums_ll(0,2)     = 0.0
    894           sums_ll(nzt+1,2) = 0.0
     893          sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
     894          sums_ll(nzt+1,1) = 0.0_wp
     895          sums_ll(0,2)     = 0.0_wp
     896          sums_ll(nzt+1,2) = 0.0_wp
    895897
    896898          DO  k = nzb+1, nzt
     
    901903          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
    902904          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
    903           sums_l(nzb,68,tn) = 0.0    ! because w* = 0 at nzb
     905          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
    904906
    905907       ENDIF
     
    907909!
    908910!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
    909        IF ( hom(nzb+1,2,57,0) /= 0.0  .OR.  hom(nzb+1,2,69,0) /= 0.0 )  THEN
     911       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )  THEN
    910912
    911913          !$OMP DO
     
    914916                DO  k = nzb_s_inner(j,i)+1, nzt
    915917
    916                    sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * (                 &
     918                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
    917919                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
    918920                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
    919                                                              ) * ddzw(k)
    920 
    921                    sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * (                 &
     921                                                                ) * ddzw(k)
     922
     923                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
    922924                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
    923                                                               )
     925                                                                )
    924926
    925927                ENDDO
     
    934936!--    Horizontal heat fluxes (subgrid, resolved, total).
    935937!--    Do it only, if profiles shall be plotted.
    936        IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
     938       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
    937939
    938940          !$OMP DO
     
    942944!
    943945!--                Subgrid horizontal heat fluxes u"pt", v"pt"
    944                    sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 *                   &
     946                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
    945947                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
    946948                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
    947949                                                 * ddx * rmask(j,i,sr)
    948                    sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 *                   &
     950                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
    949951                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
    950952                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
     
    954956                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
    955957                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
    956                                        * 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
     958                                    * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
    957959                                                 pt(k,j,i)   - hom(k,1,4,sr) )
    958                    pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
    959                                  pt(k,j,i)   - hom(k,1,4,sr) )
     960                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
     961                                    pt(k,j,i)   - hom(k,1,4,sr) )
    960962                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
    961963                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
    962                                        * 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
     964                                    * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
    963965                                                 pt(k,j,i)   - hom(k,1,4,sr) )
    964966                ENDDO
     
    967969!
    968970!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
    969           sums_l(nzb,58,tn) = 0.0
    970           sums_l(nzb,59,tn) = 0.0
    971           sums_l(nzb,60,tn) = 0.0
    972           sums_l(nzb,61,tn) = 0.0
    973           sums_l(nzb,62,tn) = 0.0
    974           sums_l(nzb,63,tn) = 0.0
     971          sums_l(nzb,58,tn) = 0.0_wp
     972          sums_l(nzb,59,tn) = 0.0_wp
     973          sums_l(nzb,60,tn) = 0.0_wp
     974          sums_l(nzb,61,tn) = 0.0_wp
     975          sums_l(nzb,62,tn) = 0.0_wp
     976          sums_l(nzb,63,tn) = 0.0_wp
    975977
    976978       ENDIF
     
    10861088       hom(:,1,37,sr) = sums(:,37)     ! w*e*
    10871089       hom(:,1,38,sr) = sums(:,38)     ! w*3
    1088        hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20 )**1.5   ! Sw
     1090       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
    10891091       hom(:,1,40,sr) = sums(:,40)     ! p
    10901092       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
     
    11151117       hom(:,1,70,sr) = sums(:,70)     ! q*2
    11161118       hom(:,1,71,sr) = sums(:,71)     ! prho
    1117        hom(:,1,72,sr) = hyp * 1E-4     ! hyp in dbar
     1119       hom(:,1,72,sr) = hyp * 1E-4_wp  ! hyp in dbar
    11181120       hom(:,1,73,sr) = sums(:,73)     ! nr
    11191121       hom(:,1,74,sr) = sums(:,74)     ! qr
     
    11431145!--    is less than 1.5 times the height where the heat flux becomes negative
    11441146!--    (positive) for the first time.
    1145        z_i(1) = 0.0
     1147       z_i(1) = 0.0_wp
    11461148       first = .TRUE.
    11471149
    11481150       IF ( ocean )  THEN
    11491151          DO  k = nzt, nzb+1, -1
    1150              IF ( first .AND. hom(k,1,18,sr) < 0.0 &
    1151                 .AND. abs(hom(k,1,18,sr)) > 1.0E-8)  THEN
     1152             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
     1153                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp)  THEN
    11521154                first = .FALSE.
    11531155                height = zw(k)
    11541156             ENDIF
    1155              IF ( hom(k,1,18,sr) < 0.0  .AND. &
    1156                   abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
     1157             IF ( hom(k,1,18,sr) < 0.0_wp  .AND.                              &
     1158                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND.                        &
    11571159                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
    1158                 IF ( zw(k) < 1.5 * height )  THEN
     1160                IF ( zw(k) < 1.5_wp * height )  THEN
    11591161                   z_i(1) = zw(k)
    11601162                ELSE
     
    11661168       ELSE
    11671169          DO  k = nzb, nzt-1
    1168              IF ( first .AND. hom(k,1,18,sr) < 0.0 &
    1169                 .AND. abs(hom(k,1,18,sr)) > 1.0E-8 )  THEN
     1170             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
     1171                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
    11701172                first = .FALSE.
    11711173                height = zw(k)
    11721174             ENDIF
    1173              IF ( hom(k,1,18,sr) < 0.0  .AND. &
    1174                   abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
     1175             IF ( hom(k,1,18,sr) < 0.0_wp  .AND.                              &
     1176                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND.                        &
    11751177                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
    1176                 IF ( zw(k) < 1.5 * height )  THEN
     1178                IF ( zw(k) < 1.5_wp * height )  THEN
    11771179                   z_i(1) = zw(k)
    11781180                ELSE
     
    11921194!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
    11931195!--             ocean case!
    1194        z_i(2) = 0.0
     1196       z_i(2) = 0.0_wp
    11951197       DO  k = nzb+1, nzt+1
    11961198          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
     
    12251227!--    the characteristic convective boundary layer temperature.
    12261228!--    The horizontal average at nzb+1 is input for the average temperature.
    1227        IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
    1228            .AND.  z_i(1) /= 0.0 )  THEN
    1229           hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
    1230                                        hom(nzb,1,18,sr) *      &
    1231                                        ABS( z_i(1) ) )**0.333333333
     1229       IF ( hom(nzb,1,18,sr) > 0.0_wp .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8_wp &
     1230           .AND.  z_i(1) /= 0.0_wp )  THEN
     1231          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) *                 &
     1232                                       hom(nzb,1,18,sr) *                      &
     1233                                       ABS( z_i(1) ) )**0.333333333_wp
    12321234!--       so far this only works if Prandtl layer is used
    12331235          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
    12341236       ELSE
    1235           hom(nzb+8,1,pr_palm,sr)  = 0.0
    1236           hom(nzb+11,1,pr_palm,sr) = 0.0
     1237          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
     1238          hom(nzb+11,1,pr_palm,sr) = 0.0_wp
    12371239       ENDIF
    12381240
     
    12611263       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
    12621264
    1263        IF ( ts_value(5,sr) /= 0.0 )  THEN
     1265       IF ( ts_value(5,sr) /= 0.0_wp )  THEN
    12641266          ts_value(22,sr) = ts_value(4,sr)**2 / &
    12651267                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
    12661268       ELSE
    1267           ts_value(22,sr) = 10000.0
     1269          ts_value(22,sr) = 10000.0_wp
    12681270       ENDIF
    12691271
     
    12781280!-- If required, sum up horizontal averages for subsequent time averaging
    12791281    IF ( do_sum )  THEN
    1280        IF ( average_count_pr == 0 )  hom_sum = 0.0
     1282       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
    12811283       hom_sum = hom_sum + hom(:,1,:,:)
    12821284       average_count_pr = average_count_pr + 1
     
    13891391!
    13901392!--    Initialize (local) summation array
    1391        sums_l = 0.0
     1393       sums_l = 0.0_wp
    13921394
    13931395!
     
    14211423             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
    14221424             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
    1423              sums_l(:,34,i) = sums_l(:,34,i) + 0.5 *                        &
    1424                               ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +   &
     1425             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        &
     1426                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
    14251427                                sums_ws2_ws_l(:,i) )    ! e*
    14261428             DO  k = nzb, nzt
    1427                 sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5 * ( &
    1428                                                       sums_us2_ws_l(k,i) +  &
    1429                                                       sums_vs2_ws_l(k,i) +  &
    1430                                                       sums_ws2_ws_l(k,i) )
     1429                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5_wp * ( &
     1430                                                      sums_us2_ws_l(k,i) +     &
     1431                                                      sums_vs2_ws_l(k,i) +     &
     1432                                                      sums_ws2_ws_l(k,i)     )
    14311433             ENDDO
    14321434          ENDDO
     
    15931595       !$acc update host( sums_l )
    15941596       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1595        CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
     1597       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL,  &
    15961598                           MPI_SUM, comm2d, ierr )
    15971599       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1598        CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
     1600       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL,  &
    15991601                           MPI_SUM, comm2d, ierr )
    16001602       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1601        CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
     1603       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL,  &
    16021604                           MPI_SUM, comm2d, ierr )
    16031605       IF ( ocean )  THEN
    16041606          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1605           CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, &
     1607          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb,       &
    16061608                              MPI_REAL, MPI_SUM, comm2d, ierr )
    16071609       ENDIF
    16081610       IF ( humidity ) THEN
    16091611          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1610           CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
     1612          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb,       &
    16111613                              MPI_REAL, MPI_SUM, comm2d, ierr )
    16121614          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1613           CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
     1615          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
    16141616                              MPI_REAL, MPI_SUM, comm2d, ierr )
    16151617          IF ( cloud_physics ) THEN
    16161618             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1617              CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
     1619             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb,    &
    16181620                                 MPI_REAL, MPI_SUM, comm2d, ierr )
    16191621             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1620              CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
     1622             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb,    &
    16211623                                 MPI_REAL, MPI_SUM, comm2d, ierr )
    16221624          ENDIF
     
    16251627       IF ( passive_scalar )  THEN
    16261628          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1627           CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
     1629          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
    16281630                              MPI_REAL, MPI_SUM, comm2d, ierr )
    16291631       ENDIF
     
    17811783          DO  j =  nys, nyn
    17821784             DO  k = nzb, nzt+1
    1783                 s1 = s1 + 0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 + w(k,j,i)**2 ) * &
     1785                s1 = s1 + 0.5_wp *                                             &
     1786                          ( u(k,j,i)**2 + v(k,j,i)**2 + w(k,j,i)**2 ) *        &
    17841787                          rmask(j,i,sr) * rflags_invers(j,i,k+1)
    17851788             ENDDO
     
    18421845!
    18431846!--                Perturbation energy
    1844                    s4 = s4 + 0.5 * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) * &
     1847                   s4 = s4 + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) *  &
    18451848                             rflags_invers(j,i,k+1)
    18461849                ENDDO
     
    18891892!
    18901893!--             Momentum flux w"u"
    1891                 s1 = s1 - 0.25 * (                   &
     1894                s1 = s1 - 0.25_wp * (                                          &
    18921895                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
    18931896                                                           ) * (               &
     
    18981901!
    18991902!--             Momentum flux w"v"
    1900                 s2 = s2 - 0.25 * (                   &
     1903                s2 = s2 - 0.25_wp * (                                          &
    19011904                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
    19021905                                                           ) * (               &
     
    19071910!
    19081911!--             Heat flux w"pt"
    1909                 s3 = s3 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
    1910                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
    1911                               * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     1912                s3 = s3 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                 &
     1913                                 * ( pt(k+1,j,i) - pt(k,j,i) )                 &
     1914                                 * ddzu(k+1) * rmask(j,i,sr)                   &
     1915                                 * rflags_invers(j,i,k+1)
    19121916             ENDDO
    19131917          ENDDO
     
    19261930             DO  i = nxl, nxr
    19271931                DO  j = nys, nyn
    1928                    s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
    1929                                  * ( sa(k+1,j,i) - sa(k,j,i) )   &
    1930                                  * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     1932                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
     1933                                    * ( sa(k+1,j,i) - sa(k,j,i) )              &
     1934                                    * ddzu(k+1) * rmask(j,i,sr)                &
     1935                                    * rflags_invers(j,i,k+1)
    19311936                ENDDO
    19321937             ENDDO
     
    19451950             DO  i = nxl, nxr
    19461951                DO  j = nys, nyn
    1947                    s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
    1948                                  * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
    1949                                  * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    1950                    s2 = s2 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
    1951                                  * ( q(k+1,j,i) - q(k,j,i) )     &
    1952                                  * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     1952                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
     1953                                    * ( vpt(k+1,j,i) - vpt(k,j,i) )            &
     1954                                    * ddzu(k+1) * rmask(j,i,sr)
     1955                                    * rflags_invers(j,i,k+1)
     1956                   s2 = s2 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
     1957                                    * ( q(k+1,j,i) - q(k,j,i) )                &
     1958                                    * ddzu(k+1) * rmask(j,i,sr)                &
     1959                                    * rflags_invers(j,i,k+1)
    19531960                ENDDO
    19541961             ENDDO
     
    19651972                DO  i = nxl, nxr
    19661973                   DO  j = nys, nyn
    1967                       s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )    &
    1968                                     * ( ( q(k+1,j,i) - ql(k+1,j,i) ) &
    1969                                       - ( q(k,j,i) - ql(k,j,i) ) )   &
    1970                                     * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     1974                      s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )           &
     1975                                       *  ( ( q(k+1,j,i) - ql(k+1,j,i) )       &
     1976                                          - ( q(k,j,i) - ql(k,j,i) ) )         &
     1977                                       * ddzu(k+1) * rmask(j,i,sr)             &
     1978                                       * rflags_invers(j,i,k+1)
    19711979                   ENDDO
    19721980                ENDDO
     
    19871995             DO  i = nxl, nxr
    19881996                DO  j = nys, nyn
    1989                    s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
    1990                                  * ( q(k+1,j,i) - q(k,j,i) )     &
    1991                                  * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     1997                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
     1998                                    * ( q(k+1,j,i) - q(k,j,i) )                &
     1999                                    * ddzu(k+1) * rmask(j,i,sr)                &
     2000                                    * rflags_invers(j,i,k+1)
    19922001                ENDDO
    19932002             ENDDO
     
    20102019                s2 = s2 + vsws(j,i) * rmask(j,i,sr)     ! w"v"
    20112020                s3 = s3 + shf(j,i)  * rmask(j,i,sr)     ! w"pt"
    2012                 s4 = s4 + 0.0 * rmask(j,i,sr)           ! u"pt"
    2013                 s5 = s5 + 0.0 * rmask(j,i,sr)           ! v"pt"
     2021                s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
     2022                s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
    20142023             ENDDO
    20152024          ENDDO
     
    20442053                DO  j =  nys, nyn
    20452054                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
    2046                    s2 = s2 + ( ( 1.0 + 0.61 * q(nzb,j,i) ) * shf(j,i) &
    2047                                + 0.61 * pt(nzb,j,i) * qsws(j,i) )
     2055                   s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * shf(j,i)    &
     2056                               + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
    20482057                ENDDO
    20492058             ENDDO
     
    20592068                DO  i = nxl, nxr
    20602069                   DO  j =  nys, nyn
    2061                       s1 = s1 + ( ( 1.0 + 0.61 * q(nzb,j,i) - ql(nzb,j,i) ) * &
    2062                                   shf(j,i) + 0.61 * pt(nzb,j,i) * qsws(j,i) )
     2070                      s1 = s1 + ( ( 1.0_wp +                                   &
     2071                                    0.61_wp * q(nzb,j,i) - ql(nzb,j,i) ) *     &
     2072                                  shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
    20632073                   ENDDO
    20642074                ENDDO
     
    21162126                s2 = s2 + vswst(j,i) * rmask(j,i,sr)    ! w"v"
    21172127                s3 = s3 + tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
    2118                 s4 = s4 + 0.0 * rmask(j,i,sr)           ! u"pt"
    2119                 s5 = s5 + 0.0 * rmask(j,i,sr)           ! v"pt"
     2128                s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
     2129                s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
    21202130             ENDDO
    21212131          ENDDO
     
    21502160                DO  j =  nys, nyn
    21512161                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
    2152                    s2 = s2 + ( ( 1.0 + 0.61 * q(nzt,j,i) ) * tswst(j,i) + &
    2153                                0.61 * pt(nzt,j,i) * qswst(j,i) )
     2162                   s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * tswst(j,i) +&
     2163                                 0.61_wp * pt(nzt,j,i) * qswst(j,i) )
    21542164                ENDDO
    21552165             ENDDO
     
    21652175                DO  i = nxl, nxr
    21662176                   DO  j =  nys, nyn
    2167                       s1 = s1 + ( ( 1.0 + 0.61 * q(nzt,j,i) - ql(nzt,j,i) ) * &
    2168                                   tswst(j,i) + 0.61 * pt(nzt,j,i) * qswst(j,i) )
     2177                      s1 = s1 + ( ( 1.0_wp +                                   &
     2178                                    0.61_wp * q(nzt,j,i) - ql(nzt,j,i) ) *     &
     2179                                  tswst(j,i) +                                 &
     2180                                  0.61_wp * pt(nzt,j,i) * qswst(j,i) )
    21692181                   ENDDO
    21702182                ENDDO
     
    22202232          DO  i = nxl, nxr
    22212233             DO  j = nys, nyn
    2222                 ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
    2223                               u(k+1,j,i) - hom(k+1,1,1,sr) )
    2224                 vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
    2225                               v(k+1,j,i) - hom(k+1,1,2,sr) )
    2226                 pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
    2227                               pt(k+1,j,i) - hom(k+1,1,4,sr) )
     2234                ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) + &
     2235                                 u(k+1,j,i) - hom(k+1,1,1,sr) )
     2236                vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) + &
     2237                                 v(k+1,j,i) - hom(k+1,1,2,sr) )
     2238                pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) + &
     2239                                 pt(k+1,j,i) - hom(k+1,1,4,sr) )
    22282240!
    22292241!--             Higher moments
     
    22322244!
    22332245!--             Energy flux w*e* (has to be adjusted?)
    2234                 s3 = s3 + w(k,j,i) * 0.5 * ( ust**2 + vst**2 + w(k,j,i)**2 ) &
     2246                s3 = s3 + w(k,j,i) * 0.5_wp * ( ust**2 + vst**2 + w(k,j,i)**2 )&
    22352247                                   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    22362248             ENDDO
     
    22542266                DO  i = nxl, nxr
    22552267                   DO  j = nys, nyn
    2256                       s1 = s1 + 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) +   &
    2257                                         sa(k+1,j,i) - hom(k+1,1,23,sr) ) &
    2258                                     * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2268                      s1 = s1 + 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +      &
     2269                                           sa(k+1,j,i) - hom(k+1,1,23,sr) )    &
     2270                                       * w(k,j,i) * rmask(j,i,sr)              &
     2271                                       * rflags_invers(j,i,k+1)
    22592272                   ENDDO
    22602273                ENDDO
     
    22932306                DO  i = nxl, nxr
    22942307                   DO  j = nys, nyn
    2295                       s1 = s1 + 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +     &
    2296                                         vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
    2297                                       w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2308                      s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +     &
     2309                                           vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
     2310                                         w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    22982311                   ENDDO
    22992312                ENDDO
     
    23092322                   DO  i = nxl, nxr
    23102323                      DO  j = nys, nyn
    2311                          s1 = s1 + 0.5 * ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) +   &
    2312                                            ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
    2313                                        * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2324                         s1 = s1 + 0.5_wp * ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) +   &
     2325                                              ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
     2326                                          * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    23142327                      ENDDO
    23152328                   ENDDO
     
    23952408                   DO  i = nxl, nxr
    23962409                      DO  j = nys, nyn
    2397                          s1 = s1 + 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +   &
    2398                                            vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
    2399                                        * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2410                         s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +   &
     2411                                              vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
     2412                                          * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    24002413                      ENDDO
    24012414                   ENDDO
     
    24082421                !$acc parallel loop present( hom, sums_l )
    24092422                DO  k = nzb, nzt_diff
    2410                    sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
    2411                                              0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
     2423                   sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
     2424                                                0.61_wp * hom(k,1,4,sr) * sums_l(k,49,tn)
    24122425                ENDDO
    24132426                !$acc end parallel loop
     
    24272440             DO  i = nxl, nxr
    24282441                DO  j = nys, nyn
    2429                    s1 = s1 + 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) +   &
    2430                                      q(k+1,j,i) - hom(k+1,1,41,sr) ) &
    2431                                  * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2442                   s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +          &
     2443                                        q(k+1,j,i) - hom(k+1,1,41,sr) )        &
     2444                                    * w(k,j,i) * rmask(j,i,sr)                 &
     2445                                    * rflags_invers(j,i,k+1)
    24322446                ENDDO
    24332447             ENDDO
     
    24502464             DO  i = nxl, nxr
    24512465                DO  j = nys, nyn
    2452                    ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
    2453                                u(k+1,j,i) - hom(k+1,1,1,sr) )
    2454                    vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
    2455                                v(k+1,j,i) - hom(k+1,1,2,sr) )
     2466                   ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +              &
     2467                                    u(k+1,j,i) - hom(k+1,1,1,sr) )
     2468                   vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +              &
     2469                                    v(k+1,j,i) - hom(k+1,1,2,sr) )
    24562470!
    24572471!--                Momentum flux w*u*
    2458                    s1 = s1 + 0.5 * ( w(k,j,i-1) + w(k,j,i) ) &
    2459                                  * ust * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2472                   s1 = s1 + 0.5_wp * ( w(k,j,i-1) + w(k,j,i) )                &
     2473                                    * ust * rmask(j,i,sr)                      &
     2474                                    * rflags_invers(j,i,k+1)
    24602475!
    24612476!--                Momentum flux w*v*
    2462                    s2 = s2 + 0.5 * ( w(k,j-1,i) + w(k,j,i) ) &
    2463                                  * vst * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2477                   s2 = s2 + 0.5_wp * ( w(k,j-1,i) + w(k,j,i) )                &
     2478                                    * vst * rmask(j,i,sr)                      &
     2479                                    * rflags_invers(j,i,k+1)
    24642480                ENDDO
    24652481             ENDDO
     
    24812497!
    24822498!--                Vertical heat flux
    2483                    s1 = s1 + 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
    2484                                      pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
    2485                                  * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2499                   s1 = s1 + 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +          &
     2500                                        pt(k+1,j,i) - hom(k+1,1,4,sr) )        &
     2501                                    * w(k,j,i) * rmask(j,i,sr)                 &
     2502                                    * rflags_invers(j,i,k+1)
    24862503                ENDDO
    24872504             ENDDO
     
    24972514                DO  i = nxl, nxr
    24982515                   DO  j = nys, nyn
    2499                       s1 = s1 + 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) +   &
    2500                                         q(k+1,j,i) - hom(k+1,1,41,sr) ) &
    2501                                     * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2516                      s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +       &
     2517                                           q(k+1,j,i) - hom(k+1,1,41,sr) )     &
     2518                                       * w(k,j,i) * rmask(j,i,sr)              &
     2519                                       * rflags_invers(j,i,k+1)
    25022520                   ENDDO
    25032521                ENDDO
     
    25252543!--    First calculate the products, then the divergence.
    25262544!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
    2527        IF ( hom(nzb+1,2,55,0) /= 0.0  .OR.  hom(nzb+1,2,68,0) /= 0.0 )  THEN
     2545       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )  THEN
    25282546
    25292547          STOP '+++ openACC porting for vertical flux div of resolved scale TKE in flow_statistics is still missing'
    2530           sums_ll = 0.0  ! local array
     2548          sums_ll = 0.0_wp  ! local array
    25312549
    25322550          !$OMP DO
     
    25352553                DO  k = nzb_s_inner(j,i)+1, nzt
    25362554
    2537                    sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * (        &
    2538                   ( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)   &
    2539                            - 0.5 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
    2540                            ) )**2                                          &
    2541                 + ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)   &
    2542                            - 0.5 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
    2543                            ) )**2                                          &
    2544                    + w(k,j,i)**2                                  )
    2545 
    2546                    sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
     2555                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
     2556                  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)    &
     2557                              - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )   &
     2558                              ) )**2                                           &
     2559                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)    &
     2560                              - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )   &
     2561                              ) )**2                                           &
     2562                + w(k,j,i)**2                                        )
     2563
     2564                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)            &
    25472565                                               * ( p(k,j,i) + p(k+1,j,i) )
    25482566
     
    25502568             ENDDO
    25512569          ENDDO
    2552           sums_ll(0,1)     = 0.0    ! because w is zero at the bottom
    2553           sums_ll(nzt+1,1) = 0.0
    2554           sums_ll(0,2)     = 0.0
    2555           sums_ll(nzt+1,2) = 0.0
     2570          sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
     2571          sums_ll(nzt+1,1) = 0.0_wp
     2572          sums_ll(0,2)     = 0.0_wp
     2573          sums_ll(nzt+1,2) = 0.0_wp
    25562574
    25572575          DO  k = nzb+1, nzt
     
    25622580          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
    25632581          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
    2564           sums_l(nzb,68,tn) = 0.0    ! because w* = 0 at nzb
     2582          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
    25652583
    25662584       ENDIF
     
    25682586!
    25692587!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
    2570        IF ( hom(nzb+1,2,57,0) /= 0.0  .OR.  hom(nzb+1,2,69,0) /= 0.0 )  THEN
     2588       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )  THEN
    25712589
    25722590          STOP '+++ openACC porting for vertical flux div of SGS TKE in flow_statistics is still missing'
     
    25762594                DO  k = nzb_s_inner(j,i)+1, nzt
    25772595
    2578                    sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * (                 &
     2596                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
    25792597                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
    25802598                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
    2581                                                              ) * ddzw(k)
    2582 
    2583                    sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * (                 &
     2599                                                                ) * ddzw(k)
     2600
     2601                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
    25842602                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
    2585                                                               )
     2603                                                                )
    25862604
    25872605                ENDDO
     
    25962614!--    Horizontal heat fluxes (subgrid, resolved, total).
    25972615!--    Do it only, if profiles shall be plotted.
    2598        IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
     2616       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
    25992617
    26002618          STOP '+++ openACC porting for horizontal flux calculation in flow_statistics is still missing'
     
    26052623!
    26062624!--                Subgrid horizontal heat fluxes u"pt", v"pt"
    2607                    sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 *                   &
     2625                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
    26082626                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
    26092627                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
    26102628                                                 * ddx * rmask(j,i,sr)
    2611                    sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 *                   &
     2629                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
    26122630                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
    26132631                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
     
    26162634!--                Resolved horizontal heat fluxes u*pt*, v*pt*
    26172635                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
    2618                                                   ( u(k,j,i) - hom(k,1,1,sr) ) &
    2619                                        * 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
    2620                                                  pt(k,j,i)   - hom(k,1,4,sr) )
    2621                    pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
    2622                                  pt(k,j,i)   - hom(k,1,4,sr) )
     2636                                     ( u(k,j,i) - hom(k,1,1,sr) ) * 0.5_wp *  &
     2637                                     ( pt(k,j,i-1) - hom(k,1,4,sr) +          &
     2638                                       pt(k,j,i)   - hom(k,1,4,sr) )
     2639                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
     2640                                    pt(k,j,i)   - hom(k,1,4,sr) )
    26232641                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
    2624                                                   ( v(k,j,i) - hom(k,1,2,sr) ) &
    2625                                        * 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
    2626                                                  pt(k,j,i)   - hom(k,1,4,sr) )
     2642                                     ( v(k,j,i) - hom(k,1,2,sr) ) * 0.5_wp *  &
     2643                                     ( pt(k,j-1,i) - hom(k,1,4,sr) +          &
     2644                                       pt(k,j,i)   - hom(k,1,4,sr) )
    26272645                ENDDO
    26282646             ENDDO
     
    26302648!
    26312649!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
    2632           sums_l(nzb,58,tn) = 0.0
    2633           sums_l(nzb,59,tn) = 0.0
    2634           sums_l(nzb,60,tn) = 0.0
    2635           sums_l(nzb,61,tn) = 0.0
    2636           sums_l(nzb,62,tn) = 0.0
    2637           sums_l(nzb,63,tn) = 0.0
     2650          sums_l(nzb,58,tn) = 0.0_wp
     2651          sums_l(nzb,59,tn) = 0.0_wp
     2652          sums_l(nzb,60,tn) = 0.0_wp
     2653          sums_l(nzb,61,tn) = 0.0_wp
     2654          sums_l(nzb,62,tn) = 0.0_wp
     2655          sums_l(nzb,63,tn) = 0.0_wp
    26382656
    26392657       ENDIF
     
    27522770       hom(:,1,37,sr) = sums(:,37)     ! w*e*
    27532771       hom(:,1,38,sr) = sums(:,38)     ! w*3
    2754        hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20 )**1.5   ! Sw
     2772       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
    27552773       hom(:,1,40,sr) = sums(:,40)     ! p
    27562774       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
     
    27812799       hom(:,1,70,sr) = sums(:,70)     ! q*2
    27822800       hom(:,1,71,sr) = sums(:,71)     ! prho
    2783        hom(:,1,72,sr) = hyp * 1E-4     ! hyp in dbar
     2801       hom(:,1,72,sr) = hyp * 1E-4_wp     ! hyp in dbar
    27842802       hom(:,1,73,sr) = sums(:,73)     ! nr
    27852803       hom(:,1,74,sr) = sums(:,74)     ! qr
     
    28092827!--    is less than 1.5 times the height where the heat flux becomes negative
    28102828!--    (positive) for the first time.
    2811        z_i(1) = 0.0
     2829       z_i(1) = 0.0_wp
    28122830       first = .TRUE.
    28132831
    28142832       IF ( ocean )  THEN
    28152833          DO  k = nzt, nzb+1, -1
    2816              IF ( first .AND. hom(k,1,18,sr) < 0.0 &
    2817                 .AND. abs(hom(k,1,18,sr)) > 1.0E-8)  THEN
     2834             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
     2835                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
    28182836                first = .FALSE.
    28192837                height = zw(k)
    28202838             ENDIF
    2821              IF ( hom(k,1,18,sr) < 0.0  .AND. &
    2822                   abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
     2839             IF ( hom(k,1,18,sr) < 0.0_wp  .AND.                              &
     2840                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND.                        &
    28232841                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
    2824                 IF ( zw(k) < 1.5 * height )  THEN
     2842                IF ( zw(k) < 1.5_wp * height )  THEN
    28252843                   z_i(1) = zw(k)
    28262844                ELSE
     
    28322850       ELSE
    28332851          DO  k = nzb, nzt-1
    2834              IF ( first .AND. hom(k,1,18,sr) < 0.0 &
    2835                 .AND. abs(hom(k,1,18,sr)) > 1.0E-8 )  THEN
     2852             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
     2853                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
    28362854                first = .FALSE.
    28372855                height = zw(k)
    28382856             ENDIF
    28392857             IF ( hom(k,1,18,sr) < 0.0  .AND. &
    2840                   abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
     2858                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND. &
    28412859                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
    2842                 IF ( zw(k) < 1.5 * height )  THEN
     2860                IF ( zw(k) < 1.5_wp * height )  THEN
    28432861                   z_i(1) = zw(k)
    28442862                ELSE
     
    28582876!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
    28592877!--             ocean case!
    2860        z_i(2) = 0.0
     2878       z_i(2) = 0.0_wp
    28612879       DO  k = nzb+1, nzt+1
    28622880          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
     
    28912909!--    the characteristic convective boundary layer temperature.
    28922910!--    The horizontal average at nzb+1 is input for the average temperature.
    2893        IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
    2894            .AND.  z_i(1) /= 0.0 )  THEN
    2895           hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
    2896                                        hom(nzb,1,18,sr) *      &
    2897                                        ABS( z_i(1) ) )**0.333333333
     2911       IF ( hom(nzb,1,18,sr) > 0.0_wp .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8_wp &
     2912           .AND.  z_i(1) /= 0.0_wp )  THEN
     2913          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) *                 &
     2914                                       hom(nzb,1,18,sr) *                      &
     2915                                       ABS( z_i(1) ) )**0.333333333_wp
    28982916!--       so far this only works if Prandtl layer is used
    28992917          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
    29002918       ELSE
    2901           hom(nzb+8,1,pr_palm,sr)  = 0.0
    2902           hom(nzb+11,1,pr_palm,sr) = 0.0
     2919          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
     2920          hom(nzb+11,1,pr_palm,sr) = 0.0_wp
    29032921       ENDIF
    29042922
     
    29272945       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
    29282946
    2929        IF ( ts_value(5,sr) /= 0.0 )  THEN
    2930           ts_value(22,sr) = ts_value(4,sr)**2 / &
     2947       IF ( ts_value(5,sr) /= 0.0_wp )  THEN
     2948          ts_value(22,sr) = ts_value(4,sr)**2_wp / &
    29312949                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
    29322950       ELSE
    2933           ts_value(22,sr) = 10000.0
     2951          ts_value(22,sr) = 10000.0_wp
    29342952       ENDIF
    29352953
     
    29462964!-- If required, sum up horizontal averages for subsequent time averaging
    29472965    IF ( do_sum )  THEN
    2948        IF ( average_count_pr == 0 )  hom_sum = 0.0
     2966       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
    29492967       hom_sum = hom_sum + hom(:,1,:,:)
    29502968       average_count_pr = average_count_pr + 1
Note: See TracChangeset for help on using the changeset viewer.