#include "dims.h"
C
      SUBROUTINE VDRIFT2
      use cons_module,only: kmax,kmaxp1,imax,imaxp2
      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 "consdyn.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

!     call addfsech('UI_IONV',' ',' ',f(1,nui),zimxp,zkmxp,zkmxp,j)
!     call addfsech('VI_IONV',' ',' ',f(1,nvi),zimxp,zkmxp,zkmxp,j)
!     call addfsech('WI_IONV',' ',' ',f(1,nwi),zimxp,zkmxp,zkmxp,j)

      RETURN
      END
C