Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 4078)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 4079)
@@ -25,4 +25,11 @@
 ! -----------------
 ! $Id$
+! Implementation of a flux limiter according to Skamarock (2006) for the 
+! vertical scalar advection. Please note, this is only implemented for the 
+! cache-optimized version at the moment. Implementation for the vector-
+! optimized version will follow after critical issues concerning 
+! vectorization are fixed.
+! 
+! 3873 2019-04-08 15:44:30Z knoop
 ! Moved ocean_mode specific code to ocean_mod
 ! 
@@ -273,4 +280,6 @@
 !> the divergence is not sufficiently reduced, resulting in erroneous fluxes 
 !> and could lead to numerical instabilities. 
+!>
+!> @todo Implement monotonic flux limiter also for vector version.
 !------------------------------------------------------------------------------!
  MODULE advec_ws
@@ -296,5 +305,5 @@
                bc_dirichlet_s, bc_radiation_l, bc_radiation_n,                 &
                bc_radiation_r, bc_radiation_s, intermediate_timestep_count,    &
-               u_gtrans, v_gtrans
+               u_gtrans, v_gtrans, dt_3d
 
     USE indices,                                                               &
@@ -1116,5 +1125,5 @@
     SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char, swap_flux_y_local,            &
                               swap_diss_y_local, swap_flux_x_local,            &
-                              swap_diss_x_local, i_omp, tn )
+                              swap_diss_x_local, i_omp, tn, flux_limitation )
 
 
@@ -1126,23 +1135,39 @@
        INTEGER(iwp) ::  k         !< grid index along z-direction
        INTEGER(iwp) ::  k_mm      !< k-2 index in disretization, can be modified to avoid segmentation faults
+       INTEGER(iwp) ::  k_mmm     !< k-3 index in disretization, can be modified to avoid segmentation faults
        INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
        INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
        INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms  
        INTEGER(iwp) ::  tn        !< number of OpenMP thread
-       
-       REAL(wp)     ::  ibit0  !< flag indicating 1st-order scheme along x-direction
-       REAL(wp)     ::  ibit1  !< flag indicating 3rd-order scheme along x-direction
-       REAL(wp)     ::  ibit2  !< flag indicating 5th-order scheme along x-direction
-       REAL(wp)     ::  ibit3  !< flag indicating 1st-order scheme along y-direction
-       REAL(wp)     ::  ibit4  !< flag indicating 3rd-order scheme along y-direction
-       REAL(wp)     ::  ibit5  !< flag indicating 5th-order scheme along y-direction
-       REAL(wp)     ::  ibit6  !< flag indicating 1st-order scheme along z-direction
-       REAL(wp)     ::  ibit7  !< flag indicating 3rd-order scheme along z-direction
-       REAL(wp)     ::  ibit8  !< flag indicating 5th-order scheme along z-direction
-       REAL(wp)     ::  diss_d !< artificial dissipation term at grid box bottom
-       REAL(wp)     ::  div    !< diverence on scalar grid
-       REAL(wp)     ::  flux_d !< 6th-order flux at grid box bottom
-       REAL(wp)     ::  u_comp !< advection velocity along x-direction 
-       REAL(wp)     ::  v_comp !< advection velocity along y-direction 
+
+       LOGICAL, OPTIONAL ::  flux_limitation !< flag indicating flux limitation of the vertical advection
+       LOGICAL           ::  limiter         !< control flag indicating the application of flux limitation
+
+       REAL(wp) ::  diss_d        !< artificial dissipation term at grid box bottom
+       REAL(wp) ::  div           !< velocity diverence on scalar grid
+       REAL(wp) ::  div_in        !< vertical flux divergence of ingoing fluxes
+       REAL(wp) ::  div_out       !< vertical flux divergence of outgoing fluxes      
+       REAL(wp) ::  f_corr_t      !< correction flux at grid-cell top, i.e. the difference between high and low-order flux
+       REAL(wp) ::  f_corr_d      !< correction flux at grid-cell bottom, i.e. the difference between high and low-order flux
+       REAL(wp) ::  f_corr_t_in   !< correction flux of ingoing flux part at grid-cell top
+       REAL(wp) ::  f_corr_d_in   !< correction flux of ingoing flux part at grid-cell bottom
+       REAL(wp) ::  f_corr_t_out  !< correction flux of outgoing flux part at grid-cell top
+       REAL(wp) ::  f_corr_d_out  !< correction flux of outgoing flux part at grid-cell bottom
+       REAL(wp) ::  fac_correction!< factor to limit the in- and outgoing fluxes
+       REAL(wp) ::  flux_d        !< 6th-order flux at grid box bottom
+       REAL(wp) ::  ibit0         !< flag indicating 1st-order scheme along x-direction
+       REAL(wp) ::  ibit1         !< flag indicating 3rd-order scheme along x-direction
+       REAL(wp) ::  ibit2         !< flag indicating 5th-order scheme along x-direction
+       REAL(wp) ::  ibit3         !< flag indicating 1st-order scheme along y-direction
+       REAL(wp) ::  ibit4         !< flag indicating 3rd-order scheme along y-direction
+       REAL(wp) ::  ibit5         !< flag indicating 5th-order scheme along y-direction
+       REAL(wp) ::  ibit6         !< flag indicating 1st-order scheme along z-direction
+       REAL(wp) ::  ibit7         !< flag indicating 3rd-order scheme along z-direction
+       REAL(wp) ::  ibit8         !< flag indicating 5th-order scheme along z-direction
+       REAL(wp) ::  max_val       !< maximum value of the quanitity along the numerical stencil (in vertical direction)
+       REAL(wp) ::  min_val       !< maximum value of the quanitity along the numerical stencil (in vertical direction)
+       REAL(wp) ::  mon           !< monotone solution of the advection equation using 1st-order fluxes
+       REAL(wp) ::  u_comp        !< advection velocity along x-direction 
+       REAL(wp) ::  v_comp        !< advection velocity along y-direction 
 !       
 !--    sk is an array from parameter list. It should not be a pointer, because 
