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

webisu commented 7 months ago

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 commented 7 months ago

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

webisu commented 7 months ago

Alex,

I believe that I am on the latest version of develop according to the following:W

@.***:~/ip/NCEPLIBS-ip$ git status On branch develop Your branch is up to date with 'origin/develop'.

nothing to commit, working tree clean @.:~/ip/NCEPLIBS-ip$ git pull Already up to date. @.:~/ip/NCEPLIBS-ip$ cat VERSION 5.0.0 @.***:~/ip/NCEPLIBS-ip$

On Wed, Apr 17, 2024 at 4:59 PM Alex Richert @.***> wrote:

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

— Reply to this email directly, view it on GitHub https://github.com/NOAA-EMC/NCEPLIBS-ip/issues/238#issuecomment-2062319372, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEIB7ZXBY4ZEGVIDYNRXBNTY53PC5AVCNFSM6AAAAABGLKEL36VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDANRSGMYTSMZXGI . You are receiving this because you authored the thread.Message ID: @.***>

AlexanderRichert-NOAA commented 7 months ago

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

webisu commented 7 months ago

Alex, I just pushed a version of wgrib2 "ip_wne_1" which calls nceplib-ip instead of the ip2lib_d.

My test was

$ wgrib2 nam.t00z.bgrd3d00.tm00.grib2.c0 -new_grid ncep grid 2 junk

and then visualize junk using GrADS. The resulting plot looks like a crashing ocean wave rather than the NAM grid. You can get the nam file from wcoss or use the copy that I have on

/cpc/home/wd51we/grib2_examples

CPC workstations are not mounted on EMC linux machines, so you have to log on to

cpc-lw-webisuzaki.ncep.noaa.gov

using your EMC linux password.

BTW this is the grid of the input file

@.***:~/grib2_examples$ wgrib2 nam.t00z.bgrd3d00.tm00.grib2.c0 -d 1 -grid 1:0:grid_template=32769:winds(grid): I am not an Arakawa E-grid. I am rotated but have no rotation angle. I am staggered. What am I? (954 x 835) units 1e-06 input WE:SN output WE:SN res 56 lat0 -7.491000 lat-center 54.000000 dlat 0.108000 lon0 215.866000 lon-center 254.000000 dlon 0.126000 #points=796590

So far, I have tested to/from lat-lon grids, to lambert-conformal grids, from WMO rotated lat-lon grds and they all look good.

Wesley

On Thu, Apr 18, 2024 at 11:11 AM Alex Richert @.***> wrote:

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

— Reply to this email directly, view it on GitHub https://github.com/NOAA-EMC/NCEPLIBS-ip/issues/238#issuecomment-2064170491, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEIB7ZSD2G3YRKYMK3NIWO3Y57PDJAVCNFSM6AAAAABGLKEL36VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDANRUGE3TANBZGE . You are receiving this because you authored the thread.Message ID: @.***>

AlexanderRichert-NOAA commented 6 months ago

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

webisu commented 6 months ago

Alex,

Wgrib2 uses mk_gdt(..)

this routine takes the sec[] and makes a g2lib igdtmpl() array for ipolates2 library
igdt is used by NCEP grib and ipolates2 libraries

// 3.32769: Rot Lat/Lon Non-E Staggered grid (Arakawa) // {32769, 21, 0, {1,1,4,1,4,1,4,4,4,4,4,-4,4,1,-4,4,4,4,1,4,4} }, {32769, "\x01\x01\x04\x01\x04\x01\x04\x04\x04\x04\x04\x0c\x04\x01\x0c\x04\x04\x04\x01\x04\x04" }, };

The New_grid.c code calls mk_gdt,

    if (mk_gdt(new_sec, &(save->gdtnum_out), &(save->gdt_out[0]),

&(save->gdt_out_size) )) fatal_error("new_grid: encoding output gdt failed","");

The gdt* variables have the input grid information., and used in the call to ipolates_grib2_single_field(..) or ipolatev_grib2_single_field(..)

Wesley

On Mon, Apr 29, 2024 at 4:49 PM Alex Richert @.***> wrote:

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

— Reply to this email directly, view it on GitHub https://github.com/NOAA-EMC/NCEPLIBS-ip/issues/238#issuecomment-2083636420, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEIB7ZUKLPZ4BEP3XFFWT6LY72W3HAVCNFSM6AAAAABGLKEL36VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDAOBTGYZTMNBSGA . You are receiving this because you authored the thread.Message ID: @.***>

AlexanderRichert-NOAA commented 6 months ago

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 commented 6 months ago

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 commented 6 months ago

Alex, The code for the NCEP rotated lat-lon grid was part of the grib1 ip library.

In the mk_kgds.c kgds[8]= floor(dlon*1000.0+0.5); The kgds array is an integer array and stores angles to the nearest milli-degree.

The grib2 grid definition table for ncep rot lat-lon grid https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp3-32769.shtml

Stores dlon as an integer in microdegrees. This is why the values are off by 1000 between gds and gdt.

Wesley

On Wed, May 8, 2024 at 2:03 PM Alex Richert @.***> wrote:

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 https://github.com/webisu is this something that can be fixed through wgrib2, or does this still indicate a problem with the ip library?

— Reply to this email directly, view it on GitHub https://github.com/NOAA-EMC/NCEPLIBS-ip/issues/238#issuecomment-2101121162, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEIB7ZUE7ZWAIB6IXHXZV23ZBJSGZAVCNFSM6AAAAABGLKEL36VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCMBRGEZDCMJWGI . You are receiving this because you were mentioned.Message ID: @.***>

AlexanderRichert-NOAA commented 6 months ago

@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 commented 6 months ago

@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 commented 6 months ago

@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 commented 6 months ago

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 commented 6 months ago

Alex,

When wgrib2 calls iplib.v3.0.0, it translates the grib2 sec[][] to the grib1 kgds by the routine mk_kgds(..). When wgrib2 call ip2lib_d, it translates the grib2 sec[][] to the GDT by the routine mk_gdt(..)

mk_gdt translates from sec[][] to gdt using the following table

// 3.32768: Rot Lat/Lon E-grid (Arakawa) // {32768, 19, 0, {1,1,4,1,4,1,4,4,4,4,4,-4,4,1,-4,4,4,4,1} }, {32768, "\x01\x01\x04\x01\x04\x01\x04\x04\x04\x04\x04\x0c\x04\x01\x0c\x04\x04\x04\x01" },

// 3.32769: Rot Lat/Lon Non-E Staggered grid (Arakawa) // {32769, 21, 0, {1,1,4,1,4,1,4,4,4,4,4,-4,4,1,-4,4,4,4,1,4,4} }, {32769, "\x01\x01\x04\x01\x04\x01\x04\x04\x04\x04\x04\x0c\x04\x01\x0c\x04\x04\x04\x01\x04\x04" },

I took the commented lines from some unknown code in g2c. Perhaps the lines have changed since 2016.

Wesley

On Mon, May 13, 2024 at 12:52 PM Alex Richert @.***> wrote:

Thanks @GeorgeGayno-NOAA https://github.com/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.

— Reply to this email directly, view it on GitHub https://github.com/NOAA-EMC/NCEPLIBS-ip/issues/238#issuecomment-2108210446, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEIB7ZTUVXPAWF3CQYFD6ATZCDVWDAVCNFSM6AAAAABGLKEL36VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCMBYGIYTANBUGY . You are receiving this because you were mentioned.Message ID: @.***>

NOAA-EMC / NCEPLIBS-ip

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