      implicit real*8(a-h,o-z)
      common/robcum/fro21,fro41,frs21,frs41,frr21,frr41
      dimension ratlin(3000,300),temp(3000),edens(300)
C
      open(13,file="inpline")
      open(6,file="OUTline")
C
C --- Which line ratios?
C --- IRAT = 1; [N II] 6548+6583 / [N II] 5755
C --- IRAT = 2; [O III] 4959+5007 / [O III] 4363
C --- IRAT = 3; [S II] 6716 / [S II] 6731
C --- IRAT = 4; [O II] 3729 / [O II] 3726
C --- IRAT = 5; [S II] 4069+4076 / [S II] 6716+6731
C
C -- Which input ratio? Give the number as RATIOIN
C
      read(13,*)IRAT,ratx
      RATIOIN=dble(ratx)
C
      abund=1.d0
      coltoti=0.d0
      coltot=0.d0
C -- dene = electron density in cm-3
C -- te = temperature in K
C -- We calculate temperatures between 1e3 K - 1e6 K,
C -- and densities between 1e2 - 1e4 cm-3.
C -- Steps are logarithmic in both density and temperature
C
C -- loop over density (in cm-3)
      do i=1,201
      denlog=2.d0+0.01d0*dble(i-1)
      dene=1.d1**denlog
      edens(i)=dene
C -- loop over temperatures (in K)
      do j=1,601
      telog=3.d0+0.005d0*dble(j-1)
      te=1.d1**telog
      temp(j)=te
      if(IRAT.eq.1) then
      CALL CFORB(TE,2,DENE,ABUND,COLTOTI,COLTOT,F21,F32,F4T123,
     &F5T23,F6T23,F54,COOL)
      ratlin(j,i)=f4t123/f54
      elseif(IRAT.eq.2) then
      CALL CFORB(TE,3,DENE,ABUND,COLTOTI,COLTOT,F21,F32,F4T123,
     &F5T23,F6T23,F54,COOL)
      ratlin(j,i)=f4t123/f54
      elseif(IRAT.eq.3) then
      CALL NFORB(TE,4,DENE,ABUND,COLTOTI,COLTOT,F23T1,F45T1,
     &F45T23,COOL)
      ratlin(j,i)=(1.d0-frr21)/frr21
      elseif(IRAT.eq.4) then
      CALL NFORB(TE,2,DENE,ABUND,COLTOTI,COLTOT,F23T1,F45T1,
     &F45T23,COOL)
      ratlin(j,i)=frr21/(1.d0-frr21)
      elseif(IRAT.eq.5) then
      CALL NFORB(TE,4,DENE,ABUND,COLTOTI,COLTOT,F23T1,F45T1,
     &F45T23,COOL)
      ratlin(j,i)=f45t1/f23t1
      endif
C       if(telog.eq.4.d0) then
C       write(6,100)edens(i),ratlin(j,i)
C       endif
      enddo
      enddo
      iskip=0
      do i=1,201
       do j=1,600
        if(iskip.eq.1) goto 5432
        if(ratlin(j,i).lt.RATIOIN.and.ratlin(j+1,i).gt.
     &          RATIOIN) then
         qqqq=(RATIOIN-ratlin(j,i))/(ratlin(j+1,i)-ratlin(j,i))
         tesave=temp(j)*(1.d0-qqqq)+temp(j+1)*qqqq
        elseif(ratlin(j,i).gt.RATIOIN.and.ratlin(j+1,i).lt.
     &          RATIOIN) then
         qqqq=(ratlin(j,i)-RATIOIN)/(ratlin(j,i)-ratlin(j+1,i))
         tesave=temp(j)*(1.d0-qqqq)+temp(j+1)*qqqq
        else
          continue
        endif
       enddo
       if(IRAT.eq.3.and.tesave.gt.8.d5) iskip=1
       if(tesave.le.1.d3) tesave=1.d3
 5432  continue
       write(6,100)edens(i),tesave,dlog10(edens(i)),
     &             dlog10(tesave)
      enddo
