anchor.diff:         endif
   
          if ( thubern_lim ) then
  -          call vertical_limiter_thuburn( q_h(:,:,:),   q_h_pl(:,:,:),   & ! [INOUT]
  -                                         q  (:,:,:),   q_pl  (:,:,:),   & ! [IN]
  -                                         d  (:,:,:),   d_pl  (:,:,:),   & ! [IN]
  -                                         ck (:,:,:,:), ck_pl (:,:,:,:)  ) ! [IN]
  +          call vertical_limiter_thuburn( q_h(:,:,:),   q_h_pl(:,:,:),  & ! [INOUT]
  +                                         q  (:,:,:),   q_pl  (:,:,:),  & ! [IN]
  +                                         d  (:,:,:),   d_pl  (:,:,:),  & ! [IN]
  +                                         ck (:,:,:,:), ck_pl (:,:,:,:) ) ! [IN]
          endif
   
          !--- update rhogq
  
positive.code:           ij     = n
            ip1jp1 = n + 1 + ADM_gall_1d
  
            flux = 0.5D0 * ( (rhovxt(ij,TI)+rhovxt(ij,TJ)) * cinterp_HN(AIJ,1,ij,l) &
                           + (rhovyt(ij,TI)+rhovyt(ij,TJ)) * cinterp_HN(AIJ,2,ij,l) &
                           + (rhovzt(ij,TI)+rhovzt(ij,TJ)) * cinterp_HN(AIJ,3,ij,l) )
  
            flx_h(2,ij    ,k,l) =  flux * cinterp_PRA(ij    ,l) * dt
            flx_h(5,ip1jp1,k,l) = -flux * cinterp_PRA(ip1jp1,l) * dt
  
negative.code:     real(8), intent(in)    :: ck    (ADM_gall,   ADM_kall,ADM_lall   ,2)
      real(8), intent(in)    :: ck_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl,2)
  
      real(8) :: Qout_min   (ADM_gall,   ADM_kall)
      real(8) :: Qout_max   (ADM_gall,   ADM_kall)
      real(8) :: Qout_min_pl(ADM_gall_pl,ADM_kall)
      real(8) :: Qout_max_pl(ADM_gall_pl,ADM_kall)
  
      real(8) :: Qin_minL, Qin_maxL
      real(8) :: Qin_minU, Qin_maxU
      real(8) :: qnext_min, qnext_max
      real(8) :: Cin, Cout
      real(8) :: CQin_min, CQin_max
      real(8) :: inflagL, inflagU
      real(8) :: zerosw
  
      integer :: n, k, l
  