#include "dims.h" SUBROUTINE FLWV32(j,POTEN) use cons_module,only: imax implicit none C **** C **** COMPUTE ELECTRIC POTENTIAL C **** #include "params.h" 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 ! ! Args: integer,intent(in) :: j real,intent(out) :: POTEN(IMAXMP) ! ! Local: integer IHEM,n,imaxd2,i real :: PHIFN2(IMAXM),ofda,e,pih,ofdc,sinth0,sinthr1,phirc real :: 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/ DATA E/1.E-6/ 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) ! ! 8/15/00: Rather than sorting out local,save vs local vs common, ! just go ahead and do this short init at every latitude: ! ! if (istar==0) then ! istar = 1 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 ! 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. if (COLAT(I)-THETA0(IHEM) >= 0.) then ifn(i) = 3 else ifn(i) = 2 endif 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 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 ! write(6,"(/'flwv32: j=',i2,' ihem=',i2)") j,ihem ! write(6,"(' theta0(ihem)=',e12.4,' pcen(ihem)=',e12.4, ! | ' rr1(ihem)=',e12.4)") theta0(ihem),pcen(ihem),rr1(ihem) ! write(6,"(' sinth0=',e12.4,' sinthr1=',e12.4)") sinth0,sinthr1 ! write(6,"(' phifun=',/(6e12.4))") phifun ! write(6,"(' phifn2=',/(6e12.4))") phifn2 ! write(6,"(' poten =',/(6e12.4))") poten RETURN END C