!= tubreac.f90 12Apr2009 20060521 Tubular reactor, from Robinson
! ... Keywords: tubular reaction simulation Robinson
 CALL timer (0, ict); CALL mainsub; CALL timer (1, ict); STOP
 END !END!

 SUBROUTINE mainsub
! ... ..................................................................
 IMPLICIT DOUBLE PRECISION (a-h, o-z)
 ALLOCATABLE c(:), dc(:), tvec(:), z(:,:)
 DATA n/ 3 / ! components
 WRITE (*, "(' *** Tubular reaction (E. R. Robinsons) ***', T66, " // &
   "'(May-2006, MC)' /)")
 ALLOCATE (c(n), dc(n))
 CALL readat (n, r, zk1inf, zk2inf, e1, e2, c, temp, timef, &
   deltat, dt, nprint, nt, ishow)
 zk1 = zk1inf * EXP(-e1 / (r * temp)); zk2 = zk2inf * EXP(-e2 / (r * temp))
 WRITE (*, 99300) (i, i=1, n)
99300 FORMAT (/ 7X, 'Time', 2X, 5(5X, 'Conc(', I1, ')':))
 time = 0
 WRITE (*, 99500) time, c
99500 FORMAT (1X, 6F12.5)
 ALLOCATE (tvec(0:nprint), z(n,0:nprint))
! ... Perform integration by simple Euler
 tvec(0) = time; z(:,0) = c
 DO jp=1, nprint; jt0 = (jp - 1) * nt
   DO jt=1, nt; time = (jt0 + jt) * dt
      CALL eulerkinet (n, c, zk1, zk2, dt, dc)
   END DO
   tvec(jp) = time; z(:,jp) = c; WRITE (*, 99500) time, c
 END DO
 time = 10; jp = MINLOC(ABS(tvec-time), DIM=1) - 1
 WRITE (*, 99800) time, z(2,jp)
99800 FORMAT (' c(2) at t =', F5.1, ':', T26, F11.4, 3X, '(in graph)')
! ... Plot
 CALL www_graph (ishow, 0, nprint, n, tvec, z)
! ... Ref.: Robinson, E. R., 1975, "Time Dependent Chemical Processes", Applied
!   . Science Publishers, Ltd., London (UK), Programme 7.3, p 137, modified
 RETURN
 END SUBROUTINE !END!

 INCLUDE './w3g.gnuParam.f90'
 SUBROUTINE eulerkinet (n, c, zk1, zk2, dt, dc)
! ... Euler method and kinetic equations................................
!   . Bilous, O., and N. Amundson, 1956, "Optimum temperature gradients 
!   . in tubular reactors", Chem. Enging. Sci., vol. 5, 81 & 115
!   . Google: +"Chemical Engineering Science" +Journal +Elsevier
 IMPLICIT DOUBLE PRECISION (a-h, o-z)
 DIMENSION c(n), dc(n)
 dc(1) = - zk1 * c(1); dc(2) = zk1 * c(1) - zk2 * c(2); dc(3) = zk2 * c(2)
 c = c + dc * dt
 RETURN
 END SUBROUTINE !END!

 SUBROUTINE readat (n, rgas, zk1inf, zk2inf, e1, e2, c, temp, &
   timef, deltat, dt, nprint, nt, ishow)
! ... ..................................................................
 IMPLICIT DOUBLE PRECISION (a-h, o-z)
 DIMENSION c(n)
 rgas = 2; READ (*, *) zk1inf, zk2inf, e1, e2
 WRITE (*, 95100) rgas, zk1inf, zk2inf, e1, e2
95100 FORMAT (' Reactor conditions and data:' / &
   3X, 'R (gas constant),', T30, 1P, G12.5, 2X, 'kcal / kmol-K' / &
   3X, 'k1, k2 _inf,', T18, 1P, 2G12.5, 2X, 'min^-1', T65, '[t] = min' / &
   3X, 'E1, E2,', T18, 1P, 2G12.5, 2X, 'kcal / kmol')
 READ (*, *) c; WRITE (*, 95200) c
95200 FORMAT (3X, 'c,', 1P, (T18, 5G12.5))
 READ (*, *) temp; READ (*, *) timef, deltat, dt; READ (*, *) ishow
 WRITE (*, 95300) temp, timef, deltat, dt
95300 FORMAT (3X, 'Temperature,', T30, G12.5, 2X, 'K' / &
   3X, 'Time_f,', T30, G12.5, 2X, 'min', T50, 'Simple Euler integration' / &
   3X, 'Print step, dt,', T22, 2G10.3, 2X, 'min')
 nprint = NINT(timef / deltat ); nt = NINT(timef / dt / nprint)
 WRITE (*, 95600) nprint, nt
95600 FORMAT (' nPrint, nt,', T30, 2I6, 2X, 'lines printed, intermed. instants')
 IF (nt > 10000) THEN; WRITE (*, "(' *** Too many instants')"); STOP
 END IF
 WRITE (*, 95700) ishow
95700 FORMAT (' Show coordinates,', T33, I9, T44, '(0|1, no|yes)', T58, 8('..'))
 RETURN
 END SUBROUTINE !END!

 SUBROUTINE timer (k, icount)
! ... ..................................................................
 DIMENSION ival0(8), ival1(8)
 IF (k == 0) THEN
   CALL SYSTEM_CLOCK (icount); CALL DATE_AND_TIME (VALUES=ival0)
   WRITE (*, 95500) ival0
95500 FORMAT (1X, I4, 2('-', I2.2), ' UTC', SP, I5.4, SS, I3, &
      2(':', I2.2), '.', I3.3: T44, 'CPU:', F10.1, ' sec.', T77, 'End')
   RETURN
 ELSE IF (k == 1) THEN
   CALL SYSTEM_CLOCK (icount1, ic_rate); CALL DATE_AND_TIME (VALUES=ival1)
   WRITE (*, *) ' '; WRITE (*, 95500) ival0
   WRITE (*, 95500) ival1, 1.*(icount1-icount)/ic_rate
   RETURN
 END IF
 WRITE (*, "(' ""Timer"": bad argument (not 0 or 1)')")
 RETURN
 END SUBROUTINE !END!
