rad.F
C TODO -- solve for ionization of the gas, don't just assume full ionization.
#define UO  9

C deg K (2.13 kev) : Compton temp of QSO spectrum
#define TCOMP_K 2.5d7 

C Implicit integration for radiative cooling
#define SEMI_IMPLICIT 
C#define IMPLICIT 

#define DT_FACT 0.1

C PHYSICIAL CONSTANTS (c.g.s)
#define MP_CGS  1.6726216d-24
#define KB_CGS  1.3806488d-16

C TRANSLATES CODE UNITS INTO c.g.s
C SCALES: 1kpc; 1Gyr; 2.5*10^7Msun
C LU:length ; TU: time ; MU:  mass
C DU:density; EU:energy; VU:volume
#define LU 3.08568025d21
#define TU 3.15360000d16
#define MU 4.97230000d40

#define DU 1.6924060d-24
#define EU 4.76042745d50
#define VU 2.93800656d64

c=======================================================================
c
c    \\\\\\\\\\      B E G I N   S U B R O U T I N E      //////////
c    //////////        R A D I A T I V E   L O S S        \\\\\\\\\\
c
c=======================================================================

      subroutine ism_rad_loss

      use real_prec
      use param
      use field
      use bndry
      use grid

      implicit none
      LOGICAL  success
      REAL(rl) ::  enew(in,jn,kn)
      INTEGER  ::  i, j, k

      success = .false. 

C     TRY OPTICALLY-THIN APPROXIMATION FIRST (JUST 
C     FOR SAVING COMPUTING TIME). IF "NOT SUCCESS",
C     WE HAVE TO DO THE EXPLICIT SERIAL INTEGRATION
      call rad_loss_co_parallel(success, enew) 

      if (.not.success) then 
      call rad_loss_co_serial  (success, enew) 
      endif

c     UPDATE ENERGY
      do k=ks,ke
         do j=js,je
            do i=is,ie               
               e(i,j,k) = enew(i,j,k)
            enddo
         enddo
      enddo

C     MARK THE BOUNDARY OUT OF DATE
      do i = 1,6
         bvstat(i,2) = 0 
      enddo

      return
      end
c=======================================================================
c    \\\\\\\\\\     E N D   R A D I A T I V E   L O S S    //////////
c=======================================================================

      subroutine ism_rad_read

      use real_prec
      use param
      use root
#ifdef MPI_USED
      use mpiyes
#else
      use mpino
#endif
      use mpipar

      use rad

      implicit none 
      INTEGER num
      parameter (num = 6)
      REAL(rl) ::  buffer(num),  ftol
      namelist /cool/ ilowt,lowt_cexp,lowt_ct0,
     .                lowt_hexp, lowt_ht0, ftol

C     default value for this:
      ilowt = 0
      lowt_cexp = 1.5
      lowt_ct0  = 2e4
      lowt_hexp = 0.0
      lowt_ht0  = 2e4
      ftol = 0.1

      if (myid.eq.0) then 
         rewind(1)
         read (1,cool)
         rewind(1)
         write(2,cool)      

         buffer(1) = ilowt
         buffer(2) = lowt_cexp
         buffer(3) = lowt_ct0
         buffer(4) = lowt_hexp
         buffer(5) = lowt_ht0
         buffer(num) = ftol
      endif

#ifdef MPI_USED      
      call MPI_BCAST(buffer,num,MPI_FLOAT,0,comm3d,ierr)

      ilowt = buffer(1) 
      lowt_cexp = buffer(2)
      lowt_ct0  = buffer(3)
      lowt_hexp = buffer(4)
      ftol  = buffer(5)
      lowt_ht0  = buffer(num)
#endif

      rad_loss_ftol = ftol

      return      
      end !of subroutine rad_loss_co_read

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine rad_loss_co_serial(success, enew)

      use real_prec
      use param
      use grid
      use field
#ifdef MPI_USED
      use mpiyes
#else
      use mpino
#endif
      use mpipar

      use agn
      use rad

      implicit none
      LOGICAL success
      REAL(rl) ::  enew(in,jn,kn)
      INTEGER  ::  i, j, k

C     This can't fail at the moment
      success = .true.

