Index: palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
===================================================================
--- palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4756)
+++ 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 )
Index: palm/trunk/SOURCE/dynamics_mod.f90
===================================================================
--- palm/trunk/SOURCE/dynamics_mod.f90	(revision 4756)
+++ palm/trunk/SOURCE/dynamics_mod.f90	(revision 4757)
@@ -24,4 +24,7 @@
 ! -----------------
 ! $Id$
+! Implement relative humidity as diagnostic output quantity
+!
+! 4731 2020-10-07 13:25:11Z schwenkel
 ! Move exchange_horiz from time_integration to modules
 !
@@ -82,5 +85,6 @@
                u, u_1, u_2, u_init, u_p, u_m_l, u_m_n, u_m_r, u_m_s,                               &
                v, v_1, v_2, v_p, v_init, v_m_l, v_m_n, v_m_r, v_m_s,                               &
-               w, w_1, w_2, w_p, w_m_l, w_m_n, w_m_r, w_m_s
+               w, w_1, w_2, w_p, w_m_l, w_m_n, w_m_r, w_m_s,                                       &
+               zu
 
     USE basic_constants_and_equations_mod,                                                         &
@@ -518,4 +522,8 @@
     REAL(wp)     ::  q_s !< saturation mixing ratio
     REAL(wp)     ::  t_l !< actual temperature
+    REAL(wp)     ::  rh     = 9999999.9_wp !< relative humidity
+    REAL(wp)     ::  rh_min = 9999999.9_wp !< max relative humidity
+    REAL(wp)     ::  height = 9999999.9_wp !< height of supersaturated regions
+    REAL(wp)     ::  min_height = 9999999.9_wp !< height of supersaturated regions
 
 !
@@ -537,5 +545,9 @@
                 q_s = rd_d_rv * e_s / ( hyp(k) - e_s )
 
-                IF ( q(k,j,i) > 1.02_wp * q_s )  realistic_q = .FALSE.
+                IF ( q(k,j,i) > 1.02_wp * q_s )  THEN
+                   realistic_q = .FALSE.
+                   rh = q(k,j,i) / q_s * 100.0_wp
+                   height = zu(k)
+                ENDIF
              ENDDO
           ENDDO
@@ -546,13 +558,17 @@
 #if defined( __parallel )
        CALL MPI_ALLREDUCE( MPI_IN_PLACE, realistic_q, 1, MPI_LOGICAL, MPI_LAND, comm2d, ierr)
+       CALL MPI_ALLREDUCE( rh, rh_min, 1, MPI_REAL, MPI_MIN, comm2d, ierr )
+       CALL MPI_ALLREDUCE( height, min_height, 1, MPI_REAL, MPI_MIN, comm2d, ierr )
 #endif
 
        IF ( .NOT. realistic_q  .AND.                                                               &
             TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
-          message_string = 'The initial mixing ratio exceeds the saturation mixing ratio.'
+            WRITE( message_string, * ) 'The initial mixing ratio exceeds the saturation mixing' // &
+               'ratio, with rh =', rh_min, '% in a height of', min_height, 'm for the first time'
           CALL message( 'dynamic_init_checks', 'PA0697', 2, 2, 0, 6, 0 )
        ELSEIF ( .NOT. realistic_q  .AND.                                                           &
                 TRIM( initializing_actions ) == 'read_restart_data' )  THEN
-          message_string = 'The mixing ratio exceeds the saturation mixing ratio.'
+          WRITE( message_string, * ) 'The initial mixing ratio exceeds the saturation mixing' //   &
+              'ratio, with rh =', rh_min, '% in a height of', min_height, 'm for the first time'
           CALL message( 'dynamic_init_checks', 'PA0697', 0, 1, 0, 6, 0 )
        ENDIF
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 4756)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 4757)
@@ -24,4 +24,7 @@
 ! -----------------
 ! $Id$
+! Implement relative humidity as diagnostic output quantity
+!
+! 4742 2020-10-14 15:11:02Z schwenkel
 ! Implement snow and graupel (bulk microphysics)
 !
@@ -105,5 +108,5 @@
 
     USE arrays_3d,                                                                                 &
-        ONLY:  ddzu, ddzw, d_exner, e, heatflux_output_conversion, hyp, km, kh,                    &
+        ONLY:  ddzu, ddzw, d_exner, e, exner, heatflux_output_conversion, hyp, km, kh,             &
                momentumflux_output_conversion,                                                     &
                nc, ni, ng, nr, ns, p, prho, prr, pt, q, qc, qi, qg, ql, qr, qs,                    &
