INCLUDE "ChemPlugin.f90"

PROGRAM titration
    Use ChemPluginModule
    IMPLICIT NONE
    TYPE(ChemPlugin) :: cp
    INTEGER :: error, i, nsp
    REAL(8) :: deltat
    CHARACTER(LEN=255), DIMENSION(:), ALLOCATABLE :: spec
    REAL(8), Dimension(:), ALLOCATABLE :: conc
    

    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()

    nsp = cp%Report("species")
    Allocate(spec(nsp))
    Allocate(conc(nsp))
    error = cp%Report(spec, "species")

    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(conc, "concentration aqueous", "umol/kg")
        DO i = 1, nsp
           IF (conc(i) >= 10.0) THEN
              WRITE(*, '(A10," = ", F12.2, " umol/kg")') spec(i), conc(i)
           END IF
        END DO
        WRITE(*,'("")')
    END DO
END Program titration