Index: palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
===================================================================
--- palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4360)
+++ palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4431)
@@ -25,36 +25,39 @@
 ! -----------------
 ! $Id$
+! added wspeed and wdir output; bugfix: set fill_value in case of masked output
+!
+! 4360 2020-01-07 11:25:50Z suehring
 ! 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
-! 
+!
 ! 4331 2019-12-10 18:25:02Z suehring
 ! - Modularize 2-m potential temperature output
 ! - New output for 10-m wind speed
-! 
+!
 ! 4329 2019-12-10 15:46:36Z motisi
 ! Renamed wall_flags_0 to wall_flags_static_0
-! 
+!
 ! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
-! 
+!
 ! 4167 2019-08-16 11:01:48Z suehring
-! Changed behaviour of masked output over surface to follow terrain and ignore 
+! Changed behaviour of masked output over surface to follow terrain and ignore
 ! buildings (J.Resler, T.Gronemeier)
-! 
+!
 ! 4157 2019-08-14 09:19:12Z suehring
 ! Initialization restructured, in order to work also when data output during
-! spin-up is enabled. 
-! 
+! spin-up is enabled.
+!
 ! 4132 2019-08-02 12:34:17Z suehring
 ! Bugfix in masked data output
-! 
+!
 ! 4069 2019-07-01 14:05:51Z Giersch
-! Masked output running index mid has been introduced as a local variable to 
-! avoid runtime error (Loop variable has been modified) in time_integration 
-! 
+! Masked output running index mid has been introduced as a local variable to
+! avoid runtime error (Loop variable has been modified) in time_integration
+!
 ! 4039 2019-06-18 10:32:41Z suehring
 ! - Add output of uu, vv, ww to enable variance calculation according temporal
@@ -64,12 +67,12 @@
 ! - Rename subroutines
 ! - Further modularization
-! 
+!
 ! 3998 2019-05-23 13:38:11Z suehring
-! Bugfix in gathering all output strings 
-! 
+! Bugfix in gathering all output strings
+!
 ! 3995 2019-05-22 18:59:54Z suehring
 ! Avoid compiler warnings about unused variable and fix string operation which
 ! is not allowed with PGI compiler
-! 
+!
 ! 3994 2019-05-22 18:08:09Z suehring
 ! Initial revision
@@ -78,5 +81,5 @@
 ! --------
 ! @author Farah Kanani-Suehring
-! 
+!
 !
 ! Description:
@@ -85,5 +88,5 @@
 !------------------------------------------------------------------------------!
  MODULE diagnostic_output_quantities_mod
- 
+
     USE arrays_3d,                                                             &
         ONLY:  ddzu,                                                           &
@@ -97,5 +100,5 @@
 
     USE basic_constants_and_equations_mod,                                     &
-        ONLY:  kappa
+        ONLY:  kappa, pi
 
     USE control_parameters,                                                    &
@@ -104,5 +107,5 @@
                message_string,                                                 &
                varnamelength
-! 
+!
 !     USE cpulog,                                                                &
 !         ONLY:  cpu_log, log_point
@@ -114,6 +117,6 @@
         ONLY:  nbgp,                                                           &
                nxl,                                                            &
-               nxlg,                                                           & 
-               nxr,                                                            & 
+               nxlg,                                                           &
+               nxr,                                                            &
                nxrg,                                                           &
                nyn,                                                            &
@@ -139,5 +142,5 @@
 
     INTEGER(iwp) ::  timestep_number_at_prev_calc = 0  !< ...at previous diagnostic output calculation
- 
+
     LOGICAL ::  initialized_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is initialized
     LOGICAL ::  prepared_diagnostic_output_quantities = .FALSE.    !< flag indicating whether output is p
@@ -149,19 +152,27 @@
 
     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 ::  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
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti_av  !< avg. rotation(u,v,w) aka turbulence intensity
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  u_center     !< u at center of grid box
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  u_center_av  !< mean of u_center
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu           !< uu
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu_av        !< mean of uu
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wspeed       !< horizontal wind speed
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wspeed_av    !< mean of horizotal wind speed
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  v_center     !< v at center of grid box
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  v_center_av  !< mean of v_center
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv           !< vv
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv_av        !< mean of vv
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wdir         !< wind direction
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wdir_av      !< mean wind direction
+    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
 
 
@@ -178,10 +189,14 @@
            pt_2m_av,                                                           &
            ti_av,                                                              &
+           u_center_av,                                                        &
            uu_av,                                                              &
            uv_10m_av,                                                          &
+           v_center_av,                                                        &
            vv_av,                                                              &
+           wdir_av,                                                            &
+           wspeed_av,                                                          &
            ww_av
-!                                                                              
-!-- Public routines                                                            
+!
+!-- Public routines
     PUBLIC doq_3d_data_averaging,                                              &
            doq_calculate,                                                      &
@@ -198,5 +213,5 @@
     INTERFACE doq_3d_data_averaging
        MODULE PROCEDURE doq_3d_data_averaging
-    END INTERFACE doq_3d_data_averaging       
+    END INTERFACE doq_3d_data_averaging
 
     INTERFACE doq_calculate
@@ -207,21 +222,21 @@
        MODULE PROCEDURE doq_check_data_output
     END INTERFACE doq_check_data_output
-    
+
     INTERFACE doq_define_netcdf_grid
        MODULE PROCEDURE doq_define_netcdf_grid
     END INTERFACE doq_define_netcdf_grid
-    
+
     INTERFACE doq_output_2d
        MODULE PROCEDURE doq_output_2d
     END INTERFACE doq_output_2d
-    
+
     INTERFACE doq_output_3d
        MODULE PROCEDURE doq_output_3d
     END INTERFACE doq_output_3d
-    
+
     INTERFACE doq_output_mask
        MODULE PROCEDURE doq_output_mask
     END INTERFACE doq_output_mask
-      
+
     INTERFACE doq_init
        MODULE PROCEDURE doq_init
@@ -231,9 +246,9 @@
        MODULE PROCEDURE doq_prepare
     END INTERFACE doq_prepare
-    
+
 !     INTERFACE doq_rrd_local
 !        MODULE PROCEDURE doq_rrd_local
 !     END INTERFACE doq_rrd_local
-    
+
     INTERFACE doq_wrd_local
        MODULE PROCEDURE doq_wrd_local
@@ -242,5 +257,5 @@
 
  CONTAINS
- 
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -254,6 +269,6 @@
         ONLY:  average_count_3d
 
-    CHARACTER (LEN=*) ::  mode     !< 
-    CHARACTER (LEN=*) ::  variable !< 
+    CHARACTER (LEN=*) ::  mode     !<
+    CHARACTER (LEN=*) ::  variable !<
 
     INTEGER(iwp) ::  i !<
@@ -270,5 +285,5 @@
              ENDIF
              ti_av = 0.0_wp
-       
+
           CASE ( 'uu' )
              IF ( .NOT. ALLOCATED( uu_av ) )  THEN
@@ -276,5 +291,5 @@
              ENDIF
              uu_av = 0.0_wp
-                
+
           CASE ( 'vv' )
              IF ( .NOT. ALLOCATED( vv_av ) )  THEN
@@ -282,5 +297,5 @@
              ENDIF
              vv_av = 0.0_wp
-                
+
           CASE ( 'ww' )
              IF ( .NOT. ALLOCATED( ww_av ) )  THEN
@@ -288,5 +303,5 @@
              ENDIF
              ww_av = 0.0_wp
-           
+
            CASE ( 'wu' )
              IF ( .NOT. ALLOCATED( wu_av ) )  THEN
@@ -294,5 +309,5 @@
              ENDIF
              wu_av = 0.0_wp
-           
+
            CASE ( 'wv' )
              IF ( .NOT. ALLOCATED( wv_av ) )  THEN
@@ -300,5 +315,5 @@
              ENDIF
              wv_av = 0.0_wp
-                
+
            CASE ( 'wtheta' )
              IF ( .NOT. ALLOCATED( wtheta_av ) )  THEN
@@ -306,5 +321,5 @@
              ENDIF
              wtheta_av = 0.0_wp
-                
+
            CASE ( 'wq' )
              IF ( .NOT. ALLOCATED( wq_av ) )  THEN
@@ -312,5 +327,5 @@
              ENDIF
              wq_av = 0.0_wp
-             
+
           CASE ( 'theta_2m*' )
              IF ( .NOT. ALLOCATED( pt_2m_av ) )  THEN
@@ -318,5 +333,5 @@
              ENDIF
              pt_2m_av = 0.0_wp
-             
+
           CASE ( 'wspeed_10m*' )
              IF ( .NOT. ALLOCATED( uv_10m_av ) )  THEN
@@ -325,4 +340,20 @@
              uv_10m_av = 0.0_wp
 
+          CASE ( 'wspeed' )
+             IF ( .NOT. ALLOCATED( wspeed_av ) )  THEN
+                ALLOCATE( wspeed_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+             ENDIF
+             wspeed_av = 0.0_wp
+
+          CASE ( 'wdir' )
+             IF ( .NOT. ALLOCATED( u_center_av ) )  THEN
+                ALLOCATE( u_center_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+             ENDIF
+             IF ( .NOT. ALLOCATED( v_center_av ) )  THEN
+                ALLOCATE( v_center_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+             ENDIF
+             u_center_av = 0.0_wp
+             v_center_av = 0.0_wp
+
           CASE DEFAULT
              CONTINUE
@@ -333,5 +364,5 @@
 
        SELECT CASE ( TRIM( variable ) )
- 
+
           CASE ( 'ti' )
              IF ( ALLOCATED( ti_av ) ) THEN
@@ -344,5 +375,5 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'uu' )
              IF ( ALLOCATED( uu_av ) ) THEN
@@ -355,5 +386,5 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'vv' )
              IF ( ALLOCATED( vv_av ) ) THEN
@@ -366,5 +397,5 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'ww' )
              IF ( ALLOCATED( ww_av ) ) THEN
