       REAL U(100),W(100),Y(100),UDUL(100),VDUL(100),WDUL(100),
     ?UUDUL2(100),VVDUL2(100),WWDUL2(100),UVDUL2(100),UWDUL2(100),
     ?VWDUL2(100),VEL(100),RSSVW(100),RSSUV(100),UV(100),VW(100)
     ?,YDD(100),YY(100),UU(100),VV(100),WW(100),PRES(100)
     ?,XX(100),UTX(16),UTZ(16),UT(16),XXNNUU(16),RHO(16),BBETAW(8)
     ?,AALFA(8)
        COMMON/CEBE/XNU
       REAL APLUS,VKC,DELTA,VE,Z
       CHARACTER*132 TITLE
       CHARACTER*8 POINT
       CHARACTER*50 HEAD
      IC=0
      CONT=0.
       BRAND=1.
       READ(17,*)XUNK1,XUNK2,XUNK3,XUNK4,XUNK5,SKINFC,BETAW,BETAE,XUNK
     ? ,XPRESGR,ZPRESGR
       PI=3.1415926/180.
       BETAW=BETAW*PI
       BETAE=BETAE*PI
       BETAW=BETAW-BETAE
       MMM=1
       READ(MMM,100)TITLE
C      WRITE(6,107)TITLE
       DO 1 I=1,3
       READ(1,100)TITLE
   1   CONTINUE
       READ(1,'(A8,I2)')POINT,N
       READ(1,100)TITLE
       READ(1,'(A39,F7.3)')TITLE,VE
        WRITE(7,*)'VE',VE
       DO 2 I=1,2
C       WRITE(6,*)TITLE
   2   READ(1,100)TITLE
      QSTAR=SQRT(SKINFC/2)*VE
       XNU=1.56E-05
      IC=IC+1
      XXNNUU(IC)=XNU
      BBETAW(IC)=BETAW
      AALFA(IC)=BETAE
      UT(IC)=QSTAR
      RHO(IC)=1.1842
      PRES(IC)=XPRESGR/2.*RHO(IC)*16.63*16.63
      PRES(IC+5)=ZPRESGR/2.*RHO(IC)*16.63*16.63
       DO 3 I=1,N
C      READ(1,99)Y(I),Z,UDUL(I),VDUL(I),WDUL(I),UUDUL2(I),VVDUL2(I),
C    ?WWDUL2(I),UVDUL2(I),UWDUL2(I),VWDUL2(I)
       READ(1,*)Y(I),Z,UDUL(I),VDUL(I),WDUL(I),UUDUL2(I),VVDUL2(I),
     ?WWDUL2(I),UVDUL2(I),UWDUL2(I),VWDUL2(I)
       Y(I)=Y(I)/100.
3      YY(I)=Y(I)
C      WRITE(6,106)POINT,N
       DO 4 I=1,4
       READ(1,100)TITLE
4      CONTINUE
       MMMMM=2
       DO 5 I=1,7
   5   READ(MMMMM,100)TITLE
       DELTA=0.0
       DO 6 I=1,N
       U(I)=UDUL(I)*VE
       W(I)=WDUL(I)*VE
       UU(I)=U(I)
       WW(I)=W(I)
       READ(MMMMM,101)YYY,YDD(I),A,B,C,D,E,F
   6   DELTA=DELTA+YYY/YDD(I)
       DELTA=DELTA/N
C      WRITE(6,105)DELTA
5001   CONTINUE
C       WRITE(6,*)'N',N
       KK=-1
       DO 9 I=1,N-4
53     KK=KK+1
       IF(KK.GT.2)KK=2
       CALL PARFIT(I,N,U,Y,B0,B1,B2)
       DUDY=B1+2*B2*Y(I+KK)
          CALL PARFIT(I,N,W,Y,C0,C1,C2)
          DWDY=C1+2*C2*Y(I+KK)
          IF(DUDY.EQ.0.AND.DWDY.EQ.0.)DUDY=10**50
          IF(DUDY.EQ.0.AND.DWDY.EQ.0.)DWDY=10**50
C      WRITE(6,*)'QSTAR',QSTAR
      XXXX=ATAN(W(I+KK)/U(I+KK))
      Q=SQRT(U(I+KK)**2+W(I+KK)**2)
      S1=Y(I+KK)
      S2=Q*COS(BETAW-XXXX)/QSTAR
      S3=QSTAR*Y(I+KK)/XNU
      S4=U(I+KK)/COS(BETAW)/QSTAR
      S5=U(I+KK)/QSTAR
      S6=U(I+KK)/QSTAR/SQRT(COS(BETAW))
      S7=Y(I+KK)*QSTAR*SQRT(COS(BETAW))/XNU
      S8=Q/QSTAR
      S9=W(I+KK)/QSTAR/SQRT(SIN(BETAW))
      S10=Y(I+KK)*QSTAR*SQRT(SIN(BETAW))/XNU
C     WRITE(7,44)S1,S2,S3,S4,S5,S6,S7,S8,S9,S10
c     write(6,44)BRAND,
c    ?S1,S2,S3,S4,S5,S6,S7,S8,S9,S10
45    format (F2.1,10F7.2)
C44    FORMAT(6F12.5,/,5F12.5)
44    FORMAT(11F12.5)
C      WRITE(8,10)S3,
C    ?ATAN(W(I+KK)/U(I+KK))/PI,ATAN(DWDY/DUDY)/PI,ATAN
C    ?(VWDUL2(I+KK)/UVDUL2(I+KK))/PI
10     FORMAT(4F12.5)
C     **************************
C     FOR VAN DEN BERG'S MODEL
      GO TO 3001
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+5)*SIN(BBETAW(IC)+AALFA(IC))
      DPDZW=PRES(IC+5)*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+KK)*UT(IC)/XXNNUU(IC)
      UPX=(UU(I+KK)*COS(BBETAW(IC))+WW(I+KK)*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+KK)*COS(BBETAW(IC))-UU(I+KK)*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))
