diff --git a/bbb/bbb.v b/bbb/bbb.v
index 9542e53d..f3047bca 100644
--- a/bbb/bbb.v
+++ b/bbb/bbb.v
@@ -1287,6 +1287,8 @@ qe              real            /1.6022e-19/            #Elementary charge
 mu0             real            /1.25663706e-6/         #Vac. magnetic perm.
 eps0            real            /8.8542e-12/            #Vac. dielectric perm.
 rt8opi          real            /1.595769121606e0/      #sqrt(8/pi)
+log_10          real            /2.302585092994046/     #log(10)
+log10ev          real           /-18.795283272599516/     #log10(ev)
 
 ***** Comtra restart:
 #Variables that contain the transport parameters.
@@ -1593,7 +1595,7 @@ mg(1:ngsp)         _real [kg] /1.67e-27/ #gas species mass, calc. fr minu
 facmg(1:nispmx)        real /nispmx*1./  #scale factor for mg to recov old case
 znucl(1:nisp)              _integer [ ]   #tot. nucl. charge, calc. from znuclin
 ni(0:nx+1,0:ny+1,1:nisp)   _real  [m^-3]  +threadprivate #ion density in primary cell (ix,iy)
-lni(0:nx+1,0:ny+1,1:nisp)  _real  [m^-3]  #log(ion dens) in prim. cell (ix,iy)
+logni(0:nx+1,0:ny+1,1:nisp)   _real  [m^-3]   +threadprivate #log ion density in primary cell (ix,iy)
 nm(0:nx+1,0:ny+1,1:nisp)   _real [kg*m^-3] +threadprivate #mass density [nm(,,1) is sum, exclud.
                                           #gas, if nusp=1, isimpon=5] in cell
 nz2(0:nx+1,0:ny+1)         _real  [m^-3]  +threadprivate #sum of ni*zi**2 over all ion species
@@ -1641,7 +1643,6 @@ logte(0:nx+1,0:ny+1)       _real  [J]	  +threadprivate #log electron temperature
 logti(0:nx+1,0:ny+1)       _real  [J]	  +threadprivate #log ion temperature in primary cell
 ng(0:nx+1,0:ny+1,1:ngsp)   _real  [m^-3]  +threadprivate #gas density in primary cell (ix,iy)
 logng(0:nx+1,0:ny+1,1:ngsp)   _real  [m^-3]   +threadprivate #log gas density in primary cell (ix,iy)
-lng(0:nx+1,0:ny+1,1:ngsp)  _real  [m^-3]  +threadprivate #log(gas dens) in prim. cell (ix,iy)
 uug(0:nx+1,0:ny+1,1:ngsp)  _real  [m/s]   +threadprivate #ratio gas-flux/density at x-face;
                                             #if orthog mesh, poloidal gas velocity
 uuxg(0:nx+1,0:ny+1,1:ngsp) _real  [m/s]   +threadprivate #poloidal-only component of uug;
@@ -3062,8 +3063,9 @@ ix2g(0:nx+1,0:ny+1)          _integer  #ix index for sec. interm. (ix,iy) pt.
 iy2g(0:nx+1,0:ny+1)          _integer  #iy index for sec. interm. (ixo,iy) pt.
 ixv2g(0:nx+1,0:ny+1)         _integer  #ixv index for sec. interm.(ixvo,iy) pt.
 iyv2g(0:nx+1,0:ny+1)         _integer  #iyv index for sec.interm.(ixvo,iy) pt.
-nis(0:nxold+1,0:nyold+1,1:nisp) _real [m^-3] +state
-                                             #ion dens at last success. calc
+nis(0:nxold+1,0:nyold+1,1:nisp) _real [m^-3] +state          
+                                   #ion dens at last success. calc
+lognis(0:nxold+1,0:nyold+1,1:nisp) _real [m^-3] +restart #log ion dens at last success. calc
 tes(0:nxold+1,0:nyold+1)        _real [J]    #elec. temp at last success. calc
 tis(0:nxold+1,0:nyold+1)        _real [J]    #ion temp at last success. calc
 tgs(0:nxold+1,0:nyold+1,1:ngsp) _real [J]    #gas temp at last success. calc
