      SUBROUTINE READIN(MODE)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
C     MODE=1 COMPUTE A MODEL
C     MODE=2 READ A PREVIOUSLY CALCULATED MODEL FOR SOME APPLICATION
C     MODE=20 SAME AS 2 BUT ON ENCOUNTERING END RETURN WITH NRHOX=0
      COMMON /ABROSS/ABROSS(kw),TAUROS(kw)
      COMMON /CONV/DLTDLP(kw),HEATCP(kw),DLRDLT(kw),VELSND(kw),
     1             GRDADB(kw),HSCALE(kw),FLXCNV(kw),VCONV(kw),MIXLTH,
     2             OVERWT,FLXCNV0(kw),FLXCNV1(kw),IFCONV
      REAL*8 MIXLTH
      COMMON /DEPART/BHYD(kw,6),BMIN(kw),NLTEON
      COMMON /ELEM/ABUND(99),ATMASS(99),ELEM(99)
      COMMON /FILENAME/FILENAME
      CHARACTER*60 FILENAME,kapname
      COMMON /FLUX/ FLUX,FLXERR(kw),FLXDRV(kw),FLXRAD(kw)
      COMMON /FREEd/WORD(6),NUMCOL,LETCOL,LAST,MORE,IFFAIL,MAXPOW
      REAL*4 WORD
      COMMON /FRESET/FRESET(1563),RCOSET(1563),NULO,NUHI,NUMNU
      COMMON /HEIGHT/HEIGHT(kw)
      COMMON /IF/IFCORR,IFPRES,IFSURF,IFSCAT,TAUSCAT,IFMOL
      COMMON /IFOP/IFOP(20)
      COMMON /ITER/ ITER,IFPRNT(15),IFPNCH(15),NUMITS
      COMMON /JUNK/TITLE(74),FREQID(6),WLTE,XSCALE
      INTEGER*4 ITITLE(148)
      EQUIVALENCE (ITITLE(1),TITLE(1))
      REAL*4 FREQID4(12)
      EQUIVALENCE (FREQID(1),FREQID4(1))
      COMMON /MUS/ANGLE(20),SURFI(20),NMU
      COMMON /PTOTAL/PTOTAL(kw)
      COMMON /PZERO/PZERO,PCON,PRADK0,PTURB0,KNU(kw),PRADK(kw),RADEN(kw)
      REAL*8 KNU
      COMMON /RAD/ACCRAD(kw),PRAD(kw)
      COMMON /RHOX/RHOX(kw),NRHOX
      COMMON /STATE/P(kw),XNE(kw),XNATOM(kw),RHO(kw)
      COMMON /STEPLG/STEPLG,TAU1LG,KRHOX
      COMMON/TABTP/TABT(56),TABP(38),TABKAP(56,38,5),PKAP(38),
     1TABV(5),PKAPV(5,38),IFTABK,KTAB(56,38,5),NT,NP
      COMMON /TEFF/TEFF,GRAV,GLOG
      COMMON /TEMP/T(kw),TKEV(kw),TK(kw),HKT(kw),TLOG(kw),ITEMP
      COMMON /TURBPR/VTURB(kw),PTURB(kw),TRBFDG,TRBCON,TRBPOW,TRBSND,
     1               IFTURB
      COMMON /WAVEY/WBEGIN,DELTAW,IFWAVE
      COMMON /XABUND/XABUND(99),WTMOLE
      INTEGER CARD(80),CARD81,CCCC(81)
      EQUIVALENCE (CARD(1),CCCC(1)),(CARD81,CCCC(81))
	common/cards/card
      DIMENSION RHOXA(kw),DUM1(kw),DUM2(kw),DUM3(kw),DUM4(kw),DUM5(kw)
      DIMENSION DUM6(kw),DUM7(kw),DUM8(kw),TAUSTD(kw),TSMOOTH(kw)
C      EQUIVALENCE (DUM1(1),DLTDLP(1)),(DUM2(1),HEATCP(1))
C      EQUIVALENCE (DUM3(1),DLRDLT(1)),(DUM4(1),VELSND(1))
C      EQUIVALENCE (DUM5(1),GRDADB(1)),(DUM6(1),HSCALE(1))
C      EQUIVALENCE (DUM7(1),FLXCNV(1)),(RHOXA(1),VCONV(1),DUM8(1))
C      EQUIVALENCE (TAUSTD(1),XNATOM(1))
      DIMENSION IFOP1(20)
C     H1,H2PLUS,HMINUS,HRAY,HE1,HE2,HEMINUS,HERAY,COOL,LUKE,
C     HOT,ELECTRON,H2RAY,HLINES,LINES,LINESCAT,XLINES,XLSCAT,XCONT,XSCAT
      DATA IFOP1/324,609929997,579591588,429928,11165,11166,564793810,
     1 15271257,173061,636997,11527,369467847,16486929,577716835,
     2 22965179,849711626,1687220147,1687711471,45152896,45946435/
      DATA WWLTE,WWNLTE/4HLTE ,4HNLTE/
      DATA BLANK/1H /
      EXP10(X)=EXP(X*2.30258509299405E0)
C      CARD81=BLANK
      CARD81=(1H )
      LAST=81
      MAXPOW=38
   98 MORE=0
      LETCOL=1
   99 READ(5,1) CARD
    1 FORMAT(80A1)
C     (M)ACHINE (I)NDEPENDENT (A)LPHAMERIC (C)ODE
C     BASE 37    A=1, Z=26, 0=27, 9=36
      MIAC=IWORDF(CARD)
      NUMCOL=LETCOL
C     TEFF
    3 IF(MIAC.EQ.1020133)GO TO 100
C     GRAVITY
      IF(MIAC.EQ.519223721)GO TO 200
C     OPACITY
      IF(MIAC.EQ.1070201044)GO TO 300
C     KAPPA
      IF(MIAC.EQ.20688921)GO TO 400
