diff --git a/src/CQL3D_SETUP/ceez.f b/src/CQL3D_SETUP/ceez.f
index 6c513a4..b1ccade 100755
--- a/src/CQL3D_SETUP/ceez.f
+++ b/src/CQL3D_SETUP/ceez.f
@@ -239,7 +239,7 @@ subroutine curv1 (n,x,y,slp1,slpn,islpsw,yp,temp,
       yp(1) = (dx1-slpp1)/diag1
       temp(1) = sdiag1/diag1
       if (n .eq. 2) go to 6
-      do 5 i = 2,nm1
+      do  i = 2,nm1
         delx2 = x(i+1)-x(i)
         if (delx2 .le. 0.) go to 9
         dx2 = (y(i+1)-y(i))/delx2
@@ -249,15 +249,17 @@ subroutine curv1 (n,x,y,slp1,slpn,islpsw,yp,temp,
         temp(i) = sdiag2/diag
         dx1 = dx2
         diag1 = diag2
-    5   sdiag1 = sdiag2
+        sdiag1 = sdiag2
+      end do
     6 diag = diag1-sdiag1*temp(nm1)
       yp(n) = (slppn-dx1-sdiag1*yp(nm1))/diag
 c
 c perform back substitution
 c
-      do 7 i = 2,n
+      do i = 2,n
         ibak = np1-i
-    7   yp(ibak) = yp(ibak)-temp(ibak)*yp(ibak+1)
+        yp(ibak) = yp(ibak)-temp(ibak)*yp(ibak+1)
+      end do
       return
 c
 c too few points
@@ -681,7 +683,8 @@ subroutine curvp1 (n,x,y,p,yp,temp,sigma,ierr)
       diag1 = diag2
       sdiag1 = sdiag2
       if (n .eq. 2) go to 2
-      do 1 i = 2,nm1
+      
+      do i = 2,nm1
         npi = n+i
         delx2 = x(i+1)-x(i)
         if (delx2 .le. 0.) go to 8
@@ -693,29 +696,35 @@ subroutine curvp1 (n,x,y,p,yp,temp,sigma,ierr)
         temp(i) = sdiag2/diag
         dx1 = dx2
         diag1 = diag2
-    1   sdiag1 = sdiag2
+        sdiag1 = sdiag2
+      end do
+      
     2 delx2 = p-(x(n)-x(1))
       dx2 = (y(1)-y(n))/delx2
       call terms (diag2,sdiag2,sigmap,delx2)
       yp(n) = dx2-dx1
       temp(nm1) = temp(2*n-1)-temp(nm1)
+      
       if (n .eq. 2) go to 4
 c
 c perform first step of back substitution
 c
-      do 3 i = 3,n
+      do i = 3,n
         ibak = np1-i
         npibak =n+ibak
         yp(ibak) = yp(ibak)-temp(ibak)*yp(ibak+1)
-    3   temp(ibak) =temp(npibak)-temp(ibak)*temp(ibak+1)
+        temp(ibak) =temp(npibak)-temp(ibak)*temp(ibak+1)
+      end do
+      
     4 yp(n) = (yp(n)-sdiag2*yp(1)-sdiag1*yp(nm1))/
      *        (diag1+diag2+sdiag2*temp(1)+sdiag1*temp(nm1))
 c
 c perform second step of back substitution
 c
       ypn =   yp(n)
-      do 5 i = 1,nm1
-    5    yp(i) = yp(i)+temp(i)*ypn
+      do i = 1,nm1
+         yp(i) = yp(i)+temp(i)*ypn
+      end do
       return
 c
 c too few points
@@ -733,6 +742,8 @@ subroutine curvp1 (n,x,y,p,yp,temp,sigma,ierr)
     8 ierr = 3
       return
       end
+
+!--------------------------------------------------------------------------------      
       function curvp2 (t,n,x,y,p,yp,sigma)
 c
       integer:: n
@@ -3199,6 +3210,9 @@ subroutine snhcsh (sinhm,coshm,x,isw)
       if (isw .eq. 3) sinhm = (expx-1./expx)/(ax+ax)-1.
       return
       end
+
+
+!--------------------------------------------------------------------------------
       subroutine surf1 (m,n,x,y,z,iz,zx1,zxm,zy1,zyn,zxy11,
      *                  zxym1,zxy1n,zxymn,islpsw,zp,temp,
      *                  sigma,ierr)
@@ -3522,54 +3536,64 @@ subroutine surf1 (m,n,x,y,z,iz,zx1,zxm,zy1,zyn,zxy11,
       del1 = x(2)-x(1)
       if (del1 .le. 0.) go to 47
       deli = 1./del1
-      do 38 j = 1,n
+      do j = 1,n
         zp(2,j,2) = deli*(z(2,j)-z(1,j))
-   38   zp(2,j,3) = deli*(zp(2,j,1)-zp(1,j,1))
+        zp(2,j,3) = deli*(zp(2,j,1)-zp(1,j,1))
+      end do
       call terms (diag1,sdiag1,sigmax,del1)
       diagi = 1./diag1
-      do 39 j = 1,n
+      do j = 1,n
         zp(1,j,2) = diagi*(zp(2,j,2)-zp(1,j,2))
-   39   zp(1,j,3) = diagi*(zp(2,j,3)-zp(1,j,3))
+        zp(1,j,3) = diagi*(zp(2,j,3)-zp(1,j,3))
+      end do
+      
       temp(n+1) = diagi*sdiag1
       if (m  .eq. 2) go to 43
-      do 42 i = 2,mm1
+      do i = 2,mm1
         im1 = i-1
         ip1 = i+1
         npi = n+i
         del2 = x(ip1)-x(i)
         if (del2 .le. 0.) go to 47
         deli = 1./del2
-        do 40 j = 1,n
+        do j = 1,n
           zp(ip1,j,2) = deli*(z(ip1,j)-z(i,j))
-   40     zp(ip1,j,3) = deli*(zp(ip1,j,1)-zp(i,j,1))
+          zp(ip1,j,3) = deli*(zp(ip1,j,1)-zp(i,j,1))
+        end do
         call terms (diag2,sdiag2,sigmax,del2)
         diagin = 1./(diag1+diag2-sdiag1*temp(npi-1))
-        do 41 j = 1,n
+        do j = 1,n
           zp(i,j,2) = diagin*(zp(ip1,j,2)-zp(i,j,2)-
      *                        sdiag1*zp(im1,j,2))
-   41     zp(i,j,3) = diagin*(zp(ip1,j,3)-zp(i,j,3)-
-     *                        sdiag1*zp(im1,j,3))
+          zp(i,j,3) = diagin*(zp(ip1,j,3)-zp(i,j,3)-
+     *         sdiag1*zp(im1,j,3))
+        end do
         temp(npi) = diagin*sdiag2
         diag1 = diag2
-   42   sdiag1 = sdiag2
+        sdiag1 = sdiag2
+      end do
+      
    43 diagin = 1./(diag1-sdiag1*temp(npm-1))
-      do 44 j = 1,n
+      do j = 1,n
         npmpj = npm+j
         zp(m,j,2) = diagin*(temp(npmpj)-zp(m,j,2)-
      *                      sdiag1*zp(mm1,j,2))
-   44   zp(m,j,3) = diagin*(temp(j)-zp(m,j,3)-
-     *                      sdiag1*zp(mm1,j,3))
+        zp(m,j,3) = diagin*(temp(j)-zp(m,j,3)-
+     *       sdiag1*zp(mm1,j,3))
+      end do
 c
 c perform back substitution
 c
-      do 45 i = 2,m
+      do i = 2,m
         ibak = mp1-i
         ibakp1 = ibak+1
         npibak = n+ibak
         t = temp(npibak)
-        do 45 j = 1,n
+        do j = 1,n
           zp(ibak,j,2) = zp(ibak,j,2)-t*zp(ibakp1,j,2)
-   45     zp(ibak,j,3) = zp(ibak,j,3)-t*zp(ibakp1,j,3)
+          zp(ibak,j,3) = zp(ibak,j,3)-t*zp(ibakp1,j,3)
+        end do
+      end do 
       return
 c
 c too few points
@@ -3582,6 +3606,8 @@ subroutine surf1 (m,n,x,y,z,iz,zx1,zxm,zy1,zyn,zxy11,
    47 ierr = 2
       return
       end
+
+!--------------------------------------------------------------------------------      
       function surf2 (xx,yy,m,n,x,y,z,iz,zp,sigma)
 c
       integer:: m,n,iz
diff --git a/src/aorsa2dMain.F b/src/aorsa2dMain.F
index 5b42cc3..befc6cc 100755
--- a/src/aorsa2dMain.F
+++ b/src/aorsa2dMain.F
@@ -866,9 +866,9 @@ program aorsa2dMain
         real :: rho_thresh
 
 
-        complex :: aij, beta
+!        complex :: aij, beta
         real,    dimension(:,:), allocatable :: work
-        complex, dimension(:,:), allocatable :: workc
+!        complex, dimension(:,:), allocatable :: workc
 
         complex, dimension(ndfmax) :: row,rowk
 
@@ -911,7 +911,7 @@ program aorsa2dMain
      &           yy(mkdim1 : mkdim2, 1 : nymx) )
 
       allocate (work(nxmx, nymx))
-      allocate (workc(nxmx, nymx))
+!      allocate (workc(nxmx, nymx))
 
       allocate ( psi(nxmx, nymx),  rho(nxmx, nymx), theta(nxmx, nymx),
      &   rhohatx(nxmx, nymx), rhohaty(nxmx, nymx), theta0(nxmx, nymx),
@@ -2255,6 +2255,13 @@ program aorsa2dMain
                uzy(i, j) =   byn(i, j)
                uzz(i, j) =   bzn(i, j)
 
+!     Rotation matrix
+!                      
+!       1-bxn^2   
+!   U = 
+!        bxn         byn         bzn
+!
+               
                dxuxx(i,j)  = dxsqx
                dxxuxx(i,j) = dxxsqx
                dyuxx(i,j)  = dysqx
@@ -3321,7 +3328,8 @@ program aorsa2dMain
                xjz(i,j) = xjantz * shapey  * gausspsi * Iz
 
                if (myid==0) then               
-               write(*,*) "gausspsi", i,j,gausspsi,xjz(i,j)
+                  write(*,*) "gausspsi", i,j,gausspsi,xjz(i,j),psiant,
+     &               psi(i,j), dpsiant,xjantz,shapey,Iz,deltay,yant_max
                end if
             end if
 
@@ -3376,7 +3384,7 @@ program aorsa2dMain
       end if      
         
 ! Set k_phi = n_phi/R_ant where Rant is the location of the current peak
-      write(*,*) 'current location',i_cur_max,j_cur_max,capr(i_cur_max),
+      write(*,*) 'current location',i_cur_max,j_cur_max,
      &   cur_mod_max,capr(i_cur_max),y(j_cur_max)
       xkphia = nphi / capr(i_cur_max)
       npar = xkphia / xk0
@@ -6568,7 +6576,7 @@ program aorsa2dMain
                   p_amat(ipos + 1) = fgk
                   p_amat(ipos + 2) = frk
 
-                  if (icnc1 .eq. 1) then
+                  if (icnc1 .eq. 1) then  !if first global row
                      p_brhs(ipos)     = xb(i,j)  !xb,c,d = current source
                      p_brhs(ipos + 1) = xc(i,j)
                      p_brhs(ipos + 2) = xd(i,j)
@@ -7263,6 +7271,7 @@ program aorsa2dMain
          end do
       end do
 
+      !JCW Why is this being done, just for plots/ doesnt work for mirror?
 *---------------------------------------------------
 *     Interpolate back onto (R, Z) grid for plotting
 *---------------------------------------------------
@@ -11702,7 +11711,7 @@ program aorsa2dMain
       deallocate (xkperp2_slow, xkperp2_fast)
       deallocate (xkprl_a, P_a)
       deallocate (work)
-      deallocate (workc)
+!      deallocate (workc)
 !if not static
       deallocate (workn)
       deallocate (rwork)
diff --git a/src/eqdsk_plot.f90 b/src/eqdsk_plot.f90
index 58fc11f..6c68b41 100644
--- a/src/eqdsk_plot.f90
+++ b/src/eqdsk_plot.f90
@@ -784,14 +784,14 @@ subroutine eqdsk_plot
          nxmx, nymx, nlevmax, title, titx, tity, scalex)
 
 
-      title = 'bmod_mid surfaces'
+      title = 'bmod_{mid} surfaces'
       call ezconc(capr, y, bmod_mid, ff, nnodex, nnodey, numb, &
          nxmx, nymx, nlevmax, title, titx, tity, iflag, scalex)
       if (iflag .eq. 0) call boundary(capr, y, rho, ff, nnodex, &
          nnodey, numb, nxmx, nymx, nlevmax, title, titx, tity,scalex)
 
 
-      title = 'capr_{bpol-mid2} surfaces'
+      title = '(R*Bpol)_{mid2} surfaces'
       call ezconc(capr, y, capr_bpol_mid2, ff, nnodex, nnodey, numb, &
          nxmx, nymx, nlevmax, title, titx, tity, iflag, scalex)
       if (iflag .eq. 0) call boundary(capr, y, rho, ff, nnodex, &
@@ -838,15 +838,15 @@ subroutine eqdsk_plot
          nnodex, nxmx)
 
 
-      title= 'Flux average bmod_mid'
+      title= 'Flux average bmod_{mid}'
       titll= 'bmod_mid (T)'
       titlr='       '
       call ezplot1(title, titll, titlr, rhon_half, bmod_midavg, &
           nnoderho_half, nrhomax)
 
 
-      title= 'Flux average capr_bpol'
-      titll= 'r*bpol (mT)'
+      title= 'Flux average R*B_{pol}'
+      titll= '<r*bpol> (m-T)'
       titlr='       '
       call ezplot1(title, titll, titlr, rhon_half, capr_bpol_midavg, &
            nnoderho_half, nrhomax)
@@ -894,7 +894,7 @@ subroutine eqdsk_plot
 
 
 
-      title = 'grad_parallel B'
+      title = 'grad_{||} B'
       call ezconc(capr, y, gradprlb, ff, nnodex, nnodey, numb, &
          nxmx, nymx, nlevmax, title, titx, tity,iflag, scalex)
       if (iflag .eq. 0) call boundary (capr, y, rho, ff, nnodex, &
diff --git a/src/eqdsk_setup.f b/src/eqdsk_setup.f
index 4e3ab3a..4df5a22 100755
--- a/src/eqdsk_setup.f
+++ b/src/eqdsk_setup.f
@@ -15,10 +15,9 @@ subroutine eqdsk_setup(myid, eqdsk, nmodesx, nmodesy,
 
       external f, error
 
-      common/fcom/fcount, bxn, byn, bzn, bmod, bratio, nxdim, nydim, 
-     &   dx, dy,
-     &   nnodex, nnodey, rt, xwleft, sgn_vprl, modb, bratio_phi, 
-     &   dxdphi, dydphi, caprx
+      common/fcom/ bxn, byn, bzn, bmod, bratio, 
+     &   dx, dy, rt, xwleft, sgn_vprl, modb, bratio_phi, 
+     &   dxdphi, dydphi, caprx, nxdim, nydim,fcount, nnodex, nnodey
 
       common/spline_com/sigma, zbxn, zbyn, zbzn, zbmod, zbratio, 
      &   xprime, yprime
@@ -151,8 +150,8 @@ subroutine eqdsk_setup(myid, eqdsk, nmodesx, nmodesy,
 
 c*** ceez.f arrays:
 
-        real zx1(nymx), zxm(nymx), zy1(nxmx), zyn(nxmx)
-        real zxy11, zxym1, zxy1n, zxymn
+      real zx1(nymx), zxm(nymx), zy1(nxmx), zyn(nxmx)
+      real zxy11, zxym1, zxy1n, zxymn
       real zbxn(nxmx, nymx, 3)
       real zbyn(nxmx, nymx, 3)
       real zbzn(nxmx, nymx, 3)
@@ -691,7 +690,7 @@ subroutine eqdsk_setup(myid, eqdsk, nmodesx, nmodesy,
       drg = rg(2) - rg(1)
       dzg = zg(2) - zg(1)
       
-      drg32 = (rg(nxeqd) - rg(1)) / 32.
+      drg32 = (rg(nxeqd) - rg(1)) / 32.  !jcw magic number
       dzg32 = (zg(nyeqd) - zg(1)) / 32.
       
       imaxis = (rmaxis - rg(1)) / (drg) + 1  
@@ -736,9 +735,9 @@ subroutine eqdsk_setup(myid, eqdsk, nmodesx, nmodesy,
 *     -------------------------
 *     adjust boundaries outward
 *     -------------------------                                  
-      ytop_auto    = ytop_auto    + dzg32           
-      ybottom_auto = ybottom_auto - dzg32                                                       
-      rwright_auto = rwright_auto + drg32 * 2.0     
+      ytop_auto    = ytop_auto    + dzg32
+      ybottom_auto = ybottom_auto - dzg32
+      rwright_auto = rwright_auto + drg32 * 2.0
       rwleft_auto  = rwleft_auto  - drg32 * 2.0
       
       if(ytop .eq. 0.0)    ytop = ytop_auto
@@ -794,7 +793,7 @@ subroutine eqdsk_setup(myid, eqdsk, nmodesx, nmodesy,
 
       xwleft = rwleft - rt
       xwright = rwright - rt
-
+      write(*,*) 'xwleft',xwleft,rwleft,rt,rwright
 
 *--------------------------------------------
 *--   Define x mesh: x(i), xprime(i), capr(i)
@@ -803,7 +802,7 @@ subroutine eqdsk_setup(myid, eqdsk, nmodesx, nmodesy,
 c--   x(i) : -xmax / 2.0   to   xmax / 2.0
       dx = xmax / nnodex
       
-      diffmin = 1.0e+05
+      diffmin = 1.0e+05  !start with large value and reduce at each step
 
       do i = 1, nnodex
          xprime(i) = (i - 1) * dx
@@ -814,13 +813,13 @@ subroutine eqdsk_setup(myid, eqdsk, nmodesx, nmodesy,
 
          xkphi(i) = nphi / capr(i)
          
-         diff = abs(capr(i) - r0)        
-         if(diff .lt. diffmin) then
+         diff = abs(capr(i) - r0)
+         if(diff .lt. diffmin) then  !find smallest diff at i0
             diffmin = diff
             i0 = i
          end if
 
-c         write(6, 1314)i, i0, xprime(i), x(i), capr(i), diff
+         write(6, 1314)i, i0, xprime(i), x(i), capr(i), diff
 
       end do
 
@@ -874,7 +873,7 @@ subroutine eqdsk_setup(myid, eqdsk, nmodesx, nmodesy,
       rplasm = rt + aplasm
       rlim   = rt + alim
 
-
+!JCW rwleft has min of 1e-3
 
 
 *------------------------------
@@ -1917,7 +1916,7 @@ subroutine aorsa_grid(nnodex, nnodey, capr, capz, nxmx, nymx,
       real rho_tor2d(nxmx, nymx)
       real psihigh, psilow
       real bmod(nxmx, nymx)
-      real ps, psr, psz, psrr, psrz, pszz
+      real ps, psr, psz, psrr, psrz, pszz, psz2, psr2
       real f_psi, f_psi_psi, f_3psi
       real f, f_a, f_a_a, f_3a, a
       real q, q_a, q_a_a, q_3a
@@ -1960,10 +1959,11 @@ subroutine aorsa_grid(nnodex, nnodey, capr, capz, nxmx, nymx,
      &            islpsw, zpr, temp,
      &            sigma, ierr)
 
+      zx1 = 0.0
       call surf1 (mr, mz, rg, zg, psizg, nxeqdmax,
      &            zx1, zxm, zy1, zyn,
      &            zxy11, zxym1, zxy1n, zxymn,
-     &            islpsw, zpz, temp,
+     &            254, zpz, temp,
      &            sigma, ierr)
 
       call curv1 (ma, psis, fs, slp1, slpn, islpsw1,
@@ -1977,11 +1977,10 @@ subroutine aorsa_grid(nnodex, nnodey, capr, capz, nxmx, nymx,
 
       do i = 1, nnodex
          do j = 1, nnodey
-            r = capr(i)  !JCW watch out for r=0
+            r = capr(i)
             z = capz(j)
 
 
-
 c           calculate fields:
 
 
@@ -2013,9 +2012,24 @@ subroutine aorsa_grid(nnodex, nnodey, capr, capz, nxmx, nymx,
 *           -------------------------------------------------
 *           Minus sign because psig(i,j) = psilim - psig(i,j)
 *           -------------------------------------------------
-            br =  1./r * psz
-            bz = -1./r * psr
+            if ( (r < 2.e-2) ) then
+               psr2 = surf2(2.e-2, z, mr, mz, rg, zg, psirg, nxeqdmax,
+     &           zpr, sigma)
 
+               psz2 = surf2(capr(2), z, mr, mz, rg, zg, psizg, nxeqdmax,
+     &              zpz, sigma)
+               !r=0 not on mesh , capr(1)=delta r. dpsi/dr(0)=0
+               br =  1./r * psz
+               !bz = -(psr2-psr)/(capr(2)-capr(1))
+               bz = -1./2.e-2 * psr2
+
+#ifdef DEBUG               
+               write(*,*) 'br,bz,axis',br,bz,r,psz,psr,ps,a
+#endif
+            else
+               br =  1./r * psz
+               bz = -1./r * psr
+            end if
 
 *           -----------------------------
 *           To set poloidal field to zero
@@ -2034,7 +2048,7 @@ subroutine aorsa_grid(nnodex, nnodey, capr, capz, nxmx, nymx,
             by0(i, j) = bz
             bz0(i, j) = bphi
 
-            bmod(i, j) = sqrt(br**2 + bz**2 + bphi**2)
+            bmod(i, j) = sqrt(br**2 + bz**2 + bphi**2) !JCW bmod defined here
 
             bxn(i, j) = br / bmod(i, j)
             byn(i, j) = bz / bmod(i, j)
@@ -2049,12 +2063,15 @@ subroutine aorsa_grid(nnodex, nnodey, capr, capz, nxmx, nymx,
 ! JCW define b0 to be consistent with F=R*BPHI on axis regardless of what eqdsk
 ! head said. Rely on it for mirrors though, eventually replace with 
 ! d/dR grad psi
-      if (r0>=1e-4) then !not mirror careful of magic number
-         a = 1.e-04 !look close to magnetic axis
+      if (r0>=1e-2) then !not mirror careful of magic number
+         a = 1.e-03 !look close to magnetic axis
          f = curv2(a, ma, psis, fs, ypf, sigma)
          b0 = f / r0
       else
          b0 = bmod(1,jequat)    ! get |B| at R=0, Z=0 assuming a mirror JCW review
+#ifdef DEBUG
+         write(*,*) 'b0 for mirror', r0,b0
+#endif         
       end if
 
          
@@ -2573,16 +2590,15 @@ subroutine eqdsk_setup2(myid, eqdsk, nmodesx, nmodesy,
      &   bmod_mid, capr_bpol_mid2, capr_bpol_mid, rho_tor2d,
      &   i_psi, dldb_tot12, dldbavg, n_prof_flux, rhomax)
      
-      use size_mod            
+      use size_mod
       
       implicit none
 
       external f, error
 
-      common/fcom/fcount, bxn, byn, bzn, bmod, bratio, nxdim, nydim, 
-     &   dx, dy,
-     &   nnodex, nnodey, rt, xwleft, sgn_vprl, modb, bratio_phi, 
-     &   dxdphi, dydphi, caprx
+      common/fcom/ bxn, byn, bzn, bmod, bratio, 
+     &   dx, dy, rt, xwleft, sgn_vprl, modb, bratio_phi, 
+     &   dxdphi, dydphi, caprx, nxdim, nydim,fcount, nnodex, nnodey
 
       common/spline_com/sigma, zbxn, zbyn, zbzn, zbmod, zbratio, 
      &   xprime, yprime
@@ -2715,8 +2731,8 @@ subroutine eqdsk_setup2(myid, eqdsk, nmodesx, nmodesy,
 
 c*** ceez.f arrays:
 
-        real zx1(nymx), zxm(nymx), zy1(nxmx), zyn(nxmx)
-        real zxy11, zxym1, zxy1n, zxymn
+      real zx1(nymx), zxm(nymx), zy1(nxmx), zyn(nxmx)
+      real zxy11, zxym1, zxy1n, zxymn
       real zbxn(nxmx, nymx, 3)
       real zbyn(nxmx, nymx, 3)
       real zbzn(nxmx, nymx, 3)
diff --git a/src/fieldws.f b/src/fieldws.f
index bc87e2f..76be9d6 100755
--- a/src/fieldws.f
+++ b/src/fieldws.f
@@ -1034,32 +1034,32 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
       titx = '1 / zeta'
       tity = 'dK/dL'
 
-      title = 'Real z0_table1'
+      title = 'Real z0 table1'
       call ezconc(zetai_table, dKdL_table, real(z0_table1),
      &     ff, nmax, mmax, numb, ntable, mtable, nlevmax,
      &     title, titx, tity, iflag,1.0)
 
-      title = 'Imag z0_table1'
+      title = 'Imag z0 table1'
       call ezconc(zetai_table, dKdL_table, aimag(z0_table1),
      &    ff, nmax, mmax,
      &   numb, ntable, mtable, nlevmax, title, titx, tity, iflag,1.0)
 
-      title = 'Real z1_table1'
+      title = 'Real z1 table1'
       call ezconc(zetai_table, dKdL_table, real(z1_table1),
      &   ff, nmax, mmax,
      &   numb, ntable, mtable, nlevmax, title, titx, tity, iflag,1.0)
 
-      title = 'Imag z1_table1'
+      title = 'Imag z1 table1'
       call ezconc(zetai_table, dKdL_table, aimag(z1_table1),
      &    ff, nmax, mmax,
      &   numb, ntable, mtable, nlevmax, title, titx, tity, iflag,1.0)
 
-      title = 'Real z2_table1'
+      title = 'Real z2 table1'
       call ezconc(zetai_table, dKdL_table, real(z2_table1),
      &   ff, nmax, mmax,
      &   numb, ntable, mtable, nlevmax, title, titx, tity, iflag,1.0)
 
-      title = 'Imag z2_table1'
+      title = 'Imag z2 table1'
       call ezconc(zetai_table, dKdL_table, aimag(z2_table1),
      &    ff, nmax, mmax,
      &   numb, ntable, mtable, nlevmax, title, titx, tity, iflag,1.0)
@@ -1112,9 +1112,9 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
      &    nnoderho_half, nrhomax)
 
       title= 'parallel gradient of B'
-      titll= 'gradprlb_avg (T/m)'
+      titll= 'gradprlb avg (T/m)'
       titlr='       '
-
+!JCW note this avg~0 for mirrors, maybe vs Z instead?
       call ezplot0(title, titll, titlr, rhon_half, gradprlb2_avg,
      &    nnoderho_half, nrhomax)
 
@@ -1234,7 +1234,7 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
          xkti3avg(n) = xkti3avg(n) / q
       end do
 
-      title= 'Flux average electron temperature'
+      title= '<electron temperature>'
       titll= 'kTe (eV)'
       titlr='       '
       call ezplot1(title, titll, titlr, rhon_half, xkteavg,
@@ -1640,15 +1640,15 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
      &    nnoderho_half, nrhomax)
       end if
       
-      title= 'normalized viscosity <mu_hat>'
-      titll= 'mu_hat (kg/m**3/s-1)'
+      title= 'normalized viscosity <mu hat>'
+      titll= 'mu hat (kg/m**3/s-1)'
       titlr='       '
       call ezplot0(title, titll, titlr, rhon_half, muhat_avg,
      &    nnoderho_half, nrhomax)
 
 
-      title= 'collisionality (nu_star)'
-      titll= 'nu_star'
+      title= 'collisionality (nu^{*})'
+      titll= 'nu*'
       titlr='       '
       call ezplot0(title, titll, titlr, rhon_half, nu_star_avg,
      &    nnoderho_half, nrhomax)
@@ -1890,14 +1890,14 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
       end do
 
 
-      title = 'Re xkperp_cold'
+      title = 'Re xkperp cold'
       call ezconc(capr, y, freal, ff, nnodex, nnodey, numb,
      &   nxmx, nymx, nlevmax, title, titx, tity, iflag,scalex)
        if (iflag .eq. 0) call boundary(capr, y, rho, ff, nnodex,
      &     nnodey, 1, nxmx, nymx, nlevmax, title, titx, tity,
      &     scalex)
 
-      title = 'Im xkperp_cold'
+      title = 'Im xkperp cold'
       call ezconc(capr, y, fimag, ff, nnodex, nnodey, numb,
      &   nxmx, nymx, nlevmax, title, titx, tity, iflag,scalex)
        if (iflag .eq. 0) call boundary(capr, y, rho, ff, nnodex,
@@ -1913,7 +1913,7 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
          fmidim(i) = aimag(xkperp_cold(i,jmid))
       end do
 
-      title = 'Midplane xkperp_cold'
+      title = 'Midplane xkperp cold'
       titll = 'Re xkperp (1/m)'
       titlr = 'Im xkperp (1/m)'
       titlb = 'R (m)'
@@ -1929,9 +1929,9 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
          fmidim(i) = aimag(xkperp_cold2(i,jmid))
       end do
 
-      title = 'Midplane xkperp_cold2'
-      titll = 'Re xkperp2 (1/m)'
-      titlr = 'Im xkperp2 (1/m)'
+      title = 'Midplane xkperp cold2'
+      titll = 'Re xkperp^{2} (1/m)'
+      titlr = 'Im xkperp^{2} (1/m)'
       titlb = 'R (m)'
 
       call ezplot2(title, titll, titlr, titlb, capr, fmidre, fmidim,
@@ -1947,7 +1947,7 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
       end do
 
 
-      title = 'Re xkperp_cold2'
+      title = 'Re xkperp cold2'
       call ezconc(capr, y, freal, ff, nnodex, nnodey, numb,
      &   nxmx, nymx, nlevmax, title, titx, tity, iflag,scalex)
       if (iflag .eq. 0) call boundary(capr, y, rho, ff, nnodex,
@@ -2245,7 +2245,7 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
             fimag(i,j) = aimag(xjpxe(i,j))
          end do
       end do
-      title = 'Real Je_alpha'
+      title = 'Real Je_{alpha}'
       call ezconc(capr, y, freal, ff, nnodex, nnodey, numb,
      &   nxmx, nymx, nlevmax, title, titx, tity, iflag,scalex)
       if (iflag .eq. 0) call boundary (capr, y, rho, ff, nnodex,
@@ -2259,7 +2259,7 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
             fimag(i,j) = aimag(xjpye(i,j))
          end do
       end do
-      title = 'Real Je_beta'
+      title = 'Real Je_{beta}'
       call ezconc(capr, y, freal, ff, nnodex, nnodey, numb,
      &   nxmx, nymx, nlevmax, title, titx, tity, iflag,scalex)
       if (iflag .eq. 0) call boundary (capr, y, rho, ff, nnodex,
@@ -2273,7 +2273,7 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
             fimag(i,j) = aimag(xjpze(i,j))
          end do
       end do
-      title = 'Real Je_b'
+      title = 'Real Je_{b}'
       call ezconc(capr, y, freal, ff, nnodex, nnodey, numb,
      &   nxmx, nymx, nlevmax, title, titx, tity, iflag,scalex)
       if (iflag .eq. 0) call boundary (capr, y, rho, ff, nnodex,
@@ -2365,7 +2365,7 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
          mod_Eminus_mid(i) = mod_Eminus(i,jmid)
       end do
 
-      title = 'Eminus'
+      title = 'E_{-}'
       titll = 'Re Eminus (V/m)'
       titlr = 'Im Eminus (V/m)'
       titlb = 'R (m)'
@@ -3126,7 +3126,7 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
          end do
       end do
 
-      title = ' Mod(E_minus)'
+      title = ' Mod(Eminus)'
       call ezconc(capr, y, mod_Eminus, ff, nnodex, nnodey, numb,
      &   nxmx, nymx, nlevmax, title, titx, tity, iflag,scalex)
       if (iflag .eq. 0) call boundary (capr, y, rho, ff, nnodex,
@@ -4603,7 +4603,7 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
       titx = 'R (m)'
       tity = 'Z (m)'
 
-      title = 'Real Bx_wave'
+      title = 'Real Bx wave'
       call ezconc(capr, y, freal, ff, nnodex, nnodey, numb,
      &   nxmx, nymx, nlevmax, title, titx, tity, iflag, scalex)
       if (iflag .eq. 0) call boundary (capr, y, rho, ff, nnodex,
@@ -4620,7 +4620,7 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
       end do
 
 
-      title = 'Real Bz_wave'
+      title = 'Real Bz wave'
       call ezconc(capr, y, freal, ff, nnodex, nnodey, numb,
      &   nxmx, nymx, nlevmax, title, titx, tity, iflag,scalex)
       if (iflag .eq. 0) call boundary (capr, y, rho, ff, nnodex,
@@ -4724,9 +4724,9 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
          fmidim(i) = aimag(bxwave(i,jmid))
       end do
 
-      title = 'Bx_wave'
-      titll = 'Re Bx_wave (T)'
-      titlr = 'Im Bx_wave (T)'
+      title = 'Bx wave'
+      titll = 'Re Bx wave (T)'
+      titlr = 'Im Bx wave (T)'
       titlb = 'R (m)'
 
       call ezplot2q(title, titll, titlr, titlb, capr, fmidre, fmidim,
@@ -4741,9 +4741,9 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
          fmidim(i) = aimag(bzwave(i,jmid))
       end do
 
-      title = 'Bz_wave'
-      titll = 'Re Bz_wave (T)'
-      titlr = 'Im Bz_wave (T)'
+      title = 'Bz wave'
+      titll = 'Re Bz wave (T)'
+      titlr = 'Im Bz wave (T)'
       titlb = 'R (m)'
 
       call ezplot2q(title, titll, titlr, titlb, capr, fmidre, fmidim,
@@ -4811,7 +4811,7 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
       end do
 
 
-      title = 'E_minus_flux_plot'
+      title = 'E_{-} flux plot'
       call ezconc(capr, y, freal, ff, nnodex, nnodey, numb,
      &   nxmx, nymx, nlevmax, title, titx, tity, iflag,scalex)
       if (iflag .eq. 0) call boundary (capr, y, rho, ff, nnodex,
@@ -4847,7 +4847,7 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
          fmidim(i) = aimag(eplus_flux_plot(i,jmid))
       end do
 
-      title = 'eplus_flux_plot'
+      title = 'Eplus flux plot'
       titll = 'Re Eplus (V/m)'
       titlr = 'Im Eplus (V/m)'
       titlb = 'R (m)'
@@ -4898,7 +4898,7 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
       end do
 
       title= 'Mod B (theta = 0)'
-      titll= 'bmod_flux (T)'
+      titll= 'bmod flux (T)'
       titlr='       '
       call ezplot0(title, titll, titlr, rhon, bmod_plot,
      &    nnoderho, nrhomax)
@@ -4961,7 +4961,7 @@ subroutine fieldws(dfquotient, rmin_zoom, rmax_zoom)
          end do
       end do
 
-      title = 'E_minus_flux'
+      title = 'E_{-} flux'
       titx   = 'rho'
       tity = 'theta'
       call ezconc(rhon, thetam, freal, ff, nnoderho, mnodetheta, numb,
@@ -6958,7 +6958,11 @@ subroutine ezplot2q(title, titll, titlr, titlb,
       ymax = ymax * 1.1
       ymin = ymin
 #ifdef DEBUG
-      write(*,*) 'title: ', title, rhomin,rhomax,ymin,ymax
+      write(*,*) 'DEBUG ranges ezplot2q: ', rhomin,rhomax,ymin,ymax
+      if (ymax < 1e-20) then
+         write(*,*) 'ymax too small',y1,y2
+         return
+      end if
 #endif
       CALL PGPAGE
       CALL PGSVP (0.15,0.85,0.15,0.85)
@@ -7220,8 +7224,10 @@ subroutine ezplot1(title, titll, titlr, x1, y1, nr, nrmax)
 #ifdef DEBUG
       write(*,*) 'ezplot1:  title: ', title, y1min,y1max
 #endif
-      if(y1max .eq. 0.0 .and. y1min .eq. 0.0)return
-
+      if(y1max .eq. 0.0 .and. y1min .eq. 0.0) then
+         return
+      end if
+      
       ymax = y1max
       ymin = y1min
 
@@ -7510,14 +7516,14 @@ subroutine ezplot1q(title, titll, titlr, titlb, x1, y1, nr, nrmax)
       if (abs(ymin-ymax) .lt. 1e-3) ymax=ymin+1.e-3
 
 #ifdef DEBUG
-      write(*,*) 'title: ', title
+      write(*,*) 'ezplot1q title: ', title
       write(6, *)"ymax = ", ymax, "   ymin = ", ymin
 #endif
 
       iflag = 0
       if(ymax == 0.0 .and. ymin == 0.0) then
+         ymax = 1.0
          iflag = 1
-         return
       end if
       
 c Advance plotter to a new page, define coordinate range of graph and draw axes
diff --git a/src/orbit.f b/src/orbit.f
index a8a4d09..82dda3a 100755
--- a/src/orbit.f
+++ b/src/orbit.f
@@ -18,10 +18,9 @@ subroutine field_line_trace(sgn_vprl_in, i, j,
 
       external f2, error
 
-      common/fcom/fcount, bxn, byn, bzn, bmod, bratio, nxdim, nydim, 
-     .   dx, dy,
-     .   nnodex, nnodey, rt, xwleft, sgn_vprl, modb, bratio_phi, 
-     .   dxdphi, dydphi, caprx
+      common/fcom/ bxn, byn, bzn, bmod, bratio, 
+     &   dx, dy, rt, xwleft, sgn_vprl, modb, bratio_phi, 
+     &   dxdphi, dydphi, caprx, nxdim, nydim,fcount, nnodex, nnodey
 
       common/spline_com/sigma, zbxn, zbyn, zbzn, zbmod, zbratio, 
      .   xprime, yprime
@@ -382,10 +381,9 @@ subroutine f2(x, y, dy)
       real zbratio(nmodesmax, mmodesmax, 3)
       real xprimea(nmodesmax), yprimea(mmodesmax)
 
-      common/fcom/fcount, bxn, byn, bzn, bmod, bratio, nxmx, nymx, 
-     &   dr, dz,
-     &   nnodex, nnodey, rt, xwleft, sgn_vprl, modb, bratio_phi, 
-     &   dxdphi, dydphi, capr
+      common/fcom/ bxn, byn, bzn, bmod, bratio, 
+     &   dr, dz, rt, xwleft, sgn_vprl, modb, bratio_phi, 
+     &   dxdphi, dydphi, capr, nxmx, nymx,fcount, nnodex, nnodey
 
       common/spline_com/sigma, zbxn, zbyn, zbzn, zbmod, zbratio, 
      &   xprimea, yprimea
@@ -454,10 +452,9 @@ subroutine f(x, y, dy)
       real zbratio(nmodesmax, mmodesmax, 3)
       real xprimea(nmodesmax), yprimea(mmodesmax)
 
-      common/fcom/fcount, bxn, byn, bzn, bmod, bratio, nxmx, nymx, 
-     &   dr, dz,
-     &   nnodex, nnodey, rt, xwleft, sgn_vprl, modb, bratio_phi, 
-     &   dxdphi, dydphi, capr
+      common/fcom/ bxn, byn, bzn, bmod, bratio, 
+     &   dr, dz, rt, xwleft, sgn_vprl, modb, bratio_phi, 
+     &   dxdphi, dydphi, capr, nxmx, nymx,fcount, nnodex, nnodey
 
       common/spline_com/sigma, zbxn, zbyn, zbzn, zbmod, zbratio, 
      &   xprimea, yprimea
diff --git a/src/field_line_trace.f b/src/unused_field_line_trace.f
similarity index 100%
rename from src/field_line_trace.f
rename to src/unused_field_line_trace.f