subroutine opacder(d,t,opcs,dlklr,dlklt)
*
* opacities, but from Lagrange re-interpolator
* Call for each chemical species. If pure, use derivatives from
* within the routine. Otherwise, use 5 point differencing.
* See opac3 for further details.
*
      implicit double precision(a-h,o-z)

      common/je/je
      common/comp/amhyhe,amheca,x(4)
      common/opcswch/first
      logical first
      real*8 op(4),dkrx(4),ytb(4),xtb(4)
      real*8 ta(30),tb(30),ty(30,30,3) 
*
* read in the opacities if we haven't already
*
      if ( first ) then
        first = .false.
        read(25,*) iaend,jbend
        read(25,705) (ta(li),li=1,iaend)
        read(25,705) (tb(li),li=1,jbend)
705     format(8f10.5)
710     format(10f8.4)
        do 1000 lk=1,3
        do 1000 li=1,iaend
        read(25,710) (ty(li,lj,lk),lj=1,jbend) 
1000    continue
      endif
*
* locate our place on the grid.
*
      ider=1
      if(km.lt.1)then
        i = 9
        j = 7
        km = 0
      endif
65    continue
      km=km+1
*
* set temp = xa and dens = xb
*
      xa = t
      xb = d

      if( xa .lt. ta(1) ) go to 50
      if( xa .gt. ta(iaend)) go to 50
1     if ( xa - ta(i)) 4, 6, 2
2     if (xa - ta(i+1)) 6, 5, 3
3     i = i + 2
      go to 1
4     i = i - 1
      go to  1
5     i = i + 1
6     continue

      if( xb .lt. tb(1) ) go to 50
      if( xb .gt. tb(jbend)) go to 50
11    if (xb - tb(j)) 14, 16, 12
12    if(xb-tb(j+1)) 16, 15, 13
13    j = j + 2
      go to 11
14    j = j - 1
      go to 11
15    j = j + 1
16    continue
      itop=0
      it10=i-1
      if(it10 .ge. 1) go to 150
      it10=it10+1
      itop=1
150   continue
      if(ty(it10,j,je).gt.-7.99 .or. ty(it10,j+1,je).gt.-7.99) goto 160
      it10=it10+1
      itop=1
      if(it10 .gt.i) go to 400
      go to 150
160   it1f=it10+3
      if(it1f.le.iaend) goto 170
      if(itop .ne. 0) go to 400
      it10=it10-1
      goto 150
170   continue
      if(ty(it1f,j,je).gt.-7.99 .or. ty(it1f,j+1,je).gt.-7.99)goto 180
      if(itop .ne. 0) go to 400
      it10 = it10-1 
      if(it10.le. 0) go to 400
      if(it1f .le. i+1) go to 400
      go to 150
180   continue

      do 250 int=1,4
      it1=it10+int-1
      if(ty(it1,j,je).le.-7.99 .or. ty(it1,j+1,je).le.-7.99) go to 200
      js1=j-1
      if(js1 .le. 0) js1=js1+1
      if(ty(it1,js1,je) .le. -7.99) js1=js1+1
      if(js1+3 .gt. jbend) js1=js1-1
      if(ty(it1,js1+3,je) .le. -7.99) js1=js1-1
      do 190 k=1,4
      k1=js1+k-1
*
* add conductive opac. to rad. opac.
*
      ocond = ocon(tb(k1),ta(it1))
      ytb(k)=dlog10(1./(10.**(-ty(it1,k1,je))+10.**(-ocond)))
190   xtb(k)=tb(k1) 
      call fh4(xb,ytb,xtb,op(int),dkrx(int),ider) 
      go to 250

200   if(ty(it1,j,je) .le. -7.99) js1=j+1
      if(ty(it1,j+1,je) .le. -7.99) js1=j-1
*
* add conductive opac. to rad. opac.
*
      ocond1 = ocon(tb(js1+1),ta(it1))
      ocond2 = ocon(tb(js1),ta(it1))
      op1=dlog10(1./(10.**(-ty(it1,js1+1,je))+10.**(-ocond1)))
      op2=dlog10(1./(10.**(-ty(it1,js1,je))+10.**(-ocond2)))
      dkrx(int)=(op1-op2)/(tb(js1+1)-tb(js1))
      op(int)=op1+dkrx(int)*(xb-tb(js1+1))
