#!/bin/bash

########################################
#  Enter Fortran compiler here
########################################
FC=ifort

########################################
#  Create and run CFPOST script
########################################
echo "
grid FP4p5.cgd
solution FP4p5.cfl

zone 1
surface u 2
cut at z 0.1
!variable x; Cfx; delta; THETA; delta*
variable x; Cfx
genplot output FP4p5.Cf.out

clear all

zone 1
surface u 4
cut at x 1.791667
variable y;x
genplot output FP4p5.BL-x.out
clear variable
variable y;u
genplot output FP4p5.BL-u.out
clear variable
variable y;v
genplot output FP4p5.BL-v.out
clear variable
variable y;rho
genplot output FP4p5.BL-rho.out
clear variable
variable y;p
genplot output FP4p5.BL-p.out
clear variable
variable y;T
genplot output FP4p5.BL-T.out
clear variable
variable y;mu
genplot output FP4p5.BL-mu.out
clear variable
variable y;mut muinf
genplot output FP4p5.BL-mut.out

exit
" > FP4p5.cfpost.com

cfpost_prod.exe < FP4p5.cfpost.com

########################################
#  Remove unnecessary text from data files
########################################
grep -v -e 'Exported' -e 'Unstructured' -e 'x' -e 'y' -e 'Cut' FP4p5.BL-x.out > FP4p5.BL-x.out2
grep -v -e 'Exported' -e 'Unstructured' -e 'y' -e 'u' -e 'Cut' FP4p5.BL-u.out > FP4p5.BL-u.out2
grep -v -e 'Exported' -e 'Unstructured' -e 'y' -e 'v' -e 'Cut' FP4p5.BL-v.out > FP4p5.BL-v.out2
grep -v -e 'Exported' -e 'Unstructured' -e 'y' -e 'rho' -e 'Cut' FP4p5.BL-rho.out > FP4p5.BL-rho.out2
grep -v -e 'Exported' -e 'Unstructured' -e 'y' -e 'p' -e 'Cut' FP4p5.BL-p.out > FP4p5.BL-p.out2
grep -v -e 'Exported' -e 'Unstructured' -e 'y' -e 'T' -e 'Cut' FP4p5.BL-T.out > FP4p5.BL-T.out2
grep -v -e 'Exported' -e 'Unstructured' -e 'y' -e 'mu' -e 'Cut' FP4p5.BL-mu.out > FP4p5.BL-mu.out2
grep -v -e 'Exported' -e 'Unstructured' -e 'y' -e 'mut' -e 'Cut' FP4p5.BL-mut.out > FP4p5.BL-mut.out2

grep -v -e 'Exported' -e 'Unstructured' -e 'x' -e 'Cfx' -e 'Cut' FP4p5.Cf.out > FP4p5.Cf.out2

########################################
#  This Fortran code consolidates the data into Cf and BL files
########################################
echo "
c234567890123456789012345678901234567890123456789012345678901234567890
c
      Real x(1000),y(1000),u(1000),v(1000),rho(1000),p(1000),T(1000),
     &     mu(1000),mut(1000),
     &     x2(1000),Cf(1000)
      Integer k,kmax,i,imax(10),imaxsum
c
c---Consolidate BL data
11    FORMAT(E15.8,T20,E15.8,T40,E15.8,T60,E15.8,T80,E15.8,
     &       T100,E15.8,T120,E15.8,T140,E15.8,T160,E15.8)
c
      open(unit=21,file='FP4p5.BL-x.out2',form='formatted')
      open(unit=22,file='FP4p5.BL-u.out2',form='formatted')
      open(unit=23,file='FP4p5.BL-v.out2',form='formatted')
      open(unit=24,file='FP4p5.BL-rho.out2',form='formatted')
      open(unit=25,file='FP4p5.BL-p.out2',form='formatted')
      open(unit=26,file='FP4p5.BL-T.out2',form='formatted')
      open(unit=27,file='FP4p5.BL-mu.out2',form='formatted')
      open(unit=28,file='FP4p5.BL-mut.out2',form='formatted')
c
      read(21,*) kmax
      read(22,*) kmax
      read(23,*) kmax
      read(24,*) kmax
      read(25,*) kmax
      read(26,*) kmax
      read(27,*) kmax
      read(28,*) kmax
