SUBROUTINE AURHT(ALFA,FLUX,DRIZL,IM)
      implicit none
C     ****     CALCULATES HORIZONTAL VARIATION OF AURORAL HEAT SOURCE
C     ****      THIS IS THE NEW VERSION WITH UPDATED ALFA ---- 10/86
      include "params.h"
      include "cons.h"
      include "cflowv3.h"

!     include "phys.h" ! for debug only

      integer :: istar,ihem
      real :: theta0,dduumm
      COMMON /IONCR/ ISTAR,IHEM,THETA0(2),
     1  DDUUMM(IMAXMP*JMX0+IMX0*JMX0+30)
      real :: rrad,h0,rh,rroth,e0,ree,rrote,fc,alfac,fd,alfad,
     |  alfk,alf6,alf21,rrot6,rrot21,rd6,rd6v,rd21,rh6,rh21,rt6,
     |  rt21,alfa0,ralfa,alfa20,ralfa2,e20,re2
      COMMON /OVALR/ RRAD(2),H0,RH,RROTH,E0,REE,RROTE,FC,ALFAC,FD,ALFAD
     1, ALFK,ALF6,ALF21,RROT6,RROT21,RD6,RD6V,RD21,RH6,RH21,RT6,RT21
     2, ALFA0,RALFA,ALFA20,RALFA2,E20,RE2
c
c Low energy auroral proton input (see ALFALP and EFLUXLP from input):
      real :: alfa_lp,eflux_lp,qteaur
      common/aurlp/ alfa_lp(zimxp),eflux_lp(zimxp),qteaur(zimxp)
!DIR$ TASKCOMMON aurlp
!
! Args:
      real,intent(out) :: ALFA(ZIMXP,2),FLUX(ZIMXP,2),DRIZL(1)
      integer,intent(in) :: im
!
! Local:
      real :: pi,alfap1,alfap2,pe1,pe2,alfap0,e0p
      integer :: i
C
      PI = 3.1415926535898
      IF (ALFA0 .GT. 0.01) THEN
C     ****     OLD ALFA:  ALFA(I) = ALFA0*(1.-RALFA*WK1(I))
      DO 1 I=1,IM
C     ****
C     ****     WK1 = COS(LAMDA)
C     ****
      WK1(I)=COS(ATAN2(SIN(ALON(I)-RROTE),COS(ALON(I)-RROTE)))
C     ****
C     ****     WK2 = AURORAL HALF WIDTH
C     ****
      WK2(I)=H0*(1.-RH*COS(ATAN2(SIN(ALON(I)-RROTH),
     1 COS(ALON(I)-RROTH))) )
C     ****
C     ****     WK3 = COLAT - RAD
C     ****
      WK3(I) = COLAT(I) - RRAD(IHEM)
C     ****     OLD ALFA:  ALFA(I) = ALFA0*(1.-RALFA*WK1(I))

      ALFA(I,1) =  ALFA0*(1.-RALFA*WK1(I))
    1 CONTINUE
      ELSE
C     ****     NEW ALFA (10/86):
C     ****      THE EQUATORWARD DISPLACEMENT OF HARDER ENERGY AROUND
C     ****         10 MLT IS   RD1=RD1A+RD1V*COS(ALON(I))
      DO 2 I=1,IM
C     ****
C     ****     WK1 = COS(LAMDA)
C     ****
      WK1(I)=COS(ATAN2(SIN(ALON(I)-RROTE),COS(ALON(I)-RROTE)))
C     ****
C     ****     WK2 = AURORAL HALF WIDTH
C     ****
      WK2(I)=H0*(1.-RH*COS(ATAN2(SIN(ALON(I)-RROTH),
     1 COS(ALON(I)-RROTH))) )
C     ****
C     ****     WK3 = COLAT - RAD
C     ****
      WK3(I) = COLAT(I) - RRAD(IHEM)
