program interpolate
      implicit none
      open(10,file='nulin')
      open(11,file='nucube')
      call nulin
      call nucube
      call plotit
      stop
      end
      subroutine nulin
      implicit none
      real Nulam,Nuturb,Pr,Re,w,Nu
      real Re1,Re2,Nu1,Nu2
      parameter (Pr=1.0,Nulam=4.0,Re1=640,Re2=2000.)
      Nu1=Nulam
      Nu2=Nuturb(Re2,Pr)
      Re=0.0
      do 100 while (Re.le.3000.)
        if(Re.lt.Re1) then
           Nu = Nulam
        else if (Re.gt.Re2)  then
           Nu=Nuturb(Re,Pr)
        else
           w = (Re - Re1)/(Re2 - Re1)
           Nu = (1.-w)*Nu1 + w*Nu2
        endif
        write(10,*)Re,Nu
        Re=Re+10.
  100 continue
      return
      end
      subroutine nucube
      implicit none
      real Nuturb,Nulam,Pr,Re,w,SRe,Re1,Re2, Nu
      real Nu1,Nu2,dNudRe1, dNudRe2, DwdSRe1, DwdSRe2, fac, a1,a2,a3
      real NuDeriv
      parameter (Pr=1.0,Nulam=4.0,Re1=640,Re2=2000.)
      Nu2=Nuturb(Re2,Pr)
      dNudRe2=NuDeriv(Re2,Pr)
      Nu1=Nulam
      dNudRe1=0.
      fac=(Re2-Re1)/(Nu2-Nu1)
      dwdSRe1=dNudRe1*fac
      dwdSRe2=dNudRe2*fac
      a3=dwdSRe1+dwdSRe2-2
      a2=3-2*dwdSRe1-dwdSRe2
      a1=dNudRe1
      Re=0.
      do 200 while (Re.le.3000.)
        if(Re.lt.Re1) then
           Nu=Nulam
        else if (Re.gt.Re2) then
           Nu=Nuturb(Re,Pr)
        else
           SRe=(Re-Re1)/(Re2-Re1)
           w= a3*SRe**3 + a2*SRe**2 + a1*SRe
           Nu= (1.-w)*Nu1+w*Nu2
        endif
        write(11,*)Re,Nu
        Re=Re+10.
  200 continue
      return
      end
      function Nuturb(Re,Pr)
      implicit none
      real Re,Pr,Nuturb
      Nuturb=0.023*exp(log(Re)*0.8+log(Pr)*0.4)
      return
      end
      function NuDeriv(Re,Pr)
      implicit none
      real Re,Pr,NuDeriv
      NuDeriv=0.0184*exp(-log(Re)*0.2+log(Pr)*0.4)
      return
      end
      subroutine plotit
       logical fexist


      inquire (file='gp-int2',exist=fexist)
      if (.not.fexist) then
         call system('cp ~jhm/201/gp-int2 .')
         print *, ' Graphics input file copied'
      endif
      close ( 11 )
      close ( 10 )
c      call system ('gnuplot gp-int2')
      return
      end