I basically agree with the expression you propose however, as nemo use a staggered grid (arakawa c-grid), u_{rot} and v_{rot} are not located at the same place (so you cannot add them directly) and the angle \alpha will won’t be the same for u_{rot} and v_{rot}.
One solution could be to first compute u_{rot} and v_{rot} on a common T point
u_rot_T(ji,jj) = 0.5 * ( u_rot_U(ji,jj) + u_rot_U(ji-1,jj) )
v_rot_T(ji,jj) = 0.5 * ( v_rot_V(ji,jj) + v_rot_V(ji,jj-1) )
and next apply your expression
u_latlon_T = u_rot_T cos(alpha_T) - v_rot_T sin(alpha_T)
v_latlon_T = u_rot_T sin(alpha_T) + v_rot_T cos(alpha_T)
The angle cos(alpha_T) and sin(alpha_T) are computed in the subroutine angle in the file src/OCE/SBC/geo2oce.F90. It is variables gsint and gcost.
The vector rotation itself is coded in the subroutine rot_rep in the same file src/OCE/SBC/geo2oce.F90. What you want to do correspond to the cases 'ij->e’ and 'ij->n’.
You can find example of the use of rot_rep in the subroutine sbc_cpl_snd in the file src/OCE/SBC/sbccpl.F90
I hope this helps,
Sébastien