c
c------------------------------------------------------------------
c Begin file /home/sting/foster/proclib/mkuivis.f
c------------------------------------------------------------------
c
      subroutine mkuivis(ixk,sui,svi,gcmht,lu)
c
c Find ion drifts from sum of dynamo and heelis potentials.
c Return ui,vi drifts at zp index ispls(ixk)
c
      include 'tgcmparam.h'
      include 'input.h'
      include 'selgrid.h'
      include 'tiegcmfld.h'
      include 'mkuivi.h'
c
      parameter(imxm1=imx-1,jmxp1=jmx+1,jmxm1=jmx-1,imxp2=imx+2)
      parameter(rs=6.475165e+8)
      dimension sui(imx,jmx),svi(imx,jmx)
      dimension phi(imx,kmx,jmx),gcmht(imx,kmx,jmx)
      data scale/1.e6/
      data magrd/0/
c
      if (ixk.le.0.or.ixk.gt.npls) then
        write(6,"('>>> mkuivis: bad ixk=',i5)") ixk 
        ixk = -1
        return
      endif
c
c Find sum of dynamo and heelis potentials:
c
      do k=1,kmx
        phi(:,k,:) = poten(:,k,:) + hpoten(:,:)
      enddo
c
c Define work array with wrap around points:
c
      do 100 k=1,kmx
        do 105 j=1,jmx
        do 105 i=1,imxm1
          wrk(i+2,k,j+1) = phi(i,k,j) 
 105    continue
c
c Over the poles:
c
        do 110 i=1,imxm1
          wrk(i+2,k,1) = phi(1+mod(i+jmxm1,imxm1),k,1)
          wrk(i+2,k,jmxp2) = phi(1+mod(i+jmxm1,imxm1),k,jmx)
 110    continue
c
c Periodic points:
c
        do 115 j=1,jmxp2
        do 115 i=1,2
          wrk(i,k,j) = wrk(i+imxm1,k,j)
          wrk(i+imxp1,k,j) = wrk(i+2,k,j)
 115    continue
 100  continue
c
c Read magnetic coords:
c
      if (magrd.le.0) then
        rewind lu
        read(lu) bx,by,bz
        magrd = 1
      endif
c
c Grid loop:
c
      do 200 j=1,jmx
      do 200 i=1,imx
        bxyzsum = bx(i,j)**2. + by(i,j)**2. + bz(i,j)**2.
        b(i,j) = bxyzsum**0.5
        do 205 k=1,kmx
c
c Calculation of electric field:
c
          ex(k) = -(wrk(i+3,k,j+1) - wrk(i+1,k,j+1)) /
     +             (2.*dlon*dtr*(rs+gcmht(i,k,j)*1.e+5)*
     +              cos(gcmlat(j)*dtr))
          ey(k) = -(wrk(i+2,k,j+2) - wrk(i+2,k,j)) /
     +             (2.*dlat*dtr*(rs+gcmht(i,k,j)*1.e+5))
 205    continue
        do 210 k=2,kmx-1
          ez(k) = -((phi(i,k+1,j)-phi(i,k,j)) /
     +              ((gcmht(i,k+1,j)-gcmht(i,k,j))*1.e+5) +
     +              (phi(i,k,j)-phi(i,k-1,j)) /
     +              ((gcmht(i,k,j)-gcmht(i,k-1,j))*1.e+5)) / 2.
 210    continue
        ez(1) = ez(2)*2.-ez(3)
        ez(kmx) = ez(kmx-1)*2.-ez(kmx-2)
c
c Find the ion drifts (m/s) at selected pressure level:
c
        k = ixk
        sui(i,j) = (ey(ispls(k))*bz(i,j) - ez(ispls(k))*by(i,j)) * 
     +            scale / b(i,j)**2.
        svi(i,j) = (ez(ispls(k))*bx(i,j) - ex(ispls(k))*bz(i,j)) * 
     +            scale / b(i,j)**2.
 200  continue
c     
      return
      end