Index: /palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- /palm/trunk/SOURCE/advec_ws.f90	(revision 4580)
+++ /palm/trunk/SOURCE/advec_ws.f90	(revision 4581)
@@ -24,4 +24,7 @@
 ! -----------------
 ! $Id$
+! Enable output of resolved-scale vertical fluxes of chemical species.
+! 
+! 4509 2020-04-26 15:57:55Z raasch
 ! file re-formatted to follow the PALM coding standard
 !
@@ -229,16 +232,28 @@
 
     USE statistics,                                                                                &
-        ONLY:  hom, weight_substep,                                                                &
-               sums_salsa_ws_l, sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l,                       &
-               sums_wsncs_ws_l, sums_wsnrs_ws_l,                                                   &
-               sums_wspts_ws_l, sums_wsqs_ws_l, sums_wsss_ws_l,                                    &
-               sums_wsqcs_ws_l, sums_wsqrs_ws_l,                                                   &
-               sums_wsqis_ws_l, sums_wsnis_ws_l,                                                   &
-               sums_wssas_ws_l, sums_wsus_ws_l, sums_wsvs_ws_l
+        ONLY:  hom,                                                                                &
+               sums_salsa_ws_l,                                                                    &
+               sums_us2_ws_l,                                                                      &
+               sums_vs2_ws_l,                                                                      &
+               sums_ws2_ws_l,                                                                      &
+               sums_wschs_ws_l,                                                                    &
+               sums_wsncs_ws_l,                                                                    &
+               sums_wsnrs_ws_l,                                                                    &
+               sums_wspts_ws_l,                                                                    &
+               sums_wsqs_ws_l,                                                                     &
+               sums_wsss_ws_l,                                                                     &
+               sums_wsqcs_ws_l,                                                                    &
+               sums_wsqrs_ws_l,                                                                    &
+               sums_wsqis_ws_l,                                                                    &
+               sums_wsnis_ws_l,                                                                    &
+               sums_wssas_ws_l,                                                                    &
+               sums_wsus_ws_l,                                                                     &
+               sums_wsvs_ws_l,                                                                     &
+               weight_substep
 
 
 
     IMPLICIT NONE
-    
+
     REAL(wp) ::  adv_mom_1            !< 1/4 - constant used in 5th-order advection scheme for momentum advection (1st-order part)
     REAL(wp) ::  adv_mom_3            !< 1/24 - constant used in 5th-order advection scheme for momentum advection (3rd-order part)
@@ -1875,5 +1890,5 @@
           ENDDO
 
-      CASE ( 'aerosol_mass', 'aerosol_number', 'salsa_gas' )
+       CASE ( 'aerosol_mass', 'aerosol_number', 'salsa_gas' )
 
           DO  k = nzb, nzt
@@ -1886,6 +1901,13 @@
           ENDDO
 
-!       CASE ( 'kc' )
-        !kk Has to be implemented for kpp chemistry
+       CASE ( 'kc' )
+          DO  k = nzb, nzt
+             sums_wschs_ws_l(k,tn)  = sums_wschs_ws_l(k,tn) +                                      &
+                                      ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )    &
+                                                  * ( w(k,j,i) - hom(k,1,3,0)                 )    &
+                                      + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )    &
+                                                  *   ABS( w(k,j,i) - hom(k,1,3,0)            )    &
+                                      ) * weight_substep(intermediate_timestep_count)
+          ENDDO
 
     END SELECT
Index: /palm/trunk/SOURCE/chem_modules.f90
===================================================================
--- /palm/trunk/SOURCE/chem_modules.f90	(revision 4580)
+++ /palm/trunk/SOURCE/chem_modules.f90	(revision 4581)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! Enable output of vertical fluxes of chemical species.
+! 
+! 4577 2020-06-25 09:53:58Z raasch
 ! further re-formatting concerning Fortran parameter variables
 !
@@ -163,6 +166,7 @@
                                             !< 4: non-cyclic
 
-    INTEGER(iwp) ::  cs_pr_count                           = 0      !< counter for chemical species
-                                                                    !< profiles
+    INTEGER(iwp) ::  cs_pr_count_fl_res                    = 0      !< counter for vertical flux profiles of chemical species (resolved-scale)
+    INTEGER(iwp) ::  cs_pr_count_fl_sgs                    = 0      !< counter for vertical flux profiles of chemical species (SGS)
+    INTEGER(iwp) ::  cs_pr_count_sp                        = 0      !< counter for chemical species profiles
     INTEGER(iwp) ::  cs_vertical_gradient_level_ind(99,10) = -9999  !< grid index values of
                                                                     !< cs_vertical_gradient_level
