Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 4507)
+++ palm/trunk/SOURCE/Makefile	(revision 4508)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# dependencies for salsa_mod updated
+# 
+# 4495 2020-04-13 20:11:20Z raasch
 # added routines and dependencies for restart data with MPI-IO
 # 
@@ -895,5 +898,4 @@
 	pmc_particle_interface.o \
 	restart_data_mpi_io_mod.o \
-	salsa_mod.o \
 	surface_layer_fluxes_mod.o \
 	surface_data_output_mod.o \
@@ -1024,5 +1026,4 @@
 	plant_canopy_model_mod.o \
 	radiation_model_mod.o \
-	salsa_mod.o \
 	subsidence_mod.o \
 	surface_mod.o \
Index: palm/trunk/SOURCE/pmc_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 4507)
+++ palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 4508)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! salsa variable name changed
+! 
+! 4444 2020-03-05 15:59:50Z raasch
 ! bugfix: cpp-directives and variable declarations for serial mode added
 ! 
@@ -288,5 +291,5 @@
 
     USE salsa_mod,                                                             &
-        ONLY:  aerosol_mass, aerosol_number, gconc_2, ibc_salsa_b,             &
+        ONLY:  aerosol_mass, aerosol_number, gconc_2, ibc_aer_b,               &
                mconc_2, nbins_aerosol,                                         &
                ncomponents_mass, nconc_2, nesting_salsa, ngases_salsa,         &
@@ -4875,5 +4878,5 @@
 !-- Set Neumann boundary conditions for aerosols and salsa gases
     IF ( salsa  .AND.  nesting_salsa )  THEN
-       IF ( ibc_salsa_b == 1 )  THEN
+       IF ( ibc_aer_b == 1 )  THEN
           DO  m = 1, bc_h(0)%ns
              ic = bc_h(0)%i(m)
Index: palm/trunk/SOURCE/salsa_mod.f90
===================================================================
--- palm/trunk/SOURCE/salsa_mod.f90	(revision 4507)
+++ palm/trunk/SOURCE/salsa_mod.f90	(revision 4508)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! decycling replaced by explicit setting of lateral boundary conditions (Siggi)
+! 
+! 4487 2020-04-03 09:38:20Z raasch
 ! bugfix for subroutine calls that contain the decycle_salsa switches as arguments
 ! 
@@ -430,6 +433,10 @@
 !
 !-- SALSA variables:
-    CHARACTER(LEN=20)  ::  bc_salsa_b = 'neumann'                 !< bottom boundary condition
-    CHARACTER(LEN=20)  ::  bc_salsa_t = 'neumann'                 !< top boundary condition
+    CHARACTER(LEN=20)  ::  bc_aer_b = 'neumann'                   !< bottom boundary condition
+    CHARACTER(LEN=20)  ::  bc_aer_l = 'undefined'                 !< left boundary condition
+    CHARACTER(LEN=20)  ::  bc_aer_n = 'undefined'                 !< north boundary condition
+    CHARACTER(LEN=20)  ::  bc_aer_r = 'undefined'                 !< right boundary condition
+    CHARACTER(LEN=20)  ::  bc_aer_s = 'undefined'                 !< south boundary condition
+    CHARACTER(LEN=20)  ::  bc_aer_t = 'neumann'                   !< top boundary condition
     CHARACTER(LEN=20)  ::  depo_pcm_par = 'zhang2001'             !< or 'petroff2010'
     CHARACTER(LEN=20)  ::  depo_pcm_type = 'deciduous_broadleaf'  !< leaf type
@@ -439,15 +446,11 @@
     CHARACTER(LEN=20)  ::  salsa_emission_mode = 'no_emission'    !< 'no_emission', 'uniform',
                                                                   !< 'parameterized', 'read_from_file'
-
-    CHARACTER(LEN=20), DIMENSION(4) ::  decycle_method_salsa =                                     &
-                                                 (/'dirichlet','dirichlet','dirichlet','dirichlet'/)
-                                     !< Decycling method at horizontal boundaries
-                                     !< 1=left, 2=right, 3=south, 4=north
-                                     !< dirichlet = initial profiles for the ghost and first 3 layers
-                                     !< neumann = zero gradient
-
     CHARACTER(LEN=3), DIMENSION(maxspec) ::  listspec = &  !< Active aerosols
                                    (/'SO4','   ','   ','   ','   ','   ','   '/)
 
+    INTEGER(iwp) ::  communicator_salsa     !< 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) ::  depo_pcm_par_num = 1   !< parametrisation type: 1=zhang2001, 2=petroff2010
     INTEGER(iwp) ::  depo_pcm_type_num = 0  !< index for the dry deposition type on the plant canopy
@@ -456,6 +459,6 @@
     INTEGER(iwp) ::  end_subrange_2a = 1    !< last index for bin subrange 2a
     INTEGER(iwp) ::  end_subrange_2b = 1    !< last index for bin subrange 2b
-    INTEGER(iwp) ::  ibc_salsa_b            !< index for the bottom boundary condition
-    INTEGER(iwp) ::  ibc_salsa_t            !< index for the top boundary condition
+    INTEGER(iwp) ::  ibc_aer_b              !< index for the bottom boundary condition
+    INTEGER(iwp) ::  ibc_aer_t              !< index for the top boundary condition
     INTEGER(iwp) ::  index_bc  = -1         !< index for black carbon (BC)
     INTEGER(iwp) ::  index_du  = -1         !< index for dust
@@ -520,6 +523,20 @@
 !-- SALSA switches:
     LOGICAL ::  advect_particle_water   = .TRUE.   !< Advect water concentration of particles
-    LOGICAL ::  decycle_salsa_lr        = .FALSE.  !< Undo cyclic boundaries: left and right
-    LOGICAL ::  decycle_salsa_ns        = .FALSE.  !< Undo cyclic boundaries: north and south
+    LOGICAL ::  bc_dirichlet_aer_l      = .FALSE.  !< flag for indicating a dirichlet condition at
+                                                   !< the left boundary
+    LOGICAL ::  bc_dirichlet_aer_n      = .FALSE.  !< flag for indicating a dirichlet condition at
+                                                   !< the north boundary
+    LOGICAL ::  bc_dirichlet_aer_r      = .FALSE.  !< flag for indicating a dirichlet condition at
+                                                   !< the right boundary
+    LOGICAL ::  bc_dirichlet_aer_s      = .FALSE.  !< flag for indicating a dirichlet condition at
+                                                   !< the south boundary
+    LOGICAL ::  bc_radiation_aer_l      = .FALSE.  !< flag for indicating a radiation/neumann
+                                                   !< condition at the left boundary
+    LOGICAL ::  bc_radiation_aer_n      = .FALSE.  !< flag for indicating a radiation/neumann
+                                                   !< condition at the north boundary
+    LOGICAL ::  bc_radiation_aer_r      = .FALSE.  !< flag for indicating a radiation/neumann
+                                                   !< condition at the right boundary
+    LOGICAL ::  bc_radiation_aer_s      = .FALSE.  !< flag for indicating a radiation/neumann
+                                                   !< condition at the south boundary
     LOGICAL ::  include_emission        = .FALSE.  !< Include or not emissions
     LOGICAL ::  feedback_to_palm        = .FALSE.  !< Allow feedback due to condensation of H2O
@@ -531,4 +548,5 @@
     LOGICAL ::  van_der_waals_coagc     = .FALSE.  !< Include van der Waals and viscous forces in coagulation
     LOGICAL ::  write_binary_salsa      = .TRUE.   !< read binary for salsa
+
 !
 !-- Process switches: nl* is read from the NAMELIST and is NOT changed.
@@ -886,9 +904,4 @@
     END INTERFACE salsa_3d_data_averaging
 
-    INTERFACE salsa_boundary_conds
-       MODULE PROCEDURE salsa_boundary_conds
-       MODULE PROCEDURE salsa_boundary_conds_decycle
-    END INTERFACE salsa_boundary_conds
-
     INTERFACE salsa_boundary_conditions
        MODULE PROCEDURE salsa_boundary_conditions
@@ -998,5 +1011,4 @@
     PUBLIC salsa_3d_data_averaging,       &
            salsa_actions,                 &
-           salsa_boundary_conds,          &
            salsa_boundary_conditions,     &
            salsa_check_data_output,       &
@@ -1029,5 +1041,6 @@
            bc_an_t_val,           &
            bc_gt_t_val,           &
-           ibc_salsa_b,           &
+           communicator_salsa,    &
+           ibc_aer_b,             &
            init_aerosol_type,     &
            init_gases_type,       &
@@ -1074,9 +1087,10 @@
                                      aerosol_flux_sigmag,                      &
                                      advect_particle_water,                    &
-                                     bc_salsa_b,                               &
-                                     bc_salsa_t,                               &
-                                     decycle_salsa_lr,                         &
-                                     decycle_method_salsa,                     &
-                                     decycle_salsa_ns,                         &
+                                     bc_aer_b,                                 &
+                                     bc_aer_l,                                 &
+                                     bc_aer_n,                                 &
+                                     bc_aer_r,                                 &
+                                     bc_aer_s,                                 &
+                                     bc_aer_t,                                 &
                                      depo_pcm_par,                             &
                                      depo_pcm_type,                            &
@@ -1175,7 +1189,7 @@
     IF ( child_domain )  THEN
        IF ( nesting_salsa )  THEN
-          bc_salsa_t = 'nested'
+          bc_aer_t = 'nested'
        ELSE
-          bc_salsa_t = 'neumann'
+          bc_aer_t = 'neumann'
        ENDIF
     ENDIF
@@ -1184,31 +1198,31 @@
     IF ( nesting_offline )  THEN
        IF ( nesting_offline_salsa )  THEN
