Clarification on setting up MOM5 velocity data for 3D advection

[hrsdawson]

Hi parcels team. I’m trying to compute 3D advection with B-grid MOM5 velocities. I've seen some discussion of B-grid interpolation here and in this documentation example, however I’m still a bit confused about the vertical dimension setup and am looking for some clarification.

In the MOM5 data, u, v, and w have the following dimensions:

u (st_ocean, xu_ocean, yu_ocean) 
v (st_ocean, xu_ocean, yu_ocean)
w (sw_ocean, xu_ocean, yu_ocean)

The vertical velocity w and the vertical coordinate for the velocity cells (sw_ocean) are defined on the bottom interface of the tracer cell such that the uppermost w value is located at the bottom of the surface layer and there is no w data for the surface at 0 m. The sw_ocean coordinate starts at 1.08 m and extends down to ~5800 m (increasing down). From the documentation example and this section of the code, I understand that parcels wants the u, v and w velocities defined at the same point in the top corner of the grid cells. Because of this I modified sw_ocean by shifting it up one cell with a zero at the surface so that sw_ocean_mod now extends from 0 to ~5600 m. I’ve then been setting up parcels as follows:

ufiles = 'datapath/ocean_ufile.nc'
vfiles = 'datapath/ocean_vfile.nc'
wfiles = 'datapath/ocean_wfile.nc'
mesh_mask = 'datapath/Southern_Ocean_coordinates.nc'

filenames = {'U': {'lon': mesh_mask, 'lat': mesh_mask, 'depth': mesh_mask, 'data': ufiles}, 
              'V': {'lon': mesh_mask, 'lat': mesh_mask, 'depth': mesh_mask, 'data': vfiles}, 
              'W': {'lon': mesh_mask, 'lat': mesh_mask, 'depth': mesh_mask, 'data': wfiles}}

variables = {'U': 'u',
             'V': 'v',
             'W': 'wt'} 
dimensions = {'U': {'lon': 'geolon_c', 'lat':'geolat_c', 'depth':'sw_ocean_mod', 'time':'time'}, 
             'V': {'lon': 'geolon_c', 'lat': 'geolat_c', 'depth': 'sw_ocean_mod', 'time': 'time'},
             'W': {'lon': 'geolon_c', 'lat': 'geolat_c', 'depth': 'sw_ocean_mod', 'time': 'time'}} 

fieldset = FieldSet.from_b_grid_dataset(filenames, variables, dimensions, field_chunksize=(1,7,300,400)) 
fieldset.W.set_scaling_factor(-1)  # set scaling factor for vertical velocities to -1

Is this the correct way to assign the vertical coordinate for these MOM5 b-grid velocities? I’m currently having an issue with particles getting stuck at the bottom of the ocean prompting me to think there may be something wrong with my set up. This also happens if I use the original sw_ocean coordinate. Should the vertical coordinate have a value at both the ocean surface and the ocean floor? In which case the dimensions of sw_ocean would not be equal to the vertical dimension of w, and I would perhaps need to add a layer of zeros to w inside the parcels code at the surface – something like this for pop interpolation.

Ideally I want to give parcels a 3D array of depth to account for bottom shaved cells used in the model, but I first want to understand the correct set up using a Rectilinear Z grid. Any help is greatly appreciated! I’m happy to share some code and/or data if that assists at all.

[erikvansebille]

Thanks @hrsdawson, for this elaborate and very detailed comment/question! You clearly really dug into this quite a bit already. I'm not really sure how much I can help you further, our grid interpolation guru, @delandmeterp, has taken up another position so is not officially part of the Parcels team anymore. I assume you've also looked at section 2.1.3 of his GMD paper, which specifically deals with 3D B-grids?

That said, here are some thoughts:

1. Assuming that W is indeed zero at the ocean surface in MOM, it could make sense to add a layer of zeros at the surface. You could either adapt the Parcels code itself, but that is somewhat cumbersome; alternatively, it might be possible to first open the data with ds=xarray.open_mfdataset(...), add a layer of zeros in this xarray DataSet (?), and then convert the DataSet to a FieldSet using FieldSet.from_xarray_dataset(ds, ...)
2. If you follow step 1 above, be careful what to do with U and V too. Does the first level of these then correspond to the depthw=0? You'll need to think about that
3. Shaved cell support is indeed something that has been on my wish-list for a long time! Let me know how this goes

Perhaps it would be good to engage some of the MOM developers too, to further discuss them with this? @StephenGriffies, would you happen to know who we could ask for advice in figuring out using Parcels with MOM5 fields?

[hrsdawson]

Hi @erikvansebille, thanks for the quick reply and the suggestions. Looks like I might need to do some more digging. Yep I have seen that section of the paper. I should have also mentioned above that some of my confusion is arising from this contributed OFAM example, which has the same dimensions as above, but uses the vertical tracer cell coordinate (st_ocean) for the depth dimension, rather than the vertical velocity coordinate (sw_ocean). I can see that a line setting up dimensions with sw_ocean has been commented out. Should I be using the vertical tracer coordinate as in this OFAM example?

