First, you need to realize that Da uses symmetry-adapted AOs rather than AOs. Use Da_subset("AO") rather than Da to get atomic orbitals.
AOs are ordered as follows: loop over atoms, loop over quantum number l (angular), loop over quantum number n (principal quantum), loop over those by m_l in pattern 0, +1, -1, +2, -2, …
A quick but non-programmatic way to match an AO index to a physical AO is to use set print_mos true and look for the “Molecular Orbitals” section of the output file. The row labels will contain the information you’re looking for. For example, minimal basis water:
1 O1 s0 0.9967095 0.2233723 0.0000000 0.0932516 0.0000000
2 O1 s0 0.0153898 -0.8532247 0.0000000 -0.5183826 0.0000000
3 O1 p0 0.0029549 -0.1154667 0.0000000 0.7687640 0.0000000
4 O1 p+1 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000
5 O1 p-1 0.0000000 0.0000000 0.6088096 0.0000000 0.0000000
6 H2 s0 -0.0034491 -0.1510180 -0.4451103 0.2948036 0.0000000
7 H3 s0 -0.0034491 -0.1510180 0.4451103 0.2948036 0.0000000