TALK=F;RUN( 1, 1)
  
    GROUP 1. Run title and other preliminaries
TEXT(Stafford's saline-layer experiment
TITLE
  DISPLAY
  A heavier saline solution lies below a lighter fresh-water layer.
  Both have the same temperature, and are at rest.
 
  The lower fluid is heated, for a limited period. It becomes lighter 
  than the upper fluid, and therefore rises through it.
  
  However, after some time it sinks again because it loses heat 
  much faster than it loses salt.
 
  The phenomena are described in the lecture which may be seen by 
   clicking here. 
  
  The setting IDISPA=1 ensures that a parphi (or parada) file is
  created enabling the development in time to be viewed by way of
  PHOTON, with the z-axis being used for time.
 
  A PHOTON USE file is supplied.
  ENDDIS
    GROUP 2. Transience; time-step specification
STEADY=F;GRDPWR(T,200,500.0,1.0)
    GROUP 3. X-direction grid specification
nx=10
nx=20
mesg(NX = :NX: OK? If not, enter new value...
  readvdu(nx,int,nx)
GRDPWR(X,nx,0.1,1.0)
    
    
    GROUP 7. Variables stored, solved & named
ONEPHS=F
SOLVE(P1,U1,U2,R1,R2)
store(den1,den2)
terms(h1,n,p,p,p,p,p)
terms(h2,n,p,p,p,p,p)
terms(r1,p,p,n,p,p,p)
terms(r2,p,p,n,p,p,p)

    GROUP 9. Properties of the medium (or media)
RHO1=LINSCAL; RHO2=LINSCAL;IBUOYB=14;IBUOYC=16
RHO1A=1000; RHO1B=-1; RHO1C=1
RHO2A=1000; RHO2B=-1; RHO2C=1
varmin(den1)=rho1a*0.9;varmin(den2)=rho1a*0.9
varmin(r1)=1.e-6;varmin(r2)=1.e-6

mesg(RHO2 = :rho2:  OK? If not, enter new value...
  readvdu(rho2,real,rho2)
    
    GROUP 10. Inter-phase-transfer processes and properties
CFIPS=10.0
cint(h1)=0.5e-3;  cint(h2)=1.0
cint(c1)=0.5e-5;  cint(c2)=1.0
    
    GROUP 11. Initialization of variable or porosity fields
real(solfrac)
solfrac=0.5
mesg(heavy-material volume fraction = :solfrac:
mesg(OK? If not, enter new value...
  readvdu(solfrac,real,solfrac)
  
  case 1: sedimentation
  fiinit(r1)=1.0-solfrac;fiinit(r2)=solfrac
 
  case 2: heavy in top half light in bottom
solve(h1,h2,c1,c2)
fiinit(r1)=1.0; fiinit(r2)=0.0
fiinit(h1)=0.0; fiinit(h2)=0.0
fiinit(c1)=1.0; fiinit(c2)=0.0
iniadd=f

              ---------------------------------
                     fluid 2
                           h init = 0  c init = 0
              ---------------------------------
                     fluid 1
                           h init = 0  c init = 1
              ---------------------------------
patch(start,inival,nx/2+1,nx,1,1,1,1,1,1)
coval(start,r1,0.0,0.0); coval(start,r2,0.0,1.0)  
  
    GROUP 13. Boundary conditions and special sources
real(gravacc)
gravacc=9.81
mesg(gravitational acceleration = :gravacc: m/s ** 2
mesg(OK? If not, enter new value...
  readvdu(gravacc,real,gravacc)
 
PATCH(GRAVITY,PHASEM,1,NX,1,1,1,1,1,LSTEP)
COVAL(GRAVITY,u1,FIXFLU,-gravacc)
COVAL(GRAVITY,u2,FIXFLU,-gravacc)

PATCH(TOP,CELL,NX,NX,1,1,1,1,1,LSTEP)
COVAL(TOP,P1,1.0,0.0);COVAL(TOP,P2,RHO2A,0.0)

patch(heat,phasem,1,nx/2,1,ny,1,nz,1,lstep/4)
coval(heat,h1,fixflu,0.02)

    GROUP 15. Termination of sweeps
LSWEEP=100;resfac=1.e-4
    GROUP 17. Under-relaxation devices
RELAX(R2,LINRLX,0.25);RELAX(R1,LINRLX,0.25)
RELAX(P1,LINRLX,0.25)
RELAX(U1,FALSDT,0.1);RELAX(U2,FALSDT,0.1)
RELAX(H1,FALSDT,1);RELAX(H2,FALSDT,1)
RELAX(C1,FALSDT,1);RELAX(C2,FALSDT,1)
SPEDAT(SET,GXMONI,TRANSIENT,L,F)    
    GROUP 21. Print-out of variables
OUTPUT(R1,N,N,N,N,N,N)
OUTPUT(P1,Y,Y,Y,Y,Y,Y);OUTPUT(R2,Y,Y,Y,Y,Y,Y)
OUTPUT(u1,Y,Y,Y,Y,Y,Y);OUTPUT(u2,Y,Y,Y,Y,Y,Y)
ntprin=lstep/10
    
    GROUP 22. Spot-value print-out
NTPRIN=lstep/5;TSTSWP=-1;IXMON=nx/2
    GROUP 23. Field print-out and plot control
IPLTF=2;IPLTL=LSWEEP;XZPR=T;IPROF=3;ITABL=3
PATCH(PROFILE,PROFIL,1,NX,1,1,1,1,1,LSTEP)
PLOT(PROFILE,P1,0.0,0.0);PLOT(PROFILE,R1,0.0,0.0)
PLOT(PROFILE,u1,-1.0,-1.0);PLOT(PROFILE,u2,-1.0,-1.0)
PATCH(TIMEPLOT,PROFIL,NX,NX,1,1,1,1,1,LSTEP)
PLOT(TIMEPLOT,R1,0.0,0.0)
    GROUP 24. Dumps for restarts
IDISPA=1
  
    photon use
    p;parphi
   1000 1 1
 
   view y
   gr ou y 1
   con r1 y 1 fi;0.1
   msg lower-phase concentration
   msg horizontal coordinate is time, increasing to the right
   
   pause;con off;red;con h1 y 1 fi;0.1
   msg lower-phase temperature
   msg horizontal coordinate is time, increasing to the right
   
   pause;con off;red;con h2 y 1 fi;0.1
   msg higher-phase temperature
   msg horizontal coordinate is time, increasing to the right
   
   pause;con off;red;con den1 y 1 fi;0.1
   msg lower-phase density
   msg horizontal coordinate is time, increasing to the right
   
   pause;con off;red;con den2 y 1 fi;0.1
   msg higher-phase density
   msg horizontal coordinate is time, increasing to the right
   
   pause;con off;red;con u1 y 1 fi;0.1
   msg lower-phase velocity
   msg horizontal coordinate is time, increasing to the right
   
   pause;con off;red;con u2 y 1 fi;0.1
   msg higher-phase velocity
   msg horizontal coordinate is time, increasing to the right
   
   pause;con off;red;con c1 y 1 fi;0.1
   msg lower-phase c-value
   msg horizontal coordinate is time, increasing to the right
   
   pause;con off;red;con c2 y 1 fi;0.1
   msg higher-phase c-value
   msg horizontal coordinate is time, increasing to the right
   
   *pause;*con off;*red;*con p1 y 1 fi;0.1
   *msg pressure in both phases
   *msg horizontal coordinate is time, increasing to the right
 
   enduse
STOP