Index: palm/trunk/SOURCE/salsa_mod.f90
===================================================================
--- palm/trunk/SOURCE/salsa_mod.f90	(revision 4116)
+++ palm/trunk/SOURCE/salsa_mod.f90	(revision 4117)
@@ -21,5 +21,17 @@
 ! Current revisions:
 ! -----------------
-! 
+! - When Dirichlet condition is applied in decycling, the boundary conditions are
+!   only set at the ghost points and not at the prognostic grid points as done
+!   before
+! - Rename decycle_ns/lr to decycle_salsa_ns/lr and decycle_method to
+!   decycle_method_salsa
+! - Allocation and initialization of special advection flags salsa_advc_flags_s
+!   used for salsa. These are exclusively used for salsa variables to
+!   distinguish from the usually-used flags which might be different when
+!   decycling is applied in combination with cyclic boundary conditions.
+!   Moreover, salsa_advc_flags_s considers extended zones around buildings where
+!   the first-order upwind scheme is applied for the horizontal advection terms.
+!   This is doneto overcome high concentration peaks due to stationary numerical
+!   oscillations caused by horizontal advection discretization.
 ! 
 ! Former revisions:
@@ -192,5 +204,11 @@
         ONLY:  call_chem_at_all_substeps, chem_gasphase_on, chem_species
 
-    USE control_parameters
+    USE control_parameters,                                                                        &
+        ONLY:  air_chemistry, bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s,      &
+               bc_lr, bc_lr_cyc, bc_ns, bc_ns_cyc, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
+               bc_radiation_s, coupling_char, debug_output, dt_3d, intermediate_timestep_count,    &
+               intermediate_timestep_count_max, land_surface, message_string, monotonic_limiter_z, &
+               plant_canopy, pt_surface, salsa, scalar_advec, surface_pressure,                    &
+               time_since_reference_point, timestep_scheme, tsc, urban_surface, ws_scheme_sca
 
     USE indices,                                                                                   &
@@ -337,5 +355,5 @@
                                                                   !< 'parameterized', 'read_from_file'
 
-    CHARACTER(LEN=20), DIMENSION(4) ::  decycle_method =                                           &
+    CHARACTER(LEN=20), DIMENSION(4) ::  decycle_method_salsa =                                     &
                                                  (/'dirichlet','dirichlet','dirichlet','dirichlet'/)
                                      !< Decycling method at horizontal boundaries
@@ -404,9 +422,13 @@
 
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  k_topo_top  !< vertical index of the topography top
+
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE  ::  salsa_advc_flags_s !< flags used to degrade order of advection
+                                                                        !< scheme for salsa variables near walls and
+                                                                        !< lateral boundaries
 !
 !-- SALSA switches:
     LOGICAL ::  advect_particle_water   = .TRUE.   !< Advect water concentration of particles
-    LOGICAL ::  decycle_lr              = .FALSE.  !< Undo cyclic boundaries: left and right
-    LOGICAL ::  decycle_ns              = .FALSE.  !< Undo cyclic boundaries: north and south
+    LOGICAL ::  decycle_salsa_lr        = .FALSE.  !< Undo cyclic boundaries: left and right
+    LOGICAL ::  decycle_salsa_ns        = .FALSE.  !< Undo cyclic boundaries: north and south
     LOGICAL ::  include_emission        = .FALSE.  !< Include or not emissions
     LOGICAL ::  feedback_to_palm        = .FALSE.  !< Allow feedback due to condensation of H2O
@@ -870,6 +892,7 @@
     NAMELIST /salsa_parameters/      aerosol_flux_dpg, aerosol_flux_mass_fracs_a,                  &
                                      aerosol_flux_mass_fracs_b, aerosol_flux_sigmag,               &
-                                     advect_particle_water, bc_salsa_b, bc_salsa_t, decycle_lr,    &
-                                     decycle_method, decycle_ns, depo_pcm_par, depo_pcm_type,      &
+                                     advect_particle_water, bc_salsa_b, bc_salsa_t,                &
+                                     decycle_salsa_lr, decycle_method_salsa, decycle_salsa_ns,     &
+                                     depo_pcm_par, depo_pcm_type,                                  &
                                      depo_surf_par, dpg, dt_salsa, feedback_to_palm, h2so4_init,   &
                                      hno3_init, init_gases_type, init_aerosol_type, listspec,      &
@@ -909,5 +932,5 @@
 
     USE control_parameters,                                                                        &
-        ONLY:  message_string
+        ONLY:  humidity
 
     IMPLICIT NONE
@@ -1127,6 +1150,6 @@
  SUBROUTINE salsa_init_arrays
 
-    USE chem_gasphase_mod,                                                                         &
-        ONLY:  nvar
+    USE advec_ws,                                                                                  &
+        ONLY: ws_init_flags_scalar
 
     USE surface_mod,                                                                               &
@@ -1368,4 +1391,82 @@
 
     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
+!-- 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
+!-- order of the scalar advection scheme near the lateral boundaries for passive scalars with
+!-- 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_0 bit 31 instead of
+!--    cyclic boundary conditions. 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_0 bit 31 must not be
+!--    modified. Hence, store the boundary conditions directly on salsa_advc_flags_s.
+!--    salsa_advc_flags_s will be later overwritten in ws_init_flags_scalar and bit 31 won't be used
+!--    to control the numerical order.
+!--    Initialize with flag 31 only.
+       salsa_advc_flags_s = 0
+       salsa_advc_flags_s = MERGE( IBSET( salsa_advc_flags_s, 31 ), 0, BTEST( wall_flags_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
+!--    ws_init_flags_scalar. The last argument in ws_init_flags_scalar indicates that a passive
+!--    scalar is being treated and the horizontal advection terms are degraded already 2 grid points
+!--    before the lateral boundary. Also, extended degradation zones are applied, where
+!--    horizontal advection of scalars is discretised by the first-order scheme at all grid points
+!--    in the vicinity of buildings (<= 3 grid points). Even though no building is within the
+!--    numerical stencil, the first-order scheme is used. At fourth and fifth grid points, the order
+!--    of the horizontal advection scheme is successively upgraded.
+!--    These degradations of the advection scheme are done to avoid stationary numerical
+!--    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,    &
+                                  bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  decycle_salsa_ns,    &
+                                  bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  decycle_salsa_lr,    &
+                                  bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  decycle_salsa_ns,    &
+                                  salsa_advc_flags_s, .TRUE. )
+    ENDIF
+
 
  END SUBROUTINE salsa_init_arrays
@@ -2531,4 +2632,7 @@
                              nyn_on_file, nysf, nysc, nys_on_file, tmp_3d, found )
 
+    USE control_parameters,                                                                        &
+        ONLY:  length, restart_string
+
     IMPLICIT NONE
 
@@ -2602,4 +2706,7 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE salsa_wrd_local
+
+    USE control_parameters,                                                                        &
+        ONLY:  write_binary
 
     IMPLICIT NONE
@@ -3212,7 +3319,4 @@
     USE basic_constants_and_equations_mod,                                                         &
         ONLY:  barometric_formula, exner_function, ideal_gas_law_rho, magnus
-
-    USE control_parameters,                                                                        &
-        ONLY: pt_surface, surface_pressure
 
     IMPLICIT NONE
@@ -3509,6 +3613,6 @@
  SUBROUTINE deposition( paero, tk, adn, mag_u, lad, kvis, schmidt_num, vc )
 
-    USE plant_canopy_model_mod,                                                &
-        ONLY: cdc
+    USE plant_canopy_model_mod,                                                                    &
+        ONLY:  cdc
 
     IMPLICIT NONE
@@ -7604,7 +7708,4 @@
  SUBROUTINE salsa_prognostic_equations_ij( i, j, i_omp_start, tn )
 
-    USE control_parameters,                                                                        &
-        ONLY:  time_since_reference_point
-
     IMPLICIT NONE
 
@@ -7671,7 +7772,4 @@
  SUBROUTINE salsa_prognostic_equations()
 
