      PROGRAM FLPLPOST
      INCLUDE 'flplpost.inc'
C------------------------------------------------------------------------
C Program for initialization/post-processing flat plate case.
C------------------------------------------------------------------------
C These lines should be contained in "flplpost.inc"
c     PARAMETER(NX=120,NY=090)
c     REAL X(NX,NY),Y(NX,NY),Q(NX,NY,4),QT(NX,NY,3)
c     REAL RO(NX,NY),VX(NX,NY),VY(NX,NY),VT(NX,NY),Vmax(NX),
c    &     PS(NX,NY),TS(NX,NY),TT(NX,NY),EN(NX,NY),MACH(NX,NY),MU(NX,NY)
c     REAL AK(NX,NY),EPS(NX,NY),EPS2(NX,NY),OMG(NX,NY),OMGW(NX,NY),
c    &     TMU(NX,NY)
c     REAL YPLUS(NX,NY),UPLUS(NX,NY),TAUW(NX),VSHE(NX),Ufs(NX),
c    &     DEL(NX),THETA(NX),THETA1(NX),RETH(NX),REX(NX),CF(NX),
c    &     UV(NX,NY)
c     REAL Xref,Pref,Tref,AMACH,RE,Rref,Aref,Uref,AMU
c     REAL PTfs,TTfs,PSfs,TSfs,ROfs,VXfs,Machfs
c     REAL GAMMA,RGAS,TSUTH,PI,GC,C2B,C2BP,GAMI,GGM1
c     INTEGER I,J,Imax,Jmax,Ile,jBLedge(NX),NC1,TURBTYPE
c     INTEGER IU(0:2),M
C
c     COMMON /GRID/X,Y,Imax,Jmax,Ile
c     COMMON /SOLUTION/Q,QT,TURBTYPE
c     COMMON /PHYSICAL/RO,VX,VY,VT,Vmax,EN,PS,TS,TT,MACH,MU
c     COMMON /TURB/AK,EPS,EPS2,OMG,OMGW,TMU
c     COMMON /FPBL/YPLUS,UPLUS,TAUW,VSHE,Ufs,jBLedge,DEL,THETA,
c    &        THETA1,RETH,REX,CF,UV
c     COMMON /REF/Xref,Pref,Tref,AMACH,NC1,RE,Rref,Aref,Uref,AMU
c     COMMON /FS/PTfs,TTfs,PSfs,TSfs,ROfs,VXfs,Machfs
c     COMMON /CONST/GAMMA,RGAS,PI,GC,TSUTH,C2B,C2BP,GAMI,GGM1
c     COMMON /PRNT/IU
C------------------------------------------------------------------------
C Units used:
C     P    lbf/ft^2     E   lbf/ft^2              AK   k/(aref^2)
C     T    Rankine      mu  lbm/(ft sec)          EPS  eps/(aref^3/xref)
C     rho  lbm/ft^3     R   lbf ft/(lbm 'R)       OMG  omg/(aref/xref)
C     x    ft           Gc  lbm ft/(lbf sec^2)    TMU  tmu/mu
C     V    ft/sec       TAU lbf/ft^2              UV   (dim-less)
C Relations used:
C     P  = rho*R*T      a = sqrt(GAMMA*Rgas*T*Gc)
C     Re = rho*V*x/mu   E = Ps/GAMMA + 0.5*rho*V^2/Gc
C------------------------------------------------------------------------
C Files Used:
C      7  PARC Restart file                                    <input>
C      8  PARC Restart file                                    <output>
C      9  Plot3d Grid file                                     <input>
C     10  History file                                         <output>
C     11  Tecplot (DIM-LESS) Flow Values (X,Y,Rho*,P*,T*,U*,M) <output>
C           (DIM-IZED by physical units, not PARC ref values)
C     12  Tecplot Turbulence (DIM-LESS) Values                 <output>
C           (X,Y,K*,Eps*,Omg*,MUt*,y+,U+)
C           (DIM-IZED by physical units, not PARC ref values)
C     13  Plot3d  Grid file                                    <output>
C     14  Tecplot Grid file, columnar form                     <output>
C     15  Tecplot FPBL Values (X,ReX,ReTH,TH,TH1,Cf,y+)        <output>
C--Use These to Monitor Convergence (iter,I,ReX,ReTH,Rho,P,T,K,Eps,MUt)
C     16  Maximum Values at ReX= 1,000,000                     <output>
C     17  Maximum Values at ReX= 4,000,000                     <output>
C     18  Maximum Values at ReX=10,000,000                     <output>
C------------------------------------------------------------------------
C Define Reference Values
      Tref  =  530.0
      Pref  = 2116.8
      Xref  =    1.0
      AMACH =    0.2
      GAMMA =    1.4 !NOT..........Read in from RST file
      RGAS  =   53.375
      TSUTH =  198.6
      PI    =    3.14159265359
      GC    =   32.174
      Ile   =    15
C
C Begin Calculations
      CALL ReadRST                ! Read NPARC restart file
      CALL CompRef                ! Compute reference values
      CALL CompPhys               ! Convert non-dim q to dim P,T,...
      CALL CompFS                 ! Compute freestream properties
      CALL CompTurb               ! Compute turbulence info from q
      CALL CompFPBL               ! Compute flat plate boundary layer
                                  !  properties (displ,mom.thks,etc)
c     CALL WritePhys              ! Write dim props to file
      CALL WriteFPBL              ! Write flat plate boundary layer
                                  !  info to file
c     CALL WriteTurb              ! Write turbulence values to file
c     CALL WriteTurb2             ! Write turbulence profiles to file
c     i = iFindReX( 1000000.)     ! Find i-index where Rex=1000000.
c     CALL WriteMax(16,i)         ! Write max values (P,T,...) at i-index
c     CALL WriteProf(20,i)        ! Write profiles / u
c     CALL WriteProf2(30,i)       ! Write profiles / utau
c     i = iFindReX( 4000000.)
c     CALL WriteMax(17,i)
c     CALL WriteProf(21,i)
c     CALL WriteProf2(31,i)
c     i = iFindReX(10000000.)
c     CALL WriteMax(18,i)
c     CALL WriteProf(22,i)
c     CALL WriteProf2(32,i)
c
      i = iFindReX(0.191E06)
      CALL WriteProf2(35,i)
      i = iFindReX(0.106E07)
      CALL WriteProf2(35,i)
      i = iFindReX(0.270E07)
      CALL WriteProf2(35,i)
      i = iFindReX(0.500E07)
      CALL WriteProf2(35,i)
      i = iFindReX(0.762E07)
      CALL WriteProf2(35,i)
      i = iFindReX(0.103E08)
      CALL WriteProf2(35,i)
c
c     i = iFindReX(0.406E06)
c     CALL WriteProf2(36,i)
c     i = iFindReX(0.627E06)
c     CALL WriteProf2(36,i)
c
C---Set Up New Grid/Solution
C     CALL ReadGrid               ! Read plot3d grid file
C     CALL NewYGrid               ! Redistribute y-grid
C     CALL WriteGrid              ! Write plot3d grid file
C     CALL InitSol                ! Initialize solution q
C     CALL WriteRST               ! write NPARC restart file
C
      WRITE(*,*) 'History File Written to fort.10'
      WRITE(*,*)
      END
C========================================================================
      SUBROUTINE ReadRST
      INCLUDE "flplpost.inc"
C----------------------------------------------------
C     Read PARC2D restart file.
C----------------------------------------------------
      WRITE(*,*) 'Is this a BL, k-e, or k-w file (0/1/2)?  '
      READ(*,*) TURBTYPE
c-BUG READ(7) NC1,GAMMA
      READ(7) NC1
      READ(7) Imax,Jmax
      READ(7) ((X(I,J),    I=1,Imax),J=1,Jmax),
     &        ((Y(I,J),    I=1,Imax),J=1,Jmax)
      READ(7) (((Q(I,J,N), I=1,Imax),J=1,Jmax),N=1,4)
      IF (TURBTYPE .eq. 0) THEN
        DO 200 I=1,Imax
          DO 100 J=1,Jmax
            QT(I,J,1) = 0.0
            QT(I,J,2) = 0.0
            QT(I,J,3) = 0.0
100       CONTINUE
200     CONTINUE
c-out   READ(7) (((QT(I,J,N),I=1,Imax),J=1,Jmax),N=3,3)
      ELSE
        READ(7) (((QT(I,J,N),I=1,Imax),J=1,Jmax),N=1,3)
      ENDIF
C
      DO 1000 M=1,IU(0)
        WRITE(IU(M),909) Imax,Jmax,NC1,GAMMA
        WRITE(IU(M),*)
909     FORMAT(I4,' x',I4,' Restart file read in @ iteration ',I8,
     &         ' with GAMMA=',F4.2)
1000  CONTINUE
      RETURN
      END
C========================================================================
      SUBROUTINE WriteRST
      INCLUDE "flplpost.inc"
C----------------------------------------------------
C     Read PARC2D restart file.
C----------------------------------------------------
      WRITE(8) NC1,GAMMA
      WRITE(8) Imax,Jmax
      WRITE(8) ((X(I,J), I=1,Imax),J=1,Jmax),
     &         ((Y(I,J), I=1,Imax),J=1,Jmax)
      WRITE(8) (((Q(I,J,N), I=1,Imax),J=1,Jmax),N=1,4)
      WRITE(8) (((QT(I,J,N),I=1,Imax),J=1,Jmax),N=1,3)
C
      DO 1000 M=1,IU(0)
        WRITE(IU(M),909) Imax,Jmax,NC1
        WRITE(IU(M),*)
909     FORMAT(I4,' x',I4,' Restart file written @ iteration ',I8)
1000  CONTINUE
      RETURN
      END
C========================================================================
      SUBROUTINE ReadGrid
      INCLUDE "flplpost.inc"
C----------------------------------------------------
C     Read Plot3d grid file.
C----------------------------------------------------
      READ(9) NB
      READ(9) Imax,Jmax
      READ(9) ((X(I,J),    I=1,Imax),J=1,Jmax),
     &        ((Y(I,J),    I=1,Imax),J=1,Jmax)
C
      DO 1000 M=1,IU(0)
        WRITE(IU(M),909) Imax,Jmax
        WRITE(IU(M),*)
909     FORMAT(I4,' x',I4,' Grid file read.  Grid units in feet.')
1000  CONTINUE
      RETURN
      END
C========================================================================
      SUBROUTINE WriteGrid
      INCLUDE "flplpost.inc"
C-----------------------------------------------
C     Write grid (x,y) to files for viewing.
C-----------------------------------------------
C Write grid in Tecplot Format
      WRITE(14,*) 'TITLE = "2-D Flat Plate Grid"'
      WRITE(14,*) 'VARIABLES = "X","Y"'
      WRITE(14,901) Imax,Jmax
901   FORMAT('ZONE T = "grid", I=',I3,',J=',I3,', F=POINT')
      DO 200 J=1,Jmax
        DO 100 I=1,Imax
          WRITE(14,919) X(I,J),Y(I,J)
919       FORMAT(E12.6,E14.6)
100     CONTINUE
200   CONTINUE
C
C Write grid in Plot3d Format
      WRITE(13) 1
      WRITE(13) Imax,Jmax
      WRITE(13) ((X(I,J), I=1,Imax),J=1,Jmax),
     &          ((Y(I,J), I=1,Imax),J=1,Jmax)
C Update User
      DO 1000 M=1,IU(0)
        WRITE(IU(M),909) Imax,Jmax,' Tecplot Grid written to fort.9'
        WRITE(IU(M),909) Imax,Jmax,' Plot3d  Grid written to fort.13'
        WRITE(IU(M),*)
909     FORMAT(I4,' x',I4,A)
1000  CONTINUE
C     
      RETURN
      END
C========================================================================
      SUBROUTINE WritePhys
      INCLUDE "flplpost.inc"
C
C Write (NORMALIZED) Values in Tecplot Format
      WRITE(11,899) '##################################################'
      WRITE(11,899) '#  Non-Dimensional Physical Flow Variables:      #'
      WRITE(11,899) '#    Rho* = RHO/RHOfs    T* = T/Tfs              #'
      WRITE(11,899) '#    P*   = P/Pfs        U* = U/Ufs              #'
      WRITE(11,899) '##################################################'
      WRITE(11,899) 'TITLE = "2-D Flat Plate Grid"'
      WRITE(11,899) 
     &  'VARIABLES = "X","Y","Rho*","P*","T*","U*","y+","U+","Mach"'
899   FORMAT(A)
      WRITE(11,903) Imax,Jmax
903   FORMAT('ZONE T = "grid", I=',I3,',J=',I3,', F=POINT')
      DO 200 J=1,Jmax
        DO 100 I=1,Imax
          WRITE(11,919) X(I,J),Y(I,J),RO(I,J)/ROfs,PS(I,J)/PSfs,
     &      TS(I,J)/TSfs,VX(I,J)/VXfs,Yplus(I,J),Uplus(I,J),MACH(I,J)
919       FORMAT(E12.6,8E14.6)
100     CONTINUE
200   CONTINUE
      DO 500 M=1,IU(0)
        WRITE(IU(M),*) 'NORMALIZED Flow Values Written to fort.11'
        WRITE(IU(M),*)
500   CONTINUE
      RETURN
      END
C========================================================================
      SUBROUTINE WriteFPBL
      INCLUDE "flplpost.inc"
C
C Write (NORMALIZED) Flat Plate Boundary Layer Values in Tecplot Format
      WRITE(15,899) '##################################################'
      WRITE(15,899) '#  Flat Plate Boundary Layer Values:             #'
      WRITE(15,899) '##################################################'
      WRITE(15,899) 'TITLE = "2-D Flat Plate BL"'
      WRITE(15,898) 
898   FORMAT('VARIABLES = "X","ReX","ReTH","Theta","Del   ","Cf","y+",',
     &  '"jBLedge","del"')
899   FORMAT(A)
      WRITE(15,903) Imax-Ile+1
903   FORMAT('ZONE T = "grid", I=',I3,', F=POINT')
      DO 100 I=Ile,Imax
        WRITE(15,919) X(I,1),ReX(I),ReTH(I),THETA(I),DEL(I),Cf(I),
     &    Yplus(I,2),jBLedge(I),DEL(I)
919     FORMAT(E12.6,6E14.6,I4,E14.6)
100   CONTINUE
      DO 500 M=1,IU(0)
        WRITE(IU(M),*) 'FPBL Values Written to fort.15'
        WRITE(IU(M),*)
500   CONTINUE
      RETURN
      END
C========================================================================
      SUBROUTINE WriteTurb
      INCLUDE "flplpost.inc"
      REAL ZERO
      ZERO = 1.0E-30
C
C Write (DIM-AL) Turbulent Values in Tecplot Format
      WRITE(12,899) '##################################################'
      WRITE(12,899) '#  Non-Dimensional turbulence values:            #'
      WRITE(12,899) '#    K*   = K/U2     eps* = eps/(U3/d)   U=Ufs   #'
      WRITE(12,899) '#    MUt* = MUt/MU   omg* = omg/(U/d)    d=del   #'
      WRITE(12,899) '##################################################'
      WRITE(12,899) 'TITLE = "2-D Flat Plate Grid"'
      WRITE(12,899) 
     &  'VARIABLES = "X","Y","K*","Eps*","Omg*","MUt*","y+","U+"',
     &  ',"k+","eps+"'
899   FORMAT(A)
      WRITE(12,901) Imax-Ile+1,Jmax
901   FORMAT('ZONE T = "grid", I=',I3,',J=',I3,', F=POINT')
      DO 200 J=1,Jmax
        DO 100 I=Ile,Imax
          WRITE(12,919) X(I,J),Y(I,J),
     &      QT(I,J,1)*(aref/Ufs(I))**2,
     &      EPS(I,J)*((aref**3/Xref)/(Ufs(I)**3/(del(I)+ZERO))),
     &      OMG(I,J)*((aref/xref)/(Ufs(I)/(del(I)+ZERO))),
     &      QT(I,J,3),Yplus(I,J),Uplus(I,J),
     &      QT(I,J,1)*(aref/VSHE(I))**2,                        ! DAY
     &      EPS(I,J)*aref**3/xref*MU(I,J)/RO(I,J)/VSHE(I)**4    ! DAY
919       FORMAT(E12.6,9E14.6)
100     CONTINUE
200   CONTINUE
      DO 500 M=1,IU(0)
        WRITE(IU(M),*) 'NON-DIM-AL Turbulence Values Written to fort.12'
        WRITE(IU(M),*)
500   CONTINUE
      RETURN
      END
C========================================================================
      SUBROUTINE WriteTurb2
      INCLUDE "flplpost.inc"
      INTEGER INDEX(0:4)
      REAL ZERO,MMperFT
c-DAY DATA INDEX/4,84,90,99,107/     ! Fraser
      DATA INDEX/4,38,52,69, 95/     ! Flat plate
      ZERO = 1.0E-30
      MMperFT=12.*25.4
C
C Write (DIM-AL) Turbulent Values in Tecplot Format
C Designed for Fraser Diffuser Output
      DO 200 IN=1,INDEX(0)
      I=INDEX(IN)
      IO=49+IN
      WRITE(IO,899) '##################################################'
      WRITE(IO,899) '#  Non-Dimensional turbulence values:            #'
      WRITE(IO,899) '#    K*   = K/U2     eps* = eps/(U3/d)   U=Ucl   #'
      WRITE(IO,899) '#    MUt* = MUt/MU   omg* = omg/(U/d)    d=dia(x)#'
      WRITE(IO,899) '#    PROD*= PROD/MU/(U/d)**2                     #'
      WRITE(IO,899) '#  Grid Index for this profile                   #'
      WRITE(IO,897) INDEX(IN)
      WRITE(IO,898) X(INDEX(IN),1)*MMperFT
      WRITE(IO,899) '#  Should be close to                            #'
      WRITE(IO,899) '#    X=116, 236, 438, 642 mm                     #'
      WRITE(IO,899) '##################################################'
      WRITE(IO,899) 'TITLE = "Fraser Diffuser Flow"'
      WRITE(IO,899) 'VARIABLES = "y","y+","K*","Eps*","Omg*",',
     &  '"MUt*","PROD*","EPROD*","PRODE*","UV*"'
896   FORMAT('#    X(',I3,')=',F6.1,' mm                            #')
897   FORMAT('#    I=',10(I3,','),'  #')
898   FORMAT('#    X=',10(F5.0,','),' mm  #')
899   FORMAT(A)
      WRITE(IO,901) 
901   FORMAT('ZONE T = "grid", F=POINT')
        DO 100 J=1,Jmax
          PROD=0.
          IF ((1.lt.J).and.(J.lt.Jmax)) THEN
            DUDX=(VX(I+1,J)-VX(I-1,J))/(X(I+1,J)-X(I-1,J))
            DUDY=(VX(I,J+1)-VX(I,J-1))/(Y(I,J+1)-Y(I,J-1))
            DVDX=(VY(I+1,J)-VY(I-1,J))/(X(I+1,J)-X(I-1,J))
            DVDY=(VY(I,J+1)-VY(I,J-1))/(Y(I,J+1)-Y(I,J-1))
            PROD=2.*QT(I,J,3)*(DUDX*DUDX+DVDY*DVDY)+
     &        QT(I,J,3)*(DUDY+DVDX)**2.
c-DAY       PROD=PROD/(Vmax(I)/Y(I,Jmax))**2.    ! Fraser
            TAU=TMU(I,J)*MU(I,J)*DUDY
          ENDIF
          xTAUW= -MU(I,Jmax)*((VX(I,Jmax)-VX(I,Jmax-1))
     &      /(Y(I,Jmax)-Y(I,Jmax-1)+ZERO))/GC
          xVSHE= (xTAUW/RO(I,Jmax-1)*GC)**0.5
c-DAY     xYPLUS=(Y(I,Jmax)-Y(I,J))*xVSHE*RO(I,J)/MU(I,J)   ! Freaser
          xYPLUS=(          Y(I,J))*VSHE(I)*RO(I,J)/MU(I,J)   ! FLPL
          xUV=TAU/(RO(I,J)*Vmax(I)**2)
          EP=EPS(I,J)*(aref**3/Xref)
c-DAY&       /(Vmax(I)**3/(Y(I,Jmax)+ZERO))
          WRITE(IO,919) (Y(I,Jmax)-Y(I,J))*MMperFT,
     &      xYPLUS,
     &      QT(I,J,1)*(aref/Vmax(I))**2,
     &      EPS(I,J)*((aref**3/Xref)/(Vmax(I)**3/(Y(I,Jmax)+ZERO))),
     &      OMG(I,J)*((aref/xref)/(Vmax(I)/(Y(I,Jmax)+ZERO))),
     &      QT(I,J,3),
     &      PROD,
     &      EP/(PROD+ZERO),
     &      PROD/(EP+ZERO) *RO(I,J),
     &      -1.*xUV
919       FORMAT(E12.6,9E14.6)
100     CONTINUE
        DO 150 M=1,IU(0)
        WRITE(IU(M),*)
     &    'NON-DIM-AL Turbulence Values Written to fort.',io
        WRITE(IU(M),*)
150     CONTINUE
200   CONTINUE
      RETURN
      END
C========================================================================
      SUBROUTINE WriteProf(UN,I)
      INCLUDE "flplpost.inc"
      INTEGER UN
      REAL ZERO
      ZERO = 1.0E-30
C
C Write (DIM-AL) Turbulent Values in Tecplot Format
      WRITE(UN,899) '##################################################'
      WRITE(UN,899) '#  Non-Dimensional turbulence profiles:          #'
      WRITE(UN,899) '#    K*   = K/U2     eps* = eps/(U3/d)   U=Ufs   #'
      WRITE(UN,899) '#    MUt* = MUt/MU   omg* = omg/(U/d)    d=del   #'
      WRITE(UN,899) '#    -uv* = -uv/U**2                             #'
      WRITE(UN,899) '##################################################'
      WRITE(UN,899) 'TITLE = "2-D Flat Plate Grid"'
      WRITE(UN,898) 
898   FORMAT('VARIABLES= "X","Y_d","K*","Eps*","Omg*","MUt*",',
     &'"y+","U+","-uv*"')
899   FORMAT(A)
      WRITE(UN,901) Jmax
901   FORMAT('ZONE T = "grid", I=',I3,', F=POINT')
      WRITE(UN,899) 
     &'#X   Y_del   K*   Eps*   Omg*   MUt*   y+   U+   -uv*'
      DO 200 J=1,Jmax
          WRITE(UN,919) X(I,J),Y(I,J)/del(I),
     &      QT(I,J,1)*(aref/Ufs(I))**2,
     &      EPS(I,J)*((aref**3/Xref)/(Ufs(I)**3/(del(I)+ZERO))),
     &      OMG(I,J)*((aref/xref)/(Ufs(I)/(del(I)+ZERO))),
     &      QT(I,J,3),Yplus(I,J),Uplus(I,J),
     &      UV(I,J)
919       FORMAT(E12.6,8E14.6)
200   CONTINUE
      DO 500 M=1,IU(0)
        WRITE(IU(M),*) 
     &    'NON-DIM-AL Turbulence Profiles Written to fort.',UN
        WRITE(IU(M),*)
500   CONTINUE
      RETURN
      END
C========================================================================
      SUBROUTINE WriteProf2(UN,I)
      INCLUDE "flplpost.inc"
      INTEGER UN
      REAL ZERO
      ZERO = 1.0E-10
C
C Write (Normalized) Velocity & Turbulent Values in Tecplot Format
      WRITE(UN,899) '##################################################'
      WRITE(UN,899) '#  Non-Dimensional turbulence profiles:          #'
      WRITE(UN,899) '#    y+   = y*Utau*ro/mu                         #'
      WRITE(UN,899) '#    U+   = U/Utau                               #'
      WRITE(UN,899) '#    k+   = k/Utau^2                             #'
      WRITE(UN,899) '#    eps+ = eps*mu/ro/Utau^4                     #'
      WRITE(UN,899) '#    omg+ = omg/(U/d)                            #'
      WRITE(UN,899) '#    -uv+ = -uv/Utau^2                           #'
      WRITE(UN,899) '#    P_dis= Production/Dissipation (dim-less)    #'
      WRITE(UN,899) '#  Note:  epsilon is the true dissipation rate.  #'
      WRITE(UN,899) '#         eps = chien_eps + 2*mu*k/ro/y^2        #'
      WRITE(UN,899) '##################################################'
      WRITE(UN,899) 'TITLE = "2-D Flat Plate Grid"'
      WRITE(UN,898) 
898   FORMAT('VARIABLES= "y+","u+","K+","Eps+","-uv+","P_dis","D"')
899   FORMAT(A)
      WRITE(UN,901) Jmax
901   FORMAT('ZONE T = "grid", I=',I3,', F=POINT')
c-out WRITE(UN,899) 
c-out&'#X   Y_del   K*   Eps*   Omg*   MUt*   y+   U+   -uv*'
      DO 200 J=1,Jmax
        PROD=0.
        DISS=1.
        IF ((1.lt.J).and.(J.lt.Jmax)) THEN
          DUDX=(VX(I+1,J)-VX(I-1,J))/(X(I+1,J)-X(I-1,J))
          DUDY=(VX(I,J+1)-VX(I,J-1))/(Y(I,J+1)-Y(I,J-1))
          DVDX=(VY(I+1,J)-VY(I-1,J))/(X(I+1,J)-X(I-1,J))
          DVDY=(VY(I,J+1)-VY(I,J-1))/(Y(I,J+1)-Y(I,J-1))
          PROD=2.*QT(I,J,3)*(DUDX*DUDX+DVDY*DVDY)+
     &            QT(I,J,3)*(DUDY+DVDX)**2.
          PROD=PROD/(aref/xref)**2.         ! prod/(muref*aref^2/xref^2)
          DISS=Q(I,J,1)*EPS(I,J)*RE         ! diss/(muref*aref^2/xref^2)
     &        +2.*MU(I,J)/AMU/Q(I,J,1)*QT(I,J,1)
     &           *(xref/(Y(I,J)-Y(I,1)+ZERO))**2
        ENDIF
        PRODDISS=PROD/DISS
        WRITE(UN,919)
     &    Y(I,J)*VSHE(I)*RO(I,J)/MU(I,J),
     &    VX(I,J)/VSHE(I),
     &    QT(I,J,1)*(aref/VSHE(I))**2,
     &    EPS(I,J)*(aref**3/Xref)*MU(I,J)/RO(I,J)/VSHE(I)**4
     &      + 2.*QT(I,J,1)*(aref*MU(I,J)/RO(I,J)/(Y(I,J)-Y(I,1)+ZERO)
     &                     /Xref)**2/VSHE(I)**4,
     &    UV(I,J)*(Ufs(I)/VSHE(I))**2,
     &    PRODDISS,
     &        2.*QT(I,J,1)*(aref*MU(I,J)/RO(I,J)/(Y(I,J)-Y(I,1)+ZERO)
     &                     /Xref)**2/VSHE(I)**4
919     FORMAT(E12.6,6E14.6)
200   CONTINUE
      DO 500 M=1,IU(0)
        WRITE(IU(M),*) 
     &    'NORMALIZED Turbulence Profiles Written to fort.',UN
        WRITE(IU(M),*)
500   CONTINUE
      RETURN
      END
C========================================================================
      SUBROUTINE WriteMax(IU1,I1)
      INCLUDE "flplpost.inc"
      INTEGER  I1,IU1
      REAL ROmax,Pmax,Tmax,AKmax,EPSmax,TMUmax
      DATA ROmax,Pmax,Tmax,AKmax,EPSmax,TMUmax/0.0,0.0,0.0,0.0,0.0,0.0/
C------------------------------------------------------------------
C     Find the maximum values at the x-location specified by I1
C------------------------------------------------------------------
      IF (I1 .gt. 0) THEN
        DO 500 J=1,Jmax
          ROmax  = MAX(ROmax, RO(I1,J))
          Pmax   = MAX(Pmax,  PS(I1,J))
          Tmax   = MAX(Tmax,  TS(I1,J))
          AKmax  = MAX(AKmax, AK(I1,J))
          EPSmax = MAX(EPSmax,EPS(I1,J))
          TMUmax = MAX(TMUmax,TMU(I1,J))
500     CONTINUE
        WRITE(IU1,598) '# MAXIMUM VALUES'
        WRITE(IU1,598) 
     &'# iter   I     ReX     ReTH    RO    P     T    K     Eps    MUt'
598     FORMAT(A)
        WRITE(IU1,599) INT(NC1),I1,ReX(I1),ReTH(I1),ROmax,Pmax,Tmax,
     &    AKmax,EPSmax,TMUmax
599     FORMAT(I6,I5,2E10.2,F8.5,F8.1,F8.2,F10.7,F10.7,F9.2)
      ENDIF
      RETURN
      END
C========================================================================
      BlOCK DATA
      COMMON /PRNT/IU(0:2)
C------------------------------------------
C     Initialize output display files.  
C     First number tell how many files.
C------------------------------------------
      DATA IU/2,6,10/
      END
C========================================================================
      SUBROUTINE CompRef
      INCLUDE 'flplpost.inc'
C
C Compute Other Reference Values
      GAMI  = GAMMA - 1.0
      GGM1  = GAMMA*(GAMMA - 1.0)
      C2B   = TSUTH/TREF
      C2BP  = C2B + 1.0
      RREF  = PREF/(RGAS*TREF)
      AREF  = SQRT(TREF*GC*RGAS*GAMMA)
C     AMU   = 2.27*(TREF**1.5/(TREF+TSUTH))*1.0E-08   (lbf sec/ft^2)
      AMU   = 7.293E-07*(TREF**1.5/(TREF+TSUTH))
      RE    = RREF*AREF*XREF/AMU
C     UREF  = AREF
C
      DO 1000 M=1,IU(0)
        WRITE(IU(M),*) 'Reference values are:'
        WRITE(IU(M),909) '  Pref  = ',Pref  ,' lbf/ft^2'
        WRITE(IU(M),909) '  Tref  = ',Tref  ,' Rankine'
        WRITE(IU(M),909) '  Xref  = ',Xref  ,' ft'
        WRITE(IU(M),909) '  AMACH = ',AMACH ,' '
        WRITE(IU(M),909) '  ROref = ',Rref  ,' lbm/ft^3'
        WRITE(IU(M),909) '  Aref  = ',Aref  ,' ft/sec'
        WRITE(IU(M),919) '  Re    = ',Re
        WRITE(IU(M),*) 
909     FORMAT(A,F13.6,A)
919     FORMAT(A,E13.6)
1000  CONTINUE
      RETURN
      END
C========================================================================
      SUBROUTINE CompPhys
      INCLUDE "flplpost.inc"
C
C PERFORM POST-PROCESSING CALCULATIONS (COMPUTE DIM-AL VALUES)
      DO 500 I = 1,Imax
        DO 300 J = 1,Jmax
          X(I,J)    = X(I,J)*Xref
          Y(I,J)    = Y(I,J)*Xref
          RO(I,J)   = Q(I,J,1)*Rref
          VX(I,J)   = Q(I,J,2)*Aref/Q(I,J,1)
          VY(I,J)   = Q(I,J,3)*Aref/Q(I,J,1)
          VT(I,J)   = ((ABS(VX(I,J)))**2. + (ABS(VY(I,J)))**2.)** 0.5
          EN(I,J)   = Q(I,J,4)*GAMMA*Pref
          PS(I,J)   = (EN(I,J) - 0.5*RO(I,J)*VT(I,J)**2./GC)*GAMI
          TS(I,J)   = PS(I,J)/(RO(I,J)*RGAS)
          MACH(I,J) = VT(I,J)/(GAMMA*RGAS*GC*TS(I,J))**(0.5)
          TT(I,J)   = TS(I,J)*(1. + (MACH(I,J)**2.)*GAMI/2.)
          MU(I,J)   = .0000007293105*(TS(I,J)**1.5)/(198.72+TS(I,J)) 
          Vmax(I) = MAX(Vmax(I),VX(I,J))
300     CONTINUE
500   CONTINUE
C Update User
      DO 1000 M=1,IU(0)
      WRITE(IU(M),*) 'Dimensional flow properties have been calculated.'
        WRITE(IU(M),*)
1000  CONTINUE
      RETURN
      END
C========================================================================
      SUBROUTINE CompFS
      INCLUDE 'flplpost.inc'
      REAL Qfs(1:4)
      INTEGER Jfs,Jfsdef
C
      PTfs = 0.0
      TTfs = 0.0
      PSfs = 0.0
      TSfs = 0.0
      ROfs = 0.0
      VXfs = 0.0
      Machfs = 0.0
      DO 5 N=1,4
        Qfs(N) = 0.0
5     CONTINUE
      WRITE(*,*) 'Computing Freestream Values,'
10    Jfsdef = Jmax - 5
      CALL IUSERINPUT('  Enter vertical gridline: ',Jfs,Jfsdef)
      IF (Jfs .GT. Jmax) GOTO 10
C
      DO 500 I=1,Imax
        PTfs = PTfs + 
     &         (PS(I,Jfs)*(1.0 + 0.5*GAMI*Mach(I,Jfs)**2)**GGM1)/Imax
        TTfs = TTfs + TT(I,Jfs)/Imax
        PSfs = PSfs + PS(I,Jfs)/Imax
        TSfs = TSfs + TS(I,Jfs)/Imax
        ROfs = ROfs + RO(I,Jfs)/Imax
        VXfs = VXfs  + VX(I,Jfs)/Imax
        Machfs = Machfs + Mach(I,Jfs)/Imax
        DO 400 N=1,4
          Qfs(N) = Qfs(N) + Q(I,Jfs,N)/Imax
400     CONTINUE
500   CONTINUE
C
      DO 700 M=1,IU(0)
        WRITE(IU(M),609) Jfs
        WRITE(IU(M),619) 'Pt   = ',PTfs,' lbf/ft^2','Q1 = ',Qfs(1)
        WRITE(IU(M),619) 'Tt   = ',TTfs,' Rankine ','Q2 = ',Qfs(2)
        WRITE(IU(M),619) 'Ps   = ',PSfs,' lbf/ft^2','Q3 = ',Qfs(3)
        WRITE(IU(M),619) 'Ts   = ',TSfs,' Rankine ','Q4 = ',Qfs(4)
        WRITE(IU(M),629) 'Rho  = ',ROfs,' lbm/ft^3'
        WRITE(IU(M),639) 'Mach = ',Machfs
        WRITE(IU(M),*)
609     FORMAT(' Computed Freestream Values along J=',I3)
619     FORMAT(3x,A,F11.5,A,12x,A,E13.6)
629     FORMAT(3x,A,F11.5,A)
639     FORMAT(3x,A,F11.5)
700   CONTINUE
      RETURN
      END
C========================================================================
      SUBROUTINE CompTurb
C
      INCLUDE 'flplpost.inc'
      REAL ZERO
      ZERO = 1.0E-25
C
      WRITE(10,*) 'Computing Turbulent Properties.'
      DO 100 I=1,Imax
        DO 50 J=1,Jmax
          AK(I,J) = QT(I,J,1)
          TMU(I,J) = QT(I,J,3)
50      CONTINUE
100   CONTINUE
      IF (TURBTYPE .eq. 0) THEN
        WRITE(10,*) '  Computing Baldwin-Lomax Values.'
        DO 110 I=1,Imax
          DO 120 J=1,Jmax
            EPS(I,J) = 0.0
            EPS2(I,J)= 0.0
            OMG(I,J) = 0.0
            OMGW(I,J) = 0.0
120       CONTINUE
110     CONTINUE
      ELSEIF (TURBTYPE .eq. 1) THEN
        WRITE(10,*) '  Computing OMEGA from EPSILON.'
        rt   = 0.00
        beta = 0.09
        DO 200 I=1,Imax
          DO 150 J=1,Jmax
            EPS(I,J) = QT(I,J,2)
            EPS2(I,J) = QT(I,J,2)
            OMG(I,J) = EPS(I,J)/(AK(I,J)+ZERO)/beta
            OMGW(I,J)= 6.*MU(I,2)/RO(I,2)/(beta*Y(I,J)**2+ZERO) /
     &                 (aref/xref)
150       CONTINUE
200     CONTINUE
      ELSE
        WRITE(10,*) '  Computing EPSILON from Omega.'
        DO 500 I=1,Imax
          DO 450 J=1,Jmax
            OMG(I,J) = QT(I,J,2)
            U = Q(I,J,2)/Q(I,J,1)
            V = Q(I,J,3)/Q(I,J,1)
            UU= U*U + V*V
            ei  = GGM1*(Q(I,J,4)/Q(I,J,1) - UU/2.0)
            fmu = (C2BP*ei**1.5)/(C2B+ei)
            rt  = Q(I,J,1)*AK(I,J)/fmu/(OMG(I,J)+ZERO)*re
            rtcrs4 = (rt/8.)**4
            beta = 0.09*(5./18.+rtcrs4)/(1.+rtcrs4)
            EPS(I,J) = OMG(I,J)*AK(I,J)*beta
            EPS2(I,J)= omg(i,j)*ak(i,j)*0.09
            OMGW(I,J)= 6.*MU(I,2)/RO(I,2)/(beta*Y(I,J)**2+ZERO) /
     &                 (aref/xref)
450       CONTINUE
500     CONTINUE
      ENDIF
      RETURN
      END
C========================================================================
      SUBROUTINE CompFPBL
C
      INCLUDE "flplpost.inc"
      INTEGER jmom,jmomo,jstop,istop,nptsFS
C
      jmomo = 20
      DO 1000 I=Ile,Imax
C-------initialize flags
        nptsFS    = 0
        jBLedge(I)= 0
        Ufs(I)    = 0.0
        THETA(I)  = 0.0
        THETA1(I)  = 0.0
        jmom      = 0
        jstop     = 0
        istop     = 0
        DO 500 J=2,Jmax
C---------compute momentum thickness
          IF (((VX(I,J).GT.VX(I,J-1)) .AND. (jstop.LT.2)) .AND.
     &      ((jmom-jmomo).LT.2)) THEN
            DELTHE = VX(I,J)/Vmax(I)*(1.-(VX(I,J)/Vmax(I)))*
     &        (Y(I,J+1)-Y(I,J-1))/2.
            THETA(I) = THETA(I) + DELTHE
              jmom = J
          ELSE
            jmomo = jmom
            jstop = 2
          ENDIF
C---------compute boundary layer edge using linear interpolation
          IF (((VX(I,J)/Vmax(I)).GT.(0.99)) .AND. (istop.LT.2)) THEN
            DEL(I) = Y(I,J-1) +(Y(I,J)-Y(I,J-1))*
     &               (0.99*Vmax(I)-VX(I,J-1))/(VX(I,J)-VX(I,J-1))
            jBLedge(I) = J
            istop = 2
          ENDIF
500     CONTINUE
        DO 550 J=2,jBLedge(I)
          DELTHE = VX(I,J)/Vmax(I)*(1.-(VX(I,J)/Vmax(I)))*
     &      (Y(I,J+1)-Y(I,J-1))/2.
          THETA1(I) = THETA1(I) + DELTHE
550     CONTINUE
C---------compute freesteam velocity
        DO 600 J=2,Jmax
          IF ((J .GT. jBLedge(I)+10) .OR. (J .GT. 70)) THEN
            Ufs(I) = Ufs(I) + VX(I,J)
            nptsFS = nptsFS + 1
          ENDIF
600     CONTINUE
        Ufs(I) = Ufs(I)/nptsFS
C-------compute boundary layer properties
        TAUW(I)    = MU(I,1)*(VX(I,2)/Y(I,2))/GC
        VSHE(I)    = (TAUW(I)/RO(I,2)*GC)**0.5
        RETH(I)    = RO(I,50)*Vmax(I)*THETA(I)/MU(I,50)
        REX(I)     = RO(I,50)*Vmax(I)*X(I,50)/MU(I,50)
        CF(I)      = TAUW(I)/(0.5*RO(I,50)*Vmax(I)**2./GC)
        DO 900 J=1,Jmax
          YPLUS(I,J) = Y(I,J)*VSHE(I)*RO(I,J)/MU(I,J)
          UPLUS(I,J) = VX(I,J)/VSHE(I)
900     CONTINUE
        DO 950 J=2,Jmax-1
          r=(Y(I,J+1)-Y(I,J)) / (Y(I,J)-Y(I,J-1))
          DUDY=(VX(I,J+1)+(r**2-1)*VX(I,J)-r**2*VX(I,J-1))/
     &         (r*(r+1.)*(Y(I,J)-Y(I,J-1)))
          TAU=TMU(I,J)*MU(I,J)*DUDY
          UV(I,J)=TAU/(RO(I,J)*VXfs**2)
950     CONTINUE
        UV(I,1)=0.
        UV(I,Jmax)=UV(I,Jmax-1)
        write(*,*) 'i=',i
C
1000  CONTINUE
      RETURN
      END
C========================================================================
      FUNCTION iFindReX(dRex)
      INCLUDE "flplpost.inc"
      REAL dREx
C-------------------------------------------------------
C     Find the first i-location with ReX = desired ReX
C-------------------------------------------------------
      iFindReX = 0
      DO 100 I=Ile,Imax
        IF ((ReX(I) .gt. dReX) .and. (iFindReX .eq. 0)) THEN
          IF (abs(ReX(I)-dReX) .lt. abs(dRex-ReX(I-1)) ) THEN
            iFindReX = I
          ELSE
            iFindReX = I-1
          ENDIF
        ENDIF
100   CONTINUE
      DO 500 M=1,IU(0)
        IF (iFindReX .eq. 0) THEN
          WRITE(IU(M),*) 'Unable to find ReX = ',dReX
        ELSE
          WRITE(IU(M),899) ReX(iFindReX),iFindReX
899       FORMAT('ReX = ',E12.4,' at I = ',I3)
        ENDIF
500   CONTINUE
      RETURN
      END
C========================================================================
      FUNCTION iFindReTH(dReTH)
      INCLUDE "flplpost.inc"
      REAL dReTH
C-------------------------------------------------------
C     Find the first i-location with ReX = desired ReX
C-------------------------------------------------------
      iFindReTH = 0
      DO 100 I=Ile,Imax
        IF ((ReTH(I) .gt. dReTH) .and. (iFindReTH .eq. 0)) THEN
          IF (abs(ReTH(I)-dReTH) .lt. abs(dReTH-ReTH(I-1)) ) THEN
            iFindReTH = I
          ELSE
            iFindReTH = I-1
          ENDIF
        ENDIF
100   CONTINUE
      DO 500 M=1,IU(0)
        IF (iFindReTH .eq. 0) THEN
          WRITE(IU(M),*) 'Unable to find ReTH = ',dReTH
        ELSE
          WRITE(IU(M),899) ReTH(iFindReTH),iFindReTH
899       FORMAT('ReTH = ',E12.4,' at I = ',I3)
        ENDIF
500   CONTINUE
      RETURN
      END
C========================================================================
      SUBROUTINE NewYGrid
C
      INCLUDE "flplpost.inc"
      REAL Ymax,Ymin,Snew
      INTEGER Jmaxold,IsOK
C
      Ymax = Y(1,Jmax)
      Ymin = Y(1,1)
      Snew = 1.0
      Jmaxold = Jmax
10    WRITE(*,*) 'Defining new vertical grid spacing'
20      CALL IUSERINPUT('  Enter new Jmax: ',Jmax,Jmax)
        IF (Jmax .gt. NY) THEN
          WRITE(*,*) '    --Jmax too big, try again:'
          GOTO 20
        ENDIF
        CALL USERINPUT('  Enter new Ymin: ',Ymin,Ymin)
        CALL USERINPUT('  Enter new Ymax: ',Ymax,Ymax)
        CALL USERINPUT('  Enter new S:    ',Snew,Snew)
C
        DO 500 J=1,Jmax
          DEXP= (Jmax-J) / (Jmax-1.0)
          DEN = (1.0 + ((Snew+1.0)/(Snew-1.0))**DEXP)
          Y(1,J) = Ymin + (Ymax - Ymin) * (1.0 - Snew + 2.*Snew/DEN)
500     CONTINUE
C
        WRITE(*,*) 'First few points off the wall are:'
        DO 700 J=1,10
          WRITE(*,609) J,Y(1,J)
609       FORMAT('   Y(1,',I2,',) = ',F10.7)
700     CONTINUE
        WRITE(*,*) 'Is this acceptable? (1=y/0=n) '
        READ(*,*) IsOK
      IF (IsOK .eq. 0) GOTO 10
C
      WRITE(10,*) 'New vertical grid spacing defined by:'
      WRITE(10,*) '  Jmax = ',Jmax
      WRITE(10,*) '  Ymin = ',Ymin
      WRITE(10,*) '  Ymax = ',Ymax
      WRITE(10,*) '  Snew = ',Snew
      WRITE(10,*) 'First few points off the wall are:'
      DO 900 J=1,10
        WRITE(10,609) J,Y(1,J)
900   CONTINUE
      WRITE(10,*)
C
      IF (Jmax .gt. Jmaxold) THEN
        DO 1200 J=Jmaxold,Jmax
          DO 1100 I=1,Imax
            X(I,J)=X(I,Jmaxold)
1100      CONTINUE
1200    CONTINUE
      ENDIF
C
      DO 1500 I=1,Imax
        DO 1400 J=1,Jmax
          DEXP= (Jmax-J) / (Jmax-1.0)
          DEN = (1.0 + ((Snew+1.0)/(Snew-1.0))**DEXP)
          Y(I,J) = Ymin + (Ymax - Ymin) * (1.0 - Snew + 2.*Snew/DEN)
1400    CONTINUE
1500  CONTINUE
C
      RETURN
      END
C========================================================================
      SUBROUTINE InitSol
C
      INCLUDE "flplpost.inc"
      REAL PTdef,TTdef,Xdef,AMACHdef,UNITSdef
      REAL PTinit,TTinit,Xinit,AMACHinit,UNITSinit
C
      PTdef = Pref
      TTdef = Tref
      Xdef  = Xref
      AMACHdef = AMACH
      UNITSdef = 1.0
C
      WRITE(*,*) 'Please enter the following flowfield information:'
      CALL USERINPUT('   Enter PT (lbf/ft^2): ',PTinit,PTdef)
      CALL USERINPUT('   Enter TT (Rankine):  ',TTinit,TTdef)
      CALL USERINPUT('   Enter Mach Number:   ',AMACHinit,AMACHdef)
      WRITE(*,*) '   Enter MULTIPLYING factor to convert grid units'
      CALL USERINPUT('     from (current) to (FT): ',UNITSinit,UNITSdef)
C
      FMinit = 1.+(GAMI/2.)*AMACHinit**2.
      TSinit = TTinit/FMinit
      PSinit = PTinit/(FMinit**GGM1)
      ROinit = PSinit/(RGAS*TSinit)
      Uinit  = AMACHinit*SQRT(GAMMA*RGAS*TSinit*GC)
      Einit  = PSinit/GAMI + 0.5*ROinit*Uinit**2./GC
C
      DO 200 J=1,Jmax
        DO 100 I=1,Imax
          Q(I,J,1) = ROinit/Rref
          Q(I,J,2) = ROinit*Uinit/(Rref*Aref)
          Q(I,J,3) = 0.0
          Q(I,J,4) = Einit/(GAMMA*Pref)
          QT(I,J,1)= 0.0
          QT(I,J,2)= 0.0
          QT(I,J,3)= 0.0
100     CONTINUE
200   CONTINUE
C
      DO 1000 M=1,IU(0)
        WRITE(IU(M),*) 'Initializing flow properties to:',
     &    '       PARC values are:'
        WRITE(IU(M),909) 'Pt  = ',PTinit,   ' lbf/ft^2','Q1 = ',Q(1,1,1)
        WRITE(IU(M),909) 'Tt  = ',TTinit,   ' Rankine ','Q2 = ',Q(1,1,2)
        WRITE(IU(M),909) 'M   = ',AMACHinit,'         ','Q3 = ',Q(1,1,3)
        WRITE(IU(M),909) 'P   = ',PSinit,   ' lbf/ft^2','Q4 = ',Q(1,1,4)
        WRITE(IU(M),909) 'T   = ',TSinit,   ' Rankine ','QT1= ',0.0
        WRITE(IU(M),909) 'rho = ',ROinit,   ' lbm/ft^3','QT2= ',0.0
        WRITE(IU(M),919) 'e   = ',Einit,    ' lbf/ft^2','QT3= ',0.0
        WRITE(IU(M),*) 
909     FORMAT(3x,A,F13.6,A,11x,A,E13.6)
919     FORMAT(3x,A,E13.6,A,11x,A,E13.6)
1000  CONTINUE
C
      RETURN
      END
C========================================================================
      SUBROUTINE UserInput(TEXT,X,Xdef)
C
      REAL X,Xdef
      CHARACTER TEXT*30
C
      WRITE(*,*) TEXT,'  [-1 ==> DEFAULT value of:  ',Xdef,']'
      READ(*,*)  X
      IF (X .LT. 0.0) X=Xdef
      Xdef=X
      RETURN
      END
C========================================================================
      SUBROUTINE IUserInput(TEXT,I,Idef)
C
      INTEGER I,Idef
      CHARACTER TEXT*30
C
      WRITE(*,*) TEXT,'  [-1 ==> DEFAULT value of:  ',Idef,']'
      READ(*,*)  I
      IF (I .LT. 0) I=Idef
      Idef=I
      RETURN
      END
C========================================================================
