diff --git a/.gitignore b/.gitignore
index 5ce4f07..c38dcf1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,4 +2,6 @@
 *.o
 *.orig
 #*
-*.nb
\ No newline at end of file
+*.nb
+.DS_Store
+**/.DS_Store
diff --git a/PROBLEMS/PUFFY/bc.c b/PROBLEMS/PUFFY/bc.c
new file mode 100644
index 0000000..3f3dd38
--- /dev/null
+++ b/PROBLEMS/PUFFY/bc.c
@@ -0,0 +1,236 @@
+//returns problem specific BC
+//int calc_bc(int ix,int iy,int iz,ldouble t,
+//	ldouble *uu,ldouble *pp,int ifinit,int BCtype)
+
+/**********************/
+//geometries
+ldouble gdet_src,gdet_bc,Fx,Fy,Fz,ppback[NV];
+int iix,iiy,iiz,iv;  	  
+
+struct geometry geom;
+fill_geometry(ix,iy,iz,&geom);
+
+struct geometry geomBL;
+fill_geometry_arb(ix,iy,iz,&geomBL,BLCOORDS);
+
+ldouble gg[4][5],GG[4][5],ggsrc[4][5],eup[4][4],elo[4][4];
+pick_g(ix,iy,iz,gg);
+pick_G(ix,iy,iz,GG);
+
+/**********************/
+
+//radius
+if(BCtype==XBCHI) //outflow in magn, atm in rad., atm. in HD
+  {
+    iix=NX-1;
+    iiy=iy;
+    iiz=iz;
+    
+    //copying everything
+    for(iv=0;iv<NV;iv++)
+      {
+	pp[iv]=get_u(p,iv,iix,iiy,iiz);
+      }
+    
+    //!! begin rescale
+    //first transform to BL
+    trans_pmhd_coco(pp, pp, MYCOORDS,BLCOORDS, geom.xxvec,&geom,&geomBL);
+    
+    struct geometry geombdBL;
+    fill_geometry_arb(iix,iiy,iiz,&geombdBL,BLCOORDS);
+    ldouble rghost = geomBL.xx;
+    ldouble rbound = geombdBL.xx;
+    ldouble deltar = rbound - rghost;
+    ldouble scale1 =  rbound*rbound/rghost/rghost;
+    ldouble scale2 = rbound/rghost;
+
+    pp[RHO]*=scale1;
+    pp[UU] *=scale1;
+    #ifdef MAGNFIELD
+    pp[B1] *=scale1;
+    pp[B2] *=scale2;
+    pp[B3] *=scale2;
+    #endif
+    //pp[VX] *=1.;
+    pp[VY] *=scale1;
+    pp[VZ] *=scale1;
+    #ifdef RADIATION
+    pp[EE0] *=scale1;
+    //pp[FX0] *=1.;
+    pp[FY0] *=scale1;
+    pp[FZ0] *=scale1;
+    #endif
+    //transform back after rescaling
+    trans_pmhd_coco(pp, pp,BLCOORDS, MYCOORDS, geomBL.xxvec,&geomBL, &geom);
+    //!! end rescale
+
+    //checking for the gas inflow
+    ldouble ucon[4]={0.,pp[VX],pp[VY],pp[VZ]};    
+    conv_vels(ucon,ucon,VELPRIM,VEL4,geom.gg,geom.GG);
+    if(MYCOORDS!=CYLCOORDS) trans2_coco(geom.xxvec,ucon,ucon,MYCOORDS,BLCOORDS);
+    if(ucon[1]<0.) //inflow, resseting to atmosphere
+      {
+	//atmosphere in rho,uint and velocities and zero magn. field
+	//set_hdatmosphere(pp,xxvec,gg,GG,4);
+	//set_hdatmosphere(pp,geom.xxvec,geom.gg,geom.GG,0);
+
+	//!! begin reset to floor        
+	//	pp[UU] = pp[RHO]*UURHORATIOMIN*3.;
+	
+	//!! end reset to floor
+
+	ucon[1]=0.;
+	/**
+	#ifdef MAGNFIELD
+	pp[B2]=pp[B3]=0.;
+	#endif
+	**/
+	if(MYCOORDS!=CYLCOORDS) trans2_coco(geomBL.xxvec,ucon,ucon,BLCOORDS,MYCOORDS);
+	conv_vels(ucon,ucon,VEL4,VELPRIM,geom.gg,geom.GG);
+	pp[VX]=ucon[1];
+	pp[VY]=ucon[2];
+	pp[VZ]=ucon[3];//atmosphere in rho,uint and velocities and zero magn. field
+      }
+    
+#ifdef RADIATION
+    ldouble urfcon[4]={0.,pp[FX0],pp[FY0],pp[FZ0]};    
+    conv_vels(urfcon,urfcon,VELPRIMRAD,VEL4,geom.gg,geom.GG);
+    if(MYCOORDS!=CYLCOORDS) trans2_coco(geom.xxvec,urfcon,urfcon,MYCOORDS,BLCOORDS);
+    if(urfcon[1]<0.) //inflow, resseting to atmosphere
+      {
+	//!! begin reset to floor
+	//pp[EE0] = ERADATMMIN;
+	//!! end reset to floor
+
+	//atmosphere in radiation
+	//set_radatmosphere(pp,xxvec,gg,GG,0);
+	urfcon[1]=0.;
+	if(MYCOORDS!=CYLCOORDS) trans2_coco(geomBL.xxvec,urfcon,urfcon,BLCOORDS,MYCOORDS);
+	conv_vels(urfcon,urfcon,VEL4,VELPRIMRAD,geom.gg,geom.GG);
+	pp[FX0]=urfcon[1];
+	pp[FY0]=urfcon[2];
+	pp[FZ0]=urfcon[3];//atmosphere in rho,uint and velocities and zero magn. field
+      }
+#endif
+
+   
+    
+    
+    p2u(pp,uu,&geom);
+    return 0;   
+  }
+ else if(BCtype==XBCLO) //outflow near BH
+   {
+     iix=0;
+     iiy=iy;
+     iiz=iz;
+
+      for(iv=0;iv<NV;iv++)
+       {
+	 pp[iv]=get_u(p,iv,0,iiy,iiz);
+       }
+
+      if(RMIN>rhorizonBL)
+	{
+	  //checking for the gas inflow
+	  ldouble ucon[4]={0.,pp[VX],pp[VY],pp[VZ]};    
+	  conv_vels(ucon,ucon,VELPRIM,VEL4,geom.gg,geom.GG);
+	  trans2_coco(geom.xxvec,ucon,ucon,MYCOORDS,BLCOORDS);
+	  if(ucon[1]<0.) //inflow, resseting to atmosphere
+	    {
+	      //atmosphere in rho,uint and velocities and zero magn. field
+	      //set_hdatmosphere(pp,xxvec,gg,GG,4);
+	      ucon[1]=0.;
+	      trans2_coco(geomBL.xxvec,ucon,ucon,BLCOORDS,MYCOORDS);
+	      conv_vels(ucon,ucon,VEL4,VELPRIM,geom.gg,geom.GG);
+	      pp[VX]=ucon[1];
+	      pp[VY]=ucon[2];
+	      pp[VZ]=ucon[3];//atmosphere in rho,uint and velocities and zero magn. field
+	    }
+
+#ifdef RADIATION
+	  ldouble urfcon[4]={0.,pp[FX0],pp[FY0],pp[FZ0]};    
+	  conv_vels(urfcon,urfcon,VELPRIMRAD,VEL4,geom.gg,geom.GG);
+	  trans2_coco(geom.xxvec,urfcon,urfcon,MYCOORDS,BLCOORDS);
+	  if(urfcon[1]<0.) //inflow, resseting to atmosphere
+	    {
+	      //atmosphere in radiation
+	      //set_radatmosphere(pp,xxvec,gg,GG,0);
+	      urfcon[1]=0.;
+	      trans2_coco(geomBL.xxvec,urfcon,urfcon,BLCOORDS,MYCOORDS);
+	      conv_vels(urfcon,urfcon,VEL4,VELPRIMRAD,geom.gg,geom.GG);
+	      pp[FX0]=urfcon[1];
+	      pp[FY0]=urfcon[2];
+	      pp[FZ0]=urfcon[3];//atmosphere in rho,uint and velocities and zero magn. field
+	    }
+#endif
+	}
+
+     p2u(pp,uu,&geom);
+     return 0;
+   }
+
+//reflections/outflow in theta 
+//in 3D polar cells overwritten with #define CORRECT_POLARAXIS_3D
+if(BCtype==YBCLO) //upper spin axis 
+  {      
+    
+    iiy=-iy-1;
+    iiz=iz;
+    iix=ix;
+    gdet_src=get_g(g,3,4,iix,iiy,iiz);  
+    gdet_bc=get_g(g,3,4,ix,iy,iz);  
+    for(iv=0;iv<NV;iv++)
+      {
+	//v_theta
+	if(iv==VY || iv==B2 || iv==FY0)
+	  pp[iv]=-get_u(p,iv,iix,iiy,iiz);
+	else
+	  pp[iv]=get_u(p,iv,iix,iiy,iiz);
+       }
+     
+
+    p2u(pp,uu,&geom);
+    return 0;
+  }
+
+if(BCtype==YBCHI) //lower spin axis
+  {
+    iiy=NY-(iy-NY)-1;
+    iiz=iz;
+    iix=ix;
+    gdet_src=get_g(g,3,4,iix,iiy,iiz);  
+    gdet_bc=get_g(g,3,4,ix,iy,iz);  
+  	  
+    for(iv=0;iv<NV;iv++)
+      {
+	if(iv==VY || iv==B2 || iv==FY0)
+	  pp[iv]=-get_u(p,iv,iix,iiy,iiz);
+	else
+	  pp[iv]=get_u(p,iv,iix,iiy,iiz);
+	
+      }
+ 
+    p2u(pp,uu,&geom); 
+    return 0; 
+  }
+   
+//periodic in phi:
+iiz=iz;
+iiy=iy;
+iix=ix;
+if(BCtype==ZBCLO) iiz=iz+NZ;
+if(BCtype==ZBCHI) iiz=iz-NZ;
+
+for(iv=0;iv<NV;iv++)
+  {
+    uu[iv]=get_u(u,iv,iix,iiy,iiz);
+    pp[iv]=get_u(p,iv,iix,iiy,iiz);      
+  }
+
+//testing if interpolated primitives make sense
+//check_floors_hd(pp,VELPRIM,gg,GG);
+//end of floor section
+ 
+return 0;
+
diff --git a/PROBLEMS/PUFFY/define.h b/PROBLEMS/PUFFY/define.h
new file mode 100644
index 0000000..6de9292
--- /dev/null
+++ b/PROBLEMS/PUFFY/define.h
@@ -0,0 +1,244 @@
+//unfortunately this is what it was
+#define MU_GAS 1.
+#define MU_I 2.
+#define MU_E 2.
+
+
+/************************************/
+//general
+/************************************/
+#define BHDISK_PROBLEMTYPE
+
+/************************************/
+//restart
+/************************************/
+#define RESTART
+#define RESTARTGENERALINDICES
+#define RESTARTNUM 0
+
+/************************************/
+//radiation choices
+/************************************/
+#define RADIATION
+// does this work? 
+#define BALANCEENTROPYWITHRADIATION
+
+#define COMPTONIZATION
+#define U2PCONV 1.e-10
+#define RADIMPLICITDAMPINGFACTOR 3.
+#define RADIMPLICITMAXENCHANGEDOWN 10.
+#define RADIMPLICIMAXENCHANGEUP 10.
+#define MAXRADIMPDAMPING 1.e-3
+#define RADIMPCONV 1.e-10
+#define RADIMPEPS 1.e-6
+#define RADIMPMAXITER 40
+#define RADIMPCONVREL 1.e-8
+#define RADIMPCONVRELERR 1.e-4
+#define RADIMPCONVRELENTR 1.e-4
+#define RADIMPCONVRELENTRERR 1.e-2
+// how does the opacity damping work?
+#define OPDAMPINIMPLICIT 0 
+
+//no bremsstrahlung in AGN
+#define BREMSSTRAHLUNG
+#define SYNCHROTRON
+//Should be used, the bridge function is adding a NR component to synchotron opacity
+#define NO_SYNCHROTRON_BRIDGE_FUNCTIONS
+#define KLEINNISHINA
+
+// try - from Brandon
+#define RADIMP_START_WITH_BISECT
+#define NPH_START_WITH_BISECT
+#define ALLOWRADCEILINGINIMPLICIT
+#define ALLOWFORENTRINF4DPRIM
+
+#define RADIMPLICITMAXNPHCHANGE 100 // default is 1e10 - goes to nph floors??? strange, check
+
+//#define SKIPFANCYOPACITIES
+
+//#define EVOLVEPHOTONNUMBER
+//we need this for evolvephotonnumber
+#define SCALE_JACOBIAN
+
+/************************************/
+//magnetic choices
+/************************************/
+#define MIMICDYNAMO
+#define CALCHRONTHEGO
+#define THETAANGLE 0.25
+#define ALPHAFLIPSSIGN                                                        
+
+#define ALPHADYNAMO (2.*0.314)
+#define DAMPBETA
+#define BETASATURATED 0.1
+#define ALPHABETA (2.*6.28)
+#define MAGNFIELD
+#define GDETIN 1
+#define VECPOTGIVEN
+#define MAXBETA .01 //target pmag/pgas int the midplane of the initial torus, overwritten later
+
+/************************************/
+//reconstruction / Courant
+/************************************/
+#define INT_ORDER 2 // 1 - lets try better (bfrom Brandon's problem)
+#define TIMESTEPPING RK2IMEX
+#define TSTEPLIM .5
+#define FLUXLIMITER 0
+#define MINMOD_THETA 1.5
+#define DOFIXUPS 1
+#define DOU2PRADFIXUPS 0
+#define DOU2PMHDFIXUPS 1
+#define DORADIMPFIXUPS 0 // 1 let's try without ... 
+
+/************************************/
+//viscosity choices
+/************************************/
+#define RADVISCOSITY SHEARVISCOSITY
+//#define ACCELRADVISCOSITY // not great ? better calc. each step - if the rad. velocities are not changed much since the last step, it will use the old values
+#define RADVISCMFPSPH
+#define RADVISCNUDAMP
+#define RADVISCMAXVELDAMP
+#define ALPHARADVISC 0.1
+#define MAXRADVISCVEL 0.1
+
+/************************************/
+//rmhd floors
+/************************************/
+#define CORRECT_POLARAXIS
+//#define POLARAXISAVGIN3D
+#define NCCORRECTPOLAR 2
+
+// we are hitting 1e-8
+#define UURHORATIOMIN 1.e-8 // 1.e-10 new from Brandon
+#define UURHORATIOMAX 1.e0 // 1.e2 new from Brandon
+
+#define EERHORATIOMIN 1.e-20
+#define EERHORATIOMAX 1e4 //1.e20 - can go even lower after regridding to force it to relax
+
+#define EEUURATIOMIN 1.e-20
+#define EEUURATIOMAX 1.e20
+
+
+#define B2UURATIOMAX 50. //around the same as uu to rho
+// maybe try this (after it is implemented correctly) to avoid breaking of the simulations - can be also usefull after regridding 
+// when we want to force it to run until it relaxes a bit (eg after regridding)
+//#define B2RHOFLOOR_BACKUP_FFFRAME
+#define B2RHOFLOORFRAME DRIFTFRAME
+//#define DUINTFROMDRHO - in ZAMOFRAME, maybe better than duint=0
+
+
+#define B2RHORATIOMAX 50.
+
+#define GAMMAMAXRAD 10.
+#define GAMMAMAXHD 10.
+//maybe? 
+#define RHOFLOOR 1.e-30
+// #define UUEERATIOMAX 10. // David - trying to remove the extra uu from ehat
+
+/************************************/
+//blackhole
+/************************************/
+#define MASS 10.
+#define BHSPIN 0.
+
+/************************************/
+//coordinates / resolution
+/************************************/
+#define myMKS2COORDS
+#define MKSR0 0.1
+#define MKSH0 0.9
+#define MKSMY1 0.001
+#define MKSMY2 0.2
+#define MKSMP0 1.5
+#define METRICAXISYMMETRIC
+#define RMIN 1.85
+#define RMAX 500.
+
+#ifdef myMKS2COORDS //modified Kerr-Shild
+#define MYCOORDS MKS2COORDS
+#define MINX (log(RMIN-MKSR0))
+#define MAXX (log(RMAX-MKSR0))
+#define MINY (0.001)
+#define MAXY (1.-0.001)
+#endif
+
+#ifdef myMKS3COORDS //modified Kerr-Shild further from axis
+#define METRICNUMERIC
+#define MYCOORDS MKS3COORDS
+#define MINX (log(RMIN-MKSR0))
+#define MAXX (log(RMAX-MKSR0))
+#define MINY 0.
+#define MAXY 1.
+#endif
+
+#define PHIWEDGE (M_PI/2.)
+#define MINZ (-PHIWEDGE/2.)
+#define MAXZ (PHIWEDGE/2.)
+
+//total resolution
+#define TNX 384 //32*10
+#define TNY 360 //64*5
+#define TNZ 1 //32 //2*8
+//number of tiles
+#define NTX 32
+#define NTY 30
+#define NTZ 1
+
+#define SPECIFIC_BC
+#define PERIODIC_ZBC
+
+/************************************/
+//output
+/************************************/
+#define DTOUT1 100.
+#define DTOUT2 500.
+
+
+#define BOXOUTPUT 0
+
+
+#define OUTCOORDS BLCOORDS                                                                    
+#define OUTVEL VEL4
+#define ALLSTEPSOUTPUT 0
+#define NOUTSTOP 5000
+#define SILOOUTPUT 0
+#define OUTOUTPUT 0
+#define COORDOUTPUT 0
+#define SIMOUTPUT 0
+#define RADOUTPUT 0 
+#define SCAOUTPUT 1
+#define AVGOUTPUT 0 
+#define THOUTPUT 0
+#define THPROFRAD1US 30
+#define SILO2D_XZPLANE
+
+#define PRINT_FIXUPS_TO_SILO
+
+/************************************/
+//common physics / torus / atmosphere
+/************************************/
+#define RHOATMMIN  1.e-24
+#define UINTATMMIN  (calc_PEQ_ufromTrho(1.e10,RHOATMMIN,0,0,0))
+#define ATMTRADINIT 3.e6
+//define ERADATMMIN  (calc_LTE_EfromT(3.e6)/10) //(calc_LTE_EfromT(ATMTRADINIT)/10)
+// From Brandon's collapsing sims
+#define ERADATMMIN  (calc_LTE_EfromT(3.e5)/10*6.62/MASS) //If this is set too high, the torus is not in hydrostatic equilibrium since torus will take on atm value
+
+#define GAMMA (5./3.)
+
+#define QUADLOOPS
+
+#define RESTORETORUS
+//#define RESCALEDENSITY 0.9
+//#define RESCALETORUS 0.8
+
+#define LT_KAPPA 6.e1
+#define EXPECTEDHR 0.3
+#define LT_XI 0.995
+#define LT_R1 20.
+#define LT_R2 350.
+#define LT_GAMMA 4./3.
+#define LT_RIN 35.
+#define BETANORMFULL
+#undef MAXBETA
+#define MAXBETA (1./20.) 
diff --git a/PROBLEMS/PUFFY/dump.c b/PROBLEMS/PUFFY/dump.c
new file mode 100644
index 0000000..893bf25
--- /dev/null
+++ b/PROBLEMS/PUFFY/dump.c
@@ -0,0 +1,34 @@
+//v1 - 24
+//..
+//v7 - 30
+
+ldouble bcond[4],bcovd[4],bsqd;
+
+
+#ifdef MAGNFIELD
+calc_bcon_prim(pp,bcond,&geom);
+indices_21(bcond,bcovd,gg); 
+bsqd = dot(bcond,bcovd);
+v1=bsqd/2.;
+v2=calc_divB(ix,iy,iz);
+#endif
+
+#ifdef RADIATION
+ldouble Rtt,Ehat,ucon[4],prad;
+calc_ff_Rtt(pp,&Rtt,ucon,&geom);
+Ehat=-Rtt; 
+prad=Ehat/3.;
+v3=prad;
+
+//tau in radius
+ldouble tau[3];
+calc_tautot(pp,xxvec,dx,tau);
+v4=tau[1];
+calc_tauabs(pp,xxvec,dx,tau);
+v5=tau[1];
+#endif
+
+//a'la plasma beta : pmag/(pgas+prad)
+v6=v1/(GAMMAM1*pp[UU] + v3);
+
+
diff --git a/PROBLEMS/PUFFY/finger.c b/PROBLEMS/PUFFY/finger.c
new file mode 100644
index 0000000..e5991d9
--- /dev/null
+++ b/PROBLEMS/PUFFY/finger.c
@@ -0,0 +1,77 @@
+//drives the density towards initial state
+int ix,iy,iz,ii,iv;
+/**************************/
+
+#pragma omp parallel for private(ix,iy,iz,iv) schedule (dynamic)
+for(ii=0;ii<Nloop_0;ii++) //domain only
+  {
+    ldouble pp[NV],pptorus[NV],uu[NV];
+    struct geometry geom;
+    ix=loop_0[ii][0];
+    iy=loop_0[ii][1];
+    iz=loop_0[ii][2];
+    fill_geometry(ix,iy,iz,&geom);
+
+    ldouble xxBL[4];
+    coco_N(geom.xxvec,xxBL,MYCOORDS, BLCOORDS);
+    
+    //Keplerian orbital period
+    ldouble Omk = 1./(BHSPIN+sqrt(xxBL[1]*xxBL[1]*xxBL[1]));
+    ldouble Pk = 2.*M_PI/Omk;
+    
+    for(iv=0;iv<NV;iv++)
+      {
+	pp[iv]=get_u(p,iv,ix,iy,iz);
+	pptorus[iv]=get_u(ptemp1,iv,ix,iy,iz);
+      }
+
+    //density driven toward initial 
+    //temperature fixed
+    ldouble drho=0.;
+    ldouble duint=0;
+    ldouble dvr=0.;
+    ldouble dvth=0.;
+    ldouble dvph=0.;
+
+    #ifdef RESTORETORUS
+
+    ldouble effscale=1.; //how fast restores
+    if(pptorus[ENTR]>0.) //inside initial torus
+      {
+	drho =  - effscale
+	  * global_dt / Pk 
+	  * (pp[RHO]-pptorus[RHO]);
+
+	//not to overshoot 0.
+	while (drho+pp[RHO]<0.)
+	  {
+	    drho/=2.;
+	  }
+
+	//drive towards torus parameters
+	duint = (pp[UU]*pp[RHO]+pptorus[UU]*drho)/(pp[RHO]+drho) - pp[UU];
+	dvr = (pp[VX]*pp[RHO]+pptorus[VX]*drho)/(pp[RHO]+drho) - pp[VX];
+	dvth = (pp[VY]*pp[RHO]+ pptorus[VY]*drho)/(pp[RHO]+drho) - pp[VY];
+	dvph = (pp[VZ]*pp[RHO] + pptorus[VZ]*drho)/(pp[RHO]+drho) - pp[VZ];
+      }
+
+    //if(iy==NY/2 && ix==NX/2)
+    //printf("%e %e %e\n",drho,pp[RHO],pptorus[RHO]);
+
+#endif
+
+    //modifying only density, leaving temperature and velocities and magn. field intact
+    pp[RHO]+=drho;
+    pp[UU]+=duint;
+    pp[VX]+=dvr;
+    pp[VY]+=dvth;
+    pp[VZ]+=dvph;
+
+    p2u(pp,uu,&geom);
+    
+    for(iv=0;iv<NV;iv++)
+      {
+	set_u(u,iv,ix,iy,iz,uu[iv]);
+	set_u(p,iv,ix,iy,iz,pp[iv]);
+      }
+}
diff --git a/PROBLEMS/PUFFY/init.c b/PROBLEMS/PUFFY/init.c
new file mode 100644
index 0000000..6203b90
--- /dev/null
+++ b/PROBLEMS/PUFFY/init.c
@@ -0,0 +1,23 @@
+ldouble pp[NV], uu[NV];
+  
+//geometries
+
+struct geometry geom;
+fill_geometry(ix,iy,iz,&geom);
+
+int iv;
+for(iv=0;iv<NV;iv++)
+  pp[iv]=get_u(ptemp1,iv,ix,iy,iz);
+pp[ENTR]=calc_Sfromu(pp[RHO],pp[UU],ix,iy,iz);
+
+p2u(pp,uu,&geom);
+
+for(iv=0;iv<NV;iv++)
+  {
+    set_u(u,iv,ix,iy,iz,uu[iv]);
+    set_u(p,iv,ix,iy,iz,pp[iv]);
+  }
+
+//entropy
+//update_entropy(ix,iy,iz,0);
+//set_cflag(0,ix,iy,iz,0);
diff --git a/PROBLEMS/PUFFY/kappa.c b/PROBLEMS/PUFFY/kappa.c
new file mode 100644
index 0000000..e6d9acb
--- /dev/null
+++ b/PROBLEMS/PUFFY/kappa.c
@@ -0,0 +1,13 @@
+//absorption
+
+ldouble rhocgs=rhoGU2CGS(rho);
+ldouble Tcgs=tempGU2CGS(Tgas);
+ldouble kappaffcgs=6.4e22*rhocgs/Tcgs/Tcgs/Tcgs/sqrt(Tcgs);
+
+//test
+//return 0.;
+kappa= kappaCGS2GU(kappaffcgs)*rho;
+
+//ldouble kappabfcgs=4.8e-24/1.67262158e-24/1.67262158e-24*rhocgs/Tcgs/Tcgs/Tcgs/sqrt(Tcgs)*ZZsun;
+//return (kappaCGS2GU(kappaffcgs)*rho+kappaCGS2GU(kappabfcgs)*rho);
+
diff --git a/PROBLEMS/PUFFY/kappaes.c b/PROBLEMS/PUFFY/kappaes.c
new file mode 100644
index 0000000..7dfd178
--- /dev/null
+++ b/PROBLEMS/PUFFY/kappaes.c
@@ -0,0 +1,16 @@
+//scattering
+  
+ldouble OpacCorrection = 1;
+
+
+
+#ifdef KLEINNISHINA
+ldouble Tkn;
+Tkn = Trad;
+OpacCorrection = 1/(1 + pow(Tkn/4.5e8,0.86) );
+#endif
+
+//return  kappaCGS2GU(0.34)*rho; //original one
+
+return  kappaCGS2GU(0.2*(1. + HFRAC)*OpacCorrection)*rho;
+
diff --git a/PROBLEMS/PUFFY/limotorus.c b/PROBLEMS/PUFFY/limotorus.c
new file mode 100644
index 0000000..701d98f
--- /dev/null
+++ b/PROBLEMS/PUFFY/limotorus.c
@@ -0,0 +1,368 @@
+#include <stdio.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_integration.h>
+
+#define FTYPE double
+
+#define LT_BHSPIN 0.9
+#define LT_KAPPA 2.e3
+#define LT_XI 0.975
+#define LT_R1 30.
+#define LT_R2 200.
+#define LT_GAMMA 4./3.
+#define LT_RIN 22.
+
+void compute_gd( FTYPE r, FTYPE th, FTYPE a, FTYPE *gdtt, FTYPE *gdtp, FTYPE *gdpp ) {
+   FTYPE Sigma, tmp;
+
+   Sigma = (pow(a*cos(th),2) + pow(r,2));
+   tmp = 2*r*pow(sin(th),2)/Sigma;
+
+   //metric (expressions taken from limotorus4.nb):
+   *gdtt = -1 + 2*r/Sigma;
+   *gdtp = -a*tmp;
+   *gdpp = (pow(r,2) + pow(a,2)*(1+tmp)) * pow(sin(th),2);
+}
+
+// Keplerian, equatorial angular momentum density: l = u_phi / u_t
+FTYPE lK(FTYPE r, FTYPE a) {
+   FTYPE curlyF, curlyG;
+
+   curlyF = 1 - 2*a/pow(r, 1.5) + pow(a/r, 2);
+   curlyG = 1 - 2/r + a/pow(r, 1.5);
+   return( pow(r,0.5) * curlyF / curlyG );
+}
+
+FTYPE l3d(FTYPE lam, FTYPE a, FTYPE lambreak1, FTYPE lambreak2, FTYPE xi) {
+   return ( xi * lK( lam<=lambreak1 ? lambreak1 : lam>=lambreak2 ? lambreak2 : lam , a) );
+}
+
+FTYPE rtbis(FTYPE (*func)(FTYPE,FTYPE*), FTYPE *parms, FTYPE x1, FTYPE x2, FTYPE xacc)
+//Taken from HARM:nrutil.c
+{
+   int j;
+   FTYPE dx,f,fmid,xmid,rtb;
+   f=(*func)(x1, parms);
+   fmid=(*func)(x2, parms);
+   if (f*fmid >= 0.0) {
+      printf( "f(%g)=%g f(%g)=%g\n", x1, f, x2, fmid );
+      printf("Root must be bracketed for bisection in rtbis\n");
+   }
+   rtb = (f < 0.0) ? (dx=x2-x1,x1) : (dx=x1-x2,x2); //Orient the search so that f>0 lies at x+dx.
+   for (j=1;j<=100;j++) {
+      fmid=(*func)(xmid=rtb+(dx *= 0.5),parms); //Bisection loop.
+      if (fmid <= 0.0) {
+         rtb=xmid;
+      }
+      if (fabs(dx) < xacc || fmid == 0.0) {
+         return rtb;
+      }
+   }
+   printf("Too many bisections in rtbis\n");
+   return 0.0; //Never get here.
+}
+
+FTYPE lamBL_func(FTYPE lam, FTYPE *parms) {
+   FTYPE gdtt, gdtp, gdpp, a, lambreak1, lambreak2, xi, l;
+
+   gdtt      = parms[0];
+   gdtp      = parms[1];
+   gdpp      = parms[2];
+   a         = parms[3];
+   lambreak1 = parms[4];
+   lambreak2 = parms[5];
+   xi        = parms[6];
+
+   l = l3d(lam, a, lambreak1, lambreak2, xi);
+
+   return ( lam*lam + l * (l*gdtp + gdpp) / (l*gdtt + gdtp) );
+}
+
+// von Zeipel cylinder radius (needs to be calculated iteratively for nonzero spin)
+FTYPE lamBL(FTYPE R, FTYPE gdtt, FTYPE gdtp, FTYPE gdpp, FTYPE a, FTYPE lambreak1, FTYPE lambreak2, FTYPE xi) {
+   // R = r*sin(th), used as initial guess for lamBL
+
+   FTYPE parms[7];
+
+   //store params in an array before function call
+   parms[0] = gdtt;
+   parms[1] = gdtp;
+   parms[2] = gdpp;
+   parms[3] = a;
+   parms[4] = lambreak1;
+   parms[5] = lambreak2;
+   parms[6] = xi;
+
+   //solve for lin using bisection, specify large enough root search range, (1e-3, 1e3)
+   //demand accuracy 5x machine prec.
+   //in non-rel limit rml = r , use 10x that as the upper limit:
+   return( rtbis( &lamBL_func, parms, R, 10*R, 5.*DBL_EPSILON ) );
+
+}
+
+FTYPE omega3d( FTYPE l, FTYPE gdtt, FTYPE gdtp, FTYPE gdpp ) {
+   return( -(gdtt*l + gdtp)*pow(gdpp + gdtp*l,-1) );
+}
+
+FTYPE compute_Agrav( FTYPE om, FTYPE gdtt, FTYPE gdtp, FTYPE gdpp ){
+   return (sqrt(fabs(1./ ( gdtt + 2*om*gdtp + pow(om,2)*gdpp ) )));
+}
+
+FTYPE rmidlam( FTYPE x, FTYPE *parms ) {
+   FTYPE lamsq, gdtt, gdtp, gdpp, a, lambreak1, lambreak2, xi;
+   FTYPE lam_x, ans;
+
+   lamsq     = parms[0];   // square of target lambda
+   gdtt      = parms[1];
+   gdtp      = parms[2];
+   gdpp      = parms[3];
+   a         = parms[4];
+   lambreak1 = parms[5];
+   lambreak2 = parms[6];
+   xi        = parms[7];
+
+   compute_gd(x, M_PI_2, a, &gdtt, &gdtp, &gdpp);
+   lam_x = lamBL(x, gdtt, gdtp, gdpp, a, lambreak1, lambreak2, xi);   // lambda at current value of x
+
+   ans = lamsq - lam_x*lam_x;
+
+   return(ans);
+}
+
+FTYPE limotorus_findrml(FTYPE lam, FTYPE gdtt, FTYPE gdtp, FTYPE gdpp, FTYPE a, FTYPE lambreak1, FTYPE lambreak2, FTYPE xi) {
+   FTYPE parms[8];
+   FTYPE rml;
+
+   //store params in an array before function call
+   parms[0] = lam*lam;
+   parms[1] = gdtt;
+   parms[2] = gdtp;
+   parms[3] = gdpp;
+   parms[4] = a;
+   parms[5] = lambreak1;
+   parms[6] = lambreak2;
+   parms[7] = xi;
+
+   //solve for lin using bisection, specify large enough root search range, (1e-3, 1e3)
+   //demand accuracy 5x machine prec.
+   //in non-rel limit rml = r , use 10x that as the upper limit:
+   rml = rtbis( &rmidlam, parms, 6, 10*parms[0], 5.*DBL_EPSILON );
+
+   return( rml );
+}
+
+struct lnf_int_params {
+   FTYPE a, lambreak1, lambreak2, xi;
+};
+
+// Integrand in expression for lnf. This function is called
+// by the GSL quadrature routine.
+FTYPE lnf_integrand(FTYPE r, void *params) {
+   struct lnf_int_params *pars = (struct lnf_int_params *) params;
+   FTYPE a = pars->a, lambreak1 = pars->lambreak1, lambreak2 = pars->lambreak2, xi = pars->xi;
+   FTYPE r2, r3, a2, lam, lam2, lamroot;
+   FTYPE gdtt, gdtp, gdpp;
+   FTYPE l, om, dl_dr, dom_dr, integrand;
+   FTYPE term1, term2, term3, term4, D, E;
+   FTYPE oneplusx, dx_dlam, dlK_dlam, om_numerator, om_denominator;
+
+   // all values below are midplane values
+   r2 = pow(r,2);
+   r3 = pow(r,3);
+   a2 = a*a;
+   compute_gd(r, M_PI_2, a, &gdtt, &gdtp, &gdpp);
+   lam = lamBL(r, gdtt, gdtp, gdpp, a, lambreak1, lambreak2, xi);
+   lam2 = lam*lam;
+
+   l = l3d(lam, a, lambreak1, lambreak2, xi);
+
+   term1 = (r-2) * l;
+   term2 = 2*a;
+   term3 = r3 + a2*(r+2);
+   term4 = term2 * l;
+
+   om_numerator = term1 + term2;
+   om_denominator = term3 - term4;
+   om = om_numerator / om_denominator;
+
+   // derivatives
+   if (lam <= lambreak1 || lam >= lambreak2) {
+      dl_dr = 0;
+   } else {
+      oneplusx = 1 - 2/lam + a*pow(lam,-1.5);
+      dx_dlam = 2*pow(lam,-2) - 1.5*a*pow(lam,-2.5);
+      lamroot = sqrt(lam);
+      dlK_dlam = (oneplusx*0.5*pow(lamroot,-1) + (a-lamroot)*dx_dlam) / pow(oneplusx,2);
+      D = term3 - 2*term4 - lam2 * (r-2);
+      E = l * (3*r2 + a2 - lam2);
+      dl_dr = E / ( 2*lam*om_numerator / (xi*dlK_dlam) - D);
+   }
+
+   dom_dr = ( om_denominator * (l + (r-2)*dl_dr) - om_numerator * (3*r2+a2+2*a*dl_dr) )
+            / pow(om_denominator, 2);
+
+   integrand = -l/(1-om*l) * dom_dr;
+
+   return( integrand );
+}
+
+int init_dsandvels_limotorus(FTYPE r, FTYPE th, FTYPE a, FTYPE *rhoout, FTYPE *uuout, FTYPE *ell)
+{
+   FTYPE kappa;
+   FTYPE hh, eps;
+   FTYPE rho, u, ur, uh, up;
+   int pl;
+   FTYPE R, rbreak1, rbreak2, xi;
+   FTYPE lambreak1, lambreak2;
+   FTYPE lam, rml, l, om, Agrav, f3d;
+   FTYPE lamin, lin, omin, Agravin, f3din;
+   FTYPE lnf3d, lnferr;
+   FTYPE gdtt, gdtp, gdpp;
+
+   // GSL quadrature stuff
+   gsl_integration_workspace *intwork;
+   int worksize = 1000;   // workspace size
+   gsl_function integrand;
+   int intstatus;
+   
+   R = r*sin(th);
+   if (R < LT_RIN) {*rhoout = 0.; return(0);}
+
+   ///
+   /// Parameters
+   ///
+
+   kappa = LT_KAPPA;   // AKMARK: entropy constant that appears in EOS
+   xi = LT_XI;   // AKMARK: omega is set to this fraction of Keplerian omega
+   rbreak1 = LT_R1; // AKMARK: locations of breaks in torus angular momentum profile
+   rbreak2 = LT_R2;
+
+   //a=-0.7
+   //rbreak1 = 44.481;
+   //rbreak2 = 1000.;
+
+
+   ///
+   /// Computations at break radii
+   ///
+
+   // AKMARK: lamBL can calculate lambda at an arbitrary location, but it needs lambreak1,2;
+   // how to calculate lambreak1,2 themselves?
+   // Solution: note that lambreak1,2 are only needed by "l3d" in order to determine which region of the torus we are in.
+   // But lambreak1,2 can be considered to be in region 2, so just need a way to make "l3d" believe that we are in region 2.
+   // This can be done by setting lambreak1,2 to very small and very large values respectively.
+   compute_gd(rbreak1, M_PI_2, a, &gdtt, &gdtp, &gdpp);
+   lambreak1 = lamBL(rbreak1, gdtt, gdtp, gdpp, a, 0, 200000, xi);
+   compute_gd(rbreak2, M_PI_2, a, &gdtt, &gdtp, &gdpp);
+   lambreak2 = lamBL(rbreak2, gdtt, gdtp, gdpp, a, lambreak1, 200000, xi);
+
+   ///
+   /// Computations at torus inner edge
+   ///
+
+   compute_gd(LT_RIN, M_PI_2, a, &gdtt, &gdtp, &gdpp);
+   lamin = lamBL(LT_RIN, gdtt, gdtp, gdpp, a, lambreak1, lambreak2, xi);
+   lin = l3d(lamin, a, lambreak1, lambreak2, xi);
+   omin = omega3d(lin, gdtt, gdtp, gdpp);
+   Agravin = compute_Agrav(omin, gdtt, gdtp, gdpp);
+   f3din = 1.;   // the way f3d is defined, it equals 1 at r=rin
+
+   ///
+   /// Computations at current point: r, th
+   ///
+
+   compute_gd(r, th, a, &gdtt, &gdtp, &gdpp);
+   lam = lamBL(R, gdtt, gdtp, gdpp, a, lambreak1, lambreak2, xi);
+   l = l3d(lam, a, lambreak1, lambreak2, xi);
+   om = omega3d(l, gdtt, gdtp, gdpp);
+   Agrav = compute_Agrav(om, gdtt, gdtp, gdpp);
+
+   rml = limotorus_findrml( lam, gdtt, gdtp, gdpp, a, lambreak1, lambreak2, xi );
+
+   // f3d requires numerical evaluation of an integral
+
+   // First, set up things that GSL quadrature routine requires
+   struct lnf_int_params pars = {a, lambreak1, lambreak2, xi};
+   integrand.function = &lnf_integrand;   // integrand: function
+   integrand.params = &pars;   // other parameters to pass to integrand (besides integration variable)
+   intwork = gsl_integration_workspace_alloc(worksize);   // integration workspace
+
+   // then perform integration to obtain ln(f3d) ...
+   intstatus = gsl_integration_qags(&integrand, LT_RIN, rml, 0., 1.e-8, worksize, intwork, &lnf3d, &lnferr);
+   gsl_integration_workspace_free( intwork );
+   if (intstatus != GSL_SUCCESS) {
+      printf("GSL integration failed during limotorus setup at r=%21.15g, th=%21.15g; setting density to zero at this point", r, th);
+      lnf3d = GSL_NEGINF;   // cause density to be 0 at the current point
+   }
+
+   // ... and finally, f3d
+   f3d = exp(-lnf3d);
+
+   hh = f3din*Agrav / (f3d*Agravin);
+   eps = (-1 + hh)*pow(LT_GAMMA,-1);
+
+   if (eps < 0) rho = 0; else rho = pow((-1 + LT_GAMMA)*eps*pow(kappa,-1),pow(-1 + LT_GAMMA,-1));
+
+   *rhoout = rho;
+   *ell=l;
+   *uuout = kappa * pow(rho, LT_GAMMA) / (LT_GAMMA - 1.);
+
+   return(0);
+
+}
+
+
+
+int main() {
+
+   int nr = 60, nth = 36;
+   FTYPE Rin = 2., Rout = 500.;   // a = 0
+   //FTYPE Rin = 1.364, Rout = 2000.;   // a = 0.9
+   FTYPE th1 = 0., th2 = M_PI_2;
+   FTYPE dr = (Rout - Rin) / nr, dth = (th2 - th1) / nth;
+   FTYPE factor;
+   FTYPE r, th, rho, uu, ell;
+   int i, j;
+   FILE *outfile, *radfile;
+
+//    FTYPE R0 = 1.05, startx1 = -0.6736, dx1 = 0.0372985;
+
+   FTYPE a = LT_BHSPIN;
+   factor = log(Rout/Rin);
+
+   outfile = fopen("slice.dat", "w");
+   radfile = fopen("radslice.dat","w");
+
+   FTYPE Sigma;
+
+   for (i=0; i<=nr; i++) 
+     {
+       r = Rin*exp((FTYPE)i/nr*factor);
+       printf("r: %.2f\n",r);
+
+       Sigma=0.;
+       for (j=0; j<=nth; j++) 
+	 {
+	   th = th1 + j*dth;
+	   if(th>M_PI) th-=0.01;
+
+	   init_dsandvels_limotorus(r, th, a, &rho, &uu, &ell);
+
+	   Sigma+=rho*r*dth;
+
+	   fprintf(outfile, "%g\t%g\t%g\t%g\t%g\n", r*sin(th), r*cos(th), rho, uu,ell);
+	 }
+      fprintf(outfile, "\n");
+      fprintf(radfile,"%g %g %g\n",r,Sigma,ell);
+
+     }
+
+   fclose(outfile);
+   fclose(radfile);
+
+   return(0);
+
+}
+
+//    intstatus = gsl_integration_qag(&integrand, lamin, lam, 1.e-7, 1.e-7, worksize, GSL_INTEG_GAUSS61, intwork, &lnf3d, &lnferr);
diff --git a/PROBLEMS/PUFFY/postinit.c b/PROBLEMS/PUFFY/postinit.c
new file mode 100644
index 0000000..e7c48e6
--- /dev/null
+++ b/PROBLEMS/PUFFY/postinit.c
@@ -0,0 +1,163 @@
+//scales magnetic pressure so MAXBETA = pmag/(pgas+prad)
+//used only intitially
+int ix,iy,iz;
+
+#ifdef MAGNFIELD
+if(PROCID==0) {printf("Renormalizing magnetic field... ");fflush(stdout);}
+ldouble maxbeta=0.;
+#pragma omp parallel for private(ix,iy,iz) schedule (dynamic)
+for(iz=0;iz<NZ;iz++)
+  {
+    for(iy=0;iy<NY;iy++)
+      {
+	for(ix=0;ix<NX;ix++)
+	  {
+	    /***********************************************/
+	    ldouble pp[NV],uu[NV];
+	    struct geometry geom;
+	    fill_geometry(ix,iy,iz,&geom);
+
+	    struct geometry geomBL;
+	    fill_geometry_arb(ix,iy,iz,&geomBL,KERRCOORDS);
+ 
+	    int iv;
+	    PLOOP(iv)
+	      pp[iv]=get_u(p,iv,ix,iy,iz);
+
+	    ldouble bcon[4],bcov[4],bsq,pmag;
+	    calc_bcon_prim(pp,bcon,&geom);
+	    indices_21(bcon,bcov,geom.gg); 
+	    bsq = dot(bcon,bcov);
+	    pmag = bsq/2.;
+	    
+	    ldouble pgas=GAMMAM1*pp[UU];
+	    ldouble ptot=pgas;
+
+#ifdef RADIATION
+
+		ldouble Rtt,Ehat,ucon[4],prad;
+		calc_ff_Rtt(pp,&Rtt,ucon,&geom);
+		Ehat=-Rtt; 
+		prad=Ehat/3.;
+		ptot+=prad;
+		
+#endif
+
+#ifdef BETANORMFULL
+		//normalizing wrt everywhere
+#pragma omp critical
+		if(pmag/ptot>maxbeta) 
+		  {
+		    maxbeta=pmag/ptot;
+		    //printf("%d %d > %e %e %e\n",ix,iy,pmag,ptot,maxbeta);
+		  }
+		    
+#else //normalizing wrt to the equatorial plane
+		if(geom.iy==NY/2)
+		  {
+#pragma omp critical
+		    if(pmag/ptot>maxbeta) 
+		      maxbeta=pmag/ptot;
+		  }
+#endif
+	  }
+      }
+  }
+
+//choosing maximal maxbeta from the whole domain
+#ifdef MPI
+  MPI_Barrier(MPI_COMM_WORLD);
+  ldouble global_maxbeta;
+  
+  MPI_Allreduce(&maxbeta, &global_maxbeta, 1, MPI_DOUBLE, MPI_MAX,
+                MPI_COMM_WORLD);
+
+  maxbeta=global_maxbeta;
+#endif
+
+ldouble fac=sqrt(MAXBETA/maxbeta);
+
+//manual normalization - verify beta!
+#ifdef BETANORMFACTOR
+fac=BETANORMFACTOR;
+#endif
+printf("rescaling magn.fields by %e (%e)\n",fac,maxbeta);
+
+#pragma omp parallel for private(ix,iy,iz) schedule (dynamic)
+for(iz=0;iz<NZ;iz++)
+  {
+    for(iy=0;iy<NY;iy++)
+      {
+	for(ix=0;ix<NX;ix++)
+	  {
+	    /***********************************************/
+	    ldouble pp[NV],uu[NV];
+	    struct geometry geom;
+	    fill_geometry(ix,iy,iz,&geom);
+	    int iv;
+
+	    PLOOP(iv)
+	      pp[iv]=get_u(p,iv,ix,iy,iz);
+	    
+	    pp[B1]*=fac;
+	    pp[B2]*=fac;
+	    pp[B3]*=fac;
+
+	    p2u(pp,uu,&geom);
+
+	    PLOOP(iv)
+	    {
+	      set_u(p,iv,ix,iy,iz,pp[iv]);
+	      set_u(u,iv,ix,iy,iz,uu[iv]);
+	    }
+	  }
+      }
+  }
+
+/*
+printf("adjusting gas pressure...\n",fac);
+
+#pragma omp parallel for private(ix,iy,iz) schedule (dynamic)
+for(iz=0;iz<NZ;iz++)
+  {
+    for(iy=0;iy<NY;iy++)
+      {
+	for(ix=0;ix<NX;ix++)
+	  {
+
+	    ldouble pp[NV],uu[NV];
+	    struct geometry geom;
+	    fill_geometry(ix,iy,iz,&geom);
+	    int iv;
+
+	    PLOOP(iv)
+	      pp[iv]=get_u(p,iv,ix,iy,iz);
+	    
+	    ldouble pgas=GAMMAM1*pp[UU];
+
+	    ldouble bcon[4],bcov[4],bsq,pmag;
+	    calc_bcon_prim(pp,bcon,&geom);
+	    indices_21(bcon,bcov,geom.gg); 
+	    bsq = dot(bcon,bcov);
+	    pmag = bsq/2.;
+
+	    pgas-=pmag;
+
+	    pp[UU]=pgas;
+
+	    p2u(pp,uu,&geom);
+
+	    PLOOP(iv)
+	    {
+	      set_u(p,iv,ix,iy,iz,pp[iv]);
+	      set_u(u,iv,ix,iy,iz,uu[iv]);
+	    }
+
+	  }
+      }
+  }
+*/
+
+if(PROCID==0) {printf("done!\n");fflush(stdout);}
+#endif
+	    
diff --git a/PROBLEMS/PUFFY/prepinit.c b/PROBLEMS/PUFFY/prepinit.c
new file mode 100644
index 0000000..96fa681
--- /dev/null
+++ b/PROBLEMS/PUFFY/prepinit.c
@@ -0,0 +1,161 @@
+int init_dsandvels_limotorus(FTYPE r, FTYPE th, FTYPE a, FTYPE *rhoout, FTYPE *uuout, FTYPE *ell);
+
+//precalculates the initial state
+//done here to be able to refer here even after restart
+//B1-B3 in ptemp1 contains only the initial vector potential!
+//ENTR positive if inside torus, negative if outside
+
+if(PROCID==0) {fflush(stdout);printf("Precalculating torus and the atmosphere...\n");}
+int pr_iix,pr_iiy,pr_iiz;
+#pragma omp parallel for private(pr_iix,pr_iiy,pr_iiz) schedule (dynamic)
+for(pr_iiz=0;pr_iiz<NZ;pr_iiz++){
+	for(pr_iiy=0;pr_iiy<NY;pr_iiy++){
+		for(pr_iix=0;pr_iix<NX;pr_iix++){
+
+	    	ldouble rho,mx,my,mz,m,E,uint,pgas,Fx,Fy,Fz,pLTE,ell;  
+	    	ldouble uu[NV], pp[NV],ppback[NV],T,uintorg;
+	    	ldouble Vphi,Vr;
+	    	ldouble D,W,eps,uT,uphi,uPhi;
+
+	    	// metric in code coords
+	    	struct geometry geom;
+	    	fill_geometry(pr_iix,pr_iiy,pr_iiz,&geom);
+
+			// BL metric
+	    	struct geometry geomBL;
+	    	fill_geometry_arb(pr_iix,pr_iiy,pr_iiz,&geomBL,KERRCOORDS);
+
+	    	ldouble r=geomBL.xx;
+	    	ldouble th=geomBL.yy;
+
+	    	init_dsandvels_limotorus(r, th, BHSPIN, &rho, &uint, &ell);
+	    	uintorg=uint;
+
+	    	if(rho<0.){ //outside donut
+	      		set_hdatmosphere(pp,geom.xxvec,geom.gg,geom.GG,0);
+				#ifdef RADIATION
+					set_radatmosphere(pp,geom.xxvec,geom.gg,geom.GG,0);
+					#ifdef EVOLVEPHOTONNUMBER
+						pp[NF]=1./2.70118/K_BOLTZ * pp[EE] / ATMTRADINIT;
+					#endif
+				#endif
+			
+				pp[ENTR]=-1.;
+	    		}
+	    
+			else{ //inside donut
+	      
+				pp[ENTR]=1.;
+				
+				// DEBORA - shouldn't this be in BL?
+				set_hdatmosphere(ppback,geom.xxvec,geom.gg,geom.GG,0);
+				#ifdef RADIATION
+					set_radatmosphere(ppback,geom.xxvec,geom.gg,geom.GG,0);
+				#endif
+
+				uint=LT_KAPPA * pow(rho, LT_GAMMA) / (LT_GAMMA - 1.);
+				// DEBORA - shouldn't we use LT_GAMMA ?? 
+				pgas = GAMMAM1 * uint;
+				ell*=-1.;
+
+				ldouble ult,ulph,ucov[4],ucon[4];
+				ulph = sqrt(-1./(geomBL.GG[0][0]/ell/ell + 2./ell*geomBL.GG[0][3] + geomBL.GG[3][3]));
+				ult = ulph / ell;
+
+				ucov[0]=ult;
+				ucov[1]=0.;
+				ucov[2]=0.;
+				ucov[3]=ulph;
+    
+				indices_12(ucov,ucon,geomBL.GG);
+
+				conv_vels_ut(ucon,ucon,VEL4,VELPRIM,geomBL.gg,geomBL.GG);
+				
+				pp[0]=my_max(rho,ppback[0]); 
+				pp[1]=my_max(uint,ppback[1]);
+				pp[2]=ucon[1]; 
+				pp[3]=ucon[2];
+				pp[4]=ucon[3];
+
+				#ifdef MAGNFIELD//setting them zero not to break the following coordinate transformation
+						pp[B1]=pp[B2]=pp[B3]=0.; 
+				#endif
+
+				#ifdef RADIATION
+					//distributing pressure
+					ldouble P,aaa,bbb;
+					// DEBORA - again, LT_GAMMA ?? 
+					P=GAMMAM1*uint;
+					//solving for T satisfying P=pgas+prad=bbb T + aaa T^4
+					aaa=4.*SIGMA_RAD/3.;
+					bbb=K_BOLTZ*rho/MU_GAS/M_PROTON;
+					ldouble naw1=cbrt(9*aaa*Power(bbb,2) - Sqrt(3)*Sqrt(27*Power(aaa,2)*Power(bbb,4) + 256*Power(aaa,3)*Power(P,3)));
+					ldouble T4=-Sqrt((-4*Power(0.6666666666666666,0.3333333333333333)*P)/naw1 + naw1/(Power(2,0.3333333333333333)*Power(3,0.6666666666666666)*aaa))/2. + Sqrt((4*Power(0.6666666666666666,0.3333333333333333)*P)/naw1 - naw1/(Power(2,0.3333333333333333)*Power(3,0.6666666666666666)*aaa) + (2*bbb)/(aaa*Sqrt((-4*Power(0.6666666666666666,0.3333333333333333)*P)/naw1 + naw1/(Power(2,0.3333333333333333)*Power(3,0.6666666666666666)*aaa))))/2.;
+
+					E=calc_LTE_EfromT(T4);
+					Fx=Fy=Fz=0.;
+					uint=calc_PEQ_ufromTrho(T4,rho,pr_iix,pr_iiy,pr_iiz);
+
+					pp[UU]=my_max(uint,ppback[1]);
+					pp[EE0]=my_max(E,ppback[EE0]);
+
+
+
+					pp[FX0]=Fx;
+					pp[FY0]=Fy;
+					pp[FZ0]=Fz;
+
+					//transforming from BL lab radiative primitives to code non-ortonormal primitives
+					prad_ff2lab(pp,pp,&geomBL);
+
+				#endif
+
+				#ifdef EVOLVEPHOTONNUMBER
+					pp[NF0]=calc_NFfromE(pp[EE0]);
+				#endif
+
+				//transforming primitives from BL to MYCOORDS
+				trans_pall_coco(pp, pp, KERRCOORDS, MYCOORDS,geomBL.xxvec,&geomBL,&geom);
+    
+				#ifdef MAGNFIELD 
+					//MYCOORDS vector potential to calculate B's
+					ldouble Acov[4];
+					Acov[0]=Acov[1]=Acov[2]=0.;
+
+					Acov[3]=my_max(pow(pp[RHO]*geomBL.xx/4.e-22,2.)-0.02,0.)*sqrt(1.e-23);
+
+					#ifdef QUADLOOPS
+						Acov[3]*=sin((M_PI/2.-geomBL.yy)/0.1);
+					#endif
+
+					#ifdef MULTIPLELOOPS
+						Acov[3]*=sin(geomBL.xx/3.);
+					#endif
+
+					#ifdef SINGLELOOP //standard single poloidal loop
+						Acov[3]=my_max(pow(pp[RHO]*geomBL.xx*geomBL.xx/4.e-20,2.)-0.02,0.)*sqrt(1.e-23)*pow(sin(fabs(geomBL.yy)),4.);
+					#endif
+
+					pp[B1]=Acov[1];
+					pp[B2]=Acov[2];
+					pp[B3]=Acov[3];
+				#endif
+
+	    	} //end of inside torus
+
+	    
+	   		//to conserved
+	    	p2u(pp,uu,&geom);
+
+	    	int iv;
+
+	    	for(iv=0;iv<NV;iv++){
+				set_u(ptemp1,iv,pr_iix,pr_iiy,pr_iiz,pp[iv]);
+	    	}
+
+		}
+    }
+}
+
+
+
diff --git a/PROBLEMS/PUFFY/tools.c b/PROBLEMS/PUFFY/tools.c
new file mode 100644
index 0000000..4f8da7d
--- /dev/null
+++ b/PROBLEMS/PUFFY/tools.c
@@ -0,0 +1,350 @@
+void compute_gd( FTYPE r, FTYPE th, FTYPE a, FTYPE *gdtt, FTYPE *gdtp, FTYPE *gdpp ) {
+   FTYPE Sigma, tmp;
+
+   Sigma = (pow(a*cos(th),2) + pow(r,2));
+   tmp = 2*r*pow(sin(th),2)/Sigma;
+
+   //metric (expressions taken from limotorus4.nb):
+   *gdtt = -1 + 2*r/Sigma;
+   *gdtp = -a*tmp;
+   *gdpp = (pow(r,2) + pow(a,2)*(1+tmp)) * pow(sin(th),2);
+}
+
+// Keplerian, equatorial angular momentum density: l = u_phi / u_t
+FTYPE lK(FTYPE r, FTYPE a) {
+   FTYPE curlyF, curlyG;
+
+   curlyF = 1 - 2*a/pow(r, 1.5) + pow(a/r, 2);
+   curlyG = 1 - 2/r + a/pow(r, 1.5);
+   return( pow(r,0.5) * curlyF / curlyG );
+}
+
+FTYPE l3d(FTYPE lam, FTYPE a, FTYPE lambreak1, FTYPE lambreak2, FTYPE xi) {
+   return ( xi * lK( lam<=lambreak1 ? lambreak1 : lam>=lambreak2 ? lambreak2 : lam , a) );
+}
+
+FTYPE rtbis(FTYPE (*func)(FTYPE,FTYPE*), FTYPE *parms, FTYPE x1, FTYPE x2, FTYPE xacc)
+//Taken from HARM:nrutil.c
+{
+   int j;
+   FTYPE dx,f,fmid,xmid,rtb;
+   f=(*func)(x1, parms);
+   fmid=(*func)(x2, parms);
+   if (f*fmid >= 0.0) {
+      printf( "f(%g)=%g f(%g)=%g\n", x1, f, x2, fmid );
+      printf("Root must be bracketed for bisection in rtbis\n");
+   }
+   rtb = (f < 0.0) ? (dx=x2-x1,x1) : (dx=x1-x2,x2); //Orient the search so that f>0 lies at x+dx.
+   for (j=1;j<=100;j++) {
+      fmid=(*func)(xmid=rtb+(dx *= 0.5),parms); //Bisection loop.
+      if (fmid <= 0.0) {
+         rtb=xmid;
+      }
+      if (fabs(dx) < xacc || fmid == 0.0) {
+         return rtb;
+      }
+   }
+   printf("Too many bisections in rtbis\n");
+   return 0.0; //Never get here.
+}
+
+FTYPE lamBL_func(FTYPE lam, FTYPE *parms) {
+   FTYPE gdtt, gdtp, gdpp, a, lambreak1, lambreak2, xi, l;
+
+   gdtt      = parms[0];
+   gdtp      = parms[1];
+   gdpp      = parms[2];
+   a         = parms[3];
+   lambreak1 = parms[4];
+   lambreak2 = parms[5];
+   xi        = parms[6];
+
+   l = l3d(lam, a, lambreak1, lambreak2, xi);
+
+   return ( lam*lam + l * (l*gdtp + gdpp) / (l*gdtt + gdtp) );
+}
+
+// von Zeipel cylinder radius (needs to be calculated iteratively for nonzero spin)
+FTYPE lamBL(FTYPE R, FTYPE gdtt, FTYPE gdtp, FTYPE gdpp, FTYPE a, FTYPE lambreak1, FTYPE lambreak2, FTYPE xi) {
+   // R = r*sin(th), used as initial guess for lamBL
+
+   FTYPE parms[7];
+
+   //store params in an array before function call
+   parms[0] = gdtt;
+   parms[1] = gdtp;
+   parms[2] = gdpp;
+   parms[3] = a;
+   parms[4] = lambreak1;
+   parms[5] = lambreak2;
+   parms[6] = xi;
+
+   //solve for lin using bisection, specify large enough root search range, (1e-3, 1e3)
+   //demand accuracy 5x machine prec.
+   //in non-rel limit rml = r , use 10x that as the upper limit:
+   return( rtbis( &lamBL_func, parms, R, 10*R, 5.*DBL_EPSILON ) );
+
+}
+
+FTYPE omega3d( FTYPE l, FTYPE gdtt, FTYPE gdtp, FTYPE gdpp ) {
+   return( -(gdtt*l + gdtp)*pow(gdpp + gdtp*l,-1) );
+}
+
+FTYPE compute_Agrav( FTYPE om, FTYPE gdtt, FTYPE gdtp, FTYPE gdpp ){
+   return (sqrt(fabs(1./ ( gdtt + 2*om*gdtp + pow(om,2)*gdpp ) )));
+}
+
+FTYPE rmidlam( FTYPE x, FTYPE *parms ) {
+   FTYPE lamsq, gdtt, gdtp, gdpp, a, lambreak1, lambreak2, xi;
+   FTYPE lam_x, ans;
+
+   lamsq     = parms[0];   // square of target lambda
+   gdtt      = parms[1];
+   gdtp      = parms[2];
+   gdpp      = parms[3];
+   a         = parms[4];
+   lambreak1 = parms[5];
+   lambreak2 = parms[6];
+   xi        = parms[7];
+
+   compute_gd(x, M_PI_2, a, &gdtt, &gdtp, &gdpp);
+   lam_x = lamBL(x, gdtt, gdtp, gdpp, a, lambreak1, lambreak2, xi);   // lambda at current value of x
+
+   ans = lamsq - lam_x*lam_x;
+
+   return(ans);
+}
+
+FTYPE limotorus_findrml(FTYPE lam, FTYPE gdtt, FTYPE gdtp, FTYPE gdpp, FTYPE a, FTYPE lambreak1, FTYPE lambreak2, FTYPE xi) {
+   FTYPE parms[8];
+   FTYPE rml;
+
+   //store params in an array before function call
+   parms[0] = lam*lam;
+   parms[1] = gdtt;
+   parms[2] = gdtp;
+   parms[3] = gdpp;
+   parms[4] = a;
+   parms[5] = lambreak1;
+   parms[6] = lambreak2;
+   parms[7] = xi;
+
+   //solve for lin using bisection, specify large enough root search range, (1e-3, 1e3)
+   //demand accuracy 5x machine prec.
+   //in non-rel limit rml = r , use 10x that as the upper limit:
+   rml = rtbis( &rmidlam, parms, 6, 10*parms[0], 5.*DBL_EPSILON );
+
+   return( rml );
+}
+
+struct lnf_int_params {
+   FTYPE a, lambreak1, lambreak2, xi;
+};
+
+// Integrand in expression for lnf. This function is called
+// by the GSL quadrature routine.
+FTYPE lnf_integrand(FTYPE r, void *params) {
+   struct lnf_int_params *pars = (struct lnf_int_params *) params;
+   FTYPE a = pars->a, lambreak1 = pars->lambreak1, lambreak2 = pars->lambreak2, xi = pars->xi;
+   FTYPE r2, r3, a2, lam, lam2, lamroot;
+   FTYPE gdtt, gdtp, gdpp;
+   FTYPE l, om, dl_dr, dom_dr, integrand;
+   FTYPE term1, term2, term3, term4, D, E;
+   FTYPE oneplusx, dx_dlam, dlK_dlam, om_numerator, om_denominator;
+
+   // all values below are midplane values
+   r2 = pow(r,2);
+   r3 = pow(r,3);
+   a2 = a*a;
+   compute_gd(r, M_PI_2, a, &gdtt, &gdtp, &gdpp);
+   lam = lamBL(r, gdtt, gdtp, gdpp, a, lambreak1, lambreak2, xi);
+   lam2 = lam*lam;
+
+   l = l3d(lam, a, lambreak1, lambreak2, xi);
+
+   term1 = (r-2) * l;
+   term2 = 2*a;
+   term3 = r3 + a2*(r+2);
+   term4 = term2 * l;
+
+   om_numerator = term1 + term2;
+   om_denominator = term3 - term4;
+   om = om_numerator / om_denominator;
+
+   // derivatives
+   if (lam <= lambreak1 || lam >= lambreak2) {
+      dl_dr = 0;
+   } else {
+      oneplusx = 1 - 2/lam + a*pow(lam,-1.5);
+      dx_dlam = 2*pow(lam,-2) - 1.5*a*pow(lam,-2.5);
+      lamroot = sqrt(lam);
+      dlK_dlam = (oneplusx*0.5*pow(lamroot,-1) + (a-lamroot)*dx_dlam) / pow(oneplusx,2);
+      D = term3 - 2*term4 - lam2 * (r-2);
+      E = l * (3*r2 + a2 - lam2);
+      dl_dr = E / ( 2*lam*om_numerator / (xi*dlK_dlam) - D);
+   }
+
+   dom_dr = ( om_denominator * (l + (r-2)*dl_dr) - om_numerator * (3*r2+a2+2*a*dl_dr) )
+            / pow(om_denominator, 2);
+
+   integrand = -l/(1-om*l) * dom_dr;
+
+   return( integrand );
+}
+
+int init_dsandvels_limotorus(FTYPE r, FTYPE th, FTYPE a, FTYPE *rhoout, FTYPE *uuout, FTYPE *ell)
+{
+   FTYPE kappa;
+   FTYPE hh, eps;
+   FTYPE rho, u, ur, uh, up;
+   int pl;
+   FTYPE R, rbreak1, rbreak2, xi;
+   FTYPE lambreak1, lambreak2;
+   FTYPE lam, rml, l, om, Agrav, f3d;
+   FTYPE lamin, lin, omin, Agravin, f3din;
+   FTYPE lnf3d, lnferr;
+   FTYPE gdtt, gdtp, gdpp;
+
+   // GSL quadrature stuff
+   gsl_integration_workspace *intwork;
+   int worksize = 1000;   // workspace size
+   gsl_function integrand;
+   int intstatus;
+   R = r*sin(th);
+   if (R < LT_RIN) {*rhoout = -1.; return(0);}
+
+   ///
+   /// Parameters
+   ///
+
+   kappa = LT_KAPPA;   // AKMARK: entropy constant that appears in EOS
+   xi = LT_XI;   // AKMARK: omega is set to this fraction of Keplerian omega
+   rbreak1 = LT_R1; // AKMARK: locations of breaks in torus angular momentum profile
+   rbreak2 = LT_R2;
+
+   //a=-0.7
+   //rbreak1 = 44.481;
+   //rbreak2 = 1000.;
+
+
+   ///
+   /// Computations at break radii
+   ///
+
+   // AKMARK: lamBL can calculate lambda at an arbitrary location, but it needs lambreak1,2;
+   // how to calculate lambreak1,2 themselves?
+   // Solution: note that lambreak1,2 are only needed by "l3d" in order to determine which region of the torus we are in.
+   // But lambreak1,2 can be considered to be in region 2, so just need a way to make "l3d" believe that we are in region 2.
+   // This can be done by setting lambreak1,2 to very small and very large values respectively.
+   compute_gd(rbreak1, M_PI_2, a, &gdtt, &gdtp, &gdpp);
+   lambreak1 = lamBL(rbreak1, gdtt, gdtp, gdpp, a, 0, 200000, xi);
+   compute_gd(rbreak2, M_PI_2, a, &gdtt, &gdtp, &gdpp);
+   lambreak2 = lamBL(rbreak2, gdtt, gdtp, gdpp, a, lambreak1, 200000, xi);
+
+   ///
+   /// Computations at torus inner edge
+   ///
+
+   compute_gd(LT_RIN, M_PI_2, a, &gdtt, &gdtp, &gdpp);
+   lamin = lamBL(LT_RIN, gdtt, gdtp, gdpp, a, lambreak1, lambreak2, xi);
+   lin = l3d(lamin, a, lambreak1, lambreak2, xi);
+   omin = omega3d(lin, gdtt, gdtp, gdpp);
+   Agravin = compute_Agrav(omin, gdtt, gdtp, gdpp);
+   f3din = 1.;   // the way f3d is defined, it equals 1 at r=rin
+
+   ///
+   /// Computations at current point: r, th
+   ///
+
+   compute_gd(r, th, a, &gdtt, &gdtp, &gdpp);
+   lam = lamBL(R, gdtt, gdtp, gdpp, a, lambreak1, lambreak2, xi);
+   l = l3d(lam, a, lambreak1, lambreak2, xi);
+   om = omega3d(l, gdtt, gdtp, gdpp);
+   Agrav = compute_Agrav(om, gdtt, gdtp, gdpp);
+
+   rml = limotorus_findrml( lam, gdtt, gdtp, gdpp, a, lambreak1, lambreak2, xi );
+
+   // f3d requires numerical evaluation of an integral
+
+   // First, set up things that GSL quadrature routine requires
+   struct lnf_int_params pars = {a, lambreak1, lambreak2, xi};
+   integrand.function = &lnf_integrand;   // integrand: function
+   integrand.params = &pars;   // other parameters to pass to integrand (besides integration variable)
+   intwork = gsl_integration_workspace_alloc(worksize);   // integration workspace
+
+   // then perform integration to obtain ln(f3d) ...
+   intstatus = gsl_integration_qags(&integrand, LT_RIN, rml, 0., 1.e-8, worksize, intwork, &lnf3d, &lnferr);
+   gsl_integration_workspace_free( intwork );
+   if (intstatus != GSL_SUCCESS) {
+      printf("GSL integration failed during limotorus setup at r=%21.15g, th=%21.15g; setting density to zero at this point", r, th);
+      lnf3d = GSL_NEGINF;   // cause density to be 0 at the current point
+   }
+
+   // ... and finally, f3d
+   f3d = exp(-lnf3d);
+
+   hh = f3din*Agrav / (f3d*Agravin);
+   eps = (-1 + hh)*pow(LT_GAMMA,-1);
+
+   if (eps < 0) rho = -1; else rho = pow((-1 + LT_GAMMA)*eps*pow(kappa,-1),pow(-1 + LT_GAMMA,-1));
+
+   *rhoout = rho;
+   *ell=l;
+   *uuout = kappa * pow(rho, LT_GAMMA) / (LT_GAMMA - 1.);
+
+   //my_err("wfsg");
+   //printf("%e %e %e\n",rho,kappa,*uuout);
+   //exit(0);   getch();
+
+   return(0);
+
+}
+
+
+/*
+int main() {
+
+   int nr = 48, nth = 16;
+   FTYPE Rin = 1.6, Rout = 200.;   // a = 0
+   //FTYPE Rin = 1.364, Rout = 2000.;   // a = 0.9
+   FTYPE th1 = 0., th2 = M_PI_2;
+   FTYPE dr = (Rout - Rin) / nr, dth = (th2 - th1) / nth;
+   FTYPE factor;
+   FTYPE r, th, rho, uu, ell;
+   int i, j;
+   FILE*outfile;
+
+//    FTYPE R0 = 1.05, startx1 = -0.6736, dx1 = 0.0372985;
+
+   FTYPE a = 0.9;
+   factor = log(Rout/Rin);
+
+   outfile = fopen("slice.dat", "w");
+//    fprintf(outfile, "variables = x, z, rho\n");
+//    fprintf(outfile, "zone t=rho, i=%i, j=%i, F=POINT\n", nr+1, nth+1);
+   for (j=0; j<=nth; j++) {
+     th = th1 + j*dth;
+     if(th>M_PI) th-=0.01;
+//       th = 1.;
+      for (i=0; i<=nr; i++) {
+//          r = Rin + i*dr;
+         r = Rin*exp((FTYPE)i/nr*factor);
+//          r = R0 + exp(startx1 + (i+0.5)*dx1);
+//          r = 70.;
+         init_dsandvels_limotorus(r, th, a, &rho, &uu, &ell);
+//          fprintf(outfile, "%g\t%g\t%g\n", r*sin(th), r*cos(th), rho);
+         fprintf(outfile, "%g\t%g\t%g\t%g\t%g\n", r, th, rho, uu,ell);
+      }
+      fprintf(outfile, "\n");
+   }
+   fclose(outfile);
+
+//    th = M_PI_2*0.3;
+//    r = 350.;
+//    init_dsandvels_limotorus(r, th, a, &rho);
+//    printf("%f\n", rho);
+
+   return(0);
+
+}
+*/
diff --git a/ana.c b/ana.c
index 3585175..c81517e 100644
--- a/ana.c
+++ b/ana.c
@@ -9,6 +9,10 @@
 #define INCLUDE_WRITE_OUTPUT
 #endif
 
