#include "dims.h" ! ! am_02/02: calculate the 2D and 3D electric drift velocities e_d1 and e_d2 ! C SUBROUTINE THREED use cons_module,only: kmaxp1,pi 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 "consdyn.h" #include "dynphi.h" ! ! 4/00: make these local (s.a., vdrift.f): real :: thetam(imaxmp),pslot(imaxmp),qslot(imaxmp) integer :: islot(imaxmp),jslot(imaxmp) ! #include "cterp.h" ! ! Global: common/ns3d/ sigma1m3d(imaxmp,jmaxm,-2:zkmxp), | sigma2m3d(imaxmp,jmaxm,-2:zkmxp), | bmodm3d(imaxmp,jmaxm), | adotv1m3d(imaxmp,jmaxm,-2:zkmxp), | adotv2m3d(imaxmp,jmaxm,-2:zkmxp), | a1a2m3d(imaxmp,jmaxm), | adotam3d(imaxmp,jmaxm), | sinim3d(imaxmp,jmaxm), | ed13d(imaxmp,jmaxm,-2:zkmxp),ed23d(imaxmp,jmaxm,-2:zkmxp) real :: sigma1m3d,sigma2m3d,bmodm3d,adotv1m3d,adotv2m3d, | a1a2m3d,adotam3d,sinim3d,ed13d,ed23d ! Local: integer :: j,k,i,ip1f,ip2f,ip3f real :: cosltm,phims,phimn,ctm character(len=80) :: title ! External: real,external :: sddot ! in util.F !ADRb real :: csth0, sym real :: ed1(imaxmp,jmaxm),ed2(imaxmp,jmaxm) ! C Calculate Ed1, Ed2 components of electric field DO J = 2,JMAXM-1 CSTH0 = COS(-pi/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)/(R0*1.e-2) ed2(I,JMAXMH+1) = (4.*PHIM(I,JMAXMH+2)-PHIM(I,JMAXMH+3) 1 -3.*PHIM(I,JMAXMH+1))/(2.*DLATM)/(R0*1.e-2) ed2(I,JMAXMH-1) = (4.*PHIM(I,JMAXMH-2)-PHIM(I,JMAXMH-3) 1 -3.*PHIM(I,JMAXMH-1))/(2.*DLATM)/(R0*1.e-2) ENDDO ! I ! 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)) sym = 1. if (J.lt.JMAXMH) sym = -1. DO 1 I = 1,IMAXM C **** THETAM(I) = THETAS FOR EACH VALUE OF I c Correction ctm = SQRT(R0/(RE+ZZM(I,J,K)))*COSLTM THETAM(I) = ACOS(ctm) 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)) ! Make THETAM negative for southern hemisphere 2 *sym 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)) 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)) 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 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) 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) ed13d(IMAXMP,J,K) = ed13d(1,J,K) ed23d(IMAXMP,J,K) = ed23d(1,J,K) 2 CONTINUE ! ! Save magnetic potential on secondary histories: 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 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