      REAL YY(220),UU(220),VV(220),WW(220),CC1(100),CC2(100),CC3(100)   3DS00010
     ?,CC2C3(100),CC1C3(100),CC1C2(100),S3(220),Q(220),PRES(100)        3DS00020
     ?,XX(200),UTX(16),UTZ(16),UT(16),XXNNUU(16),RHO(16),BBETAW(8)      3DS00030
     ?,TEMP(12),uuun(220),wwun(220)            
      CHARACTER*70 AAAA,BBBB                                              3DS00050
      open (unit=4,file='semgr.data')
      open (unit=1,file='ldv1.data')
      open (unit=3,file='semvdb1.out')
      open (unit=18,file='semvdb2.out')
      open (unit=19,file='semvdb3.out')
      open (unit=20,file='semvdb4.out')
      open (unit=7,file='semwlc1.out')
      open (unit=15,file='semwlc2.out')
      open (unit=16,file='semwlc3.out')
      open (unit=17,file='semwlc4.out')
      open (unit=8,file='res1.out')
      open (unit=9,file='res2.out')
      open (unit=10,file='res3.out')
      open (unit=11,file='res4.out')
      open (unit=12,file='res5.out')
      open (unit=13,file='sempj1.out')
      open (unit=21,file='sempj2.out')
      open (unit=14,file='semeh.out')    
      print *,'here 1'
      s9=0.
      s10=0.
      brand=0.
      i=0
      ic=0
      DO 36 JCONT=1,12                                                  3DS00110
 36   READ(4,*)XX(JCONT),PRES(JCONT) 
C     DO 37 JCONT1=1,8                                                  3DS00130
C37   READ(5,*)XX(JCONT1),UTX(JCONT1),UTZ(JCONT1)                       3DS00140
      CONT=0.                                                           3DS00150
100   READ(1,*,END=300)y,u,v,w,uun,vun,wun,u2,v2,w2,uv,uw,vw,squun,
     ?sqvun,sqwun,uvun,uwun,vwun
      read(1,*)xj1,xj2,xj3,xj4,xj5
c     print *,Y  
C      WRITE(5,*)'WRITING Y VALUE',Y                                    3DS00180
      IF(Y.NE.0.)GO TO 121                                              3DS00190
      BRAND=BRAND+1.                                                    3DS00200
C      WRITE(6,*)BRAND                                                  3DS00210
      IC=IC+1                                                           3DS00220
      READ(1,*)NST,UREF,UTAU,TEM,BETAW,PRESS,BFREE        
c     print *,nst                                                       3DS00230
      PI=3.1415926/180.                                                 3DS00240
      RO=PRESS/287./(273.15+TEM)                                        3DS00250
      XNU=(.1716*((273.15+TEM)/273.15)**(1.5)*383.7*1E-04/(273.15+TEM+  3DS00260
     ?110.6))/RO                                                        3DS00270
      PRES(IC)=PRES(IC)*1/2.*RO*27*27*100                               3DS00280
      PRES(IC+6)=PRES(IC+6)*1/2.*RO*27*27*100                           3DS00290
      RE=UREF/XNU                                                       3DS00300
C      WRITE(5,*)NST,UREF,UTAU,TEM,BETAW,PRESS,XNU                      3DS00310
      BETAW=BETAW*PI                                                    3DS00320
      BFREE=BFREE*PI                                                    3DS00330
      TEMP(IC)=PRES(IC)*COS(BFREE)+PRES(IC+6)*SIN(BFREE)                3DS00340
      PRES(IC+6)=PRES(IC+6)*COS(BFREE)-PRES(IC)*SIN(BFREE)              3DS00350
      PRES(IC)=TEMP(IC)                                                 3DS00360
      XXNNUU(IC)=XNU                                                    3DS00370
      BBETAW(IC)=BETAW                                                  3DS00380
C     UTAU=UTAU+.5*UTAU                                                 3DS00390
      UT(IC)=UTAU                                                       3DS00400
      RHO(IC)=RO                                                        3DS00410
C     WRITE(6,*)'UT,RHO,XNU',UT(IC),RHO(IC),XXNNUU(IC)                  3DS00420
      IF(Y.EQ.0..AND.BRAND.NE.1.) GO TO 200                             3DS00430
      IF(Y.EQ.0.)GO TO 100                                              3DS00440
