NWP | ECMWF | IFS | Open data

Open Dataset

ECMWF now provides a much larger open dataset to the public, representing weather forecasts at a higher resolution and a reduction in some release times, after a first set of data was made available in early 2022.

The changes can be summarised as follows:

  • Medium-range data of the Integrated Forecasting System (IFS) are provided at a resolution of 0.25 x 0.25 degrees (28 x 28 km), compared to 0.4 x 0.4 degrees before, and some additional parameters will be added.
  • Medium-range data of the latest beta version of the Artificial Intelligence/Integrated Forecasting System (AIFS) are also available at a resolution of 0.25 x 0.25 degrees.
  • ECMWF’s contribution to the sub-seasonal to seasonal (S2S) dataset, which is part of a global scheme, is now available with a 48-hour delay rather than the previous three-week delay.
  • All parameters of the EU’s Copernicus Atmosphere Monitoring Service (CAMS) that used to a have a six-day embargo are now available in real time.

An overview of open data is given on the ECMWF website.

歐洲中期天氣預報中心 (ECMWF) 的 IFS 模式現在可以作為 0.25 度解析度的開放資料取得。從 2024 年 2 月 1 日 開始，IFS 開放資料分佈的解析度從 0.4 度提高到 0.25 度。這項改進顯著提高了數據精度，尤其是在沿海地區和山區。改進後的分辨率受到了廣泛歡迎，因為之前的 0.4 度分辨率對於許多應用來說太粗糙。

資料下載連結：

https://data.ecmwf.int/forecasts/

規格如下：

  • 解析度從 0.4° 提升至 0.25°
  • 0.25° x 0.25° lat/lon grid or any multiple thereof (global or sub-area)
  • On model (Reduced Gaussian) n320 grid
  • 地表變因：2 公尺氣溫、地表溫度、降水量、氣壓、10 公尺風速
  • 9 個大氣壓力層位提供溫度、濕度或風速數據
  • 10 天預報，時間間隔 3 小時
  • 15 天集合預報，包含 51 個集合成員
  • 每天更新 4 次 (00/06/12/18)
  • 最新測試版的人工智慧/綜合預測系統（AIFS）的中尺度資料也可提供0.25 x 0.25度的解析度
  • 歐洲中期天氣預報中心（ECMWF）對次季節到季節尺度（S2S）資料集的貢獻，作為全球計劃的一部分，現在可以在48小時延遲後獲取，而不是之前的三週延遲
    • 歐洲中期天氣預報中心(ECMWF)在由世界氣象研究計畫(WWRP)和世界氣候研究計畫(WCRP)組織的預測計畫中，對次季節到季節尺度天氣預報也進行了大大減少的延遲。 ECMWF是此計畫的眾多預報中心之一。我們的預報現在只有48小時的延遲，而不是以前的三週延遲。計劃是鼓勵所有次季節性到季節尺度資料提供者同意減少他們的貢獻延遲。更多資訊可以在次季節性到季節尺度存檔頁面上找到（https://confluence.ecmwf.int/display/S2S/S2S+archive）
  • 歐盟的哥白尼大氣監測服務（CAMS）的所有參數，原本需要六天的禁止發布期，現在都可以即時取得
  • 升級後的解析度已無縫整合到 Open-Meteo 天氣 API 中。現在您可以比較 18 個不同天氣模型的預報，評估預報的不確定性並選擇最符合您需求的模型。
  • 0.25° 解析度的 51 個集合成員可以透過 Ensemble API 存取。有關 ECMWF IFS 開放資料計劃的更多詳細信息，請參閱 ECMWF 網站。

憑藉 0.25° 的分辨率，新的 IFS 開放資料模型適用於 Google GraphCast 機器學習模型的實際操作。先前，GraphCast 依賴歐洲地球觀測計畫的 ERA5 資料集，該資料集存在 5 天的延遲，使其不適用於即時預報應用。開放資料 IFS 0.25° 的發布現在允許您每 6 小時執行 GraphCast 並建立機器學習預報。唯一需要注意的是，高效運行 GraphCast 需要一個大型 GPU，但與超級電腦相比，這已經非常划算了。

  • 這些將實現有限區域建模和機器學習初始化。有關即時預測的文件可以在開放資料頁面找到。
  • 最新測試版AIFS的中期數據也以0.25 x 0.25度的解析度提供。這是目前AIFS運作的解析度。有關AIFS的更多信息，包括訪問數據的鏈接，請參閱AIFS機器學習數據頁面（https://www.ecmwf.int/en/forecasts/dataset/aifs-machine-learning-data）。

Snapshot

  • ---Only few days Available--- !!!
  • File format
    • The files are in GRIB edition 2 format, except for trajectories which are in BUFR edition 4 format.
Grib Code Mars Abbreviation Long Name Available model runs
1025 oper Atmospheric model 00z, 12z
1026 scda Atmospheric model (short cutoff) 06z, 18z
1035 enfo Ensemble prediction system 00z, 06z, 12z, 18z
1045 wave Wave model 00z, 12z
1081 waef Wave ensemble forecast

from https://codes.ecmwf.int/grib/format/mars/stream/

Time available (UTC)
Time available (UTC)

ECMWF open data: real-time forecasts from IFS and AIFS (Feb 26, 2025)

Currently, ECMWF open real-time data are available from these locations:

File format

The files are in GRIB edition 2 format, except for trajectories which are in BUFR edition 4 format. We recommend using ecCodes version 2.38.0 or newer to manipulate the GRIB and BUFR files.

File-naming convention

The files are provided with the following naming convention:

1
[ROOT]/[yyyymmdd]/[HH]z/[model]/[resol]/[stream]/[yyyymmdd][HH]0000-[step][U]-[stream]-[type].[format]
  • [ROOT] is the top-level URL of one of the sites hosting the data. See the above for possible values.
  • [yyyymmdd] is the reference date of the forecasts (base date).
  • [HH] is the reference time of the forecasts. Values are 00, 06 , 12 and 18.
  • [model] is the production model (IFS or AIFS). Note: IFS and AIFS have different options, so please be review dataset pages to see what is available.
  • [resol] is the horizontal resolution of the data. Options include: 0p25
  • [stream] is the forecasting system that produces the data. Values are:
    • oper - high-resolution forecast, atmospheric fields
    • enfo - ensemble forecast, atmospheric fields (not applicable for AIFS model)
    • waef - ensemble forecast, ocean wave fields, (not applicable for AIFS model)
    • wave - wave model, (not applicable for AIFS model)
    • scda - short cut-off high-resolution forecast, atmospheric fields (also known as "high-frequency products") (not applicable for AIFS model)
    • scwv - short cut-off high-resolution forecast, ocean wave fields (also known as "high-frequency products") (not applicable for AIFS model) and
    • mmsf - multi-model seasonal forecasts fields from the ECMWF model only (not applicable for AIFS model).
  • [step] is the forecast time step expressed in units U
  • [U] is the unit used for the time step. Values are h for hours and m for month. The latter is only valid for seasonal forecasts (mmsf).
  • [type] is once of fc (forecast), ef (ensemble forecast), ep (ensemble probabilities) or tf (trajectory forecast for tropical cyclone tracks).
  • [format] is grib2 for all fields, and bufr for the trajectories.

The valid combinations of the above are:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
format=bufr, type=tf 
   HH=00/12 
  	  stream=enfo/oper, step=240h 
   HH=06/18 
      stream=enfo, step=144h 
      stream=scda, step=90h 
format=grib2 
   HH=00/12 
      stream=enfo/waef 
        type=ef, step=0h to 144h by 3h, 144h to 360h by 6h 
        type=ep, step=240h/360h 
      stream=oper, wave 
        type=fc, step=0h to 144h by 3h, 144h to 240h by 6h 
   HH=06/18 
      stream=enfo/waef 
        type=ef, step=0h to 144h by 3h 
      stream= scda /scwv 
        type=fc, step=0h to 90h by 3h 
   HH=00 
   	  stream=mmsf, type=fc, u=m, step=1m to 7m

via Herbie

Variables

Single level - forecast

Short Name ID Long Name Units Additional information
z 129 Geopotential m^2 s^-2
msl 151 Mean sea level pressure Pa
sdor 160 Standard deviation of sub-gridscale orography m
slor 163 Slope of sub-gridscale orography Numeric
10u 165 10 metre U wind component m s^-1
10v 166 10 metre V wind component m s^-1
2t 167 2 metre temperature K
2d 168 2 metre dewpoint temperature K
ssrd 169 Surface short-wave (solar) radiation downwards J m^-2
lsm 172 Land-sea mask (0-1)
strd 175 Surface long-wave (thermal) radiation downwards J m^-2
lcc 3073 Low cloud cover %
mcc 3074 Medium cloud cover %
hcc 3075 High cloud cover %
rowe 231002 Runoff water equivalent (surface plus subsurface) kg m^-2
cp 228143 Convective precipitation kg m^-2
sf 228144 Snowfall kg m^-2
tcc 228164 Total cloud cover %
tp 228228 Total precipitation kg m^-2
100u 228246 100 metre U wind component m s^-1
100v 228247 100 metre V wind component m s^-1

Pressure level - forecast

  • All parameters are available at levels 1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100, 50 hPa
SHORT NAME ID LONG NAME UNITS ADDITIONAL INFORMATION
z 129 Geopotential m^2 s^-2
t 130 Temperature K
u 131 U component of wind m s^-1
v 132 V component of wind m s^-1
q 133 Specific humidity kg kg^-1
w 135 Vertical velocity Pa s^-1

Soil level - forecast

  • All parameters are available at levels 1, 2
SHORT NAME ID LONG NAME UNITS ADDITIONAL INFORMATION
vsw 260199 Volumetric soil moisture m^3 m^-3
sot 260360 Soil temperature K

open-data variables

IFS data High-resolution products:

Steps:

  • For times 00z &12z: 0 to 144 by 3, 150 to 360 by 6.
  • For times 06z & 18z: 0 to 144 by 3.
  • Type of level data:

Surface (sfc), Surface Other level (sol) and Pressure Level (pl)

Pressure levels available (hPa): 1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100, 50

Only a single ID is used for Volumetric soil water and Soil Temperature with layer 1 having level=1, layer 2 having level=2 etc. Layers 1-4 are available.

