Commit 6e5254ca authored by Johannes Spazier's avatar Johannes Spazier

Added support for processing grid ranges that wrap around +/-180 degree longitude.

parent 3691f96e
......@@ -91,11 +91,12 @@ extern int Jmin;
extern int Jmax;
#define idx(j,i) ((i-1)*NLat+j-1)
#define getLon(i) (LonMin+(i-1)*DLon)
extern float getLon( int i );
#define getLat(j) (LatMin+(j-1)*DLat)
#define getI(lon) round( ( (lon-LonMin) / DLon) )
#define getJ(lat) round( ( (lat-LatMin) / DLat) )
extern int getI( float lon );
#define getJ(lat) (int)( (lat - LatMin) / DLat ) + 1;
int ewLoadBathymetry();
int ewParam( int argc, char **argv );
......
......@@ -38,6 +38,23 @@ float *C2;
float *C3;
float *C4;
float getLon( int i ) {
float lon = LonMin + (i-1) * DLon;
if( lon > 180.0f )
lon -= 360.0f;
return lon;
}
int getI( float lon ) {
if( lon < LonMin )
lon += 360.0f;
return (int)( (lon - LonMin) / DLon ) + 1;
}
int ewLoadBathymetry()
{
......@@ -92,6 +109,8 @@ int ewLoadBathymetry()
ierr = fscanf( fp, " %*s %*s " ); // zmin, zmax
}
printf_v("LonMin, LonMax, LatMin, LatMax: %lf, %lf, %lf, %lf\n", LonMin, LonMax, LatMin, LatMax);
DLon = (LonMax - LonMin)/(NLon - 1); // in degrees
DLat = (LatMax - LatMin)/(NLat - 1);
......@@ -105,44 +124,61 @@ int ewLoadBathymetry()
int rgL, rgR, rgT, rgB;
/* TODO: How do we know how many steps are necessary? Assume one grid extension each 10 steps for now. */
int numIter = ceil( Par.timeMax / Par.dt / 8.0 );
int numIter = ceil( Par.timeMax / 8.0 );
rgL = My_max( 0, Imin - numIter );
rgR = My_min( Imax + numIter, NLon - 1 );
int off1, off2, len, len1, len2;
rgL = off1 = Imin - numIter;
rgR = Imax + numIter;
rgT = My_max( 0, Jmin - numIter );
rgB = My_min( Jmax + numIter, NLat - 1 );
len1 = len = rgR - rgL + 1;
len2 = off2 = 0;
double LonMinSec = getLon( rgL + 1 );
double LonMaxSec = getLon( rgR + 1 );
double LatMinSec = getLat( rgT + 1 );
double LatMaxSec = getLat( rgB + 1 );
/* Use the following lines for a static range that equals "gebco_08_atlantic_2.grd". */
/*double LonMinSec = -100.004167;
double LonMaxSec = 19.995833;
double LatMinSec = -89.995833;
double LatMaxSec = 89.987500;
if( rgL < 0 ) {
off1 = 0;
len1 = len + rgL;
rgL = getI( LonMinSec );
rgR = getI( LonMaxSec );
rgT = getJ( LatMinSec );
rgB = getJ( LatMaxSec );*/
off2 = NLon + rgL - len1;
len2 = len - len1;
LonMinSec += 360.0f;
LonMaxSec += 360.0f;
} else if( rgR > NLon - 1 ) {
off1 = 0;
len1 = rgR - NLon + 1;
off2 = rgL - len1;
len2 = len - len1;
LonMaxSec += 360.0f;
}
LonMin = LonMinSec;
LonMax = LonMaxSec;
LatMin = LatMinSec;
LatMax = LatMaxSec;
printf_v("range: %lf %lf %lf %lf\n", LonMinSec, LonMaxSec, LatMinSec, LatMaxSec );
printf_v("range: %lf %lf %lf %lf\n", LonMin, (LonMax > 180.0f ? LonMax - 360.0f : LonMax), LatMin, LatMax );
/* upon here NLot and NLat are set according to the desired section */
int NLonBase = NLon;
int NLatBase = NLat;
NLon = rgR - rgL + 1;
NLat = rgB - rgT + 1;
rgL = My_max( 0, Imin - numIter );
rgR = My_min( Imax + numIter, NLon - 1 );
printf_v("NLon, NLat: %u %u\n", NLon, NLat );
NLon = len;
NLat = rgB - rgT + 1;
if( Node.mallocMem() ) return Err.post( Err.msgAllocateMem() );
......@@ -153,12 +189,21 @@ int ewLoadBathymetry()
fseek( fp, rgT * NLonBase * sizeof(float), SEEK_CUR );
float *p = buf;
for( j = 1; j <= NLat; j++ ) {
fseek ( fp, rgL * sizeof(float), SEEK_CUR );
ierr = fread( buf + (j-1) * NLon, sizeof(float), NLon, fp );
fseek ( fp, ( (NLonBase - 1) - rgR) * sizeof(float), SEEK_CUR );
fseek( fp, off1 * sizeof(float), SEEK_CUR );
ierr = fread( p + len2, sizeof(float), len1, fp );
if( len2 > 0 ) {
fseek( fp, off2 * sizeof(float), SEEK_CUR );
ierr = fread( p, sizeof(float), len2, fp );
} else {
fseek( fp, (NLonBase - 1 - rgR) * sizeof(float), SEEK_CUR );
}
p += NLon;
}
double zmin = 0.0;
......@@ -199,10 +244,10 @@ int ewLoadBathymetry()
shval = NLat;
fwrite( &shval, sizeof(unsigned short), 1, fout );
fwrite( &LonMinSec, sizeof(double), 1, fout );
fwrite( &LonMaxSec, sizeof(double), 1, fout );
fwrite( &LatMinSec, sizeof(double), 1, fout );
fwrite( &LatMaxSec, sizeof(double), 1, fout );
fwrite( &LonMin, sizeof(double), 1, fout );
fwrite( &LonMax, sizeof(double), 1, fout );
fwrite( &LatMin, sizeof(double), 1, fout );
fwrite( &LatMax, sizeof(double), 1, fout );
fwrite( &zmin, sizeof(double), 1, fout );
fwrite( &zmax, sizeof(double), 1, fout );
......
......@@ -73,6 +73,9 @@ int ewOut2D()
nOutJ = Jmax-Jmin+1;
latOutMin = getLat(Jmin); latOutMax = getLat(Jmax);
if( lonOutMin < LonMin ) lonOutMin += 360.0f;
if( lonOutMax < LonMin ) lonOutMax += 360.0f;
// write ssh
sprintf( record, "%s.2D.%5.5d.ssh", Par.modelName, Par.time );
fp = fopen( record, "wb" );
......@@ -142,6 +145,9 @@ int ewDump2D()
nOutJ = Jmax-Jmin+1;
latOutMin = getLat(Jmin); latOutMax = getLat(Jmax);
if( lonOutMin < LonMin ) lonOutMin += 360.0f;
if( lonOutMax < LonMin ) lonOutMax += 360.0f;
// write ssh max
sprintf( record, "%s.2D.sshmax", Par.modelName );
fp = fopen( record, "wb" );
......
......@@ -116,8 +116,8 @@ int ewLoadPOIs()
}
// find the closest water grid node. Local distances could be treated as cartesian (2 min cell distortion at 60 degrees is only about 2 meters or 0.2%)
i0 = (int)((lon - LonMin)/DLon) + 1;
j0 = (int)((lat - LatMin)/DLat) + 1;
i0 = getI( lon );
j0 = getJ( lat );
if( i0<1 || i0>NLon || j0<1 || j0>NLat ) {
Log.print( "!POI out of grid: %s", record );
if( Par.poiReport ) fprintf( fpRej, "%s\n", record );
......
......@@ -223,9 +223,6 @@ int setWaveHeights() {
lon = getLon(i);
lat = getLat(j);
//lon = (LonMinSec+(i-1)*DLon);
//lat = (LatMinSec+(j-1)*DLat);
if( Node(idx(j,i), iD) != 0. )
dz = Node(idx(j,i), iH) = uZ.getVal( lon,lat );
else
......
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