+#ifdef PRINT_FIXUPS_TO_SILO
+  printf("PRINT_FIXUPS_TO_SILO works only on the run!\n\n");
+#endif
+
 int 
 main(int argc, char **argv)
 {  
diff --git a/avg.c b/avg.c
index 9bc35f9..f0d6999 100644
--- a/avg.c
+++ b/avg.c
@@ -13,6 +13,12 @@
 int 
 main(int argc, char **argv)
 {  
+
+  #ifdef PRINT_FIXUPS_TO_SILO
+    printf("PRINT_FIXUPS_TO_SILO works only on the run!\n\n");
+  
+  #endif
+
   #ifdef MPI
   printf("avg works on shared memory only, do not use MPI, please\n");
   exit(-1);
diff --git a/mnemonics.h b/mnemonics.h
index 41205c7..8f96483 100644
--- a/mnemonics.h
+++ b/mnemonics.h
@@ -1,7 +1,6 @@
 /*! \file ko.h
  \brief mnemonics for various indices
 */
-#pragma once
 
 //rad or hydro
 #define RAD 1
@@ -167,6 +166,7 @@
 //frames
 #define ZAMOFRAME 0
 #define FFFRAME 1
+#define DRIFTFRAME 2
 
 //avg quantities
 #define AVGBSQ (NV+0)
diff --git a/problem.h b/problem.h
index bce2c54..f264277 100644
--- a/problem.h
+++ b/problem.h
@@ -1,6 +1,5 @@
 //KORAL - problem.h
 //choice of the problem plus some definitions
-#pragma once
 
 //available problems:
 
@@ -144,13 +143,27 @@
 //140 RADSURVEY -- radiative parameter survey
 //141 KEPINF -- INFDISK modified for injecting keplerian 
 //142 PARTIALTDE -- INFDISK modified for partial TDE binding energy distribution non-constant
+//143 KKDISC
+// 144 PUFFY
 
 //ANDREW -- I've gone through problems 100-133 and undefined PR_KAPPA where appropriate
 //If you want to use default calc_opacities_from_state, make sure PR_KAPPA is  undefined!
 //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 140
+#define PROBLEM 144
+
+#if(PROBLEM==144)
+
+#define PR_DEFINE "PROBLEMS/PUFFY/define.h"
+#define PR_BC "PROBLEMS/PUFFY/bc.c"
+#define PR_INIT "PROBLEMS/PUFFY/init.c"
+#define PR_PREPINIT "PROBLEMS/PUFFY/prepinit.c"
+#define PR_POSTINIT "PROBLEMS/PUFFY/postinit.c"
+#define PR_KAPPAES "PROBLEMS/PUFFY/kappaes.c"
+#define PR_TOOLS "PROBLEMS/PUFFY/tools.c"
+
+#endif
 
 #if(PROBLEM==142)
 
@@ -750,7 +763,7 @@
 #define PR_DEFINE "PROBLEMS/THINDISK/define.h"
 #define PR_BC "PROBLEMS/THINDISK/bc.c"
 #define PR_INIT "PROBLEMS/THINDISK/init.c"
-#define PR_KAPPA "PROBLEMS/THINDISK/kappa.c"
+//#define PR_KAPPA "PROBLEMS/THINDISK/kappa.c"
 #define PR_KAPPAES "PROBLEMS/THINDISK/kappaes.c"
 #define PR_DUMP "PROBLEMS/THINDISK/dump.c"
 #define PR_TOOLS "PROBLEMS/THINDISK/tools.c"
diff --git a/regrid.c b/regrid.c
old mode 100644
new mode 100755
index 8d0acc9..c94b6f8
--- a/regrid.c
+++ b/regrid.c
@@ -31,9 +31,9 @@ main(int argc, char **argv)
   omp_myinit();  
   #endif
 
-  if (argc!=18 && argc!=12)
+  if (argc!=19 && argc!=12)
     {
-      printf("usage: ./regrid NX1 NY1 NZ1 coord1.dat res0001.dat NX2 NY2 NZ2 coord2.dat res0002.dat NV [regridinx=1 regridiny=1 regridinz=1 periodicphi=1 pert=0.0 interpolate=0]\n");
+      printf("usage: ./regrid NX1 NY1 NZ1 coord1.dat res0001.dat NX2 NY2 NZ2 coord2.dat res0002.dat NV [regridinx=1 regridiny=1 regridinz=1 periodicphi=1 pert=0.0 interpolate=0 divBclean=1]\n");
       exit(1);
     }
 
@@ -64,6 +64,7 @@ main(int argc, char **argv)
   char file_coord1[100],file_coord2[100],file_res1[100],file_res2[100];
   int periodicphi,regridx,regridy,regridz,interpolate;
   double pert=0.0;
+  int divBclean = 1;
 
   //random number gen. initialization
   srand ( time(NULL) );
@@ -84,7 +85,7 @@ main(int argc, char **argv)
 
   //NV=atoi(argv[11]);
 
-  if(argc==18)
+  if(argc==19)
     {
       regridx=atoi(argv[12]);
       regridy=atoi(argv[13]);
@@ -94,6 +95,8 @@ main(int argc, char **argv)
       pert=atof(argv[16]);
 
       interpolate=atoi(argv[17]);
+
+      divBclean = atoi(argv[18]);
     }
   else
     {
@@ -103,6 +106,8 @@ main(int argc, char **argv)
       periodicphi=1;
       pert=0.0;
       interpolate=0;
+      divBclean = 1;
+      
     }
 
   printf("regrid - projects snap shot file onto a new grid\n");
@@ -115,6 +120,7 @@ main(int argc, char **argv)
   if(periodicphi) printf("periodic in phi\n");
   if(pert!=0.0) printf("pert=%f\n",pert);
   if(interpolate) printf("interpolating\n");
+  if(divBclean) printf("Divergence of B cleaning\n");
 
   //allocating memory
   double *x1,**y1,*z1,*x2,**y2,*z2;
@@ -723,22 +729,22 @@ main(int argc, char **argv)
 
            if(i1==0)// && i2==0)
            { 
-             printf("Filling xxvec...\n");
+             //printf("Filling xxvec...\n");
              get_xx(ix,iy,iz,xxvec);
-             printf("%f %f %f %f\n",xxvec[0],xxvec[1],xxvec[2],xxvec[3]);
+             //printf("%f %f %f %f\n",xxvec[0],xxvec[1],xxvec[2],xxvec[3]);
 
-             printf("Filling geomout\n");
+             //printf("Filling geomout\n");
 	     fill_geometry(ix,iy,iz,&geomout);
-             printf("Filling geomin\n");
+             //printf("Filling geomin\n");
 	     fill_geometry_arb(ix,iy,iz,&geomin,BLCOORDS);//MKS2COORDS);
-             printf("%f %f \n",geomin.gdet,geomout.gdet);
-             printf("%f %f %f \n",geomin.xx,geomin.yy,geomin.zz);
+             //printf("%f %f \n",geomin.gdet,geomout.gdet);
+             //printf("%f %f %f \n",geomin.xx,geomin.yy,geomin.zz);
            }
 
 	   if(pert>0.0) dataout[i1][i2][i3][VZ]*=(1.+((double)rand()/(double)RAND_MAX-0.5)*2.*pert);
 	 }
 
-   //printf("Interpolated/copied data. \n");
+   printf("Interpolated/copied data. \n");
 
    //writing to the new file
    FILE *fout = fopen(file_res2,"wb");
@@ -761,6 +767,7 @@ main(int argc, char **argv)
   //perform divergence cleaning (only in 2D) if simulation has B-field
   // (this is non-relativistic, may have to expand code to handle this properly
   #ifdef MAGNFIELD //Only use for MHD
+  if(divBclean == 1){
   for(ix=0;ix<NX2;ix++)
     for(iy=0;iy<NY2;iy++)
       for(iz=0;iz<NZ2;iz++)
@@ -839,6 +846,7 @@ main(int argc, char **argv)
           new_Eb += (Bxx*Bxx + Byy*Byy + Bzz*Bzz)*dxx*dyy*dzz*pick_gdet(ix,iy,iz);
 	}
    #endif
+  }
 
   //compute factor for rescaling B-field to maintain same total magnetic energy in domain
   alpha = 1.;
