anchor.diff:              if( m4.gt.0 )then
   c
                   hadd(m,2,0,i)  = hadd(m,2,0,i)
  -     &                  + (1.d0/6.d0)*dtdx1d(i-1)*dtdy
  +     &                  + (1.d0/6.d0)*dtdx1d(i)*dtdy
        &                  * (bpcpapdq(m,i) - bpcmapdq(m,i))
                   hadd(m,1,0,i)  = hadd(m,1,0,i)
  -     &                  + (1.d0/6.d0)*dtdx1d(i-1)*dtdy
  +     &                  + (1.d0/6.d0)*dtdx1d(i)*dtdy
        &                  * (bmcpapdq(m,i) - bmcmapdq(m,i))
   
   
                   hadd(m,2,1,i)  = hadd(m,2,1,i)
  -     &                         - (1.d0/6.d0)*dtdx1d(i-1)*dtdy
  +     &                         - (1.d0/6.d0)*dtdx1d(i)*dtdy
        &                         * bpcpapdq(m,i)
                   hadd(m,1,1,i)  = hadd(m,1,1,i)
  -     &                         - (1.d0/6.d0)*dtdx1d(i-1)*dtdy
  +     &                         - (1.d0/6.d0)*dtdx1d(i)*dtdy
        &                         * bmcpapdq(m,i)
                   hadd(m,2,-1,i) = hadd(m,2,-1,i)
  -     &                         + (1.d0/6.d0)*dtdx1d(i-1)*dtdy
  +     &                         + (1.d0/6.d0)*dtdx1d(i)*dtdy
        &                         * bpcmapdq(m,i)
                   hadd(m,1,-1,i) = hadd(m,1,-1,i)
  -     &                         + (1.d0/6.d0)*dtdx1d(i-1)*dtdy
  +     &                         + (1.d0/6.d0)*dtdx1d(i)*dtdy
        &                         * bmcmapdq(m,i)
   c
                   hadd(m,2,0,i-1)  = hadd(m,2,0,i-1)
  
positive.code: c
  
                  gadd(m,2,0,i) = gadd(m,2,0,i)
       &                  + (1.d0/6.d0)*dtdx1d(i-1)*dtdz
       &                  * (bpcpapdq(m,i) - bpcmapdq(m,i))
                  gadd(m,1,0,i) = gadd(m,1,0,i)
       &                  + (1.d0/6.d0)*dtdx1d(i-1)*dtdz
       &                  * (bmcpapdq(m,i) - bmcmapdq(m,i))
  
  
                  gadd(m,2,1,i) = gadd(m,2,1,i)
       &                          - (1.d0/6.d0)*dtdx1d(i-1)*dtdz
       &                          * bpcpapdq(m,i)
                  gadd(m,1,1,i) = gadd(m,1,1,i)
       &                          - (1.d0/6.d0)*dtdx1d(i-1)*dtdz
       &                          * bmcpapdq(m,i)
                  gadd(m,2,-1,i) = gadd(m,2,-1,i)
       &                          + (1.d0/6.d0)*dtdx1d(i-1)*dtdz
       &                          * bpcmapdq(m,i)
                  gadd(m,1,-1,i) = gadd(m,1,-1,i)
       &                          + (1.d0/6.d0)*dtdx1d(i-1)*dtdz
       &                          * bmcmapdq(m,i)
  c
                  gadd(m,2,0,i-1) = gadd(m,2,0,i-1)
  
negative.code: c
              if(m3.lt.2) go to 180
                 gadd(m,2,0,i)   = gadd(m,2,0,i)
       &                         + dtdx1d(i)*bpcqxxp(m,i)
                 gadd(m,1,0,i)   = gadd(m,1,0,i)
       &                         + dtdx1d(i)*bmcqxxp(m,i)
                 gadd(m,2,0,i-1) = gadd(m,2,0,i-1)
       &                         - dtdx1d(i-1)*bpcqxxm(m,i)
                 gadd(m,1,0,i-1) = gadd(m,1,0,i-1)
  