@@ -174,5 +178,5 @@
     INTEGER(iwp) ::  main_street_id                        = 0      !< namelist parameter: lower bound of main street IDs
                                                                     !< (OpenStreetMaps) for PARAMETERIZED mode
-    INTEGER(iwp) ::  max_pr_cs                             = 0      !<
+    INTEGER(iwp) ::  max_pr_cs                             = 0      !< number of chemistry profiles in output
     INTEGER(iwp) ::  max_street_id                         = 0      !< namelist parameter: upper bound of main street IDs
                                                                     !< (OpenStreetMaps) for PARAMETERIZED mode
@@ -182,11 +186,13 @@
                                                                     !< (OpenStreetMaps) for PARAMETERIZED mode
 
-    INTEGER(iwp), DIMENSION(99) ::  cs_pr_index  = 0      !< index for chemical species profiles
-    INTEGER(iwp), DIMENSION(:)  ::  match_spec_nox(1:2)   !< results of matching the input and
-                                                          !< model's NOx
-    INTEGER(iwp), DIMENSION(:)  ::  match_spec_pm(1:3)    !< results of matching the input and
-                                                          !< model's PMs
-    INTEGER(iwp), DIMENSION(:)  ::  match_spec_sox(1:2)   !< results of matching the input and
-                                                          !< model's SOx!
+    INTEGER(iwp), DIMENSION(99) ::  cs_pr_index_sp      = 0 !< index for chemical species mean profiles
+    INTEGER(iwp), DIMENSION(99) ::  cs_pr_index_fl_res  = 0 !< index for chemical species sgs-flux profiles
+    INTEGER(iwp), DIMENSION(99) ::  cs_pr_index_fl_sgs  = 0 !< index for chemical species resolved-scale flux profiles
+    INTEGER(iwp), DIMENSION(99) ::  hom_index_fl_res    = 0 !< index of the resolved-scale flux profile with respect to the hom array
+    INTEGER(iwp), DIMENSION(99) ::  hom_index_fl_sgs    = 0 !< index of the SGS flux profile with respect to the hom array
+    INTEGER(iwp), DIMENSION(99) ::  hom_index_spec      = 0 !< index of the profile with respect to the hom array
+    INTEGER(iwp), DIMENSION(:)  ::  match_spec_nox(1:2)     !< results of matching the input and model's NOx
+    INTEGER(iwp), DIMENSION(:)  ::  match_spec_pm(1:3)      !< results of matching the input and model's PMs
+    INTEGER(iwp), DIMENSION(:)  ::  match_spec_sox(1:2)     !< results of matching the input and model's SOx!
 
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  match_spec_input      !< index of input chem species
@@ -262,6 +268,7 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  css          !< scaling parameter for chem species
 
-    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  emis_distribution  !< emissions final values
-                                                                     !< (main module output) ???
+    REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE, TARGET ::  sums_ws_l  !< subdomain sum of vertical chemistry flux w'ch'
+
+    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  emis_distribution  !< emissions final values (main module output) ???
 
     REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  cs_1  !< pointer for swapping of
Index: /palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 4580)
+++ /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 4581)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! Enable output of resolved-scale vertical fluxes of chemical species.
+! 
+! 4577 2020-06-25 09:53:58Z raasch
 ! further re-formatting to follow the PALM coding standard
 !
@@ -436,68 +439,68 @@
 !
 !-- Interface section
-    INTERFACE chem_3d_data_averaging
-       MODULE PROCEDURE chem_3d_data_averaging
-    END INTERFACE chem_3d_data_averaging
-
-    INTERFACE chem_boundary_conditions
-       MODULE PROCEDURE chem_boundary_conditions
-    END INTERFACE chem_boundary_conditions
-
-    INTERFACE chem_check_data_output
-       MODULE PROCEDURE chem_check_data_output
-    END INTERFACE chem_check_data_output
-
-    INTERFACE chem_data_output_2d
-       MODULE PROCEDURE chem_data_output_2d
-    END INTERFACE chem_data_output_2d
-
-    INTERFACE chem_data_output_3d
-       MODULE PROCEDURE chem_data_output_3d
-    END INTERFACE chem_data_output_3d
-
-    INTERFACE chem_data_output_mask
-       MODULE PROCEDURE chem_data_output_mask
-    END INTERFACE chem_data_output_mask
-
-    INTERFACE chem_check_data_output_pr
-       MODULE PROCEDURE chem_check_data_output_pr
-    END INTERFACE chem_check_data_output_pr
-
-    INTERFACE chem_check_parameters
-       MODULE PROCEDURE chem_check_parameters
-    END INTERFACE chem_check_parameters
-
-    INTERFACE chem_define_netcdf_grid
-       MODULE PROCEDURE chem_define_netcdf_grid
-    END INTERFACE chem_define_netcdf_grid
-
-    INTERFACE chem_header
-       MODULE PROCEDURE chem_header
-    END INTERFACE chem_header
-
-    INTERFACE chem_init_arrays
-       MODULE PROCEDURE chem_init_arrays
-    END INTERFACE chem_init_arrays
-
-    INTERFACE chem_init
-       MODULE PROCEDURE chem_init
-    END INTERFACE chem_init
-
-    INTERFACE chem_init_profiles
-       MODULE PROCEDURE chem_init_profiles
-    END INTERFACE chem_init_profiles
-
-    INTERFACE chem_integrate
-       MODULE PROCEDURE chem_integrate_ij
-    END INTERFACE chem_integrate
-
-    INTERFACE chem_parin
-       MODULE PROCEDURE chem_parin
-    END INTERFACE chem_parin
-
     INTERFACE chem_actions
        MODULE PROCEDURE chem_actions
        MODULE PROCEDURE chem_actions_ij
     END INTERFACE chem_actions
+
+    INTERFACE chem_3d_data_averaging
+       MODULE PROCEDURE chem_3d_data_averaging
+    END INTERFACE chem_3d_data_averaging
+
+    INTERFACE chem_boundary_conditions
+       MODULE PROCEDURE chem_boundary_conditions
+    END INTERFACE chem_boundary_conditions
+
+    INTERFACE chem_check_data_output
+       MODULE PROCEDURE chem_check_data_output
+    END INTERFACE chem_check_data_output
+
+    INTERFACE chem_data_output_2d
+       MODULE PROCEDURE chem_data_output_2d
+    END INTERFACE chem_data_output_2d
+
+    INTERFACE chem_data_output_3d
+       MODULE PROCEDURE chem_data_output_3d
+    END INTERFACE chem_data_output_3d
+
+    INTERFACE chem_data_output_mask
+       MODULE PROCEDURE chem_data_output_mask
+    END INTERFACE chem_data_output_mask
+
+    INTERFACE chem_check_data_output_pr
+       MODULE PROCEDURE chem_check_data_output_pr
+    END INTERFACE chem_check_data_output_pr
+
+    INTERFACE chem_check_parameters
+       MODULE PROCEDURE chem_check_parameters
+    END INTERFACE chem_check_parameters
+
+    INTERFACE chem_define_netcdf_grid
+       MODULE PROCEDURE chem_define_netcdf_grid
+    END INTERFACE chem_define_netcdf_grid
+
+    INTERFACE chem_header
+       MODULE PROCEDURE chem_header
+    END INTERFACE chem_header
+
+    INTERFACE chem_init_arrays
+       MODULE PROCEDURE chem_init_arrays
+    END INTERFACE chem_init_arrays
+
+    INTERFACE chem_init
+       MODULE PROCEDURE chem_init
+    END INTERFACE chem_init
+
+    INTERFACE chem_init_profiles
+       MODULE PROCEDURE chem_init_profiles
+    END INTERFACE chem_init_profiles
+
+    INTERFACE chem_integrate
+       MODULE PROCEDURE chem_integrate_ij
+    END INTERFACE chem_integrate
+
+    INTERFACE chem_parin
+       MODULE PROCEDURE chem_parin
+    END INTERFACE chem_parin
 
     INTERFACE chem_non_advective_processes
