Scheduled Downtime
On Friday 21 April 2023 @ 5pm MT, this website will be down for maintenance and expected to return online the morning of 24 April 2023 at the latest

help with ungrib.exe when applied to ecmwf data at T1279 resolution

S

satjemkes

Member
Although i have posted a replz to the calvosa post, this might not receive required attention.

my problem is that i can not find a solution to what i believe must be solvable.
I retrieve data from ecmwf in two different spectral truncations T799 and T1279. The T799 can be processed, the T1279 not by the ungrib routine.
This refers to the surface data. atmospheric data seems to be fine

There is a post indicating to prescribe the record length. This is not working as the data originates from ECMWF mars directly. That means that in the
rd_grib1.F routine there is a switch which indicates that if the (map%nx.eq.65535) (which is true) then only data from the NCEP grid IDs is 37 to 44 can be processed.
This is data from ECMWF (id 98) with grid id 255. So it will not pass the if statement. (around line 530 in rd_grib1.F)

Is there a fix? is there a user who succcessfully uses ECMWF data at T1279 with ungrib.

A solution would be to use the data at ncar, however, as i want to generate an ensemble, i would like to use the output from ecmwf ensemble system and it is not clear if this lives outside ecmwf mars. Or i could use forecast data )as opposed to analysis) again not clear to me if ncar has this either.
If there is a user who has solved this problem with the data from ecmwf and is willing to write a line with a solution, this will be much appreciated as right now we are stuck.
 
Hi,
the solution to the problem is not to use global ECMWF data. I have downloaded the T1279 data for a region which is a little bigger than my target d01 domain, and then it runs smoothly without any modifications.
 
Would you please clarify the following issues:
(1) where did you download the T1279 data ?
(2) did you use Vtable.ECMWF to ungrib both surface and upper air data? Is there any issue with the landsea mask?
(3) T1279 data is model-level data or pressur level data?
Thanks !
 
Dear Ming,
See my other post as well.
T1279 is from ECMWF MARS directly at the 137 model level (which are hybrid sigma).
I did not use the Vtable. ECMWF but i have a special Vtable for the surface and atmosphere. Do not know what are the differences. atmosphere is grib2 surface is grib1, so you need two separate runs with ungrib as you know.
i have not paid any attention to the landsea mask issue reported by my colleagues. Not even sure i request from MARS the lsmask. i believe is that not provided by the geos ancilliary data provided by WPS? this is at 10 m resolution i believe

I could try to get the data at 91 model levels or at fixed pressure levels. But as i indicated in the other post. T1279 surface data can not be processed by current ungrib routine even after providing the record length information in the namelist.

This is a real pity as ecmwf data at 9x9 km resolution would be nice to use for the WRF.


stephen
 
Allow me to summarize my current situation.
Currently i am using ECMWF from IFS at 799 truncation, all 137 levels. The atmospheric component can not be processed by the ungrib routine. For this one need to use the pywinter module. This is very useful tool and i would like to thank mr. Suarez for making this available.
the surface component can also be processed with pywinter, but i prefer to use the ungrib routine.

It took quite a while to figure this out. It would be appreciate if in the user guide ponts the users to the pywinter application.

I have not tried to run the system with the T1279 truncated ECMWF data.
 
Thank you for the update. We are aware of the ungrib issue when it is employed to process the high-resolution IFS data at 137 levels. We haven't used the pywinter package to process the data. It will be great if you can provide a brief description how to process IFS data using pywinter. However, I completely understand if you cannot for some reasons. Whatsoever, thank you very much for the kind information.
 
I have no problems sharing code, with the notion that i am not a programmer. The code works but i will not get the Nobel price for clean - consistent python coding. it is all in a python module which is attached. Some explanation is below. the relevant attribute here is __ungrib_oper_set__ or the __ungrib_elda_set__ with the difference that the elda (=ensemble) data is used as initial condition and the oper (=regular deterministif forecast) used as lbc.

also i have not tried to run the T1279 spectral truncated data. It should be possible with pywinter (without the ungrib for the surface component as i do now) but have no time to figure this out now.

the way i do it is a combination of the ungrib and pywinter. This is causing confusion as the output filename convention by ungrib is not consistent. And this is a problem for the calculation of the pressure grid as the wps code expects the intermediate surface file name (which has the surface pressure) in a certain format

