Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 4330)
+++ /palm/trunk/SOURCE/Makefile	(revision 4331)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# Move diagnostic surface output to diagnostic_output_quantities
+# 
+# 4309 2019-11-26 18:49:59Z suehring
 # Add dependency to parallel random generator for synthetic turbulence generator
 # 
@@ -519,5 +522,6 @@
 diagnostic_output_quantities_mod.o: \
 	mod_kinds.o \
-	modules.o
+	modules.o \
+	surface_layer_fluxes_mod.o
 diffusion_s.o: \
 	mod_kinds.o \
@@ -1128,4 +1132,5 @@
 	basic_constants_and_equations_mod.o \
 	cpulog_mod.o \
+	diagnostic_output_quantities_mod.o \
 	land_surface_model_mod.o \
 	mod_kinds.o \
Index: /palm/trunk/SOURCE/average_3d_data.f90
===================================================================
--- /palm/trunk/SOURCE/average_3d_data.f90	(revision 4330)
+++ /palm/trunk/SOURCE/average_3d_data.f90	(revision 4331)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Move 2-m potential temperature output to diagnostic_output_quantities
+! 
+! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
 ! 
@@ -315,14 +318,4 @@
              ENDIF
 
-         CASE ( 'theta_2m*' )
-             IF ( ALLOCATED( pt_2m_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      pt_2m_av(j,i) = pt_2m_av(j,i) / REAL( average_count_3d, KIND=wp )
-                   ENDDO
-                ENDDO
-                CALL exchange_horiz_2d( pt_2m_av, nbgp )
-             ENDIF
-
           CASE ( 't*' )
              IF ( ALLOCATED( ts_av ) ) THEN
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 4330)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 4331)
@@ -25,4 +25,7 @@
 ! -----------------
 ! 4172 2019-08-20 11:55:33Z oliver.maas
+! Move 2-m potential temperature output to diagnostic_output_quantities
+! 
+! 11:55:33Z oliver.maas
 ! removed message PA0421, concerning old parameter recycling_yshift
 ! 
@@ -2385,5 +2388,5 @@
 
           CASE ( 'ghf*', 'lwp*', 'ol*', 'qsws*', 'r_a*',                       &
-                 'shf*', 'ssws*', 't*', 'theta_2m*', 'tsurf*', 'us*',          &
+                 'shf*', 'ssws*', 't*', 'tsurf*', 'us*',                       &
                  'z0*', 'z0h*', 'z0q*' )
              IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
@@ -2424,10 +2427,4 @@
                                  'res passive_scalar = .TRUE.'
                 CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
-             ENDIF
-!
-!--          Activate calculation of 2m temperature if output is requested
-             IF ( TRIM( var ) == 'theta_2m*' )  THEN
-                do_output_at_2m = .TRUE.
-                unit = 'K'
              ENDIF
 
Index: /palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_2d.f90	(revision 4330)
+++ /palm/trunk/SOURCE/data_output_2d.f90	(revision 4331)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Move 2-m potential temperature output to diagnostic_output_quantities
+! 
+! 4329 2019-12-10 15:46:36Z motisi
 ! Renamed wall_flags_0 to wall_flags_static_0
 ! 
@@ -1001,36 +1004,4 @@
                 ENDIF
                 IF ( mode == 'xy' )  level_z = zu
-
-             CASE ( 'theta_2m*_xy' )        ! 2d-array
-                IF ( av == 0 )  THEN
-                   DO  m = 1, surf_def_h(0)%ns
-                      i = surf_def_h(0)%i(m)
-                      j = surf_def_h(0)%j(m)
-                      local_pf(i,j,nzb+1) = surf_def_h(0)%pt_2m(m)
-                   ENDDO
-                   DO  m = 1, surf_lsm_h%ns
-                      i = surf_lsm_h%i(m)
-                      j = surf_lsm_h%j(m)
-                      local_pf(i,j,nzb+1) = surf_lsm_h%pt_2m(m)
-                   ENDDO
-                   DO  m = 1, surf_usm_h%ns
-                      i = surf_usm_h%i(m)
-                      j = surf_usm_h%j(m)
-                      local_pf(i,j,nzb+1) = surf_usm_h%pt_2m(m)
-                   ENDDO
-                ELSE
-                   IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN
-                      ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
-                      pt_2m_av = REAL( fill_value, KIND = wp )
-                   ENDIF
-                   DO  i = nxl, nxr
-                      DO  j = nys, nyn
-                         local_pf(i,j,nzb+1) = pt_2m_av(j,i)
-                      ENDDO
-                   ENDDO
-                ENDIF
-                resorted = .TRUE.
-                two_d = .TRUE.
-                level_z(nzb+1) = zu(nzb+1)
 
              CASE ( 'w_xy', 'w_xz', 'w_yz' )
