DISPLAY
  This case is similar to case 159 except that the jet discharge
  is supersonic and the flow is three-dimensional because the
  discharge nozzle is square. 
  
  The option to simulate a 2d round jet is also provided. 
  The calculation is made for an inlet Mach number of 2 and a 
  jet-to-ambient static pressure ratio of 1.8. 
  
  The free=stream Mach number is essentially zero.
  ENDDIS
  PHOTON USE
   P
   PARPHI
 
 
   VEC X 1 SH
   pau;cl
   con mach x 1 fi;.1
   pau;cl
   con p1 x 1 fi;.1
   pau cl;
   con tmp1 x 1 fi;.1
  ENDUSE
 
REAL(AIN,CP,GAM,GM1,PTOT,HTOT,RHTOT,PAMB,HAMB,MIN,PIN,HIN,RHOIN)
REAL(WIN,EPSIN,TKEIN,DTF,DN,RAD,PRAT,RCON,RGAM,TTOT,TIN,UAC,RHAMB)
REAL(FLOWIN,AOIN,UREF,DYGDZ,DXGDZ,TAMB,WAMB,MAMB)
CHAR(CTURB);BOOLEAN(GEXPAN)
   **  GEXPAN=T activates a linear y-grid expansion with z
       Not yet considered or tested in this Q1
GEXPAN=F
 
   **  Gas properties
GAM=1.4;GM1=GAM-1.;RCON=1.0;CP=RCON*GAM/GM1;RGAM=1./GAM
 
   **  INLET CONDITIONS & NOZZLE DIAMETER
PTOT=1.0;RHTOT=1.0;TTOT=1.0;MIN=2.0;DN=1.0
HTOT=CP*TTOT
 
RAD=0.5*DN
PIN=PTOT/(1.+GM1*MIN*MIN/2.)**(GAM/GM1)
RHOIN=RHTOT*(PIN/PTOT)**RGAM;WIN=MIN*(GAM*PIN/RHOIN)**0.5
TIN=PIN/(RCON*RHOIN);HIN=CP*TIN
AOIN=(GAM*PTOT/RHTOT)*0.5;UREF=AOIN/GAM**0.5
    ** Set static pressure ratio, i.e. PIN/PAMB
PRAT=1.8
    ** Set free-stream Mach number
MAMB=1.E-3
    ** Ambient conditions
PAMB=PIN/PRAT;HAMB=HTOT
TAMB=HAMB/CP;RHAMB=PAMB/(RCON*TAMB)
WAMB=MAMB*(GAM*PAMB/RHAMB)**0.5
   ** SJET=T for a 3d square jet
          =F for a 2d round jet
