Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 4038)
+++ palm/trunk/SOURCE/Makefile	(revision 4039)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# add dependency for diagnostic output quantities
+# 
+# 4017 2019-06-06 12:16:46Z schwenkel
 # add dependency for init_pegrid
 # 
@@ -770,5 +773,4 @@
 average_3d_data.o: \
 	cpulog_mod.o \
-	diagnostic_output_quantities_mod.o \
 	exchange_horiz_2d.o \
 	mod_kinds.o \
@@ -934,5 +936,4 @@
 	bulk_cloud_model_mod.o \
 	cpulog_mod.o \
-	diagnostic_output_quantities_mod.o \
 	land_surface_model_mod.o \
 	mod_kinds.o \
@@ -947,5 +948,4 @@
 	bulk_cloud_model_mod.o \
 	cpulog_mod.o \
-	diagnostic_output_quantities_mod.o \
 	land_surface_model_mod.o \
 	mod_kinds.o \
@@ -962,5 +962,6 @@
 diagnostic_output_quantities_mod.o: \
 	mod_kinds.o \
-	modules.o
+	modules.o \
+	surface_mod.o
 diffusion_s.o: \
 	mod_kinds.o \
@@ -1179,4 +1180,5 @@
 	bulk_cloud_model_mod.o \
 	chemistry_model_mod.o \
+	diagnostic_output_quantities_mod.o \
 	gust_mod.o \
 	indoor_model_mod.o \
@@ -1224,4 +1226,5 @@
 	biometeorology_mod.o \
 	chemistry_model_mod.o \
+	diagnostic_output_quantities_mod.o \
 	gust_mod.o \
 	indoor_model_mod.o \
@@ -1472,5 +1475,4 @@
 	bulk_cloud_model_mod.o \
 	cpulog_mod.o \
-	diagnostic_output_quantities_mod.o \
 	mod_kinds.o \
 	module_interface.o \
@@ -1714,5 +1716,4 @@
 write_restart_data_mod.o: \
 	date_and_time_mod.o \
-	diagnostic_output_quantities_mod.o \
 	mod_kinds.o \
 	model_1d_mod.o \
Index: palm/trunk/SOURCE/average_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/average_3d_data.f90	(revision 4038)
+++ palm/trunk/SOURCE/average_3d_data.f90	(revision 4039)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Modularize diagnostic output
+! 
+! 3994 2019-05-22 18:08:09Z suehring
 ! output of turbulence intensity added
 ! 
@@ -191,7 +194,4 @@
         ONLY:  cpu_log, log_point
 
-    USE diagnostic_output_quantities_mod,                                      &
-        ONLY:  ti_av
-
     USE indices,                                                               &
         ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
@@ -521,16 +521,4 @@
              ENDIF
 
-
-          CASE ( 'ti' )
-             IF ( ALLOCATED( ti_av ) ) THEN
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      DO  k = nzb, nzt+1
-                         ti_av(k,j,i) = ti_av(k,j,i) / REAL( average_count_3d, KIND=wp )
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
           CASE ( 'w' )
              IF ( ALLOCATED( w_av ) ) THEN
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 4038)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 4039)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Modularize diagnostic output
+! 
+! 4017 2019-06-06 12:16:46Z schwenkel
 ! output of turbulence intensity added
 ! 
@@ -3061,7 +3064,4 @@
              CONTINUE
 
-          CASE ( 'ti' )
-             unit = '1/s'
-
           CASE ( 'ghf*', 'lwp*', 'ol*', 'qsws*', 'r_a*',                       &
                  'shf*', 'ssws*', 't*', 'theta_2m*', 'tsurf*', 'us*',          &
Index: palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_2d.f90	(revision 4038)
+++ palm/trunk/SOURCE/data_output_2d.f90	(revision 4039)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! modularize diagnostic output
+! 
+! 3994 2019-05-22 18:08:09Z suehring
 ! output of turbulence intensity added
 ! 
@@ -324,7 +327,4 @@
     USE cpulog,                                                                &
         ONLY:  cpu_log, log_point
-
-    USE diagnostic_output_quantities_mod,                                      &
-        ONLY:  ti, ti_av
 
     USE indices,                                                               &
@@ -1171,5 +1171,5 @@
                    IF ( ibc_uv_b == 0 )  local_pf(:,:,nzb) = 0.0_wp
                 ENDIF
-
+                
              CASE ( 'us*_xy' )        ! 2d-array
                 IF ( av == 0 )  THEN
@@ -1266,16 +1266,4 @@
                 two_d = .TRUE.
                 level_z(nzb+1) = zu(nzb+1)
-                
-             CASE ( 'ti_xy', 'ti_xz', 'ti_yz' )
-                IF ( av == 0 )  THEN
-                   to_be_resorted => ti
-                ELSE
-                   IF ( .NOT. ALLOCATED( ti_av ) ) THEN
-                      ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
-                      ti_av = REAL( fill_value, KIND = wp )
-                   ENDIF
-                   to_be_resorted => ti_av
-                ENDIF
-                IF ( mode == 'xy' )  level_z = zu
 
              CASE ( 'w_xy', 'w_xz', 'w_yz' )
Index: palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_3d.f90	(revision 4038)
+++ palm/trunk/SOURCE/data_output_3d.f90	(revision 4039)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! modularize diagnostic output
+! 
+! 3994 2019-05-22 18:08:09Z suehring
 ! output of turbulence intensity added
 ! 
@@ -275,7 +278,4 @@
     USE cpulog,                                                                &
         ONLY:  log_point, cpu_log
-
-    USE diagnostic_output_quantities_mod,                                      &
-        ONLY:  ti, ti_av
 
 #if defined( __parallel )
@@ -703,15 +703,4 @@
                 ENDIF
                 to_be_resorted => vpt_av
-             ENDIF
-
-          CASE ( 'ti' )
-             IF ( av == 0 )  THEN
-                to_be_resorted => ti
-             ELSE
-                IF ( .NOT. ALLOCATED( ti_av ) ) THEN
-                   ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
-                   ti_av = REAL( fill_value, KIND = wp )
-                ENDIF
-                to_be_resorted => ti_av
              ENDIF
 
Index: palm/trunk/SOURCE/data_output_mask.f90
===================================================================
--- palm/trunk/SOURCE/data_output_mask.f90	(revision 4038)
+++ palm/trunk/SOURCE/data_output_mask.f90	(revision 4039)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Modularize diagnostic output
+! 
+! 3994 2019-05-22 18:08:09Z suehring
 ! output of turbulence intensity added
 ! 
@@ -180,9 +183,9 @@
                time_since_reference_point
 
+    USE diagnostic_output_quantities_mod,                                      & 
+        ONLY:  doq_output_mask
+               
     USE cpulog,                                                                &
         ONLY:  cpu_log, log_point
-
-    USE diagnostic_output_quantities_mod,                                      &
-        ONLY:  ti, ti_av
 
     USE indices,                                                               &
@@ -669,11 +672,4 @@
              ENDIF
 
-          CASE ( 'ti' )
-             IF ( av == 0 )  THEN
-                to_be_resorted => ti
-             ELSE
-                to_be_resorted => ti_av
-             ENDIF
-
           CASE ( 'w' )
              grid = 'w'
@@ -697,4 +693,7 @@
                                            local_pf )
              ENDIF
+!
+!--          Check for diagnostic quantities
+             CALL doq_output_mask( av, domask(mid,av,ivar), found, local_pf )
 !
 !--          SALSA quantities
