anchor.diff:          real,dimension(1:nx,1:nz,1:ny-1),intent(inout) :: v
           real,dimension(1:nx,1:nz,1:ny),  intent(inout) :: w
           integer, intent(in) :: nx,nz,ny
  -
  +        
           integer :: i,j,k,n
           real, dimension(:), pointer :: q1, U2, l, f
           ! q1 = q after applying previous iteration advection
  
positive.code:                     l=q(:,j,i-1)
                      ! In MPDATA papers A = (r-l)/(r+l)
                      ! compute the denomenator first so we can check that it is not zero
                      denom=(r + l)
                      where(denom==0) denom=1e-10
                      ! U2 is the diffusive pseudo-velocity
                      v2(:,j,i-1) = abs(v(:,j,i-1)) - v(:,j,i-1)**2
                      v2(:,j,i-1) = v2(:,j,i-1) * (r-l) / denom
  
  
                      ! -----------------------
                      ! compute the w component
                      ! -----------------------
                      if (i==ny) then
                          w2(:,j,i)=0
                      else
                          if (j<nz) then
                              r=q(:,j+1,i)
                              l=q(:,j,i)
                              ! In MPDATA papers A = (r-l)/(r+l)
                              ! compute the denomenator first so we can check that it is not zero
                              denom=(r + l)
                              where(denom==0) denom=1e-10
                              ! U2 is the diffusive pseudo-velocity
                              w2(:,j,i) = abs(w(:,j,i)) - w(:,j,i)**2
                              w2(:,j,i) = w2(:,j,i) * (r-l) / denom
                          else
                              w2(:,j,i) = 0
                          endif
                      endif
  
                  endif
              end do
          end do
          !$omp end do
          !$omp end parallel
  
      end subroutine mpdata_fluxes
  
      subroutine flux_limiter(q, q2, u,v,w, nx,nz,ny)
  
negative.code:         real, dimension(nz),   target :: q1z,lz
          real, dimension(nz-1), target :: fz, U2z
          logical :: flux_is_w
          
          real :: qmax_i,qmin_i,qmax_i2,qmin_i2
          real :: beta_in_i, beta_out_i, beta_in_i2, beta_out_i2
          real :: fin_i, fout_i, fin_i2, fout_i2
          
          ! NOTE: before inclusion of FCT_core the following variables must be setup: 
          ! q1 and l (l=q0)
          !$omp parallel shared(q2,q,u,v,w) firstprivate(nx,ny,nz) default(private)
          !$omp do schedule(static)
  