Short name Long name ID Type of level
z Geopotential (step 0) 129 sfc
sdor Standard deviation of sub-gridscale orography (step 0) 160 sfc
slor Slope of sub-gridscale orography (step 0) 163 sfc
10u 10 metre U wind component 165 sfc
10v 10 metre V wind component 166 sfc
100u 100 metre U wind component 228246 sfc
100v 100 metre V wind component 228247 sfc
10fg Maximum 10 metre wind gust since previous post-processing - (appears as 10fg3 for steps 3 to 144) 49 / 228028 sfc
2t 2 metre temperature 167 sfc
2d 2 metre dewpoint temperature 168 sfc
msl Mean sea level pressure 151 sfc
mp2 Mean zero-crossing wave period 140221 sfc
mwd Mean wave direction 140230 sfc
mwp Mean wave period 140232 sfc
pp1d Peak wave period 140231 sfc
swh Significant height of combined wind waves and swell 140229 sfc
ro Runoff 205 sfc
tp Total Precipitation 228 sfc
sp Surface pressure 134 sfc
st Soil temperature 228139 sfc
tcwv Total column vertically-integrated water vapour 137 sfc
lsm Land Sea Mask 172 sfc
vsw Volumetric soil water (layers 1-4) 260199 sol
sot Soil temperature (layers 1-4) 260360 sol
mucape Most-unstable convective available potential energy 228235 sfc
asn Snow albedo 32 sfc
mn2t3 Minimum temperature at 2 metres in the last 3 hours 228027 sfc
mx2t3 Maximum temperature at 2 metres in the last 3 hours 228026 sfc
mn2t6 Minimum temperature at 2 metres in the last 6 hours 122 sfc
mx2t6 Maximum temperature at 2 metres in the last 6 hours 121 sfc
tprate Total precipitation rate 260048 sfc
ptype Precipitation type 260015 sfc
ttr Top net long-wave (thermal) radiation 179 sfc
str Surface net long-wave (thermal) radiation 177 sfc
ssr Surface net short-wave (solar) radiation 176 sfc
ssrd Surface net short-wave (solar) radiation downwards 169 sfc
strd Surface net long-wave (thermal) radiation downwards 175 sfc
nsss Time-integrated northward turbulent surface stress 181 sfc
ewss Time-integrated eastward turbulent surface stress 180 sfc
t20d Depth of 20C isotherm - currently unavailable (11/03/2024) sfc
sav300 Average salinity in the upper 300m - currently unavailable (11/03/2024) sfc
ocu Eastward sea water velocity - currently unavailable (11/03/2024) sfc
ocv Northward sea water velocity - currently unavailable (11/03/2024) sfc
sithick Sea ice thickness 262000 sfc
zos Sea surface height 262124 sfc
d Divergence 155 pl
gh Geopotential height 156 pl
q Specific humidity 133 pl
r Relative humidity 157 pl
t Temperature 130 pl
u U component of wind 131 pl
v V component of wind 132 pl
w Vertical velocity 135 pl
vo Vorticity (relative) 138 pl

rda.ucar: ECMWF IFS High-Resolution Operational Forecasts

From rda.ucar (Not use)

  • https://rda.ucar.edu/datasets/d113001/#
    • ECMWF has implemented a significant resolution upgrade and methodology for high-resolution forecasts (HRES) and ensemble forecasts (ENS) beginning January of 2016. HRES is now performed via a transform grid with a nominal grid point spacing of 9 kilometers (0.08 degrees), and is carried out with IFS (Integrated Forecast System).
    • ECMWF Operational 6-hourly atmospheric surface forecast
    • ECMWF Operational 6-hourly atmospheric isobaric analysis

Herbie

NWP | Herbie | Download NWP model output (grib2) | MPAS | WRF

Download script

  • wget [ROOT]/20240301/00z/ifs/0p25/oper/20240301060000-24h-oper-fc.grib2
  • wget [ROOT]/20240301/06z/ifs/0p25/scda/20240301060000-24h-scda-fc.grib2
  • wget [ROOT]/20240301/12z/ifs/0p25/wave/20240301060000-12h-wave-fc.grib2
  • wget [ROOT]/20240301/06z/ifs/0p25/scwv/20240301060000-18h-scwv-fc.grib2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#!/bin/sh
#################################################################
# https://data.ecmwf.int/forecasts/20241117/00z/ifs/0p25/oper/
# https://data.ecmwf.int/forecasts/20241117/00z/ifs/0p25/oper/20241117000000-114h-oper-fc.grib2
#################################################################

yyyy='2024'
mm='11'
dd='17'
hh='00'
yyyymm=${yyyy}${mm}
yyyymmddhh=${yyyy}${mm}${dd}${hh}

tmp_dir="${yyyy}${mm}/${yyyy}${mm}${dd}"
mkdir -p ${tmp_dir}
cd ${tmp_dir}
echo "Current ${tmp_dir}"

limit_cpu_core=15
count_cpu_core=1

for ifsfile in $(seq 0 3 6); do # Forecast time = 12, 24, 36, 48, 60, 72, 84, 96, 108, 130, 142, 154

    tmp_ifs_fname="${yyyymmddhh}0000-${ifsfile}h-oper-fc.grib2"
	echo "Now, ${tmp_ifs_fname} ... "
    
		echo "Downloading ... ${tmp_ifs_fname}"
		wget -c --no-check-certificate https://data.ecmwf.int/forecasts/${yyyy}${mm}${dd}/${hh}z/ifs/0p25/oper/${tmp_ifs_fname}
    
	res=$(( ${count_cpu_core} % ${limit_cpu_core} ))
	if [[ ${res} == 0 ]]; then
		echo "please wait ${tmp_ifs_fname}, ${count_cpu_core}, ${res}."
		wait
	fi
    
done

read_grib2.py

1
2
3
4
5
6
7
#!/bin/python
import pygrib

grib_file = "/home/wpsze/IFS/202411/20241117/20241117000000-0h-oper-fc.grib2"
grbs = pygrib.open(grib_file)
for grb in grbs:
    print(grb) # long list

Run as ICBC

I recognised that the Vtable.ECMWF you have mentioned was not suitable for the IFS forecast data, due to the data form of the IFS was grib2. But the Vtable.ECMWF was for grib.

https://github.com/wrf-model/WPS/blob/master/ungrib/Variable_Tables/Vtable.ECMWF

MPAS/WRF

ungrib

Solution:

  • Using Newer ECMWF Data | kwerner | Sep 28, 2024
    • The new EC grib2 data requires an additional step to convert from EC grib2 data to grib2 data WPS supports. The conversion can be done by using the convert command in the ecCodes library which is available from ECMWF. The conversion step follows the following format: 
      1
      2
      $ grib_set -r -w packingType=grid_ccsds -s packingType=grid_simple input.grib2 output.grib2
      $ grib_set -r -w packingType=grid_ccsds -s packingType=grid_simple 20250326000000-0h-oper-fc.grib2 ifs_20250326000000-0h-oper-fc.grib2
    • The updated Vtable.ECMWF in WPS v4.6 should work with the converted grib2 data.
    • I think this could be related to CCSDS compression. You can change the packing type by using eccodes for instance.
    • Hi Everyone, One of the problems with the EC data is that it's not all standardized like GFS data. Depending on where/how the data are obtained, they can be different. We've found that, for e.g., for one user, ungrib is unable to recognize the soil fields. A possible solution for this is to use GFS soil data, in addition to the EC data. We expect the output to still be reasonable, especially because the WRF LSMs are more like that used by GFS. We hope to have a better solution for this issue in the future. For additional information, see this thread in our GitHub code repository (I believe one of you is the one who was having the discussion with my colleague).
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
GRIB1| Level| From |  To  | metgrid  | metgrid  | metgrid                                  |GRIB2|GRIB2|GRIB2|GRIB2|
Param| Type |Level1|Level2| Name     | Units    | Description                              |Discp|Catgy|Param|Level|
-----+------+------+------+----------+----------+------------------------------------------+-----------------------+
 129 | 100  |   *  |      | GEOPT    | m2 s-2   |                                          |  0  |  0  |     | 100 |
 156 | 100  |   *  |      | HGT      | m        | Height                                   |  0  |  3  |  5  | 100 |  
 130 | 100  |   *  |      | TT       | K        | Temperature                              |  0  |  0  |  0  | 100 |
 131 | 100  |   *  |      | UU       | m s-1    | U                                        |  0  |  2  |  2  | 100 |
 132 | 100  |   *  |      | VV       | m s-1    | V                                        |  0  |  2  |  3  | 100 |
 157 | 100  |   *  |      | RH       | %        | Relative Humidity                        |  0  |  1  |  1  | 100 |
 165 |  1   |   0  |      | UU       | m s-1    | U                    At 10 m             |  0  |  2  |  2  | 103 | 
 166 |  1   |   0  |      | VV       | m s-1    | V                    At 10 m             |  0  |  2  |  3  | 103 | 
 167 |  1   |   0  |      | TT       | K        | Temperature          At  2 m             |  0  |  0  |  0  | 103 |
 168 |  1   |   0  |      | DEWPT    | K        |                                          |  0  |  0  |  6  | 103 |
     |  1   |   0  |      | RH       | %        | Relative Humidity    At  2 m             |  0  |  0  |     | 103 | 
 172 |  1   |   0  |      | LANDSEA  | 0/1 Flag | Land/Sea flag                            |  2  |  0  |  0  |  1  | 
 129 |  1   |   0  |      | SOILGEO  | m2 s-2   |                                          |  0  |  0  |     | 103 | 
 156 |  1   |   0  |      | SOILHGT  | m        | Terrain field of source analysis         |  0  |  0  |     | 106 | 
 134 |  1   |   0  |      | PSFC     | Pa       | Surface Pressure                         |  0  |  3  |  0  |  1  |
 151 |  1   |   0  |      | PMSL     | Pa       | Sea-level Pressure                       |  0  |  3  |  0  | 101 |
 235 |  1   |   0  |      | SKINTEMP | K        | Sea-Surface Temperature                  |  0  |  3  |     | 101 |
  31 |  1   |   0  |      | SEAICE   | 0/1 Flag | Sea-Ice-Flag                             |  0  |  3  |     | 101 |
  34 |  1   |   0  |      | SST      | K        | Sea-Surface Temperature                  |  0  |  3  |     | 101 |
 141 |  1   |   0  |      | SNOW_EC  | m        |                                          |  0  |  3  |     | 101 |
     |  1   |   0  |      | SNOW     | kg m-2   |Water Equivalent of Accumulated Snow Depth|  0  |  3  |     | 101 |
 139 | 112  |   0  |   7  | ST000007 | K        | T of 0-7 cm ground layer                 |  2  |  0  |  2  | 106 | 
 170 | 112  |   7  |  28  | ST007028 | K        | T of 7-28 cm ground layer                | 192 | 128 | 170 | 106 | 
 183 | 112  |  28  | 100  | ST028100 | K        | T of 28-100 cm ground layer              | 192 | 128 | 183 | 106 | 
 236 | 112  | 100  | 255  | ST100289 | K        | T of 100-289 cm ground layer             | 192 | 128 | 236 | 106 | 
  39 | 112  |   0  |   7  | SM000007 | fraction | Soil moisture of 0-7 cm ground layer     | 192 | 128 | 39  | 106 | 
  40 | 112  |   7  |  28  | SM007028 | fraction | Soil moisture of 7-28 cm ground layer    | 192 | 128 | 40  | 106 | 
  41 | 112  |  28  | 100  | SM028100 | fraction | Soil moisture of 28-100 cm ground layer  | 192 | 128 | 41  | 106 | 
  42 | 112  | 100  | 255  | SM100289 | fraction | Soil moisture of 100-289 cm ground layer | 192 | 128 | 42  | 106 | 
