c#######################################################################
c#
c#                Cadez Differencing Coefficients
c#
c#######################################################################
      subroutine cadez_diffcoef(Xi,dphi,JXzc,JXzf,JYzc,JYzf,
     .                          RHOzc,RHOxf,RHOyf,Zzc,CJACzc,CJACyf,
     .                          LapW,LapX,LapY,ext,int,fil,XIsing,Alpha)
      implicit none
c-----------------------------------------------------------------------
      real*8 dphi
      integer ext(6),int(6),fil(6)
      real*8 Xi(ext(2):ext(5))
      real*8 JXzc(ext(2):ext(5)),JXzf(ext(2):ext(5))
      real*8 JYzc(ext(3):ext(6)),JYzf(ext(3):ext(6))
      real*8 RHOzc(ext(2):ext(5),ext(3):ext(6))
      real*8 RHOxf(ext(2):ext(5),ext(3):ext(6))
      real*8 RHOyf(ext(2):ext(5),ext(3):ext(6))
      real*8 Zzc(ext(2):ext(5),ext(3):ext(6))
      real*8 CJACzc(ext(2):ext(5),ext(3):ext(6))
      real*8 CJACyf(ext(2):ext(5),ext(3):ext(6))
      real*8 LapW(ext(2):ext(5),ext(3):ext(6))
      real*8 LapX(ext(2):ext(5),ext(3):ext(6))
      real*8 LapY(ext(2):ext(5),ext(3):ext(6))
      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,tmp1
      real*8 ZERO,HALF,ONE,TWO,FOUR
      parameter(ZERO = 0.D0,HALF = 0.5D0,ONE = 1.D0,
     .          TWO = 2.D0,FOUR = 4.D0)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
      imini = int(2)
      imaxi = int(5)
      jmini = int(3)
      jmaxi = int(6)
      iminf = fil(2)
      imaxf = fil(5)
      jminf = fil(3)
      jmaxf = fil(6)
      imine = ext(2)
      imaxe = ext(5)
      jmine = ext(3)
      jmaxe = ext(6)
c
c-----------------------------------------------------------------------
c- Set Laplacian on W-Zone Faces
c-----------------------------------------------------------------------
      tmp1 = ONE/dphi
      do i = iminf,imaxf
        do j = jminf,jmaxf
          LapW(i,j) = tmp1/(RHOzc(i,j)*CJACzc(i,j)*JXzc(i)*JYzc(j))
c         LapW(i,j) = ZERO
        end do
      end do
c
c-----------------------------------------------------------------------
c- Set Laplacian on X-Zone Faces
c-----------------------------------------------------------------------
      do i = imine,imaxe-1
        do j = jminf,jmaxe
          LapX(i,j) = dphi*JXzf(i)*RHOxf(i,j)/JYzc(j)
        end do
      end do
c
c-----------------------------------------------------------------------
c- Set Laplacian on Y-Zone Faces
c-----------------------------------------------------------------------
      do i = iminf,imaxf
        do j = jmine,jmaxe-1
          LapY(i,j) = dphi*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
