
!#include "REAL.H"
#include "CO_Cool_globals.H"

module CO_CoolGlobal

  
  implicit none
  
  integer, parameter :: NNEQ = 1
  integer, parameter :: eint = 1
  
  integer, parameter :: ie = 1
  integer, parameter :: iH = 2
  integer, parameter :: iH2 = 3
  integer, parameter :: iHe = 4
  integer, parameter :: iHplus = 5
  integer, parameter :: iCO = 6
  integer, parameter :: iO = 7
  integer, parameter :: iC = 8
  integer, parameter :: iCplus = 9
  integer, parameter :: NTOT = 9

  integer, parameter :: NPRATES = 100
  integer, parameter :: NCRATES = 500
  integer, parameter :: NTRATES = 500
  integer, parameter :: NTEMPR = 2000
  
  real :: temprv(NTEMPR)
  real :: tratesv(NTRATES,NTEMPR), tratesd(NTRATES,NTEMPR)
  real :: tcoeff

  ! min and max temperatures for tabulated thermal rates
  real, parameter :: tmin = 5.0
  real, parameter :: tmax = 1.0e5

  real :: ATOL(NNEQ)
  
  
  !integer, parameter :: nrpar = 30
  !integer, parameter :: nipar = 30
  
  
  real, parameter :: SAFETY  =  0.05
  
  real, parameter ::  YHELIUM =  0.08
  real, parameter ::  YOXYGEN =  3.2e-4
  real, parameter ::  YCARBON =  1.4e-4
  
  real, parameter ::  zeta   =   1.0e-16
  real, parameter ::  alpha_gd = 3.2e-34
  
  ! should always be false
  logical, parameter :: evolve_temp = .false.

  !! do we include line opacity for fine-structure line cooling?
  !! potentially begins to become important above n~10^3 cm^-3
  logical, parameter :: fine_structure_line_opacity = .false.
  
  !! if true, compute CO cooling under LVG approximation
  !! using NCO_tilde = n_co / |divergence(v)| as shielding
  !! parameter
  !! if false, CO cooling uses CO column density
  logical, parameter :: CO_LVG = .true.
  
  !! if true we compute the shielding length as rho / |(grad(rho))|
  !! and compute column densities from that.
  !! if false, column densities are passed in directly in rpar(2:4)
  logical, parameter :: density_gradient_shielding_length = .true.

  !! if simple_cooling = .true.
  !! we use the exact cooling function from Kim & Ostriker 2013.
  !! shielding is *not* included
  logical, parameter :: simple_cooling = .true.

  !! include a CMB temperature floor?
  logical, parameter :: temp_floor = .true.
  
  real, parameter ::  TCMB  =    2.725
  
  real, parameter ::  kb = 1.38e-16
  real, parameter ::  mh = 1.67e-24

  real :: ei0, rpar0(nrpar), dt0

end module CO_CoolGlobal


! integrate chemical abundances and an internal energy (y) over a given timestep
subroutine COOL_STEP(dt, ei, rpar)
  
  use CO_CoolGlobal
  
  implicit none
  
  real, intent(INOUT) :: ei, rpar(nrpar)
  real, intent(IN) :: dt
  
  integer :: ipar(nipar)
  real :: yneq(NNEQ)
  real :: yneq_previous(NNEQ)
  
  ! MF = 22
  integer, parameter :: LRW = 22 + 9*NNEQ + 2*NNEQ**2
  integer, parameter :: LIW = 30 + NNEQ
  ! MF = 10
  !integer, parameter :: LRW = 20 + 16*NNEQ
  !integer, parameter :: LIW = 30
  
  integer :: ITOL, RTOL(1)
  integer :: ITASK, ISTATE, IOPT, MF, IWORK(LIW)
  real :: RWORK(LRW), TIN
  
  real :: h, t
  integer :: nsub
  logical :: reached_the_end

  integer :: i, n, nback
  
  logical :: backtrack

  external coolderiv, jac
  
  !! rpar parameters:
  !!
  !! rpar(1)  = nh
  !! rpar(2)  = total H column density
  !! rpar(3)  = total H2 column density
  !! rpar(4)  = total CO column density
  !! rpar(5)  = unattenutated radiation field in units of draine field (G0) 
  !! rpar(6)  = dust temperature
  !! rpar(7)  = magnitude of density gradient
  !! rpar(8)  = absolute value of velocity divergence
  !! rpar(9)  = 
  !! rpar(10) = sign for error
  !! rpar(11) = 
  !! rpar(12) = 
  !! rpar(13) = 

  rpar0 = rpar
  ei0 = ei
  dt0 = dt

  yneq(1) = ei

  rpar(10) = 0.0
  t = 0.0
  nsub = 0
  reached_the_end = .false.
  do while(.not. reached_the_end)
     
     nsub = nsub + 1
     nback = 0
     
     call set_timestep(yneq, rpar, h)
     
     !! reduce timestep if nearing end
     reached_the_end = (t + h) >= dt
     if(reached_the_end) h = dt - t
    
     ITOL = 2
     RTOL(1) = 1.0e-3
     ITASK = 1
     ISTATE = 1
     IOPT = 0
     MF = 22
     ! non stiff method flag:
     !MF = 10
   
     ATOL(eint) = 1.0e-14
     
     backtrack = .false.
     
     do
      
        yneq_previous = yneq
     	TIN = 0.0
        call dvode(coolderiv, NNEQ, yneq, TIN, h, ITOL, RTOL, ATOL, ITASK, ISTATE, &
             IOPT, RWORK, LRW, IWORK, LIW, jac, MF, rpar, ipar)

        if(ISTATE .lt. 0 .or. rpar(10) .gt. 0.0) then
           backtrack = .true.
        endif
        
        if(.not. backtrack) exit
        
        if(backtrack) then
           !print*, "backtracking...", nback, ISTATE, rpar(10), yneq(1)
           ISTATE = 1
           rpar(10) = 0.0
           backtrack = .false.
           yneq = yneq_previous
           h = 0.1 * h
           reached_the_end = .false.
           nback = nback + 1
        endif
        
        if(nback .gt. 5) then
           print*, "too many backtrackings!", nback
           print*, "original ei = ", ei0
           print*, "original rpar = ", rpar0
           print*, "global timestep = ", dt0
           stop
        endif
        
     enddo
     
     do n=1, NNEQ
        if(yneq(n) .lt. 0.0 .or. yneq(n) .ne. yneq(n)) then
           print*, "co_cooling: dvode produced negative / nan temperature. not right..."
           print*, "original ei = ", ei0
           print*, "original rpar = ", rpar0
           print*, "global timestep = ", dt0
           stop
        endif
     enddo
     
     t = t + h
     
  enddo
  
  ei = yneq(1) 
  
       
  return
end subroutine COOL_STEP


subroutine coolderiv(N,t,yneq,ydot,rpar,ipar)

  use CO_CoolGlobal

  implicit none


  integer :: N
  real :: t, yneq(NNEQ), ydot(NNEQ)
  real :: y(NTOT), ei
  real :: rpar(nrpar)
  integer :: ipar(nipar)
  real :: nh, tdust, temperature
  real :: heating, cooling
  real :: co_frac, c_frac
  real :: temp_min, temp_max, temp
  integer :: ntemp, i

  ei = yneq(1)
  nh = rpar(1)
  
  if(ei .le. 0.0) then
     rpar(10) = 1.0
     ydot(eint) = 0.0
     return
  endif
  
  
  !! assuming fully molecular H2 gas (very low neutral H abundance), hard code the abundances here:
  !! assume no ionized fraction to compute first temperature:
  y(iH) = 1.0e-5
  y(ie) = 6.3e-7 * sqrt(zeta / 1.0e-16)  / sqrt(nh / 1.0e2)
  y(iH2) = 0.5 * (1.0 - y(iH))
  y(iHe) = YHELIUM
  call calculate_temperature(y, ei, nh, temperature)

  !! as for carbon, let's assume a transition from C+ to CO smoothly
  !! occurring between 50 and 150 K

  call smoothstep(15.0,40.0,temperature,c_frac)
  co_frac = 1.0-c_frac

 
  y(iCO) = min(YOXYGEN, YCARBON) * co_frac
  
  y(iCplus) = YCARBON - y(iCO)
  
  y(iO) = YOXYGEN - y(iCO)
 
  if(y(iCO) .lt. 0.0) then
     print*, "negative CO!", y(iCO), temperature
     PRINT*, "yoxygen = ", YOXYGEN
     PRINT*, "ycarbon = ", YCARBON
     PRINT*, "co_frac = ", co_frac
     PRINT*, "original ei = ", ei0
     PRINT*, "original rpar = ", rpar0
     PRINT*, "global timestep = ", dt0
     stop
  endif

 if(y(iCplus) .lt. 0.0) then
     print*, "negative Cplus!", y(iCplus), temperature
     stop
  endif

 if(y(iO) .lt. 0.0) then
     print*, "negative O!", y(iO), temperature
     stop
  endif

  y(iHplus) = y(ie) - y(iCplus)
  
  y(iH2) = 0.5 * (1.0 - y(iHplus) - y(iH) )
  
  y(iHe) = YHELIUM

  call get_tdust(y, nh, rpar, tdust)
  
  call thermal(y, temperature, tdust, rpar, cooling, heating)
  ydot(eint) = heating - cooling

  return