[erikvansebille]

Ah, that may be a bug in https://github.com/OceanParcels/parcels_contributions/blob/master/Examples/example_OFAM3.py. @CKehl, do you remember why you used 'st_ocean', rather than 'sw_ocean'?

[StephenGriffies]

Hi,

Here are some comments.

1/ wt is not zero at the top model grid cell in MOM5. The model uses an explicit free surface that is fully nonlinear. And when using z* vertical coordinates, all grid cells have a time dependent thickness.

2/ The wt field in MOM5 is registered with the following axes

id_wt = register_diag_field ('ocean_model', 'wt', Grd%vel_axes_wt(1:3), Time%model_time, &
'dia-surface velocity T-points', 'm/sec', missing_value=missing_value, range=(/-10.e4,10.e4/))

where vel_axes_wt places the wt velocity at the center of the tracer cell top/bottom. I am unsure of the st_ocean and sw_ocean grid definitions you refer to. These are static. The correct static analog to the wt grid is indeed the sw_ocean grid.

3/ Veronica Tamsitt v.tamsitt@unsw.edu.au is actively using Parcels with MOM5 output. I suggest connecting with her for further details.

[adele-morrison]

Thanks Erik and Steve. @hrsdawson is doing her PhD with Matt England, myself and @vtamsitt. Veronica is still using CMS with MOM5, not Parcels yet. We're trying to switch everyone over to Parcels, starting with @hrsdawson.

A couple of points in reply to the above:

1/ To start with we are trying to get Parcels working with the approximation of stationary z layers, rather than z*. This is what we did for CM2.6 with CMS @StephenGriffies. Do we know how large the errors are for this approximation? I guess we will also have problems with particles trying to leave through the surface due to this approximation. We could in time work on setting up Parcels to represent the time-varying z* layers using the Curvilinear S Grid implementation, if we think it is worth it.

2/ Regarding what w values to insert at the surface, perhaps the best option would be to make a copy of w from the bottom of the top z layer (i.e. the uppermost w that we have saved), rather than zeros?
Erik, I don't think your point 2. above about matching the correct u/v to this new surface w is an issue. As I understand it for b-grids in Parcels, it is expected that there will be a w value at the uppermost surface (i.e. half a grid cell above u/v). But we are currently missing that uppermost w value.

[StephenGriffies]

Thanks @adele157 for clarification: My mistake forgetting that CMS was still in use.

I also wish to clarify my comment 2/ above. The wt velocity is indeed defined at sw_ocean. This setting is specified in

https://github.com/mom-ocean/MOM5/blob/master/src/mom5/ocean_core/ocean_grids.F90

which is the latest MOM5 Github repo. I see no substantial changes in that file since 2014, though not sure if there was a change made before that effects this discussion.

[StephenGriffies]

In regards to @adele157 comment about stationary z-levels rather than z*. I am not sure what is meant. Is it related to the following?

Namely, the wt diagnosed in a z* configuration of MOM5 measures the transport through the z* coordinate surface. This transport is distinct (but only by a bit) from the transport through a z-coordinate surface.

Let's write wt* for the vertical velocity diagnosed in z*-coordinate configuration and write wt for the transport diagnosed in a z-coordinate configuration. Are you concerned about the difference between wt* and wt for use in particles? That difference is very small in general, and it is a difference I never bothered with. But do you think it is a difference that has qualitative impact on particle trajectories? Or perhaps there is an issue at the boundaries?

[erikvansebille]

Hi @adele157

1/ To start with we are trying to get Parcels working with the approximation of stationary z layers, rather than z*. This is what we did for CM2.6 with CMS @StephenGriffies. Do we know how large the errors are for this approximation?

There's two types of errors here. The one that @StephenGriffies described above (the difference between wt* and wt), and another that the interpolation of U and V depends on the particle's position in the grid, so that if the depth layers move up and down, the particle's relative position in the grid moves too, even if the particle's depth doesn't change. But again, I expect that this difference is very small?

In any case, it might be a good idea at some point to compare particle trajectories assuming a fixed depth grid to particle trajectories using the time-varying cS-grid implementation. I'd be keen to see how big the difference is!

[CKehl]

Ah, that may be a bug in https://github.com/OceanParcels/parcels_contributions/blob/master/Examples/example_OFAM3.py. @CKehl, do you remember why you used 'st_ocean', rather than 'sw_ocean'?

That was based (a) of the use case we chose at that point and (b) of the variable declarations in the NetCDF file. It can be chunked either way, it just depends on how the dimensions are linked as variable parameters in NetCDF.

[eavellashaw]

Hi everyone!
I add something here as it might be interesting for the MOM5 community.
We found that when creating the fieldset with from_mom5(), it was not accounting for the grid rotation (especially important in the Arctic). We found a way to solve that by modifying the Kernel AdvectionRK4_4D to include the curvilinearity of the grid. This has been solved in #1509, don't hesitate to check it out if you are working in the Arctic!