#include "dims.h"
      SUBROUTINE ADHEEL
      implicit none
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"
#include "ceee.h"
#include "cterp.h"
#include "dynphi.h"
#include "consdyn.h"
#include "ioncr.h"
!
! Locals:
      integer :: ncc,imx,jmx,n
!
! For print_cee only:
      real :: C0(IMX0,JMX0,10),C1(IMX1,JMX1,9),C2(IMX2,JMX2,9),
     1  C3(IMX3,JMX3,9),C4(IMX4,JMX4,9)
      EQUIVALENCE (CEE,C0),(CEE(NC1),C1),(CEE(NC2),C2),(CEE(NC3),C3),
     1  (CEE(NC4),C4)
!
      crit(1) = 0.261799387
      crit(2) = 0.523598775
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     ****     MODIFY STENCILS
C     ****
!     call print_cee(c0,c1,c2,c3,c4) ! cee3

      NCC = 1
      IMX = IMX0
      JMX = JMX0
      DO 5 N = 1,5
        CALL STENMOD(IMX,JMX,CEE(NCC),PHIHM(1,JMX0),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

!     call print_cee(c0,c1,c2,c3,c4) ! cee3

      RETURN
      END
C