SUBROUTINE TRMDER(ID,T,PG,PRAD,TAU,HEATCP,DLRDLT,GRDADB,RHO)
C     ============================================================
C
C     Thermodynamic derivatives
C     Evaluation similar as in Kurucz's ATLAS
C
C     Input:  T     -  temperature
C             PG    -  gas pressure
C             PRAD  -  radiation pressure
C     Output: DEDT  -  d(energy)/d(T)
C             DRDT  -  d(rho)/d(T)
C             DEDPG -  d(energy)/d(PG)
C             DRDPG -  d(rho)/d(PG)
C             RHO   -  density
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      common/derdif/dift,difp
      common/terden/rhoter,anta,entrp
      common/adiaba/grdad0,itgrad
      DIMENSION TT(5),PP(5),RHON(5),ENER(5),entr(5)
C
C   NUMERICAL EVALUATION OF THERMODYNAMIC DERIVATIVES
C
c     prad=0.
      P=PG
      TT(1)=T*(UN+DIFT)
      TT(2)=T*(UN-DIFT)
      TT(3)=T
      TT(4)=T
      TT(5)=T
      PP(1)=P
      PP(2)=P
      PP(3)=P*(UN+DIFP)
      PP(4)=P*(UN-DIFT)
      PP(5)=P
      tmoli0=tmolim
      if(t.lt.tmolim) tmolim=t*(un+dift+0.001)
      if(t.ge.tmolim) tmolim=t*(in-dift-0.001)
      DO I=1,5
         TE=TT(I)
         TKN=TE*1.38054D-16
         ANT=PP(I)/TKN
         CALL ELDENS(ID,TE,ANT,ANE,ENRG,ENT,WM,0)
         RHON(I)=rhoter
         ENER(I)=(1.5D0*PP(I)+ENRG+3.D0*PRAD*(TE/T)**4)/RHON(I)
         entr(i)=ent/RHON(I)
      END DO
      tmolim=tmoli0
      entrp=entr(5)
c
      DRDT=(RHON(1)-RHON(2))/(2.*T*DIFT)
      DRDPG=(RHON(3)-RHON(4))/(2.*P*DIFP)
      RHO=RHON(5)
      DPDPG=1.
      dpdt=0.
      if(tau.lt.50.) DPDT=4.*PRAD/T*(un-exp(-tau))
      DLRDLT=T/RHO*(DRDT-DRDPG*DPDT/DPDPG)
      ptot=pg+prad
c
      if(ifentr.le.0) then
         DEDT= (ENER(1)-ENER(2))/(2.*T*DIFT)
         DEDPG=(ENER(3)-ENER(4))/(2.*P*DIFP)
         HEATCV=DEDT-DEDPG*DRDT/DRDPG
         HEATCP=DEDT-DEDPG*DPDT/DPDPG-
     *          PTOT/RHO/RHO*(DRDT-DRDPG*DPDT/DPDPG)
         GRDADB=-PTOT/RHO/T*DLRDLT/HEATCP
       else
         DSDT= (ENTR(1)-ENTR(2))/(2.*T*DIFT)
         DSDP= (ENTR(3)-ENTR(4))/(2.*P*DIFP)
         grdadb=-dsdp/dsdt*pg/t
c        heatcp=-PTOT/RHO/T*DLRDLT/grdadb
         heatcp=t*dsdt
      end if
c
      if(iter.le.itgrad.and.grdad0.gt.0) 
     *   grdadb=grdad0*0.4+(un-grdad0)*grdadb
      RETURN
      END