Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 4510)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 4511)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! call of chem_boundary_conds removed (respective settings are now done in the chemistry module)
+! 
+! 4495 2020-04-13 20:11:20Z raasch
 ! check new restart_data_format parameters
 ! 
@@ -131,7 +134,4 @@
 
     USE chem_modules
-
-    USE chemistry_model_mod,                                                   &
-        ONLY:  chem_boundary_conds
 
     USE control_parameters
@@ -1554,8 +1554,4 @@
 
 !
-!-- Boundary conditions for chemical species
-    IF ( air_chemistry )  CALL chem_boundary_conds( 'init' )
-
-!
 !-- Boundary conditions for horizontal components of wind speed
     IF ( bc_uv_b == 'dirichlet' )  THEN
Index: /palm/trunk/SOURCE/chem_modules.f90
===================================================================
--- /palm/trunk/SOURCE/chem_modules.f90	(revision 4510)
+++ /palm/trunk/SOURCE/chem_modules.f90	(revision 4511)
@@ -22,9 +22,12 @@
 ! Current revisions:
 ! -----------------
-!
+! 
 ! 
 ! Former revisions:
 ! -----------------
 ! $Id$
+! new variables for explicit settings of lateral boundary conditions introduced
+!
+! 4481 2020-03-31 18:55:54Z maronga
 ! added namelist flag 'emiss_read_legacy_mode' to allow concurrent
 ! functioning of new emission read mode under development (ECC)
@@ -92,4 +95,8 @@
 
     CHARACTER (LEN=20) ::  bc_cs_b        = 'dirichlet'         !< namelist parameter: surface boundary condition for concentration
+    CHARACTER (LEN=20) ::  bc_cs_l        = 'undefined'         !< left boundary condition
+    CHARACTER (LEN=20) ::  bc_cs_n        = 'undefined'         !< north boundary condition
+    CHARACTER (LEN=20) ::  bc_cs_r        = 'undefined'         !< right boundary condition
+    CHARACTER (LEN=20) ::  bc_cs_s        = 'undefined'         !< south boundary condition
     CHARACTER (LEN=20) ::  bc_cs_t        = 'initial_gradient'  !< namelist parameter: top boudary condition for concentration
     CHARACTER (LEN=30) ::  chem_mechanism = 'phstatp'           !< namelist parameter: chemistry mechanism
@@ -105,4 +112,9 @@
     CHARACTER (LEN=11), DIMENSION(99) ::  data_output_pr_cs   = 'novalue'  !< namelist parameter: tbc...???
     CHARACTER (LEN=11), DIMENSION(99) ::  surface_csflux_name = 'novalue'  !< namelist parameter: tbc...???
+
+    INTEGER(iwp) ::  communicator_chem      !< stores the number of the MPI communicator to be used
+                                            !< for ghost layer data exchange
+                                            !< 1: cyclic, 2: cyclic along x, 3: cyclic along y,
+                                            !< 4: non-cyclic
 
     INTEGER(iwp) ::  cs_pr_count                           = 0      !< counter for chemical species profiles
@@ -131,4 +143,20 @@
                                                                      !< chemical species near walls and lateral boundaries
 
+    LOGICAL ::  bc_dirichlet_cs_l         = .FALSE.  !< flag for indicating a dirichlet condition at
+                                                     !< the left boundary
+    LOGICAL ::  bc_dirichlet_cs_n         = .FALSE.  !< flag for indicating a dirichlet condition at
+                                                     !< the north boundary
+    LOGICAL ::  bc_dirichlet_cs_r         = .FALSE.  !< flag for indicating a dirichlet condition at
+                                                     !< the right boundary
+    LOGICAL ::  bc_dirichlet_cs_s         = .FALSE.  !< flag for indicating a dirichlet condition at
+                                                     !< the south boundary
+    LOGICAL ::  bc_radiation_cs_l         = .FALSE.  !< flag for indicating a radiation/neumann
+                                                     !< condition at the left boundary
+    LOGICAL ::  bc_radiation_cs_n         = .FALSE.  !< flag for indicating a radiation/neumann
+                                                     !< condition at the north boundary
+    LOGICAL ::  bc_radiation_cs_r         = .FALSE.  !< flag for indicating a radiation/neumann
+                                                     !< condition at the right boundary
+    LOGICAL ::  bc_radiation_cs_s         = .FALSE.  !< flag for indicating a radiation/neumann
+                                                     !< condition at the south boundary
     LOGICAL ::  constant_top_csflux(99)   = .TRUE.   !< internal flag, set to .FALSE. if no top_csflux is prescribed
     LOGICAL ::  constant_csflux(99)       = .TRUE.   !< internal flag, set to .FALSE. if no surface_csflux is prescribed
Index: /palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 4510)
+++ /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 4511)
@@ -27,6 +27,9 @@
 ! -----------------
 ! $Id$
+! decycling replaced by explicit setting of lateral boundary conditions
+!
+! 4487 2020-04-03 09:38:20Z raasch
 ! bugfix for subroutine calls that contain the decycle_chem switches as arguments
-! 
+!
 ! 4481 2020-03-31 18:55:54Z maronga
 ! use statement for exchange horiz added,
@@ -347,20 +350,4 @@
     REAL(kind=wp), DIMENSION(nkppctrl)                     ::  rcntrl       !< 20 real parameters for fine tuning of KPP code
                                                                             !< (e.g starting internal timestep of solver)
-!
-!-- Decycling of chemistry variables: Dirichlet BCs with cyclic is frequently not
-!-- approproate for chemicals compounds since they may accumulate too much. 
-!-- If no proper boundary conditions from a DYNAMIC input file are available,
-!-- de-cycling applies the initial profiles at the inflow boundaries for 
-!-- Dirichlet boundary conditions
-    LOGICAL ::  decycle_chem_lr    = .FALSE.    !< switch for setting decycling in left-right direction
-    LOGICAL ::  decycle_chem_ns    = .FALSE.    !< switch for setting decycling in south-norht direction
-    CHARACTER (LEN=20), DIMENSION(4) ::  decycle_method = &
-         (/'dirichlet','dirichlet','dirichlet','dirichlet'/)
-                              !< Decycling method at horizontal boundaries, 
-                              !< 1=left, 2=right, 3=south, 4=north
-                              !< dirichlet = initial size distribution and 
-                              !< chemical composition set for the ghost and 
-                              !< first three layers
-                              !< neumann = zero gradient
 
     REAL(kind=wp), PUBLIC ::  cs_time_step = 0._wp
@@ -438,9 +425,4 @@
     END INTERFACE chem_3d_data_averaging
 
