print.c

Go to the documentation of this file.
00001 /*******************************************************************************
00002  * File        : print.c                                                       *
00003  * Function    : prints the output parameters in either txt or netcdf format   *
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  *                                                                             *
00011  * This code is free for users with research purpose only, if appropriate      *
00012  * citation is refered. However, there is no warranty in any format for this   *
00013  * product.                                                                    *
00014  *                                                                             *
00015  * For questions or comments, please contact the authors of the reference.     *
00016  * One who want to use it for other consideration may also contact Dr.Duffy    *
00017  * at cxd11@psu.edu.                                                           *
00018  *******************************************************************************/
00019 
00021 
00022 /*    Header File        */
00023 #include <stdio.h>
00024 #include <stdlib.h>
00025 #include <math.h>
00026 #include <string.h>
00027 
00028 /*    NETCDF Header File    */
00029 #include <netcdf.h>
00030 
00031 /*    SUNDIALS Header File    */
00032 #include "sundials_types.h"
00033 #include "cvode.h"
00034 #include "cvode_dense.h"
00035 #include "nvector_serial.h"
00036 
00037 /*    PIHM Header Files    */
00038 #include "pihm.h"
00039 #include "calib.h"
00040 #include "print.h"
00041 
00042 
00043 #define NDIMS    2    
00045 #define ERR(e) {printf("Error: %s\n", nc_strerror(e)); return;}
00046 
00047 
00048 /*    File Pointer for TXT files    */
00049 FILE *isStatePtr, *satStatePtr, *usatStatePtr, *surfStatePtr;
00050 FILE *et0Ptr, *et1Ptr, *et2Ptr, *netPrecipPtr, *infilPtr, *rechargePtr;
00051 FILE *rivHeadPtr;
00052 FILE *rivFlowPtr, *rivBasePtr, *rivSurfPtr;
00053 
00054 /*    Strings to hold File Names    */
00055 char *isStateFile, *satStateFile, *usatStateFile, *surfStateFile;
00056 char *et0File, *et1File, *et2File, *netPrecipFile, *infilFile, *rechargeFile;
00057 char *rivHeadFile;
00058 char *rivFlowFile, *rivBaseFile, *rivSurfFile;
00059 
00060 FILE *initPtr;     
00061 char *initFile;    
00063 /*    variables to hold mean values */
00064 static double *tempIS,  *tempSatState, *tempUsatState, *tempSurfState;
00065 static double *tempET0, *tempET1, *tempET2, *tempNetPpt, *tempInfil, *tempRecharge;
00066 static double *tempFlow, *tempBase, *tempSurf, *tempHead;
00067 
00068 
00069 int NUMELE;        
00070 int NUMRIV;        
00072 /* Error handling. */
00073 int retval;        
00076 /* ncIDs for CDF (.nc) Files    */
00077 int isStateID, satStateID, usatStateID, surfStateID;
00078 int et0ID, et1ID, et2ID, netPrecipID, infilID, rechargeID;
00079 int rivHeadID;
00080 int rivFlowID, rivBaseID, rivSurfID;
00081 
00082 int ele_dimid;                           
00083 int rec_dimid;                           
00084 int dimids[NDIMS];                       
00085 int startEle[NDIMS];                     
00086 int countEle[NDIMS];                     
00087 int startRiv[NDIMS];                     
00088 int countRiv[NDIMS];                     
00090 /* Variable IDs for CDF (.nc) Files    */
00091 int isState_varid, satState_varid, usatState_varid, surfState_varid;
00092 int et0_varid, et1_varid, et2_varid, netPrecip_varid, infil_varid, recharge_varid;
00093 int rivHead_varid;
00094 int rivFlow_varid, rivBase_varid, rivSurf_varid;
00095 
00096 
00097 /********************************************************************
00098     This function calls different fuction depending on the Output File Mode
00099     and simulated variables user wants to print as declared in print.h file
00100 *********************************************************************/
00101 void FPrint(Model_Data mData, N_Vector CV_Y, realtype t)
00103 
00107 {
00108     /***************    TXT FILE MODE : START    ****************/
00109     if(FPRINT_MODE==TXT){
00110         if(ISState==YEA){
00111             printIS(mData, isStatePtr, t);
00112         }
00113         if(SatState==YEA){
00114             printSatState(mData, CV_Y, satStatePtr, t);
00115         }
00116         if(UsatState==YEA){
00117             printUsatState(mData, CV_Y, usatStatePtr, t);
00118         }
00119         if(SurfState==YEA){
00120             printSurfState(mData, CV_Y, surfStatePtr, t);
00121         }
00122 
00123         if(ET0==YEA){
00124             printET0(mData, et0Ptr, t);
00125         }
00126         if(ET1==YEA){
00127             printET1(mData, et1Ptr, t);
00128         }
00129         if(ET2==YEA){
00130             printET2(mData, et2Ptr, t);
00131         }
00132         if(NetPpt==YEA){
00133             printNetPpt(mData, netPrecipPtr, t);
00134         }
00135         if(Infil==YEA){
00136             printInfil(mData, infilPtr, t);
00137         }
00138         if(RECHARGE==YEA){
00139             printRecharge(mData, rechargePtr, t);
00140         }
00141 
00142         if(RivFlow==YEA){
00143             printRiverFlow(mData, CV_Y, rivFlowPtr, t);
00144         }
00145         if(RivBase==YEA){
00146             printRiverBase(mData, rivBasePtr, t);
00147         }
00148         if(RivSurf==YEA){
00149             printRiverSurf(mData, rivSurfPtr, t);
00150         }
00151         if(RivHead==YEA){
00152             printRiverHead(mData, CV_Y, rivHeadPtr, t);
00153         }
00154     }
00155     /***************    TXT FILE MODE : END    ****************/
00156 
00157 
00158     /***************    CDF FILE MODE : START    ****************/
00159     if(FPRINT_MODE==CDF){
00160         if(ISState==YEA){
00161             printIScdf(mData, isStateID, isState_varid, t);
00162         }
00163         if(SatState==YEA){
00164             printSatStatecdf(mData, CV_Y, satStateID, satState_varid, t);
00165         }
00166         if(UsatState==YEA){
00167             printUsatStatecdf(mData, CV_Y, usatStateID, usatState_varid, t);
00168         }
00169         if(SurfState==YEA){
00170             printSurfStatecdf(mData, CV_Y, surfStateID, surfState_varid, t);
00171         }
00172 
00173         if(ET0==YEA){
00174             printET0cdf(mData, et0ID, et0_varid, t);
00175         }
00176         if(ET1==YEA){
00177             printET1cdf(mData, et1ID, et1_varid, t);
00178         }
00179         if(ET2==YEA){
00180             printET2cdf(mData, et2ID, et2_varid, t);
00181         }
00182         if(NetPpt==YEA){
00183             printNetPptcdf(mData, netPrecipID, netPrecip_varid, t);
00184         }
00185         if(Infil==YEA){
00186             printInfilcdf(mData, infilID, infil_varid, t);
00187         }
00188         if(RECHARGE==YEA){
00189             printRechargecdf(mData, rechargeID, recharge_varid, t);
00190         }
00191 
00192 
00193         if(RivFlow==YEA){
00194             printRiverFlowcdf(mData, CV_Y, rivFlowID, rivFlow_varid, t);
00195         }
00196         if(RivBase==YEA){
00197             printRiverBasecdf(mData, rivBaseID, rivBase_varid, t);
00198         }
00199         if(RivSurf==YEA){
00200             printRiverSurfcdf(mData, rivSurfID, rivSurf_varid, t);
00201         }
00202         if(RivHead==YEA){
00203             printRiverHeadcdf(mData, CV_Y, rivHeadID, rivHead_varid, t);
00204         }
00205     }
00206     /***************    CDF FILE MODE : END    ****************/
00207 
00208 }
00209 
00210 /****************************************************************
00211 ROUTINE TO PRINT NEW INIT FILE AT THE END OF SIMULATION
00212 $$ IT ENABLES USER TO START SIMULATION FROM THERE ONWARDS
00213 *****************************************************************/
00214 
00215 void FPrintInitFile(Model_Data mData, Control_Data cData, N_Vector CV_Y, int i)
00217 
00222 {
00223     int j;
00224     char tmpFileName[100];
00225     setFileName(tmpFileName);
00226     initFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00227     strcpy(initFile, tmpFileName);
00228     strcat(initFile, ".init.end");    //TODO: replace 'end' with 'time value'
00229     initPtr=fopen(initFile, "w");
00230 
00231     fprintf(initPtr,"%lf\n",cData.Tout[i]);
00232     for(j=0; j<mData->NumEle; j++)
00233     {
00234         fprintf(initPtr,"%lf\t%lf\t%lf\t%lf\t%lf\n",mData->EleIS[j],mData->EleSnow[j],NV_Ith_S(CV_Y,j), NV_Ith_S(CV_Y,j+mData->NumEle), NV_Ith_S(CV_Y,j+2*mData->NumEle));
00235     }
00236     for(j=0; j<mData->NumRiv; j++)
00237     {
00238         fprintf(initPtr,"%lf\n",NV_Ith_S(CV_Y, j+3*mData->NumEle));
00239     }
00240     fflush(initPtr);
00241 }
00242 
00243 /****************************************************************
00244 Open the files depending on the flag values (YEA/NAY) and
00245 mode of output (TXT/NETCDF)
00246 *****************************************************************/
00247 void FPrintInit(Model_Data mData)
00249 
00251 {
00252     int i;
00253     char tmpFileName[100];
00254     setFileName(tmpFileName);
00255 
00256     NUMELE=mData->NumEle;
00257     NUMRIV=mData->NumRiv;
00258 
00259     countEle[0]=1;
00260     countEle[1]=NUMELE;
00261     startEle[1]=0;
00262     countRiv[0]=1;
00263     countRiv[1]=NUMRIV;
00264     startEle[1]=0;
00265     /***************    TXT FILE MODE : START    ****************/
00266     if(FPRINT_MODE==TXT){
00267         /* STATES */
00268         if(ISState==YEA){
00269             isStateFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00270             strcpy(isStateFile, tmpFileName);
00271             strcat(isStateFile, ".is.txt");
00272             isStatePtr=fopen(isStateFile, "w");
00273 
00274             tempIS=(double *)malloc(mData->NumEle * sizeof(double));
00275             for(i=0; i<mData->NumEle; i++)
00276                 tempIS[i]=0.0;
00277         }
00278         if(SatState==YEA){
00279             satStateFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00280             strcpy(satStateFile, tmpFileName);
00281             strcat(satStateFile, ".sat.txt");
00282             satStatePtr=fopen(satStateFile, "w");
00283 
00284             tempSatState=(double *)malloc(mData->NumEle * sizeof(double));
00285             for(i=0; i<mData->NumEle; i++)
00286                 tempSatState[i]=0.0;
00287         }
00288         if(UsatState==YEA){
00289             usatStateFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00290             strcpy(usatStateFile, tmpFileName);
00291             strcat(usatStateFile, ".usat.txt");
00292             usatStatePtr=fopen(usatStateFile, "w");
00293 
00294             tempUsatState=(double *)malloc(mData->NumEle * sizeof(double));
00295             for(i=0; i<mData->NumEle; i++)
00296                 tempUsatState[i]=0.0;
00297         }
00298         if(SurfState==YEA){
00299             surfStateFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00300             strcpy(surfStateFile, tmpFileName);
00301             strcat(surfStateFile, ".surf.txt");
00302             surfStatePtr=fopen(surfStateFile, "w");
00303 
00304             tempSurfState=(double *)malloc(mData->NumEle * sizeof(double));
00305             for(i=0; i<mData->NumEle; i++)
00306                 tempSurfState[i]=0.0;
00307         }
00308 
00309         /* FLUXES */
00310         if(ET0==YEA){
00311             et0File = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00312             strcpy(et0File, tmpFileName);
00313             strcat(et0File, ".et0.txt");
00314             et0Ptr=fopen(et0File, "w");
00315 
00316             tempET0=(double *)malloc(mData->NumEle * sizeof(double));
00317             for(i=0; i<mData->NumEle; i++)
00318                 tempET0[i]=0.0;
00319         }
00320         if(ET1==YEA){
00321             et1File = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00322             strcpy(et1File, tmpFileName);
00323             strcat(et1File, ".et1.txt");
00324             et1Ptr=fopen(et1File, "w");
00325 
00326             tempET1=(double *)malloc(mData->NumEle * sizeof(double));
00327             for(i=0; i<mData->NumEle; i++)
00328                 tempET1[i]=0.0;
00329         }
00330         if(ET2==YEA){
00331             et2File = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00332             strcpy(et2File, tmpFileName);
00333             strcat(et2File, ".et2.txt");
00334             et2Ptr=fopen(et2File, "w");
00335 
00336             tempET2=(double *)malloc(mData->NumEle * sizeof(double));
00337             for(i=0; i<mData->NumEle; i++)
00338                 tempET2[i]=0.0;
00339         }
00340         if(NetPpt==YEA){
00341             netPrecipFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00342             strcpy(netPrecipFile, tmpFileName);
00343             strcat(netPrecipFile, ".netPrecip.txt");
00344             netPrecipPtr=fopen(netPrecipFile, "w");
00345 
00346             tempNetPpt=(double *)malloc(mData->NumEle * sizeof(double));
00347             for(i=0; i<mData->NumEle; i++)
00348                 tempNetPpt[i]=0.0;
00349         }
00350         if(Infil==YEA){
00351             infilFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00352             strcpy(infilFile, tmpFileName);
00353             strcat(infilFile, ".infil.txt");
00354             infilPtr=fopen(infilFile, "w");
00355 
00356             tempInfil=(double *)malloc(mData->NumEle * sizeof(double));
00357             for(i=0; i<mData->NumEle; i++)
00358                 tempInfil[i]=0.0;
00359         }
00360         if(RECHARGE==YEA){
00361             rechargeFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00362             strcpy(rechargeFile, tmpFileName);
00363             strcat(rechargeFile, ".recharge.txt");
00364             rechargePtr=fopen(rechargeFile, "w");
00365 
00366             tempRecharge=(double *)malloc(mData->NumEle * sizeof(double));
00367             for(i=0; i<mData->NumEle; i++)
00368                 tempRecharge[i]=0.0;
00369         }
00370 
00371         /* River States */
00372         if(RivHead==YEA){
00373             rivHeadFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00374             strcpy(rivHeadFile, tmpFileName);
00375             strcat(rivHeadFile, ".rivHead.txt");
00376             rivHeadPtr=fopen(rivHeadFile, "w");
00377 
00378             tempHead=(double *)malloc(mData->NumRiv * sizeof(double));
00379             for(i=0; i<mData->NumRiv; i++)
00380                 tempHead[i]=0.0;
00381         }
00382 
00383         /* River Fluxes */
00384         if(RivFlow==YEA){
00385             rivFlowFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00386             strcpy(rivFlowFile, tmpFileName);
00387             strcat(rivFlowFile, ".rivFlow.txt");
00388             rivFlowPtr=fopen(rivFlowFile, "w");
00389 
00390             tempFlow=(double *)malloc(mData->NumRiv * sizeof(double));
00391             for(i=0; i<mData->NumRiv; i++)
00392                 tempFlow[i]=0.0;
00393         }
00394         if(RivBase==YEA){
00395             rivBaseFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00396             strcpy(rivBaseFile, tmpFileName);
00397             strcat(rivBaseFile, ".rivBase.txt");
00398             rivBasePtr=fopen(rivBaseFile, "w");
00399 
00400             tempBase=(double *)malloc(mData->NumRiv * sizeof(double));
00401             for(i=0; i<mData->NumRiv; i++)
00402                 tempBase[i]=0.0;
00403         }
00404         if(RivSurf==YEA){
00405             rivSurfFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00406             strcpy(rivSurfFile, tmpFileName);
00407             strcat(rivSurfFile, ".rivSurf.txt");
00408             rivSurfPtr=fopen(rivSurfFile, "w");
00409 
00410             tempSurf=(double *)malloc(mData->NumRiv * sizeof(double));
00411             for(i=0; i<mData->NumRiv; i++)
00412                 tempSurf[i]=0.0;
00413         }
00414     }
00415     /***************    TXT FILE MODE : END    ****************/
00416 
00417 
00418     /***************    CDF FILE MODE : START    ****************/
00419     if(FPRINT_MODE==CDF){
00420         /* STATES */
00421         if(ISState==YEA){
00422             isStateFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00423             strcpy(isStateFile, tmpFileName);
00424             strcat(isStateFile, ".is.nc");
00425             if ((retval = nc_create(isStateFile, NC_CLOBBER, &isStateID)))
00426                   ERR(retval);
00427                if ((retval = nc_def_dim(isStateID, "Elements" , NUMELE, &ele_dimid)))
00428                   ERR(retval);
00429                if ((retval = nc_def_dim(isStateID, "Time", NC_UNLIMITED, &rec_dimid)))
00430                   ERR(retval);
00431               dimids[0]=rec_dimid;
00432               dimids[1]=ele_dimid;
00433             if((retval = nc_def_var(isStateID, "Interception_Storage_State", NC_DOUBLE, NDIMS, dimids, &isState_varid)))
00434                 ERR(retval);
00435             if ((retval = nc_enddef(isStateID)))
00436                 ERR(retval);
00437 
00438             tempIS=(double *)malloc(mData->NumEle * sizeof(double));
00439             for(i=0; i<mData->NumEle; i++)
00440                 tempIS[i]=0.0;
00441         }
00442         if(SatState==YEA){
00443             satStateFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00444             strcpy(satStateFile, tmpFileName);
00445             strcat(satStateFile, ".sat.nc");
00446             if ((retval = nc_create(satStateFile, NC_CLOBBER, &satStateID)))
00447                   ERR(retval);
00448                if ((retval = nc_def_dim(satStateID, "Elements" , NUMELE, &ele_dimid)))
00449                   ERR(retval);
00450                if ((retval = nc_def_dim(satStateID, "time", NC_UNLIMITED, &rec_dimid)))
00451                   ERR(retval);
00452               dimids[0]=rec_dimid;
00453               dimids[1]=ele_dimid;
00454             if((retval = nc_def_var(satStateID, "Saturated_Zone_State", NC_DOUBLE, NDIMS, dimids, &satState_varid)))
00455                 ERR(retval);
00456             if ((retval = nc_enddef(satStateID)))
00457                 ERR(retval);
00458 
00459             tempSatState=(double *)malloc(mData->NumEle * sizeof(double));
00460             for(i=0; i<mData->NumEle; i++)
00461                 tempSatState[i]=0.0;
00462         }
00463         if(UsatState==YEA){
00464             usatStateFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00465             strcpy(usatStateFile, tmpFileName);
00466             strcat(usatStateFile, ".usat.nc");
00467             if ((retval = nc_create(usatStateFile, NC_CLOBBER, &usatStateID)))
00468                   ERR(retval);
00469                if ((retval = nc_def_dim(usatStateID, "Elements" , NUMELE, &ele_dimid)))
00470                   ERR(retval);
00471                if ((retval = nc_def_dim(usatStateID, "time", NC_UNLIMITED, &rec_dimid)))
00472                   ERR(retval);
00473               dimids[0]=rec_dimid;
00474               dimids[1]=ele_dimid;
00475             if((retval = nc_def_var(usatStateID, "Unsaturated_Zone_State", NC_DOUBLE, NDIMS, dimids, &usatState_varid)))
00476                 ERR(retval);
00477             if ((retval = nc_enddef(usatStateID)))
00478                 ERR(retval);
00479 
00480             tempUsatState=(double *)malloc(mData->NumEle * sizeof(double));
00481             for(i=0; i<mData->NumEle; i++)
00482                 tempUsatState[i]=0.0;
00483         }
00484         if(SurfState==YEA){
00485             surfStateFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00486             strcpy(surfStateFile, tmpFileName);
00487             strcat(surfStateFile, ".surf.nc");
00488             if ((retval = nc_create(surfStateFile, NC_CLOBBER, &surfStateID)))
00489                   ERR(retval);
00490                if ((retval = nc_def_dim(surfStateID, "Elements" , NUMELE, &ele_dimid)))
00491                   ERR(retval);
00492                if ((retval = nc_def_dim(surfStateID, "time", NC_UNLIMITED, &rec_dimid)))
00493                   ERR(retval);
00494               dimids[0]=rec_dimid;
00495               dimids[1]=ele_dimid;
00496             if((retval = nc_def_var(surfStateID, "Surface_Flow_State", NC_DOUBLE, NDIMS, dimids, &surfState_varid)))
00497                 ERR(retval);
00498             if ((retval = nc_enddef(surfStateID)))
00499                 ERR(retval);
00500 
00501             tempSurfState=(double *)malloc(mData->NumEle * sizeof(double));
00502             for(i=0; i<mData->NumEle; i++)
00503                 tempSurfState[i]=0.0;
00504         }
00505 
00506         /* FLUXES */
00507         if(ET0==YEA){
00508             et0File = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00509             strcpy(et0File, tmpFileName);
00510             strcat(et0File, ".et0.nc");
00511             if ((retval = nc_create(et0File, NC_CLOBBER, &et0ID)))
00512                   ERR(retval);
00513                if ((retval = nc_def_dim(et0ID, "Elements" , NUMELE, &ele_dimid)))
00514                   ERR(retval);
00515                if ((retval = nc_def_dim(et0ID, "time", NC_UNLIMITED, &rec_dimid)))
00516                   ERR(retval);
00517               dimids[0]=rec_dimid;
00518               dimids[1]=ele_dimid;
00519             if((retval = nc_def_var(et0ID, "ET0", NC_DOUBLE, NDIMS, dimids, &et0_varid)))
00520                 ERR(retval);
00521             if ((retval = nc_enddef(et0ID)))
00522                 ERR(retval);
00523 
00524             tempET0=(double *)malloc(mData->NumEle * sizeof(double));
00525             for(i=0; i<mData->NumEle; i++)
00526                 tempET0[i]=0.0;
00527         }
00528         if(ET1==YEA){
00529             et1File = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00530             strcpy(et1File, tmpFileName);
00531             strcat(et1File, ".et1.nc");
00532             if ((retval = nc_create(et1File, NC_CLOBBER, &et1ID)))
00533                   ERR(retval);
00534                if ((retval = nc_def_dim(et1ID, "Elements" , NUMELE, &ele_dimid)))
00535                   ERR(retval);
00536                if ((retval = nc_def_dim(et1ID, "time", NC_UNLIMITED, &rec_dimid)))
00537                   ERR(retval);
00538               dimids[0]=rec_dimid;
00539               dimids[1]=ele_dimid;
00540             if((retval = nc_def_var(et1ID, "ET1", NC_DOUBLE, NDIMS, dimids, &et1_varid)))
00541                 ERR(retval);
00542             if ((retval = nc_enddef(et1ID)))
00543                 ERR(retval);
00544 
00545             tempET1=(double *)malloc(mData->NumEle * sizeof(double));
00546             for(i=0; i<mData->NumEle; i++)
00547                 tempET1[i]=0.0;
00548         }
00549         if(ET2==YEA){
00550             et2File = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00551             strcpy(et2File, tmpFileName);
00552             strcat(et2File, ".et2.nc");
00553             if ((retval = nc_create(et2File, NC_CLOBBER, &et2ID)))
00554                   ERR(retval);
00555                if ((retval = nc_def_dim(et2ID, "Elements" , NUMELE, &ele_dimid)))
00556                   ERR(retval);
00557                if ((retval = nc_def_dim(et2ID, "time", NC_UNLIMITED, &rec_dimid)))
00558                   ERR(retval);
00559               dimids[0]=rec_dimid;
00560               dimids[1]=ele_dimid;
00561             if((retval = nc_def_var(et2ID, "ET2", NC_DOUBLE, NDIMS, dimids, &et2_varid)))
00562                 ERR(retval);
00563             if ((retval = nc_enddef(et2ID)))
00564                 ERR(retval);
00565 
00566             tempET2=(double *)malloc(mData->NumEle * sizeof(double));
00567             for(i=0; i<mData->NumEle; i++)
00568                 tempET2[i]=0.0;
00569         }
00570         if(NetPpt==YEA){
00571             netPrecipFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00572             strcpy(netPrecipFile, tmpFileName);
00573             strcat(netPrecipFile, ".netPrecip.nc");
00574             if ((retval = nc_create(netPrecipFile, NC_CLOBBER, &netPrecipID)))
00575                   ERR(retval);
00576                if ((retval = nc_def_dim(netPrecipID, "Elements" , NUMELE, &ele_dimid)))
00577                   ERR(retval);
00578                if ((retval = nc_def_dim(netPrecipID, "time", NC_UNLIMITED, &rec_dimid)))
00579                   ERR(retval);
00580               dimids[0]=rec_dimid;
00581               dimids[1]=ele_dimid;
00582             if((retval = nc_def_var(netPrecipID, "Net_Precipitation", NC_DOUBLE, NDIMS, dimids, &netPrecip_varid)))
00583                 ERR(retval);
00584             if ((retval = nc_enddef(netPrecipID)))
00585                 ERR(retval);
00586 
00587             tempNetPpt=(double *)malloc(mData->NumEle * sizeof(double));
00588             for(i=0; i<mData->NumEle; i++)
00589                 tempNetPpt[i]=0.0;
00590         }
00591         if(Infil==YEA){
00592             infilFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00593             strcpy(infilFile, tmpFileName);
00594             strcat(infilFile, ".infil.nc");
00595             if ((retval = nc_create(infilFile, NC_CLOBBER, &infilID)))
00596                   ERR(retval);
00597                if ((retval = nc_def_dim(infilID, "Elements" , NUMELE, &ele_dimid)))
00598                   ERR(retval);
00599                if ((retval = nc_def_dim(infilID, "time", NC_UNLIMITED, &rec_dimid)))
00600                   ERR(retval);
00601               dimids[0]=rec_dimid;
00602               dimids[1]=ele_dimid;
00603             if((retval = nc_def_var(infilID, "Infiltration", NC_DOUBLE, NDIMS, dimids, &infil_varid)))
00604                 ERR(retval);
00605             if ((retval = nc_enddef(infilID)))
00606                 ERR(retval);
00607 
00608             tempInfil=(double *)malloc(mData->NumEle * sizeof(double));
00609             for(i=0; i<mData->NumEle; i++)
00610                 tempInfil[i]=0.0;
00611         }
00612         if(RECHARGE==YEA){
00613             rechargeFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00614             strcpy(rechargeFile, tmpFileName);
00615             strcat(rechargeFile, ".recharge.nc");
00616             if ((retval = nc_create(rechargeFile, NC_CLOBBER, &rechargeID)))
00617                   ERR(retval);
00618                if ((retval = nc_def_dim(rechargeID, "Elements" , NUMELE, &ele_dimid)))
00619                   ERR(retval);
00620                if ((retval = nc_def_dim(rechargeID, "time", NC_UNLIMITED, &rec_dimid)))
00621                   ERR(retval);
00622               dimids[0]=rec_dimid;
00623               dimids[1]=ele_dimid;
00624             if((retval = nc_def_var(rechargeID, "Recharge2GW", NC_DOUBLE, NDIMS, dimids, &recharge_varid)))
00625                 ERR(retval);
00626             if ((retval = nc_enddef(rechargeID)))
00627                 ERR(retval);
00628 
00629             tempRecharge=(double *)malloc(mData->NumEle * sizeof(double));
00630             for(i=0; i<mData->NumEle; i++)
00631                 tempRecharge[i]=0.0;
00632         }
00633 
00634         /* River States */
00635         if(RivHead==YEA){
00636             rivHeadFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00637             strcpy(rivHeadFile, tmpFileName);
00638             strcat(rivHeadFile, ".rivHead.nc");
00639             if ((retval = nc_create(rivHeadFile, NC_CLOBBER, &rivHeadID)))
00640                   ERR(retval);
00641                if ((retval = nc_def_dim(rivHeadID, "RiverSegments" , NUMRIV, &ele_dimid)))
00642                   ERR(retval);
00643                if ((retval = nc_def_dim(rivHeadID, "time", NC_UNLIMITED, &rec_dimid)))
00644                   ERR(retval);
00645               dimids[0]=rec_dimid;
00646               dimids[1]=ele_dimid;
00647             if((retval = nc_def_var(rivHeadID, "RivState", NC_DOUBLE, NDIMS, dimids, &rivHead_varid)))
00648                 ERR(retval);
00649             if ((retval = nc_enddef(rivHeadID)))
00650                 ERR(retval);
00651 
00652             tempHead=(double *)malloc(mData->NumRiv * sizeof(double));
00653             for(i=0; i<mData->NumRiv; i++)
00654                 tempHead[i]=0.0;
00655         }
00656 
00657         /* River Fluxes */
00658         if(RivFlow==YEA){
00659             rivFlowFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00660             strcpy(rivFlowFile, tmpFileName);
00661             strcat(rivFlowFile, ".rivFlow.nc");
00662             if ((retval = nc_create(rivFlowFile, NC_CLOBBER, &rivFlowID)))
00663                   ERR(retval);
00664                if ((retval = nc_def_dim(rivFlowID, "RiverSegments" , NUMRIV, &ele_dimid)))
00665                   ERR(retval);
00666                if ((retval = nc_def_dim(rivFlowID, "time", NC_UNLIMITED, &rec_dimid)))
00667                   ERR(retval);
00668               dimids[0]=rec_dimid;
00669               dimids[1]=ele_dimid;
00670             if((retval = nc_def_var(rivFlowID, "RivFlow", NC_DOUBLE, NDIMS, dimids, &rivFlow_varid)))
00671                 ERR(retval);
00672             if ((retval = nc_enddef(rivFlowID)))
00673                 ERR(retval);
00674 
00675             tempFlow=(double *)malloc(mData->NumRiv * sizeof(double));
00676             for(i=0; i<mData->NumRiv; i++)
00677                 tempFlow[i]=0.0;
00678         }
00679         if(RivBase==YEA){
00680             rivBaseFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00681             strcpy(rivBaseFile, tmpFileName);
00682             strcat(rivBaseFile, ".rivBase.nc");
00683             if ((retval = nc_create(rivBaseFile, NC_CLOBBER, &rivBaseID)))
00684                   ERR(retval);
00685                if ((retval = nc_def_dim(rivBaseID, "RiverSegments" , NUMRIV, &ele_dimid)))
00686                   ERR(retval);
00687                if ((retval = nc_def_dim(rivBaseID, "time", NC_UNLIMITED, &rec_dimid)))
00688                   ERR(retval);
00689               dimids[0]=rec_dimid;
00690               dimids[1]=ele_dimid;
00691             if((retval = nc_def_var(rivBaseID, "Base2Riv", NC_DOUBLE, NDIMS, dimids, &rivBase_varid)))
00692                 ERR(retval);
00693             if ((retval = nc_enddef(rivBaseID)))
00694                 ERR(retval);
00695 
00696             tempBase=(double *)malloc(mData->NumRiv * sizeof(double));
00697             for(i=0; i<mData->NumRiv; i++)
00698                 tempBase[i]=0.0;
00699         }
00700         if(RivSurf==YEA){
00701             rivSurfFile = (char *)malloc(sizeof(char)*(20+strlen(tmpFileName)));
00702             strcpy(rivSurfFile, tmpFileName);
00703             strcat(rivSurfFile, ".rivSurf.nc");
00704             if ((retval = nc_create(rivSurfFile, NC_CLOBBER, &rivSurfID)))
00705                   ERR(retval);
00706                if ((retval = nc_def_dim(rivSurfID, "RiverSegments" , NUMRIV, &ele_dimid)))
00707                   ERR(retval);
00708                if ((retval = nc_def_dim(rivSurfID, "time", NC_UNLIMITED, &rec_dimid)))
00709                   ERR(retval);
00710               dimids[0]=rec_dimid;
00711               dimids[1]=ele_dimid;
00712             if((retval = nc_def_var(rivSurfID, "Over2Riv", NC_DOUBLE, NDIMS, dimids, &rivSurf_varid)))
00713                 ERR(retval);
00714             if ((retval = nc_enddef(rivSurfID)))
00715                 ERR(retval);
00716 
00717             tempSurf=(double *)malloc(mData->NumRiv * sizeof(double));
00718             for(i=0; i<mData->NumRiv; i++)
00719                 tempSurf[i]=0.0;
00720         }
00721 
00722     }
00723     /***************    CDF FILE MODE : END    ****************/
00724 }
00725 
00726 /*********************************************
00727     Close All the files those were opened
00728     depending on the File Mode
00729 *********************************************/
00730 void FPrintCloseAll(void)
00732 {
00733     /* if File Mode is TXT    */
00734     if(FPRINT_MODE==TXT){
00735         if(ISState==YEA)
00736             fclose(isStatePtr);
00737         if(SatState==YEA)
00738             fclose(satStatePtr);
00739         if(UsatState==YEA)
00740             fclose(usatStatePtr);
00741         if(SurfState==YEA)
00742             fclose(surfStatePtr);
00743         if(ET0==YEA)
00744             fclose(et0Ptr);
00745         if(ET1==YEA)
00746             fclose(et1Ptr);
00747         if(ET2==YEA)
00748             fclose(et2Ptr);
00749         if(NetPpt==YEA)
00750             fclose(netPrecipPtr);
00751         if(Infil==YEA)
00752             fclose(infilPtr);
00753         if(RECHARGE==YEA)
00754             fclose(rechargePtr);
00755         if(RivHead==YEA)
00756             fclose(rivHeadPtr);
00757         if(RivFlow==YEA)
00758             fclose(rivFlowPtr);
00759         if(RivBase==YEA)
00760             fclose(rivBasePtr);
00761         if(RivSurf==YEA)
00762             fclose(rivSurfPtr);
00763     }
00764 
00765     /* if File Mode is CDF    */
00766     if(FPRINT_MODE==CDF){
00767         if(ISState==YEA)
00768             ncclose(isStateID);
00769         if(SatState==YEA)
00770             ncclose(satStateID);
00771         if(UsatState==YEA)
00772             ncclose(usatStateID);
00773         if(SurfState==YEA)
00774             ncclose(surfStateID);
00775         if(ET0==YEA)
00776             ncclose(et0ID);
00777         if(ET1==YEA)
00778             ncclose(et1ID);
00779         if(ET2==YEA)
00780             ncclose(et2ID);
00781         if(NetPpt==YEA)
00782             ncclose(netPrecipID);
00783         if(Infil==YEA)
00784             ncclose(infilID);
00785         if(RECHARGE==YEA)
00786             ncclose(rechargeID);
00787         if(RivHead==YEA)
00788             ncclose(rivHeadID);
00789         if(RivFlow==YEA)
00790             ncclose(rivFlowID);
00791         if(RivBase==YEA)
00792             ncclose(rivBaseID);
00793         if(RivSurf==YEA)
00794             ncclose(rivSurfID);
00795     }
00796 
00797 }
00798 
00799 realtype FPrint_CS_AreaOrPerem(int rivOrder, realtype rivDepth, realtype rivCoeff, realtype a_pBool)
00801 
00806 {
00807     realtype rivArea, rivPerem, eq_Wid, EPSILON=0.05;
00808     switch(rivOrder)
00809     {
00810         case 1:
00811             rivArea = rivDepth*rivCoeff;
00812             rivPerem= 2.0*rivDepth+rivCoeff;
00813             eq_Wid=rivCoeff;
00814             return (a_pBool==1?rivArea:(a_pBool==2?rivPerem:eq_Wid)); //returnVal1(rivArea, rivPerem, eq_Wid, a_pBool);
00815         case 2:
00816             rivArea = pow(rivDepth,2)/rivCoeff;
00817             rivPerem = 2.0*rivDepth*pow(1+pow(rivCoeff,2),0.5)/rivCoeff;
00818             eq_Wid=2.0*pow(rivDepth+EPSILON,1/(rivOrder-1))/pow(rivCoeff,1/(rivOrder-1));
00819             return (a_pBool==1?rivArea:(a_pBool==2?rivPerem:eq_Wid)); //returnVal1(rivArea, rivPerem, eq_Wid, a_pBool);
00820         case 3:
00821             rivArea = 4*pow(rivDepth,1.5)/(3*pow(rivCoeff,0.5));
00822             rivPerem =(pow(rivDepth*(1+4*rivCoeff*rivDepth)/rivCoeff,0.5))+(log(2*pow(rivCoeff*rivDepth,0.5)+pow(1+4*rivCoeff*rivDepth,0.5))/(2*rivCoeff));
00823             eq_Wid=2.0*pow(rivDepth+EPSILON,1/(rivOrder-1))/pow(rivCoeff,1/(rivOrder-1));
00824             return (a_pBool==1?rivArea:(a_pBool==2?rivPerem:eq_Wid)); //returnVal1(rivArea, rivPerem, eq_Wid, a_pBool);
00825         case 4:
00826             rivArea = 3*pow(rivDepth,4.0/3.0)/(2*pow(rivCoeff,1.0/3.0));
00827             rivPerem = 2*((pow(rivDepth*(1+9*pow(rivCoeff,2.0/3.0)*rivDepth),0.5)/3)+(log(3*pow(rivCoeff,1.0/3.0)*pow(rivDepth,0.5)+pow(1+9*pow(rivCoeff,2.0/3.0)*rivDepth,0.5))/(9*pow(rivCoeff,1.0/3.0))));
00828             eq_Wid=2.0*pow(rivDepth+EPSILON,1/(rivOrder-1))/pow(rivCoeff,1/(rivOrder-1));
00829             return (a_pBool==1?rivArea:(a_pBool==2?rivPerem:eq_Wid)); //returnVal1(rivArea, rivPerem, eq_Wid, a_pBool);
00830         default:
00831             printf("\n Relevant Values entered are wrong");
00832             printf("\n Depth: %lf\tCoeff: %lf\tOrder: %d\t");
00833             return 0;
00834     }
00835 }
00836 
00837 realtype FPrint_OverlandFlow(int loci, int locj, int surfmode, realtype avg_y, realtype grad_y, realtype avg_sf, realtype alfa, realtype beta, realtype crossA, realtype avg_rough, int eletypeBool, realtype avg_perem) 
00838 
00839 
00852 {
00853     realtype flux;
00854     int locBool;
00855     float hydRadius;
00856     /* if surface gradient is not enough to overcome the friction */
00857     if(fabs(grad_y) <= avg_sf)
00858     {
00859          flux = 0;
00860     }
00861     else
00862     {
00863          if(grad_y > 0)
00864          {
00865               locBool=1;
00866          }
00867          else if(grad_y < 0)
00868          {
00869              locBool=-1;
00870          }
00871          switch(surfmode)
00872          {
00873                 case 1:
00874                   if(eletypeBool==1)
00875                   {
00876                        /* Kinematic Wave Approximation constitutive relationship: Manning Equation */
00877                        alfa = sqrt(locBool*grad_y)/avg_rough;
00878                        beta = pow(avg_y, 2.0/3.0);
00879                        flux = locBool*alfa*beta*crossA;
00880                        break;
00881                   }
00882                   else
00883                   {
00884                        /*alfa = sqrt(locBool*grad_y)/(avg_rough*pow((avg_perem>0?avg_perem:0), 2.0/3.0));
00885                        beta = 5.0/3.0;
00886                        *flux = locBool*alfa*pow(crossA, beta);
00887                        */
00888                        hydRadius = (avg_perem>0?crossA/avg_perem:0);
00889                        flux = locBool*sqrt(locBool*grad_y)*crossA*pow(hydRadius,2.0/3.0)/avg_rough;
00890                        break;
00891                   }
00892                 case 2:
00893                   if(eletypeBool==1)
00894                   {
00895                        /* Diffusion Wave Approximation constitutive relationship: Gottardi & Venutelli, 1993 */
00896                        alfa = pow(pow(avg_y, 1.0/3.0),2)/(1.0*avg_rough);
00897                        beta = alfa;
00898                        flux = locBool*crossA*beta*sqrt(locBool*grad_y);
00899                        break;
00900                   }
00901                   else
00902                   {
00903                        /*alfa = pow(pow(avg_y, 1.0/3.0),2)/(1.0*avg_rough);
00904                        beta = alfa;
00905                        flux = locBool*crossA*beta*sqrt(locBool*grad_y);
00906                        */
00907                        hydRadius = (avg_perem>0?crossA/avg_perem:0);
00908                        flux = locBool*sqrt(locBool*grad_y)*crossA*pow(hydRadius,2.0/3.0)/(1.0*avg_rough);
00909                        break;
00910                   }
00911                 default:
00912                   if(eletypeBool==1)
00913                   {
00914                           printf("Fatal Error: Surface Overland Mode Type Is Wrong!");
00915                   }
00916                   else
00917                   {
00918                        printf("Fatal Error: River Routing Mode Type Is Wrong!");
00919                   }
00920                   exit(1);
00921         }
00922     }
00923     return flux;
00924 }
00925 
00926 /*    Function to print River Flow in TXT format    */
00927 void  printRiverFlow(Model_Data mData, N_Vector CV_Y, FILE *flow_file, realtype t)
00929 
00934 {
00935     int i;
00936     realtype TotalY_Riv, Perem, TotalY_Riv_down, Perem_down, Avg_Perem, Avg_Y_Riv, Avg_Rough, Distance, Dif_Y_Riv, Avg_Sf, CrossA, Alfa, Beta;
00937     realtype Flux;
00938     for(i=0; i<mData->NumRiv; i++)
00939     {
00940         TotalY_Riv = NV_Ith_S(CV_Y, i + 3*mData->NumEle) + mData->Riv[i].zmin;
00941         Perem = FPrint_CS_AreaOrPerem(mData->Riv_Shape[mData->Riv[i].shape - 1].interpOrd,NV_Ith_S(CV_Y, i + 3*mData->NumEle),mData->Riv_Shape[mData->Riv[i].shape - 1].coeff,2);
00942         if(mData->Riv[i].down > 0)
00943         {
00944             TotalY_Riv_down = NV_Ith_S(CV_Y, mData->Riv[i].down - 1 + 3*mData->NumEle)  + mData->Riv[mData->Riv[i].down - 1].zmin;
00945             Perem_down = FPrint_CS_AreaOrPerem(mData->Riv_Shape[mData->Riv[mData->Riv[i].down - 1].shape - 1].interpOrd,NV_Ith_S(CV_Y, mData->Riv[i].down - 1 + 3*mData->NumEle),mData->Riv_Shape[mData->Riv[mData->Riv[i].down - 1].shape - 1].coeff,2);
00946             Avg_Perem = (Perem + Perem_down)/2.0;    
00947             if(mData->Riv[mData->Riv[i].down - 1].zmin>mData->Riv[i].zmin)
00948             {
00949                 if(mData->Riv[mData->Riv[i].down - 1].zmin > TotalY_Riv)
00950                 {
00951                     Avg_Y_Riv=NV_Ith_S(CV_Y, mData->Riv[i].down - 1 + 3*mData->NumEle)/2;
00952                 }
00953                 else
00954                 {
00955                     Avg_Y_Riv=(TotalY_Riv-mData->Riv[mData->Riv[i].down - 1].zmin+NV_Ith_S(CV_Y, mData->Riv[i].down - 1 + 3*mData->NumEle))/2;
00956                 }
00957             }
00958             else
00959             {
00960                 if(mData->Riv[i].zmin>TotalY_Riv_down)
00961                 {
00962                     Avg_Y_Riv=NV_Ith_S(CV_Y, i + 3*mData->NumEle)/2;
00963                 }
00964                 else
00965                 {
00966                     Avg_Y_Riv=(NV_Ith_S(CV_Y, i + 3*mData->NumEle)+TotalY_Riv_down-mData->Riv[i].zmin)/2;
00967                 }
00968             }
00969             Avg_Rough = (mData->Riv_Mat[mData->Riv[i].material - 1].Rough + mData->Riv_Mat[mData->Riv[mData->Riv[i].down - 1].material-1].Rough)/2.0;
00970 
00971             Distance = (mData->Riv[i].Length+mData->Riv[mData->Riv[i].down - 1].Length)/2;
00972 
00973             Dif_Y_Riv = (TotalY_Riv - TotalY_Riv_down)/Distance;
00974 
00975             Avg_Sf = (mData->Riv_Mat[mData->Riv[i].material - 1].Sf + mData->Riv_Mat[mData->Riv[mData->Riv[i].down - 1].material-1].Sf)/2.0;
00976             /*CrossA = 0.5*(FPrint_CS_AreaOrPerem(mData->Riv_Shape[mData->Riv[i].shape - 1].interpOrd,DummyY[i + 3*MD->NumEle],MD->Riv_Shape[MD->Riv[i].shape - 1].coeff,1)+CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[MD->Riv[i].down - 1].shape - 1].interpOrd,DummyY[MD->Riv[i].down - 1 + 3*MD->NumEle],MD->Riv_Shape[MD->Riv[MD->Riv[i].down - 1].shape - 1].coeff,1));*/
00977             CrossA = FPrint_CS_AreaOrPerem(mData->Riv_Shape[mData->Riv[i].shape - 1].interpOrd,Avg_Y_Riv,mData->Riv_Shape[mData->Riv[i].shape - 1].coeff,1);
00978             Flux = FPrint_OverlandFlow(i,1,mData->RivMode, Avg_Y_Riv,Dif_Y_Riv,Avg_Sf,Alfa,Beta,CrossA,Avg_Rough,0,Avg_Perem);
00979             /* Correction is being done in flux terms which can be > 0 even when there is no source water level present */
00980             if(NV_Ith_S(CV_Y, i + 3*mData->NumEle) <= 0 && Flux > 0)
00981             {
00982                 Flux = 0.0;
00983             }
00984             else if(NV_Ith_S(CV_Y, mData->Riv[i].down - 1 + 3*mData->NumEle) <= 0 && Flux < 0)
00985             {
00986                 Flux = 0.0;
00987             }
00988         }
00989         else
00990         {
00991             Flux = 0.0;
00992         }
00993 
00994         tempFlow[i]+=Flux/RivFlowT;
00995         if(((int) t)%RivFlowT==0){
00996             fprintf(flow_file, "%lf\t", tempFlow[i]);
00997             tempFlow[i]=0.0;
00998         }
00999     }
01000 
01001     if(((int) t)%RivFlowT==0){
01002         fprintf(flow_file, "\n");
01003     }
01004 }
01005 
01006 /*    Function to print River Flow in CDF format    */
01007 void  printRiverFlowcdf(Model_Data mData, N_Vector CV_Y, int ncid, int data_varid, realtype t)
01009 
01015 {
01016     int i;
01017     static int call=0;
01018     float *data;
01019     realtype TotalY_Riv, Perem, TotalY_Riv_down, Perem_down, Avg_Perem, Avg_Y_Riv, Avg_Rough, Distance, Dif_Y_Riv, Avg_Sf, CrossA, Alfa, Beta;
01020     realtype Flux;
01021     for(i=0; i<mData->NumRiv; i++)
01022     {
01023         TotalY_Riv = NV_Ith_S(CV_Y, i + 3*mData->NumEle) + mData->Riv[i].zmin;
01024         Perem = FPrint_CS_AreaOrPerem(mData->Riv_Shape[mData->Riv[i].shape - 1].interpOrd,NV_Ith_S(CV_Y, i + 3*mData->NumEle),mData->Riv_Shape[mData->Riv[i].shape - 1].coeff,2);
01025         if(mData->Riv[i].down > 0)
01026         {
01027             TotalY_Riv_down = NV_Ith_S(CV_Y, mData->Riv[i].down - 1 + 3*mData->NumEle)  + mData->Riv[mData->Riv[i].down - 1].zmin;
01028             Perem_down = FPrint_CS_AreaOrPerem(mData->Riv_Shape[mData->Riv[mData->Riv[i].down - 1].shape - 1].interpOrd,NV_Ith_S(CV_Y, mData->Riv[i].down - 1 + 3*mData->NumEle),mData->Riv_Shape[mData->Riv[mData->Riv[i].down - 1].shape - 1].coeff,2);
01029             Avg_Perem = (Perem + Perem_down)/2.0;    
01030             if(mData->Riv[mData->Riv[i].down - 1].zmin>mData->Riv[i].zmin)
01031             {
01032                 if(mData->Riv[mData->Riv[i].down - 1].zmin > TotalY_Riv)
01033                 {
01034                     Avg_Y_Riv=NV_Ith_S(CV_Y, mData->Riv[i].down - 1 + 3*mData->NumEle)/2;
01035                 }
01036                 else
01037                 {
01038                     Avg_Y_Riv=(TotalY_Riv-mData->Riv[mData->Riv[i].down - 1].zmin+NV_Ith_S(CV_Y, mData->Riv[i].down - 1 + 3*mData->NumEle))/2;
01039                 }
01040             }
01041             else
01042             {
01043                 if(mData->Riv[i].zmin>TotalY_Riv_down)
01044                 {
01045                     Avg_Y_Riv=NV_Ith_S(CV_Y, i + 3*mData->NumEle)/2;
01046                 }
01047                 else
01048                 {
01049                     Avg_Y_Riv=(NV_Ith_S(CV_Y, i + 3*mData->NumEle)+TotalY_Riv_down-mData->Riv[i].zmin)/2;
01050                 }
01051             }
01052             Avg_Rough = (mData->Riv_Mat[mData->Riv[i].material - 1].Rough + mData->Riv_Mat[mData->Riv[mData->Riv[i].down - 1].material-1].Rough)/2.0;
01053 
01054             Distance = (mData->Riv[i].Length+mData->Riv[mData->Riv[i].down - 1].Length)/2;
01055 
01056             Dif_Y_Riv = (TotalY_Riv - TotalY_Riv_down)/Distance;
01057 
01058             Avg_Sf = (mData->Riv_Mat[mData->Riv[i].material - 1].Sf + mData->Riv_Mat[mData->Riv[mData->Riv[i].down - 1].material-1].Sf)/2.0;
01059             /*CrossA = 0.5*(FPrint_CS_AreaOrPerem(mData->Riv_Shape[mData->Riv[i].shape - 1].interpOrd,DummyY[i + 3*MD->NumEle],MD->Riv_Shape[MD->Riv[i].shape - 1].coeff,1)+CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[MD->Riv[i].down - 1].shape - 1].interpOrd,DummyY[MD->Riv[i].down - 1 + 3*MD->NumEle],MD->Riv_Shape[MD->Riv[MD->Riv[i].down - 1].shape - 1].coeff,1));*/
01060             CrossA = FPrint_CS_AreaOrPerem(mData->Riv_Shape[mData->Riv[i].shape - 1].interpOrd,Avg_Y_Riv,mData->Riv_Shape[mData->Riv[i].shape - 1].coeff,1);
01061             Flux = FPrint_OverlandFlow(i,1,mData->RivMode, Avg_Y_Riv,Dif_Y_Riv,Avg_Sf,Alfa,Beta,CrossA,Avg_Rough,0,Avg_Perem);
01062             /* Correction is being done in flux terms which can be > 0 even when there is no source water level present */
01063             if(NV_Ith_S(CV_Y, i + 3*mData->NumEle) <= 0 && Flux > 0)
01064             {
01065                 Flux = 0.0;
01066             }
01067             else if(NV_Ith_S(CV_Y, mData->Riv[i].down - 1 + 3*mData->NumEle) <= 0 && Flux < 0)
01068             {
01069                 Flux = 0.0;
01070             }
01071         }
01072         else
01073         {
01074             Flux = 0.0;
01075         }
01076 
01077         tempFlow[i]+=Flux/RivFlowT;
01078 
01079     }
01080 
01081     if(((int) t)%RivFlowT==0){
01082         startRiv[0]=call++;
01083         if((retval = nc_put_vara_double(ncid, data_varid, startRiv, countRiv, &tempFlow[0])))
01084             ERR(retval);
01085         for(i=0; i<mData->NumRiv; i++)
01086             tempFlow[i]=0.0;
01087     }
01088 }
01089 
01090 /*    Function to print Baseflow to River in TXT format    */
01091 void printRiverBase(Model_Data mData, FILE *rivBaseFile, realtype t)
01093 
01097 {
01098         int i;
01099     for(i=0; i<mData->NumRiv; i++){
01100         tempBase[i]+=(mData->FluxRiv[i][4]+mData->FluxRiv[i][5])/RivBaseT;
01101         if(((int) t)%RivBaseT==0){
01102             fprintf(rivBaseFile, "%lf\t", tempBase[i]);
01103             tempBase[i]=0.0;
01104         }
01105     }
01106     if(((int) t)%RivBaseT==0){
01107         fprintf(rivBaseFile, "\n");
01108     }
01109 }
01110 
01111 /*    Function to print Base flow to River in CDF format    */
01112 void printRiverBasecdf(Model_Data mData, int ncid, int data_varid, realtype t)
01114 
01119 {
01120         int i;
01121     static int call=0;
01122     float *data;
01123     for(i=0; i<mData->NumRiv; i++){
01124         tempBase[i]+=(mData->FluxRiv[i][4]+mData->FluxRiv[i][5])/RivBaseT;
01125     }
01126     if(((int) t)%RivBaseT==0){
01127         startRiv[0]=call++;
01128         if((retval = nc_put_vara_double(ncid, data_varid, startRiv, countRiv, &tempBase[0])))
01129             ERR(retval);
01130         for(i=0; i<mData->NumRiv; i++)
01131             tempBase[i]=0.0;
01132     }
01133 }
01134 
01135 /*    Function to print Surfaceflow to River in TXT format    */
01136 void printRiverSurf(Model_Data mData, FILE *rivSurfFile, realtype t)
01138 
01142 {
01143     int i;
01144     for(i=0; i<mData->NumRiv; i++){
01145         tempSurf[i]+=(mData->FluxRiv[i][2]+mData->FluxRiv[i][3])/RivSurfT;
01146         if(((int) t)%RivSurfT==0){
01147             fprintf(rivSurfFile, "%lf\t", tempSurf[i]);
01148             tempSurf[i]=0.0;
01149         }
01150     }
01151     if(((int) t)%RivSurfT==0){
01152         fprintf(rivSurfFile, "\n");
01153     }
01154 }
01155 
01156 /*    Function to print Surfaceflow to River in CDF format    */
01157 void printRiverSurfcdf(Model_Data mData, int ncid, int data_varid, realtype t)
01159 
01164 {
01165     int i;
01166     static int call=0;
01167     float *data;
01168     for(i=0; i<mData->NumRiv; i++){
01169         tempSurf[i]+=(mData->FluxRiv[i][2]+mData->FluxRiv[i][3])/RivSurfT;
01170     }
01171     if(((int) t)%RivSurfT==0){
01172         startRiv[0]=call++;
01173         if((retval = nc_put_vara_double(ncid, data_varid, startRiv, countRiv, &tempSurf[0])))
01174             ERR(retval);
01175         for(i=0; i<mData->NumRiv; i++)
01176             tempSurf[i]=0.0;
01177     }
01178 }
01179 
01180 /*    Function to print River State (head) in TXT format    */
01181 void printRiverHead(Model_Data mData, N_Vector CV_Y, FILE *rivHeadFile, realtype t)
01183 
01188 {
01189     int i;
01190     for(i=0; i<mData->NumRiv; i++){
01191         tempHead[i]+=NV_Ith_S(CV_Y, 3*mData->NumEle + i)/RivHeadT;
01192         if(((int) t)%RivHeadT==0){
01193             fprintf(rivHeadFile, "%lf\t", tempHead[i]);
01194             tempHead[i]=0.0;
01195         }
01196     }
01197     if(((int) t)%RivHeadT==0){
01198         fprintf(rivHeadFile, "\n");
01199     }
01200 }
01201 
01202 /*    Function to print River State (head) in CDF format    */
01203 void printRiverHeadcdf(Model_Data mData, N_Vector CV_Y, int ncid, int data_varid, realtype t)
01205 
01211 {
01212     int i;
01213     static int call=0;
01214     float *data;
01215     for(i=0; i<mData->NumRiv; i++){
01216         tempHead[i]+=NV_Ith_S(CV_Y, 3*mData->NumEle + i)/RivHeadT;
01217     }
01218     if(((int) t)%RivHeadT==0){
01219         startRiv[0]=call++;
01220         if((retval = nc_put_vara_double(ncid, data_varid, startRiv, countRiv, &tempHead[0])))
01221             ERR(retval);
01222         for(i=0; i<mData->NumRiv; i++)
01223             tempHead[i]=0.0;
01224     }
01225 }
01226 
01227 /*    Function to print Interception Storage in TXT format    */
01228 void printIS(Model_Data mData, FILE *isFile, realtype t)
01230 
01234 {
01235     int i;
01236     for(i=0; i<mData->NumEle; i++){
01237         tempIS[i]+=mData->EleIS[i]/ISStateT;
01238         if(((int) t)%ISStateT==0){
01239             fprintf(isFile, "%lf\t", tempIS[i]);
01240             tempIS[i]=0.0;
01241         }
01242     }
01243     if(((int) t)%ISStateT==0){
01244         fprintf(isFile, "\n");
01245     }
01246 }
01247 
01248 /*    Function to print Interception Storage in CDF format    */
01249 void printIScdf(Model_Data mData, int ncid, int data_varid, realtype t)
01251 
01256 {
01257     int i;
01258     static int call=0;
01259     float *data;
01260     for(i=0; i<mData->NumEle; i++){
01261         tempIS[i]+=mData->EleIS[i]/ISStateT;
01262     }
01263     if(((int) t)%ISStateT==0){
01264         startEle[0]=call++;
01265         if((retval = nc_put_vara_double(ncid, data_varid, startEle, countEle, &tempIS[0])))
01266             ERR(retval);
01267         for(i=0; i<mData->NumEle; i++)
01268             tempIS[i]=0.0;
01269     }
01270 }
01271 
01272 /*    Function to print Saturated State (head) in TXT format    */
01273 void printSatState(Model_Data mData, N_Vector CV_Y, FILE *file, realtype t)
01275 
01280 {
01281     int i;
01282     for(i=0; i<mData->NumEle; i++){
01283         tempSatState[i]+=NV_Ith_S(CV_Y, 2*mData->NumEle + i)/SatStateT;
01284         if(((int) t)%SatStateT==0){
01285             fprintf(file, "%lf\t", tempSatState[i]);
01286             tempSatState[i]=0.0;
01287         }
01288     }
01289     if(((int) t)%SatStateT==0){
01290         fprintf(file, "\n");
01291     }
01292 }
01293 
01294 /*    Function to print Saturated State (head) in CDF format    */
01295 void printSatStatecdf(Model_Data mData, N_Vector CV_Y, int ncid, int data_varid, realtype t)
01297 
01303 {
01304     int i;
01305     static int call=0;
01306     float *data;
01307     for(i=0; i<mData->NumEle; i++){
01308         tempSatState[i]+=NV_Ith_S(CV_Y, 2*mData->NumEle + i)/SatStateT;
01309     }
01310     if(((int) t)%SatStateT==0){
01311         startEle[0]=call++;
01312         if((retval = nc_put_vara_double(ncid, data_varid, startEle, countEle, &tempSatState[0])))
01313             ERR(retval);
01314         for(i=0; i<mData->NumEle; i++)
01315             tempSatState[i]=0.0;
01316     }
01317 }
01318 
01319 /*    Function to print Unsaturated State (head) in TXT format    */
01320 void printUsatState(Model_Data mData, N_Vector CV_Y, FILE *file, realtype t)
01322 
01327 {
01328     int i;
01329     for(i=0; i<mData->NumEle; i++){
01330         tempUsatState[i]+=NV_Ith_S(CV_Y, 1*mData->NumEle + i)/UsatStateT;
01331         if(((int) t)%UsatStateT==0){
01332             fprintf(file, "%lf\t", tempUsatState[i]);
01333             tempUsatState[i]=0.0;
01334         }
01335     }
01336     if(((int) t)%UsatStateT==0){
01337         fprintf(file, "\n");
01338     }
01339 }
01340 
01341 /*    Function to print Unsaturated State (head) in CDF format    */
01342 void printUsatStatecdf(Model_Data mData, N_Vector CV_Y, int ncid, int data_varid, realtype t)
01344 
01350 {
01351     int i;
01352     static int call=0;
01353     float *data;
01354     for(i=0; i<mData->NumEle; i++){
01355         tempUsatState[i]+=NV_Ith_S(CV_Y, 1*mData->NumEle + i)/UsatStateT;
01356     }
01357     if(((int) t)%UsatStateT==0){
01358         startEle[0]=call++;
01359         if((retval = nc_put_vara_double(ncid, data_varid, startEle, countEle, &tempUsatState[0])))
01360             ERR(retval);
01361         for(i=0; i<mData->NumEle; i++)
01362             tempUsatState[i]=0.0;
01363     }
01364 }
01365 
01366 /*    Function to print Surface Flow State (head) in TXT format    */
01367 void printSurfState(Model_Data mData, N_Vector CV_Y, FILE *file, realtype t)
01369 
01374 {
01375     int i;
01376     for(i=0; i<mData->NumEle; i++){
01377         tempSurfState[i]+=NV_Ith_S(CV_Y, i)/SurfStateT;
01378         if(((int) t)%SurfStateT==0){
01379             fprintf(file, "%lf\t", tempSurfState[i]);
01380             tempSurfState[i]=0.0;
01381         }
01382     }
01383     if(((int) t)%SurfStateT==0){
01384         fprintf(file, "\n");
01385     }
01386 }
01387 
01388 /*    Function to print Surface Flow State (head) in CDF format    */
01389 void printSurfStatecdf(Model_Data mData, N_Vector CV_Y, int ncid, int data_varid, realtype t)
01391 
01397 {
01398     int i;
01399     static int call=0;
01400     float *data;
01401     for(i=0; i<mData->NumEle; i++){
01402         tempSurfState[i]+=NV_Ith_S(CV_Y, i)/SurfStateT;
01403     }
01404     if(((int) t)%SurfStateT==0){
01405         startEle[0]=call++;
01406         if((retval = nc_put_vara_double(ncid, data_varid, startEle, countEle, &tempSurfState[0])))
01407             ERR(retval);
01408         for(i=0; i<mData->NumEle; i++)
01409             tempSurfState[i]=0.0;
01410     }
01411 }
01412 
01413 /*    Function to print ET0 in TXT format    */
01414 void printET0(Model_Data mData, FILE *file, realtype t)
01416 
01420 {
01421     int i;
01422     for(i=0; i<mData->NumEle; i++){
01423         tempET0[i]+=(mData->EleET[i][0]*mData->Ele[i].VegFrac)/ET0T;
01424         if(((int) t)%ET0T==0){
01425             fprintf(file, "%lf\t", tempET0[i]);
01426             tempET0[i]=0.0;
01427         }
01428     }
01429     if(((int) t)%ET0T==0){
01430         fprintf(file, "\n");
01431     }
01432 }
01433 
01434 /*    Function to print ET0 in CDF format    */
01435 void printET0cdf(Model_Data mData, int ncid, int data_varid, realtype t)
01437 
01442 {
01443     int i;
01444     static int call=0;
01445     float *data;
01446     for(i=0; i<mData->NumEle; i++){
01447         tempET0[i]+=(mData->EleET[i][0]*mData->Ele[i].VegFrac)/ET0T;
01448     }
01449     if(((int) t)%ET0T==0){
01450         startEle[0]=call++;
01451         if((retval = nc_put_vara_double(ncid, data_varid, startEle, countEle, &tempET0[0])))
01452             ERR(retval);
01453         for(i=0; i<mData->NumEle; i++)
01454             tempET0[i]=0.0;
01455     }
01456 }
01457 
01458 /*    Function to print ET1 in TXT format    */
01459 void printET1(Model_Data mData, FILE *file, realtype t)
01461 
01465 {
01466     int i;
01467     for(i=0; i<mData->NumEle; i++){
01468         tempET1[i]+=mData->EleET[i][1]/ET1T;
01469         if(((int) t)%ET1T==0){
01470             fprintf(file, "%lf\t", tempET1[i]);
01471             tempET1[i]=0.0;
01472         }
01473     }
01474     if(((int) t)%ET1T==0){
01475         fprintf(file, "\n");
01476     }
01477 }
01478 
01479 /*    Function to print ET1 in CDF format    */
01480 void printET1cdf(Model_Data mData, int ncid, int data_varid, realtype t)
01482 
01487 {
01488     int i;
01489     static int call=0;
01490     float *data;
01491     data = (float *)malloc(mData->NumEle * sizeof(float));
01492     for(i=0; i<mData->NumEle; i++){
01493         tempET1[i]+=mData->EleET[i][1]/ET1T;
01494     }
01495     if(((int) t)%ET1T==0){
01496         startEle[0]=call++;
01497         if((retval = nc_put_vara_double(ncid, data_varid, startEle, countEle, &tempET1[0])))
01498             ERR(retval);
01499         for(i=0; i<mData->NumEle; i++)
01500             tempET1[i]=0.0;
01501     }
01502 }
01503 
01504 /*    Function to print ET2 in TXT format    */
01505 void printET2(Model_Data mData, FILE *file, realtype t)
01507 
01511 {
01512     int i;
01513     for(i=0; i<mData->NumEle; i++){
01514         tempET2[i]+=mData->EleET[i][2]/ET2T;
01515         if(((int) t)%ET2T==0){
01516             fprintf(file, "%lf\t", tempET2[i]);
01517             tempET2[i]=0.0;
01518         }
01519     }
01520     if(((int) t)%ET2T==0){
01521         fprintf(file, "\n");
01522     }
01523 }
01524 
01525 /*    Function to print ET2 in CDF format    */
01526 void printET2cdf(Model_Data mData, int ncid, int data_varid, realtype t)
01528 
01533 {
01534     int i;
01535     static int call=0;
01536     float *data;
01537     for(i=0; i<mData->NumEle; i++){
01538         tempET2[i]+=mData->EleET[i][2]/ET2T;
01539     }
01540     if(((int) t)%ET2T==0){
01541         startEle[0]=call++;
01542         if((retval = nc_put_vara_double(ncid, data_varid, startEle, countEle, &tempET2[0])))
01543             ERR(retval);
01544         for(i=0; i<mData->NumEle; i++)
01545             tempET2[i]=0.0;
01546     }
01547 }
01548 
01549 /*    Function to print Net Precipitation in TXT format    */
01550 void printNetPpt(Model_Data mData, FILE *file, realtype t)
01552 
01556 {
01557     int i;
01558     for(i=0; i<mData->NumEle; i++){
01559         tempNetPpt[i] += mData->EleNetPrep[i] / NetPptT;
01560         if(((int) t)%NetPptT==0){
01561             fprintf(file, "%lf\t", tempNetPpt[i]);
01562             tempNetPpt[i]=0.0;
01563         }
01564     }
01565     if(((int) t)%NetPptT==0){
01566         fprintf(file, "\n");
01567     }
01568 }
01569 
01570 /*    Function to print Net Precipitation in CDF format    */
01571 void printNetPptcdf(Model_Data mData, int ncid, int data_varid, realtype t)
01573 
01578 {
01579     int i;
01580     static int call=0;
01581     float *data;
01582     for(i=0; i<mData->NumEle; i++){
01583         tempNetPpt[i] += mData->EleNetPrep[i] / NetPptT;
01584     }
01585     if(((int) t)%NetPptT==0){
01586         startEle[0]=call++;
01587         if((retval = nc_put_vara_double(ncid, data_varid, startEle, countEle, &tempNetPpt[0])))
01588             ERR(retval);
01589         for(i=0; i<mData->NumEle; i++)
01590             tempNetPpt[i]=0.0;
01591     }
01592 }
01593 
01594 /*    Function to print Variable Infiltration in TXT format    */
01595 void printInfil(Model_Data mData, FILE *file, realtype t)
01597 
01601 {
01602     int i;
01603     for(i=0; i<mData->NumEle; i++){
01604         tempInfil[i] += mData->EleVic[i] / InfilT;
01605         if(((int) t)%InfilT==0){
01606             fprintf(file, "%lf\t", tempInfil[i]);
01607             tempInfil[i]=0.0;
01608         }
01609     }
01610     if(((int) t)%InfilT==0){
01611         fprintf(file, "\n");
01612     }
01613 }
01614 
01615 /*    Function to print Variable Infiltration in CDF format    */
01616 void printInfilcdf(Model_Data mData, int ncid, int data_varid, realtype t)
01618 
01623 {
01624     int i;
01625     static int call=0;
01626     float *data;
01627     for(i=0; i<mData->NumEle; i++){
01628         tempInfil[i] += mData->EleVic[i] / InfilT;
01629     }
01630     if(((int) t)%InfilT==0){
01631         startEle[0]=call++;
01632         if((retval = nc_put_vara_double(ncid, data_varid, startEle, countEle, &tempInfil[0])))
01633             ERR(retval);
01634         for(i=0; i<mData->NumEle; i++)
01635             tempInfil[i]=0.0;
01636     }
01637 }
01638 
01639 /*    Function to print Recharge to GW in TXT format    */
01640 void printRecharge(Model_Data mData, FILE *file, realtype t)
01642 
01646 {
01647     int i;
01648     for(i=0; i<mData->NumEle; i++){
01649         tempRecharge[i] += mData->Recharge[i] / RECHARGET;
01650         if(((int) t)%RECHARGET==0){
01651             fprintf(file, "%lf\t", tempRecharge[i]);
01652             tempRecharge[i]=0.0;
01653         }
01654     }
01655     if(((int) t)%RECHARGET==0){
01656         fprintf(file, "\n");
01657     }
01658 }
01659 
01660 /*    Function to print Recharge to GW in CDF format    */
01661 void printRechargecdf(Model_Data mData, int ncid, int data_varid, realtype t)
01663 
01668 {
01669     int i;
01670     static int call=0;
01671     float *data;
01672     for(i=0; i<mData->NumEle; i++){
01673         tempRecharge[i] += mData->Recharge[i] / RECHARGET;
01674     }
01675     if(((int) t)%RECHARGET==0){
01676         startEle[0]=call++;
01677         if((retval = nc_put_vara_double(ncid, data_varid, startEle, countEle, &tempRecharge[0])))
01678             ERR(retval);
01679         for(i=0; i<mData->NumEle; i++)
01680             tempRecharge[i]=0.0;
01681     }
01682 }
Generated on Thu Jul 12 14:34:19 2007 for PIHM by  doxygen 1.5.2