      program fivel
c
c     /////////////////////////////////////////////////////////////////
c     ////     The LICK 5-level atom program - FIVEL               ////
c     ////                                                         ////
c     ////     Language: FORTRAN 77                                ////
c     ////                                                         ////
c     ////     Certain astrophysically common ions have emission-  ////
c     ////     line intensity ratios which are sensitive to either ////
c     ////     temperature (T) or number density (N). From an      ////
c     ////     observed line ratio (ROBS), then, given a tem-      ////
c     ////     perature, the density can be calculated; or, given  ////
c     ////     the density, the temperature can be calculated.     ////
c     ////     Similarly, given the temperature and density,       ////
c     ////     the emission-line emissivities can be calculated.   ////
c     ////                                                         ////
c     ////     Program Structure:                                  ////
c     ////     1. "Physical conditions"                            ////
c     ////         (a) Specify 'T' option to use observed line     ////
c     ////              ratio and density to find temperature      ////
c     ////         (b) Specify 'N' option to use observed line     ////
c     ////             ratio and temperature to find density       ////
c     ////                                                         ////
c     ////     2. "Line emissivities"                              ////
c     ////         Supply temperature and density to determine     ////
c     ////         the level populations, critical densities,      ////
c     ////         and line emissivities                           ////
c     ////                                                         ////
c     ////     MAIN Program:                                       ////
c     ////          Select option 1 or 2 and ion to investigate.   ////
c     ////                                                         ////
c     ////     Subroutines:                                        ////
c     ////        INPUT:  Supplies data for calculation            ////
c     ////        OUTPUT: Displays final results                   ////
c     ////        PARI:   Contains atomic data independent of T    ////
c     ////        PARII:  Contains atomic data dependent on T      ////
c     ////        SOLVE:  Performs matrix algebra. Finds level     ////
c     ////                populations (POP) and line emissivities  ////
c     ////                (JL) as a function of N and T, and       ////
c     ////                critical densities (NCRIT).              ////
c     ////        TEMDEN: Employs an iterative technique (geo-     ////
c     ////                metric bisection) to find 'T' or 'N'.    ////
c     ////                Used only when T or N is unknown.        ////
c     ////                                                         ////
c     ////     Written May 1987/ M. M. De Robertis, R. J. Dufour,  ////
c     ////     and R. W. Hunt.                                     ////
c     ////---------------------------------------------------------////
c     //// ** The default logical units may have to be changed de- ////
c     //// pending on the computer/operating system. Currently,    ////
c     //// unit '*' is used for input and "interactive output",    ////
c     //// while '6' is used for writing out the calculations.     ////
c     ////                                                         ////
c     /////////////////////////////////////////////////////////////////
c
c     _________________________________________________________________
c     |                                                               |
c     |                        MAIN PROGRAM                           |
c     |                                                               |
c     |    IONUM:   integer code for ion of interest                  |
c     |    ICTG:    category (population or physical conditions)      |
c     |    PHYG:    logical for available "physical condition" ions   |
c     |    POPG:    logical for the available population ions         |
c     |    A:       transition probability matrix                     |
c     |    C:       collision strength matrix                         |
c     |    Q:       collision transition probability matrix           |
c     |    JL:      volume emissivity matrix                          |
c     |    JHB:     H-beta volume emissivity                          |
c     |    NCRIT:   critical density of upper level                   |
c     |    POP:     relative level populations                        |
c     |    WEIGHT:  statistical weights of levels                     |
c     |    T:       energy above ground state matrix (1/Angstroms)    |
c     -----------------------------------------------------------------
c
      real a(5,5),c(5,5),q(5,5),jl(5,5),t(5),pop(5)
      integer kbin,ttyout,filout
      logical phyg,popg

      common /tty/ kbin,ttyout,filout
      data kbin/5/,ttyout/6/,filout/1/

      open(filout,file='fiveout.dat')
      rewind filout

      write(ttyout,10)
  100 write(ttyout,15)
      read(kbin,*,end=100,err=100) ictg
      write(ttyout,*) ' '

c
c Physical conditions
c
      if(ictg .eq. 1) then
c
c Ion code: Atomic number//ionization stage//N(=1),T(=2)
c
        write(ttyout,*)'LICK 5-LEVEL ATOM PHYSICAL CONDITIONS CATEGORY'
        write(ttyout,*)'       *** AVAILABLE IONS ***'
        write(ttyout,*)'   7001: [N I]    N;   7012: [N II]   T;'
        write(ttyout,*)'   8002: [O I]    T;   8011: [O II]   N;'
        write(ttyout,*)'   8012: [O II]   T;   8022: [O III]  T;'
        write(ttyout,*)'  10022: [Ne III] T;  10031: [Ne IV]  N;'
        write(ttyout,*)'  10042: [Ne V]   T;  16011: [S II]   N;'
        write(ttyout,*)'  16012: [S II]   T;  16022: [S III]  T;'
        write(ttyout,*)'  17021: [Cl III] N;  17032: [Cl IV]  T;'
        write(ttyout,*)'  18022: [Ar III] T;  18031: [Ar IV]  N;'
        write(ttyout,*)'  18032: [Ar IV]  T;  18042: [Ar V]   T;'
  110   write(ttyout,20)
        read(kbin,*,end=110,err=110) ionum
c
c For new ions add appropriate parameters in IONUM and PHYG block
c
        phyg=(ionum .eq. 7001).or.(ionum .eq. 7012).or.
     1                            (ionum .eq. 8011).or.
     2       (ionum .eq. 8012).or.(ionum .eq. 8022).or.
     3       (ionum .eq.10022).or.(ionum .eq.10031).or.
     4       (ionum .eq.10042).or.(ionum .eq.16011).or.
     5       (ionum .eq.16012).or.(ionum .eq.16022).or.
     6       (ionum .eq.17021).or.(ionum .eq.17032).or.
     7       (ionum .eq.18022).or.(ionum .eq.18031).or.
     8       (ionum .eq.18032).or.(ionum .eq.18042)

        if(phyg) then
          write(filout,*)' '
          write(filout,*)'Physical conditions calculations 
     1                    for ion number: ',ionum
          write(filout,*)' '
          call input(ionum,den,tem,robs)
          call pari(ionum)
          call temden(ionum,den,tem,robs,pop,jl)
          if(robs .ne. 0.) call output(ionum,den,tem,pop,jl)
        else
          write(filout,*)'No parameters for that ion; enter another'
        endif

c
c Level populations, line ratios
c
      else if(ictg .eq. 2) then
        popg=.true.
c
c Ion code: Atomic number//ionization stage//0(=emissivities)
c
        write(ttyout,*)'LICK 5-LEVEL ATOM-POPULATIONS CATEGORY'
        write(ttyout,*)'       *** AVAILABLE IONS ***'
        write(ttyout,*)'   7000: [N I];       7010: [N II];'
        write(ttyout,*)'   8000: [O I];       8010: [O II];'
        write(ttyout,*)'   8020: [O III];    10020: [Ne III];'
        write(ttyout,*)'  10030: [Ne IV];    10040: [Ne V];  '
        write(ttyout,*)'  16010: [S II];     16020: [S III]; '
        write(ttyout,*)'  17010: [Cl II];    17020: [Cl III];'
        write(ttyout,*)'  17030: [Cl IV];    18020: [Ar III];'
        write(ttyout,*)'  18030: [Ar IV];    18040: [Ar V];  '
  120   write(ttyout,20)
        read(kbin,*,err=120,end=120) ionum