-    USE control_parameters,                                                                        &
-        ONLY:  time_since_reference_point
-
     IMPLICIT NONE
 
@@ -7743,8 +7841,5 @@
 
     USE indices,                                                                                   &
-        ONLY:  advc_flags_s, wall_flags_0
-
-    USE pegrid,                                                                                    &
-        ONLY:  threads_per_task
+        ONLY:  wall_flags_0
 
     USE surface_mod,                                                                               &
@@ -7786,11 +7881,9 @@
     IF ( timestep_scheme(1:5) == 'runge' )  THEN
        IF ( ws_scheme_sca )  THEN
-          CALL advec_s_ws( 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,                                   &
+          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,                  &
                            bc_dirichlet_n  .OR.  bc_radiation_n,                                   &
                            bc_dirichlet_r  .OR.  bc_radiation_r,                                   &
-                           bc_dirichlet_s  .OR.  bc_radiation_s,                                   &
-                           monotonic_limiter_z )
+                           bc_dirichlet_s  .OR.  bc_radiation_s, monotonic_limiter_z )
        ELSE
           CALL advec_s_pw( i, j, rs )
@@ -7892,5 +7985,5 @@
         ONLY:  diffusion_s
     USE indices,                                                                                   &
-        ONLY:  advc_flags_s, wall_flags_0
+        ONLY:  wall_flags_0
     USE surface_mod,                                                                               &
         ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
@@ -7923,6 +8016,5 @@
     IF ( timestep_scheme(1:5) == 'runge' )  THEN
        IF ( ws_scheme_sca )  THEN
-          CALL advec_s_ws( advc_flags_s, rs, id,                                                   &
-                           bc_dirichlet_l  .OR.  bc_radiation_l,                                   &
+          CALL advec_s_ws( salsa_advc_flags_s, rs, id, bc_dirichlet_l  .OR.  bc_radiation_l,       &
                            bc_dirichlet_n  .OR.  bc_radiation_n,                                   &
                            bc_dirichlet_r  .OR.  bc_radiation_r,                                   &
@@ -8248,9 +8340,9 @@
 !
 !-- Left and right boundaries
-    IF ( decycle_lr  .AND.  ( bc_lr_cyc  .OR. bc_lr == 'nested' ) )  THEN
+    IF ( decycle_salsa_lr  .AND.  ( bc_lr_cyc  .OR. bc_lr == 'nested' ) )  THEN
 
        DO  boundary = 1, 2
 
-          IF ( decycle_method(boundary) == 'dirichlet' )  THEN
+          IF ( decycle_method_salsa(boundary) == 'dirichlet' )  THEN
 !
 !--          Initial profile is copied to ghost and first three layers
@@ -8259,7 +8351,7 @@
              IF ( boundary == 1  .AND.  nxl == 0 )  THEN
                 ss = nxlg
-                ee = nxl+2
+                ee = nxl-1
              ELSEIF ( boundary == 2  .AND.  nxr == nx )  THEN
-                ss = nxr-2
+                ss = nxr+1
                 ee = nxrg
              ENDIF
@@ -8274,5 +8366,5 @@
              ENDDO
 
-          ELSEIF ( decycle_method(boundary) == 'neumann' )  THEN
+          ELSEIF ( decycle_method_salsa(boundary) == 'neumann' )  THEN
 !
 !--          The value at the boundary is copied to the ghost layers to simulate an outlet with
@@ -8300,6 +8392,6 @@
 
           ELSE
-             WRITE(message_string,*) 'unknown decycling method: decycle_method (', boundary,       &
-                                     ') ="' // TRIM( decycle_method(boundary) ) // '"'
+             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 )
           ENDIF
@@ -8309,9 +8401,9 @@
 !
 !-- South and north boundaries
-     IF ( decycle_ns  .AND.  ( bc_ns_cyc  .OR. bc_ns == 'nested' ) )  THEN
+     IF ( decycle_salsa_ns  .AND.  ( bc_ns_cyc  .OR. bc_ns == 'nested' ) )  THEN
 
        DO  boundary = 3, 4
 
-          IF ( decycle_method(boundary) == 'dirichlet' )  THEN
+          IF ( decycle_method_salsa(boundary) == 'dirichlet' )  THEN
 !
 !--          Initial profile is copied to ghost and first three layers
@@ -8320,7 +8412,7 @@
              IF ( boundary == 3  .AND.  nys == 0 )  THEN
                 ss = nysg
-                ee = nys+2
+                ee = nys-1
              ELSEIF ( boundary == 4  .AND.  nyn == ny )  THEN
-                ss = nyn-2
+                ss = nyn+1
                 ee = nyng
              ENDIF
@@ -8335,5 +8427,5 @@
              ENDDO
 
-          ELSEIF ( decycle_method(boundary) == 'neumann' )  THEN
+          ELSEIF ( decycle_method_salsa(boundary) == 'neumann' )  THEN
 !
 !--          The value at the boundary is copied to the ghost layers to simulate an outlet with 
@@ -8361,6 +8453,6 @@
 
           ELSE
-             WRITE(message_string,*) 'unknown decycling method: decycle_method (', boundary,       &
-                                     ') ="' // TRIM( decycle_method(boundary) ) // '"'
+             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 )
           ENDIF