-          bc_salsa_t = 'nesting_offline'
+          bc_aer_t = 'nesting_offline'
        ELSE
-          bc_salsa_t = 'neumann'
+          bc_aer_t = 'neumann'
        ENDIF
     ENDIF
 !
 !-- Set bottom boundary condition flag
-    IF ( bc_salsa_b == 'dirichlet' )  THEN
-       ibc_salsa_b = 0
-    ELSEIF ( bc_salsa_b == 'neumann' )  THEN
-       ibc_salsa_b = 1
+    IF ( bc_aer_b == 'dirichlet' )  THEN
+       ibc_aer_b = 0
+    ELSEIF ( bc_aer_b == 'neumann' )  THEN
+       ibc_aer_b = 1
     ELSE
-       message_string = 'unknown boundary condition: bc_salsa_b = "' // TRIM( bc_salsa_t ) // '"'
+       message_string = 'unknown boundary condition: bc_aer_b = "' // TRIM( bc_aer_t ) // '"'
        CALL message( 'salsa_check_parameters', 'PA0595', 1, 2, 0, 6, 0 )
     ENDIF
 !
 !-- Set top boundary conditions flag
-    IF ( bc_salsa_t == 'dirichlet' )  THEN
-       ibc_salsa_t = 0
-    ELSEIF ( bc_salsa_t == 'neumann' )  THEN
-       ibc_salsa_t = 1
-    ELSEIF ( bc_salsa_t == 'initial_gradient' )  THEN
-       ibc_salsa_t = 2
-    ELSEIF ( bc_salsa_t == 'nested'  .OR.  bc_salsa_t == 'nesting_offline' )  THEN
-       ibc_salsa_t = 3
+    IF ( bc_aer_t == 'dirichlet' )  THEN
+       ibc_aer_t = 0
+    ELSEIF ( bc_aer_t == 'neumann' )  THEN
+       ibc_aer_t = 1
+    ELSEIF ( bc_aer_t == 'initial_gradient' )  THEN
+       ibc_aer_t = 2
+    ELSEIF ( bc_aer_t == 'nested'  .OR.  bc_aer_t == 'nesting_offline' )  THEN
+       ibc_aer_t = 3
     ELSE
-       message_string = 'unknown boundary condition: bc_salsa_t = "' // TRIM( bc_salsa_t ) // '"'
+       message_string = 'unknown boundary condition: bc_aer_t = "' // TRIM( bc_aer_t ) // '"'
        CALL message( 'salsa_check_parameters', 'PA0596', 1, 2, 0, 6, 0 )
     ENDIF
@@ -1221,7 +1235,102 @@
 !
 !-- Check bottom boundary condition in case of surface emissions
-    IF ( salsa_emission_mode /= 'no_emission'  .AND.  ibc_salsa_b  == 0 ) THEN
-       message_string = 'salsa_emission_mode /= "no_emission" requires bc_salsa_b = "Neumann"'
+    IF ( salsa_emission_mode /= 'no_emission'  .AND.  ibc_aer_b  == 0 ) THEN
+       message_string = 'salsa_emission_mode /= "no_emission" requires bc_aer_b = "Neumann"'
        CALL message( 'salsa_check_parameters','PA0598', 1, 2, 0, 6, 0 )
+    ENDIF
+!
+!-- Check left and right boundary conditions. First set default value if not set by user.
+    IF ( bc_aer_l == 'undefined' )  THEN
+       IF ( bc_lr == 'cyclic' )  THEN
+          bc_aer_l = 'cyclic'
+       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
+          bc_aer_l = 'dirichlet'
+       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
+          bc_aer_l = 'neumann'
+       ENDIF
+    ENDIF
+    IF ( bc_aer_r == 'undefined' )  THEN
+       IF ( bc_lr == 'cyclic' )  THEN
+          bc_aer_r = 'cyclic'
+       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
+          bc_aer_r = 'neumann'
+       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
+          bc_aer_r = 'dirichlet'
+       ENDIF
+    ENDIF
+    IF ( bc_aer_l /= 'dirichlet'  .AND.  bc_aer_l /= 'neumann'  .AND.  bc_aer_l /= 'cyclic' )  THEN
+       message_string = 'unknown boundary condition: bc_aer_l = "' // TRIM( bc_aer_l ) // '"'
+       CALL message( 'salsa_check_parameters','PA0626', 1, 2, 0, 6, 0 )
+    ENDIF
+    IF ( bc_aer_r /= 'dirichlet'  .AND.  bc_aer_r /= 'neumann'  .AND.  bc_aer_r /= 'cyclic' )  THEN
+       message_string = 'unknown boundary condition: bc_aer_r = "' // TRIM( bc_aer_r ) // '"'
+       CALL message( 'salsa_check_parameters','PA0627', 1, 2, 0, 6, 0 )
+    ENDIF
+!
+!-- Check north and south boundary conditions. First set default value if not set by user.
+    IF ( bc_aer_n == 'undefined' )  THEN
+       IF ( bc_ns == 'cyclic' )  THEN
+          bc_aer_n = 'cyclic'
+       ELSEIF ( bc_ns == 'dirichlet/radiation' )  THEN
+          bc_aer_n = 'dirichlet'
+       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
+          bc_aer_n = 'neumann'
+       ENDIF
+    ENDIF
+    IF ( bc_aer_s == 'undefined' )  THEN
+       IF ( bc_ns == 'cyclic' )  THEN
+          bc_aer_s = 'cyclic'
+       ELSEIF ( bc_ns == 'dirichlet/radiation' )  THEN
+          bc_aer_s = 'neumann'
+       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
+          bc_aer_s = 'dirichlet'
+       ENDIF
+    ENDIF
+    IF ( bc_aer_n /= 'dirichlet'  .AND.  bc_aer_n /= 'neumann'  .AND.  bc_aer_n /= 'cyclic' )  THEN
+       message_string = 'unknown boundary condition: bc_aer_n = "' // TRIM( bc_aer_n ) // '"'
+       CALL message( 'salsa_check_parameters','PA0709', 1, 2, 0, 6, 0 )
+    ENDIF
+    IF ( bc_aer_s /= 'dirichlet'  .AND.  bc_aer_s /= 'neumann'  .AND.  bc_aer_s /= 'cyclic' )  THEN
+       message_string = 'unknown boundary condition: bc_aer_s = "' // TRIM( bc_aer_s ) // '"'
+       CALL message( 'salsa_check_parameters','PA0711', 1, 2, 0, 6, 0 )
+    ENDIF
+!
+!-- Cyclic conditions must be set identically at opposing boundaries
+    IF ( ( bc_aer_l == 'cyclic' .AND. bc_aer_r /= 'cyclic' )  .OR.                                 &
+         ( bc_aer_r == 'cyclic' .AND. bc_aer_l /= 'cyclic' ) )  THEN
+       message_string = 'boundary conditions bc_aer_l and bc_aer_r must both be cyclic or non-cyclic'
+       CALL message( 'salsa_check_parameters','PA0712', 1, 2, 0, 6, 0 )
+    ENDIF
+    IF ( ( bc_aer_n == 'cyclic' .AND. bc_aer_s /= 'cyclic' )  .OR.                                 &
+         ( bc_aer_s == 'cyclic' .AND. bc_aer_n /= 'cyclic' ) )  THEN
+       message_string = 'boundary conditions bc_aer_n and bc_aer_s must both be cyclic or non-cyclic'
+       CALL message( 'salsa_check_parameters','PA0713', 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_aer_n == 'dirichlet'  .AND.  nyn == ny )  bc_dirichlet_aer_n = .TRUE.
+    IF ( bc_aer_n == 'neumann'    .AND.  nyn == ny )  bc_radiation_aer_n = .TRUE.
+    IF ( bc_aer_s == 'dirichlet'  .AND.  nys ==  0 )  bc_dirichlet_aer_s = .TRUE.
+    IF ( bc_aer_s == 'neumann'    .AND.  nys ==  0 )  bc_radiation_aer_s = .TRUE.
+    IF ( bc_aer_l == 'dirichlet'  .AND.  nxl ==  0 )  bc_dirichlet_aer_l = .TRUE.
+    IF ( bc_aer_l == 'neumann'    .AND.  nxl ==  0 )  bc_radiation_aer_l = .TRUE.
+    IF ( bc_aer_r == 'dirichlet'  .AND.  nxr == nx )  bc_dirichlet_aer_r = .TRUE.
+    IF ( bc_aer_r == 'neumann'    .AND.  nxr == nx )  bc_radiation_aer_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_aer_l == 'cyclic' )  THEN
+       IF ( bc_aer_s == 'cyclic' )  THEN
+          communicator_salsa = 1
+       ELSE
+          communicator_salsa = 2
+       ENDIF
+    ELSE
+       IF ( bc_aer_s == 'cyclic' )  THEN
+          communicator_salsa = 3
+       ELSE
+          communicator_salsa = 4
+       ENDIF
     ENDIF
 !
@@ -1336,53 +1445,59 @@
        WRITE( io, 24 )
     ENDIF
+    WRITE( io, 26 )  TRIM( bc_aer_b ), TRIM( bc_aer_t ), TRIM( bc_aer_s ), TRIM( bc_aer_n ),       &
+                     TRIM( bc_aer_l ), TRIM( bc_aer_r )
 
 1   FORMAT (//' SALSA information:'/                                                               &
               ' ------------------------------'/)
