      REAL YY(100),UU(100),VV(100),WW(100),CC1(100),CC2(100),CC3(100)
     ?,CC2C3(100),CC1C3(100),CC1C2(100),S3(100),Q(100),PRES(100)
     ?,XX(100),UTX(16),UTZ(16),UT(16),XXNNUU(16),RHO(16),BBETAW(8)
     ?,UUS2(100),VVS2(100),WWS2(100)
     ?,VVWS(100),UUWS(100),UUVS(100),AALFA(100)
      IC=0
       DO 36 JCONT=1,12
 36    READ(4,*)XX(JCONT),PRES(JCONT)
*      DO 37 JCONT1=1,8
*37    READ(5,*)XX(JCONT1),UTX(JCONT1),UTZ(JCONT1)
      CONT=0.
      BRAND=0.
      I=0
100   READ(1,*,END=300)Y,U,V,W,US2,VS2,WS2,VWS,UWS,UVS
      IF(Y.NE.0.)GO TO 121
      BRAND=BRAND+1
      IC=IC+1
      IF(Y.EQ.0.)READ(1,3)NST,XDL,ZDL,QEDUREF,ALFA,CF,BETAW
      IF(Y.EQ.0.)WRITE(5,*)'NST,XDL,ZDL,QEDUREF,ALFA,CF,BETAW',
     ? NST,BRAND
      PI=3.1415926/180.
      XNU=1.551245E-05
      ALFA=ALFA*PI
      BETAW=BETAW*PI
      UREF=30.25
      QE=QEDUREF*UREF
      QSTAR=SQRT(CF/2)*QE
      XXNNUU(IC)=XNU
      BBETAW(IC)=BETAW
      AALFA(IC)=ALFA
      UT(IC)=QSTAR
      RHO(IC)=1.1842
      PRES(IC+6)=PRES(IC+6)/2.*RHO(IC)*30.25*30.25
      PRES(IC)=PRES(IC)/2.*RHO(IC)*30.25*30.25
      WRITE(5,*)'QSTAR,',QSTAR,QSTAR/UREF
      IF(Y.EQ.0..AND.BRAND.NE.1.) GO TO 200
      IF(Y.EQ.0.)GO TO 100
3     FORMAT(11X,I7,4X,F5.4,9X,F5.4,13X,F5.4,11X,F4.1,8X,F7.6,11X,F4.1
     ?)
121   I=I+1
      Q(I)=SQRT(U*U+W*W)*QE
      YY(I)=Y
      UU(I)=U*QE
      VV(I)=V*QE
      WW(I)=W*QE
      UUS2(I)=US2
      VVS2(I)=VS2
      WWS2(I)=WS2
      VVWS(I)=VWS
      UUWS(I)=UWS
      UUVS(I)=UVS
      XXXX=ATAN(W/U)
      S1=Y
      S2=Q(I)*COS(BETAW-XXXX)/QSTAR
      S3(I)=QSTAR*Y/XNU
      S4=U*QE/COS(BETAW)/QSTAR
      S5=U*QE/QSTAR
      S6=U*QE/QSTAR/SQRT(COS(BETAW))
      S7=Y*QSTAR*SQRT(COS(BETAW))/XNU
      S8=Q(I)/QSTAR
      IF(BETAW.EQ.0.) GO TO 111
      S9=W*QE/QSTAR/SQRT(SIN(BETAW))
      S10=Y*QSTAR*SQRT(SIN(BETAW))/XNU
C     **************************
C     FOR VAN DEN BERG'S MODEL
       GO TO 111
