ADCIRCtoF5.cpp

Convert ADCIRC ascii grid files and accompanying simulation output into F5. The levee information contained in the grid file is put into a specific triangular mesh for visualization purposes.

/*! Time-stamp: <@(#)ObjtoF5.c  19.10.2002 - 16:15:16  Werner Benger>
 *  $Id: ADCIRCtoF5.cpp,v 1.9 2007/06/15 22:36:50 georg Exp $
 *  $Log: ADCIRCtoF5.cpp,v $
 *  Revision 1.9  2007/06/15 22:36:50  georg
 *  added more type handling (skipping) and added pointer to doc on the net
 * 
 *  adjustment to new type in ADCIRC format 13.8.2010 marcel
 *
 *  added fort 64 support 11.9.2010 marcel
 */
 
 
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <assert.h>
#include <hdf5.h>
 
#include <F5/F5X.h>
#include <F5/F5F.h>
#include <F5/F5B.h>
#include <F5/F5R.h>
#include <F5/F5coordinates.h>
#include <F5/F5surface.h>
#include <F5/F5defs.h>
 
 
#include <vector>
#include <string>
#include <iostream>
 
using namespace std;
 
/*
  TODO: Should really implement command line argument parsing,
        instead of reading hard-coded file names in the current
        directory.
 
        Data format can be found here:
        http://www.adcirc.org/documentv46/Output_file_descriptions.html
        .14 and .63 are interesting...
        .14 is given on cmd line it contains vertices and connectivity and more
        for the types see fort.14 description.
 
        Added number of time steps to be read as second a command line argument.
 */
 
