The loading of the input grid was optimized, but there are still some things that can be done.

The memory layout should change to map the grid contents directly.
The pre-calculation of some constant arrays could be done once for each grid.

cpu/trunk/ewGrid.cpp

@@ -74,44 +74,56 @@ int ewLoadBathymetry()
   Dx = Re * g2r( DLon );     // in m along the equator
   Dy = Re * g2r( DLat );
 
-  for( j=1; j<=NLat; j++ ) {
-    for( i=1; i<=NLon; i++ ) {
+  /* 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 )
-        ierr = fread( &fval, sizeof(float), 1, fp );
+        fval = buf[ (j-1) * NLon + (i-1) ];
+    	//ierr = fread( &fval, sizeof(float), 1, fp );
       else
         ierr = fscanf( fp, " %f ", &fval );
 
       Node(m, iTopo) = fval;
-
-      for( k=0; k<MAX_VARS_PER_NODE-1; k++ )
-        Node(m, k) = 0;
       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<MAX_VARS_PER_NODE-2; k++ ) {
+	  Node.initMemory( k, 0 );
+  }
+
   fclose( fp );
 
+  if( !Par.dt ) { // time step not explicitly defined
 
-  // Make bathymetry from topography. Compute stable time step.
-  double dtLoc=RealMax;
-  /* FIXME: change loops */
-  for( i=1; i<=NLon; i++ ) {
-    for( j=1; j<=NLat; j++ ) {
+	// Make bathymetry from topography. Compute stable time step.
+	double dtLoc=RealMax;
 
-      m = idx(j,i);
-      Node(m, iD) = -Node(m, iTopo);
-      if( Node(m, iD) < 0 ) {
-        Node(m, iD) = 0.0;
-        continue;
-      }
+	for( i=1; i<=NLon; i++ ) {
+	  for( j=1; j<=NLat; j++ ) {
+		  m = idx(j,i);
+		  if( Node(m, iD) == 0.0f ) continue;
+		  dtLoc = My_min( dtLoc, 0.8 * (Dx*cosdeg(getLat(j))) / sqrt(Gravity*Node(m, iD)) );
+	  }
+	}
 
-      if( Node(m, iD) < Par.dmin ) Node(m, iD) = Par.dmin;
-      dtLoc = My_min( dtLoc, 0.8 * (Dx*cosdeg(getLat(j))) / sqrt(Gravity*Node(m, iD)) );
-    }
-  }
-  if( !Par.dt ) { // time step not explicitly defined
     Log.print("Stable CFL time step: %g sec", dtLoc);
     if( dtLoc > 15 ) Par.dt = 15;
     else if( dtLoc > 10 ) Par.dt = 10;
@@ -121,7 +133,6 @@ int ewLoadBathymetry()
     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.;
@@ -138,7 +149,6 @@ int ewLoadBathymetry()
     R6[j] = cosdeg( LatMin + (j-0.5)*DLat );
   }
 
-  /* FIXME: change loops */
   for( i=1; i<=NLon; i++ ) {
     for( j=1; j<=NLat; j++ ) {
 

cpu/trunk/ewNode.h

@@ -2,6 +2,8 @@
 #define EW_NODE_H
 
 #include <stdlib.h>
+#include <string.h>
+#include "easywave.h"
 
 #define CHKRET( x ) if( (x) == NULL ) return 1;
 
@@ -19,6 +21,8 @@ public:
 	virtual int copyPOIs() = 0;
 	virtual int freeMem() = 0;
 	virtual int run() = 0;
+
+	virtual void initMemory( int index, int val ) = 0;
 };
 
 class CStructNode : public CNode {
@@ -32,6 +36,17 @@ public:
 		return node[idx1][idx2];
 	}
 
+	void initMemory( int index, int val ) {
+
+		int m;
+		for( int i=1; i<=NLon; i++ ) {
+		  for( int j=1; j<=NLat; j++ ) {
+				m = idx(j,i);
+				this->operator ()(m, index) = val;
+		  }
+		}
+	}
+
 	int mallocMem() {
 
 		CHKRET( this->node = (Float*) malloc( sizeof(Float) * NLon * NLat) );
@@ -96,6 +111,14 @@ public:
 		return ((float**)&d)[idx2][idx1];
 	}
 
+	void *getBuf( int idx ) { return ((float**)&d)[idx]; }
+
+	virtual void initMemory( int index, int val ) {
+
+		memset( getBuf(index), 0, NLat * NLon * sizeof(float) );
+
+	}
+
 	virtual int mallocMem() {
 
 		CHKRET( this->d = (float*) malloc( sizeof(float) * NLon * NLat ) );