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

gridToCurvilinear: lambert grid support

parent d2e0c5ba
......@@ -113,6 +113,7 @@ src/griblib.c -text
src/grid.c -text
src/grid.h -text
src/grid_gme.c -text
src/grid_lambert.c -text
src/grid_rot.c -text
src/ieg.h -text
src/ieglib.c -text
......
......@@ -155,7 +155,7 @@ static void printInfo(int gridtype, int date, int time, int code, double level,
#define CDI_BIGENDIAN 0 /* Data type BIGENDIAN */
#define CDI_LITTLEENDIAN 1 /* Data type LITTLEENDIAN */
void printFiletype(int streamID)
void printFiletype(int streamID, int vlistID)
{
int filetype;
......@@ -200,9 +200,7 @@ void printFiletype(int streamID)
if ( filetype == FILETYPE_GRB )
{
int vlistID, nvars, varID;
vlistID = streamInqVlist(streamID);
int nvars, varID;
nvars = vlistNvars(vlistID);
......@@ -220,6 +218,184 @@ void printFiletype(int streamID)
}
static void printGridInfo(int vlistID)
{
static char func[] = "printGridInfo";
int ngrids, index;
int gridID, gridtype, trunc, gridsize, xsize, ysize;
int nbyte0;
ngrids = vlistNgrids(vlistID);
for ( index = 0; index < ngrids; index++ )
{
gridID = vlistGrid(vlistID, index);
gridtype = gridInqType(gridID);
trunc = gridInqTrunc(gridID);
gridsize = gridInqSize(gridID);
xsize = gridInqXsize(gridID);
ysize = gridInqYsize(gridID);
/* nbyte0 = fprintf(stdout, " %4d : %-23s : ",*/
nbyte0 = fprintf(stdout, " %4d : %-12s > ",
gridID+1, gridNamePtr(gridtype));
if ( gridtype == GRID_LONLAT ||
gridtype == GRID_GAUSSIAN ||
gridtype == GRID_GAUSSIAN_REDUCED )
{
double lonfirst = 0.0, lonlast = 0.0;
double latfirst = 0.0, latlast = 0.0;
double loninc = 0.0, latinc = 0.0;
latfirst = gridInqYval(gridID, 0);
latlast = gridInqYval(gridID, ysize-1);
latinc = gridInqYinc(gridID);
if ( gridtype == GRID_GAUSSIAN_REDUCED )
{
fprintf(stdout, "size : dim = %d nlat = %d\n", gridsize, ysize);
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "longitude : reduced\n");
}
else
{
lonfirst = gridInqXval(gridID, 0);
lonlast = gridInqXval(gridID, xsize-1);
loninc = gridInqXinc(gridID);
fprintf(stdout, "size : dim = %d nlon = %d nlat = %d\n", gridsize, xsize, ysize);
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "longitude : first = %.9g last = %.9g", lonfirst, lonlast);
if ( !DBL_IS_EQUAL(loninc, 0) )
fprintf(stdout, " inc = %.9g", loninc);
if ( gridIsCyclic(gridID) )
fprintf(stdout, " circular");
fprintf(stdout, "\n");
}
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "latitude : first = %.9g last = %.9g", latfirst, latlast);
if ( !DBL_IS_EQUAL(latinc, 0) && gridtype == GRID_LONLAT )
fprintf(stdout, " inc = %.9g", latinc);
fprintf(stdout, "\n");
if ( gridIsRotated(gridID) )
{
double lonpole, latpole;
lonpole = gridInqXpole(gridID);
latpole = gridInqYpole(gridID);
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "northpole : lon = %.9g lat = %.9g\n", lonpole, latpole);
}
if ( gridInqXbounds(gridID, NULL) || gridInqYbounds(gridID, NULL) )
{
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "available :");
if ( gridInqXbounds(gridID, NULL) ) fprintf(stdout, " xbounds");
if ( gridInqYbounds(gridID, NULL) ) fprintf(stdout, " ybounds");
fprintf(stdout, "\n");
}
}
else if ( gridtype == GRID_SPECTRAL )
{
fprintf(stdout, "size : dim = %d truncation = %d spc = %d\n",
gridsize, trunc, gridsize/2);
}
else if ( gridtype == GRID_GME )
{
int ni, nd;
ni = gridInqGMEni(gridID);
nd = gridInqGMEnd(gridID);
fprintf(stdout, "size : dim = %d nd = %d ni = %d\n", gridsize, nd, ni);
}
else if ( gridtype == GRID_CURVILINEAR )
{
fprintf(stdout, "size : dim = %d nx = %d ny = %d\n", gridsize, xsize, ysize);
if ( gridInqXvals(gridID, NULL) && gridInqYvals(gridID, NULL) )
{
int i;
double *xvals, *yvals;
double xfirst, xlast, yfirst, ylast;
xvals = (double *) malloc(gridsize*sizeof(double));
yvals = (double *) malloc(gridsize*sizeof(double));
gridInqXvals(gridID, xvals);
gridInqYvals(gridID, yvals);
xfirst = xvals[0];
xlast = xvals[0];
yfirst = yvals[0];
ylast = yvals[0];
for ( i = 1; i < gridsize; i++ )
{
if ( xvals[i] < xfirst ) xfirst = xvals[i];
if ( xvals[i] > xlast ) xlast = xvals[i];
if ( yvals[i] < yfirst ) yfirst = yvals[i];
if ( yvals[i] > ylast ) ylast = yvals[i];
}
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "longitude : first = %.9g last = %.9g", xfirst, xlast);
if ( gridIsCyclic(gridID) )
fprintf(stdout, " circular");
fprintf(stdout, "\n");
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "latitude : first = %.9g last = %.9g\n", yfirst, ylast);
free(xvals);
free(yvals);
}
}
else if ( gridtype == GRID_CELL )
{
fprintf(stdout, "size : dim = %d nvertex = %d\n", gridsize, gridInqNvertex(gridID));
}
else if ( gridtype == GRID_LAMBERT )
{
double originLon, originLat, lonParY, lat1, lat2, xincm, yincm;
gridInqLambert(gridID, &originLon, &originLat, &lonParY, &lat1, &lat2, &xincm, &yincm);
fprintf(stdout, "size : dim = %d nx = %d ny = %d\n", gridsize, xsize, ysize);
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, " originLon = %g originLat = %g lonParY = %g\n",
originLon, originLat, lonParY);
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, " lat1 = %g lat2 = %g xinc = %gm yinc = %gm\n",
lat1, lat2, xincm, yincm);
}
else /* if ( gridtype == GRID_GENERIC ) */
{
if ( ysize == 0 )
fprintf(stdout, "size : dim = %d\n", gridsize);
else
{
fprintf(stdout, "size : dim = %d nx = %d ny = %d\n", gridsize, xsize, ysize);
if ( gridIsCyclic(gridID) )
{
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "longitude : circular\n");
}
}
}
if ( gridtype == GRID_CURVILINEAR || gridtype == GRID_CELL ||
gridtype == GRID_GENERIC || gridtype == GRID_LAMBERT )
{
if ( gridInqXvals(gridID, NULL) || gridInqYvals(gridID, NULL) || gridHasArea(gridID) ||
gridInqXbounds(gridID, NULL) || gridInqYbounds(gridID, NULL) )
{
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "available :");
if ( gridInqXvals(gridID, NULL) ) fprintf(stdout, " xvals");
if ( gridInqYvals(gridID, NULL) ) fprintf(stdout, " yvals");
if ( gridInqXbounds(gridID, NULL) ) fprintf(stdout, " xbounds");
if ( gridInqYbounds(gridID, NULL) ) fprintf(stdout, " ybounds");
if ( gridHasArea(gridID) ) fprintf(stdout, " area");
fprintf(stdout, "\n");
}
}
}
}
static void printShortinfo(int streamID, int vlistID, int vardis)
{
int varID;
......@@ -227,12 +403,11 @@ static void printShortinfo(int streamID, int vlistID, int vardis)
int gridID, zaxisID, code;
int zaxistype, ltype;
int vdate, vtime;
int nrecs, nvars, ngrids, nzaxis, ntsteps;
int nrecs, nvars, nzaxis, ntsteps;
int levelID, levelsize;
int tsID, ntimeout;
int xsize, ysize, trunc;
int timeID, taxisID;
int gridtype, nbyte, nbyte0;
int nbyte, nbyte0;
int index;
char varname[128];
char longname[128];
......@@ -244,7 +419,7 @@ static void printShortinfo(int streamID, int vlistID, int vardis)
char pstr[4];
printf(" File format: ");
printFiletype(streamID);
printFiletype(streamID, vlistID);
if ( vardis )
fprintf(stdout,
......@@ -321,145 +496,8 @@ static void printShortinfo(int streamID, int vlistID, int vardis)
fprintf(stdout, "\n");
}
ngrids = vlistNgrids(vlistID);
fprintf(stdout, " Horizontal grids :\n");
for ( index = 0; index < ngrids; index++ )
{
gridID = vlistGrid(vlistID, index);
gridtype = gridInqType(gridID);
trunc = gridInqTrunc(gridID);
gridsize = gridInqSize(gridID);
xsize = gridInqXsize(gridID);
ysize = gridInqYsize(gridID);
/* nbyte0 = fprintf(stdout, " %4d : %-23s : ",*/
nbyte0 = fprintf(stdout, " %4d : %-12s > ",
gridID+1, gridNamePtr(gridtype));
if ( gridtype == GRID_LONLAT ||
gridtype == GRID_GAUSSIAN ||
gridtype == GRID_GAUSSIAN_REDUCED )
{
double lonfirst = 0.0, lonlast = 0.0;
double latfirst = 0.0, latlast = 0.0;
double loninc = 0.0, latinc = 0.0;
latfirst = gridInqYval(gridID, 0);
latlast = gridInqYval(gridID, ysize-1);
latinc = gridInqYinc(gridID);
if ( gridtype == GRID_GAUSSIAN_REDUCED )
{
fprintf(stdout, "size : dim = %d nlat = %d\n", gridsize, ysize);
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "longitude : reduced\n");
}
else
{
lonfirst = gridInqXval(gridID, 0);
lonlast = gridInqXval(gridID, xsize-1);
loninc = gridInqXinc(gridID);
fprintf(stdout, "size : dim = %d nlon = %d nlat = %d\n", gridsize, xsize, ysize);
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "longitude : first = %.9g last = %.9g", lonfirst, lonlast);
if ( !DBL_IS_EQUAL(loninc, 0) )
fprintf(stdout, " inc = %.9g", loninc);
if ( gridIsCyclic(gridID) )
fprintf(stdout, " cyclic");
fprintf(stdout, "\n");
}
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "latitude : first = %.9g last = %.9g", latfirst, latlast);
if ( !DBL_IS_EQUAL(latinc, 0) && gridtype == GRID_LONLAT )
fprintf(stdout, " inc = %.9g", latinc);
fprintf(stdout, "\n");
if ( gridIsRotated(gridID) )
{
double lonpole, latpole;
lonpole = gridInqXpole(gridID);
latpole = gridInqYpole(gridID);
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "northpole : lon = %.9g lat = %.9g\n", lonpole, latpole);
}
if ( gridInqXbounds(gridID, NULL) || gridInqYbounds(gridID, NULL) )
{
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "available :");
if ( gridInqXbounds(gridID, NULL) ) fprintf(stdout, " xbounds");
if ( gridInqYbounds(gridID, NULL) ) fprintf(stdout, " ybounds");
fprintf(stdout, "\n");
}
}
else if ( gridtype == GRID_SPECTRAL )
{
fprintf(stdout, "size : dim = %d truncation = %d spc = %d\n",
gridsize, trunc, gridsize/2);
}
else if ( gridtype == GRID_GME )
{
int ni, nd;
ni = gridInqGMEni(gridID);
nd = gridInqGMEnd(gridID);
fprintf(stdout, "size : dim = %d nd = %d ni = %d\n", gridsize, nd, ni);
}
else if ( gridtype == GRID_CURVILINEAR )
{
fprintf(stdout, "size : dim = %d nx = %d ny = %d\n", gridsize, xsize, ysize);
if ( gridIsCyclic(gridID) )
{
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "longitude : cyclic\n");
}
}
else if ( gridtype == GRID_CELL )
{
fprintf(stdout, "size : dim = %d nvertex = %d\n", gridsize, gridInqNvertex(gridID));
}
else if ( gridtype == GRID_LAMBERT )
{
double originLon, originLat, lonParY, lat1, lat2, xincm, yincm;
gridInqLambert(gridID, &originLon, &originLat, &lonParY, &lat1, &lat2, &xincm, &yincm);
fprintf(stdout, "size : dim = %d nx = %d ny = %d\n", gridsize, xsize, ysize);
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, " originLon = %g originLat = %g lonParY = %g\n",
originLon, originLat, lonParY);
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, " lat1 = %g lat2 = %g xinc = %gm yinc = %gm\n",
lat1, lat2, xincm, yincm);
}
else /* if ( gridtype == GRID_GENERIC ) */
{
if ( ysize == 0 )
fprintf(stdout, "size : dim = %d\n", gridsize);
else
{
fprintf(stdout, "size : dim = %d nx = %d ny = %d\n", gridsize, xsize, ysize);
if ( gridIsCyclic(gridID) )
{
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "longitude : cyclic\n");
}
}
}
if ( gridtype == GRID_CURVILINEAR || gridtype == GRID_CELL ||
gridtype == GRID_GENERIC || gridtype == GRID_LAMBERT )
{
if ( gridInqXvals(gridID, NULL) || gridInqYvals(gridID, NULL) || gridHasArea(gridID) ||
gridInqXbounds(gridID, NULL) || gridInqYbounds(gridID, NULL) )
{
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "available :");
if ( gridInqXvals(gridID, NULL) ) fprintf(stdout, " xvals");
if ( gridInqYvals(gridID, NULL) ) fprintf(stdout, " yvals");
if ( gridInqXbounds(gridID, NULL) ) fprintf(stdout, " xbounds");
if ( gridInqYbounds(gridID, NULL) ) fprintf(stdout, " ybounds");
if ( gridHasArea(gridID) ) fprintf(stdout, " area");
fprintf(stdout, "\n");
}
}
}
printGridInfo(vlistID);
nzaxis = vlistNzaxis(vlistID);
fprintf(stdout, " Vertical grids :\n");
......
......@@ -31,6 +31,7 @@ libcdi_a_SOURCES = \
grid.c \
grid_rot.c \
grid_gme.c \
grid_lambert.c \
zaxis.c \
stream.c \
stream_var.c \
......
......@@ -113,6 +113,7 @@ libcdi_a_SOURCES = \
grid.c \
grid_rot.c \
grid_gme.c \
grid_lambert.c \
zaxis.c \
stream.c \
stream_var.c \
......@@ -180,11 +181,12 @@ am_libcdi_a_OBJECTS = cdiFortran.$(OBJEXT) cdi_error.$(OBJEXT) \
stream_history.$(OBJEXT) stream_grb.$(OBJEXT) \
stream_cdf.$(OBJEXT) stream_srv.$(OBJEXT) stream_ext.$(OBJEXT) \
stream_ieg.$(OBJEXT) grid.$(OBJEXT) grid_rot.$(OBJEXT) \
grid_gme.$(OBJEXT) zaxis.$(OBJEXT) stream.$(OBJEXT) \
stream_var.$(OBJEXT) stream_record.$(OBJEXT) tsteps.$(OBJEXT) \
stream_int.$(OBJEXT) servicelib.$(OBJEXT) extralib.$(OBJEXT) \
ieglib.$(OBJEXT) cdf.$(OBJEXT) cdf_int.$(OBJEXT) file.$(OBJEXT) \
binary.$(OBJEXT) swap.$(OBJEXT) griblib.$(OBJEXT)
grid_gme.$(OBJEXT) grid_lambert.$(OBJEXT) zaxis.$(OBJEXT) \
stream.$(OBJEXT) stream_var.$(OBJEXT) stream_record.$(OBJEXT) \
tsteps.$(OBJEXT) stream_int.$(OBJEXT) servicelib.$(OBJEXT) \
extralib.$(OBJEXT) ieglib.$(OBJEXT) cdf.$(OBJEXT) \
cdf_int.$(OBJEXT) file.$(OBJEXT) binary.$(OBJEXT) \
swap.$(OBJEXT) griblib.$(OBJEXT)
libcdi_a_OBJECTS = $(am_libcdi_a_OBJECTS)
DEFS = @DEFS@
......@@ -200,11 +202,11 @@ am__depfiles_maybe = depfiles
@AMDEP_TRUE@ ./$(DEPDIR)/error.Po ./$(DEPDIR)/extralib.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/file.Po ./$(DEPDIR)/griblib.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/grid.Po ./$(DEPDIR)/grid_gme.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/grid_rot.Po ./$(DEPDIR)/ieglib.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/institution.Po ./$(DEPDIR)/model.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/servicelib.Po ./$(DEPDIR)/stream.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/stream_cdf.Po ./$(DEPDIR)/stream_ext.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/stream_grb.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/grid_lambert.Po ./$(DEPDIR)/grid_rot.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/ieglib.Po ./$(DEPDIR)/institution.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/model.Po ./$(DEPDIR)/servicelib.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/stream.Po ./$(DEPDIR)/stream_cdf.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/stream_ext.Po ./$(DEPDIR)/stream_grb.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/stream_history.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/stream_ieg.Po ./$(DEPDIR)/stream_int.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/stream_record.Po \
......@@ -308,6 +310,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/griblib.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/grid.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/grid_gme.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/grid_lambert.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/grid_rot.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ieglib.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/institution.Po@am__quote@
......
......@@ -3373,9 +3373,10 @@ int gridToCurvilinear(int gridID1)
{
case GRID_LONLAT:
case GRID_GAUSSIAN:
case GRID_LAMBERT:
{
int i, j, nx, ny;
double *xvals, *yvals;
double *xvals = NULL, *yvals = NULL;
double *xvals2D, *yvals2D;
double *xbounds = NULL, *ybounds = NULL;
double *xbounds2D, *ybounds2D;
......@@ -3386,40 +3387,73 @@ int gridToCurvilinear(int gridID1)
gridDefXsize(gridID2, nx);
gridDefYsize(gridID2, ny);
if ( ! (gridInqXvals(gridID1, NULL) && gridInqYvals(gridID1, NULL)) )
Error(func, "Grid has no values");
xvals = (double *) malloc(nx*sizeof(double));
yvals = (double *) malloc(ny*sizeof(double));
xvals2D = (double *) malloc(gridsize*sizeof(double));
yvals2D = (double *) malloc(gridsize*sizeof(double));
gridInqXvals(gridID1, xvals);
gridInqYvals(gridID1, yvals);
if ( gridIsRotated(gridID1) )
if ( gridtype == GRID_LAMBERT )
{
double xpole, ypole;
xpole = gridInqXpole(gridID1);
ypole = gridInqYpole(gridID1);
double xi, xj;
double originLon, originLat, lonParY, lat1, lat2, xincm, yincm;
double zlat, zlon;
int status;
gridInqLambert(gridID1, &originLon, &originLat, &lonParY, &lat1, &lat2, &xincm, &yincm);
/*
while ( originLon < 0 ) originLon += 360;
while ( lonParY < 0 ) lonParY += 360;
*/
if ( ! DBL_IS_EQUAL(xincm, yincm) )
Warning(func, "X and Y increment must be equal on Lambert grid (Xinc = %g, Yinc = %g)\n",
xincm, yincm);
if ( ! DBL_IS_EQUAL(lat1, lat2) )
Warning(func, "Lat1 and Lat2 must be equal on Lambert grid (Lat1 = %g, Lat2 = %g)\n",
lat1, lat2);
for ( j = 0; j < ny; j++ )
for ( i = 0; i < nx; i++ )
{
xvals2D[j*nx+i] = rls_to_rl(yvals[j], xvals[i], ypole, xpole);
yvals2D[j*nx+i] = phs_to_ph(yvals[j], xvals[i], ypole);
}
xi = i+1;
xj = j+1;
status = W3FB12(xi, xj, originLat, originLon, xincm, lonParY, lat1, &zlat,&zlon);
xvals2D[j*nx+i] = zlon;
yvals2D[j*nx+i] = zlat;
}
}
else
{
for ( j = 0; j < ny; j++ )
for ( i = 0; i < nx; i++ )
{
xvals2D[j*nx+i] = xvals[i];
yvals2D[j*nx+i] = yvals[j];
}
if ( ! (gridInqXvals(gridID1, NULL) && gridInqYvals(gridID1, NULL)) )
Error(func, "Grid has no values");
xvals = (double *) malloc(nx*sizeof(double));
yvals = (double *) malloc(ny*sizeof(double));
gridInqXvals(gridID1, xvals);
gridInqYvals(gridID1, yvals);
if ( gridIsRotated(gridID1) )
{
double xpole, ypole;
xpole = gridInqXpole(gridID1);
ypole = gridInqYpole(gridID1);
for ( j = 0; j < ny; j++ )
for ( i = 0; i < nx; i++ )
{
xvals2D[j*nx+i] = rls_to_rl(yvals[j], xvals[i], ypole, xpole);
yvals2D[j*nx+i] = phs_to_ph(yvals[j], xvals[i], ypole);
}
}
else
{
for ( j = 0; j < ny; j++ )
for ( i = 0; i < nx; i++ )
{
xvals2D[j*nx+i] = xvals[i];
yvals2D[j*nx+i] = yvals[j];
}
}
}
gridDefXvals(gridID2, xvals2D);
......@@ -3428,54 +3462,106 @@ int gridToCurvilinear(int gridID1)
free(xvals2D);
free(yvals2D);
if ( gridInqXboundsPtr(gridID1) )
{
xbounds = (double *) malloc(2*nx*sizeof(double));
gridInqXbounds(gridID1, xbounds);
}
else if ( nx > 1 )
{
xbounds = (double *) malloc(2*nx*sizeof(double));
gridGenXbounds(nx, xvals, xbounds);
}
if ( gridInqYboundsPtr(gridID1) )
{
ybounds = (double *) malloc(2*ny*sizeof(double));
gridInqYbounds(gridID1, ybounds);
}
else if ( ny > 1 )
{
ybounds = (double *) malloc(2*ny*sizeof(double));
gridGenYbounds(ny, yvals, ybounds);
}
if ( gridtype == GRID_LAMBERT )
{
double xi, xj;
double originLon, originLat, lonParY, lat1, lat2, xincm, yincm;
double zlat, zlon;
int status;
int index;
free(xvals);
free(yvals);
gridInqLambert(gridID1, &originLon, &originLat, &lonParY, &lat1, &lat2, &xincm, &yincm);
if ( xbounds && ybounds )
{
xbounds2D = (double *) malloc(4*gridsize*sizeof(double));
ybounds2D = (double *) malloc(4*gridsize*sizeof(double));
if ( gridIsRotated(gridID1) )
{
gridGenRotBounds(gridID1, nx, ny, xbounds, ybounds, xbounds2D, ybounds2D);
}
else
{
gridGenXbounds2D(nx, ny, xbounds, xbounds2D);
gridGenYbounds2D(nx, ny, ybounds, ybounds2D);
}
for ( j = 0; j < ny; j++ )
for ( i = 0; i < nx; i++ )
{
index = j*4*nx + 4*i;
xi = i+1.5;
xj = j+1.5;
status = W3FB12(xi, xj, originLat, originLon, xincm, lonParY, lat1, &zlat,&zlon);
xbounds2D[index+0] = zlon;
ybounds2D[index+0] = zlat;