-2   FORMAT   ('    Starts at: skip_time_do_salsa = ', F10.2, '  s')
-3   FORMAT  (/'    Timestep: dt_salsa = ', F6.2, '  s')
-4   FORMAT  (/'    Array shape (z,y,x,bins):'/                                                     &
-              '       aerosol_number:  ', 4(I5)) 
-5   FORMAT  (/'       aerosol_mass:    ', 4(I5),/                                                  &
-              '       (advect_particle_water = ', L1, ')')
-6   FORMAT   ('       salsa_gas: ', 4(I5),/                                                        &
-              '       (salsa_gases_from_chem = ', L1, ')')
-7   FORMAT  (/'    Aerosol dynamic processes included: ')
-8   FORMAT  (/'       nucleation (scheme = ', I1, ' and J3 parametrization = ', I1, ')')
-9   FORMAT  (/'       coagulation')
-10  FORMAT  (/'       condensation (of precursor gases = ', L1, ' and water vapour = ', L1, ')' )
-11  FORMAT  (/'       dissolutional growth by HNO3 and NH3')
-12  FORMAT  (/'       dry deposition (on vegetation = ', L1, ' and on topography = ', L1, ')')
-13  FORMAT  (/'    Aerosol bin subrange limits (in metres): ',  3(ES10.2E3), /                     &
-              '    Number of size bins for each aerosol subrange: ', 2I3,/                         &
-              '    Aerosol bin lower limits (in metres): ', 12(ES10.2E3))
-25  FORMAT  (/'    Bin geometric mean diameters (in metres): ', 12(ES10.2E3))
-14  FORMAT   ('    Initial number concentration in bins at the lowest level (#/m**3):', 9(ES10.2E3))
-15  FORMAT  (/'    Number of chemical components used: ', I1,/                                     &
-              '       Species: ',7(A6),/                                                           &
-              '    Initial relative contribution of each species to particle volume in:',/         &
-              '       a-bins: ', 7(F6.3),/                                                         &
-              '       b-bins: ', 7(F6.3))
-16  FORMAT  (/'    Number of gaseous tracers used: ', I1,/                                         &
-              '    Initial gas concentrations:',/                                                  &
-              '       H2SO4: ',ES12.4E3, ' #/m**3',/                                               &
-              '       HNO3:  ',ES12.4E3, ' #/m**3',/                                               &
-              '       NH3:   ',ES12.4E3, ' #/m**3',/                                               &
-              '       OCNV:  ',ES12.4E3, ' #/m**3',/                                               &
-              '       OCSV:  ',ES12.4E3, ' #/m**3')
-17   FORMAT (/'   Initialising concentrations: ', /                                                &
-              '      Aerosol size distribution: init_aerosol_type = ', I1,/                        &
-              '      Gas concentrations: init_gases_type = ', I1 )
-18   FORMAT ( '      Mode diametres: dpg(nmod) = ', 7(F7.3), ' (m)', /                             &
-              '      Standard deviation: sigmag(nmod) = ', 7(F7.2),/                               &
-              '      Number concentration: n_lognorm(nmod) = ', 7(ES12.4E3), ' (#/m3)' )
-19   FORMAT (/'      Size distribution read from a file.')
-20   FORMAT (/'   Nesting for salsa variables: ', L1 )
-21   FORMAT (/'   Offline nesting for salsa variables: ', L1 )
-22   FORMAT (/'   Emissions: salsa_emission_mode = ', A )
-23   FORMAT (/'      surface_aerosol_flux = ', ES12.4E3, ' #/m**2/s', /                            &
-              '      aerosol_flux_dpg     =  ', 7(F7.3), ' (m)', /                                 &
-              '      aerosol_flux_sigmag  =  ', 7(F7.2), /                                         &
-              '      aerosol_mass_fracs_a =  ', 7(ES12.4E3) )
-24   FORMAT (/'      (currently all emissions are soluble!)')
+2   FORMAT ('    Starts at: skip_time_do_salsa = ', F10.2, '  s')
+3   FORMAT (/'    Timestep: dt_salsa = ', F6.2, '  s')
+4   FORMAT (/'    Array shape (z,y,x,bins):'/                                                      &
+             '       aerosol_number:  ', 4(I5))
+5   FORMAT (/'       aerosol_mass:    ', 4(I5),/                                                   &
+             '       (advect_particle_water = ', L1, ')')
+6   FORMAT ( '       salsa_gas: ', 4(I5),/                                                         &
+             '       (salsa_gases_from_chem = ', L1, ')')
+7   FORMAT (/'    Aerosol dynamic processes included: ')
+8   FORMAT (/'       nucleation (scheme = ', I1, ' and J3 parametrization = ', I1, ')')
+9   FORMAT (/'       coagulation')
+10  FORMAT (/'       condensation (of precursor gases = ', L1, ' and water vapour = ', L1, ')' )
+11  FORMAT (/'       dissolutional growth by HNO3 and NH3')
+12  FORMAT (/'       dry deposition (on vegetation = ', L1, ' and on topography = ', L1, ')')
+13  FORMAT (/'    Aerosol bin subrange limits (in metres): ',  3(ES10.2E3), /                      &
+             '    Number of size bins for each aerosol subrange: ', 2I3,/                          &
+             '    Aerosol bin lower limits (in metres): ', 12(ES10.2E3))
+25  FORMAT (/'    Bin geometric mean diameters (in metres): ', 12(ES10.2E3))
+14  FORMAT  ('    Initial number concentration in bins at the lowest level (#/m**3):', 9(ES10.2E3))
+15  FORMAT (/'    Number of chemical components used: ', I1,/                                      &
+             '       Species: ',7(A6),/                                                            &
+             '    Initial relative contribution of each species to particle volume in:',/          &
+             '       a-bins: ', 7(F6.3),/                                                          &
+             '       b-bins: ', 7(F6.3))
+16  FORMAT (/'    Number of gaseous tracers used: ', I1,/                                          &
+             '    Initial gas concentrations:',/                                                   &
+             '       H2SO4: ',ES12.4E3, ' #/m**3',/                                                &
+             '       HNO3:  ',ES12.4E3, ' #/m**3',/                                                &
+             '       NH3:   ',ES12.4E3, ' #/m**3',/                                                &
+             '       OCNV:  ',ES12.4E3, ' #/m**3',/                                                &
+             '       OCSV:  ',ES12.4E3, ' #/m**3')
+17  FORMAT (/'   Initialising concentrations: ', /                                                 &
+             '      Aerosol size distribution: init_aerosol_type = ', I1,/                         &
+             '      Gas concentrations: init_gases_type = ', I1 )
+18  FORMAT ( '      Mode diametres: dpg(nmod) = ', 7(F7.3), ' (m)', /                              &
+             '      Standard deviation: sigmag(nmod) = ', 7(F7.2),/                                &
+             '      Number concentration: n_lognorm(nmod) = ', 7(ES12.4E3), ' (#/m3)' )
+19  FORMAT (/'      Size distribution read from a file.')
+20  FORMAT (/'   Nesting for salsa variables: ', L1 )
+21  FORMAT (/'   Offline nesting for salsa variables: ', L1 )
+22  FORMAT (/'   Emissions: salsa_emission_mode = ', A )
+23  FORMAT (/'      surface_aerosol_flux = ', ES12.4E3, ' #/m**2/s', /                             &
+             '      aerosol_flux_dpg     =  ', 7(F7.3), ' (m)', /                                  &
+             '      aerosol_flux_sigmag  =  ', 7(F7.2), /                                          &
+             '      aerosol_mass_fracs_a =  ', 7(ES12.4E3) )
+24  FORMAT (/'      (currently all emissions are soluble!)')
+26  FORMAT (/'   Boundary conditions for aerosols:', /                                             &
+             '      bottom/top:   ',A10,' / ',A10, /                                               &
+             '      north/south:  ',A10,' / ',A10, /                                               &
+             '      left/right:   ',A10,' / ',A10)
 
  END SUBROUTINE salsa_header
@@ -1680,13 +1795,7 @@
 
     ENDIF
-!
-!-- 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_salsa_lr = MERGE( decycle_salsa_lr, .FALSE., nxl == 0  .OR.  nxr == nx )
-    decycle_salsa_ns = MERGE( decycle_salsa_ns, .FALSE., nys == 0  .OR.  nyn == ny )
-!
-!-- Decycling can be applied separately for aerosol variables, while wind and other scalars may have
-!-- cyclic or nested boundary conditions. However, large gradients near the boundaries may produce
+
+!
+!-- Large gradients near the boundaries may produce
 !-- stationary numerical oscillations near the lateral boundaries when a higher-order scheme is
 !-- applied near these boundaries. To get rid-off this, set-up additional flags that control the
@@ -1694,8 +1803,8 @@
 !-- decycling.
     IF ( scalar_advec == 'ws-scheme' )  THEN
