TALK=T;RUN(1,1) inclined plate-flow:b103

  PHOTON USE
  p;;;;;;                                                                               
                                                                                
                                                                                
  vi 1 2 3
  gr x m                                                                          
  gr y m                                                                          
  gr z m 
  pause
  vi z
  gr z m
  dump gr
  pause
  gr off
  red
  set prop off
  con prps z 1 fi;.001  
  dump prps
  pause
  con off
  red
  vec z 1 sh
  dump vecs
  pause
  vec off
  red
  con u1 z 1 fi;.001                                                                       
  dump conu1
  pause
  con off
  red
  con p1 z 1 fi;.001 
  dump conp1                                                                      
  
  ENDUSE
TEXT(Parameterised inclined-plated flow :B103

  Display
  Air flows steadily between eight slightly-inclined steel plates, 
  which do not quite touch at their their ends.
  Body-fitted coordinates are used; and the turbulence model is 
  LVEL.
  Some problem-defining parameters can be set interactively, via 
  readvdu commands
  Enddis
  Groups 2 to 6. Geometry and grid
  Declarations
REAL(XLEN,XCO,XI,X0,XDX,XDY,XDZ,FNX,LPERCENT,DUMMY)
REAL(YLEN,YCO,YI,Y0,YDX,YDY,YDZ,FNY,YDISP)
REAL(ZLEN,ZCO,ZI,Z0,ZDX,ZDY,ZDZ,FNZ)
REAL(INVEL)
INTEGER(IXP,IYP,IZP,IX1,IX2,IY1,IY2,IFINE)
BOOLEAN(READGRID)

  grid-creation setting
READGRID=T  
  
  inlet velocity setting
INVEL=1.0
  
  geometry settings
XLEN=1.0; YLEN=0.1; ZLEN=1.0 ; YDISP=YLEN
LPERCENT=100
  grid fineness setting
IFINE=1
  
  Interactive setting of parameters
if(iqalib.ne.0) then
+ ans=y
endif
label ask
mesgm(Steady flow of air between slightly-inclined steel plates
              ----------- messages to  displau current options
mesg(length in flow direction is :xlen: m
mesg(thickness of air space plus plates is :ylen: m
mesg(displacement of plate ends normal to flow is :ydisp: m
mesg(length per cent is :lpercent:
mesg(inlet velocity of air is :invel:

if(ifine.eq.1) then
mesg(grid fineness is coarse
endif
if(ifine.eq.2) then
mesg(grid fineness is medium
endif
if(ifine.eq.3) then
mesg(grid fineness is finest
endif
mesg(readgrid option is :readgrid:
                ---------------------------- end of messages 
                ----------------- start of interactive input                     
mesgm(data OK ?(Y/n)
readvdu(ans,char,y)
if(:ans:.eq.y) then
goto continue
endif            ! all options will now be presented, in turn
mesg(insert new flow length?
readvdu(xlen,real,xlen)
mesg(insert new thickness?
readvdu(ylen,real,ylen)
mesg(insert new displacement?
readvdu(ydisp,real,ydisp)
mesg(insert new length percentage?
readvdu(lpercent,real,lpercent)
mesg(insert new inlet velocity?
readvdu(invel,real,invel)
mesg(change grid fineness? 1 coarse, 2 medium, 3 finest
readvdu(ifine,int,ifine)
mesg(change grid setting method? (y/N)
readvdu(ans,char,n)
if(:ans:.eq.y) then
readgrid=.not.:readgrid:
endif
goto ask
                ----------------- end of interactive input                     
label continue
                ----------------- start of execution                     
if(ifine.eq.1) then
NX=40; NY=20 ; NZ=1       ! coarse grid
endif
if(ifine.eq.2) then
NX=80;NY=40;NZ=1          ! medium grid
endif
if(ifine.eq.3) then
NX=160;NY=80;NZ=1         ! finest grid
endif
FNX=NX;FNY=NY;FNZ=NZ

BFC=T
if(readgrid) then
if(ifine.eq.1) then
readco(*/../d_earth/d_opt/d_bfc/inplib/103_1)
endif
if(ifine.eq.2) then
readco(*/../d_earth/d_opt/d_bfc/inplib/103_2)
endif
if(ifine.eq.3) then
readco(*/../d_earth/d_opt/d_bfc/inplib/103_3)
endif

else

mesg(creating bfc grid; Please wait  

x0 = 0.; y0 = 0.; z0 = 0.

xdx = xlen/fnx; ydy = ylen/fny; zdz = zlen/fnz

ydx = ydisp/fnx
  goto next
do ixx =0,nx
xi=ixx
if(ixx.lt.(nx/2+1)) then   ! upstream half
if(ixx.lt.(nx/4+1)) then   ! ist quarter
 y0 = y0 + ydx
else                       ! 2nd quarter
 y0=  y0 - ydx
endif
else
if(ixx.lt.(3*nx/4+1)) then ! downstream half
 y0 = y0 + ydx             ! third quarter
else
 y0=  y0 - ydx             ! fourth quarter
endif
endif
do izz =0,nz
zi=izz

do iyy =0,ny
yi=iyy

xco = x0 + xdx*xi     ! corner coordinates
yco = y0 + ydy*yi 
zco = z0 + zdz*zi
ixp=ixx+1
iyp=iyy+1
izp=izz+1
  xc(ixp,iyp,izp)=xco
  yc(ixp,iyp,izp)=yco
  zc(ixp,iyp,izp)=zco
setpt(ixp,iyp,izp,xco,yco,zco) ! equivalent to the above 3 settings

enddo
enddo
enddo

endif
label next
  Group 7
store(prps,enut,wdis)  
#solvel
turmod(lvel)
gcv=t
  Group 11
fiinit(u1)=2.0*invel
fiinit(prps)=air20


dummy=1.0 - lpercent*0.01   ! fraction cut off
dummy
dummy=dummy*fnx*0.25  ! proportion of cells cut off
dummy
ixp=dummy         ! number of cells cut off
ix1=1+ixp; ix2 = nx/4 -ixp
iy1=1; iy2=ny/4
patch(slab1,inival,ix1,ix2,iy1,iy2,1,nz,1,1)
coval(slab1,prps,0.0,steel)
iy1=3*ny/4+1; iy2=ny
patch(slab2,inival,ix1,ix2,iy1,iy2,1,nz,1,1)
coval(slab2,prps,0.0,steel)

ix1=nx/4+1+ixp; ix2 = nx/2 -ixp
iy1=1; iy2=ny/4
patch(slab3,inival,ix1,ix2,iy1,iy2,1,nz,1,1)
coval(slab3,prps,0.0,steel)
iy1=3*ny/4+1; iy2=ny
patch(slab4,inival,ix1,ix2,iy1,iy2,1,nz,1,1)
coval(slab4,prps,0.0,steel)

ix1=nx/2+1+ixp; ix2 = 3*nx/4-ixp
iy1=1; iy2=ny/4
patch(slab5,inival,ix1,ix2,iy1,iy2,1,nz,1,1)
coval(slab5,prps,0.0,steel)
iy1=3*ny/4+1; iy2=ny
patch(slab6,inival,ix1,ix2,iy1,iy2,1,nz,1,1)
coval(slab6,prps,0.0,steel)

ix1=3*nx/4+1+ixp; ix2 = nx -ixp
iy1=1; iy2=ny/4
patch(slab7,inival,ix1,ix2,iy1,iy2,1,nz,1,1)
coval(slab7,prps,0.0,steel)
iy1=3*ny/4+1; iy2=ny
patch(slab8,inival,ix1,ix2,iy1,iy2,1,nz,1,1)
coval(slab8,prps,0.0,steel)
  
  Group 13
patch(westin,west,1,1,1,ny,1,nz,1,1)
coval(westin,p1,fixflu,invel*1.189)  
coval(westin,u1,onlyms,invel)

patch(eastout,east,nx,nx,1,ny,1,nz,1,1)
coval(eastout,p1,fixp,0.0) 
if(lpercent.ge.100) then
patch(westin,west,1,1,1,ny/2,1,nz,1,1)
patch(eastout,east,nx,nx,1,ny/2,1,nz,1,1)
endif
#conprom
varmin(u1)=-1.e11;varmax(u1)=0.1*invel
varmin(v1)=-1.e11;varmax(v1)=0.1*invel
#maxmin
#endpause
tstswp=-1
lsweep=200 
nxprin=1;nyprin=1;ixprl=10;iyprl=10
                ----------------- end of execution                 
stop