$Debug
C     PROGRAM PENT(D,SAIDBO,TAPE5=D,TAPE6=SAIDBO,OUT1,TAPE7=OUT1,
C    1OUT2,TAPE10=OUT2,GR,TAPE8=GR)
C**** ENTWURF EINER ZWEIDIMENSIONALEN ODER ROTATIONSSYMETRISCHEN
C**** UEBERSCHALLWINDKANALDUESE, LFD-NAME : CHARA 14.01.76

      character*128 :: Char
      
      INTEGER RR,PNF,PN0,PNKL
      REAL MF,MD,MC,MS,MIS,MA,MAX,MREF,MISO,MSEST,MPM,M0,MCX,
     *MYREF,NYREF,NUEA,NUE,MKL,MEX
      EXTERNAL SPLINT
      COMMON /A/PNF(150,2),XF(150,2),YF(150,2),THF(150,2),MF(150,2),
     *NI(150,2)
     */C/DEL1(200)
     */D/XC(200),YC(200),YCX(200)
     */F/MC(200),TC(200),PC(200),CPC(200)
     */K/SPX(200),SPF(200)
     */I/DF(200)
      COMMON /G/A(200),B(200),C(200),D(200)
      COMMON /H/D121,D101,D241,D221,D201,D202,D203,D361,D341,D321,D301,
     *D342,D322,D302,D323,D303,D304
      COMMON /J/A131,A111,A251,A231,A211,A232,A212,A371,A351,A331,A311,
     *A352,A332,A312,A333,A313,A314
      COMMON /XYB/XB(150),YB(150),IKL
      COMMON /FTE/YXE(150),AMY(150),THY(150),IMTY
      COMMON /Q/ASPL(150,3),XSG(150),YSG(150),ISG
      COMMON /XX/PNKL(150,50),XKL(150,50),YKL(150,50),MKL(150,50),
     .THKL(150,50),NKL(50)
      COMMON /THTH/TH01,TH02,TH11,TH12,TH13,TH14,TH15,TH21,TH22,TH23,
     *TH24,TH25,TH26,TH27,TH28,TH29,TH20
      COMMON /AXAX/AX01,AX02,AX03,AX04,AX05,AX06,AX07,AX08,AX09
      COMMON /BMBM/BM01,BM02,BM03,BM04,BM05,BM06
      COMMON /ZMZM/ZM01,ZM02,ZM03,ZM04,ZM05,ZM06,ZM07,ZM08,ZM09,
     1ZM10,ZM11,ZM12,ZM13,ZM14,ZM15,ZM16,ZM17
      COMMON /QMQM/QM01,QM02,QM03,QM04,QM05,QM06,QM07,QM08,QM09,
     1QM10,QM11,QM12,QM13,QM14,QM15,QM16,QM17
      COMMON /TMTM/TM01,TM02,TM03,TM04,TM05,TM06,TM07,TM08,TM09,
     1TM10,TM11,TM12,TM13,TM14,TM15
      DIMENSION XIS(40),YIS(40),THIS(40),THISD(40),MIS(40),
     *XS(40),YS(40),THS(40),THSD(40),MS(40),XA(200),MA(200),MAX(200)
      DIMENSION SPXPA(200),SPFPA(200),DELPA1(200),NIKL(150,50)
      CPCV(T)=1.+0.4/(1.+0.4*(TH/T)**2*EXP(TH/T)/(EXP(TH/T)-1.)**2)
      AA(OM)=((FGAM1+OM**2)/FGAM2)**(FGAM2/2.)/OM
      FTFF(TFF)=39.413E-6*(TFF/873.15)**1.53*983.15/(110.+TFF)
      RLIN(XI)=FKTX1+(FKTX2-FKTX1)*(XI-XFKT1)/(XFKT2-XFKT1)
      FRC(GAM)=SQRT((1.+ZETA)*RC/(GAM+1.))
      FKAP(A0,B0,C0,D0)=(A0*GAM**2.+B0*GAM+C0)/D0
      OPEN(5,FILE='BOCALE',STATUS='OLD')
      OPEN(6,FILE='BOCALS',STATUS='UNKNOWN')
      OPEN(7,FILE='OUT1',STATUS='UNKNOWN')
      OPEN(8,FILE='OUT3',STATUS='UNKNOWN')
      OPEN(9,FILE='BOCAL.DAT',STATUS='UNKNOWN')
      OPEN(10,FILE='OUT2',STATUS='UNKNOWN')
      DATA PI/3.14159265/,GR/9.80665/,MREF/1./,PAM/1.015E5/,TAM/290./,
     *TH/3055.555/,HKM/.9/,NPM/600/,MPM/6./,IC/100/,DYC/.1/,DXC/2./,
     *DEG/57.2957795/,SIGMA/1.0/,ETA/1.8/
 1099 FORMAT(1H0)
C*****
C*****EINLESEN VON DATEN
C*****
      READ(5,*)IDELTA,IY,LC,IPLOT,IPUNCH,IDUESE,ITRAN,ISAUER
      
      write(*,*) IDELTA,IY,LC,IPLOT,IPUNCH,IDUESE,ITRAN,ISAUER
      
 1101 FORMAT(5I5)
      IF(LC.LE.0) LC=1
      IF(LC.GE.IY) LC=IY-1
      READ(5,*)RGAS,TTOT,GAMMA,PTOTB,MD,RC,MEX
      
      write(*,*) RGAS,TTOT,GAMMA,PTOTB,MD,RC,MEX
 1102 FORMAT(7F10.4)
      READ(5,*,iostat=i) CHAR

      READ(5,*) Y0M,YMINM,YEM,DEL2M,DEL3M,BDM,PAMB
      
      write(*,*) Y0M,YMINM,YEM,DEL2M,DEL3M,BDM,PAMB
      pause

!     READ(5,*)CELTH,NPTO,NISG,XSEGG,YSEGG,XSEGD,YSEGD,ALF1D
C     READ(5,*)(XSG(I),YSG(I),I=1,ISG)
      XYMD=MD
C
      IYM=0
      IF(YMINM.EQ.0) IYM=1