+
        ALLOCATE( salsa_advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 !
-!--    In case of decycling, set Neuman boundary conditions for wall_flags_total_0 bit 31 instead of
-!--    cyclic boundary conditions. Bit 31 is used to identify extended degradation zones (please see
+!--    Bit 31 is used to identify extended degradation zones (please see
 !--    the following comment). Note, since several also other modules may access this bit but may
 !--    have other boundary conditions, the original value of wall_flags_total_0 bit 31 must not be
@@ -1707,36 +1816,4 @@
        salsa_advc_flags_s = MERGE( IBSET( salsa_advc_flags_s, 31 ), 0, BTEST( wall_flags_total_0, 31 ) )
 
-       IF ( decycle_salsa_ns )  THEN
-          IF ( nys == 0 )  THEN
-             DO  i = 1, nbgp
-                salsa_advc_flags_s(:,nys-i,:) = MERGE( IBSET( salsa_advc_flags_s(:,nys,:), 31 ),   &
-                                                       IBCLR( salsa_advc_flags_s(:,nys,:), 31 ),   &
-                                                       BTEST( salsa_advc_flags_s(:,nys,:), 31 ) )
-             ENDDO
-          ENDIF
-          IF ( nyn == ny )  THEN
-             DO  i = 1, nbgp
-                salsa_advc_flags_s(:,nyn+i,:) = MERGE( IBSET( salsa_advc_flags_s(:,nyn,:), 31 ),   &
-                                                       IBCLR( salsa_advc_flags_s(:,nyn,:), 31 ),   &
-                                                       BTEST( salsa_advc_flags_s(:,nyn,:), 31 ) )
-             ENDDO
-          ENDIF
-       ENDIF
-       IF ( decycle_salsa_lr )  THEN
-          IF ( nxl == 0 )  THEN
-             DO  i = 1, nbgp
-                salsa_advc_flags_s(:,:,nxl-i) = MERGE( IBSET( salsa_advc_flags_s(:,:,nxl), 31 ),   &
-                                                       IBCLR( salsa_advc_flags_s(:,:,nxl), 31 ),   &
-                                                       BTEST( salsa_advc_flags_s(:,:,nxl), 31 ) )
-             ENDDO
-          ENDIF
-          IF ( nxr == nx )  THEN
-             DO  i = 1, nbgp
-                salsa_advc_flags_s(:,:,nxr+i) = MERGE( IBSET( salsa_advc_flags_s(:,:,nxr), 31 ),   &
-                                                       IBCLR( salsa_advc_flags_s(:,:,nxr), 31 ),   &
-                                                       BTEST( salsa_advc_flags_s(:,:,nxr), 31 ) )
-             ENDDO
-          ENDIF
-       ENDIF
 !
 !--    To initialise the advection flags appropriately, pass the boundary flags to
@@ -1751,10 +1828,9 @@
 !--    oscillations, which are responsible for high concentration maxima that may appear e.g. under
 !--    shear-free stable conditions.
-       CALL ws_init_flags_scalar(                                                                  &
-                 bc_dirichlet_l  .OR.  bc_radiation_l  .OR.  ( decycle_salsa_lr .AND. nxl == 0  ), &
-                 bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  ( decycle_salsa_ns .AND. nyn == ny ), &
-                 bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  ( decycle_salsa_lr .AND. nxr == nx ), &
-                 bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  ( decycle_salsa_ns .AND. nys == 0  ), &
-                 salsa_advc_flags_s, .TRUE. )
+       CALL ws_init_flags_scalar( bc_dirichlet_aer_l  .OR.  bc_radiation_aer_l,                    &
+                                  bc_dirichlet_aer_n  .OR.  bc_radiation_aer_n,                    &
+                                  bc_dirichlet_aer_r  .OR.  bc_radiation_aer_r,                    &
+                                  bc_dirichlet_aer_s  .OR.  bc_radiation_aer_s,                    &
+                                  salsa_advc_flags_s, .TRUE. )
     ENDIF
 
@@ -7902,20 +7978,23 @@
           CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
 !
-!--       Exchange ghost points and decycle if needed.
+!--       Exchange ghost points
           DO  ib = 1, nbins_aerosol
-             CALL exchange_horiz( aerosol_number(ib)%conc, nbgp )
-             CALL salsa_boundary_conds( aerosol_number(ib)%conc, aerosol_number(ib)%init )
+             CALL exchange_horiz( aerosol_number(ib)%conc, nbgp,                                   &
+                                  alternative_communicator = communicator_salsa )
              DO  ic = 1, ncomponents_mass
                 icc = ( ic - 1 ) * nbins_aerosol + ib
-                CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp )
-                CALL salsa_boundary_conds( aerosol_mass(icc)%conc, aerosol_mass(icc)%init )
+                CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp,                                 &
+                                     alternative_communicator = communicator_salsa )
              ENDDO
           ENDDO
           IF ( .NOT. salsa_gases_from_chem )  THEN
              DO  ig = 1, ngases_salsa
-                CALL exchange_horiz( salsa_gas(ig)%conc, nbgp )
-                CALL salsa_boundary_conds( salsa_gas(ig)%conc, salsa_gas(ig)%init )
+                CALL exchange_horiz( salsa_gas(ig)%conc, nbgp,                                     &
+                                     alternative_communicator = communicator_salsa )
              ENDDO
           ENDIF
+!
+!--       Apply only horizontal boundary conditions
+          CALL salsa_boundary_conditions( horizontal_conditions_only = .TRUE. )
           CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
 !
@@ -8108,11 +8187,9 @@
     IF ( timestep_scheme(1:5) == 'runge' )  THEN
        IF ( ws_scheme_sca )  THEN
-          CALL advec_s_ws( salsa_advc_flags_s, i, j, rs, id,                                       &
-                 flux_s, diss_s, flux_l, diss_l, i_omp_start, tn,                                  &
-                 bc_dirichlet_l  .OR.  bc_radiation_l  .OR.  ( decycle_salsa_lr .AND. nxl == 0  ), &
-                 bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  ( decycle_salsa_ns .AND. nyn == ny ), &
-                 bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  ( decycle_salsa_lr .AND. nxr == nx ), &
-                 bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  ( decycle_salsa_ns .AND. nys == 0  ), &
-                 monotonic_limiter_z )
+          CALL advec_s_ws( salsa_advc_flags_s, i, j, rs, id, flux_s, diss_s, flux_l, diss_l,       &
+                           i_omp_start, tn, bc_dirichlet_aer_l  .OR.  bc_radiation_aer_l,          &
+                           bc_dirichlet_aer_n  .OR.  bc_radiation_aer_n,                           &
+                           bc_dirichlet_aer_r  .OR.  bc_radiation_aer_r,                           &
+                           bc_dirichlet_aer_s  .OR.  bc_radiation_aer_s, monotonic_limiter_z )
        ELSE
           CALL advec_s_pw( i, j, rs )
@@ -8246,8 +8323,8 @@
        IF ( ws_scheme_sca )  THEN
           CALL advec_s_ws( salsa_advc_flags_s, rs, id,                                             &
-                 bc_dirichlet_l  .OR.  bc_radiation_l  .OR.  ( decycle_salsa_lr .AND. nxl == 0  ), &
-                 bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  ( decycle_salsa_ns .AND. nyn == ny ), &
-                 bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  ( decycle_salsa_lr .AND. nxr == nx ), &
-                 bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  ( decycle_salsa_ns .AND. nys == 0  ) )
+                           bc_dirichlet_aer_l  .OR.  bc_radiation_aer_l,                           &
+                           bc_dirichlet_aer_n  .OR.  bc_radiation_aer_n,                           &
+                           bc_dirichlet_aer_r  .OR.  bc_radiation_aer_r,                           &
+                           bc_dirichlet_aer_s  .OR.  bc_radiation_aer_s )
        ELSE
           CALL advec_s_pw( rs )
@@ -8338,47 +8415,7 @@
 ! Description:
 ! ------------
-!> Boundary conditions for prognostic variables in SALSA from module interface
-!------------------------------------------------------------------------------!
- SUBROUTINE salsa_boundary_conditions
-
-    IMPLICIT NONE
-
-    INTEGER(iwp) ::  ib              !< index for aerosol size bins
-    INTEGER(iwp) ::  ic              !< index for aerosol mass bins
-    INTEGER(iwp) ::  icc             !< additional index for aerosol mass bins
-    INTEGER(iwp) ::  ig              !< index for salsa gases
-
-
-!
-!-- moved from boundary_conds
-    CALL salsa_boundary_conds
-!
-!-- Boundary conditions for prognostic quantitites of other modules:
-!-- Here, only decycling is carried out
-    IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
-
-       DO  ib = 1, nbins_aerosol
-          CALL salsa_boundary_conds( aerosol_number(ib)%conc_p, aerosol_number(ib)%init )
-          DO  ic = 1, ncomponents_mass
-             icc = ( ic - 1 ) * nbins_aerosol + ib
-             CALL salsa_boundary_conds( aerosol_mass(icc)%conc_p, aerosol_mass(icc)%init )
-          ENDDO
-       ENDDO
-       IF ( .NOT. salsa_gases_from_chem )  THEN
-          DO  ig = 1, ngases_salsa
-             CALL salsa_boundary_conds( salsa_gas(ig)%conc_p, salsa_gas(ig)%init )
-          ENDDO
-       ENDIF
-
-    ENDIF
-
- END SUBROUTINE salsa_boundary_conditions
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
 !> Boundary conditions for prognostic variables in SALSA
 !------------------------------------------------------------------------------!
- SUBROUTINE salsa_boundary_conds
+ SUBROUTINE salsa_boundary_conditions( horizontal_conditions_only )
 
     USE arrays_3d,                                                                                 &
@@ -8400,341 +8437,638 @@
     INTEGER(iwp) ::  m    !< running index surface elements
 
+    LOGICAL, OPTIONAL ::  horizontal_conditions_only  !< switch to set horizontal bc only
+
+    REAL(wp) ::  flag  !< flag to mask topography grid points
+
     IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