Index: /palm/trunk/SOURCE/data_output_mask.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_mask.f90	(revision 4330)
+++ /palm/trunk/SOURCE/data_output_mask.f90	(revision 4331)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Formatting adjustment
+! 
+! 4329 2019-12-10 15:46:36Z motisi
 ! Renamed wall_flags_0 to wall_flags_static_0
 ! 
@@ -603,5 +606,5 @@
              IF ( .NOT. found )  THEN
                 CALL doq_output_mask( av, domask(mid,av,ivar), found, local_pf,   &
-                                   mid)
+                                      mid)
              ENDIF
 !
Index: /palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
===================================================================
--- /palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4330)
+++ /palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4331)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! - 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
 ! 
@@ -75,16 +79,20 @@
  
     USE arrays_3d,                                                             &
-        ONLY:  ddzu, u, v, w
-! 
-!     USE averaging
-! 
-!     USE basic_constants_and_equations_mod,                                     &
-!         ONLY:  c_p, lv_d_cp, l_v
-! 
-!     USE bulk_cloud_model_mod,                                                  &
-!         ONLY:  bulk_cloud_model
-! 
+        ONLY:  ddzu,                                                           &
+               pt,                                                             &
+               u,                                                              &
+               v,                                                              &
+               w,                                                              &
+               zu,                                                             &
+               zw
+
+    USE basic_constants_and_equations_mod,                                     &
+        ONLY:  kappa
+
     USE control_parameters,                                                    &
-        ONLY:  current_timestep_number, varnamelength
+        ONLY:  current_timestep_number,                                        &
+               data_output,                                                    &
+               message_string,                                                 &
+               varnamelength
 ! 
 !     USE cpulog,                                                                &
@@ -93,37 +101,26 @@
    USE grid_variables,                                                         &
         ONLY:  ddx, ddy
-! 
+
     USE indices,                                                               &
-        ONLY:  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_static_0
-! 
+        ONLY:  nbgp,                                                           &
+               nxl,                                                            &
+               nxlg,                                                           & 
+               nxr,                                                            & 
+               nxrg,                                                           &
+               nyn,                                                            &
+               nyng,                                                           &
+               nys,                                                            &
+               nysg,                                                           &
+               nzb,                                                            &
+               nzt,                                                            &
+               wall_flags_static_0
+
     USE kinds
-!
-!     USE land_surface_model_mod,                                                &
-!         ONLY:  zs
-! 
-!     USE module_interface,                                                      &
-!         ONLY:  module_interface_data_output_2d
-! 
-! #if defined( __netcdf )
-!     USE NETCDF
-! #endif
-! 
-!     USE netcdf_interface,                                                      &
-!         ONLY:  fill_value, id_set_xy, id_set_xz, id_set_yz, id_var_do2d,       &
-!                id_var_time_xy, id_var_time_xz, id_var_time_yz, nc_stat,        &
-!                netcdf_data_format, netcdf_handle_error
-! 
-!     USE particle_attributes,                                                   &
-!         ONLY:  grid_particles, number_of_particles, particle_advection_start,  &
-!                particles, prt_count
-!     
-!     USE pegrid
-! 
-!     USE surface_mod,                                                           &
-!         ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_def_h,           &
-!                surf_lsm_h, surf_usm_h
-! 
-!     USE turbulence_closure_mod,                                                &
-!         ONLY:  tcm_data_output_2d
+
+    USE surface_mod,                                                           &
+        ONLY:  surf_def_h,                                                     &
+               surf_lsm_h,                                                     &
+               surf_type,                                                      &
+               surf_usm_h
 
 
@@ -137,13 +134,15 @@
     LOGICAL ::  prepared_diagnostic_output_quantities = .FALSE.    !< flag indicating whether output is p
 
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_2m     !< 2-m air potential temperature
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_2m_av  !< averaged 2-m air potential temperature
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  uv_10m    !< horizontal wind speed at 10m
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  uv_10m_av !< averaged horizontal wind speed at 10m
+
     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 ::  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 ::  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 ::  vv_av  !< mean of vv    
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ww     !< ww
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ww_av  !< mean of ww
@@ -160,6 +159,8 @@
            prepared_diagnostic_output_quantities,                              &
            timestep_number_at_prev_calc,                                       &
+           pt_2m_av,                                                           &
            ti_av,                                                              &
            uu_av,                                                              &
+           uv_10m_av,                                                          &
            vv_av,                                                              &
            ww_av
