C----------------------------------------------------------------------
      SUBROUTINE THERMO(PR,TE,RHO,GRAD,CSPE,PMOL)
C----------------------------------------------------------------------
C
C You have  to call  the sub  for the  first time  just to read the
C tables.
C From the second call on, you must enter with GAS PRESSURE
C (PR. Be  careful: NOT  TOTAL PRESSURE!)  in cgs,  and temperature
C (TE), and you  will get as  an output: density  (RHO),
C adiabatic gradient  (GRAD), Cp  at constant  pressure (CSPE)  and
C molecular weight (PMOL).  In COMMON /01/ you  must also enter  the
C abundances (in  mass fraction)  of: Hydrogen;  He_3; He_4;  C_12,
C C_13, N_14,  N_15, 0_16  and O_17.  Of course,  if you don't have
C some of these elements, simply enter zeros.
C
C The  decrease of  Cp below  the Debye temperature is
C already included.
C
      IMPLICIT DOUBLE PRECISION  (A-H,O-Z)
C
      DIMENSION TETA(200),DEB(200)
C
      COMMON /NUMB01/ X,Y3,Y,CA,CA13,AZ,AZ15,OX
C
      COMMON /NUMB02/ TT(49),PP(49),XX(7),DEN(7),DAD(7),CSP(7),WMO(7),
     @                RO(7,49,33),AD(7,49,33),CP(7,49,33),PM(7,49,33),
     @                IT,IP(4)
C
      SAVE IZ
      DATA IZ /0/
C
      IF(IZ.EQ.0)THEN
        OPEN(UNIT=8,FILE='thermo.tab')
        IZ=1
        READ(8,111)TETA(1),DTETA
        DO K=2,200
         TETA(K)=TETA(K-1)+DTETA
        END DO
        READ(8,111)DEB
        DO 3 K=1,7
        DO 3 J=1,49
        IF(K.LT.7.OR.J.GE.30)THEN
          READ(8,222)(RO(K,J,N),N=1,33),(AD(K,J,N),N=1,33),(CP(K,J,N),
     *    N=1,33),(PM(K,J,N),N=1,33),TT(J),PP(J),XX(K)
        ELSE
          DO 2 N=1,33
          RO(7,J,N)=RO(6,J,N)
          AD(7,J,N)=AD(6,J,N)
          CP(7,J,N)=CP(6,J,N)
    2     PM(7,J,N)=1.777778
        END IF
    3   CONTINUE
        CLOSE(UNIT=8)
        DO 15 K=1,7
        DO 15 J=1,49
        TU=10.D0**(TT(J))
        DO 15 N=1,33
        CSPU=10.D0**(CP(K,J,N))
        RHU=10.D0**(RO(K,J,N))
        PMUL=PM(K,J,N)
        TET=DLOG10(1.74D3*DSQRT(RHU)/TU)
        IF(TET.GT.-3.D0) THEN
          DO LD=2,200
           IF(TET.LT.TETA(LD))GO TO 12
          END DO
12        IF(LD.GT.200) LD=200
          LI=LD-1
          RADEB=(TET-TETA(LI))/DTETA
          DEBYE=10.D0**(DEB(LI)+RADEB*(DEB(LD)-DEB(LI)))
          CSPU=CSPU*DEBYE
        END IF
        TE3=TU*TU*TU
        CPRAD=5.0433D-14*TE3/RHU+1.2312D-20*TE3*TE3*PMUL/(1.D16*RHU*RHU)
        CSPU=CSPU+CPRAD
        CP(K,J,N)=DLOG10(CSPU)
   15   CONTINUE
        RETURN
      END IF
C
C Termine la fase de lectura inicial
C
      P=DLOG10(PR)
      T=DLOG10(TE)
C
      IF(X.GT.0.D0) THEN
       ZSTAR=1.D0-X-Y3-Y
       DO K=2,5
        IF(X.GE.XX(K)) GO TO 5
       END DO