-----+------+------+------+----------+----------+------------------------------------------+-----+-----+-----+-----+
#
#  Grib codes are from Table 128 
#  http://old.ecmwf.int/publications/manuals/d/gribapi/param/filter=grib1/order=paramId/order_type=asc/p=1/table=128/
#  
#  snow depth is converted to the proper units in rrpr.F
#
#  Tested on NCAR/RDA ds113.0 dataset. http://rda.ucar.edu/datasets/ds113.0/
#  Note that for ds113.0 there is one surface data file per day and 4 pressure-level files per day.
  • Vtable error of ECMWF,pls update it | Oct 13, 2024
    • Surface data including soil data may come from grib 1 data. You may need to use eccodes to separate out grib 2 and grib 1 data first before running ungrib. We are not supporting the new enfo data yet. Please post further issue on the Forum.
    • they are loss soil data.
    • It looks like the grib code for soil data has been messed up somehow. For the time being, you can try to use GFS soil data together with the EC atmospheric data. If you need help doing so, try reading the WPS section of the User's Guide.
    • Please read carefully this section of the User's Guide for using combined datasets.
  • Change in soil variables from ECMWF Cycle 49r1 | Nov 14, 2024
    • On Tuesday (November 12), ECMWF released the 49r1 version of the operational IFS. Despite all the improvements, they modified the soil variables (temperature and humidity) which stopped the ECMWF Vtable that had been released in the 4.6.0 version of the WPS.
    • WPS still can't find the soil temperature and moisture variables.
    • I contacted the ECMWF team and they gave me this feedback. It is necessary to make a change within Ungrib, more specifically in the rd_grib2.F code so that it works with the new IFS soil variables.
    • [Jan 3, 2025] Please take a look at this document-Change in soil variables from ECMWF Cycle 49r1 , which fixes the issues related to the ECMWF 49r1 version of the operational IFS. ECMWF 49r1 is newly released and we haven't done sufficient tests to make sure it works fine for WRF. Please post any issues you may have, which will be helpful for us to implement this data for WRF modeling study. Thanks in advance.
    • [Jan 4, 2025] Please download the modified code "rd_grib2.F" from the website update era5/ifs soil field processing by jimbresch · Pull Request #266 · wrf-model/WPS, recompile WPS, then rerun ungrib.exe.
    • IFS operational model data in the ECMWF open data repository are only available for the last 4 days.
    • [Jan 15, 2025] As requested, I am providing the corrected rd_grib2F file for WPS 4.5 and my Vtable.ECMWF_GRIB2 file. For the new WPS 4.6 version, I have not yet prepared a correction for the IFS model.
    • [Jan 16, 2025] I run WPS using your rd_grib2.F and Vtable. All work fine. Thank you! The met-em file I generated only has 14 pressure levels in the vertical, --- is this what you got?
    • [Jan 16, 2025] Yes, the IFS model provided by ECMWF has only 14 pressure levels. I made a simulation with the operational IFS model at 137 hybrid levels, which we have access to at the Institute of Meteorology and Water Management. The results are similar to the results from the IFS model, which is provided as open data. This is good news.
    • [14 Jan, 2025] I have already tested this Vtable.ECMWF with the new operational IFS data and it doesn't work either. I reported this here.

recompile WPS, then rerun ungrib.exe

  • under WPS-4.6.0/ungrib/src,
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
*****************************************************************************!
! Subroutine RD_GRIB2                                                         !
!                                                                             !
! Purpose:                                                                    !
!    Read one record from the input GRIB2 file.  Based on the information in  !
!    the GRIB2 header and the user-defined Vtable, decide whether the field in!
!    the GRIB2 record is one to process or to skip.  If the field is one we   !
!    want to keep, extract the data from the GRIB2 record, and store the data !
!    in the ungrib memory structure.                                          !
!                                                                             !
! Argument list:                                                              !
!    Input:                                                                   !
!       junit   : "Unit Number" to open and read from.  Not really a Fortran  !
!                 unit number, since we do not do Fortran I/O for the GRIB2   !
!                 files.  Nor is it a UNIX File Descriptor returned from a C  !
!                 OPEN statement.  It is really just an array index to the    !
!                 array (IUARR) where the UNIX File Descriptor values are     !
!                 stored.                                                     !
!       gribflnm     : File name to open, if it is not already open.          !
!       debug_level  : Integer for various amounts of printout.               !
!       pmin         : Minimum pressure level (Pa) to process.                !
!                                                                             !
!    Output:                                                                  !
!                                                                             !
!       hdate        : The (up to)19-character date of the field to process.  !
!       grib_edition : Version of the gribfile (1 or 2)                       !
!       ireaderr     : Error flag: 0 - no error on read from GRIB2 file.      !
!                              1 - Hit the end of the GRIB2 file.             !
!                              2 - The file GRIBFLNM we tried to open does    !
!                                  not exist.                                 !
!                                                                             !
!                                                                             !
! Author: Paula McCaslin, NOAA/FSL,   Sept 2004                               !
! Code is based on code developed by Steve Gilbert NCEP & Kevin Manning NCAR  !
! Adapted for WPS: Jim Bresch, NCAR/MMM. Sept 2006                            !
!*****************************************************************************!
      
      SUBROUTINE rd_grib2(junit, gribflnm, hdate, 
     &  grib_edition, ireaderr, debug_level, pmin)

      use grib_mod
      use params
      use table          ! Included to define g2code
      use gridinfo       ! Included to define map%
      use storage_module ! Included sub put_storage
      use module_debug

      real, allocatable, dimension(:) :: hold_array
      parameter(msk1=32000,msk2=4000)
      character(len=1),allocatable,dimension(:) :: cgrib
      integer :: listsec0(3)
      integer :: listsec1(13)
      integer year, month, day, hour, minute, second, fcst
      character(len=*)  :: gribflnm
      character(len=*)  :: hdate
      character(len=8)  :: pabbrev
      character(len=20) :: labbrev
      character(len=80) :: tabbrev
      integer :: lskip, lgrib
      integer :: junit, itot, icount, iseek
      integer :: grib_edition
      integer :: i, j, ireaderr, ith , debug_level
      integer :: currlen
      logical :: unpack, expand
      type(gribfield) :: gfld
      ! For subroutine put_storage
      real :: level
      real :: scale_factor
      integer :: iplvl
      character (len=9) :: my_field
      character (len=8) :: tmp8
      ! For subroutine output
      integer , parameter :: maxlvl = 150
      real , dimension(maxlvl) :: plvl
      integer :: nlvl
      integer , dimension(maxlvl) :: level_array
      real :: glevel1, glevel2
      logical :: first = .true.
      real :: pmin

C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  SET ARGUMENTS

      unpack=.true.
      expand=.true.
      hdate = '0000-00-00_00:00:00'
      ierr=0
      itot=0
      icount=0
      iseek=0
      lskip=0
      lgrib=0
      currlen=0
      ith=1
      scale_factor = 1e6
      call mprintf(.true.,DEBUG,"Begin rd_grib2", newline=.true.)

!/* IOS Return Codes from BACIO:  */
!/*  0    All was well                                   */
!/* -1    Tried to open read only _and_ write only       */
!/* -2    Tried to read and write in the same call       */
!/* -3    Internal failure in name processing            */
!/* -4    Failure in opening file                        */
!/* -5    Tried to read on a write-only file             */
!/* -6    Failed in read to find the 'start' location    */
!/* -7    Tried to write to a read only file             */
!/* -8    Failed in write to find the 'start' location   */
!/* -9    Error in close                                 */
!/* -10   Read or wrote fewer data than requested        */

!if ireaderr =1 we have hit the end of a file. 
!if ireaderr =2 we have hit the end of all the files. 
 

      ! Open a byte-addressable file.
      CALL BAOPENR(junit,gribflnm,IOS)
      first = .true.
      if (ios.eq.0) then 
      VERSION: do

         ! Search opend file for the next GRIB2 messege (record).
         call skgb(junit,iseek,msk1,lskip,lgrib)

         ! Check for EOF, or problem
         if (lgrib.eq.0) then
            exit 
         endif

         ! Check size, if needed allocate more memory.
         if (lgrib.gt.currlen) then
            if (allocated(cgrib)) deallocate(cgrib)
            allocate(cgrib(lgrib),stat=is)
            !print *,'G2 allocate(cgrib(lgrib)) status: ',IS
            currlen=lgrib
         endif

         ! Read a given number of bytes from unblocked file.
         call baread(junit,lskip,lgrib,lengrib,cgrib)

         call mprintf ((lgrib.ne.lengrib),ERROR,
     &    "rd_grib2: IO Error. %i .ne. %i ", newline=.true.,
     &    i1=lgrib,i2=lengrib)

         iseek=lskip+lgrib
         icount=icount+1

         call mprintf (.true.,DEBUG,
     &     "G2 GRIB MESSAGE  %i starts at %i ", newline=.true.,
     &      i1=icount,i2=lskip+1)

         ! Unpack GRIB2 field
         call gb_info(cgrib,lengrib,listsec0,listsec1,
     &                numfields,numlocal,maxlocal,ierr)
         call mprintf((ierr.ne.0),ERROR,
     &     " ERROR querying GRIB2 message = %i",newline=.true.,i1=ierr)
         itot=itot+numfields

         grib_edition=listsec0(2)
         if (grib_edition.ne.2) then
              exit VERSION
         endif
         
         ! Additional print statments for developer.
!MGD         if ( debug_level .GT. 100 ) then
!MGD     print *,'G2 SECTION 0: ',(listsec0(j),j=1,3)
!MGD         print *,'G2 SECTION 1: ',(listsec1(j),j=1,13)
!MGD         print *,'G2 Contains ',numlocal,' Local Sections ',
!MGD     &           ' and ',numfields,' data fields.'
!MGD         endif


         ! ----
         ! Once per file fill in date, model and projection values.

         if (first) then 
           first = .false.

           ! Build the 19-character date string, based on GRIB2 header date
           ! and time information, including forecast time information:

           n=1
           call gf_getfld(cgrib,lengrib,n,.FALSE.,expand,gfld,ierr)


           year  =gfld%idsect(6)     !(FOUR-DIGIT) YEAR OF THE DATA
           month =gfld%idsect(7)     ! MONTH OF THE DATA
           day   =gfld%idsect(8)     ! DAY OF THE DATA
           hour  =gfld%idsect(9)     ! HOUR OF THE DATA
           minute=gfld%idsect(10)    ! MINUTE OF THE DATA
           second=gfld%idsect(11)    ! SECOND OF THE DATA

           fcst = 0

! Extract forecast time. Assume the first field's valid time is true for all fields.
! This doesn't have to be true, but ungrib is designed to decode one time-level at
! a time.

           if ( gfld%ipdtnum .ne. 8 ) then
	     if ( gfld%ipdtmpl(8) .eq. 1 ) then       ! time units are hours
	       fcst = gfld%ipdtmpl(9)
	     else if ( gfld%ipdtmpl(8) .eq. 0 ) then  ! minutes
	       fcst = gfld%ipdtmpl(9) / 60.
	     else if ( gfld%ipdtmpl(8) .eq. 2 ) then  ! days
	       fcst = gfld%ipdtmpl(9) * 24.
	     else
               call mprintf(.true.,ERROR,
     &           "Time unit in ipdtmpl(8), %i is not suported",
     &            newline=.true.,i1=gfld%ipdtmpl(8))
	     endif
	   else
