INCLUDE "ChemPlugin.f90"

PROGRAM HeatTransfer
    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, trans, tcond, &
            veloc_in, velocity, flow, &
            time_end, delta_years, next_output, &
            deltat, dt
    CHARACTER(LEN=255) :: cmd0
    CHARACTER(:), ALLOCATABLE :: cmd
    CHARACTER :: tab=char(9)
     
    WRITE(*,fmt=('("Model heat transfer in one dimension")'))
    WRITE(*,fmt=('("")'))
   
    ! Simulation parameters
    nx = 400 ! number of instances along x
    length = 100 ! m
    deltax = length / nx ! m
    deltay = 1.0  ! m
    deltaz = 1.0  ! m
    porosity = 0.25 ! volume fraction

    tcond = 2.0 ! W/m/K
    trans = deltay * deltaz * tcond / deltax ! W/K

    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
    
    time_end = 10.0 ! years
    delta_years = time_end / 5.0 ! years
    next_output = 0.0 ! years

    !Open output file and write instance positions on the first line
    open(unit=file_unit, file="HeatTransfer.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 initliaze the inlet instance
    cp_inlet = ChemPlugin()
    error = cp_inlet%Config("T = 100 C; Na+ = 0.001 mmol/kg; Cl- = 0.001 mmol/kg")
    If (cp_inlet%Initialize() /= 0) THEN
        WRITE(*,*) "Inlet failed to initialize"
        STOP
    END IF

    ! Configure and initialize the interior instances
    cmd = "T = 20 C; span 20 C to 100 C; Na+ = 0.001 mmol/kg; Cl- = 0.001 mmol/kg;"
    WRITE(cmd0, '("volume =", F5.2, " m3;")') deltax * deltay * deltaz
    cmd = trim(cmd)//trim(cmd0)
    WRITE(cmd0, '("porosity =", F5.2, ";")') porosity
    cmd = trim(cmd)//trim(cmd0)
    WRITE(cmd0, '("time end =", F5.1, " years;")') time_end
    cmd = trim(cmd)//trim(cmd0)
  
    ALLOCATE(cp(nx))
    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(*,*) "Interior Instance 1 failed to initialize"
        STOP
    END IF
    DO i = 2, nx
       cp(i) = ChemPlugin()
       error = cp(i)%Config(cmd)
       IF (cp(i)%Initialize() /= 0) THEN
           WRITE(*,'("Interior Instance ", I1, "failed to initialize")') i
           STOP
       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%HeatTrans(trans, "W/K")
    error = link0%FlowRate(flow, "m3/s")
    DO i = 2, nx
       link0 = cp(i)%Link(cp(i-1))
       error = link0%HeatTrans(trans, "W/K")
       error = link0%FlowRate(flow, "m3/s")
    END DO
    link0 = cp(nx)%Link() ! Outlet
    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)%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_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), temp(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(temp, "temperature", "C")
                   WRITE(file_unit,fmt='(a,F8.4)',advance="no") tab, temp 
                END DO
                error = cp(nx)%Report(temp, "temperature", "C")
                WRITE(file_unit,fmt='(a,F8.4)') tab, temp 
                next_output = next_output + gap
            END IF
         
        END SUBROUTINE write_line

END PROGRAM HeatTransfer