Central East Florida Shelf Realistic Application Setup Help

[cflanary]

I'm working on my dissertation to model the cold water upwelling that occurs most every summer along the central east Florida shelf. I have a long ADCP record in the nearshore and a few summers of bottom mounted ADCP data across the shelf to the shelf break, which I will use as a diagnostic comparison to model results.

I've made myself familiar with the ROMS code and the test cases and have read all of the WikiROMS pages and frequent this forum. I'm building my realistic application and I know the basic steps that need to be done, modify the cpp options, make my own *.in file, modify the stations.in file, and construct my initial, boundary conditions and forcing files. I've constructed my grid via SeaGrid.

What I don't know is the specifics of realistic application setup. Is there documentation available with examples of taking larger model data, in my case RTOFS for the initial conditions, NAMS and the AOML Gulf Stream cable for forcing, and t_tide for tidal conditions, and modifying it to the appropriate input files???

I know MATLAB very well and I've seen Kate mention that most people use MATLAB or Python scripts rather than the older CDL and Fortran options for constructing these files, but I'm finding it difficult to sort through all the available files as for what does what.

I will start with that, so thank you for any information you can provide.

[johnluick]

I'm in a hotel room somewhere and unable to check my facts, but you have reminded me of a question so will take the opportunity to reply. As I recalled there was one (and I think only one) realistic coastal application a few years ago when I first started using Rutgers ROMS, called NJ_BIGHT or something. Not long ago I went back to check something about it and it was gone. I don't remember seeing anything in the "In Memoriams" announcing its demise. I've always wondered what happened to old Whatitsname.

[arango]

There are couple of realistic cases in the test repository. You can download the full suite of ROMS idealized and realistic test cases. These will give you a good idea how a realistic case is configured. Also in the User sub-directory of the ROMS repository we have several configurations for various applications. There you can find the NJ_BIGHT configuration. However, we don't provide the NetCDF input files for the applications in the User sub-directory.

There is a lot of software in the ROMS repository to do all what you need for configuring and analyzing an application.

[wilkin]

If you look in the left had navigation panel at http://www.myroms.org for "Applications" you'll see a set of examples there. Admittedly, few of these are really up to date, but http://www.myroms.org/applications/espresso has some descriptions of how we run our US East Coast system and this is a similar set-up to your proposed project. The opendap links to output allow you to browse the full ROMS output, and since the metadata in ROMS is very comprehensive you can see the CPP options used, the various parameter choices, and back out the grid resolution etc.

[cflanary]

Thank you for the information, the espresso was a nice example to see what they have used. My biggest problem is that there is a lot of software available and it's difficult to know what to use. I'm unfortunately the only ROMS user at my school. I will start with a few direct questions and that should get me moving forward.

Is there more documentation on the input files in the *.in file other than in the glossary, specifically the multiple initial conditions files, forward solutions field, and adjoint sensitivity functional file as seen in the ocean_upwelling.in example file?

How is the ocean upwelling test example able to run without having all the input, boundary and forcing files commented out, since they don't exist, ie. in my realistic model run, is there a cpp flag that will have the model recognize these files?

Is the OA package still used with the newest available matlab codes for creating ROMS input files?

Thank you again for your information.

[hetland]

All ROMS applications can take forcing and bathymetry (and other) information from either an analytical file (look in the Analytical directory for many, many examples) or from a netCDF file. The variables, in the end, are the same to ROMS, it is just a matter of whether they are defined by FORTRAN code inside of ROMS, or a netCDF file created by some script outside of ROMS.

So, the first step is to familiarize yourself with the variables you will need. You will usually place these variables in a file called foo_grd.nc where foo is your application name. The key variables are actually only the grid metrics. You need bathymetry, h, Coriolis parameter, f, pm and pn, the inverse grid widths, and dmde and dndx, which are curvature parameters. You can find definitions for creating these parameters in the various scripts. I use a package called octant that I developed in Python. However, it seems everybody has their own set of tools. Usually people are willing to share. You can find mine at http://octant.googlecode.com. It may not be much use to you practically, as a matlab user, but you might find some algorithms that you can translate.

In terms of actually creating a realistic grid, this is a very complicated process. There is no code that just does this (although some tools in ROMS-AGRIF might come close...). As an example, I have attached a script I use to make a grid in octant for the Merrimack River region. The script was developed on the command line, and saved for later use, and shows the general workflow. Again, this might not be useful to you practically, but perhaps it will give you an idea of the steps involved. Much of the heavy lifting is buried in octant -- actually generating the curvilinear grid (I use a program called gridgen for this), calculating the Coriolis parameter, grid metrics, creating the land mask from coastline polygons.

Good luck..

Code: Select all

import os
from glob import glob
import cPickle

import numpy as np
import pylab as plt
from mpl_toolkits.basemap import Basemap
from scikits import delaunay
from scipy.ndimage import gaussian_filter

import octant
import octant.roms

import map_tools
from map_tools import proj

shp = (515, 515)

# load in boundary data from file
bry = np.load('boundary5.pickle')
lonbry, latbry = proj(bry['x'], bry['y'], inverse=True)
beta = bry['beta']
ul_idx = 1  # index of upper left hand corner

foc = octant.grid.Focus()
foc.add_focus_x(0.40, 3, 0.7)
foc.add_focus_y(0.6, 4, 0.6)

## for building and testing the grid, some boundary_interactor calls:
# bry.gridgen_options['focus'] = foc   # lines for working with BoundaryInteractor
# bry.gridgen_options['ul_idx'] = ul_idx
# bry.gridgen_options['shp'] = (131, 131)

