read_alloc.c

Go to the documentation of this file.
00001 /*******************************************************************************
00002  * File        : read_alloc.c                                                  *
00003  * Function    : read Input Files and allocate memorey for PIHM                *
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  * In this file, 7 input files are read in and dynamaticlly allocate memory    *
00011  * based on control parameters in the input files. Please refer to format.pdf  *
00012  * for details about the format of all 7 input files. All files has the same   *
00013  * file prefix provided in the command line when starting PIHM. For example,   *
00014  * "example" will incur reading in example.mesh, example.soil, etc.            *
00015  *                                                                             *
00016  * This code is free for users with research purpose only, if appropriate      *
00017  * citation is refered. However, there is no warranty in any format for this   *
00018  * product.                                                                    *
00019  *                                                                             *
00020  * For questions or comments, please contact the authors of the reference.     *
00021  * One who want to use it for other consideration may also contact Dr.Duffy    *
00022  * at cxd11@psu.edu.                                                           *
00023  *******************************************************************************/
00024 
00026 
00027 /* C Header Files */
00028 #include <stdio.h>
00029 #include <stdlib.h>
00030 #include <math.h>
00031 #include <string.h>
00032 
00033 /* SUNDIALS Header Files */
00034 #include "sundials_types.h"
00035 #include "pihm.h"
00036 #include "calib.h"
00037 
00038 
00039 
00040 
00041 /***************************************************************
00042     Function reads all the input files
00043 ****************************************************************/
00044 void read_alloc(char *filename, Model_Data DS, Control_Data *CS)
00046 
00050 {
00051     int i, j;
00052     int tempindex;
00053 
00054     int NumTout; 
00055     char *fn[8];
00056     char tempchar[5];
00057 
00058     FILE *mesh_file;                        /*    Pointer to .mesh file    */
00059     FILE *att_file;                         /*    Pointer to .att  file    */
00060     FILE *forc_file;                        /*    Pointer to .forc file    */
00061     FILE *ibc_file;                         /*    Pointer to .ibc  file    */
00062     FILE *soil_file;                        /*    Pointer to .soil file    */
00063     FILE *lc_file;                          /*    Pointer to .lc     file  */
00064     FILE *para_file;                        /*    Pointer to .para file    */
00065     FILE *riv_file;                         /*    Pointer to .riv  file    */
00066 
00067     /* CALIBRATION Parameters    */
00068         realtype roughEle_CALIB;
00069         realtype roughRiv_CALIB;
00070         realtype rivCoeff_CALIB;
00071         realtype rivDepth_CALIB;
00072         realtype alpha_CALIB;
00073         realtype set_MP;
00074         realtype lai_CALIB;
00075         realtype vegfrac_CALIB;
00076     realtype albedo_CALIB;
00077 
00078     /*    Initialize CALIBRATION parameters    */
00079     roughEle_CALIB=setroughEle_CALIB();
00080     roughRiv_CALIB=setroughRiv_CALIB();
00081     rivCoeff_CALIB=setrivCoeff_CALIB();
00082     rivDepth_CALIB=setrivDepth_CALIB();
00083     alpha_CALIB=setalpha_CALIB();
00084     set_MP=setset_MP();
00085     lai_CALIB=setlai_CALIB();
00086     vegfrac_CALIB=setvegfrac_CALIB();
00087     albedo_CALIB=setalbedo_CALIB();
00088 
00089 
00090     printf("\nStart reading in input files ... \n");
00091 
00092     /****************************************/
00093     /*========== open *.mesh file ==========*/
00094     /****************************************/
00095     printf("\n  1) reading %s.mesh ... ", filename);
00096     fn[0] = (char *)malloc((strlen(filename)+5)*sizeof(char));
00097     strcpy(fn[0], filename);
00098     mesh_file = fopen(strcat(fn[0], ".mesh"), "r");
00099 
00100     if(mesh_file == NULL)
00101     {
00102         printf("\n  Fatal Error: %s.mesh is in use or does not exist!\n", filename);
00103         exit(1);
00104     }
00105 
00106     /* start reading mesh_file */
00107     fscanf(mesh_file,"%d %d", &DS->NumEle, &DS->NumNode);
00108 
00109     DS->Ele = (element *)malloc(DS->NumEle*sizeof(element));          /*    Allocate memory for Number of Elements     */
00110     DS->Node = (nodes *)malloc(DS->NumNode*sizeof(nodes));            /*    Allocate memory for Number of Nodes        */
00111 
00112     /* read in elements information */
00113     for (i=0; i<DS->NumEle; i++)
00114     {
00115         fscanf(mesh_file, "%d", &(DS->Ele[i].index));
00116         fscanf(mesh_file, "%d %d %d", &(DS->Ele[i].node[0]), &(DS->Ele[i].node[1]), &(DS->Ele[i].node[2]));
00117         fscanf(mesh_file, "%d %d %d", &(DS->Ele[i].nabr[0]), &(DS->Ele[i].nabr[1]), &(DS->Ele[i].nabr[2]));
00118     }
00119 
00120     /* read in nodes information */
00121     for (i=0; i<DS->NumNode; i++)
00122     {
00123         fscanf(mesh_file, "%d", &(DS->Node[i].index));
00124         fscanf(mesh_file, "%lf %lf", &(DS->Node[i].x), &(DS->Node[i].y));
00125         fscanf(mesh_file, "%lf %lf", &(DS->Node[i].zmin),&(DS->Node[i].zmax));
00126     }
00127 
00128     printf("done.\n");
00129 
00130     fclose(mesh_file);
00131     /* finish reading mesh_files */
00132 
00133     /***************************************/
00134     /*========== open *.att file ==========*/
00135     /***************************************/
00136     printf("\n  2) reading %s.att  ... ", filename);
00137     fn[1] = (char *)malloc((strlen(filename)+4)*sizeof(char));
00138     strcpy(fn[1], filename);
00139     att_file = fopen(strcat(fn[1], ".att"), "r");
00140 
00141     if(att_file == NULL)
00142     {
00143         printf("\n  Fatal Error: %s.att is in use or does not exist!\n", filename);
00144         exit(1);
00145     }
00146 
00147     /* start reading att_file */
00148     DS->Ele_IC = (element_IC *)malloc(DS->NumEle*sizeof(element_IC));
00149 
00150     for (i=0; i<DS->NumEle; i++)
00151     {
00152         fscanf(att_file, "%d", &(tempindex));
00153         fscanf(att_file, "%d %d", &(DS->Ele[i].soil), &(DS->Ele[i].LC));            /*    Read Soil and Land Cover Type               */
00154         fscanf(att_file, "%lf %lf %lf %lf %lf",&(DS->Ele_IC[i].interception),&(DS->Ele_IC[i].snow),&(DS->Ele_IC[i].surf),&(DS->Ele_IC[i].unsat),&(DS->Ele_IC[i].sat));
00155                                                                                     /*    Read Initial States                         */
00156         fscanf(att_file, "%d", &(DS->Ele[i].BC));                                   /*    Read Boundary Condition                     */
00157         //DS->Ele[i].BC=0;
00158         fscanf(att_file, "%d %d", &(DS->Ele[i].prep), &(DS->Ele[i].temp));          /*    Read Precipitation & Temperature Class      */
00159         fscanf(att_file, "%d %d", &(DS->Ele[i].humidity), &(DS->Ele[i].WindVel));   /*    Read Rel. Humidity & Wind Velocity Class    */
00160         fscanf(att_file, "%d %d", &(DS->Ele[i].Rn), &(DS->Ele[i].G));               /*    Read Solar Radiation Class    + Dummy       */
00161         fscanf(att_file, "%d %d", &(DS->Ele[i].pressure), &(DS->Ele[i].source));    /*    Read Pressure Class and Source/Sinks        */
00162     }
00163 
00164     printf("done.\n");
00165 
00166     fclose(att_file);
00167     /* finish reading mesh_files */
00168 
00169     /****************************************/
00170     /*========== open *.soil file ==========*/
00171     /****************************************/
00172     printf("\n  3) reading %s.soil ... ", filename);
00173     fn[2] = (char *)malloc((strlen(filename)+5)*sizeof(char));
00174     strcpy(fn[2], filename);
00175     soil_file = fopen(strcat(fn[2], ".soil"), "r");
00176 
00177     if(soil_file == NULL)
00178     {
00179         printf("\n  Fatal Error: %s.soil is in use or does not exist!\n", filename);
00180         exit(1);
00181     }
00182 
00183     /* start reading soil_file */
00184     fscanf(soil_file, "%d", &DS->NumSoil);
00185 
00186     DS->Soil = (soils *)malloc(DS->NumSoil*sizeof(soils));
00187 
00188     for (i=0; i<DS->NumSoil; i++)
00189     {
00190         fscanf(soil_file, "%d", &(DS->Soil[i].index));
00191         fscanf(soil_file, "%lf", &(DS->Soil[i].Ksat));                                /*    Read Saturated Hydraulic Conductivity    */
00192         fscanf(soil_file, "%lf %lf", &(DS->Soil[i].SitaS), &(DS->Soil[i].SitaR));     /*    Read Porosity and Residual Porosity      */
00193         fscanf(soil_file, "%lf %lf", &(DS->Soil[i].Alpha), &(DS->Soil[i].Beta));      /*    Read Soil Parameters Alpha & Beta        */
00194         fscanf(soil_file, "%d %lf %lf", &(DS->Soil[i].Macropore), &(DS->Soil[i].base), &(DS->Soil[i].gama));
00195                                                                                       /*    Read Macropore (0/1) Base & Gamma        */
00196         fscanf(soil_file, "%lf", &(DS->Soil[i].Sf));                                  /*    Read Soil Friction Slope                 */
00197         fscanf(soil_file, "%lf", &(DS->Soil[i].RzD));                                 /*    Read Root Zone Depth                     */
00198         fscanf(soil_file, "%d", &(DS->Soil[i].Inf));
00199 
00200         /*     CALIBRATION    */
00201         DS->Soil[i].Alpha=alpha_CALIB*DS->Soil[i].Alpha;
00202         DS->Soil[i].Macropore=set_MP;
00203     }
00204 
00205     fscanf(soil_file, "%d", &DS->NumInc);
00206 
00207     DS->TSD_Inc = (TSD *)malloc(DS->NumInc*sizeof(TSD));
00208 
00209     for(i=0; i<DS->NumInc; i++)
00210     {
00211         fscanf(soil_file, "%s %d %d", DS->TSD_Inc[i].name, &DS->TSD_Inc[i].index,
00212                                   &DS->TSD_Inc[i].length);
00213 
00214         DS->TSD_Inc[i].TS = (realtype **)malloc(DS->TSD_Inc[i].length*sizeof(realtype));
00215         for(j=0; j<DS->TSD_Inc[i].length; j++)
00216         {
00217             DS->TSD_Inc[i].TS[j] = (realtype *)malloc(2*sizeof(realtype));
00218         }
00219 
00220         for(j=0; j<DS->TSD_Inc[i].length; j++)
00221         {
00222             fscanf(soil_file, "%lf %lf", &DS->TSD_Inc[i].TS[j][0],
00223                                    &DS->TSD_Inc[i].TS[j][1]);
00224         }
00225     }
00226 
00227     fclose(soil_file);
00228     printf("done.\n");
00229     /* Finish reading soil_file */
00230 
00231 
00232     /**************************************/
00233     /*========== open *.lc file ==========*/
00234     /**************************************/
00235     printf("\n  3) reading %s.lc ... ", filename);
00236     fn[3] = (char *)malloc((strlen(filename)+5)*sizeof(char));
00237     strcpy(fn[3], filename);
00238     lc_file = fopen(strcat(fn[3], ".lc"), "r");
00239 
00240     if(lc_file == NULL)
00241     {
00242         printf("\n  Fatal Error: %s.land cover is in use or does not exist!\n", filename);
00243         exit(1);
00244     }
00245 
00246     /* start reading land cover file */
00247     fscanf(lc_file, "%d", &DS->NumLC);                                                    /*    Number of Land Cover Classes             */
00248 
00249     DS->LandC = (LC *)malloc(DS->NumLC*sizeof(LC));
00250 
00251     for (i=0; i<DS->NumLC; i++)
00252     {
00253         fscanf(lc_file, "%d", &(DS->LandC[i].index));
00254         fscanf(lc_file, "%lf", &(DS->LandC[i].LAImax));                                   /*    Read Max LAI                             */
00255         fscanf(lc_file, "%lf %lf", &(DS->LandC[i].Rmin), &(DS->LandC[i].Rs_ref));         /*    Read Min Stomal Resistance & Reference   */
00256         fscanf(lc_file, "%lf %lf", &(DS->LandC[i].Albedo), &(DS->LandC[i].VegFrac));      /*    Read Albedo and Vegitation Fraction      */
00257         fscanf(lc_file, "%lf", &(DS->LandC[i].Rough));                                    /*    Read Manning's Roughnes Coefficient      */
00258 
00259         /*    CALIBRATION    */
00260         DS->LandC[i].VegFrac=vegfrac_CALIB*DS->LandC[i].VegFrac;
00261         DS->LandC[i].Albedo=albedo_CALIB*DS->LandC[i].Albedo>1.0?1.0:albedo_CALIB*DS->LandC[i].Albedo;
00262         DS->LandC[i].Rough=roughEle_CALIB*DS->LandC[i].Rough;
00263     }
00264 
00265     fclose(lc_file);
00266     printf("done.\n");
00267     /* Finish reading land cover file */
00268 
00269     /***************************************/
00270     /*========== open *.riv file ==========*/
00271     /***************************************/
00272     printf("\n  4) reading %s.riv  ... ", filename);
00273     fn[4] = (char *)malloc((strlen(filename)+4)*sizeof(char));
00274     strcpy(fn[4], filename);
00275     riv_file =  fopen(strcat(fn[4], ".riv"), "r");
00276 
00277     if(riv_file == NULL)
00278     {
00279         printf("\n  Fatal Error: %s.riv is in use or does not exist!\n", filename);
00280         exit(1);
00281     }
00282 
00283     /* start reading .riv File */
00284     fscanf(riv_file, "%d", &DS->NumRiv);
00285 
00286     DS->Riv = (river_segment *)malloc(DS->NumRiv*sizeof(river_segment));
00287     DS->Riv_IC = (river_IC *)malloc(DS->NumRiv*sizeof(river_IC));
00288 
00289     for (i=0; i<DS->NumRiv; i++)
00290     {
00291         fscanf(riv_file, "%d", &(DS->Riv[i].index));
00292         fscanf(riv_file, "%d %d", &(DS->Riv[i].FromNode), &(DS->Riv[i].ToNode));         /*    Read From and To Nodes #              */
00293         fscanf(riv_file, "%d", &(DS->Riv[i].down));                                      /*    Read Down Segment                     */
00294         fscanf(riv_file, "%d %d", &(DS->Riv[i].LeftEle), &(DS->Riv[i].RightEle));        /*    Read Left and Right Elements #        */
00295         fscanf(riv_file, "%d %d", &(DS->Riv[i].shape), &(DS->Riv[i].material));          /*    Read Shape and Material Type          */
00296         fscanf(riv_file, "%d %d", &(DS->Riv[i].IC), &(DS->Riv[i].BC));                   /*    Read Initial and Boundary Condition   */
00297         fscanf(riv_file, "%d", &(DS->Riv[i].reservoir));                                 /*    Read Reservoir Type                   */
00298     }
00299 
00300     fscanf(riv_file, "%s %d", tempchar, &DS->NumRivShape);                               /*    Read Number of Shape Types            */
00301     DS->Riv_Shape = (river_shape *)malloc(DS->NumRivShape*sizeof(river_shape));          /*    Allocate Memory for Shape Types       */
00302 
00303     for (i=0; i<DS->NumRivShape; i++)
00304     {
00305         fscanf(riv_file, "%d %lf", &DS->Riv_Shape[i].index, &DS->Riv_Shape[i].width);    /*    Read Shape # and Width                */
00306         fscanf(riv_file, "%lf %lf", &DS->Riv_Shape[i].depth, &DS->Riv_Shape[i].bed);     /*    Read Depth and Bed Elevation          */
00307         fscanf(riv_file, "%d %lf",&DS->Riv_Shape[i].interpOrd,&DS->Riv_Shape[i].coeff);  /*    Read Interpolation Order & Coeff      */
00308 
00309         /*    CALIBRATION    */
00310         DS->Riv_Shape[i].depth = rivDepth_CALIB*DS->Riv_Shape[i].depth;
00311         DS->Riv_Shape[i].coeff = rivCoeff_CALIB*DS->Riv_Shape[i].coeff;
00312     }
00313 
00314     fscanf(riv_file, "%s %d", tempchar, &DS->NumRivMaterial);                            /* Read Number of Material Types            */
00315     DS->Riv_Mat = (river_material *)malloc(DS->NumRivMaterial*sizeof(river_material));   /* Allocate Memory for Material Types       */
00316 
00317     for (i=0; i<DS->NumRivMaterial; i++)
00318     {
00319         fscanf(riv_file, "%d %lf %lf %lf", &DS->Riv_Mat[i].index, &DS->Riv_Mat[i].Rough, &DS->Riv_Mat[i].Cwr, &DS->Riv_Mat[i].Sf);
00320                                                               /* Read Roughness, Coeff. of Discharge and Manning's Roughness Coeff   */
00321 
00322         /*    CALIBRATION    */
00323         DS->Riv_Mat[i].Rough=roughRiv_CALIB*DS->Riv_Mat[i].Rough;
00324     }
00325 
00326     fscanf(riv_file, "%s %d", tempchar, &DS->NumRivIC);                                  /* Read Number of Initial Condition Types   */
00327     DS->Riv_IC = (river_IC *)malloc(DS->NumRivIC*sizeof(river_IC));                      /* Allocate Memory for Initial Cond. Types  */
00328 
00329     for (i=0; i<DS->NumRivIC; i++)
00330     {
00331         fscanf(riv_file, "%d %lf", &DS->Riv_IC[i].index, &DS->Riv_IC[i].value);          /* Read Initial Condition Value             */
00332     }
00333 
00334     fscanf(riv_file, "%s %d", tempchar, &DS->NumRivBC);                                  /* Read Number of Boundary Conditions       */
00335     DS->TSD_Riv = (TSD *)malloc(DS->NumRivBC*sizeof(TSD));                               /* Allocated memory for Boundary Conditions */
00336 
00337     for(i=0; i<DS->NumRivBC; i++)
00338     {
00339         fscanf(riv_file, "%s %d %d", DS->TSD_Riv[i].name, &DS->TSD_Riv[i].index, &DS->TSD_Riv[i].length);
00340                                                                                          /* Read Boundary Cond. TimeSeries Lenght    */
00341 
00342         DS->TSD_Riv[i].TS = (realtype **)malloc((DS->TSD_Riv[i].length)*sizeof(realtype));    /* Allocate Memory for TS              */
00343         for(j=0; j<DS->TSD_Riv[i].length; j++)
00344         {
00345             DS->TSD_Riv[i].TS[j] = (realtype *)malloc(2*sizeof(realtype));                    /* Allocate Memory for TS Cont...      */
00346         }
00347 
00348         for(j=0; j<DS->TSD_Riv[i].length; j++)
00349         {
00350             fscanf(riv_file, "%lf %lf", &DS->TSD_Riv[i].TS[j][0], &DS->TSD_Riv[i].TS[j][1]);  /* Read Time & BC Value                */
00351         }
00352     }
00353 
00354     /* read in reservoir information    */
00355     fscanf(riv_file, "%s %d", tempchar, &DS->NumRes);
00356     if(DS->NumRes > 0)
00357     {
00358         /* YET TO BE IMPLIMENTED */
00359     }
00360     fclose(riv_file);
00361     printf("done.\n");
00362     /* Finish reading .riv File */
00363 
00364     /****************************************/
00365     /*========== open *.forc file ==========*/
00366     /****************************************/
00367     printf("\n  5) reading %s.forc ... ", filename);
00368     fn[5] = (char *)malloc((strlen(filename)+5)*sizeof(char));
00369     strcpy(fn[5], filename);
00370     forc_file = fopen(strcat(fn[5], ".forc"), "r");
00371 
00372     if(forc_file == NULL)
00373     {
00374         printf("\n  Fatal Error: %s.forc is in use or does not exist!\n", filename);
00375         exit(1);
00376     }
00377 
00378     /* start reading .forc File */
00379     fscanf(forc_file, "%d %d", &DS->NumPrep, &DS->NumTemp);                /*    Read Number of Precipitaion & Temperature TS     */
00380     fscanf(forc_file, "%d %d", &DS->NumHumidity, &DS->NumWindVel);         /*     Read Number of Rel. Humidity & Wind Vel. TS     */
00381     fscanf(forc_file, "%d %d", &DS->NumRn, &DS->NumG);                     /*    Read Number of Solar Radiation & Dummy           */
00382     fscanf(forc_file, "%d %d", &DS->NumP, &DS->NumLC);                     /*    Read Number of Pressure and LAI TS               */
00383     fscanf(forc_file, "%d", &DS->NumMeltF);                                /*     Read Number of Melting Factor TS                */
00384     fscanf(forc_file, "%d", &DS->NumSource);                               /*    Read Number of Sources/Sinks TS                  */
00385 
00386     DS->TSD_Prep = (TSD *)malloc(DS->NumPrep*sizeof(TSD));                 /*    Allocate Memory for Precipitation Time Series    */
00387     DS->TSD_Temp = (TSD *)malloc(DS->NumTemp*sizeof(TSD));                 /*    Allocate Memory for Temperature Time Series      */
00388     DS->TSD_Humidity = (TSD *)malloc(DS->NumHumidity*sizeof(TSD));         /*    Allocate Memory for Rel. Humidity Time Series    */
00389     DS->TSD_WindVel = (TSD *)malloc(DS->NumWindVel*sizeof(TSD));           /*    Allocate Memory for Wind Velocity Time Series    */
00390     DS->TSD_Rn = (TSD *)malloc(DS->NumRn*sizeof(TSD));                     /*    Allocate Memory for Solar Radiation Time Series  */
00391     DS->TSD_G = (TSD *)malloc(DS->NumG*sizeof(TSD));                       /*    Allocate Memory for DUMMY Time Series            */
00392     DS->TSD_Pressure = (TSD *)malloc(DS->NumP*sizeof(TSD));                /*    Allocate Memory for Pressure Time Series         */
00393     DS->TSD_LAI = (TSD *)malloc(DS->NumLC*sizeof(TSD));                    /*    Allocate Memory for Leaf Area Index Time Series  */
00394     DS->TSD_DH = (TSD *)malloc(DS->NumLC*sizeof(TSD));                     /*    Allocate Memory for Stomal Resist Time Series    */
00395     DS->TSD_MeltF = (TSD *)malloc(DS->NumMeltF*sizeof(TSD));               /*    Allocate Memory for Melting Factor Time Series   */
00396     DS->TSD_Source = (TSD *)malloc(DS->NumSource*sizeof(TSD));             /*    Allocate Memory for Sources/Sinks Time Series    */
00397 
00398     DS->SIFactor = (realtype *)malloc(DS->NumLC*sizeof(realtype));         /*    Allocate Memory for SI Factor Time Series        */
00399     DS->WindH = (realtype *)malloc(DS->NumWindVel*sizeof(realtype));       /*    Allocate Memory for WindVel. Height Time Series  */
00400 
00401     /*    Read All the Precipitation Time Series    */
00402     for(i=0; i<DS->NumPrep; i++)
00403     {
00404         fscanf(forc_file, "%s %d %d", DS->TSD_Prep[i].name, &DS->TSD_Prep[i].index, &DS->TSD_Prep[i].length);
00405 
00406         DS->TSD_Prep[i].TS = (realtype **)malloc((DS->TSD_Prep[i].length)*sizeof(realtype));
00407 
00408         for(j=0; j<DS->TSD_Prep[i].length; j++)
00409         {
00410             DS->TSD_Prep[i].TS[j] = (realtype *)malloc(2*sizeof(realtype));
00411         }
00412 
00413         for(j=0; j<DS->TSD_Prep[i].length; j++)
00414         {
00415             fscanf(forc_file, "%lf %lf", &DS->TSD_Prep[i].TS[j][0], &DS->TSD_Prep[i].TS[j][1]);
00416         }
00417         DS->TSD_Prep[i].iCounter=0;
00418     }
00419 
00420     /*    Read All the Temperature Time Series    */
00421     for(i=0; i<DS->NumTemp; i++)
00422     {
00423         fscanf(forc_file, "%s %d %d", DS->TSD_Temp[i].name, &DS->TSD_Temp[i].index, &DS->TSD_Temp[i].length);
00424 
00425         DS->TSD_Temp[i].TS = (realtype **)malloc((DS->TSD_Temp[i].length)*sizeof(realtype));
00426 
00427         for(j=0; j<DS->TSD_Temp[i].length; j++)
00428         {
00429             DS->TSD_Temp[i].TS[j] = (realtype *)malloc(2*sizeof(realtype));
00430         }
00431 
00432         for(j=0; j<DS->TSD_Temp[i].length; j++)
00433         {
00434             fscanf(forc_file, "%lf %lf", &DS->TSD_Temp[i].TS[j][0], &DS->TSD_Temp[i].TS[j][1]);
00435         }
00436         DS->TSD_Temp[i].iCounter=0;
00437     }
00438 
00439     /*    Read All the Rel. Humidity Time Series    */
00440     for(i=0; i<DS->NumHumidity; i++)
00441     {
00442         fscanf(forc_file, "%s %d %d", DS->TSD_Humidity[i].name, &DS->TSD_Humidity[i].index, &DS->TSD_Humidity[i].length);
00443 
00444         DS->TSD_Humidity[i].TS = (realtype **)malloc((DS->TSD_Humidity[i].length)*sizeof(realtype));
00445 
00446         for(j=0; j<DS->TSD_Humidity[i].length; j++)
00447         {
00448             DS->TSD_Humidity[i].TS[j] = (realtype *)malloc(2*sizeof(realtype));
00449         }
00450 
00451         for(j=0; j<DS->TSD_Humidity[i].length; j++)
00452         {
00453             fscanf(forc_file, "%lf %lf", &DS->TSD_Humidity[i].TS[j][0], &DS->TSD_Humidity[i].TS[j][1]);
00454             DS->TSD_Humidity[i].TS[j][1] = (DS->TSD_Humidity[i].TS[j][1] > 1.0?1.0:DS->TSD_Humidity[i].TS[j][1]);
00455         }
00456         DS->TSD_Humidity[i].iCounter=0;
00457     }
00458 
00459     /*    Read All the Wind Velocity Time Series    */
00460     for(i=0; i<DS->NumWindVel; i++)
00461     {
00462         fscanf(forc_file, "%s %d %d %lf", DS->TSD_WindVel[i].name, &DS->TSD_WindVel[i].index, &DS->TSD_WindVel[i].length, &DS->WindH[i]);
00463 
00464         DS->TSD_WindVel[i].TS = (realtype **)malloc((DS->TSD_WindVel[i].length)*sizeof(realtype));
00465 
00466         for(j=0; j<DS->TSD_WindVel[i].length; j++)
00467         {
00468             DS->TSD_WindVel[i].TS[j] = (realtype *)malloc(2*sizeof(realtype));
00469         }
00470 
00471         for(j=0; j<DS->TSD_WindVel[i].length; j++)
00472         {
00473             fscanf(forc_file, "%lf %lf", &DS->TSD_WindVel[i].TS[j][0], &DS->TSD_WindVel[i].TS[j][1]);
00474         }
00475         DS->TSD_WindVel[i].iCounter=0;
00476     }
00477 
00478     /*    Read All the Solar Radiation Time Series    */
00479     for(i=0; i<DS->NumRn; i++)
00480     {
00481         fscanf(forc_file, "%s %d %d", DS->TSD_Rn[i].name, &DS->TSD_Rn[i].index, &DS->TSD_Rn[i].length);
00482 
00483         DS->TSD_Rn[i].TS = (realtype **)malloc((DS->TSD_Rn[i].length)*sizeof(realtype));
00484 
00485         for(j=0; j<DS->TSD_Rn[i].length; j++)
00486         {
00487             DS->TSD_Rn[i].TS[j] = (realtype *)malloc(2*sizeof(realtype));
00488         }
00489 
00490         for(j=0; j<DS->TSD_Rn[i].length; j++)
00491         {
00492             fscanf(forc_file, "%lf %lf", &DS->TSD_Rn[i].TS[j][0], &DS->TSD_Rn[i].TS[j][1]);
00493         }
00494         DS->TSD_Rn[i].iCounter=0;
00495     }
00496 
00497     /*    Read All the DUMMY Time Series    */
00498     for(i=0; i<DS->NumG; i++)
00499     {
00500         fscanf(forc_file, "%s %d %d", DS->TSD_G[i].name, &DS->TSD_G[i].index, &DS->TSD_G[i].length);
00501         DS->TSD_G[i].TS = (realtype **)malloc((DS->TSD_G[i].length)*sizeof(realtype));
00502         for(j=0; j<DS->TSD_G[i].length; j++)
00503         {
00504             DS->TSD_G[i].TS[j] = (realtype *)malloc(2*sizeof(realtype));
00505         }
00506 
00507         for(j=0; j<DS->TSD_G[i].length; j++)
00508         {
00509             fscanf(forc_file, "%lf %lf", &DS->TSD_G[i].TS[j][0], &DS->TSD_G[i].TS[j][1]);
00510         }
00511         DS->TSD_G[i].iCounter=0;
00512     }
00513 
00514     /*    Read All the Vapor Pressure Time Series    */
00515     for(i=0; i<DS->NumP; i++)
00516     {
00517         fscanf(forc_file, "%s %d %d", DS->TSD_Pressure[i].name, &DS->TSD_Pressure[i].index, &DS->TSD_Pressure[i].length);
00518         DS->TSD_Pressure[i].TS = (realtype **)malloc((DS->TSD_Pressure[i].length)*sizeof(realtype));
00519         for(j=0; j<DS->TSD_Pressure[i].length; j++)
00520         {
00521             DS->TSD_Pressure[i].TS[j] = (realtype *)malloc(2*sizeof(realtype));
00522         }
00523 
00524         for(j=0; j<DS->TSD_Pressure[i].length; j++)
00525         {
00526             fscanf(forc_file, "%lf %lf", &DS->TSD_Pressure[i].TS[j][0], &DS->TSD_Pressure[i].TS[j][1]);
00527         }
00528         DS->TSD_Pressure[i].iCounter=0;
00529     }
00530 
00531     /*    Read All the LAI Time Series    */
00532     for(i=0; i<DS->NumLC; i++)
00533     {
00534         fscanf(forc_file, "%s %d %d %lf", DS->TSD_LAI[i].name, &DS->TSD_LAI[i].index, &DS->TSD_LAI[i].length, &DS->SIFactor[i]);
00535         DS->TSD_LAI[i].TS = (realtype **)malloc((DS->TSD_LAI[i].length)*sizeof(realtype));
00536 
00537         for(j=0; j<DS->TSD_LAI[i].length; j++)
00538         {
00539             DS->TSD_LAI[i].TS[j] = (realtype *)malloc(2*sizeof(realtype));
00540         }
00541 
00542         for(j=0; j<DS->TSD_LAI[i].length; j++)
00543         {
00544             fscanf(forc_file, "%lf %lf", &DS->TSD_LAI[i].TS[j][0], &DS->TSD_LAI[i].TS[j][1]);
00545 
00546             DS->TSD_LAI[i].TS[j][1]=lai_CALIB*DS->TSD_LAI[i].TS[j][1];
00547         }
00548         DS->TSD_LAI[i].iCounter=0;
00549     }
00550 
00551     /*    Read All the Displacement Height Time Series    */
00552     for(i=0; i<DS->NumLC; i++)
00553     {
00554         fscanf(forc_file, "%s %d %d", DS->TSD_DH[i].name, &DS->TSD_DH[i].index, &DS->TSD_DH[i].length);
00555 
00556         DS->TSD_DH[i].TS = (realtype **)malloc((DS->TSD_DH[i].length)*sizeof(realtype));
00557 
00558         for(j=0; j<DS->TSD_DH[i].length; j++)
00559         {
00560             DS->TSD_DH[i].TS[j] = (realtype *)malloc(2*sizeof(realtype));
00561         }
00562         for(j=0; j<DS->TSD_DH[i].length; j++)
00563         {
00564             fscanf(forc_file, "%lf %lf", &DS->TSD_DH[i].TS[j][0], &DS->TSD_DH[i].TS[j][1]);
00565         }
00566         DS->TSD_DH[i].iCounter=0;
00567     }
00568 
00569     /*    Read All the Melting Factor Time Series    */
00570     for(i=0; i<DS->NumMeltF; i++)
00571     {
00572         fscanf(forc_file, "%s %d %d", DS->TSD_MeltF[i].name, &DS->TSD_MeltF[i].index, &DS->TSD_MeltF[i].length);
00573         DS->TSD_MeltF[i].TS = (realtype **)malloc((DS->TSD_MeltF[i].length)*sizeof(realtype));
00574 
00575         for(j=0; j<DS->TSD_MeltF[i].length; j++)
00576         {
00577             DS->TSD_MeltF[i].TS[j] = (realtype *)malloc(2*sizeof(realtype));
00578         }
00579 
00580         for(j=0; j<DS->TSD_MeltF[i].length; j++)
00581         {
00582             fscanf(forc_file, "%lf %lf", &DS->TSD_MeltF[i].TS[j][0], &DS->TSD_MeltF[i].TS[j][1]);
00583         }
00584         DS->TSD_MeltF[i].iCounter=0;
00585     }
00586 
00587     /*    Read All the Sources/Sinks Time Series    */
00588     for(i=0; i<DS->NumSource; i++)
00589     {
00590         fscanf(forc_file, "%s %d %d", DS->TSD_Source[i].name, &DS->TSD_Source[i].index, &DS->TSD_Source[i].length);
00591 
00592         DS->TSD_Source[i].TS = (realtype **)malloc((DS->TSD_Source[i].length)*sizeof(realtype));
00593 
00594         for(j=0; j<DS->TSD_Source[i].length; j++)
00595         {
00596             DS->TSD_Source[i].TS[j] = (realtype *)malloc(2*sizeof(realtype));
00597         }
00598 
00599         for(j=0; j<DS->TSD_Source[i].length; j++)
00600         {
00601             fscanf(forc_file, "%lf %lf", &DS->TSD_Source[i].TS[j][0], &DS->TSD_Source[i].TS[j][1]);
00602         }
00603 
00604         DS->TSD_Source[i].iCounter=0;
00605     }
00606 
00607     fclose(forc_file);
00608     printf("done.\n");
00609     /* Finish reading .forc File */
00610 
00611     /***************************************/
00612     /*========== open *.ibc file ==========*/
00613     /***************************************/
00614     printf("\n  6) reading %s.ibc  ... ", filename);
00615     fn[6] = (char *)malloc((strlen(filename)+4)*sizeof(char));
00616     strcpy(fn[6], filename);
00617     ibc_file =  fopen(strcat(fn[6], ".ibc"), "r");
00618 
00619     if(ibc_file == NULL)
00620     {
00621         printf("\n  Fatal Error: %s.ibc is in use or does not exist!\n", filename);
00622         exit(1);
00623     }
00624 
00625     /* start reading .ibc File */
00626     fscanf(ibc_file, "%d %d", &DS->Num1BC, &DS->Num2BC);
00627 
00628     if(DS->Num1BC+DS->Num2BC > 0)
00629     {
00630         DS->TSD_EleBC = (TSD *)malloc((DS->Num1BC+DS->Num2BC)*sizeof(TSD));
00631     }
00632 
00633     if(DS->Num1BC>0)
00634     {    /* For elements with Dirichilet Boundary Conditions */
00635         for(i=0; i<DS->Num1BC; i++)
00636         {
00637             fscanf(ibc_file, "%s %d %d", DS->TSD_EleBC[i].name, &DS->TSD_EleBC[i].index,
00638                                    &DS->TSD_EleBC[i].length);
00639 
00640             DS->TSD_EleBC[i].TS = (realtype **)malloc((DS->TSD_EleBC[i].length)*sizeof(realtype));
00641 
00642             for(j=0; j<DS->TSD_EleBC[i].length; j++)
00643             {
00644                 DS->TSD_EleBC[i].TS[j] = (realtype *)malloc(2*sizeof(realtype));
00645             }
00646 
00647             for(j=0; j<DS->TSD_EleBC[i].length; j++)
00648             {
00649                 fscanf(ibc_file, "%lf %lf", &DS->TSD_EleBC[i].TS[j][0],
00650                                    &DS->TSD_EleBC[i].TS[j][1]);
00651             }
00652             DS->TSD_EleBC[i].iCounter=0;
00653         }
00654     }
00655 
00656     if(DS->Num2BC>0)
00657     {    /* For elements with Nueman (non-natural) Boundary Conditions */
00658         /* This part of code has not be tested ! */
00659         for(i=DS->Num1BC; i<DS->Num1BC+DS->Num2BC; i++)
00660         {
00661             fscanf(ibc_file, "%s %d %d", DS->TSD_EleBC[i].name, &DS->TSD_EleBC[i].index,
00662                                    &DS->TSD_EleBC[i].length);
00663 
00664             DS->TSD_EleBC[i].TS = (realtype **)malloc((DS->TSD_EleBC[i].length)*sizeof(realtype));
00665             for(j=0; j<DS->TSD_EleBC[i].length; j++)
00666             {
00667                 DS->TSD_EleBC[i].TS[j] = (realtype *)malloc(2*sizeof(realtype));
00668             }
00669             for(j=0; j<DS->TSD_EleBC[i].length; j++)
00670             {
00671                 fscanf(forc_file, "%lf %lf", &DS->TSD_EleBC[i].TS[j][0], &DS->TSD_EleBC[i].TS[j][1]);
00672             }
00673             DS->TSD_EleBC[i].iCounter=0;
00674         }
00675     }
00676     fscanf(ibc_file, "%d", &DS->NumEleIC);
00677     /*DS->Ele_IC = (element_IC *)malloc(DS->NumEleIC*sizeof(element_IC));
00678     for(i=0; i<DS->NumEleIC; i++)
00679     {
00680         fscanf(ibc_file, "%d", &DS->Ele_IC[i].index);
00681         fscanf(ibc_file, "%lf", &DS->Ele_IC[i].interception);
00682         fscanf(ibc_file, "%lf", &DS->Ele_IC[i].snow);
00683         fscanf(ibc_file, "%lf", &DS->Ele_IC[i].surf);
00684         fscanf(ibc_file, "%lf", &DS->Ele_IC[i].unsat);
00685         fscanf(ibc_file, "%lf", &DS->Ele_IC[i].sat);
00686     }*/
00687 
00688     fclose(ibc_file);
00689     printf("done.\n");
00690     /* Finish reading .ibc File */
00691 
00692     /****************************************/
00693     /*========== open *.para file ==========*/
00694     /****************************************/
00695     printf("\n  7) reading %s.para ... ", filename);
00696     fn[7] = (char *)malloc((strlen(filename)+5)*sizeof(char));
00697     strcpy(fn[7], filename);
00698     para_file = fopen(strcat(fn[7], ".para"), "r");
00699 
00700     if(para_file == NULL)
00701     {
00702         printf("\n  Fatal Error: %s.para is in use or does not exist!\n", filename);
00703         exit(1);
00704     }
00705 
00706     /* start reading .para File */
00707     fscanf(para_file, "%d %d", &CS->Verbose, &CS->Debug);
00708     fscanf(para_file, "%d", &CS->int_type);
00709     fscanf(para_file, "%d %d %d %d", &CS->res_out, &CS->flux_out, &CS->q_out, &CS->etis_out);
00710     fscanf(para_file, "%d %d %d", &DS->UnsatMode, &DS->SurfMode, &DS->RivMode);
00711     fscanf(para_file, "%d", &CS->Solver);
00712     if(CS->Solver == 2)
00713     {
00714         fscanf(para_file, "%d %d %lf", &CS->GSType, &CS->MaxK, &CS->delt);
00715     }
00716     fscanf(para_file, "%lf %lf", &CS->abstol, &CS->reltol);
00717     fscanf(para_file, "%lf %lf %lf", &CS->InitStep, &CS->MaxStep, &CS->ETStep);
00718     fscanf(para_file, "%lf %lf %d", &CS->StartTime, &CS->EndTime, &CS->outtype);
00719     if(CS->outtype == 0)
00720     {
00721         fscanf(para_file, "%lf %lf", &CS->a, &CS->b);
00722     }
00723 
00724     if(CS->a != 1.0)
00725     {
00726         NumTout = (int)(log(1 - (CS->EndTime - CS->StartTime)*(1 -  CS->a)/CS->b)/log(CS->a));
00727     }
00728     else
00729     {
00730         if((CS->EndTime - CS->StartTime)/CS->b - ((int) (CS->EndTime - CS->StartTime)/CS->b) > 0)
00731         {
00732             NumTout = (int) ((CS->EndTime - CS->StartTime)/CS->b);
00733         }
00734         else
00735         {
00736             NumTout = (int) ((CS->EndTime - CS->StartTime)/CS->b - 1);
00737         }
00738     }
00739 
00740     CS->NumSteps = NumTout + 1;
00741 
00742     CS->Tout = (realtype *)malloc((CS->NumSteps + 1)*sizeof(realtype));
00743 
00744     for(i=0; i<CS->NumSteps+1; i++)
00745     {
00746         if(i == 0)
00747         {
00748             CS->Tout[i] = CS->StartTime;
00749         }
00750         else
00751         {
00752             CS->Tout[i] = CS->Tout[i-1] + pow(CS->a, i)*CS->b;
00753         }
00754     }
00755 
00756     if(CS->Tout[CS->NumSteps] < CS->EndTime)
00757     {
00758         CS->Tout[CS->NumSteps] = CS->EndTime;
00759     }
00760 
00761     fclose(para_file);
00762     printf("done.\n");
00763     /* Finish reading .para File */
00764 
00765 }
00766
Generated on Thu Jul 12 14:34:19 2007 for PIHM by  doxygen 1.5.2