!  pdt 4.8 data are time-averaged, accumulated, or min/max fields with the 
!  ending (valid) time provided. 
	     year  =gfld%ipdtmpl(16) 
	     month =gfld%ipdtmpl(17)
	     day   =gfld%ipdtmpl(18)
	     hour  =gfld%ipdtmpl(19)
	     minute=gfld%ipdtmpl(20) 
	     second=gfld%ipdtmpl(21) 
	     fcst = 0.
	   endif

           if ( gfld%idsect(5) .eq. 2 ) fcst = 0.
           ! Compute valid time.

           !print *, 'ymd',gfld%idsect(6),gfld%idsect(7),gfld%idsect(8)
           !print *, 'hhmm  ',gfld%idsect(9),gfld%idsect(10)
   
           call build_hdate(hdate,year,month,day,hour,minute,second)
           call mprintf(.true.,DEBUG,"G2 hdate = %s ", newline=.true.,
     &                  s1=hdate)
           call geth_newdate(hdate,hdate,3600*fcst)
           call mprintf(.true.,DEBUG,"G2 hdate (fcst?) = %s ",
     &                  newline=.true., s1=hdate)

           !--

           ! Indicator of the source (center) of the data.
           icenter = gfld%idsect(1)

           ! Indicator of model (or whatever) which generated the data.
           iprocess = gfld%ipdtmpl(5)


           if (icenter.eq.7) then
             if (iprocess.eq.81) then
               map%source = 'NCEP GFS Analysis'
             elseif (iprocess.eq.82) then
               map%source = 'NCEP GFS GDAS/FNL'
             elseif (iprocess.eq.83) then
               map%source = 'NCEP HRRR Model'
             elseif (iprocess.eq.84) then
               map%source = 'NCEP MESO NAM Model'
             elseif (iprocess.eq.89) then
               map%source = 'NCEP NMM '
             elseif (iprocess.eq.96) then
               map%source = 'NCEP GFS Model'
             elseif (iprocess.eq.86 .or. iprocess.eq.100) then
               map%source = 'NCEP RUC Model'    ! 60 km
             elseif (iprocess.eq.101) then
               map%source = 'NCEP RUC Model'    ! 40 km
             elseif (iprocess.eq.105) then
               if (year .gt. 2011) then
                 map%source = 'NCEP RAP Model'
               else
                 map%source = 'NCEP RUC Model'  ! 20 km
               endif
             elseif (iprocess.eq.107) then
               map%source = 'NCEP GEFS'
             elseif (iprocess.eq.109) then
               map%source = 'NCEP RTMA'
             elseif (iprocess.eq.140) then
               map%source = 'NCEP NARR'
             elseif (iprocess.eq.44) then
               map%source = 'NCEP SST Analysis'
             elseif (iprocess.eq.70) then
               map%source = 'GFDL Hurricane Model'
             elseif (iprocess.eq.80) then
               map%source = 'NCEP GFS Ensemble'
             elseif (iprocess.eq.107) then             ! renumbered as of 23 Feb 2010
               map%source = 'NCEP GFS Ensemble'
             elseif (iprocess.eq.111) then
               map%source = 'NCEP NMMB Model'
             elseif (iprocess.eq.112) then
               map%source = 'NCEP WRF-NMM Model'
             elseif (iprocess.eq.116) then
               map%source = 'NCEP WRF-ARW Model'
             elseif (iprocess.eq.129) then
               map%source = 'NCEP GODAS'
             elseif (iprocess.eq.197) then
               map%source = 'NCEP CDAS CFSV2'
             elseif (iprocess.eq.25) then
               map%source = 'NCEP SNOW COVER ANALYSIS'
             else
               map%source = 'unknown model from NCEP'
               call mprintf(.true.,STDOUT,
     &            "unknown model from NCEP %i ",newline=.true.,
     &            i1=iprocess)
               call mprintf(.true.,LOGFILE,
     &            "unknown model from NCEP %i ",newline=.true.,
     &            i1=iprocess)
             end if
           else if (icenter .eq. 57) then
             if (iprocess .eq. 87) then
               map%source = 'AFWA AGRMET'
             else
               map%source = 'AFWA'
             endif
           else if ( icenter .eq. 58 ) then
             map%source = 'US Navy FNOC'
           else if (icenter .eq. 59) then
             if (iprocess .eq. 125) then
               map%source = 'NOAA GSD Rapid Refresh Model'
             else if (iprocess .eq. 83) then
               map%source = 'NOAA GSD HRRR Model'
             else if (iprocess .eq. 105) then
               map%source = 'NOAA GSD'
             else 
               print *,'Unknown GSD source'
               stop
             endif
           else if (icenter .eq. 60) then
             map%source = 'NCAR'
           else if (icenter .eq. 98) then
             map%source = 'ECMWF'
           else if (icenter .eq. 34) then
             map%source = 'JMA'
           else if (icenter .eq. 74 .or. icenter .eq. 75 ) then
             map%source = 'UKMO'
           else if (icenter .eq. 78 .or. icenter .eq. 79 ) then
             map%source = 'DWD'
           else
             map%source = 'unknown model and orig center'
           end if
           if (icenter.eq.7) then
             if (iprocess.eq.81 .or. iprocess.eq.82 .or.
     &           iprocess.eq.96) then
!  pmin should not be 15 or 40 hPa for gfs files. Use the next highest level.
               if (pmin .eq. 1500.) then
                 pmin = 1000.
               else if (pmin .eq. 4000.) then
                 pmin = 3000.
               endif
             endif
           endif

           call mprintf(.true.,DEBUG,"G2 source is = %s ",
     &                  newline=.true., s1=map%source)

           !--

           ! Store information about the grid containing the data.
           ! This stuff gets stored in the MAP variable, as defined in 
           ! module GRIDINFO.

           map%startloc = 'SWCORNER'
           map%grid_wind = .true.

           if (gfld%igdtnum.eq.0) then ! Lat/Lon grid aka Cylindrical Equidistant
              map%igrid = 0
              map%nx = gfld%igdtmpl(8)
              map%ny = gfld%igdtmpl(9)
              map%dx = gfld%igdtmpl(17)
              map%dy = gfld%igdtmpl(18)
              map%lat1 = gfld%igdtmpl(12)
              map%lon1 = gfld%igdtmpl(13)
              write(tmp8,'(b8.8)') gfld%igdtmpl(14)
              if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
              map%r_earth = earth_radius (gfld%igdtmpl(1),
     &                         gfld%igdtmpl(2),gfld%igdtmpl(3))