121   I=I+1                                                             3DS00450
      uuun(i)=uun              
      wwun(i)=wun                                                       3DS00480
      Q(I)=SQRT(U**2+W**2)*UREF                                         3DS00490
      Y=Y*2.54/100.                                                     3DS00500
      YY(I)=Y                                                           3DS00510
      U=U*UREF                                                          3DS00520
      W=W*UREF                                                          3DS00530
      UU(I)=U                                                           3DS00540
      WW(I)=W                                                           3DS00550
      BETA=w/u                                                          3DS00460
      write(6,*)'beta',beta
      IF(BETA.EQ.999.)I=I-1                                             3DS00470
      IF(BETA.EQ.999.)GO TO 500
C     CC1(I)=US                                                         3DS00560
C     CC2(I)=VS                                                         3DS00570
C     CC3(I)=WS                                                         3DS00580
C     CC1C2(I)=UV                                                       3DS00590
C     CC1C3(I)=UW                                                       3DS00600
C     CC2C3(I)=VW                                                       3DS00610
      S1=Y                                                              3DS00620
      S2=Q(I)*COS(BETAW-BETA)/UTAU                                      3DS00630
      S3(I)=UTAU*Y/XNU                                                  3DS00640
C     WRITE(5,*)UTAU,Y,XNU,S3(I)                                        3DS00650
      S4=U/COS(BETAW)/UTAU                                              3DS00660
      S5=U/UTAU                                                         3DS00670
      S6=U/UTAU/SQRT(COS(BETAW))                                        3DS00680
      S7=Y*UTAU*SQRT(COS(BETAW))/XNU                                    3DS00690
      S8=Q(I)/UTAU                                                      3DS00700
       IF(BETAW.EQ.0.) GO TO 111                                        3DS00710
       IF(BETAW.EQ.0.)S10=999.                                          3DS00720
       IF(BETAW.EQ.0.)S9=999.                                           3DS00730
       S9=-W/UTAU/SQRT(SIN(-BETAW))                                     3DS00740
       S10=Y*UTAU*SQRT(SIN(-BETAW))/XNU                       
c  ******************
c uncer for coles 
       d1=u/q(i)*cos(betaw-beta)/utau
       d2=w/q(i)*cos(betaw-beta)/utau
       d3=q(i)/utau*(-sin(betaw-beta))
       d4=q(i)/utau*sin(betaw-beta)
       d5=-q(i)*cos(betaw-beta)/utau**2
       delbet=sqrt((u/q(i)**2*wwun(i))**2+(-w/q(i)**2*uuun(i))**2)
       delbwall=2*3.1415926/180.
       delqstar=utau*.05
       delqp=sqrt((d1*uuun(i))**2+(d2*wwun(i))**2+(d3*delbwall)**2
     ? +(d4*delbet)**2+(d5*delqstar)**2)      
       write(8,4)s3(i),s2,delqp		
c  ************** 
c uncer for johnston
       d1=1/utau/cos(betaw)
       d2=-u/utau**2/cos(betaw)
       d3=-sin(betaw)*u/utau/(cos(betaw))**2
       delup=sqrt((d1*uuun(i))**2+(d2*delqstar)**2+(d3*delbwall)**2)
       write(9,4)s3(i),s4,delup
c  ****************
c uncer for hornung-joubert
       d1=u/q(i)/utau
       d2=w/q(i)/utau
       d3=-q(i)/utau**2
       delqp=sqrt((d1*uuun(i))**2+(d2*wwun(i))**2+(d3*delqstar)**2) 
       write(10,4)s3(i),s8,delqp 
c  ****************
c uncer for p-k and chan-swamy
       d1=1/utau/sqrt(cos(betaw))
       d2=-u/utau**2/sqrt(cos(betaw))
       d3=sin(betaw)*u/2./utau/(cos(betaw))**(1.5)
       d4=1/utau/sqrt(sin(-betaw))
       d5=-w/utau**2/sqrt(sin(-betaw))
       d6=w/utau*(-1./2.*cos(betaw)/(sin(-betaw))**1.5)
       delup=sqrt((d1*uuun(i))**2+(d2*delqstar)**2+(d3*delbwall)**2)
       delwp=sqrt((d4*wwun(i))**2+(d5*delqstar)**2+(d6*delbwall)**2)
       write(11,4)s7,s6,delup
       write(12,4)s10,s9,delwp 
C      WRITE(2,4)BRAND,S1,S2,S3(I),S4,S5,S6,S7,S8,S9,S10                3DS00760
4     FORMAT(3F12.5)                                                    3DS00820
C     **************************                                        3DS00830
C     FOR VAN DEN BERG'S MODEL                                          3DS00840
c     GO TO 3000  
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+6)*SIN(BBETAW(IC))    
      DPDZW=PRES(IC+6)*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)                          
C     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)               
      UPX=Q(I)*COS(BBETAW(IC)-BETA)/UT(IC)             
C     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)              
      UPZ=Q(I)*SIN(BBETAW(IC)-BETA)/UT(IC)                      
       d1=uu(i)/q(i)*cos(bbetaw(ic)-beta)/ut(ic)
       d2=ww(i)/q(i)*cos(bbetaw(ic)-beta)/ut(ic)
       d3=q(i)/ut(ic)*(-sin(bbetaw(ic)-beta))
       d4=q(i)/ut(ic)*sin(bbetaw(ic)-beta)
       d5=-q(i)*cos(bbetaw(ic)-beta)/ut(ic)**2
       delbet=sqrt((uu(i)/q(i)**2*wwun(i))**2
     ?+(-ww(i)/q(i)**2*uuun(i))**2)
       delbwall=2*3.1415926/180.
       delqstar=ut(ic)*.05
       delqp=sqrt((d1*uuun(i))**2+(d2*wwun(i))**2+(d3*delbwall)**2
     ? +(d4*delbet)**2+(d5*delqstar)**2)      
C      WRITE(6,555)BRAND,YP,UPX,RHSX,UPZ,RHSZ                
       WRITE(3,555)YP,UPX,delqp  
      d1=1./.41*(1./yp+1/2.*alfax) 
      d2=1./.41*(1./2.*yp)
      delyp=sqrt((yy(i)/xxnnuu(ic)*delqstar)**2)
      delpx=dpdxw*2/100.
      delpz=dpdzw*2/100.
      delalx=sqrt((-3*xxnnuu(ic)/rho(ic)/ut(ic)**4*dpdxw*delqstar)**2
     ?+(xxnnuu(ic)/rho(ic)/ut(ic)**3*delpx)**2) 
      delrhsx=sqrt((d1*delyp)**2+(d2*delalx)**2)               
       WRITE(18,555)YP,RHSX,delrhsx  
       d1=uu(i)/sqrt(uu(i)**2+ww(i)**2)*sin(bbetaw(ic)-beta)/ut(ic)
       d2=ww(i)/sqrt(uu(i)**2+ww(i)**2)*sin(bbetaw(ic)-beta)/ut(ic)
       d3=sqrt(uu(i)**2+ww(i)**2)/ut(ic)*cos(bbetaw(ic)-beta)
       d4=-sqrt(uu(i)**2+ww(i)**2)/ut(ic)*cos(bbetaw(ic)-beta)
       d5=-sqrt(uu(i)**2+ww(i)**2)/ut(ic)**2*sin(bbetaw(ic-1)-beta)
       delbet=sqrt((uu(i)/(uu(i)**2+ww(i)**2)*wwun(i))**2+
     ? (-ww(i)/(uu(i)**2+ww(i)**2)*uuun(i))**2)
       delupz=sqrt((d1*uuun(i))**2+(d2*wwun(i))**2+(d3*delbwall)**2
     ? +(d4*delbet)**2+(d5*delqstar)**2)                 
       WRITE(19,555)YP,UPZ,delupz   
      d1=1/.41*yp
      d2=1/.41*alfaz     
      delalz=sqrt((-3*xxnnuu(ic)/rho(ic)/ut(ic)**4*dpdzw*delqstar)**2
     ?+(xxnnuu(ic)/rho(ic)/ut(ic)**3*delpz)**2) 
      delrhsz=sqrt((d1*delalz)**2+(d2*delyp)**2)
      write(20,555)YP,RHSZ,delrhsz                   
555   FORMAT(6(1X,F10.3))                                
C     **************************                        
C     **************************                  
111   CONTINUE                                
500   GO TO 100                           
300   CONT=1                                                    
c     print *,CONT
      IC=IC+1                                                    
      BRAND=BRAND+1                                       
200   CONTINUE                                         
       N=I                                       
C      WRITE(6,*)N                            
       KK=-1                            
C      DO 9 I=1,N-4               
C       WRITE(6,*)'I AM GOING N'                                        3DS01190
C53    KK=KK+1                                                          3DS01200
C      IF(KK.GT.2)KK=2                                                  3DS01210
C        CALL PARFIT(I,N,UU,YY,B0,B1,B2)                                3DS01220
C        DUDY=B1+2*B2*YY(I+KK)                                          3DS01230
C          CALL PARFIT(I,N,WW,YY,C0,C1,C2)                              3DS01240
C          DWDY=C1+2*C2*YY(I+KK)                                        3DS01250
C             IF(DUDY.EQ.0.AND.DWDY.EQ.0.)DUDY=10**50                   3DS01260
C             IF(DUDY.EQ.0.AND.DWDY.EQ.0.)DWDY=10**50                   3DS01270
C     IF(CC1C2(I+KK).NE.999.AND.CC2C3(I+KK).NE.999.)GO TO 11            3DS01280
C     IF(CC1C2(I+KK).EQ.999.OR.CC2C3(I+KK).EQ.999.)GO TO 12             3DS01290
10     FORMAT(5F12.5)                                                   3DS01300
C11    CONTINUE                                                         3DS01310
C       IF(CC1C2(I+KK)/CC2C3(I+KK).LE.0.) WRITE(2,10)BRAND-1,S3(I+KK),  3DS01320
C     ?ATAN(WW(I+KK)/UU(I+KK))/PI,ATAN(DWDY/DUDY)/PI,ATAN               3DS01330
C     ?(CC2C3(I+KK)/CC1C2(I+KK))/PI                                     3DS01340
C      IF(CC1C2(I+KK)/CC2C3(I+KK).GT.0.) WRITE(2,10)BRAND-1,S3(I+KK),   3DS01350
C     ?ATAN(WW(I+KK)/UU(I+KK))/PI,ATAN(DWDY/DUDY)/PI,90.-ATAN           3DS01360
C     ?(CC1C2(I+KK)/CC2C3(I+KK))/PI+GAM                                 3DS01370
C     GO TO 13                                                          3DS01380
C12    CONTINUE                                                         3DS01390
       XXX=9999.                                                        3DS01400
C       WRITE(6,*)'I AM HERE'                                           3DS01410
C        WRITE(2,10)BRAND-1,S3(I+KK),                                   3DS01420
C     ?ATAN(WW(I+KK)/UU(I+KK))/PI,ATAN(DWDY/DUDY)/PI,XXX                3DS01430
C13    IF(KK.EQ.2)GO TO 9                                               3DS01440
C      GO TO 53                                                         3DS01450
C9     CONTINUE                                                         3DS01460
C     **********************************                                3DS01470
3000  CONTINUE                                                          3DS01480
C      W-L-C METHOD                                                     3DS01490
c      GO TO 3001                                                       3DS01500
C     **********************************                                
C     PRINT *,BBETAW(IC-1),UT(IC-1),RHO(IC-1),PRES(IC-1)               
      uext=sqrt(((uu(n)**2+ww(n)**2)+
     ?(uu(n-1)**2+ww(n-1)**2)+
     ?(uu(n-2)**2+ww(n-2)**2))
     ?/3.) 
      blcont=0
      do 253 ix=1,n
      uuuu=uu(ix)**2+ww(ix)**2
      if(blcont.eq.1) goto 253 
      if(uuuu/uext**2.gt..99)idelt1=ix-1
      if(uuuu/uext**2.gt..99)idelt2=ix
253   if(uuuu/uext**2.gt..99)blcont=1
      uuuu1=uu(idelt1)**2+ww(idelt1)**2
      uuuu2=uu(idelt2)**2+ww(idelt2)**2
      blthick=(yy(idelt2)-yy(idelt1))/(uuuu2-uuuu1)*(.99*uext**2
     ?-uuuu1)+yy(idelt1)
      write(6,*)uext,blthick,yy(idelt1),yy(idelt2)
      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)/blthick)**2       