C     Wait for BC from below
      if (coords(1).eq.0) then 
         do k=ks,ke
            do j=js,je
               leff_bh(is-1,j,k) = lum_en
            enddo
         enddo
      else

         call MPI_RECV(leff_bh(is-1,1,1),1,i_slice
     .                ,n1m,16001,comm3d,stat,ierr)

      endif

      call rad_loss_co_try(enew) 

      if (coords(1).lt.ntiles(1)-1) then 
         call MPI_SEND(leff_bh(ie  ,1,1),1,i_slice
     .                ,n1p,16001,comm3d,stat,ierr)
      endif                  

      return      
      end !of subroutine rad_loss_co_serialized

ccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine rad_loss_co_parallel(success, enew) 

      use real_prec
      use param
      use grid
      use field
#ifdef MPI_USED
      use mpiyes
#else
      use mpino
#endif
      use mpipar

      use rad
      use agn

      implicit none
      LOGICAL  ::  success,lsuccess
      REAL(rl) ::  enew(in,jn,kn)
      REAL(rl) ::  lstart, lftol
      INTEGER  ::  i, j, k

      do k=ks,ke
         do j=js,je
            leff_bh(is-1,j,k) = lum_en
         enddo
      enddo

      call rad_loss_co_try(enew) 

      lsuccess = .true.

C     Absorption is cumulative: rad_loss_ftol should set the _total_ lum
C     decrement over the simualtion, then the product of all the tiles
C     should be less than rad_loss_ftol
C     lftol = rad_loss_ftol**(one/ntiles(1)) !orginal form
      lftol = 1.d0 - (1.d0-rad_loss_ftol)**(one/ntiles(1))

      do k=ks,ke
         do j=js,je
            lstart = max(leff_bh(is,j,k), TINY)
            if (abs(leff_bh(ie,j,k)/lstart - 1).gt.lftol) then 
               lsuccess = .false.
            endif
         enddo
      enddo

#ifdef MPI_USED
      call MPI_ALLREDUCE(lsuccess,success,1,MPI_LOGICAL, 
     .         MPI_LAND, comm3d, ierr)
#else
      success = lsuccess
#endif

      return 
      end !of subroutine rad_loss_co_parallel

cccccccccccccccccccccccccccccccccccccccccccc

      subroutine rad_loss_co_try(enew)
C     Radiation cooling/heating after Ciotti + Ostriker 07.  
C     Includes line cooling, free-free cooling, as well as 
C     compton heating and line heating from an AGN.

      use real_prec
      use param
      use grid
      use field
      use root
      use bndry
#ifdef MPI_USED
      use mpiyes
#else
      use mpino
#endif
      use mpipar

      use agn
      use star
      use rad

      implicit none
C     Output
      REAL(rl) ::  enew(in,jn,kn)

C     Local
      EXTERNAL ::  rad_loss_co_code
      REAL(rl) ::  rad_loss_co_code_bisect

      REAL(rl) ::  brem,line,comp,heat,cool,line_heat
      REAL(rl) ::  lambda, lum_eff, cs_sf_min, ee
      REAL(rl) ::  radii, area, asph, vol
      REAL(rl) ::  global_dtcool
      INTEGER  ::  i, j, k
C
C     For now, do this explicitly, but likely will use implicit scheme
      dtline = HUGE
      dtbrem = HUGE
      dtcomp = HUGE
      dtcool = HUGE

C     Cooling only depends on the cell itself, don't need boundary info

C NOTE -- depend on units here.  Get cooling time for unit density gas
C     with sound speed of 16 km/s, near peak of cooling curve.  Use this
C     is the sound speed is lower than 16 km/s, because otherwise the
C     cooling function goes to zero and the cooling time goes to
C     infinity.
      cs_sf_min = 16.0

      do k=ks,ke
         do j=js,je
            do i=is,ie               

               vol  = dvl1b(i)*dvl2b(j)*dvl3b(k)
               area = g2a(i)*g31a(i)*g32b(j)*dx2b(j)*dx3b(k)
               asph = 4*pi*x1b(i)**2

