Skip to: [content] [navigation]

Description


1) Converting data to NetCDF

NetCDF files must follow the formatting conventions for the standardization of NetCDF files established by the Cooperative Ocean/Atmosphere Research Data Service (COARDS). Some examples are provided at the ferret web site. More specifically, MOM4 is finicky about time units, which must adhere to the following additional requirements:

• time axis name must be either "time" or "t"
• time axis units must be "sec","min","hour(s)", or "day(s)", all case insensitive.

The easiest way I have found to set up NetCDF files is using the combination of ferret and NCO utilites. Matlab also has a bunch of utilities.

As an example of a data set whose conversion is nontrivial, regridding the OCMIP2-biotic standard euphotic zone boundary condition for phosphate (PO4) proceeds as follows:

a) Get the 4-byte binary data (po4mapnew.dat.gz) from the website.

b) uncompress it with:
gunzip po4mapnew.dat

c) Read the description file (README.po4maps):


po4mapnew.dat.gz : PO4 3D fields (12 months, 0 to 75 m). The file is compressed by command gzip. To uncompress it, use "gunzip". readmap.f : fortran program to read the map. The file is on a regular 2deg per 2deg grid. First latitude is 90N , first longitude is Greenwitch, The standard levels of Levitus from 0 to 75 m are : 0, 10, 20, 30, 50, 75 File nutrient-doc contains all information about the creation of monthly maps.

d) Since this file is of the type access='sequential', we first need to change it to access='direct' to be correctly read into ferret. Otherwise, the data will be offset by the extra formatting. This is achieved using the supplied fortran code, after adding two lines:


after the open po4mapnew.dat line:


open(22,file ='po4mapnew_direct.dat',access='direct',form='unformatted', & recl=90*180)


and after the read po4mapnew.dat line:

write(22) a

e) In ferret, define a grid, read the data into that grid, and save out a temporary file with reversed latitude axis (NOTE: be sure to use the correct big-endian/little-endian format - e.g if the data look like garbage, try adding /SWAP to the options list in the FILE read command), save out a temporary file with reversed latitude axis and list out y axis to reverse using ncdump and ncgen:


DEFINE AXIS/X=1E:359E:2 xaxis
DEFINE AXIS/Y=89S:89N:2 yaxis
DEFINE AXIS/Z/UNITS=meters/DEPTH zaxis={0, 10, 20, 30, 50, 75}
DEFINE AXIS/T=15-jan-1991:15-dec-1991:30.4/UNITS=day timeaxis
DEFINE GRID/X=xaxis/Y=yaxis/Z=zaxis/T=timeaxis xyztgrid
FILE/VAR=PO4_ocmip/GRID=xyztgrid/FORMAT=stream/SWAP po4mapnew_direct.dat
SAVE/FILE="ocmip_po4.nc"/CLOBBER po4_ocmip
LIST/y=90s:90n/format=(10(F6.1,","))/order=y -1*y[gy=yaxis]
QUIT

f) Write out netCDF file
ncdump ocmip_po4.nc > ocmip_po4.cdl

g) Open ocmip_po4.cdl to edit it by replacing the yaxis information with the reversed listing printed out in ferret.

h) Create reversed NetCDF file from description:
ncgen -o ocmip_po4.nc ocmip_po4.cdl

An alternative method, if you're having trouble getting the NetCDF to look right, is use matlab:

a) First, we create a .txt file from the .dat file using the supplied fortran code, after adding two lines:


open(22,file ='po4mapnew.txt',form='formatted')


near the beginning, and:


write(22,3) alon(jj),alat(ii),a(jj,ii)


in the loop.

b) Start matlab. I prefer the option without the Graphical User Interface (GUI), as it is much faster and more stable:
matlab -nojvm

c) Use matlab to over-write the variable PO4_OCMIP in ocmip_po4.nc (note that matlab will only write over a variable with this version of the NetCDF tools. It cannot create a variable):