@@ -1069,34 +1072,110 @@
 
 
-    CHARACTER (LEN=*)  ::  unit     !<
-    CHARACTER (LEN=*)  ::  variable !<
-    CHARACTER (LEN=*)  ::  dopr_unit
-    CHARACTER (LEN=16) ::  spec_name
-
-    INTEGER(iwp) ::  var_count, lsp  !<
-
-
-    spec_name = TRIM( variable(4:) )
-
-       IF ( .NOT. air_chemistry )  THEN
-          message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) //             &
-                           ' is not imp' // 'lemented for air_chemistry = .FALSE.'
-          CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 )
-
-       ELSE
+    CHARACTER (LEN=*)  ::  unit      !< unit
+    CHARACTER (LEN=*)  ::  variable  !< variable name
+    CHARACTER (LEN=*)  ::  dopr_unit !< unit
+    CHARACTER (LEN=16) ::  spec_name !< species name extracted from output string
+
+    INTEGER(iwp) ::  index_start     !< start index of the species name in data-output string
+    INTEGER(iwp) ::  index_end       !< end index of the species name in data-output string
+    INTEGER(iwp) ::  var_count       !< number of data-output quantity
+    INTEGER(iwp) ::  lsp             !< running index over species
+
+    SELECT CASE ( TRIM( variable(1:3 ) ) )
+
+       CASE ( 'kc_' )
+
+          IF ( .NOT. air_chemistry )  THEN
+             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) //          &
+                              ' is not implemented for air_chemistry = .FALSE.'
+             CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 )
+
+          ENDIF
+!
+!--       Output of total fluxes is not allowed to date.
+          IF ( TRIM( variable(1:4) ) == 'kc_w' )  THEN
+             IF ( TRIM( variable(1:5) ) /= 'kc_w*'  .AND.  TRIM( variable(1:5) ) /= 'kc_w"' )  THEN
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) //       &
+                              ' is currently not implemented. Please output resolved- and '//      &
+                              'subgrid-scale fluxes individually to obtain the total flux.'
+                CALL message( 'chem_check_parameters', 'CM0498', 1, 2, 0, 6, 0 )
+             ENDIF
+          ENDIF
+!
+!--       Check for profile output of first-order moments, i.e. variable(4:) equals a species name.
+          spec_name = TRIM( variable(4:) )
           DO  lsp = 1, nspec
              IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
-                cs_pr_count = cs_pr_count+1
-                cs_pr_index(cs_pr_count) = lsp
-                dopr_index(var_count) = pr_palm+cs_pr_count
-                dopr_unit  = 'ppm'
+                cs_pr_count_sp                 = cs_pr_count_sp + 1
+                cs_pr_index_sp(cs_pr_count_sp) = lsp
+                dopr_index(var_count)          = pr_palm + cs_pr_count_sp + cs_pr_count_fl_sgs + cs_pr_count_fl_res
+                dopr_unit                      = 'ppm'
                 IF ( spec_name(1:2) == 'PM')  THEN
-                     dopr_unit = 'kg m-3'
+                   dopr_unit = 'kg m-3'
                 ENDIF
-                   hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
-                   unit = dopr_unit
+                hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
+                unit                              = dopr_unit
+
+                hom_index_spec(cs_pr_count_sp)    = dopr_index(var_count)
              ENDIF
           ENDDO
-       ENDIF
+!
+!--       Check for profile output of fluxes. variable(index_start:index_end) equals a species name.
+!--       Start with SGS components.
+          IF ( TRIM( variable(1:5) ) == 'kc_w"' )  THEN
+             DO  lsp = 1, nspec
+                index_start = INDEX( TRIM( variable ), TRIM( chem_species(lsp)%name ) )
+                index_end   = LEN( TRIM( variable ) ) - 1
+                IF ( index_start /= 0 )  THEN
+                   spec_name = TRIM( variable(index_start:index_end) )
+                   IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
+                      cs_pr_count_fl_sgs                     = cs_pr_count_fl_sgs + 1
+                      cs_pr_index_fl_sgs(cs_pr_count_fl_sgs) = lsp
+                      dopr_index(var_count)                  = pr_palm + cs_pr_count_sp +          &
+                                                               cs_pr_count_fl_sgs +                &
+                                                               cs_pr_count_fl_res
+                      dopr_unit                              = 'm ppm s-1'
+                      IF ( spec_name(1:2) == 'PM')  THEN
+                         dopr_unit = 'kg m-2 s-1'
+                      ENDIF
+                      hom(:,2, dopr_index(var_count),:)    = SPREAD( zu, 2, statistic_regions+1 )
+                      unit                                 = dopr_unit
+
+                      hom_index_fl_sgs(cs_pr_count_fl_sgs) = dopr_index(var_count)
+                   ENDIF
+                ENDIF
+             ENDDO
+          ENDIF
+!
+!--       Proceed with resolved-scale fluxes.
+          IF ( TRIM( variable(1:5) ) == 'kc_w*' )  THEN
+             spec_name = TRIM( variable(6:) )
+             DO  lsp = 1, nspec
+                index_start = INDEX( TRIM( variable ), TRIM( chem_species(lsp)%name ) )
+                index_end   = LEN( TRIM( variable ) ) - 1
+                IF ( index_start /= 0 )  THEN
+                   spec_name = TRIM( variable(index_start:index_end) )
+                   IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
+                      cs_pr_count_fl_res                     = cs_pr_count_fl_res + 1
+                      cs_pr_index_fl_res(cs_pr_count_fl_res) = lsp
+                      dopr_index(var_count)                  = pr_palm + cs_pr_count_sp +          &
+                                                               cs_pr_count_fl_sgs +                &
+                                                               cs_pr_count_fl_res
+                      dopr_unit                              = 'm ppm s-1'
+                      IF ( spec_name(1:2) == 'PM')  THEN
+                         dopr_unit = 'kg m-2 s-1'
+                      ENDIF
+                      hom(:,2, dopr_index(var_count),:)    = SPREAD( zu, 2, statistic_regions+1 )
+                      unit                                 = dopr_unit
+
+                      hom_index_fl_res(cs_pr_count_fl_res) = dopr_index(var_count)
+                   ENDIF
+                ENDIF
+             ENDDO
+          ENDIF
+       CASE DEFAULT
+          unit = 'illegal'
+
+    END SELECT
 
  END SUBROUTINE chem_check_data_output_pr
