I try to cf.read a .nc file, in which there are a variabe BC(time,lev,lat,lon), lev is atmosphere hybrid sigma pressure coordinate, and associated variables for converting vertical coordiante.

When I ncdump the .nc file, BC and parameters are stored separately (see attchment 1). However, once cf.read, all these parameters become domain ancils to BC of cf.field (see attachment 2).

My questions:

1. Is there any function to automatically convert vertical coodinate to pressure coodinate?
2. If not, how can I extract parameters, so that I can convert it myself? I've tried BC.get_domain_ancillary("long_name:vertical coordinate formula term: ap(k)"). However, I've told that 'Field' object has no attribute 'get_domain_ancillary'.

Reason for asking this question is to plot the vertical profile of the variable.

Many Thanks,
Liang

Attachment 1. ncdump -h N2000_AEROSLO_ensmean.2050-2099_BC.nc

netcdf N2000_AEROSLO_ensmean.2050-2099_BC {
dimensions:
time = UNLIMITED ; // (600 currently)
lev = 26 ;
lat = 96 ;
lon = 144 ;
bnds = 2 ;
ilev = 27 ;
variables:
float BC(time, lev, lat, lon) ;
BC:long_name = "Concentration" ;
BC:units = "kg m-3" ;
BC:cell_methods = "time: mean" ;
double P0 ;
P0:long_name = "reference pressure" ;
P0:units = "Pa" ;
float PS(time, lat, lon) ;
PS:long_name = "Surface pressure" ;
PS:units = "Pa" ;
PS:cell_methods = "time: mean" ;
double a(lev) ;
a:long_name = "vertical coordinate formula term: ap(k)" ;
a:units = "Pa" ;
double a_bnds(lev, bnds) ;
a_bnds:long_name = "vertical coordinate formula term: ap(k+1/2)" ;
a_bnds:units = "Pa" ;
double b(lev) ;
b:long_name = "vertical coordinate formula term: b(k)" ;
b:units = "1" ;
double b_bnds(lev, bnds) ;
b_bnds:long_name = "vertical coordinate formula term: b(k+1/2)" ;
b_bnds:units = "1" ;
double ilev(ilev) ;
ilev:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
ilev:axis = "Z" ;
ilev:positive = "down" ;
ilev:long_name = "hybrid sigma pressure coordinate" ;
ilev:units = "level" ;
ilev:formula_terms = "ap: ap b: b ps: ps" ;
ilev:bounds = "ilev_bnds" ;
double ilev_bnds(ilev, bnds) ;
ilev_bnds:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
ilev_bnds:units = "level" ;
ilev_bnds:formula_terms = "ap: ap_bnds b: b_bnds ps: ps" ;
double lat(lat) ;
lat:standard_name = "latitude" ;
lat:long_name = "latitude" ;
lat:units = "degrees_north" ;
lat:axis = "Y" ;
double lev(lev) ;
lev:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
lev:axis = "Z" ;
lev:positive = "down" ;
lev:long_name = "hybrid sigma pressure coordinate" ;
lev:units = "level" ;
lev:formula_terms = "a: a b: b p0: P0 ps: PS" ;
lev:bounds = "lev_bnds" ;
double lev_bnds(lev, bnds) ;
lev_bnds:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
lev_bnds:units = "level" ;
lev_bnds:formula_terms = "a: a_bnds b: b_bnds p0: P0 ps: PS" ;
double lon(lon) ;
lon:standard_name = "longitude" ;
lon:long_name = "longitude" ;
lon:units = "degrees_east" ;
lon:axis = "X" ;
double time(time) ;
time:standard_name = "time" ;
time:long_name = "time" ;
time:bounds = "time_bnds" ;
time:units = "days since 2000-01-01 00:00:00" ;
time:calendar = "365_day" ;
time:axis = "T" ;
double time_bnds(time, bnds) ;

Attachment 2. print(BC)

Field: long_name=Concentration (ncvar%BC)
Data                  : long_name=Concentration(time(600), atmosphere_hybrid_sigma_pressure_coordinate(26), latitude(96), longitude(144)) kg m-3
Cell methods
                 : time(600): mean
Dimension coords
                 : time(600) = [2050-02-01 00:00:00, ..., 2100-01-01 00:00:00] 365_day
                 : atmosphere_hybrid_sigma_pressure_coordinate(26) = [3.5446380000000097, ..., 992.5560999999998] level
                 : latitude(96) = [-90.0, ..., 90.0] degrees_north
                 : longitude(144) = [0.0, ..., 357.5] degrees_east
Coord references
                 : standard_name:atmosphere_hybrid_sigma_pressure_coordinate
Domain ancils
                 : long_name=vertical coordinate formula term: ap(k)(atmosphere_hybrid_sigma_pressure_coordinate(26)) = [0.0035446380000000097, ..., 0.0] Pa
                 : long_name=vertical coordinate formula term: b(k)(atmosphere_hybrid_sigma_pressure_coordinate(26)) = [0.0, ..., 0.9925560999999998] 1
                 : long_name=reference pressure() = 100000.0 Pa
                 : long_name=Surface pressure(time(600), latitude(96), longitude(144)) = [[[68659.515625, ..., 101580.4921875]]] Pa