subroutine geometry(dpsidrho)

c  Initialize qinp, shat, b_mag, rkperp2, omegad, bgrad for general geometry
c  from a run of eiktest.  
      
      implicit none
c 
      include 'itg.par'
      include 'itg.cmn'

      integer :: ntgrid, nperiod
      real, allocatable, dimension(:) :: theta, bmag, gbdrift, gbdrift0, 
     .     gds2, gds21, gds22, cvdrift, cvdrift0, gradpar, logb, tgradrho
      real dum(lz), dpsidrho, drhodpsin, rmaj
      character*40 char
      
      integer i, ntg, l, m, n
c
c     read in data
c
      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)      


      call inter_d_cspl(ntg,theta,logb,lz,r,dum,bgrad)
      call inter_cspl(ntg,theta,bmag,lz,r,b_mag)
      call inter_cspl(ntg,theta,gbdrift, lz,r,gbd )
      call inter_cspl(ntg,theta,gbdrift0,lz,r,gbd0)
      call inter_cspl(ntg,theta,cvdrift, lz,r,cvd )
      call inter_cspl(ntg,theta,cvdrift0,lz,r,cvd0)
      call inter_cspl(ntg,theta,gds2, lz,r,gd2 )
      call inter_cspl(ntg,theta,gds21,lz,r,gd21)
      call inter_cspl(ntg,theta,gds22,lz,r,gd22)
      call inter_cspl(ntg,theta,tgradrho,lz,r,gradrho)

      gpar = gradpar(1)
      do l=1,ld
         bgrad(l)=bgrad(l)*gpar*epsn
      enddo

c make omegad and rkperp2:
cfpp$ skip R      
      do i=1,nspecies
         do n=1,nd
            do m=1,md
               do l=1,ld
                  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))
c
c     should be krho
c
                  rkperp2(l,m,n)=rkperp2(l,m,n)/b_mag(l)**2
c
c     special cases for k_theta=0, k_r=const.     
c
                  if ((m.eq.meq).and.(md.ne.1)) then
                     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)
                     rkperp2(l,m,n)=rkperp2(l,m,n)/b_mag(l)**2
                  endif

                  if (md.eq.1 .and. mrr(m,n).eq.0) then
                     omegad(l,m,n,i)=epsn*rho(i)*vt(i)/y0
     .                    *(gbd0(l)+cvd0(l))
                     rkperp2(l,m,n)=1./y0**2
                     rkperp2(l,m,n)=rkperp2(l,m,n)/b_mag(l)**2
                  endif
               enddo
            enddo
         enddo
      enddo
cfpp$ noskip R

      return
      end