-!
-!--    Surface conditions:
-       IF ( ibc_salsa_b == 0 )  THEN   ! Dirichlet
-!
-!--       Run loop over all non-natural and natural walls. Note, in wall-datatype the k coordinate 
-!--       belongs to the atmospheric grid point, therefore, set s_p at k-1
-          DO  l = 0, 1
-             !$OMP PARALLEL PRIVATE( ib, ic, icc, ig, i, j, k )
-             !$OMP DO
-             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)
-
-                DO  ib = 1, nbins_aerosol
-                   aerosol_number(ib)%conc_p(k+bc_h(l)%koff,j,i) =             &
-                                    aerosol_number(ib)%conc(k+bc_h(l)%koff,j,i)
-                   DO  ic = 1, ncomponents_mass
-                      icc = ( ic - 1 ) * nbins_aerosol + ib
-                      aerosol_mass(icc)%conc_p(k+bc_h(l)%koff,j,i) =           &
-                                    aerosol_mass(icc)%conc(k+bc_h(l)%koff,j,i)
-                   ENDDO
-                ENDDO
-                IF ( .NOT. salsa_gases_from_chem )  THEN
-                   DO  ig = 1, ngases_salsa
-                      salsa_gas(ig)%conc_p(k+bc_h(l)%koff,j,i) =               &
-                                    salsa_gas(ig)%conc(k+bc_h(l)%koff,j,i)
-                   ENDDO
-                ENDIF
+
+       IF ( .NOT. PRESENT( horizontal_conditions_only ) )  THEN
+!
+!--       Surface conditions:
+          IF ( ibc_aer_b == 0 )  THEN   ! Dirichlet
+!
+!--          Run loop over all non-natural and natural walls. Note, in wall-datatype the k
+!--          coordinate belongs to the atmospheric grid point, therefore, set s_p at k-1
+             DO  l = 0, 1
+                !$OMP PARALLEL PRIVATE( ib, ic, icc, ig, i, j, k )
+                !$OMP DO
+                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)
+
+                   DO  ib = 1, nbins_aerosol
+                      aerosol_number(ib)%conc_p(k+bc_h(l)%koff,j,i) =                              &
+                                                       aerosol_number(ib)%conc(k+bc_h(l)%koff,j,i)
+                      DO  ic = 1, ncomponents_mass
+                         icc = ( ic - 1 ) * nbins_aerosol + ib
+                         aerosol_mass(icc)%conc_p(k+bc_h(l)%koff,j,i) =                            &
+                                                       aerosol_mass(icc)%conc(k+bc_h(l)%koff,j,i)
+                      ENDDO
+                   ENDDO
+                   IF ( .NOT. salsa_gases_from_chem )  THEN
+                      DO  ig = 1, ngases_salsa
+                         salsa_gas(ig)%conc_p(k+bc_h(l)%koff,j,i) =                                &
+                                                       salsa_gas(ig)%conc(k+bc_h(l)%koff,j,i)
+                      ENDDO
+                   ENDIF
+
+                ENDDO
+                !$OMP END PARALLEL
 
              ENDDO
-             !$OMP END PARALLEL
-
-          ENDDO
-
-       ELSE   ! Neumann
-
-          DO l = 0, 1
-             !$OMP PARALLEL PRIVATE( ib, ic, icc, ig, i, j, k )
-             !$OMP DO
-             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)
-
-                DO  ib = 1, nbins_aerosol
-                   aerosol_number(ib)%conc_p(k+bc_h(l)%koff,j,i) =             &
-                                               aerosol_number(ib)%conc_p(k,j,i)
-                   DO  ic = 1, ncomponents_mass
-                      icc = ( ic - 1 ) * nbins_aerosol + ib
-                      aerosol_mass(icc)%conc_p(k+bc_h(l)%koff,j,i) =           &
-                                               aerosol_mass(icc)%conc_p(k,j,i)
-                   ENDDO
-                ENDDO
-                IF ( .NOT. salsa_gases_from_chem ) THEN
-                   DO  ig = 1, ngases_salsa
-                      salsa_gas(ig)%conc_p(k+bc_h(l)%koff,j,i) =               &
-                                               salsa_gas(ig)%conc_p(k,j,i)
-                   ENDDO
-                ENDIF
-
+
+          ELSE   ! Neumann
+
+             DO  l = 0, 1
+                !$OMP PARALLEL PRIVATE( ib, ic, icc, ig, i, j, k )
+                !$OMP DO
+                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)
+
+                   DO  ib = 1, nbins_aerosol
+                      aerosol_number(ib)%conc_p(k+bc_h(l)%koff,j,i) =                              &
+                                                                  aerosol_number(ib)%conc_p(k,j,i)
+                      DO  ic = 1, ncomponents_mass
+                         icc = ( ic - 1 ) * nbins_aerosol + ib
+                         aerosol_mass(icc)%conc_p(k+bc_h(l)%koff,j,i) =                            &
+                                                                  aerosol_mass(icc)%conc_p(k,j,i)
+                      ENDDO
+                   ENDDO
+                   IF ( .NOT. salsa_gases_from_chem ) THEN
+                      DO  ig = 1, ngases_salsa
+                         salsa_gas(ig)%conc_p(k+bc_h(l)%koff,j,i) = salsa_gas(ig)%conc_p(k,j,i)
+                      ENDDO
+                   ENDIF
+
+                ENDDO
+                !$OMP END PARALLEL
              ENDDO
-             !$OMP END PARALLEL
-          ENDDO
-
-       ENDIF
-!
-!--   Top boundary conditions:
-       IF ( ibc_salsa_t == 0 )  THEN   ! Dirichlet
-
-          DO  ib = 1, nbins_aerosol
-             aerosol_number(ib)%conc_p(nzt+1,:,:) = aerosol_number(ib)%conc(nzt+1,:,:)
-             DO  ic = 1, ncomponents_mass
-                icc = ( ic - 1 ) * nbins_aerosol + ib
-                aerosol_mass(icc)%conc_p(nzt+1,:,:) = aerosol_mass(icc)%conc(nzt+1,:,:)
+
+          ENDIF
+!
+!--       Top boundary conditions:
+          IF ( ibc_aer_t == 0 )  THEN   ! Dirichlet
+
+             DO  ib = 1, nbins_aerosol
+                aerosol_number(ib)%conc_p(nzt+1,:,:) = aerosol_number(ib)%conc(nzt+1,:,:)
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   aerosol_mass(icc)%conc_p(nzt+1,:,:) = aerosol_mass(icc)%conc(nzt+1,:,:)
+                ENDDO
              ENDDO
-          ENDDO
-          IF ( .NOT. salsa_gases_from_chem )  THEN
-             DO  ig = 1, ngases_salsa
-                salsa_gas(ig)%conc_p(nzt+1,:,:) = salsa_gas(ig)%conc(nzt+1,:,:)
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   salsa_gas(ig)%conc_p(nzt+1,:,:) = salsa_gas(ig)%conc(nzt+1,:,:)
+                ENDDO
+             ENDIF
+
+          ELSEIF ( ibc_aer_t == 1 )  THEN   ! Neumann
+
+             DO  ib = 1, nbins_aerosol
+                aerosol_number(ib)%conc_p(nzt+1,:,:) = aerosol_number(ib)%conc_p(nzt,:,:)
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   aerosol_mass(icc)%conc_p(nzt+1,:,:) = aerosol_mass(icc)%conc_p(nzt,:,:)
+                ENDDO
              ENDDO
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   salsa_gas(ig)%conc_p(nzt+1,:,:) = salsa_gas(ig)%conc_p(nzt,:,:)
+                ENDDO
+             ENDIF
+
+          ELSEIF ( ibc_aer_t == 2 )  THEN   ! Initial gradient
+
+             DO  ib = 1, nbins_aerosol
+                aerosol_number(ib)%conc_p(nzt+1,:,:) = aerosol_number(ib)%conc_p(nzt,:,:) +        &
+                                                       bc_an_t_val(ib) * dzu(nzt+1)
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   aerosol_mass(icc)%conc_p(nzt+1,:,:) = aerosol_mass(icc)%conc_p(nzt,:,:) +       &
+                                                         bc_am_t_val(icc) * dzu(nzt+1)
+                ENDDO
+             ENDDO
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   salsa_gas(ig)%conc_p(nzt+1,:,:) = salsa_gas(ig)%conc_p(nzt,:,:) +               &
+                                                     bc_gt_t_val(ig) * dzu(nzt+1)
+                ENDDO
+             ENDIF
+
+          ENDIF  ! top conditions
+
+!
+!--       Lateral boundary conditions
+          IF ( bc_dirichlet_aer_s )  THEN
+
+             DO  ib = 1, nbins_aerosol
+                DO  i = nxlg, nxrg
+                   DO  k = nzb+1, nzt
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nys-1,i), 0 ) )
+                      aerosol_number(ib)%conc_p(k,nys-1,i) = aerosol_number(ib)%init(k) * flag
+                   ENDDO
+                ENDDO
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   DO  i = nxlg, nxrg
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nys-1,i), 0 ) )
+                         aerosol_mass(icc)%conc_p(k,nys-1,i) = aerosol_mass(icc)%init(k) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDDO
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   DO  i = nxlg, nxrg
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nys-1,i), 0 ) )
+                         salsa_gas(ig)%conc_p(k,nys-1,i) = salsa_gas(ig)%init(k) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          ELSEIF ( bc_radiation_aer_s )  THEN
+
+             DO  ib = 1, nbins_aerosol
+                DO  i = nxlg, nxrg
+                   DO  k = nzb+1, nzt
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nys-1,i), 0 ) )
+                      aerosol_number(ib)%conc_p(k,nys-1,i) = aerosol_number(ib)%conc_p(k,nys,i) * flag
+                   ENDDO
+                ENDDO
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   DO  i = nxlg, nxrg
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nys-1,i), 0 ) )
+                         aerosol_mass(icc)%conc_p(k,nys-1,i) = aerosol_mass(icc)%conc_p(k,nys,i) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDDO
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   DO  i = nxlg, nxrg
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nys-1,i), 0 ) )
+                         salsa_gas(ig)%conc_p(k,nys-1,i) = salsa_gas(ig)%conc_p(k,nys,i) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
           ENDIF
 