@@ -270,4 +271,16 @@
              ENDIF
              ww_av = 0.0_wp
+             
+          CASE ( 'theta_2m*' )
+             IF ( .NOT. ALLOCATED( pt_2m_av ) )  THEN
+                ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
+             ENDIF
+             pt_2m_av = 0.0_wp
+             
+          CASE ( 'wspeed_10m*' )
+             IF ( .NOT. ALLOCATED( uv_10m_av ) )  THEN
+                ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
+             ENDIF
+             uv_10m_av = 0.0_wp
 
           CASE DEFAULT
@@ -323,4 +336,22 @@
                 ENDDO
              ENDIF
+             
+          CASE ( 'theta_2m*' )
+             IF ( ALLOCATED( pt_2m_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      pt_2m_av(j,i) = pt_2m_av(j,i) + pt_2m(j,i)
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          CASE ( 'wspeed_10m*' )
+             IF ( ALLOCATED( uv_10m_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      uv_10m_av(j,i) = uv_10m_av(j,i) + uv_10m(j,i)
+                   ENDDO
+                ENDDO
+             ENDIF
 
           CASE DEFAULT
@@ -377,4 +408,24 @@
              ENDIF
 
+         CASE ( 'theta_2m*' )
+            IF ( ALLOCATED( pt_2m_av ) ) THEN
+               DO  i = nxlg, nxrg
+                  DO  j = nysg, nyng
+                     pt_2m_av(j,i) = pt_2m_av(j,i) / REAL( average_count_3d, KIND=wp )
+                  ENDDO
+               ENDDO
+               CALL exchange_horiz_2d( pt_2m_av, nbgp )
+            ENDIF
+
+         CASE ( 'wspeed_10m*' )
+            IF ( ALLOCATED( uv_10m_av ) ) THEN
+               DO  i = nxlg, nxrg
+                  DO  j = nysg, nyng
+                     uv_10m_av(j,i) = uv_10m_av(j,i) / REAL( average_count_3d, KIND=wp )
+                  ENDDO
+               ENDDO
+               CALL exchange_horiz_2d( uv_10m_av, nbgp )
+            ENDIF
+
        END SELECT
 
@@ -389,5 +440,5 @@
 !> Check data output for diagnostic output
 !------------------------------------------------------------------------------!
- SUBROUTINE doq_check_data_output( var, unit )
+ SUBROUTINE doq_check_data_output( var, unit, i, ilen, k )
 
     IMPLICIT NONE
@@ -395,4 +446,8 @@
     CHARACTER (LEN=*) ::  unit  !< 
     CHARACTER (LEN=*) ::  var   !<
+
+    INTEGER(iwp), OPTIONAL, INTENT(IN) ::  i     !< Current element of data_output
+    INTEGER(iwp), OPTIONAL, INTENT(IN) ::  ilen  !< Length of current entry in data_output
+    INTEGER(iwp), OPTIONAL, INTENT(IN) ::  k     !< Output is xy mode? 0 = no, 1 = yes
 
     SELECT CASE ( TRIM( var ) )
@@ -409,5 +464,19 @@
        CASE ( 'ww' )
           unit = 'm2/s2'
-             
+!
+!--    Treat horizotal cross-section output quanatities
+       CASE ( 'theta_2m*', 'wspeed_10m*' )
+!
+!--       Check if output quantity is _xy only.
+          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
+             message_string = 'illegal value for data_output: "' //            &
+                              TRIM( var ) // '" & only 2d-horizontal ' //      &
+                              'cross sections are allowed for this value'
+             CALL message( 'diagnostic_output', 'PA0111', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( TRIM( var ) == 'theta_2m*'   )  unit = 'K'
+          IF ( TRIM( var ) == 'wspeed_10m*' )  unit = 'm/s'
+
        CASE DEFAULT
           unit = 'illegal'
@@ -443,5 +512,12 @@
           grid_x = 'x'
           grid_y = 'y'
-          grid_z = 'zu'    
+          grid_z = 'zu'
+!
+!--    s grid surface variables
+       CASE ( 'theta_2m*_xy', 'wspeed_10m*' )
+
+          grid_x = 'x'
+          grid_y = 'y'
+          grid_z = 'zu'
 !
 !--    u grid
@@ -572,4 +648,48 @@
           
           IF ( mode == 'xy' )  grid = 'zw'
+          
+       CASE ( 'theta_2m*_xy' )        ! 2d-array
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   local_pf(i,j,nzb+1) = pt_2m(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN
+                ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
+                pt_2m_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   local_pf(i,j,nzb+1) = pt_2m_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+          resorted = .TRUE.
+          two_d    = .TRUE.
+          grid     = 'zu1'
+          
+       CASE ( 'wspeed_10m*_xy' )        ! 2d-array
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   local_pf(i,j,nzb+1) = uv_10m(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN
+                ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
+                uv_10m_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   local_pf(i,j,nzb+1) = uv_10m_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+          resorted = .TRUE.
+          two_d    = .TRUE.
+          grid     = 'zu1'
 
        CASE DEFAULT
@@ -776,5 +896,4 @@
           flag_nr = 3
 
-
        CASE DEFAULT
           found = .FALSE.
@@ -864,5 +983,4 @@
                 uu = 0.0_wp
              ENDIF
-
 !
 !--       Allocate array for vv
@@ -872,5 +990,4 @@
                 vv = 0.0_wp
              ENDIF
-
 !
 !--       Allocate array for ww
@@ -879,4 +996,18 @@
                 ALLOCATE( ww(nzb:nzt+1,nys:nyn,nxl:nxr) )
                 ww = 0.0_wp
+             ENDIF
+!
+!--       Allocate array for 2-m potential temperature
+          CASE ( 'theta_2m*' )
+             IF ( .NOT. ALLOCATED( pt_2m ) )  THEN
+                ALLOCATE( pt_2m(nys:nyn,nxl:nxr) )
+                pt_2m = 0.0_wp
+             ENDIF
+!
+!--       Allocate array for 10-m wind speed
+          CASE ( 'wspeed_10m*' )
+             IF ( .NOT. ALLOCATED( uv_10m ) )  THEN
+                ALLOCATE( uv_10m(nys:nyn,nxl:nxr) )
+                uv_10m = 0.0_wp
              ENDIF
 
@@ -900,6 +1031,10 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i, j, k    !< grid loop index in x-, y-, z-direction
+    INTEGER(iwp) ::  i          !< grid index x-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
 
 
@@ -973,4 +1108,34 @@
                 ENDDO
              ENDDO
+!
+!--       2-m potential temperature
+          CASE ( 'theta_2m*' )
+!
+!--          2-m potential temperature is caluclated from surface arrays. In
+!--          case the 2m level is below the first grid point, MOST is applied,
+!--          else, linear interpolation between two vertical grid levels is
+!--          applied. To access all surfaces, iterate over all horizontally-
+!--          upward facing surface types.
+             surf => surf_def_h(0)
+             CALL calc_pt_2m
+             surf => surf_lsm_h
+             CALL calc_pt_2m
+             surf => surf_usm_h
+             CALL calc_pt_2m
+!
+!--       10-m wind speed
+          CASE ( 'wspeed_10m*' )
+!
+!--          10-m wind speed is caluclated from surface arrays. In
+!--          case the 10m level is below the first grid point, MOST is applied,
+!--          else, linear interpolation between two vertical grid levels is
+!--          applied. To access all surfaces, iterate over all horizontally-
+!--          upward facing surface types.
+             surf => surf_def_h(0)
+             CALL calc_wind_10m
+             surf => surf_lsm_h
+             CALL calc_wind_10m
+             surf => surf_usm_h
+             CALL calc_wind_10m
 
        END SELECT
@@ -981,11 +1146,131 @@
 !     CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' )
 
- END SUBROUTINE doq_calculate
-
-
+!
+!-- The following block contains subroutines to calculate diagnostic 
+!-- surface quantities. 
+    CONTAINS
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> ...
+!> Calculation of 2-m potential temperature.
+!------------------------------------------------------------------------------!
+       SUBROUTINE calc_pt_2m
+
+          USE surface_layer_fluxes_mod,                                        &
+              ONLY:  psi_h
+
+          IMPLICIT NONE
+
+          INTEGER(iwp) ::  kk     !< running index along the z-dimension
+          INTEGER(iwp) ::  m      !< running index for surface elements
+ 
+          DO  m = 1, surf%ns
+
+             i = surf%i(m)
+             j = surf%j(m)
+             k = surf%k(m)
+!
+!--          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
+                pt_2m(j,i) = surf%pt_surface(m) + surf%ts(m) / kappa           &
+                                * ( LOG( 2.0_wp /  surf%z0h(m) )               &
+                                  - psi_h( 2.0_wp / surf%ol(m) )               &
+                                  + psi_h( surf%z0h(m) / surf%ol(m) ) )
+!
+!--          If 2-m level is above the first grid level, 2-m temperature
+!--          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. 
+             ELSE
+!
+!--             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. 
+                pt_2m(j,i) = pt(kk-1,j,i) +                                    &
+                              ( zw(k-1) + 2.0_wp - zu(kk-1)     ) *            &
+                              ( pt(kk,j,i)       - pt(kk-1,j,i) ) /            &
+                              ( zu(kk)           - zu(kk-1)     )
+             ENDIF
+
+          ENDDO
+
+       END SUBROUTINE calc_pt_2m
+       
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculation of 10-m wind speed.
+!------------------------------------------------------------------------------!
+       SUBROUTINE calc_wind_10m
+
+          USE surface_layer_fluxes_mod,                                        &
+              ONLY:  psi_m
+
+          IMPLICIT NONE
+
+          INTEGER(iwp) ::  kk     !< running index along the z-dimension
+          INTEGER(iwp) ::  m      !< running index for surface elements
+
+          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
+
+             i = surf%i(m)
+             j = surf%j(m)
+             k = surf%k(m)
+!
+!--          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
+                uv_10m(j,i) = surf%us(m) / kappa                               &
+                          * ( LOG( 10.0_wp /  surf%z0(m) )                     &
+                              - psi_m( 10.0_wp    / surf%ol(m) )               &
+                              + psi_m( surf%z0(m) / surf%ol(m) ) )
+!
+!--          If 10-m level is above the first grid level, 10-m wind speed
+!--          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. 
+             ELSE
+!
+!--             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 
+                ENDDO
+!
+!--             kk defines the index of the first grid level >= 10m.
+                uv_l = SQRT( ( 0.5_wp * ( u(kk-1,j,i) + u(kk-1,j,i+1) ) )**2   &
+                           + ( 0.5_wp * ( v(kk-1,j,i) + v(kk-1,j+1,i) ) )**2 )
+
+                uv_u = SQRT( ( 0.5_wp * ( u(kk,j,i)   + u(kk,j,i+1)   ) )**2   &
+                           + ( 0.5_wp * ( v(kk,j,i)   + v(kk,j+1,i)   ) )**2 )
+
+                uv_10m(j,i) = uv_l + ( zw(k-1) + 10.0_wp - zu(kk-1) ) *        &
+                                     ( uv_u              - uv_l     ) /        &
+                                     ( zu(kk)            - zu(kk-1) )
+
+             ENDIF
+
+          ENDDO
+
+       END SUBROUTINE calc_wind_10m
+
+ END SUBROUTINE doq_calculate
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Preparation of the diagnostic output, counting of the module-specific
+!> output quantities and gathering of the output names.
 !------------------------------------------------------------------------------!
  SUBROUTINE doq_prepare
@@ -1173,21 +1458,31 @@
     IMPLICIT NONE
 
+    IF ( ALLOCATED( pt_2m_av ) )  THEN
+       CALL wrd_write_string( 'pt_2m_av' )
+       WRITE ( 14 )  pt_2m_av
+    ENDIF
+
     IF ( ALLOCATED( ti_av ) )  THEN
-       CALL wrd_write_string( 'ti_av' ) 
+       CALL wrd_write_string( 'ti_av' )
        WRITE ( 14 )  ti_av
     ENDIF
-    
+
     IF ( ALLOCATED( uu_av ) )  THEN
-       CALL wrd_write_string( 'uu_av' ) 
+       CALL wrd_write_string( 'uu_av' )
        WRITE ( 14 )  uu_av
     ENDIF
-       
+
+    IF ( ALLOCATED( uv_10m_av ) )  THEN
+       CALL wrd_write_string( 'uv_10m_av' )
+       WRITE ( 14 )  uv_10m_av
+    ENDIF
+
     IF ( ALLOCATED( vv_av ) )  THEN
-       CALL wrd_write_string( 'vv_av' ) 
+       CALL wrd_write_string( 'vv_av' )
        WRITE ( 14 )  vv_av
     ENDIF
-       
+
     IF ( ALLOCATED( ww_av ) )  THEN
-       CALL wrd_write_string( 'ww_av' ) 
+       CALL wrd_write_string( 'ww_av' )
        WRITE ( 14 )  ww_av
     ENDIF
Index: /palm/trunk/SOURCE/module_interface.f90
===================================================================
--- /palm/trunk/SOURCE/module_interface.f90	(revision 4330)
+++ /palm/trunk/SOURCE/module_interface.f90	(revision 4331)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Change interface for doq_check_data_output, in order to perform further 
+! output checks.
+! 
+! 4281 2019-10-29 15:15:39Z schwenkel
 ! Added dynamics boundary conditions
 ! 
@@ -872,5 +876,5 @@
 
     IF ( unit == 'illegal' )  THEN
-       CALL doq_check_data_output( variable, unit )
+       CALL doq_check_data_output( variable, unit, i, ilen, k )
     ENDIF
 
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 4330)
+++ /palm/trunk/SOURCE/modules.f90	(revision 4331)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! - do_output_at_2m, pt_2m_av
+! 
+! 4329 2019-12-10 15:46:36Z motisi
 ! Renamed wall_flags_0 to wall_flags_static_0
 ! 
@@ -413,5 +416,4 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lwp_av                 !< avg. liquid water path
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ol_av                  !< avg. Obukhov length
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_2m_av               !< avg. 2m- air potential temperature 
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  qsws_av                !< avg. surface moisture flux
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  r_a_av                 !< avg. resistance
@@ -707,5 +709,4 @@
     LOGICAL ::  do2d_at_begin = .FALSE.                          !< namelist parameter
     LOGICAL ::  do3d_at_begin = .FALSE.                          !< namelist parameter
-    LOGICAL ::  do_output_at_2m = .FALSE.                        !< flag for activating calculation of potential temperature at z = 2 m
     LOGICAL ::  do_sum = .FALSE.                                 !< contribute to time average of profile data?
     LOGICAL ::  dp_external = .FALSE.                            !< namelist parameter
Index: /palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4330)
+++ /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4331)
@@ -27,4 +27,7 @@
 ! -----------------
 ! $Id$
+! Typo corrected
+! 
+! 4329 2019-12-10 15:46:36Z motisi
 ! Renamed wall_flags_0 to wall_flags_static_0
 ! 
@@ -1135,5 +1138,5 @@
                                         'defined on top of an artificially '// &
                                         'created building grid point ' //      &
-                                        '(emerged from the filering) - ' //    &
+                                        '(emerged from the filtering) - ' //   &
                                         'LAD profile is omitted at this ' //   &
                                         'grid point: (i,j) = ', i, j
Index: /palm/trunk/SOURCE/read_restart_data_mod.f90
===================================================================
--- /palm/trunk/SOURCE/read_restart_data_mod.f90	(revision 4330)
+++ /palm/trunk/SOURCE/read_restart_data_mod.f90	(revision 4331)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Enable restart data for 2-m potential temperature output
+! 
+! 4301 2019-11-22 12:09:09Z oliver.maas
 ! removed recycling_yshift
 ! 
@@ -101,5 +104,10 @@
 
     USE diagnostic_output_quantities_mod,                                      &
-        ONLY:  ti_av, uu_av, vv_av, ww_av
+        ONLY:  pt_2m_av,                                                       &
+               ti_av,                                                          &
+               uu_av,                                                          &
+               uv_10m_av,                                                      &
+               vv_av,                                                          &
+               ww_av
 
     USE grid_variables,                                                        &
@@ -1525,4 +1533,12 @@
                    uu_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
 
+                CASE ( 'uv_10m_av' )
+                   IF ( .NOT. ALLOCATED( uv_10m_av ) )  THEN
+                      ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   uv_10m_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =                           &
+                      tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
                 CASE ( 'u_m_l' )
                    IF ( k == 1 )  THEN
Index: /palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- /palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 4330)
+++ /palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 4331)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Move 2-m potential temperature output to diagnostic_output_quantities
+! 
+! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
 ! 
@@ -64,11 +67,11 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  dzw, d_exner, e, heatflux_output_conversion, p,    &
-               pt, q, ql, ql_c, ql_v, s, u, v, vpt, w,                 &
+        ONLY:  dzw, d_exner, e, heatflux_output_conversion, p,                 &
+               pt, q, ql, ql_c, ql_v, s, u, v, vpt, w,                         &
                waterflux_output_conversion
 
     USE averaging,                                                             &
         ONLY:  e_av, ghf_av, lpt_av, lwp_av, ol_av, p_av, pc_av, pr_av, pt_av, &
-               pt_2m_av, q_av, ql_av, ql_c_av, ql_v_av, ql_vp_av, qsws_av,     &
+               q_av, ql_av, ql_c_av, ql_v_av, ql_vp_av, qsws_av,               &
                qv_av, r_a_av, s_av, shf_av, ssws_av, ts_av, tsurf_av, u_av,    &
                us_av, v_av, vpt_av, w_av, z0_av, z0h_av, z0q_av
@@ -296,10 +299,4 @@
                 ENDIF
                 vpt_av = 0.0_wp
-
-             CASE ( 'theta_2m*' )
-                IF ( .NOT. ALLOCATED( pt_2m_av ) )  THEN
-                   ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                pt_2m_av = 0.0_wp
 
              CASE ( 'w' )
@@ -730,33 +727,4 @@
              ENDIF
 
-          CASE ( 'theta_2m*' )
-             IF ( ALLOCATED( pt_2m_av ) ) THEN    
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      match_def = surf_def_h(0)%start_index(j,i) <=            &
-                                  surf_def_h(0)%end_index(j,i)
-                      match_lsm = surf_lsm_h%start_index(j,i) <=               &
-                                  surf_lsm_h%end_index(j,i)
-                      match_usm = surf_usm_h%start_index(j,i) <=               &
-                                  surf_usm_h%end_index(j,i)
-
-                      IF ( match_def )  THEN
-                         m = surf_def_h(0)%end_index(j,i)
-                         pt_2m_av(j,i) = pt_2m_av(j,i) +                       &
-                                         surf_def_h(0)%pt_2m(m)
-                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
-                         m = surf_lsm_h%end_index(j,i)
-                         pt_2m_av(j,i) = pt_2m_av(j,i) +                       &
-                                         surf_lsm_h%pt_2m(m)
-                      ELSEIF ( match_usm )  THEN
-                         m = surf_usm_h%end_index(j,i)
-                         pt_2m_av(j,i) = pt_2m_av(j,i) +                       &
-                                         surf_usm_h%pt_2m(m)
-                      ENDIF
-                   ENDDO
-                ENDDO
-             ENDIF
-             
-             
           CASE ( 't*' )
              IF ( ALLOCATED( ts_av ) ) THEN 
Index: /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 4330)
+++ /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 4331)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! Calculation of diagnostic-only 2-m potential temperature moved to 
+! diagnostic_output_quantities.
+! 
+! 4298 2019-11-21 15:59:16Z suehring
 ! Calculation of 2-m temperature adjusted to the case the 2-m level is above 
 ! the first grid point.
