**==hline.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      SUBROUTINE HLINE(Rhy,A,T,St,Beta,Alter)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL A , ALFly , ALFlyc , Alter , Beta , cc , chi , ex , EXIt , 
     &     HALf , HALfc , HBEt , HBEtc , om , P , rh , Rhy , SEATON , 
     &     SECONT , St
      REAL sum , T , tp , TWOp , y
      INTEGER k
C*** End of declarations inserted by SPAG
      COMMON /HLN   / EXIt(11) , HALf , HALfc , ALFly , ALFlyc , HBEt , 
     &                HBEtc , TWOp
      rh = Rhy/(1.+Rhy)
      y = .75*157890./T
      cc = ALOG((y+1.)/y) - 0.4/(y+1.)**2
C     CONTINUUM TO N .GE. 3
      EXIt(1) = 3.46*rh/St
C     CONTINUUM TO N .EQ. 2
      EXIt(2) = (45.0/St)*SECONT(Beta,2.)*rh/8.
C     CONTINUUM TO GROUND
      EXIt(3) = 45.0*SECONT(Beta,1.)*rh/St
      IF ( Alter.GE.1. ) EXIt(3) = 0.
C     TOTAL RECOMBINATION LINES
      chi = (.735+ALOG(Beta)+.3333/Beta)/2.
      IF ( Beta.LE..5 ) chi = 1.2*Beta + (.55+1.04*ALOG(Beta))
     &                        *Beta*Beta + (1.01*ALOG(Beta)-0.43)
     &                        *Beta**3
      chi = chi - P(Beta,1.)
      ex = .4288 + ALOG(Beta)/2. + 0.469*Beta**(-.333) - SEATON(Beta,1.)
      EXIt(4) = (2.06*1.607*13.6/St)
     &          *(ex+chi/Beta-SECONT(Beta,2.)/8.-.077)*rh
C     HBETA FROM RECOMBINATION
      EXIt(5) = .01275*EXP(-410./T)*rh/(.0001*T)**.937
C     LYMAN ALPHA FROM RECOMBINATION
      EXIt(6) = 2.06*1.6027*13.6*.75*ex*rh/St
C     TOTAL FROM COLLISIONAL EXCITATION
C     STILL USING OLD LYMAN ALPHA EXCITATION
      sum = .416*EXP(-.75*Beta)*1.195 + .07912*EXP(-.889*Beta)
     &      *1.076 + .02899*EXP(-.9375*Beta)
     &      *1.041 + .01394*EXP(-.96*Beta)
     &      *1.026 + .007799*EXP(-.972*Beta)*1.018
      EXIt(7) = 4.5E-6*1.6027E-12*(1.E23/SQRT(13.6))
     &          *sum/((1.+Rhy)*Beta**.12)
C     HBETA FROM COLLISIONAL EXCITATION : OLD H BETA
      sum = (.02899*.739*3/16.)*EXP(-.9375*Beta)
     &      *1.110 + .0745*.01394*EXP(-.96*Beta)
     &      *1.068 + .007799*.079*EXP(-.972*Beta)*1.047
      EXIt(8) = (4.5E-6*1.6027E11/SQRT(13.6))*sum/((1+Rhy)*Beta**.12)
C     LYMAN ALPHA FROM COLLISIONAL EXCITATION
      sum = .0791*EXP(-.889*Beta)*1.210 + .02899*EXP(-.9375*Beta)
     &      *1.110 + .01394*EXP(-.96*Beta)
     &      *1.068 + .007799*EXP(-.972*Beta)*1.047
      EXIt(9) = (4.5E-6*1.6027E-12*1.E23/SQRT(13.6))
     &          *sum/(Beta**.12*(1.+Rhy))
C  AGGARWAL FITS TO MCDOWELL AND CALLAWAY WITH OLD CASCADES
      tp = AMIN1(T,5.0E5)
      om = .5*(.335+1.45E-5*tp+1.39E-10*tp*tp-5.66E-15*tp*tp*tp)
      IF ( tp.GE.25000. ) om = .5*(.328+1.43E-5*tp-6.55E-12*tp*tp-
     &                         2.69E-18*tp*tp*tp)
      IF ( T.GE.5.0E5 ) om = 8.05*.276*cc
      EXIt(9) = EXIt(9) + 8.63E-6*om*10.2*1.6E11*EXP(-118000./T)
     &          /(SQRT(T)*(1.+Rhy))
C  H ALPHA EXCITATION; AGGARWAL ; CASE B, NO CASCADES
      om = (.175-1.31E-7*tp-4.08E-11*tp*tp+3.10E-15*tp*tp*tp)
      IF ( tp.GE.25000. ) om = (.138+2.50E-6*tp-4.43E-12*tp*tp+3.59E-18*
     &                         tp*tp*tp)
      EXIt(10) = 8.63E-6*om*1.89*1.6E11*EXP(-140300./T)
     &           /(SQRT(T)*(1.+Rhy))
C  TWO PHOTON FROM AGGARWAL;  NO CASCADES INCLUDED
      om = .5*(.195+1.28E-5*tp-4.69E-10*tp*tp+6.19E-15*tp*tp*tp)
      IF ( tp.GE.25000. ) om = .5*(.308+3.58E-7*tp+6.15E-13*tp*tp-
     &                         1.08E-18*tp*tp*tp)
      EXIt(11) = 8.63E-6*om*10.2*1.6E11*EXP(-118000./T)
     &           /(SQRT(T)*(1.+Rhy))
      PRINT 99001 , (EXIt(k),k=1,11)
      HALf = HALf + EXIt(5)*A*(2.654+1730./T)
      HALfc = HALfc + EXIt(10)*A
      HBEt = HBEt + EXIt(5)*A
      HBEtc = HBEtc + EXIt(8)*A
      ALFly = ALFly + EXIt(6)*A
      ALFlyc = ALFlyc + EXIt(9)*A
      TWOp = TWOp + EXIt(11)*A
      RETURN
