CDO

Sorry. The command is like this:

cdo setzaxis,myzaxisinvert.dat infile.nc temp1.nc
cdo ml2pl,100000,97500,95000,92500,90000,87500,85000,82500,80000,77500,75000,70000,67500,65000,62500,60000,57500,55000,52500,50000,47500,45000,42500,40000,37500,35000,32500,30000,27500,25000,22500,20000,17500,15000,12500,10000,7000,5000,3000 temp1.nc temp2.nc

Use setmisstoc,1.e10 to convert the missing values to the value 1.e10:

cdo -setmisstoc,1.e10 -ml2pl,100000,97500,.... temp1.nc temp2.nc

Hi Uwe,

Thanks for your quick reply. I just tried your method, but the output netcdf file still have this "-," as missing values. Can I upload the file for you? It's 3.5GB for a single file.
Thanks again.

Yes, please upload an example file. Try to reduce the size by selecting only a few timesteps. Here is an example to select the first 3 timesteps with the NCO tool ncks:

ncks -d time,0,2 ifile ofile

Hi Uwe,

Here is the new netcdf file and my z-axis file.

I need to do this:
cdo setzaxis,myzaxisinvert.dat output_for_wrf_199901.for.cdo temp1.nc
cdo ml2pl,100000,97500,95000,92500,90000,87500,85000,82500,80000,77500,75000,70000,67500,65000,62500,60000,57500,55000,52500,50000,47500,45000,42500,40000,37500,35000,32500,30000,27500,25000,22500,20000,17500,15000,12500,10000,7000,5000,3000 temp1.nc temp2.nc

You see I didn't open the "EXTRAPOLATE" option for "ml2pl", because when I open that, there will be interpolated values for Geopotential Height over the mountains regions, like west mountain region of US and Himalaya, where we want the values to be blank or infinite positive (1.oe+10). I go to the source code and didn't find how cdo extrapolate these values and so I think it's better for me to keep them blank.

Thanks for your help!

Hi Can you see the netcdf files? I have used different browers to upload them but I didn't the attachment. The file is only 25mb.

Try to upload the file again.

Thanks for the example file!
The missing value in this file is already 1.0e+10. The CDO operator setmisstoc can't be used to change a missing value to the same value as a constant. You can either use a different value:

cdo setmisstoc,1.e20 -ml2pl,100000,97500,.... temp1.nc temp2.nc

cdo setmisstoc,1.e10 -setmissval,-9e33 -ml2pl,100000,97500,.... temp1.nc temp2.nc

Hi Uwe,

In the example netcdf file, there are no missing values. I only defined the "missing_value" attribute in creating the file.

So with my example file,
after applying these two command:

cdo setzaxis,myzaxisinvert.dat output_for_wrf_199901.for.cdo temp1.nc
cdo -setmisstoc,1.e10 -ml2pl,100000,97500,.... temp1.nc temp2.nc

The new temp2.nc files will contain lots of "-," for missing values.

I have uploaded my result for you.
Thanks for your time.

Yes your original file has not missing value. But the missing_value attribute is used for the resulting missing values of ml2pl. Please try it with:

cdo setmisstoc,1.e10 -setmissval,-9e33 -ml2pl,100000,97500,.... temp1.nc temp2.nc

Hi Uwe,

Great! It works this time!
Thanks for your time and patience!
Appreciate!
BTW, do you know where can I find document about what's the method CDO was using when performing "EXTRAPOLATE" with "ml2pl". With this option, I found that the missing values for the Geopotential Height is filled with surface Height. So I want to make sure about that.

Thanks again!

Unfortunately the extrapolation method is not documented. The code used in CDO is from our ECHAM postprocessor called afterburner (https://code.zmaw.de/projects/afterburner/wiki). This is a very old tool and came without documentation.
What I saw from the code is that the extrapolated values depends on the pressure, surface height and temperature.
Here is the code:

double extra_Z(double pres, double halfp, double fullp, double geop, double temp)
{
  double alpha, tstar, tmsl, zalp, zalpal;
  double zrg;
  double zlapse = 0.0065;
  double ztlim = 290.5;

  zrg   = 1.0 / Grav;
  alpha = RD * zlapse * zrg;
  tstar = (1.0 + alpha * (halfp/fullp - 1.0)) * temp;

  if ( tstar < 255.0 ) tstar = 0.5 * (255.0 + tstar);

  tmsl = tstar + zlapse * zrg * geop;

  if ( tmsl > ztlim && tstar > ztlim )
    {
      tstar = 0.5 * (ztlim + tstar);
      tmsl  = tstar;
    }

  if ( tmsl > ztlim && tstar <= ztlim ) tmsl = ztlim;

  if ( tmsl-tstar < 0.000001 && tstar-tmsl < 0.000001 )
    alpha = 0.0;
  else if ( geop > 0.0001 || geop < -0.0001 )
    alpha = RD * (tmsl-tstar) / geop;

  zalp   = log(pres/halfp);
  zalpal = zalp * alpha;

  return ((geop - RD*tstar*zalp*(1.0 + zalpal*(0.5 + zalpal/6.0)))*zrg);
}