end subroutine coolderiv



subroutine thermal_rates(temp,trates)

  use CO_CoolGlobal
  
  implicit none
  
  real, intent(IN) :: temp
  real, intent(OUT) :: trates(NTRATES)
  


    real :: XH1, XH2, XH3, XH4, XH5
    real :: T2, T3, T5, T6, temp_var, sqrt_temp, sqrt_T5, temp_sq
    real :: gff, logtemp3, log_10_temp
    real :: aux_var, aux_var_1, aux_var_2
    real :: ln_temp, ln_temp_2, ln_temp_3, ln_temp_4
    real :: templog, templog_sq, thd
    real :: opratio, fortho, fpara
    real :: ortho, para
    
    real :: ln_ten, sqrt_minus_5
    
    
    ln_ten = log(10.0)
    sqrt_minus_5 = sqrt(1.0e-5)
 

    T2 = 1.0e-2 * temp
    T3 = 1.0e-3 * temp
    T5 = 1.0e-5 * temp
    T6 = 1.0e-6 * temp

    if (temp.le.0.0) then
       print *,'CALC_THERMO_RATE: less than zero!!!!! ', temp
    end if

    ln_temp = log(temp)
    
    ln_temp_2 = ln_temp * ln_temp
    ln_temp_3 = ln_temp * ln_temp_2
    ln_temp_4 = ln_temp * ln_temp_3

    log_10_temp = log10(temp)
 
    sqrt_temp = sqrt(temp)
 
    temp_sq = temp**2

    sqrt_T5 = sqrt_temp * sqrt_minus_5

    XH2 = 1.0 / (1.0 + sqrt_T5)
    XH1 = sqrt_temp * XH2

    aux_var = temp**(-0.397) 
    XH3 = XH2 * aux_var
    aux_var = T6**0.7
    XH4 = 1.0 / (1.0 + aux_var)
    
    aux_var_1 = exp(-470000.0 / temp)
    aux_var_2 = exp(-94000.0 / temp)
    XH5 = aux_var_1 * (1.0 + 0.3 * aux_var_2)

    T3 = 0.001 * temp
 
    templog = log10(temp)
    templog_sq = templog * templog
    
    opratio = 2.4
    fortho = opratio / (1.0 + opratio)
    fpara = 1.0 / (1.0 + opratio)
    
    logtemp3 = log10(temp/1.0e3)
 
    
    !!!!!!!!!!!!!!!!
    ! Atomic processes
    !!!!!!!!!!!!!!!!
    
    
    trates(:) = 0.0
    
! Bremsstrahlung cooling

    temp_var = 5.5 - log_10_temp
  
    aux_var = exp(- temp_var**2 / 3.0)
    gff = 1.1 + 0.34 * aux_var
     
    trates(1) =  1.42e-27 * gff * sqrt_temp
    
    
! Collisional ionization cooling of HI
    if (temp > 2000.0) then
       trates(2) = 1.27e-21 * XH1 * exp(-157809.1 / temp) 
    endif



! Collisional ionization cooling of HeI
    if (temp > 3000.0) then
       trates(3) = 9.38e-22 * XH1 * exp(-285335.4 / temp)
    endif



    ! Collisional ionization cooling of HeII
    if (temp > 4000.0) then
       trates(4) = 4.95e-22 * XH1 * exp(-631515.0 / temp)
    endif


    
    ! Recombination cooling of HII
    
    trates(5) = 8.7e-27 * XH4 * sqrt_temp * T3**(-0.2) 
    

    
    ! Recombination cooling of HeII
    !    if (temp > 5000.0) then
    
    trates(6) = 1.55e-26 * temp**0.3647
    
    

    
    ! Recombination cooling of HeIII
    !    if (temp > 5000.0) then
    
    trates(7) = 3.48e-26 * XH4 * sqrt_temp * T3**(-0.2)
    

    ! Dielectronic recombination cooling
    !    if (temp > 5000.0) then
    if (temp > 3000.0) then
       
       trates(8) = 1.24e-13 * XH5 / (temp * sqrt_temp)
       
    end if


    ! Collisional excitation cooling of HI
    !    if (temp > 5000.0) then
    if (temp > 900.0) then
       
       trates(9) = 7.5e-19 * XH2 * exp(-118348.0 / temp)
       
    endif


    ! Collisional excitation cooling of HeII
    !    if (temp > 8000.0) then
    if (temp > 3000.0) then
       
       trates(10) = 5.54e-17 * XH3 * exp(-473638.0 / temp)
       
    endif
    



    !!!!!!!!!!!!!!!!
    ! Molecular processes
    !!!!!!!!!!!!!!!!
    

    ! H2  (milos's code)
    
    
    ! low density H2 cooling function:
    
    trates(11) = 10.0**( -103.0 + 97.59 * templog - 48.05 * templog_sq + &
            10.8 * templog_sq * templog - 0.9032 * templog_sq * templog_sq )

    ! high density (LTE) H2 cooling function:

    ! factor_1
    trates(12) = 9.5e-22 * T3**3.76
    
    ! factor_2 (dens dependance)
    trates(13) = exp(- (0.13 / T3)**3)
    
    ! factor_3
    trates(14) = (1.0 + 0.12 * T3**2.1)
    
    ! factor_4 (dens dependance)
    trates(15) = 3.0e-24 * exp(- 0.51 / T3)

    ! LVLTE 
    trates(16) = 6.7e-19 * exp(-5.86 / T3) + 1.6e-18 * exp(-11.7 / T3)
    

    ! HD (Lipovka implementation):
    !! Not currently being used in thermo_rates...
    
    if(temp .gt. 2.0e4) then
       thd = log10(2.0e4)
    else
       thd = templog
    end if
    
   trates(17) =  -42.57688 & 
         + 21.93385 * thd & 
         - 10.19097 * thd**2 &
         + 2.19906  * thd**3 &
         - 0.17334  * thd**4 
    
    trates(18) =  0.92433 & 
         + 0.77952 * thd &
         - 0.54263 * thd**2 &
         + 0.11711 * thd**3 &
         - 8.35d-3 * thd**4
    
    trates(19) =  0.54962 & 
         - 1.06447 * thd &
         + 0.62343 * thd**2 &
         - 0.13768 * thd**3 &
         + 1.06d-2 * thd**4 
    
    trates(20) = - 7.676d-2 &
         + 0.11864d0 * thd &
         - 7.366d-2  * thd**2 &
         + 1.759d-2  * thd**3 &
         - 1.482d-3  * thd**4
    
    trates(21) =  2.75d-3 &
         - 3.66d-3   * thd &
         + 2.514d-3  * thd**2 &
         - 6.6631d-4 * thd**3 &
         + 6.1926d-5 * thd**4



    !!!!!!!!!!!!!!!!
    ! Metal de-excitation rates.
    ! Mostly from Glover & Jappsen (2007) :
    !!!!!!!!!!!!!!!!
    
    
    ! CII:
    
    ! collisions with e-
    if(temp .le. 2.0e3) then
       trates(22) = 3.86e-7 * T2**(-0.5)
    else
       trates(22) = 2.43e-7 * T2**(-0.345)
    end if

    ! collisions with H
    if(temp .le. 2.0e3) then
       trates(23) = 8.0e-10 * T2**(0.07) 
    else
       trates(23) = 3.1e-10 * T2**(0.385)
    end if
    
    ! collisions with H2
    if(temp .le. 250.0) then
       ortho = (4.7e-10 + 4.6e-13 * temp) * fortho
       para = 2.5e-10 * temp**0.12 * fpara
    else
       !ortho = 5.85e-10 * temp**0.07 * fortho
       !para = 4.85e-10 * temp**0.07 * fpara
       ! to make rates match at 250:
       ortho = 3.975e-10 * temp**0.07 * fortho
       para = 3.295e-10 * temp**0.07 * fpara
    end if
    trates(24) = ortho + para
    

!!!!!!!!!!!!!
!!!!!!!!!!!!!
!!!!!!!!!!!!!


    ! OI:
    
    ! collisions with H2
    
    
    ! 1->0
    ortho = 2.7e-11 * temp**0.362 * fortho
    para = 3.46e-11 * temp**0.316 * fpara
    trates(25) = ortho + para
    
    ! 2->0
    ortho = 5.49e-11 * temp**0.317 * fortho 
    para = 7.07e-11 * temp**0.268 * fpara
    trates(26) = ortho + para

    ! 2->1
    ortho = 2.74e-14 * temp**1.060 * fortho
    para = 3.33e-15 * temp**1.360 * fpara
    trates(27) = ortho + para
    

    ! collisions with H

    ! 1->0
    trates(28) = 9.2e-11 * T2**0.67
    
    ! 2->0
    trates(29) = 4.3e-11 * T2**0.80

    ! 2->1
    trates(30) = 1.1e-10 * T2**0.44
    
    
    ! collisions with e-

    ! 1->0
    trates(31) = 5.12e-10 * temp**(-0.075)
    
    ! 2->0
    trates(32) = 4.86e-10 * temp**(-0.026)

    ! 2->1
    trates(33) = 1.08e-14 * temp**(0.926)


    ! collisions with H+

    ! 1->0
    if(temp .le. 194.0) then
       trates(34) = 6.38e-11 * temp**0.4
    else if(temp .le. 3686.0) then
       trates(34) = 7.75e-12 * temp**0.8 
    else
       trates(34) = 2.65e-10 * temp**0.37
    end if
    
    ! 2->0
    if(temp .le. 511.0) then
       trates(35) = 6.10e-13 * temp**1.1
    else if(temp .le. 7150.0) then
       trates(35) = 2.12e-12 * temp**0.9
    else
       trates(35) = 49e-10 * temp**0.3
    end if
 

    ! 2->1
    if(temp .le. 2090.0) then
       trates(36) = 2.03e-11 * temp**0.56
    else
       trates(36) = 3.43e-10 * temp**0.19
    end if




!!!!!!!!!!!!!
!!!!!!!!!!!!!
!!!!!!!!!!!!!


    ! SiII:
    
    ! collisions with e-
    
    trates(37) = 1.2e-6 * T2**(-0.5)

    ! collisions with H
    
    trates(38) = 4.95e-10 * T2**(0.24)

    
!!!!!!!!!!!!!
!!!!!!!!!!!!!
!!!!!!!!!!!!!


    ! CI:
    
    ! collisions with H2
    
    ! 1->0
    ortho = 8.7e-11 - 6.6e-11 * exp(-temp / 218.3) + 6.6e-11 * exp(-2.0 * temp / 218.3) * fortho
    para = 7.9e-11 - 8.7e-11 * exp(-temp / 126.4) + 1.3e-10 * exp(-2.0 * temp / 126.4) * fpara
    trates(39) = ortho + para
    
    ! 2->0
    ortho = 1.2e-10 - 6.1e-11 * exp(-temp / 387.3) * fortho
    para = 1.1e-10 - 8.6e-11 * exp(-temp / 223.0) + 8.7e-11 * exp(-2.0 * temp / 223.0) * fpara
    trates(40) = ortho + para

    ! 2->1
    ortho = 2.9e-10 - 1.9e-10 * exp(-temp / 348.9) * fortho
    para = 2.7e-10 - 2.6e-10 * exp(-temp / 250.7) + 1.8e-10 * exp(-2.0 * temp / 250.7) * fpara
    trates(41) = ortho + para
    

    ! collisions with H

    ! 1->0
    trates(42) = 1.6e-10 * T2**0.14
    
    ! 2->0
    trates(43) = 9.2e-11 * T2**0.26

    ! 2->1
    trates(44) = 2.9e-10 * T2**0.26
    
    
    ! collisions with e-

    ! 1->0
    if(temp .le. 1000.0) then
       trates(45) = 2.88e-6 * temp**(-0.5) * exp(-9.25141 - 7.73782e-1 * ln_temp + & 
            3.61184e-1 * ln_temp**2 - 1.50892e-2 * ln_temp**3 - 6.56325e-4 * ln_temp**4)
    else
       trates(45) = 2.88e-6 * temp**(-0.5) * exp(-4.446e-2 - 2.27913e2 * ln_temp + & 
            4.2595e1 * ln_temp**2 - 3.47620 * ln_temp**3 + 1.0508e-1 * ln_temp**4)
    end if
    
    ! 2->0
    if(temp .le. 1000.0) then
       trates(46) = 1.73e-6 * temp**(-0.5) * exp(-7.69735 - 1.30743 * ln_temp + & 
            0.697638 * ln_temp**2 - 0.111338 * ln_temp**3 + 0.705277e-2 * ln_temp**4)
    else
       trates(46) = 1.73e-6 * temp**(-0.5) * exp(3.50609e2 - 1.87474e2 * ln_temp + & 
            3.61803e1 * ln_temp**2 - 3.03283 * ln_temp**3 + 9.38138e-2 * ln_temp**4)
    end if
    
    
    ! 2->1
    
    if(temp .le. 1000.0) then
       trates(47) = 1.73e-6 * temp**(-0.5) * exp(-7.4387 - 0.57443 * ln_temp + & 
            0.358264 * ln_temp**2 - 4.18166e-2 * ln_temp**3 + 2.35272e-3 * ln_temp**4)
    else
       trates(47) = 1.73e-6 * temp**(-0.5) * exp(3.86186e2 - 2.02193e2 * ln_temp + & 
            3.85049e1 * ln_temp**2 - 3.19268 * ln_temp**3 + 9.78573e-2 * ln_temp**4)
    end if


    ! collisions with H+

    ! 1->0
    if(temp .le. 5000.0) then
       trates(48) = temp**0.45 * (9.6e-11 - 1.8e-14 * temp + 1.9e-18 * temp**2)
    else
       trates(48) = 8.9e-10 * temp**0.117
    end if
    
    ! 2->0
    if(temp .le. 5000.0) then
       trates(49) = temp * (3.1e-12 - 6.0e-16 * temp + 3.9e-20 * temp**2)
    else
       trates(49) = 2.3e-9 * temp**0.0965
    end if
 

    ! 2->1
    if(temp .le. 5000.0) then
       trates(50) = temp**0.7 * (1.0e-10 - 2.2e-14 * temp + 1.7e-18 * temp**2)
    else
       trates(50) = 9.2e-9 * temp**0.0535
    end if


    !!!!!!!!!!!!!

    ! More accurate low density H2 cooling rate, Glover and Abel (2008):
    ! Assumes 3:1 ortho:para ratio
    
    ! Collisions with H
    ! Assume above 10^4 K, rate remains constant
    
    if(temp .ge. 10.0 .and. temp .le. 100.0) then
       trates(51) = -16.818342 + 37.383713 * logtemp3 + & 
            58.145166 * logtemp3**2 + 48.656103 * logtemp3**3 + & 
            20.159831 * logtemp3**4 + 3.8479610 * logtemp3**5
       trates(51) = 10.0**(trates(51))
    elseif (temp .le. 1000.0) then
       trates(51) = -24.311209 + 3.5692468 * logtemp3 + & 
            -11.332860 * logtemp3**2 - 27.850082 * logtemp3**3 & 
            -21.328264 * logtemp3**4 - 4.2519023 * logtemp3**5
       trates(51) = 10.0**(trates(51))
    elseif (temp .le. 1.0e4) then
       trates(51) = -24.311209 + 4.6450521 * logtemp3 + & 
          -3.7209846 * logtemp3**2 + 5.9369081 * logtemp3**3 & 
          -5.5108047 * logtemp3**4 + 1.5538288 * logtemp3**5
       trates(51) = 10.0**(trates(51))
    else
       trates(51) = -24.311209 + 4.6450521 + & 
            -3.7209846 + 5.9369081 + & 
            -5.5108047 + 1.5538288
       trates(51) = 10.0**(trates(51))
    end if
    

    ! collisions with H2
    ! Assume constant rate above 10^4 K
    
    if(temp .gt. 100.0 .and. temp .le. 1.0e4) then
       trates(52) = -23.962112 + 2.09433740 * logtemp3 + & 
            -0.77151436 * logtemp3**2 + 0.43693353 * logtemp3**3 & 
            -0.14913216 * logtemp3**4 - 0.033638326 * logtemp3**5
       trates(52) = 10.0**(trates(52))
    else
       trates(52) = -23.962112 + 2.09433740 + & 
            -0.77151436 + 0.43693353 + & 
            -0.14913216 - 0.033638326
       trates(52) = 10.0**(trates(52))
    end if


    ! collisions with He
    ! Assume constant rate above 10^4 K
    if(temp .gt. 10.0 .and. temp .le. 1.0e4) then
       trates(53) = -23.689237 + 2.1892372 * logtemp3 + & 
            -0.81520438 * logtemp3**2 + 0.29036281 * logtemp3**3 & 
            -0.16596184 * logtemp3**4 + 0.19191375 * logtemp3**5
       trates(53) = 10.0**(trates(53))
    else
       trates(53) = -23.689237 + 2.1892372 & 
            -0.81520438 + 0.29036281 & 
            -0.16596184 + 0.19191375
       trates(53) = 10.0**(trates(53))
    end if


    ! collisions with H+
    ! Assume constant rate above 10^4 K
    if(temp .gt. 10.0 .and. temp .le. 10000.0) then
       trates(54) = -21.716699 + 1.3865783 * logtemp3 + & 
            -0.37915285 * logtemp3**2 + 0.11453688 * logtemp3**3 & 
            -0.23214154 * logtemp3**4 + 0.058538864 * logtemp3**5
       trates(54) = 10.0**(trates(54))
    else
       trates(54) = -21.716699 + 1.3865783 + & 
            -0.37915285 + 0.11453688 + & 
            -0.23214154 + 0.058538864
       trates(54) = 10.0**(trates(54))
    end if
    
    ! collisions with e-
    if(temp .gt. 10.0 .and. temp .le. 200.0) then
       trates(55) = -34.286155 - 48.537163 * logtemp3 + & 
            -77.121176 * logtemp3**2 + -51.352459 * logtemp3**3 & 
            -15.169160 * logtemp3**4 - 0.98120322 * logtemp3**5
       trates(55) = 10.0**(trates(55))
    elseif (temp .le. 1.0e4) then
       trates(55) = -22.190316 + 1.5728955 * logtemp3 + & 
            -0.21335100 * logtemp3**2 + 0.96149759 * logtemp3**3 & 
            -0.91023195 * logtemp3**4 + 0.13749749 * logtemp3**5
       trates(55) = 10.0**(trates(55))
    else
       trates(55) = -22.190316 + 1.5728955 + & 
            -0.21335100 + 0.96149759 + & 
            -0.91023195 + 0.13749749
       trates(55) = 10.0**(trates(55))
    end if

    ! Collisions with H+, J = 0 -> 1 and vise versa

    trates(56) = 4.76e-24 * (9.0 * exp(-170.5 / temp) * 0.25 - 0.75)


  
  return
end subroutine thermal_rates



subroutine thermal(y, temp, tdust, rpar, cooling, heating)

  use CO_CoolGlobal
  
  implicit none

  real, intent(INOUT) :: y(NTOT), temp, tdust, rpar(nrpar)
  real, intent(OUT) :: cooling, heating
  
  real, parameter :: ev_in_erg = 1.602e-12
  real :: nh, G0, ne, NHtot

  real :: lambda_fs, lambda_lya, lambda_gd
  real :: gamma_grain_pe, gamma_cr, gamma_h2form
  real :: beta, G0_eff, fac, rho, ye, eff, cooling_cmb
  real :: qion, qion_HI, qion_H2, lognh, nhsq, ncr, lshield
  
  integer :: i, npoint
  real :: t_min, t_max
  

  nh = rpar(1)
  NHtot = rpar(2)
  G0 = rpar(5)
  ne = y(ie) * nh
  lognh = log10(nh)
  nhsq = nh * nh
  
  if(density_gradient_shielding_length) then
     rho = nh * mh * (1.0 + 4.0*YHELIUM)
     lshield = rho / rpar(7)
     NHtot = nh * lshield
     rpar(2) = NHtot
  endif

  
  
  if(simple_cooling) then
     
     !! Kim, Ostriker, & Kim 2013:
     lambda_lya = 2.0e-19 * exp(-1.184e5 / (temp + 1000.0)) * nhsq
     lambda_fs = 2.8e-28 * sqrt(temp) * exp(-92.0 / temp) * nhsq
     gamma_grain_pe = 2.0e-26 * (G0 + 0.0024) * nh

     cooling = lambda_fs + lambda_lya
     heating = gamma_grain_pe
     
     return
  endif
  

  call cooling_processes(y,temp,rpar,cooling)
  if(temp_floor) then
     call cooling_processes(y,TCMB,rpar,cooling_cmb)
     cooling = cooling - cooling_cmb
  endif
  
  !! Heating processes:

  !! include dust-gas energy transfer, since it should not be included as CMB cooling
  !! lambda_gd = positive for gas cooling (Tdust  < Tgas)
  !!           = negative for gas heating (Tdust > Tgas)
  lambda_gd = alpha_gd * sqrt(temp) * (temp - Tdust) * nhsq
  cooling = cooling + lambda_gd
  
  !! Grain photoelectric heating.
  !! See Omukai et al. 2008 and Krumholz 2014:

  G0_eff = G0 * exp(-0.5 * NHtot * 1.0e-21)
  beta = 0.735 / temp**0.068
  fac = G0_eff * sqrt(temp) / ne
  eff = 4.9e-2 / (1.0 + 4.0e-3 * (fac)**0.73) & 
      + 3.7e-2 * (temp / 1.0e4)**0.7 / (1.0 + 2.0e-4 * fac)
  gamma_grain_pe = 1.0e-24 * eff * G0_eff * nh & 
               - 4.65e-30 * temp**0.94 * (fac)**beta * ne * nh

  
  !! cosmic ray ionization heating rate:
  !! see Krumholz 2014
  
  qion_HI = 1.04e-11 + 4.23e-11 * sqrt(y(ie) / (y(ie) + 0.07))
 
  if(lognh .le. 2.0) then
     qion_H2 = 10.0 * ev_in_erg
  else if(lognh .gt. 2.0 .and. lognh .lt. 4.0) then
     qion_H2 = (10.0 + 3.0 * (lognh - 2.0) / 2.0) * ev_in_erg
  else if(lognh .ge. 4.0 .and. lognh .lt. 7.0) then
     qion_H2 = (13.0 + 4.0 * (lognh - 4.0) / 3.0) * ev_in_erg
  else if(lognh .ge. 7.0 .and. lognh .lt. 10.0) then
     qion_H2 = (17.0 + (lognh - 7.0) / 3.0) * ev_in_erg
  else
     qion_H2 = 18.0 * ev_in_erg
  endif

  qion = y(iH) * qion_HI + 2.0 * y(iH2) * qion_H2
  gamma_cr = zeta * qion * nh


  !! H2 formation heating on dust grains:
  !! see Hollenbach & Mckee 1979
  !! do not include h2 formation heating in this pure radiative cooling formulation
  
  !ncr = 1.0e6 / sqrt(temp) / (1.6 * y(iH) * exp( -(400.0 / temp)**2) + 1.4 * y(iH2) * exp(-(1.2e4 / (temp + 1.2e3))))
  !gamma_h2form = ev_in_erg * kh2 * nhsq * y(iH) * (0.2 + 4.2 / (1.0 + ncr / nh))
  gamma_h2form = 0.0
  
  
  heating = gamma_grain_pe + gamma_cr + gamma_h2form 
  
 
  return
end subroutine thermal




subroutine cooling_processes(y,temp,rpar,cooling)
  
  use CO_CoolGlobal
  implicit none
  
  real, intent(IN) :: y(NTOT), temp, rpar(nrpar)
  real, intent(OUT) :: cooling


  real :: lambda_CO, lambda_CO_neu, eff
  real :: lambda_fs, lambda_lya, lambda_gd, lambda_ff
  real :: trates(NTRATES), nhsq, ye

  nhsq = rpar(1) * rpar(1)
  
  !call thermal_rates_interp(temp,trates)
  call thermal_rates(temp,trates)
  
  !! tabulated Despotic CO cooling
  ! call CO_cooling(y, temp, lambda_CO, rpar)
  !! CO Cooling function from Neufeld et al. 1992, 1995
  call CO_cooling_neufeld(y,temp,lambda_CO,rpar)
  

  !! fine structure line cooling for CI, CII, and OI:  
  call fine_structure_cooling(y, temp, lambda_fs, trates, rpar)
  
  !! ly-alpha cooling (HI excitation)
  !lambda_lya = trates(9) * y(iH) * y(ie) * nhsq
  ! here we will use the exact expression from Kim and Ostriker
  lambda_lya = 2.0e-19 * exp(-1.184e5 / (temp + 1000.0)) * nhsq
 
  
  !! free-free cooling:
  lambda_ff = 0.0
  if(temp .gt. 1.0e5) then
     ! only include above 10^5 K when we can assume gas is fully ionized
     ye = 1.08
     lambda_ff = trates(1) * ye**2 * nhsq 
  endif
 
  cooling = lambda_CO + lambda_fs + lambda_lya + lambda_gd + lambda_ff
  

  return
end subroutine cooling_processes





subroutine fine_structure_cooling(y, temp, lambda_fs, trates, rpar)

  use CO_CoolGlobal
  
  implicit none
  
  real, intent(IN) :: y(NTOT), temp, rpar(nrpar), trates(NTRATES)
  real, intent(OUT) :: lambda_fs
  
  !! CII parameters
  real, parameter :: A10_CII = 2.3e-6
  real, parameter :: E10_CII = 1.27e-14
  real, parameter :: g1_CII = 4.0
  real, parameter :: g0_CII = 2.0
  
  !! OI parameters
  real, parameter :: A10_OI = 8.9e-5
  real, parameter :: A20_OI = 1.3e-10
  real, parameter :: A21_OI = 1.8e-5
  real, parameter :: E10_OI = 3.17e-14
  real, parameter :: E20_OI = 4.55e-14
  real, parameter :: E21_OI = 1.35e-14
  real, parameter :: g0_OI = 5.0
  real, parameter :: g1_OI = 3.0
  real, parameter :: g2_OI = 1.0
  
  !! CI parameters
  real, parameter :: A10_CI = 7.9e-8
  real, parameter :: A20_CI = 2.1e-14
  real, parameter :: A21_CI = 2.7e-7
  real, parameter :: E10_CI = 3.31e-15
  real, parameter :: E20_CI = 8.69e-15
  real, parameter :: E21_CI = 5.38e-15
  real, parameter :: g0_CI = 1.0
  real, parameter :: g1_CI = 3.0
  real, parameter :: g2_CI = 5.0

  real :: nh, f0, f1, f2, w
  real :: beta, ngamma
  real :: C10, C21, C20
  real :: lambda_Cplus, Lambda_C, lambda_O
  
  real :: ng20,ng21,ng10,t20,t21,t10
  real :: b20,b21,b10
  
  
  real :: f0_old, f1_old, f2_old
  integer, parameter :: maxiter = 50
  real, parameter :: tolerance = 1.0e-4
  logical :: converged
  integer :: i
  
  
  if(temp  .le. 0.0) then
     print*, "negative/zero temp in fs cooling. not ok."
     print*, temp
     stop
  endif
  
  
  nh = rpar(1)

  !! C+, two level system:
  
  lambda_Cplus = 0.0
  if(y(iCplus) .gt. 0.0) then
     beta = 1.0 
     ng10 = 1.0 / (exp(E10_CII / kb / TCMB) - 1.0)
     
     C10 = (  y(ie)  * trates(22)  & 
          + y(iH)  * trates(23) & 
          + y(iH2) * trates(24)) * nh
     
     call two_level_system(temp,A10_CII,C10,beta,ng10,g1_CII/g0_CII,E10_CII,f1,f0) 
     
     if(fine_structure_line_opacity) then
     do i=1,maxiter
        
        f1_old = f1
        f0_old = f0
        
        call calculate_beta(y(iCplus),12.0,temp,f1,f0,E10_CII,A10_CII,g1_CII/g0_CII,rpar,beta)
        call two_level_system(temp,A10_CII,C10,beta,ng10,g1_CII/g0_CII,E10_CII,f1,f0) 
        
        ! agressively force convergence:
        w = 1.0 - real(i) / real(maxiter)
        f0 = (1.0 - w) * f0_old + w * f0
        f1 = (1.0 - w) * f1_old + w * f1
        
        ! check for convergence
        converged = .false.
        converged = (abs(f1_old - f1) / (f1 * tolerance) .lt. 1.0) & 
             .and.  (abs(f0_old - f0) / (f0 * tolerance) .lt. 1.0)
        
        if(converged) exit
        
        
     enddo
     
     if(.not. converged) then
        print*, "C+ did not converge! wtf", temp
        print*, y
        stop
     endif
     endif
     
     lambda_Cplus = beta * ((1.0 + ng10) * f1 - (g1_CII / g0_CII) * ng10 * f0) * A10_CII * E10_CII * y(iCplus)
  endif
  

  ! neutral C - three level system:

  C10 = (y(ie) * trates(45) + & 
       y(iH) * trates(42) + & 
       y(iHplus) * trates(48) + & 
       y(iH2) * trates(39)) * nh
  C20 = (y(ie) * trates(46) + & 
       y(iH) * trates(43) + & 
       y(iHplus) * trates(49) + & 
       y(iH2) * trates(40)) * nh
  C21 = (y(ie) * trates(47) + & 
       y(iH) * trates(44) + & 
       y(iHplus) * trates(50) + & 
       y(iH2) * trates(41)) * nh
  
  b20 = 1.0
  b21 = 1.0
  b10 = 1.0

  ng20 = 1.0 / (exp(E20_CI / kb / TCMB) - 1.0)
  ng21 = 1.0 / (exp(E21_CI / kb / TCMB) - 1.0)
  ng10 = 1.0 / (exp(E10_CI / kb / TCMB) - 1.0)
  
  t20 = exp(-E20_CI / (kb * temp))
  t21 = exp(-E21_CI / (kb * temp))
  t10 = exp(-E10_CI / (kb * temp))
  

  call three_level_system(A20_CI,A21_CI,A10_CI,C20,C21,C10,b20,b21,b10,&
       ng20,ng21,ng10,g2_CI/g0_CI,g2_CI/g1_CI,g1_CI/g0_CI,t20,t21,t10,f0,f1,f2)


  if(fine_structure_line_opacity) then
  do i=1,maxiter
     
     f2_old = f2
     f1_old = f1
     f0_old = f0
     
     call calculate_beta(y(iC),12.0,temp,f2,f0,E20_CI,A20_CI,g2_CI/g0_CI,rpar,b20)
     call calculate_beta(y(iC),12.0,temp,f2,f1,E21_CI,A21_CI,g2_CI/g1_CI,rpar,b21)
     call calculate_beta(y(iC),12.0,temp,f1,f0,E10_CI,A10_CI,g1_CI/g0_CI,rpar,b10)
     call three_level_system(A20_CI,A21_CI,A10_CI,C20,C21,C10,b20,b21,b10,&
          ng20,ng21,ng10,g2_CI/g0_CI,g2_CI/g1_CI,g1_CI/g0_CI,t20,t21,t10,f0,f1,f2)

          
     ! agressively force convergence:
     w = 1.0 - real(i) / real(maxiter)
     f0 = (1.0 - w) * f0_old + w * f0
     f1 = (1.0 - w) * f1_old + w * f1
     f2 = (1.0 - w) * f2_old + w * f2

     ! check for convergence
     converged = .false.
     converged = (abs(f1_old - f1) / (f1 * tolerance) .lt. 1.0) & 
          .and.  (abs(f0_old - f0) / (f0 * tolerance) .lt. 1.0) & 
          .and.  (abs(f2_old - f2) / (f2 * tolerance) .lt. 1.0) 
     
     if(converged) exit
     
     !print*, i
     
  enddo

  if(.not. converged) then
     print*, "CI did not converge! wtf"
     stop
  endif
  endif

  lambda_C = (b20 * ((1.0 + ng20) * f2 - (g2_CI / g0_CI) * ng20 * f0) * A20_CI * E20_CI & 
            + b21 * ((1.0 + ng21) * f2 - (g2_CI / g1_CI) * ng21 * f1) * A21_CI * E21_CI  &
            + b10 * ((1.0 + ng10) * f1 - (g1_CI / g0_CI) * ng10 * f0) * A10_CI * E10_CI) * y(iC)
    




  ! neutral O - three level system:
  
  
  C10 = (y(ie) * trates(31) + & 
       y(iH) * trates(28) + & 
       y(iHplus) * trates(34) + & 
       y(iH2) * trates(25)) * nh
  C20 = (y(ie) * trates(32) + & 
       y(iH) * trates(29) + & 
       y(iHplus) * trates(35) + & 
       y(iH2) * trates(26)) * nh
  C21 = (y(ie) * trates(33) + & 
       y(iH) * trates(30) + & 
       y(iHplus) * trates(36) + & 
       y(iH2) * trates(27)) * nh
  
  b20 = 1.0
  b21 = 1.0
  b10 = 1.0

  ng20 = 1.0 / (exp(E20_OI / kb / TCMB) - 1.0)
  ng21 = 1.0 / (exp(E21_OI / kb / TCMB) - 1.0)
  ng10 = 1.0 / (exp(E10_OI / kb / TCMB) - 1.0)
  
  t20 = exp(-E20_OI / (kb * temp))
  t21 = exp(-E21_OI / (kb * temp))
  t10 = exp(-E10_OI / (kb * temp))
  
  call three_level_system(A20_OI,A21_OI,A10_OI,C20,C21,C10,b20,b21,b10,&
       ng20,ng21,ng10,g2_OI/g0_OI,g2_OI/g1_OI,g1_OI/g0_OI,t20,t21,t10,f0,f1,f2)

  if(fine_structure_line_opacity) then
  do i=1,maxiter
     
     f2_old = f2
     f1_old = f1
     f0_old = f0
     
     call calculate_beta(y(iO),16.0,temp,f2,f0,E20_OI,A20_OI,g2_OI/g0_OI,rpar,b20)
     call calculate_beta(y(iO),16.0,temp,f2,f1,E21_OI,A21_OI,g2_OI/g1_OI,rpar,b21)
     call calculate_beta(y(iO),16.0,temp,f1,f0,E10_OI,A10_OI,g1_OI/g0_OI,rpar,b10)
     call three_level_system(A20_OI,A21_OI,A10_OI,C20,C21,C10,b20,b21,b10,&
          ng20,ng21,ng10,g2_OI/g0_OI,g2_OI/g1_OI,g1_OI/g0_OI,t20,t21,t10,f0,f1,f2)

          
     ! agressively force convergence:
     w = 1.0 - real(i) / real(maxiter)
     f0 = (1.0 - w) * f0_old + w * f0
     f1 = (1.0 - w) * f1_old + w * f1
     f2 = (1.0 - w) * f2_old + w * f2

     ! check for convergence
     converged = .false.
     converged = (abs(f1_old - f1) / (f1 * tolerance) .lt. 1.0) & 
          .and.  (abs(f0_old - f0) / (f0 * tolerance) .lt. 1.0) & 
          .and.  (abs(f2_old - f2) / (f2 * tolerance) .lt. 1.0) 
     
     if(converged) exit
     
     !print*, i
     
  enddo

  if(.not. converged) then
     print*, "OI did not converge! wtf"
     stop
  endif
endif

  lambda_O = (b20 * ((1.0 + ng20) * f2 - (g2_OI / g0_OI) * ng20 * f0) * A20_OI * E20_OI  & 
           + b21 * ((1.0 + ng21) * f2 - (g2_OI / g1_OI) * ng21 * f1) * A21_OI * E21_OI  &
           + b10 * ((1.0 + ng10) * f1 - (g1_OI / g0_OI) * ng10 * f0) * A10_OI * E10_OI) * y(iO)
 
  
  lambda_fs = (lambda_Cplus + lambda_C + lambda_O) * nh
  
 
 return
end subroutine fine_structure_cooling



subroutine calculate_beta(y,mu,temp,f1,f0,e,a,gfac,rpar,beta)

  use CO_CoolGlobal
  
  implicit none
  real, intent(IN) :: y, mu, temp, f1, f0, e, a, gfac, rpar(nrpar)
  real, intent(OUT) :: beta

  real, parameter :: hc = 1.99e-16
  
  real :: tauline, NHtot, cs, sigma_tot, lambda
  
  NHtot = rpar(2)
  
  if(NHtot .le. 0.0) then
     beta = 0.0
     return
  endif
  
  lambda = hc / e
  
  ! assume no non-thermal component
  cs = sqrt(kb * temp / (mh * mu))
  sigma_tot = cs
  
  tauline = gfac * a * lambda**3 / sigma_tot  * y * NHtot * f0 * (1.0 - f1 / f0 / gfac) / 63.0
  
  beta = (1.0 - exp(-3.0 * tauline)) / (3.0 * tauline)
  
  return
end subroutine calculate_beta





subroutine two_level_system(t,a,c,b,ng,gfac,e,f1,f0)

  use CO_CoolGlobal

  implicit none

  real, intent(IN) :: t,a,c,b,gfac,e,ng
  real, intent(OUT) :: f1,f0
  
  real :: alpha, gamma
  
  alpha = a * b * (1.0 + ng) + c
  gamma = gfac * (c * exp(-e / (kb * t)) + a * ng * b)
  
  f1 = gamma / (alpha + gamma)
  f0 = 1.0 - f1
  
  if(f0 .lt. 0.0) then
     print*, "problem in two level system. negative level pop", f1, f0
     stop
  endif
     
  return
end subroutine two_level_system


! use cramers rule to solve three level systems
subroutine three_level_system(a20,a21,a10,c20,c21,c10,b20,b21,b10, & 
                    ng20,ng21,ng10,g20,g21,g10,t20,t21,t10,f0,f1,f2)



  real, intent(IN) :: a20,a21,a10,c20,c21,c10,b20,b21,b10
  real, intent(IN) :: ng20,ng21,ng10,g20,g21,g10,t20,t21,t10
  real, intent(OUT) :: f0,f1,f2
  
  real :: q1, q2, q3, q4, q5, q6
  real :: d, dx, dy, dz
  
  
  q1 = c21 + a20 * b20 * (1.0 + ng20) + a21 * b21 * (1.0 + ng21)
  q2 = c21 * g21 * t21 + g21 * a21 * ng21 * b21
  q3 = c20 * g20 * t20 + g20 * a20 * ng20 * b20
  q4 = c10 + c21 * g21 * t21 + b10 * a10 * (1.0 + ng10) + g21 * a21 * ng21 * b21
  q5 = b21 * a21 * c21 + a21 * ng21 * b21
  q6 = g10 * a10 * ng10 * b10 + g10 * c10 * t10
  
  d =   q3 * (-q4 - q5) & 
      - q2 * (q6 - q5) & 
      - q1 * (q6 + q4)
  
  dx = - q2 * (-q5) & 
       - q1 * (q4)
  
  dy = q3 * (-q5) & 
       -q1 * (q6)
  
  dz = -q3 * q4 & 
       -q2 * q6
  
  f0 = dx/d
  f1 = dy/d
  f2 = dz/d

  
  if(f0 .lt. 0.0 .or. f1 .lt. 0.0 .or. f2 .lt. 0.0) then
     print*, "problem in three level system. negative level pop", f0, f1, f2, a20
     print*, "dx, dy, dz, d", dx, dy, dz, d
     print*, "the q's:", q1, q2, q3, q4, q5, q6
     stop
  endif
  
  
  return
end subroutine three_level_system


!! CO cooling function from neufeld
!! only used for comparative purposes
subroutine CO_cooling_neufeld(y, temperature, lambda, rpar)

  use CO_CoolGlobal

  implicit none


  real, intent(IN) :: y(NTOT), temperature, rpar(nrpar)
  real, intent(OUT) :: lambda
  
  integer, parameter :: ntemp = 11
  integer, parameter :: nnco = 10
  
   
  real, dimension(ntemp) :: logtempv = (/ 1.000, 1.301, 1.477, 1.699, 1.903, 2.0, &
                                          2.477, 2.778, 3.000, 3.176, 3.301/)
  
  real, dimension(nnco) :: logncov = (/ 14.5, 15.0, 15.5, 16.0, 16.5, &
       17.0, 17.5, 18.0, 18.5, 19.0 /)
   
  real, dimension(ntemp) :: logL0v = (/ 24.96, 24.56, 24.36, 24.12, 23.88, 23.80, & 
                                    23.40, 23.07, 22.81, 22.61, 22.47  /)
  
  real, dimension(ntemp,nnco) :: logLLTEv = RESHAPE([21.15, 20.39, 19.97, 19.47, 19.02, 18.81, 17.82, 17.23, 16.82, 16.52, 16.33, &
                                               21.25, 20.45, 20.02, 19.50, 19.04, 18.83, 17.83, 17.24, 16.82, 16.52, 16.33, & 
                                               21.45, 20.58, 20.12, 19.57, 19.09, 18.87, 17.84, 17.24, 16.83, 16.52, 16.33, &
                                               21.73, 20.81, 20.31, 19.71, 19.20, 18.96, 17.88, 17.26, 16.84, 16.53, 16.33, &
                                               22.08, 21.11, 20.58, 19.95, 19.39, 19.14, 17.98, 17.32, 16.88, 16.56, 16.35, &
                                               22.47, 21.46, 20.91, 20.25, 19.67, 19.40, 18.15, 17.44, 16.96, 16.62, 16.40, &
                                               22.89, 21.85, 21.28, 20.60, 19.99, 19.71, 18.41, 17.65, 17.13, 16.76, 16.52, &
                                               23.32, 22.25, 21.67, 20.97, 20.35, 20.07, 18.72, 17.93, 17.38, 16.98, 16.72, &
                                               23.76, 22.68, 22.08, 21.37, 20.74, 20.45, 19.07, 18.26, 17.69, 17.27, 17.00, &
                                               24.21, 23.11, 22.51, 21.78, 21.14, 20.85, 19.45, 18.62, 18.03, 17.61, 17.33 ], &
                                               [ntemp,nnco])

  real, dimension(ntemp,nnco) :: lognhalfv = RESHAPE([3.31, 3.72, 3.96, 4.24, 4.39, 4.51, 4.95, 5.24, 5.43, 5.57, 5.68, &
                                                3.11, 3.59, 3.87, 4.17, 4.35, 4.47, 4.93, 5.23, 5.42, 5.57, 5.68, &
                                                2.80, 3.34, 3.64, 4.01, 4.23, 4.36, 4.88, 5.20, 5.42, 5.56, 5.68, &
                                                2.37, 2.95, 3.31, 3.73, 3.99, 4.13, 4.75, 5.13, 5.36, 5.54, 5.65, &
                                                1.90, 2.50, 2.87, 3.33, 3.61, 3.78, 4.50, 4.95, 5.24, 5.44, 5.58, &
                                                1.40, 2.02, 2.40, 2.86, 3.17, 3.34, 4.12, 4.65, 5.00, 5.25, 5.41, &
                                                0.90, 1.52, 1.90, 2.38, 2.69, 2.87, 3.67, 4.25, 4.63, 4.92, 5.13, &
                                                0.40, 1.02, 1.41, 1.89, 2.19, 2.38, 3.20, 3.77, 4.18, 4.15, 4.73, &
                                                -0.09, 0.52, 0.91, 1.39, 1.69, 1.88, 2.69, 3.29, 3.70, 4.03, 4.26, &
                                                -0.59, 0.03, 0.41, 0.88, 1.19, 1.38, 2.20, 2.79, 3.20, 3.54, 3.78], &
                                                [ntemp,nnco])

  
  real, dimension(ntemp,nnco) :: alphav = RESHAPE( [0.37, 0.41, 0.42, 0.41, 0.38, 0.37, 0.38, 0.38, 0.37, 0.36, 0.35, &
                                             0.38, 0.42, 0.44, 0.42, 0.39, 0.37, 0.37, 0.37, 0.37, 0.36, 0.35, &
                                             0.38, 0.43, 0.44, 0.43, 0.40, 0.39, 0.36, 0.36, 0.36, 0.35, 0.34, &
                                             0.36, 0.41, 0.43, 0.43, 0.40, 0.39, 0.35, 0.35, 0.35, 0.34, 0.33, &
                                             0.34, 0.37, 0.39, 0.39, 0.39, 0.39, 0.35, 0.34, 0.34, 0.33, 0.32, &
                                             0.31, 0.33, 0.35, 0.35, 0.35, 0.34, 0.37, 0.35, 0.33, 0.32, 0.31, &
                                             0.28, 0.31, 0.32, 0.32, 0.31, 0.31, 0.40, 0.37, 0.35, 0.32, 0.31, &
                                             0.26, 0.28, 0.29, 0.28, 0.28, 0.28, 0.43, 0.40, 0.38, 0.35, 0.33, &
                                             0.25, 0.26, 0.26, 0.26, 0.26, 0.26, 0.46, 0.43, 0.41, 0.38, 0.36, &
                                             0.23, 0.25, 0.25, 0.24, 0.25, 0.25, 0.48, 0.46, 0.44, 0.41, 0.39],&
                                             [ntemp,nnco])
  
  
  real :: LLTE, L0, nH2, nh, NCO, sigh2, ve
  real :: nhalf, alpha, tff, log_nco, LCO
  real :: log_temp, n_co, dvdz, dvdr
  real :: t, u, y1, y2, y3, y4
  integer :: temp_index, nco_index, i
  real :: dx, slope
  
  nh = rpar(1)
 
!! effective nH2 - see Glover et al. 2010
  sigh2 = 3.3e-16 / (temperature / 1000.0)**0.25
  ve = 1.03e4 * sqrt(temperature)
  nH2 = (y(iH2) + 1.414 * (2.3e-15 / sigh2) + (1.3e-8 / sigh2 / ve) * y(ie)) * nh

  log_temp = log10(temperature)
  n_co = y(iCO) * nh
  
  if(CO_LVG) then
     dvdr = rpar(8)
     NCO = y(iCO) * nh / dvdr 
  else
     NCO = rpar(2) * y(iCO)
  endif

  dvdr = rpar(8)
  NCO = y(iCO) * nh / dvdr
  ! convert from cm^-3 s^-1 to cm^-2 km^-1 s^-1
  NCO = NCO * 1.0e5
  log_nco = log10(NCO)
  
  call locate(logtempv, ntemp, log_temp, temp_index) 
  if(temp_index .le. 0) then
     temp_index = 1
     !log_temp = logtempv(temp_index)
  endif
  if(temp_index .ge. ntemp) then
     temp_index = ntemp-1
     !log_temp = logtempv(ntemp)
  endif
  call locate(logncov, nnco, log_nco, nco_index)
  if(nco_index .le. 0) then
     nco_index = 1
     log_nco = logncov(nco_index)
  endif
  if(nco_index .ge. nnco) then
     nco_index = nnco-1
     log_nco = logncov(nnco)
  endif
  t = (log_temp - logtempv(temp_index)) / ( logtempv(temp_index+1) - logtempv(temp_index))
  u = (log_nco - logncov(nco_index)) / ( logncov(nco_index+1) - logncov(nco_index))


  y1 = logLLTEv(temp_index, nco_index)
  y2 = logLLTEv(temp_index+1,nco_index)
  y3 = logLLTEv(temp_index+1,nco_index+1)
  y4 = logLLTEv(temp_index,nco_index+1)
  LLTE = 10.0**(-((1.0-t)*(1.0-u)*y1 + t*(1.0-u)*y2  + t*u*y3 + (1.0-t)*u*y4))


  y1 = alphav(temp_index, nco_index)
  y2 = alphav(temp_index+1,nco_index)
  y3 = alphav(temp_index+1,nco_index+1)
  y4 = alphav(temp_index,nco_index+1)
  alpha = (1.0-t)*(1.0-u)*y1 + t*(1.0-u)*y2 + t*u*y3 + (1.0-t)*u*y4

  y1 = lognhalfv(temp_index, nco_index)
  y2 = lognhalfv(temp_index+1,nco_index)
  y3 = lognhalfv(temp_index+1,nco_index+1)
  y4 = lognhalfv(temp_index,nco_index+1)
  nhalf = 10.0**((1.0-t)*(1.0-u)*y1 + t*(1.0-u)*y2  + t*u*y3 + (1.0-t)*u*y4)

  slope = (logL0v(temp_index+1) - logL0v(temp_index)) / ( logtempv(temp_index+1) - logtempv(temp_index))
  dx = log_temp - logtempv(temp_index)
  L0 = 10.0**(-(logL0v(temp_index) + slope * dx))

  
  LCO = 1.0 / (1.0 / L0 + nH2 / LLTE + (nH2 / nhalf)**alpha * (1.0 - nhalf*L0/LLTE) / L0)
  lambda =  nH2 * y(iCO) * nh * LCO
  
  return
end subroutine CO_cooling_neufeld



subroutine CO_cooling(y, temperature, lambda, rpar)

  use CO_CoolGlobal

  implicit none



  real, intent(IN) :: y(NTOT), temperature, rpar(nrpar)
  real, intent(OUT) :: lambda

  character(len=128) :: path, dummys
  real :: dummyr, val, val1, val2, val3, val4

  integer, save :: ntemp, nnh2, nnco
  real, save, allocatable :: log_t_table(:)
  real, save, allocatable :: log_nh2_table(:)
  real, save, allocatable :: log_nco_table(:)
  real, save, allocatable :: log_lambda_table(:,:,:)
  real, save :: logtemp_min, logtemp_max
  real, save :: lognh2_min, lognh2_max
  real, save :: lognco_min, lognco_max
  
  integer :: i, j, k
  integer :: ktemp, knh2, knco
  
  real :: NCO, nH2, nh
  real :: log_NCO, log_nH2, log_temp
  integer :: temp_index, NCO_index, nH2_index
  integer :: order
  real :: c00, c01, c10, c11, c0, c1, dy
  real :: xd, yd, zd
  real :: log_lambda
  real :: sigh2, ve, dvdr
  
  logical, save :: first_call = .true.
  
  
  if(first_call) then
     
     !! read in CO cooling table
     !! if CO_LVG is true, nnco will be the NCO_tilde
     !! parameter used in Glover et al. CO paper
     !! NCO_tilde = n_co / |dvdr|
     !!
     !! if CO_LVG is false, nco will be regular old
     !! CO column density
     
     path = adjustl(trim("CO.dat"))
     if(CO_LVG) then
        path = adjustl(trim("CO_LVG.dat"))
     endif
     open(11, file=path, form='formatted', action='read')
     read(11,*), dummys
     read(11,*), dummys
     read(11,*), dummys
     read(11,*), ntemp
     read(11,*), nnh2
     read(11,*), nnco
 
     read(11,*), dummys
    
     allocate(log_t_table(ntemp))
     allocate(log_nh2_table(nnh2))
     allocate(log_nco_table(nnco))
     allocate(log_lambda_table(ntemp,nnh2,nnco))
     
     
     
     do i=1,ntemp
        do j=1,nnh2
           do k=1,nnco
              read(11,*) val1, val2, val3, val4
              log_t_table(i) = val1
              log_nh2_table(j) = val2
              log_nco_table(k) = val3
              log_lambda_table(i,j,k) = val4
           enddo
        enddo
     enddo
     
     close(11)
     
     lognco_min = minval(log_nco_table)
     lognco_max = maxval(log_nco_table)
     lognh2_min = minval(log_nh2_table)
     lognh2_max = maxval(log_nh2_table)
     logtemp_min = minval(log_t_table)
     logtemp_max = maxval(log_t_table)
     
     
     first_call = .false.

  endif
  
  lambda = 0.0
  
  nh = rpar(1)
  if(CO_LVG) then
     dvdr = rpar(8)
     NCO = y(iCO) * nh / dvdr 
  else
     NCO = rpar(2) * y(iCO)
  endif
  
  
!! effective nH2 - see Glover et al. 2010
  sigh2 = 3.3e-16 / (temperature / 1000.0)**0.25
  ve = 1.03e4 * sqrt(temperature)
  nH2 = (y(iH2) + 1.414 * (2.3e-15 / sigh2) + (1.3e-8 / sigh2 / ve) * y(ie)) * nh
  
  !! it's possible that NCO is really, really small. Make sure it doesn't mess stuff up.
  log_NCO = log10(NCO)
  log_temp = log10(temperature)
  log_nH2 = log10(nH2)
  if(log_NCO .lt. log_nco_table(1)) then
     log_NCO = log_nco_table(1)
  endif

 
  call locate(log_t_table, ntemp, log_temp, temp_index) 
  call locate(log_nh2_table, nnh2, log_nH2, nH2_index)
  call locate(log_nco_table, nnco, log_NCO, NCO_index)
  
  if(temp_index .lt. 1) temp_index = 1
  if(temp_index .eq. ntemp) temp_index = ntemp-1

  if(nH2_index .lt. 1) nH2_index = 1
  if(nH2_index .eq. nnh2) nH2_index = nnh2-1

  if(NCO_index .lt. 1) then
     NCO_index = 1
     !! here and below, we don't want to alter the value of NCO
     !log_NCO = log_nco_table(NCO_index)
  endif
  if(NCO_index .eq. nnco) then
     NCO_index = nnco-1
     !log_NCO = log_nco_table(nnco)
  endif

  
  !! trilinear interpolation for CO cooling rate:
  
  xd = (log_temp - log_t_table(temp_index)) / (log_t_table(temp_index+1) - log_t_table(temp_index))
  yd = (log_nH2 - log_nh2_table(nH2_index)) / (log_nh2_table(nH2_index+1) - log_nh2_table(nH2_index))
  zd = (log_NCO - log_nco_table(NCO_index)) / (log_nco_table(NCO_index+1) - log_nco_table(NCO_index))
  

  c00 = log_lambda_table(temp_index,nH2_index,NCO_index) * (1.0 - xd) & 
      + log_lambda_table(temp_index+1,nH2_index,NCO_index) * xd
  c10 = log_lambda_table(temp_index,nH2_index+1,NCO_index) * (1.0 - xd) & 
      + log_lambda_table(temp_index+1,nH2_index+1,NCO_index) * xd
  c01 = log_lambda_table(temp_index,nH2_index,NCO_index+1) * (1.0 - xd) & 
      + log_lambda_table(temp_index+1,nH2_index,NCO_index+1) * xd
  c11 = log_lambda_table(temp_index,nH2_index+1,NCO_index+1) * (1.0 - xd) & 
       + log_lambda_table(temp_index+1,nH2_index+1,NCO_index+1) * xd
  
  c0 = c00 * (1.0 - yd) + c10 * yd
  c1 = c01 * (1.0 - yd) + c11 * yd
  
  log_lambda = c0 * (1.0 - zd) + c1 * zd
  
  lambda = 10.0**log_lambda * nh * y(iCO) * nh
  
  return
end subroutine CO_cooling



subroutine calculate_temperature(y,ei,nh,temp)

  use CO_CoolGlobal
  
  implicit none

  real, intent(IN) :: y(NTOT), nh, ei
  real, intent(OUT) :: temp
  real :: ntotal, gamma_t, gamma_b, gamma
  
  integer :: i, j

  real :: val1, val2, val3, val4
  integer, save :: nxh2, neint
  real, save, allocatable :: logeintv(:)
  real, save, allocatable :: logxh2v(:)
  real, save, allocatable :: gammav(:,:)
  character(len=128) :: path, dummys
  
  
  
  if(evolve_temp) then
     temp = ei
     !print*, "ctss: shouldn't be here..."
  else
     ntotal = nh * (YHELIUM + 1.0 + y(ie) - y(iH2))
     gamma = 1.66666666667
     temp = (gamma - 1.0) * ei / (ntotal * kb)
  endif
  
  
  return
  
end subroutine calculate_temperature



subroutine get_tdust(y,nh,rpar,tdust)

  use CO_CoolGlobal
  
  implicit none
  
  real, intent(IN) :: y(NTOT), nh, rpar(nrpar)
  real, intent(OUT) :: tdust
  
  real :: temp
  
  tdust = rpar(6)
  
  return
end subroutine get_tdust




subroutine set_timestep(yneq, rpar, h)

  use CO_CoolGlobal

  implicit none

  real, intent(IN) :: yneq(NNEQ), rpar(nrpar)
  real, intent(OUT) :: h

  real :: y(NTOT), YDOT(NNEQ), htry
  integer :: n, ipar(nipar)
  logical :: limiter(NNEQ)
  
  
  call coolderiv(NNEQ,1.0,yneq,ydot,rpar,ipar)
   
  h = 1.0e20

  !timestep only based on eint
  limiter(eint)   = .true.
  
  do n=1, NNEQ
     if(  (abs(ydot(n)) .gt. 0.0) .and. & 
          (yneq(n) .gt. 0.0)      .and. & 
           limiter(n)) then
        htry = abs(yneq(n) / ydot(n))
        h = min(h,htry)
     endif
  enddo
  
  h = h * SAFETY
  
  return
end subroutine set_timestep






subroutine jac
  implicit none
  print*, "dummy function, should not be here"
  return
end subroutine jac


subroutine smoothstep(xl, xh, x0, frac)
  implicit none
  real, intent(IN) :: xl, xh, x0
  real, intent(OUT) :: frac
  real :: x

  x = (x0-xl)/(xh-xl)
  if(x < 0.0) x = 0.0
  if(x > 1.0) x = 1.0
  frac = x*x*(3.0-2.0*x)
  !frac =  x*x*x*(x*(x*6.0 - 15.0) + 10.0)

return
end subroutine smoothstep
