Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 1919)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 1960)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Separate humidity and passive scalar
 !
 ! Former revisions:
@@ -197,7 +197,8 @@
     USE arrays_3d,                                                             &
         ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, ol, p, prho, prr, pt, q, qc, ql,&
-               qr, qs, qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt,  &
-               td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us,& 
-               usws, uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
+               qr, qs, qsws, qswst, rho, s, sa, ss, ssws, sswst, saswsb,       &
+               saswst, shf, td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q,        &
+               time_vert, ts, tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, &
+               vswst, w, w_subs, zw
         
     USE cloud_parameters,                                                      &
@@ -360,8 +361,8 @@
 
           DO  i = 0, threads_per_task-1
-             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
-             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
-             IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
-                                                   sums_wsqs_ws_l(:,i) !w*q*
+             sums_l(:,17,i)                        = sums_wspts_ws_l(:,i) ! w*pt* 
+             IF ( ocean          ) sums_l(:,66,i)  = sums_wssas_ws_l(:,i) ! w*sa*
+             IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)  ! w*q*
+             IF ( passive_scalar ) sums_l(:,116,i) = sums_wsss_ws_l(:,i)  ! w*s*
           ENDDO
 
@@ -440,5 +441,5 @@
              DO  j =  nys, nyn
                 DO  k = nzb_s_inner(j,i), nzt+1
-                   sums_l(k,41,tn)  = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)
+                   sums_l(k,117,tn)  = sums_l(k,117,tn) + s(k,j,i) * rmask(j,i,sr)
                 ENDDO
              ENDDO
@@ -465,5 +466,5 @@
              ENDIF
              IF ( passive_scalar )  THEN
-                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
+                sums_l(:,117,0) = sums_l(:,117,0) + sums_l(:,117,i)
              ENDIF
           ENDDO
@@ -506,5 +507,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,117,0), sums(nzb,117), nzt+2-nzb,       &
                               MPI_REAL, MPI_SUM, comm2d, ierr )
        ENDIF
@@ -522,5 +523,5 @@
           ENDIF
        ENDIF
-       IF ( passive_scalar )  sums(:,41) = sums_l(:,41,0)
+       IF ( passive_scalar )  sums(:,117) = sums_l(:,117,0)
 #endif
 
@@ -560,6 +561,6 @@
 !
 !--    Passive scalar
-       IF ( passive_scalar )  hom(:,1,41,sr) = sums(:,41) /                    &
-            ngp_2dh_s_inner(:,sr)                    ! s (q)
+       IF ( passive_scalar )  hom(:,1,117,sr) = sums(:,117) /                  &
+            ngp_2dh_s_inner(:,sr)                    ! s 
 
 !
@@ -598,5 +599,8 @@
                                   ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
                 ENDIF
-
+                IF ( passive_scalar )  THEN
+                   sums_l(k,118,tn) = sums_l(k,118,tn) + &
+                                  ( s(k,j,i)-hom(k,1,117,sr) )**2 * rmask(j,i,sr)
+                ENDIF
 !
 !--             Higher moments
@@ -627,4 +631,8 @@
                 sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
                                             qs(j,i)   * rmask(j,i,sr)
+             ENDIF
+             IF ( passive_scalar )  THEN
+                sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +        &
+                                            ss(j,i)   * rmask(j,i,sr)
              ENDIF
           ENDDO
@@ -738,7 +746,7 @@
 !--             Passive scalar flux
                 IF ( passive_scalar )  THEN
-                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
+                   sums_l(k,119,tn) = sums_l(k,119,tn)                         &
                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
-                                               * ( q(k+1,j,i) - q(k,j,i) )     &
+                                               * ( s(k+1,j,i) - s(k,j,i) )     &
                                                * ddzu(k+1) * rmask(j,i,sr)
                 ENDIF
@@ -784,6 +792,6 @@
                 ENDIF
                 IF ( passive_scalar )  THEN
-                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
-                                       qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
+                   sums_l(nzb,119,tn) = sums_l(nzb,119,tn) +                     &
+                                       ssws(j,i) * rmask(j,i,sr) ! w"s"
                 ENDIF
              ENDIF
@@ -863,6 +871,6 @@
                 ENDIF
                 IF ( passive_scalar )  THEN
-                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
-                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
+                   sums_l(nzt,119,tn) = sums_l(nzt,119,tn) + &
+                                       sswst(j,i) * rmask(j,i,sr) ! w"s"
                 ENDIF
              ENDIF
@@ -953,7 +961,7 @@
                 IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca                &
                      .OR. sr /= 0 ) )  THEN
-                   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) *        &
+                   pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +             &
+                                    s(k+1,j,i) - hom(k+1,1,117,sr) )
+                   sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *      &
                                                        rmask(j,i,sr)
                 ENDIF
@@ -1011,4 +1019,10 @@
                                       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
