
c     FRICTION v3

c     skin friction and form drag estimation program

c     by W. Mason,   October 21, 1989
c     Department of Aerospace and Ocean Engineering
c     Virginia Tech, Blacksburg, VA 24061
c     mason@aoe.vt.edu

c     version 3: mod to include laminar and Re as input, Feb. 24, 1992
c                minor mods: Nov.  7, 1992
c                            Feb. 26, 1994
c                            Sep. 13, 1996
c                            Dec. 14, 2000 by Andy Ko
c                            Mar. 14, 2001 by Andy Ko
c    						   Jan. 28, 2008 (Torenbeek form factor adopted)
c     Mod notes:
c     1. All output from the screen is also written to the output file. Also,
C        a pause statement was included in the end to prevent command window
C        from closing when program terminates. - Andy Ko 12/14/00
C     2. User is now able to change the working directory. Default is Z:\
C        which is the network drive. This can be changed to c:\ or any other
C        directory if desired. Also, a title screen was added. - Andy Ko 3/14/01 

      USE DFPORT
      DIMENSION     TITLE(15),   COMP(15,4),  SWET(15),    REFL(15),
     1              FF(15),      TC(15),      ICODE(15),   xmetbl(200),
     2              CFI(15),     R(15),       CDFCAS(200), CFSW(15),
     3              CFSWFF(15),  TRANS(15),   RNCAS(200),  CDICAS(200),
     4              alttbl(200)

      character*4   TITLE,COMP
      character*15  datset,outputfile
	CHARACTER*50 newdir
      INTEGER*4 istatus

      A      = 8.33333E-08
      IWRIT  = 6
      IREAD  = 1
	IFILE  = 8

C  Title screen - Andy Ko 3/14/01

      WRITE (*,*) '****************************************************'
	WRITE (*,*) '     VT Aerospace Aircraft Design Software Series'
	WRITE (*,*) '****************************************************'
	WRITE (*,*)
      WRITE (*,*) 'FRICTION v3.1'
      WRITE (*,*)
      WRITE (*,*) 'skin friction and form drag estimation program'
      WRITE (*,*)
      WRITE (*,*) 'by W. Mason,   October 21, 1989'
      WRITE (*,*) 'Department of Aerospace and Ocean Engineering'
      WRITE (*,*) 'Virginia Tech, Blacksburg, VA 24061'
      WRITE (*,*) 'whmason@vt.edu'
      WRITE (*,*) 
	WRITE (*,*) '****************************************************'
      WRITE (*,*)

C  Set the input/output directory. Default is Z:\ - Andy Ko 3/14/01
C	Change the above comment to make default the current directory-Lowe 2/1/2006
      
	istatus=enoent+1
C		WRITE(*,*) enoent	Commented by Lowe 2/1/2006
	IF (istatus .eq. enoent) THEN
	     WRITE(*,*) 'Cannot change directory to :\'
	     WRITE(*,*) 'Please log into the network'
	     Pause 'Press the ENTER key to exit'
	     STOP
      ENDIF
C		Added by Lowe 2/1/2006:
      WRITE (*,*) 'Current directory is the default'
	WRITE (*,*) 'If desired, please enter a new directory name. '
	WRITE (*,*) 'Type "default" for the current directory' 
	READ (*,*) newdir
	IF (newdir .eq. 'default') THEN
	   WRITE (*,*) 'Current directory selected'
      ELSE
	   istatus = CHDIR(newdir)
         IF (istatus .eq. enoent) THEN
            WRITE(*,*) 'The directory does not exist'
	      WRITE(*,*) 'The directory is set to the default'
         ELSEIF (istatus .eq. enotdir) THEN
            WRITE(*,*) 'Directory could not be changed'
	      WRITE(*,*) 'The directory is set to the default'
         ELSE
	      WRITE(*,*) 'directory has been changed to ',newdir
         ENDIF
      ENDIF