@@ -854,6 +862,7 @@ main(int argc, char **argv)
       for(iz=0;iz<NZ2;iz++)
 	{ 
           #ifdef MAGNFIELD
+          if(divBclean){
           //rescale magnetic field
           dataout[ix][iy][iz][B1] *= alpha;
           dataout[ix][iy][iz][B2] *= alpha;
@@ -875,13 +884,14 @@ main(int argc, char **argv)
           Bzz = dataout[ix][iy][iz][B3];
           newnew_Eb += (Bxx*Bxx + Byy*Byy + Bzz*Bzz)*dxx*dyy*dzz*pick_gdet(ix,iy,iz);
           #endif
-
+          }
 	  fwrite(dataout[ix][iy][iz],sizeof(double),NV,fout);
 	}
 
   #ifdef MAGNFIELD
-  //printf("old_Eb new_Eb newnew_Eb a a^2 : %e %e %e %e %e \n",old_Eb,new_Eb,newnew_Eb,alpha,alpha*alpha);
+  printf("old_Eb new_Eb newnew_Eb a a^2 : %e %e %e %e %e \n",old_Eb,new_Eb,newnew_Eb,alpha,alpha*alpha);
   #endif
+  
 
   fclose(fout);
 
diff --git a/silo.c b/silo.c
index e77763e..f99f500 100644
--- a/silo.c
+++ b/silo.c
@@ -198,7 +198,16 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
   ldouble *entrrad = (ldouble*)malloc(nx*ny*nz*sizeof(double));
   #endif //RADIATION
 