Index: palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
===================================================================
--- palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4038)
+++ palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4039)
@@ -25,4 +25,12 @@
 ! -----------------
 ! $Id$
+! - Add output of uu, vv, ww to enable variance calculation according temporal
+!   EC method
+! - Allocate arrays only when they are required
+! - Formatting adjustment
+! - Rename subroutines
+! - Further modularization
+! 
+! 3998 2019-05-23 13:38:11Z suehring
 ! Bugfix in gathering all output strings 
 ! 
@@ -43,9 +51,6 @@
  MODULE diagnostic_output_quantities_mod
  
-
-!     USE arrays_3d,                                                             &
-!         ONLY:  dzw, e, heatflux_output_conversion, nc, nr, p, pt,              &
-!                precipitation_amount, prr, q, qc, ql, ql_c, ql_v, qr, s, tend,  &
-!                u, v, vpt, w, zu, zw, waterflux_output_conversion, hyrho, d_exner
+    USE arrays_3d,                                                             &
+        ONLY:  ddzu, u, v, w
 ! 
 !     USE averaging
@@ -62,8 +67,10 @@
 !     USE cpulog,                                                                &
 !         ONLY:  cpu_log, log_point
-! 
-!     USE indices,                                                               &
-!         ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,       &
-!                nzb, nzt, wall_flags_0
+
+   USE grid_variables,                                                         &
+        ONLY:  ddx, ddy
+! 
+    USE indices,                                                               &
+        ONLY:  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
 ! 
     USE kinds
@@ -103,10 +110,19 @@
 
     INTEGER(iwp) ::  timestep_number_at_prev_calc = 0  !< ...at previous diagnostic output calculation
-
-    LOGICAL ::  initialized_diagnostic_output_quantities = .FALSE.
-    LOGICAL ::  prepared_diagnostic_output_quantities = .FALSE.
+ 
+    LOGICAL ::  initialized_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is initialized
+    LOGICAL ::  prepared_diagnostic_output_quantities = .FALSE.    !< flag indicating whether output is p
 
     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
 
 
@@ -117,43 +133,679 @@
 !
 !-- Public variables
-    PUBLIC do_all,                                                                                 &
-           initialized_diagnostic_output_quantities,                                               &
-           prepared_diagnostic_output_quantities,                                                  &
-           ti, ti_av,                                                                              &
-           timestep_number_at_prev_calc
-!
-!-- Public routines
-    PUBLIC diagnostic_output_quantities_calculate,                                                 &
-           diagnostic_output_quantities_init,                                                      &
-           diagnostic_output_quantities_prepare
-
-
-    INTERFACE diagnostic_output_quantities_init
-       MODULE PROCEDURE diagnostic_output_quantities_init
-    END INTERFACE diagnostic_output_quantities_init
-
-    INTERFACE diagnostic_output_quantities_calculate
-       MODULE PROCEDURE diagnostic_output_quantities_calculate
-    END INTERFACE diagnostic_output_quantities_calculate
-
-    INTERFACE diagnostic_output_quantities_prepare
-       MODULE PROCEDURE diagnostic_output_quantities_prepare
-    END INTERFACE diagnostic_output_quantities_prepare
+    PUBLIC do_all,                                                             &
+           initialized_diagnostic_output_quantities,                           &
+           prepared_diagnostic_output_quantities,                              &
+           timestep_number_at_prev_calc,                                       &
+           ti_av,                                                              &
+           uu_av,                                                              &
+           vv_av,                                                              &
+           ww_av
+!                                                                              
+!-- Public routines                                                            
+    PUBLIC doq_3d_data_averaging,                                              &
+           doq_calculate,                                                      &
+           doq_check_data_output,                                              &
+           doq_define_netcdf_grid,                                             &
+           doq_output_2d,                                                      &
+           doq_output_3d,                                                      &
+           doq_output_mask,                                                    &
+           doq_init,                                                           &
+           doq_prepare,                                                        &
+           doq_wrd_local
+!          doq_rrd_local,                                                      &
+
+
+    INTERFACE doq_3d_data_averaging
+       MODULE PROCEDURE doq_3d_data_averaging
+    END INTERFACE doq_3d_data_averaging       
+
+    INTERFACE doq_calculate
+       MODULE PROCEDURE doq_calculate
+    END INTERFACE doq_calculate
+
+    INTERFACE doq_check_data_output
+       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
+    END INTERFACE doq_init
+
+    INTERFACE doq_prepare
+       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
+    END INTERFACE doq_wrd_local
 
 
  CONTAINS
-
+ 
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> ...
-!------------------------------------------------------------------------------!
- SUBROUTINE diagnostic_output_quantities_init
-
-
-    USE indices,                                                                                   &
-        ONLY:  nxl, nxr, nyn, nys, nzb, nzt
+!> Sum up and time-average diagnostic output quantities as well as allocate
+!> the array necessary for storing the average.
+!------------------------------------------------------------------------------!
+ SUBROUTINE doq_3d_data_averaging( mode, variable )
+
+    USE control_parameters,                                                    &
+        ONLY:  average_count_3d
+
+    CHARACTER (LEN=*) ::  mode    !< 
+    CHARACTER (LEN=*) :: variable !< 
+
+    INTEGER(iwp) ::  i !<
+    INTEGER(iwp) ::  j !<
+    INTEGER(iwp) ::  k !<
+
+    IF ( mode == 'allocate' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+
+          CASE ( 'ti' )
+             IF ( .NOT. ALLOCATED( ti_av ) )  THEN
+                ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+             ENDIF
+             ti_av = 0.0_wp
+       
+          CASE ( 'uu' )
+             IF ( .NOT. ALLOCATED( uu_av ) )  THEN
+                ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+             ENDIF
+             uu_av = 0.0_wp
+                
+          CASE ( 'vv' )
+             IF ( .NOT. ALLOCATED( vv_av ) )  THEN
+                ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+             ENDIF
+             vv_av = 0.0_wp
+                
+          CASE ( 'ww' )
+             IF ( .NOT. ALLOCATED( ww_av ) )  THEN
+                ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+             ENDIF
+             ww_av = 0.0_wp
+
+          CASE DEFAULT
+             CONTINUE
+
+       END SELECT
+
+    ELSEIF ( mode == 'sum' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+ 
+          CASE ( 'ti' )
+             IF ( ALLOCATED( ti_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         ti_av(k,j,i) = ti_av(k,j,i) + ti(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             
+          CASE ( 'uu' )
+             IF ( ALLOCATED( uu_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         uu_av(k,j,i) = uu_av(k,j,i) + uu(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             
+          CASE ( 'vv' )
+             IF ( ALLOCATED( vv_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         vv_av(k,j,i) = vv_av(k,j,i) + vv(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             
+          CASE ( 'ww' )
+             IF ( ALLOCATED( ww_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         ww_av(k,j,i) = ww_av(k,j,i) + ww(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          CASE DEFAULT
+             CONTINUE
+
+       END SELECT
+
+    ELSEIF ( mode == 'average' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+
+          CASE ( 'ti' )
+             IF ( ALLOCATED( ti_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         ti_av(k,j,i) = ti_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+       
+          CASE ( 'uu' )
+             IF ( ALLOCATED( uu_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         uu_av(k,j,i) = uu_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             
+          CASE ( 'vv' )
+             IF ( ALLOCATED( vv_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         vv_av(k,j,i) = vv_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             
+          CASE ( 'ww' )
+             IF ( ALLOCATED( ww_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      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
+
+       END SELECT
+
+    ENDIF
+
+
+ END SUBROUTINE doq_3d_data_averaging 
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check data output for diagnostic output
+!------------------------------------------------------------------------------!
+ SUBROUTINE doq_check_data_output( var, unit )
 
     IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  unit  !< 
+    CHARACTER (LEN=*) ::  var   !<
+
+    SELECT CASE ( TRIM( var ) )
+
+       CASE ( 'ti' )
+          unit = '1/s'
+             
+       CASE ( 'uu' )
+          unit = 'm2/s2'
+             
+       CASE ( 'vv' )
+          unit = 'm2/s2'
+             
+       CASE ( 'ww' )
+          unit = 'm2/s2'
+             
+       CASE DEFAULT
+          unit = 'illegal'
+
+    END SELECT
+
+
+ END SUBROUTINE doq_check_data_output
+ 
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> 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      !< 
+
+    found  = .TRUE.
+    
+    SELECT CASE ( TRIM( variable ) )
+!
+!--    s grid
+       CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz'  )
+
+          grid_x = 'x'
+          grid_y = 'y'
+          grid_z = 'zu'    
+!
+!--    u grid
+       CASE ( 'uu', 'uu_xy', 'uu_xz', 'uu_yz' )
+
+          grid_x = 'xu'
+          grid_y = 'y'
+          grid_z = 'zu'
+!
+!--    v grid
+       CASE ( 'vv', 'vv_xy', 'vv_xz', 'vv_yz'  )
+
+          grid_x = 'x'
+          grid_y = 'yv'
+          grid_z = 'zu'
+!
+!--    w grid
+       CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz'  )
+
+          grid_x = 'x'
+          grid_y = 'y'
+          grid_z = 'zw'
+
+       CASE DEFAULT
+          found  = .FALSE.
+          grid_x = 'none'
+          grid_y = 'none'
+          grid_z = 'none'
+
+    END SELECT
+
+
+ END SUBROUTINE doq_define_netcdf_grid
+ 
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine defining 2D output variables
+!------------------------------------------------------------------------------!
+ SUBROUTINE doq_output_2d( av, variable, found, grid,                          &
+                           mode, local_pf, two_d, nzb_do, nzt_do, fill_value )
+
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  grid     !< 
+    CHARACTER (LEN=*) ::  mode     !< 
+    CHARACTER (LEN=*) ::  variable !< 
+
+    INTEGER(iwp) ::  av       !< value indicating averaged or non-averaged output
+    INTEGER(iwp) ::  flag_nr  !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
+    INTEGER(iwp) ::  i        !< grid index x-direction
+    INTEGER(iwp) ::  j        !< grid index y-direction
+    INTEGER(iwp) ::  k        !< grid index z-direction
+    INTEGER(iwp) ::  nzb_do   !< 
+    INTEGER(iwp) ::  nzt_do   !< 
+
+    LOGICAL ::  found             !< true if variable is in list
+    LOGICAL ::  resorted          !< true if array is resorted
+    LOGICAL ::  two_d             !< flag parameter that indicates 2D variables (horizontal cross sections)
+
+    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(:,:,:), POINTER ::                 to_be_resorted  !< points to array which needs to be resorted for output
+    
+    flag_nr  = 0
+    found    = .TRUE.
+    resorted = .FALSE.
+    two_d    = .FALSE.
+
+    SELECT CASE ( TRIM( variable ) )
+
+       CASE ( 'ti_xy', 'ti_xz', 'ti_yz' )
+           IF ( av == 0 )  THEN
+              to_be_resorted => ti
+           ELSE
+              IF ( .NOT. ALLOCATED( ti_av ) ) THEN
+                 ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                 ti_av = REAL( fill_value, KIND = wp )
+              ENDIF
+              to_be_resorted => ti_av
+           ENDIF
+           flag_nr = 0
+           
+           IF ( mode == 'xy' )  grid = 'zu'
+    
+       CASE ( 'uu_xy', 'uu_xz', 'uu_yz' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => uu
+          ELSE
+             IF ( .NOT. ALLOCATED( uu_av ) ) THEN
+                ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                uu_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => uu_av
+          ENDIF
+          flag_nr = 1
+          
+          IF ( mode == 'xy' )  grid = 'zu'
+          
+       CASE ( 'vv_xy', 'vv_xz', 'vv_yz' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => vv
+          ELSE
+             IF ( .NOT. ALLOCATED( vv_av ) ) THEN
+                ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                vv_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => vv_av
+          ENDIF
+          flag_nr = 2
+          
+          IF ( mode == 'xy' )  grid = 'zu'
+                
+       CASE ( 'ww_xy', 'ww_xz', 'ww_yz' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => ww
+          ELSE
+             IF ( .NOT. ALLOCATED( ww_av ) ) THEN
+                ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                ww_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => ww_av
+          ENDIF
+          flag_nr = 3
+          
+          IF ( mode == 'xy' )  grid = 'zw'
+
+       CASE DEFAULT
+          found = .FALSE.
+          grid  = 'none'
+
+    END SELECT
+    
+    IF ( found  .AND.  .NOT. resorted )  THEN     
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_do, nzt_do
+                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
+                                         REAL( fill_value, KIND = wp ),        &
+                                         BTEST( wall_flags_0(k,j,i), flag_nr ) ) 
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDIF
+ 
+ END SUBROUTINE doq_output_2d
+ 
+ 
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine defining 3D output variables
+!------------------------------------------------------------------------------!
+ SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do,  &
+                           nzt_do )
+ 
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  variable !< 
+
+    INTEGER(iwp) ::  av       !< index indicating averaged or instantaneous output
+    INTEGER(iwp) ::  flag_nr  !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
+    INTEGER(iwp) ::  i        !< index variable along x-direction
+    INTEGER(iwp) ::  j        !< index variable along y-direction
+    INTEGER(iwp) ::  k        !< index variable along z-direction
+    INTEGER(iwp) ::  nzb_do   !< lower limit of the data output (usually 0)
+    INTEGER(iwp) ::  nzt_do   !< vertical upper limit of the data output (usually nz_do3d)
+
+    LOGICAL ::  found             !< true if variable is in list
+    LOGICAL ::  resorted          !< true if array is resorted
+
+    REAL(wp) ::  fill_value       !< value for the _FillValue attribute
+
+    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
+
+    flag_nr  = 0
+    found    = .TRUE.
+    resorted = .FALSE.
+    
+    SELECT CASE ( TRIM( variable ) )
+
+       CASE ( 'ti' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => ti
+          ELSE
+             IF ( .NOT. ALLOCATED( ti_av ) ) THEN
+                ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                ti_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => ti_av
+          ENDIF
+          flag_nr = 0
+    
+       CASE ( 'uu' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => uu
+          ELSE
+             IF ( .NOT. ALLOCATED( uu_av ) ) THEN
+                ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                uu_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => uu_av
+          ENDIF
+          flag_nr = 1
+             
+       CASE ( 'vv' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => vv
+          ELSE
+             IF ( .NOT. ALLOCATED( vv_av ) ) THEN
+                ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                vv_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => vv_av
+          ENDIF
+          flag_nr = 2
+             
+       CASE ( 'ww' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => ww
+          ELSE
+             IF ( .NOT. ALLOCATED( ww_av ) ) THEN
+                ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                ww_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => ww_av
+          ENDIF
+          flag_nr = 3
+
+       CASE DEFAULT
+          found = .FALSE.
+
+    END SELECT
+    
+    IF ( found  .AND.  .NOT. resorted )  THEN     
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_do, nzt_do
+                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
+                                         REAL( fill_value, KIND = wp ),        &
+                                         BTEST( wall_flags_0(k,j,i), flag_nr ) ) 
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDIF
+
+ END SUBROUTINE doq_output_3d
+ 
+! Description:
+! ------------
+!> Resorts the user-defined output quantity with indices (k,j,i) to a
+!> temporary array with indices (i,j,k) for masked data output.
+!------------------------------------------------------------------------------!
+ SUBROUTINE doq_output_mask( av, variable, found, local_pf )
+ 
+    USE control_parameters
+        
+    USE indices
+    
+    USE surface_mod,                                                           &
+        ONLY:  get_topography_top_index_ji
+ 
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  variable  !< 
+    CHARACTER (LEN=5) ::  grid      !< flag to distinquish between staggered grids
+
+    INTEGER(iwp) ::  av           !< index indicating averaged or instantaneous output
+    INTEGER(iwp) ::  flag_nr      !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
+    INTEGER(iwp) ::  i            !< index variable along x-direction
+    INTEGER(iwp) ::  j            !< index variable along y-direction
+    INTEGER(iwp) ::  k            !< index variable along z-direction
+    INTEGER(iwp) ::  topo_top_ind !< k index of highest horizontal surface
+
+    LOGICAL ::  found             !< true if variable is in list
+    LOGICAL ::  resorted          !< true if array is resorted
+
+    REAL(wp),                                                                  &
+       DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
+          local_pf   !< 
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to array which needs to be resorted for output
+
+    flag_nr  = 0
+    found    = .TRUE.
+    resorted = .FALSE.
+    grid     = 's'
+
+    SELECT CASE ( TRIM( variable ) )
+
+       CASE ( 'ti' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => ti
+          ELSE
+             to_be_resorted => ti_av
+          ENDIF
+          grid = 's'
+          flag_nr = 0
+    
+       CASE ( 'uu' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => uu
+          ELSE
+             to_be_resorted => uu_av
+          ENDIF
+          grid = 'u'
+          flag_nr = 1
+    
+       CASE ( 'vv' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => vv
+          ELSE
+             to_be_resorted => vv_av
+          ENDIF
+          grid = 'v'
+          flag_nr = 2
+    
+       CASE ( 'ww' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => ww
+          ELSE
+             to_be_resorted => ww_av
+          ENDIF
+          grid = 'w'
+          flag_nr = 3
+
+
+       CASE DEFAULT
+          found = .FALSE.
+
+    END SELECT
+    
+    IF ( found  .AND.  .NOT. resorted )  THEN
+       IF ( .NOT. mask_surface(mid) )  THEN
+!
+!--       Default masked output
+          DO  i = 1, mask_size_l(mid,1)
+             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))
+                ENDDO
+             ENDDO
+          ENDDO
+
+       ELSE
+!
+!--       Terrain-following masked output
+          DO  i = 1, mask_size_l(mid,1)
+             DO  j = 1, mask_size_l(mid,2)
+!
+!--             Get k index of highest horizontal surface
+                topo_top_ind = get_topography_top_index_ji( mask_j(mid,j),     &
+                                                            mask_i(mid,i),     &
+                                                            grid )
+!
+!--             Save output array
+                DO  k = 1, mask_size_l(mid,3)
+                   local_pf(i,j,k) = to_be_resorted(                           &
+                                          MIN( topo_top_ind+mask_k(mid,k),     &
+                                               nzt+1 ),                        &
+                                          mask_j(mid,j),                       &
+                                          mask_i(mid,i)                     )
+                ENDDO
+             ENDDO
+          ENDDO
+
+       ENDIF
+    ENDIF
+    
+ END SUBROUTINE doq_output_mask
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Allocate required arrays
+!------------------------------------------------------------------------------!
+ SUBROUTINE doq_init
+
+    IMPLICIT NONE
+    
+    INTEGER(iwp) ::  ivar   !< loop index over all 2d/3d/mask output quantities
 !
 !-- Next line is to avoid compiler warnings about unused variables
@@ -161,13 +813,49 @@
 
     initialized_diagnostic_output_quantities = .FALSE.
-
-    IF ( .NOT. ALLOCATED( ti ) )  THEN
-       ALLOCATE( ti(nzb:nzt+1,nys:nyn,nxl:nxr) )
-       ti = 0.0_wp
-    ENDIF
+    
+    ivar = 1
+    
+    DO  WHILE ( ivar <= SIZE( do_all ) ) 
+
+       SELECT CASE ( TRIM( do_all(ivar) ) )
+!
+!--       Allocate array for 'turbulence intensity'
+          CASE ( 'ti' )
+             IF ( .NOT. ALLOCATED( ti ) )  THEN
+                ALLOCATE( ti(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                ti = 0.0_wp
+             ENDIF
+!
+!--       Allocate array for uu
+          CASE ( 'uu' )
+             IF ( .NOT. ALLOCATED( uu ) )  THEN
+                ALLOCATE( uu(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                uu = 0.0_wp
+             ENDIF
+
+!
+!--       Allocate array for vv
+          CASE ( 'vv' )
+             IF ( .NOT. ALLOCATED( vv ) )  THEN
+                ALLOCATE( vv(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                vv = 0.0_wp
+             ENDIF
+
+!
+!--       Allocate array for ww
+          CASE ( 'ww' )
+             IF ( .NOT. ALLOCATED( ww ) )  THEN
+                ALLOCATE( ww(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                ww = 0.0_wp
+             ENDIF
+
+       END SELECT
+
+       ivar = ivar + 1
+    ENDDO
 
     initialized_diagnostic_output_quantities = .TRUE.
-
- END SUBROUTINE diagnostic_output_quantities_init
+    
+ END SUBROUTINE doq_init
 
 
@@ -175,22 +863,12 @@
 ! Description:
 ! ------------
-!> ...
+!> Calculate diagnostic quantities
 !--------------------------------------------------------------------------------------------------!
- SUBROUTINE diagnostic_output_quantities_calculate
-
-
-    USE arrays_3d,                                                                                 &
-        ONLY:  ddzu, u, v, w
-
-    USE grid_variables,                                                                            &
-        ONLY:  ddx, ddy
-
-    USE indices,                                                                                   &
-        ONLY:  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
+ SUBROUTINE doq_calculate
 
     IMPLICIT NONE
 
     INTEGER(iwp) ::  i, j, k    !< grid loop index in x-, y-, z-direction
-    INTEGER(iwp) ::  ivar_all   !< loop index over all 2d/3d/mask output quantities
+    INTEGER(iwp) ::  ivar       !< loop index over all 2d/3d/mask output quantities
 
 
@@ -198,18 +876,17 @@
 !
 !-- Preparatory steps and initialization of output arrays
-    IF ( .NOT.  prepared_diagnostic_output_quantities )     CALL diagnostic_output_quantities_prepare
-    IF ( .NOT.  initialized_diagnostic_output_quantities )  CALL diagnostic_output_quantities_init
-!
-!-- Save timestep number to check in time_integration if diagnostic_output_quantities_calculate 
+    IF ( .NOT.  prepared_diagnostic_output_quantities )     CALL doq_prepare
+    IF ( .NOT.  initialized_diagnostic_output_quantities )  CALL doq_init
+!
+!-- 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.
     timestep_number_at_prev_calc = current_timestep_number
 
-
-    ivar_all = 1
-
-    DO  WHILE ( do_all(ivar_all)(1:1) /= ' ' )
-
-       SELECT CASE ( TRIM( do_all(ivar_all) ) )
+    ivar = 1
+
+    DO  WHILE ( ivar <= SIZE( do_all ) ) 
+
+       SELECT CASE ( TRIM( do_all(ivar) ) )
 !
 !--       Calculate 'turbulence intensity' from rot[(u,v,w)] at scalar grid point
@@ -218,31 +895,63 @@
                 DO  j = nys, nyn
                    DO  k = nzb+1, nzt
-
-                      ti(k,j,i) = 0.25_wp * SQRT(                                                  &
-                        (   (   w(k,j+1,i) + w(k-1,j+1,i)                                          &
-                              - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy                                  &
-                          - (   v(k+1,j,i) + v(k+1,j+1,i)                                          &
-                              - v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2                         &
-                      + (   (   u(k+1,j,i) + u(k+1,j,i+1)                                          &
-                              - u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k)                              &
-                          - (   w(k,j,i+1) + w(k-1,j,i+1)                                          &
-                              - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx     )**2                         &
-                      + (   (   v(k,j,i+1) + v(k,j+1,i+1)                                          &
-                              - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx                                  &
-                          - (   u(k,j+1,i) + u(k,j+1,i+1)                                          &
-                              - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy     )**2  )                      &
-                                          * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0) )
+                      ti(k,j,i) = 0.25_wp * SQRT(                              &
+                        (   (   w(k,j+1,i) + w(k-1,j+1,i)                      &
+                              - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy              &
+                          - (   v(k+1,j,i) + v(k+1,j+1,i)                      &
+                              - v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2     &
+                      + (   (   u(k+1,j,i) + u(k+1,j,i+1)                      &
+                              - u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k)          &
+                          - (   w(k,j,i+1) + w(k-1,j,i+1)                      &
+                              - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx     )**2     &
+                      + (   (   v(k,j,i+1) + v(k,j+1,i+1)                      &
+                              - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx              &
+                          - (   u(k,j+1,i) + u(k,j+1,i+1)                      &
+                              - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy     )**2  )  &
+                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0) )
                    ENDDO
                 ENDDO
+             ENDDO           
+!
+!--       uu
+          CASE ( 'uu' )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   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_0(k,j,i), 1) )
+                   ENDDO
+                ENDDO
              ENDDO
-
-          END SELECT
-
-          ivar_all = ivar_all + 1
+!
+!--       vv
+          CASE ( 'vv' )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   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_0(k,j,i), 2) )
+                   ENDDO
+                ENDDO
+             ENDDO
+!
+!--       ww
+          CASE ( 'ww' )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   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_0(k,j,i), 3) )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+       END SELECT
+
+       ivar = ivar + 1
     ENDDO
 
 !     CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' )
 
- END SUBROUTINE diagnostic_output_quantities_calculate
+ END SUBROUTINE doq_calculate
 
 
@@ -252,8 +961,8 @@
 !> ...
 !------------------------------------------------------------------------------!
- SUBROUTINE diagnostic_output_quantities_prepare
-
-
-    USE control_parameters,                                                                        &
+ SUBROUTINE doq_prepare
+
+
+    USE control_parameters,                                                    &
         ONLY:  do2d, do3d, domask, masks, mid
 
@@ -268,5 +977,4 @@
     INTEGER(iwp) ::  l          !< index for cutting string
 
-
     prepared_diagnostic_output_quantities = .FALSE.
 
@@ -280,5 +988,7 @@
        do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
 !
-!--    Gather 2d output quantity names, check for double occurrence of output quantity
+!--    Gather 2d output quantity names.
+!--    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
@@ -293,8 +1003,7 @@
        ivar = 1
 !
-!--    Gather 3d output quantity names, check for double occurrence of output quantity
+!--    Gather 3d output quantity names
        DO  WHILE ( do3d(av,ivar)(1:1) /= ' ' )
           do_all(ivar_all) = do3d(av,ivar)
-
           ivar = ivar + 1
           ivar_all = ivar_all + 1
@@ -303,8 +1012,11 @@
        ivar = 1
 !
-!--    Gather masked output quantity names, check for double occurrence of output quantity
+!--    Gather masked output quantity names. Also check for double output
+!--    e.g. by different masks.
        DO  mid = 1, masks
           DO  WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
-             do_all(ivar_all) = domask(mid,av,ivar)
+             IF ( .NOT.  ANY( do_all == domask(mid,av,ivar) ) )  THEN
+                do_all(ivar_all) = do2d_var(av,ivar)
+             ENDIF
 
              ivar = ivar + 1
@@ -318,5 +1030,142 @@
     prepared_diagnostic_output_quantities = .TRUE.
 
- END SUBROUTINE diagnostic_output_quantities_prepare
+ END SUBROUTINE doq_prepare
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine reads local (subdomain) restart data
+!> Note: With the current structure reading of non-standard array is not 
+!> possible
+!------------------------------------------------------------------------------!
+!  SUBROUTINE doq_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,             &
+!                            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     !< 
+! 
+!     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, 
+! !-- 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
+!              ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+!           ENDIF
+!           IF ( k == 1 )  THEN
+!              ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,  &
+!                                            nxl_on_file:nxr_on_file) )
+!              READ ( 13 )  tmp_3d_non_standard
+!           ENDIF
+!           ti_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
+!     
+!        CASE ( 'uu_av' )
+!           IF ( .NOT. ALLOCATED( uu_av ) )  THEN
+!              ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+!           ENDIF
+!           IF ( k == 1 )  THEN
+!              ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,  &
+!                                            nxl_on_file:nxr_on_file) )
+!              READ ( 13 )  tmp_3d_non_standard
+!           ENDIF
+!           uu_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
+!                    
+!        CASE ( 'vv_av' )
+!           IF ( .NOT. ALLOCATED( vv_av ) )  THEN
+!              ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+!           ENDIF
+!           IF ( k == 1 )  THEN
+!              ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,  &
+!                                            nxl_on_file:nxr_on_file) )
+!              READ ( 13 )  tmp_3d_non_standard
+!           ENDIF
+!           vv_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
+!                    
+!        CASE ( 'ww_av' )
+!           IF ( .NOT. ALLOCATED( ww_av ) )  THEN
+!              ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+!           ENDIF
+!           IF ( k == 1 )  THEN
+!              ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,  &
+!                                            nxl_on_file:nxr_on_file) )
+!              READ ( 13 )  tmp_3d_non_standard
+!           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:
+! ------------
+!> Subroutine writes local (subdomain) restart data
+!------------------------------------------------------------------------------!
+ SUBROUTINE doq_wrd_local
+
+
+    IMPLICIT NONE
+
+    IF ( ALLOCATED( ti_av ) )  THEN
+       CALL wrd_write_string( 'ti_av' ) 
+       WRITE ( 14 )  ti_av
+    ENDIF
+    
+    IF ( ALLOCATED( uu_av ) )  THEN
+       CALL wrd_write_string( 'uu_av' ) 
+       WRITE ( 14 )  uu_av
+    ENDIF
+       
+    IF ( ALLOCATED( vv_av ) )  THEN
+       CALL wrd_write_string( 'vv_av' ) 
+       WRITE ( 14 )  vv_av
+    ENDIF
+       
+    IF ( ALLOCATED( ww_av ) )  THEN
+       CALL wrd_write_string( 'ww_av' ) 
+       WRITE ( 14 )  ww_av
+    ENDIF
+
+
+ END SUBROUTINE doq_wrd_local
+ 
+ 
 
  END MODULE diagnostic_output_quantities_mod
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 4038)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 4039)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Correct conversion to kinematic scalar fluxes in case of pw-scheme and 
+! statistic regions
+! 
+! 3828 2019-03-27 19:36:23Z raasch
 ! unused variables removed
 ! 
@@ -1497,4 +1501,5 @@
                                     vpt(k+1,j,i) - hom(k+1,1,44,sr) )
                       sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
+                                               rho_air_zw(k) *                 &
                                                heatflux_output_conversion(k) * &
                                                           rmask(j,i,sr) * flag
@@ -1509,4 +1514,5 @@
                               hom(k+1,1,42,sr) )
                          sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
+                                             rho_air_zw(k) *                   &
                                              waterflux_output_conversion(k) *  &
                                                              rmask(j,i,sr)  *  &
@@ -1538,4 +1544,5 @@
                                           vpt(k+1,j,i) - hom(k+1,1,44,sr) )
                          sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
+                                              rho_air_zw(k) *                  &
                                               heatflux_output_conversion(k) *  &
                                                              rmask(j,i,sr)  *  &
@@ -1685,4 +1692,5 @@
                            ( pt(k,j,i)   - hom(k,1,4,sr) +                     &
                              pt(k+1,j,i) - hom(k+1,1,4,sr) )                   &
+                           * rho_air_zw(k)                                     &
                            * heatflux_output_conversion(k)                     &
                            * w(k,j,i) * rmask(j,i,sr) * flag
@@ -1691,4 +1699,5 @@
                                       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) *     &
+                                       rho_air_zw(k) *                         &
                                        waterflux_output_conversion(k) *        &
                                        rmask(j,i,sr) * flag
Index: palm/trunk/SOURCE/module_interface.f90
===================================================================
--- palm/trunk/SOURCE/module_interface.f90	(revision 4038)
+++ palm/trunk/SOURCE/module_interface.f90	(revision 4039)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Introduce diagnostic output
+! 
+! 4028 2019-06-13 12:21:37Z schwenkel
 ! Further modularization of particle code components
 ! 
@@ -210,4 +213,13 @@
               chem_rrd_local,                                                  &
               chem_wrd_local
+     
+    USE diagnostic_output_quantities_mod,                                      & 
+        ONLY:  doq_3d_data_averaging,                                          &
+               doq_check_data_output,                                          &
+               doq_define_netcdf_grid,                                         &
+               doq_output_2d,                                                  &
+               doq_output_3d,                                                  &
+               doq_wrd_local
+!                doq_rrd_local,                                                  &
 
     USE flight_mod,                                                            &
@@ -751,4 +763,8 @@
        CALL chem_check_data_output( variable, unit, i, ilen, k )
     ENDIF
+    
+    IF ( unit == 'illegal' )  THEN
+       CALL doq_check_data_output( variable, unit )
+    ENDIF
 
     IF ( unit == 'illegal'  .AND.  gust_module_enabled  )  THEN
@@ -1179,4 +1195,5 @@
     IF ( bulk_cloud_model    )  CALL bcm_3d_data_averaging( mode, variable )
     IF ( air_chemistry       )  CALL chem_3d_data_averaging( mode, variable )
+    CALL doq_3d_data_averaging( mode, variable )
     IF ( gust_module_enabled )  CALL gust_3d_data_averaging( mode, variable )
     IF ( land_surface        )  CALL lsm_3d_data_averaging( mode, variable )
@@ -1234,4 +1251,10 @@
             )
     ENDIF
+    
+    IF ( .NOT. found )  THEN
+       CALL doq_output_2d(                                                     &
+               av, variable, found, grid, mode, local_pf, two_d,               &
+               nzb_do, nzt_do, fill_value )
+    ENDIF
 
     IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
@@ -1314,5 +1337,10 @@
        resorted = .TRUE.
     ENDIF
-
+    
+    IF ( .NOT. found )  THEN
+       CALL doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
+       resorted = .TRUE.
+    ENDIF
+    
     IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
        CALL gust_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
@@ -1466,15 +1494,19 @@
     LOGICAL,      INTENT(INOUT) ::  found        !< flag if variable was found
 
-    REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp), INTENT(OUT) :: tmp_2d   !<
-    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp), INTENT(OUT) :: tmp_3d   !<
+    REAL(wp), &
+       DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp), &
+       INTENT(OUT) :: tmp_2d   !<
+    REAL(wp), &
+       DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp), &
+       INTENT(OUT) :: tmp_3d   !<
 
 
     IF ( debug_output )  CALL debug_message( 'module-specific read local restart data', 'start' )
 
-    IF ( .NOT. found ) CALL bio_rrd_local(                                     &
+    IF ( .NOT. found )  CALL bio_rrd_local(                                    &
                                found                                           &
                             )
 
-    IF ( .NOT. found ) CALL bcm_rrd_local(                                     &
+    IF ( .NOT. found )  CALL bcm_rrd_local(                                    &
                                map_index,                                      &
                                nxlf, nxlc, nxl_on_file,                        &
@@ -1485,5 +1517,5 @@
                             ) ! ToDo: change interface to pass variable
 
-    IF ( .NOT. found ) CALL chem_rrd_local(                                    &
+    IF ( .NOT. found )  CALL chem_rrd_local(                                   &
                                map_index,                                      &
                                nxlf, nxlc, nxl_on_file,                        &
@@ -1493,6 +1525,15 @@
                                tmp_3d, found                                   &
                             ) ! ToDo: change interface to pass variable
-
-    IF ( .NOT. found ) CALL gust_rrd_local(                                    &
+                            
+!     IF ( .NOT. found )  CALL doq_rrd_local(                                    &
+!                                map_index,                                      &
+!                                nxlf, nxlc, nxl_on_file,                        &
+!                                nxrf, nxrc, nxr_on_file,                        &
+!                                nynf, nync, nyn_on_file,                        &
+!                                nysf, nysc, nys_on_file,                        &
+!                                tmp_3d_non_standard, found                      &
+!                             ) ! ToDo: change interface to pass variable CALL doq_wrd_local
+
+    IF ( .NOT. found )  CALL gust_rrd_local(                                   &
                                map_index,                                      &
                                nxlf, nxlc, nxl_on_file,                        &
@@ -1503,5 +1544,5 @@
                             ) ! ToDo: change interface to pass variable
 
-    IF ( .NOT. found ) CALL lpm_rrd_local(                                     &
+    IF ( .NOT. found )  CALL lpm_rrd_local(                                    &
                                map_index,                                      &
                                nxlf, nxlc, nxl_on_file,                        &
@@ -1512,5 +1553,5 @@
                             ) ! ToDo: change interface to pass variable
 
-    IF ( .NOT. found ) CALL lsm_rrd_local(                                     &
+    IF ( .NOT. found )  CALL lsm_rrd_local(                                    &
                                map_index,                                      &
                                nxlf, nxlc, nxl_on_file,                        &
@@ -1521,5 +1562,5 @@
                             ) ! ToDo: change interface to pass variable
 
-    IF ( .NOT. found ) CALL ocean_rrd_local(                                   &
+    IF ( .NOT. found )  CALL ocean_rrd_local(                                  &
                                map_index,                                      &
                                nxlf, nxlc, nxl_on_file,                        &
@@ -1530,5 +1571,5 @@
                             ) ! ToDo: change interface to pass variable
 
-    IF ( .NOT. found ) CALL radiation_rrd_local(                               &
+    IF ( .NOT. found )  CALL radiation_rrd_local(                              &
                                map_index,                                      &
                                nxlf, nxlc, nxl_on_file,                        &
@@ -1539,5 +1580,5 @@
                             ) ! ToDo: change interface to pass variable
 
-    IF ( .NOT. found ) CALL salsa_rrd_local(                                   &
+    IF ( .NOT. found )  CALL salsa_rrd_local(                                  &
                                map_index,                                      &
                                nxlf, nxlc, nxl_on_file,                        &
@@ -1548,5 +1589,5 @@
                             ) ! ToDo: change interface to pass variable
 
-    IF ( .NOT. found ) CALL usm_rrd_local(                                     &
+    IF ( .NOT. found )  CALL usm_rrd_local(                                    &
                                map_index,                                      &
                                nxlf, nxlc, nxl_on_file,                        &
@@ -1558,7 +1599,7 @@
 !
 !-- Surface data do not need overlap data, so do not pass these information.
-    IF ( .NOT. found ) CALL surface_data_output_rrd_local( found )
+    IF ( .NOT. found )  CALL surface_data_output_rrd_local( found )
                            
-    IF ( .NOT. found ) CALL user_rrd_local(                                    &
+    IF ( .NOT. found )  CALL user_rrd_local(                                   &
                                map_index,                                      &
                                nxlf, nxlc, nxl_on_file,                        &
@@ -1588,4 +1629,5 @@
     IF ( bulk_cloud_model )     CALL bcm_wrd_local
     IF ( air_chemistry )        CALL chem_wrd_local
+    CALL doq_wrd_local
     IF ( gust_module_enabled )  CALL gust_wrd_local
     IF ( particle_advection )   CALL lpm_wrd_local    
Index: palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 4038)
+++ palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 4039)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Rename subroutines in module for diagnostic quantities
+! 
+! 4029 2019-06-14 14:04:35Z raasch
 ! netcdf variable NF90_NOFILL is used as argument instead of "1" in calls to NF90_DEF_VAR_FILL
 ! 
@@ -698,4 +701,7 @@
                skip_time_do3d, topography, num_leg, num_var_fl,                &
                urban_surface
+
+    USE diagnostic_output_quantities_mod,                                      & 
+        ONLY:  doq_define_netcdf_grid           
 
     USE grid_variables,                                                        &
@@ -1124,5 +1130,5 @@
                 CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',                &
                        'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
-                       's', 'theta', 'thetal', 'thetav', 'ti' )
+                       's', 'theta', 'thetal', 'thetav' )
 
                    grid_x = 'x'
@@ -1169,4 +1175,8 @@
                    ENDIF
 
+                   IF ( .NOT. found )                                          &
+                      CALL doq_define_netcdf_grid( domask(mid,av,i), found,    &
+                                                   grid_x, grid_y, grid_z )
+                   
                    IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
                       CALL gust_define_netcdf_grid( domask(mid,av,i), found,   &
@@ -1870,5 +1880,5 @@
                 CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',   &
                        'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
-                       's', 'theta', 'thetal', 'thetav', 'ti' )
+                       's', 'theta', 'thetal', 'thetav' )
 
                    grid_x = 'x'
@@ -1974,4 +1984,8 @@
                                                     grid_y, grid_z )
                    ENDIF
+                   
+                   IF ( .NOT. found )                                          &
+                      CALL doq_define_netcdf_grid( do3d(av,i), found, grid_x,  &
+                                                   grid_y, grid_z        )
                                                  
                    IF ( .NOT. found )  THEN
@@ -2858,5 +2872,5 @@
                              'ql_vp_xy', 'qr_xy', 'qv_xy',                     &
                              's_xy',                                           &
-                             'theta_xy', 'thetal_xy', 'thetav_xy', 'ti_xy' )
+                             'theta_xy', 'thetal_xy', 'thetav_xy' )
 
                          grid_x = 'x'
@@ -2944,6 +2958,9 @@
                                                           grid_z )
                          ENDIF
-
-
+                         
+                         IF ( .NOT. found )                                    &
+                            CALL doq_define_netcdf_grid(                       &
+                                                    do2d(av,i), found, grid_x, &
+                                                    grid_y, grid_z              )
 !
 !--                      Check for user-defined quantities
@@ -3762,5 +3779,5 @@
                           'ql_xz', 'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz',  &
                           'qv_xz', 's_xz',                                     &
-                          'theta_xz', 'thetal_xz', 'thetav_xz', 'ti_xz' )
+                          'theta_xz', 'thetal_xz', 'thetav_xz' )
 
                       grid_x = 'x'
@@ -3837,4 +3854,8 @@
                                                        grid_z )
                       ENDIF
+                      
+                      IF ( .NOT. found )                                       &
+                         CALL doq_define_netcdf_grid( do2d(av,i), found,       &
+                                                      grid_x, grid_y, grid_z )
 
 !
@@ -4695,4 +4716,9 @@
                                                        grid_z )
                       ENDIF
+                      
+                      IF ( .NOT. found )                                       &
+                         CALL doq_define_netcdf_grid(                          &
+                                                    do2d(av,i), found, grid_x, &
+                                                    grid_y, grid_z           )
 !
 !--                   Check for user-defined quantities
@@ -6000,5 +6026,5 @@
                        'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
                        'rho_sea_water', 's', 'sa', &
-                       'theta', 'thetal', 'thetav', 'ti', 'u', 'v' )
+                       'theta', 'thetal', 'thetav', 'u', 'v' )
 
                    grid_z = 'zu'
Index: palm/trunk/SOURCE/palm.f90
===================================================================
--- palm/trunk/SOURCE/palm.f90	(revision 4038)
+++ palm/trunk/SOURCE/palm.f90	(revision 4039)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Rename subroutines in module for diagnostic quantities
+! 
+! 4017 2019-06-06 12:16:46Z schwenkel
 ! new module for calculation and output of diagnostic quantities added
 ! 
@@ -305,5 +308,5 @@
 
     USE diagnostic_output_quantities_mod,                                      &
-        ONLY:  diagnostic_output_quantities_calculate
+        ONLY:  doq_calculate
 
     USE indices,                                                               &
@@ -559,5 +562,5 @@
 !-- If required, output of initial arrays
     IF ( do2d_at_begin )  THEN
-       CALL diagnostic_output_quantities_calculate    !TODO, will be called twice
+       CALL doq_calculate    !TODO, will be called twice
 
        CALL data_output_2d( 'xy', 0 )
@@ -567,5 +570,5 @@
 
     IF ( do3d_at_begin )  THEN
-       CALL diagnostic_output_quantities_calculate    !TODO, will be called twice
+       CALL doq_calculate    !TODO, will be called twice
 
        CALL data_output_3d( 0 )
Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4038)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4039)
@@ -28,4 +28,7 @@
 ! -----------------
 ! $Id$
+! Bugfix for masked data output
+! 
+! 4008 2019-05-30 09:50:11Z moh.hefny
 ! Bugfix in check variable when a variable's string is less than 3
 ! characters is processed. All variables now are checked if they
@@ -11197,5 +11200,5 @@
 !
 !-- Resort the array to be output, if not done above
-    IF ( .NOT. resorted )  THEN
+    IF ( found  .AND.  .NOT. resorted )  THEN
        IF ( .NOT. mask_surface(mid) )  THEN
 !
Index: palm/trunk/SOURCE/read_restart_data_mod.f90
===================================================================
--- palm/trunk/SOURCE/read_restart_data_mod.f90	(revision 4038)
+++ palm/trunk/SOURCE/read_restart_data_mod.f90	(revision 4039)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! input of uu_av, vv_av, ww_av added
+! 
+! 4017 2019-06-06 12:16:46Z schwenkel
 ! bugfix for r3998, allocation of 3d temporary arrays of various dimensions revised
 ! 
@@ -111,4 +114,6 @@
 !>
 !> @todo: Revise max_pr_cs (profiles for chemistry)
+!> @todo: Modularize reading of restart data for diagnostic quantities, which
+!>        is not possible with the current module-interface structure
 !------------------------------------------------------------------------------!
  MODULE read_restart_data_mod
@@ -135,5 +140,5 @@
 
     USE diagnostic_output_quantities_mod,                                      &
-        ONLY:  ti_av
+        ONLY:  ti_av, uu_av, vv_av, ww_av
 
     USE grid_variables,                                                        &
@@ -1553,4 +1558,15 @@
                    u_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                               &
                       tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                      
+                CASE ( 'uu_av' )
+                   IF ( .NOT. ALLOCATED( uu_av ) )  THEN
+                      ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                   ENDIF
+                   IF ( k == 1 )  THEN
+                      ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,             &
+                                                    nxl_on_file:nxr_on_file) )
+                      READ ( 13 )  tmp_3d_non_standard
+                   ENDIF
+                   uu_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
 
                 CASE ( 'u_m_l' )
@@ -1614,4 +1630,15 @@
                    v_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                               &
                       tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                      
+                CASE ( 'vv_av' )
+                   IF ( .NOT. ALLOCATED( vv_av ) )  THEN
+                      ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                   ENDIF
+                   IF ( k == 1 )  THEN
+                      ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,             &
+                                                    nxl_on_file:nxr_on_file) )
+                      READ ( 13 )  tmp_3d_non_standard
+                   ENDIF
+                   vv_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
 
                 CASE ( 'v_m_l' )
@@ -1680,4 +1707,15 @@
                    w_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                               &
                       tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                      
+                CASE ( 'ww_av' )
+                   IF ( .NOT. ALLOCATED( ww_av ) )  THEN
+                      ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                   ENDIF
+                   IF ( k == 1 )  THEN
+                      ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,             &
+                                                    nxl_on_file:nxr_on_file) )
+                      READ ( 13 )  tmp_3d_non_standard
+                   ENDIF
+                   ww_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
 
                 CASE ( 'w_m_l' )
Index: palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 4038)
+++ palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 4039)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Modularize diagnostic output
+! 
+! 3994 2019-05-22 18:08:09Z suehring
 ! output of turbulence intensity added
 ! 
@@ -286,7 +289,4 @@
     USE cpulog,                                                                &
         ONLY:  cpu_log, log_point
-
-    USE diagnostic_output_quantities_mod,                                      &
-        ONLY:  ti, ti_av
 
     USE indices,                                                               &
@@ -478,10 +478,4 @@
                 ENDIF
                 tsurf_av = 0.0_wp
-
-             CASE ( 'ti' )
-                IF ( .NOT. ALLOCATED( ti_av ) )  THEN
-                   ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
-                ENDIF
-                ti_av = 0.0_wp
 
              CASE ( 'u' )
@@ -1031,15 +1025,4 @@
              ENDIF
 
-          CASE ( 'ti' )
-             IF ( ALLOCATED( ti_av ) ) THEN
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      DO  k = nzb, nzt+1
-                         ti_av(k,j,i) = ti_av(k,j,i) + ti(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
           CASE ( 'u' )
              IF ( ALLOCATED( u_av ) ) THEN
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 4038)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 4039)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Rename subroutines in module for diagnostic quantities
+! 
+! 4029 2019-06-14 14:04:35Z raasch
 ! exchange of ghost points and boundary conditions separated for chemical species and SALSA module,
 ! bugfix: decycling of chemistry species after nesting data transfer
@@ -578,5 +581,5 @@
 
     USE diagnostic_output_quantities_mod,                                                          &
-        ONLY:  diagnostic_output_quantities_calculate,                                             &
+        ONLY:  doq_calculate,                                                                      &
                timestep_number_at_prev_calc
 
@@ -1730,5 +1733,5 @@
           IF ( time_do_sla >= dt_averaging_input )  THEN
              IF ( current_timestep_number > timestep_number_at_prev_calc )                         &
-                CALL diagnostic_output_quantities_calculate
+                CALL doq_calculate
 
              CALL sum_up_3d_data
@@ -1811,5 +1814,5 @@
        IF ( time_do2d_xy >= dt_do2d_xy  .AND.  time_since_reference_point >= skip_time_do2d_xy )  THEN
           IF ( current_timestep_number > timestep_number_at_prev_calc )                            &
-             CALL diagnostic_output_quantities_calculate
+             CALL doq_calculate
 
           CALL data_output_2d( 'xy', 0 )
@@ -1819,5 +1822,5 @@
           IF ( current_timestep_number > timestep_number_at_prev_calc )                            &
 
-             CALL diagnostic_output_quantities_calculate
+             CALL doq_calculate
           CALL data_output_2d( 'xz', 0 )
           time_do2d_xz = MOD( time_do2d_xz, MAX( dt_do2d_xz, dt_3d ) )
@@ -1825,5 +1828,5 @@
        IF ( time_do2d_yz >= dt_do2d_yz  .AND.  time_since_reference_point >= skip_time_do2d_yz )  THEN
           IF ( current_timestep_number > timestep_number_at_prev_calc )                            &
-             CALL diagnostic_output_quantities_calculate
+             CALL doq_calculate
 
           CALL data_output_2d( 'yz', 0 )
@@ -1835,5 +1838,5 @@
        IF ( time_do3d >= dt_do3d  .AND.  time_since_reference_point >= skip_time_do3d )  THEN
           IF ( current_timestep_number > timestep_number_at_prev_calc )                            &
-             CALL diagnostic_output_quantities_calculate
+             CALL doq_calculate
 
           CALL data_output_3d( 0 )
@@ -1847,5 +1850,5 @@
                .AND.  time_since_reference_point >= skip_time_domask(mid) )  THEN
              IF ( current_timestep_number > timestep_number_at_prev_calc )                         &
-                CALL diagnostic_output_quantities_calculate
+                CALL doq_calculate
 
              CALL data_output_mask( 0 )
Index: palm/trunk/SOURCE/write_restart_data_mod.f90
===================================================================
--- palm/trunk/SOURCE/write_restart_data_mod.f90	(revision 4038)
+++ palm/trunk/SOURCE/write_restart_data_mod.f90	(revision 4039)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Modularize diagnostic output
+! 
+! 4017 2019-06-06 12:16:46Z schwenkel
 ! output of turbulence intensity added
 ! 
@@ -95,7 +98,4 @@
         ONLY:  day_of_year_init, time_utc_init
 
-    USE diagnostic_output_quantities_mod,                                      &
-        ONLY:  ti_av
-
     USE grid_variables,                                                        &
         ONLY:  dx, dy
@@ -1067,9 +1067,4 @@
        ENDIF
 
-       IF ( ALLOCATED( ti_av ) )  THEN
-          CALL wrd_write_string( 'ti_av' ) 
-          WRITE ( 14 )  ti_av
-       ENDIF
-
        IF ( ALLOCATED( ts_av ) )  THEN
           CALL wrd_write_string( 'ts_av' ) 
