
c       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        subroutine COLCAT (temp, xnhmp, xnhmo)
c       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        use constantes
        use auxiliaries
        USE CHEM_DATA
        USE STR_DATA
        USE SP_DATA
        USE PROFIL

        implicit double precision (a-h,o-z)

        REAL (kind=DP), DIMENSION (nlvcat,nlvcat) :: ec_H, ec_He
        REAL (kind=DP), DIMENSION (nlvcat,nlvcat) :: ec_Hp, ec_e
        REAL (kind=DP), DIMENSION (nlvcat,nlvcat) :: ec_orthoH2
        REAL (kind=DP), DIMENSION (nlvcat,nlvcat) :: ec_paraH2, ec

c    a-collisions avec H
c Fine structure: probably Launay & Roueff 1977 (check)

        ec_H(1,2) = 1.01e-10_DP * temp**0.117_DP
        ec_H(1,3) = 4.49e-11_DP * temp**0.194_DP
        ec_H(2,3) = 1.06e-10_DP * temp**0.234_DP

c 1D - 3P
c DRF : Do not use O I + H rates

c  b-collisions avec H2
c fit from Schroder, Staemmler, Smith, Flower & Jaquet (1991)

        ec_paraH2(1,2)  = 0.80e-10_DP
        ec_orthoH2(1,2) = 0.75e-10_DP
        ec_paraH2(1,3)  = 0.90e-10_DP
        ec_orthoH2(1,3) = 3.54e-11_DP * temp**0.167_DP
        ec_paraH2(2,3)  = 2.00e-10_DP
        ec_orthoH2(2,3) = 5.25e-11_DP * temp**0.244_DP

c     c-collisions avec He
c Fine structure: Staemmler & Flower (1991)

        ec_He(1,2) = 8.38e-12_DP * temp**0.159_DP
        ec_He(1,3) = 5.98e-11_DP * temp**0.078_DP
        ec_He(2,3) = 3.68e-11_DP * temp**0.041_DP

c collisions avec H+
c from AA 236, 515 (1990)  Roueff & Le Bourlot (fit DRF)

        ec_Hp(1,2) = 10._DP**(-10.359_DP + 0.7959_DP*log10(temp)
     1             - 0.08748_DP*(log10(temp))**2._DP)
        ec_Hp(1,3) = 10._DP**(-13.232_DP + 2.4171_DP*log10(temp)
     1             - 0.29151_DP*(log10(temp))**2._DP)
        ec_Hp(2,3) = 10._DP**(-11.290_DP + 1.7915_DP*log10(temp)
     1             - 0.23010_DP*(log10(temp))**2._DP)

c     d-collisions with electrons
c Should we use Pequignot 1990...
c From Pequignot & Aldrovandi, A&A 50, 141, 1976 as compiled by Mendoza (1983)

        auxt4 = 1.83e-4_DP * temp ** 0.444_DP
        ec_e(1,4) = 8.629e-6_DP * auxt4 / gstcat(4) * gstcat(1) / 9.0_DP
        ec_e(2,4) = 8.629e-6_DP * auxt4 / gstcat(4) * gstcat(2) / 9.0_DP
        ec_e(3,4) = 8.629e-6_DP * auxt4 / gstcat(4) * gstcat(3) / 9.0_DP
        auxt4 = 9.86e-5_DP * temp ** 0.343_DP
        ec_e(1,5) = 8.629e-6_DP * auxt4 / gstcat(5) * gstcat(1) / 9.0_DP
        ec_e(2,5) = 8.629e-6_DP * auxt4 / gstcat(5) * gstcat(2) / 9.0_DP
        ec_e(3,5) = 8.629e-6_DP * auxt4 / gstcat(5) * gstcat(3) / 9.0_DP
        auxt4 = 2.77e-3_DP
        ec_e(4,5) = 8.629e-6_DP *auxt4 / gstcat(5)

