Incorrect result units for some differential operators #807

When the axis of differentiation has units (equivalent to) radians, the differential operators curl_xy, div_xy, grad_xy, and laplacian_xy display an inapproriate radian-1 or radian-2 component to the returned Field units. E.g.

>>> import cf
>>> f = cf.example_field(0)
>>> f.grad_xy(radius='earth')  # WRONG
[<CF Field: long_name=X gradient of specific_humidity(latitude(5), longitude(8)) m-1.rad-1>,
 <CF Field: long_name=Y gradient of specific_humidity(latitude(5), longitude(8)) m-1.rad-1>]
>>> # SHOULD BE
#[<CF Field: long_name=X gradient of specific_humidity(latitude(5), longitude(8)) m-1>,
# <CF Field: long_name=Y gradient of specific_humidity(latitude(5), longitude(8)) m-1>]

Units of radian are dimensionless, and in the context of these operators it is misleading for them to appear in the output. For instance, the curl of the horizontal wind fields in spherical polar corodinates should have units of s-1 (i.e. relative vorticity), not "s-1.radian-1".

The agreed solution is to provide an option to Field.derivative that just pretends that the coordinates along the axis of differentiation are dimensionless (i.e. have units of 1) when instructed to do so with a keyword argument.