2     WRITE(6,5)
      READ (5,7) datset
      OPEN (unit = IREAD,file=datset,STATUS='old',ERR=500)
	WRITE(6,6)
	READ (5,7) outputfile

	OPEN (unit = IFILE,file=outputfile)

      READ (IREAD,10) TITLE
      READ (IREAD,15) SREF,SCALE,FNCOMP,Finmd
      ncomp  = fncomp
      inmd   = finmd
   
      ASCALE = SCALE*SCALE
      SREF   = SREF/ASCALE

      WRITE(iwrit,20) TITLE,SREF,SCALE,NCOMP,inmd
	WRITE(ifile,20) TITLE,SREF,SCALE,NCOMP,inmd
      WRITE(iwrit,25)
	WRITE(ifile,25)

c     read in the component geometry and lam/turb. switch

      sum = 0.0

      DO 40 I   = 1,NCOMP
      READ (IREAD,50) (COMP(I,J),J=1,4),
     1                 SWET(I),REFL(I),TC(I),FICODE,Ftrans
      ICODE(I)  = FICODE
      TRANS(I)  = Ftrans
      SWET(I)   = SWET(I)/ASCALE
      sum       = sum + swet(i)
      REFL(I)   = REFL(I)/SCALE
      FF(I)     = 1.00

cwhm	Torenbeek form factor adopted, Jan. 28. 2006
      IF(ICODE(I).EQ.0) then
                        FF(I) = 1.0 + 2.7*TC(I) + 100.*TC(I)**4
                        else
                        FF(I) = 1.0 + 1.5*TC(I)**1.5 + 7.*TC(I)**3
                        endif
      WRITE(iwrit,60) (COMP(I,J),J=1,4),SWET(I),REFL(I),TC(I),
     1                 ICODE(I),FF(I),trans(i)
	WRITE(ifile,60) (COMP(I,J),J=1,4),SWET(I),REFL(I),TC(I),
     1                 ICODE(I),FF(I),trans(i)
   40 CONTINUE

      write(iwrit,170) sum
	write(ifile,170) sum

      TWTAW     = 1.0
      XME       = 0.2
      alt       = 0.0
      JJ        =   0

c     read in the flight condition

   65 READ (IREAD,70, end = 200) XME,xinput
      if(xme .le. 0.0) go to 200

      if(inmd .eq. 0) then 
               alt = 1000.0*xinput
               kd  = 1
               call stdatm(ALT,t,p,rho,a,xmu,ts,rr,pp,rmr,qm,kd,kk)
               if (kk.ne.0) then 
			    write(iwrit,80)
	            write(ifile,80)
	         endif
               RN  = rmr*XME

                     else

               rn  = xinput*10.**6

                     endif

c     compute the skin friction

      JJ        = JJ + 1

      DO 100 I  = 1, NCOMP
      R(I)      = RN*REFL(I)
      Rec       = R(I)*trans(i)
      CFLAM     = 0.0
      CFTURBC   = 0.0

      call turbcf(R(I),XME,TWTAW,CFTURBL)

      if(trans(i) .gt. 0.0) then
                            call turbcf(Rec,XME,TWTAW,CFTURBC)
                            call lamcf (Rec,XME,TWTAW,CFLAM)
                            endif

      CFI(I)    = CFTURBL - trans(i)*(CFTURBC - CFLAM)

   90 CFSW(I)   = CFI(I)*SWET(I)
      CFSWFF(I) = CFI(I)*SWET(I)*FF(I)
  100 CONTINUE
 
      sum       = 0.
      sum2      = 0.
      sum3      = 0.
      if (inmd .eq. 0) then
	   WRITE (iwrit,120) RN,alt,xme
         WRITE (ifile,120) RN,alt,xme
	endif         
      if (inmd .eq. 1) then 
	   write (iwrit,122) RN,xme
         write (ifile,122) RN,xme
      endif
      WRITE (iwrit,125)
	WRITE (ifile,125)

      DO 150 I  = 1, NCOMP
      CDCOMP    = CFSWFF(I)/SREF
      sum       = sum  + CFSW(I)
      sum2      = sum2 + CFSWFF(I)
      sum3      = sum3 + CDCOMP
      WRITE(iwrit,160) (COMP(I,J),J=1,4),R(I),
     1                  CFI(I),CFSW(I),CFSWFF(I),CDCOMP
  150	WRITE(ifile,160) (COMP(I,J),J=1,4),R(I),
     1                  CFI(I),CFSW(I),CFSWFF(I),CDCOMP

      WRITE(iwrit,180) sum,sum2,sum3
	WRITE(ifile,180) sum,sum2,sum3

      CD         = SUM/SREF
      CDFORM     = (SUM2-SUM)/SREF
      WRITE (iwrit,190) CD,CDFORM
	WRITE (ifile,190) CD,CDFORM

      CDICAS(JJ) = CD
      CDFCAS(JJ) = CDFORM
      RNCAS(JJ)  = RN
      xmetbl(jj) = xme
      alttbl(jj) = alt

      GO TO 65
 
c     print out a summary of the results
 
  200 CONTINUE

      if (inmd .eq. 0) then 
	   WRITE (iwrit,210)
	   WRITE (ifile,210)
	endif
      if (inmd .eq. 1) then 
	   write (iwrit,212)
	   write (ifile,212)
	endif

      DO 400 J = 1,JJ
      CDSUM = CDICAS(J)+CDFCAS(J)
      if (inmd .eq. 0) then
                WRITE(iwrit,420) J,XMETBL(J),ALTTBL(J),RNCAS(J),
     1                           CDICAS(J),CDFCAS(J),CDSUM
	          WRITE(ifile,420) J,XMETBL(J),ALTTBL(J),RNCAS(J),
     1                           CDICAS(J),CDFCAS(J),CDSUM
                       else
                WRITE(iwrit,422) J,XMETBL(J),RNCAS(J),
     1                           CDICAS(J),CDFCAS(J),CDSUM
	          WRITE(ifile,422) J,XMETBL(J),RNCAS(J),
     1                           CDICAS(J),CDFCAS(J),CDSUM
                       endif
  400 continue

      WRITE(iwrit,430)
	WRITE(ifile,430)

	close(unit=ifile)

	Pause 'Press ENTER key to terminate'
 
    5 FORMAT(' Enter name of data set:')
    6 FORMAT(' Enter name of output file')
    7 FORMAT(a15)
   10 FORMAT(15A4)
   15 FORMAT(4F10.5,I2)    
   20 FORMAT(//5X,'FRICTION - Skin Friction and Form Drag Program'/
     1 5x,'W.H. Mason, Department of Aerospace and Ocean Engineering'/
     2 5x,'Virginia Tech, Blacksburg, VA 24061 email:whmason@vt.edu'/
     1 5x,'version:   Jan. 28, 2006' //
     1 5X,'CASE TITLE:',15A4//
     2 5X,'SREF =',F12.5,2X,'MODEL SCALE =',F7.3,2X,
     3 'NO. OF COMPONENTS =',I2/
     4 5x,'input mode = ',i2,2x,
     5 '(mode=0: input M,h;  mode=1: input M, Re/L)'/)
   25 FORMAT(3x,' COMPONENT TITLE',4X,'SWET (FT2)',2X,'REFL(FT) ',
     1 2X,'TC',2X,'ICODE',1X,'FRM FCTR',1x,'FTRANS')
   50 FORMAT (4A4,4X,7F10.5,I1)
   60 FORMAT(4X,4A4,1X,F11.4,1x,f8.3,F8.3,2X,I2,2X,F7.4,2x,f7.4)
   70 FORMAT(4f10.3)
   80 format(/4x,'bad return from stdatm - check data'/)
  120 FORMAT(/5X,'REYNOLDS NO./FT =',E9.3,3x,'Altitude = ',f9.2,
     1 3x,'XME =',f7.3/)
  122 FORMAT(/5X,'REYNOLDS NO./FT =',E9.3,3x,'XME =',f7.3/)
  125 FORMAT(7X,'COMPONENT',8X,'RN',6X,'  CF ',4X,
     1          'CF*SWET',3X,'CF*SWET*FF',3X,'CDCOMP')
  160 FORMAT(4X,4A4,E10.3,2X,F7.5,F9.5,2X,F10.5,F10.5)
  170 FORMAT(/7X,'TOTAL SWET = ',F12.4)
  180 FORMAT(35X,'SUM =',F8.5,2X,F10.5,f10.5)
  190 FORMAT (/4X,'FRICTION DRAG: CDF = ',F7.5,
     1   11X,'FORM DRAG: CDFORM = ',F7.5/)
  210 FORMAT(/5X,'SUMMARY'//6X,'J',4x,'XME',5x,'Altitude',5X,
     1 ' RE/FT ',5X,' CDF  ', 2X,' CDFORM',3X,'CDF+CDFORM')
  212 FORMAT(//5X,'SUMMARY'//6X,'J',3x,'XME',5x,
     1 ' RE/FT ',5X,' CDF  ', 2X,' CDFORM',3X,'CDF+CDFORM')
  420 FORMAT(5X,I2,1X,F7.3,2x,E10.3,3x,E10.3,3F10.5)
  422 FORMAT(5X,I2,1X,F7.3,2x,E10.3,3F10.5)
  430 FORMAT(/5X,'END OF CASE'//)

      STOP
C
C   Error in opening file statement - Andy Ko 3/14/01
  500 write(*,*) "File cannot be opened or does not exist. Please try ag
     &ain."
	GOTO 2

      END

      subroutine lamcf(Rex,Xme,TwTaw,CF)
c
c     flat plate laminar skin friction routine
c
c     uses the Eckert Refereence Temperatrure Method
c     as described in White, Viscous Fluid Flow, 1974 ed., pg 589-590.  
c
c     see Boundary Layer Analysis Methods in Aerocal Methods
c         Pak #4, pages 4-39 to 4-41
c 
c     coded in FORTRAN, Feb. 24, 1992
c     by W.H. Mason 
c
c     Input:
c                  Rex   - Reynolds number
c                  Xme   - Mach Number
c                  TwTaw - Tw/Taw, the wall temperature ratio
c
c     values of gamma, Pr, and Te are set in the program
c
c     Output:
c                CF    - total skin friction coefficient

      G     = 1.4
      Pr    = .72
      R     = sqrt(Pr)
      TE    = 390.
      TK    = 200.

      TwTe  = TwTaw*(1. + R*(g - 1.)/2.*Xme**2)
      TstTe = 0.5 + 0.039*Xme**2 + 0.5*TwTe

      Cstar = sqrt(TstTe)*(1. + Tk/Te)/(TstTe + Tk/Te)

      CF  = 2.*.664*sqrt(cstar)/sqrt(Rex)

      RETURN

      END
 
      subroutine turbcf(Rex,xme,TwTaw,CF)
c
c     flat plate turbulent skin friction routine
c     
c     uses the Van Driest II Method, NASA TN D-6945
c 
c     coded in FORTRAN, October 21, 1989, by W.H. Mason
c
c     Input:       Rex   - Reynolds number (based on length)
c                  Xme   - Mach Number
c                  TwTaw - Tw/Taw, the wall temperature ratio
c
c     values of gamma, recovery factor and Te are set in the program
c
c     Output:      CF    - total skin friction coefficient

      epsmax = 0.2e-8
      G      = 1.4
      r      = 0.88
      Te     = 222.0

      xm    = (G - 1.)/2*xme**2
      TawTe = 1. + r*xm
      F     = TwTaw*TawTe
      Tw    = F * Te
      A     = SQRT(r*xm/F)
      B     = (1. + r*xm - F)/F
      denom = SQRT(4.*A**2 + B**2)
      Alpha = (2.*A**2 - B)/denom
      Beta  = B/denom
      Fc    = ((1.0 + SQRT(F))/2.0)**2
      If ( xme .gt. 0.1) Fc = r*xm/(ASIN(Alpha) + ASIN(Beta))**2

      Xnum   = (1. + 122./Tw*10**(-5/Tw))
      Denom  = (1. + 122./Te*10**(-5/Te))
      Ftheta = SQRT(1./F)*(Xnum/Denom)
      Fx     = Ftheta/Fc

      RexBar = Fx * Rex

      Cfb    = 0.074/RexBar**0.20

      iter   = 0
   10 iter   = iter + 1
      if (iter .gt. 200) write(6,100) eps, epsmax
	if (iter .gt. 200) go to 20
      Cfo    = Cfb
      Xnum   = 0.242 - SQRT(Cfb)*ALOG10(RexBar*Cfb)
      Denom  = 0.121 + SQRT(Cfb)/ALOG(10.)
      Cfb    = Cfb*(1.0 + Xnum/Denom)
      eps    = abs(Cfb - Cfo)
      if (eps .gt. epsmax) go to 10
   20	continue
      CF     = Cfb/Fc

  100 format(/3x,'did not converge in turbcf'/)

       RETURN
       END

      subroutine stdatm(z,t,p,r,a,mu,ts,rr,pp,rm,qm,kd,kk)
c
c   *********** 1976 STANDARD ATMOSPHERE SUBROUTINE **********
c
c     Mason's BASIC program, converted to FORTRAN - Sept. 1, 1989
c
c     kd -   = 0 - metric units
c           <> 0 - English units
c
c     kk - 0 - good return
c          1 - error: altitude out of table,
c              do not use output
c
c     z  - input altitude, in feet or meters (depending on kd)
c
c     output:
c
c     t  - temp.
c     p  - pressure
c     r  - density
c     a  - speed of sound
c     mu - viscosity
c     
c     ts - t/t at sea level
c     rr - rho/rho at sea level
c     pp - p/p at sea level
c
c     rm - Reynolds number per Mach per unit of length
c     qm - dynamic pressure/Mach^2
c
      real k, h, mu, ml
      KK = 0
      K  = 34.163195
      C1 = 3.048E-04
      IF (KD .eq. 0) goto 1240
      TL = 518.67
      PL = 2116.22
      RL = .0023769
      AL = 1116.45
      ML = 3.7373E-07
      BT = 3.0450963E-08
      GOTO 1260
 1240 TL = 288.15
      PL = 101325
      RL = 1.225
      C1 = .001
      AL = 340.294
      ML = 1.7894E-05
      BT = 1.458E-06
 1260 H = C1 * Z / (1 + C1 * Z / 6356.766)
      IF (H .gt. 11.0) goto 1290
      T = 288.15 - 6.5 * H
      PP = (288.15 / T) ** ( - K / 6.5)
      GOTO 1420
 1290 IF (H .gt. 20.0) goto 1310
      T = 216.65
      PP = .22336 *  EXP ( - K * (H - 11) / 216.65)
      GOTO 1420
1310  IF (H .gt. 32.0) goto 1330
      T = 216.65 + (H - 20)
      PP = .054032 * (216.65 / T) ** K
      GOTO 1420
1330  IF (H .gt. 47.0) goto 1350
      T = 228.65 + 2.8 * (H - 32)
      PP = .0085666 * (228.65 / T) ** (K / 2.8)
      GOTO 1420
1350  IF( H .gt. 51.0) goto 1370
      T = 270.65
      PP = .0010945 *  EXP ( - K * (H - 47) / 270.65)
      GOTO 1420
1370  IF (H .gt. 71.) goto 1390
      T = 270.65 - 2.8 * (H - 51)
      PP = .00066063 * (270.65 / T) ** ( - K / 2.8)
      GOTO 1420
1390  IF (H .gt. 84.852) THEN 
                              kk = 1
                              write(6,200) H
                              return
                         END IF
      T = 214.65 - 2 * (H - 71)
      PP = 3.9046E-05 * (214.65 / T) ** ( - K / 2)
1420  RR = PP / (T / 288.15)
      MU = BT * T**1.5 / (T + 110.4)
      TS = T / 288.15
      A  = AL *  SQRT (TS)
      T  = TL * TS
      R  = RL * RR
      P  = PL * PP
      RM = R * A / MU
      QM = .7 * P
c
  200 format('   Out of Table in StdAtm- too high !'//
     1        4x,'H =',f12.3,'  > 84.852 km'/)
c
      return
      end