C     ****     NEW ALFA (10/86):
C     ****      THE EQUATORWARD DISPLACEMENT OF HARDER ENERGY AROUND
C     ****         10 MLT IS   RD1=RD1A+RD1V*COS(ALON(I))
C     ALFA(I,1) = ALFK+ALF1*EXP(-((WK3(I)-(RD1A+RD1V*COS(ALON(I))))
C    1  / RAH1)**2)
C    2 *EXP(-( (ABS(ALON(I)-RAROT1)-IFIX(ABS(ALON(I)-RAROT1)/180.)*360.)
C    3 /RAHT1)**2)     +     ALF2 * EXP(-( (WK3(I)-RD2) / RAH2)**2)
C    4 *EXP(-( (ABS(ALON(I)-RAROT2)-IFIX(ABS(ALON(I)-RAROT2)/180.)*360.)
C    5 / RAHT2)**2)
C     ****
C     ****     ANOTHER NEW ALFA (3/89):
C     ****     ALON IS BETWEEN -180 (0 LT) AND +180.  ROT6 IS AROUND
C     ****      -90. AND ROT21 IS AROUND 135.  NEED DIFFERENCE 0-180.
C     ****
      ALFA(I,1) = ALFK
     1  +  ALF6 * EXP(-( (WK3(I)-RD6-RD6V*COS(ALON(I)) )/RH6)**2) *
     2 EXP(-( ATAN2(SIN(ALON(I)-RROT6),COS(ALON(I)-RROT6))/RT6 )**2 )
     3  + ALF21 * EXP(-( (WK3(I)-RD21)/RH21)**2) *
     4 EXP(-( ATAN2(SIN(ALON(I)-RROT21),COS(ALON(I)-RROT21))/RT21 )**2 )
    2 CONTINUE
      ENDIF

!     write(6,"('aurht: j=',i3,' e0=',e12.4)") j,e0

      DO 3 I=1,IM
      FLUX(I,1)=E0*(1.-REE*WK1(I))*EXP(-(    WK3(I)    /WK2(I))**2)
     1 / (2. * ALFA(I,1) * 1.602E-9)
      DRIZL(I)=EXP(-((     WK3(I) +ABS(    WK3(I)  ))/(2.*H0))**2)
      ALFA(I,2) = ALFA20*(1.-RALFA2*WK1(I))
      FLUX(I,2)=E20*(1.-RE2*WK1(I))*EXP(-(    WK3(I)    /WK2(I))**2)
     | / (2. * ALFA(I,2) * 1.602E-9)
C  The above line was commented out until 8/92, effectively making FLUX2
c
c Low-energy proton inputs:
      ALFAP1 = 5.
      ALFAP2 = 15.
C     PE1 = 0.15
C     PE2 = 0.4
      PE1 = 1.E-20
      PE2 = 1.E-20
      ALFAP0 = 0.5*(ALFAP1+ALFAP2)
      alfa_lp(i+2) = ALFAP0*(1-RALFA2*WK1(I))
      E0P = 0.5*(PE1+PE2)
      eflux_lp(i+2) = E0P*(1.-REE*WK1(I))*exp(-(wk3(i)/wk2(i))**2)
c
      qteaur(i+2) = -7.e+8*exp(-(wk3(i)/wk2(i))**2)
    3 CONTINUE
c
c Periodic points for alfa_lp and eflux_lp:
      alfa_lp(1)     = alfa_lp(im+1)
      alfa_lp(2)     = alfa_lp(im+2)
      alfa_lp(im+3)  = alfa_lp(3)
      alfa_lp(im+4)  = alfa_lp(4)
      eflux_lp(1)    = eflux_lp(im+1)
      eflux_lp(2)    = eflux_lp(im+2)
      eflux_lp(im+3) = eflux_lp(3)
      eflux_lp(im+4) = eflux_lp(4)
      qteaur(1)    = qteaur(im+1)
      qteaur(2)    = qteaur(im+2)
      qteaur(im+3) = qteaur(3)
      qteaur(im+4) = qteaur(4)
      RETURN
      END
C