5      L=K
       IF(L.GT.5) L=5
       LL=L-1
       RATCC=(X-XX(L))/(XX(LL)-XX(L))
      ELSE
       ZSTAR=1.D0-Y-OX
       IF(Y.GT.0.D0) THEN
         LL=5
         L=7
        ELSE
         LL=6
         L=7
       END IF
      END IF
C
      DO K=2,49
       IF(T.LT.TT(K)) GO TO 7
      END DO
    7 IT=K-2
      IF(IT.LT.1)IT=1
      IF(IT.GT.46)IT=46
C
      DO K=1,4
       ITK= IT+K-1
       DO J=1,33
        PB=PP(ITK)+.5D0*DFLOAT(J-1)
        IF(P.LE.PB) GO TO 9
       END DO
9      IP(K)=J-2
       IF(IP(K).LT.1)  IP(K)=1
       IF(IP(K).GT.30) IP(K)=30
      END DO
C
      CALL TERSER(LL,L,RHO,GRAD,CSPE,PMOL,P,T,ZSTAR,RATCC)
C
111   FORMAT(10F8.4)
222   FORMAT(8F10.6)
C
      RETURN
      END
C----------------------------------------------------------------------
      SUBROUTINE TERSER(LL,L,RHO,GRAD,CSPE,PMOL,P,T,ZSTAR,RATCC)
C----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION  (A-H,O-Z)
C
      PARAMETER (UM= 2.302585093D0)
C
      COMMON /NUMB01/ X,Y3,Y,CA,CA13,AZ,AZ15,OX
C
      COMMON /NUMB02/ TT(49),PP(49),XX(7),DEN(7),DAD(7),CSP(7),WMO(7),
     @                RO(7,49,33),AD(7,49,33),CP(7,49,33),PM(7,49,33),
     @                IT,IP(4)
C
      COMMON /NUMB03/ ACUB(4),BCUB(4)
C
C  Ojo este DIMENSION no estaba en lo que nos mando Mazzitelli !!!!
C
      DIMENSION VAR(4)
C
      PMOL=0.D0
c     LLEGA=3
      LLEGA=4  ! calcula el peso molecular
C
      DO 22 M=1,LLEGA
      DO 17 K=LL,L
C
      GO TO(1,5,9,13),M
C
1      DO N=1,4
       ITN1=IT+N-1
       DO J=1,4
        IPNJ= IP(N)+J-1
        ACUB(J)=RO(K,ITN1,IPNJ)
        BCUB(J)=PP(ITN1)+.5D0*DFLOAT(IPNJ-1)
       END DO
       CALL CUBICA(ALFA,BETA,GAMMA,DELTA)
       VAR(N)=((ALFA*P+BETA)*P+GAMMA)*P+DELTA
      END DO
      DO J=1,4
       ACUB(J)=VAR(J)
       BCUB(J)=TT(IT+J-1)
      END DO
      CALL CUBICA(ALFA,BETA,GAMMA,DELTA)
      DEN(K)=((ALFA*T+BETA)*T+GAMMA)*T+DELTA
      GO TO 17
C
5     DO N=1,4
       ITN1=IT+N-1
       DO J=1,4
        IPNJ= IP(N)+J-1
        ACUB(J)=AD(K,ITN1,IPNJ)
        BCUB(J)=PP(ITN1)+.5D0*DFLOAT(IPNJ-1)
       END DO
       CALL CUBICA(ALFA,BETA,GAMMA,DELTA)
       VAR(N)=((ALFA*P+BETA)*P+GAMMA)*P+DELTA
      END DO
      DO J=1,4
       ACUB(J)=VAR(J)
       BCUB(J)=TT(IT+J-1)
      END DO
      CALL CUBICA(ALFA,BETA,GAMMA,DELTA)
      DAD(K)=((ALFA*T+BETA)*T+GAMMA)*T+DELTA
      GO TO 17