-       ELSEIF ( ibc_salsa_t == 1 )  THEN   ! Neumann
-
-          DO  ib = 1, nbins_aerosol
-             aerosol_number(ib)%conc_p(nzt+1,:,:) = aerosol_number(ib)%conc_p(nzt,:,:)
-             DO  ic = 1, ncomponents_mass
-                icc = ( ic - 1 ) * nbins_aerosol + ib
-                aerosol_mass(icc)%conc_p(nzt+1,:,:) = aerosol_mass(icc)%conc_p(nzt,:,:)
+          IF ( bc_dirichlet_aer_n )  THEN
+
+             DO  ib = 1, nbins_aerosol
+                DO  i = nxlg, nxrg
+                   DO  k = nzb+1, nzt
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 0 ) )
+                      aerosol_number(ib)%conc_p(k,nyn+1,i) = aerosol_number(ib)%init(k) * flag
+                   ENDDO
+                ENDDO
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   DO  i = nxlg, nxrg
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 0 ) )
+                         aerosol_mass(icc)%conc_p(k,nyn+1,i) = aerosol_mass(icc)%init(k) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
              ENDDO
-          ENDDO
-          IF ( .NOT. salsa_gases_from_chem )  THEN
-             DO  ig = 1, ngases_salsa
-                salsa_gas(ig)%conc_p(nzt+1,:,:) = salsa_gas(ig)%conc_p(nzt,:,:)
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   DO  i = nxlg, nxrg
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 0 ) )
+                         salsa_gas(ig)%conc_p(k,nyn+1,i) = salsa_gas(ig)%init(k) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          ELSEIF ( bc_radiation_aer_n )  THEN
+
+             DO  ib = 1, nbins_aerosol
+                DO  i = nxlg, nxrg
+                   DO  k = nzb+1, nzt
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 0 ) )
+                      aerosol_number(ib)%conc_p(k,nyn+1,i) = aerosol_number(ib)%conc_p(k,nyn,i) * flag
+                   ENDDO
+                ENDDO
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   DO  i = nxlg, nxrg
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 0 ) )
+                         aerosol_mass(icc)%conc_p(k,nyn+1,i) = aerosol_mass(icc)%conc_p(k,nyn,i) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
              ENDDO
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   DO  i = nxlg, nxrg
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 0 ) )
+                         salsa_gas(ig)%conc_p(k,nyn+1,i) = salsa_gas(ig)%conc_p(k,nyn,i) * flag
+                     ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
           ENDIF
 
-       ELSEIF ( ibc_salsa_t == 2 )  THEN   ! Initial gradient
-
-          DO  ib = 1, nbins_aerosol
-             aerosol_number(ib)%conc_p(nzt+1,:,:) = aerosol_number(ib)%conc_p(nzt,:,:) +           &
-                                                    bc_an_t_val(ib) * dzu(nzt+1)
-             DO  ic = 1, ncomponents_mass
-                icc = ( ic - 1 ) * nbins_aerosol + ib
-                aerosol_mass(icc)%conc_p(nzt+1,:,:) = aerosol_mass(icc)%conc_p(nzt,:,:) +          &
-                                                      bc_am_t_val(icc) * dzu(nzt+1)
-             ENDDO
-          ENDDO
-          IF ( .NOT. salsa_gases_from_chem )  THEN
-             DO  ig = 1, ngases_salsa
-                salsa_gas(ig)%conc_p(nzt+1,:,:) = salsa_gas(ig)%conc_p(nzt,:,:) +                  &
-                                                  bc_gt_t_val(ig) * dzu(nzt+1)
-             ENDDO
-          ENDIF
-
-       ENDIF
-!
-!--    Lateral boundary conditions at the outflow
-       IF ( bc_radiation_s )  THEN
-          DO  ib = 1, nbins_aerosol
-             aerosol_number(ib)%conc_p(:,nys-1,:) = aerosol_number(ib)%conc_p(:,nys,:)
-             DO  ic = 1, ncomponents_mass
-                icc = ( ic - 1 ) * nbins_aerosol + ib
-                aerosol_mass(icc)%conc_p(:,nys-1,:) = aerosol_mass(icc)%conc_p(:,nys,:)
-             ENDDO
-          ENDDO
-          IF ( .NOT. salsa_gases_from_chem )  THEN
-             DO  ig = 1, ngases_salsa
-                salsa_gas(ig)%conc_p(:,nys-1,:) = salsa_gas(ig)%conc_p(:,nys,:)
-             ENDDO
-          ENDIF
-
-       ELSEIF ( bc_radiation_n )  THEN
-          DO  ib = 1, nbins_aerosol
-             aerosol_number(ib)%conc_p(:,nyn+1,:) = aerosol_number(ib)%conc_p(:,nyn,:)
-             DO  ic = 1, ncomponents_mass
-                icc = ( ic - 1 ) * nbins_aerosol + ib
-                aerosol_mass(icc)%conc_p(:,nyn+1,:) = aerosol_mass(icc)%conc_p(:,nyn,:)
-             ENDDO
-          ENDDO
-          IF ( .NOT. salsa_gases_from_chem )  THEN
-             DO  ig = 1, ngases_salsa
-                salsa_gas(ig)%conc_p(:,nyn+1,:) = salsa_gas(ig)%conc_p(:,nyn,:)
-             ENDDO
-          ENDIF
-
-       ELSEIF ( bc_radiation_l )  THEN
-          DO  ib = 1, nbins_aerosol
-             aerosol_number(ib)%conc_p(:,:,nxl-1) = aerosol_number(ib)%conc_p(:,:,nxl)
-             DO  ic = 1, ncomponents_mass
-                icc = ( ic - 1 ) * nbins_aerosol + ib
-                aerosol_mass(icc)%conc_p(:,:,nxl-1) = aerosol_mass(icc)%conc_p(:,:,nxl)
-             ENDDO
-          ENDDO
-          IF ( .NOT. salsa_gases_from_chem )  THEN
-             DO  ig = 1, ngases_salsa
-                salsa_gas(ig)%conc_p(:,:,nxl-1) = salsa_gas(ig)%conc_p(:,:,nxl)
-             ENDDO
-          ENDIF
-
-       ELSEIF ( bc_radiation_r )  THEN
-          DO  ib = 1, nbins_aerosol
-             aerosol_number(ib)%conc_p(:,:,nxr+1) = aerosol_number(ib)%conc_p(:,:,nxr)
-             DO  ic = 1, ncomponents_mass
-                icc = ( ic - 1 ) * nbins_aerosol + ib
-                aerosol_mass(icc)%conc_p(:,:,nxr+1) = aerosol_mass(icc)%conc_p(:,:,nxr)
-             ENDDO
-          ENDDO
-          IF ( .NOT. salsa_gases_from_chem )  THEN
-             DO  ig = 1, ngases_salsa
-                salsa_gas(ig)%conc_p(:,:,nxr+1) = salsa_gas(ig)%conc_p(:,:,nxr)
-             ENDDO
-          ENDIF
-
-       ENDIF
-
-    ENDIF
-
- END SUBROUTINE salsa_boundary_conds
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-! Undoing of the previously done cyclic boundary conditions.
-!------------------------------------------------------------------------------!
- SUBROUTINE salsa_boundary_conds_decycle ( sq, sq_init )
-
-    USE control_parameters,                                                                        &
-        ONLY:  nesting_offline
-
-    IMPLICIT NONE
-
-    INTEGER(iwp) ::  boundary  !<
-    INTEGER(iwp) ::  ee        !<
-    INTEGER(iwp) ::  copied    !<
-    INTEGER(iwp) ::  i         !<
-    INTEGER(iwp) ::  j         !<
-    INTEGER(iwp) ::  k         !<
-    INTEGER(iwp) ::  ss        !<
-
-    REAL(wp) ::  flag  !< flag to mask topography grid points
-
-    REAL(wp), DIMENSION(nzb:nzt+1) ::  sq_init  !< initial concentration profile
-
-    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sq  !< concentration array
-
-    flag = 0.0_wp
-!
-!-- Skip input if forcing from a larger-scale models is applied.
-    IF ( nesting_offline  .AND.  nesting_offline_salsa )  RETURN
-!
-!-- Left and right boundaries
-    IF ( decycle_salsa_lr  .AND.  ( bc_lr_cyc  .OR. bc_lr == 'nested' ) )  THEN
-
-       DO  boundary = 1, 2
-
-          IF ( decycle_method_salsa(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
+          IF ( bc_dirichlet_aer_l )  THEN
+
+             DO  ib = 1, nbins_aerosol
                 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 ) )
-                      sq(k,j,i) = sq_init(k) * flag
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxl-1), 0 ) )
+                      aerosol_number(ib)%conc_p(k,j,nxl-1) = aerosol_number(ib)%init(k) * flag
+                   ENDDO
+                ENDDO
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   DO  j = nysg, nyng
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxl-1), 0 ) )
+                         aerosol_mass(icc)%conc_p(k,j,nxl-1) = aerosol_mass(icc)%init(k) * flag
+                      ENDDO
                    ENDDO
                 ENDDO
              ENDDO
-
-          ELSEIF ( decycle_method_salsa(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
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   DO  j = nysg, nyng
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxl-1), 0 ) )
+                         salsa_gas(ig)%conc_p(k,j,nxl-1) = salsa_gas(ig)%init(k) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
              ENDIF
 
-              DO  i = ss, ee
+          ELSEIF ( bc_radiation_aer_l )  THEN
+
+             DO  ib = 1, nbins_aerosol
                 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 ) )
-                      sq(k,j,i) = sq(k,j,copied) * flag
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxl-1), 0 ) )
+                      aerosol_number(ib)%conc_p(k,j,nxl-1) = aerosol_number(ib)%conc_p(k,j,nxl) * flag
+                   ENDDO
+                ENDDO
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   DO  j = nysg, nyng
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxl-1), 0 ) )
+                         aerosol_mass(icc)%conc_p(k,j,nxl-1) = aerosol_mass(icc)%conc_p(k,j,nxl) * flag
+                      ENDDO
                    ENDDO
                 ENDDO
              ENDDO