! Fix for NCEP 1/12 degree grids (e.g. rtgsst)
              if (icenter .eq. 7 .and. map%dx .eq. 83000. .and. map%nx 
     &              .eq. 4320) then
                map%lat1 = 89958333.
                map%lon1 = 41667.
                map%dx = 83333.333 * sign(1.0,map%dx)
                map%dy = 83333.333 * sign(1.0,map%dy)
              endif

              if ((gfld%igdtmpl(10) .eq. 0).OR.
     &            (gfld%igdtmpl(10) .eq. 255)) THEN
          ! Scale lat/lon values to 0-180, default range is 1e6.
                map%lat1 = map%lat1/scale_factor
                map%lon1 = map%lon1/scale_factor
          ! Scale dx/dy values to degrees, default range is 1e6.
                map%dx = map%dx/scale_factor
                map%dy = map%dy/scale_factor
              else
          ! Basic angle and subdivisions are non-zero (not tested)
                map%lat1 = map%lat1 *
     &                     (gfld%igdtmpl(11)/gfld%igdtmpl(10))
                map%lon1 = map%lon1 *
     &                     (gfld%igdtmpl(11)/gfld%igdtmpl(10))
                map%dx = map%dx * 
     &                     (gfld%igdtmpl(11)/gfld%igdtmpl(10))
                map%dy = map%dy * 
     &                     (gfld%igdtmpl(11)/gfld%igdtmpl(10))
               call mprintf(.true.,STDOUT,"WARNING - Basic angle option 
     &has not been tested, continuing anyway")
               call mprintf(.true.,LOGFILE,"WARNING - Basic angle option 
     & has not been tested, continuing anyway")
              endif


! The following is needed for NCEP GFS, 0.5 degree output. The j-scan is in the -y direction.
! In WPS, the sign of dy indicates the direction of the scan.
              write(tmp8,'(b8.8)') gfld%igdtmpl(19)
              read(tmp8,'(1x,i1)') jscan
              if ( jscan .eq. 0 .and. map%dy .gt. 0. ) then
                map%dy = -1. * map%dy
              endif
!             if ( map%lat1 .gt. gfld%igdtmpl(15) .and. 
!    &               map%dy .gt. 0. ) then
!               map%dy = -1. * map%dy
!               write(6,*) 'Resetting map%dy for iprocess = ',iprocess
!             endif

           elseif (gfld%igdtnum.eq.10) then ! Mercator Grid.
              map%igrid = 1
              map%nx = gfld%igdtmpl(8)
              map%ny = gfld%igdtmpl(9)
              map%lov = 0.
              map%truelat1 = gfld%igdtmpl(13) / scale_factor
              map%truelat2 = 0.
              map%dx = gfld%igdtmpl(18) / scale_factor
              map%dy = gfld%igdtmpl(19) / scale_factor
              map%lat1 = gfld%igdtmpl(10) / scale_factor
              map%lon1 = gfld%igdtmpl(11) / scale_factor
              write(tmp8,'(b8.8)') gfld%igdtmpl(12)
              if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
              map%r_earth = earth_radius (gfld%igdtmpl(1),
     &                       gfld%igdtmpl(2),gfld%igdtmpl(3))

           elseif (gfld%igdtnum.eq.20) then ! Polar-Stereographic Grid.
              map%igrid = 5
              map%nx = gfld%igdtmpl(8)
              map%ny = gfld%igdtmpl(9)
              map%lov = gfld%igdtmpl(14) / scale_factor
              map%truelat1 = 60.
              map%truelat2 = 91.
              map%dx = gfld%igdtmpl(15) / scale_factor
              map%dy = gfld%igdtmpl(16) / scale_factor
              map%lat1 = gfld%igdtmpl(10) / scale_factor
              map%lon1 = gfld%igdtmpl(11) / scale_factor
              write(tmp8,'(b8.8)') gfld%igdtmpl(12)
              if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
              map%r_earth = earth_radius (gfld%igdtmpl(1),
     &                       gfld%igdtmpl(2),gfld%igdtmpl(3))

           elseif (gfld%igdtnum.eq.30) then ! Lambert Conformal Grid
              map%igrid = 3
              map%nx = gfld%igdtmpl(8)
              map%ny = gfld%igdtmpl(9)
              map%lov = gfld%igdtmpl(14) / scale_factor
              map%truelat1 = gfld%igdtmpl(19) / scale_factor
              map%truelat2 = gfld%igdtmpl(20) / scale_factor
              map%dx = gfld%igdtmpl(15) / scale_factor
              map%dy = gfld%igdtmpl(16) / scale_factor
              map%lat1 = gfld%igdtmpl(10) / scale_factor
              map%lon1 = gfld%igdtmpl(11) / scale_factor
              write(tmp8,'(b8.8)') gfld%igdtmpl(12)
              if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
              map%r_earth = earth_radius (gfld%igdtmpl(1),
     &                       gfld%igdtmpl(2),gfld%igdtmpl(3))

           elseif(gfld%igdtnum.eq.40) then ! Gaussian Grid (we will call it lat/lon)
              map%igrid = 4
              map%nx = gfld%igdtmpl(8)     ! Ni - # of points along a parallel
              map%ny = gfld%igdtmpl(9)     ! Nj - # of points along meridian
              map%dx = gfld%igdtmpl(17)    ! Di - i direction increment
              map%dy = gfld%igdtmpl(18)    ! N - # of parallels between pole and equator
              map%lat1 = gfld%igdtmpl(12)  ! La1 - lat of 1st grid point
              map%lon1 = gfld%igdtmpl(13)  ! Lo1 - lon of 1st grid point
              write(tmp8,'(b8.8)') gfld%igdtmpl(14)  ! resolution/component flag
              if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
              map%r_earth = earth_radius (gfld%igdtmpl(1),
     &                       gfld%igdtmpl(2),gfld%igdtmpl(3))

              ! Scale dx/dy values to degrees, default range is 1e6.
              if (map%dx.gt.10000) then 
                 map%dx = map%dx/scale_factor
              endif
              if (map%dy.gt.10000) then 
                 map%dy = (map%dy/scale_factor)*(-1)
              endif

              ! Fix for zonal shift in CFSR data, following a similar fix 
              ! for global lat-lon data in rd_grib1.F
              if ( ABS(map%nx * map%dx - 360.0) < 1.0 ) then
                 if (ABS(map%dx - (360./real(map%nx))) > 0.00001) then
                    write(0,*) 'CFSR fix: recomputing delta-longitude'
                    map%dx = 360./real(map%nx)
                 endif
              endif

              ! Scale lat/lon values to 0-180, default range is 1e6.
              if (abs(map%lat1).ge.scale_factor) then
                 map%lat1 = map%lat1/scale_factor
              endif
              if (map%lon1.ge.scale_factor) then 
                 map%lon1 = map%lon1/scale_factor
              endif
              if ( debug_level .gt. 2 ) then
              call mprintf(.true.,DEBUG,
     &     "Gaussian Grid: Dx,Dy,lat,lon,nlats %f %f %f %f %i ",
     &     newline=.true.,f1=map%dx,f2=map%dy,f3=map%lat1,f4=map%lon1,
     &     i1=nint(map%dy))
              end if
           elseif (gfld%igdtnum.eq.32769) then ! Arakawa Non-E Staggered grid.
              map%igrid = 6               
              map%nx = gfld%igdtmpl(8)        
              map%ny = gfld%igdtmpl(9)        
              map%dx = gfld%igdtmpl(17)        
              map%dy = gfld%igdtmpl(18)       
              map%lat1 = gfld%igdtmpl(12)        
              map%lon1 = gfld%igdtmpl(13)        
              map%lat0 = gfld%igdtmpl(15)        
              map%lon0 = gfld%igdtmpl(16)        
              map%r_earth = earth_radius (gfld%igdtmpl(1),
     &                         gfld%igdtmpl(2),gfld%igdtmpl(3))
              if ((gfld%igdtmpl(10) .eq. 0).OR.
     &            (gfld%igdtmpl(10) .eq. 255)) THEN
                map%lat1 = map%lat1/scale_factor
                map%lon1 = map%lon1/scale_factor
                map%lat0 = map%lat0/scale_factor
                map%lon0 = map%lon0/scale_factor
                map%dx = map%dx/scale_factor/1.e3
                map%dy = map%dy/scale_factor/1.e3
              else
          ! Basic angle and subdivisions are non-zero (not tested)
                map%lat1 = map%lat1 *
     &                     (gfld%igdtmpl(11)/gfld%igdtmpl(10))
                map%lon1 = map%lon1 *
     &                     (gfld%igdtmpl(11)/gfld%igdtmpl(10))
                map%dx = map%dx * 
     &                     (gfld%igdtmpl(11)/gfld%igdtmpl(10))
                map%dy = map%dy * 
     &                     (gfld%igdtmpl(11)/gfld%igdtmpl(10))
               call mprintf(.true.,STDOUT,"WARNING - Basic angle option 
     &has not been tested, continuing anyway")
               call mprintf(.true.,LOGFILE,"WARNING - Basic angle option 
     & has not been tested, continuing anyway")
              endif

           else
              call mprintf(.true.,STDOUT,"GRIB2 Unknown Projection: %i",
     &          newline=.true.,i1=gfld%igdtnum)
              call mprintf(.true.,STDOUT,
     &          "ungrib understands projections 0, 20, 30, and 40", 
     &          newline=.true.)
              call mprintf(.true.,LOGFILE,
     &          "GRIB2 Unknown Projection: %i",
     &          newline=.true.,i1=gfld%igdtnum)
              call mprintf(.true.,LOGFILE,
     &          "ungrib understands projections 0, 10, 20, 30, and 40", 
     &          newline=.true.)
       ! If the projection is not known, then it can't be processed by metgrid/plotfmt
                stop 'Stop in rd_grib2'
           endif

           call mprintf(.true.,DEBUG,"G2 igrid = %i ,  dx = %f ,  dy = %
     &f ", newline=.true., i1 = map%igrid, f1=map%dx, f2=map%dy)
         
           if (icenter.eq.7) then
             call ncep_grid_num (gfld%igdtnum)
           endif

           ! Deallocate arrays decoding GRIB2 record.
           call gf_free(gfld)

         endif   ! "first" if-block

         ! ----

         ! Continue to unpack GRIB2 field.
         NUM_FIELDS: do n = 1, numfields 
	   ! e.g. U and V would =2, otherwise its usually =1
           call gf_getfld(cgrib,lengrib,n,.FALSE.,expand,gfld,ierr)
           if (ierr.ne.0) then
             write(*,*) ' ERROR extracting field gf_getfld = ',ierr
             cycle
           endif

! The JMA GSM has two different grids in the same GRIB file, so we need
! to process the map info for each field separately. If any other centers do
! this, then processing will need to be added here, too.

            if (icenter .eq. 34 .and. gfld%igdtnum.eq.0) then
              map%nx = gfld%igdtmpl(8)
              map%ny = gfld%igdtmpl(9)
              map%dx = gfld%igdtmpl(17)
              map%dy = gfld%igdtmpl(18)
              ! Scale dx/dy values to degrees, default range is 1e6.
              if (map%dx.gt.10000) then
                 map%dx = map%dx/scale_factor
              endif
              if (map%dy.gt.10000) then
                 map%dy = map%dy/scale_factor
              endif
              write(tmp8,'(b8.8)') gfld%igdtmpl(19)
              read(tmp8,'(1x,i1)') jscan
              write(0,*) 'gfld%igdtmpl(19) = ',gfld%igdtmpl(19),
     &   ' jscan = ',jscan
              if ( jscan .eq. 0 .and. map%dy .gt. 0. ) then
                map%dy = -1. * map%dy
              endif
            endif             ! JMA spectral

! ------------------------------------
         ! Additional print information for developer.
         if ( debug_level .GT. 1000 ) then
!MGD           print *
!MGD           print *,'G2 FIELD ',n
!MGD           if (n==1) then
!MGD            print *,'G2 SECTION 0: ',gfld%discipline,gfld%version
!MGD            print *,'G2 SECTION 1: ',(gfld%idsect(j),j=1,gfld%idsectlen)
!MGD           endif
!MGD           if ( associated(gfld%local).AND.gfld%locallen.gt.0 ) then
!MGD              print *,'G2 SECTION 2: ',(gfld%local(j),j=1,gfld%locallen)
!MGD           endif
!MGD           print *,'G2 SECTION 3: ',gfld%griddef,gfld%ngrdpts,
!MGD     &                            gfld%numoct_opt,gfld%interp_opt,
!MGD     &                            gfld%igdtnum
!MGD           print *,'G2 GRID TEMPLATE 3.',gfld%igdtnum,': ',
!MGD     &            (gfld%igdtmpl(j),j=1,gfld%igdtlen)
!MGD           if ( gfld%num_opt .eq. 0 ) then
!MGD             print *,'G2 NO Section 3 List Defining No. of Data Points.'
!MGD           else
!MGD             print *,'G2 Section 3 Optional List: ',
!MGD     &                (gfld%list_opt(j),j=1,gfld%num_opt)
!MGD           endif
!MGD           print *,'G2 PRODUCT TEMPLATE 4.',gfld%ipdtnum,': ',
!MGD     &          (gfld%ipdtmpl(j),j=1,gfld%ipdtlen)

           pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),
     &                              gfld%ipdtmpl(2))
           !call prlevel(gfld%ipdtnum,gfld%ipdtmpl,labbrev)
           !call prvtime(gfld%ipdtnum,gfld%ipdtmpl,listsec1,tabbrev)
!MGD            print *,'G2 TEXT: ',pabbrev,trim(labbrev)," ",trim(tabbrev)

!MGD           if ( gfld%num_coord .eq. 0 ) then
!MGD             print *,'G2 NO Optional Vertical Coordinate List.'
!MGD           else
!MGD             print *,'G2 Section 4 Optional Coordinates: ',
!MGD     &             (gfld%coord_list(j),j=1,gfld%num_coord)
!MGD           endif
           if ( gfld%ibmap .ne. 255 ) then
              call mprintf(.true.,DEBUG, 
     &             'G2 Num. of Data Points = %i with BIT-MAP %i', 
     &             newline=.true., i1=gfld%ndpts, i2=gfld%ibmap)
           else
              call mprintf(.true.,DEBUG, 
     &             'G2 Num. of Data Points = %i NO BIT-MAP', 
     &             newline=.true., i1=gfld%ndpts)
           endif
!MGD           print *,'G2 DRS TEMPLATE 5.',gfld%idrtnum,': ',
!MGD     &          (gfld%idrtmpl(j),j=1,gfld%idrtlen)
         endif ! Additional Print information 
! ------------------------------------

!         do i = 1, maxvar
!           write(6,'(a10,4i8)') namvar(i),(g2code(j,i),j=1,4)
!         enddo
!MGD      if (debug_level .gt. 50) then
!MGD          write(6,*) 'looking for ',gfld%discipline,gfld%ipdtmpl(1),
!MGD     &       gfld%ipdtmpl(2),gfld%ipdtmpl(10)
!MGD          endif
           call mprintf(.true.,DEBUG,"G2 Searching the g2code array (Vta
     &ble) for this grib field %i %i %i %i %i %i ", newline=.true.,
     & i1 = gfld%discipline, i2 = gfld%ipdtmpl(1),
     & i3 = gfld%ipdtmpl(2), i4 = gfld%ipdtmpl(10),
     & i5 = gfld%ipdtmpl(12), i6 = gfld%ipdtnum )


         ! Test this data record against list of desired variables 
         ! found in Vtable.
         ! ----
         MATCH_LOOP: do i=1,maxvar ! Max variables found in Vtable,
                                   ! maxvar is defined in table.mod

           if ( gfld%discipline .eq. g2code(1,i) .and.   !Discipline 
     &          gfld%ipdtmpl(1) .eq. g2code(2,i) .and.   !Category
     &          gfld%ipdtmpl(2) .eq. g2code(3,i) .and.   !Parameter
     &          gfld%ipdtmpl(10) .eq. g2code(4,i) .and.  !Elevation
     &          gfld%ipdtnum    .eq. g2code(5,i)) then   !Template

            call gf_free(gfld)
            call gf_getfld(cgrib,lengrib,n,.TRUE.,expand,gfld,ierr)
            pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),
     &                               gfld%ipdtmpl(2))

              !my_field (e.g. RH, TMP, similar to, but not the same as pabbrev)
              my_field=namvar(i) 

!MGD    if (debug_level .gt. 50) then
!MGD     write(6,*) 'namvar(i) = ',namvar(i),' pabbrev = ',pabbrev
!MGD     write(6,*) 'Parameter = ',gfld%ipdtmpl(2)
!MGD    endif
!  The following if-block is commented out since equivalent info can be obtained from g2print
!       if (debug_level .gt. 1000) then
!          fldmax=gfld%fld(1)
!          fldmin=gfld%fld(1)
!          sum=gfld%fld(1)
!          do j=2,gfld%ndpts
!            if (gfld%fld(j).gt.fldmax) fldmax=gfld%fld(j)
!            if (gfld%fld(j).lt.fldmin) fldmin=gfld%fld(j)
!            sum=sum+gfld%fld(j)
!          enddo ! gfld%ndpts
!          call mprintf(.true.,DEBUG,'G2 FIELD=%s MIN=%f AVG=%f MAX=%f',
!    &         newline=.true., s1=pabbrev, f1=fldmin, f2=sum/gfld%ndpts,
!    &         f3=fldmax)
!       endif

