pressure data from model levels to heights

Hi Yasmine,

I've had a quick look at your script and data.

Is your data directly downloaded or something was done to it? Because all I see is zero values for all the fields at all times.

I have downloaded the t,q,z and lnsp for one date and time and made a small example that works. 

In order to try my code out, please download t and q for one date and time (just so you don't spend too much time reading a big file and filtering why getting things to work) and save it to file tq_ml.grib. Then download z and lnsp for ml=1 for the same date and time and save it to another file called zlnsp_ml.grib.

My code looks something like this:

import metview as mv

tq = mv.read('tq_ml.grib')

t_one = tq.select(shortName='t')
q_one = tq.select(shortName='q')

zlnsp = mv.read('zlnsp_ml.grib')

lnsp_one = zlnsp.select(shortName='lnsp')

z_one = zlnsp.select(shortName='z')

#create z on model levels
z_ml = mv.mvl_geopotential_on_ml(t_one,q_one,lnsp_one,z_one)

#calculate pressure on model levels
p = mv.unipressure(lnsp_one)

# calculate pressure on 100m height above ground
hlevs = [100]
p_blh_ground = mv.ml_to_hl(p, z_ml, z_one, hlevs, "ground", "linear")

# calculate pressure on 100m height above sea level
p_blh_sea = mv.ml_to_hl(p, z_ml, None, hlevs, "sea", "linear")

mv.write('p_blh_ground.grib',p_blh_ground)
mv.write('p_blh_sea.grib',p_blh_sea)

Please note that in your script you put hlevs=10 and wanted to calculate above sea level ("sea" option). This is very low so almost anywhere above the land it will have missing values.

Your code also had an error. If you interpolate above sea level then you don't need z at ml=1 (orography) in the ml_to_hl, just put None, as I did.

You need z at ml=1 for calculating above the ground. You can see both examples in my code.

I hope this helps!

Milana

Dear Milana,

Thank you for your answer. I didn't touch the data. Please find attached the scripts used to download t,q,lnsp and z.

In fact, I noticed two things:

- z in z_era5_hourly_201701.grb. has the same value (not in the points) in all the timesteps which is weird

- q (specific humidity) is 0 everywhere

I re-downloaded the data today, and the same pattern remains.

main_z_ml1_2017_hourly_around_theni.py

main_q_ml_2017_hourly_around_theni.py

Hi Yasmine,

I see now why I didn't see anything in your data! The area is so tiny!! Now when I zoom to this point in India, I see data!

z on ml=1 will always have the same value in all time stamps because that is the model orography and doesn't change. 

So, it looks like we discovered why you had NAN values. 

You are trying to interpolate to height 10m above sea level, but all your land is higher than that, so basically you are trying to get the data where it is below the ground.

I think in this case you maybe want to use 'ground' option?

Best regards,
Milana

Dear Milena,

Thank you so much for helping me with this.

I selected a small area as a mock-up because it takes forever to download era5-complete data.

Does it mean please that I need to expand my domain in order to not have the issue of q=0 everywhere?

Can you please run the script with my data because I have NAN for the output of mv.ml_to_h even with the option ground but I still bump into the same issue:
p_blh = mv.ml_to_hl(p, z, zs, hlevs, "ground", "linear")
t_blh = mv.ml_to_hl(t, z, zs, hlevs, "ground", "linear")

PLEASE HELP!

Hi Yasmine,

To answer your questions.

Does it mean please that I need to expand my domain in order to not have the issue of q=0 everywhere?

No, I saw 0s because I was using Metview and it showed 0s where there is no data.... And the area was small with global zoom that I didn't notice that there is anything there.. But look there is... I think the data is fine..

I did a test with your data and indeed got NAN output for 10m above ground. But for 12m I got some values. 

This is the comparison between surface pressure and pressure at 12m. So I am guessing that at 10m it would be lower and that is why the values are missing. But I'm no expert on this.
I guess this happens because the unipressure() function does the approximation of pressure on ML from lnsp so the values won't be so amazingly perfect...

Milana,

Thank you so much for your help. This is fantastic, you saved my day.

Two more questions please,

1- What tool are you using to visualize the data?

2- Any recommendation on the correct way to save the data with mv.write() because I didn't manage to store the results in a grib file yet.

Please have a quick glance atsrc.zip the script attached.

Hi Yasmine,

1 - I am using Metview to visualise data. It works great with GRIB data, especially from ECMWF. Here you can find many examples that can help you get started (Gallery)

2- I've had a look at your script, but other than you write to the file with the same name every time in the loop, I can't see much more wrong. Try making a simple script that will do the thing just for one date/time and saving and then build on it. But check often that it still saves.

You can make an empty fieldset:   f = mv.Fieldset()

And then append to it: f.append(out)

mv.write('10m_pblh_era5_hourly_20170101.grb', f)

I hope this helps!

Dear Milana,

Can't thank you enough for your help. Problems solved! Thank you!

No worries. Glad I could help 

Dear Milana,

I have one more question, please. My script works but I have a problem with the speed of the process as the mv.ml_to_hl takes about 250 seconds for each timestep. I'm working with hourly data for a whole year and it has been running since a few days now.

I was wondering if you had any tip to improve the efficiency of the process, I'm attaching the final version of the script to this message.

Thank you again!pressure_from_ml_to_hl_with_metview.py

Hi Yasmine,

250seconds seems a bit long. Is this for the whole process or just ml_to_hl?
Are you calculating it for the whole globe or just a small area?
The only thing I see that you can try here is to see if it is a bit faster if you select parameters before you call the function.
Instead of:
fs_in_i = fs_in.select(**my_params[i])
Something like:
t_i = fs_in.select(**my_params[i],shortName='t')
q_i = fs_in.select(**my_params[i],shortName='q')
lnsp = fs_in_i["lnsp"]
and then use those when you call mv_routine()

Breaking your huge yearly file into daily or monthly files would also help (you can use metview for this, I think, just filter days and save it before you run your big script)

Dear Milana,

Thank you again for your prompt and extremely useful answer. It takes 250 seconds only for ml_to_hl, so it takes even longer for the whole process.

I will follow your recommendation to make it faster. I was wondering if there was any recent update to metview to increase the speed?

Kind regards.

Hi Yasmine,

I don't think there was any recent update to increase the speed. A lot of time can take in slicing the big grib, but since you are already using the select() function, you are most up to date.

Good luck!

Dear Milana,

I used metview function to interpolate temperature from model to height levels. To do so, I had to download t,q for all the 137 model levels, and it took a very long time with era5 complete. I tried to do the same using only the 4 first levels to boost time efficiency and the metview function for the interpolation didn't work.

Can you please tell me whether metview vertical interpolations only work when taking all the model levels as an input? 

Dear Yasmine,

Unfortunately the metview function needs all the model levels as an input to work.

Hi Yasmine,

I read your conversation, I have a question about Lnsp(Logarithm of surface pressure) only changes with time and latitude and longitude, not with height. And it only exists at the topmost level of model（ml=1)?