C     ITERATIONS
      IF(MIAC.EQ.661856797)GO TO 500
C     MOLECULES
      IF(MIAC.EQ.930198669)GO TO 600
C     CALCULATE
      IF(MIAC.EQ.210518764)GO TO 700
C     ABUNDANCE
      IF(MIAC.EQ.74175307)GO TO 800
C     PRINT
      IF(MIAC.EQ.30911189)GO TO 900
C     PUNCH
      IF(MIAC.EQ.31069574)GO TO 1000
C     READ
      IF(MIAC.EQ.918640)GO TO 1100
C     LTE
      IF(MIAC.EQ.17173)GO TO 1200
C     NLTE
      IF(MIAC.EQ.726315)GO TO 1300
C     BEGIN
      IF(MIAC.EQ.4011517)GO TO 1500
C     SCATTERING
      IF(MIAC.EQ.1323236444)GO TO 1600
C     END
      IF(MIAC.EQ.7367)GO TO 1700
C     TITLE
      IF(MIAC.EQ.37966926)GO TO 1800
C     CONVECTION
      IF(MIAC.EQ.236883734)GO TO 1900
C     TURBULENCE
      IF(MIAC.EQ.1427151802)GO TO 2000
C     CHANGE RHOX
      IF(MIAC.EQ.223095242)GO TO 2100
C     FREQUENCIES
      IF(MIAC.EQ.450075960)GO TO 2200
C     SURFACE
      IF(MIAC.EQ.1357812572)GO TO 2300
C     PRESSURE
      IF(MIAC.EQ.1143518210)GO TO 2400
C     CORRECTION
      IF(MIAC.EQ.237080870)GO TO 2500
C     WAVELENGTH
      IF(MIAC.EQ.1597906832)GO TO 2600
C     SCALE MODEL
      IF(MIAC.EQ.35762836)GO TO 2700
C     CALL
      IF(MIAC.EQ.153784)GO TO 2800
C     VTURB
      IF(MIAC.EQ.42274019)GO TO 2900
C     SMOOTH
      IF(MIAC.EQ.1342680354)GO TO 3000
C
 9000 WRITE(6,2) CARD
    2 FORMAT(21H I DO NOT UNDERSTAND 80A1)
      stop
   97 LETCOL=MAX0(LETCOL,NUMCOL)
      MORE=1
      MIAC=IWORDF(CARD)
      IF(IFFAIL.EQ.1)GO TO 98
      MORE=0
      GO TO 3
C************
  100 TEFF=FREEFF(CARD)
      FLUX=5.6697E-5/12.5664*TEFF**4
      GO TO 97
C************
  200 GRAV=FREEFF(CARD)
      IF(GRAV.LT.10.)GRAV=EXP10(GRAV)
      GLOG=DLOG10(GRAV)
      GO TO 97
C************
  300 MIAC=IWORDF(CARD)
C     ON
      IF(MIAC.EQ.569)GO TO 380
C     OFF
      IF(MIAC.EQ.20763)GO TO 390
C     IFOP
      IF(MIAC.EQ.464662)GO TO 370
      GO TO 9000
  370 NUMCOL=LETCOL
      DO 371 I=1,20
  371 IFOP(I)=FREEFF(CARD)
      GO TO 98
C     ON
  380 ISWCH=1
      GO TO 391
C     OFF
  390 ISWCH=0
  391 MORE=1
  395 MIAC=IWORDF(CARD)
      IF(IFFAIL.EQ.1)GO TO 97
      DO 392 I=1,20
      II=I
      IF(MIAC.EQ.IFOP1(I))GO TO 393
  392 CONTINUE
      GO TO 9000
  393 IFOP(II)=ISWCH
      GO TO 395
C************
  400 GO TO 9000
C************
  500 NUMITS=FREEFF(CARD)
      DO 501 I=1,15
  501 IFPNCH(I)=0
      IFPNCH(NUMITS)=1
      GO TO 97
C************
  600 MIAC=IWORDF(CARD)
C     ON
      IF(MIAC.EQ.569)GO TO 610
C     OFF
      IF(MIAC.EQ.20763)GO TO 620
      GO TO 9000
  610 IFMOL=1
      GO TO 97
  620 IFMOL=0
      GO TO 97
C************
  700 NRHOX=FREEFF(CARD)
      TAU1LG=FREEFF(CARD)
      STEPLG=FREEFF(CARD)
      DO 701 J=1,NRHOX
      TAUROS(J)=EXP10(TAU1LG+DFLOAT(J-1)*STEPLG)
  701 T(J)=TEFF*(.75*(.710+TAUROS(J)-.1331*EXP(-3.4488*TAUROS(J))))**.25
  702 DO 703 J=1,NRHOX
      XNE(J)=0.
      PRADK(J)=2.521E-15*DMAX1(T(J)**4,TEFF**4/2.)
      VTURB(J)=0.
  703 PTURB(J)=0.
      PRADK0=PRADK(1)
      PCON=0.
      PTURB0=0.
      PZERO=PRADK0
      DO 704 J=1,NRHOX
  704 PRAD(J)=PRADK(J)-PRADK0
      CALL TTAUP(T,TAUROS,ABROSS,PTOTAL,P,PRAD,PTURB,
     1VTURB,GRAV,NRHOX)
      DO 705 J=1,NRHOX
      RHOX(J)=PTOTAL(J)/GRAV
  705 PTOTAL(J)=PTOTAL(J)+PZERO
      GO TO 97
C************
  800 MIAC=IWORDF(CARD)
C     SCALE
      IF(MIAC.EQ.35762836)GO TO 810
C     CHANGE  (SOLAR ABUNDANCE BEFORE SCALING)
      IF(MIAC.EQ.223095242)GO TO 820
