diff --git a/vlaplex.f90 b/vlaplex.f90
index 869ce81..6053586 100644
--- a/vlaplex.f90
+++ b/vlaplex.f90
@@ -1,33 +1,33 @@
-! module eos
-!   use constantParameters, only: dp
-!   implicit none
-!
-!   private
-!   public:: T_to_Z
-!
-!   contains
-!     pure function T_to_Z(T) result(Z)
-!       use constantParameters, only: eV_to_K
-!       use referenceValues,    only: Temp_ref
-!       implicit none
-!
-!       real(dp), intent(in):: T
-!       real(dp):: Z
-!
-!       ! Z = (Temp_ref * T / eV_to_K)**0.6
-!       ! Z = max(Z, 1.0_dp)
-!       ! Z = min(Z, 22.0_dp)
-!       Z = 12.0_dp
-!
-!     end function T_to_Z
-!
-! end module eos
+ module eos
+   use constantParameters, only: dp
+   implicit none
+
+   private
+   public:: T_to_Z
+
+   contains
+     pure function T_to_Z(T) result(Z)
+       use constantParameters, only: eV_to_K
+       use referenceValues,    only: Temp_ref
+       implicit none
+
+       real(dp), intent(in):: T
+       real(dp):: Z
+
+      Z = 22.5 * (Temp_ref * T / eV_to_K /100.0)**0.6
+       ! Z = max(Z, 1.0_dp)
+       ! Z = min(Z, 22.0_dp)
+       !Z = 12.0_dp
+
+     end function T_to_Z
+
+end module eos
 
 program VlaPlEx
   use constantParameters, only: dp, kb, qe, eps_0, ev_to_K, cm3_to_m3, PI
   use output
   use referenceValues
-  ! use eos,                only: T_to_Z
+  use eos,                only: T_to_Z
   use moduleTableBC
   use omp_lib
   implicit none
@@ -96,7 +96,7 @@ program VlaPlEx
 
   ! Set input parameters (remember these have to be in non-dimensional units)
   c_s      = sqrt(11.0_dp * gamma_i * 1.0_dp)
-  bc_file = 'bc.csv'
+  bc_file = 'bc_80ns_T60Z16.csv'
   call boundaryConditions%init(bc_file)
 
   ! Set domain boundaries (non-dimensional units)
@@ -105,6 +105,7 @@ program VlaPlEx
   dr =   1.0e-6_dp / L_ref
   nr = nint((rf - r0) / dr) + 1
   dr =      (rf - r0) / float(nr-1)
+  print *, '#R gridpoints: ', nr 
   allocate(r(1:nr))
   do i = 1, nr
     r(i) = dr * float(i-1) + r0
@@ -121,6 +122,7 @@ program VlaPlEx
   dv = 1.0e-1_dp
   nv = nint((vf - v0) / dv) + 1
   dv =      (vf - v0) / float(nv-1)
+  print *, '#V gridpoints: ', nv
   allocate(v(1:nv))
   do j = 1, nv
     v(j) = dv * float(j-1) + v0
@@ -135,9 +137,10 @@ program VlaPlEx
   t0 = 0.0_dp
   tf = 2.0e-7_dp / t_ref
   ! tf = 1.0e1_dp * (rf - r0) / c_s
-  dt = 1.0e-2_dp*dr/c_s
+  dt = 7.0e2_dp*dr/c_s
   nt = nint((tf - t0) / dt) 
   dt =      (tf - t0) / float(nt)
+  print *, '#timesteps: ', nt 
 
   everyOutput = nint(1.0e-9_dp/t_ref/dt)
   if (everyOutput  == 0) then
@@ -236,6 +239,12 @@ program VlaPlEx
     f0(j0:nv) = 1.0_dp / sqrt(PI*Temp_bc) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc)
     f0 = f0 * n_bc / (sum(f0)*dv)
     T_e = Temp_bc
+    print *, 'Time: ', time * t_ref
+    print *, 'Temp_bc: ', Temp_bc
+    print *, 'Zave_bc: ', Zave_bc
+    print *, 'TtoZ: ', T_to_Z(Temp_bc)
+    print *, '-------------------------'
+
 
     ! Boundary conditions
     ! r = r0, v>0