-    INTERFACE chem_boundary_conds
-       MODULE PROCEDURE chem_boundary_conds
-       MODULE PROCEDURE chem_boundary_conds_decycle
-    END INTERFACE chem_boundary_conds
-
     INTERFACE chem_boundary_conditions
        MODULE PROCEDURE chem_boundary_conditions
@@ -590,8 +572,4 @@
     END INTERFACE rc_rctot
 
-!    INTERFACE rc_comp_point_rc_eff
-!       MODULE PROCEDURE rc_comp_point_rc_eff 
-!    END INTERFACE rc_comp_point_rc_eff
-
     INTERFACE drydepo_aero_zhang_vd
        MODULE PROCEDURE drydepo_aero_zhang_vd 
@@ -604,13 +582,11 @@
 
 
-    PUBLIC chem_3d_data_averaging, chem_boundary_conds, chem_boundary_conditions, &
-            chem_boundary_conds_decycle, chem_check_data_output,              &
-         chem_check_data_output_pr, chem_check_parameters,                    &
-         chem_data_output_2d, chem_data_output_3d, chem_data_output_mask,     &
-         chem_define_netcdf_grid, chem_header, chem_init, chem_init_arrays,   &
-         chem_init_profiles, chem_integrate, chem_parin,                      &
-         chem_actions, chem_prognostic_equations, chem_rrd_local,             &
-         chem_statistics, chem_swap_timelevel, chem_wrd_local, chem_depo,     &
-         chem_non_advective_processes, chem_exchange_horiz_bounds
+    PUBLIC chem_3d_data_averaging, chem_boundary_conditions, chem_check_data_output,               &
+           chem_check_data_output_pr, chem_check_parameters, chem_data_output_2d,                  &
+           chem_data_output_3d, chem_data_output_mask, chem_define_netcdf_grid, chem_header,       &
+           chem_init, chem_init_arrays, chem_init_profiles, chem_integrate, chem_parin,            &
+           chem_actions, chem_prognostic_equations, chem_rrd_local, chem_statistics,               &
+           chem_swap_timelevel, chem_wrd_local, chem_depo, chem_non_advective_processes,           &
+           chem_exchange_horiz_bounds
 
  CONTAINS
@@ -742,16 +718,12 @@
 !> Subroutine to initialize and set all boundary conditions for chemical species
 !------------------------------------------------------------------------------!
- SUBROUTINE chem_boundary_conds( mode )                                           
-
-    USE control_parameters,                                                    &  
-        ONLY:  bc_radiation_l, bc_radiation_n, bc_radiation_r, bc_radiation_s
-
-    USE arrays_3d,                                                             &     
-        ONLY:  dzu                                                
+ SUBROUTINE chem_boundary_conditions( horizontal_conditions_only )
+
+    USE arrays_3d,                                                             &
+        ONLY:  dzu
 
     USE surface_mod,                                                           &
         ONLY:  bc_h                                                           
 
-    CHARACTER (LEN=*), INTENT(IN) ::  mode
     INTEGER(iwp) ::  i                            !< grid index x direction.
     INTEGER(iwp) ::  j                            !< grid index y direction.
@@ -761,309 +733,245 @@
     INTEGER(iwp) ::  lsp                          !< running index for chem spcs.
 
-
-    SELECT CASE ( TRIM( mode ) )        
-       CASE ( 'init' )
-
-          IF ( bc_cs_b == 'dirichlet' )  THEN
-             ibc_cs_b = 0 
-          ELSEIF ( bc_cs_b == 'neumann' )  THEN
-             ibc_cs_b = 1 
-          ELSE
-             message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"'  
-             CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 )
-          ENDIF                                                                  
-!
-!--       Set Integer flags and check for possible erroneous settings for top
-!--       boundary condition. 
-          IF ( bc_cs_t == 'dirichlet' )  THEN
-             ibc_cs_t = 0 
-          ELSEIF ( bc_cs_t == 'neumann' )  THEN
-             ibc_cs_t = 1
-          ELSEIF ( bc_cs_t == 'initial_gradient' )  THEN
-             ibc_cs_t = 2
-          ELSEIF ( bc_cs_t == 'nested' )  THEN          
-             ibc_cs_t = 3
-          ELSE
-             message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"'      
-             CALL message( 'check_parameters', 'CM0430', 1, 2, 0, 6, 0 )
-          ENDIF
-
-       CASE ( 'set_bc_bottomtop' )                   
-!
-!--       Boundary condtions for chemical species at horizontal walls      
-          DO  lsp = 1, nspec                                                      
-             IF ( ibc_cs_b == 0 )  THEN
-                DO  l = 0, 1 
-                   !$OMP PARALLEL DO PRIVATE( i, j, k )
-                   DO  m = 1, bc_h(l)%ns
-                       i = bc_h(l)%i(m)            
-                       j = bc_h(l)%j(m)
-                       k = bc_h(l)%k(m)
-                      chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =           &
+    LOGICAL, OPTIONAL ::  horizontal_conditions_only  !< switch to set horizontal bc only
+
+
+    IF ( .NOT. PRESENT( horizontal_conditions_only ) )  THEN
+!
+!--    Boundary condtions for chemical species at horizontal walls
+       DO  lsp = 1, nspec
+
+          IF ( ibc_cs_b == 0 )  THEN
+             DO  l = 0, 1
+                !$OMP PARALLEL DO PRIVATE( i, j, k )
+                DO  m = 1, bc_h(l)%ns
+                    i = bc_h(l)%i(m)
+                    j = bc_h(l)%j(m)
+                    k = bc_h(l)%k(m)
+                   chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =           &
                                       chem_species(lsp)%conc(k+bc_h(l)%koff,j,i) 
-                   ENDDO                                        
-                ENDDO                                       
-
-             ELSEIF ( ibc_cs_b == 1 )  THEN
-!
-!--             in boundary_conds there is som extra loop over m here for passive tracer
-                DO  l = 0, 1
-                   !$OMP PARALLEL DO PRIVATE( i, j, k )                                           
-                   DO m = 1, bc_h(l)%ns
-                      i = bc_h(l)%i(m)            
-                      j = bc_h(l)%j(m)
-                      k = bc_h(l)%k(m)
-                      chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =           &
+                ENDDO
+             ENDDO
+
+          ELSEIF ( ibc_cs_b == 1 )  THEN
+!
+!--          In boundary_conds there is som extra loop over m here for passive tracer
+!>           TODO: clarify the meaning of the above comment. Explain in more detail or remove it. (Siggi)
+             DO  l = 0, 1
+                !$OMP PARALLEL DO PRIVATE( i, j, k )
+                DO m = 1, bc_h(l)%ns
+                   i = bc_h(l)%i(m)
+                   j = bc_h(l)%j(m)
+                   k = bc_h(l)%k(m)
+                   chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =           &
                                          chem_species(lsp)%conc_p(k,j,i)
