
      subroutine mkionvel(fout,kmx,glat,glon,gmlat,gmlon,rdip,rdec,
     +  day,ut,u,v,w,xo2,xo,xn2,iefield,spv)
c
c Calculate ion drifts with or without efield (iefield flag), and return
c in fout(kmx,3) (ui,vi,wi):
c
      dimension fout(kmx,3),u(kmx),v(kmx),w(kmx),xo2(kmx),
     +  xo(kmx),xn2(kmx)
      data iseasav/0/, iutav/0/, bgaus/.4/, xmaso/32./, colfac/1.5/,
     +     xlag/0./
      save iseasav,iutav,bgaus,xmaxo,colfac,xlag
c
      cosd = cos(rdec)
      sind = sin(rdec)
      si = sin(rdip)
      ci = cos(rdip)
      if (iefield.gt.0) then
	timeut=ut+xlag
	call efield(gmlat,gmlon,day,ut+xlag,iseasav,iutav,pot,vu,ve)
        uit = ve
        vut = vu
      else
        uit = 0.
        vut = 0.
      endif
      do k=1,kmx
        if (xo(k).ne.spv.and.xn2(k).ne.spv.and.xo2(k).ne.spv.and.
     +      u(k).ne.spv.and.v(k).ne.spv.and.w(k).ne.spv) then
          vinp = 3.e-10 * (xo(k)+xn2(k)+xo2(k))
          rhoi = vinp / (9.57946e+3 * bgaus / xmaso)
          coi = 1. / (1. + rhoi*rhoi)
          unmag = u(k)*cosd - v(k)*sind
          vnmag = u(k)*sind + v(k)*cosd
          fout(k,1) = coi*(rhoi*rhoi*unmag+rhoi*(-vnmag*si-w(k)*ci)) +
     +                coi*(uit+rhoi*vut)
          fout(k,2) = coi*(rhoi*rhoi*vnmag+rhoi*unmag*si+
     +                vnmag*ci*ci-w(k)*si*ci)+coi*(vut*si-rhoi*uit*si)
          fout(k,3) = coi*(rhoi*rhoi*w(k)+rhoi*unmag*ci-vnmag*si*ci+
     +                w(k)*si*si)+coi*(vut*ci-rhoi*uit*ci)
        else
          fout(k,:) = spv
        endif
      enddo
      return
      end
