      SUBROUTINE FLWV32(POTEN)
C     ****
C     ****     COMPUTE ELECTRIC POTENTIAL
C     ****
      include "params.h"
      include "cons.h"
C     ****
C     ****     WORK SPACE
C     ****
      COMMON/VSCR/IFLAG(IMAXM), DLAT(IMAXM), DLON(IMAXM), RATIO(IMAXM),
     1  PSI(8), PHI(IMAXM,8), IFN(IMAXM), PHIFUN(IMAXM), COLAT(IMAXM),
     2  ALON(IMAXM), SINLAT(IMAXM),
     3  COSLAT(IMAXM), SINLON(IMAXM), COSLON(IMAXM), WK1(IMAXM),
     4  WK2(IMAXM),WK3(IMAXM)
      DIMENSION PHIFN2(IMAXM)
      COMMON /IONCR/ ISTAR,IHEM,THETA0(2),OFFA(2),OFFC(2),DSKOFA(2),
     1  DSKOFC(2),PHID(2),PHIN(2),PHIDP0(2),PHIDM0(2),PHINP0(2),
     2  PHINM0(2),PSIM(2),PSIE(2),PCEN(2),R1(2)
     3  ,DDUUMM(IMAXMP*JMX0+IMX0*JMX0+2)
C    3  ,THETAC(2),R2(2)       !   FOR OLD HEELIS POLAR CAP POTENTIAL
      COMMON /OVALPOS/ COSOFA(2),SINOFA(2),ASLONA(2),PI,PI2
      DIMENSION POTEN(IMAXMP)
      DIMENSION PHDPMX(2),PHDMMX(2),PHNPMX(2),PHNMMX(2),
     1  COSOFC(2),SINOFC(2),ASLONC(2)
      DATA E/1.E-6/,PI/3.1415927/,PI2/6.2831853/,PIH/1.5707963/
C     ****
C     ****     CONVERT PARAMETERS TO RADIANS ETC.
C     ****
      IF(ISTAR.EQ.0)THEN
        ISTAR=1
        PI=4.*ATAN(1.)
        PI2=2.*PI
        PIH=.5*PI
        DO 50 N=1,2
	IF (ABS(OFFA(N)+DSKOFA(N)) .GT. 0.01) THEN
        OFDA = SQRT(OFFA(N)**2 + DSKOFA(N)**2)
        COSOFA(N) = COS(OFDA)
        SINOFA(N) = SIN(OFDA)
        ASLONA(N) = ASIN(DSKOFA(N)/OFDA)
	ELSE
	OFDC = 0.
	COSOFA(N) = 1.
	SINOFA(N) = 0.
	ASLONA(N) = 0.
	ENDIF
	IF (ABS(OFFC(N)+DSKOFC(N)) .GT. 0.01) THEN
        OFDC = SQRT(OFFC(N)**2 + DSKOFC(N)**2)
        COSOFC(N) = COS(OFDC)
        SINOFC(N) = SIN(OFDC)
        ASLONC(N) = ASIN(DSKOFC(N)/OFDC)
	ELSE
	OFDC = 0.
	COSOFC(N) = 1.
	SINOFC(N) = 0.
	ASLONC(N) = 0.
	ENDIF
        IF(PHIN(N).LT.PHID(N))PHIN(N)=PHIN(N)+PI2
        PHDPMX(N)=.5*AMIN1(PI,(PHIN(N)-PHID(N)))
        PHNPMX(N)=.5*AMIN1(PI,(PHID(N)-PHIN(N)+PI2))
        PHNMMX(N)=PHDPMX(N)
   50   PHDMMX(N)=PHNPMX(N)
      ENDIF
C     ****               SET IHEM = 1,2  FOR  S,N HEMISPHERE
      IMAXD2 = MAX0(1,IMAX/2)
      IHEM = IFIX(DLAT(IMAXD2)*2./3.1416 + 2.)
      SINTH0=SIN(THETA0(IHEM))
C     ****
C     ****     AVERAGE AMIE RESULTS SHOW R1=-2.6 FOR 11.3 DEGREES
C     ****       (0.1972 RAD) BEYOND THETA0.
C     ****
      SINTHR1 = SIN(THETA0(IHEM)+0.1972)
C     G0=(SIN(THETAC(IHEM))/SINTH0)**R2(IHEM)    !   FOR OLD HEELIS
C     A2=1./((SIN(THETA0(IHEM)+THETAC(IHEM))/SINTH0)**R2(IHEM)-G0)
      PSI(1)=PSIE(IHEM)
      PSI(3)=PSIM(IHEM)
      DO 60 N=2,4,2
        PSI(N)=PSI(N-1)
   60 CONTINUE
      DO 70 N=1,4
        PSI(N+4)=PSI(N)
   70 CONTINUE