+                  IF ( passive_scalar )  THEN
+                     pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +            &
+                                      s(k+1,j,i) - hom(k+1,1,117,sr) )
+                     sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *      &
                                        rmask(j,i,sr)
                   ENDIF
@@ -1279,4 +1293,6 @@
           sums(k,81:88)         = sums(k,81:88)         / ngp_2dh(sr)
           sums(k,89:114)        = sums(k,89:114)        / ngp_2dh(sr)
+          sums(k,116)           = sums(k,116)           / ngp_2dh(sr)
+          sums(k,119)           = sums(k,119)           / ngp_2dh(sr)
           IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
              sums(k,8:11)          = sums(k,8:11)          / ngp_2dh_s_inner(k,sr)
@@ -1287,5 +1303,6 @@
              sums(k,64)            = sums(k,64)            / ngp_2dh_s_inner(k,sr)
              sums(k,70:80)         = sums(k,70:80)         / ngp_2dh_s_inner(k,sr)
-             sums(k,115:pr_palm-2) = sums(k,115:pr_palm-2) / ngp_2dh_s_inner(k,sr)
+             sums(k,118)           = sums(k,118)           / ngp_2dh_s_inner(k,sr)
+             sums(k,120:pr_palm-2) = sums(k,120:pr_palm-2) / ngp_2dh_s_inner(k,sr)
           ENDIF
        ENDDO
@@ -1298,4 +1315,6 @@
                                     ngp_2dh(sr)
        sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
+                                    ngp_2dh(sr)
+       sums(nzb+13,pr_palm)       = sums(nzb+13,pr_palm)       / &    ! ss
                                     ngp_2dh(sr)
 !--    eges, e*
@@ -1443,4 +1462,10 @@
        hom(:,1,114,sr) = sums(:,114)            !: L
 
+       IF ( passive_scalar )  THEN
+          hom(:,1,119,sr) = sums(:,119)     ! w"s"
+          hom(:,1,116,sr) = sums(:,116)     ! w*s*
+          hom(:,1,120,sr) = sums(:,119) + sums(:,116)    ! ws 
+       ENDIF
+
        hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
                                        ! u*, w'u', w'v', t* (in last profile)
@@ -1589,4 +1614,9 @@
        ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
 
+       IF ( passive_scalar )  THEN
+          ts_value(24,sr) = hom(nzb+13,1,119,sr)       ! w"s" ( to do ! )
+          ts_value(25,sr) = hom(nzb+13,1,pr_palm,sr)   ! s*
+       ENDIF
+
 !
 !--    Collect land surface model timeseries
@@ -1660,7 +1690,8 @@
     USE arrays_3d,                                                             &
         ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
-               qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q,&
-               td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
-               uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
+               qsws, qswst, rho, s, sa, saswsb, saswst, shf, ss, ssws, sswst,  &
+               td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts,      &
+               tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, vswst, w,      &
+               w_subs, zw
                  
         
@@ -1833,6 +1864,6 @@
              sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
              IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
-             IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
-                                                   sums_wsqs_ws_l(:,i) !w*q*
+             IF ( humidity       )  sums_l(:,49,i)  = sums_wsqs_ws_l(:,i) !w*q*
+             IF ( passive_scalar )  sums_l(:,116,i) = sums_wsss_ws_l(:,i) !w*s*
           ENDDO
 
@@ -1938,5 +1969,5 @@
        IF ( passive_scalar )  THEN
           !$OMP DO
-          !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l ) create( s1 )
+          !$acc parallel loop gang present( s, rflags_invers, rmask, sums_l ) create( s1 )
           DO  k = nzb, nzt+1
              s1 = 0
@@ -1944,8 +1975,8 @@
              DO  i = nxl, nxr
                 DO  j =  nys, nyn
-                   s1 = s1 + q(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
-                ENDDO
-             ENDDO
-             sums_l(k,41,tn) = s1
+                   s1 = s1 + s(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                ENDDO
+             ENDDO
+             sums_l(k,117,tn) = s1
           ENDDO
           !$acc end parallel loop
@@ -1981,5 +2012,5 @@
              IF ( passive_scalar )  THEN
                 !$acc parallel present( sums_l )
-                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
+                sums_l(:,117,0) = sums_l(:,117,0) + sums_l(:,117,i)
                 !$acc end parallel
              ENDIF
@@ -2024,5 +2055,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,117,0), sums(nzb,117), nzt+2-nzb,     &
                               MPI_REAL, MPI_SUM, comm2d, ierr )
        ENDIF
@@ -2053,5 +2084,5 @@
        IF ( passive_scalar )  THEN
           !$acc parallel present( sums, sums_l )
-          sums(:,41) = sums_l(:,41,0)
+          sums(:,117) = sums_l(:,117,0)
           !$acc end parallel
        ENDIF
