c     ============================================
      subroutine setaux(maxmx,maxmy,mbc,mx,my,xlower,ylower,dxc,dyc,
     &                  maux,aux)
c     ============================================
c
c     # set auxiliary arrays for advection on a curvilinear grid
c
c     # aux(i,j,1) cell centered velocity in the x-direction
c     # aux(i,j,2) cell centered velocity in the y-direction 
c     # aux(i,j,3) metric term
c
c     # we denote the different sections on the ccomputational grid by
c     # 1, 2, 3, 4, 1*, 2*, 3*, 4*
c 
c     |------------|------------|
c     |\    2*   / |\    2    / |
c     |  \     /   |  \     /   |
c     |1*  \ /  3* | 3  \ /   1 |
c     |    / \     |    / \     |
c     |  /     \   |  /     \   |
c     |/   4*    \ |/    4    \ |
c     |------------|------------|
c
c
      implicit double precision (a-h,o-z)

      parameter (maxm3 = 1005)
      dimension aux(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc, maux)
      dimension xccorn(5),yccorn(5),xpcorn(5),ypcorn(5)
      dimension xc(-3:maxm3), yc(-3:maxm3)
      dimension xp(-3:maxm3,-3:maxm3), yp(-3:maxm3,-3:maxm3),
     &          zp(-3:maxm3,-3:maxm3)
      dimension theta(-3:maxm3,-3:maxm3), phi(-3:maxm3,-3:maxm3)

      external Dxyz

      common /dmetric/ dmetric(-1:300+2,-1:150+2,5)
c
c
      if(mx.gt.300.or.my.gt.150) then
        write(6,*) "need to increase dmetric array in setaux and src2"
        stop
      endif
  
      pi2 = 8.d0*datan(1.d0)  !# = 2 * pi
      pi = 0.5*pi2

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
      do 13 j=1-mbc,my+mbc+1
         do 13 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))
 13      continue
c
c     set aux(.,.,1) and aux(.,.,2): advection speeds for solid body rotation
      do 14 j=1-mbc, my+mbc
         do 14 i=1-mbc, mx+mbc !# use for instance factor 0.025  for 
c                              !# non-zero solid body rotation
             aux(i,j,1)=0.0*(stream(xp(i,j+1),yp(i,j+1),zp(i,j+1))
     &                    - stream(xp(i,j),yp(i,j),zp(i,j))) / dyc

             aux(i,j,2)=-0.0*(stream(xp(i+1,j),yp(i+1,j),zp(i+1,j))
     &                     - stream(xp(i,j),yp(i,j),zp(i,j))) / dxc
 14        continue
 

       
      do 20 j=1,my
         do 20 i=1,mx
           xi = xlower+(i-0.5)*dxc
           yj = ylower+(j-0.5)*dyc

           if(xi.ge.-1.d0) then
             xx = xi
             istar = 1  !# upper part
           else
             xx = -2.d0-xi
             istar = -1  !# lower part
           endif
c          case 1 or 1*:
           if(xx.gt.0.d0.and.dabs(yj).le.xx) then
             dj11 = istar*Dxyz(1,xx,yj)
             dj12 = Dxyz(2,xx,yj)
             dj21 = istar*Dxyz(3,xx,yj)
             dj22 = Dxyz(4,xx,yj)
             dj31 = istar*Dxyz(5,xx,yj)
             dj32 = Dxyz(6,xx,yj)
             goto 15
           endif         
c          case 2 or 2*:
           if(yj.ge.0.d0.and.dabs(xx).le.yj) then
             dj11 = istar*Dxyz(4,yj,xx)
             dj12 = Dxyz(3,yj,xx)
             dj21 = istar*Dxyz(2,yj,xx)
             dj22 = Dxyz(1,yj,xx)
             dj31 = istar*Dxyz(6,yj,xx)
             dj32 = Dxyz(5,yj,xx)
             goto 15   
           endif  
c          case 3 or 3*:
           if(xx.le.0.d0.and.dabs(yj).le.dabs(xx)) then
             dj11 = istar*Dxyz(7,xx,yj)
             dj12 = Dxyz(8,xx,yj)
             dj21 = istar*Dxyz(9,xx,yj)
             dj22 = Dxyz(10,xx,yj)
             dj31 = istar*Dxyz(11,xx,yj)
             dj32 = Dxyz(12,xx,yj)
             goto 15
           endif