CB
C
      IF(IDELTA.EQ.1)GO TO 11
      WRITE(6,1011)
 1011 FORMAT(1H1,'ENTWURF EINER ZWEIDIMENSIONALEN UEBERSCHALLWINDKANAL
     *DUESE'///)
      GO TO 12
   11 WRITE(6,1012)
 1012 FORMAT(1H1,'ENTWURF EINER ROTATIOSSYMMETRISCHEN UEBERSCHALLWIND
     *KANALDUESE'///)
C****
C**** STROEMUNGSMECHANISCHE DATEN
C****
   12 GAM=GAMMA
      FGAM0=GAM/(GAM-1)
      FGAM1=2./(GAM-1)
      FGAM2=(GAM+1)/(GAM-1)
      FGAM3=SQRT(FGAM2)
      FGAM4=SQRT(FGAM1)
C**** DUMMY LINE
C**** DUMMY LINE
      IF (ITRAN.EQ.0) MISO=1.+FGAM2/FGAM1*(0.25/RC-FKAP(0.,1.,57.,
     1576.)/(RC**2)+FKAP(554.,915.,1449.,27648.)/(RC**3))
      IF (ITRAN.EQ.1) MISO=1.+FGAM2/FGAM1*(0.25/(RC+1.)-FKAP(0.,1.,
     1-87.,576.)/((RC+1.)**2)+FKAP(554.,819.,2889.,27648.)/((RC+1.)**3))
      IF (ISAUER.EQ.1) MISO=1.+FGAM2/FGAM1*.25/RC
      GAMTT=CPCV(TTOT)
      TD=TTOT*FGAM1/(FGAM1+MD**2)
      PTOT=1.E5*PTOTB
      PAMB=PAMB*1.E5
      PDB=PTOTB*(TD/TTOT)**FGAM0
      TREF=TTOT*FGAM1/(FGAM1+MREF**2)
      PREF=1.E5*PTOTB*(TREF/TTOT)**FGAM0
      MYREF=FTFF(TREF)
      NYREF=MYREF*RGAS*TREF/PREF
      AREF=SQRT(GAM*RGAS*TREF)
      TAUE=AA(MD)
      IF(IDELTA.EQ.1)GO TO 3
      IF(IYM.EQ.0) YEM=YMINM*TAUE
      IF(IYM.EQ.1) YMINM=YEM/TAUE
      AS=2.E-6*BDM*Y0M
      AE=2.E-6*BDM*YEM
      GO TO 4
    3 IF(IYM.EQ.0) YEM=YMINM*SQRT(TAUE)
      IF(IYM.EQ.1) YMINM=YEM/SQRT(TAUE)
      AS=1.E-6*PI*Y0M**2
      AE=1.E-6*PI*YEM**2
    4 TAUS=AS/AE*TAUE
      MSEST=0.5
      DMS=0.5
    5 DMS=DMS/2.
      TAUEST=AA(MSEST)
      DFM=(TAUEST-TAUS)/TAUS
      IF(ABS(DFM).LT.1.E-3) GO TO 6
      IF(DFM.LT.0.)GO TO 9
      MSEST=MSEST+DMS
      GO TO 5
    9 MSEST=MSEST-DMS
      GO TO 5
    6 PD=1.E5*PDB
      RHOTOT=PTOT/RGAS/TTOT
      TSTT=FGAM1/(FGAM1+MSEST**2)
      WS=MSEST*SQRT(GAM*RGAS*TSTT*TTOT)
      WD=MD*SQRT(GAM*RGAS*TD)
      RHOS=RHOTOT*TSTT**(1./(GAM-1.))                                   0
      RHOAM=PAM/RGAS/TAM
      RHOD=PD/RGAS/TD
      DMDT=AE*WD*RHOD
      DMDTV=DMDT/RHOAM
      PS=RHOS*RGAS*TSTT*TTOT
      FORCE=PS*AS-PAM*(AS-AE)-PD*AE-DMDT*(WD-WS)
      FORT=FORCE/GR/1000.
      stop
C****
C**** EINGABEDATEN AUSDRUCKEN
C****
      WRITE(6,1069)
 1069 FORMAT(1H0,'VORGEGEBENE DATEN'//)
      WRITE(6,1002)SIGMA
 1002 FORMAT(1HO,'FORMFAKTOR DER HALL-KONTUR            SIGMA =',F12.4)
      WRITE(6,1001)RC
 1001 FORMAT(1H0,'KRUEMMUNGSRADIUS DER HALL-KONTUR         RC =',F12.4)
      WRITE(6,1003)MISO
 1003 FORMAT(1H0,'MACHZAHL DER ISOTACHE                  MISO =',F12.4)
      WRITE(6,1004)IY
 1004 FORMAT(1H0,'ANZAHL DER ISOTACHENPUNKTE               IY =',I12)
      WRITE(6,1016)LC
 1016 FORMAT(1H0,'SCHALLNAHE KONTUR BEI                    LC =',I12)
      WRITE(6,1005)RGAS
 1005 FORMAT(1H0,'GASKONSTANTE                           RGAS =',F12.4,
     *' J/KG KELVIN')
      WRITE(6,1044)GAMMA
 1044 FORMAT(1H0,'CP/CV                                 GAMMA =',F12.4)
      WRITE(6,1006)TTOT
 1006 FORMAT(1H0,'RUHETEMPERATUR                         TTOT =',F12.4,
     *' KELVIN     ')
      WRITE(6,1007)PTOTB
 1007 FORMAT(1H0,'RUHEDRUCK                             PTOTB =',F12.4,
     *' BAR         ')
      WRITE(6,1008)MD
 1008 FORMAT(1H0,'ENTWURFSMACHZAHL                         MD =',F12.4)
      WRITE(6,1009)DEL2M
 1009 FORMAT(1H0,'IMPULSVERLUSTDICKE BEI M = 1          DEL2M =',F12.4,
     *' MM         ')
      WRITE(6,1010)DEL3M
 1010 FORMAT(1H0,'ENERGIEVERLUSTDICKE BEI M = 1         DEL3M =',F12.4,
     *' MM         ')
      WRITE(6,1091)Y0M
 1091 FORMAT(1H0,'HALBE DUESENHOEHE AM EINTRITT           Y0M =',F12.4,
     *' MM         ')
      IF(IYM.EQ.0) WRITE(6,1092)YMINM
 1092 FORMAT(1H0,'KLEINSTE HALBE DUESENHOEHE            YMINM =',F12.4,
     *' MM         ')
      IF(IYM.EQ.1) WRITE(6,1014)YEM
 1014 FORMAT(1H0,'HALBE DUESENHOEHE AM AUSTRITT           YEM =',F12.4,
     *' MM         ')
      IF(IDELTA.EQ.1) GO TO 2
      WRITE(6,1013)BDM
 1013 FORMAT(1H0,'DUESENBREITE                            BDM =',F12.4,
     *' MM         ')
    2 IF(IDUESE.EQ.1)GO TO 1051
 2002 WRITE(6,1068)
 1068 FORMAT(1H0//' STROEMUNGSMECHANISCHE DATEN'//)
      WRITE(6,1045)GAMTT
 1045 FORMAT(1H0,'CP/CV DER RUHETEMPERATUR              GAMTT =',F12.4)
      WRITE(6,1093)TD
 1093 FORMAT(1H0,'ENDTEMPERATUR                            TD =',F12.4,
     *' KELVIN     ')
      WRITE(6,1094)PDB
 1094 FORMAT(1H0,'ENDDRUCK                                PDB =',F12.4,
     *' BAR       ')
      WRITE(6,1057)TAUE,GAM
 1057 FORMAT(1H0,'FLAECHENVERHAELTNIS                    TAUE =',F12.4,
     *' FUER CP/CV =',F8.4)
      WRITE(6,1050)DMDT,DMDTV
 1050 FORMAT(1H0,'DURCHSATZ                              DMDT =',E12.4,
     *' KG/S        =',E12.4,' MM**3/S   ')
      WRITE(6,1052)FORCE,FORT
 1052 FORMAT(1H0,'AM DUESENFLANSCH WIRKENDE ZUGKRAFT    FORCE =',E12.4,
     *' NEWTON      =',E12.4,' TO        '//)
      IF(IDUESE.EQ.1)GO TO 2397
 1051 HIY=IY-1
C****
C**** SELECIONA A SERIE DO TRANSONICO
C****
      IF (ITRAN.EQ.0) GO TO 1212
      RC=RC+1.
      CNVR=(1.-.5/RC-1./8./(RC**2.)-1./16./(RC)**3.)
      WRITE(6,1211)
 1211 FORMAT(1H0,'SERIE DE KLIEGEL')
      GO TO 1214
 1212 WRITE(6,1213)
      CNVR=1.
 1213 FORMAT(1H0,'SERIE DE HALL')
      IF (ISAUER.EQ.1) WRITE(6,1215)
1215  FORMAT(1H0,'SERIE DE SAUER')
 1214 CONTINUE
C****
C**** SCHNITTPUNKTE ISOTACHEN/HALL-KONTUR
C****
      RC2=RC**2
      RC3=RC**3
      QQ=0.
      IF(IDELTA.EQ.1) GO TO 25
      ZETA=0.
      ZX=FRC(GAM)
      CALL FGAMM0(GAM)
      GO TO 29
   25 ZETA=1.
      ZX=FRC(GAM)
      CALL FGAMM1(GAM,ITRAN)
      CALL FGAMM3(GAM)
   29 Y=1
   23 Y2=Y**2
      Y4=Y**4
      Y6=Y**6
      IF (ITRAN.EQ.0) CALL HALL FI (QQ,Y2,Y4,Y6,F0,F1,F2)
      IF (ITRAN.EQ.1) CALL HALL FM(QQ,Y2,Y4,Y6,F0,F1,F2)
      Z=(F0+F1/RC+F2/RC2)
      IF (ISAUER.EQ.1) Z=F0
      X=Z/ZX
      X2=X**2
      X4=X**4
      IF (ITRAN.EQ.1) YC0=1.+X2/2./(RC-1.)+SIGMA*X4/8./((RC-1.)**3)     2260
      IF (ITRAN.EQ.0) YC0=1.+X2/2./RC+SIGMA*X4/8./RC3                   70
      IF (ABS((YC0-Y)/Y).LT.1.E-4) GO TO 22
      Y=YC0
      GO TO 23
C****
C**** SCHALLLINIE UND ISOTACHE
C****
   22 DELYS=YC0/HIY
      DYTT=1./HIY
      YTT=-DYTT
      YSS=-DELYS
      DO 20 IZ=1,IY
      YTT=YTT+DYTT
      YSS=YSS+DELYS
      Y=YSS
      Y2=Y**2
      Y3=Y**3
      Y4=Y**4
      Y5=Y**5
      Y6=Y**6
      Y7=Y**7
      IF (ITRAN.EQ.0) CALL HALL FI (QQ,Y2,Y4,Y6,F0,F1,F2)
      IF (ITRAN.EQ.1) CALL HALL FM(QQ,Y2,Y4,Y6,F0,F1,F2)
      ZS=(F0+F1/RC+F2/RC2)*CNVR
      IF (ISAUER.EQ.1) ZS=F0*CNVR
      ZS2=ZS**2
      ZS3=ZS**3
      IF (ITRAN.EQ.0) CALL HALL TI (F0,F1,F2,Y,Y3,Y5,Y7,TH1,TH2,TH3)
      IF (ITRAN.EQ.1) CALL HALL TM (QQ,Y,Y3,Y5,Y7,TH1,TH2,TH3)
      IF (ITRAN.EQ.0) CALL HALL QI (ZS,Y2,Y4,Y6,Q1,Q2,Q3)
      IF (ITRAN.EQ.1) CALL HALL QM(ZS,Y2,Y4,Y6,Q1,Q2,Q3)
      QS2=(1.+Q1/RC+Q2/RC2+Q3/RC3)**2
      IF (ISAUER.EQ.1) QS2=(1.+Q1/RC)**2.
      IF(QQ.NE.0.)GO TO 24
      XS(IZ)=ZS/ZX/CNVR
      YS(IZ)=YSS
      THS(IZ)=(TH1/RC+TH2/RC2+TH3/RC3)/ZX
      IF (ISAUER.EQ.1) THS(IZ)=TH1/RC/ZX
      THSD(IZ)=DEG*THS(IZ)
      MS(IZ)=SQRT(FGAM1*QS2/(FGAM2-QS2))
      GO TO 20
C****
C**** ESCOLHE A PARTIDA DA CARACTERISTICA INICIAL
C****
  24  XIS(IZ)=ZS/ZX/CNVR
      YIS(IZ)=YSS
      THIS(IZ)=(TH1/RC+TH2/RC2+TH3/RC3)/ZX
      IF (ISAUER.EQ.1) THIS(IZ)=TH1/RC/ZX
      THISD(IZ)=DEG*THIS(IZ)
      MIS(IZ)=SQRT(FGAM1*QS2/(FGAM2-QS2))
   20 CONTINUE
      IF(QQ.NE.0.)GO TO 7
      QISO=SQRT(FGAM2*MISO**2/(FGAM1+MISO**2))
      QQ=RC*(QISO-1)
      GO TO 29
    7 WRITE(6,1097)
 1097 FORMAT(1H1,'SCHALLLINIE'/)
      WRITE(6,1098)
 1098 FORMAT(1H ,4X,1HN,6X,1HX,9X,1HY,8X,'THETA',6X,1HM,/,32X,'GRAD'/)
      WRITE(6,1015)(IZ,XS(IZ),YS(IZ),THSD(IZ),MS(IZ),IZ=1,IY)
 1015 FORMAT(1H ,I5,4F10.4)
      IF(IPLOT.GT.0) WRITE(10,1032)IDELTA,IY
 1032 FORMAT(2I10)
      IF(IPLOT.GT.0) WRITE(10,1035)(XS(IZ),YS(IZ),IZ=1,IY)
 1035 FORMAT(6F10.4)
      WRITE(6,1017)
 1017 FORMAT(1H0////1X,'ISOTACHE'/)
      WRITE(6,1098)
C1098 FORMAT(1H ,4X,1HN,9X,1HX,9X,1HY,5X,'THETA',9X,1HM,/32X,'GRAD')
      WRITE(6,1015)(IZ,XIS(IZ),YIS(IZ),THISD(IZ),MIS(IZ),IZ=1,IY)
C
      IF(IPLOT.GT.0) WRITE(10,1035)(XIS(IZ),YIS(IZ),IZ=1,IY)
C
C****
C**** ANFANGSREIHE
C****
      NCHAA=IY+1
      NCHA=NCHAA
      NDIM=NCHAA
      XF(NCHA,1)=XIS(IY)
      YF(NCHA,1)=YIS(IY)
      THF(NCHA,1)=THIS(IY)
      MF(NCHA,1)=MIS(IY)
      NI(NCHA,1)=0
      DO 60 IZ=1,IY
      RR=NCHA+1-IZ
      LL=IY-IZ
      NR=RR-1
      IF(LL-1) 54,56,52
   52 CALL CHARA F(XIS(LL),YIS(LL),THIS(LL),MIS(LL),
     *XF(RR,1),YF(RR,1),THF(RR,1),MF(RR,1),
     *XF(NR,1),YF(NR,1),THF(NR,1),MF(NR,1),NIT,IDELTA,GAM)
   57 NI(NR,1)=20-NIT
      GO TO 60
   56 CALL CHARA A (XIS(1),MIS(1),XF(RR,1),YF(RR,1),THF(RR,1),MF(RR,1),
     *XF(NR,1),YF(NR,1),THF(NR,1),MF(NR,1),NIT,IDELTA,GAM)
      GO TO 57
   60 CONTINUE
   54 CALL CHARA R(XF(2,1),YF(2,1),THF(2,1),MF(2,1),
     *XF(1,1),MF(1,1),NIT,IDELTA,GAM)
      YF(1,1)=0.
      THF(1,1)=0.
      NI(1,1)=20-NIT
      FMA=MF(1,1)
C****
C**** ACHSENPUNKTE NACH HALL
C****
   49 Y2=0.
      Y4=0.
      Y6=0.
      DXA=(XF(1,1)-XS(IY))/HIY
      X=XS(IY)-DXA
      DO 50 NA=1,IY
      X=X+DXA
      XA(NA)=X
      ZA=ZX*XA(NA)*CNVR
      IF (ITRAN.EQ.0) CALL HALL QI (ZA,Y2,Y4,Y6,Q1,Q2,Q3)
      IF (ITRAN.EQ.1) CALL HALL QM(ZA,Y2,Y4,Y6,Q1,Q2,Q3)
      QA=1.+Q1/RC+Q2/RC2+Q3/RC3
      IF (ISAUER.EQ.1) QA=1.+Q1/RC
      CALL AX(ZX*CNVR,ZA,QA,D202,D203,D302,D303,D304,RC,RC2,RC3,FGAM0,
     *FGAM1,FGAM2,MA(NA),MAX(NA),ISAUER)
   50 CONTINUE
      NA=IY
      NAA=NA
      DELMA=100.*(FMA-MA(NAA))/FMA
C****
C**** ANFANGSREIHE AUSDRUCKEN
C****
      WRITE(6,1020)
 1020 FORMAT(1H0////1X,'ANFANGSREIHE'/)
      WRITE(6,1021)
 1021 FORMAT(1H ,3X,2HPN,9X,1HX,9X,1HY,5X,5HTHETA,9X,1HM,2X,
     *3HNIT,1X,4HNCHA,8X,2HMA,7X,3HMAX,5X,5HDELMA/32X,4HGRAD,
     *47X,3H0/0/)
      PN0=-97
      DO 40 L=1,NCHAA
      PN0=PN0+99
      PNF(L,1)=PN0
      THFD=DEG*THF(L,1)
      WRITE(6,1040)PNF(L,1),XF(L,1),YF(L,1),THFD,MF(L,1),NI(L,1)
 1040 FORMAT(1H ,I5,4F10.4,I5)
      IF(L.EQ.1) WRITE(6,1041)NCHA,MA(NAA),MAX(NAA),DELMA
 1041 FORMAT(1H+,50X,I5,3F10.4)
   40 CONTINUE
      IF(IPLOT.GT.0) WRITE(10,1033)NCHA
 1033 FORMAT(I5)
      IF(IPLOT.GT.0) WRITE(10,1031)(PNF(L,1),XF(L,1),YF(L,1),L=1,NCHA)
 1031 FORMAT(3(I5,2F9.4))
C****
C**** KONTURPUNKTE IN SCHALLNAEHE
C****
      LA=LC+1
      XC(1)=XF(LA,1)
      XN=XC(1)
      YC(1)=YF(LA,1)
      YCX(1)=TAN(THF(LA,1))
      JS=1
   71 JS=JS+1
      IF(YS(JS).LT.YF(LA,1))GO TO 71
      XS1=XS(JS)+(XS(JS-1)-XS(JS))*(YF(LA,1)-YS(JS))/(YS(JS-1)-YS(JS))
      DELX=(XS1-XF(LA,1))/IC
      DO 240 J=1,IC
      IT=-1
      XN=XN+DELX
      XC(J+1)=XN
      Z=ZX*XN*CNVR
      TAN1=YCX(J)
   16 DELY=DELX*TAN1
      YC(J+1)=YC(J)+DELY
      Y1=YC(J+1)
      Y3=Y1**3
      Y5=Y1**5
      Y7=Y1**7
      CALL HALL TH(Z,Y1,Y3,Y5,Y7,TH1,TH2,TH3)
      THH=(TH1/RC+TH2/RC2+TH3/RC3)/ZX
      IF (ISAUER.EQ.1) THH=TH1/RC/ZX
      TANH=TAN(THH)
      TAN2=(YCX(J)+TANH)/2.
      IT=IT+1
      IF(ABS(TAN1-TAN2).LT.1.E-8.OR.IT.EQ.50)GO TO 18
      TAN1=TAN2
      GO TO 16
   18 YCX(J+1)=TAN1
  240 CONTINUE
      XFKT1=XS1
      XFKT2=XF(LA,1)
      FKTX1=1.
      FKTX2=MF(LA,1)
      NC=0
      ICC=IC+1
      DO 250 J=1,ICC,10
      L=ICC+1-J
      NC=NC+1
      IF(J.EQ.1)GO TO 44
      XC(NC)=XC(L)
      YC(NC)=YC(L)
      YCX(NC)=YCX(L)
   44 MC(NC)=RLIN(XC(L))
  250 CALL TECEPE(FGAM0,FGAM1,TREF,MREF,PREF,NC,GAM)
      XC(1)=XC(ICC)
      YC(1)=YC(ICC)
      YCX(1)=YCX(ICC)
      NCA=NC
      IF(IDUESE.NE.1)GO TO 1036
      IKL=0
      ISG=NC
      DO 251 J=1,NC
      XSG(J)=XC(J)
      YSG(J)=ATAN(YCX(J))*57.29
  251 CONTINUE
      CALL SPLNE1
      CALL RAIZ(SPLINT,XC(1),XC(NC),XGRT,1.E-4,IER)
      DO 252 J=1,NC
      YSG(J)=YC(J)
  252 CONTINUE
      CALL SPLNE1
      YGRT=SPLINT(XGRT)
      YBZG=YMINM/YGRT
      IGRT=IGT(NC,XGRT,XC)
      NC=NC+1
      DO 253 J=IGRT,NC
      K=NC+IGRT-J
      XC(K)=XC(K-1)
      YC(K)=YC(K-1)
      YCX(K)=YCX(K-1)
      MC(K)=MC(K-1)
  253 CONTINUE
      XC(IGRT)=XGRT
      YC(IGRT)=YGRT
      YCX(IGRT)=0.
      MC(IGRT)=RLIN(XC(IGRT))
      CALL TECEPE(FGAM0,FGAM1,TREF,MREF,PREF,NC,GAM)
      NCA=NC
C****
C**** STROEMUNGSFELD DES EXPANSIONSFAECHERS
C****
 1036 WRITE(6,1037)
 1037 FORMAT(1H0,'STROEMUNGSFELD DES EXPANSIONSFAECHERS'/)
      WRITE(6,1084)
 1084 FORMAT(1H0,'   PN',6X,1HX,9X,1HY,8X,'THETA',6X,1HM,4X,
     *'  NIT',5H NCHA,16X,2HPN,7X,1HX,9X,1HY,8X,'THETA',6X,1HM,
     *'       NIT'/32X,4HGRAD,66X,4HGRAD)
      MF(NCHA,2)=MF(NCHA,1)
      CALL PRAMEY(FGAM2,FGAM3,NCHA,NUEA)
      ALPHA=THF(NCHA,1)-NUEA
      PN0=PNF(NCHA,1)
      XF0=XF(NCHA,1)
      YF0=YF(NCHA,1)
      M0=MF(NCHA,1)
      HNPM=NPM
      HKTR=HKM*PI/2./(MPM-M0)
      FMTR=TAN(HKM*PI/2.)/HNPM
      ITER=1
      DO 70 NY=1,1000
      NNY=0
   92 HNY=FLOAT(NY)
      NNY=NNY+1.
      IF(NY.EQ.1.AND.NNY.EQ.1) HNY=HNY/3.
      IF(NY.EQ.1.AND.NNY.EQ.2) HNY=2.*HNY/3.
      IF(NY.EQ.2.AND.NNY.EQ.1) HNY=3.*HNY/4.
      NCHA=NCHA+1
      IF(NCHA.GT.150)GO TO 111
      RR=NCHA
      NA=NA+1
      IF(NA.GT.150)GO TO 112
      NC=NC+1
      IF(NC.GT.170)GO TO 113
      PN0=PN0+101
      PNF(RR,2)=PN0
      XF(RR,2)=XF0
      YF(RR,2)=YF0
   98 TMTR=HNY*FMTR
      MF(RR,2)=M0+ATAN(TMTR)/HKTR
      CALL PRAMEY(FGAM2,FGAM3,NCHA,NUE)
      THF(RR,2)=ALPHA+NUE
      NI(RR,2)=0
      RR=RR+1
   72 RR=RR-1
      LL=RR-2
      NR=RR-1
      IF(LL-1) 74,75,76
   76 CALL CHARA F(XF(LL,1),YF(LL,1),THF(LL,1),MF(LL,1),
     *XF(RR,2),YF(RR,2),THF(RR,2),MF(RR,2),
     *XF(NR,2),YF(NR,2),THF(NR,2),MF(NR,2),NIT,IDELTA,GAM)
   73 PNF(NR,2)=PNF(RR,2)-99
      NI(NR,2)=20-NIT
      GO TO 72
   75 CALL CHARA A(XF(LL,1),MF(LL,1),
     *XF(RR,2),YF(RR,2),THF(RR,2),MF(RR,2),
     *XF(NR,2),YF(NR,2),THF(NR,2),MF(NR,2),NIT,IDELTA,GAM)
      GO TO 73
   74 CALL CHARA R(XF(2,2),YF(2,2),THF(2,2),MF(2,2),
     *XF(1,2),MF(1,2),NIT,IDELTA,GAM)
      PNF(1,2)=PNF(2,2)-99
      YF(1,2)=0.
      THF(1,2)=0.
      NI(1,2)=20-NIT
      DELM=MD-MF(1,2)
      IF(ABS(DELM).LT.1.E-3)GO TO 96
      IF(ITER.LT.0)GO TO 95
      IF(DELM.GT.0.AND.DELM.LE.0.05)GO TO 51
      IF(DELM.LT.0)GO TO 78
   96 IF(IDUESE.NE.1)GO TO 196
C     IF(MF(1,2).LE.MEX)GO TO 196
      IKL=IKL+1
      IF(IKL.GT.1)GO TO 194
      NDIM0=NDIM
      NA0=NA
  194 NKL(IKL)=NCHA
      DO 195 I=1,NCHA
      IW=NCHA+1-I
      PNKL(IW,IKL)=PNF(I,2)
      XKL(IW,IKL)=XF(I,2)
      YKL(IW,IKL)=YF(I,2)
      MKL(IW,IKL)=MF(I,2)
      THKL(IW,IKL)=THF(I,2)
      NIKL(IW,IKL)=NI(I,2)
  195 CONTINUE
      DO 495 I=1,NCHA
      PNF(I,1)=PNF(I,2)
      XF(I,1)=XF(I,2)
      YF(I,1)=YF(I,2)
      THF(I,1)=THF(I,2)
      MF(I,1)=MF(I,2)
      NI(I,1)=NI(I,2)
  495 CONTINUE
      GO TO 197
  196 CALL RESULT(NCHA,IPLOT)
      NDIM=NDIM+NCHA
      XA(NA)=XF(1,2)
      MA(NA)=MF(1,2)
  197 NSA=NC-1
      NSB=NC
      NSJ=NCHA+1
   77 NSJ=NSJ-1
      NSK=NSJ-1
      IF(XF(NSK,2).GT.XC(NSA))GO TO 46
      GO TO 77
   46 CALL STREAM(IST,XC(NSA),YC(NSA),YCX(NSA),XF(NSJ,2),YF(NSJ,2),THF
     *(NSJ,2),XF(NSK,2),YF(NSK,2),THF(NSK,2),XC(NSB),YC(NSB),YCX(NSB))
      IF(IST.LT.0)GO TO 77
      IF(IDUESE.NE.1)GO TO 151
      IF(IKL.EQ.0)GO TO 151
      XB(IKL)=XC(NSB)
      YB(IKL)=YC(NSB)
  151 CALL LAGRAN(XC(NSB),MC(NSB),NSJ)
      CALL TECEPE(FGAM0,FGAM1,TREF,MREF,PREF,NC,GAM)
      IF(ITER.LT.0)GO TO 99
      IF(NY.EQ.1.AND.NNY.LE.2.OR.NY.EQ.2.AND.NNY.EQ.1)GO TO 92
      GO TO 70
   51 DHNY=2.
      GO TO 95
   78 DHNY=1.
   95 ITER=-1
      RR=NCHA
      DHNY=DHNY/2.
      IF(DELM.GT.0.) HNY=HNY+DHNY
      IF(DELM.LT.0.) HNY=HNY-DHNY
      GO TO 98
   70 CONTINUE
   99 IF(IDUESE.NE.1)GO TO 199
      CALL PUNKT D(MEX,PAMB,PTOT,TTOT,RGAS,GAM,YBZG,FORCE,DMDT,XD,YD,YE,
     .XBP,YBP,IKN)
      FORT=FORCE/GR/1000.
      DMDTV=DMDT/RHOAM
      YEM=YBZG*YE
      TAUE=(YEM/YMINM)**2
      DO 397 I=1,IMTY
      YXE(I)=YXE(I)*YE
  397 CONTINUE
      I=IGT(NC,XBP,XC)
      YCX(I)=AMED(YCX(I),YCX(I-1),XC(I),XC(I-1),XBP)
      NC=I
      XC(NC)=XBP
      YC(NC)=YBP
      NA=NA0-1
      NDIM=NDIM0
      DO 798 I=1,IKN
      NA=NA+1
      XA(NA)=XKL(NKL(I),I)
      MA(NA)=MKL(NKL(I),I)
      NDIM=NDIM+NKL(I)
      KL=1
      WRITE(6,2111)
 2111 FORMAT(1H )
      DO 198 IW=1,NKL(I)
      KL=-KL
      J=NKL(I)+1-IW
      IF(KL.EQ.1)GO TO 299
      WRITE(6,2101)PNKL(J,I),XKL(J,I),YKL(J,I),THKL(J,I)*57.2957795,
     .MKL(J,I),NIKL(J,I)
 2101 FORMAT(1H ,I5,4F10.4,I5)
      GO TO 298
  299 WRITE(6,2106)PNKL(J,I),XKL(J,I),YKL(J,I),THKL(J,I)*57.2957795,
     .MKL(J,I),NIKL(J,I)
 2106 FORMAT(1H+,70X,I5,4F10.4,I5)
  298 IF(IW.EQ.1)WRITE(6,2107)NKL(I)
 2107 FORMAT(1H+,50X,I5)
 198  CONTINUE
      IF(IPLOT.NE.1)GO TO 798
      WRITE(10,3114)NKL(I)
 3114 FORMAT(I5)
      WRITE(10,2105)(PNKL(J,I),XKL(J,I),YKL(J,I),J=NKL(I),1,-1)
 2105 FORMAT(3(I5,2F9.4))
  798 CONTINUE
      NA=NA-1
C****
C**** ABSCHLIESSENDES STROEMUNGSFELD
C****
  199 NAD=NA
      NCD=NC
      XD=XF(1,1)
      WRITE(6,1039)
 1039 FORMAT(1H0,'ABSCHLIESSENDES STROEMUNGSFELD'/)
      WRITE(6,1084)
C
C
C
      ISE=1
      NSF=2
      DELY=YF(3,1)-YF(2,1)
      IF(DELY.LT.0.)DELY=-DELY
      DELX=XF(1,1)-XF(2,1)
      IF(DELY.LT.DELX)DELY=DELX
      FMD=SQRT(MF(1,1)**2-1.)
      YM=0.
      IF(IDUESE.NE.1)GO TO 81
      YM=YD
      K=IGT(IMTY,YD,YXE)
      TH0=AMED(THY(K),THY(K-1),YXE(K),YXE(K-1),YD)
      AM0=AMED(AMY(K),AMY(K-1),YXE(K),YXE(K-1),YD)
      AUX=1./AM0
      PHY0=TH0+ASIN(AUX)
   81 PN0=PNF(1,1)
      RR=1
      LL=2
      NR=2
      NC=NC+1
      IF(NC.GT.170)GO TO 113
      NSA=NC-1
      NSB=NC
      NCHA=1
      YM=YM+DELY
      PN0=PN0+101
      PNF(RR,2)=PN0
      YF(RR,2)=YM
      NI(RR,2)=0
      IF(IDUESE.EQ.1)GO TO 187
      XF(RR,2)=XD+YM*FMD
      THF(RR,2)=0.
      MF(RR,2)=MF(1,1)
      GO TO 87
  187 I=IGT(IMTY,YM,YXE)
      MF(RR,2)=AMED(AMY(I),AMY(I-1),YXE(I),YXE(I-1),YM)
      THF(RR,2)=AMED(THY(I),THY(I-1),YXE(I),YXE(I-1),YM)
      AUX=1./MF(RR,2)
      PHY=THF(RR,2)+ASIN(AUX)
      XF(RR,2)=XF(1,1)+DELY/TAN(.5*(PHY+PHY0))
      PHY0=PHY
   87 CALL CHARA F(XF(LL,1),YF(LL,1),THF(LL,1),MF(LL,1),
     *XF(RR,2),YF(RR,2),THF(RR,2),MF(RR,2),
     *XF(NR,2),YF(NR,2),THF(NR,2),MF(NR,2),NIT,IDELTA,GAM)
      PNF(NR,2)=PNF(RR,2)+99
      NI(NR,2)=20-NIT
      NCHA=NCHA+1
      IF(NCHA.GT.150)GO TO 111
      IF(XF(NR,2).GT.XC(NSA))GO TO 89
   85 NSJ=NR
      NSK=NSJ-1
      IF(NSK.EQ.0)GO TO 86
   88 CALL STREAM(IST,XC(NSA),YC(NSA),YCX(NSA),XF(NSJ,2),YF(NSJ,2),
     *THF(NSJ,2),XF(NSK,NSF),YF(NSK,NSF),THF(NSK,NSF),XC(NSB),
     *YC(NSB),YCX(NSB))
      IF(ISE.LT.0)GO TO 91
      IF(IST.LT.0)GO TO 79
   91 CALL LAGRAN(XC(NSB),MC(NSB),NSJ)
      IF(ISE.LT.0) MC(NSB)=MF(1,1)
      CALL TECEPE(FGAM0,FGAM1,TREF,MREF,PREF,NC,GAM)
      CALL RESULT(NCHA,IPLOT)
      NDIM=NDIM+NCHA
      IF(ISE.LT.0)GO TO 93
      GO TO 81
   86 ISE=-1
      NSJ=1
      NSK=1
      NSF=1
      NCHA=1
      GO TO 88
   89 RR=RR+1
      LL=LL+1
      NR=NR+1
      GO TO 87
   79 RR=RR-1
      LL=LL-1
      NR=NR-1
      GO TO 85
  111 WRITE(6,1081)NCHA
 1081 FORMAT(1H0,'COMMON/A/-DIMENSIONIERUNG UEBERSCHRITTEN NCHA',I8)
      STOP
  112 WRITE(6,1082)NA
 1082 FORMAT(1H0,'XA-DIMENSIONIERUNG UEBERSCHRITTEN   NA =',I8)
      STOP
  113 WRITE(6,1083)NC
 1083 FORMAT(1H0,'COMMON/D/-DIMENSIONIERUNG UEBERSCHRITTEN    NC =',I8)
      STOP
   93 NCE=NC
      WRITE(6,1090) NDIM
 1090 FORMAT(1H0/' ANZAHL DER CHARAKTERISTIKENPUNKTE    NDIM =',I10//)
C****
C**** MACHZAHLGRADIENT AUF DER ACHSE
C****
      IF(IDUESE.EQ.1)GO TO 2002
 2397 WRITE(6,1054)
 1054 FORMAT(1H0/' MACHZAHLVERLAUF AUF DER ACHSE'/)
      WRITE(6,1029)
 1029 FORMAT(1H0,4(2X,1HN,7X,2HXA,7X,2HMA,7X,3HMAX,1X)/)
      WRITE(6,1030)(N,XA(N),MA(N),MAX(N),N=1,NAA)
 1030 FORMAT(4(1H ,I4,3F9.4))
      IF(IPLOT.GT.0) WRITE(10,1023)NAA,NAD,NDIM
 1023 FORMAT(3I10)
      IF(IPLOT.GT.0) WRITE(10,1025)(XA(N),MA(N),MAX(N),N=1,NAA)
 1025 FORMAT(6F10.4)
      WRITE(6,1099)
C1099 FORMAT(1H0)
      MA(NAA)=FMA
      N1=1
      N2=NAD
      DO 30 J=N1,N2
      SPX(J)=XA(J)
   30 SPF(J)=MA(J)
      CALL SPLINE(N1,N2)
      CALL TAYLOR(N2,B(N2))
      WRITE(6,1030)(N,XA(N),MA(N),B(N),N=NAA,N2)
C1030 FORMAT(4(1H ,I5,3F10.4))
      IF(IPLOT.GT.0) WRITE(10,1025)(XA(N),MA(N),B(N),N=NAA,N2)
C1025 FORMAT(6F10.4)
C****
C**** MACHZAHLGRADIENT AN DER KONTUR
C****
      TAUEFF=YC(NCE)/YC(1)
      IF(IDUESE.EQ.1) TAUEFF=YC(NCE)/YC(IGRT)
      IF(IDELTA.EQ.1) TAUEFF=TAUEFF**2
      WRITE(6,1058)
 1058 FORMAT(1H0/' DATEN DER POTENTIALKONTUR (OHNE UNTERSCHALLKONTUR)')
      WRITE(6,1042)TAUEFF
 1042 FORMAT(1H0,'EFFEKTIVES FLAECHENVERHAELTNIS   TAUEFF =',F8.4/)
      WRITE(6,1028)
 1028 FORMAT(1H0,4X,1HN,6X,2HXC,8X,2HYC,7X,3HTHC,8X,2HMC,7X,
     *3HMCX,8X,2HTC,           7X,3HCPC,8X,2HPC/31X,4HGRAD,24X,6HKELVIN,
     *16X,3HBAR/)
      IF(IPLOT.GT.0) WRITE(10,1046)NCA,NCD,NCE
 1046 FORMAT(3I10)
      N1=1
      N2=NCA
      DO 120 L=1,3
      IF(L.EQ.2) N2=NCD
      IF(L.EQ.3) N2=NCE
      DO 110 N=N1,N2
      SPX(N)=XC(N)
  110 SPF(N)=MC(N)
      CALL SPLINE(N1,N2)
      CALL TAYLOR(N2,B(N2))
      DO 90 K=N1,N2
      IF(IPLOT.GT.0) WRITE(10,1026)XC(K),YC(K),YCX(K),MC(K),B(K)
 1026 FORMAT(5F10.4)
      THC=DEG*ATAN(YCX(K))
   90 WRITE(6,1027)K,XC(K),YC(K),THC,MC(K),B(K),TC(K),CPC(K),PC(K)
 1027 FORMAT(1H ,I5,8F10.4)
      N1=N2
  120 WRITE(6,1099)
C1099 FORMAT(1H )
C****
C**** KRUEMMUNGSVERLAUF DER POTENTIALKONTUR
C****
      WRITE(6,1018)
 1018 FORMAT(1H1,'KRUEMMUNGSVERLAUF DER POTENTIALKONTUR (OHNE UNTERSCHAL
     *KONTUR)'//4H   N,6X,1HX,9X,1HY,'      KRUEMMUNG',25X,'GRAPHISCHE D
     *ARSTELLUNG DES KRUEMMUNGSVERLAUFS',/,10X,2HRE,8X,2HRE,8X,4H1/RE/)
      N1=1
      N2=NCA
      DO 200 L=1,3
      IF(L.EQ.2) N2=NCD
      IF(L.EQ.3) N2=NCE
      DO 190 J=N1,N2
      SPX(J)=XC(J)
  190 SPF(J)=YC(J)
      CALL SPLINE(N1,N2)
      CALL CURVAT(N1,N2)
  200 N1=N2
C****
C**** GRENZSCHICHT MIT VERDRAENGUNGSDICKE
C****
      ELREF=XC(NCE)-XC(1)
      R0=10**6
      DO 140 J=1,NCE
      SPX(J)=XC(J)/ELREF
      IF(IDELTA.EQ.0) SPF(J)=R0
      IF(IDELTA.EQ.1) SPF(J)=-YC(J)/ELREF
  140 CONTINUE
      EM=YEM/YC(NCE)
      ELREFM=EM*ELREF
      EML=ELREFM
      REREF=MREF*AREF*ELREFM/NYREF/1000
      PREFB=1.E-5*PREF
      DEL2=DEL2M/EML
      DEL3=DEL3M/EML
      WRITE(6,1064)MREF,TREF,PREFB,ELREFM,REREF,AREF,NYREF
 1064 FORMAT(1H0//' TURBULENTE GRENZSCHICHT',14X,'BEZUGSGROESSEN'//
     *' MACHZAHL',31X,'MREF =',F10.4//' TEMPERATUR',29X,'TREF =',F10.4,
     *'  KELVIN'//' DRUCK',34X,'PREF =',F10.4,'  BAR'//' LAENGE',31X,
     *'ELREFM =',F10.4,'   MM'//' REYNOLDSZAHL',26X,'REREF =',E10.4//
     *' SCHALLGESCHWINDIGKEIT',18X,'AREF =',F10.4,'   M/S'//
     *' KINEMATISCHE ZAEHIGKEIT',15X,'NYREF =',E10.4,' M**2/S'//)
      IF(IDUESE.NE.0)GO TO 994
      CALL ENSAIO(NCE,CELTH,NPTO,ELREFM,IDELTA)
      CALL DIFSUP(YEM,YSEGG,XSEGG,NISG,XYMD,GAM,TREF,MREF,PREF,NCE,ELREF
     #M,NDIF,XSEGD,YSEGD,ALF1D)
      NCE=NCE+NDIF
      DO 993 I=1,NCE
      WRITE(7,*)I,SPX(I),SPF(I),PC(I)
 993  CONTINUE
 994  CALL B LAYER(DEL2,DEL3,NCE,MREF,TREF,REREF,GAM)
      IF(IDUESE.NE.0)GO TO 498
C.....MODIFICACAO P/ CALCULO DA CAM. DE ENSAIO & DIFUSOR SUPERSONICO
      NCE=NCE-NPTO-NDIF
      DO 497 I=NCE+1,NCE+NPTO+NDIF
      K=I-NCE
      SPXPA(K)=EML*SPX(I)
      SPFPA(K)=YEM
      IF(I.LE.NCE+NPTO)GO TO 494
      IZ=IGT(60,SPXPA(K),XSG)
      DX=SPXPA(K)-XSG(IZ-1)
      YTS=YSG(IZ-1)
      IF(SPXPA(K).LT.XSG(IZ))GO TO 493
      YTS=YSG(IZ)
      DX=SPXPA(K)-XSG(IZ)
 493  SPFPA(K)=YTS+DX*(ASPL(IZ-1,1)+DX*(ASPL(IZ-1,2)+DX*ASPL(IZ-1,3)))
 494  DELPA1(K)=EML*DEL1(I)
      SPFPA(K)=SPFPA(K)+DELPA1(K)
      SPX(I)=0.
      SPF(I)=0.
      MC(I)=0.
      TC(I)=0.
      PC(I)=0.
      CPC(I)=0.
      DEL1(I)=0.
 497  CONTINUE
C......................................................................
C****
C**** GLAETTUNG DER VERDRAENGTEN KONTUR
C****
  498 WRITE(6,1073)
 1073 FORMAT(1H1,'GLAETTUNG DER VERDRAENGTEN KONTUR'/)
      WRITE(6,1067)DYC
 1067 FORMAT(1H0,'VORGEGEBENE FEHLERABWEICHUNG DF =',F5.2,'  0/00'//)
      WRITE(6,1070)
 1070 FORMAT(1H0,'  N',9X,1HX,9X,1HY,3X,'Y GLATT',3X,'GLAETTUNG'/
     *43X,'0/0'/)
      DY=DYC/1000
      DO 150 J=1,NCE
      IF(IDELTA.EQ.0) SPF(J)=YC(J)/ELREF+DEL1(J)
      IF(IDELTA.EQ.1) SPF(J)=DEL1(J)-SPF(J)
  150 DF(J)=DY*SPF(J)
      N1=1
      N2=NCE
      S=N2-N1+1
      CALL SMOOTH(N1,N2,S)
      IF(IYM.EQ.0) EML=YMINM/A(1)
      IF(IYM.EQ.1) EML=YEM/A(NCE)
      IF(IDUESE.EQ.1) EML=(YMINM+DEL1(IGRT))/A(IGRT)
C     IF(IDUESE.EQ.1) EML=YMINM/A(IGRT)
      DO 100 J=1,NCE
      XC(J)=EML*SPX(J)
      YC(J)=EML*A(J)
      YCC=EML*SPF(J)
      DEL1(J)=EML*DEL1(J)
      SPX(J)=XC(J)
      SPF(J)=YC(J)
      DELYC=100*(YCC-YC(J))/YCC
  100 WRITE(6,1075)J,XC(J),YCC,YC(J),DELYC
 1075 FORMAT(1H ,I3,3F10.4,E12.4)
C****
C**** UNTERSCHALLKONTUR
C****
      CALL SUBCON(NCE,Y0M,ETA,ICS)
      NE=NCE
      NCS=ICS+1
      NCA=NCA+ICS
      NCD=NCD+ICS
      NCE=NCE+ICS
      DO 80 J=1,NCE
      XC(J)=SPX(J)
   80 YC(J)=SPF(J)
C****
C**** KRUEMMUNGSVERLAUF DER VERDRAENGTEN GEGLAETTETEN KONTUR
C****
      WRITE(6,1048)
 1048 FORMAT(1H1,'KRUEMMUNGSVERLAUF DER VERDRAENGTEN GEGLAETTETEN KONTUR
     *'//'   N',1HX,9X,1HY,'   KRUEMMUNG',25X,'GRAPHISCHE DARSTELLUNG DE
     *S KRUEMMUNGSVERLAUFS'/12X,2HMM,8X,4H1/MM/)
      N1=1
      N2=NCS
      DO 170 L=1,4
      IF(L.EQ.2) N2=NCA
      IF(L.EQ.3) N2=NCD
      IF(L.EQ.4) N2=NCE
      DO 160 J=N1,N2
      SPX(J)=XC(J)
  160 SPF(J)=YC(J)
      CALL SPLINE(N1,N2)
      CALL CURVAT(N1,N2)
      N1=N2
  170 WRITE(6,1099)
C1099 FORMAT(1H0)
C****
C**** INTERPOLATION DER GESAMTEN KONTUR UND STANZEN DER KONTURKOORDENATE
C****
      N1=1
      N2=NCE
      IF(IPUNCH.GT.0) CALL CONDA(N1,N2,DXC)
C****
C**** AUSDRUCKEN DER KONTURDATEN
C****
      DO 220 NTAB=1,1
      WRITE(6,1061)
 1061 FORMAT(1H1,'DATEN DER GESAMTEN KONTUR OHNE INTERPOLIERTE KOORDINAT
     *EN'/)
      WRITE(6,1049)
 1049 FORMAT(1H0,'  N',5X,2HXC,6X,7HYC+DEL1,4X,4HDEL1/9X,2HMM,8X,2HMM,
     *9X,2HMM/)
      IF(IDUESE.NE.0)GO TO 338
C.....MODIFICACAO P/ CALCULO DA CAMARA DE ENSAIO & DIF. SUPERSONICO
      DO 337 I=NCE+1,NCE+NPTO+NDIF
      K=I-NCE
      J=I-ICS
      XC(I)=SPXPA(K)
      YC(I)=SPFPA(K)-DEL1(NE)
      DEL1(J)=DELPA1(K)
 337  CONTINUE
      NCE=NCE+NPTO+NDIF
      NE=NE+NPTO+NDIF
C......................................................................
  338 DO 220 K=1,NCE
      J=NCE+1-K
      L=NE+1-K
      X=XC(J)-XC(1)
      Y=YC(J)
      YP1=Y0M-YC(J)
      IF(L.LT.1)GO TO 61
      D1=DEL1(L)
      WRITE(8,*)L,X,Y-D1,MC(L)
      IF(IDUESE.EQ.0)WRITE(6,1055)K,X,Y,YP1,D1
      IF(IDUESE.NE.0)THEN
	 WRITE(6,2054)K,X,Y,D1
       WRITE(9,2054)K,X,Y,-Y,D1
      ENDIF
 2054 FORMAT(1H ,I3,3F10.4,E12.4)
 1055 FORMAT(1H ,I3,3F10.4,E12.4)
      GO TO 220
   61 IF(IDUESE.EQ.0)WRITE(6,1056)K,X,Y,YP1
      IF(IDUESE.NE.0)THEN
	 WRITE(6,2057)K,X,Y
       WRITE(9,2057)K,X,Y,-Y
      ENDIF
      WRITE(8,*)K,X,Y,MC(J)
 1056 FORMAT(1H ,I3,3F10.4)
 2057 FORMAT(1H ,I3,3F10.4)
  220 CONTINUE
C**** DUMMY LINE
  555 STOP
      END
C**** HALL-KOEFFIZIENTEN FI
      SUBROUTINE HALL FI(Q,Y2,Y4,Y6,F0,F1,F2)
      COMMON /H/D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13,D14,D15,D16,
     *D17
      E1=D1*Y2+D2
      E2=D3*Y4+D4*Y2+D5
      E3=Y2+D6
      E4=D7
      E5=D8*Y6+D9*Y4+D10*Y2+D11
      E6=D12*Y4+D13*Y2+D14
      E7=D15*Y2+D16
      E8=D17
      F0=Q-E1
      F1=-(E2+F0*(E3+E4*F0))
      F2=-(E5+F0*(E6+F0*(E7+E8*F0))+F1*(E3+2*E4*F0))
      RETURN
      END
C**** COEFICIENTES FI MODIFICADOS:
      SUBROUTINE HALL FM(QQ,Y2,Y4,Y6,F0,F1,F2)
      COMMON /ZMZM/ZM01,ZM02,ZM03,ZM04,ZM05,ZM06,ZM07,ZM08,ZM09,
     1ZM10,ZM11,ZM12,ZM13,ZM14,ZM15,ZM16,ZM17
      F0=ZM01*Y2+ZM02+QQ
      F1=ZM03*Y2+ZM04+(ZM05*Y2+ZM06)*QQ+ZM07*QQ*QQ
      F2=ZM08*Y6+ZM09*Y4+ZM10*Y2+ZM11+(ZM12*Y4+ZM13*Y2+ZM14)*QQ+
     1(ZM15*Y2+ZM16)*QQ*QQ+ZM17*QQ**3
      RETURN
      END
C**** HALL-KOEFFIZIENTEN TI
      SUBROUTINE HALL TI(F0,F1,F2,Y,Y3,Y5,Y7,T1,T2,T3)
      COMMON /J/A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12,A13,A14,A15,A16,
     *A17
      B1=A1*Y3+A2*Y
      B2=Y
      B3=A3*Y5+A4*Y3+A5*Y
      B4=A6*Y3+A7*Y
      B5=A8*Y7+A9*Y5+A10*Y3+A11*Y
      B6=A12*Y5+A13*Y3+A14*Y
      B7=A15*Y3+A16*Y
      B8=A17*Y
      T1=B1+B2*F0
      T2=B2*F1+B3+B4*F0
      T3=B2*F2+B4*F1+B5+F0*(B6+F0*(B7+B8*F0))
      RETURN
      END
C**** HALL-KOEFFIZIENTEN QI
      SUBROUTINE HALL QI(Z,Y2,Y4,Y6,Q1,Q2,Q3)
      COMMON /H/D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13,D14,D15,D16,
     *D17
      E1=D1*Y2+D2
      E2=D3*Y4+D4*Y2+D5
      E3=Y2+D6
      E4=D7
      E5=D8*Y6+D9*Y4+D10*Y2+D11
      E6=D12*Y4+D13*Y2+D14
      E7=D15*Y2+D16
      E8=D17
      Q1=E1+Z
      Q2=E2+E3*Z+E4*Z**2
      Q3=E5+E6*Z+E7*Z**2+E8*Z**3
      RETURN
      END
C**** MODIFICACAO DA HALL QI
      SUBROUTINE HALL QM(ZS,Y2,Y4,Y6,Q1,Q2,Q3)
      COMMON /QMQM/QM01,QM02,QM03,QM04,QM05,QM06,QM07,QM08,QM09,QM10,
     1QM11,QM12,QM13,QM14,QM15,QM16,QM17
      Q1=QM01*Y2+QM02+ZS
      Q2=QM03*Y4+QM04*Y2+QM05+(Y2+QM06)*ZS+QM07*ZS*ZS
      Q3=QM08*Y6+QM09*Y4+QM10*Y2+QM11+(QM12*Y4+QM13*Y2+QM14)*ZS+
     1(QM15*Y2+QM16)*ZS*ZS+QM17*ZS**3
      RETURN
      END
      SUBROUTINE HALL TH(Z,Y,Y3,Y5,Y7,TH1,TH2,TH3)
      COMMON /J/T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15,T16,
     *T17
      TH1=T1*Y3+(T2+Z)*Y
      TH2=T3*Y5+T4*Y3+T5*Y+(T6*Y3+T7*Y)*Z
      TH3=T8*Y7+T9*Y5+T10*Y3+T11*Y+(T12*Y5+T13*Y3+T14*Y)*Z+(T15*Y3+T16*
     *Y)*Z**2+T17*Y*Z**3
      RETURN
      END
      SUBROUTINE FGAMM0(GAM)
      COMMON /H/D121,D101,D241,D221,D201,D202,D203,D361,D341,D321,D301,
     *D342,D322,D302,D323,D303,D304
      COMMON /J/A131,A111,A251,A231,A211,A232,A212,A371,A351,A331,A311,
     *A352,A332,A312,A333,A313,A314
      FKAP(A0,B0,C0,D0)=(A0*GAM**2+B0*GAM+C0)/D0
      D121=1./2.
      D101=-1./6.
      D241=FKAP(0.,1.,6.,18.)
      D221=-FKAP(0.,2.,9.,18.)
      D201=FKAP(0.,1.,30.,270.)
      D202=-1./2.
      D203=-FKAP(0.,2.,-3.,6.)
      D361=FKAP(362.,1629.,3357.,12960.)
      D341=-FKAP(194.,909.,1737.,2592.)
      D321=FKAP(854.,3867.,6939.,12960.)
      D301=-FKAP(782.,5523.,22887.,272160.)
      D342=FKAP(26.,75.,285.,288.)
      D322=-FKAP(26.,75.,213.,144.)
      D302=FKAP(134.,429.,1743.,4320.)
      D323=-FKAP(0.,3.,-2.,4.)
      D303=FKAP(0.,7.,-18.,36.)
      D304=FKAP(2.,-33.,9.,72.)
      A131=1./6.
      A111=-1./6.
      A251=FKAP(0.,22.,45.,360.)
      A231=-FKAP(0.,5.,15.,54.)
      A211=FKAP(0.,34.,165.,1080.)
      A232=FKAP(0.,2.,6.,9.)
      A212=-FKAP(0.,2.,6.,9.)
      A371=FKAP(6574.,19257.,18639.,181440.)
      A351=-FKAP(2254.,7929.,8739.,25920.)
      A331=FKAP(5026.,21639.,28917.,77760.)
      A311=-FKAP(7570.,42735.,74817.,544320.)
      A352=FKAP(362.,957.,1107.,2160.)
      A332=-FKAP(194.,609.,711.,648.)
      A312=FKAP(854.,3219.,3789.,6480.)
      A333=FKAP(26.,3.,-15.,144.)
      A313=-FKAP(26.,27.,9.,144.)
      A314=-FKAP(0.,1.,1.,2.)
      RETURN
      END
C**** ERMITTLUNG EINES CHARAKTERISTIKENPUNKTES N AUS DEN FELDPUNKTEN
C**** L UND R                   C H A R A  F
      SUBROUTINE CHARA F (XL,YL,THL,ML,XR,YR,THR,MR,XN,YN,THN,MN,NIT,ID,
     *G)
      REAL ML,MR,MN,MUEL,MUER,MUEN
      DATA EPSIT/1.E-10/
      IS=1
      NIT=5
      AUX=1./ML
      MUEL=ASIN(AUX)
      AUX=1./MR
      MUER=ASIN(AUX)
      A1L=TAN(THL+MUEL)
      A1R=TAN(THR-MUER)
      A2L=SQRT(ML**2-1.)/ML/(1.+(G-1.)/2.*ML**2)                        030
      A2R=SQRT(MR**2-1.)/MR/(1.+(G-1.)/2.*MR**2)                        040
      IF(ID.EQ.0)GO TO 1
      A3L=SIN(THL)/ML/YL/COS(THL+MUEL)
      A3R=SIN(THR)/MR/YR/COS(THR-MUER)
      GO TO 2
    1 A3L=0.
      A3R=0.
    2 A1LM=A1L
      A1RM=A1R
      A2LM=A2L
      A2RM=A2R
      A3LM=A3L
      A3RM=A3R
    6 XN=(YR-YL+A1LM*XL-A1RM*XR)/(A1LM-A1RM)
      YN1=YL+A1LM*(XN-XL)
      YN2=YR+A1RM*(XN-XR)
      YN=(YN1+YN2)/2.
      MN=(THR-THL+A2LM*ML+A2RM*MR+A3LM*(XN-XL)+A3RM*(XN-XR))/(A2LM+A2RM)
      IF(MN.LE.1)GO TO 7
      THN1=THL+A2LM*(MN-ML)-A3LM*(XN-XL)
      THN2=THR-A2RM*(MN-MR)+A3RM*(XN-XR)
      THN=(THN1+THN2)/2.
      IF(IS.LT.0)GO TO 12
      IS=-1
      GO TO 14
   12 IF(ABS(THNN-THN).LT.EPSIT.OR.NIT.EQ.0) RETURN
   14 NIT=NIT-1
      THNN=THN
      AUX=1./MN
      MUEN=ASIN(AUX)
      A1NL=TAN(THN+MUEN)
      A1NR=TAN(THN-MUEN)
      A2N=SQRT(MN**2-1.)/MN/(1.+(G-1.)/2.*MN**2)                        350
      IF(ID.EQ.0)GO TO 3
      ZN=SIN(THN)/MN/YN
      A3NL=ZN/COS(THN+MUEN)
      A3NR=ZN/COS(THN-MUEN)
      GO TO 4
    3 A3NL=0.
      A3NR=0.
    4 A1LM=(A1L+A1NL)/2.
      A1RM=(A1R+A1NR)/2.
      A2LM=(A2L+A2N)/2.
      A2RM=(A2R+A2N)/2.
      A3LM=(A3L+A3NL)/2.
      A3RM=(A3R+A3NR)/2.
      GO TO 6
    7 WRITE(6,101)
  101 FORMAT(1H0,'UNTERSCHALL IN CHARA F'//)
      STOP
      END
C**** ERMITTLUNG EINES CHARAKTERISTIKENPUNKTES N AUS EINEM FELDPUNKT R
C**** UND EINEM ACHSENPUNKT L  C H A R A  A
      SUBROUTINE CHARA A (XL,ML,XR,YR,THR,MR,XN,YN,THN,MN,NIT,ID,G)
      REAL ML,MR,MN,MUEL,MUER,MUEN
      DATA EPSIT/1.E-8/
      IS=1
      NIT=5
      AUX=1./ML
      MUEL=ASIN(AUX)
      AUX=1./MR
      MUER=ASIN(AUX)
      A1L=1./SQRT(ML**2-1.)                                             0
      A1R=TAN(THR-MUER)
      A2L=1./A1L/ML/(1.+(G-1.)/2.*ML**2)                                650
      A2R=SQRT(MR**2-1.)/MR/(1.+(G-1.)/2.*MR**2)                        660
      IF(ID.EQ.0)GO TO 1
      A3R=SIN(THR)/MR/YR/COS(THR-MUER)
      A3L=A3R
      GO TO 2
    1 A3L=0.
      A3R=0.
    2 A1LM=A1L
      A1RM=A1R
      A2LM=A2L
      A2RM=A2R
      A3LM=A3L
      A3RM=A3R
    6 XN=(YR+A1LM*XL-A1RM*XR)/(A1LM-A1RM)
      YN1=A1LM*(XN-XL)
      YN2=YR+A1RM*(XN-XR)
      YN=(YN1+YN2)/2.
      MN=(THR+A2LM*ML+A2RM*MR+A3LM*(XN-XL)+A3RM*(XN-XR))/(A2LM+A2RM)
      IF(MN.LE.1.)GO TO 7
      THN1=A2LM*(MN-ML)-A3LM*(XN-XL)
      THN2=THR-A2RM*(MN-MR)+A3RM*(XN-XR)
      THN=(THN1+THN2)/2.
      IF (IS.LT.0) GO TO 12
      IS=-1
      GO TO 14
   12 IF(ABS(THNN-THN).LT.EPSIT.OR.NIT.EQ.0) RETURN
   14 NIT=NIT-1
      THNN=THN
      AUX=1./MN
      MUEN=ASIN(AUX)
      A1NL=TAN(THN+MUEN)
      A1NR=TAN(THN-MUEN)
      A2N=SQRT(MN**2-1.)/MN/(1.+(G-1.)/2.*MN**2)                        970
      IF(ID.EQ.0)GO TO 3
      ZN=SIN(THN)/YN
      A3LM=4.*ZN/(ML+MN)/(COS(MUEL)+COS(THN+MUEN))
      A3NR=ZN/MN/COS(THN-MUEN)
      GO TO 4
    3 A3NL=0.
      A3NR=0.
    4 A1LM=(A1L+A1NL)/2.
      A1RM=(A1R+A1NR)/2.
      A2LM=(A2L+A2N)/2.
      A2RM=(A2R+A2N)/2.
      A3RM=(A3R+A3NR)/2.
      GO TO 6
    7 WRITE(6,101)
  101 FORMAT(1H0,'UNTERSCHALL IN CHARA A'//)
      STOP
      END
C**** ERMITTLUNG EINES CHARAKTERISTIKENPUNKTES N AUF DER SYMMETRIEACHSE
C**** AUS EINEM FELDPUNKT R    C H A R A  R
      SUBROUTINE CHARA R(XR,YR,THR,MR,XN,MN,NIT,ID,G)
      REAL MR,MN,MNN,MUER,MUEN
      DATA EPSIT/1E-8/
      IS=1
      NIT=5
      AUX=1./MR
      MUER=ASIN(AUX)
      A1R=TAN(THR-MUER)
      A2R=SQRT(MR**2-1.)/MR/(1.+(G-1.)/2.*MR**2)                        240
      IF(ID.EQ.0)GO TO 1
      A3R=SIN(THR)/MR/YR/COS(THR-MUER)
      GO TO 2
    1 A3R=0
    2 A1RM=A1R
      A2RM=A2R
      A3RM=A3R
    6 XN=XR-YR/A1RM
      MN=(THR+A2RM*MR+A3RM*(XN-XR))/A2RM
      IF(MN.LT.1.)GO TO 7
      IF(IS.LT.0)GO TO 12
      IS=-1
      GO TO 14
   12 IF(ABS(MNN-MN).LT.EPSIT.OR.NIT.EQ.0) RETURN
   14 NIT=NIT-1
      MNN=MN
      AUX=1./MN
      MUEN=ASIN(AUX)
      COMUEN=SQRT(MN**2-1.)
      A1NR=-1./COMUEN
      A2N=COMUEN/MN/(1.+(G-1.)/2.*MN**2)                                40
      A1RM=(A1NR+A1R)/2.
      A2RM=(A2N+A2R)/2.
      IF(ID.EQ.0)GO TO 6
      A3RM=4.*SIN(THR)/YR/(MR+MN)/(COS(THR-MUER)+COS(MUEN))
      GO TO 6
    7 WRITE(6,101)
  101 FORMAT(1H0,'UNTERSCHALL IN CHARA R'//)
      STOP
      END
C**** STROMLINIENPUNKT
      SUBROUTINE STREAM(IS,XA,YA,YAX,XJ,YJ,THJ,XK,YK,THK,XB,YB,YBX)
      IS=1
      YJX=TAN(THJ)
      YKX=TAN(THK)
      DXK=XK-XJ
      DYK=YK-YJ
      QY=DYK/DXK
      QYX=(YKX-YJX)/DXK
      S=YAX+YJX-QYX*(XJ+XA)-2.*QY
      T2=2.*(YA-YJ+XJ*QY)
      IF(ABS(QYX).LT.1.E-12)GO TO 2
      T1=XA*(QYX*XJ-YJX-YAX)
      T=T1+T2
      DIS=S/2./QYX
      SQ=SQRT(DIS**2-T/QYX)
      XB=SQ-DIS
      IF(DIS.LT.0.) XB=-SQ-DIS
      GO TO 5
    2 XB=(XA*(YAX+YJX)-T2)/(YAX+YJX-2.*QY)
      YBX=YJX
      GO TO 4
    5 YBX=YJX+QYX*(XB-XJ)
    4 YB=YJ+QY*(XB-XJ)
      DXB=XB-XJ
      DYB=YB-YJ
      IF(SQRT(DXB**2+DYB**2).LE.SQRT(DXK**2+DYK**2))GO TO 6
      GO TO 7
    6 IF(DXK.GE.0.AND.DYK.LT.0.AND.DXB.GE.0.AND.DYB.LT.0)GO TO 8
      IF(DXK.GE.0.AND.DYK.GE.0.AND.DXB.GE.0.AND.DYB.GE.0)GO TO 8
      IF(DXK.LT.0.AND.DYK.GE.0.AND.DXB.LT.0.AND.DYB.GE.0)GO TO 8
      IF(DXK.LT.0.AND.DYK.LT.0.AND.DXB.LT.0.AND.DYB.LT.0)GO TO 8
    7 IS=-1
    8 RETURN
      END
C***  MACHZAHLVERLAUF UND DESSEN ABLEITUNGEN AUF DER ACHSE
      SUBROUTINE AX(ZX,ZA,QA,D1,D2,D3,D4,D5,RC,RC2,RC3,F0,F1,F2,MA,
     *MAX,ISAUER)
      REAL MA,MAX,MQ
      ZA2=ZA**2
      ZA3=ZA**3
      QA2=QA**2
      MA=SQRT(F1*QA2/(F2-QA2))
      Q1Z=1.
      Q2Z=D1+2.*D2*ZA
      Q3Z=D3+2.*D4*ZA+3.*D5*ZA2                                         0
      QAX=(Q1Z/RC+Q2Z/RC2+Q3Z/RC3)*ZX
      IF (ISAUER.EQ.1) QAX=Q1Z/RC*ZX
      MQ=MA/QA*F2/(F2-QA2)
      MAX=MQ*QAX
      RETURN
      END
C**** TEMPERATUREN, DRUCKE UND CP-WERTE DER KONTUR
      SUBROUTINE TECEPE(FG0,FG1,TR,MR,PR,N,G)
      REAL M,MR
      COMMON /F/M(200),T(200),P(200),C(200)
      FM1=FG1+MR**2
      FM2=2./G/MR**2
      TTR=FM1/(FG1+M(N)**2)
      PPR=TTR**FG0
      T(N)=TTR*TR
      P(N)=1.E-5*PPR*PR
      C(N)=FM2*(PPR-1.)
      RETURN
      END
C**** KOEFFIZIENTEN DER INTERPOLIERENDEN SPLINE-FUNKTION
      SUBROUTINE SPLINE(N1,N2)
      COMMON /K/X(200),F(200)
      COMMON /G/A(200),B(200),C(200),D(200)
      M1=N1+1
      M2=N2-1
      S=0.
      DO 10 I=N1,M2
      D(I)=X(I+1)-X(I)
      R=(F(I+1)-F(I))/D(I)
      C(I)=R-S
   10 S=R
      S=0.
      R=0.
      C(N1)=0.
      C(N2)=0.
      DO 20 I=M1,M2
      C(I)=C(I)+R*C(I-1)
      B(I)=2.*(X(I-1)-X(I+1))-R*S
      S=D(I)
   20 R=S/B(I)
      K=N2-M1
      DO 30 J=1,K
      I=N2-J
   30 C(I)=(D(I)*C(I+1)-C(I))/B(I)
      DO 40 I=N1,M2
      S=D(I)
      R=C(I+1)-C(I)
      D(I)=R/S
      C(I)=3.*C(I)
      B(I)=(F(I+1)-F(I))/S-S*(C(I)+R)
   40 A(I)=F(I)
      A(N2)=F(N2)
      RETURN
      END
C***  KOEFFIZIENTEN DER GLAETTENDEN SPLINE-FUNKTION
      SUBROUTINE SMOOTH(N1,N2,S)
      COMMON /K/X(200),F(200)
      COMMON /G/A(200),B(200),C(200),D(200)
      COMMON /I/DF(200)
      DIMENSION R(200),R1(200),R2(200),T(200),T1(200),U(200),V(200)
      I=N2+1
      R(N1)=0.
      R(N2)=0.
      R(I)=0.
      R1(N1)=0.
      R1(N2)=0.
      R2(N1)=0.
      U(N1)=0.
      U(N2)=0.
      U(I)=0.
      P=0.
      M1=N1+1
      M2=N2-1
      H=X(M1)-X(N1)
      F1=(F(M1)-F(N1))/H
      DO 10 I=M1,M2
      G=H
      H=X(I+1)-X(I)
      E=F1
      F1=(F(I+1)-F(I))/H
      A(I)=F1-E
      T(I)=2.*(G+H)/3.
      T1(I)=H/3.
      R(I)=DF(I-1)/G
      R2(I)=DF(I+1)/H
   10 R1(I)=-DF(I)*(1./G+1./H)                                          0
      DO 20 I=M1,M2
      B(I)=R(I)**2+R1(I)**2+R2(I)**2
      C(I)=R1(I)*R(I+1)+R2(I)*R1(I+1)
   20 D(I)=R2(I)*R(I+2)
      F2=-S
    1 DO 30 I=M1,M2
      IF(I.EQ.M1)GO TO 3
      R(I)=1./SQRT(P*B(I)+T(I)-R1(I-1)**2-R2(I-2)**2)
      GO TO 4
    3 R(I)=1./SQRT(P*B(I)+T(I)-R1(I-1)**2)
    4 R1(I)=(P*C(I)+T1(I)-R1(I-1)*R2(I-1))*R(I)
   30 R2(I)=P*D(I)*R(I)
      DO 40 I=M1,M2
      IF(I.EQ.M1)GO TO 5
      U(I)=(A(I)-R1(I-1)*U(I-1)-R2(I-2)*U(I-2))*R(I)
      GO TO 40
    5 U(I)=(A(I)-R1(I-1)*U(I-1))*R(I)
   40 CONTINUE
      K=N2-M1
      DO 50 J=1,K
      I=N2-J
   50 U(I)=(U(I)-R1(I)*U(I+1)-R2(I)*U(I+2))*R(I)
      E=0.
      H=0.
      DO 60 I=N1,M2
      G=H
      H=(U(I+1)-U(I))/(X(I+1)-X(I))
      H1=H-G
      V(I)=H1*DF(I)**2
   60 E=E+V(I)*H1
      V(N2)=-H*DF(N2)**2
      G=V(N2)
      E=E-G*H
      G=F2
      F2=E*P**2
      IF(F2.GE.S.OR.F2.LE.G)GO TO 2
      F1=0.
      H=(V(M1)-V(N1))/(X(M1)-X(N1))
      DO 70 I=M1,M2
      G=H
      H=(V(I+1)-V(I))/(X(I+1)-X(I))
      IF(I.EQ.M1)GO TO 6
      R(I)=(H-G-R1(I-1)*R(I-1)-R2(I-2)*R(I-2))*R(I)
      GO TO 70
    6 R(I)=(H-G-R1(I-1)*R(I-1))*R(I)
   70 F1=F1+R(I)**2
      H=E-P*F1
      IF(H.LE.0.)GO TO 2
      P=P+(S-F2)/((SQRT(S/E)+P)*H)
      GO TO 1
    2 DO 80 I=N1,N2
      A(I)=F(I)-P*V(I)
   80 C(I)=U(I)
      DO 90 I=N1,M2
      H=X(I+1)-X(I)
      D(I)=(C(I+1)-C(I))/(3.*H)
   90 B(I)=(A(I+1)-A(I))/H-(H*D(I)+C(I))*H
      RETURN
      END
C**** LAGRANGE'SCHE INTERPOLATION
      SUBROUTINE LAGRAN(XE,F,N)
      REAL MF
      COMMON /A/PNF(150,2),XF(150,2),YF(150,2),THF(150,2),MF(150,2),
     *NI(150,2)
      DIMENSION EL(150)
      J=N
      IF(N.GT.4) J=4
      L=1+N-J
      F=0.
      DO 10 I=L,N
      EL(I)=1.
      DO 20 K=L,N
      IF(I.EQ.K)GO TO 20
      EL(I)=EL(I)*(XE-XF(K,2))/(XF(I,2)-XF(K,2))
   20 CONTINUE
   10 F=F+EL(I)*MF(I,2)
      RETURN
      END
C**** FORTSETZUNGSPROBLEM
      SUBROUTINE TAYLOR(N,G)
      COMMON /K/A(200),E(200)
      HN1=A(N)-A(N-1)
      HN2=A(N)-A(N-2)
      HN3=A(N)-A(N-3)
      FN=E(N)
      F1=E(N-1)
      F2=E(N-2)
      F3=E(N-3)
      G=((HN2*HN3)**2*(FN-F1)+(HN3*HN1)**2*(FN-F2)-2.*(HN1*HN2)**2*
     *(FN-F3))/(HN1*HN2*HN3*(HN2*HN3+HN3*HN1-2.*HN1*HN2))
      RETURN
      END
C*****PRANDTL-MEYER WINKEL
      SUBROUTINE PRAMEY(F2,F3,NR,N)
      COMMON /A/PNF(150,2),XF(150,2),YF(150,2),THF(150,2),MF(150,2),
     *NI(150,2)
      REAL MF,N
      FM=MF(NR,2)**2-1.
      N=F3*ATAN(SQRT(FM/F2))-ATAN(SQRT(FM))
      RETURN
      END
C*****AUSDRUCKEN DER FELDPUNKTE
      SUBROUTINE RESULT(N,IPL)
      INTEGER PNF
      REAL MF
      COMMON /A/PNF(150,2),XF(150,2),YF(150,2),THF(150,2),MF(150,2),
     *NI(150,2)
      DATA D/57.2957795/
      WRITE(6,111)
  111 FORMAT(1H0)
      KL=1
      DO 10 L=1,N
      PNF(L,1)=PNF(L,2)
      XF(L,1)=XF(L,2)
      YF(L,1)=YF(L,2)
      THF(L,1)=THF(L,2)
      THFD=D*THF(L,1)
      MF(L,1)=MF(L,2)
      NI(L,1)=NI(L,2)
      KL=-KL
      IF(KL.EQ.1)GO TO 3
      WRITE(6,101)PNF(L,1),XF(L,1),YF(L,1),THFD,MF(L,1),NI(L,1)
  101 FORMAT(1H ,I5,4F10.4,I5)
      GO TO 4
    3 WRITE(6,106)PNF(L,1),XF(L,1),YF(L,1),THFD,MF(L,1),NI(L,1)
  106 FORMAT(1H+,70X,I5,4F10.4,I5)
    4 IF(L.EQ.1)GO TO 1
      GO TO 10
    1 WRITE(6,107)N
  107 FORMAT(1H+,50X,I5)
   10 CONTINUE
      IF(IPL.LE.0) RETURN
      WRITE(10,102)N
  102 FORMAT(I5)
      WRITE(10,105)(PNF(L,1),XF(L,1),YF(L,1),L=1,N)
  105 FORMAT(3(I5,2F9.4))
      RETURN
      END
C*****INTERPOLATION UND STANZEN DER KONTUROORDINATEN
      SUBROUTINE CONDA(N1,N2,DXC)
      COMMON /D/XC(200),YC(200),YCX(200)
      COMMON /G/A(200),B(200),C(200),D(200)
      DX=0.2*DXC
      L=0
      M=0
      DO 10 K=N1,N2
      M=M+1
      J=N2+1-M
      L=L+1
      XIN=XC(N2)-XC(J)
      YIN=YC(J)
      IS=1
      GO TO 1
    3 IF(K.EQ.N2)GO TO 5
    2 X=XIN+DXC
      XX=X+DX
      IF(XX.GE.XC(N2)-XC(J-1))GO TO 10
      L=L+1
      XIN=X
      H=XC(N2)-XC(J-1)-X
      YIN=((D(J-1)*H+C(J-1))*H+B(J-1))*H+A(J-1)
      IS=2
    1 IF(L.EQ.1) X1=XIN
      IF(L.EQ.1) YI=YIN
      IF(L.EQ.2) X2=XIN
      IF(L.EQ.2) Y2=YIN
      IF(L.EQ.3)GO TO 12
      IF(IS-2) 3,2,2
   12 L=0
      WRITE(10,101)X1,Y1,X2,Y2,XIN,YIN
  101 FORMAT(6F10.4)
      IF(IS-2) 3,2,4
    5 IF(L.EQ.0) RETURN
      IS=3
      IF(L-2) 8,9,12
    8 WRITE(10,101)X1,Y1
      RETURN
    9 WRITE(10,101)X1,Y1,X2,Y2
    4 RETURN
   10 CONTINUE
      END
C*****UNTERSCHALLKONTUR
      SUBROUTINE SUBCON(N,Y0,E,I)
      COMMON /K/X(200),Y(200)
      COMMON /G/A(200),B(200),C(200),D(200)
      I=N/3
      IF(I.GT.30) I=30
      DLX=E*(Y0-Y(1))
      DELX=X(N)-X(1)
      XT=AINT(DELX+DLX)
      X0=X(N)-XT
      DO 10 K=1,N
      J=N+1-K
      L=J+I
      X(L)=X(J)
   10 Y(L)=Y(J)
      X(I)=X0
      Y(I)=Y0
      N1=I
      N2=I+6
      CALL SPLINE(N1,N2)
      X(1)=X0
      Y(1)=Y0
      DX=ABS((X(I+1)-X0)/I)
      XS=X0
      DO 20 K=2,I
      XS=XS+DX
      X(K)=XS
      H=X(K)-X(1)
   20 Y(K)=((D(I)*H+C(I))*H+B(I))*H+A(I)
      RETURN
      END
C*****KRUEMMUNGSVERLAUF
      SUBROUTINE CURVAT(N1,N2)
      COMMON /D/XC(200),YC(200),YCX(200)
      COMMON /G/A(200),B(200),C(200),D(200)
      COMMON /K/X(200),FX(200)
      DIMENSION CUR(200),ISI(200),COL(10)
      DATA IBLANC,ISTAR,IPONT/1H ,1H*,1H./,DELI/.01/,LI/25/
      M1=N1+1
      M2=N2-1
      DO 10 IZ=N1,M2
      QN=SQRT(1+B(IZ)**2)
   10 CUR(IZ)=2*C(IZ)/QN**3
      CUR(N1)=CUR(M1)
      IF(N2-N1.LT.6)GO TO 7
      NP=N1+2
      DO 60 LP=1,3
      X(LP)=XC(NP+LP)
   60 FX(LP)=CUR(NP+LP)
      CALL PAREX(3,A0,A1,A2)
      DO 70 LP=1,3
      KP=N1+3-LP
   70 CUR(KP)=(A2*XC(KP)+A1)*XC(KP)+A0
    7 CUR(N2)=CUR(M2)
      IF(N2-N1.LT.6)GO TO 8
      NP=N2-6
      DO 80 LP=1,3
      X(LP)=XC(NP+LP)
   80 FX(LP)=CUR(NP+LP)
      CALL PAREX(3,A0,A1,A2)
      DO 90 LP=1,3
      KP=N2-3+LP
   90 CUR(KP)=(A2*XC(KP)+A1)*XC(KP)+A0
    8 CV=ABS(CUR(N1))
      DO 100 J=M1,N2
      FX(J)=ABS(CUR(J))
      IF(FX(J).LE.CV)GO TO 100
      CV=FX(J)
  100 CONTINUE
      DLI=0.
    6 DLI=DLI+DELI
      IF(DLI.GE.CV)GO TO 2
      GO TO 6
    2 DLI=DLI/2
      CO=-3.*DLI
      DO 50 I=1,6
      CO=CO+DLI
   50 COL(I)=CO
      WRITE(6,102)(COL(N),N=1,5)
  102 FORMAT(1H0,20X,5(14X,F8.4)//)
      Q=DLI/LI
      DO 20 K=N1,N2
      N=1
      DO 30 J=1,118
      IF(J.EQ.N)GO TO 4
      ISI(J)=IBLANC
      GO TO 30
    4 ISI(J)=IPONT
      N=N+LI
   30 CONTINUE
      CU=CUR(K)
      DO 40 L=1,101
      HL=L
      IF(CU.LT.0.) HL=-L
      CC=HL*Q
      IF(CU.LT.0.)GO TO 3
      IF(CC.LT.CU)GO TO 40
      M=50+L
      GO TO 5
    3 IF(CC.GT.CU)GO TO 40
      M=51-L
      GO TO 5
   40 CONTINUE
    5 ISI(M)=ISTAR
      WRITE(6,104)K,XC(K),YC(K),CUR(K)
  104 FORMAT(1H ,I3,2F10.4,E12.4,5X)
   20 CONTINUE
      RETURN
      END
C*****PARABELEXTRAPOLATION
      SUBROUTINE PAREX(N,A0,A1,A2)
      COMMON /K/X(200),Y(200)
      X12=X(N-2)-X(N-1)
      X31=X(N)-X(N-2)
      XX21=X(N-1)**2-X(N-2)**2
      XX13=X(N-2)**2-X(N)**2
      Y12=Y(N-2)-Y(N-1)
      Y13=Y(N-2)-Y(N)
      A2=(Y13+Y12*X31/X12)/(XX13-XX21*X31/X12)
      A1=(Y12+A2*XX21)/X12
      A0=Y(N)-A1*X(N)-A2*X(N)**2
      RETURN
      END
C BERECHNUNG DER LAMINAREN UND TURBULENTEN GRENZSCHICHT
C IN KOMPRESSIBLER STROEMUNG OHNE WAERMEUEBERGANG
C INTEGRALSAETZE FUER IMPULS- UND KINETISCHE ENERGIE
C VERFASSER H.SENKBEIL
C BESCHREIBUNG IN: J.C.ROTTA, FORTRAN IV - RECHENPROGRAMM FUER GRENZ-
C SCHICHTEN BEI KOMPRESSIBLEN EBENEN UND ACHSENSYMMETRISCHEN
C STROEMUNGEN. DLR-FB  71-51  DFVLR-AVA,  INSTITUT FUER STROEMUNGS-
C MECHANIK, GOETTINGEN 1971
      SUBROUTINE B LAYER(D2,D3,M,EMAREF,TREF,REREF,GAMM)
      REAL MC
      COMMON /B/THETA,OMEGA,S,EMAI(200),PA(10),PA1(10)
      COMMON /C/D1(200)
      COMMON /K/XI(200),RI(200)
      COMMON /F/MC(200),TC(200),PC(200),CPI(200)
      IFALK=0
      IF(D2.EQ.0.) IFALK=1
      KU=1
      X0=XI(1)-(XI(M)-XI(1))/100
      XU=(X0+XI(1))/2
      XEND=XI(M)
      TH0=0.
      D21=0.
      D22=0.
      D23=0.
      D31=0.
      D32=0.
      D33=0.
      DN1=0.
      DN2=0.
      DX=0.
      DN3=0.
      FN1=0.
      FN2=0.
      FN3=0.
      DELT22=0
      DELT23=0
      DO 999  I=1,10
      PA(I)=0.
  999 PA1(I)=0.
      NI=1
C STOFFGROESSEN FUER LUFT
      THETA=110.
      OMEGA=1.5
      P0=1.
      MZW=16**5
      MGT=M
      S=THETA/TREF
      RHOU=1.
      T=1.
      TD=1.
      EMA=0.
      CALL HCFL (GAMM,D2,D3,REREF,X,RHOU,EMA,TD,R,H12,H32,H12U,CF,CD,
     *NAB,NRL)
      CALL HCFT (GAMM,D2,D3,REREF,X,RHOU,EMA,TD,R,H12,H32,H12U,CF,CD,
     *NAB,NR)
      CALL UKR(KU,X,XU,D2,D3,REREF,RHOU,EMA,T,R,H12,CD,NAB,XE,GAMM,NRU)
      WRITE(6,24)    GAMM,OMEGA,THETA,P0,NRL,NR
      WRITE(6,1024) X0,TH0,XEND,MZW,NRU,XU,KU
      I=1
      HN=0.
      HN1=0.
      X=XI(I)
      XE=X
      R=RI(I)
      CP=CPI(I)
      CALL IS(CP,EMAREF,Q,EMA,RHO,U,T,RHOU,P0,GAMM)
      IW=0
      IF(D2.LE.0.)GO TO 8000
      IF(X.LT.XU.OR.KU.NE.1)GO TO 25
 1004 CALL HCFT (GAMM,D2,D3,REREF,X,RHOU,EMA,T,R,H12,H32,H12U,CF,CD,
     *NAB,NR)
      XU=XU-ABS(XEND)
      WRITE(6,5000)
      GO TO 26
   25 CALL HCFL (GAMM,D2,D3,REREF,X,RHOU,EMA,TD,R,H12,H32,H12U,CF,CD,
     *NAB,NR)
      CALL UKR(KU,X,XU,D2,D3,REREF,RHOU,EMA,T,R,H12,CD,NAB,XE,GAMM,NRU)
      IF(KU.EQ.1.AND.XU.LE.X)GO TO 1004
      WRITE(6,5001)
   26 CQ=CF/2.*Q*R
      IF(H12.EQ.0.)GO TO 8001
      FN=CQ*U**H12
      DN=CD*Q*R*U*RHO**(1.-GAMM)
C DATEN BIS ZUR NAECHSTEN UNSTETIGKEIT ODER ENDE
 2000 IUN=I-1
      MGT=MGT-IUN
      DO 6001  J=1,MGT
      JIUN=J+IUN
      XI(J)=XI(JIUN)
      RI(J)=RI(JIUN)
      CPI(J)=CPI(JIUN)
      IF(EMAREF.LT.0.1) EMAI(J)=SQRT(1.-CPI(J))
      IF(EMAREF.GE.0.1)CALL IS(CPI(J),EMAREF,QIS,EMA1,RHOIS,UIS,TIS,
     *RHOUIS,P0,GAMM)
      IF(EMAREF.GE.0.1) EMAI(J)=EMA1/EMAREF
      IF(J.EQ.1)GO TO 6001
      IF(XI(J).EQ.XI(J-1))GO TO 6002
      IF(XI(J).LT.XI(J-1))GO TO 8111
 6001 CONTINUE
      M=MGT
      GO TO 6543
 6002 M=J-1
 6543 I=1
C HAUPTROUTINE
   29 D1(I)=H12*D2
      WRITE(6,27) I,XE,D2,D3,CF,H12,H12U,H32,CP,EMA,R,D1(I)
      IF(NAB.GT.0)GO TO 1766
      GO TO 77
 1766 WRITE(6,5005)
   77 IF(X.NE.XU.OR.KU.NE.1)GO TO 28
      CALL HCFT (GAMM,D2,D3,REREF,X,RHOU,EMA,T,R,H12,H32,H12U,CF,CD,
     *NAB,NR)
      IF(NAB.GT.0.OR.H12.EQ.0.)GO TO 2471
      IF(XU.LT.XI(I)) I=I-1
      XU=XU-ABS(XEND)
      GO TO 29
   28 IF(NAB.GT.0.AND.KU.NE.1)GO TO 1002
      IF(X.GE.XEND)GO TO 1001
C NEUER WERT  XE
      XE1=XE
   34 IF(XE1.GE.XU.OR.XU.GE.XI(I+1))GO TO 31
      XE=XU
      CALL INPOL(M,CPE,RE,XE,I,NI,EMAREF,P0,GAMM)
      GO TO 32
 1001 WRITE(6,5002)
      IF(X.LE.XU.OR.KU.EQ.1) RETURN
 1002 WRITE(6,1003)
      RETURN
   31 I=I+1
      IF(I.LE.MGT)GO TO 33
      XE=XEND
      CALL INPOL(M,CPE,RE,XE,I,NI,EMAREF,P0,GAMM)
      GO TO 32
   33 XE=XI(I)
      CPE=CPI(I)
      RE=RI(I)
      IF(DX.EQ.0..AND.XE-XE1.EQ.0.)GO TO 6003
   32 DX=XE-XE1
      IF(DX.NE.0.)GO TO 6003
C UNSTETIGKEIT IN CP
      IF(EMA.LE.1..OR.CPE.LE.CPI(I-1))GO TO 6004
      P2P1=(GAMM/2.*EMAREF**2*CPE+1.)/(GAMM/2.*EMAREF**2*CPI(I-1)+1.)
      EM12=P2P1*(GAMM+1.)+GAMM-1.
      IF(EM12.LE.GAMM*2.*EMA**2)GO TO 9010
      P2P1=(2.*GAMM*EMA**2-(GAMM-1.))/(GAMM+1.)
      CPE=(P2P1*(GAMM/2.*EMAREF**2*CPI(I-1)+1.)-1.)*2./(GAMM*EMAREF**2)
      CPI(I)=CPE
 9010 PT2PT1=(((GAMM+1.)*P2P1+GAMM-1.)/((GAMM-1.)*P2P1+GAMM+1.))
     ***(GAMM/(GAMM-1.))/P2P1**(1./(GAMM-1.))
      P0=P0*PT2PT1
 6004 WRITE(6,6005) P0
C NEUER SCHRITT
 6003 N=1
      K=1
C UMSPEICHERN
   53 X1=X
      EN=N
      HN2=HN1
      HN1=HN
      D23=D22
      D22=D21
      D21=D2
      D33=D32
      D32=D31
      D31=D3
      CQ1=CQ
      H121=H12
      Q1=Q
      R1=R
      RHOU1=RHOU
      RHO1=RHO
      U1=U
      DN3=DN2
      DN2=DN1
      DN1=DN
      FN3=FN2
      FN2=FN1
      FN1=FN
      DO 6006 IK=1,10
 6006 PA1(IK)=PA(IK)
      IW=IW+1
      IF(DX.EQ.0.) IW=1
      IF(K.LT.N)GO TO 35
      X=XE
      CP=CPE
      R=RE
      GO TO 36
C SCHRITTUNTERTEILUNG
   42 N=N+N
      IF(N.GT.16**5)GO TO 7001
      K=K+K-1
   35 EK=K
      EN=N
      IF(DX.EQ.0.)GO TO 6007
      X=XE1+EK/EN*DX
      CALL INPOL(M,CP,R,X,I,NI,EMAREF,P0,GAMM)
   36 HN=DX/EN
      CALL IS(CP,EMAREF,Q,EMA,RHO,U,T,RHOU,P0,GAMM)
      IF(IW.EQ.1)GO TO 37
      IF(IW.EQ.2)GO TO 38
C 3. UND FOLGENDER SCHRITT
      BETA0=(HN*(HN+2.*HN1+HN2))/(HN1*(HN1+HN2))
      BETA1=-((HN*(HN+HN1+HN2))/(HN2*HN1))
      BETA2=((HN+HN1)*HN)/(HN2*(HN1+HN2))
   40 A0=(2.*HN+3.*HN1)/(HN+HN1)
      A1=(HN+3.*HN1)/HN1
      A2=-(HN**2/((HN+HN1)*HN1))
      GO TO 41
 6007 CP=CPI(I-1)+EK/EN*(CPE-CPI(I-1))
      GO TO 36
C ERSTER SCHRITT
   37 A0=3.
      A1=3.
      A2=0.
      BETA0=0.
      BETA1=0.
      BETA2=0.
      GO TO 41
C ZWEITER SCHRITT
   38 BETA0=HN/HN1
      BETA1=-(HN/HN1)
      BETA2=0.
      GO TO 40
C EXTRAPOLATION FUER ERSTE NAEHERUNG
   41 D2S=D21*BETA0+D22*BETA1+D23*BETA2+D21
      D3S=D31*BETA0+D32*BETA1+D33*BETA2+D31
      J=0
C BEDINGUNGEN F. VERKLEINERUNG DER SCHRITTWEITE (A)
   48 CONTINUE
      IF(J.GT.8.OR.D2S.LE.0.OR.D3S.LE.0)GO TO 42
C FORMPARAMETER'  CF' CD
      IF(XU.LT.X.AND.KU.EQ.1)GO TO 44
      CALL HCFL(GAMM,D2S,D3S,REREF,X,RHOU,EMA,T,R,H12,H32,H12U,CF,CD,
     *NAB,NR)
      GO TO 45
   44 CALL HCFT(GAMM,D2S,D3S,REREF,X,RHOU,EMA,T,R,H12,H32,H12U,CF,CD,
     *NAB,NR)
C BEDINGUNG F. VERKLEINERUNG DER SCHRITTWEITE (B)
   45 CONTINUE
      IF(IFALK.EQ.1.AND.IW.EQ.1)GO TO 66
      IF(ABS(H12-H121).GT.0.05.OR.ABS((D2S-D21)/D21).GT.0.05)GO TO 42
C QUADRATUR
   66 D2Q1=D21*R1*Q1
      D3U1=D31*Q1*U1*RHO1**(1.-GAMM) *R1
      H12S=(H12+H121)/2.
      CQ=CF/2.*Q*R
      DN=CD*Q*R*U*RHO**(1.-GAMM)
      FN=CQ*U**H12S
      FI=HN/6.*(FN*A0+FN1*A1+FN2*A2)
      DI=HN/6.*(DN*A0+DN1*A1+DN2*A2)
      D2=(D2Q1*U1**H12S+FI)/(Q*R*U**H12S)
      D3=(D3U1+DI)/(Q*R*U*RHO**(1.-GAMM))
      DELT12=DELT22
      DELT13=DELT23
      DELT22=D2-D2S
      DELT23=D3-D3S
      AA=DELT22
      BB=DELT23
      CC=D2*0.0001
      DD=D3*0.0001
      IF(ABS(AA).LT.CC.AND.ABS(BB).LT.DD)GO TO 46
      DK12=D2S
      DK13=D3S
      D2S=D2
      D3S=D3
      IF(J.GT.0)GO TO 47
   49 J=J+1
      GO TO 48
C BEDINGUNGEN F. VERKLEINERUNG DER SCHRITTWEITE
   47 IF(ABS(DELT12-DELT22).LT.ABS(CC))GO TO 65
      D2S=D2+(D2-DK12)/(DELT12-DELT22)*DELT22
   65 IF(ABS(DELT13-DELT23).LT.ABS(DD))GO TO 49
      D3S=D3+(D3-DK13)/(DELT13-DELT23)*DELT23
      GO TO 49
   46 CALL UKR(KU,X,XU,D2,D3,REREF,RHOU,EMA,T,R,H12,CD,NAB,XE,GAMM,NRU)
C KRITERIUM FUER ABLOESUNG
      IF(NAB.EQ.0)GO TO 51
      DEL1=H12*D2
      WRITE(6,55) X,D2,D3,CF,H12,H12U,H32,CP,EMA,R,DEL1
      IF(X.GT.XU.AND.KU.EQ.1)GO TO 2471
      XU=X
      XE=X
      GO TO 29
 2471 D1(I)=H12*D2
      WRITE(6,55) X,D2,D3,CF,H12,H12U,H32,CP,EMA,R,D1(I)
      WRITE(6,5006)
      GO TO 1000
   51 IF(XU.EQ.X)GO TO 29
      IF(K.LT.N)GO TO 52
      IF(DX.NE.0.)GO TO 29
      IW=0
      GO TO 2000
   52 IF(DX.EQ.0.)GO TO 5553
C DRUCKEN DER ZWISCHENERGEBNISSE
      NZW=1
      IF(MZW.GE.0.AND.MZW.LT.N) NZW=N/(MZW+1)
      KA=K
 5554 KA=KA-NZW
      IF(KA.EQ.0)GO TO 5552
      IF(KA.GT.0)GO TO 5554
      GO TO 5555
 5552 DEL1=H12*D2
      WRITE(6,55) X,D2,D3,CF,H12,H12U,H32,CP,EMA,R,DEL1
C SCHRITTVERGROESSERUNG
 5555 IF((-1.)**K.GT.0.)GO TO 7000                                      0
 5553 K=K+1
      GO TO 53
 7000 N=N/2
      K=K/2
      GO TO 5553
C BERECHNUNG DER ANFANGSWERTE
 8000 IF(XI(1).LE.X0)GO TO 8002
      BETA=TH0/90.
      ZETA=1.
      IF(RI(1).GT.10.*XI(1))GO TO 8003
      BETA=(0.1+0.4*TH0/90.)*BETA
      ZETA=SQRT(1./3.)
 8003 IF(EMA.GT.1.) BETA=0.
      ZETA=ZETA*SQRT(2.-BETA)
      DELT2=0.47-0.3*BETA+0.12*BETA**2
      H32=1.625-0.0171*(1.-BETA)-0.0359*(1.-BETA)**3
      YM=T**OMEGA*(1.+S)/(T+S)
      D2=DELT2*ZETA*SQRT((XI(1)-X0)/RHOU*YM)/SQRT(REREF)
      D3=H32*D2
      IF(XU.LT.XI(1))  XU=XI(1)
      GO TO 25
C FEHLERMELDUNGEN
 8110 WRITE(6,9000)     CPMAX,CPMIN
      GO TO 1000
 8111 JIUN1=JIUN-1
      WRITE(6,8005) JIUN1,XI(JIUN1),JIUN,XI(JIUN)
      GO TO 1000
 8002 WRITE(6,8009) X0
      GO TO 1000
 8001 WRITE(6,8010) H32
      GO TO 1000
 7001 WRITE(6,7002)
 1000 RETURN
C EINLESE FORMAT-STATEMENTS
   18 FORMAT(6E11.4,I1)
   19 FORMAT(F7.4,E11.4,F7.4,F7.1)
 6000 FORMAT(F4.1,F6.1,F4.1,I3)
   20 FORMAT(3E11.4)
C AUSGABE-FORMAT-STATEMENTS
 5000 FORMAT(1H0,'TURBULENTE GRENZSCHICHT'/)
 5001 FORMAT(1H0,'LAMINARE GRENZSCHICHT'/)
 5002 FORMAT(1H0,'KEINE ABLOESUNG'/)
 5005 FORMAT(1H0,'LAMINARE ABLOESUNG'/)
 5006 FORMAT(1H0,'TURBULENTE ABLOESUNG'/)
 6005 FORMAT(1H0,'UNSTETIGKEIT P0=',E11.4)
 9000 FORMAT(1H0,5X,'UNZULAESSIGER WERT CP'/5X,'CPMAX=',E15.4,
     *2X,'CPMIN=',E15.4/)
 8005 FORMAT(1H0,'FALSCHE REIHENFOLGE',5X,'X(',I4,')=',E12.5,5X,
     *'X(',I4,')=',E12.5/)
 8009 FORMAT(1H0,'UNZULAESSIGER WERT X0=',E12.5/)
 8010 FORMAT(1H0,'UNZULAESSIGER ANFANGSWERT  H32=',E12.5/)
 7002 FORMAT(1H0,'RECHNUNG ABGEBROCHEN WEGEN ZU WEIT GEHENDER SCHRITTVER
     *KLEINERUNG'/)
 1003 FORMAT(1H0,'UMSCHLAGKRITERIUM (UKR) ERLAUBT KEINE TURBULENTE GRENZ
     *SCHICHT'/)
   55 FORMAT(1H ,14X,F9.5,2X,E10.4,2X,E10.4,2X,E10.4,6F10.5,2X,E10.4)
   27 FORMAT(1H ,I3,2X,F9.5,11X,E10.4,2X,E10.4,2X,E10.4,6F10.5,2X,E10.4)
   24 FORMAT(1H0,16X,'BERECHNUNG DER LAMINAREN UND TURBULENTEN GRENZSCHI
     *CHT IN KOMPRESSIBLER STROEMUNG OHNE WAERMEUEBERGANG'/41X,'INTEGRAL
     *SAETZE FUER IMPULS UND KINETISCHE ENERGIE'/////90X,'PROG.NR. 11020
     *019',3X,//10X,'ADIABATENEXPONENT =',F4.1,5X,'ZAEHIGKEIT OMEGA =',
     *F4.1,5X,'THETA=',F7.2,1X,'GRAD',5X,'RUHEDRUCKVERH. P0 =',F7.4//10X
     *,'UNTERPROGRAMM HCFL NRL.',I6,5X,'HCFT NR.',I6/)
 1024 FORMAT(1H0,'ANFANG X0=',F9.5,5X,'OEFFNUNGSWINKEL=',F5.1,1X,
     *'GRAD',5X,'ENDE X=',F9.5,5X,'MZW=',I4//10X,'UMSCLAG LAM./TURBULEN
     *T UKR NR.',I6,5X,'EINGABE  X=',F9.5,5X,'KU=',I2//3X,'N',6X,'XE',
     *8X,'X',9X,'D2',10X,'D3',10X,'CF',9X,'H12',7X,'H12U',6X,'H32',
     *8X,'CP',7X,'MA',9X,'R',9X,'D1'/)
      END
C     TURBULEMTE GRENZSCHICHT
C     UNTERPROGRAMM ZUR GRENZSCHICHTRECHNUNG VOM 8. 12. 67
C     FORMELN NACH DLR FB 65-52
C     REIBUNGSBEIWERT NACH AVA BER. 68A04
C     DISSIPATIONSWERT F. GLEICHGEWICHTSGRENZSCHICHT (NASH) MIT
C     VERZOEGERUNG
C     VERFASSER I. WOELKE
C     SUBROUTINE HCFT NR. 110768 FORTRAN 1. 7. 68 GOETINGEN
C     LAETZTE AENDERUNG JUNI 1971
      SUBROUTINE HCFT(GAMM,D2,D3,REREF,X,RHOU,EMA,TD,R,H12,H32,H12U,CF,
     *CD,NAB,NR)
      COMMON /B/THETA,OMEGA,S,EMAI(200),PA(10),PA1(10)
      DIMENSION XC(20),CDE(20),EL(20)
      NR=010768
      AL=4.
      J=1
      IF(PA(1).GT.0.5)GO TO 16
      PA(1)=0.9
      DO 600 I=1,20
      XC(I)=0.
      CDE(I)=0.
      EL(I)=0.
  600 CONTINUE
      RETURN
   16 NAB=0
      H32=D3/D2
      IF(H32.GT.1.9)  H32=1.9
      B=0.15*TANH(0.22*EMA)**2
      IF(B.GE..0005)GO TO 12
      H32U=H32/(1.+2.*B)
      GO TO 13
   12 H32U=(1.+2.*B)/(2.*B)*(1.-SQRT(1.-(4.*B*H32)/((1.+2.*B)**2)))
   13 IF(H32U.LE.1.34)GO TO 500
      H12U=H32U/(3.*H32U-4.)
      D2UD2=1.+.89*(GAMM-1)/2.*EMA**2*H32*(2.-H32U)*(1.-.043*EMA)
      D=D2*(H12U+1.)*H12U/(H12U-1.)*D2UD2
      DX=D*AL
      TW=(1.+.89*(GAMM-1.)/2.*EMA**2)*TD
      YMW=TW**OMEGA*(1.+S)/(TW+S)
      RE2=REREF*(RHOU*D2*D2UD2)/YMW*TD/TW
      CFO=.246*10.**(-0.678*H12U)*RE2**(-.268)
      CF=CFO*TD/TW
      GI=(H12U-1.)/(H12U*SQRT(CFO/2.))
      CD=CF*H32/2.*(1.+(H12U-1.)/H12U*((GI+1.7)**2/37.21-1.8))
      RER=1.+(DX/D2)*CFO/2.*(1.+(H12U+1.)/H12U*((GI+1.7)**2/37.21-1.8))
      IF(RER.LT.0.1) RER=0.1
C
C     SPEICHERUNG VON CD
   24 IF(PA(1).GT.1.)GO TO 14
      J=1
      XC(1)=X
      CDE(1)=CD
      PA(1)=2.
      GO TO 26
   14 IF(X.LE.X1)  J=J-1
      IF(J.LT.6)GO TO 26
      IF((X+DX).LT.(XC(J)+D))GO TO 77
      IF(J.LT.20)GO TO 26
      DO 18 J=2,20
      XC(J-1)=XC(J)
   18 CDE(J-1)=CDE(J)
      J=20
      GO TO 15
   26 J=J+1
   15 XC(J)=X+DX
      CDE(J)=CD*RER**(-.268)
      X1=X
      GO TO 17
   77 IF(X.LE.X1)   J=J+1
C
C     INTERPOLATION
   17 CD=CDE(1)
      IF(X.LE.XC(1))GO TO 21
      DO 19 JE=1,J
      IF(XC(JE).GT.X)GO TO 20
   19 CONTINUE
      JE=J
   20 M1=JE-2
      IF(M1.LT.1)    M1=1
      CD=0.
      M2=JE+1
      IF(M2.GT.J)   M2=J
      DO 22 JX=M1,M2
      EL(JX)=1.
      DO 23 JY=M1,M2
      IF(JY.EQ.JX)GO TO 23
      EL(JX)=(X-XC(JY))/(XC(JX)-XC(JY))*EL(JX)
   23 CONTINUE
   22 CD=EL(JX)*CDE(JX)+CD
   21 H32=D3/D2
      H12=H12U*D2UD2+(GAMM-1.)/2.*EMA**2*H32
C     KENNZEICHEN FUER ABLOESUNG
      IF(H12U.GE.4.)  NAB=1
      RETURN
  500 H12=0.
      RETURN
      END
C     LAMINARE GRENZSCHICHT
C     UNTERPROGRAMM ZUR GRENZSCHICHTRECHNUNG VOM 8. 12. 67
C     VERFASSER H. SENKBEIL
C     SUBROUTINE HCFL NR. 121267 FORTRAN 30. 9. 67 GOETTINGEN
C     STEWARTSON-ILLINGWORTH-TRANSFORMATION
C     HARTREE-PROFILE, NACH R. EPPLER, ING-ARCH. 32, 1963
      SUBROUTINE HCFL(GAMM,D2,D3,REREF,X,RHOU,EMA,TD,R,H12,H32,H12U,CF,
     *CD,NAB,NR)
      COMMON /B/THETA,OMEGA,S,EMAI(200),PA(10),PA1(10)
      NR=121267
      IF(D2.EQ.0.)RETURN
      NAB=0
      H12U=100000.
      H32=D3/D2
      IF(H32.LE.1.51509)GO TO 501
      IF(H32.LE.1.57258)GO TO 502
      IF(H32.GT.1.8)  H32=1.8
      HSTERN=79.870845-89.582142*H32+25.715786*H32**2
      IF(HSTERN.LT.1.) HSTERN=1.
      EPS=1.372391-4.226253*H32+2.221687*H32**2
      GO TO 503
  502 HSTERN=4.02922-(583.60182-724.55916*H32+227.18220*H32**2)
     **SQRT(H32-1.51509)
      EPS=2.512589-1.686095*HSTERN+0.391541*HSTERN**2-0.031729*HSTERN**3
  503 D=7.853976-10.260551*H32+3.418898*H32**2
      TW=(1.+0.85*(GAMM-1.)/2.*EMA**2)*TD
      YW=TW**OMEGA*(1.+S)/(TW+S)
      RE2=REREF*(RHOU*D2)/YW*TD/TW
      CF=2.*(TD/TW)**2*EPS/RE2
      CD=2./RE2*(TD/TW)**2*D
      H12=TW*HSTERN+(GAMM-1.)/2.*EMA**2
      H32=D3/D2
C     KENNZEICHEN FUER ABLOESUNG
      IF(H32.LE.1.51512)  NAB=1
      RETURN
  501 H12=0.
      RETURN
      END
C     ISENTROPE STROEMUNG
C     UNTERPROGRAMM ZUR GRENZSCHICHTRECHNUNG VOM 8. 12. 1967
C     VERFASSER H. SENKBEIL
C     SUBROUTINE IS  FORTRAN 8. 12. 1967  GOETTINGEN
      SUBROUTINE IS(CP,EMAREF,Q,EMA,RHO,U,T,RHOU,P0,GAMM)
      P=1.+CP*GAMM/2.*EMAREF**2
      RHO=P0**((GAMM-1.)/GAMM)*P**(1./GAMM)
      IF(EMAREF.GE..1)GO TO 10
      Q=1.-CP
      GO TO 11
   10 Q=2./((GAMM-1.)*EMAREF**2)*(RHO*(1.+(GAMM-1.)/2.*EMAREF**2)-P)
   11 EMA=EMAREF*SQRT(Q/P)
      U=SQRT(Q/RHO)
      RHOU=RHO*U
      T=(P/P0)**((GAMM-1.)/GAMM)
      RETURN
      END
C     DISSIPATIONSKRITERIUM FUER DEN UMSCHLAG LAMINAR/TURBULENT
C     UNTERPROGRAMM ZUR GRENZSCHICHTRECHNUNG VOM 8. 12. 1967
C     VERFASSER H. SENKBEIL
C     SUBROUTINE UKR NR 240268 FORTRAN 26. 2. 68  GOETTINGEN
      SUBROUTINE UKR(KU,X,XU,D2,D3,REREF,RHOU,EMA,T,R,H12,CD,NAB,XE,
     *GAMM,NRU)
      NRU=240268
      IF(KU.EQ.1)RETURN
      IF(X.LT.XU.AND.NAB.EQ.0.)RETURN
      KU=1
      CALL HCFT (GAMM,D2,D3,REREF,X,RHOU,EMA,TD,R,H12,H32,H12U,CF,CD,
     *NAB,NR)
      IF(H12.EQ.0..OR.H.EQ.0.)GO TO 2
      IF(CD.GT.CDT)GO TO 1
      XU=X
      XE=X
      RETURN
    1 KU=0
      RETURN
    2 WRITE(6,3)
    3 FORMAT(10X,'DISSIPATIONSKRITERIUM NICHT PRUEFBAR')
      RETURN
      END
C     LAGRANGESCHE INTERPOLATIONSFORMEL
C     UNTERPROGRAMM ZUR GRENZSCHITRECHNUNG VOM 8. 12. 1967
C     VERFASSER  H. SENKBEIL
C     SUBROUTINE INPOL FORTRAN 8. 12. 1967  GOETTINGEN
      SUBROUTINE INPOL(M,CPE,RE,XE,IE,NI,EMAREF,P0,GAMM)
      COMMON /B/THETA,OMEGA,S,EMAI(200),PA(10),PA1(10)
      COMMON /K/XI(200),RI(200)
      DIMENSION EL(200)
      CPE=0.
      RE=0.
      M1=IE-1-NI
      IF(M1.LT.1) M1=1
      M2=IE+NI
      IF(M2.GT.M) M2=M
      DO 2 J=M1,M2
      EL(J)=1.
      DO 3 I=M1,M2
      IF(I.EQ.J)GO TO 3
    1 EL(J)=(XE-XI(I))/(XI(J)-XI(I))*EL(J)
    3 CONTINUE
      CPE=EL(J)*EMAI(J)+CPE
    2 RE=EL(J)*RI(J)+RE
      IF(EMAREF.LT.0.1)  CPE=1.-CPE**2
      IF(EMAREF.GE.0.1)  CPE=2./(GAMM*EMAREF**2)*(P0*((1.+(GAMM-1.)/2.
     *                      *EMAREF**2)/(1.+(GAMM-1.)/2.*EMAREF**2*CPE
     *                        **2))**(GAMM/(GAMM-1.))-1.)
      RETURN
      END
C......................................................................
      SUBROUTINE ENSAIO(NCE,AL,NP,ELREFM,IDELTA)
      COMMON /K/XI(200),RI(200)
      COMMON /F/MC(200),TC(200),PC(200),CPI(200)
      R0=10.**6
      DL=AL/NP
      DO 10 I=NCE+1,NP+NCE
      K=I-NCE
      XI(I)=XI(NCE)+K*DL/ELREFM
      IF(IDELTA.EQ.0)RI(I)=R0
      IF(IDELTA.EQ.1)RI(I)=RI(NCE)
      MC(I)=MC(NCE)
      TC(I)=TC(NCE)
      PC(I)=PC(NCE)
      CPI(I)=CPI(NCE)
      WRITE(6,1027)I,XI(I),MC(I),TC(I),CPI(I),PC(I)
 1027 FORMAT(1H ,I5,5F10.4)
 10   CONTINUE
      NCE=NCE+NP
      RETURN
      END
      SUBROUTINE DIFSUP(YEM,YS,XS,N,MD,G,TREF,MREF,PREF,NCE,ELREFM,NDIF,
     #XS1,YS1,ALF1)
      REAL ML,MR,MN,MREF,MD
      COMMON /K/XC(200),RI(200)
      COMMON /F/YM(200),TC(200),PC(200),CPC(200)
      COMMON /Q/A(150,3),XSG(150),YSG(150),ISG
      COMMON /A/PNF(150,2),X(150,2),Y(150,2),TET(150,2),XM(150,2),
     #NI(150,2)
      DATA I1/1/,Y0/0./,TH0/0./
C     GERACAO DOS PONTOS DO DIFUSOR SUPERSONICO(EQ. CUBICA)
      AA1=3.*(YS-YEM)/(XS**2.)
      AA2=-2.*(YS-YEM)/(XS**3.)
      XSG(1)=0.
      YSG(1)=YEM
      DO 47 I=2,49
      XSG(I)=XSG(I-1)+XS/49.
      YSG(I)=YEM+AA1*(XSG(I)**2.)+AA2*(XSG(I)**3.)
 47   CONTINUE
      ISG=I-1
      AA4=2*(YS-YS1)/(XS1**3.)+TAN(ALF1)/(XS1**2.)
      AA3=0.5*TAN(ALF1)/XS1-1.5*AA4*XS1
      X0=0.
      DO 91 I=ISG+1,49+ISG
      X0=X0+XS1/49
      XSG(I)=XSG(I-1)+XS1/49
      YSG(I)=YS+AA3*(X0**2.)+AA4*(X0**3.)
 91   CONTINUE
      ISG=I-1
C
      DO 49 I=1,ISG
 49   XSG(I)=XSG(I)+XC(NCE)*ELREFM
      CALL SPLNE1
      XINIT=XSG(1)
      XFIM=XSG(ISG)
      ID=0
C
      FG0=G/(G-1)
      FG1=2./(G-1)
      IN=NCE
      XC(IN)=XINIT/ELREFM
      DO 1 I=1,N
      X(I,1)=XINIT
      Y(I,1)=(I-1)*YEM/(N-1)
      TET(I,1)=0.0
      XM(I,1)=MD
1     CONTINUE
      J=2
C
 4    DO 2 I=1,N-1
      XL=X(I,1)
      YL=Y(I,1)
      ML=XM(I,1)
      THL=TET(I,1)
      XR=X(I+1,1)
      YR=Y(I+1,1)
      MR=XM(I+1,1)
      THR=TET(I+1,1)
      CALL CHARA F(XL,YL,THL,ML,XR,YR,THR,MR,XN,YN,THN,MN,NIT,ID,G)
      IF(XN-XFIM)20,21,21
 20   X(I,2)=XN
      Y(I,2)=YN
      TET(I,2)=THN
      XM(I,2)=MN
C     WRITE (6,*)I,J,XN,YN,THN,MN
 2    CONTINUE
      DO 19 I=1,N-1
      X(I,1)=X(I,2)
      Y(I,1)=Y(I,2)
      TET(I,1)=TET(I,2)
      XM(I,1)=XM(I,2)
 19   CONTINUE
C
      J=J+1
      XR=X(1,1)
      YR=Y(1,1)
      THR=TET(1,1)
      MR=XM(1,1)
      CALL CHARA R(XR,YR,THR,MR,XN,MN,NIT,ID,G)
      IF(XN-XFIM)22,21,21
 22   X(1,2)=XN
      Y(1,2)=0.0
      TET(1,2)=0.0
      XM(1,2)=MN
      WRITE(6,*)I1,J,XN,Y0,TH0,MN
      DO 3 I=1,N-2
      XL=X(I,1)
      YL=Y(I,1)
      THL=TET(I,1)
      ML=XM(I,1)
      XR=X(I+1,1)
      YR=Y(I+1,1)
      THR=TET(I+1,1)
      MR=XM(I+1,1)
      CALL CHARA F(XL,YL,THL,ML,XR,YR,THR,MR,XN,YN,THN,MN,NIT,ID,G)
      X(I+1,2)=XN
      Y(I+1,2)=YN
      TET(I+1,2)=THN
      XM(I+1,2)=MN
      WRITE(6,*)I+1,J,XN,YN,THN,MN
 3    CONTINUE
      XL=X(N-1,1)
      YL=Y(N-1,1)
      ML=XM(N-1,1)
      THL=TET(N-1,1)
      CALL CHARA C(XL,YL,THL,ML,XN,YN,THN,MN,G,MD,XINIT,XFIM)
      X(N,2)=XN
      Y(N,2)=YN
      TET(N,2)=THN
      XM(N,2)=MN
      WRITE(6,*)N,J,XN,YN,THN,MN
      XC(IN+1)=XC(IN)+((X(N,2)-X(N,1))/COS((TET(N,2)+TET(N,1))/2))/ELREF
     *M
      YM(IN+1)=XM(N,2)
      CALL TECEPE(FG0,FG1,TREF,MREF,PREF,IN+1,G)
      IN=IN+1
      WRITE(6,1027)IN,XC(IN),YM(IN),TC(IN),CPC(IN),PC(IN)
1027  FORMAT(1H ,I5,5F10.4)
      DO 23 I=1,N
      X(I,1)=X(I,2)
      Y(I,1)=Y(I,2)
      TET(I,1)=TET(I,2)
      XM(I,1)=XM(I,2)
 23   CONTINUE
      J=J+1
      GO TO 4
 21   NDIF=IN-1-NCE
      DO 50 I=NCE,NCE+NDIF
      RI(I)=10.**6
 50   CONTINUE
      RETURN
      END
C     .........................................................
C
C     CALCULA A SOLUCAO DE F(XRAIZ)=0, BISSECCAO
C     .................................................................
      SUBROUTINE RAIZ(F,AX,BX,XRAIZ,ERROR,IRAIZ)
      EXTERNAL F
      X1=AX
      X2=BX
      F1=F(X1)
      IF(F1.NE.0.0E0)GO TO 1
      XRAIZ=X1
      IRAIZ=1
      RETURN
    1 F2=F(X2)
      IF(F1*F2)3,9,3
    2 RETURN
    3 XHOLD=X1
      XOLD=1.0E0
      DO 7 J=1,200
      XRAIZ=(X1*F2-X2*F1)/(F2-F1)
      FRAIZ=F(XRAIZ)
      IF(F1*FRAIZ)4,8,5
    4 X2=XRAIZ
      F2=FRAIZ
      GO TO 6
    5 X1=XRAIZ
      F1=FRAIZ
    6 IF(ABS(1.0E0-XRAIZ/XOLD).LT.ERROR.AND.ABS(FRAIZ).LT.ERROR)
     *GO TO 8
      XOLD=XRAIZ
    7 CONTINUE
      WRITE(6,10)
    8 X1=XHOLD
      IRAIZ=1
      RETURN
    9 XRAIZ=X2
      IRAIZ=2
      RETURN
   10 FORMAT(1X,10('+'),' ITERATION ON XRAIZ FAILS TO CONVERGE IN 200 IT
     *ERATIONS'/1X,11X,'SUBROUTINE RAIZ RETURNS THE 200TH ITERATED VALUE
     * OF XRAIZ')
      END
C......................................................................
      SUBROUTINE SPLNE1
      COMMON /Q/A(150,3),X(150),Y(150),N
      DIMENSION DX(150),D2YDX2(150)
      NM1=N-1
      DO 1 I=1,NM1
      DX(I)=X(I+1)-X(I)
      DYDXI=(Y(I+1)-Y(I))/DX(I)
      IF(I.NE.1)A(I,1)=DYDXI-DYDX
    1 DYDX=DYDXI
      A(2,2)=2.*(DX(1)+DX(2))
      DO 2 I=3,NM1
      A(I,2)=2.*(X(I+1)-X(I-1))-DX(I-1)**2/A(I-1,2)
    2 A(I,1)=A(I,1)-DX(I-1)*A(I-1,1)/A(I-1,2)
      D2YDX2(N)=0.
      D2YDX2(1)=D2YDX2(N)
      D2YDX2(N-1)=A(N-1,1)/A(N-1,2)
      NM2=N-2
      DO 3 I=2,NM2
      J=N-I
    3 D2YDX2(J)=(A(J,1)-DX(J)*D2YDX2(J+1))/A(J,2)
      DO 4 I=1,NM1
      A(I,1)=(Y(I+1)-Y(I))/DX(I)-(D2YDX2(I+1)+2.*D2YDX2(I))*DX(I)
      A(I,2)=3.*D2YDX2(I)
    4 A(I,3)=(D2YDX2(I+1)-D2YDX2(I))/DX(I)
      RETURN
      END
C     .............................................................
C
C     CALCULA A POSICAO (IGT) DE X DENTRO DO VETOR Y, DE IX
C     POSICOES. Y(IGT-1)<=X<=Y(IGT)
C....................................................................
      FUNCTION IGT(IX,X,Y)
      DIMENSION Y(IX)
      IGT=0
      IF(X.LT.Y(1))GO TO 2
      DO 1 I=2,IX
 1    IF(X.GE.Y(I-1).AND.X.LE.Y(I))IGT=I
      IF(IGT.EQ.0)IGT=IX
      GO TO 3
    2 IGT=2
    3 RETURN
      END
C...................................................................
      FUNCTION FCC(X)
      COMMON /P/A,B,THX
      COMMON /Q/C(150,3),XC(150),YC(150),IX
      I=IGT(IX,X,XC)
      ITS=I-1
      DX=X-XC(ITS)
      YITS=YC(ITS)
      IF(X.LT.XC(I))GO TO 1
      YITS=YC(I)
      DX=X-XC(I)
 1    FCC=YITS+DX*(C(ITS,1)+DX*(C(ITS,2)+DX*C(ITS,3)))-(A*X+B)
      Y=YITS+DX*(C(ITS,1)+DX*(C(ITS,2)+DX*C(ITS,3)))
      IF(ITS.EQ.1)GO TO 2
      IF(I.EQ.IX)GO TO 3
      DYDX=((Y-YITS)/DX+(YC(I)-Y)/(XC(I)-X))/2.
      GO TO 4
 2    DYDX=(YC(I)-Y)/(XC(I)-X)
      GO TO 4
 3    IF(DX.NE.0.)GO TO 5
      DX=0.5*(XC(I)-XC(I-1))
      Y=YITS+DX*(C(ITS,1)+DX*(C(ITS,2)+DX*C(ITS,3)))
      DYDX=(YC(I)-Y)/DX
      GO TO 4
 5    DYDX=(Y-YITS)/DX
 4    THX=ATAN(DYDX)
      RETURN
      END
C.....................................................................
      SUBROUTINE CHARA C(XL,YL,THL,ML,XN,YN,THN,MN,G,MD,XI,XFIM)
      REAL ML,MN,MD
      COMMON /R/XNIN,GDD
      COMMON /P/A,B,THX
      EXTERNAL FCC,FPM
      AUX=1./ML
      XMI=ASIN(AUX)
      A=TAN(THL+XMI)
      B=YL-A*XL
      CALL RAIZ(FCC,XI,XFIM,XN,1.E-3,IRAIZ)
      YN=A*XN+B
      WRITE(7,*)XN,YN,THX
      GDD=SQRT((G+1.)/(G-1.))
      XNI=GDD*ATAN(SQRT(ML**2-1.)/GDD)-ATAN(SQRT(ML**2-1.))
      XNIN=XNI+THN-THL
      CALL RAIZ(FPM,1.,MD,MN,1.E-3,IRAIZ)
      THN=THX
      RETURN
      END
C...................................................................
      FUNCTION FPM(X)
      COMMON /R/XNIN,GDD
      FPM=GDD*ATAN(SQRT(X**2-1)/GDD)-ATAN(SQRT(X**2-1))-XNIN
      RETURN
      END
C......................................................................
C......................................................................
C***********   INICIO DAS MODIFICACOES   ******************************
C......................................................................
      SUBROUTINE FGAMM1(GAM,ITRAN)
      COMMON /H/D121,D101,D241,D221,D201,D202,D203,D361,D341,D321,D301,
     *D342,D322,D302,D323,D303,D304
      COMMON /J/A131,A111,A251,A231,A211,A232,A212,A371,A351,A331,A311,
     *A352,A332,A312,A333,A313,A314
      FKAP(A0,B0,C0,D0)=(A0*GAM**2+B0*GAM+C0)/D0
      IF (ITRAN.EQ.1) GO TO 1
      D121=1./2.
      D101=-1./4.
      D241=FKAP(0.,2.,9.,24.)
      D221=-FKAP(0.,4.,15.,24.)
      D201=FKAP(0.,10.,57.,288.)
      D202=-5./8.
      D203=-FKAP(0.,2.,-3.,6.)
      D361=FKAP(556.,1899.,3231.,10308.)
      D341=-FKAP(388.,1233.,1953.,2304.)
      D321=FKAP(304.,858.,1269.,1728.)
      D301=-FKAP(2708.,7839.,14211.,82944.)
      D342=FKAP(52.,99.,375.,384.)
      D322=-FKAP(52.,99.,303.,192.)
      D302=FKAP(200.,72.,639.,1152.)
      D323=-FKAP(0.,5.,-5.,8.)
      D303=FKAP(0.,13.,-27.,48.)
      D304=FKAP(4.,-57.,27.,144.)
      A131=1./4.
      A111=-1./4.
      A251=FKAP(0.,8.,15.,72.)
      A231=-FKAP(0.,20.,45.,96.)
      A211=FKAP(0.,28.,75.,288.)
      A232=FKAP(0.,4.,9.,12.)
      A212=-FKAP(0.,4.,9.,12.)
      A371=FKAP(6836.,16695.,14211.,82944.)
      A351=-FKAP(3380.,8703.,7875.,13824.)
      A331=FKAP(3748.,8859.,8964.,13824.)
      A311=-FKAP(9044.,17631.,20745.,82944.)
      A352=FKAP(556.,1113.,981.,1728.)
      A332=-FKAP(388.,801.,693.,576.)
      A312=FKAP(304.,645.,549.,864.)
      A333=FKAP(52.,3.,-33.,192.)
      A313=-FKAP(52.,27.,-9.,192.)
      A314=-FKAP(0.,1.,1.,4.)
      GO TO 2
1     D121=1./2.
      D101=-1./4.
      D241=FKAP(0.,2.,9.,24.)
      D221=FKAP(0.,4.,3.,-24.)
      D201=FKAP(0.,10.,-15.,288.)
      D202=3./8.
      D203=FKAP(0.,2.,-3.,-6.)
      D361=FKAP(556.,1899.,3231.,10368.)
      D341=FKAP(388.,849.,225.,-2304.)
      D321=FKAP(304.,282.,-27.,1728.)
      D301=FKAP(2708.,2079.,2115.,-82944.)
      D342=FKAP(52.,99.,375.,384.)
      D322=FKAP(52.,99.,-81.,-192.)
      D302=FKAP(92.,180.,351.,1152.)
      D323=FKAP(0.,5.,-5.,-8.)
      D303=FKAP(0.,19.,-21.,-48.)
      D304=FKAP(4.,-57.,27.,144.)
      A131=1./4.
      A111=-1./4.
      A251=FKAP(0.,8.,15.,72.)
      A231=FKAP(0.,20.,9.,-96.)
      A211=FKAP(0.,28.,-33.,288.)
      A232=FKAP(0.,4.,9.,12.)
      A212=FKAP(0.,4.,-9.,-12.)
      A371=FKAP(6836.,16695.,14211.,82944.)
      A351=FKAP(3380.,4863.,675.,-13824.)
      A331=FKAP(3424.,1983.,-756.,13824.)
      A311=FKAP(7100.,-585.,5625.,-82944.)
      A352=FKAP(556.,1113.,981.,1728.)
      A332=FKAP(388.,321.,-387.,-576.)
      A312=FKAP(304.,-75.,549.,864.)
      A333=FKAP(52.,3.,-33.,192.)
      A313=FKAP(52.,27.,-9.,-192.)
      A314=FKAP(0.,1.,1.,-4.)
2     RETURN
      END
C     ..................................................................
C**** SUBROTINA FGAMM3
      SUBROUTINE FGAMM3(GAM)
      COMMON /ZMZM/ZM01,ZM02,ZM03,ZM04,ZM05,ZM06,ZM07,ZM08,ZM09,
     1ZM10,ZM11,ZM12,ZM13,ZM14,ZM15,ZM16,ZM17
      COMMON /QMQM/QM01,QM02,QM03,QM04,QM05,QM06,QM07,QM08,QM09,
     1QM10,QM11,QM12,QM13,QM14,QM15,QM16,QM17
      COMMON /TMTM/TM01,TM02,TM03,TM04,TM05,TM06,TM07,TM08,TM09,
     1TM10,TM11,TM12,TM13,TM14,TM15
      FKAP(A0,B0,C0,D0)=(A0*GAM**2+B0*GAM+C0)/D0
      ZM01=-1./2.
      ZM02=1./4.
      ZM03=FKAP(0.,4.,-3.,48.)
      ZM04=-FKAP(0.,4.,-15.,288.)
      ZM05=-FKAP(0.,2.,3.,6.)
      ZM06=FKAP(0.,4.,-3.,24.)
      ZM07=FKAP(0.,2.,-3.,6.)
      ZM08=FKAP(364.,1089.,909.,20736.)
      ZM09=-FKAP(52.,63.,27.,1536.)
      ZM10=-FKAP(652.,15.,333.,13824.)
      ZM11=FKAP(412.,270.,909.,41472.)
      ZM12=-FKAP(52.,291.,279.,1152.)
      ZM13=FKAP(68.,75.,9.,288.)
      ZM14=-FKAP(460.,303.,225.,6912.)
      ZM15=-FKAP(28.,-3.,-39.,96.)
      ZM16=FKAP(28.,-3.,21.,192.)
      ZM17=FKAP(28.,-39.,45.,144.)
      QM01=1./2.
      QM02=-1./4.
      QM03=FKAP(0.,2.,9.,24.)
      QM04=-FKAP(0.,4.,3.,24.)
      QM05=FKAP(0.,10.,-15.,288.)
      QM06=3./8.
      QM07=-FKAP(0.,2.,-3.,6.)
      QM08=FKAP(556.,1899.,3231.,10368.)
      QM09=-FKAP(388.,849.,225.,2304.)
      QM10=FKAP(304.,282.,-27.,1728.)
      QM11=-FKAP(2708.,2079.,2115.,82944.)
      QM12=FKAP(52.,99.,375.,384.)
      QM13=-FKAP(52.,99.,-81.,192.)
      QM14=FKAP(92.,180.,351.,1152.)
      QM15=-FKAP(0.,5.,-5.,8.)
      QM16=-FKAP(0.,19.,-21.,48.)
      QM17=FKAP(4.,-57.,27.,144.)
      TM01=-1./4.
      TM02=-FKAP(0.,1.,3.,18.)
      TM03=FKAP(0.,4.,-3.,32.)
      TM04=1./4.
      TM05=-FKAP(0.,4.,-3.,24.)
      TM06=FKAP(0.,2.,-3.,6.)
      TM07=FKAP(188.,-915.,-2223.,27648.)
      TM08=FKAP(12.,33.,-43.,384.)
      TM09=-FKAP(844.,-273.,441.,9216.)
      TM10=-FKAP(364.,789.,225.,3456.)
      TM11=FKAP(52.,15.,63.,384.)
      TM12=FKAP(652.,15.,333.,6912.)
      TM13=FKAP(52.,291.,135.,576.)
      TM14=-FKAP(68.,27.,45.,288.)
      TM15=FKAP(28.,-75.,9.,144.)
      RETURN
      END
C     .................................................................
C     CALCULA OS TETAS, SERIE DE KLIEGEL
      SUBROUTINE HALL TM(QQ,Y,Y3,Y5,Y7,TH1,TH2,TH3)
      COMMON /TMTM/TM01,TM02,TM03,TM04,TM05,TM06,TM07,TM08,TM09,
     1TM10,TM11,TM12,TM13,TM14,TM15
      TH1=TM01*Y3+Y*QQ
      TH2=TM02*Y5+TM03*Y3+(TM04*Y3+TM05*Y)*QQ+TM06*Y*QQ*QQ
      TH3=TM07*Y7+TM08*Y5+TM09*Y3+(TM10*Y5+TM11*Y3+TM12*Y)*QQ+
     1(TM13*Y3+TM14*Y)*QQ**2+TM15*Y*QQ**3
      RETURN
      END
C     ................................................................
C
C      CALCULO DOS PARAMETROS DA TUBEIRA DE EMPUXO MAXIMO
C      ...............................................................
      SUBROUTINE PUNKT D(ME,PAB,PTOT,TTOT,R,GAM,YBZG,TRACAO,DMDT,XD,YD,
     .YE,XBP,YBP,IKN)
      INTEGER PNKL
      REAL LBD2,LBD3,ME,MD,L,MKL,MKX,MF
      EXTERNAL FLBD2,FLBD3,FXYD,FCOT,FXYB,FINTD,FINTE,FCHRCE,FTRAC,FXO
      COMMON /YH/YHLF(150)
      COMMON /NI/N,I,XHLF(150)
      COMMON /XYB/XB(150),YB(150),IXYB
      COMMON /FLB/LBD2,LBD3,GAM2,MD,THD
      COMMON /FI/RHOTOT,TOT
      COMMON /FO/AXO,BXO,XMA,XMB
      COMMON /DR/XR,YR,THO,AMO,YO
      COMMON /GM/AG,GAM1,RGAM
      COMMON /FTE/YXE(150),AM(150),THP(150),IMTY
      COMMON /XYD/YYE(150),YED(150),IXYD
      COMMON /A/PNF(150,2),XF(150,2),YF(150,2),THF(150,2),MF(150,2),
     .NI(150,2)
      COMMON /XX/PNKL(150,50),XKL(150,50),YKL(150,50),MKL(150,50),
     .THKL(150,50),NKL(50)
      DIMENSION XD1(150),YD1(150),ICH(150),THX(2),MKX(2),YX(2)
      IKL=IXYB
      BERR=1.E-4
      ISTP=0
      IXYD=0
      IMTY=0
      PI=3.1416
      GAM1=GAM-1
      GAM2=GAM+1
      AG=-GAM/GAM1
      TOT=TTOT
      RGAM=SQRT(R*GAM)
      RHOTOT=PTOT/(R*TOT)
      AUX=1./ME
      ALFAE=ASIN(AUX)
      PE=PTOT*(1.+.5*GAM1*ME**2)**AG
      RHOE=RHOTOT*(1.+.5*GAM1*ME**2)**(-1./GAM1)
      TE=TTOT/(1.+.5*GAM1*ME**2)
      WE=SQRT(TE)*RGAM*ME
      AUX=(PE-PAB)/(.5*RHOE*TAN(ALFAE)*WE**2)
      THE=.5*ASIN(AUX)
      LBD2=SQRT(GAM2/(GAM1+2./ME**2))*COS(THE-ALFAE)/COS(ALFAE)
      LBD3=ME**2*(1.+.5*GAM1*ME**2)**AG*TAN(ALFAE)*SIN(THE)**2
      DTHF=.5*PI/180.
      YED(1)=THE
C**** CALCULO DE YE
      DO 1 I=1,1000
      THD=YED(I)
      CALL RAIZ(FLBD3,1.5,10.,XD1(I),BERR,IER)
      AUX=1./XD1(I)
      ALFA=ASIN(AUX)
      YD1(I)=LBD3/(XD1(I)**2*(1.+.5*GAM1*XD1(I)**2)**AG*TAN(ALFA)
     .*SIN(YED(I))**2)
      YED(I+1)=YED(I)+DTHF
      IF (I.GT.1) THEN
      IF (YD1(I).GE.YD1(I-1)) GO TO 2
      ENDIF
      IMTY=IMTY+1
    1 CONTINUE
    2 DO 3 I=1,IMTY
      J=IMTY+1-I
      YXE(J)=YD1(I)
      AM(J)=XD1(I)
      THP(J)=YED(I)
    3 CONTINUE
C**** CALCULO DE D1
      DO 5 KN=1,IKL
      I=IKL+1-KN
      N=NKL(I)
      DO 4 J=1,N
      XHLF(J)=XKL(J,I)
    4 CONTINUE
      write(*,*)'n=',n,' i=',i,' ikl=',ikl
      IF(MKL(N,I).LT.ME.OR.ISTP.EQ.1)GO TO 6
      CALL RAIZ(FLBD2,XHLF(1),XHLF(N),XDI,BERR,IER)
      AUX=1./MD
      YYEI=LBD3/(MD**2*(1.+.5*GAM1*MD**2)**AG*TAN(ASIN(AUX))*
     .SIN(THD)**2)
      IF(YYEI.GE.1.)GO TO 5
      IXYD=IXYD+1
      ICH(IXYD)=I
      XD1(IXYD)=XDI
      YYE(IXYD)=YYEI
      K=IGT(N,XDI,XHLF)
      YD1(IXYD)=AMED(YKL(K,I),YKL(K-1,I),XHLF(K),XHLF(K-1),XDI)
      XINT=FINTL(FINTD,XB(I),XDI,100)
      YINT=FINTL(FINTE,YYEI,1.,100)
      YED(IXYD)=SQRT(XINT/YINT)-YD1(IXYD)/YYEI
      IF(YED(IXYD)/YED(1).GE.0.)GO TO 5
      ISTP=ISTP+1
    5 CONTINUE
    6 CALL RAIZ(FXYD,YYE(1),YYE(IXYD),YYEI,BERR,IER)
      K=IGT(IXYD,YYEI,YYE)
      XD=AMED(XD1(K),XD1(K-1),YYE(K),YYE(K-1),YYEI)
      YD=AMED(YD1(K),YD1(K-1),YYE(K),YYE(K-1),YYEI)
      YE=YD/YYEI
      L=XD+FINTL(FCOT,YD,YE,100)
      DX=(XD-XHLF(1))/49.
      XF0=XD+DX
C**** CALCULO DA"CARACTERISTICA MEDIA" DIREITA, DE D PARA B
      DO 10 I=1,50
      XF(I,1)=XF0-DX
      DO 8 IK=1,2
      IX=ICH(K+1-IK)
      DO 7 J=1,NKL(IX)
      XHLF(J)=XKL(J,IX)
    7 CONTINUE
      KX=IGT(NKL(IX),XF(I,1),XHLF)
      THX(IK)=AMED(THKL(KX,IX),THKL(KX-1,IX),XHLF(KX),XHLF(KX-1),
     .XF(I,1))
      MKX(IK)=AMED(MKL(KX,IX),MKL(KX-1,IX),XHLF(KX),XHLF(KX-1),XF(I,1))
      YX(IK)=AMED(YKL(KX,IX),YKL(KX-1,IX),XHLF(KX),XHLF(KX-1),XF(I,1))
    8 CONTINUE
      IF(I.EQ.1)CTE=(YX(2)-YD)/(YD-YX(1))
      YF(I,1)=(YX(2)+CTE*YX(1))/(1.+CTE)
      XF0=XF(I,1)
      IF(YX(1).EQ.YX(2))GO TO 9
      THF(I,1)=AMED(THX(1),THX(2),YX(1),YX(2),YF(I,1))
      MF(I,1)=AMED(MKX(1),MKX(2),YX(1),YX(2),YF(I,1))
      GO TO 10
    9 THF(I,1)=.5*(THX(1)+THX(2))
      MF(I,1)=.5*(MKX(1)+MKX(2))
   10 CONTINUE
      DO 11 I=1,50
      J=51-I
      XHLF(J)=XF(I,1)
      YHLF(J)=YF(I,1)
   11 CONTINUE
C**** CALCULO DA POSICAO DE B NO CONTORNO INICIAL
      CALL RAIZ(FXYB,XB(1),XB(IKL),XBP,BERR,IER)
      I=IGT(IKL,XBP,XB)
      YBP=AMED(YB(I),YB(I-1),XB(I),XB(I-1),XBP)
      IKN=IKL-K+1
      DO 12 K=1,50
      J=51-K
      NI(J,1)=0.
      XKL(J,IKN)=XF(K,1)
      YKL(J,IKN)=YF(K,1)
      MKL(J,IKN)=MF(K,1)
      THKL(J,IKN)=THF(K,1)
   12 CONTINUE
      NKL(IKN)=50
      XR=XD
      YR=YD
      YRE=YR/YE
      K=IGT(IMTY,YRE,YXE)
      AMR=AMED(AM(K),AM(K-1),YXE(K),YXE(K-1),YRE)
      THR=AMED(THP(K),THP(K-1),YXE(K),YXE(K-1),YRE)
      WRITE(7,*)IMTY,K,(YXE(I),I=1,IMTY)
      IMTY=IMTY-K+2
      DO 13 I=2,IMTY
      J=K+I-2
      YXE(I)=YXE(J)
      AM(I)=AM(J)
      THP(I)=THP(J)
   13 CONTINUE
      YXE(1)=YRE
      AM(1)=AMR
      THP(1)=THR
      DO 14 I=1,IMTY
      J=IMTY+1-I
      YD1(J)=YXE(I)
      XD1(J)=AM(I)
      YED(J)=THP(I)
   14 CONTINUE
      IKN1=IKN-1
      DO 16 KN=1,IKN1
      I=IKN1+1-KN
      N=NKL(I)
      DO 15 J=1,N
      XHLF(J)=XKL(J,I)
   15 CONTINUE
      CALL RAIZ(FCHRCE,XHLF(1),XHLF(N),XO,BERR,IERR)
      IMTY=IMTY+1
      IF(YO.LT.0.)GO TO 17
      YD1(IMTY)=YO/YE
      XD1(IMTY)=AMO
      YED(IMTY)=THO
      XR=XO
      YR=YO
      WRITE(7,*)I,XO,YO/YE,AMO,THO*57.29,IMTY
   16 CONTINUE
      IMTY=IMTY+1
   17 M=NKL(I-1)
      YD1(IMTY)=0.
      YED(IMTY)=0.
      AXO=XKL(N,I)
      BXO=XKL(M,I-1)
      XMA=MKL(N,I)
      XMB=MKL(M,I-1)
      CALL RAIZ(FXO,AXO,BXO,XO,BERR,IER)
      XD1(IMTY)=AMO
      WRITE(7,*)I,XO,YD1(IMTY),AMO,YED(IMTY),IMTY
      DO 18 K=1,IMTY
      J=IMTY+1-K
      YXE(J)=YD1(K)
      AM(J)=XD1(K)
      THP(J)=YED(K)
   18 CONTINUE
      YXE(IMTY+1)=PAB
      WRITE(7,*)IMTY+1,(YXE(I),I=1,IMTY),YXE(IMTY+1)
      WRITE(7,*)IMTY,(AM(I),I=1,IMTY)
      WRITE(7,*)IMTY,(THP(I)*57.29,I=1,IMTY)
      TRACAO=2.*PI*FINTL(FTRAC,0.,1.,100)*(YBZG*YE/1000.)**2            0
      YXE(IMTY+1)=-1000.
      DMDT=2.*PI*FINTL(FTRAC,0.,1.,100)*(YBZG*YE/1000.)**2              0
      WRITE(7,*)TRACAO,DMDT
      RETURN
      END
C......................................................................
C
C     CALCULO DO MACH NO PONTO D1
C     ............................................................
      FUNCTION FLBD2(X)
      REAL LBD2,LBD3
      INTEGER PNKL
      COMMON /NI/NK,IK,XHLF(150)
      COMMON /FLB/LBD2,LBD3,GAM2,AMX,THX
      COMMON /GM/AG,GAM1,RGAM
      COMMON /XX/PNKL(150,50),XKL(150,50),YKL(150,50),AMKL(150,50),
     .THKL(150,50),NKL(50)
      I=IGT(NK,X,XHLF)
      THX=AMED(THKL(I,IK),THKL(I-1,IK),XHLF(I),XHLF(I-1),X)
      AMX=AMED(AMKL(I,IK),AMKL(I-1,IK),XHLF(I),XHLF(I-1),X)
      AUX=1./AMX
      ALFAX=ASIN(AUX)
      FLBD2=LBD2-SQRT(GAM2/(GAM1+2./AMX**2))*COS(THX-ALFAX)/COS(ALFAX)
      RETURN
      END
C.......................................................................
C
C      INTEGRAL SOBRE A CARACTERISTICA DIREITA, DE B A D
C     .............................................................
      FUNCTION FINTD(X)
      INTEGER PNKL
      COMMON /NI/NK,IK,XHLF(150)
      COMMON /FI/RHOTOT,TTOT
      COMMON /GM/AG,GAM1,RGAM
      COMMON /XX/PNKL(150,50),XKL(150,50),YKL(150,50),AMKL(150,50),
     .THKL(150,50),NKL(50)
      I=IGT(NK,X,XHLF)
      AMX=AMED(AMKL(I,IK),AMKL(I-1,IK),XHLF(I),XHLF(I-1),X)
      YX=AMED(YKL(I,IK),YKL(I-1,IK),XHLF(I),XHLF(I-1),X)
      THX=AMED(THKL(I,IK),THKL(I-1,IK),XHLF(I),XHLF(I-1),X)
      AUX=1./AMX
      ALFAX=ASIN(AUX)
      RHO=RHOTOT*(1.+.5*GAM1*AMX**2)**(-1./GAM1)
      T=TTOT/(1.+.5*GAM1*AMX**2)
      W=SQRT(T)*AMX*RGAM
      FINTD=RHO*W*YX*SIN(ALFAX)/COS(THX-ALFAX)
      RETURN
      END
C.......................................................................
C
C     INTEGRAL SOBRE A CARACTERISTICA ESQUERDA, DE D A E
C      ...............................................................
      FUNCTION FINTE(Y)
      EXTERNAL FLBD3
      COMMON /FTE/YXE(150),AM(150),THF(150),IMTY
      COMMON /FI/RHOTOT,TTOT
      COMMON /GM/AG,GAM1,RGAM
      K=IGT(IMTY,Y,YXE)
      AMX=AMED(AM(K),AM(K-1),YXE(K),YXE(K-1),Y)
      THX=AMED(THF(K),THF(K-1),YXE(K),YXE(K-1),Y)
      AUX=1./AMX
      ALFAX=ASIN(AUX)
      RHO=RHOTOT*(1.+.5*GAM1*AMX**2)**(-1./GAM1)
      T=TTOT/(1.+.5*GAM1*AMX**2)
      FINTE=RHO*AMX*SQRT(T)*RGAM*SIN(ALFAX)*Y/SIN(THX+ALFAX)
      RETURN
      END
C.......................................................................
C
C     CALCULA O MACH NO PONTO D
C     ...............................................................
      FUNCTION FLBD3(X)
      REAL LBD2,LBD3
      COMMON /FLB/LBD2,LBD3,GAM2,AMD,THX
      COMMON /GM/AG,GAM1,RGAM
      AUX=1./X
      ALFAX=ASIN(AUX)
      FLBD3=LBD2-SQRT(GAM2/(GAM1+2./X**2))*COS(THX-ALFAX)/COS(ALFAX)
      RETURN
      END
C.......................................................................
C
C     INTERPOLA PARA ACHAR O YE
C      ...............................................................
      FUNCTION FXYD(X)
      COMMON /XYD/YYE(150),YED(150),IXYD
      K=IGT(IXYD,X,YYE)
      FXYD=AMED(YED(K),YED(K-1),YYE(K),YYE(K-1),X)
      RETURN
      END
C.......................................................................
      FUNCTION FCOT(Y)
      EXTERNAL FLBD3
      COMMON /FTE/YXE(150),AM(150),THF(150),IMTY
      K=IGT(IMTY,Y,YXE)
      AMX=AMED(AM(K),AM(K-1),YXE(K),YXE(K-1),Y)
      AUX=1./AMX
      ALFAX=ASIN(AUX)
      THX=AMED(THF(K),THF(K-1),YXE(K),YXE(K-1),Y)
      FCOT=1./TAN(THX+ALFAX)
      RETURN
      END
C.......................................................................
C
C     ***** INTERPOLACAO LINEAR *****
C     ..............................................................
      FUNCTION AMED(Y2,Y1,X2,X1,X)
      AMED=Y2-(X2-X)*(Y2-Y1)/(X2-X1)
      RETURN
      END
C......................................................................
      FUNCTION SPLINT(XM)
      COMMON /Q/A(150,3),X(150),Y(150),N
      I=IGT(N,XM,X)
      ITS=I-1
      DX=XM-X(ITS)
      YITS=Y(ITS)
      IF(XM.LT.X(I))GO TO 1
      YITS=Y(I)
      DX=XM-X(I)
    1 SPLINT=YITS+DX*(A(ITS,1)+DX*(A(ITS,2)+DX*A(ITS,3)))
      RETURN
      END
C......................................................................
      FUNCTION FXYB(X)
      COMMON /YH/YHLF(150)
      COMMON /NI/NK,IK,XHLF(150)
      COMMON /XYB/XB(150),YB(150),IKL
      COMMON /A/P(150,2),XF(150,2),Y(150,2),T(150,2),Z(150,2),N(150,2)
      I=IGT(IKL,X,XB)
      Y1=AMED(YB(I),YB(I-1),XB(I),XB(I-1),X)
      I=IGT(50,X,XHLF)
      Y2=AMED(YHLF(I),YHLF(I-1),XHLF(I),XHLF(I-1),X)
      FXYB=Y2-Y1
      RETURN
      END
C......................................................................
C
C      INTEGRAL POR SIMPSON (MODIFICADO)
C     ............................................................
      FUNCTION FINTL(F,A,B,N)
      DX=(B-A)/2./N
      X=A
      S=F(X)
      DO 1 I=1,N
      X=X+DX
      S=S+4.*F(X)
      X=X+DX
      FX=F(X)
      S=S+2.*FX
    1 CONTINUE
      FINTL=(S-FX)*DX/3.
      RETURN
      END
C......................................................................
C      TERMINA A CARACTERISTICA A ESQUERDA (QUE PASSA POR D)
C      DENTRO DO LEQUE DE EXPANSAO
C-------------------------------------------------------------------
      FUNCTION FCHRCE(X)
      INTEGER PNKL
      COMMON /NI/NK,IK,XHLF(150)
      COMMON /FTE/YXE(150),AM(150),THF(150),IMTY
      COMMON /XX/PNKL(150,50),XKL(150,50),YKL(150,50),AMKL(150,50),
     .THKL(150,50),NKL(50)
      COMMON /DR/XR,YR,THX,AMX,YX
      K=IGT(NK,X,XHLF)
      THX=AMED(THKL(K,IK),THKL(K-1,IK),XHLF(K),XHLF(K-1),X)
      AMX=AMED(AMKL(K,IK),AMKL(K-1,IK),XHLF(K),XHLF(K-1),X)
      YX=AMED(YKL(K,IK),YKL(K-1,IK),XHLF(K),XHLF(K-1),X)
      AUX=1./AMX
      ALFAX=ASIN(AUX)
      FCHRCE=YR-YX-(XR-X)*TAN(THX+ALFAX)
      RETURN
      END
C......................................................................
C
C     CALCULO DO EMPUXO FORNECIDO PELA TUBEIRA
C     .................................................................
      FUNCTION FTRAC(Y)
      COMMON /FTE/YXE(150),AM(150),THF(150),IMTY
      COMMON /FI/RHOTOT,TTOT
      COMMON /GM/AG,GAM1,RGAM
      K=IGT(IMTY,Y,YXE)
      AMX=AMED(AM(K),AM(K-1),YXE(K),YXE(K-1),Y)
      THX=AMED(THF(K),THF(K-1),YXE(K),YXE(K-1),Y)
      AUX=1./AMX
      ALFAX=ASIN(AUX)
      RHO=RHOTOT*(1.+.5*GAM1*AMX**2)**(-1./GAM1)
      T=TTOT/(1.+.5*GAM1*AMX**2)
      PHI=THX+ALFAX
      W=AMX*SQRT(T)*RGAM
      IF(YXE(IMTY+1).LT.0.)GO TO 1
      P=RHO*T*RGAM**2/(GAM1+1)
      FTRAC=((P-YXE(IMTY+1))+RHO*W**2*SIN(PHI-THX)*COS(THX)/SIN(PHI))*Y
      RETURN
    1 FTRAC=Y*RHO*W*SIN(PHI-THX)/SIN(PHI)
      RETURN
      END
C.....................................................................
      FUNCTION FXO(X)
      COMMON /FO/AXO,BXO,XMA,XMB
      COMMON /DR/XR,YR,THX,AMX,YX
      AMX=AMED(XMA,XMB,AXO,BXO,X)
      AUX=1./AMX
      ALFAX=ASIN(AUX)
      FXO=YR-(XR-X)*TAN(ALFAX)
      RETURN
      END

