<jupyter_start><jupyter_text>Code to check c-grid definition NEMO vs MITgcmMake three fieldsets (radial rotation example):- NEMO: f-points in northeast, u-points in the east, v-points in the north- MIT: f-points in the southwest, u-points in the west, v-points in the south- corrected: like MIT, but u shifted in j-direction, v shifted in i-direction<jupyter_code>#make fieldsets based on NEMO c-grid, MITgcm c-grid and an MITgcm c-grid that is corrected

xdim=ydim=201
a = b = 20000  # domain size
lon = np.linspace(-a/2, a/2, xdim, dtype=np.float32)
lat = np.linspace(-b/2, b/2, ydim, dtype=np.float32)
dx, dy = lon[2]-lon[1], lat[2]-lat[1]
omega = 2 * np.pi / delta(days=1).total_seconds()

def calc_r_phi(ln, lt):
    return np.sqrt(ln**2 + lt**2), np.arctan2(ln, lt)

def calculate_UVR(lat,lon,dx,dy,omega,cgrid='nemo'):
    U = np.zeros((lat.size, lon.size), dtype=np.float32)
    V = np.zeros((lat.size, lon.size), dtype=np.float32)
    R = np.zeros((lat.size, lon.size), dtype=np.float32)
    
    for i in range(lon.size):
        for j in range(lat.size):
            r, phi = calc_r_phi(lon[i], lat[j])
            R[j, i] = r
            if cgrid=='nemo':
                r, phi = calc_r_phi(lon[i]-dx/2, lat[j])
                V[j, i] = -omega * r * np.sin(phi)
                r, phi = calc_r_phi(lon[i], lat[j]-dy/2)
                U[j, i] = omega * r * np.cos(phi)
            if cgrid=='mit':
                r, phi = calc_r_phi(lon[i]+dx/2, lat[j])
                V[j, i] = -omega * r * np.sin(phi)
                r, phi = calc_r_phi(lon[i], lat[j]+dy/2)
                U[j, i] = omega * r * np.cos(phi)
    
    return U, V, R

#save NEMO fieldset 
U,V,R = calculate_UVR(lat,lon,dx,dy,omega,cgrid='nemo')
data = {'U': U, 'V': V, 'R': R}
dimensions = {'lon': lon, 'lat': lat}
fieldset_NEMOstyle = FieldSet.from_data(data, dimensions, mesh='flat')
fieldset_NEMOstyle.U.interp_method = 'cgrid_velocity'
fieldset_NEMOstyle.V.interp_method = 'cgrid_velocity'
fieldset_NEMOstyle.U.gridindexingtype = 'nemo'
fieldset_NEMOstyle.V.gridindexingtype = 'nemo'

#save MITgcm fieldset 
U,V,R = calculate_UVR(lat,lon,dx,dy,omega,cgrid='mit')
data = {'U': U, 'V': V, 'R': R}
dimensions = {'lon': lon, 'lat': lat}
fieldset_MITstyle = FieldSet.from_data(data, dimensions, mesh='flat')
fieldset_MITstyle.U.interp_method = 'cgrid_velocity'
fieldset_MITstyle.V.interp_method = 'cgrid_velocity'
fieldset_MITstyle.U.gridindexingtype = 'mitgcm'
fieldset_MITstyle.V.gridindexingtype = 'mitgcm'

#corrected MITgcm fieldset
UU = np.zeros((lat.size, lon.size), dtype=np.float32)
VV = np.zeros((lat.size, lon.size), dtype=np.float32)
UU[1:,:]=U[0:-1,:] #shift grids of U and V by adding row/column of zeros 
VV[:,1:]=V[:,0:-1] #(U in j-direction, V in i-direction)
data = {'U': UU, 'V': VV, 'R': R}
dimensions = {'lon': lon, 'lat': lat}
fieldset = FieldSet.from_data(data, dimensions, mesh='flat')
fieldset_corrected = FieldSet.from_data(data, dimensions, mesh='flat')
fieldset_corrected.U.interp_method = 'cgrid_velocity'
fieldset_corrected.V.interp_method = 'cgrid_velocity'
fieldset_corrected.U.gridindexingtype = 'nemo'
fieldset_corrected.V.gridindexingtype = 'nemo'<jupyter_output><empty_output><jupyter_text>Run same particle in three different fieldsSecond figure shows that without changing UV grid, c-grid definition of the MITgcm gives different results.<jupyter_code>def UpdateR(particle, fieldset, time):
    particle.radius = fieldset.R[time, particle.depth, particle.lat, particle.lon]

class MyParticle(ScipyParticle):
    radius = Variable('radius', dtype=np.float32, initial=0.)
    radius_start = Variable('radius_start', dtype=np.float32, initial=fieldset_NEMOstyle.R)

pset = ParticleSet(fieldset_NEMOstyle, pclass=MyParticle, lon=0, lat=4e3, time=0)

output = pset.ParticleFile(name='Cgrid_NEMOstyle.nc', outputdt=delta(hours=1))
pset.execute(pset.Kernel(UpdateR) + AdvectionRK4,
             runtime=delta(hours=12),
             dt=delta(minutes=5),  
             output_file=output)
output.close()

class MyParticle(ScipyParticle):
    radius = Variable('radius', dtype=np.float32, initial=0.)
    radius_start = Variable('radius_start', dtype=np.float32, initial=fieldset_MITstyle.R)

pset = ParticleSet(fieldset_MITstyle, pclass=MyParticle, lon=0, lat=4e3, time=0)

output = pset.ParticleFile(name='Cgrid_MITstyle.nc', outputdt=delta(hours=1))
pset.execute(pset.Kernel(UpdateR) + AdvectionRK4,
             runtime=delta(hours=12),
             dt=delta(minutes=5),  
             output_file=output)
output.close()

class MyParticle(ScipyParticle):
    radius = Variable('radius', dtype=np.float32, initial=0.)
    radius_start = Variable('radius_start', dtype=np.float32, initial=fieldset_corrected.R)

pset = ParticleSet(fieldset_corrected, pclass=MyParticle, lon=0, lat=4e3, time=0)

output = pset.ParticleFile(name='Cgrid_corrected.nc', outputdt=delta(hours=1))
pset.execute(pset.Kernel(UpdateR) + AdvectionRK4,
             runtime=delta(hours=12),
             dt=delta(minutes=5),  
             output_file=output)
output.close()
plotTrajectoriesFile('Cgrid_NEMOstyle.nc')
plotTrajectoriesFile('Cgrid_MITstyle.nc')
plotTrajectoriesFile('Cgrid_corrected.nc')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.149099271751629, 'fraction_numerical': 0.016673054810272134, 'mean_word_length': 4.2770475227502525, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '59921', 'license': 'MIT', 'n_tokens_mistral': 2037, 'n_tokens_neox': 1957, 'n_words': 392, 'path': 'documentation/20.08_CgridTest/Test_Cgrid.ipynb', 'repo_name': 'SLYpma/GalapagosAnalysis', 'size': '5218'}