-
-                   ENDDO
                 ENDDO
-             ENDIF
-       ENDDO    ! end lsp loop  
-!
-!--    Top boundary conditions for chemical species - Should this not be done for all species?
-          IF ( ibc_cs_t == 0 )  THEN                    
-             DO  lsp = 1, nspec
-                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)        
-             ENDDO
-          ELSEIF ( ibc_cs_t == 1 )  THEN
-             DO  lsp = 1, nspec
-                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:)
-             ENDDO
-          ELSEIF ( ibc_cs_t == 2 )  THEN
-             DO  lsp = 1, nspec
-                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1)
              ENDDO
           ENDIF
 
-       CASE ( 'set_bc_lateral' )                       
-!
-!--             Lateral boundary conditions for chem species at inflow boundary
-!--             are automatically set when chem_species concentration is
-!--             initialized. The initially set value at the inflow boundary is not
-!--             touched during time integration, hence, this boundary value remains
-!--             at a constant value, which is the concentration that flows into the
-!--             domain.                                                           
-!--             Lateral boundary conditions for chem species at outflow boundary 
-
-          IF ( bc_radiation_s )  THEN
-             DO  lsp = 1, nspec
-                chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:)
-             ENDDO
-          ELSEIF ( bc_radiation_n )  THEN
-             DO  lsp = 1, nspec
-                chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:)
-             ENDDO
-          ELSEIF ( bc_radiation_l )  THEN
-             DO  lsp = 1, nspec
-                chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl)
-             ENDDO
-          ELSEIF ( bc_radiation_r )  THEN
-             DO  lsp = 1, nspec
-                chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr)
-             ENDDO
-          ENDIF
-
-    END SELECT
-
- END SUBROUTINE chem_boundary_conds
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Subroutine for boundary conditions
-!------------------------------------------------------------------------------!
- SUBROUTINE chem_boundary_conditions
-
-    IMPLICIT NONE
-
-    INTEGER(iwp) ::  lsp             !<
-    INTEGER(iwp) ::  lsp_usr         !<
-
-!
-!-- Top/bottom boundary conditions for chemical species
-    CALL chem_boundary_conds( 'set_bc_bottomtop' )
-!
-!-- Lateral boundary conditions for chemical species
-    CALL chem_boundary_conds( 'set_bc_lateral' )
-
-!
-!--  Boundary conditions for prognostic quantitites of other modules:
-!--  Here, only decycling is carried out
-     DO  lsp = 1, nvar
-        lsp_usr = 1
-        DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
-           IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN
-              CALL chem_boundary_conds( chem_species(lsp)%conc_p,                          &
-                                        chem_species(lsp)%conc_pr_init )
-           ENDIF
-           lsp_usr = lsp_usr + 1
-        ENDDO
-     ENDDO
-
-
- END SUBROUTINE chem_boundary_conditions
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Boundary conditions for prognostic variables in chemistry: decycling in the 
-!> x-direction-
-!> Decycling of chemistry variables: Dirichlet BCs with cyclic is frequently not
-!> approproate for chemicals compounds since they may accumulate too much.
-!> If no proper boundary conditions from a DYNAMIC input file are available,
-!> de-cycling applies the initial profiles at the inflow boundaries for
-!> Dirichlet boundary conditions
-!------------------------------------------------------------------------------!
- SUBROUTINE chem_boundary_conds_decycle( cs_3d, cs_pr_init )
-
-    USE control_parameters,                                                                         &
-        ONLY:  nesting_offline
-
-    INTEGER(iwp) ::  boundary  !< 
-    INTEGER(iwp) ::  ee        !<
-    INTEGER(iwp) ::  copied    !< 
-    INTEGER(iwp) ::  i         !< 
-    INTEGER(iwp) ::  j         !< 
-    INTEGER(iwp) ::  k         !< 
-    INTEGER(iwp) ::  ss        !<
-
-    REAL(wp), DIMENSION(nzb:nzt+1) ::  cs_pr_init
-    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  cs_3d
-    REAL(wp) ::  flag          !< flag to mask topography grid points
-
-
-    flag = 0.0_wp
-!
-!-- Skip input if forcing from a larger-scale model is applied
-    IF ( nesting_offline  .AND.  nesting_offline_chem )  RETURN
-!
-!-- Left and right boundaries
-    IF ( decycle_chem_lr  .AND.  bc_lr_cyc )  THEN
-
-       DO  boundary = 1, 2
-
-          IF ( decycle_method(boundary) == 'dirichlet' )  THEN
-!
-!--          Initial profile is copied to ghost and first three layers
-             ss = 1
-             ee = 0
-             IF ( boundary == 1  .AND.  nxl == 0 )  THEN
-                ss = nxlg
-                ee = nxl-1
-             ELSEIF ( boundary == 2  .AND.  nxr == nx )  THEN
-                ss = nxr+1
-                ee = nxrg
-             ENDIF
-
-             DO  i = ss, ee
-                DO  j = nysg, nyng
-                   DO  k = nzb+1, nzt
-                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
-                      cs_3d(k,j,i) = cs_pr_init(k) * flag
-                   ENDDO
-                ENDDO
-             ENDDO
-
-          ELSEIF ( decycle_method(boundary) == 'neumann' )  THEN
-!
-!--          The value at the boundary is copied to the ghost layers to simulate
-!--          an outlet with zero gradient
-             ss = 1
-             ee = 0
-             IF ( boundary == 1  .AND.  nxl == 0 )  THEN
-                ss = nxlg
-                ee = nxl-1
-                copied = nxl
-             ELSEIF ( boundary == 2  .AND.  nxr == nx )  THEN
-                ss = nxr+1
-                ee = nxrg
-                copied = nxr
-             ENDIF
-
-             DO  i = ss, ee
-                DO  j = nysg, nyng
-                   DO  k = nzb+1, nzt
-                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
-                      cs_3d(k,j,i) = cs_3d(k,j,copied) * flag
-                   ENDDO
-                ENDDO
-             ENDDO
-
-          ELSE
-             WRITE(message_string,*)                                           &
-                  'unknown decycling method: decycle_method (', &
-                  boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
-             CALL message( 'chem_boundary_conds_decycle', 'CM0431',           &
-                  1, 2, 0, 6, 0 )
+       ENDDO    ! end lsp loop  
+
+!
+!--    Top boundary conditions for chemical species - Should this not be done for all species?
+!>     TODO: This question also needs to be clarified. I guess it can be removed because the loop
+!>           already runs over all species? (Siggi)
+       DO  lsp = 1, nspec
+          IF ( ibc_cs_t == 0 )  THEN                    
+             chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)
+          ELSEIF ( ibc_cs_t == 1 )  THEN
+             chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:)
+          ELSEIF ( ibc_cs_t == 2 )  THEN
+             chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) +             &
+                                                   bc_cs_t_val(lsp) * dzu(nzt+1)
           ENDIF
        ENDDO