C      CALL KFORB(TE,1,DENE,ABUND,COLTOTI,COLTOT,F21,F31,F41,
C     &F51,F45T23,F6T45,COOL)
C      CALL BEFORB(TE,1,DENE,ABUND,COLTOTI,COLTOT,F31,F41,F51,
C     &F6T43,F65,COOL)
C      CALL NFORB(TE,4,DENE,ABUND,COLTOTI,COLTOT,F23T1,F45T1,
C     &F45T23,COOL)
C      call OFORB(TE,1,DENE,ABUND,COLTOTI,COLTOT,F21,F32,F4T123,
C     &F5T12,F54,COOL)
C      fffa1=frr21*f23t1
C      fffa2=(1.d0-frr21)*f23t1
C      fffa3=frr41*f45t1
C      fffa3=frr41*f45t23
C      fffa4=(1.d0-frr41)*f45t1
C      fffa4=(1.d0-frr41)*f45t23
C      write(6,100)dene,F4T123/F6T23
C      write(6,100)dene,frr21/(1.d0-frr21),fffa3/f23t1,fffa4/f23t1,
C     &f45t1/f23t1
C      write(6,100)denlog,dene,fffa2/fffa1,te
C     &,fffa3/fffa1,fffa4/fffa1,f45t1/fffa1
CCC      ratlin(i,ijji)=f4t123/f54
C -- write out the ( [O III] 4959+5007 / [O III 4363) - ratio
C      write(6,100)te,f4t123/f54
C      ratlin(i,ijji)=F23T1/F45T1
CCC      temp(i)=te
C      enddo
C      enddo
CCC      write(6,100)abund,(edens(ijji),ijji=1,36)
CCC      do i=1,2001
CCC        write(6,100)temp(i),(ratlin(i,ijji),ijji=1,36)
CCC      enddo
C      write(6,100)denlog,(1.d0-frr21)/frr21
C      enddo
C      enddo
  100 format(1x,40e13.5)
      stop
      end
c*********************************************************************   
c*********************************************************************   
      DOUBLE PRECISION FUNCTION ESCAPE(TAU)                             
C                                                                       
C     (Valid for Doppler profiles)                                      
C     Exact evaluation for TAU<1.E12                                    
C     Asymptotic expression for TAU>1.E12                               
C************************************************************           
C                                                                       
      IMPLICIT REAL*8(A-H,O-Z)                                          
C                                                                       
      COMMON/CONSTS/PI,PISQRT,CLIGHT                                    
C                                                                       
      DIMENSION X(11),F(11),G(11),T(11),XA(8),XB(8),TT(8)               
C                                                                       
      DATA XA/0.0D0,0.2D0,0.9D0,1.3D0,2.0D0,2.9D0,3.9D0,4.9D0/       
      DATA XB/3.3D0,3.4D0,3.5D0,3.6D0,3.9D0,4.5D0,5.0D0,6.0D0/          
      DATA TT/1.D0,1.3D0,1.6D0,2.D0,3.D0,5.D0,8.D0,12.D0/               
C                                                                       
      SUM=0.0D0                                                         
      if(tau.le.1.d-5) pesc=1.d0
      if(tau.le.1.d-5) goto 9999
      TTLOG=DLOG10(TAU)                                                 
      IF(TTLOG.LT.12.D0) GOTO 1000                                      
      PESC=1.D0/(2.D0*TAU*DSQRT(DLOG(TAU/PISQRT)))                      
      GOTO 9999                                                         
C                                                                       
 1000 CONTINUE                                                          
      IXLOOP=11                                                         
      ITLOOP=8                                                          
C                                                                       
      TTLOG=DLOG10(TAU)                                                 
      IF(TTLOG.GE.TT(1)) GOTO 40                                        
      XMIN=XA(1)                                                        
      XMAX=XB(1)                                                        
      GOTO 99                                                           
   40 CONTINUE                                                          
      IF(TTLOG.GE.TT(ITLOOP)) ISAVE=ITLOOP-1                            
      IF(TTLOG.GE.TT(ITLOOP)) GOTO 88                                   
      DO 50 I=1,ITLOOP-1                                                
      IF(TTLOG.LT.TT(I+1).AND.TTLOG.GE.TT(I)) GOTO 51                   
      GOTO 50                                                           
   51 ISAVE=I                                                           
      GOTO 88                                                           
   50 CONTINUE                                                          
   88 CC=(TTLOG-TT(ISAVE))/(TT(ISAVE+1)-TT(ISAVE))                      
      XMIN=XA(ISAVE)+(XA(ISAVE+1)-XA(ISAVE))*CC                         
      XMAX=XB(ISAVE)+(XB(ISAVE+1)-XB(ISAVE))*CC                         
   99 DX=(XMAX-XMIN)/DBLE(IXLOOP-1)                                     
C                                                                       
      DO 100 I=1,IXLOOP                                                 
      X(I)=XMIN+DX*DBLE(I-1)                                            
  100 CONTINUE                                                          
C                                                                       
      DO 200 I=1,IXLOOP                                                 
      F(I)=FIFUNC(X(I))                                                 
      G(I)=TAU*F(I)                                                     
  200 T(I)=F(I)*EXPINT(2,G(I))                                          
C                                                                       
      ISUM=(IXLOOP-1)/2                                                 
C                                                                       
      DO 350 I=1,ISUM                                                   
      SUM=SUM+(DX/3.D0)*(T(2*I-1)+T(2*I+1)+4.D0*T(2*I))                 
  350 CONTINUE                                                          
      PESC=2.D0*SUM                                                     
C                                                                       
 9999 ESCAPE=PESC                                                       
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
C************************************************************           
C************************************************************           
C                                                                       
      DOUBLE PRECISION FUNCTION FIFUNC(X)                               
C                                                                       
      IMPLICIT REAL*8(A-H,O-Z)                                          
C                                                                       
      COMMON/CONSTS/PI,PISQRT,CLIGHT                                    
C                                                                       
      FIFUNC=EXPFN(X**2)/PISQRT                                       
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
C************************************************************           
C************************************************************           
C                                                                       
      DOUBLE PRECISION FUNCTION EXPINT(N,X)                             
C                                                                       
C  COMPUTES THE N'TH EXPONENTIAL INTEGRAL OF X                          
C  INPUT - X,  INDEPENDENT VARIABLE (-100. .LE. X .LE. +100.)           
C         N,  ORDER OF DESIRED EXPONENTIAL INTEGRAL (1 .LE. N .LE. 8)   
C  OUTPUT - EXPINT,  THE DESIRED RESULT                                 
C           EX,  EXPF(-X)                                               
C  NOTE   RETURNS WITH E1(0)=0, (NOT INFINITY).                         
C  RUDOLF LOESER, 3 MAR 66                                              
C-----GENERAL COMPILATION OF 1 AUGUST 1967.                             
C                                                                       
      IMPLICIT REAL*8(A-H,O-Z)                                          
      DIMENSION TAB(20),XINT(7)                                         
      DATA XINT/1.,2.,3.,4.,5.,6.,7./                                   
      DATA TAB /.2707662555,.2131473101,.1746297218,.1477309984,        
     1.1280843565,.1131470205,.1014028126,.0919145454,.0840790292,      
     1.0774922515,.0718735405,.0670215610,.0627878642,.0590604044,      
     1.0557529077,.0527977953,.0501413386,.0477402600,.0455592945,      
     1.0435694088/                                                      
C                                                                       
      XSAVE=0.0D0                                                       
      U=X                                                               
      IF(U)603,602,603                                                  
  602 EX=1.                                                             
      IF(N-1)800,800,801                                                
  800 EXPINT=0.                                                         
      GOTO 777                                                          
  801 EXPINT=1./XINT(N-1)                                               
      GOTO 777                                                          
  603 IF(U-XSAVE)604,503,604                                            
  604 XSAVE=U                                                           
      XM=-U                                                             
      EMX=DEXP(XM)                                                      
C                                                                       
C  SELECT METHOD FOR COMPUTING EI(XM)                                   
C  COPIED FROM THE SHARE ROUTINE NUEXPI, WRITTEN BY J. W. COOLEY,       
C  COURANT INSTITUTE OF MATHEMATICAL SCIENCES, NEW YORK UNIVERSITY      
C                                                                       
      IF(XM-24.5)501,400,400                                            
  501 IF(XM-5.)502,300,300                                              
  502 IF(XM+1.)100,200,200                                              
  503 EISAVE=-ARG                                                       
      EXSAVE=EMX                                                        
C                                                                       
C  NOW RECURSE TO HIGHER ORDERS                                         
C                                                                       
      IF(N-1)507,507,505                                                
  505 DO 506 I=2,N                                                      
        EISAVE=(U*EISAVE-EXSAVE)/(-XINT(I-1))                           
  506 CONTINUE                                                          
  507 EXPINT=EISAVE                                                     
      EX=EXSAVE                                                         
  777 RETURN                                                            
C                                                                       
C  EI(XM) FOR XM .LT. -1.0                                              
C  HASTINGS POLYNOMIAL APPROXIMATION                                    
C                                                                       
  100 ARG=((((((U+8.573328740 )*U+18.05901697  )*U+8.634760893 )*U      
     *+.2677737343)/XM)*EMX)/((((U+9.573322345 )*U+25.63295615  )*U     
     *+21.09965308  )*U+3.958496923 )                                   
      GOTO 503                                                          
C     EI(XM) FOR -1. .LE. XM .LT. 5.0                                   
C     POWER SERIES EXPANSION ABOUT ZERO                                 
  200 ARG=DLOG(DABS(XM))                                                
      ARG=((((((((((((((((.41159050D-14*XM+.71745406D-13)*XM+.76404637D-
     *12)*XM+.11395905D-10)*XM+.17540077D-9)*XM+.23002666D-8)*XM+.275360
     *18D-7)*XM+.30588626D-6)*XM+.31003842D-5)*XM+.28346991D-4)*XM+.2314
     *8057D-3)*XM+.0016666574)*XM+.010416668)*XM+.055555572)*XM+.25)*XM+
     *.99999999)*XM+.57721566)+ARG                                      
      GOTO 503                                                          
C                                                                       
C  EI(XM) FOR 5.0 .LE. XM .LT. 24.5                                     
C  TABLE LOOK-UP AND INTERPOLATION                                      
C                                                                       
  300 I=XM+.5                                                           
      XZERO=I                                                           
      DELTA=XZERO-XM                                                    
      ARG=TAB(I-4)                                                      
      IF(DELTA)303,305,303                                              
  303 Y=ARG                                                             
      DELTAX=DELTA/XZERO                                                
      POWER=1./DELTAX                                                   
      DO 304 I=1,7                                                      
        POWER=POWER*DELTAX                                              
        Y=((Y-POWER/XZERO)*DELTA)/XINT(I)                               
        ARG=ARG+Y                                                       
        IF(DABS(Y/ARG)-1.D-8)305,304,304                                
  304 CONTINUE                                                          
  305 ARG=EMX*ARG                                                       
      GOTO 503                                                          
C     EI(XM) FOR 24.5 .LE. XM                                           
C     TRUNCATED CONTINUED FRACTION                                      
  400 ARG=((((XM-15.)*XM+58.)*XM-50.)*EMX)/((((XM-16.)*XM+72.)*XM-96.)  
     **XM+24.)                                                          
      GOTO 503                                                          
      END                                                               
C                                                                       
C************************************************************           
C************************************************************
      SUBROUTINE ESCACF(ENE,A21,ADAMP,G1,G2,TE,AMASS,TCORR,
     &COLTOTI,COLTOT,BE)
C********************************************************
C     Evaluate the escape probability (BE) for photons esca-
C     ping inwards and outwards. Take a mean value of these.
C     ENE in eV, AMASS in atomic mass units, TCORR correction
C     in opacity due to stimulated emission.
C********************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      parameter (md=105)
C
      if(tcorr.le.0.0d0) tcorr=0.0d0
      vterm=1.285d4*dsqrt(te/amass)
      wl=1.2398d4/ene
      dnud=1.d8*vterm/wl
      f12=1.499d-16*g2*a21*wl**2./g1
      cap=1.4974d-2*f12*tcorr/dnud
      t0=cap*coltoti
      t0t=cap*(coltot-coltoti)
      if(t0.le.0.0d0) then
       be1=1.d0
      else
       CALL ESCAP(T0,WL,ADAMP,VTERM,BE1)
      endif
      if(t0t.le.0.0d0) then
       be2=1.d0
      else
       CALL ESCAP(T0T,WL,ADAMP,VTERM,BE2)
      endif
      BE=0.5d0*(BE1+BE2)
C
      RETURN
      END
C************************************************************
C************************************************************
      SUBROUTINE ESCAP(T0,WL,ADAMP,VTERM,BE)
C********************************************************
C     Escape probability (from Ferland, coded by C.F.)
C     (Does not apply to lines trapped by the continuum.)
C********************************************************
      IMPLICIT REAL*8(A-H,O-Z)
C
      IF(VTERM.LE.0.) VTERM=1.285E6        
      A=7.958E-10*WL*ADAMP/VTERM
      AT=A*T0
      IF(AT.GT.1.) THEN
        B=1.6 + 3./((2.*A)**.12*(1.+AT))
      ELSEIF(AT.GT.0.) THEN
        B=1.6 + 3./((2.*A)**.12*(1.+1./SQRT(AT)))
      ELSE
        B=1.6
      ENDIF
      IF(B.GT.5) THEN
        B=5.
      ENDIF
      BE = 1./(1. + B * T0)
C
      RETURN
      END
C************************************************************           
c*********************************************************************   
      SUBROUTINE ESC(taul,tauc,cl,pesc,plin,pcont,dtau)
c*********************************************************************   
c     Calculates probabilities that a line photon escapes (pesc),
c     is self-absorbed (plin), or is trapped in the continuum (pcont).
c     For plin/pcont, it is assumed that tauc/taul=cl, i.e.
c     local conditions govern plin/pcont. This routine is good
c     for lines with Voigt parameter  a < 1E-3. Both absorption
c     and emission profiles are assumed to be Voigt profiles.
c     The routine also works for zero continuum absorption.
c*********************************************************************   
      implicit real*8(a-h,o-z)
c
      dimension acl(9),fcl(9),tl(13),gtl(13)
c
      data acl/-10.,-7.,-5.,-4.,-3.,-2.,-1.,0.,1./
      data fcl/2.064,1.310,.8797,.8122,.7372,.6385,.5172,.4265,.3997/
      data tl/-1.,-0.5,0.,0.5,1.,1.5,2.,3.,4.,5.,7.,9.,11./
      data gtl/-.01768,-.05425,-.1677,-.4882,-1.1204,-1.7481,-2.3201,
     &         -3.4121,-4.4541,-5.3275,-6.0770,-7.9987,-93.876/
c
      if(cl.lt.1.d-10) then
       q=1.159d2*cl
      elseif(cl.ge.1.d-10.and.cl.lt.1.d1) then
       cllog=dlog10(cl)
        do 100 icl=1,8
        if(cllog.ge.acl(icl).and.cllog.lt.acl(icl+1)) then
         y=(cllog-acl(icl))/(acl(icl+1)-acl(icl))
         fqlog=(1.d0-y)*fcl(icl)+y*fcl(icl+1)
         q=cl*1.d1**fqlog
        endif
  100   continue
      elseif(cl.ge.1.d1) then
       q=2.51d0*cl
      endif
c
      if(taul.lt.0.1d0) then
       g=1.d0-0.39894d0*taul
      elseif(taul.ge.0.1d0.and.taul.lt.1.d11) then
       tlog=dlog10(taul)
        do 101 itl=1,12
        if(tlog.ge.tl(itl).and.tlog.lt.tl(itl+1)) then
         x=(tlog-tl(itl))/(tl(itl+1)-tl(itl))
         fglog=(1.d0-x)*gtl(itl)+x*gtl(itl+1)
         g=1.d1**fglog
        endif
  101   continue
      elseif(taul.ge.1.d11) then
       g=1.33d-94
      endif
c
      if(dtau.lt.1.d-2) then
       pesc=g*expfn(tauc)
      else
       pesc=(g*(1.d0-expfn(dtau))*expfn(tauc))/dtau
      endif
      ptrap=1.d0-pesc
      plin=ptrap/(1.d0+q)
      pcont=(q*ptrap)/(1.d0+q)
c
      return
      end
c*********************************************************************   
C************************************************************           
      SUBROUTINE MATRIX(N)                                              
      IMPLICIT REAL*8(A-H,O-Z)                                          
C
      parameter (mmatr=80)
C
      COMMON/MSTUFF/B(mmatr),A(mmatr,mmatr),EPS,NERR
C
      DIMENSION LASTN(mmatr),ASAVE(mmatr)                                     
C                                                                       
C -- Find the column index of the last non-zero element in each row --  
C -- Speeds up the solution of loose matrices considerably --           
C                                                                       
      DO 100 I=1,N                                                      
      DO 101 J=N,I,-1                                                   
      JSAVE=J                                                           
      IF(A(I,J).NE.0.0D0) GOTO 102                                      
  101 CONTINUE                                                          
  102 LASTN(I)=JSAVE                                                    
  100 CONTINUE                                                          
C                                                                       
C -- ************************** --                                      
C -- Forward elimination scheme --                                      
C -- ************************** --                                      
C                                                                       
      DO 200 I=1,N-1                                                    
C                                                                       
C -- Partial pivoting routine --                                        
C                                                                       
      AI=0.0D0                                                          
      DO 201 L=I,N                                                      
      ASAVE(L)=0.0D0                                                    
      IF(DABS(A(L,I)).GT.AI) LSAVE=L                                    
      IF(DABS(A(L,I)).GT.AI) AI=DABS(A(L,I))                            
  201 CONTINUE                                                          
      IF(LSAVE.EQ.I) GOTO 202                                           
C                                                                       
C -- Interchange rows A(I,?) and A(LSAVE,?) if LSAVE.NE.I --            
C                                                                       
      LAST=MAX0(LASTN(I),LASTN(LSAVE))                                  
      DO 203 L=I,LAST                                                   
      ASAVE(L)=A(I,L)                                                   
      A(I,L)=A(LSAVE,L)                                                 
  203 A(LSAVE,L)=ASAVE(L)                                               
      BSAVE=B(I)                                                        
      B(I)=B(LSAVE)                                                     
      B(LSAVE)=BSAVE                                                    
      LASTNI=LASTN(I)                                                   
      LASTN(I)=LASTN(LSAVE)                                             
      LASTN(LSAVE)=LASTNI                                               
  202 CONTINUE                                                          
C                                                                       
C -- Elimination routine --                                             
C                                                                       
      BI=B(I)                                                           
      AII=A(I,I)                                                        
      DO 301 J=I+1,N                                                    
      IF(A(J,I).EQ.0.0D0.OR.(I+1).GT.LASTN(I)) GOTO 301                 
      AJIAII=A(J,I)/A(I,I)                                              
      DO 302 K=I+1,LASTN(I)                                             
  302 A(J,K)=A(J,K)-AJIAII*A(I,K)                                       
      B(J)=B(J)-AJIAII*BI                                               
      LASTN(J)=MAX0(LASTN(J),LASTN(I))                                  
  301 CONTINUE                                                          
C                                                                       
C -- ************************** --                                      
C -- End of forward elimination --                                      
C -- ************************** --                                      
C                                                                       
  200 CONTINUE                                                          
C                                                                       
C -- *************** --                                                 
C -- Back-substitute --                                                 
C -- *************** --                                                 
C                                                                       
      DO 400 K=N,1,-1                                                   
      BK=B(K)                                                           
      IF((K+1).GT.N) GOTO 401                                           
      DO 402 L=K+1,N                                                    
      BK=BK-A(K,L)*B(L)                                                 
  402 CONTINUE                                                          
  401 B(K)=BK/A(K,K)                                                    
  400 CONTINUE                                                          
C                                                                       
      RETURN                                                            
      END                                                               
C************************************************************           
C***************************************************                    
      DOUBLE PRECISION FUNCTION SIMUL(N,A,B,X)                          
      IMPLICIT REAL*8(A-H,O-Z)                                          
      parameter(nf=40,nfp1=41,nf2=80,nf21=81)                           
      dimension a(nf21,nf21),b(nf21),x(nf21)                            
      DIMENSION LASTN(nf21),ASAVE(nf21),BS(nf21)                        
C                                                                       
C -- Find the column index of the last non-zero element in each row --  
C -- Speeds up the solution of loose matrices considerably --           
C                                                                       
      DO 100 I=1,N                                                      
      BS(I)=B(I)                                                        
      DO 101 J=N,I,-1                                                   
      JSAVE=J                                                           
      IF(A(I,J).NE.0.0D0) GOTO 102                                      
  101 CONTINUE                                                          
  102 LASTN(I)=JSAVE                                                    
  100 CONTINUE                                                          
C                                                                       
C -- ************************** --                                      
C -- Forward elimination scheme --                                      
C -- ************************** --                                      
C                                                                       
      DO 200 I=1,N-1                                                    
C                                                                       
C -- Partial pivoting routine --                                        
C                                                                       
      AI=0.0D0                                                          
      DO 201 L=I,N                                                      
      ASAVE(L)=0.0D0                                                    
      IF(DABS(A(L,I)).GT.AI) LSAVE=L                                    
      IF(DABS(A(L,I)).GT.AI) AI=DABS(A(L,I))                            
  201 CONTINUE                                                          
      IF(LSAVE.EQ.I) GOTO 202                                           
C                                                                       
C -- Interchange rows A(I,?) and A(LSAVE,?) if LSAVE.NE.I --            
C                                                                       
      LAST=MAX0(LASTN(I),LASTN(LSAVE))                                  
      DO 203 L=I,LAST                                                   
      ASAVE(L)=A(I,L)                                                   
      A(I,L)=A(LSAVE,L)                                                 
  203 A(LSAVE,L)=ASAVE(L)                                               
      BSAVE=B(I)                                                        
      B(I)=B(LSAVE)                                                     
      B(LSAVE)=BSAVE                                                    
      LASTNI=LASTN(I)                                                   
      LASTN(I)=LASTN(LSAVE)                                             
      LASTN(LSAVE)=LASTNI                                               
  202 CONTINUE                                                          
C                                                                       
C -- Elimination routine --                                             
C                                                                       
      BI=B(I)                                                           
      AII=A(I,I)                                                        
      DO 301 J=I+1,N                                                    
      IF(A(J,I).EQ.0.0D0.OR.(I+1).GT.LASTN(I)) GOTO 301                 
      AJIAII=A(J,I)/A(I,I)                                              
      DO 302 K=I+1,LASTN(I)                                             
  302 A(J,K)=A(J,K)-AJIAII*A(I,K)                                       
      B(J)=B(J)-AJIAII*BI                                               
      LASTN(J)=MAX0(LASTN(J),LASTN(I))                                  
  301 CONTINUE                                                          
C                                                                       
C -- ************************** --                                      
C -- End of forward elimination --                                      
C -- ************************** --                                      
C                                                                       
  200 CONTINUE                                                          
C                                                                       
C -- *************** --                                                 
C -- Back-substitute --                                                 
C -- *************** --                                                 
C                                                                       
      DO 400 K=N,1,-1                                                   
      BK=B(K)                                                           
      IF((K+1).GT.N) GOTO 401                                           
      DO 402 L=K+1,N                                                    
      BK=BK-A(K,L)*B(L)                                                 
  402 CONTINUE                                                          
  401 B(K)=BK/A(K,K)                                                    
  400 CONTINUE                                                          
      DO 9 I=1,N                                                        
      X(I)=B(I)                                                         
9     B(I)=BS(I)                                                        
C                                                                       
      simul=0.                                                          
      RETURN                                                            
      END                                                               
C************************************************************           
C************************************************************
      double precision function expfn(x)
C
      IMPLICIT REAL*8(A-H,O-Z)
      if(x.gt.2.d2) then
       expfn=1.384d-87
      elseif(dabs(x).lt.1.d-15) then
       expfn=1.d0
      else
       expfn=dexp(-x)
      endif
      return
      end
C************************************************************
C*********************************************************************
      SUBROUTINE BEFORB(TE,I,DENE,ABUND,COLTOTI,COLTOT,F31,F41,F51,
     &F6T43,F65,COOL)
C*********************************************************************
C  Forbidden lines (Be I-like). Enumeration of the ions are: (= I)
C
C  1 = C III  2 = Mg I  3 = Al II  4 = Si III
C
C  The atom is treated as a 6-level atom.
C  (Returned to RAD when E21/Tev > 20).
C
C  A-values for C III from Allard et al. (AAS 84, 563, '90) (superseded
C  by Verner, Verner & Ferland '96), 
C  Fleming et al. (Phys.Scr. 49, 316, '94; MNRAS, 279, 1289) and 
C  Nussbaumer & Storey (A&A 64, 139, '78). 
C  Energies from Mendoza ('83) and Daresbury ('86).
C  Omegas from Seaton (J.Phys.B. 20, 6431, '87), and Keenan et al.
C  (ApJ 389, 443, '92).
C  A-values for Al II from Johnson et al. (ApJ 308, 1013, '86) 
C  and Aggarwal & Keenan (J.Phys.B 27, 2343, '94). Energies from 
C  the ADS NIST Database ('95). Omegas from Keenan et al. (ApJ 385,
C  375, '92) and Aggarwal & Keenan (J.Phys.B 27, 5321, '94). 
C  Energies for Mg I and Si III from the ADS NIST Database ('95). 
C  (Note that for Mg I there are 4s levels (not included) intervening 
C  the model atom below. Note also that the 3PE level is not included 
C  for C III. This lies below 1D, which is included.)
C*********************************************************************
      implicit real*8(a-h,o-z)
C
      parameter (mmatr=80)
C
      COMMON/MSTUFF/X(mmatr),A(mmatr,mmatr),EPS,NERR
c
      dimension ene2(4),ene3(4),ene4(4),ene5(4),ene6(4),
     &a31(4),a32(4),a41(4),a42(4),a43(4),a51(4),a63(4),
     &a64(4),a65(4),om1t234(4),om15(4),om16(4),om23(4),om24(4),
     &om34(4),om234t5(4),om234t6(4),om56(4),al1t234(4),al15(4),
     &al16(4),al23(4),al24(4),al34(4),al234t5(4),al234t6(4),
     &al56(4),g(6),ama(4)
c
      data ene2/6.4924,2.7092,4.6364,6.5373/
      data ene3/6.4954,2.7117,4.6439,6.5533/
      data ene4/6.5024,2.7168,4.6593,6.5857/
      data ene5/12.690,4.3460,7.4210,10.2768/
      data ene6/18.090,5.7535,10.5988,15.1534/
      data a31/1.038d2,0.,3.33d3,0./
      data a32/2.39d-7,1.45d-7,3.93d-6,3.82d-5/
      data a41/5.19d-3,0.,3.33d-3,0./
      data a42/4*0./
      data a43/2.41d-6,9.10d-7,2.57d-5,2.42d-4/
      data a51/1.79d9,0.,1.44d9,0./
      data a63/4.85d3,0.,7.75d3,0./
      data a64/2.97d4,0.,2.03d4,0./
      data a65/1.40d8,0.,2.63d8,0./
      data om1t234/1.095d0,0.,3.167d0,0./
      data om15/4.35d0,0.,7.382d0,0./
      data om16/0.415d0,0.,1.326d0,0./
      data om23/0.95d0,0.,1.67d0,0./
      data om24/0.695d0,0.,2.00d0,0./
      data om34/2.75d0,0.,6.54d0,0./
      data om234t5/3.408d0,0.,7.635d0,0./
      data om234t6/1.425d0,0.,11.63d0,0./
      data om56/10.48d0,0.,5.324d0,0./
      data al1t234/-0.0734d0,0.,0.189d0,0./
      data al15/0.063d0,0.,0.0835d0,0./
      data al16/-0.0419d0,0.,-0.126d0,0./
      data al23/0.18d0,0.,0.039d0,0./
      data al24/0.27d0,0.,0.076d0,0./
      data al34/0.27d0,0.,0.0656d0,0./
      data al234t5/-0.163d0,0.,-0.186d0,0./
      data al234t6/-0.0665d0,0.,-0.364d0,0./
      data al56/0.0490d0,0.,-0.27d0,0./
      data g/1.,1.,3.,5.,3.,5./
      data ama/12.,24.305,26.982,28.086/
C
      cool=0.0d0
      f31=0.0d0
      f41=0.0d0
      f51=0.0d0
      f6t43=0.0d0
      f65=0.0d0
      if(abund.ge.1.d-10) then
c
      tev=te/1.1609d4
      t4=te/1.d4
      cons=8.63d-6*dene/dsqrt(te)
c
      e21=ene2(i)
        IF((E21/TEV).lt.2.d1) THEN
      e31=ene3(i)
      e41=ene4(i)
      e51=ene5(i)
      e61=ene6(i)
      e32=ene3(i)-ene2(i)
      e42=ene4(i)-ene2(i)
      e52=ene5(i)-ene2(i)
      e62=ene6(i)-ene2(i)
      e43=ene4(i)-ene3(i)
      e53=ene5(i)-ene3(i)
      e63=ene6(i)-ene3(i)
      e54=ene5(i)-ene4(i)
      e64=ene6(i)-ene4(i)
      e65=ene6(i)-ene5(i)
C
      o1t234=om1t234(i)*t4**al1t234(i)
      o15=om15(i)*t4**al15(i)
      o16=om16(i)*t4**al16(i)
      o23=om23(i)*t4**al23(i)
      o24=om24(i)*t4**al24(i)
      o34=om34(i)*t4**al34(i)
      o234t5=om234t5(i)*t4**al234t5(i)
      o234t6=om234t6(i)*t4**al234t6(i)
      o56=om56(i)*t4**al56(i)
      if(i.eq.3) then
	if(te.gt.1.5d4) then
          o1t234=3.419d0*(t4/1.5d0)**(-.169)
	endif
	if(te.gt.4.d4) then
          o16=1.113d0*(t4/4.d0)**.185
	endif
      endif
      gsum=g(2)+g(3)+g(4)
      o12=(g(2)*o1t234)/gsum
      o13=(g(3)*o1t234)/gsum
      o14=(g(4)*o1t234)/gsum
      o25=(g(2)*o234t5)/gsum
      o26=(g(2)*o234t6)/gsum
      o35=(g(3)*o234t5)/gsum
      o36=(g(3)*o234t6)/gsum
      o45=(g(4)*o234t5)/gsum
      o46=(g(4)*o234t6)/gsum
      c21=cons*o12/g(2)
      c31=cons*o13/g(3)
      c41=cons*o14/g(4)
      c51=cons*o15/g(5)
      c61=cons*o16/g(6)
      c32=cons*o23/g(3)
      c42=cons*o24/g(4)
      c52=cons*o25/g(5)
      c62=cons*o26/g(6)
      c43=cons*o34/g(4)
      c53=cons*o35/g(5)
      c63=cons*o36/g(6)
      c54=cons*o45/g(5)
      c64=cons*o46/g(6)
      c65=cons*o56/g(6)
      c12=c21*expfn(e21/tev)*g(2)/g(1)
      c13=c31*expfn(e31/tev)*g(3)/g(1)
      c14=c41*expfn(e41/tev)*g(4)/g(1)
      c15=c51*expfn(e51/tev)*g(5)/g(1)
      c16=c61*expfn(e61/tev)*g(6)/g(1)
      c23=c32*expfn(e32/tev)*g(3)/g(2)
      c24=c42*expfn(e42/tev)*g(4)/g(2)
      c25=c52*expfn(e52/tev)*g(5)/g(2)
      c26=c62*expfn(e62/tev)*g(6)/g(2)
      c34=c43*expfn(e43/tev)*g(4)/g(3)
      c35=c53*expfn(e53/tev)*g(5)/g(3)
      c36=c63*expfn(e63/tev)*g(6)/g(3)
      c45=c54*expfn(e54/tev)*g(5)/g(4)
      c46=c64*expfn(e64/tev)*g(6)/g(4)
      c56=c65*expfn(e65/tev)*g(6)/g(5)
c
      amass=ama(i)
      m=0
 9988 continue
      m=m+1
      adamp=a31(i)+a32(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(3)*g(1))/(x(1)*g(3)))
      call escacf(e31,a31(i),adamp,g(1),g(3),te,amass,tcorr,
     &coltoti,coltot,p31)
      adamp=a41(i)+a43(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(4)*g(1))/(x(1)*g(4)))
      call escacf(e41,a41(i),adamp,g(1),g(4),te,amass,tcorr,
     &coltoti,coltot,p41)
      adamp=a51(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(5)*g(1))/(x(1)*g(5)))
      call escacf(e51,a51(i),adamp,g(1),g(5),te,amass,tcorr,
     &coltoti,coltot,p51)
      adamp=a31(i)+a32(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(3)*g(2))/(x(2)*g(3)))
      call escacf(e32,a32(i),adamp,g(2),g(3),te,amass,tcorr,
     &coltoti,coltot,p32)
      adamp=a41(i)+a43(i)+a32(i)+a31(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(3)*(1.d0-(x(4)*g(3))/(x(3)*g(4)))
      call escacf(e43,a43(i),adamp,g(3),g(4),te,amass,tcorr,
     &coltoti,coltot,p43)
      adamp=a63(i)+a64(i)+a65(i)+a31(i)+a32(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(3)*(1.d0-(x(6)*g(3))/(x(3)*g(6)))
      call escacf(e63,a63(i),adamp,g(3),g(6),te,amass,tcorr,
     &coltoti,coltot,p63)
      adamp=a63(i)+a64(i)+a65(i)+a41(i)+a43(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(4)*(1.d0-(x(6)*g(4))/(x(4)*g(6)))
      call escacf(e64,a64(i),adamp,g(4),g(6),te,amass,tcorr,
     &coltoti,coltot,p64)
      adamp=a63(i)+a64(i)+a65(i)+a51(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(5)*(1.d0-(x(6)*g(5))/(x(5)*g(6)))
      call escacf(e65,a65(i),adamp,g(5),g(6),te,amass,tcorr,
     &coltoti,coltot,p65)
      a(1,1)=c12+c13+c14+c15+c16
      a(1,2)=-c21
      a(1,3)=-(c31+p31*a31(i))
      a(1,4)=-(c41+p41*a41(i))
      a(1,5)=-(c51+p51*a51(i))
      a(1,6)=-c61
      a(2,1)=-c12
      a(2,2)=c21+c23+c24+c25+c26
      a(2,3)=-(c32+p32*a32(i))
      a(2,4)=-c42
      a(2,5)=-c52
      a(2,6)=-c62
      do 100 j=1,6
      x(j)=0.d0
      a(3,j)=1.d0
  100 continue
      x(3)=1.d0
      a(4,1)=-c14
      a(4,2)=-c24
      a(4,3)=-c34
      a(4,4)=c41+c42+c43+c45+c46+p41*a41(i)+p43*a43(i)
      a(4,5)=-c54
      a(4,6)=-(c64+p64*a64(i))
      a(5,1)=-c15
      a(5,2)=-c25
      a(5,3)=-c35
      a(5,4)=-c45
      a(5,5)=c51+c52+c53+c54+c56+p51*a51(i)
      a(5,6)=-(c65+p65*a65(i))
      a(6,1)=-c16
      a(6,2)=-c26
      a(6,3)=-c36
      a(6,4)=-c46
      a(6,5)=-c56
      a(6,6)=c61+c62+c63+c64+c64+p63*a63(i)+p64*a64(i)+p65*a65(i)
c
      call matrix(6)
      if(m.lt.2) goto 9988
c
      elab=(1.6022d-12*abund)/(12.57d0*dene)
      cool=elab*(e21*(x(1)*c12-x(2)*c21)+e31*
     &(x(1)*c13-x(3)*c31)+e41*(x(1)*c14-x(4)*c41)+e51*(x(1)*c15-
     &x(5)*c51)+e61*(x(1)*c16-x(6)*c61)+e32*(x(2)*c23-x(3)*c32)+
     &e42*(x(2)*c24-x(4)*c42)+e52*(x(2)*c25-x(5)*c52)+e62*(x(2)*
     &c26-x(6)*c62)+e43*(x(3)*c34-x(4)*c43)+e53*(x(3)*c35-x(5)*
     &c53)+e63*(x(3)*c36-x(6)*c63)+e54*(x(4)*c45-x(5)*c54)+
     &e64*(x(4)*c46-x(6)*c64)+e65*(x(5)*c56-x(6)*c65))
      cool=12.57d0*cool
      f31=elab*x(3)*a31(i)*e31*p31
      f41=elab*x(4)*a41(i)*e41*p41
      f51=elab*x(5)*a51(i)*e51*p51
      f6t43=elab*x(6)*(p64*e64*a64(i)+p63*e63*a63(i))
      f65=elab*x(6)*e65*a65(i)*p65
c
        ENDIF
      endif
C
      RETURN
      END
C*********************************************************************
C*********************************************************************
      SUBROUTINE BFORB(TE,I,DENE,DENH,ABUND,COLTOTI,COLTOT,F21,F345T12,
     &F43,F54,COOL)
C*********************************************************************
C  Forbidden lines (B I-like). Enumeration of the ions are: (= I)
C
C  1 = C II  2 = N III  3 = O IV  4 = S IV  5 = Si II  6 = Ne VI
C  (Ar VI is mainly in the EUV/UV.)
C
C  The atom is treated as a 5-level atom.
C  (Reduced to a 2-level atom for E32/Tev > 20)
C
C  A-values not mentioned specifically are from Mendoza ('83).
C  A-values for C II from Biemont et al. (J. Phys.B, 27, 5841, '94.) 
C  A-values for N III from Brage et al. (ApJ, 445, 457, '95). 
C  (These A-values are somewhat lower than those of Biemont et al, 
C  J.Phys.B, '94.) A-values for O IV from Pradhan's compilation 
C  (Jan. '96; Merkelis et al. '94.), updated with Brage, Judge & Brekke
C  (ApJ, 464, 1030, '96). Si II from Calamai et al, 
C  (ApJ, 415, L59, '93). A43 and A54 are isoelectronically guessed for 
C  N III. For S IV, Si II and Ne VI they are put equal to the values for 
C  O IV. A21 and om21 for Si II from Kafatos & Lynch, ApJS 42, 611, 1980. 
C  A21 of Ne VI from Mendoza ('83) and for intercombination lines from
C  Dankwort & Trefftz (A & A, 65, 93, '78).
C  For C and O we also use Opacity Project data (Blum & Pradhan, ApJS 
C  80, 425, 1992). For O IV and Ne VI we use the collisional strengths
C  of Zhang et al. (A & A, 283, 319, '93). N III is updated with Stafford, 
C  Bell, & Hibbert, MN 266, 715 (1994). om34, om35 and om45 isoelectronically 
C  guessed for Si at 1E4 K (use S IV's T-dependence for Si II).
C  Collisions with H for fine structure lines of C II and Si II are
C  included (E.Roueff, A&A 234, 567, 1990).
C  Energy levels of C II, N III & O IV slightly changed by Peng & Pradhan
C  (AAS, 112, 151, '95.)
C*********************************************************************
      implicit real*8(a-h,o-z)
C
      parameter (mmatr=80)
C
      COMMON/MSTUFF/X(mmatr),A(mmatr,mmatr),EPS,NERR
c
      dimension ene2(6),ene3(6),ene4(6),ene5(6),a21(6),a31(6),
     &a41(6),a32(6),a42(6),a52(6),a43(6),a54(6),om12(6),om13(6),
     &om14(6),om15(6),om23(6),om24(6),om25(6),om34(6),om35(6),
     &om45(6),al12(6),al13(6),al14(6),al15(6),al23(6),al24(6),
     &al25(6),al34(6),al35(6),al45(6),g(5),ama(6)
c
      data ene2/.007932,.02163,.04785,.1177,.03558,.16182/
      data ene3/5.3283,7.0869,8.8574,8.7861,5.3095,12.4231/
      data ene4/5.3310,7.0943,8.8736,8.8361,5.3230,12.4775/
      data ene5/5.3346,7.1044,8.8966,8.8590,5.3447,12.5578/
      data a21/2.29d-6,4.77d-5,5.20d-4,7.73d-3,2.17d-4,2.02d-2/
      data a31/7.27d1,3.61d2,1.47d3,5.50d4,5.20d3,1.27d4/
      data a41/1.88,9.11,38.4,140.,13.2,301./
      data a32/82.8,372.,1.43d3,3.39d4,4.41d3,1.11d4/
      data a42/10.3,65.1,294.,1.95d4,1.22d3,2.80d3/
      data a52/50.9,282.,1.17d3,3.95d4,2.46d3,9.87d3/
      data a43/2.39d-7,3.5d-6,4*5.07d-5/
      data a54/3.67d-7,6.1d-6,4*1.02d-4/
      data om12/0.,1.44,2.60,6.42,3.78,1.48/
      data om13/.242,.198,.139,.51,.406,.058/
      data om14/.362,.298,.217,.87,.692,.09/
      data om15/.235,.201,.169,.95,.752,.074/
      data om23/.177,.151,.123,.66,.545,.054/
      data om24/.477,.399,.309,1.47,1.17,.131/
      data om25/1.02,.844,.619,2.53,2.01,.259/
      data om34/.824,1.10,1.20,3.04,2.00,.587/
      data om35/.853,.668,.822,2.53,1.85,.401/
      data om45/1.98,2.04,2.35,5.85,4.50,.586/
      data al12/0.,.182,-.398,2*0.,-.569/
      data al13/.029,.063,-.309,2*-.16,-.284/
      data al14/.035,.083,-.328,2*-.133,-.409/
      data al15/.062,.185,.091,2*-.121,-.509/
      data al23/.057,.16,-.377,2*-.122,-.466/
      data al24/.046,.126,-.357,2*-.134,-.451/
      data al25/.036,.089,-.33,2*-.123,-.396/
      data al34/.365,.072,-.557,2*-.174,-.434/
      data al35/.186,.08,-.275,2*-.141,-.399/
      data al45/.269,.08,-.287,2*-.197,-.419/
      data g/2.,4.,2.,4.,6./
      data ama/12.,14.,16.,32.1,28.09,20.179/
C
      cool=0.0d0
      f21=0.0d0
      f345t12=0.0d0
      f43=0.0d0
      f54=0.0d0
      if(abund.ge.1.d-10) then
c
      tmx=2.d5
      tev=te/1.1609d4
      t4=te/1.d4
      cons=8.63d-6*dene/dsqrt(te)
c
      e21=ene2(i)
      e31=ene3(i)
      e41=ene4(i)
      e51=ene5(i)
      e32=ene3(i)-ene2(i)
      e42=ene4(i)-ene2(i)
      e52=ene5(i)-ene2(i)
      e43=ene4(i)-ene3(i)
      e53=ene5(i)-ene3(i)
      e54=ene5(i)-ene4(i)
c
      if(i.eq.1) then
       o12=dmax1(1.58d0,1.58d0*(te/1.d3)**0.11d0)
      elseif(i.eq.2) then
       o12=cheb(tmx,te,0.724d0,0.2349d0,-0.0682d0,-0.0517d0)
      elseif(i.eq.3) then
        if(t4.gt.0.6d0) then
           o12=2.70d0*(te/6.d3)**(-.25)
        else
           o12=2.70d0*(te/6.d3)**(.13)
        endif
      elseif(i.ge.4) then
        o12=om12(i)*t4**al12(i)
      endif
      c21=cons*o12/g(2)
C - H-atom collisions for C II and Si II
      if(i.eq.1) c21=c21+2.05d-9*denh*(te/2.d3)**.0759
      if(i.eq.5) c21=c21+1.1d-9*denh*(te/2.d3)**.345
      c12=c21*expfn(e21/tev)*g(2)/g(1)
       IF((E32/TEV).lt.2.d1) then
      if(i.ne.2.or.i.ne.3) then
       o13=om13(i)*t4**al13(i)
       o14=om14(i)*t4**al14(i)
       o15=om15(i)*t4**al15(i)
       o23=om23(i)*t4**al23(i)
       o24=om24(i)*t4**al24(i)
       o25=om25(i)*t4**al25(i)
       o34=om34(i)*t4**al34(i)
       o35=om35(i)*t4**al35(i)
       o45=om45(i)*t4**al45(i)
      elseif(i.eq.2) then
       o13=cheb(tmx,te,-3.6087d0,-0.0656d0,-0.0831d0,-0.0165d0)
       o14=cheb(tmx,te,-2.7917d0,-0.0501d0,-0.0914d0,-0.0197d0)
       o15=cheb(tmx,te,-3.5774d0,0.0274d0,-0.1330d0,-0.0337d0)
       o23=cheb(tmx,te,-4.1572d0,0.0119d0,-0.1247d0,-0.0311d0)
       o24=cheb(tmx,te,-2.2053d0,-0.0191d0,-0.1080d0,-0.0257d0)
       o25=cheb(tmx,te,-0.7075d0,-0.0465d0,-0.0934d0,-0.0204d0)
       o34=cheb(tmx,te,-0.0975d0,0.0837d0,-0.0437d0,-0.0195d0)
       o35=cheb(tmx,te,-0.7808d0,0.1475d0,-0.0418d0,-0.0211d0)
       o45=cheb(tmx,te,1.2787d0,0.1128d0,-0.0426d0,-0.0203d0)
      elseif(i.eq.3) then
        if(t4.gt.0.4d0) then
           o13=0.154d0*(te/4.d3)**(-.24)
        else
           o13=0.154d0*(te/4.d3)**(.12)
        endif
        if(t4.gt.0.4d0) then
           o14=0.24d0*(te/4.d3)**(-.25)
        else
           o14=0.24d0*(te/4.d3)**(.14)
        endif
        if(t4.gt.0.5d0) then
           o15=0.189d0*(te/5.d3)**(-.32)
        else
           o15=0.189d0*(te/5.d3)**(.22)
        endif
        if(t4.gt.0.4d0) then
           o23=0.137d0*(te/4.d3)**(-.28)
        else
           o23=0.137d0*(te/4.d3)**(.23)
        endif
        if(t4.gt.0.4d0) then
           o24=0.342d0*(te/4.d3)**(-.27)
        else
           o24=0.342d0*(te/4.d3)**(.18)
        endif
        if(t4.gt.0.4d0) then
           o25=0.189d0*(te/4.d3)**(-.25)
        else
           o25=0.189d0*(te/4.d3)**(.15)
        endif
        if(t4.gt.0.7d0) then
           o34=1.23d0*(te/7.d3)**(-.25)
        else
           o34=1.23d0*(te/7.d3)**(.064)
        endif
        if(t4.gt.0.7d0) then
           o35=0.836d0*(te/7.d3)**(-.23)
        else
           o35=0.836d0*(te/7.d3)**(.055)
        endif
        if(t4.gt.0.7d0) then
           o45=2.40d0*(te/7.d3)**(-.25)
        else
           o45=2.40d0*(te/7.d3)**(.060)
        endif
      endif
      c31=cons*o13/g(3)
      c41=cons*o14/g(4)
      c51=cons*o15/g(5)
      c32=cons*o23/g(3)
      c42=cons*o24/g(4)
      c52=cons*o25/g(5)
      c43=cons*o34/g(4)
      c53=cons*o35/g(5)
      c54=cons*o45/g(5)
      c13=c31*expfn(e31/tev)*g(3)/g(1)
      c14=c41*expfn(e41/tev)*g(4)/g(1)
      c15=c51*expfn(e51/tev)*g(5)/g(1)
      c23=c32*expfn(e32/tev)*g(3)/g(2)
      c24=c42*expfn(e42/tev)*g(4)/g(2)
      c25=c52*expfn(e52/tev)*g(5)/g(2)
      c34=c43*expfn(e43/tev)*g(4)/g(3)
      c35=c53*expfn(e53/tev)*g(5)/g(3)
      c45=c54*expfn(e54/tev)*g(5)/g(4)
       ENDIF
c
      amass=ama(i)
      m=0
 9988 continue
      m=m+1
      adamp=a21(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(2)*g(1))/(x(1)*g(2)))
      call escacf(e21,a21(i),adamp,g(1),g(2),te,amass,tcorr,
     &coltoti,coltot,p21)
       IF((E32/TEV).lt.2.d1) then
      adamp=a31(i)+a32(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(3)*g(1))/(x(1)*g(3)))
      call escacf(e31,a31(i),adamp,g(1),g(3),te,amass,tcorr,
     &coltoti,coltot,p31)
      adamp=a41(i)+a42(i)+a43(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(4)*g(1))/(x(1)*g(4)))
      call escacf(e41,a41(i),adamp,g(1),g(4),te,amass,tcorr,
     &coltoti,coltot,p41)
      adamp=a21(i)+a32(i)+a31(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(3)*g(2))/(x(2)*g(3)))
      call escacf(e32,a32(i),adamp,g(2),g(3),te,amass,tcorr,
     &coltoti,coltot,p32)
      adamp=a21(i)+a41(i)+a42(i)+a43(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(4)*g(2))/(x(2)*g(4)))
      call escacf(e42,a42(i),adamp,g(2),g(4),te,amass,tcorr,
     &coltoti,coltot,p42)
      adamp=a21(i)+a52(i)+a54(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(5)*g(2))/(x(2)*g(5)))
      call escacf(e52,a52(i),adamp,g(2),g(5),te,amass,tcorr,
     &coltoti,coltot,p52)
      adamp=a31(i)+a32(i)+a41(i)+a42(i)+a43(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(3)*(1.d0-(x(4)*g(3))/(x(3)*g(4)))
      call escacf(e43,a43(i),adamp,g(3),g(4),te,amass,tcorr,
     &coltoti,coltot,p43)
      adamp=a52(i)+a54(i)+a41(i)+a42(i)+a43(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(4)*(1.d0-(x(5)*g(4))/(x(4)*g(5)))
      call escacf(e54,a54(i),adamp,g(4),g(5),te,amass,tcorr,
     &coltoti,coltot,p54)
      do 100 j=1,5
      x(j)=0.d0
      a(1,j)=1.d0
  100 continue
      x(1)=1.d0
      a(2,1)=-c12
      a(2,2)=c21+c23+c24+c25+p21*a21(i)
      a(2,3)=-(c32+p32*a32(i))
      a(2,4)=-(c42+p42*a42(i))
      a(2,5)=-(c52+p52*a52(i))
      a(3,1)=-c13
      a(3,2)=-c23
      a(3,3)=c31+c32+c34+c35+p31*a31(i)+p32*a32(i)
      a(3,4)=-(c43+p43*a43(i))
      a(3,5)=-c53
      a(4,1)=-c14
      a(4,2)=-c24
      a(4,3)=-c34
      a(4,4)=c41+c42+c43+c45+p41*a41(i)+p42*a42(i)+p43*a43(i)
      a(4,5)=-(c54+p54*a54(i))
      a(5,1)=-c15
      a(5,2)=-c25
      a(5,3)=-c35
      a(5,4)=-c45
      a(5,5)=c51+c52+c53+c54+p52*a52(i)+p54*a54(i)
c
      call matrix(5)
       ELSE
      a(1,1)=1.d0
      a(1,2)=1.d0
      a(2,1)=-c12
      a(2,2)=c21+p21*a21(i)
      x(1)=1.d0
      x(2)=0.0d0
      call matrix(2)
       ENDIF
      if(m.lt.2) goto 9988
c
      elab=(1.6022d-12*abund)/(12.57d0*dene)
       IF((E32/TEV).lt.2.d1) then
      cool=elab*(e21*(x(1)*c12-x(2)*c21)+e31*
     &(x(1)*c13-x(3)*c31)+e41*(x(1)*c14-x(4)*c41)+e51*(x(1)*c15-
     &x(5)*c51)+e32*(x(2)*c23-x(3)*c32)+
     &e42*(x(2)*c24-x(4)*c42)+e52*(x(2)*c25-x(5)*c52)+
     &e43*(x(3)*c34-x(4)*c43)+e53*(x(3)*c35-x(5)*c53)+
     &e54*(x(4)*c45-x(5)*c54))
      f345t12=elab*(x(4)*(p41*e41*a41(i)+p42*e42*a42(i))+
     &x(5)*p52*e52*a52(i)+x(3)*(p31*e31*a31(i)+p32*e32*a32(i)))
      f43=elab*x(4)*a43(i)*e43*p43
      f54=elab*x(5)*a54(i)*e54*p54
       ELSE
      cool=elab*(e21*(x(1)*c12-x(2)*c21))
       ENDIF
      cool=12.57d0*cool
      f21=elab*x(2)*a21(i)*e21*p21
c
      endif
C
      RETURN
      END
C*********************************************************************
C*********************************************************************
      SUBROUTINE CFORB(TE,I,DENE,ABUND,COLTOTI,COLTOT,F21,F32,F4T123,
     &F5T23,F6T23,F54,COOL)
C*********************************************************************
C  Forbidden lines (C I-like). Enumeration of the ions are: (= I)
C
C  [last updated May 16, 1997]
C
C  1 = C I  2 = N II  3 = O III  4 = Ne V  5 = S III  6 = Ar V  
C  7 = Ca VII  8 = Mg VII  9 = Fe XIII
C
C  The atom is treated as a 6-level atom.
C  (Reduced to a 3-level atom for E43/Tev > 20)
C
C  WARNING !! Note that coll.strengths may be for multiplets. om23t5 and
C      om23t6 are for ALL transitions from levels 1,2 & 3 to 5 & 6.
C
C  Energy levels for O III from Bashkin & Stoner '75 (experimental)
C  Ne V from Lennon & Burke (MNRAS 251, 628; 1991 & AAS 103, 273; 1994).
C  (Ne V Omegas for fine structure lines may be wrong by up to a factor
C  of 3; Oliva et al. 1996, A&A 305, L21.)
C  S III 5S2 A-values from Hayes 1986, S III Coll.strengths from Tayal
C  (ApJ 1997), also from S III Galavis et al. (AAS, 111, 347, '95). Ar V and 
C  Ca VII from Galavis et al. Ca VII multiplet levels from Nahar & Pradhan
C  (J.Phys.B, 26, 1103, '93.). Fine structure levels from Bashkin & 
C  Stoner ('78). A-values for Ca VII put equal to O III A-values 
C  (which should give a lower limit to the cooling.)
C  A-values A62,A63 of N II from Brage et al. (ApJ, 478, 423, '97).
C  S, Ar, Ca & Mg  essentially treated as a 5-level atom. 
C  Mg Omegas from Lennon & Burke (1994). A-values and level energies 
C  from Nussbaumer & Rusca (A&A 72, 129, '79).
C  Fe XIII energy levels from Corliss & Sugar ('82) and Fawcett 
C  (J.Phys.B. '71). A-values from J.Phys.Chem.Ref. Data ('88), 
C  and Omegas from Tayal (ApJ 446, 895, '95). (Omegas are probably
C  valid only for T > 1E5 K.)
C  C I coll. strengths from Johnson, Burke & Kingston, J.Phys.B 20, 2553 
C  (1987). O III from Aggarwal, ApJ Suppl. 85, 197 (1993) and N II from 
C  Stafford et al., MNRAS 268, 816 (1994)
C  A-values of C, N, O, Ne and Mg updated with Galavis et al. 1997, A&AS,
C  123, 159
C*********************************************************************
      implicit real*8(a-h,o-z)
C
      parameter (mmatr=80)
C
      COMMON/MSTUFF/X(mmatr),A(mmatr,mmatr),EPS,NERR
c
      dimension ene2(9),ene3(9),ene4(9),ene5(9),ene6(9),a21(9),
     &a31(9),a32(9),a41(9),a42(9),a43(9),a52(9),a53(9),a54(9),
     &a62(9),a63(9),om12(9),om13(9),om23(9),om123t4(9),om23t5(9),
     &om45(9),om23t6(9),al12(9),al13(9),al23(9),al123t4(9),
     &al23t5(9),al45(9),al23t6(9),g(6),ama(9),a0(3,2),a1(3,2),
     &a2(3,2),a3(3,2),a4(3,2),yy(3)
c
      data ene2/.002033,.006038,.01403,.05113,.03685,.09471,
     &.0202,.1398,1.1531/
      data ene3/.005381,.01622,.03796,.1376,.1032,.2516,.0505,
     &.3645,2.3008/
      data ene4/1.2637,1.8989,2.4857,3.7555,1.4035,2.0208,2.394,
     &5.080,5.9559/
      data ene5/2.6839,4.0528,5.3250,7.9240,3.3677,4.3492,5.755,
     &10.562,11.338/
      data ene6/4.1825,5.8004,7.4490,10.9553,7.200,3*12.0,35.589/
C      data a21/7.93d-8,2.08d-6,2.66d-5,1.27d-3,4.72d-4,8.0d-3,
C     &2.66d-5,2.51d-2,14.d0/
C      data a31/2.05d-14,1.13d-12,3.09d-11,4.97d-9,4.61d-8,
C     &1.2d-6,3.09d-11,2.32d-7,6.3d-3/
C      data a32/2.65d-7,7.46d-6,9.70d-5,4.59d-3,2.07d-3,2.7d-2,
C     &9.70d-5,8.05d-2,9.86d0/
C      data a41/5.92d-8,3.55d-7,1.69d-6,1.94d-5,5.82d-6,3.5d-5,
C     &1.69d-6,1.16d-4,0.d0/
C      data a42/7.48d-5,1.02d-3,6.99d-3,1.25d-1,2.21d-2,2.0d-1,
C     &6.99d-3,1.19d0,63.04d0/
C      data a43/2.23d-4,3.00d-3,2.04d-2,3.50d-1,5.76d-2,4.8d-1,
C     &2.04d-2,3.11d0,75.07d0/
C      data a52/2.38d-3,3.30d-2,2.27d-1,3.99d0,7.96d-1,6.6d0,
C     &2.27d-1,35.8d0,1.0d3/
C      data a53/2.12d-5,1.32d-4,6.09d-4,6.29d-3,1.05d-2,5.7d-2,
C     &6.09d-4,3.61d-2,3.8d0/
C      data a54/6.42d-1,1.02d0,1.56d0,2.83d0,2.22d0,3.3d0,
C     &1.56d0,3.96d0,8.1d0/
      data a21/7.93d-8,2.1d-6,2.6d-5,1.27d-3,4.72d-4,8.0d-3,
     &2.66d-5,2.51d-2,14.d0/
      data a31/2.05d-14,1.2d-12,3.0d-11,4.97d-9,4.61d-8,
     &1.2d-6,3.09d-11,2.32d-7,6.3d-3/
      data a32/2.65d-7,7.5d-6,9.8d-5,4.59d-3,2.07d-3,2.7d-2,
     &9.70d-5,8.05d-2,9.86d0/
      data a41/5.92d-8,5.4d-7,2.7d-6,1.94d-5,5.82d-6,3.5d-5,
     &1.69d-6,1.16d-4,0.d0/
      data a42/7.48d-5,1.0d-3,6.7d-3,1.25d-1,2.21d-2,2.0d-1,
     &6.99d-3,1.19d0,63.04d0/
      data a43/2.23d-4,3.0d-3,2.0d-2,3.50d-1,5.76d-2,4.8d-1,
     &2.04d-2,3.11d0,75.07d0/
      data a52/2.38d-3,3.4d-2,2.2d-1,3.99d0,7.96d-1,6.6d0,
     &2.27d-1,35.8d0,1.0d3/
      data a53/2.12d-5,1.5d-4,7.8d-4,6.29d-3,1.05d-2,5.7d-2,
     &6.09d-4,3.61d-2,3.8d0/
      data a54/6.42d-1,1.1d0,1.8d0,2.83d0,2.22d0,3.3d0,
     &1.56d0,3.96d0,8.1d0/
C     data a62/6.94d0,4.8d1,2.12d2,2.37d3,5*0./
C     data a63/1.56d1,1.07d2,5.22d2,6.06d3,5*0./
      data a62/6.94d0,5.36d1,2.12d2,2.37d3,6.22d3,4*0./
      data a63/1.56d1,1.306d2,5.22d2,6.06d3,1.70d4,4*0./
      data om12/0.,0.401,0.545,1.408,2.33,2.94,1.15,0.34,10.3/
      data om13/0.,0.279,0.271,1.81,1.11,1.84,1.00,0.30,3.96/
      data om23/0.,1.13,1.29,5.832,5.41,7.81,3.70,1.08,20.3/
      data om123t4/1.14,2.68,2.29,2.1,7.95,3.21,4.25,0.86,
     &0.0/
      data om23t5/0.252,0.352,0.292,0.246,1.15,0.559,0.86,
     &0.18,0.0/
      data om45/0.277,0.411,0.582,0.65,1.30,1.70,0.63,0.4,
     &0.668/
      data om23t6/0.671,1.28,1.18,1.43,5*1.d-10/
      data al12/3*0.,-.328,-.132,-.24,.142,0.11,-.684/
      data al13/2*0.,.044,-.494,0.07,-.155,.161,0.11,-.544/
      data al23/2*0.,.068,-.443,0.,-.194,.153,0.10,-.598/
      data al123t4/.969,0.,.129,0.,0.,.142,-.082,0.095,0./
      data al23t5/.896,0.,.287,2*0.,.103,-.233,-.074,0./
      data al45/.500,.035,.085,0.,.275,0.,.279,-.13,-.101/
      data al23t6/.616,0.,.120,6*0./
      data a0/-9.25141,-7.69735,-7.4387,444.6,350.609,386.186/
      data a1/-.773782,-1.30743,-.57443,-227.913,-187.474,
     &-202.192/
      data a2/.361184,.697638,.358264,42.5952,36.1803,38.5049/
      data a3/-.0150892,-.111338,-.0418166,-3.4762,-3.03283,
     &-3.19268/
      data a4/-.000656325,-.00705277,-.00235272,.105085,
     &.0938138,.0978573/
      data g/1.,3.,5.,5.,1.,5./
      data ama/12.,14.,16.,20.,32.1,39.948,40.08,24.305,55.847/
C
      cool=0.0d0
      f21=0.0d0
      f32=0.0d0
      f4t123=0.0d0
      f5t23=0.0d0
      f6t23=0.0d0
      f54=0.0d0
      if(abund.ge.1.d-10) then
c
      tev=te/1.1609d4
      t4=te/1.d4
      telog=dlog10(te)
      th=dmax1(telog,3.4d0)
      th=dmin1(5.3d0,th)
      tmx=1.25d5
      cons=8.63d-6*dene/dsqrt(te)
c
      e21=ene2(i)
      e31=ene3(i)
      e41=ene4(i)
      e51=ene5(i)
      e61=ene6(i)
      e32=ene3(i)-ene2(i)
      e42=ene4(i)-ene2(i)
      e52=ene5(i)-ene2(i)
      e62=ene6(i)-ene2(i)
      e43=ene4(i)-ene3(i)
      e53=ene5(i)-ene3(i)
      e63=ene6(i)-ene3(i)
      e54=ene5(i)-ene4(i)
      if(i.eq.1) then
       y=dmin1(dlog(te),9.2103d0)
       ipc=1
       if(te.gt.1.d3) ipc=2
        do 5001 ip=1,3
        yy(ip)=a0(ip,ipc)+a1(ip,ipc)*y+a2(ip,ipc)*y**2+
     &  a3(ip,ipc)*y**3+a4(ip,ipc)*y**4
 5001   continue
       o12=dexp(yy(1))
       o13=dexp(yy(2))
       o23=dexp(yy(3))
      elseif(i.eq.2) then
       o12=cheb1(tmx,te,-1.534d0,0.1436d0,0.0107d0,-0.0076d0)
       o13=cheb1(tmx,te,-2.4004d0,0.1962d0,-0.044d0,-0.0296d0)
       o23=cheb1(tmx,te,0.4597d0,0.1717d0,-0.0189d0,-0.0198d0)
      else
       o12=om12(i)*t4**al12(i)
       o13=om13(i)*t4**al13(i)
       o23=om23(i)*t4**al23(i)
      endif
      c21=cons*o12/g(2)
      c31=cons*o13/g(3)
      c32=cons*o23/g(3)
      c12=c21*expfn(e21/tev)*g(2)/g(1)
      c13=c31*expfn(e31/tev)*g(3)/g(1)
      c23=c32*expfn(e32/tev)*g(3)/g(2)
C
       IF((E43/TEV).lt.2.d1) then
      gsum=g(1)+g(2)+g(3)
      if(i.eq.3) then
       fc=(g(2)+g(3))/gsum
       o123t4=35.674-25.2742*th+6.2041*th**2-0.4944*th**3
       o23t5=fc*(5.9018-4.1547*th+1.0034*th**2-0.0788*th**3)
       o23t6=fc*(14.1717-10.5248*th+2.7476*th**2-0.233*th**3)
       o45=42.233-40.2293*th+14.2334*th**2-2.196*th**3+0.125*
     & th**4
      elseif(i.eq.5) then
       fc=(g(2)+g(3))/gsum
       if(te.le.3.d4) then
         o123t4=7.37d0*(te/3.d4)**0.0146
         o23t5=fc*1.44d0*(te/3.d4)**0.0614
         o23t6=fc*3.35d0*(te/3.d4)**0.112
       else
	 o123t4=7.37d0*(te/3.d4)**(-0.176)
         o23t5=fc*1.44d0*(te/3.d4)**(-0.124)
         o23t6=fc*3.35d0*(te/3.d4)**(-0.158)
       endif
       if(te.le.6.d4) then
         o45=1.88d0*(te/6.d4)**0.116
       else
         o45=1.88d0*(te/6.d4)**(-0.085)
       endif
      elseif(i.eq.2) then
       o14=cheb1(tmx,te,-2.1882d0,-0.0023d0,0.0082d0,-0.0005d0)
       o15=cheb1(tmx,te,-6.3792d0,-0.0076d0,0.019d0,-0.0013d0)
       o16=cheb1(tmx,te,-4.1428d0,-0.0268d0,-0.0218d0,-0.0125d0)
       o24=cheb1(tmx,te,0.009d0,-0.0023d0,0.0082d0,-0.0006d0)
       o25=cheb1(tmx,te,-4.182d0,-0.0076d0,0.019d0,-0.0013d0)
       o26=cheb1(tmx,te,-1.9489d0,-0.0296d0,-0.0239d0,-0.0139d0)
       o34=cheb1(tmx,te,1.0307d0,-0.0023d0,0.0082d0,-0.0006d0)
       o35=cheb1(tmx,te,-3.1604d0,-0.0076d0,0.019d0,-0.0013d0)
       o36=cheb1(tmx,te,-0.9239d0,-0.0268d0,-0.0217d0,-0.0125d0)
       o45=cheb1(tmx,te,-1.2463d0,0.0961d0,0.0291d0,-0.0024d0)
      elseif(i.eq.9) then
       o14=4.11d0*t4**(-.367)
       o15=0.32d0*t4**(-.705)
       o16=1.d-10
       o24=2.84d0*t4**(-.594)
       o25=0.36d0*t4**(-.655)
       o26=1.d-10
       o34=10.5d0*t4**(-.589)
       o35=0.20d0*t4**(-.266)
       o36=1.d-10
       o45=0.668d0*t4**(-.101)
      else
       fc=(g(2)+g(3))/gsum
       o123t4=om123t4(i)*t4**al123t4(i)
       o23t5=fc*om23t5(i)*t4**al23t5(i)
       o23t6=fc*om23t6(i)*t4**al23t6(i)
       o45=om45(i)*t4**al45(i)
      endif
c
      if(i.eq.1.or.i.ge.3) then
       if(i.le.8) then
       o14=g(1)*o123t4/gsum
       o24=g(2)*o14/g(1)
       o34=g(3)*o14/g(1)
       o25=g(2)*o23t5/(g(2)+g(3))
       o15=g(1)*o25/g(2)
       o35=g(3)*o25/g(2)
       o26=g(2)*o23t6/(g(2)+g(3))
       o16=g(1)*o26/g(2)
       o36=g(3)*o26/g(2)
       endif
      endif
      c41=cons*o14/g(4)
      c51=cons*o15/g(5)
      c61=cons*o16/g(6)
      c42=cons*o24/g(4)
      c52=cons*o25/g(5)
      c62=cons*o26/g(6)
      c43=cons*o34/g(4)
      c53=cons*o35/g(5)
      c63=cons*o36/g(6)
      c54=cons*o45/g(5)
      c14=c41*expfn(e41/tev)*g(4)/g(1)
      c15=c51*expfn(e51/tev)*g(5)/g(1)
      c16=c61*expfn(e61/tev)*g(6)/g(1)
      c24=c42*expfn(e42/tev)*g(4)/g(2)
      c25=c52*expfn(e52/tev)*g(5)/g(2)
      c26=c62*expfn(e62/tev)*g(6)/g(2)
      c34=c43*expfn(e43/tev)*g(4)/g(3)
      c35=c53*expfn(e53/tev)*g(5)/g(3)
      c36=c63*expfn(e63/tev)*g(6)/g(3)
      c45=c54*expfn(e54/tev)*g(5)/g(4)
       ENDIF
c
C      write(6,6485)te,(o14+o24+o34),(o15+o25+o35),o45
 6485 format(1x,6e12.5)
      amass=ama(i)
      m=0
 9988 continue
      m=m+1
      adamp=a21(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(2)*g(1))/(x(1)*g(2)))
      call escacf(e21,a21(i),adamp,g(1),g(2),te,amass,tcorr,
     &coltoti,coltot,p21)
      adamp=a31(i)+a32(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(3)*g(1))/(x(1)*g(3)))
      call escacf(e31,a31(i),adamp,g(1),g(3),te,amass,tcorr,
     &coltoti,coltot,p31)
      adamp=a21(i)+a32(i)+a31(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(3)*g(2))/(x(2)*g(3)))
      call escacf(e32,a32(i),adamp,g(2),g(3),te,amass,tcorr,
     &coltoti,coltot,p32)
       IF((E43/TEV).lt.2.d1) then
      adamp=a41(i)+a42(i)+a43(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(4)*g(1))/(x(1)*g(4)))
      call escacf(e41,a41(i),adamp,g(1),g(4),te,amass,tcorr,
     &coltoti,coltot,p41)
      adamp=a21(i)+a41(i)+a42(i)+a43(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(4)*g(2))/(x(2)*g(4)))
      call escacf(e42,a42(i),adamp,g(2),g(4),te,amass,tcorr,
     &coltoti,coltot,p42)
      adamp=a21(i)+a52(i)+a54(i)+a53(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(5)*g(2))/(x(2)*g(5)))
      call escacf(e52,a52(i),adamp,g(2),g(5),te,amass,tcorr,
     &coltoti,coltot,p52)
      adamp=a31(i)+a32(i)+a41(i)+a42(i)+a43(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(3)*(1.d0-(x(4)*g(3))/(x(3)*g(4)))
      call escacf(e43,a43(i),adamp,g(3),g(4),te,amass,tcorr,
     &coltoti,coltot,p43)
      adamp=a52(i)+a54(i)+a53(i)+a31(i)+a32(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(3)*(1.d0-(x(5)*g(3))/(x(3)*g(5)))
      call escacf(e53,a53(i),adamp,g(3),g(5),te,amass,tcorr,
     &coltoti,coltot,p53)
      adamp=a52(i)+a54(i)+a53(i)+a41(i)+a42(i)+a43(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(4)*(1.d0-(x(5)*g(4))/(x(4)*g(5)))
      call escacf(e54,a54(i),adamp,g(4),g(5),te,amass,tcorr,
     &coltoti,coltot,p54)
      adamp=a62(i)+a63(i)+a21(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(6)*g(2))/(x(2)*g(6)))
      call escacf(e62,a62(i),adamp,g(2),g(6),te,amass,tcorr,
     &coltoti,coltot,p62)
      adamp=a62(i)+a63(i)+a31(i)+a32(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(3)*(1.d0-(x(6)*g(3))/(x(3)*g(6)))
      call escacf(e63,a63(i),adamp,g(3),g(6),te,amass,tcorr,
     &coltoti,coltot,p63)
      a(1,1)=c12+c13+c14+c15+c16
      a(1,2)=-(c21+p21*a21(i))
      a(1,3)=-(c31+p31*a31(i))
      a(1,4)=-(c41+p41*a41(i))
      a(1,5)=-c51
      a(1,6)=-c61
      a(2,1)=-c12
      a(2,2)=c21+c23+c24+c25+c26+p21*a21(i)
      a(2,3)=-(c32+p32*a32(i))
      a(2,4)=-(c42+p42*a42(i))
      a(2,5)=-(c52+p52*a52(i))
      a(2,6)=-(c62+p62*a62(i))
      do 100 j=1,6
      x(j)=0.d0
      a(3,j)=1.d0
  100 continue
      x(3)=1.d0
      a(4,1)=-c14
      a(4,2)=-c24
      a(4,3)=-c34
      a(4,4)=c41+c42+c43+c45+p41*a41(i)+p42*a42(i)+p43*a43(i)
      a(4,5)=-(c54+p54*a54(i))
      a(4,6)=0.0d0
      a(5,1)=-c15
      a(5,2)=-c25
      a(5,3)=-c35
      a(5,4)=-c45
      a(5,5)=c51+c52+c53+c54+p52*a52(i)+p53*a53(i)+p54*a54(i)
      a(5,6)=0.0d0
      a(6,1)=-c16
      a(6,2)=-c26
      a(6,3)=-c36
      a(6,4)=0.0d0
      a(6,5)=0.0d0
      a(6,6)=c61+c62+c63+p62*a62(i)+p63*a63(i)
c
      call matrix(6)
       ELSE
      a(1,1)=c12+c13
      a(1,2)=-(c21+p21*a21(i))
      a(1,3)=-(c31+p31*a31(i))
      a(2,1)=-c12
      a(2,2)=c21+c23+p21*a21(i)
      a(2,3)=-(c32+p32*a32(i))
      a(3,1)=1.d0
      a(3,2)=1.d0
      a(3,3)=1.d0
      x(1)=0.0d0
      x(2)=0.0d0
      x(3)=1.d0
c
      call matrix(3)
       ENDIF
      if(m.lt.2) goto 9988
c
      elab=(1.6022d-12*abund)/(12.57d0*dene)
       IF((E43/TEV).lt.2.d1) then
      cool=elab*(e21*(x(1)*c12-x(2)*c21)+e31*
     &(x(1)*c13-x(3)*c31)+e41*(x(1)*c14-x(4)*c41)+e51*(x(1)*c15-
     &x(5)*c51)+e61*(x(1)*c16-x(6)*c61)+e32*(x(2)*c23-x(3)*c32)+
     &e42*(x(2)*c24-x(4)*c42)+e52*(x(2)*c25-x(5)*c52)+e62*(x(2)*
     &c26-x(6)*c62)+e43*(x(3)*c34-x(4)*c43)+e53*(x(3)*c35-x(5)*
     &c53)+e63*(x(3)*c36-x(6)*c63)+e54*(x(4)*c45-x(5)*c54))
      f4t123=elab*x(4)*(p41*e41*a41(i)+p42*e42*a42(i)+
     &p43*e43*a43(i))
      f5t23=elab*x(5)*(p52*e52*a52(i)+p53*e53*a53(i))
      f6t23=elab*x(6)*(p62*e62*a62(i)+p63*e63*a63(i))
      f54=elab*x(5)*e54*a54(i)*p54
       ELSE
      cool=elab*(e21*(x(1)*c12-x(2)*c21)+e31*(x(1)*c13-x(3)*c31)+
     &e32*(x(2)*c23-x(3)*c32))
       ENDIF
      cool=12.57d0*cool
      f21=elab*x(2)*a21(i)*e21*p21
      f32=elab*x(3)*a32(i)*e32*p32
C      if(i.eq.4.or.i.eq.6) then
         decrit=(a43(i)*dene)/(c41+c42+c43+c45+c46)
C         write(6,1818)i,decrit
C      endif
 1818  format(1x,i4,3e14.5)
c
      endif
C
      RETURN
      END
C*********************************************************************
C*********************************************************************
      SUBROUTINE NFORB(TE,I,DENE,ABUND,COLTOTI,COLTOT,F23T1,F45T1,
     &F45T23,COOL)
C*********************************************************************
C  Forbidden lines (N I-like). Enumeration of the ions are: (= I)
C
C  1 = N I  2 = O II  3 = Ne IV  4 = S II  5 = Na V  6 = Ar IV
C  7 = Mg VI  8 = Si VIII
C
C  The atom is treated as a 5-level atom.
C  No levels are calculated if E21/TEV > 20.
C
C  Note that some of the levels have inverted energy levels !
C  N I Omegas Mendoza ('83) and A-values from Kaufman & Sugar (J.Phys.
C  Chem. '86). O II, Ne IV & Ar IV from Osterbrock ('89). 
C  O II & Ne IV checked against the original papers, and small changes
C  may occur. O II Omegas updated with McLaughlin & Bell (J.Phys.B '93).
C  New A-values for S II (Keenan et al., Phys. Scripta 48, 129 (1993). 
C  S II (Omegas) updated with Ramsbottom, Bell & Stafford (1996) [Different
C  from Cai & Pradhan, ApjS 88, 329 (1993)]. 
C  Na V from Mendoza ('83). A-values, energies and 10Ryd-Omegas for 
C  Mg VI and Si VIII from Bhatia & Manson (MNRAS 190, 925, '80).
C  (From isoel.scaling this actually gives rather good Omegas.)
C*********************************************************************
      implicit real*8(a-h,o-z)
C
      parameter (mmatr=80)
C
      COMMON/MSTUFF/X(mmatr),A(mmatr,mmatr),EPS,NERR
      common/robcum/fro21,fro41,frs21,frs41,frr21,frr41
c
      dimension ene2(8),ene3(8),ene4(8),ene5(8),a21(8),a31(8),
     &a32(8),a41(8),a42(8),a43(8),a51(8),a52(8),a53(8),a54(8),
     &om12(8),om13(8),om14(8),om15(8),om23(8),om24(8),om25(8),
     &om34(8),om35(8),om45(8),al12(8),al13(8),al14(8),al15(8),
     &al23(8),al24(8),al25(8),al34(8),al35(8),al45(8),g(5),ama(8)
c
      data ene2/2.381916,3.321846,5.108967,1.840287,5.9860,
     &2.6131,7.2166,8.9326/
      data ene3/2.382994,3.324262,5.114530,1.844227,5.9917,
     &2.6291,7.2171,8.9725/
      data ene4/3.573140,5.013923,7.735647,3.038635,9.0697,
     &4.3186,10.233,12.878/
      data ene5/3.573189,5.014060,7.736477,3.044446,9.0740,
     &4.3405,10.243,12.943/
      data a21/6.92d-6,3.54d-5,4.84d-4,8.9d-4,1.46d-3,2.23d-2,
     &0.145d0,1.96d0/
      data a31/1.62d-5,1.80d-4,5.54d-3,2.73d-4,2.69d-2,1.77d-3,
     &3.11d-3,3.75d-2/
      data a41/2.46d-3,5.75d-2,5.21d-1,9.06d-2,1.76d0,8.62d-1,
     &4.92d0,29.9d0/
      data a51/6.18d-3,2.36d-2,1.27d0,2.25d-1,4.27d0,2.6d0,
     &12.d0,71.2d0/
      data a32/1.07d-8,1.31d-7,1.48d-6,3.35d-7,1.55d-6,2.3d-5,
     &3.5d-5,3.54d-4/
      data a42/3.45d-2,1.07d-1,1.15d-1,1.41d-1,1.41d-1,6.03d-1,
     &1.85d0,11.4d0/
      data a52/5.48d-2,5.63d-2,4.00d-1,1.21d-1,9.18d-1,7.89d-1,
     &2.99d0,19.d0/
      data a43/4.71d-2,5.79d-2,3.93d-1,6.27d-2,9.55d-1,1.19d-1,
     &7.67d-2,0.119d0/
      data a53/2.47d-2,9.40d-2,4.37d-1,1.52d-1,1.29d0,5.98d1,
     &1.69d0,10.5d0/
      data a54/5.17d-13,1.35d-10,2.68d-9,1.03d-6,3.64d-7,
     &4.94d-5,4.36d-6,1.29d-3/
      data om12/.290,.825,.838,2.76,.551,1.30,.19,.15/
      data om13/.194,.55,.559,4.14,.368,1.94,.22,.18/
      data om14/.0567,.276,.156,1.17,.12,.147,.058,.048/
      data om15/.113,.138,.313,2.35,.239,.293,.10,.084/
      data om23/.269,1.17,1.36,7.47,.696,6.13,.28,.22/
      data om24/.109,.882,.368,1.79,.201,1.67,.15,.10/
      data om25/.266,.356,.900,3.00,.502,2.47,.19,.13/
      data om34/.097,.493,.336,2.20,.19,1.79,.14,.10/
      data om35/.147,.332,.509,4.99,.279,4.44,.37,.26/
      data om45/.071,.287,.343,2.71,.438,2.33,.15,.12/
      data al12/.904,.038,0.,-.161,4*0./
      data al13/.913,.038,0.,-.160,4*0./
      data al14/.928,.0738,0.,-.104,4*0./
      data al15/.921,.0738,0.,-.101,4*0./
      data al23/1.07,2*0.,-.135,4*0./
      data al24/.800,.1357,.054,-.153,4*0./
      data al25/.715,.1357,0.,-.120,4*0./
      data al34/.691,.1357,0.,-.115,4*0./
      data al35/.780,.1357,0.,-.141,4*0./
      data al45/1.110,.049,.109,-.158,4*0./
      data ama/14.,16.,20.,32.1,22.99,39.948,24.305,28.086/
C
      cool=0.0d0
      f23t1=0.0d0
      f45t1=0.0d0
      f45t23=0.0d0
      frr21=0.0d0
      frr41=0.0d0
      if(abund.ge.1.d-10) then
c
      tev=te/1.1609d4
      t4=te/1.d4
      telog=dlog10(te)
      th=telog
      if(th.le.3.d0) th=3.d0
      if(th.ge.4.6d0) th=4.6d0
      cons=8.63d-6*dene/dsqrt(te)
c
      g(1)=4.d0
      g(2)=4.d0
      g(3)=6.d0
      if(i.eq.1.or.i.eq.2) then
       g(2)=6.d0
       g(3)=4.d0
      endif
      if(i.eq.3.or.i.eq.5) then
       g(2)=6.d0
       g(3)=4.d0
      endif
      if(i.ne.2) then
       g(4)=2.d0
       g(5)=4.d0
      else
       g(4)=4.d0
       g(5)=2.d0
      endif
      e21=ene2(i)
       IF((E21/TEV).lt.2.d1) THEN
      e31=ene3(i)
      e41=ene4(i)
      e51=ene5(i)
      e32=ene3(i)-ene2(i)
      e42=ene4(i)-ene2(i)
      e52=ene5(i)-ene2(i)
      e43=ene4(i)-ene3(i)
      e53=ene5(i)-ene3(i)
      e54=ene5(i)-ene4(i)
      if(i.ne.4) then
       o12=om12(i)*t4**al12(i)
       o13=om13(i)*t4**al13(i)
       o14=om14(i)*t4**al14(i)
       o15=om15(i)*t4**al15(i)
       o23=om23(i)*t4**al23(i)
       o24=om24(i)*t4**al24(i)
       o25=om25(i)*t4**al25(i)
       o34=om34(i)*t4**al34(i)
       o35=om35(i)*t4**al35(i)
       o45=om45(i)*t4**al45(i)
      else
       o12=3.9949-1.1021*th+0.5487*th**2-0.0823*th**3
       o13=6.1758-1.836*th+0.8798*th**2-0.1291*th**3
       o14=13.1138-10.7869*th+2.9634*th**2-0.2608*th**3
       o15=26.7492-22.005*th+6.0445*th**2-0.5322*th**3
       o23=3.2171+3.0844*th-0.3551*th**2-0.038*th**3
       o24=7.0174-5.7352*th+1.8524*th**2-0.187*th**3
       o25=13.3032-10.9238*th+3.4266*th**2-0.3392*th**3
       o34=10.6418-8.6917*th+2.6787*th**2-0.2618*th**3
       o35=23.5642-19.3646*th+6.0722*th**2-0.602*th**3
       o45=37.3916-31.082*th+8.6052*th**2-0.7627*th**3
      endif
      c21=cons*o12/g(2)
      c31=cons*o13/g(3)
      c41=cons*o14/g(4)
      c51=cons*o15/g(5)
      c32=cons*o23/g(3)
      c42=cons*o24/g(4)
      c52=cons*o25/g(5)
      c43=cons*o34/g(4)
      c53=cons*o35/g(5)
      c54=cons*o45/g(5)
      c12=c21*expfn(e21/tev)*g(2)/g(1)
      c13=c31*expfn(e31/tev)*g(3)/g(1)
      c14=c41*expfn(e41/tev)*g(4)/g(1)
      c15=c51*expfn(e51/tev)*g(5)/g(1)
      c23=c32*expfn(e32/tev)*g(3)/g(2)
      c24=c42*expfn(e42/tev)*g(4)/g(2)
      c25=c52*expfn(e52/tev)*g(5)/g(2)
      c34=c43*expfn(e43/tev)*g(4)/g(3)
      c35=c53*expfn(e53/tev)*g(5)/g(3)
      c45=c54*expfn(e54/tev)*g(5)/g(4)
c
      amass=ama(i)
      m=0
 9988 continue
      m=m+1
      adamp=a21(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(2)*g(1))/(x(1)*g(2)))
      call escacf(e21,a21(i),adamp,g(1),g(2),te,amass,tcorr,
     &coltoti,coltot,p21)
      adamp=a31(i)+a32(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(3)*g(1))/(x(1)*g(3)))
      call escacf(e31,a31(i),adamp,g(1),g(3),te,amass,tcorr,
     &coltoti,coltot,p31)
      adamp=a41(i)+a42(i)+a43(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(4)*g(1))/(x(1)*g(4)))
      call escacf(e41,a41(i),adamp,g(1),g(4),te,amass,tcorr,
     &coltoti,coltot,p41)
      adamp=a41(i)+a42(i)+a43(i)+a51(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(5)*g(1))/(x(1)*g(5)))
      call escacf(e51,a51(i),adamp,g(1),g(5),te,amass,tcorr,
     &coltoti,coltot,p51)
      adamp=a21(i)+a32(i)+a31(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(3)*g(2))/(x(2)*g(3)))
      call escacf(e32,a32(i),adamp,g(2),g(3),te,amass,tcorr,
     &coltoti,coltot,p32)
      adamp=a21(i)+a41(i)+a42(i)+a43(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(4)*g(2))/(x(2)*g(4)))
      call escacf(e42,a42(i),adamp,g(2),g(4),te,amass,tcorr,
     &coltoti,coltot,p42)
      adamp=a21(i)+a51(i)+a52(i)+a54(i)+a53(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(5)*g(2))/(x(2)*g(5)))
      call escacf(e52,a52(i),adamp,g(2),g(5),te,amass,tcorr,
     &coltoti,coltot,p52)
      adamp=a31(i)+a32(i)+a41(i)+a42(i)+a43(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(3)*(1.d0-(x(4)*g(3))/(x(3)*g(4)))
      call escacf(e43,a43(i),adamp,g(3),g(4),te,amass,tcorr,
     &coltoti,coltot,p43)
      adamp=a51(i)+a52(i)+a54(i)+a53(i)+a31(i)+a32(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(3)*(1.d0-(x(5)*g(3))/(x(3)*g(5)))
      call escacf(e53,a53(i),adamp,g(3),g(5),te,amass,tcorr,
     &coltoti,coltot,p53)
      adamp=a51(i)+a52(i)+a54(i)+a53(i)+a41(i)+a42(i)+a43(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(4)*(1.d0-(x(5)*g(4))/(x(4)*g(5)))
      call escacf(e54,a54(i),adamp,g(4),g(5),te,amass,tcorr,
     &coltoti,coltot,p54)
c
      a(1,1)=c12+c13+c14+c15
      a(1,2)=-(c21+p21*a21(i))
      a(1,3)=-(c31+p31*a31(i))
      a(1,4)=-(c41+p41*a41(i))
      a(1,5)=-(c51+p51*a51(i))
      a(2,1)=-c12
      a(2,2)=c21+c23+c24+c25+p21*a21(i)
      a(2,3)=-(c32+p32*a32(i))
      a(2,4)=-(c42+p42*a42(i))
      a(2,5)=-(c52+p52*a52(i))
      do 100 j=1,5
      x(j)=0.d0
      a(3,j)=1.d0
  100 continue
      x(3)=1.d0
      a(4,1)=-c14
      a(4,2)=-c24
      a(4,3)=-c34
      a(4,4)=c41+c42+c43+c45+p41*a41(i)+p42*a42(i)+p43*a43(i)
      a(4,5)=-(c54+p54*a54(i))
      a(5,1)=-c15
      a(5,2)=-c25
      a(5,3)=-c35
      a(5,4)=-c45
      a(5,5)=c51+c52+c53+c54+p51*a51(i)+p52*a52(i)+p53*a53(i)+
     &p54*a54(i)
c
      call matrix(5)
      if(m.lt.2) goto 9988
c
      elab=(1.6022d-12*abund)/(12.57d0*dene)
      cool=elab*(e21*(x(1)*c12-x(2)*c21)+e31*
     &(x(1)*c13-x(3)*c31)+e41*(x(1)*c14-x(4)*c41)+e51*(x(1)*c15-
     &x(5)*c51)+e32*(x(2)*c23-x(3)*c32)+
     &e42*(x(2)*c24-x(4)*c42)+e52*(x(2)*c25-x(5)*c52)+
     &e43*(x(3)*c34-x(4)*c43)+e53*(x(3)*c35-x(5)*c53)+
     &e54*(x(4)*c45-x(5)*c54))
      cool=12.57d0*cool
      f23t1=elab*(x(2)*p21*a21(i)*e21+x(3)*p31*a31(i)*e31)
      f45t1=elab*(x(4)*p41*a41(i)*e41+x(5)*p51*a51(i)*e51)
      f45t23=elab*(x(4)*(e42*a42(i)*p42+e43*a43(i)*p43)+x(5)*
     &(e52*a52(i)*p52+e53*a53(i)*p53))
C - line ratios -
      frr21=(elab*x(2)*p21*a21(i)*e21)/f23t1
      if(i.ne.2) then
       frr41=(elab*x(4)*p41*a41(i)*e41)/f45t1
      else
       frr41=(elab/f45t23)*(x(4)*e43*a43(i)*p43+x(5)*e53*a53(i)*p53)
      endif
c
       ENDIF
      endif
C
      RETURN
      END
C*********************************************************************
C*********************************************************************
      SUBROUTINE OFORB(TE,I,DENE,ABUND,COLTOTI,COLTOT,F21,F32,F4T123,
     &F5T12,F54,COOL)
C*********************************************************************
C  Forbidden lines (O I-like). Enumeration of the ions are: (= I)
C
C  1 = O I  2 = Ne III  3 = S I  4 = Na IV  5 = Mg V  6 = Ar III  
C  7 = Ca V
C
C  The atom is treated as a 5-level atom. (The energy levels of lower
C  3P-state are inverted compared to C I-like states.)
C  (The atom is reduced to a 3-level atom for E43/Tev > 20).
C
C  ene4 & ene5 of Ar III from Hansen & Persson ('87). Other energies 
C  from Mendoza ('83). A-values and Omegas from Mendoza ('83) and
C  Osterbrock ('89), except Omegas of Ar III and Ca V which are from
C  Galavis et al. (AAS, 111, 347, '95), and of O I (Berrington, 
C  J.Phys.B 21, 1083, '87). S I Omegas are assumed to follow those of 
C  O I, which is probably an underestimate.
C
C  A-values of O, Ne, Na and Mg updated with Galavis et al. 1997, A&AS,
C  123, 159.
C*********************************************************************
      implicit real*8(a-h,o-z)
C
      parameter (mmatr=80)
C
      COMMON/MSTUFF/X(mmatr),A(mmatr,mmatr),EPS,NERR
c
      dimension ene2(7),ene3(7),ene4(7),ene5(7),a21(7),a31(7),
     &a41(7),a51(7),a32(7),a42(7),a52(7),a43(7),a54(7),
     &om12(7),om13(7),om23(7),om123t4(7),om123t5(7),om45(7),
     &al12(7),al13(7),al23(7),al123t4(7),al123t5(7),al45(7),g(5),
     &ama(7)
c
      data ene2/.01965,.07966,.04908,.13707,.22093,.13779,.29794/
      data ene3/.02809,.1141,.07108,.19522,.31245,.19455,.40585/
      data ene4/1.9660,3.2017,1.1447,3.8211,4.4520,1.6685,2.3331/
      data ene5/4.1869,6.9075,2.7481,8.2389,9.5757,4.0455,5.4313/
      data a21/8.86d-5,5.97d-3,1.39d-3,3.05d-2,.127d0,3.08d-2,
     &.31d0/
      data a31/1.28d-10,2.08d-8,6.71d-8,1.60d-7,9.68d-7,2.37d-6,
     &3.67d-5/
      data a32/1.77d-5,1.16d-3,3.02d-4,5.56d-3,2.17d-2,5.17d-3,
     &3.54d-2/
      data a41/6.51d-3,.173d0,2.78d-2,.61d0,1.87d0,0.314d0,1.9d0/
      data a42/2.11d-3,5.34d-2,8.16d-3,.184d0,.535d0,8.23d-2,
     &.426d0/
      data a43/6.39d-7,8.27d-6,3.84d-6,2.20d-5,5.15d-5,2.21d-5,
     &8.42d-5/
      data a51/2.94d-4,3.98d-3,8.23d-3,1.07d-2,2.50d-2,4.17d-2,
     &.145d0/
      data a52/7.91d-2,2.03,.35d0,7.15d0,21.6d0,3.91d0,23.1d0/
      data a54/1.12d0,2.56d0,1.53d0,3.32d0,4.09d0,2.59d0,3.73d0/
      data om12/0.,.527,0.,.471,.4,3.1,2.3/
      data om13/0.,.131,0.,.111,.091,.671,.648/
      data om23/0.,.185,0.,.177,.156,1.261,.671/
      data om123t4/0.,1.34,0.,1.17,1.02,4.825,3.067/
      data om123t5/0.,.151,0.,.163,.146,.841,.522/
      data om45/0.,.236,0.,.157,.129,1.219,1.348/
      data al12/5*0.,.034,.163/
      data al13/5*0.,.16,.176/
      data al23/5*0.,-.098,.141/
      data al123t4/5*0.,-.065,.104/
      data al123t5/5*0.,-.099,.185/
      data al45/0.,.184,3*0.,-.152,0./
      data g/5.,3.,1.,5.,1./
      data ama/16.,20.,32.1,22.99,24.305,39.948,40.08/
C
      cool=0.0d0
      f21=0.0d0
      f32=0.0d0
      f4t123=0.0d0
      f5t12=0.0d0
      f54=0.0d0
      if(abund.ge.1.d-10) then
c
      tev=te/1.1609d4
      t4=te/1.d4
      telog=dlog10(te)
      th=dmax1(telog,3.4d0)
      th=dmin1(5.d0,th)
      cons=8.63d-6*dene/dsqrt(te)
c
      e21=ene2(i)
      e31=ene3(i)
      e41=ene4(i)
      e51=ene5(i)
      e32=ene3(i)-ene2(i)
      e42=ene4(i)-ene2(i)
      e52=ene5(i)-ene2(i)
      e43=ene4(i)-ene3(i)
      e53=ene5(i)-ene3(i)
      e54=ene5(i)-ene4(i)
c
      gsum=g(1)+g(2)+g(3)
C -- S I omegas assumed to follow those of O I (probably an
C    underestimate !) --
c -- O I from Berrington, J.Phys.B 21, 1083 (1987) --
      if(i.eq.1.or.i.eq.3) then
       o45=1.05d-1*t4**0.520d0
       if(te.le.1.d3) then
        o12=9.72d-3*(te/1.d3)**.773d0
        o13=4.11d-3*(te/1.d3)**.676d0
        o23=6.36d-4*(te/1.d3)**1.301d0
        o123t4=1.51d-2*(te/1.d3)**1.380d0
        o123t5=1.84d-3*(te/1.d3)**1.501d0
       else 
        o12=9.72d-3*(te/1.d3)**1.038d0
        o13=4.11d-3*(te/1.d3)**.893d0
        o23=6.36d-4*(te/1.d3)**1.648d0
        o123t4=2.66d-1*t4**1.246d0
        o123t5=3.24d-2*t4**1.246d0
       endif
      else
       o12=om12(i)*t4**al12(i)
       o13=om13(i)*t4**al13(i)
       o23=om23(i)*t4**al23(i)
       o123t4=om123t4(i)*t4**al123t4(i)
       o123t5=om123t5(i)*t4**al123t5(i)
       o45=om45(i)*t4**al45(i)
      endif
       if(i.eq.7.and.o123t5.gt.0.8d0) o123t5=0.8d0
      c21=cons*o12/g(2)
      c31=cons*o13/g(3)
      c32=cons*o23/g(3)
      c12=c21*expfn(e21/tev)*g(2)/g(1)
      c13=c31*expfn(e31/tev)*g(3)/g(1)
      c23=c32*expfn(e32/tev)*g(3)/g(2)
       IF((E43/TEV).lt.2.d1) then
      o14=g(1)*o123t4/gsum
      o24=g(2)*o14/g(1)
      o34=g(3)*o14/g(1)
      o25=g(2)*o123t5/gsum
      o15=g(1)*o25/g(2)
      o35=g(3)*o25/g(2)
      c41=cons*o14/g(4)
      c51=cons*o15/g(5)
      c42=cons*o24/g(4)
      c52=cons*o25/g(5)
      c43=cons*o34/g(4)
      c53=cons*o35/g(5)
      c54=cons*o45/g(5)
      c14=c41*expfn(e41/tev)*g(4)/g(1)
      c15=c51*expfn(e51/tev)*g(5)/g(1)
      c24=c42*expfn(e42/tev)*g(4)/g(2)
      c25=c52*expfn(e52/tev)*g(5)/g(2)
      c34=c43*expfn(e43/tev)*g(4)/g(3)
      c35=c53*expfn(e53/tev)*g(5)/g(3)
      c45=c54*expfn(e54/tev)*g(5)/g(4)
       ENDIF
c
      amass=ama(i)
      m=0
 9988 continue
      m=m+1
      adamp=a21(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(2)*g(1))/(x(1)*g(2)))
      call escacf(e21,a21(i),adamp,g(1),g(2),te,amass,tcorr,
     &coltoti,coltot,p21)
      adamp=a31(i)+a32(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(3)*g(1))/(x(1)*g(3)))
      call escacf(e31,a31(i),adamp,g(1),g(3),te,amass,tcorr,
     &coltoti,coltot,p31)
      adamp=a21(i)+a32(i)+a31(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(3)*g(2))/(x(2)*g(3)))
      call escacf(e32,a32(i),adamp,g(2),g(3),te,amass,tcorr,
     &coltoti,coltot,p32)
       IF((E43/TEV).lt.2.d1) then
      adamp=a41(i)+a42(i)+a43(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(4)*g(1))/(x(1)*g(4)))
      call escacf(e41,a41(i),adamp,g(1),g(4),te,amass,tcorr,
     &coltoti,coltot,p41)
      adamp=a51(i)+a52(i)+a54(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(5)*g(1))/(x(1)*g(5)))
      call escacf(e51,a51(i),adamp,g(1),g(5),te,amass,tcorr,
     &coltoti,coltot,p51)
      adamp=a21(i)+a41(i)+a42(i)+a43(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(4)*g(2))/(x(2)*g(4)))
      call escacf(e42,a42(i),adamp,g(2),g(4),te,amass,tcorr,
     &coltoti,coltot,p42)
      adamp=a21(i)+a52(i)+a54(i)+a51(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(5)*g(2))/(x(2)*g(5)))
      call escacf(e52,a52(i),adamp,g(2),g(5),te,amass,tcorr,
     &coltoti,coltot,p52)
      adamp=a31(i)+a32(i)+a41(i)+a42(i)+a43(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(3)*(1.d0-(x(4)*g(3))/(x(3)*g(4)))
      call escacf(e43,a43(i),adamp,g(3),g(4),te,amass,tcorr,
     &coltoti,coltot,p43)
      adamp=a52(i)+a54(i)+a51(i)+a41(i)+a42(i)+a43(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(4)*(1.d0-(x(5)*g(4))/(x(4)*g(5)))
      call escacf(e54,a54(i),adamp,g(4),g(5),te,amass,tcorr,
     &coltoti,coltot,p54)
c
      a(1,1)=c12+c13+c14+c15
      a(1,2)=-(c21+p21*a21(i))
      a(1,3)=-(c31+p31*a31(i))
      a(1,4)=-(c41+p41*a41(i))
      a(1,5)=-(c51+p51*a51(i))
      a(2,1)=-c12
      a(2,2)=c21+c23+c24+c25+p21*a21(i)
      a(2,3)=-(c32+p32*a32(i))
      a(2,4)=-(c42+p42*a42(i))
      a(2,5)=-(c52+p52*a52(i))
      do 100 j=1,5
      x(j)=0.d0
      a(3,j)=1.d0
  100 continue
      x(3)=1.d0
      a(4,1)=-c14
      a(4,2)=-c24
      a(4,3)=-c34
      a(4,4)=c41+c42+c43+c45+p41*a41(i)+p42*a42(i)+p43*a43(i)
      a(4,5)=-(c54+p54*a54(i))
      a(5,1)=-c15
      a(5,2)=-c25
      a(5,3)=-c35
      a(5,4)=-c45
      a(5,5)=c51+c52+c53+c54+p51*a51(i)+p52*a52(i)+p54*a54(i)
c
      call matrix(5)
       ELSE
      a(1,1)=c12+c13
      a(1,2)=-(c21+p21*a21(i))
      a(1,3)=-(c31+p31*a31(i))
      a(2,1)=-c12
      a(2,2)=c21+c23+p21*a21(i)
      a(2,3)=-(c32+p32*a32(i))
      a(3,1)=1.d0
      a(3,2)=1.d0
      a(3,3)=1.d0
      x(1)=0.0d0
      x(2)=0.0d0
      x(3)=1.d0
c
      call matrix(3)
       ENDIF
      if(m.lt.2) goto 9988
c
      elab=(1.6022d-12*abund)/(12.57d0*dene)
       IF((E43/TEV).lt.2.d1) then
      cool=elab*(e21*(x(1)*c12-x(2)*c21)+e31*
     &(x(1)*c13-x(3)*c31)+e41*(x(1)*c14-x(4)*c41)+e51*(x(1)*c15-
     &x(5)*c51)+e32*(x(2)*c23-x(3)*c32)+
     &e42*(x(2)*c24-x(4)*c42)+e52*(x(2)*c25-x(5)*c52)+
     &e43*(x(3)*c34-x(4)*c43)+e53*(x(3)*c35-x(5)*c53)+
     &e54*(x(4)*c45-x(5)*c54))
      f4t123=elab*x(4)*(e41*p41*a41(i)+e42*p42*a42(i)+
     &e43*p43*a43(i))
      if(i.eq.2) then
        f4t123=elab*x(4)*e41*p41*a41(i)
      endif
      f5t12=elab*x(5)*(e51*p51*a51(i)+e52*p52*a52(i))
      f54=elab*x(5)*e54*p54*a54(i)
       ELSE
      cool=elab*(e21*(x(1)*c12-x(2)*c21)+e31*(x(1)*c13-x(3)*c31)+
     &e32*(x(2)*c23-x(3)*c32))
       ENDIF
      cool=12.57d0*cool
      f21=elab*x(2)*p21*a21(i)*e21
      f32=elab*x(3)*p32*a32(i)*e32
c
      endif
C
      RETURN
      END
C*********************************************************************
C*********************************************************************
      SUBROUTINE KFORB(TE,I,DENE,ABUND,COLTOTI,COLTOT,F21,F31,F41,
     &F51,F45T23,F6T45,COOL)
C*********************************************************************
C  Lines (K I-like). Enumeration of the ions are: (= I)
C
C  1 = Ca II
C
C  The atom is treated as a 6-level atom.
C  (Returned to RAD when E21/Tev > 20).
C
C  A-values of Ca II from Ali & Kim (Phys.Rev.A 38, 3992, '88) and 
C  Zeippen (A & A 229, 248, '90) for 4s-3d. (Smaller than Osterbrock's 
C  value, 1.3, from '51 by 30 %.) 4p-5s and 4s-4p from Wiese et al.
C  ('69) scaled by the multiplet calculations of Burgess et al.
C  (A & A 300, 627, '95). 3d-4p from Wiese et al. ('69).
C  A31 from Ritter & Eichmann (J.Phys.B 30, L141 '97). Energies
C  are from the ADS NIST Database ('95). Omegas for most transitions
C  from Burgess et al. ('95), except between the sublevels of 3d which 
C  are from Li & McCray (ApJ 405, 730, '93; based on the results of 
C  Shine ('73)).
C*********************************************************************
      implicit real*8(a-h,o-z)
C
      parameter (mmatr=80)
C
      COMMON/MSTUFF/X(mmatr),A(mmatr,mmatr),EPS,NERR
c
      dimension ene2(1),ene3(1),ene4(1),ene5(1),ene6(1),
     &a21(1),a31(1),a32(1),a41(1),a42(1),a51(1),a52(1),a53(1),
     &a64(1),a65(1),om1t23(1),om1t45(1),om16(1),om23(1),om23t45(1),
     &om23t6(1),om45t6(1),om45(1),al1t23(1),al1t45(1),al16(1),
     &al23(1),al23t45(1),al23t6(1),al45t6(1),al45(1),g(6),ama(1)
c
      data ene2/1.6926d0/
      data ene3/1.7001d0/
      data ene4/3.1237d0/
      data ene5/3.1514d0/
      data ene6/6.4687d0/
      data a21/0.97d0/
C     data a31/1.00d0/
      data a31/1.03d0/
      data a41/1.36d8/
      data a51/1.40d8/
      data a32/2.45d-6/
      data a42/9.45d6/
      data a52/9.90d5/
      data a53/8.83d6/
      data a64/8.8d7/
      data a65/1.7d8/
      data om1t23/0./
      data om1t45/0./
      data om16/0./
      data om23/21.1d0/
      data om23t45/0./
      data om23t6/0./
      data om45/9.85d0/
      data om45t6/0./
      data al1t23/0./
      data al1t45/0./
      data al16/0./
      data al23/0.01d0/
      data al23t45/0./
      data al23t6/0./
      data al45t6/0./
      data al45/-0.01d0/
      data g/2,4,6,2,4,2/
      data ama/40.08d0/
C
      cool=0.0d0
      f21=0.0d0
      f31=0.0d0
      f41=0.0d0
      f51=0.0d0
      f45t23=0.0d0
      f6t45=0.0d0
      if(abund.ge.1.d-10.and.te.ge.1.d3) then
c
      tev=te/1.1609d4
      t4=te/1.d4
      cons=8.63d-6*dene/dsqrt(te)
c
      e21=ene2(i)
        IF((E21/TEV).lt.2.d1) THEN
      e31=ene3(i)
      e41=ene4(i)
      e51=ene5(i)
      e61=ene6(i)
      e32=ene3(i)-ene2(i)
      e42=ene4(i)-ene2(i)
      e52=ene5(i)-ene2(i)
      e62=ene6(i)-ene2(i)
      e43=ene4(i)-ene3(i)
      e53=ene5(i)-ene3(i)
      e63=ene6(i)-ene3(i)
      e54=ene5(i)-ene4(i)
      e64=ene6(i)-ene4(i)
      e65=ene6(i)-ene5(i)
C
      if(i.eq.1) then
        o23=om23(i)*t4**al23(i)
        o45=om45(i)*t4**al45(i)
        o1t23=upsil(2,0.1248d0,1.4d0,2.56d0,2.954d0,3.762d0,
     &	            4.339d0,5.023d0,te)
        o16=upsil(2,0.4757d0,1.5d0,2.601d0,1.778d0,1.981d0,
     &  	  2.261d0,3.382d0,te)
        o23t6=upsil(2,0.3509d0,4.4d0,0.9722d0,0.7578d0,0.6744d0,
     &  	  0.5929d0,0.4160d0,te)
        o1t45=upsil(1,0.2311d0,2.9d0,2.196d1,2.953d1,3.235d1,
     &	            3.743d1,3.756d1,te)
        o23t45=upsil(1,0.1063d0,5.6d0,4.342d1,3.318d1,2.932d1,
     &	            2.786d1,2.418d1,te)
        o45t6=upsil(1,0.2446d0,2.6d0,7.667d0,5.636d0,7.14d0,
     &	            1.149d1,1.708d1,te)
      endif
      gsum1=g(2)+g(3)
      gsum2=g(4)+g(5)
      gmult=gsum1*gsum2
      o12=(g(2)*o1t23)/gsum1
      o13=(g(3)*o1t23)/gsum1
      o14=(g(4)*o1t45)/gsum2
      o15=(g(5)*o1t45)/gsum2
      o24=(g(2)*g(4)*o23t45)/gmult
      o25=(o24*g(5))/g(4)
      o26=(g(2)*o23t6)/gsum1
      o34=(g(3)*g(4)*o23t45)/gmult
      o35=(o34*g(5))/g(4)
      o36=(g(3)*o23t6)/gsum1
      o46=(g(4)*o45t6)/gsum2
      o56=(g(5)*o45t6)/gsum2
      c21=cons*o12/g(2)
      c31=cons*o13/g(3)
      c41=cons*o14/g(4)
      c51=cons*o15/g(5)
      c61=cons*o16/g(6)
      c32=cons*o23/g(3)
      c42=cons*o24/g(4)
      c52=cons*o25/g(5)
      c62=cons*o26/g(6)
      c43=cons*o34/g(4)
      c53=cons*o35/g(5)
      c63=cons*o36/g(6)
      c54=cons*o45/g(5)
      c64=cons*o46/g(6)
      c65=cons*o56/g(6)
      c12=c21*expfn(e21/tev)*g(2)/g(1)
      c13=c31*expfn(e31/tev)*g(3)/g(1)
      c14=c41*expfn(e41/tev)*g(4)/g(1)
      c15=c51*expfn(e51/tev)*g(5)/g(1)
      c16=c61*expfn(e61/tev)*g(6)/g(1)
      c23=c32*expfn(e32/tev)*g(3)/g(2)
      c24=c42*expfn(e42/tev)*g(4)/g(2)
      c25=c52*expfn(e52/tev)*g(5)/g(2)
      c26=c62*expfn(e62/tev)*g(6)/g(2)
      c34=c43*expfn(e43/tev)*g(4)/g(3)
      c35=c53*expfn(e53/tev)*g(5)/g(3)
      c36=c63*expfn(e63/tev)*g(6)/g(3)
      c45=c54*expfn(e54/tev)*g(5)/g(4)
      c46=c64*expfn(e64/tev)*g(6)/g(4)
      c56=c65*expfn(e65/tev)*g(6)/g(5)
c
      amass=ama(i)
      m=0
 9988 continue
      m=m+1
      adamp=a21(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(2)*g(1))/(x(1)*g(2)))
      call escacf(e21,a21(i),adamp,g(1),g(2),te,amass,tcorr,
     &coltoti,coltot,p21)
      adamp=a31(i)+a32(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(3)*g(1))/(x(1)*g(3)))
      call escacf(e31,a31(i),adamp,g(1),g(3),te,amass,tcorr,
     &coltoti,coltot,p31)
      adamp=a41(i)+a42(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(4)*g(1))/(x(1)*g(4)))
      call escacf(e41,a41(i),adamp,g(1),g(4),te,amass,tcorr,
     &coltoti,coltot,p41)
      adamp=a51(i)+a52(i)+a53(i)
      if(m.eq.1) tcorr=1.d0
      if(m.ne.1) tcorr=x(1)*(1.d0-(x(5)*g(1))/(x(1)*g(5)))
      call escacf(e51,a51(i),adamp,g(1),g(5),te,amass,tcorr,
     &coltoti,coltot,p51)
      adamp=a31(i)+a32(i)+a21(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(3)*g(2))/(x(2)*g(3)))
      call escacf(e32,a32(i),adamp,g(2),g(3),te,amass,tcorr,
     &coltoti,coltot,p32)
      adamp=a21(i)+a41(i)+a42(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(4)*g(2))/(x(2)*g(4)))
      call escacf(e42,a42(i),adamp,g(2),g(4),te,amass,tcorr,
     &coltoti,coltot,p42)
      adamp=a21(i)+a51(i)+a52(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(2)*(1.d0-(x(5)*g(2))/(x(2)*g(5)))
      call escacf(e52,a52(i),adamp,g(2),g(5),te,amass,tcorr,
     &coltoti,coltot,p52)
      adamp=a31(i)+a32(i)+a51(i)+a52(i)+a53(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(3)*(1.d0-(x(5)*g(3))/(x(3)*g(5)))
      call escacf(e53,a53(i),adamp,g(3),g(5),te,amass,tcorr,
     &coltoti,coltot,p53)
      adamp=a41(i)+a42(i)+a64(i)+a65(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(4)*(1.d0-(x(6)*g(4))/(x(4)*g(6)))
      call escacf(e64,a64(i),adamp,g(4),g(6),te,amass,tcorr,
     &coltoti,coltot,p64)
      adamp=a51(i)+a52(i)+a53(i)+a64(i)+a65(i)
      if(m.eq.1) tcorr=0.d0
      if(m.ne.1) tcorr=x(5)*(1.d0-(x(6)*g(5))/(x(5)*g(6)))
      call escacf(e65,a65(i),adamp,g(5),g(6),te,amass,tcorr,
     &coltoti,coltot,p65)
      a(1,1)=c12+c13+c14+c15+c16
      a(1,2)=-(c21+p21*a31(i))
      a(1,3)=-(c31+p31*a31(i))
      a(1,4)=-(c41+p41*a41(i))
      a(1,5)=-(c51+p51*a51(i))
      a(1,6)=-c61
      a(2,1)=-c12
      a(2,2)=c21+c23+c24+c25+c26+p21*a21(i)
      a(2,3)=-(c32+p32*a32(i))
      a(2,4)=-(c42+p42*a42(i))
      a(2,5)=-(c52+p52*a52(i))
      a(2,6)=-c62
      do 100 j=1,6
      x(j)=0.d0
      a(3,j)=1.d0
  100 continue
      x(3)=1.d0
      a(4,1)=-c14
      a(4,2)=-c24
      a(4,3)=-c34
      a(4,4)=c41+c42+c43+c45+c46+p41*a41(i)+p42*a42(i)
      a(4,5)=-c54
      a(4,6)=-(c64+p64*a64(i))
      a(5,1)=-c15
      a(5,2)=-c25
      a(5,3)=-c35
      a(5,4)=-c45
      a(5,5)=c51+c52+c53+c54+c56+p51*a51(i)+p52*a52(i)+
     &       p53*a53(i)
      a(5,6)=-(c65+p65*a65(i))
      a(6,1)=-c16
      a(6,2)=-c26
      a(6,3)=-c36
      a(6,4)=-c46
      a(6,5)=-c56
      a(6,6)=c61+c62+c63+c64+c64+p64*a64(i)+p65*a65(i)
c
      call matrix(6)
      if(m.lt.2) goto 9988
c
      elab=(1.6022d-12*abund)/(12.57d0*dene)
      cool=elab*(e21*(x(1)*c12-x(2)*c21)+e31*
     &(x(1)*c13-x(3)*c31)+e41*(x(1)*c14-x(4)*c41)+e51*(x(1)*c15-
     &x(5)*c51)+e61*(x(1)*c16-x(6)*c61)+e32*(x(2)*c23-x(3)*c32)+
     &e42*(x(2)*c24-x(4)*c42)+e52*(x(2)*c25-x(5)*c52)+e62*(x(2)*
     &c26-x(6)*c62)+e43*(x(3)*c34-x(4)*c43)+e53*(x(3)*c35-x(5)*
     &c53)+e63*(x(3)*c36-x(6)*c63)+e54*(x(4)*c45-x(5)*c54)+
     &e64*(x(4)*c46-x(6)*c64)+e65*(x(5)*c56-x(6)*c65))
      cool=12.57d0*cool
      f21=elab*x(2)*a21(i)*e21*p21
      f31=elab*x(3)*a31(i)*e31*p31
      f41=elab*x(4)*a41(i)*e41*p41
      f51=elab*x(5)*a51(i)*e51*p51
      f45t23=elab*(x(4)*e42*a42(i)*p42+x(5)*(e52*a52(i)*p52+
     &e53*a53(i)*p53))
      f6t45=elab*x(6)*(p64*e64*a64(i)+p65*e65*a65(i))
c
        ENDIF
      endif
C
      RETURN
      END
C*********************************************************************
c*********************************************************************   
      double precision function spline(p1,p2,p3,p4,p5,x)
C
C     Devised by Burgess & Tully (A & A 254, 436, '92).
C************************************************************           
      implicit real*8(a-h,o-z)
C
      s=1.d0/30.d0
      s2=32.d0*s*(19.d0*p1-43.d0*p2+30.d0*p3-7.d0*p4+p5)
      s3=160.d0*s*(-p1+7.d0*p2-12.d0*p3+7.d0*p4-p5)
      s4=32.d0*s*(p1-7.d0*p2+30.d0*p3-43.d0*p4+19.d0*p5)
      if(x.gt.0.25d0) goto 1
      x0=x-0.125d0
      t3=0.0d0
      t2=0.5d0*s2
      t1=4.d0*(p2-p1)
      t0=0.5d0*(p1+p2)-0.015625d0*t2
      goto 4
    1 if(x.gt.0.5d0) goto 2
      x0=x-0.375d0
      t3=20.d0*s*(s3-s2)
      t2=0.25d0*(s2+s3)
      t1=4.0d0*(p3-p2)-0.015625d0*t3
      t0=0.5d0*(p2+p3)-0.015625d0*t2
      goto 4
    2 if(x.gt.0.75d0) goto 3
      x0=x-0.625d0
      t3=20.d0*s*(s4-s3)
      t2=0.25d0*(s4+s3)
      t1=4.0d0*(p4-p3)-0.015625d0*t3
      t0=0.5d0*(p4+p3)-0.015625d0*t2
      goto 4
    3 x0=x-0.875d0
      t3=0.0d0
      t2=0.5d0*s4
      t1=4.d0*(p5-p4)
      t0=0.5d0*(p5+p4)-0.015625d0*t2
    4 spline=t0+x0*(t1+x0*(t2+x0*t3))
C  
      return
      end
c*********************************************************************   
c*********************************************************************   
      double precision function upsil(k,eij,c,p1,p2,p3,p4,p5,t)
C
C     Devised by Burgess & Tully (A & A 254, 436, '92).
C************************************************************           
      implicit real*8(a-h,o-z)
C
      teij=1.57888d5*eij
      e=dabs(t/teij)
      if(k.eq.1.or.k.eq.4) x=dlog((e+c)/c)/dlog(e+c)
      if(k.eq.2.or.k.eq.3) x=e/(e+c)
      y=spline(p1,p2,p3,p4,p5,x)
      if(k.eq.1) y=y*dlog(e+2.71828d0)
      if(k.eq.3) y=y/(e+1)
      if(k.eq.4) y=y*dlog(e+c)
      upsil=y
C
      return
      end
c*********************************************************************   
C*********************************************************************
      double precision function cheb(tmx,t,a0,a1,a2,a3)
      implicit real*8(a-h,o-z)
C
      y=t
      if(y.lt.5.d3) y=5.d3
      if(y.gt.tmx) y=tmx
      x=0.6199646d0*dlog(y)-6.2803580d0
      t1=x
      t2=2.d0*x**2-1.d0
      t3=x*(4.d0*x**2-3.d0)
      omlog=0.5d0*a0+a1*t1+a2*t2+a3*t3
      cheb=expfn(-omlog)
C      cheb=dexp(omlog)
C
      return
      end
C*********************************************************************
C*********************************************************************
      double precision function cheb1(tmx,t,a0,a1,a2,a3)
      implicit real*8(a-h,o-z)
C
      y=t
      if(y.lt.5.d3) y=5.d3
      if(y.gt.tmx) y=tmx
      x=0.542519d0*dlog(y)-5.620738d0
      t1=x
      t2=2.d0*x**2-1.d0
      t3=x*(4.d0*x**2-3.d0)
      omlog=0.5d0*a0+a1*t1+a2*t2+a3*t3
      cheb1=dexp(omlog)
C
      return
      end
C*********************************************************************


