program prandtline
      parameter(jxx=201,lxx=101)
      real c(jxx),g(jxx),dg(jxx),y(jxx),eta(jxx)
      real w(jxx),t(jxx),dem(jxx)
      real a0(jxx),a1(jxx),a2(jxx),b0(jxx),b1(jxx)
      real l(jxx),d(jxx),q(jxx),at(jxx)
      real cx(lxx),cz(lxx),cq(lxx),inc(lxx)
      real cmf(jxx),cmt(jxx),fz(jxx)
      integer m(jxx)
      integer kxtrm(lxx)
      character*4 title(18)
      data a0/jxx*0./a1/jxx*0./a2/jxx*0./
      data m/jxx*2/
      open(unit=13,file='polarbl.dat',form='formatted')
      open(unit=14,file='polarbl.cdcl',form='formatted')
      open(unit=15,file='prandtline.data',form='formatted')
      open(unit=16,file='prandtline.in',form='formatted')
      open(unit=17,file='prandtline.ycdt',form='formatted')
      open(unit=18,file='prandtline.ygw',form='formatted')
      open(unit=19,file='prandtline.yg',form='formatted')
      open(unit=20,file='prandtline.yw',form='formatted')
      open(unit=21,file='prandtline.ycl',form='formatted')
      open(unit=22,file='prandtline.yat',form='formatted')
      open(unit=23,file='prandtline.cdcl',form='formatted')
      open(unit=24,file='prandtline.out',form='formatted')
      open(unit=25,file='prandtline.clcdcq',form='formatted')
      open(unit=26,file='prandtline.ycmf',form='formatted')
      open(unit=27,file='prandtline.yfz',form='formatted')
      open(unit=28,file='prandtline.ycmt',form='formatted')
c*****constants
      eps=1.e-7
      pi=2.*asin(1.)
      degrad=pi/180.
c*****polar data
      write(6,*)
      write(6,*)'******do you want to use polar data? Y/N=1/0'
      read(5,*)ipolar
      if(ipolar.ne.1)goto 4
      read(13,1000)title
      write(6,1000)title
      read(13,1000)title
      write(6,1000)title
      read(13,1000)title
      write(6,1000)title
      read(13,1000)title
      write(6,1000)title
      read(13,1000)title
      write(6,1000)title
      read(13,1000)title
      write(6,1000)title
      read(13,1000)title
      write(6,1000)title
      read(13,1000)title
      write(6,1000)title
      read(13,1000)title
      write(6,1000)title
      read(13,1000)title
      write(6,1000)title
      read(13,1000)title
      write(6,1000)title
      read(13,1000)title
      write(6,1000)title
      write(6,*)
      write(6,*)'*****extrema of the Cl(alpha) function:'
      prod=1.
      km=1
      kfirst=0
      do 1 k=1,lxx
         read(13,*,end=2)inc(k),cz(k),cx(k),dum,cq(k)
         cq(k)=cq(k)
         kxtrm(k)=0
         dcz=cz(k)-cz(km)+eps
         prod=prod*dcz
         if(prod.lt.-eps)then
            write(6,*)'kxtrm=',km,' cz(kxtrm)=',cz(km)
            kxtrm(km)=km
            if(kfirst.le.0)then
               mxtrm=km
               kfirst=1
            endif
         endif
         prod=sign(1.,dcz)
         km=k
 1    continue
 2    continue
      kx=k-1
      write(6,*)
      write(6,*)'*************profile data from Xfoil:'
      write(6,*)'  k=        inc=         Cz=         Cx=      Cm,ac= '
      do 3 k=1,kx
         kp=k+1
         km=k-1
         prod=1.
         if(k.gt.1.and.k.lt.kx)then
         dcxm=((cx(k)-cx(km))*(cz(kp)-cz(k))*(cz(k)+cz(kp)-2.*cz(km))
     &        -(cx(kp)-cx(k))*(cz(k)-cz(km))**2)
     &        /((cz(kp)-cz(k))*(cz(kp)-cz(km))*(cz(k)-cz(km)))
         dcxp=((cx(k)-cx(kp))*(cz(km)-cz(k))*(cz(k)+cz(km)-2.*cz(kp))
     &        -(cx(km)-cx(k))*(cz(k)-cz(kp))**2)
     &        /((cz(km)-cz(k))*(cz(km)-cz(kp))*(cz(k)-cz(kp)))
         prod=dcxm*dcxp
         endif
         if(prod.lt.-eps.and.(kxtrm(km).ne.0.or.kxtrm(kp).ne.0))then
            write(6,*)'bad data distribution: interpolate a new data'
         endif
         dinc=inc(k)
         inc(k)=degrad*inc(k)
         write(6,1001)k,dinc,cz(k),cx(k),cq(k)
         write(14,*)cx(k),cz(k),cq(k),dinc
 3    continue
      write(6,*)
      write(6,*)'******extrema pointer:'
      write(6,*)'kxtrm(k)=',(kxtrm(k),k=1,kx)
      write(6,*)
 4    continue
c*****dimensionless variables
      write(6,*)'************************'
      write(6,*)'dimensionless variables:'
      write(6,*)'                      Y=0.5*B*y'
      write(6,*)'                      C=0.5*B*c'
      write(6,*)'                      A=B*C*am'
      write(6,*)'                      D=C*d'
      write(6,*)'                      W=U*w'
      write(6,*)'                   GAMA=0.5*U*B*g'
      write(6,*)'                   LIFT=0.5*RHO*U**2*A*Cl'
      write(6,*)'                   DRAG=0.5*RHO*U**2*A*Cd'
      write(6,*)'                 MOMENT=0.5*RHO*U**2*A*Cav*Cmo'
      write(6,*)'               REYNOLDS=RHO*U*C/AMU'
      write(6,*)'                     Fz=0.5*RHO*U**2*A*fz'
      write(6,*)'                     Mf=0.25*RHO*U**2*A*B*cmf'
      write(6,*)'                     Mt=0.5*RHO*U**2*A*Cav*cmt'
c*****read in data
      read(15,*)jx
      if(jx.gt.jxx)then
         write(6,*)'jx=',jx,'>jxx=',jxx,', change dimension:exiting!'
         stop
      endif
      jx2=jx/2
      is=mod(jx,2)
      si=is
      read(15,*)omega
      read(15,*)avis
      read(15,*)B
      read(15,*)cxm
      read(15,*)dm
      read(15,*)tmd
      read(15,*)irect
      read(15,*)am
      read(15,*)alphad
      read(15,*)Rho
      read(15,*)Vinf
      read(15,*)Amu
      tm=degrad*tmd
      arm=2./(am*cxm)
      am=2.*cxm*am
      cavm=cxm
      if(irect.eq.0)then
         cavm=8.*cxm/(3.*pi)
      endif
      alpha=degrad*alphad
      Dyn=0.5*Rho*Vinf**2
      Re=0.5*Rho*Vinf*B*cxm/Amu
c*****initialization
c*****mesh, geometry and flow
      write(6,*)
      write(6,*)'******point distribution, geometry and flow:'
      write(6,*)'read file y/n=1/0 ?'
      read(5,*)inprandtl
      write(6,1002)
      dtet=pi/(jx-1)
      dum=eps
      dumav=0.
      jm=1
      do 5 j=1,jx
         tetj=(j-1)*dtet
         yj=-cos(tetj)
         y(j)=yj
         etaj=-cos(tetj+.5*dtet)
         eta(j)=etaj
         c(j)=cxm*sin(tetj)
c         if(abs(y(j)).gt.0.63)then
c            c(j)=cxm*(1.5331-1.1331*abs(y(j)))
c         else
c            c(j)=cxm*(1.-0.3091*abs(y(j)))
c         endif
         if(irect.eq.1)then
            c(j)=cxm
         endif
         dum=dum+0.5*(c(jm)+c(j))*(y(j)-y(jm))
         dumav=dumav+(0.5*(c(jm)+c(j)))**2*(y(j)-y(jm))
         dem(j)=dm
         t(j)=tm
         g(j)=0.
         w(j)=0.
         at(j)=0.
         a0(j)=-pi*dem(j)
         a1(j)=0.
         b0(j)=2.*pi*(2.*dem(j))
         b1(j)=2.*pi
         if(inprandtl.eq.1)then
            read(16,*)y(j),eta(j),c(j),t(j),dem(j),g(j),w(j)
            at(j)=alpha+atan2(w(j),1.)+t(j)
         endif
         write(6,1003)y(j),eta(j),c(j),t(j),dem(j),g(j),w(j)
         cle=-0.25*c(j)-0.
         cte=0.75*c(j)-0.
         write(17,*)y(j),cle,cte,dem(j),t(j)
         jm=j
 5    continue
      eta(jx)=eta(jx-1)
      am=dum
      cavm=dumav/dum
      write(6,*)'***************'
      write(6,*)'numerical data:'
      write(6,*)'         number of points jx=',jx
      write(6,*)'                       omega=',omega
      write(6,*)'  viscosity coefficient avis=',avis
      write(6,*)'main wing data:'
      write(6,*)'                 wing span B=',B,' (m)'
      write(6,*)'               maximum chord=',cxm,' (ref. B/2)'
      write(6,*)'   average aerodynamic chord=',cavm,' (ref. B/2)'
      write(6,*)'      relative camber height=',dm,' (ref. C)'
      write(6,*)'          wing setting angle=',tm,' (rd) ='
     &     ,tmd,' (deg)'
      write(6,*)'  rectangular wing y/n(=1/0)=',irect
      write(6,*)'                   wing area=',am,' (ref. B**2/4)'
      write(6,*)'           wing aspect ratio=',arm
      write(6,*)'   geometric angle of attack=',alpha,' (rd) ='
     &     ,alphad,' (deg)'
      write(6,*)'air data:'
      write(6,*)'             air density Rho=',Rho,' (kg/m**3)'
      write(6,*)'          wind velocity Vinf=',Vinf,' (m/s)'
      write(6,*)'       dynamic viscosity Amu=',Amu,' (kg/(m*s))'
      write(6,*)'        dynamic pressure Dyn=',Dyn,' (Pa)'
      write(6,*)'       reference Reynolds Re=',Re
c*****iterations
      write(6,*)'title: airfoil, thickness, camber, Re <72 char.'
      read(5,1000)title
      write(25,1000)title
 100  continue
      write(6,*)'enter first integer alfa value=? or -100(exit)' 
     &     ,' or 360(given alpha)'
      read(5,*)ialfi
      if(ialfi.le.-100)goto 500
      write(6,*)'enter last integer alfa value=? or same as previous'
      read(5,*)ialff
      write(6,*)'enter integer alfa increment=? or 1'
      read(5,*)idalf
      if(idalf.eq.0)idalf=1
      do 400 ial=ialfi,ialff,idalf
      ialfd=ial
      if(ialfi.ne.360)alpha=degrad*ialfd
      alphad=alpha/degrad
      iter=0
      do 200 it=1,2001
      iter=iter+1
      if(ipolar.ne.1)goto 9
c     search for point on polar and polar coefficients
      do 8 j=1,jx
         atj=at(j)
         mj=1
         prod=atj+.5*pi
         do 6 k=2,kx-1
            prod=prod*(atj-inc(k))
            if(prod.le.eps)goto 7
            prod=1.
            mj=k
 6      continue
 7      continue
         b1(j)=(cz(mj+1)-cz(mj))/(inc(mj+1)-inc(mj))
         b0(j)=cz(mj)-b1(j)*inc(mj)
         czj=b0(j)+b1(j)*atj
         if(kxtrm(mj).ne.0)mj=mj+1
         if(mj.lt.2)mj=2
         if(mj.gt.kx-1)mj=kx-1
         m(j)=mj
         dcxdcz=(cx(mj)-cx(mj-1))/(cz(mj)-cz(mj-1))
         d2cxdcz2=((cx(mj+1)-cx(mj))*(cz(mj)-cz(mj-1))
     &        -(cx(mj)-cx(mj-1))*(cz(mj+1)-cz(mj)))
     &     /((cz(mj+1)-cz(mj))*(cz(mj+1)-cz(mj-1))*(cz(mj)-cz(mj-1)))
         cxj=cx(mj-1)+dcxdcz*(czj-cz(mj-1))
     &        +d2cxdcz2*(czj-cz(mj-1))*(czj-cz(mj))
c         write(6,*)'czj=',czj,' cxj=',cxj
         a0(j)=cx(mj-1)-dcxdcz*cz(mj-1)+d2cxdcz2*cz(mj-1)*cz(mj)
         a1(j)=dcxdcz-d2cxdcz2*(cz(mj-1)+cz(mj))
         a2(j)=d2cxdcz2
         l(j)=czj
         d(j)=cxj
         q(j)=cq(mj)+(cq(mj+1)-cq(mj))*(czj-cz(mj))/(cz(mj+1)-cz(mj))
 8    continue
 9    continue
      jdx=0
      dgx=0.
      do 11 j=2,jx-1
         sum=0.
         do 10 k=1,jx-1
            sum=sum+(g(k+1)-g(k))/(y(j)-eta(k))
 10      continue
         wj=-sum/(4.*pi)
         atj=alpha+atan2(wj,1.)
         attj=atj+t(j)
         czj=b0(j)+b1(j)*attj
         if(ipolar.eq.1)then
            vis=0.
            if(m(j).ge.mxtrm)then
               vis=-c(j)*b1(j)
     &              *(1./(y(j)-eta(j-1))-1./(y(j)-eta(j)))/(16.*pi)
               if(vis.lt.eps)vis=0.
            endif
         endif
         dg(j)=omega*(0.5*c(j)*czj-g(j)
     &        +(avis+vis)*(g(j+1)-2.*g(j)+g(j-1)))
     &        /(1.+c(j)*b1(j)*(1./(y(j)-eta(j-1))-1./(y(j)-eta(j)))
     &        /(8.*pi)+2.*(avis+vis))
         w(j)=wj
         at(j)=attj
         g(j)=g(j)+dg(j)
         if(abs(dgx).lt.abs(dg(j)))then
            dgx=dg(j)
            jdx=j
         endif
 11   continue
      w(1)=w(2)
     &     +(w(3)-w(2))*(y(1)-y(2))/(y(3)-y(2))
      w(jx)=w(jx-1)
     &     +(w(jx-2)-w(jx-1))*(y(jx)-y(jx-1))/(y(jx-2)-y(jx-1))
      at(1)=at(2)
     &     +(at(3)-at(2))*(y(1)-y(2))/(y(3)-y(2))
      at(jx)=at(jx-1)
     &     +(at(jx-2)-at(jx-1))*(y(jx)-y(jx-1))/(y(jx-2)-y(jx-1))
      if(abs(dgx).lt.eps)goto 300
 200  continue
      write(6,*)'************NOT CONVERGED!!!!!'
 300  continue
      cl=0.
      cmo=0.
      fz(1)=0.
      fz(jx)=0.
      cmf(1)=0.
      cmf(2)=0.
      cmf(jx)=0.
      cmt(1)=0.
      cmt(jx)=0.
      do 12 j=2,jx-1
         cl=cl+g(j)*(eta(j)-eta(j-1))
         if(ipolar.ne.1)then
            q(j)=a0(j)+a1(j)*at(j)
            cmo=cmo+(c(j)**2*q(j)
     &           -0.25*cxm*c(j)*(b0(j)+b1(j)*at(j)))*(eta(j)-eta(j-1))
            if(j.eq.jx2+1)then
               cmt(j)=-cmt(j-1)+(1.-si)*c(j)**2*q(j)*(eta(j)-eta(j-1))
            else
               cmt(j)=cmt(j-1)+c(j)**2*q(j)*(eta(j)-eta(j-1))
            endif
         else
            cmo=cmo
     &           +(c(j)**2*q(j)-0.25*cxm*c(j)*l(j))*(eta(j)-eta(j-1))
         endif
         if(j.eq.jx2+1)then
            fz(j)=-fz(j-1)+(1.-si)*g(j)*(eta(j)-eta(j-1))
            cmt(j)=-cmt(j-1)+(1.-si)*c(j)**2*q(j)*(eta(j)-eta(j-1))
         else
            fz(j)=fz(j-1)+g(j)*(eta(j)-eta(j-1))
            cmt(j)=cmt(j-1)+c(j)**2*q(j)*(eta(j)-eta(j-1))
         endif
         cmf(j+1)=cmf(j)-fz(j)*(eta(j+1)-eta(j))
 12   continue
      cl=0.5*arm*cl
      cmo=0.25*arm*cmo/cavm
      sum=0.
      sum0=0.
      sum1=0.
      sum2=0.
      do 13 j=1,jx
         jm=j-1
         if(j.eq.1)jm=1
         czj=b0(j)+b1(j)*at(j)
         sum=sum+g(j)*w(j)*(eta(j)-eta(jm))
         if(ipolar.ne.1)then
            rey=c(j)*Re
            if(rey.lt.1000.)rey=1000.
            Cd0=1.328/sqrt(rey)
            if(rey.gt.1.0E5)then
               Cd0=.072/rey**.2
            endif
            Cd0=2.*Cd0
            sum0=sum0+2.*Cd0*(eta(j)-eta(jm))/arm
            l(j)=b0(j)+b1(j)*at(j)
            d(j)=Cd0
         else
            sum0=sum0+a0(j)*c(j)*(eta(j)-eta(jm))
            sum1=sum1+a1(j)*czj*c(j)*(eta(j)-eta(jm))
            sum2=sum2+a2(j)*czj*czj*c(j)*(eta(j)-eta(jm))
         endif
 13   continue
      cdi=-0.5*arm*sum
      cdv=0.25*arm*(sum0+sum1+sum2)
      em=cl*cl/(pi*arm*cdi)
      cd=cdi+cdv
      write(6,*)'alphad=',alphad,' iter=',iter,' dgx=',dgx,' jdx=',jdx
      write(6,*)'CL=',cl,' CD=',cd,' CM,o=',cmo,' eta=',em
      write(6,*)'m(j)=',(m(j),j=1,jx)
      write(25,1011)alphad,cl,cd,cmo
 400  continue
      goto 100
 500  continue
c*****distributions
      write(6,1004)
      do 14 j=1,jx
         tetj=(j-1)*dtet
         ge=4.*(0.00675*sin(tetj)+0.0*sin(2.*tetj))
         ge=4.*0.1666666*alpha*sin(tetj)
         ge=0.5*cl*c(j)
         we=-0.00675
         we=-0.16666666*alpha
         we=-cl*cxm/8.
         write(6,1003)y(j),c(j),t(j),dem(j),g(j),w(j),l(j),d(j)
         write(18,*)y(j),g(j),w(j)
         write(19,*)y(j),g(j),ge
         write(20,*)y(j),w(j),we
         write(21,*)y(j),l(j)
         write(22,*)y(j),at(j)
         write(23,*)d(j),l(j)
         write(24,*)y(j),eta(j),c(j),t(j),dem(j),g(j),w(j)
 14   continue
c*****force and moment
      do 15 j=1,jx
         fz(j)=0.5*arm*fz(j)
         cmf(j)=0.5*arm*cmf(j)
         cmt(j)=0.25*arm*cmt(j)/cavm
         write(26,*)y(j),cmf(j)
         write(27,*)eta(j),fz(j)
         write(28,*)eta(j),cmt(j)
 15   continue
      write(6,*)
      write(6,*)'********'
      write(6,*)'results:'
      write(6,1005)cdi,em
      write(6,1006)cdv
      write(6,1007)cd,cl
      write(6,1008)cmo
      write(6,1009)-cmf(jx2+1)
      write(6,1010)-cmt(jx2+1)
      write(6,*)
c      write(6,*)'g(j)=',(g(j),j=1,jx)
c      write(6,*)'w(j)=',(w(j),j=1,jx)
      write(6,*)'at(j)=',(at(j),j=1,jx)
      write(6,*)
c*****files
      write(6,*)'******data file:'
      write(6,*)'polarbl.dat           :profile polar from Xfoil'
      write(6,*)'prandtline.data       :geometric data for wing'
      write(6,*)'prandtline.in         :restart file'
      write(6,*)'******output files:'
      write(6,*)'polarbl.cdcl          :profile polar cd,cl,cq,inc'
      write(6,*)'prandtline.ycdt       :y,c(y),d(y) and t(y)'
      write(6,*)'prandtline.ygw        :y,g(y) and w(y)'
      write(6,*)'prandtline.yg         :y,g(y) circulation gama'
      write(6,*)'prandtline.yw         :y,w(y) downwash'
      write(6,*)'prandtline.ycl        :y,cl(y)local lift coefficient' 
      write(6,*)'prandtline.yat        :y,at(y) angle of attack'
      write(6,*)'prandtline.cdcl       :cd(y),cl(y) polar'
      write(6,*)'prandtline.clcdcq     :alpha,CL,CD,CM,o'
      write(6,*)'prandtline.out        :save results for restart'
      write(6,*)'prandtline.ycmf       :y,cmf(y)'
      write(6,*)'prandtline.yfz        :y,fz(y)'
      write(6,*)'prandtline.ycmt       :y,cmt(y)'
 1000 format(18a4)
 1001 format(1x,i4,4x,f8.4,4x,f8.4,4x,f8.4,4x,f8.4)
 1002 format(7x,'y(j)=',5x,'eta(j)=',7x,'c(j)=',7x,'t(j)='
     &     ,7x,'d(j)=',7x,'g(j)=',7x,'w(j)=')
 1003 format(8f12.4)
 1004 format(7x,'y(j)=',7x,'c(j)=',7x,'t(j)=',7x,'d(j)=',7x,'g(j)=',
     &     7x,'w(j)=',7x,'  Cl=',7x,'  Cd=')
 1005 format(' inviscid contribution: CDi=',f10.6
     &     ,'    Oswald efficiency e=',f10.4)
 1006 format('  viscous contribution: CDv=',f10.6)
 1007 format('        global results:  CD=',f10.6
     &     ,'    lift coefficient CL=',f10.4)
 1008 format('                             '
     &     ,'pitching moment coefficient CM,o=',f10.4)
 1009 format('                             '    
     &     ,'  root bending moment coef. CM,x=',f10.4)
 1010 format('                             '
     &     ,'  root torsion moment coef. CM,y=',f10.4)
 1011 format(1x,f10.4,4x,f10.4,4x,f10.5,4x,f10.4)
c
      read(5,1000)press
c
 
      end