INCLUDE "ChemPlugin.f90"

PROGRAM Advection
    Use ChemPluginModule
    IMPLICIT NONE
    TYPE(ChemPlugin) :: cp_inlet
    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, &
            deltat, dt
    CHARACTER(LEN=255) :: cmd1
    CHARACTER(:), allocatable :: cmd
    CHARACTER :: tab=char(9)
 
     
    WRITE(*,fmt=('("Model advection-dispersion in one dimension")'))
    WRITE(*,fmt=('("")'))
   
    ! Simulation parameters
    nx = 400 ! number of instances along x
    length = 100 ! m
    deltax = length / nx ! m
    deltay = 1.0
    deltaz = 1.0
    porosity = 0.25

    WRITE(*,'("Please enter fluid velocity in m/yr: ")',advance="no") 
    READ(*,*) veloc_in
    velocity = veloc_in / 31557600 ! m/s
    flow = deltay * deltaz * porosity * velocity ! m3/s

    diffcoef = 1e-10 ! m2/s
    dispersivity = 1.0 ! m
    dispcoef = velocity * dispersivity + diffcoef ! m2/s
    trans = deltay * deltaz * porosity * dispcoef / deltax ! m3/s

    time_end = 10.0 ! years
    delta_years = time_end / 5.0 ! years
    next_output = 0.0
     
    !Open output file and write instance positions on the first line
    OPEN(unit=file_unit, file="Advection.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 inlet instance
    cp_inlet = ChemPlugin()
    error = cp_inlet%Config("Na+= 1 mmol/kg; Cl- = 1 mmol/kg")
    If (cp_inlet%Initialize() /= 0) THEN
        WRITE(*,*) "Inlet failed to initialize"
    END IF
    error = cp_inlet%Console("stdout")

   ! Configure and initialize interior instances
    ALLOCATE(cp(nx))
    cmd = "Na+ = 0.001 mmol/kg; Cl- = 0.001 mmol/kg;"
    WRITE(cmd1,'("volume =", G10.3, " m3;")') deltax * deltay * deltaz
    cmd = trim(cmd)//trim(cmd1)
    WRITE(cmd1,'("porosity =", G10.3, ";")') porosity
    cmd = trim(cmd)//trim(cmd1)
    WRITE(cmd1, '("time end =", G10.1, " year;")') time_end
    cmd = trim(cmd)//trim(cmd1)

    cp(1) = ChemPlugin()
    error = cp(1)%Console("stdout")
    error = cp(1)%Config("pluses = banner")
    error = cp(1)%Config(cmd)
    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(cmd)
       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
    link0 = cp(1)%Link(cp_inlet)
    error = link0%Transmissivity(trans, "m3/s")
    error = link0%FlowRate(flow, "m3/s")

    DO i = 2, nx
       link0 = cp(i)%Link(cp(i-1))
       error = link0%Transmissivity(trans, "m3/s")
       error = link0%FlowRate(flow, "m3/s")
    END DO

    ! Outlet link
    link0 = cp(nx)%Link()
    error = link0%FlowRate(-flow, "m3/s")

    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 Advection