!+
PROGRAM SkinFrictionTables
! ---------------------------------------------------------------------------
! PURPOSE - A little test program to check that you have the function
!  TurbulentSkinFriction copied and compiled correctly.
! If you run it you should get the following output

!            0.0     1.0     2.0     3.0     5.0    10.0
!1.0E+05 0.00718 0.00676 0.00583 0.00489 0.00357 0.00235
!2.0E+05 0.00609 0.00571 0.00488 0.00403 0.00285 0.00170
!5.0E+05 0.00499 0.00466 0.00393 0.00321 0.00219 0.00119
!1.0E+06 0.00435 0.00405 0.00339 0.00274 0.00183 0.00094
!2.0E+06 0.00382 0.00355 0.00296 0.00237 0.00155 0.00077
!5.0E+06 0.00326 0.00302 0.00250 0.00198 0.00128 0.00060
!1.0E+07 0.00291 0.00269 0.00222 0.00175 0.00111 0.00051
!2.0E+07 0.00262 0.00241 0.00198 0.00155 0.00098 0.00043
!5.0E+07 0.00229 0.00211 0.00172 0.00134 0.00083 0.00036
!1.0E+08 0.00208 0.00192 0.00156 0.00121 0.00075 0.00031
!2.0E+08 0.00190 0.00175 0.00142 0.00110 0.00067 0.00028
!5.0E+08 0.00170 0.00156 0.00126 0.00097 0.00059 0.00024
!1.0E+09 0.00156 0.00143 0.00116 0.00089 0.00053 0.00021

!             50     100     150     200     250     300
!1.0E+05 0.00684 0.00681 0.00679 0.00678 0.00677 0.00676
!2.0E+05 0.00577 0.00575 0.00573 0.00572 0.00571 0.00571
!5.0E+05 0.00470 0.00469 0.00468 0.00467 0.00466 0.00466
!1.0E+06 0.00408 0.00407 0.00406 0.00405 0.00405 0.00405
!2.0E+06 0.00358 0.00356 0.00356 0.00355 0.00355 0.00355
!5.0E+06 0.00304 0.00303 0.00302 0.00302 0.00302 0.00302
!1.0E+07 0.00271 0.00270 0.00270 0.00269 0.00269 0.00269
!2.0E+07 0.00243 0.00242 0.00242 0.00242 0.00242 0.00241
!5.0E+07 0.00212 0.00212 0.00212 0.00211 0.00211 0.00211
!1.0E+08 0.00193 0.00192 0.00192 0.00192 0.00192 0.00192
!2.0E+08 0.00176 0.00176 0.00175 0.00175 0.00175 0.00175
!5.0E+08 0.00157 0.00156 0.00156 0.00156 0.00156 0.00156
!1.0E+09 0.00144 0.00144 0.00144 0.00144 0.00144 0.00143

!1.0E+05 0.00718
!2.0E+05 0.00609
!5.0E+05 0.00499
!1.0E+06 0.00435
!2.0E+06 0.00382
!5.0E+06 0.00326
!1.0E+07 0.00291
!2.0E+07 0.00262
!5.0E+07 0.00229
!1.0E+08 0.00208
!2.0E+08 0.00190
!5.0E+08 0.00170
!1.0E+09 0.00156
!Normal stop in SkinFrictionTables


  IMPLICIT NONE

  REAL,PARAMETER,DIMENSION(13):: RNTAB = (/ 1E5,2E5,5E5, 1E6,2E6,5E6, &
                                            1E7,2E7,5E7, 1E8,2E8,5E8, 1E9 /) 
  REAL,PARAMETER,DIMENSION(6):: MACHTAB = (/ 0,1,2,3,5,10 /)                                           
  
  INTEGER:: i,j
  REAL,DIMENSION(6):: cf
!----------------------------------------------------------------------------

  WRITE(*,'(8X,6F8.1)' ) MACHTAB
  DO j=1,SIZE(RNTAB)  ! all using default of T=250
    DO i=1,SIZE(MACHTAB)
      cf(i)=TurbulentSkinFriction(RNTAB(j), MACHTAB(i))
    END DO
    WRITE(*, '(ES8.1, 6F8.5)' ) RNTAB(j), cf
  END DO
  WRITE(*,*)

  WRITE(*,*) "             50     100     150     200     250     300"
  DO j=1,SIZE(RNTAB)
    DO i=1,SIZE(cf)
      cf(i)=TurbulentSkinFriction(RNTAB(j), 1.0, REAL(50*i))
    END DO
    WRITE(*, '(ES8.1, 6F8.5)' ) RNTAB(j), cf
  END DO
  WRITE(*,*)

  DO j=1,SIZE(RNTAB)  ! default values of Mach=0 and T=250
    cf(1)=TurbulentSkinFriction(RNTAB(j))
    WRITE(*, '(ES8.1, 6F8.5)' ) RNTAB(j), cf(1)
  END DO

  WRITE(*,*) "Normal stop in SkinFrictionTables"
  STOP

CONTAINS


!+
FUNCTION TurbulentSkinFriction(reynoldsNumber, machNumber, temperature) &
                               RESULT(cf) 
! ---------------------------------------------------------------------------
! PURPOSE - Compute the value of total flat-plate all-turbulent skin
!    friction using the method of Sommer and Short (NACA TM xxxx).
!    Although this semi-empirical method is somewhat more complicated
!    than other formulas that are given in textbooks of aerodynamics,
!    it does exhibit good correlation in the upper hypersonic regime
!    where temperature is a significant parameter.
!    At lower Mach Numbers, the result is quite insensitive to the
!    value of temperature.
  IMPLICIT NONE
  REAL,INTENT(IN):: reynoldsNumber   ! based on length of plate
  REAL,INTENT(IN),OPTIONAL:: machNumber
  REAL,INTENT(IN),OPTIONAL:: temperature          ! deg Kelvin

  INTEGER,PARAMETER:: DP = SELECTED_REAL_KIND(10,50)  
  REAL(DP):: cf
  REAL(DP),PARAMETER:: ONE = 1.0_DP
  REAL(DP),PARAMETER:: RF = 0.115_DP   ! recovery factor
  REAL(DP),PARAMETER:: SP1 = 0.088_DP   ! Sievells and Payne factor1
  REAL(DP),PARAMETER:: SP2 = 1.5_DP     ! Sievells and Payne factor2
  REAL(DP),PARAMETER:: SUTH = 110.4
  REAL(DP),PARAMETER:: TEMP_STRATOSPHERE = 300.0_DP  ! kelvins
  REAL(DP):: reynolds,mach,tempK,xx,cfi,z,denom,rstar
!----------------------------------------------------------------------------
  reynolds=reynoldsNumber   ! makes it DP
  IF (PRESENT(machNumber)) THEN
    mach=machNumber     ! makes it DP
  ELSE
    mach=0.0
  END IF
  IF (PRESENT(temperature)) THEN
    tempK=temperature   ! makes it DP
  ELSE
    tempK=TEMP_STRATOSPHERE
  END IF
  xx=LOG10(reynolds)-SP2
  cfi=SP1/(xx*xx)    ! Sievells and Payne
  z=ONE+RF*mach*mach
  denom=(tempK+SUTH)/(z*tempK+SUTH)
  rstar=reynolds/(denom*(z**2.5_DP))
  xx=xx/(LOG10(rstar)-SP2)
  cf=xx*xx*(cfi/z)
  RETURN
END Function TurbulentSkinFriction ! ----------------------------------------

END Program SkinFrictionTables   ! ==========================================