SUBROUTINE iorad(DL,TL,X,Y,Z,C12,X16,FKRAD)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Compute radiative opacity using Iben's 1975 fits to
C the Cox+Stewart opacities, in conjunction with the Christy opacity
c fits
C
C INPUT:   DL = log density (g/cm**3)
c          TL = log temperature (K)
c           X = hydrogen mass fraction
c           Y = He4 mass fraction
c           Z = mass fraction of all "metals"
c         C12 = C12 mass fraction
c         X14 = N14 mass fraction
c         X16 = O16 mass fraction
c
c OUTPUT: fkrad = radiative opacity
c        dlkrro = d(log kappa rad)/d(log rho) at constant T
c         dlkrt = d(log kappa rad)/d(log T) at constant rho
c
      dimension Xsaha(5),Ysaha(11,5)
      FLT6=TL-6.d0
      T6=10.d0**(FLT6)
      T=T6*1.d6
      RHO=10.d0**(DL)
      Zresid=1.d0-X-Y-C12-X16
      FMU = X+Y+3.d0*C12+4.d0*X16+(4.528d0)*Zresid - 1.d0
      ZB=C12+X16
c Electron scattering opacity
      D=0.05d0*(FLT6-1.7d0)
      FKE=(0.2d0-D-DSQRT(D*D+0.0004d0))*(1.d0+X)
c
c Iben Equations (B3)-(B9)
      A=(1.25d0+0.488d0*dsqrt(FMU)+0.092d0*FMU)/0.67d0
      B=3.86d0+0.252d0*dsqrt(FMU)+0.018d0*FMU
      C=(2.019d-4*RHO/(T6**1.7d0))
      C=C**2.425d0
      AP=1.d0+C*(1.d0+C/24.55d0)
      FLDB=-A+B*FLT6
c Iben Equation (B2)
      FLK1=0.67d0*(DL-FLDB)+dlog10(AP)
      FK1=10.d0**FLK1
c
c Iben Equations (B12)-(B20)
c
c note that FLKRHO below is from Iben's code directly and doesn't
c appear in the paper
      DLT=0.15d0*(0.13d0+0.212d0*X+FLT6)
      FLKRHO=-0.174d0-DLT+DSQRT(DLT*DLT+0.002809d0)
      FKRHO=10.d0**FLKRHO
      FLD1=-2.3d0+4.833d0*FLT6
      FLD2=-(1.86d0+0.3125d0*X)+3.86d0*FLT6
      FLDB=DMIN1(FLD1,FLD2)
      IF (T6.GT.1.d0) THEN
        FLD0=-(1.68d0+0.35d0*X)+1.8d0*FLT6
      ELSE
        FLD0=-(1.68d0+0.35d0*X)+(3.42d0-0.52d0*X)*FLT6
      ENDIF
      FAC=FLT6-0.22d0+0.1375d0*X
      A1=1.22d0*dexp(-(1.74d0-0.755d0*X)*(FAC*FAC))
      FAC=FLT6-0.18d0+0.1625d0*X
      W=4.05d0*dexp(-(0.306d0-0.04125d0*X)*(FAC*FAC))
      AZ0=A1*dexp(-2.76d0*(DL-FLD0)*(DL-FLD0)/(W*W))
      AZ=AZ0*(Z+ZB)/0.02d0
c Iben Equation (B12)
      FLK2=FKRHO*(DL-FLDB)+AZ
      FK2=10.d0**FLK2
c
c Christy radiative opacity. Note: need PE for this part.  Here we
c use the GETEOS routine to get it...
c   The Christie opacity begins to suck at densities of 
c   about 1.d-3 and above when T is of order 1e5.
c   Can be off by a factor of two before that spike...
c
      IF (T6.le.1.5d0) then
        Xsaha(1)=X
        Xsaha(2)=Y
        V=1.d0/RHO
        call GETEOS(T,V,5.5d0,x,y,z,p,pe,pv,pt,et,e,totln)
        v12=dsqrt(v)
        v14=dsqrt(v12)
c-- The commented lines here are Icko's original version-------
c        t1=t/1.d4
c        ses=(5.4d-13)*v/t1
c        trt=sqrt(t1)
c        t2=t1*t1
c        t3=t2*t1
c        t4=t3*t1
c        d1=20.+5.*t3+t4
c        sz=z/(trt*d1)
c        t46=t4*t2
c        d1=(1.4d3)*t1
c        d2=d1+t46
c        d3=1.d6+0.1d0*t46
c        sy=1./d2+1.5/d3
c        sy=y*sy
c        t10=t46*t4
c        d1=2.d6+2.1*t10
c        s1x=trt*t4/d1
c        t4v14=t4*v14
c        d1=1.d0+0.05d0*t4v14
c        d2=4.5*d1*t3
c        d3=d2+250.d0
c        s2x=d1/(t3*d3)
c        sx=(s1x+s2x)*x
c        fkch=ses+sx+sy+sz
c        fkch=pe*fkch
c ----------------------------------------------------------------
        T4=T/1.d4
        FLNT4=dlog(T4)
        T44=dexp(4.d0*FLNT4)
        T45=dexp(5.d0*FLNT4)
        T46=dexp(6.d0*FLNT4)
c
        TEMPX=4.0d-3/T44+2.0d-4*v14
        TEMPX=1.d0/(4.5d0*T46+1.d0/T4/TEMPX)
        TEMPX=X*(TEMPX+dsqrt(T4)/(2.0d6/T44+2.1d0*T46))
        TEMPY=1.d0/(1.4d3*T4+T46)+1.5d0/(1.0d6+0.1d0*T46)
        TEMPY=Y*TEMPY
        TEMPZ=Z*dsqrt(T4)/(20.d0*T4+5.d0*T44+T45)
        FKCH=PE*(5.4d-13*V/T4+TEMPX+TEMPY+TEMPZ)
      ELSE
        FKCH=0.d0
      ENDIF
c
      IF (X.GT.0.D0.AND.T6.GE.1.5D0) THEN
c        print *," Case 1"
        FKRAD=FKE+FK2
      ELSEIF (X.GT.0.d0.AND.T6.LE.1.d0) THEN
c        print *," Case 2"
        FKRAD=FKCH
      ELSEIF (X.GT.0.d0.AND.T6.LT.1.5d0.AND.T6.GT.1.d0) THEN
c        print *," Case 3"
        FKRAD=2.d0*FKCH*(1.5d0-T6)+2.d0*(FK2+FKE)*(T6-1.d0)
      ELSEIF (X.EQ.0.d0.AND.ZB.GT.Z) THEN
c        print *," Case 4"
        FKRAD=FKE+FK1
      ELSEIF (X.EQ.0.d0.AND.ZB.LE.Z.AND.T6.GE.1.d0) THEN
c        print *," Case 5"
        FKRAD=(FKE+FK1)*ZB/Z+(FKE+FK2)*(1.d0-ZB/Z)
      ELSE
        print *,' Tl=',TL,' Dl=',dl
        print *,' X=',X,'Y=',Y,'C=',C12,'O=',O16
        PAUSE "Impossible logic in Iben-style Opacity"
      ENDIF
      RETURN
      END 
c