Index: palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
===================================================================
--- palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4671)
+++ palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4757)
@@ -24,4 +24,7 @@
 ! -----------------
 ! $Id$
+! Implement relative humidity as diagnostic output quantity
+!
+! 4671 2020-09-09 20:27:58Z pavelkrc
 ! Implementation of downward facing USM and LSM surfaces
 ! 
@@ -110,4 +113,6 @@
     USE arrays_3d,                                                                                 &
         ONLY:  ddzu,                                                                               &
+               exner,                                                                              &
+               hyp,                                                                                &
                pt,                                                                                 &
                q,                                                                                  &
@@ -119,9 +124,14 @@
 
     USE basic_constants_and_equations_mod,                                                         &
-        ONLY:  kappa, pi
+        ONLY:  kappa,                                                                              &
+               pi,                                                                                 &
+               magnus,                                                                             &
+               rd_d_rv
 
     USE control_parameters,                                                                        &
         ONLY:  current_timestep_number,                                                            &
                data_output,                                                                        &
+               data_output_pr,                                                                     &
+               humidity,                                                                           &
                message_string,                                                                     &
                restart_data_format_output,                                                         &
@@ -153,8 +163,18 @@
     USE kinds
 
+    USE profil_parameter,                                                                          &
+        ONLY:  dopr_index
+
     USE restart_data_mpi_io_mod,                                                                   &
         ONLY:  rd_mpi_io_check_array,                                                              &
                rrd_mpi_io,                                                                         &
                wrd_mpi_io
+
+    USE statistics,                                                                                &
+        ONLY:  hom,                                                                                &
+               pr_palm,                                                                            &
+               rmask,                                                                              &
+               statistic_regions,                                                                  &
+               sums_l
 
     USE surface_mod,                                                                               &
@@ -164,5 +184,4 @@
                surf_usm_h
 
-
     IMPLICIT NONE
 
@@ -179,4 +198,6 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  uv_10m_av !< averaged horizontal wind speed at 10m
 
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  rh     !< relative humidity
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  rh_av  !< avg. relative humidity
     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
@@ -220,4 +241,5 @@
            doq_calculate,                                                                          &
            doq_check_data_output,                                                                  &
+           doq_check_data_output_pr,                                                               &
            doq_define_netcdf_grid,                                                                 &
            doq_init,                                                                               &
@@ -225,4 +247,5 @@
            doq_output_3d,                                                                          &
            doq_output_mask,                                                                        &
+           doq_statistics,                                                                         &
            doq_rrd_local,                                                                          &
            doq_wrd_local
@@ -241,4 +264,8 @@
     END INTERFACE doq_check_data_output
 
+    INTERFACE doq_check_data_output_pr
+       MODULE PROCEDURE doq_check_data_output_pr
+    END INTERFACE doq_check_data_output_pr
+
     INTERFACE doq_define_netcdf_grid
        MODULE PROCEDURE doq_define_netcdf_grid
@@ -260,4 +287,8 @@
        MODULE PROCEDURE doq_init
     END INTERFACE doq_init
+
+    INTERFACE doq_statistics
+       MODULE PROCEDURE doq_statistics
+    END INTERFACE doq_statistics
 
     INTERFACE doq_prepare
@@ -274,5 +305,4 @@
     END INTERFACE doq_wrd_local
 
-
  CONTAINS
 
@@ -300,4 +330,10 @@
        SELECT CASE ( TRIM( variable ) )
 
+          CASE ( 'rh' )
+             IF ( .NOT. ALLOCATED( rh_av ) )  THEN
+                ALLOCATE( rh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             ENDIF
+             rh_av = 0.0_wp
+
           CASE ( 'ti' )
              IF ( .NOT. ALLOCATED( ti_av ) )  THEN
@@ -384,4 +420,15 @@
 
        SELECT CASE ( TRIM( variable ) )
+
+          CASE ( 'rh' )
+             IF ( ALLOCATED( rh_av ) )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         rh_av(k,j,i) = rh_av(k,j,i) + rh(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
 
           CASE ( 'ti' )
@@ -523,4 +570,15 @@
        SELECT CASE ( TRIM( variable ) )
 
+          CASE ( 'rh' )
+             IF ( ALLOCATED( rh_av ) )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         rh_av(k,j,i) = rh_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
           CASE ( 'ti' )
              IF ( ALLOCATED( ti_av ) )  THEN
@@ -688,4 +746,7 @@
     SELECT CASE ( TRIM( var ) )
 
+       CASE ( 'rh' )
+          unit = '%'
+
        CASE ( 'ti' )
           unit = '1/s'
@@ -740,4 +801,44 @@
  END SUBROUTINE doq_check_data_output
 
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Set the unit of user defined profile output quantities. For those variables not recognized by the
+!> user, the parameter unit is set to "illegal", which tells the calling routine that the
+!> output variable is not defined and leads to a program abort.
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE doq_check_data_output_pr( variable, var_count, unit, dopr_unit )
+
+    CHARACTER (LEN=*) ::  unit      !<
+    CHARACTER (LEN=*) ::  variable  !<
+    CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
+
+    INTEGER(iwp) ::  pr_index       !<
+    INTEGER(iwp) ::  var_count      !<
+
+    SELECT CASE ( TRIM( variable ) )
+
+       CASE ( 'rh' )
+          IF ( .NOT.  humidity )  THEN
+             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) //       &
+                              ' requires humidity'
+             CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
+          ENDIF
+          pr_index = 130
+          dopr_index(var_count) = pr_index
+          dopr_unit     = '%'
+          unit = dopr_unit
+          hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
+
+       CASE DEFAULT
+          unit = 'illegal'
+
+    END SELECT
+
+
+ END SUBROUTINE doq_check_data_output_pr
+
+
 !--------------------------------------------------------------------------------------------------!
 !
@@ -763,5 +864,6 @@
 !
 !--    s grid
-       CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz',                                                     &
+       CASE ( 'rh', 'rh_xy', 'rh_xz', 'rh_yz',                                                     &
+              'ti', 'ti_xy', 'ti_xz', 'ti_yz',                                                     &
               'wspeed', 'wspeed_xy', 'wspeed_xz', 'wspeed_yz',                                     &
               'wdir', 'wdir_xy', 'wdir_xz', 'wdir_yz' )
@@ -854,4 +956,18 @@
 
     SELECT CASE ( TRIM( variable ) )
+
+       CASE ( 'rh_xy', 'rh_xz', 'rh_yz' )
+           IF ( av == 0 )  THEN
+              to_be_resorted => rh
+           ELSE
+              IF ( .NOT. ALLOCATED( rh_av ) )  THEN
+                 ALLOCATE( rh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                 rh_av = REAL( fill_value, KIND = wp )
+              ENDIF
+              to_be_resorted => rh_av
+           ENDIF
+           flag_nr = 0
+
+           IF ( mode == 'xy' )  grid = 'zu'
 
        CASE ( 'ti_xy', 'ti_xz', 'ti_yz' )
@@ -1095,4 +1211,16 @@
     SELECT CASE ( TRIM( variable ) )
 
+       CASE ( 'rh' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => rh
+          ELSE
+             IF ( .NOT. ALLOCATED( rh ) )  THEN
+                ALLOCATE( rh(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                rh = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => rh
+          ENDIF
+          flag_nr = 0
+
        CASE ( 'ti' )
           IF ( av == 0 )  THEN
@@ -1280,4 +1408,13 @@
 
     SELECT CASE ( TRIM( variable ) )
+
+       CASE ( 'rh' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => rh
+          ELSE
+             to_be_resorted => rh_av
+          ENDIF
+          grid = 's'
+          flag_nr = 0
 
        CASE ( 'ti' )
@@ -1450,4 +1587,11 @@
        SELECT CASE ( TRIM( do_all(ivar) ) )
 !
+!--       Allocate array for 'relative humidity'
+          CASE ( 'rh' )
+             IF ( .NOT. ALLOCATED( rh ) )  THEN
+                ALLOCATE( rh(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                rh = 0.0_wp
+             ENDIF
+!
 !--       Allocate array for 'turbulence intensity'
           CASE ( 'ti' )
@@ -1556,4 +1700,33 @@
 ! Description:
 ! ------------
+!> Calculation of module-specific statistics, i.e. horizontally averaged profiles and time series.
+!> This is called for every statistic region sr, but at least for the region "total domain" (sr=0).
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE doq_statistics( mode, sr, tn )
+
+
+    CHARACTER (LEN=*) ::  mode   !<
+
+!     INTEGER(iwp) ::  i    !<
+!     INTEGER(iwp) ::  j    !<
+!     INTEGER(iwp) ::  k    !<
+    INTEGER(iwp) ::  sr   !<
+    INTEGER(iwp) ::  tn   !<
+!
+!-- Next line is to avoid compiler warning about unused variables. Please remove.
+    IF ( sr == 0  .OR.  tn == 0 )  CONTINUE
+
+    IF ( mode == 'profiles' )  THEN
+
+    ELSEIF ( mode == 'time_series' )  THEN
+
+    ENDIF
+
+ END SUBROUTINE doq_statistics
+
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Calculate diagnostic quantities
 !--------------------------------------------------------------------------------------------------!
@@ -1567,4 +1740,8 @@
     INTEGER(iwp) ::  ivar       !< loop index over all 2d/3d/mask output quantities
 
+    REAL(wp) ::  temp           !< temperature
+    REAL(wp) ::  e_s            !< saturation vapor pressure
+    REAL(wp) ::  q_s            !< saturation mixing ratio
+
     TYPE(surf_type), POINTER ::  surf     !< surf-type array, used to generalize subroutines
 
@@ -1583,4 +1760,20 @@
 
        SELECT CASE ( TRIM( do_all(ivar) ) )
+!
+!--       rh (relative humidity)
+          CASE ( 'rh' )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb+1, nzt
+                      IF ( humidity )  THEN
+                         temp = exner(k) * pt(k,j,i)
+                         e_s = magnus( temp )
+                         q_s = rd_d_rv * e_s / ( hyp(k) - e_s )
+                         rh(k,j,i) = q(k,j,i) / q_s * 100.0_wp                                     &
+                                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDDO
 !
 !--       Calculate 'turbulence intensity' from rot[(u,v,w)] at scalar grid point
@@ -1996,4 +2189,12 @@
                       tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
+       CASE ( 'rh_av' )
+          IF ( .NOT. ALLOCATED( rh_av ) )  THEN
+             ALLOCATE( rh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( k == 1 )  READ ( 13 )  tmp_3d
+          rh_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
@@ -2123,4 +2324,10 @@
     ENDIF
 
+    CALL rd_mpi_io_check_array( 'rh_av' , found = array_found )
+    IF ( array_found )  THEN
+       IF ( .NOT. ALLOCATED( rh_av ) )  ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+       CALL rrd_mpi_io( 'rh_av', rh_av )
+    ENDIF
+
     CALL rd_mpi_io_check_array( 'u_center_av' , found = array_found )
     IF ( array_found )  THEN
@@ -2208,4 +2415,9 @@
        ENDIF
 
+       IF ( ALLOCATED( rh_av ) )  THEN
+          CALL wrd_write_string( 'rh_av' )
+          WRITE ( 14 )  rh_av
+       ENDIF
+
        IF ( ALLOCATED( ti_av ) )  THEN
           CALL wrd_write_string( 'ti_av' )
@@ -2271,4 +2483,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( 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 )