@@ -1153,10 +1178,11 @@
        REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !<  advected scalar
 
-       REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_r !< discretized artificial dissipation at rightward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_t !< discretized artificial dissipation at top of the grid box
-       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_n !< discretized 6th-order flux at northward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_r !< discretized 6th-order flux at rightward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_t !< discretized 6th-order flux at top of the grid box
+       REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_n     !< discretized artificial dissipation at northward-side of the grid box
+       REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_r     !< discretized artificial dissipation at rightward-side of the grid box
+       REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_t     !< discretized artificial dissipation at top of the grid box
+       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_n     !< discretized 6th-order flux at northward-side of the grid box
+       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_r     !< discretized 6th-order flux at rightward-side of the grid box
+       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_t     !< discretized 6th-order flux at top of the grid box
+       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_t_1st !< discretized 1st-order flux at top of the grid box
        
        REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) ::  swap_diss_y_local !< discretized artificial dissipation at southward-side of the grid box 
@@ -1178,4 +1204,8 @@
           nzb_max_l = nzb_max
        END IF
+!
+!--    Set control flag for flux limiter 
+       limiter = .FALSE.
+       IF ( PRESENT( flux_limitation) )  limiter = flux_limitation
 !
 !--    Compute southside fluxes of the respective PE bounds.
@@ -1231,9 +1261,9 @@
                                   +           ( sk(k,j+2,i) + sk(k,j-3,i) )   &
                                                 ) * adv_sca_5
-              swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                    &
+             swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                     &
                                     10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )   &
                                   -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )   &
                                   +             sk(k,j+2,i) - sk(k,j-3,i)     &
-                                                         ) * adv_sca_5
+                                                        ) * adv_sca_5
 
           ENDDO
@@ -1266,5 +1296,5 @@
                                                   )
 
-              swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                  &
+             swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                   &
                                                ( 10.0_wp * ibit2 * adv_sca_5  &
                                             +     3.0_wp * ibit1 * adv_sca_3  &
@@ -1279,5 +1309,5 @@
                                                ) *                            &
                                             ( sk(k,j,i+2) - sk(k,j,i-3)    )  &
-                                                           )
+                                                          )
 
           ENDDO
@@ -1535,5 +1565,108 @@
                              ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
                                                          )
