anchor.diff:           do j = 1,my/2
               jj1 = 2*(j-1) + 1     !! indexing starts at 1
               jj2 = 2*(j-1) + 2
  -            deltam = (fm0(jj1,mq) + fm0(jj2,mq) +
  -     &               eff0(jj1,mq) + eff0(jj2,mq))
  -            deltap = (fp1(jj1,mq) + fp1(jj2,mq) +
  -     &               eff1(jj1,mq) + eff1(jj2,mq))
  +            deltam = fmfine0(jj1,mq) + fmfine0(jj2,mq) +
  +     &               eff0(jj1,mq,1) + eff0(jj2,mq,1)
  +            deltap = fpfine1(jj1,mq) + fpfine1(jj2,mq) +
  +     &               eff1(jj1,mq,1) + eff1(jj2,mq,1)
   
               do jj = jj1,jj2
                  do ii = 1,mbc
  
positive.code:                if (.not. skip_this_grid) then
                    fineval = qfine_dummy(i2(1),j2(1),mq)
                    if (iface_coarse .eq. 2) then
                       deltac = fineval + gpcoarse2(ic,mq) + efc2(ic,mq)
                       areac = area2(ic)
                    else
                       deltac = fineval + gmcoarse3(ic,mq) + efc3(ic,mq)
                       areac = area3(ic)
                    endif
                    qcoarse(ic,jc,mq) = qcoarse(ic,jc,mq) + deltac/areac
  
negative.code:                if (.not. skip_this_grid) then
                    fineval = qfine_dummy(i2(1),j2(1),mq)
                    if (iface_coarse .eq. 0) then
                       deltac = fineval + fpcoarse0(jc,mq) + efc0(jc,mq,0)
                       areac = area0(jc)
                    else
                       deltac = fineval + fmcoarse1(jc,mq) + efc1(jc,mq,0)
                       areac = area1(jc)
                    endif
                    qcoarse(ic,jc,mq) = qcoarse(ic,jc,mq) + deltac/areac
  