@@ -1976,4 +2055,10 @@
     spec_conc_3 (:,:,:,:) = 0.0_wp
 
+!
+!-- Allocate array to store locally summed-up resolved-scale vertical fluxes.
+    IF ( scalar_advec == 'ws-scheme' )  THEN
+       ALLOCATE( sums_ws_l(nzb:nzt+1,0:threads_per_task-1,nspec) )
+       sums_ws_l = 0.0_wp
+    ENDIF
 
     DO  lsp = 1, nspec
@@ -2070,4 +2155,5 @@
              ENDDO
           ENDIF
+
        ENDIF
 !
@@ -2620,9 +2706,10 @@
 !
 !-- Determine the number of user-defined profiles and append them to the standard data output
-!>> (data_output_pr)
+!-- (data_output_pr)
     max_pr_cs_tmp = 0
     i = 1
 
     DO  WHILE ( data_output_pr(i)  /= ' '  .AND.  i <= 100 )
+
        IF ( TRIM( data_output_pr(i)(1:3) ) == 'kc_' )  THEN
           max_pr_cs_tmp = max_pr_cs_tmp+1
@@ -2708,4 +2795,11 @@
           ENDIF
 
+       CASE ( 'before_timestep' )
+!
+!--       Set array used to sum-up resolved scale fluxes to zero.
+          IF ( ws_scheme_sca )  THEN
+             sums_ws_l = 0.0_wp
+          ENDIF
+
        CASE DEFAULT
           CONTINUE
@@ -2876,4 +2970,5 @@
        IF ( timestep_scheme(1:5) == 'runge' )  THEN
           IF ( ws_scheme_sca )  THEN
+             sums_wschs_ws_l(nzb:,0:) => sums_ws_l(:,:,ilsp)
              CALL advec_s_ws( cs_advc_flags_s, chem_species(ilsp)%conc, 'kc',                      &
                               bc_dirichlet_cs_l  .OR.  bc_radiation_cs_l,                          &
@@ -2986,4 +3081,5 @@
        IF ( timestep_scheme(1:5) == 'runge' )  THEN
           IF ( ws_scheme_sca )  THEN
+             sums_wschs_ws_l(nzb:,0:) => sums_ws_l(nzb:nzt+1,0:threads_per_task-1,ilsp)
              CALL advec_s_ws( cs_advc_flags_s,                                                     &
                               i,                                                                   &
@@ -3029,6 +3125,5 @@
                * ( chem_species(ilsp)%conc(k,j,i) - chem_species(ilsp)%conc_pr_init(k) )           &
                )                                                                                   &
-               * MERGE( 1.0_wp, 0.0_wp,                                                            &
-               BTEST( wall_flags_total_0(k,j,i), 0 )                                               &
+               * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 )                      &
                )
 
@@ -3170,36 +3265,101 @@
     CHARACTER (LEN=*) ::  mode   !<
 
-    INTEGER(iwp) ::  i    !< running index on x-axis
-    INTEGER(iwp) ::  j    !< running index on y-axis
-    INTEGER(iwp) ::  k    !< vertical index counter
-    INTEGER(iwp) ::  lpr  !< running index chem spcs
-    INTEGER(iwp) ::  sr   !< statistical region
-    INTEGER(iwp) ::  tn   !< thread number
+    INTEGER(iwp) ::  i       !< running index on x-axis
+    INTEGER(iwp) ::  j       !< running index on y-axis
+    INTEGER(iwp) ::  k       !< vertical index counter
+    INTEGER(iwp) ::  lpr     !< running index chem spcs
+    INTEGER(iwp) ::  m       !< running index for surface elements
+    INTEGER(iwp) ::  sr      !< statistical region
+    INTEGER(iwp) ::  surf_e  !< end surface index
+    INTEGER(iwp) ::  surf_s  !< start surface index
+    INTEGER(iwp) ::  tn      !< thread number
+
+    REAL(wp)                                            ::  flag     !< topography masking flag
+    REAL(wp), DIMENSION(nzb:nzt+1,0:threads_per_task-1) ::  sums_tmp !< temporary array used to sum-up profiles
 
     IF ( mode == 'profiles' )  THEN
