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