Processing and using NCAR ocean current data from processed CESM4 ocean model run, from NetCDF file.

The important key to using these data is the recognition that whatever projection is used for the model, the latitude and longitude of the center of each grid cell are stored in the NetCDF file. The goal is to extract the relevant data for each month and then georeference to a geographic latlong location. The native grid is xy, but represents the entire world.

This example uses an NCAR NetCDF file for ancient ocean currents

Key variables
UVEL is eastward current vector
VVEL is northward current vector
ULAT is the latitude of velocity grid points (TLAT is latitude for temperature and atmospheric aerosols)
ULONG is the longitude of velocity grid points (TLONG is latitude for temperature and atmospheric aerosols)
NETCDF_DIMENSION_z_t is the midpoint depth below surface of a band of water in which the current vector is modeled
NETCDF_z_t_units=centimeters are the units in which z_t is measured
Note also that time units are also important if you want to do a time series. These include time units (usually days) and start point. Equally important for long series is the year used. Commonly, a 365 day year us used. This means that it loses a day every 4 years. So this has to be taken into account. In some cases a 360 day year or other year is used. This is in the NetCDF header metadata.

For UVEL and VVEL files, there are 60 bands of different depths below surface for each month. I used the shallowest one, centered at 500 cm below surface for surface currents.

See also <http://www.ncl.ucar.edu/Applications/pop2lat.shtml> for doing this in NCL

1. Extract needed UVEL and VVEL files: bands 1, 61, 121, 181, 241, 301, 361, 421, 481, 541, 601, 661 for the first 12 months of the file

[Must import these into an xy location because they use an arbitrary (and rather weird 3 pole "grid" = POP Grid)]

EXAMPLE: r.in.gdal input=NETCDF:"/Users/cmbarton/Desktop/paleoclimate_data/b40.lm850-1850.1deg.001.pop.h.UVEL.100001-100912.nc":UVEL output=CCSM4_UVEL_1000feb band=61

2. Create points file from each relevant UVEL and VVEL band

[First set region to match one of the global grid files to make sure that all the points are extacted properly]

EXAMPLE: r.to.vect input=CCSM4_UVEL_1000jan output=CCSM4_UVEL_1000jan_g7 type=point column=currvect

Note that I'm using GRASS 7 for this. In GRASS 6, you cannot specify a column to hold the current vector data. You will need to add a column for the current vector and update the column as described below

3. Add new east and north fields to each of the points files (or current vector, east, and north if in GRASS 6)

EXAMPLE: v.db.addcolumn map=CCSM4_UVEL_1000jan_g7 layer=1 columns='lon double precision,lat double precision'

For GRASS 6:
EXAMPLE: v.db.addcolumn map=CCSM4_UVEL_1000jan_g7 layer=1 columns='currvect double precision,lon double precision,lat double precision'

4 Update points file east and north fields with lon and lat values from ULAT and ULONG subsets grid files  (or current vector, east, and north if in GRASS 6)

EXAMPLE: v.what.rast map=CCSM4_UVEL_1000jan_g7 layer=1 raster=VVEL.ULONG column=lon
EXAMPLE: v.what.rast map=CCSM4_UVEL_1000jan_g7 layer=1 raster=VVEL.ULAT column=lat

For GRASS 6, also:
EXAMPLE: v.what.rast map=CCSM4_UVEL_1000jan_g7 layer=1 raster=CCSM4_UVEL_1000jan column=currvect

5. Export points file to ascii

EXAMPLE: v.out.ogr input=CCSM4_UVEL_1000jan_g7 type=point layer=1 dsn=CCSM4_UVEL_1000jan format=CSV

6. Reimport points file into latlon region using the lat and lon field values as y and x respectively

[First switch to a latlon location and set region as needed]

EXAMPLE: v.in.ascii input=/Users/cmbarton/CCSM4_paleocurrent_AD1000/CCSM4_UVEL_1000jan_g7.csv output=CCSM4_UVEL_1000jan_g7 fs=, skip=1 columns='cat int,windvect double precision,lon double precision,lat double precision' x=3 y=4 cat=1

7. Reinterpolate to ETOPO2 resolution for N Atlantic

EXAMPLE: v.surf.bspline input=CCSM4_UVEL_1000jan_g7 raster=CCSM4_UVEL_1000jan.interp sie=0.8603 sin=0.8603 method=bicubic layer=1 column=windvect

8. Convert to current direction and speed

EXAMPLE FOR CURRENT VELOCITY: r.mapcalc expression='CCSM4_1000jun.windspeed=sqrt(pow(CCSM4_UVEL_1000jan.interp,2) + pow(CCSM4_VVEL_1000jan.interp,2))'

EXAMPLE FOR CURRENT VELOCITY: r.mapcalc expression='CCSM4_1000jun.winddir=atan(CCSM4_VVEL_1000jan.interp, CCSM4_UVEL_1000jan.interp)'