void AppendSurfaceData( hid_t fileID, const string&gridname, 
                        const vector<F5_vec3_point_t>&Coords,
                        const vector<F5_triangle_t>&Triangles,
                        int time_steps_to_load) 
{
char    buf[2048];
 
double  time = 0; 
//FILE*kat63 = fopen("fortkat22.63", "rt"), 
//    *kat74 = fopen("fortkat22.74", "rt"); 
FILE*kat63 = fopen("fort.63", "rt"),
    *kat64 = fopen("fort.64", "rt"),
    *kat74 = fopen("fort.74", "rt"); 
FILE*elev = 0;
 
#ifdef  R_OK
        if (access("fort.63.gz", R_OK) == 0)
        {
                elev = popen("zcat fort.63.gz", "r"); 
                if (elev==0)
                        perror("Error in reading fort.63"); 
                else
                        fputs("Fine, could open fort63 file.\n", stderr);
        } 
#endif
 
        if (!(kat63 || kat74 || kat64|| elev))
        {
                fputs("No additional field information available.\n", stderr);
                return; 
        }
 
 
unsigned int    Timesteps = 0, Nodes = 0;
 
 
        // Skip the first 3 lines, ctj 
        if (kat63)
        {
                fgets(buf, sizeof(buf), kat63); 
                fprintf(stderr, buf ); 
 
                fgets(buf, sizeof(buf), kat63);
                sscanf(buf, "%d %d", &Timesteps, &Nodes);
 
                fgets(buf, sizeof(buf), kat63); 
                fprintf(stderr, buf );
        }
        if (kat64)
        {       
                fgets(buf, sizeof(buf), kat64); 
                fprintf(stderr, buf ); 
 
                fgets(buf, sizeof(buf), kat64); 
                sscanf(buf, "%d %d", &Timesteps, &Nodes);
 
                fgets(buf, sizeof(buf), kat64); 
                fprintf(stderr, buf );
        } 
        if (kat74)
        {       
                fgets(buf, sizeof(buf), kat74); 
                fprintf(stderr, buf ); 
 
                fgets(buf, sizeof(buf), kat74); 
                sscanf(buf, "%d %d", &Timesteps, &Nodes);
 
                fgets(buf, sizeof(buf), kat74); 
                fprintf(stderr, buf );
        } 
        if (elev)
        {       
                fgets(buf, sizeof(buf), elev); 
                fprintf(stderr, buf ); 
 
                fgets(buf, sizeof(buf), elev); 
                sscanf(buf, "  %d %d", &Timesteps, &Nodes); 
                //      fprintf(stderr, "----> %s (%d: t=%d, n=%d)\n",
                //      buf,i, Timesteps, Nodes ); 
 
                fgets(buf, sizeof(buf), elev); 
                fprintf(stderr, buf );
        } 
 
 
        fprintf(stderr,"ADCIRC Surface Data: %d time steps, %d vertices\n", Timesteps, Nodes); 
        if (Nodes != Coords.size() )
        {
                fputs("Vertices and Coordinates mismatch.", stderr); 
                return;
        } 
 
        if(time_steps_to_load > 0 && (unsigned)time_steps_to_load < Timesteps )
                Timesteps = time_steps_to_load;
 
 
vector<float>   Elevation; 
vector<F5_vec3f_t> Wind; 
vector<F5_vec3f_t> DepthAverageVelocity;
 
        if (kat63 || elev) Elevation.resize( Coords.size() ); 
        if (kat74)         Wind.resize( Coords.size() ); 
        if (kat64)  DepthAverageVelocity.resize( Coords.size() );
 
F5Path*FirstSurface = 0;
        for(unsigned tstep=0; tstep < Timesteps; tstep++)
        {
                // Time unit is minutes since... driving MM5 data start.
                // one time step is 30 mins
        double  Time = time + 30*tstep; 
                fprintf(stderr, "Timestep %d -> t=%lg\n", tstep, Time); 
 
        F5Path *myGrid;
                if (!FirstSurface)
                {
                        FirstSurface = F5write_triangular_surface(fileID, Time, 
                                                                  gridname.c_str(),
                                                                  &Coords[0],  Coords.size(),
                                                                  &Triangles[0], Triangles.size() ); 
                        myGrid = FirstSurface;
                } 
                else
                        myGrid = F5Flink_triangular_surface(FirstSurface, Time,
                                                            gridname.c_str(),
                                                            0, Coords.size(), // share vertices
                                                            0, Triangles.size() // share triangles
                                                            );
 
                if (kat63)
                {
                        for(unsigned int i=0; i<Coords.size(); i++)
                        {
                        int     num;
                                fgets(buf, sizeof(buf), kat63); 
                                if ( sscanf(buf, "%d %g\n", &num, &Elevation[i] ) != 2 )
                                {
                                        fprintf(stderr, "Error can't find scalar value for vertex %d\n", i ); 
                                        return;
                                } 
 
                                // ??? No idea what this is good for ???
                                if ( Elevation[i] < -3333 ) Elevation[i] = 0.0; 
                        } 
                        fgets(buf, sizeof(buf), kat63); 
 
                        F5Fwrite_1D(myGrid, "height", Coords.size(),
                                    H5T_NATIVE_FLOAT, H5T_NATIVE_FLOAT,
                                    &Elevation[0], F5P_DEFAULT);
                }
 
                if (elev)
                {
                        for(unsigned i=0; i<Coords.size(); i++)
                        {unsigned
                        int     num;
                                fgets(buf, sizeof(buf), elev); 
                                if ( sscanf(buf, "%d %g\n", &num, &Elevation[i] ) != 2 )
                                {
                                        fprintf(stderr, "Error can't find scalar value for vertex %d\n", i ); 
                                        return;
                                }
 
                                // ??? No idea what this is good for ???
                                //if ( Elevation[i] < -3333 ) Elevation[i] = 0.0; 
                                if (num != i+1)
                                        fprintf(stderr, "Index mismatch: %s\n", buf);
                        } 
                        fgets(buf, sizeof(buf), elev); 
 
                        F5Fwrite_1D(myGrid, "elevation", Coords.size(),
                                    H5T_NATIVE_FLOAT, H5T_NATIVE_FLOAT,
                                    &Elevation[0], F5P_DEFAULT);
                }
 
 
                if (kat74)
                {
                        for(unsigned int i=0; i<Coords.size(); i++)
                        {
                        int     num;
                                fgets(buf, sizeof(buf), kat74); 
                                if ( sscanf(buf, "%d %g %g\n", &num, &Wind[i].x, &Wind[i].y ) != 3 )
                                {
                                        fprintf(stderr, "Error can't find vector value for vertex %d\n", i ); 
                                        return;
                                } 
                                Wind[i].z = 0.0;
                        } 
                        fgets(buf, sizeof(buf), kat74); 
 
                        F5Fwrite_1D(myGrid, "wind", Coords.size(),
                                    F5T_VEC3_FLOAT, F5T_VEC3_FLOAT,
                                    &Wind[0], F5P_DEFAULT);
                } 
 
                if (kat64)
                {
                        for(unsigned int i=0; i<Coords.size(); i++)
                        {
                        int     num;
                                fgets(buf, sizeof(buf), kat64); 
                                if ( sscanf(buf, "%d %g %g\n", &num, 
                                            &DepthAverageVelocity[i].x, 
                                            &DepthAverageVelocity[i].y ) != 3 )
                                {
                                        fprintf(stderr, "Error can't find vector value for vertex %d\n", i ); 
                                        return;
                                } 
                                DepthAverageVelocity[i].z = 0.0;
                        } 
                        fgets(buf, sizeof(buf), kat64); 
 
                        F5Fwrite_1D(myGrid, "depth average velocity", Coords.size(),
                                    F5T_VEC3_FLOAT, F5T_VEC3_FLOAT,
                                    &DepthAverageVelocity[0], F5P_DEFAULT);
                } 
 
 
                if (myGrid != FirstSurface)
                        F5Fclose( myGrid );
        } 
        F5Fclose(FirstSurface);
 
        if (kat63) fclose(kat63); 
        if (kat74) fclose(kat74); 
        if (kat64) fclose(kat64);
 
        if (elev) pclose(elev);
}
 
 
void flip(const vector<F5_vec3_point_t>&Coords, F5_triangle_t &T)
{
F5_vec3_point_t
        A = Coords[T.i],
        B = Coords[T.j],
        C = Coords[T.k],
 
        BA = { B.x - A.x, B.y - A.y, B.z - A.z },
        CA = { C.x - A.x, C.y - A.y, C.z - A.z };
 
double  zDir = BA.x*CA.y - CA.x*BA.y;
 
        if (zDir < 0)
        {
        int     tmp = T.j;
                T.j = T.k;
                T.k = tmp;
        }
}
 
