runPIHM/pihm/initialize.cpp

Go to the documentation of this file.
00001 /*******************************************************************************
00002  * File        : initialize.c                                                  *
00003  * Function    : initialization of data structue                               *
00004  * Programmer  : Yizhong Qu @ Pennsylvania State Univeristy                    *
00005  * Version     : May, 2004 (1.0)                                               *
00006  *-----------------------------------------------------------------------------*
00007  *                                                                             *
00008  * Part of initialization work is done in read_alloc, and the rest is done here*
00009  *                                                                             *
00010  * This code is free for users with research purpose only, if appropriate      *
00011  * citation is refered. However, there is no warranty in any format for this   *
00012  * product.                                                                    *
00013  *                                                                             *
00014  * For questions or comments, please contact the authors of the reference.     *
00015  * One who want to use it for other consideration may also contact Dr.Duffy    *
00016  * at cxd11@psu.edu.                                                           *
00017  *******************************************************************************/
00018 
00019 #include <stdio.h>
00020 #include <stdlib.h>
00021 #include <math.h>
00022 #include <string.h>
00023 
00024 #include "sundials_types.h"
00025 #include "nvector_serial.h"
00026 #include "ihm.h"
00027 #include "calib.h"
00028 
00029 //#define satD_CALIB    0.0
00030 //#define br_CALIB      4.0
00031 //#define poros_CALIB   0.35
00032 //#define icsat_CALIB   0.56
00033 //#define rivEle_CALIB  0.62
00034 
00035 realtype satD_CALIB;
00036 realtype br_CALIB;
00037 realtype poros_CALIB;
00038 realtype icsat_CALIB;
00039 realtype rivEle_CALIB;
00040 
00041 int lbool;
00042 
00043 void initialize(char *filename, Model_Data DS, Control_Data *CS, N_Vector CV_Y)
00044 {
00045         int i,j,counterMin,counterMax,MINCONST,domcounter;
00046         realtype a_x, a_y, b_x, b_y, c_x, c_y,MAXCONST;
00047         realtype a_zmin, a_zmax, b_zmin, b_zmax, c_zmin, c_zmax;
00048         realtype tempvalue;
00049         FILE *int_file;
00050         char *fn;
00051         realtype *zmin_cor;
00052 
00053         satD_CALIB = setsatD_CALIB();
00054         br_CALIB   = setbr_CALIB();
00055         poros_CALIB = setporos_CALIB();
00056         icsat_CALIB = seticsat_CALIB();
00057         rivEle_CALIB = setrivEle_CALIB();
00058 
00059         zmin_cor=(realtype *)malloc(DS->NumEle*sizeof(realtype));
00060 
00061         printf("\nInitializing data structure ... ");
00062 
00063         for(i=0; i<DS->NumEle; i++)
00064         {
00065         a_x = DS->Node[DS->Ele[i].node[0]-1].x;
00066         b_x = DS->Node[DS->Ele[i].node[1]-1].x;
00067         c_x = DS->Node[DS->Ele[i].node[2]-1].x;
00068         a_y = DS->Node[DS->Ele[i].node[0]-1].y;
00069         b_y = DS->Node[DS->Ele[i].node[1]-1].y;
00070         c_y = DS->Node[DS->Ele[i].node[2]-1].y;
00071 
00072         a_zmin = DS->Node[DS->Ele[i].node[0]-1].zmin;
00073         b_zmin = DS->Node[DS->Ele[i].node[1]-1].zmin;
00074         c_zmin = DS->Node[DS->Ele[i].node[2]-1].zmin;
00075         a_zmax = DS->Node[DS->Ele[i].node[0]-1].zmax;
00076         b_zmax = DS->Node[DS->Ele[i].node[1]-1].zmax;
00077         c_zmax = DS->Node[DS->Ele[i].node[2]-1].zmax;
00078 
00079         MINCONST=10000000;
00080         MAXCONST=-1000000;
00081         for(j=0;j<3;j++)
00082         {
00083                 if(DS->Node[DS->Ele[i].node[j]-1].zmin<MINCONST)
00084                 {
00085                         MINCONST=DS->Node[DS->Ele[i].node[j]-1].zmin;
00086                         counterMin=j;
00087                 }
00088                 if(DS->Node[DS->Ele[i].node[j]-1].zmax>MAXCONST)
00089                 {
00090                         MAXCONST=DS->Node[DS->Ele[i].node[j]-1].zmax;
00091                         counterMax=j;
00092                 }
00093         }
00094         DS->Ele[i].NodeZmin= DS->Node[DS->Ele[i].node[counterMin]-1].zmin;
00095         DS->Ele[i].NodeZmax= DS->Node[DS->Ele[i].node[counterMax]-1].zmax;
00096         DS->Ele[i].NodeDist= sqrt(pow(DS->Node[DS->Ele[i].node[counterMin]-1].x-DS->Node[DS->Ele[i].node[counterMax]-1].x,2)+pow(DS->Node[DS->Ele[i].node[counterMin]-1].y-DS->Node[DS->Ele[i].node[counterMax]-1].y,2));
00097 
00098         DS->Ele[i].area = 0.5*((b_x - a_x)*(c_y - a_y) - (b_y - a_y)*(c_x - a_x));
00099                 //printf("\n%lf",DS->Ele[i].area);
00100         DS->Ele[i].zmax = (a_zmax + b_zmax + c_zmax)/3.0;
00101         DS->Ele[i].zmin = (a_zmin + b_zmin + c_zmin)/3.0;
00102         //DS->Ele[i].zmin =DS->Ele[i].zmax-br_CALIB;
00103                 /*MAXCONST=-1000000;
00104                 domcounter=0;
00105         for(j=0;j<3;j++)
00106         {
00107                         if(fabs(DS->Ele[i].zmax-DS->Ele[DS->Ele[i].nabr[j]-1].zmax)>MAXCONST)
00108                         {
00109                                 MAXCONST=fabs(DS->Ele[i].zmax-DS->Ele[DS->Ele[i].nabr[j]-1].zmax);
00110                                 domcounter=DS->Ele[i].zmax-DS->Ele[DS->Ele[i].nabr[j]-1].zmax<0?j:j+3;
00111                         }
00112         }
00113                 printf("\n%d %d",i,domcounter);
00114                 */
00115         DS->Ele[i].edge[0] = pow((b_x - c_x), 2) + pow((b_y - c_y), 2);
00116         DS->Ele[i].edge[1] = pow((c_x - a_x), 2) + pow((c_y - a_y), 2);
00117         DS->Ele[i].edge[2] = pow((a_x - b_x), 2) + pow((a_y - b_y), 2);
00118         /* calculate centroid of triangle */
00119 
00120         DS->Ele[i].x = (a_x + b_x + c_x)/3.0;
00121         DS->Ele[i].y = (a_y + b_y + c_y)/3.0;
00122 
00123 
00124         /* calculate circumcenter of triangle */
00125                 /* DS->Ele[i].x = a_x - ((b_y - a_y)*DS->Ele[i].edge[2] - (c_y - a_y)*DS->Ele[i].edge[0])/(4*DS->Ele[i].area);
00126         DS->Ele[i].y = a_y + ((b_x - a_x)*DS->Ele[i].edge[2] - (c_x - a_x)*DS->Ele[i].edge[0])/(4*DS->Ele[i].area);
00127         */
00128         DS->Ele[i].edge[0] = sqrt(DS->Ele[i].edge[0]);
00129         DS->Ele[i].edge[1] = sqrt(DS->Ele[i].edge[1]);
00130         DS->Ele[i].edge[2] = sqrt(DS->Ele[i].edge[2]);
00131         /***************** Temporary ET inputs ***********************************/
00132         DS->Ele[i].windH = 10;
00133         /*************************************************************************/
00134 
00135         }
00136         /* Routine to Find Sinks 
00137         for(i=0; i<DS->NumEle; i++){
00138                 lbool=0;
00139                 for(j=0; j<3; j++){
00140                         if(DS->Ele[i].nabr[j]>0){
00141                                 if(DS->Ele[i].zmax+1.0>DS->Ele[DS->Ele[i].nabr[j]-1].zmax)
00142                                         lbool=1;
00143                         }
00144                 }
00145                 if(lbool==0){
00146                         printf("%d is sink\n",i+1);
00147                 }
00148         }
00149         */      
00150         //getchar();
00151         /* allocate flux */
00152         DS->FluxSurf = (realtype **)malloc(DS->NumEle*sizeof(realtype));
00153         DS->FluxSub = (realtype **)malloc(DS->NumEle*sizeof(realtype));
00154         DS->FluxRiv = (realtype **)malloc(DS->NumRiv*sizeof(realtype));
00155         DS->EleET = (realtype **)malloc(DS->NumEle*sizeof(realtype));
00156 
00157         for(i=0; i<DS->NumEle; i++)
00158         {
00159         DS->FluxSurf[i] = (realtype *)malloc(3*sizeof(realtype));
00160         DS->FluxSub[i] = (realtype *)malloc(3*sizeof(realtype));
00161         DS->EleET[i] = (realtype *)malloc(4*sizeof(realtype));
00162         }
00163 
00164         for(i=0; i<DS->NumRiv; i++)
00165         {
00166         DS->FluxRiv[i] = (realtype *)malloc(6*sizeof(realtype));
00167         }
00168 
00169         DS->ElePrep = (realtype *)malloc(DS->NumEle*sizeof(realtype));
00170         DS->EleVic = (realtype *)malloc(DS->NumEle*sizeof(realtype));
00171         DS->Recharge = (realtype *)malloc(DS->NumEle*sizeof(realtype));
00172         DS->EleIS = (realtype *)malloc(DS->NumEle*sizeof(realtype));
00173         DS->EleISmax = (realtype *)malloc(DS->NumEle*sizeof(realtype));
00174         DS->EleSnow = (realtype *)malloc(DS->NumEle*sizeof(realtype));
00175         DS->EleTF = (realtype *)malloc(DS->NumEle*sizeof(realtype));
00176         DS->Ele2IS = (realtype *)malloc(DS->NumEle*sizeof(realtype));
00177         DS->EleNetPrep = (realtype *)malloc(DS->NumEle*sizeof(realtype));
00178 
00179         for(i=0; i<DS->NumEle; i++)
00180         {
00181         DS->Ele[i].Ksat = DS->Soil[(DS->Ele[i].soil-1)].Ksat;
00182         DS->Ele[i].Porosity = DS->Soil[(DS->Ele[i].soil-1)].SitaS -
00183                           DS->Soil[(DS->Ele[i].soil-1)].SitaR;
00184         DS->Ele[i].Porosity = poros_CALIB*DS->Ele[i].Porosity;
00185         DS->Ele[i].Alpha = DS->Soil[(DS->Ele[i].soil-1)].Alpha;
00186         DS->Ele[i].Beta = DS->Soil[(DS->Ele[i].soil-1)].Beta;
00187         DS->Ele[i].Sf = DS->Soil[(DS->Ele[i].soil-1)].Sf;
00188         DS->Ele[i].RzD=DS->Soil[DS->Ele[i].soil-1].RzD>(DS->Ele[i].zmax-DS->Ele[i].zmin)?0.5*(DS->Ele[i].zmax-DS->Ele[i].zmin):DS->Soil[DS->Ele[i].soil-1].RzD;  /* Revisit this modification */
00189 
00190         DS->Ele[i].LAImax = DS->LandC[DS->Ele[i].LC-1].LAImax;
00191         DS->Ele[i].Rmin = DS->LandC[DS->Ele[i].LC-1].Rmin;
00192         DS->Ele[i].Rs_ref = DS->LandC[DS->Ele[i].LC-1].Rs_ref;
00193         DS->Ele[i].Albedo = DS->LandC[DS->Ele[i].LC-1].Albedo;
00194         DS->Ele[i].VegFrac = DS->LandC[DS->Ele[i].LC-1].VegFrac;
00195         DS->Ele[i].Rough = DS->LandC[DS->Ele[i].LC-1].Rough;
00196         }
00197 
00198         for(i=0; i<DS->NumRiv; i++)
00199         {
00200         DS->Riv[i].x = (DS->Node[DS->Riv[i].FromNode-1].x +
00201                     DS->Node[DS->Riv[i].ToNode-1].x)/2;
00202         DS->Riv[i].y = (DS->Node[DS->Riv[i].FromNode-1].y +
00203                     DS->Node[DS->Riv[i].ToNode-1].y)/2;
00204         DS->Riv[i].zmax = (DS->Node[DS->Riv[i].FromNode-1].zmax + DS->Node[DS->Riv[i].ToNode-1].zmax)/2;
00205         DS->Riv[i].depth = DS->Riv_Shape[DS->Riv[i].shape-1].depth;
00206         DS->Riv[i].zmin = DS->Riv[i].zmax - DS->Riv[i].depth;
00207 
00208         DS->Riv[i].Length = sqrt(pow(DS->Node[DS->Riv[i].FromNode-1].x -
00209                                  DS->Node[DS->Riv[i].ToNode-1].x, 2) +
00210                              pow(DS->Node[DS->Riv[i].FromNode-1].y -
00211                                  DS->Node[DS->Riv[i].ToNode-1].y, 2));
00212         }
00213         /*
00214         for(i=0; i<DS->NumRiv; i++)
00215         {
00216                 if(DS->Riv[i].down>0)
00217                 {
00218                         if(DS->Riv[i].zmin-DS->Riv[DS->Riv[i].down-1].zmin<=0.0)
00219                         {
00220                                 printf("\n%d\t%lf",i,DS->Riv[i].zmin-DS->Riv[DS->Riv[i].down-1].zmin);//DS->Ele[DS->Riv[i].LeftEle-1].zmax,DS->Ele[DS->Riv[i].RightEle-1].zmax);
00221                         }
00222                 }
00223         }
00224         getchar();
00225 */
00226 
00227         /* initialize state varible */
00228         /* relax cases */
00229         if (CS->int_type == 0)
00230         {
00231         for(i=0; i<DS->NumEle; i++)
00232         {
00233                 DS->EleIS[i] = 0;
00234                 NV_Ith_S(CV_Y, i) = 0;
00235                 NV_Ith_S(CV_Y, i + DS->NumEle) = (1/DS->Ele[i].Alpha)*(1-exp(-DS->Ele[i].Alpha*(0.1+0.05)));
00236                 NV_Ith_S(CV_Y, i + 2*DS->NumEle) = DS->Ele[i].zmax - DS->Ele[i].zmin - 0.1;
00237         }
00238 
00239         for(i=0; i<DS->NumRiv; i++)
00240         {
00241                 NV_Ith_S(CV_Y, i + 3*DS->NumEle) = 0;
00242         }
00243         }
00244         /* type mode */
00245         else if (CS->int_type == 1)
00246         {
00247                 if(DS->UnsatMode ==1)
00248                 {
00249                 for(i=0; i<DS->NumEle; i++)
00250                 {
00251                         DS->EleIS[i] = DS->Ele_IC[i].interception;
00252                         DS->EleSnow[i]=DS->Ele_IC[i].snow;
00253                         NV_Ith_S(CV_Y, i) = DS->Ele_IC[i].surf;
00254                                 /*NV_Ith_S(CV_Y, i + DS->NumEle) = DS->Ele_IC[DS->Ele[i].IC-1].unsat; */
00255                         NV_Ith_S(CV_Y, i + DS->NumEle) = DS->Ele_IC[i].sat;
00256 
00257                 }
00258 
00259                 for(i=0; i<DS->NumRiv; i++)
00260                 {
00261                         NV_Ith_S(CV_Y, i + 2*DS->NumEle) = DS->Riv_IC[DS->Riv[i].IC-1].value;
00262                 }
00263         }
00264                 if(DS->UnsatMode ==2)
00265                 {
00266                 for(i=0; i<DS->NumEle; i++)
00267                 {
00268                         /*DS->Ele_IC[i].sat= (DS->Ele_IC[i].sat>(DS->Ele[i].zmax - DS->Ele[i].zmin+zmin_cor[i]))?(DS->Ele[i].zmax - DS->Ele[i].zmin)-0.01: (DS->Ele_IC[i].sat-zmin_cor[i]>0?DS->Ele_IC[i].sat-zmin_cor[i]:0);*/
00269                         //DS->Ele_IC[i].sat=  (DS->Ele[i].zmax - DS->Ele[i].zmin)*icsat_CALIB;
00270                         DS->Ele_IC[i].sat= (DS->Ele[i].zmax-(br_CALIB*(1.0-icsat_CALIB)))>DS->Ele[i].zmin?(DS->Ele[i].zmax-DS->Ele[i].zmin-(br_CALIB*(1.0-icsat_CALIB))):(DS->Ele[i].zmax-DS->Ele[i].zmin)/2.0;
00271                         DS->EleIS[i] = DS->Ele_IC[i].interception;
00272                         NV_Ith_S(CV_Y, i) = DS->Ele_IC[i].surf;
00273 
00274                         DS->Ele_IC[i].sat=(DS->Ele[i].zmax - DS->Ele[i].zmin)-(DS->Ele_IC[i].sat+satD_CALIB)>0?((DS->Ele_IC[i].sat+satD_CALIB)<0?0.0:(DS->Ele_IC[i].sat+satD_CALIB)):(DS->Ele[i].zmax - DS->Ele[i].zmin)-0.01;
00275 
00276                         NV_Ith_S(CV_Y, i + DS->NumEle) = DS->Ele_IC[i].unsat+(1-exp(-1.0*DS->Ele[i].Alpha*((DS->Ele[i].zmax - DS->Ele[i].zmin)-DS->Ele_IC[i].sat)))/DS->Ele[i].Alpha; /* delete the later part: Just for Juniata */
00277                         NV_Ith_S(CV_Y, i + 2*DS->NumEle) = DS->Ele_IC[i].sat;
00278 
00279 
00280                         if ((NV_Ith_S(CV_Y, i + DS->NumEle) + NV_Ith_S(CV_Y, i + 2*DS->NumEle)) >= (DS->Ele[i].zmax - DS->Ele[i].zmin))
00281                         {
00282                                 NV_Ith_S(CV_Y, i + DS->NumEle) = ((DS->Ele[i].zmax - DS->Ele[i].zmin) - NV_Ith_S(CV_Y, i + 2*DS->NumEle))*0.9;
00283                                 // printf("\n here %d %e %e %e",i,NV_Ith_S(CV_Y, i + DS->NumEle),NV_Ith_S(CV_Y, i + 2*DS->NumEle),(DS->Ele[i].zmax - DS->Ele[i].zmin));
00284                                 if (NV_Ith_S(CV_Y, i + DS->NumEle) < 0)
00285                                 {
00286                                                 NV_Ith_S(CV_Y, i + DS->NumEle) = 0;
00287                                         }
00288                         }
00289                 }
00290 
00291                 for(i=0; i<DS->NumRiv; i++)
00292                 {
00293                         NV_Ith_S(CV_Y, i + 3*DS->NumEle) = DS->Riv_IC[DS->Riv[i].IC-1].value;
00294                         //NV_Ith_S(CV_Y, DS->Riv[i].LeftEle-1 + 2*DS->NumEle) = (DS->Ele[DS->Riv[i].LeftEle-1].zmax - DS->Ele[DS->Riv[i].LeftEle-1].zmin)*rivEle_CALIB;
00295                         //NV_Ith_S(CV_Y, DS->Riv[i].RightEle-1 + 2*DS->NumEle) = (DS->Ele[DS->Riv[i].RightEle-1].zmax - DS->Ele[DS->Riv[i].RightEle-1].zmin)*rivEle_CALIB;
00296                         //NV_Ith_S(CV_Y, DS->Riv[i].LeftEle-1 + DS->NumEle) = (1-exp(-1.0*DS->Ele[i].Alpha*((DS->Ele[i].zmax - DS->Ele[i].zmin)-DS->Ele_IC[i].sat)))/DS->Ele[i].Alpha;
00297                         //NV_Ith_S(CV_Y, DS->Riv[i].RightEle-1 + DS->NumEle) = (1-exp(-1.0*DS->Ele[i].Alpha*((DS->Ele[i].zmax - DS->Ele[i].zmin)-DS->Ele_IC[i].sat)))/DS->Ele[i].Alpha;
00298                 }
00299         }
00300         }
00301         /* hot start mode */
00302         else if(CS->int_type == 2)
00303         {
00304         fn = (char *)malloc((strlen(filename)+4)*sizeof(char));
00305         strcpy(fn, filename);
00306         int_file = fopen(strcat(fn, ".int"), "r");
00307 
00308         if(int_file == NULL)
00309         {
00310                 printf("\n  Fatal Error: %s.int is in use or does not exist!\n", filename);
00311                 exit(1);
00312         }
00313         else
00314         {
00315                 for(i=0; i<DS->NumEle; i++)
00316                 {
00317                         fscanf(int_file, "%lf", &tempvalue);
00318                         if(tempvalue <= 0) {tempvalue = 0.01;}
00319                         NV_Ith_S(CV_Y, i + DS->NumEle) = tempvalue;
00320                 }
00321 
00322                 for(i=0; i<DS->NumEle; i++)
00323                 {
00324                         fscanf(int_file, "%lf", &tempvalue);
00325                         if(tempvalue <= 0) {tempvalue = 0.01;}
00326                         if(tempvalue >= (DS->Ele[i].zmax - DS->Ele[i].zmin)) {tempvalue = (DS->Ele[i].zmax - DS->Ele[i].zmin) - 0.01;}
00327                         NV_Ith_S(CV_Y, i + 2*DS->NumEle) = tempvalue;
00328                 }
00329 
00330                 for(i=0; i<DS->NumEle; i++)
00331                 {
00332                         DS->EleIS[i] = 0;
00333                         NV_Ith_S(CV_Y, i) = 0;
00334                 }
00335 
00336                 for(i=0; i<DS->NumRiv; i++)
00337                 {
00338                         NV_Ith_S(CV_Y, i + 3*DS->NumEle) = 0;
00339                 }
00340         }
00341         fclose(int_file);
00342         }
00343         /**************************************************/
00344         /* Adding routine to initialize state from a file */
00345         /* with initCondition at any given time           */
00346         /**************************************************/
00347         else if(CS->int_type == 3)
00348         {
00349                 fn = (char *)malloc((strlen(filename)+5)*sizeof(char));
00350                 strcpy(fn, filename);
00351                 int_file = fopen(strcat(fn, ".init"), "r");
00352 
00353                 if(int_file == NULL)
00354                 {
00355                         printf("\n  Fatal Error: %s.int is in use or does not exist!\n", filename);
00356                     exit(1);
00357                 }
00358                 else
00359                 {
00360                         fscanf(int_file, "%lf", &tempvalue);
00361                         if(tempvalue != CS->StartTime){
00362                                 printf("\n  Fatal Error: Initial time in .init file does not match start time\n");
00363                                 exit(1);
00364                         }
00365 
00366                         if(DS->UnsatMode == 1)
00367                         {
00368                                 for(i=0; i<DS->NumEle; i++)
00369                             {
00370                                 fscanf(int_file, "%lf", &tempvalue);
00371                                 DS->EleIS[i] = tempvalue;
00372                                 fscanf(int_file, "%lf", &tempvalue);
00373                                 DS->EleSnow[i]=tempvalue;
00374                                 fscanf(int_file, "%lf", &tempvalue);
00375                                 NV_Ith_S(CV_Y, i) = tempvalue;
00376                                 /*NV_Ith_S(CV_Y, i + DS->NumEle) = DS->Ele_IC[DS->Ele[i].IC-1].unsat; */
00377                                 fscanf(int_file, "%lf", &tempvalue);
00378                                 NV_Ith_S(CV_Y, i + DS->NumEle) = tempvalue;
00379                             }
00380 
00381                             for(i=0; i<DS->NumRiv; i++)
00382                             {
00383                                 fscanf(int_file, "%lf", &tempvalue);
00384                                 NV_Ith_S(CV_Y, i + 2*DS->NumEle) = tempvalue;
00385                             }
00386                         }
00387                         if(DS->UnsatMode == 2)
00388                         {
00389                                 for(i=0; i<DS->NumEle; i++)
00390                                 {
00391                                         fscanf(int_file, "%lf", &tempvalue);
00392                                         DS->EleIS[i] = tempvalue;
00393                                         fscanf(int_file, "%lf", &tempvalue);
00394                                         DS->EleSnow[i]=tempvalue;
00395                                         fscanf(int_file, "%lf", &tempvalue);
00396                                         NV_Ith_S(CV_Y, i) = tempvalue;
00397                                         fscanf(int_file, "%lf", &tempvalue);
00398                                         NV_Ith_S(CV_Y, i + DS->NumEle) = tempvalue;
00399                                         fscanf(int_file, "%lf", &tempvalue);
00400                                         NV_Ith_S(CV_Y, i + 2*DS->NumEle) = tempvalue;
00401                                 }
00402 
00403                                 for(i=0; i<DS->NumRiv; i++)
00404                                 {
00405                                         fscanf(int_file, "%lf", &tempvalue);
00406                                         NV_Ith_S(CV_Y, i + 3*DS->NumEle) = tempvalue;
00407                                 }
00408                         }
00409                 }
00410     }
00411 
00412     printf("done.\n");
00413 }
00414
Generated on Sun Aug 5 17:33:59 2007 for PIHMgis by  doxygen 1.5.2