C
C PURPOSE: GENERAL ROUTINE FOR COMPUTING COLLISIONAL RATES
C     IN THE SYSNTHESIS CODE FOR M1 LINES, THIS SHOULD BE USED 
C     ONLY FOR DETERMINING IONIZATION BALANCE
C
C INPUTS:
C
C OUTPUTS:
C
C COMMON:
C     CATOM CATMOS CATMO2 CCONST CINPUT CLU
C COMMENTS: October 6, 1999, P. JUDGE
C 
C
C*********************************************************************
C
      SUBROUTINE GENCOL
C
C  GENERAL ROUTINE FOR COMPUTING COLLISIONAL RATES (VERSION MAR '89)
C  FORMAT IS EITHER GIVEN IN TERMS OF A COLLISION STRENGTH OMEGA
C         OR                          A CE(TE) RATE
C
C  OMEGAS ARE USED FOR CHARGED IONS IN GENERAL, 
C  SINCE THEY REMAIN ROUGHLY CONSTANT WITH TEMPERATURE.  
C  RATE IS PROP. TO 1/SQRT(TE) * EXP(DELTA-E)
C
C  CE(TE) VALUES ARE USED LARGELY FOR NEUTRALS FOR SAME REASON.
C  RATE IS PROP. TO SQRT(TE) * EXP(DELTA-E)
C
C  CI(TE) VALUES ARE USED FOR IONIZATION RATES: THESE HAVE A
C         SQRT(TE) * EXP(DELTA-E) DEPENDENCE AS FOR CE(TE)
C
C  ORIGINALLY CODED BY P.G. JUDGE,  APRIL 2ND, 1987
C  MODIFIED FOR FORTRAN-77 STANDARD BY P.G. JUDGE, MAY 21ST 1987
C:
C: GENCOL 87-05-21  NEW ROUTINE: (PHILIP JUDGE)
C:        GENERAL COLLISONAL ROUTINE READING TABLES
C:
C:        89-03-24  MODFICATIONS: (PHILIP JUDGE)
C:        NEW INPUT 'CP' AND 'CH' ADDED FOR INELASTIC COLLISIONS
C:        WITH NEUTRAL AND IONIZED HYDROGEN.
C:        'SEMI ' CHANGED SO THAT THE ABSORPTION OSCILLATOR STRENGTH 
C:        IS REQUIRED FOR INPUT.
C:
C:        89-06-07  MODIFICATIONS: (MATS CARLSSON)
C:        DO LOOP VARIABLE IT CHANGED TO JT IN ORDER NOT TO CONFLICT
C:        WITH THE MAIN ITERATION VARIABLE IT IN ITER
C:
C:        90-04-10  MODIFICATIONS: (PHILIP JUDGE)
C:        LOW TEMPERATURE DIELECTRONIC RECOMBINATION OPTION ADDED
C:        FOLLOWING NUSSBAUMER AND STOREY (A+A  )
C:        NOTE THAT ONLY THE RATE DOWNWARD IS ADDED TO THE FIXED RATE
C:
C:        90-08-27  MODIFICATIONS: (MATS CARLSSON)
C:        DO LOOP 1 REPLACED BY A GOTO TO REMOVE RESTRICTION ON
C:        NUMBER OF LINES WITH COLLISIONAL DATA
C:
C:        92-08-10  MODIFICATIONS (MATS CARLSSON)
C:        INTEGER ARRAYS WITH ELEMENT NAME CHANGED TO CHARACTER ARRAYS
C:
C:        92-09-07  MODIFICATIONS (MATS CARLSSON)
C:        ADDED OPTIONS RECOMB AND ALPHA
C:
C:        94-01-24  MODIFICATIONS (PHILIP JUDGE)
C:        AR85 OPTIONS ADDED.  ARNAUD AND ROTHENFLUG 1985
C:        94-02-22  MODIFICATIONS (PHILIP JUDGE)
C:        AR85-CEA FORMALISM ADDED
C:        94-03-08  MODIFICATIONS (PHILIP JUDGE)
C:        AR85-CH, AR85-CH+, AR85-CHE, AR85-CHE+, FORMALISM ADDED
C
      INCLUDE 'PREC'
      PARAMETER (MTGRD=50,MSHELL=5)
      INCLUDE 'PARAM'
      INCLUDE 'CATOM'
      INCLUDE 'CATMOS'
      INCLUDE 'CATMO2'
      INCLUDE 'CCONST'
      INCLUDE 'CINPUT'
      INCLUDE 'CLU'
C
C  LOCAL VARIABLES: UP TO MTGRD POINTS IN TEMPERATURE GRID ALLOWED
C  THIS CAN BE INCREASED
C
      DIMENSION TGRID(MTGRD),CGRID(MTGRD)
      DIMENSION CDI(5,MSHELL)
      CHARACTER *20 KEY
      CNST=0.0
C
C  OUTPUT FOR IWATOM .GT. 1
C
      IF(IWATOM .GT. 1)WRITE(LOUT,800)
  800   FORMAT('1 COLLISION RATES FROM SUBROUTINE  G E N C O L'//
     *  1X,'PARAMETERS ARE AS FOLLOWS:'//
     *  1X,'TEMP -   TEMPERATURE GRID GIVEN IN FILE ATOM'/
     *  1X,'         (SPLINE INTEPOLATIONS ARE MADE)'//
     *  1X,'CE   -   NEUTRAL B-B COLLISION RATE WHERE'/
     *  1X,'         R(U-L) = CE * NE * G(L) * SQRT(TEMP) / G(U)'//
     *  1X,'OHM  -   ION B-B COLLISION RATE WHERE'/
     *  1X,'         R(U-L) = OHM * NE * 8.63E-6 / SQRT(TEMP) / G(U)'//
     *  1X,'SEMI -   B-B COLLISIONAL RATE USING THE SEMI-EMPIRICAL')
      IF(IWATOM .GT. 1) WRITE(LOUT,801)
  801   FORMAT(1X,'         FORMULAE OF VAN REGEMORTER AND SHEVELKO'//
     *  1X,'CP   -   B-B COLLISIONAL RATE WITH PROTONS WHERE'/
     *  1X,'         R(U-L) = CP * NP '//
     *  1X,'CH   -   B-B COLLISIONAL RATE WITH NEUTRAL HYDROGEN WHERE'/
     *  1X,'         R(U-L) = CH * N(H I) '//
     *  1X,'CI   -   B-F IONIZATION RATE WHERE'/
     *  1X,'         R(L-U) = CI * NE * EXP(-DELTAE/KT) * SQRT(T)'//
     *  1X,'CALP -   COLLISIONS WITH ALPHA PARTICLES'/
     *  1X,'         R(I-J) = CALP * N(NK)'///
     *  1X,'***********************************************'/
     *  1X,'*** ALL THE ABOVE HAVE REVERSE RATES IN LTE ***'/
     *  1X,'*** THE FOLLOWING HAVE ZERO REVERSE RATES   ***'/
     *  1X,'*** UNLESS SPECIFIED EXPLICITLY             ***'/
     *  1X,'***********************************************'//
     *  1X,'CH0  -   CHARGE TRANSFER RATE WITH NEUTRAL HYDROGEN'/
     *  1X,'         R(I-J) = CH0 * NH(LEVEL1)'//
     *  1X,'CH+  -   CHARGE TRANSFER RATE WITH IONIZED HYDROGEN'/
     *  1X,'         R(I-J) = CH+ * N(H II)'///
     *  1X,'RECO -   RECOMBINATION RATE'/
     *  1X,'         R(J-I) = RECO * NE'///
     *  1X,'CONS -   CONSTANT RATE'/
     *  1X,'         R(J-I) = CONST '///
     *  1X,'KEY ','FROM  TO  PARAMETERS... '/)
C
C  ASSUME INITIALLY THAT DATA ARE GIVEN INDEPENDENT OF TEMPERATURE:
C  I.E. THAT NTEMP IS ORIGINALLY SET = 1
C
      NTEMP=1
C
C  READ THE KEYWORD 'KEY' AND ASSOCIATED PARAMETERS
C
   50 CONTINUE
        READ(LDUMS,100,END=301)KEY
  100   FORMAT(A)
        CALL LJUST(KEY)
        IF(KEY(1:1) .EQ. ' ') GOTO 50
        IF (KEY(1:3) .EQ. 'END')THEN
          GOTO 300
        ELSE IF (KEY(1:4) .EQ. 'TEMP')THEN
          READ(LDUMS,*)NTEMP,(TGRID(JT),JT=1,MIN(NTEMP,MTGRD))
          IF(NTEMP .GT. MTGRD)THEN
            CALL STOP(' GENCOL: WORK ARRAYS (TGRID) TOO SMALL')
          ENDIF
          IF (IWATOM .GT. 1)THEN
            WRITE(LOUT,803)KEY(1:4),(TGRID(JT),JT=1,NTEMP)
  803       FORMAT(1X,A,8X,10(1P,1X,E10.3))
          ENDIF
          GOTO 50
        ELSE IF (KEY(1:4) .EQ. 'SEMI') THEN
          READ(LDUMS, *) IL,IH, FAB
          IF(IWATOM .GT. 1)WRITE(LOUT,8004)KEY(1:4),IL,IH
 8004     FORMAT(1X,A,2(2X,I2))
        ELSE IF (KEY(1:4) .EQ. 'LTDR') THEN
          READ(LDUMS, *) IL,IH, ALTDR,BLTDR,CLTDR,DLTDR,FLTDR
          IF(IWATOM .GT. 1)WRITE(LOUT,8005)KEY(1:4),IL,IH,ALTDR,
     *      BLTDR,CLTDR,DLTDR,FLTDR
 8005     FORMAT(1X,A,2(2X,I2),5(2X,F11.4))
C
C  ADDITION P. JUDGE 24-JAN-1994: BEGIN
C
        ELSE IF (KEY(1:7) .EQ. 'AR85-RR') THEN
          READ(LDUMS, *) IL,IH,ARADAR,ETAAR
        ELSE IF (KEY(1:8) .EQ. 'AR85-CDI') THEN
          READ(LDUMS, *) IL,IH,NCDI
          IF(NCDI .GT. MSHELL) CALL STOP('GENCOL: NCDI .GT. MSHELL')
          READ(LDUMS,*) ((CDI(J,I),J=1,5),I=1,NCDI)
C
C  ADDITION P. JUDGE 24-JAN-1994: END
C
C:
C: 22-FEB-1994 P.G.JUDGE MODIFICATIONS START:
C:
        ELSE IF (KEY(1:8) .EQ. 'AR85-CEA') THEN
          READ(LDUMS, *) IL,IH,CEAFAK
          ILO=MIN( IL, IH )
          IHI=MAX( IL, IH )
          CALL AR85CEA(ILO,IHI,CEA)
            CEA=CEA*CEAFAK
C:
C: 22-FEB-1994 P.G.JUDGE MODIFICATIONS END:
C:
C:
C: 08-MAR-1994 P.G.JUDGE MODIFICATIONS START:
C:  ALL FOUR CHARGE TRANSFER KEYWORDS HERE
        ELSE IF (KEY(1:7) .EQ. 'AR85-CH') THEN
          READ(LDUMS, *) IL,IH
          READ(LDUMS,*) AR85T1,AR85T2,AR85A,AR85B,AR85C,AR85D
C:
C: 08-MAR-1994 P.G.JUDGE MODIFICATIONS END:
C:
        ELSE IF (KEY(1:7) .EQ. 'SHULL82') THEN
          READ(LDUMS, *) IL,IH,ACOLSH,TCOLSH,ARADSH,XRADSH,ADISH,
     *      BDISH,T0SH,T1SH
          IF(IWATOM .GT. 1)WRITE(LOUT,8006)KEY(1:6),IL,IH,ACOLSH,
     *    TCOLSH,ARADSH,XRADSH,ADISH,BDISH,T0SH,T1SH
 8006     FORMAT(1X,A,2(2X,I2),8(2X,1P,E11.4))
        ELSE IF (KEY(1:6) .EQ. 'CORONA') THEN
          READ(LDUMS, *) IL,IH
          IF(IWATOM .GT. 1)WRITE(LOUT,8004)KEY(1:6),IL,IH
        ELSE IF (KEY(1:8) .EQ. 'CONSTANT') THEN
          READ(LDUMS, *) IL,IH, CNST
          IF(IWATOM .GT. 1)WRITE(LOUT,8014)KEY(1:4),IL,IH, CT
 8014     FORMAT(1X,A,2(2X,I2),2X,1P,E11.4)
        ELSE IF (KEY(1:3) .EQ. 'OHM'
     *  .OR.     KEY(1:2) .EQ. 'CE'
     *  .OR.     KEY(1:2) .EQ. 'CP'
     *  .OR.     KEY(1:2) .EQ. 'CH'
     *  .OR.     KEY(1:4) .EQ. 'CALP'
     *  .OR.     KEY(1:4) .EQ. 'RECO'
     *  .OR.     KEY(1:3) .EQ. 'CH0'
     *  .OR.     KEY(1:3) .EQ. 'CH+'
     *  .OR.     KEY(1:2) .EQ. 'CI'  )THEN
          READ(LDUMS,*)IL,IH,(CGRID(JT),JT=1,NTEMP)
          NINTEP=NTEMP
C:
C: 23-FEB-1994 P.G.JUDGE MODIFICATIONS START:
C: SPECIAL CASE OF RAPID TE DEPENDENCE- INTERPOLATE LINEARLY
          IF(KEY(1:3) .EQ. 'CH+' .OR. KEY(1:3) .EQ. 'CH0') THEN
            NINTEP=2
c            WRITE(LJOBLO,*)'LINEAR INTERPOLATION OF LOG10 ',KEY(1:3)
c     *      ,  ' BETWEEN LEVELS', IL, IH
            DO 333 JT=1,NTEMP
              CGRID(JT)=LOG10(CGRID(JT))
  333       CONTINUE
          ENDIF
C:
C: 23-FEB-1994 P.G.JUDGE MODIFICATIONS END:
C:
          IF(IWATOM .GT. 1)THEN
            WRITE(LOUT,802)KEY(1:4), IL, IH, (CGRID(JT),JT=1,NTEMP)
  802       FORMAT(1X,A,2(2X,I2),10(1P,1X,E10.3))
          ENDIF
        ELSE
          WRITE(LJOBLO,191)KEY(1:4)
  191     FORMAT('  GENCOL: KEYWORD(1:4) IS  (',A,')')
          CALL STOP(' GENCOL:  UNKNOWN KEYWORD IN ATOMIC DATA')
        ENDIF
C
C  IDENTIFY THE UPPER AND LOWER LEVELS:
C
        ILO=MIN( IL, IH )
        IHI=MAX( IL, IH )
C
C  CHECK TO SEE THAT THE LEVELS ARE BETWEEN 1 AND NK
C
        IF(ILO .LT. 1 .OR. IHI .GT. NK)
     *  CALL STOP(' GENCOL:  LEVEL INDEX OUTSIDE RANGE (1,NK)')
C
C  CHECK TO SEE IF THERE ARE LEVELS WITH ENERGIES NOT INCREASING WITH
C  THE LEVEL LABEL (THIS AFFECTS CI ETC.)
C
        IF(EV(IHI) .LT. EV(ILO) )THEN
          WRITE(LJOBLO,1001)KEY(1:4),IHI,ILO
 1001     FORMAT('  GENCOL:  WARNING- KEYWORD ',A,2(2X,I2),
     *           ' ENERGY OF UPPER LEVEL IS .LT. LOWER LEVEL'/
     *           '           WARNING- CHECK RESULTS CAREFULLY')
        ENDIF
C
C  DEPTH DO-LOOP:
C
C
C  SEMI-EMPIRICAL COLLISIONS
C
          IF (KEY(1:4)  .EQ.  'SEMI')THEN
            EUPCM=EV(IHI)*EE/CC/HH
            ELOCM=EV(ILO)*EE/CC/HH
            FEM = FAB*G(ILO)/G(IHI)
            CT=SEMIC(ION(ILO),EUPCM,ELOCM,FEM,TEMP,III)
            CDN=CT*NE
            CUP = CDN * NSTAR(IHI) / NSTAR(ILO)
            C(IHI,ILO) = CDN + C(IHI,ILO)
            C(ILO,IHI) = CUP + C(ILO,IHI)
C
C  CORONAL ION BALANCE DATA (FROM DATA READ IN OPACITY PACKAGE)
C
          ELSE IF (KEY(1:6)  .EQ.  'CORONA')THEN
            PE=BK*NE*TEMP
            CALL CORONR(ATOMID(1:3),ION(ILO),TEMP,PE,RECRAT,RIRAT)
            C(IHI,ILO) =  C(IHI,ILO) + RECRAT*NE
            C(ILO,IHI) =  C(ILO,IHI) + RIRAT*NE
          ELSE
            IF(TEMP .LT. TGRID(1)) THEN
              CT=CGRID(1)
            ELSE IF (TEMP .GT. TGRID(NTEMP)) THEN
              CT=CGRID(NTEMP)
            ELSE
              CT=SPLIN(TEMP,TGRID,CGRID,NTEMP,NINTEP)
            ENDIF
          ENDIF
          IF(KEY(1:3) .EQ. 'CH+' .OR. KEY(1:3) .EQ. 'CH0') CT=10.**CT
          IF(CT .LT. 0.) THEN
            WRITE(LJOBLO,*)'INTERPOLATED COLLISION PARAMETER < 0.',
     *        ' AT TEMPERATURE ', TEMP
            WRITE(LJOBLO,*)'FOR TRANSITION ', KEY,' ',IL,' ', IH
            WRITE(LJOBLO,*)(TGRID(JT),JT=1,NTEMP)
            WRITE(LJOBLO,*)(CGRID(JT),JT=1,NTEMP)
            CALL STOP('NEGATIVE COLLISION RATES ENCOUNTERED')
          ENDIF
C
C  OMEGAS ARE GIVEN (+VE IONS)
C
          IF(KEY(1:3) .EQ. 'OHM') THEN
            CDN = 8.63E-06 * CT * NE / ( G(IHI)*SQRT(TEMP) )
            CUP = CDN * NSTAR(IHI) / NSTAR(ILO)
            C(IHI,ILO) = CDN + C(IHI,ILO)
            C(ILO,IHI) = CUP + C(ILO,IHI)
C
C  CE VALUES ARE GIVEN (NEUTRALS)
C
          ELSE IF (KEY(1:2) .EQ. 'CE')THEN
            CDN = NE * CT * G(ILO) * SQRT(TEMP) / G(IHI)
            CUP= CDN * NSTAR(IHI) / NSTAR(ILO)
            C(IHI,ILO) = CDN + C(IHI,ILO)
            C(ILO,IHI) = CUP + C(ILO,IHI)
C
C  CP VALUES ARE GIVEN (B-B COLLISIONS WITH PROTONS)
C
          ELSE IF (KEY(1:2) .EQ. 'CP')THEN
            CDN = NH(6) * CT
            CUP= CDN * NSTAR(IHI) / NSTAR(ILO)
            C(IHI,ILO) = CDN + C(IHI,ILO)
            C(ILO,IHI) = CUP + C(ILO,IHI)
C
C  CH VALUES ARE GIVEN (COLLISIONS WITH NEUTRAL HYDROGEN)
C
          ELSE IF (KEY(1:3) .EQ. 'CH ')THEN
            CDN = NH(1) * CT 
            CUP= CDN * NSTAR(IHI) / NSTAR(ILO)
            C(IHI,ILO) = CDN + C(IHI,ILO)
            C(ILO,IHI) = CUP + C(ILO,IHI)
C
C  CI VALUES ARE GIVEN
C
          ELSE IF (KEY(1:2) .EQ. 'CI')THEN
            DEKT=  (EV(IHI)-EV(ILO)) * EK / TEMP
            CUP = NE * CT * EXP(-DEKT) * SQRT(TEMP)
            CDN = CUP * NSTAR(ILO) / NSTAR(IHI)
            C(IHI,ILO) = CDN + C(IHI,ILO)
            C(ILO,IHI) = CUP + C(ILO,IHI)
C
C  CH OR CH+ VALUES GIVEN:  CHARGE TRANSFER COLLISIONS WITH 
C            HYDROGEN ATOMS/IONS  NB:  ORDERING IMPORTANT HERE
C  EXAMPLE INPUT:
C
C  CH0
C  1  2  1.E-9 1.E-9 1.E-9
C
C  THIS ASSUMES THAT THE RATE 1.E-9*NH(1) WILL BE ADDED TO THE
C  COLLISION RATE FROM LEVEL 1 TO LEVEL 2
C
          ELSE IF (KEY(1:3) .EQ. 'CH0')  THEN
            C(IL,IH)= NH(1) * CT  + C(IL,IH)
          ELSE IF (KEY(1:3) .EQ. 'CH+') THEN
            C(IL,IH)= NH(6) * CT  + C(IL,IH)
C
C  LOW-TEMPERATURE DIELECTRONIC RECOMBINATION COEFFS. ARE GIVEN
C  FORMULA (9) OF NUSSBAUMER AND STOREY (A+A PAPER II)
C
          ELSE IF(KEY(1:4) .EQ. 'LTDR') THEN
            T4LTDR=TEMP/1.E4
            CDN = 1.E-12*NE* (ALTDR/T4LTDR + BLTDR +
     *             CLTDR*T4LTDR +DLTDR*T4LTDR*T4LTDR) *
     *             EXP(-FLTDR/T4LTDR)/T4LTDR/SQRT(T4LTDR)
            C(IHI,ILO) = CDN + C(IHI,ILO)
C
C  ADDITION P. JUDGE 24-JAN-1994: BEGIN
C
C  PARAMETRIC FITS TO DATA OF ARNAUD AND ROTHENFLUG 1985,
C  APJS 60, 425.
C  CDI IS 5 BY MSHELL ELEMENT ARRAY, CONTAINING PARAMETERS
C  I, A,B,C,D OF TABLE 1 OF ARNAUD AND ROTHENFLUG
          ELSE IF(KEY(1:7) .EQ. 'AR85-RR') THEN
            CDN= ARADAR*(TEMP/1.E4)**ETAAR
            C(IHI,ILO) = CDN + C(IHI,ILO)
          ELSE IF(KEY(1:8) .EQ. 'AR85-CDI') THEN
            CUP=0.
            DO 321 J=1,NCDI
              XJ=CDI(1,J)*EE/BK/TEMP
              FXJ=EXP(-XJ)*SQRT(XJ) *(CDI(2,J)+CDI(3,J)*(1.+XJ) +
     *         (CDI(4,J)-XJ*(CDI(2,J)+CDI(3,J)*(2.+XJ)))*FONE(XJ) +
     *         CDI(5,J)*XJ*FTWO(XJ))
              CUP = CUP+ 6.69E-07 / CDI(1,J) / SQRT(CDI(1,J))*FXJ
  321       CONTINUE
            IF(CUP .LT. 0.) CALL STOP('GENCOL: CDI-CUP .LT.0')
            C(ILO,IHI) = CUP*NE + C(ILO,IHI)
            C(IHI,ILO) = CUP*NE*NSTAR(ILO)/NSTAR(IHI)
     *          + C(IHI,ILO)
C
C  ADDITION P. JUDGE 24-JAN-1994: END
C
C:
C: 22-FEB-1994 P.G.JUDGE MODIFICATIONS START:
C:
          ELSE IF(KEY(1:8) .EQ. 'AR85-CEA') THEN
            C(ILO,IHI) = CEA*NE + C(ILO,IHI)
C
C  THIS IS INCORRECT SINCE POPULATION OF UPPER LEVEL IS THE UPPER LEVEL
C  OF THE DOUBLY EXCITED STATE, **, I.E. THIS SHOULD BE MULTIPLIED BY 
C  G(**)/G(IHI) EXP(E(**) - E(IHI)).  A BETTER APPROXIMATION WOULD SEEM TO
C  BE 0. TIMES THIS TO AVOID PROBLEMS. HENCE COMMENT THIS LINE OUT
C  
C            C(IHI,ILO) = CEA*NE*NSTAR(ILO)/NSTAR(IHI)
C     *          + C(IHI,ILO)
C:
C: 22-FEB-1994 P.G.JUDGE MODIFICATIONS END:
C:
C:
C: 08-MAR-1994 P.G.JUDGE MODIFICATIONS START:
C:
C:  AR85-CH    CHARGE TRANSFER RECOMBINATION WITH NEUTRAL HYDROGEN
C:
          ELSE IF(KEY(1:7) .EQ. 'AR85-CH' .AND. KEY(1:8) .NE. 
     *  'AR85-CH+' .AND. KEY(1:8) .NE. 'AR85-CHE') THEN
            IF(TEMP .GE. AR85T1 .AND. TEMP .LE. AR85T2) THEN 
              T4=TEMP/1.E4
              CUP = AR85A * 1.E-9 * T4**AR85B * (1. + AR85C 
     *               * EXP(AR85D*T4))*NH(1)           
              C(IL,IH) = CUP + C(IL,IH)
            ENDIF
C:
C:  AR85-CH+   CHARGE TRANSFER WITH IONIZED HYDROGEN
C:
          ELSE IF(KEY(1:8) .EQ. 'AR85-CH+') THEN 
            IF(TEMP .GE. AR85T1 .AND. TEMP .LE. AR85T2) THEN 
              T4=TEMP/1.E4
              CUP = AR85A * 1.E-9 * T4**AR85B * EXP(-AR85C*T4) 
     *               * EXP(-AR85D*EE/BK/TEMP)*NH(6) 
              C(IL,IH) = CUP + C(IL,IH)
            ENDIF
C: 
C:  AR85-CHE   CHARGE TRANSFER WITH NEUTRAL HELIUM
C:
C:
C: 08-MAR-1994 P.G.JUDGE MODIFICATIONS END:
C:
C
C  COMPUTATIONS FROM TABLES OF SHULL & VAN STEENBERG, 
C  AP J SUPPL  48 P95
C
          ELSE IF(KEY(1:7) .EQ. 'SHULL82') THEN
            CDN=ARADSH*(TEMP/1.E4)**(-XRADSH)  +  
     *        ADISH /TEMP/SQRT(TEMP) * EXP(-T0SH/TEMP)*  
     *        (1.E0+BDISH * (EXP(- T1SH/TEMP)))
         CDN=CDN*NE
         CUP=ACOLSH * SQRT(TEMP) * EXP( -TCOLSH / TEMP) 
     *         / (1.E0 + 0.1 * TEMP / TCOLSH) 
         CUP=CUP*NE
C*
C*  3-BODY RECOMBINATION (HIGH DENSITY LIMIT)
C*
         CDN=CDN+CUP*NSTAR(ILO)/NSTAR(IHI)
C*  
         C(IHI,ILO) = CDN + C(IHI,ILO)
         C(ILO,IHI) = CUP + C(ILO,IHI)
C
C  CALP VALUES ARE GIVEN (COLLISIONS WITH ALPHA PARTICLES)
C  ATOM HAS TO BE HE
C  INITIA IS CALLED FOR FIRST RATE TO READ POPULATIONS
C  FROM RSTRT
C
C  RECOMBINATION COEFFICIENTS ARE GIVEN
C  C(J-I)=NE*CT
C
          ELSE IF (KEY(1:4) .EQ. 'RECO') THEN
            C(IHI,ILO)= NE*CT + C(IHI,ILO)
C
C  SIMPLE ADDITION OF CONSTANT, FOR EXAMPLE FOR PHOTOIONIZATION
C  C(J-I)=CT
C
          ELSE IF (KEY(1:8) .EQ. 'CONSTANT') THEN
            C(IL,IH)= CNST + C(IL,IH)
C
C  MORE KEYWORDS CAN BE ADDED HERE
C
          ENDIF
      GOTO 50
C
 301  CALL STOP(' GENCOL: KEYWORD ''END'' NOT FOUND IN COLLISION FILE')
 300  CONTINUE
      RETURN
      END
C
C***********************************************************************
C
      FUNCTION SEMIC(Z,EUP,ELO,FEM,TEMP,IFLAG)
C
C  RETURNS COLLISION RATE DOWNWARDS USING SEMI-EMPIRICAL GF'S
C
C  P.G. JUDGE, NOVEMBER 1987.
C
C  INPUT:
C         Z      - CHARGE OF ION+1 (E.G. FOR C II Z=2)
C         EUP    - ENERGY OF UPPER LEVEL IN CM-1
C         ELO    - ENERGY OF LOWER LEVEL IN CM-1
C         FEM    - EMISSION OSCILLATOR STRENGTH ( = GF/G(UPPER))
C         TEMP   - TEMPERATURE IN DEGREES KELVIN
C
C  OUTPUT:
C         IFLAG  -  = 1 FOR VAN REGEMORTER APPROXIMATION
C         IFLAG  -  = 2 FOR SHEVELKO       APPROXIMATION
C         SEMIC-  DOWNWARD COLLISION RATE IN CM3 /S
C
      INCLUDE 'PREC'
      INTEGER Z
      DATA C1,RY /1.43882,109737.312/
C
C INITIALIZE
C
      DELTE=EUP-ELO
      BETA=DELTE*C1/TEMP
      SEMIC=0.0
C
C  IFLAG=1 => VAN REGEMORTER APPROXN.
C
      IF(BETA .GE. 0.01)THEN
        IFLAG=1
        SEMIC = 3.2E-7*FEM*(RY/DELTE)**1.5*SQRT(BETA)*P(Z,BETA)
C
C  IFLAG=2 => SHEVELKO APPROX
C
      ELSE IF(BETA .LT. 0.01)THEN
        IFLAG=2
        SEMIC=1.74E-07*FEM*EXP(BETA)/SQRT(BETA*DELTE/RY)
     +  *LOG( ELO / DELTE * SQRT(RY / C1 / TEMP))
      ENDIF
C
      RETURN
      END
C
C **********************************************************************
C
      FUNCTION P(Z,B)
C
C  THERMAL P-FUNCTION FOR SEMIC EMPIRICAL COLLISION RATES
C  REFERENCE: SOBEL'MAN - 'ATOMIC PHYSICS II.  '
C
      INCLUDE 'PREC'
      PARAMETER (N1=10)
      DIMENSION  BETREF(N1),PNREF(N1),PCREF(N1)
      INTEGER Z
      DATA BETREF/-2.0,-1.699,-1.398,-1.0,-0.699,-0.398,0.0
     +  ,0.301,0.602,1.0/
      DATA PNREF/1.160,0.956,0.758,0.493,0.331,0.209
     +  ,0.100,0.063,0.040,0.023/
      DATA PCREF/1.160,0.977,0.788,0.554,0.403
     +  ,0.290,0.214,0.201,0.200,0.200/
C
C  LIMIT OF HIGH TEMPERATURE: KT >> E
C
      IF (B .LT. 0.01)THEN
        P = 0.27566 * (0.577 + LOG(B))
      RETURN
C
C  INTERMEDIATE TEMPS (MOST IMPORTANT FOR EQM PLASMAS)
C  SPLINE INTERPOLATION ONTO LOGB GRID:
C
      ELSE IF(B .GT. 0.01 .AND. B .LT. 10.0)THEN
        BB=LOG10(B)
        IF (Z .EQ. 1)P = SPLIN(BB,BETREF,PNREF,N1,N1)
        IF (Z .GT. 1)P = SPLIN(BB,BETREF,PCREF,N1,N1)
      RETURN
C
C  LIMIT OF LOW TEMPERATURE: KT << E
C
      ELSE  IF(B .GT. 10.0)THEN
        IF(Z .EQ. 1)P = 0.066 / SQRT(B)
        IF(Z .GT. 1)P = 0.200
        RETURN
      ENDIF
      END
C
C **********************************************************************
C
      SUBROUTINE CORONR(CEL,ION,TEMP,PRESS,RRATE,RIRATE)
C
C  CORONAL ION BALANCE CALCULATION
C  NEW ROUTINE FOR USE WITH V2.0 OPACITY PACKAGE AND GENCOL
C
C  INPUT:  ELEMENT IDENTIFIER    CEL (CHARACTER)
C          ION STAGE IDENTIFIER  ION (INTEGER)
C          TEMP   (ELECTRON TEMPERATURE IN K       (RL)
C          PRESS  (ELECTRON PRESSURE IN DYNE/CM2)  (RL)
C
C  OUTPUT: RRATE    RECOMBINATION RATE (CM3/S)
C          RIRATE   IONIZATION RATE    (CM3/S)
C  
C  EXAMPLE:
C          CEL='C  ', ION=1, WILL RETURN 
C          THE RECOMBINATION RATE PER C II ATOM PER ELECTRON
C          AND THE IONIZATION RATE PER C I ATOM PER ELECTRON
C
C  NOTES:
C          COMPUTATIONS FROM TABLES OF SHULL & VAN STEENBERG, 
C          AP J SUPPL  48 P95
C          NO PE-DEPENDENCE, THIS COULD BE ADDED WHEN COMPUTATIONS BECOME
C          AVAILABLE
C          IF THE ION FRACTION IS NOT FOUND THEN THE ROUTINE WILL STOP
C          
C  89-03-28 NEW ROUTINE (PHILIP JUDGE)
C
C: CORONR 92-08-10  MODIFICATIONS (MATS CARLSSON)
C:        INTEGER ARRAYS WITH ELEMENT NAME CHANGED TO CHARACTER ARRAYS
C:
C:        92-09-09  MODIFICATIONS (MATS CARLSSON)
C:        INTEGER IEL IN WRITE STATEMENT CHANGED TO CHARACTER
C:
      INCLUDE 'PREC'
      PARAMETER (MCOR=30)
      INTEGER NCOR,IONCOR(MCOR)
      DOUBLE PRECISION ARADC(MCOR),XRADC(MCOR),ADIC(MCOR),T0C(MCOR),
     * BDIC(MCOR),T1C(MCOR),ACOLC(MCOR),TCOLC(MCOR),TECOR,DECOR
      COMMON /CCORON/ TECOR,DECOR,
     * ARADC,XRADC,ADIC,T0C,BDIC,T1C,ACOLC,TCOLC,NCOR,IONCOR
      CHARACTER*3 CCOR
      COMMON /CCCOR/ CCOR(MCOR)
      COMMON /CLU/ LINPUT,LATOM,LATOM2,LATMOS,LDSCAL,LABUND,LOUT,
     * LTIME,LRSTRT,LDSCA2,LWMAT,LNIIT,LDUMS,LDUMI,LDUMC,LOPC,
     * LJNY,LINIT,LPHI,LJOBLO,LATHSE
      CHARACTER*(*) CEL
C
C  IDENTIFY THE ELEMENT AND IONIZATION STAGE TO BE COMPUTED
C
      DO 100 N=1,NCOR
        IF( CEL .EQ. CCOR(N) .AND. ION .EQ. IONCOR(N) )GOTO 99
  100 CONTINUE
      WRITE(LJOBLO,1024) CEL,ION
 1024  FORMAT(' CORONR: SPECIES',A3,I3,' NOT FOUND IN ABSDAT FILE')
      CALL STOP(' ')
C
C  ELEMENT HAS BEEN COMPUTED- NOW RETURN THE RATES
C
   99 RRATE=ARADC(N)*(TEMP/1.E4)**(-XRADC(N))  +     
     *  ADIC(N) * (TEMP**(-1.5)) * EXP(-T0C(N)/TEMP)*
     *  (1.E0+BDIC(N) * (EXP(- T1C(N)/TEMP)))
      RIRATE=ACOLC(N) * SQRT(TEMP) * EXP( -TCOLC(N) / TEMP)
     *    / (1.E0 + 0.1 * TEMP / TCOLC(N))
      RETURN
      END
