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

Added function gaussianLatitudesInDegrees().

parent 3640d519
......@@ -233,6 +233,27 @@ gen_grid_gme(GridDesciption &grid, const char *pline)
}
}
static void
gaussianLatitudesInDegrees(std::vector<double> &lats, std::vector<double> &lat_bounds, size_t nlat)
{
// lats(nlat)
// lat_bounds(nlat+1)
std::vector<double> latw(nlat), latw_cumsum(nlat);
gaussianLatitudes(&lats[0], &latw[0], nlat);
for (size_t j = 0; j < nlat; ++j) lats[j] = RAD2DEG*std::asin(lats[j]);
latw_cumsum[0] = latw[0];
for (size_t j = 1; j < nlat; ++j) latw_cumsum[j] = latw_cumsum[j-1] + latw[j];
lat_bounds[0] = 1.0;
for (size_t j = 1; j < nlat; ++j) lat_bounds[j] = 1.0 - latw_cumsum[j-1];
lat_bounds[nlat] = -1.0;
for (size_t j = 0; j < nlat+1; ++j) lat_bounds[j] = RAD2DEG*std::asin(lat_bounds[j]);
}
static void
gen_grid_gea(GridDesciption &grid, const char *pline)
{
......@@ -252,31 +273,14 @@ gen_grid_gea(GridDesciption &grid, const char *pline)
size_t nlat = 0.5 * polar_circumference / dy;
if (nlat&2) nlat = nlat + 1;
std::valarray<double> lats(nlat), latw(nlat), latw_cumsum(nlat), lat_bounds(nlat+1), dlat(nlat), cell_height(nlat);
std::vector<double> lats(nlat), lat_bounds(nlat+1);
gaussianLatitudesInDegrees(lats, lat_bounds, nlat);
gaussianLatitudes(&lats[0], &latw[0], nlat);
lats = RAD2DEG*std::asin(lats);
latw_cumsum[0] = latw[0];
for (size_t j = 1; j < nlat; ++j) latw_cumsum[j] = latw_cumsum[j-1] + latw[j];
std::vector<double> cell_height(nlat);
for (size_t j = 0; j < nlat; ++j) cell_height[j] = 0.25 * polar_circumference * (lat_bounds[j]-lat_bounds[j+1]) / 90.0;
lat_bounds[0] = 1.0;
for (size_t j = 1; j < nlat; ++j) lat_bounds[j] = 1.0 - latw_cumsum[j-1];
lat_bounds[nlat] = -1.0;
lat_bounds = RAD2DEG*std::asin(lat_bounds);
for (size_t j = 0; j < nlat; ++j) dlat[j] = lat_bounds[j]-lat_bounds[j+1];
cell_height = 0.25 * polar_circumference * dlat / 90.0;
/*
for ( size_t j = 0; j < nlat; ++j )
printf("%zu %10.3f %10.3f %10.3f %10.3f %10.3f %g\n",
j+1, lat_bounds[j], lats[j], lat_bounds[j+1], dlat[j], cell_height[j], latw[j]);
*/
size_t nlone = equator_circumference / dx;
if (nlone%2) nlone = nlone+1;
if (nlone%2) nlone++;
std::vector<int> rowlon(nlat);
size_t ncells = 0;
......@@ -287,7 +291,7 @@ gen_grid_gea(GridDesciption &grid, const char *pline)
const double dx_to_use = dx * dy / cell_height[j];
size_t nlon = std::max((int)std::lround(circumference / dx_to_use), 1);
if (nlon%2) nlon = nlon + 1;
if (nlon%2) nlon++;
rowlon[j] = nlon;
ncells += nlon;
......@@ -295,8 +299,8 @@ gen_grid_gea(GridDesciption &grid, const char *pline)
// printf("%zu %zu %zu %zu %g\n", ncells, nlone, nlat, nlone*nlat, 100.0*ncells/(nlone*nlat));
size_t ij = 0;
std::vector<double> lons(ncells);
size_t ij = 0;
for (size_t j = 0; j < nlat; ++j)
{
const size_t nlon = rowlon[j];
......@@ -310,10 +314,13 @@ gen_grid_gea(GridDesciption &grid, const char *pline)
grid.np = nlat/2;
grid.xvals.resize(ncells);
grid.yvals.resize(nlat);
grid.ybounds.resize(nlat*2);
grid.rowlon.resize(nlat);
for ( size_t i = 0; i < ncells; ++i ) grid.xvals[i] = lons[i];
for ( size_t j = 0; j < nlat; ++j ) grid.yvals[j] = lats[j];
for ( size_t j = 0; j < nlat; ++j ) grid.rowlon[j] = rowlon[j];
for (size_t i = 0; i < ncells; ++i) grid.xvals[i] = lons[i];
for (size_t j = 0; j < nlat; ++j) grid.yvals[j] = lats[j];
for (size_t j = 0; j < nlat; ++j) grid.ybounds[j*2+1] = lat_bounds[j];
for (size_t j = 0; j < nlat; ++j) grid.ybounds[j*2] = lat_bounds[j+1];
for (size_t j = 0; j < nlat; ++j) grid.rowlon[j] = rowlon[j];
}
}
}
......@@ -424,6 +431,16 @@ gridFromName(const char *gridnameptr)
grid.def_xfirst = true;
grid.def_yfirst = true;
grid.yvals.resize(grid.ysize);
grid.ybounds.resize(grid.ysize*2);
size_t nlat = grid.ysize;
std::vector<double> lats(nlat), lat_bounds(nlat+1);
gaussianLatitudesInDegrees(lats, lat_bounds, nlat);
for (size_t j = 0; j < nlat; ++j) grid.yvals[j] = lats[j];
for (size_t j = 0; j < nlat; ++j) grid.ybounds[j*2+1] = lat_bounds[j];
for (size_t j = 0; j < nlat; ++j) grid.ybounds[j*2] = lat_bounds[j+1];
}
}
}
......
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