u/v re-mapping consistency

[dgwyther]

Hi folks

I'm looking for a bit of guidance concerning re-mapping between the rho and u/v grids in preprocess/load_ocean_data_roms.py.

I note that the other drivers (e.g. nemo) essentially load LEVOFU etc from the grid:

   # Ocean fraction at each level:
   LEVOFT = maskT*100.0
   LEVOFU = maskU*100.0
   LEVOFV = maskV*100.0

whereas ROMS calculates it quite elaborately:

DEPTHO = xr.DataArray( np.zeros((my,mx))*np.nan, dims=['eta_rho', 'xi_rho'] ) # bathymetry
DEPTHOT = DEPTHO.isel(xi_rho=slice(imin,imax+1),eta_rho=slice(jmin,jmax+1)) # subset bathymetry
DEPTHOU = xr.DataArray( 0.5 * (DEPTHOT.values+DEPTHOT.shift(xi_rho=1).values), dims=['eta_u', 'xi_u'] ) #subset bathymetry *roughly* (not exactly) on the u grid

ROMS has a mask_u variable in the grid file, which is ==1 for ocean and ==0 for land - is this equivalent to maskU in NEMO?

DEPTHOU is used a few lines below:

   LEVOFU = LEVOFU.where( ( (DepTUV==0.) & (DEPFLFU==0.) ) | ( (DEPTHOU>DepTUV) & (DEPFLFU<DepTUV) ), 0.e0 ) \
            * ncM.mask_u.isel(xi_u=slice(imin,imax+1),eta_u=slice(jmin,jmax+1))

I don't really understand what LEVOFU (depth of the u points?) is doing.

The reason I'm raising this possible issue, is that this line:

   LEVOFU = LEVOFU.where( ( (DepTUV==0.) & (DEPFLFU==0.) ) | ( (DEPTHOU>DepTUV) & (DEPFLFU<DepTUV) ), 0.e0 ) \
            * ncM.mask_u.isel(xi_u=slice(imin,imax+1),eta_u=slice(jmin,jmax+1))

