#include "dims.h"
      SUBROUTINE POLHT(QP,IM,ihem,j)
      implicit none
C     ****     CALCULATES HORIZONTAL VARIATION OF POLAR CUSP HEATING
#include "params.h"
#include "ioncr.h"
#include "cflowv3.h"
      real,parameter :: S5=.08726646, S20=.34906585, PI=3.1415927
!
! Args:
      real,intent(out) :: QP(1)
      integer,intent(in) :: im,ihem,j
!
! Local:
      integer :: i
      real :: rad
!
C     ****    ON RECOMMENDATION OF ROD HEELIS, PLACE THE CUSP AT THE
C     ****     REVERSAL BOUNDARY WITH A RADIUS OF 5 DEGREES
!
      RAD=THETA0(IHEM)
!
! QP, defined here, is CUSP, passed from heelis as T5:
!
      DO 1 I=1,IM
      QP(I)=(EXP(-((RAD-COLAT(I))/S5)**2)+EXP(-((PI-RAD-COLAT(I))/S5)
     1**2))*EXP(-(ATAN2(SIN(ALON(I)-PHID(IHEM)),
     2 COS(ALON(I)-PHID(IHEM)))/S20)**2)

!         write(6,"('cusp: lat=',i3,' i=',i3,' ihem=',i2,' theta0=',
!    |      e12.4,' colat=',e12.4,' alon=',e12.4,' phid=',e12.4,
!    |      ' cusp=',e12.4)") j,i,ihem,theta0(ihem),colat(i), 
!    |      alon(i),phid(ihem),qp(i) 

    1 CONTINUE

!     write(6,"('aurora_cusp: lat=',i2,' cusp=',/,(6e12.4))") j,qp(1:im)

      RETURN
      END
C