! need to match up soil levels with those requested.
! For the Vtable levels, -88 = all levels, -99 = missing. The units
! vary depending on the level code (e.g. 106 = cm, 103 = m).
! The grib2 standard allows scaling of the units, so make sure the soil level
! units are in cm (as used in the Vtable).
!
!              if ( gfld%ipdtmpl(10) .eq. 106 ) then
              if (( gfld%ipdtmpl(10) .eq. 106 ) .or. 
     &           ( gfld%ipdtmpl(10) .eq. 151 )) then              !! Added by Prof. Mariusz Figurski, 20250110
                if ( ( gfld%ipdtmpl(14) .EQ. -1*(2**07-1) ) .AND.
!    &               ( gfld%ipdtmpl(15) .EQ. -1*(2**31-1) ) ) THEN ! Some compilers cannot
                                                                   ! handle the initial 2**31
                                                                   ! part of the computation, 
                                                                   ! which is an arithmetic 
                                                                   ! overflow on 32 bit signed ints
     &               ( gfld%ipdtmpl(15) .EQ. -2147483647  ) ) THEN
! special UM grib2
                   glevel1 = gfld%ipdtmpl(12)
                   glevel2 = gfld%ipdtmpl(11)
                else
                   glevel1 = 100. * gfld%ipdtmpl(12)*
     &                         (10.**(-1.*gfld%ipdtmpl(11)))
                   glevel2 = 100. * gfld%ipdtmpl(15)*
     &                         (10.**(-1.*gfld%ipdtmpl(14)))
!The following if-block accounts for the 2024 changes to IFS soil fields. They use a local table. 10.01.2024 
!! Added by Prof. Mariusz Figurski, 20250110
                  if (icenter .eq. 98 .and.
     &                gfld%discipline .eq. 2 .and.     !! Added by Prof. Mariusz Figurski, 20250110
     &                (gfld%ipdtmpl(1) .eq. 3 .or.      !! Added by Prof. Mariusz Figurski, 20250110
     &                gfld%ipdtmpl(1) .eq. 0) .and.      !! Added by Prof. Mariusz Figurski, 20250110
     &                gfld%ipdtmpl(10).eq. 151) then    !! Added by Prof. Mariusz Figurski, 20250110
!     &                gfld%discipline .eq. 192 .and.
!     &                gfld%ipdtmpl(1) .eq. 128 .and.
!     &                gfld%ipdtmpl(10).eq. 106) then
                      glevel1 = gfld%ipdtmpl(12)
                      glevel2 = gfld%ipdtmpl(15)                      
!!    Added by Prof. Mariusz Figurski, 20250110  start special for IFS grib2
                      if ( gfld%ipdtmpl(12) .eq. 0) then
                          glevel1 = 0
                      elseif ( gfld%ipdtmpl(12) .eq. 1) then
                          glevel1 = 7
                      elseif ( gfld%ipdtmpl(12) .eq. 2) then
                          glevel1 = 28
                      elseif ( gfld%ipdtmpl(12) .eq. 3) then
                          glevel1 = 100
                      elseif ( gfld%ipdtmpl(12) .eq. 4) then
                          glevel1 = 289
                      endif

                      if ( gfld%ipdtmpl(15) .eq. 0) then
                          glevel2 = 0
                      elseif ( gfld%ipdtmpl(15) .eq. 1) then
                          glevel2 = 7
                      elseif ( gfld%ipdtmpl(15) .eq. 2) then
                          glevel2 = 28
                      elseif ( gfld%ipdtmpl(15) .eq. 3) then
                          glevel2 = 100
                      elseif ( gfld%ipdtmpl(15) .eq. 4) then
                          glevel2 = 289
                      endif
!!   end special block for IFS grib2

                      if ( glevel1 .eq. 100. .and.
     &                     glevel2 .lt. -1.e+9 ) then
                        glevel2 = 255.
                      endif
                  endif
                end if
                TMP8LOOP: do j = 1, maxvar
                  if (((g2code(4,j) .eq. 106) .or. 
     &               (g2code(4,j) .eq. 151)) .and.              !! Added by Prof. Mariusz Figurski, 20250110
     &               (gfld%ipdtmpl(1) .eq. g2code(2,j)) .and.   !! Added by Dr. Haiqing SONG, 20181201
     &               (gfld%ipdtmpl(2) .eq. g2code(3,j)) .and.
     &               (glevel1 .eq. level1(j)) .and.
     &               ((glevel2 .eq. level2(j)) .or.
     &                                   (level2(j) .le. -88))) then
                    my_field = namvar(j)
                    exit TMP8LOOP
                  endif
                enddo TMP8LOOP
                if (j .gt. maxvar ) then
                  write(6,'(a,i6,a)') 'Subsoil level ',
     &               gfld%ipdtmpl(12), 
     &           ' in the GRIB2 file, was not found in the Vtable'
                  cycle MATCH_LOOP
                endif
!MGD         if (debug_level .gt. 50) write(6,*) 'my_field is now ',my_field
              endif

              ! Level (eg. 10000 mb)
              if(gfld%ipdtmpl(10).eq.100) then
                 ! Pressure level (range from 1000mb to 0mb)
                 level=gfld%ipdtmpl(12) *
     &                           (10. ** (-1. * gfld%ipdtmpl(11)))
		 if ( level .lt. pmin ) cycle MATCH_LOOP
              elseif((gfld%ipdtmpl(10).eq.105).or.
     &               (gfld%ipdtmpl(10).eq.118).or.
     &               (gfld%ipdtmpl(10).eq.150))then
                 ! Hybrid level (range from 1 to N)
                 level=gfld%ipdtmpl(12)
              elseif(gfld%ipdtmpl(10).eq.104) then
                 ! Sigma level (range from 10000 to 0)
                 level=gfld%ipdtmpl(12)
              elseif(gfld%ipdtmpl(10).eq.101) then
                 ! MSL
                 level=201300.
              elseif(gfld%ipdtmpl(10).eq.103) then
	         ! Height above ground (m)
		 if (gfld%ipdtmpl(12) .eq. 2. .or. 
     &               gfld%ipdtmpl(12) .eq. 1000. .or.     ! temp fix for hrrr maxref
     &               gfld%ipdtmpl(12) .eq. 10. ) then
                   level=200100.
		 else
                   cycle MATCH_LOOP
		 endif
              elseif((gfld%ipdtmpl(10).ge.206 .and.
     &               gfld%ipdtmpl(10).le.234) .or.
     &              (gfld%ipdtmpl(10).ge.242 .and.     
     &               gfld%ipdtmpl(10).le.254) .or.
     &              (gfld%ipdtmpl(10).eq.200) .or.
     &              (gfld%ipdtmpl(10).eq.10) ) then
                 ! NCEP cloud layers used for plotting
                   level=200100.
              elseif(gfld%ipdtmpl(10).eq.106 .or. 
     &               gfld%ipdtmpl(10).eq.151 .or.       !! Added by Prof. Mariusz Figurski, 20250110
     &               gfld%ipdtmpl(10).eq.1) then
                 ! Misc near ground/surface levels
                 level=200100.
              elseif(gfld%ipdtmpl(10).eq.6) then
                 ! Level of Max wind
                 level=200100.  
              elseif(gfld%ipdtmpl(10).eq.7) then
                 ! Tropopause
                 level=200100.
              else
                 ! If we are here then the Vtable contains a level code
                 ! which we cannot handle. Write an info message and skip it.
                 call mprintf(.true.,INFORM,"Rd_grib2 does not know abou
     &t level code %i (field = %s). Skipping this field. If you want thi
     &s level, rd_grib2.F must be modified", i1 = gfld%ipdtmpl(10),
     & s1 =  my_field )
                 cycle MATCH_LOOP
              endif
              iplvl = int(level)

              ! Store the unpacked 2D slab from the GRIB2 record
              allocate(hold_array(gfld%ngrdpts))
              do j=1,gfld%ngrdpts
                 hold_array(j)=gfld%fld(j)
              enddo

!   Some grids need to be reordered. Until we get an example, this is
!   a placeholder
!             call reorder_it (hold_array, map%nx, map%ny, map%dx, 
!    &                 map%dy, iorder)

              ! When we have reached this point, we have a data array ARRAY 
              ! which has some data we want to save, with field name FIELD 
              ! at pressure level LEVEL (Pa).  Dimensions of this data are 
              ! map%nx and map%ny.  Put this data into storage.

              !print *,'call put_storage',iplvl,my_field,hold_array(55),ith
              !e.g. call put_storage(200100, 'RH', my_field, 1, ith)
!             call mprintf(.true.,DEBUG,"Calling put_storage for
!    &level = %i , field = %s , g2level = %i ", newline=.true.,
!    & i1 = iplvl, s1 = my_field, i2 = gfld%ipdtmpl(12) )

              call put_storage(iplvl,my_field,
     &           reshape(hold_array(1:map%nx*map%ny),
     &           (/map%nx, map%ny/)), map%nx,map%ny)
              deallocate(hold_array)

              ! If Specific Humidity is present on hybrid levels AND 
              ! upper-air RH is missing, see if we can compute RH from 
              ! Specific Humidity.
              if (.not. is_there(iplvl, 'RH') .and.
     &            is_there(iplvl, 'SH') .and.
     &            is_there(iplvl, 'TT') .and.
     &            is_there(iplvl, 'P')) then
                  call g2_compute_rh_spechumd_upa(map%nx,map%ny,iplvl)
                 !call llstor_remove(iplvl, 'SH') !We are done with SH
              endif

              ! If Specific Humidity is present on hybrid levels AND 
              ! upper-air RH is missing, see if we can compute RH from 
              ! Specific Humidity - v2
              if (.not. is_there(iplvl, 'RH') .and.
     &            is_there(iplvl, 'SPECHUMD') .and.
     &            is_there(iplvl, 'THETA') .and.
     &            is_there(iplvl, 'TT')) then
                  call g2_compute_rh_spechumd_upa2(map%nx,map%ny,iplvl)
              endif

              ! If Temperature and Theta are present on hybrid levels AND
              ! upper-air PRESSURE is missing, see if we can compute PRESSURE from
              ! Temperature and Theta
              if (.not. is_there(iplvl, 'PRESSURE') .and.
     &            is_there(iplvl, 'THETA') .and.
     &            is_there(iplvl, 'TT')) then
                  call g2_compute_pressure_tth_upa(map%nx,map%ny,iplvl)
              endif

              ith=ith+1
              exit MATCH_LOOP

           endif ! Selected param.


         enddo MATCH_LOOP

         ! Deallocate arrays decoding GRIB2 record.
         call gf_free(gfld)

         enddo NUM_FIELDS


      enddo VERSION ! skgb


       if ( debug_level .gt. 100 ) then
         call mprintf (.true.,DEBUG,
     &   "G2 total number of fields found = %i ",newline=.true.,i1=itot)
       end if

       CALL BACLOSE(junit,IOS)

       nullify(gfld%local)            ! must be nullified before opening next file
       ireaderr=1
      else 
        call mprintf (.true.,DEBUG,"open status failed because %i ",
     &                newline=.true., i1=ios)
        hdate = '9999-99-99_99:99:99'
        ireaderr=2
      endif ! ireaderr check 

      END subroutine rd_grib2

