#include "dims.h" C SUBROUTINE VDRIFT use cons_module,only: kmax,kmaxp1,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 **** ! real :: phi3d,emx,emy,emz,thetam,qslot,pslot ! integer :: islot,jslot ! 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) ! ! 4/00: make these local (some locals of same name also in threed) real :: phi3d(0:imaxmp2,0:jmaxmp,zkmxp), emx(imaxmp,jmaxm,zkmxp), | emy(imaxmp,jmaxm,zkmxp), emz(imaxmp,jmaxm,zkmxp), thetam(imaxg), | qslot(imaxg), pslot(imaxg) integer :: islot(imaxg),jslot(imaxg) ! ! Local: integer :: i,j,k,imaxmh,imxm3q,imxm1q real :: csth0 real,dimension(zimxp,zkmxp) :: extmp,eytmp,eztmp C **** C **** Copy PHIM3D to PHI3D adding extra points for C **** differentiation. C **** ! do j=1,jmaxm ! do k=1,zkmxp,zkmxp-1 ! write(6,"('vdrift: j=',i3,' k=',i3,' phim3d(:,j,k)=',/, ! | (6e12.4))") j,k,phim3d(:,j,k) ! enddo ! enddo 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 **** ! write(6,"('efield: dlatm=',e12.4,' dt0dts(nmlat)=',/,(6e12.4))") ! | dlatm,dt0dts 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) ! write(6,"('efield: k=',i3,' j=',i3,' i=',i3, ! | ' phi3d(i,j+1,k)=',e12.4,' phi3d(i,j-1,k)=',e12.4, ! | ' emy=',e12.4)") ! | k,j,i,phi3d(i,j+1,k),phi3d(i,j-1,k),emy(i,j,k) 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(-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) ! write(6,"('efield: k=',i3,' j=',i3,' i=',i3, ! | ' csth0=',e12.4,' phi3d(i+1,j,k)=',e12.4, ! | ' phi3d(i-1,j,k)=',e12.4,' emx=',e12.4)") ! | k,j,i,csth0,phi3d(i+1,j,k),phi3d(i-1,j,k),emx(i,j,k) 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) ! write(6,"('efield: k=',i3,' i=',i3,' emx(i,1,k)=',e12.4, ! | ' emx(i,nmlat,k)=',e12.4)") k,i,emx(i,1,k),emx(i,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) ! write(6,"('efield: j=',i2,' k=',i2,' i=',i2,' p,qslot=', ! | 2e12.4,' i,jslot=',2i3,' ex(i,j,k)=',e12.4)") ! | j,k,i,pslot(i),qslot(i),islot(i),jslot(i),ex(i,j,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) ! write(6,"('efield: j=',i2,' k=',i2,' i=',i2,' p,qslot=', ! | 2e12.4,' i,jslot=',2i3,' ey(i,j,k)=',e12.4)") ! | j,k,i,pslot(i),qslot(i),islot(i),jslot(i),ey(i,j,k) ! print *,'ey=',ey(i,j,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 ! ! Save e-field on geographic grid to secondary history. ! dynphi.h: EX(IMAXGP,JMAXG,ZKMXP), EY(IMAXGP,JMAXG,ZKMXP), ! local: real,dimension(zimxp,zkmxp) :: extmp,eytmp,eztmp ! do j=1,zjmx ! zjmx==jmaxg extmp(3:zimx+2,:) = ex(1:imaxg,j,:) ! 3:74 <- 1:72 eytmp(3:zimx+2,:) = ey(1:imaxg,j,:) eztmp(3:zimx+2,:) = ez(1:imaxg,j,:) do i=1,2 extmp(i,:) = extmp(zimx+i,:) ! 1:2 <- 73:74 eytmp(i,:) = eytmp(zimx+i,:) eztmp(i,:) = eztmp(zimx+i,:) extmp(zimx+2+i,:) = extmp(i+2,:) ! 75:76 <- 3:4 eytmp(zimx+2+i,:) = eytmp(i+2,:) eztmp(zimx+2+i,:) = eztmp(i+2,:) enddo ! call addfsech('EX',' ',' ',extmp,zimxp,zkmxp,zkmx,j) ! call addfsech('EY',' ',' ',eytmp,zimxp,zkmxp,zkmx,j) ! call addfsech('EZ',' ',' ',eztmp,zimxp,zkmxp,zkmx,j) enddo RETURN END