subroutine recomb(flam,t,rho,z,qrec) implicit double precision (a-h,o-z) zeta=1.57d5*z*z/t temp=1.850d-4*(z**6)*((2.d0*z)**(-2))*rho*rho*flam*flam qrec=dexp(-zeta-(2.428d-5*(rho**(0.66667d0))/flam))*temp return end ************************************************************************