diff --git a/cgyro/src/cgyro_equilibrium.F90 b/cgyro/src/cgyro_equilibrium.F90
index 453ea1681..55b2493bf 100644
--- a/cgyro/src/cgyro_equilibrium.F90
+++ b/cgyro/src/cgyro_equilibrium.F90
@@ -203,6 +203,25 @@ subroutine cgyro_equilibrium
   ! 2. Compute terms required for O(M) rotation, and no rotation.
   !
 
+  ! EAB LE3 test
+  !print *, geo_f, geo_ffprime/geo_f
+  !open(unit=1,file='out.cgyro.geogk',status='replace')
+  !write (1,'(i2)') n_theta
+  !write (1,'(e16.8)') theta(:)
+  !write (1,'(e16.8)') geo_b(:)
+  !write (1,'(e16.8)') 1.0/(q*rmaj*geo_g_theta(:))
+  !write (1,'(e16.8)') (geo_dbdt(:)/geo_b(:))/(q*rmaj*geo_g_theta(:))
+  !write (1,'(e16.8)') -1.0/geo_b(:)*geo_grad_r(:)/rmaj*geo_gsin(:) 
+  !write (1,'(e16.8)') -1.0/geo_b(:)*geo_gq(:)/rmaj&
+  !     *(geo_gcos1(:)+geo_gcos2(:)+geo_captheta(:)*geo_gsin(:))
+  !write (1,'(e16.8)') 1.0/geo_b(:)*geo_gq(:)/rmaj*geo_gcos2(:)
+  !write (1,'(e16.8)') geo_grad_r(:)**2
+  !write (1,'(e16.8)') geo_gq(:)**2 * (1.0+geo_captheta(:)**2)
+  !write (1,'(e16.8)') 2.0*geo_grad_r(:)*geo_gq(:)*geo_captheta(:)
+  !write (1,'(e16.8)') geo_chi2(:)
+  !write (1,'(e16.8)') -(geo_gcos1(:) + geo_gcos2(:))/rmaj
+  !close(1)
+
   do it=1,n_theta
 
      do is=1,n_species
diff --git a/le3/src/le3_geometry.f90 b/le3/src/le3_geometry.f90
index ce1eaefd1..e6045975d 100644
--- a/le3/src/le3_geometry.f90
+++ b/le3/src/le3_geometry.f90
@@ -46,6 +46,14 @@ subroutine le3_geometry
   allocate(vdrift_x(nt,np))
   allocate(vdrift_dt(nt,np))
   allocate(vdrift_dp(nt,np))
+  allocate(vdrift_gk(nt,np,3))
+  allocate(grad_perpsq_gk(nt,np,3))
+  allocate(grad_cc(nt,np))
+  allocate(grad_tt(nt,np))
+  allocate(grad_pp(nt,np))
+  allocate(grad_pt(nt,np))
+  allocate(grad_ct(nt,np))
+  allocate(grad_cp(nt,np))
   allocate(vexb_dt(nt,np))
   allocate(vexb_dp(nt,np))
   allocate(dgdp(nt,np))
diff --git a/le3/src/le3_geometry_matrix.f90 b/le3/src/le3_geometry_matrix.f90
index f8efac05e..15d612553 100644
--- a/le3/src/le3_geometry_matrix.f90
+++ b/le3/src/le3_geometry_matrix.f90
@@ -2,7 +2,7 @@ subroutine le3_geometry_matrix
 
   use le3_globals
   implicit none
-  integer :: i,j,its,ips,kt,kp
+  integer :: i,j,k,its,ips,kt,kp
   
   real, dimension(:), allocatable :: vec_vdriftx
   real, dimension(:), allocatable :: vec_flux
@@ -17,29 +17,107 @@ subroutine le3_geometry_matrix
   real :: J_boozer, I_boozer
   real, dimension(:), allocatable :: vec_boozer_p, vec_boozer_t, vec_boozer
   real, dimension(:,:), allocatable :: mat_boozer, t_boozer, p_boozer
+  real, dimension(:,:), allocatable :: xtemp
+  real, dimension(:), allocatable :: tmap
   real :: xsum
   integer, dimension(:), allocatable :: i_piv
 
