SUBROUTINE OPACTD(IJ)
C     =====================
C
C     Absorption and emission coefficients, and their derivatives
C
C     This procedure is very similar to OPACT1, the only difference is
C     the evaluation of derivatives
C
C     Input: 
C      IJ    - depth index
C     Output:
C      ABSO1 - array of absorption coefficient
C      EMIS1 - array of emission coefficient
C      SCAT1 - array of scattering coefficient
C      Dxxy  - array of derivatives of xx (=AB for absorption, =EM for
C              emission, =SC for scattering) coefficient 
C              wrt y (=T for temperature, =N for density)
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ALIPAR.FOR'
      INCLUDE 'ARRAY1.FOR'
      INCLUDE 'ITERAT.FOR'
      common/rhoder/drhodt(mdepth)
      PARAMETER (DELT=1.D-3, DELR=1.D-3)
      common/dsctva/dsct1(mdepth),dscn1(mdepth)
      common/hmolab/anh2(mdepth),anhm(mdepth)
C
      imodf=0
      if(ifryb.gt.0) imodf=1
      FR=FREQ(IJ)
      DO ID=1,ND
         T=TEMP(ID)
         T1=T*(UN+DELT)
         RHO=DENS(ID)
         RHO1=RHO*(UN+DELR)
         XKF(ID)=EXP(-HKT1(ID)*FR)
         XKF1(ID)=UN-XKF(ID)
         XKFB(ID)=XKF(ID)*BNUE(IJ)
         PLAN=XKFB(ID)/XKF1(ID)
         DPLAN=PLAN/XKF1(ID)*HKT1(ID)*FR/T
         CALL OPCTAB(FR,IJ,ID,T,RHO,AB,SC,SCT,imodf)
         CALL OPCTAB(FR,IJ,ID,T1,RHO,AB1,SC1,SCT1,imodf)
         CALL OPCTAB(FR,IJ,ID,T,RHO1,AB2,SC2,SCT2,imodf)
         ABSO1(ID)=AB+sct
         SCAT1(ID)=SCT
         EMIS1(ID)=AB*PLAN
         ABSOT(ID)=ABSO1(ID)
         IF(IMODF.EQ.0) ABSOT(ID)=ABSO1(ID)/DENS(ID)
c
c        derivatives w.r.t. temperature
c
         DABT1(ID)=(AB1-AB)/T/DELT
         DEMT1(ID)=AB*DPLAN+DABT1(ID)*PLAN
         DSCT1(ID)=(SCT1-SCT)/T/DELT
         dabt1(id)=dabt1(id)+dsct1(id)
         if(ifryb.gt.0) then
c
c        derivatives w.r.t. density    
c
         DABN1(ID)=(AB2-AB)/RHO/DELR
         DEMN1(ID)=DABN1(ID)*PLAN
         DSCN1(ID)=(SCT2-SCT)/RHO/DELR
         dabn1(id)=dabn1(id)+dscn1(id)
c
c        modify derivatives  in case density is not a state parameter
c
         IF(INHE.LE.0) THEN 
            DABT1(ID)=DABT1(ID)+DABN1(ID)*DRHODT(ID)
            DEMT1(ID)=DEMT1(ID)+DEMN1(ID)*DRHODT(ID)
            DSCT1(ID)=DSCT1(ID)+DSCN1(ID)*DRHODT(ID)
            DABN1(ID)=0.
            DEMN1(ID)=0.
            DSCN1(ID)=0.
         END IF
         end if
      END DO
C
C     store quantities for explicit (linearized) frequencies
C
      IF(IJEX(IJ).LE.0.OR.IFRYB.GT.0) RETURN
      IJE=IJEX(IJ)
      DO ID=1,ND
         ABSOEX(IJE,ID)=ABSO1(ID)
         SCATEX(IJE,ID)=SCAT1(ID)
         EMISEX(IJE,ID)=EMIS1(ID)
         DABTEX(IJE,ID)=DABT1(ID)
         DEMTEX(IJE,ID)=DEMT1(ID)
         DABNEX(IJE,ID)=DABN1(ID)
         DEMNEX(IJE,ID)=DEMN1(ID)
      END DO
C
      RETURN
      END