c
      do k=1,kmax
        read(21,*) imax(k)
        read(22,*) imax(k)
        read(23,*) imax(k)
        read(24,*) imax(k)
        read(25,*) imax(k)
        read(26,*) imax(k)
        read(27,*) imax(k)
        read(28,*) imax(k)
        if (k .EQ. 1) then
	  imaxsum=0
	endif
        do i=1,imax(k)
	  read(21,*) y(i+imaxsum),x(i+imaxsum)
	  read(22,*) y(i+imaxsum),u(i+imaxsum)
	  read(23,*) y(i+imaxsum),v(i+imaxsum)
	  read(24,*) y(i+imaxsum),rho(i+imaxsum)
	  read(25,*) y(i+imaxsum),p(i+imaxsum)
	  read(26,*) y(i+imaxsum),T(i+imaxsum)
	  read(27,*) y(i+imaxsum),mu(i+imaxsum)
	  read(28,*) y(i+imaxsum),mut(i+imaxsum)
        enddo
	imaxsum=imaxsum+imax(k)
      enddo
c
      close(21)
      close(22)
      close(23)
      close(24)
      close(25)
      close(26)
      close(27)
      close(28)
c
      open(unit=41,file='FP4p5.BL.out3',form='formatted')
      do i=1,imaxsum
        write(41,11) y(i),x(i),u(i),v(i),rho(i),p(i),T(i),mu(i),mut(i)
      enddo
      close(41)
c
c---Consolidate Cf data
12    FORMAT(E15.8,T20,E15.8)
      open(unit=31,file='FP4p5.Cf.out2',form='formatted')
      read(31,*) kmax
      do k=1,kmax
        read(31,*) imax(k)
        if (k .EQ. 1) then
	  imaxsum=0
	endif
        do i=1,imax(k)
	  read(31,*) x2(i+imaxsum),Cf(i+imaxsum)
        enddo
	imaxsum=imaxsum+imax(k)
      enddo
      close(31)
      open(unit=42,file='FP4p5.Cf.out3',form='formatted')
      do i=1,imaxsum
        write(42,12) x2(i),Cf(i)
      enddo
      close(42)
c
      end
" > FP4p5consolidate.f

$FC -o FP4p5consolidate.exe FP4p5consolidate.f
./FP4p5consolidate.exe

########################################
# Sort the data in the Cf and BL files
########################################
sort -g -u FP4p5.BL.out3 > FP4p5.BL.out4

sort -g -u FP4p5.Cf.out3 > FP4p5.Cf.out4

########################################
#  This Fortran code computes the BL thickness
########################################
echo "
c234567890123456789012345678901234567890123456789012345678901234567890
c
      Real x(1000),y(1000),u(1000),v(1000),rho(1000),p(1000),T(1000),
     &     mu(1000),mut(1000),Vs(1000),
     &     Uinf,BLfactor,nuinf,
     &     ydelta,xdelta,Rexdelta,udelta,vdelta,Vsdelta,rhodelta,
     &     theta,deltastar,sumy,sumrhoV,sumrhoV2
      Integer i,imax,idelta
c
11    FORMAT(E15.8,T20,E15.8,T40,E15.8,T60,E15.8,T80,E15.8)
c
      Uinf=2307.38
      nuinf=(0.814118E-07)/(0.751033E-04)
      BLfactor=0.99
c
      open(unit=21,file='FP4p5.BL.out4',form='formatted')
c
      imax=81
      do i=1,imax
        read(21,*) y(i),x(i),u(i),v(i),rho(i),p(i),T(i),mu(i),mut(i)
	Vs(i)=sqrt(u(i)**2+v(i)**2)
      enddo
      udelta=BLfactor*Uinf
      do i=1,(imax-1)
        if (u(i) .LE. udelta) then
	  if (u(i+1) .GT. udelta) then
	    ydelta=((y(i+1)-y(i))/(u(i+1)-u(i)))*(udelta-u(i))+y(i)
            xdelta=((x(i+1)-x(i))/(y(i+1)-y(i)))*(ydelta-y(i))+x(i)
	    Rexdelta=xdelta*Uinf/nuinf
            vdelta=((v(i+1)-v(i))/(y(i+1)-y(i)))*(ydelta-y(i))+v(i)
            Vsdelta=((Vs(i+1)-Vs(i))/(y(i+1)-y(i)))*(ydelta-y(i))+Vs(i)
            rhodelta=((rho(i+1)-rho(i))/(y(i+1)-y(i)))*
     &               (ydelta-y(i))+rho(i)
	    idelta=i
	  endif
	endif
      enddo
