
c
c
c     =====================================================
      subroutine rpn2(ixy,maxm,meqn,mwaves,mbc,mx,ql,qr,auxl,auxr,
     &			wave,s,amdq,apdq)
c     =====================================================
c
c     # Roe-solver for the 2D shallow water equations
c     #  on the sphere, using 3d Cartesian representation of velocities
c
c     # solve Riemann problems along one slice of data.
c
c     # On input, ql contains the state vector at the left edge of each cell
c     #           qr contains the state vector at the right edge of each cell
c
c     # This data is along a slice in the x-direction if ixy=1
c     #                            or the y-direction if ixy=2.
c     # On output, wave contains the waves, s the speeds,
c     # and amdq, apdq the decomposition of the flux difference
c     #   f(qr(i-1)) - f(ql(i))
c     # into leftgoing and rightgoing parts respectively.
c     # With the Roe solver we have
c     #    amdq  =  A^- \Delta q    and    apdq  =  A^+ \Delta q
c     # where A is the Roe matrix.  An entropy fix can also be incorporated
c     # into the flux differences.
c
c     # Note that the i'th Riemann problem has left state qr(i-1,:)
c     #                                    and right state ql(i,:)
c     # From the basic clawpack routines, this routine is called with ql = qr
c
c
      implicit double precision (a-h,o-z)
c
      dimension wave(1-mbc:maxm+mbc, meqn, mwaves)
      dimension    s(1-mbc:maxm+mbc, mwaves)
      dimension   ql(1-mbc:maxm+mbc, meqn)
      dimension   qr(1-mbc:maxm+mbc, meqn)
      dimension  apdq(1-mbc:maxm+mbc, meqn)
      dimension  amdq(1-mbc:maxm+mbc, meqn)
      dimension auxl(1-mbc:maxm+mbc, 16)
      dimension auxr(1-mbc:maxm+mbc, 16)
c
c     local arrays -- common block comroe is passed to rpt2
c     ------------
      parameter (maxm2 = 1002)  !# assumes at most 1000x1000 grid with mbc=2
      dimension delta(3)
      logical efix

      common /sw/  g
      common /comroe/ u(-1:maxm2),v(-1:maxm2),a(-1:maxm2),h(-1:maxm2)
      common /comxyt/ dtcom,dxcom,dycom,tcom,icom,jcom
c
      data efix /.true./    !# use entropy fix for transonic rarefactions
c
      if (-1.gt.1-mbc .or. maxm2 .lt. maxm+mbc) then
	 write(6,*) 'need to increase maxm2 in rpA'
	 stop
	 endif
      
      if(ixy.eq.1) then
        dy = dycom
      else
        dy = dxcom
      endif 

c     The aux array has the following elements:
c         1  kappa = ratio of cell area to dxc*dyc
c         2  enx = x-component of normal vector to left edge in tangent plane
c         3  eny = y-component of normal vector to left edge in tangent plane
c         4  enz = z-component of normal vector to left edge in tangent plane
c         5  etx = x-component of tangent vector to left edge in tangent plane
c         6  ety = y-component of tangent vector to left edge in tangent plane
c         7  etz = z-component of tangent vector to left edge in tangent plane
c         8  enx = x-component of normal vector to bottom edge in tangent plane
c         9  eny = y-component of normal vector to bottom edge in tangent plane
c        10  enz = z-component of normal vector to bottom edge in tangent plane
c        11  etx = x-component of tangent vector to bottom edge in tangent plane
c        12  ety = y-component of tangent vector to bottom edge in tangent plane
c        13  etz = z-component of tangent vector to bottom edge in tangent plane
c        14  erx = x-component of unit vector in radial direction at cell ctr
c        15  ery = y-component of unit vector in radial direction at cell ctr
c        16  erz = z-component of unit vector in radial direction at cell ctr
c
c     # offset to index into aux array for enx, eny, etx, ety, gamma 
c     #    depends on whether ixy=1 (left edge) or ixy=2 (bottom edge).
      ioff = 6*(ixy-1) + 1
c
c
c     # find a1 thru a3, the coefficients of the 3 eigenvectors:

      do 20 i = 2-mbc, mx+mbc

         enx =   auxl(i,ioff+1)
         eny =   auxl(i,ioff+2)
         enz =   auxl(i,ioff+3)
         etx =   auxl(i,ioff+4)
         ety =   auxl(i,ioff+5)
         etz =   auxl(i,ioff+6)
         gamma = dsqrt(etx**2 + ety**2 + etz**2)
         etx =   etx / gamma
         ety =   ety / gamma
         etz =   etz / gamma