3001  CONTINUE
C     **************************
C     **************************
       IF(KK.EQ.2)GO TO 9
       GO TO 53
9      CONTINUE
C      WRITE(7,*)VE
C99    FORMAT(1X,F5.3,1X,F5.2,3F6.2,6F8.4)
100    FORMAT(A)
101    FORMAT(F5.2,7(2X,F8.4))
102    FORMAT(3E13.4)
103    FORMAT(A21,F7.3)
104    FORMAT(3X,'Y cm',9X,'uv/Ul-2',6X,'vw/Ul-2',6X,'uv/Ul-2',6X,
     ?'vw/Ul-2',6X,'meas uv/',5X,'meas vw/',/,42X,'comp',9X,'comp',
     ?9X,'comp uv',6X,'comp vw')
105    FORMAT('BOUNDARY LAYER THICKNESS=',F7.4,'cm')
106    FORMAT(18X,A7,I3)
107    FORMAT(A72)
108    FORMAT(18X,'VE=',F7.3)
200    FORMAT(A7)
109    FORMAT(A50)
110    FORMAT('ANISOTROPY CONSTANT T=',F5.3)
111    FORMAT(6F10.4)
C     **********************************
C      W-L-C METHOD
      GO TO 3002
C     **********************************
      WLCCONT=0.
      CALL BLT (N,UU,WW,YY,BLTHK)
      WRITE(7,*)YY(N),BLTHK
      DO 23 IWLC=1,N
      UTWLC=SQRT(UT(IC)**2*COS(BBETAW(IC)))
      ZETA=TAN(BBETAW(IC))
      UP=UU(IWLC)/UTWLC
      WP=WW(IWLC)/UTWLC
      YP=YY(IWLC)/XXNNUU(IC)*UTWLC
      DPDXS=PRES(IC)*COS(AALFA(IC))
     ? +PRES(IC+5)*SIN(AALFA(IC))
      ALFA=XXNNUU(IC)/RHO(IC)/UTWLC**3*DPDXS
      IF((ALFA*YP).LT.-1)GO TO 123
      S=SQRT(1.+ALFA*YP)
      S0=SQRT(1.+.1108*ALFA)
      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,YP,UP,RHSU,WP,RHSW
556   FORMAT(6(1X,F10.3))
      IF (YY(IWLC)/BLTHK.EQ.1.)WLCCONT=1.
123   CONTINUE
23    CONTINUE
C     **********************************
 3002  CONTINUE
C     **********************************
C     **********************************
C     P-J MODEL
C     GO TO 2008
C     **********************************
      DPDXW=PRES(IC)*COS(BBETAW(IC)+AALFA(IC))
     ?+PRES(IC+5)*SIN(BBETAW(IC)+AALFA(IC))
      DPDZW=PRES(IC+5)*COS(BBETAW(IC)+AALFA(IC))
     ?-PRES(IC)*SIN(BBETAW(IC)+AALFA(IC))
      TETA=(DPDXW*UT(IC))/
     ?(SQRT(DPDXW*DPDXW+DPDZW*DPDZW)*SQRT(UT(IC)*UT(IC)
     ?))
C     PRINT *,'TETA',ACOS(TETA)/PI
      ALFA=SQRT(DPDXW*DPDXW+DPDZW*DPDZW)/RHO(IC)
C     PRINT *,'ALFA',ALFA
      SUM=0.
      SUM1=0.
      YPRE=1.
C     RHSPRE=1/.41*(1/1.*(1-2*TETA*(ALFA*XXNNUU(IC)/
C    ?(UT(IC))**3)*1.+
C    ?(ALFA*XXNNUU(IC)/
C    ?(UT(IC))**3)**2*1.**2)**(1./4.))
      DO 49 IPJ=1,N
      YP=YY(IPJ)*UT(IC)/XXNNUU(IC)
C     RHS=1/.41*(1/YP*(1-2*TETA*(ALFA*XXNNUU(IC)/
C    ?(UT(IC))**3)*YP+
C    ?(ALFA*XXNNUU(IC)/
C    ?(UT(IC))**3)**2*YP**2)**(1./4.))
      IF(IPJ.EQ.1)RL=SQRT(WW(1)*WW(1)+UU(1)*UU(1))/(UT(IC))
      IF(IPJ.NE.1)RL=SQRT((WW(IPJ)-WW(IPJ-1))**2+(UU(IPJ)-UU(IPJ-1))**2)
     ?/(UT(IC))
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)/
     ?(UT(IC))**3)*YP+
     ?(ALFA*XXNNUU(IC)/
     ?(UT(IC))**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     **********************************
2008  CONTINUE
C     **********************************
      GO TO 3000
C     **********************************
C     EAST-HOXEY MODEL
C     **********************************
      UBL=SQRT(UU(N)**2+WW(N)**2)
      U0=16.63
      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
      BETA0=ASIN(1/(XK*UT(IC)/U0)*SIN(GAMA))-GAMA
      WRITE(7,*)BETA0/PI,BETAW/PI
      DO 64 IEH=1,N
      XLHS=UU(IEH)/COS(BETA0)/UT(IC)
      YP=YY(IEH)*UT(IC)/XXNNUU(IC)
      WRITE(6,558)BRAND-1,YP,XLHS
558   FORMAT(3(1X,F10.4))
64    CONTINUE
C     **************************
3000  CONTINUE
C     **************************
C     **************************
       STOP
       END
