diff --git a/PROBLEMS/MONOPOLE_2D/define.h b/PROBLEMS/MONOPOLE_2D/define.h
index 80501a9..f1adbde 100755
--- a/PROBLEMS/MONOPOLE_2D/define.h
+++ b/PROBLEMS/MONOPOLE_2D/define.h
@@ -36,18 +36,18 @@
 //rmhd floors
 /************************************/
 
-#define B2RHORATIOMAXINIT 100
-#define B2UURATIOMAXINIT 100
+#define B2RHORATIOMAXINIT 100//100
+#define B2UURATIOMAXINIT 100//100
 
 #if defined(FORCEFREE)
 
+//#define ENFORCEENTROPY
 #define HYBRID_FORCEFREE
-#define HYBRID_FORCEFREE_SIGMACUT 50
+#define HYBRID_FORCEFREE_SIGMACUT 25
 #define FORCEFREE_SOLVE_PARALLEL
-//#define FORCEFREE_SOLVE_PARALLEL_OLD
+#//define FORCEFREE_SOLVE_PARALLEL_OLD
 
 //#define SKIPALLFLOORS // TODO seems critical for SOLVE_PARALLEL? 
-#define B2RHOFLOORFRAME DRIFTFRAME //ZAMOFRAME //DRIFTFRAME
 //#define CORRECT_POLARAXIS
 //#define NCCORRECTPOLAR 2
 
@@ -72,8 +72,8 @@
 #define B2UURATIOMIN 0.
 #define B2UURATIOMAX 2500.
 #define B2RHORATIOMIN 0.
-#define B2RHORATIOMAX 200.
-#define GAMMAMAXHD 50.
+#define B2RHORATIOMAX 500.
+#define GAMMAMAXHD 500.
 #endif
 
 /************************************/
@@ -145,7 +145,7 @@
 #define OUTVEL VEL4
 #define ALLSTEPSOUTPUT 0
 #define NSTEPSTOP 1.e10
-#define NOUTSTOP 102
+#define NOUTSTOP 110
 #define SILOOUTPUT 1
 #define PRIMOUTPUT 1
 //#define PRIMOUTPUTINMYCOORDS
diff --git a/finite.c b/finite.c
index cabdd6e..c21a66c 100644
--- a/finite.c
+++ b/finite.c
@@ -15,6 +15,7 @@ int
 reduce_order_check(ldouble *pm2,ldouble *pm1,ldouble *p0,ldouble *pp1,ldouble *pp2,int ix,int iy,int iz)
 {
   int reconstrpar=0;
+
 #ifdef REDUCEORDERTEMP
   ldouble t0,tp1,tm1,tmin;
   t0 =calc_PEQ_Tfromurho(p0[UU], p0[RHO]);
@@ -31,7 +32,22 @@ reduce_order_check(ldouble *pm2,ldouble *pm1,ldouble *p0,ldouble *pp1,ldouble *p
   if(xxBL[1]<REDUCEORDERRADIUS)
     reconstrpar=1;
 #endif
-  
+
+  // ANDREW TODO -- different behavior on boundary vs interior of FF region? 
+#if defined(FORCEFREE) && defined(REDUCEORDERFF)
+  int ffflag;
+  ffflag = get_cflag(FFINVFLAG,ix,iy,iz);
+  if(ffflag>0)
+    reconstrpar=1;
+#endif
+
+#ifdef REDUCEORDERATBH
+  ldouble xxBL[4];
+  get_xx_arb(ix,iy,iz,xxBL,BLCOORDS);
+  if(xxBL[1]<rhorizonBL)
+    reconstrpar=1;
+#endif
+
   return reconstrpar;
 }
 
@@ -93,20 +109,19 @@ avg2point(ldouble *um2,ldouble *um1,ldouble *u0,ldouble *up1,ldouble *up2,
 {
   ldouble r0[NV],rm1[NV],rp1[NV];
 
-  if(param!=0) //overrule the standard reconstruction
+  int int_order_local = INT_ORDER;
+  
+  if(param!=0) //reduce integration order
   {
-    int i;
-    if(param==1) //DONOR CELL
-    {
-      for(i=0;i<NV;i++)
-      {
-        ur[i]=u0[i];
-        ul[i]=u0[i];
-      }
-    }
+    //int_order_local = INT_ORDER-param;
+    int_order_local=0; // ANDREW TODO reduce directly to 0
   }  // if(param!=0)
+
+  // default to donor cell
+  if(int_order_local<0)
+    int_order_local = 0;
   
-  else if(INT_ORDER==0) //DONOR CELL
+  if(int_order_local==0) //DONOR CELL
   {
     int i;
     
@@ -117,7 +132,7 @@ avg2point(ldouble *um2,ldouble *um1,ldouble *u0,ldouble *up1,ldouble *up2,
     }
   }  // else if(INT_ORDER==0)
   
-  else if(INT_ORDER==1) //linear
+  else if(int_order_local==1) //linear
   {
     ldouble diffpar=theta;
     int i;
@@ -201,7 +216,7 @@ avg2point(ldouble *um2,ldouble *um1,ldouble *u0,ldouble *up1,ldouble *up2,
     } // for(i=0;i<NV;i++)
   }  // else if(INT_ORDER==1)
 
-  else if(INT_ORDER==2) //parabolic PPM
+  else if(int_order_local==2) //parabolic PPM
   {
     //The following is based on Colella & Woodward (J. Comp. Phys. 54, 174, 1984).
     //It uses five points: m2, m1, 0, p1, p2.
@@ -1876,7 +1891,7 @@ int
 set_grid(ldouble *mindx,ldouble *mindy, ldouble *mindz, ldouble *maxdtfac)
 {
   int i1,i2,ix,iy,iz;
-  ldouble mdx,mdy,mdz,dx,dy,dz,gloc[4][5],xx[4];
+  ldouble mdx,mdy,mdz,dx,dy,dz,gloc[4][5],xx[4],xxBL[4];
   mdx=mdy=mdz=-1;
   ldouble maxdt=-1;
 
@@ -1908,14 +1923,6 @@ set_grid(ldouble *mindx,ldouble *mindy, ldouble *mindz, ldouble *maxdtfac)
       if(i1>-NG) set_x(i1-1,2,.5*(get_xb(i1,2)+get_xb(i1-1,2)));
     }
 
-  //consistency check
-#if(MYCOORDS==BLCOORDS)
-  if(get_x(-1,0)<=rhorizonBL)
-    {
-      printf("ix %d > %f\n",-1,get_x(-1,0));
-      my_err("-1 cell inside horizon\n");
-    }
-#endif
 
   // what is the minimal cell size
   for(ix=ix1;ix<ix2;ix++)
@@ -2467,23 +2474,14 @@ int
 get_xx_arb(int ix,int iy,int iz,ldouble *xx,int COORDSOUT)
 {
   
-#ifdef PRECOMPUTE_MY2OUT // use precomputed coordinates if COORDS == OUTCOORDS
-  if(COORDSOUT == OUTCOORDS)
-  {
+#if defined(PRECOMPUTE_MY2OUT) && (COORDSOUT==OUTCOORDS) // use precomputed coordinates if COORDS == OUTCOORDS
     get_xxout(ix, iy, iz, xx); // time will be nonsense! seems ok everywhere this is used
-  }
-  else
-  {
-    ldouble xx0[4];
-    get_xx(ix,iy,iz,xx0);
-    coco_N(xx0,xx,MYCOORDS,COORDSOUT);
-  }
 #else
     ldouble xx0[4];
     get_xx(ix,iy,iz,xx0);
     coco_N(xx0,xx,MYCOORDS,COORDSOUT);
-#endif  
-  return 0;
+#endif
+    return 0;
 }
 
 
@@ -5187,8 +5185,7 @@ cell_fixup(int type)
 	  (get_cflag(RADIMPFIXUPFLAG,ix,iy,iz)!=0 && type==FIXUP_RADIMP)) && is_cell_active(ix,iy,iz)
 	)
 	{
-	  //ANDREW - I think maybe this should be here, to set the flag back to its default? 
-	  //this should not be here, should it?
+
 	  /*
 	  if(type==FIXUP_U2PMHD) set_cflag(HDFIXUPFLAG,ix,iy,iz,0); //try only once
 	  if(type==FIXUP_U2PRAD) set_cflag(RADFIXUPFLAG,ix,iy,iz,0); //try only once
diff --git a/frames.c b/frames.c
index f1810fb..ee65542 100644
--- a/frames.c
+++ b/frames.c
@@ -95,8 +95,9 @@ trans_pmhd_coco(ldouble *ppin, ldouble *ppout, int CO1,int CO2, ldouble *xxvec,
     }
 
 #ifdef FORCEFREE
-    pp2[UUFF]=pp1[UUFF];
-      
+
+
+    // perpindicular
     ldouble uconff[4];
     uconff[0]=0.;
     uconff[1]=pp1[VXFF];
@@ -108,6 +109,13 @@ trans_pmhd_coco(ldouble *ppin, ldouble *ppout, int CO1,int CO2, ldouble *xxvec,
     pp2[VXFF]=uconff[1];
     pp2[VYFF]=uconff[2];
     pp2[VZFF]=uconff[3];
+
+    //parallel
+    //ANDREW TODO this is NOT right
+    //but it doesn't matter for now as long as FF prims are not used in evolution
+    //FF conserveds are all derived from normal prims
+    pp2[UUFF]=pp1[UUFF];
+    
 #endif
     
 #endif
diff --git a/ko.c b/ko.c
index e3e09ab..6c3f044 100644
--- a/ko.c
+++ b/ko.c
@@ -96,7 +96,8 @@ main(int argc, char **argv)
 
   //precalculates metric, Christoffels, etc.
   calc_metric();
-  
+  calc_cells_under_horiz();
+ 
   //save coordinate file
   #ifdef COORDOUTPUT
   fprint_coordfile(folder,"coord");
@@ -112,12 +113,6 @@ main(int argc, char **argv)
 #endif 
 
   //print scalings GU->CGS
-  #ifdef FORCEFREE
-  if(PROCID==0) printf("NVHD %d NVMHD %d NV %d\n",NVHD,NVMHD,NV);
-  if(PROCID==0) printf("%d %d %d %d %d %d %d %d %d %d %d %d %d\n",
-		       RHO,UU,VX,VY,VZ,ENTR,B1,B2,B3,UUFF,VXFF,VYFF,VZFF);
-  #endif
-
   if(PROCID==0) print_scalings();
 
   //**************
@@ -177,6 +172,11 @@ main(int argc, char **argv)
         printf("Sending initial data... ");
 	fflush(stdout);
       }
+
+#ifdef FORCEFREE
+      fill_ffprims(); // make force-free primitives consistent
+#endif
+
       //exchange initial state
       mpi_exchangedata();  
       calc_avgs_throughout();
@@ -284,6 +284,7 @@ main(int argc, char **argv)
   int iix,iiy,iiz;
   struct geometry geom;
   ldouble PERTURBATION;
+
   
   // reset Te so that ue/ugas is constant
 #ifdef RESETELECTRONTEMPERATURETOUEUGASRATIO
diff --git a/ko.h b/ko.h
index d114005..95f1f73 100644
--- a/ko.h
+++ b/ko.h
@@ -177,6 +177,7 @@ ldouble min_dx,min_dy,min_dz;
 
 //precalculated metric parameters
 ldouble rhorizonBL,rISCOBL,rmboundBL,rphotonBL,etaNT;
+int cells_under_horizon;
 
 // jet coordinate specific
 #if (MYCOORDS==JETCOORDS)
@@ -676,6 +677,7 @@ int my_finger(ldouble t);
 ///////////////////////////////////////////////////////////////
 
 void initialize_constants();
+int calc_cells_under_horiz();
 int print_scalings();
 int set_initial_profile();
 void am_i_sane();
@@ -1003,15 +1005,16 @@ int calc_primitives(int ix,int iy,int iz,int type,int setflags);
 int u2p(ldouble *uu0, ldouble *pp,void *ggg,int corrected[3],int fixups[2],int type);
 int check_floors_mhd(ldouble *pp, int whichvel,void *ggg);
 
-int u2p_solver(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose);
+int u2p_solver_mhd(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose);
 int u2p_solver_nonrel(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose);
 int u2p_solver_Wp(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose);
 int u2p_solver_W(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose);
 int u2p_solver_Bonly(ldouble *uu, ldouble *pp, void *ggg);
-int u2p_solver_ff(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose);
+int u2p_solver_ff(ldouble *uu, ldouble *pp, void *ggg,int verbose);
 int fill_ffprims();
 int fill_ffprims_cell(ldouble *pp, void *gg);
-
+ldouble get_driftvel_velr(ldouble *pp,ldouble *velff,void* ggg);
+  
 int count_entropy(int *n, int *n2);
 int copy_entropycount();
 int update_entropy();
@@ -1334,6 +1337,8 @@ ldouble pick_ViscousHeating(int ix,int iy,int iz);
 int test_gammagas();
 int test_calcgamma();
 
+int calc_faraday(ldouble F[4][4],int ix,int iy,int iz,int when);
+int calc_current(void* ggg,ldouble jcon[4],int *derdir);
 ///////////////////////////////////////////////////////////////
 // nonthermal.c ///////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////
diff --git a/misc.c b/misc.c
index f78b6fc..5f8973a 100644
--- a/misc.c
+++ b/misc.c
@@ -190,10 +190,55 @@ void initialize_constants()
   */
   
   #endif
+
   
   return;
 }
 
+//**********************************************************************
+/*! \fn int calc_cells_under_horiz()
+ \brief calculate the number of cells beneath the horizon for bhdisk problems
+*/
+//**********************************************************************
+
+int
+calc_cells_under_horiz()
+{
+    //consistency check -- how many cell centers are under horizon? 
+#ifdef BHDISK_PROBLEMTYPE
+  ldouble xx[4],xxBL[4];
+  int ix;
+  if(TOI==0 && TOJ==0 && TOK==0) // only do on 0th tile
+  {
+    cells_under_horizon=0;
+    for(ix=0;ix<NX;ix++)
+    {
+      get_xx(ix,0,0,xx); //BH problem types should have r coord independent of theta,phi
+      coco_N(xx,xxBL,MYCOORDS,BLCOORDS);
+      if(xxBL[1]>rhorizonBL)
+	break;
+      else
+	cells_under_horizon+=1;
+    }
+    
+    if(cells_under_horizon<3)
+    {
+      printf("There are only %d cells under horizon at rh=%.2f! increase to at least 4\n",
+	     cells_under_horizon,rhorizonBL);
+      exit(-1);
+    }
+    if(cells_under_horizon==NX)
+    {
+      printf("All cells on inner tile are under horizon! \n");
+      exit(-1);
+    }
+
+    printf("There are %d cells under the horizon at rh=%.2f\n",
+	   cells_under_horizon,rhorizonBL);
+  }
+#endif
+  return 0;
+}
 
 //**********************************************************************
 /*! \fn int print_scalings()
@@ -358,7 +403,6 @@ if (NTZ % 2 != 0)
 #endif
 
 #ifdef FORCEFREE
-  printf("WARNING! FORCEFREE is untested and does not account for current sheets! use only with monopole!\n");
   #ifdef NONRELMHD
   printf("FORCEFREE does not work with NONRELMHD!\n");
   exit(-1);
@@ -539,6 +583,8 @@ if (NTZ % 2 != 0)
   if(PROCID==0) printf("RESTARTFROMMHD\n");
   if(PROCID==0) printf("urad/uu: %e ue/uu: %e\n", INITURADFRAC, INITUEFRAC);
 #endif
+
+
   return;
 }
 
@@ -1476,7 +1522,7 @@ print_primitives(ldouble *p)
   printf("B^2 = %.15e\n",p[B2]);
   printf("B^3 = %.15e\n",p[B3]);
 #ifdef FORCEFREE
-  printf("uu (ff) = %.15e\n", p[UUFF]);
+  printf("upar (ff) = %.15e\n", p[UUFF]);
   printf("u^1 (ff) = %.15e\n",p[VXFF]);
   printf("u^2 (ff) = %.15e\n",p[VYFF]);
   printf("u^3 (ff) = %.15e\n",p[VZFF]);
@@ -1525,6 +1571,7 @@ print_conserved(ldouble *u)
   printf("B^2 = %.15e\n",u[B2]);
   printf("B^3 = %.15e\n",u[B3]);
 #ifdef FORCEFREE
+  printf("mu b^0 (ff) = %.15e\n",u[UUFF]);
   printf("T^t_1 (ff) = %.15e\n",u[VXFF]);
   printf("T^t_2 (ff) = %.15e\n",u[VYFF]);
   printf("T^t_3 (ff) = %.15e\n",u[VZFF]);
diff --git a/mnemonics.h b/mnemonics.h
index 90af3df..1f2442e 100644
--- a/mnemonics.h
+++ b/mnemonics.h
@@ -109,7 +109,7 @@
 #define ZBCHI 6
 
 //cell flags
-#define NFLAGS 7
+#define NFLAGS 8
 
 //boolean
 #define ENTROPYFLAG 0
@@ -119,6 +119,7 @@
 #define ENTROPYFLAG2 4
 #define ENTROPYFLAG3 5
 #define RADSOURCETYPEFLAG 6
+#define FFINVFLAG 7
 
 //values for RADFIXUP
 #define RADSOURCETYPEEXPLICIT 1
diff --git a/p2u.c b/p2u.c
index 3fd7bb5..27e5350 100644
--- a/p2u.c
+++ b/p2u.c
@@ -116,7 +116,7 @@ p2u_mhd(ldouble *pp, ldouble *uu, void *ggg)
   for(ib=0;ib<NRELBIN;ib++)
   {
     uu[NEREL(ib)]=gdetu*pp[NEREL(ib)]*ut;
-  }    
+  }    x
 #endif
 
   //************************************
@@ -127,32 +127,22 @@ p2u_mhd(ldouble *pp, ldouble *uu, void *ggg)
   uu[B2]=gdetu*pp[B2];
   uu[B3]=gdetu*pp[B3];
 
+  // ANDREW TODO I changed to use regular prims instead of FF prims
 #ifdef FORCEFREE
-  // FF conserved parallel enthalpy
-  ldouble ugas_ff = pp[UUFF];
-  ldouble w_s = 1 + gamma*ugas_ff/rho; // specific enthalpy
-  uu[UUFF] = gdetu * w_s * bcon[0];
-  
-  // TODO -- do we need to recompute these things?
-  // TODO -- T^\mu\nu EM should be invariant to any parallel velocity
-  ldouble vcon_ff[4],ucon_ff[4],ucov_ff[4];
-  ldouble bcon_ff[4],bcov_ff[4],bsq_ff;
-  
-  // FF stress-energy components
-  // ANDREW TODO -- necessary to recompute ucon, etc?
-  vcon_ff[0]=0.;
-  vcon_ff[1]=pp[VXFF];
-  vcon_ff[2]=pp[VYFF];
-  vcon_ff[3]=pp[VZFF];
 
-  conv_vels_both(vcon_ff,ucon_ff,ucov_ff,VELPRIM,VEL4,gg,GG);
-  calc_bcon_bcov_bsq_from_4vel(pp, ucon_ff, ucov_ff, geom, bcon_ff, bcov_ff, &bsq_ff);
+#ifdef FORCEFREE_PARALLEL_COLD // neglect pressure
+  uu[UUFF] = gdetu * bcon[0];
+#else // adiabatic, with pressure
 
-  //ldouble Tttt_ff = bsq*ucon[0]*ucov[0] - bcon[0]*bcov[0] + 0.5*bsq + rho*ucon[0];  
-  ldouble Ttr_ff  = bsq_ff*ucon_ff[0]*ucov_ff[1] - bcon_ff[0]*bcov_ff[1];
-  ldouble Ttth_ff = bsq_ff*ucon_ff[0]*ucov_ff[2] - bcon_ff[0]*bcov_ff[2];
-  ldouble Ttph_ff = bsq_ff*ucon_ff[0]*ucov_ff[3] - bcon_ff[0]*bcov_ff[3];
+  ldouble w_s = 1 + gamma*pp[UU]/pp[RHO]; // specific enthalpy
+  uu[UUFF] = gdetu * w_s * bcon[0];
+#endif
+  
+  ldouble Ttr_ff  = bsq*ucon[0]*ucov[1] - bcon[0]*bcov[1];
+  ldouble Ttth_ff = bsq*ucon[0]*ucov[2] - bcon[0]*bcov[2];
+  ldouble Ttph_ff = bsq*ucon[0]*ucov[3] - bcon[0]*bcov[3];
 
+  
   uu[VXFF] = gdetu*Ttr_ff;
   uu[VYFF] = gdetu*Ttth_ff;
   uu[VZFF] = gdetu*Ttph_ff;
diff --git a/physics.c b/physics.c
index f60a8a0..77efade 100644
--- a/physics.c
+++ b/physics.c
@@ -1266,7 +1266,7 @@ int f_flux_prime(ldouble *pp, int idim, int ix, int iy, int iz,ldouble *ff,int l
 #endif
 
   //magnetic fluxes
-  // ANDREW TODO -- need to modify for FORCEFREE? 
+
 #ifdef MAGNFIELD
   ff[B1]=gdetu*(bcon[1]*ucon[idim+1] - bcon[idim+1]*ucon[1]);
   ff[B2]=gdetu*(bcon[2]*ucon[idim+1] - bcon[idim+1]*ucon[2]);
@@ -1286,11 +1286,16 @@ int f_flux_prime(ldouble *pp, int idim, int ix, int iy, int iz,ldouble *ff,int l
 
 #ifdef FORCEFREE
 
-  // FF conserved parallel enthalpy
-  ldouble ugas_ff = pp[UUFF];
-  ldouble w_s = 1. + gamma*ugas_ff/rho; // specific enthalpy  
-  ff[UUFF]= gdetu*w_s*bcon[idim+1];
+  // ANDREW TODO modified to usue regular prims
 
+#ifdef FORCEFREE_PARALLEL_COLD // neglect pressure
+  ff[UUFF]= gdetu*bcon[idim+1];
+#else // adiabatic, with pressure
+  
+  ldouble w_s = 1. + gamma*ugas/rho; // specific enthalpy
+  ff[UUFF]= gdetu*w_s*bcon[idim+1];
+#endif
+  
   // FF magnetic fluxes
   ldouble T_ff[4][4];
   calc_Tij_ff(pp,&geom,T_ff);
@@ -1386,12 +1391,6 @@ calc_Tij(ldouble *pp, void* ggg, ldouble T[][4])
    for(i=1;i<4;i++)
      T[0][i]=T[i][0]=(T[0][0] + ptot) *ucon[i]*ucon[0] + ptot*GG[i][0] - bcon[i]*bcon[0];
 
-   /*
-#elif defined(FORCEFREE) // ANDREW get rid of this part
-  for(i=0;i<4;i++)
-    for(j=0;j<4;j++)
-      T[i][j]=bsq*ucon[i]*ucon[j] + 0.5*bsq*GG[i][j] - bcon[i]*bcon[j];
-   */
 #else //normal GRMHD
   for(i=0;i<4;i++)
     for(j=0;j<4;j++)
@@ -1418,14 +1417,14 @@ calc_Tij_ff(ldouble *pp, void* ggg, ldouble T[][4])
   ldouble utcon[4],ucon[4],ucov[4];  
   ldouble bcon[4],bcov[4],bsq=0.;
 
-  // ANDREW TODO : which 4-velocity to use here generally? 
+  // ANDREW TODO changed to regular prims
   // as long as perp components are the same, it shouldn't matter,
   // since parallel components cancel out
 
   utcon[0]=0.;
-  utcon[1]=pp[VXFF]; // TODO or pp[VX] ? 
-  utcon[2]=pp[VYFF]; // TODO or pp[VY] ? 
-  utcon[3]=pp[VZFF]; // TODO or pp[VZ] ? 
+  utcon[1]=pp[VX]; 
+  utcon[2]=pp[VY]; 
+  utcon[3]=pp[VZ]; 
 
   conv_vels_both(utcon,ucon,ucov,VELPRIM,VEL4,gg,GG);
 
@@ -1780,7 +1779,6 @@ calc_ufromS3rho(ldouble S3,ldouble rho,int type,int ix,int iy,int iz)
  return uint;
 }
 
-
 //**********************************************************************
 //updates entropy in the specified cell (p[ENTR]) basing on new primitives
 //or stays with the old one if entropy u2p solver was involved
@@ -3104,10 +3102,231 @@ pick_ViscousHeating(int ix,int iy,int iz)
 }
 
 
