Skip to content

GitLab

Commit dd50a8aa authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Added read support for hybrid sigma pressure coordinate with formula term P0.

parent be897882
Hide whitespace changes
Inline Side-by-side
Showing with 59 additions and 47 deletions
ChangeLog
View file @ dd50a8aa
2016-09-26 Uwe Schulzweida
* Added read support for hybrid sigma pressure coordinate with formula term P0
2016-09-22 Uwe Schulzweida
* GRIB: sort only unsorted pressure levels
......
src/stream_cdf_i.c
View file @ dd50a8aa
......@@ -83,9 +83,9 @@ typedef struct {
bool defmissval;
bool deffillval;
int xtype;
int ndims;
int gmapid;
int positive;
int ndims;
int dimids[8];
int dimtype[8];
int chunks[8];
......@@ -732,13 +732,11 @@ void cdf_set_dim(ncvar_t *ncvars, int ncvarid, int dimid, int dimtype)
}
static
int scan_hybrid_formulaterms(int ncid, int ncfvarid, int *apvarid, int *bvarid, int *psvarid, int *avarid, int *p0varid)
void scan_hybrid_formulaterms(int ncid, int ncfvarid, int *avarid, int *bvarid, int *psvarid, int *p0varid)
{
int status = 1;
*apvarid = -1;
*avarid = -1;
*bvarid = -1;
*psvarid = -1;
*avarid = -1;
*p0varid = -1;
char attstring[1024];
......@@ -766,7 +764,8 @@ int scan_hybrid_formulaterms(int ncid, int ncfvarid, int *apvarid, int *bvarid,
int status_nc = nc_inq_varid(ncid, varname, &dimvarid);
if ( status_nc == NC_NOERR )
{
if ( strcmp(tagname, "ap:") == 0 ) *apvarid = dimvarid;
if ( strcmp(tagname, "ap:") == 0 ) *avarid = dimvarid;
else if ( strcmp(tagname, "a:") == 0 ) *avarid = dimvarid;
else if ( strcmp(tagname, "b:") == 0 ) *bvarid = dimvarid;
else if ( strcmp(tagname, "ps:") == 0 ) *psvarid = dimvarid;
else if ( strcmp(tagname, "p0:") == 0 ) *p0varid = dimvarid;
......@@ -778,8 +777,6 @@ int scan_hybrid_formulaterms(int ncid, int ncfvarid, int *apvarid, int *bvarid,
if ( lstop ) break;
}
return status;
}
static
......@@ -794,59 +791,68 @@ bool isHybridSigmaPressureCoordinate(int ncid, int ncvarid, ncvar_t *ncvars, con
status = true;
ncvar->zaxistype = ZAXIS_HYBRID;
//int ndims = ncvar->ndims;
int dimid = ncvar->dimids[0];
size_t dimlen = ncdims[dimid].len;
int ret = 0;
int apvarid1 = -1, bvarid1 = -1, psvarid1 = -1, avarid1 = -1, p0varid1 = -1;
int avarid1 = -1, bvarid1 = -1, psvarid1 = -1, p0varid1 = -1;
int ncfvarid = ncvarid;
if ( ncvars[ncfvarid].lformulaterms )
ret = scan_hybrid_formulaterms(ncid, ncfvarid, &apvarid1, &bvarid1, &psvarid1, &avarid1, &p0varid1);
// printf("ret %d apvarid1, bvarid1, psvarid1, avarid1, p0varid1 %d %d %d %d %d\n", ret, apvarid1, bvarid1, psvarid1, avarid1, p0varid1);
if ( apvarid1 != -1 ) ncvars[apvarid1].isvar = FALSE;
if ( bvarid1 != -1 ) ncvars[bvarid1].isvar = FALSE;
if ( psvarid1 != -1 ) ncvar->psvarid = psvarid1;
scan_hybrid_formulaterms(ncid, ncfvarid, &avarid1, &bvarid1, &psvarid1, &p0varid1);
// printf("avarid1, bvarid1, psvarid1, p0varid1 %d %d %d %d\n", avarid1, bvarid1, psvarid1, p0varid1);
if ( avarid1 != -1 ) ncvars[avarid1].isvar = FALSE;
if ( bvarid1 != -1 ) ncvars[bvarid1].isvar = FALSE;
if ( psvarid1 != -1 ) ncvar->psvarid = psvarid1;
if ( p0varid1 != -1 ) ncvar->p0varid = p0varid1;
if ( ncvar->bounds != UNDEFID && ncvars[ncvar->bounds].lformulaterms )
{
ncfvarid = ncvar->bounds;
int apvarid2 = -1, bvarid2 = -1, psvarid2 = -1, avarid2 = -1, p0varid2 = -1;
ret = 0;
int avarid2 = -1, bvarid2 = -1, psvarid2 = -1, p0varid2 = -1;
if ( ncvars[ncfvarid].lformulaterms )
ret = scan_hybrid_formulaterms(ncid, ncfvarid, &apvarid2, &bvarid2, &psvarid2, &avarid2, &p0varid2);
// printf("ret %d apvarid2, bvarid2, psvarid2, avarid2, p0varid2 %d %d %d %d %d\n", ret, apvarid2, bvarid2, psvarid2, avarid2, p0varid2);
if ( ret == 1 ) avarid2 = apvarid2;
scan_hybrid_formulaterms(ncid, ncfvarid, &avarid2, &bvarid2, &psvarid2, &p0varid2);
// printf("avarid2, bvarid2, psvarid2, p0varid2 %d %d %d %d\n", avarid2, bvarid2, psvarid2, p0varid2);
if ( avarid2 != -1 && bvarid2 != -1 )
{
ncvars[avarid2].isvar = FALSE;
ncvars[bvarid2].isvar = FALSE;
if ( dimid == ncvars[avarid2].dimids[0] && ncdims[ncvars[avarid2].dimids[1]].len == 2 )
int ndims2 = ncvars[avarid2].ndims;
int dimid2 = ncvars[avarid2].dimids[0];
size_t dimlen2 = ncdims[dimid2].len;
if ( (ndims2 == 2 && dimid == ncvars[avarid2].dimids[0] ) ||
(ndims2 == 1 && dimlen == dimlen2-1 ) )
{
double px = 1;
if ( ret == 2 && p0varid1 == p0varid2 )
if ( p0varid1 != -1 && p0varid1 == p0varid2 )
cdf_get_var_double(ncid, p0varid2, &px);
double abuf[dimlen*2], bbuf[dimlen*2];
cdf_get_var_double(ncid, avarid2, abuf);
cdf_get_var_double(ncid, bvarid2, bbuf);
/*
for ( int i = 0; i < dimlen; ++i )
printf("%d %g %g %g %g\n", i, abuf[i*2], abuf[i*2+1], bbuf[i*2], bbuf[i*2+1]);
*/
size_t vctsize = (dimlen+1)*2;
double *vct = (double *) Malloc(vctsize * sizeof(double));
for ( size_t i = 0; i < dimlen; ++i )
double *vct = (double *) Malloc(vctsize*sizeof(double));
if ( ndims2 == 2 )
{
for ( size_t i = 0; i < dimlen; ++i )
{
vct[i] = abuf[i*2];
vct[i+dimlen+1] = bbuf[i*2];
}
vct[dimlen] = abuf[dimlen*2-1];
vct[dimlen*2+1] = bbuf[dimlen*2-1];
}
else
{
vct[i] = abuf[i*2];
vct[i+dimlen+1] = bbuf[i*2];
for ( size_t i = 0; i < dimlen2; ++i )
{
vct[i] = abuf[i];
vct[i+dimlen+1] = bbuf[i];
}
}
vct[dimlen] = abuf[dimlen*2-1];
vct[dimlen*2+1] = bbuf[dimlen*2-1];
if ( ret == 2 && IS_NOT_EQUAL(px, 1) )
if ( p0varid1 != -1 && IS_NOT_EQUAL(px, 1) )
for ( size_t i = 0; i < dimlen+1; ++i ) vct[i] *= px;
ncvar->vct = vct;
......
src/stream_cdf_o.c
View file @ dd50a8aa
......@@ -1597,7 +1597,15 @@ void cdf_def_zaxis_hybrid_cf(stream_t *streamptr, int type, int *ncvaridp, int z
int fileID = streamptr->fileID;
if ( streamptr->ncmode == 2 ) cdf_redef(fileID);
strcpy(axisname, "lev");
char zname[CDI_MAX_NAME]; zname[0] = 0;
char zlongname[CDI_MAX_NAME]; zlongname[0] = 0;
char zunits[CDI_MAX_NAME]; zunits[0] = 0;
cdiZaxisInqKeyStr(zaxisID, CDI_KEY_NAME, CDI_MAX_NAME, zname);
//cdiZaxisInqKeyStr(zaxisID, CDI_KEY_LONGNAME, CDI_MAX_NAME, zlongname);
cdiZaxisInqKeyStr(zaxisID, CDI_KEY_UNITS, CDI_MAX_NAME, zunits);
if ( zname[0] ) strcpy(axisname, zname);
if ( zlongname[0] == 0 ) strcpy(zlongname, "hybrid sigma pressure coordinate");
if ( zunits[0] == 0 ) strcpy(zunits, "1");
int ncvarid;
cdf_def_dim(fileID, axisname, dimlen, dimID);
......@@ -1607,24 +1615,21 @@ void cdf_def_zaxis_hybrid_cf(stream_t *streamptr, int type, int *ncvaridp, int z
{
static const char sname[] = "standard_name",
sname_v[] = "atmosphere_hybrid_sigma_pressure_coordinate",
lname[] = "long_name",
lname_v[] = "hybrid sigma pressure coordinate",
units[] = "units",
units_v[] = "1",
axis[] = "axis",
axis_v[] = "Z",
direction[] = "positive",
direction_v[] = "down";
struct attTxtTab2 tab[] = {
{ sname, sname_v, sizeof (sname_v) - 1 },
{ lname, lname_v, sizeof (lname_v) - 1 },
{ units, units_v, sizeof (units_v) - 1 },
{ axis, axis_v, sizeof (axis_v) - 1 },
{ direction, direction_v, sizeof (direction_v) - 1 },
};
enum { nAtt = sizeof (tab) / sizeof (tab[0]) };
for ( size_t i = 0; i < nAtt; ++i )
cdf_put_att_text(fileID, ncvarid, tab[i].attName, tab[i].valLen, tab[i].attVal);
cdf_put_att_text(fileID, ncvarid, "long_name", strlen(zlongname), zlongname);
cdf_put_att_text(fileID, ncvarid, "units", strlen(zunits), zunits);
}
{
char txt[CDI_MAX_NAME];
......@@ -1671,16 +1676,14 @@ void cdf_def_zaxis_hybrid_cf(stream_t *streamptr, int type, int *ncvaridp, int z
cdf_put_att_text(fileID, ncvarid, "bounds", axisnameLen, axisname);
{
static const char sname[] = "standard_name",
sname_v[] = "atmosphere_hybrid_sigma_pressure_coordinate",
units[] = "units",
units_v[] = "1";
sname_v[] = "atmosphere_hybrid_sigma_pressure_coordinate";
struct attTxtTab2 tab[] = {
{ sname, sname_v, sizeof (sname_v) - 1 },
{ units, units_v, sizeof (units_v) - 1 },
};
enum { nAtt = sizeof (tab) / sizeof (tab[0]) };
for ( size_t i = 0; i < nAtt; ++i )
cdf_put_att_text(fileID, ncbvarid, tab[i].attName, tab[i].valLen, tab[i].attVal);
cdf_put_att_text(fileID, ncbvarid, "units", strlen(zunits), zunits);
}
{
char txt[CDI_MAX_NAME];
......
src/zaxis.c
View file @ dd50a8aa
......@@ -1204,13 +1204,12 @@ void zaxisResize(int zaxisID, int size)
int zaxisDuplicate(int zaxisID)
{
int zaxisIDnew;
zaxis_t *zaxisptr = zaxis_to_pointer(zaxisID);
int zaxistype = zaxisInqType(zaxisID);
int zaxissize = zaxisInqSize(zaxisID);
zaxisIDnew = zaxisCreate(zaxistype, zaxissize);
int zaxisIDnew = zaxisCreate(zaxistype, zaxissize);
zaxis_t *zaxisptrnew = zaxis_to_pointer(zaxisIDnew);
zaxis_copy(zaxisptrnew, zaxisptr);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment