[
Date Prev][
Date Next][
Thread Prev][
Thread Next][
Date Index][
Thread Index]
Re: 20020410: gribtonc CDLs for AVN and could CDL for LAPS work
- Subject: Re: 20020410: gribtonc CDLs for AVN and could CDL for LAPS work
- Date: Tue, 23 Apr 2002 18:32:47 -0600 (MDT)
Tiziana,
I got a proto version working but there were bugs in gribtonc. I,m
attaching the 3 files plus the the cdl. Insert the code into the dir and
do a make then try it out. I started with a another cdl so could you
check for vars that are not used and then comment them out with '//'. Let
me know what you find out.
Thanks,
Robb...
===============================================================================
Robb Kambic Unidata Program Center
Software Engineer III Univ. Corp for Atmospheric Research
address@hidden WWW: http://www.unidata.ucar.edu/
===============================================================================
netcdf avn-x { // 126 Wave, 18 Layer Spectral Model Aviation Run
// on expanded quasi-regular "thinned" grids
dimensions:
record = UNLIMITED ; // (reference time, forecast time)
level = 26 ;
lat = 181 ; // latitude
lon = 360 ; // longitude
lpdg = 2 ; // boundary layer levels
fhg = 2 ; // fixed height above ground levels
fh = 3 ; // fixed height above ground levels
soil_lpdg = 6 ; // soil boundary layer levels
ls_all = 7 ; // whole atmosphere layer
sigma = 1 ; // sigma level
datetime_len = 21 ; // string length for datetime strings
nmodels = 3 ; // both AVN and SSIAVN models
accum = 2 ; // time range for accumulations
nav = 1 ; // For navigation. Variables that use
// this dimension define a mapping between
// (lat, lon) indices and (lat, lon) coords.
nav_len = 100 ; // max string length for navigation strings
variables:
double reftime(record) ; // reference time of the model
reftime:long_name = "reference time" ;
reftime:units = "hours since 1992-1-1" ;
double valtime(record) ; // forecast time ("valid" time)
valtime:long_name = "valid time" ;
valtime:units = "hours since 1992-1-1" ;
:record = "reftime, valtime" ; // "dimension attribute" -- means
// (reftime, valtime) uniquely
// determine record
char datetime(record, datetime_len) ; // derived from reftime
datetime:long_name = "reference date and time" ;
// units YYYY-MM-DD hh:mm:ssZ (ISO 8601)
float valtime_offset(record) ; // derived as valtime-reftime
valtime_offset:long_name = "hours from reference time" ;
valtime_offset:units = "hours" ;
float level(level) ;
level:long_name = "level" ;
level:units = "hectopascals" ;
// (soil_lpdg_bot, soil_lpdg_top) uniquely determines soil_lpdg
float soil_lpdg_bot(soil_lpdg) ;
soil_lpdg_bot:long_name = "bottom level of boundary layer
between 2 levels from ground to levels" ;
//soil_lpdg_bot:units = "centimeters" ;
float soil_lpdg_top(soil_lpdg) ;
soil_lpdg_top:long_name = "top level of boundary layer between
2 levels from ground to levels" ;
//soil_lpdg_top:units = "centimeters" ;
:lpdg = "lpdg_bot, lpdg_top" ; // (lpdg_bot, lpdg_top) uniquely
// determines lpdg
float lpdg_bot(lpdg) ;
lpdg_bot:long_name = "bottom level of boundary layer between 2
levels at specified pressure differences from ground to levels" ;
lpdg_bot:units = "hectopascals" ;
float lpdg_top(lpdg) ;
lpdg_top:long_name = "top level of boundary layer between 2
levels at specified pressure differences from ground to levels" ;
lpdg_top:units = "hectopascals" ;
float fhg(fhg) ; // fixed height above ground
fhg:long_name = "fixed height above ground" ;
fhg:units = "meters" ;
float fh(fh) ; // fixed height above ground
fh:long_name = "fixed height above ground" ;
fh:units = "meters" ;
:ls_all = "ls_all_bot, ls_all_top" ;
float ls_all_bot(ls_all) ;
ls_all_bot:long_name = "bottom level of atmosphere between 2
sigma levels" ;
ls_all_bot:units = "" ;
float ls_all_top(ls_all) ;
ls_all_top:long_name = "top level of atmosphere between 2 sigma
levels" ;
ls_all_top:units = "" ;
float sigma(sigma) ; // fixed height above ground
sigma:long_name = "sigma level" ;
sigma:units = "" ; // dimensionless
long model_id(nmodels) ;
model_id:long_name = "generating process ID number" ;
// The following lat and lon coordinate variables are redundant,
// since the navigation variables provide the necessary information.
// The extra information is included here for human readability.
float lat(lat) ;
lat:long_name = "latitude" ;
lat:units = "degrees_north" ;
float lon(lon) ;
lon:long_name = "longitude" ;
lon:units = "degrees_east" ;
// navigation variables all use nav dimension
char nav_model(nav, nav_len) ; // navigation parameterization
nav_model:long_name = "navigation model name" ;
int grid_type_code(nav) ;
grid_type_code:long_name = "GRIB-1 GDS data representation
type" ;
char grid_type(nav, nav_len) ;
grid_type:long_name = "GRIB-1 grid type" ;
char grid_name(nav, nav_len) ;
grid_name:long_name = "grid name" ;
int grid_center(nav) ;
grid_center:long_name = "GRIB-1 originating center ID" ;
int grid_number(nav) ;
grid_number:long_name = "GRIB-1 catalogued grid numbers" ;
grid_number:_FillValue = -9999 ;
char i_dim(nav, nav_len) ;
i_dim:long_name = "longitude dimension name" ;
char j_dim(nav, nav_len) ;
j_dim:long_name = "latitude dimension name" ;
int Ni(nav) ;
Ni:long_name = "number of points along a latitude circle" ;
int Nj(nav) ;
Nj:long_name = "number of points along a longitude circle" ;
float La1(nav) ;
La1:long_name = "latitude of first grid point" ;
La1:units = "degrees_north" ;
float Lo1(nav) ;
Lo1:long_name = "longitude of first grid point" ;
Lo1:units = "degrees_east" ;
float La2(nav) ;
La2:long_name = "latitude of last grid point" ;
La2:units = "degrees_north" ;
float Lo2(nav) ;
Lo2:long_name = "longitude of last grid point" ;
Lo2:units = "degrees_east" ;
float Di(nav) ;
Di:long_name = "Longitudinal direction increment" ;
Di:units = "degrees" ;
float Dj(nav) ;
Dj:long_name = "Latitudinal direction increment" ;
Dj:units = "degrees" ;
byte ResCompFlag(nav) ;
ResCompFlag:long_name = "resolution and component flags" ;
// end of navigation variables
float P(record, lat, lon) ;
P:long_name = "pressure" ;
P:units = "Pa" ;
P:_FillValue = -9999.f ;
P:navigation = "nav" ;
float P_maxwind(record, lat, lon) ;
P_maxwind:long_name = "pressure at maximum wind" ;
P_maxwind:units = "Pa" ;
P_maxwind:_FillValue = -9999.f ;
P_maxwind:navigation = "nav" ;
float P_trop(record, lat, lon) ;
P_trop:long_name = "pressure at tropopause" ;
P_trop:units = "Pa" ;
P_trop:_FillValue = -9999.f ;
P_trop:navigation = "nav" ;
float P_msl(record, lat, lon) ;
P_msl:long_name = "pressure reduced to MSL" ;
P_msl:units = "Pa" ;
P_msl:_FillValue = -9999.f ;
P_msl:navigation = "nav" ;
float P_sfc(record, lat, lon) ;
P_sfc:long_name = "pressure at surface" ;
P_sfc:units = "Pa" ;
P_sfc:_FillValue = -9999.f ;
P_sfc:navigation = "nav" ;
float RH(record, level, lat, lon) ;
RH:long_name = "relative humidity" ;
RH:units = "percent" ;
RH:_FillValue = -9999.f ;
RH:navigation = "nav" ;
float T(record, level, lat, lon) ;
T:long_name = "temperature" ;
T:units = "degK" ;
T:_FillValue = -9999.f ;
T:navigation = "nav" ;
float T_maxwind(record, lat, lon) ;
T_maxwind:long_name = "temperature at maxwind" ;
T_maxwind:units = "degK" ;
T_maxwind:_FillValue = -9999.f ;
T_maxwind:navigation = "nav" ;
float T_trop(record, lat, lon) ;
T_trop:long_name = "temperature at tropopause" ;
T_trop:units = "degK" ;
T_trop:_FillValue = -9999.f ;
T_trop:navigation = "nav" ;
float T_sigma(record, sigma, lat, lon) ;
T_sigma:long_name = "temperature" ;
T_sigma:units = "degK" ;
T_sigma:_FillValue = -9999.f ;
T_sigma:navigation = "nav" ;
float RH_sigma(record, sigma, lat, lon) ;
RH_sigma:long_name = "relative humidity at sigma level" ;
RH_sigma:units = "percent" ;
RH_sigma:_FillValue = -9999.f ;
RH_sigma:navigation = "nav" ;
float u_sigma(record, sigma, lat, lon ) ;
u_sigma:long_name = "u-component of wind" ;
u_sigma:units = "m/s" ;
u_sigma:_FillValue = -9999.f ;
u_sigma:navigation = "nav" ;
float v_sigma(record, sigma, lat, lon ) ;
v_sigma:long_name = "v-component of wind" ;
v_sigma:units = "m/s" ;
v_sigma:_FillValue = -9999.f ;
v_sigma:navigation = "nav" ;
float theta_sigma(record, sigma, lat, lon) ;
theta_sigma:long_name = "Potential temperature" ;
theta_sigma:units = "degK" ;
theta_sigma:_FillValue = -9999.f ;
theta_sigma:navigation = "nav" ;
float omega_sigma(record, sigma, lat, lon) ;
omega_sigma:long_name = "pressure vertical velocity" ;
omega_sigma:units = "Pa/s" ;
omega_sigma:_FillValue = -9999.f ;
omega_sigma:navigation = "nav" ; // georeference info
float Z(record, level, lat, lon) ;
Z:long_name = "geopotential height" ;
Z:units = "gp m" ;
Z:_FillValue = -9999.f ;
Z:navigation = "nav" ; // georeference info
float Z_maxwind(record, lat, lon) ;
Z_maxwind:long_name = "geopotential height at maxwind" ;
Z_maxwind:units = "gp m" ;
Z_maxwind:_FillValue = -9999.f ;
Z_maxwind:navigation = "nav" ;
float Z_trop(record, lat, lon) ;
Z_trop:long_name = "geopotential height at tropopause" ;
Z_trop:units = "gp m" ;
Z_trop:_FillValue = -9999.f ;
Z_trop:navigation = "nav" ;
float T_lpdg(record, lpdg, lat, lon) ;
T_lpdg:long_name = "temperature in boundary layer" ;
T_lpdg:units = "degK" ;
T_lpdg:_FillValue = -9999.f ;
T_lpdg:navigation = "nav" ;
float RH_lpdg(record, lpdg, lat, lon) ;
RH_lpdg:long_name = "relative humidity in boundary layer" ;
RH_lpdg:units = "percent" ;
RH_lpdg:_FillValue = -9999.f ;
RH_lpdg:navigation = "nav" ;
float cin_lpdg(record, lpdg, lat, lon ) ;
cin_lpdg:long_name = "boundary convective inhibition" ;
cin_lpdg:units = "J/kg" ;
cin_lpdg:_FillValue = -9999.f ;
cin_lpdg:navigation = "nav" ;
float spec_hum_lpdg(record, lpdg, lat, lon ) ;
spec_hum_lpdg:long_name = "specific humidity" ;
spec_hum_lpdg:units = "kg/kg" ;
spec_hum_lpdg:_FillValue = -9999.f ;
spec_hum_lpdg:navigation = "nav" ;
float cape_lpdg(record, lpdg, lat, lon ) ;
cape_lpdg:long_name = "boundary convective available potential
energy" ;
cape_lpdg:units = "J/kg" ;
cape_lpdg:_FillValue = -9999.f ;
cape_lpdg:navigation = "nav" ;
float u_lpdg(record, lpdg, lat, lon) ;
u_lpdg:long_name = "u-component of wind in boundary layer" ;
u_lpdg:units = "meters/second" ;
u_lpdg:_FillValue = -9999.f ;
u_lpdg:navigation = "nav" ;
float v_lpdg(record, lpdg, lat, lon) ;
v_lpdg:long_name = "v-component of wind in boundary layer" ;
v_lpdg:units = "meters/second" ;
v_lpdg:_FillValue = -9999.f ;
v_lpdg:navigation = "nav" ;
float u(record, level, lat, lon) ;
u:long_name = "u-component of wind" ;
u:units = "meters/second" ;
u:_FillValue = -9999.f ;
u:navigation = "nav" ;
float u_maxwind(record, lat, lon) ;
u_maxwind:long_name = "u-component of wind at max wind" ;
u_maxwind:units = "meters/second" ;
u_maxwind:_FillValue = -9999.f ;
u_maxwind:navigation = "nav" ;
float u_trop(record, lat, lon) ;
u_trop:long_name = "u-component of wind at tropopause" ;
u_trop:units = "meters/second" ;
u_trop:_FillValue = -9999.f ;
u_trop:navigation = "nav" ;
float v(record, level, lat, lon) ;
v:long_name = "v-component of wind" ;
v:units = "meters/second" ;
v:_FillValue = -9999.f ;
v:navigation = "nav" ;
float v_maxwind(record, lat, lon) ;
v_maxwind:long_name = "v-component of wind at max wind" ;
v_maxwind:units = "meters/second" ;
v_maxwind:_FillValue = -9999.f ;
v_maxwind:navigation = "nav" ;
float v_trop(record, lat, lon) ;
v_trop:long_name = "v-component of wind at tropopause" ;
v_trop:units = "meters/second" ;
v_trop:_FillValue = -9999.f ;
v_trop:navigation = "nav" ;
float u_fhg(record, fhg, lat, lon) ;
u_fhg:long_name = "u-component of wind at fixed height above
ground" ;
u_fhg:units = "meters/second" ;
u_fhg:_FillValue = -9999.f ;
u_fhg:navigation = "nav" ;
float v_fhg(record, fhg, lat, lon) ;
v_fhg:long_name = "v-component of wind at fixed height above
ground" ;
v_fhg:units = "meters/second" ;
v_fhg:_FillValue = -9999.f ;
v_fhg:navigation = "nav" ;
float RH_fhg(record, fhg, lat, lon) ;
RH_fhg:long_name = "relative humidity at fixed height above
ground" ;
RH_fhg:units = "percent" ;
RH_fhg:_FillValue = -9999.f ;
RH_fhg:navigation = "nav" ;
float spec_hum_fhg(record, fhg, lat, lon) ;
spec_hum_fhg:long_name = "specific humidity at fixed height
above ground" ;
spec_hum_fhg:units = "kg/kg" ;
spec_hum_fhg:_FillValue = -9999.f ;
spec_hum_fhg:navigation = "nav" ;
float T_fhg(record, fhg, lat, lon) ;
T_fhg:long_name = "temperature at fixed height above ground" ;
T_fhg:units = "degK" ;
T_fhg:_FillValue = -9999.f ;
T_fhg:navigation = "nav" ;
float Tmax_fhg(record, fhg, lat, lon) ;
Tmax_fhg:long_name = " maximum temperature at fixed height
above ground" ;
Tmax_fhg:units = "degK" ;
Tmax_fhg:_FillValue = -9999.f ;
Tmax_fhg:navigation = "nav" ;
float Tmin_fhg(record, fhg, lat, lon) ;
Tmin_fhg:long_name = " minimum temperature at fixed height
above ground" ;
Tmin_fhg:units = "degK" ;
Tmin_fhg:_FillValue = -9999.f ;
Tmin_fhg:navigation = "nav" ;
float u_fh(record, fh, lat, lon) ;
u_fh:long_name = "u-component of wind at fixed height above
ground" ;
u_fh:units = "meters/second" ;
u_fh:_FillValue = -9999.f ;
u_fh:navigation = "nav" ;
float v_fh(record, fh, lat, lon) ;
v_fh:long_name = "v-component of wind at fixed height above
ground" ;
v_fh:units = "meters/second" ;
v_fh:_FillValue = -9999.f ;
v_fh:navigation = "nav" ;
float RH_fh(record, fh, lat, lon) ;
RH_fh:long_name = "relative humidity at fixed height above
ground" ;
RH_fh:units = "percent" ;
RH_fh:_FillValue = -9999.f ;
RH_fh:navigation = "nav" ;
float T_fh(record, fh, lat, lon) ;
T_fh:long_name = "temperature at fixed height above ground" ;
T_fh:units = "degK" ;
T_fh:_FillValue = -9999.f ;
T_fh:navigation = "nav" ;
float RH_ls(record, ls_all, lat, lon) ;
RH_ls:long_name = "relative humidity" ;
RH_ls:units = "percent" ;
RH_ls:_FillValue = -9999.f ;
RH_ls:navigation = "nav" ;
float PRECIP(record, lat, lon) ;
PRECIP:long_name = "total precipitation over accumulation
interval" ;
PRECIP:units = "kg/m2" ;
PRECIP:_FillValue = -9999.f ;
PRECIP:navigation = "nav" ;
float PRECIP_accum_times(record, accum) ;
PRECIP_accum_times:long_name = "precipitation accumulation
interval" ;
PRECIP_accum_times:units = "hours" ;
PRECIP_accum_times:_FillValue = -9999.f ;
float precip_cn(record, lat, lon) ;
precip_cn:long_name = "convective precipitation over
accumulation interval" ;
precip_cn:units = "kg/m2" ;
precip_cn:_FillValue = -9999.f ;
precip_cn:navigation = "nav" ;
float precip_cn_accum_times(record, accum) ;
precip_cn_accum_times:long_name = "convective precipitation
accumulation interval" ;
precip_cn_accum_times:units = "hours" ;
precip_cn_accum_times:_FillValue = -9999.f ;
float precip_rt(record, lat, lon ) ;
precip_rt:long_name = "precipitation rate" ;
precip_rt:units = "kg/m2/s" ;
precip_rt:_FillValue = -9999.f ;
precip_rt:navigation = "nav" ;
float watr(record, lat, lon ) ;
watr:long_name = "water runoff" ;
watr:units = "kg/m2" ;
watr:_FillValue = -9999.f ;
watr:navigation = "nav" ;
float pr_water_atm(record, lat, lon ) ; // entire atmosphere as
single layer
pr_water_atm:long_name = "precipitable water" ;
pr_water_atm:units = "kg/m2" ;
pr_water_atm:_FillValue = -9999.f ;
pr_water_atm:navigation = "nav" ;
float cprat(record, lat, lon) ;
cprat:long_name = "Convective precipitation rate" ;
cprat:units = "kg/m2/sec" ;
cprat:_FillValue = -9999.f ;
cprat:navigation = "nav" ;
float crain(record, lat, lon ) ;
crain:long_name = "Categorical rain" ;
crain:_FillValue = -9999.f ;
crain:navigation = "nav" ;
float cfrzrn(record, lat, lon ) ;
cfrzrn:long_name = "Categorical freezing rain" ;
cfrzrn:_FillValue = -9999.f ;
cfrzrn:navigation = "nav" ;
float csnow(record, lat, lon ) ;
csnow:long_name = "Categorical snow" ;
csnow:_FillValue = -9999.f ;
csnow:navigation = "nav" ;
float ice_conc(record, lat, lon ) ;
ice_conc:long_name = "Ice concentration" ;
ice_conc:_FillValue = -9999.f ;
ice_conc:navigation = "nav" ;
float LI(record, lat, lon ) ;
LI:long_name = "lifted index" ;
LI:units = "degK" ;
LI:_FillValue = -9999.f ; // To fill grid corners
LI:navigation = "nav" ;
float T_sfc(record, lat, lon) ;
T_sfc:long_name = "surface temperature" ;
T_sfc:units = "degK" ;
T_sfc:_FillValue = -9999.f ;
T_sfc:navigation = "nav" ;
float Z_sfc(record, lat, lon) ;
Z_sfc:long_name = "terrain" ;
Z_sfc:units = "gp m" ;
Z_sfc:_FillValue = -9999.f ;
Z_sfc:navigation = "nav" ;
float sen_ht_sfc(record, lat, lon ) ;
sen_ht_sfc:long_name = "Sensible heat net flux" ;
sen_ht_sfc:units = "W / m2" ;
sen_ht_sfc:_FillValue = -9999.f ;
sen_ht_sfc:navigation = "nav" ;
float cin_sfc(record, lat, lon ) ;
cin_sfc:long_name = "surface convective inhibition" ;
cin_sfc:units = "J/kg" ;
cin_sfc:_FillValue = -9999.f ;
cin_sfc:navigation = "nav" ;
float uswrf_sfc(record, lat, lon ) ;
uswrf_sfc:long_name = "Upward short wave rad. flux" ;
uswrf_sfc:units = "W/m2" ;
uswrf_sfc:_FillValue = -9999.f ;
uswrf_sfc:navigation = "nav" ;
float dswrf_sfc(record, lat, lon ) ;
dswrf_sfc:long_name = "Downward short wave rad. flux" ;
dswrf_sfc:units = "W/m2" ;
dswrf_sfc:_FillValue = -9999.f ;
dswrf_sfc:navigation = "nav" ;
float ulwrf_sfc(record, lat, lon ) ;
ulwrf_sfc:long_name = "Upward long wave rad. flux" ;
ulwrf_sfc:units = "W / m2" ;
ulwrf_sfc:_FillValue = -9999.f ;
ulwrf_sfc:navigation = "nav" ;
float dlwrf_sfc(record, lat, lon ) ;
dlwrf_sfc:long_name = "Downward long wave rad. flux" ;
dlwrf_sfc:units = "W / m2" ;
dlwrf_sfc:_FillValue = -9999.f ;
dlwrf_sfc:navigation = "nav" ;
float land_mask_sfc(record, lat, lon ) ;
land_mask_sfc:long_name = "Land-Sea mask" ;
land_mask_sfc:units = "bit" ;
land_mask_sfc:_FillValue = -9999.f ;
land_mask_sfc:navigation = "nav" ;
float albedo_sfc(record, lat, lon ) ;
albedo_sfc:long_name = "Albedo" ;
albedo_sfc:_FillValue = -9999.f ;
albedo_sfc:navigation = "nav" ;
// Latent heat net flux
float lat_ht_sfc(record, lat, lon ) ;
lat_ht_sfc:long_name = "Latent heat net flux" ;
lat_ht_sfc:units = "W / m2" ;
lat_ht_sfc:_FillValue = -9999.f ;
lat_ht_sfc:navigation = "nav" ;
float LI4_sfc(record, lat, lon ) ;
LI4_sfc:long_name = "Best 4 layer lift index" ;
LI4_sfc:units = "K" ;
LI4_sfc:_FillValue = -9999.f ;
LI4_sfc:navigation = "nav" ;
float cape_sfc(record, lat, lon ) ;
cape_sfc:long_name = "surface convective available potential
energy" ;
cape_sfc:units = "J/kg" ;
cape_sfc:_FillValue = -9999.f ;
cape_sfc:navigation = "nav" ;
float u_flx_sfc(record, lat, lon ) ;
u_flx_sfc:long_name = "Momentum flux, u componet" ;
u_flx_sfc:units = "N/m2" ;
u_flx_sfc:_FillValue = -9999.f ;
u_flx_sfc:navigation = "nav" ;
float v_flx_sfc(record, lat, lon ) ;
v_flx_sfc:long_name = "Momentum flux, v componet" ;
v_flx_sfc:units = "N/m2" ;
v_flx_sfc:_FillValue = -9999.f ;
v_flx_sfc:navigation = "nav" ;
// Planetary boundary layer height
float hpbl_sfc(record) ;
hpbl_sfc:long_name = "Planetary boundary layer height" ;
hpbl_sfc:units = "m" ;
hpbl_sfc:_FillValue = -9999.f ;
hpbl_sfc:navigation = "nav" ;
float RH_atm(record, lat, lon) ;
RH_atm:long_name = "relative humidity entire atmosphere" ;
RH_atm:units = "percent" ;
RH_atm:_FillValue = -9999.f ;
RH_atm:navigation = "nav" ;
float totoz_atm(record, lat, lon) ;
totoz_atm:long_name = "Total ozone entire atmosphere" ;
//totoz_atm:units = "Dobson" ;
totoz_atm:_FillValue = -9999.f ;
totoz_atm:navigation = "nav" ;
float RH_frzlvl(record, lat, lon ) ;
RH_frzlvl:long_name = "relative humidity at 0 degree isotherm" ;
RH_frzlvl:units = "percent" ;
RH_frzlvl:_FillValue = -9999.f ;
RH_frzlvl:navigation = "nav" ;
float Z_frzlvl(record, lat, lon ) ;
Z_frzlvl:long_name = "geopotential height at 0 isotherm" ;
Z_frzlvl:units = "gp m" ;
Z_frzlvl:_FillValue = -9999.f ;
Z_frzlvl:navigation = "nav" ;
float T_hctl(record, lat, lon) ;
T_hctl:long_name = "temperature at high cloud top level" ;
T_hctl:units = "degK" ;
T_hctl:_FillValue = -9999.f ;
T_hctl:navigation = "nav" ;
float P_hctl(record, lat, lon) ;
P_hctl:long_name = "Pressure at high cloud top level" ;
P_hctl:units = "Pa" ;
P_hctl:_FillValue = -9999.f ;
P_hctl:navigation = "nav" ;
float T_mctl(record, lat, lon) ;
T_mctl:long_name = "temperature at middle cloud top level" ;
T_mctl:units = "degK" ;
T_mctl:_FillValue = -9999.f ;
T_mctl:navigation = "nav" ;
float P_mctl(record, lat, lon) ;
P_mctl:long_name = "Pressure at middle cloud top level" ;
P_mctl:units = "Pa" ;
P_mctl:_FillValue = -9999.f ;
P_mctl:navigation = "nav" ;
float T_lctl(record, lat, lon) ;
T_lctl:long_name = "temperature at low cloud top level" ;
T_lctl:units = "degK" ;
T_lctl:_FillValue = -9999.f ;
T_lctl:navigation = "nav" ;
float P_lctl(record, lat, lon) ;
P_lctl:long_name = "Pressure at low cloud top level" ;
P_lctl:units = "Pa" ;
P_lctl:_FillValue = -9999.f ;
P_lctl:navigation = "nav" ;
float P_mcbl(record, lat, lon) ;
P_mcbl:long_name = "Pressure at middle cloud bottom level" ;
P_mcbl:units = "Pa" ;
P_mcbl:_FillValue = -9999.f ;
P_mcbl:navigation = "nav" ;
float P_cctl(record, lat, lon) ;
P_cctl:long_name = "Pressure at convective cloud top layer" ;
P_cctl:units = "Pa" ;
P_cctl:_FillValue = -9999.f ;
P_cctl:navigation = "nav" ;
float P_ccbl(record, lat, lon) ;
P_ccbl:long_name = "Pressure at convective cloud bottom layer" ;
P_ccbl:units = "Pa" ;
P_ccbl:_FillValue = -9999.f ;
P_ccbl:navigation = "nav" ;
float uswrf_topa(record, lat, lon) ;
uswrf_topa:long_name = "Upward short wave rad.flux" ;
uswrf_topa:units = "W / m2" ;
uswrf_topa:_FillValue = -9999.f ;
uswrf_topa:navigation = "nav" ;
float ulwrf_topa(record, lat, lon) ;
ulwrf_topa:long_name = "Upward long wave rad.flux" ;
ulwrf_topa:units = "W / m2" ;
ulwrf_topa:_FillValue = -9999.f ;
ulwrf_topa:navigation = "nav" ;
float N(record, lat, lon) ;
N:long_name = "Total cloud cover" ;
N:units = "percent" ;
N:_FillValue = -9999.f ;
N:navigation = "nav" ;
float N_hcy(record, lat, lon) ;
N_hcy:long_name = "Total cloud cover, high cloud layer" ;
N_hcy:_FillValue = -9999.f ;
N_hcy:navigation = "nav" ;
float N_mcy(record, lat, lon) ;
N_mcy:long_name = "Total cloud cover, middle cloud layer" ;
N_mcy:_FillValue = -9999.f ;
N_mcy:navigation = "nav" ;
float N_lcy(record, lat, lon) ;
N_lcy:long_name = "Total cloud cover, low cloud layer" ;
N_lcy:_FillValue = -9999.f ;
N_lcy:navigation = "nav" ;
float N_bcy(record, lat, lon) ;
N_bcy:long_name = "Total cloud cover, boundary layer cloud
layer" ;
N_bcy:_FillValue = -9999.f ;
N_bcy:navigation = "nav" ;
float N_ccy(record, lat, lon) ;
N_ccy:long_name = "Total cloud cover, covective cloud layer" ;
N_ccy:_FillValue = -9999.f ;
N_ccy:navigation = "nav" ;
float N_atm(record, lat, lon) ;
N_atm:long_name = "Total cloud cover entire atmosphere" ;
N_atm:_FillValue = -9999.f ;
N_atm:navigation = "nav" ;
float P_hcbl(record, lat, lon) ;
P_hcbl:long_name = "Pressure, high cloud bottom level" ;
P_hcbl:units = "Pa" ;
P_hcbl:_FillValue = -9999.f ;
P_hcbl:navigation = "nav" ;
float P_lcbl(record, lat, lon) ;
P_lcbl:long_name = "Pressure, low cloud bottom level" ;
P_lcbl:units = "Pa" ;
P_lcbl:_FillValue = -9999.f ;
P_lcbl:navigation = "nav" ;
float omega(record, level, lat, lon) ;
omega:long_name = "pressure vertical velocity" ;
omega:units = "Pa/s" ;
omega:_FillValue = -9999.f ;
omega:navigation = "nav" ; // georeference info
float absvor(record, level, lat, lon) ;
absvor:long_name = "absolute vorticity" ;
absvor:units = "1/s" ;
absvor:_FillValue = -9999.f ;
absvor:navigation = "nav" ;
// Cloud water
float clwmr(record, level, lat, lon ) ;
clwmr:long_name = "Cloud water" ;
clwmr:units = "kg / kg" ;
clwmr:_FillValue = -9999.f ;
clwmr:navigation = "nav" ;
float cloud_wat_atm(record, lat, lon ) ;
cloud_wat_atm:long_name = "Cloud water" ;
cloud_wat_atm:units = "kg / m2" ;
cloud_wat_atm:_FillValue = -9999.f ;
cloud_wat_atm:navigation = "nav" ;
float snow_wat(record, lat, lon ) ;
snow_wat:long_name = "Water equiv. of accumulated snow depth" ;
snow_wat:units = "kg / m2" ;
snow_wat:_FillValue = -9999.f ;
snow_wat:navigation = "nav" ;
float cicepl(record, lat, lon ) ;
cicepl:long_name = "Categorical ice pellets" ;
cicepl:_FillValue = -9999.f ;
cicepl:navigation = "nav" ;
// Ozone mixing ratio
float o3mr(record, level, lat, lon ) ;
o3mr:long_name = "Ozone mixing ratio" ;
o3mr:units = "kg / kg" ;
o3mr:_FillValue = -9999.f ;
o3mr:navigation = "nav" ;
float Zdev(record, level, lat, lon ) ;
Zdev:long_name = "Geopotential height anomaly" ;
Zdev:units = "gp m" ;
Zdev:_FillValue = -9999.f ;
Zdev:navigation = "nav" ;
float gpt_hgt5(record, level, lat, lon ) ;
gpt_hgt5:long_name = "5-wave Geopotential height" ;
gpt_hgt5:units = "gp m" ;
gpt_hgt5:_FillValue = -9999.f ;
gpt_hgt5:navigation = "nav" ;
float gflux(record, lat, lon ) ;
gflux:long_name = "Ground heat flux" ;
gflux:units = "W / m2" ;
gflux:_FillValue = -9999.f ;
gflux:navigation = "nav" ;
float vert_sshr_trop(record, lat, lon ) ;
vert_sshr_trop:long_name = "vertical speed shear" ;
vert_sshr_trop:units = "1/s" ;
vert_sshr_trop:_FillValue = -9999.f ;
vert_sshr_trop:navigation = "nav" ;
float T_lbls(record, soil_lpdg, lat, lon ) ;
T_lbls:long_name = "Temperature layer between 2 depth below
surface" ;
T_lbls:units = "K" ;
T_lbls:_FillValue = -9999.f ;
T_lbls:navigation = "nav" ;
float Z_htfl(record, lat, lon ) ;
Z_htfl:long_name = "geopotential height" ;
Z_htfl:units = "gp m" ;
Z_htfl:_FillValue = -9999.f ;
Z_htfl:navigation = "nav" ;
float RH_htfl(record, lat, lon ) ;
RH_htfl:long_name = "relative humidity" ;
RH_htfl:units = "percent" ;
RH_htfl:_FillValue = -9999.f ;
RH_htfl:navigation = "nav" ;
float reserved(record, level, lat, lon ) ;
reserved:long_name = "" ;
reserved:_FillValue = -9999.f ;
reserved:navigation = "nav" ;
float reserved_lbls(record, soil_lpdg, lat, lon ) ;
reserved_lbls:long_name = "Volumetric soil moisture content" ;
reserved_lbls:_FillValue = -9999.f ;
reserved_lbls:navigation = "nav" ;
float reserved_atm(record, lat, lon ) ;
reserved_atm:long_name = "Cloud workfunction" ;
reserved_atm:_FillValue = -9999.f ;
reserved_atm:navigation = "nav" ;
float reserved_sfc(record, lat, lon ) ;
reserved_sfc:long_name = "Meridional flux of gravity wave
stress" ;
reserved_sfc:_FillValue = -9999.f ;
reserved_sfc:navigation = "nav" ;
// global attributes:
:history = "created by gribtonc from HRS broadcast" ;
:title = "NMC Global Product Set" ;
:Conventions = "NUWG" ;
:version = 0.0 ; // still just a draft
data:
level = 1000, 975, 950, 925, 900, 850, 800, 750, 700, 650, 600, 550, 500, 450,
400, 350, 300, 250, 200, 150, 100, 70, 50, 30, 20, 10 ;
soil_lpdg_top = 0, 0, 5, 10, 60, 150 ;
soil_lpdg_bot = 5, 10, 30, 200, 90, 180 ;
lpdg_bot = 0, 0 ;
lpdg_top = 30, 180 ;
fhg = 2, 10 ;
fh = 1829, 2743, 3658 ;
ls_all_top = 0.0, 0.33, 0.44, 0.44, 0.47, 0.84, 0.72 ;
ls_all_bot = 1.0, 1.0, 1.0, 0.72, 1.0, 0.98, 0.94 ;
sigma = 0.9950 ;
model_id = 77, 81, 96;
// Navigation
nav_model = "GRIB1" ;
grid_type_code = 0 ;
grid_type = "Latitude/Longitude" ;
grid_name = "Global 1.0 x 1.0 degree grid" ;
grid_center = 7 ; // NMC
grid_number = 3 ;
i_dim = "lon" ;
j_dim = "lat" ;
Ni = 360 ;
Nj = 181 ;
La1 = 90.0 ;
Lo1 = 0.0 ;
La2 = -90.0 ;
Lo2 = 359.0 ;
Di = 1.0 ;
Dj = 1.0 ;
ResCompFlag = 0x80 ;
lon = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
30, 31, 32, 33, 34, 35, 36, 37, 38, 39,
40, 41, 42, 43, 44, 45, 46, 47, 48, 49,
50, 51, 52, 53, 54, 55, 56, 57, 58, 59,
60, 61, 62, 63, 64, 65, 66, 67, 68, 69,
70, 71, 72, 73, 74, 75, 76, 77, 78, 79,
80, 81, 82, 83, 84, 85, 86, 87, 88, 89,
90, 91, 92, 93, 94, 95, 96, 97, 98, 99,
100, 101, 102, 103, 104, 105, 106, 107, 108, 109,
110, 111, 112, 113, 114, 115, 116, 117, 118, 119,
120, 121, 122, 123, 124, 125, 126, 127, 128, 129,
130, 131, 132, 133, 134, 135, 136, 137, 138, 139,
140, 141, 142, 143, 144, 145, 146, 147, 148, 149,
150, 151, 152, 153, 154, 155, 156, 157, 158, 159,
160, 161, 162, 163, 164, 165, 166, 167, 168, 169,
170, 171, 172, 173, 174, 175, 176, 177, 178, 179,
180, 181, 182, 183, 184, 185, 186, 187, 188, 189,
190, 191, 192, 193, 194, 195, 196, 197, 198, 199,
200, 201, 202, 203, 204, 205, 206, 207, 208, 209,
210, 211, 212, 213, 214, 215, 216, 217, 218, 219,
220, 221, 222, 223, 224, 225, 226, 227, 228, 229,
230, 231, 232, 233, 234, 235, 236, 237, 238, 239,
240, 241, 242, 243, 244, 245, 246, 247, 248, 249,
250, 251, 252, 253, 254, 255, 256, 257, 258, 259,
260, 261, 262, 263, 264, 265, 266, 267, 268, 269,
270, 271, 272, 273, 274, 275, 276, 277, 278, 279,
280, 281, 282, 283, 284, 285, 286, 287, 288, 289,
290, 291, 292, 293, 294, 295, 296, 297, 298, 299,
300, 301, 302, 303, 304, 305, 306, 307, 308, 309,
310, 311, 312, 313, 314, 315, 316, 317, 318, 319,
320, 321, 322, 323, 324, 325, 326, 327, 328, 329,
330, 331, 332, 333, 334, 335, 336, 337, 338, 339,
340, 341, 342, 343, 344, 345, 346, 347, 348, 349,
350, 351, 352, 353, 354, 355, 356, 357, 358, 359 ;
lat = 90, 89, 88, 87, 86, 85, 84, 83, 82, 81,
80, 79, 78, 77, 76, 75, 74, 73, 72, 71,
70, 69, 68, 67, 66, 65, 64, 63, 62, 61,
60, 59, 58, 57, 56, 55, 54, 53, 52, 51,
50, 49, 48, 47, 46, 45, 44, 43, 42, 41,
40, 39, 38, 37, 36, 35, 34, 33, 32, 31,
30, 29, 28, 27, 26, 25, 24, 23, 22, 21,
20, 19, 18, 17, 16, 15, 14, 13, 12, 11,
10, 9, 8, 7, 6, 5, 4, 3, 2, 1,
0, -1, -2, -3, -4, -5, -6, -7, -8, -9,
-10, -11, -12, -13, -14, -15, -16, -17, -18, -19,
-20, -21, -22, -23, -24, -25, -26, -27, -28, -29,
-30, -31, -32, -33, -34, -35, -36, -37, -38, -39,
-40, -41, -42, -43, -44, -45, -46, -47, -48, -49,
-50, -51, -52, -53, -54, -55, -56, -57, -58, -59,
-60, -61, -62, -63, -64, -65, -66, -67, -68, -69,
-70, -71, -72, -73, -74, -75, -76, -77, -78, -79,
-80, -81, -82, -83, -84, -85, -86, -87, -88, -89,
-90 ;
}
/*
* Copyright 1995, University Corporation for Atmospheric Research.
* See ../COPYRIGHT file for copying and redistribution conditions.
*/
/* $Id: levels.h,v 1.12 2000/02/15 18:25:39 rkambic Exp $ */
/* GRIB levels */
#ifndef LEVELS_H_
#define LEVELS_H_
#define LEVEL_SURFACE 1 /* surface of the Earth */
#define LEVEL_CLOUD_BASE 2 /* cloud base level */
#define LEVEL_CLOUD_TOP 3 /* cloud top level */
#define LEVEL_ISOTHERM 4 /* 0 degree isotherm level */
#define LEVEL_ADIABAT 5 /* adiabatic condensation level */
#define LEVEL_MAX_WIND 6 /* maximium wind speed level */
#define LEVEL_TROP 7 /* at the tropopause */
#define LEVEL_TOP 8 /* nominal top of atmosphere */
#define LEVEL_SEABOT 9 /* sea bottom */
#define LEVEL_TMPL 20 /* temperature in 1/100 K (2 octets) */
#define LEVEL_ISOBARIC 100 /* isobaric level */
#define LEVEL_LISO 101 /* layer between two isobaric levels */
#define LEVEL_MEAN_SEA 102 /* mean sea level */
#define LEVEL_FH 103 /* fixed height level */
#define LEVEL_LFHM 104 /* layer between 2 height levels above MSL */
#define LEVEL_FHG 105 /* fixed height above ground */
#define LEVEL_LFHG 106 /* layer between 2 height levels above ground */
#define LEVEL_SIGMA 107 /* sigma level */
#define LEVEL_LS 108 /* layer between 2 sigma levels */
#define LEVEL_HY 109 /* Hybrid level */
#define LEVEL_LHY 110 /* Layer between 2 hybrid levels */
#define LEVEL_Bls 111 /* Depth below land surface */
#define LEVEL_LBls 112 /* Layer between 2 depths below land surface */
#define LEVEL_ISEN 113 /* Isentropic (theta) level */
#define LEVEL_LISEN 114 /* Layer between 2 isentropic (theta) levels */
#define LEVEL_PDG 115 /* level at specified pressure difference from
ground */
#define LEVEL_LPDG 116 /* layer between levels at specif. pressure
diffs from ground */
#define LEVEL_PV 117 /* potential vorticity */
#define LEVEL_ETAL 119 /* ETA value in 1/10000 */
#define LEVEL_LETA 120 /* layer between two ETA levels */
#define LEVEL_LISH 121 /* layer between 2 isobaric surfaces (high
precision) */
#define LEVEL_FHGH 125 /* height level above ground (high precision) */
#define LEVEL_LSH 128 /* layer between 2 sigma levels (high
precision) */
#define LEVEL_LISM 141 /* layer between 2 isobaric surfaces (mixed
precision) */
#define LEVEL_DBS 160 /* depth below sea level */
#define LEVEL_ATM 200 /* entire atmosphere considered as a single
layer */
#define LEVEL_OCEAN 201 /* entire ocean considered as a single layer */
#define LEVEL_HTFL 204 /* Highest tropospheric freezing level */
#define LEVEL_BCBL 209 /* Boundary layer cloud bottom level */
#define LEVEL_BCTL 210 /* Boundary layer cloud top level */
#define LEVEL_BCY 211 /* Boundary layer cloud layer */
#define LEVEL_LCBL 212 /* Low cloud bottom level */
#define LEVEL_LCTL 213 /* Low cloud top level */
#define LEVEL_LCY 214 /* Low cloud layer */
#define LEVEL_MCBL 222 /* Middle cloud bottom level */
#define LEVEL_MCTL 223 /* Middle cloud top level */
#define LEVEL_MCY 224 /* Middle cloud layer */
#define LEVEL_HCBL 232 /* High cloud bottom level */
#define LEVEL_HCTL 233 /* High cloud top level */
#define LEVEL_HCY 234 /* Highcloud layer */
#define LEVEL_CCBL 242 /* convective cloud bottom level */
#define LEVEL_CCTL 243 /* convective cloud top level */
#define LEVEL_CCY 244 /* convective cloud layer */
#define LEVEL_MTHE 246 /* maximum equivalent potential temperature
level */
#define LEVEL_EHLT 247 /* equilibrium level */
#define LEVEL_FL 9999 /* FAA flight level */
#ifdef __cplusplus
extern "C" double mblev(int* levels);
extern "C" char* levelname(int level);
extern "C" char* levelsuffix(int level);
extern "C" char* levelunits(int level);
extern "C" long level_index(double level, float* levels, long nlevels);
extern "C" long layer_index(double top, double bot, float* tops,
float* bots, long nlayers);
extern "C" int level1(int flag, int* levels);
extern "C" int level2(int flag, int* levels);
#elif defined(__STDC__)
extern double mblev(int* levels);
extern char* levelname(int level);
extern char* levelsuffix(int level);
extern char* levelunits(int level);
extern long level_index(double level, float* levels, long nlevels);
extern long layer_index(double top, double bot, float* tops,
float* bots, long nlayers);
extern int level1(int flag, int* levels);
extern int level2(int flag, int* levels);
#else
extern double mblev ( /* int* levels */ );
#define const
extern char* levelname ( /* int level */ );
extern char* levelsuffix ( /* int level */ );
extern char* levelunits ( /* int level */ );
extern long level_index ( /* double level, float* levels, long nlevels */ );
extern long layer_index( /* double top, double bot, float* tops,
float* bots, long nlayers */ );
extern int level1( /* int flag, int* levels */ );
extern int level2( /* int flag, int* levels */ );
#endif
#endif /* LEVELS_H_ */
/*
* Copyright 1996 University Corporation for Atmospheric Research
* See ../COPYRIGHT file for copying and redistribution conditions.
*/
/* $Id: levels.c,v 1.16 2000/02/15 18:25:17 rkambic Exp $ */
#include "ulog.h"
#include "levels.h"
#include <math.h>
#define float_near(x,y) ((float)((y) + 0.1*fabs((x)-(y))) == (float)(y)) /*
true if x is "close to" y */
/*
* Atmospheric level in mb from two 8-bit integers in GRIB product
*/
double
mblev(levels)
int* levels;
{
return 256.*levels[0] + levels[1];
}
/*
* Search in table for specified level. Returns index
* of value in table, or -1 if the value was not found.
*/
long
level_index(level, levels, nlevels)
double level; /* level sought */
float* levels; /* table of levels */
long nlevels; /* number of levels in table */
{
int nn = 0;
while (nlevels--) {
float ll;
ll = *levels++;
if (float_near(level, ll)) /* true if level "is close to" ll */
break;
nn++;
}
return nlevels==-1 ? -1 : nn;
}
/*
* Search in table for specified layer. Returns index
* of value in table, or -1 if the value was not found.
*/
long
layer_index(top, bot, tops, bots, nlayers)
double top; /* top of layer sought */
double bot; /* bottom of layer sought */
float* tops; /* table of layer tops */
float* bots; /* table of layer bottoms */
long nlayers; /* number of layers in table */
{
int nn = 0;
while (nlayers--) {
if (float_near(top, *tops) && float_near(bot, *bots))
break;
tops++;
bots++;
nn++;
}
return nlayers==-1 ? -1 : nn;
}
/*
* Return name for GRIB level, given GRIB level code.
*/
char *
levelname(ii)
int ii;
{
switch(ii){
case LEVEL_SURFACE:
return "Surface";
case LEVEL_CLOUD_BASE:
return "Cloud Base";
case LEVEL_CLOUD_TOP:
return "Cloud Top";
case LEVEL_ISOTHERM:
return "0 Isotherm";
case LEVEL_ADIABAT:
return "Adiabatic Condensation";
case LEVEL_MAX_WIND:
return "Maximum Wind";
case LEVEL_TROP:
return "Tropopause";
case LEVEL_TOP:
return "Top of Atmosphere";
case LEVEL_SEABOT:
return "Sea Bottom";
case LEVEL_TMPL:
return "Temperature in 1/100 K";
case LEVEL_ISOBARIC:
return "Isobaric";
case LEVEL_LISO:
return "Layer Between Two Isobaric";
case LEVEL_MEAN_SEA:
return "Mean Sea";
case LEVEL_FH:
return "Fixed Height";
case LEVEL_LFHM:
return "Layer Between Two Heights Above MSL";
case LEVEL_FHG:
return "Fixed Height Above Ground";
case LEVEL_LFHG:
return "Layer Between Two Fixed Heights Above Ground";
case LEVEL_SIGMA:
return "Sigma";
case LEVEL_LS:
return "Layer Between Two Sigma";
case LEVEL_HY:
return "Hybrid level";
case LEVEL_LHY:
return "Layer between 2 hybrid levels";
case LEVEL_Bls:
return "Below Land Surface";
case LEVEL_LBls:
return "Layer Between Two Depths Below Land Surface";
case LEVEL_ISEN:
return "Isentropic (theta) level";
case LEVEL_LISEN:
return "Layer between 2 isentropic (theta) levels";
case LEVEL_PDG:
return "level at specified pressure difference from ground to level";
case LEVEL_LPDG:
return "layer between 2 levels at specified pressure differences from
ground to levels";
case LEVEL_PV:
return "potential vorticity";
case LEVEL_ETAL:
return "ETA level";
case LEVEL_LETA:
return "Layer between two ETA levels";
case LEVEL_LISH:
return "Layer Between Two Isobaric Surfaces, High Precision";
case LEVEL_FHGH:
return "Height level above ground (high precision)";
case LEVEL_LSH:
return "Layer Between Two Sigma Levels, High Precision";
case LEVEL_LISM:
return "Layer Between Two Isobaric Surfaces, Mixed Precision";
case LEVEL_DBS:
return "Depth Below Sea";
case LEVEL_ATM:
return "Entire atmosphere considered as a single layer";
case LEVEL_OCEAN:
return "Entire ocean considered as a single layer";
case LEVEL_HTFL:
return "Highest tropospheric freezing level";
case LEVEL_BCBL:
return "Boundary layer cloud bottom level";
case LEVEL_BCTL:
return "Boundary layer cloud top level";
case LEVEL_BCY:
return "Boundary layer cloud layer";
case LEVEL_LCBL:
return "Low cloud bottom level";
case LEVEL_LCTL:
return "Low cloud top level";
case LEVEL_LCY:
return "Low cloud layer";
case LEVEL_MCBL:
return "Middle cloud bottom level";
case LEVEL_MCTL:
return "Middle cloud top level";
case LEVEL_MCY:
return "Middle cloud layer";
case LEVEL_HCBL:
return "High cloud bottom level";
case LEVEL_HCTL:
return "High cloud top level";
case LEVEL_HCY:
return "Highcloud layer";
case LEVEL_CCBL:
return "Convective cloud bottom level";
case LEVEL_CCTL:
return "Convective cloud top level";
case LEVEL_CCY:
return "Convective cloud layer";
case LEVEL_MTHE:
return "maximum equivalent potential temperature level";
case LEVEL_EHLT:
return "equilibrium level";
case LEVEL_FL:
return "flight_level";
}
/* default */
uerror("unknown level: %d", ii);
return "reserved or unknown";
}
/*
* Returns level suffix, used in netCDF names for variables on special
* levels and in "gribdump -b" level abbreviations.
*/
char *
levelsuffix(lev)
int lev;
{
/* Note: If any suffixes are added or
changed, they must be added or changed in
the function grib_pcode() as well. */
switch(lev) {
case LEVEL_SURFACE: return "sfc"; /* surface of the Earth */
case LEVEL_CLOUD_BASE: return "clbs"; /* cloud base level */
case LEVEL_CLOUD_TOP: return "cltp"; /* cloud top level */
case LEVEL_ISOTHERM: return "frzlvl"; /* 0 degree isotherm level */
case LEVEL_ADIABAT: return "adcn"; /* adiabatic condensation level */
case LEVEL_MAX_WIND: return "maxwind"; /* maximium wind speed level */
case LEVEL_TROP: return "trop"; /* at the tropopause */
case LEVEL_TOP: return "topa"; /* nominal top of atmosphere */
case LEVEL_SEABOT: return "sbot"; /* sea bottom */
case LEVEL_TMPL: return "tmpl"; /* temperature in 1/100 K */
case LEVEL_ISOBARIC: return ""; /* isobaric level */
case LEVEL_LISO: return "liso"; /* layer between two isobaric levels */
case LEVEL_MEAN_SEA: return "msl"; /* mean sea level */
case LEVEL_FH: return "fh"; /* fixed height level */
case LEVEL_LFHM: return "lfhm"; /* layer between 2 height levels above MSL
*/
case LEVEL_FHG: return "fhg"; /* fixed height above ground */
case LEVEL_LFHG: return "lfhg"; /* layer between 2 height levels above
ground */
case LEVEL_SIGMA: return "sigma"; /* sigma level */
case LEVEL_LS: return "ls"; /* layer between 2 sigma levels */
case LEVEL_HY: return "hybr"; /* Hybrid level */
case LEVEL_LHY: return "lhyb"; /* Layer between 2 hybrid levels */
case LEVEL_Bls: return "bls"; /* Depth below land surface */
case LEVEL_LBls: return "lbls"; /* Layer between 2 depths below land
surface */
case LEVEL_ISEN: return "isen"; /* Isentropic (theta) level */
case LEVEL_LISEN: return "lisn"; /* Layer between 2 isentropic (theta)
levels */
case LEVEL_PDG: return "pdg"; /* level at specified pressure difference
from ground */
case LEVEL_LPDG: return "lpdg"; /* layer between levels at specif. pressure
diffs from ground */
case LEVEL_PV: return "pv"; /* level of specified potential vorticity */
case LEVEL_ETAL: return "etal"; /* ETA level */
case LEVEL_LETA: return "leta"; /* layer between 2 ETA levels */
case LEVEL_LISH: return "lish"; /* layer between 2 isobaric surfaces (high
precision) */
case LEVEL_FHGH: return "fhgh"; /* height level above ground (high
precision) */
case LEVEL_LSH: return "lsh"; /* layer between 2 sigma levels (high
precision) */
case LEVEL_LISM: return "lism"; /* layer between 2 isobaric surfaces (mixed
precision) */
case LEVEL_DBS: return "dbs"; /* depth below sea level */
case LEVEL_ATM: return "atm"; /* entire atmosphere considered as a single
layer */
case LEVEL_OCEAN: return "ocn"; /* entire ocean considered as a single
layer */
case LEVEL_HTFL: return "htfl"; /* Highest tropospheric freezing level */
case LEVEL_BCBL: return "bcbl"; /* Boundary layer cloud bottom level */
case LEVEL_BCTL: return "bctl"; /* Boundary layer cloud top level */
case LEVEL_BCY: return "bcy"; /* Boundary layer cloud layer */
case LEVEL_LCBL: return "lcbl"; /* Low cloud bottom level */
case LEVEL_LCTL: return "lctl"; /* Low cloud top level */
case LEVEL_LCY: return "lcy"; /* Low cloud layer */
case LEVEL_MCBL: return "mcbl"; /* Middle cloud bottom level */
case LEVEL_MCTL: return "mctl"; /* Middle cloud top level */
case LEVEL_MCY: return "mcy"; /* Middle cloud layer */
case LEVEL_HCBL: return "hcbl"; /* High cloud bottom level */
case LEVEL_HCTL: return "hctl"; /* High cloud top level */
case LEVEL_HCY: return "hcy"; /* Highcloud layer */
case LEVEL_CCBL: return "ccbl"; /* convective cloud bottom height */
case LEVEL_CCTL: return "cctl"; /* convective cloud top height */
case LEVEL_CCY: return "ccy"; /* convective cloud layer */
case LEVEL_MTHE: return "mthe"; /* maximum equivalent potential temperature
pressure */
case LEVEL_EHLT: return "ehlt"; /* equilibrium level height */
case LEVEL_FL: return "fl"; /* FAA flight level */
}
/* default: */
uerror("bad level flag: %d", lev);
return "";
}
/*
* Returns int for first level (if 2 levels) or level (if only one level)
*/
int
level1(flag, ii)
int flag; /* GRIB level flag */
int *ii; /* GRIB levels */
{
switch(flag){
case LEVEL_SURFACE:
case LEVEL_CLOUD_BASE:
case LEVEL_CLOUD_TOP:
case LEVEL_ISOTHERM:
case LEVEL_ADIABAT:
case LEVEL_MAX_WIND:
case LEVEL_TROP:
case LEVEL_TOP:
case LEVEL_SEABOT:
case LEVEL_MEAN_SEA:
case LEVEL_ATM:
case LEVEL_OCEAN:
case LEVEL_HTFL:
case LEVEL_LCBL:
case LEVEL_LCTL:
case LEVEL_LCY:
case LEVEL_MCBL:
case LEVEL_MCTL:
case LEVEL_MCY:
case LEVEL_HCBL:
case LEVEL_HCTL:
case LEVEL_HCY:
case LEVEL_FL:
return 0;
case LEVEL_TMPL:
case LEVEL_ISOBARIC:
case LEVEL_FH:
case LEVEL_FHG:
case LEVEL_SIGMA:
case LEVEL_HY:
case LEVEL_Bls:
case LEVEL_ISEN:
case LEVEL_PDG:
case LEVEL_PV:
case LEVEL_ETAL:
case LEVEL_FHGH:
case LEVEL_DBS:
return 256*ii[0]+ii[1]; /* 2-octet level */
case LEVEL_LISO:
case LEVEL_LFHM:
case LEVEL_LFHG:
case LEVEL_LS:
case LEVEL_LHY:
case LEVEL_LBls:
case LEVEL_LISEN:
case LEVEL_LPDG:
case LEVEL_LETA:
case LEVEL_LISH:
case LEVEL_LSH:
case LEVEL_LISM:
return ii[0]; /* 1-octet level */
}
/* default */
uerror("unknown level: %d", ii);
return -1;
}
/*
* Returns int for second level (if 2 levels) or 0 (if only one level)
*/
int
level2(flag, ii)
int flag; /* GRIB level flag */
int *ii; /* GRIB levels */
{
switch(flag){
case LEVEL_SURFACE:
case LEVEL_CLOUD_BASE:
case LEVEL_CLOUD_TOP:
case LEVEL_ISOTHERM:
case LEVEL_ADIABAT:
case LEVEL_MAX_WIND:
case LEVEL_TROP:
case LEVEL_TOP:
case LEVEL_SEABOT:
case LEVEL_MEAN_SEA:
case LEVEL_ATM:
case LEVEL_OCEAN:
case LEVEL_HTFL:
case LEVEL_LCBL:
case LEVEL_LCTL:
case LEVEL_LCY:
case LEVEL_MCBL:
case LEVEL_MCTL:
case LEVEL_MCY:
case LEVEL_HCBL:
case LEVEL_HCTL:
case LEVEL_HCY:
case LEVEL_FL:
case LEVEL_TMPL:
case LEVEL_ISOBARIC:
case LEVEL_FH:
case LEVEL_FHG:
case LEVEL_SIGMA:
case LEVEL_HY:
case LEVEL_Bls:
case LEVEL_ISEN:
case LEVEL_PDG:
case LEVEL_PV:
case LEVEL_ETAL:
case LEVEL_FHGH:
case LEVEL_DBS:
return 0;
case LEVEL_LISO:
case LEVEL_LFHM:
case LEVEL_LFHG:
case LEVEL_LS:
case LEVEL_LHY:
case LEVEL_LBls:
case LEVEL_LISEN:
case LEVEL_LPDG:
case LEVEL_LETA:
case LEVEL_LISH:
case LEVEL_LSH:
case LEVEL_LISM:
return ii[1];
}
/* default */
uerror("unknown level: %d", ii);
return -1;
}
/*
* Return GRIB units (as a string) for various kinds of levels.
*/
char *
levelunits(ii)
{
switch(ii){
case LEVEL_SURFACE:
case LEVEL_CLOUD_BASE:
case LEVEL_CLOUD_TOP:
case LEVEL_ISOTHERM:
case LEVEL_ADIABAT:
case LEVEL_MAX_WIND:
case LEVEL_TROP:
case LEVEL_TOP:
case LEVEL_SEABOT:
case LEVEL_MEAN_SEA:
case LEVEL_HY:
case LEVEL_LHY:
case LEVEL_ATM:
case LEVEL_OCEAN:
case LEVEL_HTFL:
case LEVEL_LCBL:
case LEVEL_LCTL:
case LEVEL_LCY:
case LEVEL_MCBL:
case LEVEL_MCTL:
case LEVEL_MCY:
case LEVEL_HCBL:
case LEVEL_HCTL:
case LEVEL_HCY:
case LEVEL_FL:
return "" ;
case LEVEL_ISOBARIC:
case LEVEL_PDG:
case LEVEL_LPDG:
case LEVEL_LISH:
return "hPa";
case LEVEL_LISO:
case LEVEL_LISM:
return "kPa";
case LEVEL_FH:
case LEVEL_FHG:
case LEVEL_DBS:
return "meters";
case LEVEL_LFHM:
case LEVEL_LFHG:
return "hm" ;
case LEVEL_SIGMA:
return ".0001"; /* dimensionless */
case LEVEL_LS:
case LEVEL_LETA:
return ".01"; /* dimensionless */
case LEVEL_Bls:
case LEVEL_LBls:
case LEVEL_FHGH:
return "cm";
case LEVEL_ISEN:
case LEVEL_LISEN:
return "degK";
case LEVEL_TMPL:
return ".01 degK";
case LEVEL_LSH:
return ".001";
case LEVEL_PV:
return ".000001 K m2/kg/sec";
case LEVEL_ETAL:
return ".0001";
}
/* default */
return "unknown" ;
}
/*
* Copyright 1996 University Corporation for Atmospheric Research
* See ../COPYRIGHT file for copying and redistribution conditions.
*/
/* $Id: nc.c,v 1.33 2000/04/11 22:51:34 rkambic Exp $ */
#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <math.h>
#include <netcdf.h>
#include "ulog.h"
#include "mkdirs_open.h"
#include "nc.h"
#include "nuwg.h"
#include "emalloc.h"
#include "params.h"
#include "units.h"
#include "levels.h"
#include "timeunits.h"
#include "gbds.h" /* only for FILL_VAL */
#include "gdes.h"
#include "recs.h"
#ifndef FILL_NAME
#define FILL_NAME "_FillValue"
#endif
typedef struct levels_table {
int id; /* dimension id of netCDF dimension */
float *vals; /* levels */
long num; /* number of levels */
utUnit *bunitp; /* level units */
} levels_table;
typedef struct layers_table {
int id; /* dimension id of netCDF dimension */
float *tops; /* top levels of layers */
float *bots; /* bottom levels of layers */
long num; /* number of layers */
utUnit *bunitp; /* top (and bottom) level units */
} layers_table;
typedef struct tris_table {
int id; /* dimension id of netCDF dimension */
float *starts; /* starts of time ranges */
float *ends; /* ends of time ranges */
long num; /* number of ranges */
} tris_table;
#ifdef __STDC__
static ncdim* new_dim(int dimid);
static ncvar* new_var(char* ncname, int varid);
static void free_var(ncvar* var);
static void free_dim(ncdim* dim);
static char* parmname(ncfile* nc, int parm, int level);
static int make_ncfile(char* ncname, ncfile* out);
static void free_ncfile(ncfile* np);
static levels_table* getlevtab(ncfile* nc, ncvar* var);
static layers_table* getlaytab(ncfile* nc, ncvar* var);
static long getlev(product_data* pp, ncfile* nc,
ncvar* var);
static long gettri(product_data* pp, ncfile* nc,
ncvar* var);
static int make_var(char* ncname, int varid, ncvar* out);
static int var_as_int(ncfile* nc, enum ncpart comp, int* val);
static int var_as_float(ncfile* nc, enum ncpart comp, float* val);
static int var_as_lset(ncfile* nc, enum ncpart comp, lset* list);
static void varerr(ncfile* nc, enum ncpart comp);
static int make_navgrid(ncfile* nc, navinfo* nav);
static int make_navinfo(ncfile* nc, navinfo* nav);
static void free_navinfo(navinfo* np);
static navinfo* new_navinfo(ncfile* nc);
static int gd_fne_err(product_data* pp, ncfile* nc, enum ncpart comp, double
pval, double nval);
static int gd_ine_err(product_data* pp, ncfile* nc, enum ncpart comp, int pval,
int nval);
static int gd_igt_err(product_data* pp, ncfile* nc, enum ncpart comp, int pval,
int nval);
static int check_gd(product_data* pp, ncfile* nc);
static int check_pp(product_data* pp, ncfile* nc);
static int subgrid(ncfile* nc, product_data* pp, long* ix0, long* ix1);
#endif
int ncid; /* netCDF id of open netCDF output file.
* This is at file scope so routine registered
* with atexit() can get it to close file. */
static ncdim *
new_dim(dimid)
int dimid;
{
char dimname[MAX_NC_NAME];
long size;
ncdim *out = (ncdim *)emalloc(sizeof(ncdim));
if (ncdiminq(ncid, dimid, dimname, &size) == -1) {
return 0;
}
out->name = estrdup(dimname);
out->id = dimid;
/* We don't cache size because it may change,
e.g. adding records, and we don't want to
have to keep cached value up to date. To
get size, use ncdiminq(). */
return out;
}
static int
make_var(ncname, varid, out)
char *ncname; /* netCDF pathanme, only used in error msg */
int varid; /* variable ID */
ncvar *out; /* place to put constructed ncvar */
{
char varname[MAX_NC_NAME];
nc_type type;
int ndims;
int dims[MAX_VAR_DIMS];
int id;
if(ncvarinq(ncid, varid, varname, &type, &ndims, dims, (int *)0) == -1) {
return -1;
}
out->id = varid;
out->name = estrdup(varname);
out->type = type;
out->ndims = ndims;
out->dims = (int *)emalloc(ndims * sizeof(int));
for (id = 0; id < ndims; id++)
out->dims[id] = dims[id];
/* get value of _FillValue attribute, if any,
as a float */
{
nc_type atttype;
int attlen;
if (ncattinq(ncid, varid, FILL_NAME, &atttype, &attlen) == -1) {
out->fillval = 0; /* no fill-value attribute */
} else {
union {
char c;
unsigned char u;
short s;
long l;
float f;
double d;
}fillval;
out->fillval = (float *) emalloc(sizeof(float));
if(ncattget(ncid, varid, FILL_NAME, (void *)&fillval) == -1) {
uerror("%s: can't get % attribute for variable %s",
ncname, FILL_NAME, varname);
free(out->dims);
return -1;
}
switch(atttype) {
case NC_CHAR:
*out->fillval = fillval.c;
break;
case NC_BYTE:
*out->fillval = fillval.u;
break;
case NC_SHORT:
*out->fillval = fillval.s;
break;
case NC_LONG:
*out->fillval = fillval.l;
break;
case NC_FLOAT:
*out->fillval = fillval.f;
break;
case NC_DOUBLE:
*out->fillval = fillval.d;
break;
}
}
}
if (get_units(ncid, varid, &out->bunitp) == -1) {
uerror("%s: can't get units attribute for variable %s",
ncname, varname);
return -1;
}
if (out->bunitp && grib_pcode(varname) != -1) {
out->uc = uconv(varname, out->bunitp);
} else {
out->uc = 0;
}
return 0;
}
static void
free_var(var)
ncvar *var;
{
if (var) {
if(var->name)
free(var->name);
if(var->dims)
free(var->dims);
if(var->bunitp)
free(var->bunitp);
if(var->fillval)
free(var->fillval);
free(var);
}
}
/*
* Creates a new ncvar structure and fills it in with the information from
* the open netCDF file whose handle is ncid. A pointer to the structure is
* returned, or 0 on failure.
*/
static ncvar *
new_var(ncname, varid)
char *ncname;
{
ncvar *out;
if (varid == -1) /* handle common failure case with message
at higher level */
return 0;
out = (ncvar *)emalloc(sizeof(ncvar));
if (make_var(ncname, varid, out) != 0) {
free_var(out);
return 0;
}
return out;
}
static void
free_dim(dim)
ncdim *dim;
{
if (dim) {
if(dim->name)
free(dim->name);
free(dim);
}
}
/*
* Name of environment variable containing directory in which to search for
* CDL files if not found relative to directory from which this process is
* invoked.
*/
#define LDM_ETCDIR "LDM_ETCDIR"
/*
* Checks to see if netCDF file with specified name exists. If not,
* makes netCDF file from CDL template file.
* Returns netCDF file ID on success, or -1 on error.
*/
int
cdl_netcdf (cdlname, ncname)
char *cdlname; /* CDL file specifying netCDF structure */
char *ncname; /* filename of netcdf file to be created */
{
char cmnd[2*_POSIX_PATH_MAX+20];
char *cdlfile = cdlname;
static char *cdldir = 0;
char envcdl[_POSIX_PATH_MAX];
ncopts = 0; /* turn off netCDF error messages from
library, means we will have to check all
netCDF status returns and interpret */
if (access(ncname, (R_OK|W_OK)) != 0) { /* no writable ncname exists */
if (cdlfile == 0) {
uerror("%s doesn't exist, and didn't specify a CDL filename",
ncname);
return -1;
}
/* Try to create any directories in output path that don't exist */
if (diraccess(ncname, (R_OK | W_OK), !0) == -1) {
serror("can't access directories leading to %s", ncname);
return -1;
}
/* If CDL file not found, look in environment variable LDM_ETCDIR */
if (access(cdlname, R_OK) == -1) { /* no CDL file or unreadable */
if (cdldir == 0) { /* executed once, when cdldir first needed */
char *ldm_etcdir = getenv(LDM_ETCDIR);
int slen;
if (ldm_etcdir == 0) {
uerror("CDL file %s not found & LDM_ETCDIR not in
environment",
cdlname);
return -1;
}
slen = strlen(ldm_etcdir);
cdldir = (char *)emalloc((slen+2) * sizeof(char));
strcpy(cdldir, ldm_etcdir);
if (cdldir[slen-1] != '/') { /* append "/" to dir name */
strcat(cdldir, "/");
}
}
strcat(envcdl,cdldir);
strcat(envcdl,cdlname);
if (access(envcdl, R_OK) == -1) {
uerror("can't find CDL file %s, or unreadable", envcdl);
return -1;
}
cdlfile = envcdl;
}
(void) strcpy(cmnd, "ncgen -o ");
(void) strcat(cmnd, ncname);
(void) strcat(cmnd , " ");
(void) strcat(cmnd, cdlfile);
if (system(cmnd) != 0) {
uerror("can't run \"%s\"", cmnd);
return -1;
}
}
return ncopen(ncname, NC_WRITE);
}
void
setncid(id)
int id;
{
ncid = id;
}
int
getncid()
{
return ncid;
}
/* Close open netCDF file, if any */
void
nccleanup()
{
ncclose(ncid);
}
/*
* Creates and returns a pointer to a one-dimensional table of levels
* from the currently-open netCDF file (handle ncid). Returns 0 on
* failure.
*/
static levels_table*
getlevtab(nc, var)
ncfile *nc; /* currently open netCDF file */
ncvar *var; /* level variable */
{
int did;
levdim *lp = nc->levdims;
/* See if appropriate levels table already exists for this file */
if (var->ndims < 4) {
uerror("%s: variable %s has too few dimensions for a level",
nc->ncname, var->name);
return 0;
}
did = var->dims[1]; /* level dimension */
while (lp) {
if(did == lp->levtab->id) { /* found it */
return lp->levtab;
}
lp = lp->next;
}
/* Not, there, so we must create it */
lp = (levdim *)emalloc(sizeof(levdim));
lp->levtab = (levels_table *)emalloc(sizeof(levels_table));
lp->levtab->id = did;
lp->next = nc->levdims;
/* Initialize array of levels */
{
levels_table *out = lp->levtab;
char levname[MAX_NC_NAME];
int levvarid; /* variable id of level variable */
ncvar *lev; /* level variable */
long start = 0;
/* get number of levels */
if (ncdiminq(ncid, out->id, levname, &out->num) == -1) {
uerror("%s: can't get number of %s levels", nc->ncname, var->name);
return 0;
}
out->vals = (float *) emalloc(out->num * sizeof(float));
levvarid = ncvarid(ncid, levname);
if(levvarid == -1) {
uerror("%s: No %s coordinate variable for %s level",
nc->ncname, levname, var->name);
return 0;
}
lev = nc->vars[levvarid];
/* Check consistency of lev variable */
if (strcmp(lev->name, levname) != 0 ||
lev->type != NC_FLOAT ||
lev->ndims != 1 ||
lev->dims[0] != out->id) {
uerror("%s: variable %s must be float %s(%s)", nc->ncname,
levname, lev->name, lev->name);
return 0;
}
if(get_units(ncid, levvarid, &out->bunitp) == -1) {
uerror("%s: error getting units attribute for %s",
nc->ncname, levname);
return 0;
}
if(ncvarget(ncid, levvarid, &start, &out->num,
(void *)out->vals) ==-1) {
uerror("%s: no %s variable for level",
nc->ncname, levname);
return 0;
}
}
nc->levdims = lp; /* if all goes well */
return lp->levtab;
}
/*
* Creates and returns a pointer to a one-dimensional table of layers
* from the currently-open netCDF file (handle ncid). Returns 0 on
* failure.
*/
static layers_table*
getlaytab(nc, var)
ncfile *nc; /* currently open netCDF file */
ncvar *var; /* layer variable */
{
int did;
laydim *lp = nc->laydims;
/* See if appropriate layers table already exists for this file */
if (var->ndims < 4) {
uerror("%s: variable %s has too few dimensions for a layer",
nc->ncname, var->name);
return 0;
}
did = var->dims[1]; /* layer dimension */
while (lp) {
if(did == lp->laytab->id) { /* found it */
return lp->laytab;
}
lp = lp->next;
}
/* Not, there, so we must create it */
lp = (laydim *)emalloc(sizeof(laydim));
lp->laytab = (layers_table *)emalloc(sizeof(layers_table));
lp->laytab->id = did;
lp->next = nc->laydims;
/* Initialize array of layers */
{
layers_table *out = lp->laytab;
char layname[MAX_NC_NAME];
char topname[MAX_NC_NAME];
char botname[MAX_NC_NAME];
int topvarid; /* variable id of layer top variable */
int botvarid; /* variable id of layer bottom variable */
ncvar *top; /* layer top variable */
ncvar *bot; /* layer bottom variable */
long start = 0;
/* get number of layers */
if (ncdiminq(ncid, out->id, layname, &out->num) == -1) {
uerror("%s: can't get number of %s layers", nc->ncname, var->name);
return 0;
}
if (strlen(layname) + strlen("_top") > (size_t) MAX_NC_NAME) {
uerror("%s: name of layer dimension too long (%s)",
nc->ncname, layname);
return 0;
}
out->tops = (float *) emalloc(out->num * sizeof(float));
strcpy(topname, layname);
strcat(topname, "_top");
topvarid = ncvarid(ncid, topname);
if(topvarid == -1) {
uerror("%s: No %s coordinate variable for %s layer top",
nc->ncname, layname, var->name);
return 0;
}
top = nc->vars[topvarid];
/* Check consistency of top variable */
if (strcmp(top->name, topname) != 0 ||
top->type != NC_FLOAT ||
top->ndims != 1 ||
top->dims[0] != out->id) {
uerror("%s: variable %s must be float %s(%s)", nc->ncname,
layname, top->name, top->name);
return 0;
}
if(get_units(ncid, topvarid, &out->bunitp) == -1) {
uerror("%s: error getting units attribute for %s",
nc->ncname, topname);
return 0;
}
if(ncvarget(ncid, topvarid, &start, &out->num,
(void *)out->tops) ==-1) {
uerror("%s: no %s variable for top of layer",
nc->ncname, topname);
return 0;
}
out->bots = (float *) emalloc(out->num * sizeof(float));
strcpy(botname, layname);
strcat(botname, "_bot");
botvarid = ncvarid(ncid, botname);
if(botvarid == -1) {
uerror("%s: No %s coordinate variable for %s layer bot",
nc->ncname, layname, var->name);
return 0;
}
bot = nc->vars[botvarid];
/* Check consistency of bot variable */
if (strcmp(bot->name, botname) != 0 ||
bot->type != NC_FLOAT ||
bot->ndims != 1 ||
bot->dims[0] != out->id) {
uerror("%s: variable %s must be float %s(%s)", nc->ncname,
layname, bot->name, bot->name);
return 0;
}
if(ncvarget(ncid, botvarid, &start, &out->num,
(void *)out->bots) ==-1) {
uerror("%s: no %s variable for bottom of layer",
nc->ncname, botname);
return 0;
}
}
nc->laydims = lp; /* if all goes well */
return lp->laytab;
}
/*
* Handle levels
*/
static long
levaux(pp, nc, var)
product_data *pp; /* decoded GRIB data to be written */
ncfile *nc; /* netCDF file to be written */
ncvar *var; /* netCDF variable to be written */
{
int level_flg = pp->level_flg;
char *varname = var->name;
/*
Assumes level dimension is second dimension of a variable, and that
there is a coordinate variable associated with a level, values of
which get stored in the level table.
*/
int levdim;
double lev;
long levix;
levels_table *levtab = getlevtab(nc, var);
if (levtab == 0) { /* initialize table of levels */
return -1;
}
lev = level1(pp->level_flg, pp->level);
/* Must convert to units of level table */
{
utUnit bfunit;
char *funits = levelunits(level_flg);
double slope=1.;
double intercept=0.;
if(utScan(funits, &bfunit) != 0) { /* "from" unit */
uerror("Error parsing unit `%s' for level %s", funits, varname);
return -1;
}
if (levtab->bunitp) {
if(utConvert(&bfunit, levtab->bunitp, &slope, &intercept) ==
UT_ECONVERT) {
uerror("units `%s' not conformable with variable %s:units",
funits, varname);
return -1;
}
}
lev = slope * lev + intercept;
}
levix = level_index(lev, levtab->vals, levtab->num);
if (levix == -1) {
uerror("%s: In %s, no %f level for %s",
pp->header, nc->ncname, lev, varname);
return -1;
}
return levix;
}
/*
* Handle layers
*/
static long
layaux(pp, nc, var)
product_data *pp; /* decoded GRIB data to be written */
ncfile *nc; /* netCDF file to be written */
ncvar *var; /* netCDF variable to be written */
{
int layer_flg = pp->level_flg;
char *varname = var->name;
/*
Assumes layer dimension is second dimension of a variable, and that
there are _top and _bot levels associated with each layer, values of
which which get stored in the layer table.
*/
int laydim;
double top;
double bot;
long layix;
layers_table *laytab = getlaytab(nc, var);
if (laytab == 0) {
return -1;
}
top = pp->level[0];
bot = pp->level[1];
/* Must convert top,bot to units of layer tables */
{
utUnit bfunit;
char *funits = levelunits(layer_flg);
double slope=1.;
double intercept=0.;
if(utScan(funits, &bfunit) != 0) { /* "from" unit */
uerror("Error parsing unit `%s' for level %s", funits, varname);
return -1;
}
if(laytab->bunitp) {
if(utConvert(&bfunit, laytab->bunitp, &slope, &intercept) ==
UT_ECONVERT) {
uerror("units `%s' not conformable with variable %s:units",
funits, varname);
return -1;
}
}
top = slope * top + intercept;
bot = slope * bot + intercept;
}
layix = layer_index(top, bot, laytab->tops, laytab->bots, laytab->num);
if (layix == -1) {
uerror("%s: In %s, no (%g,%g) level for %s",
pp->header, nc->ncname, top, bot, varname);
return -1;
}
return layix;
}
/*
* Return netCDF level dimension index appropriate for decoded GRIB
* product. Returns -2 if no level dimension appropriate (e.g. for surface
* variables) or -1 in case of failure.
*/
static long
getlev(pp, nc, var)
product_data *pp; /* decoded GRIB data to be written */
ncfile *nc; /* netCDF file to be written */
ncvar *var; /* netCDF variable to be written */
{
switch (pp->level_flg) {
/* Levels */
case LEVEL_ISOBARIC:
case LEVEL_FHG:
case LEVEL_SIGMA:
case LEVEL_HY:
case LEVEL_FH:
case LEVEL_Bls:
case LEVEL_ISEN:
case LEVEL_PDG:
case LEVEL_FHGH:
case LEVEL_DBS:
case LEVEL_FL:
return levaux(pp, nc, var);
/* Layers */
case LEVEL_LBls:
case LEVEL_LFHG:
case LEVEL_LFHM:
case LEVEL_LHY:
case LEVEL_LISEN:
case LEVEL_LISH:
case LEVEL_LISM:
case LEVEL_LISO:
case LEVEL_LPDG:
case LEVEL_LS:
case LEVEL_LSH:
return layaux(pp, nc, var);
/* Special levels, just one so no dimension needed */
case LEVEL_SURFACE:
case LEVEL_CLOUD_BASE:
case LEVEL_CLOUD_TOP:
case LEVEL_ISOTHERM:
case LEVEL_ADIABAT:
case LEVEL_MAX_WIND:
case LEVEL_TROP:
case LEVEL_TOP:
case LEVEL_SEABOT:
case LEVEL_MEAN_SEA:
case LEVEL_ATM:
case LEVEL_OCEAN:
return -2;
}
/* default: */
return -1;
}
/*
* Get conventional netcdf variable name, with level indicator appended if
* appropriate, and see if it exists in the open netCDF file. If so return
* name. If not, return 0. The name is a pointer to a static string, so
* should be copied if needed beyond the next call to parmname.
*/
static char*
parmname(nc, parm, level)
ncfile *nc; /* netCDF file */
int parm; /* parameter code from GRIB product */
int level; /* level flag from GRIB product */
{
char *varname = grib_pname(parm); /* netcdf variable base name */
int ncopts_save = ncopts;
char *suffix;
static char string[MAX_NC_NAME];
char *name = string;
if (!varname) {
uerror("unrecognized GRIB parameter code %d", parm);
return 0;
}
/* Add level modifier to name, if appropriate */
suffix = "";
suffix = levelsuffix(level);
strcpy(name, varname);
if( suffix == "" )
uerror("var name : %s", name);
/* The "_sfc", "_msl", and "_liso" suffixes are redundant for some
parameters, so we explicitly exclude those here */
if((level != LEVEL_SURFACE || !sfcparam(parm)) &&
(level != LEVEL_MEAN_SEA || !mslparam(parm)) &&
(level != LEVEL_LISO || !lisoparam(parm))) {
if (suffix[0] != '\0') {
strcat(name, "_");
strcat(name, suffix);
}
}
ncopts = 0;
if (ncvarid(ncid, name) == -1) {
uerror("%s: no variable %s in netCDF file", nc->ncname, name);
name = 0;
}
ncopts = ncopts_save;
return name;
}
/*
* Stores value of a netCDF variable identified by a NUWG conventional id.
* In case of failure, returns -1. The value is converted from
* whatever type is used for the netCDF variable.
*/
static int
var_as_int(nc, comp, val)
ncfile *nc;
enum ncpart comp;
int *val; /* where to store the resulting value */
{
ncvar *var = (ncvar *)emalloc(sizeof(ncvar));
long start[] = {0};
long count[] = {1};
double buf[1]; /* generic data buffer */
if(make_var(nc->ncname, nuwg_getvar(ncid, comp), var) == -1) {
free_var(var);
return -1;
}
if (ncvarget(ncid, var->id, start, count, (void *)buf) == -1) {
free_var(var);
return -1;
}
switch (var->type) { /* return the value as an int, no
matter how it is stored */
case NC_BYTE:
*val = *(unsigned char *) buf;
break;
case NC_CHAR:
*val = *(char *) buf;
break;
case NC_SHORT:
*val = *(short *) buf;
break;
case NC_LONG:
*val = *(nclong *) buf;
break;
case NC_FLOAT:
*val = *(float *) buf;
break;
case NC_DOUBLE:
*val = *(double *) buf;
break;
}
free_var(var);
return 0;
}
/*
* Stores value of a netCDF variable identified by a NUWG conventional id.
* In case of failure, returns -1. The value is converted from
* whatever type is used for the netCDF variable.
*/
static int
var_as_float(nc, comp, val)
ncfile *nc;
enum ncpart comp;
float *val;
{
ncvar *var = (ncvar *)emalloc(sizeof(ncvar));
long start[] = {0};
long count[] = {1};
double buf[1]; /* generic data buffer */
if(make_var(nc->ncname, nuwg_getvar(ncid, comp), var) == -1) {
return -1;
}
if (ncvarget(ncid, var->id, start, count, (void *)buf) == -1) {
free_var(var);
return -1;
}
switch (var->type) { /* return the value as a float, no
matter how it is stored */
case NC_BYTE:
*val = *(unsigned char *) buf;
break;
case NC_CHAR:
*val = *(char *) buf;
break;
case NC_SHORT:
*val = *(short *) buf;
break;
case NC_LONG:
*val = *(nclong *) buf;
break;
case NC_FLOAT:
*val = *(float *) buf;
break;
case NC_DOUBLE:
*val = *(double *) buf;
break;
}
free_var(var);
return 0;
}
/*
* Stores values of a netCDF variable (of longs) identified by a NUWG
* conventional id. Values are just stored in a list of longs, which can be
* used as a set in which values are looked up. In case of failure, returns
* -1.
*/
static int
var_as_lset(nc, comp, list)
ncfile *nc;
enum ncpart comp;
lset *list; /* where to store the resulting list */
{
ncvar *var = (ncvar *)emalloc(sizeof(ncvar));
static long start[MAX_VAR_DIMS];
static long count[MAX_VAR_DIMS];
long prod;
int i;
if(make_var(nc->ncname, nuwg_getvar(ncid, comp), var) == -1) {
return -1;
}
if (var->type != NC_LONG) {
uerror("%s: variable %s must be of type long", nc->ncname,
nuwg_name(comp));
free_var(var);
return -1;
}
prod=1;
for (i=0; i<var->ndims; i++) {
start[i] = 0;
if (ncdiminq(ncid, var->dims[i], (char *)0, &count[i]) == -1) {
uerror("%s: can't get size of dimension for %s", nc->ncname,
nuwg_name(comp));
free_var(var);
return -1;
}
prod *= count[i];
}
list->n = prod;
list->vals = (nclong *)emalloc(sizeof(nclong) * prod);
if (ncvarget(ncid, var->id, start, count, (void *)list->vals) == -1) {
uerror("%s: can't get values for %s", nc->ncname, nuwg_name(comp));
free_var(var);
free(list->vals);
return -1;
}
free_var(var);
return 0;
}
static void
varerr(nc,comp)
ncfile *nc;
enum ncpart comp;
{
uerror("%s: no variable for %s", nc->ncname, nuwg_name(comp));
}
/*
* Returns 0 on success, -1 on failure.
*/
static int
make_navgrid(nc, nav )
ncfile *nc; /* netCDF file */
navinfo *nav; /* where to put the navinfo */
{
int *ip;
float *fp;
switch(nav->grid_type_code) {
case GRID_LL:
case GRID_RLL:
case GRID_SLL:
case GRID_SRLL:
{
gdes_ll *gg = &nav->grid.ll;
if (var_as_int(nc, VAR_NI, &gg->ni) == -1) {
varerr(nc, VAR_NI);
return -1;
}
if (var_as_int(nc, VAR_NJ, &gg->nj) == -1) {
varerr(nc, VAR_NJ);
return -1;
}
if (var_as_float(nc, VAR_LA1, &gg->la1) == -1) {
varerr(nc, VAR_LA1);
return -1;
}
if (var_as_float(nc, VAR_LO1, &gg->lo1) == -1) {
varerr(nc, VAR_LO1);
return -1;
}
if (var_as_float(nc, VAR_LA2, &gg->la2) == -1) {
varerr(nc, VAR_LA2);
return -1;
}
if (var_as_float(nc, VAR_LO2, &gg->lo2) == -1) {
varerr(nc, VAR_LO2);
return -1;
}
if (var_as_float(nc, VAR_DI, &gg->di) == -1) {
varerr(nc, VAR_DI);
return -1;
}
if (var_as_float(nc, VAR_DJ, &gg->dj) == -1) {
varerr(nc, VAR_DJ);
return -1;
}
if (nav->grid_type_code == GRID_RLL || nav->grid_type_code ==
GRID_SRLL) {
gg->rot = (rotated *)emalloc(sizeof(rotated));
if (var_as_float(nc, VAR_ROTLAT, &gg->rot->lat) == -1) {
varerr(nc, VAR_ROTLAT);
return -1;
}
if (var_as_float(nc, VAR_ROTLON, &gg->rot->lon) == -1) {
varerr(nc, VAR_ROTLON);
return -1;
}
if (var_as_float(nc, VAR_ROTANGLE, &gg->rot->angle) == -1) {
varerr(nc, VAR_ROTANGLE);
return -1;
}
}
if (nav->grid_type_code == GRID_SLL || nav->grid_type_code ==
GRID_SRLL) {
gg->strch = (stretched *)emalloc(sizeof(stretched));
if (var_as_float(nc, VAR_STRETCHLAT, &gg->strch->lat) == -1) {
varerr(nc, VAR_STRETCHLAT);
return -1;
}
if (var_as_float(nc, VAR_STRETCHLON, &gg->strch->lon) == -1) {
varerr(nc, VAR_STRETCHLON);
return -1;
}
if (var_as_float(nc, VAR_STRETCHFACTOR, &gg->strch->factor) == -1) {
varerr(nc, VAR_STRETCHFACTOR);
return -1;
}
}
}
break;
case GRID_GAU:
case GRID_RGAU:
case GRID_SGAU:
case GRID_SRGAU:
{
gdes_gau *gg = &nav->grid.gau;
if (var_as_int(nc, VAR_NI, &gg->ni) == -1) {
varerr(nc, VAR_NI);
return -1;
}
if (var_as_int(nc, VAR_NJ, &gg->nj) == -1) {
varerr(nc, VAR_NJ);
return -1;
}
if (var_as_float(nc, VAR_LA1, &gg->la1) == -1) {
varerr(nc, VAR_LA1);
return -1;
}
if (var_as_float(nc, VAR_LO1, &gg->lo1) == -1) {
varerr(nc, VAR_LO1);
return -1;
}
if (var_as_float(nc, VAR_LA2, &gg->la2) == -1) {
varerr(nc, VAR_LA2);
return -1;
}
if (var_as_float(nc, VAR_LO2, &gg->lo2) == -1) {
varerr(nc, VAR_LO2);
return -1;
}
if (var_as_float(nc, VAR_DI, &gg->di) == -1) {
varerr(nc, VAR_DI);
return -1;
}
if (var_as_int(nc, VAR_N, &gg->n) == -1) {
varerr(nc, VAR_N);
return -1;
}
if (nav->grid_type_code == GRID_RLL || nav->grid_type_code ==
GRID_SRLL) {
gg->rot = (rotated *)emalloc(sizeof(rotated));
if (var_as_float(nc, VAR_ROTLAT, &gg->rot->lat) == -1) {
varerr(nc, VAR_ROTLAT);
return -1;
}
if (var_as_float(nc, VAR_ROTLON, &gg->rot->lon) == -1) {
varerr(nc, VAR_ROTLON);
return -1;
}
if (var_as_float(nc, VAR_ROTANGLE, &gg->rot->angle) == -1) {
varerr(nc, VAR_ROTANGLE);
return -1;
}
}
if (nav->grid_type_code == GRID_SLL || nav->grid_type_code ==
GRID_SRLL) {
gg->strch = (stretched *)emalloc(sizeof(stretched));
if (var_as_float(nc, VAR_STRETCHLAT, &gg->strch->lat) == -1) {
varerr(nc, VAR_STRETCHLAT);
return -1;
}
if (var_as_float(nc, VAR_STRETCHLON, &gg->strch->lon) == -1) {
varerr(nc, VAR_STRETCHLON);
return -1;
}
if (var_as_float(nc, VAR_STRETCHFACTOR, &gg->strch->factor) == -1) {
varerr(nc, VAR_STRETCHFACTOR);
return -1;
}
}
}
break;
case GRID_SPH:
case GRID_RSPH:
case GRID_SSPH:
case GRID_SRSPH:
{
gdes_sph *gg = &nav->grid.sph;
if (var_as_int(nc, VAR_J, &gg->j) == -1) {
varerr(nc, VAR_J);
return -1;
}
if (var_as_int(nc, VAR_K, &gg->k) == -1) {
varerr(nc, VAR_K);
return -1;
}
if (var_as_int(nc, VAR_M, &gg->m) == -1) {
varerr(nc, VAR_M);
return -1;
}
if (var_as_int(nc, VAR_TYPE, &gg->type) == -1) {
varerr(nc, VAR_TYPE);
return -1;
}
if (var_as_int(nc, VAR_MODE, &gg->mode) == -1) {
varerr(nc, VAR_MODE);
return -1;
}
if (nav->grid_type_code == GRID_RLL || nav->grid_type_code ==
GRID_SRLL) {
gg->rot = (rotated *)emalloc(sizeof(rotated));
if (var_as_float(nc, VAR_ROTLAT, &gg->rot->lat) == -1) {
varerr(nc, VAR_ROTLAT);
return -1;
}
if (var_as_float(nc, VAR_ROTLON, &gg->rot->lon) == -1) {
varerr(nc, VAR_ROTLON);
return -1;
}
if (var_as_float(nc, VAR_ROTANGLE, &gg->rot->angle) == -1) {
varerr(nc, VAR_ROTANGLE);
return -1;
}
}
if (nav->grid_type_code == GRID_SLL || nav->grid_type_code ==
GRID_SRLL) {
gg->strch = (stretched *)emalloc(sizeof(stretched));
if (var_as_float(nc, VAR_STRETCHLAT, &gg->strch->lat) == -1) {
varerr(nc, VAR_STRETCHLAT);
return -1;
}
if (var_as_float(nc, VAR_STRETCHLON, &gg->strch->lon) == -1) {
varerr(nc, VAR_STRETCHLON);
return -1;
}
if (var_as_float(nc, VAR_STRETCHFACTOR, &gg->strch->factor) == -1) {
varerr(nc, VAR_STRETCHFACTOR);
return -1;
}
}
}
break;
case GRID_MERCAT:
{
gdes_mercator *gg = &nav->grid.mercator;
if (var_as_int(nc, VAR_NI, &gg->ni) == -1) {
varerr(nc, VAR_NI);
return -1;
}
if (var_as_int(nc, VAR_NJ, &gg->nj) == -1) {
varerr(nc, VAR_NJ);
return -1;
}
if (var_as_float(nc, VAR_LA1, &gg->la1) == -1) {
varerr(nc, VAR_LA1);
return -1;
}
if (var_as_float(nc, VAR_LO1, &gg->lo1) == -1) {
varerr(nc, VAR_LO1);
return -1;
}
if (var_as_float(nc, VAR_LA2, &gg->la2) == -1) {
varerr(nc, VAR_LA2);
return -1;
}
if (var_as_float(nc, VAR_LO2, &gg->lo2) == -1) {
varerr(nc, VAR_LO2);
return -1;
}
if (var_as_float(nc, VAR_LATIN, &gg->latin) == -1) {
varerr(nc, VAR_LATIN);
return -1;
}
if (var_as_float(nc, VAR_DI, &gg->di) == -1) {
varerr(nc, VAR_DI);
return -1;
}
if (var_as_float(nc, VAR_DJ, &gg->dj) == -1) {
varerr(nc, VAR_DJ);
return -1;
}
}
break;
case GRID_GNOMON: /* fall through */
case GRID_POLARS:
{
gdes_polars *gg = &nav->grid.polars;
if (var_as_int(nc, VAR_NX, &gg->nx) == -1) {
varerr(nc, VAR_NX);
return -1;
}
if (var_as_int(nc, VAR_NY, &gg->ny) == -1) {
varerr(nc, VAR_NY);
return -1;
}
if (var_as_float(nc, VAR_LA1, &gg->la1) == -1) {
varerr(nc, VAR_LA1);
return -1;
}
if (var_as_float(nc, VAR_LO1, &gg->lo1) == -1) {
varerr(nc, VAR_LO1);
return -1;
}
if (var_as_float(nc, VAR_LOV, &gg->lov) == -1) {
varerr(nc, VAR_LOV);
return -1;
}
if (var_as_float(nc, VAR_DX, &gg->dx) == -1) {
varerr(nc, VAR_DX);
return -1;
}
if (var_as_float(nc, VAR_DY, &gg->dy) == -1) {
varerr(nc, VAR_DY);
return -1;
}
if (var_as_int(nc, VAR_PROJFLAG, &gg->pole) == -1) {
varerr(nc, VAR_PROJFLAG);
return -1;
}
}
break;
case GRID_LAMBERT:
{
gdes_lambert *gg = &nav->grid.lambert;
if (var_as_int(nc, VAR_NX, &gg->nx) == -1) {
varerr(nc, VAR_NX);
return -1;
}
if (var_as_int(nc, VAR_NY, &gg->ny) == -1) {
varerr(nc, VAR_NY);
return -1;
}
if (var_as_float(nc, VAR_LA1, &gg->la1) == -1) {
varerr(nc, VAR_LA1);
return -1;
}
if (var_as_float(nc, VAR_LO1, &gg->lo1) == -1) {
varerr(nc, VAR_LO1);
return -1;
}
if (var_as_float(nc, VAR_LOV, &gg->lov) == -1) {
varerr(nc, VAR_LOV);
return -1;
}
if (var_as_float(nc, VAR_DX, &gg->dx) == -1) {
varerr(nc, VAR_DX);
return -1;
}
if (var_as_float(nc, VAR_DY, &gg->dy) == -1) {
varerr(nc, VAR_DY);
return -1;
}
if (var_as_int(nc, VAR_PROJFLAG, &gg->pole) == -1) {
varerr(nc, VAR_PROJFLAG);
return -1;
}
if (var_as_float(nc, VAR_LATIN1, &gg->latin1) == -1) {
varerr(nc, VAR_LATIN1);
return -1;
}
if (var_as_float(nc, VAR_LATIN2, &gg->latin2) == -1) {
varerr(nc, VAR_LATIN2);
return -1;
}
if (var_as_float(nc, VAR_SPLAT, &gg->splat) == -1) {
varerr(nc, VAR_SPLAT);
return -1;
}
if (var_as_float(nc, VAR_SPLON, &gg->splon) == -1) {
varerr(nc, VAR_SPLON);
return -1;
}
}
break;
case GRID_SPACEV:
{
gdes_spacev *gg = &nav->grid.spacev;
if (var_as_int(nc, VAR_NX, &gg->nx) == -1) {
varerr(nc, VAR_NX);
return -1;
}
if (var_as_int(nc, VAR_NY, &gg->ny) == -1) {
varerr(nc, VAR_NY);
return -1;
}
if (var_as_float(nc, VAR_LAP, &gg->lap) == -1) {
varerr(nc, VAR_LAP);
return -1;
}
if (var_as_float(nc, VAR_LOP, &gg->lop) == -1) {
varerr(nc, VAR_LOP);
return -1;
}
if (var_as_float(nc, VAR_DX, &gg->dx) == -1) {
varerr(nc, VAR_DX);
return -1;
}
if (var_as_float(nc, VAR_DY, &gg->dy) == -1) {
varerr(nc, VAR_DY);
return -1;
}
if (var_as_float(nc, VAR_XP, &gg->xp) == -1) {
varerr(nc, VAR_XP);
return -1;
}
if (var_as_float(nc, VAR_YP, &gg->yp) == -1) {
varerr(nc, VAR_YP);
return -1;
}
if (var_as_float(nc, VAR_ORIENTATION, &gg->orient) == -1) {
varerr(nc, VAR_ORIENTATION);
return -1;
}
if (var_as_float(nc, VAR_NR, &gg->nr) == -1) {
varerr(nc, VAR_NR);
return -1;
}
if (var_as_float(nc, VAR_XO, &gg->xo) == -1) {
varerr(nc, VAR_XO);
return -1;
}
if (var_as_float(nc, VAR_YO, &gg->yo) == -1) {
varerr(nc, VAR_YO);
return -1;
}
}
break;
case GRID_ALBERS:
case GRID_OLAMBERT:
case GRID_UTM:
case GRID_SIMPOL:
case GRID_MILLER:
default:
uerror("%s: can't handle %s grids", nc->ncname,
gds_typename(nav->grid_type_code));
return -1;
}
return 0;
}
/*
* Make an in-memory structure for netCDF navigation information and
* initialize it from specified file. We cache this information because
* some of it must be consulted for every decoded product, and it won't
* change. Returns 0 on success, -1 on failure.
*/
static int
make_navinfo(nc, nav)
ncfile *nc; /* netCDF file */
navinfo *nav; /* where to put the navinfo */
{
char *cp;
int *ip;
float *fp;
nav->navid = nuwg_getdim(ncid, DIM_NAV);
if (var_as_int(nc, VAR_GRID_TYPE_CODE, &nav->grid_type_code) == -1) {
varerr(nc, VAR_GRID_TYPE_CODE);
return -1;
}
if (var_as_lset(nc, VAR_GRID_CENTER, &nav->grid_center) == -1) {
varerr(nc, VAR_GRID_CENTER);
return -1;
}
/* Multiple grid numbers allowed, they get stitched together */
if (var_as_lset(nc, VAR_GRID_NUMBER, &nav->grid_numbers) == -1) {
varerr(nc, VAR_GRID_NUMBER);
return -1;
}
/* GRIB resolution and component flags */
if (var_as_int(nc, VAR_RESCOMP, &nav->rescomp) == -1) {
varerr(nc, VAR_RESCOMP);
return -1;
}
if (make_navgrid(nc, nav) == -1) {
return -1;
}
return 0;
}
static void
free_navinfo(np)
navinfo* np;
{
if (np) {
if(np->nav_model)
free(np->nav_model);
if(np->grid_numbers.vals)
free(np->grid_numbers.vals);
free(np);
}
}
static navinfo *
new_navinfo(nc)
ncfile *nc;
{
navinfo *out = (navinfo *)emalloc(sizeof(navinfo));
if (make_navinfo(nc, out) != 0) {
free_navinfo(out);
return 0;
}
return out;
}
/*
* Make an in-memory structure for netCDF information and initialize it from
* specified file. We cache this information because some of it must be
* consulted for every decoded product, and it won't change.
*/
static int
make_ncfile(ncname, out)
char *ncname;
ncfile *out;
{
int ndims; /* number of dimensions */
int nvars; /* number of variables */
int recid;
long nrecs; /* number of records */
int dimid;
int varid;
int *ip;
out->ncname = estrdup(ncname);
out->ncid = ncid;
if (ncinquire(ncid, &ndims, &nvars, (int *)0, &recid) == -1) {
uerror("%s: ncinquire() failed", ncname);
return -1;
}
out->ndims = ndims;
out->nvars = nvars;
out->dims = (ncdim **)emalloc(ndims * sizeof(ncdim *));
for (dimid = 0; dimid < ndims; dimid++) {
out->dims[dimid] = new_dim(dimid);
}
out->vars = (ncvar **)emalloc(nvars * sizeof(ncvar *));
for (varid = 0; varid < nvars; varid++) {
out->vars[varid] = new_var(ncname, varid);
}
if (recid == -1) {
uerror("%s: no record dimension", ncname);
return -1;
}
out->recid = recid;
out->reftimeid = nuwg_getvar(ncid, VAR_REFTIME);
out->valtimeid = nuwg_getvar(ncid, VAR_VALTIME);
if(out->reftimeid == -1) {
uerror("%s: no reftime variable", ncname);
return -1;
}
if(out->valtimeid == -1) {
uerror("%s: no valtime variable", ncname);
return -1;
}
out->datetimeid = nuwg_getvar(ncid, VAR_DATETIME);
out->valoffsetid = nuwg_getvar(ncid, VAR_VALOFFSET);
if(out->datetimeid == -1) {
uerror("%s: no datetimeid variable", ncname);
return -1;
}
if(out->valoffsetid == -1) {
uerror("%s: no valoffsetid variable", ncname);
return -1;
}
/* initialize reftime,valtime,record table */
if (new_recs(out) == -1) {
uerror("%s: can't initialize reftime,valtime table", ncname);
return -1;
}
/* Multiple model numbers allowed, e.g. for initialization */
if (var_as_lset(out, VAR_MODELID, &out->models) == -1) {
varerr(out, VAR_MODELID);
return -1;
}
out->nav = new_navinfo(out);
if (!out->nav) { /* get navigation information */
uerror("%s: can't get navigation information", out->ncname);
return -1;
}
out->levdims = 0; /* only add level dimensions as needed */
out->laydims = 0; /* only add layer dimensions as needed */
return 0;
}
static void
free_ncfile(np)
ncfile* np;
{
if (np) {
if(np->ncname)
free(np->ncname);
if(np->models.vals)
free(np->models.vals);
if(np->dims) {
int dimid;
for (dimid = 0; dimid < np->ndims; dimid++) {
free_dim(np->dims[dimid]);
}
free(np->dims);
}
if(np->vars) {
int varid;
for (varid = 0; varid < np->nvars; varid++) {
free_var(np->vars[varid]);
}
free(np->vars);
}
if(np->rt)
free_recs(np->rt);
free(np);
if(np->nav)
free_navinfo(np->nav);
if(np->levdims)
free(np->levdims);
if(np->laydims)
free(np->laydims);
}
}
/*
* Creates a new ncfile structure and fills it in with the information from
* the open netCDF file whose handle is ncid. A pointer to the structure is
* returned, or 0 on failure.
*/
ncfile *
new_ncfile(ncname)
char *ncname;
{
ncfile *out = (ncfile *)emalloc(sizeof(ncfile));
if (make_ncfile(ncname, out) != 0) {
uerror("make_ncfile failed");
free_ncfile(out);
return 0;
}
return out;
}
/*
* Print an error message if two specified floating-point values are not
* sufficiently close. Returns 1 if not close, zero otherwise.
*/
static int
gd_fne_err(pp, nc, comp, pval, nval)
product_data *pp; /* decoded GRIB product */
ncfile *nc; /* netCDF file */
enum ncpart comp; /* name-independent id of component */
float pval; /* first value, from GRIB product */
float nval; /* second value, from netCDF file */
{
#define float_near(x,y) ((float)((y) + 0.1*fabs((x)-(y))) == (float)(y)) /*
true if x is "close to" y */
if (! float_near(pval, nval)) {
uerror("%s, %s nav. mismatch %s: %g != %g",
pp->header, nc->ncname, nuwg_name(comp), pval, nval);
return 1;
}
return 0;
}
/*
* Print an error message if two specified floating-point values are not
* integer close. Returns 1 if not close, zero otherwise.
*/
static int
gd_fnei_err(pp, nc, comp, pval, nval)
product_data *pp; /* decoded GRIB product */
ncfile *nc; /* netCDF file */
enum ncpart comp; /* name-independent id of component */
float pval; /* first value, from GRIB product */
float nval; /* second value, from netCDF file */
{
if ((int) (pval+0.5) != (int) (nval+0.5)) {
uerror("%s, %s nav. mismatch %s: %d != %d",
pp->header, nc->ncname, nuwg_name(comp), (int) (pval+0.5),
(int) (nval+0.5));
return 1;
}
return 0;
}
/*
* Print an error message if two specified int values are not
* equal. Return 1 if not equal, zero otherwise.
*/
static int
gd_ine_err(pp, nc, comp, pval, nval)
product_data *pp; /* decoded GRIB product */
ncfile *nc; /* netCDF file */
enum ncpart comp; /* name-independent id of component */
int pval; /* first value, from GRIB product */
int nval; /* second value, from netCDF file */
{
if (pval != nval) {
uerror("%s, %s nav. mismatch %s: %d != %d",
pp->header, nc->ncname, nuwg_name(comp), pval, nval);
return 1;
}
return 0;
}
/*
* Print an error message if first int value is greater than second.
* Return 1 if not equal, zero otherwise.
*/
static int
gd_igt_err(pp, nc, comp, pval, nval)
product_data *pp; /* decoded GRIB product */
ncfile *nc; /* netCDF file */
enum ncpart comp; /* name-independent id of component */
int pval; /* first value, from GRIB product */
int nval; /* second value, from netCDF file */
{
if (pval > nval) {
uerror("%s, %s nav. mismatch %s: %d > %d",
pp->header, nc->ncname, nuwg_name(comp), pval, nval);
return 1;
}
return 0;
}
/*
* Check consistency of decoded GRIB product grid description section
* information with the navigation information in the netCDF file. Return
* -1 if inconsistency found, 0 otherwise.
*/
static int
check_gd(pp, nc)
product_data *pp; /* decoded GRIB product */
ncfile *nc; /* netCDF file */
{
gengrid *gp = &pp->gd->grid; /* decoded Grid Description Section */
navinfo *nav = nc->nav;
gengrid *gn; /* Grid Description Section info in netCDF */
int errs = 0;
if (gp == 0) {
uerror("%s: no grid description in product?", pp->header);
return -1;
}
if (nav == 0) {
uerror("%s: no navigation info in netCDF file?", nc->ncname);
return -1;
}
gn = &nc->nav->grid; /* GDS info in netCDF file */
switch(nav->grid_type_code) {
case GRID_LL:
case GRID_RLL:
case GRID_SLL:
case GRID_SRLL:
{
gdes_ll *gng = &gn->ll;
gdes_ll *gpg = &gp->ll;
errs += gd_igt_err(pp, nc, VAR_NI, gpg->ni, gng->ni);
errs += gd_igt_err(pp, nc, VAR_NJ, gpg->nj, gng->nj);
errs += gd_fne_err(pp, nc, VAR_DI, gpg->di, gng->di);
errs += gd_fne_err(pp, nc, VAR_DJ, gpg->dj, gng->dj);
if (nav->grid_type_code == GRID_RLL || nav->grid_type_code ==
GRID_SRLL) {
if (gpg->rot == 0) {
uerror("%s: rotated grid, but no rotation info?",
pp->header);
return -1;
}
if (gng->rot == 0) {
uerror("%s: no rotation info for rotated grid in netCDF",
nc->ncname);
return -1;
}
errs += gd_fne_err(pp, nc, VAR_ROTLAT, gpg->rot->lat,
gng->rot->lat);
errs += gd_fne_err(pp, nc, VAR_ROTLON, gpg->rot->lon,
gng->rot->lon);
errs += gd_fne_err(pp, nc, VAR_ROTANGLE, gpg->rot->angle,
gng->rot->angle);
}
if (nav->grid_type_code == GRID_SLL || nav->grid_type_code ==
GRID_SRLL) {
if (gpg->strch == 0) {
uerror("%s: stretched grid, but no stretch info?",
pp->header);
return -1;
}
if (gng->strch == 0) {
uerror("%s: no stretch info for stretched grid in netCDF",
nc->ncname);
return -1;
}
errs += gd_fne_err(pp, nc, VAR_STRETCHLAT, gpg->strch->lat,
gng->strch->lat);
errs += gd_fne_err(pp, nc, VAR_STRETCHLON, gpg->strch->lon,
gng->strch->lon);
errs += gd_fne_err(pp, nc, VAR_STRETCHFACTOR, gpg->strch->factor,
gng->strch->factor);
}
}
break;
case GRID_GAU:
case GRID_RGAU:
case GRID_SGAU:
case GRID_SRGAU:
{
gdes_gau *gng = &gn->gau;
gdes_gau *gpg = &gp->gau;
errs += gd_ine_err(pp, nc, VAR_NI, gpg->ni, gng->ni);
errs += gd_ine_err(pp, nc, VAR_NJ, gpg->nj, gng->nj);
errs += gd_ine_err(pp, nc, VAR_N, gpg->n, gng->n);
errs += gd_fne_err(pp, nc, VAR_LA1, gpg->la1, gng->la1);
errs += gd_fne_err(pp, nc, VAR_LO1, gpg->lo1, gng->lo1);
errs += gd_fne_err(pp, nc, VAR_LA2, gpg->la2, gng->la2);
errs += gd_fne_err(pp, nc, VAR_LO2, gpg->lo2, gng->lo2);
errs += gd_fne_err(pp, nc, VAR_DI, gpg->di, gng->di);
if (nav->grid_type_code == GRID_RLL || nav->grid_type_code ==
GRID_SRLL) {
if (gpg->rot == 0) {
uerror("%s: rotated grid, but no rotation info?",
pp->header);
return -1;
}
if (gng->rot == 0) {
uerror("%s: no rotation info for rotated grid in netCDF",
nc->ncname);
return -1;
}
errs += gd_fne_err(pp, nc, VAR_ROTLAT, gpg->rot->lat,
gng->rot->lat);
errs += gd_fne_err(pp, nc, VAR_ROTLON, gpg->rot->lon,
gng->rot->lon);
errs += gd_fne_err(pp, nc, VAR_ROTANGLE, gpg->rot->angle,
gng->rot->angle);
}
if (nav->grid_type_code == GRID_SLL || nav->grid_type_code ==
GRID_SRLL) {
if (gpg->strch == 0) {
uerror("%s: stretched grid, but no stretch info?",
pp->header);
return -1;
}
if (gng->strch == 0) {
uerror("%s: no stretch info for stretched grid in netCDF",
nc->ncname);
return -1;
}
errs += gd_fne_err(pp, nc, VAR_STRETCHLAT, gpg->strch->lat,
gng->strch->lat);
errs += gd_fne_err(pp, nc, VAR_STRETCHLON, gpg->strch->lon,
gng->strch->lon);
errs += gd_fne_err(pp, nc, VAR_STRETCHFACTOR, gpg->strch->factor,
gng->strch->factor);
}
}
break;
case GRID_SPH:
case GRID_RSPH:
case GRID_SSPH:
case GRID_SRSPH:
{
gdes_sph *gng = &gn->sph;
gdes_sph *gpg = &gp->sph;
errs += gd_ine_err(pp, nc, VAR_J, gpg->j, gng->j);
errs += gd_ine_err(pp, nc, VAR_K, gpg->k, gng->k);
errs += gd_ine_err(pp, nc, VAR_M, gpg->m, gng->m);
errs += gd_ine_err(pp, nc, VAR_TYPE, gpg->type, gng->type);
errs += gd_ine_err(pp, nc, VAR_MODE, gpg->mode, gng->mode);
if (nav->grid_type_code == GRID_RLL || nav->grid_type_code ==
GRID_SRLL) {
if (gpg->rot == 0) {
uerror("%s: rotated grid, but no rotation info?",
pp->header);
return -1;
}
if (gng->rot == 0) {
uerror("%s: no rotation info for rotated grid in netCDF",
nc->ncname);
return -1;
}
errs += gd_fne_err(pp, nc, VAR_ROTLAT, gpg->rot->lat,
gng->rot->lat);
errs += gd_fne_err(pp, nc, VAR_ROTLON, gpg->rot->lon,
gng->rot->lon);
errs += gd_fne_err(pp, nc, VAR_ROTANGLE, gpg->rot->angle,
gng->rot->angle);
}
if (nav->grid_type_code == GRID_SLL || nav->grid_type_code ==
GRID_SRLL) {
if (gpg->strch == 0) {
uerror("%s: stretched grid, but no stretch info?",
pp->header);
return -1;
}
if (gng->strch == 0) {
uerror("%s: no stretch info for stretched grid in netCDF",
nc->ncname);
return -1;
}
errs += gd_fne_err(pp, nc, VAR_STRETCHLAT, gpg->strch->lat,
gng->strch->lat);
errs += gd_fne_err(pp, nc, VAR_STRETCHLON, gpg->strch->lon,
gng->strch->lon);
errs += gd_fne_err(pp, nc, VAR_STRETCHFACTOR, gpg->strch->factor,
gng->strch->factor);
}
}
break;
case GRID_MERCAT:
{
gdes_mercator *gng = &gn->mercator;
gdes_mercator *gpg = &gp->mercator;
errs += gd_ine_err(pp, nc, VAR_NI, gpg->ni, gng->ni);
errs += gd_ine_err(pp, nc, VAR_NJ, gpg->nj, gng->nj);
errs += gd_fne_err(pp, nc, VAR_LA1, gpg->la1, gng->la1);
errs += gd_fne_err(pp, nc, VAR_LO1, gpg->lo1, gng->lo1);
errs += gd_fne_err(pp, nc, VAR_LA2, gpg->la2, gng->la2);
errs += gd_fne_err(pp, nc, VAR_LO2, gpg->lo2, gng->lo2);
errs += gd_fne_err(pp, nc, VAR_LATIN, gpg->latin, gng->latin);
errs += gd_fne_err(pp, nc, VAR_DI, gpg->di, gng->di);
errs += gd_fne_err(pp, nc, VAR_DJ, gpg->dj, gng->dj);
}
break;
case GRID_GNOMON: /* fall through */
case GRID_POLARS:
{
gdes_polars *gng = &gn->polars;
gdes_polars *gpg = &gp->polars;
errs += gd_ine_err(pp, nc, VAR_NX, gpg->nx, gng->nx);
errs += gd_ine_err(pp, nc, VAR_NY, gpg->ny, gng->ny);
errs += gd_fne_err(pp, nc, VAR_LA1, gpg->la1, gng->la1);
errs += gd_fne_err(pp, nc, VAR_LO1, gpg->lo1, gng->lo1);
errs += gd_fne_err(pp, nc, VAR_LOV, gpg->lov, gng->lov);
errs += gd_fnei_err(pp, nc, VAR_DX, gpg->dx, gng->dx);
errs += gd_fnei_err(pp, nc, VAR_DY, gpg->dy, gng->dy);
errs += gd_ine_err(pp, nc, VAR_PROJFLAG, gpg->pole, gng->pole);
}
break;
case GRID_LAMBERT:
{
gdes_lambert *gng = &gn->lambert;
gdes_lambert *gpg = &gp->lambert;
errs += gd_ine_err(pp, nc, VAR_NX, gpg->nx, gng->nx);
errs += gd_ine_err(pp, nc, VAR_NY, gpg->ny, gng->ny);
errs += gd_fne_err(pp, nc, VAR_LA1, gpg->la1, gng->la1);
errs += gd_fne_err(pp, nc, VAR_LO1, gpg->lo1, gng->lo1);
errs += gd_fne_err(pp, nc, VAR_LOV, gpg->lov, gng->lov);
errs += gd_fnei_err(pp, nc, VAR_DX, gpg->dx, gng->dx);
errs += gd_fnei_err(pp, nc, VAR_DY, gpg->dy, gng->dy);
errs += gd_ine_err(pp, nc, VAR_PROJFLAG, gpg->pole, gng->pole);
errs += gd_fne_err(pp, nc, VAR_LATIN1, gpg->latin1, gng->latin1);
errs += gd_fne_err(pp, nc, VAR_LATIN2, gpg->latin2, gng->latin2);
errs += gd_fne_err(pp, nc, VAR_SPLAT, gpg->splat, gng->splat);
errs += gd_fne_err(pp, nc, VAR_SPLON, gpg->splon, gng->splon);
}
break;
case GRID_SPACEV:
{
gdes_spacev *gng = &gn->spacev;
gdes_spacev *gpg = &gp->spacev;
errs += gd_ine_err(pp, nc, VAR_NX, gpg->nx, gng->nx);
errs += gd_ine_err(pp, nc, VAR_NY, gpg->ny, gng->ny);
errs += gd_fne_err(pp, nc, VAR_LAP, gpg->lap, gng->lap);
errs += gd_fne_err(pp, nc, VAR_LOP, gpg->lop, gng->lop);
errs += gd_fnei_err(pp, nc, VAR_DX, gpg->dx, gng->dx);
errs += gd_fnei_err(pp, nc, VAR_DY, gpg->dy, gng->dy);
errs += gd_fne_err(pp, nc, VAR_XP, gpg->xp, gng->xp);
errs += gd_fne_err(pp, nc, VAR_YP, gpg->yp, gng->yp);
errs += gd_fne_err(pp, nc, VAR_ORIENTATION, gpg->orient, gng->orient);
errs += gd_ine_err(pp, nc, VAR_NR, gpg->nr, gng->nr);
errs += gd_fne_err(pp, nc, VAR_XO, gpg->xo, gng->xo);
errs += gd_fne_err(pp, nc, VAR_YO, gpg->yo, gng->yo);
}
break;
case GRID_ALBERS:
case GRID_OLAMBERT:
case GRID_UTM:
case GRID_SIMPOL:
case GRID_MILLER:
default:
uerror("%s: can't handle %s grids", nc->ncname,
gds_typename(nav->grid_type_code));
return -1;
}
if (errs > 0)
return -1;
return 0;
}
/*
* Consistency check between decoded product data and netCDF file to which
* it will be written. Returns -1 if a serious inconsistency was found.
*
* If a value is missing in the netCDF file, it should get filled in here.
* ( *** not implemented yet *** ).
*
* The netCDF information is read from cached info rather than from the file
* for speed, since this check takes place on every product.
*/
static int
check_pp(pp, nc)
product_data *pp; /* decoded GRIB product */
ncfile *nc; /* netCDF file */
{
navinfo* nav = nc->nav;
int found;
int i;
found = 0;
for (i = 0; i < nc->models.n; i++) {
if (pp->model == nc->models.vals[i]) {
found = 1;
break;
}
}
if (!found) {
uerror("%s model %d not in modelid list in %s", pp->header,
pp->model, nc->ncname);
return -1;
}
found = 0;
for (i = 0; i < nav->grid_center.n; i++) {
if (pp->center == nav->grid_center.vals[i]) {
found = 1;
break;
}
}
if (!found) {
uerror("%s center %d not in grid_center list in %s", pp->header,
pp->center, nc->ncname);
return -1;
}
found = 0;
for (i = 0; i < nav->grid_numbers.n; i++) {
if (pp->grid == nav->grid_numbers.vals[i]) {
found = 1;
break;
}
}
if (!found) {
uerror("%s grid %d not in grid_number list in %s", pp->header,
pp->grid, nc->ncname);
return -1;
}
/* gd->res_flags should match what's in nc, but if not, rewrite it. */
if (pp->gd->res_flags != nav->rescomp) {
udebug("%s: rescompflag in %s doesn't match, %d, %d",
pp->header, nc->ncname, (int) pp->gd->res_flags, nav->rescomp);
}
if(check_gd(pp, nc) == -1) {
return -1;
}
return 0;
}
/*
* Figure out subgrid location from netCDF navigation information and
* product Grid Description Section information. Returns 0 on success, -1
* on failure.
*/
static int
subgrid(nc, pp, ix0, ix1)
ncfile *nc;
product_data *pp;
long *ix0;
long *ix1;
{
navinfo* nav = nc->nav;
gdes_ll *gll;
float plon;
*ix0 = 0;
*ix1 = 0;
switch (pp->gd->type) {
case GRID_LL:
case GRID_RLL:
case GRID_SLL:
case GRID_SRLL:
if(pp->gd->scan_mode & 0x20 == 1 || /*adjacent points in j direction */
pp->gd->scan_mode & 0x80 == 1 ) { /* points scan in -i direction */
uerror("%s: can't handle scan mode of %x ",
pp->header,pp->gd->scan_mode);
return -1;
}
/* If scanning mode flag indicates North to South scan, we
reverse the rows so we can always assume South to North scan */
if( (pp->gd->scan_mode != 0) &&
(pp->gd->scan_mode & 0x40) == 0) { /* north to south */
int row;
int nrows = pp->gd->nrows;
int ncols = pp->gd->ncols;
float tmp;
for (row = 0; row < nrows/2; row++) {
int col;
float *upper = pp->data + row * ncols;
float *lower = pp->data + (nrows - 1 - row) * ncols;
/* swap row (upper) and nrows-1-row (lower) */
for (col = 0; col < ncols; col++) {
tmp = *upper;
*upper = *lower;
*lower = tmp;
upper++;
lower++;
}
}
gll = &pp->gd->grid.ll;
/*
uerror( "gll->la1=%f, gll->la2=%f,gll->lo1=%f, gll->lo2=%f",
gll->la1, gll->la2, gll->lo1, gll->lo2 );
*/
tmp = gll->la1;
gll->la1 = gll->la2;
gll->la2 = tmp;
pp->gd->scan_mode |= 0x40;
}
/* Compare pp->gdes->grid.ll->la1 to value of La1, La2 netCDF
variables and similarly for lo1, lo2. */
gll = &pp->gd->grid.ll;
/*
uerror( "gll->la1=%f, gll->la2=%f,gll->lo1=%f, gll->lo2=%f",
gll->la1, gll->la2, gll->lo1, gll->lo2 );
*/
plon = gll->lo1;
while (plon < nc->nav->grid.ll.lo1)
plon += 360.0 ;
/* handle case where right edge is same as left edge */
if (plon == nc->nav->grid.ll.lo2 &&
nc->nav->grid.ll.lo1 + 360 == nc->nav->grid.ll.lo2) {
plon = nc->nav->grid.ll.lo1;
}
/*
uerror("gll->la1=%f, nc->nav->grid.ll.la1=%f, nc->nav->grid.ll.dj=%f",
gll->la1, nc->nav->grid.ll.la1, nc->nav->grid.ll.dj);
*/
*ix0 = (gll->la1 - nc->nav->grid.ll.la1) /nc->nav->grid.ll.dj + 0.5;
*ix1 = (plon - nc->nav->grid.ll.lo1) /nc->nav->grid.ll.di + 0.5;
/*
uerror("*ix0=%f, *ix1=%f", *ix0, *ix1);
*/
}
/* default: */
return 0;
}
/*
* Handle writing extra time-range indicator information, if any, in
* auxilliary variables. For example, accumulation interval for
* precipitation variables would be handled here.
* For now, we do this in an ad hoc way, until we have a mechanism for
* writing general auxilliary GRIB info associated with a variable.
*/
static int
triaux(pp, nc, var, start)
product_data *pp; /* decoded GRIB data to be written */
ncfile *nc; /* netCDF file to be written */
ncvar *var; /* netCDF variable to be written */
long *start; /* index where variable to be written */
{
char tri_name[MAX_NC_NAME];
char *suf;
ncvar *trivar;
int trivarid;
long ix[2];
long count[2];
float trivals[2];
if(pp->tr_flg == TRI_P1 || pp->tr_flg == TRI_LP1)
return 0; /* usually time range info is in valid_time variable */
/* check if auxilliary time-range variable exists */
suf = trisuffix(pp->tr_flg);
if (strlen(var->name) + 1 + strlen(suf) > (size_t) MAX_NC_NAME) {
uerror("%s: name of %s TRI variable too long (%s)",
nc->ncname, suf, var->name);
return -1;
}
strcpy(tri_name, var->name);
strcat(tri_name, "_");
strcat(tri_name, suf);
trivarid = ncvarid(ncid, tri_name);
if(trivarid == -1) {
return 0; /* not an error, since optional whether file
* has auxilliary time-range variable */
}
/* *** should check units of _accum_len variable and convert to
those, use float_nc() *** */
trivar = nc->vars[trivarid];
switch (trinum(pp->tr_flg)) {
case 2:
ix[0]=start[0]; /* record number */
ix[1]=0;
count[0] = 1;
count[1] = 2; /* two values for time range */
trivals[0] = pp->tr[0];
trivals[1] = pp->tr[1];
if (ncvarput(ncid, trivarid, ix, count, &trivals) == -1) {
uerror("%s: can't write accum_len variable for (%s)",
nc->ncname, var->name);
return -1;
}
break;
case 1:
ix[0]=start[0];
trivals[0] = frcst_time(pp);
if (ncvarput1(ncid, trivarid, ix, &trivals) == -1) {
uerror("%s: can't write accum_len variable for (%s)",
nc->ncname, var->name);
return -1;
}
break;
default:
uerror("%s: can't handle time flag %d for variable (%s)",
nc->ncname, pp->tr_flg, var->name);
return -1;
}
return 0;
}
/*
* Writes decoded GRIB product to netCDF file open for writing.
* Returns 0 on success, -1 on failure.
*/
int
nc_write(pp, nc)
product_data *pp; /* decoded GRIB product to be written */
ncfile *nc; /* netCDF file to write */
{
double reftime, valtime;
long rec;
char *varname;
ncvar *var;
int varid;
long lev; /* level to write, if level dimension */
char *cp = parmname(nc, pp->param, pp->level_flg); /* var to write */
int dim=0; /* which dimension we are dealing with */
#define MAX_PARM_DIMS 4 /* Maximum dimensions for a parameter.
X(rec,lev,lat,lon) This would be 4 for a
variable with a level dimension, 3 if no
level dimension. */
long start[MAX_PARM_DIMS];
long count[MAX_PARM_DIMS];
/* Check consistency of product with
* netCDF file, e.g. model_id, gdes vs.
* netCDF navigation */
if(check_pp(pp, nc) == -1) {
/* any messages about inconsistency should have already been logged */
return -1;
}
if (!cp) {
uerror("%s: unrecognized (param,level_flg) combination (%d,%d)",
pp->header, pp->param, pp->level_flg);
return -1;
}
varname = estrdup(cp);
/* locate variable in output netCDF file */
varid = ncvarid(ncid,varname);
if (varid == -1) {
uerror("%s: no variable %s in %s", pp->header, varname,
nc->ncname);
return -1;
}
var = nc->vars[varid];
if (var->dims[0] != nc->recid) { /* no record dimension */
dim--;
} else { /* handle record dimension */
humtime ht;
if (rvhours(pp, nc, &reftime, &valtime, &ht) != 0) {
uerror("%s: can't get reftime,valtime", pp->header);
return -1;
}
rec = getrec(nc, reftime, valtime, &ht); /* which record to write */
if (rec == -1) {
uerror("%s: can't get record number for variable %s in %s",
pp->header, varname, nc->ncname);
return -1;
}
start[dim] = rec;
count[dim] = 1;
}
/* Handle auxilliary time-range indicator information, if any */
if (triaux(pp, nc, var, start) == -1)
return -1;
/* handle level dimension, if any */
lev = getlev(pp, nc, var);
if (lev == -1)
return -1;
if(lev >= 0) {
dim++;
start[dim] = lev;
count[dim] = 1;
}
/* Write varname(rec,lev) or varname(rec) hyperslab, but start in the
right place if this is only a subgrid of the netCDF file domain */
{
double slope, intercept;
long ix0; /* if subgridded, lat index, otherwise 0 */
long ix1; /* if subgridded, lon index, otherwise 0 */
if (subgrid(nc, pp, &ix0, &ix1) == -1) {
uerror("%s: can't fit %s variable subgrid into %s",
pp->header, varname, nc->ncname);
return -1;
}
dim++;
start[dim] = ix0;
count[dim] = pp->gd->nrows;
dim++;
start[dim] = ix1;
count[dim] = pp->gd->ncols;
if (var->uc) { /* units conversion */
slope = var->uc->slope;
intercept = var->uc->intercept;
} else {
slope = 1.0;
intercept = 0.0;
}
if (float_nc(ncid, varid, start, count,
pp->data, slope, intercept, FILL_VAL) == -1) {
uerror("%s: error writing %s variable in %s",
pp->header, varname, nc->ncname);
return -1;
}
}
if (ncsync(ncid) == -1) {
uerror("%s: ncsync failed after writing %s variable in %s",
pp->header, varname, nc->ncname);
return -1;
}
if (var->dims[0] != nc->recid) { /* no record dimension */
if (lev < 0) {
uinfo("%s: wrote %s(*,*) to %s", pp->header, varname,
nc->ncname);
} else if (lev >= 0) {
uinfo("%s: wrote %s(%ld,*,*) to %s", pp->header, varname,
lev, nc->ncname);
}
} else {
if (lev < 0) {
uinfo("%s: wrote %s(%ld,*,*) to %s", pp->header, varname,
rec, nc->ncname);
} else if (lev >= 0) {
uinfo("%s: wrote %s(%ld,%ld,*,*) to %s", pp->header, varname,
rec, lev, nc->ncname);
}
}
free(varname);
return 0;
}
