#include "dims.h"
      SUBROUTINE POLHT(QP,IM,ihem)
      implicit none
C     ****     CALCULATES HORIZONTAL VARIATION OF POLAR CUSP HEATING
#include "params.h"
#include "ioncr.h"
#include "cflowv3.h"
!
! Args:
      integer,intent(in) :: im,ihem
      real,intent(out) :: QP(1)
!
! Local:
      real,parameter :: S5=.08726646, S20=.34906585, PI=3.1415927
      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 passed as CUSP through flowv2, flowv1, flowx0,
! flowxx, and finally as T5 in heelis, which is ref'd by orora. 
!
      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)
    1 CONTINUE
      RETURN
      END
C