@@ -377,5 +408,5 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'wu' )
              IF ( ALLOCATED( wu_av ) ) THEN
@@ -388,5 +419,5 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'wv' )
              IF ( ALLOCATED( wv_av ) ) THEN
@@ -399,5 +430,5 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'wtheta' )
              IF ( ALLOCATED( wtheta_av ) ) THEN
@@ -410,5 +441,5 @@
                 ENDDO
              ENDIF
-                          
+
           CASE ( 'wq' )
              IF ( ALLOCATED( wq_av ) ) THEN
@@ -421,5 +452,5 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'theta_2m*' )
              IF ( ALLOCATED( pt_2m_av ) ) THEN
@@ -436,4 +467,27 @@
                    DO  j = nys, nyn
                       uv_10m_av(j,i) = uv_10m_av(j,i) + uv_10m(j,i)
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          CASE ( 'wspeed' )
+            IF ( ALLOCATED( wspeed_av ) ) THEN
+               DO  i = nxl, nxr
+                  DO  j = nys, nyn
+                     DO  k = nzb, nzt+1
+                         wspeed_av(k,j,i) = wspeed_av(k,j,i) + wspeed(k,j,i)
+                     ENDDO
+                  ENDDO
+               ENDDO
+            ENDIF
+
+          CASE ( 'wdir' )
+             IF ( ALLOCATED( u_center_av )  .AND.  ALLOCATED( v_center_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                        u_center_av(k,j,i) = u_center_av(k,j,i) + u_center(k,j,i)
+                        v_center_av(k,j,i) = v_center_av(k,j,i) + v_center(k,j,i)
+                      ENDDO
                    ENDDO
                 ENDDO
@@ -459,5 +513,5 @@
                 ENDDO
              ENDIF
-       
+
           CASE ( 'uu' )
              IF ( ALLOCATED( uu_av ) ) THEN
@@ -470,5 +524,5 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'vv' )
              IF ( ALLOCATED( vv_av ) ) THEN
@@ -481,5 +535,5 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'ww' )
              IF ( ALLOCATED( ww_av ) ) THEN
@@ -503,5 +557,5 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'wv' )
              IF ( ALLOCATED( wv_av ) ) THEN
@@ -514,5 +568,5 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'wtheta' )
              IF ( ALLOCATED( wtheta_av ) ) THEN
@@ -525,5 +579,5 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'wq' )
              IF ( ALLOCATED( wq_av ) ) THEN
@@ -557,4 +611,35 @@
             ENDIF
 