c **************************************
c uncer for white-lessmann-christoph  for r.h.s. of  u
      d1=1/2.*yp/sqrt(1+alfa*yp)
      d2=1/2.*alfa/sqrt(1+alfa*yp)
      d3=1/2.*.1108/sqrt(1+.1108*alfa)
      d4=xxnnuu(ic-1)/rho(ic-1)*pres(ic-1)*(-3)/utwlc**4
      d5=3/2.*alfa*sin(bbetaw(ic-1))/utwlc/sqrt(cos(bbetaw(ic-1))) 
      d6=alfa/pres(ic-1)
      d7=yy(iwlc)/xxnnuu(ic-1)*cos(bbetaw(ic-1))
      delqstar=ut(ic-1)*5/100.
      delbetw=2*3.1415926/180.
      delpx=pres(ic-1)*2/100.
      delalfa=sqrt((d4*delqstar)**2+(d5*delbetw)**2+(d6*delpx)**2)
      delyp=sqrt((d7*delqstar)**2)
      dels0=sqrt((d3*delalfa)**2)
      dels=sqrt((d1*delalfa)**2+(d2*delyp)**2)
      d8=1/.41*(2+2/((s-1)*(s+1)))
      d9=-1/.41*(2+2/((s0+1)*(s0-1)))
      delrhs=sqrt((d8*dels)**2+(d9*dels0)**2) 
      WRITE(15,556)YP,RHSU,delrhs
c uncer for left hand side of u
       d1=1/ut(ic-1)/sqrt(cos(bbetaw(ic-1)))
       d2=-uu(iwlc)/ut(ic-1)**2/sqrt(cos(bbetaw(ic-1)))
       d3=sin(bbetaw(ic-1))*uu(iwlc)
     ?/2./ut(ic-1)/(cos(bbetaw(ic-1)))**(1.5)
       delup=sqrt((d1*uuun(iwlc))**2+(d2*delqstar)**2+(d3*delbetw)**2)
      WRITE(7,556)YP,UP,delup  
      d1=1/ut(ic-1)/sqrt(cos(bbetaw(ic-1)))
      d2=-ww(iwlc)/ut(ic-1)**2/sqrt(cos(bbetaw(ic-1)))
      d3=ww(iwlc)/2./ut(ic-1)*sin(bbetaw(ic-1))/
     ?(cos(bbetaw(ic-1)))**1.5
      delwp=sqrt((d1*wwun(iwlc))**2+(d2*delqstar)**2+(d3*delbetw)**2)
      WRITE(16,556)YP,WP,delwp
      d1=tan(bbetaw(ic-1))*(1-yy(iwlc)/yy(n))**2
      d2=uu(iwlc)/ut(ic-1)/sqrt(cos(bbetaw(ic-1)))*(1+tan(bbetaw(ic-1))
     ?**2)*(1-yy(iwlc)/yy(n))**2
      delwp=sqrt((d1*delup)**2+(d2*delbetw)**2)
      WRITE(17,556)YP,RHSW,delwp                    
556   FORMAT(3(1X,F10.4))                                    
123   CONTINUE                                          
23    CONTINUE                        
C     **********************************                      
3001  CONTINUE                                       
C     **********************************         
C     *******************************
C     P-J MODEL                             
c     GO TO 3002                    
C     **********************************        
      DPDXW=PRES(IC-1)*COS(BBETAW(IC-1))                                
     ?+PRES(IC-1+6)*SIN(BBETAW(IC-1))                           
      DPDZW=PRES(IC-1+6)*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',teta                        
c     print *,'TA',ACOS(TETA)/PI          
      ALFA=SQRT(DPDXW*DPDXW+DPDZW*DPDZW)/RHO(IC-1)                   
