INCLUDE "ChemPlugin.f90"

PROGRAM Diffusion
    Use ChemPluginModule
    IMPLICIT NONE
    TYPE(ChemPlugin),dimension(:), allocatable :: cp
    TYPE(CpiLink) :: link0
    INTEGER :: nx, error, i,j, file_unit=20
    REAL(8) :: length, deltax, deltay, deltaz, &
            porosity, veloc_in, velocity, flow, &
            diffcoef, dispersivity, dispcoef, trans, &
            time_end, delta_years, next_output, &
            xstable, deltat, dt
    CHARACTER(LEN=255) :: cmd, cmd1, cmd2
    CHARACTER(:), ALLOCATABLE :: cmd0
    CHARACTER :: tab=char(9)
     
    WRITE(*,fmt=('("Model diffusion in one dimension")'))
    WRITE(*,fmt=('("")'))
   
    ! Simulation parameters
    nx = 100 ! number of instances along x
    length = 100 ! cm
    deltax = length / nx ! cm
    deltay = 1.0  ! cm
    deltaz = 1.0  ! cm
    porosity = 0.25 ! volume fraction

    diffcoef = 1e-6 ! cm2/s
    trans = deltay * deltaz * porosity * diffcoef / deltax ! cm3/s

    xstable = 0.9
    time_end = 15.0 ! years
    delta_years = time_end / 3.0 ! years
    next_output = 0.0 ! years

     
    !Open output file and write instance positions on the first line
    open(unit=file_unit, file="Diffusion.txt", action='write')
    WRITE(file_unit,'("years")',advance="no") 
    Do i = 1, nx-1
       WRITE(file_unit,fmt='(a,F8.4)',advance="no") tab, (i - 0.5) * deltax  
    END DO
    WRITE(file_unit,fmt='(a,F8.4)') tab, (nx - 0.5) * deltax 
   

    ! Configure and initialize the instances
    
    WRITE(cmd, '("volume =", F4.1, " cm3;")') deltax * deltay * deltaz
    cmd0 = cmd
    WRITE(cmd, '("porosity =", F5.2, ";")') porosity
    cmd0 = trim(cmd0)//trim(cmd)
    WRITE(cmd, '("time end =", F5.1, " years;")') time_end
    cmd0 = trim(cmd0)//trim(cmd)
    WRITE(cmd, '("Xstable =", F5.2)') xstable
    cmd0 = trim(cmd0)//trim(cmd)  
    cmd1 = "Na+ = 1 mmol/kg; Cl- = 1 mmol/kg"
    cmd2 = "Na+ = 0.001 mmol/kg; Cl- = 0.001 mmol/kg"

    ALLOCATE(cp(nx))
    cp(1) = ChemPlugin()
    error = cp(1)%Console("stdout")
    error = cp(1)%Config("pluses = banner")
    error = cp(1)%Config(cmd0)
    error = cp(1)%Config(cmd1)
    IF (cp(1)%Initialize() /= 0) THEN
        WRITE(*,*) "Instance 1 failed to initialize"
    END IF
    DO i = 2, nx
       cp(i) = ChemPlugin()
       error = cp(i)%Config(cmd0)
       IF (i < nx/2 + 1) THEN 
          error = cp(i)%Config(cmd1)
       ELSE
          error = cp(i)%Config(cmd2)
       END IF
       IF (cp(i)%Initialize() /= 0) THEN
           WRITE(*,'("Instance ", I1, "failed to initialize")') i
       END IF
    END DO

    call write_line(file_unit, cp, nx, delta_years, next_output)
    
    ! Link the instances
    DO i = 2, nx
       link0 = cp(i)%Link(cp(i-1))
       error = link0%Transmissivity(trans, "cm3/s")
    END DO

    DO WHILE(.TRUE.) 
       deltat = 1e9
       DO i = 1, nx
          dt = cp(i)%ReportTimeStep()
          deltat = min(deltat, dt)
       END DO
       DO i = 1, nx
          IF (cp(i)%AdvanceTimeStep(deltat) /= 0) THEN
              READ(*,*)
              STOP
          END IF
       END DO
       DO i = 1, nx
          IF (cp(i)%AdvanceTransport() /= 0) THEN
              STOP
          END IF
       END DO   
       DO i = 1, nx
          IF (cp(i)%AdvanceChemical() /= 0) THEN
              STOP
          END IF
       END DO
       call write_line(file_unit, cp, nx, delta_years, next_output)

    END DO 
       

    CONTAINS 
        SUBROUTINE write_line(file_unit, cp, nx, gap, next_output) 
            IMPLICIT NONE 
            INTEGER, INTENT(in) :: file_unit
            TYPE(ChemPlugin),dimension(:), INTENT(in) :: cp
            INTEGER, INTENT(in) :: nx
            REAL(8), INTENT(in) :: gap
            REAL(8), INTENT(inout) :: next_output
            INTEGER :: i , error
            REAL(8) :: now(1), conc(1)
            !REAL ::  now(1)

            error = cp(1)%Report(now, "Time", "years")
            
            IF ((next_output - now(1)) <= gap/1d4) THEN
                WRITE(file_unit,fmt='(F8.4)',advance="no") now
                Do i = 1, nx-1
                   error = cp(i)%Report(conc, "concentration Na+", "mmol/kg")
                   WRITE(file_unit,fmt='(a,F8.4)',advance="no") tab, conc 
                END DO
                error = cp(nx)%Report(conc, "concentration Na+", "mmol/kg")
                WRITE(file_unit,fmt='(a,F8.4)') tab, conc 
                next_output = next_output + gap
            END IF
         
        END SUBROUTINE write_line
END PROGRAM Diffusion