-       ENDDO 
+       ENDDO
+
+       IF ( limiter )  THEN
+!
+!--       Compute monotone first-order fluxes which are required for mononte
+!--       flux limitation.
+          flux_t_1st(nzb) = 0.0_wp
+          DO  k = nzb+1, nzb_max_l
+             flux_t_1st(k) = ( w(k,j,i)   * ( sk(k+1,j,i)  + sk(k,j,i) )       &
+                       -  ABS( w(k,j,i) ) * ( sk(k+1,j,i)  - sk(k,j,i) ) )     &
+                           * rho_air_zw(k) * adv_sca_1
+!
+!--          In flux limitation the total flux will be corrected. For the sake
+!--          of cleariness the higher-order advective and disspative fluxes 
+!--          will be merged onto flux_t. 
+             flux_t(k) = flux_t(k) + diss_t(k)
+             diss_t(k) = 0.0_wp
+          ENDDO
+!
+!--       Flux limitation of vertical fluxes according to Skamarock (2006). 
+!--       Please note, as flux limitation implies linear dependencies of fluxes,
+!--       flux limitation is only made for the vertical advection term. Limitation
+!--       of the horizontal terms cannot be parallelized. 
+!--       Due to the linear dependency, the following loop will not be vectorized.
+!--       Further, note that the flux limiter is only applied within the urban
+!--       layer, i.e up to the topography top.  
+          DO  k = nzb+1, nzb_max_l
+!
+!--          Compute one-dimensional divergence along the vertical direction,
+!--          which is used to correct the advection discretization. This is 
+!--          necessary as in one-dimensional space the advection velocity
+!--          should actually be constant. 
+             div = ( w(k,j,i)   * rho_air_zw(k)                                &
+                   - w(k-1,j,i) * rho_air_zw(k-1)                              &      
+                   ) * drho_air(k) * ddzw(k)
+!
+!--          Compute monotone solution of the advection equation from 
+!--          1st-order fluxes. Please note, the advection equation is corrected
+!--          by the divergence term (in 1D the advective flow should be divergence
+!--          free). Moreover, please note, as time-increment the full timestep 
+!--          is used, even though a Runge-Kutta scheme will be used. However, 
+!--          the length of the actual time increment is not important at all 
+!--          since it cancels out later when the fluxes are limited.   
+             mon = sk(k,j,i) + ( - ( flux_t_1st(k) - flux_t_1st(k-1) )         &
+                             * drho_air(k) * ddzw(k)                           &
+                             + div * sk(k,j,i)                                 &
+                               ) * dt_3d 
+!
+!--          Determine minimum and maximum values along the numerical stencil.
+             k_mmm = MAX( k - 3, nzb + 1 )
+             k_ppp = MIN( k + 3, nzt + 1 ) 
+
+             min_val = MINVAL( sk(k_mmm:k_ppp,j,i) )
+             max_val = MAXVAL( sk(k_mmm:k_ppp,j,i) )
+!
+!--          Compute difference between high- and low-order fluxes, which may 
+!--          act as correction fluxes
+             f_corr_t = flux_t(k)   - flux_t_1st(k)
+             f_corr_d = flux_t(k-1) - flux_t_1st(k-1)
+!
+!--          Determine outgoing fluxes, i.e. the part of the fluxes which can 
+!--          decrease the value within the grid box
+             f_corr_t_out = MAX( 0.0_wp, f_corr_t )
+             f_corr_d_out = MIN( 0.0_wp, f_corr_d )
+!
+!--          Determine ingoing fluxes, i.e. the part of the fluxes which can 
+!--          increase the value within the grid box 
+             f_corr_t_in = MIN( 0.0_wp, f_corr_t)
+             f_corr_d_in = MAX( 0.0_wp, f_corr_d)
+!
+!--          Compute divergence of outgoing correction fluxes
+             div_out = - ( f_corr_t_out - f_corr_d_out ) * drho_air(k)         &
+                                                         * ddzw(k) * dt_3d
+!
+!--          Compute divergence of ingoing correction fluxes
+             div_in = - ( f_corr_t_in - f_corr_d_in )    * drho_air(k)         &
+                                                         * ddzw(k) * dt_3d
+!
+!--          Check if outgoing fluxes can lead to undershoots, i.e. values smaller
+!--          than the minimum value within the numerical stencil. If so, limit
+!--          them. 
+             IF ( mon - min_val < - div_out  .AND.  ABS( div_out ) > 0.0_wp )  &
+             THEN
+                fac_correction = ( mon - min_val ) / ( - div_out )
+                f_corr_t_out = f_corr_t_out * fac_correction
+                f_corr_d_out = f_corr_d_out * fac_correction
+             ENDIF
+!
+!--          Check if ingoing fluxes can lead to overshoots, i.e. values larger
+!--          than the maximum value within the numerical stencil. If so, limit
+!--          them. 
+             IF ( mon - max_val > - div_in  .AND.  ABS( div_in ) > 0.0_wp )    &
+             THEN
+                fac_correction = ( mon - max_val ) / ( - div_in )
+                f_corr_t_in = f_corr_t_in * fac_correction
+                f_corr_d_in = f_corr_d_in * fac_correction
+             ENDIF
+!
+!--          Finally add the limited fluxes to the original ones. If no 
+!--          flux limitation was done, the fluxes equal the original ones. 
+             flux_t(k)   = flux_t_1st(k)   + f_corr_t_out + f_corr_t_in
+             flux_t(k-1) = flux_t_1st(k-1) + f_corr_d_out + f_corr_d_in
+          ENDDO
+       ENDIF
 
        DO  k = nzb+1, nzb_max_l