-
-          ELSE
-             WRITE(message_string,*) 'unknown decycling method: decycle_method_salsa (', boundary, &
-                                     ') ="' // TRIM( decycle_method_salsa(boundary) ) // '"'
-             CALL message( 'salsa_boundary_conds_decycle', 'PA0626', 1, 2, 0, 6, 0 )
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   DO  j = nysg, nyng
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxl-1), 0 ) )
+                         salsa_gas(ig)%conc_p(k,j,nxl-1) = salsa_gas(ig)%conc_p(k,j,nxl) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
           ENDIF
-       ENDDO
-    ENDIF
-
-!
-!-- South and north boundaries
-     IF ( decycle_salsa_ns  .AND.  ( bc_ns_cyc  .OR. bc_ns == 'nested' ) )  THEN
-
-       DO  boundary = 3, 4
-
-          IF ( decycle_method_salsa(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
+
+          IF ( bc_dirichlet_aer_r )  THEN
+
+             DO  ib = 1, nbins_aerosol
+                DO  j = nysg, nyng
+                   DO  k = nzb+1, nzt
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 0 ) )
+                      aerosol_number(ib)%conc_p(k,j,nxr+1) = aerosol_number(ib)%init(k) * flag
+                   ENDDO
+                ENDDO
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   DO  j = nysg, nyng
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 0 ) )
+                         aerosol_mass(icc)%conc_p(k,j,nxr+1) = aerosol_mass(icc)%init(k) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDDO
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   DO  j = nysg, nyng
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 0 ) )
+                         salsa_gas(ig)%conc_p(k,j,nxr+1) = salsa_gas(ig)%init(k) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
              ENDIF
 
-             DO  i = nxlg, nxrg
-                DO  j = ss, ee
+          ELSEIF ( bc_radiation_aer_r )  THEN
+
+             DO  ib = 1, nbins_aerosol
+                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 ) )
-                      sq(k,j,i) = sq_init(k) * flag
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 0 ) )
+                      aerosol_number(ib)%conc_p(k,j,nxr+1) = aerosol_number(ib)%conc_p(k,j,nxr) * flag
+                   ENDDO
+                ENDDO
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   DO  j = nysg, nyng
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 0 ) )
+                         aerosol_mass(icc)%conc_p(k,j,nxr+1) = aerosol_mass(icc)%conc_p(k,j,nxr) * flag
+                      ENDDO
                    ENDDO
                 ENDDO
              ENDDO
-
-          ELSEIF ( decycle_method_salsa(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
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   DO  j = nysg, nyng
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 0 ) )
+                         salsa_gas(ig)%conc_p(k,j,nxr+1) = salsa_gas(ig)%conc_p(k,j,nxr) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
              ENDIF
 
-              DO  i = nxlg, nxrg
-                DO  j = ss, ee
+          ENDIF
+
+       ELSE
+!
+!--       Conditions for timelevel t (if called from salsa_exchange_horiz_bounds).
+!--       Only horizontal boundary conditions are set here.
+          IF ( bc_dirichlet_aer_s )  THEN
+
+             DO  ib = 1, nbins_aerosol
+                DO  i = nxlg, nxrg
                    DO  k = nzb+1, nzt
-                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
-                      sq(k,j,i) = sq(k,copied,i) * flag
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nys-1,i), 0 ) )
+                      aerosol_number(ib)%conc(k,nys-1,i) = aerosol_number(ib)%init(k) * flag
+                   ENDDO
+                ENDDO
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   DO  i = nxlg, nxrg
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nys-1,i), 0 ) )
+                         aerosol_mass(icc)%conc(k,nys-1,i) = aerosol_mass(icc)%init(k) * flag
+                      ENDDO
                    ENDDO
                 ENDDO
              ENDDO
-
-          ELSE
-             WRITE(message_string,*) 'unknown decycling method: decycle_method_salsa (', boundary, &
-                                     ') ="' // TRIM( decycle_method_salsa(boundary) ) // '"'
-             CALL message( 'salsa_boundary_conds_decycle', 'PA0627', 1, 2, 0, 6, 0 )
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   DO  i = nxlg, nxrg
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nys-1,i), 0 ) )
+                         salsa_gas(ig)%conc(k,nys-1,i) = salsa_gas(ig)%init(k) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          ELSEIF ( bc_radiation_aer_s )  THEN
+
+             DO  ib = 1, nbins_aerosol
+                DO  i = nxlg, nxrg
+                   DO  k = nzb+1, nzt
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nys-1,i), 0 ) )
+                      aerosol_number(ib)%conc(k,nys-1,i) = aerosol_number(ib)%conc(k,nys,i) * flag
+                   ENDDO
+                ENDDO
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   DO  i = nxlg, nxrg
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nys-1,i), 0 ) )
+                         aerosol_mass(icc)%conc(k,nys-1,i) = aerosol_mass(icc)%conc(k,nys,i) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDDO
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   DO  i = nxlg, nxrg
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nys-1,i), 0 ) )
+                         salsa_gas(ig)%conc(k,nys-1,i) = salsa_gas(ig)%conc(k,nys,i) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
           ENDIF
-       ENDDO
-    ENDIF
-
- END SUBROUTINE salsa_boundary_conds_decycle
+
+          IF ( bc_dirichlet_aer_n )  THEN
+
+             DO  ib = 1, nbins_aerosol
+                DO  i = nxlg, nxrg
+                   DO  k = nzb+1, nzt
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 0 ) )
+                      aerosol_number(ib)%conc(k,nyn+1,i) = aerosol_number(ib)%init(k) * flag
+                   ENDDO
+                ENDDO
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   DO  i = nxlg, nxrg
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 0 ) )
+                         aerosol_mass(icc)%conc(k,nyn+1,i) = aerosol_mass(icc)%init(k) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDDO
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   DO  i = nxlg, nxrg
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 0 ) )
+                         salsa_gas(ig)%conc(k,nyn+1,i) = salsa_gas(ig)%init(k) * flag
+                     ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          ELSEIF ( bc_radiation_aer_n )  THEN
+
+             DO  ib = 1, nbins_aerosol
+                DO  i = nxlg, nxrg
+                   DO  k = nzb+1, nzt
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 0 ) )
+                      aerosol_number(ib)%conc(k,nyn+1,i) = aerosol_number(ib)%conc(k,nyn,i) * flag
+                   ENDDO
+                ENDDO
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   DO  i = nxlg, nxrg
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 0 ) )
+                         aerosol_mass(icc)%conc(k,nyn+1,i) = aerosol_mass(icc)%conc(k,nyn,i) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDDO
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   DO  i = nxlg, nxrg
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 0 ) )
+                         salsa_gas(ig)%conc(k,nyn+1,i) = salsa_gas(ig)%conc(k,nyn,i) * flag
+                     ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          ENDIF
+
+          IF ( bc_dirichlet_aer_l )  THEN
+
+             DO  ib = 1, nbins_aerosol
+                DO  j = nysg, nyng
+                   DO  k = nzb+1, nzt
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxl-1), 0 ) )
+                      aerosol_number(ib)%conc(k,j,nxl-1) = aerosol_number(ib)%init(k) * flag
+                   ENDDO
+                ENDDO
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   DO  j = nysg, nyng
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxl-1), 0 ) )
+                         aerosol_mass(icc)%conc(k,j,nxl-1) = aerosol_mass(icc)%init(k) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDDO
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   DO  j = nysg, nyng
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxl-1), 0 ) )
+                         salsa_gas(ig)%conc(k,j,nxl-1) = salsa_gas(ig)%init(k) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          ELSEIF ( bc_radiation_aer_l )  THEN
+
+             DO  ib = 1, nbins_aerosol
+                DO  j = nysg, nyng
+                   DO  k = nzb+1, nzt
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxl-1), 0 ) )
+                      aerosol_number(ib)%conc(k,j,nxl-1) = aerosol_number(ib)%conc(k,j,nxl) * flag
+                   ENDDO
+                ENDDO
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   DO  j = nysg, nyng
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxl-1), 0 ) )
+                         aerosol_mass(icc)%conc(k,j,nxl-1) = aerosol_mass(icc)%conc(k,j,nxl) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDDO
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   DO  j = nysg, nyng
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxl-1), 0 ) )
+                         salsa_gas(ig)%conc(k,j,nxl-1) = salsa_gas(ig)%conc(k,j,nxl) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          ENDIF
+
+          IF ( bc_dirichlet_aer_r )  THEN
+
+             DO  ib = 1, nbins_aerosol
+                DO  j = nysg, nyng
+                   DO  k = nzb+1, nzt
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 0 ) )
+                      aerosol_number(ib)%conc(k,j,nxr+1) = aerosol_number(ib)%init(k) * flag
+                   ENDDO
+                ENDDO
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   DO  j = nysg, nyng
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 0 ) )
+                         aerosol_mass(icc)%conc(k,j,nxr+1) = aerosol_mass(icc)%init(k) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDDO
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   DO  j = nysg, nyng
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 0 ) )
+                         salsa_gas(ig)%conc(k,j,nxr+1) = salsa_gas(ig)%init(k) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          ELSEIF ( bc_radiation_aer_r )  THEN
+
+             DO  ib = 1, nbins_aerosol
+                DO  j = nysg, nyng
+                   DO  k = nzb+1, nzt
+                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 0 ) )
+                      aerosol_number(ib)%conc(k,j,nxr+1) = aerosol_number(ib)%conc(k,j,nxr) * flag
+                   ENDDO
+                ENDDO
+                DO  ic = 1, ncomponents_mass
+                   icc = ( ic - 1 ) * nbins_aerosol + ib
+                   DO  j = nysg, nyng
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 0 ) )
+                         aerosol_mass(icc)%conc(k,j,nxr+1) = aerosol_mass(icc)%conc(k,j,nxr) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDDO
+             IF ( .NOT. salsa_gases_from_chem )  THEN
+                DO  ig = 1, ngases_salsa
+                   DO  j = nysg, nyng
+                      DO  k = nzb+1, nzt
+                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 0 ) )
+                         salsa_gas(ig)%conc(k,j,nxr+1) = salsa_gas(ig)%conc(k,j,nxr) * flag
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          ENDIF
+
+       ENDIF  ! PRESENT( horizontal conditions only )
+
+    ENDIF
+
+ END SUBROUTINE salsa_boundary_conditions
+
 
 !------------------------------------------------------------------------------!
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 4507)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 4508)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! salsa decycling replaced by explicit setting of lateral boundary conditions
+!
+! 4502 2020-04-17 16:14:16Z schwenkel
 ! Implementation of ice microphysics
 !
