netcdf
Lee Byerle
labyerle at MET.UTAH.EDU
Sat Apr 2 15:34:43 EST 2005
Rajesh,
One way to do this is with the lats4d.gs script in the grads script library
at http://grads.iges.org/grads/gadoc/library.html (I attached it here).
Here is an example (assuming you have grads-readable data, in this case using a
control (ctl) file.
lats4d -v -ftype ctl -i u700.ctl -o newnc
where:
u700.ctl is the current control file (input)
newnc is the name of the new .nc (netcdf file) that will be created.
Good luck,
Lee
On Sat, 2 Apr 2005, Rajesh J wrote:
> hello,
> is it possible to write the output of grads to netcdf?
> thanks,
> rajesh
>
-------------- next part --------------
* file: lats4d.gs
*
* A minimum fuss gs script for writing NetCDF, HDF-SDS or GRIB
* files from GrADS using the PCMDI (http://www-pcmdi.llnl.gov) LATS
* interface. This script requires GrADS Version 1.7.beta.9 or later.
*
* See usage() for documentation; from GrADS enter "lats4d -h"
*
* REVISION HISTORY:
*
* 21dec1998 da Silva First crack.
* 24dec1998 da Silva Further development: options -lat -lon -cal -mean,
* better error trapping.
* 31dec1998 da Silva First beta version; automatic generation of
* LATS parameter table; now variables with 1 vertical
* level are identified as sfc variables.
* 04jan1999 da Silva Added -ftype option, automatic detection of
* file types, removed -sdf/-xdf options.
* 05jan1999 da Silva Cosmetic changes in man page.
* 06jan1999 da Silva Added -precision option.
* 08jan1999 da Silva Added -precision as valid option in
* parsecmd(), small other changes
* 11jan1999 da Silva now file name can have 3 tokens, so SDFopen
* can use templates.
* 17jan1999 da Silva Added -func and -freq options, fixed bug when
* there was only 1 time step on file; -func is
* currently undocumented.
* 22Jan1999 da Silva Fixed bug with -time option: now 0Z15dec1991
* works.
* 09Feb1999 da Silva Added -grid option; added detection of invalid
* time frequency in minutes (1.0.Beta.8)
* 17Mar1999 da Silva Fixed parsing bug in -freq; fixed bug in the
* evaluation of -func arguments; increased
* long name column size in parm table to 70
* to accomodate files going to ECS; removed FUNCTION
* from variable long name - this information is
* now part of the global attribute "comments"
* (1.0.Beta.9)
* 18Mar1999 da Silva Fixed another -func related bug in getgid().
* 19Mar1999 da Silva Added -func to options to avoid parsing error;
* revised documentation + small mod in -freq.
* 29Mar1999 da Silva Added -zrev option
* 28Apr1999 da Silva Added "-de" option (1.0.beta.13); now we
* can regrid; also added undocumented options
* -geos1x1, -geos1x1a, -geos1x1s, and same for
* 2x2.5, 4x5; these are DAO specific.
* 04May1999 da Silva Fixed small bug in -levs; now -levs 0.4 is
* possible.
* 30Jun1999 da Silva Added CPTEC to list of centers, fixed -vars
* to allow variables starting with 'e'.
* 13Jul1999 da Silva Fixed "-i" bug: now file names can start with "e".
* 16Jul1999 da Silva Changed order of upper-air var/level. This
* way the data is written in the same order they
* come in a standard IEEE grads file. Implemented
* stream (fwrite) output; no automatic ctl generation
* yet.
* 23Jul1999 da Silva Added -title option.
* 22Nov1999 da Silva Made v1.0.beta.18 the first stable version: v1.0.
* 23Feb2000 da Silva Special handling for grib output and p<1; now it
* uses hybrid level number for vertical coordinate
* in such cases.
* 01Jun2000 da Silva Added -xvars option.
* 16Jun2000 da Silva Added -xsfc/xupper options.
* 10Oct2001 Owens Added -tmean option to allow the user to specify
* a date for mean file other than tmid
* 17dec2002 da Silva Bug fix in getgid for pre-projected data; change
* in openf() to detect URLs.
*
*......................................................................
*
* main(args) Main driver checking config, handling error conditions, etc.
*
function main ( args )
_myname = 'lats4d: '
_version = '1.4 of 19 Dec 2002'
* Parse command line
* ------------------
rc = parsecmd(args)
if ( rc != 0 )
if ( _quit = 1 )
say _myname 'exiting from GrADS...'
'quit'
endif
return rc
endif
* Check GrADS configuration
* -------------------------
rc = chkcfg()
if ( rc != 0 )
if ( _quit = 1 )
say _myname 'exiting from GrADS...'
'quit'
endif
return rc
endif
* Output file message
* -------------------
if ( _format = 'stream' )
_fmsg = _format ' file ' _ofname '.bin'
else
if ( _format = 'coards' )
_fmsg = _format ' file ' _ofname '.nc'
else
_fmsg = _format ' file ' _ofname '.grb'
endif
endif
* OK, ready for the real work
* ---------------------------
rc = lats4d ( args )
if ( rc != 0 )
say _myname 'error creating ' _fmsg
else
say _myname 'created ' _fmsg
endif
if ( _quit = 1 )
say _myname 'exiting from GrADS...'
'quit'
endif
return rc
*............................................................................
*
* lats4d(args) The main lats4d entry point.
*
function lats4d ( args )
if ( _verb = 1 )
say _myname 'Version ' _version
endif
* Open input file and reset dimension environment if desired...
* Assumes global grid, i.e., longitudinal wrapping
* -------------------------------------------------------------
if ( _ifname != '' )
'reinit'
rc=openf(_ifname,_ftype)
if ( rc != 0 )
say _myname 'cannot open file "'_ifname'"'
return rc
endif
rc = xyrange()
'set x ' _xmin ' ' _xmax
'set y ' _ymin ' ' _ymax
endif
* Get file number
* ---------------
'q file'
if ( rc != 0 )
say _myname 'no files open; specify "-i" option or open a file'
return rc
endif
_datf = subwrd(result,2)
* Open dimension environment file if so desired,
* otherwise let dim env file be the same as the current data file
* ---------------------------------------------------------------
rc = setdimf()
if ( rc != 0 )
say _myname 'could not open dimension env file: ' _dimenv
return rc
endif
* File information
* ----------------
if ( _verb = 1 )
'q file ' _datf
say _myname 'Data file is '
say result
if ( _datf != _dimf )
'q file ' _dimf
say _myname 'Dimension environment file is '
say result
else
say _myname 'Dimension environment file same as data file'
endif
endif
* Metadata describing origin of the dataset
*------------------------------------------
'set lats model ' _model
'set lats center ' _center
* Parameter table
* ---------------
rc = getvars()
ch = substr(_table,1,1)
* generate parameter table
if ( ch = '@' )
_table = substr(_table,2,256)
rc = mkptab()
if ( rc != 0 ); return 1; endif
endif
* use internal table
if ( ch = '=' )
_table = ''
endif
if ( _table != '' )
'set lats parmtab ' _table
id = subwrd(result,5)
if ( id = 0 ); return 1; endif
endif
* LATS comments (Title for GRIB files)
* ------------------------------------
if ( _title = '' )
if ( _func = '' )
comment = 'File created by GrADS using LATS4D available from '
else
comment = 'File created by GrADS using LATS4D '
comment = comment '(-func ' _func ') available from '
endif
url = 'http://dao.gsfc.nasa.gov/software/grads/lats4d/'
'set lats comment ' comment url
else
'set lats comment ' _title
endif
* Metadata describing conventions/format
* --------------------------------------
if ( _format != 'stream' )
'set lats convention ' _format
if ( rc != 0 )
say _myname 'not really, we stop right here'
return 1
endif
endif
* Get (z,t) range on file
* -----------------------
rc = ztrange()
* Metadata describing time frequency of grids
* -------------------------------------------
if ( _freq = '' )
rc = getfreq()
if ( rc != 0 )
say _myname 'could not determine time frequency; specify -freq option'
return rc
endif
endif
'set lats calendar ' _cal
if ( rc != 0 ); return rc; endif
'set lats frequency ' _freq
if ( rc != 0 ); return rc; endif
if ( _tinc='' ); _tinc=1; endif
odeltat = _tinc * _deltat
'set lats deltat ' odeltat
if ( rc != 0 ); return rc; endif
* Redefine time range (_tmin,_tmax) if user wants to
* --------------------------------------------------
rc = gettim()
if ( rc = 1 )
say _myname 'invalid time range ' _time
return
endif
if ( _verb = 1 )
'q time'
t1 = subwrd(result,3)
t2 = subwrd(result,5)
say _myname 'time range: ' t1 ' ' t2 ' by ' _tinc ', delta t: ' _deltat ' ' _freq
endif
* Define the vertical dimension
* -----------------------------
if ( _levels = '' )
rc = getlevs()
endif
if ( _verb = 1 )
say _myname 'vertical levels: ' _levels
endif
rc = setlevs()
if ( rc = 1 ); return 1; endif
* Set LATS grid type
* ------------------
'set lats gridtype ' _gridtype
if ( rc != 0 ); return rc; endif
* Instruct LATS to use the dimension environment to set the LATS time
* -------------------------------------------------------------------
'set lats timeoption dim_env'
if ( rc != 0 ); return rc; endif
* Define the horizontal grid (based on dim env file)
* --------------------------------------------------
'set dfile ' _dimf
rc = setgrid()
if ( rc != 0 ); return rc; endif
if ( _verb = 1 )
say _myname 'latitudinal range: ' _lat
say _myname 'longitudinal range: ' _lon
endif
_gid = getgid()
if ( _gid < 1 ); return 1; endif
'set dfile ' _datf
* Define the output file name
* ---------------------------
if ( _format = 'stream' )
'set fwrite '_ofname'.bin'
else
'set lats create ' _ofname
_fid = subwrd(result,5)
if ( _fid < 1 ); return 1; endif
endif
* Define all the variables to be written to the LATS file
* -------------------------------------------------------
rc = defvars()
if ( rc != 0 )
say _myname 'fatal error defining variables'
return rc
endif
if ( _verb = 1 )
if ( _func !='' ); say _myname 'Function expression: ' _func; endif
if ( _svars !='' ); say _myname 'surface variables: ' _svars ; endif
if ( _uvars !='' ); say _myname 'upper air variables: ' _uvars ; endif
endif
if ( _svars='' & _uvars='' )
say _myname 'All invalid sfc/upper variables: ' _vars
return 1
endif
* Put GrADS in LATS output mode
* -----------------------------
if ( _format = 'stream' )
'set gxout fwrite'
else
'set gxout latsdata'
if ( rc != 0 ); return rc; endif
endif
* Now, if user wants time mean then redefine time environment
* -----------------------------------------------------------
if ( _mean = 1 )
_tbeg = _tmin
_tend = _tmax
if ( _tinc = '' ); _tinc = 1; endif
_tmid = ( _tmin + _tmax ) / 2
if ( _tmean != '' ); _tmid = _tmean; endif
_tmin = _tmid
_tmax = _tmid
if ( _verb = 1 )
'set t ' _tmid
'q time'
t1 = subwrd(result,3)
say _myname 'time average grid valid at ' t1
say _myname 'averaging increment: ' _tinc ', delta t: ' _deltat ' ' _freq
endif
endif
* The LATS interface only allows one horizontal slice of data
* to be written at a time (z and t dimensions must be fixed),
* so we setup 2 loops
* --------------------------------------------------------------
* Loop over time...
* -----------------
_t = _tmin
while ( _t <= _tmax )
'set t ' _t
if ( rc != 0 ); return rc; endif
'set z ' _zmin
if ( rc != 0 ); return rc; endif
if ( _verb = 1 )
'q time'
time = subwrd(result,3)
say _myname 'writing to ' _fmsg ' on ' time
endif
* Write surface variables
* -----------------------
n = 1
while ( n <= _nsvar )
var = subwrd(_svars,n)
if ( _dimenv != '' )
var = var '.' _datf '(t=' _t ')'
endif
var = subst(var,_func)
if ( _mean = 1 )
var = 'ave('var',t='_tbeg',t='_tend','_tinc')'
endif
vid = subwrd(_svids,n)
if ( _format != 'stream' )
'set lats write ' _fid ' ' vid
if ( rc != 0 ); return rc; endif
endif
'set dfile ' _dimf
'display ' var
if ( rc != 0 ); return rc; endif
'set dfile ' _datf
n = n + 1
endwhile
*
* Loop over upper air variables...
* --------------------------------
n = 1
while ( n <= _nuvar )
var = subwrd(_uvars,n)
if ( _dimenv != '' )
var = var '.' _datf '(t=' _t ')'
endif
var = subst(var,_func)
if ( _mean = 1 )
var = 'ave('var',t='_tbeg',t='_tend','_tinc')'
endif
vid = subwrd(_uvids,n)
* Loop over levels...
* -------------------
k = 1
_level = subwrd(_levels,k)
while ( _level != '' )
_latlevel = subwrd(_latlevels,k)
'set lev ' _level
if ( _format != 'stream' )
'set lats write ' _fid ' ' vid ' ' _latlevel
if ( rc != 0 ); return rc; endif
endif
'set dfile ' _dimf
'display ' var
if ( rc != 0 ); return rc; endif
'set dfile ' _datf
k = k + 1
_level = subwrd(_levels,k)
endwhile
n = n + 1
endwhile
_t = _t + _tinc
endwhile
rc = restdim()
* Close LATS file
* ---------------
if ( _format = 'stream' )
'disable fwrite'
else
'set lats close ' _fid
if ( rc != 0 ); return rc; endif
endif
* All done
* --------
return 0
*.........................................................................
function usage(flag)
say ''
say 'NAME'
say ''
say ' lats4d - LATS for Dummies (Version ' _version ')'
say ''
say 'SYNOPSIS'
say ''
say ' lats4d [-i fn] [-o fn] [-cal calendar] [-center ctr] [-de fn]'
say ' [-format fmt] [-ftype ctl|sdf|xdf] [-freq ...] '
say ' [-func expr] [-h] [-grid type]'
say ' [-lat y1 y2] [-levs ...] [-lon x1 x2] '
say ' [-model mod] [-mean] [-precision nbits] [-table tab] '
say ' [-time t1 t2 [tincr]] [-title ...]'
say ' [-v] [-vars ...] [-xvars] [-zrev] [-q] '
say ''
if ( flag != 1 )
say ' Enter "lats4d -h" for additional information'
say ''
return
endif
say 'DESCRIPTION'
say ''
say ' A minimum fuss gs script for writing NetCDF, HDF-SDS or '
say ' GRIB files from GrADS using the PCMDI LATS interface '
say ' (http://www-pcmdi.llnl.gov). This script can serve as a'
say ' general purpose file conversion and subsetting utility.'
say ' Any GrADS readable file (GrADS IEEE, GSFC Phoenix, GRIB, '
say ' NetCDF or HDF-SDS) can be subset and converted to GRIB, '
say ' NetCDF, HDF-SDS or flat binary (direct access) using a'
say ' single command line.'
say ' '
say ' When invoked without arguments this script will create a'
say ' COARDS compliant NetCDF or HDF-SDS file named '
say ' "grads.lats.nc", with all the contents of the default '
say ' file (all variables, levels, times). The file name and '
say ' several other attributes can be customized at the command'
say ' line, see OPTIONS below.'
say ' '
say ' NetCDF files are obtained by running this script under the'
say ' executable "gradsnc". HDF-SDS files can be produced with'
say ' the "gradshdf" executable. Notice that the classic version'
say ' of grads, "gradsc", does not include support for LATS and'
say ' therefore cannot be used with lats4d. This script requires'
say ' GrADS Version 1.7.beta.9 or later.'
say ''
say 'OPTIONS'
say ''
say ' -i fn input file name; it can be any'
say ' of the following:'
say ' - an ASCII control (ctl) file'
say ' used for GRIB and IEEE files'
say ' - a binary NetCDF file/template'
say ' - a binary HDF-SDS file/template'
say ' - an ASCII data descriptor file (ddf)'
say ' used for non-COARDS compliant'
say ' NetCDF/HDF-SDS files through'
say ' the "xdfopen" command'
say ' If the option "-ftype" is not'
say ' specified lats4d attempts to'
say ' determine the file type using'
say ' a heuristic algorithm.'
say ' NOTE: When the option "-i" is '
say ' specified a GrADS "reinit" is '
say ' issued before the file is opened.'
say ' For NetCDF/HDF-SDS templates in'
say ' GrADS consult the SDFopen home'
say ' page listed under SEE ALSO'
say ''
say ' -o fn output (base) file name; default: '
say ' "grads.lats"'
say ''
say ' -cal calendar calendar type: "standard", "noleap", '
say ' "clim", or "climleap"; default: '
say ' "standard"'
say ''
say ' -center ctr center, e.g., PCMDI, GSFC, NCEP, etc'
say ''
say ' -de fn Dimension environment file name;'
say ' defaut: same as "-i" argument.'
say ' This option is useful for using'
say ' lats4d with the user defined function'
say ' (udf) regrid2. See REGRIDDING below'
say ' for more information.'
say ''
say ' -format fmt LATS file format: coards, grib,'
say ' grads_grib or stream; specify '
say ' "grads_grib" instead of "grib" for'
say ' getting ctl and gribmap files as well.'
say ' NOTE: The option "stream" creates'
say ' a flat binary file using the GrADS'
say ' command "set gxout fwrite" which'
say ' is not part of LATS.'
say ' '
say ''
say ' -ftype ctl|sdf|xdf Specifies the input file type:'
say ' ctl standard GrADS control (ctl)'
say ' file used for IEEE and GRIB '
say ' files'
say ' sdf COARDS compliant NetCDF/HDF-SDS'
say ' binary data file'
say ' xdf data descriptor file (ddf)'
say ' used for non-COARDS compliant'
say ' NetCDF/HDF-SDS files through'
say ' the "xdfopen" command'
say ' By default lats4d attempts to '
say ' determine the file type using a'
say ' heuristic algorithm; use this'
say ' option if lats4d fails to properly'
say ' detect the input file type'
say ' '
say ' -freq [n] unit Time frequency of the input file.'
say ' LATS4D usually detects this from'
say ' the GrADS metadata, but sometimes'
say ' it fails with an error message.'
say ' In such cases use this option.'
say ' Example: -freq 6 hourly '
say ' NOTE: unlike GrADS, LATS does not'
say ' support time frequency in minutes'
say ' Default: n=1, e.g., -freq daily'
say ' '
say ' -func expr Evaluates the expression "expr"'
say ' before writing to the output'
say ' file. The character "@" is used'
say ' to denote the variable name in'
say ' "expr". Example:'
say ' '
say ' -func ave(@,t-1,t+1)'
say ' '
say ' will replace "@" with each '
say ' variable name and produce a file'
say ' with running means. Default:'
say ' expr = @'
say ' '
say ' -grid type Grid type: linear, gaussian or'
say ' generic; default: linear'
say ' '
say ' -h displays this man page'
say ''
say ' -lat y1 y2 latitude range, e.g., "-30 30" for '
say ' 30S thru 30N; default: latitude '
say ' dimension environment'
say ''
say ' -levs lev1 ... levN list of levels; default: all levels'
say ''
say ' -lon x1 x2 longitude range, e.g., "-50 20" for '
say ' 50W thru 20E; default: longitude '
say ' dimension environment'
say ''
say ' -mean saves time mean to file; the actual'
say ' averaging period is specified with'
say ' the "-time" option; the "tincr" '
say ' parameter is the time increment'
say ' for the average (see GrADS ave()'
say ' function)'
say ''
say ' -model mod model name, e.g., GEOS/DAS'
say ''
say ' -precision nbits specify the number of bits of'
say ' precision when storing in GRIB.'
say ' This option is only used when'
say ' lats4d automatically generates'
say ' a parameter table file (see option'
say ' -table below), and the output'
say ' format is "grib" or "grads_grib".'
say ' Default: nbits = 16'
say ''
say ' -table tab LATS parameter table file, e.g., '
say ' "dao.lats.table". If the table name'
say ' starts with "@" (e.g., @my.table)'
say ' then lats4d automatically generates'
say ' a LATS parameter table appropriate'
say ' for the current file and saves it '
say ' to a file; the file name in this'
say ' case is the same as "tab" with the'
say ' @ character removed (e.g., my.table).'
say ' Specify tab as "=" for using the'
say ' internal LATS parameter table.'
say ' See below for additional info on'
say ' parameter tables.'
say ' Default: @.grads.lats.table'
say ''
say ''
say ' -time t1 t2 [tincr] time range and time increment in '
say ' units of the "delta t" in the'
say ' input file; "tincr" is optional;'
say ' Example: "0z1dec91 18z31dec91 2"'
say ' to write every other '
say ' time step'
say ' Defaults: (t1,t2) is taken from '
say ' the time dimension environment,'
say ' and tincr=1. Note: enter "= ="'
say ' for keeping the default values'
say ' for (t1,t2) while specifying tincr'
say ''
say ' -title text output dataset TITLE for GRIB files,'
say ' COMMENTS for NetCDF/HDF files'
say ''
say ' -v verbose mode'
say ''
say ' -vars var1 ... varN list of variables; default: all '
say ' variables on the current file will'
say ' be written to the output file'
say ''
say ' -xsfc exclude all surface variables'
say ''
say ' -xvars var1 ... varN list of variables to exclude;'
say ' default: none'
say ''
say ' -xupper exclude all upper air variables'
say ''
say ' -zrev reverse order of vertical levels'
say ''
say ' -q quits GrADS upon return'
say ' '
say 'LATS PARAMETER TABLES'
say ' '
say ' LATS maintains an internal parameter table that prescribes'
say ' variable names, description, units, datatype, basic'
say ' structure (e.g., upper air or surface), and compression'
say ' (GRIB options). These descriptors are inferred from the'
say ' parameter name only, and thus most of the metadata needed'
say ' to write GRIB and/or netCDF data are located in the'
say ' parameter table and need not be specified at the command'
say ' line. The option "-table" is provided to override the '
say ' internal table with an external parameter file. For'
say ' additional information on LATS parameter tables'
say ' consult http://www-pcmdi.llnl.gov/software/lats/.'
say ''
say ' The only inconvenience of this approach is that variables'
say ' names being written to file must match those defined in '
say ' this internal parameter table (which is essentially the '
say ' same as the "AMIPS2" LATS table, see URL above).'
say ' To circumvent this problem lats4d can automatically'
say ' generate a parameter table based on the current file'
say ' metadata. Since GrADS keeps no units or GRIB packing'
say ' information, this parameter file sets the units entry'
say ' to blank and uses defaults for the GRIB packing parameters.'
say ' The default GRIB packing algorithm is "16-bit fixed width '
say ' compression" and produces GRIB files which are about half'
say ' the size of NetCDF/HDF-SDS files. The option "-precision"'
say ' allows the user to define the number of bits of precision'
say ' at the command line; see EXAMPLES ex2a,b,c below.'
say ' If you care about having proper metadata written to'
say ' your file or need more efficient GRIB packing then you can '
say ' either change your variable names to match those in the '
say ' internal LATS table, or customize an existing LATS parameter'
say ' table; see URL above for sample parameter tables.'
say ' '
say 'LATS QUALITY CONTROL WARNINGS'
say ' '
say ' Quality control (QC) information is included in some '
say ' LATS parameter tables to help the user ensure that their'
say ' data is being written properly. In such cases, if LATS'
say ' detects suspect data it writes a warning message to the'
say ' screen and saves additional information in a log file.'
say ' Consult the LATS home page for additional information.'
say ' '
say 'REGRIDDING'
say ' '
say ' This script can be used with Mike Fiorino s user'
say ' defined function (udf) regrid2(). This combination'
say ' allows you to convert any GrADS redable file to any'
say ' other horizontal resolution/domain of your choice. '
say ' Here is a quick roadmap:'
say ' 1. Start by installing regrid2() available from'
say ' ftp://grads.iges.org/grads/sprite/udf/regrid2beta.tar'
say ' 2. If you already have a sample file at the desired new'
say ' resolution, great! Otherwise you can get one by creating '
say ' a fake GrADS control file. There are a few samples'
say ' on the last4d home page: geos1x1.ctl, geos4x5.ctl and'
say ' geos2x25.ctl. This file is used to define the dimension'
say ' environment at the new desired resolution through the'
say ' "-de" option.'
say ' 3. Here is an example which converts the sample model.???'
say ' data file from 4x5 (latxlon) resolution to 1x1:'
say ' '
say ' lats4d -i model -de geos1x1 -func regrid2(@,1,1,bs_p1,-180,-90)'
say ' '
say ' The resulting "grads.lats.nc" file is at 1x1 degree'
say ' resolution.'
say ' '
say 'EXAMPLES'
say ''
say ' Download files "model.ctl", "model.gmp" and "model.grb"'
say ' from http://dao.gsfc.nasa.gov/software/grads/lats4d/.'
say ' Then start "gradsnc" or "gradshdf" and try these,'
say ' carefully examining the files produced:'
say ' '
say ' lats4d -h'
say ' lats4d -v -q -i model -o ex1 '
say ' lats4d -v -q -i model -o ex2a -format grads_grib'
say ' lats4d -v -q -i model -o ex2b -format grads_grib -precision 8'
* say ' lats4d -v -q -i model -o ex2c -format grads_grib -precision 32'
say ' lats4d -v -q -i model -o ex3 -levs 700 500 -vars ua va'
say ' lats4d -v -q -i model -o ex4 -time 1jan1987 3jan1987'
say ' lats4d -v -q -i model -o ex5 -time = = 2'
say ' lats4d -v -q -i model -o ex6 -mean'
say ' lats4d -v -q -i model -o ex7 -mean -time = = 2'
say ' lats4d -v -q -i model -o ex8 -lat 20 70 -lon -140 -60'
say ''
say ' Note: the "-q" option above is to make sure you'
say ' restart GrADS; see BUGS below. You may want to'
say ' enter these from your OS shell, e.g.,'
say ''
say ' % gradsnc -blc "lats4d -v -q -i model -o ex1"'
say ''
say ' The sh(1) script "lats4d" allows you to enter lats4d'
say ' options directly from the Unix command line, e.g.,'
say ''
say ' % lats4d -v -i model -o ex1 '
say ''
say 'BUGS'
say ''
say ' Sometimes lats4d will only work if you exit and'
say ' restart GrADS.'
say ''
say ' The option "-precision 32" does not quite work. This'
say ' appears to be a LATS bug.'
say ''
say ' Because of a limitation in the GRIB format, "grib" or '
say ' "grads_grib" output cannot have levels where p<1.'
say ' To circumvent this problem, a hybrid level number is'
say ' is used in such cases.'
say ''
say 'SEE ALSO '
say ''
say ' GrADS http://grads.iges.org/grads/ '
say ' LATS http://www-pcmdi.llnl.gov/software/lats'
say ' LATS4D http://dao.gsfc.nasa.gov/software/grads/lats4d'
say ' SDFopen http://www.cdc.noaa.gov/~hoop/grads.html'
say ' XDFopen http://www.cdc.noaa.gov/~hoop/xdfopen.shtml'
say ' NetCDF http://www.unidata.ucar.edu/packages/netcdf/'
say ' HDF http://hdf.ncsa.uiuc.edu/'
say ' GRIB ftp://ncardata.ucar.edu/docs/grib/prev-vers/guide.txt'
say ' http://www.wmo.ch/web/www/reports/Guide-binary-2.html'
say ' '
say 'COPYRIGHT'
say ''
say ' Copyright (c) 1998-1999 A. da Silva'
say ' Permission is granted to copy, modify and distribute this'
say ' software provided it is not sold for profit, and provided '
say ' this notice is included. '
say ''
say 'NO WARRANTY'
say ' '
say ' Because lats4d is provided free of charge, it is provided'
say ' "as is" WITHOUT WARRANTY OF ANY KIND, either expressed or'
say ' implied, including, but not limited to, the implied'
say ' warranties of merchantability and fitness for a particular'
say ' purpose. USE AT YOUR OWN RISK.'
say ' '
say 'CREDITS'
say ''
say ' Arlindo da Silva (NASA/GSFC) wrote the lats4d.gs script. '
say ' Mike Fiorino (PCMDI/LLNL) wrote the LATS interface to'
say ' GrADS. Robert Drach, Mike Fiorino and Peter Gleckler'
say ' (PCMDI/LLNL) wrote the LATS library.'
say ''
return 1
*.........................................................................
*
* parsecmd() Parse command line arguments
*
function parsecmd(args)
*
* Note: customize defaults for your site
*
_ifname = ''
_ofname = 'grads.lats'
_format = 'coards'
_model = 'geos/das'
_center = 'gsfc'
_table = '@.grads.lats.table'
_precision = 16
_vars = ''
_xvars = ''
_xsfc = ''
_xupper = ''
_func = ''
_dimenv = ''
_title = ''
_lat = ''
_lon = ''
_levels = ''
_time = ''
_tmean = ''
_tinc = ''
_freq = ''
_ftype = ''
_gridtype = 'linear'
_help = 0
_verb = 0
_quit = 0
_mean = 0
_zrev = 0
_cal = 'standard'
options = '-model -center -table -v -q -i -o -format -vars -levs -time'
options = options ' -h -help -lat -lon -cal -mean -precision -ftype -tmean'
options = options ' -grid -func -zrev -de -title -xvars -xsfc -xupper'
options = options ' -geos1x1 -geos4x5 -geos2x25'
options = options ' -geos1x1a -geos4x5a -geos2x25a'
options = options ' -geos1x1s -geos4x5s -geos2x25s'
i = 1
token = subwrd(args,i)
while ( token != '' )
* Handle each option separately ...
if ( token = '-help' | token = '-h' )
_help = 1
i = i + 1
token = subwrd(args,i)
endif
if ( token = '-v' )
_verb = 1
i = i + 1
token = subwrd(args,i)
endif
if ( token = '-q' )
_quit = 1
i = i + 1
token = subwrd(args,i)
endif
if ( token = '-mean' )
_mean = 1
i = i + 1
token = subwrd(args,i)
endif
if ( token = '-xsfc' )
_xsfc = 1
i = i + 1
token = subwrd(args,i)
endif
if ( token = '-xupper' )
_xupper = 1
i = i + 1
token = subwrd(args,i)
endif
if ( token = '-zrev' )
_zrev = 1
i = i + 1
token = subwrd(args,i)
endif
if ( token = '-i' )
i = i + 1
token = subwrd(args,i)
first = x '' substr(token,1,1)
while ( first != 'x-' & first != 'x' )
_ifname= _ifname ' ' token
i = i + 1
token = subwrd(args,i)
first = x '' substr(token,1,1)
endwhile
* allows sdfopen templates
if ( critique(1,3,'-i',_ifname) = 1 ); return 1; endif
endif
if ( token = '-o' )
_ofname = subwrd(args,i+1)
i = i + 2
token = subwrd(args,i)
if ( critique(1,1,'-o',_ofname) = 1 ); return 1; endif
endif
if ( token = '-cal' )
_cal = subwrd(args,i+1)
i = i + 2
token = subwrd(args,i)
if ( critique(1,1,'-cal',_cal) = 1 ); return 1; endif
endif
if ( token = '-format' )
_format = subwrd(args,i+1)
i = i + 2
token = subwrd(args,i)
if ( critique(1,1,'-format',_format) = 1 ); return 1; endif
if ( _format = 'fwrite' ); _format = 'stream'; endif
endif
if ( token = '-grid' )
_gridtype = subwrd(args,i+1)
i = i + 2
token = subwrd(args,i)
if ( critique(1,1,'-grid',_gridtype) = 1 ); return 1; endif
endif
if ( token = '-ftype' )
_ftype= subwrd(args,i+1)
i = i + 2
token = subwrd(args,i)
if ( critique(1,1,'-ftype',_format) = 1 ); return 1; endif
if ( _ftype='nc'|_ftype='netcdf'|_ftype='hdf'|_ftype='hdf-sds'|_ftype='gfio' )
_ftype = 'sdf'
endif
if ( _ftype='grib'|_ftype='grads_grib'|_ftype='ieee'|_ftype='sequential')
_ftype = 'ctl'
endif
endif
if ( token = '-model' )
_model = subwrd(args,i+1)
i = i + 2
token = subwrd(args,i)
if ( critique(1,1,'-model',_model) = 1 ); return 1; endif
endif
if ( token = '-center' )
_center = subwrd(args,i+1)
i = i + 2
token = subwrd(args,i)
if ( critique(1,1,'-center',_center) = 1 ); return 1; endif
endif
if ( token = '-de' )
_dimenv = subwrd(args,i+1)
i = i + 2
token = subwrd(args,i)
if ( critique(1,1,'-de',_dimenv) = 1 ); return 1; endif
endif
if ( token = '-table' )
_table = subwrd(args,i+1)
i = i + 2
token = subwrd(args,i)
if ( critique(1,1,'-table',_table) = 1 ); return 1; endif
endif
if ( token = '-precision' )
_precision = subwrd(args,i+1)
i = i + 2
token = subwrd(args,i)
if ( critique(1,1,'-precision',_precision) = 1 ); return 1; endif
endif
if ( token = '-func' )
_func= subwrd(args,i+1)
i = i + 2
token = subwrd(args,i)
if ( critique(1,1,'-func',_func) = 1 ); return 1; endif
endif
if ( token = '-vars' )
i = i + 1
token = subwrd(args,i)
first = x '' substr(token,1,1)
while ( first != 'x-' & first != 'x' )
_vars = _vars ' ' token
i = i + 1
token = subwrd(args,i)
first = x '' substr(token,1,1)
endwhile
if ( critique(1,999,'-vars',_vars) = 1 ); return 1; endif
endif
if ( token = '-xvars' )
i = i + 1
token = subwrd(args,i)
first = x '' substr(token,1,1)
while ( first != 'x-' & first != 'x' )
_xvars = _xvars ' ' token
i = i + 1
token = subwrd(args,i)
first = x '' substr(token,1,1)
endwhile
if ( critique(1,999,'-xvars',_xvars) = 1 ); return 1; endif
endif
if ( token = '-title' )
i = i + 1
token = subwrd(args,i)
first = x '' substr(token,1,1)
while ( first != 'x-' & first != 'x' )
_title = _title ' ' token
i = i + 1
token = subwrd(args,i)
first = x '' substr(token,1,1)
endwhile
if ( critique(1,999,'-title',_title) = 1 ); return 1; endif
endif
*
* Note: we force -lon/-lat to have 2 args because lats do
* not allow single point grids.
*
if ( token = '-lat' )
tok1 = subwrd(args,i+1)
tok2 = subwrd(args,i+2)
_lat = tok1 ' ' tok2
i = i + 3
token = subwrd(args,i)
if ( critique(2,2,'-lat',_lat) = 1 ); return 1; endif
endif
if ( token = '-lon' )
tok1 = subwrd(args,i+1)
tok2 = subwrd(args,i+2)
_lon = tok1 ' ' tok2
i = i + 3
token = subwrd(args,i)
if ( critique(2,2,'-lon',_lon) = 1 ); return 1; endif
endif
if ( token = '-levs' )
i = i + 1
token = subwrd(args,i)
first = x '' substr(token,1,1)
while ( first != 'x-' & first != 'x' )
_levels = _levels ' ' token
i = i + 1
token = subwrd(args,i)
first = x '' substr(token,1,1)
endwhile
endif
if ( token = '-time' )
i = i + 1
token = subwrd(args,i)
first = x '' substr(token,1,1)
while ( first != 'x-' & first != 'x' )
_time = _time ' ' token
i = i + 1
token = subwrd(args,i)
first = x '' substr(token,1,1)
endwhile
if ( critique(1,3,'-time',_time) = 1 ); return 1; endif
_tinc = subwrd(_time,3)
t1 = subwrd(_time,1)
t2 = subwrd(_time,2)
if ( _tinc != '' )
_time = t1 ' ' t2
else
_tinc = 1
endif
if ( t1 = '=' | t2 = '=' ); _time = ''; endif
if ( _tinc < 1 )
rc = usage()
say _myname 'invalid time increment ' _tinc
return 1
endif
endif
if ( token = '-freq' )
_deltat = subwrd(args,i+1)
_freq = subwrd(args,i+2)
token = _freq ' ' _deltat
if ( critique(1,2,'-freq',token) = 1 ); return 1; endif
first = substr(_freq,1,1)
if ( first = '-' | first = '' )
_freq = _deltat
_deltat = 1
i = i + 2
else
i = i + 3
endif
token = subwrd(args,i)
endif
* Undocumented options (DAO specific)
* -----------------------------------
if ( token = '-geos1x1' )
rc = mkgeosf('geos1x1.ctl',1,1,360,181,'bl_p1')
if ( rc != 0 ); return 1; endif
i = i + 1
token = subwrd(args,i)
endif
if ( token = '-geos1x1a' )
rc = mkgeosf('geos1x1.ctl',1,1,360,181,'ba_p1')
if ( rc != 0 ); return 1; endif
i = i + 1
token = subwrd(args,i)
endif
if ( token = '-geos1x1s' )
rc = mkgeosf('geos1x1.ctl',1,1,360,181,'bs_p1')
if ( rc != 0 ); return 1; endif
i = i + 1
token = subwrd(args,i)
endif
if ( token = '-geos2x2.5' )
rc = mkgeosf('geos2x25.ctl',2.5,2,144,91,'bl_p1')
if ( rc != 0 ); return 1; endif
i = i + 1
token = subwrd(args,i)
endif
if ( token = '-geos2x2.5a' )
rc = mkgeosf('geos2x25.ctl',2.5,2,144,91,'ba_p1')
if ( rc != 0 ); return 1; endif
i = i + 1
token = subwrd(args,i)
endif
if ( token = '-geos2x2.5s' )
rc = mkgeosf('geos2x25.ctl',2.5,2,144,91,'bs_p1')
if ( rc != 0 ); return 1; endif
i = i + 1
token = subwrd(args,i)
endif
if ( token = '-geos4x5' )
rc = mkgeosf('geos4x5.ctl',5,4,72,46,'bl_p1')
if ( rc != 0 ); return 1; endif
i = i + 1
token = subwrd(args,i)
endif
if ( token = '-geos4x5a' )
rc = mkgeosf('geos4x5.ctl',5,4,72,46,'ba_p1')
if ( rc != 0 ); return 1; endif
i = i + 1
token = subwrd(args,i)
endif
if ( token = '-geos4x5s' )
rc = mkgeosf('geos4x5.ctl',5,4,72,46,'bs_p1')
if ( rc != 0 ); return 1; endif
i = i + 1
token = subwrd(args,i)
endif
if ( token = '-tmean' )
_tmean = subwrd(args,i+1)
i = i + 2
token = subwrd(args,i)
endif
* Validate option
* ---------------
j = 1
opt = subwrd(options,j)
valid = 0
while ( opt != '' )
if ( token = opt ); valid = 1; endif
j = j + 1
opt = subwrd(options,j)
endwhile
if ( valid = 0 & token != '' )
rc = usage()
say _myname 'invalid option "'token'"'
return 1
endif
endwhile
if ( _help = 1 )
rc = usage(1)
return 1
endif
return 0
function printopt()
say ' Input fname: ' _ifname
say ' Output fname: ' _ofname
say ' Format: ' _format
say ' Variables: ' _vars
say ' Time Range: ' _time
say ' Levels: ' _levels
return
*.......................................................................
*
* openf(fname,ftype) Opens a file according to the file type (ftype).
* If ftype is not specified it attempts to determine it
* by a heuristic algorithm;
* ftype can be 'ctl', 'sdf' or 'xdf'
*
function openf(fname,ftype)
* Determine file type
* -------------------
if ( ftype = '' | ftype ='ftype' )
* fname may be a template...
* filen is always a file name
filen = subwrd(fname,1)
http = substr(filen,1,7)
if ( http = 'http://' )
ftype = 'sdf'
else
buf = read(filen)
rc = sublin(buf,1)
if ( rc != 0 )
buf = read(filen'.ctl')
rc = sublin(buf,1)
endif
if ( rc != 0 )
say _myname 'cannot read file ' filen ' or ' filen '.ctl'
return rc
endif
rec = sublin(buf,2)
tok = subwrd(rec,1)
if ( tok = 'dset' | tok='DSET' )
is_xdf = 0
i = 1
tok = substr(filen,i,4)
while ( tok != '' )
if ( tok='.ddf' | tok='.DDF' ); is_xdr = 1; endif
i = i + 1
tok = substr(filen,i,4)
endwhile
if ( is_xdr = 1 )
ftype = 'xdf'
else
ftype = 'ctl'
endif
else
ftype = 'sdf'
endif
endif
endif
* Open according to file type
* ---------------------------
if ( ftype = 'ctl' )
'open ' fname
_result = result
return rc
endif
if ( ftype = 'sdf' )
'sdfopen ' fname
_result = result
return rc
endif
if ( ftype = 'xdf' )
'xdfopen ' fname
_result = result
return rc
endif
say _myname 'cannot handle file type "' ftype '"'
return 1
*
* setdimf - Open dimension env file
*
function setdimf()
_dimf = 0
if ( _dimenv != '' )
rc=openf(_dimenv,'')
if ( rc != 0 )
say _myname 'cannot open file "'_dimenv'"'
return rc
endif
result = _result
i = 0
while ( token != '' )
i = i + 1
lin = sublin(result,i)
token = subwrd(lin,1)
word = subwrd(lin,2)
if ( token = 'Data' & word = 'file' )
_dimf = subwrd(lin,8)
endif
endwhile
if ( _dimf = 0 ); return 1; endif
'set dfile ' _dimf
rc = xyrange()
'set x ' _xmin ' ' _xmax
'set y ' _ymin ' ' _ymax
'set dfile ' _datf
else
_dimf = _datf
endif
return 0
*
* gettim() Redefines time range according to user defined _time
*
function gettim()
if ( _time = '' ); return 0; endif
tbeg = subwrd(_time,1)
tend = subwrd(_time,2)
'set time ' tbeg ' ' tend
if ( rc = 1 ); return 1; endif
rc = savedim()
if ( _t1save < _tmin | _t2save > _tmax )
return 1
else
_tmin = _t1save
_tmax = _t2save
endif
return 0
*
* getfreq() Examines two subsequent times and determine
* time increment. This is heuristic and not
* guaranteed to always work.
*
function getfreq()
* Special case: only 1 time step on file
if ( _tmax = 1 )
_freq = hourly
_deltat = 6
say _myname" only 1 time step on file, assuming frequency is " _deltat ' ' _freq
return 0
endif
*
* Get year-month-day-hour for 2 time intervals
*
'q time'
tb = subwrd(result,3)
te = subwrd(result,5)
'set t 1 2'
'q time'
t1 = subwrd(result,3)
t2 = subwrd(result,5)
'set time ' tb ' ' te
rc = split( t1 )
y1=_y; m1=_m; d1=_d; h1=_h
rc = split( t2 )
y2=_y; m2=_m; d2=_d; h2=_h
*
* Minutes: not supported
*
if ( y1=y2 & m1=m2 & d1=d2 & h1=h2 )
say _myname 'LATS does not support time frequency in minutes'
return 1
endif
*
* Hourly
*
if ( y1=y2 & m1=m2 & d1=d2 )
_freq = hourly
_deltat = h2 - h1
if ( _deltat <= 0 ); return 1; endif
return 0
endif
*
* Daily
*
if ( y1=y2 & m1=m2 )
_freq = daily
_deltat = d2 - d1
if ( _deltat <= 0 ); return 1; endif
return 0
endif
*
* Monthly
*
if ( y1=y2 )
_freq = monthly
_deltat = m2 - m1
if ( _deltat <= 0 ); return 1; endif
return 0
endif
*
* Yearly
*
if ( y1 < y2 )
_freq = yearly
_deltat = y2 - y1
return 0
endif
return 1
*
* split() Returns year, month, day & hour from date of the form
* 00Z02JAN1987
*
function split ( t )
ch = substr(t,3,1)
if ( ch = ':' )
off = 3
else
off = 0
endif
_h = substr(t,1,2)
_d = substr(t,4+off,2)
m3 = substr(t,6+off,3)
_y = substr(t,9+off,4)
mons = 'JAN FEB MAR APR MAY JUN JUL AUG SEP OCT NOV DEC'
i = 1
while(i<=12)
mon = subwrd(mons,i)
if ( m3 = mon )
_m = i
return
endif
i = i + 1
endwhile
return
*
* ztrange() Get (z,t) range on file
*
function ztrange()
'q file'
tmp = sublin ( result, 5 )
_zmin = 1
_zmax = subwrd(tmp,9)
_tmin = 1
_tmax = subwrd(tmp,12)
return
function xyrange()
'q file'
tmp = sublin ( result, 5 )
_xmin = 1
_xmax = subwrd(tmp,3)
_ymin = 1
_ymax = subwrd(tmp,6)
return
*
* getlevs() Get all levels from dimension environment
*
function getlevs()
_levels = ''
if ( _zrev = 1 )
z = _zmax
while ( z >= _zmin )
'set z ' z
lev = subwrd(result,4)
_levels = _levels ' ' lev
z = z - 1
endwhile
else
z = _zmin
while ( z <= _zmax )
'set z ' z
lev = subwrd(result,4)
_levels = _levels ' ' lev
z = z + 1
endwhile
endif
return
*
* setlevs() Set LATS vertical dimension.
*
function setlevs()
* GRIB does not allow pressure < 1 hPa
* ------------------------------------
plev = 1
zlevs = ''
if ( _format = 'grads_grib' | _format = 'grib' )
k = 1
lev = subwrd(_levels,k)
while ( lev != '' )
if ( lev < 1 )
say _myname 'invalid plev ' lev ' for GRIB output'
plev = 0
endif
'set lev ' lev
'q dims'
tmp = sublin(result,4)
zlevs = zlevs ' ' subwrd(tmp,9)
k = k + 1
lev = subwrd(_levels,k)
endwhile
endif
* Set pressure levels
* -------------------
if ( plev = 1 )
if ( _verb = 1 )
say _myname 'using PRESSURE for vertical coordinate'
endif
_latlevels = _levels
'set lats vertdim plev ' _latlevels
_lid = subwrd(result,5)
if ( _lid < 1 ); return 1; endif
_sid = 0
* Set hybrid levels
* -----------------
else
if ( _verb = 1 )
say _myname 'using HYBRID level number for vertical coordinate'
endif
_latlevels = zlevs
'set lats vertdim hybrid ' _latlevels
_lid = subwrd(result,5)
if ( _lid < 1 ); return 1; endif
_sid = 0
endif
return 0
*
* setgrid() Sets the horizontal grid if user specified
* -lat or -lon at the command line
*
function setgrid()
rc = savedim()
if ( _lon != '' )
'set lon ' _lon
if ( rc!=0 ); return rc; endif
else
_lon = _lonmin ' ' _lonmax
endif
if ( _lat != '' )
'set lat ' _lat
if ( rc!=0 ); return rc; endif
else
_lat = _latmin ' ' _latmax
endif
return 0
*
* getgid() Get horizontal grid id
*
function getgid()
'q file'
tmp = sublin(result,7)
var = subwrd(tmp,1)
rc = savedim()
'set z 1'
'set t 1'
'set gxout latsgrid'
'd ' var
bufr = sublin(result,1)
grid = subwrd(bufr,2)
if ( grid != 'GRID' )
bufr = sublin(result,2)
grid = subwrd(bufr,2)
if ( grid != 'GRID' ); return -1; endif
endif
gid = subwrd(bufr,5)
rc = restdim()
return gid
*
* savedim() Save dimension environment
*
function savedim()
'q dims'
tmp = sublin(result,2)
type = subwrd(tmp,3)
if ( type = 'varying' )
_x1save = subwrd(tmp,11)
_x2save = subwrd(tmp,13)
_lonmin = subwrd(tmp,6)
_lonmax = subwrd(tmp,8)
else
_x1save = subwrd(tmp,9)
_x2save = _x1save
_lonmin = subwrd(tmp,6)
_lonmax = _lonmin
endif
tmp = sublin(result,3)
type = subwrd(tmp,3)
if ( type = 'varying' )
_y1save = subwrd(tmp,11)
_y2save = subwrd(tmp,13)
_latmin = subwrd(tmp,6)
_latmax = subwrd(tmp,8)
else
_y1save = subwrd(tmp,9)
_y2save = _y1save
_latmin = subwrd(tmp,6)
_latmax = _latmin
endif
tmp = sublin(result,4)
type = subwrd(tmp,3)
if ( type = 'varying' )
_z1save = subwrd(tmp,11)
_z2save = subwrd(tmp,13)
else
_z1save = subwrd(tmp,9)
_z2save = _z1save
endif
tmp = sublin(result,5)
type = subwrd(tmp,3)
if ( type = 'varying' )
_t1save = subwrd(tmp,11)
_t2save = subwrd(tmp,13)
else
_t1save = subwrd(tmp,9)
_t2save = _t1save
endif
return
*
* restdim() Restore saved dimension environment
*
function restdim()
if ( _x1save = '_x1save' )
say 'restdim: dimensions not saved'
return 1
endif
'set x ' _x1save ' ' _x2save
'set y ' _y1save ' ' _y2save
'set z ' _z1save ' ' _z2save
'set t ' _t1save ' ' _t2save
return
*
* getvars() Get all variables from current file
*
function getvars()
'q file'
tmp = sublin(result,6)
nvars = subwrd(tmp,5)
n = 1
_svars = ''; _slong=''; _nsvar = 0
_uvars = ''; _ulong=''; _nuvar = 0
while ( n <= nvars )
tmp = sublin(result,6+n)
var = subwrd(tmp,1)
long = ''
i = 4
token = subwrd(tmp,i)
while ( token != '' )
long = long ' ' token
i = i + 1
token = subwrd(tmp,i)
endwhile
nlevs = subwrd(tmp,2)
rc = validate(var,_vars)
if ( rc = 1 )
rc = xvalidat(var,_xvars)
if ( rc = 0 )
if ( _verb = 1 )
say _myname 'excluding variable ' var
endif
endif
endif
if ( rc = 1 )
if ( nlevs < 2 )
_nsvar = _nsvar + 1
_svars = _svars ' ' var
j = _nsvar; _slong.j = long
if ( _verb = 1 )
*** say _myname 'selecting sfc var ' var ': ' _slong.j
endif
else
_nuvar = _nuvar + 1
_uvars = _uvars ' ' var
j = _nuvar; _ulong.j = long
if ( _verb = 1 )
*** say _myname 'selecting up var ' var ': ' _ulong.j
endif
endif
endif
n = n + 1
endwhile
* Exclude all surface variables if user so chooses
* ------------------------------------------------
if ( _xsfc = 1 )
_nsvar = 0
_svars = ''
endif
* Exclude all upper air variables if user so chooses
* --------------------------------------------------
if ( _xupper = 1 )
_nuvar = 0
_uvars = ''
endif
return
*
* mkptab() Creates a LATS parameter table from variable list.
*
function mkptab()
*
* Temporary file
*
fname = _table
if ( _verb = 1 )
say _myname 'creating LATS PARAMETER TABLE file ' fname
endif
*
* Identify ourselves...
*
rc=write(fname,'# LATS PARAMETER TABLE automatically generated by lats4d')
rc=subwrd(rc,1)
if ( rc > 0 )
say rc
say _myname 'cannot create LATS table file ' fname ' on current directory'
return rc
endif
rc=write(fname,'# Note: This table lacks the UNITS and GRIB decimal_scale_factor columns.', append)
rc=write(fname,'# It also does not include QC information. For additional', append)
rc=write(fname,'# information on LATS parameter tables consult:' , append)
rc=write(fname,'# http://www-pcmdi.llnl.gov/software/lats/' , append)
*
* General info
*
rc=write(fname,'#---------------------------------------------------------------------------------------------------',append)
rc=write(fname,'#',append)
#rc=write(fname,'# A parameter file is divided into sections, indicated by #! comments.',append)
rc=write(fname,'# The sections may appear in any order. The 'center' section is only required for GRIB output.',append)
rc=write(fname,'#',append)
rc=write(fname,'# #!variable',append)
rc=write(fname,'#',append)
rc=write(fname,'# Variable table: defines variable-specific parameters',append)
rc=write(fname,'#',append)
rc=write(fname,'# #!vert',append)
rc=write(fname,'#',append)
rc=write(fname,'# Vertical dimension type table: defines categories of vertical dimensions',append)
rc=write(fname,'#',append)
rc=write(fname,'# #!center',append)
rc=write(fname,'#',append)
rc=write(fname,'# Center table: defines GRIB parameters which identify the originating process, center, and subcenter.',append)
rc=write(fname,'#',append)
rc=write(fname,'# #!qc',append)
rc=write(fname,'#',append)
rc=write(fname,'# Quality control marks table: defines the values controlling the quality control routines.',append)
rc=write(fname,'# ',append)
rc=write(fname,'#---------------------------------------------------------------------------------------------------',append)
*
* Variable list
*
rc=write(fname,'#!variable',append)
rc=write(fname,'#',append)
rc=write(fname,'# Variables',append)
rc=write(fname,'# (max number of entries = LATS_MAX_PARMS in lats.h)',append)
rc=write(fname,'#',append)
rc=write(fname,'# The format of each record is:',append)
rc=write(fname,'# name | id | title | units | datatype | surface | decimal_scale_factor | precision | comments_1 | comments_2',append)
rc=write(fname,'#',append)
rc=write(fname,'# name = variable name (no blanks)',append)
rc=write(fname,'# id = GRIB parameter number (>127 => AMIP-2 specific)',append)
rc=write(fname,'# title = long name (description)',append)
rc=write(fname,'# units = variable units',append)
rc=write(fname,'# datatype = 'float' or 'int'',append)
rc=write(fname,'# level_type = level_type in vertical dimension table, or blank if values must be defined via lats_vert_dim',append)
rc=write(fname,'# decimal_scale_factor = GRIB decimal scale factor, or -999 if no decimal scaling',append)
rc=write(fname,'# precision = number of bits of precision if stored in GRIB,',append)
rc=write(fname,'# or -999 for level-dependent bit length (ignored if decimal_scale_factor is set)',append)
rc=write(fname,'# comments_1 = comments, ignored by LATS',append)
rc=write(fname,'# comments_2 = comments, ignored by LATS',append)
rc=write(fname,'#',append)
rc=write(fname,'#---------------------------------------------------------------------------------------------------',append)
* Fake these for now
units = ' '
scale_factor = -999
precision = _precision
dattype = 'float'
* Sfc variables
id = 1
i = 1
levtype = 'sfc'
while(i <= _nsvar )
name = subwrd(_svars,i) ' '
title = _slong.i ' '
name = substr(name,1,8)
title = substr(title,1,70)
idc = substr(id' ',1,2)
record = name '| ' idc ' |' title ' | ' units ' | ' dattype ' | ' levtype ' | ' scale_factor ' | ' precision ' | | |'
rc = write(fname,record,append)
i = i + 1
id = id + 1
endwhile
* Upper air variables
i = 1
levtype = ' '
while(i <= _nuvar )
name = subwrd(_uvars,i) ' '
title = _ulong.i ' '
name = substr(name,1,8)
title = substr(title,1,70)
idc = substr(id' ',1,2)
record = name '| ' idc ' |' title ' | ' units ' | ' dattype ' | ' levtype ' | ' scale_factor ' | ' precision ' | | |'
rc = write(fname,record,append)
i = i + 1
id = id + 1
endwhile
*
* Vertical dimension
*
rc=write(fname,'#---------------------------------------------------------------------------------------------------',append)
rc=write(fname,'#! vert',append)
rc=write(fname,'# Vertical dimension types',append)
rc=write(fname,'# (max number of entries = LATS_MAX_VERT_TYPES in lats.h)',append)
rc=write(fname,'#',append)
rc=write(fname,'# The format of each record is:',append)
rc=write(fname,'# level_type | description | units | verticality | positive | default | GRIB_id | GRIB_p1 | GRIB_p2 | GRIB_p3',append)
rc=write(fname,'#',append)
rc=write(fname,'# level_type = level type',append)
rc=write(fname,'# description = level description',append)
rc=write(fname,'# units = units for this level type',append)
rc=write(fname,'# verticality = 'single' (single surface) or 'multi' (variable can have multiple levels of this type)',append)
rc=write(fname,'# positive = 'up' (increasing values point up) or 'down' (increasing values point down)',append)
rc=write(fname,'# GRIB_id = GRIB level type indicator (PDS octet 10)',append)
rc=write(fname,'# GRIB_p1 = GRIB PDS octet 11',append)
rc=write(fname,'# GRIB_p2 = GRIB PDS octet 12',append)
rc=write(fname,'# GRIB_p3 = combined GRIB octets 11, 12 - overrides values of GRIB_p1, GRIB_p2 if specified',append)
rc=write(fname,'',append)
rc=write(fname,'0degiso | 0 deg isotherm | hPa | single | up | 4 | 0 | 0 | 0',append)
rc=write(fname,'atm | Atmosphere (entire) | | single | up | 200 | 0 | 0 | 0 ',append)
rc=write(fname,'ocn | Ocean (entire depth) | | single | up | 201 | 0 | 0 | 0 ',append)
rc=write(fname,'clhbot | High Cloud Bottom Level | hPa | single | up | 232 | 0 | 0 | 0',append)
rc=write(fname,'clhlay | High Cloud Top Layer | | single | up | 234 | 0 | 0 | 0',append)
rc=write(fname,'clhtop | High Cloud Top Level | hPa | single | up | 233 | 0 | 0 | 0',append)
rc=write(fname,'cllbot | Low Cloud Bottom Level | hPa | single | up | 212 | 0 | 0 | 0',append)
rc=write(fname,'clllay | Low Cloud Top Layer | | single | up | 214 | 0 | 0 | 0',append)
rc=write(fname,'clltop | Low Cloud Top Level | hPa | single | up | 213 | 0 | 0 | 0',append)
rc=write(fname,'clmbot | Mid Cloud Bottom Level | hPa | single | up | 222 | 0 | 0 | 0',append)
rc=write(fname,'clmlay | Mid Cloud Top Layer | | single | up | 224 | 0 | 0 | 0',append)
rc=write(fname,'clmtop | Mid Cloud Top Level | hPa | single | up | 223 | 0 | 0 | 0',append)
rc=write(fname,'cltbot | Cloud base level | hPa | single | up | 2 | 0 | 0 | 0',append)
rc=write(fname,'cltlay | Total Cloud layer | | single | up | 3 | 0 | 0 | 0',append)
rc=write(fname,'cltmax | Highest Cloud height | m | single | up | 105 | 0 | 0 | 0',append)
rc=write(fname,'landd | Below ground, 10 to 200 cm| | single | up | 112 |10 |200 | 0',append)
rc=write(fname,'lands | Below ground, 0 to 10 cm | | single | up | 112 | 0 | 10 | 0',append)
rc=write(fname,'landt | Below ground, 0 to 200 cm| | single | up | 112 | 0 |200 | 0',append)
rc=write(fname,'lcl | Adiabatic cond level | hPa | single | up | 5 | 0 | 0 | 0',append)
rc=write(fname,'maxwnd | Maximum wind speed level | hPa | single | up | 6 | 0 | 0 | 0',append)
rc=write(fname,'msl | Mean Sea Level | | single | up | 102 | 0 | 0 | 0',append)
rc=write(fname,'ocnbot | Ocean bottom | | single | up | 9 | 0 | 0 | 0',append)
rc=write(fname,'plev | Pressure level | hPa | multi | down | 100 | 0 | 0 | 0',append)
rc=write(fname,'pbltop | Top of PBL | | single | up | 21 | 0 | 0 | 0',append)
rc=write(fname,'sfc | Earth surface | | single | up | 1 | 0 | 0 | 0',append)
rc=write(fname,'sfclo | Sfc Layer Ocean | | single | up | 112 | 0 |300 | 0',append)
rc=write(fname,'sfc10m | 10 meters above earth surface| m | single | up | 105 | 0 | 0 | 10',append)
rc=write(fname,'sfc2m | 2 meters above earth surface| m | single | up | 105 | 0 | 0 | 2',append)
rc=write(fname,'toa | Top of atmosphere | | single | up | 8 | 0 | 0 | 0',append)
rc=write(fname,'modtop | Top of Model | | single | up | 20 | 0 | 0 | 0',append)
rc=write(fname,'toasat | TOA satellite | | single | up | 22 | 0 | 0 | 0',append)
rc=write(fname,'troplev | Tropopause level | hPa | single | up | 7 | 0 | 0 | 0',append)
rc=write(fname,'theta | Isentropic Level | K | multi | up | 113 | 0 | 0 | 0',append)
rc=write(fname,'sigma | Sigma level | | multi | down | 107 | 0 | 0 | 0',append)
rc=write(fname,'hybrid | Hybrid Model level number | | multi | up | 109 | 0 | 0 | 0',append)
rc=write(fname,'zocean | Depth below sea level | m | multi | down | 160 | 0 | 0 | 0',append)
rc=write(fname,'',append)
rc=write(fname,'#---------------------------------------------------------------------------------------------------',append)
rc=write(fname,'#! Center',append)
rc=write(fname,'# Modeling centers (GRIB only)',append)
rc=write(fname,'# (max number of entries = LATS_MAX_CENTERS in lats.h)',append)
rc=write(fname,'#',append)
rc=write(fname,'# The format of each record is:',append)
rc=write(fname,'# center | GRIB_id | GRIB_center | GRIB_subcenter',append)
rc=write(fname,'#',append)
rc=write(fname,'# center = mnemonic for the center',append)
rc=write(fname,'# GRIB_id = GRIB generating process id (PDS octet 6)',append)
rc=write(fname,'# GRIB_center = the id of center managing the data (for AMIP II this is PCMDI) - see GRIB Table 0',append)
rc=write(fname,'# GRIB_subcenter = the id of the subcenter',append)
rc=write(fname,'# ',append)
rc=write(fname,'#',append)
rc=write(fname,'# Acronym AMIP Group Location',append)
rc=write(fname,'# ------- ---------- --------',append)
rc=write(fname,'#',append)
rc=write(fname,'# bmrc Bureau of Meteorology Research Centre Melbourne, Australia',append)
rc=write(fname,'# ccc Canadian Centre for Climate Modelling and Analysis Victoria, Canada',append)
rc=write(fname,'# ccsr Center for Climate System Research Tokyo, Japan',append)
rc=write(fname,'# cnrm Centre National de Recherches Meteorologiques Toulouse, France',append)
rc=write(fname,'# cola Center for Ocean-Land-Atmosphere Studies Calverton, Maryland',append)
rc=write(fname,'# csiro Commonwealth Scientific & Industrial Research Organization Mordialloc, Australia',append)
rc=write(fname,'# csu Colorado State University Fort Collins, Colorado',append)
rc=write(fname,'# derf Dynamical Extended Range Forecasting (at GFDL) Princeton, New Jersey',append)
rc=write(fname,'# dnm Department of Numerical Mathematics Moscow, Russia',append)
rc=write(fname,'# ecmwf European Centre for Medium-Range Weather Forecasts Reading, England',append)
rc=write(fname,'# gfdl Geophysical Fluid Dynamics Laboratory Princeton, New Jersey',append)
rc=write(fname,'# giss Goddard Institute for Space Studies New York, New York',append)
rc=write(fname,'# gla Goddard Laboratory for Atmospheres Greenbelt, Maryland',append)
rc=write(fname,'# gsfc Goddard Space Flight Center Greenbelt, Maryland',append)
rc=write(fname,'# iap Institute of Atmospheric Physics Beijing, China',append)
rc=write(fname,'# jma Japan Meteorological Agency Tokyo, Japan',append)
rc=write(fname,'# llnl Lawrence Livermore National Laboratory Livermore, California',append)
rc=write(fname,'# lmd Laboratoire de Meteorologie Dynamique Paris, France',append)
rc=write(fname,'# mgo Main Geophysical Observatory St. Petersburg, Russia',append)
rc=write(fname,'# mpi Max-Planck-Institut fur Meteorologie Hamburg, Germany',append)
rc=write(fname,'# mri Meteorological Research Institute Ibaraki-ken, Japan',append)
rc=write(fname,'# ncar National Center for Atmospheric Research Boulder, Colorado',append)
rc=write(fname,'# nmc National Meteorological Center Suitland, Maryland',append)
rc=write(fname,'# nrl Naval Research Laboratory Monterey, California',append)
rc=write(fname,'# ntu National Taiwan University Taipei, Taiwan',append)
rc=write(fname,'# pcmdi Program for Climate Model Diagnosis and Intercomparison, LLNL Livermore, California',append)
rc=write(fname,'# rpn Recherche en Privision Numerique Dorval, Canada',append)
rc=write(fname,'# sunya State University of New York at Albany Albany, New York',append)
rc=write(fname,'# sunya/ncar State University of New York at Albany/NCAR Albany, New York/Boulder, Colorado',append)
rc=write(fname,'# ucla University of California at Los Angeles Los Angeles, California',append)
rc=write(fname,'# ugamp The UK Universities Global Atmospheric Modelling Programme Reading, England',append)
rc=write(fname,'# uiuc University of Illinois at Urbana-Champaign Urbana, Illinois',append)
rc=write(fname,'# ukmo United Kingdom Meteorological Office Bracknell, UK',append)
rc=write(fname,'# yonu Yonsei University Seoul, Korea',append)
rc=write(fname,'#',append)
rc=write(fname,'',append)
rc=write(fname,'bmrc | 1 | 100 | 2',append)
rc=write(fname,'ccc | 2 | 100 | 2',append)
rc=write(fname,'cnrm | 3 | 100 | 2',append)
rc=write(fname,'cola | 4 | 100 | 2',append)
rc=write(fname,'csiro | 5 | 100 | 2',append)
rc=write(fname,'csu | 6 | 100 | 2',append)
rc=write(fname,'dnm | 7 | 100 | 2',append)
rc=write(fname,'ecmwf | 8 | 100 | 2',append)
rc=write(fname,'gfdl | 9 | 100 | 2',append)
rc=write(fname,'derf | 10 | 100 | 2',append)
rc=write(fname,'giss | 11 | 100 | 2',append)
rc=write(fname,'gla | 12 | 100 | 2',append)
rc=write(fname,'gsfc | 13 | 100 | 2',append)
rc=write(fname,'iap | 14 | 100 | 2',append)
rc=write(fname,'jma | 15 | 100 | 2',append)
rc=write(fname,'lmd | 16 | 100 | 2',append)
rc=write(fname,'mgo | 17 | 100 | 2',append)
rc=write(fname,'mpi | 18 | 100 | 2',append)
rc=write(fname,'mri | 19 | 100 | 2',append)
rc=write(fname,'ncar | 20 | 100 | 2',append)
rc=write(fname,'ncep | 21 | 100 | 2',append)
rc=write(fname,'nrl | 22 | 100 | 2',append)
rc=write(fname,'rpn | 23 | 100 | 2',append)
rc=write(fname,'sunya | 24 | 100 | 2',append)
rc=write(fname,'sunya/ncar| 25 | 100 | 2',append)
rc=write(fname,'ucla | 26 | 100 | 2',append)
rc=write(fname,'ugamp | 27 | 100 | 2',append)
rc=write(fname,'uiuc | 28 | 100 | 2',append)
rc=write(fname,'ukmo | 29 | 100 | 2',append)
rc=write(fname,'yonu | 30 | 100 | 2',append)
rc=write(fname,'ccsr | 31 | 100 | 2',append)
rc=write(fname,'llnl | 32 | 100 | 2',append)
rc=write(fname,'ntu | 33 | 100 | 2',append)
rc=write(fname,'cptec | 46 | 100 | 2',append)
rc=write(fname,'pcmdi | 100 | 100 | 2',append)
rc=write(fname,'#---------------------------------------------------------------------------------------------------',append)
rc=write(fname,'#!qc',append)
rc=write(fname,'# Quality control marks',append)
rc=write(fname,'# (no limit on number of entries)',append)
rc=write(fname,'#',append)
rc=write(fname,'# The format of each record is:',append)
rc=write(fname,'# variable | level_type | level | mean | std | tolerance | range | rangetol',append)
rc=write(fname,'#',append)
rc=write(fname,'# variable = variable name',append)
rc=write(fname,'# level_type = type of level, as defined in the leveltypes section, or blank if no associated level',append)
rc=write(fname,'# level = level value, or blank if no specified level',append)
rc=write(fname,'# mean = observed mean at specified level',append)
rc=write(fname,'# std = observed standard deviation at specified level',append)
rc=write(fname,'# tolerance = number of standard deviations about mean',append)
rc=write(fname,'# (if abs(calculated_mean - mean) > (std * tolerance), flag the value as 'mean out of range')',append)
rc=write(fname,'# range = observed (maximum - minimum)',append)
rc=write(fname,'# rangetol = range tolerance:',append)
rc=write(fname,'# (if calculated_range > (rangetol * range), flag as 'range is too large')',append)
rc=write(fname,'',append)
rc=write(fname,'# NOTE: no QC table yet',append)
rc=write(fname,'#',append)
rc=subwrd(rc,1)
if(rc!=0)
say _myname 'problems creating LATS table file ' fname
return rc
endif
rc=close(fname)
rc=subwrd(rc,1)
if(rc>0)
say _myname 'problems closing LATS table file ' fname
return rc
endif
return 0
*
* mkgeosf() Creates a template GEOS ctl for regridding;
* also sets _dimenv, _func
*
function mkgeosf(fname,dlon,dlat,nlon,nlat,algo)
*
* Temporary file
*
if ( _verb = 1 )
say _myname 'creating GEOS template file ' fname
endif
rc=write(fname,'DSET nofile')
rc=subwrd(rc,1)
if ( rc > 0 )
say rc
say _myname 'cannot create GEOS template file ' fname ' on current directory'
return rc
endif
xdef = 'xdef ' nlon ' linear -180 ' dlon
ydef = 'ydef ' nlat ' linear -90 ' dlat
rc=write(fname,'title Template GrADS for regridding',append)
rc=write(fname,'options template',append)
rc=write(fname,'undef 1e+20',append)
rc=write(fname,xdef,append)
rc=write(fname,ydef,append)
rc=write(fname,'zdef 1 levels 1000 ',append)
rc=write(fname,'tdef 1 linear 0Z1jan1900 1dy',append)
rc=write(fname,'vars 1',append)
rc=write(fname,'var 0 0 generic sfc variable',append)
rc=write(fname,'endvars',append)
rc=close(fname)
rc=subwrd(rc,1)
if(rc>0)
say _myname 'problems closing GEOS template file ' fname
return rc
endif
_dimenv = fname
_func = 'regrid2(@,' dlon ',' dlat ',' algo ',-180,-90)'
return 0
*
* validate(var,vars) See whether token var is in the list vars
*
function validate(var,vars)
if ( vars = '' ); return 1; endif
valid = 0
i = 1
v = subwrd(vars,i)
while ( v != '' )
if ( var = v ); valid = 1; endif
i = i + 1
v = subwrd(vars,i)
endwhile
return valid
*
* xvalidate(var,vars) See whether token var is NOT in the list vars
*
function xvalidat(var,vars)
if ( vars = '' ); return 1; endif
valid = 1
i = 1
v = subwrd(vars,i)
while ( v != '' )
if ( var = v ); valid = 0; endif
i = i + 1
v = subwrd(vars,i)
endwhile
return valid
*
* defvars() Define variables in output LATS file
*
function defvars()
* BUG: for now, cannot handle "accum" variables, oh well...
* ---------------------------------------------------------
if ( _mean = 1 )
timestat = 'average'
else
timestat = 'instant'
endif
* Surface variables
* ------------------
n = 1; _svids = ''
while ( n <= _nsvar )
var = subwrd(_svars,n)
if ( _format != stream )
'set lats var ' _fid ' ' var ' ' timestat ' ' _gid ' ' _sid
endif
vid = subwrd(result,5)
if ( vid = 0 ); return 1; endif
_svids = _svids ' ' vid
n = n + 1
endwhile
* Upper air variables
* -------------------
n = 1; _uvids = ''
while ( n <= _nuvar )
var = subwrd(_uvars,n)
if ( _format != stream )
'set lats var ' _fid ' ' var ' ' timestat ' ' _gid ' ' _lid
endif
vid = subwrd(result,5)
_uvids = _uvids ' ' vid
n = n + 1
endwhile
return 0
*............................................................................
* subst(var,expr) Substitute all the occurences of '@' in expr with the
* contents of var. Example:
*
* expr = cos(lat*pi/180)*@*hgrad(@)
* var = ua
*
* results in cos(lat*pi/180)*ua*hgrad(ua).
*
* If expr='' it returns var.
*
function subst(var,expr)
if ( expr ='' ); return var; endif
str = ''
i = 1
ch = substr(expr,i,1)
while( ch != '' )
if ( ch = '@' )
str = str '' var
else
str = str '' ch
endif
i = i + 1
ch = substr(expr,i,1)
endwhile
return str
*............................................................................
*
* chkcfg() Check whether grads version is powerful enough to run lats4d
*
function chkcfg()
'q config'
if ( rc != 0 )
say _myname 'GrADS version appears earlier than 1.7Beta'
return 1
endif
* Check for LATS
* --------------
lats = 0
cfg = sublin(result,1)
i = 1
version = subwrd(cfg,2)
token = subwrd(cfg,i)
while ( token != '' )
if ( token = 'lats' ); lats = 1; endif
i = i + 1
token = subwrd(cfg,i)
endwhile
if ( lats = 0 )
say _myname 'this build of GrADS ' version ' does not include the LATS option'
say _myname 'please try "gradsnc" or "gradshdf", version 1.7beta9 or later'
return 1
endif
return 0
*............................................................................
*
* crtique(nmin,nmax,opt,string) Returns 0 if the number of words in "string"
* is in between nmin and nmax, returns 1 otherwise.
*
function critique ( nmin, nmax, opt, string)
n = 1
while ( subwrd(string,n) != '' ); n = n + 1; endwhile
n = n - 1
if ( n < nmin )
rc = usage()
say _myname 'you specified "' opt ' ' string '"'
say _myname 'option "' opt '" requires at least ' nmin ' argument(s)'
return 1
endif
if ( n > nmax )
rc = usage()
say _myname 'you specified "' opt ' ' string '"'
say _myname 'option "' opt '" requires at most ' nmax ' argument(s)'
return 1
endif
return 0
More information about the gradsusr
mailing list