Analytical grid generation with pyroms

[bperfect]

Hello,

I am using ROMS to do some modeling with an idealized domain. At the behest of a mentor, I have installed pyroms (which appears to be working fine) in order to overcome some of the limitations of the ana_* files. My question is, how would I use pyroms to generate a full set of netCDF files (ocean_grd, ocean_clm, ocean_ini, ocean_bry) from scratch? Does anyone have any code that deals with this? Everything I can find in the examples folder either starts from formatted data files, or pre-existing netCDF files.

[fagundesmo]

I am not sure I'm the best person to give advices since I'm only an undergraduate student but what I did was I took the function to interpolate from z to sigma levels from pyroms and the horizontal interpolation I did using another function and it worked. The ocean_grd I had some trouble with the mask especially because my domain is the African and Brazilian continents and since pyroms uses the coastsegs generated by python it didn't mask all the continents correctly,the thing that I and other colleague did was:

After, I set the resolution of my ocean_grd I used to use the following code:

Code: Select all

# Grid dimension
Mm=len(lats1) #latitude
Lm=len(lons1) #Longitude

lon0=lons1.min() ; lat0 =lats1.max()
lon1=lons1.min() ; lat1 =lats1.min()
lon2 =lons1.max() ; lat2 =lats1.min()
lon3 =lons1.max() ; lat3 =lats1.max()

map = Basemap(projection='merc', llcrnrlon=lons1.min(), llcrnrlat=lats1.min(),urcrnrlon=lons1.max(), urcrnrlat=lats1.max(), resolution='f')

lonp = np.array([lon0, lon1, lon2, lon3])
latp = np.array([lat0, lat1, lat2, lat3])
beta = np.array([1., 1., 1., 1.])
hgrd = pyroms.grid.Gridgen(lonp, latp, beta, (Mm,Lm),proj=map)

lonv, latv = map(hgrd.x_vert, hgrd.y_vert, inverse=True)

hgrd = pyroms.grid.CGrid_geo(lonv, latv, map)
for verts in map.coastsegs:
    hgrd.mask_polygon(verts)

coast = pyroms.utility.get_coast_from_map(map)
pyroms.grid.edit_mask_mesh_ij(hgrd, coast=coast)

But then I and my other colleague realized that we could use the function np.where from python to do the job and we modified from the previous to:

Code: Select all

# Grid dimension
Mm=len(lats1) #latitude
Lm=len(lons1) #Longitude

lon0=lons1.min() ; lat0 =lats1.max()
lon1=lons1.min() ; lat1 =lats1.min()
lon2 =lons1.max() ; lat2 =lats1.min()
lon3 =lons1.max() ; lat3 =lats1.max()

map = Basemap(projection='merc', llcrnrlon=lons1.min(), llcrnrlat=lats1.min(),urcrnrlon=lons1.max(), urcrnrlat=lats1.max(), resolution='f')

lonp = np.array([lon0, lon1, lon2, lon3])
latp = np.array([lat0, lat1, lat2, lat3])
beta = np.array([1., 1., 1., 1.])

hgrd = pyroms.grid.Gridgen(lonp, latp, beta, (Mm,Lm),proj=map)

lonv, latv = map(hgrd.x_vert, hgrd.y_vert, inverse=True)

hgrd = pyroms.grid.CGrid_geo(lonv, latv, map)
    
mask_rho=np.where(bath1[:-1,:-1]>-0.1,0.,hgrd.mask_rho)
mask_rho=np.where(bath1[:-1,:-1]<-0.1,1.,mask_rho)
hgrd.mask_rho=mask_rho
coast = pyroms.utility.get_coast_from_map(map)
pyroms.grid.edit_mask_mesh_ij(hgrd, coast=coast)

where:
bath1 is the bathymetry already interpolated using the resolution of my domain, bath1 is bigger than the mask_rho so I adjusted that.

It saved me some time otherwise I'd have either make some loops to mask the land or click with the pyroms.grid.edit_mask_mesh_ij command.

I used SODA to the ocean and NCEP/NCAR to the atmosphere. For the atmosphere you could think as it was the SSH data so you could either interpolate to your domain following .cdl files or you could put all the data in a .nc file, but again I am not the best person to answer questions related to pyroms since I've used this toolbox only for my final project since January and if I said anything wrong about the commands and etc I'd appreciate the help too.

I hope it can help you at least a little.

Matheus