Combine ip and ip2 libraries

[kgerheiser]

I've been working on this for the past two weeks and I finally have a combined ip library. All the grib1/grib2 tests pass in this one repo. This is a Draft PR while I add some final touches. The history is based of off ip2, but I plan to make the actual PR into NCEPLIBS-ip.

Instead of operating on the raw grib1/grib2 descriptors the code has been refactored to make use of some new objects.

To start with there's the abstract type ip_grid_descriptor in ip_grid_descriptor_mod.f90, and there's two subclasses grib1_descriptor and grib2_descriptor which basically just wrap the grib1/grib2 kgds and igdtmpl arrays. To create a descriptor you call init_descriptor and pass the grib arrays.

Then, there are the grid classes which are created by calling init_grid with an ip_grid_descriptor object. At the top-level there's the abstract type ip_grid which contains common variables for each grid like im, jm, rerth, etc. Each concrete subclass of ip_grid (e.g. ip_equid_cylind_grid) contains grid-specific variables and must implement the routines init_grib1, init_grib2 and gdswzd. The init routines are responsible for handling grib1/grib2 specific initializations and setting the proper variables.

The top-level ipolates(v) routines retain their interfaces and then create descriptors, which are used to create the grid object, which is then passed to the actual interpolation routines.

With that framework, all interpolation routines and gdswzd routines now take grid objects instead of passing the descriptors around.

Spectral interpolation has not been merged together yet, and I just copied the grib1/grib2 routines into the module with an interface. There are a lot of specific checks using the grib descriptor arrays that are difficult to generalize. I'm still looking into this.

I am working on cleaning up the code, adding documentation, and doing some final touches.

[GeorgeGayno-NOAA]

@kgerheiser When you are ready, please compile it on one of our officially supported machines. I would like to test it with UFS_UTILS.

[kgerheiser]

@GeorgeGayno-NOAA

I installed in on Hera.

Set export ip_ROOT=/home/Kyle.Gerheiser/NCEPLIBS-ip2/build/install to use it with CMake in UFS_UTILS.

You'll have to add use ip_mod wherever you use the library. That module contains all the public facing routines: ipolates, ipolatev, and gdswzd.

It's possible you might also have to adjust your calls because now that everything is in a module you have to match the interface exactly. Like, in the test code ibi was being passed as a scalar when the interface says it's supposed to be an array. So, I changed it to an array of size one. Or rlat/rlon were declared as 2d arrays, but the polate routines take 1-d arrays. Don't know if you do anything like that.

And the C routines were changed. C doesn't have function overloading so the grib1/grib2 C versions of gdswzd can't both be called gdswzd. I left the ip2 version alone and named the ip1 version gdswzd_grib1.

[kgerheiser]

Building UFS_UTILS I ran into what I thought might happen. In snowdat.f in UFS_UTILS rlat and rlon are 2-d arrays, but the interface wants a 1-d array. I think this is because of the interface gdswzd and it can't find a matching interface that takes 2-d arrays.

I think this can be fixed by adding a gdswzd_2d_array_grib1 in the same way there's a grib2 gdswzd_2d_array in ip2

[kgerheiser]

I was able to get UFS_UTILS to build. You can the changes I had to make:

https://github.com/kgerheiser/UFS_UTILS/commit/a7a8d84879eeec850dafaebc7f295de6c38b893b

There's a lot of inconsistency with the shape of various arrays (sometimes 1-d, sometimes 2-d). Using an interface is a no-go because of that, and I didn't want to change anything. So, I used the specific call directly ipolates_grib1 because that allows you to pass 1-d arrays to 2-d arguments or vice-versa.

I also had to change some intent(in) because in ip they're intent(inout) and you can't pass an intent(in) argument to a routine that's intent(inout).

[GeorgeGayno-NOAA]

I was able to get UFS_UTILS to build. You can the changes I had to make:

kgerheiser/UFS_UTILS@a7a8d84

There's a lot of inconsistency with the shape of various arrays (sometimes 1-d, sometimes 2-d). Using an interface is a no-go because of that, and I didn't want to change anything. So, I used the specific call directly ipolates_grib1 because that allows you to pass 1-d arrays to 2-d arguments or vice-versa.

I also had to change some intent(in) because in ip they're intent(inout) and you can't pass an intent(in) argument to a routine that's intent(inout).

Thanks. I will open an issue and create a branch with your updates.

[GeorgeGayno-NOAA]

@GeorgeGayno-NOAA

I installed in on Hera.

Set export ip_ROOT=/home/Kyle.Gerheiser/NCEPLIBS-ip2/build/install to use it with CMake in UFS_UTILS.

