anchor.diff:         ix = modulo(floor(solvent%pos(1,i)/solvent_cells%a), L(1)) + 1
          iy = modulo(floor(solvent%pos(2,i)/solvent_cells%a), L(2)) + 1
          iz = modulo(floor(solvent%pos(3,i)/solvent_cells%a), L(3)) + 1
  +       v_xyz_count(iz, iy, ix) = v_xyz_count(iz, iy, ix) + 1
  +       v_xyz(:, iz, iy, ix) = v_xyz(:, iz, iy, ix) + solvent%vel(:, i)
          if ( (iy == L(2)/2) .or. (iy == 1+L(2)/2) ) then
             v_xz_count(iz, ix) = v_xz_count(iz, ix) + 1
             v_xz(1, iz, ix) = v_xz(1, iz, ix) + solvent%vel(1, i)
  
positive.code:   double precision :: tau, dt , T, alpha
    double precision :: skin, co_max, so_max
    integer :: N_MD_steps, N_loop
    integer :: N_therm
    integer :: n_extra_sorting
    double precision :: kin_e, temperature
  
negative.code:         rhoz% data = rhoz% data / rhoz% dx
          call elem_rhoz% append(rhoz% data, i, i*tau)
          rhoz% data = 0
          call compute_vxz_and_vxyz
          if (modulo(i-N_therm, vxz_interval)==0) then
             call div_vxz()
             call v_xz_el%append(v_xz, i, i*tau)
             v_xz_count = 0
             v_xz = 0
  