This archive contains answers to questions sent to Unidata support through mid-2025. Note that the archive is no longer being updated. We provide the archive for reference; many of the answers presented here remain technically correct, even if somewhat outdated. For the most up-to-date information on the use of NSF Unidata software and data services, please consult the Software Documentation first.
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; }