-    
+	
+	// Debora - print out where are the cells which needs fixups (any)
+	#ifdef PRINT_FIXUPS_TO_SILO
+		int *fixups_u2pmhd = (int*)malloc(nx*ny*nz*sizeof(int));
+		int *fixups_u2prad = (int*)malloc(nx*ny*nz*sizeof(int));
+		int *fixups_radimp = (int*)malloc(nx*ny*nz*sizeof(int));
+		int *fixups = (int*)malloc(nx*ny*nz*sizeof(int));
+	#endif
+
+
   //first fill coordinates on nodes
 #pragma omp parallel for private(ix,iy,iz,iv,imx,imy,imz,i,j,pp,uu,xxvec,xxveccar,xxvecsph,xx1,xx2) schedule (static)
 
@@ -369,6 +378,8 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
 	      ldouble th=xxvecsph[2];
 	      ldouble ph=xxvecsph[3];
 
+
+
 	      //gdet and coordinates of cells +- 1 in radius
 	      ldouble gdet1,gdet2,gdet,gdetu;
 	      gdet=geomout.gdet;
@@ -438,6 +449,9 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
 	      //magnetic fields
               #ifdef MAGNFIELD
 	      ldouble bcon[4],bcov[4];
+
+
+
 	      if(doingavg==0)
 		{
 		  calc_bcon_prim(pp,bcon,&geomout);
@@ -599,10 +613,38 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
 	      #ifdef PRINTVISCHEATINGTOSILO
 	      deltae[zonalindex]=calc_ViscousElectronHeatingFraction(&get_u(p,0,ix,iy,iz),&geomout);
 	      vischeat[zonalindex]=get_u_scalar(vischeating,ix,iy,iz);
-	      vischeatnege[zonalindex]=get_u_scalar(vischeatingnegebalance,ix,iy,iz);; 
-	      vischeatnegi[zonalindex]=get_u_scalar(vischeatingnegibalance,ix,iy,iz);; 
+	      vischeatnege[zonalindex]=get_u_scalar(vischeatingnegebalance,ix,iy,iz); 
+	      vischeatnegi[zonalindex]=get_u_scalar(vischeatingnegibalance,ix,iy,iz);
               dtauarr[zonalindex]=-1.;
 	      #endif
+
+		#ifdef PRINT_FIXUPS_TO_SILO
+			fixups_u2pmhd[zonalindex] = 0;
+			fixups_u2prad[zonalindex] = 0;
+			fixups_radimp[zonalindex] = 0;
+			fixups[zonalindex] = 0;
+			int nfixMHD = 0;
+			int nfixRAD = 0;
+			int nfixRADimp = 0;
+			
+			nfixMHD = get_cflag(HDFIXUPFLAG,ix,iy,iz);
+			if (nfixMHD!=0){
+				fixups_u2pmhd[zonalindex]+=nfixMHD;
+				fixups++;
+			}
+			nfixRAD = get_cflag(RADFIXUPFLAG,ix,iy,iz);
+			if (nfixRAD!=0){
+				fixups_u2prad[zonalindex]+=nfixRAD;
+				fixups++;
+			}  
+			nfixRADimp = get_cflag(RADIMPFIXUPFLAG,ix,iy,iz);
+			if (nfixRADimp!=0){
+				fixups_radimp[zonalindex]+=nfixRADimp;
+				fixups++;
+			}
+		#endif	
+	
+	
 		}
 	      
 	      else //using averaged data
@@ -730,6 +772,8 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix)
 		  deltae[zonalindex]=-1.;
 		  dtauarr[zonalindex]=-1.;
 #endif
+
+
 		} //doingavg
 
 
@@ -1783,7 +1827,18 @@ DBPutQuadvar1(file, "Gtff","mesh1", Gtff,
 		DB_DOUBLE, DB_ZONECENT, optList);
 
   #endif
-  #endif
+  #endif // MAGFIELD
+
+	#ifdef PRINT_FIXUPS_TO_SILO
+		DBPutQuadvar1(file, "fixups_u2pmhd","mesh1", fixups_u2pmhd,
+  		dimensions, ndim, NULL, 0, DB_INT, DB_ZONECENT, optList);
+		DBPutQuadvar1(file, "fixups_u2prad","mesh1", fixups_u2prad,
+  		dimensions, ndim, NULL, 0, DB_INT, DB_ZONECENT, optList);
+		DBPutQuadvar1(file, "fixups_radimp","mesh1", fixups_radimp,
+  		dimensions, ndim, NULL, 0, DB_INT, DB_ZONECENT, optList);
+		DBPutQuadvar1(file, "fixups","mesh1", fixups,
+  		dimensions, ndim, NULL, 0, DB_INT, DB_ZONECENT, optList);
+	#endif  		
 
 
   // Write vectors to file
