Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix support for rotated lat/lon in ip2 code to support wgrib2 #238

Closed
webisu opened this issue Apr 17, 2024 · 14 comments · Fixed by #241
Closed

fix support for rotated lat/lon in ip2 code to support wgrib2 #238

webisu opened this issue Apr 17, 2024 · 14 comments · Fixed by #241
Assignees
Labels
bug Something isn't working

Comments

@webisu
Copy link

webisu commented Apr 17, 2024

Using NCEPLIBS-ip 5.0.0
Interpolating a from COSMO (grid template 1) rot lat-lon grid to NCEP grid 2 looks good.
Interpolating from a NAM (grid template 32769) rot lat-lon grid to NCEP grid 2 looks bad.
The grid should be the NAM grid, and it is completely different.

These results are from testing with wgrib2 with preliminary NCEPLIBS-ip support.

@AlexanderRichert-NOAA
Copy link
Contributor

Can you try with the latest copy of develop, or in any case something at least as recent as commit 30beaa0? The changes we made for Arakawa grids are not yet in a release (though they can be if needed).

@webisu
Copy link
Author

webisu commented Apr 18, 2024 via email

@AlexanderRichert-NOAA
Copy link
Contributor

Would you be able to provide a test case (input file+call to whatever ip subroutine)?

@webisu
Copy link
Author

webisu commented Apr 19, 2024 via email

@AlexanderRichert-NOAA
Copy link
Contributor

How does wgrib2 determine the input grid type? Is it reading 32769 from the file and passing that to ip?

@webisu
Copy link
Author

webisu commented Apr 30, 2024 via email

@AlexanderRichert-NOAA
Copy link
Contributor

AlexanderRichert-NOAA commented Apr 30, 2024

Updates:

  • I can reproduce the problem using wgrib2 3.1.1. Using ipolates=1, I can process the NAM file using wgrib2 nam.t00z.bgrd3d00.tm00.grib2 -new_grid ncep grid 2 regridded_ncep2, and the plots in grads looks reasonable, i.e., the contours are in the right part of the map (I'm plotting soil temp since it should match up cleanly with continent boundaries). With wgrib2 'ip_wne_1' branch (commit 4024a36), it's obvious that the coordinates are totally wrong. I see what you mean now about crashing ocean waves...
  • When I run wgrib2 with ip@5, it appears from the gprof output that the relevant ip routines are not being called (or any ip routines as far as I can tell), whereas in the gprof output for my correctly working copy of wgrib2, I can see, for instance, gdswzdcd_mod:gdswzdcd_vect_rot getting called many times over. I can see all the ip routines in the nm output of my wgrib2+ip@5 executable, so hopefully it's just a matter of a missing function call, grid number mismatch, whatever, as opposed to something "deeper" with the arithmetic in ip. not sure about this-- it appears to be doing something with ip...

I'll update when I know more.

@AlexanderRichert-NOAA
Copy link
Contributor

I can find at least one source of trouble here. When I run the correctly working version of wgrib2+internal ip, it uses KGDS(8) to derive rlon0. The ip library, in src/ip_rot_equid_cylind_grid_mod.F90, believes than rlon0 is to be derived from igdtmpl(21), however the equivalent value of kgds(8) (equivalent meaning scaled by a factor of 1000) is being passed in under igdtmpl(16) by wgrib2. @webisu is this something that can be fixed through wgrib2, or does this still indicate a problem with the ip library?

@webisu
Copy link
Author

webisu commented May 8, 2024 via email

@AlexanderRichert-NOAA
Copy link
Contributor

@webisu please take a look at https://github.com/NOAA-EMC/NCEPLIBS-ip/tree/arakawa_grib2_fix and see if that works for you. What I'm now working to understand is why the logic for the grib1 and grib2 grid setups are so different.

6f9dcfb

@AlexanderRichert-NOAA
Copy link
Contributor

@GeorgeGayno-NOAA I know the code is from a while ago, but I was wondering if you might have any insight into why ip's logic for the non-E Arakawa grid setup is so different between init_grib1 and init_grib2? For instance, in init_grib2, it's adding 90 degrees to RLAT0, and the [NESW]BD variables are derived very differently. It's only by copying and pasting the logic from init_grib1 into init_grib2 that I can make the above test file look right (the file is a NAM output file, nam.t00z.bgrd3d00.tm00.grib2).

@GeorgeGayno-NOAA
Copy link
Contributor

@GeorgeGayno-NOAA I know the code is from a while ago, but I was wondering if you might have any insight into why ip's logic for the non-E Arakawa grid setup is so different between init_grib1 and init_grib2? For instance, in init_grib2, it's adding 90 degrees to RLAT0, and the [NESW]BD variables are derived very differently. It's only by copying and pasting the logic from init_grib1 into init_grib2 that I can make the above test file look right (the file is a NAM output file, nam.t00z.bgrd3d00.tm00.grib2).

Why are init_grib1 and init_grib2 different? You will need to ask the WMO. They set the standard.

In GRIB2, RLAT0 is specified as the south pole of the projection. In GRIB1, it is the center latitude of the grid. [NESW]BD are computed from the corner points. In GRIB1, the corners are in rotated space. They need to be unrotated. In GRIB2, they are already in unrotated space.

https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp3-1.shtml
https://www.nco.ncep.noaa.gov/pmb/docs/on388/tabled.html

@AlexanderRichert-NOAA
Copy link
Contributor

Thanks @GeorgeGayno-NOAA for the info.

@webisu-- Here is the current state of my understanding of this issue.

  • When working correctly, wgrib2 takes a grib2 message with GDTN=32769 and passes it to the built-in iplib.v3.0.0, wherein it translates the GDT parameters into the corresponding grib1 GDT parameters but otherwise applies the same formulas in deriving the grid mapping. This is consistent with how the NAM generates its grib2 outputs via ncep_post (grib2_module.f). In other words, iplib.v3.0.0 and ncep_post agree on the meaning of each grib2 GDT element and therefore operate self consistently.
  • On the other hand, the current ip code (post-ip2 merge, I suspect) and the WMO have different ideas from that. Specifically, the definitions of several of the GDT parameters are different, and the math that the currect ip library uses to derive the grid is very different (see George's comment above re: parameter definitions/usage).

Since of course I want to avoid breaking existing functionality, I'm hesitant to remove/replace the existing grid derivation for GDTN=32769 in the ip library unless there's a specific and solid reason to think that it's wrong. That said, one possible approach would be to add an optional parameter to the non-E version of init_grib2() and its various calling subroutines that would use a grid definition consistent with what's in ncep_post. Does that make sense? In what I'm envisioning, you would add something like ncep_post=true to the subroutine call to access that code path in ip.

@webisu
Copy link
Author

webisu commented May 20, 2024 via email

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment