aboutsummaryrefslogtreecommitdiff
path: root/orrs/src/osrun.f
diff options
context:
space:
mode:
authorDimitri Sokolyuk <demon@dim13.org>2009-05-11 00:27:49 +0000
committerDimitri Sokolyuk <demon@dim13.org>2009-05-11 00:27:49 +0000
commit0d4f43d355de79178b1142e9735902cf641670b6 (patch)
tree2ced2323f6351db2a51090b3fd13eb11f69ff53f /orrs/src/osrun.f
Xfoil 6.97
Diffstat (limited to 'orrs/src/osrun.f')
-rwxr-xr-xorrs/src/osrun.f528
1 files changed, 528 insertions, 0 deletions
diff --git a/orrs/src/osrun.f b/orrs/src/osrun.f
new file mode 100755
index 0000000..8b8ee33
--- /dev/null
+++ b/orrs/src/osrun.f
@@ -0,0 +1,528 @@
+
+ PROGRAM OSRUN
+C---------------------------------------------------------------------------
+C Program for executing and displaying Orr-Sommerfeld solution
+C Usage:
+C
+C % osrun [ inputfile ]
+C
+C The inputfile contains mean-flow profile data:
+C
+C y1 U1 S1
+C y2 U2 S1
+C . . .
+C . . .
+C
+C where S = dU/dy .
+C If the optional argument "inputfile" is missing, then
+C a Falkner-Skan profile will be generated, for a specified
+C H or betaU = x/ue due/dx
+C
+C In either case, the user is prompted for R_theta and w_r values.
+C
+C---------------------------------------------------------------------------
+C
+ PARAMETER (NMAX=2001)
+ DIMENSION ETA(NMAX), F(NMAX), U(NMAX), S(NMAX)
+ DIMENSION UTR(NMAX), UTI(NMAX), UT(NMAX),
+ & VTR(NMAX), VTI(NMAX), VT(NMAX),
+ & WTR(NMAX), WTI(NMAX), WT(NMAX),
+ & CTR(NMAX), CTI(NMAX), CT(NMAX),
+ & PTR(NMAX), PTI(NMAX), PT(NMAX)
+ DIMENSION UU(NMAX), VV(NMAX), UV(NMAX), QQ(NMAX)
+ DIMENSION TVT(NMAX), TV(NMAX)
+ CHARACTER*1 ANS
+ CHARACTER*80 FNAME, ARGP1
+ DIMENSION XLIN(2), YLIN(2)
+C
+ DIMENSION AINPUT(10)
+ LOGICAL ERROR
+C
+ DATA CV1 / 7.1 /
+C
+ LST = 1
+ LRE = 1
+C
+ IDEV = 1
+ IDEVRP = 2
+ IPSLU = 0
+C
+ SIZE = 6.0
+ PAR = 0.75
+C
+ CALL PLINITIALIZE
+C
+C
+ N = 1001
+ ETAE = 16.0
+ GEO = 1.01
+C
+ CH = 0.021
+C
+ IF(N.GT.NMAX) STOP 'TEST: Array overflow.'
+C
+ CALL GETARG0(1,ARGP1)
+ IF(ARGP1(1:1).EQ.' ') GO TO 50
+C
+ FNAME = ARGP1
+C
+C---- try formatted read first
+ OPEN(1,FILE=FNAME,STATUS='OLD',ERR=50)
+ DO I=1, NMAX
+ READ(1,*,ERR=30,END=25) ETA(I), U(I), S(I)
+ ENDDO
+C
+ 25 CLOSE(1)
+ N = I-1
+ GO TO 80
+C
+C---- now try unformatted read
+ 30 CONTINUE
+ OPEN(19,FILE=FNAME,STATUS='OLD',FORM='UNFORMATTED',ERR=50)
+ DO I=1, NMAX
+ READ(19,ERR=50,END=35) ETA(I), U(I), S(I)
+ ENDDO
+C
+ 35 CLOSE(19)
+ N = I-1
+ GO TO 80
+C
+C------------------------------------------------------------------------
+C---- no argument specified or read error... get Falkner-Skan parameter
+ 50 CONTINUE
+ WRITE(*,*) 'Enter Falkner-Skan parameter Beta (or H)'
+ READ (*,*) PARM
+C
+ IF(PARM .GT. 1.0) THEN
+ write(*,*) 'Enter max y/theta, GEO'
+ read (*,*) etae, geo
+ H = PARM
+ CALL FS(3,2,BU,H,N,ETAE,GEO,ETA,F,U,S,DELTA)
+C
+ ELSE
+ write(*,*) 'Enter max y/theta, GEO'
+ read (*,*) etae, geo
+ BU = PARM / (2.0 - PARM)
+ CALL FS(3,1,BU,H,N,ETAE,GEO,ETA,F,U,S,DELTA)
+C
+ ENDIF
+C
+ WRITE(*,*) 'BetaU, H =', BU, H
+ccc GO TO 90
+C
+C------------------------------------------------------------------------
+C---- normalize input profiles
+ 80 CONTINUE
+ DS = 0.
+ TH = 0.
+ DO I = 1, N-1
+ DY = ETA(I+1) - ETA(I)
+ UA = (U(I+1) + U(I))*0.5 / U(N)
+ DS = DS + (1.0 - UA) *DY
+ TH = TH + (1.0 - UA)*UA*DY
+ ENDDO
+C
+ H = DS/TH
+ DO I = 1, N
+ ETA(I) = ETA(I) / TH
+ S(I) = S(I) * TH
+ ENDDO
+ ETAE = ETA(N)
+C
+ WRITE(*,*) 'H =', H
+C
+C---------------------
+C
+ 90 CONTINUE
+
+
+ do i = 1, n
+ write(*,'(1x,i5,3g14.6)') i, eta(i), u(i), s(i)
+ enddo
+
+
+ EWT = 1.0/ETAE
+C
+ DETA = 1.0
+ IF(ETAE .GT. 16.01) DETA = 2.0
+ IF(ETAE .GT. 30.01) DETA = 5.0
+ ETAE = DETA * AINT( ETAE/DETA + 0.51 )
+C
+ CALL PLOPEN(0.7,IPSLU,IDEV)
+ CALL PLOTABS(0.5,1.0,-3)
+C
+ CALL NEWFACTOR(SIZE)
+ CALL GETCOLOR(ICOL0)
+C
+ CALL PLOT(0.5*PAR,0.0,-3)
+C
+ CALL NEWPEN(1)
+ CALL XAXIS(0.0,0.0,-PAR, 0.2*PAR,0.0, 0.2,0.7*CH,1)
+ CALL YAXIS(0.0,0.0,-1.0,DETA*EWT,0.0,DETA,0.7*CH,1)
+C
+ CALL NEWPEN(2)
+ XL = -4.0*CH
+ YL = (ETAE-1.5*DETA)*EWT - 0.5*CH
+ CALL PLCHAR(XL,YL,CH,'y/ ',0.0,3)
+ CALL PLMATH(XL,YL,CH,' q',0.0,3)
+C
+C
+ UWT = PAR
+ CALL NEWPEN(4)
+ CALL XYLINE(N,U,ETA,0.0,UWT,0.0,EWT,1)
+ CALL NEWPEN(3)
+ CALL XYSYMB(N,U,ETA,0.0,UWT,0.0,EWT,0.125*CH,1)
+C
+ CALL PLFLUSH
+C
+ RE = 100.0
+ OMEGAR = 0.1
+ ALPHAR = 2.0*OMEGAR
+ ALPHAI = 0.
+C
+ DO 100 IPASS=1, 50
+ ITMAX = 20
+C
+ 95 AINPUT(1) = RE
+ AINPUT(2) = OMEGAR
+ WRITE(*,2100) AINPUT(1), AINPUT(2)
+ 2100 FORMAT(1X,' Enter Rtheta, Wreal:', F9.1, F10.5)
+ CALL READR(2,AINPUT,ERROR)
+ IF(ERROR) GO TO 95
+C
+ RE = AINPUT(1)
+ OMEGAR = AINPUT(2)
+C
+ IF(RE .EQ. 0.0) GO TO 101
+C
+c RD = RE*H
+c WR = OMEGAR/RE
+c WRITE(*,*) ' '
+c WRITE(*,*) 'Rd* =', RD, ' Wr/Rth =', WR
+C
+ OMEGAI = 0.0
+C
+ 97 AINPUT(1) = ALPHAR
+ AINPUT(2) = ALPHAI
+ WRITE(*,2200) AINPUT(1), AINPUT(2)
+ 2200 FORMAT(1X,' Enter initial ar, ai:', 2F10.5)
+ CALL READR(2,AINPUT,ERROR)
+ IF(ERROR) GO TO 97
+C
+ ALPHAR = AINPUT(1)
+ ALPHAI = AINPUT(2)
+C
+ ITLIM = ITMAX
+ CALL ORRS(LST,LRE,N,ETA,U,S, RE, ITLIM,
+ & ALPHAR,ALPHAI, OMEGAR,OMEGAI,
+ & UTR,UTI,VTR,VTI,WTR,WTI,CTR,CTI, DELMAX)
+C
+ CALL OSPRES(N,ETA,U, ALPHAR,ALPHAI, VTR,VTI, PTR,PTI )
+C
+ DO I=1, N
+ UT(I) = SQRT(UTR(I)**2 + UTI(I)**2)
+ VT(I) = SQRT(VTR(I)**2 + VTI(I)**2)
+ PT(I) = SQRT(PTR(I)**2 + PTI(I)**2)
+ UU(I) = 0.5*(UTR(I)*UTR(I) + UTI(I)*UTI(I))
+ VV(I) = 0.5*(VTR(I)*VTR(I) + VTI(I)*VTI(I))
+ UV(I) = 0.5*(UTR(I)*VTR(I) + UTI(I)*VTI(I))
+ QQ(I) = UU(I) + VV(I)
+ SLIM = MAX( S(I) , 1.0E-5 )
+ TVT(I) = ABS( -UV(I) / SLIM )
+C
+ TV(I) = SQRT( SQRT( TVT(I)*CV1**3/RE**3 ) )
+ IF(I .GE. 3) TV(I) = TV(I-1)
+C
+ IF(TV(I) .GT. 0.0) THEN
+ DO ITER = 1, 10
+ CHI = RE*TV(I)
+ CHI3 = CHI**3
+ CHI3_N = 3.0*CHI**2 * RE
+ RES = TVT(I)*(CHI3 + CV1**3) - TV(I)*CHI3
+ RES_N = TVT(I)*CHI3_N - TV(I)*CHI3_N - CHI3
+ DN = -RES/RES_N
+c write(*,'(1x,2g13.6,2e12.4)') tvt(i), tv(i), res, dn
+ TV(I) = TV(I) + DN
+c pause
+ ENDDO
+ ENDIF
+C
+ ENDDO
+C
+ QTHIK = 0.
+ DQTDX = 0.
+ UPRES = 0.
+ TWORK = 0.
+ DISS1 = 0.
+ DISS2 = 0.
+ DISS3 = 0.
+ PQINT = 0.
+ DO I = 2, N
+ UA = (U(I) + U(I-1))*0.5
+ DU = U(I) - U(I-1)
+ DY = ETA(I) - ETA(I-1)
+C
+ URA = (UTR(I) + UTR(I-1))*0.5
+ UIA = (UTI(I) + UTI(I-1))*0.5
+ VRA = (VTR(I) + VTR(I-1))*0.5
+ VIA = (VTI(I) + VTI(I-1))*0.5
+ WRA = (WTR(I) + WTR(I-1))*0.5
+ WIA = (WTI(I) + WTI(I-1))*0.5
+ PRA = (PTR(I) + PTR(I-1))*0.5
+ PIA = (PTI(I) + PTI(I-1))*0.5
+C
+ DUR = UTR(I) - UTR(I-1)
+ DUI = UTI(I) - UTI(I-1)
+ DVR = VTR(I) - VTR(I-1)
+ DVI = VTI(I) - VTI(I-1)
+ DWR = WTR(I) - WTR(I-1)
+ DWI = WTI(I) - WTI(I-1)
+C
+ QTHIK = QTHIK + 0.25*(UU(I)+UU(I-1)
+ & +VV(I)+VV(I-1))*UA*DY
+C
+ UDUDX = - (ALPHAI*URA + ALPHAR*UIA)*URA
+ & + (ALPHAR*URA - ALPHAI*UIA)*UIA
+ VDVDX = - (ALPHAI*VRA + ALPHAR*VIA)*VRA
+ & + (ALPHAR*VRA - ALPHAI*VIA)*VIA
+C
+ PDUDX = - (ALPHAI*URA + ALPHAR*UIA)*PRA
+ & + (ALPHAR*URA - ALPHAI*UIA)*PIA
+ UDPDX = - (ALPHAI*PRA + ALPHAR*PIA)*URA
+ & + (ALPHAR*PRA - ALPHAI*PIA)*UIA
+C
+ DQTDX = DQTDX + 0.5*(UDUDX + VDVDX)*DY * UA
+C
+ UPRES = UPRES - 0.5*(UDPDX + PDUDX)*DY
+C
+ TWORK = TWORK - 0.50*(UV(I)+UV(I-1))*DU
+C
+ DISS1 = DISS1 + ( ALPHAI*URA + ALPHAR*UIA )**2 * DY
+ & + ( ALPHAR*URA - ALPHAI*UIA )**2 * DY
+C
+ DISS2 = DISS2 + ( DVR**2 + DVI**2 ) / DY
+C
+ DISS3 = DISS3
+ & + 0.5 * ( DUR/DY - ALPHAI*VRA - ALPHAR*VIA )**2 * DY
+ & + 0.5 * ( DUI/DY + ALPHAR*VRA - ALPHAI*VIA )**2 * DY
+C
+ PQINT = PQINT
+ & - 0.5*URA*DWR - 0.5*(ALPHAI*WRA + ALPHAR*WIA)*VRA * DY
+ & - 0.5*UIA*DWI + 0.5*(ALPHAR*WRA - ALPHAI*WIA)*VIA * DY
+ ENDDO
+C
+ DISS1 = DISS1 / RE
+ DISS2 = DISS2 / RE
+ DISS3 = DISS3 / RE
+ PQINT = PQINT / RE
+C
+C
+ DQTDX = DQTDX / QTHIK
+ UPRES = UPRES / QTHIK
+ TWORK = TWORK / QTHIK
+ DISS1 = DISS1 / QTHIK
+ DISS2 = DISS2 / QTHIK
+ DISS3 = DISS3 / QTHIK
+ PQINT = PQINT / QTHIK
+C
+ DISS = DISS1 + DISS2 + DISS3
+C
+ WRITE(*,*)
+ WRITE(*,*) 'dEdx, P+Dx+D :',DQTDX,TWORK+UPRES+PQINT
+ WRITE(*,*) 'P Dx D e:', TWORK, UPRES, PQINT, -DISS
+ WRITE(*,*)
+C
+ IF(IPASS.EQ.1) THEN
+ CALL SCALIT(N,UT,0.0,USF,ANN,NANN)
+ UWT = PAR*USF
+C
+ CALL SCALIT(N,VT,0.0,VSF,ANN,NANN)
+ VWT = PAR*VSF
+C
+ CALL SCALIT(N,PT,0.0,PSF,ANN,NANN)
+ PWT = PAR*PSF
+C
+ CALL SCALIT(N,QQ,0.0,TSF,ANN,NANN)
+ TWT = PAR*TSF
+C
+ EOFF = 0.
+ UOFF = 0.
+ POFF = 0.
+ TOFF = 0.
+C
+ PWT = UWT
+ ENDIF
+C
+C
+ CALL NEWPEN(3)
+C
+ XL = PAR + 5.0*CH
+C
+ YL = ETAE*EWT
+ CALL PLCHAR(XL ,YL , CH,'H = ',0.0, 8)
+ CALL PLNUMB(XL+8.0*CH,YL , CH, H ,0.0, 3)
+C
+ YL = YL - 3.5*CH
+ CALL PLCHAR(XL ,YL , CH,'Re = ',0.0, 8)
+ CALL PLMATH(XL+1.9*CH,YL-0.4*CH,0.8*CH, 'q' ,0.0, 1)
+ CALL PLNUMB(XL+8.0*CH,YL , CH, RE ,0.0,-1)
+C
+ YL = YL - 2.5*CH
+ CALL PLMATH(XL ,YL , CH,'w q/ = ',0.0, 8)
+ CALL PLCHAR(XL ,YL , CH,' U ',0.0, 8)
+ CALL PLCHAR(XL+0.9*CH,YL-0.3*CH,0.8*CH, 'r' ,0.0, 1)
+ CALL PLNUMB(XL+8.0*CH,YL , CH, OMEGAR ,0.0, 5)
+C
+ YL = YL - 3.5*CH
+ CALL PLMATH(XL ,YL , CH,'a q = ',0.0, 8)
+ CALL PLCHAR(XL+0.9*CH,YL-0.3*CH,0.8*CH, 'r' ,0.0, 1)
+ CALL PLNUMB(XL+8.0*CH,YL , CH, ALPHAR ,0.0, 5)
+C
+ YL = YL - 2.5*CH
+ CALL PLMATH(XL ,YL , CH,'a q = ',0.0, 8)
+ CALL PLCHAR(XL+0.9*CH,YL-0.3*CH,0.8*CH, 'i' ,0.0, 1)
+ CALL PLNUMB(XL+8.0*CH,YL , CH, ALPHAI ,0.0, 5)
+C
+C
+C
+ XLIN(1) = -7.5*CH
+ XLIN(2) = -1.5*CH
+ YLIN(1) = 0.5*CH
+ YLIN(2) = 0.5*CH
+C
+ CALL NEWPEN(2)
+C
+ XL = PAR + 12.0*CH
+ YL = 0.50*ETAE*EWT
+C
+ CALL NEWCOLORNAME('red')
+ CALL XYLINE(N,UTR,ETA,UOFF,UWT,EOFF,EWT,2)
+ CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,2)
+ CALL PLCHAR(XL ,YL , CH,'u /U',0.0,4)
+ CALL PLMATH(XL ,YL , CH,' ` ',0.0,4)
+ CALL PLCHAR(XL+0.9*CH,YL-0.3*CH,0.8*CH, 'r' ,0.0,1)
+C
+ YL = YL - 2.5*CH
+C
+ CALL NEWCOLORNAME('orange')
+ CALL XYLINE(N,UTI,ETA,UOFF,UWT,EOFF,EWT,3)
+ CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,3)
+ CALL PLCHAR(XL ,YL , CH,'u /U',0.0,4)
+ CALL PLMATH(XL ,YL , CH,' ` ',0.0,4)
+ CALL PLCHAR(XL+0.9*CH,YL-0.3*CH,0.8*CH, 'i' ,0.0,1)
+C
+ YL = YL - 2.5*CH
+C
+ CALL NEWCOLORNAME('yellow')
+ CALL XYLINE(N,UT ,ETA,UOFF,UWT,EOFF,EWT,1)
+ CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,1)
+ CALL PLCHAR(XL ,YL , CH,'u /U',0.0,4)
+ CALL PLMATH(XL ,YL , CH,' ` ',0.0,4)
+ CALL PLCHAR(XL-0.6*CH,YL ,0.9*CH,'| |' ,0.0,3)
+C
+C
+ YL = YL - 3.5*CH
+C
+ CALL NEWCOLORNAME('violet')
+ CALL XYLINE(N,PTR,ETA,POFF,PWT,EOFF,EWT,5)
+ CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,5)
+ CALL PLCHAR(XL ,YL , CH,'p / U ',0.0,6)
+ CALL PLMATH(XL ,YL , CH,' ` r 2',0.0,6)
+ CALL PLCHAR(XL+0.9*CH,YL-0.3*CH,0.8*CH, 'r' ,0.0,1)
+C
+ YL = YL - 2.5*CH
+C
+ CALL NEWCOLORNAME('blue')
+ CALL XYLINE(N,PTI,ETA,POFF,PWT,EOFF,EWT,6)
+ CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,6)
+ CALL PLCHAR(XL ,YL , CH,'p / U ',0.0,6)
+ CALL PLMATH(XL ,YL , CH,' ` r 2',0.0,6)
+ CALL PLCHAR(XL+0.9*CH,YL-0.3*CH,0.8*CH, 'i' ,0.0,1)
+C
+ YL = YL - 2.5*CH
+C
+ CALL NEWCOLORNAME('cyan')
+ CALL XYLINE(N,PT ,ETA,POFF,PWT,EOFF,EWT,1)
+ CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,1)
+ CALL PLCHAR(XL ,YL , CH,'p / U ',0.0,6)
+ CALL PLMATH(XL ,YL , CH,' ` r 2',0.0,6)
+ CALL PLCHAR(XL-0.6*CH,YL ,0.9*CH,'| |' ,0.0,3)
+C
+ YL = YL - 3.5*CH
+C
+ CALL NEWCOLORNAME('green')
+ CALL XYLINE(N,QQ ,ETA,TOFF,TWT,EOFF,EWT,2)
+ CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,2)
+ CALL PLMATH(XL ,YL+0.2*CH, CH,' ___ ' ,0.0,5)
+ CALL PLMATH(XL ,YL , CH,' ` ` 2',0.0,8)
+ CALL PLCHAR(XL ,YL , CH,' q q /U ',0.0,8)
+C
+ YL = YL - 2.5*CH
+C
+ CALL NEWCOLORNAME('green')
+ CALL XYLINE(N,UV ,ETA,TOFF,-10.0*TWT,EOFF,EWT,1)
+ CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,1)
+ CALL PLMATH(XL ,YL+0.2*CH, CH,' ___ ',0.0,12)
+ CALL PLMATH(XL ,YL , CH,' ` ` 2 # ',0.0,12)
+ CALL PLCHAR(XL ,YL , CH,'-u v /U 10',0.0,12)
+C
+ YL = YL - 3.5*CH
+C
+ CALL NEWCOLORNAME('magenta')
+ CALL XYLINE(N,TVT,ETA,TOFF,TWT,EOFF,EWT,1)
+ CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,1)
+ CALL PLMATH(XL ,YL+0.2*CH, CH,'n 2 ',0.0,9)
+ CALL PLCHAR(XL ,YL , CH,' /U ',0.0,9)
+ CALL PLCHAR(XL+0.9*CH,YL-0.3*CH,0.8*CH, 't' ,0.0,1)
+C
+ YL = YL - 2.5*CH
+C
+ CALL NEWCOLORNAME('magenta')
+ CALL XYLINE(N,TV,ETA,TOFF,TWT,EOFF,EWT,2)
+ CALL XYLINE(2,XLIN,YLIN,-XL,1.0,-YL,1.0,2)
+ CALL PLMATH(XL ,YL+0.2*CH, CH,'n 2',0.0,5)
+ CALL PLMATH(XL ,YL+0.2*CH, CH,'~' ,0.0,1)
+ CALL PLCHAR(XL ,YL , CH,' /U ',0.0,5)
+C
+ CALL NEWCOLOR(ICOL0)
+C
+ CALL PLFLUSH
+C
+ 99 CONTINUE
+ WRITE(*,*) 'Zoom, Unzoom, Annotate, Hardcopy, Dump ?'
+ READ (*,1000) ANS
+ 1000 FORMAT(A)
+C
+ IF (INDEX('Zz',ANS) .NE. 0) THEN
+ CALL USETZOOM(.FALSE.,.TRUE.)
+ CALL REPLOT(IDEV)
+ GO TO 99
+ ELSEIF(INDEX('Uu',ANS) .NE. 0) THEN
+ CALL CLRZOOM
+ CALL REPLOT(IDEV)
+ GO TO 99
+ ELSEIF(INDEX('Aa',ANS) .NE. 0) THEN
+ CALL ANNOT(CH)
+ GO TO 99
+ ELSEIF(INDEX('Hh',ANS) .NE. 0) THEN
+ CALL REPLOT(IDEVRP)
+ GO TO 99
+ ELSEIF(INDEX('Dd',ANS) .NE. 0) THEN
+ LU = 9
+ WRITE(LU,9922) ALPHAR,ALPHAI, OMEGAR,OMEGAI
+ DO I = 1, N
+ 9922 FORMAT(1X, 8E16.7)
+ WRITE(LU,9922)
+ & UTR(I), UTI(I), VTR(I), VTI(I), PTR(I), PTI(I)
+ ENDDO
+ WRITE(*,*) 'Written to fort.9'
+C
+ ENDIF
+C
+ 100 CONTINUE
+ 101 CONTINUE
+C
+ CALL PLCLOSE
+ STOP
+ END
+