#ifdef SEMI_IMPLICIT
C     Try an explicit step, if it's too big do an implicit step
               radii = sqrt(r2b(i,j,k)) 
               call rad_loss_co_code(d(i,j,k),e(i,j,k),leff_bh(i-1,j,k),
     .                      radii, line,brem,comp,line_heat, heat, cool)
               lambda = brem + comp + line               

               if (dt.lt.abs(DT_FACT*e(i,j,k)/lambda)) then 
                  enew(i,j,k) = e(i,j,k) + dt*lambda 

                  leff_bh(i,j,k) = leff_bh(i-1,j,k) -
     $                 (max(line,zro) + max(comp,zro))*asph*vol/area
                  if (leff_bh(i,j,k).lt.zro) leff_bh(i,j,k) = zro

               else
                  enew(i,j,k) = rad_loss_co_code_bisect(  d(i,j,k),  
     1                                 e(i,j,k),  leff_bh(i-1,j,k), 
     2                                 radii,line_heat,heat,cool  )

                  leff_bh(i,j,k) = leff_bh(i-1,j,k) -
     $                 (max(line,zro) + max(comp,zro))*asph*vol/area
                  if (leff_bh(i,j,k).lt.zro) leff_bh(i,j,k) = zro

               endif

#ifdef STARFORM
               if (lambda.gt.zro) then 
                  cooling_time(i,j,k) = huge

               else if ( ilowt.ge.1 ) then 
C              Low temp cooling: compute cooling time self-consistently
                  cooling_time(i,j,k) = -e(i,j,k)/lambda

               else if ( ilowt.eq.0 .and. (cs_sf_min**2).lt. 
     .                  (gamma*gamm1*e(i,j,k)/d(i,j,k)) ) then
C              No low temp cooling, but gas temp above min temp
                  cooling_time(i,j,k) = -e(i,j,k)/lambda

               else 
C              No low temp cooling, but gas temp below min temp,
C              evaluate cooling terms again
                  ee = d(i,j,k)*cs_sf_min**2/(gamma*gamm1)
                  call rad_loss_co_code( d(i,j,k), ee, leff_bh(i-1,j,k),
     .                     radii, line,brem,comp, line_heat, heat,cool )

                  lambda = brem + comp + line               
                  if (lambda.gt.zro) then 
                     cooling_time(i,j,k) = huge
                  else
C                 Set cooling time given the actual energy density 
                     cooling_time(i,j,k) = -e(i,j,k)/lambda
C                 Set cooling time as though gas has all the energy 
C                 it would have if it were at 16 km/s
C                     cooling_time(i,j,k) = -ee/lambda
                  endif
               endif
#endif /* STARFORM */          

#else /* ! SEMI_IMPLICIT */

#ifdef IMPLICIT
C     USING IMPLICIT SCHEME
               write(UO,*) "Implement SF cooling time tracking!"
               stop
C              in the implicit scheme, dtcool is never reset from HUGE
C              so doesn't influence timestep
               radii = sqrt(r2b(i,j,k))
               enew(i,j,k) = rad_loss_co_code_bisect(  d(i,j,k),
     .                              e(i,j,k),  leff_bh(i-1,j,k), 
     .                              radii,line_heat,heat,cool  )

               leff_bh(i,j,k) = leff_bh(i-1,j,k) -
     .              (max(line,zro) + max(comp,zro))*asph*vol/area
               if (leff_bh(i,j,k).lt.zro) leff_bh(i,j,k) = zro

#else /* ! IMPLICIT */
C     USING EXPLICIT SCHEME
               write(UO,*) "Implement SF cooling time tracking!"
               stop
C              Get cooling rates separately for each proces
               radii = sqrt(r2b(i,j,k))
               call rad_loss_co_code(d(i,j,k),e(i,j,k),leff_bh(i-1,j,k),
     .                       radii, line,brem,comp,line_heat,heat,cool )

               lambda = brem + comp + line
               enew(i,j,k) = e(i,j,k) + dt*lambda               

               leff_bh(i,j,k) = leff_bh(i-1,j,k) -
     .              (max(line,zro) + max(comp,zro))*asph*vol/area
               if (leff_bh(i,j,k).lt.zro) leff_bh(i,j,k) = zro

C              Get cooling timescales to set timestep
               dtcool = min(dtcool, abs(DT_FACT*e(i,j,k)/lambda))
