Author's explanation 01.04.09
  -----------------------------
  This case is a development of case 290 which concerns turbulent
  flow over and downstream of a backward-facing step but uses a grid 
  which is far too coarse for accuracy.
  
  It has therefore been modified so as to enable the effects of the  
  grid-defining parameters NX, NY, XPOWER and YPOWER to be easily varied.
  This is done by changing the name of the STEP Patch to .STEP . 
  
  The '.' which starts its name is a signal to EARTH that locating 
  arguments are geometric rather than indicial, albeit non-
  dimensionalised by reference to the domain dimensions.
  A multiplier of 1000 is used (and matched by division on arrival 
  in EARTH) because, although the arguments are real numbers, the
  SATELLITE reduces them to whole numbers.
  
  Although it is not necessary, because they lie at domain 
  boundaries, the inlet, outlet and wall boundary conditions are
  also expressed via dot-patches so as to test and illustrate the 
  technique.
  
  DBS
  -----------------------------------------------------------------
  Setting the grid fineness and distribution interactively
NX=20; ny=15     ! defaults
real(xpower,ypower)
xpower=1.0; ypower=1.0
if(iqalib.ne.0) then
+ ans=y
endif
label message
mesgm(nx=:nx: ny=:ny: xpower=:xpower: ypower=:ypower: If ok, RETURN
readvdu(ans,char,y)
if(:ans:.eq.y) then
 goto begin
else
 mesg(insert required NX, or RETURN
 readvdu(nx,int,nx)  
 mesg(insert required NY, or RETURN 
 readvdu(ny,int,ny)  
 mesg(insert required XPOWER, or RETURN 
 readvdu(xpower,real,xpower)  
 mesg(insert required YPOWER, or RETURN 
 readvdu(ypower,real,ypower)  
 goto message
endif 
label begin
  
  PHOTON USE
  p;;;;
 
  use patgeo
  gr ou z 1
  msg contours of velocity, and velocity vectors
  con u1 z 1 fi;0.01;vec z 1
  msg press RETURN for length scale
  pause;con off;vec off;red
  msg contours of turbulence length scale
  con len1 z 1 fi;0.01;pause;con off;red
  msg contours of effective viscosity
  con enut z 1 fi;0.01
  enduse
    GROUP 1. Run title and other preliminaries
TEXT(Backward-Facing Step; K-E or LVEL 
TITLE
  DISPLAY
  This case involves steady, turbulent flow over a backward-facing
  step, of height h, in a two-dimensional channel of width 3*h.
  The calculations are started at a distance 4h upstream of the step
  and terminated at a distance 16h downstream.
 
    //////////////////////// wall ///////////////////////
    -----------------------------------------------------
                                                           Pressure
  Inlet  ------->                              --------->  fixed at
                                                           zero
   /________________
  | ////////////////|                                  Exit
  |       wall     /| Recirculation
  |                /|  <----                ----->
   \ constant      /|____________________________________
  prescribed mass  //////////////////////////////////////
  inflow rate                        wall
  ENDDIS
REAL(HEIGHT,WIDTH,CLEN,SLEN,UIN,TKEIN,EPSIN,GENUT)
UIN=1.0
     ** Calculation of domain specifications
HEIGHT=0.0381;
WIDTH=3.*HEIGHT
SLEN=4.*HEIGHT;
CLEN=20.*HEIGHT

XULAST=CLEN
YVLAST=WIDTH
GRDPWR(X,NX,XULAST,XPOWER)
GRDPWR(Y,NY,YVLAST,YPOWER)
    GROUP 7. Variables stored, solved & named
SOLVE(P1,U1,V1);SOLUTN(P1,Y,Y,Y,N,N,N);STORE(ENUT,LEN1)
mesg(Turbulence model is KEMODL. Use LVEL instead? (y/n)
CHAR(char1)
READVDU(ANS,CHAR,N)
IF(:ANS:.EQ.Y) THEN
 TURMOD(LVEL)
 char1=LVEL
 store(wgap)
ELSE
 TURMOD(KEMODL)
 char1=KE-EPS
 KELIN=3
 mesgb(kelin=3; no under-relaxation is used for ke or ep0
ENDIF
    GROUP 9. Properties of the medium (or media)
    GROUP 11. Initialization of variable or porosity fields
     ** Calculation of KE (where fric=0.018)...
TKEIN=0.25*UIN*UIN*0.018
     ** Calculation of EP (where lmix=0.09 x h)...
EPSIN=TKEIN**1.5*0.1643/(0.09*height)
     ** Initial values
FIINIT(U1)=UIN;FIINIT(P1)=1.3E-4;FIINIT(KE)=TKEIN;FIINIT(EP)=EPSIN
     ** Initialization of variables in blocked region
 EGWF=T;STORE(PRPS)
   declarations and settings of dot-patch arguments
 REAL(RXF,RXL,RYF,RYL,RZF,RZL)
 RXF=0.0; RXL=SLEN/XULAST*1000
 RYF=0.0; RYL=HEIGHT/YVLAST*1000
 RZF=0.0; RZL=ZWLAST/ZWLAST*1000
 PATCH(.STEP,INIVAL,RXF,RXL,RYF,RYL,RZF,RZL,1,1) ! Important warning
    ! Although it might seem reasonable to do so, DO NOT set the 
    ! last but one arguments above to 0 because doing so confuses the 
    ! indexing system of the satellite, with unforeseeable results.
    ! Note also that 1,1000 is satisfactory for the last two 
    ! arguments of the group 13 dot-patches; but the equally 
    ! reasonable 0,1000, although OK for the x, y and z arguments,
    ! causes the aforesaid  damage.
 COVAL(.STEP,PRPS,0.0,198.0)
TEXT( :char1: model;  EGWF= :EGWF:
TITLE
   GROUP 13. Boundary conditions and special sources
  **  Inlet
  PATCH(INLET,WEST,1,1,1,NY,1,1,1,1) ! the no-dot patch
PATCH(.INLET,WEST,0,0,0,1000,0,1000,1,1000)
COVAL(.INLET,P1,FIXFLU,RHO1*UIN)
COVAL(.INLET,U1,ONLYMS,UIN)
COVAL(.INLET,KE,ONLYMS,TKEIN)
COVAL(.INLET,EP,ONLYMS,EPSIN)
  **  Exit
  PATCH(OUTLET,EAST,NX,NX,1,NY,1,1,1,1) ! the no-dot patch
PATCH(.OUTLET,EAST,1000,1000,0,1000,0,1000,1,1000)
COVAL(.OUTLET,P1,1.0E5,0.0)
  **    N-wall
  WALL (WFNN,NORTH,1,NX,NY,NY,1,1,1,1)  ! the no-dot patch
WALL (.WFNN,NORTH,0,1000,1000,1000,0,1000,1,1000)
  **    S-wall
WALL (.WFNS,SOUTH,0,1000,0,0,0,1000,1,1000)
    GROUP 15. Termination of sweeps
LSWEEP=1000;SELREF=T;RESFAC=0.01
    GROUP 16. Termination of iterations
LITER(U1)=20;LITER(V1)=20
    GROUP 17. Under-relaxation devices
RELAX(U1,FALSDT,SLEN/UIN);RELAX(V1,FALSDT,SLEN/UIN)
relax(ke,linrlx,0.5);relax(ep,linrlx,0.5)
    GROUP 22. Monitor print-out
IYMON=NY/2;IXMON=NX/2;NPRMON=100;TSTSWP=12;UWATCH=T
    GROUP 23. Field print-out and plot control
NPLT=1;IPLTL=LSWEEP;ITABL=1
PATCH(MAP,CONTUR,1,NX,1,NY,1,1,1,1)
PLOT(MAP,P1,0.0,10.0);PLOT(MAP,U1,0.0,10.0);PLOT(MAP,V1,0.0,10.0)
tstswp=-1
  lsweep=1
  debug=t
  dbindx=t
  search=t
  idbf0=414
  flag=t
stop