anchor.diff:      real (kind=dp_t) :: beta_xm(nprim), beta_xp(nprim)
   
       real (kind=dp_t) :: dq0, dqp
  -    real (kind=dp_t) :: Iminus(nprim), Iplus(nprim)
  +    real (kind=dp_t) :: Iminus(nprim,nprim), Iplus(nprim,nprim)
   
       type(gridedgevar_t) :: Qminus, Qplus
       type(gridvar_t) :: Q6
  
positive.code:                     (Q%data(i,n) - &
                             0.5_dp_t*(Qminus%data(i,n) + Qplus%data(i,n))) ) then
  
               Qplus%data(i,n) = 3.0_dp_t*Q%data(i,n) - 2.0_dp_t*Qminus%data(i,n)
  
            endif
  
negative.code: 
            else if ( (Qplus%data(i,n) - Qminus%data(i,n)) * &
                      (Q%data(i,n) - &
                        0.5_dp_t*(Qminus%data(i,n) + Qplus%data(i,n))) > &
                     (Qplus%data(i,n) - Qminus%data(i,n))**2/6.0_dp_t ) then
  
            !else if (abs(Qminus%data(i,n) - Q%data(i,n)) >= &
            !     2.0*abs(Qplus%data(i,n) - Q%data(i,n))) then
               Qminus%data(i,n) = 3.0_dp_t*Q%data(i,n) - 2.0_dp_t*Qplus%data(i,n)
  
            else if (-(Qplus%data(i,n) - Qminus%data(i,n))**2/6.0_dp_t > &
  