C     ABSOLUTE (FINAL ABUNDANCE AFTER SCALING)
      IF(MIAC.EQ.74075686)GO TO 830
      GO TO 9000
  810 NUMCOL=LETCOL
      XSCALE=FREEFF(CARD)
      GO TO 97
  820 MORE=1
  821 IZ=FREEFF(CARD)
      IF(IFFAIL.EQ.1)GO TO 98
      ABUND(IZ)=FREEFF(CARD)
      GO TO 821
  830 MORE=1
  831 IZ=FREEFF(CARD)
      IF(IFFAIL.EQ.1)GO TO 98
      ABUND(IZ)=FREEFF(CARD)
      IF(IZ.GT.2)ABUND(IZ)=ABUND(IZ)/XSCALE
      GO TO 831
C************
  900 DO 901 I=1,NUMITS
  901 IFPRNT(I)=FREEFF(CARD)
      GO TO 97
C************
 1000 DO 1001 I=1,NUMITS
 1001 IFPNCH(I)=FREEFF(CARD)
      GO TO 97
C************
 1100 MIAC=IWORDF(CARD)
      NUMCOL=LETCOL
C     FREQUENCIES
      IF(MIAC.EQ.450075960)GO TO 1110
C     DEPARTURE COEFFICIENTS
      IF(MIAC.EQ.287559136)GO TO 1120
C     STARTING T-TAU
      IF(MIAC.EQ.1355094447)GO TO 1130
C     DECK
      IF(MIAC.EQ.209579)GO TO 1140
C     DECK6
      IF(MIAC.EQ.7754456)GO TO 1140
C     KAPPA
      IF(MIAC.EQ.20688921)GO TO 1150
C     MOLECULES
      IF(MIAC.EQ.930198669)GO TO 1160
C     PUNCH
      IF(MIAC.EQ.31069574)GO TO 1170
      GO TO 9000
C     FREQUENCIES
 1110 NUM=FREEFF(CARD)
      NULO=FREEFF(CARD)
      NUHI=FREEFF(CARD)
      NUMNU=NUM
      LETCOL=NUMCOL
      NDUMMY=IWORDF(CARD)
      DO 1111 I=1,6
 1111 FREQID4(2*I-1)=WORD(I)
      NUMCOL=LETCOL
      DO 1112 I=1,NUMNU
      NU=FREEFR(CARD)
      FRESET(NU)=FREEFF(CARD)
C     PROVISION FOR READING WAVELENGTHS
      IF(FRESET(NU).LT.1.E7)FRESET(NU)=2.997925E17/FRESET(NU)
C     PROVISION FOR READING WAVENUMBERS SCALED BY 1.E25
      IF(FRESET(NU).GT.1.E20)FRESET(NU)=FRESET(NU)*(2.997925E10/1.E25)
 1112 RCOSET(NU)=FREEFF(CARD)
      GO TO 98
C     DEPARTURE COEFFICIENTS
 1120 NRHOX=FREEFF(CARD)
      DO 1122 J=1,NRHOX
      NUMCOL=1
      READ(5,1)CARD
      DUMMY=FREEFF(CARD)
      DO 1121 I=1,6
 1121 BHYD(J,I)=FREEFF(CARD)
 1122 BMIN(J)=FREEFF(CARD)
      WLTE=WWNLTE
      NLTEON=1
      GO TO 98
 1130 NRHOX=FREEFF(CARD)
      DO 1131 J=1,NRHOX
      NUMCOL=1
      READ(5,1)CARD
      TAUROS(J)=FREEFF(CARD)
 1131 T(J)=FREEFF(CARD)
      IF(TAUROS(1).GT.0.)GO TO 702
      DO 1132 J=1,NRHOX
 1132 TAUROS(J)=EXP10(TAUROS(J))
      GO TO 702
 1140 NRHOX=FREEFF(CARD)
      DO 1141 J=1,NRHOX
      NUMCOL=1
      READ(5,1)CARD
      RHOX(J)=FREEFF(CARD)
      T(J)=FREEFF(CARD)
      MORE=1
      P(J)=FREEFF(CARD)
      XNE(J)=FREEFF(CARD)
      ABROSS(J)=FREEFF(CARD)
      PRAD(J)=FREEFF(CARD)
      VTURB(J)=FREEFF(CARD)
 1141 MORE=0
      IF(RHOX(1).GE.0.)GO TO 1143
      DO 1142 J=1,NRHOX
 1142 RHOX(J)=EXP10(RHOX(J))
 1143 PRADK0=0.
      PTURB0=PTURB(1)
      PCON=0.
      PZER0=PCON+PRADK0+PTURB0
      CALL INTEG(RHOX,ABROSS,TAUROS,NRHOX,ABROSS(1)*RHOX(1))
      IF(MIAC.NE.7754456)GO TO 98
      READ(5,1)CARD
      NUMCOL=1
      PRADK0=FREEFF(CARD)
      DO 1144 J=1,NRHOX
 1144 ACCRAD(J)=PRAD(J)
      CALL INTEG(RHOX,ACCRAD,PRAD,NRHOX,ACCRAD(1)*RHOX(1))
      DO 1145 J=1,NRHOX
 1145 PRADK(J)=PRAD(J)+PRADK0
      GO TO 98

1150	read(5,'((a))')kapname
        OPEN(1,file=kapname,status='old')
      READ(1,1)CARD
c      WRITE(6,1151)CARD
 1151 FORMAT(80A1)
      READ(1,1)
      DO 1154 IT=1,56
      DO 1154 IP=1,38
      READ(1,1152)(TABKAP(IT,IP,IV),IV=1,5)
 1152 FORMAT(5X,5X,5F7.3)
 1154 CONTINUE
      IFTABK=1
      CLOSE(UNIT=1)
      GO TO 98
 1160 CALL READMOL
      GO TO 98