@@ -2066,6 +2121,13 @@ DBPutQuadvar1(file, "Gtff","mesh1", Gtff,
   free(entrrad);
   #endif
 
+	#ifdef PRINT_FIXUPS_TO_SILO
+		free(fixups_u2pmhd);
+  		free(fixups_u2prad);
+  		free(fixups_radimp);
+		free(fixups);
+	#endif  
+
   return (0);
 }
 
diff --git a/u2p.c b/u2p.c
index 9df57f0..bf1ee57 100644
--- a/u2p.c
+++ b/u2p.c
@@ -453,37 +453,66 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg)
   ldouble ucond[4],ucovd[4];
   ldouble bcond[4],bcovd[4],bsq,magpre;
   ldouble etacon[4],etarel[4];
-  for(iv=1;iv<4;iv++)
+  ldouble f = 1.;
+  ldouble fuu = 1.;
+  int rhofloored = 0;
+  int uufloored = 0;
+
+  for(iv=1;iv<4;iv++){
     ucond[iv]=pp[1+iv];
+  }
+  
   calc_ucon_ucov_from_prims(pp, geom, ucond, ucovd);
   calc_bcon_bcov_bsq_from_4vel(pp, ucond, ucovd, geom, bcond, bcovd, &bsq);
   //magpre = 0.5 * bsq;
   magpre = bsq; // ANDREW: B2RHORATIOMAX and B2UURATIOMAX make more sense with bsq, not .5bsq
-  
-  calc_normalobs_ncon(GG, geom->alpha, etacon);
-  conv_vels_ut(etacon,etarel,VEL4,VELPRIM,gg,GG);
 
-  if(magpre>B2RHORATIOMAX*pp[RHO]) 
-  {
+
+  // check density vs bsq
+  if(magpre>B2RHORATIOMAX*pp[RHO]) {
     if(verbose) printf("mag_floors CASE 2 at (%d,%d,%d): %e %e\n",geom->ix+TOI,geom->iy+TOJ,geom->iz,pp[RHO],magpre);
-    ldouble f=magpre/(B2RHORATIOMAX*pp[RHO]);
+    f = magpre/(B2RHORATIOMAX*pp[RHO]);
+    rhofloored = 1;
+  }
 
-    for(iv=0;iv<NVMHD;iv++)
-    {
+  // check ugas vs bsq
+  if(magpre>B2UURATIOMAX*pp[UU]) {
+    if(verbose) printf("mag_floors CASE 3 at (%d,%d,%d): %e %e\n",geom->ix+TOI,geom->iy+TOJ,geom->iz,pp[UU],magpre);
+    fuu=magpre/(B2UURATIOMAX*pp[UU]); //increase factor, fuu>1   
+    uufloored = 1;
+  }
+
+   // floors were activated, apply them
+  if((rhofloored==1 || uufloored==1)){
+    // save backups
+    for(iv=0;iv<NVMHD;iv++){
       pporg[iv]=pp[iv];
     }
+    // Debora: updated to be consistent with driftframe
+    if (B2RHOFLOORFRAME==ZAMOFRAME){ //new mass in ZAMO
+      if(verbose) printf ("mag_floors in zamo frame\n");
+      ldouble dpp[NV],duu[NV];
 
-#if (B2RHOFLOORFRAME==ZAMOFRAME) //new mass in ZAMO
+      ldouble drho =pp[RHO]*(f-1.);
+         
+      for(iv=0;iv<NVMHD;iv++){
+        dpp[iv]=0.0;
+      }
 
-      ldouble dpp[NV],duu[NV];
-      ldouble drho=pp[RHO]*(f-1.);
-   
-      for(iv=0;iv<NVMHD;iv++)
-	dpp[iv]=0.0;
+      calc_normalobs_ncon(GG, geom->alpha, etacon);
+      conv_vels_ut(etacon,etarel,VEL4,VELPRIM,gg,GG);
 
-      dpp[RHO]=drho;
+      dpp[RHO] = drho;
       //do not inject energy - just density
-      dpp[UU]=0.;
+      dpp[UU] = 0;
+      
+      #ifdef DUINTFROMDRHO
+        if (fuu!=1){
+          ldouble duint = pp[UU]*drho/pp[RHO]; //Debora - from the acient thindisk works version
+          dpp[UU] = duint;
+        }
+      #endif
+      
       dpp[VX] = etarel[1];
       dpp[VY] = etarel[2];
       dpp[VZ] = etarel[3];
@@ -492,42 +521,112 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg)
 
       p2u_mhd(dpp,duu,geom);
  
-      for(iv=0;iv<NVMHD;iv++)
-      {
-	uu[iv]+=duu[iv];
+      for(iv=0;iv<NVMHD;iv++){
+	      uu[iv]+=duu[iv];
       }
 
       int rettemp=0;
       rettemp=u2p_solver(uu,pp,geom,U2P_HOT,0); 
       if(rettemp<0)
-	rettemp=u2p_solver(uu,pp,geom,U2P_ENTROPY,0); 
+	      rettemp=u2p_solver(uu,pp,geom,U2P_ENTROPY,0); 
       
-      if(rettemp<0) 
-      {
-#ifdef BHDISK_PROBLEMTYPE
-	//if(geom->ix+TOI>5) //report only outside horizon
-#endif
-	  printf("u2p failed after imposing bsq over rho floors at %d %d %d with gamma=%f\n",geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK,get_u_scalar(gammagas,geom->ix,geom->iy,geom->iz));
-
-#ifdef B2RHOFLOOR_BACKUP_FFFRAME
-       // Backup bsq/rho floor -- if zamo frame fails, do fluid frame instead of crashing 
-       for(iv=0;iv<NVMHD;iv++)
-         pp[iv]=pporg[iv];
-       pp[RHO]*=f;
-       pp[UU]*=f;
-#else
-       print_primitives(pp);
-       exit(-1);
-#endif
+      if(rettemp<0) {
+	      printf("u2p failed after imposing bsq over rho floors at %d %d %d with gamma=%f\n",geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK,get_u_scalar(gammagas,geom->ix,geom->iy,geom->iz));
+
+        #ifdef B2RHOFLOOR_BACKUP_FFFRAME
+          // Backup bsq/rho floor -- if zamo frame fails, do fluid frame instead of crashing 
+          for(iv=0;iv<NVMHD;iv++)
+            pp[iv]=pporg[iv];
+          pp[RHO]*=f;
+          pp[UU]*=fuu;
+        #else
+          print_primitives(pp);
+          exit(-1);
+        #endif
       }
+    }
     
-#elif(B2RHOFLOORFRAME==FFFRAME) //new mass in fluid frame
+    else if(B2RHOFLOORFRAME==FFFRAME){ //new mass in fluid frame
+      if(verbose) printf ("mag_floors in fluid frame\n");
       pp[RHO]*=f;
-      pp[UU]*=f;
+      pp[UU]*=fuu;
+    }
 
-#endif //B2RHOFLOORFRAME==ZAMOFRAME
+    //Debora new from Andrew's force free version, similar to HARM
+    else if(B2RHOFLOORFRAME==DRIFTFRAME){ // new mass in drift frame, Ressler+2017
+      if(verbose) printf ("mag_floors in drift frame\n");
+      ldouble betapar,betasq,betasqmax,gammapar;
+      ldouble udotB,QdotB,Bsq,Bmag,wold,wnew;
+      ldouble xx,vpar;
+      ldouble Bcon[4],Bcov[4],ucondr[4],vcon[4],ucont[4];
+      ldouble pgamma = GAMMA;
+      #ifdef CONSISTENTGAMMA
+        pgamma=pick_gammagas(geom->ix,geom->iy,geom->iz);
+      #endif
+      ldouble pgammam1=pgamma-1.;
+
+      // old enthalpy
+      wold = pporg[RHO] + pporg[UU]*pgamma;
+
+      // apply the floors to the scalars
+      pp[RHO] = pporg[RHO]*f;
+      //pp[UU] = pporg[UU]*fuu;
+
+      // new enthalpy
+      wnew = pp[RHO] + pp[UU]*pgamma;
+
+      // parallel velocity and Lorentz¨
+      betapar = -bcond[0]/((bsq)*ucond[0]);
+      betasq = betapar*betapar*bsq;
+      betasqmax = 1. - 1./(GAMMAMAXHD*GAMMAMAXHD);
+      if(betasq>betasqmax){
+        betasq=betasqmax;
+        printf("floored parallel velocity\n");
+      }
+      gammapar = 1./sqrt(1.-betasq);
+
+      // drift frame velocity
+      for(iv=0;iv<4;iv++)
+        ucondr[iv] = gammapar*(ucond[iv] + betapar*bcond[iv]);
+
+      // magnetic field 
+      Bcon[0]=0.;
+      Bcon[1]=pp[B1];
+      Bcon[2]=pp[B2];
+      Bcon[3]=pp[B3];
+      indices_21(Bcon,Bcov,gg);
+      Bsq = dot(Bcon,Bcov);
+      Bmag = sqrt(Bsq);
+
+      // conserved parallel momentum
+      udotB = dot(ucond,Bcov);
+      QdotB = udotB*wold*ucond[0];
+
+      // get new parallel velocity
+      xx = 2.*QdotB/(Bmag*wnew*ucondr[0]);
+      vpar = xx / (ucondr[0]*(1.+sqrt(1.+xx*xx)));
+
+      // new three velocity
+      vcon[0]=1.;
+      for(iv=1;iv<4;iv++){
+        vcon[iv]=vpar*Bcon[iv]/(Bmag) + ucondr[iv]/ucondr[0];
+      }
 
-#ifdef EVOLVEELECTRONS
+      // to VELPRIM
+      conv_vels(vcon,ucont,VEL3,VELPRIM,gg,GG);
+    
+      pp[VX] = ucont[1];
+      pp[VY] = ucont[2];
+      pp[VZ] = ucont[3];
+
+    }
+  
+
+    // ANDREW
+    // TODO this seems wrong for species?
+    // TODO if uint changes above, we will not have ue + ui = uint using this method!
+    // TODO multiply both ui,ue by fuu? 
+    #ifdef EVOLVEELECTRONS
 
       //keep energy density in ions and electrons fixed after modifying B
       ldouble Tg,Te,Ti,ptot,uint,theta;
@@ -538,20 +637,20 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg)
       ldouble ne=calc_thermal_ne(pporg); //thermal only
       ldouble pe=K_BOLTZ*ne*Te;
       ldouble gammae=GAMMAE;
-#ifdef CONSISTENTGAMMA
-#ifndef FIXEDGAMMASPECIES
-      gammae=calc_gammaintfromtemp(Te,ELECTRONS);
-#endif
-#endif
+      #ifdef CONSISTENTGAMMA
+      #ifndef FIXEDGAMMASPECIES
+            gammae=calc_gammaintfromtemp(Te,ELECTRONS);
+      #endif
+      #endif
       ldouble ue=pe/(gammae-1.);  
       ldouble ni=pporg[RHO]/MU_I/M_PROTON;
       ldouble pi=K_BOLTZ*ni*Ti;
       ldouble gammai=GAMMAI;
-#ifdef CONSISTENTGAMMA
-#ifndef FIXEDGAMMASPECIES
-      gammai=calc_gammaintfromtemp(Ti,IONS);
-#endif
-#endif
+      #ifdef CONSISTENTGAMMA
+      #ifndef FIXEDGAMMASPECIES
+            gammai=calc_gammaintfromtemp(Ti,IONS);
+      #endif
+      #endif
       ldouble ui=pi/(gammai-1.);
 
       //calculate new entropy using pp[]
@@ -562,13 +661,13 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg)
       mass = M_ELECTR;
       gamma = GAMMAE;
       n = calc_thermal_ne(pp);
-#ifdef CONSISTENTGAMMA
-      Tenew=solve_Teifromnmu(n, mass, ue,ELECTRONS); //solves in parallel for gamma and temperature
-      theta=K_BOLTZ*Tenew/mass;  
-#ifndef FIXEDGAMMASPECIES
-      gamma=calc_gammaintfromtheta(theta); //the same gamma as just solved     
-#endif
-#endif
+      #ifdef CONSISTENTGAMMA
+            Tenew=solve_Teifromnmu(n, mass, ue,ELECTRONS); //solves in parallel for gamma and temperature
+            theta=K_BOLTZ*Tenew/mass;  
+      #ifndef FIXEDGAMMASPECIES
+            gamma=calc_gammaintfromtheta(theta); //the same gamma as just solved     
+      #endif
+      #endif
       pnew=(ue)*(gamma-1.);
       Tenew=pnew/K_BOLTZ/n;
 
@@ -581,31 +680,24 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg)
       mass = M_PROTON;
       gamma = GAMMAI;
       n = calc_thermal_ne(pp);
-#ifdef CONSISTENTGAMMA
-      Tinew=solve_Teifromnmu(n, mass, ui,IONS); //solves in parallel for gamma and temperature
-      theta=K_BOLTZ*Tinew/mass;
-#ifndef FIXEDGAMMASPECIES
-      gamma=calc_gammaintfromtheta(theta); //the same gamma as just solved     
-#endif
-#endif
+      #ifdef CONSISTENTGAMMA
+            Tinew=solve_Teifromnmu(n, mass, ui,IONS); //solves in parallel for gamma and temperature
+            theta=K_BOLTZ*Tinew/mass;
+      #ifndef FIXEDGAMMASPECIES
+            gamma=calc_gammaintfromtheta(theta); //the same gamma as just solved     
+      #endif
+      #endif
       
       pnew=(ui)*(gamma-1.);
       Tinew=pnew/K_BOLTZ/n;
       Sinew=calc_SefromrhoT(pp[RHO],Tinew,IONS);
       
       pp[ENTRI]=Sinew;
-#endif  //EVOLVEELECTRONS
+    #endif  //EVOLVEELECTRONS
 
-      ret=-1;      
+    ret=-1;      
   } //if(magpre>B2RHORATIOMAX*pp[RHO]) 
 
-  //independent check on ugas vs 
-  if(magpre>B2UURATIOMAX*pp[UU]) 
-  {
-      //if(verbose) printf("mag_floors CASE 3 at (%d,%d,%d): %e %e\n",geom->ix+TOI,geom->iy+TOJ,geom->iz,pp[UU],magpre);
-      pp[UU]*=magpre/(B2UURATIOMAX*pp[UU]);
-      ret=-1;      
-  }
 
 #endif //MAGNFIELD
 