!*****************************************************************************!
! Subroutine edition_num                                                      !
!                                                                             !
! Purpose:                                                                    !
!    Read one record from the input GRIB2 file.  Based on the information in  !
!    the GRIB2 header and the user-defined Vtable, decide whether the field in!
!    the GRIB2 record is one to process or to skip.  If the field is one we   !
!    want to keep, extract the data from the GRIB2 record, and pass the data  !
!    back to the calling routine.                                             !
!                                                                             !
! Argument list:                                                              !
!    Input:                                                                   !
!       JUNIT   : "Unit Number" to open and read from.  Not really a Fortran  !
!                 unit number, since we do not do Fortran I/O for the GRIB2   !
!                 files.  Nor is it a UNIX File Descriptor returned from a C  !
!                 OPEN statement.  It is really just an array index to the    !
!                 array (IUARR) where the UNIX File Descriptor values are     !
!                 stored.                                                     !
!       GRIB2FILE: File name to open, if it is not already open.              !
!                                                                             !
!    Output:                                                                  !
!       GRIB_EDITION: Set to 1 for GRIB and set to 2 for GRIB2                ! 
!       IERR     : Error flag: 0 - no error on read from GRIB2 file.          !
!                              1 - Hit the end of the GRIB2 file.             !
!                              2 - The file GRIBFLNM we tried to open does    !
!                                  not exist.                                 !
! Author: Paula McCaslin                                                      !
! NOAA/FSL                                                                    !
! Sept 2004                                                                   !
!*****************************************************************************!
      
      SUBROUTINE edition_num(junit, gribflnm, 
     &  grib_edition, ireaderr)

      use grib_mod
      use params
      use module_debug

      parameter(msk1=32000,msk2=4000)
      character(len=1),allocatable,dimension(:) :: cgrib
      integer :: listsec0(3)
      integer :: listsec1(13)
      character(len=*)  :: gribflnm
      integer :: lskip, lgrib
      integer :: junit
      integer :: grib_edition
      integer :: i, j, ireaderr
      integer :: currlen

      character(len=4) :: ctemp
      character(len=4),parameter :: grib='GRIB',c7777='7777'

C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  SET ARGUMENTS

      itot=0
      icount=0
      iseek=0
      lskip=0
      lgrib=0
      currlen=0

!/* IOS Return Codes from BACIO:  */
!/*  0    All was well                                   */
!/* -1    Tried to open read only _and_ write only       */
!/* -2    Tried to read and write in the same call       */
!/* -3    Internal failure in name processing            */
!/* -4    Failure in opening file                        */
!/* -5    Tried to read on a write-only file             */
!/* -6    Failed in read to find the 'start' location    */
!/* -7    Tried to write to a read only file             */
!/* -8    Failed in write to find the 'start' location   */
!/* -9    Error in close                                 */
!/* -10   Read or wrote fewer data than requested        */

!if ireaderr =1 we have hit the end of a file. 
!if ireaderr =2 we have hit the end of all the files. 
!if ireaderr =3 beginning characters 'GRIB' not found

      ! Open a byte-addressable file.
      CALL BAOPENR(junit,gribflnm,IOS)
      if (ios.eq.0) then 

         ! Search opend file for the next GRIB2 messege (record).
         call skgb(junit,iseek,msk1,lskip,lgrib)

         ! Check for EOF, or problem
         call mprintf((lgrib.eq.0),ERROR,
     &     "Grib2 file or date problem, stopping in edition_num.",
     &     newline=.true.)
 
         ! Check size, if needed allocate more memory.
         if (lgrib.gt.currlen) then
            if (allocated(cgrib)) deallocate(cgrib)
            allocate(cgrib(lgrib),stat=is)
            currlen=lgrib
         endif

         ! Read a given number of bytes from unblocked file.
         call baread(junit,lskip,lgrib,lengrib,cgrib)

         ! Check for beginning of GRIB message in the first 100 bytes
         istart=0
         do j=1,100
            ctemp=cgrib(j)//cgrib(j+1)//cgrib(j+2)//cgrib(j+3)
            if (ctemp.eq.grib ) then
              istart=j
              exit
            endif
         enddo
         if (istart.eq.0) then
            ireaderr=3
            print*, "The beginning 4 characters >GRIB< were not found."
         endif
   
         ! Unpack Section 0 - Indicator Section to extract GRIB edition field
         iofst=8*(istart+5)
         call gbyte(cgrib,discipline,iofst,8)     ! Discipline
         iofst=iofst+8
         call gbyte(cgrib,grib_edition,iofst,8)   ! GRIB edition number

!        print *, 'ungrib - grib edition num',  grib_edition
         CALL BACLOSE(junit,IOS)
         ireaderr=1
      else if (ios .eq. -4) then
        call mprintf(.true.,ERROR, 
     &    "edition_num: unable to open %s",newline=.true.,s1=gribflnm)
      else 
         print *,'edition_num: open status failed because',ios,gribflnm
         ireaderr=2
      endif ! ireaderr check 

      END subroutine edition_num