C     READS A FILE THAT HAS BEEN PUNCHED PREVIOUSLY EXCEPT FOR THE BEGIN CARD
1170	read(5,'((a))')filename
	open(3,file=filename,status='old')
      READ(3,1171) TEFF,GLOG,WLTE,TITLE,IFOP,A,MIXLTH,B,TRBFDG,
     1TRBPOW,TRBSND,TRBCON,XSCALE,(IZ,ABUND(IZ),IZ=1,99)
 1171 FORMAT(5HTEFF ,F7.0,9H  GRAVITY,F8.5,1X,A4/6HTITLE ,74A1/
     1 13H OPACITY IFOP,20I2/
     2 12H CONVECTION ,A3,F6.2,12H TURBULENCE, A3,4F6.2/
     3 16HABUNDANCE SCALE ,F9.5,17H ABUNDANCE CHANGE,2(I2,F8.5)/
     4(17H ABUNDANCE CHANGE,6(I3,F7.2)))
      GRAV=10.**GLOG
      FLUX=5.6697E-5/12.5664*TEFF**4
      IF(WLTE.EQ.4HLTE )NLTEON=0
      IF(WLTE.EQ.4HNLTE)NLTEON=1
      IF(A.EQ.3HOFF)IFCONV=0
      IF(A.EQ.3HON )IFCONV=1
      IF(B.EQ.3HOFF)IFTURB=0
      IF(B.EQ.3HON )IFTURB=1
      READ(3,1172)NRHOX,(RHOX(J),T(J),P(J),XNE(J),ABROSS(J),ACCRAD(J),
     1VTURB(J),J=1,NRHOX)
 1172 FORMAT(10HREAD DECK6,I3,33H RHOX,T,P,XNE,ABROSS,ACCRAD,VTURB/
     1(1PE15.8,0PF9.1,1P5E10.3))
      READ(3,1173)PRADK0
 1173 FORMAT(5HPRADK,1PE11.4)
      CALL INTEG(RHOX,ACCRAD,PRAD,NRHOX,ACCRAD(1)*RHOX(1))
      DO 1177 J=1,NRHOX
 1177 PRADK(J)=PRAD(J)+PRADK0
      IF(NLTEON.EQ.1)THEN
      READ(3,1174)NRHOX,(RHOX(J),(BHYD(J,I),I=1,6),BMIN(J),J=1,NRHOX)
 1174 FORMAT(27HREAD DEPARTURE COEFFICIENTSI3,21H RHOX  BHYD 1-6  BMIN/
     1(1PE11.4,0P7F9.4))
      ENDIF
      READ(3,1175)ITERPUN
 1175 FORMAT(5HBEGIN,20X,10HITERATION I3,10H COMPLETED )
      CLOSE(UNIT=3)
      GO TO 98
C************
 1200 NLTEON=0
      WLTE=WWLTE
      DO 1202 J=1,kw
      DO 1201 I=1,6
 1201 BHYD(J,I)=1.
 1202 BMIN(J)=1.
      GO TO 97
C************
 1300 NLTEON=1
      WLTE=WWNLTE
      GO TO 97
C************
 1500 IF(MODE.NE.1)GO TO 1510
      IF(NUMITS.EQ.0)stop
      IF(NRHOX.EQ.0)stop
      IF(NUMNU.EQ.0)stop
      IF(TEFF.EQ.0.)stop
      IF(GRAV.EQ.0.)stop
 1510 CONTINUE
      IF(ABUND(1).LT.0.)ABUND(1)=EXP10(ABUND(1))
      IF(ABUND(2).LT.0.)ABUND(2)=EXP10(ABUND(2))
      DO 1511 IZ=3,99
      IF(ABUND(IZ).GT.0.)ABUND(IZ)=DLOG10(ABUND(IZ))
 1511 CONTINUE
      WRITE(6,1512)TEFF,GLOG,WLTE,TITLE,XSCALE,
     1(ELEM(IZ),ABUND(IZ),IZ=1,99)
 1512 FORMAT(1H1/////5H TEFF,F7.0,8H   LOG G,F8.4,3X,A4/
     17H0TITLE ,74A1/7H0XSCALE,F10.6,2(3X,A2,F8.5)/(10(1X,A2,F6.2)))
      DO 1513 IZ=3,99
 1513 XABUND(IZ)=EXP10(ABUND(IZ))*XSCALE
      XABUND(1)=ABUND(1)
      XABUND(2)=ABUND(2)
      WTMOLE=0.
      DO 1514 IZ=1,99
 1514 WTMOLE=WTMOLE+XABUND(IZ)*ATMASS(IZ)
      DO 1516 J=1,NRHOX
      TK(J)=1.38054E-16*T(J)
      HKT(J)=6.6256E-27/TK(J)
      TKEV(J)=8.6171E-5*T(J)
      TLOG(J)=DLOG(T(J))
      XNATOM(J)=P(J)/TK(J)-XNE(J)
      RHO(J)=XNATOM(J)*WTMOLE*1.660E-24
      IF(IFTURB.GT.0)PTURB(J)=.5*RHO(J)*VTURB(J)**2
 1516 CONTINUE
      WRITE(6,1517)IFOP
 1517 FORMAT(3H0H1I2,7H H2PLUSI2,7H HMINUSI2,5H HRAYI2,4H HE1I2,
     1 4H HE2I2,8H HEMINUSI2,6H HERAYI2,5H COOLI2,5H LUKEI2/
     2 4H HOTI2,9H ELECTRONI2,6H H2RAYI2,7H HLINESI2,6H LINESI2,
     3 9H LINESCATI2,7H XLINESI2,7H XLSCATI2,6H XCONTI2,6H XSCATI2)
      WRITE(6,1518)IFCORR,IFPRES,IFSURF,IFSCAT,TAUSCAT,IFCONV,MIXLTH,
     1IFMOL,IFTURB,TRBFDG,TRBPOW,TRBSND,TRBCON
 1518 FORMAT(7H0IFCORR,I2,8H  IFPRES,I2,8H  IFSURF,I2,8H  IFSCAT,I2,
     1 9H  TAUSCAT,1PE10.2,
     2 8H  IFCONV,I2,8H  MIXLTH,0PF6.2,7H  IFMOL,I2/7H IFTURB,I2,
     3 8H  TRBFDG,F6.2,8H  TRBPOW,F6.2,8H  TRBSND,F6.2,8H  TRBCON,F6.2)
      IF(MODE.NE.1)GO TO 1575
      WRITE(6,1521)NUMITS,IFPRNT,IFPNCH
 1521 FORMAT(7H NUMITS,I3,8H  IFPRNT,15I2,8H  IFPNCH,15I2)
      IF(IFWAVE.EQ.0)GO TO 1560
      WRITE(6,1536)WBEGIN,DELTAW,NUMNU
 1536 FORMAT(7H0WBEGIN,F11.4,9H   DELTAW,F7.4,8H   NUMNU,I5)
      GO TO 1575
 1560 WRITE(6,1561)FREQID,NUMNU,NULO,NUHI
 1561 FORMAT(8H0FREQID ,6A1,8H   NUMNU,I4,7H   NULO,I4,7H   NUHI,I4)
      NN=(NUMNU+3)/4
      NNN=NUMNU-NN*3
      IF(NNN.LT.1)NN=1
      IF(NNN.LT.1)NNN=1
      WRITE(6,1563)((NU,FRESET(NU),RCOSET(NU),NU=I,NUMNU,NN),I=1,NNN)
 1563 FORMAT((4(I6,1P2E13.6)))
      IF(NN.EQ.NNN)GO TO 1575
      NNN=NNN+1
      WRITE(6,1564)((NU,FRESET(NU),RCOSET(NU),NU=I,NUMNU,NN),I=NNN,NN)
 1564 FORMAT((3(I6,1P2E13.6)))
      SUM=0.
      DO 9500 NU=NULO,NUHI
 9500 SUM=SUM+RCOSET(NU)