-  ! bhat dot grad = bdotgrad * (iota d/dt + d/dp)  
+  ! anorm* bhat dot grad = bdotgrad * (iota d/dt + d/dp)  
   bdotgrad(:,:) = 1.0/(bmag * g)
 
-  ! (bhat dot grad B)/B
+  ! anorm * (bhat dot grad B)/B
   bdotgradB_overB(:,:) = bdotgrad * (iota * dbdt + dbdp) / bmag
 
-  ! bhat cross grad B dot grad r / B^2
-  vdrift_x(:,:) = 1/(rmin*bmag * g**2) &
+  ! bhat cross grad B dot grad chi / (B^2 * r)
+  vdrift_x(:,:) = 1/(rmin*bmag* g**2) &
        * (-dbdt * (gpp + iota * gpt) + dbdp * (gpt + iota*gtt)) / bmag**2
 
   ! bhat cross grad B dot grad theta / B^2
   vdrift_dt(:,:) = 1/g**2 &
        * ((b1*bmag/chi1 - dtheta*dbdt)*(gpp + iota * gpt) &
-       - dbdp*gc) / bmag**2
+       - dbdp*gc) / bmag**3
 
   ! bhat cross grad B dot grad phi / B^2
   vdrift_dp(:,:) = 1/g**2 &
        * (-(b1*bmag/chi1 - dtheta*dbdt)*(gpt + iota * gtt) &
-       + dbdt*gc) / bmag**2
+       + dbdt*gc) / bmag**3
+ 
+  ! (a/cs)*vdrift_gradB = 1/(cs*anorm*Omega_ca_unit)*(vperp^2/2+vpar^2)*[(1) d/d(r/a) + (2) (i k_theta a)]
+  ! Remap theta: 0..2pi -> -pi..pi
+
+  allocate(xtemp(nt,np))
+  allocate(tmap(nt))
+  do i=1,nt/2
+     k=nt/2+i
+     tmap(i) = t(k)-2*pi
+     tmap(k)   = t(i)
+  enddo
+
+  vdrift_gk(:,:,1) = vdrift_x(:,:)
+  call remap_theta(vdrift_gk(:,:,1),xtemp(:,:),0.0)
+  vdrift_gk(:,:,1) = xtemp(:,:)
+
+  vdrift_gk(:,:,2)  = (-iota*vdrift_dp(:,:) + vdrift_dt(:,:))*rmin
+  call remap_theta(vdrift_gk(:,:,2),xtemp(:,:),0.0)
+  vdrift_gk(:,:,2) = xtemp(:,:)
+  do i=1,nt
+     do j=1,np
+        vdrift_gk(i,j,2) = vdrift_gk(i,j,2) &
+             - vdrift_gk(i,j,1) * rmin**2 *(iota_p/iota)*tmap(i)
+     enddo
+  enddo
 
+  ! (a/cs)*vdrift_gradp =  1/(cs*anorm*Omega_ca_unit)*(vpar^2) [(3) (i k_theta a)]
+  vdrift_gk(:,:,3) = (-0.5*beta_star)/bmag(:,:)**2
+  call remap_theta(vdrift_gk(:,:,3),xtemp(:,:),0.0)
+  vdrift_gk(:,:,3) = xtemp(:,:)
+  
+  ! grad_perpsq: (1) d^2/d(r/a)^2 + (2) (k_theta a)^2 + (3) (i k_theta a) d/d(r/a)
+  grad_cc =  1.0/g**2 * (gtt * gpp - gpt**2)
+  grad_tt =  1.0/g**2 * (gpp * gcc - gcp**2)
+  grad_pp =  1.0/g**2 * (gtt * gcc - gct**2)
+  grad_pt =  1.0/g**2 * (gcp * gct - gpt*gcc)
+  grad_ct =  1.0/g**2 * (gpt * gcp - gpp*gct)
+  grad_cp =  1.0/g**2 * (gpt * gct - gtt*gcp)
+  
+  grad_perpsq_gk(:,:,1) = 1.0/rmin**2 * grad_cc(:,:)
+  call remap_theta(grad_perpsq_gk(:,:,1),xtemp(:,:),0.0)
+  grad_perpsq_gk(:,:,1) = xtemp(:,:)
+  
+  grad_perpsq_gk(:,:,2) = grad_pp + 1.0/iota**2*grad_tt - 2.0/iota*grad_pt
+  call remap_theta(grad_perpsq_gk(:,:,2),xtemp(:,:),0.0)
+  grad_perpsq_gk(:,:,2) = xtemp(:,:)
+  call remap_theta(grad_cc(:,:),xtemp(:,:),0.0)
+  do i=1,nt
+     do j=1,np
+        grad_perpsq_gk(i,j,2) = grad_perpsq_gk(i,j,2) + (iota_p/iota**2*tmap(i))**2 * xtemp(i,j)
+     enddo
+  enddo
+  call remap_theta(grad_cp(:,:),xtemp(:,:),0.0)
+    do i=1,nt
+     do j=1,np
+        grad_perpsq_gk(i,j,2) = grad_perpsq_gk(i,j,2) + 2.0*iota_p/iota**2*tmap(i)*xtemp(i,j)
+     enddo
+  enddo
+  call remap_theta(grad_ct(:,:),xtemp(:,:),0.0)
+    do i=1,nt
+     do j=1,np
+        grad_perpsq_gk(i,j,2) = grad_perpsq_gk(i,j,2) - 2.0*iota_p/iota**3*tmap(i)*xtemp(i,j)
+     enddo
+  enddo
+  grad_perpsq_gk(:,:,2) = grad_perpsq_gk(:,:,2)*(iota*rmin)**2
+  
+  grad_perpsq_gk(:,:,3) = grad_cp - 1.0/iota*grad_ct
+  call remap_theta(grad_perpsq_gk(:,:,3),xtemp(:,:),0.0)
+  grad_perpsq_gk(:,:,3) = xtemp(:,:)
+  call remap_theta(grad_cc(:,:),xtemp(:,:),0.0)
+  do i=1,nt
+     do j=1,np
+        grad_perpsq_gk(i,j,3) =  grad_perpsq_gk(i,j,3) + iota_p/iota**2*tmap(i)* xtemp(i,j)
+     enddo
+  enddo
+  grad_perpsq_gk(:,:,3) =  grad_perpsq_gk(:,:,3)*(-2.0*iota)
+  
   ! flux-surface d volume / dr
   vprime = 0.0
   do i=1,nt
@@ -64,7 +142,7 @@ subroutine le3_geometry_matrix
   vexb_dp(:,:) = -1/(rmin*bmag * g**2) & 
        * (gpt + iota*gtt) * bmag/bsq_avg
 
-  ! construct the geo collocation matices
+  ! construct the geo collocation matices for NEO-3D
 
   allocate(mat_stream_dt(matsize,matsize))
   allocate(mat_stream_dp(matsize,matsize))
@@ -174,7 +252,7 @@ subroutine le3_geometry_matrix
   enddo
   close(1)
 
-  ! Construct the geo vectors
+  ! Construct the geo vectors for NEO-3D
 
   vec_thetabar(:) = vec_thetabar(:) / (nt*np)
   vec_vdriftx(:)  = vec_vdriftx(:) / (nt*np)
@@ -212,7 +290,7 @@ subroutine le3_geometry_matrix
      write (1,'(i2,i2,i2)') itype(i),m_indx(i), n_indx(i)
   enddo
   close(1)
-
+  
   ! Map to Boozer coordinates
 
   ! J: B ~ J(psi) grad phi
@@ -412,6 +490,67 @@ subroutine le3_geometry_matrix
 
   call le3_compute_theory
 
