if(:ask:.eq.n) then
 goto proceed
endif
 
mesg(Defaults are LSTEP=:lstep:, NX=:nx:, NY=:ny:, NZ=:nz:
mesg(GAMA=:gama:, TLAST=:tlast:, uniform expansion,
mesgb(impervious walls. Are these settings OK? (Y/n)
readvdu(ans,char,y)
if(:ans:.eq.n) then
mesg(Change grid? (Y/n)
readvdu(ans,char,n)
if(:ans:.eq.y) then
MESG(LSTEP = :lstep: OK ?  If not, insert desired value
READVDU(LSTEP,INT,LSTEP)
LSTEP
MESG(NX= :NX: OK ?  If not, insert desired value
READVDU(NX,INT,NX)
NX
MESG(NY= :NY: OK ?  If not, insert desired value
READVDU(NY,INT,NY)
NY
MESG(NZ= :NZ: OK ?  If not, insert desired value
READVDU(NZ,INT,nz)
NZ
endif
mesg(The specific-heat ratio, GAMA, = :gama: OK?
mesg( If not, insert the desired value
READVDU(gama,real,gama)
GAMA
 
mesg( TLAST = :TLAST: OK?
mesg( If not, insert the desired value
READVDU(TLAST,real,TLAST)
TLAST
mesg(Change expansion rates, eg to compression? (y/N)
readvdu(ans,char,n)
if(:ans:.eq.y) then
mesg(Expansion ratios per time step in x, y and z directions are
mesg(currently :rsg8:, :rsg10:, and :rsg12:
mesg(Values smaller than 1 dictate compression. Enter three values.
readvdu(rsg8,real,rsg8)
rsg8
readvdu(rsg10,real,rsg10)
rsg10
readvdu(rsg12,real,rsg12)
rsg12
endif
mesg(Open the expanding boundaries ? (y/N)
readvdu(ans,char,n)
if(:ans:.eq.y) then
open=t
open
endif
endif
label proceed
    GROUP 2. Transience; time-step specification
GRDPWR(T,LSTEP,TLAST,1.0)
 
    GROUP 6. Body-fitted coordinates or grid distortion
    GROUP 7. Variables stored, solved & named
 
char(movebfc);movebfc=$B005
#movebfc
    The grid will be specified by way of coding in GROUND,
    and will vary with time.
 
    GROUP 9. Properties of the medium (or media)
RHO1=COMPRESS;RHO1B=1./GAMA;PRESS0=1.E5;RHO1A=1./PRESS0**RHO1B
DRH1DP=COMPRESS
    GROUP 11. Initialization of variable or porosity fields
FIINIT(P1)=0.0;FIINIT(U1)=0.0;FIINIT(V1)=0.0;FIINIT(W1)=0.0
    GROUP 15. Termination of sweeps
LSWEEP=50;selref=t;resfac=0.01
    GROUP 16. Termination of iterations
LITER(P1)=50
 
    GROUP 17. Under-relaxation devices
real(dtf); dtf=0.01*tlast
RELAX(U1,FALSDT,dtf);RELAX(V1,FALSDT,dtf);RELAX(W1,FALSDT,dtf)
RELAX(P1,LINRLX,0.5)
 
    GROUP 19. Data communicated by satellite to GROUND
 
   
    Set the initial box size:
SPEDAT(GRIDS, CASE,  C, BOX)
SPEDAT(GRIDS, XCMIN, R, 0.0)
SPEDAT(GRIDS, XCMAX, R, 1.0)
SPEDAT(GRIDS, YCMIN, R, 0.0)
SPEDAT(GRIDS, YCMAX, R, 2.0)
SPEDAT(GRIDS, ZCMIN, R, 0.0)
SPEDAT(GRIDS, ZCMAX, R, 3.0)
    Set the time factors
if(nx.gt.1) then
 SPEDAT(GRIDS, TFCXMX, R, 1.1)
 SPEDAT(GRIDS, TFCXMN, R, 1.0)
else
 SPEDAT(GRIDS, TFCXMX, R, 1.0)
 SPEDAT(GRIDS, TFCXMN, R, 1.0)
endif
 
if(nx.gt.1) then
 SPEDAT(GRIDS, TFCYMX, R, 1.1)
 SPEDAT(GRIDS, TFCYMN, R, 1.0)
else
 SPEDAT(GRIDS, TFCYMX, R, 1.0)
 SPEDAT(GRIDS, TFCYMN, R, 1.0)
endif
 
if(nx.gt.1) then
 SPEDAT(GRIDS, TFCZMX, R, 1.1)
 SPEDAT(GRIDS, TFCZMN, R, 1.0)
else
 SPEDAT(GRIDS, TFCZMX, R, 1.0)
 SPEDAT(GRIDS, TFCZMN, R, 1.0)
endif
    Set the slopes
SPEDAT(GRIDS, DXCDZC, R, 0.0)
SPEDAT(GRIDS, DYCDXC, R, 0.0)
SPEDAT(GRIDS, DYCDZC, R, 0.0)
 
    GROUP 22. Spot-value print-out
TSTSWP=-1;IXMON=NX/2;IYMON=NY/2;IZMON=NZ/2
    GROUP 23. Field print-out and plot control
NPLT=1
PATCH(LONGPLOT,PROFIL,1,nx,1,1,1,1,1,LSTEP)
PLOT(LONGPLOT,P1,0.0,0.0);PLOT(LONGPLOT,U1,0.0,0.0)
 
CSG1=PHI;CSG2=XYZ;IDISPA=1
 
IPLTF=3;XZPR=T;itabl=1
IF(:OPEN:.EQ.T) THEN
 IF(NX.GT.1) THEN
+ PATCH(E,CELL,NX,NX,1,NY,1,NZ,1,LSTEP)
+ COVAL(E,P1,FIXVAL,0.0)
 ENDIF
 IF(NY.GT.1) THEN
+ PATCH(N,CELL,1,NX,NY,NY,1,NZ,1,LSTEP)
+ COVAL(N,P1,FIXVAL,0.0)
ENDIF
 IF(NZ.GT.1) THEN
+ PATCH(H,CELL,1,NX,1,NY,NZ,NZ,1,LSTEP)
+ COVAL(H,P1,FIXVAL,0.0)
 ENDIF
ENDIF