+
+!
+!--    Lateral boundary conditions.
+!--    Dirichlet conditions have been already set when chem_species concentration is initialized.
+!--    The initially set value is not touched during time integration, hence, this boundary value
+!--    remains at a constant value.
+!--    Neumann conditions:
+       IF ( bc_radiation_cs_s )  THEN
+          DO  lsp = 1, nspec
+             chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:)
+          ENDDO
+       ENDIF
+       IF ( bc_radiation_cs_n )  THEN
+          DO  lsp = 1, nspec
+             chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:)
+          ENDDO
+       ENDIF
+       IF ( bc_radiation_cs_l )  THEN
+          DO  lsp = 1, nspec
+             chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl)
+          ENDDO
+       ENDIF
+       IF ( bc_radiation_cs_r )  THEN
+          DO  lsp = 1, nspec
+             chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr)
+          ENDDO
+       ENDIF
+
+!--    For testing: set Dirichlet conditions for all boundaries
+!      IF ( bc_dirichlet_cs_s )  THEN
+!         DO  lsp = 1, nspec
+!            DO  i = nxlg, nxrg
+!               DO  j = nysg, nys-1
+!                  DO  k = nzb, nzt
+!                     IF ( k /= nzb )  THEN
+!                        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!                     ELSE
+!                        flag = 1.0
+!                     ENDIF
+!                     chem_species(lsp)%conc_p(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag
+!                  ENDDO
+!               ENDDO
+!            ENDDO
+!         ENDDO
+!      ENDIF
+!      IF ( bc_dirichlet_cs_n )  THEN
+!         DO  lsp = 1, nspec
+!            DO  i = nxlg, nxrg
+!               DO  j = nyn+1, nyng
+!                  DO  k = nzb, nzt
+!                     IF ( k /= nzb )  THEN
+!                        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!                     ELSE
+!                        flag = 1.0
+!                     ENDIF
+!                     chem_species(lsp)%conc_p(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag
+!                  ENDDO
+!               ENDDO
+!            ENDDO
+!         ENDDO
+!      ENDIF
+!      IF ( bc_dirichlet_cs_l )  THEN
+!         DO  lsp = 1, nspec
+!            DO  i = nxlg, nxl-1
+!               DO  j = nysg, nyng
+!                  DO  k = nzb, nzt
+!                     IF ( k /= nzb )  THEN
+!                        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!                     ELSE
+!                        flag = 1.0
+!                     ENDIF
+!                     chem_species(lsp)%conc_p(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag
+!                  ENDDO
+!               ENDDO
+!            ENDDO
+!         ENDDO
+!      ENDIF
+!      IF ( bc_dirichlet_cs_r )  THEN
+!         DO  lsp = 1, nspec
+!            DO  i = nxr+1, nxrg
+!               DO  j = nysg, nyng
+!                  DO  k = nzb, nzt
+!                     IF ( k /= nzb )  THEN
+!                        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!                     ELSE
+!                        flag = 1.0
+!                     ENDIF
+!                     chem_species(lsp)%conc_p(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag
+!                  ENDDO
+!               ENDDO
+!            ENDDO
+!         ENDDO
+!      ENDIF
+
+    ELSE
+!
+!--    Lateral Neumann booundary conditions for timelevel t.
+!--    This branch is executed when routine is called after the non-advective processes /
+!--    before the prognostic equations.
+       IF ( bc_radiation_cs_s )  THEN
+          DO  lsp = 1, nspec
+             chem_species(lsp)%conc(:,nys-1,:) = chem_species(lsp)%conc(:,nys,:)
+          ENDDO
+       ENDIF
+       IF ( bc_radiation_cs_n )  THEN
+          DO  lsp = 1, nspec
+             chem_species(lsp)%conc(:,nyn+1,:) = chem_species(lsp)%conc(:,nyn,:)
+          ENDDO
+       ENDIF
+       IF ( bc_radiation_cs_l )  THEN
+          DO  lsp = 1, nspec
+             chem_species(lsp)%conc(:,:,nxl-1) = chem_species(lsp)%conc(:,:,nxl)
+          ENDDO
+       ENDIF
+       IF ( bc_radiation_cs_r )  THEN
+          DO  lsp = 1, nspec
+             chem_species(lsp)%conc(:,:,nxr+1) = chem_species(lsp)%conc(:,:,nxr)
+          ENDDO
+       ENDIF
+
+!--    For testing: set Dirichlet conditions for all boundaries
+!      IF ( bc_dirichlet_cs_s )  THEN
+!         DO  lsp = 1, nspec
+!            DO  i = nxlg, nxrg
+!               DO  j = nysg, nys-1
+!                  DO  k = nzb, nzt
+!                     IF ( k /= nzb )  THEN
+!                        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!                     ELSE
+!                        flag = 1.0
+!                     ENDIF
+!                     chem_species(lsp)%conc(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag
+!                  ENDDO
+!               ENDDO
+!            ENDDO
+!         ENDDO
+!      ENDIF
+!      IF ( bc_dirichlet_cs_n )  THEN
+!         DO  lsp = 1, nspec
+!            DO  i = nxlg, nxrg
+!               DO  j = nyn+1, nyng
+!                  DO  k = nzb, nzt
+!                     IF ( k /= nzb )  THEN
+!                        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!                     ELSE
+!                        flag = 1.0
+!                     ENDIF
+!                     chem_species(lsp)%conc(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag
+!                  ENDDO
+!               ENDDO
+!            ENDDO
+!         ENDDO
+!      ENDIF
+!      IF ( bc_dirichlet_cs_l )  THEN
+!         DO  lsp = 1, nspec
+!            DO  i = nxlg, nxl-1
+!               DO  j = nysg, nyng
+!                  DO  k = nzb, nzt
+!                     IF ( k /= nzb )  THEN
+!                        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!                     ELSE
+!                        flag = 1.0
+!                     ENDIF
+!                     chem_species(lsp)%conc(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag
+!                  ENDDO
+!               ENDDO
+!            ENDDO
+!         ENDDO
+!      ENDIF
+!      IF ( bc_dirichlet_cs_r )  THEN
+!         DO  lsp = 1, nspec
+!            DO  i = nxr+1, nxrg
+!               DO  j = nysg, nyng
+!                  DO  k = nzb, nzt
+!                     IF ( k /= nzb )  THEN
+!                        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!                     ELSE
+!                        flag = 1.0
+!                     ENDIF
+!                     chem_species(lsp)%conc(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag
+!                  ENDDO
+!               ENDDO
+!            ENDDO
+!         ENDDO
+!      ENDIF
+
     ENDIF
-!
-!-- South and north boundaries
-    IF ( decycle_chem_ns  .AND.  bc_ns_cyc )  THEN
-
-       DO  boundary = 3, 4
-
-          IF ( decycle_method(boundary) == 'dirichlet' )  THEN
-!
-!--          Initial profile is copied to ghost and first three layers
-             ss = 1
-             ee = 0
-             IF ( boundary == 3  .AND.  nys == 0 )  THEN
-                ss = nysg
-                ee = nys-1
-             ELSEIF ( boundary == 4  .AND.  nyn == ny )  THEN
-                ss = nyn+1
-                ee = nyng
-             ENDIF
-
-             DO  i = nxlg, nxrg
-                DO  j = ss, ee
-                   DO  k = nzb+1, nzt
-                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
-                      cs_3d(k,j,i) = cs_pr_init(k) * flag
-                   ENDDO
-                ENDDO
-             ENDDO
-
-
-          ELSEIF ( decycle_method(boundary) == 'neumann' )  THEN
-!
-!--          The value at the boundary is copied to the ghost layers to simulate
-!--          an outlet with zero gradient
-             ss = 1
-             ee = 0
-             IF ( boundary == 3  .AND.  nys == 0 )  THEN
-                ss = nysg
-                ee = nys-1
-                copied = nys
-             ELSEIF ( boundary == 4  .AND.  nyn == ny )  THEN
-                ss = nyn+1
-                ee = nyng
-                copied = nyn
-             ENDIF
-
-             DO  i = nxlg, nxrg
-                DO  j = ss, ee
-                   DO  k = nzb+1, nzt
-                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
-                      cs_3d(k,j,i) = cs_3d(k,copied,i) * flag
-                   ENDDO
-                ENDDO
-             ENDDO
-
-          ELSE
-             WRITE(message_string,*)                                           &
-                  'unknown decycling method: decycle_method (', &
-                  boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
-             CALL message( 'chem_boundary_conds_decycle', 'CM0432',           &
-                  1, 2, 0, 6, 0 )
-          ENDIF
-       ENDDO
-    ENDIF
-
-
- END SUBROUTINE chem_boundary_conds_decycle
+
+ END SUBROUTINE chem_boundary_conditions
 
 
@@ -1194,4 +1102,6 @@
  SUBROUTINE chem_check_parameters
 
+    USE control_parameters,                                                                        &
+        ONLY:  bc_lr, bc_ns
 
     LOGICAL  ::  found
@@ -1223,11 +1133,4 @@
     ENDIF
 !
-!-- check for  decycling of chem species
-    IF ( (.NOT. any(decycle_method == 'neumann') ) .AND. (.NOT. any(decycle_method == 'dirichlet') ) )  THEN
-       message_string = 'Incorrect boundary conditions. Only neumann or '   &
-                // 'dirichlet &available for decycling chemical species '
-       CALL message( 'chem_check_parameters', 'CM0426', 1, 2, 0, 6, 0 )
-    ENDIF
-!
 !-- check for chemical mechanism used 
     CALL get_mechanism_name
@@ -1236,4 +1139,30 @@
        CALL message( 'chem_check_parameters', 'CM0462', 1, 2, 0, 6, 0 )
     ENDIF
+
+!
+!-- Check bottom boundary condition and set internal steering parameter
+    IF ( bc_cs_b == 'dirichlet' )  THEN
+       ibc_cs_b = 0
+    ELSEIF ( bc_cs_b == 'neumann' )  THEN
+       ibc_cs_b = 1
+    ELSE
+       message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"'
+       CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 )
+    ENDIF
+!
+!-- Check top boundary condition and set internal steering parameter
+    IF ( bc_cs_t == 'dirichlet' )  THEN
+       ibc_cs_t = 0
+    ELSEIF ( bc_cs_t == 'neumann' )  THEN
+       ibc_cs_t = 1
+    ELSEIF ( bc_cs_t == 'initial_gradient' )  THEN
+       ibc_cs_t = 2
+    ELSEIF ( bc_cs_t == 'nested' )  THEN
+       ibc_cs_t = 3
+    ELSE
+       message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"'
+       CALL message( 'check_parameters', 'CM0430', 1, 2, 0, 6, 0 )
+    ENDIF
+
 !
 !-- If nesting_chem = .F., set top boundary condition to its default value
@@ -1242,9 +1171,107 @@
        bc_cs_t = 'initial_gradient'
     ENDIF
+
+!
+!-- Check left and right boundary conditions. First set default value if not set by user.
+    IF ( bc_cs_l == 'undefined' )  THEN
+       IF ( bc_lr == 'cyclic' )  THEN
+          bc_cs_l = 'cyclic'
+       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
+          bc_cs_l = 'dirichlet'
+       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
+          bc_cs_l = 'neumann'
+       ENDIF
+    ENDIF
+    IF ( bc_cs_r == 'undefined' )  THEN
+       IF ( bc_lr == 'cyclic' )  THEN
+          bc_cs_r = 'cyclic'
+       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
+          bc_cs_r = 'neumann'
+       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
+          bc_cs_r = 'dirichlet'
+       ENDIF
+    ENDIF
+    IF ( bc_cs_l /= 'dirichlet'  .AND.  bc_cs_l /= 'neumann'  .AND.  bc_cs_l /= 'cyclic' )  THEN
+       message_string = 'unknown boundary condition: bc_cs_l = "' // TRIM( bc_cs_l ) // '"'
+       CALL message( 'chem_check_parameters','PA0505', 1, 2, 0, 6, 0 )
+    ENDIF
+    IF ( bc_cs_r /= 'dirichlet'  .AND.  bc_cs_r /= 'neumann'  .AND.  bc_cs_r /= 'cyclic' )  THEN
+       message_string = 'unknown boundary condition: bc_cs_r = "' // TRIM( bc_cs_r ) // '"'
+       CALL message( 'chem_check_parameters','PA0551', 1, 2, 0, 6, 0 )
+    ENDIF
+!
+!-- Check north and south boundary conditions. First set default value if not set by user.
+    IF ( bc_cs_n == 'undefined' )  THEN
+       IF ( bc_ns == 'cyclic' )  THEN
+          bc_cs_n = 'cyclic'
+       ELSEIF ( bc_ns == 'dirichlet/radiation' )  THEN
+          bc_cs_n = 'dirichlet'
+       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
+          bc_cs_n = 'neumann'
+       ENDIF
+    ENDIF
+    IF ( bc_cs_s == 'undefined' )  THEN
+       IF ( bc_ns == 'cyclic' )  THEN
+          bc_cs_s = 'cyclic'
+       ELSEIF ( bc_ns == 'dirichlet/radiation' )  THEN
+          bc_cs_s = 'neumann'
+       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
+          bc_cs_s = 'dirichlet'
+       ENDIF
+    ENDIF
+    IF ( bc_cs_n /= 'dirichlet'  .AND.  bc_cs_n /= 'neumann'  .AND.  bc_cs_n /= 'cyclic' )  THEN
+       message_string = 'unknown boundary condition: bc_cs_n = "' // TRIM( bc_cs_n ) // '"'
+       CALL message( 'chem_check_parameters','PA0714', 1, 2, 0, 6, 0 )
+    ENDIF
+    IF ( bc_cs_s /= 'dirichlet'  .AND.  bc_cs_s /= 'neumann'  .AND.  bc_cs_s /= 'cyclic' )  THEN
+       message_string = 'unknown boundary condition: bc_cs_s = "' // TRIM( bc_cs_s ) // '"'
+       CALL message( 'chem_check_parameters','PA0715', 1, 2, 0, 6, 0 )
+    ENDIF
+!
+!-- Cyclic conditions must be set identically at opposing boundaries
+    IF ( ( bc_cs_l == 'cyclic' .AND. bc_cs_r /= 'cyclic' )  .OR.                                 &
+         ( bc_cs_r == 'cyclic' .AND. bc_cs_l /= 'cyclic' ) )  THEN
+       message_string = 'boundary conditions bc_cs_l and bc_cs_r must both be cyclic or non-cyclic'
+       CALL message( 'chem_check_parameters','PA0716', 1, 2, 0, 6, 0 )
+    ENDIF
+    IF ( ( bc_cs_n == 'cyclic' .AND. bc_cs_s /= 'cyclic' )  .OR.                                 &
+         ( bc_cs_s == 'cyclic' .AND. bc_cs_n /= 'cyclic' ) )  THEN
+       message_string = 'boundary conditions bc_cs_n and bc_cs_s must both be cyclic or non-cyclic'
+       CALL message( 'chem_check_parameters','PA0717', 1, 2, 0, 6, 0 )
+    ENDIF
+!
+!-- Set the switches that control application of horizontal boundary conditions at the boundaries
+!-- of the total domain
+    IF ( bc_cs_n == 'dirichlet'  .AND.  nyn == ny )  bc_dirichlet_cs_n = .TRUE.
+    IF ( bc_cs_n == 'neumann'    .AND.  nyn == ny )  bc_radiation_cs_n = .TRUE.
+    IF ( bc_cs_s == 'dirichlet'  .AND.  nys ==  0 )  bc_dirichlet_cs_s = .TRUE.
+    IF ( bc_cs_s == 'neumann'    .AND.  nys ==  0 )  bc_radiation_cs_s = .TRUE.
+    IF ( bc_cs_l == 'dirichlet'  .AND.  nxl ==  0 )  bc_dirichlet_cs_l = .TRUE.
+    IF ( bc_cs_l == 'neumann'    .AND.  nxl ==  0 )  bc_radiation_cs_l = .TRUE.
+    IF ( bc_cs_r == 'dirichlet'  .AND.  nxr == nx )  bc_dirichlet_cs_r = .TRUE.
+    IF ( bc_cs_r == 'neumann'    .AND.  nxr == nx )  bc_radiation_cs_r = .TRUE.
+
+!
+!-- Set the communicator to be used for ghost layer data exchange
+!-- 1: cyclic, 2: cyclic along x, 3: cyclic along y, 4: non-cyclic
+    IF ( bc_cs_l == 'cyclic' )  THEN
+       IF ( bc_cs_s == 'cyclic' )  THEN
+          communicator_chem = 1
+       ELSE
+          communicator_chem = 2
+       ENDIF
+    ELSE
+       IF ( bc_cs_s == 'cyclic' )  THEN
+          communicator_chem = 3
+       ELSE
+          communicator_chem = 4
+       ENDIF
+    ENDIF
+
 !
 !-- chem_check_parameters is called before the array chem_species is allocated!
 !-- temporary switch of this part of the check
 !    RETURN                !bK commented
-
+!> TODO: this workaround definitely needs to be removed from here!!!
     CALL chem_init_internal
 !
@@ -1732,5 +1759,4 @@
 !
 !-- initializatoin of Surface and profile chemical species
-
     lsp = 1
     docsinit_chr ='Chemical species for initial surface and profile emissions: ' 
@@ -1750,8 +1776,11 @@
     IF ( nesting_chem )  WRITE( io, 12 )  nesting_chem
     IF ( nesting_offline_chem )  WRITE( io, 13 )  nesting_offline_chem
+
+    WRITE( io, 14 )  TRIM( bc_cs_b ), TRIM( bc_cs_t ), TRIM( bc_cs_s ), TRIM( bc_cs_n ),           &
+                     TRIM( bc_cs_l ), TRIM( bc_cs_r )
+
 !
 !-- number of variable and fix chemical species and number of reactions
     cs_fixed = nspec - nvar
-
     WRITE ( io, * ) '   --> Chemical Mechanism        : ', cs_mech 
     WRITE ( io, * ) '   --> Chemical species, variable: ', nvar
@@ -1774,6 +1803,9 @@
 12  FORMAT (/'   Nesting for chemistry variables: ', L1 )
 13  FORMAT (/'   Offline nesting for chemistry variables: ', L1 )
-!
-!
+14  FORMAT (/'   Boundary conditions for chemical species:', /                                     &
+             '      bottom/top:   ',A10,' / ',A10, /                                               &
+             '      north/south:  ',A10,' / ',A10, /                                               &
+             '      left/right:   ',A10,' / ',A10)
+
  END SUBROUTINE chem_header
 
@@ -1884,4 +1916,5 @@
     INTEGER(iwp) ::  lpr_lev           !< running index for chem spcs profile level
 
+    REAL(wp)     ::  flag              !< flag for masking topography/building grid points
 !
 !-- 20200203 ECC
@@ -1948,13 +1981,5 @@
 
     ENDDO
-!
-!-- Set control flags for decycling only at lateral boundary cores, within the 
-!-- inner cores the decycle flag is set to .False.. Even though it does not 
-!-- affect the setting of chemistry boundary conditions, this flag is used to
-!-- set advection control flags appropriately. 
-    decycle_chem_lr = MERGE( decycle_chem_lr, .FALSE.,                         &
-                             nxl == 0  .OR.  nxr == nx )
-    decycle_chem_ns = MERGE( decycle_chem_ns, .FALSE.,                         &
-                             nys == 0  .OR.  nyn == ny )
+
 !
 !-- For some passive scalars decycling may be enabled. This case, the lateral 
@@ -1982,10 +2007,9 @@
 !--    Initialize with flag 31 only. 
        cs_advc_flags_s = 0
-       cs_advc_flags_s = MERGE( IBSET( cs_advc_flags_s, 31 ), 0,               &
-                                BTEST( wall_flags_total_0, 31 ) )
-
-       IF ( decycle_chem_ns )  THEN
+       cs_advc_flags_s = MERGE( IBSET( cs_advc_flags_s, 31 ), 0, BTEST( wall_flags_total_0, 31 ) )
+
+       IF ( bc_dirichlet_cs_n .OR. bc_dirichlet_cs_s )  THEN
           IF ( nys == 0  )  THEN
-             DO  i = 1, nbgp     
+             DO  i = 1, nbgp
                 cs_advc_flags_s(:,nys-i,:) = MERGE(                            &
                                         IBSET( cs_advc_flags_s(:,nys,:), 31 ), &
@@ -1996,5 +2020,5 @@
           ENDIF
           IF ( nyn == ny )  THEN
-             DO  i = 1, nbgp  
+             DO  i = 1, nbgp
                 cs_advc_flags_s(:,nyn+i,:) = MERGE(                            &
                                         IBSET( cs_advc_flags_s(:,nyn,:), 31 ), &
@@ -2005,7 +2029,7 @@
           ENDIF
        ENDIF
-       IF ( decycle_chem_lr )  THEN
+       IF ( bc_dirichlet_cs_l .OR. bc_dirichlet_cs_r )  THEN
           IF ( nxl == 0  )  THEN
-             DO  i = 1, nbgp   
+             DO  i = 1, nbgp
                 cs_advc_flags_s(:,:,nxl-i) = MERGE(                            &
                                         IBSET( cs_advc_flags_s(:,:,nxl), 31 ), &
@@ -2015,6 +2039,6 @@
              ENDDO
           ENDIF
-          IF ( nxr == nx )  THEN 
-             DO  i = 1, nbgp   
+          IF ( nxr == nx )  THEN
+             DO  i = 1, nbgp
                 cs_advc_flags_s(:,:,nxr+i) = MERGE(                            &
                                         IBSET( cs_advc_flags_s(:,:,nxr), 31 ), &
@@ -2023,5 +2047,5 @@
                                                   )
              ENDDO
-          ENDIF      
+          ENDIF
        ENDIF
 !
@@ -2039,10 +2063,9 @@
 !--    oscillations, which are responsible for high concentration maxima that may
 !--    appear under shear-free stable conditions. 
-       CALL ws_init_flags_scalar(                                              &
-                  bc_dirichlet_l  .OR.  bc_radiation_l  .OR.  ( decycle_chem_lr .AND. nxl == 0  ), &
-                  bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  ( decycle_chem_ns .AND. nyn == ny ), &
-                  bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  ( decycle_chem_lr .AND. nxr == nx ), &
-                  bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  ( decycle_chem_ns .AND. nys == 0  ), &
-                  cs_advc_flags_s, .TRUE. )
+       CALL ws_init_flags_scalar( bc_dirichlet_cs_l  .OR.  bc_radiation_cs_l,                      &
+                                  bc_dirichlet_cs_n  .OR.  bc_radiation_cs_n,                      &
+                                  bc_dirichlet_cs_r  .OR.  bc_radiation_cs_r,                      &
+                                  bc_dirichlet_cs_s  .OR.  bc_radiation_cs_s,                      &
+                                  cs_advc_flags_s, .TRUE. )
     ENDIF
 !
@@ -2074,9 +2097,12 @@
 !
 !--       Transfer initial profiles to the arrays of the 3D model
+!--       Concentrations within buildings are set to zero.
           DO  lsp = 1, nspec
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
                    DO lpr_lev = 1, nz + 1
-                      chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(lpr_lev,j,i), 0 ) )
+                      chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)&
+                                                            * flag
                    ENDDO
                 ENDDO
@@ -2090,5 +2116,9 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   chem_species(lsp)%conc(:,j,i) = chem_species(lsp)%conc_pr_init    
+                   DO  lpr_lev = nzb, nz+1
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(lpr_lev,j,i), 0 ) )
+                      chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)&
+                                                            * flag
+                   ENDDO
                 ENDDO
              ENDDO
@@ -2100,6 +2130,11 @@
        IF ( cs_surface_initial_change(1) /= 0.0_wp )  THEN            
           DO  lsp = 1, nspec 
-             chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) +  &
-                                               cs_surface_initial_change(lsp)
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb,j,i), 0 ) )
+                   chem_species(lsp)%conc(nzb,j,i) = chem_species(lsp)%conc(nzb,j,i) +  &
+                                                     cs_surface_initial_change(lsp) * flag
+                ENDDO
+             ENDDO
           ENDDO
        ENDIF 
@@ -2346,4 +2381,8 @@
 
     NAMELIST /chemistry_parameters/  bc_cs_b,                          &
+         bc_cs_l,                          &
+         bc_cs_n,                          &
+         bc_cs_r,                          &
+         bc_cs_s,                          &
          bc_cs_t,                          &
          call_chem_at_all_substeps,        &
@@ -2360,7 +2399,4 @@
          cs_vertical_gradient_level,       &
          daytype_mdh,                      &
-         decycle_chem_lr,                  &
-         decycle_chem_ns,                  &            
-         decycle_method,                   &
          deposition_dry,                   &
          emissions_anthropogenic,          &
@@ -2590,5 +2626,5 @@
        solver_type = 'Rang3'
     ELSE
-       message_string = 'illegal solver type'
+       message_string = 'illegal Rosenbrock-solver type'
        CALL message( 'chem_parin', 'PA0506', 1, 2, 0, 6, 0 )
     END IF
@@ -2774,5 +2810,4 @@
 
    INTEGER(iwp) ::  lsp       !<
-   INTEGER(iwp) ::  lsp_usr   !<
  
 !
@@ -2780,21 +2815,10 @@
        CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'start' )
        DO  lsp = 1, nvar
-          CALL exchange_horiz( chem_species(lsp)%conc, nbgp )    
-          lsp_usr = 1  
-          DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
-             IF ( TRIM(chem_species(lsp)%name) == TRIM(cs_name(lsp_usr)) )  THEN
-!
-!--             As chem_exchange_horiz_bounds is called at the beginning
-!--             of prognostic_equations, boundary conditions are set on 
-!--             %conc. 
-                CALL chem_boundary_conds( chem_species(lsp)%conc,              &
-                                          chem_species(lsp)%conc_pr_init ) 
-             
-             ENDIF
-             lsp_usr = lsp_usr + 1
-          ENDDO
-
-
+          CALL exchange_horiz( chem_species(lsp)%conc, nbgp,                                       &
+                               alternative_communicator = communicator_chem )
        ENDDO
+
+       CALL chem_boundary_conditions( horizontal_conditions_only = .TRUE. )
+
        CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'stop' )
 
@@ -2832,8 +2856,8 @@
           IF ( ws_scheme_sca )  THEN
              CALL advec_s_ws( cs_advc_flags_s, chem_species(ilsp)%conc, 'kc',                      &
-                  bc_dirichlet_l  .OR.  bc_radiation_l  .OR.  ( decycle_chem_lr .AND. nxl == 0  ), &
-                  bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  ( decycle_chem_ns .AND. nyn == ny ), &
-                  bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  ( decycle_chem_lr .AND. nxr == nx ), &
-                  bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  ( decycle_chem_ns .AND. nys == 0  ) )
+                              bc_dirichlet_cs_l  .OR.  bc_radiation_cs_l,                          &
+                              bc_dirichlet_cs_n  .OR.  bc_radiation_cs_n,                          &
+                              bc_dirichlet_cs_r  .OR.  bc_radiation_cs_r,                          &
+                              bc_dirichlet_cs_s  .OR.  bc_radiation_cs_s )
           ELSE
              CALL advec_s_pw( chem_species(ilsp)%conc )
@@ -2942,5 +2966,9 @@
        IF ( timestep_scheme(1:5) == 'runge' )  THEN
           IF ( ws_scheme_sca )  THEN
