#include "dims.h" C SUBROUTINE VDRIFT use cons_module,only: kmax,kmaxp1,imax,pi implicit none 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" integer,parameter :: imaxmp2=imaxm+2, jmaxmp=jmaxm+1 #include "consdyn.h" #include "dynphi.h" #include "fieldz.h" C **** C **** Work space in / / C **** ! 4/98: changed this from blank common to local space to avoid ! conflict w/ F array, which is also blank common. ! ! 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) ! real :: 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) integer :: ISLOT(IMAXG), JSLOT(IMAXG) integer :: i,j,k,imaxmh,imxm3q,imxm1q real :: csth0 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 **** 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 **** distorted geomagnetic grid to standard apex latitudes. C **** DO K = 1,KMAXP1 DO J = 2,JMAXM-1 CSTH0 = COS(-pi/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./pi+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)+pi)/DLONM + 1. ISLOT(I) = PSLOT(I) PSLOT(I) = PSLOT(I)-FLOAT(ISLOT(I)) QSLOT(I) = (THETAM(I)+pi/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