@@ -108,5 +112,5 @@
                constant_waterflux, coupling_mode,                              &
                debug_output_timestep,                                          &
-               do_output_at_2m, humidity,                                      &
+               humidity,                                                       &
                ibc_e_b, ibc_pt_b, indoor_model,                                &
                land_surface, large_scale_forcing, lsf_surf, message_string,    &
@@ -162,5 +166,9 @@
     PRIVATE
 
-    PUBLIC init_surface_layer_fluxes, phi_m, surface_layer_fluxes
+    PUBLIC init_surface_layer_fluxes,                                          &
+           phi_m,                                                              &
+           psi_h,                                                              &
+           psi_m,                                                              &
+           surface_layer_fluxes
 
     INTERFACE init_surface_layer_fluxes
@@ -171,4 +179,12 @@
        MODULE PROCEDURE phi_m
     END INTERFACE phi_m
+
+    INTERFACE psi_h
+       MODULE PROCEDURE psi_h
+    END INTERFACE psi_h
+
+    INTERFACE psi_m
+       MODULE PROCEDURE psi_m
+    END INTERFACE psi_m
 
     INTERFACE surface_layer_fluxes
@@ -259,7 +275,4 @@
           CALL calc_scaling_parameters
           CALL calc_surface_fluxes
-          IF ( do_output_at_2m )  THEN
-             CALL calc_pt_near_surface ( '2m' )
-          ENDIF
        ENDIF
 !
