CDO

This topic: https://code.mpimet.mpg.de/boards/1/topics/167
provides relevant input netcdf file where cdo's ml2pl would work and I have tried replicating the netcdf structure but still unsuccessful.

ncatted -a bounds,lev,o,c,lev_bnds ta0_short1.nc ta3.nc # change attribute 
ncks -C -O -x -v lev_bnds,time_bnds ta3.nc ta4.nc # remove variable
ncks -C -O -x -v lev_bounds,time_bounds ta4.nc ta5.nc # remove variable
ncap2 -O -s 'lon_bnds=make_bounds(lon,$bnds);lat_bnds=make_bounds(lat,$bnds);time_bnds=make_bounds(time,$bnds);lev_bnds=make_bounds(lev,$bnds)' ta5.nc ta6.nc # add variable bounds
ncatted -O -a formula,lev_bnds,c,c,'p = ap + b*ps' -a units,ap,c,c,Pa -a units,ap_bnds,c,c,Pa -a formula_term,lev,d,, ta6.nc ta7.nc # missing/incorrect attributes based on good that that was properly handled by ml2pl
ncatted -O \
-a online_operation,ap,d,, \
-a coordinates,ap,d,, \
-a units,ap,d,, \
-a online_operation,ap_bnds,d,, \
-a coordinates,ap_bnds,d,, \
-a units,ap_bnds,d,, \
-a online_operation,b,d,, \
-a coordinates,b,d,, \
-a online_operation,b_bnds,d,, \
-a coordinates,b_bnds,d,, \
-a online_operation,ta,d,, \
-a interval_operation,ta,d,, \
-a interval_write,ta,d,, \
-a coordinates,ta,d,, \
-a description,ta,d,, \
-a bounds,lat,c,c,"lat_bnds" \
-a axis,lat_bnds,d,, \
-a long_name,lat_bnds,d,, \
-a standard_name,lat_bnds,d,, \
-a units,lat_bnds,d,, \
-a formula_terms,lev,c,c,"ap: ap b: b ps: ps" \
-a axis,lev,c,c,Z \
-a name,lev,d,, \
-a units,lev,d,, \
-a bounds,lev_bnds,d,, \
-a formula_term,lev_bnds,d,, \
-a name,lev_bnds,d,, \
-a positive,lev_bnds,d,, \
-a formula_terms,lev_bnds,c,c,"ap: ap_bnds b: b_bnds ps: ps" \
-a long_name,lev_bnds,o,c,"atmospheric model level bounds" \
-a units,lev_bnds,o,c," " \
-a bounds,lon,c,c,"lon_bnds" \
-a axis,lon_bnds,d,, \
-a long_name,lon_bnds,d,, \
-a standard_name,lon_bnds,d,, \
-a units,lon_bnds,d,, \
-a long_name,ps,o,c,"Surface Air Pressure" \
-a online_operation,ps,d,, \
-a cell_methods,ps,d,, \
-a interval_operation,ps,d,, \
-a interval_write,ps,d,, \
-a _FillValue,ps,d,, \
-a missing_value,ps,d,, \
-a coordinates,ps,d,, \
-a standard_name,ps,d,, \
-a description,ps,d,, \
-a history,ps,d,, \
-a cell_measures,ps,d,, \
-a bounds,time,o,c,"time_bnds" \
-a long_name,time,o,c,"time" \
-a axis,time_bnds,d,, \
-a calendar,time_bnds,d,, \
-a long_name,time_bnds,d,, \
-a standard_name,time_bnds,d,, \
-a time_origin,time_bnds,d,, \
-a units,time_bnds,d,, \
-a bounds,time_bnds,d,, \
ta7.nc ta8.nc
cdo ml2pl,100000,85000,40000,15000,1000,700,500,300,200,100 -selname,ta ta8.nc ta9.nc # convert model to pressure levels

The relevant input file ncdump -h after applying the code above but before applying the last command, ml2pl

netcdf ta8 {
dimensions:
    lon = 256 ;
    bnds = 2 ;
    lat = 128 ;
    time = UNLIMITED ; // (2 currently)
    lev = 91 ;
variables:
    double lev_bnds(lev, bnds) ;
        lev_bnds:formula = "p = ap + b*ps" ;
        lev_bnds:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
        lev_bnds:units = " " ;
        lev_bnds:long_name = "atmospheric model level bounds" ;
        lev_bnds:formula_terms = "ap: ap_bnds b: b_bnds ps: ps" ;
    double ap(lev) ;
        ap:long_name = "vertical coordinate formula term: ap(k)" ;
    double ap_bnds(bnds, lev) ;
        ap_bnds:long_name = "vertical coordinate formula term: ap(k+1/2)" ;
    double b(lev) ;
        b:long_name = "vertical coordinate formula term: b(k)" ;
    double b_bnds(bnds, lev) ;
        b_bnds:long_name = "vertical coordinate formula term: b(k+1/2)" ;
    double lev(lev) ;
        lev:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
        lev:long_name = "atmospheric model level" ;
        lev:positive = "down" ;
        lev:formula = "p = ap + b*ps" ;
        lev:bounds = "lev_bnds" ;
        lev:formula_terms = "ap: ap b: b ps: ps" ;
        lev:axis = "Z" ;

The relevant input file from the other topic:

netcdf tmp_sig_repair {
dimensions:
    lev = 26 ;
    bnds = 2 ;
    time = UNLIMITED ; // (1 currently)
    lat = 80 ;
    lon = 180 ;
variables:
    double ap(lev) ;
        ap:long_name = "vertical coordinate formula term: ap(k)" ;
    double ap_bnds(lev, bnds) ;
        ap_bnds:long_name = "vertical coordinate formula term: ap(k+1/2)" ;
    double b(lev) ;
        b:long_name = "vertical coordinate formula term: b(k)" ;
    double b_bnds(lev, bnds) ;
        b_bnds:long_name = "vertical coordinate formula term: b(k+1/2)" ;
    double lev(lev) ;
        lev:bounds = "lev_bnds" ;
        lev:axis = "Z" ;
        lev:positive = "down" ;
        lev:long_name = "atmospheric model level" ;
        lev:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
        lev:formula = "p = ap + b*ps" ;
        lev:formula_terms = "ap: ap b: b ps: ps" ;
    double lev_bnds(lev, bnds) ;
        lev_bnds:formula = "p = ap + b*ps" ;
        lev_bnds:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
        lev_bnds:units = "" ;
        lev_bnds:formula_terms = "ap: ap_bnds b: b_bnds ps: ps" ;
        lev_bnds:long_name = "atmospheric model level bounds" ;

Did you find a solution to this? I wonder if support team can help?? I have the same problem.

Andrew Williams wrote in RE: ml2pl bounds issue with CNRM-CM6-1:

Did you find a solution to this? I wonder if support team can help?? I have the same problem.

Not with CDO no.
You can have a look at a similar method from NCAR's geocat: https://geocat-comp.readthedocs.io/en/latest/user_api/generated/geocat.comp.interpolation.interp_hybrid_to_pressure.html#geocat.comp.interpolation.interp_hybrid_to_pressure
It isn't as efficient as CDO but does the work well.
I had trouble installing Geocat, but copying the relevant functions from the source code and slightly modifying it worked for me.

Thanks Waren!

Geocat looks like a great option