c      PRINT 9501,SUM
 9501 FORMAT(1PE13.6)
C
 1575 CONTINUE
      WRITE(6,1576)(J,RHOX(J),T(J),P(J),XNE(J),ABROSS(J),PRAD(J),
     1VTURB(J),(BHYD(J,I),I=1,6),BMIN(J),J=1,NRHOX)
 1576 FORMAT(1H1/////9X,4HRHOX,9X,1HT,8X,1HP,8X,3HXNE,6X,6HABROSS,
     1 5X,4HPRAD,6X,5HVTURB,23X,4HBHYD,25X,4HBMIN/
     2(I4,1PE13.6,0PF9.1,1P5E10.3,0P7F8.4))
C
      RETURN
C************
 1600 MIAC=IWORDF(CARD)
C     ON
      IF(MIAC.EQ.569)GO TO 1610
C     OFF
      IF(MIAC.EQ.20763)GO TO 1620
C     LINES
      IF(MIAC.EQ.22965179)GO TO 1630
      GO TO 9000
 1610 IFSCAT=1
      GO TO 97
 1620 IFSCAT=0
      GO TO 97
 1630 IFSCAT=1
      NUMCOL=LETCOL
      TAUSCAT=FREEFF(CARD)
      GO TO 97
C************
 1700 IF(MODE.NE.20)stop
      NRHOX=0
      RETURN
C************
 1800 DO 1801 I=1,74
      III=CARD(I+6)
      ITITLE(2*I)=(1H )
 1801 ITITLE(2*I-1)=III
C 1801 TITLE(I)=CARD(I+6)
      GO TO 98
C************
 1900 MIAC=IWORDF(CARD)
C     ON
      IF(MIAC.EQ.569)GO TO 1910
C     OFF
      IF(MIAC.EQ.20763)GO TO 1920
C     OVER
      IF(MIAC.EQ.790116)GO TO 1930
      GO TO 9000
 1910 IFCONV=1
      NUMCOL=LETCOL
      MIXLTH=FREEFF(CARD)
      GO TO 97
 1920 IFCONV=0
      MIXLTH=1.
      DO 1921 J=1,NRHOX
      DLTDLP(J)=0.
      HEATCP(J)=0.
      DLRDLT(J)=0.
      VELSND(J)=0.
      GRDADB(J)=0.
      HSCALE(J)=0.
      FLXCNV(J)=0.
 1921 VCONV(J)=0.
      GO TO 97
 1930 IFCONV=1
      NUMCOL=LETCOL
      MIXLTH=FREEFF(CARD)
      OVERWT=FREEFF(CARD)
      GO TO 97
C************
 2000 MIAC=IWORDF(CARD)
C     ON
      IF(MIAC.EQ.569)GO TO 2010
C     OFF
      IF(MIAC.EQ.20763)GO TO 2020
      GO TO 9000
 2010 IFTURB=1
      NUMCOL=LETCOL
      TRBFDG=FREEFF(CARD)
      TRBPOW=FREEFF(CARD)
      TRBSND=FREEFF(CARD)
      TRBCON=FREEFF(CARD)
      GO TO 97
 2020 IFTURB=0
      TRBFDG=0.
      TRBPOW=0.
      TRBSND=0.
      TRBCON=0.
      GO TO 97
C************
 2100 NNEW=FREEFF(CARD)
      DO 2101 J=1,NNEW
 2101 RHOXA(J)=FREEFR(CARD)
      IDUM=MAP1(RHOX,T,NRHOX,RHOXA,DUM1,NNEW)
      IDUM=MAP1(RHOX,P,NRHOX,RHOXA,DUM2,NNEW)
      IDUM=MAP1(RHOX,XNE,NRHOX,RHOXA,DUM3,NNEW)
      IDUM=MAP1(RHOX,ABROSS,NRHOX,RHOXA,DUM4,NNEW)
      IDUM=MAP1(RHOX,VTURB,NRHOX,RHOXA,DUM5,NNEW)
      IDUM=MAP1(RHOX,PRAD,NRHOX,RHOXA,DUM6,NNEW)
      IDUM=MAP1(RHOX,BMIN,NRHOX,RHOXA,DUM7,NNEW)
      DO 2102 J=1,NNEW
      T(J)=DUM1(J)
      P(J)=DUM2(J)
      XNE(J)=DUM3(J)
      ABROSS(J)=DUM4(J)
      VTURB(J)=DUM5(J)
      PRAD(J)=DUM6(J)
      PRADK(J)=PRAD(J)+PRADK0
 2102 BMIN(J)=DUM7(J)
      DO 2105 I=1,6
      IDUM=MAP1(RHOX,BHYD(1,I),NRHOX,RHOXA,DUM1,NNEW)
      DO 2104 J=1,NNEW
 2104 BHYD(J,I)=DUM1(J)
 2105 CONTINUE
      NRHOX=NNEW
      DO 2106 J=1,NRHOX
 2106 RHOX(J)=RHOXA(J)
      GO TO 97