addpath ~jpd/matlab6tools ~jpd/matlab6tools/toolbox ~jpd/matlab6tools/toolbox/netcdf ~jpd/matlab6tools/toolbox/netcdf/ncutility/ load po4mapnew.txt lat=reshape(po4mapnew(:,1),[180 90 6 12]);
lon=reshape(po4mapnew(:,2),[180 90 6 12]);
po4=reshape(po4mapnew(:,3),[180 90 6 12]);
lat_ocmip=zeros([12 6 90 180]);
lon_ocmip=zeros([12 6 90 180]);
PO4_OCMIP=zeros([12 6 90 180]);
for t=1:12
for k=1:6
for j=1:90
for i=1:180
lat_ocmip(t,k,91-j,i)=lat(i,j,k,t);
lon_ocmip(t,k,91-j,i)=lon(i,j,k,t);
PO4_OCMIP(t,k,91-j,i)=po4(i,j,k,t);
end
end
end
end
ncsave('ocmip_po4.nc',PO4_OCMIP)

2) Extrapolating the data vertically over the MOM4 grid domain

MOM4 requires that all depths in a model input have a value. For example, tracers in the OM2 version of MOM4 are defined from 5 meters to 5316.4 meters. For convenience - since I already have gone through the trouble of getting World Ocean Atlas data into NetCDF, I use the Levitus grid for this extrapolation, which goes from 0 meters to 5500 meters.


a) In order to run ferret on the Analysis Cluster (AC), you must have something equivalent to the following in your ~/.cshrc file:


if ( `hostname` =~ anc[1-2] ) then

# put aliases unique to AC graphics software here
set ferret_path_file = /home/mh2/.ferret_paths_i686
set mh2_ext_fcns = /home/mh2/.ferret_paths_ac
if (-e ) then

source
source

endif
setenv FMS_BIN_DIR /net2/mh2/fms/bin

endif

b) Log onto the AC:
qlogin -pe ac.inter 1 exec tcsh cd ~/jakarta/scripts mk_om1_ocmip2_biotic
c) Start ferret, read in the fields, and create a new variable which is OCMIP PO4 in the upper 75 meters, and WOA01 everywhere else:


USE "/home/jpd/input_fields/ocmip/ocmip_po4.nc"
USE "/home/jpd/input_fields/woa01_phosphate_annual.nc"
LET po4_full=if (z[g=phosphate[d=2]] lt 80) then po4_ocmip[d=1,gz=GBP1] else phosphate[d=2,gz=GBP1,l=1]


c-alternate) Alternatively, one could create an arbitrary vertical axis and then extrapolate automatically, for example:


DEFINE AXIS/Z/UNITS=meters/DEPTH zaxis2={0, 10, 20, 30, 50, 75, 100, 1000, 2000, 3000, 4000, 5000, 6000}
DEFINE GRID/X=xaxis/Y=yaxis/Z=zaxis2/T=timeaxis xyztgrid2
LET po4_full=phosphate[d=2,gz=zaxis2,z=@fnr,i=@fnr,j=@fnr]

3) Extrapolation of the data horizontally over the globe:

Continue the ferret session with Matt Harrison's extrapolation function:


let po4=extrap(po4_full,500,0.001,0.6)
SET MEMORY/SIZE=100
SAVE/FILE="ocmip_po4_monthly_om2_extrap.nc"/CLOBBER po4

4) Regridding extrapolated data onto OM1 grid

Grid Generation and re-gridding c-shell scripts have been developed by Zhi Liang and are documented in the MOM4 guide. This software is automatically included in the source code of the model, and can be found in:


/home/[USERNAME]/jakarta/om1_ocmip2_biotic/src/preprocessing/regrid_3d/regrid_3d.csh

At the top of regrid_3d.csh, you will need to edit the path/file information in the following lines to suit your particular needs:


set output_dir = /preprocessing/regrid_3d
set src_file = /archive/fms/mom4/input_data/levitus_ewg.nc
set dest_grid = /home/z1l/mom4_data_set/test3/grid_orig.nc
set dest_file = regrid_3d


In addition, lower down in regrid_3d.csh, you will need to change the variable name information:


numfields = 2
src_field_name = 'temp', 'salt'


In the AC, run the shell script simply by:
cd /home/[USERNAME]/jakarta/om1_ocmip2_biotic/src/preprocessing/regrid_3d/ regrid_3d.csh

5) making modifications to the final NetCDF header/description

If you need to further edit the netCDF file in order to conform to formatting requirements, you may want to edit the file. To do this, you must first dump the file to ascii format and then probably (if it's big) separate the file into parts after finding out how many lines there are:

ncdump woa01_sio4_om1_bc.nc > woa01_sio4_om1_bc.cdl wc -l woa01_sio4_om1_bc.cdl head --lines=12630 woa01_sio4_om1_bc.cdl > woa01_sio4_om1_bc.head tail --lines=2699284 woa01_sio4_om1_bc.cdl > woa01_sio4_om1_bc.tail

Make changes to head to conform to formatting... see example in /archive/jpd/postinchon/input_fields/woa01_po4_om1_bc.nc, which has the following header information (upon ncdump -h [FILENEME]):


netcdf woa01_po4_om1_bc {
dimensions:
grid_x_T = 96 ;
grid_y_T = 40 ;
zt = 24 ;
TIME = UNLIMITED ; // (12 currently)
variables:
float grid_x_T(grid_x_T) ;
grid_x_T:cartesian_axis = "X" ;
grid_x_T:units = "degree_east" ;
float grid_y_T(grid_y_T) ;
grid_y_T:cartesian_axis = "Y" ;
grid_y_T:units = "degree_north" ;
float zt(zt) ;
zt:long_name = "zt" ;
zt:units = "meters" ;
zt:cartesian_axis = "z" ;
zt:positive = "down" ;
double TIME(TIME) ;
TIME:units = "seconds since 0000-01-01 00:00:00" ;
TIME:time_origin = "01-JAN-0000 00:00:00" ;
TIME:modulo = "T" ;
double x_T(grid_y_T, grid_x_T) ;
x_T:long_name = "noname" ;
x_T:units = "nounits" ;
double y_T(grid_y_T, grid_x_T) ;
y_T:long_name = "noname" ;
y_T:units = "nounits" ;
double po4(TIME, zt, grid_y_T, grid_x_T) ;
po4:long_name = "World Ocean Atlas 2001 Silicate" ;
po4:units = "mol PO4 m-3" ;

for a monthly climatology in units of seconds, the time axis should have the attributes of:

double TIME(TIME) ;
TIME:units = "SEC since 0000-01-01 00:00:00" ;
TIME:time_origin = "01-JAN-0000 00:00:00" ;
TIME:modulo = "y" ;

and have time values of:

time = 1296000, 3974400, 6480000, 9158400, 11750400, 14428800, 17020800, 19699200, 22377600, 24969600, 27648000, 30240000;

Make the necessary edits, then recombine the head and tail and convert back to NetCDF binary with:

cat woa01_sio4_om1_bc.head.head woa01_sio4_om1_bc.tail > woa01_sio4_om1_bc.cdl ncgen -o woa01_sio4_om2_bc_173.nc woa01_sio4_om1_bc.cdl

The file is now ready to be used as a boundary condition!

• About us
  • Welcome
  • Visiting GFDL
  • Overview
  • Tour
  • Gallery
  • People
  • Employment
  • Help
  • Notices
• Research
  • Atmospheric Physics and Chemistry
  • Biospheric Processes
  • Climate Diagnostics
  • Climate Dynamics and Prediction
  • Oceans and Climate
  • Weather and Atmospheric Dynamics
• Products and Services
  • Flexible Modeling System
  • Our Data Portal
  • Our GForge
  • Technical Services
  • Visualization
• Reference
  • Our Online Bibliography
  • Our Meetings and Seminars
  • Our Library
• GFDL Only
  • External forms and links
  • External web mail
  • Help Desk
  • Internal Administrative Services
  • Internal Technical Services
  • NOAA's Hot Items
  • webTA Login

smaller bigger reset