diff --git a/.github/workflows/idefix-ci-jobs.yml b/.github/workflows/idefix-ci-jobs.yml
index dfb89e0f..b95385e4 100644
--- a/.github/workflows/idefix-ci-jobs.yml
+++ b/.github/workflows/idefix-ci-jobs.yml
@@ -177,6 +177,8 @@ jobs:
         run: scripts/ci/run-tests $IDEFIX_DIR/test/MHD/sphBragTDiffusion -all $TESTME_OPTIONS
       - name: Spherical anisotropic viscosity
         run: scripts/ci/run-tests $IDEFIX_DIR/test/MHD/sphBragViscosity -all $TESTME_OPTIONS
+      - name: Collisionless thermal conduction
+        run: scripts/ci/run-tests $IDEFIX_DIR/test/MHD/clessTDiffusion -all $TESTME_OPTIONS
 
   Examples:
     needs: [Fargo, Dust, Planet, ShearingBox, SelfGravity]
diff --git a/reference b/reference
index 5ecb9588..c3289420 160000
--- a/reference
+++ b/reference
@@ -1 +1 @@
-Subproject commit 5ecb95884f845a23fa0baf45492770c30e22b4fc
+Subproject commit c328942083b618a85b84eb16ce9fd35e67c00597
diff --git a/test/MHD/clessTDiffusion/setup.cpp b/test/MHD/clessTDiffusion/setup.cpp
index 651a81f5..c540c0de 100644
--- a/test/MHD/clessTDiffusion/setup.cpp
+++ b/test/MHD/clessTDiffusion/setup.cpp
@@ -45,7 +45,7 @@ void UserDefBoundary(Hydro *hydro, int dir, BoundarySide side, real t) {
     IdefixArray4D<real> Vs = hydro->Vs;
     IdefixArray1D<real> x1 = data->x[IDIR];
 
-    real rc,vc,vwind0;
+    real rc,vwind0;
     real cs=cs_vescGlob*sqrt(2.);
     real va_vesc = va_vescGlob;
     real mu = va_vesc * sqrt(2.);
@@ -53,7 +53,6 @@ void UserDefBoundary(Hydro *hydro, int dir, BoundarySide side, real t) {
 
     PonRho = cs*cs;
     rc = 0.25 / (cs_vescGlob*cs_vescGlob);
-    vc = cs;
     vwind0 = ParkerWind(1./rc) * cs;
 
     hydro->boundary->BoundaryFor("UserDefBoundary", dir, side,
@@ -96,8 +95,8 @@ void Setup::InitFlow(DataBlock &data) {
   // Create a host copy
   DataBlockHost d(data);
 
-  real r,th,rl;
-  real PonRho, vwind0, rc, vc;
+  real r,rl;
+  real PonRho, vwind0, rc;
   real cs=cs_vescGlob*sqrt(2.);
 
 
@@ -110,7 +109,7 @@ void Setup::InitFlow(DataBlock &data) {
     for(int j = 0; j < d.np_tot[JDIR] ; j++) {
       for(int i = 0; i < d.np_tot[IDIR] ; i++) {
         r=d.x[IDIR](i);
-        th=d.x[JDIR](j);
+
         real vwind;
 
         vwind = ParkerWind(r/rc) * cs;