C************
 2200 NUM=FREEFF(CARD)
      NULO=FREEFF(CARD)
      NUHI=FREEFF(CARD)
      LETCOL=NUMCOL
      MIAC=IWORDF(CARD)
      NUMNU=NUM
      DO 2201 I=1,6
 2201 FREQID4(2*I-1)=WORD(I)
C     BIG
      IF(MIAC.EQ.3078)GO TO 2210
C     LITTLE
      IF(MIAC.EQ.850035822)GO TO 2220
      GO TO 9000
C     FREQUENCY SET BIG
C     NUMNU=337 IS CURRENT, NUMNU=331 IS OLD AND WILL WORK
 2210 DO 2211 NU=1,NUMNU
 2211 FRESET(NU)=2.997925E17/FRESET(NU+5)
      DO 2212 NU=2,NUMNU-1
 2212 RCOSET(NU)=(FRESET(NU-1)-FRESET(NU+1))/2.
      RCOSET(1)=(2.997925E17/8.97666-FRESET(2))/2.
      RCOSET(NUMNU)=FRESET(NUMNU-1)/2.
      GO TO 98
C     FREQUENCY SET LITTLE
C     NUMNU=1221 IS CURRENT, NUMNU=1215 IS OLD AND WILL WORK
 2220 DO 2221 NU=1,NUMNU
 2221 FRESET(NU)=2.997925E17/FRESET(NU+342)
      DO 2222 NU=2,NUMNU-1
 2222 RCOSET(NU)=(FRESET(NU-1)-FRESET(NU+1))/2.
      RCOSET(1)=(2.997925E17/8.97666-FRESET(2))/2.
      RCOSET(NUMNU)=FRESET(NUMNU-1)/2.
      GO TO 98
C************
 2300 MIAC=IWORDF(CARD)
C     INTENSITY
      IF(MIAC.EQ.651354309)GO TO 2310
C     FLUX
      IF(MIAC.EQ.321147)GO TO 2320
C     OFF
      IF(MIAC.EQ.20763)GO TO 2330
      GO TO 9000
 2310 NMU=FREEFF(CARD)
      DO 2311 MU=1,NMU
 2311 ANGLE(MU)=FREEFR(CARD)
      IFSURF=2
      GO TO 97
 2320 IFSURF=1
      GO TO 97
 2330 IFSURF=0
      GO TO 97
C************
 2400 MIAC=IWORDF(CARD)
C     ON
      IF(MIAC.EQ.569)GO TO 2410
C     OFF
      IF(MIAC.EQ.20763)GO TO 2420
      GO TO 9000
 2410 IFPRES=1
      GO TO 97
 2420 IFPRES=0
      GO TO 97
C************
 2500 MIAC=IWORDF(CARD)
C     ON
      IF(MIAC.EQ.569)GO TO 2510
C     OFF
      IF(MIAC.EQ.20763)GO TO 2520
      GO TO 9000
 2510 IFCORR=1
      GO TO 97
 2520 IFCORR=0
      DO 2521 J=1,kw
      FLXERR(J)=0.
 2521 FLXDRV(J)=0.
      GO TO 97
C************
 2600 WBEGIN=FREEFF(CARD)
      DELTAW=FREEFF(CARD)
      WEND=FREEFF(CARD)
      IFWAVE=1
      NULO=1
      NUHI= INT(ABS(WEND-WBEGIN)/ABS(DELTAW)+.5)+1
      NUMNU=NUHI
      IF(WBEGIN.LT.1.E7)GO TO 97
      IF(WBEGIN.GT.1.E20)GO TO 2610
C     FREQUENCY STEPS
      WBEGIN=2.997925E17/WBEGIN
      DELTAW=2.997925E17/DELTAW
      WEND=2.997925E17/WEND
      GO TO 97
C     WAVENUMBER STEPS SCALE BY 1.E25
 2610 WBEGIN=1.E7/(WBEGIN/1.E25)
      DELTAW=1.E7/(DELTAW/1.E25)
      WEND=1.E7/(WEND/1.E25)
      GO TO 97
