PHOTON USE
  p;;;;;
 
  msg 1. CLDA tests
  gr ou z 1
  vec z 1
  msg contours for the upwind difference scheme
  con upwc z 1 fi;0.02
  msg contours for the upwind difference scheme.
  msg Press  to continue
  pause;con off;vec off;red
  msg contours for the CLDA
  con cnor z 1 fi;0.02
  msg            -
  msg Press  to continue
  pause
  msg            -
  msg Press e to END
  enduse
 
   GROUP 1. Run title and other preliminaries
TEXT(GXCLDA tests :               N301
TITLE
  DISPLAY
  CONSERVATIVE LOW-DISPERSION ALGORITHM (CLDA) - flow at various
  angles 2-dimensional (x-y), Cartesian, steady/unsteady, elliptic
  simulation.This demonstration of the CLDA (conservative low-
  dispersion algorithm) concerns the flow of two fluids, "red"
  and "blue" in a 2D domain. Physical diffusion is absent, so that
  the apparent mixing is a consequence of numerical diffusion alone.
 
  The subroutine GXCLDA.F id used.
 
  Several different runs can be made, including a transient.
  If the transient is chosen, PHOTON-readable files P1, P2, etc
  will be created; but the PHOTON USE commands supplied at the top
  of this Q1 call in only that for the last time step. To see the
  others type p;P1 , P;P2, etc
 
  ENDDIS
char(ans)
mesg(Press RETURN to continue
readvdu(ans,char,y)
 
  The 3D-unsteady case remains to be coded.
 
  No BFC testing has been done.
 
char(exam)
REAL(UVEL,VVEL)
boolean(edgeso,cornso,shuttop)
 
  data:
  Input data **************************************
LSTEP=10;NX=10;NY=10
STEADY=t
UVEL=1.0;VVEL=1.0
xulast=1.0;yvlast=1.0
edgeso=t;shuttop=t
mesg(Please enter a number (1-9)
mesg( 1. flow diagonal to the grid
mesg( 2. vertical velocity =0.5 * horizontal velocity
mesg( 3. Point source in corner
mesg( 4. vertical vel. =0.5 * horiz. vel. with corner source
mesg( 5. vertical flow with corner source
mesg( 6. diagonal flow with corner source and shuttop
mesg( 7. as for 6, but finner grids
mesg( 8. Fluid flow is steady; scalar transport is transient
mesg( 9. A limited duration source released from the origin
 
readvdu(exam,char,1)
 
if(:exam:.eq.1) then
 steady=t;uvel=1.0;edgeso=t;shuttop=f
endif
if(:exam:.eq.2) then
 steady=t;uvel=0.5;edgeso=t;shuttop=f
endif
if(:exam:.eq.3) then
 steady=t;uvel=1.0;edgeso=f;shuttop=f
endif
if(:exam:.eq.4) then
 steady=t;uvel=0.5;edgeso=f;shuttop=f
endif
if(:exam:.eq.5) then
 steady=t;uvel=0.0;edgeso=f;shuttop=f
endif
if(:exam:.eq.6) then
 steady=t;uvel=1.0;edgeso=f;shuttop=t
endif
if(:exam:.eq.7) then
 steady=t;uvel=1.0;edgeso=f;shuttop=t
 nx=20;ny=20
endif
if(:exam:.eq.8) then
 steady=f;uvel=1.0;edgeso=t;shuttop=f
 SPEDAT(SET,GXMONI,TRANSIENT,L,F)
endif
if(:exam:.eq.9) then
 steady=f;uvel=1.0;edgeso=f;shuttop=f
 SPEDAT(SET,GXMONI,TRANSIENT,L,F)
endif
if(:exam:.eq.10) then
 steady=f;uvel=1.0;edgeso=f;shuttop=t
 SPEDAT(SET,GXMONI,TRANSIENT,L,F)
endif
 
 
   Input data **************************************
 
if(edgeso) then
 cornso=f
else
 cornso=t
endif
if(steady) then
 lstep=1
endif
mesg( Data for this run are:-
mesg( exam = :exam:
mesg( nx=:nx:, ny=:ny:, steady=:steady:, lstep=:lstep:
mesg( uvel=:uvel:, vvel=:vvel:, xulast=:xulast:, yvlast=:yvlast:
mesg( edgeso = :edgeso:, cornso = :cornso:, shuttop = :shuttop:
mesg(
 
    GROUP 2. Transience; time-step specification
GRDPWR(T,LSTEP,1.0,1.0)
 
    GROUP 3. X-direction grid specification
GRDPWR(X,NX,xulast,1.0)
 
    GROUP 4. Y-direction grid specification
GRDPWR(Y,NY,yvlast,1.0)
    GROUP 7. Variables stored, solved & named
Solve(p1)
if(nx.gt.1) then
 solve(u1)
endif
if(ny.gt.1) then
 solve(v1)
endif
solve(CNEW,CWES,CEAS,CSOU,CNOR,COLD,UPWC)
solutn(cnew,p,p,n,y,p,p)
solutn(cwes,p,p,n,y,p,p)
solutn(ceas,p,p,n,y,p,p)
solutn(cnor,p,p,n,y,p,p)
solutn(csou,p,p,n,y,p,p)
solutn(cold,p,p,n,y,p,p)
solutn(upwc,y,y,n,n,p,p)
if(nx.eq.1) then
 solutn(cwes,p,n,p,p,p,p)
 solutn(ceas,p,n,p,p,p,p)
endif
if(ny.eq.1) then
 solutn(csou,p,n,p,p,p,p)
 solutn(cnor,p,n,p,p,p,p)
endif
if(steady) then
 solutn(cnew,p,n,p,p,p,p)
 solutn(cold,p,n,p,p,p,p)
endif
if(.not.steady) then
 store(CNEO)
endif
 
    GROUP 8. Terms (in differential equations) & devices
TERMS(CNEW,N,N,N,n,Y,Y)
TERMS(CWES,N,N,N,n,Y,Y)
TERMS(CEAS,N,N,N,n,Y,Y)
TERMS(CSOU,N,N,N,n,Y,Y)
TERMS(CNOR,N,N,N,n,Y,Y)
TERMS(COLD,N,N,N,n,Y,Y)
 
    GROUP 11. Initialization of variable or porosity fields
FIINIT(U1)=UVEL;FIINIT(V1)=VVEL
FIINIT(CNOR)=0.0;FIINIT(CSOU)=0.0
FIINIT(CWES)=0.0;FIINIT(CEAS)=0.0
if(.not.steady) then
 fiinit(CNEO)=0.0
 FIINIT(COLD)=0.0
 FIINIT(CNEW)=0.0
endif
    GROUP 13. Boundary conditions and special sources
patch(westedge,west,1,1,1,ny,1,1,1,lstep)
coval(westedge,p1,fixflu,uvel)
coval(westedge,u1,onlyms,uvel)
coval(westedge,v1,onlyms,vvel)
if(edgeso) then
 coval(westedge,cwes,fixval,0.0)
endif
if(cornso) then
 patch(corn,cell,1,1,1,1,1,1,1,1)
 coval(corn,upwc,fixval,1.0)
 coval(corn,cwes,fixval,1.0)
 coval(corn,cold,fixval,1.0)
 coval(corn,cnew,fixval,1.0)
 coval(corn,ceas,fixval,1.0)
 coval(corn,csou,fixval,1.0)
 
 if(.not.steady) then
+ patch(corn3,cell,1,1,1,1,1,1,2,lstep)
+ coval(corn3,cold,fixval,0.0)
+ if(nx.gt.1) then
+   patch(corn4,cell,2,nx,1,1,1,1,1,lstep)
+   coval(corn4,cold,fixval,0.0)
+ endif
+ if(ny.gt.1) then
+   patch(corn5,cell,1,1,2,ny,1,1,1,lstep)
+   coval(corn5,cold,fixval,0.0)
+ endif
 endif
 if(ny.gt.1.and.uvel.gt.0.0) then
+ patch(westedg2,cell,1,1,2,ny,1,1,1,lstep)
+ coval(westedg2,cwes,fixval,0.0)
 endif
 patch(soutedg2,cell,2,nx,1,1,1,1,1,lstep)
 coval(soutedg2,csou,fixval,0.0)
endif
if(ny.gt.1) then
 PATCH(SOUTEDGE,SOUTH,1,NX,1,1,1,1,1,LSTEP)
 coval(soutedge,p1,fixflu,vvel)
 coval(soutedge,v1,onlyms,vvel)
 coval(soutedge,u1,onlyms,uvel)
 if(edgeso) then
+ COVAL(SOUTEDGE,CSOU,fixval,1.0)
 endif
endif
patch(eastedge,east,nx,nx,1,ny,1,1,1,lstep)
coval(eastedge,p1,fixp*uvel,0.0)
 
if(ny.gt.1) then
 if(.not.shuttop) then
+ patch(nortedge,north,1,nx,ny,ny,1,1,1,lstep)
+ coval(nortedge,p1,fixp*vvel,0.0)
 endif
endif
 
    The variable: "upwc"
 
if(edgeso) then
 COVAL(SOUTEDGE,upwc,onlyms,1.0)
 coval(westedge,upwc,onlyms,0.0)
endif
 
    The clda patch
patch(CLDA,cell,1,nx,1,ny,1,nz,1,lstep)
if(ny.gt.1) then
 coval(CLDA,cnor,grnd,grnd)
 coval(CLDA,csou,grnd,grnd)
endif
if(nx.gt.1) then
 coval(CLDA,ceas,grnd,grnd)
 coval(CLDA,cwes,grnd,grnd)
endif
if(.not.steady) then
 coval(CLDA,cnew,grnd,grnd)
 coval(CLDA,cold,grnd,grnd)
endif
 
    GROUP 19. Data communicated by satellite to GROUND
READQ1=T
 
    GROUP 21. Print-out of variables
itabl=1
output(p1,y,n,n,n,n,n)
output(u1,y,n,n,n,n,n)
output(v1,y,n,n,n,n,n)
 
nxprin=nx/5;nyprin=ny/5
NTPRIN=1
if(nx.eq.1) then
 output(u1,n,n,n,n,n,n);output(cwes,n,n,n,n,n,n)
 output(ceas,n,n,n,n,n,n)
endif
if(ny.eq.1) then
 output(v1,n,n,n,n,n,n);output(csou,n,n,n,n,n,n)
 output(cnor,n,n,n,n,n,n)
endif
if(steady) then
 OUTPUT(cold,N,N,N,N,N,N);OUTPUT(cnew,N,N,N,N,N,N)
 output(cnew,n,n,n,n,n,n)
endif
 
    GROUP 22. Spot-value print-out
IXMON=nx/2;IYMON=ny/2;UWATCH=T
 
    GROUP 23. Field print-out and plot control
orsiz=0.2
if(nx.eq.1) then
 patch(nxeq1,profil,1,1,1,ny,1,1,1,lstep)
 plot(nxeq1,cnew,0.0,1.0)
 plot(nxeq1,upwc,0.0,1.0)
endif
if(ny.eq.1) then
 patch(nyeq1,profil,1,nx,1,1,1,1,1,lstep)
 plot(nyeq1,cnew,0.0,1.0)
 plot(nyeq1,upwc,0.0,1.0)
endif
if(nx.gt.1.and.ny.gt.1) then
 PATCH(MAP,CONTUR,1,NX,1,NY,1,1,1,LSTEP)
 if(.not.steady) then
+ PLOT(MAP,CNEW,0.0,10);PLOT(MAP,COLD,0.0,10)
 endif
 PLOT(MAP,CWES,0.0,10);PLOT(MAP,CEAS,0.0,10)
 PLOT(MAP,CSOU,0.0,10);PLOT(MAP,CNOR,0.0,10)
 PLOT(MAP,UPWC,0.0,10)
endif
if(nx.gt.1.and.ny.gt.1.and..not.steady) then
 CSG1=PHI;CSG2=XYZ;IDISPA=2
endif
RELAX(U1,FALSDT,1.0);RELAX(V1,FALSDT,1.0)
LSWEEP=ny+1
FIINIT(P1)=1.0;liter(p1)=10
tstswp=-1