#include #include #include #include #include "easywave.h" int NLon,NLat; double LonMin,LonMax,LatMin,LatMax; double DLon,DLat; // steps in grad double Dx,Dy; // steps in m, dx must be multiplied by cos(y) before use float *R6; float *C1; float *C2; float *C3; float *C4; int ewLoadBathymetry() { FILE *fp; char fileLabel[5]; unsigned short shval; int ierr,isBin,i,j,m,k; float fval; double dval; CNode& Node = *gNode; Log.print( "Loading bathymetry from %s", Par.fileBathymetry ); // check if bathymetry file is in ascii or binary format if( (fp=fopen(Par.fileBathymetry,"rb")) == NULL ) return Err.post( Err.msgOpenFile(Par.fileBathymetry) ); memset( fileLabel, 0, 5 ); ierr = fread( fileLabel, 4, 1, fp ); if( !strcmp( fileLabel,"DSAA" ) ) isBin = 0; else if( !strcmp( fileLabel,"DSBB" ) ) isBin = 1; else return Err.post( "%s: not GRD-file!", Par.fileBathymetry ); fclose(fp); if( isBin ) { fp = fopen( Par.fileBathymetry, "rb" ); ierr = fread( fileLabel, 4, 1, fp ); ierr = fread( &shval, sizeof(unsigned short), 1, fp ); NLon = shval; ierr = fread( &shval, sizeof(unsigned short), 1, fp ); NLat = shval; } else { fp = fopen( Par.fileBathymetry, "rt" ); ierr = fscanf( fp, "%s", fileLabel ); ierr = fscanf( fp, " %d %d ", &NLon, &NLat ); } // try to allocate memory for GRIDNODE structure and for caching arrays printf_v("Size: %d %d %luMB\n", NLon, NLat, sizeof(float)*MAX_VARS_PER_NODE*NLon*NLat/1024/1024); if( Node.mallocMem() ) return Err.post( Err.msgAllocateMem() ); if( isBin ) { ierr = fread( &LonMin, sizeof(double), 1, fp ); ierr = fread( &LonMax, sizeof(double), 1, fp ); ierr = fread( &LatMin, sizeof(double), 1, fp ); ierr = fread( &LatMax, sizeof(double), 1, fp ); ierr = fread( &dval, sizeof(double), 1, fp ); ierr = fread( &dval, sizeof(double), 1, fp ); // zmin zmax } else { ierr = fscanf( fp, " %lf %lf ", &LonMin, &LonMax ); ierr = fscanf( fp, " %lf %lf ", &LatMin, &LatMax ); ierr = fscanf( fp, " %*s %*s " ); // zmin, zmax } DLon = (LonMax - LonMin)/(NLon - 1); // in degrees DLat = (LatMax - LatMin)/(NLat - 1); Dx = Re * g2r( DLon ); // in m along the equator Dy = Re * g2r( DLat ); /* NOTE: optimal would be reading everything in one step, but that does not work because rows and columns are transposed * (only possible with binary data at all) - use temporary buffer for now (consumes additional memory!) */ float *buf = new float[ NLat*NLon ]; ierr = fread( buf, sizeof(float), NLat*NLon, fp ); for( i=1; i<=NLon; i++ ) { for( j=1; j<=NLat; j++ ) { m = idx(j,i); if( isBin ) fval = buf[ (j-1) * NLon + (i-1) ]; //ierr = fread( &fval, sizeof(float), 1, fp ); else ierr = fscanf( fp, " %f ", &fval ); Node(m, iTopo) = fval; Node(m, iTime) = -1; Node(m, iD) = -fval; if( Node(m, iD) < 0 ) { Node(m, iD) = 0.0f; } else if( Node(m, iD) < Par.dmin ) { Node(m, iD) = Par.dmin; } } } delete[] buf; for( k=1; k 15 ) Par.dt = 15; else if( dtLoc > 10 ) Par.dt = 10; else if( dtLoc > 5 ) Par.dt = 5; else if( dtLoc > 2 ) Par.dt = 2; else if( dtLoc > 1 ) Par.dt = 1; else return Err.post("Bathymetry requires too small time step (<1sec)"); } // Correct bathymetry for edge artefacts for( i=1; i<=NLon; i++ ) { if( Node(idx(1,i), iD) != 0 && Node(idx(2,i), iD) == 0 ) Node(idx(1,i), iD) = 0.; if( Node(idx(NLat,i), iD) != 0 && Node(idx(NLat-1,i), iD) == 0 ) Node(idx(NLat,i), iD) = 0.; } for( j=1; j<=NLat; j++ ) { if( Node(idx(j,1), iD) != 0 && Node(idx(j,2), iD) == 0 ) Node(idx(j,1), iD) = 0.; if( Node(idx(j,NLon), iD) != 0 && Node(idx(j,NLon-1), iD) == 0 ) Node(idx(j,NLon), iD) = 0.; } // Calculate caching grid parameters for speedup for( j=1; j<=NLat; j++ ) { R6[j] = cosdeg( LatMin + (j-0.5)*DLat ); } for( i=1; i<=NLon; i++ ) { for( j=1; j<=NLat; j++ ) { m = idx(j,i); if( Node(m, iD) == 0 ) continue; Node(m, iR1) = Par.dt/Dy/R6[j]; if( i != NLon ) { if( Node(m+NLat, iD) != 0 ) { Node(m, iR2) = 0.5*Gravity*Par.dt/Dy/R6[j]*(Node(m, iD)+Node(m+NLat, iD)); Node(m, iR3) = 0.5*Par.dt*Omega*sindeg( LatMin + (j-0.5)*DLat ); } } else { Node(m, iR2) = 0.5*Gravity*Par.dt/Dy/R6[j]*Node(m, iD)*2; Node(m, iR3) = 0.5*Par.dt*Omega*sindeg( LatMin + (j-0.5)*DLat ); } if( j != NLat ) { if( Node(m+1, iD) != 0 ) { Node(m, iR4) = 0.5*Gravity*Par.dt/Dy*(Node(m, iD)+Node(m+1, iD)); Node(m, iR5) = 0.5*Par.dt*Omega*sindeg( LatMin + j*DLat ); } } /* FIXME: Bug? */ else { Node(m, iR2) = 0.5*Gravity*Par.dt/Dy*Node(m, iD)*2; Node(m, iR3) = 0.5*Par.dt*Omega*sindeg( LatMin + j*DLat ); } } } for( i=1; i<=NLon; i++ ) { C1[i] = 0; if( Node(idx(1,i), iD) != 0 ) C1[i] = 1./sqrt(Gravity*Node(idx(1,i), iD)); C3[i] = 0; if( Node(idx(NLat,i), iD) != 0 ) C3[i] = 1./sqrt(Gravity*Node(idx(NLat,i), iD)); } for( j=1; j<=NLat; j++ ) { C2[j] = 0; if( Node(idx(j,1), iD) != 0 ) C2[j] = 1./sqrt(Gravity*Node(idx(j,1), iD)); C4[j] = 0; if( Node(idx(j,NLon), iD) != 0 ) C4[j] = 1./sqrt(Gravity*Node(idx(j,NLon), iD)); } return 0; }