You'll have to add use ip_mod wherever you use the library. That module contains all the public facing routines: ipolates, ipolatev, and gdswzd.

It's possible you might also have to adjust your calls because now that everything is in a module you have to match the interface exactly. Like, in the test code ibi was being passed as a scalar when the interface says it's supposed to be an array. So, I changed it to an array of size one. Or rlat/rlon were declared as 2d arrays, but the polate routines take 1-d arrays. Don't know if you do anything like that.

And the C routines were changed. C doesn't have function overloading so the grib1/grib2 C versions of gdswzd can't both be called gdswzd. I left the ip2 version alone and named the ip1 version gdswzd_grib1.

@kgerheiser Can you open up the permissions on your build directory?

[kgerheiser]

Try it now

[GeorgeGayno-NOAA]

The UFS_UTILS snow2mdl regression test is failing with a seg fault. It happens in routine snowdat.f at line 1257. I tried two ways to call gdswzd.

call gdswzd_grib1(kgds,(-1),npts,fill,xpts,ypts,rlon,rlat,nret) and call gdswzd(kgds,(-1),npts,fill,xpts,ypts,rlon,rlat,nret)

In both cases, the routine returns nret=0. And xpts(1) and ypts(1) = 9999. The input arrays, kgds, rlat, rlon are correct. You can see the log files on Hera: /scratch2/NCEPDEV/stmp1/George.Gayno/reg_tests.snow2mdl.seg.fault.

[kgerheiser]

What kind of grid is it running on?

How can I run just that test and do it interactively?

[kgerheiser]

I think I have it working. All I need is that fort.41 file? I turned on debug flags and see an out-of-bounds error.

[kgerheiser]

I think I found the bug (or at least a bug):

In gdswzd_05 there's this line: IF(H.EQ.-1)ORIENT=ORIENT+180.

But in ip2 there isn't, and I guess I missed it when I was comparing the two files.

[kgerheiser]

I ran the tests myself, and it still doesn't pass, but I plotted it and compared it to the original, and it's close. Here is the relative error between them (array2 - array1) / array1

I'll see if I can find what the issue is.

[kgerheiser]

I am going to use the base UFS_UTILS and change one call at a time while testing at each stage and link in both the new library and the old library to track where the difference comes from.

[GeorgeGayno-NOAA]

@kgerheiser I will try your fix at 9aab5fe.

[edwardhartnett]

Every time you make a fix, you should also ensure there is a test that will fail without that fix.

If you discover some problems with external, manual testing, then you must update the internal, automatic tests, so that they will find the problem next time.

Although it might seem that every bug is unique, careful study has shown that programmers make the same few mistakes over and over again. So putting a test in place every time that happens is a great way to build up a set of test which will be valuable in the future, catching the mistakes immediately.

[kgerheiser]

I can add a new test that uses that particular grid. There's a polar stereo grid test, but it doesn't test the elliptical part.

@GeorgeGayno-NOAA would I be able to edit the output descriptor of the current lat-lon -> polar stereo grid test to and just make it elliptical? Then, I could test it against the existing ip and ip2 library with the same settings.

[GeorgeGayno-NOAA]

I can add a new test that uses that particular grid. There's a polar stereo grid test, but it doesn't test the elliptical part.

@GeorgeGayno-NOAA would I be able to edit the output descriptor of the current lat-lon -> polar stereo grid test to and just make it elliptical? Then, I could test it against the existing ip and ip2 library with the same settings.

@kgerheiser Point me to the source code file for the lat-lon to polar grid test.

[kgerheiser]

Sure, it's the tests you created, but very slightly modified and combined with CMake:

https://github.com/kgerheiser/NCEPLIBS-ip2/blob/feature/combined-ip/tests/interp_mod_grib1.f90

