CCM Test: Laminar flow through a gradual expansion.
  **************************************************************
  DISPLAY
  ----------------------------------------------------------
   This  case  concerns  flow  inside  a channel with smooth
   expansion. The  geometry for  it, proposed  by P.  Roache
   (P.  Roache,  "Scaling  of  High  Reynolds  Number Weakly
   Separated  Channel  Flows",  in  Numerical  and  Physical
   Aspects  of  Aerodynamic  Flows,  Spriger-Verlag,  1981),
   depends on the value of the Reynolds number (REUNO).  The
   channel becomes longer and straighter as Reynolds  number
   increases.
 
   NOTE, that the simulation for Re = 10 needs smaller value
   of linear relaxation for P1 (i.e 0.05), than default.
 
   User can switch from the default colocated  computational
   algorithm (CCM)  to the  staggered one  (STAG) by setting
   LCCM= F.
  ----------------------------------------------------------
  ENDDIS
L(PAUSE
  **************************************************************
BOOLEAN(LCCM,LUNIF);  LCCM= T;  LUNIF= F
  **************************************************************
  PHOTON USE
   p ; ; ; ; ;
 
   msg Computational Domain:
   gr i 1
   msg Press Any Key to Continue...
   pause
   cl
   set vec av off
   msg Velocity Vectors:
   vec i 1 sh
   msg Press Any Key to Continue...
   pause
   cl
   msg Contours of Pressure:
   con p1 i 1 fi;0.0001
   msg Press E  to exit PHOTON ...
  ENDUSE
  **************************************************************
    GROUP 1. Run title and other preliminaries
IF(LCCM) THEN
+ TEXT(CCM : 2D-duct with F(Re)-expansion.
ELSE
+ TEXT(STAG: 2D-duct with F(Re)-expansion.
+ NONORT= T
ENDIF
TITLE
INTEGER(NY1,NY2)
REAL(REYNO,UIN,LZ,LY,DZZ,DY,DZ1,DZ2,TAN2,ZCUR,YCUR,WCR,TARG,YWL)
  ** Problem definition:
REYNO= 15.0;   UIN = 1.0;        LZ  = REYNO/3.0;
NY1  = 8;      NY2 = 8;          NZ  = 32;     NY= NY1+NY2
DZZ  = LZ/NZ;  DZ1 = DZZ*0.3;    DZ2 = DZZ*2
ZCUR = LZ;     TAN2= TANH(2.0);  TARG=-8.0
LY   = 1.0-(TANH(TARG)-TAN2)/2.0
    GROUP 6. Body-fitted coordinates or grid distortion
BFC= T; GSET(D,NX,NY,NZ,1.0,LY,LZ)
GSET(P,P1,0.0,-1.0,   0.0);  GSET(P,P2,0.0,1.0,   0.0)
GSET(P,P3,0.0,-1.0,   DZ1);  GSET(P,P4,0.0,1.0,   DZ1)
GSET(P,P5,0.0, -LY,LZ-DZ2);  GSET(P,P6,0.0, LY,LZ-DZ2)
GSET(P,P7,0.0, -LY,    LZ);  GSET(P,P8,0.0, LY,    LZ)
GSET(V,YWT,S,P4,SPLINE)
DO II =1,5
+ ZCUR= (LZ/2-DZ1)/5*II + DZ1; TARG= 2.0-10.*ZCUR/LZ
+ YCUR= 1.0-(TANH(TARG)-TAN2)/2.0
+ GSET(V,0.0,YCUR,ZCUR)
ENDDO
GSET(V,YWT,E,P6,0.0,0.0,1.0,0.0)
GSET(L,LWT,P4,P6,NZ-2,1.2CRVYWT)
GSET(V,YWB,S,P3,SPLINE)
DO II =1,5
+ ZCUR= (LZ/2-DZZ)/5*II + DZ1; TARG= 2.0-10.*ZCUR/LZ
+ YCUR= -(1.0-(TANH(TARG)-TAN2)/2.0)
+ GSET(V,0.0,YCUR,ZCUR)
ENDDO
GSET(V,YWB,E,P5,0.0,0.0,1.0,0.0)
GSET(L,LWB,P3,P5,NZ-2,1.2CRVYWB)
GSET(L,L12,P1,P2,NY);  GSET(L,L24,P2,P4, 1)
GSET(L,L43,P4,P3,NY);  GSET(L,L31,P3,P1, 1)
GSET(L,L65,P6,P5,NY);  GSET(L,L87,P8,P7,NY)
GSET(L,L66,P6,P8, 1);  GSET(L,L75,P7,P5, 1)
GSET(F,F1,P1,-,P2,-,P4,-,P3,-); GSET(M,F1,+J+K,1,1, 1)
GSET(F,F2,P3,-,P4,-,P6,-,P5,-); GSET(M,F2,+J+K,1,1, 2)
GSET(F,F3,P5,-,P6,-,P8,-,P7,-); GSET(M,F3,+J+K,1,1,NZ)
GSET(C,I:NX+1:,F,I1,1,NY,1,NZ,+,1.0,0.0,0.0,INC,1.0)
GVIEW(X); VIEW
    GROUP 7. Variables stored, solved & named
SOLVE(P1,V1,W1); SOLUTN(P1,Y,Y,Y,N,N,N)
IF(LCCM) THEN
L($F150)
ENDIF
    GROUP 9. Properties of the medium (or media)
ENUL = 1./REYNO
    GROUP 11. Initialization of variable or porosity fields
INIADD=F
IF(LCCM) THEN
+ FIINIT(WC1)= UIN; FIINIT(VC1)= 1.E-6
ELSE
+ FIINIT(W1) = UIN; FIINIT(V1) = 1.E-6
ENDIF
    GROUP 13. Boundary conditions and special sources
DO II=1,NY
+ YCUR=-1+(2*II-1)/NY
+ IF(LUNIF) THEN
+  WCR =UIN
+ ELSE
+  WCR =3./2.*(1.0-YCUR**2)
+ ENDIF
+ INLET(INL:II:,LOW,1,NX,II,II,1,1,1,LSTEP)
+ VALUE(INL:II:,P1,WCR*RHO1)
+ IF(LCCM) THEN
+  VALUE(INL:II:,VC1,0.0); VALUE(INL:II:,WC1,WCR)
+ ELSE
+  VALUE(INL:II:,V1,0.0); VALUE(INL:II:,W1,WCR)
+ ENDIF
ENDDO
    ** Walls.
PATCH(WN,NWALL,1,NX,NY,NY,1,NZ,1,LSTEP)
PATCH(WS,SWALL,1,NX, 1, 1,1,NZ,1,LSTEP)
IF(LCCM) THEN
+COVAL(WN,VC1,1.0,0.0); COVAL(WS,VC1,1.0,0.0)
+COVAL(WN,WC1,1.0,0.0); COVAL(WS,WC1,1.0,0.0)
ELSE
+COVAL(WN, W1,1.0,0.0); COVAL(WS, W1,1.0,0.0)
ENDIF
    ** Outlet.
PATCH(OUT,HIGH,1,NX,1,NY,NZ,NZ,1,LSTEP)
COVAL(OUT,P1,1000.0,0.0)
IF(LCCM) THEN
+ COVAL(OUT,WC1,ONLYMS,SAME); COVAL(OUT,VC1,ONLYMS,SAME)
ELSE
+ COVAL(OUT, W1,ONLYMS,SAME); COVAL(OUT, V1,ONLYMS,SAME)
ENDIF
    GROUP 15. Termination of sweeps
LSWEEP = 200; TSTSWP = -1
    GROUP 16. Termination of iterations
SELREF = T;   RESFAC = 1.E-3
    GROUP 17. Under-relaxation devices
RELAX(P1,LINRLX,0.25)
IF(.NOT.LCCM) THEN
+ RELAX(W1,FALSDT,0.5); RELAX(V1,FALSDT,0.5)
ENDIF
    GROUP 19. Data communicated by satellite to GROUND
IF(LCCM) THEN
    * LSG4= T activates nonorthogonality treatment in CCM;
    * LSG7= T permits CCM-solver to use higher order schemes.
+ LSG4= T;  LSG7= T
  SCHMBEGIN
      VARNAM VC1  SCHEME SUPERB
      VARNAM WC1  SCHEME SUPERB
  SCHMEND
ENDIF
    GROUP 22. Spot-value print-out
IXMON= 1;  IYMON = NY/2+1;  IZMON= NZ/2+1