INCLUDE "ChemPlugin.f90" PROGRAM HeatConduction 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, trans, tcond, & time_end, delta_years, next_output, & deltat, dt CHARACTER(LEN=255) :: cmd0 CHARACTER(:), ALLOCATABLE :: cmd CHARACTER :: tab=char(9) WRITE(*,fmt=('("Model heat conduction in one dimension")')) WRITE(*,fmt=('("")')) ! Simulation parameters nx = 100 ! 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 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="HeatConduction.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 cmd = "span 20 C to 100 C; Na+ = 0.001 mmol/kg; Cl- = 0.001 mmol/kg;" WRITE(cmd0, '("volume =", F4.1, " 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) error = cp(1)%Config("T = 100 C") 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 (i < nx/2 + 1) THEN error = cp(i)%Config("T = 100 C") ELSE error = cp(i)%Config("T = 20 C") 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%HeatTrans(trans, "W/K") 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 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 HeatConduction