      SUBROUTINE STENCIL(ZIGM,CS,IMX0,JMX0,SYM,CEE,ARRAY,NCOEF,WK)
      implicit none
      SAVE
C     ****
C     ****     CALCULATES CONTRIBUTION TO 3 BY 3 STENCIL FROM
C     ****       COEFFICIENT ,ZIGM, AT EACH GRID POINT AND LEVEL
C     ****
C     ****     INPUT:
C     ****       ZIGM((IMX0+1)/2,(JMX0+1)/2)   ARRAY OF COEFFICIENT
C     ****         VALUES.
C     ****       CS(-15:JMX0) ARRAY OF COSINES OF LATITUDE VALUES
C     ****       IMX0, JMX0  LONGITUDE AND LATITUDE DIMENSIONS
C     ****         OF FINEST GRID LEVEL
C     ****       SYM = 1.  IF ZIGM SYMMETRIC W.R.T. EQUATOR
C     ****           = -1.  FOR ANTISYMMETRY
C     ****       NCOEF = NUMBER IDENTIFYING COEFFICIENT
C     ****         NCOEF = 1 FOR ZIGM11
C     ****         NCOEF = 2 FOR ZIGM12 (=ZIGMC+ZIGM2)
C     ****         NCOEF = 3 FOR ZIGM21 (=ZIGMC-ZIGM2)
C     ****         NCOEF = 4 FOR ZIGM22
C     ****
C     ****     OUTPUT:
C     ****       CEE    STENCIL ARRAY CONSISTING OF
C     ****         C0(IMX0,JMX0,10), C1(IMX1,JMX1,9), C2(IMX2,JMX2,9),
C     ****         C3(IMX3,JMX3,9), C4(IMX4,JMX4,9)
C     ****           WHERE:   IMX1 = (IMX0+1)/2, JMX1 = (JMX0+1)/2,
C     ****                    IMX2 = (IMX1+1)/2, JMX2 = (JMX1+1)/2,
C     ****              AND SO ON
C     ****
C     ****     WORK SPACE:
C     ****       ARRAY(-15:IMX0+16,-15:JMX0+2)
C     ****         USED FOR INTERPOLATED AND EXTENDED ARRAY OF
C     ****         COEFFICIENT ZIGM
C     ****
!
! Args:
      integer,intent(in) :: imx0,jmx0,ncoef
      real,intent(in) :: ZIGM(IMX0,JMX0), CS(1),CEE(1),sym,
     1  ARRAY(IMX0+32,1),wk(1)
!
! Local:
      integer :: nc,imx,jmx,n
C     ****
C     ****     PERFORM HALF WAY INTERPOLATION AND EXTEND ZIGM IN ARRAY
C     ****
      CALL HTRPEX(ZIGM,ARRAY,IMX0,JMX0,SYM)
C     ****
C     ****     CALCULATE CONTRIBUTION TO STENCIL FOR EACH GRID POINT
C     ****       AND LEVEL
C     ****
      NC = 1
      IMX = IMX0
      JMX = JMX0
      DO 1 N = 1,5
        CALL CNM(ARRAY,CS,IMX0,JMX0,IMX,JMX,CEE(NC),NCOEF,WK)
        NC = NC+9*IMX*JMX
        IF(N.EQ.1)NC = NC+IMX*JMX
        IMX = (IMX+1)/2
        JMX = (JMX+1)/2
    1 CONTINUE
      RETURN
      END
