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) ?,AALFA(100) BRAND=0. I=0 IC=0 DO 36 JCONT=1,14 36 READ(4,*)XX(JCONT),PRES(JCONT) DO 37 JJCC=15,20 37 PRES(JJCC)=1.00 CONT=0. 100 READ(1,*,END=300)Y,GAM,BET,SU,U,V,W,C1,C2,C3,C1C2,C1C3,C2C3 IF(Y.NE.0.)GO TO 121 IC=IC+1 BRAND=BRAND+1 IF(Y.EQ.0.)READ(1,*)NST,NX,Z,ALFA,UDELTA,UTAU,TOTAL,TEMP,NPRESS C WRITE(7,*)UTAU PRESS=NPRESS PI=3.1415926/180. XNU=(.1716*((273.15+TEMP)/273.15)**(1.5)*383.7*1E-04/(273.15+TEMP+ ?110.6))/(PRESS*100./287./(273.15+TEMP)) BETAW=(TOTAL-ALFA)/1.11111111 ALFA=ALFA/1.1111111111 AALFA(IC)=ALFA*PI XXNNUU(IC)=XNU BBETAW(IC)=BETAW*PI UT(IC)=UTAU RHO(IC)=PRESS*100./(287.*(TEMP+273.15)) PRES(IC)=PRES(IC)/2.*RHO(1)*22.67*22.67 PRES(IC+7)=PRES(IC+7)/2.*RHO(1)*22.67*22.67 C WRITE(7,*)PRES(IC),PRES(IC+7) IF(Y.EQ.0..AND.BRAND.NE.1.) GO TO 200 IF(Y.EQ.0.)GO TO 100 121 I=I+1 Y=Y/1000. Q(I)=SQRT(U*U+W*W)*UDELTA YY(I)=Y UU(I)=U*UDELTA VV(I)=V*UDELTA WW(I)=W*UDELTA CC1(I)=C1 CC2(I)=C2 CC3(I)=C3 CC1C2(I)=C1C2 CC1C3(I)=C1C3 CC2C3(I)=C2C3 S1=Y S2=Q(I)*COS((BETAW-GAM)*PI)/UTAU S3(I)=UTAU*Y/XNU S4=U*UDELTA/COS(BETAW*PI)/UTAU S5=U*UDELTA/UTAU S6=U*UDELTA/UTAU/SQRT(COS(BETAW*PI)) S7=Y*UTAU*SQRT(COS(BETAW*PI))/XNU S8=Q(I)/UTAU IF(BETAW.EQ.0.) GO TO 1111 S9=W*UDELTA/UTAU/SQRT(SIN(BETAW*PI)) S10=Y*UTAU*SQRT(SIN(BETAW*PI))/XNU C ************************** C FOR VAN DEN BERG'S MODEL 1111 CONTINUE GO TO 111 C ************************** C BETAX=XXNNUU(IC)/UT(IC)**2*UTX(IC) C BETAZ=XXNNUU(IC)/UT(IC)**2*UTZ(IC) BETAX=0. BETAZ=0. DPDXW=PRES(IC)*COS(BBETAW(IC)+AALFA(IC)) ? +PRES(IC+7)*SIN(BBETAW(IC)+AALFA(IC)) DPDZW=PRES(IC+7)*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 C WRITE(7,*)'ALFAX,ALFAZ,DPDXW,DPDZW',ALFAX,ALFAZ,DPDXW,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,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) GO TO 100 300 CONT=1 IC=IC+1 BRAND=BRAND+1 200 CONTINUE N=I 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 IF(CC1C2(I+KK).NE.0.AND.CC2C3(I+KK).NE.0.)GO TO 11 IF(CC1C2(I+KK).EQ.0.OR.CC2C3(I+KK).EQ.0.)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+GAM 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(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 DPDXS=PRES(IC-1)*COS(AALFA(IC-1)) ?+PRES(IC-1+7)*SIN(AALFA(IC-1)) ALFA=XXNNUU(IC-1)/RHO(IC-1)/UTWLC**3*DPDXS C PRINT *,'ALFA',ALFA*YP,IC-1 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 4001 C ********************************** DPDXW=PRES(IC-1)*COS(BBETAW(IC-1)+AALFA(IC-1)) ?+PRES(IC-1+7)*SIN(BBETAW(IC-1)+AALFA(IC-1)) DPDZW=PRES(IC-1+7)*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 ********************************** 4001 CONTINUE C ********************************** GO TO 3000 C ********************************** C EAST-HOXEY MODEL C ********************************** UBL=SQRT(UU(N)**2+WW(N)**2) IF((IC-1).EQ.1)U0=UBL C WRITE(7,*)U0 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 STOP END