double zDir(const vector<F5_vec3_point_t>&Coords, F5_triangle_t &T)
{
F5_vec3_point_t
        A = Coords[T.i],
        B = Coords[T.j],
        C = Coords[T.k],
 
        BA = { B.x - A.x, B.y - A.y, B.z - A.z },
        CA = { C.x - A.x, C.y - A.y, C.z - A.z };
 
double  zDir = BA.x*CA.y - CA.x*BA.y; 
 
        return zDir;
}
 
void flip(F5_triangle_t &T)
{
int     tmp = T.j; 
        T.j = T.k; 
        T.k = tmp;
}
 
 
/**
 Conversion application that reads an ADCIRC .grd file and
 writes it as a surface in F5.
 */
int     main(int argc, char*argv[])
{
        double  time = 0.0;
        hid_t   fileID;
        string  gridname;
        char    buf[1024]; 
        int     nCoords = 0, 
                nTriangles = 0;
        FILE*what; 
 
        int time_steps_to_load = -1;
 
        if (argc<2)
        {
                puts("Usage: ADCIRCtoF5 <inputfile.grd> <number of time steps>"); 
                puts("Time step argument may be ommitted to convert all time steps");
                return 1;
        } 
 
        if (argc > 2)
        {
                time_steps_to_load = atoi(argv[2]); 
        }
 
/*      putenv("F5_VERBOSITY=30");*/
 
        strcpy(buf, argv[1]);
        {
        char*c = strrchr(buf, '.');
                if (c)
                {
                        *c = 0;
                        strcat(c, ".f5");
                        fileID = H5Fcreate(buf,       /* The FileName */
                                           H5F_ACC_TRUNC,    /* Creation Flags */
                                           H5P_DEFAULT,      /* File creation property list identifier */
                                           H5P_DEFAULT);     /* File access property list identifier */
                }
                else
                {
                        puts("invalid filename.");
                        return 3;
                }
                *c = 0;
                gridname = buf;
        }
 
 
        what = fopen(argv[1], "rt");
        if (!what)
        {
                printf("Cannot open '%s'. That's not good.\n", argv[1] );
                return 2;
        }
        fgets(buf, sizeof(buf), what); 
 
        /* get no of points and no of triangles here*/ 
 
        if (fscanf(what, "%d %d\n",&nTriangles, &nCoords) !=2) 
        {
                fprintf(stderr,"Error cant find no of points and triangles in header\n"); 
                return 5;
        } 
        fprintf(stderr,"Read no of coords %d and triangles %d\n",nCoords, nTriangles); 
        fflush(stderr);
 
        vector<F5_vec3_point_t> Coords; 
        Coords.resize( nCoords );
 
        /* VERTEX Information */ 
        {for (int i=0; i<nCoords; i++)
        {
        int  index; 
        double  x,y,z;
                if (fscanf(what, "%d %lg %lg %lg\n",&index, &x, &y, &z) !=4) {
                        fprintf(stderr,"Error cant find 3 coords for index %d\n",i); 
                        return 5;
                } 
                assert(index == i+ 1); 
                Coords[i].x = x; 
                Coords[i].y = y; 
                Coords[i].z = z; 
                //need to scale down to see something .. it's difficult with 
                //meters/degrees ratio :(
 
                //some scaling... which was done before. removed by marcel
                //Coords[i].z /= 15000.0; 
                //Coords[i].z *= -1; 
                //fprintf(stderr,"Reading coords idx:%d (%f, %f, %f)\n",i,Coords[i].x,Coords[i].y,Coords[i].z);
        }} 
 
        /* allocate for the triangle buffer */ 
        vector<F5_triangle_t> Triangles; 
        Triangles.resize( nTriangles ); 
 
        vector<F5_triangle_t> Levees; 
        vector<vector<int> >  LeveeComplex;
 
        {for (int i=0; i<nTriangles; i++) 
        {
        int index; 
        int nodes_nr; 
        int idx[3]; 
                if (fscanf(what, "%d %d %d %d %d\n",&index, &nodes_nr, &idx[0], &idx[1], &idx[2]) !=5) 
                {
                        fprintf(stderr,"Error cant find 3 vertices for index %d\n",i); 
                        return 5;
                } 
 
                assert(index == i + 1); 
                assert(nodes_nr == 3); 
                Triangles[i].i = idx[0]-1; 
                Triangles[i].j = idx[1]-1; 
                Triangles[i].k = idx[2]-1; 
                //fprintf(stderr,"Reading triangles idx:%d (%d, %d, %d)\n",i,Triangles[i].i,Triangles[i].j, Triangles[i].k);
        }} 
 
 
        /*
          see http://www.adcirc.org/documentv46/fort_14.html
         */
 
        vector<vector<int> > WeirComplex;
 
int     DoID = -1;
char    *ADCIRC_ID = getenv("ADCIRC_ID");
        if (ADCIRC_ID)
                DoID = atoi( ADCIRC_ID );
 
        for(int N=0; N<2; N++)
        {
        int     Number;
                fgets(buf, sizeof(buf), what); 
                if (sscanf(buf, "%d", &Number) < 1)
                {
                        puts("ERROR in decoding object number"); 
                        break;
                } 
                printf("Got Object number line: %s", buf);
 
        int TotalNodeNumber; 
                fgets(buf, sizeof(buf), what); 
                if (sscanf(buf, "%d", &TotalNodeNumber) < 1)
                {
                        puts("ERROR in getting the total number of nodes"); 
                        break;
                }
                printf("Got Object total node number: %s", buf);
 
                for(int ID=0; ID<Number; ID++)
                {
                int nodes, type;
 
                        fgets(buf, sizeof(buf), what); 
                int got = sscanf(buf, "%d %d", &nodes, &type); 
                        //      printf("Got Object line: %s", buf);
 
                        if (got < 2)
                        {
                                if (got<1)
                                {
                                        puts("ERROR in decoding boundary type"); 
                                        printf("LINE: %s\n", buf); 
                                        return 1;
                                }
                                // open boundary 
                                for(int i=0; i<nodes; i++)
                                {
                                int     I;
                                        fgets(buf, sizeof(buf), what); 
                                        sscanf(buf, "%d\n", &I); 
                                } 
                                puts("open boundary done.");
                        } 
                        // NBVV(k,j) include this line only if IBTYPE(k) = 0, 1, 2, 10, 11, 12, 20, 21, 22, 30
                        else if ( type==0 || type==1 || type==2 || type==10 ||  type==11 ||type==11 ||
                                  type==12 || type==20 || type==21 ||  type==22  
                                  || type == 52 // new flux boundary??
                                )    // have a boundary
                        {
                        int     I; 
                                for(int i=0; i<nodes; i++)
                                {
                                        fgets(buf, sizeof(buf), what); 
                                        if (sscanf(buf, "%d\n", &I)<1)
                                        {
                                                puts("ERROR: cannot decode boundary"); 
                                                break;
                                        }
                                } 
                                printf("Boundary No. %d done\n", ID); fflush(stdout);
                        } 
                        //  NBVV(k,j), BARLANHT(k,j), BARLANCFSP(k,j ) - include this line only if IBTYPE(k) = 3, 13, 23 
                        else if ( type==3 || type==13 || type==23 )
                        {
                        int I; 
                        double height, CoeffA; 
                                WeirComplex.push_back( vector<int>() ); 
                                LeveeComplex.push_back( vector<int>() );
 
                                fgets(buf, sizeof(buf), what); 
                                if (sscanf(buf, "%d %lg %lg\n", &I, &height, &CoeffA )<3)
                                {
                                        puts("ERROR: cannot decode weirs"); 
                                        break;
                                }
                                I--; // Convert from FORTRAN to C indexing
 
                        double  HeightScale = 0.001;
 
 
                                height *= HeightScale;
                        int     VI; 
                                VI = Coords.size(); 
                        F5_vec3_point_t VIc = Coords[I]; 
                                VIc.z = height; 
                                Coords.push_back( VIc ); 
 
//                      int     VJ; 
//                              VJ = Coords.size(); 
//                      F5_vec3_point_t VJc = Coords[J]; 
//                              VJc.z = height; 
//                              Coords.push_back( VJc ); 
 
                                std::cerr << "Weir " << ID << "found - adding " 
                                          << 2*nodes << "triangles (now at "
                                          << Triangles.size() << ")" << std::endl; 
 
                                //fprintf(stderr, "Weir %d found - adding %d triangles (now at %d)\n", ID, 2*nodes, (int)Triangles.size() ); 
                                fflush(stderr); 
 
                                //bool  flipit = false; 
 
                                for(int i=0; i<nodes-1; i++)
                                {
                                int     Inew;
                                        fgets(buf, sizeof(buf), what); 
                                        if (sscanf(buf, "%d  %lg %lg\n", &Inew,  &height, &CoeffA)<3)
                                        {
                                                puts("ERROR: cannot decode weirs"); 
                                                break;
                                        } 
                                        Inew--;  // Convert from FORTRAN to C indexing 
/*
 
                                        if (DoID>0 && ID != DoID) continue;
 
                                        // Close weirs 
                                        {
                                        F5_triangle_t left, right; 
                                                left.i = I;  left.j = J;   left.k = Inew; 
                                                right.i=Inew; right.j = J; right.k = Jnew; 
 
                                                if (i==0)
                                                {
                                                        flipit = (zDir( Coords, left) <= 0.0)
                                                                 //||
                                                                 //(zDir( Coords, right) <= 0.0)
                                                                 ;
                                                }
 
                                        //if (i>1) continue;
                                        //              printf("LeftZDir: %lg\n", zDir( Coords, left) );
                                        //              printf("RightZDir: %lg\n", zDir( Coords, right) );
 
                                                WeirComplex.back().push_back( Triangles.size() );
                                                if (flipit) flip(left);
                                                Triangles.push_back( left ); 
 
                                                WeirComplex.back().push_back( Triangles.size() ); 
                                                if (flipit) flip(right);
                                                Triangles.push_back( right ); 
                                        }
*/ 
/*
                                        // Construct Levee 
                                        {
                                                height *= HeightScale;
 
                                        int     VInew; 
                                                VInew = Coords.size(); 
                                        F5_vec3_point_t VIc = Coords[Inew]; 
                                                VIc.z = height; 
                                                Coords.push_back( VIc ); 
 
                                        int     VJnew; 
                                                VJnew = Coords.size(); 
                                        F5_vec3_point_t VJc = Coords[Jnew]; 
                                                VJc.z = height; 
                                                Coords.push_back( VJc ); 
 
                                        F5_triangle_t Trngl;
 
                                                Trngl.i = I; Trngl.j = VI; Trngl.k = Inew;
                                                if (flipit) flip(Trngl);
                                                LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl );
 
                                                Trngl.i = Inew; Trngl.j = VI; Trngl.k = VInew;
                                                if (flipit) flip(Trngl);
                                                LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl );
 
 
                                                Trngl.i = VI; Trngl.j = VJ; Trngl.k = VInew;
                                                if (flipit) flip(Trngl);
                                                LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl );
 
                                                Trngl.i = VJ; Trngl.j = VJnew; Trngl.k = VInew;
                                                if (flipit) flip(Trngl);
                                                LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl );
 
 
                                                Trngl.i = VJnew; Trngl.j = J; Trngl.k = Jnew;
                                                if (flipit) flip(Trngl);
                                                LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl );
 
                                                Trngl.i = VJnew; Trngl.j = VJ; Trngl.k = J;
                                                if (flipit) flip(Trngl);
                                                LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl );
 
 
                                                VI = VInew; 
                                                VJ = VJnew;
                                        }
*/                              
 
                                        // Reuse indices
                                        I = Inew; 
                                } 
                        } 
                        // NBVV(k,j), IBCONN(k,j), BARINHT(k,j), BARINCFSB(k,j), BARINCFSP(k,j) - include this line only if IBTYPE(k) = 4, 24
                        else if (type==24 || type== 4) // have a weir!
                        {
                        int I, J; 
                        double height, CoeffA, CoeffB; 
                                WeirComplex.push_back( vector<int>() ); 
                                LeveeComplex.push_back( vector<int>() );
 
                                fgets(buf, sizeof(buf), what); 
                                if (sscanf(buf, "%d %d %lg %lg %lg\n", &I, &J, &height, &CoeffA, &CoeffB)<5)
                                {
                                        puts("ERROR: cannot decode weirs"); 
                                        break;
                                }
                                I--; J--; // Convert from FORTRAN to C indexing
 
                        double  HeightScale = 0.001;
 
 
                                height *= HeightScale;
                        int     VI; 
                                VI = Coords.size(); 
                        F5_vec3_point_t VIc = Coords[I]; 
                                VIc.z = height; 
                                Coords.push_back( VIc );
 
                        int     VJ;
                                VJ = Coords.size(); 
                        F5_vec3_point_t VJc = Coords[J]; 
                                VJc.z = height; 
                                Coords.push_back( VJc ); 
 
                                std::cerr << "Weir " << ID << " found - adding " 
                                          << 2*nodes << " triangles (now at "
                                          << Triangles.size() << ")" << std::endl; 
                                //fprintf(stderr, "Weir %d found - adding %d triangles (now at %d)\n", ID, 2*nodes, (int)Triangles.size() ); 
                                fflush(stderr);
 
                        bool    flipit = false; 
 
                                for(int i=0; i<nodes-1; i++)
                                {
                                int     Inew, Jnew;
                                        fgets(buf, sizeof(buf), what); 
                                        if (sscanf(buf, "%d %d %lg %lg %lg\n", &Inew, &Jnew, &height, &CoeffA, &CoeffB)<5)
                                        {
                                                puts("ERROR: cannot decode weirs"); 
                                                break;
                                        } 
                                        Inew--; Jnew--; // Convert from FORTRAN to C indexing
 
 
                                        if (DoID>0 && ID != DoID) continue;
 
                                        // Close weirs 
                                        {
                                        F5_triangle_t left, right; 
                                                left.i = I;  left.j = J;   left.k = Inew; 
                                                right.i=Inew; right.j = J; right.k = Jnew; 
 
                                                if (i==0)
                                                {
                                                        flipit = (zDir( Coords, left) <= 0.0)
                                                                 //||
                                                                 //(zDir( Coords, right) <= 0.0)
                                                                 ;
                                                }
 
                                        //if (i>1) continue;
                                        //              printf("LeftZDir: %lg\n", zDir( Coords, left) );
                                        //              printf("RightZDir: %lg\n", zDir( Coords, right) );
 
                                                WeirComplex.back().push_back( Triangles.size() );
                                                if (flipit) flip(left);
                                                Triangles.push_back( left ); 
 
                                                WeirComplex.back().push_back( Triangles.size() ); 
                                                if (flipit) flip(right);
                                                Triangles.push_back( right ); 
                                        }
 
                                        // Construct Levee 
                                        {
                                                height *= HeightScale;
 
                                        int     VInew; 
                                                VInew = Coords.size(); 
                                        F5_vec3_point_t VIc = Coords[Inew]; 
                                                VIc.z = height; 
                                                Coords.push_back( VIc ); 
 
                                        int     VJnew; 
                                                VJnew = Coords.size(); 
                                        F5_vec3_point_t VJc = Coords[Jnew]; 
                                                VJc.z = height; 
                                                Coords.push_back( VJc ); 
 
                                        F5_triangle_t Trngl;
 
                                                Trngl.i = I; Trngl.j = VI; Trngl.k = Inew;
                                                if (flipit) flip(Trngl);
                                                LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl );
 
                                                Trngl.i = Inew; Trngl.j = VI; Trngl.k = VInew;
                                                if (flipit) flip(Trngl);
                                                LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl );
 
 
                                                Trngl.i = VI; Trngl.j = VJ; Trngl.k = VInew;
                                                if (flipit) flip(Trngl);
                                                LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl );
 
                                                Trngl.i = VJ; Trngl.j = VJnew; Trngl.k = VInew;
                                                if (flipit) flip(Trngl);
                                                LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl );
 
 
                                                Trngl.i = VJnew; Trngl.j = J; Trngl.k = Jnew;
                                                if (flipit) flip(Trngl);
                                                LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl );
 
                                                Trngl.i = VJnew; Trngl.j = VJ; Trngl.k = J;
                                                if (flipit) flip(Trngl);
                                                LeveeComplex.back().push_back( Levees.size() ); Levees.push_back( Trngl );
 
 
                                                VI = VInew; 
                                                VJ = VJnew;
                                        }
                                        
 
                                        // Reuse indices
                                        I = Inew; 
                                        J = Jnew;
                                } 
                        } 
                        
                        else
                        {
                                puts("ERROR in decoding ADCIRC file"); 
                                printf("ERROR LINE: %s\n", buf); 
                                printf("Unsupported type = %d", type);
                                return 1;
                                break;
                        }
                } 
                printf("Boundary Complex %d done.\n", N);
        } 
 
        std::cerr << "Triangular surface: " << nCoords << " vertices, "
                  << Triangles.size() << " triangles - " 
                  << Triangles.size() - nTriangles << " for weirs" << std::endl; 
 