This descriptor specifically (there's also interp_mod_grib2.f90):

https://github.com/kgerheiser/NCEPLIBS-ip2/blob/9aab5fefb03a8dd8994ef8e03ea2107b95cd5706/tests/interp_mod_grib1.f90#L85

Here is where the top-level executable is called:

https://github.com/kgerheiser/NCEPLIBS-ip2/blob/9aab5fefb03a8dd8994ef8e03ea2107b95cd5706/tests/CMakeLists.txt#L38

[GeorgeGayno-NOAA]

@kgerheiser I will try your fix at 9aab5fe.

@kgerheiser I reran my UFS_UTILS snow2mdl regression test. It failed, but it is much closer to the baseline. The depth data is determined from AFWA data using nearest neighbor interpolation. That field appears to match. The snow cover data is determined from IMS data using budget interpolation. That field has slight differences along the snow edge. The differences are slight enough that I think the radius of the earth is no longer the problem. The budget interpolation has lots of options stored in the 'ipopt' array. Can you check to see if the 'ipopt' array is being used correctly. Here is what I use for the IMS interpolation: https://github.com/NOAA-EMC/UFS_UTILS/blob/develop/sorc/emcsfc_snow2mdl.fd/snow2mdl.f#L250 I can help setup the snow test for you if you like.

[kgerheiser]

Yep, I see that something is off in there.

[GeorgeGayno-NOAA]

Sure, it's the tests you created, but very slightly modified and combined with CMake:

https://github.com/kgerheiser/NCEPLIBS-ip2/blob/feature/combined-ip/tests/interp_mod_grib1.f90

This descriptor specifically (there's also interp_mod_grib2.f90):

https://github.com/kgerheiser/NCEPLIBS-ip2/blob/9aab5fefb03a8dd8994ef8e03ea2107b95cd5706/tests/interp_mod_grib1.f90#L85

Here is where the top-level executable is called:

https://github.com/kgerheiser/NCEPLIBS-ip2/blob/9aab5fefb03a8dd8994ef8e03ea2107b95cd5706/tests/CMakeLists.txt#L38

For an elliptical earth grid, the sixth element (GDS Octet 17) of the 'grd212' array should be changed from '8' to '40'. https://www.nco.ncep.noaa.gov/pmb/docs/on388/tabled.html https://www.nco.ncep.noaa.gov/pmb/docs/on388/table7.html

But I would retain the 212 test and add another '212' test that is the same except for an elliptical earth. Give it any number you want.

[kgerheiser]

I found a problem with the negative number sub-grid thing. I was calling gdswzd on the original output grid object (which has a negative grid number and implements gdswzd with an empty call), instead of on the new grid object created with 255 + grid_number.

[GeorgeGayno-NOAA]

I found a problem with the negative number sub-grid thing. I was calling gdswzd on the original output grid object (which has a negative grid number and implements gdswzd with an empty call), instead of on the new grid object created with 255 + grid_number.

Great. I will rerun my test.

[kgerheiser]

I was able to generate test data for a non-spherical Earth pretty easily. Now, I would also like a test for the sub grid option. How would craft the input to do that?

[GeorgeGayno-NOAA]

I was able to generate test data for a non-spherical Earth pretty easily. Now, I would also like a test for the sub grid option. How would craft the input to do that?

You don't modify the input data. Rather you interpolate the input data to a subset of the grid. A subset is just a subsection of the input grid. You can try a single point to start. Allocate output_rlat/rlon, output_bitmap and output_data for one point. Set the lat/lon to a value in the middle of North America. Set 'no' to one. Subtract kgds from 255. Then call ipolates. I think that will work. You can expand the test to do the 'left' half of the grid. You would need to know the lat/lons for each point on the sub grid (read them from a file?) I use this option so I can just pass in the land points as I don't care about snow at water points.

[kgerheiser]

@GeorgeGayno-NOAA could you try this again to double check me? It seems to have worked for me on Hera.

export ip_ROOT=/home/Kyle.Gerheiser/NCEPLIBS-ip2/build/install

[GeorgeGayno-NOAA]

@GeorgeGayno-NOAA could you try this again to double check me? It seems to have worked for me on Hera.

export ip_ROOT=/home/Kyle.Gerheiser/NCEPLIBS-ip2/build/install

Yes, it works now. Thanks.

[edwardhartnett]

Why do we need environment variable ip_ROOT?

[kgerheiser]

You don't, but that's one of the ways CMake finds libraries so George could build UFS_UTILS with my build, the same way as appending to CMAKE_PREFIX_PATH.

[aerorahul]

You don't, but that's one of the ways CMake finds libraries so George could build UFS_UTILS with my build, the same way as appending to CMAKE_PREFIX_PATH.

@edwardhartnett @kgerheiser One could also avoid both pkg_ROOT and CMAKE_PREFIX_PATH by passing -Dpkg_DIR=/path/to/pkg/lib/cmake as arguments to cmake at configuration time. /path/to/pkg/lib/cmake should contain pkg-config.cmake. Ofcourse it is easier to set pkg_ROOT environment variable in the module that loads the pkg.

[edwardhartnett]

Is this ready to merge? Lot's of changes here! ;-)

[kgerheiser]

Well, I broke something yesterday, but otherwise it works (with UFS_UTILS, at least). I was going to finish up converting the documentation to Doxygen. I also think it should be tested with some other application that uses ip (weather model, UPP) before it's merged.

[edwardhartnett]

Can you save the doxygen changes for a separate PR?

[kgerheiser]

Yes, that can be separate.

[edwardhartnett]

@kgerheiser this can be close now, right?