@@ -112,5 +115,5 @@
 
     USE basic_constants_and_equations_mod,                                                         &
-        ONLY:  g, lv_d_cp
+        ONLY:  g, lv_d_cp, magnus, rd_d_rv
 
     USE bulk_cloud_model_mod,                                                                      &
@@ -220,4 +223,9 @@
     REAL(wp) ::  v2               !<
     REAL(wp) ::  w2               !<
+
+    REAL(wp) ::  temp           !< temperature
+    REAL(wp) ::  e_s            !< saturation vapor pressure
+    REAL(wp) ::  q_s            !< saturation mixing ratio
+    REAL(wp) ::  rh_tmp         !< local value of relative humidity
 
     REAL(wp) ::  dptdz(nzb+1:nzt+1)    !<
@@ -1390,4 +1398,9 @@
                       END IF
                    END IF
+                   temp = exner(k) * pt(k,j,i)
+                   e_s = magnus( temp )
+                   q_s = rd_d_rv * e_s / ( hyp(k) - e_s )
+                   rh_tmp = q(k,j,i) / q_s * 100.0_wp
+                   sums_l(k,130,tn) = sums_l(k,130,tn) + rh_tmp * rmask(j,i,sr) * flag
                 ENDIF
 !
@@ -1747,6 +1760,6 @@
                 IF ( max_pr_cs > 0 )  THEN
                      sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+ max_pr_cs,0) =          &
-                               sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs,0) + &
-                               sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs,i)
+                               sums_l(:,pr_palm+max_pr_user+1:pr_palm+max_pr_user+max_pr_cs,0) + &
+                               sums_l(:,pr_palm+max_pr_user+1:pr_palm+max_pr_user+max_pr_cs,i)
 
                 ENDIF
@@ -1832,5 +1845,5 @@
              sums(k,116)           = sums(k,116)           / ngp_2dh_s_inner(k,sr)
              sums(k,118:pr_palm-2) = sums(k,118:pr_palm-2) / ngp_2dh_s_inner(k,sr)
-             sums(k,123:129)       = sums(k,123:129) * ngp_2dh_s_inner(k,sr)  / ngp_2dh(sr)
+             sums(k,123:130)       = sums(k,123:130) * ngp_2dh_s_inner(k,sr)  / ngp_2dh(sr)
           ENDIF
        ENDDO
@@ -1946,4 +1959,5 @@
        hom(:,1,128,sr) = sums(:,128)   ! ns
        hom(:,1,129,sr) = sums(:,129)   ! qs
+       hom(:,1,130,sr) = sums(:,130)   ! rh
        hom(:,1,73,sr) = sums(:,73)     ! nr
        hom(:,1,74,sr) = sums(:,74)     ! qr
Index: palm/trunk/SOURCE/module_interface.f90
===================================================================
--- palm/trunk/SOURCE/module_interface.f90	(revision 4756)
+++ palm/trunk/SOURCE/module_interface.f90	(revision 4757)
@@ -24,4 +24,7 @@
 ! -----------------
 ! $Id$
+! Implement doq calls
+!
+! 4753 2020-10-21 14:55:41Z raasch
 ! file re-formatted to follow the PALM coding standard
 !
@@ -342,8 +345,10 @@
         ONLY:  doq_3d_data_averaging,                                                              &
                doq_check_data_output,                                                              &
+               doq_check_data_output_pr,                                                           &
                doq_define_netcdf_grid,                                                             &
                doq_init,                                                                           &
                doq_output_2d,                                                                      &
                doq_output_3d,                                                                      &
+               doq_statistics,                                                                     &
                doq_rrd_local,                                                                      &
                doq_wrd_local
@@ -870,4 +875,5 @@
 
     CALL dynamics_check_data_output_pr( variable, var_count, unit, dopr_unit )
+    CALL doq_check_data_output_pr( variable, var_count, unit, dopr_unit )
 
     IF ( unit == 'illegal'  .AND.  bulk_cloud_model )  THEN
@@ -1701,4 +1707,5 @@
 
     CALL dynamics_statistics( mode, sr, tn )
+    CALL doq_statistics( mode, sr, tn )
 
     IF ( gust_module_enabled )  CALL gust_statistics( mode, sr, tn, dots_max )
