#include "dims.h"
      SUBROUTINE AMIEPA (VELNS,VELEW,POTEN,CUSP,ALFA,FLUX,DRIZL,WION)
      use cons_module,only: imax
      implicit none
C
C  GET ION VELOCITIES, POTENTIAL AND AURORAL OVAL FROM AMIE RESULTS
C
#include "params.h"
#include "cflowv3.h"
#include "phys.h"
!
! Args:
      real,intent(out) :: velns(zimx),velew(zimx),poten(zimx),
     |  cusp(zimx),alfa(zimxp,2),flux(zimxp,2),drizl(zimx),wion(zimx)
      real :: rlatm,rlonm,dumid,rlatmp,rlonmp,dlons
      COMMON/TRGM/RLATM(ZIMXP,ZJMX),RLONM(ZIMXP,ZJMX),
     | DUMID(ZIMXP,ZJMX,6),RLATMP(ZJMX),RLONMP(ZJMX),DLONS(ZJMX)
!
#include "ovalr.h"
!
      real :: alfa3,flux3
      COMMON /AURXTRA/ ALFA3(ZIMXP),FLUX3(ZIMXP)
!$OMP THREADPRIVATE (/aurxtra/)
!DIR$ TASKCOMMON AURXTRA
!
      integer :: iamie,iuts1,iuts2,nspsec,ntims,ntst
      real :: sectgcm,crad,phida,secuta,hpsha,hpnha,cpsha,cpnha,
     |  cuspsta,cuspnta,cuspsla,cuspnla,potkva1,vixyza1,ekeva1,
     |  eflxa1,potkva2,vixyza2,ekeva2,eflxa2
      COMMON/AMIE/IAMIE,IUTS1,IUTS2,NSPSEC,NTIMS,NTST,SECTGCM,CRAD(2),
     1 PHIDA(2),SECUTA(2),HPSHA(2),HPNHA(2),CPSHA(2),CPNHA(2),
     2 CUSPSTA(2),CUSPNTA(2),CUSPSLA(2),CUSPNLA(2),
     3 POTKVA1(73,36),VIXYZA1(73,36,3),EKEVA1(73,36),EFLXA1(73,36),
     4 POTKVA2(73,36),VIXYZA2(73,36,3),EKEVA2(73,36),EFLXA2(73,36)
!
! Local:
      real :: dlat(zimx),dlon(zimx),seca,f1,f2
      real,parameter :: S5=0.08726646, S10=0.174532925, S20=0.34906585
      real,parameter :: PI=3.14159265358979
      integer :: ih,i,im
!
      IH = 1
      IF (J .GE. ZJMX/2) IH = 2
C  GET AURORAL OVAL
      DO 1000 I=1,IMAX
      DLAT(I) = RLATM(I+2,J)
 1000 DLON(I) = RLONM(I+2,J)-DLONS(J)
C  CALCULATE COLAT AND ALON FOR CUSP AND DRIZZLE
C  VALUES IN OVALPOS WILL ALREADY BE DETERMINED FROM CALL TO FLOWXX FOR
C   S. H.
      DO 2000 I=1,IMAX
      SINLAT(I) = SIN(ABS(DLAT(I)))
      COSLAT(I) = COS(DLAT(I))
      SINLON(I) = SIN(DLON(I))
      COSLON(I) = COS(DLON(I))
      COLAT(I) = ACOS( SINLAT(I) )
 2000 ALON(I) = AMOD(ATAN2(+SINLON(I)*COSLAT(I),COSLAT(I)*COSLON(I))
     | +3.*PI,2.*PI) - PI
C  SUBSTITUTE IN AMIE VALUES
C  POTKV IN KV, POTEN IN V
C  VIXYZA IN M/S, VELEW,NS,WION IN CM/S
C  EKEVA (2*ALFA) IN KEV, ALFA IN KEV
C  EFLXA IN W/M2, NFLUX ASSUMES IT IS IN ERG/CM2-S
      SECA = AMAX1(SECUTA(1),SECTGCM)
      SECA = AMIN1(SECUTA(2),SECA)
      F1 = SECA - SECUTA(1)
      F2 = SECUTA(2) - SECA
      DO 3000 I=1,IMAX
      POTEN(I) = (F2*POTKVA1(I,J)+F1*POTKVA2(I,J))*1000./FLOAT(NSPSEC)
      VELEW(I) = (F2*VIXYZA1(I,J,1)+F1*VIXYZA2(I,J,1)) * 100. /
     | FLOAT(NSPSEC)
      VELNS(I) = (F2*VIXYZA1(I,J,2)+F1*VIXYZA2(I,J,2)) * 100. /
     | FLOAT(NSPSEC)
      WION(I) = (F2*VIXYZA1(I,J,3)+F1*VIXYZA2(I,J,3)) * 100. /
     | FLOAT(NSPSEC)
      ALFA(I,1) = (F2*EKEVA1(I,J)+F1*EKEVA2(I,J))*0.5/FLOAT(NSPSEC)
C  WARNING --- THIS WORKS ONLY BECAUSE ALFA21=ALFA22=ALFA20
      ALFA(I,2) = ALFA20
      ALFA3(I+2) = ALFA30
C  CUSP AND DRIZL ARE FRACTIONS, NUMBERS BETWEEN 0 AND 1.
C  ON RECOMMENDATION OFROD HEELIS, PLACE THE CUSP AT THE REVERSAL
C   BOUNDARY WITH A RADIUS OF 5 DEGREES
      CUSP(I)=(EXP(-((CRAD(IH)-COLAT(I))/S5)**2)+EXP(-((PI-CRAD(IH)-
     |COLAT(I))/S5)**2))*EXP(-(ATAN2(SIN(ALON(I)-PHIDA(IH)),COS(ALON(I)
     | - PHIDA(IH)))/S20)**2)
 3000 DRIZL(I) = EXP(-(( COLAT(I)-CRAD(IH) + ABS(COLAT(I)-CRAD(IH)) )
     | / S10)**2)
      DO 4000 I=1,IMAX
      FLUX3(I+2) = 1.E-20
      FLUX(I,1) = (F2*EFLXA1(I,J)+F1*EFLXA2(I,J))*1000. / (2.*ALFA(I,1)
     | * 1.602E-9 * FLOAT(NSPSEC) )
C  7/92  The following flux is much too large, and immediately develops
C        large Te's and low Ne's in the auroral region
C     FLUX(I,2) = (F2*EFLXA1(I,J)+F1*EFLXA2(I,J))*1000. / (2.*ALFA20
C    | * 1.602E-9 * FLOAT(NSPSEC) )
C  WK1 = COS(LAMDA)
      WK1(I) = COS(ATAN2(SIN(ALON(I)-RROTE),COS(ALON(I)-RROTE)))
C  WK2 = auroral half-width
      WK2(I) = H0*(1.-RH*COS(ATAN2(SIN(ALON(I)-RROTH),COS(ALON(I)-
     | RROTH))) )
C  WK3 = colat - rad
      WK3(I) = COLAT(I) - CRAD(IH)
 4000 FLUX(I,2) = E20*(1.-RE2*WK1(I)) * EXP(-(WK3(I)/WK2(I))**2)
     | / (2. * ALFA20 * 1.602E-9)
C  Insert periodic points for ALFA3 and FLUX3
      im = imax ! added 10/10/00 (btf)
      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