[编译器] Powerstation上遇到的浮点数错误

PROGRAM TTMEHL 
DIMENSION X(0:120),D(0:120,0:120),BSA(-5:-1,-1:1),BSB(11:15,-1:1),SAB(0:120),P(0:120),H(0:120),TEMF(0:120,-5:15),TEMP(0:120,0:10),DSABDX(0:120),Z(0:10),ROR(0:120,0:10),U(0:120,0:10),DUDZ(0:120,0:10),TAU(0:120,0:10),AITA(0:120,0:10),DPDX(0:120),DPDT(0:120),RHAL(0:120,0:10),RHA(0:120,0:10),ACOE(-1:1,119),BCOE(-1:1,119),TEML(0:120,-5:15),PL(0:120),HL(0:120),C(119,119),SKK(0:120),DENE(0:120),DENL(0:120),ZA(-6:0),ZB(10:16),EPSI(0:120),DENS(0:120),PJUG(0:120),FRIC(64),HAVE(64),HH0(64),PREM(0:120,64),P16(0:120),P32(0:120),P48(0:120),H16(0:120),H32(0:120),H48(0:120),S16(0:120),S32(0:120),S48(0:120),T16(0:120,-5:15),T32(0:120,-5:15),T48(0:120,-5:15)

  OPEN(UNIT=15,FILE='RES.DAT',STATUS='UNKNOWN')

  PI=4.0*ATAN(1.0)
  N=120
  M=110

  SLIP=1.2
  AA=0.10
  AB=0.05
  ALG=0.4
  BLG=0.2

  W0=3.0E-5
  U0=1.5E-11
  G0=5000.0
  R0=0.027
  T0=313.0
  TAU0=1.0E+7
  AITA0=0.08
  Z0=0.6
  S0=1.1
  ALFA=2.19E-8

  HK0=0.38
  CALL MESH(X,N)
  DZ=0.1
  DO 20 J=0,10
  Z(J)=DZ*J
 20 CONTINUE

  E1=G0/ALFA
  PH=E1/4.0*SQRT(8.0*W0/PI)
  BH=R0*SQRT(8.0*W0/PI)
  LSTEP=64

  WRITE(15,40)
 40 FORMAT(1X,15('*'), 'Eyring Transient Micro-EHL Solution',15('*')/)
  WRITE(15,60) W0,U0,G0
 60 FORMAT(5X,'W0=',E13.6,10X,'U0=',E13.6,10X,'G0=',F8.2)
  WRITE(15,70) PH,BH,SLIP
 70 FORMAT(5X,'PH=',E13.6,10X,'BH=',E13.6,10X,'SLIP=',F6.3)

  TAU00=TAU0/PH
  ALAMDA=0.75*U0*PI*PI/(W0*W0)
  AL1=LOG(AITA0)+9.67
  AL2=5.1*1E-9*PH
  AL3=1.0/(1.0-138.0/T0)
  AL4=(138.0/T0)/(1.0-138.0/T0)
  CL1=0.6*1.0E-9*PH
  CL2=1.7*1.0E-9*PH
  CL3=0.00065*T0

  H0=BH*BH/R0
  DZA=0.05
  DZB=0.05

  UA=(2.0+SLIP)/2.0
  UB=(2.0-SLIP)/2.0
  TLG=BLG/UB
  DT=TLG/LSTEP
  WRITE(15,75)AA,AB,ALG,BLG
 75 FORMAT(5X,'AA=',F10.5,5X,'AB=',F10.5,5X,'ALG=',F7.3,5X,'BLG=',F7.3/)

  CF0=2000.0
  CFA=470.0
  CFB=470.0
  DEN0=870.0
  DENA=7850.0
  DENB=7850.0
  FK0=0.14
  FKA=46.0
  FKB=46.0

  UR=U0*E1*R0/AITA0
  PRDL=CF0*AITA0/FK0
  REYN=DEN0*UR*BH/AITA0
  HOB=H0/BH
  RESTAR=REYN*HOB**2
  EC=UR*UR/(CF0*T0)
  YT=PH*CF0*BH/(UR*FK0)

  CNA=BH*UR*CFA*DENA/FKA
  CNB=BH*UR*CFB*DENB/FKB
  CMA=(FK0/FKA)/HOB
  CMB=(FK0/FKB)/HOB
  CH=(PH*H0)/(AITA0*UR)
  CNS1=RESTAR*PRDL
  CNS2=CL3*EC*YT*HOB*HOB
  CNS3=EC*YT*HOB

  CALL DEFORM(X,D,PI,N)
  CALL COEFCL(ACOE,BCOE,SKK,X,N)
  CALL TCOEFF(BSA,BSB,ZA,ZB)

  DO 90 I=1,N-1
  DO 80        J=1,N-1
  C(I,J)=ACOE(-1,I)*D(I-1,J)+ACOE(0,I)*D(I,J)+ACOE(1,I)*D(I+1,J)
 80 CONTINUE
 90 CONTINUE

  TIME=0.0
  CALL SABCAL(SAB,DSABDX,X,AA,AB,ALG,BLG,UA,UB,TIME,PI,N)
  CALL PREGET(PL,X,SAB,SKK,PI,HK0,N)
  CALL TEMGET(TEML,PL,X,Z,ZA,ZB,SLIP,N)

  DO 110 I=0,N
  P(I)=PL(I)
  DO 100 J=0,10
  TEMP(I,J)=TEML(I,J)
 100 CONTINUE
 110 CONTINUE

  DO 125 I=0,N
  DO 115 J=-5,15
  TEMF(I,J)=TEML(I,J)
 115 CONTINUE
 125 CONTINUE

  DE0=0.0
  DO 130 J=0,N
  DE0=DE0+D(70,J)*PL(J)
 130 CONTINUE
  H00=HK0-DE0

  CALL HCAL(D,PL,X,HL,SAB,H00,N)
  CALL ARCAL(AITA,ROR,TEMP,PL,AL1,AL2,AL3,AL4,CL1,CL2,CL3,Z0,S0,N)
  CALL PARTFX(PL,DPDX,ACOE,X,N)
  CALL TAUCAL(DPDX,HL,Z,TAU,AITA,TAU00,UA,UB,CH,HOB,DZ,N)
  CALL RHEOCL(EPSI,DENS,DENL,TAU,HL,AITA,ROR,Z,TAU00,DZ,UA,UB,ALAMDA,N)

  DO 150 I=0,N
  RHAL(I,0)=0.0
  DO 140 J=1,10
  RHAL(I,J)=RHAL(I,J-1)+0.5*DZ*(ROR(I,J-1)+ROR(I,J))
 140 CONTINUE
 150 CONTINUE
  DO 170 I=0,N
  DO 160 J=0,10
  RHAL(I,J)=RHAL(I,J)*HL(I)
 160 CONTINUE
 170 CONTINUE

  IJUDGE=0
  DO 2000 LARGE=1,20

  WRITE(*,*)
  WRITE(*,260) LARGE
 260 FORMAT(//5X,30('*'),1X,'LARGE LOOP',I3,1X,30('*')//)

  DO 280 I=0,N
  PJUG(I)=PL(I)
 280 CONTINUE

  DO 1000 LTIME=1,LSTEP

  WRITE(*,*) 'LTIME=',LTIME
  TIME=DT*LTIME
  CALL SABCAL(SAB,DSABDX,X,AA,AB,ALG,BLG,UA,UB,TIME,PI,N)

  IF (IJUDGE.EQ.0) THEN

  CALL PRESS(P,H,HL,TEMP,D,X,SAB,AITA,ROR,DPDX,ACOE,BCOE,EPSI,Z,TAU,DENE,DENS,DENL,DSABDX,SKK,C,M,H00,AL1,AL2,AL3,AL4,CL1,CL2,CL3,Z0,S0,TAU00,UA,UB,CH,HOB,DZ,DT,ALAMDA,PI,N)

  DO 360 I=0,N
  PREM(I,LTIME)=P(I)
 360 CONTINUE

  CALL PARTFX(P,DPDX,ACOE,X,N)
  DO 370 I=0,N
  DPDT(I)=(P(I)-PL(I))/DT
 370 CONTINUE
  CALL HCAL(D,P,X,H,SAB,H00,N)

  CALL TEMPER(TEMF,TEMP,BSA,BSB,ROR,U,DUDZ,P,H,TAU,AITA,DPDX,DPDT,RHAL,RHA,ACOE,TEML,Z,X,CNA,CNB,CMA,CMB,CNS1,CNS2,CNS3,DT,DZA,DZB,UA,UB,CH,HOB,DZ,AL1,AL2,AL3,AL4,CL1,CL2,CL3,Z0,S0,TAU00,N)

  ELSE

  DO 380 I=0,N
  P(I)=PREM(I,LTIME)
 380 CONTINUE

  CALL PARTFX(P,DPDX,ACOE,X,N)
  DO 390 I=0,N
  DPDT(I)=(P(I)-PL(I))/DT
 390 CONTINUE
  CALL HCAL(D,P,X,H,SAB,H00,N)

  CALL TEMPER(TEMF,TEMP,BSA,BSB,ROR,U,DUDZ,P,H,TAU,AITA,DPDX,DPDT,RHAL,RHA,ACOE,TEML,Z,X,CNA,CNB,CMA,CMB,CNS1,CNS2,CNS3,DT,DZA,DZB,UA,UB,CH,HOB,DZ,AL1,AL2,AL3,AL4,CL1,CL2,CL3,Z0,S0,TAU00,N)
  CALL PRESS(P,H,HL,TEMP,D,X,SAB,AITA,ROR,DPDX,ACOE,BCOE,EPSI,Z,TAU,DENE,DENS,DENL,DSABDX,SKK,C,M,H00,AL1,AL2,AL3,AL4,CL1,CL2,CL3,Z0,S0,TAU00,UA,UB,CH,HOB,DZ,DT,ALAMDA,PI,N)

  DO 400 I=0,N
  PREM(I,LTIME)=P(I)
 400 CONTINUE

  END IF

  CALL HAVFRC(H,TAU,SKK,PI,HAV,FRC,N)
  HAVE(LTIME)=HAV
  FRIC(LTIME)=FRC
  HH0(LTIME)=H00

  DO 410 I=0,N
  PL(I)=P(I)
  HL(I)=H(I)
 410 CONTINUE

  DO 440 I=0,N
  DO 430 J=-5,15
  TEML(I,J)=TEMF(I,J)
 430 CONTINUE
 440 CONTINUE

  DO 460 I=0,N
  DO 450 J=0,10
  RHAL(I,J)=RHA(I,J)
 450 CONTINUE
 460 CONTINUE

  IF (LTIME.EQ.16) THEN
  DO 520 I=0,N
  P16(I)=P(I)
  H16(I)=H(I)
  S16(I)=SAB(I)
  DO 510 J=-5,15
  T16(I,J)=TEMF(I,J)
 510 CONTINUE
 520 CONTINUE
  ELSE IF (LTIME.EQ.32) THEN
  DO 540 I=0,N
  P32(I)=P(I)
  H32(I)=H(I)
  S32(I)=SAB(I)
  DO 530 J=-5,15
  T32(I,J)=TEMF(I,J)
 530 CONTINUE
 540 CONTINUE
  ELSE IF (LTIME.EQ.48) THEN
  DO 560 I=0,N
  P48(I)=P(I)
  H48(I)=H(I)
  S48(I)=SAB(I)
  DO 550 J=-5,15
  T48(I,J)=TEMF(I,J)
 550 CONTINUE
 560 CONTINUE
  END IF

  CALL ARCAL(AITA,ROR,TEMP,PL,AL1,AL2,AL3,AL4,CL1,CL2,CL3,Z0,S0,N)
  CALL TAUCAL(DPDX,HL, Z,TAU,AITA,TAU00,UA,UB,CH,HOB,DZ,N)
  CALL RHEOCL(EPSI,DENS,DENL,TAU,HL,AITA,ROR,Z,TAU00,DZ,UA,UB,ALAMDA,N)

 1000 CONTINUE

  ERROR=0.0
  SS=0.0
  DO 1020 I=0,N
  ERROR=ERROR+ABS(P(I)-PJUG(I))
  SS=SS+P(I)
 1020 CONTINUE
  ERROR=ERROR/SS
  
WRITE(*,*)
  WRITE(*,*) 'LARGE=',LARGE,'ERROR=',ERROR
  WRITE(15,*) 'LARGE=',LARGE,'ERROR=',ERROR

  IF (ERROR.LT.0.05) IJUDGE=1
  IF ((ERROR.LT.0.002).AND.(LARGE.GT.10)) GOTO 2010
 2000 CONTINUE

 2010 WRITE(15,2020)
 2020 FORMAT(//3X,11('-'),'The Final Results',10('-')/)
      WRITE (15,2040)
 2040 FORMAT (/3X,10('*'),'Friction Coefficient',10('*')/)
      WRITE (15,'(6(1X,I4,F8.5))') (I,FRIC(I),I=1,LSTEP)
          WRITE(15,2050)
 2050 FORMAT(/3X,10('*'),'Average Film Thickness',10('*')/)
      WRITE(15,'(6(1X,I4,F8.5))') (I,HAVE(I),I=1,LSTEP)
          WRITE(15,2060)
 2060 FORMAT(/3X,10('*'),'Rigid Surface Distance H00',10('*')/)
      WRITE(15,'(6(1X,I4,F8.4))') (I,HH0(I),I=1,LSTEP)

      WRITE(15,2080)
 2080 FORMAT(//1X,20('*'),'Results at Time Step 16',21('*')/)
      CALL OUTRES(P16,H16,S16,T16,X,N,T0)
          WRITE(15,2100)
 2100 FORMAT(//1X,20('*'),'Results at Time Step 32',21('*')/)
      CALL OUTRES(P32,H32,S32,T32,X,N,T0)
          WRITE(15,2120)
 2120 FORMAT(//1X,20('*'),'Results at Time Step 48',21('*')/)
      CALL OUTRES(P48,H48,S48,T48,X,N,T0)
          WRITE(15,2150)
 2150 FORMAT(//1X,20('*'),'Results at Time Step 64',21('*')/)
      CALL OUTRES(P,H,SAB,TEMF,X,N,T0)
          WRITE(15,2200)
 2200 FORMAT(//30X,'end of output file')
      END


      SUBROUTINE PRESS(P,H,HL,TEMP,D,X,SAB,AITA,ROR,DPDX,ACOE,BCOE,EPSI,Z,TAU,DENE,DENS,DENL,DSABDX,SKK,C,M,H00,AL1,AL2,AL3,AL4,CL1,CL2,CL3,Z0,S0,TAU00,UA,UB,CH,HOB,DZ,DT,ALAMDA,PI,N)

      DIMENSION P(0:120),POLD(0:120),H(0:120),C(119,119),HL(0:120),TEMP(0:120,0:10),D(0:120,0:120),X(0:120),SAB(0:120),AITA(0:120,0:10),ROR(0:120,0:10),DPDX(0:120),ACOE(-1:1,119),BCOE(-1:1,119),Z(0:10),TAU(0:120,0:10),EPSI(0:120),DENS(0:120),DENE(0:120),DENL(0:120),DXEPSI(0:120),DXDENS(0:120),E(-1:1,119),F(0:120),SK(0:120,0:120),SY(0:120),DSABDX(0:120),SKK(0:120),PRH0(0:120),SA(0:120,0:120),SB(0:120),DP(0:120),IIMX(0:120),JJMX(0:120),ZZR(0:120),ZR(0:120,0:120),POLDK(0:120)

      DO 600 KKK=1,10

      DO 20 I=0,N
          POLDK(I)=P(I)
 20 CONTINUE
      ERRPNL=1.0

      DO 500 NNN=1,10

     DO 30 I=0,N
          POLD(I)=P(I)
 30 CONTINUE
      HOLD=H00

      CALL HCAL(D,P,X,H,SAB,H00,N)
          CALL ARCAL(AITA,ROR,TEMP,P,AL1,AL2,AL3,AL4,CL1,CL2,CL3,Z0,S0,N)
          CALL PARTFX(P,DPDX,ACOE,X,N)
          CALL TAUCAL(DPDX,H,Z,TAU,AITA,TAU00,UA,UB,CH,HOB,DZ,N)
          CALL RHEOCL(EPSI,DENS,DENE,TAU,H,AITA,ROR,Z,TAU00,DZ,UA,UB,ALAMDA,N)
          CALL PARTFX(EPSI,DXEPSI,ACOE,X,N)
          CALL PARTFX(DENS,DXDENS,ACOE,X,N)

      DO 50 I=1,M-1
          E(-1,I)=DXEPSI(I)*ACOE(-1,I)+EPSI(I)*BCOE(-1,I)
          E(0,I)=DXEPSI(I)*ACOE(0,I)+EPSI(I)*BCOE(0,I)
          E(1,I)=DXEPSI(I)*ACOE(1,I)+EPSI(I)*BCOE(1,I)
 50 CONTINUE

      DO 60 I=1,M-1
      F(I)=DXDENS(I)+(2.0*DENE(I)-DENL(I))/DT
 60 CONTINUE
      DO 80 I=1,M-1
          DO 70 J=1,M-1
          SK(I,J)=-(F(I)*D(I,J)+DENS(I)*C(I,J))
 70 CONTINUE
 80 CONTINUE

      DO 90 I=1,M-1
          SK(I,I-1)=SK(I,I-1)+E(-1,I)
          SK(I,I)=SK(I,I)+E(0,I)
          SK(I,I+1)=SK(I,I+1)+E(1,I)
 90 CONTINUE

      DO 100 I=1,M-1
          SK(I,M)=-F(I)
          SY(I)=F(I)*(0.5*X(I)**2-SAB(I))+DENS(I)*(X(I)-DSABDX(I))-DENE(I)*HL(I)/DT 
 100 CONTINUE

      DO 130 J=1,M-1
          SK(M,J)=SKK(J)
 130 CONTINUE
      SK(M,M)=0.0
          SY(M)=PI/2.0

      DO 140 I=0,M-1
          PRH0(I)=P(I)
 140 CONTINUE
      PRH0(M)=H00
          DO 150 I=M+1,N
          PRH0(I)=0.0
 150 CONTINUE

      IF (NNN.EQ.1) THEN
          DO 170 I=1,M
          DO 160 J=1,M
          SA(I,J)=SK(I,J)
 160 CONTINUE
 170 CONTINUE
      END IF

      DO 200 I=1,M
          SB(I)=SY(I)
          DO 190 J=1,M
          SB(I)=SB(I)-SK(I,J)*PRH0(J)
 190 CONTINUE
 200 CONTINUE

      IF (NNN.EQ.1) THEN
          CALL GAUSSP(SA,SB,DP,IIMX,JJMX,ZZR,ZR,M)
          ELSE
          CALL GAUSSP(SA,SB,DP,IIMX,JJMX,ZZR,ZR,M)
          END IF
     DO 210 I=1,M
     PRH0(I)=PRH0(I)+DP(I)
 210 CONTINUE

      SSS=0.0
          SER=0.0
          DO 270 I=1,M-1
          SER=SER+ABS(DP(I))
          SSS=SSS+PRH0(I)
 270 CONTINUE
      ERRPN=SER/SSS

      OMIGA=0.2
          DO 290 I=1,M-1
          P(I)=POLD(I)+OMIGA*DP(I)
 290 CONTINUE
      H00=HOLD+OMIGA*DP(M)

      IF (ERRPN.GT.ERRPNL) THEN
          WRITE(*,*) 'BAD! NNN=',NNN,'ERRPN=',ERRPN
          GOTO 510
          END IF

      DO 310 I=0,N
          IF (P(I).LT.0.0) P(I)=0.0
 310 CONTINUE
 
      IF ((ERRPN.LT.0.001).AND.(NNN.GT.5)) GOTO 510
          ERRPNL=ERRPN

 500 CONTINUE

 510  SSK=0.0
      SEK=0.0
          DO 520 I=1,M-1
          SSK=SSK+P(I)
          SEK=SEK+ABS(P(I)-POLDK(I))
 520 CONTINUE
      ERRPK=SEK/SSK
          CALL GETM(P,M,N)
          WRITE(*,*) 'KKK=',KKK,'ERRPK=',ERRPK,'M=',M
          IF ((ERRPK.LT.0.001).AND.(KKK.GT.1)) GOTO 610
 600 CONTINUE
 610 RETURN
     END


      SUBROUTINE TEMPER(TEMF,TEMP,BSA,BSB,ROR,U,DUDZ,P,H,TAU,AITA,DPDX,DPDT,RHAL,RHA,ACOE,TEML,Z,X,CNA,CNB,CMA,CMB,CNS1,CNS2,CNS3,DT,DZA,DZB,UA,UB,CH,HOB,DZ,AL1,AL2,AL3,AL4,CL1,CL2,CL3,Z0,S0,TAU00,N)
          DIMENSION TEMF(0:120,-5:15),TEMP(0:120,0:10),TK(21),BSA(-5:-1,-1:1),BSB(11:15,-1:1),P(0:120),H(0:120),ROR(0:120,0:10),U(0:120,0:10),Z(0:10),DUDZ(0:120,0:10),TAU(0:120,0:10),X(0:120),AITA(0:120,0:10),DPDX(0:120),RHAL(0:120,0:10),RHA(0:120,0:10),RH(0:120,0:10),ACOE(-1:1,119),TOLD(0:120,-5:15),STA(-6:16,-6:16),STB(-6:16),SA(21,21),SB(21),TEML(0:120,-5:15),DPDT(0:120)
          J1=5
          DO 500 LMN=1,20

      DO 60 I=0,N
          DO 50 J=-5,15
          TOLD(I,J)=TEMF(I,J)
 50 CONTINUE
 60 CONTINUE
    
      DO 80 I=0,N
          DO 70 J=0,10
          TEMP(I,J)=TEMF(I,J)
 70 CONTINUE
 80 CONTINUE

      CALL ARCAL(AITA,ROR,TEMP,P,AL1,AL2,AL3,AL4,CL1,CL2,CL3,Z0,S0,N)
          CALL TAUCAL(DPDX,H,Z,TAU,AITA,TAU00,UA,UB,CH,HOB,DZ,N)
          CALL UCAL(U,DUDZ,TAU,AITA,H,CH,UA,UB,TAU00,DZ,N)
          CALL RHABCL(RH,RHA,RHAL,H,ROR,U,ACOE,DZ,DT,N)

      DO 300 K=1,N-1

      DXP=X(K)-X(K-1)
          DXM=X(K+1)-X(K)
          HDZ=H(K)*DZ

      DO 110 I=-6,16
          STB(I)=0.0
          DO 100 J=-6,16
          STA(I,J)=0.0
 100 CONTINUE
 110 CONTINUE

      DO 125 J=-5,-1
          STA(J,J-1)=-BSA(J,-1)
          STA(J,J)=CNA/DT+CNA*UA/DXP-BSA(J,0)
          STA(J,J+1)=-BSA(J,1)
          STB(J)=CNA/DT*TEML(K,J)+CNA*UA/DXP*TEMF(K-1,J)
 125 CONTINUE

      J=0
          STA(J,J-1)=1.0/DZA
          STA(J,J)=-(CMA/HDZ+1.0/DZA)
          STA(J,J+1)=CMA/HDZ
          STB(J)=0.0

      DO 130 J=1,9
          IF (U(K,J).GT.0.0) THEN
          STA(J,J-1)=-(1.0/HDZ**2-0.5*CNS1*RH(K,J)/HDZ)
          STA(J,J)=CNS1*ROR(K,J)/DT+CNS1*ROR(K,J)*U(K,J)/DXP+2.0/HDZ**2
          STA(J,J+1)=-(1.0/HDZ**2+0.5*CNS1*RH(K,J)/HDZ)          
          STB(J)=CNS3*TAU(K,J)/H(K)*DUDZ(K,J)+CNS1*ROR(K,J)/DT*TEML(K,J)+CNS1/DXP*ROR(K,J)*U(K,J)*TEMF(K-1,J)+CNS2/ROR(K,J)*(DPDT(K)+U(K,J)*DPDX(K))*TEMF(K,J)
          ELSE
          STA(J,J-1)=-(1.0/HDZ**2-0.5*CNS1*RH(K,J)/HDZ)
          STA(J,J)=CNS1*ROR(K,J)/DT-CNS1*ROR(K,J)*U(K,J)/DXM+2.0/HDZ**2
          STA(J,J+1)=-(1.0/HDZ**2+0.5*CNS1*RH(K,J)/HDZ)
          STB(J)=CNS3*TAU(K,J)/H(K)*DUDZ(K,J)+CNS1*ROR(K,J)/DT*TEML(K,J)-CNS1/DXM*ROR(K,J)*U(K,J)*TEMF(K+1,J)+CNS2/ROR(K,J)*(DPDT(K)+U(K,J)*DPDX(K))*TEMF(K,J)
          END IF
 130 CONTINUE 

      J=10
          STA(J,J-1)=CMB/HDZ
          STA(J,J)=-(CMB/HDZ+1.0/DZB)
          STA(J,J+1)=1.0/DZB
          STB(J)=0.0

      DO 140 J=11,15
          STA(J,J-1)=-BSB(J,-1)
          STA(J,J)=CNB/DT+CNB*UB/DXP-BSB(J,0)
          STA(J,J+1)=-BSB(J,1)
          STB(J)=CNB/DT*TEML(K,J)+CNB*UB/DXP*TEMF(K-1,J)
 140 CONTINUE

      J=-5
          ABCD=STA(J,J-1)*1.0
          STB(J)=STB(J)-ABCD
          J=15
          ABCD=STA(J,J+1)*1.0
          STB(J)=STB(J)-ABCD

      DO 170 I=1,21
          SB(I)=STB(I-6)
          DO 160 J=1,21
          SA(I,J)=STA(I-6,J-6)
 160 CONTINUE
 170 CONTINUE

      CALL CATCHT(SA,SB,TK)
          DO 180 J=1,21
          IF(TK(J).LT.1.0) TK(J)=1.0
          TEMF(K,J-6)=TK(J)
 180 CONTINUE

 300 CONTINUE

      DO 310 J=-5,15
          IF ((J.GT.0).AND.(J.LT.10)) THEN
           IF (U(0,J).GE.0.0) THEN
           TEMF(0,J)=1.0
           ELSE
           TEMF(0,J)=TEMF(1,J)
           END IF
          ELSE
          TEMF(0,J)=1.0
          END IF
          TEMF(N,J)=TEMF(N-1,J)
 310 CONTINUE

      TERR=0.0
          TSSS=0.0
          DO 320 I=1,N-1
          TERR=TERR+ABS(TEMF(I,J1)-TOLD(I,J1))
          TSSS=TSSS+ABS(TEMF(I,J1))
 320 CONTINUE
      TERR=TERR/TSSS
           WRITE(*,*)'---LMN=',LMN,'---TERR=',TERR

      OMIT=0.4
          DO 340 I=1,N
          DO 330 J=-5,15
          TEMF(I,J)=TOLD(I,J)+OMIT*(TEMF(I,J)-TOLD(I,J))
          IF (TEMF(I,J).LT.1.0) TEMF(I,J)=1.0
 330 CONTINUE
 340 CONTINUE

      DO 360 I=0,N
          DO 350 J=0,10
          TEMP(I,J)=TEMF(I,J)
 350 CONTINUE
 360 CONTINUE

     IF ((TERR.LT.0.0001).AND.(LMN.GT.4)) GOTO 600

 500 CONTINUE
 600  CALL ARCAL(AITA,ROR,TEMP,P,AL1,AL2,AL3,AL4,CL1,CL2,CL3,Z0,S0,N)
      CALL TAUCAL(DPDX,H,Z,TAU,AITA,TAU00,UA,UB,CH,HOB,DZ,N)
          CALL UCAL(U,DUDZ,TAU,AITA,H,CH,UA,UB,TAU00,DZ,N)
          CALL RHABCL(RH,RHA,RHAL,H,ROR,U,ACOE,DZ,DT,N)
          END


      SUBROUTINE GAUSSP(SA,SB,P,IIMX,JJMX,ZZR,ZR,N)
          DIMENSION SA(0:120,0:120),SB(0:120),P(0:120),IIMX(0:120),JJMX(0:120),ZZR(0:120),ZR(0:120,0:120)
          M=N
          DO 50 K=1,M
          JJMX(K)=K
 50 CONTINUE

      DO 300 K=1,M-1
          ABSMX=0.0
          IMAX=K
          JMAX=K
          DO 125 J=K,M
          DO 110 I=K,M
          ABSA=ABS(SA(I,J))
           IF (ABSA.GT.ABSMX) THEN
           ABSMX=ABSA
           IMAX=I
           JMAX=J
           END IF
 110 CONTINUE
 125 CONTINUE
      IIMX(K)=IMAX

      DO 130 J=K,M
          TEMP=SA(K,J)
          SA(K,J)=SA(IMAX,J)
          SA(IMAX,J)=TEMP
 130 CONTINUE
      TEMP=SB(K)
          SB(K)=SB(IMAX)
          SB(IMAX)=TEMP
          DO 135 I=1,M
          TEMP=SA(I,K)
          SA(I,K)=SA(I,JMAX)
          SA(I,JMAX)=TEMP
 135 CONTINUE
      IT=JJMX(K)
          JJMX(K)=JJMX(JMAX)
          JJMX(JMAX)=IT

      ZZ=SA(K,K)
          ZZR(K)=ZZ
          DO 140 J=K,M
          SA(K,J)=SA(K,J)/ZZ
 140 CONTINUE
      SB(K)=SB(K)/ZZ
          DO 200 I=K+1,M
          Z=SA(I,K)
          ZR(K,I)=Z
          DO 150 J=K,M
          SA(I,J)=SA(I,J)-SA(K,J)*Z
 150 CONTINUE
      SB(I)=SB(I)-SB(K)*Z
 200 CONTINUE
 300 CONTINUE

      P(M)=SB(M)/SA(M,M)
          DO 500 K=M-1,1,-1
          SUM=SB(K)
          DO 400 J=K+1,M
          SUM=SUM-P(J)*SA(K,J)
 400 CONTINUE
      P(K)=SUM
 500 CONTINUE
      DO 700 J=1,M
          ITRUE=JJMX(J)
          SB(ITRUE)=P(J)
 700 CONTINUE
      DO 800 J=1,M
          P(J)=SB(J)
 800 CONTINUE
      END


      SUBROUTINE GAUSSR(SA,SB,P,IIMX,JJMX,ZZR,ZR,N)
          DIMENSION SA(0:120,0:120),SB(0:120),P(0:120),IIMX(0:120),JJMX(0:120),ZZR(0:120),ZR(0:120,0:120)
          M=N
          DO 300 K=1,M-1
          TEMP=SB(K)
          SB(K)=SB(IIMX(K))
          SB(IIMX(K))=TEMP
          SB(K)=SB(K)/ZZR(K)
          DO 200 I=K+1,M
          SB(I)=SB(I)-SB(K)*ZR(K,I)
 200 CONTINUE
 300 CONTINUE
     P(M)=SB(M)/SA(M,M)
         DO 500 K=M-1,1,-1
         SUM=SB(K)
         DO 400 J=K+1,M
         SUM=SUM-P(J)*SA(K,J)
 400 CONTINUE
     P(K)=SUM
 500 CONTINUE
     DO 700 J=1,M
         ITRUE=JJMX(J)
         SB(ITRUE)=P(J)
 700 CONTINUE
     DO 800 J=1,M
         P(J)=SB(J)
 800 CONTINUE
     END


     SUBROUTINE GETM(P,M,N)
         DIMENSION P(0:120)
         M2=M-1
          IF (P(M2).GT.0.003) THEN
          M=M+1
          ELSE IF (P(M2).LE.1.0E-8) THEN
          M=M-1
          ELSE
          M=M
          END IF
         IF (M.GT.N) M=N
         END


     SUBROUTINE ARCAL(AITA,ROR,TEMP,P,AL1,AL2,AL3,AL4,CL1,CL2,CL3,Z0,S0,N)
         DIMENSION AITA(0:120,0:10),ROR(0:120,0:10),TEMP(0:120,0:10),P(0:120)
         DO 50 I=0,N
         EA=(1.0+AL2*P(I))**Z0
         FA=CL1*P(I)/(1.0+CL2*P(I))
         DO 40 J=0,10
         EB=(AL3*TEMP(I,J)-AL4)**(-S0)
         EC=AL1*(-1.0+EA*EB)
         FB=CL3*(TEMP(I,J)-1.0)
         AITA(I,J)=EXP(EC)
         ROR(I,J)=1.0+FA-FB
 40 CONTINUE
 50 CONTINUE
     END


     SUBROUTINE TAUCAL(DPDX,H,Z,TAU,AITA,TAU00,UA,UB,CH,HOB,DZ,N)
         DIMENSION DPDX(0:120),H(0:120),TAU(0:120,0:10),Z(0:10),AITA(0:120,0:10),VIS(0:10)
         UUAB=UA-UB
         DO 100 K=0,N
         PXT0=DPDX(K)*H(K)*HOB/TAU00

     DO 20 J=0,10
         VIS(J)=AITA(K,J)
 20 CONTINUE

     FF1=0.0
         FF2=0.0
         DO 30 J=0,9
         FF1=FF1+0.5*DZ*(TAU00/VIS(J)*COSH(Z(J)*PXT0)+TAU00/VIS(J+1)*COSH(Z(J+1)*PXT0))
         FF2=FF2+0.5*DZ*(TAU00/VIS(J)*SINH(Z(J)*PXT0)+TAU00/VIS(J+1)*SINH(Z(J+1)*PXT0))
 30 CONTINUE

     FF1=FF1*H(K)
         FF2=FF2*H(K)
         RRL=SQRT(UUAB**2+CH**2*(FF1**2-FF2**2))
         RRL0=(RRL-UUAB)/(CH*(FF1+FF2))
          IF (RRL0.LE.0.0) THEN
          TAUA0=-2.0
          ELSE
          TAUA0=LOG(RRL0)
          END IF
          TAU(K,0)=TAUA0*TAU00
          DO 50 J=1,10
          TAU(K,J)=TAU(K,0)+HOB*H(K)*Z(J)*DPDX(K)
 50 CONTINUE
 100 CONTINUE
      END


     SUBROUTINE UCAL(U,DUDZ,TAU,AITA,H,CH,UA,UB,TAU00,DZ,N)
         DIMENSION U(0:120,0:10),TAU(0:120,0:10),DUDZ(0:120,0:10),H(0:120),AITA(0:120,0:10)
         DO 20 I=0,N
         DO 10 J=0,10
         DUDZ(I,J)=CH*H(I)*(TAU00/AITA(I,J))* SINH(TAU(I,J)/TAU00)
 10 CONTINUE
 20 CONTINUE
     DO 80 I=0,N
         U(I,0)=UA
         U(I,10)=UB
         DO 70 J=1,9
         U(I,J)=U(I,J-1)+0.5*DZ*(DUDZ(I,J-1)+DUDZ(I,J))
 70 CONTINUE
 80 CONTINUE
     END


     SUBROUTINE RHABCL(RH,RHA,RHAL,H,ROR,U,ACOE,DZ,DT,N)
         DIMENSION RH(0:120,0:10),RHA(0:120,0:10),RHB(0:120,0:10),RHAL(0:120,0:10),H(0:120),ROR(0:120,0:10),U(0:120,0:10),ACOE(-1:1,119)
         DO 200 I=0,N
         RHA(I,0)=0.0
         RHB(I,0)=0.0
         DO 100 J=1,10
         RHA(I,J)=RHA(I,J-1)+0.5*DZ*(ROR(I,J-1)+ROR(I,J))
         RHB(I,J)=RHB(I,J-1)+0.5*DZ*(ROR(I,J-1)*U(I,J-1)+ROR(I,J)*U(I,J))
 100 CONTINUE
 200 CONTINUE
     DO 400 I=0,N
         DO 300 J=0,10
         RHA(I,J)=RHA(I,J)*H(I)
         RHB(I,J)=RHB(I,J)*H(I)
 300 CONTINUE
 400 CONTINUE
     DO 600 I=1,N-1
         DO 500 J=0,10
         RH(I,J)=(RHA(I,J)-RHAL(I,J))/DT+ACOE(-1,I)*RHB(I-1,J)+ACOE(0,I)*RHB(I,J)+ACOE(1,I)*RHB(I+1,J)
 500 CONTINUE
 600 CONTINUE
     DO 800 J=0,10
         RH(0,J)=RH(1,J)
         RH(N,J)=RH(N-1,J)
 800 CONTINUE
     END


     SUBROUTINE CATCHT(A,B,T)
         DIMENSION A(21,21),B(21),T(21),C(21),D(21)
         C(1)=A(1,2)/A(1,1)
         D(1)=B(1)/A(1,1)
         DO 100 K=2,20
         C(K)=A(K,K+1)/(A(K,K)-A(K,K-1)*C(K-1))
         D(K)=(B(K)-A(K,K-1)*D(K-1))/(A(K,K)-A(K,K-1)*C(K-1))
 100 CONTINUE
     T(21)=(B(21)-A(21,20)*D(20))/(A(21,21)-A(21,20)*C(20))
         DO 200 K=20,1,-1
         T(K)=D(K)-C(K)*T(K+1)
 200 CONTINUE
     END


     SUBROUTINE RHEOCL(EPSI,DENS,DENE,TAU,H,AITA,ROR,Z,TAU00,DZ,UA,UB,ALAMDA,N)
         DIMENSION EPSI(0:120),DENS(0:120),DENE(0:120),H(0:120),TAU(0:120,0:10),AITA(0:120,0:10),ROR(0:120,0:10),Z(0:10),VA(0:10),VA0(0:10),VA1(0:10)
         DO 100 K=0,N
         DO 15 J=0,10
         TT0=TAU(K,J)/TAU00
          IF (ABS(TT0).LT.0.001) THEN
          VA(J)=1.0/AITA(K,J)
          ELSE
          VA(J)=SINH(TT0)/(AITA(K,J)*TT0)
          END IF
 15 CONTINUE
    VA0(0)=0.0
        VA1(0)=0.0
        DO 40 I=1,10
        VA0(I)=0.0
        VA1(I)=0.0
        DO 30 J=0,I-1
        VA0(I)=VA0(I)+0.5*DZ*(VA(J)+VA(J+1))
        VA1(I)=VA1(I)+0.5*DZ*(VA(J)*Z(J)+VA(J+1)*Z(J+1))
 30 CONTINUE
 40 CONTINUE

    AIE=1.0/VA0(10)
        AIE1=1.0/VA1(10)

    DEN0=0.0
        DEN1=0.0
        DEN2=0.0
        DO 50 J=0,9
        DEN0=DEN0+0.5*DZ*(ROR(K,J)+ROR(K,J+1))
        DEN1=DEN1+0.5*DZ*(ROR(K,J)*VA0(J)+ROR(K,J+1)*VA0(J+1))
        DEN2=DEN2+0.5*DZ*(ROR(K,J)*VA1(J)+ROR(K,J+1)*VA1(J+1))
 50 CONTINUE
    EPSI(K)=12.0*(AIE*DEN1/AIE1-DEN2)*H(K)**3/ALAMDA
        DENS(K)=2.0*(DEN1*AIE*(UB-UA)+DEN0*UA)/(UA+UB)
        DENE(K)=DEN0
 100 CONTINUE
    END


    SUBROUTINE HCAL(D,P,X,H,SAB,H00,N)
        DIMENSION D(0:120,0:120),P(0:120),X(0:120),H(0:120),SAB(0:120),DE(0:120)
        DO 100 I=0,N
        DE(I)=0.0
        DO 90 J=0,N
        DE(I)=DE(I)+D(I,J)*P(J)
 90 CONTINUE
 100 CONTINUE
    DO 130 I=0,N
        H(I)=H00+X(I)*X(I)/2.0+DE(I)-SAB(I)
 130 CONTINUE
    END


    SUBROUTINE PARTFX(F,DFDX,ACOE,X,N)
        DIMENSION F(0:120),DFDX(0:120),ACOE(-1:1,119),X(0:120)
        DO 50 I=1,N-1
        DFDX(I)=ACOE(-1,I)*F(I-1)+ACOE(0,I)*F(I)+ACOE(1,I)*F(I+1)
 50 CONTINUE
    DFDX(0)=(F(1)-F(0))/(X(1)-X(0))
        DFDX(N)=(F(N)-F(N-1))/(X(N)-X(N-1))
        END


    SUBROUTINE TCOEFF(BSA,BSB,ZA,ZB)
        DIMENSION ZA(-6:0),ZB(10:16),BSA(-5:-1,-1:1),BSB(11:15,-1:1)
        ZA(-6)=-3.15
        ZA(-5)=-1.55
        ZA(-4)=-0.75
        ZA(-3)=-0.35
        ZA(-2)=-0.15
        ZA(-1)=-0.05
        ZA(0)=0.0
        ZB(10)=0.0
        ZB(11)=0.05
        ZB(12)=0.15
        ZB(13)=0.35
        ZB(14)=0.75
        ZB(15)=1.55
        ZB(16)=3.15
        DO 100 J=-5,-1
        BSA(J,-1)=2.0/((ZA(J-1)-ZA(J))*(ZA(J-1)-ZA(J+1)))
        BSA(J,0)=2.0/((ZA(J)-ZA(J-1))*(ZA(J)-ZA(J+1)))
        BSA(J,1)=2.0/((ZA(J+1)-ZA(J-1))*(ZA(J+1)-ZA(J)))
 100 CONTINUE
    DO 200 J=11,15
        BSB(J,-1)=2.0/((ZB(J-1)-ZB(J))*(ZB(J-1)-ZB(J+1)))
        BSB(J,0)=2.0/((ZB(J)-ZB(J-1))*(ZB(J)-ZB(J+1)))
        BSB(J,1)=2.0/((ZB(J+1)-ZB(J-1))*(ZB(J+1)-ZB(J)))
 200 CONTINUE
    END


    SUBROUTINE DEFORM(X,D,PI,N)
        DIMENSION X(0:120),D(0:120,0:120)
        DP=-0.5/PI
        DO 20 I=0,N
        DO 10 J=0,N
        D(I,J)=0.0
 10 CONTINUE
 20 CONTINUE

    DO 500 I=0,N
        DO 400 JJ=1,N/2
        J=2*JJ-1
        D1=(X(J-1)-X(J))*(X(J-1)-X(J+1))
        D2=(X(J)-X(J-1))*(X(J)-X(J+1))
        D3=(X(J+1)-X(J-1))*(X(J+1)-X(J))
        A1=2.0/D1
        A2=-(X(J)+X(J+1))/D1+A1*X(J)
        A3=2.0/D2
        A4=-(X(J+1)+X(J-1))/D2+A3*X(J)
        A5=2.0/D3
        A6=-(X(J)+X(J-1))/D3+A5*X(J)

    ZMIN=X(I)-X(J-1)
        ZMAX=X(I)-X(J+1)

    IF (ABS(ZMIN).LT.1.0E-6) THEN
        FM1=0.0
        FN2=0.0
        ELSE
        FM1=ZMIN**2*(LOG(ZMIN**2)-3.0)
        FN2=ZMIN**3*(LOG(ABS(ZMIN**3))-4.0)
        END IF

    IF (ABS(ZMAX).LT.1.0E-6) THEN
        FM2=0.0
        FN1=0.0
        ELSE
        FM2=ZMAX**2*(LOG(ZMAX**2)-3.0)
        FN1=ZMAX**3*(LOG(ABS(ZMAX**3))-4.0)
        END IF

    FM=FM1-FM2
        FN=FN1-FN2
        FK=FM*(X(I)-X(J))/2.0+2.0/9.0*FN

    DIJM1=DP*(A1*FK+A2*FM/2.0)
        DIJ=DP*(A3*FK+A4*FM/2.0)
        DIJP1=DP*(A5*FK+A6*FM/2.0)
        D(I,J-1)=D(I,J-1)+DIJM1
        D(I,J)=D(I,J)+DIJ
        D(I,J+1)=D(I,J+1)+DIJP1
 400 CONTINUE
 500 CONTINUE
     END


    SUBROUTINE OUTRES(P,H,SAB,TEMF,X,N,T0)
        DIMENSION P(0:120),H(0:120),TEMF(0:120,-5:15),X(0:120),T(0:120,-5:15),SAB(0:120)
        DO 30 I=0,N
        DO 20 J=-5,15
        T(I,J)=T0*(TEMF(I,J)-1.0)
 20 CONTINUE
 30 CONTINUE
    DO 50 I=0,N
        WRITE(15,40) I,X(I),P(I),H(I),SAB(I)
 40 FORMAT(1X,'I=',I4,4X,'X=',F8.4,4X,'P=',F8.4,4X,'H=',F8.4,4X,'Sa+Sb=',F8.4)
 50 CONTINUE
    WRITE(15,60)
 60 FORMAT(/1X,'**X**',18X,'**Film Temperature Rises (^C)**'/)
    DO 80 I=0,N
        WRITE(15,70) X(I),(T(I,J),J=0,10)
 70 FORMAT(1X,F7.4,2X,11F6.2)
 80 CONTINUE
    WRITE(15,100)
 100 FORMAT(/1X,'**X**',5X,'**Temp.Rises of solid A**',7X,'**Temp.Rises of Solid B**'/)
    DO 120 I=0,N
        WRITE(15,110) X(I),(T(I,J),J=-5,-1),(T(I,J),J=11,15)
 110 FORMAT(1X,F7.4,3X,5F6.2,'---',5F6.2)
 120 CONTINUE
     END


    SUBROUTINE HAVFRC(H,TAU,SKK,PI,HAV,FRC,N)
        DIMENSION H(0:120),TAU(0:120,0:10),SKK(0:120)
        HAV=0.0
        DO 20 I=30,109
        HAV=HAV+0.0125*(H(I)+H(I+1))
 20 CONTINUE
    HAV=HAV/2.0
        FRC=0.0
        DO 40 I=0,N
        FRC=FRC+SKK(I)*TAU(I,0)
 40 CONTINUE
    FRC=-FRC/(PI/2.0)
        END


    SUBROUTINE MESH(X,N)
        DIMENSION X(0:120)
        DO 10,I=0,19
        X(I)=0.025*(1.0-1.16**(20-I))/0.16-1.25
 10 CONTINUE
    DO 20 I=20,N
        X(I)=0.025*(I-70)
 20 CONTINUE
    END


    SUBROUTINE SABCAL(SAB,DSABDX,X,AA,AB,ALG,BLG,UA,UB,TIME,PI,N)
        DIMENSION SAB(0:120),DSABDX(0:120),X(0:120)
        DO 10 I=0,N
        SAB(I)=AA*COS(2.0*PI*(X(I)-UA*TIME)/ALG)+AB*COS(2.0*PI*(X(I)-UB*TIME)/BLG)
        DSABDX(I)=-2.0*PI/ALG*AA*SIN(2.0*PI*(X(I)-UA*TIME)/ALG)-2.0*PI/BLG*AB*SIN(2.0*PI*(X(I)-UB*TIME)/BLG)
 10 CONTINUE
    END


    SUBROUTINE COEFCL(ACOE,BCOE,SKK,X,N)
        DIMENSION ACOE(-1:1,119),BCOE(-1:1,119),SKK(0:120),X(0:120),DL(0:119)
        DO 10 I=0,N-1
        DL(I)=X(I+1)-X(I)
 10 CONTINUE
    DO 20 I=1,N-1
        AB=DL(I-1)*(DL(I-1)+DL(I))
        ACOE(-1,I)=-DL(I)/AB
        BCOE(-1,I)=2.0/AB
        AB=DL(I-1)*DL(I)
        ACOE(0,I)=(DL(I)-DL(I-1))/AB
        BCOE(0,I)=-2.0/AB
        AB=DL(I)*(DL(I-1)+DL(I))
        ACOE(1,I)=DL(I-1)/AB
        BCOE(1,I)=2.0/AB
 20 CONTINUE
    SKK(0)=DL(0)/2.0
        SKK(N)=DL(N-1)/2.0
        DO 30 I=1,N-1
        SKK(I)=(DL(I)+DL(I-1))/2.0
 30 CONTINUE
    END


    SUBROUTINE PREGET(P,X,SAB,SKK,PI,HK0,N)
        DIMENSION P(0:120),X(0:120),SAB(0:120),SKK(0:120)
        X0=-1.4
        X1=-0.8
        X2=0.8
        X3=1.0
        PX1=SQRT(1.0-X1*X1)
        PX2=SQRT(1.0-X2*X2)
        DO 10 I=0,N
         IF (X(I).LE.X0) THEN
         P(I)=0.0
         ELSE IF(X(I).LE.X1) THEN
         P(I)=((X(I)-X0)/(X1-X0))**2*PX1
         ELSE IF (X(I).LE.X2) THEN
         P(I)=SQRT(1.0-X(I)*X(I))
         ELSE IF (X(I).LE.X3) THEN
         P(I)=((X(I)-X3)/(X2-X3))**2*PX2
         ELSE 
         P(I)=0.0
         END IF
 10 CONTINUE
    FC=0.75
        DO 20 I=0,N
        P(I)=P(I)*(1.0+FC*SAB(I)/HK0)
 20 CONTINUE
    W=0.0
        DO 30 I=0,N
        W=W+SKK(I)*P(I)
 30 CONTINUE
    WFC=0.5*PI/W
        DO 40 I=0,N
        P(I)=P(I)*WFC
 40 CONTINUE
    END


    SUBROUTINE TEMGET(T,P,X,Z,ZA,ZB,SLIP,N)
        DIMENSION T(0:120,-5:15),P(0:120),X(0:120),Z(0:10),ZA(-6:0),ZB(10:16),TM(0:120)
        X1=-1.4
        X2=X(N)
        ZM=0.5+0.1*SLIP**2
        E=0.2
        TOUT=0.06
        C=0.4
        FK=TOUT/(X2-X1)
        DO 20 I=0,N
        IF (X(I).LE.X1) THEN
        TM(I)=1.0
        ELSE
        TM(I)=E*SLIP*(P(I)+FK*(X(I)-X1))+1.0
        END IF
 20 CONTINUE
    DO 50 I=0,N
        IF (X(I).LE.X1) THEN
        DO 30 J=0,10
        T(I,J)=1.0
 30 CONTINUE
    ELSE
        A=C*SLIP*(TM(I)-1.0)/((C*SLIP+1.0)*ZM*ZM-(1.0-ZM)**2)
        DO 40 J=0,10
        T(I,J)=TM(I)-A*(Z(J)-ZM)**2
 40 CONTINUE
    END IF 
 50 CONTINUE

    DO 80 I=0,N
        G=(T(I,0)-1.0)/ZA(-5)**2
        DO 60 J=-5,-1
        T(I,J)=G*(ZA(J)-ZA(-5))**2+1.0
 60 CONTINUE
    G=(T(I,10)-1.0)/ZB(15)**2
        DO 70 J=11,15
        T(I,J)=G*(ZB(J)-ZB(15))**2+1.0
 70 CONTINUE
 80 CONTINUE
    END