+  ! Write out quantitites for GK-3D
+  ! remap theta: 0..2pi -> -pi..pi
+  
+  open(unit=1,file='out.le3.geogk',status='replace')
+  write (1,'(i4)') nt
+  write (1,'(i4)') np
+  write (1,'(e16.8)') tmap
+  write (1,'(e16.8)') p
+
+  ! theta_bar(theta,phi) and derivatives
+  call remap_theta(tb(:,:),xtemp(:,:),1.0)
+  write (1,'(e16.8)') xtemp(:,:)
+  call remap_theta(dtbdt(:,:),xtemp(:,:),0.0)
+  write (1,'(e16.8)') xtemp(:,:)
+  call remap_theta(dtbdp(:,:),xtemp(:,:),0.0)
+  write (1,'(e16.8)') xtemp(:,:)
+  
+  ! B/Bunit
+  call remap_theta(bmag(:,:),xtemp(:,:),0.0)
+  write (1,'(e16.8)') xtemp(:,:)
+  
+  ! anorm*(bhat dot grad) = bdotgrad * (iota d/dt + d/dp) = (1) d/dt + (2) d/dp
+  call remap_theta(bdotgrad(:,:),xtemp(:,:),0.0)
+  write (1,'(e16.8)') xtemp(:,:)*iota
+  write (1,'(e16.8)') xtemp(:,:)
+  
+  ! anorm*(bhat dot grad B)/B
+  call remap_theta(bdotgradB_overB(:,:),xtemp(:,:),0.0)
+  write (1,'(e16.8)') xtemp(:,:)
+  
+  ! (a/cs)*vdrift_gradB = 1/(cs*anorm*Omega_ca_unit)*(vperp^2/2+vpar^2)*[(1) d/d(r/a) + (2) (i k_theta a)]
+  write (1,'(e16.8)') vdrift_gk(:,:,1)
+  write (1,'(e16.8)') vdrift_gk(:,:,2)
+  
+  ! (a/cs)*vdrift_gradp =  1/(cs*anorm*Omega_ca_unit)*(par^2)*(i k_theta a) [(3) (i ktheta a)]
+  write (1,'(e16.8)') vdrift_gk(:,:,3)
+  
+  ! grad_perpsq: (1) d^2/d(r/a)^2 + (2) (k_theta a)^2 + (3) (i k_theta a) d/d(r/a)
+  write (1,'(e16.8)') grad_perpsq_gk(:,:,1)
+  write (1,'(e16.8)') grad_perpsq_gk(:,:,2)
+  write (1,'(e16.8)') grad_perpsq_gk(:,:,3)
+
+  call remap_theta(b1(:,:),xtemp(:,:),0.0)
+  write (1,'(e16.8)') xtemp(:,:)
+
+  call remap_theta(chi1(:,:),xtemp(:,:),0.0)
+  write (1,'(e16.8)') xtemp(:,:)
+
+  call remap_theta(dchi(:,:),xtemp(:,:),0.0)
+  write (1,'(e16.8)') xtemp(:,:)
+
+  call remap_theta(dtheta(:,:),xtemp(:,:),0.0)
+  write (1,'(e16.8)') xtemp(:,:)
+
+  call remap_theta(dtheta_t(:,:),xtemp(:,:),0.0)
+  write (1,'(e16.8)') xtemp(:,:)
+  
+  close(1)
+  deallocate(xtemp)
+  deallocate(tmap)
+  
   deallocate(itype)
   deallocate(m_indx)
   deallocate(n_indx)
@@ -423,6 +562,9 @@ subroutine le3_geometry_matrix
   deallocate(gpp)
   deallocate(gtt)
   deallocate(gpt)
+  deallocate(gct)
+  deallocate(gcp)
+  deallocate(gcc)
   deallocate(cosu)
   deallocate(btor)
   deallocate(bpol)
@@ -440,6 +582,14 @@ subroutine le3_geometry_matrix
   deallocate(vexb_dp)
   deallocate(vdrift_dt)
   deallocate(vdrift_dp)
+  deallocate(vdrift_gk)
+  deallocate(grad_perpsq_gk)
+  deallocate(grad_cc)
+  deallocate(grad_tt)
+  deallocate(grad_pp)
+  deallocate(grad_pt)
+  deallocate(grad_ct)
+  deallocate(grad_cp)
   deallocate(mat_stream_dt)
   deallocate(mat_stream_dp)
   deallocate(mat_trap)
@@ -458,3 +608,18 @@ subroutine le3_geometry_matrix
   deallocate(vec_ntv)
 
 end subroutine le3_geometry_matrix
+
+subroutine remap_theta(xold,xnew,fac)
+  use le3_globals
+  implicit none
+  real, intent(inout), dimension(nt,np) :: xold,xnew
+  real, intent(in) :: fac
+  integer :: i,j,k
+  do i=1,nt/2
+     k=nt/2+i
+     do j=1,np
+        xnew(i,j)   = xold(k,j) - fac*2*pi
+        xnew(k,j)   = xold(i,j)
+     enddo
+  enddo
+end subroutine remap_theta
diff --git a/le3/src/le3_geometry_rho.f90 b/le3/src/le3_geometry_rho.f90
index 0a2db3f3a..87f6b044c 100644
--- a/le3/src/le3_geometry_rho.f90
+++ b/le3/src/le3_geometry_rho.f90
@@ -13,7 +13,8 @@ subroutine le3_geometry_rho
   real :: ycosuv
   real :: ysinuv
   real :: up,pp,pt