+         CASE ( 'wspeed' )
+             IF ( ALLOCATED( wspeed_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         wspeed_av(k,j,i) = wspeed_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          CASE ( 'wdir' )
+             IF ( ALLOCATED( u_center_av )  .AND.  ALLOCATED( v_center_av ) ) THEN
+
+                IF ( .NOT. ALLOCATED( wdir_av ) )  THEN
+                   ALLOCATE( wdir_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                ENDIF
+                wdir_av = 0.0_wp
+
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         u_center_av(k,j,i) = u_center_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                         v_center_av(k,j,i) = v_center_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                         wdir_av(k,j,i) = ATAN2( u_center_av(k,j,i), v_center_av(k,j,i) ) &
+                                        / pi * 180.0_wp + 180.0_wp
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
        END SELECT
 
@@ -562,6 +647,6 @@
 
 
- END SUBROUTINE doq_3d_data_averaging 
- 
+ END SUBROUTINE doq_3d_data_averaging
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -573,5 +658,5 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  unit  !< 
+    CHARACTER (LEN=*) ::  unit  !<
     CHARACTER (LEN=*) ::  var   !<
 
@@ -584,25 +669,31 @@
        CASE ( 'ti' )
           unit = '1/s'
-             
+
        CASE ( 'uu' )
           unit = 'm2/s2'
-             
+
        CASE ( 'vv' )
           unit = 'm2/s2'
-             
+
        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'
+
+       CASE ( 'wspeed' )
+          unit = 'm/s'
+
+       CASE ( 'wdir' )
+          unit = 'degree'
 !
 !--    Treat horizotal cross-section output quanatities
@@ -627,5 +718,5 @@
 
  END SUBROUTINE doq_check_data_output
- 
+
 !------------------------------------------------------------------------------!
 !
@@ -633,21 +724,23 @@
 ! ------------
 !> Subroutine defining appropriate grid for netcdf variables.
-!------------------------------------------------------------------------------! 
+!------------------------------------------------------------------------------!
  SUBROUTINE doq_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
-    
+
     IMPLICIT NONE
 
     CHARACTER (LEN=*), INTENT(IN)  ::  variable    !<
-    LOGICAL, INTENT(OUT)           ::  found       !< 
-    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !< 
-    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !< 
-    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !< 
+    LOGICAL, INTENT(OUT)           ::  found       !<
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
 
     found  = .TRUE.
-    
+
     SELECT CASE ( TRIM( variable ) )
 !
 !--    s grid
        CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz',                                 &
+              'wspeed', 'wspeed_xy', 'wspeed_xz', 'wspeed_yz',                 &
+              'wdir', 'wdir_xy', 'wdir_xz', 'wdir_yz',                         &
               'wu', 'wu_xy', 'wu_xz', 'wu_yz',                                 &
               'wv', 'wv_xy', 'wv_xz', 'wv_yz',                                 &
@@ -697,5 +790,5 @@
 
  END SUBROUTINE doq_define_netcdf_grid
- 
+
 !------------------------------------------------------------------------------!
 !
@@ -710,7 +803,7 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  grid     !< 
-    CHARACTER (LEN=*) ::  mode     !< 
-    CHARACTER (LEN=*) ::  variable !< 
+    CHARACTER (LEN=*) ::  grid     !<
+    CHARACTER (LEN=*) ::  mode     !<
+    CHARACTER (LEN=*) ::  variable !<
 
     INTEGER(iwp) ::  av       !< value indicating averaged or non-averaged output
@@ -719,6 +812,6 @@
     INTEGER(iwp) ::  j        !< grid index y-direction
     INTEGER(iwp) ::  k        !< grid index z-direction
-    INTEGER(iwp) ::  nzb_do   !< 
-    INTEGER(iwp) ::  nzt_do   !< 
+    INTEGER(iwp) ::  nzb_do   !<
+    INTEGER(iwp) ::  nzt_do   !<
 
     LOGICAL ::  found             !< true if variable is in list
@@ -727,8 +820,8 @@
 
     REAL(wp) ::  fill_value       !< value for the _FillValue attribute
-    
-    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< 
+
+    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
     REAL(wp), DIMENSION(:,:,:), POINTER ::                 to_be_resorted  !< points to array which needs to be resorted for output
-    
+
     flag_nr  = 0
     found    = .TRUE.
@@ -749,7 +842,7 @@
            ENDIF
            flag_nr = 0
-           
+
            IF ( mode == 'xy' )  grid = 'zu'
-    
+
        CASE ( 'uu_xy', 'uu_xz', 'uu_yz' )
           IF ( av == 0 )  THEN
@@ -763,7 +856,7 @@
           ENDIF
           flag_nr = 1
-          
+
           IF ( mode == 'xy' )  grid = 'zu'
-          
+
        CASE ( 'vv_xy', 'vv_xz', 'vv_yz' )
           IF ( av == 0 )  THEN
@@ -777,7 +870,7 @@
           ENDIF
           flag_nr = 2
-          
+
           IF ( mode == 'xy' )  grid = 'zu'
-                
+
        CASE ( 'ww_xy', 'ww_xz', 'ww_yz' )
           IF ( av == 0 )  THEN
@@ -791,7 +884,7 @@
           ENDIF
           flag_nr = 3
-          
+
           IF ( mode == 'xy' )  grid = 'zw'
-          
+
        CASE ( 'wu_xy', 'wu_xz', 'wu_yz' )
           IF ( av == 0 )  THEN
@@ -805,7 +898,7 @@
           ENDIF
           flag_nr = 0
-          
+
           IF ( mode == 'xy' )  grid = 'zw'
-          
+
        CASE ( 'wv_xy', 'wv_xz', 'wv_yz' )
           IF ( av == 0 )  THEN
@@ -819,7 +912,7 @@
           ENDIF
           flag_nr = 0
-          
+
           IF ( mode == 'xy' )  grid = 'zw'
-          
+
        CASE ( 'wtheta_xy', 'wtheta_xz', 'wtheta_yz' )
           IF ( av == 0 )  THEN
@@ -833,7 +926,7 @@
           ENDIF
           flag_nr = 0
-          
+
           IF ( mode == 'xy' )  grid = 'zw'
-          
+
        CASE ( 'wq_xy', 'wq_xz', 'wq_yz' )
           IF ( av == 0 )  THEN
@@ -847,7 +940,7 @@
           ENDIF
           flag_nr = 0
-          
+
           IF ( mode == 'xy' )  grid = 'zw'
-          
+
        CASE ( 'theta_2m*_xy' )        ! 2d-array
           IF ( av == 0 )  THEN
@@ -871,5 +964,5 @@
           two_d    = .TRUE.
           grid     = 'zu1'
-          
+
        CASE ( 'wspeed_10m*_xy' )        ! 2d-array
           IF ( av == 0 )  THEN
@@ -894,4 +987,32 @@
           grid     = 'zu1'
 
+       CASE ( 'wspeed_xy', 'wspeed_xz', 'wspeed_yz' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wspeed
+          ELSE
+             IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN
+                ALLOCATE( wspeed_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wspeed_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => wspeed_av
+          ENDIF
+          flag_nr = 0
+
+          IF ( mode == 'xy' )  grid = 'zu'
+
+       CASE ( 'wdir_xy', 'wdir_xz', 'wdir_yz' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wdir
+          ELSE
+             IF ( .NOT. ALLOCATED( wdir_av ) ) THEN
+                ALLOCATE( wdir_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wdir_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => wdir_av
+          ENDIF
+          flag_nr = 0
+
+          IF ( mode == 'xy' )  grid = 'zu'
+
        CASE DEFAULT
           found = .FALSE.
@@ -899,6 +1020,6 @@
 
     END SELECT
-    
-    IF ( found  .AND.  .NOT. resorted )  THEN     
+
+    IF ( found  .AND.  .NOT. resorted )  THEN
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -906,13 +1027,13 @@
                 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
                                      REAL( fill_value, KIND = wp ),            &
-                                     BTEST( wall_flags_total_0(k,j,i), flag_nr ) ) 
+                                     BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
              ENDDO
           ENDDO
        ENDDO
     ENDIF
- 
+
  END SUBROUTINE doq_output_2d
- 
- 
+
+
 !------------------------------------------------------------------------------!
 !
@@ -923,8 +1044,8 @@
  SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do,  &
                            nzt_do )
- 
+
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  variable !< 
+    CHARACTER (LEN=*) ::  variable !<
 
     INTEGER(iwp) ::  av       !< index indicating averaged or instantaneous output
@@ -941,5 +1062,5 @@
     REAL(wp) ::  fill_value       !< value for the _FillValue attribute
 
-    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf        !< 
+    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf        !<
     REAL(wp), DIMENSION(:,:,:), POINTER ::                 to_be_resorted  !< points to array which needs to be resorted for output
 
@@ -947,5 +1068,5 @@
     found    = .TRUE.
     resorted = .FALSE.
-    
+
     SELECT CASE ( TRIM( variable ) )
 
@@ -961,5 +1082,5 @@
           ENDIF
           flag_nr = 0
-    
+
        CASE ( 'uu' )
           IF ( av == 0 )  THEN
@@ -973,5 +1094,5 @@
           ENDIF
           flag_nr = 1
-             
+
        CASE ( 'vv' )
           IF ( av == 0 )  THEN
@@ -985,5 +1106,5 @@
           ENDIF
           flag_nr = 2
-             
+
        CASE ( 'ww' )
           IF ( av == 0 )  THEN
@@ -1021,5 +1142,5 @@
           ENDIF
           flag_nr = 0
-          
+
        CASE ( 'wtheta' )
           IF ( av == 0 )  THEN
@@ -1033,5 +1154,5 @@
           ENDIF
           flag_nr = 0
-          
+
        CASE ( 'wq' )
           IF ( av == 0 )  THEN
@@ -1046,10 +1167,34 @@
           flag_nr = 0
 
+       CASE ( 'wspeed' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wspeed
+          ELSE
+             IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN
+                ALLOCATE( wspeed_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wspeed_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => wspeed_av
+          ENDIF
+          flag_nr = 0
+
+       CASE ( 'wdir' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wdir
+          ELSE
+             IF ( .NOT. ALLOCATED( wdir_av ) ) THEN
+                ALLOCATE( wdir_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wdir_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => wdir_av
+          ENDIF
+          flag_nr = 0
+
        CASE DEFAULT
           found = .FALSE.
 
     END SELECT
-    
-    IF ( found  .AND.  .NOT. resorted )  THEN     
+
+    IF ( found  .AND.  .NOT. resorted )  THEN
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -1057,5 +1202,5 @@
                 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
                                      REAL( fill_value, KIND = wp ),            &
-                                     BTEST( wall_flags_total_0(k,j,i), flag_nr ) ) 
+                                     BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
              ENDDO
           ENDDO
@@ -1064,5 +1209,5 @@
 
  END SUBROUTINE doq_output_3d
- 
+
 ! Description:
 ! ------------
@@ -1071,7 +1216,7 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE doq_output_mask( av, variable, found, local_pf, mid )
- 
+
     USE control_parameters
-        
+
     USE indices
 
@@ -1097,5 +1242,5 @@
     REAL(wp),                                                                  &
        DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
-          local_pf   !< 
+          local_pf   !<
     REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to array which needs to be resorted for output
 
@@ -1117,5 +1262,5 @@
           grid = 's'
           flag_nr = 0
-    
+
        CASE ( 'uu' )
           IF ( av == 0 )  THEN
@@ -1126,5 +1271,5 @@
           grid = 'u'
           flag_nr = 1
-    
+
        CASE ( 'vv' )
           IF ( av == 0 )  THEN
@@ -1135,5 +1280,5 @@
           grid = 'v'
           flag_nr = 2
-    
+
        CASE ( 'ww' )
           IF ( av == 0 )  THEN
@@ -1144,5 +1289,5 @@
           grid = 'w'
           flag_nr = 3
-    
+
        CASE ( 'wu' )
           IF ( av == 0 )  THEN
@@ -1153,5 +1298,5 @@
           grid = 's'
           flag_nr = 0
-    
+
        CASE ( 'wv' )
           IF ( av == 0 )  THEN
@@ -1162,5 +1307,5 @@
           grid = 's'
           flag_nr = 0
-              
+
        CASE ( 'wtheta' )
           IF ( av == 0 )  THEN
@@ -1171,5 +1316,5 @@
           grid = 's'
           flag_nr = 0
-               
+
        CASE ( 'wq' )
           IF ( av == 0 )  THEN
@@ -1181,9 +1326,27 @@
           flag_nr = 0
 
+       CASE ( 'wspeed' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wspeed
+          ELSE
+             to_be_resorted => wspeed_av
+          ENDIF
+          grid = 's'
+          flag_nr = 0
+
+       CASE ( 'wdir' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => wdir
+          ELSE
+             to_be_resorted => wdir_av
+          ENDIF
+          grid = 's'
+          flag_nr = 0
+
        CASE DEFAULT
           found = .FALSE.
 
     END SELECT
-    
+
     IF ( found  .AND.  .NOT. resorted )  THEN
        IF ( .NOT. mask_surface(mid) )  THEN
@@ -1193,7 +1356,13 @@
              DO  j = 1, mask_size_l(mid,2)
                 DO  k = 1, mask_size_l(mid,3)
-                   local_pf(i,j,k) =  to_be_resorted(mask_k(mid,k),            &
-                                                     mask_j(mid,j),            &
-                                                     mask_i(mid,i))
+                   local_pf(i,j,k) = MERGE( to_be_resorted(mask_k(mid,k),  &
+                                                           mask_j(mid,j),  &
+                                                           mask_i(mid,i)), &
+                                            REAL( fill_value, KIND = wp ), &
+                                            BTEST( wall_flags_total_0(     &
+                                                           mask_k(mid,k),  &
+                                                           mask_j(mid,j),  &
+                                                           mask_i(mid,i)), &
+                                                   flag_nr ) )
                 ENDDO
              ENDDO
@@ -1226,5 +1395,5 @@
        ENDIF
     ENDIF
-    
+
  END SUBROUTINE doq_output_mask
 
@@ -1237,5 +1406,5 @@
 
     IMPLICIT NONE
-    
+
     INTEGER(iwp) ::  ivar   !< loop index over all 2d/3d/mask output quantities
 
@@ -1248,8 +1417,8 @@
 
     initialized_diagnostic_output_quantities = .FALSE.
-    
+
     ivar = 1
-    
-    DO  WHILE ( ivar <= SIZE( do_all ) ) 
+
+    DO  WHILE ( ivar <= SIZE( do_all ) )
 
        SELECT CASE ( TRIM( do_all(ivar) ) )
@@ -1324,4 +1493,27 @@
                 uv_10m = 0.0_wp
              ENDIF
+!
+!--       Allocate array for wspeed
+          CASE ( 'wspeed' )
+             IF ( .NOT. ALLOCATED( wspeed ) )  THEN
+                ALLOCATE( wspeed(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wspeed = 0.0_wp
+             ENDIF
+
+!
+!--       Allocate array for wdir
+          CASE ( 'wdir' )
+             IF ( .NOT. ALLOCATED( u_center ) )  THEN
+                ALLOCATE( u_center(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                u_center = 0.0_wp
+             ENDIF
+             IF ( .NOT. ALLOCATED( v_center ) )  THEN
+                ALLOCATE( v_center(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                v_center = 0.0_wp
+             ENDIF
+             IF ( .NOT. ALLOCATED( wdir ) )  THEN
+                ALLOCATE( wdir(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                wdir = 0.0_wp
+             ENDIF
 
        END SELECT
@@ -1331,5 +1523,5 @@
 
     initialized_diagnostic_output_quantities = .TRUE.
-    
+
  END SUBROUTINE doq_init
 
@@ -1345,8 +1537,8 @@
 
     INTEGER(iwp) ::  i          !< grid index x-dimension
-    INTEGER(iwp) ::  j          !< grid index y-dimension 
+    INTEGER(iwp) ::  j          !< grid index y-dimension
     INTEGER(iwp) ::  k          !< grid index z-dimension
     INTEGER(iwp) ::  ivar       !< loop index over all 2d/3d/mask output quantities
-    
+
     TYPE(surf_type), POINTER ::  surf     !< surf-type array, used to generalize subroutines
 
@@ -1355,5 +1547,5 @@
 
 !
-!-- Save timestep number to check in time_integration if doq_calculate 
+!-- Save timestep number to check in time_integration if doq_calculate
 !-- has been called already, since the CALL occurs at two locations, but the calculations need to be
 !-- done only once per timestep.
@@ -1362,5 +1554,5 @@
     ivar = 1
 
-    DO  WHILE ( ivar <= SIZE( do_all ) ) 
+    DO  WHILE ( ivar <= SIZE( do_all ) )
 
        SELECT CASE ( TRIM( do_all(ivar) ) )
@@ -1387,5 +1579,5 @@
                    ENDDO
                 ENDDO
-             ENDDO           
+             ENDDO
 !
 !--       uu
@@ -1502,4 +1694,31 @@
              surf => surf_usm_h
              CALL calc_wind_10m
+!
+!--       horizontal wind speed
+          CASE ( 'wspeed' )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb, nzt+1
+                      wspeed(k,j,i) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2             &
+                                          + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 )           &
+                                    * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0) )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+!
+!--       horizontal wind direction
+          CASE ( 'wdir' )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb, nzt+1
+                      u_center(k,j,i) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
+                      v_center(k,j,i) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
+
+                      wdir(k,j,i) = ATAN2( u_center(k,j,i), v_center(k,j,i) ) &
+                                  / pi * 180.0_wp + 180.0_wp
+                   ENDDO
+                ENDDO
+             ENDDO
 
        END SELECT
@@ -1511,6 +1730,6 @@
 
 !
-!-- The following block contains subroutines to calculate diagnostic 
-!-- surface quantities. 
+!-- The following block contains subroutines to calculate diagnostic
+!-- surface quantities.
     CONTAINS
 !------------------------------------------------------------------------------!
@@ -1528,5 +1747,5 @@
           INTEGER(iwp) ::  kk     !< running index along the z-dimension
           INTEGER(iwp) ::  m      !< running index for surface elements
- 
+
           DO  m = 1, surf%ns
 
@@ -1535,5 +1754,5 @@
              k = surf%k(m)
 !
-!--          If 2-m level is below the first grid level, MOST is 
+!--          If 2-m level is below the first grid level, MOST is
 !--          used for calculation of 2-m temperature.
              IF ( surf%z_mo(m) > 2.0_wp )  THEN
@@ -1546,15 +1765,15 @@
 !--          is linearly interpolated between the two nearest vertical grid
 !--          levels. Note, since 2-m temperature is only computed for
-!--          horizontal upward-facing surfaces, only a vertical 
-!--          interpolation is necessary. 
+!--          horizontal upward-facing surfaces, only a vertical
+!--          interpolation is necessary.
              ELSE
 !
-!--             zw(k-1) defines the height of the surface. 
+!--             zw(k-1) defines the height of the surface.
                 kk = k
                 DO WHILE ( zu(kk) - zw(k-1) < 2.0_wp  .AND.  kk <= nzt )
-                   kk = kk + 1 
-                ENDDO
-!
-!--             kk defines the index of the first grid level >= 2m. 
+                   kk = kk + 1
+                ENDDO
+!
+!--             kk defines the index of the first grid level >= 2m.
                 pt_2m(j,i) = pt(kk-1,j,i) +                                    &
                               ( zw(k-1) + 2.0_wp - zu(kk-1)     ) *            &
@@ -1566,5 +1785,5 @@
 
        END SUBROUTINE calc_pt_2m
-       
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -1584,5 +1803,5 @@
           REAL(wp) ::  uv_l !< wind speed at lower grid point
           REAL(wp) ::  uv_u !< wind speed at upper grid point
- 
+
           DO  m = 1, surf%ns
 
@@ -1591,5 +1810,5 @@
              k = surf%k(m)
 !
-!--          If 10-m level is below the first grid level, MOST is 
+!--          If 10-m level is below the first grid level, MOST is
 !--          used for calculation of 10-m temperature.
              IF ( surf%z_mo(m) > 10.0_wp )  THEN
@@ -1602,12 +1821,12 @@
 !--          is linearly interpolated between the two nearest vertical grid
 !--          levels. Note, since 10-m temperature is only computed for
-!--          horizontal upward-facing surfaces, only a vertical 
-!--          interpolation is necessary. 
+!--          horizontal upward-facing surfaces, only a vertical
+!--          interpolation is necessary.
              ELSE
 !
-!--             zw(k-1) defines the height of the surface. 
+!--             zw(k-1) defines the height of the surface.
                 kk = k
                 DO WHILE ( zu(kk) - zw(k-1) < 10.0_wp  .AND.  kk <= nzt )
-                   kk = kk + 1 
+                   kk = kk + 1
                 ENDDO
 !
@@ -1653,5 +1872,5 @@
     INTEGER(iwp) ::  ivar_all   !< loop index
     INTEGER(iwp) ::  l          !< index for cutting string
-    INTEGER(iwp) ::  mid          !< masked output running index 
+    INTEGER(iwp) ::  mid          !< masked output running index
 
     prepared_diagnostic_output_quantities = .FALSE.
@@ -1667,6 +1886,6 @@
 !
 !--    Gather 2d output quantity names.
-!--    Check for double occurrence of output quantity, e.g. by _xy, 
-!--    _yz, _xz. 
+!--    Check for double occurrence of output quantity, e.g. by _xy,
+!--    _yz, _xz.
        DO  WHILE ( do2d_var(av,ivar)(1:1) /= ' ' )
           IF ( .NOT.  ANY( do_all == do2d_var(av,ivar) ) )  THEN
@@ -1709,10 +1928,10 @@
 
  END SUBROUTINE doq_prepare
- 
+
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Subroutine reads local (subdomain) restart data
-!> Note: With the current structure reading of non-standard array is not 
+!> Note: With the current structure reading of non-standard array is not
 !> possible
 !------------------------------------------------------------------------------!
@@ -1720,43 +1939,43 @@
 !                            nxr_on_file, nynf, nync, nyn_on_file, nysf,         &
 !                            nysc, nys_on_file, tmp_3d_non_standard, found )
-!  
-! 
+!
+!
 !     USE control_parameters
-!         
+!
 !     USE indices
-!     
+!
 !     USE kinds
-!     
+!
 !     USE pegrid
-! 
-! 
+!
+!
 !     IMPLICIT NONE
-! 
-!     INTEGER(iwp) ::  k               !< 
-!     INTEGER(iwp) ::  nxlc            !< 
-!     INTEGER(iwp) ::  nxlf            !< 
-!     INTEGER(iwp) ::  nxl_on_file     !< 
-!     INTEGER(iwp) ::  nxrc            !< 
-!     INTEGER(iwp) ::  nxrf            !< 
-!     INTEGER(iwp) ::  nxr_on_file     !< 
-!     INTEGER(iwp) ::  nync            !< 
-!     INTEGER(iwp) ::  nynf            !< 
-!     INTEGER(iwp) ::  nyn_on_file     !< 
-!     INTEGER(iwp) ::  nysc            !< 
-!     INTEGER(iwp) ::  nysf            !< 
-!     INTEGER(iwp) ::  nys_on_file     !< 
-! 
+!
+!     INTEGER(iwp) ::  k               !<
+!     INTEGER(iwp) ::  nxlc            !<
+!     INTEGER(iwp) ::  nxlf            !<
+!     INTEGER(iwp) ::  nxl_on_file     !<
+!     INTEGER(iwp) ::  nxrc            !<
+!     INTEGER(iwp) ::  nxrf            !<
+!     INTEGER(iwp) ::  nxr_on_file     !<
+!     INTEGER(iwp) ::  nync            !<
+!     INTEGER(iwp) ::  nynf            !<
+!     INTEGER(iwp) ::  nyn_on_file     !<
+!     INTEGER(iwp) ::  nysc            !<
+!     INTEGER(iwp) ::  nysf            !<
+!     INTEGER(iwp) ::  nys_on_file     !<
+!
 !     LOGICAL, INTENT(OUT)  :: found
-!     
+!
 !     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE  ::  tmp_3d_non_standard !< temporary array for storing 3D data with non standard dimensions
 ! !
-! !-- If temporary non-standard array for reading is already allocated, 
+! !-- If temporary non-standard array for reading is already allocated,
 ! !-- deallocate it.
 !     IF ( ALLOCATED( tmp_3d_non_standard ) )  DEALLOCATE( tmp_3d_non_standard )
-!     
+!
 !     found = .TRUE.
-! 
+!
 !     SELECT CASE ( restart_string(1:length) )
-! 
+!
 !        CASE ( 'ti_av' )
 !           IF ( .NOT. ALLOCATED( ti_av ) )  THEN
@@ -1769,5 +1988,5 @@
 !           ENDIF
 !           ti_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
-!     
+!
 !        CASE ( 'uu_av' )
 !           IF ( .NOT. ALLOCATED( uu_av ) )  THEN
@@ -1780,5 +1999,5 @@
 !           ENDIF
 !           uu_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
-!                    
+!
 !        CASE ( 'vv_av' )
 !           IF ( .NOT. ALLOCATED( vv_av ) )  THEN
@@ -1791,5 +2010,5 @@
 !           ENDIF
 !           vv_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
-!                    
+!
 !        CASE ( 'ww_av' )
 !           IF ( .NOT. ALLOCATED( ww_av ) )  THEN
@@ -1802,14 +2021,14 @@
 !           ENDIF
 !           ww_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
-!                         
-! 
+!
+!
 !        CASE DEFAULT
-! 
+!
 !           found = .FALSE.
-! 
+!
 !     END SELECT
-! 
+!
 !  END SUBROUTINE doq_rrd_local
- 
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -1861,10 +2080,10 @@
        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' )
@@ -1872,7 +2091,22 @@
     ENDIF
 
+    IF ( ALLOCATED( wspeed_av ) )  THEN
+       CALL wrd_write_string( 'wspeed_av' )
+       WRITE ( 14 )  wspeed_av
+    ENDIF
+
+    IF ( ALLOCATED( u_center_av ) )  THEN
+       CALL wrd_write_string( 'u_center_av' )
+       WRITE ( 14 )  u_center_av
+    ENDIF
+
+    IF ( ALLOCATED( v_center_av ) )  THEN
+       CALL wrd_write_string( 'v_center_av' )
+       WRITE ( 14 )  v_center_av
+    ENDIF
+
  END SUBROUTINE doq_wrd_local
- 
- 
+
+
 
  END MODULE diagnostic_output_quantities_mod
