CCM Test: Laminar/turbulent flow in a U-turned duct.
  **************************************************************
  DISPLAY
   This case concerns plane or axisymmetric (LCYL= T),  two-
   dimensional,  incompressible  laminar/turbulent (LTURB=T)
   flow through a  180 degree turnaround  duct. The size  of
   the inlet cross-section differs from that at the  outlet.
   The  flow  exhibits  large  streamline curvature together
   with flow separation near the bend exit next to the inner
   surface of the duct.
 
   Calculations of the turbulent flow could be made  whether
   with high-Re K-E models (standard (LKEMOD=T); RNG-derived
   (LRNG=T); 2-scale (LTWOSC=T)), or with low-Re  extensions
   (LOWRE=T) of K-E  models (standard; Chen-Kim  (LCHEK=T)).
   Note  that  low-Re  models  needs further grid refinement
   near walls, as compared to the grid set in the Q1-file.
 
   User can switch from the default colocated  computational
   algorithm (CCM)  to the  staggered one  (STAG) by setting
   LCCM= F.
   ---------------------------------------------------------
  ENDDIS
L(PAUSE
  **************************************************************
BOOLEAN(LCCM,LTURB,LCYL,LOWRE,LKEMOD,LRNG,LTWOSC,LCHEK)
LCCM = T;  LTURB = F;  LCYL= T
LKEMOD= T;  LRNG = F;  LTWOSC= F
LOWRE = F;  LCHEK= 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.001
   msg Press E  to exit PHOTON ...
  ENDUSE
  **************************************************************
    GROUP 1. Run title and other preliminaries
IF(LCCM) THEN
+IF(LTURB) THEN
+ IF(LOWRE) THEN
+  TEXT(CCM: U-duct (K-E + Lam-Brem).
+  IF(LCHEK) THEN
+   TEXT(CCM: U-duct (Chen-Kim K-E + Lam-Brem).
+  ENDIF
+ ELSE
+  TEXT(CCM: U-duct (K-E).
+  IF(LRNG) THEN
+   TEXT(CCM: U-duct (RNG K-E).
+  ENDIF
+  IF(LTWOSC) THEN
+   TEXT(CCM: U-duct (2-scale K-E).
+  ENDIF
+ ENDIF
+ELSE
+ TEXT(CCM: U-duct (Y-Z plane, Re=1000.).
+ENDIF
ELSE
+NONORT= T
+IF(LTURB) THEN
+ IF(LOWRE) THEN
+  TEXT(STAG: U-duct (K-E + Lam-Brem).
+  IF(LCHEK) THEN
+   TEXT(STAG: U-duct (Chen-Kim K-E + Lam-Brem).
+  ENDIF
+ ELSE
+  TEXT(STAG: U-duct (K-E).
+  IF(LRNG) THEN
+   TEXT(STAG: U-duct (RNG K-E).
+  ENDIF
+  IF(LTWOSC) THEN
+   TEXT(STAG: U-duct (2-scale K-E).
+  ENDIF
+ ENDIF
+ELSE
+ TEXT(STAG: U-duct (Y-Z plane, Re=1000.).
+ENDIF
ENDIF
TITLE
REAL(PI,REYNO,WIN,HP1,HP2,LPIP,RP1,RP2,DTHYD,TKEIN,EPSIN)
INTEGER(NZP1,NZC,NZP2,IC1)
PI = 3.1415
  ** Problem definition:
HP1 = 0.045;    HP2 = 0.0188;       LPIP = 0.1095
RP1 = 0.5*HP2;  RP2 = (HP1+2*RP1+HP2)/2
IF(LTURB) THEN
+ REYNO= 1.0E5
ELSE
+ REYNO= 1000.0
ENDIF
WIN = 26.25;  ENUL= WIN*HP1/REYNO
NX  = 1;      NY  = 12
NZP1= 8;      NZC = 16;  NZP2= 8;  NZ= NZP1+NZC+NZP2
    GROUP 6. Body-fitted coordinates or grid distortion
BFC = T; VUP = T; GSET(D,NX,NY,NZ,0.1,HP1,LPIP)
GSET(P,P1,0.0,          0.0, 0.0); GSET(P,P2,0.0,          0.0,LPIP)
GSET(P,P3,0.0,          HP1,LPIP); GSET(P,P4,0.0,          HP1, 0.0)
GSET(P,P5,0.0,HP1+2*RP1+HP2,LPIP); GSET(P,P6,0.0,HP1+2*RP1+HP2, 0.0)
GSET(P,P7,0.0,    HP1+2*RP1, 0.0); GSET(P,P8,0.0,    HP1+2*RP1,LPIP)
GSET(P,A2,0.0,          0.0,LPIP); GSET(P,A3,0.0,          HP1,LPIP)
GSET(P,A5,0.0,HP1+2*RP1+HP2,LPIP); GSET(P,A8,0.0,    HP1+2*RP1,LPIP)
GSET(L,L12,P1,P2,NZP1,1.0);        GSET(L,L23,P2,P3,NY,S1.5)
GSET(L,L34,P3,P4,NZP1,1.0);        GSET(L,L41,P4,P1,NY,S1.5)
GSET(F,F1,P1,-,P2,-,P3,-,P4,-);    GSET(M,F1,+K+J,1,1,1)
GSET(L,L25,A2,A5,NZC,-1.2,ARC,0.0,    RP2,LPIP+RP2)
GSET(L,L58,A5,A8,NY,S1.5)
GSET(L,L83,A8,A3,NZC, 1.0,ARC,0.0,HP1+RP1,LPIP+RP1)
GSET(L,L32,A3,A2,NY,S1.5)
GSET(F,F2,A2,-,A5,-,A8,-,A3,-);    GSET(M,F2,+K+J,1,1,NZP1+1)
GSET(L,L56,P5,P6,NZP2,1.0);        GSET(L,L67,P6,P7,NY,S1.5)
GSET(L,L78,P7,P8,NZP2,1.0);        GSET(L,L85,P8,P5,NY,S1.5)
GSET(F,F3,P5,-,P6,-,P7,-,P8,-); GSET(M,F3,+K+J,1,1,NZP1+NZC+1)
IF(LCYL) THEN
+ GSET(C,I:NX+1:,F,I1,1,NY,1,NZ,RZ,-0.1,0.0,0.0,INC,1.0)
ELSE
+ GSET(C,I:NX+1:,F,I1,1,NY,1,NZ, +, 0.1,0.0,0.0,INC,1.0)
ENDIF
GVIEW(X); VIEW
    GROUP 7. Variables stored, solved & named
SOLVE(P1,V1,W1)
IF(LTURB) THEN
+ IF(LTWOSC) THEN
+  TURMOD(TSKEMO)
+ ELSE
+  IF(LOWRE) THEN
+   IF(LCHEK) THEN
+    TURMOD(KECHEN-LOWRE)
+   ELSE
+    TURMOD(KEMODL-LOWRE)
+   ENDIF
+  ELSE
+   IF(LRNG) THEN
+    TURMOD(KERNG)
+   ELSE
+    TURMOD(KEMODL)
+   ENDIF
+  ENDIF
+ ENDIF
+ STORE(ENUT,LEN1)
ENDIF
IF(LCCM) THEN
L($F150)
ENDIF
    GROUP 11. Initialization of variable or porosity fields
INIADD = F
IF(LCCM) THEN
+ FIINIT(VC1) = 1.E-3; FIINIT(WC1) = 1.E-3
ELSE
+ FIINIT(V1)  = 1.E-3; FIINIT(W1)  = 1.E-3
ENDIF
IF(LTURB) THEN
+ TKEIN= (0.05*WIN)**2; EPSIN     = TKEIN**1.5*0.1643/0.09*HP1
+ FIINIT(P1)= 1.3E-4;   FIINIT(KE)= TKEIN;  FIINIT(EP)= EPSIN
ENDIF
    GROUP 13. Boundary conditions and special sources
    ** Inlet.
INLET(IN,LOW,1,NX,1,NY,1,1,1,LSTEP)
COVAL(IN,P1,FIXFLU,RHO1*WIN)
IF(LCCM) THEN
+ COVAL(IN,VC1,ONLYMS,0.0); COVAL(IN,WC1,ONLYMS,WIN)
ELSE
+ COVAL(IN,V1, ONLYMS,0.0); COVAL(IN,W1, ONLYMS,WIN)
ENDIF
    ** Outlet.
PATCH(OUT,HIGH,1,NX,1,NY,NZ,NZ,1,LSTEP); COVAL(OUT,P1,FIXP,0.0)
    ** Walls.
PATCH(WS1,SWALL,1,NX,1, 1,      1,NZP1,1,LSTEP)
PATCH(WS2,SWALL,1,NX,1, 1, NZP1+1,  NZ,1,LSTEP)
PATCH( WN,NWALL,1,NX,NY,NY,     1,  NZ,1,LSTEP)
IF(LTURB) THEN
+ COVAL(IN, KE, ONLYMS,TKEIN); COVAL(IN, EP,ONLYMS,EPSIN)
+ IF(LCCM) THEN
+  COVAL(WS2,VC1,LOGLAW,0.0); COVAL(WN,VC1,LOGLAW,0.0)
+  COVAL(WS2,WC1,LOGLAW,0.0); COVAL(WN,WC1,LOGLAW,0.0)
+ ELSE
+  COVAL(WS2,W1, LOGLAW,0.0); COVAL(WN, W1,LOGLAW,0.0)
+ ENDIF
+ IF(LOWRE) THEN
+  COVAL(WS2,  KE,1.0,0.0);  COVAL(WN,  KE,1.0,0.0)
+  COVAL(WS2,LTLS,1.0,0.0);  COVAL(WN,LTLS,1.0,0.0)
+  IF(.NOT.LCYL) THEN
+   COVAL(WS1,LTLS,1.0,0.0); COVAL(WS1, KE,1.0,0.0)
+  ENDIF
+ ELSE
+  IF(LTWOSC) THEN
+   COVAL(WS2,KP,LOGLAW,LOGLAW); COVAL(WN,KP,LOGLAW,LOGLAW)
+   COVAL(WS2,ET,LOGLAW,LOGLAW); COVAL(WN,ET,LOGLAW,LOGLAW)
+   COVAL(WS2,KT,LOGLAW,LOGLAW); COVAL(WN,KT,LOGLAW,LOGLAW)
+   IF(.NOT.LCYL) THEN
+    COVAL(WS1,KP,LOGLAW,LOGLAW);COVAL(WS1,ET,LOGLAW,LOGLAW)
+    COVAL(WS1,KT,LOGLAW,LOGLAW)
+   ENDIF
+  ELSE
+   COVAL(WS2,KE,LOGLAW,LOGLAW); COVAL(WN,KE,LOGLAW,LOGLAW)
+   IF(.NOT.LCYL) THEN
+    COVAL(WS1,KE,LOGLAW,LOGLAW)
+   ENDIF
+  ENDIF
+  COVAL(WS2,EP,LOGLAW,LOGLAW); COVAL(WN,EP,LOGLAW,LOGLAW)
+  IF(.NOT.LCYL) THEN
+   COVAL(WS1,EP,LOGLAW,LOGLAW)
+  ENDIF
+ ENDIF
+ COVAL(WS2,C3,LOGLAW, 0.0);   COVAL(WN, C3, LOGLAW, 0.0)
+ COVAL(OUT,KE,ONLYMS,0.0);   COVAL(OUT,EP, ONLYMS,0.0)
+ IF(.NOT.LCYL) THEN
+  COVAL(WS1,C3,LOGLAW, 0.0)
+  IF(LCCM) THEN
+   COVAL(WS1,VC1,LOGLAW, 0.0); COVAL(WS1,WC1,LOGLAW, 0.0)
+  ELSE
+   COVAL(WS1,W1, LOGLAW, 0.0)
+  ENDIF
+ ENDIF
ELSE
+ IF(LCCM) THEN
+  COVAL(WS2,VC1,1.0,0.0);   COVAL(WN, VC1,1.0,0.0)
+  COVAL(WS2,WC1,1.0,0.0);   COVAL(WN, WC1,1.0,0.0)
+ ELSE
+  COVAL(WS2,W1, 1.0,0.0);   COVAL(WN, W1, 1.0,0.0)
+ ENDIF
+ IF(.NOT.LCYL) THEN
+  IF(LCCM) THEN
+   COVAL(WS1,VC1,1.0,0.0);  COVAL(WS1,WC1,1.0,0.0)
+  ELSE
+   COVAL(WS1,W1, 1.0,0.0)
+  ENDIF
+ ENDIF
ENDIF
    GROUP 15. Termination of sweeps
LSWEEP = 500; TSTSWP = -1
    GROUP 16. Termination of iterations
SELREF = T;   RESFAC = 1.E-3
    GROUP 17. Under-relaxation devices
DTHYD = 10.*LPIP/NZP1/WIN; RELAX(P1,LINRLX,0.5)
IF(LTURB) THEN
+ KELIN= 1;  RELAX(EP, FALSDT,DTHYD)
+ IF(LTWOSC) THEN
+  RELAX(KP,FALSDT,DTHYD);RELAX(ET,FALSDT,DTHYD)
+  RELAX(KT,FALSDT,DTHYD)
+ ELSE
+  RELAX(KE,FALSDT,DTHYD)
+ ENDIF
ENDIF
IF(LCCM) THEN
+ RELAX(VC1,FALSDT,DTHYD);RELAX(WC1,FALSDT,DTHYD)
ENDIF
    GROUP 19. Data communicated by satellite to GROUND
IF(LCCM) THEN
    * LSG4 activates non-orthogonality treatment in CCM/MBFGE.
+ LSG4= T
ENDIF
    GROUP 22. Spot-value print-out
IXMON= 1;  IYMON= NY/2+1;  IZMON= NZ/2+1