c
c     =====================================================
      subroutine qcor(ixy,i,m,aux,q,maxm,meqn,mbc,qc)
c     =====================================================
c
c     # maintain conservation on sphere
c
      implicit double precision (a-h,o-z)

      dimension aux(1-mbc:maxm+mbc, *)
      dimension  q(1-mbc:maxm+mbc, meqn)
      dimension qc(4)

      common /comxyt/ dtcom,dxcom,dycom,tcom,icom,jcom
      common /sw/  g
c
c
      if(ixy.eq.1) then
        in = 2
        dy = dycom
      else
        in = 8
        dy = dxcom 
      endif
c     # left interface 
      etxl = aux(i,in+3)
      etyl = aux(i,in+4)
      etzl = aux(i,in+5)
      gammal = dsqrt(etxl**2.d0 + etyl**2.d0 + etzl**2.d0) / dy
      enxl = aux(i,in) * gammal
      enyl = aux(i,in+1) * gammal
      enzl = aux(i,in+2) * gammal
    

c     # right interface
      etxr = aux(i+1,in+3)
      etyr = aux(i+1,in+4)
      etzr = aux(i+1,in+5)        
      gammar = dsqrt(etxr**2.d0 + etyr**2.d0 + etzr**2.d0) / dy
      enxr = aux(i+1,in) * gammar
      enyr = aux(i+1,in+1) * gammar
      enzr = aux(i+1,in+2)  * gammar
c
      qc(1) = (enxr-enxl)*q(i,2)
     &       +(enyr-enyl)*q(i,3)
     &       +(enzr-enzl)*q(i,4)

      qc(2) = (enxr-enxl)*(q(i,2)**2.d0/q(i,1)+0.5d0*g*q(i,1)**2.d0)
     &       +(enyr-enyl)*(q(i,2)*q(i,3)/q(i,1))
     &       +(enzr-enzl)*(q(i,2)*q(i,4)/q(i,1))
    
      qc(3) = (enxr-enxl)*(q(i,2)*q(i,3)/q(i,1))
     &       +(enyr-enyl)*(q(i,3)**2.d0/q(i,1)+0.5d0*g*q(i,1)**2.d0)
     &       +(enzr-enzl)*(q(i,3)*q(i,4)/q(i,1))
   
      qc(4) = (enxr-enxl)*(q(i,2)*q(i,4)/q(i,1))
     &       +(enyr-enyl)*(q(i,3)*q(i,4)/q(i,1))
     &       +(enzr-enzl)*(q(i,4)**2.d0/q(i,1)+0.5d0*g*q(i,1)**2.d0)

c     # project qc to tangent plane:
      erx = aux(i,14)
      ery = aux(i,15)
      erz = aux(i,16)

      qcn = erx*qc(2) + ery*qc(3) + erz*qc(4)
      qc(2) = qc(2) - qcn*erx
      qc(3) = qc(3) - qcn*ery
      qc(4) = qc(4) - qcn*erz
      return
      end
