Index: palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
===================================================================
--- palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4860)
+++ palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4861)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! ------------------
-!
+! Implemented air temperature as diagnostic output quantity
 !
 ! Former revisions:
@@ -31,5 +31,5 @@
 ! 4671 2020-09-09 20:27:58Z pavelkrc
 ! Implementation of downward facing USM and LSM surfaces
-! 
+!
 ! 4583 2020-06-29 12:36:47Z raasch
 ! file re-formatted to follow the PALM coding standard
@@ -127,5 +127,6 @@
 
     USE basic_constants_and_equations_mod,                                                         &
-        ONLY:  kappa,                                                                              &
+        ONLY:  degc_to_k,                                                                          &
+               kappa,                                                                              &
                pi,                                                                                 &
                magnus,                                                                             &
@@ -203,4 +204,6 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  rh     !< relative humidity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  rh_av  !< avg. relative humidity
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ta     !< air temperature
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ta_av  !< avg. air temperature
     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
@@ -339,4 +342,10 @@
              rh_av = 0.0_wp
 
+          CASE ( 'ta' )
+             IF ( .NOT. ALLOCATED( ta_av ) )  THEN
+                ALLOCATE( ta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             ENDIF
+             ta_av = 0.0_wp
+
           CASE ( 'ti' )
              IF ( .NOT. ALLOCATED( ti_av ) )  THEN
@@ -430,4 +439,15 @@
                       DO  k = nzb, nzt+1
                          rh_av(k,j,i) = rh_av(k,j,i) + rh(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          CASE ( 'ta' )
+             IF ( ALLOCATED( ta_av ) )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         ta_av(k,j,i) = ta_av(k,j,i) + ta(k,j,i)
                       ENDDO
                    ENDDO
@@ -584,4 +604,15 @@
              ENDIF
 
+          CASE ( 'ta' )
+             IF ( ALLOCATED( ta_av ) )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         ta_av(k,j,i) = ta_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
           CASE ( 'ti' )
              IF ( ALLOCATED( ti_av ) )  THEN
@@ -752,4 +783,7 @@
           unit = '%'
 
+       CASE ( 'ta' )
+          unit = 'degree_C'
+
        CASE ( 'ti' )
           unit = '1/s'
@@ -868,4 +902,5 @@
 !--    s grid
        CASE ( 'rh', 'rh_xy', 'rh_xz', 'rh_yz',                                                     &
+              'ta', 'ta_xy', 'ta_xz', 'ta_yz',                                                     &
               'ti', 'ti_xy', 'ti_xz', 'ti_yz',                                                     &
               'wspeed', 'wspeed_xy', 'wspeed_xz', 'wspeed_yz',                                     &
@@ -974,4 +1009,18 @@
            IF ( mode == 'xy' )  grid = 'zu'
 
+       CASE ( 'ta_xy', 'ta_xz', 'ta_yz' )
+           IF ( av == 0 )  THEN
+              to_be_resorted => ta
+           ELSE
+              IF ( .NOT. ALLOCATED( ta_av ) )  THEN
+                 ALLOCATE( ta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                 ta_av = REAL( fill_value, KIND = wp )
+              ENDIF
+              to_be_resorted => ta_av
+           ENDIF
+           flag_nr = 0
+
+           IF ( mode == 'xy' )  grid = 'zu'
+
        CASE ( 'ti_xy', 'ti_xz', 'ti_yz' )
            IF ( av == 0 )  THEN
@@ -1226,4 +1275,16 @@
           flag_nr = 0
 
+       CASE ( 'ta' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => ta
+          ELSE
+             IF ( .NOT. ALLOCATED( ta ) )  THEN
+                ALLOCATE( ta(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ta = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => ta
+          ENDIF
+          flag_nr = 0
+
        CASE ( 'ti' )
           IF ( av == 0 )  THEN
@@ -1417,4 +1478,13 @@
           ELSE
              to_be_resorted => rh_av
+          ENDIF
+          grid = 's'
+          flag_nr = 0
+
+       CASE ( 'ta' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => ta
+          ELSE
+             to_be_resorted => ta_av
           ENDIF
           grid = 's'
@@ -1597,4 +1667,11 @@
              ENDIF
 !
+!--       Allocate array for 'air temperature'
+          CASE ( 'ta' )
+             IF ( .NOT. ALLOCATED( ta ) )  THEN
+                ALLOCATE( ta(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ta = 0.0_wp
+             ENDIF
+!
 !--       Allocate array for 'turbulence intensity'
           CASE ( 'ti' )
@@ -1776,4 +1853,14 @@
                                      * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
                       ENDIF
+                   ENDDO
+                ENDDO
+             ENDDO
+!
+!--       ta (air temperature)
+          CASE ( 'ta' )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb+1, nzt
+                      ta(k,j,i) = exner(k) * pt(k,j,i) - degc_to_k
                    ENDDO
                 ENDDO
@@ -2200,4 +2287,12 @@
                       tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
+       CASE ( 'ta_av' )
+          IF ( .NOT. ALLOCATED( ta_av ) )  THEN
+             ALLOCATE( ta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( k == 1 )  READ ( 13 )  tmp_3d
+          ta_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
+                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
        CASE ( 'ti_av' )
           IF ( .NOT. ALLOCATED( ti_av ) )  THEN
@@ -2321,4 +2416,10 @@
     ENDIF
 
+    CALL rd_mpi_io_check_array( 'ta_av' , found = array_found )
+    IF ( array_found )  THEN
+       IF ( .NOT. ALLOCATED( ta_av ) )  ALLOCATE( ta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+       CALL rrd_mpi_io( 'ta_av', ta_av )
+    ENDIF
+
     CALL rd_mpi_io_check_array( 'ti_av' , found = array_found )
     IF ( array_found )  THEN
@@ -2423,4 +2524,9 @@
        ENDIF
 
+       IF ( ALLOCATED( ta_av ) )  THEN
+          CALL wrd_write_string( 'ta_av' )
+          WRITE ( 14 )  ta_av
+       ENDIF
+
        IF ( ALLOCATED( ti_av ) )  THEN
           CALL wrd_write_string( 'ti_av' )
@@ -2487,4 +2593,5 @@
        IF ( ALLOCATED( pt_2m_av ) )     CALL wrd_mpi_io( 'pt_2m_av', pt_2m_av )
        IF ( ALLOCATED( rh_av ) )        CALL wrd_mpi_io( 'rh_av', rh_av )
+       IF ( ALLOCATED( ta_av ) )        CALL wrd_mpi_io( 'ta_av', ta_av )
        IF ( ALLOCATED( ti_av ) )        CALL wrd_mpi_io( 'ti_av', ti_av )
        IF ( ALLOCATED( u_center_av ) )  CALL wrd_mpi_io( 'u_center_av', u_center_av )
