aboutsummaryrefslogtreecommitdiff
path: root/src/xbl.f
diff options
context:
space:
mode:
Diffstat (limited to 'src/xbl.f')
-rw-r--r--src/xbl.f1581
1 files changed, 1581 insertions, 0 deletions
diff --git a/src/xbl.f b/src/xbl.f
new file mode 100644
index 0000000..235870f
--- /dev/null
+++ b/src/xbl.f
@@ -0,0 +1,1581 @@
+C***********************************************************************
+C Module: xbl.f
+C
+C Copyright (C) 2000 Mark Drela
+C
+C This program is free software; you can redistribute it and/or modify
+C it under the terms of the GNU General Public License as published by
+C the Free Software Foundation; either version 2 of the License, or
+C (at your option) any later version.
+C
+C This program is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with this program; if not, write to the Free Software
+C Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+C***********************************************************************
+C
+ SUBROUTINE SETBL
+C-------------------------------------------------
+C Sets up the BL Newton system coefficients
+C for the current BL variables and the edge
+C velocities received from SETUP. The local
+C BL system coefficients are then
+C incorporated into the global Newton system.
+C-------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ INCLUDE 'XBL.INC'
+ REAL USAV(IVX,2)
+ REAL U1_M(2*IVX), U2_M(2*IVX)
+ REAL D1_M(2*IVX), D2_M(2*IVX)
+ REAL ULE1_M(2*IVX), ULE2_M(2*IVX)
+ REAL UTE1_M(2*IVX), UTE2_M(2*IVX)
+ REAL MA_CLMR, MSQ_CLMR, MDI
+C
+C---- set the CL used to define Mach, Reynolds numbers
+ IF(LALFA) THEN
+ CLMR = CL
+ ELSE
+ CLMR = CLSPEC
+ ENDIF
+C
+C---- set current MINF(CL)
+ CALL MRCL(CLMR,MA_CLMR,RE_CLMR)
+ MSQ_CLMR = 2.0*MINF*MA_CLMR
+C
+C---- set compressibility parameter TKLAM and derivative TK_MSQ
+ CALL COMSET
+C
+C---- set gas constant (= Cp/Cv)
+ GAMBL = GAMMA
+ GM1BL = GAMM1
+C
+C---- set parameters for compressibility correction
+ QINFBL = QINF
+ TKBL = TKLAM
+ TKBL_MS = TKL_MSQ
+C
+C---- stagnation density and 1/enthalpy
+ RSTBL = (1.0 + 0.5*GM1BL*MINF**2) ** (1.0/GM1BL)
+ RSTBL_MS = 0.5*RSTBL/(1.0 + 0.5*GM1BL*MINF**2)
+C
+ HSTINV = GM1BL*(MINF/QINFBL)**2 / (1.0 + 0.5*GM1BL*MINF**2)
+ HSTINV_MS = GM1BL*( 1.0/QINFBL)**2 / (1.0 + 0.5*GM1BL*MINF**2)
+ & - 0.5*GM1BL*HSTINV / (1.0 + 0.5*GM1BL*MINF**2)
+C
+C---- set Reynolds number based on freestream density, velocity, viscosity
+ HERAT = 1.0 - 0.5*QINFBL**2*HSTINV
+ HERAT_MS = - 0.5*QINFBL**2*HSTINV_MS
+C
+ REYBL = REINF * SQRT(HERAT**3) * (1.0+HVRAT)/(HERAT+HVRAT)
+ REYBL_RE = SQRT(HERAT**3) * (1.0+HVRAT)/(HERAT+HVRAT)
+ REYBL_MS = REYBL * (1.5/HERAT - 1.0/(HERAT+HVRAT))*HERAT_MS
+C
+ AMCRIT = ACRIT
+ IDAMPV = IDAMP
+C
+C---- save TE thickness
+ DWTE = WGAP(1)
+C
+ IF(.NOT.LBLINI) THEN
+C----- initialize BL by marching with Ue (fudge at separation)
+ WRITE(*,*)
+ WRITE(*,*) 'Initializing BL ...'
+ CALL MRCHUE
+ LBLINI = .TRUE.
+ ENDIF
+C
+ WRITE(*,*)
+C
+C---- march BL with current Ue and Ds to establish transition
+ CALL MRCHDU
+C
+ DO 5 IS=1, 2
+ DO 6 IBL=2, NBL(IS)
+ USAV(IBL,IS) = UEDG(IBL,IS)
+ 6 CONTINUE
+ 5 CONTINUE
+C
+ CALL UESET
+C
+ DO 7 IS=1, 2
+ DO 8 IBL=2, NBL(IS)
+ TEMP = USAV(IBL,IS)
+ USAV(IBL,IS) = UEDG(IBL,IS)
+ UEDG(IBL,IS) = TEMP
+ 8 CONTINUE
+ 7 CONTINUE
+C
+ ILE1 = IPAN(2,1)
+ ILE2 = IPAN(2,2)
+ ITE1 = IPAN(IBLTE(1),1)
+ ITE2 = IPAN(IBLTE(2),2)
+C
+ JVTE1 = ISYS(IBLTE(1),1)
+ JVTE2 = ISYS(IBLTE(2),2)
+C
+ DULE1 = UEDG(2,1) - USAV(2,1)
+ DULE2 = UEDG(2,2) - USAV(2,2)
+C
+C---- set LE and TE Ue sensitivities wrt all m values
+ DO 10 JS=1, 2
+ DO 110 JBL=2, NBL(JS)
+ J = IPAN(JBL,JS)
+ JV = ISYS(JBL,JS)
+ ULE1_M(JV) = -VTI( 2,1)*VTI(JBL,JS)*DIJ(ILE1,J)
+ ULE2_M(JV) = -VTI( 2,2)*VTI(JBL,JS)*DIJ(ILE2,J)
+ UTE1_M(JV) = -VTI(IBLTE(1),1)*VTI(JBL,JS)*DIJ(ITE1,J)
+ UTE2_M(JV) = -VTI(IBLTE(2),2)*VTI(JBL,JS)*DIJ(ITE2,J)
+ 110 CONTINUE
+ 10 CONTINUE
+C
+ ULE1_A = UINV_A(2,1)
+ ULE2_A = UINV_A(2,2)
+C
+C**** Go over each boundary layer/wake
+ DO 2000 IS=1, 2
+C
+C---- there is no station "1" at similarity, so zero everything out
+ DO 20 JS=1, 2
+ DO 210 JBL=2, NBL(JS)
+ JV = ISYS(JBL,JS)
+ U1_M(JV) = 0.
+ D1_M(JV) = 0.
+ 210 CONTINUE
+ 20 CONTINUE
+ U1_A = 0.
+ D1_A = 0.
+C
+ DUE1 = 0.
+ DDS1 = 0.
+C
+C---- similarity station pressure gradient parameter x/u du/dx
+ IBL = 2
+ BULE = 1.0
+C
+C---- set forced transition arc length position
+ CALL XIFSET(IS)
+C
+ TRAN = .FALSE.
+ TURB = .FALSE.
+C
+C**** Sweep downstream setting up BL equation linearizations
+ DO 1000 IBL=2, NBL(IS)
+C
+ IV = ISYS(IBL,IS)
+C
+ SIMI = IBL.EQ.2
+ WAKE = IBL.GT.IBLTE(IS)
+ TRAN = IBL.EQ.ITRAN(IS)
+ TURB = IBL.GT.ITRAN(IS)
+C
+ I = IPAN(IBL,IS)
+C
+C---- set primary variables for current station
+ XSI = XSSI(IBL,IS)
+ IF(IBL.LT.ITRAN(IS)) AMI = CTAU(IBL,IS)
+ IF(IBL.GE.ITRAN(IS)) CTI = CTAU(IBL,IS)
+ UEI = UEDG(IBL,IS)
+ THI = THET(IBL,IS)
+ MDI = MASS(IBL,IS)
+C
+ DSI = MDI/UEI
+C
+ IF(WAKE) THEN
+ IW = IBL - IBLTE(IS)
+ DSWAKI = WGAP(IW)
+ ELSE
+ DSWAKI = 0.
+ ENDIF
+C
+C---- set derivatives of DSI (= D2)
+ D2_M2 = 1.0/UEI
+ D2_U2 = -DSI/UEI
+C
+ DO 30 JS=1, 2
+ DO 310 JBL=2, NBL(JS)
+ J = IPAN(JBL,JS)
+ JV = ISYS(JBL,JS)
+ U2_M(JV) = -VTI(IBL,IS)*VTI(JBL,JS)*DIJ(I,J)
+ D2_M(JV) = D2_U2*U2_M(JV)
+ 310 CONTINUE
+ 30 CONTINUE
+ D2_M(IV) = D2_M(IV) + D2_M2
+C
+ U2_A = UINV_A(IBL,IS)
+ D2_A = D2_U2*U2_A
+C
+C---- "forced" changes due to mismatch between UEDG and USAV=UINV+dij*MASS
+ DUE2 = UEDG(IBL,IS) - USAV(IBL,IS)
+ DDS2 = D2_U2*DUE2
+C
+ CALL BLPRV(XSI,AMI,CTI,THI,DSI,DSWAKI,UEI)
+ CALL BLKIN
+C
+C---- check for transition and set TRAN, XT, etc. if found
+ IF(TRAN) THEN
+ CALL TRCHEK
+ AMI = AMPL2
+ ENDIF
+ IF(IBL.EQ.ITRAN(IS) .AND. .NOT.TRAN) THEN
+ WRITE(*,*) 'SETBL: Xtr??? n1 n2: ', AMPL1, AMPL2
+ ENDIF
+C
+C---- assemble 10x4 linearized system for dCtau, dTh, dDs, dUe, dXi
+C at the previous "1" station and the current "2" station
+C
+ IF(IBL.EQ.IBLTE(IS)+1) THEN
+C
+C----- define quantities at start of wake, adding TE base thickness to Dstar
+ TTE = THET(IBLTE(1),1) + THET(IBLTE(2),2)
+ DTE = DSTR(IBLTE(1),1) + DSTR(IBLTE(2),2) + ANTE
+ CTE = ( CTAU(IBLTE(1),1)*THET(IBLTE(1),1)
+ & + CTAU(IBLTE(2),2)*THET(IBLTE(2),2) ) / TTE
+ CALL TESYS(CTE,TTE,DTE)
+C
+ TTE_TTE1 = 1.0
+ TTE_TTE2 = 1.0
+ DTE_MTE1 = 1.0 / UEDG(IBLTE(1),1)
+ DTE_UTE1 = -DSTR(IBLTE(1),1) / UEDG(IBLTE(1),1)
+ DTE_MTE2 = 1.0 / UEDG(IBLTE(2),2)
+ DTE_UTE2 = -DSTR(IBLTE(2),2) / UEDG(IBLTE(2),2)
+ CTE_CTE1 = THET(IBLTE(1),1)/TTE
+ CTE_CTE2 = THET(IBLTE(2),2)/TTE
+ CTE_TTE1 = (CTAU(IBLTE(1),1) - CTE)/TTE
+ CTE_TTE2 = (CTAU(IBLTE(2),2) - CTE)/TTE
+C
+C----- re-define D1 sensitivities wrt m since D1 depends on both TE Ds values
+ DO 35 JS=1, 2
+ DO 350 JBL=2, NBL(JS)
+ J = IPAN(JBL,JS)
+ JV = ISYS(JBL,JS)
+ D1_M(JV) = DTE_UTE1*UTE1_M(JV) + DTE_UTE2*UTE2_M(JV)
+ 350 CONTINUE
+ 35 CONTINUE
+ D1_M(JVTE1) = D1_M(JVTE1) + DTE_MTE1
+ D1_M(JVTE2) = D1_M(JVTE2) + DTE_MTE2
+C
+C----- "forced" changes from UEDG --- USAV=UINV+dij*MASS mismatch
+ DUE1 = 0.
+ DDS1 = DTE_UTE1*(UEDG(IBLTE(1),1) - USAV(IBLTE(1),1))
+ & + DTE_UTE2*(UEDG(IBLTE(2),2) - USAV(IBLTE(2),2))
+C
+ ELSE
+C
+ CALL BLSYS
+C
+ ENDIF
+C
+C
+C---- Save wall shear and equil. max shear coefficient for plotting output
+ TAU(IBL,IS) = 0.5*R2*U2*U2*CF2
+ DIS(IBL,IS) = R2*U2*U2*U2*DI2*HS2*0.5
+ CTQ(IBL,IS) = CQ2
+ DELT(IBL,IS) = DE2
+ USLP(IBL,IS) = 1.60/(1.0+US2)
+C
+C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+c IF(WAKE) THEN
+c ALD = DLCON
+c ELSE
+c ALD = 1.0
+c ENDIF
+cC
+c IF(TURB .AND. .NOT.WAKE) THEN
+c GCC = GCCON
+c HKC = HK2 - 1.0 - GCC/RT2
+c IF(HKC .LT. 0.01) THEN
+c HKC = 0.01
+c ENDIF
+c ELSE
+c HKC = HK2 - 1.0
+c ENDIF
+cC
+c HR = HKC / (GACON*ALD*HK2)
+c UQ = (0.5*CF2 - HR**2) / (GBCON*D2)
+cC
+c IF(TURB) THEN
+c IBLP = MIN(IBL+1,NBL(IS))
+c IBLM = MAX(IBL-1,2 )
+c DXSSI = XSSI(IBLP,IS) - XSSI(IBLM,IS)
+c IF(DXXSI.EQ.0.0) DXSSI = 1.0
+c GUXD(IBL,IS) = -LOG(UEDG(IBLP,IS)/UEDG(IBLM,IS)) / DXSSI
+c GUXQ(IBL,IS) = -UQ
+c ELSE
+c GUXD(IBL,IS) = 0.0
+c GUXQ(IBL,IS) = 0.0
+c ENDIF
+C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+C
+C---- set XI sensitivities wrt LE Ue changes
+ IF(IS.EQ.1) THEN
+ XI_ULE1 = SST_GO
+ XI_ULE2 = -SST_GP
+ ELSE
+ XI_ULE1 = -SST_GO
+ XI_ULE2 = SST_GP
+ ENDIF
+C
+C---- stuff BL system coefficients into main Jacobian matrix
+C
+ DO 40 JV=1, NSYS
+ VM(1,JV,IV) = VS1(1,3)*D1_M(JV) + VS1(1,4)*U1_M(JV)
+ & + VS2(1,3)*D2_M(JV) + VS2(1,4)*U2_M(JV)
+ & + (VS1(1,5) + VS2(1,5) + VSX(1))
+ & *(XI_ULE1*ULE1_M(JV) + XI_ULE2*ULE2_M(JV))
+ 40 CONTINUE
+C
+ VB(1,1,IV) = VS1(1,1)
+ VB(1,2,IV) = VS1(1,2)
+C
+ VA(1,1,IV) = VS2(1,1)
+ VA(1,2,IV) = VS2(1,2)
+C
+ IF(LALFA) THEN
+ VDEL(1,2,IV) = VSR(1)*RE_CLMR + VSM(1)*MSQ_CLMR
+ ELSE
+ VDEL(1,2,IV) =
+ & (VS1(1,4)*U1_A + VS1(1,3)*D1_A)
+ & + (VS2(1,4)*U2_A + VS2(1,3)*D2_A)
+ & + (VS1(1,5) + VS2(1,5) + VSX(1))
+ & *(XI_ULE1*ULE1_A + XI_ULE2*ULE2_A)
+ ENDIF
+C
+ VDEL(1,1,IV) = VSREZ(1)
+ & + (VS1(1,4)*DUE1 + VS1(1,3)*DDS1)
+ & + (VS2(1,4)*DUE2 + VS2(1,3)*DDS2)
+ & + (VS1(1,5) + VS2(1,5) + VSX(1))
+ & *(XI_ULE1*DULE1 + XI_ULE2*DULE2)
+C
+C
+ DO 50 JV=1, NSYS
+ VM(2,JV,IV) = VS1(2,3)*D1_M(JV) + VS1(2,4)*U1_M(JV)
+ & + VS2(2,3)*D2_M(JV) + VS2(2,4)*U2_M(JV)
+ & + (VS1(2,5) + VS2(2,5) + VSX(2))
+ & *(XI_ULE1*ULE1_M(JV) + XI_ULE2*ULE2_M(JV))
+ 50 CONTINUE
+C
+ VB(2,1,IV) = VS1(2,1)
+ VB(2,2,IV) = VS1(2,2)
+C
+ VA(2,1,IV) = VS2(2,1)
+ VA(2,2,IV) = VS2(2,2)
+C
+ IF(LALFA) THEN
+ VDEL(2,2,IV) = VSR(2)*RE_CLMR + VSM(2)*MSQ_CLMR
+ ELSE
+ VDEL(2,2,IV) =
+ & (VS1(2,4)*U1_A + VS1(2,3)*D1_A)
+ & + (VS2(2,4)*U2_A + VS2(2,3)*D2_A)
+ & + (VS1(2,5) + VS2(2,5) + VSX(2))
+ & *(XI_ULE1*ULE1_A + XI_ULE2*ULE2_A)
+ ENDIF
+C
+ VDEL(2,1,IV) = VSREZ(2)
+ & + (VS1(2,4)*DUE1 + VS1(2,3)*DDS1)
+ & + (VS2(2,4)*DUE2 + VS2(2,3)*DDS2)
+ & + (VS1(2,5) + VS2(2,5) + VSX(2))
+ & *(XI_ULE1*DULE1 + XI_ULE2*DULE2)
+C
+C
+ DO 60 JV=1, NSYS
+ VM(3,JV,IV) = VS1(3,3)*D1_M(JV) + VS1(3,4)*U1_M(JV)
+ & + VS2(3,3)*D2_M(JV) + VS2(3,4)*U2_M(JV)
+ & + (VS1(3,5) + VS2(3,5) + VSX(3))
+ & *(XI_ULE1*ULE1_M(JV) + XI_ULE2*ULE2_M(JV))
+ 60 CONTINUE
+C
+ VB(3,1,IV) = VS1(3,1)
+ VB(3,2,IV) = VS1(3,2)
+C
+ VA(3,1,IV) = VS2(3,1)
+ VA(3,2,IV) = VS2(3,2)
+C
+ IF(LALFA) THEN
+ VDEL(3,2,IV) = VSR(3)*RE_CLMR + VSM(3)*MSQ_CLMR
+ ELSE
+ VDEL(3,2,IV) =
+ & (VS1(3,4)*U1_A + VS1(3,3)*D1_A)
+ & + (VS2(3,4)*U2_A + VS2(3,3)*D2_A)
+ & + (VS1(3,5) + VS2(3,5) + VSX(3))
+ & *(XI_ULE1*ULE1_A + XI_ULE2*ULE2_A)
+ ENDIF
+C
+ VDEL(3,1,IV) = VSREZ(3)
+ & + (VS1(3,4)*DUE1 + VS1(3,3)*DDS1)
+ & + (VS2(3,4)*DUE2 + VS2(3,3)*DDS2)
+ & + (VS1(3,5) + VS2(3,5) + VSX(3))
+ & *(XI_ULE1*DULE1 + XI_ULE2*DULE2)
+C
+C
+ IF(IBL.EQ.IBLTE(IS)+1) THEN
+C
+C----- redefine coefficients for TTE, DTE, etc
+ VZ(1,1) = VS1(1,1)*CTE_CTE1
+ VZ(1,2) = VS1(1,1)*CTE_TTE1 + VS1(1,2)*TTE_TTE1
+ VB(1,1,IV) = VS1(1,1)*CTE_CTE2
+ VB(1,2,IV) = VS1(1,1)*CTE_TTE2 + VS1(1,2)*TTE_TTE2
+C
+ VZ(2,1) = VS1(2,1)*CTE_CTE1
+ VZ(2,2) = VS1(2,1)*CTE_TTE1 + VS1(2,2)*TTE_TTE1
+ VB(2,1,IV) = VS1(2,1)*CTE_CTE2
+ VB(2,2,IV) = VS1(2,1)*CTE_TTE2 + VS1(2,2)*TTE_TTE2
+C
+ VZ(3,1) = VS1(3,1)*CTE_CTE1
+ VZ(3,2) = VS1(3,1)*CTE_TTE1 + VS1(3,2)*TTE_TTE1
+ VB(3,1,IV) = VS1(3,1)*CTE_CTE2
+ VB(3,2,IV) = VS1(3,1)*CTE_TTE2 + VS1(3,2)*TTE_TTE2
+C
+ ENDIF
+C
+C---- turbulent intervals will follow if currently at transition interval
+ IF(TRAN) THEN
+ TURB = .TRUE.
+C
+C------ save transition location
+ ITRAN(IS) = IBL
+ TFORCE(IS) = TRFORC
+ XSSITR(IS) = XT
+C
+C------ interpolate airfoil geometry to find transition x/c
+C- (for user output)
+ IF(IS.EQ.1) THEN
+ STR = SST - XT
+ ELSE
+ STR = SST + XT
+ ENDIF
+ CHX = XTE - XLE
+ CHY = YTE - YLE
+ CHSQ = CHX**2 + CHY**2
+ XTR = SEVAL(STR,X,XP,S,N)
+ YTR = SEVAL(STR,Y,YP,S,N)
+ XOCTR(IS) = ((XTR-XLE)*CHX + (YTR-YLE)*CHY)/CHSQ
+ YOCTR(IS) = ((YTR-YLE)*CHX - (XTR-XLE)*CHY)/CHSQ
+ ENDIF
+C
+ TRAN = .FALSE.
+C
+ IF(IBL.EQ.IBLTE(IS)) THEN
+C----- set "2" variables at TE to wake correlations for next station
+C
+ TURB = .TRUE.
+ WAKE = .TRUE.
+ CALL BLVAR(3)
+ CALL BLMID(3)
+ ENDIF
+C
+ DO 80 JS=1, 2
+ DO 810 JBL=2, NBL(JS)
+ JV = ISYS(JBL,JS)
+ U1_M(JV) = U2_M(JV)
+ D1_M(JV) = D2_M(JV)
+ 810 CONTINUE
+ 80 CONTINUE
+C
+ U1_A = U2_A
+ D1_A = D2_A
+C
+ DUE1 = DUE2
+ DDS1 = DDS2
+C
+C---- set BL variables for next station
+ DO 190 ICOM=1, NCOM
+ COM1(ICOM) = COM2(ICOM)
+ 190 CONTINUE
+C
+C---- next streamwise station
+ 1000 CONTINUE
+C
+ IF(TFORCE(IS)) THEN
+ WRITE(*,9100) IS,XOCTR(IS),ITRAN(IS)
+ 9100 FORMAT(1X,'Side',I2,' forced transition at x/c = ',F7.4,I5)
+ ELSE
+ WRITE(*,9200) IS,XOCTR(IS),ITRAN(IS)
+ 9200 FORMAT(1X,'Side',I2,' free transition at x/c = ',F7.4,I5)
+ ENDIF
+C
+C---- next airfoil side
+ 2000 CONTINUE
+C
+ RETURN
+ END
+
+
+ SUBROUTINE IBLSYS
+C---------------------------------------------
+C Sets the BL Newton system line number
+C corresponding to each BL station.
+C---------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ INCLUDE 'XBL.INC'
+C
+ IV = 0
+ DO 10 IS=1, 2
+ DO 110 IBL=2, NBL(IS)
+ IV = IV+1
+ ISYS(IBL,IS) = IV
+ 110 CONTINUE
+ 10 CONTINUE
+C
+ NSYS = IV
+ IF(NSYS.GT.2*IVX) STOP '*** IBLSYS: BL system array overflow. ***'
+C
+ RETURN
+ END
+
+
+ SUBROUTINE MRCHUE
+C----------------------------------------------------
+C Marches the BLs and wake in direct mode using
+C the UEDG array. If separation is encountered,
+C a plausible value of Hk extrapolated from
+C upstream is prescribed instead. Continuous
+C checking of transition onset is performed.
+C----------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ INCLUDE 'XBL.INC'
+ LOGICAL DIRECT
+ REAL MSQ
+C
+C---- shape parameters for separation criteria
+ HLMAX = 3.8
+ HTMAX = 2.5
+C
+ DO 2000 IS=1, 2
+C
+ WRITE(*,*) ' side ', IS, ' ...'
+C
+C---- set forced transition arc length position
+ CALL XIFSET(IS)
+C
+C---- initialize similarity station with Thwaites' formula
+ IBL = 2
+ XSI = XSSI(IBL,IS)
+ UEI = UEDG(IBL,IS)
+C BULE = LOG(UEDG(IBL+1,IS)/UEI) / LOG(XSSI(IBL+1,IS)/XSI)
+C BULE = MAX( -.08 , BULE )
+ BULE = 1.0
+ UCON = UEI/XSI**BULE
+ TSQ = 0.45/(UCON*(5.0*BULE+1.0)*REYBL) * XSI**(1.0-BULE)
+ THI = SQRT(TSQ)
+ DSI = 2.2*THI
+ AMI = 0.0
+C
+C---- initialize Ctau for first turbulent station
+ CTI = 0.03
+C
+ TRAN = .FALSE.
+ TURB = .FALSE.
+ ITRAN(IS) = IBLTE(IS)
+C
+C---- march downstream
+ DO 1000 IBL=2, NBL(IS)
+ IBM = IBL-1
+C
+ IW = IBL - IBLTE(IS)
+C
+ SIMI = IBL.EQ.2
+ WAKE = IBL.GT.IBLTE(IS)
+C
+C------ prescribed quantities
+ XSI = XSSI(IBL,IS)
+ UEI = UEDG(IBL,IS)
+C
+ IF(WAKE) THEN
+ IW = IBL - IBLTE(IS)
+ DSWAKI = WGAP(IW)
+ ELSE
+ DSWAKI = 0.
+ ENDIF
+C
+ DIRECT = .TRUE.
+C
+C------ Newton iteration loop for current station
+ DO 100 ITBL=1, 25
+C
+C-------- assemble 10x3 linearized system for dCtau, dTh, dDs, dUe, dXi
+C at the previous "1" station and the current "2" station
+C (the "1" station coefficients will be ignored)
+C
+C
+ CALL BLPRV(XSI,AMI,CTI,THI,DSI,DSWAKI,UEI)
+ CALL BLKIN
+C
+C-------- check for transition and set appropriate flags and things
+ IF((.NOT.SIMI) .AND. (.NOT.TURB)) THEN
+ CALL TRCHEK
+ AMI = AMPL2
+C
+ IF(TRAN) THEN
+ ITRAN(IS) = IBL
+ IF(CTI.LE.0.0) THEN
+ CTI = 0.03
+ S2 = CTI
+ ENDIF
+ ELSE
+ ITRAN(IS) = IBL+2
+ ENDIF
+C
+C
+ ENDIF
+C
+ IF(IBL.EQ.IBLTE(IS)+1) THEN
+ TTE = THET(IBLTE(1),1) + THET(IBLTE(2),2)
+ DTE = DSTR(IBLTE(1),1) + DSTR(IBLTE(2),2) + ANTE
+ CTE = ( CTAU(IBLTE(1),1)*THET(IBLTE(1),1)
+ & + CTAU(IBLTE(2),2)*THET(IBLTE(2),2) ) / TTE
+ CALL TESYS(CTE,TTE,DTE)
+ ELSE
+ CALL BLSYS
+ ENDIF
+C
+ IF(DIRECT) THEN
+C
+C--------- try direct mode (set dUe = 0 in currently empty 4th line)
+ VS2(4,1) = 0.
+ VS2(4,2) = 0.
+ VS2(4,3) = 0.
+ VS2(4,4) = 1.0
+ VSREZ(4) = 0.
+C
+C--------- solve Newton system for current "2" station
+ CALL GAUSS(4,4,VS2,VSREZ,1)
+C
+C--------- determine max changes and underrelax if necessary
+ DMAX = MAX( ABS(VSREZ(2)/THI),
+ & ABS(VSREZ(3)/DSI) )
+ IF(IBL.LT.ITRAN(IS)) DMAX = MAX(DMAX,ABS(VSREZ(1)/10.0))
+ IF(IBL.GE.ITRAN(IS)) DMAX = MAX(DMAX,ABS(VSREZ(1)/CTI ))
+C
+ RLX = 1.0
+ IF(DMAX.GT.0.3) RLX = 0.3/DMAX
+C
+C--------- see if direct mode is not applicable
+ IF(IBL .NE. IBLTE(IS)+1) THEN
+C
+C---------- calculate resulting kinematic shape parameter Hk
+ MSQ = UEI*UEI*HSTINV / (GM1BL*(1.0 - 0.5*UEI*UEI*HSTINV))
+ HTEST = (DSI + RLX*VSREZ(3)) / (THI + RLX*VSREZ(2))
+ CALL HKIN( HTEST, MSQ, HKTEST, DUMMY, DUMMY)
+C
+C---------- decide whether to do direct or inverse problem based on Hk
+ IF(IBL.LT.ITRAN(IS)) HMAX = HLMAX
+ IF(IBL.GE.ITRAN(IS)) HMAX = HTMAX
+ DIRECT = HKTEST.LT.HMAX
+ ENDIF
+C
+ IF(DIRECT) THEN
+C---------- update as usual
+ccc IF(IBL.LT.ITRAN(IS)) AMI = AMI + RLX*VSREZ(1)
+ IF(IBL.GE.ITRAN(IS)) CTI = CTI + RLX*VSREZ(1)
+ THI = THI + RLX*VSREZ(2)
+ DSI = DSI + RLX*VSREZ(3)
+ ELSE
+C---------- set prescribed Hk for inverse calculation at the current station
+ IF(IBL.LT.ITRAN(IS)) THEN
+C----------- laminar case: relatively slow increase in Hk downstream
+ HTARG = HK1 + 0.03*(X2-X1)/T1
+ ELSE IF(IBL.EQ.ITRAN(IS)) THEN
+C----------- transition interval: weighted laminar and turbulent case
+ HTARG = HK1 + (0.03*(XT-X1) - 0.15*(X2-XT))/T1
+ ELSE IF(WAKE) THEN
+C----------- turbulent wake case:
+C- asymptotic wake behavior with approximate Backward Euler
+ CONST = 0.03*(X2-X1)/T1
+ HK2 = HK1
+ HK2 = HK2 - (HK2 + CONST*(HK2-1.0)**3 - HK1)
+ & /(1.0 + 3.0*CONST*(HK2-1.0)**2)
+ HK2 = HK2 - (HK2 + CONST*(HK2-1.0)**3 - HK1)
+ & /(1.0 + 3.0*CONST*(HK2-1.0)**2)
+ HK2 = HK2 - (HK2 + CONST*(HK2-1.0)**3 - HK1)
+ & /(1.0 + 3.0*CONST*(HK2-1.0)**2)
+ HTARG = HK2
+ ELSE
+C----------- turbulent case: relatively fast decrease in Hk downstream
+ HTARG = HK1 - 0.15*(X2-X1)/T1
+ ENDIF
+C
+C---------- limit specified Hk to something reasonable
+ IF(WAKE) THEN
+ HTARG = MAX( HTARG , 1.01 )
+ ELSE
+ HTARG = MAX( HTARG , HMAX )
+ ENDIF
+C
+ WRITE(*,1300) IBL, HTARG
+ 1300 FORMAT(' MRCHUE: Inverse mode at', I4, ' Hk =', F8.3)
+C
+C---------- try again with prescribed Hk
+ GO TO 100
+C
+ ENDIF
+C
+ ELSE
+C
+C-------- inverse mode (force Hk to prescribed value HTARG)
+ VS2(4,1) = 0.
+ VS2(4,2) = HK2_T2
+ VS2(4,3) = HK2_D2
+ VS2(4,4) = HK2_U2
+ VSREZ(4) = HTARG - HK2
+C
+ CALL GAUSS(4,4,VS2,VSREZ,1)
+C
+C--------- added Ue clamp MD 3 Apr 03
+ DMAX = MAX( ABS(VSREZ(2)/THI),
+ & ABS(VSREZ(3)/DSI),
+ & ABS(VSREZ(4)/UEI) )
+ IF(IBL.GE.ITRAN(IS)) DMAX = MAX( DMAX , ABS(VSREZ(1)/CTI))
+C
+ RLX = 1.0
+ IF(DMAX.GT.0.3) RLX = 0.3/DMAX
+C
+C--------- update variables
+ccc IF(IBL.LT.ITRAN(IS)) AMI = AMI + RLX*VSREZ(1)
+ IF(IBL.GE.ITRAN(IS)) CTI = CTI + RLX*VSREZ(1)
+ THI = THI + RLX*VSREZ(2)
+ DSI = DSI + RLX*VSREZ(3)
+ UEI = UEI + RLX*VSREZ(4)
+C
+ ENDIF
+C
+C-------- eliminate absurd transients
+ IF(IBL.GE.ITRAN(IS)) THEN
+ CTI = MIN(CTI , 0.30 )
+ CTI = MAX(CTI , 0.0000001 )
+ ENDIF
+C
+ IF(IBL.LE.IBLTE(IS)) THEN
+ HKLIM = 1.02
+ ELSE
+ HKLIM = 1.00005
+ ENDIF
+ MSQ = UEI*UEI*HSTINV / (GM1BL*(1.0 - 0.5*UEI*UEI*HSTINV))
+ DSW = DSI - DSWAKI
+ CALL DSLIM(DSW,THI,UEI,MSQ,HKLIM)
+ DSI = DSW + DSWAKI
+C
+ IF(DMAX.LE.1.0E-5) GO TO 110
+C
+ 100 CONTINUE
+ WRITE(*,1350) IBL, IS, DMAX
+ 1350 FORMAT(' MRCHUE: Convergence failed at',I4,' side',I2,
+ & ' Res =', E12.4)
+C
+C------ the current unconverged solution might still be reasonable...
+CCC IF(DMAX .LE. 0.1) GO TO 110
+ IF(DMAX .LE. 0.1) GO TO 109
+C
+C------- the current solution is garbage --> extrapolate values instead
+ IF(IBL.GT.3) THEN
+ IF(IBL.LE.IBLTE(IS)) THEN
+ THI = THET(IBM,IS) * (XSSI(IBL,IS)/XSSI(IBM,IS))**0.5
+ DSI = DSTR(IBM,IS) * (XSSI(IBL,IS)/XSSI(IBM,IS))**0.5
+ ELSE IF(IBL.EQ.IBLTE(IS)+1) THEN
+ CTI = CTE
+ THI = TTE
+ DSI = DTE
+ ELSE
+ THI = THET(IBM,IS)
+ RATLEN = (XSSI(IBL,IS)-XSSI(IBM,IS)) / (10.0*DSTR(IBM,IS))
+ DSI = (DSTR(IBM,IS) + THI*RATLEN) / (1.0 + RATLEN)
+ ENDIF
+ IF(IBL.EQ.ITRAN(IS)) CTI = 0.05
+ IF(IBL.GT.ITRAN(IS)) CTI = CTAU(IBM,IS)
+C
+ UEI = UEDG(IBL,IS)
+ IF(IBL.GT.2 .AND. IBL.LT.NBL(IS))
+ & UEI = 0.5*(UEDG(IBL-1,IS) + UEDG(IBL+1,IS))
+ ENDIF
+C
+ 109 CALL BLPRV(XSI,AMI,CTI,THI,DSI,DSWAKI,UEI)
+ CALL BLKIN
+C
+C------- check for transition and set appropriate flags and things
+ IF((.NOT.SIMI) .AND. (.NOT.TURB)) THEN
+ CALL TRCHEK
+ AMI = AMPL2
+ IF( TRAN) ITRAN(IS) = IBL
+ IF(.NOT.TRAN) ITRAN(IS) = IBL+2
+ ENDIF
+C
+C------- set all other extrapolated values for current station
+ IF(IBL.LT.ITRAN(IS)) CALL BLVAR(1)
+ IF(IBL.GE.ITRAN(IS)) CALL BLVAR(2)
+ IF(WAKE) CALL BLVAR(3)
+C
+ IF(IBL.LT.ITRAN(IS)) CALL BLMID(1)
+ IF(IBL.GE.ITRAN(IS)) CALL BLMID(2)
+ IF(WAKE) CALL BLMID(3)
+C
+C------ pick up here after the Newton iterations
+ 110 CONTINUE
+C
+C------ store primary variables
+ IF(IBL.LT.ITRAN(IS)) CTAU(IBL,IS) = AMI
+ IF(IBL.GE.ITRAN(IS)) CTAU(IBL,IS) = CTI
+ THET(IBL,IS) = THI
+ DSTR(IBL,IS) = DSI
+ UEDG(IBL,IS) = UEI
+ MASS(IBL,IS) = DSI*UEI
+ TAU(IBL,IS) = 0.5*R2*U2*U2*CF2
+ DIS(IBL,IS) = R2*U2*U2*U2*DI2*HS2*0.5
+ CTQ(IBL,IS) = CQ2
+ DELT(IBL,IS) = DE2
+ TSTR(IBL,IS) = HS2*T2
+C
+C------ set "1" variables to "2" variables for next streamwise station
+ CALL BLPRV(XSI,AMI,CTI,THI,DSI,DSWAKI,UEI)
+ CALL BLKIN
+ DO 310 ICOM=1, NCOM
+ COM1(ICOM) = COM2(ICOM)
+ 310 CONTINUE
+C
+C------ turbulent intervals will follow transition interval or TE
+ IF(TRAN .OR. IBL.EQ.IBLTE(IS)) THEN
+ TURB = .TRUE.
+C
+C------- save transition location
+ TFORCE(IS) = TRFORC
+ XSSITR(IS) = XT
+ ENDIF
+C
+ TRAN = .FALSE.
+C
+ IF(IBL.EQ.IBLTE(IS)) THEN
+ THI = THET(IBLTE(1),1) + THET(IBLTE(2),2)
+ DSI = DSTR(IBLTE(1),1) + DSTR(IBLTE(2),2) + ANTE
+ ENDIF
+C
+ 1000 CONTINUE
+ 2000 CONTINUE
+C
+ RETURN
+ END
+
+
+ SUBROUTINE MRCHDU
+C----------------------------------------------------
+C Marches the BLs and wake in mixed mode using
+C the current Ue and Hk. The calculated Ue
+C and Hk lie along a line quasi-normal to the
+C natural Ue-Hk characteristic line of the
+C current BL so that the Goldstein or Levy-Lees
+C singularity is never encountered. Continuous
+C checking of transition onset is performed.
+C----------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ INCLUDE 'XBL.INC'
+ REAL VTMP(4,5), VZTMP(4)
+ REAL MSQ
+ccc REAL MDI
+C
+ DATA DEPS / 5.0E-6 /
+C
+C---- constant controlling how far Hk is allowed to deviate
+C- from the specified value.
+ SENSWT = 1000.0
+C
+ DO 2000 IS=1, 2
+C
+C---- set forced transition arc length position
+ CALL XIFSET(IS)
+C
+C---- set leading edge pressure gradient parameter x/u du/dx
+ IBL = 2
+ XSI = XSSI(IBL,IS)
+ UEI = UEDG(IBL,IS)
+CCC BULE = LOG(UEDG(IBL+1,IS)/UEI) / LOG(XSSI(IBL+1,IS)/XSI)
+CCC BULE = MAX( -.08 , BULE )
+ BULE = 1.0
+C
+C---- old transition station
+ ITROLD = ITRAN(IS)
+C
+ TRAN = .FALSE.
+ TURB = .FALSE.
+ ITRAN(IS) = IBLTE(IS)
+C
+C---- march downstream
+ DO 1000 IBL=2, NBL(IS)
+ IBM = IBL-1
+C
+ SIMI = IBL.EQ.2
+ WAKE = IBL.GT.IBLTE(IS)
+C
+C------ initialize current station to existing variables
+ XSI = XSSI(IBL,IS)
+ UEI = UEDG(IBL,IS)
+ THI = THET(IBL,IS)
+ DSI = DSTR(IBL,IS)
+
+CCC MDI = MASS(IBL,IS)
+C
+C------ fixed BUG MD 7 June 99
+ IF(IBL.LT.ITROLD) THEN
+ AMI = CTAU(IBL,IS)
+ CTI = 0.03
+ ELSE
+ CTI = CTAU(IBL,IS)
+ IF(CTI.LE.0.0) CTI = 0.03
+ ENDIF
+C
+CCC DSI = MDI/UEI
+C
+ IF(WAKE) THEN
+ IW = IBL - IBLTE(IS)
+ DSWAKI = WGAP(IW)
+ ELSE
+ DSWAKI = 0.
+ ENDIF
+C
+ IF(IBL.LE.IBLTE(IS)) DSI = MAX(DSI-DSWAKI,1.02000*THI) + DSWAKI
+ IF(IBL.GT.IBLTE(IS)) DSI = MAX(DSI-DSWAKI,1.00005*THI) + DSWAKI
+C
+C------ Newton iteration loop for current station
+ DO 100 ITBL=1, 25
+C
+C-------- assemble 10x3 linearized system for dCtau, dTh, dDs, dUe, dXi
+C at the previous "1" station and the current "2" station
+C (the "1" station coefficients will be ignored)
+C
+ CALL BLPRV(XSI,AMI,CTI,THI,DSI,DSWAKI,UEI)
+ CALL BLKIN
+C
+C-------- check for transition and set appropriate flags and things
+ IF((.NOT.SIMI) .AND. (.NOT.TURB)) THEN
+ CALL TRCHEK
+ AMI = AMPL2
+ IF( TRAN) ITRAN(IS) = IBL
+ IF(.NOT.TRAN) ITRAN(IS) = IBL+2
+ ENDIF
+C
+ IF(IBL.EQ.IBLTE(IS)+1) THEN
+ TTE = THET(IBLTE(1),1) + THET(IBLTE(2),2)
+ DTE = DSTR(IBLTE(1),1) + DSTR(IBLTE(2),2) + ANTE
+ CTE = ( CTAU(IBLTE(1),1)*THET(IBLTE(1),1)
+ & + CTAU(IBLTE(2),2)*THET(IBLTE(2),2) ) / TTE
+ CALL TESYS(CTE,TTE,DTE)
+ ELSE
+ CALL BLSYS
+ ENDIF
+C
+C-------- set stuff at first iteration...
+ IF(ITBL.EQ.1) THEN
+C
+C--------- set "baseline" Ue and Hk for forming Ue(Hk) relation
+ UEREF = U2
+ HKREF = HK2
+C
+C--------- if current point IBL was turbulent and is now laminar, then...
+ IF(IBL.LT.ITRAN(IS) .AND. IBL.GE.ITROLD ) THEN
+C---------- extrapolate baseline Hk
+ UEM = UEDG(IBL-1,IS)
+ DSM = DSTR(IBL-1,IS)
+ THM = THET(IBL-1,IS)
+ MSQ = UEM*UEM*HSTINV / (GM1BL*(1.0 - 0.5*UEM*UEM*HSTINV))
+ CALL HKIN( DSM/THM, MSQ, HKREF, DUMMY, DUMMY )
+ ENDIF
+C
+C--------- if current point IBL was laminar, then...
+ IF(IBL.LT.ITROLD) THEN
+C---------- reinitialize or extrapolate Ctau if it's now turbulent
+ IF(TRAN) CTAU(IBL,IS) = 0.03
+ IF(TURB) CTAU(IBL,IS) = CTAU(IBL-1,IS)
+ IF(TRAN .OR. TURB) THEN
+ CTI = CTAU(IBL,IS)
+ S2 = CTI
+ ENDIF
+ ENDIF
+C
+ ENDIF
+C
+C
+ IF(SIMI .OR. IBL.EQ.IBLTE(IS)+1) THEN
+C
+C--------- for similarity station or first wake point, prescribe Ue
+ VS2(4,1) = 0.
+ VS2(4,2) = 0.
+ VS2(4,3) = 0.
+ VS2(4,4) = U2_UEI
+ VSREZ(4) = UEREF - U2
+C
+ ELSE
+C
+C********* calculate Ue-Hk characteristic slope
+C
+ DO 20 K=1, 4
+ VZTMP(K) = VSREZ(K)
+ DO 201 L=1, 5
+ VTMP(K,L) = VS2(K,L)
+ 201 CONTINUE
+ 20 CONTINUE
+C
+C--------- set unit dHk
+ VTMP(4,1) = 0.
+ VTMP(4,2) = HK2_T2
+ VTMP(4,3) = HK2_D2
+ VTMP(4,4) = HK2_U2*U2_UEI
+ VZTMP(4) = 1.0
+C
+C--------- calculate dUe response
+ CALL GAUSS(4,4,VTMP,VZTMP,1)
+C
+C--------- set SENSWT * (normalized dUe/dHk)
+ SENNEW = SENSWT * VZTMP(4) * HKREF/UEREF
+ IF(ITBL.LE.5) THEN
+ SENS = SENNEW
+ ELSE IF(ITBL.LE.15) THEN
+ SENS = 0.5*(SENS + SENNEW)
+ ENDIF
+C
+C--------- set prescribed Ue-Hk combination
+ VS2(4,1) = 0.
+ VS2(4,2) = HK2_T2 * HKREF
+ VS2(4,3) = HK2_D2 * HKREF
+ VS2(4,4) =( HK2_U2 * HKREF + SENS/UEREF )*U2_UEI
+ VSREZ(4) = -(HKREF**2)*(HK2 / HKREF - 1.0)
+ & - SENS*(U2 / UEREF - 1.0)
+C
+ ENDIF
+C
+C-------- solve Newton system for current "2" station
+ CALL GAUSS(4,4,VS2,VSREZ,1)
+C
+C-------- determine max changes and underrelax if necessary
+C-------- (added Ue clamp MD 3 Apr 03)
+ DMAX = MAX( ABS(VSREZ(2)/THI),
+ & ABS(VSREZ(3)/DSI),
+ & ABS(VSREZ(4)/UEI) )
+ IF(IBL.GE.ITRAN(IS)) DMAX = MAX(DMAX,ABS(VSREZ(1)/(10.0*CTI)))
+C
+ RLX = 1.0
+ IF(DMAX.GT.0.3) RLX = 0.3/DMAX
+C
+C-------- update as usual
+ IF(IBL.LT.ITRAN(IS)) AMI = AMI + RLX*VSREZ(1)
+ IF(IBL.GE.ITRAN(IS)) CTI = CTI + RLX*VSREZ(1)
+ THI = THI + RLX*VSREZ(2)
+ DSI = DSI + RLX*VSREZ(3)
+ UEI = UEI + RLX*VSREZ(4)
+C
+C-------- eliminate absurd transients
+ IF(IBL.GE.ITRAN(IS)) THEN
+ CTI = MIN(CTI , 0.30 )
+ CTI = MAX(CTI , 0.0000001 )
+ ENDIF
+C
+ IF(IBL.LE.IBLTE(IS)) THEN
+ HKLIM = 1.02
+ ELSE
+ HKLIM = 1.00005
+ ENDIF
+ MSQ = UEI*UEI*HSTINV / (GM1BL*(1.0 - 0.5*UEI*UEI*HSTINV))
+ DSW = DSI - DSWAKI
+ CALL DSLIM(DSW,THI,UEI,MSQ,HKLIM)
+ DSI = DSW + DSWAKI
+C
+ IF(DMAX.LE.DEPS) GO TO 110
+C
+ 100 CONTINUE
+C
+ WRITE(*,1350) IBL, IS, DMAX
+ 1350 FORMAT(' MRCHDU: Convergence failed at',I4,' side',I2,
+ & ' Res =', E12.4)
+C
+C------ the current unconverged solution might still be reasonable...
+CCC IF(DMAX .LE. 0.1) GO TO 110
+ IF(DMAX .LE. 0.1) GO TO 109
+C
+C------- the current solution is garbage --> extrapolate values instead
+ IF(IBL.GT.3) THEN
+ IF(IBL.LE.IBLTE(IS)) THEN
+ THI = THET(IBM,IS) * (XSSI(IBL,IS)/XSSI(IBM,IS))**0.5
+ DSI = DSTR(IBM,IS) * (XSSI(IBL,IS)/XSSI(IBM,IS))**0.5
+ UEI = UEDG(IBM,IS)
+ ELSE IF(IBL.EQ.IBLTE(IS)+1) THEN
+ CTI = CTE
+ THI = TTE
+ DSI = DTE
+ UEI = UEDG(IBM,IS)
+ ELSE
+ THI = THET(IBM,IS)
+ RATLEN = (XSSI(IBL,IS)-XSSI(IBM,IS)) / (10.0*DSTR(IBM,IS))
+ DSI = (DSTR(IBM,IS) + THI*RATLEN) / (1.0 + RATLEN)
+ UEI = UEDG(IBM,IS)
+ ENDIF
+ IF(IBL.EQ.ITRAN(IS)) CTI = 0.05
+ IF(IBL.GT.ITRAN(IS)) CTI = CTAU(IBM,IS)
+ ENDIF
+C
+ 109 CALL BLPRV(XSI,AMI,CTI,THI,DSI,DSWAKI,UEI)
+ CALL BLKIN
+C
+C------- check for transition and set appropriate flags and things
+ IF((.NOT.SIMI) .AND. (.NOT.TURB)) THEN
+ CALL TRCHEK
+ AMI = AMPL2
+ IF( TRAN) ITRAN(IS) = IBL
+ IF(.NOT.TRAN) ITRAN(IS) = IBL+2
+ ENDIF
+C
+C------- set all other extrapolated values for current station
+ IF(IBL.LT.ITRAN(IS)) CALL BLVAR(1)
+ IF(IBL.GE.ITRAN(IS)) CALL BLVAR(2)
+ IF(WAKE) CALL BLVAR(3)
+C
+ IF(IBL.LT.ITRAN(IS)) CALL BLMID(1)
+ IF(IBL.GE.ITRAN(IS)) CALL BLMID(2)
+ IF(WAKE) CALL BLMID(3)
+C
+C------ pick up here after the Newton iterations
+ 110 CONTINUE
+C
+ SENS = SENNEW
+C
+C------ store primary variables
+ IF(IBL.LT.ITRAN(IS)) CTAU(IBL,IS) = AMI
+ IF(IBL.GE.ITRAN(IS)) CTAU(IBL,IS) = CTI
+ THET(IBL,IS) = THI
+ DSTR(IBL,IS) = DSI
+ UEDG(IBL,IS) = UEI
+ MASS(IBL,IS) = DSI*UEI
+ TAU(IBL,IS) = 0.5*R2*U2*U2*CF2
+ DIS(IBL,IS) = R2*U2*U2*U2*DI2*HS2*0.5
+ CTQ(IBL,IS) = CQ2
+ DELT(IBL,IS) = DE2
+ TSTR(IBL,IS) = HS2*T2
+C
+C------ set "1" variables to "2" variables for next streamwise station
+ CALL BLPRV(XSI,AMI,CTI,THI,DSI,DSWAKI,UEI)
+ CALL BLKIN
+ DO 310 ICOM=1, NCOM
+ COM1(ICOM) = COM2(ICOM)
+ 310 CONTINUE
+C
+C
+C------ turbulent intervals will follow transition interval or TE
+ IF(TRAN .OR. IBL.EQ.IBLTE(IS)) THEN
+ TURB = .TRUE.
+C
+C------- save transition location
+ TFORCE(IS) = TRFORC
+ XSSITR(IS) = XT
+ ENDIF
+C
+ TRAN = .FALSE.
+C
+ 1000 CONTINUE
+C
+ 2000 CONTINUE
+C
+ RETURN
+ END
+
+
+ SUBROUTINE XIFSET(IS)
+C-----------------------------------------------------
+C Sets forced-transition BL coordinate locations.
+C-----------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ INCLUDE 'XBL.INC'
+C
+ IF(XSTRIP(IS).GE.1.0) THEN
+ XIFORC = XSSI(IBLTE(IS),IS)
+ RETURN
+ ENDIF
+C
+ CHX = XTE - XLE
+ CHY = YTE - YLE
+ CHSQ = CHX**2 + CHY**2
+C
+C---- calculate chord-based x/c, y/c
+ DO 10 I=1, N
+ W1(I) = ((X(I)-XLE)*CHX + (Y(I)-YLE)*CHY) / CHSQ
+ W2(I) = ((Y(I)-YLE)*CHX - (X(I)-XLE)*CHY) / CHSQ
+ 10 CONTINUE
+C
+ CALL SPLIND(W1,W3,S,N,-999.0,-999.0)
+ CALL SPLIND(W2,W4,S,N,-999.0,-999.0)
+C
+ IF(IS.EQ.1) THEN
+C
+C----- set approximate arc length of forced transition point for SINVRT
+ STR = SLE + (S(1)-SLE)*XSTRIP(IS)
+C
+C----- calculate actual arc length
+ CALL SINVRT(STR,XSTRIP(IS),W1,W3,S,N)
+C
+C----- set BL coordinate value
+ XIFORC = MIN( (SST - STR) , XSSI(IBLTE(IS),IS) )
+C
+ ELSE
+C----- same for bottom side
+C
+ STR = SLE + (S(N)-SLE)*XSTRIP(IS)
+ CALL SINVRT(STR,XSTRIP(IS),W1,W3,S,N)
+ XIFORC = MIN( (STR - SST) , XSSI(IBLTE(IS),IS) )
+C
+ ENDIF
+C
+ IF(XIFORC .LT. 0.0) THEN
+ WRITE(*,1000) IS
+ 1000 FORMAT(/' *** Stagnation point is past trip on side',I2,' ***')
+ XIFORC = XSSI(IBLTE(IS),IS)
+ ENDIF
+C
+ RETURN
+ END
+
+
+
+
+ SUBROUTINE UPDATE
+C------------------------------------------------------------------
+C Adds on Newton deltas to boundary layer variables.
+C Checks for excessive changes and underrelaxes if necessary.
+C Calculates max and rms changes.
+C Also calculates the change in the global variable "AC".
+C If LALFA=.TRUE. , "AC" is CL
+C If LALFA=.FALSE., "AC" is alpha
+C------------------------------------------------------------------
+ INCLUDE 'XFOIL.INC'
+ REAL UNEW(IVX,2), U_AC(IVX,2)
+ REAL QNEW(IQX), Q_AC(IQX)
+ EQUIVALENCE (VA(1,1,1), UNEW(1,1)) ,
+ & (VB(1,1,1), QNEW(1) )
+ EQUIVALENCE (VA(1,1,IVX), U_AC(1,1)) ,
+ & (VB(1,1,IVX), Q_AC(1) )
+ REAL MSQ
+C
+C---- max allowable alpha changes per iteration
+ DALMAX = 0.5*DTOR
+ DALMIN = -0.5*DTOR
+C
+C---- max allowable CL change per iteration
+ DCLMAX = 0.5
+ DCLMIN = -0.5
+ IF(MATYP.NE.1) DCLMIN = MAX(-0.5 , -0.9*CL)
+C
+ HSTINV = GAMM1*(MINF/QINF)**2 / (1.0 + 0.5*GAMM1*MINF**2)
+C
+C---- calculate new Ue distribution assuming no under-relaxation
+C- also set the sensitivity of Ue wrt to alpha or Re
+ DO 1 IS=1, 2
+ DO 10 IBL=2, NBL(IS)
+ I = IPAN(IBL,IS)
+C
+ DUI = 0.
+ DUI_AC = 0.
+ DO 100 JS=1, 2
+ DO 1000 JBL=2, NBL(JS)
+ J = IPAN(JBL,JS)
+ JV = ISYS(JBL,JS)
+ UE_M = -VTI(IBL,IS)*VTI(JBL,JS)*DIJ(I,J)
+ DUI = DUI + UE_M*(MASS(JBL,JS)+VDEL(3,1,JV))
+ DUI_AC = DUI_AC + UE_M*( -VDEL(3,2,JV))
+ 1000 CONTINUE
+ 100 CONTINUE
+C
+C-------- UINV depends on "AC" only if "AC" is alpha
+ IF(LALFA) THEN
+ UINV_AC = 0.
+ ELSE
+ UINV_AC = UINV_A(IBL,IS)
+ ENDIF
+C
+ UNEW(IBL,IS) = UINV(IBL,IS) + DUI
+ U_AC(IBL,IS) = UINV_AC + DUI_AC
+C
+ 10 CONTINUE
+ 1 CONTINUE
+C
+C---- set new Qtan from new Ue with appropriate sign change
+ DO 2 IS=1, 2
+ DO 20 IBL=2, IBLTE(IS)
+ I = IPAN(IBL,IS)
+ QNEW(I) = VTI(IBL,IS)*UNEW(IBL,IS)
+ Q_AC(I) = VTI(IBL,IS)*U_AC(IBL,IS)
+ 20 CONTINUE
+ 2 CONTINUE
+C
+C---- calculate new CL from this new Qtan
+ SA = SIN(ALFA)
+ CA = COS(ALFA)
+C
+ BETA = SQRT(1.0 - MINF**2)
+ BETA_MSQ = -0.5/BETA
+C
+ BFAC = 0.5*MINF**2 / (1.0 + BETA)
+ BFAC_MSQ = 0.5 / (1.0 + BETA)
+ & - BFAC / (1.0 + BETA) * BETA_MSQ
+C
+ CLNEW = 0.
+ CL_A = 0.
+ CL_MS = 0.
+ CL_AC = 0.
+C
+ I = 1
+ CGINC = 1.0 - (QNEW(I)/QINF)**2
+ CPG1 = CGINC / (BETA + BFAC*CGINC)
+ CPG1_MS = -CPG1/(BETA + BFAC*CGINC)*(BETA_MSQ + BFAC_MSQ*CGINC)
+C
+ CPI_Q = -2.0*QNEW(I)/QINF**2
+ CPC_CPI = (1.0 - BFAC*CPG1)/ (BETA + BFAC*CGINC)
+ CPG1_AC = CPC_CPI*CPI_Q*Q_AC(I)
+C
+ DO 3 I=1, N
+ IP = I+1
+ IF(I.EQ.N) IP = 1
+C
+ CGINC = 1.0 - (QNEW(IP)/QINF)**2
+ CPG2 = CGINC / (BETA + BFAC*CGINC)
+ CPG2_MS = -CPG2/(BETA + BFAC*CGINC)*(BETA_MSQ + BFAC_MSQ*CGINC)
+C
+ CPI_Q = -2.0*QNEW(IP)/QINF**2
+ CPC_CPI = (1.0 - BFAC*CPG2)/ (BETA + BFAC*CGINC)
+ CPG2_AC = CPC_CPI*CPI_Q*Q_AC(IP)
+C
+ DX = (X(IP) - X(I))*CA + (Y(IP) - Y(I))*SA
+ DX_A = -(X(IP) - X(I))*SA + (Y(IP) - Y(I))*CA
+C
+ AG = 0.5*(CPG2 + CPG1 )
+ AG_MS = 0.5*(CPG2_MS + CPG1_MS)
+ AG_AC = 0.5*(CPG2_AC + CPG1_AC)
+C
+ CLNEW = CLNEW + DX *AG
+ CL_A = CL_A + DX_A*AG
+ CL_MS = CL_MS + DX *AG_MS
+ CL_AC = CL_AC + DX *AG_AC
+C
+ CPG1 = CPG2
+ CPG1_MS = CPG2_MS
+ CPG1_AC = CPG2_AC
+ 3 CONTINUE
+C
+C---- initialize under-relaxation factor
+ RLX = 1.0
+C
+ IF(LALFA) THEN
+C===== alpha is prescribed: AC is CL
+C
+C----- set change in Re to account for CL changing, since Re = Re(CL)
+ DAC = (CLNEW - CL) / (1.0 - CL_AC - CL_MS*2.0*MINF*MINF_CL)
+C
+C----- set under-relaxation factor if Re change is too large
+ IF(RLX*DAC .GT. DCLMAX) RLX = DCLMAX/DAC
+ IF(RLX*DAC .LT. DCLMIN) RLX = DCLMIN/DAC
+C
+ ELSE
+C===== CL is prescribed: AC is alpha
+C
+C----- set change in alpha to drive CL to prescribed value
+ DAC = (CLNEW - CLSPEC) / (0.0 - CL_AC - CL_A)
+C
+C----- set under-relaxation factor if alpha change is too large
+ IF(RLX*DAC .GT. DALMAX) RLX = DALMAX/DAC
+ IF(RLX*DAC .LT. DALMIN) RLX = DALMIN/DAC
+C
+ ENDIF
+C
+ RMSBL = 0.
+ RMXBL = 0.
+C
+ DHI = 1.5
+ DLO = -.5
+C
+C---- calculate changes in BL variables and under-relaxation if needed
+ DO 4 IS=1, 2
+ DO 40 IBL=2, NBL(IS)
+ IV = ISYS(IBL,IS)
+C
+
+
+C-------- set changes without underrelaxation
+ DCTAU = VDEL(1,1,IV) - DAC*VDEL(1,2,IV)
+ DTHET = VDEL(2,1,IV) - DAC*VDEL(2,2,IV)
+ DMASS = VDEL(3,1,IV) - DAC*VDEL(3,2,IV)
+ DUEDG = UNEW(IBL,IS) + DAC*U_AC(IBL,IS) - UEDG(IBL,IS)
+ DDSTR = (DMASS - DSTR(IBL,IS)*DUEDG)/UEDG(IBL,IS)
+C
+C-------- normalize changes
+ IF(IBL.LT.ITRAN(IS)) DN1 = DCTAU / 10.0
+ IF(IBL.GE.ITRAN(IS)) DN1 = DCTAU / CTAU(IBL,IS)
+ DN2 = DTHET / THET(IBL,IS)
+ DN3 = DDSTR / DSTR(IBL,IS)
+ DN4 = ABS(DUEDG)/0.25
+C
+C-------- accumulate for rms change
+ RMSBL = RMSBL + DN1**2 + DN2**2 + DN3**2 + DN4**2
+C
+C-------- see if Ctau needs underrelaxation
+ RDN1 = RLX*DN1
+ IF(ABS(DN1) .GT. ABS(RMXBL)) THEN
+ RMXBL = DN1
+ IF(IBL.LT.ITRAN(IS)) VMXBL = 'n'
+ IF(IBL.GE.ITRAN(IS)) VMXBL = 'C'
+ IMXBL = IBL
+ ISMXBL = IS
+ ENDIF
+ IF(RDN1 .GT. DHI) RLX = DHI/DN1
+ IF(RDN1 .LT. DLO) RLX = DLO/DN1
+C
+C-------- see if Theta needs underrelaxation
+ RDN2 = RLX*DN2
+ IF(ABS(DN2) .GT. ABS(RMXBL)) THEN
+ RMXBL = DN2
+ VMXBL = 'T'
+ IMXBL = IBL
+ ISMXBL = IS
+ ENDIF
+ IF(RDN2 .GT. DHI) RLX = DHI/DN2
+ IF(RDN2 .LT. DLO) RLX = DLO/DN2
+C
+C-------- see if Dstar needs underrelaxation
+ RDN3 = RLX*DN3
+ IF(ABS(DN3) .GT. ABS(RMXBL)) THEN
+ RMXBL = DN3
+ VMXBL = 'D'
+ IMXBL = IBL
+ ISMXBL = IS
+ ENDIF
+ IF(RDN3 .GT. DHI) RLX = DHI/DN3
+ IF(RDN3 .LT. DLO) RLX = DLO/DN3
+C
+C-------- see if Ue needs underrelaxation
+ RDN4 = RLX*DN4
+ IF(ABS(DN4) .GT. ABS(RMXBL)) THEN
+ RMXBL = DUEDG
+ VMXBL = 'U'
+ IMXBL = IBL
+ ISMXBL = IS
+ ENDIF
+ IF(RDN4 .GT. DHI) RLX = DHI/DN4
+ IF(RDN4 .LT. DLO) RLX = DLO/DN4
+C
+ 40 CONTINUE
+ 4 CONTINUE
+C
+C---- set true rms change
+ RMSBL = SQRT( RMSBL / (4.0*FLOAT( NBL(1)+NBL(2) )) )
+C
+C
+ IF(LALFA) THEN
+C----- set underrelaxed change in Reynolds number from change in lift
+ CL = CL + RLX*DAC
+ ELSE
+C----- set underrelaxed change in alpha
+ ALFA = ALFA + RLX*DAC
+ ADEG = ALFA/DTOR
+ ENDIF
+C
+C---- update BL variables with underrelaxed changes
+ DO 5 IS=1, 2
+ DO 50 IBL=2, NBL(IS)
+ IV = ISYS(IBL,IS)
+C
+ DCTAU = VDEL(1,1,IV) - DAC*VDEL(1,2,IV)
+ DTHET = VDEL(2,1,IV) - DAC*VDEL(2,2,IV)
+ DMASS = VDEL(3,1,IV) - DAC*VDEL(3,2,IV)
+ DUEDG = UNEW(IBL,IS) + DAC*U_AC(IBL,IS) - UEDG(IBL,IS)
+ DDSTR = (DMASS - DSTR(IBL,IS)*DUEDG)/UEDG(IBL,IS)
+C
+ CTAU(IBL,IS) = CTAU(IBL,IS) + RLX*DCTAU
+ THET(IBL,IS) = THET(IBL,IS) + RLX*DTHET
+ DSTR(IBL,IS) = DSTR(IBL,IS) + RLX*DDSTR
+ UEDG(IBL,IS) = UEDG(IBL,IS) + RLX*DUEDG
+C
+ IF(IBL.GT.IBLTE(IS)) THEN
+ IW = IBL - IBLTE(IS)
+ DSWAKI = WGAP(IW)
+ ELSE
+ DSWAKI = 0.
+ ENDIF
+C
+C-------- eliminate absurd transients
+ IF(IBL.GE.ITRAN(IS))
+ & CTAU(IBL,IS) = MIN( CTAU(IBL,IS) , 0.25 )
+C
+ IF(IBL.LE.IBLTE(IS)) THEN
+ HKLIM = 1.02
+ ELSE
+ HKLIM = 1.00005
+ ENDIF
+ MSQ = UEDG(IBL,IS)**2*HSTINV
+ & / (GAMM1*(1.0 - 0.5*UEDG(IBL,IS)**2*HSTINV))
+ DSW = DSTR(IBL,IS) - DSWAKI
+ CALL DSLIM(DSW,THET(IBL,IS),UEDG(IBL,IS),MSQ,HKLIM)
+ DSTR(IBL,IS) = DSW + DSWAKI
+C
+C-------- set new mass defect (nonlinear update)
+ MASS(IBL,IS) = DSTR(IBL,IS) * UEDG(IBL,IS)
+C
+ 50 CONTINUE
+C
+C------ make sure there are no "islands" of negative Ue
+ DO IBL = 3, IBLTE(IS)
+ IF(UEDG(IBL-1,IS) .GT. 0.0 .AND.
+ & UEDG(IBL ,IS) .LE. 0.0 ) THEN
+ UEDG(IBL,IS) = UEDG(IBL-1,IS)
+ MASS(IBL,IS) = DSTR(IBL,IS) * UEDG(IBL,IS)
+ ENDIF
+ ENDDO
+ 5 CONTINUE
+C
+C
+C---- equate upper wake arrays to lower wake arrays
+ DO 6 KBL=1, NBL(2)-IBLTE(2)
+ CTAU(IBLTE(1)+KBL,1) = CTAU(IBLTE(2)+KBL,2)
+ THET(IBLTE(1)+KBL,1) = THET(IBLTE(2)+KBL,2)
+ DSTR(IBLTE(1)+KBL,1) = DSTR(IBLTE(2)+KBL,2)
+ UEDG(IBLTE(1)+KBL,1) = UEDG(IBLTE(2)+KBL,2)
+ TAU(IBLTE(1)+KBL,1) = TAU(IBLTE(2)+KBL,2)
+ DIS(IBLTE(1)+KBL,1) = DIS(IBLTE(2)+KBL,2)
+ CTQ(IBLTE(1)+KBL,1) = CTQ(IBLTE(2)+KBL,2)
+ DELT(IBLTE(1)+KBL,1) = DELT(IBLTE(2)+KBL,2)
+ TSTR(IBLTE(1)+KBL,1) = TSTR(IBLTE(2)+KBL,2)
+ 6 CONTINUE
+C
+ RETURN
+ END
+
+
+
+ SUBROUTINE DSLIM(DSTR,THET,UEDG,MSQ,HKLIM)
+ IMPLICIT REAL (A-H,M,O-Z)
+C
+ H = DSTR/THET
+ CALL HKIN(H,MSQ,HK,HK_H,HK_M)
+C
+ DH = MAX( 0.0 , HKLIM-HK ) / HK_H
+ DSTR = DSTR + DH*THET
+C
+ RETURN
+ END
+
+
+
+ SUBROUTINE BLPINI
+ INCLUDE 'BLPAR.INC'
+C
+ SCCON = 5.6
+ GACON = 6.70
+ GBCON = 0.75
+ GCCON = 18.0
+ DLCON = 0.9
+C
+ CTRCON = 1.8
+ CTRCEX = 3.3
+C
+ DUXCON = 1.0
+C
+ CTCON = 0.5/(GACON**2 * GBCON)
+C
+ CFFAC = 1.0
+C
+ RETURN
+ END
+