Tew at al Source code of: A Stirling engine Computer model, written in Fortran by Tew at al
The engine of the future project
The following is a ASCII version of the FORTRAN 5 program, run on a UNIVAC 1110.
I found it in A_Stirling_Engine_Computer_Model_for_Per.pdf
DOE/NASA/1011-78/24
NASA TM-78884
http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19780022051.pdf
and typed it in a ASCII file that can be used as a FORTRAN program.
I published a version in Quick Basic on http://jordaan.info/greengasoline/stirlingtewbas.html
This version is the original fortran program.
How do you download the program? In fact it is already on your computer.
You can select the part between the ===== lines and copy it in a simple editor and save it with a name with the extention .f
Download a FORTRAN compiler and compile the program.
Warning: There will still be errors in this program.
Please e-mail me when you find one. (hvd@haghen.nl)

Now I'm working on a better version.
This page is ment for progammers.
For the html-javascript version of the Isothermal second order calculation Program (Written by Wiliam R. Martini in Fortran) Click here
The program is in the source code.
A more advanced version of this program will be published in a Quick Basic version on this site soon.
Click here for HOME
'====================================================================
C PAGE 47 (53 IN MANUAL PDF)----------------
C MAIN PROGRAM
C DOES INITIALIZATION. UPDATES PISTON POSITIONS--EXPANSION, COMPRESSION
C AND BUFFER SPACE VOLUMES. INTEGRATES P*DELTA(V) TO DETERMINE WORK.
C DETERMINES WHEN CYCLES IS COMPLETE.
      DIMENSION X(2),V(13),REYNO(13),        QR(5),TMAX(12)
      DIMENSION TG(12),TGA(13),TM(13),TMAVG(13),DTM(13)
      COMMON DTIME,F(12),WG(13),DELP(13),P,PE,PC
      COMMON /CYC/ QEXP,QHEATR,QREGEN,QCOOLR,QCOMP,ENFRTH,ENFHTR,
     1ENFCTR,ENFRTC,WRKP,WRKEXP,
     2WRKCMP,WRKLEX,WRKLCM,WRKLP,WRKH,SUMNUM(13),SUMDEN(13),ENTH,
     3WDISP,WPIST,WTPIST,QEXPN,QEXPP,QHEATN,QHEATP,
     4QCOOLN,QCOOLP,QCOMPN,QCOMPP,QREGN,QREGP,WDISPP,WDISPN,WPISTP,
     5WPISTN,QCNDRI,QCNDRO,QCNDCL,QCNDD,QCNDCN, QSHTL,TGEXPA,TGCMPA
      COMMON /DELTAP/AD,AP,		HA(13)
      COMMON  /TAVCYC/TGCYC(12),TGACYC(13),TMCYC(13)
      COMMON/CTIT/NCURV,ISET1,ISET2,ISET3
      COMMON /JNDEX/JIP,NOCYC
      NAMELIST /STRLNG/P,REALGS,FIPCYC,IPRINT,ITMPS,COEFF,T3,P3,OMEGA,
     *FACT1,FACT2,NOCYC,NSTRT,NOEND,MWGAS,TG,TGA,TM,RHCFAC,JIP
      DATA		AD,AP,AR,APISTR /
     1			6.0,6.0,0.11,0.60/
      DATA RCRANK,E,RODL,V3CL,DIAMD,DIAMP,DIAMDR,DIAMPR/
     10.543,0.82,1.81,24.84,2.75,2.75,0.375,0.875/
      DATA   V/0.0,2*1.5948,.933,5*0.624,3*0.240,0.0/
   10 CONTINUE
C READ IN INITIAL PRESSURES, TEMPERATURES AND OTHER PARAMETERS WHICH
C DEFINE THE NATURE OF THE RUN TO BE MADE
      READ (5,STRLNG,END=1000)
      DTIME=1./(OMEGA*FIPCYC)
      WRITE (6,STRLNG)
C     INITIALIZATION
      FMULT=1.
      ITER=0
      TIME=0.0
      PSI=0.0
      PSIDEG=0.0
      SAVET=0.0
      IPLOT=0
      PRSUM=0.0
      NOITPC=0
      DO 11 I=1,13
      IF(I.EQ.13) GO TO 12
      IGCYC(I)=TG(I)
   12 TGACYC(I)=TGA(I)
   11 TMCYC(I)=TM(I)
      CPCYC=1./(OMEGA*DTIME)
C R,CP,CV IN UNITS OF IN-LBF/ LBM-DEG. R
      GAMMA=1.394
      R=9197.
      CP=32557.
      CV=23360.
      B=.02358
      IF (MWGAS.NE.4) GO TO 15
      B=.01613