c  Bring'em together

       do jl=1,nlvcat-1
         do ju=jl+1,nlvcat
           ec(jl,ju) = ec_H(jl,ju) * ab(i_h)
     1               + ec_Hp(jl,ju) * ab(i_hp)
     2               + ec_paraH2(jl,ju) * xnhmp
     3               + ec_orthoH2(jl,ju) * xnhmo
     4               + ec_He(jl,ju) * ab(i_he)
     5               + ec_e(jl,ju) * ab(nspec)
         end do
       end do

       do ju=2,nlvcat
         do jl = 1,ju-1
           ec(ju,jl) = ec(jl,ju) * (gstcat(ju)/gstcat(jl))
     1               * exp(-(elkcat(ju)-elkcat(jl)) / temp)
         end do
       end do

       eqcol = zero

       do jl=1,nlvcat-1
         do ju=jl+1,nlvcat
           eqcol(jl,ju) = eqcol(jl,ju) - ec(jl,ju)
           eqcol(ju,ju) = eqcol(ju,ju) + ec(jl,ju)
         end do
       end do

       do ju=2,nlvcat
         do jl = 1,ju-1
           eqcol(ju,jl) = eqcol(ju,jl) - ec(ju,jl)
           eqcol(jl,jl) = eqcol(jl,jl) + ec(ju,jl)
         end do
       end do

       return
       end

c       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        subroutine COLOAT (temp, xnhmp, xnhmo)
c       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        use constantes
        use auxiliaries
        USE CHEM_DATA
        USE STR_DATA
        USE SP_DATA
        USE PROFIL

        implicit double precision (a-h,o-z)

        REAL (kind=DP), DIMENSION (nlvoat,nlvoat) :: ec_H, ec_He
        REAL (kind=DP), DIMENSION (nlvoat,nlvoat) :: ec_Hp, ec_e
        REAL (kind=DP), DIMENSION (nlvoat,nlvoat) :: ec_orthoH2
        REAL (kind=DP), DIMENSION (nlvoat,nlvoat) :: ec_paraH2, ec

c    collisions avec H
c  fit refait d apres Launay & Roueff

       ec_H(1,2) = 4.37d-12 * temp**0.66d0
       ec_H(1,3) = 1.06d-12 * temp**0.80d0
       ec_H(2,3) = 1.35d-11 * temp**0.45d0

c 1D - 3P
c fit by JLB, from datas of Federman & Shipsey, Ap J, 1983, 269, 791 (using potential V2)

       ec_H(1,4) = 9.85e-14_DP * temp**0.245_DP
       ec_H(2,4) = 9.85e-14_DP * temp**0.245_DP
       ec_H(3,4) = 9.85e-14_DP * temp**0.245_DP

c  collisions avec H2
c  fit from Jaquet, Staemmler, Smith & Flower, J.Phys.B (1992)

       ec_paraH2(1,2) = 3.46d-11 * temp**0.316d0
       ec_orthoH2(1,2) = 2.70d-11 * temp**0.362d0
       ec_paraH2(1,3) = 7.07d-11 * temp**0.268d0
       ec_orthoH2(1,3) = 5.49d-11 * temp**0.317d0
       ec_paraH2(2,3) = 1.44d-14 * temp**1.109d0
       ec_orthoH2(2,3) = 4.64d-14 * temp**0.976d0

c     collisions avec He
c fit from Monteiro & Flower, MNRAS 228, 101 (1987)

       ec_He(1,2) = 1.55d-12 * temp**0.694d0
       ec_He(1,3) = 2.52d-12 * temp**0.681d0
       ec_He(2,3) = 2.00d-15 * temp**1.428d0

c    collisions avec H+
c Pequignot A&A 231, 499 (1990) + erratum A&A 313, 1026 (1996)

        if (temp < 194.664_dp) then
          toto = 1.40e-3_dp * (temp*1.0e-2_dp)**0.90_dp
        else if (temp < 3686.2414_dp) then
          toto = 2.14e-2_dp * (temp*1.0e-3_dp)**1.30_dp
        else
          toto = 2.78e-1_dp * (temp*1.0e-4_dp)**0.87_dp
        endif
        ec_HP(1,2) = toto * 8.629e-6_dp / (sqrt(temp) * gstoat(2))

        if (temp < 511.9563_dp) then
          toto = 1.12e-4_dp * (temp*1.0e-2_dp)**1.60_dp
        else if (temp < 7510.155_dp) then
          toto = 3.90e-3_dp * (temp*1.0e-3_dp)**1.40_dp
        else
          toto = 8.25e-2_dp * (temp*1.0e-4_dp)**0.80_dp
        endif
        ec_Hp(1,3) = toto * 8.629e-6_dp / sqrt(temp)

        if (temp < 2090.2558_dp) then
          toto = 3.10e-4_dp * (temp*1.0e-2_dp)**1.06_dp
        else
          toto = 2.29e-2_dp * (temp*1.0e-4_dp)**0.69_dp
        endif
        ec_Hp(2,3) = toto * 8.629e-6_dp / sqrt(temp)