@@ -213,11 +216,10 @@
 
     USE arrays_3d,                                                                                 &
-        ONLY:  diss, diss_p, dzu, e_p, nc_p, nr_p, prho, pt, pt_p, pt_init, q, qc_p, qr_p, q_init, &
-               q_p, ref_state, rho_ocean, sa_p, s_p, tend, u, u_p, v, vpt, v_p, w_p,               &
-               qi_p, ni_p
+        ONLY:  diss, diss_p, dzu, e_p, nc_p, ni_p, nr_p, prho, pt, pt_p, pt_init, q, qc_p, qr_p,   &
+               q_init, q_p, qi_p, ref_state, rho_ocean, sa_p, s_p, tend, u, u_p, v, vpt, v_p, w_p
 
 #if defined( __parallel )  &&  ! defined( _OPENACC )
     USE arrays_3d,                                                                                 &
-        ONLY:  e, nc, nr, qc, qr, s, w, qi, ni
+        ONLY:  e, nc, ni, nr, qc, qi, qr, s, w
 #endif
 
@@ -228,5 +230,5 @@
     USE bulk_cloud_model_mod,                                                                      &
         ONLY: bulk_cloud_model, calc_liquid_water_content, collision_turbulence,                   &
-              microphysics_morrison, microphysics_seifert, microphysics_ice_extension
+              microphysics_ice_extension, microphysics_morrison, microphysics_seifert
 
     USE calc_mean_profile_mod,                                                                     &
@@ -274,7 +276,7 @@
                time_dopr_listing, time_dopts, time_dosp, time_dosp_av, time_dots, time_do_av,      &
                time_do_sla, time_disturb, time_run_control, time_since_reference_point,            &
-               turbulent_inflow, turbulent_outflow, urban_surface,                                 &
+               timestep_count, turbulent_inflow, turbulent_outflow, urban_surface,                 &
                use_initial_profile_as_reference, use_single_reference_value, u_gtrans, v_gtrans,   &
-               virtual_flight, virtual_measurement, ws_scheme_mom, ws_scheme_sca, timestep_count
+               virtual_flight, virtual_measurement, ws_scheme_mom, ws_scheme_sca
 
 #if defined( __parallel )
@@ -365,6 +367,7 @@
     USE salsa_mod,                                                                                 &
         ONLY: aerosol_number, aerosol_mass, bc_am_t_val, bc_an_t_val, bc_gt_t_val,                 &
-              nbins_aerosol, ncomponents_mass, ngases_salsa, salsa_boundary_conds,                 &
-              salsa_emission_update, salsa_gas, salsa_gases_from_chem, skip_time_do_salsa
+              communicator_salsa, nbins_aerosol, ncomponents_mass, ngases_salsa,                   &
+              salsa_boundary_conditions, salsa_emission_update, salsa_gas, salsa_gases_from_chem,  &
+              skip_time_do_salsa
 
     USE spectra_mod,                                                                               &
@@ -434,7 +437,7 @@
                flux_s_w,                                                                           &
                heatflux_output_conversion,                                                         &
-               kh, km, momentumflux_output_conversion, nc, nr, p, ptdf_x, ptdf_y, qc, qr, rdf,     &
+               kh, km, momentumflux_output_conversion, nc, ni, nr, p, ptdf_x, ptdf_y, qc, qi, qr, rdf,     &
                rdf_sc, rho_air, rho_air_zw, s, tdiss_m, te_m, tpt_m, tu_m, tv_m, tw_m, ug, u_init, &
-               u_stokes_zu, vg, v_init, v_stokes_zu, w, zu, qi, ni
+               u_stokes_zu, vg, v_init, v_stokes_zu, w, zu
 
     USE control_parameters,                                                                        &
@@ -448,6 +451,6 @@
                sums_wsus_ws_l, sums_vs2_ws_l, sums_wsvs_ws_l, sums_ws2_ws_l, sums_wspts_ws_l,      &
                sums_wsqs_ws_l, sums_wssas_ws_l, sums_wsqcs_ws_l, sums_wsqrs_ws_l, sums_wsncs_ws_l, &
-               sums_wsnrs_ws_l, sums_wsss_ws_l, weight_substep, sums_salsa_ws_l,                   &
-               sums_wsqis_ws_l, sums_wsnis_ws_l
+               sums_wsnrs_ws_l, sums_wsss_ws_l, weight_substep, sums_salsa_ws_l, sums_wsqis_ws_l,  &
+               sums_wsnis_ws_l
 
     USE surface_mod,                                                                               &
@@ -703,4 +706,5 @@
        CALL module_interface_actions( 'before_timestep' )
 
+!
 !--    Start of intermediate step loop
        intermediate_timestep_count = 0
@@ -770,4 +774,5 @@
              CALL prognostic_equations_vector
           ENDIF
+
 !
 !--       Movement of agents in multi agent system
@@ -826,16 +831,20 @@
           IF ( salsa  .AND.  time_since_reference_point >= skip_time_do_salsa )  THEN
              DO  ib = 1, nbins_aerosol
-                CALL exchange_horiz( aerosol_number(ib)%conc_p, nbgp )
+                CALL exchange_horiz( aerosol_number(ib)%conc_p, nbgp,                              &
+                                     alternative_communicator = communicator_salsa )
                 DO  ic = 1, ncomponents_mass
                    icc = ( ic - 1 ) * nbins_aerosol + ib
-                   CALL exchange_horiz( aerosol_mass(icc)%conc_p, nbgp )
+                   CALL exchange_horiz( aerosol_mass(icc)%conc_p, nbgp,                            &
+                                        alternative_communicator = communicator_salsa )
                 ENDDO
              ENDDO
              IF ( .NOT. salsa_gases_from_chem )  THEN
                 DO  ig = 1, ngases_salsa
-                   CALL exchange_horiz( salsa_gas(ig)%conc_p, nbgp )
+                   CALL exchange_horiz( salsa_gas(ig)%conc_p, nbgp,                                &
+                                        alternative_communicator = communicator_salsa )
                 ENDDO
              ENDIF
           ENDIF
+
           CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
 
@@ -845,4 +854,5 @@
 !--       boundary conditions for module-specific variables
           CALL module_interface_boundary_conditions
+
 !
 !--       Incrementing timestep counter
@@ -906,7 +916,8 @@
                    ENDIF
                    IF ( bulk_cloud_model  .AND.  microphysics_ice_extension )  THEN
-                       CALL exchange_horiz( qi, nbgp )
-                       CALL exchange_horiz( ni, nbgp )
+                      CALL exchange_horiz( qi, nbgp )
+                      CALL exchange_horiz( ni, nbgp )
                    ENDIF
+
                 ENDIF
 
@@ -927,13 +938,16 @@
                 IF ( salsa  .AND. time_since_reference_point >= skip_time_do_salsa )  THEN
                    DO  ib = 1, nbins_aerosol
-                      CALL exchange_horiz( aerosol_number(ib)%conc, nbgp )
+                      CALL exchange_horiz( aerosol_number(ib)%conc, nbgp,                          &
+                                           alternative_communicator = communicator_salsa )
                       DO  ic = 1, ncomponents_mass
                          icc = ( ic - 1 ) * nbins_aerosol + ib
-                         CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp )
+                         CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp,                        &
+                                              alternative_communicator = communicator_salsa )
                       ENDDO
                    ENDDO
                    IF ( .NOT. salsa_gases_from_chem )  THEN
                       DO  ig = 1, ngases_salsa
-                         CALL exchange_horiz( salsa_gas(ig)%conc, nbgp )
+                         CALL exchange_horiz( salsa_gas(ig)%conc, nbgp,                            &
+                                              alternative_communicator = communicator_salsa )
                       ENDDO
                    ENDIF
@@ -960,21 +974,4 @@
                    ENDDO
                 ENDDO
-             ENDIF
-
-!
-!--          Set SALSA boundary conditions (decycling)
-             IF ( salsa  .AND. time_since_reference_point >= skip_time_do_salsa )  THEN
-                DO  ib = 1, nbins_aerosol
-                   CALL salsa_boundary_conds( aerosol_number(ib)%conc, aerosol_number(ib)%init )
-                   DO  ic = 1, ncomponents_mass
-                      icc = ( ic - 1 ) * nbins_aerosol + ib
-                      CALL salsa_boundary_conds( aerosol_mass(icc)%conc, aerosol_mass(icc)%init )
-                   ENDDO
-                ENDDO
-                IF ( .NOT. salsa_gases_from_chem )  THEN
-                   DO  ig = 1, ngases_salsa
-                      CALL salsa_boundary_conds( salsa_gas(ig)%conc, salsa_gas(ig)%init )
-                   ENDDO
-                ENDIF
              ENDIF
 
