anchor.diff:          rhoz% data = rhoz% data / rhoz% dx
           call elem_rhoz% append(rhoz% data, i, i*tau)
           rhoz% data = 0
  +        call compute_vxz
  +        call v_xz_el%append(v_xz, i, i*tau)
           call varia%tac()
   
           call sphere_io%position%append(colloids%pos)
  
positive.code:   type(profile_t) :: vx
    type(h5md_element_t) :: elem_tz, elem_tz_count
    type(h5md_element_t) :: elem_rhoz
  
    integer :: rho
    integer :: N
  
negative.code:   call tz% init(0.d0, solvent_cells% edges(3), L(3))
    call rhoz% init(0.d0, solvent_cells% edges(3), L(3))
  
    allocate(v_xz_count(L(3), L(1)))
    allocate(v_xz(2, L(3), L(1)))
  
    call h5gcreate_f(hfile%id, 'fields', fields_group, error)
    call vx_el%create_time(fields_group, 'vx', vx%data, ior(H5MD_LINEAR,H5MD_STORE_TIME), &
         step=N_MD_steps, time=N_MD_steps*dt)
    call elem_tz% create_time(fields_group, 'tz', tz% data, ior(H5MD_TIME, H5MD_STORE_TIME))
    call elem_tz_count% create_time(fields_group, 'tz_count', tz% count, ior(H5MD_TIME, H5MD_STORE_TIME))
    call elem_rhoz% create_time(fields_group, 'rhoz', rhoz% data, ior(H5MD_TIME, H5MD_STORE_TIME))
    call v_xz_el%create_time(fields_group, 'v_xz', v_xz, ior(H5MD_TIME, H5MD_STORE_TIME))
    call h5gclose_f(fields_group, error)
  
    call h5gcreate_f(sphere_io%group, 'box', box_group, error)
  