c     ============================================
      subroutine setaux(maxmx,maxmy,mbc,mx,my,xlower,ylower,dxc,dyc,
     &                  maux,aux)
c     ============================================
c
c     # set auxiliary arrays for shallow water equations on the sphere
c
c     # on input, (xc(i),yc(j)) gives uniformly spaced computational grid.
c     # on output, 
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     
      implicit double precision (a-h,o-z)
      parameter (maxm3 = 1005)
      dimension xc(-3:maxm3), yc(-3:maxm3)
      dimension xp(-3:maxm3,-3:maxm3), yp(-3:maxm3,-3:maxm3)
      dimension zp(-3:maxm3,-3:maxm3)
      dimension theta(-3:maxm3,-3:maxm3), phi(-3:maxm3,-3:maxm3)
      dimension aux(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,maux)
      common /comsphere/ Rsphere
 
      pi = 4.d0*datan(1.d0)
c
      if (mbc .gt. 4) then
	 write(6,*)'***  increase size of local arrays in setaux ***'
	 stop
	 endif
      if (mx+mbc+1.gt.maxm3 .or. my+mbc+1.gt.maxm3) then
	 write(6,*)'***  increase size of maxm3 in setaux ***'
	 stop
	 endif
c
c     # Set xc and yc so that (xc(i),yc(j)) is the 
c     # lower left corner of (i,j) cell in computational space:
c
      do 10 i=1-mbc,mx+mbc+1
         xc(i) = xlower + (i-1.d0) * dxc
   10    continue
c
      do 12 j=1-mbc,my+mbc+1
         yc(j) = ylower + (j-1.d0) * dyc
   12    continue

c     # compute cell corners on sphere and angles phi, theta
c     # related to latitude and longitude
c
      do 15 j=1-mbc,my+mbc+1
         do 15 i=1-mbc,mx+mbc+1

c           # map computational point to (xp,yp,zp) on sphere:
            call mapc2m(xc(i),yc(j),xp(i,j),yp(i,j),zp(i,j))

c           # compute longitude theta from positive x axis:
            r = dsqrt(xp(i,j)**2 + yp(i,j)**2)
            theta(i,j) = atan2(xp(i,j),yp(i,j))
c           # compute phi, angle down from north pole:
            rrat = dmax1(dmin1(r/Rsphere,1.d0), -1.d0)
            if (zp(i,j) .gt. 0.d0) then
                phi(i,j) = pi/2.d0 - dacos(rrat)  
              else
                phi(i,j) = pi/2.d0 + dacos(rrat)  
              endif
c           write(17,1701) i,j,r,phi(i,j),theta(i,j)
   15	    continue

c
      areamax = 0.d0
      areamin = 1.d3
      areatotal = 0.d0 
      do 20 j=1-mbc,my+mbc
         do 20 i=1-mbc,mx+mbc
c
c           # compute normal and tangent vectors to left edge (in tangent plane)
c
c           # tangent vector is edge vector:
            etx = xp(i,j+1) - xp(i,j)
            ety = yp(i,j+1) - yp(i,j)
            etz = zp(i,j+1) - zp(i,j)
            aux(i,j,5) = etx
            aux(i,j,6) = ety
            aux(i,j,7) = etz

c           # normal to sphere in radial direction at midpoint of edge:
            erx = 0.5d0*(xp(i,j) + xp(i,j+1))
            ery = 0.5d0*(yp(i,j) + yp(i,j+1))
            erz = 0.5d0*(zp(i,j) + zp(i,j+1))
c
c           # normal to edge in tangent plane is cross product of et and er:
            enx = ety*erz - etz*ery
            eny = etz*erx - etx*erz
            enz = etx*ery - ety*erx
            ennorm = dsqrt(enx**2 + eny**2 + enz**2)
            aux(i,j,2) = enx/ennorm
            aux(i,j,3) = eny/ennorm
            aux(i,j,4) = enz/ennorm


c           # compute normal and tangent vectors to bottom edge (in tang. pl.)
c
c           # tangent vector is edge vector:
            etx = xp(i+1,j) - xp(i,j)
            ety = yp(i+1,j) - yp(i,j)
            etz = zp(i+1,j) - zp(i,j)
            aux(i,j,11) = etx
            aux(i,j,12) = ety
            aux(i,j,13) = etz

c           # normal to sphere in radial direction at midpoint of edge:
            erx = 0.5d0*(xp(i,j) + xp(i+1,j))
            ery = 0.5d0*(yp(i,j) + yp(i+1,j))
            erz = 0.5d0*(zp(i,j) + zp(i+1,j))