-       !
-!
-!--    Sample on how to calculate horizontally averaged profiles of user-defined quantities. Each
-!--    quantity is identified by the index "pr_palm+#" where "#" is an integer starting from 1.
-!--    These user-profile-numbers must also be assigned to the respective strings given by
-!--    data_output_pr_cs in routine user_check_data_output_pr.
-!--    hom(:,:,:,:) =  dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g.
-!--                     w*pt*), dim-4 = statistical region.
-
-!$OMP DO
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb, nzt+1
-                DO lpr = 1, cs_pr_count
-
-                   sums_l(k,pr_palm+max_pr_user+lpr,tn) = sums_l(k,pr_palm+max_pr_user+lpr,tn) +   &
-                        chem_species(cs_pr_index(lpr))%conc(k,j,i) *                               &
+!
+!--    Sum-up profiles for the species
+       tn = 0
+       !$OMP PARALLEL PRIVATE( i, j, k, tn, lpr, sums_tmp )
+       !$ tn = omp_get_thread_num()
+       !$OMP DO
+       DO  lpr = 1, cs_pr_count_sp
+          sums_tmp(:,tn) = 0.0_wp
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb, nzt+1
+                   sums_tmp(k,tn) = sums_tmp(k,tn) +                                               &
+                        chem_species(cs_pr_index_sp(lpr))%conc(k,j,i) *                            &
                         rmask(j,i,sr)  *                                                           &
-                        MERGE( 1.0_wp, 0.0_wp,                                                     &
-                        BTEST( wall_flags_total_0(k,j,i), 22 ) )
+                        MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 22 ) )
                 ENDDO
              ENDDO
           ENDDO
+          sums_l(nzb:nzt+1,hom_index_spec(lpr),tn) = sums_tmp(nzb:nzt+1,tn)
        ENDDO
+!
+!--    Sum-up profiles for vertical fluxes of the the species. Note, in case of WS5 scheme the
+!--    profiles of resolved-scale fluxes have been already summed-up, while resolved-scale fluxes
+!--    need to be calculated in case of PW scheme. 
+!--    For summation employ a temporary array.
+       !$OMP PARALLEL PRIVATE( i, j, k, tn, lpr, sums_tmp, flag )
+       !$ tn = omp_get_thread_num()
+       !$OMP DO
+       DO  lpr = 1, cs_pr_count_fl_sgs
+          sums_tmp(:,tn) = 0.0_wp
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb, nzt
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 23 ) ) *        &
+                          MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 9  ) )
+                   sums_tmp(k,tn) = sums_tmp(k,tn) -                                               &
+                        0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                                       &
+                               * ( chem_species(cs_pr_index_fl_sgs(lpr))%conc(k+1,j,i) -           &
+                                   chem_species(cs_pr_index_fl_sgs(lpr))%conc(k,j,i) ) *           &
+                        ddzu(k+1) * rmask(j,i,sr)  * flag
+                ENDDO
+!
+!--             Add surface fluxes
+                surf_s = surf_def_h(0)%start_index(j,i)
+                surf_e = surf_def_h(0)%end_index(j,i)
+                DO  m = surf_s, surf_e
+                   k   = surf_def_h(0)%k(m) + surf_def_h(0)%koff
+                   sums_tmp(k,tn) = sums_tmp(k,tn) + surf_def_h(0)%cssws(cs_pr_index_fl_sgs(lpr),m)
+                ENDDO
+                surf_s = surf_def_h(1)%start_index(j,i)
+                surf_e = surf_def_h(1)%end_index(j,i)
+                DO  m = surf_s, surf_e
+                   k   = surf_def_h(1)%k(m)  + surf_def_h(1)%koff
+                   sums_tmp(k,tn) = sums_tmp(k,tn) + surf_def_h(1)%cssws(cs_pr_index_fl_sgs(lpr),m)
+                ENDDO
+                surf_s = surf_lsm_h%start_index(j,i)
+                surf_e = surf_lsm_h%end_index(j,i)
+                DO  m = surf_s, surf_e
+                   k   = surf_lsm_h%k(m) + surf_lsm_h%koff
+                   sums_tmp(k,tn) = sums_tmp(k,tn) + surf_lsm_h%cssws(cs_pr_index_fl_sgs(lpr),m)
+                ENDDO
+                surf_s = surf_usm_h%start_index(j,i)
+                surf_e = surf_usm_h%end_index(j,i)
+                DO  m = surf_s, surf_e
+                   k   = surf_usm_h%k(m) + surf_usm_h%koff
+                   sums_tmp(k,tn) = sums_tmp(k,tn) + surf_usm_h%cssws(cs_pr_index_fl_sgs(lpr),m)
+                ENDDO
+             ENDDO
+          ENDDO
+          sums_l(nzb:nzt+1,hom_index_fl_sgs(lpr),tn) = sums_tmp(nzb:nzt+1,tn)
+       ENDDO
+!
+!--    Resolved-scale fluxes from the WS5 scheme
+       IF ( ws_scheme_sca )  THEN
+          !$OMP PARALLEL PRIVATE( tn, lpr )
+          !$ tn = omp_get_thread_num()
+          !$OMP DO
+          DO  lpr = 1, cs_pr_count_fl_res
+             sums_l(nzb:nzt+1,hom_index_fl_res(lpr),tn) =                                           &
+                                                      sums_ws_l(nzb:nzt+1,tn,cs_pr_index_fl_res(lpr))
+          ENDDO
+       ENDIF
+
     ELSEIF ( mode == 'time_series' )  THEN
 !      @todo
