**==aa0001.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      IMPLICIT NONE
c**==AA0001.spg  processed by SPAG 3.09I  at 14:01 on  1 Mar 1993
C          Program to compute lines and cooling from
c             a hot diffuse gas.
C
c      (units 2,7,8,9,10,11 reserved for I/O)
c      First a few common blocks.
 
      INCLUDE 'rayspec.inc'
 
      INTEGER NBIn , NJ(12)
      REAL ABUnd , ABUnj(12) , BINmin , BINsyz
      COMMON /PARAMS/ NJ , ABUnj , ABUnd , BINmin , BINsyz , NBIn
 
      INTEGER NBLn
      REAL BLN(MXBINS) , BLNmin , BLNsyz
      COMMON /BLN   / BLN , BLNmin , BLNsyz , NBLn
 
      REAL CONce(30) , DNE , GNDrec(30) , HENeut , HEPlus , PCOol , 
     &     PM(4) , POT , POU , POWer(220) , RE , RHY , TU
      COMMON /RESULT/ CONce , GNDrec , POWer , RHY , HENeut , HEPlus , 
     &                DNE , PCOol , POU , POT , RE , TU , PM
 
      REAL RF(500) , TAU(150) , TAUhe(150) , TPLus(150)
      COMMON /PT    / RF , TAU , TAUhe , TPLus
 
      REAL CNC(12,30) , PTOt(12,220) , ABIn(MXBINS) , BIN(MXBINS)
      COMMON /COM   / CNC , PTOt , ABIn , BIN
 
c Local variables
 
      REAL be , CION , cionh , coolrate , crec , dene , t
      INTEGER i , icx , idens , ier , iphot , iprint , jprint , k , 
     &        ncmd , ncmdwl , num
      CHARACTER*1 irep , iwl , iplot , char
      CHARACTER*70 cmd(10) , cmdwl(10)
      CHARACTER cmdtmp*70
      LOGICAL*4 flag , backflag , repflag
 
      DATA ABUnj/10.93 , 8.52 , 7.96 , 8.82 , 7.92 , 7.42 , 7.52 , 7.2 , 
     &     6.9 , 6.3 , 7.5 , 6.3/
      DATA RF/500*1.0/
 
      WRITE (*,'(a)') 
     &               ' This program computes equilibrium optically thin'
     &               , ' X-ray and EUV emission spectra.' , 
     &               ' More documentation can be found in ray.doc'
c
      repflag = .TRUE.
      irep = 'n'
      NBIn = 100
      BINmin = 100.
      BINsyz = 10.
      iwl = 'n'
      num = 12
      idens = 1
      icx = 1
      iprint = 0
      jprint = 0
      t = 1E7
      dene = 1.
      iplot = 'y'
c
      DO WHILE ( repflag )
c
 50      PRINT * , ' No. of bins for emission spectrum. Must be <= ' , 
     &         MXBINS
         PRINT * , ' 0 will just output lines. (' , NBIn , ')'
         READ (*,*)  , NBIn
 100     PRINT * , 'Energy in eV of lowest energy and binsize (' , 
     &         BINmin , ',' , BINsyz , ')'
         READ (*,*,END=50)  , BINmin , BINsyz
c
 150     NBLn = 0
         CALL QUESTCHAR('Do you want wavelength bins',flag,backflag,iwl)
         IF ( backflag ) GOTO 100
         IF ( flag ) CALL WAVEBINS(NBIn,BINmin,BINsyz,NBLn,BLNmin,
     &                             BLNsyz)
c
 200     char = 'y'
         CALL QUESTCHAR('Do you want to skip the next few questions',
     &                  flag,backflag,char)
         IF ( backflag ) GOTO 150
         IF ( .NOT.flag ) THEN
 220        PRINT * , ' No. of elements other than Hydrogen? (' , num , 
     &            ')'
            READ (*,*,END=200)  , num
            iphot = 0
 240        PRINT * , 
     &            ' If you want photoionisation then put Iphot .ne. 0 ('
     &            , iphot , ')'
            READ (*,*,END=220)  , iphot
            DO WHILE ( iphot.NE.0 )
               PRINT * , 
     &' In which case you had better look at the subroutine    radrd bec
     &ause it reads in some data from unit 15 which does     not appear 
     &to exist.'
               READ (*,*,END=220)  , iphot
            ENDDO
