      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 AAAA*71,BBBB*33,CCCC*68
      S9=0.0
      S10=0.0
      BRAND=0.
      I=0
      IC=0
      DO 36 JCONT=1,16
36    READ(4,*)XX(JCONT),PRES(JCONT)
      DO 37 JCONT1=1,8
37    READ(5,*)XX(JCONT1),UTX(JCONT1),UTZ(JCONT1)
      CONT=0.
100   READ(1,*,END=300)Y,U,BETA,US,VS,WS,UV,VW,UW
C     WRITE(7,*)Y,U,BETA,US,VS,WS,UV,VW,UW
      IF(Y.NE.0.)GO TO 121
      BRAND=BRAND+1.
C      WRITE(6,*)BRAND
      READ(1,88)BBBB
      WRITE(7,88)BBBB
      WRITE(7,*)
      READ(1,88)CCCC
      WRITE(7,88)CCCC
      WRITE(7,*)
88    FORMAT(A)
      IC=IC+1
      READ(1,*)NST,X,UD,Z,UTAU,XNU,BETAW,RHO(IC)
      WRITE (7,*)NST,X,UD,Z,UTAU,XNU,BETAW,RHO(IC)
      READ(1,88)AAAA
      WRITE(7,88)AAAA
      WRITE(7,*)
      PI=3.1415926/180.
      ZERO=0.
C      WRITE(3,*)X,UTAU,BETAW,ZERO
      BETAW=BETAW*PI
      XXNNUU(IC)=XNU
      BBETAW(IC)=BETAW
      UT(IC)=UTAU
C     WRITE(6,*)'UT,RHO,XNU',UT(IC),RHO(IC),XXNNUU(IC)
      IF(Y.EQ.0..AND.BRAND.NE.1.) GO TO 200
      IF(Y.EQ.0.)GO TO 100
121   I=I+1
      IF(BETA.EQ.999.)I=I-1
      IF(BETA.EQ.999.)GO TO 500
      BETA=BETA*PI/1.111111111
      Q(I)=U
      W=U*SIN(BETA)
      U=U*COS(BETA)
      YY(I)=Y
      UU(I)=U
      WW(I)=W
      CC1(I)=US
      CC2(I)=VS
      CC3(I)=WS
      CC1C2(I)=UV
      CC1C3(I)=UW
      CC2C3(I)=VW
      S1=Y
      S2=Q(I)*COS(BETAW-BETA)/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
c       IF(BETAW.EQ.0.) GO TO 111
       IF(BETAW.EQ.0.)S10=999.
       IF(BETAW.EQ.0.)S9=999.
       if(betaw.eq.0)go to 234
       S9=-W/UTAU/SQRT(SIN(-BETAW))
       S10=Y*UTAU*SQRT(SIN(-BETAW))/XNU
234   CONTINUE
C     WRITE(6,4)BRAND,S1,S2,S3(I),S4,S5,S6,S7,S8,S9,S10
4     FORMAT(11F12.5)
C      WRITE(2,4)BRAND,S1,S2,S3(I),S4,S5,S6,S7,S8,S9,S10
c4     FORMAT(6F12.5,/,5F12.5)
C     **************************
C     FOR VAN DEN BERG'S MODEL
        GO TO 111
C     **************************
C     BETAX=XXNNUU(IC)/UT(IC)**2*(UTX(IC)+TAN(BBETAW(IC))*UTZ(IC))
C     BETAZ=XXNNUU(IC)/UT(IC)**2*(UTZ(IC)-TAN(BBETAW(IC))*UTX(IC))
      BETAX=0.0
      BETAZ=0.0
      DPDXW=PRES(IC)*COS(BBETAW(IC))+PRES(IC+8)*SIN(BBETAW(IC))
      DPDZW=PRES(IC+8)*COS(BBETAW(IC))-PRES(IC)*SIN(BBETAW(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
500   GO TO 100
300   CONT=1
      IC=IC+1
      BRAND=BRAND+1
200   CONTINUE
       N=I
C      WRITE(6,*)N
       KK=-1
       DO 9 I=1,N-4
C       WRITE(6,*)'I AM GOING N'
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
      IF(CC1C2(I+KK).NE.999.AND.CC2C3(I+KK).NE.999.)GO TO 11
      IF(CC1C2(I+KK).EQ.999.OR.CC2C3(I+KK).EQ.999.)GO TO 12
10     FORMAT(5F12.5)
11     CONTINUE
C       IF(CC1C2(I+KK)/CC2C3(I+KK).LE.0.) WRITE(2,10)BRAND-1,S3(I+KK),
C     ?ATAN(WW(I+KK)/UU(I+KK))/PI,ATAN(DWDY/DUDY)/PI,ATAN
C     ?(CC2C3(I+KK)/CC1C2(I+KK))/PI
C      IF(CC1C2(I+KK)/CC2C3(I+KK).GT.0.) WRITE(2,10)BRAND-1,S3(I+KK),
C     ?ATAN(WW(I+KK)/UU(I+KK))/PI,ATAN(DWDY/DUDY)/PI,90.-ATAN
C     ?(CC1C2(I+KK)/CC2C3(I+KK))/PI+GAM
      GO TO 13
12     CONTINUE
       XXX=9999.
C       WRITE(6,*)'I AM HERE'
C        WRITE(2,10)BRAND-1,S3(I+KK),
C     ?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
C     **********************************
C      W-L-C METHOD
       GO TO 3001
C     **********************************
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
      ALFA=XXNNUU(IC-1)/RHO(IC-1)/UTWLC**3*PRES(IC-1)
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
C     **********************************
3001  CONTINUE
C     **********************************
C     **********************************
C     P-J MODEL
C     GO TO 3009
C     **********************************
      DPDXW=PRES(IC-1)*COS(BBETAW(IC-1))
     ?+PRES(IC-1+8)*SIN(BBETAW(IC-1))
      DPDZW=PRES(IC-1+8)*COS(BBETAW(IC-1))
     ?-PRES(IC-1)*SIN(BBETAW(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 91 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 91
      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))
91    CONTINUE
87    CONTINUE
C     **********************************
C     **********************************
C     **********************************
C     **********************************
C     EAST-HOXEY MODEL
3009  CONTINUE
      GO TO 3000
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
      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   CONTINUE
      STOP
      END