Index: /palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- /palm/trunk/SOURCE/flow_statistics.f90	(revision 4580)
+++ /palm/trunk/SOURCE/flow_statistics.f90	(revision 4581)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Formatting adjustment
+! 
+! 4551 2020-06-02 10:22:25Z suehring
 ! Bugfix in summation for statistical regions
 ! 
@@ -1883,18 +1886,18 @@
        IF ( air_chemistry  .AND.  max_pr_cs > 0 )  THEN
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-             DO  i = 1, max_pr_cs
-                CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+max_pr_user+i,0),       &
-                                    sums(nzb,pr_palm+max_pr_user+i),           &
-                                    nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr )
-             ENDDO
+          DO  i = 1, max_pr_cs
+             CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+max_pr_user+i,0),          &
+                                 sums(nzb,pr_palm+max_pr_user+i),              &
+                                 nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr )
+          ENDDO
        ENDIF
 
        IF ( salsa  .AND.  max_pr_salsa > 0 )  THEN
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-             DO  i = 1, max_pr_salsa
-                CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+max_pr_user+max_pr_cs+i,0),                 &
-                                    sums(nzb,pr_palm+max_pr_user+max_pr_user+i),                   &
-                                    nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr )
-             ENDDO
+          DO  i = 1, max_pr_salsa
+             CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+max_pr_user+max_pr_cs+i,0),                    &
+                                 sums(nzb,pr_palm+max_pr_user+max_pr_user+i),                      &
+                                 nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr )
+          ENDDO
        ENDIF
 
Index: /palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 4580)
+++ /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 4581)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Minor formatting of error message
+! 
+! 4535 2020-05-15 12:07:23Z raasch
 ! bugfix for restart data format query
 ! 
@@ -1038,5 +1041,5 @@
           IF (  .NOT.  land_surface )  THEN
              message_string = 'data_output_pr = ' //                           &
-                              TRIM( data_output_pr(var_count) ) // ' is' //    &
+                              TRIM( data_output_pr(var_count) ) // ' is ' //   &
                               'not implemented for land_surface = .FALSE.'
              CALL message( 'lsm_check_data_output_pr', 'PA0402', 1, 2, 0, 6, 0 )
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 4580)
+++ /palm/trunk/SOURCE/modules.f90	(revision 4581)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! +sums_wschs_ws_l
+!
+! 4574 2020-06-24 16:33:32Z pavelkrc
 !
 ! 4572 2020-06-24 10:09:41Z gronemeier
@@ -1478,4 +1481,6 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_ls_l        !< subdomain sum of large scale forcing and nudging tendencies
 
+    REAL(wp), DIMENSION(:,:), POINTER     ::  sums_wschs_ws_l  !< subdomain sum of vertical chemistry flux w'ch'
+
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  hom_sum             !< sum array for horizontal mean
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  rmask               !< REAL flag array (0.0 or 1.0) for statistic regions
Index: /palm/trunk/SOURCE/nesting_offl_mod.f90
===================================================================
--- /palm/trunk/SOURCE/nesting_offl_mod.f90	(revision 4580)
+++ /palm/trunk/SOURCE/nesting_offl_mod.f90	(revision 4581)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Omit explicit pressure forcing via geostrophic wind components in case of
+! mesoscale nesting.
+! 
+! 4561 2020-06-12 07:05:56Z suehring
 ! Adapt mass-flux correction also for the anelastic approximation
 ! 
