Getting NaN valus in initial and boundary files

[09CL6010]

Hi...

I am trying to create initial and boundary netcdf files from daily data for interannual simulation in ROMS. For this purpose I have created initial and boundary file from Mercator Ocean data using my grid file. But after interpolation from Mercator data I found some NaN values (for temperature, salinity) exist in both initial and boundary files. For simulation there should not be any NaN values in those files. Is this a problem of grid I have generated or a problem of interpolation from Mercator data??

[kate]

Could be either. How did you go about building the grid and how did you go about doing the interpolation? I beat on the pyroms scripts until I could interpolate an initial file from Mercator and even then it had NaN at the North Pole which I just made go away.

[09CL6010]

Dear Kate, Thanks a lot for the reply.
I have created grid file using Seagrid from ROMS matlab toolbox for my simulation in North Indian Ocean and used etopo1 as bathymetry data. For initial and boundary condition files, I have used the matlab codes for Mercator data available in ROMS matlab toolbox. The codes are running successfully but few NaNs are exist always. Should I remove NaNs with some values (like 0)??

[kate]

Do you know where the NaNs are and why they should be there? Are they on land or at sea? At all depths? For my North Pole problem I didn't want to replace NaN by zero because that does weird things to the density field via salinity. I simply averaged some neighboring values together and replaced the two NaN points by that. You need to figure out what would be sensible in your case.

[wilkin]

Mercator-Ocean data are reported with netcdf _FillValue assigned to the land area. When interpolating these data to a ROMS grid, points on your ROMS grid near the coast can fall into cells in Mercator-Ocean that are enclosed by at least one _FillValue - which Matlab will interpret as NaN - and propagate the _FillValue to that cell in ROMS.

If you are using Matlab, then you should read the 'help' on griddedInterpolant and scatteredInterpolant because options to extrapolate outside the convex hull of the data can help you avoid some of these issues, for example by invoking a 'nearest' neighbor extrapolation. To use scatteredInterpolant you would want to exclude all Mercator-Ocean land data from the input.

However, an island in Mercator-Ocean will not be a region outside the convex hull of the data and you need to be wary of interpolating across an island - and potentially across quite different water masses. I made this mistake with Mercator-Ocean output in the Philippine archipelago and got some pretty crazy values in some deep basins.]

Inevitably, you are going to have to carefully locate the values in your interpolated ROMS fields - using the 'find(isnan(....' function in Matlab - and then make intelligent choices on how to fill those values from the data at hand.

Be aware that Mercator-Ocean data are reported on standard z-levels, and therefore have a different land mask at every depth.

[kate]

If the NaNs are all over at the land-sea edge, I would see that as something that needs to be fixed by the tool-builders. We do a "flood-fill", extrapolating ocean values into the land for each z-level before doing any kind of interpolation. It's such a pervasive problem that Raphael even does it for the winds and other atmospheric fields.

[wilkin]

The "flood-fill" into land is what I am meaning by careful extrapolation, but in my experience is it tricky to build a very robust tool for this. Caution is always required.

For me, in the Philippine archipelago a lateral flood fill caused absurd values in the deep Sulu Sea because my deepest ROMS point was deeper than the deepest valid Mercator-Ocean point (I had more detailed bathymetry) and the flood fill propagated across deep T/S values from the Celebes Sea, which is a completely different water mass. In this instance I had to concoct a tailored combination of vertical extrapolation and lateral extrapolation that respected the separation of water bodies by peninsulas and submarine sills.

If you have a robust code set for doing this it would be great to share it here. This is a problem we all grapple with in down-scaling basin/global models to regional applications.

[kate]

I'll attach a Fortran code I got from Raphael to do the flooding. We compile to a shared library, then call this from Python. I've used it for HYCOM files. It looks like I copied from an older set of files for the Mercator flooding - that's all in Python.