                     dg::blas1::pointwiseDot( transferH2dX, volX2d, transferH2dX); //multiply by sqrt(g)
                     poloidal_average( transferH2dX, t1d, false); //average over eta
                     dg::blas1::scal( t1d, 4*M_PI*M_PI*f0); //
-                    dg::blas1::pointwiseDivide( t1d, dvdpsip, fsa1d );
-                    if( record_name[0] == 'j')
-                        dg::blas1::pointwiseDot( fsa1d, dvdpsip, fsa1d );
+                    dg::blas1::copy( 0., fsa1d); //get rid of previous nan in fsa1d (nasty bug)
+                    if( record_name[0] != 'j')
+                        dg::blas1::pointwiseDivide( t1d, dvdpsip, fsa1d );
+                    else
+                        dg::blas1::copy( t1d, fsa1d);
                     //3. Interpolate fsa on 2d plane : <f>
                     dg::blas2::gemv(fsa2rzmatrix, fsa1d, transferH2d); //fsa on RZ grid
                 }

                err = nc_put_vara_double( ncid_out, id2d.at(record_name+"_fsa2d"),
                    start2d_out, count2d, transferH2d.data() );
                if( record_name[0] == 'j')
                    record_name[1] = 's';
                err = nc_put_vara_double( ncid_out, id2d.at(record_name+"_cta2d"),
                    start2d_out, count2d, t2d_mp.data() ); //this is still js...
                if( record_name[0] == 'j')
                    record_name[1] = 'v';
                //4. Read 2d variable and compute fluctuations
                available = true;
                try{

            long_name.data());

        name = record_name + "_cta2d";
        long_name = record.long_name + " (Convoluted toroidal average on 2d plane.)";
        err = nc_def_var( ncid_out, name.data(), NC_DOUBLE, 3, dim_ids,
            &id2d[name]);