!*****************************************************************************!

      SUBROUTINE g2_compute_rh_spechumd_upa(ix, jx, iiplvl)
      ! Compute relative humidity from specific humidity in the upper air.
      use storage_module
      implicit none
      integer :: ix, jx
      integer :: iiplvl
      real :: lat1, lon1, dx, dy
      real, dimension(ix,jx) :: T, P, RH, Q
    
      real, parameter :: svp1=611.2
      real, parameter :: svp2=17.67
      real, parameter :: svp3=29.65
      real, parameter :: svpt0=273.15
      real, parameter :: eps = 0.622
    
      real startlat, startlon, deltalat, deltalon

      call get_storage(iiplvl, 'P', P, ix, jx)
      call get_storage(iiplvl, 'TT', T, ix, jx)
      call get_storage(iiplvl, 'SH', Q, ix, jx)
    
      rh=1.E2*(p*q/(q*(1.-eps)+eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
     
      call put_storage(iiplvl, 'RH', rh, ix, jx)
    
      end subroutine g2_compute_rh_spechumd_upa

!*****************************************************************************!

      SUBROUTINE g2_compute_rh_spechumd_upa2(ix, jx, iiplvl)
      ! Compute relative humidity from specific humidity in the upper air.
      use storage_module
      implicit none
      integer :: ix, jx
      integer :: iiplvl
      real :: lat1, lon1, dx, dy
      real, dimension(ix,jx) :: T, TH, RH, Q, P
    
      real, parameter :: svp1=611.2
      real, parameter :: svp2=17.67
      real, parameter :: svp3=29.65
      real, parameter :: svpt0=273.15
      real, parameter :: eps = 0.622
    
      real startlat, startlon, deltalat, deltalon

      call get_storage(iiplvl, 'THETA', TH, ix, jx)
      call get_storage(iiplvl, 'TT', T, ix, jx)
      call get_storage(iiplvl, 'SPECHUMD', Q, ix, jx)
    
      p=1.e5*(t/th)**(1005/287.05)
     
      rh=1.E2*(p*q/(q*(1.-eps)+eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
     
      call put_storage(iiplvl, 'RH', rh, ix, jx)
    
      end subroutine g2_compute_rh_spechumd_upa2

!*****************************************************************************!

      SUBROUTINE g2_compute_pressure_tth_upa(ix, jx, iiplvl)
      ! Compute relative humidity from specific humidity in the upper air.
      use storage_module
      implicit none
      integer :: ix, jx
      integer :: iiplvl
      real :: lat1, lon1, dx, dy
      real, dimension(ix,jx) :: T, TH, P
    
      real, parameter :: svp1=611.2
      real, parameter :: svp2=17.67
      real, parameter :: svp3=29.65
      real, parameter :: svpt0=273.15
      real, parameter :: eps = 0.622
    
      real startlat, startlon, deltalat, deltalon

      call get_storage(iiplvl, 'THETA', TH, ix, jx)
      call get_storage(iiplvl, 'TT', T, ix, jx)
    
      p=1.e5*(t/th)**(1005/287.05)
     
      call put_storage(iiplvl, 'PRESSURE', p, ix, jx)
    
      end subroutine g2_compute_pressure_tth_upa

!*****************************************************************************!

      subroutine ncep_grid_num (pnum)
!
!  Find the grib number for descriptive  labelling.
!  Grib2 doesn't have a grid-number entry, so we have to figure it out
!  from the parameters. 
!
      use gridinfo       ! Included to define map%
      integer :: pnum
      real, parameter :: eps = .01
      character (len=8) :: tmp8

!     write(6,*) 'begin ncep_grid_num'
!     write(6,*) 'dx = ',map%dx,' pnum = ',pnum,' nx = ',map%nx
      tmp8 = '        '
      if (pnum .eq. 30) then            ! lambert conformal
        if ( abs(map%dx - 12.19058) .lt. eps .and. map%nx .eq. 614) then
          write(tmp8,'("GRID 218")') 
        else if (abs(map%dx - 40.63525) .lt. eps 
     &     .and. map%nx .eq. 185) then
          write(tmp8,'("GRID 212")') 
        else if (abs(map%dx - 40.63525) .lt. eps 
     &     .and. map%nx .eq. 151) then
          write(tmp8,'("GRID 236")') 
        else if (abs(map%dx - 81.2705) .lt. eps 
     &     .and. map%nx .eq. 93) then
          write(tmp8,'("GRID 211")') 
        else if (abs (map%dx - 32.46341) .lt. eps 
     &     .and. map%nx .eq. 349) then
          write(tmp8,'("GRID 221")') 
        else if (abs(map%dx - 20.317625) .lt. eps 
     &     .and. map%nx .eq. 301) then
          write(tmp8,'("GRID 252")') 
        else if (abs(map%dx - 13.545087) .lt. eps 
     &     .and. map%nx .eq. 451) then
          write(tmp8,'("GRID 130")') 
        endif
      else if (pnum .eq. 20) then     ! polar stereographic
        if (abs(map%dx - 15.0) .lt. eps) then
          write(tmp8,'("GRID  88")') 
        endif
      else if (pnum .eq. 0) then      ! lat/lon
        if (abs(map%dx - 1.) .lt. eps .and. map%nx .eq. 360) then
          write(tmp8,'("GRID   3")') 
        else if (abs(map%dx - 0.5) .lt. eps .and. map%nx .eq. 720) then
          write(tmp8,'("GRID   4")') 
        endif
      endif
      map%source(25:32) = tmp8
!     write(6,*) 'map%source = ',map%source
      end subroutine ncep_grid_num
!*****************************************************************************!

      function earth_radius (icode, iscale, irad_m)
! Grib2 Code Table 3.2. Returns the spherical earth's radius in km.
      use module_debug
      real :: earth_radius
      integer :: icode
      integer :: iscale, irad_m
      if ( icode .eq. 0 ) then
        earth_radius = 6367470. * .001
      else if ( icode .eq. 1) then
        earth_radius = 0.001 * float(irad_m) / 10**iscale
      else if ( icode .eq. 6 ) then
        earth_radius = 6371229. * .001
      else if ( icode .eq. 8 ) then
        earth_radius = 6371200. * .001
      else
        call mprintf(.true.,ERROR,
     &    "unknown earth radius for code %i",newline=.true.,i1=icode)
      endif
      end function earth_radius
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
GRIB1| Level| From |  To  | metgrid  | metgrid  | metgrid                                  |GRIB2|GRIB2|GRIB2|GRIB2|
Param| Type |Level1|Level2| Name     | Units    | Description                              |Discp|Catgy|Param|Level|
-----+------+------+------+----------+----------+------------------------------------------+-----------------------+
 129 | 100  |   *  |      | HGT      | m        | Height                                   |  0  |  3  |  5  | 100 |
 130 | 100  |   *  |      | TT       | K        | Temperature                              |  0  |  0  |  0  | 100 |
 131 | 100  |   *  |      | UU       | m s-1    | U                                        |  0  |  2  |  2  | 100 |
 132 | 100  |   *  |      | VV       | m s-1    | V                                        |  0  |  2  |  3  | 100 |
 157 | 100  |   *  |      | RH       | %        | Relative Humidity                        |  0  |  1  |  1  | 100 |
 165 |  1   |   0  |      | UU       | m s-1    | U                 at 10 m                |  0  |  2  |  2  | 103 |
 166 |  1   |   0  |      | VV       | m s-1    | V                 at 10 m                |  0  |  2  |  3  | 103 |
 167 |  1   |   0  |      | TT       | K        | Temperature       at  2 m                |  0  |  0  |  0  | 103 |
 168 |  1   |   0  |      | DEWPT    | K        |                                          |  0  |  0  |  6  | 103 |
     |  1   |   0  |      | RH       | %        | Relative Humidity at 2 m                 |  0  |  1  |  1  | 103 |
 172 |  1   |   0  |      | LANDSEA  | 0/1 Flag | Land/Sea flag (1=land, 0 or 2=sea)       |  2  |  0  |  0  |   1 |
 134 |  1   |   0  |      | PSFC     | Pa       | Surface Pressure                         |  0  |  3  |  0  |   1 |
 151 |  1   |   0  |      | PMSL     | Pa       | Sea-level Pressure                       |  0  |  3  |  0  | 101 |
 235 |  1   |   0  |      | SKINTEMP | K        | Skin temperature                         |  0  |  0  | 17  |   1 |
  18 | 112  |   0  |   7  | ST000007 | K        | T of 0-7 cm  ground layer                |  2  |  3  |  18 | 151 | 
  18 | 112  |   7  |  28  | ST007028 | K        | T of 7-28 cm ground layer                |  2  |  3  |  18 | 151 | 
  18 | 112  |  28  | 100  | ST028100 | K        | T of 28-100 cm ground layer              |  2  |  3  |  18 | 151 | 
  18 | 112  | 100  | 289  | ST100289 | K        | T of 100-289 cm ground layer             |  2  |  3  |  18 | 151 | 
  25 | 112  |   0  |   7  | SM000007 | m3 m-3   | Soil moisture of 0-7 cm ground layer     |  2  |  0  |  25 | 151 | 
  25 | 112  |   7  |  28  | SM007028 | m3 m-3   | Soil moisture of 7-28 cm ground layer    |  2  |  0  |  25 | 151 | 
  25 | 112  |  28  | 100  | SM028100 | m3 m-3   | Soil moisture of 28-100 cm ground layer  |  2  |  0  |  25 | 151 | 
  25 | 112  | 100  | 289  | SM100289 | m3 m-3   | Soil moisture of 100-289 cm ground layer |  2  |  0  |  25 | 151 | 
-----+------+------+------+----------+---------+-----------------------------------------+-------------------------+
#

Vtable.ECMWF_GRIB2/Vtable.ECMWF

1
2
-rw-rw-r-- 1 wpsze wpsze 321M Apr  1 16:27 ERA5:2025-03-26_00-ECMWF_GRIB2
-rw-rw-r-- 1 wpsze wpsze 286M Apr  1 14:26 ERA5:2025-03-26_00-Vtable.ECMWF
  • The intermediate file (ERA5.xxx) by Vtable.ECMWF_GRIB2 is larger file size than that of Vtable.ECMWF
  • ERA5.xxx by Vtable.ECMWF fails in init_atmosphere_model step.

init_atmosphere_model

Error

  • MPASv7.3 init_atmosphere_model
1
2
3
4
ERROR: extrap_type == 2 not implemented for target_z >= zf(1,nz)
ERROR: *****************************************************************
ERROR: Error in interpolation of t(k,iCell) for k=  46, iCell=         1
CRITICAL ERROR: ********************************************************
  • set 'linear',
1
2
3
&interpolation_control
    config_extrap_airtemp = 'linear'
/
1
CRITICAL ERROR: Error in interpolation of st_fg to MPAS grid: num_st = 27244
  • Error creating initial conditions with NCEP FNL ds083.2 | May 8, 2024
    • I'm trying to create the initial conditions file for MPAS (v8.1) using data obtained from the NCEP FNL ds083.2
    • Just had to change config_nfgsoillevels to 2 and config_extrap_airtemp to linear in namelist.init_atmosphere. Once those two modifications were made init_atmosphere_model and atmosphere_model ran without issues.
    • Note that FNL data only contains 2-level of soil information, which is why we need to set config_nfgsoillevels=2. The top of the FNL data is 10 hpa, and therefore extrapolation is needed if MPAS top level is above 10 hpa.
    • For the option of config_extrap_airtemp, would you please let me know what is the model top of your MPAS run? Thanks.
  • Initial condition CRITICAL ERROR: Error in interpolation of st_fg to MPAS | Feb 16, 2024
  • 'Error in interpolation of st_fg to MPAS grid' with variable-resolution mesh simulation | Jul 23, 2020
    • The error message generally indicates that there are grid cells that have a soil temperature value less than or equal to zero.
    • I solved the problem by downloading the "Vtable from the RDR page/WRF" and using WPS version 4.0 or higher.
  • CRITICAL ERROR: Error in interpolation of st_fg to MPAS grid | Start dateMay 3, 2019
    • This error message is an indication that there were one or more points in the MPAS mesh that received a zero or negative value for soil temperature.
    • Can you also verify that there are soil moisture and soil temperature fields in your WPS intermediate file, and that the number of soil levels in your intermediate files matches the value of the config_nfgsoillevels in your namelist.init_atmosphere file?
  • CRITICAL ERROR: Error in interpolation of st_fg to MPAS grid: num_st = 18642 | Apr 12, 2023
    • MPAS can only process soil data at specific levels with specific field names. If your input soil data are on levels or with names unrecognized by MPAS, then it won't work.
    • Please take a look at the code "mpas_init_atm_cases.F", and the lines 3908 - 4118 process soil temperature data. This might give you some idea how MPAS works.
  • Trying to Use RAP as IC/BC | Start dateJun 26, 2019
    • The message "Error in interpolation of st_fg to MPAS grid: num_st = XXXX" is generally indicative of an issue with the soil temperature. Would it be possible to use some of the utilities that come with the WRF Pre-processing System (WPS) to check that there are valid soil temperature data in your intermediate file? Specifically, you could use the "rd_intermediate" utility to check that there are fields named, e.g., ST000010, and you could use the "int2nc" utility to convert the intermediate file to netCDF format, where it would be easier to view the soil fields with, e.g., ncview.
  • Use WPSv4.6 with above renewed Vtable and rd_grib2.F

One more issue, 資料集可用性 2024-11-12T12:00:00Z–2025-04-16T12:00:00Z

  • 資料集可用性 2024-11-12T12:00:00Z–2025-04-16T12:00:00Z
  • The grib files from ecmwf.ini and from aws are different. It will make different/debug on the ungrib step and init_atmosphere_model step.
    • Not all different.

It is because,

parameter "config_nfglevels"

  • configuration problem for parameter "config_nfglevels" | Sep 2, 2021
    • If you have 37 isobaric levels of atmospheric data, plus surface fields (u10, v10, t2, etc.), that would give a total of 38 levels from the perspective of the init_atmosphere_model program. As the error message states, config_nfglevels must therefore be set to at least 38 -- so value of 38, 39, 40, 41, etc. should all work.
    • Internally, the 'config_nfglevels' namelist option is used by the init_atmosphere_model program to allocate 3-d arrays to hold the horizontally interpolated first-guess ("fg" in the config_nfglevels name) fields. Allocating too few vertical levels in these 3-d arrays is problematic, but allocating extra levels -- which will simply go unused -- isn't a problem, though it does waste memory.

Done

1
MPAS-120km_2022063000.init.nc

and MPASv7.3 runs (./atmosphere_model) with IFS successfully as well (both with Vtable.ECMWF_GRIB2 above).

cronjob

For 00z,

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#!/bin/sh
#################################################################
# https://data.ecmwf.int/ecpds/home/opendata/20250402/00z/ifs/0p25/oper/20250402000000-0h-oper-fc.grib2
#################################################################
source /home/wpsze/micromamba/etc/profile.d/micromamba.sh
micromamba activate venv

BASE_DIR=$(pwd)
echo "Base Dir = ${BASE_DIR}"

for ifsday in $(seq 0 1 3); do # day-shift = 0, 1, 2, 3

    YYYYMMDDHH=$(date -d "today - ${ifsday} day" +%Y%m%d00)

    yyyy=${YYYYMMDDHH:0:4}
    mm=${YYYYMMDDHH:4:2}
    dd=${YYYYMMDDHH:6:2}
    hh=${YYYYMMDDHH:8:2}

    yyyymm=${yyyy}${mm}
    yyyymmdd=${yyyy}${mm}${dd}

    tmp_dir="./${yyyy}${mm}/${yyyy}${mm}${dd}"
    mkdir -p ${tmp_dir}
    cd ${tmp_dir}
    echo "Current Dir: $(pwd)"

    limit_cpu_core=15
    count_cpu_core=1

    #============== 00z ===================================
    tmp_ifs_fname="${YYYYMMDDHH}0000-0h-oper-fc.grib2"

    echo "Now, ${tmp_ifs_fname} ... "

	    echo "Downloading ... ${tmp_ifs_fname}"
	    #wget -c --no-check-certificate https://data.ecmwf.int/forecasts/${yyyy}${mm}${dd}/00z/ifs/0p25/oper/${tmp_ifs_fname}
	    wget -c --no-check-certificate https://data.ecmwf.int/ecpds/home/opendata/${yyyy}${mm}${dd}/00z/ifs/0p25/oper/${tmp_ifs_fname}


    res=$(( ${count_cpu_core} % ${limit_cpu_core} ))
    if [[ ${res} == 0 ]]; then
	    echo "please wait ${tmp_ifs_fname}, ${count_cpu_core}, ${res}."
	    wait
    fi
    
    cd ${BASE_DIR}
done

Dissemination schedule

Time available (UTC)
Time available (UTC)
  • HKT = UTC + 8
  • How to set the crontab: $ crontab -e
1
2
$ crontab -l
30 17 * * * cd /home/wpsze/mpas/IFS; sh wpsze_download_IFS-cronjob-00z.sh

References

  1. ECMWF releases a much larger open dataset | 29 February 2024
  2. Open data | ECMWF real-time forecast data
  3. ECMWF IFS 免费开放预报升级到 0.25° 分辨率 | 2024年02月26日
  4. 重大升级！欧洲中期天气预报中心AI加持下免费开放了更高分辨率、更多参数的气象数据！| 2024年03月01日
NWP
#WRF #WPS #MPAS #IC/BC #ECMWF #IFS
NWP | ECMWF | IFS | Open data
https://waipangsze.github.io/2024/11/12/ECMWF-IFS/
Author
wpsze
Posted on
November 12, 2024
Updated on
October 27, 2025
Licensed under