@@ -272,7 +285,4 @@
           CALL calc_scaling_parameters
           CALL calc_surface_fluxes
-          IF ( do_output_at_2m )  THEN
-             CALL calc_pt_near_surface ( '2m' )
-          ENDIF
        ENDIF
 !
@@ -285,7 +295,4 @@
           CALL calc_scaling_parameters
           CALL calc_surface_fluxes
-          IF ( do_output_at_2m )  THEN
-             CALL calc_pt_near_surface ( '2m' )
-          ENDIF
 !
 !--       Calculate 10cm temperature, required in indoor model
@@ -1863,5 +1870,4 @@
        INTEGER(iwp)                    :: j      !< grid index y-dimension 
        INTEGER(iwp)                    :: k      !< grid index z-dimension
-       INTEGER(iwp)                    :: kk     !< running index along the z-dimension
        INTEGER(iwp)                    :: m      !< running index for surface elements
 
@@ -1882,49 +1888,7 @@
                                      - psi_h( 0.1_wp / surf%ol(m) )            &
                                      + psi_h( surf%z0h(m) / surf%ol(m) ) )
-                                    
-             ENDDO
-
-!
-!--       2-m temperature. Note, this is only calculated for output reasons at
-!--       horizontal upward-facing surfaces. 
-          CASE ( '2m' )
-      
-             DO  m = 1, surf%ns
-
-                i = surf%i(m)
-                j = surf%j(m)
-                k = surf%k(m)
-!
-!--             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
-                   surf%pt_2m(m) = surf%pt_surface(m) + surf%ts(m) / kappa     &
-                                   * ( LOG( 2.0_wp /  surf%z0h(m) )            &
-                                     - psi_h( 2.0_wp / surf%ol(m) )            &
-                                     + psi_h( surf%z0h(m) / surf%ol(m) ) )
-!
-!--             If 2-m level is above the first grid level, 2-m temperature
-!--             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. 
-                ELSE
-!
-!--                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. 
-                   surf%pt_2m(m) = pt(kk-1,j,i) +                              &
-                                 ( zw(k-1) + 2.0_wp - zu(kk-1)     ) *         &
-                                 ( pt(kk,j,i)       - pt(kk-1,j,i) ) /         &
-                                 ( zu(kk)           - zu(kk-1)     )
-                ENDIF
-
-             ENDDO
-          
-       
+
+             ENDDO
+
        END SELECT
 