C     PRINT *,'ALFA',ALFA                                       
c     go to 5001
      SUM=0.                                                  
      SUM1=0. 
      sum2=0.
      YPRE=0.                                                     
      RHSPRE=0.                                              
      pre1=0.
      pre2=0.
      pre3=0.
      srhs1=0.
      srhs2=0.
      srhs3=0.
      DO 49 IPJ=1,N                                     
      YP=YY(IPJ)*UT(IC-1)/XXNNUU(IC-1)              
      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.))                      
      drhsdt=1/4./.41/yy(ipj)*(2*sin(teta)*alfa*yy(ipj)/ut(ic-1)**2)
     ?*(1-2*cos(teta)*(alfa*yy(ipj)/ut(ic-1)**2)+(alfa/ut(ic-1)**2)**2
     ?*yy(ipj)**2)**-.75
      drhsda=1/.41/4./yy(ipj)*(-2.*cos(teta)*yy(ipj)/ut(ic-1)**2
     ?+2*(alfa/ut(ic-1)**2)*yy(ipj)**2/ut(ic-1)**2)
     ?*(1-2*cos(teta)*(alfa*yy(ipj)/ut(ic-1)**2)+(alfa/ut(ic-1)**2)**2
     ?*yy(ic-1)**2)-.75
      drhsdq=1/.41/4./yy(ipj)*(1-2*cos(teta)*(alfa/ut(ic-1)**2)*yy(ipj)
     ?+(alfa/ut(ic-1)**2)**2*yy(ipj)**2)**-.75*
     ?(-2*cos(teta)*alfa*yy(ipj)*(-2/ut(ic-1)**3)+2*(alfa*yy(ipj)**2/
     ?ut(ic-1)**2)*(-2*alfa/ut(ic-1)**3))
      deltet=2*3.1415926/180.
      delalf=alfa*.02
      delqstar=ut(ic-1)*.05
      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 
      delrl=sqrt(((((ww(ipj)-ww(ipj-1))*wwun(ipj))**2
     ?+((ww(ipj)-ww(ipj-1))*
     ?wwun(ipj-1))**2+((uu(ipj)-uu(ipj-1))*uuun(ipj))**2+
     ?((uu(ipj)-uu(ipj-1))*uuun(ipj-1))**2)/ut(ic-1)**2)+
     ?(rl/ut(ic-1)*delqstar)**2)
      sum2=sum2+rl+delrl 
      SUM=SUM+(RHS+RHSPRE)/2.*(YP-YPRE) 
      srhs1=srhs1+(drhsdt+pre1)/2.*(yy(ipj)-yy(ipj-1))
      srhs2=srhs2+(drhsda+pre2)/2.*(yy(ipj)-yy(ipj-1))
      srhs3=srhs3+(drhsdq+pre3)/2.*(yy(ipj)-yy(ipj-1))   
      YPRE=YP                                                     
      RHSPRE=RHS                                      
      pre1=drhsdt
      pre2=drhsda
      pre3=drhsdq                  
      delrhs=sqrt(srhs1**2*deltet**2+srhs2**2*delalf**2+srhs3**2*
     ?delqstar**2)
      WRITE(13,557)YP,SUM1,sum2-sum1                              
c     WRITE(21,557)YP,SUM+5.61,delrhs                              
557   FORMAT(4(1X,F10.4))                                               
49    CONTINUE                                                        
C     **********************************                                        
5001  continue
      SUM1=0. 
      sum2=0.
      pre1=0.
      pre2=0.
      pre3=0.
      srhs1=0.
      srhs2=0.
      srhs3=0.
      delrhs=0
      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.))                                        
      sum1=sum
      SUM=SUM+(RHS+RHSPRE)/2.*(YP-YPRE)                                         
      drhsdt=1/4./.41/yp*(2*sin(teta)*alfa*xxnnuu(ic-1)*yp/ut(ic-1)**3)
     ?*(1-2*cos(teta)*(alfa*xxnnuu(ic-1)*yp/ut(ic-1)**3)+(alfa/ut(ic-1)*
     ?*3*xxnnuu(ic-1))**2
     ?*yp**2)**-.75
      drhsda=1/.41/4./yp*(-2.*cos(teta)*xxnnuu(ic-1)*yp/ut(ic-1)**3
     ?+2*(alfa*xxnnuu(ic-1)/ut(ic-1)**3)*yp**2*xxnnuu(ic-1)/ut(ic-1)**3)
     ?*(1-2*cos(teta)*(alfa*xxnnuu(ic-1)*yp/ut(ic-1)**3)+(alfa/ut(ic-1)*
     ?*3*xxnnuu(ic-1))**2
     ?*yp**2)**-.75
      drhsdq=(sum-sum1)/(yp-ypre)*yp/ut(ic-1)
      deltet=2*3.1415926/180.
      delalf=alfa*.02
      delqstar=ut(ic-1)*.05
      srhs1=srhs1+(drhsdt+pre1)/2.*(yp-ypre)
      srhs2=srhs2+(drhsda+pre2)/2.*(yp-ypre)
      srhs3=srhs3+(drhsdq+pre3)/2.*(yp-ypre)   
      delrhs=sqrt(srhs1**2*deltet**2+srhs2**2*delalf**2+srhs3**2*
     ?delqstar**2)
