Index: palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
===================================================================
--- palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4329)
+++ 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