c     d-collisions with electrons
c      From Pequignot, 1990, A&A, 231, 499 and Berrington, 1998, MNRAS, 293, L83

       auxt4 = 4.28e-1_DP * (temp/1.e4_DP) ** 1.43_DP
     1       / (1.00_DP + 6.05e-1_DP * (temp/1.e4_DP) ** 1.105_DP)
       ec_e(1,4) = 8.629e-6_DP * auxt4 / gstoat(4)
     1       / sqrt(temp) * gstoat(1) / 9.0_DP
       ec_e(2,4) = 8.629e-6_DP * auxt4 / gstoat(4)
     1       / sqrt(temp) * gstoat(2) / 9.0_DP
       ec_e(3,4) = 8.629e-6_DP * auxt4 / gstoat(4)
     1       / sqrt(temp) * gstoat(3) / 9.0_DP
       auxt4 = 5.88e-2_DP * (temp/1.e4_DP) ** 1.50_DP
     1       / (1.00_DP + 8.00e-1_DP * (temp/1.e4_DP) ** 1.125_DP)
       ec_e(1,5) = 8.629e-6_DP * auxt4 / gstoat(5)
     1       / sqrt(temp) * gstoat(1) / 9.0_DP
       ec_e(2,5) = 8.629e-6_DP * auxt4 / gstoat(5)
     1       / sqrt(temp) * gstoat(2) / 9.0_DP
       ec_e(3,5) = 8.629e-6_DP * auxt4 / gstoat(5)
     1       / sqrt(temp) * gstoat(3) / 9.0_DP
       auxt4 = 1.16e-1_DP * (temp/1.e4_DP) ** 0.53_DP
     1       / (1.00_DP + 1.11e-1_DP * (temp/1.e4_DP) ** 0.160_DP)
       ec_e(4,5) = 8.629e-6_DP *auxt4 / gstoat(5) / sqrt(temp)

       auxt4 = 1.49e-4_DP * temp ** 0.4565_DP
       ec_e(1,2) = 8.629e-6_DP * auxt4 / gstoat(2) / sqrt(temp)
       auxt4 = 4.98e-5_DP * temp ** 0.4955_DP
       ec_e(1,3) = 8.629e-6_DP * auxt4 / gstoat(3) / sqrt(temp)
       auxt4 = 1.83e-9_DP * temp ** 1.347_DP
       ec_e(2,3) = 8.629e-6_DP * auxt4 / gstoat(3) / sqrt(temp)

c  Bring'em together

       do jl=1,nlvoat-1
         do ju=jl+1,nlvoat
           ec(jl,ju) = ec_H(jl,ju) * ab(i_h)
     1               + ec_Hp(jl,ju) * ab(i_hp)
     2               + ec_paraH2(jl,ju) * xnhmp
     3               + ec_orthoH2(jl,ju) * xnhmo
     4               + ec_He(jl,ju) * ab(i_he)
     5               + ec_e(jl,ju) * ab(nspec)
         end do
       end do

       do ju=2,nlvoat
         do jl = 1,ju-1
           ec(ju,jl) = ec(jl,ju) * (gstoat(ju)/gstoat(jl))
     1               * exp(-(elkoat(ju)-elkoat(jl)) / temp)
         end do
       end do

       eqcol = zero

       do jl=1,nlvoat-1
         do ju=jl+1,nlvoat
           eqcol(jl,ju) = eqcol(jl,ju) - ec(jl,ju)
           eqcol(ju,ju) = eqcol(ju,ju) + ec(jl,ju)
         end do
       end do

       do ju=2,nlvoat
         do jl = 1,ju-1
           eqcol(ju,jl) = eqcol(ju,jl) - ec(ju,jl)
           eqcol(jl,jl) = eqcol(jl,jl) + ec(ju,jl)
         end do
       end do

       return
       end

