  Gravitational wave astronomy pipelines rely on template waveform models for
searches and parameter estimation purposes. For coalescing binary neutron stars
(BNS), such models need to accurately reproduce numerical relativity (NR) up to
merger, in order to provide robust estimate of the stars' equation of state -
dependent parameters. In this work we present an improved version of the
Effective One Body (EOB) model $\tt TEOBResumS$ for gravitational waves from
BNS systems. Building upon recent post-Newtonian calculations, we include
subleading order tidal terms in the waveform multipoles and EOB metric
potentials, as well as add up to 5.5PN terms in the gyro-gravitomagnetic
functions entering the spin-orbit sector of the model. In order to further
improve the EOB-NR agreement in the last few orbital cycles before merger, we
introduce next-to-quasicircular corrections in the waveform -- informed by a
large number of BNS NR simulations -- and introduce a new NR-informed parameter
entering the tidal sector of our conservative dynamics. The performance of our
model is then validated against 14 new eccentricity reduced simulations of
unequal mass, spinning binaries with varying equation of state. A time-domain
phasing analysis and mismatch computations demonstrate that the new model
overall improves over the previous version of $\tt TEOBResumS$. Finally, we
present a closed-form frequency domain representation of the (tidal) amplitude
and phase of the new model. This representation accounts for mass-ratio,
aligned spin and (resummed) spin-quadrupole effects in the tidal phase and --
within the calibration region -- it is faithful to the original model.
