Discrepancies / bugs in grid generation

[NoraLoose]

As I am trying to understand and cleaning up the code in grid.py (as well as writing tests for it), I found some discrepancies between the original matlab code and the python code, specifically within the function _make_initial_lon_lat_ds:

1. The matlab code generates a domain of (nx+1, ny+1) points, whereas the python code a domain of (nx+2, ny+2) points. Is this intended? @sdbachman @TomNicholas

2. This loop will do exactly the same thing 100 times: https://github.com/CWorthy-ocean/roms-tools/blob/ba7762b14f7edbcf4d9a27b4727bcda76b95faa7/roms_tools/setup/grid.py#L240-L269

because in line 269 it updates latitude_array_1d_in_degrees rather than domain_width_in_degrees_latitude, where the latter is what is (correctly) updated by the corresponding matlab code:

dlat = width/r_earth;
for it = 1:100
    y1 = log(tan(pi/4-dlat/4));
    y2 = log(tan(pi/4+dlat/4));
    y = (y2-y1)*[-0.5:1:nw+0.5]/nw + y1;  
    lat1d = atan(sinh(y));
    dlat_cen = 0.5*(lat1d(round(nw/2)+1)-lat1d(round(nw/2)-1));
    dlon_cen = dlon/nl;
    mul = dlat_cen/dlon_cen*length/width*nw/nl;
    dlat = dlat/mul;
end

I will put in a PR to fix the issue described in 2.

[TomNicholas]

Neither of those seem deliberate! Good catches @NoraLoose

[NoraLoose]

Sure that the first one is not deliberate, though? The README (documentation) shows

grid = Grid(
    nx=100,          # number of points in the x-direction (not including 2 boundary cells on either end)
    ny=100,          # number of points in the y-direction (not including 2 boundary cells on either end)
    size_x=1800,     # size of the domain in the x-direction (in km)
    size_y=2400,     # size of the domain in the y-direction (in km)
    center_lon=-21,  # longitude of the center of the domain
    center_lat=61,   # latitude of the center of the domain
    rot=20,          # rotation of the grid's x-direction from lines of constant longitude, with positive values being a counter-clockwise rotation
)

so it mentions 2 (not 1) boundary cells on either end.

[TomNicholas]

1 boundary cell on either end would give nx+2 no? So the README of the python code is currently inconsistent with the current behaviour of the python code?

[NoraLoose]

The README of the python code is consistent with the current behavior of the python code. However, the behavior of the python code is inconsistent with the behavior of the original matlab code. But this may be intentional?

[TomNicholas]

The correct behaviour is to generate a grid which matches whatever ROMS writes out. If you look at some ROMS output files you'll see they have larger chunks at the end, like {10, 8, 8, 8, 10} etc. or {9, 8, 8, 8, 9}. I can't remember which of these two it is, but by looking at some output files you should be able to see.

[NoraLoose]

I am confused now. I thought the purpose of our python code is to generate a netcdf file that serves as an input file for ROMS. Or does the ROMS fortran code not require a netcdf file with the following info

but is happy with just inputs like

    nx=100,          # number of points in the x-direction (not including 2 boundary cells on either end)
    ny=100,          # number of points in the y-direction (not including 2 boundary cells on either end)
    size_x=1800,     # size of the domain in the x-direction (in km)
    size_y=2400,     # size of the domain in the y-direction (in km)
    center_lon=-21,  # longitude of the center of the domain
    center_lat=61,   # latitude of the center of the domain
    rot=20,          # rotation of the grid's x-direction from lines of constant longitude, with positive values being a counter-clockwise rotation

?

[TomNicholas]

I thought the purpose of our python code is to generate a netcdf file that serves as an input file for ROMS.

It is. But IIRC that input grid file has a similar structure to the output data files that ROMS will write out. So I'm suggesting looking at examples of those files to deduce what exactly the files written by the python script need to look like.

[NoraLoose]

Closing this because #19 fixed this