#include "dims.h" SUBROUTINE FLWV32(POTEN) use cons_module,only: imax implicit none C **** C **** COMPUTE ELECTRIC POTENTIAL C **** #include "params.h" ! ! Args: real,intent(out) :: POTEN(IMAXMP) C **** C **** WORK SPACE C **** ! iflag, dlat, dlon, and ratio were set in potm. ! integer :: iflag,ifn real :: dlat,dlon,ratio,psi,phi,phifun,colat,alon,sinlat, | coslat,sinlon,coslon,wk1,wk2,wk3 COMMON/flwv32_com/IFLAG(IMAXM), DLAT(IMAXM), DLON(IMAXM), + RATIO(IMAXM), PSI(8), PHI(IMAXM,8), IFN(IMAXM), PHIFUN(IMAXM), + COLAT(IMAXM), ALON(IMAXM), SINLAT(IMAXM), + COSLAT(IMAXM), SINLON(IMAXM), COSLON(IMAXM), WK1(IMAXM), + WK2(IMAXM),WK3(IMAXM) #include "ioncr.h" real :: cosofa,sinofa,aslona,pi,pi2 COMMON /OVALPOS/ COSOFA(2),SINOFA(2),ASLONA(2),PI,PI2 ! ! Local: real :: PHIFN2(IMAXM) real :: PHDPMX(2),PHDMMX(2),PHNPMX(2),PHNMMX(2), 1 COSOFC(2),SINOFC(2),ASLONC(2) integer IHEM real OFDA,ofdc,sinth0,sinthr1,phirc,pih ! DATA E/1.E-6/,PI/3.1415927/,PI2/6.2831853/,PIH/1.5707963/ real,parameter :: E=1.E-6 integer :: n,imaxd2,i C **** C **** CONVERT PARAMETERS TO RADIANS ETC. C **** ! ! Once per iteration init. Istar is in reg common (ioncr.h) ! (is reset to 0 by adheel before this routine) ! (note istar is used for once-per-run init in heelis.f) ! PI=4.*ATAN(1.) PI2=2.*PI PIH=.5*PI DO 50 N=1,2 OFDA = SQRT(OFFA(N)**2 + DSKOFA(N)**2) COSOFA(N) = COS(OFDA) SINOFA(N) = SIN(OFDA) ASLONA(N) = ASIN(DSKOFA(N)/OFDA) OFDC = SQRT(OFFC(N)**2 + DSKOFC(N)**2) COSOFC(N) = COS(OFDC) SINOFC(N) = SIN(OFDC) ASLONC(N) = ASIN(DSKOFC(N)/OFDC) 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) PHDMMX(N)=PHNPMX(N) 50 CONTINUE 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.) ifn(i) = 3 if (colat(i) < theta0(ihem)) ifn(i) = 2 ! IFLAG(I) = merge(IFN(I),IFLAG(I),IFLAG(I)-1 == 0) if (iflag(i)==1) iflag(i) = ifn(i) 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)**RR1(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)**RR1(IHEM)*EXP(7.*(1.-AMAX1(SIN(COLAT(I)), ! 6 SINTHR1)/SINTHR1) ), IFLAG(I)-2 == 0) ! 6 CONTINUE ! DO 6 I=1,IMAXM if (iflag(i)==2) then POTEN(I) = (2.*(PCEN(IHEM)-PHIFUN(I))+(PHIFUN(I)-PHIFN2(I))* | 0.75)*(COLAT(I)/THETA0(IHEM))**3 + | (1.5*(PHIFUN(I)+PHIFN2(I))-3.*PCEN(IHEM))*(COLAT(I)/ | THETA0(IHEM))**2 + 0.75*(PHIFUN(I)-PHIFN2(I))*(COLAT(I)/ | THETA0(IHEM)) + PCEN(IHEM) else POTEN(I) = PHIFUN(I)*(AMAX1(SIN(COLAT(I)), | SINTH0)/SINTH0)**RR1(IHEM)*EXP(7.*(1.-AMAX1(SIN(COLAT(I)), | SINTHR1)/SINTHR1)) endif 6 CONTINUE RETURN END C