C     **************************
      BETAX=XXNNUU(IC)/UT(IC)**2*UTX(IC)
      BETAZ=XXNNUU(IC)/UT(IC)**2*UTZ(IC)
      BETAX=0.
      BETAZ=0.
      DPDXW=PRES(IC)*COS(BBETAW(IC)+AALFA(IC))
     ?+PRES(IC+6)*SIN(BBETAW(IC)+AALFA(IC))
      DPDZW=PRES(IC+6)*COS(BBETAW(IC)+AALFA(IC))
     ?-PRES(IC)*SIN(BBETAW(IC)+AALFA(IC))
      ALFAX=XXNNUU(IC)/RHO(IC)/UT(IC)**3*DPDXW
      ALFAZ=XXNNUU(IC)/RHO(IC)/UT(IC)**3*DPDZW
      YP=YY(I)*UT(IC)/XXNNUU(IC)
      UPX=(UU(I)*COS(BBETAW(IC))+WW(I)*SIN(BBETAW(IC)))/UT(IC)
      RHSX=1/.41*(LOG(YP)+1/2.*ALFAX*YP
     ?+1/2.*BETAX*(LOG(YP))**2/(.41)**2*YP)
      UPZ=(WW(I)*COS(BBETAW(IC))-UU(I)*SIN(BBETAW(IC)))/UT(IC)
      RHSZ=1/.41*(ALFAZ*YP+BETAZ*(LOG(YP))**2/(.41)**2*YP)
      WRITE(6,555)BRAND,YP,UPX,RHSX,UPZ,RHSZ
555   FORMAT(6(1X,F10.3))
C     **************************
C     **************************
111   CONTINUE
C      WRITE(6,4)BRAND,S3(I),S4
C     WRITE(6,4)BRAND,S1,S2,S3(I),S4,S5,S6,S7,S8,S9,S10
C4     FORMAT(6F12.5,/,5F12.5)
4     FORMAT(11F12.5)
      GO TO 100
300   CONT=1
      IC=IC+1
      BRAND=BRAND+1
200   CONTINUE
       N=I
       WRITE(5,*)N
       KK=-1
       DO 9 I=1,N-4
53     KK=KK+1
       IF(KK.GT.2)KK=2
         CALL PARFIT(I,N,UU,YY,B0,B1,B2)
         DUDY=B1+2*B2*YY(I+KK)
           CALL PARFIT(I,N,WW,YY,C0,C1,C2)
           DWDY=C1+2*C2*YY(I+KK)
              IF(DUDY.EQ.0.AND.DWDY.EQ.0.)DUDY=10**50
              IF(DUDY.EQ.0.AND.DWDY.EQ.0.)DWDY=10**50
                 CALL PARFIT(I,N,Q,YY,D0,D1,D2)
                 DQDY=D1+2*D2*YY(I+KK)
          TEMP=(Q(I+KK)*DQDY-UU(I+KK)*DUDY)/WW(I+KK)
C     WRITE(5,*)'TEMP/DWDY',TEMP/DWDY
      IF(UUVS(I).NE.0.)GO TO 11
      IF(UUVS(I).EQ.0.)GO TO 12
C       WRITE(5,*)YY(I),((WW(I)-C0+C1*YY(I)
C     ?+C2*YY(I)**2)**2+(WW(I+1)-
C     ?C0+C1*YY(I+1)+C2*YY(I+1)**2)**2+(WW(I+2)-C0+C1*YY(I+2)+C2*YY(I+2)
C     ?**2)**2+(WW(I+3)-C0+C1*YY(I+3)+C2*YY(I+3)**2)**2+
C     ?(WW(I+4)-C0+C1*YY(I+4)+C2*YY(I+4)**2)**2)
C     ?/(WW(I)**2+WW(I+1)**2+WW(I+2)**2+WW(I+3)**2+WW(I+4)**2)
C     ?,((UU(I)-B0+B1*YY(I)
C     ?+B2*YY(I)**2)**2+(UU(I+1)-
C     ?B0+B1*YY(I+1)+B2*YY(I+1)**2)**2+(UU(I+2)-B0+B1*YY(I+2)+B2*YY(I+2)
C     ?**2)**2+(UU(I+3)-B0+B1*YY(I+3)+B2*YY(I+3)**2)**2+
C     ?(UU(I+4)-B0+B1*YY(I+4)+B2*YY(I+4)**2)**2)
C     ?/(UU(I)**2+UU(I+1)**2+UU(I+2)**2+UU(I+3)**2+UU(I+4)**2)
10     FORMAT(5F12.5)
11     CONTINUE
C      WRITE(2,10)BRAND-1,S3(I+KK),
C    ?ATAN(WW(I+KK)/UU(I+KK))/PI,ATAN(DWDY/DUDY)/PI,ATAN
C    ?(VVWS(I+KK)/UUVS(I+KK))/PI
      GO TO 13
