pihm/f.cpp

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