#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
      real,dimension(zimxp,zkmxp) :: extmp,eytmp,eztmp

!     write(6,"('ionvel: r00=',e12.4)") r00

      extmp(3:zimx+2,:) = ex(1:zimx,j,:) ! 3:74 <- 1:72
      eytmp(3:zimx+2,:) = ey(1:zimx,j,:)
      eztmp(3:zimx+2,:) = ez(1:zimx,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_IONV',' ',' ',extmp,zimxp,zkmxp,zkmxp,j)
!     call addfsech('EY_IONV',' ',' ',eytmp,zimxp,zkmxp,zkmxp,j)
!     call addfsech('EZ_IONV',' ',' ',eztmp,zimxp,zkmxp,zkmxp,j)

      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))

!         write(6,"('lat=',i2,' k=',i2,' i=',i2,' rjac=',2e12.4,
!    |      ' ex=',e12.4,' ey=',e12.4,' z=',e12.4,' eex=',e12.4)")
!    |      j,k,i,rjac(i,j,1,1),rjac(i,j,2,1),ex(i,j,k),
!    |      ey(i,j,k),f(i+2,nzk),eex(i,k)

          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     ****

! these are ok
!     call addfsech('Z_IONV',' ',' ',f(1,nj+nz),zimxp,zkmxp,zkmxp,j)

      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))

!         write(6,"('ionvel: lat=',i2,' k=',i2,' i=',i2,
!    |      ' ez(i,lat,k)=',e14.7,' z(k+1)=',e14.7,' z(k-1)=', 
!    |      e14.7)") j,k,i,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

      extmp(3:zimx+2,:) = eex(1:zimx,:) ! 3:74 <- 1:72
      eytmp(3:zimx+2,:) = eey(1:zimx,:)
      eztmp(3:zimx+2,:) = eez(1:zimx,:)
      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('EEX',' ',' ',extmp,zimxp,zkmxp,zkmxp,j)
!     call addfsech('EEY',' ',' ',eytmp,zimxp,zkmxp,zkmxp,j)
!     call addfsech('EEZ',' ',' ',eztmp,zimxp,zkmxp,zkmxp,j)

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

!         write(6,"('ionvel: lat=',i2,' k=',i2,' i=',i2, 
!    |      ' eey=',e14.7,' eez=',e14.7)") j,k,i,eey(i,k),eez(i,k)
!         write(6,"('  xb=',e14.7,' zb=',e14.7)")
!    |      xb(i,j),zzb(i,j) 
!         write(6,"('  bmod=',e14.7,' ui=',e14.7)") 
!    |      bmod(i,j),f(i+2,nuik)

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

!     call addfsech('UI_VEL',' ',' ',f(1,nui),zimxp,zkmxp,zkmxp,j)
!     call addfsech('VI_VEL',' ',' ',f(1,nvi),zimxp,zkmxp,zkmxp,j)
!     call addfsech('WI_VEL',' ',' ',f(1,nwi),zimxp,zkmxp,zkmxp,j)

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_VEL',' ',' ',f(1,nui),zimxp,zkmxp,zkmxp,j)
!     call addfsech('VI_VEL',' ',' ',f(1,nvi),zimxp,zkmxp,zkmxp,j)
!     call addfsech('WI_VEL',' ',' ',f(1,nwi),zimxp,zkmxp,zkmxp,j)

      RETURN
      END
C