c       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        subroutine COLCPL (temp, xnhmp, xnhmo)
c       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        use constantes
        use auxiliaries
        USE STR_DATA
        USE CHEM_DATA
        USE SP_DATA
        USE PROFIL

        implicit double precision (a-h,o-z)

        REAL (kind=DP), DIMENSION (nlvcpl,nlvcpl) :: ec_H, ec_He
        REAL (kind=DP), DIMENSION (nlvcpl,nlvcpl) :: ec_Hp, ec_e
        REAL (kind=DP), DIMENSION (nlvcpl,nlvcpl) :: ec_orthoH2
        REAL (kind=DP), DIMENSION (nlvcpl,nlvcpl) :: ec_paraH2, ec


c    a-collisions avec H

        ec_H(1,2) = 8.86e-10_dp

c  b-collisions avec H2

        ec_paraH2(1,2) = 3.94e-10_dp
        ec_orthoH2(1,2) = 4.92e-10_dp

c     c-collisions avec electrons

c k21 : Hayes & Nussbaumer, 1984, A&A 134, 193
c Dec 02, CM & DRF
c k21 changed to power law fit to Wilson & Bell, 2002, MNRAS, 337, 1027
c other : Lennon et al., 1985, Ap J 294, 200

        toto = 0.65582_dp * temp**(0.13244_DP)
        ec_e(1,2) = 8.629e-6_DP * min(toto,2.5_dp)
     1            / (4.0_DP*sqrt(temp))
        ec_e(1,3) = 8.629e-6_DP * (0.288_DP - 5.89e-7_DP*temp)
     1            / (2.0_DP * sqrt(temp))
        ec_e(1,4) = 8.629e-6_DP * (0.424_DP - 8.35e-7_DP*temp)
     1            / (4.0_DP * sqrt(temp))
        ec_e(1,5) = 8.629e-6_DP * (0.257_DP - 3.95e-7_DP*temp)
     1            / (6.0_DP * sqrt(temp))
        ec_e(2,3) = 8.629e-6_DP * (0.197_DP - 3.15e-7_DP*temp)
     1            / (2.0_DP * sqrt(temp))
        ec_e(2,4) = 8.629e-6_DP * (0.544_DP - 9.84e-7_DP*temp)
     1            / (4.0_DP * sqrt(temp))
        ec_e(2,5) = 8.629e-6_DP * (1.200_DP - 2.40e-6_DP*temp)
     1            / (6.0_DP * sqrt(temp))

c  Bring'em together

       do jl=1,nlvcpl-1
         do ju=jl+1,nlvcpl
           ec(jl,ju) = ec_H(jl,ju) * ab(i_h)
     1               + ec_Hp(jl,ju) * ab(i_hp)
     2               + ec_paraH2(jl,ju) * xnhmp
     3               + ec_orthoH2(jl,ju) * xnhmo
     4               + ec_He(jl,ju) * ab(i_he)
     5               + ec_e(jl,ju) * ab(nspec)
         end do
       end do

       do ju=2,nlvcpl
         do jl = 1,ju-1
           ec(ju,jl) = ec(jl,ju) * (gstcpl(ju)/gstcpl(jl))
     1               * exp(-(elkcpl(ju)-elkcpl(jl)) / temp)
         end do
       end do

       eqcol = zero

       do jl=1,nlvcpl-1
         do ju=jl+1,nlvcpl
           eqcol(jl,ju) = eqcol(jl,ju) - ec(jl,ju)
           eqcol(ju,ju) = eqcol(ju,ju) + ec(jl,ju)
         end do
       end do

       do ju=2,nlvcpl
         do jl = 1,ju-1
           eqcol(ju,jl) = eqcol(ju,jl) - ec(ju,jl)
           eqcol(jl,jl) = eqcol(jl,jl) + ec(ju,jl)
         end do
       end do

        return
        end