Index: /palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_mod.f90	(revision 4330)
+++ /palm/trunk/SOURCE/surface_mod.f90	(revision 4331)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! -pt_2m - array is moved to diagnostic_output_quantities
+! 
+! 4329 2019-12-10 15:46:36Z motisi
 ! Renamed wall_flags_0 to wall_flags_static_0
 ! 
@@ -289,5 +292,4 @@
        
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_10cm             !< near surface air potential temperature at distance 10 cm from the surface (K)
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_2m               !< near surface air potential temperature at distance 2 m from the surface (K)     
        
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  alpha_vg          !< coef. of Van Genuchten
@@ -1281,7 +1283,4 @@
 !--    Salinity surface flux
        IF ( ocean_mode )  DEALLOCATE ( surfaces%sasws )
-!
-!--    2-m potential temperature (for output quantity theta_2m*)
-       IF ( do_output_at_2m )  DEALLOCATE ( surfaces%pt_2m )
 
     END SUBROUTINE deallocate_surface_attributes_h
@@ -1412,11 +1411,4 @@
 !--    Salinity surface flux
        IF ( ocean_mode )  ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
-!
-!--    2-m potential temperature (for output quantity theta_2m*)
-       IF ( do_output_at_2m )  THEN
-          ALLOCATE ( surfaces%pt_2m(1:surfaces%ns) )
-          surfaces%pt_2m = -9999.0_wp  !< output array (for theta_2m*) must be initialized here,
-                                       !< otherwise simulation crash at do2d_at_begin with spinup=.F.
-       ENDIF
 
     END SUBROUTINE allocate_surface_attributes_h
