anchor.diff:  c                    # interior cell. 
                        delta = neighborval + gpthis2(i1,mq) 
        &                          - efthis2(i1,mq,0)
  -                     areathis = area2(i1)
  +                     area = area2(i1)
  +
  +                     write(6,100) this_iface, i1, neighborval, 
  +     &                       gpthis2(i1,mq), 
  +     &                       efthis2(i1,mq,0), delta                     
   
   c                    # Reset these, since they will no longer be needed                     
                        gpthis2(i1,mq) = 0
                        efthis2(i1,mq,0) = 0
  -                     efthis2(i1,mq,1) = 0
  +c                     efthis2(i1,mq,1) = 0
                     else
   c                    # Coarse grid is bottom grid; fine is top grid   
   c                    # efc3(0) is flux stored in the coarse grid 
   c                    # interior cell.                   
                        delta = neighborval + gmthis3(i1,mq) + 
        &                            efthis3(i1,mq,0)
  -                     areathis = area3(i1)
  +                     area = area3(i1)
  +
  +                     write(6,100) this_iface, i1, neighborval, 
  +     &                       gmthis3(i1,mq), 
  +     &                       efthis3(i1,mq,0), delta                     
  +
   
   c                    # Reset these, since they will no longer be needed                     
                        gmthis3(i1,mq) = 0
                        efthis3(i1,mq,0) = 0
  -                     efthis3(i1,mq,1) = 0
  +c                     efthis3(i1,mq,1) = 0
                     endif
  -                  qthis(i1,j1,mq) = qthis(i1,j1,mq) + delta/areathis                 
  -               enddo
  -            enddo
  -         endif
  -      enddo
  +                  qthis(i1,j1,mq) = qthis(i1,j1,mq) + 0.5*delta/area                 
  +              enddo
  +          enddo
  +      endif
         end
   
   
  
positive.code:      
              deltap = (fpfine1(jj1,mq) + fpfine1(jj2,mq)) -
       &               (eff1(jj1,mq,1) + eff1(jj2,mq,1))
  
              if (j .eq. 1) then
  c                write(6,200) fmfine0(jj1,mq),eff0(jj1,mq,1)
              endif
   200        format(2F16.8)
  
  c           # Put the same value in four ghost cells;  grab the first
  c           # one that shows up in the transformation.          
  
negative.code: c           # NOTE : Fluxes have already been scaled by gamma, e.g. 
  c           # x-face-length/dy or y-face-length/dx
  
              fp_left(j,m)  = fp_left(j,m) - dt*dy*fp(1,j,m)
              fm_left(j,m)  = fm_left(j,m) + dt*dy*fm(1,j,m)
  
              fp_right(j,m) = fp_right(j,m)-dt*dy*fp(mx+1,j,m)
              fm_right(j,m) = fm_right(j,m)+dt*dy*fm(mx+1,j,m)
  