99001 FORMAT (
     &' EXITS   REC3      REC2      REC1      RECT      RH B    RLYA    
     &  EXCT      EBET      ELYA      EH A      TWOP',/,5X,11E10.2)
      END
**==setpar.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      SUBROUTINE SETPAR(Vshock,Colcut)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL ABUnd , ABUnj , AH , AM , AN , BINmin , BINsyz , bz , BZEro , 
     &     Colcut , DENlim , DENz , DYNam , EERg , HEZi , HEZii , HMAx , 
     &     PARt , PB , PREss
      REAL RF1 , RFAct , RMAx , SI , SS , TINt , TLOw , tz , TZEro , 
     &     vs , Vshock , XMAx
      INTEGER ICOnt , ICX , IHE , ISTar , ITFin , NBIn , NI , NJ , NUM
C*** End of declarations inserted by SPAG
C  READS IN PARAMETERS FOR SHOCK RUN FROM SHPARAM.DAT
C  INITIALIZES VARIOUS CONSTANTS
C
      CHARACTER*80 LABel
      COMMON /LAB   / LABel
      COMMON /SLV   / DENz , PREss , PB , DYNam
      COMMON /SHPAR / BZEro , TINt , HEZi , HEZii , TLOw , RF1 , SS , 
     &                SI , RMAx , RFAct , NUM , ITFin , ICOnt , IHE , 
     &                ICX , NI , AH , AN , AM , XMAx , HMAx , DENlim , 
     &                EERg , ISTar
      COMMON /PARAMS/ NJ(12) , ABUnj(12) , ABUnd , BINmin , BINsyz , 
     &                NBIn
      COMMON /TZ    / TZEro , PARt
C
      READ (17,99003) LABel
      READ (17,99001) vs , DENz
      Vshock = vs*1.0E5
      DENz = 10.**DENz
      READ (17,99001) tz , bz
      TZEro = 10.**tz
      BZEro = .000001*bz
      READ (17,99001) TINt , PARt
      READ (17,99001) HEZi , HEZii
      READ (17,99001) BINmin , BINsyz
      READ (17,99001) TLOw , RF1
      READ (17,99001) SS , SI
      READ (17,99001) RMAx , RFAct
C  CONTROL INTEGERS
      READ (17,99002) NUM , ITFin
      READ (17,99002) ICOnt , IHE
      READ (17,99002) ICX , NI
      READ (17,99002) NBIn , ISTar
C  READ ELEMENTAL ABUNDANCES
      READ (17,99001) ABUnj(1) , ABUnj(2)
      READ (17,99001) ABUnj(3) , ABUnj(4)
      READ (17,99001) ABUnj(5) , ABUnj(6)
      READ (17,99001) ABUnj(7) , ABUnj(8)
      READ (17,99001) ABUnj(9) , ABUnj(10)
      READ (17,99001) ABUnj(11) , ABUnj(12)
      READ (17,99001) Colcut
      Colcut = 10.**Colcut
C  COMPUTE SHOCK CONSTANTS FOR JUMP CONDITIONS
      AH = 10.**(ABUnj(1)-12.)
      AN = 1. + AH
      AM = (1.+4.*AH)/AN
      XMAx = SQRT
     &       (1.27*DENz*Vshock*Vshock*8.*3.1416/(BZEro*BZEro*6.0225E23))
      PB = BZEro*BZEro*.125/3.1416
      DYNam = AM*DENz*Vshock*Vshock/6.0223E23
      PREss = DYNam + PB + DENz*(1.+PARt/AN)*1.38E-16*TZEro
      HMAx = (0.5*DYNam+2.*PB)/DENz + 2.50*(1.+PARt/AN)*1.38E-16*TZEro + 
     &       (PARt/AN)*13.595*1.602E-12 + (AH*HEZii/AN)*24.4*1.602E-12 + 
     &       (AH*(1.-HEZii-HEZi)/AN)*54.4*1.602E-12
      DENlim = XMAx*DENz
      EERg = (13.6*PARt+24.4*AH*HEZii+54.4*AH*(1.-HEZi-HEZii))*AN
C  PRINT OUT INPUT PARAMETERS
      PRINT 99004 , LABel
      PRINT 99005 , vs , DENz , TZEro , bz
      PRINT 99006 , PARt , HEZi , HEZii , EERg
      PRINT 99007 , NUM , ICOnt , ICX , ITFin , BINmin , BINsyz , NBIn
      PRINT 99008 , TINt , NI , SI , SS , RMAx , RFAct
      PRINT 99009 , TLOw , RF1 , IHE , ISTar
      PRINT 99010 , DENlim , AN , AH , AM
      PRINT 99011 , ABUnj
      PRINT 99012 , XMAx , PB , DYNam , HMAx , PREss
      PRINT 99013 , Colcut
      RETURN
C
99001 FORMAT (2F10.3)
99002 FORMAT (2I5)
99003 FORMAT (A80)
99004 FORMAT ('1',50X,'++++  SHOCK  ++++',//,A80)
99005 FORMAT (/,F10.1,' KM/S  ;   DENZ= ',E12.3,'  CM-3  ; TZ= ',E12.3,
     &        '  ; BZERO = ',E12.2,' MICROGAUSS')
99006 FORMAT (/,' IONIZED H, NEUTRAL HE, HE+',3F10.5,' EERG=',F5.2)
99007 FORMAT (/,' NUM,ICONT,ICX,ITFIN,BINMIN,BINSYZ,NBIN ',4I5,2F5.2,I5)
99008 FORMAT (/,' TINT,NI,SI,SS,RMAX,RFACT  ',F5.2,I5,4F10.3)
99009 FORMAT (/,' TLOW,RF1,IHE,ISTAR   ',2F5.2,2I5)
99010 FORMAT (/,' DENLIM,AN,AH,AM ',E12.2,3X,3F10.2)
99011 FORMAT (
     &'0   HE     C      N      O      NE     MG     SI     S    AR     
     &CA     FE     NI',/12F7.2)
99012 FORMAT (/,' XMAX=',E12.3,' PB=',E12.3,' DYNAM=',E12.3,' HMAX=',
     &        E12.3,' PRESS=',E12.3)
99013 FORMAT ('  CUTOFF AT LIMITING HYDROGEN COLUMN DENSITY ',E12.3)
      END
**==adder.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      SUBROUTINE ADDER(A1,A2,It,Itfin,Num,Denz,Vshock)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL A1 , A2 , ABIn , ABUnd , ABUnj , ah , ALFj , an , BIN , 
     &     BINmin , BINsyz , C1J , C2J , CJ , CNC , COLumn , CONce , 
     &     Denz , E3J , EJ
      REAL ESJ , ETWO , FCOol , FJ , FOUt , FTOt , FWV , GNDrec , 
     &     HENeut , HEPlus , hnu , PHOt2 , PHOt3 , PHOt4 , POWer , 
     &     PTOt , PWR , RF , RHY , SIGj
      REAL SJ , t1 , t2 , t3 , ta , TAU , TAUhe , TAUP , tb , TBIn , 
     &     TPLus , Vshock , w , WJ
      INTEGER i , It , Itfin , j , k , l , LLJ , NBIn , NJ , no , Num
C*** End of declarations inserted by SPAG
C  ADDS EMISSION, OPTICAL DEPTH AND COLUMN DENSITIES FOR ONE STEP TO TOT
      INCLUDE 'rayspec.inc'
      COMMON /ADDER / TBIn(500) , COLumn(10,6) , PHOt2 , PHOt3 , PHOt4
      COMMON /FEL   / WJ(12,220) , FJ(12,220) , E3J(12,220) , 
     &                PWR(12,220) , C1J(12,30) , C2J(12,30) , EJ(12,30)
     &                , SJ(12,30) , CJ(12,30) , LLJ(12,30) , SIGj(12,30)
     &                , ALFj(12,30) , ESJ(12,30)
      COMMON /PARAMS/ NJ(12) , ABUnj(12) , ABUnd , BINmin , BINsyz , 
     &                NBIn
      COMMON /COM   / CNC(12,30) , PTOt(12,220) , ABIn(MXBINS) , 
     &                BIN(MXBINS)
      COMMON /FLN   / FWV(320) , FTOt(320) , FOUt(320) , FCOol
      COMMON /RESULT/ CONce(30) , GNDrec(30) , POWer(220) , RHY , 
     &                HENeut , HEPlus
      COMMON /PT    / RF(500) , TAU(150) , TAUhe(150) , TPLus(150)
C
      ah = 10.**(ABUnj(1)-12.)
      an = 1. + ah
C
C  COLUMN DENSITIES
C
      DO 100 no = 1 , 10
         DO 50 j = 1 , 6
            COLumn(no,j) = COLumn(no,j) + A1*10.**(ABUnj(no)-12.)
     &                     *CNC(no,j)
 50      CONTINUE
 100  CONTINUE
C
C  OPTICAL DEPTHS
C
      t1 = 6.3E-18*A1/(an*(1.+RHY))
      t2 = ah*A1*HENeut*8.32E-18/(1.+ah)
      t3 = ah*A1*HEPlus*1.75E-18/(1.+ah)
      DO 200 k = It , Itfin
         TAU(k) = TAU(k) + t1
         TAUhe(k) = TAUhe(k) + t2
         TPLus(k) = TPLus(k) + t3
 200  CONTINUE
C
C  BINS, PHOT2,PHOT2,PHOT4
C
      DO 300 i = 1 , NBIn
         hnu = BINmin + (i-.05)*BINsyz
         w = 12399./hnu
         ta = TAUP(w,It)
         tb = TAUP(w,Itfin)
         ABIn(i) = ABIn(i)*EXP(-ta) + A2*BIN(i)*ETWO(ta)
         PHOt3 = PHOt3 + BIN(i)*A2*(1.-EXP(-tb))
     &           /(hnu*1.6E-12*Denz*Vshock)
         IF ( hnu.GE.24.4 ) PHOt4 = PHOt4 + BIN(i)*A2*(1.-EXP(-tb))
     &                              /(hnu*1.6E-12*Denz*Vshock)
         TBIn(i) = TBIn(i) + BIN(i)*A2*EXP(-tb)
 300  CONTINUE
C
C  PERMITTED LINES
C
      DO 400 no = 1 , Num
         DO 350 l = 1 , 220
            w = WJ(no,l)
            ta = 0.
            IF ( w.GT.0. ) ta = TAUP(w,Itfin)
            PTOt(no,l) = PTOt(no,l) + PWR(no,l)*A2*EXP(-ta)
 350     CONTINUE
 400  CONTINUE
C
C  FORBIDDEN LINES
C
      DO 500 i = 1 , 320
         ta = 0.
         IF ( FWV(i).GT.0. ) ta = TAUP(FWV(i),Itfin)
         FTOt(i) = FTOt(i) + A2*FOUt(i)*EXP(-ta)
 500  CONTINUE
      RETURN
      END
**==secont.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      FUNCTION SECONT(B,Q)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL a , B , GAMMA , Q , r , SECONT
C*** End of declarations inserted by SPAG
      a = B/(Q*Q)
      r = Q**.6666667
      IF ( a.LE.50.0 ) THEN
         SECONT = 1.0 + (.1728/r)*((.33333-2.*a)*GAMMA(.33333,a)+1.0)
     &            - (.0496/(r*r))
     &            *(.66667*(1.-a-3.*a*a)*GAMMA(.66667,a)+(1.+2.*a))
      ELSE
         SECONT = (1.0-.1728/r) - 0.0496/(r*r)
      ENDIF
      RETURN
      END
**==lnprt.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      SUBROUTINE LNPRT(Num)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL ABIn , ABUnd , ABUnj , ALFj , ALFly , ALFlyc , BIN , BINmin , 
     &     BINsyz , c , C1J , C2J , CJ , CNC , E3J , EJ , ESJ , EXIt , 
     &     FCOol , FJ
      REAL FOUt , FTOt , FWV , HALf , HALfc , HBEt , HBEtc , PDUm , 
     &     PTOt , PWR , SIGj , SJ , TWOp , WDUm , WJ
      INTEGER i , il , ix , iy , j , JDUm , JMX , jt , k , l , LDUm , 
     &        LLJ , loc , n , NBIn , NF , nhigh , NJ , nlow , no
      INTEGER Num
C*** End of declarations inserted by SPAG
      INCLUDE 'rayspec.inc'
      COMMON /HLN   / EXIt(11) , HALf , HALfc , ALFly , ALFlyc , HBEt , 
     &                HBEtc , TWOp
      COMMON /PARAMS/ NJ(12) , ABUnj(12) , ABUnd , BINmin , BINsyz , 
     &                NBIn
      COMMON /FLN   / FWV(320) , FTOt(320) , FOUt(320) , FCOol
      COMMON /FEL   / WJ(12,220) , FJ(12,220) , E3J(12,220) , 
     &                PWR(12,220) , C1J(12,30) , C2J(12,30) , EJ(12,30)
     &                , SJ(12,30) , CJ(12,30) , LLJ(12,30) , ALFj(12,30)
     &                , ESJ(12,30) , SIGj(12,30)
      COMMON /COM   / CNC(12,30) , PTOt(12,220) , ABIn(MXBINS) , 
     &                BIN(MXBINS)
      INTEGER bb
      COMMON /LNPR  / JMX(11) , NF(11) , LDUm(5) , JDUm(5) , WDUm(5) , 
     &                PDUm(5)
      PRINT 99001 , HBEt , HBEtc , TWOp
      c = 100./(HBEt+HBEtc)
      DO 100 no = 1 , Num
         n = NJ(no)
         PRINT 99002 , n
         loc = 1
         ix = 0
         jt = AMIN0(JMX(no),n)
         DO 50 j = 1 , jt
            iy = LLJ(no,j)*3
            IF ( iy.GT.0 ) THEN
               DO 10 l = 1 , iy
                  bb = ix + l
                  IF ( WJ(no,bb).GT.0 ) THEN
                     LDUm(loc) = l
                     JDUm(loc) = j
                     WDUm(loc) = WJ(no,bb)
                     PDUm(loc) = PTOt(no,bb)*c
                     IF ( loc.EQ.5 ) PRINT 99003 , n , 
     &                    ((JDUm(k),LDUm(k),WDUm(k),PDUm(k)),k=1,5)
                     loc = MOD(loc,5) + 1
                  ENDIF
 10            CONTINUE
            ENDIF
            ix = ix + iy
 50      CONTINUE
         il = loc - 1
         IF ( loc.NE.1 ) PRINT 99003 , n , 
     &                         ((JDUm(k),LDUm(k),WDUm(k),PDUm(k)),k=1,
     &                         il)
         IF ( n.NE.2 ) THEN
            PRINT 99004 , n
            nlow = NF(no-1) + 1
            nhigh = NF(no)
            loc = 1
            DO 60 i = nlow , nhigh
               IF ( FWV(i)*FTOt(i).GT.0. ) THEN
                  LDUm(loc) = i
                  WDUm(loc) = FWV(i)
                  PDUm(loc) = FTOt(i)*c
                  IF ( loc.EQ.5 ) PRINT 99005 , n , 
     &                                  ((LDUm(k),WDUm(k),PDUm(k)),k=1,
     &                                  5)
                  loc = MOD(loc,5) + 1
               ENDIF
 60         CONTINUE
            il = loc - 1
            IF ( loc.NE.1 ) PRINT 99005 , n , 
     &                            ((LDUm(k),WDUm(k),PDUm(k)),k=1,il)
         ENDIF
 100  CONTINUE
      RETURN
99001 FORMAT ('    HBET, HBETC, TWOP    ',3E10.3)
99002 FORMAT (' LINES FOR ELEMENT',I5)
99003 FORMAT (I4,5(I4,I3,F8.2,E10.3))
99004 FORMAT ('  ***** FORBIDDEN LINES ***** N=',I3)
99005 FORMAT (1X,I2,5(I6,F9.1,E10.3))
      END
**==solvx.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      FUNCTION SOLVX(F,E,H)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL ax0 , ax1 , ax2 , DENz , DYNam , E , F , H , PB , plink , 
     &     PREss , qx , rx , SOLVX , thetax , x
      INTEGER k
C*** End of declarations inserted by SPAG
      COMMON /SLV   / DENz , PREss , PB , DYNam
C     CALCULATE COMPRESSION
      ax2 = 2.*DENz*(H-E)/PB
      ax1 = -5.*PREss/PB
      ax0 = 4.*DYNam/PB
      qx = ax2*ax2/9. - ax1/3.
      rx = (3.*ax0-ax1*ax2)/6. + (ax2**3)/27.
      plink = rx/qx**1.5
      IF ( plink.LE.0.997 ) THEN
         thetax = ACOS(plink)/3.
         SOLVX = SQRT(qx)*(COS(thetax)+1.732*SIN(thetax)) - ax2/3.
         RETURN
      ENDIF
      x = (-ax1+SQRT(ax1*ax1-4.*ax2*ax0))/(2.*ax2)
      DO 100 k = 1 , 4
         x = (-ax1+SQRT(ax1*ax1-4.*ax2*(ax0+x**3)))/(2.*ax2)
 100  CONTINUE
      SOLVX = x
      RETURN
      END
**==hcool.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      FUNCTION HCOOL(T,Rhy,Alter)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL ALFly , ALFlyc , Alter , beta , EXIt , HALf , HALfc , HBEt , 
     &     HBEtc , HCOOL , Rhy , st , T , TWOp
C*** End of declarations inserted by SPAG
      COMMON /HLN   / EXIt(11) , HALf , HALfc , ALFly , ALFlyc , HBEt , 
     &                HBEtc , TWOp
      beta = 157890./T
      st = SQRT(T)
      CALL HLINE(Rhy,0,T,st,beta,Alter)
      HCOOL = EXIt(1) + EXIt(2) + EXIt(3) + EXIt(4) + EXIt(8) + EXIt(9)
     &        + EXIt(10) + EXIt(11)
      RETURN
      END
**==taup.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      FUNCTION TAUP(W,It)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      INTEGER It
      REAL RF , TAU , TAUhe , TAUP , TPLus , W
C*** End of declarations inserted by SPAG
      COMMON /PT    / RF(500) , TAU(150) , TAUhe(150) , TPLus(150)
      TAUP = 1.E-8
      IF ( W.GT.912. ) RETURN
      TAUP = TAUP + TAU(It)*(W/912.)**3
      IF ( W.LT.504.6 ) TAUP = TAUP + TAUhe(It)
     &                         *(.763*(W/504.6)**1.99+.237*(W/504.6)
     &                         **2.99)
      IF ( W.LT.228. ) TAUP = TAUP + TPLus(It)*(W/228.)**3
      RETURN
      END
**==expf.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      FUNCTION EXPF(X)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL EXPF , X
C*** End of declarations inserted by SPAG
      EXPF = 1.0E-37
      IF ( X.GE.-75. .AND. X.LE.75. ) EXPF = EXP(X)
      IF ( X.GT.75. ) EXPF = 1.0E37
      RETURN
      END
**==gamma.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      FUNCTION GAMMA(A,X)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL A , fn , fn1 , fn2 , GAMMA , X
      INTEGER in , mn , n , nn
C*** End of declarations inserted by SPAG
      REAL*8 can1 , can2 , cbn1 , cbn2 , an , bn , can , cbn
C     EXP(-X) * X ** A   TAKEN OUTSIDE
      can2 = 1.
      can1 = 0.
      cbn2 = 0.
      cbn1 = 1.
      an = 1.
      bn = X
      n = 1
      can = bn*can1 + an*can2
      cbn = bn*cbn1 + an*cbn2
      fn1 = can/cbn
      DO 100 n = 2 , 1000
         IF ( can.GT.1.E30 ) THEN
            can1 = can1*1.E-30
            can = can*1.E-30
            cbn1 = cbn1*1.E-30
            cbn = cbn*1.E-30
         ENDIF
         in = n
         can2 = can1
         cbn2 = cbn1
         can1 = can
         cbn1 = cbn
         mn = MOD(n,2)
         bn = mn*X + (1.-mn)
         nn = n/2
         an = nn - (1.-mn)*A
         can = bn*can1 + an*can2
         cbn = bn*cbn1 + an*cbn2
         fn = can/cbn
         IF ( ABS((fn-fn1)/fn1).LE..0001 .AND. n.GE.20 ) GOTO 200
         fn2 = fn1
         fn1 = fn
 100  CONTINUE
      GAMMA = (fn+fn2)/2.
      RETURN
 200  GAMMA = fn
      RETURN
      END
**==cycle.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      SUBROUTINE CYCLE(N,Dene,Elapse)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL ageinv , apop , blatz , C1 , C2 , CA , CTH , d1dwn , decay , 
     &     delta , Dene , DNE , DONce , down , dp , dt , du , dup , dw , 
     &     Elapse
      REAL EXPF , FABLE , fast , HENeut , HEPlus , p1 , p2 , pa , 
     &     PCOol , PM , POT , POU , POWer , RE , RHY , stay , TU , 
     &     turn , up
      INTEGER iterm , j , jot , k , lm , N
C*** End of declarations inserted by SPAG
      COMMON /RESULT/ DONce(30) , CA(30) , POWer(220) , RHY , HENeut , 
     &                HEPlus , DNE , PCOol , POU , POT , RE , TU , PM(4)
      COMMON /RATES / C1(30) , C2(30) , CTH
      DIMENSION p1(30) , p2(30) , pa(30)
      DIMENSION delta(30) , up(30,4) , stay(30) , down(30,2) , 
     &          ageinv(30) , turn(30) , decay(30) , dup(30,2) , 
     &          d1dwn(30) , du(30) , dp(30) , dw(30) , apop(30)
      jot = N + 1
      C2(1) = 0
      C1(jot) = 0
      p2(N+2) = 0.
      fast = 0
      CA(jot) = 0
      CA(jot-1) = 0
      ageinv(N+2) = 0
      ageinv(N+3) = 0
      DO 100 j = 1 , jot
         ageinv(j) = (C1(j)+C2(j)+CA(j))*Dene
         blatz = DONce(j)*ageinv(j)
         IF ( ageinv(j).LE.0 ) THEN
            p1(j) = 0
            p2(j) = 0
            pa(j) = 0
         ELSE
            p1(j) = C1(j)*Dene/ageinv(j)
            p2(j) = C2(j)*Dene/ageinv(j)
            pa(j) = CA(j)*Dene/ageinv(j)
         ENDIF
         IF ( DONce(j).LE..0001 ) blatz = 0
         fast = AMAX1(fast,blatz)
 100  CONTINUE
      iterm = 100
      IF ( 5.*fast*Elapse.LE.100. ) iterm = 5*fast*Elapse
      IF ( N.LE.1 ) THEN
         DONce(1) = ageinv(2)/(ageinv(1)+ageinv(2))
     &              + (DONce(1)-ageinv(2)/(ageinv(1)+ageinv(2)))
     &              *EXPF(-Elapse*(ageinv(1)+ageinv(2)))
         DONce(2) = 1.0 - DONce(1)
      ELSE
         iterm = MIN0(iterm,100)
         IF ( iterm.LE.6 ) iterm = 6
         dt = Elapse/iterm
         DO 150 j = 1 , jot
            IF ( j.LE.1 ) THEN
               d1dwn(1) = 0
            ELSE
               d1dwn(j) = C2(j)*Dene*FABLE(ageinv(j-1),ageinv(j),dt)
            ENDIF
            turn(j) = ageinv(j)*dt
            decay(j) = EXPF(-turn(j))
            dup(j,1) = C1(j)*Dene*FABLE(ageinv(j+1),ageinv(j),dt)
            dup(j,2) = CA(j)*Dene*FABLE(ageinv(j+2),ageinv(j),dt)
            du(j) = p1(j)*(1.-decay(j)) - dup(j,1)
            dp(j) = p2(j)*(1.-decay(j)) - d1dwn(j)
            dw(j) = pa(j)*(1.-decay(j)) - dup(j,2)
 150     CONTINUE
         DO 250 j = 1 , jot
            IF ( j.EQ.jot ) THEN
               stay(j) = decay(j) + p1(j-1)*dp(j)
               DO 160 k = 1 , 4
                  up(j,k) = 0
 160           CONTINUE
            ELSEIF ( j.EQ.(jot-1) ) THEN
               IF ( j.EQ.1 ) THEN
                  up(1,1) = dup(1,1)
               ELSE
                  up(j,1) = dup(j,1) + pa(j-1)*dp(j)
               ENDIF
               DO 170 k = 2 , 4
                  up(j,k) = 0
 170           CONTINUE
            ELSEIF ( j.EQ.(jot-2) ) THEN
               DO 180 k = 3 , 4
                  up(j,k) = 0
 180           CONTINUE
               GOTO 200
            ELSEIF ( j.EQ.(jot-3) ) THEN
               up(j,4) = 0
               up(j,3) = p1(j+2)*dw(j) + pa(j+1)*du(j)
               GOTO 200
            ELSE
               IF ( j.EQ.1 ) THEN
                  up(j,1) = dup(j,1) + p2(j+2)*dw(j)
               ELSE
                  up(j,1) = dup(j,1) + pa(j-1)*dp(j) + p2(j+2)*dw(j)
               ENDIF
               up(j,4) = pa(j+2)*dw(j)
               up(j,3) = p1(j+2)*dw(j) + pa(j+1)*du(j)
               up(j,2) = dup(j,2) + p1(j+1)*du(j)
            ENDIF
            GOTO 220
 200        IF ( j.EQ.1 ) THEN
               up(1,1) = dup(1,1) + p2(j+2)*dw(j)
            ELSE
               up(j,1) = dup(j,1) + pa(j-1)*dp(j) + p2(j+2)*dw(j)
            ENDIF
            up(j,2) = dup(j,2) + p1(j+1)*du(j)
 220        IF ( j.EQ.1 ) THEN
               stay(j) = decay(j) + p2(j+1)*du(j)
            ELSE
               stay(j) = decay(j) + p1(j-1)*dp(j) + p2(j+1)*du(j)
            ENDIF
            IF ( j.EQ.1 ) THEN
               down(j,1) = 0
               down(j,2) = 0
            ELSEIF ( j.EQ.2 ) THEN
               down(j,2) = 0
               down(j,1) = d1dwn(j)
            ELSE
               down(j,1) = d1dwn(j)
               down(j,2) = p2(j-1)*dp(j)
            ENDIF
 250     CONTINUE
         DO 300 lm = 1 , 30
            IF ( lm.LE.(N+1) ) THEN
               apop(lm) = DONce(lm)
            ELSE
               apop(lm) = 0
               down(lm,1) = 0
               down(lm,2) = 0
            ENDIF
 300     CONTINUE
         DO 350 k = 1 , iterm
            DO 320 j = 1 , jot
               delta(j) = apop(j)*stay(j) + apop(j+1)*down(j+1,1)
     &                    + apop(j+2)*down(j+2,2)
               IF ( j.NE.1 ) THEN
                  IF ( j.EQ.2 ) THEN
                     delta(j) = delta(j) + apop(j-1)*up(j-1,1)
                  ELSEIF ( j.EQ.3 ) THEN
                     delta(j) = delta(j) + apop(j-1)*up(j-1,1)
     &                          + apop(j-2)*up(j-2,2)
                  ELSEIF ( j.EQ.4 ) THEN
                     delta(j) = delta(j) + apop(j-1)*up(j-1,1)
     &                          + apop(j-2)*up(j-2,2) + apop(j-3)
     &                          *up(j-3,3)
                  ELSE
                     delta(j) = delta(j) + apop(j-4)*up(j-4,4)
     &                          + apop(j-3)*up(j-3,3) + apop(j-2)
     &                          *up(j-2,2) + apop(j-1)*up(j-1,1)
                  ENDIF
               ENDIF
 320        CONTINUE
            DO 340 j = 1 , jot
               apop(j) = delta(j)
               delta(j) = 0
 340        CONTINUE
 350     CONTINUE
         DO 400 j = 1 , jot
            DONce(j) = apop(j)
 400     CONTINUE
      ENDIF
C	PRINT 998,DONCE
      RETURN
C	PRINT 999,N,DENE,ELAPSE
C	PRINT 998,DONCE
C	PRINT 998,C1
C	PRINT 998,C2
C	PRINT 998,CA
99001 FORMAT (I5,2E12.3)
99002 FORMAT (10E12.3)
      END
**==fable.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      FUNCTION FABLE(X,Y,Z)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL duck , EXPF , FABLE , X , Y , Z
C*** End of declarations inserted by SPAG
      IF ( X.NE.Y ) THEN
         duck = (EXPF(-X*Z)-EXPF(-Y*Z))/(Y-X)
      ELSE
         duck = Z*EXP(-X*Z)
      ENDIF
      FABLE = duck
      RETURN
      END
**==fsolv.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      SUBROUTINE FSOLV(Nl,T,Dene,Lev,Ind)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL A21 , A31 , a4 , A41 , A42 , A43 , a5 , A51 , A52 , A53 , 
     &     AMG , amp , au1 , au2 , au3 , AUG , AUM , aump , CON , d2
      REAL d3 , deg1 , deg2 , deg3 , deg4 , deg5 , degg , degm , degu , 
     &     delta , Dene , EM , EN , EN1 , EN2 , EN3 , EN4 , EN5 , EU , 
     &     FCOol
      REAL FCT , FOUt , FTOt , FWV , g23 , g32 , g45 , OH23 , OH45 , 
     &     OHM , OHMg , OHUg , OHUm , omp , pg , pm , pu , qdown , qgm , 
     &     qgu
      REAL qmg , qmm , qmu , quem , queu , qug , qum , qup , quu , r2 , 
     &     r3 , r4 , r5 , rm , ru , s , SQ , T
      INTEGER Ind , JDUm , k , Lev , NDUm , Nl
C*** End of declarations inserted by SPAG
      COMMON /F2    / NDUm , JDUm , OHM , EN
      COMMON /F     / OHUm , OHUg , OHMg , AUM , AUG , AMG , EU , EM , 
     &                SQ , CON , FCT , EN1 , EN2 , EN3 , EN4 , EN5 , 
     &                A53 , A52 , A51 , A43 , A42 , A41 , A31 , A21 , 
     &                OH45 , OH23
      COMMON /FLN   / FWV(320) , FTOt(320) , FOUt(320) , FCOol
      DIMENSION d2(2,6) , d3(3,8)
      DATA d2/2. , 4. , 4. , 2. , 6. , 12. , 9. , 5. , 1. , 9. , 2. , 
     &     10./
      DATA d3/9. , 5. , 1. , 1. , 3. , 5. , 5. , 3. , 1. , 10. , 8. , 
     &     6. , 9. , 7. , 5. , 9. , 9. , 33. , 4. , 6. , 8. , 4. , 12. , 
     &     22./
      IF ( Nl.EQ.2 ) THEN
         qdown = 8.63E-6*OHM/(d2(2,Ind)*SQ)
         qup = (d2(2,Ind)/d2(1,Ind))*qdown*EXP(-1.4388*EN/T)
         ru = qup/(A21/Dene+qup+qdown)
         FWV(Lev) = 1.E+8/EN
         FOUt(Lev) = ru*A21*EN*FCT*CON
         FCOol = FCOol + FOUt(Lev)
         Lev = Lev + 1
         RETURN
      ELSEIF ( Nl.EQ.3 ) THEN
C     THREE LEVEL PACKAGE
         degg = d3(1,Ind)
         degm = d3(2,Ind)
         degu = d3(3,Ind)
         IF ( Ind.EQ.3 .AND. NDUm-JDUm.EQ.7 .AND. NDUm.GE.26 ) THEN
            degm = 1.
            degu = 3.
         ENDIF
      ELSEIF ( Nl.EQ.4 ) THEN
         RETURN
      ELSEIF ( Nl.EQ.5 ) THEN
C     FIVE LEVEL PACKAGE
         deg1 = 4.
         deg2 = 4.
         deg3 = 6.
         deg4 = 2.
         deg5 = 4.
         IF ( NDUm-JDUm.EQ.6 ) deg2 = 6.
         IF ( NDUm-JDUm.EQ.6 ) deg3 = 4.
         IF ( Ind.EQ.4 ) THEN
            deg1 = 21.
            deg2 = 2.
            deg3 = 3.
            deg4 = 9.
            deg5 = 9.
         ENDIF
         EM = (EN2*deg2+EN3*deg3)/(deg1+deg2)
         degu = deg4 + deg5
         degm = deg2 + deg3
         degg = deg1
         a4 = A43 + A42 + A41
         a5 = A53 + A52 + A51
         qmg = 8.63E-6*OHMg/(degm*SQ)
         queu = 8.63E-6*(OHUm+OHUg)/(degu*SQ)
         quem = qmg + 8.63E-6*OHUm*EXP(-1.4388*(EU-EM)/T)/(degm*SQ)
         quu = (degu/(deg4*deg5))*8.63E-6*OH45/SQ
         qmm = degm/(deg2*deg3)*8.63E-6*OH23/SQ
         g45 = (deg4/deg5)*(a5+Dene*(queu+quu))/(a4+Dene*(queu+quu))
         g23 = (deg2/deg3)*(A31+Dene*(quem+qmm))/(A21+Dene*(quem+qmm))
         g32 = (Dene*deg3*qmm/degm)/(A31+Dene*(quem+qmm))
         au3 = (A53+g45*A43)/(1.+g45)
         au2 = (A52+g45*A42)/(1.+g45)
         au1 = (A51+g45*A41)/(1.+g45)
         AMG = (A31+g23*A21)/(1.+g23)
         amp = (A21+g32*A31)/(1.+g32)
         omp = (A21+Dene*quem+(A31+Dene*quem)*g32)/(1.+g32)
         aump = au2 - deg2*au3/deg3
         AUM = degm*au3/deg3
C        CAN NOW BE DONE AS THREE LEVEL PROBLEM
         AUG = au1 + aump*(amp+Dene*qmg)/omp
      ELSE
         RETURN
      ENDIF
      qum = 8.63E-6*OHUm/(degu*SQ)
      qmu = (degu/degm)*qum*EXP(-(EU-EM)*1.4388/T)
      qug = 8.63E-6*OHUg/(degu*SQ)
      qgu = qug*EXP(-EU*1.4388/T)*degu/degg
      qmg = 8.63E-6*OHMg/(degm*SQ)
      qgm = qmg*EXP(-EM*1.4388/T)*degm/degg
      pu = (qgm*qmu*Dene**2+(AMG+qmu*Dene+qmg*Dene)*qgu*Dene)
      pm = (qgu*Dene*(AUM+qum*Dene)+qgm*Dene*(AUM+AUG+qum*Dene+qug*Dene)
     &     )
      pg = ((AUM+AUG+(qum+qug)*Dene)*(AMG+qmu*Dene+qmg*Dene)-(qmu*Dene)
     &     *(AUM+qum*Dene))
      s = pu + pm + pg
      ru = pu/s
      rm = pm/s
      IF ( Nl.GE.4 ) THEN
         delta = ru*aump/omp
         ru = ru/(1.+delta)
         rm = rm/(1.+delta)
         r5 = ru/(1.+g45)
         r4 = g45*r5
         r3 = rm/(1.+g23) + g32*delta/((1.+delta)*(1.+g32))
         r2 = g23*rm/(1.+g23) + delta/((1.+delta)*(1.+g23))
         IF ( Ind.NE.4 ) THEN
            FWV(Lev) = 1.E+8/EU
            FOUt(Lev) = (r5*A51+r4*A41)*EU*FCT*CON
            FWV(Lev+1) = 1.E+8/EN3
            FOUt(Lev+1) = r3*A31*EN3*FCT*CON
            FWV(Lev+2) = 1.E+8/EN2
            FOUt(Lev+2) = r2*A21*EN2*FCT*CON
            FWV(Lev+3) = 1.E+8/(EU-EN2)
            FOUt(Lev+3) = (r5*A52+r4*A42)*(EU-EN2)*FCT*CON
            FWV(Lev+4) = 1.E+8/(EU-EN3)
            FOUt(Lev+4) = (r5*A53+r4*A43)*(EU-EN3)*FCT*CON
            DO 20 k = 1 , 5
               FCOol = FCOol + FOUt(Lev-1+k)
 20         CONTINUE
            Lev = Lev + 5
            RETURN
         ENDIF
      ELSE
         FWV(Lev) = 1.E+8/(EU-EM)
         FOUt(Lev) = ru*(EU-EM)*AUM*FCT*CON
         IF ( Ind.EQ.6 .OR. Ind.EQ.8 ) FWV(Lev) = 1.E+8/EU
         IF ( Ind.EQ.6 .OR. Ind.EQ.8 ) FOUt(Lev) = ru*EU*AUG*FCT*CON
         FWV(Lev+1) = 1.E+8/EM
         FOUt(Lev+1) = rm*EM*AMG*FCT*CON
         FCOol = FCOol + FOUt(Lev)
         FCOol = FCOol + FOUt(Lev+1)
         Lev = Lev + 2
         IF ( Ind.GE.2 ) RETURN
         FWV(Lev) = 1.E+8/EU
         FOUt(Lev) = ru*EU*AUG*FCT*CON
         FCOol = FCOol + FOUt(Lev)
         Lev = Lev + 1
         RETURN
      ENDIF
      FWV(Lev) = 1.E+8/EN5
C   FE VII
      FWV(Lev+1) = 1.E+8/EN4
      FWV(Lev+2) = 1.E+8/(EN4-EN2)
      FWV(Lev+3) = 1.E+8/EN2
      FOUt(Lev) = r5*A51*EN5*FCT*CON
      FOUt(Lev+1) = r4*A41*EN4*FCT*CON
      FOUt(Lev+2) = r4*(A42+A43)*(EN4-EN2)*FCT*CON
      FOUt(Lev+3) = (r3+r2)*A21*EN2*FCT*CON
      DO 100 k = 1 , 4
         FCOol = FCOol + FOUt(Lev-1+k)
 100  CONTINUE
      Lev = Lev + 4
      RETURN
      END
**==f12.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      FUNCTION F12(Wi,W1)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL F12 , W1 , Wi
C*** End of declarations inserted by SPAG
      F12 = 0.70
      IF ( W1.LE.Wi ) F12 = 0.13
      IF ( 1.8*W1.LE.Wi ) F12 = 0.049
      IF ( 2.4*W1.LE.Wi ) F12 = 0.023
      RETURN
      END
**==fsf.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      SUBROUTINE FSF(T,Dene,Elohy)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL A21 , A31 , A41 , A42 , A43 , A51 , A52 , A53 , aalf , ABIn , 
     &     ABUnd , ABUnj , acarb , acinf , acl , aco , aff , AMG , an , 
     &     aoi
      REAL aoo , ap , asi , asii , asio , aso , AUG , AUM , BIN , 
     &     BINmin , BINsyz , CNC , CON , Dene , dne , e3d , e3p2 , 
     &     ealf , eb , ecinf
      REAL ecl , ecopt , eff , efn , Elohy , EM , EN , EN1 , EN2 , EN3 , 
     &     EN4 , EN5 , eoinf , eoopt , eosf , ep , esf , esi , esii , 
     &     esio
      REAL eso , EU , EXINT1 , FCOol , FCT , FIR , FOUt , FTOt , FWV , 
     &     ocinf , ocl , oco , OH23 , OH45 , OHM , OHMg , ohn , OHUg , 
     &     OHUm , omalf
      REAL ombf , omff , ompal , ompb , ompbsl , ompcl , ompcls , 
     &     ompeff , omsf , ooi , ooo , op , osi , osii , osiopt , oso , 
     &     POWer , realn , SQ , T
      REAL t4 , x
      INTEGER i , j , JDUm , jtop , lev , levno , n , NBIn , NDUm , NJ , 
     &        nn , no
C*** End of declarations inserted by SPAG
C O II AND S II FROM PRADHAN, MENDOZA & ZEIPPEN PAPERS
C  O I AND C I FROM PEQ AND ALD; N I FROM BERRINGTON&BURKE
C  NE V FROM GILES
C
C  WON'T DO NI BUT NOW DOES ARGON
C
      INCLUDE 'rayspec.inc'
      COMMON /COM   / CNC(12,30) , POWer(12,220) , ABIn(MXBINS) , 
     &                BIN(MXBINS)
      COMMON /FLN   / FWV(320) , FTOt(320) , FOUt(320) , FCOol
      COMMON /PARAMS/ NJ(12) , ABUnj(12) , ABUnd , BINmin , BINsyz , 
     &                NBIn
      COMMON /F     / OHUm , OHUg , OHMg , AUM , AUG , AMG , EU , EM , 
     &                SQ , CON , FCT , EN1 , EN2 , EN3 , EN4 , EN5 , 
     &                A53 , A52 , A51 , A43 , A42 , A41 , A31 , A21 , 
     &                OH45 , OH23
      COMMON /F2    / NDUm , JDUm , OHM , EN
      DIMENSION ombf(10) , eb(10) , ecopt(2,10) , aco(3,10) , oco(3,10)
     &          , ecinf(2,10) , acinf(3,10) , ocinf(3,10) , an(8,9) , 
     &          ohn(5,9) , efn(3,9) , eosf(8) , eoopt(2,8) , aoo(3,8) , 
     &          ooo(3,8) , eoinf(2,8) , aoi(3,8) , ooi(3,8) , aff(7) , 
     &          eff(7) , omff(7) , esf(7) , omsf(7) , aalf(5) , ealf(5)
     &          , omalf(5) , asio(3,5) , osiopt(3,5) , esio(2,5) , 
     &          esii(2,5) , asii(3,5) , osii(3,5) , op(5,4) , ep(3,4) , 
     &          aso(3,3) , oso(3,3) , eso(2,3) , asi(3,3) , osi(3,3) , 
     &          esi(2,3)
      DIMENSION acl(3) , ecl(3) , ocl(3)
      DIMENSION ompal(5) , ompb(10) , ompbsl(10) , ompcl(3) , ompcls(3)
      DIMENSION acarb(10) , levno(12)
      DATA levno/0 , 1 , 10 , 24 , 43 , 67 , 91 , 124 , 167 , 213 , 
     &     272 , 317/
      DIMENSION ap(8,4)
      DATA ap/77.8 , 184. , 354. , .33 , 69.9 , 186. , 2.22 , 53. , 
     &     2.88 , 4.25 , 12. , .39 , 2.94 , 4.8 , .0196 , 2.68 , .79 , 
     &     .60 , 2.11 , .12 , .60 , .86 , .0018 , .022 , .179 , .133 , 
     &     .225 , .0779 , .163 , .0906 , .00026 , .000882/
      DATA acarb/9000. , 300. , 80. , 21. , 3.07 , .324 , .0202 , 
     &     .000518 , .0000477 , .00000236/
      REAL ln
      DATA ombf/.033 , .066 , .10 , .134 , .184 , .272 , .43 , .80 , 
     &     1.09 , 1.44/
      DATA eb/118200. , 39200. , 22657. , 13266. , 6990. , 3304. , 
     &     1310. , 386.5 , 174.5 , 64.0/
      DATA ecopt/261000. , 223000. , 172100. , 109000. , 158100. , 
     &     85840. , 121500. , 67500. , 107780. , 52960. , 85647. , 
     &     41429. , 63900. , 30294. , 43183. , 20271. , 32687. , 
     &     15351. , 21648. , 10193.7/
      DATA aco/6.48 , 91800. , 19700. , 6.49 , 8500. , 1630. , 5.5 , 
     &     3030. , 662. , 5.71 , 824. , 103. , 5.4 , 210. , 26.9 , 
     &     4.18 , 37. , 4.69 , 2.6 , 4.20 , .52 , 1.6 , .23 , .028 , 
     &     1.08 , .0342 , .0041 , .50 , .0026 , .00031/
      DATA oco/.018 , .015 , .1 , .036 , .028 , .191 , .048 , .041 , 
     &     .27 , .059 , .054 , .351 , .0773 , .0845 , .558 , .123 , 
     &     .128 , .58 , .60 , .24 , 1.92 , .42 , .36 , 2.39 , .42 , 
     &     .36 , 2.72 , 0. , 0. , 0./
      DATA ecinf/120000. , 66489. , 35911. , 17556. , 21835. , 10040. , 
     &     11800. , 5340. , 6460. , 2590. , 2939. , 1127. , 1112. , 
     &     414. , 306.8 , 113.4 , 131.3 , 49.1 , 43.5 , 16.4/
      DATA acinf/1600. , .862 , 4780. , 79.3 , .0057 , 94.7 , 21.5 , 
     &     .0008 , 17.9 , 4.74 , .000077 , 2.2 , .779 , 6.1E-6 , .311 , 
     &     .008 , 2.67E-7 , .0258 , .00459 , 5.1E-9 , .00128 , 9.76E-5 , 
     &     3.0E-11 , 2.6E-5 , 7.46E-6 , 1.16E-12 , 2.08E-6 , 2.7E-7 , 
     &     1.7E-14 , 7.9E-8/
      DATA ocinf/.05 , .009 , .02 , .089 , .019 , .037 , .12 , .026 , 
     &     .053 , .16 , .033 , .068 , .187 , .048 , .063 , .282 , .073 , 
     &     .095 , .57 , .12 , .24 , .95 , .212 , .378 , 1.13 , .28 , 
     &     .405 , .5 , .1 , .2/
      DATA an/1300. , 9000. , 500000. , 2000. , .65 , 770. , 32. , .88 , 
     &     180. , 710. , 2200. , 100. , .35 , 123. , 2.0 , .12 , 175. , 
     &     368. , 993. , 42. , 186. , 483. , 2.73 , 114. , 33. , 50. , 
     &     300. , 18. , .3 , 45. , 2.0 , .056 , 16. , 9. , 80. , 9.2 , 
     &     .18 , 21. , .07 , .023 , 3.8 , 2.4 , 13. , 2.5 , .15 , 5.3 , 
     &     .12 , .0054 , .437 , .40 , 1.27 , .393 , .115 , .52 , .0055 , 
     &     .00048 , .0614 , .117 , .0564 , .102 , .0615 , .0232 , 
     &     .00017 , 3.8E-5 , .0276 , .061 , .0066 , .053 , .035 , .027 , 
     &     2.02E-5 , 7.27E-6/
      DATA ohn/.12 , .035 , .040 , .025 , .035 , .24 , .06 , .15 , 
     &     .046 , .128 , .36 , .08 , .22 , .071 , .18 , .481 , .124 , 
     &     .282 , .094 , .231 , .529 , .193 , .468 , .161 , .298 , 
     &     .824 , .289 , .710 , .250 , .468 , 2.14 , .464 , 1.37 , .37 , 
     &     1.33 , 1.77 , .427 , 1.36 , .30 , 1.14 , 5*0./
      DATA efn/298000. , 160657. , 123174. , 156000. , 96899. , 82919. , 
     &     151405. , 97924. , 94589. , 123200. , 81699. , 76628. , 
     &     103600. , 67420. , 67140. , 82790. , 54171. , 54150. , 
     &     62153. , 40995. , 40950. , 40468. , 26830. , 26811. , 
     &     28840. , 19233. , 19224./
      DATA eosf/6.4E+6 , 5.6E+6 , 4.E+6 , 1.09E+6 , 8.05E+5 , 5.5E+5 , 
     &     314148. , 73768./
      DATA eoopt/290857. , 121294. , 167200. , 67660. , 150265. , 
     &     73668. , 122300. , 58000. , 99780. , 47000. , 77712. , 
     &     36348. , 55747. , 25841. , 33792. , 15868./
      DATA aoo/11. , 1.E+5 , 20000. , 7.4 , 8201. , 961.3 , 6.69 , 
     &     2450. , 78. , 6.9 , 710. , 79.2 , 5.5 , 148. , 16.2 , 4.23 , 
     &     21.4 , 2.39 , 2.7 , 2.0 , .223 , 1.22 , .0735 , .00845/
      DATA ooo/.023 , .025 , .12 , .039 , .036 , .238 , .055 , .056 , 
     &     .36 , .071 , .076 , .493 , .082 , .084 , .664 , .129 , .145 , 
     &     1.02 , .284 , .157 , 1.32 , 3*0./
      DATA eoinf/89316. , 73823. , 28660. , 24464. , 19645. , 16217. , 
     &     10630. , 7970. , 5570. , 4030. , 2519. , 1780. , 927. , 
     &     642.9 , 226.5 , 158.5/
      DATA aoi/220. , 16000. , .52 , 3.84 , .0111 , 317. , 69. , .0015 , 
     &     60. , .95 , .00022 , 11.4 , .194 , 1.82E-5 , 1.45 , .0217 , 
     &     1.01E-6 , .127 , .00115 , 2.18E-8 , .00597 , 1.74E-5 , 0. , 
     &     8.92E-5/
      DATA ooi/.029 , .01 , .054 , .045 , .025 , .114 , .068 , .035 , 
     &     .166 , .092 , .045 , .218 , .12 , .065 , .102 , .156 , .09 , 
     &     .40 , .19 , .13 , .53 , .016 , .0107 , .0231/
      DATA aff/2000. , 486. , 106. , 9. , 1.8 , .2 , .00859/
      DATA eff/102600. , 30028. , 18059. , 10130. , 5100. , 2226. , 
     &     782./
      DATA omff/.056 , .094 , .134 , .174 , .240 , .302 , .246/
      DATA aalf/190.8 , 60.05 , .716 , .0077 , .000214/
      DATA ealf/27762. , 18852. , 4305. , 950.2 , 287.32/
      DATA omalf/.17 , .19 , .51 , 1.6 , 6.0/
      DATA asio/5.14 , 901.7 , 165.1 , 4.4 , 33.5 , 3.82 , 3.29 , 6.60 , 
     &     .204 , 2.22 , .89 , .797 , 1.14 , .0031 , .030/
      DATA osiopt/.34 , .02 , .17 , .77 , .06 , .57 , 1.23 , 1.03 , 
     &     3.42 , 1.278 , 1.205 , 4.92 , 1. , 1. , 1./
      DATA esio/75220. , 43600. , 68588. , 34656. , 44999. , 18389. , 
     &     27163. , 11320. , 15394. , 6298./
      DATA esii/18561. , 9302.5 , 4087. , 1627. , 2029.2 , 763.9 , 
     &     832.5 , 297.2 , 223.16 , 77.11/
      DATA asii/9.55 , .007 , 13.74 , .199 , .000019 , .0771 , .027 , 
     &     1.24E-6 , .008 , .002 , 4.69E-8 , 4.32E-4 , 4.2E-5 , 0. , 
     &     8.75E-6/
      DATA osii/.26 , .07 , .04 , .65 , .21 , .14 , 1.04 , .32 , .26 , 
     &     2.7 , .56 , 1.2 , 1. , 1. , 1./
      DATA op/2.03 , .15 , .72 , .25 , .57 , 2.56 , .16 , .52 , .17 , 
     &     .41 , 5.65 , .635 , 1.42 , .6 , 1.35 , 12.1 , 2.31 , 6.67 , 
     &     2.33 , 7.11/
      DATA ep/78375. , 46083. , 41560. , 45125. , 27417. , 27000. , 
     &     35975. , 21219. , 21090. , 24540. , 14884.8 , 14853./
      DATA aso/12.52 , 1055. , 89.2 , 4.6 , 24.2 , 2.37 , 2.59 , 3.95 , 
     &     .396/
      DATA oso/.32 , .025 , .20 , .79 , .11 , .90 , .823 , .680 , 4.74/
      DATA eso/74473. , 31916. , 42681. , 17664. , 33266. , 12898./
      DATA asi/.36 , .012 , 34.08 , .0354 , .000045 , .311 , .031 , 
     &     2.4E-6 , .0052/
      DATA osi/.055 , .080 , .26 , .19 , .22 , .75 , 2.24 , .53 , 1.18/
      DATA esi/14440. , 12668. , 3276. , 2406. , 1570. , 458.1/
      DATA acl/69.15 , .543 , .0527/
      DATA ecl/15683. , 3115. , 1431.2/
      DATA ocl/.28 , 1.05 , .64/
      DATA ompal/.094 , .138 , .290 , .408 , .274/
      DATA ompb/.065 , .060 , .053 , .045 , .036 , .061 , .086 , .070 , 
     &     -.109 , -.225/
      DATA ompbsl/1.0E-8 , 1.2E-8 , 3.8E-8 , 8.7E-8 , 1.55E-7 , 3.2E-7 , 
     &     7.6E-7 , 2.40E-6 , 5.0E-6 , 1.28E-5/
      DATA ompcl/.081 , .379 , .50/      ! FIND PROTON RATE FOR AR
      DATA ompcls/3.96E-7 , 1.50E-7 , 4.57E-6/
      FIR(x) = 1. - x + x*x*EXP(x)*EXINT1(x,1)
      FCOol = 0.
      SQ = SQRT(T)
      DO 200 no = 2 , 10
C   GOES ONLY THROUGH CA
         n = NJ(no)
         lev = levno(no)
         NDUm = n
         realn = n
         ln = ALOG10(realn)
         nn = 12 - no
         FCT = 10.**(ABUnj(no)-5.0)*1.986/Dene
         jtop = n - 4
         DO 150 j = 1 , jtop
            JDUm = j
            i = n - j + 1
            IF ( i.GT.20 ) GOTO 150
            CON = CNC(no,j)
            IF ( CON.LT..0001 ) THEN
               IF ( i.EQ.1 ) THEN
               ELSEIF ( i.EQ.2 ) THEN
               ELSEIF ( i.EQ.3 ) THEN
               ELSEIF ( i.EQ.4 ) THEN
               ELSEIF ( i.EQ.5 .OR. i.EQ.9 .OR. i.EQ.13 .OR. i.EQ.17 )
     &                  THEN
                  GOTO 100
               ELSEIF ( i.EQ.6 .OR. i.EQ.7 .OR. i.EQ.14 .OR. 
     &                  i.EQ.15 .OR. i.EQ.16 ) THEN
                  GOTO 60
               ELSEIF ( i.EQ.8 ) THEN
                  FOUt(lev) = 0.
                  lev = lev + 1
                  GOTO 60
               ELSEIF ( i.EQ.10 ) THEN
               ELSEIF ( i.EQ.11 ) THEN
               ELSEIF ( i.EQ.12 ) THEN
               ELSEIF ( i.EQ.18 ) THEN
               ELSEIF ( i.EQ.19 ) THEN
                  GOTO 40
               ELSEIF ( i.EQ.20 ) THEN
                  GOTO 80
               ELSE
                  GOTO 20
               ENDIF
               GOTO 150
            ENDIF
 20         IF ( i.EQ.1 ) GOTO 150
            IF ( i.EQ.2 ) GOTO 150
            IF ( i.EQ.3 ) GOTO 150
            IF ( i.EQ.4 ) GOTO 150
            IF ( i.EQ.5 ) THEN
C   BE ISO INTERCOMBINATION LINE DONE WITH PERMITTED
               OHM = ombf(nn)
C  B ISO FORBIDDEN LINE; ] LINE DONE WITH PERMITTED
               IF ( n.EQ.6 ) OHM = OHM + .012*(T/100.)**.8/Elohy
               EN = eb(nn)
               A21 = acarb(nn)
               OHM = OHM + AMAX1(.83-.58*(ln-1),.83-15.7*(ln-1))
     &               *EXP((196000.-56000.*n)/T)
               ompeff = -ompb(nn) + ompbsl(nn)*T
               OHM = OHM + AMAX1(ompeff,0.)*(Dene-dne)/Dene
               CALL FSOLV(2,T,Dene,lev,1)
               GOTO 150
            ELSEIF ( i.EQ.6 ) THEN
               EU = ecopt(1,nn)
               EM = ecopt(2,nn)
               AUM = aco(1,nn)
               AUG = aco(2,nn)
               AMG = aco(3,nn)
               OHUm = oco(1,nn)
               OHUg = oco(2,nn)
               OHMg = oco(3,nn)
               IF ( n.EQ.6 ) THEN
                  t4 = AMIN1(.0001*T,3.)
                  OHMg = 1.156*t4*(1.085-.07507*t4-.0215*t4*t4)
                  OHUg = .2482*t4*(1.411-.4545*t4+.05848*t4*t4)
                  OHUm = .2764*SQRT(t4)*1.003
               ENDIF
               CALL FSOLV(3,T,Dene,lev,1)
               EU = ecinf(1,nn)
               EM = ecinf(2,nn)
               AUM = acinf(1,nn)
               AUG = acinf(2,nn)
               AMG = acinf(3,nn)
               OHUm = ocinf(1,nn)
               OHUg = ocinf(2,nn)
               OHMg = ocinf(3,nn)
               CALL FSOLV(3,T,Dene,lev,2)
               GOTO 150
            ELSEIF ( i.EQ.7 ) THEN
               A53 = an(1,nn)
               A52 = an(2,nn)
               A51 = an(3,nn)
               A43 = an(4,nn)
               A42 = an(5,nn)
               A41 = an(6,nn)
               A31 = an(7,nn)
               A21 = an(8,nn)
               OHUm = ohn(1,nn)
               OHUg = ohn(2,nn)
               OHMg = ohn(3,nn)
               OH45 = ohn(4,nn)
               OH23 = ohn(5,nn)
               IF ( n.EQ.7 ) THEN
                  t4 = .0001*T
                  OHUm = .82*t4**.59
                  OHMg = 2.53/(1.+4.23/t4)
                  OHUg = .802/(1.+3.72/t4)
               ENDIF
               EU = efn(1,nn)
               EN3 = efn(2,nn)
               EN2 = efn(3,nn)
               CALL FSOLV(5,T,Dene,lev,5)
               GOTO 150
            ELSEIF ( i.EQ.8 ) THEN
               A21 = 10.**(10.8-100./n)
C  OXYGEN ISO INTERCOMBINATION LINE ; O I SKELTON AND SHINE
               EN = eosf(nn)
               OHM = 444000.*SQRT(8./n)/EN
               IF ( n.EQ.8 ) OHM = 5.5E-5*T
               CALL FSOLV(2,T,Dene,lev,4)
               EU = eoopt(1,nn)
C   OXYGEN ISO OPTICAL LINES
               EM = eoopt(2,nn)
               AUM = aoo(1,nn)
               AUG = aoo(2,nn)
               AMG = aoo(3,nn)
               OHUm = ooo(1,nn)
               OHUg = ooo(2,nn)
               OHMg = ooo(3,nn)
               IF ( n.EQ.8 ) THEN
                  t4 = AMIN1(.0001*T,3.)
                  OHUg = .040*t4*(.97872+.05849*t4-.04724*t4*t4)
                  OHMg = .2676*t4*(1.0087+.01676*t4-.02952*t4*t4)
                  OHUm = .105*SQRT(t4)*(1.0483-.04759*t4)
               ENDIF
               CALL FSOLV(3,T,Dene,lev,1)
               EU = eoinf(1,nn)
               EM = eoinf(2,nn)
               AUM = aoi(1,nn)
               AUG = aoi(2,nn)
               AMG = aoi(3,nn)
               IF ( n.EQ.8 ) THEN
                  OHUm = .016*(T/1320.)*FIR(1320./T)
                  OHUg = .0107*(T/1320.)*FIR(1320./T)
     &                   + 4.32E-4*(T/1000.)**.8/Elohy
                  OHMg = .0231*(T/740.)*FIR(740./T) + 2.95E-5*(T/1000.)
     &                   **.8/Elohy
               ENDIF
               CALL FSOLV(3,T,Dene,lev,3)
               GOTO 150
            ELSEIF ( i.EQ.9 ) THEN
               A21 = aff(nn)
C  FLUORINE ISO FORBIDDEN LINE
               EN = eff(nn)
               OHM = omff(nn)
               OHM = OHM + AMIN1(-.11+5.24/n,1.13-9.6/n)
     &               *EXP((361000.-67300.*n)/T)
     &               + AMIN1(.24,-.315+11.3/n)
     &               *EXP(-10.**(7.939-24.4/n)/T)
               ompeff = -(.202-.92/n) + 51.55*T/realn**6.947
               OHM = OHM + AMAX1(ompeff,0.)*(Dene-dne)/Dene
               CALL FSOLV(2,T,Dene,lev,2)
               GOTO 150
            ELSEIF ( i.EQ.10 ) THEN
               GOTO 150
            ELSEIF ( i.EQ.11 ) THEN
               GOTO 150
            ELSEIF ( i.EQ.12 ) THEN
               GOTO 150
            ELSEIF ( i.EQ.13 ) THEN
               A21 = aalf(nn)
C   ALUMINUM ISO FORBIDDEN LINE
               EN = ealf(nn)
               OHM = omalf(nn)
               e3p2 = -394000. + 36000.*n
               e3d = -530000. + 46000.*n
               OHM = OHM + (1.173-4.5/n)*(2.72-.996/(1.+.126*T/e3p2))
     &               *EXP(-e3p2/T) + (.733+7.2/n)
     &               *(1.03-.404/(1.+.126*T/e3d))*EXP(-e3d/T)
               ompeff = -ompal(nn) + 56053.*T/realn**8.0673
               OHM = OHM + AMAX1(ompeff,0.)*(Dene-dne)/Dene
               CALL FSOLV(2,T,Dene,lev,1)
               GOTO 150
            ELSEIF ( i.EQ.14 ) THEN
               AUM = asio(1,nn)
               AUG = asio(2,nn)
               AMG = asio(3,nn)
               OHUm = osiopt(1,nn)
               OHUg = osiopt(2,nn)
               OHMg = osiopt(3,nn)
               EU = esio(1,nn)
               EM = esio(2,nn)
               CALL FSOLV(3,T,Dene,lev,1)
               AUM = asii(1,nn)
               AUG = asii(2,nn)
               AMG = asii(3,nn)
               OHUm = osii(1,nn)
               OHUg = osii(2,nn)
               OHMg = osii(3,nn)
               EU = esii(1,nn)
               EM = esii(2,nn)
               CALL FSOLV(3,T,Dene,lev,2)
               GOTO 150
            ELSEIF ( i.EQ.15 ) THEN
               A53 = ap(1,nn)
               A52 = ap(2,nn)
               A51 = ap(3,nn)
               A43 = ap(4,nn)
               A42 = ap(5,nn)
               A41 = ap(6,nn)
               A31 = ap(7,nn)
               A21 = ap(8,nn)
               OHUm = op(1,nn)
               OHUg = op(2,nn)
               OHMg = op(3,nn)
               OH45 = op(4,nn)
               OH23 = op(5,nn)
               EU = ep(1,nn)
               EN3 = ep(2,nn)
               EN2 = ep(3,nn)
               CALL FSOLV(5,T,Dene,lev,5)
               GOTO 150
            ELSEIF ( i.EQ.16 ) THEN
               IF ( n.EQ.16 ) GOTO 60
               AUM = aso(1,nn)
               AUG = aso(2,nn)
               AMG = aso(3,nn)
               OHUm = oso(1,nn)
               OHUg = oso(2,nn)
               OHMg = oso(3,nn)
               EU = eso(1,nn)
               EM = eso(2,nn)
               CALL FSOLV(3,T,Dene,lev,1)
               AUM = asi(1,nn)
               AUG = asi(2,nn)
               AMG = asi(3,nn)
               OHUm = osi(1,nn)
               OHUg = osi(2,nn)
               OHMg = osi(3,nn)
               EU = esi(1,nn)
               EM = esi(2,nn)
               CALL FSOLV(3,T,Dene,lev,3)
               GOTO 150
            ELSEIF ( i.EQ.17 ) THEN
               A21 = acl(nn)
               EN = ecl(nn)
               OHM = ocl(nn)
               OHM = OHM + (26/n)
     &               *(1.25*EXP((159000.-18800.*n)/T)+3.2*EXP
     &               ((1320000.-82300.*n)/T))
               ompeff = -ompcl(nn) + ompcls(nn)*T
               CALL FSOLV(2,T,Dene,lev,2)
               GOTO 150
            ELSEIF ( i.EQ.18 ) THEN
               GOTO 150
            ELSEIF ( i.NE.19 ) THEN
               GOTO 150
            ENDIF
 40         IF ( n.EQ.20 ) THEN
               EN = 13686.
               A21 = 1.01
               OHM = 2.3
               CALL FSOLV(2,T,Dene,lev,6)
            ENDIF
            GOTO 150
 60         FOUt(lev) = 0.
            lev = lev + 1
 80         IF ( n.NE.26 .AND. i.EQ.20 ) GOTO 150
            FOUt(lev) = 0.
            lev = lev + 1
            FOUt(lev) = 0.
            lev = lev + 1
            FOUt(lev) = 0.
            lev = lev + 1
 100        FOUt(lev) = 0.
            lev = lev + 1
            GOTO 150
            lev = lev + 2
            lev = lev + 2
 150     CONTINUE
 200  CONTINUE
      RETURN
      END
**==fprint.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      SUBROUTINE FPRINT(Num)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL ABIn , ABUnd , ABUnj , BIN , BINmin , BINsyz , CNC , FCOol , 
     &     FOUt , FTOt , FWV , PDUm , PTOt , WDUm
      INTEGER i , ibn , il , JDUm , JMX , k , LDUm , loc , n , NBIn , 
     &        NF , nhigh , NJ , nlow , no , Num
C*** End of declarations inserted by SPAG
C  PRINTS OUT FORBIDDEN LINE EMISSIVITIES
      INCLUDE 'rayspec.inc'
      COMMON /FLN   / FWV(320) , FTOt(320) , FOUt(320) , FCOol
      COMMON /LNPR  / JMX(11) , NF(11) , LDUm(5) , JDUm(5) , WDUm(5) , 
     &                PDUm(5)
      COMMON /COM   / CNC(12,30) , PTOt(12,220) , ABIn(MXBINS) , 
     &                BIN(MXBINS)
      COMMON /PARAMS/ NJ(12) , ABUnj(12) , ABUnd , BINmin , BINsyz , 
     &                NBIn
C
      DO 100 no = 2 , Num
         n = NJ(no)
         loc = 1
         nlow = NF(no-1) + 1
         nhigh = NF(no)
         DO 50 i = nlow , nhigh
            IF ( FWV(i)*FOUt(i).GT.1. ) THEN
               LDUm(loc) = i
               WDUm(loc) = FWV(i)
               PDUm(loc) = FOUt(i)
               IF ( loc.EQ.5 ) PRINT 99001 , n , 
     &                               ((LDUm(k),WDUm(k),PDUm(k)),k=1,5)
               ibn = (12399./FWV(i)-BINmin)/BINsyz
               IF ( ibn.LE.NBIn .AND. ibn.GE.1 ) BIN(ibn) = BIN(ibn)
     &              + FOUt(i)
               loc = MOD(loc,5) + 1
            ENDIF
 50      CONTINUE
         il = loc - 1
         IF ( loc.NE.1 ) PRINT 99001 , n , 
     &                         ((LDUm(k),WDUm(k),PDUm(k)),k=1,il)
 100  CONTINUE
      RETURN
99001 FORMAT (1X,I2,5(I6,F9.1,E10.3))
      END
**==fecool.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      FUNCTION FECOOL(Dene,T,Ipr)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL ABIn , ABUnd , ABUnj , BIN , BINmin , BINsyz , c2 , c3 , 
     &     CNC , Dene , fc , fct , FECOOL , FEForb , FEIii , FEIiio , 
     &     FEInf , FEPerm , PTOt , st
      REAL T
      INTEGER Ipr , NBIn , NJ
C*** End of declarations inserted by SPAG
C  COOLING DUE TO FE II AND FE III ASSUMING THAT CONCE=CONCE(MG)
C  IF IRON IS NOT EXPLICITLY CALCULATED.
C  NUSSBAUMER AND STOREY FE II OMEGA'S AND ROBB FE III
C
      INCLUDE 'rayspec.inc'
      COMMON /PARAMS/ NJ(12) , ABUnj(12) , ABUnd , BINmin , BINsyz , 
     &                NBIn
      COMMON /COM   / CNC(12,30) , PTOt(12,220) , ABIn(MXBINS) , 
     &                BIN(MXBINS)
      COMMON /FE    / FEInf , FEPerm , FEForb , FEIii , FEIiio
C
      st = SQRT(T)
      c2 = CNC(6,2)
      c3 = CNC(6,3)
      IF ( NJ(11).EQ.26 ) c2 = CNC(11,2)
      IF ( NJ(11).EQ.26 ) c3 = CNC(11,3)
      fct = 1.6E11*10.**(ABUnj(11)-12.)*8.63E-6*c2/st
      FEInf = fct*(.31*.048*EXP(-553./T)+.42*.99*EXP(-11400./T))
      FEInf = FEInf*580.*st/(Dene+580.*st)
      FEPerm = fct*3.4*4.47*EXP(-55000./T)
      FEForb = fct*0.1*2.9*EXP(-33000./T)*1.2E6*st/(Dene+1.2E6*st)
      fct = 1.6E11*10.**(ABUnj(11)-12.)*8.63E-6*c3/st
      FEIii = fct*.30*.054*EXP(-625./T)*1590.*st/(Dene+1590.*st)
      FEIiio = fct*.43*2.66*EXP(-30840./T)*521000.*st/(Dene+521000.*st)
C
      fc = FEInf + FEPerm + FEForb + FEIii + FEIiio
      FECOOL = fc
      IF ( Ipr.NE.0 ) PRINT 99001 , fc , FEInf , FEPerm , FEForb , 
     &                      FEIii , FEIiio
      RETURN
99001 FORMAT (
     &   ' IRON COOLING; FC,II IR, II UV, II OPT, III IR, III      OPT '
     &   ,6F9.6)
      END
**==etwo.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      FUNCTION ETWO(X)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL ETWO , EXINT1 , X
C*** End of declarations inserted by SPAG
      ETWO = EXP(-X) - X*EXINT1(X,1)
      RETURN
      END
