anchor.diff:                          iff = (i-1)*refratio + ii
                           jf = (j-1)*refratio + jj
                           qf = qc + shiftx*gradx + shifty*grady -
  -     &                        aux_fine(19,iff,jf)
  +     &                        aux_fine(mbathy,iff,jf)
                           qfine(mq,iff,jf) = qf
                        enddo
                     enddo
  
positive.code:                if (mq .eq. 1) then
  c                 # Interpolate sea surface height rather than just the
  c                 # water column height.
                    qc = qcoarse(mq,ic,jc) + aux_coarse(19,ic,jc)
  
                    sl = (qc - qcoarse(mq,ic-1,jc) -
       &                  aux_coarse(19,ic-1,jc))
                    sr = (qcoarse(mq,ic+1,jc)+aux_coarse(19,ic+1,jc) - qc)
                    gradx = compute_slopes(sl,sr,mth)
  
                    sl = (qc - qcoarse(mq,ic,jc-1) -
       &                  aux_coarse(ic,jc-1,19))
                    sr = (qcoarse(mq,ic,jc+1)+aux_coarse(19,ic,jc+1) - qc)
                    grady = compute_slopes(sl,sr,mth)
  
  c                 # Fill in fine grid values from coarse grid cell (ic,jc)
  
negative.code:       subroutine geoclaw_interpolate2fine(mx,my,mbc,meqn,qcoarse,
       &      qfine,maux,aux_coarse,aux_fine,mbathy,
       &      p4est_refineFactor,refratio,igrid)
        implicit none
  
        integer mx,my,mbc,meqn,p4est_refineFactor,refratio
        integer igrid, maux,mbathy
        double precision qcoarse(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)
        double precision qfine(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)
        double precision aux_coarse(maux,1-mbc:mx+mbc,1-mbc:my+mbc)
  