INCLUDE "ChemPlugin.f90"

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

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

    cp = ChemPlugin("stdout")
    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")

    error = cp%Initialize()
    if (error /= 0) then
        Write(*,*) "Instance failed to initialize"
    end if

    error = cp%Report(Xi, "Xi")
    error = cp%Report(pH, "pH")
    WRITE(*,'(" Xi = ", F4.2, " pH = ", F4.2)') Xi, pH

    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%Report(Xi, "Xi")
        error = cp%Report(pH, "pH")
        WRITE(*,'(" Xi = ", F4.2, " pH = ", F4.2)') Xi, pH
    END DO
    
END Program titration