PIHM: initialize.c Source File

00001 /*******************************************************************************
00002  * File        : initialize.c                                                  *
00003  * Function    : initialization of data structue                               *
00004  * Programmers : Yizhong Qu   @ Pennsylvania State Univeristy                  *
00005  *               Mukesh Kumar @ Pennsylvania State Univeristy                  *
00006  *               Gopal Bhatt  @ Pennsylvania State Univeristy                  *
00007  * Version     : 2.0 (July 10, 2007)                                           *
00008  *-----------------------------------------------------------------------------*
00009  *                                                                             *
00010  * Part of initialization work is done in read_alloc, and the rest is done here*
00011  *                                                                             *
00012  * This code is free for users with research purpose only, if appropriate      *
00013  * citation is refered. However, there is no warranty in any format for this   *
00014  * product.                                                                    *
00015  *                                                                             *
00016  * For questions or comments, please contact the authors of the reference.     *
00017  * One who want to use it for other consideration may also contact Dr.Duffy    *
00018  * at cxd11@psu.edu.                                                           *
00019  *******************************************************************************/
00020 
00022 
00023 /* C Header Files */
00024 #include <stdio.h>
00025 #include <stdlib.h>
00026 #include <math.h>
00027 #include <string.h>
00028 
00029 /* SUNDIALS Header Files */
00030 #include "sundials_types.h"
00031 #include "nvector_serial.h"
00032 
00033 /* PIHM Header Files */
00034 #include "pihm.h"
00035 #include "calib.h"
00036 
00037 /* Calibration Parameters */
00038 realtype satD_CALIB;
00039 realtype br_CALIB;
00040 realtype poros_CALIB;
00041 realtype icsat_CALIB;
00042 realtype rivEle_CALIB;
00043 
00044 int lbool;    
00046 /*******************************************************************************
00047 *    Initialize Model & Control Data
00048 ********************************************************************************/
00049 void initialize(char *filename, Model_Data DS, Control_Data *CS, N_Vector CV_Y)
00051 
00056 {
00057       int i,j,counterMin, counterMax, MINCONST, domcounter;
00058       realtype a_x, a_y, b_x, b_y, c_x, c_y, MAXCONST;
00059       realtype a_zmin, a_zmax, b_zmin, b_zmax, c_zmin, c_zmax;
00060       realtype tempvalue;
00061       FILE *int_file;
00062       char *fn;
00063 
00064       realtype *zmin_cor;
00065 
00066     /* initializing calibration parameters    */
00067       satD_CALIB = setsatD_CALIB();
00068       br_CALIB   = setbr_CALIB();
00069       poros_CALIB = setporos_CALIB();
00070       icsat_CALIB = seticsat_CALIB();
00071       rivEle_CALIB = setrivEle_CALIB();
00072 
00073       zmin_cor=(realtype *)malloc(DS->NumEle*sizeof(realtype));
00074 
00075       printf("\nInitializing data structure ... ");
00076 
00077       for(i=0; i<DS->NumEle; i++)
00078       {
00079         a_x = DS->Node[DS->Ele[i].node[0]-1].x;          /* x coordinate of node 1       */
00080         b_x = DS->Node[DS->Ele[i].node[1]-1].x;          /* x coordinate of node 2       */
00081         c_x = DS->Node[DS->Ele[i].node[2]-1].x;          /* x coordinate of node 3       */
00082         a_y = DS->Node[DS->Ele[i].node[0]-1].y;          /* y coordinate of node 1       */
00083         b_y = DS->Node[DS->Ele[i].node[1]-1].y;          /* y coordinate of node 2       */
00084         c_y = DS->Node[DS->Ele[i].node[2]-1].y;          /* y coordinate of node 3       */
00085 
00086         a_zmin = DS->Node[DS->Ele[i].node[0]-1].zmin;    /* Bed elevation of node 1      */
00087         b_zmin = DS->Node[DS->Ele[i].node[1]-1].zmin;    /* Bed elevation of node 2      */
00088         c_zmin = DS->Node[DS->Ele[i].node[2]-1].zmin;    /* Bed elevation of node 3      */
00089         a_zmax = DS->Node[DS->Ele[i].node[0]-1].zmax;    /* Surface elevation of node 1  */
00090         b_zmax = DS->Node[DS->Ele[i].node[1]-1].zmax;    /* Surface elevation of node 2  */
00091         c_zmax = DS->Node[DS->Ele[i].node[2]-1].zmax;    /* Surface elevation of node 3  */
00092 
00093         /* Finding Lowest and Highest Elevation of an Element and Distance between them    */
00094         MINCONST =  10000000;
00095         MAXCONST = -1000000;
00096         for(j=0; j<3; j++)
00097         {
00098             if(DS->Node[DS->Ele[i].node[j]-1].zmin<MINCONST)
00099             {
00100                 MINCONST=DS->Node[DS->Ele[i].node[j]-1].zmin;
00101                 counterMin=j;
00102             }
00103             if(DS->Node[DS->Ele[i].node[j]-1].zmax>MAXCONST)
00104             {
00105                 MAXCONST=DS->Node[DS->Ele[i].node[j]-1].zmax;
00106                 counterMax=j;
00107             }
00108         }
00109         DS->Ele[i].NodeZmin= DS->Node[DS->Ele[i].node[counterMin]-1].zmin;            /*    Lowest Elevation      */
00110         DS->Ele[i].NodeZmax= DS->Node[DS->Ele[i].node[counterMax]-1].zmax;            /*    Highest Elevation     */
00111         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));
00112                                                               /*    Distance between Lowest & Highest Ele Nodes    */
00113         DS->Ele[i].area = 0.5*((b_x - a_x)*(c_y - a_y) - (b_y - a_y)*(c_x - a_x));    /*    Area of the Element    */
00114         //printf("\n%lf",DS->Ele[i].area);
00115         DS->Ele[i].zmax = (a_zmax + b_zmax + c_zmax)/3.0;                             /*    Mean Surface Elevation of an element    */
00116         DS->Ele[i].zmin = (a_zmin + b_zmin + c_zmin)/3.0;                             /*    Mean Bed Elevation of an element        */
00117 
00118          //DS->Ele[i].zmin =DS->Ele[i].zmax-br_CALIB;
00119 
00120 
00121         /* Calculate centroid of triangle */
00122         DS->Ele[i].x = (a_x + b_x + c_x)/3.0;
00123         DS->Ele[i].y = (a_y + b_y + c_y)/3.0;
00124 
00125 
00126         /*    Calculate Edge Lengths    */
00127         DS->Ele[i].edge[0] = pow((b_x - c_x), 2) + pow((b_y - c_y), 2);                /*    Length of Edge 1 of an element    */
00128         DS->Ele[i].edge[1] = pow((c_x - a_x), 2) + pow((c_y - a_y), 2);                /*    Length of Edge 2 of an element    */
00129         DS->Ele[i].edge[2] = pow((a_x - b_x), 2) + pow((a_y - b_y), 2);                /*    Length of Edge 3 of an element    */
00130 
00131 
00132         /* calculate circumcenter of triangle */
00133           /*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);
00134         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);
00135         */
00136         DS->Ele[i].edge[0] = sqrt(DS->Ele[i].edge[0]);
00137         DS->Ele[i].edge[1] = sqrt(DS->Ele[i].edge[1]);
00138         DS->Ele[i].edge[2] = sqrt(DS->Ele[i].edge[2]);
00139 
00140         /***************** Temporary ET inputs ***********************************/
00141         DS->Ele[i].windH = DS->WindH[DS->Ele[i].WindVel - 1]; //10;
00142         /*************************************************************************/
00143 
00144       }
00145       /*********************************************/
00146     /*    Optional :: Routine to Find Sinks         */
00147     /*********************************************/
00148     /*for(i=0; i<DS->NumEle; i++){
00149         lbool=0;
00150         for(j=0; j<3; j++){
00151             if(DS->Ele[i].nabr[j]>0){
00152                 if(DS->Ele[i].zmax+1.0>DS->Ele[DS->Ele[i].nabr[j]-1].zmax)
00153                     lbool=1;
00154             }
00155         }
00156         if(lbool==0){
00157             printf("%d is sink\n",i+1);
00158         }
00159     }
00160     //getchar();
00161     */
00162 
00163       /*    Memory allocation for flux terms     */
00164       DS->FluxSurf = (realtype **)malloc(DS->NumEle*sizeof(realtype));       /* Memory allocation for Surface flux                    */
00165       DS->FluxSub = (realtype **)malloc(DS->NumEle*sizeof(realtype));        /* Memory allocation for Sursurface flux                 */
00166       DS->FluxRiv = (realtype **)malloc(DS->NumRiv*sizeof(realtype));        /* Memory allocation for River Flux                      */
00167       DS->EleET = (realtype **)malloc(DS->NumEle*sizeof(realtype));          /* Memory allocation for Evapotranspiration              */
00168 
00169       for(i=0; i<DS->NumEle; i++)
00170       {
00171         DS->FluxSurf[i] = (realtype *)malloc(3*sizeof(realtype));            /* Memory allocation for Surface flux: Cont...           */
00172         DS->FluxSub[i] = (realtype *)malloc(3*sizeof(realtype));             /* Memory allocation for Subsurface flux: Cont...        */
00173         DS->EleET[i] = (realtype *)malloc(4*sizeof(realtype));               /* Memory allocation for Evapotranspiration: Cont...     */
00174       }
00175 
00176       for(i=0; i<DS->NumRiv; i++)
00177       {
00178         DS->FluxRiv[i] = (realtype *)malloc(6*sizeof(realtype));             /* Memory allocation for River flux: Cont...             */
00179       }
00180 
00181       DS->ElePrep = (realtype *)malloc(DS->NumEle*sizeof(realtype));         /* Memory allocation for Precipitation to the Element    */
00182       DS->EleVic = (realtype *)malloc(DS->NumEle*sizeof(realtype));          /* Memory allocation for Infiltration to the Element     */
00183       DS->Recharge = (realtype *)malloc(DS->NumEle*sizeof(realtype));        /* Memory allocation for Recharge to GW of a kernel      */
00184       DS->EleIS = (realtype *)malloc(DS->NumEle*sizeof(realtype));           /* Memory allocation for Interception Storage of Element */
00185       DS->EleISmax = (realtype *)malloc(DS->NumEle*sizeof(realtype));        /* Memory allocation for Maximum Interception Storage    */
00186       DS->EleSnow = (realtype *)malloc(DS->NumEle*sizeof(realtype));         /* Memory allocation for Snow accumulation of Element    */
00187       DS->EleTF = (realtype *)malloc(DS->NumEle*sizeof(realtype));           /* Memory allocation for Throughfall of an Element       */
00188       DS->Ele2IS = (realtype *)malloc(DS->NumEle*sizeof(realtype));          /* Memory allocation for Interception Storage Rate       */
00189       DS->EleNetPrep = (realtype *)malloc(DS->NumEle*sizeof(realtype));      /* Memory allocation for Net Precipitation to Element    */
00190 
00191       for(i=0; i<DS->NumEle; i++)
00192       {
00193         DS->Ele[i].Ksat = DS->Soil[(DS->Ele[i].soil-1)].Ksat;                /* Saturation Hydraulic Conductivity of an Element       */
00194         DS->Ele[i].Porosity = DS->Soil[(DS->Ele[i].soil-1)].SitaS -
00195                           DS->Soil[(DS->Ele[i].soil-1)].SitaR;               /* Porosity of an Element                                */
00196         DS->Ele[i].Porosity = poros_CALIB*DS->Ele[i].Porosity;               /* CALIBRATION                                           */
00197         DS->Ele[i].Alpha = DS->Soil[(DS->Ele[i].soil-1)].Alpha;              /* Soil Parameter Alpha of an Element                    */
00198         DS->Ele[i].Beta = DS->Soil[(DS->Ele[i].soil-1)].Beta;                /* Soil Parameter Beta of an Element                     */
00199         DS->Ele[i].Sf = DS->Soil[(DS->Ele[i].soil-1)].Sf;                    /* Slope Fricition Factor of an Element                  */
00200         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 */
00201                                                                              /* Root Zone Depth of the Element                        */
00202         DS->Ele[i].LAImax = DS->LandC[DS->Ele[i].LC-1].LAImax;               /* Maximum Leaf Area Index of the element                */
00203         DS->Ele[i].Rmin = DS->LandC[DS->Ele[i].LC-1].Rmin;                   /* Minimum Stomatal Resistance of the element            */
00204         DS->Ele[i].Rs_ref = DS->LandC[DS->Ele[i].LC-1].Rs_ref;               /* Reference ??        */
00205         DS->Ele[i].Albedo = DS->LandC[DS->Ele[i].LC-1].Albedo;               /* Albedo of an element                                  */
00206         DS->Ele[i].VegFrac = DS->LandC[DS->Ele[i].LC-1].VegFrac;             /* Vegetation Fraction/Cover of an Element               */
00207         DS->Ele[i].Rough = DS->LandC[DS->Ele[i].LC-1].Rough;                 /* Manning's Roughness Coefficient of an Element         */
00208       }
00209 
00210       for(i=0; i<DS->NumRiv; i++)
00211       {
00212         DS->Riv[i].x = (DS->Node[DS->Riv[i].FromNode-1].x +
00213                     DS->Node[DS->Riv[i].ToNode-1].x)/2;                      /* x-centroid of a river segment                         */
00214         DS->Riv[i].y = (DS->Node[DS->Riv[i].FromNode-1].y +
00215                     DS->Node[DS->Riv[i].ToNode-1].y)/2;                      /* y-centroid of a river segment                         */
00216         DS->Riv[i].zmax = (DS->Node[DS->Riv[i].FromNode-1].zmax + DS->Node[DS->Riv[i].ToNode-1].zmax)/2;
00217                                                                              /* Bank Elevation a river segment                        */
00218         DS->Riv[i].depth = DS->Riv_Shape[DS->Riv[i].shape-1].depth;          /* Depth of a river segment                              */
00219         DS->Riv[i].zmin = DS->Riv[i].zmax - DS->Riv[i].depth;                /* Bed Elevation of a river segment                      */
00220 
00221         DS->Riv[i].Length = sqrt(pow(DS->Node[DS->Riv[i].FromNode-1].x -
00222                                  DS->Node[DS->Riv[i].ToNode-1].x, 2) +
00223                              pow(DS->Node[DS->Riv[i].FromNode-1].y -
00224                                  DS->Node[DS->Riv[i].ToNode-1].y, 2));       /* Length of a River Segment                             */
00225       }
00226       /************************************************/
00227     /* Optional : Routine to see River Bed Slope    */
00228       /************************************************/
00229       /*for(i=0; i<DS->NumRiv; i++)
00230       {
00231         if(DS->Riv[i].down>0)
00232         {
00233             if(DS->Riv[i].zmin-DS->Riv[DS->Riv[i].down-1].zmin<=0.0)
00234             {
00235                 printf("\n%d\t%lf",i,DS->Riv[i].zmin-DS->Riv[DS->Riv[i].down-1].zmin);
00236                 //DS->Ele[DS->Riv[i].LeftEle-1].zmax,DS->Ele[DS->Riv[i].RightEle-1].zmax);
00237             }
00238         }
00239     }
00240     getchar();
00241     */
00242 
00243       /* Initialize State Variables */
00244       if (CS->int_type == 0)      /* relax cases */
00245       {
00246         for(i=0; i<DS->NumEle; i++)
00247         {
00248               DS->EleIS[i] = 0;                             /* Initialize Interception Storage         */
00249               NV_Ith_S(CV_Y, i) = 0;                        /* Initialize Surface State                */
00250               NV_Ith_S(CV_Y, i + DS->NumEle) = (1/DS->Ele[i].Alpha)*(1-exp(-DS->Ele[i].Alpha*(0.1+0.05)));
00251                                                             /* Initialize Unsaturated Zone State       */
00252               NV_Ith_S(CV_Y, i + 2*DS->NumEle) = DS->Ele[i].zmax - DS->Ele[i].zmin - 0.1;
00253                                                             /* Initialize Saturated Zone State         */
00254         }
00255 
00256         for(i=0; i<DS->NumRiv; i++)
00257         {
00258               NV_Ith_S(CV_Y, i + 3*DS->NumEle) = 0;         /* Initialize River State                  */
00259         }
00260       }
00261       /* type mode */
00262       else if (CS->int_type == 1)
00263       {
00264         if(DS->UnsatMode ==1)
00265         {
00266             for(i=0; i<DS->NumEle; i++)
00267             {
00268                   DS->EleIS[i] = DS->Ele_IC[i].interception;            /* Initialize Interception Storage        */
00269                   DS->EleSnow[i]=DS->Ele_IC[i].snow;                    /* Initialize Snow Storage                */
00270                   NV_Ith_S(CV_Y, i) = DS->Ele_IC[i].surf;               /* Initialize Surface State               */
00271                   NV_Ith_S(CV_Y, i + DS->NumEle) = DS->Ele_IC[i].sat;   /* Initialize SubSurface State            */
00272             }
00273 
00274             for(i=0; i<DS->NumRiv; i++)
00275             {
00276                   NV_Ith_S(CV_Y, i + 2*DS->NumEle) = DS->Riv_IC[DS->Riv[i].IC-1].value;
00277                                                                         /* Initialize River State                 */
00278             }
00279         }
00280         if(DS->UnsatMode ==2)
00281         {
00282             for(i=0; i<DS->NumEle; i++)
00283             {
00284                 /*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);*/
00285                   //DS->Ele_IC[i].sat=  (DS->Ele[i].zmax - DS->Ele[i].zmin)*icsat_CALIB;
00286                 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;
00287                   DS->EleIS[i] = DS->Ele_IC[i].interception;                 /* Initialize Interception Storage         */
00288                   NV_Ith_S(CV_Y, i) = DS->Ele_IC[i].surf;                    /* Initialize Surface State                */
00289 
00290                   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;
00291 
00292                   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 */
00293                                                                              /* Initialize Unsaturated Zone State    */
00294                   NV_Ith_S(CV_Y, i + 2*DS->NumEle) = DS->Ele_IC[i].sat;      /* Initialize Unsaturated Zone State    */
00295 
00296                 /* If Sat + UnSat > Bed/Soil Thickness    */
00297                   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))
00298                   {
00299                     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;
00300                        // 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));
00301                     if (NV_Ith_S(CV_Y, i + DS->NumEle) < 0)
00302                     {
00303                         NV_Ith_S(CV_Y, i + DS->NumEle) = 0;
00304                     }
00305                   }
00306             }
00307 
00308             for(i=0; i<DS->NumRiv; i++)
00309             {
00310                   NV_Ith_S(CV_Y, i + 3*DS->NumEle) = DS->Riv_IC[DS->Riv[i].IC-1].value;    /* Initialize River State      */
00311                   //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;
00312                   //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;
00313                   //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;
00314                   //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;
00315             }
00316         }
00317       }
00318       /* hot start mode */
00319       else if(CS->int_type == 2)
00320       {
00321         fn = (char *)malloc((strlen(filename)+4)*sizeof(char));
00322         strcpy(fn, filename);
00323         int_file = fopen(strcat(fn, ".int"), "r");
00324 
00325         if(int_file == NULL)
00326         {
00327               printf("\n  Fatal Error: %s.int is in use or does not exist!\n", filename);
00328               exit(1);
00329         }
00330         else
00331         {
00332               for(i=0; i<DS->NumEle; i++)
00333               {
00334                 fscanf(int_file, "%lf", &tempvalue);
00335                 if(tempvalue <= 0) {tempvalue = 0.01;}
00336                 NV_Ith_S(CV_Y, i + DS->NumEle) = tempvalue;
00337               }
00338 
00339               for(i=0; i<DS->NumEle; i++)
00340               {
00341                 fscanf(int_file, "%lf", &tempvalue);
00342                 if(tempvalue <= 0) {tempvalue = 0.01;}
00343                 if(tempvalue >= (DS->Ele[i].zmax - DS->Ele[i].zmin)) {tempvalue = (DS->Ele[i].zmax - DS->Ele[i].zmin) - 0.01;}
00344                 NV_Ith_S(CV_Y, i + 2*DS->NumEle) = tempvalue;
00345               }
00346 
00347               for(i=0; i<DS->NumEle; i++)
00348               {
00349                 DS->EleIS[i] = 0;
00350                 NV_Ith_S(CV_Y, i) = 0;
00351               }
00352 
00353               for(i=0; i<DS->NumRiv; i++)
00354               {
00355                 NV_Ith_S(CV_Y, i + 3*DS->NumEle) = 0;
00356               }
00357         }
00358         fclose(int_file);
00359       }
00360       /**************************************************/
00361       /* Adding routine to initialize state from a file */
00362       /* with initCondition/States at any given time    */
00363       /**************************************************/
00364       else if(CS->int_type == 3)
00365       {
00366         /* Open .init File for read    */
00367         fn = (char *)malloc((strlen(filename)+5)*sizeof(char));
00368         strcpy(fn, filename);
00369         int_file = fopen(strcat(fn, ".init"), "r");
00370 
00371         if(int_file == NULL)    /* Exit if couldn't open the file    */
00372         {
00373             printf("\n  Fatal Error: %s.int is in use or does not exist!\n", filename);
00374             exit(1);
00375         }
00376         else
00377         {
00378             fscanf(int_file, "%lf", &tempvalue);
00379             /* Check if the expected Start Time matches with that of the .init file    */
00380             if(tempvalue != CS->StartTime){
00381                 printf("\n  Fatal Error: Initial time in .init file does not match start time\n");
00382                 exit(1);
00383             }
00384 
00385             if(DS->UnsatMode == 1)
00386             {
00387                 for(i=0; i<DS->NumEle; i++)
00388                 {
00389                     fscanf(int_file, "%lf", &tempvalue);      /* Read Interception Storage State    */
00390                       DS->EleIS[i] = tempvalue;
00391                       fscanf(int_file, "%lf", &tempvalue);    /* Read Snow Storage State            */
00392                       DS->EleSnow[i]=tempvalue;
00393                       fscanf(int_file, "%lf", &tempvalue);    /* Read Surface Flow State            */
00394                       NV_Ith_S(CV_Y, i) = tempvalue;
00395                      /*NV_Ith_S(CV_Y, i + DS->NumEle) = DS->Ele_IC[DS->Ele[i].IC-1].unsat; */
00396                      fscanf(int_file, "%lf", &tempvalue);     /* Read SubSurface Flow State         */
00397                       NV_Ith_S(CV_Y, i + DS->NumEle) = tempvalue;
00398                 }
00399 
00400                 for(i=0; i<DS->NumRiv; i++)
00401                 {
00402                     fscanf(int_file, "%lf", &tempvalue);      /* Read River Flow State              */
00403                     NV_Ith_S(CV_Y, i + 2*DS->NumEle) = tempvalue;
00404                 }
00405             }
00406             if(DS->UnsatMode == 2)
00407             {
00408                 for(i=0; i<DS->NumEle; i++)
00409                 {
00410                     fscanf(int_file, "%lf", &tempvalue);      /* Read Interception Storage State    */
00411                     DS->EleIS[i] = tempvalue;
00412                     fscanf(int_file, "%lf", &tempvalue);      /* Read Snow Storage State            */
00413                     DS->EleSnow[i]=tempvalue;
00414                     fscanf(int_file, "%lf", &tempvalue);      /* Read Surface Flow State            */
00415                     NV_Ith_S(CV_Y, i) = tempvalue;
00416                     fscanf(int_file, "%lf", &tempvalue);      /* Read Unsaturated Zone State        */
00417                     NV_Ith_S(CV_Y, i + DS->NumEle) = tempvalue;
00418                     fscanf(int_file, "%lf", &tempvalue);      /* Read Saturated Zone State          */
00419                     NV_Ith_S(CV_Y, i + 2*DS->NumEle) = tempvalue;
00420                 }
00421 
00422                 for(i=0; i<DS->NumRiv; i++)
00423                 {
00424                     fscanf(int_file, "%lf", &tempvalue);       /* Read River Flow State             */
00425                     NV_Ith_S(CV_Y, i + 3*DS->NumEle) = tempvalue;
00426                 }
00427             }
00428         }
00429     }
00430 
00431     printf("done.\n");
00432 }
00433