250   continue
      do 260 k=1,4
260   xtb(k)=ta(it10+k-1)
      call fh4(xa,op,xtb,opc,dlklt,ider)
      opcs=opc
      if(ider .ne. 1)then
        return
      endif
      call fh4(xa,dkrx,xtb,dlklr,dumfh,0)
      return
400   continue

      do 500 k=1,2
      it1=i+k-1
      if(ty(it1,j,je).le.-7.99 .and. ty(it1,j+1,je).le.-7.99) go to 450
      js1=j
      if(ty(it1,j,je).le.-7.99) js1=j+1
      if(ty(it1,j+1,je).le.-7.99) js1=j-1
      frxs=(xb-tb(js1+1))/(tb(js1+1)-tb(js1))
*
* add conductive opac. to rad. opac.
*
      ocond1 = ocon(tb(js1+1),ta(it1))
      ocond2 = ocon(tb(js1),ta(it1))
      op1=dlog10(1./(10.**(-ty(it1,js1+1,je))+10.**(-ocond1)))
      op2=dlog10(1./(10.**(-ty(it1,js1,je))+10.**(-ocond2)))
      op(k)=op1+(op1-op2)*frxs
      if(ider .ne. 1) go to 500
*
* add conductive opac. to rad. opac.
*
      dop1=(op1-op2)/(tb(js1+1)-tb(js1))
      dop2=dop1
      if(js1-1 .le. 0) go to 420
      if(ty(it1,js1-1,je) .gt. -7.99) then
*
* add conductive opac. to rad. opac.
*
      ocond1 = ocon(tb(js1-1),ta(it1))
      ocond2 = ocon(tb(js1),ta(it1))
      op1=dlog10(1./(10.**(-ty(it1,js1-1,je))+10.**(-ocond1)))
      op2=dlog10(1./(10.**(-ty(it1,js1,je))+10.**(-ocond2)))
      dop1=0.5*(dop1+(op2-op1)/(tb(js1)-tb(js1-1)))
      endif
420   if(js1+2 .gt. jbend) go to 430
      if(ty(it1,js1+2,je).gt.-7.99) then
*
* add conductive opac. to rad. opac.
*
      ocond1 = ocon(tb(js1+1),ta(it1))
      ocond2 = ocon(tb(js1+2),ta(it1))
      op1=dlog10(1./(10.**(-ty(it1,js1+1,je))+10.**(-ocond1)))
      op2=dlog10(1./(10.**(-ty(it1,js1+2,je))+10.**(-ocond2)))
      dop2=0.5*(dop2+(op3-op1)/(tb(js1+2)-tb(js1+1)))
      endif
430   dkrx(k)=dop2+(dop2-dop1)/(tb(js1+1)-tb(js1))*(xb-tb(js1+1))
      go to 500

450   continue
      jc=1
      jco=0
      js1=j
      if(j.gt. 11) jc=-1
455   do 460 ljs=1,jbend
      js1=js1+jc
      if(js1.ge.jbend .or. js1.le.1) go to 465
      if(ty(it1,js1,je) .gt. -7.99) go to 470
460   continue
465   if(jco .ne. 0) go to 50 
      jco=jco+1
      jc=-jc
      go to 455
470   if(jc .eq. -1) js1=js1-1
*
* add conductive opac. to rad. opac.
*
      ocond1 = ocon(tb(js1+1),ta(it1))
      ocond2 = ocon(tb(js1),ta(it1))
      op1=dlog10(1./(10.**(-ty(it1,js1+1,je))+10.**(-ocond1)))
      op2=dlog10(1./(10.**(-ty(it1,js1,je))+10.**(-ocond2)))
      dkrx(k)=(op1-op2)/(tb(js1+1)-tb(js1))
      op(k)=op1+dkrx(k)*(xb-tb(js1+1))

500   continue
      dlklt=(op(2)-op(1))/(ta(i+1)-ta(i))
      opcs=op(2)+dlklt*(xa-ta(i+1))
      if(ider .eq. 0) return
      dlklr=dkrx(2)+(dkrx(2)-dkrx(1))*(xa-ta(i+1))/(ta(i+1)-ta(i))
      return

50    continue
      stop
      end 

************************************************************************