C PAGE 48 (54 IN MANUAL PDF)-----------------------
      GAMMA=1.667
      R=4632.
      CP=11580.
      CV=6948.
   15 CONTINUE
      IP=IPRINT
      IPRNT2=IPRINT/25
      ADR=AD-AR
      APR=AP-APISTR
      X1MIN=SQRT((RODL-RCRANK)**2-E**2)
      RODYMX=SQRT(RODL**2-(E-RCRANK)**2)
      RODYMN=SQRT(RODL**2-(E+RCRANK)**2)
      STROKV=2.*(RODYMX-RODYMN)
      RODY=RODYMX
      X(1)=RODYMX
      X(2)=-RODYMX
      V(1)=AD*(X(1)-X1MIN)+0.965
      V(13)=2.*ADR*(RODYMX-RODY)+0.353
      V3=APR*(-X(2)-X1MIN)+V3CL
      V3STAR=V3
      VOVRT=0.0
      DO 30 I=1,13
   30 VOVRT=VOVRT+V(I)/(TGA(I)+B*P*REALGS)
      W3=P3*V3STAR/(R*(T3+B*P3*REALGS))
      ALPHA=W3*R/VOVRT
      SIG3=ALPHA/(P+ALPHA)
      W=W3/SIG3
      PE=P
      PC=P
      DO 50 I=1,13
   50 WG(I)=P*V(I)/(R*(TGA(I)+B*P*REALGS))
      ISTART=1
      J=1
  100 CONTINUE
      IF (J .LT. ISTART) GO TO 140
  125 IF (JIP.GT.0) GO TO 140
      IF (J.LT.(NOCYC-1)) WRITE (6,130)
  130 FORMAT ('2  TIME   ANGLE    X(1)   X(2)   T4',7X,'T9',7X,'T1',4X,'T2
     1', 5X,' PE      PC     PRIN    PROUT ',3X,'F(4)',5X,'F(9) ')
      IF (J.EQ.(NOCYC-1)) WRITE (6,132)
  132 FORMAT ('  TIME    ANGLE    RE1      RE2     RE3     RE4      RE5'
     1'     RE6     RE7     RE8     RE9      RE10     RE11     RE12    '
     2'RE13 ')
C 'FORGOTTEN IN LINE AFTER RE12?
C     1'     RE6     RE7     RE8     RE9      RE10     RE11     RE12     
      IF (J.EQ.NOCYC) WRITE (6,133)
  133 FORMAT ('  TIME    ANGLE     F1       F2       F3     F4       F5'
     1'      F6      F7      F8      F9       F10       F11      F12')
C     BEGINNING OF LOOP
  140 ITER=ITER+1
      IF(J.GE.(NOCYC-4)) IPRINT=IPRNT2
      IF  ((TIME+SAVET).GT.5.0) GO TO 700
      IP=IP+1
      VTOTL=0.0
      DO 142 IJK=1,13
  142 VTOTL=VTOTL+V(IJK)
      IF (JIP.GT.0) GO TO 170
      IF (IP.LT.IPRINT) GO TO 180
      IF (J .LT. ISTART) GO TO 170
      IF (J.LT.(NOCYC-1)) WRITE (6,160) TIME,PSIDEG,X,TG(4),TG(9),
     1TGA(1),TGA(13),PE,PC,PRIN,PROUT,F(4),F(9)
  160 FORMAT (1X,F7.4,F8.1,2F6.2,1X,2F8.2,F8.0,F6.0,1X,4F8.2,2F8.5)
      IF (J.EQ.(NOCYC-1)) WRITE (6,162) TIME,PSIDEG,REYNO
  162 FORMAT (1X,F7.4,F8.1,13F8.0)
      IF (J.EQ.NOCYC) WRITE (6,163) TIME,PSIDEG,F
  163 FORMAT (1X,F7.4,F8.1,12F8.5)
      IF (ITMPS.EQ.1) WRITE (6,660) TG,TGA,TM
  660 FORMAT ('  TG=',4X,12F8.0/,' TGA=',13F8.0/,'  TM=',13F8.0/)
C PAGE 49 (55 IN MANUAL PDF)--------------------------------
  170 IP=0
  180 TIME=DTIME*FLOAT(ITER)
      IF(ITER.EQ.IFIX(.3/DTIME) .AND. J.LT.6) GO TO 700
C CALL HEAT EXCHANGER SUBROUTINE TO UPDATE GAS LUMP MASSES, FLOW RATES,
C TEMPERATURES, PRESSURE + REGENERATOR METAL TEMPERATURES.
      CALL HEATX (TIME,TG,TGA,TM,REYNO,QR,X,RODY,J,DIAMD,X1MIN,
     1RODYMX,V,SIG3,W,MWGAS,RHCFAC,REALGS,CONDID,NOCYC)
C CONDID is not used ??????????
      PRSUM=PRSUM+P
      NOITPC=NOITPC+1
C UPDATE EXPANSION AND COMPRESSION SPACE PRESSURES. (LBF/IN2)
      PEOLD=PE
      PCOLD=PC
      POLD=P
      IF (J.LT.(NOCYC-5)) GO TO 240
      PRIN=(DELP(5)+DELP(6)+0.5*DELP(7))*FMULT+P
      PROUT=P-(0.5*DELP(7)+DELP(8)+DELP(9))*FMULT
      PE=(DELP(2)+DELP(3)+DELP(4))*FMULT+PRIN
      PC=PROUT-(DELP(10)+DELP(11)+DELP(12))*FMULT
      GO TO 245
  240 PE=P
      PC=P
      PRIN=P
      PROUT=P
  245 CONTINUE
      AE=(PE+PEOLD)/2.
      AC=(PC+PCOLD)/2.
      A=(P+POLD)/2.
      DELPRG=PRIN-PROUT
      DELPHT=PE-PRIN
      DELPCL=PROUT-PC
C CALCULATE INTERNAL ENERGY OF WORKING SPACE GAS. (FT-LBF)
      UTOTAL=0.0
      DO 250 I=1,13
  250 UTOTAL=WG(I)*CV*TGA(I)/12.+UTOTAL
  260 CONTINUE
C BUFFER SPACE BLEED FLOW,FLOTOB (LBM/SEC)
      IF ((PC-P3) .GE. 0.) GO TO 270
      SIG3P=-COEFF*(P3-PC)**.5/W
      GO TO 280
  270 SIG3P=COEFF*(PC-P3)**.5/W
  280 CONTINUE
      FLOTOB=SIG3P*W
      IF (FLOTOB.GE.0) ENTHWB=CP*TGA(13)*FLOTOB*DTIME/12.
      IF (FLOTOB.LT.0) ENTHWB=CP*T3*FLOTOB*DTIME/12.
      ENTH=ENTH+ENTHWB
      SIG3=SIG3+SIG3P*DTIME
      WG(13)=WG(13)-FLOTOB*DTIME
      IF (FLOTOB.LT.0.0) TGA(13)=TGA(13)-FLOTOB*DTIME*
     1(T3-TGA(13))/WG(13)
  520 P3OLD=P3
      P3=T3/(V3/(SIG3*W*R)-B*REALGS)
      A3=(P3OLD+P3)/2.
C UPDATE CRANK ANGLE, PSI(RADIANS), PSIDEG(DEGREES)
      PSIOLD=PSI
      PSI=PSI+OMEGA*DTIME*6.28318
      IF (PSI.GE.6.28318) PSI=PSI-6.28318
      PSIDEG=PSI*360./6.28318
C UPDATE PISTON POSITIONS (IN)
      RODY=SQRT(RODL**2-(E-RCRANK*COS(PSI))**2)
      X1OLD=X(1)
      X2OLD=X(2)
      X(1)=RCRANK*SIN(PSI)+RODY
      X(2)=RCRANK*SIN(PSI)-RODY
C UPDATE EXPANSION, COMPRESSION, AND BUFFER SPACE VOLUMES. (IN3)
      V1OLD=V(1)
      V13OLD=V(13)
C PAGE 50 (56 IN MANUAL PDF)-----------------------------------      
      V(1)=AD*(X(1)-X1MIN)+0.965
      V(13)=2.*ADR*(RODYMX-RODY)+0.353
      V3=APR*(-X(2)-X1MIN)+V3CL
C CALCULATE WORK (FT-LBF)--EXPANSION SPACE WORK,WRKEXP; COMPRESSION
C SPACE WORK,WRKCMP; TOTAL WORK; DISPLACER WORK,WDISP; POWER PISTON
C WORK,WPIST; PRESSURE DROP WORK LOSS IN EXPANSION AND COMPRESSION
C SPACE, WRKLEX AND WRKLCM
      WRKP=WRKP+AC*(V(13)-V13OLD)/12.
      DWEXP=AE*(V(1)-V1OLD)/12.
      WRKEXP=WRKEXP+DWEXP
      IF (DWEXP.GT.0.0) WEXPP=WEXPP+DWEXP
      IF (DWEXP.LT.0.0) WEXPN=WEXPN+DWEXP
      DWCMP=AC*(V(13)-V13OLD)/12.
      WRKCMP=WRKCMP+DWCMP
      IF (DWCMP.GT.0.0) WCMPP=WCMPP+DWCMP
      IF (DWCMP.LT.0.0) WCMPN=WCMPN+DWCMP
      WRKLEX=WRKLEX+(A-AE)*(V(1)-V1OLD)/12.
      WRKLCM=WRKLCM+(A-AC)*(V(13)-V13OLD)/12.
      WRKLP=WRKLP+(A-AC)*(V(13)-V13OLD)/12.
      WRKH=WRKH+(AE*(V(1)-V1OLD)+AC*(V(13)-V13OLD))/12.
      DWDISP=(AE*AD-AC*ADR)*(X(1)-X1OLD)/12.
      WDISP=WDISP+DWDISP
      IF (DWDISP.GT.0.0) WDISPP=WDISPP+DWDISP
      IF (DWDISP.LT.0.0) WDISPN=WDISPN+DWDISP
      DWPIST=AC*ADR*(X(2)-X2OLD)/12.
      WPIST=WPIST+DWPIST
      IF (DWPIST.GT.0.0) WPISTP=WPISTP+DWPIST
      IF (DWPIST.LT.0.0) WPISTN=WPISTN+DWPIST
      WTPIST=WTPIST+AC*ADR*(X(2)-X2OLD)/12.
      DO 530 I=1,13
      IF(I.EQ.13) GO TO 531
      TGCYC(I)=TGCYC(I)+TG(I)
  531 TGACYC(I)=TGACYC(I)+TGA(I)
  530 TMCYC(I)=TMCYC(I)+TM(I)
      IF (PSIOLD.GT.3.14159 .AND. PSI.LT.3.14159) GO TO 650
      GO TO 140
C     END OF LOOP
  650 CONTINUE
C CALCULATE AVERAGE GAS TEMPERATURES OVER THE CYCLE FOR EACH CONTROL VOLUME
      DO 533 I=1,13
      IF(I.EQ.13) GO TO 534
      TGCYC(I)=TGCYC(I)/CPCYC
  534 TGACYC(I)=TGACYC(I)/CPCYC
  533 TMCYC(I)=TMCYC(I)/CPCYC
      TGEXPA=TGACYC(1)
      TGCMPA=TGACYC(13)
      DO 667 II=1,6
      JJ=II+6
      TMAX(II)=0.0
  667 TMAX(JJ)=5000.
      IF (JIP.GT.0.AND.J.NE.(NOCYC-1)) GO TO 668
      IF (J.LT.(NOCYC-1)) WRITE (6,160) TIME,PSIDEG,X,P,P3,TGA(1),
     1TGA(13),PE,PC,PRIN,PROUT,V(1),V(13)
      IF (JIP.EQ.0) WRITE (6,660) TG,TGA,TM
  668 CONTINUE
      IF(J.LT.NSTRT .OR. J.GT.NOEND) GO TO 675
C CONVERGENCE SCHEME FOR REGENERATOR METAL TEMPERATURES.
      DO 670 JI=5.9
  670 DTM(JI)=FACT1*DTM(JI)+FACT2*SUMNUM(JI)/SUMDEN(JI)
      TM(5)=TM(5)-(5.*DTM(5)+4.*DTM(6)+3.*DTM(7)+2.*DTM(8)+   DTM(9))/3.
      TM(6)=TM(6)-(4.*DTM(5)+8.*DTM(6)+6.*DTM(7)+4.*DTM(8)+2.*DTM(9))/3.
      TM(7)=TM(7)-(1.*DTM(5)+2.*DTM(6)+3.*DTM(7)+2.*DTM(8)+   DTM(9))
      TM(8)=TM(8)-(2.*DTM(5)+4.*DTM(6)+6.*DTM(7)+8.*DTM(8)+4.*DTM(9))/3.
      TM(9)=TM(9)-(   DTM(5)+2.*DTM(6)+3.*DTM(7)+4.*DTM(8)+5.*DTM(9))/3.
      GO TO 675
  675 CONTINUE
 C PAGE 51 (57 IN MANUAL PDF)-------------
      DO 673 I=1,5
  673 QR(I)=0.0
      AVGWSP=PRSUM/FLOAT(NOITPC)
      CALL CYCL (TIME,UTOTAL,STROKV,AVGWSP)
      DO 532 I=1,13
      IF (I.EQ.13) GO TO 535
      TGCYC(I)=0.0
  535 TGACYC(I)=0.0
  532 TMCYC(I)=0.0
      PRSUM=0.0
      NOITPC=0
      J=J+1
      IF (J.NE.(NOCYC-1)) GO TO 680
      IPLOT=0
      ITER=0
      SAVET=TIME
      TIME=0.
  680 IF(J.LE.NOCYC) GO TO 125
  700 CONTINUE
      GO TO 10
 1000 CONTINUE
      STOP
      END
      
C SUBROUTINE HEATEX ---------
C THIS SUBROUTINE CALCALATES HEAT TRANSFER AND UPDATES TEMPERATURES AND
C PRESSURE
      SUBROUTINE HEATX (TIME,TG,TGA,TMA,REYNO,QR,X,RODY,JCYCLE,DIAMD,
     1X1MIN,RODYMX,V,SIG3,W,MWGAS,RHCFAC,REALGS,CONDID,NOCYC)
C CONDID is not used ??????????
      COMMON DTIME,F(12),WG(13),DELP(13),P,PE,PC
      COMMON /CYC/ QECP,QHEATR,QREGEN,QCOOLR,QCOMP,ENFRTH,ENFHTR,
     1ENFCTR,ENFRTC,WRKP,WRKEXP,
     2WRKCMP,WRKLEX,WRKLCM,WRKLP,WRKH,SUMNUM(13),SUMDEN(13),ENTH,
     3WDISP,WPIST,WTPIST,QEXPN,QEXPP,QHEATN,QHEATP,
     4QCOOLN,QCOOLP,QCOMPN,QCOMPP,QREGN,QREGP,WDISPP,WDISPN,WPISTP,
     5WPISTN,QCNDRI,QCNDRO,QCNDCL,QCNDD,QCNDCN, QSHTL,TGEXPA,TGCMPA
      COMMON /DELTAP/AD,AP,        HA(13)
      DIMENSION COEFX(40),DPX(40)
      DIMENSION TGA(13),       TMA(13),TGAO(13),TMAO(13), X(2)
      DIMENSION CPM(13),       Q(13),TMIX(13),TG(12),DH(13),ACS(13),
     1AHT(13),COND(13),VIS(13),REYNO(13),                    QR(5)
      DIMENSION V(13),WGOLD(13),XL(13),FAVG(13),DNSTY(13),FRICT(13)
      DATA DH/1.437,3*0.1190,5*0.003800,3*0.0409,1.487/,
     1    ACS/3.250,3*0.05561,5*0.4380,3*0.05124,3.150/,
     1    AHT/0.0,2*5.720,3.922,5*82.56,3*2.288,0.0/,
     1     CPM/4*10.0,5*.011,               4*10.0/,
     1     XL/0.0,2*3.06,2.098,5*0.178,3*0.587,0.0/,
     1 IDEX/1/
C JH(I)=HYDRAULIC DIAMETER,IN.; ACS(I)=FLOW AREA, IN2; AHT=HEAT TRANSFER,IN2;
C CPM(I)=METAL HEAT CAPACITY, BTU/LBM-DEG,R; XL(I)=LENGTHS FOR CALCULATING
C PRESSURE DROP,IN
C CP,CV ARE IN BTU/LBM-DEG.R;GAS CONSTANT, R, IS
C IN IN-LBF/(LBM-DEG R)
      B=.02358
      CP=3.484
      GAMMA=1.394
C PAGE 52 (58 IN MANUAL PDF)---------------------
      R=9197.
      IF (MWGAS.NE.4) GO TO 5
      B=.01613
      CP=1.239
      R=4632.
      GAMMA=1.667
    5 CONTINUE
      CV=CP/GAMMA
      DO 10 I=1,13
      IF (MWGAS.EQ.4) GO TO 6
C CONDUCTIVITY, COND(I) HAS UNITS BTU/IN-SEC-DEG.R
C VISCOSITY,VIS(I) HAS UNITS LBM/IN-SEC
      COND(I)=2.481E-9*TGA(I)+1.263E-6
      VIS(I)=4.683E-10*TGA(I)+2.619E-7
      GO TO 8
    6 COND(I)=1.962E-9*TGA(I)+1.115E-6
      VIS(I)=1.058E-9*TGA(I)+5.919E-7
    8 TMAO(I)=TMA(I)
   10 TGAO(I)=TGA(I)
      POLD=P
C CALCULATE PRESSURE,LBF/IN2
      SUMWT=0.0
      SUMV=0.0
      DO 14 I=1,13
      SUMWT=SUMWT+WG(I)*TGA(I)
   14 SUMV=SUMV+V(I)
      P=R*SUMWT/(SUMV-B*R*W*(1.-SIG3)*REALGS)
      FGAMMA=(P/POLD)**((GAMMA-1.)/GAMMA)
      DO 20 J=1,13
      TGAO(J)=TGAO(J)*FGAMMA
   20 TMIX(J)=TGAO(J)
C CALCULATE MASSES, WG(I), IN EACH GAS LUMP. ,LBM
      SUM=0.0
      DO 200 I=1,13
  200 SUM=SUM+P*V(I)/(R*(TMIX(I)+B*P*REALGS))
      DO 220 I=1,13
  220 WGOLD(I)=WG(I)
      DO 240 I=1,13
  240 WG(I)=W*(1.-SIG3)*P*V(I)/(R*SUM*(TMIX(I)+B*P*REALGS))
      IF(MDEX.EQ.2) GO TO 243
      DO 242 I=1,13
  242 WGOLD(I)=WG(I)
C CALCULATES FLOW RATES, F(I) IN LBM/SEC, BETWEEN GAS LUMPS.
  243 MDEX=2
      F(1)=(WGOLD(1)-WG(1))/DTIME
      DO 300 I=2,12
  300 F(I)=(WGOLD(I)-WG(I))/DTIME+F(I-1)
      DO 24 I=1,12
      IF(F(I).LT.0.) GO TO 22
      TG(I)=TGAO(I)
      IF (I.GT.4 .AND. I.LT.10) TG(I)=TGAO(I)-(TMA(5)-TMA(9))/8.+
     1F(I)*DTIME*(TMA(5)-TMA(9))/8./WGOLD(I)
      GO TO 24
   22 TG(I)=TGAO(I+1)
      IF(I.GT.3 .AND. I.LT.9) TG(I)=TGAO(I+1)+(TMA(5)-TMA(9))/8.+
     1F(I)*DTIME*(TMA(5)-TMA(9))/8./WGOLD(I+1)
   24 CONTINUE
C CALCULATE AVERAGE FLOW RATES FOR EACH GAS LUMP.
      FAVG(1)=F(1)
      DO 30 I=2,12
   30 FAVG(I)=(F(I-1)+F(I))/2.
      FAVG(13)=F(13)
C MUST BE A MISTAKE, COMMON F(12) IN HEATX !!!!!!!!!!!!!!!
C REVISE GAS TEMPERATURES TO ACCOUNT FOR MIXING.
      IF (F(1).LT.0.0) TMIX(1)=(WGOLD(1)*TMIX(1)-F(1)*DTIME*
     1TG(1))/WG(1)
      DO 100 I=2,12
C PAGE 53 (59 IN MANUAL PDF)---------------------------
  100 TMIX(I)=(WGOLD(I)*TMIX(I)+(F(I-1)*TG(I-1)-F(I)*TG(I))*
     1DTIME)/WG(I)
      IF (F(12).GT.0.0) TMIX(13)=(WGOLD(13)*TMIX(13)+F(12)*
     1DTIME*TG(12))/WG(13)
  110 CONTINUE
C CALCULATE REYNOLDS NO S, REYNO(I); HEAT TRANSFER AREAS, AHT(I);
C HEAT TRANSFER COEFFICIENTS, HA(I) BETWEEN GAS + METAL LUMPS.
C HA(I) HAS UNITS BTU/SEC-DEG.R
      REYNO(1)=DH(1)*ABS(F(1))/(VIS(1)*ACS(1))+.0001
      AHT(1)=3.1416*DIAMD*(X(1)-X1MIN)+16.6
      HA(1)=COND(1)/DH(1)*(1.871E-3*DH(1)*ABS(F(1))/(VIS(1)*
     1ACS(1))+15.)*AHT(1)
C TO MAKE EXPANSION SPACE PROCES ADIABATIC SET HA 1 =0.0
C     HA(1)=  0.0
C
C  QUENCHING IN EXP. SPACE--RADIATION .70 EFFECTIVE--HCONV=5/3600/144 MI
      I=1
      XL(1)=X(1)-X1MIN+.01
      XL(13)=V(13)/AP+.01
      FOA=0.7
      IF(ABS(TMIX(1)-TMA(1)).LE.1.0)  GO TO 42
      QRAD=0.173*FOA*((TMIX(1)/100.)**4-(TMA(1)/100.)**4)
      HRAD=QRAD/(TMIX(1)-TMA(1))/3600./144.
      GO TO 43
   42 HRAD=0.0
   43 CONTINUE
      PR=CP*VIS(I)/COND(I)
      HCONV=0.023* REYNO(I)**0.8*PR**0.4*COND(I)/DH(I)
      IF(REYNO(I).GT.2100.0.AND.REYNO(I).LT.10000.) HCONV=HCONV*(1.+(DH(
     1I)/XL(I))**0.7)
      GRAETZ=REYNO(I)*PR*DH(I)/XL(I)
      IF(REYNO(I).LE.2100.0.AND.GRAETZ.GT.10.) HCONV=1.86*GRAETZ**0.333*
     1COND(I)/DH(I)
      IF(REYNO(I).LE.2100.0.AND.GRAETZ.LE.10.)  HCONV=5.0/3600./144.
      IF(I.EQ.13)  GO TO 44
C ACCOUNT FOR HEAT TRANSFER TO INSULATED PORTION OF HEATER TUBES
C (ADJACENT TO EXPANSION SPACE) ALONG WITH HEAT TRANSFER TO EXPANSION
C SPACE WALL
      XL1=1.321
      REYNO1=DH(2)*ABS(F(1))/8.0/(VIS(1)*ACS(2))+.0001
      AHT1=14.955*XL1
      HCONV1=0.023*REYNO1**0.8*PR**0.4*COND(1)/DH(2)
      HA1=(HRAD+HCONV1)*AHT1
      HA(1)=(HRAD+HCONV)*AHT(1)+HA1
C
      DO 40 I=2,12
      REYNO(I)=DH(I)*ABS(F(I-1)+F(I))/2./8./(VIS(I)*ACS(I))+.0001
      IF (I.LT.5.OR.I.GT.9) HA(I)=COND(I)/DH(I)*(1.871E-3*
     1DH(I)*ABS((F(I-1)+F(I))/2.)/8./
     2(VIS(I)*ACS(I))+15.)*AHT(I)*8.
      IF (I.GE.5.AND.I.LE.9) HA(I)=COND(I)/DH(I)*(0.06071*
     1DH(I)*ABS((F(I-1)+F(I))/2.)/8./
     2(VIS(I)*ACS(I))+3.7)*AHT(I)*8.
   40 CONTINUE
      HA(4)=0.0
      DO 41 I=5,9
   41 HA(I)=RHCFAC*HA(I)
      REYNO(13)=DH(13)*ABS(F(12))/(VIS(13)*ACS(13))+.0001
      AHT(13)=6.2832*DIAMD*(RODYMX-RODY)+5.5
      HA(13)=COND(13)/DH(13)*(1.871E-3*DH(13)*ABS(F(12))/
     1(VIS(13)*ACS(13))+15.)*AHT(13)
C TO MAKE COMPRESSION SPACE PROCES ADIABATIC SET HA 13 =0.0
C     HA(13)=  0.0
      I=13
      IF(I.EQ.13) GO TO 42
   44 HA(13)=HCONV*AHT(13)
C PAGE 54 (60 IN MANUAL PDF)-----------------
C CALCULATE RATES OF HEAT TRANSFER BETWEEN GAS AND METAL REVISE
C GAS TEMPERATURES TO ACCOUNT FOR HEAT TRANSFER TO OR FROM METAL.
      DO 120 N=1,13
      HAWC=-DTIME*HA(N)/(WG(N)*CP)
      FCTR=1.-EXP(HAWC)
      DTGA=(TMAO(N)-TMIX(N))*FCTR
      Q(N)=-DTGA*WG(N)*CP/DTIME
      WFCTR=FCTR*WG(N)
      WDTGA=DTGA*WG(N)
      SUMNUM(N)=SUMNUM(N)+WDTGA
      SUMDEN(N)=SUMDEN(N)+WFCTR
  120 TGA(N)=TMIX(N)+DTGA
      IF (JCYCLE.LT.(NOCYC-5)) GO TO 64
C CALCULATE DENSITY,FRICTION FACTORS, AND DELTA PRESSURES
C FOR EACH LUMP.
      DO 130 I=1,13
  130 DNSTY(I)=P/(R*(TGA(I)+B*P*REALGS))
      DO 140 I=2,12
      IF (REYNO(I).LT.1500.AND.I.LT.5.OR.I.GT.9) FRICT(I)=16./REYNO(I)
      IF (REYNO(I).GE.1500.AND.I.LT.5.OR.I.GT.9) FRICT(I)=0.046/
     1REYNO(I)**0.2
C CHANGED. tHE LINE BELOW HAS A SMALL RISK OF FAILING.
C      IF (REYNO(I).GE.1500..AND.I.LT.5.OR.I.GT.9) FRICT(I)=0.046/
      IF (I.GE.5.AND.I.LE.9) FRICT(I)=0.96*EXP(-0.0190*REYNO(I))+0.54
  140 DELP(I)=FRICT(I)*FAVG(I)*ABS(FAVG(I))*XL(I)/
     1(32.2*DH(I)*DNSTY(I)*ACS(I)**2*6.*64.)
C PRESSURE DROP CALCULATIONS
      DO 141 I=1,40
  141 DPX(I)=GAMMA
      ACSEXP=V(1)/AD*3.1416*2.75/8.
      ACSCOM=V(13)/AP*3.1416*2.75/8.
      J=1
      IF(F(1).LE.0.0)  J=2
      Z=P
      F01=F(1)
      F02=F(2)
      F04=F(4)
      F09=F(9)
      F12=F(12)
      CALL XDELP(J,ACSEXP, ACS(02),DH(02),XL(02),DNSTY(02),VIS(02),F01
     X,Z,COEFX(01),DPX(01),1)
C
      CALL XDELP(3,ACS(02),ACS(02),DH(02),XL(02),DNSTY(02),VIS(02),F01
     X,Z,COEFX(02),DPX(02),2)
      CALL XDELP(4,ACS(02),ACS(02),DH(02),XL(02),DNSTY(03),VIS(03),F02
     X,Z,COEFX(03),DPX(03),3)
C
      DPFRIC=0.0
      XL2=XL(2)
      DO 142 I=2,12
      IF(I.EQ.2)  XL(2)=XL2+1.321
      IF(I.LT.5.OR.I.GT.9) KTYPE=10
      IF(I.GE.5.AND.I.LE.9) KTYPE=11      
      CALL XDELP(KTYPE,ACS(I),ACS(I),DH(I),XL(I),DNSTY(I),VIS(I),FAVG(I
     X),Z,COEFX(I+2),DPX(I+02),(I+2))
C    X,Z,COEFX(I+2),DPX(I+02),(I+2))  MUST BE WRONG (HARRIE) !!!!!!!!
      XL(2)=XL2
      DELP(I)=DPX(I+2)
      DPFRIC=DPFRIC+DELP(I)
  142 CONTINUE
