c cinter2.f c c
      program interpolate
      implicit none
c
c     Demonstation of 2 simple interpolation methods for
c     smoothly joining results from 2 disjoint regions
c
c     John Mahaffy 2/3/96
c
c    Open up two files in your current directory to receive results
c
      open(10,file='nulin')
      open(11,file='nucube')
c
c   First a linear interpolation
c
      call nulin
c
c   Next a cubic interpolation
c
      call nucube
c
c   If you are in Hammond Lab, change this subroutine to plot results
c
      call plotit
c
      stop
      end
c
      subroutine nulin
      implicit none
      real Nulam,Nuturb,Pr,Re,w,Nu
      real Re1,Re2,Nu1,Nu2
c
c   Demonstration of simple linear end-point interpolation
c   Calculation of Nusselt Number
c
c
c     John Mahaffy 2/3/96
c
c   Nu   -  Nusselt Number
c   Nulam - laminar Nusselt Number
c   Nuturb - Turbulent Nusselt Number
c   Re  -  Reynolds number
c   w   -  interpolation weighting function
c
      parameter (Pr=1.0,Nulam=4.0,Re1=640,Re2=2000.)
c
      Nu1=Nulam
      Nu2=Nuturb(Re2,Pr)
      Re=0.0
c
c    Evaluate the Nusselt number over the range 0.<=Re<=3000.
c
      do 100 while (Re.le.3000.)
        if(Re.lt.Re1) then
           Nu = Nulam
c
        else if (Re.gt.Re2)  then
c
           Nu=Nuturb(Re,Pr)
c
        else
c
           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
c
c   Demonstration of simple cubic end-point interpolation between two
c   correlations for the Nusselt number
c
c
c     John Mahaffy 2/3/96
c
c   Nu   -  Nusselt Number
c   Nulam - laminar Nusselt Number
c   Nuturb - Turbulent Nusselt Number
c   Re  -  Reynolds number
c   SRe -  Scaled Reynolds number, ranges from 0 to 1 over interpolation range
c   w   -  interpolation weighting function
c
      parameter (Pr=1.0,Nulam=4.0,Re1=640,Re2=2000.)
c
c
c     Prepare some constants necessary to generate polynomial coefficients
c     giving continuity with the bounding functions, first derivatives of
c     those functions.
c
      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
c
c    Next Evaluate the Nusselt number over the range 0.<=Re<=3000.
c
      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)
c
c    Calculate a Turbulent heat transfer coefficient
c    obtained from a Dittus-Boelter correlation
c
c     John Mahaffy   2/2/96
c
      implicit none
      real Re,Pr,Nuturb
c
c   Nuturb - Turbulent Nusselt number (Dittus-Boelter correlation)
c   Re  -  Reynolds number
c   Pr  -  Prandl number
c
c     Nuturb=0.023*Re**.8*Pr**0.4
c
      Nuturb=0.023*exp(log(Re)*0.8+log(Pr)*0.4)
      return
      end
      function NuDeriv(Re,Pr)
c
c    Calculate the Derivative of Turbulent Nusselt Number with respect
c    to Reynolds number.
c    Based on the Function Nuturb
c
c     John Mahaffy   2/2/96
c
      implicit none
      real Re,Pr,NuDeriv
c
c   NuDeriv - Derivative of Nusselt number
c   Re  -  Reynolds number
c   Pr  -  Prandl number
c
c     NuDeriv=0.0184*Re**(-.2)*Pr**0.4
c
c 
      NuDeriv=0.0184*exp(-log(Re)*0.2+log(Pr)*0.4)
c 
      return
      end
      subroutine plotit
c
c    This subroutine is for those of you who are getting bored.  At this
c    stage of the class you will not be held responsible new material here.
       logical fexist

c      fexist -  flag indicating whether or not the gnuin file exists

c   First check to see if a file containing commands to the graphic program
c   exists in your local directory
c
      inquire (file='gp-int2',exist=fexist)
c
c   If it doesn't then issue a unix command to copy it from my directory
c   the subroutine "system" is not a standard part of Fortran, but is
c   provided by some Unix systems to give a crude form of what is called
c   a "system call".  Any computer will supply an equivalent command, and
c   other functions or subroutines that give more direct access to system
c   support.
c
      if (.not.fexist) then
         call system('cp ~jhm/201/gp-int2 .')
         print *, ' Graphics input file copied'
      endif
c
c   Close the data files so they are ready for reading by gnuplot
c
      close ( 11 )
      close ( 10 )
c
c   If you are sitting at a Hammond lab work station remove the "c" in
c   column 1 of the following "call system" to plot your results.  If you
c   are running from NCSA/BYU Telnet on a Mac, you can see plots by
c   making the same change, and by following instructions in "gp-int2"
c   to change that file for Tektronics emulation in gnuplot.
c
c      call system ('gnuplot gp-int2')
c
      return
      end
c
c c