diff --git a/bbb/boundary.m b/bbb/boundary.m
index ab6140b0..ce53bcab 100755
--- a/bbb/boundary.m
+++ b/bbb/boundary.m
@@ -4827,7 +4827,7 @@ c_mpi        call mpi_allreduce(sumarr,gsumarr,2,MPI_DOUBLE_PRECISION,MPI_SUM,
 c   Compute fluxes along left (inner divertor - need to change sign)
       do jx = 1, nxpt  #number of x-points; =1 for single null
         do iy = 1, ny
-          sxcos = sx(ixlb(jx),iy)*cos(angfx(ixlb(jx),iy))
+          sxcos = sx(ixlb(jx),iy)*cosangfx(ixlb(jx),iy)
           ue_ion_flux_sum = 0.
           do ifld = 1, 1  #placeholder; needs extension for impurities
             ue_part_fluxh2p1lb(iy,jx) = -fnix(ixlb(jx),iy,1)/sxcos
@@ -4858,7 +4858,7 @@ c   Compute fluxes along right (outer divertor)
         nxt = ixrb(jx)-1   #effec nx
         nxt1 = ixrb(jx)    #effec nx+1
         do iy = 1, ny
-          sxcos = sx(nxt,iy)*cos(angfx(nxt,iy))
+          sxcos = sx(nxt,iy)*cosangfx(nxt,iy)
           ue_ion_flux_sum = 0.
           do ifld = 1, 1  #placeholder; needs extension for impurities
             ue_part_fluxh2p1rb(iy,jx) = fnix(nxt,iy,1)/sxcos
diff --git a/bbb/convert.m b/bbb/convert.m
index 1dd8f7ed..b684b97e 100755
--- a/bbb/convert.m
+++ b/bbb/convert.m
@@ -107,7 +107,7 @@
             if(ineudif .ne. 3) then
               yl(iv) = ng(ix,iy,igsp)/n0g(igsp)
             elseif(ineudif .eq. 3) then
-              yl(iv) = lng(ix,iy,igsp)
+              yl(iv) = logng(ix,iy,igsp)
             endif
             rtol(iv) = rtolv(igrid)*bfac
             atol(iv) = cngatol*rtol(iv)*bfac*abs(yl(iv))
@@ -261,7 +261,8 @@ c_mpi      integer status(MPI_STATUS_SIZE)
                ne(ix,iy) = ne(ix,iy) + zi(ifld)*ni(ix,iy,ifld)
                if (isupgon(1).eq.1 .and. zi(ifld).eq.0) then
                   ng(ix,iy,1) = ni(ix,iy,ifld)
-                  if (ineudif .eq. 3) lng(ix,iy,1)=log(ng(ix,iy,1))
+                  logng(ix,iy,1)=log(abs(ng(ix,iy,1)))
+                  if (ineudif .eq. 3) logng(ix,iy,1)=log(ng(ix,iy,1))
                else
                   nit(ix,iy) = nit(ix,iy) + ni(ix,iy,ifld)
                   if (isimpon.ge.5 .and. nusp_imp.eq.0) 
@@ -280,6 +281,7 @@ c_mpi      integer status(MPI_STATUS_SIZE)
 	    if(isteonxy(ix,iy) .eq. 1) then
                te(ix,iy)=yl(idxte(ix,iy))*ennorm/(1.5*ntemp)
                te(ix,iy) = max(te(ix,iy), temin*ev)  #NEW Feb4,2018
+               logte(ix,iy)=log(abs(te(ix,iy)))
             endif
             do 65 igsp =1, ngsp
 	       if(isngonxy(ix,iy,igsp) .eq. 1) then
@@ -292,9 +294,10 @@ c_mpi      integer status(MPI_STATUS_SIZE)
 		     igspneg = igsp
                    endif
                  elseif(ineudif .eq. 3) then
-                   lng(ix,iy,igsp) = yl(idxg(ix,iy,igsp))
-                   ng(ix,iy,igsp) = exp(lng(ix,iy,igsp))
+                   logng(ix,iy,igsp) = yl(idxg(ix,iy,igsp))
+                   ng(ix,iy,igsp) = exp(logng(ix,iy,igsp))
                  endif
+                 logng(ix,iy,igsp)=log(abs(ng(ix,iy,igsp)))
                endif
 	       if(istgonxy(ix,iy,igsp) .eq. 1) then
                  ntemp = ng(ix,iy,igsp)
@@ -303,6 +306,7 @@ c_mpi      integer status(MPI_STATUS_SIZE)
      .                                                  (1.5*ntemp)
                  tg(ix,iy,igsp) = max(tg(ix,iy,igsp), tgmin*ev) 
                endif
+               logtg(ix,iy,igsp)=log(abs(tg(ix,iy,igsp)))
  65         continue
             ntemp = nit(ix,iy) + cngtgx(1)*ng(ix,iy,1)
             if(isflxvar .eq. 0) ntemp = nnorm
@@ -515,8 +519,11 @@ subroutine convsr_aux (ixl, iyl)
          do 12 iy = js, je
             do 11 ix = is, ie
                pri(ix,iy,ifld) = ni(ix,iy,ifld) * ti(ix,iy)
-               if (ifld.eq.iigsp .and. istgon(1).eq.1) pri(ix,iy,ifld) = 
-     .                                     ni(ix,iy,ifld) * tg(ix,iy,1)
+               logpri(ix,iy,ifld) = logni(ix,iy,ifld) + logti(ix,iy)
+               if (ifld.eq.iigsp .and. istgon(1).eq.1) then
+                    pri(ix,iy,ifld) = ni(ix,iy,ifld) * tg(ix,iy,1)
+                    logpri(ix,iy,ifld) = logni(ix,iy,ifld) + logtg(ix,iy,1)
+               endif
                if (zi(ifld).ne.0.) then
                   pr(ix,iy) = pr(ix,iy) + pri(ix,iy,ifld)
                   zeff(ix,iy)=zeff(ix,iy)+zi(ifld)**2*ni(ix,iy,ifld)
@@ -552,6 +559,7 @@ ccc               ne(ix,iy) = rnec (ni(ix,iy,1), nzsp, ni(ix,iy,nhsp+1),
      .                       (1-istgcon(igsp))*rtg2ti(igsp)*ti(ix,iy) + 
      .                          istgcon(igsp)*tgas(igsp)*ev
 	       pg(ix,iy,igsp) = ng(ix,iy,igsp)*tg(ix,iy,igsp)
+	       logpg(ix,iy,igsp) = logng(ix,iy,igsp) + logtg(ix,iy,igsp)
            enddo
    15    continue
  16   continue
diff --git a/bbb/geometry.m b/bbb/geometry.m
index e136b025..cb06dee3 100755
--- a/bbb/geometry.m
+++ b/bbb/geometry.m
@@ -868,13 +868,13 @@ call s2fill (nx+2, ny+2, 0., fxpyv(0:nx+1,0:ny+1,iu), 1, nx+2)
             call s2fill (nx+2, ny+2, 0., fymxv(0:nx+1,0:ny+1,iu), 1, nx+2)
             call s2fill (nx+2, ny+2, 0., fypxv(0:nx+1,0:ny+1,iu), 1, nx+2)
          enddo
-      if (isnonog .ge. 1) then
-         call nonorthg   # Sets dist btwn interp pts as dxnog, dynog
+         if (isnonog .ge. 1) then
+            call nonorthg   # Sets dist btwn interp pts as dxnog, dynog
                          # Also sets nonorth stencils:fx0,fxm,fy0,fym etc
          do ix = 0, nx+1    # Bdry value set: matters little except vygtan
 ccc            gxfn(ix,0) = gxfn(ix,1)
             dxnog(ix,0) = dxnog(ix,1)
-	    dxnog(ix,ny+1) = dxnog(ix,ny)
+            dxnog(ix,ny+1) = dxnog(ix,ny)
          enddo
          do iy = 0, ny+1  # likewise, dxnog(nx+1,) not relevant, but avoid 0
            dxnog(nx+1,iy) = 0.1*dxnog(nx,iy)
@@ -906,10 +906,10 @@ ccc            gxfn(ix,0) = gxfn(ix,1)
              dx(ix,iy) = max(dx(ix,iy), dxmin)
 cc      for nonorthogonal mesh, dy gets reduced by cosine(average_angle)
             if ((islimon .ne. 0) .and. (ix .eq. ix_lim+1)) then
-               dy(ix,iy) = dy(ix,iy)*cos( angfx(ix,iy) )
+               dy(ix,iy) = dy(ix,iy)*cos(angfx(ix,iy))
             elseif (nxpt==2 .and. ix==ixlb(nxpt)) then  # full double null
 c              Effectively, use angfx(ix-1,iy)=angfx(ix,iy) at left boundaries
-               dy(ix,iy) = dy(ix,iy)*cos( angfx(ix,iy) )
+               dy(ix,iy) = dy(ix,iy)*cos(angfx(ix,iy))
             else
                dy(ix,iy) = dy(ix,iy)*cos( 0.5*
      .                          (angfx(ix,iy) + angfx(ixm1(ix,iy),iy)) )
@@ -1957,9 +1957,11 @@ c     of cell (ix,iy) and vertex 1 of cell (ix+1,iy+1).
          iym1 = max(0, iy-1)
          do ix = 0, nx+1
             angfx(ix,iy) = 0.5*(vtag(ix,iy)+vtag(ix,iym1))
+            cosangfx(ix,iy) = COS(angfx(ix,iy))
             if (islimon .eq. 1) then
                if (iy.eq.iy_lims) then
                   angfx(ix,iy) = vtag(ix,iy)
+                  cosangfx(ix,iy) = COS(angfx(ix,iy))
                endif
                if (iy .eq. iy_lims-1) then
                   angfx(ix,iy) = vtag(ix,iym1)
@@ -1976,6 +1978,7 @@ c     of cell (ix,iy) and vertex 1 of cell (ix+1,iy+1).
       if(redopltvtag == 1) then
        do iy = 0, ny+1
          angfx(0,iy) = 2*angfx(1,iy) - angfx(2,iy)
+         cosangfx(0,iy) = COS(angfx(0,iy))
          vtag(0,iy) = 2*vtag(1,iy) - vtag(2,iy)
          if (nxomit .gt. 0) then
             angfx(0,iy) = angfx(1,iy)
@@ -1983,23 +1986,31 @@ c     of cell (ix,iy) and vertex 1 of cell (ix+1,iy+1).
          endif
          if (islimon .ne. 0) then
             angfx(ix_lim-1,iy) = 2*angfx(ix_lim-2,iy)-angfx(ix_lim-3,iy)
+            cosangfx(ix_lim-1,iy) = COS(angfx(ix_lim-1,iy))
             vtag(ix_lim-1,iy) = 2*vtag(ix_lim-2,iy)-vtag(ix_lim-3,iy)
             angfx(ix_lim,iy) = angfx(ix_lim-1,iy)
+            cosangfx(ix_lim,iy) = COS(angfx(ix_lim,iy))
             vtag(ix_lim,iy) = vtag(ix_lim-1,iy)
             angfx(ix_lim+1,iy) = 2*angfx(ix_lim+2,iy)-angfx(ix_lim+3,iy)
+            cosangfx(ix_lim+1,iy) = COS(angfx(ix_lim+1,iy))
             vtag(ix_lim+1,iy) = 2*vtag(ix_lim+2,iy)-vtag(ix_lim+3,iy)
          endif
          if (nxpt==2) then                # full double null
            angfx(ixrb(1),iy) = 2*angfx(ixrb(1)-1,iy) - angfx(ixrb(1)-2,iy)
+           cosangfx(ixrb(1),iy) = COS(angfx(ixrb(1),iy))
            vtag(ixrb(1),iy) = 2*vtag(ixrb(1)-1,iy) - vtag(ixrb(1)-2,iy)
            angfx(ixrb(1)+1,iy) = angfx(ixrb(1),iy)
+           cosangfx(ixrb(1)+1,iy) = COS(angfx(ixrb(1)+1,iy))
            vtag(ixrb(1)+1,iy) = vtag(ixrb(1),iy)
            angfx(ixlb(2),iy) = 2*angfx(ixlb(2)+1,iy) - angfx(ixlb(2)+2,iy)
+           cosangfx(ixlb(2),iy) = COS(angfx(ixlb(2),iy))
            vtag(ixlb(2),iy) = 2*vtag(ixlb(2)+1,iy) - vtag(ixlb(2)+2,iy)
          endif
          angfx(nx,iy) = 2*angfx(nx-1,iy) - angfx(nx-2,iy)
+         cosangfx(nx,iy) = COS(angfx(nx,iy))
          vtag(nx,iy) = 2*vtag(nx-1,iy) - vtag(nx-2,iy)
          angfx(nx+1,iy) = angfx(nx,iy)    # not really used
+         cosangfx(nx+1,iy) = COS(angfx(nx+1,iy))
          vtag(nx+1,iy) = vtag(nx,iy)      # not really used
        enddo
       endif
@@ -2007,7 +2018,9 @@ c     of cell (ix,iy) and vertex 1 of cell (ix+1,iy+1).
 c...  Set angfx and vtag in the y-guard cells to adjacent values
       do ix = 0, nx+1
          angfx(ix,0) = angfx(ix,1)
+         cosangfx(ix,0) = COS(angfx(ix,0))
          angfx(ix,ny+1) = angfx(ix,ny)
+         cosangfx(ix,ny+1) = COS(angfx(ix,ny+1))
          vtag(ix,0) = vtag(ix,1)
          vtag(ix,ny+1) = vtag(ix,ny)
       enddo
@@ -2017,13 +2030,19 @@ c     of cell (ix,iy) and vertex 1 of cell (ix+1,iy+1).
       do jx = 1, nxpt
          if (ixpt1(jx) .gt. 0) then
             angfx(ixpt1(jx),iysptrx1(jx)-1) = 0.
+            cosangfx(ixpt1(jx),iysptrx1(jx)-1) = COS(0.)
             angfx(ixpt1(jx),iysptrx1(jx)  ) = 0.
+            cosangfx(ixpt1(jx),iysptrx1(jx)) = COS(0.)
             angfx(ixpt1(jx),iysptrx1(jx)+1) = 0.
+            cosangfx(ixpt1(jx),iysptrx1(jx)+1) = COS(0.)
          endif
          if (ixpt2(jx) .gt. 0) then
             angfx(ixpt2(jx),iysptrx2(jx)-1) = 0.
+            cosangfx(ixpt2(jx),iysptrx2(jx)-1) = COS(0.)
             angfx(ixpt2(jx),iysptrx2(jx)  ) = 0.
+            cosangfx(ixpt2(jx),iysptrx2(jx)) = COS(0.)
             angfx(ixpt2(jx),iysptrx2(jx)+1) = 0.
+            cosangfx(ixpt2(jx),iysptrx2(jx)+1) = COS(0.)
          endif
       enddo
 
@@ -2031,11 +2050,13 @@ c     of cell (ix,iy) and vertex 1 of cell (ix+1,iy+1).
       if (isfixlb(1).eq.2 .and. ixpt2(1).gt.0) then
          do iy = 0, iysptrx1(1)
             angfx(ixpt2(1),iy) = 0.
+            cosangfx(ixpt2(1),iy) = COS(0.)
           enddo
        endif
        if (isfixrb(1).eq.2 .and. ixpt1(1).gt.0) then
          do iy = 0, iysptrx1(1)
             angfx(ixpt1(1),iy) = 0.
+            cosangfx(ixpt1(1),iy) = COS(0.)
           enddo
        endif
         
@@ -2044,8 +2065,11 @@ c     of cell (ix,iy) and vertex 1 of cell (ix+1,iy+1).
       if ((isudsym==1.or.(geometry.eq.'dnXtarget')) .and. nxc.gt.0) then
          do iy = 0, ny+1
             angfx(nxc-1,iy) = 0.
+            cosangfx(nxc-1,iy) = COS(0.)
             angfx(nxc  ,iy) = 0.
+            cosangfx(nxc,iy) = COS(0.)
             angfx(nxc+1,iy) = 0.
+            cosangfx(nxc+1,iy) = COS(0.)
             vtag(nxc-1,iy) = 0.
             vtag(nxc  ,iy) = 0.
             vtag(nxc+1,iy) = 0.
diff --git a/bbb/oderhs.m b/bbb/oderhs.m
index 626831db..9cb3b834 100755
--- a/bbb/oderhs.m
+++ b/bbb/oderhs.m
@@ -1449,7 +1449,7 @@ cc              if (difnit(ifld) .gt. 1.e-20 .and. zi(ifld) .eq. 1.
      .                     fypx(ix,iy,0)*log(ni(ix6,iy2,ifld)) ) ) /
      .                                                      dxnog(ix,iy)
             vytan(ix,iy,ifld)=(fcdif*difni(ifld) + dif_use(ix,iy,ifld)) *
-     .                                      (grdnv/cos(angfx(ix,iy)) - 
+     .                                      (grdnv/cosangfx(ix,iy) - 
      .                       (log(ni(ix2,iy,ifld)) - log(ni(ix,iy,ifld)))
      .                                                 * gxf(ix,iy) )
             if (islimon.eq.1.and. ix.eq.ix_lim.and. iy.ge.iy_lims) then
@@ -4345,7 +4345,7 @@ c               if (ifld .ne. iigsp) then
      .                         (log(te(ix2,iy)) + log(te(ix,iy))) )* 
      .                               (fcdif*kye+kye_use(ix,iy))*0.5*
      .                                       (ne(ix2,iy)+ne(ix,iy))*
-     .                                     (grdnv/cos(angfx(ix,iy)) - 
+     .                                     (grdnv/cosangfx(ix,iy) - 
      .                         (log(te(ix2,iy)) - log(te(ix,iy)))* 
      .                                         gxf(ix,iy))*sx(ix,iy)
 
@@ -4375,7 +4375,7 @@ c               if (ifld .ne. iigsp) then
      .                                     (nit(ix2,iy)+nit(ix,iy))
      .          + cftiexclg*cfneut*cfneutsor_ei*0.25*(hcyn(ix ,iy)+hcyn(ix ,iy1)
      .                              +hcyn(ix2,iy)+hcyn(ix4,iy1)) ) *
-     .                                 (  grdnv/cos(angfx(ix,iy))
+     .                                 (  grdnv/cosangfx(ix,iy)
      .                         - (log(ti(ix2,iy)) - log(ti(ix,iy)))*
      .                                        gxf(ix,iy) )*sx(ix,iy)
 c...  Flux limit with flalftxt even though hcys have parallel FL built in
@@ -5661,22 +5661,22 @@ c  If we are solving for the neutral parallel momentum (isupgon=1)
             endif
             vygtan(ix,iy,igsp) = 0.  # vygtan is grad(T) rad vel on x-face
             if (isnonog .eq. 1 .and. iy .le. ny) then
-                  grdnv =(    ( fym (ix,iy,1)*log(tg(ix2,iy1,igsp)) +
-     .                          fy0 (ix,iy,1)*log(tg(ix2,iy ,igsp)) +
-     .                          fyp (ix,iy,1)*log(tg(ix2,iy2,igsp)) +
-     .                          fymx(ix,iy,1)*log(tg(ix ,iy1,igsp)) +
-     .                          fypx(ix,iy,1)*log(tg(ix, iy2,igsp)) )
-     .                       -( fym (ix,iy,0)*log(tg(ix ,iy1,igsp)) +
-     .                          fy0 (ix,iy,0)*log(tg(ix ,iy ,igsp)) +
-     .                          fyp (ix,iy,0)*log(tg(ix ,iy2,igsp)) + 
-     .                          fymx(ix,iy,0)*log(tg(ix4,iy1,igsp)) +
-     .                          fypx(ix,iy,0)*log(tg(ix6,iy2,igsp)) ) )/ 
+                  grdnv =(    ( fym (ix,iy,1)*logtg(ix2,iy1,igsp) +
+     .                          fy0 (ix,iy,1)*logtg(ix2,iy ,igsp) +
+     .                          fyp (ix,iy,1)*logtg(ix2,iy2,igsp) +
+     .                          fymx(ix,iy,1)*logtg(ix ,iy1,igsp) +
+     .                          fypx(ix,iy,1)*logtg(ix, iy2,igsp) )
+     .                       -( fym (ix,iy,0)*logtg(ix ,iy1,igsp) +
+     .                          fy0 (ix,iy,0)*logtg(ix ,iy ,igsp) +
+     .                          fyp (ix,iy,0)*logtg(ix ,iy2,igsp) + 
+     .                          fymx(ix,iy,0)*logtg(ix4,iy1,igsp) +
+     .                          fypx(ix,iy,0)*logtg(ix6,iy2,igsp) ) )/ 
      .                                                      dxnog(ix,iy)
                vygtan(ix,iy,igsp) = exp( 0.5*
-     .                     (log(tg(ix2,iy,igsp))+log(tg(ix,iy,igsp))) )*
+     .                     (logtg(ix2,iy,igsp)+logtg(ix,iy,igsp)) )*
      .                      ( cngfx(igsp) / (mg(igsp)*0.5*(nu1+nu2)) ) *
-     .                                     ( grdnv/cos(angfx(ix,iy)) - 
-     .                       (log(tg(ix2,iy,igsp)) - log(tg(ix,iy,igsp)))
+     .                                     ( grdnv/cosangfx(ix,iy) - 
+     .                       (logtg(ix2,iy,igsp) - logtg(ix,iy,igsp))
      .                                                 * gxf(ix,iy) ) 
              if (islimon.eq.1.and. ix.eq.ix_lim.and. iy.ge.iy_lims) then
                vygtan(ix,iy,igsp) = 0.
@@ -5756,16 +5756,16 @@ c  If we are solving for the neutral parallel momentum (isupgon=1)
      .                fxmy(ix,iy,1)*tg(ixm1(ix,iy)  ,iy  ,igsp) +
      .                fxpy(ix,iy,1)*tg(ixp1(ix,iy)  ,iy  ,igsp)
               elseif (isintlog .eq. 1) then
-                ty0=exp(fxm (ix,iy,0)*log(tg(ixm1(ix,iy)  ,iy  ,igsp)) + 
-     .                  fx0 (ix,iy,0)*log(tg(ix           ,iy  ,igsp)) +
-     .                  fxp (ix,iy,0)*log(tg(ixp1(ix,iy)  ,iy  ,igsp)) +
-     .                  fxmy(ix,iy,0)*log(tg(ixm1(ix,iy+1),iy+1,igsp)) +
-     .                  fxpy(ix,iy,0)*log(tg(ixp1(ix,iy+1),iy+1,igsp)) )
-                ty1=exp(fxm (ix,iy,1)*log(tg(ixm1(ix,iy+1),iy+1,igsp)) + 
-     .                  fx0 (ix,iy,1)*log(tg(ix           ,iy+1,igsp)) +
-     .                  fxp (ix,iy,1)*log(tg(ixp1(ix,iy+1),iy+1,igsp)) +
-     .                  fxmy(ix,iy,1)*log(tg(ixm1(ix,iy)  ,iy  ,igsp)) +
-     .                  fxpy(ix,iy,1)*log(tg(ixp1(ix,iy)  ,iy  ,igsp)) )
+                ty0=exp(fxm (ix,iy,0)*logtg(ixm1(ix,iy)  ,iy  ,igsp) + 
+     .                  fx0 (ix,iy,0)*logtg(ix           ,iy  ,igsp) +
+     .                  fxp (ix,iy,0)*logtg(ixp1(ix,iy)  ,iy  ,igsp) +
+     .                  fxmy(ix,iy,0)*logtg(ixm1(ix,iy+1),iy+1,igsp) +
+     .                  fxpy(ix,iy,0)*logtg(ixp1(ix,iy+1),iy+1,igsp) )
+                ty1=exp(fxm (ix,iy,1)*logtg(ixm1(ix,iy+1),iy+1,igsp) + 
+     .                  fx0 (ix,iy,1)*logtg(ix           ,iy+1,igsp) +
+     .                  fxp (ix,iy,1)*logtg(ixp1(ix,iy+1),iy+1,igsp) +
+     .                  fxmy(ix,iy,1)*logtg(ixm1(ix,iy)  ,iy  ,igsp) +
+     .                  fxpy(ix,iy,1)*logtg(ixp1(ix,iy)  ,iy  ,igsp) )
               endif
               qtgf = cngfy(igsp) * fgtdy(iy)* sy(ix,iy) * 
      .                           ave(gy(ix,iy)/nu1, gy(ix,iy+1)/nu2) *
@@ -5832,16 +5832,16 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg,
      .                 (ix==ixrb(jx).and.ixmxbcl==1) ) isxyfl = .false.
                enddo
                if (methgx .eq. 6) then  # log interpolation
-               grdnv =(   ( fym (ix,iy,1)*log(ng(ix2,iy1 ,igsp)) +
-     .                      fy0 (ix,iy,1)*log(ng(ix2,iy  ,igsp)) +
-     .                      fyp (ix,iy,1)*log(ng(ix2,iy+1,igsp)) + 
-     .                      fymx(ix,iy,1)*log(ng(ix ,iy1 ,igsp)) +
-     .                      fypx(ix,iy,1)*log(ng(ix, iy+1,igsp)) )
-     .                   -( fym (ix,iy,0)*log(ng(ix ,iy1 ,igsp)) +
-     .                      fy0 (ix,iy,0)*log(ng(ix ,iy  ,igsp)) +
-     .                      fyp (ix,iy,0)*log(ng(ix ,iy+1,igsp)) +
-     .                      fymx(ix,iy,0)*log(ng(ix4,iy1 ,igsp)) + 
-     .                      fypx(ix,iy,0)*log(ng(ix6,iy+1,igsp)) ) )/ 
+               grdnv =(   ( fym (ix,iy,1)*logng(ix2,iy1 ,igsp) +
+     .                      fy0 (ix,iy,1)*logng(ix2,iy  ,igsp) +
+     .                      fyp (ix,iy,1)*logng(ix2,iy+1,igsp) + 
+     .                      fymx(ix,iy,1)*logng(ix ,iy1 ,igsp) +
+     .                      fypx(ix,iy,1)*logng(ix, iy+1,igsp) )
+     .                   -( fym (ix,iy,0)*logng(ix ,iy1 ,igsp) +
+     .                      fy0 (ix,iy,0)*logng(ix ,iy  ,igsp) +
+     .                      fyp (ix,iy,0)*logng(ix ,iy+1,igsp) +
+     .                      fymx(ix,iy,0)*logng(ix4,iy1 ,igsp) + 
+     .                      fypx(ix,iy,0)*logng(ix6,iy+1,igsp) ) )/ 
      .                                                  dxnog(ix,iy)
                elseif (methgx .eq. 7) then  # inverse interpolation
                grdnv =( 1/(fym (ix,iy,1)/ng(ix2,iy1 ,igsp) + 
@@ -5874,12 +5874,12 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg,
      .                           0.5*(nuiz(ix,iy,igsp)+nuiz(ix2,iy,igsp))
               if (methgx .eq. 6) then
                fngxy(ix,iy,igsp) =  exp( 0.5*
-     .                     (log(ng(ix2,iy,igsp))+log(ng(ix,iy,igsp))) )*
-     .                               difgx2*(grdnv/cos(angfx(ix,iy)) -
-     .                     (log(ng(ix2,iy,igsp)) - log(ng(ix,iy,igsp)))*
+     .                     (logng(ix2,iy,igsp)+logng(ix,iy,igsp)) )*
+     .                               difgx2*(grdnv/cosangfx(ix,iy) -
+     .                     (logng(ix2,iy,igsp) - logng(ix,iy,igsp))*
      .                                 gxf(ix,iy) ) * sx(ix,iy)
               else
-               fngxy(ix,iy,igsp) = difgx2*( grdnv/cos(angfx(ix,iy)) -
+               fngxy(ix,iy,igsp) = difgx2*( grdnv/cosangfx(ix,iy) -
      .                             (ng(ix2,iy,igsp) - ng(ix,iy,igsp))*
      .                                 gxf(ix,iy) ) * sx(ix,iy)
               endif
@@ -6159,22 +6159,22 @@ c      write (*,*) "neudifpg"
             vygtan(ix,iy,igsp) = 0.  # vygtan is grad(T) rad vel on x-face
                                      # vygtan) only from thermal force
             if (isnonog .eq. 1 .and. iy .le. ny) then
-                  grdnv =(    ( fym (ix,iy,1)*log(tg(ix2,iy1,igsp)) +
-     .                          fy0 (ix,iy,1)*log(tg(ix2,iy ,igsp)) +
-     .                          fyp (ix,iy,1)*log(tg(ix2,iy2,igsp)) +
-     .                          fymx(ix,iy,1)*log(tg(ix ,iy1,igsp)) +
-     .                          fypx(ix,iy,1)*log(tg(ix, iy2,igsp)) )
-     .                       -( fym (ix,iy,0)*log(tg(ix ,iy1,igsp)) +
-     .                          fy0 (ix,iy,0)*log(tg(ix ,iy ,igsp)) +
-     .                          fyp (ix,iy,0)*log(tg(ix ,iy2,igsp)) + 
-     .                          fymx(ix,iy,0)*log(tg(ix4,iy1,igsp)) +
-     .                          fypx(ix,iy,0)*log(tg(ix6,iy2,igsp)) ) )/ 
+                  grdnv =(    ( fym (ix,iy,1)*logtg(ix2,iy1,igsp) +
+     .                          fy0 (ix,iy,1)*logtg(ix2,iy ,igsp) +
+     .                          fyp (ix,iy,1)*logtg(ix2,iy2,igsp) +
+     .                          fymx(ix,iy,1)*logtg(ix ,iy1,igsp) +
+     .                          fypx(ix,iy,1)*logtg(ix, iy2,igsp) )
+     .                       -( fym (ix,iy,0)*logtg(ix ,iy1,igsp) +
+     .                          fy0 (ix,iy,0)*logtg(ix ,iy ,igsp) +
+     .                          fyp (ix,iy,0)*logtg(ix ,iy2,igsp) + 
+     .                          fymx(ix,iy,0)*logtg(ix4,iy1,igsp) +
+     .                          fypx(ix,iy,0)*logtg(ix6,iy2,igsp) ) )/ 
      .                                                      dxnog(ix,iy)
                vygtan(ix,iy,igsp) = exp( 0.5*
-     .                     (log(tg(ix2,iy,igsp))+log(tg(ix,iy,igsp))) )*
+     .                     (logtg(ix2,iy,igsp)+logtg(ix,iy,igsp)) )*
      .                      ( alftng / (mg(igsp)*0.5*(nu1+nu2)) ) *
-     .                                     ( grdnv/cos(angfx(ix,iy)) - 
-     .                       (log(tg(ix2,iy,igsp)) - log(tg(ix,iy,igsp)))
+     .                                     ( grdnv/cosangfx(ix,iy) - 
+     .                       (logtg(ix2,iy,igsp) - logtg(ix,iy,igsp))
      .                                                 * gxf(ix,iy) ) 
              if (islimon.eq.1.and. ix.eq.ix_lim.and. iy.ge.iy_lims) then
                vygtan(ix,iy,igsp) = 0.
@@ -6279,16 +6279,16 @@ c      write (*,*) "neudifpg"
      .                fxmy(ix,iy,1)*tg(ixm1(ix,iy)  ,iy  ,igsp) +
      .                fxpy(ix,iy,1)*tg(ixp1(ix,iy)  ,iy  ,igsp)
               elseif (isintlog .eq. 1) then
-                ty0=exp(fxm (ix,iy,0)*log(tg(ixm1(ix,iy)  ,iy  ,igsp)) + 
-     .                  fx0 (ix,iy,0)*log(tg(ix           ,iy  ,igsp)) +
-     .                  fxp (ix,iy,0)*log(tg(ixp1(ix,iy)  ,iy  ,igsp)) +
-     .                  fxmy(ix,iy,0)*log(tg(ixm1(ix,iy+1),iy+1,igsp)) +
-     .                  fxpy(ix,iy,0)*log(tg(ixp1(ix,iy+1),iy+1,igsp)) )
-                ty1=exp(fxm (ix,iy,1)*log(tg(ixm1(ix,iy+1),iy+1,igsp)) + 
-     .                  fx0 (ix,iy,1)*log(tg(ix           ,iy+1,igsp)) +
-     .                  fxp (ix,iy,1)*log(tg(ixp1(ix,iy+1),iy+1,igsp)) +
-     .                  fxmy(ix,iy,1)*log(tg(ixm1(ix,iy)  ,iy  ,igsp)) +
-     .                  fxpy(ix,iy,1)*log(tg(ixp1(ix,iy)  ,iy  ,igsp)) )
+                ty0=exp(fxm (ix,iy,0)*logtg(ixm1(ix,iy)  ,iy  ,igsp) + 
+     .                  fx0 (ix,iy,0)*logtg(ix           ,iy  ,igsp) +
+     .                  fxp (ix,iy,0)*logtg(ixp1(ix,iy)  ,iy  ,igsp) +
+     .                  fxmy(ix,iy,0)*logtg(ixm1(ix,iy+1),iy+1,igsp) +
+     .                  fxpy(ix,iy,0)*logtg(ixp1(ix,iy+1),iy+1,igsp) )
+                ty1=exp(fxm (ix,iy,1)*logtg(ixm1(ix,iy+1),iy+1,igsp) + 
+     .                  fx0 (ix,iy,1)*logtg(ix           ,iy+1,igsp) +
+     .                  fxp (ix,iy,1)*logtg(ixp1(ix,iy+1),iy+1,igsp) +
+     .                  fxmy(ix,iy,1)*logtg(ixm1(ix,iy)  ,iy  ,igsp) +
+     .                  fxpy(ix,iy,1)*logtg(ixp1(ix,iy)  ,iy  ,igsp) )
               endif
               qtgf = alftng * fgtdy(iy)* sy(ix,iy) * 
      .                           ave(gy(ix,iy)/nu1, gy(ix,iy+1)/nu2) *
@@ -6414,11 +6414,11 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg,
               if (methgx .eq. 6) then
                fngxy(ix,iy,igsp) =  exp( 0.5*
      .                     (log(pg(ix2,iy,igsp))+log(pg(ix,iy,igsp))) )*
-     .                               difgx2*(grdnv/cos(angfx(ix,iy)) -
+     .                               difgx2*(grdnv/cosangfx(ix,iy) -
      .                     (log(pg(ix2,iy,igsp)) - log(pg(ix,iy,igsp)))*
      .                                 gxf(ix,iy) ) * sx(ix,iy)
               else
-               fngxy(ix,iy,igsp) = difgx2*(grdnv/cos(angfx(ix,iy)) -
+               fngxy(ix,iy,igsp) = difgx2*(grdnv/cosangfx(ix,iy) -
      .                             (pg(ix2,iy,igsp) - pg(ix,iy,igsp))*
      .                                 gxf(ix,iy) ) * sx(ix,iy)
               endif
@@ -6735,23 +6735,23 @@ c   physically meaningful gas variables (flux, velocity, etc)
             endif
             vygtan(ix,iy,igsp) = 0.  # vygtan is grad(T) rad vel on x-face
             if (isnonog .eq. 1 .and. iy .le. ny) then
-                  grdnv =(    ( fym (ix,iy,1)*log(tg(ix2,iy1,igsp)) +
-     .                          fy0 (ix,iy,1)*log(tg(ix2,iy ,igsp)) +
-     .                          fyp (ix,iy,1)*log(tg(ix2,iy2,igsp)) +
-     .                          fymx(ix,iy,1)*log(tg(ix ,iy1,igsp)) +
-     .                          fypx(ix,iy,1)*log(tg(ix, iy2,igsp)) )
-     .                       -( fym (ix,iy,0)*log(tg(ix ,iy1,igsp)) +
-     .                          fy0 (ix,iy,0)*log(tg(ix ,iy ,igsp)) +
-     .                          fyp (ix,iy,0)*log(tg(ix ,iy2,igsp)) + 
-     .                          fymx(ix,iy,0)*log(tg(ix4,iy1,igsp)) +
-     .                          fypx(ix,iy,0)*log(tg(ix6,iy2,igsp)) ) )/ 
+                  grdnv =(    ( fym (ix,iy,1)*logtg(ix2,iy1,igsp) +
+     .                          fy0 (ix,iy,1)*logtg(ix2,iy ,igsp) +
+     .                          fyp (ix,iy,1)*logtg(ix2,iy2,igsp) +
+     .                          fymx(ix,iy,1)*logtg(ix ,iy1,igsp) +
+     .                          fypx(ix,iy,1)*logtg(ix, iy2,igsp) )
+     .                       -( fym (ix,iy,0)*logtg(ix ,iy1,igsp) +
+     .                          fy0 (ix,iy,0)*logtg(ix ,iy ,igsp) +
+     .                          fyp (ix,iy,0)*logtg(ix ,iy2,igsp) + 
+     .                          fymx(ix,iy,0)*logtg(ix4,iy1,igsp) +
+     .                          fypx(ix,iy,0)*logtg(ix6,iy2,igsp) ) )/ 
      .                                                      dxnog(ix,iy)
                vygtan(ix,iy,igsp) = exp( 0.5*
-     .                     (log(tg(ix2,iy,igsp))+log(tg(ix,iy,igsp))) )*
+     .                     (logtg(ix2,iy,igsp)+logtg(ix,iy,igsp)) )*
      .                                  ( cngfx(igsp) / (mg(igsp)*0.5*
      .                         (nuix(ix,iy,igsp)+nuix(ix2,iy,igsp))) ) *
-     .                             ( grdnv/cos(angfx(ix,iy)) - 
-     .                       (log(tg(ix2,iy,igsp)) - log(tg(ix,iy,igsp)))
+     .                             ( grdnv/cosangfx(ix,iy) - 
+     .                       (logtg(ix2,iy,igsp) - logtg(ix,iy,igsp))
      .                                                 * gxf(ix,iy) ) 
                if (islimon.eq.1.and. ix.eq.ix_lim.and. iy.ge.iy_lims) then
                   vygtan(ix,iy,igsp) = 0.
@@ -6766,7 +6766,7 @@ c   physically meaningful gas variables (flux, velocity, etc)
             if (isupgon(igsp) .eq. 1) then
               qtgf = qtgf + rrv(ix,iy)*up(ix,iy,iigsp)*sx(ix,iy)
             endif
-            qsh = csh * (lng(ix,iy,igsp)-lng(ix2,iy,igsp)) + qtgf
+            qsh = csh * (logng(ix,iy,igsp)-logng(ix2,iy,igsp)) + qtgf
             qr = abs(qsh/qfl)
 c...  Because guard-cell values may be distorted from B.C., possibly omit terms on
 c...  boundary face - shouldnt matter(just set BC) except for guard-cell values
@@ -6785,7 +6785,7 @@ c   physically meaningful gas variables (flux, velocity, etc)
 c...  now add the convective velocity for charge-exchange neutrals
          if(igsp .eq. 1) floxg(ix,iy) = 
      .              floxg(ix,iy) + cngflox(1)*sx(ix,iy)*uu(ix,iy,1)
-         floxg(ix,iy) = floxg(ix,iy)*2/(lng(ix,iy,igsp)+lng(ix2,iy,igsp))
+         floxg(ix,iy) = floxg(ix,iy)*2/(logng(ix,iy,igsp)+logng(ix2,iy,igsp))
 
   887    continue
          conxg(nx+1,iy) = 0
@@ -6825,16 +6825,16 @@ c   physically meaningful gas variables (flux, velocity, etc)
      .                fxmy(ix,iy,1)*tg(ixm1(ix,iy)  ,iy  ,igsp) +
      .                fxpy(ix,iy,1)*tg(ixp1(ix,iy)  ,iy  ,igsp)
               elseif (isintlog .eq. 1) then
-                ty0=exp(fxm (ix,iy,0)*log(tg(ixm1(ix,iy)  ,iy  ,igsp)) + 
-     .                  fx0 (ix,iy,0)*log(tg(ix           ,iy  ,igsp)) +
-     .                  fxp (ix,iy,0)*log(tg(ixp1(ix,iy)  ,iy  ,igsp)) +
-     .                  fxmy(ix,iy,0)*log(tg(ixm1(ix,iy+1),iy+1,igsp)) +
-     .                  fxpy(ix,iy,0)*log(tg(ixp1(ix,iy+1),iy+1,igsp)) )
-                ty1=exp(fxm (ix,iy,1)*log(tg(ixm1(ix,iy+1),iy+1,igsp)) + 
-     .                  fx0 (ix,iy,1)*log(tg(ix           ,iy+1,igsp)) +
-     .                  fxp (ix,iy,1)*log(tg(ixp1(ix,iy+1),iy+1,igsp)) +
-     .                  fxmy(ix,iy,1)*log(tg(ixm1(ix,iy)  ,iy  ,igsp)) +
-     .                  fxpy(ix,iy,1)*log(tg(ixp1(ix,iy)  ,iy  ,igsp )) )
+                ty0=exp(fxm (ix,iy,0)*logtg(ixm1(ix,iy)  ,iy  ,igsp) + 
+     .                  fx0 (ix,iy,0)*logtg(ix           ,iy  ,igsp) +
+     .                  fxp (ix,iy,0)*logtg(ixp1(ix,iy)  ,iy  ,igsp) +
+     .                  fxmy(ix,iy,0)*logtg(ixm1(ix,iy+1),iy+1,igsp) +
+     .                  fxpy(ix,iy,0)*logtg(ixp1(ix,iy+1),iy+1,igsp) )
+                ty1=exp(fxm (ix,iy,1)*logtg(ixm1(ix,iy+1),iy+1,igsp) + 
+     .                  fx0 (ix,iy,1)*logtg(ix           ,iy+1,igsp) +
+     .                  fxp (ix,iy,1)*logtg(ixp1(ix,iy+1),iy+1,igsp) +
+     .                  fxmy(ix,iy,1)*logtg(ixm1(ix,iy)  ,iy  ,igsp) +
+     .                  fxpy(ix,iy,1)*logtg(ixp1(ix,iy)  ,iy  ,igsp ) )
               endif
               qtgf = cngfy(igsp) * fgtdy(iy) * sy(ix,iy) * 
      .                      ave( gy(ix,iy)/nuix(ix,iy,igsp) ,
@@ -6847,7 +6847,7 @@ c   physically meaningful gas variables (flux, velocity, etc)
             if(methgy.ne.2) nconv =
      .                         ngy0(ix,iy,igsp)*0.5*(1+sign(1.,qtgf)) +
      .                         ngy1(ix,iy,igsp)*0.5*(1-sign(1.,qtgf))
-            qsh = csh * (lng(ix,iy,igsp)-lng(ix,iy+1,igsp)) + qtgf
+            qsh = csh * (logng(ix,iy,igsp)-logng(ix,iy+1,igsp)) + qtgf
             qr = abs(qsh/qfl)
             if(iy.eq.0 .and. iymnbcl.eq.1) then
                qr = gcfacgy*qr
@@ -6866,7 +6866,7 @@ c   physically meaningful gas variables (flux, velocity, etc)
          if(igsp .eq. 1) floyg(ix,iy) =  
      .             floyg(ix,iy)+cngfloy(1)*sy(ix,iy)*vy(ix,iy,1)
 
-         floyg(ix,iy)=floyg(ix,iy)*2/(lng(ix,iy,igsp)+lng(ix,iy+1,igsp))
+         floyg(ix,iy)=floyg(ix,iy)*2/(logng(ix,iy,igsp)+logng(ix,iy+1,igsp))
   889    continue
   890 continue
 
@@ -6875,7 +6875,7 @@ c   physically meaningful gas variables (flux, velocity, etc)
 *  --------------------------------------------------------------------
 
       call fd2tra (nx,ny,floxg,floyg,conxg,conyg,
-     .             lng(0:nx+1,0:ny+1,igsp),flngx(0:nx+1,0:ny+1,igsp),
+     .             logng(0:nx+1,0:ny+1,igsp),flngx(0:nx+1,0:ny+1,igsp),
      .             flngy(0:nx+1,0:ny+1,igsp),0,methg)
 
 c...  Addition for nonorthogonal mesh
@@ -6896,16 +6896,16 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg,
                   if ( (ix==ixlb(jx).and.ixmnbcl==1) .or.
      .                 (ix==ixrb(jx).and.ixmxbcl==1) )isxyfl = .false.
                enddo
-               grdnv =( (fym (ix,iy,1)*lng(ix2,iy1 ,igsp) + 
-     .                   fy0 (ix,iy,1)*lng(ix2,iy  ,igsp) +
-     .                   fyp (ix,iy,1)*lng(ix2,iy+1,igsp) + 
-     .                   fymx(ix,iy,1)*lng(ix ,iy1 ,igsp) +
-     .                   fypx(ix,iy,1)*lng(ix, iy+1,igsp))
-     .                - (fym (ix,iy,0)*lng(ix ,iy1 ,igsp) +
-     .                   fy0 (ix,iy,0)*lng(ix ,iy  ,igsp) +
-     .                   fyp (ix,iy,0)*lng(ix ,iy+1,igsp) +
-     .                   fymx(ix,iy,0)*lng(ix4,iy1 ,igsp) + 
-     .                   fypx(ix,iy,0)*lng(ix6,iy+1,igsp)) ) / 
+               grdnv =( (fym (ix,iy,1)*logng(ix2,iy1 ,igsp) + 
+     .                   fy0 (ix,iy,1)*logng(ix2,iy  ,igsp) +
+     .                   fyp (ix,iy,1)*logng(ix2,iy+1,igsp) + 
+     .                   fymx(ix,iy,1)*logng(ix ,iy1 ,igsp) +
+     .                   fypx(ix,iy,1)*logng(ix, iy+1,igsp))
+     .                - (fym (ix,iy,0)*logng(ix ,iy1 ,igsp) +
+     .                   fy0 (ix,iy,0)*logng(ix ,iy  ,igsp) +
+     .                   fyp (ix,iy,0)*logng(ix ,iy+1,igsp) +
+     .                   fymx(ix,iy,0)*logng(ix4,iy1 ,igsp) + 
+     .                   fypx(ix,iy,0)*logng(ix6,iy+1,igsp)) ) / 
      .                                                  dxnog(ix,iy)
 
                difgx2 = ave( tg(ix ,iy,igsp)/nuix(ix ,iy,igsp),
@@ -6913,8 +6913,8 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg,
      .                         + rld2dxg(igsp)**2*(1/gxf(ix,iy)**2)*
      .                           0.5*(nuiz(ix,iy,igsp)+nuiz(ix2,iy,igsp))
 
-               flngxy(ix,iy,igsp) = difgx2*(grdnv/cos(angfx(ix,iy)) -
-     .                             (lng(ix2,iy,igsp) - lng(ix,iy,igsp))*
+               flngxy(ix,iy,igsp) = difgx2*(grdnv/cosangfx(ix,iy) -
+     .                             (logng(ix2,iy,igsp) - logng(ix,iy,igsp))*
      .                                 gxf(ix,iy) ) * sx(ix,iy)
 
 c...  Now flux limit with flalfgxy
@@ -6977,14 +6977,14 @@ cc               uu(ix,iy,iigsp) = uug(ix,iy,igsp)
          do ix = i1, i5
             ix2 = ixp1(ix,iy)
             fngx(ix,iy,igsp) = flngx(ix,iy,igsp) *
-     .                     exp(0.5*(lng(ix,iy,igsp)+lng(ix2,iy,igsp)))
+     .                     exp(0.5*(logng(ix,iy,igsp)+logng(ix2,iy,igsp)))
          enddo
       enddo
 c ...   now do fniy
       do iy = j1, j5    # same loop ranges as for fngy in fd2tra
          do ix = i4, i8
             fngy(ix,iy,igsp) = flngy(ix,iy,igsp) *
-     .                     exp(0.5*(lng(ix,iy,igsp)+lng(ix,iy+1,igsp)))
+     .                     exp(0.5*(logng(ix,iy,igsp)+logng(ix,iy+1,igsp)))
          enddo
       enddo
 
@@ -7147,21 +7147,21 @@ cc               uu(ix,iy,iigsp) = uug(ix,iy,igsp)
      .                     fymx(ix,iy,0)*tg(ix4,iy1,igsp) -
      .                     fypx(ix,iy,0)*tg(ix6,iy2,igsp) )/dxnog(ix,iy)
                elseif (isintlog .eq. 1) then
-                  grdnv =( exp( fym (ix,iy,1)*log(tg(ix2,iy1,igsp)) +
-     .                          fy0 (ix,iy,1)*log(tg(ix2,iy ,igsp)) +
-     .                          fyp (ix,iy,1)*log(tg(ix2,iy2,igsp)) +
-     .                          fymx(ix,iy,1)*log(tg(ix ,iy1,igsp)) +
-     .                          fypx(ix,iy,1)*log(tg(ix, iy2,igsp)) )
-     .                    -exp( fym (ix,iy,0)*log(tg(ix ,iy1,igsp)) +
-     .                          fy0 (ix,iy,0)*log(tg(ix ,iy ,igsp)) +
-     .                          fyp (ix,iy,0)*log(tg(ix ,iy2,igsp)) + 
-     .                          fymx(ix,iy,0)*log(tg(ix4,iy1,igsp)) +
-     .                          fypx(ix,iy,0)*log(tg(ix6,iy2,igsp)) ) )/ 
+                  grdnv =( exp( fym (ix,iy,1)*logtg(ix2,iy1,igsp) +
+     .                          fy0 (ix,iy,1)*logtg(ix2,iy ,igsp) +
+     .                          fyp (ix,iy,1)*logtg(ix2,iy2,igsp) +
+     .                          fymx(ix,iy,1)*logtg(ix ,iy1,igsp) +
+     .                          fypx(ix,iy,1)*logtg(ix, iy2,igsp) )
+     .                    -exp( fym (ix,iy,0)*logtg(ix ,iy1,igsp) +
+     .                          fy0 (ix,iy,0)*logtg(ix ,iy ,igsp) +
+     .                          fyp (ix,iy,0)*logtg(ix ,iy2,igsp) + 
+     .                          fymx(ix,iy,0)*logtg(ix4,iy1,igsp) +
+     .                          fypx(ix,iy,0)*logtg(ix6,iy2,igsp) ) )/ 
      .                                                      dxnog(ix,iy)
                endif
                vygtan(ix,iy,igsp) = ( cngfx(igsp) / (mg(igsp)*0.5*
      .                         (nuix(ix,iy,igsp)+nuix(ix2,iy,igsp))) ) *
-     .                             ( grdnv/cos(angfx(ix,iy)) - 
+     .                             ( grdnv/cosangfx(ix,iy) - 
      .                             (tg(ix2,iy,igsp) - tg(ix,iy,igsp))
      .                                                 * gxf(ix,iy) ) 
                if (islimon.eq.1.and. ix.eq.ix_lim.and. iy.ge.iy_lims) then
@@ -7241,16 +7241,16 @@ cc               uu(ix,iy,iigsp) = uug(ix,iy,igsp)
      .                fxmy(ix,iy,1)*tg(ixm1(ix,iy)  ,iy  ,igsp) +
      .                fxpy(ix,iy,1)*tg(ixp1(ix,iy)  ,iy  ,igsp)
               elseif (isintlog .eq. 1) then
-                ty0=exp(fxm (ix,iy,0)*log(tg(ixm1(ix,iy)  ,iy  ,igsp)) + 
-     .                  fx0 (ix,iy,0)*log(tg(ix           ,iy  ,igsp)) +
-     .                  fxp (ix,iy,0)*log(tg(ixp1(ix,iy)  ,iy  ,igsp)) +
-     .                  fxmy(ix,iy,0)*log(tg(ixm1(ix,iy+1),iy+1,igsp)) +
-     .                  fxpy(ix,iy,0)*log(tg(ixp1(ix,iy+1),iy+1,igsp)) )
-                ty1=exp(fxm (ix,iy,1)*log(tg(ixm1(ix,iy+1),iy+1,igsp)) + 
-     .                  fx0 (ix,iy,1)*log(tg(ix           ,iy+1,igsp)) +
-     .                  fxp (ix,iy,1)*log(tg(ixp1(ix,iy+1),iy+1,igsp)) +
-     .                  fxmy(ix,iy,1)*log(tg(ixm1(ix,iy)  ,iy  ,igsp)) +
-     .                  fxpy(ix,iy,1)*log(tg(ixp1(ix,iy)  ,iy  ,igsp)) )
+                ty0=exp(fxm (ix,iy,0)*logtg(ixm1(ix,iy)  ,iy  ,igsp) + 
+     .                  fx0 (ix,iy,0)*logtg(ix           ,iy  ,igsp) +
+     .                  fxp (ix,iy,0)*logtg(ixp1(ix,iy)  ,iy  ,igsp) +
+     .                  fxmy(ix,iy,0)*logtg(ixm1(ix,iy+1),iy+1,igsp) +
+     .                  fxpy(ix,iy,0)*logtg(ixp1(ix,iy+1),iy+1,igsp) )
+                ty1=exp(fxm (ix,iy,1)*logtg(ixm1(ix,iy+1),iy+1,igsp) + 
+     .                  fx0 (ix,iy,1)*logtg(ix           ,iy+1,igsp) +
+     .                  fxp (ix,iy,1)*logtg(ixp1(ix,iy+1),iy+1,igsp) +
+     .                  fxmy(ix,iy,1)*logtg(ixm1(ix,iy)  ,iy  ,igsp) +
+     .                  fxpy(ix,iy,1)*logtg(ixp1(ix,iy)  ,iy  ,igsp) )
               endif
               qtgf = cngfy(igsp) * fgtdy(iy) * sy(ix,iy) * 
      .                      ave( gy(ix,iy)/nuix(ix,iy,igsp) ,
@@ -7340,16 +7340,16 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg,
                ix5 = ixm1(ix,iy+1)
                ix6 = ixp1(ix,iy+1) 
                if (methgx .eq. 6) then  # log interpolation
-               grdnv =( exp(fym (ix,iy,1)*log(ng(ix2,iy1 ,igsp)) + 
-     .                      fy0 (ix,iy,1)*log(ng(ix2,iy  ,igsp)) +
-     .                      fyp (ix,iy,1)*log(ng(ix2,iy+1,igsp)) + 
-     .                      fymx(ix,iy,1)*log(ng(ix ,iy1 ,igsp)) +
-     .                      fypx(ix,iy,1)*log(ng(ix, iy+1,igsp)))
-     .                - exp(fym (ix,iy,0)*log(ng(ix ,iy1 ,igsp)) +
-     .                      fy0 (ix,iy,0)*log(ng(ix ,iy  ,igsp)) +
-     .                      fyp (ix,iy,0)*log(ng(ix ,iy+1,igsp)) +
-     .                      fymx(ix,iy,0)*log(ng(ix4,iy1 ,igsp)) + 
-     .                      fypx(ix,iy,0)*log(ng(ix6,iy+1,igsp))) ) / 
+               grdnv =( exp(fym (ix,iy,1)*logng(ix2,iy1 ,igsp) + 
+     .                      fy0 (ix,iy,1)*logng(ix2,iy  ,igsp) +
+     .                      fyp (ix,iy,1)*logng(ix2,iy+1,igsp) + 
+     .                      fymx(ix,iy,1)*logng(ix ,iy1 ,igsp) +
+     .                      fypx(ix,iy,1)*logng(ix, iy+1,igsp) )
+     .                - exp(fym (ix,iy,0)*logng(ix ,iy1 ,igsp) +
+     .                      fy0 (ix,iy,0)*logng(ix ,iy  ,igsp) +
+     .                      fyp (ix,iy,0)*logng(ix ,iy+1,igsp) +
+     .                      fymx(ix,iy,0)*logng(ix4,iy1 ,igsp) + 
+     .                      fypx(ix,iy,0)*logng(ix6,iy+1,igsp)) ) / 
      .                                                  dxnog(ix,iy)
                elseif (methgx .eq. 7) then  # inverse interpolation
                grdnv =( 1/(fym (ix,iy,1)/ng(ix2,iy1 ,igsp) + 
@@ -7380,7 +7380,7 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg,
      .                       tg(ix2,iy,igsp)/nuix(ix2,iy,igsp) )/mg(igsp)
      .                         + rld2dxg(igsp)**2*(1/gxf(ix,iy)**2)*
      .                           0.5*(nuiz(ix,iy,igsp)+nuiz(ix2,iy,igsp))
-               fngxy(ix,iy,igsp) = difgx2*(grdnv/cos(angfx(ix,iy)) -
+               fngxy(ix,iy,igsp) = difgx2*(grdnv/cosangfx(ix,iy) -
      .                             (ng(ix2,iy,igsp) - ng(ix,iy,igsp))*
      .                                 gxf(ix,iy) ) * sx(ix,iy)
  8904       continue
@@ -7741,16 +7741,16 @@ call fd2tra (nx,ny,floxge(0:nx+1,0:ny+1,igsp),
 c --- Note: this four-point average results in not getting the full Jac. for
 c --- a nonorthogonal mesh because of ngy1,0 - see def. of hcyn
 
-                grdnv =( ( fym (ix,iy,1)*log(tg(ix2,iy1 ,igsp)) +  
-     .                     fy0 (ix,iy,1)*log(tg(ix2,iy  ,igsp)) +
-     .                     fyp (ix,iy,1)*log(tg(ix2,iy+1,igsp)) +  
-     .                     fymx(ix,iy,1)*log(tg(ix ,iy1 ,igsp)) +
-     .                     fypx(ix,iy,1)*log(tg(ix ,iy+1,igsp)) ) 
-     .                  -( fym (ix,iy,0)*log(tg(ix ,iy1 ,igsp)) +
-     .                     fy0 (ix,iy,0)*log(tg(ix ,iy  ,igsp)) +
-     .                     fyp (ix,iy,0)*log(tg(ix ,iy+1,igsp)) +
-     .                     fymx(ix,iy,0)*log(tg(ix4,iy1 ,igsp)) +  
-     .                     fypx(ix,iy,0)*log(tg(ix6,iy+1,igsp)) ) ) / 
+                grdnv =( ( fym (ix,iy,1)*logtg(ix2,iy1 ,igsp) +  
+     .                     fy0 (ix,iy,1)*logtg(ix2,iy  ,igsp) +
+     .                     fyp (ix,iy,1)*logtg(ix2,iy+1,igsp) +  
+     .                     fymx(ix,iy,1)*logtg(ix ,iy1 ,igsp) +
+     .                     fypx(ix,iy,1)*logtg(ix ,iy+1,igsp) ) 
+     .                  -( fym (ix,iy,0)*logtg(ix ,iy1 ,igsp) +
+     .                     fy0 (ix,iy,0)*logtg(ix ,iy  ,igsp) +
+     .                     fyp (ix,iy,0)*logtg(ix ,iy+1,igsp) +
+     .                     fymx(ix,iy,0)*logtg(ix4,iy1 ,igsp) +  
+     .                     fypx(ix,iy,0)*logtg(ix6,iy+1,igsp) ) ) / 
      .                                                  dxnog(ix,iy)  
                difgx2 = ave( tg(ix ,iy,igsp)/nu1,
      .                       tg(ix2,iy,igsp)/nu2 )/mg(igsp)
@@ -7758,10 +7758,10 @@ call fd2tra (nx,ny,floxge(0:nx+1,0:ny+1,igsp),
      .                           0.5*(nuiz(ix,iy,igsp)+nuiz(ix2,iy,igsp))
 
                fegxy(ix,iy,igsp) = cfegxy*exp( 0.5*
-     .                 (log(tg(ix2,iy,igsp))+log(tg(ix,iy,igsp))) )*
+     .                 (logtg(ix2,iy,igsp)+logtg(ix,iy,igsp)) )*
      .                      difgx2*ave(ng(ix2,iy,igsp),ng(ix,iy,igsp))*
-     .                                 ( grdnv/cos(angfx(ix,iy))
-     .                  - (log(tg(ix2,iy,igsp)) - log(tg(ix,iy,igsp)))*
+     .                                 ( grdnv/cosangfx(ix,iy)
+     .                  - (logtg(ix2,iy,igsp) - logtg(ix,iy,igsp))*
      .                                        gxf(ix,iy) )*sx(ix,iy)
 c...  Flux limit with flalftxt even though hcys have parallel FL built in
                t0 = max(tg(ix,iy,igsp),tgmin*ev)
@@ -10586,7 +10586,7 @@ call s2fill (nx+2, ny+2, 0., prdu, 1, nx+2)
      .                            rm(ixv+nj,iyv,3)-rm(ix+nj,iy,0) )
               dthgy = abs(theta_ray1-theta_ray2)
               frth = min(dthgy, 2*pi-dthgy)/(2*pi)  # frac.; need angle < pi
-              sxo = sx(ixv,iyv)/(cos(angfx(ixv,iyv)))
+              sxo = sx(ixv,iyv)/(cosangfx(ixv,iyv))
               pwr_pltz(iyv,ip) = pwr_pltz(iyv,ip) + 
      .                                 prdu(ix,iy)*vol(ix,iy)*frth/sxo
               pwr_plth(iyv,ip) = pwr_plth(iyv,ip) + (
@@ -10707,8 +10707,8 @@ real sxi(0:ny+1,1:nxpt),sxo(0:ny+1,1:nxpt)
         ixi=ixlb(jx)	# ixi=0
         ixo=ixrb(jx)	# ixo=nx
         do iy=1,ny+1
-          sxo(iy,jx) = sx(ixo,iy)/(cos(angfx(ixo,iy)))
-          sxi(iy,jx) = sx(ixi,iy)/(cos(angfx(ixi,iy)))
+          sxo(iy,jx) = sx(ixo,iy)/(cosangfx(ixo,iy))
+          sxi(iy,jx) = sx(ixi,iy)/(cosangfx(ixi,iy))
           do id = 1, nfsp
 	    gdilb(iy,id,jx) = -fnix(ixi,iy,id)/sxi(iy,jx)
 	    gdirb(iy,id,jx) =  fnix(ixo,iy,id)/sxo(iy,jx)
@@ -11970,149 +11970,6 @@ real ps_tmp(0:n1+1,0:n2+1), ps(0:n1+1,0:n2+1)
       return
       end
 c ======================================================================
-c
-      subroutine coneq
-
-c ... This subroutine calculates the fluxes needed for the ion continuity
-c ... equations
-
-      implicit none
-
-*  -- local variables
-      integer methnx,methny,ifld
-      #Former Aux module variables
-      integer ix,iy,igsp,iv,iv1,iv2,iv3,ix2,ix3,ix4,ix5,ix6
-      real t,t0,t1,t2,a
-      Use(Dim)      # nx,ny,nhsp,nisp,ngsp
-      Use(Share)    # nxpt,geometry,nxc,cutlo
-      Use(Xpoint_indices)      # ixpt1,ixpt2,iysptrx
-      Use(UEpar)    # methn,nlimix,nlimiy,nlimiy
-      Use(Coefeq)   # cnfx,cnfy
-      Use(Selec)    # i1,i2,i3,i4,i5,i6,i7,i8,j1,j2,j3,j4,j5,j6,j7,j8
-                    # xlinc,xrinc,yinc,ixm1,ixp1
-      Use(Comgeo)   # vol, gx, gy, sx ,xy
-      Use(Compla)   # mi, zi, ni, uu, up, v2, v2ce, vygtan, mg
-      Use(Comflo)   # flnix,flniy
-      Use(Indices_domain_dcl)   # ixmxbcl
-	  
-c ------------------
-
-      do 104 ifld = 1, nfsp
-*  ---------------------------------------------------------------------
-*     compute flux, residual
-*     The residual is: res := snic + sniv * ni - outflow(ni).
-*  ---------------------------------------------------------------------
-
-*  -- compute flnix, flox, conx --
-
-         methnx = mod(methn, 10)
-         methny = methn/10
-         do 81 iy = j4, j8
-            do 80 ix = i1, i5
-              if (zi(ifld).eq.0. .and. ineudif <= 2) then
-                 flnix(ix,iy,ifld) = fngx(ix,iy,1)
-              else
-               ix2 = ixp1(ix,iy)
-
-               if (methnx .eq. 2) then   # central differencing
-
-                  t2 = ( lni(ix, iy,ifld) + lni(ix2,iy,ifld) ) / 2
-
-               elseif (methnx .eq. 3) then   # upwind differencing
-
-                  if( uu(ix,iy,ifld) .ge. 0.) then
-                     t2 = lni(ix,iy,ifld)
-                  else
-                     t2 = lni(ix2,iy,ifld)
-                  endif
-
-               else   # interp. ave or harmonic ave depending on wind*grad
-
-                  t0 = ( lni(ix, iy,ifld)*gx(ix, iy) +
-     .                   lni(ix2,iy,ifld)*gx(ix2,iy) ) / 
-     .                                      ( gx(ix,iy)+gx(ix2,iy) )
-                  t1 = ( gx(ix,iy)+gx(ix2,iy) ) * lni(ix,iy,ifld) *
-     .                   lni(ix2,iy,ifld) / ( cutlo + lni(ix,iy,ifld)*
-     .                   gx(ix2,iy) + lni(ix2,iy,ifld)*gx(ix,iy) )
-                  if( uu(ix,iy,ifld)*(lni(ix,iy,ifld)-lni(ix2,iy,ifld)) 
-     .                                                     .ge. 0.) then
-                     t2 = t0
-                  else
-                     t2 = t1
-                  endif
-
-               endif
-
-               flnix(ix,iy,ifld) = cnfx*uu(ix,iy,ifld) * sx(ix,iy) * t2
-                  flnix(ix,iy,ifld) = flnix(ix,iy,ifld)/( 1 +
-     .                             (nlimix(ifld)/lni(ix2,iy,ifld))**2 +
-     .                             (nlimix(ifld)/lni(ix ,iy,ifld))**2 )
-               endif
-   80      continue
-           if ((isudsym==1.or.geometry.eq.'dnXtarget') 
-     .                               .and. nxc > 1) flnix(nxc,iy,ifld)=0.
-           if (islimon.ne.0 .and. iy.ge.iy_lims) flnix(ix_lim,iy,ifld)=0.
-           if (nxpt==2.and.ixmxbcl==1) flnix(ixrb(1)+1,iy,ifld)=0.
-   81    continue
-
-
-*  -- compute flniy, floy, cony --
-
-         do 83 iy = j1, j5
-            do 82 ix = i4, i8
-               if (zi(ifld).eq.0.) then
-                  flniy(ix,iy,ifld) = fngy(ix,iy,1)
-               else
-                  if (methny .eq. 2) then   # central differencing
-
-                     t2 = ( niy0(ix,iy,ifld) + niy1(ix,iy,ifld) ) / 2
-
-                  elseif (methny .eq. 3) then    # upwind differencing
-
-                     if( vy(ix,iy,ifld) .ge. 0.) then
-                        t2 = niy0(ix,iy,ifld)
-                     else
-                        t2 = niy1(ix,iy,ifld)
-                     endif
-
-                  else    # interp. ave or harmonic ave depending on wind*grad
-
-                     t0 = ( niy0(ix,iy,ifld)*gy(ix,iy  ) +
-     .                    niy1(ix,iy,ifld)*gy(ix,iy+1) ) / 
-     .                    ( gy(ix,iy)+gy(ix,iy+1) )
-                     t1 = ( gy(ix,iy)+gy(ix,iy+1) ) * niy0(ix,iy,ifld)*
-     .                    niy1(ix,iy,ifld) / ( cutlo + niy0(ix,iy,ifld)*
-     .                    gy(ix,iy+1) + niy1(ix,iy,ifld)*gy(ix,iy) )
-                     if( (niy0(ix,iy,ifld)-niy1(ix,iy,ifld))*
-     .                    vy(ix,iy,ifld) .ge. 0.) then
-                        t2 = t0
-                     else
-                        t2 = t1
-                     endif
-                  
-                  endif
-               
-                  flniy(ix,iy,ifld) = cnfy*vy(ix,iy,ifld)*sy(ix,iy)*t2
-                  if (vy(ix,iy,ifld)*(lni(ix,iy,ifld)-lni(ix,iy+1,ifld))
-     .                                                      .lt. 0.) then
-                     flniy(ix,iy,ifld) = flniy(ix,iy,ifld)/( 1 +
-     .                               (nlimiy(ifld)/lni(ix,iy+1,ifld))**2 +
-     .                               (nlimiy(ifld)/lni(ix,iy  ,ifld))**2 )
-                  endif
-               endif
- 82         continue
- 83      continue
-         
-         do ix = i4, i8
-            flniy(ix,ny+1,ifld) = 0.0e0
-         enddo
-
- 104  continue
-
-      return
-      end
-c ***** End of subroutine coneq **********
-c ======================================================================
 c
       subroutine upvisneo
 
diff --git a/bbb/odesetup.m b/bbb/odesetup.m
index 2641db6f..8b2fed81 100755
--- a/bbb/odesetup.m
+++ b/bbb/odesetup.m
@@ -1606,8 +1606,9 @@ call remark('*** Error: ngs <= 0 ***')
                       write (*,*) 'Error at ix=', ix,'  iy=',iy
                       call xerrab("")
                    endif
-                   lng(ix,iy,igsp) = log(ng(ix,iy,igsp))
+                   logng(ix,iy,igsp)=LOG(ABS(ng(ix,iy,igsp)))
                    tg(ix,iy,igsp) = tgs(ix,iy,igsp)
+                   logtg(ix,iy,igsp) = LOG(tg(ix,iy,igsp))
                 enddo
                 te(ix,iy)      = tes(ix,iy)
                 ti(ix,iy)      = tis(ix,iy)
@@ -1772,7 +1773,7 @@ call intpvar (tgs(0:,0:,igsp), tg(0:,0:,igsp), 0, nxold, nyold)
                   ng(ix,iy,igsp) = max(ng(ix,iy,igsp),
      .                                 1.0e-01*ngbackg(igsp))
                   if (ineudif .eq. 2) then
-                    lng(ix,iy,igsp) = log(ng(ix,iy,igsp))
+                    logng(ix,iy,igsp) = log(ng(ix,iy,igsp))
                   endif
                endif
              enddo
diff --git a/com/com.v b/com/com.v
index 13bf3388..476e584c 100755
--- a/com/com.v
+++ b/com/com.v
@@ -455,6 +455,7 @@ lconeg(0:nx+1,0:ny+1) _real [m]   #glob banana-wid corrected elec conn leng
 vtag(0:nx+1,0:ny+1)    _real  [ ]  #angle at upper-right vertex; 0 rad. is orthg
                                    #clockwise rot. of face => vtag > 0
 angfx(0:nx+1,0:ny+1)   _real  [ ]  #angle on x-face;=.5*(vtag(,iy)+vtag(,iy-1))
+cosangfx(0:nx+1,0:ny+1) _real [ ]  # cos of angfx
 fxm(0:nx+1,0:ny+1,0:1) _real  [ ]  #frac. of (ix-1,iy+k) pt used for y deriv,ave
                                    # for (ix,iy) cell where k is third index
 fx0(0:nx+1,0:ny+1,0:1) _real  [ ]  #frac. of (ix,iy+k) pt used for y deriv,ave
diff --git a/ppp/debug_parallel.F90 b/ppp/debug_parallel.F90
index 82a647d1..4f32443b 100644
--- a/ppp/debug_parallel.F90
+++ b/ppp/debug_parallel.F90
@@ -443,8 +443,8 @@ subroutine DebugHelper(FileName)
       call WriteArrayReal(kye_use,size(kye_use),iunit)
       write(iunit,*) "kyi_use"
       call WriteArrayReal(kyi_use,size(kyi_use),iunit)
-      write(iunit,*) "lng"
-      call WriteArrayReal(lng,size(lng),iunit)
+      write(iunit,*) "logng"
+      call WriteArrayReal(logng,size(logng),iunit)
       write(iunit,*) "loglambda"
       call WriteArrayReal(loglambda,size(loglambda),iunit)
       write(iunit,*) "msor"