C              These are diagnostic, not dynamical
               dtcomp = min(dtcomp, abs(DT_FACT*e(i,j,k)/comp))
               dtline = min(dtline, abs(DT_FACT*e(i,j,k)/line))
               dtbrem = min(dtbrem, abs(DT_FACT*e(i,j,k)/brem))       

#endif /* IMPLICIT */
#endif /* SEMI_IMPLICIT */

#if defined RADCOOL && defined ACCRETION
C     Star formation code wants to be able to calculate cooling times by
C     calling rad_loss_co_code, even if cooling is not enabled.
C     However, if cooling is not enabled, these global variables are not
C     defined, so the code won't compile.  Thus the enclosing ifdef.
               tot_cool(i,j,k) = cool
               tot_heat(i,j,k) = heat
               tot_line_heat(i,j,k) = line_heat
#endif /* RADCOOL && ACCRETION */

            enddo
         enddo
      enddo

#ifndef    SEMI_IMPLICIT 
#ifndef      IMPLICIT  
C    IF USING EXPLICIT SCHEME, GET THE
C    GLOBAL COOLING-TIME SCALE::dtcool
      call MPI_ALLREDUCE ( dtcool,global_dtcool,1,
     .           MPI_FLOAT, MPI_MIN, comm3d, ierr )
      dtcool = global_dtcool
#endif /*    !IMPLICIT    */
#endif /* EXPLICIT SCHEME */

      return
      end !of subroutine rad_loss_co_try

cccccccccccccccccccccccccccccccccccccccc

      REAL*8 function rad_loss_co_code_bisect
     .                (dens,en_dens,lum_eff,rad,line_heat,heat,cool)
C     IMPLICIT SCHEME FOR RADIATIVE COOLING INTEGRATION, WHICH SOLVES
C     THE EQUATION:   (enew -dt*lambda(enew))/eold -1 == 0

      use real_prec
      use param
      use root

C     Find root by bisection      
      implicit none
      REAL(rl) ::  fderiv, ftol
      INTEGER  ::  maxiter, niter
      parameter   (fderiv=0.01,ftol=0.01,maxiter=100)

      REAL(rl) ::  dens,en_dens,lum_eff,rad
      REAL(rl) ::  line_heat, heat, cool
      REAL(rl) ::  eh, e0, el, fh,f0, fl
      REAL(rl) ::  lambda, ll, bb, cc

C     Function to find root is (enew-dt*lambda)/eold-1=0
C     First identify the appropriate interval to search
      eh = en_dens
      el = en_dens
      call rad_loss_co_code(dens,en_dens,lum_eff,rad,
     .                ll,bb,cc, line_heat, heat,cool)
      lambda = ll + bb + cc

      fh = (eh-lambda*dt)/en_dens - 1
      fl = (el-lambda*dt)/en_dens - 1

      if (lambda.lt.0.0) then 
C     Cooling::  expand the interval downward, 
C     seek spot where function turns negiative
         niter = 0
         do while (niter.lt.maxiter.and.fl.gt.0.0)
            el = 0.9*el
            call rad_loss_co_code(dens,el,lum_eff,rad, 
     .                   ll,bb,cc,line_heat,heat,cool)
            lambda = ll + bb + cc

            fl = (el-lambda*dt)/en_dens - 1
            niter = niter + 1
         enddo

      else
C     Heating::  expand the interval upward, 
C     seek spot where function turns positive
         niter = 0
         do while (niter.lt.maxiter.and.fh.lt.0.0)
            eh = 1.1*eh
            call rad_loss_co_code(dens,eh,lum_eff,rad,
     .                   ll,bb,cc,line_heat,heat,cool)
            lambda = ll + bb + cc

            fh = (eh-lambda*dt)/en_dens - 1
            niter = niter + 1
         enddo
      endif

      if (fl*fh.le.0) then 
C     Check for landing on a root      
         if (fl.eq.0) eh=el
         if (fh.eq.0) el=eh

         f0 = HUGE
         niter = 0
         do while (niter.lt.maxiter.and.abs(f0).gt.ftol) 
            e0 = 0.5*(eh + el)
            call rad_loss_co_code(dens,e0,lum_eff,rad,
     .                   ll,bb,cc,line_heat,heat,cool)      
            lambda = ll + bb + cc

            f0 = (e0-lambda*dt)/en_dens - 1

            if (f0*fl.gt.0) then 
               el = e0
               fl = f0
            endif
            if (f0*fh.gt.0) then 
               eh = e0
               fh = f0
            endif
            niter = niter + 1
         enddo
      endif

      if (niter.eq.maxiter.or.fl*fh.gt.0)  then         
         if (niter.eq.maxiter)  then         
