[ecco-support] MITgcm on tiles: vorticity and interpolation coded

[Fabien Roquet]

Hello,
I have coded a function to compute the vorticity on the G-point and a function to interpolate a tracer from G-point to C-point. Are you interested? Also I have a function to turn the orientation of tiles to plot velocities. If you are interested, can you check it works well and maybe add it to the ecco module?
Cheers,
Fabien



def turn_tiles(oceX,oceY):
    oceE=oceX.copy()
    oceE[6:13,:,:]=oceY[6:13,:,:]
    oceE[6,:,:]=np.nan
    oceN=oceY.copy()
    oceN[6:13,:,:]=-oceX[6:13,:,:]
    oceN[6,:,:]=np.nan
    return oceE,oceN

def compute_curl(fieldU,fieldV,ds_grid):
    # compute the curl of a velocity field on the MITgcm grid. Result on the G-point.

    fieldU *= ds_grid.dxC
    fieldV *= ds_grid.dyC

    dUdY = xr.zeros_like(ds_grid.XG)
    dUdY.values[:,1:,:] = fieldU.values[:,1:,:] - fieldU.values[:,:-1,:]
    dUdY.values[[1,2,4,5,6],0,:] = fieldU.values[[1,2,4,5,6],0,:] - fieldU.values[[0,1,3,4,5],-1,:]
    dUdY.values[10:,0,:] = fieldU.values[10:,0,:] - fieldU.values[7:10,-1,:]
    dUdY.values[7:10,0,1:] = fieldU.values[7:10,0,1:] + fieldV.values[[5,4,3],1:,-1][:,::-1]
    dUdY.values[7:10,0,0] = fieldU.values[7:10,0,0] + fieldV.values[[6,5,4],0,-1]

    dVdX = xr.zeros_like(ds_grid.XG)
    dVdX.values[:,:,1:] = fieldV.values[:,:,1:] - fieldV.values[:,:,:-1]
    dVdX.values[[3,4,5,7,8,9,11,12],:,0] = fieldV.values[[3,4,5,7,8,9,11,12],:,0] - fieldV.values[[0,1,2,6,7,8,10,11],:,-1]
    dVdX.values[[0,1,2,6,10],1:,0] = fieldV.values[[0,1,2,6,10],1:,0] + fieldU.values[[12,11,10,2,6],-1,1:][:,::-1]
    dVdX.values[[1,2,6,10],0,0] = fieldV.values[[1,2,6,10],0,0] + fieldU.values[[12,11,5,7],-1,0]

    zeta = (-dUdY + dVdX)/ecco_ds.rAz

    return zeta

def interp_G_to_C(fieldG,ds_grid):
    # interpolate a tracer field from G-point to C-point

    fieldU = grid.interp(fieldG,['Y'],boundary='fill')

    fieldC = xr.zeros_like(ds_grid.XC)
    fieldC.values[:,:,:-1] = (fieldU.values[:,:,:-1] + fieldU.values[:,:,1:]) / 2
    fieldC.values[[0,1,2,6,7,8,10,11],:,-1] = (fieldU.values[[0,1,2,6,7,8,10,11],:,-1] + fieldU.values[[3,4,5,7,8,9,11,12],:,0]) / 2
    fieldC.values[[5,4,3],1:,-1] = (fieldU.values[[5,4,3],1:,-1] + fieldU.values[7:10,0,1:][:,::-1] ) / 2
    fieldC.values[[5,4],0,-1] = (fieldU.values[[5,4],0,-1] + fieldU.values[[8,9],0,0] ) / 2

    return fieldC




-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.mit.edu/pipermail/ecco-support/attachments/20211017/2c248199/attachment.html