-  real :: d0,eta
+  real :: d0, eta, d1, gpp_1, gtt_1, gpt_1
+  real :: sum1, sum2, sum3, sum4
 
   real, dimension(:,:), allocatable :: s1a,s1b
   real, dimension(:,:), allocatable :: s2a,s2b
@@ -42,9 +43,13 @@ subroutine le3_geometry_rho
   ! Global second order quantities
   allocate(dtheta(nt,np))
   allocate(dchi(nt,np))
-  allocate(dthetap(nt,np))
+  allocate(dtheta_t(nt,np))
+  allocate(dtheta_p(nt,np))
   allocate(b1(nt,np))
   allocate(gc(nt,np))
+  allocate(gct(nt,np))
+  allocate(gcp(nt,np))
+  allocate(gcc(nt,np))
 
   do i=1,nt
      do j=1,np
@@ -62,7 +67,7 @@ subroutine le3_geometry_rho
         pt = (x/rc(i,j))*ycosuv-chi1t(i,j)/chi1(i,j)*ysinuv
 
         eta = 1.0/rc(i,j)+cosu(i,j)/r(i,j)
-
+        
         c_m0 = iota*x**2+x*ycosuv
         c_n0 = r(i,j)**2+y**2+iota*x*ycosuv
 
@@ -86,7 +91,6 @@ subroutine le3_geometry_rho
         a21a(i,j) = (-2*c_m0+iota*x**2)/d0
         a21b(i,j) = x**2/d0
 
-
      enddo
   enddo
 
@@ -180,17 +184,18 @@ subroutine le3_geometry_rho
 
   dtheta  = 0.0
   dchi    = 0.0
-  dthetap = 0.0
+  dtheta_t = 0.0
+  dtheta_p = 0.0
   do k=1,matsize
      call le3_basis(itype(k),m_indx(k),n_indx(k),bk(:,:),'d0')
      call le3_basis(itype(k),m_indx(k),n_indx(k),bk_t(:,:),'dt')
+     call le3_basis(itype(k),m_indx(k),n_indx(k),bk_p(:,:),'dp')
      do j=1,np
         do i=1,nt
-
-           dtheta(i,j) = dtheta(i,j) + sys_b(k)*bk(i,j)
-           dchi(i,j)   = dchi(i,j)   + sys_b(k+matsize)*bk(i,j)
-
-           dthetap(i,j) = dthetap(i,j) + sys_b(k)*bk_t(i,j)
+           dchi(i,j)     = dchi(i,j)     + sys_b(k+matsize)*bk(i,j)
+           dtheta(i,j)   = dtheta(i,j)   + sys_b(k)*bk(i,j)
+           dtheta_t(i,j) = dtheta_t(i,j) + sys_b(k)*bk_t(i,j)
+           dtheta_p(i,j) = dtheta_p(i,j) + sys_b(k)*bk_p(i,j)
 
         enddo
      enddo
@@ -224,19 +229,53 @@ subroutine le3_geometry_rho
         c_dm = x*up+pt+iota*x**2*2.0/rc(i,j)+chi1(i,j)*iota_p*x**2-c_m0*eta 
         c_dn = 2*r(i,j)*cosu(i,j)+iota*x*up+2*pp+iota*pt+chi1(i,j)*iota_p*x*ycosuv-c_n0*eta
 
+        ! D0 = sqrt(g_s)
         d0 = g(i,j)
+
+        ! D1/D0
+        d1 = eta - chi1(i,j)*(dchi(i,j)+dtheta_t(i,j))
+
+        gpp_1 = 2.0*(r(i,j)*cosu(i,j) - chi1p(i,j)/chi1(i,j)*ysinuv &
+             + (up - chi1(i,j)*x*dtheta_p(i,j))*ycosuv)
+
+        gtt_1 = 2.0 * x**2 * (1.0/rc(i,j)- chi1(i,j)*dtheta_t(i,j))
+
+        gpt_1 = x*(up - chi1(i,j)*x*dtheta_p(i,j)) &
+             + (x/rc(i,j) - chi1(i,j)*x*dtheta_t(i,j))*ycosuv &
+             - chi1t(i,j)/chi1(i,j)*ysinuv
+        
         !-------------------------------------------------------------
 
         ! B1/Bs