xbry, ybry = proj(lonbry, latbry)
# running pyroms.gridgen returns a pyroms.Grid object

print '...running gridgen.'
if os.path.exists('gridgen.pickle'):
    print "   ### loading gridgen results from 'gridgen.pickle'"
    grd = np.load('gridgen.pickle')
else:
    grd = octant.grid.Gridgen(xbry, ybry, beta, shp, \
                    focus=foc, ul_idx=ul_idx, verbose=True)
    grd.focus = None  # can't pickle functions...
    f = open('gridgen.pickle', 'wb')
    cPickle.dump(grd, f, protocol=-1)
    f.close()

# custom set coreolis paremeter (otherwise determined from grd.lat)
grd.f = 1.0e-4 * np.ones(grd.x_rho.shape, dtype='d')

# load in the bathymetry
lonb, latb, bathy = np.load('new_bathy.npy').T
xb, yb = proj(lonb, latb)

# interpolate the bathymetry to the grid
idx = np.where((xb<grd.x_rho.max()) & (xb>grd.x_rho.min()) & \
               (yb<grd.y_rho.max()) & (yb>grd.y_rho.min()))
xb = xb[idx]
yb = yb[idx]
bathy = bathy[idx]
xb, idx = np.unique1d(xb, return_index=True)
yb = yb[idx]
bathy = bathy[idx]
h = np.zeros_like(grd.x_rho)
h[grd.mask==1] = delaunay.Triangulation(xb, yb)\
                    .nn_extrapolator(bathy)(grd.x_rho[grd.mask==1], 
                                            grd.y_rho[grd.mask==1])


h_min = 2.0
h_max = 50.0
h = h.clip(min=h_min, max=h_max)

print '...smooth the bathymetry and straighten out the edges'
hf_n = gaussian_filter(h, 1.0)       # modestly smoothed bathy in the north (old val 0.75)
hf_s = gaussian_filter(h, 2.0)        # strongly smoothed bathy in the south

ii, jj = np.mgrid[0:514, 0:514]
# weighting for blending bathy.  One in the north, zero in the south.
weight = ((ii - 200.0) / 50.0).clip(min=0.0, max=1.0)
hf = weight * hf_n + (1-weight) * hf_s

print '...modifying min/max depth and depth array at boundaries'



print '...masking land from maskfiles'
for maskfile in glob('*.mask'):
    lon, lat = plt.load(maskfile, unpack=True)
    x, y = proj(lon, lat)
    grd.mask_polygon(zip(x, y))


print '...fix mask near south jetty'
lon = np.array([-70.809356641927451,
                -70.807724158705241,
                -70.80674047240781 ,
                -70.807340045005887,
                -70.809260613502616])
lat = np.array([42.811496540734552,
                42.811989712803737,
                42.816295717852896,
                42.816240386352455,
                42.815277526598322])
x, y = proj(lon, lat)
grd.mask_polygon(zip(x, y))


h_north = hf[-88, :]
h_north[3:-3]
for j in range(1, 88):
    hf[-j, :] = hf[-88,:]
    grd.mask_rho[-j, :] = grd.mask_rho[-88,:]

for i in range(2):
    hf[:, i] = hf[:,2]
    grd.mask_rho[:, i] = grd.mask_rho[:,2]

# for j in range(10):
#     hf[j, :] = hf[10,:]
#     grd.mask_rho[j, :] = grd.mask_rho[10, :]

grd.h = hf.copy()

grd.mask[13, 288] = 1
grd.mask[14, 281] = 1
grd.mask[15, 282] = 1
grd.mask[17, 274] = 0
grd.mask[17, 275] = 0
grd.mask[17, 276] = 0

print '...pickling grid'
grd.focus = None  # can't pickle functions...
f = open('grd.pickle', 'wb')
cPickle.dump(grd, f, protocol=-1)
f.close()

grd.dx = grd.dx.filled( grd.dx.mean() )
grd.dy = grd.dy.filled( grd.dy.mean() )
grd.dndx = grd.dndx.filled( grd.dndx.mean() )
grd.dmde = grd.dmde.filled( grd.dmde.mean() )
grd.angle = grd.angle.filled( grd.angle.mean() )

print '...writing netcdf file'
octant.roms.write_grd(grd, 'midfield_3_grd.nc', verbose=True)

[ruoying]

cflanary, you might want to take a look at our ROMS setup that covers both the Gulf of Mexico and South Atlantic Bight (encompassing the east Florida shelf). This model is running everyday providing daily regional ocean nowcast and 72-hour forecast.

http://omgsrv1.meas.ncsu.edu:8080/ocean-circulation/

Regarding upwelling on the east Florida shelf you mentioned, we have a paper on this topic, focusing on summer 2003:

Hyun, K. H. and He, R. (2010), Coastal Upwelling in the South Atlantic Bight: A Revisit of the 2003 Cold Event Using Long Term Observations and Model Hindcast Solutions, Journal of Marine Systems, v83, 1-13.

The paper has more info on how the model was set up.

Ruoying

[maxblu33]

Dear hetland
I have struggle to make a grid with pyroms(from hedstrom), when I try to interpolate the bathimetry using griddata from matplotlib.mlab it gives me nothing and the the process is killed after 15 minutes of a blank screen, I saw how you do your bathimetry interpolation in the previous scrip example and try to do it with matplotlib dealunay but it did not work, with an error regarding
" popitem():dictionary is empty"

can you give any suggestions about scikits.dealuney, because I'm using Anaconda python distribution and there are many different version of scikits for biology, learning etc
Another question:
I try to install ocatant but when I run python setup.py install there is an error about missing numpy but I'm already using it with pyroms
thanks in advance
Luis Sorinas-Morales