12     CONTINUE
       XXX=9999.
       WRITE(2,10)BRAND-1,S3(I+KK),
     ?ATAN(WW(I+KK)/UU(I+KK))/PI,ATAN(DWDY/DUDY)/PI,XXX
13     IF(KK.EQ.2)GO TO 9
       GO TO 53
9      CONTINUE
       N=I
C     **********************************
C      W-L-C METHOD
C     **********************************
      GO TO 3003
C     PRINT *,BBETAW(IC-1),UT(IC-1),RHO(IC-1),PRES(IC-1)
      CALL BLT (N,UU,WW,YY,BLTHK)
      WLCCONT=0.
      DO 23 IWLC=1,N
      UTWLC=SQRT(UT(IC-1)**2*COS(BBETAW(IC-1)))
      ZETA=TAN(BBETAW(IC-1))
      UP=UU(IWLC)/UTWLC
      WP=WW(IWLC)/UTWLC
      YP=YY(IWLC)/XXNNUU(IC-1)*UTWLC
      DPDXS=PRES(IC-1)*COS(AALFA(IC-1))
     ?+PRES(IC-1+6)*SIN(AALFA(IC-1))
      ALFA=XXNNUU(IC-1)/RHO(IC-1)/UTWLC**3*DPDXS
C     PRINT *,'ALFA',ALFA,YP
      IF((ALFA*YP).LT.-1)GO TO 123
      S=SQRT(1.+ALFA*YP)
      S0=SQRT(1.+.1108*ALFA)
C     PRINT *,(S-1)/(S+1)*(S0+1)/(S0-1)
      IF (((S-1)/(S+1)*(S0+1)/(S0-1)).LE.0.)GO TO 123
      RHSU=1/.41*(2*(S-S0)+LOG((S-1)/(S+1)*(S0+1)/(S0-1)))
      RHSW=RHSU*ZETA*(1-YY(IWLC)/BLTHK)**2
      IF(WLCCONT.EQ.1.)GO TO 123
      WRITE(6,556)BRAND-1,YP,UP,RHSU,WP,RHSW
556   FORMAT(6(1X,F10.3))
      IF(YY(IWLC)/BLTHK.EQ.1.)WLCCONT=1.
123   CONTINUE
23    CONTINUE
3003  CONTINUE
C     **********************************
C     **********************************
C     **********************************
C     P-J MODEL
C     GO TO 3001
C     **********************************
      DPDXW=PRES(IC-1)*COS(BBETAW(IC-1)+AALFA(IC-1))
     ?+PRES(IC-1+6)*SIN(BBETAW(IC-1)+AALFA(IC-1))
      DPDZW=PRES(IC-1+6)*COS(BBETAW(IC-1)+AALFA(IC-1))
     ?-PRES(IC-1)*SIN(BBETAW(IC-1)+AALFA(IC-1))
      TETA=(DPDXW*UT(IC-1))/
     ?(SQRT(DPDXW*DPDXW+DPDZW*DPDZW)*SQRT(UT(IC-1)*UT(IC-1)
     ?))
C     PRINT *,'TETA',ACOS(TETA)/PI
      ALFA=SQRT(DPDXW*DPDXW+DPDZW*DPDZW)/RHO(IC-1)
C     PRINT *,'ALFA',ALFA
      SUM=0.
      SUM1=0.
      YPRE=1.
C     RHSPRE=1/.41*(1/1.*(1-2*TETA*(ALFA*XXNNUU(IC-1)/
C    ?(UT(IC-1))**3)*1.+
C    ?(ALFA*XXNNUU(IC-1)/
C    ?(UT(IC-1))**3)**2*1.**2)**(1./4.))
      DO 49 IPJ=1,N
      YP=YY(IPJ)*UT(IC-1)/XXNNUU(IC-1)