C     ****
C     ****     TRANSFORM TO AURORAL CIRCLE CORDINATES
C     ****
      DO 4 I=1,IMAXM
        SINLAT(I)=SIN(ABS(DLAT(I)))
        COSLAT(I)=COS(DLAT(I))
        SINLON(I)=SIN(DLON(I)+ASLONC(IHEM))
        COSLON(I)=COS(DLON(I)+ASLONC(IHEM))
        COLAT(I) = COSOFC(IHEM)*SINLAT(I) - SINOFC(IHEM)*COSLAT(I)*
     1    COSLON(I)
        ALON(I) = AMOD(ATAN2(+SINLON(I)*COSLAT(I),SINLAT(I)*
     1    SINOFC(IHEM)+COSOFC(IHEM)*COSLAT(I)*COSLON(I)) -
     2    ASLONC(IHEM)+3.*PI,PI2) - PI
        COLAT(I)=ACOS(COLAT(I))*SQRT(RATIO(I))
C       ****     CALCULATE BOUNDARIES FOR LONGITUDINAL FUNCTION
        WK1(I)=((COLAT(I)-THETA0(IHEM))/THETA0(IHEM))**2
        PHI(I,4)=PHID(IHEM)+E-AMIN1(PHIDM0(IHEM)+WK1(I)*
     1    (PIH-PHIDM0(IHEM)),PHDMMX(IHEM))
        PHI(I,5)=PHID(IHEM)-E+AMIN1(PHIDP0(IHEM)+WK1(I)*
     1    (PIH-PHIDP0(IHEM)),PHDPMX(IHEM))
        PHI(I,6)=PHIN(IHEM)+E-AMIN1(PHINM0(IHEM)+WK1(I)*
     1    (PIH-PHINM0(IHEM)),PHNMMX(IHEM))
        PHI(I,7)=PHIN(IHEM)-E+AMIN1(PHINP0(IHEM)+WK1(I)*
     1    (PIH-PHINP0(IHEM)),PHNPMX(IHEM))
        PHI(I,1)=PHI(I,5)-PI2
        PHI(I,2)=PHI(I,6)-PI2
        PHI(I,3)=PHI(I,7)-PI2
        PHI(I,8)=PHI(I,4)+PI2
        PHIFUN(I)=0.
        PHIFN2(I) = 0.
        IFN(I) = merge(3,2,COLAT(I)-THETA0(IHEM)>=0.)
        IFLAG(I) = merge(IFN(I),IFLAG(I),IFLAG(I)-1==0.)
C       ****    ADD RING CURRENT ROTATION TO POTENTIAL (PHIRC)
        PHIRC = 0.
        WK2(I) = AMOD(ALON(I)+PHIRC+2.*PI2+PI, PI2) - PI
        WK3(I) = AMOD(ALON(I)+PHIRC+3.*PI2,PI2) - PI
    4 CONTINUE
C     ****     CALCULATE LONGITUDINAL VARIATION
      DO 5 N=1,7
C     DO 5 N=1,5
      DO 5 I=1,IMAXM
      PHIFUN(I)=PHIFUN(I)+.25*(PSI(N)+PSI(N+1)+(PSI(N)-
     1  PSI(N+1))*COS(AMOD(PI*(WK2(I)-PHI(I,N))/(PHI(I,N+1)-
     2  PHI(I,N)),PI2)))*(1.-SIGN(1.,(WK2(I)-PHI(I,N))*(WK2(I)-
     3  PHI(I,N+1))))
        PHIFN2(I)=PHIFN2(I)+.25*(PSI(N)+PSI(N+1)+(PSI(N)-
     1  PSI(N+1))*COS(AMOD(PI*(WK3(I)-PHI(I,N))/(PHI(I,N+1)-
     2  PHI(I,N)),PI2)))*(1.-SIGN(1.,(WK3(I)-PHI(I,N))*(WK3(I)-
     3  PHI(I,N+1))))
    5 CONTINUE
C     ****
C     ****     EVALUATE TOTAL POTENTIAL
C     ****
      DO 6 I=1,IMAXM
C     ****     OLD HEELIS MODEL (JGR, 1982)
C     POTEN(I)=PHIFUN(I)*merge(A2*((SIN(COLAT(I)+THETAC(IHEM))/SINTH0)**
C    1  R2(IHEM)-G0),(AMAX1(SIN(COLAT(I)),SINTH0)/SINTH0)**R1(IHEM),
C    2  IFLAG(I)-2==0.)
C     ****     NEW HEELIS MODEL (PERSONAL COMMUNICATION, 1987)
       POTEN(I)=merge( (2.*(PCEN(IHEM)-PHIFUN(I))+(PHIFUN(I)-PHIFN2(I))*
     1  0.75)*(COLAT(I)/THETA0(IHEM))**3 +
     2  (1.5*(PHIFUN(I)+PHIFN2(I))-3.*PCEN(IHEM))*(COLAT(I)/
     3  THETA0(IHEM))**2  +  0.75*(PHIFUN(I)-PHIFN2(I))*(COLAT(I)/
     4  THETA0(IHEM)) + PCEN(IHEM), PHIFUN(I)*(AMAX1(SIN(COLAT(I)),
     5  SINTH0)/SINTH0)**R1(IHEM)*EXP(7.*(1.-AMAX1(SIN(COLAT(I)),
     6  SINTHR1)/SINTHR1) ), IFLAG(I)-2==0.)
    6 CONTINUE
      RETURN
      END
C
