APPENDIX IV. MODIFIED "CONP.FOR" SOURCE CODE USED IN NUMERICAL MODEL

PROGRAM CONP
C
C     Integration of adiabatic, constant pressure kinetics problems
C     Modified by Javier Molero. July 1995
C     Modified February 1997 for input of temperatures and residence times
C
C*****precision > double
      IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER(I-N)
C*****END precision > double
C
      PARAMETER (LENIWK=5000, LENRWK=13000, LENCWK=500, NK=5, NLMAX=55,
     1           LIN=5, LOUT=6, LINKCK=25, KMAX=61, ITOL=1, IOPT=0,
     2           RTOL=1.0E-6, ITASK=1, ATOL=1.0E-15)
C
      DIMENSION IWORK(LENIWK), RWORK(LENRWK), X(KMAX), Z(KMAX)
      CHARACTER CWORK(LENCWK)*16, KSYM(KMAX)*16, LINE*80
      LOGICAL KERR, IERR
      EXTERNAL FUN
C
      COMMON /RCONS/ P, RU, TT2, T
      COMMON /ICONS/ KK, NWT, NH, NWDOT
C
      DATA KERR/.FALSE./, X/KMAX*0.0/, KSYM/KMAX*' '/
C
C*****open statements included (cif)
      OPEN (LIN, FILE='INP', STATUS='OLD')
      OPEN (LOUT, FILE='CONP.OUT', STATUS='UNKNOWN')
C*****end of open statements !! (cif)
      WRITE (LOUT, 15)
   15 FORMAT (
     1'  CONP:      CHEMKINII  Vers 1.2  Aug. 1992',
C*****precision > double
     2'        DOUBLE PRECISION')
C*****END precision > double
C
C     Open the CHEMKIN LINK file
C
C*****OPEN statement > unix
      OPEN (LINKCK, FORM='UNFORMATTED', file='cklink')
C*****END OPEN statement > unix
C
C*****Initialize CHEMKIN
C
      CALL CKLEN  (LINKCK, LOUT, LENI, LENR, LENC)
      CALL CKINIT (LENIWK, LENRWK, LENCWK, LINKCK, LOUT, IWORK,
     1             RWORK, CWORK)
      CALL CKINDX (IWORK, RWORK, MM, KK, II, NFIT)
C
      NEQ   = KK + 1
      LRW   = 22 + 9*NEQ + 2*NEQ**2
      NVODE = LENR + 1
      NWT   = NVODE + LRW
      NH    = NWT  + KK
      NWDOT = NH   + KK
      NTOT  = NWDOT+ KK - 1
C
      LIW   = 30 + NEQ
      IVODE = LENI + 1
      ITOT  = IVODE + LIW - 1
C
      IF (KK .GT. KMAX) THEN
       WRITE (LOUT, *)
     1   ' Error...KMAX too small...must be at least ',KK
       KERR = .TRUE.
      ENDIF
C
      IF (LENRWK .LT. NTOT) THEN
       KERR = .TRUE.
       WRITE (LOUT, *)
     1   ' Error...LENRWK too small...must be at least', NTOT
      ENDIF
C
      IF (LENIWK .LT. ITOT) THEN
       KERR = .TRUE.
       WRITE (LOUT, *)
     1   ' Error...LENIWK too small...must be at least', ITOT
      ENDIF
C
      IF (KERR) STOP
C
      CALL CKSYMS (CWORK, LOUT, KSYM, IERR)
      IF (IERR) KERR = .TRUE.
      CALL CKWT   (IWORK, RWORK, RWORK(NWT))
      CALL CKRP   (IWORK, RWORK, RU, RUC, PATM)
C
C*****Read initial non-zero moles from input file
C
   40 CONTINUE
      LINE = ' '
      READ  (LIN, '(A)', END=45) LINE
      ILEN = INDEX (LINE, '!')
      IF (ILEN .EQ. 1) GO TO 40
C
      ILEN = ILEN - 1
      IF (ILEN .LE. 0) ILEN = LEN(LINE)
      IF (INDEX(LINE(:ILEN), 'END') .EQ. 0) THEN
       IF (LINE(:ILEN) .NE. ' ') THEN
          CALL CKSNUM (LINE(:ILEN), 1, LOUT, KSYM, KK, KNUM,
     1                   NVAL, VAL, IERR)
          IF (IERR) THEN
             WRITE (LOUT,*) ' Error reading moles...'
             KERR = .TRUE.
          ELSE
             X(KNUM) = VAL
          ENDIF
       ENDIF
       GO TO 40
      ENDIF