Index: palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 4078)
+++ palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 4079)
@@ -27,4 +27,9 @@
 ! -----------------
 ! $Id$
+! Application of monotonic flux limiter for the vertical scalar advection 
+! up to the topography top (only for the cache-optimized version at the 
+! moment).
+! 
+! 4069 2019-07-01 14:05:51Z Giersch
 ! Masked output running index mid has been introduced as a local variable to 
 ! avoid runtime error (Loop variable has been modified) in time_integration
@@ -354,5 +359,7 @@
                 dt_3d, humidity, initializing_actions, message_string,                             &
                 omega, tsc, intermediate_timestep_count, intermediate_timestep_count_max,          &
-         max_pr_user, timestep_scheme, use_prescribed_profile_data, ws_scheme_sca, air_chemistry
+                max_pr_user,                                                                       &
+                monotonic_limiter_z,                                                               &
+                timestep_scheme, use_prescribed_profile_data, ws_scheme_sca, air_chemistry
 
     USE arrays_3d,                                                                                 &
@@ -2750,5 +2757,5 @@
              CALL advec_s_ws( i, j, chem_species(ilsp)%conc, 'kc', chem_species(ilsp)%flux_s_cs,   &
                               chem_species(ilsp)%diss_s_cs, chem_species(ilsp)%flux_l_cs,          &
-                              chem_species(ilsp)%diss_l_cs, i_omp_start, tn )
+                              chem_species(ilsp)%diss_l_cs, i_omp_start, tn, monotonic_limiter_z )
           ELSE
              CALL advec_s_pw( i, j, chem_species(ilsp)%conc )
@@ -2782,8 +2789,4 @@
                BTEST( wall_flags_0(k,j,i), 0 )                                                     &
                )
-
-          IF ( chem_species(ilsp)%conc_p(k,j,i) < 0.0_wp )  THEN
-             chem_species(ilsp)%conc_p(k,j,i) = 0.1_wp * chem_species(ilsp)%conc(k,j,i)    !FKS6
-          ENDIF
        ENDDO
 !
@@ -3315,4 +3318,5 @@
        r_aero_surf = surf_lsm_h%r_a(m)
        solar_rad   = surf_lsm_h%rad_sw_dir(m) + surf_lsm_h%rad_sw_dif(m)
+       
        lai = surf_lsm_h%lai(m)
        sai = lai + 1
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 4078)
+++ palm/trunk/SOURCE/modules.f90	(revision 4079)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! + monotonic_limiter_z
+! 
+! 4069 2019-07-01 14:05:51Z Giersch
 ! Masked output running index mid has been introduced as a local variable to 
 ! avoid runtime error (Loop variable has been modified) in time_integration
@@ -1360,4 +1363,5 @@
     LOGICAL ::  masking_method = .FALSE.                         !< namelist parameter
     LOGICAL ::  mg_switch_to_pe0 = .FALSE.                       !< internal multigrid switch for steering the ghost point exchange in case that data has been collected on PE0
