diff --git a/vlaplex.f90 b/vlaplex.f90
index a804b59..26dab14 100644
--- a/vlaplex.f90
+++ b/vlaplex.f90
@@ -47,7 +47,7 @@ program VlaPlEx
   real(dp):: time
   real(dp):: dr, dv, dt
   integer::  nr, nv, nt, nz
-  integer::   i, iz,  j,  t
+  integer::   i, iz,  j,  t, z_inj
   integer:: j0 ! First integer of positive velocity
 
   real(dp):: Temp_bc ! Temperature
@@ -241,18 +241,33 @@ program VlaPlEx
     time = t * dt + t0
     call boundaryConditions%get(time, n_bc, u_bc, Temp_bc, Zave_bc)
     call writeOutputBoundary(t, dt, n_bc, u_bc, Temp_bc, Zave_bc)
+    z_inj = minloc(abs(Z_list - T_to_Z(Temp_bc)),1)
+    Zave_bc = Z_list(z_inj)
     u_bc = sqrt(Zave_bc * Temp_bc)
-    do iz = 1, nz
-      ! f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc)
-      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)
 
-      ! Boundary conditions
-      ! r = r0, v>0
-      f_i_old(iz,1,j0:nv) = f0
-      f_i(iz,1,j0:nv) = f_i_old(iz,1,j0:nv)
-      T_i(iz,1) = Temp_bc
-    end do
+    ! f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc)
+    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)
+
+    !f_i_old(:,1,j0:nv) = 0
+    !f_i(:,1,j0:nv) = 0
+    !T_i(:,1) = 0
+
+    f_i_old(z_inj,1,j0:nv) = f0
+    f_i(z_inj,1,j0:nv) = f_i_old(z_inj,1,j0:nv)
+    T_i(z_inj,1) = Temp_bc
+
+    !do iz = 1, nz
+    !  ! f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc)
+    !  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)
+
+    !  ! Boundary conditions
+    !  ! r = r0, v>0
+    !  f_i_old(iz,1,j0:nv) = f0
+    !  f_i(iz,1,j0:nv) = f_i_old(iz,1,j0:nv)
+    !  T_i(iz,1) = Temp_bc
+    !end do
     T_e = Temp_bc
     print *, 'Time: ', time * t_ref
     print *, 'Temp_bc: ', Temp_bc