Add curl computation function

[bradyrx]

I think it would be nice to have a function that computed wind stress curl (or perhaps vector curl agnostic to the variable) over the POP grid.

In past years, I had a lot of trouble getting this done effectively in python, and was under the impression that curl/divergence computations in python weren't supported too well. I turned to NCL for computing it (https://www.ncl.ucar.edu/Document/Functions/Built-in/uv2vrF-1.shtml). However, that function depends on a fixed (non-POP) grid.

Here's a reference that might be helpful for python: https://auckland.figshare.com/articles/Curl_of_a_vector_field/5768280/1.

cc @matt-long @andersy005. Any thoughts? Is this better for pop-tools or esmlab? If the latter, this should be planned to handle all grid types. If the former, just POP. I'd be happy to lead this effort with some feedback, but #12 is higher priority when I can put it to the top of the docket.

This of course should be straight forward following:

[matt-long]

I think xgcm has a curl operator.

[bradyrx]

It might be nice then to add a utility here that converts raw POP output to an xGCM object, or xGCM compatible object, as we discussed before. Ryan Abernathy added a code snippet for this in #12, and it shows that it's not straight forward to a new user.

[matt-long]

Agreed. I think a pop2xgcm_grid method would be great. I could also see adding making xgcm a dependency and building documentation about how to do things with it.

[andersy005]

Agreed. I think a pop2xgcm_grid method would be great. I could also see adding making xgcm a dependency and building documentation about how to do things with it.

Now that we have to_xgcm_grid_dataset(), it would be nice to revisit this issue. I recently received an email asking me how to compute gradient, curl and Laplacian on POP original grids via pop-tools, and I didn't have any examples to point to.

[bradyrx]

Yes this should be fairly trivial with the to_xgcm_grid_dataset(). I'd be interested to see the implementation. That'll help guide #12 once I get to it at the end of the semester. :) I think now that the xgcm infrastructure is here I can go with @rabernat's suggestion to use xgcm there. I'll probably do it with and without xgcm for a sanity check on answers.

[matt-long]

@dcherian may have perspective on how to move forward here.

[DapengLi17]

Hi All:
I am trying to compute wind stress curl over the POP original grids. After using to_xgcm_grid_dataset() to make CESM POP grids output compatible to xgcm, I still have difficulty using the xgcm gradient operator. Both TAUX and TAUY are on the Ugrids. The dimensions of dTAUX/dy and dTAUY/dx are ('nlat_t', 'nlon_u') and ('nlat_u', 'nlon_t'). When you subtract the two, the dimension becomes ('nlat_t', 'nlon_u', 'nlat_u', 'nlon_t'), 4D array. Please see the attached codes. I used jupyter notebook and the GitHub only allows certain file extensions. I saved the codes as .ipynb and changed the extension to .txt. Pleae replace .txt to .ipynb
to open it. I would really appreciate it if anyone could help me to resolve this problem.
Thanks so much for your time and work!
Regards!
Dapeng

ComputeCurlTau_PleaseRenametxtToipynbToOpen.txt

[dcherian]

Sorry for the late response; I was out of town.

@ALDepp and I are working on this. Currently I have a rough implementation

def x_divh(U, V, grid, boundary=None):

    dy = grid.get_metric(U, "Y")
    dzu = grid.get_metric(U, "Z")
    dx = grid.get_metric(V, "X")
    dzv = grid.get_metric(V, "Z")
    
    UT = U * dy
    VT = V * dx
    
    if grid.arakawa == "C":
        # c grid needs the dz multiply
        # but we want to avoid this for b grid so that 
        # we don't create unnecessary tasks
        UT = UT * dzu
        VT = V * dzv
        
    elif grid.arakawa == "B":
        # B grid needs an interp step here
        #    U in Y; V in X
        UT = grid.interp(UT, "Y", boundary="extend")
        VT = grid.interp(VT, "X", boundary="extend")
    
    div = grid.diff(UT, "X", boundary="extend") + grid.diff(VT, "Y", boundary="extend")
    area = grid.get_metric(div, "XY")
    div = div / area
    
    # TODO: for C grid need (div / DZ) to get dimensionality right: ∂w/∂z
    #       but the xgcm example notebook does not do that.
        
    return div


def x_curlz(U, V, grid):
    
    Udx = U * grid.get_metric(U, "X")
    Vdy = V * grid.get_metric(V, "Y")
    
    if grid.arakawa == "B":
        # B grid needs an interp step here
        #    Udx in X; Vdy in Y
        Udx = grid.interp(Udx, "X", boundary="extend")
        Vdy = grid.interp(Vdy, "Y", boundary="extend")

    der = -grid.diff(Udx, "Y", boundary="extend") + grid.diff(Vdy, "X", boundary="extend")
    area = grid.get_metric(der, "XY")
    curlz = der / area
    
    return curlz

This leverages the new metrics functionality in xgcm: https://xgcm.readthedocs.io/en/latest/grid_metrics.html. The metrics setup could be added to to_xgcm_grid_dataset and then things would be really easy.

I am envisioning the following aspiration syntax that would work regardless of model:

vorticity = grid.curlz(ds.UVEL, ds.VVEL)
diveregence = grid.divh(ds.UVEL, ds.VVEL)

We will iterate on this and start a PR/discussion over at xgcm on adding these higher order operators to xgcm soon.

[DapengLi17]

Hi @dcherian
Thanks so much for your reply. I managed to compute curl on the POP original grids using xgcm following https://github.com/pangeo-data/pangeo-ocean-examples/blob/master/cesm-pop-highres-ocean.ipynb Now I have difficulty plotting the 0.1-deg CESM POP original grids output. I tried different Python cartopy projections and the pictures look strange. Please see the attached platecarre projection plot. Can you please let me know how to plot the CESM POP original grids output?
Thanks so much for your help!

[matt-long]

@DapengLi17 here's a gist that demonstrates one way to plot the POP grid using a pop_add_cyclic function:
https://gist.github.com/matt-long/50433da346da8ac17cde926eec90a87c

[DapengLi17]

Hi @matt-long
Thanks so much for your reply. I followed your gist and replaced "POP_gx1v7" in your codes with ''POP_tx0.1v3'' for my 0.1-deg resolution grids. It seems the pic is still not correct. The Asia and America continents seem to be out of shaped. Please see the attached plot and my jupyter notebook. https://github.com/DapengLi17/JupyterNotebook/blob/master/pltPOPgrids_2020Mar21.ipynb I am not sure if it is ok to plot without specifying the cartopy projections. Any comments are highly appreciated.
Regards!
Dapeng

[matt-long]

pop_add_cyclic just rearranges the grid for plotting with a projection. Looks to me like you are just visualizing the array itself with no projection.

[DapengLi17]

Hi @matt-long
Just figured out the projection issue. Now it works fine. Thanks so much for your kind help.
Out of curiosity, your util.pop_add_cyclic works on the grid file (e.g. POP_tx0.1v3) not the CESM POP output large nc data files. Can we just rearrange lat and lon in the CESM POP original output nc files? Then we do not need the grid file anymore since lat and lon are available in the CESM POP output nc files.

[matt-long]

You can apply pop_add_cyclic to any dataset that has TLAT and TLONG.

[DapengLi17]

Hi @matt-long
Thanks so much for your reply. I guess the reason why my laptop crashed when I applied pop_add_cyclic to the original CESM POP nc files is because of the nc file size and my laptop memory.
Thanks again for your kind help!
Regards!
Dapeng

[bradyrx]

@andersy005, @dcherian I think this can be closed in favor of xgcm/xgcm#187.