-             CALL advec_s_ws( cs_advc_flags_s, i, j, chem_species(ilsp)%conc, 'kc',                &
+             CALL advec_s_ws( cs_advc_flags_s,                                                     &
+                              i,                                                                   &
+                              j,                                                                   &
+                              chem_species(ilsp)%conc,                                             &
+                              'kc',                                                                &
                               chem_species(ilsp)%flux_s_cs,                                        &
                               chem_species(ilsp)%diss_s_cs,                                        &
@@ -2949,9 +2977,9 @@
                               i_omp_start,                                                         &
                               tn,                                                                  &
-                  bc_dirichlet_l  .OR.  bc_radiation_l  .OR.  ( decycle_chem_lr .AND. nxl == 0  ), &
-                  bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  ( decycle_chem_ns .AND. nyn == ny ), &
-                  bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  ( decycle_chem_lr .AND. nxr == nx ), &
-                  bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  ( decycle_chem_ns .AND. nys == 0  ), &
-                  monotonic_limiter_z )
+                              bc_dirichlet_cs_l  .OR.  bc_radiation_cs_l,                          &
+                              bc_dirichlet_cs_n  .OR.  bc_radiation_cs_n,                          &
+                              bc_dirichlet_cs_r  .OR.  bc_radiation_cs_r,                          &
+                              bc_dirichlet_cs_s  .OR.  bc_radiation_cs_s,                          &
+                              monotonic_limiter_z )
           ELSE
              CALL advec_s_pw( i, j, chem_species(ilsp)%conc )
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 4510)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 4511)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! chemistry decycling replaced by explicit setting of lateral boundary conditions
+! 
+! 4508 2020-04-24 13:32:20Z raasch
 ! salsa decycling replaced by explicit setting of lateral boundary conditions
 !
