TALK=T;RUN( 1, 1)
************************************************************
Group 1
************************************************************
save1begin
Notes of experiences, by dbs:
1. The inadvertent erasure of the following line
> DOM, SIZE, :domlong:,:domwide: , :domhigh:
caused the floor object to be ignored, and unexpected objects
to appear.
I recommend that satellte should print a clear warning message
reporting why it is ignoring objects which are represented
correctly in the Q1 file.
2. The pipe VR object, de-activated below led to convergence
difficulties, even without the presence of the wall and the
release. Its size was such that, with the grid prescribed, it
created y-direction lines of four adjacent cut cells.
It has therefore been replaced by two in-form LINE object, to
leave space for the release object between their ends.
3. However, it emerged that the line object, at least as
implemented, did not stretch across the domain. An ellipsoid
object will be tried instead.
4. An old 'quick fix' in e2eaio.for switched off conwiz (i.e.
lsg57) when in-form objects were present. I have deactivated
this, because conwiz-related features, such as the control of
dvdp's, are need to procude convergence.
Without them very large (~1.e9) dvdps were found for v1 and w1
near the wall.
Whether the VR-pipe convergence difficulties would disappear
with conwiz=t has not been investigated.
DISPLAY
This parameterised relational Q1 was created by dbs in Nov. 2008
It represents the release of gas from a hole in a pipeline.
Of interest is what influence a nearby wall has on the spread of
gas into the atmosphere.
__________________ north boundary ________________
� � �
� � �
�< indist >< walldist >�< outdist >�
� � �
� * � �
� � �
�_________________ south boundary _______________�_
west (inlet) release wall east (outlet)
------)--------
- ) wind angle
-
-
The effect of non-zero wind angle can be represented in two ways,
namely:
(1) grid turns with wind, and
(2) flow occurs through north, south, east and west boundaries.
In the latter case,in order to avoid switching from inlet to
outlet type, but also for physical realism, it is best to fix
velocities rather than mass flows, and to make the top boundary
an open constant-pressure one. This may just as well be done for
method (1) also.
However, Method (2) is the only one represented below.
ENDDIS
!-------------------------------- declarations
! primary parameters, real
real(wallhigh,wallthck,walldist,indist,outdist)
real(domhigh,domwide)
real(pipezpos,relhigh,floorco)
real(pipediam,reldiam,relangl,gaspr)
real(windvel,windangl) ! wind parameters
! secondary parameters
real(domlong,molwt,scale)
real(uwind,vwind) ! horizontal wind components
real(xlen,xtot)
real(xp,yp,zp,xs,ys,zs,al,be,ga)
real(xpw,ypw,zpw,xsw,ysw,zsw) ! wall parameters for use in BOUND
real(rad)
real(valu)
! primary parameters, integer
integer(nxgrid,nygrid,nzgrid)
! secondary parameters
integer(nrx,ixwal1,ixwal2)
integer(intger)
integer(secno)
integer(nrz)
integer(io) ! index of in-form objects
integer(iii)
! primary parameters, character
char(gas) ! what gas is released
char(tmodel,profl) ! turbulence model; wind profile
! secondary parameters
char(uprofl,vprofl,gasden,gasvel,gasflo)
char(item,ch)
! primary parameters, boolean
boolean(pipe,wal,releas,diagnos,uniform,restart)
! secondary parameters
boolean(bool,subwal)
!----------------------------------------- settings
! reals
! geometrical settings
wallhigh=2.5 ! used as a general length-scale parameter
scale=wallhigh ! by means of this setting of scale
wallthck=0.2*scale ! which therefore influence all the rest
domhigh=3.5*scale
domwide= 12*scale
pipezpos=0.4*scale
indist= 2.5*scale
walldist=7.5*scale
outdist =14*scale
reldiam=2.0*0.0254 ! release diameter specified in inches
relhigh= 0.4*scale
relangl= 90 *3.14159/180 ! release angle; 0 = +ve x direction
relangl= 0. *3.14159/180 ! release angle; 0 = +ve x direction
! 90 = +ve z direction
floorco=0.01
floorco=0.005
pipediam=0.16*scale
windvel = 2
windangl=-45 ! wind from south-west
windangl=45 ! wind from north-west
windangl= 0 ! wind from west
windangl=180 ! wind from east
windangl=-90 ! wind from south
windangl= 90 ! wind from north
! integers, defining the grid
nxgrid=50
nxgrid=5
nygrid=30
nygrid=1
nzgrid=35
nzgrid=5
! characters
gas=h2s ! hydrogen sulphide
gas=propane ! other possibilities are methane, ethane and propane
gaspr = 1.e5 ! i.e. one atmosphere
gaspr = .0
gaspr
profl=(zg/:zwlast:)^(1./7.) ! one-seventh power law
profl=1.0 ! uniform velocity
tmodel=lvel
tmodel=ke
! booleans
! these logicals enable features of the simulation to be
! switched on or off, which permits their individual
! influences to be studied
releas=t
releas=f
wal=t
wal=f
pipe=t
pipe=f
diagnos=f
uniform=f
restart=f
restart=t
! ----------------------- end of settings -----------------
! ------- Interactive inputs Part 1 Information ------------
label input ! print to screen and satlog.txt
mesgm(
mesga(Summary of parameters which can be changed interactively
mesg(1. geometry settings in meters:
mesg(wall height=:wallhigh: ;wall thick=:wallthck:
mesg(domwide=:domwide:; domhigh=:domhigh:
mesg(indist=:indist:, walldist=:walldist:;
domlong =indist+walldist+outdist
mesg(dist. from wall to east=:outdist:; so domlong=:domlong:
mesga(2. wind parameters
mesg(windvel=:windvel:; angle=:windangl:; profile=:profl:
mesga(3. logicals
mesg(releas=:releas:;restart=:restart:; wall =:wal:; pipe=:pipe:
mesg(uniform=:uniform:; diagnos=:diagnos:
mesga(4. release parameters
mesg(gas=:gas:; excess pressure in atm,=:gaspr:
mesg(relhigh=:relhigh:; relangl=:relangl:; reldiam=:reldiam:
mesg(pipediam=:pipediam:; pipezpos=:pipezpos:
mesga(5. grid parameters
mesg(nxgrid=:nxgrid:; nygrid=:nygrid:; nzgrid=:nzgrid:
! -------- Interactive inputs Part 2 Question and answer --------
mesgm(Are above settings OK?. If not enter section number to change
readvdu(secno,int,0)
if(secno.eq.0) then ! section number 0
goto continue
else
if(secno.eq.1) then ! section number 1
! Note that the messages in these sections are copies of those
! which have appeared above. The repetition could be avoided by
! use of the PIL 'array' feature.
mesg(1. geometry settings in meters:
mesg(wall height=:wallhigh: ;wall thick=:wallthck:
mesg(domwide=:domwide:; domhigh=:domhigh:
mesg(dist. of release from west, to wall=:indist:, walldist
mesg(dist. from wall to east=:outdist:; so domlong=:domlong:
valu=wallhigh ;item='wallhigh'
CALL CHREAL
valu=wallthck; item='wallthck'
CALL CHREAL
valu=domwide ; item='domwide'
CALL CHREAL
valu=domhigh ; item='domhigh'
CALL CHREAL
valu=indist; item='indist'
CALL CHREAL
valu=walldist; item='walldist'
CALL CHREAL
valu=outdist ; item='outdist'
CALL CHREAL
endif
if(secno.eq.2) then ! section number 2
mesga(2. wind parameters
mesg(speed=:windvel:; angle=:windangl:; profile=:profl:
valu=windvel ; item='windvel'
CALL CHREAL
valu=windangl; item='windangl'
CALL CHREAL
ch=profl ;item='profl'
CALL CHCHAR
endif
if(secno.eq.3) then ! section number 3
mesga(3. logicals
mesg(releas=:releas:; restart=:restart:; wall =:wal:; pipe=:pipe:
mesg(uniform=:uniform:; dignos=:diagnos
bool=restart; item='restart'
CALL CHBOOL
bool=wal; item='wal'
CALL CHBOOL
bool=pipe; item='pipe'
CALL CHBOOL
bool=uniform; item='uniform'
CALL CHBOOL
bool=diagnos; item='diagnos'
CALL CHBOOL
endif
if(secno.eq.4) then ! section number 4
mesga(4. release parameters
mesg(gas=:gas:; excess pressure in atm,=:gaspr:
mesg(relhigh=:relhigh:; relangl=:relangl:; reldiam=:reldiam:
mesg(pipediam=:pipediam:; pipezpos=:pipezpos:
ch=:gas: ; item='gas'
CALL CHCHAR
valu=relhigh ; item='relhigh'
CALL CHREAL
valu=relangl ; item='relangl'
CALL CHREAL
valu=reldiam ; item='reldiam'
CALL CHREAL
valu=pipediam ; item='pipediam'
CALL CHREAL
valu=pipezpos ; item='pipezpos'
CALL CHREAL
endif
if(secno.eq.5) then ! section number 5
endif
mesga(5. grid parameters
mesg(nxgrid=:nxgrid:; nygrid=:nygrid:; nzgrid=:nzgrid:
mesg(uniform=:uniform:
intger=nxgrid; item='nxgrid'
CALL CHINT
mesg(uniform=:uniform:
intger=nygrid; item='nygrid'
CALL CHINT
mesg(uniform=:uniform:
intger=nzgrid; item='nzgrid'
CALL CHINT
endif
label continue ! end of interactive sections
! --------------- Interactive inputs End --------
!-------------- consequential settings -------------------
mesgm(consequences of settings
!----------------------- consequences GROUP 1 ------------
molwt=28.9 ! i.e. that of air, the default gas
if(:gas:.eq.h2s) then
molwt=34.08
endif
if(:gas:.eq.methane) then
molwt=16
endif
if(:gas:.eq.ethane) then
molwt=30
endif
if(:gas:.eq.propane) then
molwt=44
endif
mesg(gas=:gas:; molwt=:molwt:
domlong =indist+walldist+outdist ! the implied domain length
windangl=windangl*3.14159/180 ! convert angle to radians
uwind=windvel*cos(windangl) ! deduce x- & y-direction
vwind=-windvel*sin(windangl) ! wind components
uprofl=:uwind:*:profl:
vprofl=:vwind:*:profl:
if(releas) then
gasden=1.189*molwt/28.9 ! air density * molecular-weight ratio
gasvel = (2*:gaspr:/:gasden:)^.5 ! an approximation
gasflo = :gasvel:*:gasden:*:relDIAM:^2*3.14159/4 ! vel*den*area
else
gasvel=0.0
gasflo=0.0
endif
gasvel ! print to satlog for checking
gasflo
TEXT(gas-release with wind angle = :windangl:
libref=162
save1end
!---------------- consequences for further groups ---------
************************************************************
Group 2. Transience
save2begin
STEADY = T
save2end
************************************************************
********************* non-uniform-grid settings **********
Group 3
save3begin
(a name="unigrid">
if(uniform) then ! activate the following if uniform grid is ! chosen
nx=nxgrid
ny=nygrid
nz=nzgrid
xulast=domlong
yvlast=domwide
zwlast=domhigh
#unigrid
mesg(go to next
goto next ! skip the next settings if grid is uniform
endif
(a name="finegrid">
********************* x grid *****************************
ixwal1=0
mesgm(grid information
nregx=5 ! 5 regions, of length 1. indist-pipediam
2. 6.*pipediam
3. walldist- 5.pipediam-wallthck
4. 3.*wallthck
5. outdist-2.&wallthck
xlen=indist-pipediam
xtot=xlen
nrx=nxgrid/10
ixwal1=ixwal1+nrx
iregx=1;grdpwr(x,-nrx,-xlen,-1.15) ! from west to release
mesg(west to release; iregx=:iregx: nrx=:nrx: xlen=:xlen:
xlen=6.*pipediam
xtot=xtot+xlen
nrx=nxgrid/5
ixwal1=ixwal1+nrx
iregx=2;grdpwr(x,nrx,-xlen,1.0) ! vicinity of release
mesg(around release; iregx=:iregx: nrx=:nrx: xlen=:xlen:
xlen=walldist-wallthck-5.*pipediam
xtot=xtot+xlen
nrx=nxgrid/5
ixwal1=ixwal1+nrx
iregx=3;grdpwr(x,-nrx,-xlen,1.15) ! from release to wall
mesg(release to wall; iregx=:iregx: nrx=:nrx: xlen=:xlen:
xlen=6.0*wallthck
xtot=xtot+xlen
nrx=nxgrid/10
ixwal2=ixwal1+nrx
iregx=4;grdpwr(x,nrx,-xlen,1.0) ! vicinity of wall
mesg(around wall; iregx=:iregx: nrx=:nrx: xlen=:xlen:
xlen=outdist-4.*wallthck
xtot=xtot+xlen
mesg(total length is :xtot:)
nrx=2*nxgrid/5
iregx=5;grdpwr(x,nrx,-xlen,1.05) ! from wall to east
mesg(wall to east; iregx=:iregx: nrx=:nrx: xlen=:xlen:
Group 4
*************************** y grid *********************
char(form)
form=-nygrid,-domwide,0.9
form=-nygrid,-domwide,1.0
grdpwr(y,:form:)
Group 5
*************************** z grid *********************
nregz=2
xlen=domhigh/2
nrz=5*nzgrid/7
iregz=1;grdpwr(z,nrz,-xlen,1.0)
nrz=2*nzgrid/7
iregz=2;grdpwr(z,nrz,-xlen,1.1)
nogrid=t ! setting needed to tell the Satellite not to modify
! the just-specified grid in any way.
label next ! end of grid settings
view(k,1) ! display grid (not yet active)
save3end
************************************************************
Group 7. Variables: STOREd,SOLVEd,NAMEd
************************************************************
Echo InForm settings for Group 7
save7begin
store (p1,u1,v1,w1,prps,vabs)
solutn (u1, y, y, y, n, n, y )
varmax(u1) = 1.0 ;varmin(u1) =-1.0e+11
output(u1,y,n,y,y,y,y)
solutn (v1, y, y, y, n, n, y )
varmax(v1) = 1.0 ;varmin(v1) =-1.0e+11
output(v1,y,n,y,y,y,y)
solutn (w1, y, y, y, n, n, y )
varmax(w1) = 1.0 ;varmin(w1) =-1.0e+11
output(w1,y,n,y,y,y,y)
solutn (p1, y, y, y, n, n, y )
output(p1,y,n,y,y,y,y)
if(releas) then
store (gas)
store (g.25)
solutn (gas, y, y, y, n, n, y )
output(gas,y,n,y,y,y,y)
varmin(gas)=0.0
(stored var g.25 is gas^0.25) ! g.25, i.e. the fourth root of gas,
! is computed because it is more
! visible than gas when contours are
! plotted in the viewer or photon
output(g.25,y,n,y,y,y,y)
endif
store (mark,epor,enut,den1)
if(diagnos) then ! quantities assisting diagnosis
store (imb1,pcor,v1sl,v1ad) ! if needed
fiinit(du1p)=0
fiinit(dv1p)=0
fiinit(dw1p)=0
relax(du1p,linrlx,0.1)
relax(dv1p,linrlx,0.1)
relax(dw1p,linrlx,0.1)
varmax(du1p)=100
varmax(dv1p)=100
varmax(dw1p)=100
endif
save7end
************************************************************
Group 8. Terms & Devices
************************************************************
Group 9. Properties
************************************************************
Echo InForm settings for Group 9
save9begin
! properties are taken as those of air at 20 degK, except that
! density increases with mass fraction of gas in proportion to
! the molecular-weight difference
if(releas) then
(property den1 is 1.189*(1.+gas*(molwt/28.9-1.)))
else
(property den1 is 1.189)
endif
(property enul is 1.544e-5)
if(:tmodel:.eq.lvel) then
turmod(lvel)
liter(ltls) = 1000; endit(ltls) = 0.0
output(wdis.y,n,n,n,n,n)
endif
if(:tmodel:.eq.ke) then
turmod(kemodl)
solutn (ke, y, y, y, n, n, y )
solutn (ep, y, y, y, n, n, y )
output (ke, y, n,y, y, y, y)
output (ep, y, n,y, y, y, y)
varmin(ke)=0.005*windvel**2
varmin(ep)=10.*varmin(ke)**1.5
endif
save9end
************************************************************
Group 10.Inter-Phase Transfer Processes
************************************************************
Group 11.Initialise Var/Porosity Fields
************************************************************
Echo InForm settings for Group 11
save11Begin
patch(whole,cell,1,nx,1,ny,1,nz,1,lstep)
(initial of u1 at whole is :uprofl:)
mesgm(u and v initial profiles
(initial of v1 at whole is :vprofl:)
if(releas) then
(initial of gas at whole is 0.0)
endif
if(restart) then
restrt(all)
endif
save11end
************************************************************
Group 13. Boundary & Special Sources
************************************************************
Echo InForm settings for Group 13
save13begin
************************ patches *********************************
patch(zerov,cell,1,nx,1,ny,1,nz,1,1) ! a patch used during
! testing
coval(zerov,v1,1.e6,0.0)
do ixx=1,nx
+ patch(zerow:ixx:,cell,:ixx:,:ixx:,1,ny,1,nz,1,1)
! used during testing
+ coval(zerow:ixx:,w1,1.e9,0.0)
enddo
if(releas) then
patch(gravity,phasem,1,nx,1,ny,1,nz,1,lstep) ! when density varies
(source of w1 at gravity is coval(fixflu,-9.81*(den1-1.189)))
endif
patch(ceiling,high,1,nx,1,ny,nz,nz,1,1) ! fixed u and v at ceilinq
coval(ceiling,p1,fixp,0.0)
(source of u1 at ceiling is coval(fixval,:uprofl:))
(source of v1 at ceiling is coval(fixval,:vprofl:))
(source of gas at ceiling is coval(fixval,0.0))
patch(ceiling2,high,1,nx,1,ny,nz-1,nz-1,1,1) ! some in and outflow
(source of w1 at ceiling2 is coval(1.e3,0.0)) ! alllowed by ceiling2
(source of gas at ceiling2 is coval(fixval,0.0))
patch(floor,low,1,nx,1,ny,1,1,1,1) ! a rough surface, with friction
! factor floorco
(source of u1 at floor is coval(den1*vabs*floorco,0.0))
(source of v1 at floor is coval(den1*vabs*floorco,0.0))
********************** In-Form objects **************************
! the use of In-Form objects
! note that infob_1, infob_2, infob_3 must appear in numerical
! order.
! otherwise unexpected things happen to solution and print-out
! The use of io, and io=io+1 for each activated object, allows
! this. Note however that this use if io as a variable which can
! appear on both side of the = sign means that it is not a PIL
! variable in the sense understood by PRELUDE.
!---------------------------------------------------------------
io=0
io=io+1 ! in-form object 1
patch(infob:io:,west,1,1,1,ny,1,nz,1,1) ! infob_:io:, the west
! boundary
xp=0
yp=0
zp=0
xs=domlong * xfrac(1) ! to make the box one cell thick
ys=domwide
zs=domhigh
if(wal) then
real(xpw,ypw,zpw,xsw,ysw,zsw) ! wall parameters for use in BOUND
xpw=indist+walldist-1 ! in order that the -box feature
xsw=wallthck+2 ! may be used to avoid enforcing
ypw=-domwide ! boundary conditions at the ends
ysw=domwide *3 ! of the wall1
zpw=0.0
zsw=wallhigh
endif
subwal=f
CALL BOUND
io=io+1 ! next infob
patch(infob:io:,east,nx-2,nx,1,ny,1,nz,1,1) ! infob_:io:, the
! east boundary
iii=nx-2
xp=xfrac(iii)
xp=domlong*xp
yp=0
zp=0
iii=nx-1
xs=domlong-xp
ys=domwide
zs=domhigh
subwal=f
CALL BOUND
if(ny.gt.1.or.windangl.eq.0.0) then ! to allow testng with ny=1
io=io+1 ! next in-form object, the south boundary
patch(infob:io:,south,1,nx,1,1,1,nz,1,1)
xp=0
yp=0
zp=0
xs=domlong
ys=domwide * yfrac(1) ! to make the box one cell thick
zs=domhigh
subwal=wal ! subtract wall box if wal
CALL BOUND ! note that PIL subroutine names must be upper-case
io=io+1
patch(infob:io:,north,1,nx,ny-1,ny,1,nz,1,1) ! infob_:io:, the
! north boundary
iii=ny-2
xp=0
yp=yfrac(iii)
yp=domwide*yp
zp=0
xs=domlong
ys=domwide-yp
zs=domhigh
subwal-wal ! subtract wall box if wal
CALL BOUND
endif ! end of if(ny.gt.1)
if(wal) then ! next infob the wall
io=io+1
patch(infob:io:,cell,ixwal1-1,ixwal2+1,1,ny,1,nz/2,1,1)
xp=indist+walldist
xs=wallthck
yp=-domwide
ys=domwide *3
zp=0.0
zs=wallhigh
(infob at infob:io: is box(xp,yp,zp,xs,ys,zs,al,be,ga) with infob_:$
io:)
CALL VELZERO
endif
if(releas) then
io=io+1 ! next infob, the gas release
mesg(at release, infob no is :io: gasflo=:gasflo:
patch(infob:io:,cell,1,nx/2,ny/4,3*ny/4+1,1,nz/2,1,1)
patch(infob:io:,cell,1,nx,1,ny/4,1,nz,1,1)
(infob at infob:io: is point(:indist:,:domwide/2:,:pipezpos+pipedi$
am/2:,0.01) with infob_:io:)
(source of p1 at infob:io: is :gasflo: with infob_:io:)
(source of p1 at infob:io: is 1.e6 with infob_:io:)
(source of u1 at infob:io: is :gasvel:*:cos(relangl): with infob_:i$
o:!only)
(source of w1 at infob:io: is :gasvel:*:sin(relangl): with infob_:i$
o:!only)
(source of gas at infob:io: is 1.0 with infob_:io:!only)
endif
if(pipe) then
io=io+1 ! next infob, the southern half of the gas-carrying
! pipe
patch(infob:io:,phasem,1,nx/2,1,ny/2-2,1,nz/2,1,1)
real(rad)
rad=:pipediam/2:
(infob at infob:io: is ellpsd(:indist:,:domwide/2:,:pipezpos:+:rad:$
,:rad:,100.0,:rad:,0.,0.,0.) with infob_:io:)
CALL VELZERO
io=io+1 ! next infob, the northern part of the gas-carrying pipe
patch(infob:io:,phasem,1,nx/2,ny/2+1,ny,1,nz/2,1,1)
(infob at infob:io: is ellpsd(:indist:,:domwide/2:,:pipezpos:+:rad:$
,:rad:,100.0,:rad:,0.,0.,0.) with infob_:io:)
CALL VELZERO
endif
save13end
************************************************************
Group 15. Terminate Sweeps
************************************************************
Echo InForm settings for Group 15
save15begin
LSWEEP = 50
RESFAC = 1.000000E-03
save15end
************************************************************
Group 16. Terminate Iterations
************************************************************
Echo InForm settings for Group 16
save16begin
liter(p1) = 200 ;liter(gas ) = 200
endit(gas) = 0.0
save16end
************************************************************
Group 17. Relaxation
************************************************************
Group 18. Limits
save18begin
save18end
************************************************************
Group 19. EARTH Calls To GROUND Station
save19begin
conwiz = t ! use convergence wizard
calfor = t ! calculate forces
isg50 = 1 ! enforce endpause
isg52 = 1 ! plot max and min values as convergence indicators
save19end
************************************************************
Group 20. Preliminary Printout
************************************************************
Group 21. Print-out of Variables
************************************************************
Group 22. Monitor Print-Out
************************************************************
Group 23.Field Print-Out & Plot Control
************************************************************
Echo InForm settings for Group 23
save23begin
xzpr = t
tstswp=-1
save23end
************************************************************
Group 24. Dumps For Restarts
************************************************************
save25begin
GVIEW(P,0.000000E+00,-1.000000E+00,0.000000E+00)
GVIEW(UP,0.000000E+00,0.000000E+00,1.000000E+00)
Monitor MONITOROBJECT
> DOM, MONIT, 5, 15, 21
> DOM, SCALE, 1.000000E+00, 1.000000E+00, 1.000000E+00
> DOM, SNAPSIZE, 1.000000E-02
> DOM, SIZE, :domlong:,:domwide: , :domhigh:
save25end
debug=t
flag=t
izdb2=nz
iswdb2=lsweep
ENDMAIN
! start of subroutines (upper-case names needed)
SUBROUTINE CHINT
mesgm(:item:=:intger: press n if not ok
readvdu(ans,char,y)
if(:ans:.ne.y) then
label int
mesg(insert desired integer value for :item:
readvdu(:item:,int,12345)
if(:item:.eq.12345) then
goto int
endif
goto input
endif
ENDSUB
SUBROUTINE CHREAL
mesgm(:item:=:valu: press n if not ok
readvdu(ans,char,y)
if(:ans:.ne.y) then
label real
mesg(insert desired real value for :item:
readvdu(:item:,real,.12345)
if(:item:.eq..12345) then
goto real
endif
goto input
endif
ENDSUB
SUBROUTINE CHBOOL
mesgm(:item:=:bool: press n if not ok
readvdu(ans,char,y)
if(:ans:.eq.n) then
goto input
endif
ENDSUB
SUBROUTINE CHCHAR
mesgm(:item:=:ch: press n if not ok
readvdu(ans,char,y)
if(:ans:.ne.y) then
label char
mesg(insert desired character value for :item:
readvdu(:item:,char,pqrst)
if(:item:.eq.pqrst) then
goto char
endif
goto input
endif
mesgm(:item:=:ch: return if ok, else insert characters
readvdu(:item:,char,ch)
if(:item:.ne..pqrs) then
goto input
endif
ENDSUB
SUBROUTINE BOUND
al=0.0
be=0.0
ga=0.0
if(subwal) then ! subtract wall box for north and south
(infob at infob:io: is box(xp,yp,zp,xs,ys,zs,al,be,ga)-$
box(xpw,ypw,zpw,xsw,ysw,zsw,0.,0.,0.) with infob_:io:)
else
(infob at infob:io: is box(xp,yp,zp,xs,ys,zs,al,be,ga)$
with infob_:io:)
endif
(source of P1 at infob:io: is -1.e10*p1 with infob_:io:!line)
(source of U1 at infob:io: is :uprofl: with infob_:io:!fixv)
(source of V1 at infob:io: is :vprofl: with infob_:io:!fixv)
(source of w1 at infob:io: is 0. with infob_:io:!fixv)
if(releas) then
(source of gas at infob:io: is 0.0 with infob_:io:!only)
endif
if(tmodel.eq.ke) then
+ (source of ke at infob:io: is same with infob_:io:!only)
+ (source of ep at infob:io: is same with infob_:io:!only)
endif
(initial of mark at infob:io: is :io:+1 with infob_:io:)
ENDSUB
SUBROUTINE VELZERO
(source of u1 at infob:io: is 1.e9*(0.0-u1) with infob_:io:!line)
! why did coval not work
(source of v1 at infob:io: is 1.e9*(0.0-v1) with infob_:io:!line)
! why did coval not work ?
(source of w1 at infob:io: is 1.e9*(0.0-w1) with infob_:io:!line)
if(:tmodel:.eq.lvel) then
(source of ltls at infob:io: is 1.e9*(0.0-ltls) with infob_:io:!lin$
e)
endif
(initial of mark at infob:io: is 1.0 with infob_:io:)
ENDSUB
STOP
VRV USE
* Start of frame
* Setting display switches
TEXT CLEAR
AXIS ON
CELPOS OFF
CONTOUR SCALE ON
POSITION CONTOURKEY 2.252252E-02 8.389262E-02
TEXT ON
POSITION TITLE 2.590090E-01 8.993289E-01
POSITION PROBE 7.972973E-01 8.389262E-02
GRID OFF
* Domain scaling factors
SCALE 1.000000E+00 1.000000E+00 1.000000E+00
* Settings for current slice
PROBE 5.468750E+00 1.450000E+01 5.000000E-02; PROBE ON
SLICE Y
SLICE OUTLINE ON
* View and up directions
VIEW -6.841546E-02 -9.966097E-01 -4.569499E-02
UP -4.894052E-02 -4.239444E-02 9.979016E-01
* View centre
VIEW CENTRE 2.900758E+01 1.717579E+01 3.993689E+00
* View size
VIEW SIZE 1.704600E+01
* View perspective
VIEW DEPTH 3.000000E+00;VIEW TILT 0.8
DOMAIN ON
* Setting object visibility and painting status
OBJECT SHOW TYPE BLOCKAGE
OBJECT PAINT TYPE BLOCKAGE OFF
OBJECT WIREFRAME TYPE BLOCKAGE OFF
OBJECT SHOW TYPE INLET
OBJECT PAINT TYPE INLET OFF
OBJECT WIREFRAME TYPE INLET OFF
OBJECT SHOW TYPE OUTLET
OBJECT PAINT TYPE OUTLET OFF
OBJECT WIREFRAME TYPE OUTLET OFF
OBJECT SHOW TYPE USER_DEFINED
OBJECT PAINT TYPE USER_DEFINED OFF
OBJECT WIREFRAME TYPE USER_DEFINED OFF
VARIABLE u1
VECTOR ON
VECTOR COLOUR MULTI
VECTOR TYPE TOTAL
CONTOUR ON
CONTOUR BLANK ON
CONTOUR AVERAGE ON
SURFACE OFF
MINMAX OFF
CONTOUR OPAQUENESS 100
PAUSE
VRVEND