c
c------------------------------------------------------------------
c Begin file /home/sting/foster/proclib/mkuivi.f
c------------------------------------------------------------------
c
      subroutine mkuivi(phi,gcmht,ui,vi,wi,lu)
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,kmx,jmx),gcmht(imx,kmx,jmx)
      dimension ui(imx,kmx,jmx),vi(imx,kmx,jmx),wi(imx,kmx,jmx)
      data scale/1.e6/
      data magrd/0/
c
c     write(6,"('mkuivi: ncalls=',i3,' lu=',i3)") ncalls,lu
c     do j=1,jmx,6
c     do k=1,kmx,5
c       write(6,"('k j=',2i3,' gcmht=',/(6e12.4))") k,j,
c    +    (gcmht(i,k,j),i=1,imx,2)
c     enddo
c     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):
c
        do 215 k=1,kmx
          ui(i,k,j) = (ey(k)*bz(i,j) - ez(k)*by(i,j)) * scale / 
     +                b(i,j)**2
          vi(i,k,j) = (ez(k)*bx(i,j) - ex(k)*bz(i,j)) * scale / 
     +                b(i,j)**2
          wi(i,k,j) = (ex(k)*by(i,j) - ey(k)*bx(i,j)) * scale / 
     +                b(i,j)**2
 215    continue
 200  continue
c     
      return
      end