Index: /palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
===================================================================
--- /palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4865)
+++ /palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4866)
@@ -19,9 +19,13 @@
 ! Current revisions:
 ! ------------------
-! Implemented air temperature as diagnostic output quantity
+!
 !
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Implemented vertical passive scalar flux
+! Implemented air temperature as diagnostic output quantity (previous commit)
+!
+! 4861 2021-02-01 10:51:45Z raasch
 ! Enable 3D data output also with 64-bit precision
 !
@@ -120,4 +124,5 @@
                pt,                                                                                 &
                q,                                                                                  &
+               s,                                                                                  &
                u,                                                                                  &
                v,                                                                                  &
@@ -139,4 +144,5 @@
                humidity,                                                                           &
                message_string,                                                                     &
+               passive_scalar,                                                                     &
                restart_data_format_output,                                                         &
                varnamelength
@@ -230,4 +236,6 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wq           !< wq
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wq_av        !< mean of wq
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ws           !< ws
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ws_av        !< mean of ws
 
 
@@ -396,4 +404,10 @@
              wq_av = 0.0_wp
 
+          CASE ( 'ws' )
+             IF ( .NOT. ALLOCATED( ws_av ) )  THEN
+                ALLOCATE( ws_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             ENDIF
+             ws_av = 0.0_wp
+
           CASE ( 'theta_2m*' )
              IF ( .NOT. ALLOCATED( pt_2m_av ) )  THEN
@@ -538,4 +552,15 @@
                       DO  k = nzb, nzt+1
                          wq_av(k,j,i) = wq_av(k,j,i) + wq(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          CASE ( 'ws' )
+             IF ( ALLOCATED( ws_av ) )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         ws_av(k,j,i) = ws_av(k,j,i) + ws(k,j,i)
                       ENDDO
                    ENDDO
@@ -703,4 +728,15 @@
              ENDIF
 
+          CASE ( 'ws' )
+             IF ( ALLOCATED( ws_av ) )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         ws_av(k,j,i) = ws_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
          CASE ( 'theta_2m*' )
             IF ( ALLOCATED( pt_2m_av ) )  THEN
@@ -810,4 +846,12 @@
           unit = 'm/s'
 
+       CASE ( 'ws' )
+          IF ( .NOT.  passive_scalar )  THEN
+             message_string = 'data_output = ' // TRIM( var ) //                                   &
+                              ' is not implemented for passive_scalar = .FALSE.'
+             CALL message( 'diagnostic_output', 'PA0093', 1, 2, 0, 6, 0 )
+          ENDIF
+          unit = 'm/s'
+
        CASE ( 'wspeed' )
           unit = 'm/s'
@@ -937,5 +981,6 @@
               'wv', 'wv_xy', 'wv_xz', 'wv_yz',                                                     &
               'wtheta', 'wtheta_xy', 'wtheta_xz', 'wtheta_yz',                                     &
-              'wq', 'wq_xy', 'wq_xz', 'wq_yz' )
+              'wq', 'wq_xy', 'wq_xz', 'wq_yz',                                                     &
+              'ws', 'ws_xy', 'ws_xz', 'ws_yz' )
 
           grid_x = 'x'
@@ -1135,4 +1180,18 @@
           IF ( mode == 'xy' )  grid = 'zw'
 
+       CASE ( 'ws_xy', 'ws_xz', 'ws_yz' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => ws
+          ELSE
+             IF ( .NOT. ALLOCATED( ws_av ) )  THEN
+                ALLOCATE( ws_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ws_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => ws_av
+          ENDIF
+          flag_nr = 0
+
+          IF ( mode == 'xy' )  grid = 'zw'
+
        CASE ( 'theta_2m*_xy' )        ! 2d-array
           IF ( av == 0 )  THEN
@@ -1383,4 +1442,16 @@
           flag_nr = 0
 
+       CASE ( 'ws' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => ws
+          ELSE
+             IF ( .NOT. ALLOCATED( ws_av ) )  THEN
+                ALLOCATE( ws_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ws_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => ws_av
+          ENDIF
+          flag_nr = 0
+
        CASE ( 'wspeed' )
           IF ( av == 0 )  THEN
@@ -1559,4 +1630,13 @@
           ELSE
              to_be_resorted => wq_av
+          ENDIF
+          grid = 's'
+          flag_nr = 0
+
+       CASE ( 'ws' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => ws
+          ELSE
+             to_be_resorted => ws_av
           ENDIF
           grid = 's'
@@ -1728,4 +1808,11 @@
                 ALLOCATE( wq(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 wq = 0.0_wp
+             ENDIF
+!
+!--       Allocate array for ws
+          CASE ( 'ws' )
+             IF ( .NOT. ALLOCATED( ws ) )  THEN
+                ALLOCATE( ws(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ws = 0.0_wp
              ENDIF
 !
@@ -1970,4 +2057,15 @@
              ENDDO
 !
+!--       ws
+          CASE ( 'ws' )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb+1, nzt-1
+                      ws(k,j,i) = w(k,j,i) * 0.5_wp * ( s(k,j,i) + s(k+1,j,i) )                    &
+                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
+                   ENDDO
+                ENDDO
+             ENDDO
+!
 !--       2-m potential temperature
           CASE ( 'theta_2m*' )
@@ -2359,4 +2457,12 @@
                       tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
+       CASE ( 'ws_av' )
+          IF ( .NOT. ALLOCATED( ws_av ) )  THEN
+             ALLOCATE( ws_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( k == 1 )  READ ( 13 )  tmp_3d
+          ws_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
+                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
        CASE ( 'wspeed_av' )
           IF ( .NOT. ALLOCATED( wspeed_av ) )  THEN
@@ -2494,4 +2600,10 @@
     ENDIF
 
+    CALL rd_mpi_io_check_array( 'ws_av' , found = array_found )
+    IF ( array_found )  THEN
+       IF ( .NOT. ALLOCATED( ws_av ) )  ALLOCATE( ws_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+       CALL rrd_mpi_io( 'ws_av', ws_av )
+    ENDIF
+
     CALL rd_mpi_io_check_array( 'wspeed_av' , found = array_found )
     IF ( array_found )  THEN
@@ -2582,4 +2694,9 @@
           CALL wrd_write_string( 'wq_av' )
           WRITE ( 14 )  wq_av
+       ENDIF
+
+       IF ( ALLOCATED( ws_av ) )  THEN
+          CALL wrd_write_string( 'ws_av' )
+          WRITE ( 14 )  ws_av
        ENDIF
 
@@ -2602,4 +2719,5 @@
        IF ( ALLOCATED( wtheta_av ) )    CALL wrd_mpi_io( 'wtheta_av', wtheta_av )
        IF ( ALLOCATED( wq_av ) )        CALL wrd_mpi_io( 'wq_av', wq_av )
+       IF ( ALLOCATED( ws_av ) )        CALL wrd_mpi_io( 'ws_av', ws_av )
        IF ( ALLOCATED( wspeed_av ) )    CALL wrd_mpi_io( 'wspeed_av', wspeed_av )
        IF ( ALLOCATED( wu_av ) )        CALL wrd_mpi_io( 'wu_av', wu_av )