To be more specific: i use ungrib for the surface component. I could use pywinter for this as well, but have not figured out exactly how to handle the multiple soil layers. So i use ungrib for this. The problem is that sometimes ungrib generates an output file with the :00:00 (surface_file = datetime.datetime.strftime(validity_time, "FILE_SFC:%Y-%m-%d_%H:00:00")) and sometimes not. I have no time to dig through the code to understand why.
So i make sure that both files exists (dst = datetime.datetime.strftime(validity_time, "FILE_SFC:%Y-%m-%d_%H").

pywinter allows me to have more control over the output file of the atmosphere part.

step 1: take the atmosphere.grib file from ECMWF IFS and convert to netcdf using the eccodes

cmd = "grib_to_netcdf {} -o {}".format(grib_file, atmosphere_netcdf_file)

this cmd is executed in a ksh script

then read the netcdf file into the main python script and write the intermediate file using pywinter

#[ [4.0] Convert the atmosphere file
ec_atmos_data = Dataset(os.path.join(ungrib_work_directory, atmosphere_netcdf_file), 'r')

lat = ec_atmos_data.variables['latitude'][:] # degrees north
lon = ec_atmos_data.variables['longitude'][:] # degrees east

dlat = numpy.abs(lat[1] - lat[0])
dlon = numpy.abs(lon[1] - lon[0])
nlats = ec_atmos_data.variables['u'].shape[2]//2
# create winter Geo-information for cylindrical equidistant projection
winter_geo = pyw.Geo4(lat[0], lon[0], nlats, dlon, False)

winter_u = pyw.V3d('UU', ec_atmos_data.variables['u'][0, ::-1, :, :])
winter_v = pyw.V3d('VV', ec_atmos_data.variables['v'][0, ::-1, :, :])
winter_t = pyw.V3d('TT', ec_atmos_data.variables['t'][0, ::-1, :, :])
winter_q = pyw.V3d('SPECHUMD', ec_atmos_data.variables['q'][0, ::-1, :, :])
print("Shape info u {}".format(ec_atmos_data.variables['u'].shape))
print("Shape info v {}".format(ec_atmos_data.variables['v'].shape))
print("Shape info t {}".format(ec_atmos_data.variables['t'].shape))
print("Shape info q {}".format(ec_atmos_data.variables['q'].shape))

total_fields = [winter_t,
winter_q,
winter_u,
winter_v]

# Write intermediate file
winter_date = datetime.datetime.strftime(validity_time, "%Y-%m-%d_%H")
pyw.cinter('FILE_ML', winter_date, winter_geo, total_fields, ungrib_work_directory)

=================================

Calculate the pressure coordinates of the 137 level atmospheric file
==============================================
# [5.0]
# we need the a / b coefficients for the conversion sigma => p
src = os.path.join(
self.configuration['file system']["maatdaf.support"],
self.configuration['file system']["code repository"]["wps"], "ecmwf_coeffs")
dst = os.path.join(ungrib_work_directory, "ecmwf_coeffs")
if os.path.isfile(dst):
os.remove(dst)
shutil.copy(src, dst)
script_file = "calc_ecmwf_p.ksh"
# flags
self.execution_mode = 0o744
fp = open(script_file, 'w')
fp.write("#! /bin/ksh\n")
cmd = "source ~/.muprc"
fp.write("{}\n".format(cmd))
fp.write("ln -sf {} {}\n".format(
os.path.join(
self.configuration['file system']["maatdaf.support"],
self.configuration['file system']["code repository"]["wps"], "util/calc_ecmwf_p.exe"),
"calc_ecmwf_p.exe")
)
fp.write("./{}\n".format("calc_ecmwf_p.exe")) # code partition
fp.close()
os.chmod(script_file, self.execution_mode)
cmd = ["ksh", "./{}".format(script_file)]
self.logger.debug("issued {}".format(cmd))
subprocess.call(cmd)
===========================================

Finally combine all into the d1 file
==========================================
#
# the metgrid combines the pressure, surface and atmosphere file into a format suitable for wrf
# output is a met_em.d01 file with %Y-%m-%d_%H:%M:00.nc format
script_file = "metgrid.ksh"
# flags
self.execution_mode = 0o744
fp = open(script_file, 'w')
fp.write("#! /bin/ksh\n")
cmd = "source ~/.muprc"
fp.write("{}\n".format(cmd))
fp.write("ln -sf {} {}\n".format(
os.path.join(
self.configuration['file system']["maatdaf.support"],
self.configuration['file system']["code repository"]["wps"], "metgrid.exe"),
"metgrid.exe")
)
fp.write("./{}\n".format("metgrid.exe")) # code partition
fp.close()
os.chmod(script_file, self.execution_mode)
cmd = ["ksh", "./{}".format(script_file)]
self.logger.debug("issued {}".format(cmd))
subprocess.call(cmd)

=====================================================
 
You must log in or register to reply here.
Top