c
c For new ions, add appropriate parameters in IONUM and POPG block.
c
        phyg=(ionum .eq. 7000).or.(ionum .eq. 7010).or.
     1       (ionum .eq. 8000).or.(ionum .eq. 8010).or.
     2       (ionum .eq. 8020).or.(ionum .eq.10020).or.
     3       (ionum .eq.10030).or.(ionum .eq.10040).or.
     4       (ionum .eq.16010).or.(ionum .eq.16020).or.
     5       (ionum .eq.17010).or.(ionum .eq.17020).or.
     6       (ionum .eq.17030).or.(ionum .eq.18020).or.
     7       (ionum .eq.18030).or.(ionum .eq.18040)

        if(popg) then
          write(filout,*) 'Line emissivity calculations',
     1                    'for ion number ',ionum
          write(filout,*) ' '
          call input(ionum,den,tem,robs)
          call pari(ionum)
          call parii(ionum,tem)
          call solve(den,tem,pop,jl)
          call output(ionum,den,tem,pop,jl)
        else
          write(filout,*) 'No parameters for that ion; try another'
        endif
c
c Stop
c
      else if(ictg .eq. 3) then
        stop
      else
c
c Try again
c 
        write(ttyout,*) 'Enter 1, 2, or 3 only'
      endif
      go to 100
  200 continue

   10 format(//,20x,'<<< LICK 5-LEVEL ATOM PROGRAM >>>',//)
   15 format(/,' Calculate PHYSICAL CONDITIONS(1), LINE RATIOS(2),',
     1' or STOP(3): ')
   20 format(' Select one of the above by typing number: ')
      end
c-----------------------------------------------------------------------
      subroutine input(ionum,den,tem,robs)
c
c     _________________________________________________________________
c     |                                                               |
c     |                        INPUT ROUTINE                          |
c     |                                                               |
c     |   ROBS: the observed line ratio                               |
c     |   DEN:  the electron number density of the gas                |
c     |   TEM:  the electron temperature of the gas                   |
c     |   PHYS: logical to decide to read ROBS                        |
c     |                                                               |
c     -----------------------------------------------------------------
c
      integer kbin,ttyout,filout
      logical phys
      common /tty/ kbin,ttyout,filout

      phys = .true.
c
c  Enter relevant physical parameters as required
c
      if(mod(ionum,10) .eq. 0) then
  100   write(ttyout,10)
        read(kbin,*,err=100,end=100) den,tem
        write(filout,50) int(den),int(tem)
        if(tem .lt. 1.e3 .or. tem .gt. 3.6e4) then
          write(ttyout,*)'Indicated temperature outside permitted range'
          stop
        endif
      else if(mod(ionum,10) .eq. 2) then
  110   write(ttyout,20)
        read(kbin,*,err=110,end=110) den
      else
  120   write(ttyout,30) 
        read(kbin,*,err=120,end=120) tem
        if(tem .lt. 1.e3 .or. tem .gt. 3.6e4) then
          write(ttyout,*)'Indicated temperature outside permitted range'
          stop
        endif
      endif
c
c  For new ion, line ratio must be included in IF block below
c
  130 if(ionum .eq. 7001) then
        write(ttyout,*) '[N I] N calculation'
        write(ttyout,*) 'I(5200)/I(5198) = '
      else if(ionum .eq. 7012) then
        write(ttyout,*) '[N II] T calculation'
        write(ttyout,*) 'I(6548+6583)/I(5755) = '
      else if(ionum .eq. 8011) then
        write(ttyout,*) '[O II] N calculation'
        write(ttyout,*) 'I(3726)/I(3729) = '
      else if(ionum .eq. 8012) then
        write(ttyout,*) '[O II] N calculation'
        write(ttyout,*) 'I(3727)/I(7324) = '
      else if(ionum .eq. 8022) then
        write(ttyout,*) '[O III] T calculation'
        write(ttyout,*) 'I(4959+5007)/I(4363) = '
      else if(ionum .eq. 10022) then
        write(ttyout,*) '[Ne III] T calculation'
        write(ttyout,*) 'I(3869+3968)/I(3343) = '
      else if(ionum .eq. 10031) then
        write(ttyout,*) '[Ne IV] N calculation'
        write(ttyout,*) 'I(2423)/I(2425) = '
      else if(ionum .eq. 10042) then
        write(ttyout,*) '[Ne V] T calculation'
        write(ttyout,*) 'I(3426+3346)/I(2975) = '
      else if(ionum .eq. 16011) then
        write(ttyout,*) '[S II] N calculation'
        write(ttyout,*) 'I(6717)/I(6731) = '
      else if(ionum .eq. 16012) then
        write(ttyout,*) '[S II] T calculation'
        write(ttyout,*) 'I(6717+6731)/I(4069+4076) = '
      else if(ionum .eq. 16022) then
        write(ttyout,*) '[S III] T calculation'
        write(ttyout,*) 'I(9069+9532)/I(6312) = '
      else if(ionum .eq. 17021) then
        write(ttyout,*) '[Cl III] N calculation'
        write(ttyout,*) 'I(5518)/I(5538) = '
      else if(ionum .eq. 17032) then
        write(ttyout,*) '[Cl IV] T calculation'
        write(ttyout,*) 'I(7530+8045)/I(5323) = '
      else if(ionum .eq. 18022) then
        write(ttyout,*) '[Ar III] T calculation'
        write(ttyout,*) 'I(7136+7751)/I(5192) = '
      else if(ionum .eq. 18031) then
        write(ttyout,*) '[Ar IV] N calculation'
        write(ttyout,*) 'I(4711)/I(4740) = '
      else if(ionum .eq. 18032) then
        write(ttyout,*) '[Ar IV] T calculation'
        write(ttyout,*) 'I(4711+4740)/I(2854+2868) = '
      else if(ionum .eq. 18042) then
        write(ttyout,*) '[Ar V] T calculation'
        write(ttyout,*) 'I(6435+7006)/I(4626) = '
      else
        phys = .false.
      endif
      if(phys) then
        read(kbin,*,err=130,end=130) robs
        write(filout,40) robs
      endif
   10 format(' Enter density and temperature to find ratio: ')
   20 format(' Enter fixed density to find temperature: ')
   30 format(' Enter fixed temperature to find density: ')
   40 format(' Observed line ratio = ',f8.2)
   50 format(' Electron density = ', i7,' Electron temperature = ', 
     .       i6)
      return
      end
c----------------------------------------------------------------------
      subroutine output(ionum,den,tem,pop,jl)
c
c     ________________________________________________________________
c     |                                                              |
c     |                           OUTPUT                             |
c     |                                                              |
c     |     POP:      fraction of electrons in each level            |
c     |     JL:       volume emission coefficient                    |
c     |     JHB:      volume emissivity of H-beta (case B)           |
c     |               (using Brocklehurst's interpolation formula    |
c     |     NCRIT:    critical density                               |
c     |     XL,XI:    wavelength (Angstroms), intensity              |
c     |               (N.B.- wavelengths > 1e5 Angstroms in microns) |
c     |                                                              |
c     ----------------------------------------------------------------
c
      real jl(5,5), pop(5), ncrit, xl(5), xi(5), jhb
      integer weight(5)
      integer kbin,ttyout,filout
      common/atom/c(5,5),a(5,5),ncrit(5),t(5),weight
      common/tty/kbin,ttyout,filout
c
c  Write level populations and critical densities
c
      write(filout,*) ' '
      write(filout,10) pop(1),(pop(j),ncrit(j), j=2,5)
      write(filout,*) ' '
c
c  Pause statement for screen convenience
c
   50 if(mod(ionum,10) .eq. 0) then
      write(ttyout,*) 'Hit 0 to continue'
      read(kbin,*,err=50,end=50) dummy
c
c  Transform wavelength into microns if > 100000 Angstroms
c  If emissivity is < 1e-35 erg/s, set it to zero.
c
      do 100 i=2,5
         do 200 j=1,4
            if(j .lt. i) then
              xl(j)=1./(t(i)-t(j))
              if(xl(j) .ge. 1.e5) xl(j)=xl(j)*1.e-4
              xi(j)=jl(i,j)/den
              if(xi(j) .lt. 1.e-35) xi(j)=0.
            endif
  200    continue
c
         if(i .eq. 2) then
           write(filout,11) (xl(k),k=1,i-1),(i,k,k=1,i-1),
     1                      (xi(k),k=1,i-1)
         elseif(i .eq. 3) then
           write(filout,12) (xl(k),k=1,i-1),(i,k,k=1,i-1),
     1                      (xi(k),k=1,i-1)
         elseif(i .eq. 4) then
           write(filout,13) (xl(k),k=1,i-1),(i,k,k=1,i-1),
     1                      (xi(k),k=1,i-1)
         elseif(i .eq. 5) then
           write(filout,14) (xl(k),k=1,i-1),(i,k,k=1,i-1),
     1                      (xi(k),k=1,i-1)
         endif
  100 continue
c
c  Brocklehurst's interpolation formula for H-beta emissivity
c  Accurate for DEN <= 1.e4/cm**3
c
      t4=tem/1.e4
      e42=1.387/t4**0.983/10.**(0.0424/t4)
      jhb=e42*10.**-25
      write(filout,15) jhb
c
c  Output final results for TEM/DEN ions.
c
      elseif(ionum .eq. 7001) then
        write(filout,*) 'DENS(N0) = ',den
      elseif(ionum .eq. 7012) then
        write(filout,*) 'TEMP(N+) = ',tem
      elseif(ionum .eq. 8011) then
        write(filout,*) 'DENS(O+) = ',den
      elseif(ionum .eq. 8012) then
        write(filout,*) 'TEMP(O+) = ',tem
      elseif(ionum .eq. 8022) then
        write(filout,*) 'TEMP(O++) = ',tem
      elseif(ionum .eq. 10022) then
        write(filout,*) 'TEMP(Ne++) = ',tem
      elseif(ionum .eq. 10031) then
        write(filout,*) 'DENS(Ne+3) = ',den
      elseif(ionum .eq. 10042) then
        write(filout,*) 'TEMP(Ne+4) = ',tem
      elseif(ionum .eq. 16011) then
        write(filout,*) 'DENS(S+) = ',den
      elseif(ionum .eq. 16012) then
        write(filout,*) 'TEMP(S+) = ',tem
      elseif(ionum .eq. 16022) then
        write(filout,*) 'TEMP(S++) = ',tem
      elseif(ionum .eq. 17021) then
        write(filout,*) 'DENS(Cl++) = ',den
      elseif(ionum .eq. 17032) then
        write(filout,*) 'TEMP(Cl+3) = ',tem
      elseif(ionum .eq. 18022) then
        write(filout,*) 'TEMP(Ar++) = ',tem
      elseif(ionum .eq. 18031) then
        write(filout,*) 'DENS(Ar+3) = ',den
      elseif(ionum .eq. 18042) then
        write(filout,*) 'TEMP(Ar+4) = ',tem
      endif
c
      write(filout,*) 'Log10 x = ',alog10(den/sqrt(tem))
c
   10 format(10x,'Level Populations - Critical Densities (cm**-3):',/,
     110x,'Level 1: ',1pe8.2,/,
     210x,'Level 2: ',1pe8.2,13x,1pe8.2,/,
     310x,'Level 3: ',1pe8.2,13x,1pe8.2,/,
     410x,'Level 4: ',1pe8.2,13x,1pe8.2,/,
     510x,'Level 5: ',1pe8.2,13x,1pe8.2)
   11 format(5x,'Wavelength',5x,f11.2,/,
     1       1x,'Upper--Lower Levels',5x,'(',i1,'--',i1,') ',/,
     2       2x,'Volume Emissivity ',1pe11.3,/)
   12 format(20x,2f11.2,/,20x,2(5x,'(',i1,'--',i1,')'),/,
     1       20x,2(1pe11.3),/)
   13 format(20x,3f11.2,/,20x,3(5x,'(',i1,'--',i1,')'),/,
     1       20x,3(1pe11.3),/)
   14 format(20x,4f11.2,/,20x,4(5x,'(',i1,'--',i1,')'),/,
     1       20x,4(1pe11.3),/)
   15 format(1x,'H-beta volume emissivity: ',1pe9.2,' N(H+)Ne',
     1       ' erg/s',/)
      return
      end
c-----------------------------------------------------------------------
      subroutine pari(ionum)
c
c     _________________________________________________________________
c     |                                                               |     
c     |                           PAR I                               |
c     |                                                               |
c     |  T:      energy-level separation from ground in 1/Angstroms   |
c     |  WEIGHT: statistical weights for each level                   |
c     |  A:      radiative transition probabilities                   |
c     |                                                               |
c     -----------------------------------------------------------------

      real ncrit
      integer weight(5)
      integer kbin,ttyout,filout
      common /atom/ c(5,5),a(5,5),ncrit(5),t(5),weight
      common /tty/ kbin,ttyout,filout

c
c  Zero arrays
c
      do 100 i=1,5
         do 200 j=1,5
            a(i,j)=0.0
            c(i,j)=0.0
  200    continue
  100 continue
c
      t(1)=0.0
      if((ionum .eq. 7000) .or. (ionum .eq. 7001)
     1                     .or. (ionum .eq. 7002)) then
c
c  N I ion parameters
c
        a(2,1)=7.27e-6
        a(3,1)=2.02e-5
        a(4,1)=2.71e-3
        a(5,1)=6.58e-3
        a(3,2)=1.27e-8
        a(4,2)=3.45e-2
        a(5,2)=6.14e-2
        a(4,3)=5.29e-2
        a(5,3)=2.76e-2
        t(2)=1.92245e-4
        t(3)=1.92332e-4
        t(4)=2.88389e-4
        t(5)=2.88393e-4
        weight(1)=4
        weight(2)=6
        weight(3)=4
        weight(4)=2
        weight(5)=4
      else if((ionum .eq. 7010) .or. (ionum .eq. 7011)
     1                          .or. (ionum .eq. 7012)) then
c
c  N II ion parameters
c
        a(2,1)=2.08e-6
        a(3,1)=1.16e-12
        a(4,1)=5.35e-7
        a(3,2)=7.46e-6
        a(4,2)=1.01e-3
        a(5,2)=3.38e-2
        a(4,3)=2.99e-3
        a(5,3)=1.51e-4
        a(5,4)=1.12   
        t(2)=4.87e-7
        t(3)=1.308e-6
        t(4)=1.53162e-4
        t(5)=3.26888e-4
        weight(1)=1
        weight(2)=3
        weight(3)=5
        weight(4)=5
        weight(5)=1
      else if((ionum .eq. 8000) .or. (ionum .eq. 8001)
     1                          .or. (ionum .eq. 8002)) then
c
c  O I ion parameters
c
        a(2,1)=8.92e-5
        a(4,1)=6.34e-3
        a(5,1)=2.88e-4
        a(3,2)=1.74e-5
        a(4,2)=2.11e-3
        a(5,2)=7.32e-2
        a(4,3)=7.23e-7
        a(5,4)=1.22   
        t(2)=1.583e-6
        t(3)=2.27e-6
        t(4)=1.58679e-4
        t(5)=3.37926e-4
        weight(1)=5
        weight(2)=3
        weight(3)=1
        weight(4)=5
        weight(5)=1

      else if((ionum .eq. 8010) .or. (ionum .eq. 8011)
     1                          .or. (ionum .eq. 8012)) then
c
c  O II ion parameters
c
        a(2,1)=3.82e-5
        a(3,1)=1.65e-4
        a(4,1)=5.64e-2
        a(5,1)=2.32e-2
        a(3,2)=1.2e-7
        a(4,2)=1.17e-1
        a(5,2)=6.15e-2
        a(4,3)=6.14e-2
        a(5,3)=1.02e-1
        a(5,4)=2.08e-11
        t(2)=2.68107e-4
        t(3)=2.68302e-4
        t(4)=4.04675e-4
        t(5)=4.04686e-4
        weight(1)=4
        weight(2)=6
        weight(3)=4
        weight(4)=4
        weight(5)=2
c
      else if((ionum .eq. 8020) .or. (ionum .eq. 8021)
     1                          .or. (ionum .eq. 8022)) then
c
c  O III ion parameters
c
        a(2,1)=2.62e-5
        a(3,1)=3.02e-11
        a(4,1)=2.74e-6
        a(3,2)=9.76e-5
        a(4,2)=6.74e-3
        a(5,2)=2.23e-1
        a(4,3)=1.96e-2
        a(5,3)=7.85e-4
        a(5,4)=1.78
        t(2)=1.132e-6
        t(3)=3.062e-6
        t(4)=2.02733e-4
        t(5)=4.318577e-4
        weight(1)=1
        weight(2)=3
        weight(3)=5
        weight(4)=5
        weight(5)=1
c
      else if((ionum .eq. 10020) .or. (ionum .eq. 10021)
     1                           .or. (ionum .eq. 10022)) then
c
c  Ne III ion parameters
c
        a(2,1)=5.97e-3
        a(3,1)=2.18e-8
        a(4,1)=1.71e-1
        a(5,1)=3.94e-3
        a(3,2)=1.15e-3
        a(4,2)=5.42e-2
        a(5,2)=2.0
        a(4,3)=8.51e-6
        a(5,4)=2.71
        t(2)=6.429e-6
        t(3)=9.205e-6
        t(4)=2.58408e-4
        t(5)=5.57506e-4
        weight(1)=5
        weight(2)=3
        weight(3)=1
        weight(4)=5
        weight(5)=1
c
      else if((ionum .eq. 10030) .or. (ionum .eq. 10031)
     1                           .or. (ionum .eq. 10032)) then
c
c  Ne IV ion parameters
c
        a(2,1)=4.84e-4
        a(3,1)=5.54e-3
        a(4,1)=5.21e-1
        a(5,1)=1.27
        a(3,2)=1.48e-6
        a(4,2)=1.15e-1
        a(5,2)=4.0e-1
        a(4,3)=3.93e-1
        a(5,3)=4.37e-1
        a(5,4)=2.68e-9
        t(2)=4.12346e-4
        t(3)=4.12795e-4
        t(4)=6.24346e-4
        t(5)=6.24413e-4
        weight(1)=4
        weight(2)=6
        weight(3)=4
        weight(4)=2
        weight(5)=4
c
      else if((ionum .eq. 10040) .or. (ionum .eq. 10041)
     1                           .or. (ionum .eq. 10042)) then
c
c  Ne V ion parameters
c
        a(2,1)=1.28e-3
        a(3,1)=5.08e-9
        a(4,1)=2.37e-5
        a(3,2)=4.59e-3
        a(4,2)=1.31e-1
        a(5,2)=4.21
        a(4,3)=3.65e-1
        a(5,3)=6.69e-3
        a(5,4)=2.85
        t(2)=4.124e-6
        t(3)=1.1101e-5
        t(4)=3.02915e-4
        t(5)=6.39136e-4
        weight(1)=1
        weight(2)=3
        weight(3)=5
        weight(4)=5
        weight(5)=1
c
      else if((ionum .eq. 16010) .or. (ionum .eq. 16011)
     1                           .or. (ionum .eq. 16012)) then
c
c  S II ion parameters
c
        a(2,1)=8.82e-4
        a(3,1)=2.60e-4
        a(4,1)=9.06e-2
        a(5,1)=2.25e-1
        a(3,2)=3.35e-7
        a(4,2)=1.63e-1
        a(5,2)=1.33e-1
        a(4,3)=7.79e-2
        a(5,3)=1.79e-1
        a(5,4)=1.03e-6
        t(2)=1.48530e-4
        t(3)=1.48848e-4
        t(4)=2.45249e-4
        t(5)=2.45718e-4
        weight(1)=4
        weight(2)=4
        weight(3)=6
        weight(4)=2
        weight(5)=4
c
      else if((ionum .eq. 16020) .or. (ionum .eq. 16021)
     1                           .or. (ionum .eq. 16022)) then
c
c  S III ion parameters
c
        a(2,1)=4.72e-4
        a(3,1)=4.61e-8
        a(4,1)=5.82e-6
        a(3,2)=2.07e-3
        a(4,2)=2.21e-2
        a(5,2)=7.96e-1
        a(4,3)=5.76e-2
        a(5,3)=1.05e-2
        a(5,4)=2.22
        t(2)=2.972e-6
        t(3)=8.325e-6
        t(4)=1.1320e-4
        t(5)=2.7163e-4
        weight(1)=1
        weight(2)=3
        weight(3)=5
        weight(4)=5
        weight(5)=1
c
      else if((ionum .eq. 17010) .or. (ionum .eq. 17011)
     1                           .or. (ionum .eq. 17012)) then
c
c  Cl II ion parameters
c
        a(2,1)=7.57e-3
        a(3,1)=4.57e-7
        a(4,1)=1.04e-1
        a(5,1)=1.97e-2
        a(3,2)=1.46e-3
        a(4,2)=2.92e-2
        a(5,2)=1.31
        a(4,3)=9.82e-6
        a(5,4)=2.06
        t(2)=6.96e-6
        t(3)=9.965e-6
        t(4)=1.06571e-4
        t(5)=2.78780e-4
        weight(1)=5
        weight(2)=3
        weight(3)=1
        weight(4)=5
        weight(5)=1
c
      else if((ionum .eq. 17020) .or. (ionum .eq. 17021)
     1                           .or. (ionum .eq. 17022)) then
c
c  Cl III ion parameters
c 
        a(2,1)=4.83e-3
        a(3,1)=7.04e-4
        a(4,1)=3.05e-1
        a(5,1)=7.54e-1
        a(3,2)=3.22e-6
        a(4,2)=3.03e-1
        a(5,2)=3.23e-1
        a(4,3)=1.0e-1
        a(5,3)=3.16e-1
        a(5,4)=7.65e-6
        t(2)=1.8053e-4
        t(3)=1.81186e-4
        t(4)=2.9812e-4
        t(5)=2.9907e-4
        weight(1)=5
        weight(2)=3
        weight(3)=1
        weight(4)=5
        weight(5)=1
c
      else if((ionum .eq. 17030) .or. (ionum .eq. 17031)
     1                           .or. (ionum .eq. 17032)) then
c
c  Cl IV ion parameters
c
        a(2,1)=2.14e-3
        a(3,1)=2.70e-7
        a(4,1)=1.54e-5
        a(3,2)=8.25e-3
        a(4,2)=7.23e-2
        a(5,2)=2.47
        a(4,3)=1.79e-1
        a(5,3)=2.62e-2
        a(5,4)=2.80
        t(2)=4.92e-6
        t(3)=1.3419e-5
        t(4)=1.37676e-4
        t(5)=3.2547e-4
        weight(1)=1
        weight(2)=3
        weight(3)=5
        weight(4)=5
        weight(5)=1
c
      else if((ionum .eq. 18020) .or. (ionum .eq. 18021)
     1                           .or. (ionum .eq. 18022)) then
c
c  Ar III ion parameters
c
        a(2,1)=3.08e-2
        a(3,1)=2.37e-6
        a(4,1)=3.14e-1
        a(5,1)=4.17e-2
        a(3,2)=5.17e-3
        a(4,2)=8.32e-2
        a(5,2)=3.91
        a(4,3)=2.21e-5
        a(5,3)=0.0
        a(5,4)=2.59
        t(2)=1.1121e-5
        t(3)=1.5702e-5
        t(4)=1.40100e-4
        t(5)=3.32657e-4
        weight(1)=5
        weight(2)=3
        weight(3)=1
        weight(4)=5
        weight(5)=1
c
      else if((ionum .eq. 18030) .or. (ionum .eq. 18031)
     1                           .or. (ionum .eq. 18032)) then
c
c  Ar IV ion parameters
c
        a(2,1)=2.23e-2
        a(3,1)=1.77e-3
        a(4,1)=8.62e-1
        a(5,1)=2.11
        a(3,2)=2.3e-5
        a(4,2)=6.03e-1
        a(5,2)=7.89e-1
        a(4,3)=0.119
        a(5,3)=0.598
        a(5,4)=4.94e-5
        t(2)=2.10904e-4
        t(3)=2.12193e-4
        t(4)=3.48555e-4
        t(5)=3.50326e-4
        weight(1)=4
        weight(2)=4
        weight(3)=6
        weight(4)=2
        weight(5)=4
c
      else if((ionum .eq. 18040) .or. (ionum .eq. 18041)
     1                           .or. (ionum .eq. 18042)) then
c
c  Ar V ion parameters
c
        a(2,1)=7.99e-3
        a(3,1)=1.24e-6
        a(4,1)=3.50e-5
        a(3,2)=2.72e-2
        a(4,2)=0.204
        a(5,2)=6.55
        a(4,3)=0.476
        a(5,3)=5.69e-2
        a(5,4)=3.29
        t(2)=7.639e-6
        t(3)=2.0292e-5
        t(4)=1.62994e-4
        t(5)=3.79125e-4
        weight(1)=1
        weight(2)=3
        weight(3)=5
        weight(4)=5
        weight(5)=1
c
      endif
c
      return
      end
c-----------------------------------------------------------------------
      subroutine parii(ionum,tem)
c
c     ________________________________________________________________
c     |                                                              |
c     |                           PAR II                             |
c     |                                                              |
c     |     Calculate temperature-dependent atomic parameters        |
c     |                                                              |
c     |     C:  collision strengths                                  |
c     |     CJ: coefficients for quadratic fit in T4(T/10000)        |
c     |                                                              |
c     ----------------------------------------------------------------
c
      real cj(3),ncrit
      integer weight(5)
      integer kbin,ttyout,filout
      common /atom/ c(5,5),a(5,5),ncrit(5),t(5),weight
      common /tty/ kbin,ttyout,filout
c
      t4=1.0e-4*tem
      cj(1)=0.0
      cj(2)=0.0
      cj(3)=0.0
c
c  For new ions, appropriate parameters must be included in IF block
c
      if((ionum .eq. 7000) .or. (ionum .eq. 7001)
     1                     .or. (ionum .eq. 7002)) then
c
c  N I ion parameters
c
        c(2,1)=-1.2124e-2+t4*(0.3616-t4*5.88e-2)
        c(3,1)=-8.386e-3+t4*(0.2419-t4*3.94e-2)
        c(4,1)=-2.601e-3+t4*(0.07003-t4*1.07e-2)
        c(5,1)=-5.074e-3+t4*(0.1395-t4*0.02126)
        c(3,2)=-4.166e-2+t4*(0.368-t4*0.05733)
        c(4,2)=1.227e-2+t4*(0.1046-t4*7.868e-3)
        c(5,2)=4.600e-2+t4*(0.244-t4*0.0240)
        c(4,3)=1.8601e-2+t4*(8.76e-2-t4*0.0092)
        c(5,3)=1.8265e-2+t4*(0.1406-t4*1.187e-2)
        c(5,4)=-3.265e-3+t4*(0.0704-t4*0.00387)
c
      else if((ionum .eq. 7010) .or. (ionum .eq. 7011)
     1                          .or. (ionum .eq. 7012)) then
c
c  N II ion parameters
c
        cj(1)=2.577+t4*(0.137-t4*0.03)
        cj(2)=0.353+t4*(-0.005807+t4*0.006003)
        c(2,1)=0.401
        c(3,1)=0.279
        c(4,1)=cj(1)/9
        c(5,1)=cj(2)/9
        c(3,2)=1.128
        c(4,2)=cj(1)*3/9
        c(5,2)=cj(2)*3/9
        c(4,3)=cj(1)*5/9
        c(5,3)=cj(2)*5/9
        c(5,4)=0.3993+t4*(0.0109+t4*0.001)
c
      else if((ionum .eq. 8000) .or. (ionum .eq. 8001)
     1                          .or. (ionum .eq. 8002)) then
c
c  O I ion parameters
c
        cj(1)=0.016656+t4*(0.3004-t4*0.02067)
        cj(2)=0.00201+t4*(0.03686-t4*0.002742)
        c(2,1)=-0.00214+t4*(0.09715+t4*0.003707)
        c(3,1)=-0.00109+t4*(0.0332-t4*0.00295)
        c(4,1)=cj(1)*5/9
        c(5,1)=cj(2)*5/9
        c(3,2)=-0.000128+t4*(0.0186+t4*0.00807)
        c(4,2)=cj(1)*3/9
        c(5,2)=cj(2)*3/9
        c(4,3)=cj(1)/9
        c(5,3)=cj(2)/9
        c(5,4)=0.0218+t4*(0.1076-t4*0.02234)
c
      else if((ionum .eq. 8010) .or. (ionum .eq. 8011)
     1                          .or. (ionum .eq. 8012)) then
c
c  O II ion parameters
c
        c(2,1)=0.7894+t4*(0.0098+t4*0.002282)
        c(3,1)=0.5269+t4*(5.204e-3+t4*1.923e-3)
        c(4,1)=0.26+t4*(1.001e-2-t4*2.915e-6)
        c(5,1)=0.1311+t4*(3.445e-3+t4*4.99e-4)
        c(3,2)=1.283+t4*(-0.1389+t4*0.0262)
        c(4,2)=0.7061+t4*(0.0235+t4*4.9e-4)
        c(5,2)=0.285+t4*0.01
        c(4,3)=0.3948+t4*(0.0123+t4*6.262e-4)
        c(5,3)=0.2644+t4*(0.0117-t4*9.16e-4)
        c(5,4)=0.2731+t4*(0.014-t4*2.919e-4)
c
      else if((ionum .eq. 8020) .or. (ionum .eq. 8021)
     1                          .or. (ionum .eq. 8022)) then
c
c  O III ion parameters
c
        cj(1)=1.835+t4*(0.3981-t4*0.06)
        cj(2)=0.2127+t4*(0.0767-0.013*t4)
        c(2,1)=0.4825+t4*(0.0806-t4*0.022)
        c(3,1)=0.2397+t4*(0.0381-t4*0.007)
        c(4,1)=cj(1)/9
        c(5,1)=cj(2)/9
        c(3,2)=1.1325+t4*(0.203-t4*0.05)
        c(4,2)=cj(1)/3
        c(5,2)=cj(2)/3
        c(4,3)=cj(1)*5/9
        c(5,3)=cj(2)*5/9
        c(5,4)=0.3763+t4*(0.3375-t4*0.105)
c
      else if((ionum .eq. 10020) .or. (ionum .eq. 10021)
     1                           .or. (ionum .eq. 10022)) then
c
c  Ne III ion parameters
c
        cj(1)=1.6028+t4*(0.07967-t4*0.03282)
        cj(2)=0.12956+t4*(0.05331-t4*0.01487)
        c(2,1)=1.0314+t4*(0.1544-t4*0.05638)
        c(3,1)=0.2910+t4*(0.02765-t4*0.012639)
        c(4,1)=cj(1)*5/9
        c(5,1)=cj(2)*5/9
        c(3,2)=0.3080+t4*(0.06120-t4*0.02096)
        c(4,2)=cj(1)/3
        c(5,2)=cj(2)/3
        c(4,3)=cj(1)/9
        c(5,3)=cj(2)/9
        c(5,4)=0.1739+t4*(0.05962-t4*0.00862)
c
      else if((ionum .eq. 10030) .or. (ionum .eq. 10031)
     1                           .or. (ionum .eq. 10032)) then
c
c  Ne IV ion parameters
c
        c(2,1)=0.8473+t4*(-0.005832-t4*0.002886)
        c(3,1)=0.5652+t4*(-0.00452-t4*0.001535)
        c(4,1)=0.1496+t4*(9.539e-3-t4*3.437e-3)
        c(5,1)=0.2976+t4*(0.0232-t4*0.0088)
        c(3,2)=1.362+t4*(0.0217-t4*0.0186)
        c(4,2)=0.3057+t4*(0.0859-t4*0.026)
        c(5,2)=0.8014+t4*(0.1364-t4*0.0415)
        c(4,3)=0.308+t4*(0.0391-t4*0.0119)
        c(5,3)=0.4291+t4*(0.1103-t4*0.0336)
        c(5,4)=0.2883+t4*(0.0662-t4*0.0127)
c
      else if((ionum .eq. 10040) .or. (ionum .eq. 10041)
     1                           .or. (ionum .eq. 10042)) then
c
c  Ne V ion parameters
c
        cj(1)=1.6175+t4*(0.171-t4*0.01)
        cj(2)=0.3315+t4*(-0.1142+t4*0.034)
        c(2,1)=0.244
        c(3,1)=0.122
        c(4,1)=cj(1)/9
        c(5,1)=cj(2)/9
        c(3,2)=0.578
        c(4,2)=cj(1)/3
        c(5,2)=cj(2)/3
        c(4,3)=cj(1)*5/9
        c(5,3)=cj(2)*5/9
        c(5,4)=1.27-0.02*t4
c
      else if((ionum .eq. 16010) .or. (ionum .eq. 16011)
     1                           .or. (ionum .eq. 16012)) then
c
c  S II ion parameters
c
        c(2,1)=3.06+t4*(-0.29+t4*0.02)
        c(3,1)=4.592+t4*(-0.437+t4*0.03)
        c(4,1)=0.78+t4*(-0.0242+t4*0.01)
        c(5,1)=1.33+t4*(0.264-t4*0.08)
        c(3,2)=8.79+t4*(-1.36+t4*0.16)
        c(4,2)=1.4575+t4*(0.111-t4*0.05)
        c(5,2)=3.282+t4*(0.223-t4*0.13)
        c(4,3)=2.502+t4*(0.1431-t4*0.09)
        c(5,3)=4.613+t4*(0.343-t4*0.17)
        c(5,4)=2.383+t4*(0.043-t4*0.05)
c
      else if((ionum .eq. 16020) .or. (ionum .eq. 16021)
     1                           .or. (ionum .eq. 16022)) then
c
c  S III ion parameters
c
        cj(1)=9.903+t4*(-2.017+t4*0.59)
        cj(2)=1.135+t4*0.0522
        c(2,1)=2.6725+t4*(0.01897-t4*0.13)
        c(3,1)=1.053+t4*(0.143-t4*0.05)
        c(4,1)=cj(1)/9
        c(5,1)=cj(2)/9
        c(3,2)=5.71+t4*(0.3181-t4*0.26)
        c(4,2)=cj(1)/3
        c(5,2)=cj(2)/3
        c(4,3)=cj(1)*5/9
        c(5,3)=cj(2)*5/9
        c(5,4)=0.82+t4*(1.424-t4*0.4)
c
      else if((ionum .eq. 17010) .or. (ionum .eq. 17011)
     1                           .or. (ionum .eq. 17012)) then
c
c  Cl II ion parameters
c
        c(2,1)=2.17
        c(3,1)=0.443
        c(4,1)=3.86*5/9
        c(5,1)=0.456*5/9
        c(3,2)=0.933
        c(4,2)=3.86/3
        c(5,2)=0.456/3
        c(4,3)=3.86/9
        c(5,3)=0.456/9
        c(5,4)=1.15
c
      else if((ionum .eq. 17020) .or. (ionum .eq. 17021)
     1                           .or. (ionum .eq. 17022)) then
c
c  Cl III ion parameters
c
        c(2,1)=1.1120+t4*(0.3837-t4*0.13750)
        c(3,1)=1.6660+t4*(0.5886-t4*0.21062)
        c(4,1)=0.3912+t4*(0.0085+t4*0.01505)
        c(5,1)=0.7810+t4*(0.0213+t4*0.02784)
        c(3,2)=4.0507+t4*(0.7741-t4*0.29264)
        c(4,2)=1.2051+t4*(0.6197-t4*0.18408)
        c(5,2)=1.8324+t4*(0.4803-t4*0.13531)
        c(4,3)=1.3373+t4*(0.2975-t4*0.08182)
        c(5,3)=3.2157+t4*(1.3672-t4*0.40935)
        c(5,4)=1.7478-t4*(0.0450-t4*0.05217)
c
      else if((ionum .eq. 17030) .or. (ionum .eq. 17031)
     1                           .or. (ionum .eq. 17032)) then
c
c  Cl IV ion parameters
c
        cj(1)=4.702+t4*(0.771-t4*0.01)
        cj(2)=1.712+t4*(0.791-t4*0.25)
        c(2,1)=0.475
        c(3,1)=0.4
        c(4,1)=cj(1)/9
        c(5,1)=cj(2)/9
        c(3,2)=1.5
        c(4,2)=cj(1)/3
        c(5,2)=cj(2)/3
        c(4,3)=cj(1)*5/9
        c(5,3)=cj(2)*5/9
        c(5,4)=0.3388+t4*(1.3214-t4*0.265)
c
      else if((ionum .eq. 18020) .or. (ionum .eq. 18021)
     1                           .or. (ionum .eq. 18022)) then
c
c  Ar III ion parameters
c
        c(2,1)=2.24
        c(3,1)=0.531
        c(4,1)=4.74*5/9
        c(5,1)=0.68*5/9
        c(3,2)=1.18
        c(4,2)=4.74/3
        c(5,2)=0.68/3
        c(4,3)=4.74/9
        c(5,3)=0.68/9
        c(5,4)=0.823
c
      else if((ionum .eq. 18030) .or. (ionum .eq. 18031)
     1                           .or. (ionum .eq. 18032)) then
c
c  Ar IV ion parameters
c
        c(2,1)=2.2661-t4*(1.2805-t4*0.32167)
        c(3,1)=3.3993-t4*(1.9217-t4*0.48293)
        c(4,1)=0.1637-t4*(0.0351-t4*0.01790)
        c(5,1)=0.3356-t4*(0.0817-t4*0.03930)
        c(3,2)=6.4696-t4*(0.3631-t4*0.04479)
        c(4,2)=1.4523+t4*(0.4098-t4*0.15965)
        c(5,2)=2.3424+t4*(0.2396-t4*0.11250)
        c(4,3)=1.7193+t4*(0.1391-t4*0.07235)
        c(5,3)=3.9465+t4*(0.7872-t4*0.30578)
        c(5,4)=2.1475+t4*(0.1047+t4*0.09440)
c
      else if((ionum .eq. 18040) .or. (ionum .eq. 18041)
     1                           .or. (ionum .eq. 18042)) then
c
c  Ar V ion parameters
c
        cj(1)=5.2075+t4*(-1.985+t4*0.55)
        cj(2)=1.133+t4*(0.127-t4*0.09)
        c(2,1)=0.257
        c(3,1)=0.32
        c(4,1)=cj(1)/9
        c(5,1)=cj(2)/9
        c(3,2)=1.04
        c(4,2)=cj(1)/3
        c(5,2)=cj(2)/3
        c(4,3)=cj(1)*5/9
        c(5,3)=cj(2)*5/9
        c(5,4)=1.27+t4*(-0.02+t4*2.4e-5)
c
      endif
      return
      end
c------------------------------------------------------------------------
      subroutine solve(den,tem,pop,jl)
c
c     _________________________________________________________________
c     |                                                               |
c     |                            SOLVE                              |
c     |                                                               |
c     |     TS:       square root of temperature                      |
c     |     POP1:     storage array for level populations             |
c     |     POP:      level population array                          |
c     |     A:        radiative transition probabilities              |
c     |     C:        collision strengths                             |
c     |     Q:        collision transition probabilities              |
c     |     E:        relative energy separation                      |
c     |     KT,HC,CQ: "constants"                                     |
c     |                                                               |
c     -----------------------------------------------------------------
c
      real e(5,5), q(5,5), a1(5,5), a2(5,5), jl(5,5)
      real pop(5), b1(5), b2(5), dtl(5), pop1(5), ncrit, kt, hc
      integer weight(5)
      integer kbin,ttyout,filout
c
      common /atom/c(5,5),a(5,5),ncrit(5),t(5),weight
      common /tty/kbin,ttyout,filout
c
      kt=(1.3807e-16)*tem
      hc=1.9865e-08
      cq=8.629e-06
      ts=sqrt(tem)
c
c  Transition energy differences
c
      do 100 i=1,5
         do 110 j=i+1,5
            e(j,i)=hc*(t(j)-t(i))
  110    continue
  100 continue
c
c  Collision transition probabilities (if i=j, Q=0)
c
      do 200 i=1,5
         do 210 j=1,5
            if(j .gt. i) then
              q(i,j)=cq*c(j,i)*exp(-e(j,i)/kt)/(weight(i)*ts)
            else if(j .eq. i) then
              q(i,j)=0.0
            else
              q(i,j)=cq*c(i,j)/(weight(i)*ts)
            endif
  210    continue
  200 continue
c
c  Critical density calculation
c
      do 300 i=2,5
         asum=0.0
         qsum=0.0
         do 400 j=1,i-1
            asum=asum+a(i,j)
            qsum=qsum+q(i,j)
  400    continue
         ncrit(i)=asum/qsum
  300 continue
c
c  Matrix manipulation
c
      do 500 i=1,5
         do 600 j=1,5
            if(j .eq. i) then
              a1(i,j)=0.0
            else if(j .gt. i) then
              a1(i,j)=den*q(i,j)
            else
              a1(i,j)=den*q(i,j)+a(i,j)
            endif
  600    continue
  500 continue
c
      do 700 i=1,5
         b1(i)=0.0
         do 800 j=1,5 
            b1(i)=b1(i)+a1(i,j)
  800    continue
  700 continue
c
      do 900 i=1,5
         b2(i)=-a1(1,i)
         a1(i,i)=-b1(i)
  900 continue
c
      do 1000 i=1,5
         do 1100 j=1,5
            do 1200 k=1,5
               a2(j,k)=a1(j,k)
 1200       continue
 1100    continue
         do 1300 j=1,5
            a2(i,j)=b2(j)
 1300    continue
         dtl(i)=1
         do 1400 j=2,4
            t1=a2(j,j)
            if(t1 .eq. 0.0) then
              write(filout,*)'Divide by 0 during T1 normalization
     1 in SOLVE'
              stop
            endif
            dtl(i)=dtl(i)*t1
            do 1500 k=j,5
               a2(j,k)=a2(j,k)/t1
 1500       continue
            do 1600 k=j+1,5
               t2=a2(k,j)
               do 1700 l=j,5 
                  a2(k,l)=a2(k,l)-t2*a2(j,l)
 1700          continue
 1600       continue
 1400    continue
c
         dtl(i)=dtl(i)*a2(5,5)
 1000 continue
c
      pop1(1)=1
      do 1800 i=2,5
         pop1(i)=dtl(i)/dtl(1)
 1800 continue
c
      t1=0
      do 1900 i=1,5
         t1=t1+pop1(i)
 1900 continue
c
c  Finally, level populations, emissivities
c
      do 2000 i=1,5
         pop(i)=pop1(i)/t1
         do 2100 j=1,5
            jl(i,j)=pop(i)*a(i,j)*e(i,j)
            if(jl(i,j) .lt. 1.0e-34) jl(i,j)=1.0e-34
 2100    continue
 2000 continue
c
      return
      end
c------------------------------------------------------------------------
      subroutine temden(ionum,den,tem,robs,pop,jl)
c
c     _________________________________________________________________
c     |                                                               |
c     |                           TEMDEN                              |
c     |                                                               |
c     |     ITER:     iterative index used to find TEM or DEN         |
c     |     RCAL:     'characteristic' ratio calculated in SOLVE      |
c     |     X1,X2:    initial (absolute) limits on TEM or DEN         |
c     |     XLO,XHI:  current range on iterated TEM or DEN            |
c     |     RLO,RHI:  ratios calculated at XLO, XHI                   |
c     |     X:        current 'best guess' for TEM or DEN             |
c     |     ERR:      percentage tolerance in X                       |
c     |                                                               |
c     -----------------------------------------------------------------
c
      real a(5,5), c(5,5), jl(5,5), cj(3), pop(5), t(5)
      integer kbin,ttyout,filout

      common /tty/ kbin,ttyout,filout

      err = 0.001
c
c  Temperature limits
c
      if(mod(ionum,10) .eq. 2) then
        x1=1.0e3
        x2=3.6e4
        write(filout,20)
c
c  Density limits
c
      elseif(mod(ionum,10) .eq. 1) then
        x1=1.0
        x2=1.0e8
        write(filout,10)
      endif
c
c  Permit up to 100 iterations
c
      do 100 iter=1,100
         if(iter .eq. 1) then
           xlo=x1
           xhi=x2
           x=xlo
         else if(iter .eq. 2) then
           xhi=x2
           x=xhi
         else
           x=sqrt(xlo*xhi)
         endif
c
         if(mod(ionum,10) .eq. 2) then
           tem=x
         else if(mod(ionum,10) .eq. 1) then
           den=x
         endif
c
         call parii(ionum,tem)
         call solve(den,tem,pop,jl)
c
c  For new ions, enter appropriate ratio in IF block below
c
         if(ionum .eq. 7001) then
           rcal=jl(2,1)/jl(3,1)
         else if(ionum .eq. 7012) then
           rcal=(jl(4,3)+jl(4,2))/jl(5,4)
         else if(ionum .eq. 8011) then
           rcal=jl(3,1)/jl(2,1)
         else if(ionum .eq. 8012) then
           rcal=(jl(2,1)+jl(3,1))/(jl(5,3)+jl(5,2)+jl(4,3)+jl(4,2))
         else if(ionum .eq. 8022) then
           rcal=(jl(4,3)+jl(4,2))/jl(5,4)
         else if(ionum .eq. 10022) then
           rcal=(jl(4,1)+jl(4,2))/jl(5,4)
         else if(ionum .eq. 10031) then
           rcal=jl(3,1)/jl(2,1)
         else if(ionum .eq. 10042) then
           rcal=(jl(4,3)+jl(4,2))/jl(5,4)
         else if(ionum .eq. 16011) then
           rcal=jl(3,1)/jl(2,1)
         else if(ionum .eq. 16012) then
           rcal=(jl(2,1)+jl(3,1))/(jl(4,1)+jl(5,1))
         else if(ionum .eq. 16022) then
           rcal=(jl(4,3)+jl(4,2))/jl(5,4)
         else if(ionum .eq. 17021) then
           rcal=jl(3,1)/jl(2,1)
         else if(ionum .eq. 17032) then
           rcal=(jl(4,2)+jl(4,3))/jl(5,4)
         else if(ionum .eq. 18022) then
           rcal=(jl(4,1)+jl(4,2))/jl(5,4)
         else if(ionum .eq. 18031) then
           rcal=jl(3,1)/jl(2,1)
         else if(ionum .eq. 18032) then
           rcal=(jl(2,1)+jl(3,1))/(jl(4,1)+jl(5,1))
         else if(ionum .eq. 18042) then
           rcal=(jl(4,2)+jl(4,3))/jl(5,4)
         endif
c
c  Output intermediate iterations
c
         if(mod(ionum,10) .eq. 1) then
           write(filout,30) iter,den,xlo,tem,xhi,(rcal-robs)/robs
         else if(mod(ionum,10) .eq. 2) then
           write(filout,30) iter,tem,xlo,den,xhi,(rcal-robs)/robs
         endif
c
c  Test for convergence
c
         if(iter .eq. 1) then
           rlo=rcal
         else if((iter .eq. 2) .and. ((rcal-robs)/(rlo-robs) .lt. 
     . 0.0)) then
           rhi=rcal
         else if((iter .eq. 2) .and. ((rcal-robs)/(rlo-robs) .gt.
     . 0.0)) then
           write(filout,40) amin1(rlo,rcal),amax1(rlo,rcal)
           write(filout,50) rcal,robs,rlo
           robs=0.0
           return
         else if(abs(0.5*(xhi-xlo)/(xhi+xlo)) .lt. err) then
           return
         else if((rcal-robs)/(rhi-robs) .lt. 0.0) then
           rlo=rcal
           xlo=x
         else if((rcal-robs)/(rhi-robs) .gt. 0.0) then
           rhi=rcal
           xhi=x
         endif
  100 continue
c
      write(ttyout,*) 'Did not converge after 100 iterations; stopping'
   10 format(1x,'Iter#',6x,'Density',8x,'XLO',2x,'Temperature',7x,
     1'XHI',3x,' (RCAL-ROBS)/ROBS')
   20 format(1x,'Iter#',2x,'Temperature',8x,'XLO',6x,'Density',7x,
     1'XHI',3x,' (RCAL-ROBS)/ROBS')
   30 format(3x,i3,2(3x,1pe10.3,1x,1pe10.3),9x,1pe9.2)
   40 format(' (RCAL-ROBS)/(RLO-ROBS) is > 0, so the entered RATIO',/,
     1    ' predicts unreasonable conditions. Please try again.',/,
     2    ' The allowable RATIO range is:',1pe9.2,'--',1pe9.2)
   50 format(' RCAL = ',f8.3,' ROBS = ',f8.3,' RLO = ',f8.3)
c
      return
      end