c
c           # normal to edge in tangent plane is cross product of er and et:
            enx = ery*etz - erz*ety
            eny = erz*etx - erx*etz
            enz = erx*ety - ery*etx
            ennorm = dsqrt(enx**2 + eny**2 + enz**2)
            aux(i,j,8) = enx/ennorm
            aux(i,j,9) = eny/ennorm
            aux(i,j,10) = enz/ennorm

c           # normal to sphere in radial direction at cell center:
            erx = 0.25d0*(xp(i,j) + xp(i+1,j) 
     &                            + xp(i,j+1) + xp(i+1,j+1))
            ery = 0.25d0*(yp(i,j) + yp(i+1,j) 
     &                            + yp(i,j+1) + yp(i+1,j+1))
            erz = 0.25d0*(zp(i,j) + zp(i+1,j) 
     &                            + zp(i,j+1) + zp(i+1,j+1))
            ernorm = dsqrt(erx**2 + ery**2 + erz**2)
            aux(i,j,14) = erx/ernorm
            aux(i,j,15) = ery/ernorm
            aux(i,j,16) = erz/ernorm

c           # compute area of physical cell from four corners:
c           # find area on the sphere of two spherical triangles obtained
c           # by subdividing rectangle.  See
c           #    http://mathforum.org/library/drmath/view/65316.html

c           # corners are labeled  
c           #      1: (i,j)      2: (i+1,j)    3: (i,j+1)    4: (i+1,j+1)

            beta12 = dsin(phi(i,j))*dsin(phi(i+1,j))*
     &                  dcos(theta(i,j) - theta(i+1,j))
     &                + dcos(phi(i,j))*dcos(phi(i+1,j))
            beta23 = dsin(phi(i,j+1))*dsin(phi(i+1,j))*
     &                  dcos(theta(i,j+1) - theta(i+1,j))
     &                 + dcos(phi(i,j+1))*dcos(phi(i+1,j))
            beta13 = dsin(phi(i,j+1))*dsin(phi(i,j))*
     &                  dcos(theta(i,j+1) - theta(i,j))
     &                 + dcos(phi(i,j+1))*dcos(phi(i,j))
            beta24 = dsin(phi(i+1,j+1))*dsin(phi(i+1,j))*
     &                  dcos(theta(i+1,j+1) - theta(i+1,j))
     &                 + dcos(phi(i+1,j+1))*dcos(phi(i+1,j))
            beta34 = dsin(phi(i+1,j+1))*dsin(phi(i,j+1))*
     &                  dcos(theta(i+1,j+1) - theta(i,j+1))
     &                 + dcos(phi(i+1,j+1))*dcos(phi(i,j+1))

c           # great circles distances between corners:
            d12 = Rsphere * dacos(beta12)
            d23 = Rsphere * dacos(beta23)
            d13 = Rsphere * dacos(beta13)
            d24 = Rsphere * dacos(beta24)
            d34 = Rsphere * dacos(beta34)

            s123 = 0.5d0 * (d12 + d23 + d13)
            s234 = 0.5d0 * (d23 + d34 + d24)

c           # spherical excess for each triangle:
            t123 = dtan(s123/2.d0)*dtan((s123-d12)/2.d0)
     &             *dtan((s123-d23)/2.d0)*dtan((s123-d13)/2.d0)
            t123 = dmax1(t123,0.d0)
            E123 = 4.d0*datan(sqrt(t123))

            t234 = dtan(s234/2.d0)*dtan((s234-d23)/2.d0)
     &             *dtan((s234-d34)/2.d0)*dtan((s234-d24)/2.d0)
            t234 = dmax1(t234,0.d0)
            E234 = 4.d0*datan(sqrt(t234))

c           area = Rsphere**2 * (E123 + E234)
            area = (E123 + E234)
c
c           # capacity kappa:
	    aux(i,j,1) = area / (dxc*dyc)
            if (i.gt.0 .and. i.le.mx .and. j.gt.0 .and. j.le.my) then
               areamax = dmax1(areamax,area)
               areamin = dmin1(areamin,area)
               areatotal = areatotal+aux(i,j,1)*dxc*dyc
               endif
c           write(17,1701) i,j,area,phi(i,j),theta(i,j),dxc,areamax
 1701       format(2i5,5d16.6)
   20       continue
c      write(6,*) 'areamax = ',areamax,'  areamin = ',areamin
c      write(6,*) 'area ratio = ',areamax/areamin
c      write(6,*) 'areatotal = ', areatotal, '  |4pi-areatotal| = ',
c    &       dabs(4.d0*pi-areatotal)

       return


       end