c
      sumy=0.0
      sumrhoV=0.0
      sumrhoV2=0.0
      do i=1,(idelta-1)
        sumy=sumy+(y(i+1)-y(i))*(1+1)/2
        sumrhoV=sumrhoV+(y(i+1)-y(i))*(rho(i+1)*Vs(i+1)+rho(i)*Vs(i))/2
        sumrhoV2=sumrhoV2+(y(i+1)-y(i))*
     &           (rho(i+1)*Vs(i+1)**2+rho(i)*Vs(i)**2)/2
      enddo
      sumy=sumy+(ydelta-y(idelta))*(1+1)/2
      sumrhoV=sumrhoV+(ydelta-y(idelta))*
     &        (rhodelta*Vsdelta+rho(idelta)*Vs(idelta))/2
      sumrhoV2=sumrhoV2+(ydelta-y(idelta))*
     &         (rhodelta*Vsdelta**2+rho(idelta)*Vs(idelta)**2)/2
      theta=sumy-(1/(rhodelta*Vsdelta))*sumrhoV
      deltastar=(1/(rhodelta*Vsdelta))*sumrhoV-
     &          (1/(rhodelta*Vsdelta**2))*sumrhoV2
c
      open(unit=31,file='FP4p5.delta.out',form='formatted')
      write(31,11) xdelta,Rexdelta,ydelta,theta,deltastar
      close(31)
c
      end
" > FP4p5delta.f

$FC -o FP4p5delta.exe FP4p5delta.f
./FP4p5delta.exe

########################################
#  This Fortran code computes Cf data
########################################
echo "
c234567890123456789012345678901234567890123456789012345678901234567890
c
      Real x(1000),Cf(1000),Rex(1000),xdelta,
     &     Rexdelta,delta,theta,deltastar,Cfdelta,
     &     Uinf,nuinf
      Integer i,imax
      Character*45 charline1
      Character*43 charline2
      Character*45 charline3
      Character*9 charline4
      Character*9 charline5
      
c
11    FORMAT(E15.8,T20,E15.8,T40,E15.8,T60,E15.8,T80,E15.8,
     &       T100,E15.8)
c
c      imax=141
      imax=31
      Uinf=2307.38
      nuinf=(0.814118E-07)/(0.751033E-04)
c
      open(unit=21,file='FP4p5.Cf.out4',form='formatted')
      open(unit=22,file='FP4p5.delta.out',form='formatted')
c
      open(unit=31,file='FP4p5.Cf.dat',form='formatted')
      charline1 = 'TITLE     = \"Mach 4.5 Flow Over a Flat Plate:'
      charline2 = 'Skin Friction and Integrated BL Quantities\"'
      write(31,'(A,A)') charline1,charline2
      write(31,'(A)') 'VARIABLES = \"x [ft]\"'
      write(31,'(A)') '\"Rex\"'
      write(31,'(A)') '\"Cfx\"'
      write(31,'(A)') '\"delta [ft]\"'
      write(31,'(A)') '\"theta [ft]\"'
      write(31,'(A)') '\"delta* [ft]\"'
      write(31,'(A)') ''
      charline3 = 'ZONE T=\"Structured Grid, Unstructured Solver,'
      charline4 = ' SST\", I='
      charline5 = ', F=POINT'
      write(31,'(A,A,I4,A)') charline3,charline4,(imax+1),charline5
c
      do i=1,imax
        read(21,*) x(i),Cf(i)
	Rex(i)=x(i)*Uinf/nuinf
      enddo
      read(22,*) xdelta,Rexdelta,delta,theta,deltastar
      close(21)
      close(22)
c
      do i=1,imax
        if (x(i+1) .GT. xdelta) then
	  if (x(i) .LE. xdelta) then
	    write(31,11) x(i),Rex(i),Cf(i),0,0,0
	    Cfdelta=((Cf(i+1)-Cf(i))/(x(i+1)-x(i)))*(xdelta-x(i))+Cf(i)
	    write(31,11) xdelta,Rexdelta,Cfdelta,delta,theta,deltastar
	  else
	    write(31,11) x(i),Rex(i),Cf(i),0,0,0
	  endif	    
	else
	  write(31,11) x(i),Rex(i),Cf(i),0,0,0
	endif
      enddo
      close(23)
c
      end
" > Cfsort.f

$FC -o Cfsort.exe Cfsort.f
./Cfsort.exe

########################################
#  This Fortran code computes BL data
########################################
echo "
c234567890123456789012345678901234567890123456789012345678901234567890
c
      Real x(1000),y(1000),u(1000),v(1000),rho(1000),p(1000),T(1000),
     &     mu(1000),mut(1000)
      Integer i,imax
      Character*45 charline1
      Character*42 charline2
      Character*45 charline3
      Character*9 charline4
      Character*9 charline5
