c
c------------------------------------------------------------------
c Begin file /home/sting/foster/tigcm/mkuivih.f
c------------------------------------------------------------------
c
      subroutine mkuivih(phi,gcmht,ui,vi,wi,lu,ixk)
c
c Find heelis drifts from heelis potential (height independent) 
c
      include 'tgcmparam.h'
      include 'input.h'
      include 'mkuivi.h'
c
c Calculate ui,vi,wi from electric potential (input phi), and from
c magnetic field coords bx,by,bz (read from file attached to unit lu)
c
      parameter(imxm1=imx-1,jmxp1=jmx+1,
     +          jmxm1=jmx-1,imxp2=imx+2)
      parameter(rs=6.475165e+8)
      dimension phi(imx,jmx),gcmht(imx,kmx,jmx)
      dimension ui(imx,jmx),vi(imx,jmx),wi(imx,jmx)
      data scale/1.e6/
      data magrd/0/
c
c Define work array with wrap around points:
c
      do 105 j=1,jmx
      do 105 i=1,imxm1
        wrk(i+2,ixk,j+1) = phi(i,j) 
 105  continue
c
c Over the poles:
c
      do 110 i=1,imxm1
        wrk(i+2,ixk,1) = phi(1+mod(i+jmxm1,imxm1),1)
        wrk(i+2,ixk,jmxp2) = phi(1+mod(i+jmxm1,imxm1),jmx)
 110  continue
c
c Periodic points:
c
      do 115 j=1,jmxp2
      do 115 i=1,2
        wrk(i,ixk,j) = wrk(i+imxm1,ixk,j)
        wrk(i+imxp1,ixk,j) = wrk(i+2,ixk,j)
 115  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
c
c Set drifts to 0 for latitudes below 27.5:
c (j=13 is lat -27.5; j=24 is lat +27.5)
c
        if (j.gt.13.and.j.lt.24) then
          ui(:,j) = 0.
          vi(:,j) = 0.
          wi(:,j) = 0.
          goto 200
        endif
        do 205 i=1,imx
          bxyzsum = bx(i,j)**2 + by(i,j)**2 + bz(i,j)**2
          b(i,j) = bxyzsum**0.5
c
c Calculation of electric field:
c
          ex(ixk) = -(wrk(i+3,ixk,j+1) - wrk(i+1,ixk,j+1)) /
     +             (2.*dlon*dtr*(rs+gcmht(i,ixk,j)*1.e+5)*
     +              cos(gcmlat(j)*dtr))
          ey(ixk) = -(wrk(i+2,ixk,j+2) - wrk(i+2,ixk,j)) /
     +             (2.*dlat*dtr*(rs+gcmht(i,ixk,j)*1.e+5))
          ez(ixk) = -(ex(ixk)*bx(i,j) + ey(ixk)*by(i,j)) / bz(i,j) 
c
c Find the ion drifts (m/s):
c
          ui(i,j) = (ey(ixk)*bz(i,j) - ez(ixk)*by(i,j)) * scale / 
     +                b(i,j)**2
          vi(i,j) = (ez(ixk)*bx(i,j) - ex(ixk)*bz(i,j)) * scale / 
     +                b(i,j)**2
          wi(i,j) = (ex(ixk)*by(i,j) - ey(ixk)*bx(i,j)) * scale / 
     +                b(i,j)**2
 205  continue
 200  continue
c     
      return
      end