c          case 4 or 4*:
           if(yj.le.0.d0.and.dabs(xx).le.dabs(yj)) then
             dj11 = istar*Dxyz(10,yj,xx)
             dj12 = Dxyz(9,yj,xx)
             dj21 = istar*Dxyz(8,yj,xx)
             dj22 = Dxyz(7,yj,xx)
             dj31 = istar*Dxyz(12,yj,xx)
             dj32 = Dxyz(11,yj,xx)
             goto 15
           endif

 15        continue
c          along the 'diagonals' we calculate metric terms by averaging
c          between the two possible formulas
           i_diag_fix = 0  !# assume grid cell are not on 'diagonal'
c
c          fix along diagonal between region 1 and 2 or 1* and 2*,respectively:
           if(xx.ge.0.d0.and.yj.ge.0.d0.and.dabs(xx-yj).le.1.d-10) then
             i_diag_fix = 1
             dj11_1 = istar*Dxyz(1,xx,yj)
             dj12_1 = Dxyz(2,xx,yj)
             dj21_1 = istar*Dxyz(3,xx,yj)
             dj22_1 = Dxyz(4,xx,yj)
             dj31_1 = istar*Dxyz(5,xx,yj)
             dj32_1 = Dxyz(6,xx,yj)
c
             dj11_2 = istar*Dxyz(4,yj,xx)
             dj12_2 = Dxyz(3,yj,xx)
             dj21_2 = istar*Dxyz(2,yj,xx)
             dj22_2 = Dxyz(1,yj,xx)
             dj31_2 = istar*Dxyz(6,yj,xx)
             dj32_2 = Dxyz(5,yj,xx)

           endif   
c          fix along diagonal between region 2 and 3 or 2* and 3*,respectively:
           if(xx.ge.-1.d0.and.xx.le.0.d0.and.yj.gt.0.d0.and.
     &        dabs(-xx-yj).le.1.d-10) then
             i_diag_fix = 1
             dj11_1 = istar*Dxyz(4,yj,xx)
             dj12_1 = Dxyz(3,yj,xx)
             dj21_1 = istar*Dxyz(2,yj,xx)
             dj22_1 = Dxyz(1,yj,xx)
             dj31_1 = istar*Dxyz(6,yj,xx)
             dj32_1 = Dxyz(5,yj,xx)
c
             dj11_2 = istar*Dxyz(7,xx,yj)
             dj12_2 = Dxyz(8,xx,yj)
             dj21_2 = istar*Dxyz(9,xx,yj)
             dj22_2 = Dxyz(10,xx,yj)
             dj31_2 = istar*Dxyz(11,xx,yj)
             dj32_2 = Dxyz(12,xx,yj)          
           endif
c          fix along diagonal between region 3 and 4 or 3* and 4*,respectively:
           if(xx.ge.-1.d0.and.xx.lt.0.d0.and.yj.lt.0.d0.and.
     &        dabs(xx-yj).le.1.d-10) then  
             i_diag_fix = 1
             dj11_1 = istar*Dxyz(7,xx,yj)
             dj12_1 = Dxyz(8,xx,yj)
             dj21_1 = istar*Dxyz(9,xx,yj)
             dj22_1 = Dxyz(10,xx,yj)
             dj31_1 = istar*Dxyz(11,xx,yj)
             dj32_1 = Dxyz(12,xx,yj)
c
             dj11_2 = istar*Dxyz(10,yj,xx)
             dj12_2 = Dxyz(9,yj,xx)
             dj21_2 = istar*Dxyz(8,yj,xx)
             dj22_2 = Dxyz(7,yj,xx)
             dj31_2 = istar*Dxyz(12,yj,xx)
             dj32_2 = Dxyz(11,yj,xx)              
           endif   
c          fix along diagonal between region 4 and 1 or 4* and 1*,respectively:
           if(xx.ge.0.d0.and.yj.le.0.d0.and.dabs(xx+yj).le.1.d-10) then 
             i_diag_fix = 1
             dj11_1 = istar*Dxyz(10,yj,xx)
             dj12_1 = Dxyz(9,yj,xx)
             dj21_1 = istar*Dxyz(8,yj,xx)
             dj22_1 = Dxyz(7,yj,xx)
             dj31_1 = istar*Dxyz(12,yj,xx)
             dj32_1 = Dxyz(11,yj,xx)                   
c        
             dj11_2 = istar*Dxyz(1,xx,yj)
             dj12_2 = Dxyz(2,xx,yj)
             dj21_2 = istar*Dxyz(3,xx,yj)
             dj22_2 = Dxyz(4,xx,yj)
             dj31_2 = istar*Dxyz(5,xx,yj)
             dj32_2 = Dxyz(6,xx,yj) 
           endif    
