`cf.Field.subspace` - allow multiple subspacing conditions for the same axis

cf.Field.subspace has, to date, always disallowed setting two or more conditions for a single axis, e.g.

>>> import cf
>>> f = cf.example_field(0)
>>> f.subspace(Y=cf.lt(-10), latitude=cf.gt(20))
ValueError: Error: Can't specify 2 conditions for 1 axis: (<CF Query: (lt -10)>, <CF Query: (gt 20)>). Consider applying the conditions separately.

To get around this you could either apply the conditions separately:

>>> f.subspace(Y=cf.lt(-10)).subspace(Y=cf.gt(20))

or roll the conditions into one:

>>> f.subspace(Y=cf.lt(-10) & cf.gt(20))

There is, I think, is no need for this restriction though, and removing it will make subspacing Fields with a discrete axis (such as DSGs, UGRID, etc) a more intuitive process. For such fields, we typically have have multiple 1-d coordinates (with different identities) for the same axis, and we'd like to be able to do this:

>>> u = cf.example_field(8)
>>> print(u)
Field: air_temperature (ncvar%ta)
---------------------------------
Data            : air_temperature(time(2), ncdim%nMesh2_face(3)) K
Cell methods    : time(2): point (interval: 3600 s)
Dimension coords: time(2) = [2016-01-02 01:00:00, 2016-01-02 11:00:00] gregorian
Auxiliary coords: longitude(ncdim%nMesh2_face(3)) = [-44.0, -44.0, -42.0] degrees_east
                : latitude(ncdim%nMesh2_face(3)) = [34.0, 32.0, 34.0] degrees_north
Topologies      : cell:face(ncdim%nMesh2_face(3), 4) = [[2, ..., --]]
Connectivities  : connectivity:edge(ncdim%nMesh2_face(3), 5) = [[0, ..., --]]

>>> print(u.subspace(X=cf.gt(-43), Y=cf.wi(30, 40)))  # X and Y refer to the same discrete axis
Field: air_temperature (ncvar%ta)
---------------------------------
Data            : air_temperature(time(2), ncdim%nMesh2_face(1)) K
Cell methods    : time(2): point (interval: 3600 s)
Dimension coords: time(2) = [2016-01-02 01:00:00, 2016-01-02 11:00:00] gregorian
Auxiliary coords: longitude(ncdim%nMesh2_face(1)) = [-42.0] degrees_east
                : latitude(ncdim%nMesh2_face(1)) = [34.0] degrees_north
Topologies      : cell:face(ncdim%nMesh2_face(1), 4) = [[1, ..., --]]
Connectivities  : connectivity:edge(ncdim%nMesh2_face(1), 5) = [[2, ..., --]]

Implementing this will be "straight forward", as we essentially just have to introduce a loop over multiple conditions, AND-ing the selected indices from the application of each condition to its corresponding coordinates.

This change will be wholly backwards compatible - no existing code will change its behaviour, but new API patterns will become possible.