BOOLEAN(SJET)
MESG( Enter T for a 3d square jet (default)
MESG( Enter F for a 2d round jet
READVDU(SJET,BOOLEAN,T)
    GROUP 1. Run title and other preliminaries
    GROUP 2. Transience; time-step specification
PARAB=T
    GROUP 3. X-direction grid specification
IF(SJET) THEN
+ TEXT(3D Supersonic Underexpanded Square Jet
+ CARTES=T;AIN=DN*RAD
+ NREGX=2
+ IREGX=1;GRDPWR(X, 5,  RAD,1.0)
+ IREGX=2;GRDPWR(X,10,   DN,1.0)
IF(GEXPAN) THEN
+ DXGDZ=0.0875;AZXU=1.0;ZWADD=(DN+RAD)/DXGDZ
ENDIF
ELSE
+ TEXT(2D SUPERSONIC UNDEREXPANDED ROUND JET:160
+ CARTES=F;NX=1;XULAST=0.1;AIN=0.5*XULAST*RAD*RAD
ENDIF
TITLE
    GROUP 4. Y-direction grid specification
IF(SJET) THEN
+ NREGY=3
+ IREGY=1;GRDPWR(Y,10, DN,1.0)
+ IREGY=2;GRDPWR(Y,10, DN,1.0)
+ IREGY=3;GRDPWR(Y,10, DN,1.0)
ELSE
+ NREGY=2
+ IREGY=1;GRDPWR(Y,10,RAD,1.0)
+ IREGY=2;GRDPWR(Y,20,2*RAD,1.0)
IF(GEXPAN) THEN
+ DYGDZ=0.069927;AZYV=1.0;ZWADD=DN/DYGDZ
ENDIF
ENDIF
    GROUP 5. Z-direction grid specification
NZ=500;AZDZ=PROPY;DZW1=0.01
    GROUP 7. Variables stored, solved & named
SOLVE(P1,V1,W1)
IF(NX.GT.1) THEN
+ SOLVE(U1);SOLUTN(U1,P,P,Y,Y,P,P);TERMS(U1,Y,Y,N,P,P,P)
ENDIF
SOLUTN(V1,P,P,Y,Y,P,P);SOLUTN(W1,P,P,Y,Y,P,P)
  ** Provide storage for the density.
STORE(RHO1,MACH,ENUT,LEN1,H1,TMP1)
IPARAB=5;STORE(MACZ)
   ** RMACHZ=1 unless enut=0.
RMACHZ=10.0
 
MESG( Enter the required turbulence model:
MESG(  KE   -  Standard k-e model (Default)
MESG(  LAM  -  Laminar model
MESG(
READVDU(CTURB,CHAR,KE)
CASE :CTURB: OF
WHEN KE,2
+ MESG(Standard k-e model
+ TURMOD(KEMODL);KELIN=3
WHEN LAM,3
+ MESG(Laminar model
+ ENUT=0.
  + RMACHZ=1.E3
ENDCASE
    GROUP 8. Terms (in differential equations) & devices
TERMS(V1,Y,Y,N,P,P,P)
DIFCUT=0;DENPCO=T
    GROUP 9. Properties of the medium (or media)
RHO1=IDEALGAS;RHO1B=1./RCON;PRESS0=0.0
DRH1DP=IDEALGAS;RHO1C=RGAM
TMP1=VARSTAGH; CP1=CP
ENUL=1.8E-5/(UREF*DN)
    GROUP 11. Initialization of variable or porosity fields
FIINIT(U1)=0.0; FIINIT(V1)=0.0
FIINIT(W1)=WIN; FIINIT(MACZ)= MIN
FIINIT(RHO1)=RHOIN; FIINIT(P1)=PIN
FIINIT(H1)=HTOT; FIINIT(TMP1)=TIN
TKEIN=(0.05*WIN)**2;EPSIN=0.1643*TKEIN**1.5/(0.1*RAD/DN)
FIINIT(EP)=EPSIN; FIINIT(KE)=TKEIN
INIADD=F
IF(SJET) THEN
+ PATCH(INITFS1,INIVAL,1,NX,#1,#1,1,1,1,1)
+ COVAL(INITFS1,W1,0.0,WAMB);COVAL(INITFS1,P1,0.0,PAMB)
+ COVAL(INITFS1,TMP1,0.0,TAMB);COVAL(INITFS1,RHO1,0.0,RHAMB)
+ COVAL(INITFS1,MACZ,0.0,MAMB)
 
+ PATCH(INITFS3,INIVAL,1,NX,#3,#3,1,1,1,1)
+ COVAL(INITFS3,W1,0.0,WAMB);COVAL(INITFS3,P1,0.0,PAMB)
+ COVAL(INITFS3,TMP1,0.0,TAMB);COVAL(INITFS3,RHO1,0.0,RHAMB)
+ COVAL(INITFS3,MACZ,0.0,MAMB)
 
+ PATCH(INITFS4,INIVAL,#2,#2,#1,#3,1,1,1,1)
+ COVAL(INITFS4,W1,0.0,WAMB);COVAL(INITFS4,P1,0.0,PAMB)
+ COVAL(INITFS4,TMP1,0.0,TAMB);COVAL(INITFS4,RHO1,0.0,RHAMB)
+ COVAL(INITFS4,MACZ,0.0,MAMB)
+ COVAL(INITFS4,p1,0.0,PAMB)
ELSE
+ PATCH(INITFS3,INIVAL,1,NX,#2,#2,1,1,1,1)
+ COVAL(INITFS3,W1,0.0,WAMB);COVAL(INITFS3,P1,0.0,PAMB)
+ COVAL(INITFS3,TMP1,0.0,TAMB);COVAL(INITFS3,RHO1,0.0,RHAMB)
+ COVAL(INITFS3,MACZ,0.0,MAMB)
ENDIF
    GROUP 13. Boundary conditions and special sources
IF(SJET) THEN
+ INLET(IN1,LOW,1,NX,#1,#1,1,1,1,1)
+ VALUE(IN1,P1,RHAMB*WAMB)
+ VALUE(IN1,W1,WAMB)
 
+ INLET(IN2,LOW,#1,#1,#2,#2,1,1,1,1)
+ VALUE(IN2,P1,RHOIN*WIN)
+ VALUE(IN2,W1,WIN)
+ VALUE(IN2,KE,TKEIN);VALUE(IN2,EP,EPSIN)
 
+ INLET(IN3,LOW,1,NX,#3,#3,1,1,1,1)
+ VALUE(IN3,P1,RHAMB*WAMB)
+ VALUE(IN3,W1,WAMB)
 
+ INLET(IN4,LOW,#2,#2,#2,#2,1,1,1,1)
+ VALUE(IN4,P1,RHAMB*WAMB)
+ VALUE(IN4,W1,WAMB)
 
+ PATCH(EB,EAST,NX,NX,1,NY,1,NZ,1,1)
+ COVAL(EB,P1,FIXVAL,PAMB);COVAL(EB,W1,ONLYMS,WAMB)
 
+ PATCH(SB,SOUTH,1,NX,1,1,1,NZ,1,1)
+ COVAL(SB,P1,FIXVAL,PAMB);COVAL(SB,W1,ONLYMS,WAMB)
ELSE
+ INLET(IN,LOW,1,NX,#1,#1,1,1,1,1)
+ VALUE(IN,P1,RHOIN*WIN);VALUE(IN,W1,WIN)
+ VALUE(IN,KE,TKEIN);VALUE(IN,EP,EPSIN)
 
+ INLET(INFS,LOW,1,NX,#2,#2,1,1,1,1)
+ VALUE(INFS,P1,RHAMB*WAMB);VALUE(INFS,W1,WAMB)
+ VALUE(INFS,KE,TKEIN);VALUE(INFS,EP,EPSIN)
ENDIF
 
PATCH(NB,NORTH,1,NX,NY,NY,1,NZ,1,1)
COVAL(NB,P1,FIXVAL,PAMB);COVAL(NB,W1,ONLYMS,WAMB)
    GROUP 16. Termination of iterations
IF(SJET) THEN
+ LITHYD=15
ELSE
+ LITHYD=10
ENDIF
TSTSWP=-1;IPLTL=LITHYD
LITER(U1)=30; LITER(V1)=30; LITER(W1)=30; LITER(P1)=100
LITER(KE)=30; LITER(EP)=30
 
FLOWIN=RHOIN*WIN*AIN
SELREF=F;RESREF(P1)=1.E-12*FLOWIN
RESREF(W1)=1.E-12*FLOWIN*WIN
RESREF(KE)=1.E-12*FLOWIN*TKEIN;RESREF(EP)=1.E-12*FLOWIN*EPSIN
RESREF(V1)=RESREF(W1);RESREF(U1)=RESREF(W1)
 
    GROUP 17. Under-relaxation devices
DTF=10.*DZW1*YVLAST/(WIN+UAC)
RELAX(U1,LINRLX,0.5); RELAX(V1,LINRLX,0.5)
RELAX(W1,LINRLX,0.5); RELAX(P1,LINRLX,0.5)
RELAX(H1,LINRLX,0.5); RELAX(TMP1,LINRLX,0.5)
RELAX(KE,LINRLX,0.3); RELAX(EP,LINRLX,0.3)
RELAX(MACZ,LINRLX,0.8)
    GROUP 18. Limits on variables or increments to them
VARMIN(P1)=0.0001*PTOT
VARMIN(RHO1)=0.01*RHTOT;VARMAX(RHO1)=5.0*RHTOT
VARMAX(P1)=5.0*PTOT
    GROUP 19. Data communicated by satellite to GROUND
    GROUP 20. Preliminary print-out
OUTPUT(MACH,P,P,P,P,Y,P)
    GROUP 21. Print-out of variables
PATCH(constz,profil,1,1,1,ny,1,1,1,1)
COVAL(constz,p1,0.0,0.0);COVAL(constz,w1,0.0,0.0)
COVAL(constz,rho1,0.0,0.0)
PATCH(constz2,profil,1,1,1,ny,1,1,1,1)
COVAL(constz2,macz,0.0,0.0);COVAL(constz2,v1,0.0,0.0)
    GROUP 22. Spot-value print-out
IXMON=1;IYMON=2*NY/3;IZMON=1
    GROUP 23. Field print-out and plot control
ITABL=3;NPLT=1;NYPRIN=1;NZPRIN=NZ/5
IF(NZ.GT.1) THEN
+ IDISPA=nz/20;IDISPB=1;IDISPC=NZ
ENDIF