SUBROUTINE VDRIFT 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 "cons.h" include "consts.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 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 **** 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(-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 C SUBROUTINE VDRIFT2 implicit none C **** C **** Copy ion drift velocities for latitude J from arrays UI, C **** VI and WI in /DYNPHI/ to slots NUI, NVI and NWI in the C **** F-array. C **** include "params.h" include "fcom.h" include "buff.h" include "cons.h" include "consts.h" include "dynphi.h" include "fieldz.h" include "index.h" include "phys.h" real :: eex,eey,eez COMMON/VSCR/EEX(IMAXGP,ZKMXP), EEY(IMAXGP,ZKMXP), 1 EEZ(IMAXGP,ZKMXP) !$OMP THREADPRIVATE (/vscr/) !DIR$ TASKCOMMON vscr ! ! Local: integer :: nzk,k,i,nuik,nvik,nwik character(len=80) :: title NZK = NJ+NZ-1 DO K = 1,KMAXP1 NZK = NZK+1 ! write(6,"('vdrift2: j=',i3,' k=',i3,' f(:,nzk)=',/,(6e12.4))") ! | j,k,f(:,nzk) DO I = 1,IMAX+1 C **** C **** For latitude J, rotate EX and EY to geographic C **** orientation using Jacobian. Divide by distance from C **** center of earth. C **** C **** EEX = rotated EX C **** EEY = rotated EY C **** EEX(I,K) = (RJAC(I,J,1,1)*EX(I,J,K) + 1 RJAC(I,J,2,1)*EY(I,J,K))/(R00+F(I+2,NZK)) EEY(I,K) = (RJAC(I,J,1,2)*EX(I,J,K) + 1 RJAC(I,J,2,2)*EY(I,J,K))/(R00+F(I+2,NZK)) ENDDO ENDDO C **** C **** For K = 2,KMAX divide EZ by (Z(K+1) - Z(K-1)) C **** NZK = NJ+NZ DO K = 2,KMAX NZK = NZK+1 DO I = 1,IMAX+1 EEZ(I,K) = EZ(I,J,K)/(F(I+2,NZK+1)-F(I+2,NZK-1)) ENDDO ENDDO C **** C **** Extrapolate for values of EEZ at K = 1 and KMAXP1 C **** DO I = 1,IMAX+1 EEZ(I,1) = 2.*EEZ(I,2)-EEZ(I,3) EEZ(I,KMAXP1) = 2.*EEZ(I,KMAX)-EEZ(I,KMAX-1) ENDDO C **** C **** VI = (E X B/B**2) C **** Multiply by 1.E6 for results in m/sec C **** NUIK = NUI-1 NVIK = NVI-1 NWIK = NWI-1 DO K = 1,KMAXP1 NUIK = NUIK+1 NVIK = NVIK+1 NWIK = NWIK+1 DO I = 1,IMAX+1 C **** C **** ui = x-component of ion drift velocity C **** F(I+2,NUIK) = -(EEY(I,K)*ZZB(I,J)+EEZ(I,K)*XB(I,J))*1.E6/ 1 BMOD(I,J)**2 C **** C **** vi = y-component of ion drift velocity C **** F(I+2,NVIK) = (EEZ(I,K)*YB(I,J)+EEX(I,K)*ZZB(I,J))*1.E6/ 1 BMOD(I,J)**2 C **** C **** wi = Z-component of ion drift velocity C **** F(I+2,NWIK) = (EEX(I,K)*XB(I,J)-EEY(I,K)*YB(I,J))*1.E6/ 1 BMOD(I,J)**2 ENDDO ENDDO C **** C **** Periodic points C **** NUIK = NUI-1 NVIK = NVI-1 NWIK = NWI-1 DO I = 1,2 DO K = 1,KMAXP1 F(I,NUIK+K) = F(I+IMAX,NUIK+K) F(I,NVIK+K) = F(I+IMAX,NVIK+K) F(I,NWIK+K) = F(I+IMAX,NWIK+K) F(I+IMAXP2,NUIK+K) = F(I+2,NUIK+K) F(I+IMAXP2,NVIK+K) = F(I+2,NVIK+K) F(I+IMAXP2,NWIK+K) = F(I+2,NWIK+K) ENDDO ENDDO RETURN END C