c#######################################################################
c#
c#                Cadez Differencing Coefficients
c#
c#######################################################################
      subroutine cadez_diffcoef(Xi,JXzc,JXzf,JYzc,JYzf,
     .                          RHOzc,RHOxf,RHOyf,Zzc,CJACyf,
     .                          LapX,LapY,ext,int,fil,XIsing,Alpha)
      implicit none
c-----------------------------------------------------------------------
      integer ext(4),int(4),fil(4)
      real*8 Xi(ext(1):ext(3))
      real*8 JXzc(ext(1):ext(3)),JXzf(ext(1):ext(3))
      real*8 JYzc(ext(2):ext(4)),JYzf(ext(2):ext(4))
      real*8 RHOzc(ext(1):ext(3),ext(2):ext(4))
      real*8 RHOxf(ext(1):ext(3),ext(2):ext(4))
      real*8 RHOyf(ext(1):ext(3),ext(2):ext(4))
      real*8 Zzc(ext(1):ext(3),ext(2):ext(4))
      real*8 CJACyf(ext(1):ext(3),ext(2):ext(4))
      real*8 LapX(ext(1):ext(3),ext(2):ext(4))
      real*8 LapY(ext(1):ext(3),ext(2):ext(4))
      real*8 XIsing,Alpha
c-----------------------------------------------------------------------
      integer i,j
      integer imini,imaxi,jmini,jmaxi
      integer iminf,imaxf,jminf,jmaxf
      integer imine,imaxe,jmine,jmaxe
      real*8 radfac
      real*8 ZERO,HALF,ONE,TWO,FOUR,PI,TWOPI
      parameter(ZERO = 0.D0,HALF = 0.5D0,ONE = 1.D0,
     .          TWO = 2.D0,FOUR = 4.D0)
      parameter(PI = 3.1415926535897931160D0, TWOPI = 2.D0*PI)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
      imini = int(1)
      imaxi = int(3)
      jmini = int(2)
      jmaxi = int(4)
      iminf = fil(1)
      imaxf = fil(3)
      jminf = fil(2)
      jmaxf = fil(4)
      imine = ext(1)
      imaxe = ext(3)
      jmine = ext(2)
      jmaxe = ext(4)
c
c-----------------------------------------------------------------------
c- Set Laplacian on X-Zone Faces
c-----------------------------------------------------------------------
      do i = imine,imaxe-1
        do j = jminf,jmaxe
          LapX(i,j) = TWOPI*JXzf(i)*RHOxf(i,j)/JYzc(j)
        end do
      end do
c
c-----------------------------------------------------------------------
c- Set rho on Y-Zone Faces
c-----------------------------------------------------------------------
      do i = iminf,imaxf
        do j = jmine,jmaxe-1
          LapY(i,j) = TWOPI*JYzf(j)*RHOyf(i,j)/JXzc(i)
        end do
      end do
c
c-----------------------------------------------------------------------
c- Correct the Inner Boundary differencing coefficients
c- Save Inner Boundary "Psi"-correction into LapX
c-----------------------------------------------------------------------
      if(jmini.eq.jminf) then
        if(Xi(imini).lt.XIsing) then
          radfac = ONE
        else
          radfac = Alpha
        end if
        j = jmini - 1
        do i = iminf,imaxf
          LapX(i,j) = -TWO*radfac/
     .                        (FOUR*sqrt(CJACyf(i,j))*JYzf(j) - radfac)
        end do
      end if
c
c-----------------------------------------------------------------------
c- Correct the Outer Boundary differencing coefficients
c- Save Outer Boundary "1 + A/r"-correction into LapX
c-----------------------------------------------------------------------
      if(jmaxi.eq.jmaxf) then
        j = jmaxi
        do i = iminf,imaxf
          LapX(i,j+1) = ONE -
     .             sqrt((RHOzc(i,j)**2   + Zzc(i,j)**2  )/
     .                  (RHOzc(i,j+1)**2 + Zzc(i,j+1)**2))
          LapY(i,j) = LapY(i,j)*LapX(i,j+1)
        end do
      end if
c
      return
      end
