      SUBROUTINE ADHEEL
C     ****
C     ****     MODIFY STENCILS AND RHS SO AS TO INSERT HEELIS
C     ****       POTENTIAL FOR GEOMAGNETIC LATITUDE .GE. 60.0 DEG.
C     ****
C     ****
      include "params.h"
      PARAMETER(IMAXG=ZIMX,JMAXG=ZJMX,IMAXGP=IMAXG+1,JMAXGP=JMAXG+1)
      include "cons.h"
      include "ceee.h"
      include "cterp.h"
      include "dynphi.h"
      include "consts.h"
      COMMON /IONCR/ ISTAR,IHEM,THETA0(2),OFFA(2),OFFC(2),DSKOFA(2),
     1  DSKOFC(2),PHID(2),PHIN(2),PHIDP0(2),PHIDM0(2),PHINP0(2),
     2  PHINM0(2),PSIM(2),PSIE(2),PCEN(2),RR1(2),COLATC(IMAXMP,JMX0),
     3  P(IMX0,JMX0),CRIT(2)
C     ****
C     ****     FILL ARRAY COLATC WITH NORTHERN AURORAL COLATITUDES
C     ****       CORRESPONDING TO EACH NORTHERN GEOMAGNETIC GRID
C     ****       POINT.  VALUES IN COLATC(IMX0,JMX0).
C     ****       ALSO CALCULATE FRACTIONAL PRESENCE OF DYNAMO
C     ****       EQUATION GIVEN CRITICAL COLATITUDES CRIT(2).
C     ****       VALUES IN P(IMX0,JMX0).
C     ****
      CALL COLATH
C     ****
C     ****     CALCULATE PHIHM, THE HEELIS POTENTIAL IN
C     ****       GEOMAGNETIC COORDINATES.
C     ****
      ISTAR = 0
      CALL POTM
C     ****
C     ****     Calculate the symmetric and antisymmetric parts of PHIHM
C     ****     Put the symmetric part in as if it were the NH part (equ-pole)
C     ****     (IMX0=IMAXMP=IMAXM+1=80+1=81, JMXO=(JMAXM+1)/2=49, JMAXM=97)
C     ****
      DO 100 I=1,IMAXMP
      PHISYM(I,1)=PHIHM(I,JMX0)
      PHIANT(I,JMX0) = 0.
      DO 100 J=1,JMX0-1
      JJ = JMAXM + 1 - J
      J0 = JMX0 - J + 1
      PHISYM(I,J0) = 0.5 * (PHIHM(I,J) + PHIHM(I,JJ))
      PHIANT(I,J) = 0.5 * (PHIHM(I,J) + PHIHM(I,JJ)) - PHIHM(I,J)
100   PHIANT(I,JJ) = 0.5 * (PHIHM(I,J) + PHIHM(I,JJ)) - PHIHM(I,JJ)
C  Plot of symmetric hi-lat potential
C     CALL EZCNTR (PHISYM,IMX0,JMX0)
C  Plot of anti-symmetric hi-lat potential
C     CALL EZCNTR (PHIANT,IMX0,JMAXM)
C     ****
C     ****     MODIFY STENCILS
C     ****
      NCC = 1
      IMX = IMX0
      JMX = JMX0
      DO 5 N = 1,5
	CALL STENMOD(IMX,JMX,CEE(NCC),PHISYM,P,COLATC)
        NCC = NCC+9*IMX*JMX
        IF(N.EQ.1)NCC = NCC+IMX*JMX
        IMX = (IMX+1)/2
        JMX = (JMX+1)/2
    5 CONTINUE
C  SEEMS TO BE ADDED BACK WRONG!  elim for now
C  1/95:  Alter the anti-symmetric potential by the stencil P
C    P=0 from pole to crit1, P=1 from crit2 to equator
      DO 200 I=1,IMAXMP
      DO 200 J=1,JMX0-1
      JJ = JMAXM + 1 - J
      J0 = JMX0 - J + 1
      PHIANT(I,J) = PHIANT(I,J) * (1. - P(I,J0))
 200  PHIANT(I,JJ) = PHIANT(I,JJ) * (1. - P(I,J0))
C  Plot of anti-symmetric hi-lat potential which will be added back
C     CALL EZCNTR (PHIANT,IMX0,JMAXM)
C  Plot of stencil P used (P=0 from pole to crit1, P=1 from crit2 to equator)
C     CALL EZCNTR (P,IMX0,JMX0)
      RETURN
      END
C