@@ -8417,7 +8509,4 @@
  SUBROUTINE salsa_emission_update
 
-    USE control_parameters,                                                                        &
-        ONLY:  time_since_reference_point
-
     IMPLICIT NONE
 
@@ -8450,7 +8539,4 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE salsa_emission_setup( init )
-
-    USE control_parameters,                                                                        &
-        ONLY:  time_since_reference_point
 
     USE date_and_time_mod,                                                                         &
@@ -9139,7 +9225,4 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE salsa_gas_emission_setup( init )
-
-    USE control_parameters,                                                                        &
-        ONLY:  time_since_reference_point
 
     USE date_and_time_mod,                                                                         &
@@ -9489,7 +9572,4 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE salsa_check_data_output( var, unit )
- 
-    USE control_parameters,                                                                        &
-        ONLY:  message_string
 
     IMPLICIT NONE
@@ -9551,9 +9631,6 @@
  SUBROUTINE salsa_3d_data_averaging( mode, variable )
 
-    USE control_parameters
-
-    USE indices
-
-    USE kinds
+    USE control_parameters,                                                                        &
+        ONLY:  average_count_3d
 
     IMPLICIT NONE
@@ -10746,5 +10823,5 @@
 
     USE control_parameters,                                                                        &
-        ONLY:  mask_size_l, mask_surface
+        ONLY:  mask_i, mask_j, mask_k, mask_size_l, mask_surface, nz_do3d
 
     USE surface_mod,                                                                               &
Index: palm/trunk/TESTS/cases/urban_environment_salsa/INPUT/urban_environment_salsa_p3d
===================================================================
--- palm/trunk/TESTS/cases/urban_environment_salsa/INPUT/urban_environment_salsa_p3d	(revision 4116)
+++ palm/trunk/TESTS/cases/urban_environment_salsa/INPUT/urban_environment_salsa_p3d	(revision 4117)
@@ -20,5 +20,5 @@
 !------------------------------------------------------------------------------
     fft_method           = 'fftw',
-
+    monotonic_limiter_z  = .T., 
 !
 !-- Mode
@@ -243,7 +243,7 @@
 !-- Decycle at the north & south boundaries: initial concentration  
 !--   is copied to ghost and first three layers at both boundaries
-    decycle_lr     = .T.,         
-    decycle_ns     = .F.,
-    decycle_method = 'dirichlet','dirichlet','dirichlet','dirichlet',  
+    decycle_salsa_lr     = .T.,         
+    decycle_salsa_ns     = .F.,
+    decycle_method_salsa = 'dirichlet','dirichlet','dirichlet','dirichlet',  
     bc_salsa_b     = 'neumann', !'dirichlet',  ! surface flux requires 'neumann'
     bc_salsa_t     = 'dirichlet',  ! top
