Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 3873)
+++ 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
