anchor.diff:  c
               climit = zero
   c            climit = one
  -c            climit = 1.e-1
  +CC            climit = 1.e-1
               alpha_LF(iel) = climit * max(abs(dot_product(u0(iel,:)-um0(iel,:),nv0(iel,:))-c0(iel)),
        &                  abs(dot_product(u1(iel,:)-um1(iel,:),nv1(iel,:))-c1(iel)))
               alpha = alpha_LF(iel)
  
positive.code: c
        deallocate(tmpmu0,tmpmu1)
  c
              call kinematic_condition(ri0,y0,y1)
              call kinematic_condition(ri1,y1,y0)
  c
  c...LHS calculations...
  c
  
negative.code:       endif
  #endif
  c
  c...UPWIND????
  c
  C      ri0(iel,16:20) = ri0(iel,16:20) + f1n0(1:5)
  