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