C
9     DO N=1,4
       ITN1=IT+N-1
       DO J=1,4
        IPNJ= IP(N)+J-1
        ACUB(J)=CP(K,ITN1,IPNJ)
        BCUB(J)=PP(ITN1)+.5D0*DFLOAT(IPNJ-1)
       END DO
       CALL CUBICA(ALFA,BETA,GAMMA,DELTA)
       VAR(N)=((ALFA*P+BETA)*P+GAMMA)*P+DELTA
      END DO
      DO J=1,4
       ACUB(J)=VAR(J)
       BCUB(J)=TT(IT+J-1)
      END DO
      CALL CUBICA(ALFA,BETA,GAMMA,DELTA)
      CSP(K)=((ALFA*T+BETA)*T+GAMMA)*T+DELTA
      GO TO 17
C
13    DO N=1,4
       ITN1=IT+N-1
       DO J=1,4
        IPNJ= IP(N)+J-1
        ACUB(J)=PM(K,ITN1,IPNJ)
        BCUB(J)=PP(ITN1)+.5D0*DFLOAT(IPNJ-1)
       END DO
       CALL CUBICA(ALFA,BETA,GAMMA,DELTA)
       VAR(N)=((ALFA*P+BETA)*P+GAMMA)*P+DELTA
      END DO
      DO J=1,4
       ACUB(J)=VAR(J)
       BCUB(J)=TT(IT+J-1)
      END DO
      CALL CUBICA(ALFA,BETA,GAMMA,DELTA)
      WMO(K)=((ALFA*T+BETA)*T+GAMMA)*T+DELTA
   17 CONTINUE
C
      GO TO(18,19,20,21),M
C
18    IF(L.NE.7) THEN
        A=1.D0/DEXP(UM*DEN(L))
        RHO=1.D0/(A+(1.D0/DEXP(UM*DEN(LL))-A)*RATCC)
       ELSE
        RHO=ZSTAR/DEXP(UM*DEN(6))+OX/DEXP(UM*DEN(7))
        IF(LL.EQ.5) RHO=RHO+Y/DEXP(UM*DEN(5))
        RHO=1.D0/RHO
      END IF
      GO TO 22
C
19    IF(L.NE.7) THEN
        GRAD=DAD(L)*(1.D0-RATCC)+DAD(LL)*RATCC
       ELSE
        GRAD=ZSTAR*DAD(6)+OX*DAD(7)
        IF(LL.EQ.5) GRAD=GRAD+Y*DAD(5)
      END IF
      GO TO 22
C
20    IF(L.NE.7) THEN
        A=DEXP(UM*CSP(L))
        CSPE=DEXP(UM*(CSP(L)*(1.D0-RATCC)+RATCC*CSP(LL)))
       ELSE
        CSPE=ZSTAR*DEXP(UM*CSP(6))+OX*DEXP(UM*CSP(7))
        IF(LL.EQ.5) CSPE=CSPE+Y*DEXP(UM*CSP(5))
      END IF
       GO TO 22
C
21    IF(L.NE.7) THEN
        PMOL=WMO(L)*(1.D0-RATCC)+WMO(LL)*RATCC
       ELSE
        PMOL=ZSTAR*WMO(6)+OX*WMO(7)
        IF(LL.EQ.5) PMOL=PMOL+Y*WMO(5)
      END IF
C
22    CONTINUE
C
      RETURN
      END
C----------------------------------------------------------------------
      SUBROUTINE CUBICA(ALFA,BETA,GAMMA,DELTA)
C----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /NUMB03/ A1,A2,A3,A4,B1,B2,B3,B4
C
      D1=B2*(B2+B1)
      C1=B3*(B3+B1)-D1
      C2=B3-B2
      D2=(A2-A1)/(B2-B1)
      C3=(A3-A1)/(B3-B1)-D2
      C4=B4*(B4+B1)-D1
      C5=B4-B2
      C6=(A4-A1)/(B4-B1)-D2
      D3=C1*C5-C2*C4
      ALFA=(C3*C5-C2*C6)/D3
      BETA=(C1*C6-C4*C3)/D3
      GAMMA=D2-ALFA*(D1+B1*B1)-BETA*(B2+B1)
      DELTA=A1-((ALFA*B1+BETA)*B1+GAMMA)*B1
C
      RETURN
      END
