From 0d4f43d355de79178b1142e9735902cf641670b6 Mon Sep 17 00:00:00 2001 From: Dimitri Sokolyuk Date: Mon, 11 May 2009 00:27:49 +0000 Subject: Xfoil 6.97 --- orrs/src/ospres.f | 246 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 246 insertions(+) create mode 100644 orrs/src/ospres.f (limited to 'orrs/src/ospres.f') diff --git a/orrs/src/ospres.f b/orrs/src/ospres.f new file mode 100644 index 0000000..c4632f0 --- /dev/null +++ b/orrs/src/ospres.f @@ -0,0 +1,246 @@ + + SUBROUTINE OSPRES(NI,YI,UI, ALPHAR,ALPHAI, VTR,VTI, + & PTR,PTI ) + DIMENSION YI(NI), UI(NI) + DIMENSION VTR(NI), VTI(NI) + DIMENSION PTR(NI), PTI(NI) +C--------------------------------------------------------------------- +C Routine for calculating the Orr-Sommerfeld pressure profile. +C +C Input: +C ------ +C NI total number of points in profiles +C YI normal BL coordinate array +C UI mean flow u(y) profile +C ALPHAR real part of complex wavenumber +C ALPHAI imag. part of complex wavenumber +C VTR real part of perturbation y-velocity profile +C VTI imag. part of perturbation y-velocity profile +C +C Output: +C ------- +C PTR real part of perturbation pressure profile +C PTI imag. part of perturbation pressure profile +C--------------------------------------------------------------------- +C + INCLUDE 'OSPRES.INC' +C +C---- convergence tolerance + DATA EPS / 1.0E-4 / +C + IF(NI.GT.NMAX) STOP 'OSPRES: Array overflow.' +C + N = NI + DO 5 I=1, N + Y(I) = YI(I) + U(I) = UI(I) + VT(I) = CMPLX( VTR(I) , VTI(I) ) + 5 CONTINUE +C + ALPHA = CMPLX(ALPHAR,ALPHAI) +C +C---- set number of righthand sides + NRHS = 1 +C + DO I=1, N + F0(I) = 0. + F1(I) = 0. + ENDDO + ISOL = 0 +C + CALL SETUP_P + CALL SOLVE_P + CALL UPDATE_P +C + DO 200 I=1, N + PTR(I) = REAL(F0(I)) + PTI(I) = IMAG(F0(I)) + 200 CONTINUE +C + RETURN + END ! OSPRES + + + SUBROUTINE SETUP_P + INCLUDE 'OSPRES.INC' + COMPLEX VTA +C +C---- zero out A,B,C blocks and righthand sides R + DO 20 I=1, N + DO 201 J=1, 2 + DO 2001 K=1, 2 + A(J,K,I) = (0.0,0.0) + B(J,K,I) = (0.0,0.0) + C(J,K,I) = (0.0,0.0) + 2001 CONTINUE + DO 2002 K=1, NRMAX + R(J,K,I) = (0.0,0.0) + 2002 CONTINUE + 201 CONTINUE + 20 CONTINUE +C + I = 1 +C +C---- set 1st wall BC + R(2,1,I) = F1(I) + A(2,2,I) = 1.0 +C +C---- set interior equations + DO 50 I=1,N-1 +C + DY = Y(I+1) - Y(I) + DU = U(I+1) - U(I) +C +C--------------------------------------------------------------- +C + R(1,1,I) = F0(I+1) - F0(I) - 0.5*DY*(F1(I+1)+F1(I)) + A(1,1,I) = -1.0 + C(1,1,I) = 1.0 + A(1,2,I) = -0.5*DY + C(1,2,I) = -0.5*DY +C--------------------------------------------------------------- +C + R(2,1,I+1) = F1(I+1) - F1(I) - 0.5*DY*(F0(I+1)+F0(I))*ALPHA**2 + & + (0.0,1.0)*ALPHA*DU*(VT(I+1) + VT(I)) + B(2,1,I+1) = -0.5*DY*ALPHA**2 + A(2,1,I+1) = -0.5*DY*ALPHA**2 + B(2,2,I+1) = -1.0 + A(2,2,I+1) = 1.0 +C--------------------------------------------------------------- +C + 50 CONTINUE +C +C---- set asymptotic regularity conditions at outer edge +C + R(1,1,N) = F1(N) + F0(N)*ALPHA + A(1,1,N) = ALPHA + A(1,2,N) = 1.0 +C + RETURN + END ! SETUP + + + SUBROUTINE SOLVE_P + INCLUDE 'OSPRES.INC' + COMPLEX PIVOT, TEMP +C--------------------------------------------------- +C 2x2 complex tridiagonal block solver. +C--------------------------------------------------- +C +CCC** Forward sweep: Elimination of lower block diagonal (B's). + DO 1 I=1, N +C + IM = I-1 +C +C------ don't eliminate B1 block because it doesn't exist + IF(I.EQ.1) GO TO 12 +C +C------ eliminate Ci block, thus modifying Ai and Ri blocks + DO 111 L=1, 2 + K = 1 + A(K,L,I) = A(K,L,I) + & - B(K,1,I)*C(1,L,IM) + & - B(K,2,I)*C(2,L,IM) + K = 2 + A(K,L,I) = A(K,L,I) + & - B(K,1,I)*C(1,L,IM) + & - B(K,2,I)*C(2,L,IM) + 111 CONTINUE + DO 112 L=1, NRHS + K = 1 + R(K,L,I) = R(K,L,I) + & - B(K,1,I)*R(1,L,IM) + & - B(K,2,I)*R(2,L,IM) + K = 2 + R(K,L,I) = R(K,L,I) + & - B(K,1,I)*R(1,L,IM) + & - B(K,2,I)*R(2,L,IM) + 112 CONTINUE +C +C -1 +CCC---- multiply Ci block and righthand side Ri vectors by (Ai) +C using Gaussian elimination. +C + 12 CONTINUE +C + DO 13 KPIV=1, 2 +C + KP1 = KPIV+1 +C + PIVOT = 1.0/A(KPIV,KPIV,I) +C +C-------- normalize pivot row + DO 132 L=KP1, 2 + A(KPIV,L,I) = A(KPIV,L,I)*PIVOT + 132 CONTINUE +C + C(KPIV,1,I) = C(KPIV,1,I)*PIVOT + C(KPIV,2,I) = C(KPIV,2,I)*PIVOT +C + DO 134 L=1, NRHS + R(KPIV,L,I) = R(KPIV,L,I)*PIVOT + 134 CONTINUE +C +C-------- eliminate lower off-diagonal elements in Ai block + DO 135 K=KP1, 2 + TEMP = A(K,KPIV,I) + DO 1351 L=KP1, 2 + A(K,L,I) = A(K,L,I) - TEMP*A(KPIV,L,I) + 1351 CONTINUE + C(K,1,I) = C(K,1,I) - TEMP*C(KPIV,1,I) + C(K,2,I) = C(K,2,I) - TEMP*C(KPIV,2,I) + DO 1352 L=1, NRHS + R(K,L,I) = R(K,L,I) - TEMP*R(KPIV,L,I) + 1352 CONTINUE + 135 CONTINUE +C + 13 CONTINUE +C +C------ back substitute everything + DO 15 KPIV=1, 1, -1 + KP1 = KPIV+1 + DO 151 K=KP1, 2 + C(KPIV,1,I) = C(KPIV,1,I) - A(KPIV,K,I)*C(K,1,I) + C(KPIV,2,I) = C(KPIV,2,I) - A(KPIV,K,I)*C(K,2,I) + DO 1511 L=1, NRHS + R(KPIV,L,I) = R(KPIV,L,I) - A(KPIV,K,I)*R(K,L,I) + 1511 CONTINUE + 151 CONTINUE + 15 CONTINUE +C + 1 CONTINUE +C +CCC** Backward sweep: Back substitution using upper block diagonal (Ci's). + DO 2 I=N-1, 1, -1 + IP = I+1 + DO 21 L=1, NRHS + DO 211 K=1, 2 + R(K,L,I) = R(K,L,I) + & - (R(1,L,IP)*C(K,1,I) + R(2,L,IP)*C(K,2,I)) + 211 CONTINUE + 21 CONTINUE + 2 CONTINUE +C + RETURN + END ! SOLVE + + + SUBROUTINE UPDATE_P + INCLUDE 'OSPRES.INC' + COMPLEX DF0,DF1 +C + RLX = 1.0 +C +C---- perform Newton update on modes + DO 50 I=1, N + DF0 = -R(1,1,I) + DF1 = -R(2,1,I) +C + F0(I) = F0(I) + RLX*DF0 + F1(I) = F1(I) + RLX*DF1 +C + 50 CONTINUE +C + RETURN + END ! UPDATE + -- cgit v1.2.3