@@ -2102,6 +2133,6 @@
        IF ( passive_scalar )  THEN
           !$acc parallel present( hom, ngp_2dh_s_inner, sums )
-          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
-          hom(:,1,41,sr) = sums(:,41)                ! s (q)
+          sums(:,117)     = sums(:,117) / ngp_2dh_s_inner(:,sr)
+          hom(:,1,117,sr) = sums(:,117)                ! s
           !$acc end parallel
        ENDIF
@@ -2230,4 +2261,17 @@
           ENDDO
           sums_l(nzb+12,pr_palm,tn) = s1
+          !$acc end parallel
+       ENDIF
+
+       IF ( passive_scalar )  THEN
+          !$acc parallel present( ss, rmask, sums_l ) create( s1 )
+          s1 = 0
+          !$acc loop vector collapse( 2 ) reduction( +: s1 )
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+                s1 = s1 + ss(j,i) * rmask(j,i,sr)
+             ENDDO
+          ENDDO
+          sums_l(nzb+13,pr_palm,tn) = s1
           !$acc end parallel
        ENDIF
@@ -2411,5 +2455,5 @@
        IF ( passive_scalar )  THEN
 
-          !$acc parallel loop gang present( ddzu, kh, q, rflags_invers, rmask, sums_l ) create( s1 )
+          !$acc parallel loop gang present( ddzu, kh, s, rflags_invers, rmask, sums_l ) create( s1 )
           DO  k = nzb, nzt_diff
              s1 = 0
@@ -2418,10 +2462,10 @@
                 DO  j = nys, nyn
                    s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
-                                    * ( q(k+1,j,i) - q(k,j,i) )                &
+                                    * ( s(k+1,j,i) - s(k,j,i) )                &
                                     * ddzu(k+1) * rmask(j,i,sr)                &
                                     * rflags_invers(j,i,k+1)
                 ENDDO
              ENDDO
-             sums_l(k,48,tn) = s1
+             sums_l(k,119,tn) = s1
           ENDDO
           !$acc end parallel loop
@@ -2532,13 +2576,13 @@
 
              !$OMP DO
-             !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
+             !$acc parallel present( ssws, rmask, sums_l ) create( s1 )
              s1 = 0
              !$acc loop vector collapse( 2 ) reduction( +: s1 )
              DO  i = nxl, nxr
                 DO  j =  nys, nyn
-                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
-                ENDDO
-             ENDDO
-             sums_l(nzb,48,tn) = s1
+                   s1 = s1 + ssws(j,i) * rmask(j,i,sr)  ! w"s"
+                ENDDO
+             ENDDO
+             sums_l(nzb,119,tn) = s1
              !$acc end parallel
 
@@ -2651,13 +2695,13 @@
 
              !$OMP DO
-             !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
+             !$acc parallel present( sswst, rmask, sums_l ) create( s1 )
              s1 = 0
              !$acc loop vector collapse( 2 ) reduction( +: s1 )
              DO  i = nxl, nxr
                 DO  j =  nys, nyn
-                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
-                ENDDO
-             ENDDO
-             sums_l(nzt,48,tn) = s1
+                   s1 = s1 + sswst(j,i) * rmask(j,i,sr) ! w"s"
+                ENDDO
+             ENDDO
+             sums_l(nzt,119,tn) = s1
              !$acc end parallel
 
@@ -2890,5 +2934,5 @@
        IF ( passive_scalar  .AND.  ( .NOT. ws_scheme_sca  .OR.  sr /= 0 ) )  THEN
 
-          !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
+          !$acc parallel loop gang present( hom, s, rflags_invers, rmask, sums_l, w ) create( s1 )
           DO  k = nzb, nzt_diff
              s1 = 0
@@ -2896,6 +2940,6 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   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) )        &
+                   s1 = s1 + 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +          &
+                                        s(k+1,j,i) - hom(k+1,1,117,sr) )        &
                                     * w(k,j,i) * rmask(j,i,sr)                 &
                                     * rflags_invers(j,i,k+1)
@@ -2981,4 +3025,24 @@
                 ENDDO
                 sums_l(k,49,tn) = s1
+             ENDDO
+             !$acc end parallel loop
+
+          ENDIF
+
+          IF ( passive_scalar )  THEN
+
+             !$acc parallel loop gang present( hom, s, rflags_invers, rmask, sums_l, w ) create( s1 )
+             DO  k = nzb, nzt_diff
+                s1 = 0
+                !$acc loop vector collapse( 2 ) reduction( +: s1 )
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      s1 = s1 + 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +      &
+                                           s(k+1,j,i) - hom(k+1,1,117,sr) )    &
+                                       * w(k,j,i) * rmask(j,i,sr)              &
+                                       * rflags_invers(j,i,k+1)
+                   ENDDO
+                ENDDO
+                sums_l(k,116,tn) = s1
              ENDDO
              !$acc end parallel loop