+    LOGICAL ::  monotonic_limiter_z = .FALSE.                    !< use monotonic flux limiter for vertical scalar advection
     LOGICAL ::  nesting_offline = .FALSE.                        !< flag controlling offline nesting in COSMO model  
     LOGICAL ::  neutral = .FALSE.                                !< namelist parameter
Index: palm/trunk/SOURCE/nesting_offl_mod.f90
===================================================================
--- palm/trunk/SOURCE/nesting_offl_mod.f90	(revision 4078)
+++ palm/trunk/SOURCE/nesting_offl_mod.f90	(revision 4079)
@@ -25,4 +25,12 @@
 ! -----------------
 ! $Id$
+! - Set boundary condition for w at nzt+1 at the lateral boundaries, even 
+!   though these won't enter the numerical solution. However, due to the mass
+!   conservation these values might some up to very large values which will
+!   occur in the run-control file
+! - Bugfix in offline nesting of chemical species
+! - Do not set Neumann conditions for TKE and passive scalar
+! 
+! 4022 2019-06-12 11:52:39Z suehring
 ! Detection of boundary-layer depth in stable boundary layer on basis of 
 ! boundary data improved 
@@ -353,4 +361,5 @@
                                    BTEST( wall_flags_0(k,j,-1), 3 ) )
              ENDDO
+             w(nzt,j,-1) = w(nzt-1,j,-1)
           ENDDO
 
@@ -427,4 +436,5 @@
                                     BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
              ENDDO
+             w(nzt,j,nxr+1) = w(nzt-1,j,nxr+1)
           ENDDO
 
@@ -504,4 +514,5 @@
                                   BTEST( wall_flags_0(k,-1,i), 3 ) )
              ENDDO
+             w(nzt,-1,i) = w(nzt-1,-1,i)
           ENDDO
 
@@ -580,4 +591,5 @@
                                       BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
              ENDDO
+             w(nzt,nyn+1,i) = w(nzt-1,nyn+1,i)
           ENDDO
 
@@ -710,6 +722,6 @@
                    DO  j = nys, nyn
                       chem_species(n)%conc(nzt+1,j,i) = interpolate_in_time(   &
-                                              nest_offl%chem_north(0,j,i,n),   &
-                                              nest_offl%chem_north(1,j,i,n),   &
+                                              nest_offl%chem_top(0,j,i,n),   &
+                                              nest_offl%chem_top(1,j,i,n),   &
                                               fac_dt                       )
                    ENDDO
@@ -727,17 +739,17 @@
           IF (  bc_dirichlet_n )  diss(:,nyn+1,:) = diss(:,nyn,:)
        ENDIF
-       IF ( .NOT. constant_diffusion )  THEN
-          IF (  bc_dirichlet_l )  e(:,:,nxl-1) = e(:,:,nxl)
-          IF (  bc_dirichlet_r )  e(:,:,nxr+1) = e(:,:,nxr)
-          IF (  bc_dirichlet_s )  e(:,nys-1,:) = e(:,nys,:)
-          IF (  bc_dirichlet_n )  e(:,nyn+1,:) = e(:,nyn,:)
-          e(nzt+1,:,:) = e(nzt,:,:)
-       ENDIF
-       IF ( passive_scalar )  THEN
-          IF (  bc_dirichlet_l )  s(:,:,nxl-1) = s(:,:,nxl)
-          IF (  bc_dirichlet_r )  s(:,:,nxr+1) = s(:,:,nxr)
-          IF (  bc_dirichlet_s )  s(:,nys-1,:) = s(:,nys,:)
-          IF (  bc_dirichlet_n )  s(:,nyn+1,:) = s(:,nyn,:)
-       ENDIF
+!        IF ( .NOT. constant_diffusion )  THEN
+!           IF (  bc_dirichlet_l )  e(:,:,nxl-1) = e(:,:,nxl)
+!           IF (  bc_dirichlet_r )  e(:,:,nxr+1) = e(:,:,nxr)
+!           IF (  bc_dirichlet_s )  e(:,nys-1,:) = e(:,nys,:)
+!           IF (  bc_dirichlet_n )  e(:,nyn+1,:) = e(:,nyn,:)
+!           e(nzt+1,:,:) = e(nzt,:,:)
+!        ENDIF
+!        IF ( passive_scalar )  THEN
+!           IF (  bc_dirichlet_l )  s(:,:,nxl-1) = s(:,:,nxl)
+!           IF (  bc_dirichlet_r )  s(:,:,nxr+1) = s(:,:,nxr)
+!           IF (  bc_dirichlet_s )  s(:,nys-1,:) = s(:,nys,:)
+!           IF (  bc_dirichlet_n )  s(:,nyn+1,:) = s(:,nyn,:)
+!        ENDIF
 
        CALL exchange_horiz( u, nbgp )
