Test: Laminar flow in a 2D-duct with gradual expansion.
  *******************************************************
  DISPLAY
    The problem considered is that defined by the International
  Association for Hydraulic Research Working Group on Refined
  Modelling of Flows for their 5th meeting. The report of the
  workshop is presented as a paper entitled 'Laminar flow in a
  complex geometry: a comparison', Int. j. numer. methods fluids,
  vol.5, 667-683(1985).
    The main features of the problem are:
  * laminar flow, Re=10.0 ( and 100.0);
  * the north and south boundaries of the channel are given by the
    following analytical expression :
           y(z) = +/- [tanh(2-30*z/Re)-tanh(2)]/2
           for  0 < z< Re/3
  * a parabolic profile for the axial velocity is specified at the
    inlet,
           w=3*(y-y*y/2) for z=0 and -1 < y < 1,
           v=0.
  ENDDIS
  PHOTON USE
   p ; ; ; ; ;
 
   msg Computational Domain:
   gr k 1
   pause
   cl
   set vec av off
   msg Velocity Vectors:
   vec k 1 sh
   msg Press Any Key to Continue...
   pause
   cl
   msg Contours of UCRT-component:
   con ucrt k 1 fi;0.01
   msg Press Any Key to Continue...
   pause
   cl
   msg Contours of p1
   con p1 k 1 fi;0.01
   msg Press E  to exit PHOTON ...
  ENDUSE
BOOLEAN(LUNIF); LUNIF= F
  **************************************************************
    GROUP 1. Run title and other preliminaries
INTEGER(NY1,NY2)
REAL(REYNO,UIN,LX,LY,DX,DY,DX1,DX2)
REAL(TAN2,XCUR,YCUR,UCR,TARG,YWL)
  ** Problem definition:
MESG(Select Reynolds Number:
MESG(   A - Re = 10
MESG(   B - Re = 100
MESG(Enter choice (A - Re=10)
READVDU(ANS,CHAR,A)
REYNO=10
IF(:ANS:.EQ.2) THEN
REYNO=100
ENDIF
UIN= 1.0
NY1  = 32;    NY2= NY1; NX  = 32;       NY  = NY1+NY2
LX   = REYNO/3.0;      DX  = LX/NX;     DX1 = DX*0.3; DX2= DX*2
XCUR = LX ;            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,LX,LY,1.0)
NONORT=T
GSET(P,P1,   0.0,-1.0,0.0);  GSET(P,P2,   0.0,1.0,0.0)
GSET(P,P3,   DX1,-1.0,0.0);  GSET(P,P4,   DX1,1.0,0.0)
GSET(P,P5,LX-DX2, -LY,0.0);  GSET(P,P6,LX-DX2, LY,0.0)
GSET(P,P7,    LX, -LY,0.0);  GSET(P,P8,    LX, LY,0.0)
GSET(V,YWT,S,P4,SPLINE)
DO II =1,5
+ XCUR= (LX/2-DX1)/5*II + DX1; TARG= 2.0-10.*XCUR/LX
+ YCUR= 1.0-(TANH(TARG)-TAN2)/2.0
+ GSET(V,XCUR,YCUR,0.0)
ENDDO
GSET(V,YWT,E,P6,1.0,0.0,0.0,0.0)
GSET(L,LWT,P4,P6,NX-2,1.2CRVYWT)
GSET(V,YWB,S,P3,SPLINE)
DO II =1,5
+ XCUR= (LX/2-DX)/5*II + DX1; TARG= 2.0-10.*XCUR/LX
+ YCUR= -(1.0-(TANH(TARG)-TAN2)/2.0)
+ GSET(V,XCUR,YCUR,0.)
ENDDO
GSET(V,YWB,E,P5,1.0,0.0,0.0,0.0)
GSET(L,LWB,P3,P5,NX-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+I, 1,1,1)
GSET(F,F2,P3,-,P4,-,P6,-,P5,-); GSET(M,F2,+J+I, 2,1,1,LAP1)
GSET(F,F3,P5,-,P6,-,P8,-,P7,-); GSET(M,F3,+J+I,NX,1,1)
GSET(C,K:NZ+1:,F,K1,1,NX,1,NY,+,0.0,0.0,1.0,INC,1.0)
MESG(Show grid? (N/y)
READVDU(ANS,CHAR,N)
IF(:ANS:.EQ.Y) THEN
VIEW(K,1)
ENDIF
    GROUP 7. Variables stored, solved & named
SOLVE(P1,U1,V1)
    GROUP 9. Properties of the medium (or media)
ENUL = 1./REYNO
    GROUP 11. Initialization of variable or porosity fields
INIADD=F
FIINIT(U1)  = UIN; FIINIT(V1)  = 1.E-3
    GROUP 12. Unused
    GROUP 13. Boundary conditions and special sources
DO II=1,NY
+ YCUR=-1+(2*II-1)/NY
+ IF(LUNIF) THEN
+  UCR =UIN
+ ELSE
+  UCR =3./2.*(1.0-YCUR**2)
+ ENDIF
+ INLET(BFCINL:II:,WEST,1,1,II,II,1,NZ,1,LSTEP)
+ VALUE(BFCINL:II:,P1,GRND1)
+ VALUE(BFCINL:II:,U1,GRND1)
+ VALUE(BFCINL:II:,V1,GRND1)
+ VALUE(BFCINL:II:,UCRT,UCR)
ENDDO
BFCA=RHO1
    ** Walls.
PATCH(WN,NWALL,1,NX,NY,NY,1,NZ,1,LSTEP)
PATCH(WS,SWALL,1,NX, 1, 1,1,NZ,1,LSTEP)
COVAL(WN, U1,1.0,0.0); COVAL(WS, U1,1.0,0.0)
    ** Outlet.
PATCH(OUT,EAST,NX,NX,1,NY,1,NZ,1,LSTEP)
COVAL(OUT,P1,1000.0,0.0)
COVAL(OUT, U1,ONLYMS,0.0); COVAL(OUT, V1,ONLYMS,0.0)
    GROUP 15. Termination of sweeps
LSWEEP = 200; TSTSWP = -1
    GROUP 16. Termination of iterations
SELREF = T;   RESFAC = 1.E-3
    GROUP 19. Data communicated by satellite to GROUND
MESG(Use GCV solver? (Y/n)
READVDU(ANS,CHAR,Y)
GCV=:ANS:.EQ.Y
IF(GCV) THEN
 LSG9=T; LSG8=T
 TEXT(IAHR Smooth Channel, GCV, Re = :REYNO:
ELSE
 SYMBFC=T
 RELAX(P1,LINRLX,0.6)
 RELAX(U1,FALSDT,0.5); RELAX(V1,FALSDT,0.5)
 LSWEEP=400
 TEXT(IAHR Smooth Channel, Re = :REYNO:
ENDIF
TITLE
    GROUP 22. Spot-value print-out
IXMON= NX/2+1;  IYMON= NY/2+1;  IZMON= 1
    GROUP 24. Dumps for restarts