INCLUDE "ChemPlugin.f90"

PROGRAM titration
    Use ChemPluginModule
    IMPLICIT NONE
    TYPE(ChemPlugin) :: cp
    INTEGER :: error
    REAL(8) :: deltat
    

    WRITE(*,'("ChemPlugin example -- pH titration")')
    WRITE(*,'("")')

    ! Create ChemPlugin instance
    cp = ChemPlugin("stdout")

    ! Configure the instance
    error = cp%Config("Ca++ = 1 mmol/kg; Na+ = 1 mmol/kg")
    error = cp%Config("Cl- = 3 mmol/kg; HCO3- = 2 mmol/kg; pH = 4")
    error = cp%Config("react 3 mmol/kg NaOH; delxi = 0.1")
    
    ! Initialize the instance
    error = cp%Initialize()
    error = cp%PlotHeader("myPlot.gtp", "char")
    error = cp%PlotBlock()
    
    DO WHILE(.True.)
        deltat = cp%ReportTimeStep() 
        IF (cp%AdvanceTimeStep(deltat) /= 0) THEN
            EXIT
        END IF
        IF (cp%AdvanceChemical() /= 0) then
            EXIT
        END IF 
        error = cp%PlotBlock()  
    END DO
    error = cp%PlotTrailer() ! Not neccessary

END Program titration