f.c

Go to the documentation of this file.
00001 /*******************************************************************************
00002  * File        : f.c                                                           *
00003  * Function    : calculate fluxes, right hand side and construct ODE system    *
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  * This is the model kernel. Idealy, Users can modify this kernel based on     *
00011  * hydrological conceptual model. Also, flexible constutive relationships can  *
00012  * apply. This is a multiple scale model in that                               *
00013  *   1. All equations solved are in integrated form, which make possible to    *
00014  *      use larger elements;                                                   *
00015  *   2. Flexible constitutive relationships can be applied;                    *
00016  *   3. Users can modify this kernel based on appropriate conceptual model.    *
00017  * We have prove that in finer limit, this model works as good as any other    *
00018  * tranditional finite element/difference model. However, for larger scale     *
00019  * application, this model is possible to be applied due to above features.    *
00020  *                                                                             *
00021  * This code is free for users with research purpose only, if appropriate      *
00022  * citation is refered. However, there is no warranty in any format for this   *
00023  * product.                                                                    *
00024  *                                                                             *
00025  * For questions or comments, please contact the authors of the reference.     *
00026  * One who want to use it for other consideration may also contact Dr.Duffy    *
00027  * at cxd11@psu.edu.                                                           *
00028  *******************************************************************************/
00029 
00031 
00032 /*    C Header Files    */
00033 #include <stdio.h>
00034 #include <stdlib.h>
00035 #include <math.h>
00036 #include <string.h>
00037 
00038 /*    SUNDIALS Header Files    */
00039 #include "nvector_serial.h"
00040 #include "sundials_types.h"
00041 
00042 /*    PIHM Header Files    */
00043 #include "pihm.h"
00044 #include "calib.h"
00045 
00046 #define EPSILON 0.05
00047 #define THRESH  0.00
00048 
00049 
00050 #define GRAV  73231257600.0     
00052 #define ABS_TOL    1E-4         
00054 realtype Vic_CALIB;
00055 realtype rivK_CALIB;
00056 realtype Kh_CALIB;
00057 realtype Rec_CALIB;
00058 realtype et2_CALIB;
00059 realtype et1_CALIB;
00060 realtype sat_THRESH;
00061 realtype mp_MULTFH;
00062 realtype mp_MULTFV;
00063 realtype mpArea_CALIB;
00064 realtype ovl_THRESH_H;
00065 realtype ovl_THRESH_V;
00066 realtype rzd_CALIB;
00067 
00068 
00069 /*    Function Declarations    */
00070 realtype Interpolation(TSD *Data, realtype t);
00071 realtype returnVal(realtype rArea, realtype rPerem, realtype eqWid, realtype ap_Bool);
00072 realtype CS_AreaOrPerem(int rivOrder, realtype rivDepth, realtype rivCoeff, realtype a_pBool);
00073 void OverlandFlow(realtype **flux, 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);
00074 void OLflowFromEleToRiv(realtype sideEle_y,realtype sideEle_zmax,realtype rivX,realtype sideEleX,realtype rivY,realtype sideEleY,realtype cwr,realtype rivZmax,realtype loc_yriver,realtype **fluxriv,int loc_i,int loc_j,realtype length);
00075 void GWflowFromEleToRiv(realtype sideEle_y,realtype sideEle_zmax, realtype sideEle_zmin,realtype rivX,realtype sideEleX,realtype rivY,realtype sideEleY,int loc_McPore,realtype loc_yriver,realtype loc_totyriver,realtype **fluxriv,int loc_i,int loc_j,realtype length, realtype loc_base,realtype loc_gama, realtype loc_perem,realtype loc_ksat,realtype ele_Thresh);
00076 
00077 
00078 
00079 
00080 
00081 
00082 
00083 int f(realtype t, N_Vector CV_Y, N_Vector CV_Ydot, void *DS)
00085 
00090 {
00091     int i, j;
00092 
00093     realtype Delta, Gamma;
00094     realtype Rn, T, Vel, RH, VP,P,LAI,zero_dh,cnpy_h,rl,r_a,r_s,alpha_r,f_r,eta_s,beta_s,Rmax;
00095 
00096 
00097     realtype Avg_Y_Surf, Dif_Y_Surf;
00098     realtype Distance;
00099     realtype Grad_Y_Surf, Avg_Sf;
00100 
00101     realtype TotalY_Riv, TotalY_Riv_down;
00102     realtype Perem, Perem_down, loc_perem;
00103     realtype Avg_Perem;
00104     realtype Avg_Y_Riv, Avg_Rough;
00105     realtype Left_Ele_Y, Left_Ele_YH;
00106     realtype Right_Ele_Y, Right_Ele_YH;
00107     realtype Dif_Y_Riv;
00108     realtype Avg_Y_Sub, Dif_Y_Sub;
00109     realtype Avg_Ksat, Grad_Y_Sub;
00110     realtype mp_factor,m_factor;
00111     realtype G, GI;
00112     realtype Cwr, RivPrep;
00113     realtype temp1, temp2;
00114 
00115     realtype Alfa, Beta, CrossA;
00116     realtype bank_ele;
00117     realtype AquiferDepth, Deficit, PH,elemSatn,eleSatn;
00118     realtype loc_bcEle, Avg_BedDepth;
00119 
00120     realtype *Y, *DY,*DummyY,*DummyDY;
00121     Model_Data MD;
00122 
00123     Vic_CALIB=setVic_CALIB();
00124     rivK_CALIB=setrivK_CALIB();
00125     Kh_CALIB=setKh_CALIB();
00126     Rec_CALIB=setRec_CALIB();
00127     et2_CALIB=setet2_CALIB();
00128     et1_CALIB=setet1_CALIB();
00129     sat_THRESH=setsat_THRESH();
00130     mp_MULTFH=setmp_MULTFH();
00131     mp_MULTFV=setmp_MULTFV();
00132     mpArea_CALIB=setmpArea_CALIB();
00133     ovl_THRESH_H=setovl_THRESH_H();
00134     ovl_THRESH_V=setovl_THRESH_V();
00135     rzd_CALIB=setrzd_CALIB();
00136 
00137 
00138     Y = NV_DATA_S(CV_Y);
00139     DY = NV_DATA_S(CV_Ydot);
00140     MD = (Model_Data) DS;
00141 
00142     DummyY=(realtype *)malloc((3*MD->NumEle+MD->NumRiv)*sizeof(realtype));
00143     DummyDY=(realtype *)malloc((3*MD->NumEle+MD->NumRiv)*sizeof(realtype));
00144 
00145 
00146     /* Lateral Flux Calculation Follows */
00147     for(i=0; i<3*MD->NumEle+MD->NumRiv; i++)
00148     {
00149 
00150         //printf("%.35lf %d %lf\n",t, i, Y[i]);
00151 
00152          DummyDY[i]=0.0;
00153            if(Y[i]<=ABS_TOL)
00154            {
00155               DummyY[i]=0;
00156          }
00157            else
00158            {
00159                 DummyY[i]=Y[i];
00160          }
00161     }
00162     for(i=MD->NumEle; i<2*MD->NumEle; i++)
00163     {
00164           //      if((i>=MD->NumEle)&&(i<2*MD->NumEle)){
00165         if(Y[i]>(MD->Ele[i-MD->NumEle].zmax-MD->Ele[i-MD->NumEle].zmin))
00166         {
00167             DummyY[i]=MD->Ele[i-MD->NumEle].zmax - MD->Ele[i-MD->NumEle].zmin;
00168         }
00169          if(Y[i+MD->NumEle]>(MD->Ele[i-MD->NumEle].zmax-MD->Ele[i-MD->NumEle].zmin))
00170          {
00171             DummyY[i+MD->NumEle]=MD->Ele[i-MD->NumEle].zmax - MD->Ele[i-MD->NumEle].zmin;
00172         }
00173         if(DummyY[i+MD->NumEle]+DummyY[i]>(MD->Ele[i-MD->NumEle].zmax-MD->Ele[i-MD->NumEle].zmin))
00174         {
00175             if((DummyY[i]<(MD->Ele[i-MD->NumEle].zmax-MD->Ele[i-MD->NumEle].zmin)))
00176             {
00177                 DummyY[i]=1.0*((MD->Ele[i-MD->NumEle].zmax-MD->Ele[i-MD->NumEle].zmin)-DummyY[i+MD->NumEle]);
00178             }
00179             else //Only possible if DummyY UnSat = Zmax-Zmin : Check to see if this really happens in natural conditions?
00180             {
00181                 //DummyY[i]=MD->Ele[i-MD->NumEle].zmax - MD->Ele[i-MD->NumEle].zmin;
00182                 DummyY[i+MD->NumEle]=0.0;
00183             }
00184         }
00185       }
00186 
00187     /* Lateral Flux Calculation between Triangular elements Follows  */
00188     for(i=0; i<MD->NumEle; i++){
00189          for(j=0; j<3; j++){
00190               if(MD->Ele[i].nabr[j] > 0){
00191               /* Subsurface Lateral Flux Calculation between Triangular elements Follows */
00192                    if(MD->Ele[MD->Ele[i].nabr[j]-1].zmin>MD->Ele[i].zmin){
00193                         if(MD->Ele[MD->Ele[i].nabr[j]-1].zmin>MD->Ele[i].zmin+DummyY[i+2*MD->NumEle]){
00194                              Avg_Y_Sub=DummyY[MD->Ele[i].nabr[j]-1 + 2*MD->NumEle]/2;
00195                         }
00196                         else{
00197                              Avg_Y_Sub=(DummyY[i+2*MD->NumEle]+MD->Ele[i].zmin-MD->Ele[MD->Ele[i].nabr[j]-1].zmin+DummyY[MD->Ele[i].nabr[j]-1 + 2*MD->NumEle])/2;
00198                         }
00199                    }
00200                    else{
00201                         if(MD->Ele[i].zmin>MD->Ele[MD->Ele[i].nabr[j]-1].zmin+DummyY[MD->Ele[i].nabr[j]-1 + 2*MD->NumEle]){
00202                              Avg_Y_Sub=DummyY[i+2*MD->NumEle]/2;
00203                         }
00204                         else{
00205                              Avg_Y_Sub=(DummyY[i+2*MD->NumEle]+DummyY[MD->Ele[i].nabr[j]-1 + 2*MD->NumEle]+MD->Ele[MD->Ele[i].nabr[j]-1].zmin-MD->Ele[i].zmin)/2;
00206                         }
00207                    }
00208                    Dif_Y_Sub = (DummyY[i+2*MD->NumEle] + MD->Ele[i].zmin) - (DummyY[MD->Ele[i].nabr[j]-1 + 2*MD->NumEle] + MD->Ele[MD->Ele[i].nabr[j]-1].zmin);
00209                    Distance = sqrt(pow((MD->Ele[i].x - MD->Ele[MD->Ele[i].nabr[j] - 1].x), 2) + pow((MD->Ele[i].y - MD->Ele[MD->Ele[i].nabr[j] - 1].y), 2));
00210                    Avg_Ksat = (MD->Ele[i].Ksat + MD->Ele[MD->Ele[i].nabr[j] - 1].Ksat)/2.0;
00211                    Grad_Y_Sub = Dif_Y_Sub/Distance;
00212                       /* #? see the source of -1 in the nabr[j]-1 */
00213                    /* take care of macropore effect */
00214                    if((MD->Ele[i].zmax-MD->Ele[i].zmin)-DummyY[i+2*MD->NumEle]-DummyY[i+MD->NumEle]<=0){
00215                         elemSatn=1.0;
00216                    }
00217                    else{
00218                         elemSatn = DummyY[i+2*MD->NumEle]/(MD->Ele[i].zmax-MD->Ele[i].zmin);   /*  Will have to change this for other formulation */
00219                    }
00220                    if (MD->Soil[(MD->Ele[i].soil-1)].Macropore == 0){
00221                         mp_factor = 1;
00222                    }
00223                    else if (MD->Soil[(MD->Ele[i].soil-1)].Macropore == 1){
00224                         if((elemSatn>=sat_THRESH)&&(DummyY[i]>ovl_THRESH_H)){
00225                              temp1=1.0+mp_MULTFH*(elemSatn-sat_THRESH)/(1-sat_THRESH);
00226                              temp1=temp1*mpArea_CALIB+1*(1-mpArea_CALIB);
00227                         }
00228                         else{
00229                              temp1 = 1.0;
00230                         }
00231 
00232                         if((MD->Ele[MD->Ele[i].nabr[j] - 1].zmax-MD->Ele[MD->Ele[i].nabr[j] - 1].zmin)-DummyY[MD->Ele[i].nabr[j] - 1+2*MD->NumEle]-DummyY[MD->Ele[i].nabr[j] - 1+MD->NumEle]<=0){
00233                              elemSatn=1.0;
00234                         }
00235                         else{
00236                              elemSatn = DummyY[MD->Ele[i].nabr[j] - 1+2*MD->NumEle]/(MD->Ele[MD->Ele[i].nabr[j] - 1].zmax-MD->Ele[MD->Ele[i].nabr[j] - 1].zmin);   /*  Will have to change this for other formulation */
00237                         }
00238                         if((elemSatn>=sat_THRESH)&&(DummyY[MD->Ele[i].nabr[j] - 1]>ovl_THRESH_H)){
00239                              temp2=1.0+mp_MULTFH*(elemSatn-sat_THRESH)/(1-sat_THRESH);
00240                              temp2=temp2*mpArea_CALIB+1*(1-mpArea_CALIB);
00241                         }
00242                         else{
00243                              temp2 = 1.0;
00244                         }
00245 
00246                         mp_factor = (temp1 + temp2)/2.0;
00247                    }
00248 
00249                    /* groundwater flow modeled by Darcy's law */
00250                    MD->FluxSub[i][j] = mp_factor*Kh_CALIB*Avg_Ksat*Grad_Y_Sub*Avg_Y_Sub*MD->Ele[i].edge[j];
00251 
00252                    /*     Correction is being done in flux terms which can be > 0 even when there is no source water level present */
00253 
00254                    if(DummyY[i + 2*MD->NumEle] <= 0 && MD->FluxSub[i][j] > 0){
00255                         MD->FluxSub[i][j] = 0;
00256                    }
00257                    if(DummyY[MD->Ele[i].nabr[j] - 1 + 2*MD->NumEle] <= 0 && MD->FluxSub[i][j] < 0){
00258                         MD->FluxSub[i][j] = 0;
00259                    }
00260 
00261                    /* Saturation check */
00262 
00263                    if((DummyY[i + 2*MD->NumEle] >= (MD->Ele[i].zmax - MD->Ele[i].zmin)) && MD->FluxSub[i][j] < 0){
00264                         MD->FluxSub[i][j] = 0;
00265                    }
00266                    if((DummyY[MD->Ele[i].nabr[j] - 1 + 2*MD->NumEle] >= (MD->Ele[MD->Ele[i].nabr[j] - 1].zmax - MD->Ele[MD->Ele[i].nabr[j] - 1].zmin)) && MD->FluxSub[i][j] >0){
00267                         MD->FluxSub[i][j] = 0;
00268                    }
00269 
00270                    /* Surface Lateral Flux Calculation between Triangular elements Follows */
00271                    if(MD->Ele[MD->Ele[i].nabr[j] - 1].zmax>MD->Ele[i].zmax){
00272                         if(MD->Ele[MD->Ele[i].nabr[j] - 1].zmax>MD->Ele[i].zmax+DummyY[i]){
00273                              Avg_Y_Surf=DummyY[MD->Ele[i].nabr[j]-1]/2;
00274                         }
00275                         else{
00276                              Avg_Y_Surf=(DummyY[i]+MD->Ele[i].zmax-MD->Ele[MD->Ele[i].nabr[j] - 1].zmax+DummyY[MD->Ele[i].nabr[j]-1])/2;
00277                         }
00278                    }
00279                    else{
00280                         if(MD->Ele[i].zmax>MD->Ele[MD->Ele[i].nabr[j] - 1].zmax+DummyY[MD->Ele[i].nabr[j]-1]){
00281                              Avg_Y_Surf=DummyY[i]/2;
00282                         }
00283                         else{
00284                              Avg_Y_Surf=(DummyY[i]+DummyY[MD->Ele[i].nabr[j]-1]+MD->Ele[MD->Ele[i].nabr[j] - 1].zmax-MD->Ele[i].zmax)/2;
00285                         }
00286                    }
00287                    Dif_Y_Surf = (DummyY[i] + MD->Ele[i].zmax) - (DummyY[MD->Ele[i].nabr[j] - 1] + MD->Ele[MD->Ele[i].nabr[j] - 1].zmax);
00288                    Grad_Y_Surf = Dif_Y_Surf/Distance;
00289                    Avg_Sf = (MD->Ele[i].Sf + MD->Ele[MD->Ele[i].nabr[j] - 1].Sf)/2.0;
00290                    Avg_Rough = 0.5*(MD->Ele[i].Rough + MD->Ele[MD->Ele[i].nabr[j] - 1].Rough);
00291                    CrossA = Avg_Y_Surf*MD->Ele[i].edge[j];
00292 
00293                    OverlandFlow(MD->FluxSurf,i,j,MD->SurfMode, Avg_Y_Surf,Grad_Y_Surf,Avg_Sf,Alfa,Beta,CrossA,Avg_Rough,1,1);
00294 
00295                     if(isnan(MD->FluxSurf[i][j])==1){
00296                         printf("\n1: %f %d %d %lf %lf %lf",t,MD->Ele[i].index,MD->Ele[MD->Ele[i].nabr[j]-1].index,DummyY[i],DummyY[MD->Ele[i].nabr[j]-1],MD->FluxSurf[i][j]);
00297                         getchar();
00298                    }
00299                    /*     Correction is being done in flux terms which can be > 0 even when there is no source water level present */
00300 
00301                    if(DummyY[i] <= 0 && MD->FluxSurf[i][j] > 0){
00302                         MD->FluxSurf[i][j] = 0;
00303                    }
00304                    if(DummyY[MD->Ele[i].nabr[j] - 1] <= 0 && MD->FluxSurf[i][j] < 0){
00305                         MD->FluxSurf[i][j] = 0;
00306                    }
00307               }
00308               else{
00309                    /* Handle boundary conditions for elements. No flow (natural) boundary condition is default */
00310                    if(MD->Ele[i].BC == 0){
00311                         MD->FluxSurf[i][j] = 0;
00312                         MD->FluxSub[i][j] = 0;
00313                    }
00314                    else{
00315                         //MD->FluxSurf[i][j] = 0;
00316 
00317                         /* this part of code has not been tested! */
00318                         if(MD->Ele[i].BC > 0){         // Dirichlet BC
00319                 loc_bcEle = Interpolation(&MD->TSD_EleBC[(MD->Ele[i].BC)-1], t);
00320 
00321                                 if( loc_bcEle < MD->Ele[i].zmin ){
00322                                         Avg_Y_Surf = (DummyY[i])/2;
00323                                         Avg_Y_Sub  = (DummyY[i+2*MD->NumEle])/2;
00324                                 }
00325                                 else if(loc_bcEle < MD->Ele[i].zmax){
00326                                         Avg_Y_Surf = (DummyY[i])/2;
00327                                         Avg_Y_Sub  = (loc_bcEle - MD->Ele[i].zmin + DummyY[i+2*MD->NumEle])/2.0;
00328                                 }
00329                                 else{
00330                                         Avg_Y_Surf = (DummyY[i] + loc_bcEle - MD->Ele[i].zmax)/2;
00331                                         Avg_Y_Sub  = MD->Ele[i].zmax - MD->Ele[i].zmin;
00332                                 }
00333                                 Distance = sqrt(pow(MD->Ele[i].edge[0]*MD->Ele[i].edge[1]*MD->Ele[i].edge[2]/(4*MD->Ele[i].area), 2) - pow(MD->Ele[i].edge[j]/2, 2));
00334 
00335                                 //Sub
00336                                 Dif_Y_Sub = DummyY[i+2*MD->NumEle] + MD->Ele[i].zmin - loc_bcEle;
00337                                 Avg_Ksat = MD->Ele[i].Ksat;
00338                                 Grad_Y_Sub = Dif_Y_Sub/Distance;
00339                                 MD->FluxSub[i][j] = Avg_Ksat * Grad_Y_Sub * Avg_Y_Sub * MD->Ele[i].edge[j];
00340 
00341                                 //Surf
00342                                 Dif_Y_Surf = DummyY[i] + MD->Ele[i].zmax - loc_bcEle;
00343                                 Grad_Y_Surf = Dif_Y_Surf / Distance;
00344                                 MD->FluxSurf[i][j] = 0.1*(Grad_Y_Surf>0?1:-1)* (Avg_Y_Sub * MD->Ele[i].edge[j] ) * (pow(pow(Avg_Y_Surf, 1.0/3.0),2)/(MD->Ele[i].Rough)) * sqrt((Grad_Y_Surf>0?1:-1)*Grad_Y_Surf);
00345 
00346                                 MD->FluxSurf[i][j] = sqrt(DummyY[i]/Distance)*pow(DummyY[i],2.0/3.0)*(DummyY[i]/2)*MD->Ele[i].edge[j]/MD->Ele[i].Rough;
00347                         }
00348 
00349                                  /* Note the distance calculated here is the distance between circumcenter of the triangle and the edge on which BDD. condition is defined
00350                              Distance = sqrt(pow(MD->Ele[i].edge[0]*MD->Ele[i].edge[1]*MD->Ele[i].edge[2]/(4*MD->Ele[i].area), 2) - pow(MD->Ele[i].edge[j]/2, 2));
00351                              Avg_Ksat = MD->Ele[i].Ksat;
00352                              Grad_Y_Sub = Dif_Y_Sub/Distance;
00353 
00354                              MD->FluxSub[i][j] = Avg_Ksat*Grad_Y_Sub*Avg_Y_Sub*MD->Ele[i].edge[j];
00355                         } */
00356 
00357 
00358                         else{                          // Nuemann BC
00359                              MD->FluxSub[i][j] = Interpolation(&MD->TSD_EleBC[(-MD->Ele[i].BC)-1+MD->Num1BC], t);
00360                         }
00361             /*
00362             Distance = sqrt(pow(MD->Ele[i].edge[0]*MD->Ele[i].edge[1]*MD->Ele[i].edge[2]/(4*MD->Ele[i].area), 2) - pow(MD->Ele[i].edge[j]/2, 2));
00363                         MD->FluxSub[i][j]=MD->Ele[i].Ksat*(DummyY[i+2*MD->NumEle]/Distance)*(DummyY[i+2*MD->NumEle]/2)*MD->Ele[i].edge[j];
00364             */
00365                         /* zero-depth-gradient boundary conditions */
00366             /*
00367                             Dif_Y_Surf=DummyY[i]/Distance;
00368                         Avg_Rough = MD->Ele[i].Rough;
00369                         Avg_Y_Riv = DummyY[i]/2;
00370                         MD->FluxSurf[i][j] = sqrt(Dif_Y_Surf)*pow(DummyY[i],2.0/3.0)*Avg_Y_Riv*MD->Ele[i].edge[j]/Avg_Rough;
00371             */
00372                    }
00373 
00374                    /*     Correction is being done in flux terms which can be > 0 even when there is no source water level present */
00375                    if(DummyY[i + 2*MD->NumEle] <= 0 && MD->FluxSub[i][j] > 0){
00376                         MD->FluxSub[i][j] = 0;
00377                    }
00378               }
00379          }
00380          /************************************************************************************************/
00381 
00382          /**************************************Evaporation from ground*********************************************************/
00383          if((MD->Ele[i].zmax-MD->Ele[i].zmin)-DummyY[i+2*MD->NumEle]-DummyY[i+MD->NumEle]<=0){
00384               elemSatn=1.0;
00385             }
00386          else{
00387               //elemSatn = 0.5*(1-cos(3.14*(DummyY[i+MD->NumEle]/(MD->Ele[i].zmax-MD->Ele[i].zmin-DummyY[i+2*MD->NumEle]))));    /*  Will have to change this for other formulation */
00388               elemSatn = (DummyY[i+2*MD->NumEle]+DummyY[i+MD->NumEle]>(MD->Ele[i].zmax-MD->Ele[i].zmin)-MD->Ele[i].RzD-rzd_CALIB)? 0.5*(1-cos(3.14*(DummyY[i+2*MD->NumEle]/(MD->Ele[i].zmax-MD->Ele[i].zmin)))):0;
00389          }
00390          Rn = Interpolation(&MD->TSD_Rn[MD->Ele[i].Rn-1], t);
00391          //G = Interpolation(&MD->TSD_G[MD->Ele[i].G-1], t);
00392          T = Interpolation(&MD->TSD_Temp[MD->Ele[i].temp-1], t);
00393          Vel = Interpolation(&MD->TSD_WindVel[MD->Ele[i].WindVel-1], t);
00394          RH = Interpolation(&MD->TSD_Humidity[MD->Ele[i].humidity-1], t);
00395          VP = Interpolation(&MD->TSD_Pressure[MD->Ele[i].pressure-1], t);
00396          P = 101.325*pow(10,3)*pow((293-0.0065*MD->Ele[i].zmax)/293,5.26);
00397          Delta = 2503*pow(10,3)*exp(17.27*T/(T+237.3))/(pow(237.3 + T, 2));
00398          Gamma = P*1.0035*0.92/(0.622*2441);
00399          LAI = Interpolation(&MD->TSD_LAI[MD->Ele[i].LC-1], t);
00400          zero_dh=Interpolation(&MD->TSD_DH[MD->Ele[i].LC-1], t);
00401          cnpy_h = zero_dh/(1.1*(0.0000001+log(1+pow(0.007*LAI,0.25))));
00402 
00403          /*if(LAI<2.85)
00404         {
00405         rl= 0.0002 + 0.3*cnpy_h*pow(0.07*LAI,0.5);
00406         }
00407         else
00408         {
00409         rl= 0.3*cnpy_h*(1-(zero_dh/cnpy_h));
00410         }
00411         */
00412         rl=Interpolation(&MD->TSD_DH[MD->Ele[i].LC-1], t);
00413         r_a = log(MD->Ele[i].windH/rl)*log(10*MD->Ele[i].windH/rl)/(Vel*0.16);
00414 
00415  /* $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ DELETE 0.01 BELOW $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ */
00416 
00417         MD->EleET[i][2] = (1-MD->Ele[i].VegFrac)*(Rn*(1-MD->Ele[i].Albedo)*Delta+(1.2*1003.5*((VP/RH)-VP)/r_a))/(1000.0*2441000.0*(Delta+Gamma));
00418         MD->EleET[i][2]=et2_CALIB*MD->EleET[i][2];
00419 
00420         if(LAI>0.0){
00421              Rmax = 5000.0/(24*3600);        /* Unit day_per_m */
00422              f_r= 1.1*Rn*(1-exp(-LAI))/(MD->Ele[i].Rs_ref*LAI);
00423              alpha_r= (1+f_r)/(1+(MD->Ele[i].Rmin/Rmax));
00424              eta_s= 1- 0.0016*(pow((24.85-T),2));
00425              beta_s= elemSatn>EPSILON/1000.0?elemSatn:EPSILON/1000.0;
00426 
00427              r_s=(MD->Ele[i].Rmin*alpha_r/(beta_s*LAI*pow(eta_s,4)))> Rmax?Rmax:(MD->Ele[i].Rmin*alpha_r/(beta_s*LAI*pow(eta_s,4)));
00428 
00429              MD->EleET[i][1] = (LAI/MD->Ele[i].LAImax)*MD->Ele[i].VegFrac*(1-pow((MD->EleIS[i]<0?0:MD->EleIS[i])/MD->EleISmax[i],2.0/3))*(Rn*(1-MD->Ele[i].Albedo)*Delta+(1.2*1003.5*((VP/RH)-VP)/r_a))/(1000*2441000.0*(Delta+Gamma*(1+r_s/r_a)));
00430         }
00431         else{
00432              MD->EleET[i][1] =0.0;
00433         }
00434 
00435 
00436         /**********************************************************************************************************************/
00437 
00438         MD->EleET[i][1]=et1_CALIB*MD->EleET[i][1];
00439 
00440 
00441         AquiferDepth = MD->Ele[i].zmax - MD->Ele[i].zmin;
00442 
00443         if(DummyY[i+2*MD->NumEle]+DummyY[i+MD->NumEle] >= AquiferDepth){
00444              Deficit = 0;
00445              MD->EleVic[i] = 0;
00446         }
00447         else if(DummyY[i+2*MD->NumEle]+DummyY[i+MD->NumEle]>AquiferDepth-MD->Ele[i].RzD){
00448                Deficit = AquiferDepth - DummyY[i+2*MD->NumEle];
00449                MD->EleVic[i] = MD->Ele[i].Ksat*(1+(DummyY[i]/MD->Ele[i].RzD));
00450         }
00451         else{
00452              Deficit = AquiferDepth - DummyY[i+2*MD->NumEle];
00453              MD->EleVic[i] = MD->Ele[i].Ksat*(1+(DummyY[i]-log((DummyY[i+MD->NumEle]+EPSILON/10000.0)/(Deficit-MD->Ele[i].RzD+EPSILON/10000.0))/MD->Ele[i].Alpha)/(MD->Ele[i].RzD));/* Interpolation(&MD->TSD_Inc[MD->Soil[(MD->Ele[i].soil-1)].Inf-1], t);*/    /* ## Take care of this using saturation rate instead of */
00454         }
00455 
00456         MD->EleVic[i]=MD->EleVic[i]/Vic_CALIB;
00457 
00458         if(DummyY[i+MD->NumEle] < Deficit){                                     /*  ## redo this condition with only phi and z */
00459              if(DummyY[i] > 0){
00460              /***************************Addn***************************/
00461                   if((MD->Ele[i].zmax-MD->Ele[i].zmin)-DummyY[i+2*MD->NumEle]-DummyY[i+MD->NumEle]<=0){
00462                        elemSatn=1.0;
00463                   }
00464                   else{
00465                        elemSatn = DummyY[i+2*MD->NumEle]/(MD->Ele[i].zmax-MD->Ele[i].zmin);   /*  Will have to change this for other formulation */
00466                   }
00467 
00468                   if (MD->Soil[(MD->Ele[i].soil-1)].Macropore == 0){
00469                        mp_factor = 1.0;
00470                   }
00471                   else if (MD->Soil[(MD->Ele[i].soil-1)].Macropore == 1){
00472                        if((elemSatn>=sat_THRESH)&&(DummyY[i]>ovl_THRESH_V)){
00473                             mp_factor=1.0+mp_MULTFV*(elemSatn-sat_THRESH)/(1-sat_THRESH);
00474                             mp_factor=mp_factor*mpArea_CALIB+1*(1-mpArea_CALIB);
00475                        }
00476                        else{
00477                             mp_factor = 1.0;
00478                        }
00479                   }
00480 
00481                   /*********************************************************/
00482                   if((MD->Ele[i].zmax-MD->Ele[i].zmin)-DummyY[i+2*MD->NumEle]-DummyY[i+MD->NumEle]<=0){
00483                        elemSatn=1.0;
00484                   }
00485                   else{
00486                        elemSatn = (DummyY[i+2*MD->NumEle]+DummyY[i+MD->NumEle]>(MD->Ele[i].zmax-MD->Ele[i].zmin)-MD->Ele[i].RzD-rzd_CALIB)? 0.5*(1-cos(3.14*(DummyY[i+2*MD->NumEle]/(MD->Ele[i].zmax-MD->Ele[i].zmin)))):0;
00487                   }
00488 
00489                   DummyDY[i] = MD->EleNetPrep[i] - mp_factor*MD->EleVic[i] - MD->EleET[i][2];
00490                   DummyDY[i+MD->NumEle] = (1.0-mpArea_CALIB)*MD->EleVic[i];                    /* ##Clean DY variables to some dummy for defining right hand side */
00491                   DummyDY[i+2*MD->NumEle] = DummyDY[i+2*MD->NumEle]+(mp_factor-1*(1-mpArea_CALIB))*MD->EleVic[i];
00492              }
00493              else if((MD->EleNetPrep[i] > MD->EleVic[i]+MD->EleET[i][2])){
00494                   DummyDY[i] = MD->EleNetPrep[i] - (MD->EleVic[i]+MD->EleET[i][2]);
00495                   DummyDY[i+MD->NumEle] = MD->EleVic[i];
00496              }
00497              else if((MD->EleNetPrep[i] < MD->EleVic[i]+MD->EleET[i][2])){       /* ##Think abt putting <= here */
00498 
00499                   if(MD->EleNetPrep[i] < MD->EleVic[i]){
00500                        DummyDY[i] = 0;
00501                          DummyDY[i+MD->NumEle] = MD->EleNetPrep[i]-elemSatn*MD->EleET[i][2];
00502                 //BHATT
00503                 MD->EleET[i][2]=elemSatn*MD->EleET[i][2];
00504                   }
00505                   else{
00506                          DummyDY[i] = 0;
00507                          DummyDY[i+MD->NumEle] = MD->EleVic[i];
00508                   }
00509              }
00510              else{
00511                   DummyDY[i] = 0;
00512                   DummyDY[i+MD->NumEle] = MD->EleNetPrep[i]-elemSatn*MD->EleET[i][2];
00513             //BHATT
00514             MD->EleET[i][2]=elemSatn*MD->EleET[i][2];
00515              }
00516         }
00517         else{                                    /*  ## redo this condition */
00518 
00519              /* The reason of this happening is not clear, therefore the treatment is not secure
00520              Fortranately, this will not happen in most cases. */
00521              if(DummyY[i]>0){
00522                   DummyDY[i] = MD->EleNetPrep[i]-MD->EleET[i][2];
00523                    DummyDY[i+2*MD->NumEle] = 0;
00524              }
00525               else{
00526                   if(MD->EleNetPrep[i]>0){
00527                        if(MD->EleNetPrep[i]>MD->EleET[i][2]){
00528                             DummyDY[i] =MD->EleNetPrep[i]-MD->EleET[i][2];
00529                          }
00530                           else{
00531                                DummyDY[i] = 0;
00532                           }
00533                           DummyDY[i+2*MD->NumEle] = 0;
00534                   }
00535                    else{
00536                         DummyDY[i] = 0;
00537                         DummyDY[i+2*MD->NumEle] = 0-elemSatn*MD->EleET[i][2];
00538                     //BHATT
00539                     MD->EleET[i][2]=elemSatn*MD->EleET[i][2];
00540                   }
00541              }
00542         }
00543 
00544         if(DummyY[i+2*MD->NumEle]>=AquiferDepth-MD->Ele[i].RzD){
00545              DummyDY[i+2*MD->NumEle]=DummyDY[i+2*MD->NumEle] - MD->EleET[i][1];
00546         }
00547         else{
00548               DummyDY[i+MD->NumEle]=DummyDY[i+MD->NumEle] - MD->EleET[i][1];
00549         }
00550 
00551         /*for(j=0; j<3; j++)
00552         {
00553         DummyDY[i] =  DummyDY[i] - MD->FluxSurf[i][j]/MD->Ele[i].area;
00554         }
00555 
00556         for(j=0; j<3; j++)
00557         {
00558          DummyDY[i+2*MD->NumEle] = DummyDY[i+2*MD->NumEle] - MD->FluxSub[i][j]/MD->Ele[i].area;
00559         } */
00560          /************************************************************************************************/
00561     }
00562 
00563 
00564     /* initialize river flux */
00565     for(i=0; i<MD->NumRiv; i++){
00566          for(j=0; j<6; j++){
00567               MD->FluxRiv[i][j] = 0;
00568          }
00569     }
00570 
00571     /* Lateral Flux Calculation between River-River and River-Triangular elements Follows */
00572     for(i=0; i<MD->NumRiv; i++){
00573 
00574          /* Note: the ordering of river segment in input file has to be: from UP to DOWN */
00575          TotalY_Riv = DummyY[i + 3*MD->NumEle] + MD->Riv[i].zmin;
00576          Perem = CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[i].shape - 1].interpOrd,DummyY[i + 3*MD->NumEle],MD->Riv_Shape[MD->Riv[i].shape - 1].coeff,2);
00577          /*    if(DummyY[10 + 3*MD->NumEle]>0)
00578          {
00579                printf("\n%lf %e %e %e %e area",t,DummyY[10 + 3*MD->NumEle],MD->Riv_Shape[MD->Riv[10].shape - 1].coeff,Perem,CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[10].shape - 1].interpOrd,DummyY[10 + 3*MD->NumEle],MD->Riv_Shape[MD->Riv[10].shape - 1].coeff,1));
00580                getchar();
00581          }
00582          */
00583          /* Lateral Flux Calculation between River-River element Follows */
00584          if(MD->Riv[i].down > 0){
00585               TotalY_Riv_down = DummyY[MD->Riv[i].down - 1 + 3*MD->NumEle] + MD->Riv[MD->Riv[i].down - 1].zmin;
00586               Perem_down = 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,2);
00587               Avg_Perem = (Perem + Perem_down)/2.0;    /* Avg perimeter */
00588               if(MD->Riv[MD->Riv[i].down - 1].zmin>MD->Riv[i].zmin){
00589                    if(MD->Riv[MD->Riv[i].down - 1].zmin>MD->Riv[i].zmin+DummyY[i + 3*MD->NumEle]){
00590                         Avg_Y_Riv=DummyY[MD->Riv[i].down - 1 + 3*MD->NumEle]/2;
00591                    }
00592                    else{
00593                         Avg_Y_Riv=(DummyY[i + 3*MD->NumEle]+MD->Riv[i].zmin-MD->Riv[MD->Riv[i].down - 1].zmin+DummyY[MD->Riv[i].down - 1 + 3*MD->NumEle])/2;
00594                    }
00595               }
00596               else{
00597                    if(MD->Riv[i].zmin>MD->Riv[MD->Riv[i].down - 1].zmin+DummyY[MD->Riv[i].down - 1 + 3*MD->NumEle]){
00598                         Avg_Y_Riv=DummyY[i + 3*MD->NumEle]/2;
00599                    }
00600                    else{
00601                         Avg_Y_Riv=(DummyY[i + 3*MD->NumEle]+DummyY[MD->Riv[i].down - 1 + 3*MD->NumEle]+MD->Riv[MD->Riv[i].down - 1].zmin-MD->Riv[i].zmin)/2;
00602                    }
00603               }
00604              Avg_Rough = (MD->Riv_Mat[MD->Riv[i].material - 1].Rough + MD->Riv_Mat[MD->Riv[MD->Riv[i].down - 1].material-1].Rough)/2.0;
00605              //Distance = sqrt(pow(MD->Riv[i].x - MD->Riv[MD->Riv[i].down - 1].x, 2) + pow(MD->Riv[i].y - MD->Riv[MD->Riv[i].down - 1].y, 2));
00606              Distance = (MD->Riv[i].Length+MD->Riv[MD->Riv[i].down - 1].Length)/2;
00607 
00608              Dif_Y_Riv = (TotalY_Riv - TotalY_Riv_down)/Distance;
00609              Avg_Sf = (MD->Riv_Mat[MD->Riv[i].material - 1].Sf + MD->Riv_Mat[MD->Riv[MD->Riv[i].down - 1].material-1].Sf)/2.0;
00610              /*CrossA = 0.5*(CS_AreaOrPerem(MD->Riv_Shape[MD->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));
00611              */
00612              CrossA = CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[i].shape - 1].interpOrd,Avg_Y_Riv,MD->Riv_Shape[MD->Riv[i].shape - 1].coeff,1);
00613              OverlandFlow(MD->FluxRiv,i,1,MD->RivMode, Avg_Y_Riv,Dif_Y_Riv,Avg_Sf,Alfa,Beta,CrossA,Avg_Rough,0,Avg_Perem);
00614 
00615              /* Correction is being done in flux terms which can be > 0 even when there is no source water level present */
00616              if(DummyY[i + 3*MD->NumEle] <= 0 && MD->FluxRiv[i][1] > 0){
00617                   MD->FluxRiv[i][1] = 0;
00618              }
00619              else if(DummyY[MD->Riv[i].down - 1 + 3*MD->NumEle] <= 0 && MD->FluxRiv[i][1] < 0){
00620                   MD->FluxRiv[i][1] = 0;
00621              }
00622 
00623              /* accumulate to get in-flow for down segments: [0] for inflow, [1] for outflow */
00624              MD->FluxRiv[MD->Riv[i].down - 1][0] = MD->FluxRiv[MD->Riv[i].down - 1][0] + MD->FluxRiv[i][1];
00625          }
00626          /* Here the if statement ends for if there is a downstream river segment or not */
00627          else{
00628               switch(MD->Riv[i].down){
00629                    case -1:
00630 
00631                         /* Dirichlet boundary condition */
00632                         TotalY_Riv_down = Interpolation(&MD->TSD_Riv[(MD->Riv[i].BC)-1], t) + MD->Node[MD->Riv[i].ToNode-1].zmin + MD->Riv_Shape[MD->Riv[i].shape-1].bed;
00633                         Distance = (MD->Riv[i].Length)*0.5;
00634                         Dif_Y_Riv = (TotalY_Riv - TotalY_Riv_down)/Distance;
00635                         Avg_Sf = MD->Riv_Mat[MD->Riv[i].material - 1].Sf;
00636                         Avg_Rough = MD->Riv_Mat[MD->Riv[i].material-1].Rough;
00637                         Avg_Y_Riv = DummyY[i + 3*MD->NumEle];
00638                         Avg_Perem = Perem;
00639                         CrossA = CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[i].shape - 1].interpOrd,DummyY[i + 3*MD->NumEle],MD->Riv_Shape[MD->Riv[i].shape - 1].coeff,1);
00640 
00641                         OverlandFlow(MD->FluxRiv,i,1,MD->RivMode, Avg_Y_Riv,Dif_Y_Riv,Avg_Sf,Alfa,Beta,CrossA,Avg_Rough,0,Avg_Perem);
00642 
00643                         break;
00644 
00645                    case -2:
00646 
00647                         /* Neumann boundary condition */
00648                         MD->FluxRiv[i][1] = Interpolation(&MD->TSD_Riv[MD->Riv[i].BC-1], t);
00649                         break;
00650 
00651                    case -3:
00652 
00653                         /* zero-depth-gradient boundary conditions */
00654                         Distance = (MD->Riv[i].Length)*0.5;
00655                         //Dif_Y_Riv = (MD->Riv[i].zmin - (MD->Node[MD->Riv[i].ToNode-1].zmax -MD->Riv[i].depth+ MD->Riv_Shape[MD->Riv[i].shape-1].bed))/Distance;
00656                           Dif_Y_Riv=0.1/Distance;
00657                         Avg_Rough = MD->Riv_Mat[MD->Riv[i].material-1].Rough;
00658                         Avg_Y_Riv = DummyY[i + 3*MD->NumEle];
00659                         Avg_Perem = Perem;
00660                         CrossA = CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[i].shape - 1].interpOrd,DummyY[i + 3*MD->NumEle],MD->Riv_Shape[MD->Riv[i].shape - 1].coeff,1);
00661                         MD->FluxRiv[i][1] = sqrt(Dif_Y_Riv)*CrossA*(Perem>0?pow(CrossA/Perem,2.0/3.0):0)/Avg_Rough;
00662                         break;
00663                         /* #? How is critical dept being defined */
00664                    case -4:
00665 
00666                         /* Critical Depth boundary conditions */
00667                         CrossA = CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[i].shape - 1].interpOrd,DummyY[i + 3*MD->NumEle],MD->Riv_Shape[MD->Riv[i].shape - 1].coeff,1);
00668                         MD->FluxRiv[i][1] = CrossA*sqrt(GRAV*DummyY[i + 3*MD->NumEle]);
00669                         break;
00670 
00671                    default:
00672                         printf("Fatal Error: River Routing Boundary Condition Type Is Wrong!");
00673                         exit(1);
00674               }
00675 
00676               /*     Correction is being done in flux terms which can be > 0 even when there is no source water level present */
00677               if(DummyY[i + 3*MD->NumEle] <= 0 && MD->FluxRiv[i][1] > 0){
00678                    MD->FluxRiv[i][1] = 0;
00679               }
00680               /*      else if(DummyY[MD->Riv[i].down - 1 + 2*MD->NumEle] <= 0 && MD->FluxRiv[i][1] < 0)
00681               {
00682                       MD->FluxRiv[i][1] = 0;
00683                        }  */            /* ######## Does this need to be commented */
00684 
00685               /* #? Have to take care of this before it can be implemented on ny basin other than ShaleHIlls */
00686               /* need some work to take care of multiple output Q */
00687               MD->Q = MD->FluxRiv[i][1];
00688 
00689          }
00690 
00691          /* Lateral Surface Flux Calculation between River-Triangular element Follows */
00692          if (MD->Riv[i].LeftEle > 0){
00693               OLflowFromEleToRiv(DummyY[MD->Riv[i].LeftEle - 1],MD->Ele[MD->Riv[i].LeftEle - 1].zmax,MD->Riv[i].x,MD->Ele[MD->Riv[i].LeftEle - 1].x,MD->Riv[i].y,MD->Ele[MD->Riv[i].LeftEle - 1].y,MD->Riv_Mat[MD->Riv[i].material-1].Cwr, MD->Riv[i].zmax,TotalY_Riv,MD->FluxRiv,i,2,MD->Riv[i].Length);
00694 
00695               /*     Correction is being done in flux terms which can be > 0 even when there is no source water level present */
00696               if(DummyY[i + 3*MD->NumEle] <= 0 && MD->FluxRiv[i][2] > 0){
00697                    MD->FluxRiv[i][2] = 0;
00698               }
00699 
00700               if(DummyY[MD->Riv[i].LeftEle - 1] <= 0 && MD->FluxRiv[i][2] < 0){
00701                    MD->FluxRiv[i][2] = 0;
00702               }
00703 
00704               /* replace overland flux item */
00705               for(j=0; j < 3; j++){
00706                    /* this may cause trouble when there are 2 boundary side and one is not neutral.*/
00707                    if(MD->Ele[MD->Riv[i].LeftEle - 1].nabr[j] == MD->Riv[i].RightEle){
00708                         MD->FluxSurf[MD->Riv[i].LeftEle - 1][j] = -MD->FluxRiv[i][2];
00709                    }
00710               }
00711          }
00712 
00713          if (MD->Riv[i].RightEle > 0){
00714               OLflowFromEleToRiv(DummyY[MD->Riv[i].RightEle - 1],MD->Ele[MD->Riv[i].RightEle - 1].zmax,MD->Riv[i].x,MD->Ele[MD->Riv[i].RightEle - 1].x,MD->Riv[i].y,MD->Ele[MD->Riv[i].RightEle - 1].y,MD->Riv_Mat[MD->Riv[i].material-1].Cwr, MD->Riv[i].zmax,TotalY_Riv,MD->FluxRiv,i,3,MD->Riv[i].Length);
00715 
00716               /*     Correction is being done in flux terms which can be > 0 even when there is no source water level present */
00717               if(DummyY[i + 3*MD->NumEle] <= 0 && MD->FluxRiv[i][3] > 0){
00718                    MD->FluxRiv[i][3] = 0;
00719               }
00720 
00721               if(DummyY[MD->Riv[i].RightEle - 1] <= 0 && MD->FluxRiv[i][3] < 0){
00722                    MD->FluxRiv[i][3] = 0;
00723               }
00724 
00725               /* replace overland flux item */
00726               for(j=0; j < 3; j++){
00727                    if(MD->Ele[MD->Riv[i].RightEle - 1].nabr[j] == MD->Riv[i].LeftEle){
00728                         MD->FluxSurf[MD->Riv[i].RightEle - 1][j] = -MD->FluxRiv[i][3];
00729                         break;
00730                    }
00731               }
00732          }
00733          /* Lateral Sub-surface Flux Calculation between River-Triangular element Follows */
00734          if (MD->Riv[i].LeftEle > 0){
00735               /* Left Neighbor Groundwater Head */
00736               /*********************TRANSFER THIS INSIDE FUNC*****************************/
00737               if((MD->Ele[MD->Riv[i].LeftEle - 1].zmax-MD->Ele[MD->Riv[i].LeftEle - 1].zmin)-DummyY[MD->Riv[i].LeftEle - 1+2*MD->NumEle]-DummyY[MD->Riv[i].LeftEle - 1+MD->NumEle]<=0){
00738                    elemSatn=1.0;
00739               }
00740               else{
00741                    elemSatn = DummyY[MD->Riv[i].LeftEle - 1+2*MD->NumEle]/(MD->Ele[MD->Riv[i].LeftEle - 1].zmax-MD->Ele[MD->Riv[i].LeftEle - 1].zmin);   /*  Will have to change this for other formulation */
00742               }
00743 
00744               if (MD->Soil[(MD->Ele[MD->Riv[i].LeftEle - 1].soil-1)].Macropore == 0){
00745                    mp_factor = 1;
00746               }
00747               else if (MD->Soil[(MD->Ele[MD->Riv[i].LeftEle - 1].soil-1)].Macropore == 1){
00748                    if((elemSatn>=sat_THRESH)&&(DummyY[MD->Riv[i].LeftEle - 1]>ovl_THRESH_H)){
00749                         mp_factor=1.0+mp_MULTFH*(elemSatn-sat_THRESH)/(1-sat_THRESH);
00750                         mp_factor=mp_factor*mpArea_CALIB+1*(1-mpArea_CALIB);
00751                    }
00752                    else{
00753                         mp_factor = 1.0;
00754                    }
00755               }
00756               /**************************************************************************/
00757 
00758 
00759               if( (MD->Riv[i].zmin + DummyY[i+3*MD->NumEle] - (DummyY[MD->Riv[i].LeftEle - 1 + 2*MD->NumEle]+MD->Ele[MD->Riv[i].LeftEle - 1].zmin))>0 ){
00760                       loc_perem=Perem;
00761               }
00762               else{
00763                     if(MD->Ele[MD->Riv[i].LeftEle - 1].zmin < MD->Riv[i].zmin){
00764                           if( (DummyY[MD->Riv[i].LeftEle - 1 + 2*MD->NumEle] + MD->Ele[MD->Riv[i].LeftEle - 1].zmin - MD->Riv[i].zmin) > MD->Riv[i].depth ){
00765                               loc_perem = CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[i].shape - 1].interpOrd,MD->Riv[i].depth,MD->Riv_Shape[MD->Riv[i].shape - 1].coeff,2);
00766 
00767                         }
00768                           else{
00769                               loc_perem = CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[i].shape - 1].interpOrd,(DummyY[MD->Riv[i].LeftEle - 1 + 2*MD->NumEle] + MD->Ele[MD->Riv[i].LeftEle - 1].zmin - MD->Riv[i].zmin),MD->Riv_Shape[MD->Riv[i].shape - 1].coeff,2);
00770                         }
00771                     }
00772                     else{
00773                         if( DummyY[MD->Riv[i].LeftEle - 1 + 2*MD->NumEle] > MD->Riv[i].depth){
00774                             loc_perem = CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[i].shape - 1].interpOrd,MD->Riv[i].depth,MD->Riv_Shape[MD->Riv[i].shape - 1].coeff,2);
00775                         }
00776                         else{
00777                             loc_perem = CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[i].shape - 1].interpOrd,DummyY[MD->Riv[i].LeftEle - 1 + 2*MD->NumEle],MD->Riv_Shape[MD->Riv[i].shape - 1].coeff,2);
00778                         }
00779                     }
00780               }
00781 
00782               GWflowFromEleToRiv(DummyY[MD->Riv[i].LeftEle - 1 + 2*MD->NumEle],MD->Ele[MD->Riv[i].LeftEle - 1].zmax,MD->Ele[MD->Riv[i].LeftEle - 1].zmin,MD->Riv[i].x,MD->Ele[MD->Riv[i].LeftEle - 1].x,MD->Riv[i].y,MD->Ele[MD->Riv[i].LeftEle - 1].y,MD->Soil[(MD->Ele[MD->Riv[i].LeftEle - 1].soil-1)].Macropore,DummyY[i+3*MD->NumEle],TotalY_Riv,MD->FluxRiv,i,4,MD->Riv[i].Length,MD->Soil[(MD->Ele[MD->Riv[i].LeftEle - 1].soil-1)].base,mp_factor,loc_perem,MD->Ele[MD->Riv[i].LeftEle - 1].Ksat,MD->Ele[MD->Riv[i].LeftEle-1].RzD); /* delete replace Wid by 0.5*avg_perim */
00783 
00784               /* Saturation check */
00785               if((DummyY[MD->Riv[i].LeftEle - 1 + 2*MD->NumEle] >= MD->Ele[MD->Riv[i].LeftEle - 1 + 2*MD->NumEle].zmax-MD->Ele[MD->Riv[i].LeftEle - 1 + 2*MD->NumEle].zmin) && MD->FluxRiv[i][4] > 0){
00786                    MD->FluxRiv[i][4]  = 0;
00787               }
00788 
00789 
00790               /*     Correction is being done in flux terms which can be > 0 even when there is no source water level present */
00791               if(DummyY[i+3*MD->NumEle] <= 0 && MD->FluxRiv[i][4] > 0){
00792                    MD->FluxRiv[i][4] = 0;
00793               }
00794               if(DummyY[MD->Riv[i].LeftEle - 1 + 2*MD->NumEle] <= 0 && MD->FluxRiv[i][4] < 0){
00795                    MD->FluxRiv[i][4] = 0;
00796               }
00797 
00798               /* modify groundwater flux item */
00799               for(j=0; j < 3; j++){
00800                    if(MD->Ele[MD->Riv[i].LeftEle - 1].nabr[j] == MD->Riv[i].RightEle){
00801                         if(MD->Ele[MD->Ele[MD->Riv[i].LeftEle - 1].nabr[j]-1].zmin>MD->Ele[MD->Riv[i].LeftEle - 1].zmin){
00802                                 if(MD->Ele[MD->Ele[MD->Riv[i].LeftEle - 1].nabr[j]-1].zmin>MD->Ele[MD->Riv[i].LeftEle - 1].zmin+DummyY[MD->Riv[i].LeftEle - 1+2*MD->NumEle]){
00803                                    Avg_Y_Sub=DummyY[MD->Ele[MD->Riv[i].LeftEle - 1].nabr[j]-1 + 2*MD->NumEle]/2;
00804                              }
00805                                 else{
00806                                       Avg_Y_Sub=(DummyY[MD->Riv[i].LeftEle - 1+2*MD->NumEle]+MD->Ele[MD->Riv[i].LeftEle - 1].zmin-MD->Ele[MD->Ele[MD->Riv[i].LeftEle - 1].nabr[j]-1].zmin+DummyY[MD->Ele[MD->Riv[i].LeftEle - 1].nabr[j]-1 + 2*MD->NumEle])/2;
00807                              }
00808                         }
00809                            else{
00810                                 if(MD->Ele[MD->Riv[i].LeftEle - 1].zmin>MD->Ele[MD->Ele[MD->Riv[i].LeftEle - 1].nabr[j]-1].zmin+DummyY[MD->Ele[MD->Riv[i].LeftEle - 1].nabr[j]-1 + 2*MD->NumEle]){
00811                                    Avg_Y_Sub=DummyY[MD->Riv[i].LeftEle - 1+2*MD->NumEle]/2;
00812                                 }
00813                                 else{
00814                                       Avg_Y_Sub=(DummyY[MD->Riv[i].LeftEle - 1+2*MD->NumEle]+DummyY[MD->Ele[MD->Riv[i].LeftEle - 1].nabr[j]-1 + 2*MD->NumEle]+MD->Ele[MD->Ele[MD->Riv[i].LeftEle - 1].nabr[j]-1].zmin-MD->Ele[MD->Riv[i].LeftEle - 1].zmin)/2;
00815                                 }
00816                            }
00817                         Avg_BedDepth = (MD->Ele[MD->Riv[i].RightEle-1].zmax-MD->Ele[MD->Riv[i].RightEle-1].zmin + MD->Ele[MD->Riv[i].LeftEle-1].zmax-MD->Ele[MD->Riv[i].LeftEle-1].zmin)*0.5;
00818                         if(Avg_BedDepth - MD->Riv[i].depth <= 0)
00819                         {
00820                             MD->FluxSub[MD->Riv[i].LeftEle - 1][j] = 0.0;
00821                         }
00822                         else{
00823                             if (Avg_Y_Sub > Avg_BedDepth - MD->Riv[i].depth){
00824                                 MD->FluxSub[MD->Riv[i].LeftEle - 1][j] = MD->FluxSub[MD->Riv[i].LeftEle - 1][j] * (Avg_BedDepth - MD->Riv[i].depth) / Avg_Y_Sub;
00825                             }
00826                         }
00827                         break;
00828                    }
00829               }
00830          }
00831 
00832          if (MD->Riv[i].RightEle > 0){
00833          /************************TRANSFER THIS INSIDE FUNC*******************/
00834               if((MD->Ele[MD->Riv[i].RightEle - 1].zmax-MD->Ele[MD->Riv[i].RightEle - 1].zmin)-DummyY[MD->Riv[i].RightEle - 1+2*MD->NumEle]-DummyY[MD->Riv[i].RightEle - 1+MD->NumEle]<=0){
00835                    elemSatn=1.0;
00836               }
00837               else{
00838                    elemSatn = DummyY[MD->Riv[i].RightEle - 1+2*MD->NumEle]/(MD->Ele[MD->Riv[i].RightEle - 1].zmax-MD->Ele[MD->Riv[i].RightEle - 1].zmin);   /*  Will have to change this for other formulation */
00839               }
00840 
00841               if (MD->Soil[(MD->Ele[MD->Riv[i].RightEle - 1].soil-1)].Macropore == 0){
00842                    mp_factor = 1;
00843               }
00844               else if (MD->Soil[(MD->Ele[MD->Riv[i].RightEle - 1].soil-1)].Macropore == 1){
00845                    if((elemSatn>=sat_THRESH)&&(DummyY[MD->Riv[i].RightEle - 1]>ovl_THRESH_H)){
00846                         mp_factor=1.0+mp_MULTFH*(elemSatn-sat_THRESH)/(1-sat_THRESH);
00847                         mp_factor=mp_factor*mpArea_CALIB+1*(1-mpArea_CALIB);
00848                    }
00849                    else{
00850                         mp_factor = 1.0;
00851                    }
00852               }
00853               /*******************************************************************/
00854               if( (MD->Riv[i].zmin + DummyY[i+3*MD->NumEle] - (DummyY[MD->Riv[i].RightEle - 1 + 2*MD->NumEle]+MD->Ele[MD->Riv[i].RightEle - 1].zmin))>0 ){
00855                       loc_perem=Perem;
00856               }
00857               else{
00858                     if(MD->Ele[MD->Riv[i].RightEle - 1].zmin < MD->Riv[i].zmin){
00859                           if( (DummyY[MD->Riv[i].RightEle - 1 + 2*MD->NumEle] + MD->Ele[MD->Riv[i].RightEle - 1].zmin - MD->Riv[i].zmin) > MD->Riv[i].depth ){
00860                               loc_perem = CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[i].shape - 1].interpOrd,MD->Riv[i].depth,MD->Riv_Shape[MD->Riv[i].shape - 1].coeff,2);
00861 
00862                         }
00863                           else{
00864                               loc_perem = CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[i].shape - 1].interpOrd,(DummyY[MD->Riv[i].RightEle - 1 + 2*MD->NumEle] + MD->Ele[MD->Riv[i].RightEle - 1].zmin - MD->Riv[i].zmin),MD->Riv_Shape[MD->Riv[i].shape - 1].coeff,2);
00865                         }
00866                     }
00867                     else{
00868                         if( DummyY[MD->Riv[i].RightEle - 1 + 2*MD->NumEle] > MD->Riv[i].depth){
00869                             loc_perem = CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[i].shape - 1].interpOrd,MD->Riv[i].depth,MD->Riv_Shape[MD->Riv[i].shape - 1].coeff,2);
00870                         }
00871                         else{
00872                             loc_perem = CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[i].shape - 1].interpOrd,DummyY[MD->Riv[i].RightEle - 1 + 2*MD->NumEle],MD->Riv_Shape[MD->Riv[i].shape - 1].coeff,2);
00873                         }
00874                     }
00875               }
00876 
00877               GWflowFromEleToRiv(DummyY[MD->Riv[i].RightEle - 1 + 2*MD->NumEle],MD->Ele[MD->Riv[i].RightEle - 1].zmax,MD->Ele[MD->Riv[i].RightEle - 1].zmin,MD->Riv[i].x,MD->Ele[MD->Riv[i].RightEle - 1].x,MD->Riv[i].y,MD->Ele[MD->Riv[i].RightEle - 1].y,MD->Soil[(MD->Ele[MD->Riv[i].RightEle - 1].soil-1)].Macropore,DummyY[i+3*MD->NumEle],TotalY_Riv,MD->FluxRiv,i,5,MD->Riv[i].Length,MD->Soil[(MD->Ele[MD->Riv[i].RightEle - 1].soil-1)].base,mp_factor,loc_perem,MD->Ele[MD->Riv[i].RightEle - 1].Ksat,MD->Ele[MD->Riv[i].RightEle-1].RzD); /* delete replace Wid by 0.5*avg_perim */
00878 
00879               /* Saturation check */
00880               if((DummyY[MD->Riv[i].RightEle - 1 + 2*MD->NumEle] >= MD->Ele[MD->Riv[i].RightEle - 1 + 2*MD->NumEle].zmax-MD->Ele[MD->Riv[i].RightEle - 1 + 2*MD->NumEle].zmin) && MD->FluxRiv[i][5] > 0){
00881                    MD->FluxRiv[i][5]  = 0;
00882               }
00883 
00884               /*     Correction is being done in flux terms which can be > 0 even when there is no source water level present */
00885               if(DummyY[i + 3*MD->NumEle] <= 0 && MD->FluxRiv[i][5] > 0){
00886                    MD->FluxRiv[i][5] = 0;
00887               }
00888               if(DummyY[MD->Riv[i].RightEle - 1 + 2*MD->NumEle] <= 0 && MD->FluxRiv[i][5] < 0){
00889                    MD->FluxRiv[i][5] = 0;
00890               }
00891 
00892               /* modify groundwater flux item */
00893               for(j=0; j < 3; j++){
00894                   if(MD->Ele[MD->Riv[i].RightEle - 1].nabr[j] == MD->Riv[i].LeftEle){
00895                         if(MD->Ele[MD->Ele[MD->Riv[i].RightEle - 1].nabr[j]-1].zmin>MD->Ele[MD->Riv[i].RightEle - 1].zmin){
00896                                 if(MD->Ele[MD->Ele[MD->Riv[i].RightEle - 1].nabr[j]-1].zmin>MD->Ele[MD->Riv[i].RightEle - 1].zmin+DummyY[MD->Riv[i].RightEle - 1+2*MD->NumEle]){
00897                                    Avg_Y_Sub=DummyY[MD->Ele[MD->Riv[i].RightEle - 1].nabr[j]-1 + 2*MD->NumEle]/2;
00898                              }
00899                                 else{
00900                                       Avg_Y_Sub=(DummyY[MD->Riv[i].RightEle - 1+2*MD->NumEle]+MD->Ele[MD->Riv[i].RightEle - 1].zmin-MD->Ele[MD->Ele[MD->Riv[i].RightEle - 1].nabr[j]-1].zmin+DummyY[MD->Ele[MD->Riv[i].RightEle - 1].nabr[j]-1 + 2*MD->NumEle])/2;
00901                              }
00902                         }
00903                            else{
00904                                 if(MD->Ele[MD->Riv[i].RightEle - 1].zmin>MD->Ele[MD->Ele[MD->Riv[i].RightEle - 1].nabr[j]-1].zmin+DummyY[MD->Ele[MD->Riv[i].RightEle - 1].nabr[j]-1 + 2*MD->NumEle]){
00905                                    Avg_Y_Sub=DummyY[MD->Riv[i].RightEle - 1+2*MD->NumEle]/2;
00906                                 }
00907                                 else{
00908                                       Avg_Y_Sub=(DummyY[MD->Riv[i].RightEle - 1+2*MD->NumEle]+DummyY[MD->Ele[MD->Riv[i].RightEle - 1].nabr[j]-1 + 2*MD->NumEle]+MD->Ele[MD->Ele[MD->Riv[i].RightEle - 1].nabr[j]-1].zmin-MD->Ele[MD->Riv[i].RightEle - 1].zmin)/2;
00909                                 }
00910                            }
00911                           Avg_BedDepth = (MD->Ele[MD->Riv[i].LeftEle-1].zmax-MD->Ele[MD->Riv[i].LeftEle-1].zmin + MD->Ele[MD->Riv[i].RightEle-1].zmax-MD->Ele[MD->Riv[i].RightEle-1].zmin)*0.5;
00912                         if(Avg_BedDepth - MD->Riv[i].depth <= 0){
00913                             MD->FluxSub[MD->Riv[i].RightEle - 1][j] = 0.0;
00914                         }
00915                         else{
00916                               if (Avg_Y_Sub > Avg_BedDepth - MD->Riv[i].depth){
00917                                   MD->FluxSub[MD->Riv[i].RightEle - 1][j] = MD->FluxSub[MD->Riv[i].RightEle - 1][j] * (Avg_BedDepth - MD->Riv[i].depth) / Avg_Y_Sub;
00918                             }
00919                           }
00920                         break;
00921                   }
00922               }
00923          }
00924          /**************************************************************************************************/
00925          if (MD->Riv[i].LeftEle > 0){
00926               DummyDY[MD->Riv[i].LeftEle-1 + 2*MD->NumEle] = DummyDY[MD->Riv[i].LeftEle-1 + 2*MD->NumEle] + MD->FluxRiv[i][4]/MD->Ele[MD->Riv[i].LeftEle-1].area;
00927          }
00928          if (MD->Riv[i].RightEle > 0){
00929               DummyDY[MD->Riv[i].RightEle-1 + 2*MD->NumEle] = DummyDY[MD->Riv[i].RightEle-1 + 2*MD->NumEle]
00930                                  + MD->FluxRiv[i][5]/MD->Ele[MD->Riv[i].RightEle-1].area;
00931          }
00932 
00933          /**************************************************************************************************/
00934           if((isnan(MD->FluxRiv[i][0])==1)||(isnan(MD->FluxRiv[i][1])==1)||(isnan(MD->FluxRiv[i][2])==1)||(isnan(MD->FluxRiv[i][3])==1)||(isnan(MD->FluxRiv[i][4])==1)||(isnan(MD->FluxRiv[i][5])==1)){
00935               printf("\n2: %d %d %d %d :%lf %lf: %lf %lf: %lf %lf: %lf %lf %lf %lf %lf %lf",MD->Riv[i].index,MD->Riv[MD->Riv[i].down-1].index,MD->Riv[i].LeftEle,MD->Riv[i].RightEle,DummyY[i+3*MD->NumEle],DummyY[MD->Riv[i].down-1+3*MD->NumEle],DummyY[MD->Riv[i].LeftEle-1 + 2*MD->NumEle],DummyY[MD->Riv[i].LeftEle-1 + 0*MD->NumEle],DummyY[MD->Riv[i].RightEle-1 + 2*MD->NumEle],DummyY[MD->Riv[i].RightEle-1 + 0*MD->NumEle],MD->FluxRiv[i][0],MD->FluxRiv[i][1],MD->FluxRiv[i][2],MD->FluxRiv[i][3],MD->FluxRiv[i][4],MD->FluxRiv[i][5]);
00936               getchar();
00937          }
00938     }
00939 
00940 
00941     for(i=0; i<MD->NumEle; i++){
00942          for(j=0; j<3; j++){
00943               DummyDY[i] =  DummyDY[i] - MD->FluxSurf[i][j]/MD->Ele[i].area;
00944               DummyDY[i+2*MD->NumEle] = DummyDY[i+2*MD->NumEle] - MD->FluxSub[i][j]/MD->Ele[i].area;
00945          }
00946          /*
00947          if(i==144){
00948               //printf("\nIn Flux %f %f %f %f %f %f %f",t,DummyY[i],DummyY[i+MD->NumEle],DummyY[i+2*MD->NumEle],DummyDY[i],DummyDY[i+MD->NumEle],DummyDY[i+2*MD->NumEle]);
00949          }
00950          */
00951 
00952         //      AquiferDepth = MD->Ele[i].zmax - MD->Ele[i].zmin;
00953         //      Deficit = AquiferDepth - DummyY[i+MD->NumEle];
00954         //      G = MD->Ele[i].Porosity*(1-pow( 1+pow(MD->Ele[i].Alpha*(Deficit+EPSILON), MD->Ele[i].Beta),-(MD->Ele[i].Beta+1)/MD->Ele[i].Beta));
00955 
00956 
00957         /* for test case 1 only */
00958         /*
00959         G = MD->Ele[i].Porosity;
00960         GI = 0;
00961         */
00962         /*        if(i==160)
00963         {
00964                   printf("\nt= %f y =%.9f sy= %.9f",t,DummyY[i],DummyY[i+MD->NumEle]);
00965         }
00966         */
00967         //DummyDY[i+MD->NumEle] = DummyDY[i+MD->NumEle]/MD->Ele[i].Porosity;        /* ## Explicitly define porosity */
00968 
00969         AquiferDepth = MD->Ele[i].zmax - MD->Ele[i].zmin;
00970         Deficit = AquiferDepth - DummyY[i+2*MD->NumEle];
00971         //      PH = 1 - exp(-MD->Ele[i].Alpha*Deficit);
00972 
00973         //     elemSatn = Deficit-DummyY[i+MD->NumEle]>0?DummyY[i+MD->NumEle]/Deficit:(Deficit==0?1.0:0.0);
00974         if((MD->Ele[i].zmax-MD->Ele[i].zmin)-DummyY[i+2*MD->NumEle]-DummyY[i+MD->NumEle]<=0){
00975                 elemSatn=1.0;
00976         }
00977         else{
00978              elemSatn = 0.5*(1-cos(3.14*(DummyY[i+MD->NumEle]/(MD->Ele[i].zmax-MD->Ele[i].zmin-DummyY[i+2*MD->NumEle]))));    /*  Will have to change this for other formulation */
00979         }
00980 
00981         MD->Recharge[i] = elemSatn==0.0?0:Rec_CALIB*(-MD->Ele[i].Ksat*elemSatn*AquiferDepth*(1-(2*log(elemSatn)/(AquiferDepth*MD->Ele[i].Alpha)))/((AquiferDepth-DummyY[i+2*MD->NumEle])+DummyY[i+2*MD->NumEle]*elemSatn));
00982 
00983         if(DummyY[i+MD->NumEle]<=0 && MD->Recharge[i] < 0){
00984              MD->Recharge[i] = 0;
00985         }
00986         if(DummyY[i+2*MD->NumEle]<=0 && MD->Recharge[i] > 0){
00987              MD->Recharge[i] = 0;
00988         }
00989 
00990         DummyDY[i+MD->NumEle] = DummyDY[i+MD->NumEle] + MD->Recharge[i];
00991         DummyDY[i+MD->NumEle] = DummyDY[i+MD->NumEle]/MD->Ele[i].Porosity;
00992 
00993         /* Source check */
00994         /*      if(DummyY[i+MD->NumEle]>=Deficit-EPSILON && DummyDY[i+MD->NumEle]>0)
00995         {
00996         DummyDY[i+MD->NumEle] = 0;
00997         }
00998         if(DummyY[i+MD->NumEle]<=0 && DummyDY[i+MD->NumEle]<0)
00999         {
01000         DummyDY[i+MD->NumEle] = 0;
01001         }
01002         */
01003         DummyDY[i+2*MD->NumEle] = DummyDY[i+2*MD->NumEle]-MD->Recharge[i];
01004         DummyDY[i+2*MD->NumEle] = DummyDY[i+2*MD->NumEle]/MD->Ele[i].Porosity;
01005 
01006         if(DummyY[i+2*MD->NumEle]>=AquiferDepth && DummyDY[i+2*MD->NumEle]>0){
01007              DummyDY[i+2*MD->NumEle] = 0;
01008         }
01009         if(DummyY[i+2*MD->NumEle]<=0 && DummyDY[i+2*MD->NumEle]<0){
01010              DummyDY[i+2*MD->NumEle] = 0;
01011         }
01012         if(DummyY[i+MD->NumEle]>=AquiferDepth && DummyDY[i+MD->NumEle]>0){  // WHAT?
01013              DummyDY[i+MD->NumEle] = 0;
01014         }
01015         if(DummyY[i+MD->NumEle]<=0 && DummyDY[i+MD->NumEle]<0){
01016              DummyDY[i+MD->NumEle] = 0;
01017         }
01018         if(DummyY[i]<=0 && DummyDY[i]<0){
01019              DummyDY[i] = 0;
01020         }
01021         // if there is a source (well) in it
01022         /*      if(MD->Ele[i].source > 0)
01023         {
01024         DummyDY[i+MD->NumEle] = DummyDY[i+MD->NumEle] - Interpolation(&MD->TSD_Source[MD->Ele[i].source - 1], t)/(MD->Ele[i].Porosity*MD->Ele[i].area);
01025         }
01026         */
01027         // check if it is out of bound
01028 
01029         /*      if(Y[i+MD->NumEle]>AquiferDepth && DummyDY[i+MD->NumEle]>0)
01030         {
01031         DummyDY[i+MD->NumEle] = 0;
01032         }
01033         if(Y[i+MD->NumEle]<0 && DummyDY[i+MD->NumEle]<0)
01034         {
01035         DummyDY[i+MD->NumEle] = 0;
01036         }
01037         */
01038     }
01039     for(i=0; i<MD->NumRiv; i++){
01040          DummyDY[i+3*MD->NumEle] = MD->FluxRiv[i][0] - MD->FluxRiv[i][1];
01041          DummyDY[i+3*MD->NumEle] = DummyDY[i+3*MD->NumEle] - MD->FluxRiv[i][2] - MD->FluxRiv[i][3];
01042          DummyDY[i+3*MD->NumEle] = DummyDY[i+3*MD->NumEle] - MD->FluxRiv[i][4] - MD->FluxRiv[i][5];
01043          DummyDY[i+3*MD->NumEle] = DummyDY[i+3*MD->NumEle]/(MD->Riv[i].Length*CS_AreaOrPerem(MD->Riv_Shape[MD->Riv[i].shape - 1].interpOrd,DummyY[i + 3*MD->NumEle],MD->Riv_Shape[MD->Riv[i].shape - 1].coeff,3)); /* delete derive denominator to be replace by volume */
01044          if(DummyY[i+3*MD->NumEle]<=0 && DummyDY[i+3*MD->NumEle]<0){
01045               DummyDY[i+3*MD->NumEle] = 0;
01046          }
01047     }
01048     /*  for(i=0;i<MD->NumEle;i++)
01049         {
01050     printf("{");
01051     for(j=0;j<3;j++)
01052         {
01053         printf("%lf ",MD->FluxSub[i][j]);
01054         }
01055     printf("},");
01056     }
01057     for(i=0;i<MD->NumEle;i++)
01058     {
01059     printf("{");
01060     for(j=0;j<3;j++)
01061         {
01062         printf("%lf ",MD->FluxSurf[i][j]);
01063         }
01064     printf("},");
01065     }
01066     */
01067 
01068     //printf("Flux: %lf %lf %lf %lf %lf %lf\n", MD->FluxSub[5067-2*MD->NumEle][0], MD->FluxSub[5067-2*MD->NumEle][1], MD->FluxSub[5067-2*MD->NumEle][2], MD->Recharge[5067-2*MD->NumEle],MD->FluxRiv[386][4], MD->FluxRiv[386][5]);
01069     for(i=0; i<3*MD->NumEle+MD->NumRiv; i++){
01070          /* Y[i]=DummyY[i]; */
01071            DY[i]=DummyDY[i]/(60.0*24.0);
01072     //printf("%d\t%lf\n", i, DY[i]);
01073     }
01074 
01075     free(DummyY);
01076     free(DummyDY);
01077     return(0);
01078 }
01079 
01080 
01081 realtype Interpolation(TSD *Data, realtype t)
01083 
01086 {
01087     int i, success;
01088     realtype result;
01089 
01090     //i=0;
01091     i=Data->iCounter;
01092     success = 0;
01093     t=t/(24.0*60.0);
01094     while(i<Data->length && t>Data->TS[i][0]){
01095          i++;
01096     }
01097 
01098     if(i==0){
01099          /* t is smaller than the 1st node */
01100          result = Data->TS[i][1];
01101     }
01102     else if(i >= Data->length){
01103          result = Data->TS[i-1][1];
01104     }
01105     else{
01106          result = ((Data->TS[i][0]-t)*Data->TS[i-1][1] + (t-Data->TS[i-1][0])*Data->TS[i][1])/(Data->TS[i][0]-Data->TS[i-1][0]);
01107          success = 1;
01108     }
01109 
01110     if(success == 0){
01111          /*
01112          printf("\nWarning:  Exterpolation is used ...\n");
01113          */
01114     }
01115 
01116     return result;
01117 }
01118 
01119 
01120 realtype returnVal(realtype rArea, realtype rPerem, realtype eqWid, realtype ap_Bool)
01122 
01127 {
01128     if(ap_Bool==1)
01129     {
01130         return rArea;
01131     }
01132     else if(ap_Bool==2)
01133     {
01134         return rPerem;
01135     }
01136     else
01137     {
01138         return eqWid;
01139     }
01140 }
01141 
01142 /*    Area or Peremeter of River Segment's cross-section    */
01143 realtype CS_AreaOrPerem(int rivOrder, realtype rivDepth, realtype rivCoeff, realtype a_pBool)
01145 
01150 {
01151     realtype rivArea, rivPerem, eq_Wid;
01152     switch(rivOrder)
01153     {
01154         case 1:
01155             rivArea = rivDepth*rivCoeff;
01156             rivPerem= 2.0*rivDepth+rivCoeff;
01157             eq_Wid=rivCoeff;
01158             return returnVal(rivArea, rivPerem, eq_Wid, a_pBool);
01159         case 2:
01160             rivArea = pow(rivDepth,2)/rivCoeff;
01161             rivPerem = 2.0*rivDepth*pow(1+pow(rivCoeff,2),0.5)/rivCoeff;
01162             eq_Wid=2.0*pow(rivDepth+EPSILON,1/(rivOrder-1))/pow(rivCoeff,1/(rivOrder-1));
01163             return returnVal(rivArea, rivPerem, eq_Wid, a_pBool);
01164         case 3:
01165             rivArea = 4*pow(rivDepth,1.5)/(3*pow(rivCoeff,0.5));
01166             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));
01167             eq_Wid=2.0*pow(rivDepth+EPSILON,1/(rivOrder-1))/pow(rivCoeff,1/(rivOrder-1));
01168             return returnVal(rivArea, rivPerem, eq_Wid, a_pBool);
01169         case 4:
01170             rivArea = 3*pow(rivDepth,4.0/3.0)/(2*pow(rivCoeff,1.0/3.0));
01171             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))));
01172             eq_Wid=2.0*pow(rivDepth+EPSILON,1/(rivOrder-1))/pow(rivCoeff,1/(rivOrder-1));
01173             return returnVal(rivArea, rivPerem, eq_Wid, a_pBool);
01174         default:
01175             printf("\n Relevant Values entered are wrong");
01176             printf("\n Depth: %lf\tCoeff: %lf\tOrder: %d\t");
01177             return 0;
01178     }
01179 }
01180 
01181 /*    Function to compute surface-flux between elements    */
01182 void OverlandFlow(realtype **flux, 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) /* delete perimeter should be passed */
01184 
01198 {
01199     int locBool;
01200     float hydRadius;
01201     /* if surface gradient is not enough to overcome the friction */
01202     if(fabs(grad_y) <= avg_sf)
01203     {
01204          flux[loci][locj] = 0;
01205     }
01206     else
01207     {
01208          if(grad_y > 0)
01209          {
01210               locBool=1;
01211          }
01212          else if(grad_y < 0)
01213          {
01214               locBool=-1;
01215          }
01216          switch(surfmode)
01217          {
01218                 case 1:
01219                   if(eletypeBool==1)
01220                   {
01221                        /* Kinematic Wave Approximation constitutive relationship: Manning Equation */
01222                        alfa = sqrt(locBool*grad_y)/avg_rough;
01223                        beta = pow(avg_y, 2.0/3.0);
01224                        flux[loci][locj] = locBool*alfa*beta*crossA;
01225                        break;
01226                   }
01227                   else
01228                   {
01229                        /*alfa = sqrt(locBool*grad_y)/(avg_rough*pow((avg_perem>0?avg_perem:0), 2.0/3.0));
01230                        beta = 5.0/3.0;
01231                        flux[loci][1] = locBool*alfa*pow(crossA, beta);
01232                        */
01233                        hydRadius = (avg_perem>0?crossA/avg_perem:0);
01234                        flux[loci][1] = locBool*sqrt(locBool*grad_y)*crossA*pow(hydRadius,2.0/3.0)/avg_rough;
01235                        break;
01236                   }
01237                 case 2:
01238                   if(eletypeBool==1)
01239                   {
01240                        /* Diffusion Wave Approximation constitutive relationship: Gottardi & Venutelli, 1993 */
01241                        alfa = pow(pow(avg_y, 1.0/3.0),2)/(1.0*avg_rough);
01242                        beta = alfa;
01243                        flux[loci][locj] = locBool*crossA*beta*sqrt(locBool*grad_y);
01244                        break;
01245                   }
01246                   else
01247                   {
01248                        /*alfa = pow(pow(avg_y, 1.0/3.0),2)/(1.0*avg_rough);
01249                        beta = alfa;
01250                        flux[loci][1] = locBool*crossA*beta*sqrt(locBool*grad_y);
01251                        */
01252                        hydRadius = (avg_perem>0?crossA/avg_perem:0);
01253                        flux[loci][1] = locBool*sqrt(locBool*grad_y)*crossA*pow(hydRadius,2.0/3.0)/(1.0*avg_rough);
01254                        break;
01255                   }
01256                 default:
01257                   if(eletypeBool==1)
01258                   {
01259                           printf("Fatal Error: Surface Overland Mode Type Is Wrong!");
01260                   }
01261                   else
01262                   {
01263                        printf("Fatal Error: River Routing Mode Type Is Wrong!");
01264                   }
01265                   exit(1);
01266         }
01267     }
01268 }
01269 
01270 
01271 /*    Surface Interaction between Element and River Segment    */
01272 void OLflowFromEleToRiv(realtype sideEle_y,realtype sideEle_zmax,realtype rivX,realtype sideEleX,realtype rivY,realtype sideEleY,realtype cwr,realtype rivZmax,realtype loc_yriver,realtype **fluxriv,int loc_i,int loc_j,realtype length)
01274 
01289 {
01290       realtype loc_cwr,dist,ele_YH,ele_Y,loc_bele;
01291       ele_Y = sideEle_y;
01292       ele_YH = sideEle_y + sideEle_zmax;
01293       dist = sqrt(pow(rivX - sideEleX, 2) + pow(rivY - sideEleY, 2));
01294       loc_cwr = cwr;
01295       if (rivZmax < sideEle_zmax)
01296       {
01297            loc_bele = sideEle_zmax;
01298       }
01299       else
01300       {
01301            loc_bele = rivZmax;
01302       }
01303 
01304       /*  #? where is the reference to this */
01305       /* ## i think its in Sorab PAnday's paper */
01306 
01307       if (loc_yriver > ele_YH)
01308       {
01309               if (ele_YH > loc_bele)
01310             {
01311                   fluxriv[loc_i][loc_j] = loc_cwr*2.0*sqrt(2*GRAV)*length*sqrt(loc_yriver - ele_YH)*(loc_yriver - loc_bele)/3.0;
01312             }
01313             else
01314             {
01315                 if(loc_bele<loc_yriver)
01316                 {
01317                       fluxriv[loc_i][loc_j] = loc_cwr*2.0*sqrt(2*GRAV)*length*sqrt(loc_yriver - loc_bele)*(loc_yriver - loc_bele)/3.0;
01318                      }
01319                 else
01320                   {
01321                       fluxriv[loc_i][loc_j]=0.0;
01322                   }
01323               }
01324 
01325             /*fluxriv[loc_i][loc_j]=0.0; */ //uncommented last if
01326       }
01327       else{
01328            if (loc_yriver > loc_bele)
01329            {
01330                 fluxriv[loc_i][loc_j] = -loc_cwr*2.0*sqrt(2*GRAV)*length*sqrt(ele_YH - loc_yriver)*(ele_YH - loc_bele)/3.0;
01331            }
01332            else
01333            {
01334                 /* This is basicaly a lump representation */
01335                 if(loc_bele<ele_YH)
01336                 {
01337                      fluxriv[loc_i][loc_j] = -loc_cwr*2.0*sqrt(2*GRAV)*length*sqrt(ele_YH - loc_bele)*(ele_YH - loc_bele)/3.0;
01338                 }
01339                 else
01340                 {
01341                      fluxriv[loc_i][loc_j]=0.0;
01342                 }
01343            }
01344       }
01345 }
01346 
01347 /*    SubSurface Interaction between Element and River Segment    */
01348 void GWflowFromEleToRiv(realtype sideEle_y,realtype sideEle_zmax, realtype sideEle_zmin,realtype rivX,realtype sideEleX,realtype rivY,realtype sideEleY,int loc_McPore,realtype loc_yriver,realtype loc_totyriver,realtype **fluxriv,int loc_i,int loc_j,realtype length, realtype loc_base,realtype loc_gama, realtype loc_perem,realtype loc_ksat,realtype ele_Thresh) /* delete 0.5 perimeter */
01350 
01370 {
01371      realtype loc_mpfactor,dist,ele_YH,ele_Y;
01372 
01373     realtype loc_sat, loc_rivK_CALIB, mp_Rzd=0.8;
01374     mp_Rzd = ((sideEle_zmax - sideEle_zmin - mp_Rzd) < 0)?(sideEle_zmax-sideEle_zmin):mp_Rzd;
01375     loc_sat = (sideEle_y-(sideEle_zmax-sideEle_zmin-mp_Rzd))/mp_Rzd;
01376     loc_rivK_CALIB = loc_sat>=0?(1.0+loc_sat*rivK_CALIB):1.0;
01377 
01378      ele_Y = sideEle_y;
01379      ele_YH = sideEle_y + sideEle_zmin;
01380      dist = sqrt(pow(rivX - sideEleX, 2) + pow(rivY - sideEleY, 2));
01381 
01382 
01383      if (loc_McPore == 0){
01384           loc_mpfactor = 1;
01385      }
01386      else if (loc_McPore == 1){
01387           loc_mpfactor = loc_gama;
01388      }
01389 
01390      //fluxriv[loc_i][loc_j] = length*(0.5*loc_perem)* loc_ksat *rivK_CALIB* ((loc_totyriver-ele_YH)>0? (loc_yriver>0? ((loc_totyriver-loc_yriver)>(ele_YH+ele_Thresh)?ele_Thresh:(loc_totyriver-ele_YH)):0): (ele_Y>0?((sideEle_zmin-loc_totyriver)>0?(-1.0*ele_Y):(loc_totyriver- ele_YH)):0))/dist; /* delete 0.5* perimeter */
01391         fluxriv[loc_i][loc_j] = length*(0.5*loc_perem)* loc_ksat *loc_rivK_CALIB* ((loc_totyriver-ele_YH)>0? (loc_yriver>0? ((loc_totyriver-loc_yriver)>(ele_YH+ele_Thresh)?ele_Thresh:(loc_totyriver-ele_YH)):0): (ele_Y>0?((sideEle_zmin-loc_totyriver)>0?(loc_totyriver- ele_YH):(loc_totyriver- ele_YH)):0))/dist;
01392         fluxriv[loc_i][loc_j] = fluxriv[loc_i][loc_j]>0?fluxriv[loc_i][loc_j]:fluxriv[loc_i][loc_j]*loc_mpfactor;
01393 }
Generated on Thu Jul 12 14:34:18 2007 for PIHM by  doxygen 1.5.2