A.1: Codes for Vadose Zone Leaching and Saturated Zone Mixing A.2 Codes of Groundwater Flow Appendix A. SOURCE CODES FOR VG MODELA.1 Codes for Vadose Zone Leaching and Saturated Zone Mixing
//************************************** // Main.cpp, By Sam Sang Lee (URI), 1998 // ************************************* #include "VadClass.h" void main(int x, char **FName){ vadoseClass V; V.VadoseOutputFiles(); }//end main //************************************** // VadClass.cpp //************************************** #include "VadClass.h" float Abs(float Num) { if (Num < 0) Num = Num * -1; return (Num); }//end Abs vadoseClass::vadoseClass(){ Vad.ReadInput("Data.inp"); }//end default constructor vadoseClass::vadoseClass(char *FileName){ Vad.ReadInput(FileName); }//end default constructor void vadoseClass::VadoseOutputFiles(){ OpenFiles(); Prm.Output(Vad); MainLoop(); CloseFiles(); }//end VadoseOutputFiles void vadoseClass::OpenFiles(){ char FileName[255]; int Len, Cut; strcpy(FileName, Vad.OutFileName); Len = strlen(FileName); Cut = Len - 3; FileName[Cut] = '\0'; strcat(FileName, "prm"); Prm.OpenFile(FileName); FileName[Cut] = '\0'; strcat(FileName, "out"); Out.OpenFile(FileName); FileName[Cut] = '\0'; strcat(FileName, "prf"); Prf.OpenFile(FileName); FileName[Cut] = '\0'; strcat(FileName, "sim"); Sim.OpenFile(FileName); FileName[Cut] = '\0'; strcat(FileName, "gim"); Gim.OpenFile(FileName); FileName[Cut] = '\0'; strcat(FileName, "cmx"); Cmx.OpenFile(FileName); }//end OpenFiles void vadoseClass::CloseFiles(){ Prm.CloseFile(); Out.CloseFile(); Prf.CloseFile(); Sim.CloseFile(); Gim.CloseFile(); Cmx.CloseFile(); }//end CloseFiles //BEGINNING OF MAIN CALCULATIONS AND //OUTPUT void vadoseClass::MainLoop(){ int i,j; //loop controls Out.SimTitle(Vad.SimTitle); for (i=0; i < Vad.Ntime; ++i){ Gtotal[i] = 0; }//end for i Cmx.Title(Vad.Angl); for (i=0; i < Vad.Nscol; ++i){ if (Vad.Angl == 'T') Cmx.ColInfo(Vad.Angl, i+1); VadMain(i); for(j=0; j < Vad.Ntime; ++j){ Gwim[i][j] = Gwimp[j]; Gtot[i][j] = Gtotal[j]; }//end forj }//end for i Gcum = 0; Out.TotalGWTitle(); for(i=0; i < Vad.Ntime; ++i){ Gwimtl = 0; Gtottl = 0; for(j=0; j < Vad.Nscol; ++j) { Gwimtl += Gwim[j][i]; Gtottl += Gtot[j][i]; }//end for j Gcum += Gtottl; Out.TotalGWImpact(Vad.Ntime, (i+1)*Vad.Ptime, Gwimtl, Gcum); }//end for i if (Vad.Angl == 'P'){ Cmx.ColInfo(Vad.Angl, 1); Qsum[0] = Vad.Volfxc[0] + Vad.Pdepth[0] * Vad.Width[0] * Vad.Flxaqf[0]; for(j=0; j <= Vad.Ntimest; ++j) { Time = (j+1) * Vad.Delt; Cmixe[0][j] = (Vad.Pdepth[0] * Vad.Width[0] * Vad.Flxaqf[0]* Vad.Cinpaq[0] + Vad.Volfxc[0] * Cliqend[0][j])/Qsum[0]; Cmx.PAngleOutput_1(Time, Cliqend[0][j] / CONVERT_2, Vad.Cinpaq[0] / CONVERT_2, Cmixe[0][j] / CONVERT_2); }//end for j for(i=1; i < Vad.Nscol; ++i) { Cmx.ColInfo(Vad.Angl, (i+1)); Qsum[i]=Vad.Volfxc[i]+ Vad.Pdepth[i] * Vad.Width[i] * Vad.Flxaqf[i] + Qsum[i-1]; for(j=0; j <= Vad.Ntimest; ++j) { Time = (j+1) * Vad.Delt; Cmixe[i][j] = (Vad.Pdepth[i] * Vad.Width[i]*Vad.Flxaqf[i]*Vad.Cinpaq[i] +Vad.Volfxc[i]*Cliqend[i][j]+Cmixe[i-1] [j] * Qsum[i-1]) / Qsum[i]; Cmx.PAngleOutput_2(Time, Cliqend[i] [j]/CONVERT_2, Vad.Cinpaq[i] / CONVERT_2, Cmixe[i-1][j] / CONVERT_2, Cmixe[i][j] / CONVERT_2); }//end for j }//end for i }//end if P Sim.Terminate(); Gim.Terminate(); }//end MainLoop void vadoseClass::VadMain(int Col){ int i,j, Ncel; int count = 0; Nceldm = Vad.Ncell[Col]; Ncel = Vad.Ncell[Col] - 1; for(i=0; i <= Ncel; ++i){ Cgas[i] = 0; Cliq[i] = 0; Cs[i] = 0; }//end for i for(i=0; i<Vad.GrNum[Col]; ++i){ for(j=(Vad.J1[Col][i]-1); j < Vad.J2[Col][i]; ++j) { if(Vad.Kd[Col][j] == 0) { Cliq[j] = Vad.Xcon[Col][i] * CONVERT_2; Cs[j] = 0; }//end if else { Cs[j] = float(Vad.Xcon[Col][i] * 0.000001); Cliq[j] = Cs[j] / Vad.Kd[Col][j]; }//end else }//end for j }//end for i if (Vad.Typtbc[Col] == 1){ Cs[0]=Vad.Vatbc[Col] * Vad.Kd[Col][0]; Cliq[0] = Vad.Vatbc[Col]; }//end if //Call Equil function VadEquil(Col); Out.SoilTitle(Vad.SoilTitle[Col]); if (Iflag > 0) { Out.Warn(0, Iflag); }//end if Out.AtTime((Col+1), 0, Mtotal, Mgtot, Mltot, Mstot); Mt0 = Mtotal; Mtp = Mtotal; MdfSum = 0; MatCum = 0; MatInt = 0; MabCum = 0; MabInt = 0; MdtCum = 0; MdtInt = 0; MdbCum = 0; MdbInt = 0; MdcCum = 0; MdcInt = 0; Itime = 0; Prf.Title(Vad.SoilTitle[Col], 0); for(i=0; i < Nceldm; ++i){ Depthz = i * Vad.Delz[Col]; Prf.Output(Depthz, Cgas[i]/CONVERT_2, Cliq[i]/CONVERT_2, float(Cs[i]*(1e+6))); }//end for i Cofm=Vad.Delt/(2*Vad.Delz[Col]*Vad.Delz[Col]); Cofn = Vad.Delt / (4*Vad.Delz[Col]); if (Vad.Decy1[Col]) Li = Vad.Delt / 2*Vad.Decy1[Col]; else Li = 0; Mi[0]=Cofm/Vad.Theta[Col][0]*Vad.Disp[Col][0]; Ni[0]=Cofn/Vad.Theta[Col][0]*Vad.Q[Col][0]; Mip[0] = 0; Nip[0] = 0; for(i=1; i<(Nceldm-1); ++i){ Cofmi = Cofm / Vad.Theta[Col][i]; Cofni = Cofn / Vad.Theta[Col][i]; Mi[i] = Cofmi * Vad.Disp[Col][i]; Ni[i] = Cofni * Vad.Q[Col][i]; Mip[i] = Cofmi*( Vad.Disp[Col][i+1] - Vad.Disp[Col][i-1] ) / 4.0; Nip[i] = Cofni*( Vad.Q[Col][i+1] - Vad.Q[Col][i-1] ); }//end for i Mi[Ncel] = Cofm / Vad.Theta[Col][Ncel] * Vad.Disp[Col][Ncel]; Ni[Ncel] = Cofn / Vad.Theta[Col][Ncel] * Vad.Q[Col][Ncel]; Mip[Ncel] = 0; Nip[Ncel] = 0; for(i=1; i<(Nceldm-1); ++i) { A[i] = -Mi[i] + Mip[i] - Ni[i]; B[i] = 1 + 2*Mi[i] + Nip[i] + Li; C[i] = -Mi[i] - Mip[i] + Ni[i]; }//end for i Areas = Vad.Width[Col] * Vad.Pdepth[Col]; if (Vad.Plt[Col] == 'Y') Gim.Title(Col); Iptime = 1; Ipltime = 1; Iprtim = 1; //Loop Begins do { ++Itime; Time = Itime * Vad.Delt; if (Time > Vad.Timpl) Vad.Vatbc[Col] = 0; //Top boundary conditions B[0] = 1.; Ddzmpn1 = Vad.Disp[Col][0] / ( Vad.Delz[Col]*( Mi[0] + Ni[0]) ); Aaa1=Ddzmpn1*(2.*Mi[0]+Li+1)/ 4.; Bbb1 = Ddzmpn1 * Mi[0] / 2; Ccc1=Ddzmpn1*(2.*Mi[0]+Li-1.)/4.; Phi = Aaa1 + Vad.Q[Col][0] / 2.; Psi = Bbb1; Omg = Ccc1 + Vad.Q[Col][0] / 2.; if(Vad.Typtbc[Col] == 1) { C[0] = 0.; Rhs[0] = Vad.Vatbc[Col] * float(exp(-Vad.Decy2[Col] * Time)); } else if(Vad.Typtbc[Col] == 3) { C[0] = -Psi / Phi; Rhs[0] = (-Omg*Cliq[0] +Psi*Cliq[1]+Vad.Q[Col][0]*Vad.Vatbc[Col]*float(exp(-Vad.Decy2[Col]*Time)))/Phi; } else cout << "Check type of top boundary condition !!!\n"; //Bottom boundary conditions type 2: Cn - //Cn-1 = 0 Valbbc = 0; Rhs[Ncel] = Valbbc; B[Ncel] = 1.; A[Ncel] = -1.; //Values between the top and bottom for(i=1; i<Ncel; ++i) { Rhs[i] = -A[i] * Cliq[i-1] + ( 2. - B[i] ) * Cliq[i] - C[i]*Cliq[i+1]; }//end for i //Flux calculations Ddzmpnn = Vad.Disp[Col][Ncel] / ( Vad.Delz[Col] *( Mi[Ncel] + Ni[Ncel] ) ); Aaan = Ddzmpnn*(2.*Mi[Ncel]+Li+1.) / 4.; Bbbn = Ddzmpnn*Mi[Ncel] / 2.; Cccn = Ddzmpnn*(2.*Mi[Ncel]+Li-1.) / 4.; Mat00 = Vad.Q[Col][0] * Cliq[0]; Mdt00 = Ccc1*Cliq[0] - Bbb1*Cliq[1]; Mab00 = Vad.Q[Col][Ncel] * Cliq[Ncel]; Mdb00=Cccn*Cliq[Ncel-1]-Bbbn*Cliq[Ncel]; //Call to Tridia VadTridia(A, B, C, Rhs, Cliq, Ncel); //Top b' dry flux calculations by //Adevection Mat=Vad.Delt*(Vad.Q[Col][0]*Cliq[0]+ Mat00) / 2.; //Top b' dry flux calculations by Mdt=Vad.Delt*(Mdt00+Aaa1*Cliq[0]-Bbb1 *Cliq[1]); //Bottom b' dry flux calculations by //Adevection Mab=-Vad.Delt*(Vad.Q[Col][Ncel] *Cliq[Ncel] + Mab00) / 2.; //Bottom b' dry flux calculations by //Dispersion Mdb = -Vad.Delt * (Mdb00 + Aaan *Cliq[Ncel-1] - Bbbn*Cliq[Ncel]); //Mass equilibration step Mttold = Mtotal; VadEquil(Col); if (Iflag > 0) { Out.Warn(0, Iflag); }//end if Mdiff = Mtotal - Mttold; MdfSum += Mdiff; //Mass Decay in the soil column Mdc=-Vad.Decy1[Col] * Mtotal * Vad.Delt; //Mass Balance calculations MatCum += Mat; MabCum += Mab; MdtCum += Mdt; MdbCum += Mdb; MdcCum += Mdc; MatInt += Mat; MabInt += Mab; MdtInt += Mdt; MdbInt += Mdb; MdcInt += Mdc; //Write Groundwater impact data to plot //file if (Vad.Plt[Col] == 'Y') { Gim.NumOutput(Time, (-Vad.Areac[Col] * (Mab + Mdb) / Vad.Delt)); }//end if //Write soil concentration - depth data //to plot file Pltim = Vad.Pltime[Col] * Ipltime; if ((Vad.Plt[Col] == 'Y') && (Abs(Time-Pltim) < (1e-6))) { ++Ipltime; Sim.Title((Col+1), Time); for(i=0; i < Nceldm; ++i) { Sim.Output(float(Cliq[i] / CONVERT_2), (-(i) * Vad.Delz[Col])); }//end for i }//end if Ptime1 = int(Vad.Ptime * Iptime); if (Abs(Time - Ptime1) < (1e-6)) { Gwimp[count] = -Vad.Areac[Col] * (Mab + Mdb) / Vad.Delt; Gtotal[count] = -Vad.Areac[Col] * (MabInt + MdbInt); //Call to Output VadOutput(Col); ++Iptime; ++count; MdfSum = 0; }//end if //Write vert. concent profiles Prtim = int(Vad.Prtime * Iprtim); if (Abs(Time-Prtim) < 1e-6) { ++Iprtim; Prf.Title(Vad.SoilTitle[Col], Time); for(i=0; i < Nceldm; ++i) { Depthz = i * Vad.Delz[Col]; Prf.Output(Depthz,Cgas[i]/CONVERT_2, Cliq[i]/ CONVERT_2, float(Cs[i] * (1e+6))); }//end for i }//end if if (Vad.Angl == 'T') { Qaqf = Vad.Flxaqf[Col] * Areas; Cmixed = Qaqf*Vad.Cinaqf[Col]+Cliq [Ncel]*Vad.Areac[Col]*Vad.Q[Col][Ncel])/(Qaqf+Vad.Areac[Col]*Vad.Q[Col][Ncel]); Cmx.RAngleOutput(Time,Cliq[Ncel] /CONVERT_2,Vad.Cinaqf[Col]/CONVERT_2, Cmixed / CONVERT_2); }//end if Angl==T else { Cliqend[Col][Itime-1] = Cliq[Ncel]; }//end else }while(Time < Vad.Stime); //Loop Ends Out.GWImpact((Col+1), Vad.Ntime, Vad.Ptime, Gwimp, Vad.Areac[Col]); }//end VadMain void vadoseClass::VadOutput(int Col){ float Del1, Del2; Out.AtTime((Col+1), Time, Mtotal, Mgtot, Mltot, Mstot); Out.SinceLast(Time - Vad.Ptime); Del1 = Mtotal - Mtp; Del2 = MatInt + MdtInt + MabInt + MdcInt + MdbInt; Out.Change(Del1, MatInt + MdtInt, MabInt + MdbInt, MdcInt, Del2); Out.SinceBegin(0); Del1 = Mtotal - Mt0; Del2 = MatCum + MdtCum + MabCum + MdbCum + MdcCum; Out.Change(Del1, MatCum + MdtCum, MabCum + MdbCum, MdcCum, Del2); MatInt = 0; MdtInt = 0; MabInt = 0; MdbInt = 0; MdcInt = 0; Mtp = Mtotal; }//end VadOutput void vadoseClass::VadEquil(int Col){ Mgtot = 0; Mltot = 0; Mstot = 0; Mtotal = 0; Iflag = 0; for(int i=0; i < Nceldm; ++i) { Ml = Cliq[i] * Vad.Delz[Col] * Vad.Thetaw[Col][i]; Cgas[i] = Cliq[i] * Vad.Kh; Mg = Cgas[i] * Vad.Delz[Col] * Vad.Pora[Col][i]; Cs[i] = Cliq[i] * Vad.Kd[Col][i]; Ms=Cs[i]*Vad.Delz[Col]* Vad.Rhob[Col] [i]; Mgtot += Mg; Mltot += Ml; Mstot += Ms; Mt = Mg + Ml + Ms; if (i==1) Mtotal1 = Mt; Mtotal += Mt; if (Cliq[i] > Vad.Cmax) ++Iflag; }//end for i }//end VadEquil void vadoseClass::VadTridia(vadNumType A[], vadNumType B[], vadNumType C[], vadNumType F[], vadNumType X[], int N){ //Variables vadNumType Alpha, Beta[1000], Y[1000]; int i; //Calculating beta and y arrays Alpha = B[0]; Y[0] = F[0] / Alpha; for(i=1; i<=N; ++i) { Beta[i-1] = C[i-1] / Alpha; Alpha = B[i] - A[i]*Beta[i-1]; Y[i]=(F[i]-A[i] * Y[i-1] ) / Alpha; }//end for i //Back substitution to calculate x array X[N] = Y[N]; for(i=(N-1); i>=0; --i) X[i] = Y[i] - Beta[i]*X[i+1]; }//end VadTridia //************************************** // VadInput.cpp //************************************** #include "VadInput.h" vadoseInput::vadoseInput(){ }//end default constructor vadoseInput::vadoseInput(char *FileName){ ReadInput(FileName); }//end constructor void vadoseInput::ReadInput(char *FileName){ int i, j, k; char ch; float fl1,fl2; vadNumType Qj, Rhobij, Rhobj, Porj, Poraj, Thetawj, Thetaj, Focj, Dispj, Kdj, Darj; ifstream In(FileName); In.getline(OutFileName, MAX_CHAR-1, '\n'); In.getline(SimTitle, MAX_CHAR - 1, '\n'); In >> Nscol; In >> Delt >> Stime >> Timpl >> Ptime >> Prtime >> Angl; In >> Koci >> Kh >> Cmaxi >> Dairi; Koc = Koci / CONVERT_1; Cmax = Cmaxi * CONVERT_2; Dair = Dairi * CONVERT_3; Ntime = int(Stime/Ptime); Ntimest = int(Stime/Delt); //WRITING IN SOIL COLUMN DATA FOR EACH //NSCOL for(i=0; i < Nscol; ++i) { In.get(ch); In.get(ch); In.getline(SoilTitle[i], MAX_CHAR - 1, '\n'); In>>Ncell[i]>>Delz[i]>>Width[i] >>Length[i]>>Typtbc[i]>>Valtbc[i] >>Aldisp[i]; In >> Plt[i] >> Pltime[i] >> Decy1[i] >> Decy2[i]; In >> GrNum[i]; Areac[i] = Width[i] * Length[i]; Vatbc[i] = Valtbc[i] * CONVERT_2; //WRITING IN AQUIFER DATA FOR EACH NCELL for(j=0; j < GrNum[i]; ++j) { In >> J1[i][j] >> J2[i][j] >> Qj >> Rhobij >> Porj >> Thetawj >> Focj >>Xcon[i][j]; Rhobj = Rhobij * CONVERT_1; Kdj = Koc * Focj; Poraj = Porj - Thetawj; Darj=float(Dair*pow(Poraj,7./3.)/ (Porj*Porj)); Thetaj=Thetawj+Poraj*Kh+Rhobj*Kdj; Dispj = Aldisp[i]*Qj + Poraj*Darj*Kh; for(k=(J1[i][j]-1); k<J2[i][j]; ++k) { Rhob[i][k] = Rhobj; Kd[i][k] = Kdj; Por[i][k] = Porj; Pora[i][k] = Poraj; Foc[i][k] = Focj; Thetaw[i][k] = Thetawj; Dar[i][k] = Darj; Theta[i][k] = Thetaj; Disp[i][k] = Dispj; Q[i][k] = Qj; }//end for k }//end for j In >> Flxaqf[i] >> Cinaqfi[i] >> Pdepth[i] >> ch >> fl1 >> fl2; Cinaqf[i] = Cinaqfi[i] * CONVERT_2; Cinpaq[i] = Cinaqf[i]; Volfxc[i]=Q[i][Ncell[i]-1]*Areac[i]; }//end for i In.close(); }//end ReadInput //************************************** // OutOutput.cpp //************************************** #include "OutOutput.h" outClass::outClass(){ }//end default constructor outClass::outClass(char *FileName){ OpenFile(FileName); }//end constructor void outClass::OpenFile(char *FileName){ Out.open(FileName); }//end OpenFile void outClass::CloseFile(){ Out.close(); }//end CloseFile void outClass::SimTitle(char *Title){ Out << Title << endl; }//end SimTitle void outClass::Warn(float Time, float Iflag){ Out<< "WARNING!!! At time " << Time <<"aqueous solubility was exceeded in" << Iflag << " cells." << endl; }//end Warn void outClass::SoilTitle(char *Title){ Out << endl << endl; Out << Title << endl; }//end SoilTitle void outClass::AtTime(int SoilCol, float Time, float Mtotal, float Mgtot, float Mltot, float Mstot){ Out << endl << endl; Out << "Soil Column " << SoilCol << endl; Out << "At time = " << Time <<", total mass in vadose zone =" << Mtotal << " g/ft^2" << endl; Out << "Mass in gas phase = " << Mgtot << " g/ft^2" << endl; Out << "Mass in liquid phase = " << Mltot << " g/ft^2" << endl; Out << "Mass sorbed = " << Mstot << " g/ft^2" << endl; }//end AtTime void outClass::SinceBegin(float Time){ Out << endl << endl; Out << "Since beginning of run at time = " << Time << endl; }//end SinceBegin void outClass::SinceLast(float Time){ Out << endl << endl; Out <<"Since last printout at time =" << Time << endl; }//end SinceLast void outClass::Change(float Del1, float Mt, float Mb, float Mc, float Del2){ Out << " Change in Total mass = " << Del1 << " g/ft^2" << endl; Out << " Mass Flux Across Top Boundary =" << Mt << " g/ft^2" << endl; Out << " Mass Flux Across Bottom Boundary = " << Mb << " g/ft^2" << endl; Out << " Mass Loss due to Decay in Column = " << Mc << " g/ft^2" << endl; Out << " Net Mass Flux Through Boundaries = " << Del2 << " g/ft^2" << endl; Out << " Mass discrepancy = " << Del1 - Del2 << " g/ft^2" << endl; }//end Change void outClass::GWImpact(int SoilCol, int Ntime, float Ptime, float Gwimp[], float Areac){ Out << endl << endl; Out << "******************************************************************************" << endl; Out << "Ground-water IMPACT OF soil column" << SoilCol << endl; Out << endl; Out << "Time Mass flux (g/yr/ft^2) Total Mass (g/yr)" << endl; for(int i=0; i < Ntime; ++i) { Out << setprecision(3) << setw(9) << (i+1) * Ptime << setprecision(8) << setw(22) << Gwimp[i] / Areac << setprecision(5) << setw(31) << Gwimp[i] << endl;}//end for t Out << "******************************************************************************" << endl; //Spaces Out << endl << endl; }//end GWImpact void outClass::TotalGWTitle() { Out << "******************************************************************************" << endl; Out << "TOTAL GROUND WATER IMPACT" << endl; Out << "******************************************************************************" << endl; Out << endl; Out << "Time (yr) Mass (g/yr) Cumulative Mass (g)" << endl; }//end TotalGWTitle void outClass::TotalGWImpact(int Ntime, float Ptime, float Gwimtl, float Gcum){ Out << setprecision(3) << setw(8) << Ptime << setprecision(3) << setw(18) << Gwimtl << setprecision(3) << setw(24) << Gcum << endl; }//end TotalGWImpact //************************************** // PrfOutput.cpp //************************************** #include "PrfOutput.h" prfClass::prfClass() { }//end default constructor prfClass::prfClass(char *FileName){ OpenFile(FileName); }//end constructor void prfClass::OpenFile(char *FileName){ Out.open(FileName); }//end OpenFile void prfClass::CloseFile() { Out.close(); }//end CloseFile void prfClass::Title(char *Title, float Time) { Out << endl << endl; Out << Title << endl; Out << "Time : " << Time << endl; Out << "Dist(ft) Cgas(mg/L) Cliq(mg/L) Cs(mg/Kg)" << endl; }//end Title void prfClass::Output(float Depthz, float Cgas, float Cliq, float Cs) { Out << setprecision(3) << setw(9) << Depthz << setprecision(8) << setw(17) << Cgas << setprecision(8) << setw(18) << Cliq << setprecision(8) << setw(17) << Cs << endl; }//end Output void prfClass::Spaces() { Out << endl << endl; }//end Spaces //************************************** // PrmOutput.cpp //************************************** #include "PrmOutput.h" prmClass::prmClass() { }//end default condtructor prmClass::prmClass(char *FileName) { OpenFile(FileName); }//end condtructor void prmClass::OpenFile(char *FileName){ Out.open(FileName); }//end OutputFile void prmClass::CloseFile(){ Out.close(); }//end Close void prmClass::Output(vadoseInput &Vad){ int i,j,J1; //Space between each Out statement //represent different lines of Output Out << Vad.SimTitle << endl; Out << "There are " << Vad.Nscol << " soil Columns." << endl; Out<<"Timestep="<<Vad.Delt<<"years. Simulation length="<<Vad.Stime<<"years."<<endl; Out<<"Printout every"<<Vad.Ptime<<"years.Vertical Profile stored every"<< Vad.Prtime<<" years." << endl; Out << "Koc = " << Vad.Koci << " ml/g "; Out << Vad.Koc << " ft/g" << endl; Out << "Kh = " << Vad.Kh << " (dimensionless)" << endl; Out <<"Aqueous solulbility = " <<Vad.Cmaxi<<" ml/g, " <<Vad.Cmax<<" g/ft^3"<<endl; Out<<"Free air diffusion coeficient =" <<Vad.Dairi<<"cm^2/sec"<<Vad.Dair <<" ft^2/g"<<endl; for(i=0; i < Vad.Nscol; ++i){ //Spaces Out << endl << endl; Out << "SOIL COLUMN : " << (i+1) << endl << endl; Out << Vad.SoilTitle[i] << endl; Out << " Column Area = " << Vad.Areac[i] << " ft^2" << endl; Out << " "<<Vad.Ncell[i] <<"cells, each cell" <<Vad.Delz[i]<<"ft. thick." <<endl; Out<<"Longitudinal dispersivity (aL)in vadose zone=" <<Vad.Aldisp[i] <<"ft."<< endl; Out<<"Contaminant Decay Rate in the vadose zone ="<<Vad.Decy1[i]<<" 1/yr" << endl; Out << " Decay Rate in the Source = " << Vad.Decy2[i] << " 1/yr" << endl; //Spaces Out << endl << endl; for(j=0; j<Vad.GrNum[i]; ++j) { J1 = Vad.J1[i][j] - 1; Out<<"*Soil Properties*(Cell No." <<Vad.J1[i][j]<<"-" <<Vad.J2[i][j]<<")"<<endl; Out<<"Bulk density =" <<Vad.Rhob[i][J1]/CONVERT_1<<"g/ml," <<Vad.Rhob[i][J1] << " g/ft^3" << endl; Out<<"Porosity="<<Vad.Por[i][J1] <<"Volumetric water content=" << Vad.Thetaw[i][J1] << endl; Out << "Organic carbon content = " << Vad.Foc[i][J1] << endl; Out << "Recharge Rate = " << Vad.Q[i][J1] << " ft/yr" << endl; if (Vad.Typtbc[i] == 1) { Out << "Fixed Conc. at Cell #1 (1st"; }//end if Tybtc=1 else { Out << "Fixed Flux at Cell #1 (3rd"; }//end else Out<<"type b.c.)=" <<Vad.Valtbc[i] <<"mg/l,="<<Vad.Vatbc[i] <<"g/ft^3"<<endl; //Spaces Out << endl << endl; }//end for j Out << "**** Aquifer Properties ****" << endl; Out << "Horizontal GW Flux (q=ki) = " << Vad.Flxaqf[i] << " ft/yr" << endl; Out << "Background Conc. of GW Flux = " << Vad.Cinaqfi[i] << " mg/L" << endl; Out << "Leachate Penetration Depth = " << Vad.Pdepth[i] << " ft" << endl; }//end for i }//end Output //************************************** // SimOutput.cpp //************************************** #include "SimOutput.h" simClass::simClass() { }//end default constrructor simClass::simClass(char *FileName){ OpenFile(FileName); }//end constructor void simClass::OpenFile(char *FileName){ Out.open(FileName); }//end OpenFile void simClass::CloseFile() { Out.close(); }//CloseFile void simClass::Title(int SoilCol, float Time) { Out << "Soil Column " << SoilCol << " at t = " << Time << endl; Out << SoilCol << endl; }//end Title void simClass::Output(float Cs, float Delz) { Out << Cs << ", " << Delz << endl; }//end Output void simClass::Terminate() { Out << "9999" << endl; }//Terminate //************************************** // CmxOutput.cpp //************************************** #include "CmxOutput.h" cmxClass::cmxClass() { }//end default constructor cmxClass::cmxClass(char *FileName){ OpenFile(FileName); }//end constructor void cmxClass::OpenFile(char *FileName){ Out.open(FileName); }//end OpenFile void cmxClass::CloseFile() { Out.close(); }//end CloseFile void cmxClass::Title(char Angl) { Out << "SATURATED GW CONCENTRATION (underneath columns) DUE TO MIXING" << endl; if(Angl=='T')Out<<"SOIL COLUMNS ARE ARRANGED 90 DEGREE TO GROUNDWATER DIRECTION" << endl; else Out<<"SOIL COLUMNS ARE ARRANGED PARALLEL TO GROUNDWATER DIRECTION"<<endl; Out << " Liquid Phase Concentrations [mg/L] " << endl; }//end Title void cmxClass::ColInfo(char Angl, int SoilCol) { Out << endl << endl; Out << "Soil Column " << SoilCol << " (mg/L)" << endl; if (Angl == 'T') Out<<"TIME Cin(vadose) Cin(aquifer) Cout(mixed)" << endl; Else { if (SoilCol == 1) Out << " TIME Cw,vad(in) Caq,aqui(in) Cmx(out)" << endl; else Out << " TIME Cw,vad(in) Caq,aqui(in) Cmx(in) Cmx(out)" << endl; }//end else }//end ColINfo void cmxClass::RAngleOutput(float Time, float Cliq, float Cinaqf, float Cmixed){ Out.setf(ios::dec); Out << setiosflags(ios::dec) << setprecision(3) << setw(8) << Time; Out << setiosflags(ios::scientific) << setw(17) << Cliq << setiosflags(ios::scientific) << setw(20) << Cinaqf << setiosflags(ios::scientific) << setw(19) << Cmixed << endl; }//end RAngleOutput void cmxClass::PAngleOutput_1(float Time, float Cliqend, float Cinpaq, float Cmixe) { Out.unsetf(Flag); Out << setiosflags(ios::dec) << setprecision(3) << setw(8) << Time << setiosflags(ios::scientific) << setw(18) << Cliqend << setiosflags(ios::scientific) << setw(20) << Cinpaq << setiosflags(ios::scientific) << setw(16) << Cmixe << endl; }//end PAngleOutput_1 void cmxClass::PAngleOutput_2(float Time, float Cliqend, float Cinpaq, float CmixeIn, float CmixeOut) { Out.unsetf(Flag); Out << setiosflags(ios::dec) << setprecision(3) << setw(8) << Time << setiosflags(ios::scientific) << setw(18) << Cliqend << setiosflags(ios::scientific) << setw(20) << Cinpaq << setiosflags(ios::scientific) << setw(15) << CmixeIn << setiosflags(ios::scientific) << setw(16) << CmixeOut << endl; }//end PAngleOutput_2 //************************************** // GimOutput.cpp //************************************** #include "GimOutput.h" gimClass::gimClass() { }//end default constructor gimClass::gimClass(char *FileName) { OpenFile(FileName); }//end constructor void gimClass::OpenFile(char *FileName) { Out.open(FileName); }//end OpenFile void gimClass::CloseFile() { Out.close(); }//end CloseFile void gimClass::Title(int SoilCol) { Out << "Soil Column " << SoilCol+1 << endl; }//end Title void gimClass::NumOutput(float Time, float Num) { Out << Time << ", " << Num << endl; }//end NumOutput void gimClass::Terminate() {
Last modified: Oct 15, 1999 VG Model / Samuel Lee / VADOSE.NET |