NLON, ELAT not computed when TLAT, TLON, ANGLET on grid file
Closed this issue · 14 comments
This was raised in the CICE forum. We should fix this before the next release. I can look into it.
Dave
How should we fix this? Should we NOT write the N/E coords with B grid or should we initialize N/E coords with B grid? I can create a PR unless someone else wants to take it. Will try to fix today.
I think I'll go with not writing the N/E coords with B grid as suggested by @TillRasmussen. Let me know if others prefer another solution.
Why aren't our tests catching this?? Do we need to add some tests or some other kind of IC file to our test suite? @phil-blain?
I think I'll go with not writing the N/E coords with B grid as suggested by @TillRasmussen. Let me know if others prefer another solution.
I agree with this approach. It also makes the history files a little less voluminous, which is good.
Why aren't our tests catching this?? Do we need to add some tests or some other kind of IC file to our test suite? @phil-blain?
As I wrote in the forum, I think this is because all our grids have only ULON, ULAT and ANGLE, and so the corresponding TLON, TLAT and ANGLET are computed by the code (in ice_grid::Tlatlon).
If the grid has ANGLET, then the code assumes that it also has TLON, TLAT and Tlatlon is not called. Since the N- and E- locations lat/lon arrays are also computed in that subroutine, then they are not set before the history code writes them to the history files.
What if you have an input grid that has TLON, TLAT, and ANGLET but then you are running on the C-grid. I think this is not a B grid issue. I think we have to compute the E/N coords if ANGLET is included on the grid file. Am I missing something?
This is my understanding as well. I would retitle the issue
NLON, ELAT, etc. not computed under l_readCenter=.true.
We should add a grid that has TLON, TLAT, ANGLET to our test suite. Anyone have a grid we could use? If someone can provide a file, I'll add some test cases. We could also create one from the gx3, maybe I'll try to do that.
I think it should be easy to create one for gx3 from the history output.
Thanks @phil-blain, that's what I'm trying to do now.
If I understand it correct then latitudes are used for Coriolis calulations and longitudes/latitudes at velocity grid are used for netcdf coordinates. Center points (Tlon, Tlat) are used for coupling.
Therefore would it be more correct to read in ulon,ulat on B grid and never calculate nlon, nlat, elat and elon. Then on c and cd grids one would read in nlon, nlat, elat and elon and never calculate b grid.
Then skip the calculations of lon/lat at the velocity points in the TLATLON function.
Then some logic is required when the coordinates for the netcdf file is created.
I think I introduced the NLON, ELAT etc. when we first started with the CD-grid implementation. I think it might be good to have these for generality. I believe currently the binary grid files have ULON, ULAT, ANGLE(U) and then TLON, TLAT, and ANGLET are computed internally. However, I think the netCDF grid files also have the T location of these as well. Honestly, we are moving to the MOM6 model and it uses a "supergrid" with all of this on a higher resolution grid. I need to implement something in CICE that understands this.
The binary files are their own thing. In our test suites, we have only the gx3, gx1, and tx1 binary files. I am going to create a gx3 netcdf file with the T data on that file, so we can test the code.
Regarding @TillRasmussen's comments. On the B-grid, we read U grid variables and compute T grid variables if they are not on the grid file. We don't need to compute N, E grid variables for B, but might anyway. On the C/CD grid, we read in the same (U grid based) file. Yes, this may not be ideal but allows users to use the same and existing grid file for C-grid runs. That is convenient. It might make sense to add a new grid format that would have the N and E grids on them, but we'd have to think about what is best and I think we still need to provide ability to read the current grid files and run C-grid. It would surprise me if the C-grid never needs to know the T lon and lat. You could be right that the physics doesn't care. It's an interesting question.
Ok, one more place. tlon, tlat are also used for the zenith angle of the sun
OK, I spent several hours this afternoon trying to clean up the CICE N/E grid stuff. There are many challenges.
First, the compilers I'm testing on don't abort in the default case when NLON is allocated but not initialized then written. On my compiler, it writes 0s. So, on other compilers or with other compiler flags, this could abort, but it won't all the time. We still have a task to debug with Nan initialization, but we're not there yet. If we had that, it probably would fail like the user.
So, for our testing to truly trap the error, the best thing to do is NOT allocate the N/E variables at all. Also, it's not just NLON, it's N/E LON, LAT, area, mask, and so forth. The fact this failed on NLON for the user is just the tip of the iceberg. OK, maybe a small iceberg, but still...
So I moved all the N/E grid variable allocation into an if "C-grid" block. Then I had to add a bunch of if statements in the ice_grid initialization. Then I had to modify ice_history_write.F90 (note that there are two of them to do netcdf and pio) to provide increased flexibility in the implementation to write only T and U or to write T, U, N, and E. As implemented, somethings are hardwired with parameter statements to write all 4 T, U, N, and E locations, so that had to be refactored. And ice_history_write is one of the files with the worst indentation, so I fixed that too because I was struggling to read the code.
Then, I finally got to a point where it got thru initialization and initial history writing and I then realize uarea and narea are needed by the dynamics for locate triangles. I think this has nothing to do with the C-grid. I suspect some old code that was doing the earea/narea computation locally before the C-grid was added was changed to use the new C-grid variables for convenience/performance. So now I'm adding back in earea and narea to the B-grid implementation and that also requires that d[x,y][E,N] also exist. So now I'm adding back in a some of the N/E grid arrays that I removed before. But I also don't love that we allocate some N/E grid variables but not others. This seems confusing to me. I prefer all or none.
So, circling back again. Maybe I'll leave all the grid N/E allocation and initialization in for B-grid, but try to just remove the N/E variables in the history file (is this really worth it). Again, the problem is that implementation will not help us trap this error in the future. And I cannot verify that the fixes I implement actually resolve this issue completely (I have to NOT allocate variables to do that). And I do wonder if some people might actually want some of the N/E grid written to the history file as part of the B-grid runs. What we'd be doing is taking it out without an easy option to turn it back on.
OK, I think I might have just talked myself into leaving it as is. N and E grid variables are allocated and initialized even when just running the B-grid. And those variables are going to be written to the history file by default for now. OK?