A.1: Codes for Vadose Zone Leaching and Saturated Zone Mixing
A.2 Codes of Groundwater Flow











 

Appendix A. SOURCE CODES FOR VG MODEL


A.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()        {

  Out << "9999" << endl; }//end Terminate

 


 



Last modified: Oct 15, 1999
VG Model / Samuel Lee / VADOSE.NET