INCLUDE "ChemPlugin.f90"

PROGRAM titration
    Use ChemPluginModule
    IMPLICIT NONE
    TYPE(ChemPlugin) :: cp
    INTEGER :: error, i, nsp
    REAL(8) :: deltat
    CHARACTER(LEN=255) :: spec
    CHARACTER(LEN=100) :: spec_keywords, conc_keywords 
    REAL(8) :: 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%Report1i("nspecies")
    
    DO WHILE(.True.)
        deltat = cp%ReportTimeStep() 
        IF (cp%AdvanceTimeStep(deltat) /= 0) THEN
            EXIT
        END IF
        IF (cp%AdvanceChemical() /= 0) then
            EXIT
        END IF 
        
        DO i = 1, nsp
           WRITE(spec_keywords, '("species", I3)') i-1
           WRITE(conc_keywords, '("concentration aqueous", I3)') i-1
           spec = cp%Report1c(spec_keywords)
           conc = cp%Report1(conc_keywords, "umol/kg")
           IF (conc >= 10.0) THEN
              WRITE(*, '(A10," = ", F12.2, " umol/kg")') spec, conc
           END IF
        END DO
        WRITE(*,'("")')
    END DO
END Program titration