Appendix A
FORTRAN Source Code for the Constantly Stirred Tank Reduction Model

The following FORTRAN code defines the Constantly Stirred Tank Reduction (CSTR) Model used in the example of how to integrate a model into FRAMES.

!     Last change:  KJC  21 Jun 2003   10:50 pm
program main
implicit none
! The following lines represent data associate with each chemical.
  CHARACTER(80) CNAME,CASID
  REAL Qco,kd,S,lambda,Ctp,TFinal,Ct
! QCo (g/yr)   Mass flux into CSTR
! kd (cm^3/g)  Soil partition coefficient
! S (mg/L)     Solubility limit
! lamda (1/yr) Decay rate of the chemical
! CTp (mg/kg)  Initial concentration on total weight
! CT  (g/cm^3) Initial concentration on total volume
! TFinal (yr)  Last year to report
! EnterCT is a flag to signify whether CT or CTp is entered.
  INTEGER NTimes
! The number of times to report


! The following line represents data associate with the site
  REAL L,W,Th,b,Ks,i,thetaFc,rhoS,Bd
! L (cm)         Length of site
! W (cm)         Width of site
! Th (cm)        Thickness of site
! b ()           Soil type exponent
! Ks (cm/s converted to cm/yr) Staturated K
! i (cm/yr)      Infiltration rate
! thetaFc ()     Field capacity
! rhoS (g/cm^3)  Particle density
! Bd (g/cm^3)    Bulk density

! The following line represents data calculated for a site.
REAL Vt,n,m,thetaCalc,theta,Q
! The following line represents data calculated for a site given a chemical.
REAL Rf,HalfLife,Ci,Mass,Lambda2,Tau1

! The following line is associated with creating the results and reading
! multiple chemicals.
INTEGER NumChem,Chem,Flag

! The following lines are associated with calculating each time.
REAL dt,arg,C,t,Vw,D,Vs,Tfp,Mfp,Fmc
LOGICAL SolLimited
INTEGER AddNTimes,ioerr

! Open files and read site information.
  OPEN(UNIT=20,  FILE="CSTRInput.txt",STATUS='OLD')
  OPEN(UNIT=21,  FILE="CSTROutput.txt")
  OPEN(UNIT=23,  FILE="CSTRMsgs.txt")
  WRITE(23,*) "Settings"
  READ(20,*)L,W,Th,b,Ks,i,thetaFc,rhoS,Bd
  WRITE(23,*)L,W,Th,b,Ks,i,thetaFc,rhoS,Bd
! Calculate additional site parameters.
  Vt=L*W*Th;
  n=1.0-Bd/rhoS
  m=2.0*b+3.0
  thetaCalc=n*(i/Ks)**(1/m)
  if (thetaCalc.gt.thetaFc) then
    theta=thetaCalc
  else
    theta=thetaFC
  end if
  Q=i*L*W
  WRITE(23,*) "Calculated Values"
  WRITE(23,*) Vt,n,m,thetaCalc,theta,Q
  READ(20,*)NumChem
  WRITE(21,*)NumChem
  WRITE(23,*)NumChem
  do chem=1,NumChem
    WRITE(23,*)"Chemical specific"
    READ(20,*) CNAME
    WRITE(23,*) CNAME
    READ(20,*) CASID
    write(23,*) CASID
    READ(20,*)Qco,kd,S,lambda,Tfinal
    READ(20,*)NTimes
!   Convert S.
    S=S/(1000.0*1000.0)
    write(23,*)Qco,kd,S,lambda,Tfinal
    WRITE(23,*)NTimes
    READ(20,*)flag
    write(23,*)flag
    if (flag.ne.0) then
      READ(20,*)Ct
      write(23,*)Ct
    else
      READ(20,*)Ctp
      write(23,*)Ctp
    end if
    if (lambda.lt.1.0E-20) then
      HalfLife=1.0E+20
    else
      HalfLife=-log(0.5)/lambda
    end if
!   Calculate additional chemical and site combined parameters.
    Rf=1.0+Bd*kd/theta
    if (flag.ne.0) then
      Ct=Ctp*1.0E-6*Bd
    end if
    Ci=Ct/(theta+Bd*kd)
    Mass=Ct*Vt
    lambda2=Qco/(Vt*theta)
    Tau1=Vt*theta/Q
    WRITE(23,*)HalfLife,Rf,Ct,Ci,Mass,Lambda2,Tau1
!   Calculate the time series.
    t=0.0
    Vs=Vt*(1.0-n)
    Vw=Vt*theta
    D=Mass/(kd*rhoS*Vs+Vw)
    SolLimited=S.lt.D
    WRITE(23,*) "Volumes"
    WRITE(23,*) Vs,Vw,D
    if (SolLimited) then
      WRITE(23,*) "Solubility limited"
      Mfp=Mass-s*(kd*rhos*Vs+Vw)
      Fmc=Q*S
      Tfp=Mfp/Fmc
      AddNTimes=1
      Ci=S
      WRITE(23,*)Mfp,Fmc,Tfp,Ci
    else
      Tfp=0.0
      AddNTimes=0
    end if
    WRITE(21,*)CName
    WRITE(21,*)CASID
    WRITE(21,*)NTimes+AddNTimes
    WRITE(21,*)"yr"," ","g/ml"," ","g/yr"
    dt=TFinal/NTimes
    arg=((1.0+Tau1*Rf*Lambda)/(Rf*Tau1))
    if (SolLimited) then
      WRITE(21,FMT='(3g8.3)')0.0,S,Fmc
    end if
    do WHILE (t+dt.le.TFinal)
      c=Tau1*lambda2/(1.0+Tau1*Rf*lambda)+(Ci-(Tau1*lambda2/(1.0+Tau1*Rf*lambda)))**(-arg*t)
      WRITE(21,FMT='(3g8.3)')t+tfp,C,Q*C
      t=t+dt
    end do
  end do
  CLOSE(UNIT=20)
  CLOSE(UNIT=21)
  CLOSE(UNIT=23)
end program

Home | Security and Privacy | Contact Us