INCLUDE "ChemPlugin.f90"

PROGRAM ReactiveTransport
    Use ChemPluginModule
    IMPLICIT NONE
    TYPE(ChemPlugin) :: cp_inlet
    TYPE(ChemPlugin),Dimension(:),Allocatable :: cp
    TYPE(CpiLink) :: link0
    INTEGER :: file_unit=20, ios=0, error, nx, i, scope = 0
    REAL(8) :: time_end, delta_years, next_output, &
               length, deltax, deltay, deltaz, &
               veloc_in, velocity, flow, porosity, &
               diffcoef, dispersivity, dispcoef, trans,&
               tcond, ttrans, dt, deltat
    CHARACTER(LEN=500) :: input_file, line, cmd

    WRITE(*,'("Model reactive transport in one dimension")')
    WRITE(*,'("")')
    nx = 100 ! number of instances along x
    length = 100 ! m
    deltax = length / nx ! m
    deltay = 1.0  ! m
    deltaz = 1.0  ! m
    time_end = 10.0  ! years
    delta_years = time_end / 5 ! years
    next_output = 0

    ! Create the inlet and interior instances
    cp_inlet = ChemPlugin()
    
    ALLOCATE(cp(nx))
     
    Do i = 1, nx
       cp(i) = ChemPlugin()
    END DO
    error = cp_inlet%Console("stdout")
    error = cp(1)%Console("stdout")
    error = cp(1)%Config("pluses = banner")

    WRITE(cmd, '("volume =", F6.2, " m3; ", "time end =", F5.1, " years;")') deltax * deltay * deltaz, time_end
    Do i = 1, nx
       error = cp(i)%Config(cmd)
    END DO
    WRITE(*,'( "Enter RTM input script: ")', ADVANCE="no")
    READ(*,'(A)') input_file

    !Open output file and write instance positions on the first line
    OPEN(unit=file_unit, file=input_file, action="read")
    Do WHILE(ios == 0) 
       READ(file_unit,'(A)',IOSTAT=ios) line
       IF (line == 'go') THEN
           EXIT
       ELSE IF (line == 'scope inlet') THEN
           scope = 1
       ELSE IF (line == 'scope initial') THEN
           scope = 2
       ELSE IF (scope == 1) THEN
           error = cp_inlet%Config(line)
       ELSE IF (scope == 2) THEN
           DO i = 1, nx
              error = cp(i)%Config(line)
           END DO
       END IF
    END DO
  
    ! Initialize the inlet and interior instances; write out initial conditions
    IF (cp_inlet%Initialize() /= 0) THEN
        WRITE(*,*) "Inlet failed to initialize"
        STOP
    END IF
    error = cp_inlet%PrintOutput("RTM.txt", "Inlet fluid")

    DO i = 1, nx
       IF (cp(i)%Initialize() /= 0) THEN
          WRITE(*, '("Instance", I2, " failed to initalize")') i
          STOP
       END IF
    END DO

    call write_results(cp, nx, delta_years, next_output)

    ! Query user for velocity; calculate flow rate and transmissivities
    WRITE(*,'(A)', Advance='no') "Please enter fluid velocity in m/yr: "
    READ(*,*) veloc_in
    WRITE(*,'()', Advance='yes')
    velocity = veloc_in / 31557600 ! m/s

    porosity = cp(1)%Report1("porosity")
    flow = deltay * deltaz * porosity * velocity ! m3/s

    diffcoef = 1d-10 ! m2/s
    dispersivity = 1.0 ! m
    dispcoef = velocity * dispersivity + diffcoef ! m2/s
    trans = deltay * deltaz * porosity * dispcoef / deltax  ! m3/s
    
    tcond = 2.0 ! W/m/K
    ttrans = deltay * deltaz * tcond / deltax

    ! Link the instances
    link0 = cp(1)%Link(cp_inlet)
    error = link0%Transmissivity(trans, "m3/s")
    error = link0%HeatTrans(ttrans, "W/K")
    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%HeatTrans(ttrans, "W/K")
       error = link0%FlowRate(flow, "m3/s")
    END DO
    link0 = cp(nx)%Link()
    error = link0%FlowRate(-flow, "m3/s")

    DO WHILE(.TRUE.)
       deltat = 1e10
       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
              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)%AdvanceHeatTransport() /= 0) THEN
              STOP
          END IF
       END DO  
       DO i = 1, nx
          IF (cp(i)%AdvanceChemical() /= 0) THEN
              STOP
          END IF
       END DO  
       call write_results(cp, nx, delta_years, next_output)
    END DO
        
CONTAINS 
        SUBROUTINE write_results(cp, nx, gap, next_output) 
            IMPLICIT NONE 
            TYPE(ChemPlugin),dimension(:), INTENT(in) :: cp
            INTEGER, INTENT(in) :: nx
            REAL(8), INTENT(in) :: gap
            REAL(8), INTENT(inout) :: next_output
            CHARACTER(LEN=11) :: label
            INTEGER :: i , error
            REAL(8) :: now

            now = cp(1)%Report1("Time", "years")
            
            IF ((next_output - now) <= gap/1d4) THEN
                Do i = 1, nx
                   WRITE(label,'("Instance ",I2)') i-1
                   error = cp(i)%PrintOutput("RTM.txt", label)
                END DO
                next_output = next_output + gap
            END IF
         
        END SUBROUTINE write_results

END PROGRAM ReactiveTransport