aspen RPLUG 的fortran动力学子程序纠错？

这是我仿照论坛上的示例修改的子程序，请高手帮忙看一下，从语法上看看有没什么错误，多谢！编译后无法运行。示例文件在这里 https://bbs.mahoupao.net/forum.php?mod=viewthread&tid=76210 C User Kinetics Subroutine for Urea synthensis, 2NH3+CO2=CARB; CARB=UREA+H2O C SUBROUTINE USURA (SOUT, NSUBS, IDXSUB, ITYPE, NINT, 2 INT, NREAL, REAL, IDS, NPO, 3 NBOPST, NIWORK, IWORK, NWORK, WORK, 4 NC, NR, STOIC, RATES, FLUXM, 5 FLUXS, XCURR, NTCAT, RATCAT, NTSSAT, 6 RATSSA, KCALL, KFAIL, KFLASH, NCOMP, 7 IDX, Y, X, X1, X2, 8 NRALL, RATALL, NUSERV, USERV, NINTR, 9 INTR, NREALR, REALR, NIWR, IWR, 1 NWR, WR) C IMPLICIT NONE C C DECLARE VARIABLES USED IN DIMENSIONING C INTEGER NSUBS, NINT, NPO, NIWORK,NWORK, + NC, NR, NTCAT, NTSSAT,NCOMP, + NRALL, NUSERV,NINTR, NREALR,NIWR, + NWR C #include "pputl_ppglob.cmn" #include "ppexec_user.cmn" EQUIVALENCE (RMISS, USER_RUMISS) EQUIVALENCE (IMISS, USER_IUMISS) C C C C C.....RCSTR... #include "rcst_rcstri.cmn" #include "rxn_rcstrr.cmn" C C.....RPLUG... #include "rplg_rplugi.cmn" #include "rplg_rplugr.cmn" EQUIVALENCE (XLEN, RPLUGR_UXLONG) EQUIVALENCE (DIAM, RPLUGR_UDIAM) C C.....RBATCH... #include "rbtc_rbati.cmn" #include "rbtc_rbatr.cmn" C C.....PRES-RELIEF... #include "rbtc_presrr.cmn" C C.....REACTOR (OR PRES-RELIEF VESSEL) PROPERTIES... #include "rxn_rprops.cmn" EQUIVALENCE (TEMP, RPROPS_UTEMP) EQUIVALENCE (PRES, RPROPS_UPRES) EQUIVALENCE (VFRAC, RPROPS_UVFRAC) EQUIVALENCE (BETA, RPROPS_UBETA) EQUIVALENCE (VVAP, RPROPS_UVVAP) EQUIVALENCE (VLIQ, RPROPS_UVLIQ) EQUIVALENCE (VLIQS, RPROPS_UVLIQS) C #include "shs_stwork.cmn" EQUIVALENCE (MKBAS, STWORK_NDUM) EQUIVALENCE (MKPHAS, STWORK_NBLM) EQUIVALENCE (MTAPP, STWORK_NCOVAR) EQUIVALENCE (MKBASS, STWORK_NWR) EQUIVALENCE (MTAPPS, STWORK_NIWR) EQUIVALENCE (SSALT, STWORK_RDUM1) EQUIVALENCE (VSALT, STWORK_RDUM2) EQUIVALENCE (FSALT, STWORK_FFSALT) C #include "dms_ncomp.cmn" #include "dms_plex.cmn" EQUIVALENCE (IB(1), B(1)) C C DECLARE ARGUMENTS C INTEGER IDXSUB(NSUBS),ITYPE(NSUBS), INT(NINT), + IDS(2),NBOPST(6,NPO),IWORK(NIWORK), + IDX(NCOMP), INTR(NINTR), IWR(NIWR), + NREAL, KCALL, KFAIL, KFLASH,I REAL*8 SOUT(1), WORK(NWORK), + STOIC(NC,NSUBS,NR), RATES(5), + FLUXM(1), FLUXS(1), RATCAT(NTCAT), + RATSSA(NTSSAT), Y(NCOMP), + X(NCOMP), X1(NCOMP), X2(NCOMP) REAL*8 RATALL(NRALL),USERV(NUSERV), + REALR(NREALR),WR(NWR), XCURR C C C DECLARE LOCAL VARIABLES C INTEGER IMISS, MKBAS, MKPHAS,MTAPP, MKBASS, + MTAPPS,LMW REAL*8 REAL(NREAL), XL(7), B(1), RMISS, XLEN, + DIAM, TEMP, PRES, VFRAC, BETA, + VVAP, VLIQ, VLIQS, SSALT, VSALT, + FSALT, RTEMP, DENLIQ,KKA, KKB, + AREA, RTA, RTB, SUMR INTEGER XMW #include "dms_ipoff1.cmn" C C INITIALIZE RATES C C C STATEMENT FUNCTIONS FOLLOW C C C BEGIN EXECUTABLE CODE C XMW(I) = LMW + I C C SET PLEX OFFSETS C LMW = IPOFF1_IPOFF1(306) DO 100 I = 1, NC RATES(I) = 0D0 XL(I) = 0D0 100 CONTINUE C DO 101 I = 1, NCOMP XL(IDX(I)) = X(I) 101 CONTINUE C C The structure in the array SOUT is as follows: C C SOUT(1) - SOUT(NCC) : Component flowrates(kg-moles/sec) C SOUT(NCC+1) : Total flowrates(kg-moles/sec) C SOUT(NCC+2) : Temperature(K) C SOUT(NCC+3) : Pressure(N/SQM) C SOUT(NCC+4) : Mass enthalpy(J/KG) C SOUT(NCC+5) : Molar vapor fraction C SOUT(NCC+6) : Molar liquid fraction C SOUT(NCC+7) : Mass entropy(J/KG-K) C SOUT(NCC+8) : Mass density(KG/CUM) C SOUT(NCC+9) : Molecular Weight C C Details of this infomation can be found in the User Guide C on page 17-10 for 8.5-6 or Reference Manual Vol. 6, User Models C for Release 9. C C Set Reactor Temperature C RTEMP = SOUT(NCOMP_NCC+2) C C Get liquid mixture density C C DENLIQ = 1.0 / STWORK_VL C C Compute Concentration for reactants: C C CONC(i) = X(i) * RHO C C CONAA = XL(1) * DENLIQ C CONACT = XL(2) * DENLIQ C C C C CALCULATE RATES OF REACTIONS C C The rates must be set for every single component exists in the C system(both reactants and products) according to stoichiometry. C C RATES - The rates for all components C C Acetone + Allyl Alcohol -> N-Propyl Propionate C XL1=UREA,XL2=CARB,XL3=NH3,XL4=CO2,XL5=H2O C FTERM = REALR(1) * EXP( -REALR(2) / PPGLOB_RGAS / RTEMP ) C STERM = conaa * conact**0.5 C RATE = FTERM * STERM * RCSTRR_VOLRC C RATES(1) = -rate C RATES(2) = -rate C RATES(3) = rate C AREA=DFLOAT(RPLUGI_NTUBE)*3.141592654/4.*DIAM**2 KKA=XL(2)/XL(4)/XL(4)/XL(3) KKB=XL(1)*XL(5)/XL(2) RTA=1.0e12*EXP(-1E07/PPGLOB_RGAS/RTEMP)*(XL(4)**2*XL(3)-XL(2)/KKA) RTB=1.5e09*EXP(-1E08/PPGLOB_RGAS/RTEMP)*(XL(2)-XL(1)*XL(5)/KKB) RATES(1)=AREA*RTB RATES(2)=AREA*(RTA-RTB) RATES(3)=AREA*RTA RATES(4)=-2*RTA*AREA RATES(5)=RTB*AREA C Checking Mass Balance for RATES Vector C SUMR = 1.0 DO 102 I=1,NCOMP_NCC SUMR = SUMR + RATES(I)*B(XMW(I)) 102 CONTINUE IF (DABS(1.0 - SUMR) .GT. 1E-6) THEN WRITE(USER_NHSTRY,*) + 'Calculated Rates Not in Mass Balance!' WRITE(USER_NHSTRY,*) + 'Errors in User Kinetics Subroutines!' ENDIF write(USER_NHSTRY,*) 'volrc = ', RCSTRR_VOLRC write(USER_NHSTRY,*) 'vfrrc = ', RCSTRR_VFRRC RETURN #undef P_NPOFF1 END 显示全部