anchor.diff:  c        # For correction terms below, need average of dtdx in cell
   c        # i-1 and i.  Compute these and overwrite dtdx1d:
   c
  -         dtdx1d(i-1) = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i))
  +         dtdxave = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i))
   c
            do 120 m=1,meqn
               cqxx(m,i) = 0.d0
               do 119 mw=1,mwaves
                  cqxx(m,i) = cqxx(m,i) + 0.5d0 * dabs(s(mw,i))
  -     &             * (1.d0 - dabs(s(mw,i))*dtdx1d(i-1)) * wave(m,mw,i)
  +     &             * (1.d0 - dabs(s(mw,i))*dtdxave) * wave(m,mw,i)
     119          continue
               fadd(m,i) = fadd(m,i) + cqxx(m,i)
     120       continue
  
positive.code: c
              if(m3.lt.2) go to 200
                 hadd(m,2,0,i)   = hadd(m,2,0,i)
       &                         + dtdx1d(i-1)*cpcqxxp(m,i)
                 hadd(m,1,0,i)   = hadd(m,1,0,i)
       &                         + dtdx1d(i-1)*cmcqxxp(m,i)
                 hadd(m,2,0,i-1) = hadd(m,2,0,i-1)
       &                         - dtdx1d(i-1)*cpcqxxm(m,i)
                 hadd(m,1,0,i-1) = hadd(m,1,0,i-1)
  
negative.code: c
  
                  gadd(m,2,0,i) = gadd(m,2,0,i)
       &                  + (1.d0/6.d0)*dtdx1d(i)*dtdz
       &                  * (bpcpapdq(m,i) - bpcmapdq(m,i))
                  gadd(m,1,0,i) = gadd(m,1,0,i)
       &                  + (1.d0/6.d0)*dtdx1d(i)*dtdz
       &                  * (bmcpapdq(m,i) - bmcmapdq(m,i))
  
  
                  gadd(m,2,1,i) = gadd(m,2,1,i)
       &                          - (1.d0/6.d0)*dtdx1d(i)*dtdz
       &                          * bpcpapdq(m,i)
                  gadd(m,1,1,i) = gadd(m,1,1,i)
       &                          - (1.d0/6.d0)*dtdx1d(i)*dtdz
       &                          * bmcpapdq(m,i)
                  gadd(m,2,-1,i) = gadd(m,2,-1,i)
       &                          + (1.d0/6.d0)*dtdx1d(i)*dtdz
       &                          * bpcmapdq(m,i)
                  gadd(m,1,-1,i) = gadd(m,1,-1,i)
       &                          + (1.d0/6.d0)*dtdx1d(i)*dtdz
       &                          * bmcmapdq(m,i)
  c
                  gadd(m,2,0,i-1) = gadd(m,2,0,i-1)
  