@@ -471,10 +475,10 @@
 !
 !--    Read geostrophic wind components
-       DO  t = nest_offl%tind, nest_offl%tind_p
-          CALL get_variable_pr( pids_id, 'ls_forcing_ug', t+1,                 &
-                                nest_offl%ug(t-nest_offl%tind,nzb+1:nzt) )
-          CALL get_variable_pr( pids_id, 'ls_forcing_vg', t+1,                 &
-                                nest_offl%vg(t-nest_offl%tind,nzb+1:nzt) )
-       ENDDO
+!        DO  t = nest_offl%tind, nest_offl%tind_p
+!           CALL get_variable_pr( pids_id, 'ls_forcing_ug', t+1,                 &
+!                                 nest_offl%ug(t-nest_offl%tind,nzb+1:nzt) )
+!           CALL get_variable_pr( pids_id, 'ls_forcing_vg', t+1,                 &
+!                                 nest_offl%vg(t-nest_offl%tind,nzb+1:nzt) )
+!        ENDDO
 !
 !--    Read data at lateral and top boundaries. Please note, at left and
@@ -1567,7 +1571,6 @@
 ! Description:
 !------------------------------------------------------------------------------!
-!>  Update of the geostrophic wind components.
-!>  @todo: update geostrophic wind also in the child domains (should be done
-!>         in the nesting.
+!>  Update of the geostrophic wind components. Note, currently this routine is
+!>  not invoked.
 !------------------------------------------------------------------------------!
     SUBROUTINE nesting_offl_geostrophic_wind
@@ -2135,16 +2138,11 @@
           ENDIF
 !
-!--       Initialize geostrophic wind components. Actually this is already done in
-!--       init_3d_model when initializing_action = 'inifor', however, in speical 
-!--       case of user-defined initialization this will be done here again, in 
-!--       order to have a consistent initialization. 
-          ug(nzb+1:nzt) = nest_offl%ug(0,nzb+1:nzt)
-          vg(nzb+1:nzt) = nest_offl%vg(0,nzb+1:nzt)
-!
-!--       Set bottom and top boundary condition for geostrophic wind components
-          ug(nzt+1) = ug(nzt)
-          vg(nzt+1) = vg(nzt)
-          ug(nzb)   = ug(nzb+1)
-          vg(nzb)   = vg(nzb+1)
+!--       In case of offline nesting the pressure forms itself based on the prescribed lateral
+!--       boundary conditions. Hence, explicit forcing by pressure gradients via geostrophic
+!--       wind components is not necessary and would be canceled out by the perturbation
+!--       pressure otherwise. For this reason, set geostrophic wind components to zero.
+          ug(nzb+1:nzt) = 0.0_wp
+          vg(nzb+1:nzt) = 0.0_wp
+
        ENDIF
 !
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 4580)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 4581)
@@ -20,9 +20,13 @@
 ! Current revisions:
 ! ------------------
-!
-!
+! 
+! 
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Omit explicit pressure forcing via geostrophic wind components in case of
+! mesoscale nesting.
+! 
+! 4578 2020-06-25 15:43:32Z gronemeier
 ! bugfix: removed unused variables
 !
@@ -343,5 +347,4 @@
     USE nesting_offl_mod,                                                                          &
         ONLY:  nesting_offl_bc,                                                                    &
-               nesting_offl_geostrophic_wind,                                                      &
                nesting_offl_input,                                                                 &
                nesting_offl_interpolation_factor,                                                  &
@@ -1006,5 +1009,5 @@
              CALL nesting_offl_interpolation_factor
              CALL nesting_offl_bc
-             CALL nesting_offl_geostrophic_wind
+!              CALL nesting_offl_geostrophic_wind
           ENDIF
 !
Index: /palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4580)
+++ /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4581)
@@ -27,4 +27,7 @@
 ! -----------------
 ! $Id$
+! Missing initialization in case of cyclic_fill runs
+! 
+! 4535 2020-05-15 12:07:23Z raasch
 ! bugfix for restart data format query
 ! 
@@ -4997,6 +5000,5 @@
 !
 !--  Initialization for restart runs
-     IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.                               &
-          TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
+     IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
 
 !
