Index: /palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
===================================================================
--- /palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4350)
+++ /palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4351)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! added output of wu, wv, wtheta and wq to enable covariance calculation
+! according to temporal EC method
+! 
+! 4346 2019-12-18 11:55:56Z motisi
 ! Introduction of wall_flags_total_0, which currently sets bits based on static
 ! topography information used in wall_flags_static_0
@@ -85,4 +89,5 @@
         ONLY:  ddzu,                                                           &
                pt,                                                             &
+               q,                                                              &
                u,                                                              &
                v,                                                              &
@@ -145,10 +150,18 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti     !< rotation(u,v,w) aka turbulence intensity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti_av  !< avg. rotation(u,v,w) aka turbulence intensity    
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu     !< uu
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu_av  !< mean of uu    
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv     !< vv
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv_av  !< mean of vv    
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ww     !< ww
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ww_av  !< mean of ww
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu         !< uu
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu_av      !< mean of uu
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv         !< vv
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv_av      !< mean of vv
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ww         !< ww
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ww_av      !< mean of ww
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wu         !< wu
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wu_av      !< mean of wu
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wv         !< wv
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wv_av      !< mean of wv
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wtheta     !< wtheta
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wtheta_av  !< mean of wtheta
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wq         !< wq
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wq_av      !< mean of wq
 
 
@@ -275,4 +288,28 @@
              ENDIF
              ww_av = 0.0_wp