C************
C     SCALING MODELS OR CHANGING RHOX SPACING TO BE UNIFORM IN TAUROS
 2700 KRHOX=FREEFF(CARD)
      TAU1LG=FREEFF(CARD)
      STEPLG=FREEFF(CARD)
      MORE=1
      TEFF1=FREEFF(CARD)
      GNEW=FREEFF(CARD)
      IF(GNEW.LT.10.)GNEW=EXP10(GNEW)
      MORE=0
      DO 2701 J=1,KRHOX
 2701 TAUSTD(J)=EXP10(TAU1LG+DFLOAT(J-1)*STEPLG)
      CALL INTEG(RHOX,ABROSS,TAUROS,NRHOX,ABROSS(1)*RHOX(1))
      TAUROS(1)=MIN(TAUROS(1),TAUSTD(1))
      IDUM=MAP1(TAUROS,RHOX,NRHOX,TAUSTD,DUM1,KRHOX)
      IDUM=MAP1(TAUROS,T,NRHOX,TAUSTD,DUM2,KRHOX)
      IDUM=MAP1(TAUROS,P,NRHOX,TAUSTD,DUM3,KRHOX)
      IDUM=MAP1(TAUROS,XNE,NRHOX,TAUSTD,DUM4,KRHOX)
      IDUM=MAP1(TAUROS,ABROSS,NRHOX,TAUSTD,DUM5,KRHOX)
      IDUM=MAP1(TAUROS,PRAD,NRHOX,TAUSTD,DUM6,KRHOX)
      IDUM=MAP1(TAUROS,VTURB,NRHOX,TAUSTD,DUM7,KRHOX)
      IDUM=MAP1(TAUROS,BMIN,NRHOX,TAUSTD,DUM8,KRHOX)
      DO 2702 J=1,KRHOX
      RHOX(J)=DUM1(J)
      T(J)=DUM2(J)
      P(J)=DUM3(J)
      XNE(J)=DUM4(J)
      ABROSS(J)=DUM5(J)
      PRAD(J)=DUM6(J)
      PRADK(J)=PRAD(J)+PRADK0
      VTURB(J)=DUM7(J)
 2702 BMIN(J)=DUM8(J)
      DO 2704 I=1,6
      IDUM=MAP1(TAUROS,BHYD(1,I),NRHOX,TAUSTD,DUM1,KRHOX)
      DO 2703 J=1,KRHOX
 2703 BHYD(J,I)=DUM1(J)
 2704 CONTINUE
      NRHOX=KRHOX
      IF(TEFF1.EQ.0.)GO TO 97
      IF(TEFF1.EQ.TEFF.AND.GNEW.EQ.GRAV)GO TO 97
      IF(TEFF1.LT.TEFF+1.AND.TEFF1.GT.TEFF-1..AND.
     1GNEW.LT.GRAV*1.001.AND.GNEW.GT.GRAV*.999)GO TO 97
      DO 2710 J=1,NRHOX
      TAUROS(J)=TAUSTD(J)
      T(J)=T(J)*TEFF1/TEFF
      PTURB(J)=0.
      PRADK(J)=PRADK(J)*(TEFF1/TEFF)**4
 2710 PRAD(J)=PRAD(J)*(TEFF1/TEFF)**4
      PRADK0=PRADK0*(TEFF1/TEFF)**4
      PZERO=PCON+PRADK0+PTURB0
      TEFF=TEFF1
      FLUX=5.6697E-5/12.5664*TEFF**4
      GRAV=GNEW
      GLOG=DLOG10(GRAV)
      CALL TTAUP(T,TAUROS,ABROSS,PTOTAL,P,PRAD,PTURB,
     1VTURB,GRAV,NRHOX)
      DO 2711 J=1,NRHOX
      RHOX(J)=PTOTAL(J)/GRAV
 2711 PTOTAL(J)=PTOTAL(J)+PZERO
C2711 RHOX(J)=PTOTAL(J)/GRAV
      GO TO 97
C************
 2800 CALL DUMMYR
      GO TO 97
C************
 2900 VNEW=FREEFF(CARD)
C     VTURB IN CM/S
      DO 2901 J=1,NRHOX
 2901 VTURB(J)=VNEW
      GO TO 97