+//***********************************************************************************
+//calculates current 4-vector
+// if when=-1, use ptm1
+//derdir[] determines the type of derivative in each dimension (left,right,centered)
+//************************************************************************************
+int
+calc_faraday(ldouble Ft[4][4],int ix,int iy,int iz,int when)
+{
+  struct geometry geom;
+  fill_geometry(ix,iy,iz,&geom);
+
+  int iv;
+  ldouble pp[NV];
+  for(iv=0;iv<NV;iv++)
+  {
+    if(when==-1) pp[iv]=get_u(ptm1,iv,ix,iy,iz);
+    else pp[iv]=get_u(p,iv,ix,iy,iz);
+  }
+	      
+  //start with computing F^\mu\nu
+
+  ldouble ucon[4],ucov[4],bcon[4],bcov[4],bsq;
+  calc_ucon_ucov_from_prims(pp, &geom, ucon, ucov);
+  calc_bcon_bcov_bsq_from_4vel(pp, ucon, ucov, &geom, bcon, bcov, &bsq);
+
+  ldouble Fstarcov[4][4];
+  int i,j,k,l;
+  for(i=0;i<4;i++)
+  {
+    for(j=0;j<4;j++)
+    {
+      Fstarcov[i][j] = bcov[i]*ucov[j] - bcov[j]*ucov[i];
+    }
+  }
+      
+  for(i=0;i<4;i++)
+  {
+    for(j=0;j<4;j++)
+    {
+      Ft[i][j]=0;
+      for(k=0;k<4;k++)
+      {
+	for(l=0;l<4;l++)
+	{
+	  int sign = epsilon_4d(i,j,k,l);
+	  //Ft[i][j] += sign*bcov[k]*ucov[l];
+	  Ft[i][j] += 0.5*sign*Fstarcov[k][l];
+	}
+      }
+    Ft[i][j] /= geom.gdet;
+    }
+  }
+
+  return 0;
+}
+
+
+int
+calc_current(void* ggg,ldouble jcon[4],int *derdir)
+{
+  int i,j,k,l,iv;
+
+  struct geometry *geomin
+    = (struct geometry *) ggg;
+
+  int ix,iy,iz;
+  ix=geomin->ix;
+  iy=geomin->iy;
+  iz=geomin->iz;
+
+  //lets get geometry again, just in case geomin is not in internal coordinates
+  struct geometry geom;
+  fill_geometry(ix,iy,iz,&geom);
+
+  ldouble dF[4];
+
+  //instead:
+  for(i=0;i<4;i++)
+  {
+    //force d/dt = 0 for now // ANDREW TODO
+    dF[i] = 0.;
+  }
+
+  ldouble Fm1[4][4],Fp1[4][4],Fc[4][4];
+  ldouble xxvec[4],xxvecm1[4],xxvecp1[4];
+  int idim;
+  calc_faraday(Fc,ix,iy,iz,0);
+  get_xx(ix,iy,iz,xxvec);
+
+  #ifdef CURRENTTIMEDERIV
+  if(!doingpostproc)
+  {
+    calc_faraday(Fm1,ix,iy,iz,-1); //Faraday tensor at previous timestep
+    for(i=0;i<4;i++)
+      dF[i] += (Fc[i][0] - Fm1[i][0])/dt; 
+  }
+  #endif
+
+  //derivatives
+  for(idim=1;idim<4;idim++)
+  {
+   
+
+      if(idim==1)
+      {
+	  get_xx(ix-1,iy,iz,xxvecm1);
+	  get_xx(ix+1,iy,iz,xxvecp1);
+	  calc_faraday(Fm1,ix-1,iy,iz,0);
+	  calc_faraday(Fp1,ix+1,iy,iz,0);
+      }
+
+      else if(idim==2)
+      {
+	  get_xx(ix,iy-1,iz,xxvecm1);
+	  get_xx(ix,iy+1,iz,xxvecp1);	  
+	  calc_faraday(Fm1,ix,iy-1,iz,0);
+	  calc_faraday(Fp1,ix,iy+1,iz,0);	    
+      }
+
+      else if(idim==3)
+      {
+	 get_xx(ix,iy,iz-1,xxvecm1);
+	 get_xx(ix,iy,iz+1,xxvecp1);
+	 calc_faraday(Fm1,ix,iy,iz-1,0);
+	 calc_faraday(Fp1,ix,iy,iz+1,0);	    
+      }
+
+     // if only one cedll
+     if(TNY==1 && idim==2) continue;
+     if(TNZ==1 && idim==3) continue;
+ 
+     ldouble dl,dr,dc;
+     
+     for(i=0;i<4;i++)
+     {
+
+	 dc=(Fp1[i][idim] -Fm1[i][idim]) / (xxvecp1[idim] - xxvecm1[idim]);
+	 dr=(Fp1[i][idim] -Fc[i][idim]) / (xxvecp1[idim] - xxvec[idim]);
+	 dl=(Fc[i][idim]  -Fm1[i][idim]) / (xxvec[idim] - xxvecm1[idim]);
+	 //TODO : whether MPI handles this properly
+	 //to avoid corners
+	 if(TNX==1 && idim==1)
+	 {
+	   dF[i] += 0;
+	 }
+	 else if(TNY==1 && idim==2)
+	 {
+	   dF[i] += 0;
+	 }
+	 else if(TNZ==1 && idim==3)
+	 {
+	   dF[i] += 0;
+	 }
+	 else if((ix<0 && iy==0 && iz==0 && idim!=1) ||
+	    (iy<0 && ix==0 && iz==0 && idim!=2) ||
+	    (iz<0 && ix==0 && iy==0 && idim!=3))
+	 {
+	     dF[i] += dr;
+	 }
+	 else if((ix<0 && iy==NY-1 && iz==NZ-1 && idim!=1) ||
+	    (iy<0 && ix==NX-1 && iz==NZ-1 && idim!=2) ||
+	    (iz<0 && ix==NX-1 && iy==NY-1 && idim!=3))
+	 {
+	     dF[i] += dl;
+	 }
+	 else if((ix>=NX && iy==0 && iz==0 && idim!=1) ||
+	    (iy>=NY && ix==0 && iz==0 && idim!=2) ||
+	    (iz>=NZ && ix==0 && iy==0 && idim!=3))
+	 {
+	     dF[i] += dr;
+	 }
+	 else if((ix>=NX && iy==NY-1 && iz==NZ-1 && idim!=1) ||
+	    (iy>=NY && ix==NX-1 && iz==NZ-1 && idim!=2) ||
+	    (iz>=NZ && ix==NX-1 && iy==NY-1 && idim!=3))
+	 {
+	     dF[i] += dl;
+	 }
+	 else
+	 {
+	     //choice of 1st order derivative
+	     if(derdir[idim-1]==0)
+	       {
+		 dF[i] += dc;
+	       }
+	     if(derdir[idim-1]==1)
+	       {
+		 dF[i] += dl;
+	       }
+	     if(derdir[idim-1]==2)
+	       {
+		 dF[i] += dr;
+	       }
+	 }
+     }  // loop over i     
+  } // loop over idim
+
+  //covariant derivative tensor du^i;j - only for expansion
+  ldouble dcu2[4][4];
+  ldouble Krsum[4];
+
+  for(i=0;i<4;i++)
+  {
+    Krsum[i]=0;
+    for(j=0;j<4;j++)
+    {
+      for(k=0;k<4;k++)
+      {
+	Krsum[i] += get_gKr(j,j,k,ix,iy,iz)*Fc[i][k];
+      }
+    }
+  }
+
+  //total current
+  for(i=0;i<4;i++)
+  {
+    jcon[i] = dF[i]+Krsum[i];
+  }
+
+  return 0;
+}
+
 
 //*********************************************************************************
 //tests
-//*********************************************************************************
+//******************************************************************x***************
 
 
 /*! int test_gammagas()
diff --git a/problem.c b/problem.c
index 1344fa0..fb94a97 100644
--- a/problem.c
+++ b/problem.c
@@ -753,7 +753,7 @@ solve_the_problem(ldouble tstart, char* folder)
             #endif
 	    
 	    // Calculate primitves
-	    // ANDREW is this excessive? Should be consistent after implicit!
+	    // ANDREW TODO is this excessive? Should be consistent after implicit!
 	    calc_u2p(0,1); //do not calculate visc. heating, do count entropy inversions
 	    
 	    // Heat species at end
@@ -783,10 +783,11 @@ solve_the_problem(ldouble tstart, char* folder)
       
       //for outputs - use what came out of 2nd implicit: 
       //ANDREW why?? try current prims instead
+      // ANDREW: try saving out current prims instead of postimplicit
       
       // Set uforget = p and p = ppostimplicit over domain
       copy_u(1.,p,uforget); //backup current primitives
-      //copy_u(1.,ppostimplicit,p); // ANDREW: try saving out current prims
+      //copy_u(1.,ppostimplicit,p); 
    
       //counting faiures and average parameters of the implicit solver
       int nfailures[3],nfailuresloc[3]={global_int_slot[GLOBALINTSLOT_NTOTALRADIMPFAILURES],
@@ -880,7 +881,7 @@ solve_the_problem(ldouble tstart, char* folder)
       
 	  //zero out avg values over domain
 
-    long long Navg = (long long) SX*SY*SZ*(NV+NAVGVARS); // RN: Apr 2, 2019
+          long long Navg = (long long) SX*SY*SZ*(NV+NAVGVARS); // RN: Apr 2, 2019
 	  copy_u_core(0., pavg, pavg, Navg);
 	  avgtime=0.;
 	  
@@ -912,6 +913,9 @@ solve_the_problem(ldouble tstart, char* folder)
 	  // Set boundary conditions on conserveds in the ghost cells
 	  set_bc(t,0);
 
+	  #ifdef FORCEFREE
+	  fill_ffprims();
+	  #endif
 	  //print restart file
 	  fprint_restartfile(t,folder);
 
@@ -952,7 +956,6 @@ solve_the_problem(ldouble tstart, char* folder)
 
 
 	  nfout1++;
-
 	  
 	  #ifdef DTOUT_LOG
           dtout = pow(10, DTOUT1_LOG_INIT + DTOUT_LOG*(nfout1-1));
diff --git a/problem.h b/problem.h
index 7654bea..7087025 100644
--- a/problem.h
+++ b/problem.h
@@ -152,8 +152,8 @@
 //if you are using a problem older than 100 with a different kappa defined in kappa.c
 //make sure your kappa.c ends with (kappa=(.....)) NOT (return kappa)!!
 
-#define PROBLEM 132
-//#define PROBLEM 144
+//#define PROBLEM 132
+#define PROBLEM 144
 
 #if(PROBLEM==144)
 
diff --git a/rad.c b/rad.c
index 793bd86..989c325 100644
--- a/rad.c
+++ b/rad.c
@@ -1394,10 +1394,10 @@ implicit_apply_constraints(ldouble *pp, ldouble *uu, ldouble *uu0, void* ggg, in
    
     //and invert back to primitives
     int rettemp=0;
-    rettemp=u2p_solver(uu,pp,geom,U2P_HOT,0); 
+    rettemp=u2p_solver_mhd(uu,pp,geom,U2P_HOT,0); 
     #ifdef ALLOWFORENTRINF4DPRIM
     if(rettemp<0)
-      rettemp=u2p_solver(uu,pp,geom,U2P_ENTROPY,0);
+      rettemp=u2p_solver_mhd(uu,pp,geom,U2P_ENTROPY,0);
     #endif
     //u2pret=rettemp;
     if(rettemp<0) u2pret=-2; 
@@ -5538,7 +5538,7 @@ test_jon_solve_implicit_lab()
       int corr[2],fixup[2],u2pret,radcor;
      
       //printf("inverting...\n");
-      u2pret=u2p_solver(uu,pp,&geom,U2P_HOT,0); //hd
+      u2pret=u2p_solver_mhd(uu,pp,&geom,U2P_HOT,0); //hd
       if(u2pret<0) printf("u2pret mhd: (%d)\n",u2pret);
       u2p_rad(uu,pp,&geom,&radcor); //rad
       if(radcor!=0) printf("u2pcor rad: (%d)\n",radcor);
@@ -5559,7 +5559,7 @@ test_jon_solve_implicit_lab()
       if(s2/s1 < 0.9 | u2pret<0.)
 	{ 
 	  printf("\n PROBLEM DETECTED IN ENTROPY OR U2P_HOT DID NOT SUCCEED!\n");
-	  u2pret=u2p_solver(uu,pp,&geom,U2P_ENTROPY,0); //hd
+	  u2pret=u2p_solver_mhd(uu,pp,&geom,U2P_ENTROPY,0); //hd
 	  if(u2pret<0) printf("u2pret mhd: (%d)\n",u2pret);
 	  printf("\n..........................\nafter u2p_ENTROPY:\n\n");
 	  print_Nvector(pp,NV);
diff --git a/silo.c b/silo.c
index 73be334..987a861 100644
--- a/silo.c
+++ b/silo.c
@@ -68,10 +68,8 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
   ldouble *rho = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *entr = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *uint = (ldouble*)malloc(nx*ny*nz*sizeof(double));
-    ldouble *uintff = (ldouble*)malloc(nx*ny*nz*sizeof(double));
+
   ldouble *temp = (ldouble*)malloc(nx*ny*nz*sizeof(double));
-  ldouble *trad = (ldouble*)malloc(nx*ny*nz*sizeof(double));
-  ldouble *tradlte = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *Omega = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *muBe = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *Be = (ldouble*)malloc(nx*ny*nz*sizeof(double));
@@ -86,13 +84,10 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
   ldouble *u3 = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *lorentz = (ldouble*)malloc(nx*ny*nz*sizeof(double));
 
-  ldouble *u0_perp = (ldouble*)malloc(nx*ny*nz*sizeof(double));
-  ldouble *u1_perp = (ldouble*)malloc(nx*ny*nz*sizeof(double));
-  ldouble *u2_perp = (ldouble*)malloc(nx*ny*nz*sizeof(double));
-  ldouble *u3_perp = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *lorentz_perp = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *lorentz_par = (ldouble*)malloc(nx*ny*nz*sizeof(double));
-  
+  int *ffinv = (int*)malloc(nx*ny*nz*sizeof(int));
+    
   ldouble *vx = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *vy = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *vz = (ldouble*)malloc(nx*ny*nz*sizeof(double));
@@ -118,6 +113,7 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
   ldouble *sigmaw = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *betainv = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *divB = (ldouble*)malloc(nx*ny*nz*sizeof(double));
+  ldouble *jdens = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *OmegaF = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *Qtheta = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *Qphi = (ldouble*)malloc(nx*ny*nz*sizeof(double));
@@ -183,6 +179,9 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
   #endif  
 
   #ifdef RADIATION
+  ldouble *trad = (ldouble*)malloc(nx*ny*nz*sizeof(double));
+  ldouble *tradlte = (ldouble*)malloc(nx*ny*nz*sizeof(double));
+
   ldouble *taucoupling = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *tausca = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   ldouble *tauabs = (ldouble*)malloc(nx*ny*nz*sizeof(double));
@@ -445,11 +444,13 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
       
       if(doingavg==0) //using snapshot data
       {
+	
 	rho[zonalindex]  = pp[RHO];
 	entr[zonalindex] = pp[ENTR];
 	uint[zonalindex] = pp[UU];
 	pregas = (gamma-1.)*pp[UU];
-	
+
+	vel[0] = 0.;
 	vel[1]=pp[VX];
 	vel[2]=pp[VY];
 	vel[3]=pp[VZ];
@@ -479,10 +480,6 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
 	sigma[zonalindex] = bsq[zonalindex]/(rho[zonalindex]);
 	sigmaw[zonalindex] = bsq[zonalindex]/(rho[zonalindex] + gamma*uint[zonalindex]);
 
-	// field line angular velocity
-	// TODO why doesn't this look right? 
-	//OmegaF[zonalindex] = vel[3]/vel[0] - (vel[1]/vel[0])*(pp[B3]/pp[B1]);
-	OmegaF[zonalindex] = (Omega[zonalindex] - (vel[2]*B3comp[zonalindex])/(vel[0]*B2comp[zonalindex]));
         #endif		  
 
 	// stress energy
@@ -564,6 +561,7 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
       }
       else //using averaged data
       {
+
 	 rho[zonalindex] = get_uavg(pavg,RHO,ix,iy,iz);
 	 entr[zonalindex] = get_uavg(pavg,ENTR,ix,iy,iz);
 	 uint[zonalindex] = get_uavg(pavg,UU,ix,iy,iz);
@@ -607,9 +605,7 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
 	 sigma[zonalindex]=bsq[zonalindex]/(rho[zonalindex] + uint[zonalindex] + pregas);
 	 sigmaw[zonalindex]=bsq[zonalindex]/(rho[zonalindex] + gamma*uint[zonalindex]);
 
-	 // field line angular velocity
-	 //OmegaF[zonalindex] = vel[3]/vel[0] - (vel[1]/vel[0])*(pp[B3]/pp[B1]);
-	 OmegaF[zonalindex] = (Omega[zonalindex] - (vel[2]*B3comp[zonalindex])/(vel[0]*B2comp[zonalindex]));
+     
          #endif		  
 	 
 	 // stress energy tensor
@@ -699,6 +695,9 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
 	 deltae[zonalindex]=-1.;
 	 dtauarr[zonalindex]=-1.;
 #endif //PRINTVISCHEATINGTOSILO	 
+
+
+	 
       } //doingavg
 
       
@@ -741,33 +740,47 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
 
       // put quantities into arrays
       expansion[zonalindex]=exploc;
-
       lorentz[zonalindex]=fabs(vel[0])/sqrt(fabs(geomout.GG[0][0]));
       u0[zonalindex]=fabs(vel[0]);
       u1[zonalindex]=vel[1];
       u2[zonalindex]=vel[2];
       u3[zonalindex]=vel[3];
 
-#ifdef FORCEFREE // must have VELPRIM=VELR
+#ifdef MAGNFIELD
+      //OmegaF[zonalindex] = Omega[zonalindex] - (u1[zonalindex]/u0[zonalindex])*(B3comp[zonalindex]/B1comp[zonalindex]);
+
+      OmegaF[zonalindex]=Omega[zonalindex];
+      if(fabs(u2[zonalindex]*B3comp[zonalindex])>1.e-8)
+	OmegaF[zonalindex] -= (u2[zonalindex]*B3comp[zonalindex])/(u0[zonalindex]*B2comp[zonalindex]);
+
       ldouble velff[4];
-      velff[0] = 0.;
-      velff[1]=pp[VXFF];
-      velff[2]=pp[VYFF];
-      velff[3]=pp[VZFF];
+      ldouble vpar;
+
+      vpar = get_driftvel_velr(pp,velff,&geomout);
       conv_vels(velff,velff,VELR,VEL4,geomout.gg,geomout.GG);
+      
       lorentz_perp[zonalindex]=fabs(velff[0])/sqrt(fabs(geomout.GG[0][0]));
       lorentz_par[zonalindex]=1./sqrt(1+pow(lorentz[zonalindex],-2)-pow(lorentz_perp[zonalindex],-2));
 
-      u0_perp[zonalindex]=velff[0];
-      u1_perp[zonalindex]=velff[1];
-      u2_perp[zonalindex]=velff[2];
-      u3_perp[zonalindex]=velff[3];
+#endif
+
+
+#ifdef FORCEFREE // must have VELPRIM=VELR
+
+      ffinv[zonalindex]=get_cflag(FFINVFLAG,ix,iy,iz);
+      
+      int derdir2[3] = {2,2,2};
+      ldouble jcon[4],jcov[4];
+      calc_current(&geom,jcon,derdir2);
+      indices_21(jcon,jcov,geom.gg); // not in lab frame
+
+      //ldouble FF[4][4];
+      //calc_faraday(FF,geom.ix,geom.iy,geom.iz,0);
+      jdens[zonalindex] = sqrt(fabs(dot(jcon,jcov)));
 
-      uintff[zonalindex]=pp[UUFF];
 #else
-      //TODO actually get || and perp velocities here
-      lorentz_perp[zonalindex]=-1;
-      lorentz_par[zonalindex]=-1;
+      jdens[zonalindex]=0.;      
+      ffinv[zonalindex]=0;
 #endif
       
       #ifdef EVOLVEELECTRONS
@@ -1527,7 +1540,7 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
   DBPutQuadvar1(file, "entropyinv","mesh1", entropyinv,
   		dimensions, ndim, NULL, 0, 
 		DB_INT, DB_ZONECENT, optList);
-
+  
   DBPutQuadvar1(file, "uint","mesh1", uint,
   		dimensions, ndim, NULL, 0, 
 		DB_DOUBLE, DB_ZONECENT, optList);
@@ -1540,9 +1553,6 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
   		dimensions, ndim, NULL, 0, 
 		DB_DOUBLE, DB_ZONECENT, optList);
 
-  DBPutQuadvar1(file, "omegaF","mesh1", OmegaF,
-  		dimensions, ndim, NULL, 0, 
-		DB_DOUBLE, DB_ZONECENT, optList);
 
   DBPutQuadvar1(file, "u0","mesh1", u0,
   		dimensions, ndim, NULL, 0, 
@@ -1557,32 +1567,20 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
   		dimensions, ndim, NULL, 0, 
 		DB_DOUBLE, DB_ZONECENT, optList);
 
-  DBPutQuadvar1(file, "lorentz","mesh1", lorentz,
-  		dimensions, ndim, NULL, 0, 
-		DB_DOUBLE, DB_ZONECENT, optList);
-
-#ifdef FORCEFREE
-  DBPutQuadvar1(file, "uint_ff","mesh1", uintff,
-  		dimensions, ndim, NULL, 0, 
-		DB_DOUBLE, DB_ZONECENT, optList);
-
-  DBPutQuadvar1(file, "u0_perp","mesh1", u0_perp,
-  		dimensions, ndim, NULL, 0, 
-		DB_DOUBLE, DB_ZONECENT, optList);
-  DBPutQuadvar1(file, "u1_perp","mesh1", u1_perp,
-  		dimensions, ndim, NULL, 0, 
-		DB_DOUBLE, DB_ZONECENT, optList);
-  DBPutQuadvar1(file, "u2_perp","mesh1", u2_perp,
+  
+  DBPutQuadvar1(file, "ffinv","mesh1", ffinv,
   		dimensions, ndim, NULL, 0, 
-		DB_DOUBLE, DB_ZONECENT, optList);
-  DBPutQuadvar1(file, "u3_perp","mesh1", u3_perp,
+		DB_INT, DB_ZONECENT, optList);
+  
+  DBPutQuadvar1(file, "lorentz","mesh1", lorentz,
   		dimensions, ndim, NULL, 0, 
 		DB_DOUBLE, DB_ZONECENT, optList);
-#endif
+  
   DBPutQuadvar1(file, "lorentz_perp","mesh1", lorentz_perp,
   		dimensions, ndim, NULL, 0, 
 		DB_DOUBLE, DB_ZONECENT, optList);
   DBPutQuadvar1(file, "lorentz_par","mesh1", lorentz_par,
+		
   		dimensions, ndim, NULL, 0, 
 		DB_DOUBLE, DB_ZONECENT, optList);
 
@@ -1742,6 +1740,7 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
 #endif
 
 #ifdef MAGNFIELD
+
   DBPutQuadvar1(file, "bsq","mesh1", bsq,
   		dimensions, ndim, NULL, 0, 
 		DB_DOUBLE, DB_ZONECENT, optList);
@@ -1755,6 +1754,10 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
   		dimensions, ndim, NULL, 0, 
 		DB_DOUBLE, DB_ZONECENT, optList);
 
+  DBPutQuadvar1(file, "OmegaF","mesh1", OmegaF,
+  		dimensions, ndim, NULL, 0, 
+		DB_DOUBLE, DB_ZONECENT, optList);
+
   DBPutQuadvar1(file, "phi","mesh1", phi,
   		dimensions, ndim, NULL, 0, 
 		DB_DOUBLE, DB_ZONECENT, optList);
@@ -1785,6 +1788,10 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
   		dimensions, ndim, NULL, 0, 
 		DB_DOUBLE, DB_ZONECENT, optList);
 
+  DBPutQuadvar1(file, "jdens","mesh1", jdens,
+  		dimensions, ndim, NULL, 0, 
+		DB_DOUBLE, DB_ZONECENT, optList);
+
   #ifdef MIMICDYNAMO
   DBPutQuadvar1(file, "phidyn","mesh1", phidyn,
   		dimensions, ndim, NULL, 0, 
@@ -1956,7 +1963,7 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
   DBPutQuadvar(file, "rad_vel","mesh1", 3, names, vels, 
   		dimensions, ndim, NULL, 0, 
 		DB_DOUBLE, DB_ZONECENT, optList);
-  #endif
+#endif
  
   // Close the Silo file and free memory
   DBClose(file);
@@ -1967,24 +1974,23 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
 
   free(rho);
   free(uint);
-  free(uintff);
   free(entr);
   free(temp);
   free(Omega);
   free(muBe);
+  free(Be);
   free(expansion);
   free(NHr);
   free(entrlnT);
-
+  free(entropyinv);
+  
   free(u0);
   free(u1);
   free(u2);
   free(u3);
   free(lorentz);
-  free(u0_perp);
-  free(u1_perp);
-  free(u2_perp);
-  free(u3_perp);
+
+  free(ffinv);
   free(lorentz_perp);
   free(lorentz_par);
   free(vx);
@@ -1994,25 +2000,6 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
   free(Edoty);
   free(Edotz);
 
-#ifdef EVOLVEELECTRONS
-  free(tempe);
-  free(tempi);
-  free(ue);
-  free(ui);
-#ifdef RELELECTRONS
-  free(gammabrk);
-  free(urelel);
-  free(nrelel); 
-  free(neth);
-  free(uratio_tot);
-  free(uratio_th);
-  free(G0relel);
-  free(G0icrelel);
-  free(G0synrelel);
-#endif
-#endif
-
-
 #ifdef MAGNFIELD
   free(Bangle);
   free(bsq);
@@ -2029,6 +2016,7 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
   free(sigmaw);
   free(betainv); 
   free(divB);
+  free(jdens);
   free(OmegaF);
   free(Qtheta);
   free(Qphi);
@@ -2064,18 +2052,37 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
 #ifdef PRINTGAMMATOSILO
   free(gammag);
 #endif
- 
+
+
+#ifdef EVOLVEELECTRONS
+  free(tempe);
+  free(tempi);
+  free(ue);
+  free(ui);
+#ifdef RELELECTRONS
+  free(gammabrk);
+  free(urelel);
+  free(nrelel); 
+  free(neth);
+  free(uratio_tot);
+  free(uratio_th);
+  free(G0relel);
+  free(G0icrelel);
+  free(G0synrelel);
+#endif
+#endif
+
+
 #ifdef RADIATION
+  free(trad);
+  free(tradlte);
+
   free(tausca);
   free(tauabs);
   free(taueff);
   free(taucoupling);
   free(Erad);
-  free(trad);
-  free(tradlte);
   free(Ehat);
-  free(Gtff);
-  free(G0icth);
   
   free(Fx);
   free(Fy);
@@ -2084,6 +2091,8 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
   free(uradx);
   free(urady);
   free(uradz);
+  free(Gtff);
+  free(G0icth);
 
   free(tauscar);
   free(tauabsr);
diff --git a/u2p.c b/u2p.c
index 400e88d..d64f9b0 100644
--- a/u2p.c
+++ b/u2p.c
@@ -155,16 +155,16 @@ u2p(ldouble *uu0, ldouble *pp, void *ggg, int corrected[3], int fixups[2], int t
   #endif
   gdetu_inv = 1. / gdetu;
   
-  int verbose=1;
+  int verbose=0;
   int hdcorr=0;
   int radcor=0;
   corrected[0]=corrected[1]=corrected[2]=0;
   fixups[0]=fixups[1]=0;
   
-  int u2pret,u2pentrret,ret;
+
   ldouble ppbak[NV];
-  for(u2pret=0;u2pret<NV;u2pret++)
-    ppbak[u2pret]=pp[u2pret];
+  for(iv=0;iv<NV;iv++)
+    ppbak[iv]=pp[iv];
 
   //************************************
   //magneto-hydro part
@@ -172,20 +172,23 @@ u2p(ldouble *uu0, ldouble *pp, void *ggg, int corrected[3], int fixups[2], int t
   ldouble u0=pp[UU];
   
   //************************************
-  //hot hydro - conserving energy
-  ret=0;
-  u2pret=-1;
+  // hot hydro - conserving energy
+  int ret=0;
+  int u2pret=-1;
     
-  //negative uu[RHO] = rho u^t -> set to floor 
-  if(uu[RHO] * gdetu_inv < 0.)
+  // negative uu[RHO] = rho u^t -> set rho to floor
+  ldouble alpha=geom->alpha;
+  //if(uu[RHO] * gdetu_inv < 0.)
+  if(uu[RHO] * alpha * gdetu_inv < 0.)    // ANDREW -- change, since uu[RHO]=rho*gamma*gdet/alpha. What if alpha<0? 
   {
     int gix,giy,giz;
     mpi_local2globalidx(geom->ix,geom->iy,geom->iz,&gix,&giy,&giz);
-    if(verbose)
+    if(1 || verbose)
       printf("%4d > %4d %4d %4d > neg rhout - requesting fixup\n",PROCID,gix,giy,giz);
 
     pp[RHO]=RHOFLOOR; 
-    uu[RHO]=RHOFLOOR*gdetu;
+    //uu[RHO]=RHOFLOOR*gdetu;
+    uu[RHO]=RHOFLOOR*gdetu/alpha; // ANDREW -- change, since uu[RHO]=rho*gamma*gdet/alpha. What if alpha<0? 
     ret=-2;  // to request fixup in all cases
     
     #ifndef SWAPPAPC
@@ -193,64 +196,109 @@ u2p(ldouble *uu0, ldouble *pp, void *ggg, int corrected[3], int fixups[2], int t
     #endif
   }
 
-  // TODO ANDREW MAKE THIS CLEANER
-#if defined(ENFORCEENTROPY) || (defined(FORCEFREE) && !defined(HYBRID_FORCEFREE))
-  printf("TESTTEST\n");
-  u2pret = -1;  //skip hot energy-conserving inversion and go straight to entropy inversion
-#else
-  u2pret = u2p_solver(uu,pp,ggg,U2P_HOT,0);  // invert using the hot energy equation    
-#endif //ENFORCEENTROPY
+  // force-free inversion
+  int ffflag=0;
+#ifdef FORCEFREE
+  ffflag = get_cflag(FFINVFLAG, geom->ix,geom->iy,geom->iz);
+#endif
 
-  if(ALLOWENTROPYU2P)  // Inversion with entropy equation -- on by default (see choices.h)
+  // use ff solver if not hybrid or if high sigma
+  if(ffflag>0)
+  {
+    u2pret = u2p_solver_ff(uu,pp,ggg,verbose);
+  }
+  else // MHD inversion
   {
-    if(u2pret<0)  // true if energy equation failed, or if energy equation was not required (because ENFORCEENTROPY is defined)
+
+    #ifdef ENFORCEENTROPY
+    u2pret = -1;  //skip hot energy-conserving inversion and go straight to entropy inversion
+    #else
+    u2pret = u2p_solver_mhd(uu,pp,ggg,U2P_HOT,0);  // invert using the hot energy equation    
+    #endif //ENFORCEENTROPY
+
+    if(ALLOWENTROPYU2P)  // Inversion with entropy equation -- on by default (see choices.h)
     {
-      ret=-1;
-      
-      if(verbose>2 )
+      if(u2pret<0)  // true if energy equation failed, or if energy equation was not required (because ENFORCEENTROPY is defined)
       {
-        printf("u2p_entr %d  > %d %d %d \n",
-	       u2pret, geom->ix + TOI, geom->iy + TOJ,geom->iz + TOK);
-      }
-      
-      //************************************
-      //entropy solver - invert using entropy equation
-      u2pret = u2p_solver(uu,pp,ggg,U2P_ENTROPY,0);
+        ret=-1;
       
-      if(u2pret<0)
-      {
-        ret=-2;
-        
-        if(verbose>1)
+        if(verbose>2 )
         {
-          printf("u2p_entr err %4d > %4d %4d %4d\n",
-		 u2pret,geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK);
+          printf("%4d > u2p_entr %d %d %d > %d %d %d \n",
+		 PROCID, u2pret, ret,ffflag,geom->ix + TOI, geom->iy + TOJ,geom->iz + TOK);
         }
+      
+        //entropy solver - invert using entropy equation
+        u2pret = u2p_solver_mhd(uu,pp,ggg,U2P_ENTROPY,0);
+      
+        if(u2pret<0)
+        {
+          ret=-2;
+        
+          if(verbose>1)
+          {
+            printf("%4d > u2p_entr ERROR %4d %d %d > %4d %4d %4d\n",
+		   PROCID,u2pret,ret,ffflag,geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK);
+          }
 	
-      } // if(u2pret<0) // second time -- entropy eqn
-    } // if(u2pret<0) // first time -- energy eqn
-  } // if(ALLOWENTROPYU2P)
+        } // if(u2pret<0) // second time -- entropy eqn
+      } // if(u2pret<0) // first time -- energy eqn
+    } // if(ALLOWENTROPYU2P)
+  }
+
+  // if entropy equation failed, check if we can try force-free
+#if defined(FORCEFREE) && defined(HYBRID_FORCEFREE)
+    if(ffflag==0 && u2pret<0)
+    {
+      ldouble ucon[4],ucov[4],bcon[4],bcov[4];
+      ldouble bsq;
+      calc_ucon_ucov_from_prims(pp, &geom, ucon, ucov);
+      calc_bcon_bcov_bsq_from_4vel(pp, ucon, ucov, &geom, bcon, bcov, &bsq);
+      if(bsq> pp[RHO]*HYBRID_FORCEFREE_SIGMACUT)
+      {
+	ffflag=1;
+	u2pret = u2p_solver_ff(uu,pp,ggg,verbose);
+        set_cflag(FFINVFLAG, geom->ix,geom->iy,geom->iz,ffflag);
+	if(u2pret>=0)
+	  ret=0;
+        if(verbose > 1)
+	  printf("%4d > MHD->FF %d %d %d > %4d %4d %4d \n",
+          PROCID,ret ,u2pret,ffflag,geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK);
+      }
+    }
+#endif
   
   if(u2pret<0)  // entropy equation also failed
   {
-    //leaving primitives unchanged
-    if(verbose > 1 || 1) // TODO ANDREW
+    
+    //leave primitives unchanged
+    if(1 || verbose) // TODO ANDREW
     {
-      printf("%4d > MHDU2PFAIL %d > %4d %4d %4d > u2p prim. unchanged \n",
-	     PROCID,u2pret,geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK);
+      printf("%4d > MHDU2PFAIL %d %d %d > %4d %4d %4d > u2p prim. unchanged \n",
+	     PROCID,u2pret,ret,ffflag,geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK);
     }
     ret=-3;
-    for(u2pret=0;u2pret<NV;u2pret++)
-      pp[u2pret]=ppbak[u2pret];
+    for(iv=0;iv<NV;iv++)
+      pp[iv]=ppbak[iv];
   }
   
-  if(ret<0)  // update conserved
+  if(ret<0) // update conserved
+  {
     hdcorr=1;  
+  }
   if(ret<-1) // request fixup when entropy failed
+  {
     fixups[0]=1;
+    if(verbose > 1)
+    {
+      printf("%4d > FIXUPS %d %d %d > %4d %4d %4d \n",
+             PROCID,u2pret,ret,ffflag,geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK);
+    }
+  }
   else
+  {
     fixups[0]=0;
-  
+  }
   //************************************
   //radiation part
   //************************************
@@ -471,7 +519,7 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg)
   
   //**********************************************************************
   //too magnetized
-#ifndef FORCEFREE //no floors for force-free evolution
+
 #ifdef MAGNFIELD
   ldouble ucon[4],ucov[4];
   ldouble bcon[4],bcov[4],bsq,magpre;
@@ -489,7 +537,6 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg)
   magpre = 0.5 * bsq;
   //magpre = bsq; // ANDREW: B2RHORATIOMAX and B2UURATIOMAX make more sense with bsq, not 0.5 bsq
 
-    
   // check density vs bsq
   if(bsq>B2RHORATIOMAX*pp[RHO]) 
   {
@@ -514,8 +561,13 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg)
     uufloored=1;
   }
 
+  int ffflag=0;
+  #ifdef FORCEFREE
+  ffflag = get_cflag(FFINVFLAG, geom->ix,geom->iy,geom->iz);
+  #endif
+    
   // floors were activated, apply them
-  if(rhofloored==1 || uufloored==1)
+  if((rhofloored==1 || uufloored==1) && ffflag==0)
   {
     // save backups
     for(iv=0;iv<NVMHD;iv++)
@@ -619,9 +671,9 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg)
 
     // find new prims after adding delta conserved in ZAMO
     int rettemp=0;
-    rettemp = u2p_solver(uu,pp,geom,U2P_HOT,0); 
+    rettemp = u2p_solver_mhd(uu,pp,geom,U2P_HOT,0); 
     if(rettemp<0)
-      rettemp = u2p_solver(uu,pp,geom,U2P_ENTROPY,0); 
+      rettemp = u2p_solver_mhd(uu,pp,geom,U2P_ENTROPY,0); 
       
     if(rettemp<0) // new mass in ZAMO frame failed
     {
@@ -724,7 +776,6 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg)
 
 
 #endif //MAGNFIELD
-#endif //ndef FORCEFREE
   
   //**********************************************************************
   //too fast
@@ -873,11 +924,6 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg)
   if(ret<0)
   {
     pp[ENTR]=calc_Sfromu(pp[RHO],pp[UU],geom->ix,geom->iy,geom->iz);
-    #ifdef FORCEFREE
-    //if(verbose>0) printf("updating ff in check_floors\n");
-    //update_ffprims_cell(pp, &geom);
-    //pp[UUFF] = pp[UU];
-    #endif
   }
 #endif //SKIPALLFLOORS
   
@@ -886,60 +932,27 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg)
 
 
 //**********************************************************************
-// solver wrapper
+// MHD solver wrapper
 //**********************************************************************
-//ANDREW TODO do we need to use function pointers?
 
 int
-u2p_solver(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose)
+u2p_solver_mhd(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose)
 {
-  struct geometry *geom
-      = (struct geometry *) ggg;
 
-#ifdef NONRELMHD
+  #ifdef NONRELMHD
   return u2p_solver_nonrel(uu,pp,ggg,Etype,verbose);
-#endif
-
-  int do_u2p_ff=0;
-#ifdef FORCEFREE
-  #ifndef HYBRID_FORCEFREE
-  do_u2p_ff=1;
-  #else
-  
-  //ANDREW TODO what should be the criterion here? local sigma from guess pp?
-  ldouble ucon0[4],ucov0[4],bcon0[4],bcov0[4];
-  ldouble bsq0,sigma0;
-  calc_ucon_ucov_from_prims(pp, geom, ucon0, ucov0);
-  calc_bcon_bcov_bsq_from_4vel(pp, ucon0, ucov0, geom, bcon0, bcov0, &bsq0);
-  sigma0 = bsq0/pp[RHO];
-  if(sigma0 > HYBRID_FORCEFREE_SIGMACUT)
-    do_u2p_ff = 1;
   #endif
-#endif //FORCEFREE
-
-  // use ff solver if not hybrid or if high sigma
-  if(do_u2p_ff)
-  {
-    return u2p_solver_ff(uu,pp,ggg,U2P_ENTROPY,verbose);
-  }
-  else
-  {
-  
-    //int (*solver)(ldouble*,ldouble*,void*,int,int);
   
-    #if (U2P_SOLVER==U2P_SOLVER_WP)
-    //solver = & u2p_solver_Wp;
-    return u2p_solver_Wp(uu,pp,ggg,Etype,verbose);
-    #endif
+  #if (U2P_SOLVER==U2P_SOLVER_WP)
+  return u2p_solver_Wp(uu,pp,ggg,Etype,verbose);
+  #endif
   
-    #if (U2P_SOLVER==U2P_SOLVER_W)  // this is the default
-    //solver = & u2p_solver_W;
-    return u2p_solver_W(uu,pp,ggg,Etype,verbose);
-    #endif
-  }
+  #if (U2P_SOLVER==U2P_SOLVER_W)  // this is the default
+  return u2p_solver_W(uu,pp,ggg,Etype,verbose);
+  #endif
+ 
+  return -200; // no solver found 
 
-  return -200;
-  //return (*solver)(uu,pp,ggg,Etype,verbose);
 } 
 
 
@@ -1257,7 +1270,8 @@ u2p_solver_W(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose)
   
   if(uint<0. || gamma2<0. || isnan(W) || !isfinite(W))
   {
-    if(verbose>0) printf("neg u in u2p_solver %e %e %e %e\n",rho,uint,gamma2,W);//getchar();
+    if(verbose>0) printf("neg u in u2p_solver %e %e %e %e\n",rho,uint,gamma2,W);
+    //getchar();
     return -104;
   }
   
@@ -2231,7 +2245,9 @@ test_inversion()
   fill_geometry(ix,iy,iz,&geom);
   ldouble alpha=geom.alpha;
   ldouble (*gg)[5];
+  ldouble (*GG)[5];
   gg=geom.gg;
+  GG=geom.GG;
   //fill_geometry_arb(ix,iy,iz,&geomBL,BLCOORDS);
 
   ldouble rho=1.3;
@@ -2291,14 +2307,14 @@ test_inversion()
 
   ldouble Sc = gamma*calc_Sfromu(rho,uint,ix,iy,iz);
   
-  printf("sigma %e\n",Bsq/(gamma_perp*gamma_perp)/rho);
+  printf("\nsigma %e\n",Bsq/(gamma_perp*gamma_perp)/rho);
   printf("Bsq %e\n",Bsq);
   printf("W %e\n",W);
   printf("D %e Y %e Z %e\n",D,Y,Z);
   printf("gammam2perp %e\n",1./gamma_perp/gamma_perp);
   printf("Sc %e \n",Sc);
   printf("check1 %e\n",Z-(p-W));
-  printf("check2 %e\n",(afac/(gamma*gamma)-1)*W - afac*D/gamma - Z);
+  printf("check2 %e\n\n",(afac/(gamma*gamma)-1)*W - afac*D/gamma - Z);
   
   // put primitives in
   pp[RHO]=rho;
@@ -2313,40 +2329,67 @@ test_inversion()
   pp[B2] = Bcon[2]/alpha;
   pp[B3] = Bcon[3]/alpha;
 
+
+  // look at electric field
+  ldouble ucon[4],ucov[4],bcon[4],bcov[4],bsq;
+  calc_ucon_ucov_from_prims(pp, &geom, ucon, ucov);
+  calc_bcon_bcov_bsq_from_4vel(pp, ucon, ucov, &geom, bcon, bcov, &bsq);
+  ldouble E1 = (bcov[2]*ucov[3]-bcov[3]*ucov[2])/geom.gdet;
+  ldouble E2 = (bcov[3]*ucov[1]-bcov[1]*ucov[3])/geom.gdet;
+  ldouble E3 = (bcov[1]*ucov[2]-bcov[2]*ucov[1])/geom.gdet;
+  printf("Efield: %e %e %e\n",E1,E2,E3);
+  
   #ifdef FORCEFREE
-  pp[UUFF] = uint;
+  pp[UUFF] = gamma_par*vpar; // ANDREW TODO ?? 
   pp[VXFF] = gamma_perp*vcon_perp[1];
   pp[VYFF] = gamma_perp*vcon_perp[2];
   pp[VZFF] = gamma_perp*vcon_perp[3];
   #endif
 
+  printf("RAISETEST\n");
+  ldouble ucon_velr[4];
+  ldouble ucov_velr[4];
+  ucon_velr[0] = 0;
+  ucon_velr[1] = pp[VX]; ucon_velr[2]=pp[VY]; ucon_velr[3]=pp[VZ];
+   
+
 
-  printf("input\n");
+  indices_21(ucon_velr,ucov_velr,gg);
+  printf("ucov: %e %e %e %e\n",ucov[0],ucov[1],ucov[2],ucov[3]);
+  printf("ucov_velr %e %e %e %e\n",ucov_velr[0],ucov_velr[1],ucov_velr[2],ucov_velr[3]);
+  printf("ucon: %e %e %e %e\n",ucon[0],ucon[1],ucon[2],ucon[3]);
+  printf("ucon_velr: %e %e %e %e\n",ucon_velr[0],ucon_velr[1],ucon_velr[2],ucon_velr[3]);
+  indices_12(ucov_velr,ucon_velr,GG);
+  printf("ucon_velr 2: %e %e %e %e\n",ucon_velr[0],ucon_velr[1],ucon_velr[2],ucon_velr[3]);
+  
+  printf("\n u2p input\n");
   print_primitives(pp);
   p2u(pp,uu,&geom);
 
+  printf("uu input\n");
+  print_conserved(uu);
+  
   PLOOP(iv) pp0[iv]=pp[iv];
 
+
   pp[RHO]*=3.1245325124;
   pp[UU]*=23.124124214421124;
-  //pp[ENTR] = calc_Sfromu(pp[RHO],pp[UU],ix,iy,iz);
+  pp[ENTR] = calc_Sfromu(pp[RHO],pp[UU],ix,iy,iz);
   
   pp[VX]*=12.3; 
   pp[VY]*=21.1;
   pp[VZ]*=.66;
   
-
   #ifdef FORCEFREE
-  pp[UUFF]=pp[UU];
+  pp[UUFF]*=.8;
   pp[VXFF]*=22.3;
   pp[VYFF]*=12.1;
   pp[VZFF]*=5.6;
   #endif
 
-  //print_conserved(uu);
-  u2p_solver_ff(uu,pp,&geom,U2P_ENTROPY,2);
-  //u2p_solver(uu,pp,&geom,U2P_HOT,0);
-
+  
+  u2p_solver_ff(uu,pp,&geom,2);
+  
   printf("solved\n");
   print_primitives(pp);
 
@@ -2359,7 +2402,7 @@ test_inversion()
   printf("projection test....\n");
   PLOOP(iv) pp[iv]=pp0[iv];
   fill_ffprims_cell(pp, &geom);
-  print_primitives(pp);
+  //print_primitives(pp);
   PLOOP(iv) pperr[iv]=(pp[iv]-pp0[iv])/pp0[iv];
   printf("projection error\n");
   print_primitives(pperr);
@@ -2411,7 +2454,7 @@ test_inversion_nonrel()
   print_conserved(uu);
   printf("gdet = %e\n",geom.gdet);
 
-  u2p_solver(uu,pp,&geom,U2P_HOT,0); 
+  u2p_solver_mhd(uu,pp,&geom,U2P_HOT,0); 
   print_primitives(pp);
 
   return 0;