#include "dims.h"
      SUBROUTINE AURHT(ALFA,FLUX,DRIZL,t7,t9,IM,IHEM)
C     ****     CALCULATES HORIZONTAL VARIATION OF AURORAL HEAT SOURCE
C     ****      THIS IS THE NEW VERSION WITH UPDATED ALFA ---- 10/86
!
! Called from heelis.f:
!                alfa  flux  drizl alfa2 flux2
!     call aurht(t6(3),t8(3),t3(3),t7(3),t9(3),imax,IHEM)
!
      implicit none
!
#include "params.h"
#include "cflowv3.h"
#include "ioncr.h"
#include "ovalr.h"
#include "phys.h"
!
      real :: alfa3,flux3
      COMMON /AURXTRA/ ALFA3(ZIMXP),FLUX3(ZIMXP)
!$OMP THREADPRIVATE (/aurxtra/)
!DIR$ TASKCOMMON AURXTRA
!
! Args:
      integer,intent(in) :: im
      real,intent(out) :: ALFA(IM),FLUX(IM),DRIZL(IM),t7(IM),t9(IM)
      integer,intent(in) :: IHEM
!
! Local:
      integer :: i
!
! In tgcm22, ALFA0 is 2.5
      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) =  ALFA0*(1.-RALFA*WK1(I))
    1   CONTINUE
!
! ALFA0 <= 0.01:
      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)*EXP(-( (ABS(ALON(I)-RAROT1)-IFIX(ABS(ALON(I)-
C    2     RAROT1)/180.)*360.)/RAHT1)**2)+ALF2 * EXP(-( (WK3(I)-RD2) / 
C    3     RAH2)**2)*EXP(-( (ABS(ALON(I)-RAROT2)-IFIX(ABS(ALON(I)-RAROT2)
C    4     /180.)*360.) / 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) = 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))/
     5     RT21)**2)
    2   CONTINUE
      ENDIF

!
      DO 3 I=1,IM
      FLUX(I)=E0*(1.-REE*WK1(I))*EXP(-(    WK3(I)    /WK2(I))**2)
     1 / (2. * ALFA(I) * 1.602E-9)
      DRIZL(I)=EXP(-((     WK3(I) +ABS(    WK3(I)  ))/(2.*H0))**2)
      t7(i) = ALFA20*(1.-RALFA2*WK1(I))
      t9(i) = E20*(1.-RE2*WK1(I))*EXP(-(    WK3(I)    /WK2(I))**2)
     | / (2. * t7(i) * 1.602E-9)
C  The above line was commented out until 12/91, effectively making FLUX
      ALFA3(I+2) = ALFA30
      FLUX3(I+2) = E30*EXP(-(WK3(I)/WK2(I))**2)/1.602E-6

    3 CONTINUE
C  Insert periodic points for ALFA3 and FLUX3
      ALFA3(1) = ALFA3(IM+1)
      ALFA3(2) = ALFA3(IM+2)
      FLUX3(1) = FLUX3(IM+1)
      FLUX3(2) = FLUX3(IM+2)
      ALFA3(IM+3) = ALFA3(3)
      ALFA3(IM+4) = ALFA3(4)
      FLUX3(IM+3) = FLUX3(3)
      FLUX3(IM+4) = FLUX3(4)
      RETURN
      END
C