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"
      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
!
! Args:
      real,intent(out) :: ALFA(ZIMXP,2),FLUX(ZIMXP,2),DRIZL(1)
      integer,intent(in) :: im
!
! Local:
      real :: pi
      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))

!     write(6,"(/'aurht: i=',i2,' ihem=',i2)") i,ihem
!     write(6,"('  alfa0=',e12.4,' alon(i)=',e12.4)") alfa0,alon(i)
!     write(6,"('  rrote=',e12.4,' h0=',e12.4,' rh=',e12.4)")
!    |  rrote,h0,rh
!     write(6,"('  rroth=',e12.4,' rrad(ihem)=',e12.4)") 
!    |  rroth,rrad(ihem)
!     write(6,"('  ralfa=',e12.4)") ralfa
 
    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
!
      DO 3 I=1,IM

!     write(6,"('e0=',e12.4,' ree=',e12.4,' alfa(i)=',e12.4)")
!    |  e0,ree,alfa(i,1)
!     write(6,"('wk1=',e12.4,' wk2=',e12.4,' wk3=',e12.4)")
!    |  wk1(i),wk2(i),wk3(i)

      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

!     write(6,"('aurht do 3: i=',i2,' flux=',e12.4,' drizl=',e12.4)")
!    |  i,flux(i,1),drizl(i)
!     write(6,"('  flux2=',e12.4,' alfa2=',e12.4)")
!    |  alfa(i,2),flux(i,2)

    3 CONTINUE
      RETURN
      END
C