c
11    FORMAT(E15.8,T20,E15.8,T40,E15.8,T60,E15.8,T80,E15.8,
     &       T100,E15.8,T120,E15.8,T140,E15.8,T160,E15.8)
c
      imax=81
c
      open(unit=21,file='FP4p5.BL.out4',form='formatted')
      do i=1,imax
        read(21,*) y(i),x(i),u(i),v(i),rho(i),p(i),T(i),mu(i),mut(i)
      enddo
      close(21)
c
      open(unit=31,file='FP4p5.BL.dat',form='formatted')
      charline1 = 'TITLE     = \"Mach 4.5 Flow Over a Flat Plate:'
      charline2 = 'Boundary Layer Profile at Rex=3.8 million\"'
      write(31,'(A,A)') charline1,charline2
      write(31,'(A)') 'VARIABLES = \"x [ft]\"'
      write(31,'(A)') '\"y [ft]\"'
      write(31,'(A)') '\"u [ft/s]\"'
      write(31,'(A)') '\"v [ft/s]\"'
      write(31,'(A)') '\"rho [slug/ft^3]\"'
      write(31,'(A)') '\"p [psi]\"'
      write(31,'(A)') '\"T [deg R]\"'
      write(31,'(A)') '\"mu [slug/(ft*s)]\"'
      write(31,'(A)') '\"mut/muinf\"'
      write(31,'(A)') ''
c
      charline3 = 'ZONE T=\"Structured Grid, Unstructured Solver,'
      charline4 = ' SST\", I='
      charline5 = ', F=POINT'
      write(31,'(A,A,I3,A)') charline3,charline4,imax,charline5
c
      do i=1,imax
        write(31,11) x(i),y(i),u(i),v(i),rho(i),p(i),T(i),mu(i),mut(i)
      enddo
c
      close(31)
c
      end
" > BLsort.f

$FC -o BLsort.exe BLsort.f
./BLsort.exe

########################################
#  Extract max N-S residuals from LIS file (assumes single zone)
########################################
grep -a RES FP4p5.lis | grep -v -e 'Gradient' | awk '{ print $2 "    " $7 "    " $8 }'  > FP4p5.NSresid.out
grep -a RES FP4p5.lis | grep -v -e 'Gradient' | awk '{ print $2 } ' | awk '{ c += length + 1 } END { print NR }' > FP4p5.imax.out

########################################
#  This Fortran code formats the N-S residuals data
########################################
echo "
c234567890123456789012345678901234567890123456789012345678901234567890
c This code formats N-S residuals data
c
      Real big(100000),l2(100000)
      Integer i,iter(100000),imax
      Character*45 charline1
      Character*24 charline2
      Character*45 charline3
      Character*9 charline4
      Character*9 charline5
c
11    FORMAT(I10,T20,E15.8,T40,E15.8)
c
      open(unit=20,file='FP4p5.imax.out',form='formatted')
      read(20,*) imax
      close(20)
c
      open(unit=21,file='FP4p5.NSresid.out',form='formatted')
      do i=1,imax
        read(21,*) iter(i),big(i),l2(i)
      enddo
      close(21)
c
      open(unit=31,file='FP4p5.NSresid.dat',form='formatted')
      charline1 = 'TITLE     = \"Mach 4.5 Flow Over a Flat Plate:'
      charline2 = 'Navier-Stokes Residuals\"'
      write(31,'(A,A)') charline1,charline2
      write(31,'(A)') 'VARIABLES = \"iteration\"'
      write(31,'(A)') '\"N-S residual\"'
      write(31,'(A)') '\"L2 residual\"'
      write(31,'(A)') ''
c
      charline3 = 'ZONE T=\"Structured Grid, Unstructured Solver,'
      charline4 = ' SST\", I='
      charline5 = ', F=POINT'
      write(31,'(A,A,I7,A)') charline3,charline4,imax,charline5
c
      do i=1,imax
        write(31,11) iter(i),big(i),l2(i)
      enddo
c
      close(31)
c
      end
" > NSresid.f

$FC -o NSresid.exe NSresid.f
./NSresid.exe

########################################
#  Clean up intermediate output files
########################################
rm -f FP4p5.BL*.out*
rm -f FP4p5.Cf.out*
rm -f FP4p5.delta.out
rm -f FP4p5.NSresid.out FP4p5.imax.out
rm -f *.jou *.f *.exe
