C SUBROUTINE THREED implicit none C **** C **** CALCULATE 3-D DYNAMO POTENTIAL IN GEOMAGNETIC COORDINATES C **** C **** PHIM(IMAXMP,JMAXM) = 2-D POTENTIAL IN GEOMAGNETIC C **** COORDINATES OUTPUT BY DYNAMO CODE. C **** C **** PHIM3D(IMAXMP,JMAXM,-2:ZKMXP) = 3-D DYNAMO POTENTIAL IN C **** GEOMAGNETIC COORDINATES GENERATED BY THIS SUBROUTINE. C **** include "params.h" ! integer,parameter :: kmax=zkmx, kmaxp=kmax+1 include "consts.h" include "dynphi.h" ! real :: thetam,pslot,qslot ! integer :: islot,jslot ! COMMON THETAM(IMAXMP),ISLOT(IMAXMP),JSLOT(IMAXMP), ! 1 PSLOT(IMAXMP),QSLOT(IMAXMP) ! ! 4/00: make these local (s.a., vdrift.f): real :: thetam(imaxmp),pslot(imaxmp),qslot(imaxmp) integer :: islot(imaxmp),jslot(imaxmp) ! include "cterp.h" ! ! Local: integer :: j,k,i,ip1f,ip2f,ip3f real :: cosltm,phims,phimn character(len=80) :: title ! External: real,external :: sddot ! in util.F !ADRb include "cons.h" ! need for C(110) (pi) real :: csth0, sym real :: ed1(IMAXMP,JMAXM),ed2(IMAXMP,JMAXM), | ed13d(imaxmp,jmaxm,-2:zkmxp), ed23d(imaxmp,jmaxm,-2:zkmxp) C Calculate Ed1, Ed2 components of electric field DO J = 2,JMAXM-1 CSTH0 = COS(-C(110)/2. + (J-1)*DLATM) DO I = 2,IMAXMP-1 ed1(I,J) = -(PHIM(I+1,J)-PHIM(I-1,J))/ 1 (2.*DLONM*CSTH0)*RCOS0S(J)/(R0*1.e-2) ENDDO ed1(1,J) = -(PHIM(2,J)-PHIM(IMAXMP-1,J))/ 1 (2.*DLONM*CSTH0)*RCOS0S(J)/(R0*1.e-2) ed1(IMAXMP,J) = ed1(1,J) ENDDO ! J ! Southern hemisphere DO J = 2,JMAXMH-1 DO I = 1,IMAXMP ed2(I,J) = -(PHIM(I,J+1)-PHIM(I,J-1))/ 1 (2.*DLATM)*DT1DTS(J)/(R0*1.e-2) ENDDO ENDDO ! J ! Northern hemisphere DO J = JMAXMH+1,JMAXM-1 DO I = 1,IMAXMP ed2(I,J) = (PHIM(I,J+1)-PHIM(I,J-1))/ 1 (2.*DLATM)*DT1DTS(J)/(R0*1.e-2) ENDDO ENDDO ! J DO I = 1,IMAXMP ! Poles ip1f = i + IMAXM/4 if (ip1f.gt.IMAXMP) ip1f = ip1f - IMAXM ip2f = i + IMAXM/2 if (ip2f.gt.IMAXMP) ip2f = ip2f - IMAXM ip3f = i + 3*IMAXM/4 if (ip3f.gt.IMAXMP) ip3f = ip3f - IMAXM ed1(I,1) = .25*(ed1(I,2) - ed1(ip2f,2) + 1 ed2(ip1f,2) - ed2(ip3f,2)) ed1(I,JMAXM) = .25*(ed1(I,JMAXM-1) - ed1(ip2f,JMAXM-1) + 1 ed2(ip1f,JMAXM-1) - ed2(ip3f,JMAXM-1)) ed2(I,1) = .25*(ed2(I,2) - ed2(ip2f,2) - 1 ed1(ip1f,2) + ed1(ip3f,2)) ed2(I,JMAXM) = .25*(ed2(I,JMAXM-1) - ed2(ip2f,JMAXM-1) - 1 ed1(ip1f,JMAXM-1) + ed1(ip3f,JMAXM-1)) ! Equator ed2(I,JMAXMH) = (4.*PHIM(I,JMAXMH+1)-PHIM(I,JMAXMH+2) 1 -3.*PHIM(I,JMAXMH))/(2.*DLATM)*DT1DTS(JMAXMH)/(R0*1.e-2) ENDDO ! I !ADRe C **** C **** FOR EACH GEOMAGNETIC GRID POINT OVER THE 3-D GRID, C **** DETERMINE THE LOCATION OF THE FOOT OF THE MAGNETIC C **** FIELD PASSING THROUGH THE GRID POINT. ASSUME THAT C **** THE FIELD LINES ARE DIPOLE IN SHAPE. RECALLING THAT C **** THE MODEL ASSUMES CONSTANT ELECTRIC POTENTIAL ALONG C **** EACH FIELD LINE, INTERPOLATION IN THE 2-D ARRAY PHIM C **** GIVES THE POTENTIAL AT THE ORIGINAL GRID POINT. C **** C **** LOOP OVER GEOMAGNETIC GRID. C **** ! kmaxp1 = zkmxp DO 1 K = -2,kmaxp1 DO 1 J = 2,JMAXM-1 C **** COSLTM = COS(THETAS) COSLTM = COS(YLATM(J)) !ADRb sym = 1. if (J.lt.JMAXMH) sym = -1. !ADRe DO 1 I = 1,IMAXM C **** THETAM(I) = THETAS FOR EACH VALUE OF I THETAM(I) = ACOS(SQRT(R0/(R00+ZZM(I,J,K)))*COSLTM) C **** C **** FIND POSITION OF THETAM IN TABLE OF C **** THETA0 VS. THETAS C **** PSLOT(I) = THETAM(I)*180./PI+1. ISLOT(I) = PSLOT(I) PSLOT(I) = PSLOT(I)-FLOAT(ISLOT(I)) C **** C **** DETERMINE THATA0 BY INTERPOLATION IN TABLE C **** THETAM(I) = (1.-PSLOT(I))*TABLE(ISLOT(I),2)+ PSLOT(I)* 1 TABLE(ISLOT(I)+1,2) !ADRb ! Make THETAM negative for southern hemisphere 2 *sym !ADRe C **** C **** NOW IDENTIFY GEOMAGNETIC GRID ELEMENT IN WHICH C **** THE POINT LIES IN PHIM C **** ISLOT(I) = I PSLOT(I) = 0. QSLOT(I) = (THETAM(I)+PI/2.)/DLATM+1. JSLOT(I) = QSLOT(I) QSLOT(I) = QSLOT(I)-FLOAT(JSLOT(I)) C **** C **** INTERPOLATE IN PHIM USING BILINEAR INTERPOLATION C **** PHIM3D(I,J,K) = (1.-QSLOT(I))*((1.-PSLOT(I))* 1 PHIM(ISLOT(I),JSLOT(I))+PSLOT(I)* 2 PHIM(ISLOT(I)+1,JSLOT(I)))+QSLOT(I)*((1.-PSLOT(I))* 3 PHIM(ISLOT(I),JSLOT(I)+1)+PSLOT(I)* 4 PHIM(ISLOT(I)+1,JSLOT(I)+1)) !ADRb ed13d(I,J,K) = (1.-QSLOT(I))*((1.-PSLOT(I))* 1 ed1(ISLOT(I),JSLOT(I))+PSLOT(I)* 2 ed1(ISLOT(I)+1,JSLOT(I)))+QSLOT(I)*((1.-PSLOT(I))* 3 ed1(ISLOT(I),JSLOT(I)+1)+PSLOT(I)* 4 ed1(ISLOT(I)+1,JSLOT(I)+1)) ed23d(I,J,K) = (1.-QSLOT(I))*((1.-PSLOT(I))* 1 ed2(ISLOT(I),JSLOT(I))+PSLOT(I)* 2 ed2(ISLOT(I)+1,JSLOT(I)))+QSLOT(I)*((1.-PSLOT(I))* 3 ed2(ISLOT(I),JSLOT(I)+1)+PSLOT(I)* 4 ed2(ISLOT(I)+1,JSLOT(I)+1)) !ADRe 1 CONTINUE C **** C **** CALCULATE VALUE OF PHIM3D AT GEOMAGNETIC POLES C **** PHIMS = sddot(IMAXM,UNIT,PHIM(1,1))/FLOAT(IMAXM) PHIMN = sddot(IMAXM,UNIT,PHIM(1,JMAXM))/FLOAT(IMAXM) DO 4 K = -2,kmaxp1 DO 4 I = 1,IMAXM PHIM3D(I,1,K) = PHIMS PHIM3D(I,JMAXM,K) = PHIMN !ADRb ed13d(I,1,K) = ed1(I,1) ed13d(I,JMAXM,K) = ed1(I,JMAXM) ed23d(I,1,K) = ed2(I,1) ed23d(I,JMAXM,K) = ed2(I,JMAXM) !ADRe 4 CONTINUE C **** C **** PERIODIC POINT C **** DO 2 K = -2,kmaxp1 DO 2 J = 1,JMAXM PHIM3D(IMAXMP,J,K) = PHIM3D(1,J,K) !ADRb ed13d(IMAXMP,J,K) = ed13d(1,J,K) ed23d(IMAXMP,J,K) = ed23d(1,J,K) !ADRe 2 CONTINUE ! do j=1,jmaxm ! do k=1,zkmxp,zkmxp-1 ! write(6,"('threed: j=',i3,' k=',i3,' phim3d(:,j,k)=',/, ! | (6e12.4))") j,k,phim3d(:,j,k) ! enddo ! enddo ! ! Save magnetic potential on secondary histories: !ADRb do j=1,jmaxm call addfsech('PHIM3D','PHIM3D','Volts', | phim3d(:,j,:),imaxmp,kmaxp1+3,kmaxp1+3,j) call addfsech('ED1M3D','ED1M3D','V/m', | ed13d(:,j,:),imaxmp,kmaxp1+3,kmaxp1+3,j) call addfsech('ED2M3D','ED2M3D','V/m', | ed23d(:,j,:),imaxmp,kmaxp1+3,kmaxp1+3,j) enddo !ADRe C **** C **** TRANSFORM PHIM3D TO GEOGRAPHIC COORDINATES IN C **** DYNPOT(IMAXGP,0:JMAXGP,kmaxp1) C **** DO 3 K = 1,kmaxp1 CALL GRDTRP(PHIM3D(1,1,K),DYNPOT(1,0,K),IM(1,0),JM(1,0), 1 DIM(1,0),DJM(1,0),IMAXGP,IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM) 3 CONTINUE C **** C **** PERIODIC POINTS C **** DO 5 K = 1,kmaxp1 DO 5 J = 0,JMAXG+1 DYNPOT(IMAXGP,J,K) = DYNPOT(1,J,K) 5 CONTINUE C **** C **** TRANSFORM PHIHM.TO GEOGRAPHIC COORDINATES IN C **** PHIH(IMAXGP,0:JMAXGP) C **** CALL GRDTRP(PHIHM(1,1),PHIH(1,0),IM(1,0),JM(1,0), 1 DIM(1,0),DJM(1,0),IMAXGP,IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM) C **** C **** PERIODIC POINTS C **** DO 6 J = 0,JMAXG+1 PHIH(IMAXGP,J) = PHIH(1,J) 6 CONTINUE RETURN END C