c            write (UO,*) "Warning, exceeded max iter! reverting!"
c            write (UO,*) dens, en_dens, lum_eff, rad
         endif
         if (fl*fh.gt.0)  then         
c            write (UO,*) "Warning, couldn't bracket root! reverting!"
c            write (UO,*) dens, en_dens, lum_eff, rad
         endif

C     IF IMPLICIT SCHEME FAILED, then perform "SAFE EXPLICIT STEP" 
C     to prevent too-large jumps:: using   "g(x)=(1-exp(-x))"  on 
C     the down side to prevent zero crossings and  "log(1+x)"  on
C     the up side to allow fast heating without being ridiculous
         call rad_loss_co_code(dens, en_dens, lum_eff, rad,
     .                         ll,bb,cc,line_heat,heat,cool)
         lambda = ll + bb + cc

         if (lambda.gt.0) then 
            e0 = e0*(1+log(1+lambda*dt/e0))
         else
            e0 = e0*exp(lambda*dt/e0)
         endif
      endif      

      rad_loss_co_code_bisect = e0

      return
      end !of function rad_loss_co_code_bisect

cccccccccccccccccccccccccccccccccccccccccccccccc

C      REAL*8 function rad_loss_co_code_newton

cccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine rad_loss_co_code(dens,en_dens,lum_eff,radii,
     .               line, brem, comp, line_heat, heat, cool)

      use real_prec
      use param
      use root

      use rad

      implicit none
C     Input: dens, en_dens, lum_eff, rad
C     Output: line, brem, comp
C     line = net line heating/cooling
C     brem = brem. cooling
C     comp = compton heating/cooling
C     line_heat = line heating only
C     heat = total heating
C     cool = total cooling

      REAL(rl) ::  mean_molecular_weight
      REAL(rl) ::  ionization_parameter_cgs
      REAL(rl) ::  rad_loss_co_cgs_comp
      REAL(rl) ::  rad_loss_co_cgs_brem
      EXTERNAL ::  rad_loss_co_cgs_line

      REAL(rl) ::  dens, en_dens, lum_eff, radii
      REAL(rl) ::  line, brem, comp,  line_heat 
      REAL(rl) ::  heat, cool, heat_cgs, cool_cgs
      REAL(rl) ::  line_heat_cgs, line_cool_cgs
      REAL(rl) ::  line_cgs, brem_cgs, comp_cgs
      REAL(rl) ::  mmw, dens_cgs, nn_cgs, cs, cs_cgs, temp_cgs
      REAL(rl) ::  lowt_temp_cgs
      REAL(rl) ::  lowt_line_heat_cgs, lowt_line_cool_cgs

      REAL(rl) ::  lum_eff_cgs, rad_cgs, zeta

C     NOTE I believe that the C+O functions want baryon density
      dens_cgs = dens*DU
      nn_cgs   = dens_cgs/MP_CGS

C     But temperature cares about mmw
      cs  = sqrt(gamma*gamm1*en_dens/dens)
      mmw = mean_molecular_weight(cs)

      cs_cgs   = cs*LU/TU
      temp_cgs = mmw*MP_CGS*cs_cgs**2/(gamma*KB_CGS)

      lum_eff_cgs = lum_eff*EU/TU
      rad_cgs = radii*LU
      zeta = ionization_parameter_cgs(nn_cgs,lum_eff_cgs,rad_cgs)

      comp_cgs = rad_loss_co_cgs_comp(nn_cgs,temp_cgs,zeta)

      if (ilowt.eq.0) then 
         call rad_loss_co_cgs_line(nn_cgs, temp_cgs, zeta,
     .               line_heat_cgs,line_cool_cgs,line_cgs)
      else
         if (lowt_ct0.eq.lowt_ht0) then 
C        Calculate cooling as though temp were lowt_t0
            lowt_temp_cgs = temp_cgs +
     $           lowt_ct0/(1+exp((temp_cgs-lowt_ct0)/(0.25*lowt_ct0)))
            call rad_loss_co_cgs_line( nn_cgs, lowt_temp_cgs, zeta,
     $              lowt_line_heat_cgs,lowt_line_cool_cgs,line_cgs )
