subroutine geo(dpsidrho)

!  Initialize qinp, shat, b_mag, rkperp2, omegad, bgrad for general geometry
!  from a run of eiktest.  
      
  use itg_data
  implicit none
  
  integer :: ntgrid
  real, dimension(:), allocatable :: theta, bmag, gbdrift, gbdrift0, &
       gds2, gds21, gds22, cvdrift, cvdrift0, gradpar, logb, tgradrho, dum
  real :: dpsidrho, drhodpsin, rmaj
  character*40 char
      
  integer i, ntg, l, m, n
!
!     read in data
!
  open(unit=21,file='eik.out',status='unknown')
  
  read(21,*) char
  read(21,*) ntgrid, i, i, drhodpsin, rmaj, shr
  dpsidrho = 1./drhodpsin
  
  ntg = 2*ntgrid+1
  allocate(theta(ntg), bmag(ntg), gbdrift(ntg), gbdrift0(ntg), &
       gds2(ntg), gds21(ntg), gds22(ntg), cvdrift(ntg), &
       cvdrift0(ntg), gradpar(ntg), logb(ntg), tgradrho(ntg))

  read(21,*) char
  do i=1,ntg
     read(21,*) gbdrift(i), gradpar(i), tgradrho(i), theta(i)
     gbdrift(i) = gbdrift(i)*rmaj/4.
     gradpar(i) = gradpar(i)*rmaj*epsn
     tgradrho(i) = tgradrho(i)/rmaj/epsn
  enddo
  
  read(21,*) char
  do i=1,ntg
     read(21,*) cvdrift(i), gds2(i), bmag(i)
     cvdrift(i) = cvdrift(i)*rmaj/4.
     logb(i) = log(bmag(i))
  enddo
  
  read(21,*) char
  do i=1,ntg
     read(21,*) gds21(i), gds22(i)
  enddo
  
  read(21,*) char
  do i=1,ntg
     read(21,*) cvdrift0(i), gbdrift0(i)
     cvdrift0(i) = cvdrift0(i)*rmaj/4.
     gbdrift0(i) = gbdrift0(i)*rmaj/4.
  enddo
  
  close(21)      
    
  allocate (dum(ld))
  call inter_d_cspl(ntg,theta,logb,ld,r,dum,bgrad)
  deallocate (dum)
  
  call inter_cspl(ntg,theta,bmag,ld,r,b_mag)
  call inter_cspl(ntg,theta,gbdrift, ld,r,gbd )
  call inter_cspl(ntg,theta,gbdrift0,ld,r,gbd0)
  call inter_cspl(ntg,theta,cvdrift, ld,r,cvd )
  call inter_cspl(ntg,theta,cvdrift0,ld,r,cvd0)
  call inter_cspl(ntg,theta,gds2, ld,r,gd2 )
  call inter_cspl(ntg,theta,gds21,ld,r,gd21)
  call inter_cspl(ntg,theta,gds22,ld,r,gd22)
  call inter_cspl(ntg,theta,tgradrho,ld,r,gradrho)
  
  gpar = gradpar(1)
  do l=1,ld
     bgrad(l)=bgrad(l)*gpar*epsn
  enddo
  
! make omegad and rkperp2:
!cfpp$ skip R      
  do i=1,nspecies
     do n=1,nd
        do m=1,md
           do l=1,ld
              rkx(l,m,n)=(mrr(m,n)/y0)*sqrt(gd22(l))/b_mag(l)
!
!     The 1/B above needs to be handled consistently throughout
!     the code.
!     
              omegad(l,m,n,i)=epsn*rky(m,n)*rho(i)*vt(i) &
                   *(gbd(l)+cvd(l)+r0(m,n)*(gbd0(l)+cvd0(l)))
              rkperp2(l,m,n)=rky2(m,n) &
                   *(gd2(l)+2.*r0(m,n)*gd21(l)+r0(m,n)**2*gd22(l))
!
!     special cases for k_theta=0, k_r=const.     
!
              if ((m == meq).and.(md /= 1)) then
                 rkx(l,m,n) = -mrr(m+1,n)/y0*r0(m+1,n)*shr/b_mag(l)
                 omegad(l,m,n,i)=epsn*rho(i)*vt(i) &
                      *rky(m+1,n)*r0(m+1,n)*(gbd0(l)+cvd0(l))
                 rkperp2(l,m,n)=(rky(m+1,n)*r0(m+1,n))**2*gd22(l)
              endif
              
              if (md == 1 .and. mrr(m,n) == 0) then
                 rkx(l,m,n)=1./y0
                 omegad(l,m,n,i)=epsn*rho(i)*vt(i)/y0*(gbd0(l)+cvd0(l))
                 rkperp2(l,m,n)=1./y0**2
              endif
!
!     should be krho
!
              rkx2(l,m,n)=rkx(l,m,n)**2
              rkperp2(l,m,n)=rkperp2(l,m,n)/b_mag(l)**2
           enddo
        enddo
     enddo
  enddo
!cfpp$ noskip R
  
  return
end subroutine geo