92    YPRE=YP                                                                   
      RHSPRE=RHS                                                                
      pre1=drhsdt
      pre2=drhsda
      pre3=drhsdq                  
      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(21,561)YP,SUM+5.61,delrhs                                          
561   FORMAT(3(1X,F10.4))                                                       
88    CONTINUE                                                                  
87    CONTINUE                                                                  
C     **********************************                                        
C     **********************************                             
3002  CONTINUE                                                   
C     **********************************                      
C     EAST-HOXEY MODEL                                       
c     GO TO 3003                                                        3DS02140
C     **********************************                                3DS02150
      UBL=SQRT(UU(N)**2+WW(N)**2)                                      
      IF((IC-1).EQ.1)U0=UBL 
      u0=27.47                                          
      XK=19.45                                                       
      YPRE=0.                                                       
      DELPRE1=0.                                                  
      DELPRE2=0.
      delpuun=0
      delpwun=0
      sumuun=0 
      sumwun=0                                                
      SUM1=0.           
      SUM2=0.         
      DO 63 IEH=1,N                                                    
      DEL1=(1-UU(IEH)/UBL)                                         
      DEL2=-WW(IEH)/Ubl 
      deuun=-1/ubl*uuun(ieh)
      dewun=-1/ubl*wwun(ieh)                                     
      Y1=YY(IEH)                                       
      SUM1=SUM1+(DEL1+DELPRE1)/2.*(Y1-YPRE)         
      SUM2=SUM2+(DEL2+DELPRE2)/2.*(Y1-YPRE)
      sumuun=sumuun+(delpuun+deuun)/2.*(y1-ypre)
      sumwun=sumwun+(delpwun+dewun)/2.*(y1-ypre)
      delpuun=deuun
      delpwun=dewun
      DELPRE1=DEL1                                                
      DELPRE2=DEL2                                             
63    YPRE=Y1       
      GAMA=-SUM2/SUM1 
      write(6,*)sum1,sum2
      beta0=asin(1/(xk*ut(ic-1)/u0)*sin(gama))-gama
      delgama=sqrt((sum1/(sum1**2+sum2**2)*sumuwn)**2
     ?+(sum2/(sum1**2+sum2**2)*sumuun)**2)
      den=sqrt((xk*ut(ic-1))**2-(u0*sin(gama))**2)
      d1=sin(gama)/den
      d2=u0/den
      d3=-u0*sin(gama)/xk/den
      d4=-u0*sin(gama)/ut(ic-1)/den
      deluz=.02*u0
      dxk=5.5/100*xk
      delqstar=5/100.*ut(ic-1)
      delbetz=sqrt((d1*deluz)**2+((d2-1)*delgama)**2+(d3*dxk)**2
     ?+(d4*delqstar)**2)
      DO 64 IEH=1,N                             
      XLHS=UU(IEH)/COS(BETA0)/UT(IC-1)        
c uncer for east-hoxey
       d1=1/ut(ic-1)/cos(beta0)
       d2=-uu(ieh)/ut(ic-1)**2/cos(beta0)
       d3=-sin(beta0)*uu(ieh)/ut(ic-1)/(cos(beta0))**2
       delup=sqrt((d1*uuun(i))**2+(d2*delqstar)**2+(d3*delbetz)**2)
      YP=YY(IEH)*UT(IC-1)/XXNNUU(IC-1)    
      WRITE(14,58)YP,XLHS,delup       
58    FORMAT(3(1X,F10.4))           
64    CONTINUE                     
C     **************************                                        3DS02410
3003  CONTINUE                                                          3DS02420
C     **************************                                        3DS02430
C     **************************                                        3DS02440
      i=0                                                               3DS02450
      IF(CONT.NE.0.)GO TO 301                                           3DS02460
      GO TO 100                                                         3DS02470
301   CONTINUE                                                          3DS02480
      STOP                                                              3DS02490
      end
