API tools faq
paste
Advertisement
rgedies

micro

rgedies
May 26th, 2023 (edited)
528
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Fortran 48.78 KB | Source Code | 0 0
raw download clone embed print report
  1.     SUBROUTINE THERMODYNAMICS(XCOMPOS,XFE,XPE,TMS,TBS,AE1,AE3,TDISS&,AE4,Tsolid,Tliquid,as1,krate)
  2.     IMPLICIT REAL*8 (A-H,O-Z)
  3.     PARAMETER(ZERO=1.0D-6)
  4.     REAL*8 XFE,XPE,CEUT,TMS,TBS,AE1,AE3,TDISS,AE4,Tsolid,Tliquid,as1
  5.     REAL*8 C,Mn,Si,Ni,Cr,Mo,Cu,W,V,Nb,Ti,Al,N,FC, PC, BC
  6.     REAL*8 XCOMPOS(*), KRATE(*)
  7.  
  8.       C = XCOMPOS(1)
  9.     Mn = XCOMPOS(2)
  10.     Si = XCOMPOS(3)
  11.     Ni = XCOMPOS(4)
  12.     Cr = XCOMPOS(5)
  13.     Mo = XCOMPOS(6)
  14.     Cu = XCOMPOS(7)
  15.     W = XCOMPOS(8)
  16.     V = XCOMPOS(9)
  17.     Nb = XCOMPOS(10)
  18.     Ti = XCOMPOS(11)
  19.     Al = XCOMPOS(12)
  20.     N = XCOMPOS(13)
  21.  
  22.     CALL THERMCALC (C,Mn,Si,Ni,Cr,Mo,Cu,W,V,Nb,Ti,Al,N,&  XFE,XPE,CEUT,AE1,AE3,TDISS,AE4,Tsolid,Tliquid,as1)
  23.     TBS=637.D0-58.D0*C-35.D0*Mn-15.D0*NI-34.D0*Cr-41.D0*Mo
  24.     TMS=539.D0-423.D0*C-30.4D0*Mn-17.7D0*Ni-12.1D0*Cr-7.5D0*Mo-7.5D0*Si
  25.     FC=EXP(-0.281D0+1.780D0*Mn+0.306D0*Si+1.123D0*Ni&  +2.720D0*Cr+4.057D0*Mo)
  26.     PC=EXP(-4.246D0+4.356D0*Mn+0.441D0*Si+1.705D0*Ni &  +3.329D0*Cr+5.194D0*SQRT(Mo))
  27.     BC=EXP(-10.230D0+10.177D0*C+0.847D0*Mn+0.547D0*Ni &   +0.902D0*Cr+0.364D0*Mo)
  28.     KRATE(1) = FC
  29.     KRATE(2) = PC
  30.     KRATE(3) = BC
  31.     RETURN
  32.     END
  33.  
  34.     SUBROUTINE THERMCALC (C,Mn,Si,Ni,Cr,Mo,Cu,W,V,Nb,Ti,Al,N,& XFE,XPE,CEUT,AE1,AE3,TDISS,AE4,Tsolid,Tliquid,as1)
  35.     IMPLICIT REAL*8 (A-H,O-Z)
  36.     REAL*8 C,Mn,Si,Ni,Cr,Mo,Cu,W,V,Nb,Ti,Al,N,X(0:10)
  37.     REAL*8 AS1,AE1,AE3,TDISS,AE4,TNJ,TNH,TAH,THJ,THB,TJE,TAB,TBC
  38.     CALL XCOMP(C,MN,SI,NI,CR,MO,CU,W,V,NB,X)
  39.     CALL AS1CALC(X,AS1,CEUT,CF)
  40.     CALL AE1CALC(X,AE1,CEUTU,CFU)
  41.     IF (AE1 .LT. AS1) THEN
  42.       TEMP=AE1
  43.       AE1=AS1
  44.       AS1=TEMP
  45.       CEUT=CEUTU
  46.       CF=CFU
  47.     ENDIF
  48.  
  49.     XFE=(CEUT-C)/(CEUT-CF)
  50.     XPE=(C-CF)/(CEUT-CF)
  51.  
  52.     IF (C .GT. CEUT) THEN
  53.        WRITE(*,*) '*** WARNING ***'
  54.        WRITE(*,*) 'Composition of the steel is out of the range of ', &              'this model!'
  55.          WRITE(*,*) 'No results available!'
  56.        STOP
  57.     ENDIF
  58.  
  59.     CALL AE3CALC(X,AE3)
  60.  
  61.     CALL TDISSCALC (C,Mn,Cr,Mo,V,Nb,Ti,&  Al,N,AE3,TDISS)
  62.  
  63.         IF (C .LT. 0.20D0) THEN
  64.        CALL TNJCALC(X,TNJ)
  65.        IF (TNJ .NE. 0.0D0) THEN
  66.           AE4=TNJ
  67.        ENDIF
  68.     ENDIF
  69.  
  70.     IF ((C .LT. 0.10D0) .AND. (TNJ .GT. 0.0D0)) THEN
  71.        CALL TNHCALC(X,TNH)
  72.        IF (TNH .NE. 0.0D0) THEN
  73.           CALL TAHCALC (X,TAH)
  74.           Tsolid=TAH
  75.           CALL TABCALC (X,TAB)
  76.           Tliquid=TAB
  77.        ENDIF
  78.     ENDIF
  79.  
  80.     IF ((TNJ .NE. 0.0D0) .AND. (TNH .EQ. 0.0D0)) THEN
  81.        CALL THJCALC(X,THJ,XCH,XCJ)
  82.        Tsolid=THJ
  83.        CALL THBCALC(X,THB)
  84.        CALL TABCALC (X,TAB)
  85.        Tliquid=TAB
  86.     ENDIF
  87.    
  88.       IF (TNJ .EQ. 0.0D0) THEN
  89.        CALL TJECALC(X,TJE)
  90.        Tsolid=TJE
  91.        AE4=TJE
  92.  
  93.        CALL TJBCALC (X,TJB)
  94.        CALL THBCALC (X,THB)
  95.        CALL TABCALC (X,TAB)
  96.        IF ((TJB .NE. 0.0D0) .AND.  &  ((THB .NE. 0.0D0) .OR. (TAB .NE. 0.0D0))) THEN
  97.             Tliquid=TAB
  98.          ELSE
  99.           CALL TBCCALC (X,TBC)
  100.           Tliquid=TBC
  101.        ENDIF  
  102.     ENDIF
  103.  
  104.     RETURN
  105.     END
  106.  
  107.  
  108.  
  109.     SUBROUTINE XCOMP (C,MN,SI,NI,CR,MO,CU,W,V,NB,X)
  110.     IMPLICIT REAL*8 (A-H,O-Z)
  111.     REAL*8 X(0:10),AW(0:10)
  112.     REAL*8 C,Mn,Si,Ni,Cr,Mo,Cu,W,V,Nb,XTOTAL
  113.  
  114.     DATA AW /55.875D0, 12.011D0, 54.938D0, 28.0855D0, 58.7D0, 51.996D0,&  95.94D0, 63.546D0, 183.85D0, 50.9415D0, 92.9064D0/
  115.  
  116.     X(0)=(100.0D0-C-Mn-Si-Ni-Cr-Mo-Cu-W-V-Nb)/AW(0)
  117.     X(1)=C/AW(1)
  118.     X(2)=Mn/AW(2)
  119.     X(3)=Si/AW(3)
  120.     X(4)=Ni/AW(4)
  121.     X(5)=Cr/AW(5)
  122.     X(6)=Mo/AW(6)
  123.     X(7)=Cu/AW(7)
  124.     X(8)=W/AW(8)
  125.     X(9)=V/AW(9)
  126.     X(10)=Nb/AW(10)
  127.  
  128.     XTOTAL=0.0D0
  129.     DO 10 J=0,10,1
  130.        XTOTAL=XTOTAL+X(J)
  131. 10  CONTINUE
  132.     DO 20 J=0,10,1
  133.        X(J)=X(J)/XTOTAL
  134. 20  CONTINUE
  135.    
  136.     RETURN
  137.     END
  138.  
  139.  
  140.     SUBROUTINE AS1CALC (X,AS1,CEUT,CF)
  141.     IMPLICIT REAL*8 (A-H,O-Z)
  142.     REAL*8 AS1,CEUT,CF
  143.     REAL*8 X(0:10),XA(0:10),XF(0:10)
  144.     REAL*8 A0(2:10),A1(2:10),Y(0:10)
  145.     REAL*8 DGAF0,DGAF(1:10),GCA(0:10)
  146.     REAL*8 EIIA(1:10),EIIF(1:10),ECIA(1:10),ECIF(1:10)
  147.     REAL*8 F(2,2),G(2,2),FAVG,GAVG
  148.     REAL*8 AW(0:10),T0,T,DT,DTT,XC,DX,DXC,Z,DTTMAX
  149.     INTEGER INDEX,ITER,I,J,K
  150.  
  151.     DATA AW /55.875D0, 12.011D0, 54.9398D0, 28.0855D0, 58.70D0, 51.996D0,& 95.94D0, 63.546D0, 183.85D0, 50.9415D0, 92.9064D0/
  152.  
  153.     T0=1000.D0
  154.     XC=0.0364D0
  155.     INDEX=0
  156.     DX=1.0D-5
  157.     DT=0.01D0
  158.     DXC=1.0D-5
  159.     DTTMAX=10.
  160.     DXCMAX=XC/2.0D0
  161.  
  162. 10  INDEX=INDEX+1
  163.     DO 200 J=1,2,1
  164.        T=T0+(J-1.5D0)*DT
  165.     DO 200 K=1,2,1
  166.        XA(1)=XC+(K-1.5D0)*DX
  167.        RT=1.987D0*T
  168.        CALL SFECALC (T,SFE)
  169.        CALL DGAF0CALC (T,DGAF0)
  170.        CALL DGAFCALC (T,DGAF)
  171.        CALL EIIFCALC (T,SFE,EIIF)
  172.        CALL EIIACALC (T,EIIA)
  173.        CALL ECIFCALC (T,ECIF)
  174.        CALL ECIACALC (T,ECIA)
  175.        CALL GCACALC (T,GCA)
  176.        CALL A01CALC (T,A0,A1)
  177.  
  178.        XF(1)=XA(1)*EXP(DGAF(1)/RT)
  179.        XF(1)=XF(1)*EXP(EIIA(1)*XA(1)-EIIF(1)*XF(1))
  180.        EXSUM=ECIA(1)*XA(1)
  181.        DO 20 I=2,10,1
  182.           XA(I)=X(I)
  183.           XF(I)=XA(I)*EXP(DGAF(I)/RT+ECIA(I)*XA(1)-ECIF(I)*XF(1))
  184.           EXSUM=EXSUM+ECIA(I)*XA(I)
  185. 20     CONTINUE
  186.  
  187.        DO 40 I=2,10,1
  188.           IF ((I .NE. 3) .AND. (I .NE. 7)) THEN
  189.              Y(I)=XA(I)*EXP(-GCA(I)/RT-A0(I)/RT+ECIA(I)*XA(1) &   +EIIA(I)*XA(I)+DLOG(XA(1))/3.0D0+EXSUM/3.0D0)
  190.           ELSE
  191.              Y(I)=0.
  192.           ENDIF
  193. 40     CONTINUE
  194.  
  195.        DO 80 ITER=1,10,1
  196.           XF(1)=XA(1)*DEXP(DGAF(1)/RT+EIIA(1)*XA(1)-EIIF(1)*XF(1))
  197.           EXSUM=ECIA(1)*XA(1)
  198.  
  199.           DO 60 I=2,10,1
  200.              IF ((I .EQ. 3) .OR. (I .EQ. 7)) THEN
  201.                 XF(I)=X(I)*(XF(1)-0.25)/(X(1)-0.25)
  202.              ELSE
  203.                 XF(I)=X(I)+(XF(I)-0.75*Y(I))*(XF(1)-X(1))&  /(XF(1)-0.25)
  204.              ENDIF
  205.              XA(I)=XF(I)*EXP(-DGAF(I)/RT+ECIF(I)*XF(1)-ECIA(I)*XA(I) &   +EIIF(I)*XF(I)-EIIA(I)*XA(I))
  206.              XF(1)=XF(1)*EXP(EIIA(I)*XA(I)-ECIF(I)*XF(I))
  207.              EXSUM=EXSUM+ECIA(I)*XA(I)
  208. 60        CONTINUE
  209.  
  210.           DO 80 I=2,10,1
  211.              IF ((I .NE. 3) .AND. (I .NE. 7)) THEN
  212.                 Y(I)=XA(I)*EXP(-GCA(I)/RT-(A0(I)+Y(I)*A1(I))&  *(1-Y(I))*(1-Y(I))/RT+ECIA(I)*XA(1)&  +EIIA(I)*XA(I)+DLOG(XA(1))/3.0D0+EXSUM/3.0D0)
  213.              ELSE
  214.                 Y(I)=0.0D0
  215.              ENDIF
  216. 80     CONTINUE
  217.  
  218.        XA(0)=1.0D0-XA(1)
  219.        XF(0)=1.0D0-XF(1)
  220.        F(J,K)=DGAF0/RT+EIIF(1)*XF(1)*XF(1)/2-EIIA(1)*XA(1)*XA(1)/2
  221.        EXSUM=ECIA(1)*XA(1)
  222.        DO 100 I=2,10,1
  223.           XA(0)=XA(0)-XA(I)
  224.           XF(0)=XF(0)-XF(I)
  225.           F(J,K)=F(J,K)+ECIF(I)*XF(1)*XF(I)-ECIA(I)*XA(1)*XA(I)& +(EIIF(I)*XF(I)*XF(I)-EIIA(I)*XA(I)*XA(I))/2.0D0
  226.           EXSUM=EXSUM+ECIA(I)*XA(I)
  227. 100    CONTINUE
  228.        F(J,K)=F(J,K)+DLOG(XA(0)/XF(0))
  229.  
  230.        Y(0)=1.0D0
  231.        G(J,K)=GCA(0)/RT-DLOG(XA(0))+EIIA(1)*XA(1)*XA(1)/2.0D0&  -DLOG(XA(1))/3.0D0-EIIA(1)*XA(1)/3.0D0
  232.        DO 120 I=2,10,1
  233.           G(J,K)=G(J,K)+Y(I)*Y(I)*(A0(I)+Y(I)*A1(I))/RT &  +EIIA(I)*XA(I)*XA(I)/2+ECIA(I)*XA(1)*XA(I) &  -ECIA(I)*XA(I)/3.0D0
  234.           Y(0)=Y(0)-Y(I)
  235. 120    CONTINUE
  236.        G(J,K)=G(J,K)+DLOG(Y(0))
  237. 200 CONTINUE
  238.  
  239.     DFX=(F(1,2)-F(1,1)+F(2,2)-F(2,1))/DX/2
  240.     DGX=(G(1,2)-G(1,1)+G(2,2)-G(2,1))/DX/2
  241.     DFT=(F(2,1)-F(1,1)+F(2,2)-F(1,2))/DT/2
  242.     DGT=(G(2,1)-G(1,1)+G(2,2)-G(1,2))/DT/2
  243.     Z=DFX*DGT-DFT*DGX
  244.     FAVG=(F(1,1)+F(1,2)+F(2,1)+F(2,2))/4
  245.     GAVG=(G(1,1)+G(1,2)+G(2,1)+G(2,2))/4
  246.     DXC=(-FAVG*DGT+GAVG*DFT)/Z
  247.     DTT=(-GAVG*DFX+FAVG*DGX)/Z
  248.  
  249.     IF ((ABS(DXC/XC) .LE. 0.01D0) .AND. (ABS(DTT) .LE. 1.0D0)) THEN
  250.        AS1=T0+DTT-273.15D0
  251.        GOTO 400
  252.     ENDIF
  253.  
  254.     IF (INDEX .GT. 40) THEN
  255.        IF ((DTTMAX .GT. DT) .OR. (DXCMAX .GT. DX)) THEN
  256.             IF (DTTMAX .GT. DT) THEN
  257.              DTTMAX=2.0D0*DTTMAX/3.0D0
  258.           ENDIF
  259.           IF (DXCMAX .GT. DX) THEN
  260.              DXCMAX=2.0D0*DXCMAX/3.0D0
  261.           ENDIF
  262.           INDEX=0
  263.           GOTO 10
  264.        ELSE
  265.           AS1=T0-273.15D0
  266.           GOTO 400
  267.        ENDIF
  268.     ENDIF
  269.  
  270.     IF (DABS(DTT) .GT. DTTMAX) THEN
  271.        DTT=SIGN(DTTMAX,DTT)
  272.     ENDIF
  273.  
  274.     IF (DABS(DXC) .GT. DXCMAX) THEN
  275.        DXC=SIGN(DXCMAX,DXC)
  276.     ENDIF
  277.  
  278.     T0=T0+DTT
  279.     XC=XC+DXC
  280.     GOTO 10
  281.    
  282. 400     XAM=0.0D0
  283.     XFM=0.0D0
  284.     DO 500 I=0,10,1
  285.        XAM=XAM+XA(I)*AW(I)
  286.        XFM=XFM+XF(I)*AW(I)
  287. 500 CONTINUE
  288.     CEUT=100.0D0*XA(1)*AW(1)/XAM
  289.     CF=100.0D0*XF(1)*AW(1)/XFM
  290.     RETURN
  291.  
  292.     END
  293.  
  294.  
  295.     SUBROUTINE AE1CALC (X,AE1,CEUT,CF)
  296.     IMPLICIT REAL*8 (A-H,O-Z)
  297.     REAL*8 AE1,CEUT,CF
  298.     REAL*8 X(0:10),XA(0:10),XF(0:10)
  299.     REAL*8 A0(2:10), A1(2:10), A(1:10),Y(0:10)
  300.     REAL*8 DGAF0,DGAF(1:10),GCA(0:10)
  301.     REAL*8 EIIA(1:10),EIIF(1:10),ECIA(1:10),ECIF(1:10)
  302.     REAL*8 AW(0:10),F(2,2),G(2,2),FAVG,GAVG
  303.     REAL*8 T0,T,DT,DTT,XC,DX,DXC,Z,DTTMAX,DXCMAX
  304.     INTEGER INDEX,ITER,I,J,K
  305.  
  306.     DATA AW /55.875D0, 12.011D0, 54.9398D0, 28.0855D0, 58.7D0, 51.996D0,&   95.94D0, 63.546D0, 183.85D0, 50.9415D0, 92.9064D0/
  307.  
  308.     T0=1000.0D0
  309.     XC=0.0364D0
  310.     INDEX=0
  311.     DT=0.01
  312.     DX=1.0D-5
  313.     DTTMAX=10.0D0
  314.     DXCMAX=1.0D-4
  315.  
  316. 10  INDEX=INDEX+1
  317.     DO 200 J=1,2,1
  318.        T=T0+(J-1.5)*DT
  319.        RT=1.987*T
  320.     DO 200 K=1,2,1
  321.        XA(1)=XC+(K-1.5D0)*DX
  322.        CALL SFECALC (T,SFE)
  323.        CALL DGAF0CALC (T,DGAF0)
  324.        CALL DGAFCALC (T,DGAF)
  325.        CALL EIIFCALC (T,SFE,EIIF)
  326.        CALL EIIACALC (T,EIIA)
  327.        CALL ECIFCALC (T,ECIF)
  328.        CALL ECIACALC (T,ECIA)
  329.        CALL GCACALC (T,GCA)
  330.        CALL A01CALC (T,A0,A1)
  331.  
  332.        XF(1)=XA(1)*EXP(DGAF(1)/RT)
  333.        XF(1)=XF(1)*EXP(ECIA(1)*XA(1)-ECIF(1)*XF(1))
  334.          B=(XA(1)-X(1))/(XA(1)-XF(1))
  335.        DO 20 I=2,10,1
  336.           A(I)=EXP(DGAF(I)/RT+ECIA(I)*XA(1)-ECIF(I)*XF(1))
  337.           XA(I)=X(I)/(1-B*(1-A(I)))
  338.           XF(I)=A(I)*XA(I)
  339.           XF(1)=XF(1)*EXP(ECIA(I)*XA(I)-ECIF(I)*XF(I))
  340. 20     CONTINUE
  341.  
  342.        DO 40 ITER=1,10,1
  343.           XF(1)=XA(1)*EXP(DGAF(1)/RT+ECIA(1)*XA(1)-ECIF(1)*XF(1))
  344.             B=(XA(1)-X(1))/(XA(1)-XF(1))
  345.        DO 40 I=2,10,1
  346.           A(I)=EXP(DGAF(I)/RT+ECIA(I)*XA(1)-ECIF(I)*XF(1)& +EIIA(I)*XA(I)-EIIF(I)*XF(I))
  347.           XA(I)=X(I)/(1-B*(1-A(I)))
  348.           XF(I)=A(I)*XA(I)
  349.           XF(1)=XF(1)*EXP(ECIA(I)*XA(I)-ECIF(I)*XF(I))
  350. 40     CONTINUE
  351.  
  352.        XA(0)=1.0-XA(1)
  353.        XF(0)=1.0-XF(1)
  354.        F(J,K)=DGAF0/RT+EIIF(1)*XF(1)*XF(1)/2-EIIA(1)*XA(1)*XA(1)/2
  355.        EXSUM=ECIA(1)*XA(1)
  356.        DO 60 I=2,10,1
  357.           XA(0)=XA(0)-XA(I)
  358.           XF(0)=XF(0)-XF(I)
  359.           F(J,K)=F(J,K)+ECIF(I)*XF(1)*XF(I)-ECIA(I)*XA(1)*XA(I)&   +(EIIF(I)*XF(I)*XF(I)-EIIA(I)*XA(I)*XA(I))/2
  360.           EXSUM=EXSUM+ECIA(I)*XA(I)
  361. 60     CONTINUE
  362.        F(J,K)=F(J,K)+DLOG(XA(0)/XF(0))
  363.  
  364.        Y(0)=1.0
  365.        Y(1)=1.0
  366.        DO 80 I=2,10,1
  367.           IF ((I .NE. 3) .AND. (I .NE. 7)) THEN
  368.              Y(I)=XA(I)*DEXP(-GCA(I)/RT-A0(I)/RT+ECIA(I)*XA(1)&    +EIIA(I)*XA(I)+DLOG(XA(1))/3.0D0+EXSUM/3.0D0)
  369.           ELSE
  370.              Y(I)=0.
  371.           ENDIF
  372.           Y(0)=Y(0)-Y(I)
  373. 80     CONTINUE
  374.  
  375.        DO 100 ITER=1,10,1
  376.        DO 100 I=2,10,1
  377.           IF ((I .NE. 3) .AND. (I .NE. 7)) THEN
  378.              Y(I)=XA(I)*DEXP(-GCA(I)/RT-(A0(I)+Y(I)*A1(I))&  *(1-Y(I))*(1-Y(I))/RT+ECIA(I)*XA(1)&  +EIIA(I)*XA(I)+DLOG(XA(1))/3.0D0+EXSUM/3.0D0)
  379.           ELSE
  380.              Y(I)=0.0
  381.           ENDIF
  382. 100    CONTINUE
  383.  
  384.        G(J,K)=GCA(0)/RT-DLOG(XA(0))+EIIA(1)*XA(1)*XA(1)/2.0D0&   -DLOG(XA(1))/3.0D0-ECIA(1)*XA(1)/3.0D0
  385.        DO 120 I=2,10,1
  386.           G(J,K)=G(J,K)+Y(I)*Y(I)*(A0(I)+Y(I)*A1(I))/RT&   +EIIA(I)*XA(I)*XA(I)/2+ECIA(I)*XA(1)*XA(I) &  -ECIA(I)*XA(I)/3.0D0
  387.           Y(0)=Y(0)-Y(I)
  388. 120    CONTINUE
  389.        G(J,K)=G(J,K)+DLOG(Y(0))
  390.  
  391. 200 CONTINUE
  392.  
  393.     DFX=(F(1,2)-F(1,1)+F(2,2)-F(2,1))/DX/2
  394.     DGX=(G(1,2)-G(1,1)+G(2,2)-G(2,1))/DX/2
  395.     DFT=(F(2,1)-F(1,1)+F(2,2)-F(1,2))/DT/2
  396.     DGT=(G(2,1)-G(1,1)+G(2,2)-G(1,2))/DT/2
  397.     Z=DFX*DGT-DFT*DGX
  398.     FAVG=(F(1,1)+F(1,2)+F(2,1)+F(2,2))/4
  399.     GAVG=(G(1,1)+G(1,2)+G(2,1)+G(2,2))/4
  400.     DXC=(-FAVG*DGT+FAVG*DFT)/Z
  401.     DTT=(-GAVG*DFX+FAVG*DGX)/Z
  402.    
  403.     IF ((DABS(DXC/XC) .LE. 0.01) .AND. (DABS(DTT) .LE. 1.0D0)) THEN
  404.        AE1=T0+DTT-273.15D0
  405.        GOTO 400
  406.     ENDIF
  407.  
  408.     IF (INDEX .GT. 40) THEN
  409.        IF ((DTTMAX .GT. DT) .OR. (DXCMAX .GT. DX)) THEN
  410.             IF (DTTMAX .GT. DT) THEN
  411.              DTTMAX=2*DTTMAX/3.0D0
  412.           ENDIF
  413.           IF (DXCMAX .GT. DX) THEN
  414.              DXCMAX=2*DXCMAX/3.0D0
  415.           ENDIF
  416.           INDEX=0
  417.           GOTO 10
  418.        ELSE
  419.           AE1=T0-273.15D0
  420.           GOTO 400
  421.        ENDIF
  422.     ENDIF
  423.  
  424.     IF (DABS(DTT) .GT. DTTMAX) THEN
  425.        DTT=SIGN(DTTMAX,DTT)
  426.     ENDIF
  427.  
  428.     IF (DABS(DXC) .GT. DXCMAX) THEN
  429.        DXC=SIGN(DXCMAX,DXC)
  430.     ENDIF
  431.  
  432.     T0=T0+DTT
  433.     XC=XC+DXC
  434.     GOTO 10
  435.  
  436. 400     XAM=0.0D0
  437.     XFM=0.0D0
  438.     DO 500 I=0,10,1
  439.        XAM=XAM+XA(I)*AW(I)
  440.        XFM=XFM+XF(I)*AW(I)
  441. 500 CONTINUE
  442.     CEUT=100.0*XA(1)*AW(1)/XAM
  443.     CF=100.0*XF(1)*AW(1)/XFM
  444.     RETURN
  445.  
  446.     END
  447.  
  448.  
  449.     SUBROUTINE AE3CALC (XA,AE3)
  450.     IMPLICIT REAL*8 (A-H,O-Z)
  451.     REAL*8 XA(0:10),XF(0:10),AE3
  452.     REAL*8 DGAF0,DGAF(1:10)
  453.     REAL*8 EIIF(1:10),EIIA(1:10),ECIF(1:10),ECIA(1:10)
  454.     REAL*8 F(1:2),FAVG
  455.     REAL*8 T0,DT,T,RT,DTT,DTTMAX
  456.     INTEGER INDEX,ITER,I
  457.  
  458.     T0=1127.5D0-4805*XA(1)
  459.     DT=0.01D0
  460.     DTTMAX=10.0D0
  461.     INDEX=0
  462.  
  463. 10  INDEX=INDEX+1
  464.     DO 100 K=1,2,1
  465.     T=T0+(K-1.5D0)*DT
  466.     RT=1.987D0*T
  467.     CALL SFECALC (T,SFE)
  468.     CALL DGAF0CALC (T,DGAF0)
  469.     CALL DGAFCALC (T,DGAF)
  470.     CALL EIIFCALC (T,SFE,EIIF)
  471.     CALL EIIACALC (T,EIIA)
  472.     CALL ECIFCALC (T,ECIF)
  473.     CALL ECIACALC (T,ECIA)
  474.  
  475.     XF(1)=XA(1)*DEXP(DGAF(1)/RT)
  476.       XF(1)=XF(1)*DEXP(ECIA(1)*XA(1)-ECIF(1)*XF(1))
  477.     DO 20 I=2,10,1
  478.        XF(I)=XA(I)*EXP(DGAF(I)/RT+ECIA(I)*XA(1)-ECIF(I)*XF(1))
  479.            XF(1)=XF(1)*EXP(ECIA(I)*XA(I)-ECIF(I)*XF(I))
  480. 20  CONTINUE
  481.  
  482.     DO 40 ITER=1,10,1
  483.        XF(1)=XA(1)*DEXP(DGAF(1)/RT+ECIA(1)*XA(1)-ECIF(1)*XF(1))
  484.     DO 40 I=2,10,1
  485.        XF(I)=XA(I)*DEXP(DGAF(I)/RT+ECIA(I)*XA(1)-ECIF(I)*XF(1)  &    +EIIA(I)*XA(I)-EIIF(I)*XF(I))
  486.        XF(3)=XA(3)*DEXP(DGAF(3)/RT+ECIA(3)*XA(1)-ECIF(3)*XF(1)  &    +EIIA(3)*XA(3)-EIIF(3)*XF(3))
  487.        XF(1)=XF(1)*DEXP(ECIA(I)*XA(I)-ECIF(I)*XF(I))
  488. 40  CONTINUE
  489.  
  490.     XF(0)=1.0
  491.     DO 60 I=1,10,1
  492.        XF(0)=XF(0)-XF(I)
  493. 60  CONTINUE
  494.  
  495.     F(K)=DGAF0/RT+DLOG(XA(0)/XF(0))+EIIF(1)*XF(1)*XF(1)/2&  -EIIA(1)*XA(1)*XA(1)/2
  496.     DO 100 I=2,10,1
  497.        F(K)=F(K)+ECIF(I)*XF(1)*XF(I)-ECIA(I)*XA(1)*XA(I) &    +EIIF(I)*XF(I)*XF(I)/2-EIIA(I)*XA(I)*XA(I)/2
  498. 100 CONTINUE
  499.  
  500.     FAVG=(F(1)+F(2))/2
  501.     DTT=-DT*FAVG/(F(2)-F(1))
  502.  
  503.     IF (DABS(DTT) .LE. 1.0D0) THEN
  504.        AE3=T0+DTT-273.15D0
  505.        RETURN
  506.     ENDIF
  507.  
  508.     IF (INDEX .GE. 40) THEN
  509.        WRITE(*,*) '****** WARNING ******'
  510.        WRITE(*,*) 'AE3 COMPUTATION FAILS TO CONVERGE',&  ' IN 40 ITERATIONS!'
  511.        STOP
  512.     ELSEIF (DABS(DTT) .GT. DTTMAX) THEN
  513.        DTT=SIGN(DTTMAX,DTT)
  514.        T0=T0+DTT
  515.        GOTO 10
  516.     ELSE
  517.        T0=T0+DTT
  518.        GOTO 10
  519.     ENDIF
  520.  
  521.     END
  522.  
  523.  
  524.     SUBROUTINE TNJCALC (XA,TNJ)
  525.     IMPLICIT REAL*8 (A-H,O-Z)
  526.     REAL*8 XA(0:10),XF(0:10),TNJ
  527.     REAL*8 DGAF0,DGAF(1:10),F(1:2),FAVG
  528.     REAL*8 EIIF(1:10),EIIA(1:10),ECIF(1:10),ECIA(1:10)
  529.     REAL*8 T0,DT,T,RT,DTT,DTTMAX
  530.  
  531.     T0=1673.0D0+10555.0D0*XA(1)
  532.     DT=0.01
  533.     DTTMAX=10.
  534.     INDEX=0
  535.  
  536. 10  INDEX=INDEX+1
  537.     DO 100 K=1,2,1
  538.     T=T0+(K-1.5D0)*DT
  539.     RT=1.987D0*T
  540.     CALL SFECALC (T,SFE)
  541.     CALL DGAF0CALC (T,DGAF0)
  542.     CALL DGAFCALC (T,DGAF)
  543.     CALL EIIFCALC (T,SFE,EIIF)
  544.     CALL EIIACALC (T,EIIA)
  545.     CALL ECIFCALC (T,ECIF)
  546.     CALL ECIACALC (T,ECIA)
  547.  
  548.       XF(1)=XA(1)*DEXP(DGAF(1)/RT)
  549.     XF(1)=XF(1)*DEXP(EIIA(1)*XA(1)-EIIF(1)*XF(1))
  550.     DO 20 I=2,10,1
  551.        XF(I)=XA(I)*DEXP(DGAF(I)/RT+ECIA(I)*XA(1)-ECIF(I)*XF(1))
  552.        XF(I)=XF(I)*DEXP(EIIA(I)*XA(I)-EIIF(I)*XF(I))
  553.        XF(1)=XF(1)*DEXP(ECIA(I)*XA(I)-ECIF(I)*XF(I))
  554. 20  CONTINUE
  555.  
  556.     DO 40 ITER=1,10,1
  557.        XF(1)=XA(1)*DEXP(DGAF(1)/RT+EIIA(1)*XA(1)-EIIF(1)*XF(1))
  558.        DO 40 I=2,10,1
  559.           XF(I)=XA(I)*DEXP(DGAF(I)/RT+ECIA(I)*XA(1)-ECIF(I)*XF(1)   &     +EIIA(I)*XA(I)-EIIF(I)*XF(I))
  560.           XF(1)=XF(1)*DEXP(ECIA(I)*XA(I)-ECIF(I)*XF(I))
  561. 40  CONTINUE
  562.  
  563.     XF(0)=1.0
  564.     DO 60 I=1,10,1
  565.        XF(0)=XF(0)-XF(I)
  566. 60  CONTINUE
  567.  
  568.     F(K)=DGAF0/RT+DLOG(XA(0)/XF(0))&     +(EIIF(1)*XF(1)*XF(1)-EIIA(1)*XA(1)*XA(1))/2
  569.     DO 80 I=2,10,1
  570.        F(K)=F(K)+ECIF(I)*XF(1)*XF(I)-ECIA(I)*XA(1)*XA(I) &        +(EIIF(I)*XF(I)*XF(I)-EIIA(I)*XA(I)*XA(I))/2
  571. 80  CONTINUE
  572. 100 CONTINUE
  573.  
  574.     FAVG=(F(1)+F(2))/2
  575.     DTT=-DT*FAVG/(F(2)-F(1))
  576.  
  577.     IF (DABS(DTT) .LE. 1.0) THEN
  578.        TNJ=T0+DTT-273.15D0
  579.        RETURN
  580.     ELSEIF ((T0+DTT) .GT. 1820.0D0) THEN
  581.        TNJ=0.0D0
  582.        RETURN
  583.     ENDIF
  584.  
  585.     IF (INDEX .GE. 40) THEN
  586.        WRITE(*,*) '****** WARNING ******'
  587.        WRITE(*,*) 'TNJ COMPUTATION FAILS TO CONVERGE',&                ' IN 40 ITERATIONS!'
  588.        STOP
  589.     ELSEIF (DABS(DTT) .GT. DTTMAX) THEN
  590.        DTT=SIGN(DTTMAX,DTT)
  591.        T0=T0+DTT
  592.        GOTO 10
  593.     ELSE
  594.        T0=T0+DTT
  595.        GOTO 10
  596.     ENDIF
  597.  
  598. 1000    RETURN
  599.     END
  600.  
  601.  
  602.  
  603.     SUBROUTINE TNHCALC (XF,TNH)
  604.     IMPLICIT REAL*8 (A-H,O-Z)
  605.     REAL*8 XA(0:10),XF(0:10),TNH
  606.     REAL*8 DGAF0,DGAF(1:10),F(1:2),FAVG
  607.     REAL*8 EIIF(1:10),EIIA(1:10),ECIF(1:10),ECIA(1:10)
  608.     REAL*8 T0,DT,T,RT,DTT,DTTMAX
  609.  
  610.     T0=1673.0D0+23750.0D0*XF(1)
  611.     DT=0.01D0
  612.     DTTMAX=10.0D0
  613.     INDEX=0
  614.  
  615. 10  INDEX=INDEX+1
  616.     DO 100 K=1,2,1
  617.     T=T0+(K-1.5D0)*DT
  618.     RT=1.987D0*T
  619.     CALL SFECALC (T,SFE)
  620.     CALL DGAF0CALC (T,DGAF0)
  621.     CALL DGAFCALC (T,DGAF)
  622.     CALL EIIFCALC (T,SFE,EIIF)
  623.     CALL EIIACALC (T,EIIA)
  624.     CALL ECIFCALC (T,ECIF)
  625.     CALL ECIACALC (T,ECIA)
  626.  
  627.       XA(1)=XF(1)*DEXP(-DGAF(1)/RT)
  628.     XA(1)=XA(1)*DEXP(ECIF(1)*XF(1)-ECIA(1)*XA(1))
  629.     DO 20 I=2,10,1
  630.        XA(I)=XF(I)*DEXP(-DGAF(I)/RT+ECIF(I)*XF(1)-ECIA(I)*XA(1))
  631.        XA(I)=XA(I)*DEXP(EIIF(I)*XF(I)-EIIA(I)*XA(I))
  632.        XA(1)=XA(1)*DEXP(ECIF(I)*XF(I)-ECIA(I)*XA(I))
  633. 20  CONTINUE
  634.  
  635.     DO 40 ITER=1,10,1
  636.        XA(1)=XF(1)*DEXP(-DGAF(1)/RT+ECIF(1)*XF(1)-ECIA(1)*XA(1))
  637.     DO 40 I=2,10,1
  638.        XA(I)=XF(I)*DEXP(-DGAF(I)/RT+ECIF(I)*XF(1)+EIIF(I)*XF(I)&                     -EIIA(I)*XA(I)-ECIA(I)*XA(1))
  639.        XA(1)=XA(1)*DEXP(ECIF(I)*XF(I)-ECIA(I)*XA(I))
  640. 40  CONTINUE
  641.  
  642.     XA(0)=1.0
  643.     DO 60 I=1,10,1
  644.        XA(0)=XA(0)-XA(I)
  645. 60  CONTINUE
  646.  
  647.     F(K)=DGAF0/RT+DLOG(XA(0)/XF(0))+EIIF(1)*XF(1)*XF(1)/2&     -EIIA(1)*XA(1)*XA(1)/2
  648.     DO 100 I=2,10,1
  649.        F(K)=F(K)+ECIF(I)*XF(1)*XF(I)-ECIA(I)*XA(1)*XA(I)&        +EIIF(I)*XF(I)*XF(I)/2-EIIA(I)*XA(I)*XA(I)/2
  650. 100 CONTINUE
  651.  
  652.     FAVG=(F(1)+F(2))/2
  653.     DTT=-DT*FAVG/(F(2)-F(1))
  654.  
  655.     IF (DABS(DTT) .LE. 1.0) THEN
  656.        TNH=T0+DTT-273.15D0
  657.        RETURN
  658.     ELSEIF ((T0+DTT) .GT. 1820.0D0) THEN
  659.        TNH=0.0D0
  660.        RETURN
  661.     ENDIF
  662.  
  663.     IF (INDEX .GE. 40) THEN
  664.        WRITE(*,*) '****** WARNING ******'
  665.        WRITE(*,*) 'TNH COMPUTATION FAILS TO CONVERGE',&                ' IN 40 ITERATIONS!'
  666.          TNH=0.0D0
  667.        RETURN
  668.     ELSEIF (DABS(DTT) .GT. DTTMAX) THEN
  669.        DTT=SIGN(DTTMAX,DTT)
  670.        T0=T0+DTT
  671.        GOTO 10
  672.     ELSE
  673.        T0=T0+DTT
  674.        GOTO 10
  675.     ENDIF
  676.  
  677. 1000    RETURN
  678.     END
  679.  
  680.  
  681.     SUBROUTINE TAHCALC (XF,TAH)
  682.     IMPLICIT REAL*8 (A-H,O-Z)
  683.     REAL*8 XL(0:10),XF(0:10),TAH
  684.     REAL*8 DGAF0,DGLF0,DGLF(1:10),F(1:2),FAVG
  685.     REAL*8 EIIF(1:10),EIIL(1:10),ECIF(1:10),ECIL(1:10)
  686.     REAL*8 T0,DT,T,RT,DTT,DTTMAX
  687.     INTEGER INDEX,ITER,I
  688.  
  689.     T0=1811.0D0-10750.0D0*XF(1)
  690.     DT=0.01D0
  691.     DTTMAX=10.0D0
  692.     INDEX=0
  693.  
  694. 10  INDEX=INDEX+1
  695.     DO 100 K=1,2,1
  696.     T=T0+(K-1.5D0)*DT
  697.     RT=1.987D0*T
  698.     CALL SFECALC (T,SFE)
  699.     CALL DGAF0CALC (T,DGAF0)
  700.     CALL DGLFCALC (T,DGLF)
  701.     CALL EIIFCALC (T,SFE,EIIF)
  702.     CALL EIILCALC (T,EIIL)
  703.     CALL ECIFCALC (T,ECIF)
  704.     CALL ECILCALC (T,ECIL)
  705.        
  706.     DGLF0=-2646.0D0+39.1693D0*T-5.27D0*T*DLOG(T)+0.001D0*T*T+DGAF0
  707.  
  708.     XL(1)=XF(1)*DEXP(-DGLF(1)/RT)
  709.       XL(1)=XL(1)*DEXP(ECIF(1)*XF(1)-ECIL(1)*XL(1))
  710.     DO 20 I=2,10,1
  711.        XL(I)=XF(I)*DEXP(-DGLF(I)/RT+ECIF(I)*XF(1)-ECIL(I)*XL(1))
  712.        XL(I)=XL(I)*DEXP(EIIF(I)*XF(I)-EIIL(I)*XL(I))
  713.        XL(1)=XL(1)*DEXP(ECIF(I)*XF(I)-ECIF(I)*XL(I))
  714. 20  CONTINUE
  715.  
  716.     DO 40 ITER=1,10,1
  717.        XL(1)=XF(1)*DEXP(-DGLF(1)/RT+ECIF(1)*XF(1)-ECIL(1)*XL(1))
  718.     DO 40 I=2,10,1
  719.        XL(I)=XF(I)*DEXP(-DGLF(I)/RT+ECIF(I)*XF(1)-ECIL(I)*XL(1) &     +EIIF(I)*XF(I)-EIIL(I)*XL(I))
  720.        XL(1)=XL(1)*DEXP(ECIF(I)*XF(I)-ECIL(I)*XL(I))
  721. 40  CONTINUE
  722.  
  723.     XL(0)=1.0D0
  724.     DO 60 I=1,10,1
  725.        XL(0)=XL(0)-XL(I)
  726. 60  CONTINUE
  727.  
  728.     F(K)=DGLF0/RT+DLOG(XL(0)/XF(0))+EIIF(1)*XF(1)*XF(1)/2 &       -EIIL(1)*XL(1)*XL(1)/2
  729.     DO 100 I=2,10,1
  730.        F(K)=F(K)+ECIF(I)*XF(1)*XF(I)-ECIL(I)*XL(1)*XL(I)&        +EIIF(I)*XF(I)*XF(I)/2-EIIL(I)*XL(I)*XL(I)/2
  731. 100 CONTINUE
  732.  
  733.     FAVG=(F(1)+F(2))/2
  734.     DTT=-DT*FAVG/(F(2)-F(1))
  735.  
  736.     IF (DABS(DTT) .LE. 1.0) THEN
  737.        TAH=T0+DTT-273.15D0
  738.        RETURN
  739.     ELSEIF ((T0+DTT) .GT. 1820.0D0) THEN
  740.        TAH=0.0D0
  741.        RETURN
  742.     ENDIF
  743.  
  744.     IF (INDEX .GE. 40) THEN
  745.        WRITE(*,*) '****** WARNING ******'
  746.        WRITE(*,*) 'TAH COMPUTATION FAILS TO CONVERGE',&                ' IN 40 ITERATIONS!'
  747.        STOP
  748.     ELSEIF (DABS(DTT) .GT. DTTMAX) THEN
  749.        DTT=SIGN(DTTMAX,DTT)
  750.        T0=T0+DTT
  751.        GOTO 10
  752.     ELSE
  753.        T0=T0+DTT
  754.        GOTO 10
  755.     ENDIF
  756.  
  757.     END
  758.  
  759.  
  760.  
  761.     SUBROUTINE THJCALC(X,THJ,XCH,XCJ)
  762.     IMPLICIT REAL*8 (A-H,O-Z)
  763.     REAL*8 X(0:10),XA(0:10),XF(0:10),XL(0:10),THJ
  764.     REAL*8 DGAF0,DGFA(1:10),DGLA0,DGLA(1:10)
  765.     REAL*8 EIIL(1:10),EIIF(1:10),EIIA(1:10)
  766.     REAL*8 ECIL(1:10),ECIF(1:10),ECIA(1:10)
  767.     REAL*8 A(1:10),B,F(2,2),G(2,2)
  768.     REAL*8 T0,T,DT,DTT,XC,DX,DXC,Z,DTTMAX
  769.     INTEGER INDEX
  770.  
  771.     T0=1765.0D0
  772.     XC=0.008D0
  773.     DTTMAX=10.0D0
  774.     INDEX=0
  775.     DX=1.0D-5
  776.     DT=0.01D0
  777.     DXC=1.0D-5
  778.     DTTMAX=10
  779.     DXCMAX=XC/2.0D0
  780.  
  781. 10  INDEX=INDEX+1
  782.     DO 100 J=1,2,1
  783.        T=T0+(J-1.5D0)*DT
  784.        RT=1.987D0*T
  785.     DO 100 K=1,2,1
  786.        XA(1)=XC+(K-1.5D0)*DX
  787.        CALL DGAF0CALC(T,DGAF0)
  788.        CALL DGFACALC(T,DGFA)
  789.        CALL DGLACALC(T,DGLA)
  790.        CALL EIIACALC(T,EIIA)
  791.        CALL EIIFCALC(T,SFE,EIIF)
  792.        CALL EIILCALC(T,EIIL)
  793.          CALL ECIACALC(T,ECIA)
  794.        CALL ECIFCALC(T,ECIF)
  795.        CALL ECILCALC(T,ECIL)
  796.  
  797.        DGLA0=-2646.0D0+39.1693D0*T-5.27D0*T*DLOG(T)+0.001D0*T*T
  798.  
  799.          XF(1)=XA(1)*DEXP(-DGFA(1)/RT)
  800.        XF(1)=XF(1)*DEXP(EIIA(1)*XA(1)-EIIF(1)*XF(1))
  801.        XL(1)=XA(1)*DEXP(-DGLA(1)/RT)
  802.        XL(1)=XL(1)*DEXP(EIIA(1)*XA(1)-EIIL(1)*XL(1))
  803.        B=(XA(1)-X(1))/(XA(1)-XF(1))
  804.        DO 20 I=2,10,1
  805.           A(I)=DEXP(-DGFA(I)/RT+ECIA(I)*XA(1)-ECIF(I)*XF(1))
  806.           XA(I)=X(I)/(1-B*(1-A(I)))
  807.           XF(I)=XA(I)*A(I)
  808.           XF(1)=XF(1)*DEXP(ECIA(I)*XA(I)-ECIF(I)*XF(I))
  809.           XL(I)=XA(I)*DEXP(-DGLA(I)/RT+ECIA(I)*XA(1)-ECIL(I)*XL(1))
  810.           XL(1)=XL(1)*DEXP(ECIA(I)*XA(I)-ECIL(I)*XL(I))
  811. 20     CONTINUE
  812.  
  813.        DO 40 ITER=1,10,1
  814.           B=(XA(1)-X(1))/(XA(1)-XF(1))
  815.           XF(1)=XA(1)*DEXP(-DGFA(1)/RT+EIIA(1)*XA(1)-EIIF(1)*XF(1))
  816.           XL(1)=XA(1)*DEXP(-DGLA(1)/RT+EIIA(1)*XA(1)-EIIL(1)*XL(1))
  817.        DO 40 I=2,10,1
  818.           A(I)=DEXP(-DGFA(I)/RT+ECIA(I)*XA(1)-ECIF(I)*XF(1) &               +EIIA(I)*XA(I)-EIIF(I)*XF(I))
  819.           XA(I)=X(I)/(1-B*(1-A(I)))
  820.           XF(I)=XA(I)*A(I)
  821.           XF(1)=XF(1)*DEXP(ECIA(I)*XA(I)-ECIF(I)*XF(I))
  822.           XL(I)=XA(I)*DEXP(-DGLA(I)/RT+ECIA(I)*XA(1)-ECIL(I)*XL(1)&                      +EIIA(I)*XA(I)-EIIL(I)*XL(I))
  823.           XL(1)=XL(1)*DEXP(ECIA(I)*XA(I)-ECIL(I)*XL(I))
  824. 40     CONTINUE
  825.              
  826.          XA(0)=1.0D0
  827.        XF(0)=1.0D0
  828.        XL(0)=1.0D0
  829.        DO 60 I=1,10,1
  830.           XA(0)=XA(0)-XA(I)
  831.           XF(0)=XF(0)-XF(I)
  832.           XL(0)=XL(0)-XL(I)
  833. 60     CONTINUE
  834.  
  835.        F(J,K)=DGAF0/RT+DLOG(XA(0)/XF(0)) &          +(EIIF(1)*XF(1)*XF(1)-EIIA(1)*XA(1)*XA(1))/2
  836.        G(J,K)=DGLA0/RT+DLOG(XL(0)/XA(0))&          +(EIIA(1)*XA(1)*XA(1)-EIIL(1)*XL(1)*XL(1))/2
  837.        DO 80 I=2,10,1
  838.           F(J,K)=F(J,K)+(EIIF(I)*XF(I)*XF(I)-EIIA(I)*XA(I)*XA(I))/2&             +ECIF(I)*XF(1)*XF(I)-ECIA(I)*XA(1)*XA(I)
  839.           G(J,K)=G(J,K)+(EIIA(I)*XA(I)*XA(I)-EIIL(I)*XL(I)*XL(I))/2 &             +ECIA(I)*XA(1)*XA(I)-ECIL(I)*XL(1)*XL(I)
  840. 80     CONTINUE
  841. 100 CONTINUE
  842.  
  843.     DFX=(F(1,2)-F(1,1)+F(2,2)-F(2,1))/DX/2
  844.     DGX=(G(1,2)-G(1,1)+G(2,2)-G(2,1))/DX/2
  845.     DFT=(F(2,1)-F(1,1)+F(2,2)-F(1,2))/DT/2
  846.     DGT=(G(2,1)-G(1,1)+G(2,2)-G(1,2))/DT/2
  847.     Z=DFX*DGT-DFT*DGX
  848.     FAVG=(F(1,1)+F(1,2)+F(2,1)+F(2,2))/4
  849.     GAVG=(G(1,1)+G(1,2)+G(2,1)+G(2,2))/4
  850.     DXC=(-FAVG*DGT+GAVG*DFT)/Z
  851.     DTT=(-GAVG*DFX+FAVG*DGX)/Z
  852.  
  853.     IF ((DABS(DXC/XC) .LE. 0.01D0) .AND. (DABS(DTT) .LE. 1.0D0)) THEN
  854.        THJ=T0+DTT-273.15D0
  855.        XCH=XF(1)
  856.        XCJ=XA(1)
  857.        RETURN
  858.     ENDIF
  859.  
  860.     IF (INDEX .GT. 40) THEN
  861.        THJ=0.0
  862.        RETURN
  863.     ENDIF
  864.  
  865.     IF (DABS(DTT) .GT. DTTMAX) THEN
  866.        DTT=SIGN(DTTMAX,DTT)
  867.     ENDIF
  868.  
  869.     IF (DABS(DXC) .GT. DXCMAX) THEN
  870.        DXC=SIGN(DXCMAX,DXC)
  871.     ENDIF
  872.  
  873.     T0=T0+DTT
  874.     XC=XC+DXC
  875.     GOTO 10
  876.    
  877.     END
  878.  
  879.  
  880.     SUBROUTINE THBCALC(X,THB)
  881.     IMPLICIT REAL*8 (A-H,O-Z)
  882.     REAL*8 X(0:10),XA(0:10),XF(0:10),XL(0:10),THB
  883.     REAL*8 DGAF0,DGFA0,DGFA(1:10),DGLA0,DGLA(1:10),DGLF(1:10)
  884.     REAL*8 EIIL(1:10),EIIF(1:10),EIIA(1:10)
  885.     REAL*8 ECIL(1:10),ECIF(1:10),ECIA(1:10)
  886.     REAL*8 A(1:10),B,F(2,2),G(2,2)
  887.     REAL*8 T0,T,DT,DTT,XC,DX,DXC,Z,DTTMAX
  888.     INTEGER INDEX
  889.  
  890.     T0=1765.0D0
  891.     XC=0.008D0
  892.     DTTMAX=10.0D0
  893.     INDEX=0
  894.     DX=1.0D-5
  895.     DT=0.01D0
  896.     DXC=1.0D-5
  897.     DTTMAX=10
  898.     DXCMAX=XC/2.0D0
  899.  
  900. 10  INDEX=INDEX+1
  901.     DO 100 J=1,2,1
  902.        T=T0+(J-1.5D0)*DT
  903.     DO 100  K=1,2,1
  904.        XA(1)=XC+(K-1.5D0)*DX
  905.        RT=1.987D0*T
  906.        CALL DGAF0CALC(T,DGAF0)
  907.        CALL DGFACALC(T,DGFA)
  908.        CALL DGLACALC(T,DGLA)
  909.        CALL DGLFCALC(T,DGLF)
  910.        CALL SFECALC(T,SFE)
  911.        CALL EIIACALC(T,EIIA)
  912.        CALL EIIFCALC(T,SFE,EIIF)
  913.        CALL EIILCALC(T,EIIL)
  914.          CALL ECIACALC(T,ECIA)
  915.        CALL ECIFCALC(T,ECIF)
  916.        CALL ECILCALC(T,ECIL)
  917.  
  918.        DGFA0=-DGAF0
  919.        DGLA0=-2646.0D0+39.1693D0*T-5.27D0*T*DLOG(T)+0.001D0*T*T
  920.  
  921.        XF(1)=XA(1)*DEXP(-DGFA(1)/RT)
  922.        XF(1)=XF(1)*DEXP(EIIA(1)*XA(1)-EIIF(1)*XF(1))
  923.        XL(1)=XA(1)*DEXP(-DGLA(1)/RT)
  924.        XL(1)=XL(1)*DEXP(EIIA(1)*XA(1)-EIIL(1)*XL(1))
  925.        B=(XL(1)-X(1))/(XL(1)-XF(1))
  926.        DO 20 I=2,10,1
  927.           A(I)=DEXP(DGLF(I)/RT+ECIL(I)*XL(1)-ECIL(I)*XF(1))
  928.           XL(I)=X(I)/(1-B*(1-A(I)))
  929.           XF(I)=XL(I)*A(I)
  930.           XA(I)=XF(I)*DEXP(DGFA(I)/RT+ECIF(I)*XF(1)-ECIA(I)*XA(1))
  931.           XF(1)=XF(1)*DEXP(ECIA(I)*XA(I)-ECIF(I)*XF(I))
  932.           XL(1)=XL(1)*DEXP(ECIA(I)*XA(I)-ECIL(I)*XL(I))
  933. 20     CONTINUE
  934.                  
  935.          XA(0)=1.0D0-XA(1)
  936.        XF(0)=1.0D0-XF(1)
  937.        XL(0)=1.0D0-XL(1)
  938.        DO 60 I=1,10,1
  939.           XA(0)=XA(0)-XA(I)
  940.           XF(0)=XF(0)-XF(I)
  941.           XL(0)=XL(0)-XL(I)
  942. 60     CONTINUE
  943.          
  944.        F(J,K)=DGFA0/RT+DLOG(XF(0)/XA(0))&          +(EIIA(1)*XA(1)*XA(1)-EIIF(1)*XF(1)*XF(1))/2
  945.        G(J,K)=DGLA0/RT+DLOG(XL(0)/XA(0))&          +(EIIA(1)*XA(1)*XA(1)-EIIL(1)*XL(1)*XL(1))/2
  946.        DO 80 I=2,10,1
  947.           F(J,K)=F(J,K)+(EIIA(I)*XA(I)*XA(I)-EIIF(I)*XF(I)*XF(I))/2. &             +ECIA(I)*XA(1)*XA(I)-ECIF(I)*XF(1)*XF(I)
  948.           G(J,K)=G(J,K)+(EIIA(I)*XA(I)*XA(I)-EIIL(I)*XL(I)*XL(I))/2.&             +ECIA(I)*XA(1)*XA(I)-ECIL(I)*XL(1)*XL(I)
  949.    
  950.  
  951. 80  CONTINUE
  952. 100     CONTINUE
  953.     DFX=(F(1,2)-F(1,1)+F(2,2)-F(2,1))/DX/2
  954.     DGX=(G(1,2)-G(1,1)+G(2,2)-G(2,1))/DX/2
  955.     DFT=(F(2,1)-F(1,1)+F(2,2)-F(1,2))/DT/2
  956.     DGT=(G(2,1)-G(1,1)+G(2,2)-G(1,2))/DT/2
  957.     Z=DFX*DGT-DFT*DGX
  958.     FAVG=(F(1,1)+F(1,2)+F(2,1)+F(2,2))/4
  959.     GAVG=(G(1,1)+G(1,2)+G(2,1)+G(2,2))/4
  960.     DXC=(-FAVG*DGT+GAVG*DFT)/Z
  961.     DTT=(-GAVG*DFX+FAVG*DGX)/Z
  962.  
  963.     IF ((ABS(DXC/XC) .LE. 0.01) .AND. (ABS(DTT) .LE. 1.)) THEN
  964.        THB=T0+DTT-273.15D0
  965.        IF (X(1) .GT. XL(1)) THEN
  966.           THB=0.0D0
  967.        ENDIF
  968.        RETURN
  969.     ENDIF
  970.  
  971.     IF (INDEX .GT. 40) THEN
  972.        THB=0.0D0
  973.        RETURN
  974.         ENDIF
  975.  
  976.     IF (((T0+DTT) .GT. 1820.D0) .OR. ((T0+DTT) .LT. 600.0D0)) THEN
  977.        THB=0.0D0
  978.        RETURN
  979.     ENDIF
  980.  
  981.     IF (DABS(DTT) .GT. DTTMAX) THEN
  982.        DTT=SIGN(DTTMAX,DTT)
  983.     ENDIF
  984.  
  985.     IF (DABS(DXC) .GT. DXCMAX) THEN
  986.        DXC=SIGN(DXCMAX,DXC)
  987.     ENDIF
  988.  
  989.     T0=T0+DTT
  990.     XC=XC+DXC
  991.     GOTO 10
  992.    
  993.     END
  994.  
  995.  
  996.  
  997.     SUBROUTINE TABCALC(XL,TAB)
  998.     IMPLICIT REAL*8 (A-H,O-Z)
  999.     REAL*8 XL(0:10),XF(0:10),TAB
  1000.     REAL*8 DGAF0,DGLA0,DGLF0,DGFA(1:10),DGLF(1:10),SFE
  1001.     REAL*8 EIIL(1:10),EIIF(1:10),ECIL(1:10),ECIF(1:10)
  1002.     REAL*8 F(1:2),T0,T,DT,DTT,DTTMAX
  1003.     INTEGER I,J,INDEX
  1004.  
  1005.     T0=1811.15D0-2000*XL(1)
  1006.     DTTMAX=10.0D0
  1007.     INDEX=0
  1008.     DT=0.01D0
  1009.     DTTMAX=10.0D0
  1010.  
  1011. 10  INDEX=INDEX+1
  1012.     DO 100 J=1,2,1
  1013.        T=T0+(J-1.5D0)*DT
  1014.        RT=1.987D0*T
  1015.        CALL DGAF0CALC(T,DGAF0)
  1016.        CALL DGLFCALC(T,DGLF)
  1017.        CALL DGFACALC(T,DGFA)
  1018.        CALL SFECALC(T,SFE)
  1019.        CALL EIIFCALC(T,SFE,EIIF)
  1020.        CALL EIILCALC(T,EIIL)
  1021.        CALL ECIFCALC(T,ECIF)
  1022.        CALL ECILCALC(T,ECIL)
  1023.  
  1024.        DGLA0=-2646.0D0+39.1693D0*T-5.27D0*T*DLOG(T)+0.001D0*T*T
  1025.        DGLF0=DGLA0+DGAF0
  1026.  
  1027.        XF(1)=XL(1)*DEXP(DGLF(1)/RT)
  1028.        XF(1)=XF(1)*DEXP(EIIL(1)*XL(1)-EIIF(1)*XF(1))
  1029.        DO 20 I=2,10,1
  1030.           XF(I)=XL(I)*DEXP(DGLF(I)/RT+ECIL(I)*XL(1)-ECIF(I)*XF(1))
  1031.           XF(I)=XF(I)*DEXP(EIIL(I)*XL(I)-EIIF(I)*XF(I))
  1032.           XF(1)=XF(1)*DEXP(ECIL(I)*XL(I)-ECIF(I)*XF(I))
  1033. 20     CONTINUE
  1034.    
  1035.        GOTO 60
  1036.        DO 40 ITER=1,10,1
  1037.           XF(1)=XL(1)*DEXP(DGLF(1)/RT+EIIL(1)*XL(1)-EIIF(1)*XF(1))
  1038.           DO 40 I=2,10,1
  1039.              XF(I)=XL(I)*DEXP(DGLF(I)/RT+ECIL(I)*XL(1)-ECIF(I)*XF(1)&                         +EIIL(I)*XL(I)-EIIF(I)*XF(I))
  1040.              XF(1)=XF(1)*DEXP(ECIL(I)*XL(I)-ECIF(I)*XF(I))
  1041. 40        CONTINUE
  1042. 60     CONTINUE
  1043.              
  1044.        XF(0)=1.0-XF(1)
  1045.        F(J)=DGLF0/RT+(EIIF(1)*XF(1)*XF(1)-EIIL(1)*XL(1)*XL(1))/2.0D0
  1046.        DO 80 I=2,10,1
  1047.           XF(0)=XF(0)-XF(I)
  1048.           F(J)=F(J)+(EIIF(I)*XF(I)*XF(I)-EIIL(I)*XL(I)*XL(I))/2.0D0&           +ECIF(I)*XF(1)*XF(I)-ECIL(I)*XL(1)*XL(I)
  1049. 80     CONTINUE
  1050.  
  1051.        F(J)=F(J)+DLOG(XL(0)/XF(0))
  1052.  
  1053. 100 CONTINUE
  1054.  
  1055.     FAVG=(F(1)+F(2))/2.0D0
  1056.     DTT=-DT*FAVG/(F(2)-F(1))
  1057.    
  1058.     IF (DABS(DTT) .LE. 1.0D0) THEN
  1059.        TAB=T0+DTT-273.15D0
  1060.        RETURN
  1061.     ENDIF
  1062.  
  1063.     IF (((T0+DTT) .GT. 1820.0D0) .OR. ((T0+DTT) .LT. 600.0D0)) THEN
  1064.        TAB=0.0D0
  1065.        RETURN
  1066.     ENDIF
  1067.  
  1068.     IF (INDEX .GT. 40) THEN
  1069.        TAB=0.0D0
  1070.          RETURN
  1071.     ENDIF
  1072.  
  1073.     IF (DABS(DTT) .GT. DTTMAX) THEN
  1074.        DTT=SIGN(DTTMAX,DTT)
  1075.     ENDIF
  1076.  
  1077.     T0=T0+DTT
  1078.     GOTO 10
  1079.    
  1080.     END
  1081.  
  1082.  
  1083.     SUBROUTINE TJECALC (XA,TJE)
  1084.     IMPLICIT REAL*8 (A-H,O-Z)
  1085.     REAL*8 XL(0:10),XA(0:10),TJE
  1086.     REAL*8 DGLA0,DGFA(1:10),DGLA(1:10),F(1:2),FAVG
  1087.     REAL*8 EIIL(1:10),EIIA(1:10),ECIL(1:10),ECIA(1:10)
  1088.     REAL*8 T0,DT,T,RT,DTT,DTTMAX
  1089.     INTEGER INDEX,I
  1090.  
  1091.     T0=1808.0D0-4480.0D0*XA(1)
  1092.     DT=0.01D0
  1093.     DTTMAX=10.
  1094.     INDEX=0
  1095.            
  1096. 10  INDEX=INDEX+1
  1097.     DO 100 K=1,2,1
  1098.     T=T0+(K-1.5D0)*DT
  1099.     RT=1.987D0*T
  1100.       CALL DGLACALC (T,DGLA)
  1101.     CALL DGFACALC (T,DGFA)
  1102.     CALL EIILCALC (T,EIIL)
  1103.     CALL EIIACALC (T,EIIA)
  1104.     CALL ECILCALC (T,ECIL)
  1105.     CALL ECIACALC (T,ECIA)
  1106.  
  1107.     DGLA0=-2646.0D0+39.1693D0*T-5.27D0*T*DLOG(T)+0.001D0*T*T
  1108.  
  1109.     XL(1)=XA(1)*DEXP(-DGLA(1)/RT)
  1110.       XL(1)=XL(1)*DEXP(ECIA(1)*XA(1)-ECIL(1)*XL(1))
  1111.     DO 20 I=2,10,1
  1112.        XL(I)=XA(I)*DEXP(-DGLA(I)/RT+ECIA(I)*XA(1)-ECIL(I)*XL(1))
  1113.        XL(I)=XL(I)*DEXP(EIIA(I)*XA(I)-EIIL(I)*XL(I))
  1114.        XL(1)=XL(1)*DEXP(ECIA(I)*XA(I)-ECIL(I)*XL(I))
  1115. 20  CONTINUE
  1116.    
  1117.     XL(0)=1.0
  1118.     DO 60 I=1,10,1
  1119.        XL(0)=XL(0)-XL(I)
  1120. 60  CONTINUE
  1121.  
  1122.     F(K)=DGLA0/RT+DLOG(XL(0)/XA(0))&     +(EIIA(1)*XA(1)*XA(1)-EIIL(1)*XL(1)*XL(1))/2
  1123.     DO 100 I=2,10,1
  1124.        F(K)=F(K)+ECIA(I)*XA(1)*XA(I)-ECIL(I)*XL(1)*XL(I)&        +(EIIA(I)*XA(I)*XA(I)-EIIL(I)*XL(I)*XL(I))/2
  1125. 100 CONTINUE
  1126.  
  1127.     FAVG=(F(1)+F(2))/2
  1128.     DTT=-DT*FAVG/(F(2)-F(1))
  1129.  
  1130.     IF (DABS(DTT) .LE. 1.0) THEN
  1131.        TJE=T0+DTT-273.15D0
  1132.        RETURN
  1133.     ENDIF
  1134.  
  1135.     IF (INDEX .GE. 40) THEN
  1136.        TJE=0.0D0
  1137.        RETURN
  1138.     ELSEIF (DABS(DTT) .GT. DTTMAX) THEN
  1139.        DTT=SIGN(DTTMAX,DTT)
  1140.        T0=T0+DTT
  1141.        GOTO 10
  1142.     ELSE
  1143.        T0=T0+DTT
  1144.        GOTO 10
  1145.     ENDIF
  1146.  
  1147.     END
  1148.  
  1149.  
  1150.     SUBROUTINE TJBCALC(X,TJB)
  1151.     IMPLICIT REAL*8 (A-H,O-Z)
  1152.     REAL*8 X(0:10),XA(0:10),XF(0:10),XL(0:10),TJB
  1153.     REAL*8 DGAF0,DGFA0,DGFA(1:10),DGLA0,DGLA(1:10)
  1154.     REAL*8 EIIL(1:10),EIIF(1:10),EIIA(1:10)
  1155.     REAL*8 ECIL(1:10),ECIF(1:10),ECIA(1:10)
  1156.     REAL*8 A(1:10),B,F(2,2),G(2,2)
  1157.     REAL*8 T0,T,DT,DTT,XC,DX,DXC,Z,DTTMAX
  1158.     INTEGER INDEX
  1159.  
  1160.     T0=1765.0D0
  1161.     XC=0.008D0
  1162.     DTTMAX=10.0D0
  1163.     INDEX=0
  1164.     DX=1.0D-5
  1165.     DT=0.01D0
  1166.     DTTMAX=10.0D0
  1167.     DXCMAX=XC/2.0D0
  1168.  
  1169. 10  INDEX=INDEX+1
  1170.     DO 100 J=1,2,1
  1171.        T=T0+(J-1.5)*DT
  1172.        RT=1.987D0*T
  1173.     DO 100 K=1,2,1
  1174.        XA(1)=XC+(K-1.5D0)*DX
  1175.        CALL DGAF0CALC(T,DGAF0)
  1176.        CALL DGFACALC(T,DGFA)
  1177.        CALL DGLACALC(T,DGLA)
  1178.        CALL EIIACALC(T,EIIA)
  1179.        CALL EIIFCALC(T,SFE,EIIF)
  1180.        CALL EIILCALC(T,EIIL)
  1181.          CALL ECIACALC(T,ECIA)
  1182.        CALL ECIFCALC(T,ECIF)
  1183.        CALL ECILCALC(T,ECIL)
  1184.        DGFA0=-DGAF0
  1185.        DGLA0=-2646.0D0+39.1693D0*T-5.27D0*T*DLOG(T)+0.001D0*T*T
  1186.  
  1187.        XF(1)=XA(1)*DEXP(-DGFA(1)/RT)
  1188.        XF(1)=XF(1)*DEXP(EIIA(1)*XA(1)-EIIF(1)*XF(1))
  1189.        XL(1)=XA(1)*DEXP(-DGLA(1)/RT)
  1190.        XL(1)=XL(1)*DEXP(EIIA(1)*XA(1)-EIIL(1)*XL(1))
  1191.        B=(X(1)-XA(1))/(XL(1)-XA(1))
  1192.        DO 20 I=2,10,1
  1193.           A(I)=DEXP(-DGLA(I)/RT+ECIA(I)*XA(1)-ECIL(I)*XL(1))
  1194.           XA(I)=X(I)*(1+B*(A(I)-1))
  1195.           XL(I)=XA(I)*A(I)
  1196.           XL(1)=XL(1)*DEXP(ECIA(I)*XA(I)-ECIL(I)*XL(I))
  1197.           XF(I)=XA(I)*DEXP(-DGLA(I)/RT+ECIA(I)*XA(1)-ECIL(I)*XL(1))
  1198.           XF(1)=XF(1)*DEXP(ECIA(I)*XA(I)-ECIF(I)*XF(I))
  1199. 20     CONTINUE
  1200.  
  1201.        DO 50 ITER=1,10,1
  1202.           B=(X(1)-XA(1))/(XL(1)-XA(1))
  1203.           XF(1)=XA(1)*DEXP(-DGFA(1)/RT+EIIA(1)*XA(1)-EIIF(1)*XF(1))
  1204.           XL(1)=XA(1)*DEXP(-DGLA(1)/RT+EIIA(1)*XA(1)-EIIL(1)*XL(1))
  1205.            DO 40 I=2,10,1
  1206.             A(I)=DEXP(-DGLA(I)/RT+ECIA(I)*XA(1)-ECIL(I)*XL(1) &                 +EIIA(I)*XA(I)-EIIL(I)*XL(I))
  1207.             XA(I)=X(I)*(1+B*(A(I)-1))
  1208.             XL(I)=XA(I)*A(I)
  1209.             XL(1)=XL(1)*DEXP(ECIA(I)*XA(I)-ECIL(I)*XL(I))
  1210.             XF(I)=XA(I)*DEXP(-DGLA(I)/RT+ECIA(I)*XA(1)-ECIL(I)*XL(1))
  1211.             XF(1)=XF(1)*DEXP(ECIA(I)*XA(I)-ECIF(I)*XF(I))
  1212. 40       CONTINUE
  1213. 50      CONTINUE
  1214.              
  1215.          XA(0)=1.0D0
  1216.        XF(0)=1.0D0
  1217.        XL(0)=1.0D0
  1218.        DO 60 I=1,10,1
  1219.           XA(0)=XA(0)-XA(I)
  1220.           XF(0)=XF(0)-XF(I)
  1221.           XL(0)=XL(0)-XL(I)
  1222. 60     CONTINUE
  1223.  
  1224.        F(J,K)=DGFA0/RT+DLOG(XF(0)/XA(0))&          +(EIIA(1)*XA(1)*XA(1)-EIIF(1)*XF(1)*XF(1))/2
  1225.        G(J,K)=DGLA0/RT+DLOG(XL(0)/XA(0)) &          +(EIIA(1)*XA(1)*XA(1)-EIIL(1)*XL(1)*XL(1))/2
  1226.        DO 80 I=2,10,1
  1227.           F(J,K)=F(J,K)+(EIIA(I)*XA(I)*XA(I)-EIIF(I)*XF(I)*XF(I))/2&             +ECIA(I)*XA(1)*XA(I)-ECIF(I)*XF(1)*XF(I)
  1228.           G(J,K)=G(J,K)+(EIIA(I)*XA(I)*XA(I)-EIIL(I)*XL(I)*XL(I))/2&             +ECIA(I)*XA(1)*XA(I)-ECIL(I)*XL(1)*XL(I)
  1229. 80     CONTINUE
  1230. 100 CONTINUE
  1231.  
  1232.     DFX=(F(1,2)-F(1,1)+F(2,2)-F(2,1))/DX/2
  1233.     DGX=(G(1,2)-G(1,1)+G(2,2)-G(2,1))/DX/2
  1234.     DFT=(F(2,1)-F(1,1)+F(2,2)-F(1,2))/DT/2
  1235.     DGT=(G(2,1)-G(1,1)+G(2,2)-G(1,2))/DT/2
  1236.     Z=DFX*DGT-DFT*DGX
  1237.     FAVG=(F(1,1)+F(1,2)+F(2,1)+F(2,2))/4
  1238.     GAVG=(G(1,1)+G(1,2)+G(2,1)+G(2,2))/4
  1239.     DXC=(-FAVG*DGT+GAVG*DFT)/Z
  1240.     DTT=(-GAVG*DFX+FAVG*DGX)/Z
  1241.  
  1242.     IF ((DABS(DXC/XC) .LE. 0.01D0) .AND. (DABS(DTT) .LE. 1.0D0)) THEN
  1243.        IF ((X(1) .GE. XF(1)) .AND. (X(1) .LE. XL(1))) THEN
  1244.           TJB=T0+DTT-273.15D0
  1245.           RETURN
  1246.        ELSE
  1247.           TJB=0.0D0
  1248.           RETURN
  1249.        ENDIF
  1250.     ENDIF
  1251.  
  1252.     IF ((T0+DTT) .GT. 1820.0D0) THEN
  1253.        TJB=0.0D0
  1254.        RETURN
  1255.     ENDIF
  1256.  
  1257.     IF (INDEX .GT. 20) THEN
  1258.        TJB=0.0
  1259.        RETURN
  1260.     ENDIF
  1261.  
  1262.     IF (DABS(DTT) .GT. DTTMAX) THEN
  1263.        DTT=SIGN(DTTMAX,DTT)
  1264.     ENDIF
  1265.  
  1266.     IF (DABS(DXC) .GT. DXCMAX) THEN
  1267.        DXC=SIGN(DXCMAX,DXC)
  1268.     ENDIF
  1269.  
  1270.     T0=T0+DTT
  1271.     XC=XC+DXC
  1272.     GOTO 10
  1273.    
  1274.     END
  1275.  
  1276.  
  1277.  
  1278.     SUBROUTINE TBCCALC (XL,TBC)
  1279.     IMPLICIT REAL*8 (A-H,O-Z)
  1280.     REAL*8 XL(0:10),XA(0:10),TBC
  1281.     REAL*8 DGLA0,DGLF(1:10),F(1:2),FAVG,DGLA(1:10),DGFA(1:10)
  1282.     REAL*8 EIIL(1:10),EIIA(1:10),ECIL(1:10),ECIA(1:10)
  1283.     REAL*8 T0,DT,T,RT,DTT,DTTMAX
  1284.     INTEGER INDEX,ITER,I
  1285.  
  1286.     T0=1818.0D0-1820.0D0*XL(1)
  1287.     DT=0.01D0
  1288.     DTTMAX=10.0D0
  1289.     INDEX=0
  1290.  
  1291. 10  INDEX=INDEX+1
  1292.     DO 100 K=1,2,1
  1293.     T=T0+(K-1.5D0)*DT
  1294.     RT=1.987D0*T
  1295.     CALL DGLFCALC (T,DGLF)
  1296.     CALL DGFACALC (T,DGFA)
  1297.     CALL DGLACALC (T,DGLA)
  1298.     CALL EIILCALC (T,EIIL)
  1299.     CALL EIIACALC (T,EIIA)
  1300.     CALL ECILCALC (T,ECIL)
  1301.     CALL ECIACALC (T,ECIA)
  1302.  
  1303.     DGLA0=-2646.0D0+39.1693D0*T-5.27D0*T*DLOG(T)+0.001D0*T*T
  1304.  
  1305.     DO 12 I=1,10,1
  1306.        DGLA(I)=DGLF(I)+DGFA(I)
  1307. 12    CONTINUE
  1308.     XA(1)=XL(1)*DEXP(DGLA(1)/RT)
  1309.       XA(1)=XA(1)*DEXP(ECIL(1)*XL(1)-ECIA(1)*XA(1))
  1310.     DO 20 I=2,10,1
  1311.        XA(I)=XL(I)*DEXP(DGLA(I)/RT+ECIL(I)*XL(1)-ECIA(I)*XA(1))
  1312.        XA(I)=XA(I)*DEXP(EIIL(I)*XL(I)-EIIA(I)*XA(I))
  1313.        XL(1)=XL(1)*DEXP(ECIL(I)*XL(I)-ECIA(I)*XA(I))
  1314. 20  CONTINUE
  1315.  
  1316.     DO 40 ITER=1,10,1
  1317.        XA(1)=XL(1)*DEXP(DGLA(1)/RT+ECIL(1)*XL(1)-ECIA(1)*XA(1))
  1318.     DO 40 I=2,10,1
  1319.        XA(I)=XL(I)*DEXP(DGLA(I)/RT+ECIL(I)*XA(1)-ECIA(I)*XA(1&                 +EIIL(I)*XL(I)-EIIA(I)*XA(I))
  1320.        XA(1)=XA(1)*DEXP(ECIL(I)*XL(I)-ECIA(I)*XA(I))
  1321. 40  CONTINUE
  1322.  
  1323. 50  XA(0)=1.0
  1324.     DO 60 I=1,10,1
  1325.        XA(0)=XA(0)-XA(I)
  1326. 60  CONTINUE
  1327.  
  1328.     F(K)=DGLA0/RT+DLOG(XL(0)/XA(0)) &     +(EIIA(1)*XA(1)*XA(1)-EIIL(1)*XL(1)*XL(1))/2
  1329.     DO 100 I=2,10,1
  1330.        F(K)=F(K)+ECIA(I)*XA(1)*XA(I)-ECIL(I)*XL(1)*XL(I)&        +EIIA(I)*XA(I)*XA(I)/2-EIIL(I)*XL(I)*XL(I)/2
  1331. 100 CONTINUE
  1332.  
  1333.     FAVG=(F(1)+F(2))/2
  1334.     DTT=-DT*FAVG/(F(2)-F(1))
  1335.  
  1336.     IF (DABS(DTT) .LE. 1.0D0) THEN
  1337.        TBC=T0+DTT-273.15D0
  1338.        RETURN
  1339.     ENDIF
  1340.  
  1341.     IF (INDEX .GE. 40) THEN
  1342.        TBC=0.0D0
  1343.        RETURN
  1344.     ELSEIF (DABS(DTT) .GT. DTTMAX) THEN
  1345.        DTT=SIGN(DTTMAX,DTT)
  1346.        T0=T0+DTT
  1347.        GOTO 10
  1348.     ELSE
  1349.        T0=T0+DTT
  1350.        GOTO 10
  1351.     ENDIF
  1352.  
  1353.     END
  1354.  
  1355.     SUBROUTINE SFECALC(TEMP,SFE)
  1356.     IMPLICIT REAL*8 (A-H,O-Z)
  1357.     REAL*8 TEMP,SFE
  1358.     REAL*8 T(15),S(15)
  1359.  
  1360.     DATA T /400.D0,500.D0,600.D0,700.D0,800.D0,900.D0,1000.D0,1100.D0,&         1200.D0,1300.D0,1400.D0,1500.D0,1600.D0,1700.D0,1800.D0/      
  1361.     DATA S/0.05D0,0.09D0,0.16D0,0.28D0,0.45D0,0.72D0,1.16D0,1.67D0, &          1.84D0,1.94D0,1.99D0,2.02D0,2.04D0,2.05D0,2.06D0/
  1362.  
  1363.     DO 50 I=1,14,1
  1364.        IF ((TEMP .GE. T(I)) .AND. (TEMP .LE. T(I+1))) THEN
  1365.           SFE=S(I)+(S(I+1)-S(I))*(TEMP-T(I))/(T(I+1)-T(I))
  1366.        GOTO 100
  1367.        ENDIF
  1368. 50  CONTINUE
  1369.  
  1370. 100 RETURN
  1371.     END
  1372.  
  1373.  
  1374.  
  1375.     SUBROUTINE DGAF0CALC (TEMP,DGAF0)
  1376.     IMPLICIT REAL*8 (A-H,O-Z)
  1377.     REAL*8 T(126),G(126)
  1378.     REAL*8 TEMP,DGAF0
  1379.  
  1380.     DATA T /600.0D0,610.0D0,620.0D0,630.0D0,640.0D0,650.0D0,660.0D0,&          670.0D0,680.0D0,690.0D0,700.0D0,710.0D0,720.0D0,730.0D0,
  1381.      &          740.0D0,750.0D0,760.0D0,770.0D0,780.0D0,790.0D0,800.0D0,&          810.0D0,820.0D0,830.0D0,840.0D0,850.0D0,860.0D0,870.0D0,
  1382.      &          880.0D0,890.0D0,900.0D0,910.0D0,920.0D0,930.0D0,940.0D0, &          950.0D0,960.0D0,970.0D0,980.0D0,990.0D0,1000.D0,1010.D0,
  1383.      &          1020.D0,1030.D0,1040.D0,1050.D0,1060.D0,1070.D0,1080.D0, &        1090.D0,1100.D0,1110.D0,1120.D0,1130.D0,1140.D0,1150.D0,&          1160.D0,1170.D0,1180.D0,1184.D0,1190.D0,1200.D0,1210.D0,&          1220.D0,1230.D0,1240.D0,1250.D0,1260.D0,1270.D0,1280.D0,&          1290.D0,1300.D0,1310.D0,1320.D0,1330.D0,1340.D0,1350.D0,&          1360.D0,1370.D0,1380.D0,1390.D0,1400.D0,1410.D0,1420.D0,&          1430.D0,1440.D0,1450.D0,1460.D0,1470.D0,1480.D0,1490.D0,&          1500.D0,1510.D0,1520.D0,1530.D0,1540.D0,1550.D0,1560.D0, &          1570.D0,1580.D0,1590.D0,1600.D0,1610.D0,1620.D0,1630.D0,&          1640.D0,1650.D0,1660.D0,1668.D0,1670.D0,1680.D0,1690.D0,&          1700.D0,1710.D0,1720.D0,1730.D0,1740.D0,1750.D0,1760.D0,&          1770.D0,1780.D0,1790.D0,1800.D0,1810.D0,1820.D0,1830.D0/
  1384.  
  1385.     DATA G /677.4D0,658.6D0,640.0D0,621.6D0,603.3D0,585.2D0,567.3D0, &          549.5D0,531.9D0,514.4D0,497.0D0,479.7D0,462.6D0,445.7D0,&          429.0D0,412.5D0,396.1D0,379.9D0,364.0D0,348.4D0,333.0D0,&          317.8D0,302.8D0,288.1D0,273.5D0,259.1D0,244.7D0,230.8D0,&          217.3D0,204.1D0,191.1D0,178.3D0,166.0D0,154.1D0,142.4D0,&          131.0D0,119.9D0,109.3D0,99.30D0,89.80D0,80.80D0,72.18D0,&          64.16D0,56.75D0,50.15D0,44.31D0,39.11D0,34.21D0,29.92D0,&          26.07D0,22.22D0,18.37D0,14.92D0,11.84D0, 9.10D0, 6.66D0,&           4.48D0, 2.52D0, 0.71D0, 0.00D0,-0.89D0,-2.40D0,-3.82D0,&          -5.15D0,-6.39D0,-7.55D0,-8.63D0,-9.63,-10.55,-11.40D0,&   -12.17D0,-12.87D0,-13.50D0,-14.06D0,-14.54D0,-14.95D0,-15.29D0,&   -15.57D0,-15.79D0,-15.95D0,-16.05D0,-16.10D0,-16.10D0,-16.06D0,&   -15.97D0,-15.83D0,-15.64D0,-15.40D0,-15.11D0,-14.77D0,-14.38D0,&   -13.95D0,-13.48D0,-12.97D0,-12.41D0,-11.80D0,-11.15D0,-10.45D0,&          -9.92D0,-9.71D0,-8.98D0,-7.20D0,-6.27D0,-5.30D0,-4.28D0,&          -3.22D0,-2.11D0,-0.96D0, 0.00D0, 0.23D0, 1.47D0, 2.75D0,&           4.08D0, 5.45D0, 6.86D0, 8.31D0, 9.81D0,11.36D0,12.97D0,&          14.63D0,16.35D0,18.13D0,19.98D0,21.83D0,23.68D0,25.53D0/
  1386.  
  1387.     DO 10 I=1,126,1
  1388.        IF ((TEMP .GT. T(I)) .AND. (TEMP .LE. T(I+1))) THEN
  1389.           DGAF0=G(I)+(TEMP-T(I))*(G(I+1)-G(I))/(T(I+1)-T(I))
  1390.           GOTO 100
  1391.        ENDIF
  1392. 10  CONTINUE
  1393.     WRITE(*,*) '****** ERROR ******'
  1394.     WRITE(*,*) 'TEMPERATURE OUT OF THE RANGE 600-1800 K'
  1395.     STOP
  1396.  
  1397. 100 RETURN
  1398.     END
  1399.  
  1400.  
  1401.     SUBROUTINE DGAFCALC (TEMP,DGAF)
  1402.     IMPLICIT REAL*8 (A-H,O-Z)
  1403.  
  1404.  
  1405.     REAL*8 TEMP,DGAF(1:10)
  1406.  
  1407.     DGAF(1)=-15323.0D0+7.686D0*TEMP
  1408.     DGAF(2)=-26650.0D0+42.69D0*TEMP-0.017D0*TEMP*TEMP
  1409.     DGAF(3)=-5954.0D0+38.7996D0*TEMP-4.7244D0*TEMP*DLOG(TEMP)
  1410.     DGAF(4)=-4545.0D0+3.233D0*TEMP
  1411.     DGAF(5)=-367.0D0-4.656D0*TEMP+0.6568D0*TEMP*DLOG(TEMP)
  1412.     DGAF(6)=565+0.15*TEMP
  1413.     DGAF(7)=-25500.0D0+41.183D0*TEMP-0.017D0*TEMP*TEMP
  1414.     DGAF(8)=1290.0D0-0.3D0*TEMP
  1415.     DGAF(9)=-8357.0D0+13.8D0*TEMP-0.0051D0*TEMP*TEMP
  1416.     DGAF(10)=5160.0D0-2.875D0*TEMP
  1417.  
  1418.     RETURN
  1419.     END
  1420.  
  1421.     SUBROUTINE DGFACALC (TEMP,DGFA)
  1422.     IMPLICIT REAL*8 (A-H,O-Z)
  1423.  
  1424.  
  1425.     REAL*8 TEMP,DGFA(1:10)
  1426.  
  1427.     DGFA(1)=15940.0D0-5.7D0*TEMP
  1428.     DGFA(2)=430.0D0-0.305D0*TEMP
  1429.     DGFA(3)=5954.0D0-38.799D0*TEMP+4.7244D0*TEMP*DLOG(TEMP)
  1430.     DGFA(4)=1330.0D0-0.26D0*TEMP
  1431.     DGFA(5)=367.0D0+4.646D0*TEMP-0.6568D0*TEMP*DLOG(TEMP)
  1432.     DGFA(6)=-565.0D0-0.15D0*TEMP
  1433.     DGFA(7)=-1450.0D0-0.8D0*TEMP
  1434.     DGFA(8)=-2500.0D0-0.15D0*TEMP
  1435.     DGFA(9)=-8357.0D0-13.8D0*TEMP+0.0051D0*TEMP*TEMP
  1436.     DGFA(10)=-5700.0D0+0.04D0*TEMP
  1437.  
  1438.     RETURN
  1439.     END
  1440.  
  1441.  
  1442.  
  1443.     SUBROUTINE DGLFCALC (TEMP,DGLF)
  1444.     IMPLICIT REAL*8 (A-H,O-Z)
  1445.  
  1446.     REAL*8 TEMP,DGLF(1:10)
  1447.  
  1448.     DGLF(1)=-21300.0D0+6.3D0*TEMP
  1449.     DGLF(2)=3500.0D0-2.308D0*TEMP
  1450.     DGLF(3)=-8200.0D0+3.9D0*TEMP
  1451.     DGLF(4)=-2120.0D0-0.38D0*TEMP
  1452.     DGLF(5)=4600.0D0+2.19D0*TEMP
  1453.     DGLF(6)=-6600.0D0+2.29D0*TEMP
  1454.     DGLF(7)=-1200.0D0
  1455.     DGLF(8)=-7500.0D0+3.65D0*TEMP
  1456.     DGLF(9)=-5100.0D0+2.3D0*TEMP
  1457.     DGLF(10)=-5500.0D0+2.3D0*TEMP
  1458.  
  1459.     RETURN
  1460.     END
  1461.  
  1462.  
  1463.     SUBROUTINE DGLACALC (TEMP,DGLA)
  1464.     IMPLICIT REAL*8 (A-H,O-Z)
  1465.  
  1466.     REAL*8 TEMP,DGLA(1:10)
  1467.  
  1468.     DGLA(1)=-5360.0D0+0.6D0*TEMP
  1469.     DGLA(2)=3930.0D0-2.613D0*TEMP
  1470.     DGLA(3)=-2246.0D0-34.899D0*TEMP+4.7244D0*TEMP*DLOG(TEMP)
  1471.     DGLA(4)=-790.0D0-0.64D0*TEMP
  1472.     DGLA(5)=-4233.0D0+6.89D0*TEMP-0.6568D0*TEMP*DLOG(TEMP)
  1473.     DGLA(6)=-7165.0D0+2.14D0*TEMP
  1474.     DGLA(7)=-2650.0D0+0.8D0*TEMP
  1475.     DGLA(8)=-8790.0D0+3.95D0*TEMP
  1476.     DGLA(9)=3257.0D0-11.5D0*TEMP+0.0051D0*TEMP*TEMP
  1477.     DGLA(10)=-10661.0D0+5.175D0*TEMP
  1478.  
  1479.     RETURN
  1480.     END
  1481.  
  1482.     SUBROUTINE EIIFCALC (TEMP,SFE,EIIF)
  1483.     IMPLICIT REAL*8 (A-H,O-Z)
  1484.     REAL*8 TEMP,SFE,EIIF(1:10)
  1485.  
  1486.     EIIF(1)=1.3D0
  1487.     EIIF(2)=-3.082D0+4679.0D0/TEMP+1509.8D0*SFE/TEMP
  1488.     EIIF(3)=-16.35D0+44829.0D0/TEMP
  1489.     EIIF(4)=2.013D0-4595.0D0/TEMP+385.5D0*SFE/TEMP
  1490.     EIIF(5)=2.819D0-6039.0D0/TEMP
  1491.     EIIF(6)=-0.219D0-4772.0D0/TEMP+402.6D0*SFE/TEMP
  1492.     EIIF(7)=0.634D0-11270.0D0/TEMP+1006.5D0*SFE/TEMP
  1493.     EIIF(8)=-9058.0D0/TEMP
  1494.     EIIF(9)=9.863D0-3849.0D0/TEMP
  1495.     EIIF(10)=9.863D0-3849.0D0/TEMP
  1496.  
  1497.     RETURN
  1498.     END
  1499.  
  1500.     SUBROUTINE EIIACALC (TEMP,EIIA)
  1501.     IMPLICIT REAL*8 (A-H,O-Z)
  1502.     REAL*8 TEMP,EIIA(1:10)
  1503.  
  1504.     EIIA(1)=4.786D0+5066.D0/TEMP
  1505.     EIIA(2)=2.406D0-175.0D0/TEMP
  1506.     EIIA(3)=26048.0D0/TEMP
  1507.     EIIA(4)=-2839.0D0/TEMP
  1508.     EIIA(5)=7.655D0-3154.0D0/TEMP-0.661D0*DLOG(TEMP)
  1509.     EIIA(6)=-2330.0D0/TEMP
  1510.     EIIA(7)=-1.376D0-5766.0D0/TEMP
  1511.     EIIA(8)=0.453D0-7851.0D0/TEMP
  1512.     EIIA(9)=11.366D0-3524.0D0/TEMP
  1513.     EIIA(10)=11.366D0-3524.0D0/TEMP
  1514.  
  1515.     RETURN
  1516.     END
  1517.  
  1518.     SUBROUTINE EIILCALC (TEMP,EIIL)
  1519.     IMPLICIT REAL*8 (A-H,O-Z)
  1520.     REAL*8 TEMP,EIIL(1:10)
  1521.  
  1522.     EIIL(1)=0.389D0+7810.0D0/TEMP
  1523.     EIIL(2)=31.4D0/TEMP
  1524.     EIIL(3)=24751.0D0/TEMP
  1525.     EIIL(4)=312.6D0/TEMP
  1526.     EIIL(5)=9.8D0/TEMP
  1527.     EIIL(6)=-1343.0D0/TEMP
  1528.     EIIL(7)=11011.0D0/TEMP
  1529.     EIIL(8)=-4289.0D0/TEMP
  1530.     EIIL(9)=-6056.0D0/TEMP
  1531.     EIIL(10)=-6056.0D0/TEMP
  1532.  
  1533.     RETURN
  1534.     END
  1535.  
  1536.     SUBROUTINE ECIFCALC (TEMP,ECIF)
  1537.     IMPLICIT REAL*8 (A-H,O-Z)
  1538.     REAL*8 TEMP,ECIF(1:10)
  1539.    
  1540.     ECIF(1)=1.3D0
  1541.     ECIF(2)=-5834.0D0/TEMP
  1542.     ECIF(3)=16205.0D0/TEMP
  1543.     ECIF(4)=5533.0D0/TEMP
  1544.     ECIF(5)=1.558-6160/TEMP
  1545.     ECIF(6)=-10714.0D0/TEMP
  1546.     ECIF(7)=6.615-5533/TEMP
  1547.     ECIF(8)=-3447.0D0/TEMP
  1548.     ECIF(9)=-21650.0D0/TEMP
  1549.     ECIF(10)=-28770.0D0/TEMP
  1550.  
  1551.     RETURN
  1552.     END
  1553.  
  1554.     SUBROUTINE ECIACALC (TEMP,ECIA)
  1555.     IMPLICIT REAL*8 (A-H,O-Z)
  1556.     REAL*8 TEMP,ECIA(1:10)
  1557.    
  1558.     ECIA(1)=4.786D0+5066.0D0/TEMP  
  1559.     ECIA(2)=-4811.0D0/TEMP
  1560.     ECIA(3)=14794.0D0/TEMP
  1561.     ECIA(4)=5533.0D0/TEMP
  1562.     ECIA(5)=14.193-30209.0D0/TEMP
  1563.     ECIA(6)=-10714.0D0/TEMP
  1564.     ECIA(7)=6.615D0-5533.0D0/TEMP
  1565.     ECIA(8)=-10342.0D0/TEMP
  1566.     ECIA(9)=-21650.0D0/TEMP
  1567.     ECIA(10)=-28770.0D0/TEMP
  1568.  
  1569.     RETURN
  1570.     END
  1571.  
  1572.     SUBROUTINE ECILCALC (TEMP,ECIL)
  1573.     IMPLICIT REAL*8 (A-H,O-Z)
  1574.     REAL*8 TEMP,ECIL(1:10)
  1575.  
  1576.     ECIL(1)=0.389D0+7810.0D0/TEMP
  1577.     ECIL(2)=-5060.0D0/TEMP
  1578.     ECIL(3)=-0.428D0+18740.0D0/TEMP
  1579.     ECIL(4)=5340.0D0/TEMP
  1580.     ECIL(5)=-9500.0D0/TEMP
  1581.     ECIL(6)=-7490.0D0/TEMP
  1582.     ECIL(7)=7586.0D0/TEMP
  1583.     ECIL(8)=-12230.0D0/TEMP
  1584.     ECIL(9)=-30100.0D0/TEMP
  1585.     ECIL(10)=-46615.0D0/TEMP
  1586.  
  1587.     RETURN
  1588.     END
  1589.    
  1590.     SUBROUTINE GCACALC (TEMP,GCA)
  1591.     IMPLICIT REAL*8 (A-H,O-Z)
  1592.     REAL*8 TEMP,GCA(0:10)
  1593.  
  1594.     GCA(0)=13320.D0-64.72D0*TEMP+7.4833D0*TEMP*DLOG(TEMP)
  1595.     GCA(1)=46150.D0-19.221D0*TEMP
  1596.     GCA(2)=-13532.+730.0D0+10.0D0*TEMP
  1597.     GCA(4)=14540.D0-2.367D0*TEMP-(-14600.D0+8800.D0)
  1598.     GCA(5)=-850.D0-14.58D0*TEMP-(10460.D0+0.628D0*TEMP)&         -(13110.D0-31.820*TEMP+2.748D0*TEMP*DLOG(TEMP))
  1599.     GCA(6)=502.D0-(10460.D0+0.628D0*TEMP)-9686.
  1600.     GCA(8)=63200.0D0-30.392*TEMP-(10460.D0+0.628D0*TEMP)&         -(32635.0D0-1.883D0*TEMP)
  1601.     GCA(9)=-12763.D0-14.345D0*TEMP-(9000.D0+3.556D0*TEMP)&         -(-21350.D0-12.75D0*TEMP)
  1602.     GCA(10)=GCA(9)
  1603.     DO 10 I=0,10,1
  1604.        GCA(I)=(GCA(I)-GCA(1)/3.0)/4.184D0
  1605. 10  CONTINUE
  1606.  
  1607.     RETURN
  1608.     END
  1609.  
  1610.     SUBROUTINE  A01CALC (TEMP,A0,A1)
  1611.     IMPLICIT REAL*8 (A-H,O-Z)
  1612.     REAL*8 TEMP,A0(2:10), A1(2:10)
  1613.     A0(2)=1996.D0-3.633D0*TEMP
  1614.     A0(4)=0.0D0
  1615.     A0(5)=427.8D0
  1616.     A0(6)=0.0D0
  1617.     A0(8)=0.0D0
  1618.     A0(9)=-7170.D0
  1619.     A0(10)=-7170.D0
  1620.     A1(2)=-956.D0
  1621.  
  1622.     RETURN
  1623.     END
  1624.  
  1625.     SUBROUTINE TDISSCALC (C,Mn,Cr,Mo,V,Nb,Ti,Al,N,&                          AE3,TDISS)
  1626.     IMPLICIT REAL*8 (A-H,O-Z)
  1627.     PARAMETER (ZERO=1.0D-6)
  1628.     REAL*8 C,Mn,Cr,Mo,V,Nb,Ti,Al,N,AE3,TDISS
  1629.     REAL*8 NbC,TiC,CrC,VC,MoC,VN,VN2,AlN,AlN2,NbN,NbN2,TiN
  1630.     IF (Nb .EQ. 0.D0) THEN
  1631.        Nb=ZERO
  1632.     ENDIF
  1633.  
  1634.     IF(Ti .EQ. 0.D0) THEN
  1635.        Ti=ZERO
  1636.     ENDIF
  1637.  
  1638.     IF (Cr .EQ. 0.D0) THEN
  1639.        Cr=ZERO
  1640.     ENDIF
  1641.  
  1642.     IF (Mo .EQ. 0.D0) THEN
  1643.        Mo=ZERO
  1644.     ENDIF
  1645.  
  1646.     IF (Al .EQ. 0.D0) THEN
  1647.        Al=ZERO
  1648.     ENDIF
  1649.  
  1650.     IF (V .EQ. 0.D0) THEN
  1651.        V=ZERO
  1652.     ENDIF
  1653.  
  1654.     IF (N .EQ. 0.D0) THEN
  1655.        N=ZERO
  1656.     ENDIF
  1657.  
  1658.     NbC=7520.D0/(3.11D0-DLOG10(Nb*C))-273.D0
  1659.     TiC=7000.D0/(2.75D0-DLOG10(Ti*C))-273.D0
  1660.     CrC=7375.D0/(5.90D0-DLOG10((Cr**23.0D0)*(C**6.0D0)))-273.D0
  1661.     VC=8000.D0/(5.36D0-DLOG10((V**4.0D0)*(C**3.0D0)))-273.D0   
  1662.     MoC=7375.D0/(5.0D0-DLOG10((Mo**2.0D0)*C))-273.D0
  1663.     VN=7070.D0/(2.27D0-DLOG10(V*N))-273.D0
  1664.     VN2=8330.D0/(3.46D0-0.012*Mn-DLOG10(V*N))-273.D0
  1665.     AlN=7750.D0/(1.8D0-DLOG10(Al*N))-273.15D0
  1666.     AlN2=6770.D0/(1.03D0-DLOG10(Al*N))-273.D0
  1667.     NbN=10230.D0/(4.04D0-DLOG10(Nb*N))-273.D0
  1668.     NbN2=10800.D0/(3.7D0-DLOG10(Nb*N))-273.D0
  1669.     TiN=8000.D0/(0.32D0-DLOG10(Ti*N))-273.D0   
  1670.     TDISS=MAX(NbC,TiC,CrC,VC,MoC,VN,VN2,AlN,AlN2,NbN,&          NbN2,TiN,AE3)
  1671.  
  1672.     RETURN
  1673.     END
  1674.  
  1675.     SUBROUTINE STARTSTOPH (TEMP1,TEMP2,TIME1,TIME2,&                       T3,T4,TEMP0,DTEMP,DTIME)
  1676.     IMPLICIT REAL*8 (A-H,O-Z)
  1677.  
  1678.     REAL*8 TEMP1,TEMP2,TIME1,TIME2,T3,T4,TEMP0,DTEMP,DTIME
  1679.     REAL*8 TIME0
  1680.  
  1681.     IF (TEMP1 .GE. T3) THEN
  1682.        TEMP0=TEMP1
  1683.        TIME0=TIME1
  1684.     ELSE
  1685.        TEMP0=T3
  1686.        TIME0=TIME1+(TIME2-TIME1)*(TEMP0-TEMP1)/(TEMP2-TEMP1)
  1687.     ENDIF
  1688.  
  1689.     IF (TEMP2 .LE. T4) THEN
  1690.        DTEMP=TEMP2-TEMP0
  1691.        DTIME=TIME2-TIME0
  1692.     ELSE
  1693.        DTEMP=T4-TEMP0
  1694.        DTIME=TIME1-TIME0+(TIME2-TIME1)*(T4-TEMP1)/(TEMP2-TEMP1)
  1695.     ENDIF    
  1696.  
  1697.     RETURN
  1698.     END
  1699.  
  1700.  
  1701.     SUBROUTINE STARTSTOPC (TEMP1,TEMP2,TIME1,TIME2, &   T3,T4,TEMP0,DTEMP,DTIME)
  1702.     IMPLICIT REAL*8 (A-H,O-Z)
  1703.     REAL*8 TEMP1,TEMP2,TIME1,TIME2,T3,T4,TEMP0,DTEMP,DTIME
  1704.     REAL*8 TIME0
  1705.     IF (TEMP1 .LT. T3) THEN
  1706.        TEMP0=TEMP1
  1707.        TIME0=TIME1
  1708.     ELSE
  1709.        TEMP0=T3
  1710.        TIME0=TIME1+(TIME2-TIME1)*(T3-TEMP1)/(TEMP2-TEMP1)
  1711.     ENDIF
  1712.  
  1713.     IF (TEMP2 .GT. T4) THEN
  1714.        DTEMP=TEMP2-TEMP0
  1715.        DTIME=TIME2-TIME0
  1716.     ELSE
  1717.        DTEMP=T4-TEMP0
  1718.        DTIME=TIME1-TIME0+(TIME2-TIME1)*(T4-TEMP1)/(TEMP2-TEMP1)
  1719.     ENDIF
  1720.  
  1721.     RETURN
  1722.     END
  1723.  
  1724.     SUBROUTINE GRAINGROWTH (TEMP,TIMEINC,GS,gsmax)
  1725.     IMPLICIT REAL*8 (A-H,O-Z)
  1726.     REAL*8 TIMEINC,TEMP,GS,gsmax
  1727.     REAL*8 DGS,QR,T
  1728.     REAL*8 N
  1729.  
  1730.       N=4.0D0
  1731.       QR=69300.0D0
  1732.     T=TEMP+273.15D0
  1733.     DGS=2.969D15*TIMEINC*DEXP(-QR/T)/N/GS**(N-1.0D0)
  1734.     GS=GS+DGS
  1735.       IF (GS .gt. gsmax) then
  1736.        GS=gsmax
  1737.     ENDIF
  1738.  
  1739.     RETURN
  1740.     END
  1741.  
  1742.  
  1743.     REAL*8 FUNCTION XTERM (X)
  1744.  
  1745.     REAL*8 X
  1746.  
  1747.     IF ( X .GT. 1.0D0 ) THEN
  1748.         X = 1.0D0
  1749.     END IF
  1750.  
  1751.     IF (X .LT. 1.0D-6) THEN
  1752.        X=1.0D-6
  1753.     ENDIF
  1754.     XTERM=X**(0.4D0*(1.0D0-X))*(1.0D0-X)**(0.4D0*X)
  1755.  
  1756.     RETURN
  1757.     END
  1758.  
  1759.  
  1760.     REAL*8 FUNCTION DIFF(TEMP)
  1761.     REAL*8 TEMP
  1762.     DIFF=DEXP(-27500.D0/1.987D0/(TEMP+273.15D0))
  1763.     RETURN
  1764.     END
  1765.  
  1766.  
  1767.     SUBROUTINE HARDNESS (TEMP2,VR,NNODE)
  1768.     IMPLICIT REAL*8 (A-H,O-Z)
  1769.     PARAMETER (NN=10000)
  1770.     real*8 C,Mn,Si,Ni,Cr,Mo,Cu,W,V,Nb,Ti,Al,N
  1771.     REAL*8 C2,Mn2,Si2,Ni2,Cr2,Mo2,Cu2,W2,V2,Nb2,Ti2,Al2,N2
  1772.     REAL*8 CE,HM,HB,HVM,HVB,HVFP
  1773.     REAL*8 VR(NN),TEMP2(NN)
  1774.     REAL*8 XA(NN),XF(NN),XP(NN),XB(NN),XM(NN),GS(NN),G(NN),HV(NN),XL(NN)
  1775.     real*8 AC1_1,AC3_1,TDISS1,AE4_1,AE3_1,AE1_1,TBS1,TMS1,as1_1&    ,Tliquid1,Tsolid1
  1776.     real*8 AC1_2,AC3_2,TDISS2,AE4_2,AE3_2,AE1_2,TBS2,TMS2,as1_2&    ,Tliquid2,Tsolid2
  1777.  
  1778.     COMMON/COMP/C,Mn,Si,Ni,Cr,Mo,Cu,W,V,Nb,Ti,Al,N
  1779.     COMMON/COMP2/C2,Mn2,Si2,Ni2,Cr2,Mo2,Cu2,W2,V2,Nb2,Ti2,Al2,N2
  1780.     COMMON/STRUCT/XA,XF,XP,XB,XM,GS,G,HV,XL
  1781.     COMMON/TEMPS/AC1_1,AC3_1,TDISS1,AE4_1,AE3_1,AE1_1,TBS1,TMS1,as1_1&  ,Tliquid1,Tsolid1
  1782.     COMMON/TEMPS2/AC1_2,AC3_2,TDISS2,AE4_2,AE3_2,AE1_2,TBS2,TMS2,as1_2& ,Tliquid2,Tsolid2
  1783.  
  1784.     DO 20 I=1,NNODE,1
  1785.          if ((HV(I) .lt. 1.0D0) .and. (VR(I) .gt. 1.0D0) &       .and. (TEMP2(I) .lt. TMS) .and. (XA(I) .le. 1.0d-6)) then
  1786.           HM=294.0D0+884.0D0*C-265.2D0*C**3.0D0
  1787.           CE=C+Si/24.D0+Mn/5.D0+Cu/10.D0+Ni/18.D0&         +Cr/5.D0+Mo/2.5D0+V/5.0D0+Nb/3.0D0
  1788.           IF (CE .LE. 0.75D0) THEN       
  1789.              HB=117.0D0+197.0D0*CE
  1790.           ELSE
  1791.              HB=145.0D0+130.D0*DTANH(2.65D0*CE-0.69D0)
  1792.           ENDIF      
  1793.           HVM=902.60D0*C+121.15D0+26.68D0*DLOG10(VR(I))
  1794.           HVB=-323.0D0+185.0D0*C+330.0D0*Si+153.0D0*Mn+66.0D0*Ni& +144.0D0*Cr+191.0D0*Mo+DLOG10(VR(I))*(89.0D0+53.D0*C&       -55.D0*Si-22.0D0*Mn-10.D0*Ni-20.D0*Cr-33.D0*Mo)
  1795.           HVFP=42.0D0+223.D0*C+53.D0*Si+30.D0*Mn+12.6D0*Ni+7.D0*Cr& +19.D0*Mo&  +DLOG10(VR(I))*(10.D0-19.D0*Si+4.D0*Ni+8.D0*Cr+130.D0*V)
  1796.           IF(HVM .GT. HM) THEN
  1797.              HVM=HM
  1798.           ENDIF
  1799.           IF(HVB .GT. HB) THEN
  1800.             HVB=HB
  1801.           ENDIF
  1802.           HV(I)=XM(I)*HVM+XB(I)*HVB+(XF(I)+XP(I))*HVFP
  1803.          endif
  1804. 20  CONTINUE
  1805.  
  1806.     RETURN
  1807.     END
  1808.     SUBROUTINE HARDNESS2(XCOMPOS,XF,XP,XB,XM,VR,HV)
  1809.     IMPLICIT REAL*8 (A-H,O-Z)
  1810.     real*8 xcompos(*)
  1811.     REAL*8 CE,HM,HB,HVM,HVB,HVFP
  1812.     REAL*8 VR,HV
  1813.     REAL*8 XF,XP,XB,XM
  1814.       C = XCOMPOS(1)
  1815.     Mn = XCOMPOS(2)
  1816.     Si = XCOMPOS(3)
  1817.     Ni = XCOMPOS(4)
  1818.     Cr = XCOMPOS(5)
  1819.     Mo = XCOMPOS(6)
  1820.     Cu = XCOMPOS(7)
  1821.     W = XCOMPOS(8)
  1822.     V = XCOMPOS(9)
  1823.     Nb = XCOMPOS(10)
  1824.     Ti = XCOMPOS(11)
  1825.     Al = XCOMPOS(12)
  1826.     N = XCOMPOS(13)
  1827.  
  1828.       HM=294.0D0+884.0D0*C-265.2D0*C**3.0D0
  1829.       CE=C+Si/24.D0+Mn/5.D0+Cu/10.D0+Ni/18.D0&         +Cr/5.D0+Mo/2.5D0+V/5.0D0+Nb/3.0D0
  1830.       IF (CE .LE. 0.75D0) THEN       
  1831.          HB=117.0D0+197.0D0*CE
  1832.       ELSE
  1833.          HB=145.0D0+130.D0*DTANH(2.65D0*CE-0.69D0)
  1834.       ENDIF      
  1835.       HVM=902.6D0*C+121.15D0+26.68D0*DLOG10(VR)
  1836.       HVB=-323.D0+185.0D0*C+330.0D0*Si+153.0D0*Mn+66.0D0*Ni&    +144.0D0*Cr+191.0D0*Mo+DLOG10(VR)*(89.0D0+53.D0*C&    -55.D0*Si-22.0D0*Mn-10.D0*Ni-20.D0*Cr-33.D0*Mo)
  1837.     HVFP=42.D0+223.D0*C+53.D0*Si+30.D0*Mn+12.6D0*Ni+7.D0*Cr+19.D0*Mo &    +DLOG10(VR)*(10.D0-19.D0*Si+4.D0*Ni+8.D0*Cr+130.D0*V)
  1838.     IF(HVM .GT. HM) THEN
  1839.        HVM=HM
  1840.     ENDIF
  1841.     IF(HVB .GT. HB) THEN
  1842.        HVB=HB
  1843.     ENDIF
  1844.     HV=XM*HVM+XB*HVB+(XF+XP)*HVFP
  1845.     RETURN
  1846.     END
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!