C           Then correct via power law scaling
            line_cool_cgs = lowt_line_cool_cgs *
     $                    (temp_cgs/lowt_temp_cgs)**lowt_cexp
            line_heat_cgs = lowt_line_heat_cgs *
     $                    (temp_cgs/lowt_temp_cgs)**lowt_hexp
C           And recompute net
            line_cgs = line_cool_cgs + line_heat_cgs
         else
C        Need two separate calls here
            lowt_temp_cgs = temp_cgs +
     $           lowt_ct0/(1+exp((temp_cgs-lowt_ct0)/(0.25*lowt_ct0)))
            call rad_loss_co_cgs_line( nn_cgs, lowt_temp_cgs, zeta,
     $              lowt_line_heat_cgs,lowt_line_cool_cgs,line_cgs )
            line_cool_cgs = lowt_line_cool_cgs *
     $                    (temp_cgs/lowt_temp_cgs)**lowt_cexp

            lowt_temp_cgs = temp_cgs +
     $           lowt_ht0/(1+exp((temp_cgs-lowt_ht0)/(0.25*lowt_ht0)))
            call rad_loss_co_cgs_line( nn_cgs, lowt_temp_cgs, zeta,
     $              lowt_line_heat_cgs,lowt_line_cool_cgs,line_cgs )
            line_heat_cgs = lowt_line_heat_cgs *
     $                    (temp_cgs/lowt_temp_cgs)**lowt_hexp
C           And recompute net
            line_cgs = line_cool_cgs + line_heat_cgs
         endif 
      endif

      brem_cgs = rad_loss_co_cgs_brem(nn_cgs,temp_cgs)

      heat_cgs = zro
      cool_cgs = zro
      if (comp_cgs.gt.zro) heat_cgs = heat_cgs + comp_cgs
      if (comp_cgs.lt.zro) cool_cgs = cool_cgs + comp_cgs
      heat_cgs = heat_cgs + line_heat_cgs
      cool_cgs = cool_cgs + line_cool_cgs
      cool_cgs = cool_cgs + brem_cgs

      comp = comp_cgs*VU*TU/EU
      line = line_cgs*VU*TU/EU
      brem = brem_cgs*VU*TU/EU
      heat = heat_cgs*VU*TU/EU
      cool = cool_cgs*VU*TU/EU

      line_heat=line_heat_cgs*VU*TU/EU

      return
      end

ccccccccccccccccccccccccccccccccccccccccc      

      REAL*8 function rad_loss_co_cgs_comp (nn, temp, zeta)
C     COMPTON HEATING/COOLING RATE DUE TO AGN ILLUMINATING

      use real_prec
      use agn

      implicit none 
      REAL(rl) ::  nn, temp, zeta, comp, C_COMP
C     Compton heating coefficient                              
      C_COMP = 4.1d-35   ! cm^2/K                               

#ifdef   VCOMPT
      comp = nn**2 * C_COMP*(Tcomp  -temp)*zeta        
#endif /*VCOMPT*/
#ifndef  VCOMPT
      comp = nn**2 * C_COMP*(TCOMP_K-temp)*zeta        
#endif

      rad_loss_co_cgs_comp = comp

      return
      end !of subroutine rad_loss_co_code

ccccccccccccccccccccccccccccccccccccccccc

      REAL*8  function rad_loss_co_cgs_brem(nn, temp)
C     COMPUTE THE LOCAL BREMSTRAULLUNG COOLING RATE

      use real_prec

      implicit none
      REAL(rl) ::  C_BREM
      REAL(rl) ::  nn, temp, brem
C     erg cm^3 / s : Bremstraullung coef                                 
C     FIXME -- Osterbrock 74 says = 1.42 sqrt(T) n_e n_z Z^2
C     C+O says that the n here refers to Hydrogen number density
C     Accounting for H and He, I get that the Brem formula is 
C     1.42e-27 sqrt(T) n^2 (1-0.5*Y)  or
C     1.42e-27 sqrt(T) n_H^2 (1-0.5*Y)/(1-Y)**2