@@ -755,5 +767,8 @@
           ENDDO
        ENDIF
-       
+!
+!--    Set top boundary condition at all horizontal grid points, also at the
+!--    lateral boundary grid points. 
+       w(nzt+1,:,:) = w(nzt,:,:)       
 !
 !--    In case of Rayleigh damping, where the profiles u_init, v_init
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 4078)
+++ palm/trunk/SOURCE/parin.f90	(revision 4079)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! +monotonic_limiter_z
+! 
+! 4022 2019-06-12 11:52:39Z suehring
 ! Change default top boundary condition for pressure to Neumann in offline
 ! nesting case
@@ -590,4 +593,5 @@
              loop_optimization, lsf_exception, masking_method, mg_cycles,      &
              mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
+             monotonic_limiter_z,                                              &
              netcdf_precision, neutral, ngsrb,                                 &
              nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega,           &
@@ -660,4 +664,5 @@
              loop_optimization, lsf_exception, masking_method, mg_cycles,      &
              mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
+             monotonic_limiter_z,                                              &
              netcdf_precision, neutral, ngsrb,                                 &
              nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega,           &
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 4078)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 4079)
@@ -25,4 +25,10 @@
 ! -----------------
 ! $Id$
+! Application of monotonic flux limiter for the vertical scalar advection 
+! up to the topography top (only for the cache-optimized version at the 
+! moment). Please note, at the moment the limiter is only applied for passive
+! scalars.
+! 
+! 4048 2019-06-21 21:00:21Z knoop
 ! Moved tcm_prognostic_equations to module_interface
 ! 
@@ -419,5 +425,7 @@
                humidity, intermediate_timestep_count,                          &
                intermediate_timestep_count_max, large_scale_forcing,           &
-               large_scale_subsidence, neutral, nudging,                       &
+               large_scale_subsidence,                                         &
+               monotonic_limiter_z,                                            &
+               neutral, nudging,                                               &
                ocean_mode, passive_scalar, plant_canopy, pt_reference,         &
                scalar_advec, scalar_advec, simulated_time, sloping_surface,    &
@@ -978,6 +986,10 @@
              THEN
                 IF ( ws_scheme_sca )  THEN
+!
+!--                For scalar advection apply monotonic flux limiter near 
+!--                topography.
                    CALL advec_s_ws( i, j, s, 's', flux_s_s, &
-                                diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn )
+                                    diss_s_s, flux_l_s, diss_l_s, i_omp_start, &
+                                    tn, monotonic_limiter_z )
                 ELSE
                    CALL advec_s_pw( i, j, s )
Index: palm/trunk/SOURCE/salsa_mod.f90
===================================================================
--- palm/trunk/SOURCE/salsa_mod.f90	(revision 4078)
+++ palm/trunk/SOURCE/salsa_mod.f90	(revision 4079)
@@ -26,4 +26,9 @@
 ! -----------------
 ! $Id$
+! Application of monotonic flux limiter for the vertical scalar advection 
+! up to the topography top (only for the cache-optimized version at the 
+! moment).
+! 
+! 4069 2019-07-01 14:05:51Z Giersch
 ! Masked output running index mid has been introduced as a local variable to 
 ! avoid runtime error (Loop variable has been modified) in time_integration
@@ -7773,5 +7778,6 @@
     IF ( timestep_scheme(1:5) == 'runge' )  THEN
        IF ( ws_scheme_sca )  THEN
-          CALL advec_s_ws( i, j, rs, id, flux_s, diss_s, flux_l, diss_l, i_omp_start, tn )
+          CALL advec_s_ws( i, j, rs, id, flux_s, diss_s, flux_l, diss_l, i_omp_start, tn,          &
+                           monotonic_limiter_z )
        ELSE
           CALL advec_s_pw( i, j, rs )