c        # compute normal and tangential momentum at cell edge:
         hunl = enx*ql(i,2) + eny*ql(i,3) + enz*ql(i,4)
         hunr = enx*qr(i-1,2) + eny*qr(i-1,3) + enz*qr(i-1,4)

         hutl = etx*ql(i,2) + ety*ql(i,3) + etz*ql(i,4)
         hutr = etx*qr(i-1,2) + ety*qr(i-1,3) + etz*qr(i-1,4)

c        # compute the Roe-averaged variables needed in the Roe solver.
c        # These are stored in the common block comroe since they are
c        # later used in routine rpt2 to do the transverse wave splitting.
c
         hl = ql(i,1)
         hr = qr(i-1,1)
         h(i) = (hl+hr)*0.50d0
         hsqr = dsqrt(hr)
         hsql = dsqrt(hl)
         hsq = hsqr + hsql
         u(i) = (hunr/hsqr + hunl/hsql) / hsq
         v(i) = (hutr/hsqr + hutl/hsql) / hsq
         a(i) = dsqrt(g*h(i))

c        # Split the jump in q at each interface into waves
         delta(1) = hl - hr
         delta(2) = hunl - hunr
         delta(3) = hutl - hutr

         a1 = ((u(i)+a(i))*delta(1) - delta(2))*(0.50d0/a(i))
         a2 = -v(i)*delta(1) + delta(3)
         a3 = (-(u(i)-a(i))*delta(1) + delta(2))*(0.50d0/a(i))

c
c        # Compute the waves.
c
         wave(i,1,1) = a1
         wave(i,2,1) = a1*(u(i)-a(i))*enx + a1*v(i)*etx
         wave(i,3,1) = a1*(u(i)-a(i))*eny + a1*v(i)*ety
         wave(i,4,1) = a1*(u(i)-a(i))*enz + a1*v(i)*etz
         s(i,1) = (u(i)-a(i)) * gamma/dy 
c
         wave(i,1,2) = 0.0d0
         wave(i,2,2) = a2*etx
         wave(i,3,2) = a2*ety
         wave(i,4,2) = a2*etz
         s(i,2) = u(i) * gamma/dy 
c
         wave(i,1,3) = a3
         wave(i,2,3) = a3*(u(i)+a(i))*enx + a3*v(i)*etx
         wave(i,3,3) = a3*(u(i)+a(i))*eny + a3*v(i)*ety
         wave(i,4,3) = a3*(u(i)+a(i))*enz + a3*v(i)*etz
         s(i,3) = (u(i)+a(i)) * gamma/dy 
  281    format(2i4,5d12.4)
  283    format(8x,5d12.4)
   20    continue
c
c
c    # compute flux differences amdq and apdq.
c    ---------------------------------------
c
      if (efix) go to 110
c
c     # no entropy fix
c     ----------------
c
c     # amdq = SUM s*wave   over left-going waves
c     # apdq = SUM s*wave   over right-going waves
c
      do 100 m=1,meqn
         do 100 i=2-mbc, mx+mbc
	    amdq(i,m) = 0.d0
	    apdq(i,m) = 0.d0
	    do 90 mw=1,mwaves
	       if (s(i,mw) .lt. 0.d0) then
		   amdq(i,m) = amdq(i,m) + s(i,mw)*wave(i,m,mw)
		 else
		   apdq(i,m) = apdq(i,m) + s(i,mw)*wave(i,m,mw)
		 endif
   90          continue
  100       continue
c
c     # project momentum components of amdq and apdq onto tangent plane:

      do i=2-mbc,mx+mbc  
            erx = auxr(i-1,14)
            ery = auxr(i-1,15)
            erz = auxr(i-1,16)
            amn = erx*amdq(i,2) + ery*amdq(i,3) + erz*amdq(i,4)
            amdq(i,2) = amdq(i,2) - amn*erx
            amdq(i,3) = amdq(i,3) - amn*ery
            amdq(i,4) = amdq(i,4) - amn*erz

            erx = auxl(i,14)
            ery = auxl(i,15)
            erz = auxl(i,16)
            apn = erx*apdq(i,2) + ery*apdq(i,3) + erz*apdq(i,4)
            apdq(i,2) = apdq(i,2) - apn*erx
            apdq(i,3) = apdq(i,3) - apn*ery
            apdq(i,4) = apdq(i,4) - apn*erz
            enddo
      go to 900
c
c-----------------------------------------------------
c
  110 continue
