Add a couple of diagnostic figures at the end of the function for sanity checking the function. Also change the way the open boundary nodes are index by using the same approach as is in the code to generate the Casename_obc.dat file. This should avoid problems where the interpolated arrays store the data in a different order from how FVCOM is expecting it (I am assuming here that FVCOM assumes the order of boundary nodes specified in Casename_obc.dat will be common across all input data which use the open boundary nodes. This seems sensible to me)
Change the way the open boundary nodes are obtained to match that which is in the code which generates Casename_obc.dat. This should stop problems with ordering values in the input arrays where the manner in which the list of open boundary nodes are generated in different ways resulting in different output orders
Flip the order of the vertical levels in the POLCOMS data to match its depth array values (i.e. go from the seabed to the surface instead of the other way around). This changes means the POLCOMS data matches the convention in the FVCOM arrays (where the surface is the first index) and in the POLCOMS depth data. Also change the 1D interpolation to use the pchip interpolation method (shape-preserving piecewise cubic interpolation) instead of linear interpolation. The original change from pchip to linear was because I thought pchip didn't do extrapolation, but since we're scaling the POLCOMS vertical profiles to the FVCOM depths, it doesn't matter. Finally, some miscellaneous changes (moving the creation of the open boundary nodes array to outside the loops and tidying up some of the whitespace/figures code