C
      IF(F(4).GT.0.0) GO TO 61
      CALL XDELP(1,ACS(04),.622000,DH(04),XL(04),DNSTY(04),VIS(04),F04
     X,Z,COEFX(15),DPX(15),15)
      GO TO 62
   61 CONTINUE
      CALL XDELP(1,.622000,ACS(10),DH(10),XL(10),DNSTY(10),VIS(10),F09
     X,Z,COEFX(16),DPX(16),16)
   62 CONTINUE

( PAGE 55, 61 IN PDF -------------------------

      CALL XDELP(5,ACS(12),.043370,DH(12),XL(12),DNSTY(12),VIS(12),F12
     X,Z,COEFX(17),DPX(17),17)
      J=2
      IF(F(12).LE.0.0)  J=1
      CALL XDELP(J,ACS(12),ACSCOM,DH(12),XL(12),DNSTY(12),VIS(12),F12
     X,Z,COEFX(18),DPX(18),18)
      DPSUM=0.0
      DO 63 I=1,17
      DPSUM=DPSUM+DPX(I)
   63 CONTINUE
      DELP(1)=DPFRIC
      DELP(13)=DPSUM
C
      DELP(02)=DELP(02)+DPX(01)+DPX(02)
      DELP(03)=DELP(03)+DPX(03)
      DELP(04)=DELP(04)+DPX(15)
      DELP(10)=DELP(10)+DPX(16)
      DELP(12)=DELP(12)+DPX(17)+DPX(18)
  174 CONTINUE
C CALCULATE HEAT CONDUCTION (CALL HEAT CONDUCTION SUBROUTINE)
      IF (NCOND.NE.3) CALL CNDCT (TMA,QCREG ,QCYL,QCONDD,QCAN,QSHTTL,
     1MWGAS)
      NCOND=3
   64 CONTINUE
C REVISE REGENERATOR METAL TEMPERATURES.
      DO 177 J=5,9
  177 TMA(J)=TMA(J)+Q(J)*DTIME/CPM(J)
C  CONVERT Q(J)'S TO F1-LBF & SUM UP FOR EACH COMPONENT.
      QTOT=0.0
      QEX=0.0
      QHEAT=0.0
      QHETN=0.0
      QHETP=0.0
      QREGE=0.0
      QRGP=0.0
      QRGN=0.0
      QCOOL=0.0
      QCOLP=0.0
      QCOLN=0.0
      QCON=0.0
      DO 179 J=1,13
      QADD=Q(J)*DTIME*777.3
      IF (J.EQ.1) QEX=QEX+QADD
      IF (J.GE.2.AND.J.LE.4) QHEAT=QHEAT+QADD
      IF (J.GE.2.AND.J.LE.4.AND.QADD.LT.0.0) QHETN=QHETN+QADD
      IF (J.GE.2.AND.J.LE.4.AND.QADD.GT.0.0) QHETP=QHETP+QADD
      IF (J.GE.5.AND.J.LE.9) QREGE=QREGE+QADD
      IF (J.GE.5.AND.J.LE.9.AND.QADD.GT.0.0) QRGP=QRGP+QADD
      IF (J.GE.5.AND.J.LE.9.AND.QADD.LT.0.0) QRGN=QRGN+QADD
      IF (J.GE.10.AND.J.LE.12) QCOOL=QCOOL+QADD
      IF (J.GE.10.AND.J.LE.12.AND.QADD.LT.0.0) QCOLN=QCOLN+QADD
      IF (J.GE.10.AND.J.LE.12.AND.QADD.GT.0.0) QCOLP=QCOLP+QADD
      IF (J.EQ.13) QCOM=QCOM+QADD
  179 QTOT=QTOT+QADD
      DO 185 J=1,5
  185 QR(J)=QR(J)+Q(J)*DTIME/CPM(J)
      QEXP=QEXP+QEX
      IF (QEX .LT.0.0) QEXPN=QEXPN+QEX
      IF (QEX .GT.0.0) QEXPP=QEXPP+QEX
      QHEATR=QHEATR+QHEAT
      QHEATN=QHEATN+QHETN
      QHEATP=QHEATP+QHETP
      QREGEN=QREGEN+QREGE
      QREGN=QREGN+QRGN
      QREGP=QREGP+QRGP
      QCOOLR=QCOOLR+QCOOL
C PAGE 56 (62 IN PDF) -----------------------------

      QCOOLP=QCOOLP+QCOLP
      QCOOLN=QCOOLN+QCOLN
      QCOMP=QCOMP+QCOM
      IF (QCOM.GT.0.0) QCOMPP=QCOMPP+QCOM
      IF (QCOM.LT.0.0) QCOMPN=QCOMPN+QCOM
      QCNDRI=QCNDRI+QCREG   *DTIME*777.3
      QCNDRO=QCNDRO+QCREG   *DTIME*777.3
      QCNDCL=QCNDCL+QCYL*DTIME*777.3
      QCNDD=QCNDD+QCONDD*DTIME*777.3
      QCNDCN=QCNDCN+QCAN*DTIME*777.3
      QSHTL=QSHTTL*DTIME*777.3+QSHTL
  350 CONTINUE
      ENFRTH=ENFRTH-AMIN1(F(4)*CP,0.0)*TG(4)*DTIME
      ENFHTR=ENFHTR+AMAX1(F(4)*CP,0.0)*TG(4)*DTIME
      ENFCTR=ENFCTR-AMIN1(F(9)*CP,0.0)*TG(9)*DTIME
      ENFRTC=ENFRTC+AMAX1(F(9)*CP,0.0)*TG(9)*DTIME
      RETURN
      END

C SUBROUTINE CYCL
      
C THIS SUBROUTINE IS CALLED JUST ONCE PER CYCLE--WHEN THE CYCLE IS COMLETE.
C IT CALCULATES--NET HEAT IN AND OUT, NET WORK PER CYCLE, INDICATED POWERS
C AND EFFICIENCY, ETC.  HEAT AND WORK QANTITIES ARE STORED FOR CALCULATION
C QUANTITIES ARE RESET TO ZERO IN PREPARATION FOR NEXT CYCLE.
      SUBROUTINE CYCL (TIME,UTOTAL,STROKV,AVGWSP)
      DIMENSION VAR(106),STORE(11,86)
      COMMON DTIME,F(12),WG(13),DELP(13),P,PE,PC
      DATA J,K,STORE,OLDTIM/2*0,947*0./
      DATA  JIP1/0/
      COMMON /CYC/QEXP,QHEATR,QREGEN,QCOOLR,QCOMP,ENFRTH,ENFHTR,
     1ENFCTR,ENFRTC,WRKP,WRKEXP,
     2WRKCMP,WRKLEX,WRKLCM,WRKLP,WRKH,SUMNUM(13),SUMDEN(13),ENTH,
     3WDISP,WPIST,WTPIST,QEXPN,QEXPP,QHEATN,QHEATP,
     4QCOOLN,QCOOLP,QCOMPN,QCOMPP,QREGN,QREGP,WDISPP,WDISPN,WPISTP,
     5WPISTN,QCNDRI,QCNDRO,QCNDCL,QCNDD,QCNDCN, QSHTL,TGEXPA,TGCMPA
      COMMON  /TAVCYC/TGCYC(12),TGACYC(13),TMCYC(13)
      COMMON /JNDEX/JIP,NOCYC
      EQUIVALENCE (VAR(1),QEXP)
      J=J+1
      K=K+1
      JIP1=JIP1+1
      IF(OLDTIM.GT.TIME) OLDTIME=0.
      DELTIM=TIME-OLDTIM
      OLDTIM=TIME
      QIN=-QEXP-QHEATR+QCNDRI+QCNDCL+QCNDD+QCNDCN+QSHTL-WRKLT/2.
      QINB=-QEXP-QHEATR
      QOUT=QCOOLR+QCOMP+QCNDRO+QCNDCL+QCNDD+QCNDCN+QSHTL
      WRKTOT=WRKEXP+WRKCMP
      WRKLT=WRKLEX+WRKLCM
      WRKBAS=WRKTOT+WRKLT
      EFFTOT=WRKTOT/QIN
      EFFP=WRKP/QIN
      REFF1=ENFRTH/ENFHTR
      REFF2=ENFCTR/ENFRTC
      PWRHP=WRKH/(550.*DELTIM)

C PAGE 57 (63 IN PDF)  ----------------------------------------

      PWRKW=PWRHP*.7457
      FREQ=1./DELTIM
      EN3PM=STROKV*FREQ*60.
C STORE VARIABLES FOR CALCULATION OF AVERAGES OVER 5 CYCLES.
      STORE(K,1)=QIN
      STORE(K,2)=QOUT
      STORE(K,3)=WRKEXP
      STORE(K,4)=WRKCMP
      STORE(K,5)=WRKTOT
      STORE(K,6)=EFFTOT
      STORE(K,7)=REFF1
      STORE(K,8)=REFF2
      STORE(K,9)=QREGEN
      STORE(K,10)=PWRHP
      STORE(K,11)=FREQ
      STORE(K,12)=WRKLT
      STORE(K,13)=ENTH
      STORE(K,14)=UTOTAL
      STORE(K,15)=EN3PM
      STORE(K,16)=AVGWSP
      STORE(K,17)=WDISP
      STORE(K,18)=WPIST
      STORE(K,19)=WTPIST
      STORE(K,20)=WRKLP
      STORE(K,21)=QEXPN
      STORE(K,22)=QEXPP
      STORE(K,23)=QHEATN
      STORE(K,24)=QHEATP
      STORE(K,25)=QCOOLN
      STORE(K,26)=QCOOLP
      STORE(K,27)=QCOMPN
      STORE(K,28)=QCOMPP
      STORE(K,29)=QREGN
      STORE(K,30)=QREGP
      STORE(K,31)=WDISPP
      STORE(K,32)=WDISPN
      STORE(K,33)=WPISTP
      STORE(K,34)=WPISTN
      STORE(K,35)=QCNDRI
      STORE(K,36)=QCNDRO
      STORE(K,37)=QCNDCL
      STORE(K,38)=QCNDD
      STORE(K,39)=QCNDCN
      STORE(K,40)=TGEXPA
      STORE(K,41)=TGCMPA
      STORE(K,42)=QSHTL
      STORE(K,43)=QINB
      STORE(K,44)=WRKBAS
      DO 90 L=1,13
      IF(L.EQ.13) GO TO 91
      M=L+48
      STORE(K,M)=TGCYC(L)
   91 M=L+60
      STORE(K,M)=TGACYC(L)
      M=L+73
   90 STORE(K,M)=TMCYC(L)
  100 CONTINUE
      IF (JIP.GT.0.AND.JIP1.NE.NOCYC) GO TO 301
      IF (JIP.GT.0) WRITE (6,202)
  202 FORMAT ('LAST CYCLE')
      WRITE (6,200) QIN,QOUT,WRKEXP,WRKCMP,WRKTOT,EFFTOT,REFF1,REFF2,
     *QREGEN,PWRHP,FREQ,WRKLT,     ENTH,UTOTAL,
     *EN3PN,AVGWSP,       WDISP,WPIST,WTPIST,WRKLP,
     *QEXPN,QEXPP,QHEATN,QHEATP,QCOOLN,QCOOLP,QCOMPN,QCOMPP,
     *QREGN,QREGP,WDISPP,WDISPN,WPISTP,WPISTN,QCNDRI,QCNDRO,
     *QCNDCL,QCNDD,QCNDCN,TGEXPA,TGCMPA,QSHTL,QINB,WRKBAS

C PAGE 58 (64 IN PDF) ----------------------------------

  200 FORMAT (/' QIN=',F8.3,1X,'QOUT=',F8.3,1X,'WRKEXP=',F8.3,1X,
     1'WRKCMP=',F8.3,1X,'WRKTOT=',F8.3,1X,
     2'EFFTOT=',F8.3,1X,'REFF1=',F8.3,1X,'REFF2=',F8.3,1X/
     3' QREGEN=',F8.3,' PWRHP=',F8.3,' FREQ=',F8.3,
     4' WRKLT=',F8.3,              'ENTH=',F8.3,' UTOTAL=',
     5F10.3/' EN3PN=',F10.3,' AVGWSP=',F9.3,
     6' WDISP=',F8.3,' WPIST=',F8.3,' WTPIST=',F8.3,' WRKLP=',F8.3/,
     7' QEXPN=',F8.3,' QEXPP=',F8.3,' QHEATN=',F8.3,' QHEATP=',F8.3,
     8' QCOOLN=',F8.3,' QCOOLP=',F8.3,' QCOMPN=',F8.3,' QCOMPP=',F8.3/,
     9' QREGN=',F9.3,' QREGP=',F8.3,' WDISPP=',F8.3,' WDISPN=',F8.3,
     1' WPISTP=',F8.3,' WPISTN=',F8.3,' QCNDRI=',F7.3,' QCNDRO=',F7.3,/
     2' QCNDCL=',F8.3,' QCNDD=',F8.3,' QCNDCN=',F8.3,' TGEXPA=',F8.3,
     3' TGCMPA=',F8.3,' QSHTL=',F8.3/,'QINB =',F8.3,' WRKBAS =',F8.3)
  301 CONTINUE
      DO 300 I=1,66
  300 VAR(I)=0.0
      IF(J.EQ.5) GO TO 400
      RETURN
  400 K=0
      J=0
      DO 450 KK=1,86
      STORE(11,KK)=STORE(1,KK)/5.
      DO 450 IK=2,5
  450 STORE(11,KK)=STORE(11,KK)+STORE(IK,KK)/5.
  575 CONTINUE
      IF (JIP.GT.0.AND.JIP1.NE.NOCYC) GO TO 601
      WRITE (6,600)
      WRITE (6,200) (STORE(11,KK),KK=1,44)
      IF (JIP.GT.0) WRITE (6,204)
  204 FORMAT (' AVERAGE TEMPERATURES OVER LAST CYCLE')
      WRITE(6,201)  TGCYC,TGACYC,TMCYC
  201 FORMAT ('  TGCYC=',4X,12F8.0/,' TGACYC=',13F8.0/,'  TMCYC=',13F8.0
     X/)
  600 FORMAT (1HD,'AVERAGE VALUES OF LAST 5 CYCLES')
  601 RETURN
      END
      
C SUBROUTINE CNDCT

      SUBROUTINE CNDCT (TMA,QCREG, QCYL,QCONDD,QCAN,QSHTTL,MWGAS)
      DIMENSION TMA(13),                     TCYL(3),TCAN(3),
     1        ACAN(2),         DCAN(2)
      COMMON /CYC/ QEXP,DHEATR,QREGEN,QCOOLR,QCOMP,ENFRTH,ENFHTR,
     1ENFCTR,ENFRTC,WRKP,WRKEXP,
     2WRKCMP,WRKLEX,WRKLCM,WRKLP,WRKH,SUMNUM(13),SUMDEN(13),ENTH,
     3WDISP,WPIST,WTPIST,QEXPN,QEXPP,QHEATN,QHEATP,
     4QCOOLN,QCOOLP,QCOMPN,QCOMPP,QREGN,QREGP,WDISPP,WDISPN,WPISTP,
     5WPISTN,QCNDRI,QCNDRO,QCNDRCL,QCNDD,QCNDCN, QSHTL,TGEXPA,TGCMPA
      DATA TR0,TR1,TR2,TCYL/1760.,1172.,585.,1760.,1255.,975./,
     1     RI2,R2,R3,R4,XL3,XL4/1.38,1.62,1.56,1.53,1.126,0.400/,
     1     RI1,R0,R1,XL1,XL2,/0.454,0.521,0.495,0.400,0.470/,
     1       ACONDD,DCONDD/0.5875,1.716/,
     1       TCAN/1400.,1280.,1100./,
     1       ACAN,DCAN/1.154,1.196,0.5,0.75/,
     1       CLEAR,DISPD,STROKE/.010,2.75,1.23/
C ACONDR,ACYL,ACONDD,ACAN---AREA S IN IN2
C DCONDR,DCYL,DCONDD,DCAN---DISTANCE IN IN

C PAGE 59 (65 IN PDF) ------------------

C TCYL,TCAN--- TEMPERATURES IN DEGREES R
C CLEAR--CLEARANCE BETEEN DISPLACER AND CYLINDER, IN.
C DISPD---DISPLACER DIAMETER, IN
C STROKE---DISPLACER STROKE, IN.

C CALCULATE HEAT CONDUCTED THROUGH CYLINDER,BTU/SEC
      TCAVG1=(TCYL(1)+TCYL(2))/2.
      TCAVG2=(TCYL(2)+TCYL(3))/2.
      THCND1=(.00387*TCAVG1+6.793)/(3600.*12.)
      IF (TCAVG1.GT.1060.) THCND1=(.00467*TCAVG1+5.932)/(3600.*12.)
      THCND2=(.00387*TCAVG2+6.793)/(3600.*12.)
      IF (TCAVG2.GT.1060.) THCND2=(.00467*TCAVG2+5.932)/(3600.*12.)
      QCYL1=6.283*THCND1*RI2*(R3-R2)*(TCYL(1)-TCYL(2))/(XL3*ALOG(((R3-
     1RI2)*(R2+RI2))/((R3+RI2)*(R2-RI2))))
      QCYL2=6.283*THCND2*RI2*(R4-R3)*(TCYL(2)-TCYL(3))/(XL4*ALOG(((R4-
     1RI2)*(R3+RI2))/((R4+RI2)*(R3-RI2))))
      QCYL=(QCYL1+QCYL2)/2.
C CALCULATE HEAT CONDUCTED THROUGH DISPLACER,BTU/SEC
      TGAVG=(TGEXPA+TGCMPA)/2.
      THCOND=(.00387*TGAVG+6.793)/(3600.*12.)
      IF (TGAVG.GT.1060.) THCOND=(.00467*TGAVG+5.932)/(3600.*12.)
      QCONDD=ACONDD*THCOND*(TGEXPA-TGCMPA)/DCONDD
C CALCULATE HEAT CONDUCTED THROUGH EXTERNAL CAN,BTU/SEC
      QCAN=0.0
      DO 300 I=1,2
      TCAVG=(TCAN(I)+TCAN(I+1))/2.
      THCOND=(.00387*TCAVG+6.793)/(3600.*12.)
      IF (TCAVG.GT.1060.) THCOND=(.00467*TCAVG+5.932)/(3600.*12.)
  300 QCAN=QCAN+ACAN(I)*THCOND*(TCAN(I)-TCAN(I+1))/(DCAN(I)*2.)
C CALCULATE HEAT CONDUCTION THROUGH REGENERATOR,BTU/SEC
      TRAVG1=(TR0+TR1)/2.
      TRAVG2=(TR1+TR2)/2.
      THCND1=(.00387*TRAVG1+6.793)/(3600.*12.)
      IF (TRAVG1.GT.1060.) THCND1=(.00467*TRAVG1+5.932)/(3600.*12.)
      THCND2=(.00387*TRAVG2+6.793)/(3600.*12.)
      IF (TRAVG2.GT.1060.) THCND2=(.00467*TRAVG2+5.932)/(3600.*12.)
      QCRIN=3.142*(R0**2-RI1**2)*THCND1 *(TR0-TR1)/XL1
      QCROUT=6.283*THCND2*RI1*(R1-R0)*(TR0-TR1)/(XL2*ALOG(((R1-RI1)*
     1(R0+RI1))/((R1+RI1)*(R0-RI1))))
      QCREG=(QCRIN+QCROUT)/2.*8.
C CALCULATE SHUTTLE HEAT TRANSFER, BTU/SEC
      THCNDG=2.481E-09*(TCYL(1)+TCYL(3))/2.+1.263E-06
      IF (MWGAS.EQ.4) THCNDG=1.962E-09*( TCYL(1)+TCYL(3))/2.+1.115E-6
      QSHTTL=THCNDG*3.1416*DISPD*STROKE**2*(TCYL(1)-TCYL(3))/(8.0*CLEAR*
     1(XL3+XL4        ))
      RETURN
      END

C SUBROUTINE XDELP
      SUBROUTINE  XDELP(K,AREAIN,AREAOT,DHX,XLGTH,ROX,VISX,FLOWIN,PXIN,
     XCOEF,DP,M)
C INDENTIFICATION OF TYPE OF LOSSES
C KTYPE=1 CONTRACTION--KTYPE=2 EXPANSION--KTYPE=3 90-TURN --KTYPE4 I
C KTYPE=5 OTHER TURN --KTYPE=6 FRICTION-TUBE--KTYPE=7 FRICTION-REGENE
C KTYPE=8 FRICTION+MOMENTUM FOR TUBE--ISOTHERMAL
C KTYPE=9 FRICTION+MOMENTUM FOR REGENRATOR--ISOTHERMAL

C PAGE 60 (66 IN PDF) ----------------------

      DIMENSION  COEFX(40),DPX(40)
      DATA  GRAV/32.2/,  PSI/0.006945/
      GAMMA=DP
      IF (ABS(FLOWIN).LE.0.0) RETURN
      PR=0.
      POUT=0.
      DPMON=0.
C DPMON COULD BE OPMON OR OPMOM OR OPNON ????
C VAR USED IN PROGRAM?
      KTYPE=K
      AIN=AREAIN/144.
      AOUT=AREAOT/144.
      DE=DHX/12.
      KLG=XLGTH/12.
      RO=ROX*144.*12.
      VISC=VISX*12
      FLOW=ABS(FLOWIN)/8.0
      PIN=PXIN
      ROMEAN=RO
      IF(AIN-AOUT)  30,30,31
   30 AMIN=AIN
      A=AIN/AOUT
      GO TO 32
   31 AMIN=AOUT
      A=AOUT/AIN
   32 CONTINUE
      IF(A.LT.0.001)  A=0.0
      VEL=FLOW/AMIN/RO
      VELHD=RO*VEL**2/(2.0*GRAV)
      RE=RO*VEL*DE/VISC
      GO TO (1,2,3,4,5,6,7,6,7,6,7),KTYPE
C
    1 CONTINUE
C KTYPE=1 CONTRACTION
      IF(RE-3000.0)  12,12,11
   11 COEF=-0.4*A+0.5
      GO TO 100
   12 COEF=-0.4*A+1.0
      GO TO 100
C
    2 CONTINUE
C KTYPE=2 EXPANSION
      IF(RE-3000.0)  21,21,22
   21 COEF=(1.0-2.6*A+1.005*A**2)*(-1.0)
      GO TO 100
   22 COEF=(1.0-2.0928*A+0.996*A**2)*(-1.0)
      GO TO 100
C
    3 CONTINUE
      COEF=0.22
      GO TO 102
C
    4 CONTINUE
      COEF=3.0*0.22
      GO TO 102
C
    5 CONTINUE
      COEF=0.11
      GO TO 102
C
    6 CONTINUE
C KTYPE=6 FRICTION-TUBE
      IF(RE.LT.1500.) COEF=16.0/RE*(XLG/(DE/4.0))
      IF(RE.GE.1500.)  COEF=0.046/RE**0.2*(XLG/(DE/4.0))
      GO TO 102
C
    7 CONTINUE
C KTYPE=7 FRICTION-REGENERATOR

C PAGE 61 (67 IN PDF) ---------------

      COEF=0.96*EXP(-0.0190*RE)+0.54
      COEF=COEF*(XLG/(DE/4.0))
      GO TO 102
  100 COEF=1.0-A**2+COEF
  102 DP=COEF*VELHD*PSI
      IF(KTYPE.LT.8) GO TO 101
C KTYPE=8 FRICTION+MOMENTUM --TUBE   KTYPE=9 SAME FOR REGENERATOR
C KTYPE=10 --TUBE      MON.+FRIC.  KTYPE=11--FOR REGENERATOR  ADI
      IF(KTYPE.GE.10) GAMMA1=GAMMA
      KMA1=VEL/SQRT(GAMMA*PIN*144.0/RO+GRAV)
      KMA2=1.0/SQRT(1.0/XMA1**2-GAMMA*COEF)
      FMA21=1.0
      IF(ABS((XMA2-XMA1)/XMA1).LE.0.001) GO TO 52
   51 XMA0=(XMA1+XMA2)*0.5
      FMA21=((GAMMA-1.0)*XMA0**2+2.0)/((GAMMA-1.0)*XMA1**2+2.0)
      IF(KTYPE.LT.10)  FMA21=1.0
      IF(KTYPE.LT.10)  GAMMA1=1.0
      XMA2=1.0/SQRT(1.0/XMA1**2-GAMMA*COEF+GAMMA*(GAMMA1+1.0)/2.0/GAMMA1
     1*ALOG((XMA1/XMA0)**2*FMA21))
      IF(ABS(XMA2-XMA0).GT.0.001)  GO TO 51
   52 TR=1.0/FMA21
      VR=XMA2/XMA1*SQRT(1.0/FMA21)
      PR=XMA1/XMA2*SQRT(1.0/FMA21)
      POUT=PR*PIN
      DP=PIN-POUT
  101 CONTINUE
      IF(FLOWIN.LE.0.0) DP=(-1.)*DP
      COEFX(M)=COEF
      DPX(M)=DP
      RETURN
      END

'==================================================================

Here you find a list of the input and output variables with their meanings.

'input variables
DATA "COEFF"
DATA "leage coefficient for flow between buffer and compression spaces"
DATA "FACT1"
DATA "FACT1 and FACT2 are convergence constants used in regenerator 
metal temperature convergence scheme" DATA "FACT2" DATA "FACT2 and FACT1 are convergence constants used in regenerator
metal temperature convergence scheme" DATA "FIPCYC" DATA "number of iterations per cycle" DATA "IPRINT" DATA "number of iterations between printouts (if JIP=0)" DATA "ITMPS" DATA "1, printout gas and metal control volume temperature at every iteration
0, do not do the above" DATA "JIP" DATA ">0, short-form printout
=0, long-form printout" DATA "MWGAS" DATA "4, use helium working gas
not 4, use hydrogen working gas" DATA "NOCYC" DATA "number of engine cycles to be run" DATA "NOEND" DATA"number of cycle at wich convergence scheme is turned on" DATA "NSTRT" DATA "number of cycle at which regenerator temperature convergence
scheme is turned on" DATA "OMEGA" DATA "drive speed, Hz" DATA "P" DATA "initial pressure, lbf/inch^2 (N/m^2) (approximately equal to mean
pressure for GPU)" DATA "P3" DATA "initial buffer space pressure, lbf/inch^2 (N/m^2)" DATA "REALGS" DATA "1, use real-gas equation of state
0, use ideal-gas equation of state" DATA "RHCFAC" DATA "regenerator heat-transfer coefficient multiplication factor" DATA "TG" DATA "12 control volume interface gas temperatures, degrees R (Kelvin)" DATA "TGA" DATA "13 control volume gas temperatures, degrees R (Kelvin)" DATA "TM" DATA "13 metal boundary temperatures, degrees R (Kelvin)" DATA "T3" DATA "buffer space temperature (assumed constant), degrees R (Kelvin)" 'output variables DATA "AVGWSP" DATA "mean working-space pressure, lbf/inch^2 (N/m^2)" DATA "EFFTOT" DATA "efficiency, WRKTOT/QIN" DATA "ENTH" DATA "net enthalpy flow from working space to buffer space
over one cycle, ft-lbf (J)" DATA "EN3PM" DATA "rate of displacement of working-space gas, inch^3/min (cm^3/min)" DATA "FREQ" DATA "rotational frequency of crankshaft, Hz" DATA "PWRHP" DATA "indicated power -WRKTOT/cycle period, hp (kW)" DATA "QCNDCL" DATA "heat conduction through cylinder, ft-lbf/cycle (J/cycle)" DATA "QCNDCN" DATA "heat conduction through insulation container, ft-lbf/cycle (J/cycle)" DATA "QCNDD" DATA "heat conduction through displacer, ft-lbf/cycle (J/cycle)" DATA "QCNDRI" DATA "heat conduction through regenerator casing, ft-lbf/cycle (J/cycle)" DATA "QCNDRO" DATA "heat conduction through regenerator casing, ft-lbf/cycle (J/cycle)" DATA "QCOMPN" DATA "heat from compression-space wall to gas, ft-lbf/cycle (J/cycle)" DATA "QCOMPP" DATA "heat from gas to compression-space wall, ft-lbf/cycle (J/cycle" DATA "QCOOLN" DATA "heat from cooler to gas, ft-lbf/cycle (J/cycle)" DATA "QCOOLP" DATA "heat from gas to cooler, ft-lbf/cycle (J/cycle)" DATA "QEXPN" DATA "heat from metal wall to gas in expansion space, ft-lbf/cycle (J/cycle)" DATA "QEXPP" DATA "heat from gas to metal wall in expansion space ft-lbf/cycle (J/cycle)" DATA "QHEATN" DATA "heat from heater to gas, ft-lbf/cycle (J/cycle)" DATA "QHEATP" DATA "heat from gas to heater, ft-lbf/cycle (J/cycle)" DATA "QIN" DATA "heat into gas in expansion space and heater tubes over one cycle,
ft-lbf (J)" DATA "QOUT" DATA "heat out of gas in cooler tubes and compression space over
one cycle, ft-lbf (J)" DATA "QREGEN" DATA "net heat flow from gas to regenerator metal over one cycle,
ft-lbf/cycle (J/cycle)" DATA "QREGN" DATA "heat from regenerator metal to gas, ft-lbf/cycle (J/cycle)" DATA "QREGP" DATA "heat from gas to regenerator metal, ft-lbf/cycle (J/cycle)" DATA "QSHTL" DATA "shuttle heat loss, ft-lbf/cycle (J/cycle)" DATA "REFF1" DATA "like REFF2 a measure of regenerator effectiveness" DATA "REFF2" DATA "like REFF1 a measure of regenerator effectiveness" DATA "TG" DATA "gas temperature at interfase of two adjacent control volumes
degrees R (kelvin)" DATA "TGA" DATA "working bulk gas temperature of a control volume, degrees R (kelvin)" DATA "TGACYC" DATA "cycle-time average temperature of TGA, degrees R (kelvin)" DATA "TGCMPA" DATA "cycle-time average temperature in compression space, degrees R (kelvin)" DATA "TGCYC" DATA "cycle time average temperature of TG, degree R (kelvin)" DATA "TGEXPA" DATA "cycle time average temperature in expansion space, degree R (kelvin)" DATA "TM" DATA "metal temperature of a metal control volume, degree R (kelvin)" DATA "TMCYC" DATA "cycke-time average temperarure of TM, degree R (kelvin)" DATA "UTOTAL" DATA "internal energy content of working-space gas, ft-lbf (J)" DATA "WDISP" DATA "net work done on displacer piston by working-space gas, ft-lbf (J)" DATA "WDISPN" DATA "work done on gas by displacer, ft-lbf/cycle (J/cycle)" DATA "WDISPP" DATA "work done on piston by gas, ft-lbf/cycle (J/cycle)" DATA "WPISTP" DATA "work done on power piston by gas, ft-lbf/cycle (J/cycle" DATA "WPIST" DATA "work done on power piston by working-space gas ft-lbf (J)" DATA "WPISTN" DATA "work done on gas by piston, ft-lbf/cycle ft-lbf (J/cycle)" DATA "WRKCMP" DATA "work done by the gas in compression space over one cycle, ft-lbf (J)" DATA "WRKEXP" DATA "worl done by the gas in compression space over one cycle, ft-lbf (J)" DATA "WRKLP" DATA "power-piston work lost due to pressure drop from regenerator
to compression space, ft-lbf/cycle J/cycle" DATA "WRKLT" DATA "work lost due to pressure drop over one cycle, ft-lbf (J)" DATA "WORKTOT" DATA "net work done by gas in working space over one cycle,
WRKEXP + WRKCMP, ft-lbf (J)"