can raise an error, depending on the test case. For example, with the ROMS_test case, it runs fine, but for a different grid size (my actual model data), this term ncM.mask_u.isel(xi_u=slice(imin,imax+1),eta_u=slice(jmin,jmax+1) has a different spatial size than DEPTHOU etc, as though the re-mapping to the u-grid by the .shift() command didn't work correctly.

I re-wrote the calculations of the DEPTHOFU and DEPTHOFV to be the averages of the rho points, and then made a new depth vector ztmp for the u, v, and rho grids. This seemed to work, but didn't work with the ROMS_test data.

So just looking for some thoughts on what might be going on. Any advice is appreciated!

[dgwyther]

A little bit more information. In preproc/load_ocean_data_roms.py there are some lines to find the average of the bathymetry and ice draft rho grid:

   DEPTHOU = xr.DataArray( 0.5 * (DEPTHOT.values+DEPTHOT.shift(xi_rho=1).values), dims=['eta_u', 'xi_u'] )
   DEPFLFU = xr.DataArray( 0.5 * (DEPFLFT.values+DEPFLFT.shift(xi_rho=1).values), dims=['eta_u', 'xi_u'] )
   DEPTHOV = xr.DataArray( 0.5 * (DEPTHOT.values+DEPTHOT.shift(eta_rho=1).values), dims=['eta_v', 'xi_v'] )
   DEPFLFV = xr.DataArray( 0.5 * (DEPFLFT.values+DEPFLFT.shift(eta_rho=1).values), dims=['eta_v', 'xi_v'] )

But they're not exactly moving the to the u/v grid, as the .shift() command is shifting in one direction and padding with a NaN value. So the resultant grid doesn't have dims 'eta_u' and 'xi_u', it still has dims 'eta_rho' and 'xi_rho'.

I could replace these lines with something like:

DEPTHOU = xr.DataArray(0.5 * (DEPTHOT.isel(xi_rho=slice(0, -1)) + DEPTHOT.isel(xi_rho=slice(1, None))).values,dims=["eta_u", "xi_u"])

Which would be the correct average from rho -> u points. This works for my setup, but with the test data, where the larger grid is subset, there is a multiple of two different fields:

LEVOFU = LEVOFU.where( ( (DepTUV==0.) & (DEPFLFU==0.) ) | ( (DEPTHOU>DepTUV) & (DEPFLFU<DepTUV) ), 0.e0 ) \
            * ncM.mask_u.isel(xi_u=slice(imin,imax+1),eta_u=slice(jmin,jmax+1))

Here, DEPTHOU and DEPFLFU will have the size of imin:imax, jmin:jmax+1 (as they were calculated by finding the center point by 1 cell in the x direction. Whereas ncM.mask_u.isel(xi_u=slice(imin,imax+1),eta_u=slice(jmin,jmax+1) obviously has the full size of imin:imax+1 and jmin:jmin+1

[dgwyther]

Just an update with more information on this issue.

1. when imin, imax, jmin and jmax are calculated, they're selected based on the lat/lon of the rho grid.

2. the same extents are used to subset the u/v grids.

3. if the entire domain is selected (i.e. if the ROMS domain is == the MISOMIP2 target domain), then the entire rho grid will be chosen by imin, imax, jmin and jmax. However, when these are then used to subset the u/v grids, the code fails, as by definition, the u,v grids are 1 smaller than the rho grid in the x, y directions.

This issue was not detected in test code, as the source ROMS grid was much larger than the MISOMIP2 target grid, and so the subset of the u,v grid didn't ask for indices outside of the source grid.

[dgwyther]

@pmathiot or @nicojourdain? Keen for any suggestions :)

[dgwyther]

I can force a reduction in the limits for u, v subsets:

   imin_u = imin
   imax_u = min(imax, ncM.sizes['xi_u'] - 1)

   jmin_u = jmin
   jmax_u = min(jmax, ncM.sizes['eta_u'] - 1)

   imin_v = imin
   imax_v = min(imax, ncM.sizes['xi_v'] - 1)

   jmin_v = jmin
   jmax_v = min(jmax, ncM.sizes['eta_v'] - 1)

But it is somewhat hiding the issue... And it runs into an issue at the end of the script where everything is being forced into 1-D vector:

   #----------
   # Create new xarray dataset including all useful variables:
   # reshaping (x,y) as 1-dimensional (sxy)

   nxy=(jmax-jmin+1)*(imax-imin+1)

   ds = xr.Dataset(
      {
       "SO":        (["time", "z", "sxy"], np.reshape( SO_z, (mtime,mz,nxy)) ),

Obviously, we can't then have some variables being shorter by 1 in 1 dimension.

So I probably need to update the u, v masks to make them rho shaped..

[pmathiot]

'ROMS has a mask_u variable in the grid file, which is ==1 for ocean and ==0 for land - is this equivalent to maskU in NEMO?' yes

[dgwyther]

Thanks @pmathiot . If NEMO is an Arakawa-C grid model, why are the u, v, grid quantities (e.g. u-velocity, lat_u, etc) still using the same coordinates "x" and "y" as the rho grid coordinates in the preproc/load_ocean_data_nemo.py:

   # ocean y-ward velocity [m s-1]
   if ( "voce" in ncV.data_vars ):
     VY = ncV.voce
   elif ( "vomecrty" in ncV.data_vars ):
     VY = ncV.vomecrty
   elif ( "vo" in ncV.data_vars ):
     VY = ncV.vo
   else:
     print('    WARNING :   No data found for VY  -->  filled with NaNs')
     VY = xr.DataArray( np.zeros((mtime,mz,my,mx))*np.nan, dims=['time', 'z', 'y', 'x'] )

or

  ds = xr.Dataset(
      {
...
       "UX":        (["time", "z", "sxy"], np.reshape( UX.isel(x=slice(imin,imax+1),y=slice(jmin,jmax+1)).values, (mtime,mz,nxy)) ),
       "VY":        (["time", "z", "sxy"], np.reshape( VY.isel(x=slice(imin,imax+1),y=slice(jmin,jmax+1)).values, (mtime,mz,nxy)) ),

Noting that UX and VY are indexed with the x and y coordinate which is the same for THETA0 etc.

I'm asking because the load_ocean_data_roms.py script is substantially more complex around remapping of u- and v-grid quantities. I think I will have to re-write it to do a re-mapping to rho points when the fields are initially loaded, and essentially put everything on the rho grid.

[nicojourdain]

Hi Dave,

As I said previously, I don't know ROMS, I just quickly coded an example of some outputs that I got from Ole because people were asking that we share this kind of tools. Your ROMS output may have a different format that is not compatible with what I did.

The outputs in NEMO have the same indices (i,j) for the rho, u and v grids, even if they do not correspond to the same longitude and latitude.

The imin, imax, jmin, jmax are just there to reduce the size of the model output, it is not something very accurate. Then the interpolation will take the closest points using the correct coordinates. In the provided code, we take more than the MISOMIP2 domain bounds, so that the interpolation is fine:

   [lonmin,lonmax,latmin,latmax] = grid_bounds_oce(region=region)
   wdeg = 1.5 * np.max([dxT.where(latT <= latmax).max(),dyT.where(latT <= latmax).max()]) / 6.37e6 * 180. / np.pi / np.cos(0.5*(latmin+latmax)*np.pi/180.)
   lonmin=lonmin - wdeg # take a bit more for interpolation
   lonmax=lonmax + wdeg
   latmin=latmin - wdeg
   latmax=latmax + wdeg

If you need one more point in ROMS, you can increase 1.5 to 2.5 or calculate specific imin/imax for the rho, u and v grids as you are suggesting.

[dgwyther]

Hi @nicojourdain

Thanks for the information. That makes sense.

The issue arises because the ROMS u, v grids have different sizes to the rho grid.

imin, imax, jmin, jmax

are determined on the rho grid, but are then reused directly to subset variables on the u, v grids:

UX = UX.isel(xi_rho=slice(imin,imax+1),eta_rho=slice(jmin,jmax+1))
VY = VY.isel(xi_rho=slice(imin,imax+1),eta_rho=slice(jmin,jmax+1))

which is only safe if the selected rho subset does not touch the domain boundary.

For example, suppose the full ROMS domain is:

rho: 0 ... 501
u:   0 ... 500

and MISOMIP2 only needs:
rho: 50 ... 450
Then:
slice(50,451)
is valid on both rho and u grids, so the bug remains hidden.

However, when the entire domain is selected:

imin = 0
imax = 501

the rho subset:
slice(0,502)
is fine. But then:
UX.isel(xi_u=slice(0,502))
is requested. xarray quietly truncates that to the available range:
xi_u = 0:501
giving 501 points.
When the quantities are combined e.g.

LEVOFU = LEVOFU.where( ( (DepTUV==0.) & (DEPFLFU==0.) ) | ( (DEPTHOU>DepTUV) & (DEPFLFU<DepTUV) ), 0.e0 ) \
            * mask_u.isel(xi_rho=slice(imin,imax+1),eta_rho=slice(jmin,jmax+1))

There will be a multiplication of a 501-long dimension by a 502-long dimension.

I have 'fixed' it, by first interpolating all u and v quantities to the rho grid at the start. That way, the resultant fields are guaranteed to have shape of xi_rho * eta_rho. The downside with this is that the fields are interpolated, so they may not be exact.

[nicojourdain]

Ok, I see. Glad that you found a way around it.

We are not planning to use the interpolated fields to calculate heat or salt budgets, so I don't think it is an issue that there are two interpolations.

On a more general note, it would be good that you create a branch with your modified code, so that we can pull it to the main repo later on.

[dgwyther]

Yep, I'm putting my code together as cleanly as I can, then I'll make a pull request.