C     RHS=1/.41*(1/YP*(1-2*TETA*(ALFA*XXNNUU(IC-1)/
C    ?(UT(IC-1))**3)*YP+
C    ?(ALFA*XXNNUU(IC-1)/
C    ?(UT(IC-1))**3)**2*YP**2)**(1./4.))
      IF(IPJ.EQ.1)RL=SQRT(WW(1)*WW(1)+UU(1)*UU(1))/(UT(IC-1))
      IF(IPJ.NE.1)RL=SQRT((WW(IPJ)-WW(IPJ-1))**2+(UU(IPJ)-UU(IPJ-1))**2)
     ?/(UT(IC-1))
C     PRINT *,RL
      SUM1=SUM1+RL
C     SUM=SUM+(RHS+RHSPRE)/2.*(YP-YPRE)
C     YPRE=YP
C     RHSPRE=RHS
       WRITE(8,557)BRAND-1,YP,SUM1
557   FORMAT(3(1X,F10.4))
49    CONTINUE
C     **********************************
C     **********************************
      SUM=0.
      YPRE=0.
      RHSPRE=0.
      DO 87 IPJ1=1,4
      XPJ3=IPJ1-1.
      DO 88 IPJ2=100,900
      XPJ4=IPJ2/100.
      YP=10**(XPJ3)*XPJ4
      IF (IPJ1.EQ.1.AND.IPJ2.EQ.100) RHS=0
      IF (IPJ1.EQ.1.AND.IPJ2.EQ.100) GO TO 92
      IF (YP.GT.4000)GO TO 88
      RHS=1/.41*(1/YP*(1-2*TETA*(ALFA*XXNNUU(IC-1)/
     ?(UT(IC-1))**3)*YP+
     ?(ALFA*XXNNUU(IC-1)/
     ?(UT(IC-1))**3)**2*YP**2)**(1./4.))
      SUM=SUM+(RHS+RHSPRE)/2.*(YP-YPRE)
92    YPRE=YP
      RHSPRE=RHS
      IF (IPJ2.EQ.100.OR.IPJ2.EQ.200.OR.IPJ2.EQ.300.
     ?OR.IPJ2.EQ.400.OR.IPJ2.EQ.500.OR.IPJ2.EQ.600.
     ?OR.IPJ2.EQ.700.OR.IPJ2.EQ.800.OR.IPJ2.EQ.900)
     ? WRITE(6,561)BRAND-1,YP,SUM+5.61
561   FORMAT(4(1X,F10.4))
88    CONTINUE
87    CONTINUE
C     **********************************
C     **********************************
3001  CONTINUE
C     **********************************
C     **********************************
      GO TO 3000
C     EAST-HOXEY MODEL
C     **********************************
      UBL=SQRT(UU(N)**2+WW(N)**2)
      IF((IC-1).EQ.1)U0=UBL
      XK=19.45
      YPRE=0.
      DELPRE1=0.
      DELPRE2=0.
      SUM1=0.
      SUM2=0.
      DO 63 IEH=1,N
      DEL1=(1-UU(IEH)/UBL)
      DEL2=-WW(IEH)/UBL
      Y1=YY(IEH)
      SUM1=SUM1+(DEL1+DELPRE1)/2.*(Y1-YPRE)
      SUM2=SUM2+(DEL2+DELPRE2)/2.*(Y1-YPRE)
      DELPRE1=DEL1
      DELPRE2=DEL2
63    YPRE=Y1
      GAMA=-SUM2/SUM1
      ZARI=(1/(XK*UT(IC-1)/U0)*SIN(GAMA))
      IF(ZARI.GT.1.)WRITE(7,*)ZARI
      IF(ZARI.GT.1.)GO TO 301
      BETA0=ASIN(1/(XK*UT(IC-1)/U0)*SIN(GAMA))-GAMA
      DO 64 IEH=1,N
      XLHS=UU(IEH)/COS(BETA0)/UT(IC-1)
      YP=YY(IEH)*UT(IC-1)/XXNNUU(IC-1)
      WRITE(6,558)BRAND-1,YP,XLHS
558   FORMAT(3(1X,F10.4))
64    CONTINUE
C     **************************
3000  CONTINUE
C     **************************
C     **************************
      I=0
      IF(CONT.NE.0.)GO TO 301
      GO TO 100
301   STOP
      END