-        b1(i,j) = 0.5*chi1(i,j)/a12a(i,j)*(iota_p*c_m0/d0+s1a(i,j)) &
-             -0.5*(eta-dchi(i,j)-chi1(i,j)*dthetap(i,j))
-
+        b1(i,j) = -d1 &
+             + 0.5/(c_n0 + iota*c_m0) * (gpp_1 + iota**2 * gtt_1 + 2.0*iota*gpt_1 &
+             + chi1(i,j) * iota_p * (2*iota*gtt(i,j) + gpt(i,j)) )
+        
         ! g_cp + i g_ct
         gc(i,j) = ysinuv/chi1(i,j)-c_m0*dtheta(i,j)
 
+        ! g_ct
+        gct(i,j) = -dtheta(i,j) * x**2
+        
+        ! g_cp
+        gcp(i,j) = ysinuv/chi1(i,j) - dtheta(i,j)*x*ycosuv
+        
+        ! g_cc
+        gcc(i,j) = 1.0/chi1(i,j)**2 + dtheta(i,j)**2 * x**2
+
      enddo
   enddo
 
+  !do i=1,nt
+  !   j=1
+  !   x = sqrt(gtt(i,j))
+     !print *, b1(i,j), -1.0/(r(i,j)**2 + iota**2 * x**2) &
+     !     * (r(i,j)*cosu(i,j) + iota**2*x**2/rc(i,j)), &
+     !     chi1(i,j)*r(i,j)/sqrt(gtt(i,j)), &
+     !     r(i,j)/sqrt(gtt(i,j))&
+     !  *(-1/rc(i,j)+cosu(i,j)/r(i,j)+chi1(i,j)*(dchi(i,j)+dtheta_t(i,j)))/iota
+  !enddo
+  
   !--------------------------------------------------------------------
   ! Check with GEO result
   !
@@ -270,7 +309,7 @@ subroutine le3_geometry_rho
 
   endif
   !--------------------------------------------------------------------
-
+  
   deallocate(sys_m)
   deallocate(sys_b)
   deallocate(i_piv)
diff --git a/le3/src/le3_globals.f90 b/le3/src/le3_globals.f90
index 52e1340fa..033ad57be 100644
--- a/le3/src/le3_globals.f90
+++ b/le3/src/le3_globals.f90
@@ -38,6 +38,7 @@ module le3_globals
   real, dimension(:,:), allocatable  :: sinn,cosn
 
   real, dimension(:,:), allocatable :: gpp,gtt,gpt
+  real, dimension(:,:), allocatable :: gct, gcp, gcc
   real, dimension(:,:), allocatable :: cosu,rc
   real, dimension(:,:), allocatable :: bpol,btor
   real, dimension(:,:), allocatable :: chi1
@@ -51,7 +52,7 @@ module le3_globals
   ! "Second order" quantities
   real, dimension(:,:), allocatable :: dtheta
   real, dimension(:,:), allocatable :: dchi
-  real, dimension(:,:), allocatable :: dthetap
+  real, dimension(:,:), allocatable :: dtheta_t, dtheta_p
   real, dimension(:,:), allocatable :: b1
   real, dimension(:,:), allocatable :: gc
 
@@ -72,6 +73,8 @@ module le3_globals
   real, dimension(:,:), allocatable :: bdotgradB_overB
   real, dimension(:,:), allocatable :: dgdp,vexb_dt,vexb_dp
   real, dimension(:,:), allocatable :: vdrift_x, vdrift_dt, vdrift_dp
+  real, dimension(:,:,:), allocatable :: grad_perpsq_gk, vdrift_gk
+  real, dimension(:,:), allocatable :: grad_cc, grad_tt, grad_pp, grad_pt, grad_ct, grad_cp
   real, dimension(:,:), allocatable :: g
   real, dimension(:,:), allocatable :: mat_stream_dt
   real, dimension(:,:), allocatable :: mat_stream_dp