C************
C     SMOOTH TEMPERTURE DISTRIBUTION
 3000 J1=FREEFF(CARD)
      J2=FREEFF(CARD)
      WTJM1=FREEFF(CARD)
      WTJ=FREEFF(CARD)
      WTJP1=FREEFF(CARD)
      DO 3001 J=J1,J2
 3001 TSMOOTH(J)=WTJM1*T(J-1)+WTJ*T(J)+WTJP1*T(J+1)
      DO 3002 J=J1,J2
 3002 T(J)=TSMOOTH(J)
      GO TO 97
      END
      SUBROUTINE DUMMYR
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
C     DUMMY INPUT ROUTINE FOR MODIFYING READIN
C
C
      RETURN
      END
      FUNCTION FREEFR(CARD)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
      COMMON /FREEd/WORD(6),NUMCOL,LETCOL,LAST,MORE,IFFAIL,MAXPOW
      INTEGER CARD(1),WORD
      MORE=1
      FREEFR=FREEFF(CARD)
      IF(IFFAIL.EQ.0)RETURN
      L=LAST-1
      READ(5,1)(CARD(I),I=1,L)
    1 FORMAT(80A1)
      NUMCOL=1
      FREEFR=FREEFF(CARD)
      RETURN
      END
      FUNCTION FREEFF(CARD)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
      COMMON /FREEd/WORD(6),NUMCOL,LETCOL,LAST,MORE,IFFAIL,MAXPOW
      DIMENSION CARD(1)
      DIMENSION A(10)
      INTEGER CARD,A,C,QPT,QCM,QMI,QE,QPL,QBL,WORD
      DATA A/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/
      DATA QPT/1H./,QCM/1H,/,QMI/1H-/,QE/1HE/,QPL/1H+/,QBL/1H /
      IFFAIL=0
      IF(NUMCOL.GT.LAST)GO TO 3002
      ANSWER=0.
      ASIGN=1.
      ISIGN=1
      NPT=0
      IF0=0
      N=0
      ASSIGN 100 TO NSWCH
      DO 1000 NCOL=NUMCOL,LAST
      C=CARD(NCOL)
      GO TO NSWCH,(100,200,300,400)
  100 IF(C.EQ.QBL)GO TO 104
      DO 101 I=1,10
      IF(C.EQ.A(I))GO TO 102
  101 CONTINUE
      IF(C.EQ.QPT)GO TO 103
      IF(C.EQ.QCM)GO TO 104
      IF(C.EQ.QMI)GO TO 105
  999 ASIGN=1.
      ANSWER=0.
      NPT=0
      IF0=0
      N=0
      ASSIGN 100 TO NSWCH
      GO TO 1000
  102 N=N+1
      ANSWER=10.*ANSWER+DFLOAT(I-1)
      GO TO 1000
  103 ASSIGN 200 TO NSWCH
      GO TO 1000
  104 IF(N.EQ.0)GO TO 999
      FREEFF=ANSWER*ASIGN
      GO TO 998
  105 IF(N.EQ.0)GO TO 106
      GO TO 999
  106 ASIGN=-1.
      GO TO 1000
  200 DO 201 I=1,10
      IF(C.EQ.A(I))GO TO 202
  201 CONTINUE
      IF(C.EQ.QE )GO TO 203
      IF(C.EQ.QMI)GO TO 304
      IF(C.EQ.QPL)GO TO 303
      IF(C.EQ.QBL)GO TO 204
      IF(C.EQ.QCM)GO TO 204
      GO TO 999
  202 N=N+1
      NPT=NPT+1
      ANSWER=10.*ANSWER+DFLOAT(I-1)
      GO TO 1000
  203 ASSIGN 300 TO NSWCH
      GO TO 1000
  204 IF(N.EQ.0)GO TO 999
      FREEFF=ANSWER*ASIGN/10.**NPT
      GO TO 998
  300 DO 301 I=1,10
      IF(C.EQ.A(I))GO TO 302
  301 CONTINUE
      IF(C.EQ.QBL)GO TO 303
      IF(C.EQ.QMI)GO TO 304
      IF(C.EQ.QPL)GO TO 303
      GO TO 999
  302 NPOWER=I-1
      IF0=1
  310 ASSIGN 400 TO NSWCH
      GO TO 1000
  303 NPOWER=0
      GO TO 310
  304 ISIGN=-1
      NPOWER=0
      GO TO 310
  400 DO 401 I=1,10
      IF(C.EQ.A(I))GO TO 402
  401 CONTINUE
      IF(C.EQ.QCM)GO TO 403
      IF(C.EQ.QBL)GO TO 403
      GO TO 999
  402 NPOWER=10*NPOWER+I-1
      IF0=1
      IF(NPOWER.GE.MAXPOW)GO TO 999
      GO TO 1000
  403 IF(IF0.EQ.0)GO TO 999
      FREEFF=ANSWER*ASIGN*10.**(ISIGN*NPOWER-NPT)
      GO TO 998
 1000 CONTINUE
      NUMCOL=LAST+1
 3002 IFFAIL=1
      IF(MORE.GT.0)GO TO 3000
      WRITE(6,3001)(CARD(I),I=1,LAST)
 3001 FORMAT(28H1FREEFF HAS READ OFF THE END/(1X80A1))
      stop
 3000 FREEFF=0.
      RETURN
  998 NUMCOL=NCOL+1
      RETURN
      END
      FUNCTION IWORDF(CARD)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
      COMMON /FREEd/WORD(6),NUMCOL,LETCOL,LAST,MORE,IFFAIL,MAXPOW
      DIMENSION CARD(1)
      DIMENSION A(36)
      EQUIVALENCE (QE,A(5))
      INTEGER CARD,A,C,WORD,QPT,QCM,QEQ,QBL,QE
      DATA A/1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM,1HN,
     11HO,1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,1HZ,1H0,1H1,1H2,1H3,
     21H4,1H5,1H6,1H7,1H8,1H9/
      DATA QPT/1H./,QCM/1H,/,QEQ/1H=/,QBL/1H /
      DO 1 I=1,6
    1 WORD(I)=QBL
      IFFAIL=0
      IF(LETCOL.GT.LAST)GO TO 4002
      N=0
C     (M)ACHINE (I)NDEPENDENT (A)LPHAMERIC (C)ODE
      MIAC=0
      ASSIGN 500 TO NSWCH
      DO 2000 NCOL=LETCOL,LAST
      C=CARD(NCOL)
      GO TO NSWCH,(500,600)
  500 IF(C.EQ.QBL)GO TO 1999
      DO 501 II=1,26
      IF(C.EQ.A(II))GO TO 502
  501 CONTINUE
 1999 MIAC=0
      N=0
      ASSIGN 500 TO NSWCH
      GO TO 2000
  502 IF(C.NE.QE )GO TO 506
      IF(NCOL.EQ.1)GO TO 506
      C=CARD(NCOL-1)
      DO 503 I=27,36
      IF(C.EQ.A(I))GO TO 504
  503 CONTINUE
      IF(C.NE.QPT)GO TO 506
  504 C=CARD(NCOL+1)
      DO 505 I=27,36
      IF(C.EQ.A(I))GO TO 1999
  505 CONTINUE
      IF(C.EQ.QBL)GO TO 1999
  506 N=N+1
      MIAC=II
      WORD(1)=A(II)
      ASSIGN 600 TO NSWCH
      GO TO 2000
  600 IF(C.EQ.QBL)GO TO 603
      IF(C.EQ.QEQ)GO TO 603
      IF(C.EQ.QCM)GO TO 603
      DO 601 II=1,36
      IF(C.EQ.A(II))GO TO 602
  601 CONTINUE
      GO TO 1999
  602 N=N+1
      IF(N.GT.6)GO TO 604
      MIAC=37*MIAC+II
      WORD(N)=A(II)
  604 GO TO 2000
  603 IWORDF=MIAC
      GO TO 998
 2000 CONTINUE
      LETCOL=LAST+1
 4002 IFFAIL=1
      IF(MORE.GT.0)GO TO 4000
      WRITE(6,4001)(CARD(I),I=1,LAST)
 4001 FORMAT(28H1IWORDF HAS READ OFF THE END/(1X80A1))
      stop
 4000 IWORDF=0
      RETURN
  998 LETCOL=NCOL+1
      RETURN
      END
	subroutine ttaup(d1,d2,d3,d4,d5,d6,d7,d8,x,n)
	implicit real*8 (a-h,o-z)
	dimension d1(*),d2(*),d3(*),d4(*),d5(*),d6(*),d7(*),d8(*)
	end