c
 260        PRINT * , 
     &' Input Idens .eq. 0 to switch off ionisation from  metastable lev
     &els and density dependent diel. recomb. (' , idens , ')'
            READ (*,*,END=240)  , idens
 280        PRINT * , 
     &            ' Input Icx .eq. 0 to switch off charge transfer (' , 
     &            icx , ')'
            READ (*,*,END=260)  , icx
         ENDIF
 300     DO WHILE ( .TRUE. )
c
c               Level of Output Required.
            WRITE (*,'(a/a,I1,a)') 
     &         ' Level of output for ionisation balance,lines etc 0-3 ?'
     &         , ' -1 for a list. (' , iprint , ')'
            READ (*,*,END=280)  , iprint
            IF ( iprint.LT.0 ) THEN
               WRITE (*,'(a)') ' 0 outputs no separate line information'
     &                         , ' 1 outputs ionisation fractions' , 
     &                         ' 2 outputs above + major lines' , 
     &           ' 3 outputs above + ionisation and recombination rates'
               iprint = 0
               GOTO 350
            ENDIF
            IF ( iprint.GT.0 ) OPEN (UNIT=8,FILE='lines.out',STATUS=
     &                               'new')
            GOTO 400
 350     ENDDO
 400     DO WHILE ( .TRUE. )
            WRITE (*,'(a/a,I1,a)') 
     &                           ' Level of output for continuum 0-2 ? '
     &                           , ' -1 for a list. (' , jprint , ')'
            READ (*,*,END=300)  , jprint
            IF ( jprint.LT.0 ) THEN
               WRITE (*,'(a)') 
     &                    ' 0 outputs no separate continuum information'
     &                    , 
     &                    ' 1 outputs recomb., brems. and 2phot. rates'
     &                    , ' 2 outputs above + major edges'
               jprint = 0
               GOTO 450
            ENDIF
            IF ( jprint.GT.0 ) CALL OPENWR(7,'continuum.out','new',' ',
     &           'l',0,0,ier)
            IF ( jprint.GT.1 ) OPEN (UNIT=11,FILE='edges.out',STATUS=
     &                               'new')
C
            CALL OPENWR(2,'atomic.data','old',' ',' ',0,1,ier)
            CALL ATREAD(num)
            CLOSE (2)
C
            IF ( iphot.NE.0 ) CALL RADRD(1)
            GOTO 500
 450     ENDDO
 500     DO WHILE ( .TRUE. )
C
c              Modify Abundances if Required.
            char = 'n'
            CALL QUESTCHAR('Do you wish to modify abundances',flag,
     &                     backflag,char)
            IF ( .NOT.(backflag) ) THEN
               IF ( flag ) CALL MODABUND(ABUnj)
c
               t = ALOG10(t)
               GOTO 550
            ENDIF
         ENDDO
 550     PRINT * , ' Log temperature. Allowable values 4 to 8(' , t , 
     &         ')'
         READ (*,*,END=500)  , t
         t = 10.**t
C COMPUTE RATIO OF IONIZED TO NEUTRAL HYDROGEN
         cionh = CION(1,1,13.6,t)
         be = 157890./t
         crec = 2.06E-11*(.4288+.5*ALOG(be)+.469*be**(-.33333))/SQRT(t)
         RHY = cionh/crec
         dene = ALOG10(dene)
         PRINT * , ' Log density. Should be .le. 12(' , dene , ')'
         READ (*,*,END=550) dene
         dene = 10.**dene
C
C
         DO 600 i = 1 , NBIn
            BIN(i) = 0.
 600     CONTINUE
         CALL FELINE(t,dene,num,iprint,jprint,0,iphot,idens,icx)
c
         CLOSE (7)
         CLOSE (8)
         CLOSE (11)
c
c              Print Out Total Spectrum Versus Energy.
         IF ( NBIn.GT.0 ) THEN
            OPEN (UNIT=9,FILE='energy.out',STATUS='new')
            WRITE (9,99001) t , BINmin , BINsyz , ABUnj
            coolrate = 0.
            DO 620 i = 1 , NBIn
               coolrate = coolrate + BIN(i)
 620        CONTINUE
            WRITE (9,99002) coolrate
            DO 640 k = 1 , NBIn
               WRITE (9,99003) (BINmin+(k-0.5)*BINsyz) , BIN(k)
 640        CONTINUE
            CLOSE (9)
c               Print Out Total Spectrum Versus Wavelength.
            IF ( NBLn.GT.0 ) THEN
               OPEN (UNIT=10,FILE='wlength.out',STATUS='new')
               WRITE (10,99001) t , BLNmin , BLNsyz , ABUnj
               WRITE (10,99002) coolrate
               DO 650 k = 1 , NBLn
                  WRITE (10,99003) (BLNmin+(k-0.5)*BLNsyz) , BLN(k)
 650           CONTINUE
               CLOSE (10)
            ENDIF
         ENDIF
c
         IF ( NBIn.GT.0 ) THEN
            CALL QUESTCHAR('Do you want an energy plot',flag,backflag,
     &                     iplot)
            IF ( flag ) THEN
c    Plot out the emissivity versus energy
 
               cmdtmp = 'xaxis lin '
               WRITE (cmdtmp(11:21),'(1pe11.3)') BINmin + 0.5*BINsyz
               WRITE (cmdtmp(22:33),'(1x,1pe11.3)') BINsyz
               cmd(1) = cmdtmp
 
               cmd(2) = 'la x energy (eV)'
               cmd(3) = 
     &'la y emissivity (10\u-23\dn\de\un\dH\uerg s\u-1\dcm\u3\dbin\u-1\d
     &)'
               cmd(4) = 'la file emissivity per energy bin'
               cmd(5) = 'time off'
               cmd(6) = 'log y'
               cmd(7) = 'LI ON'
               ncmd = 7
               CALL PLT(BIN,0,MXBINS,NBIn,1,cmd,ncmd,ier)
            ENDIF
         ENDIF
c
         IF ( NBLn.GT.0 ) THEN
            CALL QUESTCHAR('Do you want a wavelength plot',flag,
     &                     backflag,iplot)
            IF ( flag ) THEN
c   Plot out the emissivity versus wavelength
 
               cmdtmp = 'xaxis lin '
               WRITE (cmdtmp(11:21),'(1pe11.3)') BLNmin + 0.5*BLNsyz
               WRITE (cmdtmp(22:33),'(1x,1pe11.3)') BLNsyz
               cmdwl(1) = cmdtmp
 
               cmdwl(2) = 'la x wavelength (A)'
               cmdwl(3) = 
     &'la y emissivity (10\u-23\dn\de\un\dH\uerg s\u-1\dcm\u3\dbin\u-1\d
     &)'
               cmdwl(4) = 'la file emissivity per wavelength bin'
               cmdwl(5) = 'time off'
               cmdwl(6) = 'log y'
               cmdwl(7) = 'LI ON'
               ncmdwl = 7
               CALL PLT(BLN,0,MXBINS,NBLn,1,cmdwl,ncmdwl,ier)
            ENDIF
         ENDIF
         DO WHILE ( .TRUE. )
c
            CALL QUESTCHAR('Do you want another run',repflag,backflag,
     &                     irep)
            IF ( .NOT.(backflag) ) GOTO 700
         ENDDO
c
 700  ENDDO
      STOP
99001 FORMAT (' Temperature=',1pe10.2,' Binmin=',0pf5.1,' Binsize=',
     &        F5.1,/,' Abund',12F6.2)
99002 FORMAT (' Total cooling rate in this range is ',1P1e12.4)
99003 FORMAT (1P2e12.4)
      END
