INCLUDE "ChemPlugin.f90"

PROGRAM FlowThrough
    Use ChemPluginModule
    IMPLICIT NONE
    TYPE(ChemPlugin) :: cp_inlet, cp
    TYPE(CpiLink) :: link1, link2
    INTEGER :: error, i
    REAL(8) :: deltat
    CHARACTER(LEN=255),dimension(3) :: &
    cmds = [character(len=255) :: &
            "swap Calcite for Ca++; Calcite = 0.03 free m3; Cl- = 1 mmol/kg", &
            "swap CO2(g) for H+; fugacity CO2(g) = 1; balance on HCO3-", &
            "volume = 1 m3; fix f CO2(g); delxi = 0.01; pluses = banner"]
 
    WRITE(*,fmt=('("Model a flow-through reactor")'))
    WRITE(*,fmt=('("")'))

    ! Configure and initialize the inlet fluid
    cp_inlet = ChemPlugin("stdout")
    error = cp_inlet%Config("Ca++ = 1 mmol/kg; HCO3- = 1 mmol/kg;pH = 1; balance on Cl-")
    error = cp_inlet%Initialize()

    ! configure and initialize the stirred reactor
  

    cp = ChemPlugin("stdout")
    DO i = 1, 3
       error = cp%Config(cmds(i))
    END DO
    error = cp%Initialize(1.0d+0, "day")

    ! Link reactor to inlet and free outlet; set rate of flow.
    link1 = cp%Link(cp_inlet)
    link2 = cp%Outlet()
    error = link1%FlowRate(10.0d+0, "m3/day")
    error = link2%FlowRate(-10.0d+0, "m3/day")

    error = cp%PlotHeader("FlowThrough.gtp", "char")
    error = cp%PlotBlock()

    DO WHILE(.True.)
        deltat = cp%ReportTimeStep() 
        IF (cp%AdvanceTimeStep(deltat) /= 0) THEN
            EXIT
        END IF
        IF (cp%AdvanceTransport() /= 0) THEN
            EXIT
        END IF
        IF (cp%AdvanceChemical() /= 0) THEN
            EXIT
        END IF
        error = cp%PlotBlock() 
    END DO
    
END PROGRAM FlowThrough