//      fprintf(stderr, "Triangular surface: %d vertices, %d triangles - %d for weirs\n", 
//              nCoords, (int)Triangles.size(), (int)Triangles.size() - nTriangles ); 
        fflush(stderr); 
 
        F5Path *MainGrid = F5write_triangular_surface(fileID, time, gridname.c_str(),
                                                      &Coords[0],  nCoords,
                                                      &Triangles[0], Triangles.size() );
 
        if (!MainGrid)
                puts("ERROR in writing grid...");
 
 
        F5Fclose(MainGrid); 
 
 
        std::cerr << "Levee Triangular surface: " << Coords.size() << " vertices, "
                  << Levees.size() << " triangles" << std::endl; 
 
//      fprintf(stderr, "Levee Triangular surface: %d vertices, %d triangles\n", 
//              (int)Coords.size(), (int)Levees.size() ); 
        fflush(stderr); 
 
        // No levee data found, don't write
        if( Levees.size()>0 )
        {
        F5Path *LeveeGrid = F5write_triangular_surface(fileID, time, (gridname + ".levee").c_str(),
                                                       &Coords[0],  Coords.size(),
                                                       &Levees[0], Levees.size() ); 
 
        if (!LeveeGrid)
                puts("ERROR in writing grid..."); 
 
 
        F5Fclose(LeveeGrid); 
        }
 
        // Due to the levees, we have additional vertices now. 
        // In order to write surface data, we need to crop them 
        // to their original count.
        Coords.resize( nCoords ); 
        AppendSurfaceData( fileID, gridname, Coords, Triangles, time_steps_to_load ); 
        
 
        H5Fclose(fileID);       
        return 0;
}