diff --git a/u2prad.c b/u2prad.c
index 70a95e8..6d8f14b 100644
--- a/u2prad.c
+++ b/u2prad.c
@@ -1,138 +1,11 @@
-/*! \file u2prad.c
- \brief Radiative routines related to u2p inversion
- */
-
+//KORAL - u2p.c
+//radiative routines related to u2p inversion
 
 #include "ko.h"
 
 
-static int get_m1closure_gammarel2_cold(int verbose, void *ggg,
-					FTYPE *Avcon,FTYPE *Avcov,
-					FTYPE *gammarel2return, FTYPE *deltareturn,
-					FTYPE *numeratorreturn, FTYPE *divisorreturn,
-					FTYPE *Erfreturn, FTYPE *urfconrel);
-static int get_m1closure_gammarel2(int verbose,void *ggg,
-				   ldouble *Avcon, ldouble *Avcov,
-				   ldouble *gammarel2return,ldouble *deltareturn,
-				   ldouble *numeratorreturn, ldouble *divisorreturn);
-static int get_m1closure_Erf(void *ggg, ldouble *Avcon, ldouble gammarel2, ldouble *Erfreturn);
-static int get_m1closure_urfconrel(int verbose, void *ggg,
-				   ldouble *pp, ldouble *Avcon,
-				   ldouble *Avcov, ldouble gammarel2,
-				   ldouble delta, ldouble numerator,
-				   ldouble divisor, ldouble *Erfreturn,
-				   ldouble *urfconrel, int *corflag);
-
-//**********************************************************************
-//Do U2P including for photon number
-//**********************************************************************
-
-int
-u2p_rad(ldouble *uu, ldouble *pp, void *ggg, int *corrected)
-{
-  //whether primitives corrected for caps, floors etc. - if so, conserved will be updated
-  *corrected=0;
-
-  int u2pret=u2p_rad_urf(uu,pp,ggg,corrected);
-
-  #ifdef EVOLVEPHOTONNUMBER
-  struct geometry *geom
-   = (struct geometry *) ggg;
-  ldouble gdetu=geom->gdet;
-  #if (GDETIN==0) //gdet out of derivatives
-  gdetu=1.;
-  #endif
-  ldouble urfcon[4];
-  urfcon[0]=0.;
-  urfcon[1]=pp[FX];
-  urfcon[2]=pp[FY];
-  urfcon[3]=pp[FZ];
-  conv_vels(urfcon,urfcon,VELPRIMRAD,VEL4,geom->gg,geom->GG);
- 
-  pp[NF]=uu[NF]/urfcon[0]/gdetu;
-  #endif
-
-  return u2pret;
-}
-
-//**********************************************************************
-//basic conserved to primitives solver for radiation
-//uses M1 closure in arbitrary frame/metric
-//radiative primitives: (E,\tilde u^i)
-//  E - radiative energy density in the rad.rest frame
-//  u^i - relative velocity of the rad.rest frame
-//takes conserved R^t_mu in uu
-//**********************************************************************
-
-int
-u2p_rad_urf(ldouble *uu, ldouble *pp,void* ggg, int *corrected)
-{
-  struct geometry *geom
-    = (struct geometry *) ggg;
-
-  ldouble (*gg)[5],(*GG)[5],gdet,gdetu;
-  gg=geom->gg;
-  GG=geom->GG;
-  gdet=geom->gdet;
-  gdetu=gdet;
-
-  #if (GDETIN==0) //gdet out of derivatives
-  gdetu=1.;
-  #endif
-
-  int verbose=0;
-  
-  //whether primitives corrected for caps, floors etc. - if so, conserved will be updated 
-  *corrected=0;
-
-  ldouble urfcon[4],Erf;
-
-  //conserved - R^t_mu
-  ldouble Avcov[4]={uu[EE]/gdetu,uu[FX]/gdetu,uu[FY]/gdetu,uu[FZ]/gdetu};
-  ldouble Avcon[4];
-
-  //indices up - R^tmu
-  indices_12(Avcov,Avcon,GG);
-
-  ldouble gammarel2,delta,numerator,divisor;
-
-  // get \gamma^2 for relative 4-velocity
-  if(get_m1closure_gammarel2(verbose,ggg,Avcon,Avcov,&gammarel2,&delta,&numerator,&divisor)<0)
-    {
-      //printf("get_m1closure_gammarel2 failed at %d %d %d\n",geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK);
-      return -1;
-    }
-
-  // get E in radiation frame
-  get_m1closure_Erf(ggg,Avcon,gammarel2,&Erf);
-
-  // get relative 4-velocity
-  if(get_m1closure_urfconrel(verbose,ggg,pp,Avcon,Avcov,gammarel2,delta,numerator,divisor,&Erf,urfcon,corrected)<0)
-    {
-      //printf("get_m1closure_urfconrel failed at %d %d %d\n",geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK);
-      return -1;
-    }
-
-  if (VELR != VELPRIMRAD)
-  {
-    conv_vels(urfcon,urfcon,VELR,VELPRIMRAD,gg,GG);
-  }
-
-  //new primitives
-  pp[EE]=Erf;
-  pp[FX]=urfcon[1];
-  pp[FY]=urfcon[2];
-  pp[FZ]=urfcon[3];
-
-  return 0;
-}
-
-
-//**********************************************************************
-//gets gamma assuming fixed E rather than using original R^t_t that we assume is flawed near floor regions.
-//We want to preserve R^t_i (i.e conserved momentum)
-//**********************************************************************
 