+           
+           CASE ( 'wu' )
+             IF ( .NOT. ALLOCATED( wu_av ) )  THEN
+                ALLOCATE( wu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+             ENDIF
+             wu_av = 0.0_wp
+           
+           CASE ( 'wv' )
+             IF ( .NOT. ALLOCATED( wv_av ) )  THEN
+                ALLOCATE( wv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+             ENDIF
+             wv_av = 0.0_wp
+                
+           CASE ( 'wtheta' )
+             IF ( .NOT. ALLOCATED( wtheta_av ) )  THEN
+                ALLOCATE( wtheta_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+             ENDIF
+             wtheta_av = 0.0_wp
+                
+           CASE ( 'wq' )
+             IF ( .NOT. ALLOCATED( wq_av ) )  THEN
+                ALLOCATE( wq_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+             ENDIF
+             wq_av = 0.0_wp
              
           CASE ( 'theta_2m*' )
@@ -341,4 +378,48 @@
              ENDIF
              
+          CASE ( 'wu' )
+             IF ( ALLOCATED( wu_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         wu_av(k,j,i) = wu_av(k,j,i) + wu(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             
+          CASE ( 'wv' )
+             IF ( ALLOCATED( wv_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         wv_av(k,j,i) = wv_av(k,j,i) + wv(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             
+          CASE ( 'wtheta' )
+             IF ( ALLOCATED( wtheta_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         wtheta_av(k,j,i) = wtheta_av(k,j,i) + wtheta(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+                          
+          CASE ( 'wq' )
+             IF ( ALLOCATED( wq_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         wq_av(k,j,i) = wq_av(k,j,i) + wq(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             
           CASE ( 'theta_2m*' )
              IF ( ALLOCATED( pt_2m_av ) ) THEN
@@ -407,4 +488,48 @@
                       DO  k = nzb, nzt+1
                          ww_av(k,j,i) = ww_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          CASE ( 'wu' )
+             IF ( ALLOCATED( wu_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         wu_av(k,j,i) = wu_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             
+          CASE ( 'wv' )
+             IF ( ALLOCATED( wv_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         wv_av(k,j,i) = wv_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             
+          CASE ( 'wtheta' )
+             IF ( ALLOCATED( wtheta_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         wtheta_av(k,j,i) = wtheta_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             
+          CASE ( 'wq' )
+             IF ( ALLOCATED( wq_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         wq_av(k,j,i) = wq_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                       ENDDO
                    ENDDO
@@ -468,4 +593,16 @@
        CASE ( 'ww' )
           unit = 'm2/s2'
+          
+       CASE ( 'wu' )
+          unit = 'm2/s2'
+             
+       CASE ( 'wv' )
+          unit = 'm2/s2'
+             
+       CASE ( 'wtheta' )
+          unit = 'Km/s'
+             
+       CASE ( 'wq' )
+          unit = 'm/s'
 !
 !--    Treat horizotal cross-section output quanatities
@@ -512,5 +649,9 @@
 !
 !--    s grid
-       CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz'  )
+       CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz',                                 &
+              'wu', 'wu_xy', 'wu_xz', 'wu_yz',                                 &
+              'wv', 'wv_xy', 'wv_xz', 'wv_yz',                                 &
+              'wtheta', 'wtheta_xy', 'wtheta_xz', 'wtheta_yz',                 &
+              'wq', 'wq_xy', 'wq_xz', 'wq_yz')
 
           grid_x = 'x'
@@ -650,4 +791,60 @@
           ENDIF
           flag_nr = 3
+          
+          IF ( mode == 'xy' )  grid = 'zw'
+          
+       CASE ( 'wu_xy', 'wu_xz', 'wu_yz' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wu
+          ELSE
+             IF ( .NOT. ALLOCATED( wu_av ) ) THEN
+                ALLOCATE( wu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wu_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => wu_av
+          ENDIF
+          flag_nr = 0
+          
+          IF ( mode == 'xy' )  grid = 'zw'
+          
+       CASE ( 'wv_xy', 'wv_xz', 'wv_yz' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wv
+          ELSE
+             IF ( .NOT. ALLOCATED( wv_av ) ) THEN
+                ALLOCATE( wv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wv_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => wv_av
+          ENDIF
+          flag_nr = 0
+          
+          IF ( mode == 'xy' )  grid = 'zw'
+          
+       CASE ( 'wtheta_xy', 'wtheta_xz', 'wtheta_yz' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wtheta
+          ELSE
+             IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN
+                ALLOCATE( wtheta_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wtheta_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => wtheta_av
+          ENDIF
+          flag_nr = 0
+          
+          IF ( mode == 'xy' )  grid = 'zw'
+          
+       CASE ( 'wq_xy', 'wq_xz', 'wq_yz' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wq
+          ELSE
+             IF ( .NOT. ALLOCATED( wq_av ) ) THEN
+                ALLOCATE( wq_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wq_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => wq_av
+          ENDIF
+          flag_nr = 0
           
           IF ( mode == 'xy' )  grid = 'zw'
@@ -801,4 +998,52 @@
           flag_nr = 3
 
+       CASE ( 'wu' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wu
+          ELSE
+             IF ( .NOT. ALLOCATED( wu_av ) ) THEN
+                ALLOCATE( wu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wu_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => wu_av
+          ENDIF
+          flag_nr = 0
+
+       CASE ( 'wv' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wv
+          ELSE
+             IF ( .NOT. ALLOCATED( wv_av ) ) THEN
+                ALLOCATE( wv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wv_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => wv_av
+          ENDIF
+          flag_nr = 0
+          
+       CASE ( 'wtheta' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wtheta
+          ELSE
+             IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN
+                ALLOCATE( wtheta_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wtheta_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => wtheta_av
+          ENDIF
+          flag_nr = 0
+          
+       CASE ( 'wq' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wq
+          ELSE
+             IF ( .NOT. ALLOCATED( wq_av ) ) THEN
+                ALLOCATE( wq_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wq_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => wq_av
+          ENDIF
+          flag_nr = 0
+
        CASE DEFAULT
           found = .FALSE.
@@ -899,4 +1144,40 @@
           grid = 'w'
           flag_nr = 3
+    
+       CASE ( 'wu' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wu
+          ELSE
+             to_be_resorted => wu_av
+          ENDIF
+          grid = 's'
+          flag_nr = 0
+    
+       CASE ( 'wv' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wv
+          ELSE
+             to_be_resorted => wv_av
+          ENDIF
+          grid = 's'
+          flag_nr = 0
+              
+       CASE ( 'wtheta' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wtheta
+          ELSE
+             to_be_resorted => wtheta_av
+          ENDIF
+          grid = 's'
+          flag_nr = 0
+               
+       CASE ( 'wq' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wq
+          ELSE
+             to_be_resorted => wq_av
+          ENDIF
+          grid = 's'
+          flag_nr = 0
 
        CASE DEFAULT
@@ -1000,4 +1281,32 @@
                 ALLOCATE( ww(nzb:nzt+1,nys:nyn,nxl:nxr) )
                 ww = 0.0_wp
+             ENDIF
+!
+!--       Allocate array for wu
+          CASE ( 'wu' )
+             IF ( .NOT. ALLOCATED( wu ) )  THEN
+                ALLOCATE( wu(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wu = 0.0_wp
+             ENDIF
+!
+!--       Allocate array for wv
+          CASE ( 'wv' )
+             IF ( .NOT. ALLOCATED( wv ) )  THEN
+                ALLOCATE( wv(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wv = 0.0_wp
+             ENDIF
+!
+!--       Allocate array for wtheta
+          CASE ( 'wtheta' )
+             IF ( .NOT. ALLOCATED( wtheta ) )  THEN
+                ALLOCATE( wtheta(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wtheta = 0.0_wp
+             ENDIF
+!
+!--       Allocate array for wq
+          CASE ( 'wq' )
+             IF ( .NOT. ALLOCATED( wq ) )  THEN
+                ALLOCATE( wq(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wq = 0.0_wp
              ENDIF
 !
@@ -1086,5 +1395,6 @@
                    DO  k = nzb+1, nzt
                       uu(k,j,i) = u(k,j,i) * u(k,j,i)                          &
-                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1) )
+                       * MERGE( 1.0_wp, 0.0_wp,                                &
+                       BTEST( wall_flags_total_0(k,j,i), 1) )
                    ENDDO
                 ENDDO
@@ -1097,5 +1407,6 @@
                    DO  k = nzb+1, nzt
                       vv(k,j,i) = v(k,j,i) * v(k,j,i)                          &
-                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2) )
+                       * MERGE( 1.0_wp, 0.0_wp,                                &
+                       BTEST( wall_flags_total_0(k,j,i), 2) )
                    ENDDO
                 ENDDO
@@ -1108,5 +1419,54 @@
                    DO  k = nzb+1, nzt-1
                       ww(k,j,i) = w(k,j,i) * w(k,j,i)                          &
-                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3) )
+                       * MERGE( 1.0_wp, 0.0_wp,                                &
+                       BTEST( wall_flags_total_0(k,j,i), 3) )
+                   ENDDO
+                ENDDO
+             ENDDO
+!
+!--       wu
+          CASE ( 'wu' )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb+1, nzt-1
+                      wu(k,j,i) = w(k,j,i) * u(k,j,i)                          &
+                       * MERGE( 1.0_wp, 0.0_wp,                                &
+                       BTEST( wall_flags_total_0(k,j,i), 0) )
+                   ENDDO
+                ENDDO
+             ENDDO
+!
+!--       wv
+          CASE ( 'wv' )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb+1, nzt-1
+                      wv(k,j,i) = w(k,j,i) * v(k,j,i)                          &
+                       * MERGE( 1.0_wp, 0.0_wp,                                &
+                       BTEST( wall_flags_total_0(k,j,i), 0) )
+                   ENDDO
+                ENDDO
+             ENDDO
+!
+!--       wtheta
+          CASE ( 'wtheta' )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb+1, nzt-1
+                      wtheta(k,j,i) = w(k,j,i) * pt(k,j,i)                     &
+                       * MERGE( 1.0_wp, 0.0_wp,                                &
+                       BTEST( wall_flags_total_0(k,j,i), 0) )
+                   ENDDO
+                ENDDO
+             ENDDO
+!
+!--       wq
+          CASE ( 'wq' )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb+1, nzt-1
+                      wq(k,j,i) = w(k,j,i) * q(k,j,i)                          &
+                       * MERGE( 1.0_wp, 0.0_wp,                                &
+                       BTEST( wall_flags_total_0(k,j,i), 0) )
                    ENDDO
                 ENDDO
@@ -1492,4 +1852,23 @@
     ENDIF
 
+    IF ( ALLOCATED( wu_av ) )  THEN
+       CALL wrd_write_string( 'wu_av' )
+       WRITE ( 14 )  wu_av
+    ENDIF
+
+    IF ( ALLOCATED( wv_av ) )  THEN
+       CALL wrd_write_string( 'wv_av' )
+       WRITE ( 14 )  wv_av
+    ENDIF
+    
+    IF ( ALLOCATED( wtheta_av ) )  THEN
+       CALL wrd_write_string( 'wtheta_av' )
+       WRITE ( 14 )  wtheta_av
+    ENDIF
+    
+    IF ( ALLOCATED( wq_av ) )  THEN
+       CALL wrd_write_string( 'wq_av' )
+       WRITE ( 14 )  wq_av
+    ENDIF
 
  END SUBROUTINE doq_wrd_local