Index: /palm/trunk/SOURCE/time_integration_spinup.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration_spinup.f90	(revision 4330)
+++ /palm/trunk/SOURCE/time_integration_spinup.f90	(revision 4331)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Enable output of diagnostic quantities, e.g. 2-m temperature
+! 
+! 4227 2019-09-10 18:04:34Z gronemeier
 ! implement new palm_date_time_mod
 ! 
@@ -81,4 +84,7 @@
     USE cpulog,                                                                &
         ONLY:  cpu_log, log_point_s
+
+    USE diagnostic_output_quantities_mod,                                      &
+        ONLY:  doq_calculate
 
     USE indices,                                                               &
@@ -468,4 +474,5 @@
 !--       2d-data output (cross-sections)
           IF ( time_do2d_xy >= dt_do2d_xy )  THEN
+             CALL doq_calculate
              CALL data_output_2d( 'xy', 0 )
              time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
@@ -475,4 +482,5 @@
 !--       3d-data output (volume data)
           IF ( time_do3d >= dt_do3d )  THEN
+             CALL doq_calculate
              CALL data_output_3d( 0 )
              time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
@@ -560,5 +568,5 @@
 
        USE basic_constants_and_equations_mod,                                  &
-       ONLY:  pi
+           ONLY:  pi
       
        USE kinds
Index: /palm/trunk/SOURCE/write_restart_data_mod.f90
===================================================================
--- /palm/trunk/SOURCE/write_restart_data_mod.f90	(revision 4330)
+++ /palm/trunk/SOURCE/write_restart_data_mod.f90	(revision 4331)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! - Move 2-m potential temperature to diagnostic_output_quantities.
+! 
+! 4301 2019-11-22 12:09:09Z oliver.maas
 ! removed recycling_yshift
 ! 
@@ -942,9 +945,4 @@
        ENDIF
 
-       IF ( ALLOCATED( pt_2m_av ) )  THEN
-          CALL wrd_write_string( 'pt_2m_av' ) 
-          WRITE ( 14 )  pt_2m_av
-       ENDIF
-
        IF ( humidity )  THEN
 