c
c     # With entropy fix
c     ------------------
c
c    # compute flux differences amdq and apdq.
c    # First compute amdq as sum of s*wave for left going waves.
c    # Incorporate entropy fix by adding a modified fraction of wave
c    # if s should change sign.
c
         do 200 i=2-mbc,mx+mbc
            enx =   auxl(i,ioff+1)
            eny =   auxl(i,ioff+2)
            enz =   auxl(i,ioff+3)
            etx =   auxl(i,ioff+4)
            ety =   auxl(i,ioff+5)
            etz =   auxl(i,ioff+6)
            gamma = dsqrt(etx**2 + ety**2 + etz**2)
            etx =   etx / gamma
            ety =   ety / gamma
            etz =   etz / gamma
c           # compute normal and tangential momentum at cell edge:
            hunl = enx*ql(i,2) + eny*ql(i,3) + enz*ql(i,4)
            hunr = enx*qr(i-1,2) + eny*qr(i-1,3) + enz*qr(i-1,4)

c           check 1-wave
            him1 = qr(i-1,1)
            s0 =  (hunr/him1 - dsqrt(g*him1)) * gamma / dy
c           check for fully supersonic case :
            if (s0.gt.0.0d0.and.s(i,1).gt.0.0d0) then
               do 60 m=1,4
                  amdq(i,m)=0.0d0
   60          continue
               goto 200
            endif
c
            h1 = qr(i-1,1)+wave(i,1,1)
            hu1= hunr + enx*wave(i,2,1) + eny*wave(i,3,1) 
     &                                        + enz*wave(i,4,1)
            s1 = (hu1/h1 - dsqrt(g*h1))*gamma/dy  !speed just to right of 1-wave
            if (s0.lt.0.0d0.and.s1.gt.0.0d0) then
c              transonic rarefaction in 1-wave
               sfract = s0*((s1-s(i,1))/(s1-s0))
            else if (s(i,1).lt.0.0d0) then
c              1-wave is leftgoing
               sfract = s(i,1)
            else
c              1-wave is rightgoing
               sfract = 0.0d0
            endif
            do 120 m=1,4
               amdq(i,m) = sfract*wave(i,m,1)
  120       continue
c           check 2-wave
            if (s(i,2).gt.0.0d0) then
c	       #2 and 3 waves are right-going
	       go to 200 
	       endif

            do 140 m=1,4
               amdq(i,m) = amdq(i,m) + s(i,2)*wave(i,m,2)
  140       continue
c
c           check 3-wave
c
            hi = ql(i,1)
            s03 = (hunl/hi + dsqrt(g*hi)) * gamma/dy
            h3=ql(i,1)-wave(i,1,3)
            hu3=hunl-(enx*wave(i,2,3)+eny*wave(i,3,3)+enz*wave(i,4,3))
            s3=(hu3/h3 + dsqrt(g*h3)) * gamma/dy
            if (s3.lt.0.0d0.and.s03.gt.0.0d0) then
c              transonic rarefaction in 3-wave
               sfract = s3*((s03-s(i,3))/(s03-s3))
            else if (s(i,3).lt.0.0d0) then
c              3-wave is leftgoing
               sfract = s(i,3)
            else
c              3-wave is rightgoing
               goto 200
            endif
            do 160 m=1,4
               amdq(i,m) = amdq(i,m) + sfract*wave(i,m,3)
  160       continue
  200       continue
c
c           compute rightgoing flux differences :
c
            do 220 i = 2-mbc,mx+mbc  
               do 222 m=1,4
                  df = 0.0d0
                  do 210 mw=1,mwaves
                     df = df + s(i,mw)*wave(i,m,mw)
  210             continue
                  apdq(i,m)=df - amdq(i,m)
 222              continue
c
c                 project momentum components onto tangent plane
c
                  erx = auxr(i-1,14)
                  ery = auxr(i-1,15)
                  erz = auxr(i-1,16)
                  amn = erx*amdq(i,2)+ery*amdq(i,3)+erz*amdq(i,4)
                  amdq(i,2) = amdq(i,2) - amn*erx
                  amdq(i,3) = amdq(i,3) - amn*ery
                  amdq(i,4) = amdq(i,4) - amn*erz

                  erx = auxl(i,14)
                  ery = auxl(i,15)
                  erz = auxl(i,16)
                  apn = erx*apdq(i,2)+ery*apdq(i,3)+erz*apdq(i,4)
                  apdq(i,2) = apdq(i,2) - apn*erx
                  apdq(i,3) = apdq(i,3) - apn*ery
                  apdq(i,4) = apdq(i,4) - apn*erz

  220          continue

  900 continue
      return
      end