@@ -242,14 +245,6 @@
 
     USE chem_modules,                                                                              &
-        ONLY:  bc_cs_t_val, chem_species, emissions_anthropogenic, emiss_read_legacy_mode,         &
-               n_matched_vars
-
-#if defined( __parallel )
-    USE chem_modules,                                                                              &
-        ONLY:  cs_name
-#endif
-
-    USE chemistry_model_mod,                                                                       &
-        ONLY:  chem_boundary_conds
+        ONLY:  bc_cs_t_val, chem_species, communicator_chem, emissions_anthropogenic,              &
+               emiss_read_legacy_mode, n_matched_vars
 
     USE control_parameters,                                                                        &
@@ -469,10 +464,6 @@
     INTEGER(iwp) ::  icc                 !< additional index for aerosol mass bins
     INTEGER(iwp) ::  ig                  !< index for salsa gases
-    INTEGER(iwp) ::  lsp                 !<
     INTEGER(iwp) ::  mid                 !< masked output running index
-#if defined( __parallel )
-    INTEGER(iwp) ::  lsp_usr             !<
     INTEGER(iwp) ::  n                   !< loop counter for chemistry species
-#endif
 
     REAL(wp) ::  dt_3d_old  !< temporary storage of timestep to be used for
@@ -675,7 +666,7 @@
            bc_q_t_val  = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
            IF ( air_chemistry )  THEN
-              DO  lsp = 1, nvar
-                 bc_cs_t_val = (  chem_species(lsp)%conc_pr_init(nzt+1)                            &
-                                - chem_species(lsp)%conc_pr_init(nzt) )                            &
+              DO  n = 1, nvar
+                 bc_cs_t_val = (  chem_species(n)%conc_pr_init(nzt+1)                              &
+                                - chem_species(n)%conc_pr_init(nzt) )                              &
                                / dzu(nzt+1)
               ENDDO
@@ -824,6 +815,7 @@
           IF ( passive_scalar )  CALL exchange_horiz( s_p, nbgp )
           IF ( air_chemistry )  THEN
-             DO  lsp = 1, nvar
-                CALL exchange_horiz( chem_species(lsp)%conc_p, nbgp )
+             DO  n = 1, nvar
+                CALL exchange_horiz( chem_species(n)%conc_p, nbgp,                                 &
+                                     alternative_communicator = communicator_chem )
              ENDDO
           ENDIF
@@ -932,5 +924,6 @@
                 IF ( air_chemistry )  THEN
                    DO  n = 1, nvar
-                      CALL exchange_horiz( chem_species(n)%conc, nbgp )
+                      CALL exchange_horiz( chem_species(n)%conc, nbgp,                             &
+                                           alternative_communicator = communicator_chem )
                    ENDDO
                 ENDIF
@@ -960,19 +953,4 @@
 !--          Set boundary conditions again after interpolation and anterpolation.
              CALL pmci_boundary_conds
-
-!
-!--          Set chemistry boundary conditions (decycling)
-             IF ( air_chemistry )  THEN
-                DO  lsp = 1, nvar
-                   lsp_usr = 1
-                   DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
-                      IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN
-                         CALL chem_boundary_conds( chem_species(lsp)%conc,                         &
-                                                   chem_species(lsp)%conc_pr_init )
-                      ENDIF
-                      lsp_usr = lsp_usr + 1
-                   ENDDO
-                ENDDO
-             ENDIF
 
              CALL cpu_log( log_point(60), 'nesting', 'stop' )