C C+O Value
C      C_BREM = -3.8d-27
C Osterbrock w/ my correction for primordial abundance
C      C_BREM = -1.55d-27
C Normalized to Sutherland + Dopita metallicity ~ -0.5
      C_BREM = -2.0d-27
      brem = nn**2 * C_BREM * sqrt(temp)

C     Free-free doens't operate when gas isn't ionized.  On the other
C     hand, doesn't make a difference when line cooling dominates, so
C     turn it off at 10^5 K
      if (temp.lt.1e5) brem = 0

      rad_loss_co_cgs_brem = brem

      return
      end ! of function rad_loss_co_cgs_brem

cccccccccccccccccccccccccccccccccccccccccccc      

      subroutine rad_loss_co_cgs_line(nn, temp, zeta, heat, cool, net)

      use real_prec

      implicit none
C NOTE -- Note that I'm not sure what nn is supposed to be.  The paper
C says it's the hydrogen number density.  Could also be the baryon
C density, which would make the most sense.  Probably is not the
C particle number density (including electrons), but one could define it
C that way.

C Input: number density, temperature, ionization parameter, 
C all in cgs units.
C Output: volumetric cooling rate in CGS units
C
C heat = total heating
C cool = total cooling
C net = heating - cooling
      REAL(rl) ::  nn, temp, zeta, heat, cool, net
      REAL(rl) ::  C_BOUND
      REAL(rl) ::  aa, bb, cc, logtemp, zeta0, bound

C     Bound transitions coefficient                                 
      C_BOUND = 1.d-23  !erg*cm^3/s                                 

C     TODO -- log10 or loge here?
C     TODO 1.8e5 in cc, correct? 
C     TODO 1.5e3 in zeta_0

      logtemp = log10(temp)
      aa = - ( 18/exp(25 *(logtemp - 4.35)**2)
     .       + 80/exp(5.5*(logtemp - 5.2 )**2)
     .       + 17/exp(3.6*(logtemp - 6.5 )**2) )
      bb = 1.7e4*temp**(-0.7)
C     Rewritten cc and zeta0 to make float overflow less likely
      cc    = 1.1 - 1.1/exp(temp/1.8e5) + (7.952d3/temp)**4
      zeta0 = ( 1 / (1.5*temp**(-0.5) + (7.420d4/temp)**2.5)
     .      + (2.e5/temp)**2*(1.d0+80/exp((temp-1e4)/1.5e3)) )

C     Manually take a few limits to avoid generating NaN's
C     Shouldn't need this since c > 5 at T < 10^4
      if (cc.lt.5) then 
         heat = nn**2 * C_BOUND*(bb*(zeta/zeta0)**cc) 
     .        / (1 + (zeta/zeta0)**cc)
         cool = nn**2 * C_BOUND* aa
     .        / (1 + (zeta/zeta0)**cc)
      else 
         if (zeta.lt.zeta0) then 
            heat = 0.0
            cool = nn**2 * C_BOUND * aa
         else
            heat = nn**2 * C_BOUND * bb
            cool = 0.0
         endif
      endif
      net = heat + cool

      return
      end !of rad_loss_co_cgs_line

cccccccccccccccccccccccccccccccccc

      REAL*8 function ionization_parameter_cgs(nn,lum_eff,rad)
C     Compute ionization parameter given density and effective                 
C     optical luminosity as a function of radius
      use real_prec

      implicit none
      REAL(rl) ::  nn, lum_eff, rad

      ionization_parameter_cgs = lum_eff/(nn*rad**2)

      return
      end

ccccccccccccccccccccccccccccccccccc

      REAL*8  function mean_molecular_weight(cs)

      use real_prec

      implicit none
      REAL(rl) ::  cs
      REAL(rl) ::  result
      REAL(rl) ::  mul, muh, cl, ch
C     Very simple handling of mean molecular weight + ionization.
C     NOTE -- depend on units of velocity here
      mul = 1.23
      muh = 0.6
      cl = 10
      ch = 20
      if (cs.lt.cl) then 
         result = mul
      else if (cs.gt.ch) then 
         result = muh
      else
         result = (muh-mul)*(cs-cl)/(ch-cl) + mul
      endif

      mean_molecular_weight = result

      return
      end !of function mean_molecular_weight
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc