#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,jj,j0,j,i
      
!     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     ****
!     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

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    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))

      RETURN
      END
C
C--------------------------------------------------------------------------
#include "dims.h"
      SUBROUTINE ADHEELMD
!      
! am_02/02:   adheelmd: adheel for modified mudpack solver modmud
!             (isolve = 2 in dyncal)
!             subroutine stenmd: calculate modified and unmodified
!             coefficient stencil
!
      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,jj,j0,j,i
!
!     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     ****
!     NCC = 1
!     IMX = IMX0
!     JMX = JMX0
!     DO 5 N = 1,5
!      CALL STENMD(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


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    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))
      RETURN
      END
C