#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"
#include "trgm.h"
#include "ovalr.h"
#include "amie.h"
!
! Local:
      real :: VELNS(ZIMX),VELEW(ZIMX),POTEN(ZIMX),WION(ZIMX),
     | ALFA(ZIMXP,2),FLUX(ZIMXP,2),DRIZL(ZIMX),CUSP(ZIMX),dlat(ZIMX),
     | dlon(ZIMX)
      integer :: ih,i
      real :: s5,s10,s20,pi,f1,f2,seca
      DATA S5/0.08726646/, S10/0.174532925/, S20/0.34906585/
      DATA PI/3.14159265358979/
C
      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
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
      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)
      RETURN
      END
C