      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)
      CHARACTER*80 AAA
      CHARACTER*80 BBB
      BRAND=0.
      I=0
      IC=0
      DO 37 JCONT1=1,4
37    READ(5,*)XX(JCONT1),UTX(JCONT1),UTZ(JCONT1)
      CONT=0.
100   READ(1,*,END=300)Y,U,W
      IF(Y.NE.0.)GO TO 121
      BRAND=BRAND+1
      IF(Y.EQ.0.)READ(1,1001)
1001  FORMAT(/)
      IC=IC+1
      READ(1,*)BETAW,UTAU,XLOC
C     WRITE(7,*)BETAW,UTAU,XLOC
      PI=3.1415926/180.
      XNU=1.555688E-05
C     WRITE(3,*)UTAU,BETAW,XLOC*2.54/100.
      BETAW=BETAW*PI
      XXNNUU(IC)=XNU
      BBETAW(IC)=BETAW
      UT(IC)=UTAU
C      WRITE(5,*)XNU,BETAW,UTAU,COS(BETAW)
      IF(Y.EQ.0..AND.BRAND.NE.1.) GO TO 200
      IF(Y.EQ.0.)GO TO 100
121   I=I+1
      Q(I)=SQRT(U*U+W*W)
      YY(I)=Y
      UU(I)=U
      WW(I)=W
      S1=Y
      S2=Q(I)*COS(BETAW-ATAN(W/U))/UTAU
      S3(I)=UTAU*Y/XNU
      S4=U/COS(BETAW)/UTAU
      S5=U/UTAU
      S6=U/UTAU/SQRT(COS(BETAW))
      S7=Y*UTAU*SQRT(COS(BETAW))/XNU
      S8=Q(I)/UTAU
      IF(BETAW.EQ.0.) GO TO 111
      S9=W/UTAU/SQRT(SIN(BETAW))
      S10=Y*UTAU*SQRT(SIN(BETAW))/XNU
111   CONTINUE
C     **************************
C     FOR VAN DEN BERG'S MODEL
      GO TO 5001
C     **************************
      BETAX=XXNNUU(IC)/UT(IC)**2*(UTX(IC)+TAN(BBETAW(IC))*UTZ(IC))
      BETAZ=XXNNUU(IC)/UT(IC)**2*(UTZ(IC)-TAN(BBETAW(IC))*UTX(IC))
      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.*BETAX*(LOG(YP))**2/(.41)**2*YP)
      UPZ=(WW(I)*COS(BBETAW(IC))-UU(I)*SIN(BBETAW(IC)))/UT(IC)
      RHSZ=1/.41*(BETAZ*(LOG(YP))**2/(.41)**2*YP)
      WRITE(6,230)BRAND,YP,UPX,RHSX,UPZ,RHSZ
230   FORMAT(6(1X,F11.3))
C     **************************
5001  CONTINUE
C     **************************
C      WRITE(6,4)BRAND,S1,S2,S3(I),S4,S5,S6,S7,S8,S9,S10
4     FORMAT(11F12.5)
C       WRITE(6,4)BRAND,S1,S2,S3(I),S4,S5,S6,S7,S8,S9,S10
C4     FORMAT(6F12.5,/,5F12.5)
      GO TO 100
300   CONT=1
      BRAND=BRAND+1
      IC=IC+1
200   CONTINUE
       N=I
C       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)-WW(I+KK)/(.0635+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
10     FORMAT(5F12.5)
C      WRITE(2,10)BRAND-1,S3(I+KK),
C    ?ATAN(WW(I+KK)/UU(I+KK))/PI,ATAN(DWDY/DUDY)/PI
       IF(KK.EQ.2)GO TO 9
       GO TO 53
9      CONTINUE
C     **********************************
C     **********************************
C     P-J MODEL
C     GO TO 3009
C     **********************************
      SUM=0.
      SUM1=0.
      YPRE=1.
C     RHSPRE=1/.41
      DO 49 IPJ=1,N
      YP=YY(IPJ)*UT(IC-1)/XXNNUU(IC-1)
C     RHS=1/.41*(1/YP)
      IF(IPJ.EQ.1)
     ? RL=SQRT((WW(1)-18.1356)*(WW(1)-18.1356)+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))
      SUM1=SUM1+RL
C     SUM=SUM+(RHS+RHSPRE)/2.*(YP-YPRE)
C     YPRE=YP
C     RHSPRE=RHS
      WRITE(8,*)BRAND-1,YP,SUM1
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     **********************************
3009  CONTINUE
C     **********************************
C     EAST-HOXEY MODEL
      GO TO 3003
C     **********************************
      UBL=SQRT(UU(N)**2+WW(N)**2)
      IF((IC-1).EQ.1)U0=UBL
      XK=19.45
      YPRE=0.
      DELPRE1=0.
      DELPRE2=-18.1356/UBL
      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-1)/U0)*SIN(GAMA))-GAMA
      DO 64 IEH=1,N
      XLHS=UU(IEH)/COS(BETA0)/UT(IC-1)
      RHS=YY(IEH)*UT(IC-1)/XXNNUU(IC-1)
      YP=RHS
      WRITE(6,*)BRAND-1,YP,XLHS
64    CONTINUE
C     **************************
3003  CONTINUE
C     **************************
C     **********************************
C      W-L-C METHOD
      GO TO 3002
C     **********************************
C     PRINT *,BBETAW(IC),UT(IC),RHO(IC),PRES(IC)
      WLCCONT=0.
      CALL BLT (N,UU,WW,YY,BLTHK)
      UTWLC=SQRT(UT(IC-1)**2*COS(BBETAW(IC-1)))
      ZETA=TAN(BBETAW(IC-1))
      WRITE(7,*)ZETA
      DO 23 IWLC=1,N
      UP=UU(IWLC)/UTWLC
      WP=WW(IWLC)/UTWLC
      YP=YY(IWLC)/XXNNUU(IC-1)*UTWLC
      RHSU=1/.41*LOG(YP)+5.5
      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
C     **********************************
 3002  CONTINUE
C     **************************
      I=0
      IF(CONT.NE.0.)GO TO 301
      GO TO 100
301   STOP
      END