+// get's gamma assuming fixed E rather than using original R^t_t that we assume is flawed near floor regions.  We want to preserve R^t_i (i.e conserved momentum)
 static int get_m1closure_gammarel2_cold(int verbose,
 					void *ggg, 
 					FTYPE *Avcon, 
@@ -157,8 +30,7 @@ static int get_m1closure_gammarel2_cold(int verbose,
   int jj;
 
   //static
-  FTYPE gctt, gn11, gn12,  gn13,  gn14,  gn22,  gn23,  gn24,  gn33,  gn34,  gn44;
-  FTYPE Rtt,  Rtx,  Rty,  Rtz,  Rdtt,  Rdtx,  Rdty,  Rdtz;
+  FTYPE gctt, gn11, gn12,  gn13,  gn14,  gn22,  gn23,  gn24,  gn33,  gn34,  gn44,  Rtt,  Rtx,  Rty,  Rtz,  Rdtt,  Rdtx,  Rdty,  Rdtz;
   gn11=GG[0][0];
   gn12=GG[0][1];
   gn13=GG[0][2];
@@ -229,18 +101,11 @@ static int get_m1closure_gammarel2_cold(int verbose,
   return(0);
 }
 
-//**********************************************************************
-// gets gamma^2 for lab-frame gamma using Rd and gcon
-//**********************************************************************
 
-static int get_m1closure_gammarel2(int verbose,
-				   void *ggg,
-				   ldouble *Avcon,
-				   ldouble *Avcov,
-				   ldouble *gammarel2return,
-				   ldouble *deltareturn,
-				   ldouble *numeratorreturn,
-				   ldouble *divisorreturn)
+
+
+// get's gamma^2 for lab-frame gamma  using Rd and gcon
+static int get_m1closure_gammarel2(int verbose,void *ggg, ldouble *Avcon, ldouble *Avcov, ldouble *gammarel2return,ldouble *deltareturn, ldouble *numeratorreturn, ldouble *divisorreturn)
 {
   struct geometry *geom
    = (struct geometry *) ggg;
@@ -252,10 +117,9 @@ static int get_m1closure_gammarel2(int verbose,
   ldouble gamma2,gammarel2,delta,numerator,divisor;
   ldouble gamma2a,gamma2b;
 
-  // mathematica solution that avoids catastrophic cancellation when Rtt very small
-  // (otherwise above gives gamma2=1/2 oddly when gamma2=1) -- otherwise same as above
+  // mathematica solution that avoids catastrophic cancellation when Rtt very small (otherwise above gives gamma2=1/2 oddly when gamma2=1) -- otherwise same as above
+  // well, then had problems for R~1E-14 for some reason when near BH.  Couldn't quickly figure out, so use no replacement of gv11.
   // see u2p_inversion.nb
-
   //static
   ldouble gctt, gn11, gn12,  gn13,  gn14,  gn22,  gn23,  gn24,  gn33,  gn34,  gn44,  Rtt,  Rtx,  Rty,  Rtz,  Rdtt,  Rdtx,  Rdty,  Rdtz;
   gn11=GG[0][0];
@@ -293,8 +157,7 @@ static int get_m1closure_gammarel2(int verbose,
 
   if(verbose>1)     printf("gamma2a: %e\n",gamma2a);
 
-  if( gamma2a<GAMMASMALLLIMIT || !isfinite(gamma2a) )
-  {
+  if( gamma2a<GAMMASMALLLIMIT || !isfinite(gamma2a) ){
     gamma2b=(0.25*(-2.*Power(gn11,2)*Power(Rdtt,2) - 1.*gn11*(4.*gn12*Rdtt*Rdtx + gn22*Power(Rdtx,2) + 
                                                               Rdty*(4.*gn13*Rdtt + 2.*gn23*Rdtx + gn33*Rdty) + 2.*(2.*gn14*Rdtt + gn24*Rdtx + gn34*Rdty)*Rdtz + gn44*Power(Rdtz,2)) + 
                    gn11*Rdtt*Sqrt(4.*Power(gn11,2)*Power(Rdtt,2) + Power(gn12*Rdtx + gn13*Rdty + gn14*Rdtz,2) + 
@@ -310,8 +173,8 @@ static int get_m1closure_gammarel2(int verbose,
 
     if(verbose)     printf("gamma2b: %e\n",gamma2b);
   }
-  else
-  {
+  else{
+    // choose
     gamma2=gamma2a;
   }
 
@@ -321,14 +184,22 @@ static int get_m1closure_gammarel2(int verbose,
   //  if(!isnan(gamma2a) && isnan(gamma2b))
   //    gamma2=gamma2a;
 
-  //********************************************
+  ////////////////////////
+  //
   //cap on u^t
-  //********************************************
-
+  //
+  ///////////////////////
   ldouble alpha=geom->alpha;
 
+
   // get relative 4-velocity, that is always >=1 even in GR
   gammarel2 = gamma2*alpha*alpha;
+
+  //to consider it a failure - leads to instability
+  //if(isnan(gammarel2)) return -1;
+  
+  // check for machine error away from 1.0 that happens sometimes
+  //if(isnan(gammarel2) || (gammarel2>GAMMASMALLLIMIT && gammarel2<1.0)) - WTF?
   
   if((gammarel2>GAMMASMALLLIMIT && gammarel2<1.0))
   {
@@ -337,6 +208,17 @@ static int get_m1closure_gammarel2(int verbose,
     gammarel2=1.0;
   }
 
+  /*
+  if(isnan(gammarel2))
+    {
+      printf("nan in get_m1closure_gammarel2\n");
+      printf("%e %e %e\n",gamma2,gamma2a,gamma2b);
+      print_4vector(Avcon);
+      print_4vector(Avcov);
+      //      getchar();
+    }
+  */
+
   *gammarel2return=gammarel2; 
   *deltareturn=delta=0;
   *numeratorreturn=numerator=0;
@@ -344,10 +226,8 @@ static int get_m1closure_gammarel2(int verbose,
   return(0);
 }
 
-//**********************************************************************
-// get Erf
-//**********************************************************************
 
+// get Erf
 static int get_m1closure_Erf(void *ggg, ldouble *Avcon, ldouble gammarel2, ldouble *Erfreturn)
 {
   struct geometry *geom
@@ -355,18 +235,28 @@ static int get_m1closure_Erf(void *ggg, ldouble *Avcon, ldouble gammarel2, ldoub
 
   ldouble alpha=geom->alpha;
 
+  ////////////
+  //
   // get initial attempt for Erf
-  // If delta<0, then gammarel2=nan and Erf<RADLIMIT check below will fail 
+  // If delta<0, then gammarel2=nan and Erf<RADLIMIT check below will fail as good.
+  //
+  ////////////
   *Erfreturn = 3.*Avcon[0]*alpha*alpha/(4.*gammarel2-1.0);  // JCM
 
+  /*
+  if(isnan(*Erfreturn)) 
+    {
+      printf("nan in get_m1closure_Erf\n %d %d %d %e %e\n",geom->ix,geom->iy,geom->iz,*Erfreturn,gammarel2);
+      print_4vector(Avcon);
+      getchar();
+    }
+  */
+
   return(0);
 }
 
-
-//**********************************************************************
-// get contravariant rad frame 4-velocity in lab frame
-//**********************************************************************
-static int get_m1closure_urfconrel(int verbose,
+// get contravariant relative 4-velocity in lab frame
+static int get_m1closure_urfconrel(int verbose, 
 				   void *ggg, 
 				   ldouble *pp, 
 				   ldouble *Avcon, 
@@ -393,53 +283,53 @@ static int get_m1closure_urfconrel(int verbose,
   ldouble gammamax=GAMMAMAXRAD; 
   int ii,jj,kk,usingfast;
 
-  // Fix-up inversion if problem with gamma (i.e. velocity)
-  // or energy density in radiation rest-frame (i.e. Erf)
+  //////////////////////
+  //
+  // Fix-up inversion if problem with gamma (i.e. velocity) or energy density in radiation rest-frame (i.e. Erf)
+  //
+  //////////////////////
+
   // NOTE: gammarel2 just below 1.0 already fixed to be =1.0
-  int nonfailure = gammarel2>=1.0 && Erf>ERADFLOOR && gammarel2<=gammamax*gammamax/GAMMASMALLLIMIT/GAMMASMALLLIMIT;
- 
+  int nonfailure=gammarel2>=1.0 && Erf>ERADFLOOR && gammarel2<=gammamax*gammamax/GAMMASMALLLIMIT/GAMMASMALLLIMIT;
   // falilure1 : gammarel2 normal, but already Erf<ERADFLOOR (note for M1 that gammarel2>=1/4 for any reasonable chance for correct non-zero Erf
   int failure1=Avcon[0]<0.0 || (gammarel2>0.0 && gammarel2<=0.25 && delta>=0.0 && divisor!=0.0) || numerator==0.0 || (gammarel2>=1.0 && delta>=0.0 && divisor!=0.0 && Erf<ERADFLOOR);
-
+  //int failure1=Avcon[0]<0.0 || gammarel2>0.0 && gammarel2<=0.25 && delta>=0.0 && divisor!=0.0 || numerator==0.0 || gammarel2>=1.0 && delta>=0.0 && divisor!=0.0 && Erf<ERADFLOOR;
   // gamma probably around 1
   int failure2=gammarel2<1.0 && gammarel2>0.0 && delta>=0.0;
-
   // i.e. all else, so not really used below.
   int failure3=(gammarel2>gammamax*gammamax && Erf>=ERADFLOOR) || gammarel2<0.0 || delta<0.  || (divisor==0.0 && numerator==0.0) || (divisor==0.0 && numerator!=0.0);
+  //int failure3=gammarel2>gammamax*gammamax && Erf>=ERADFLOOR || gammarel2<0.0 || delta<0.  || divisor==0.0 && numerator==0.0 || divisor==0.0 && numerator!=0.0;
 
   // any failure
   int failure=!nonfailure || isinf(gammarel2) || isinf(Erf);
 
-  if(failure && (failure1==0 && failure2==0 && failure3==0))
-  {
+  if(failure && (failure1==0 && failure2==0 && failure3==0)){
     printf("Undetected failure, now considered\n");
   }
 
-  if(nonfailure)
-  {
+
+  if(nonfailure){
     // get good relative velocity
     ldouble gammarel=sqrt(gammarel2);
  
     for(ii=0;ii<4;ii++)
-    {	  
-      urfconrel[ii] = alpha * (Avcon[ii] + 1./3.*Erf*GG[0][ii]*(4.0*gammarel2-1.0) )/(4./3.*Erf*gammarel);
-    }
+      {	  
+	urfconrel[ii] = alpha * (Avcon[ii] + 1./3.*Erf*GG[0][ii]*(4.0*gammarel2-1.0) )/(4./3.*Erf*gammarel);
+      }
 
     *Erfreturn=Erf; // pass back new Erf to pointer
     *corflag=0;
     return(0);
   }
-  else //failure -- try getting cold gammarel2
-  {
+  else{
     if(verbose) 
       {
-	printf("failure: %d %d %d %d %d\n",nonfailure,failure1,failure2,failure3,failure);
-	printf("%e %e %e\n",Erf,gammarel2,Avcov[0]);
+	printf("failure: %d %d %d %d %d at %d %d %d\n",nonfailure,failure1,failure2,failure3,failure,geom->ix,geom->iy,geom->iz);
+  printf("%e %e %e\n",Erf,gammarel2,Avcov[0]);
 	print_4vector(Avcon);
 	print_4vector(Avcov);
 	if(isnan(Erf))	getchar();
       }
-    
     // get \gammarel=1 case
     ldouble gammarel2slow=pow(1.0+10.0*NUMEPSILON,2.0);
     ldouble Avconslow[4],Avcovslow[4],urfconrelslow[4];
@@ -463,10 +353,9 @@ static int get_m1closure_urfconrel(int verbose,
     get_m1closure_gammarel2_cold(verbose,ggg,Avconfast,Avcovfast,&gammarel2fast,&delta,&numerator,&divisor,&Erffast,urfconrelfast);
 
     usingfast=1;
-
-    // choose by which Avcov[0] is closest to original
-    if(fabs(Avcovslow[0]-Avcov[0])>fabs(Avcovfast[0]-Avcov[0]))
-    {
+    // choose by which Avcov[0] is closest to original&&
+    //    if( fabs(Avcovslow[0]-Avcov[0])>fabs(Avcovfast[0]-Avcov[0])){
+    if( fabs(Avcovslow[0]-Avcov[0])>fabs(Avcovfast[0]-Avcov[0])){
       usingfast=1;
       for(jj=0;jj<4;jj++)
 	{
@@ -477,8 +366,7 @@ static int get_m1closure_urfconrel(int verbose,
       gammarel2=gammarel2fast;
       Erf=Erffast;
     }
-    else
-    {
+    else{
       usingfast=0;
       for(jj=0;jj<4;jj++)
 	{
@@ -495,10 +383,14 @@ static int get_m1closure_urfconrel(int verbose,
     if(verbose) printf("CASEGEN: gammarel>gammamax (cold, usingfast=%d): gammarel2=%g Erf=%g : i=%d j=%d k=%d\n",usingfast,gammarel2,Erf,geom->ix,geom->iy,geom->iz);
   }
 
+
   *Erfreturn=Erf; // pass back new Erf to pointer
 
-  if(verbose) printf("end: %e %e %e %e\n",Erf,*Erfreturn,Erfslow,Erffast);
-      
+  if(verbose) 
+      {
+	printf("end: %e %e %e %e\n",Erf,*Erfreturn,Erfslow,Erffast);
+      }
+
   if(!isfinite(Erf) || !isfinite(gammarel2) || !isfinite(urfconrel[1])|| !isfinite(urfconrel[2])|| !isfinite(urfconrel[3]))
     {
       if(verbose)      
@@ -507,22 +399,286 @@ static int get_m1closure_urfconrel(int verbose,
 	  mpi_local2globalidx(geom->ix,geom->iy,geom->iz,&gix,&giy,&giz);
 	  printf("------------\nJONNAN: %e %e %e %e\n",Erf,*Erfreturn,Erfslow,Erffast);
 	  if(usingfast) printf("JONNAN: usingfast==1\n"); else printf("JONNAN: usingfast==0\n");
+	  
 	  printf("JONNAN: ijk=%d %d %d :  %g %g : %g %g %g %g : %d %d %d %d : %g %g %g %g\n",gix,giy,giz,Erf,gammarel2,urfconrel[0],urfconrel[1],urfconrel[2],urfconrel[3],failure1,failure2,failure3,failure,Avcon[0],Avcon[1],Avcon[2],Avcon[3]);
 	}
     
 
       return -1;
-    }
+  }
 
   return(0);
 }
 
 
+
 //**********************************************************************
 //**********************************************************************
-//checks if rad primitives make sense
+//basic conserved to primitives solver for radiation
+//uses M1 closure in arbitrary frame/metric
+//radiative primitives: (E,\tilde u^i)
+//  E - radiative energy density in the rad.rest frame
+//  u^i - relative velocity of the rad.rest frame
+//takes conserved R^t_mu in uu
+//**********************************************************************
 //**********************************************************************
+int
+u2p_rad_urf(ldouble *uu, ldouble *pp,void* ggg, int *corrected)
+{
+  struct geometry *geom
+    = (struct geometry *) ggg;
+
+  ldouble (*gg)[5],(*GG)[5],gdet,gdetu;
+  gg=geom->gg;
+  GG=geom->GG;
+  gdet=geom->gdet;gdetu=gdet;
+#if (GDETIN==0) //gdet out of derivatives
+  gdetu=1.;
+#endif
+
+  int verbose=0;
+  
+  //whether primitives corrected for caps, floors etc. - if so, conserved will be updated 
+  *corrected=0;
+
+  ldouble Rij[4][4];
+  ldouble urfcon[4],urfcov[4],Erf;
+  //conserved - R^t_mu
+  ldouble Avcov[4]={uu[EE]/gdetu,uu[FX]/gdetu,uu[FY]/gdetu,uu[FZ]/gdetu};
+  ldouble Avcon[4];
+  //indices up - R^tmu
+  indices_12(Avcov,Avcon,GG);
+
+  ldouble gammarel2,delta,numerator,divisor;
+
+  // get \gamma^2 for relative 4-velocity
+  if(get_m1closure_gammarel2(verbose,ggg,Avcon,Avcov,&gammarel2,&delta,&numerator,&divisor)<0)
+    {
+      //printf("get_m1closure_gammarel2 failed at %d %d %d\n",geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK);
+      return -1;
+    }
+
+  // get E in radiation frame
+  get_m1closure_Erf(ggg,Avcon,gammarel2,&Erf);
+
+  // get relative 4-velocity
+  if(get_m1closure_urfconrel(verbose,ggg,pp,Avcon,Avcov,gammarel2,delta,numerator,divisor,&Erf,urfcon,corrected)<0)
+    {
+      //printf("get_m1closure_urfconrel failed at %d %d %d\n",geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK);
+      return -1;
+    }
+
+  conv_vels(urfcon,urfcon,VELR,VELPRIMRAD,gg,GG);
+
+  //new primitives
+  pp[EE]=Erf;
+  pp[FX]=urfcon[1];
+  pp[FY]=urfcon[2];
+  pp[FZ]=urfcon[3];
+
+
+  return 0;
+}
+
 
+int
+u2p_rad(ldouble *uu, ldouble *pp, void *ggg, int *corrected)
+{
+  //whether primitives corrected for caps, floors etc. - if so, conserved will be updated
+  *corrected=0;
+
+  int u2pret=u2p_rad_urf(uu,pp,ggg,corrected);
+
+  #ifdef EVOLVEPHOTONNUMBER
+  struct geometry *geom
+   = (struct geometry *) ggg;
+  ldouble gdetu=geom->gdet;
+  #if (GDETIN==0) //gdet out of derivatives
+  gdetu=1.;
+  #endif
+  ldouble urfcon[4];
+  urfcon[0]=0.;
+  urfcon[1]=pp[FX];
+  urfcon[2]=pp[FY];
+  urfcon[3]=pp[FZ];
+  conv_vels(urfcon,urfcon,VELPRIMRAD,VEL4,geom->gg,geom->GG);
+ 
+  pp[NF]=uu[NF]/urfcon[0]/gdetu;
+  #endif
+
+  //M1
+  return u2pret;
+}
+
+
+//**********************************************************************
+//**********************************************************************
+//**********************************************************************
+//numerical conserved to primitives solver for radiation
+//works in ortonormal fluid frame
+//used e.g. for not-frame-invariant  Eddington apr. 
+//solves in 4dimensions using frame boosts etc.
+int f_u2prad_num(ldouble *uu,ldouble *pp, void* ggg,ldouble *f)
+{
+  /*struct geometry *geom
+    = (struct geometry *) ggg;
+
+  ldouble (*gg)[5],(*GG)[5],(*tlo)[4],(*tup)[4];
+  gg=geom->gg;
+  GG=geom->GG;
+  tlo=geom->tlo;
+  tup=geom->tup;
+
+  ldouble Rij[4][4];
+
+  calc_Rij_M1_ff(pp,Rij);  
+  trans22_on2cc(Rij,Rij,tlo);  
+  boost22_ff2lab(Rij,Rij,pp,gg,GG); 
+  indices_2221(Rij,Rij,gg);  
+
+  f[0]=-Rij[0][0]+uu[EE0];
+  f[1]=-Rij[0][1]+uu[FX0];
+  f[2]=-Rij[0][2]+uu[FY0];
+  f[3]=-Rij[0][3]+uu[FZ0];
+*/
+  return 0;
+} 
+
+int
+print_state_u2prad_num (int iter, ldouble *x, ldouble *f)
+{
+  printf ("iter = %3d x = % .3e % .3e % .3e % .3e "
+	  "f(x) = % .3e % .3e % .3e % .3e\n",
+	  iter,
+	  x[0],x[1]/x[0],x[2]/x[0],x[3]/x[0],f[0],f[1],f[2],f[3]);
+  return 0;
+}
+
+int
+u2p_rad_onff(ldouble *uu, ldouble *pp, void* ggg, int *corrected)
+{
+  struct geometry *geom
+    = (struct geometry *) ggg;
+
+  ldouble pp0[NV],pporg[NV];
+  ldouble J[4][4],iJ[4][4];
+  ldouble x[4],f1[4],f2[4],f3[4];
+  int i,j,k,iter=0;
+
+  ldouble EPS = 1.e-6;
+  ldouble CONV = 1.e-6;
+
+  int verbose=0;
+
+  for(i=EE0;i<NV;i++)
+    {
+      pporg[i]=pp[i];
+    }
+
+  //converting radiative primitives to fluid frame ortonormal
+  //ad_lab2ff(pp,pp,geom);
+
+  if(verbose!=0)   print_Nvector(uu,NV);
+  do
+    {
+      iter++;
+      for(i=EE0;i<NV;i++)
+	{
+	  pp0[i]=pp[i];
+	}
+
+      //valueas at zero state
+      f_u2prad_num(uu,pp,geom,f1);
+ 
+      //calculating approximate Jacobian
+      for(j=0;j<4;j++)
+	{
+	  pp[j+EE0]=pp[j+EE0]+EPS*pp[EE0];
+	    
+	  if(verbose>0)    print_Nvector(pp,NV);
+ 	  f_u2prad_num(uu,pp,geom,f2);
+	  if(verbose>0)    print_state_u2prad_num (iter,x,f2); 
+     
+	  for(i=0;i<4;i++)
+	    {
+	      J[i][j]=(f2[i] - f1[i])/(EPS*pp[EE0]);
+	    }
+
+	  pp[j+EE0]=pp0[j+EE0];
+	}
+	
+
+      //inversion
+      inverse_44matrix(J,iJ);
+
+      //updating x
+      for(i=0;i<4;i++)
+	{
+	  x[i]=pp0[i+EE0];
+	}
+
+      for(i=0;i<4;i++)
+	{
+	  for(j=0;j<4;j++)
+	    {
+	      x[i]-=iJ[i][j]*f1[j];
+	    }
+	}
+      if(verbose>0)    print_state_u2prad_num (iter,x,f1); 
+
+      for(i=0;i<4;i++)
+	{
+	  pp[i+EE0]=x[i];
+	}
+  
+      //test convergence
+      for(i=0;i<4;i++)
+	{
+	  f3[i]=(pp[i+EE0]-pp0[i+EE0]);
+	  f3[i]=fabs(f3[i]/pp0[EE0]);
+	}
+
+      if(f3[0]<CONV && f3[1]<CONV && f3[2]<CONV && f3[3]<CONV)
+	break;
+
+      if(iter>20)
+	{
+	  printf("iter exceeded in u2prad_num() %d %d %d\n",geom->ix,geom->iy,geom->iz);getchar();
+
+	  pp[EE0]=pporg[EE0];
+	  pp[FX0]=pp[FY0]=pp[FZ0]=0.;
+	  
+	  *corrected=1;
+	  return -1;
+
+	  break;
+	}
+     
+    }
+  while(1);
+  
+  if(pp[EE0]<EEFLOOR) 
+    {
+      printf("enegative u2prad() %d %d %d\n",geom->ix,geom->iy,geom->iz); getchar();
+      pp[EE0]=EEFLOOR;
+      pp[7]=pp[8]=pp[9]=0.;
+      *corrected=1;
+      return -1;
+    }
+  
+  //converting to lab primitives
+  //no - in EDD_APR I use fluid frame fluxes as primitives
+  //prad_ff2lab(pp,pp,geom);
+  
+  if(verbose!=0)   {print_Nvector(pp,NV);}
+  if(verbose>0)   {printf("----\n");}//getchar();}
+
+  *corrected=0;
+  return 0;
+
+}
+
+/// Debora - taken from koral_lite to count in photon number
 int
 check_floors_rad(ldouble *pp, int whichvel,void *ggg)
 {
@@ -697,3 +853,123 @@ check_floors_rad(ldouble *pp, int whichvel,void *ggg)
   return ret;
 
 }
+
+
+
+// Debora - old function  
+//**********************************************************************
+//**********************************************************************
+//**********************************************************************
+//checks if rad primitives make sense
+/*int
+check_floors_rad(ldouble *pp, int whichvel,void *ggg)
+{
+  //skip floors for some time
+  //return 0;
+
+  int verbose=0;
+  int ret=0;
+
+  struct geometry *geom
+    = (struct geometry *) ggg;
+
+  ldouble (*gg)[5],(*GG)[5];
+  gg=geom->gg;
+  GG=geom->GG;
+
+
+#ifdef RADIATION  
+  ldouble pp2[NV];
+  int iv;
+
+  //ff rad energy density
+  ldouble Rtt,Ehat,ugas[4],Eratio;
+  calc_ff_Rtt(pp,&Rtt,ugas,geom);
+  Ehat=-Rtt;
+  Eratio=pp[EE0]/Ehat; //ratio of energy densities in two frames
+
+  //absolute rad-frame EE:
+  if(pp[EE0]<ERADFLOOR) 
+    {
+      if(verbose) printf("rhd_floors CASE R0 at (%d,%d,%d): %e %e\n",geom->ix,geom->iy,geom->iz,pp[0],pp[EE0]);
+      pp[EE0]=ERADFLOOR;
+      ret=-1;
+     }
+
+  #ifndef SKIPRADSOURCE
+
+  //Ehat/rho ratios 
+  if(Ehat<EERHORATIOMIN*pp[0]) 
+    {
+      if(verbose) printf("hd_floors CASE R2 at (%d,%d,%d): %e %e\n",geom->ix,geom->iy,geom->iz,pp[0],pp[EE0]);
+      pp[EE0]=Eratio*EERHORATIOMIN*pp[0];
+      ret=-1;
+    }
+
+  if(Ehat>EERHORATIOMAX*pp[0]) 
+    {
+      if(verbose) printf("hd_floors CASE R3 at (%d,%d,%d): %e %e\n",geom->ix,geom->iy,geom->iz,pp[0],Ehat);
+      pp[0]=1./EERHORATIOMAX*Ehat;
+      ret=-1;
+    }
+  
+  //Ehat/uint ratios  
+  if(Ehat<EEUURATIOMIN*pp[1]) 
+    {
+      if(verbose) printf("hd_floors CASE R4 at (%d,%d,%d): %e %e\n",geom->ix,geom->iy,geom->iz,pp[1],Ehat);
+      pp[EE0]=Eratio*EEUURATIOMIN*pp[1];
+      ret=-1;
+    }
+
+  if(Ehat>EEUURATIOMAX*pp[1]) 
+    {
+      if(verbose) printf("hd_floors CASE R5 at (%d,%d,%d): %e %e\n",geom->ix,geom->iy,geom->iz,pp[1],Ehat);
+      pp[1]=1./EEUURATIOMAX*Ehat;
+      ret=-1;
+    }
+
+
+#ifdef SKIP_MAGNFIELD
+
+  ldouble ucond[4],ucovd[4];
+  ldouble bcond[4],bcovd[4],magpre;  
+  ldouble etacon[4],etarel[4];
+  for(iv=1;iv<4;iv++)
+    ucond[iv]=pp[1+iv];
+  conv_vels(ucond,ucond,VELPRIM,VEL4,gg,GG);
+ 
+  indices_21(ucond,ucovd,gg);
+  calc_bcon_4vel(pp,ucond,ucovd,bcond);
+  indices_21(bcond,bcovd,gg); 
+  magpre = dotB(bcond,bcovd)/2.;
+
+  //Ehat/uint ratios 
+  
+  if(magpre>B2EERATIOMAX*Ehat) 
+    {
+      if(verbose) printf("rad_floors CASE MR4 at (%d,%d,%d): %e %e\n",geom->ix,geom->iy,geom->iz,magpre,Ehat);
+      pp[EE0]=Eratio*magpre/B2EERATIOMAX;
+      ret=-1;
+    }
+  
+
+
+#endif
+
+#endif
+
+#ifdef PUTNFFLOOR
+  //floor on number of photons - does not work in general
+  if(pp[NF]<SMALL)
+    {
+      //printf("NF floor imposed at  (%d,%d,%d): %e %e\n",geom->ix,geom->iy,geom->iz,pp[NF],SMALL);
+      pp[NF]=SMALL;
+    }
+#endif
+#endif
+ 
+
+  return ret;
+
+}
+*/
\ No newline at end of file