PROGRAM DRIVER C*****precision > double IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) C*****END precision > double C*****precision > single C IMPLICIT REAL (A-H,O-Z), INTEGER (I-N) C*****END precision > single C C LENLWK allocates the logical working space C LENIWK allocates the integer working space C LENRWK allocates the real working space C LENCWK allocates the character working space C LENSYM is the length of a character string C LIN is the unit from which user input is read C LOUT is the unit to which printed output is written C LRIN is the unit from which the restart file is read C LROUT is the unit to which the save file is written C LRCRVR is the unit to which the scratch save file is written C LINKCK is unit from which the Chemkin binary file is read C LINTP is unit from which the Transport binary file is read C PARAMETER (LENLWK=1500, LENIWK=200000, LENRWK=130000000, 1 LENCWK=280, 1 LENSYM=16, LIN=8, LOUT=6, LREST=14, LSAVE=15, 2 LRCRVR=20, LINKCK=25, LINKMC=35,OTB=8.1E9) C C All storage needed by the flame code is allocated in the C following three arrays. LWORK is for logical variables, C IWORK is for integer variables, RWORK is of real variables, C CWORK is for character variables. C DIMENSION I(LENIWK), R(LENRWK) CHARACTER C(LENCWK)*(LENSYM) LOGICAL L(LENLWK) C C NMAX is the total number of grid points allowed NMAX =1000 C open the restart file OPEN (LIN, FILE='shuru.txt') c OPEN (LREST, FORM='UNFORMATTED', STATUS='UNKNOWN', C 1 FILE='rest.bin') C open the save output file C OPEN (LSAVE, FORM='UNFORMATTED', STATUS='UNKNOWN', C 1 FILE='save.bin') OPEN (LSAVE, FILE='SAVE.txt') C open the recover output file OPEN (LRCRVR, FILE='SENS.txt') C open the Chemkin binary file OPEN (LINKCK, FORM='FORMATTED', STATUS='UNKNOWN', 1 FILE='chem.asc') C open the Transport binary file OPEN (LINKMC, FORM='UNFORMATTED', STATUS='UNKNOWN', 1 FILE='tran.bin') OPEN (LREST, FILE='shuchu.txt') C CALL OPPDIF (NMAX, LIN, LOUT, LINKCK, LINKMC, LREST, LSAVE, 1 LRCRVR, LENLWK, L, LENIWK, I, LENRWK, R, 2 LENCWK, C) C STOP END C/////////////////////////////////////////////////////////////////// C C ONE DIMENSIONAL OPPOSED FLOW FLAME CODE, C C WRITTEN BY: C ROBERT J. KEE, ANDREW E. LUTZ C COMPUTATIONAL MECHANICS DIVISION C SANDIA NATIONAL LABORATORIES C LIVERMORE, CA 94550 C TEL (510) 294-3272 C FAX (510) 294-1459 C C///////////////////////////////////////////////////////////////////// C C VERSION 3.3 C 1.1 ADD CORRECTION VELOCITY TO DIFFUSION IN -MDIFV- C COMPUTE SPECIES EQN FOR LAST SPECIES C 1.2 VELOCITY DEPENDENT WINDWARD DIFFERENCING C 1.2A INITIAL G=DFDX=CONST C 1.2B BRIEF PRINT (3 SPECIES REQUESTED) C 1.2C LINEAR INITIAL PROFILE IN MASS FRACTIONS -START- C OPTION FOR COMPUTING STAGNATION PLANE LOCATION C IF XCEN = 0. C 1.2D COMPUTE EQUILIBRIUM AS GUESS FOR ENERGY SOLUTION C 1.2e solve for species with temp fixed at equil guess C 1.2G Multiply reaction rates by factor (RATE). C VERSION 2.0 C Call list for TWOPNT requires additional input; optional C use of new keywords reset default values: C 'ISTP' n - sets NINIT initial time steps before newton C (default is 0) C 'IRET' n - set retirement age IRETIR of old time step C (default 50) C 'NJAC' n - set retirement age NJAC of Jacobian during C steady state newton (default 20) C 'TJAC' n - set retirement age ITJAC of Jacobian during C time stepping (default 20) C 'DTMX' x - set maximum time step DTMAX (default 1.0E-4) C VERSION 2.1 C 1. Initial profile is plateau with equil species. C 2. Use binary files instead of data from restart solution. C VERSION 3.0 C 1. Upgrade to EQUIL.30 C 2. CALL CKRDEX instead of CKRAEX to perturb rates for C sensitivities C CHANGES FOR VERSION 3.1 (3/15/94 F. Rupley) C 1. DOS/PC compatibility effort includes adding file names to C OPEN statements, removing unused variables in CALL lists, C unusued but possibly initialized variables. C CHANGES FOR VERSION 3.2 (1/19/95 F. Rupley) C 1. Add integer error flag to CKLEN, CKINIT, MCLEN, MCINIT C call lists. C CHANGES FOR VERSION 3.3 (2/27/95 F. Rupley) C 1. Change character index "(:" to "(1:" C CHANGES FOR VERSION 3.5 (06/03/99 S. R. Lee) C 1. A optically thin radition model is adpoted. C C///////////////////////////////////////////////////////////////////// C SUBROUTINE OPPDIF (NMAX, LIN, LOUT, LINKCK, LINKMC, LREST, LSAVE, 1 LRCRVR, LENLWK, L, LENIWK, I, LENRWK, R, 2 LENCWK, C) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C DIMENSION I(*), R(*) CHARACTER C(*)*(*) LOGICAL L(*),JIA C COMMON /FLFLFL/ NCKW, NMCW, NYOX, NYFL, NWT, NFL, NPD, NOX, 1 NSCH, NX, NVIS, NCON, NTGV, NXGV, ND, NDKJ, 2 NTDR, NYV, NABV, NBLW, NBUF, NTWP, NS, NSN, 3 NFF, NFN, NDS, NSSV, NA, 4 ICKW, IMCW, ITPW, IKS, ICC, IIP, LAC, LMK, C ADDED FOR EQUIL... 5 IEQW, IATOM, NREAC, NDCS, NEQW C COMMON /LOCS/ NT, NG, NF, NH,NTH ,NYS, NY COMMON TIME_BEGIN, TIME_END COMMON/JI/NUMJIA,NKKK COMMON/HUI/OFIX,OX,OM COMMON/QIAO/ JIA COMMON/FUSHE/IH2O,ICO2 C C C CPU TIME CHECKING C---- CPU Time checking (START) -------------------------- CALL CPU_TIME(TIME_BEGIN) C C WRITE VERSION NUMBER C WRITE (LOUT, 15) 15 FORMAT( 1/' OPPDIF: Opposed-flow diffusion-flame, pressure gradient', 2/' Version 3.6, July 2001. Partly Modified by C.B.Oh', C*****precision > double 3/' DOUBLE PRECISION') C*****END precision > double C*****precision > single C 3/' SINGLE PRECISION') C*****END precision > single C C SET UP INTERNAL WORK POINTERS C c REWIND LSAVE CALL POINT (LINKCK, LINKMC, NMAX, LOUT, LSAVE, KK, II, MM, NATJ, 1 LENIWK, LENRWK, LENCWK, LENITW, LENTWP, 2 LENICK, LENRCK, LENIEQ, LENREQ, 3 LTOT, ITOT, NTOT, ICTOT, I, R, C) C C CHECK FOR ENOUGH SPACE C WRITE (LOUT, 7000) LENLWK, LTOT, LENIWK, ITOT, LENRWK, NTOT, 1 LENCWK, ICTOT 7000 FORMAT (/,' WORKING SPACE REQUIREMENTS', 1 /,' PROVIDED REQUIRED ', 2 /,' LOGICAL ' , 2I15, 3 /,' INTEGER ' , 2I15, 4 /,' REAL ' , 2I15, 5 /,' CHARACTER' , 2I15,/) C IF (LTOT.GT.LENLWK .OR. ITOT.GT.LENIWK .OR. NTOT.GT.LENRWK 1 .OR. ICTOT.GT.LENCWK) THEN WRITE (LOUT, *) ' FATAL ERROR, NOT ENOUGH WORK SPACE PROVIDED' STOP ENDIF C CALL FLDRIV (LENITW, LENTWP, LIN, LOUT, LREST, 1 LSAVE, LRCRVR, KK, II, MM, NATJ, NMAX, LENICK, LENRCK, 2 R(NCKW), R(NMCW), R(NYOX), R(NYFL), R(NWT), R(NFL), 3 R(NPD), R(NOX), R(NSCH), R(NX), R(NVIS), R(NCON), 4 R(NTGV), R(NXGV), R(ND), R(NDKJ), R(NTDR), R(NYV), 5 R(NABV), R(NBLW), R(NBUF), R(NTWP), R(NS), R(NSN), 6 R(NFF), R(NFN), R(NDS), R(NSSV), R(NA), 7 I(ICKW), I(IMCW), I(ITPW), C(IKS), C(ICC), 8 I(IIP), L(LAC), L(LMK), C EQUIL ARRAYS... 9 I(IEQW), LENIEQ, R(NEQW), LENREQ, R(NREAC), 1 R(NDCS), C(IATOM) ) C RETURN END C C---------------------------------------------------------------------- C SUBROUTINE POINT (LINKCK, LINKMC, NMAX, LOUT, LSAVE, KK, II, MM, 1 NATJ, LENIWK, LENRWK, LENCWK, LENITW, LENTWP, 2 LENICK, LENRCK, LENIEQ, LENREQ, 3 LTOT, ITOT, NTOT, ICTOT, I, R, C) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C DIMENSION I(*), R(*) CHARACTER C(*)*(*), ILINK*16 C COMMON /FLFLFL/ NCKW, NMCW, NYOX, NYFL, NWT, NFL, NPD, NOX, 1 NSCH, NX, NVIS, NCON, NTGV, NXGV, ND, NDKJ, 2 NTDR, NYV, NABV, NBLW, NBUF, NTWP, NS, NSN, 3 NFF, NFN, NDS, NSSV, NA, 4 ICKW, IMCW, ITPW, IKS, ICC, IIP, LAC, LMK, C ADDED FOR EQUIL... 5 IEQW, IATOM, NREAC, NDCS, NEQW C COMMON /LOCS/ NT, NG, NF, NH,NTH ,NYS, NY C C C SET THE POINTERS INTO THE SOLUTION VECTOR C NT = 1 NG = 2 NF = 3 NH = 4 NTH = 5 NYS = 5 NY = 6 C C CALL FOR CHEMKIN WORK ARRAY LENGTHS C c CALL CKLEN (LINKCK, LOUT, LENICK, LENRCK, LENCCK) c CALL MCLEN (LINKMC, LOUT, LENIMC, LENRMC) CALL CKLEN (LINKCK, LOUT, LENICK, LENRCK, LENCCK, IFLAG1) CALL MCLEN (LINKMC, LOUT, LENIMC, LENRMC) C IF (IFLAG1.GT.0 .OR. IFLAG2.GT.0) STOP C C real chemkin work space NCKW = 1 C real transport work space NMCW = NCKW + LENRCK NTOT = NMCW + LENRMC C C integer chemkin space ICKW = 1 C integer transport space IMCW = ICKW + LENICK ITOT = IMCW + LENIMC C C characcter chemkin space ICC = 1 ICTOT = ICC + LENCCK C IF (ITOT.LT.LENIWK .AND. NTOT.LT.LENRWK .AND. ICTOT.LT.LENCWK) 1 THEN C CALL CKINIT (LENICK, LENRCK, LENCCK, LINKCK, LOUT, I, R, C) CALL CKINIT (LENICK, LENRCK, LENCCK, LINKCK, LOUT, I, R, C, 1 IFLAG) C IF (IFLAG .GT. 0) STOP ILINK = 'CKLINK' C WRITE (LSAVE) ILINK C CALL CKSAVE (LOUT, LSAVE, I, R, C) CALL CKINDX (I, R, MM, KK, II, NFIT) C CALL MCINIT (LINKMC, LOUT, LENIMC, LENRMC, I(IMCW), R(NMCW)) C CALL MCINIT (LINKMC, LOUT, LENIMC, LENRMC, I(IMCW), R(NMCW), C 1 IFLAG) C IF (IFLAG .GT. 0) STOP ILINK = 'MCLINK' C WRITE (LSAVE) ILINK C CALL MCSAVE (LOUT, LSAVE, I(IMCW), R(NMCW)) ENDIF C C TWOPOINT WORK SPACE LENGTHS C NATJ = KK+5 LENITW = 3 * NMAX LENTWP = (7*NATJ+2) * NMAX C C STANJAN WORK SPACE LENGTHS C NCON = 0 CALL EQLEN (MM, KK, NCON, LENIEQ, LENREQ) C C APPORTION THE FLOATING POINT SPACE C NYOX = NTOT NYFL = NYOX + KK NWT = NYFL + KK NFL = NWT + KK NPD = NFL + KK NOX = NPD + KK NSCH = NOX + KK NX = NSCH + KK*5 NVIS = NX + NMAX NCON = NVIS + NMAX NTGV = NCON + NMAX NXGV = NTGV + NMAX ND = NXGV + NMAX NDKJ = ND + KK*NMAX NTDR = NDKJ + KK*KK*NMAX NYV = NTDR + KK*NMAX NABV = NYV + KK*NMAX NBLW = NABV + NATJ*NMAX NBUF = NBLW + NATJ*NMAX NTWP = NBUF + NATJ*NMAX NS = NTWP + LENTWP NSN = NS + NATJ*NMAX NFF = NSN + NATJ*NMAX NFN = NFF + NATJ*NMAX NDS = NFN + NATJ*NMAX NSSV = NDS + NMAX NA = NSSV + NMAX C NTOT = NA + (6*NATJ-2) * (NATJ*NMAX) - 1 NEQW = NA + (6*NATJ-2) * (NATJ*NMAX) - 1 NREAC = NEQW + LENREQ NDCS = NREAC + 2*KK NTOT = NDCS + KK C C APPORTION THE INTEGER SPACE C ITPW = ITOT IIP = ITPW + LENITW C ITOT = IIP + NATJ*NMAX - 1 IEQW = IIP + NATJ*NMAX - 1 ITOT = IEQW + LENIEQ C C APPORTION THE LOGICAL SPACE C LAC = 1 LMK = LAC + NATJ LTOT = LMK + NMAX - 1 C C APPORTION THE CHARACTER SPACE C IKS = ICTOT IATOM = IKS + KK ICTOT = IATOM + MM RETURN END C C---------------------------------------------------------------------- C SUBROUTINE FLDRIV (LENITW, LENTWP, LIN, LOUT, 1 LREST, LSAVE, LRCRVR, KK, II, MM, NATJ, NMAX, 2 LENICK, LENRCK, RCKWRK, RMCWRK, YOXID, YFUEL, 3 WT, FUEL, PROD, OXID, SCRCHK, X, VISC, COND, 4 TGIVEN, XGIVEN, D, DKJ, TDR, YV, ABOVE, BELOW, 5 BUFFER, TWPWK, S, SN, F, FN, DS, SSAVE, A, 6 ICKWRK, IMCWRK, ITWPWK, KSYM, CCKWRK, IP, 7 ACTIVE, MARK, C EQUIL ARRAYS... 7 IEQWRK, LENIEQ, REQWRK, LENREQ, REAC, 8 DCS, ATOM ) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C COMMON /LOCS/ NT, NG, NF, NH,NTH ,NYS, NY COMMON/JI/NUMJIA,NKKK COMMON/HUI/OFIX,OX,OM C DIMENSION ICKWRK(*), RCKWRK(*), IMCWRK(*), RMCWRK(*), YOXID(*), 1 YFUEL(*), WT(*), FUEL(*), PROD(*), OXID(*), X(*), 2 VISC(*), COND(*), D(KK,*), TDR(KK,*), YV(KK,*), 3 DKJ(KK,KK,*), S(NATJ,*), SN(NATJ,*), F(NATJ,*), 4 XGIVEN(*), TGIVEN(*), SCRCHK(KK,5), KW(6),KWW(1), 5 WDOT(KK),Q(II),QM(KK),RK(II),C(KK), HML(KK) C C EQUIL SPACE C DIMENSION IEQWRK(LENIEQ), REQWRK(LENREQ), REAC(KK,2), DCS(KK) C C DIMENSION NEWTON SPACE C DIMENSION ABOVE(NATJ,*), BELOW(NATJ,*), BUFFER(NATJ,*), 1 TWPWK(*), ITWPWK(*), A(6 * NATJ - 2, *), IP(NATJ, *), 3 FN(NATJ,*), DS(*), SSAVE(*), LEVEL(2) C CHARACTER KSYM(*)*(*), CCKWRK(*)*(*), ATOM(*)*(*) C LOGICAL LTIME, LTIME2, LENRGY, LMULTI, LTDIF, LUMESH, LRSTRT, 1 LCNTUE, LSEN, LVARMC, RSTCNT, ERROR, FUNCTN, JACOBN, 2 REENTR, SOLVE, STORE, SUCCES, ADAPT, SHOW, SAVE, UPDATE, 3 ACTIVE(*), MARK(*), ENERGY, LUSTGV, LREGRD,JIA,LFXTEQ COMMON/QIAO/ JIA COMMON/FUSHE/IH2O,ICO2 C INTEGER CALL, CALLS C CHARACTER*16 ISOLUT, ISENSI, ICHR, ICKLNK, IMCLNK, WUZHI LOGICAL IERR, LEQUIL DATA ISOLUT/'SOLUTION '/, ISENSI/'SENSITIVITY '/, 1 ICKLNK/'CKLINK '/, IMCLNK/'MCLINK '/ DATA ADAFLR /1.0E-08/, LCNTUE/.FALSE./, RFAC/1.0/ C C COMPUTE THE UNIT ROUNDOFF, AND THE RELATIVE AND ABSOLUTE C PERTURBATIONS FOR THE JACOBIAN EVALUATION. C OO=1 NNUMBER=1 OSS=0.65 OWUZHI=0.0 JIA=.TRUE. WUZHI='OH' U = 1.0 10 CONTINUE U = U*0.5 COMP = 1.0 + U IF (COMP .NE. 1.0) GO TO 10 ABSOL = SQRT(2.0*U) RELAT = SQRT(2.0*U) C RSTCNT = .FALSE. C REWIND LRCRVR C CALL CKSYMS (CCKWRK, LOUT, KSYM, IERR) CALL CKWT (ICKWRK, RCKWRK, WT) CALL CKRP (ICKWRK, RCKWRK, RU, RUC, PATM) CALL CKCRAY (WUZHI, KK, KSYM, LOUT, 1, KWW, NFD, IERR) NUMJIA=KWW(1) CALL CKCRAY ('H2O', KK, KSYM, LOUT, 1, KWW, NFD, IERR) IH2O=KWW(1) CALL CKCRAY ('CO2', KK, KSYM, LOUT, 1, KWW, NFD, IERR) ICO2=KWW(1) CALL CKCRAY ('C2H4', KK, KSYM, LOUT, 1, KWW, NFD, IERR) NKKK=KWW(1) WRITE (LOUT, *)NKKK C 200 CONTINUE C C CALL THE KEYWORD INPUT C CALL RDKEY (KK, NMAX, LIN, LOUT, KSYM, PATM, LTIME, LTIME2, 1 LUSTGV, LENRGY, LEQUIL, LMULTI, LTDIF, LUMESH, LRSTRT, 2 LCNTUE, IPROFL, IPRNT, MFILE, UFAC,DFAC,RATEF,N1CALL, 3 LSEN, VFUEL, VOXID, AFUEL, AOXID, TFUEL, TOXID, TMAX, 4 P, JJ, X,FUEL, PROD, OXID, XCEN, XEND, WMIX, 5 GRAD, CURV, NADP, IRETIR, SFLR, NTEMP, XGIVEN, TGIVEN, 6 ATOL, RTOL, ATIM, RTIM, NUMDT, DT1, NUMDT2, DT2, 7 DTMIN, WNDFAC, LREGRD, JJREGD, PCTADP, RATGTC, 8 KW, NW, LFXTEQ, NJAC, NINIT, ITJAC, DTMAX ) C C C SET THE SOLUTION BOUNDS C CALL CKXTY (OXID, ICKWRK, RCKWRK, YOXID) CALL CKXTY (FUEL, ICKWRK, RCKWRK, YFUEL) CALL CKRHOY (P, TFUEL, YFUEL, ICKWRK, RCKWRK, RHOF) CALL CKRHOY (P, TOXID, YOXID, ICKWRK, RCKWRK, RHOO) VOXID=-(VFUEL*sqrt(RHOF/RHOO)) 201 CONTINUE DO 100 J = 1, NMAX BELOW(NT,J) = 200. ABOVE(NT,J) = 6000.0E0 BELOW(NG,J) = -10000. ABOVE(NG,J) = 10000. BELOW(NH,J) = -100.E10 ABOVE(NH,J) = 100.E10 BELOW(NF,J) = -10000. ABOVE(NF,J) = 10000. BELOW(NTH,J) = 200. ABOVE(NTH,J) = 1800. DO 100 K = 1, KK BELOW (NYS+K, J) = SFLR ABOVE (NYS+K, J) = 10. 100 CONTINUE C ENERGY = LENRGY C C C SET THE VARIABLES ON WHICH TO ADAPT C DO 300 N = 1, NATJ ACTIVE(N) = .TRUE. 300 CONTINUE C IF (IPRNT .LT. 10) THEN LEVEL(1) = MIN (2, IPRNT) LEVEL(2) = LEVEL(1) ELSE LEVEL(1) = IPRNT/10 LEVEL(2) = IPRNT - 10*LEVEL(1) LEVEL(1) = MAX( LEVEL(1), LEVEL(2) ) ENDIF C IF (LRSTRT) THEN IF (.NOT. RSTCNT) THEN C C THIS IS A RESTART C NREST = 0 C320 CONTINUE C ENDIF C C LIMIT SPECIES MINIMUM ON RESTART C CC DO 326 J = 2, JJ CC DO 325 K = 1, KK CC CC*****precision > single CC S(NYS+K,J) = MAX( 1.E-15, S(NYS+K,J) ) CC*****END precision > single CC*****precision > double CC S(NYS+K,J) = MAX( 1.D-15, S(NYS+K,J) ) CC*****END precision > double CC CC325 CONTINUE CC326 CONTINUE C CALL RESTRT (KK, NMAX, NATJ, JJ, LOUT, LUSTGV, FUEL, 1 OXID, NTEMP, XGIVEN, TGIVEN, XEND, 2 ICKWRK, RCKWRK, YOXID, YFUEL, X, S) C C IF (RSTCNT .AND. LREGRD) THEN IF (LREGRD) THEN WRITE (LOUT, *) ' REGRIDDING TO ', JJREGD, ' POINTS' CALL REGRID (SN(1,2), S(1,2), NATJ, JJREGD-1, JJ-1, 1 VISC(2), X(2), COND, PCTADP, RATGTC, NT) JJ = JJREGD DO 380 J = 2, JJ X(J) = VISC(J) DO 370 N = 1, NTH S(N,J) = SN(N,J) 370 CONTINUE DO 375 K = 1, KK C C*****precision > single C S(NYS+K,J) = MAX( 1.E-15, SN(NYS+K,J) ) C*****END precision > single C*****precision > double S(NYS+K,J) = MAX( 1.D-15, SN(NYS+K,J) ) C*****END precision > double C 375 CONTINUE 380 CONTINUE C NTEMP = JJ DO 390 N = 1, NTEMP XGIVEN(N) = X(N) TGIVEN(N) = S(NT,N) 390 CONTINUE ENDIF C ELSE C C SET THE STARTING PROFILES C CALL START (KK, NMAX, NATJ, JJ, LOUT, LUMESH, IPROFL, TMAX, 1 TFUEL, TOXID, VFUEL, VOXID, P, 2 FUEL, PROD, OXID, NTEMP, XGIVEN, TGIVEN, XCEN, 3 XEND, WMIX, ICKWRK, RCKWRK, SCRCHK(1,4), YOXID, 4 YFUEL, X, S) C IF (LEQUIL .AND. (IPROFL.EQ.3) ) THEN PA = P / PATM CALL EQGUES ( JJ, KK, MM, LOUT, NATJ, S, PA, ATOM, KSYM, 2 ICKWRK, LENICK, RCKWRK, LENRCK, IEQWRK, LENIEQ, REQWRK, 3 LENREQ, REAC, DCS, LFXTEQ ) NTEMP = JJ DO 420 N = 1, NTEMP XGIVEN(N) = X(N) TGIVEN(N) = S(NT,N) 420 CONTINUE LEQUIL = .FALSE. ENDIF C ENDIF C C C HOW MANY CALLS TO TWOPNT C IF (ENERGY .AND. LEQUIL ) THEN CALLS = 3 ELSEIF (ENERGY .AND. .NOT. LEQUIL ) THEN CALLS = 2 ELSE CALLS = 1 ENDIF C C/// TOP OF THE LOOP OVER CALLS TO TWOPNT. C C DO 1100 CALL = N1CALL, CALLS C IF (CALL .EQ. 1) THEN LENRGY = .FALSE. ADAPT = .FALSE. ELSEIF (CALL .EQ. 2) THEN IF (LEQUIL) THEN PA = P / PATM CALL EQGUES 1 ( JJ, KK, MM, LOUT, NATJ, S, PA, ATOM, KSYM, 2 ICKWRK, LENICK, RCKWRK, LENRCK, IEQWRK, LENIEQ, REQWRK, 3 LENREQ, REAC, DCS, LFXTEQ ) NTEMP = JJ DO 450 N = 1, NTEMP XGIVEN(N) = X(N) TGIVEN(N) = S(NT,N) 450 CONTINUE ELSE LENRGY = ENERGY ADAPT = .TRUE. WRITE (LOUT,'(/A/)') 1 ' FLDRIV: FINISHED FIXED TEMPERATURE, ADDING ENERGY EQUATION' ENDIF ELSEIF (CALL .EQ. 3) THEN LENRGY = ENERGY ADAPT = .TRUE. WRITE (LOUT,'(/A/)') 1 ' FLDRIV: FINISHED FIXED TEMPERATURE, ADDING ENERGY EQUATION' ENDIF C IF (LENRGY) RFAC = RATEF C IF ((LTIME2) .AND. CALL.GT.1) THEN NUMDT = NUMDT2 DT1 = DT2 ENDIF C RSTCNT = .FALSE. C IPASSS = 1 500 CONTINUE C CALL TWOPNT (ERROR, LOUT, LEVEL, LENITW, ITWPWK, LENTWP, TWPWK, 1 ABOVE, ACTIVE, ADAPT, BELOW, 0, BUFFER, NATJ, 2 CONDIT, FUNCTN, JACOBN, MARK, X, IPASS, 3 IPASSS, NADP, NMAX, JJ, REENTR, IREPRT, 4 SAVE, 0, SHOW, SOLVE, ATOL, NJAC, RTOL, 5 NINIT, NUMDT, IRETIR, STORE, SUCCES, 6 ATIM, ITJAC, DFAC, RTIM, LTIME, UFAC, DTMAX, DTMIN, 7 ADAFLR, GRAD, CURV, DT1, DT, UPDATE, S) C IF (ERROR) THEN STOP ELSE IF (.NOT.SUCCES .AND. .NOT.REENTR) THEN IF (IREPRT .EQ. 2) WRITE (LOUT, *) 1 ' TWOPNT requires more mesh points, but NMAX too small' C C ELSE IF (REENTR) THEN C IF (FUNCTN) THEN C LVARMC = .TRUE. CALL FUN (KK, JJ, NATJ, LENRGY, LMULTI, LTDIF, LVARMC, 1 LTIME, P, WT, YOXID, YFUEL, DT, NTEMP, XGIVEN, 2 TGIVEN, X, SN, BUFFER, WNDFAC, SCRCHK(1,1), 3 SCRCHK(1,5), YV, SCRCHK(1,2), SCRCHK(1,3), 4 SCRCHK(1,4), TFUEL, TOXID, VFUEL, VOXID, 4 AFUEL, AOXID, 5 VISC, COND, D, DKJ, TDR, 6 ICKWRK, RCKWRK, IMCWRK, RMCWRK, F, RFAC) C C*****precision > double CALL DCOPY (NATJ*JJ, F, 1, BUFFER, 1) C*****END precision > double C C*****precision > single C CALL SCOPY (NATJ*JJ, F, 1, BUFFER, 1) C*****END precision > single C C ELSE IF (JACOBN) THEN C CALL JACOB (KK, JJ, NATJ, LENRGY, LMULTI, LTDIF, 1 LTIME, P, WT, YOXID, YFUEL, DT, NTEMP, XGIVEN, 2 TGIVEN, X, SN, BUFFER, WNDFAC, ABSOL, RELAT, 3 SCRCHK, YV, TFUEL, TOXID, VFUEL, VOXID, 3 AFUEL, AOXID, 4 VISC, COND, D, DKJ, TDR, ICKWRK, RCKWRK, 5 IMCWRK, RMCWRK, F, FN, A, DS, SSAVE, RFAC) C C C*****precision > double CALL DGBCO C*****END precision > double C*****precision > single C CALL SGBCO C*****END precision > single + (A, 6 * NATJ - 2, NATJ * JJ, 2 * NATJ - 1, 2 * NATJ - 1, + IP, RCOND, FN) C IF (RCOND .LE. 0.0E0) THEN WRITE (LOUT, *) ' FATAL ERROR, SINGULAR JACOBIAN ' STOP ENDIF CONDIT = 1.0 / RCOND C ELSE IF (SOLVE) THEN C C*****precision > double CALL DGBSL C*****END precision > double C*****precision > single C CALL SGBSL C*****END precision > single + (A, 6 * NATJ - 2, NATJ * JJ, 2 * NATJ - 1, 2 * NATJ - 1, + IP, BUFFER, 0) C ELSE IF (STORE) THEN C C*****precision > double CALL DCOPY (NATJ*JJ, BUFFER, 1, SN, 1) C*****END precision > double C DO 1040 J = 1, JJ DO 1035 K = 1, KK SN(NYS+K,J) = MAX (SN(NYS+K,J), 1.E-20) 1035 CONTINUE 1040 CONTINUE C*****precision > single C CALL SCOPY (NATJ*JJ, BUFFER, 1, SN, 1) C*****END precision > single C ELSE IF (SHOW) THEN CALL PRINT (LOUT, KK, JJ, NATJ, P, X, BUFFER, SCRCHK(1,1), 1 KSYM, ICKWRK, RCKWRK, KW, NW,LREST) C ELSE IF (SAVE) THEN C REWIND LRCRVR C WRITE (LRCRVR) ISOLUT C WRITE (LRCRVR) NATJ, JJ, P C WRITE (LRCRVR) VFUEL, VOXID C WRITE (LRCRVR) (X(J), J=1,JJ) C WRITE (LRCRVR) ((BUFFER(N,J), N=1,NATJ), J=1,JJ) C ELSE IF (UPDATE) THEN IF (.NOT. LENRGY) THEN DO 900 J = 1, JJ CALL TEMP (NTEMP, X(J), XGIVEN, TGIVEN, TI) BUFFER(NT,J) = TI 900 CONTINUE ENDIF IF(JIA)THEN DO 999 J=1,JJ-1 BUFFER(NTH,J)=BUFFER(NTH,JJ) 999 CONTINUE ENDIF C ENDIF C GO TO 500 C ENDIF C 1100 CONTINUE C C WRITE TO LSAVE WHEN SOLUTION IS COMPLETE C IF (LEVEL(2) .EQ. 0) 1 CALL PRINT (LOUT, KK, JJ, NATJ, P, X, S, SCRCHK(1,1), 2 KSYM, ICKWRK, RCKWRK, KW, NW,LREST) C C WRITE (LSAVE) ISOLUT C WRITE (LSAVE) NATJ, JJ, P C WRITE (LSAVE) VFUEL, VOXID C WRITE (LSAVE) (X(J), J=1,JJ) C WRITE (LSAVE) ((S(N,J), N=1,NATJ), J=1,JJ) C C c IF(ABS(S(NTH,1)- 1237.531).LT.1)OO=2 c LSEN= .TRUE. IF(OO.EQ.2) THEN c LDA = 6*NATJ - 2 c CALL BSEN (II, KK, JJ, NATJ, LDA, LENRGY, LMULTI, LTDIF, c 1 LRCRVR, LOUT, LVARMC, LTIME, P, WT, c 2 YOXID, YFUEL, DT, NTEMP, XGIVEN, c 3 TGIVEN, X, SN, S, WNDFAC, ABSOL, RELAT, c 4 SCRCHK, YV, VISC, COND, D, DKJ, TDR, ICKWRK, c 5 RCKWRK, IMCWRK, RMCWRK, F, FN, A, DS, SSAVE, c 6 IP, RFAC) WRITE (LOUT,'(/A/)') ' SENSITIVITY CALCULATION COMPLETE' LSEN= .TRUE. DO 7821 J=2,JJ-1 CALL CKWYP(P,S(NT,J),S(NY,J),ICKWRK,RCKWRK,WDOT) c ת³ÉĦ¶û·ÖÊý CALL CKYTX(S(NY,J),ICKWRK,RCKWRK,QM) c ÎïÖÊŨ¶È mole/cm**3 CALL CKXTCP(P,S(NT,J),QM,ICKWRK,RCKWRK,C) C ¸ø¶¨Î¶ȺÍÎïÖÊŨ¶Èϵķ´Ó¦ËÙÂÊ CALL CKQC(S(NT,J),C,ICKWRK,RCKWRK,Q) c ¶ÔijÎïÖʵÄÓй±Ï׵ķ´Ó¦ CALL CKCONT(NKKK,Q,ICKWRK,RCKWRK,RK) CALL CKHML (S(NT,J), ICKWRK, RCKWRK, HML) IF (S(NF,J) .LE. 0.) THEN WINDA = 1. WINDB = -1. WINDC = 0. ELSE WINDA = 0. WINDB = 1. WINDC = -1. ENDIF DXWIND = WINDA*X(J+1) + WINDB*X(J) + WINDC*X(J-1) C C DYK = WINDA*S(NYS+NUMJIA,J+1) + WINDB*S(NYS+NUMJIA,J) 1 + WINDC*S(NYS+NUMJIA,J-1) C RCON= 2.0*S(NF,J) * DYK / DXWIND C RAC= WDOT(NUMJIA) * WT(NUMJIA) C RDIF=-RCON + RAC c WRITE(LSAVE,816)X(J),RCON,RAC,RDIF c WRITE(LSAVE,799)X(J),Q(38) DO 7822 I=1,II CALL CKHRX(I,HML,ICKWRK,RCKWRK,HRXI) HREC=HRXI*Q(I) WRITE(LSAVE,799)X(J),Q(I),HREC,RK(I) 7822 CONTINUE 7821 CONTINUE ENDIF 799 FORMAT(2X,F6.3,3X,E16.3,3X,E16.3,3X,E16.3) c816 FORMAT(2X,'X=',3X,F6.3,2X,E10.3,2X,E10.3,2X,E10.3) C----------------------------------------------------------------------------- IF (LSEN) THEN C LDA = 6*NATJ - 2 CALL REASEN (II, KK, JJ, NATJ, LDA, LENRGY, LMULTI, LTDIF, 1 LRCRVR, LOUT, LVARMC, LTIME, P, WT, 2 YOXID, YFUEL, DT, NTEMP, XGIVEN, 3 TGIVEN, X, SN, S, WNDFAC, ABSOL, RELAT, 4 SCRCHK, YV, VISC, COND, D, DKJ, TDR, ICKWRK, 5 RCKWRK, IMCWRK, RMCWRK, F, FN, A, DS, SSAVE, 6 IP, RFAC) C WRITE (LOUT,'(/A/)') ' SENSITIVITY CALCULATION COMPLETE' C STOP ENDIF C C CHECK FOR CONTINUATION C IF (LCNTUE) THEN C WRITE (LOUT,'(/////)') DO 1210 L = 1, 5 WRITE (LOUT, *) 1 ' ////////////////// CONTINUING TO NEW PROBLEM /////////////' 1210 CONTINUE WRITE (LOUT,'(/////)') C RSTCNT = .TRUE. LRSTRT = .TRUE. C GO TO 200 ENDIF IF(OSS.LT.0.66)THEN DO 77 J =1,JJ-1 IF (S(NYS+NUMJIA,J+1).GT.OWUZHI) THEN OWUZHI=S(NYS+NUMJIA,J+1) ENDIF 77 CONTINUE ENDIF IF(OSS.LT.0.66)THEN DO 78 J = 1,JJ IF(S(NYS+NUMJIA,J).GE.(OSS*OWUZHI))THEN OX=X(J) GO TO 79 ENDIF 78 CONTINUE ENDIF 79 CONTINUE OFIX=OSS*OWUZHI OSS=OSS+0.5 c IF(NNUMBER.GT.20)OSS=OSS+2000 JIA=.FALSE. NNUMBER=NNUMBER+1 WRITE(LOUT,*)OFIX WRITE(LOUT,*)S(NTH,1) c LTDIF=.TRUE. IF(NNUMBER.LT.5000)THEN ENERGY=.TRUE. LEQUIL=.FALSE. RSTCNT =.TRUE. LRSTRT =.TRUE. N1CALL=2 GO TO 201 ENDIF STOP END C C--------------------------------------------------------------------- C SUBROUTINE EQGUES 1 ( JJ, KK, MM, LOUT, NATJ, SOLN, PRES, ATOM, KSYM, 2 ICKWRK, LENICK, RCKWRK, LENRCK, IEQWRK, LENIEQ, REQWRK, 3 LENREQ, REAC, DCS, LFXTEQ ) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C PARAMETER (NMAX=1) DIMENSION ICKWRK(LENICK), RCKWRK(LENRCK), IEQWRK(LENIEQ), 1 REQWRK(LENREQ), REAC(KK,2), DCS(*), 4 XCON(NMAX), KCON(NMAX), SOLN(NATJ,*) C CHARACTER ATOM(*)*(*), KSYM(*)*(*) C LOGICAL EQST, LFXTEQ C COMMON /LOCS/ NT, NG, NF, NH,NTH ,NYS, NY C DATA EQST/.FALSE./, KMON/0/, NCON/0/ C IF (LFXTEQ) THEN NOP = 1 ELSE NOP = 5 ENDIF C C LOOP OVER INNER GRID POINTS, COMPUTE LOCAL EQUIL C DO 100 J = 2, JJ-1 C C PASS IN SPECIES MOLE FRACTIONS C CALL CKYTX (SOLN(NY,J), ICKWRK, RCKWRK, REAC(1,1)) C TE = SOLN(NT,J) PE = PRES CALL EQUIL (LOUT, .FALSE., 0, EQST, .FALSE., ICKWRK, RCKWRK, 1 LENIEQ, IEQWRK, LENREQ, REQWRK, MM, KK, ATOM, 2 KSYM, NOP, KMON, REAC(1,1), SOLN(NT,J), TE, 3 PRES, PE, NCON, KCON, XCON) C CALL EQSOL (KK, REQWRK, REAC(1,2), DCS, SOLN(NT,J), PDUM, 1 HDUM, VDUM, SDUM, WMDUM, CDUM, CDET) C CALL CKXTY (REAC(1,2), ICKWRK, RCKWRK, SOLN(NY,J)) 100 CONTINUE C RETURN END C C--------------------------------------------------------------------- C SUBROUTINE REGRID (SNEW, SOLD, NVAR, JJNEW, JJOLD, XNEW, XOLD, 1 WORK, P, R, N ) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C DIMENSION SNEW(NVAR,*), SOLD(NVAR,*), XNEW(*), XOLD(*), WORK(*) C C interpolate the solution onto an equidistributing mesh C C C compute weight coefficients C R0 = 1. - P R1 = P *R /(R+1.) R2 = P - R1 TV1 = 0. DO 10 I = 2, JJOLD TV1 = TV1 + ABS( SOLD(N,I) - SOLD(N,I-1) ) 10 CONTINUE TV2 = 0. DO 20 I = 2, JJOLD-1 TV2 = TV2 + ABS( (SOLD(N,I+1)-SOLD(N,I))/(XOLD(I+1)-XOLD(I)) 1 - (SOLD(N,I)-SOLD(N,I-1))/(XOLD(I)-XOLD(I-1)) ) 20 CONTINUE XLEN = XOLD(JJOLD) - XOLD(1) B1 = R1 * XLEN /(R0 * TV1) B2 = R2 * XLEN /(R0 * TV2) C C compute partial sums of weight function C WORK(1) = 0. DO 50 I = 2, JJOLD DX = XOLD(I) -XOLD(I-1) WORK(I) = DX + B1*ABS( SOLD(N,I) - SOLD(N,I-1) ) + WORK(I-1) 1 + B2*ABS( (SOLD(N,I+1)-SOLD(N,I))/(XOLD(I+1)-XOLD(I)) 2 - (SOLD(N,I)-SOLD(N,I-1))/(XOLD(I)-XOLD(I-1)) ) 50 CONTINUE DO 65 I = 2, JJOLD WORK(I) = WORK(I)/WORK(JJOLD) 65 CONTINUE C C interpolate onto uniform eta grid to find new x C XNEW(1) = XOLD(1) XNEW(JJNEW) = XOLD(JJOLD) ISTART = 2 DETA = 1./FLOAT(JJNEW-1) DO 80 J = 2, JJNEW-1 ETAJ = (J-1)*DETA DO 70 I = ISTART, JJOLD IF (ETAJ .LE. WORK(I)) THEN DEL = (ETAJ-WORK(I-1))/(WORK(I)-WORK(I-1)) XNEW(J) = XOLD(I-1)+(XOLD(I)-XOLD(I-1))*DEL GO TO 80 ELSE ISTART = I ENDIF 70 CONTINUE WRITE (6, *) ' *** VALUE OF ETA NOT FOUND ***' 80 CONTINUE C C interpolate solution... C CALL INTPL8 (SOLD, SNEW, XOLD, XNEW, NVAR, JJOLD, JJNEW) C RETURN END C C------------------------------------------------------------------ SUBROUTINE INTPL8 (F1, F2, X1, X2, MVAR, N1, N2) C------------------------------------------------------------------ C C interpolate to get f2(x2) from f1(x1) C C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C DIMENSION F1(MVAR,*), F2(MVAR,*), X1(*), X2(*) DO 100 J = 2, N2-1 XVAL = X2(J) DO 50 I = 2, N1 IF(XVAL.LE.X1(I)) THEN DEL = (XVAL-X1(I-1))/(X1(I)-X1(I-1)) DO 40 L = 1, MVAR F2(L,J)=F1(L,I-1)+(F1(L,I)-F1(L,I-1))*DEL 40 CONTINUE GOTO 80 ENDIF 50 CONTINUE WRITE(7,*) ' *** STOP...INTERPOLATION ERROR.' 80 CONTINUE 100 CONTINUE C C endpoints.. C DO 110 L = 1, MVAR F2(L,1) = F1(L,1) F2(L,N2) = F1(L,N1) 110 CONTINUE END C C------------------------------------------------------------------- C SUBROUTINE JACOB (KK, POINTS, COMPS, LENRGY, LMULTI, LTDIF, 1 LTIME, P, WT, YOXID, YFUEL, DT, 2 NTEMP, XGIVEN, TGIVEN, MESH, SN, X0, WNDFAC, 3 ABSOL, RELAT, SCRTCH, YV, 4 TFUEL, TOXID, VFUEL, VOXID, AFUEL, AOXID, 4 VISC, COND, D, DKJ, 5 TDR, ICKWRK, RCKWRK, IMCWRK, RMCWRK, Y1, Y0, 6 A, PERTRB, SAVE, RFAC) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) DOUBLE PRECISION MESH C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C REAL MESH C*****END precision > single C INTEGER CLASS, COMP, COMPS, INDEX, LOWER, POINT, POINTS, ROWS, 1 UPPER LOGICAL LENRGY, LMULTI, LTDIF, LTIME C PARAMETER (LOWER = 1, UPPER = 1, MODULS = LOWER + 1 + UPPER) C DIMENSION A((2 * LOWER + UPPER + 3) * COMPS - 2, *), 1 MESH(*), PERTRB(*), SAVE(*), SCRTCH(KK,*), 2 X0(COMPS, *), Y0(COMPS, *), Y1(COMPS, *),WT(*), 3 YOXID(*),YFUEL(*),XGIVEN(*),TGIVEN(*),SN(COMPS,*), 4 YV(KK,*),VISC(*),COND(*),D(KK,*), DKJ(KK,KK,*), 5 TDR(KK,*),ICKWRK(*), RCKWRK(*),IMCWRK(*), RMCWRK(*) C KOFSET = (LOWER + UPPER + 2) * COMPS - 1 ROWS = (2 * LOWER + UPPER + 3) * COMPS - 2 C C/// ZERO THE MATRIX STORAGE SPACE. C ZERO = 0.0 C*****precision > double CALL DCOPY (ROWS * COMPS * POINTS, ZERO, 0, A, 1) C*****END precision > double C C*****precision > single C CALL SCOPY (ROWS * COMPS * POINTS, ZERO, 0, A, 1) C*****END precision > single C C C/// CALL THE FUNCTION AT X0 AND STORE IN Y0. C CALL FUN (KK, POINTS, COMPS, LENRGY, LMULTI, LTDIF, .TRUE., 1 LTIME, P, WT, YOXID, YFUEL, DT, NTEMP, XGIVEN, 2 TGIVEN, MESH, SN, X0, WNDFAC, SCRTCH(1, 1), 3 SCRTCH(1,5), YV, SCRTCH(1, 2), SCRTCH(1, 3), 4 SCRTCH(1, 4), TFUEL, TOXID, VFUEL, VOXID, 4 AFUEL, AOXID, 5 VISC, COND, D, DKJ, TDR, ICKWRK, RCKWRK, IMCWRK, 6 RMCWRK, Y0, RFAC) C C C/// TOP OF THE LOOPS OVER THE RESIDUE CLASSES AND SOLUTION COMPONENTS. C DO 0400 CLASS = 1, MODULS DO 0400 COMP = 1, COMPS C C/// FOR A GIVEN RESIDUE CLASS AND A GIVEN SOLUTION COMPONENT, C/// PERTRB THE X0 VECTOR AT POINTS IN THE SAME RESIDUE CLASS. C INDEX = 0 DO 0100 POINT = CLASS, POINTS, MODULS INDEX = INDEX + 1 SAVE(INDEX) = X0(COMP, POINT) PERTRB(INDEX) = ABS(X0(COMP, POINT)) * RELAT + ABSOL X0(COMP, POINT) = X0(COMP, POINT) + PERTRB(INDEX) 0100 CONTINUE C C/// CALL THE FUNCTION AT THE PERTRBED X0 AND STORE THE RESULT IN Y1. C CALL FUN (KK, POINTS, COMPS, LENRGY, LMULTI, LTDIF, .FALSE., 1 LTIME, P, WT, YOXID, YFUEL, DT, NTEMP, XGIVEN, 2 TGIVEN, MESH, SN, X0, WNDFAC, SCRTCH(1, 1), 3 SCRTCH(1,5), YV, SCRTCH(1, 2), SCRTCH(1, 3), 4 SCRTCH(1, 4), TFUEL, TOXID, VFUEL, VOXID, 4 AFUEL, AOXID, 5 VISC, COND, D, DKJ, TDR, ICKWRK, RCKWRK, 6 IMCWRK, RMCWRK, Y1, RFAC) C C/// RESTORE X0 TO ITS ORIGINAL VALUE. C INDEX = 0 DO 0200 POINT = CLASS, POINTS, MODULS INDEX = INDEX + 1 X0(COMP, POINT) = SAVE(INDEX) 0200 CONTINUE C C/// DIFFERENCE TO GET THE COLUMNS OF THE JACOBIAN. C INDEX = 0 DO 0300 POINT = CLASS, POINTS, MODULS INDEX = INDEX + 1 TEMP = 1.0 / PERTRB(INDEX) K = COMP + (POINT - 1) * COMPS DO 0300 J2 = MAX (POINT - LOWER, 1), + MIN (POINT + UPPER, POINTS) JOFSET = (J2 - 1) * COMPS - K + KOFSET DO 0300 J1 = 1, COMPS A(J1 + JOFSET, K) = (Y1(J1, J2) - Y0(J1, J2)) * TEMP c INUM=J1 + JOFSET c WRITE(15,1234) INUM,K,A(J1 + JOFSET, K) 0300 CONTINUE C C/// BOTTOM OF THE LOOPS OVER THE RESIDUE CLASSES AND SOLUTION C/// COMPONENTS. C 0400 CONTINUE C C WRITE(15,2345) C2345 FORMAT(/2X, 'END') C1234 FORMAT(1X,I5,2X,I5,2X,E10.3) RETURN END C C--------------------------------------------------------------------- C SUBROUTINE FUN (KK, JJ, NATJ, LENRGY, LMULTI, LTDIF, LVARMC, 1 LTIME, P, WT, YOXID, YFUEL, DT, NTEMP, XGIVEN, 2 TGIVEN, X, SN, S, WNDFAC, YAV, XAV, YV, WDOT, CP, 3 H, TFUEL, TOXID, VFUEL, VOXID, 3 AFUEL, AOXID, VISC, 4 COND, D, DKJ, TDR, ICKWRK, RCKWRK, IMCWRK, RMCWRK, 5 F, RFAC) C C This subroutine forms the function for the differential equations C Indicies: C NT : Temperature C NG : Velocity gradient C NF : Velocity C NH : Eigenvalue, pressure gradient C KK : no. of species C JJ : no. of gridpoints C NATJ : no. variables at each point C LENERGY = .TRUE. if solve energy equation C LMULTI C LTDIFF C LVARMC C LTIME C VFUEL, VOXID, TFUEL, TOXID, YFUEL, YOXID... C C--------------------------------------------------------------------- C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C COMMON /LOCS/ NT, NG, NF, NH,NTH ,NYS, NY C DIMENSION WT(*), YOXID(*), YFUEL(*), XGIVEN(*), TGIVEN(*), X(*), 1 S(NATJ,*), SN(NATJ,*), YAV(*), XAV(*), YV(KK,*), 2 WDOT(*), CP(*), H(*), VISC(*), COND(*), D(KK,*), 3 TDR(KK,*), DKJ(KK,KK,*), ICKWRK(*), RCKWRK(*), 4 IMCWRK(*), RMCWRK(*), F(NATJ,*), C Added by S. R. Lee (06/03/99) 5 XM(KK), XXM(KK), FIX(KK),QM(KK),Q(179),C(KK) C LOGICAL LENRGY, LMULTI, LTDIF, LVARMC, LTIME, JIA COMMON/QIAO/ JIA COMMON/JI/NUMJIA,NKKK COMMON/HUI/OFIX ,OX, OM COMMON/FUSHE/IH2O,ICO2 C C C timelimit C C*****unicos timelimit C call tremain (secs) C if (secs .lt. 60.0) stop C*****END unicos timelimit C C EVALUATE AND STORE THE TRANSPORT COEFFICIENTS C IF (LVARMC .AND. .NOT.LTIME) THEN CALL MTRNPR (KK, JJ, NATJ, LENRGY, LMULTI, LTDIF, P, X, S, 1 YAV, ICKWRK, RCKWRK, IMCWRK, RMCWRK, WT, 2 H, CP, XAV, VISC, COND, D, TDR, DKJ) ENDIF C C EVALUATE AND STORE THE DIFFUSION VELOCITIES C (in the following call h(*) and cp(*) are used C temporarily for scratch space) C CALL MDIFV (KK, JJ, NATJ, LMULTI, LTDIF, X, S, WT, YAV, 1 H, CP, P, D, TDR, ICKWRK, RCKWRK, YV, DKJ) C C FUEL BOUNDARY (X=0) C TAV = 0.5 * (S(NT,1) + S(NT,2)) DO 10 K = 1, KK YAV(K) = 0.5 * (S(NYS+K,1) + S(NYS+K,2)) 10 CONTINUE CALL CKRHOY (P, TAV, YAV, ICKWRK, RCKWRK, RHOP) C DO 30 K = 1, KK c F(NYS+K, 1) = RHOP*YFUEL(K)*VFUEL- c 2 RHOP*S(NYS+K,1)*VFUEL-RHOP*YV(K,1) F(NYS+K, 1) = S(NYS+K, 1) - YFUEL(K) 30 CONTINUE IF (LENRGY) THEN F(NT,1) = TFUEL - S(NT,1) ELSE CALL TEMP (NTEMP, X(1), XGIVEN, TGIVEN, TI) F(NT, 1) = S(NT,1) - TI ENDIF F(NF,1) = S(NF,1) - VFUEL*RHOP/2. F(NG,1) = AFUEL*RHOP - S(NG,1) F(NH,1) = S(NH,2) - S(NH,1) F(NTH,1) = S(NTH,2) - S(NTH,1) C C INTERIOR MESH POINTS C DO 1000 J = 2, JJ-1 C TAV = 0.5 * (S(NT,J) + S(NT,J+1)) DO 100 K = 1, KK YAV(K) = 0.5 * (S(NYS+K,J) + S(NYS+K,J+1)) 100 CONTINUE C RHOM = RHOP CALL CKRHOY (P, TAV, YAV, ICKWRK, RCKWRK, RHOP) CALL CKRHOY (P, S(NT,J), S(NY,J), ICKWRK, RCKWRK, RHOJ) CALL CKRHOY (P, S(NT,J+1), S(NY,J+1), ICKWRK, RCKWRK, RHOJP1) CALL CKRHOY (P, S(NT,J-1), S(NY,J-1), ICKWRK, RCKWRK, RHOJM1) C C FORM THE CHEMICAL RATE TERMS C CALL CKWYP (P, S(NT,J), S(NY,J), ICKWRK, RCKWRK, WDOT) CALL CKHML (S(NT,J), ICKWRK, RCKWRK, H) CALL CKCPBS (S(NT,J), S(NY,J), ICKWRK, RCKWRK, CPB) CALL CKCPMS (S(NT,J), ICKWRK, RCKWRK, CP) C C ACCELERATE RATES TO FACILITATE SOLUTION? C DO 120 K = 1, KK WDOT(K) = WDOT(K) *RFAC 120 CONTINUE C C FORM THE MESH DIFFERENCES C DXP = (X(J+1) - X(J) ) DXM = (X(J) - X(J-1)) DXAV = 0.5 * (X(J+1) - X(J-1)) DXPM = (X(J+1) - X(J-1)) C C FORM THE COEFFICIENTS FOR CENTRAL DIFFERENCES C CENDFM = - DXP / (DXM*DXPM) CENDFC = (DXP-DXM) / (DXP*DXM) CENDFP = DXM / (DXP*DXPM) C C WINDWARD DIFFERENCING FACTORS C IF (WNDFAC .EQ. 1.0) THEN IF (S(NF,J) .LE. 0.) THEN WINDA = 1. WINDB = -1. WINDC = 0. ELSE WINDA = 0. WINDB = 1. WINDC = -1. ENDIF DXWIND = WINDA*X(J+1) + WINDB*X(J) + WINDC*X(J-1) ENDIF C C SPECIES CONSERVATION EQUATION C IF (WNDFAC .EQ. 1.0) THEN DO 200 K = 1, KK DYK = WINDA*S(NYS+K,J+1) + WINDB*S(NYS+K,J) 1 + WINDC*S(NYS+K,J-1) F(NYS+K,J) = 2.0*S(NF,J) * DYK / DXWIND 1 - WDOT(K) * WT(K) 2 + (RHOP*YV(K,J) - RHOM*YV(K,J-1)) / DXAV 200 CONTINUE ELSE DO 210 K = 1, KK F(NYS+K,J) = 2.0*S(NF,J) * 1 (CENDFP*S(NYS+K,J+1) + CENDFC*S(NYS+K,J) + 1 CENDFM*S(NYS+K,J-1) ) - 2 WDOT(K) * WT(K) + 3 (RHOP*YV(K,J) - RHOM*YV(K,J-1)) / DXAV 210 CONTINUE ENDIF C C STREAMFUNCTION EQUATION C F(NF,J) = (S(NG,J)+S(NG,J-1))/2. - (S(NF,J) - S(NF,J-1))/DXM C C MOMENTUM EQUATION C BIGCP = VISC(J) BIGCM = VISC(J-1) IF (WNDFAC .EQ. 1.0) THEN DUGR = WINDA*( S(NF,J+1)*S(NG,J+1)/RHOJP1 ) 1 + WINDB*( S(NF,J )*S(NG,J )/RHOJ ) 2 + WINDC*( S(NF,J-1)*S(NG,J-1)/RHOJM1 ) F(NG,J) = - 3.0*S(NG,J)**2/RHOJ + 2 2.0 * DUGR / DXWIND 3 - ( BIGCP * (S(NG,J+1)/RHOJP1 - S(NG,J)/RHOJ) / DXP - 3 BIGCM * (S(NG,J)/RHOJ -S(NG,J-1)/RHOJM1)/DXM ) / 3 DXAV - S(NH,J) ELSE F(NG,J) = - 3.0*S(NG,J)**2/RHOJ + 2 2.0 * (CENDFP*S(NG,J+1)*S(NF,J+1)/RHOJP1 + 2 CENDFC*S(NG,J)*S(NF,J)/RHOJ + 2 CENDFM*S(NG,J-1)*S(NF,J-1)/RHOJM1 ) - 3 ( BIGCP * (S(NG,J+1)/RHOJP1 - S(NG,J)/RHOJ) / DXP - 3 BIGCM * (S(NG,J)/RHOJ -S(NG,J-1)/RHOJM1)/DXM ) / 3 DXAV - 4 S(NH,J) ENDIF IF (JIA) THEN F(NTH,J) = S(NTH,J+1) - S(NTH,J) ELSE IF (X(J).LT.OX) THEN F(NTH,J)=S(NTH,J+1)-S(NTH,J) ELSE IF(ABS(X(J)-OX).LT.1E-5)THEN F(NTH,J) =S(NYS+NUMJIA,J) - OFIX ELSE F(NTH,J) = S(NTH,J) - S(NTH,J-1) ENDIF ENDIF C C ENERGY EQUATION C C IF (LENRGY) THEN CALL CKYTX(S(NY,J),ICKWRK, RCKWRK,XM) IF (ICO2 .EQ. 0) THEN XCO2 = 0. ELSE XCO2 = MAX(0., XM(ICO2)) ENDIF IF (IH2O .EQ. 0) THEN XH2O = 0. ELSE XH2O = MAX(0., XM(IH2O)) ENDIF PKH2O = 1.10263E7 *(2.0E-7 - 6.4E-11*S(NT,J)) PKCO2 = 1.10263E7 *(4.6E-7 - 1.3E-10*S(NT,J)) ABSP = PKH2O*XH2O + PKCO2*XCO2 ABSP = ABSP*P/(100*1013250) T4 =S(NT,J)*S(NT,J)*S(NT,J)*S(NT,J) OPQFUN = 4.0 * ABSP * 5.67E-12 * (T4 - OTB)*1.E7 C C*****precision > double TDOT = DDOT (KK, WDOT, 1, H, 1) C*****END precision > double C*****precision > single C TDOT = SDOT (KK, WDOT, 1, H, 1) C*****END precision > single +OPQFUN/CPB C c TDOT=TDOT+2*HRXI*Q(99) SUM = 0.0 DO 400 K = 1, KK SUM = SUM + 0.5 * (YV(K,J) + YV(K,J-1)) * 1 CP(K) * (CENDFP*S(NT,J+1) + CENDFC*S(NT,J) + 2 CENDFM*S(NT,J-1) ) 400 CONTINUE C IF (WNDFAC .EQ. 1.0) THEN DTEMP = WINDA*S(NT,J+1) + WINDB*S(NT,J) + WINDC*S(NT,J-1) F(NT,J) = 1 2.0*S(NF,J) * DTEMP / DXWIND - 1 ( COND(J)*(S(NT,J+1)-S(NT,J))/DXP - 2 COND(J-1)*(S(NT,J)-S(NT,J-1))/DXM ) / 3 (CPB*DXAV) + 4 RHOJ * SUM / CPB + 5 TDOT / CPB+OPQFUN/CPB C +OPQFUN/CPB C==================================================================== ELSE F(NT,J) = 2.0*S(NF,J) * 1 (CENDFP*S(NT,J+1) + CENDFC*S(NT,J) + 1 CENDFM*S(NT,J-1) ) - 1 ( COND(J)*(S(NT,J+1)-S(NT,J))/DXP - 2 COND(J-1)*(S(NT,J)-S(NT,J-1))/DXM ) / 3 (CPB*DXAV) + 4 RHOJ * SUM / CPB + 5 TDOT / CPB+OPQFUN/CPB C +OPQFUN/CPB C==================================================================== ENDIF C ELSE CALL TEMP (NTEMP, X(J), XGIVEN, TGIVEN, TI) F(NT, J) = S(NT,J) - TI ENDIF C C EIGENVALUE EQUATION, H C F(NH,J) = S(NH,J+1) - S(NH,J) C 1000 CONTINUE C C OXIDIZER BOUNDARY (X=L) C CALL CKRHOY (P, S(NT,JJ), S(NY,JJ), ICKWRK, RCKWRK, RHOJ) C IF(JIA) THEN F(NTH,JJ)=TOXID-S(NTH,JJ) OM=sqrt(RHOJ)*VOXID ELSE F(NTH,JJ)=S(NTH,JJ-1)-S(NTH,JJ) VOXID=OM/sqrt(RHOJ) ENDIF DO 1200 K = 1, KK F(NYS+K, JJ) = S(NYS+K,JJ) - YOXID(K) c F(NYS+K, JJ) = YOXID(K)*RHOJ*VOXID c 1 - RHOJ*VOXID*S(NYS+K,JJ) - RHOJ*YV(K,JJ-1) 1200 CONTINUE C IF (LENRGY) THEN F(NT, JJ) = S(NT,JJ) - S(NTH,JJ) ELSE CALL TEMP (NTEMP, X(JJ), XGIVEN, TGIVEN, TI) F(NT, JJ) = S(NT,JJ) - TI ENDIF C F(NH,JJ) = RHOJ*VOXID/2.0 - S(NF,JJ) F(NG,JJ) = S(NG,JJ) - AOXID*RHOJ F(NF,JJ) = (S(NG,JJ)+S(NG,JJ-1))/2. - (S(NF,JJ) - S(NF,JJ-1))/DXP C C ADD THE TIME STEP, IF NEEDED C IF (.NOT. LTIME) RETURN C DO 2500 J = 2, JJ-1 C CALL CKRHOY (P, S(NT,J), S(NY,J), ICKWRK, RCKWRK, RHO) C DO 2200 K = 1, KK DYDT = (S(NYS+K,J) - SN(NYS+K,J)) / DT F(NYS+K,J) = F(NYS+K,J) + RHO*DYDT 2200 CONTINUE DGDT = (S(NG,J) - SN(NG,J)) / DT F(NG,J) = F(NG,J) + DGDT C IF (LENRGY) THEN DTDT = (S(NT,J) - SN(NT,J)) / DT F(NT,J) = F(NT,J) + RHO*DTDT ENDIF 2500 CONTINUE C RETURN END C-------------------------------------------------------------------- SUBROUTINE MTRNPR (KK, JJ, NATJ, LENRGY, LMULTI, LTDIF, P, X, S, 1 YAV, ICKWRK, RCKWRK, IMCWRK, RMCWRK, WT, 2 XMF, XMFP, XAV, VISC, COND, D, TDR, DKJ) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C COMMON /LOCS/ NT, NG, NF, NH,NTH ,NYS, NY C DIMENSION X(*), S(NATJ,*), YAV(*), XAV(*), XMF(*),WT(*), 1 XMFP(*), ICKWRK(*), RCKWRK(*), IMCWRK(*), RMCWRK(*), 2 VISC(*), COND(*), DKJ(KK,KK,*), D(KK,*), TDR(KK,*) LOGICAL LENRGY, LTDIF, LMULTI C DATA EPS/1.0E-30/ C IF (LMULTI) THEN C-----------------MULTICOMPONENT FORMULAS--------------------- CALL CKYTX (S(NY,1), ICKWRK, RCKWRK, XMFP) DO 200 J = 1, JJ-1 TAV = 0.5 * (S(NT,J) + S(NT,J+1)) C C DIMENSIONAL TEMPERATURE AT THE GRID POINTS C DO 100 K = 1, KK YAV(K) = 0.5 * (S(NYS+K,J) + S(NYS+K,J+1)) XMF(K) = XMFP(K) 100 CONTINUE CALL CKYTX (YAV, ICKWRK, RCKWRK, XAV) CALL CKYTX (S(NY,J+1), ICKWRK, RCKWRK, XMFP) CALL CKMMWX (XAV, ICKWRK, RCKWRK, WTMAV) CALL MCMDIF (P, TAV, XAV, KK, IMCWRK, RMCWRK, DKJ(1,1,J)) CALL MCAVIS (TAV, XAV, RMCWRK, VISC(J)) DO 75 K = 1, KK SUMN = 0.0 DO 50 L = 1, KK CC SUMN = SUMN + WT(L) * DKJ(K,L,J) * CC 1 (XMFP(L) - XMF(L)) / (X(J+1) - X(J)) SUMN = SUMN + DKJ(K,L,J) * 1 (S(NYS+L,J+1) - S(NYS+L,J)) / (X(J+1) - X(J)) 50 CONTINUE CC DENOM = - (XMFP(K) - XMF(K)) / (X(J+1) - X(J)) DENOM = - (S(NYS+K,J+1) - S(NYS+K,J)) / (X(J+1) - X(J)) D(K,J) = (SUMN + EPS) / ( WTMAV * (DENOM + EPS)) 75 CONTINUE C C DETERMINE THE MIXTURE CONDUCTIVITY AND C THERMAL DIFFUSION COEFFICIENT AT J C IF (LENRGY .OR. LTDIF) 1 CALL MCMCDT (P, TAV, XAV, IMCWRK, RMCWRK, 2 ICKWRK, RCKWRK, TDR(1,J), COND(J)) C 200 CONTINUE c RETURN c ENDIF ELSE C-----------------MIXTURE-AVERAGED FORMULAS--------------------- C DO 400 J = 1, JJ-1 TAV = 0.5 * (S(NT,J) + S(NT,J+1)) C C DIMENSIONAL TEMPERATURE AT THE GRID POINTS C DO 300 K = 1, KK YAV(K) = 0.5 * (S(NYS+K,J) + S(NYS+K,J+1)) 300 CONTINUE CALL CKYTX (YAV, ICKWRK, RCKWRK, XAV) CALL MCADIF (P, TAV, XAV, RMCWRK, D(1,J)) CALL MCAVIS (TAV, XAV, RMCWRK, VISC(J)) C C DETERMINE THE MIXTURE CONDUCTIVITY AT J C IF (LENRGY) CALL MCACON (TAV, XAV, RMCWRK, COND(J)) IF (LTDIF)THEN CALL MCATDR (TAV, XAV, IMCWRK, RMCWRK, TDR(1,J)) CALL CKRHOY (P, TAV, YAV, ICKWRK, RCKWRK, RHOAV) CALL CKMMWY (YAV, ICKWRK, RCKWRK, WTM) DO 350 K = 1, KK TDR(K,J) = D(K,J) * TDR(K,J) * RHOAV * WT(K)/WTM 350 CONTINUE ENDIF C 400 CONTINUE C C ENDIF RETURN END C-------------------------------------------------------------------- SUBROUTINE MDIFV (KK, JJ, NATJ, LMULTI, LTDIF, X, S, WT, YAV, 1 XMF, XMFP, P, D, TDR, ICKWRK, RCKWRK, YV, DKJ) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C COMMON /LOCS/ NT, NG, NF, NH,NTH ,NYS, NY C DIMENSION X(*), S(NATJ,*), WT(*), YAV(*), XMF(*), XMFP(*), 1 D(KK,*), TDR(KK,*), DKJ(KK,KK,*), YV(KK,*), 2 ICKWRK(*), RCKWRK(*) LOGICAL LTDIF, LDIRCT, LMULTI C LDIRCT = .TRUE. C CALL CKYTX (S(NY,1), ICKWRK, RCKWRK, XMFP) C C LOOP OVER ALL MESH POINTS, COMPUTING THE DIFFUSION C VELOCITY AT THE MID POINTS. THE INDEXING IS SUCH THAT C YV(K,J) IS THE DIFFUSION VELOCITY OF THE KTH SPECIES C MIDWAY BETWEEN NODES J AND J+1. C DO 1000 J = 1, JJ-1 C TMF = S(NT, J) TMFP = S(NT, J+1) TAV = 0.5 * (TMF + TMFP) DO 300 K = 1, KK YAV(K) = 0.5 * (S(NYS+K,J) + S(NYS+K,J+1)) XMF(K) = XMFP(K) 300 CONTINUE C CALL CKMMWY (YAV, ICKWRK, RCKWRK, WTM) CALL CKRHOY (P, TAV, YAV, ICKWRK, RCKWRK, RHOAV) CALL CKYTX (S(NY,J+1), ICKWRK, RCKWRK, XMFP) C IF (LMULTI. AND. LDIRCT) THEN C EVALUATE THE MULTICOMPONENT DIFFUSION VELOCITY DIRECTLY, C RATHER THAN USE THE MIXTURE-AVERAGED FORM FOR D(K,J) DO 475 K = 1, KK SUM = 0.0 DO 450 L = 1, KK SUM = SUM + WT(L) * DKJ(K,L,J) * 1 (XMFP(L)-XMF(L)) / (X(J+1)-X(J)) 450 CONTINUE YV(K,J) = (WT(K)/WTM**2) * SUM 475 CONTINUE ELSE C USE MIXTURE-AVERAGED FORM FOR FICKIAN DIFFUSION, C WHETHER WE ARE USING THE MULTICOMPONENT FORMALISM C OR MIXTURE-AVERAGED DO 500 K = 1, KK YV(K,J) = - D(K,J) * (WT(K)/WTM) * 1 (XMFP(K)-XMF(K)) / (X(J+1)-X(J)) 500 CONTINUE ENDIF C C ADD THE THERMAL DIFFUSION, IF REQUESTED C IF (LTDIF) THEN C DO 600 K = 1, KK YV(K,J) = YV(K,J) - 1 (TDR(K,J) / (TAV*RHOAV)) * 2 ( TMFP-TMF )/ (X(J+1)-X(J)) 600 CONTINUE C ENDIF C C ADD CORRECTION VELOCITY TO DIFFUSION C SUM = 0. DO 700 K = 1, KK SUM = SUM + YV(K,J) 700 CONTINUE VC = - SUM DO 800 K = 1, KK YV(K,J) = YV(K,J) + YAV(K)*VC 800 CONTINUE 1000 CONTINUE C RETURN END C-------------------------------------------------------------------- SUBROUTINE PRINT (LOUT, KK, JJ, NATJ, P, X, S, XMF, KSYM, 1 ICKWRK, RCKWRK, KW, NW,LREST) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C COMMON /LOCS/ NT, NG, NF, NH,NTH ,NYS, NY C DIMENSION X(*), S(NATJ, *), XMF(*), ICKWRK(*), RCKWRK(*), 1 KW(*) CHARACTER KSYM(*)*(*) C C DATA KPERLN /6/ C C LS = 1 K2 = MIN (KPERLN, KK) C C C PRINT THE FLOW VARIABLES C C WRITE (LOUT, 7070) WRITE (LREST, 7070) WRITE (LREST, 7021)S(NTH,1) DO 100 J = 1, JJ CALL CKRHOY (P, S(NT,J), S(NY,J), ICKWRK, RCKWRK, RHO) VEL = 2.0*S(NF,J)/RHO C WRITE (LOUT, 7020) X(J), S(NF,J), VEL, S(NG,J), S(NTH,J), C 1 S(NT,J), RHO WRITE (LREST, 7020) X(J), S(NF,J), VEL, S(NG,J), S(NTH,J), 1 S(NT,J), RHO 100 CONTINUE C C PRINT THE MOLE FRACTIONS C IF (NW .EQ. 0) THEN C DO 200 L = LS, KK, KPERLN K1 = L K2 = L+KPERLN-1 K2 = MIN (K2, KK) C C WRITE (LOUT, 7060) (KSYM(K), K=K1,K2) WRITE (LREST, 7060) (KSYM(K), K=K1,K2) DO 200 J = 1, JJ CALL CKYTX (S(NY,J), ICKWRK, RCKWRK, XMF) C WRITE (LOUT, 7020) X(J), (XMF(K), K=K1,K2) WRITE (LREST, 7020) X(J), (XMF(K), K=K1,K2) 200 CONTINUE C ELSE C C WRITE (LOUT, 7060) (KSYM(KW(K)), K=1,NW) WRITE (LREST, 7060) (KSYM(KW(K)), K=1,NW) DO 210 J = 1, JJ CALL CKYTX (S(NY,J), ICKWRK, RCKWRK, XMF) C WRITE (LOUT, 7020) X(J), (S(NYS+KW(K),J), K=1,NW) C WRITE (LOUT, 7020) X(J), (XMF(KW(K)), K=1, NW) WRITE (LREST, 7020) X(J), (XMF(KW(K)), K=1, NW) 210 CONTINUE ENDIF C RETURN C 7020 FORMAT(1X, F6.3, 6(1PE13.5)) 7021 FORMAT(1X, F10.4) 7060 FORMAT(/2X, 'X(cm)' , 3X, 10(A12, 1X)) 7070 FORMAT(/2X, 'X(cm)', 5X, 'F', 9X, 'V(cm/s)', 9X, 'G', 1 10X, 'TH', 12X, 'T(K)', 9X, 'RHO') C END C C--------------------------------------------------------------------- C SUBROUTINE START (KK, NMAX, NATJ, JJ, LOUT, LUMESH, IPROFL, TMAX, 1 TFUEL, TOXID, VFUEL, VOXID, P, 2 FUEL, PROD, OXID, NTEMP, XGIVEN, TGIVEN, 3 XCEN, XEND, WMIX, ICKWRK, RCKWRK, Y, 4 YOXID, YFUEL, X, S) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C COMMON /LOCS/ NT, NG, NF, NH,NTH ,NYS, NY C DIMENSION FUEL(*), OXID(*), PROD(*), 1 ICKWRK(*), RCKWRK(*), Y(*), YOXID(*), YFUEL(*), 2 XGIVEN(*), TGIVEN(*), X(*), S(NATJ, *) C LOGICAL LUMESH C C INITIALIZE MASS FLUX FRACTIONS AND MASS FRACTIONS TO 0. C DO 100 K = 1, KK YOXID(K) = 0.0 YFUEL(K) = 0.0 DO 100 J = 1, NMAX S(NYS+K, J) = 0.0 100 CONTINUE C C SET THE MASS FLUX FRACTION BOUNDARY CONDITIONS C CALL CKXTY (OXID, ICKWRK, RCKWRK, YOXID) CALL CKXTY (FUEL, ICKWRK, RCKWRK, YFUEL) CALL CKXTY (PROD, ICKWRK, RCKWRK, Y ) C C COMPUTE STAGNATION PLANE LOCATION C CALL CKRHOY (P, TFUEL, YFUEL, ICKWRK, RCKWRK, RHOF) CALL CKRHOY (P, TOXID, YOXID, ICKWRK, RCKWRK, RHOO) XSTAG = XEND /( 1. + RHOO/RHOF*( (VOXID/VFUEL)**2 ) ) WRITE (LOUT,'(/1X,A,1PE10.3/)') 1 'ESTIMATED STAGNATION LOCATION =', XSTAG IF (XCEN .EQ. 0.) THEN XCEN = XSTAG IF (LUMESH) THEN WMIX = 2.*MIN( XEND-XCEN, XCEN ) ELSE WMIX = 2.*MIN( X(JJ)-XCEN, XCEN-X(1) ) ENDIF ENDIF C IF (LUMESH) THEN C C SET UNIFORM X MESH COORDINATES C DX = (XEND-0.0) / FLOAT(JJ-1) DO 200 J = 1, JJ X(J) = 0.0 + DX*FLOAT(J-1) 200 CONTINUE ENDIF C IF (IPROFL .EQ. 1) THEN C C SET TRIANGLE PROFILES C DO 800 K = 1,KK DO 800 J = 1,JJ CALL TRNGL (WMIX, XCEN, YFUEL(K), Y(K), YOXID(K), 1 X(J), S(NYS+K,J) ) 800 CONTINUE C C SET THE TEMPERATURE PROFILE C DO 850 J = 1, JJ CALL TRNGL (WMIX, XCEN, TFUEL, TMAX, TOXID, 1 X(J), S(NT,J) ) 850 CONTINUE C ELSEIF (IPROFL .EQ. 2) THEN C C SET LINEAR PROFILES C DO 900 K = 1,KK DO 900 J = 1,JJ CALL LINE (WMIX, XCEN, YFUEL(K), YOXID(K), 1 X(J), S(NYS+K,J) ) 900 CONTINUE DO 950 J = 1, JJ S(NT,J) = TFUEL 950 CONTINUE C C ELSEIF (IPROFL .EQ. 3) THEN C C SET PLATEAU PROFILES C DO 960 K = 1,KK DO 960 J = 1,JJ CALL PLATOW (WMIX, XCEN, YFUEL(K), Y(K), YOXID(K), 1 X(1), X(J), X(JJ), S(NYS+K,J) ) 960 CONTINUE DO 970 J = 1,JJ CALL PLATOW (WMIX, XCEN, TFUEL, TMAX, TOXID, 1 X(1), X(J), X(JJ), S(NT,J) ) 970 CONTINUE C ENDIF C C CONVERT ALL STARTING ESTIMATES TO MASS FRACTION C NO NEED TO DO THIS IF SET PROFILE IN Y FIRST! CC C DO 1600 J = 1,JJ C CALL CKXTY (S(NY, J), ICKWRK, RCKWRK, Y) C DO 1600 K = 1, KK C S(NYS+K, J) = Y(K) C1600 CONTINUE C C SET THE 'GIVEN' TEMPERATURES TO THE COMPUTED PROFILE C NTEMP = JJ DO 1770 N = 1, NTEMP XGIVEN(N) = X(N) TGIVEN(N) = S(NT,N) 1770 CONTINUE C C CALL CKRHOY (P, S(NT,1 ), S(NY,1 ), ICKWRK, RCKWRK, RHOF) C CALL CKRHOY (P, S(NT,JJ), S(NY,JJ), ICKWRK, RCKWRK, RHOO) FFUEL = RHOF*VFUEL/2. FOXID = RHOO*VOXID/2. DFDX = (FOXID - FFUEL)/(X(JJ) - X(1)) DO 1800 J = 1, JJ C C SET A LINEAR F PROFILE C S(NF,J) = DFDX* (X(J)-X(1)) + FFUEL C C SET G AS CONSTANT (TO MATCH LINEAR PROFILE IN F) C S(NG,J) = DFDX C C SET THE PRESSURE GRADIENT EIGENVALUE C S(NH,J) = -100. S(NTH,J)= TOXID C 1800 CONTINUE C RETURN END C C------------------------------------------------------------------- C SUBROUTINE TRNGL (WMIX, XCEN, FLOORL, FMAX, FLOORR, X, FVAL) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C IF (X .LE. (XCEN-WMIX/2.) ) THEN FVAL = FLOORL RETURN ELSE IF ((X .GT. (XCEN-WMIX/2.)) .AND. 1 (X .LE. XCEN) ) THEN FVAL = (FMAX-FLOORL)/(0.5*WMIX) * (X-XCEN+WMIX*0.5) + FLOORL RETURN ELSE IF ((X .GT. XCEN) .AND. 1 (X .LE. (XCEN+WMIX/2.)) ) THEN FVAL = (FLOORR-FMAX)/(0.5*WMIX) * (X-XCEN) + FMAX RETURN ELSE IF (X .GT. (XCEN+WMIX/2.) ) THEN FVAL = FLOORR RETURN ENDIF C WRITE (6,*) ' ERROR IN SUBROUTINE TRNGL...' STOP C END C C------------------------------------------------------------------- C SUBROUTINE LINE (WMIX, XCEN, FLOORL, FLOORR, X, FVAL) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C IF (X .LE. (XCEN-WMIX/2.) ) THEN FVAL = FLOORL RETURN ELSE IF ( X .LT. (XCEN+WMIX/2.) ) THEN FVAL = (FLOORR-FLOORL)/WMIX * (X-XCEN+WMIX*0.5) + FLOORL RETURN ELSE FVAL = FLOORR RETURN ENDIF END C C------------------------------------------------------------------- C SUBROUTINE PLATOW (W, XC, FL, FP, FR, X0, X, XL, F) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C IF (X .LE. (XC-W/2.) ) THEN F = FL + (FP-FL)*(X-X0)/(XC-W/2.-X0) RETURN ELSE IF ( X .LT. (XC+W/2.) ) THEN F = FP RETURN ELSE IF ( X .LE. XL ) THEN F = FP + (FR-FP)*(X-(XC+W/2.))/(XL-(XC+W/2.)) RETURN ELSE WRITE (6,*) 'STOP, X NOT FOUND IN -PLATOW-' STOP ENDIF END C C------------------------------------------------------------------- C SUBROUTINE RESTRT (KK, NMAX, NATJ, JJ, LOUT, LUSTGV, 1 FUEL, OXID, NTEMP, XGIVEN, TGIVEN, XEND, 2 ICKWRK, RCKWRK, YOXID, YFUEL, X, S) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C COMMON /LOCS/ NT, NG, NF, NH,NTH ,NYS, NY C DIMENSION FUEL(*), OXID(*), 1 ICKWRK(*), RCKWRK(*), YOXID(*), YFUEL(*), 2 XGIVEN(*), TGIVEN(*), X(*), S(NATJ, *) C LOGICAL LUSTGV C C C INITIALIZE MASS FLUX FRACTIONS TO 0. C DO 100 K = 1, KK YOXID(K) = 0.0 YFUEL(K) = 0.0 100 CONTINUE C C IF A NEW XEND .GT. X(JJ), THEN ADD A POINT AT JJ+1 C IF A NEW XEND .LE. X(JJ), THEN REDUCE JJ, AND SET X(JJ)=XEND C IF (XEND .GT. (X(JJ)+1.E-4) ) THEN JJ = JJ + 1 IF (JJ .GT. NMAX) THEN WRITE (LOUT, *) ' ERROR...NEW XEND NEEDS TOO MANY POINTS' STOP ENDIF X(JJ) = XEND DO 400 N = 1, NATJ S(N,JJ) = S(N,JJ-1) 400 CONTINUE ELSE C DO 500 J = 1, JJ C IF (XEND .GE. X(J)) THEN C X(J) = XEND C GO TO 550 C ENDIF C500 CONTINUE C550 CONTINUE C JJ = J ENDIF C C SET THE MASS FLUX FRACTION BOUNDARY CONDITIONS C CALL CKXTY (FUEL, ICKWRK, RCKWRK, YFUEL) CALL CKXTY (OXID, ICKWRK, RCKWRK, YOXID) C C SET XGIVEN AND TGIVEN TO THE OLD SOLUTION C IF (.NOT. LUSTGV) THEN NTEMP = JJ DO 1250 N = 1, NTEMP XGIVEN(N) = X(N) TGIVEN(N) = S(NT,N) 1250 CONTINUE ENDIF C C RETURN END C C-------------------------------------------------------------------- C SUBROUTINE RDKEY (KK, NMAX, LIN, LOUT, KSYM, PATM, LTIME, 1 LTIME2, LUSTGV, LENRGY, LEQUIL, LMULTI, LTDIF, 2 LUMESH, LRSTRT, LCNTUE, IPROFL, IPRNT, MFILE, 3 UFAC, DFAC, RATEF, N1CALL, LSEN, VFUEL, VOXID, 3 AFUEL, AOXID, TFUEL, TOXID, 4 TMAX, P, NPTS, X, FUEL, 5 PROD, OXID, XCEN, XEND, WMIX, GRAD, 6 CURV, NADP, IRETIR, SFLR, NTEMP, XX, TT, ATOL, 7 RTOL, ATIM, RTIM, NUMDT, DT, NUMDT2, DT2, DTMIN, 8 WNDFAC, LREGRD, JJREGD, PCTADP, RATGTC, 9 KW, NW, LFXTEQ, NJAC, NINIT, ITJAC, DTMAX ) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C COMMON /LOCS/ NT, NG, NF, NH,NTH ,NYS, NY C DIMENSION FUEL(*), OXID(*), PROD(*), 1 XX(*), TT(*), X(*), VALUE(5), KW(*) C CHARACTER KEYWRD*4, KSYM(*)*(*), LINE*80 C LOGICAL LTIME, LTIME2, LENRGY, LMULTI, LTDIF, LUMESH, LRSTRT, 1 LCNTUE, LUSTGV, LSEN, NEC(8), NOPT(5), CNTNUD, LFIRST, 2 IERR, KERR, LREGRD, LEQUIL, LFXTEQ C C DATA NEC/8*.FALSE./, NOPT/5*.FALSE./ DATA CNTNUD/.FALSE./ C C INITIALIZE VARIABLES C KERR = .FALSE. IF (LCNTUE) THEN C LCNTUE = .FALSE. CNTNUD = .TRUE. LFIRST = .TRUE. NP = 0 C ELSE C DO 10 K = 1, KK FUEL(K) = 0. OXID(K) = 0. PROD(K) = 0. 10 CONTINUE DO 11 K = 1, 6 KW(K) = 0 11 CONTINUE C NPTS = 6 GRAD = 0.1 CURV = 0.5 NADP = 3 IRETIR = 50 ATOL = 1.0E-9 RTOL = 1.0E-4 ATIM = 1.0E-9 RTIM = 1.0E-4 SFLR = -5.E-4 NUMDT = 100 DT = 1.0E-6 NUMDT2= 100 DT2 = 1.0E-6 DTMIN = 1.E-10 NTEMP = 0 NP = 0 NW = 0 WNDFAC= 1.0 VFUEL = 0. VOXID = 0. AFUEL = 0. AOXID = 0. MFILE = 1 UFAC = 2.0 DFAC = 2.2 N1CALL= 1 P = PATM TFUEL = 300. TOXID = 300. TMAX = 2200. JJREGD= 40 PCTADP= .75 RATGTC= 1.0 IPRNT = 1 RATEF = 1.0 NJAC = 20 NINIT = 0 ITJAC = 20 DTMAX = 1.E-4 LFIRST= .TRUE. LUMESH= .TRUE. LUSTGV= .FALSE. LCNTUE= .FALSE. LRSTRT= .FALSE. LTDIF = .FALSE. LMULTI= .FALSE. LTIME = .FALSE. LTIME2= .FALSE. LENRGY= .FALSE. LSEN = .FALSE. LREGRD= .FALSE. IPROFL= 1 LEQUIL= .FALSE. LFXTEQ= .FALSE. ENDIF C C-------------------------------------------------------------- C C READ NEXT INPUT LINE C WRITE (LOUT,'(/A/)') ' KEYWORD INPUT ' C 90 CONTINUE KEYWRD = ' ' LINE = ' ' READ (LIN, 7000) KEYWRD, LINE WRITE (LOUT, 8000) KEYWRD, LINE C C IS THIS A KEYWORD COMMENT? C IF (KEYWRD(1:1) .EQ. '.' .OR. KEYWRD(1:1) .EQ. '/' 1 .OR. KEYWRD(1:1) .EQ. '!') GO TO 90 C CALL UPCASE (KEYWRD) C C--------------PROBLEM TYPE KEYWORDS-------------------- C C ENERGY EQUATION IS NOT INCLUDED C IF (KEYWRD .EQ. 'TGIV') THEN LENRGY = .FALSE. C C ENERGY EQUATION IS INCLUDED C ELSE IF (KEYWRD .EQ. 'ENRG') THEN LENRGY = .TRUE. C C FACTOR ON REACTION RATES C ELSE IF (KEYWRD .EQ. 'GFAC') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, RATEF, IERR) KERR = KERR.OR.IERR C C USE EQUILIBRIUM SOLUTION ESTIMATE C ELSE IF (KEYWRD .EQ. 'EQUL') THEN LEQUIL = .TRUE. C C FIX TEMP IN EQUIL SOLUTION ESTIMATE C ELSE IF (KEYWRD .EQ. 'FXEQ') THEN LFXTEQ = .TRUE. C C LINEAR PROFILE ON INITIAL SOLUTION C ELSE IF (KEYWRD .EQ. 'LINE') THEN IPROFL = 2 C C PROFILE ON INITIAL SOLUTION IS PLATEAU C ELSE IF (KEYWRD .EQ. 'PLAT') THEN IPROFL = 3 C C--------------METHOD OPTIONS KEYWORDS-------------------- C C C ABSOLUTE NEWTON ITERATION CONVERGENCE CRITERIA C ELSE IF (KEYWRD .EQ. 'ATOL') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, ATOL, IERR) KERR = KERR.OR.IERR C C RELATIVE NEWTON ITERATION CONVERGENCE CRITERIA C ELSE IF (KEYWRD .EQ. 'RTOL') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, RTOL, IERR) KERR = KERR.OR.IERR C C ABSOLUTE NEWTON CONVERGENCE CRITERIA FOR TIMESTEPS C ELSE IF (KEYWRD .EQ. 'ATIM') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, ATIM, IERR) KERR = KERR.OR.IERR C C RELATIVE NEWTON CONVERGENCE CRITERIA FOR TIMESTEPS C ELSE IF (KEYWRD .EQ. 'RTIM') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, RTIM, IERR) KERR = KERR.OR.IERR C C TIME STEP STARTING PROCEDURE C ELSE IF (KEYWRD .EQ. 'TIME') THEN LTIME = .TRUE. CALL CKXNUM (LINE, 2, LOUT, NVAL, VALUE, IERR) KERR = KERR.OR.IERR NUMDT = INT(VALUE(1)) DT = VALUE(2) C C TIME STEPPING, AFTER ADDING THE ENERGY EQUATION C ELSE IF (KEYWRD .EQ. 'TIM2') THEN LTIME2 = .TRUE. CALL CKXNUM (LINE, 2, LOUT, NVAL, VALUE, IERR) KERR = KERR.OR.IERR NUMDT2 = INT(VALUE(1)) DT2 = VALUE(2) C C TIMESTEP INCREASE WHEN TIMESTEP DOES NOT CHANGE SOLUTION C ELSE IF (KEYWRD .EQ. 'UFAC') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, UFAC, IERR) KERR = KERR.OR.IERR C C TIMESTEP DECREASE WHEN NEWTON FAILS CONVERGENCE ON TIMESTEP C ELSE IF (KEYWRD .EQ. 'DFAC') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, DFAC, IERR) KERR = KERR.OR.IERR C C MINIMUM TIMESTEP C ELSE IF (KEYWRD .EQ. 'DTMN') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, DTMIN, IERR) KERR = KERR.OR.IERR C C MAXIMUM TIMESTEP C ELSE IF (KEYWRD .EQ. 'DTMX') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, DTMAX, IERR) KERR = KERR.OR.IERR C C DO NOT DO THE FIXED TEMPERATURE PROBLEM C ELSE IF (KEYWRD .EQ. 'NOFT') THEN N1CALL = 2 C C WINDWARD DIFFERENCING C ELSE IF (KEYWRD .EQ. 'WDIF') THEN WNDFAC = 1.0 C C CENTRAL DIFFERENCING C ELSE IF (KEYWRD .EQ. 'CDIF') THEN WNDFAC = 0.0 C C FLOOR VALUE FOR THE SPECIES BOUNDS C ELSE IF (KEYWRD .EQ. 'SFLR') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, SFLR, IERR) KERR = KERR.OR.IERR C C OXIDIZER C ELSE IF (KEYWRD .EQ. 'KOUT') THEN CALL CKCRAY (LINE, KK, KSYM, LOUT, 6, KW, NFD, IERR) IF (IERR) THEN WRITE (LOUT,'(A)') 1 ' ERROR READING KEYWRD '//KEYWRD KERR = .TRUE. IF (NFD .GT. 6) THEN WRITE (LOUT,'(A,I4,A)') 'FOUND ', NFD,' SPECIES' NW = 6 ENDIF ELSE NW = NFD ENDIF C C--------------GRID PARAMETER KEYWORDS-------------------- C C C NUMBER OF INITIAL MESH POINTS C (THIS IS OVERWRITTEN 'GRID' INPUT) C ELSE IF (KEYWRD .EQ. 'NPTS') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, VALUE, IERR) KERR = KERR.OR.IERR NPTS = INT(VALUE(1)) C C INITIAL MESH C ELSE IF (KEYWRD .EQ. 'GRID') THEN LUMESH = .FALSE. CALL CKXNUM (LINE, 1, LOUT, NVAL, VALUE, IERR) IF (IERR .OR. NP+1.GT.NMAX) THEN KERR = .TRUE. ELSE NP = NP + 1 X(NP) = VALUE(1) ENDIF C C GRADIENT MESH ADAPTION PARAMETER C ELSE IF (KEYWRD .EQ. 'GRAD') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, GRAD, IERR) KERR = KERR.OR.IERR C C CURVATURE MESH ADAPTION PARAMETER C ELSE IF (KEYWRD .EQ. 'CURV') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, CURV, IERR) KERR = KERR.OR.IERR C C NUMBER OF MESH POINTS TO ADD AT ONCE C ELSE IF (KEYWRD .EQ. 'NADP') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, VALUE, IERR) KERR = KERR.OR.IERR NADP = INT(VALUE(1)) C C RETIREMENT PERIOD BEFORE INCREASING THE TIMESTEP C ELSE IF (KEYWRD .EQ. 'IRET') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, VALUE, IERR) KERR = KERR.OR.IERR IRETIR = INT(VALUE(1)) C C INITIAL TIME STEPS BEFORE NEWTON C ELSE IF (KEYWRD .EQ. 'ISTP') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, VALUE, IERR) KERR = KERR.OR.IERR NINIT = INT(VALUE(1)) C C RETIREMENT AGE OF JACOBIAN DURING NEWTON C ELSE IF (KEYWRD .EQ. 'NJAC') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, VALUE, IERR) KERR = KERR.OR.IERR NJAC = INT(VALUE(1)) C C RETIREMENT AGE OF JACOBIAN DURING TIME STEPPING C ELSE IF (KEYWRD .EQ. 'TJAC') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, VALUE, IERR) KERR = KERR.OR.IERR ITJAC = INT(VALUE(1)) C C CENTER OF MIXING REGION C ELSE IF (KEYWRD .EQ. 'XCEN') THEN NOPT(2) = .TRUE. CALL CKXNUM (LINE, 1, LOUT, NVAL, XCEN, IERR) KERR = KERR.OR.IERR C C DISTANCE AT WHICH END BOUNDARY CONDITION IS APPLIED C ELSE IF (KEYWRD .EQ. 'XEND') THEN NEC(1) = .TRUE. CALL CKXNUM (LINE, 1, LOUT, NVAL, XEND, IERR) KERR = KERR.OR.IERR C C WIDTH OF MIXING ZONE C ELSE IF (KEYWRD .EQ. 'WMIX') THEN NOPT(1) = .TRUE. CALL CKXNUM (LINE, 1, LOUT, NVAL, WMIX, IERR) KERR = KERR.OR.IERR C C--------------FLAME DEFINITION KEYWORDS-------------------- C C C REACTANT INLET VELOCITY (r=R) C ELSE IF (KEYWRD .EQ. 'VFUE') THEN NEC(4) = .TRUE. CALL CKXNUM (LINE, 1, LOUT, NVAL, VFUEL, IERR) KERR = KERR.OR.IERR C ELSE IF (KEYWRD .EQ. 'VOXI') THEN C NEC(5) = .TRUE. C CALL CKXNUM (LINE, 1, LOUT, NVAL, VOXID, IERR) C KERR = KERR.OR.IERR C VOXID = - VOXID C C REACTANT INLET VELOCITY SLOPE (r=R) C ELSE IF (KEYWRD .EQ. 'AFUE') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, AFUEL, IERR) KERR = KERR.OR.IERR AFUEL = - AFUEL ELSE IF (KEYWRD .EQ. 'AOXI') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, AOXID, IERR) KERR = KERR.OR.IERR AOXID = - AOXID C C INLET TEMPERATURE C ELSE IF (KEYWRD .EQ. 'TOXI') THEN NEC(2) = .TRUE. CALL CKXNUM (LINE, 1, LOUT, NVAL, TOXID, IERR) KERR = KERR.OR.IERR ELSE IF (KEYWRD .EQ. 'TFUE') THEN NEC(3) = .TRUE. CALL CKXNUM (LINE, 1, LOUT, NVAL, TFUEL, IERR) KERR = KERR.OR.IERR C C MAXIMUM TEMPERATURE C ELSE IF (KEYWRD .EQ. 'TMAX') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, TMAX, IERR) KERR = KERR.OR.IERR C C PRESSURE C ELSE IF (KEYWRD .EQ. 'PRES') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, P, IERR) KERR = KERR.OR.IERR P = P*PATM C C FUEL C ELSE IF (KEYWRD .EQ. 'FUEL') THEN IF (LFIRST) THEN LFIRST = .FALSE. DO 1100 K = 1, KK FUEL(K) = 0. 1100 CONTINUE ENDIF CALL CKSNUM (LINE, 1, LOUT, KSYM, KK, KSPEC, NVAL, 1 VALUE, IERR) IF (IERR) THEN WRITE (LOUT,'(A)') 1 ' ERROR READING KEYWRD '//KEYWRD ELSE FUEL(KSPEC) = VALUE(1) ENDIF C C PRODUCT C ELSE IF (KEYWRD .EQ. 'PROD') THEN CALL CKSNUM (LINE, 1, LOUT, KSYM, KK, KSPEC, NVAL, 1 VALUE, IERR) IF (IERR) THEN WRITE (LOUT,'(A)') 1 ' ERROR READING KEYWRD '//KEYWRD KERR = .TRUE. ELSE PROD(KSPEC) = VALUE(1) ENDIF C C OXIDIZER C ELSE IF (KEYWRD .EQ. 'OXID') THEN CALL CKSNUM (LINE, 1, LOUT, KSYM, KK, KSPEC, NVAL, 1 VALUE, IERR) IF (IERR) THEN WRITE (LOUT,'(A)') 1 ' ERROR READING KEYWRD '//KEYWRD KERR = .TRUE. ELSE OXID(KSPEC) = VALUE(1) ENDIF C C READ SPECIFIED TEMPERATURE PROFILE (X,T) PAIRS C ELSE IF (KEYWRD .EQ. 'TEMP') THEN CALL CKXNUM (LINE, 2, LOUT, NVAL, VALUE, IERR) KERR = KERR.OR.IERR IF (NTEMP+1 .GT. NMAX) THEN WRITE (LOUT, '(A,I4,A)') 1 ' ERROR... THE PROBLEM IS ONLY DIMENSIONED FOR ', 2 ' (X,T) PAIRS' ELSE NTEMP = NTEMP+1 XX(NTEMP) = VALUE(1) TT(NTEMP) = VALUE(2) ENDIF C C ON A RESTART USE GIVEN TEMPERATURE PROFILE, NOT THE ONE ON C THE RESTART FILE C ELSE IF (KEYWRD .EQ. 'USTG') THEN LUSTGV = .TRUE. C C--------------TRANSPORT OPTIONS KEYWORDS-------------------- C C C MULTICOMPONENT FORMULAS USED C ELSE IF (KEYWRD .EQ. 'MULT') THEN LMULTI = .TRUE. C C MIXTURE-AVERAGED FORMULAS USED C ELSE IF (KEYWRD(1:3) .EQ. 'MIX') THEN LMULTI = .FALSE. C C C THERMAL DIFFUSION INCLUDED C ELSE IF (KEYWRD .EQ. 'TDIF') THEN LTDIF = .TRUE. C C--------------SENSITIVITY KEYWORDS-------------------- C C C ALL REACTION SENSITIVITY C ELSE IF (KEYWRD .EQ. 'ASEN') THEN LSEN = .TRUE. C C--------------PRINTING AND RESTARTING KEYWORDS-------------------- C C C PRINT CONTROL C ELSE IF (KEYWRD .EQ. 'PRNT') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, VALUE, IERR) KERR = KERR.OR.IERR IPRNT = INT(VALUE(1)) C C RESTART SKIPS C ELSE IF (KEYWRD .EQ. 'SKIP') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, VALUE, IERR) KERR = KERR.OR.IERR MFILE = INT(VALUE(1)) + 1 C C RESTART CHECK C ELSE IF (KEYWRD .EQ. 'RSTR') THEN LRSTRT = .TRUE. C C CONTINUATION FLAG AND ASSOCIATED PARAMETERS C ELSE IF (KEYWRD .EQ. 'CNTN') THEN LCNTUE = .TRUE. C C NUMBER OF MESH POINTS IN THE REGRID C ELSE IF (KEYWRD .EQ. 'JJRG') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, VALUE, IERR) KERR = KERR.OR.IERR JJREGD = INT(VALUE(1)) LREGRD = .TRUE. C C PERCENTAGE OF REGRID POINTS DEDICATED TO ADAPTION C ELSE IF (KEYWRD .EQ. 'PCAD') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, PCTADP, IERR) KERR = KERR.OR.IERR C C RATIO OF GRADIENT REGRID POINTS TO CURVATURE POINTS C ELSE IF (KEYWRD .EQ. 'RGTC') THEN CALL CKXNUM (LINE, 1, LOUT, NVAL, RATGTC, IERR) KERR = KERR.OR.IERR C C LAST CARD C ELSE IF (KEYWRD(1:3) .EQ. 'END') THEN GO TO 6000 ELSE C C-------------THE KEYWORDS AFTER HERE ARE NOT OPERATIONAL------------ C C--------------END OF KEYWORDS-------------------- C C C TO GET HERE, AN INVALID KEYWORD WAS READ C WRITE (LOUT, *) ' ERROR...ILLEGAL KEYWORD' KERR = .TRUE. ENDIF GO TO 90 C C CHECK THE REACTANT AND PRODUCT SUMS C 6000 CONTINUE IF (KERR) STOP C SUMF = 0. SUMO = 0. SUMP = 0. DO 6100 K = 1, KK SUMF = SUMF+FUEL(K) SUMO = SUMO+OXID(K) SUMP = SUMP+PROD(K) 6100 CONTINUE C C NORMALIZE REACTANT AND PRODUCT FRACTIONS C DO 6200 K = 1, KK FUEL(K) = FUEL(K)/SUMF OXID(K) = OXID(K)/SUMO PROD(K) = PROD(K)/SUMP 6200 CONTINUE C IF (ABS(SUMF-1.0) .GT. 1.E-6) WRITE (LOUT, *) 1 ' CAUTION...FUEL FRACTIONS SUM TO ', SUMF IF (ABS(SUMO-1.0) .GT. 1.E-6) WRITE (LOUT, *) 1 ' CAUTION...OXIDIZER FRACTIONS SUM TO ', SUMO IF (ABS(SUMP-1.0) .GT. 1.E-6) WRITE (LOUT, *) 1 ' CAUTION...PRODUCT FRACTIONS SUM TO ', SUMP C C CHECK FOR NECESSARY INPUT C IF (.NOT. NEC(1)) THEN WRITE (LOUT, *) ' ERROR..."XEND" NOT SPECIFIED' STOP ENDIF C IF (.NOT. NEC(2) ) THEN WRITE (LOUT, *) ' ERROR..."TFUEL" NOT GIVEN' STOP ENDIF IF (.NOT. NEC(3) ) THEN WRITE (LOUT, *) ' ERROR..."TOXID" NOT GIVEN' STOP ENDIF C IF (.NOT. NEC(4) ) THEN WRITE (LOUT, *) ' ERROR..."VFUEL" NOT GIVEN' STOP ENDIF C IF (.NOT. NEC(5) ) THEN C WRITE (LOUT, *) ' ERROR..."VOXID" NOT GIVEN' C STOP C C C C ENDIF C C MAKE SURE THE INITIAL GRID POINTS ARE IN ORDER C IF ((.NOT. CNTNUD) .AND. (.NOT. LUMESH)) THEN NPTS = NP DO 5550 N = 2, NPTS IF (X(N-1) .GE. X(N)) THEN WRITE (LOUT, *) 1 ' ERROR...INITIAL GRID IS OUT OF ORDER' STOP ENDIF 5550 CONTINUE ENDIF C IF (NTEMP.GT.0) THEN C C MAKE SURE THE (X,T) PAIRS ARE IN ORDER C DO 5500 N = 2, NTEMP IF (XX(N-1) .GE. XX(N)) THEN WRITE (LOUT, *) 1 ' ERROR...SPECIFIED TEMPERATURES ARE OUT OF ORDER' STOP ENDIF 5500 CONTINUE C C MAKE SURE THE GIVEN TEMPERATURES SPAN THE XEND-XSTR DOMAIN C IF (.NOT.LRSTRT .OR. .NOT.CNTNUD .OR. LUSTGV) THEN IF (XX(1).GT.0.0 .OR. XX(NTEMP).LT.XEND) THEN WRITE (LOUT, *) 1 ' ERROR...GIVEN TEMPERATURE PROFILE DOES NOT SPAN XEND-XSTR' STOP ENDIF ENDIF ENDIF C C SET OPTIONAL INPUT IF NEEDED C IF (.NOT. NOPT(1)) WMIX = (XEND-0.0)*0.50 C IF (.NOT. NOPT(2)) XCEN = (XEND-0.0)*0.35 C C formats C RETURN 7000 FORMAT(A4, A) 8000 FORMAT(10X, A4, A76) END C C---------------------------------------------------------------------- C SUBROUTINE TEMP(NPTS, X, XX, TT, T) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C DIMENSION XX(*), TT(*) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C THIS SUBROUTINE USES BISECTION TO LINEARLY INTERPOLATE C AN ARRAY OF XX,TT PAIRS. GIVEN AN XX,TT PAIR THIS ROUTINE C RETURNS THE INTERPOLATED VALUE OF THE T AT THE POINT X. C C INPUT- C NPTS - NUMBER OF XX,TT PAIRS. C X - LOCATION AT WHICH INTERPOLATED T IS DESIRED. C XX - ARRAY OF X POINTS AT WHICH TT ARE GIVEN. C TT - ARRAY OF T VALUES AT THE XX LOCATIONS. C C OUTPUT- C T - INTERPOLATED T AT POINT X C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C check for x outside (1,npts) C IF (X .LE. XX(2)) THEN N = 2 S = (TT(N) - TT(N-1)) / (XX(N) - XX(N-1)) ELSEIF (X .GE. XX(NPTS-1)) THEN N = NPTS-1 S = (TT(N+1) - TT(N)) / (XX(N+1) - XX(N)) ELSE NLO = 1 NHI = NPTS S = 0.0 C C bisect interval C 50 CONTINUE N = (NLO+NHI)/2 IF (X .LT. XX(N)) THEN IF (X .LT. XX(N-1)) THEN NHI = N GO TO 50 ELSEIF (X .EQ. XX(N-1)) THEN N = N-1 ELSE S = (TT(N) - TT(N-1)) / (XX(N) - XX(N-1)) ENDIF ELSEIF (X .GT. XX(N)) THEN IF (X .GT. XX(N+1)) THEN NLO = N GO TO 50 ELSEIF (X .EQ. XX(N+1)) THEN N = N + 1 ELSE S = (TT(N+1) - TT(N)) / (XX(N+1) - XX(N)) ENDIF ENDIF ENDIF C 100 CONTINUE T = TT(N) + S * (X - XX(N)) RETURN END C C-------------------------------------------------------------------- C SUBROUTINE REASEN (II, KK, JJ, NATJ, LDA, LENRGY, LMULTI, 1 LTDIF, LRCRVR, LOUT, LVARMC, LTIME, P, 2 WT, YOXID, YFUEL, DT, NTEMP, 3 XGIVEN, TGIVEN, X, SN, S, WNDFAC, ABS, REL, 4 SCRCHK, YV, VISC, COND, D, DKJ, TDR, ICKWRK, 5 RCKWRK, IMCWRK, RMCWRK, F, FN, A, DS, SSAVE, 6 IP, RFAC) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C DIMENSION WT(*), YOXID(*), YFUEL(*), XGIVEN(*), TGIVEN(*), X(*), 1 SN(NATJ,*), S(NATJ, *), F(NATJ,*), FN(NATJ,*), 2 SCRCHK(KK,*), YV(KK,*), VISC(*), COND(*), D(KK,*), 3 TDR(KK,*), DKJ(KK,KK,*), DS(*), SSAVE(*), IP(*), 4 ICKWRK(*), RCKWRK(*), IMCWRK(*), RMCWRK(*), A(LDA, *) C COMMON /LOCS/ NT, NG, NF, NH,NTH ,NYS, NY COMMON/JI/NUMJIA,NKKK C LOGICAL LENRGY, LMULTI, LTDIF, LVARMC, LTIME C CHARACTER*18 ISENSI DATA ISENSI /'SENSITIVITY '/ C C*****precision > single C DATA ONE/1.0E0/ C*****END precision > single C*****precision > double DATA ONE/1.0D0/ C*****END precision > double LTIME = .FALSE. LVARMC = .TRUE. ML = 2*NATJ - 1 MU = 2*NATJ - 1 C C FORM THE JACOBIAN, AND AT THE SAME TIME EVALUATE THE C UNPERTURBED FUNCTION FN C CALL JACOB (KK, JJ, NATJ, LENRGY, LMULTI, LTDIF, LTIME, 1 P, WT, YOXID, YFUEL, DT, NTEMP, XGIVEN, TGIVEN, 2 X, SN, S, WNDFAC, ABS, REL, SCRCHK, YV, 3 TFUEL, TOXID, VFUEL, VOXID, AFUEL, AOXID, 3 VISC, COND, D, DKJ, 4 TDR, ICKWRK, RCKWRK, IMCWRK, RMCWRK, F, FN, A, DS, 5 SSAVE, RFAC) C C FACTOR THE JACOBIAN C C*****precision > double CALL DGBFA (A, LDA, NATJ*JJ, ML, MU, IP, INFO) C*****END precision > double C C*****precision > single C CALL SGBFA (A, LDA, NATJ*JJ, ML, MU, IP, INFO) C*****END precision > single C IF (INFO .NE. 0) THEN WRITE(LOUT,*) 'ERROR IN FACTORING THE JACOBIAN, INFO = ',INFO STOP ENDIF C C WRITE (LSAVE) ISENSI C C COMPUTE THE RAW SENSITIVITY COEFFICIENTS WITH RESPECT TO C THE RATE CONSTANTS, D(MASS FRACTION)/D(RATE CONSTANT). C WRITE(LRCRVR,*)S(NTH,1) DO 1000 I = 1, II C CALL CKRDEX (I, RCKWRK, SAVEP) DP = REL*SAVEP + ABS CALL CKRDEX (-I, RCKWRK, SAVEP+DP) C CALL FUN (KK, JJ, NATJ, LENRGY, LMULTI, LTDIF, LVARMC, LTIME, 1 P, WT, YOXID, YFUEL, DT, NTEMP, XGIVEN, TGIVEN, 2 X, SN, S, WNDFAC, SCRCHK(1,1), SCRCHK(1,5), YV, 3 SCRCHK(1,2), SCRCHK(1,3), SCRCHK(1,4), 4 TFUEL, TOXID, VFUEL, VOXID, AFUEL, AOXID, 4 VISC, COND, D, DKJ, 5 TDR, ICKWRK, RCKWRK, IMCWRK, RMCWRK, F, RFAC) C CALL CKRDEX (-I, RCKWRK, SAVEP) C DO 100 J = 1, JJ DO 100 N = 1, NATJ SN(N,J) = - (F(N,J)-FN(N,J)) / DP 100 CONTINUE C C C*****precision > double CALL DGBSL (A, LDA, NATJ*JJ, ML, MU, IP, SN, 0) C*****END precision > double C C*****precision > single C CALL SGBSL (A, LDA, NATJ*JJ, ML, MU, IP, SN, 0) C*****END precision > single C C C SN(N,J) NOW CONTAINS THE RAW SENSIVITY MATRIX DS(N,J)/DAi C C C NORMALIZE THE SENSIVITY COEFFICIENTS C DO 200 J = 1, JJ SN(NT,J) = SN(NT,J) * MAX(ONE,SAVEP) / S(NT,J) SN(NTH,J)= SN(NTH,J)*SAVEP/S(NTH,J) 200 CONTINUE C DO 500 J = 1, JJ SUM = 0.0 DO 300 L = 1, KK SUM = SUM + SN(NYS+L,J) / WT(L) 300 CONTINUE CALL CKMMWY (S(NY,J), ICKWRK, RCKWRK, WTM) WTMSUM = WTM * SUM DO 400 K = 1, KK SN(NYS+K,J) = SAVEP * (SN(NYS+K,J) / S(NYS+K,J) 2 - WTMSUM) 400 CONTINUE WRITE(LRCRVR,7020)X(J),SN(NYS+NUMJIA,J),SN(NTH,J) 500 CONTINUE C c SN(NTH,1) = SN(NTH,1) * SAVEP / S(NTH,1) C WRITE (LSAVE) I, ((SN(N,J),N=1,NATJ), J=1,JJ) C WRITE(LSAVE,7020) C 7020 FORMAT(1X,'X=',2X,F6.3,3X,E16.3,3X,E16.3) 1000 CONTINUE C WRITE(LRCRVR,*)'complet 123' RETURN END C C C----------------------------------------------------------------------- SUBROUTINE BSEN (II, KK, JJ, NATJ, LDA, LENRGY, LMULTI, 1 LTDIF, LRCRVR, LOUT, LVARMC, LTIME, P, 2 WT, YOXID, YFUEL, DT, NTEMP, 3 XGIVEN, TGIVEN, X, SN, S, WNDFAC, ABS, REL, 4 SCRCHK, YV, VISC, COND, D, DKJ, TDR, ICKWRK, 5 RCKWRK, IMCWRK, RMCWRK, F, FN, A, DS, SSAVE, 6 IP, RFAC) C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C DIMENSION WT(*), YOXID(*), YFUEL(*), XGIVEN(*), TGIVEN(*), X(*), 1 SN(NATJ,*), S(NATJ, *), F(NATJ,*), FN(NATJ,*), 2 SCRCHK(KK,*), YV(KK,*), VISC(*), COND(*), D(KK,*), 3 TDR(KK,*), DKJ(KK,KK,*), DS(*), SSAVE(*), IP(*), 4 ICKWRK(*), RCKWRK(*), IMCWRK(*), RMCWRK(*), A(LDA, *) C COMMON /LOCS/ NT, NG, NF, NH,NTH ,NYS, NY COMMON/JI/NUMJIA,NKKK C LOGICAL LENRGY, LMULTI, LTDIF, LVARMC, LTIME,LLJIA C C C*****precision > single C DATA ONE/1.0E0/ C*****END precision > single C*****precision > double DATA ONE/1.0D0/ C*****END precision > double LLJIA=.TRUE. LTIME = .FALSE. LVARMC = .TRUE. ML = 2*NATJ - 1 MU = 2*NATJ - 1 C C FORM THE JACOBIAN, AND AT THE SAME TIME EVALUATE THE C UNPERTURBED FUNCTION FN C CALL JACOB (KK, JJ, NATJ, LENRGY, LMULTI, LTDIF, LTIME, 1 P, WT, YOXID, YFUEL, DT, NTEMP, XGIVEN, TGIVEN, 2 X, SN, S, WNDFAC, ABS, REL, SCRCHK, YV, 3 TFUEL, TOXID, VFUEL, VOXID, AFUEL, AOXID, 3 VISC, COND, D, DKJ, 4 TDR, ICKWRK, RCKWRK, IMCWRK, RMCWRK, F, FN, A, DS, 5 SSAVE, RFAC) C C FACTOR THE JACOBIAN C C*****precision > double CALL DGBFA (A, LDA, NATJ*JJ, ML, MU, IP, INFO) C*****END precision > double C C*****precision > single C CALL SGBFA (A, LDA, NATJ*JJ, ML, MU, IP, INFO) C*****END precision > single C IF (INFO .NE. 0) THEN WRITE(LOUT,*) 'ERROR IN FACTORING THE JACOBIAN, INFO = ',INFO STOP ENDIF C DO 1001 K=1,KK CALL JIA(K,RCKWRK,RD,31,LLJIA,REL,ABS) CALL FUN (KK, JJ, NATJ, LENRGY, LMULTI, LTDIF, LVARMC, LTIME, 1 P, WT, YOXID, YFUEL, DT, NTEMP, XGIVEN, TGIVEN, 2 X, SN, S, WNDFAC, SCRCHK(1,1), SCRCHK(1,5), YV, 3 SCRCHK(1,2), SCRCHK(1,3), SCRCHK(1,4), 4 TFUEL, TOXID, VFUEL, VOXID, AFUEL, AOXID, 4 VISC, COND, D, DKJ, 5 TDR, ICKWRK, RCKWRK, IMCWRK, RMCWRK, F, RFAC) C C CALL JIA(-K,RCKWRK,RD,31,LLJIA,REL,ABS) C----------------------------------------------------------------------- DO 100 J = 1, JJ DO 100 N = 1, NATJ SN(N,J) = - (F(N,J)-FN(N,J)) / RD 100 CONTINUE C C C*****precision > double CALL DGBSL (A, LDA, NATJ*JJ, ML, MU, IP, SN, 0) c----------------------------------------------------------------------- SN(NTH,1)= SN(NTH,1)*RD/S(NTH,1) WRITE(LRCRVR,7025)K,SN(NTH,1) 7025 FORMAT(1X,'K=',2X,I3,3X,E16.3) 1001 CONTINUE CALL FUN (KK, JJ, NATJ, LENRGY, LMULTI, LTDIF, LVARMC, LTIME, 1 P, WT, YOXID, YFUEL, DT, NTEMP, XGIVEN, TGIVEN, 2 X, SN, S, WNDFAC, SCRCHK(1,1), SCRCHK(1,5), YV, 3 SCRCHK(1,2), SCRCHK(1,3), SCRCHK(1,4), 4 TFUEL, TOXID, VFUEL, VOXID, AFUEL, AOXID, 4 VISC, COND, D, DKJ, 5 TDR, ICKWRK, RCKWRK, IMCWRK, RMCWRK, F, RFAC) WRITE(LRCRVR,*)'complet 123' RETURN END