aboutsummaryrefslogtreecommitdiff
path: root/src/iopol.f
diff options
context:
space:
mode:
Diffstat (limited to 'src/iopol.f')
-rw-r--r--src/iopol.f746
1 files changed, 746 insertions, 0 deletions
diff --git a/src/iopol.f b/src/iopol.f
new file mode 100644
index 0000000..6cea1a9
--- /dev/null
+++ b/src/iopol.f
@@ -0,0 +1,746 @@
+C***********************************************************************
+C Module: iopol.f
+C
+C Copyright (C) 2000 Mark Drela, Harold Youngren
+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***********************************************************************
+
+ SUBROUTINE POLREAD(LU,FNPOL,ERROR,
+ & NAX,NA,CPOL,
+ & REYN1,MACH1,ACRIT,XTRIP,
+ & PTRAT,ETAP,
+ & NAME, IRETYP,IMATYP,
+ & ISX,NBL,CPOLSD,
+ & CODE,VERSION )
+ INCLUDE 'PINDEX.INC'
+ CHARACTER*(*) FNPOL, NAME
+ LOGICAL ERROR
+ CHARACTER*(*) CODE
+ REAL CPOL(NAX,IPTOT), CPOLSD(NAX,ISX,JPTOT)
+ REAL MACH1, XTRIP(ISX)
+C--------------------------------------------------------
+C Reads in polar save file
+C
+C Input:
+C LU logical unit to use for reading
+C FNPOL name of polar file to be read,
+C if FNPOL(1:1).eq.' ', unit LU will be read
+C if it is already open
+C NAX polar point array dimension
+C ISX airfoil side array dimension
+C
+C Output:
+C ERROR T if a READ error occurred
+C NA number polar points
+C CPOL polar coefficients and parameters
+C REYN1 Reynolds number for CL=1
+C MACH1 Mach number for CL=1
+C ACRIT Critical amplification ratio
+C XTRIP Trip locations
+C PTRAT Actuator disk total-pressure ratio
+C ETAP Actuator disk thermal efficiency
+C NAME airfoil name string
+C IRETYP flag giving type of Re variation with CL
+C IMATYP flag giving type of Ma variation with CL
+C NBL number of airfoil elements
+C CPOLSD airfoil side-related parameters
+C CODE code used to compute polar
+C VERSION code version
+C--------------------------------------------------------
+ CHARACTER*128 LINE
+ CHARACTER*1 DUMMY
+ REAL RINP(0:IPTOT+2*JPTOT)
+C
+ INTEGER IPOL(IPTOT), ISPOL(2,JPTOT)
+ INTEGER ITMP(IPTOT+2*JPTOT), ITMP0(IPTOT+2*JPTOT)
+ LOGICAL LOPEN, LHEAD, LDLAB
+ LOGICAL LIRE, LIMA, LINC, LJTP
+ CHARACTER*20 CPNAME
+C
+C
+ ERROR = .FALSE.
+ LHEAD = .TRUE.
+C
+ NA = 0
+ NBL = 1
+C
+c KCH = 0
+c KMC = 0
+C
+ NIPOL = 0
+ DO IP = 1, IPTOT
+ IPOL(IP) = 0
+ ENDDO
+ DO JP = 1, JPTOT
+ ISPOL(1,JP) = 0
+ ISPOL(2,JP) = 0
+ ENDDO
+C
+C---- assume Re,Mach will not be given in header
+ IRETYP = 0
+ IMATYP = 0
+C
+C---- do we have to open the file?
+ LOPEN = FNPOL .NE. ' '
+C
+ IF(LOPEN) OPEN(LU,FILE=FNPOL,STATUS='OLD',ERR=90)
+C
+C=============================================================
+C---- start data reading loop
+ 500 CONTINUE
+ READ(LU,1000,END=80) LINE
+ IF(LINE.EQ.' ') GO TO 500
+C
+ IF(LHEAD) THEN
+C----- parse to get header info
+C
+C----- assume this will be the data-label line
+ LDLAB = .TRUE.
+C
+C--------------------------------------------
+ K = INDEX(LINE,'Version')
+ IF(K.NE.0) THEN
+C------ code,version line
+ DO K1=1, 128
+ IF(LINE(K1:K1).NE.' ') GO TO 10
+ ENDDO
+C
+ 10 CONTINUE
+ IF(K.GT.K1) THEN
+ CODE = LINE(K1:K-1)
+ READ(LINE(K+7:128),*,ERR=11) VERSION
+ ENDIF
+ 11 CONTINUE
+ LDLAB = .FALSE.
+ ENDIF
+C
+C--------------------------------------------
+ KF = INDEX(LINE,'for:')
+ IF(KF.NE.0) THEN
+C------ airfoil name line
+ NAME = LINE(KF+5:128)
+ LDLAB = .FALSE.
+ ENDIF
+C
+C--------------------------------------------
+ KE = INDEX(LINE,'elements')
+ IF(KE.GT.0) THEN
+C------ element-number line
+ READ(LINE(KE-4:KE-1),*,ERR=60) NBL
+C------ truncate name line to eliminate elements #
+ NAME = LINE(KF+5:KE-4)
+C
+ IF(2*NBL .GT. ISX) THEN
+ NBL = ISX/2
+ WRITE(*,*)
+ & 'POLREAD: Number of elements set to array limit', NBL
+ ENDIF
+ LDLAB = .FALSE.
+ ENDIF
+C
+C--------------------------------------------
+ KR = INDEX(LINE,'Reynolds number')
+ KM = INDEX(LINE,'Mach number')
+C
+ IF(KR.NE.0) THEN
+C------ Re-type line
+ IF(KM.GT.KR) THEN
+ KEND = KM-1
+ ELSE
+ KEND = 128
+ ENDIF
+ IF (INDEX(LINE(KR:KEND),'fixed').NE.0) THEN
+ IRETYP = 1
+ ELSEIF(INDEX(LINE(KR:KEND),'1/sqrt(CL)').NE.0) THEN
+ IRETYP = 2
+ ELSEIF(INDEX(LINE(KR:KEND),'1/CL').NE.0) THEN
+ IRETYP = 3
+ ENDIF
+ LDLAB = .FALSE.
+ ENDIF
+C
+ IF(KM.NE.0) THEN
+C------ Ma-type line
+ IF(KR.GT.KM) THEN
+ KEND = KR-1
+ ELSE
+ KEND = 128
+ ENDIF
+ IF (INDEX(LINE(KM:KEND),'fixed').NE.0) THEN
+ IMATYP = 1
+ ELSEIF(INDEX(LINE(KM:KEND),'1/sqrt(CL)').NE.0) THEN
+ IMATYP = 2
+ ELSEIF(INDEX(LINE(KM:KEND),'1/CL').NE.0) THEN
+ IMATYP = 3
+ ENDIF
+ LDLAB = .FALSE.
+ ENDIF
+C
+C--------------------------------------------
+C---- find specified BL trip location
+ K = INDEX(LINE,'xtrf')
+ IF(K.NE.0) THEN
+C------ new style xtrip line
+ KT = INDEX(LINE,'(top)')
+ KB = INDEX(LINE,'(bottom)')
+ KE = INDEX(LINE,'element ')
+C--- check for old style trip line
+ KS = INDEX(LINE,'(suc')
+ KP = INDEX(LINE,'(pre')
+C
+ IF(KE.NE.0) THEN
+ READ(LINE(KE+7:KE+12),*,ERR=21) N
+ ELSE
+ N = 1
+ ENDIF
+ IF(N.LE.NBL) THEN
+ IS1 = 2*N-1
+ IS2 = 2*N
+ XTRIP(IS1) = 1.0
+ XTRIP(IS2) = 1.0
+ IF(KT.GT.0) READ(LINE(K+6:KT-1) ,*,ERR=21) XTRIP(IS1)
+ IF(KB.GT.KT) READ(LINE(KT+5:KB-1),*,ERR=21) XTRIP(IS2)
+ IF(KS.GT.0) READ(LINE(K+6:KS-1) ,*,ERR=21) XTRIP(IS1)
+ IF(KP.GT.KS) READ(LINE(KS+5:KP-1),*,ERR=21) XTRIP(IS2)
+ ENDIF
+ 21 CONTINUE
+ LDLAB = .FALSE.
+ ENDIF
+C
+C--------------------------------------------
+ K = INDEX(LINE,'Mach =')
+ IF(K.NE.0) THEN
+ READ(LINE(K+6:128),*,ERR=31) MACH1
+ 31 CONTINUE
+ LDLAB = .FALSE.
+ ENDIF
+C
+C--------------------------------------------
+ K = INDEX(LINE,'Re =')
+ IF(K.NE.0) THEN
+ READ(LINE(K+4:128),*,ERR=32) REYN1
+ REYN1 = REYN1 * 1.0E6
+ 32 CONTINUE
+ LDLAB = .FALSE.
+ ENDIF
+C
+C--------------------------------------------
+ K = INDEX(LINE,'Ncrit =')
+ IF(K.NE.0) THEN
+ READ(LINE(K+7:128),*,ERR=33) ACRIT
+ 33 CONTINUE
+ LDLAB = .FALSE.
+ ENDIF
+C
+C--------------------------------------------
+ K = INDEX(LINE,'pi_p =')
+ IF(K.NE.0) THEN
+ READ(LINE(K+6:128),*,ERR=34) PTRAT
+ 34 CONTINUE
+ LDLAB = .FALSE.
+ ENDIF
+C
+C--------------------------------------------
+ K = INDEX(LINE,'eta_p =')
+ IF(K.NE.0) THEN
+ READ(LINE(K+7:128),*,ERR=35) ETAP
+ 35 CONTINUE
+ LDLAB = .FALSE.
+ ENDIF
+C
+C--------------------------------------------
+ IF(LDLAB .AND. NIPOL.EQ.0) THEN
+C------ process line for possible data labels
+ DO IP = 1, IPTOT
+ CALL STRIP(CPOLNAME(IP),NNAME)
+C
+C-------- mark this parameter for reading
+ K = INDEX(LINE,CPOLNAME(IP)(1:NNAME))
+ ITMP0(IP) = K
+ ITMP(IP) = K
+ ENDDO
+C
+ DO JP = 1, JPTOT
+ CALL STRIP(CPOLSNAME(JP),NNAME)
+C
+ CPNAME = 'Top ' // CPOLSNAME(JP)
+ K1 = INDEX(LINE,CPNAME(1:NNAME+4))
+ CPNAME = 'Top_' // CPOLSNAME(JP)
+ K2 = INDEX(LINE,CPNAME(1:NNAME+4))
+ ITMP0(IPTOT+JP) = MAX(K1,K2)
+ ITMP (IPTOT+JP) = MAX(K1,K2)
+C
+ CPNAME = 'Bot ' // CPOLSNAME(JP)
+ K1 = INDEX(LINE,CPNAME(1:NNAME+4))
+ CPNAME = 'Bot_' // CPOLSNAME(JP)
+ K2 = INDEX(LINE,CPNAME(1:NNAME+4))
+ ITMP0(IPTOT+JP+JPTOT) = MAX(K1,K2)
+ ITMP (IPTOT+JP+JPTOT) = MAX(K1,K2)
+ ENDDO
+C
+C------ bubble-sort data label positions in line string
+ DO IPASS = 1, IPTOT+2*JPTOT
+ DO IP = 1, IPTOT+2*JPTOT-1
+ IF(ITMP(IP).GT.ITMP(IP+1)) THEN
+ ITMPP1 = ITMP(IP+1)
+ ITMP(IP+1) = ITMP(IP)
+ ITMP(IP) = ITMPP1
+ ENDIF
+ ENDDO
+ ENDDO
+C
+C------ assign data position to each parameter
+ DO IPT = 1, IPTOT+2*JPTOT
+ IF(ITMP(IPT).GT.0) THEN
+ NIPOL = NIPOL + 1
+ DO IP = 1, IPTOT
+ IF(ITMP(IPT).EQ.ITMP0(IP)) IPOL(IP) = NIPOL
+ ENDDO
+ DO JP = 1, JPTOT
+ IF(ITMP(IPT).EQ.ITMP0(IPTOT+JP )) ISPOL(1,JP) = NIPOL
+ IF(ITMP(IPT).EQ.ITMP0(IPTOT+JPTOT+JP)) ISPOL(2,JP) = NIPOL
+ ENDDO
+ ENDIF
+ ENDDO
+C
+ ENDIF
+C
+C--------------------------------------------
+ IF(INDEX(LINE,'-----').NE.0) THEN
+ LHEAD = .FALSE.
+ ENDIF
+C
+C--------------------------------------------------------------
+ ELSE
+C----- read polar data lines
+ IA = NA + 1
+C
+ NINP = IPTOT+2*JPTOT
+ CALL GETFLT(LINE,RINP(1),NINP,ERROR)
+ IF(ERROR) GO TO 90
+C
+ DO IP = 1, IPTOT
+ CPOL(IA,IP) = RINP(IPOL(IP))
+ ENDDO
+C
+ DO JP = 1, JPTOT
+ DO N = 1, NBL
+ IS1 = 2*N-1
+ IS2 = 2*N
+ CPOLSD(IA,IS1,JP) = RINP(ISPOL(1,JP)+2*(N-1))
+ CPOLSD(IA,IS2,JP) = RINP(ISPOL(2,JP)+2*(N-1))
+ ENDDO
+ ENDDO
+C
+ ACL = MAX( CPOL(IA,ICL) , 0.001 )
+C
+C
+C----- try to find Re, Ma, Ncrit, Xtrip in polar data
+ LIRE = .FALSE.
+ LIMA = .FALSE.
+ LINC = .FALSE.
+ LJTP = .FALSE.
+ DO KP = 1, NIPOL
+ IF(IPOL(KP) .EQ. IRE) LIRE = .TRUE.
+ IF(IPOL(KP) .EQ. IMA) LIMA = .TRUE.
+ IF(IPOL(KP) .EQ. INC) LINC = .TRUE.
+ IF(ISPOL(1,KP) .EQ. JTP) LJTP = .TRUE.
+ ENDDO
+C
+ IF(.NOT. LIRE) THEN
+C------ Re was not in polar data... set using header info
+ IF (IRETYP.EQ.1) THEN
+ CPOL(IA,IRE) = REYN1
+ ELSEIF(IRETYP.EQ.2) THEN
+ CPOL(IA,IRE) = REYN1/SQRT(ACL)
+ ELSEIF(IRETYP.EQ.3) THEN
+ CPOL(IA,IRE) = REYN1/ACL
+ ENDIF
+ ENDIF
+C
+ IF(.NOT. LIMA) THEN
+C------ Mach was not in polar data... set using header info
+ IF (IMATYP.EQ.1) THEN
+ CPOL(IA,IMA) = MACH1
+ ELSEIF(IMATYP.EQ.2) THEN
+ CPOL(IA,IMA) = MACH1/SQRT(ACL)
+ ELSEIF(IMATYP.EQ.3) THEN
+ CPOL(IA,IMA) = MACH1/ACL
+ ENDIF
+ ENDIF
+C
+ IF(.NOT. LINC) THEN
+C------ Ncrit was not in polar data... set using header info
+ CPOL(IA,INC) = ACRIT
+ ENDIF
+C
+ IF(.NOT. LJTP) THEN
+C------ set trip data using header info
+ DO IS = 1, 2*NBL
+ CPOLSD(IA,IS,JTP) = XTRIP(IS)
+ ENDDO
+ ENDIF
+C
+ NA = IA
+ ENDIF
+C
+ 60 CONTINUE
+C---- go read next line
+ GO TO 500
+C=============================================================
+C
+ 80 CONTINUE
+C---- if file was opened here, then close it
+ IF(LOPEN) CLOSE(LU)
+ RETURN
+C
+ 90 CONTINUE
+ IF(LOPEN) CLOSE(LU)
+ ERROR = .TRUE.
+ RETURN
+C
+C..........................................
+ 1000 FORMAT(A)
+ END ! POLREAD
+
+
+ SUBROUTINE POLWRIT(LU,FNPOL,ERROR, LHEAD,
+ & NAX, IA1,IA2, CPOL, IPOL,NIPOL,
+ & REYN1,MACH1,ACRIT,XTRIP,
+ & PTRAT,ETAP,
+ & NAME, IRETYP,IMATYP,
+ & ISX,NBL,CPOLSD, JPOL,NJPOL,
+ & CODE,VERSION, LQUERY )
+ INCLUDE 'PINDEX.INC'
+ CHARACTER*(*) FNPOL, NAME
+ LOGICAL ERROR, LHEAD,LQUERY
+ CHARACTER*(*) CODE
+ REAL CPOL(NAX,IPTOT), CPOLSD(NAX,ISX,JPTOT)
+ REAL MACH1, XTRIP(ISX)
+ INTEGER IPOL(IPTOT), JPOL(JPTOT)
+C--------------------------------------------------------
+C Writes polar save file
+C
+C Input:
+C LU logical unit to use for writing
+C FNPOL name of polar file to be read,
+C if FNPOL(1:1).eq.' ', unit LU is assumed
+C to be already open
+C NAX polar point array dimension
+C ISX airfoil side array dimension
+C IA1,IA2 only polar points IA1..IA2 are written
+C CPOL polar coefficients and parameters
+C IPOL(.) indices of data quantities to be written
+C NIPOL number of data quantities to be written
+C REYN1 Reynolds number for CL=1
+C MACH1 Mach number for CL=1
+C ACRIT Critical amplification ratio
+C XTRIP Trip locations
+C PTRAT Actuator disk total-pressure ratio
+C ETAP Actuator disk thermal efficiency
+C NAME airfoil name string
+C IRETYP flag giving type of Re variation with CL
+C IMATYP flag giving type of Ma variation with CL
+C NBL number of airfoil elements
+C CPOLSD airfoil side-related parameters
+C JPOL(.) indices of side data quantities to be written
+C NJPOL number of side data quantities to be written
+C LHEAD T if header and column label are to be written
+C CODE code used to compute polar
+C VERSION code version
+C LQUERY if T, asks permission to overwrite existing file
+C
+C Output:
+C ERROR T if a OPER or WRITE error occurred
+C--------------------------------------------------------
+ CHARACTER*29 LINE1, LINE2
+ CHARACTER*128 LINEL, LINED, LINEF
+ CHARACTER*1 ANS
+ LOGICAL LOPEN
+C
+ ERROR = .FALSE.
+C
+C---- do we have to open the file?
+ LOPEN = FNPOL .NE. ' '
+C
+ IF(LOPEN) THEN
+ OPEN(LU,FILE=FNPOL,STATUS='OLD',ERR=20)
+C
+ IF(LQUERY) THEN
+ WRITE(*,*)
+ WRITE(*,*) 'Output file exists. Overwrite? Y'
+ READ(*,1000) ANS
+C
+ IF(INDEX('Nn',ANS).EQ.0) GO TO 22
+C
+ CLOSE(LU)
+ WRITE(*,*) 'Polar file not saved'
+ RETURN
+ ENDIF
+C
+ 20 OPEN(LU,FILE=FNPOL,STATUS='UNKNOWN',ERR=90)
+ 22 REWIND(LU)
+ ENDIF
+C
+ IF(LHEAD) THEN
+ WRITE(LU,*) ' '
+ WRITE(LU,8000) CODE, VERSION
+ WRITE(LU,*) ' '
+ IF(NBL.EQ.1) THEN
+ WRITE(LU,9001) NAME
+ ELSE
+ WRITE(LU,9002) NAME, NBL
+ ENDIF
+C
+ IFFBC = 0
+ ISMOM = 0
+C
+ IF(IFFBC.NE.0 .AND. ISMOM.NE.0) THEN
+ IF(IFFBC.EQ.1) LINE1 = ' Solid wall far field '
+ IF(IFFBC.EQ.2) LINE1 = ' Vortex + doublet far field '
+ IF(IFFBC.EQ.3) LINE1 = ' Constant pressure far field '
+ IF(IFFBC.EQ.4) LINE1 = ' Supersonic wave far field '
+ IF(IFFBC.GE.5) LINE1 = ' '
+ IF(ISMOM.EQ.1) LINE2 = ' S-momentum conserved '
+ IF(ISMOM.EQ.2) LINE2 = ' Entropy conserved '
+ IF(ISMOM.EQ.3) LINE2 = ' Entropy conserved near LE '
+ IF(ISMOM.EQ.4) LINE2 = ' S-mom conserved at shocks '
+ IF(ISMOM.GE.5) LINE2 = ' '
+ WRITE(LU,9006) LINE1, LINE2
+ 9006 FORMAT(1X,3X,2A29)
+ ENDIF
+C
+ WRITE(LU,*) ' '
+C
+ LINE1 = ' '
+ LINE2 = ' '
+ IF(IRETYP.EQ.1) LINE1 = ' Reynolds number fixed '
+ IF(IRETYP.EQ.2) LINE1 = ' Reynolds number ~ 1/sqrt(CL)'
+ IF(IRETYP.EQ.3) LINE1 = ' Reynolds number ~ 1/CL '
+ IF(IMATYP.EQ.1) LINE2 = ' Mach number fixed '
+ IF(IMATYP.EQ.2) LINE2 = ' Mach number ~ 1/sqrt(CL) '
+ IF(IMATYP.EQ.3) LINE2 = ' Mach number ~ 1/CL '
+ WRITE(LU,9005) IRETYP, IMATYP, LINE1, LINE2
+C
+ WRITE(LU,*) ' '
+ DO N = 1, NBL
+ IS1 = 2*N-1
+ IS2 = 2*N
+ IF(NBL.EQ.1) THEN
+ WRITE(LU,9011) XTRIP(IS1), XTRIP(IS2)
+ ELSE
+ WRITE(LU,9012) XTRIP(IS1), XTRIP(IS2), N
+ ENDIF
+ ENDDO
+ WRITE(LU,9015) MACH1, REYN1/1.0E6, ACRIT
+ IF(PTRAT .NE. 0.0) WRITE(LU,9017) PTRAT, ETAP
+ WRITE(LU,*) ' '
+C
+ LINEL = ' '
+ LINED = ' '
+C
+ KL = 1
+ KD = 1
+C
+ DO 30 KP = 1, NIPOL
+ IP = IPOL(KP)
+ IF(IP.EQ.0) GO TO 30
+C
+ KDOT = INDEX(CPOLFORM(IP),'.')
+ IF(KDOT.EQ.0) KDOT = LEN(CPOLFORM(IP))
+ READ(CPOLFORM(IP)(2:KDOT-1),*,ERR=95) NFORM
+C
+ CALL STRIP(CPOLNAME(IP),NNAME)
+ NBLANK = MAX( (NFORM-NNAME+2)/2 , 0 )
+C
+ LINEL(KL+1+NBLANK:KL+NNAME+NBLANK) = CPOLNAME(IP)(1:NNAME)
+ KL = KL + NFORM
+C
+ LINED(KD+2:KD+NFORM) = '--------------------------------'
+ KD = KD + NFORM
+ 30 CONTINUE
+C
+ DO 32 KP = 1, NJPOL
+ JP = JPOL(KP)
+ IF(JP.EQ.0) GO TO 32
+C
+ KDOT = INDEX(CPOLSFORM(JP),'.')
+ IF(KDOT.EQ.0) KDOT = LEN(CPOLSFORM(JP))
+ READ(CPOLSFORM(JP)(2:KDOT-1),*,ERR=95) NFORM
+C
+ CALL STRIP(CPOLSNAME(JP),NNAME)
+ NBLANK = MAX( (NFORM-NNAME-2)/2 , 0 )
+C
+ DO N = 1, NBL
+ LINEL(KL+1+NBLANK:KL+4+NNAME+NBLANK) =
+ & 'Top_' // CPOLSNAME(JP)(1:NNAME)
+ KL = KL + NFORM
+C
+ LINED(KD+2:KD+NFORM) = '--------------------------------'
+ KD = KD + NFORM
+C
+ LINEL(KL+1+NBLANK:KL+4+NNAME+NBLANK) =
+ & 'Bot_' // CPOLSNAME(JP)(1:NNAME)
+ KL = KL + NFORM
+C
+ LINED(KD+2:KD+NFORM) = '--------------------------------'
+ KD = KD + NFORM
+ ENDDO
+ 32 CONTINUE
+C
+C
+C
+C LINEL =
+C & ' alpha CL CD CDp CM Top_Xtr Bot_Xtr'
+CCC 1234567890123456789012345678901234567890123456789012345678901234567890
+C K = 62
+C
+C
+C LINEL =
+C & ' ------- -------- --------- --------- -------- ------- -------'
+CCC 3.453 1.3750 0.00921 0.00512 -0.1450 0.9231 0.5382
+CCC 3.453 1.3750 0.00921 0.00213 -0.1450 0.9231 0.5382
+C K = 62
+
+ WRITE(LU,1000) LINEL(1:KL)
+ WRITE(LU,1000) LINED(1:KD)
+C
+ ENDIF
+C
+
+ LINEF = '(1X'
+ KF = 3
+ DO KP = 1, NIPOL
+ IP = IPOL(KP)
+ NF = LEN(CPOLFORM(IP))
+C
+ LINEF(KF+1:KF+NF+1) = ',' // CPOLFORM(IP)
+ KF = KF + NF + 1
+ ENDDO
+ DO KP = 1, NJPOL
+ JP = JPOL(KP)
+ NF = LEN(CPOLSFORM(JP))
+C
+ DO N = 1, NBL
+ LINEF(KF+1:KF+NF+1) = ',' // CPOLSFORM(JP)
+ KF = KF + NF + 1
+C
+ LINEF(KF+1:KF+NF+1) = ',' // CPOLSFORM(JP)
+ KF = KF + NF + 1
+ ENDDO
+ ENDDO
+ LINEF(KF+1:KF+1) = ')'
+ KF = KF + 1
+C
+C
+ DO 40 IA = IA1, IA2
+ WRITE(LU,LINEF)
+ & (CPOL(IA,IPOL(KP)), KP=1, NIPOL),
+ & ((CPOLSD(IA,IS,JPOL(KP)), IS=1, 2*NBL), KP=1, NJPOL)
+ 40 CONTINUE
+C
+C
+ 80 CONTINUE
+C---- if file was opened here, then close it
+ IF(LOPEN) CLOSE(LU)
+ RETURN
+C
+ 90 CONTINUE
+ ERROR = .TRUE.
+ RETURN
+C
+ 95 CONTINUE
+ WRITE(*,*) '? Bad format specification in PINDEX.INC'
+ STOP
+C
+C......................................................................
+ 1000 FORMAT(A)
+ 8000 FORMAT(7X,A,9X,'Version', F5.2)
+ 9001 FORMAT(1X,'Calculated polar for: ', A)
+ 9002 FORMAT(1X,'Calculated polar for: ', A, I4,' elements')
+ 9005 FORMAT(1X,I1,I2,2A29)
+ 9011 FORMAT(1X,
+ &'xtrf = ',F7.3,' (top) ',F9.3,' (bottom) ')
+ 9012 FORMAT(1X,
+ &'xtrf = ',F7.3,' (top) ',F9.3,' (bottom) element', I3)
+ 9015 FORMAT(1X,
+ &'Mach = ',F7.3,5X,'Re = ',F9.3,' e 6',5X,'Ncrit = ',F7.3)
+ 9017 FORMAT(1X,
+ &'pi_p = ',F7.4,5X,'eta_p = ',F9.4)
+ 9100 FORMAT(1X,F7.3,F9.4,2F10.5,F9.4,2F8.4 , F9.5)
+CCC 3.453 1.3750 0.00921 0.500 -0.1450 0.9231 0.5382 -0.00942
+CCC 3.453 1.3750 0.00921 0.500 -0.1450 0.9231 0.5382
+ END
+
+
+
+ SUBROUTINE POLREF(LU,FNREF,ERROR,
+ & NFX,NF,XYREF,LABREF )
+ INCLUDE 'PINDEX.INC'
+ CHARACTER*(*) FNREF,LABREF
+ LOGICAL ERROR
+ DIMENSION NF(4)
+ DIMENSION XYREF(NFX,2,4)
+C--------------------------------------------------------
+C Reads in polar reference data file
+C
+C Input:
+C LU logical unit to use for reading
+C FNREF name of polar file to be read,
+C if FNREF(1:1).eq.' ', unit LU is assumed
+C to be already open
+C NFX polar point array dimension
+C
+C Output:
+C ERROR T if a READ error occurred
+C NF(.) number of points in each data block
+C XYREF(...) reference polar data
+C LABREF(.) reference polar label
+C--------------------------------------------------------
+ LOGICAL LOPEN
+ CHARACTER*80 LINE
+C
+ ERROR = .FALSE.
+ LOPEN = FNREF(1:1) .NE. ' '
+ IF(LOPEN) OPEN(LU,FILE=FNREF,STATUS='OLD',ERR=900)
+C
+C---- try to read data label
+ READ(LU,1000,END=900) LINE
+ 1000 FORMAT(A)
+C
+C---- set data label if present
+ IF(LINE(1:1).EQ.'#') THEN
+ LABREF = LINE(2:80)
+ ELSE
+ LABREF = ' '
+ REWIND(LU)
+ ENDIF
+C
+ DO 100 K=1, 4
+ DO 10 I=1, NFX
+ READ(LU,*,END=11,ERR=900) XYREF(I,1,K), XYREF(I,2,K)
+ IF(XYREF(I,1,K) .EQ. 999.0) GO TO 11
+ 10 CONTINUE
+ 11 NF(K) = I-1
+ 100 CONTINUE
+ IF(LOPEN) CLOSE(LU)
+ RETURN
+C
+ 900 CONTINUE
+ ERROR = .TRUE.
+C
+ RETURN
+ END ! POLREF