c
           if(i_diag_fix.eq.0) then !# not on diagonal
            det_h = (dj11**2.d0+dj21**2.d0+dj31**2.d0)
     &       *(dj12**2.d0+dj22**2.d0+dj32**2.d0) - (dj12*dj11+
     &         dj22*dj21+dj32*dj31)**2.d0 
            h11 = (dj12**2.d0+dj22**2.d0+dj32**2.d0)/det_h
            h12 = -(dj11*dj12+dj21*dj22+dj31*dj32)/det_h
            h21 = -(dj12*dj11+dj22*dj21+dj32*dj31)/det_h 
            h22 = (dj11**2.d0+dj21**2.d0+dj31**2.d0)/det_h
           else  !# for grid cells on diagonal
            det_h1 = (dj11_1**2.d0+dj21_1**2.d0+dj31_1**2.d0)
     &       *(dj12_1**2.d0+dj22_1**2.d0+dj32_1**2.d0) - (dj12_1*dj11_1+
     &         dj22_1*dj21_1+dj32_1*dj31_1)**2.d0 
            det_h2 = (dj11_2**2.d0+dj21_2**2.d0+dj31_2**2.d0)
     &       *(dj12_2**2.d0+dj22_2**2.d0+dj32_2**2.d0) - (dj12_2*dj11_2+
     &         dj22_2*dj21_2+dj32_2*dj31_2)**2.d0
            det_h = 0.5d0*(det_h1+det_h2)
            h11_1 = (dj12_1**2.d0+dj22_1**2.d0+dj32_1**2.d0)/det_h1
            h12_1 = -(dj11_1*dj12_1+dj21_1*dj22_1+dj31_1*dj32_1)/det_h1
            h21_1 = -(dj12_1*dj11_1+dj22_1*dj21_1+dj32_1*dj31_1)/det_h1 
            h22_1 = (dj11_1**2.d0+dj21_1**2.d0+dj31_1**2.d0)/det_h1
            h11_2 = (dj12_2**2.d0+dj22_2**2.d0+dj32_2**2.d0)/det_h2
            h12_2 = -(dj11_2*dj12_2+dj21_2*dj22_2+dj31_2*dj32_2)/det_h2
            h21_2 = -(dj12_2*dj11_2+dj22_2*dj21_2+dj32_2*dj31_2)/det_h2 
            h22_2 = (dj11_2**2.d0+dj21_2**2.d0+dj31_2**2.d0)/det_h2
            h11 = 0.5*(h11_1+h11_2)
            h12 = 0.5*(h12_1+h12_2)
            h21 = 0.5*(h21_1+h21_2)
            h22 = 0.5*(h22_1+h22_2) 
           endif   
           dmetric(i,j,1) =  det_h
c          #  h^{11}
           dmetric(i,j,2) = h11 
c          #  h^{12}
           dmetric(i,j,3) = h12 
c          #  h^{21}
           dmetric(i,j,4) = h21 
c          #  h^{22} 
           dmetric(i,j,5) = h22 
c
           aux(i,j,3) = dsqrt(det_h)
 20    continue


c      # set bc's at top and bottom
       do i=1, mx
         i2 = mx+1-i
         aux(i,0,3) = aux(i2,1,3) 
         aux(i,-1,3) = aux(i2,2,3)         
         aux(i,my+1,3) = aux(i2,my,3)
         aux(i,my+2,3) = aux(i2,my-1,3)
         do k=1, 5
            dmetric(i,0,k) = dmetric(i2,1,k) 
            dmetric(i,-1,k) = dmetric(i2,2,k)         
c
            dmetric(i,my+1,k) = dmetric(i2,my,k)
            dmetric(i,my+2,k) = dmetric(i2,my-1,k)            
         enddo
       enddo
c      # set periodic  boundary conditions
       do j=1-mbc, my+mbc
         do ibc=1, mbc
             aux(1-ibc,j,3) = aux(mx+1-ibc,j,3)
             aux(mx+ibc,j,3)= aux(ibc,j,3)            
           do k=1, 5
             dmetric(1-ibc,j,k) = dmetric(mx+1-ibc,j,k)
             dmetric(mx+ibc,j,k)= dmetric(ibc,j,k)  
           enddo 
         enddo
       enddo
      
c
c      # calculate area on surface of the sphere:
       
       dsum = 0.d0
       do i=1, mx
         do j=1, my
           dsum = dsum+aux(i,j,3)*dxc*dyc      
         enddo
       enddo
       write(6,*) "error in the calculation of the surface",dsum-4.d0*pi


       return
       end

