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:             hadd(m,2,0,i-1) = hadd(m,2,0,i-1)
       &                      - 0.5d0*dtdx1d(i-1)*cpamdq(m,i)
              hadd(m,1,0,i)   = hadd(m,1,0,i)
       &                      - 0.5d0*dtdx1d(i-1)*cmapdq(m,i)
              hadd(m,2,0,i)   = hadd(m,2,0,i)
       &                      - 0.5d0*dtdx1d(i-1)*cpapdq(m,i)
  c
  c           # Transverse propagation of the increment wave (and the
  c           # correction wave if m4=2) between cells
  
negative.code:             gadd(m,2,0,i-1) = gadd(m,2,0,i-1)
       &                      - 0.5d0*dtdx1d(i-1)*bpamdq(m,i)
              gadd(m,1,0,i)   = gadd(m,1,0,i)
       &                      - 0.5d0*dtdx1d(i)*bmapdq(m,i)
              gadd(m,2,0,i)   = gadd(m,2,0,i)
       &                      - 0.5d0*dtdx1d(i)*bpapdq(m,i)
  c
  c           # Transverse propagation of the increment wave (and the
  c           # correction wave if m4=2) between cells
  