Powerstation上遇到的浮点数错误 - Fortran Coder

[Fortran] 纯文本查看 复制代码
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