C
   45 CONTINUE
C
      IF (KERR) STOP
C
C*****Normalize mole fractions
C
      XTOT = 0.00
      DO 50 K = 1, KK
       XTOT = XTOT + X(K)
   50 CONTINUE
      DO 55 K = 1, KK
       X(K) = X(K) / XTOT
   55 CONTINUE
C
C*****Set time interval for integration constant at 0.001 s
C     Inititalise accumulated residence time
C     Set system pressure constant at 1.15 bar
C
      DT = 0.001
      RTIME = 0.0
      PA = 1.15
      P = PA*PATM
C
C
C*****Read first set of conditions, T and T2
C   
   60 READ (LIN, *) T, T2
C
C*****Convert mole fractions to mass fractions
C
      CALL CKXTY (X, IWORK, RWORK, Z(2))
C
C*****Print spreadsheet heading and first set of conditions 
C      
      WRITE (LOUT, 7020) (K, K=1, KK+2)
      WRITE (LOUT, 7000) (KSYM(K)(:10), K=1, KK)
      WRITE (LOUT, 7010) RTIME, (10**6 * Z(K), K=2, KK+1), T
C
C*****Integration loop begins here
C*****Read new temperature and final time from input file
C
  100 READ (LIN, *) T, T2
C      
C*****End of computation
C      
      IF (T .EQ. 0.0 .AND. T2 .EQ. 0.0) THEN
       STOP
      ENDIF
C
C*****Sets initial conditions and mass fractions and integration
C     parameters for DVODE
C
      TT1  = 0.0
      TT2   = TT1
      Z(1) = T
      MF = 22
      ISTATE= 1
      NLINES=NLMAX + 1
C
C
  250 CONTINUE
C
C*****Calculate accumulated residence time and
C     print solution of present stage to output file
C
      IF (TT2 .GE. T2) THEN
       RTIME = RTIME + T2
       WRITE (LOUT, 7010) RTIME, (10**6 * Z(K), K=2, KK+1), T
       GOTO 100
      ENDIF
C
C*****Increase computation time
C
      TT2 = MIN(TT2 + DT, T2)
C
C*****Call the differential equation solver
C
  350 CONTINUE
C*****precision > double
      CALL DVODE
C*****END precision > double
     *           (FUN, NEQ, Z, TT1, TT2, ITOL, RTOL, ATOL, ITASK,
     1            ISTATE, IOPT, RWORK(NVODE), LRW, IWORK(IVODE),
     2            LIW, JAC, MF, RWORK, IWORK)
C
      IF (ISTATE .LE. -2) THEN
       IF (ISTATE .EQ. -1) THEN
          ISTATE = 2
          GO TO 350
       ELSE
          WRITE (LOUT,*) ' ISTATE=',ISTATE
          STOP
       ENDIF
      ENDIF
      GO TO 250
C
C         FORMATS
C
 7000 FORMAT (4X, 'Time/ss', 1X, 61(1X,A11), 2X, 'Temperat') 
 7010 FORMAT (100E12.4)
 7020 FORMAT (100I12)
      END
C
      SUBROUTINE FUN (N, TIME, Z, ZP, RPAR, IPAR)
C
C*****precision > double
      IMPLICIT DOUBLE PRECISION(A-H,O-Z), INTEGER(I-N)
C*****END precision > double
C
      COMMON /RCONS/ P, RU, TT2, T
      COMMON /ICONS/ KK, NWT, NH, NWDOT
      DIMENSION Z(*), ZP(*), RPAR(*), IPAR(*)
C
C     Variables in Z are:  Z(1) = T, Z(K+1) = Y(K)
C
C     Call CHEMKIN subroutines
C
      CALL CKRHOY (P, T, Z(2), IPAR, RPAR, RHO)
      CALL CKWYP  (P, T, Z(2), IPAR, RPAR, RPAR(NWDOT))
C
C     Form governing equation
C
      DO 100 K = 1, KK
       WDOT    = RPAR(NWDOT + K - 1)
       WT      = RPAR(NWT   + K - 1)
       ZP(K+1) = WDOT * WT / RHO
  100  CONTINUE
C
      RETURN
      END

Previous

Table of Contents