C SUBROUTINE VDRIFT C **** C **** Calculate 3-d electric field by numerical differentiation C **** of the 3-d dynamo potential in geomagnetic coordinates. C **** C **** Transform the three components of the field to geographic C **** space. C **** C **** Assemble needed common blocks: C **** include "params.h" PARAMETER (IMAXG = ZIMX, JMAXG = ZJMX, IMAXMP2 = IMAXM+2) PARAMETER (IMAXGP = IMAXG+1,JMAXGP = JMAXG+1, JMAXMP = JMAXM+1) include "cons.h" include "consts.h" include "dynphi.h" include "fieldz.h" C **** C **** Work space in / / C **** COMMON PHI3D(0:IMAXMP2,0:JMAXMP,ZKMXP), EMX(IMAXMP,JMAXM,ZKMXP), 1 EMY(IMAXMP,JMAXM,ZKMXP), EMZ(IMAXMP,JMAXM,ZKMXP), THETAM(IMAXG), 2 QSLOT(IMAXG), PSLOT(IMAXG), ISLOT(IMAXG), JSLOT(IMAXG) C **** C **** Copy PHIM3D to PHI3D adding extra points for C **** differentiation. C **** DO K = 1,KMAXP1 DO J = 1,JMAXM DO I = 1,IMAXMP PHI3D(I,J,K) = PHIM3D(I,J,K) ENDDO ENDDO ENDDO C **** C **** Add points corresponding to J = 0 and JMAX+1 performing C **** extrapolation across the poles. C **** IMAXMH = IMAXM/2 DO K = 1,KMAXP1 DO I = 1,IMAXMP PHI3D(I,0,K) = PHI3D(1+MOD(I-1+IMAXMH,IMAXM),2,K) PHI3D(I,JMAXM+1,K) = PHI3D(1+MOD(I-1+IMAXMH,IMAXM),JMAXM-1,K) ENDDO ENDDO C **** C **** Add periodic points to PHI3D C **** DO K = 1,KMAXP1 DO J = 0,JMAXM+1 PHI3D(0,J,K) = PHI3D(IMAXM,J,K) PHI3D(IMAXMP+1,J,K) = PHI3D(2,J,K) ENDDO ENDDO C CALL EZCNTR(PHIM3D(1,1,KMAXP),IMAXMP,JMAXM) C CALL MKCON(PHIM3D(1,1,KMAXP),IMAXMP,JMAXM,0.,0, C 1 'MLONG','MLAT','VDRIFT: PHIM3D') C **** C **** Calculate EMY, the meridional component of the electric C **** field, for all geomagnetic grid points. C **** C **** Note: The factor DT0DTS converts from the latitudinally C **** distorted geomagnetic grid to standard apex latitudes. C **** DO K = 1,KMAXP1 DO J = 1,JMAXM DO I = 1,IMAXMP EMY(I,J,K) = -(PHI3D(I,J+1,K)-PHI3D(I,J-1,K))/ 1 (2.*DLATM)*DT0DTS(J) ENDDO ENDDO ENDDO C **** C **** EMX = zonal electric field for all geomagnetic grid C **** points C **** C **** Note: The factor RCOS0S converts from the latitudinally C **** distortded geomagnetic grid to standard apex latitudes. C **** DO K = 1,KMAXP1 DO J = 2,JMAXM-1 CSTH0 = COS(-C(110)/2. + (J-1)*DLATM) DO I = 1,IMAXMP EMX(I,J,K) = -(PHI3D(I+1,J,K)-PHI3D(I-1,J,K))/ 1 (2.*DLONM*CSTH0)*RCOS0S(J) ENDDO ENDDO ENDDO C **** C **** Polar values of EMX C **** IMXM3Q = (3*IMAXM)/4 IMXM1Q = IMAXM/4 DO K = 1,KMAXP1 DO I = 1,IMAXMP EMX(I,1,K) = EMY(1+MOD(I-1+IMXM1Q,IMAXM),1,K) EMX(I,JMAXM,K) = EMY(1+MOD(I-1+IMXM3Q,IMAXM),JMAXM,K) ENDDO ENDDO C **** C **** Now calculate EMZ = d(PHI)/dz C **** DO K = 2,KMAX DO J = 1,JMAXM DO I = 1,IMAXMP EMZ(I,J,K) = -(PHIM3D(I,J,K+1)-PHI3D(I,J,K-1)) ENDDO ENDDO ENDDO C **** C **** Transform EMX, EMY, EMZ to geographic space. C **** C **** For each geographic grid point, determine the C **** corresponding theta0 in the distorted geomagnetic C **** latitude grid. C **** C **** ALATM(I,J) is the thetas (standard apex latitude) C **** corresponding to geographic grid point (I,J). Find the C **** slot frr this value in the table of theta0 vs.thetas, C **** which has 1 degree spacing from 0. to 90. degrees. C **** DO J = 1,JMAXG DO I = 1,IMAXG PSLOT(I) = ABS(ALATM(I,J))*180./C(110)+1. ISLOT(I) = PSLOT(I) PSLOT(I) = PSLOT(I)-FLOAT(ISLOT(I)) C **** C **** Now interpolate for theta0 in table using linear C **** interpolation. C **** THETAM(I) = SIGN((1.-PSLOT(I))*TABLE(ISLOT(I),2) + 1 PSLOT(I)*TABLE(ISLOT(I)+1,2),ALATM(I,J)) C **** C **** Locate magnetic grid element containing geographic C **** grid point (I,J) C **** PSLOT(I) = (ALONM(I,J)+C(110))/DLONM + 1. ISLOT(I) = PSLOT(I) PSLOT(I) = PSLOT(I)-FLOAT(ISLOT(I)) QSLOT(I) = (THETAM(I)+C(110)/2.)/DLATM + 1. JSLOT(I) = QSLOT(I) QSLOT(I) = QSLOT(I)-FLOAT(JSLOT(I)) ENDDO DO K = 1,KMAXP1 DO I = 1,IMAXG C **** C **** Transform EMX, EMY, EMZ to geographic frame as C **** EX, EY and EZ C **** EX(I,J,K) = (1.-PSLOT(I))*(1.-QSLOT(I)) * 1 EMX(ISLOT(I),JSLOT(I),K) + 2 PSLOT(I)*(1.-QSLOT(I)) * EMX(ISLOT(I)+1,JSLOT(I),K) + 3 PSLOT(I)*QSLOT(I) * EMX(ISLOT(I)+1,JSLOT(I)+1,K) + 4 (1.-PSLOT(I))*QSLOT(I) * EMX(ISLOT(I),JSLOT(I)+1,K) EY(I,J,K) = (1.-PSLOT(I))*(1.-QSLOT(I)) * 1 EMY(ISLOT(I),JSLOT(I),K) + 2 PSLOT(I)*(1.-QSLOT(I)) * EMY(ISLOT(I)+1,JSLOT(I),K) + 3 PSLOT(I)*QSLOT(I) * EMY(ISLOT(I)+1,JSLOT(I)+1,K) + 4 (1.-PSLOT(I))*QSLOT(I) * EMY(ISLOT(I),JSLOT(I)+1,K) ENDDO EX(IMAXG+1,J,K) = EX(1,J,K) EY(IMAXG+1,J,K) = EY(1,J,K) ENDDO DO K = 2,KMAX DO I = 1,IMAXG EZ(I,J,K) = (1.-PSLOT(I))*(1.-QSLOT(I)) * 1 EMZ(ISLOT(I),JSLOT(I),K) + 2 PSLOT(I)*(1.-QSLOT(I)) * EMZ(ISLOT(I)+1,JSLOT(I),K) + 3 PSLOT(I)*QSLOT(I) * EMZ(ISLOT(I)+1,JSLOT(I)+1,K) + 4 (1.-PSLOT(I))*QSLOT(I) * EMZ(ISLOT(I),JSLOT(I)+1,K) ENDDO EZ(IMAXG+1,J,K) = EZ(1,J,K) ENDDO ENDDO RETURN END