PHOTON USE
  p;parphi
  50 50 1
 
  set ref vec .1
  *  view z
  msg water levels at successive times. Press any key to pause
  do kk=1,m
  con p1 z kk fi;1
  enddo
  pause;cl
  msg x-direction velocity contours. Press any key to pause 
  do kk=1,m
  con u1 z kk fi;1000
  enddo
  pause;cl
  msg y-direction velocity contours. Press any key to pause 
  do kk=1,m
  con v1 z kk fi;1000
  enddo
  pause;cl
 
  msg pollutant contours. Press any key to pause
  do kk=1,m
  con pol z kk fi;01
  enddo
  pause;cl
 
  ENDUSE
#cls  
  DISPLAY
  The transient flow into and out of a harbour is simulated by means
  of the shallow-water compressible-flow analogy.
  
  This is effected by:
   
   * setting RHO2 = -rhonom, as a signal to GXRHO 
     (in file GXDEBS.HTM); and 
   * by changing PRESS0 from the value set in the SEASURF PIL 
     fragment, namely :  9.81 * ZWLAST * RHONOM
     to               :  9.81 * DZUPPER * RHONOM
     where DZUPPER is the depth of the upper cells.
  
  
  ENDDIS
#pause    
    GROUP 1. Run title and other preliminaries
TEXT(Pollutant Flow In A Tidal Harbour
TITLE
REAL(CYCLETIM, TIDERANGE)
INTEGER(NCYCLES,NSTEPS)
TIDERANGE=0.5; CYCLETIM=12*3600; NCYCLES=1; NSTEPS=200
 
    GROUP 2. time grid specification
STEADY=F; LSTEP=NSTEPS*NCYCLES
GRDPWR(T,LSTEP,NCYCLES*CYCLETIM,1.0)
XULAST=100; YVLAST=100; ZWLAST =2

NX=40; NY=40; NZ=1
  ** use macros to specify a uniform grid and solve for velocities
                                                    and pressure
#UNIGRID
#solvel

    GROUP 7. Variables stored, solved & named
STORE(RHO1,PRPS)
SOLVE(POL)

    GROUP 9. Properties of the medium (or media)
  rho1=1000 * [(p1-p0) / ( 1000 * g * zwlast)] ** 0.5
RHO2=1.0; RHO1=COMPRESS; DRH1DP=COMPRESS
REAL(RHONOM); RHONOM=1000.0
RHO1B=0.5 ; PRESS0= 9.81 * ZWLAST * RHONOM
RHO1A=RHONOM / PRESS0 ** RHO1B; RHO1C=0.0
 

    GROUP 11. Initialization of variable or porosity fields
FIINIT(RHO1)=RHONOM
  ** introduce seawall, jetty and island by way of blockages
PATCH(SEAWALL,INIVAL,NX/2,NX/2,NY/5,4*NY/5,1,1,1,LSTEP)
COVAL(SEAWALL,PRPS,0.0,198.0)

PATCH(JETTY,INIVAL,3*NX/4,3*NX/4,4*NY/5,NY,1,1,1,LSTEP)
COVAL(JETTY,PRPS,0.0,198.0)

PATCH(ISLAND,INIVAL,4*NX/5,9*NX/10,NY/5,2*NY/5,1,1,1,LSTEP)
COVAL(ISLAND,PRPS,0.0,198.0)
  
  ** Tidal boundary by way of time patch
PATCH(TIMINLET,WEST,1,1,1,NY,1,1,1,LSTEP)
COVAL(TIMINLET,P1,FIXVAL,GRND3); ITIMA=NSTEPS; ITIMB=1; ITIMC=NSTEPS
  pressure amplitude
TIMA=9.81*RHONOM*TIDERANGE

  ** bottom friction
PATCH(BOTTOM,HWALL,1,NX,1,NY,1,1,1,LSTEP)
COVAL(BOTTOM,U1,1.0, 0); COVAL(BOTTOM,V1,1.0, 0)

  ** pollutant source
PATCH(POLLUTNT,CELL,3*NX/4,3*NX/4,NY/2,NY/2+1,1,1,1,LSTEP)  
COVAL(POLLUTNT,POL,FIXFLU,100)

  GROUP 15. Termination of sweeps
LSWEEP=10
 
    GROUP 16. Termination of iterations
LITER(P1)=50
 
    GROUP 17. Under-relaxation devices
REAL(DTF);DTF=CYCLETIM*0.0001
RELAX(U1,FALSDT,DTF); RELAX(V1,FALSDT,DTF)
SPEDAT(SET,GXMONI,TRANSIENT,L,F) 
    GROUP 18. Limits on variables or increments to them
VARMAX(RHO1)=5.0*RHONOM; VARMIN(RHO1)=0.2*RHONOM
 
    GROUP 22. Spot-value print-out
IXMON=3*NX/4; IYMON=3*NY/4; IZMON=1; TSTSWP=-2; NPLT=1
 
    GROUP 23. Field print-out and plot control
  ** line-prinnter plots
PATCH(TIMPLT1,PROFIL, 1,1, 1,1, 1,1, 1,LSTEP)  
COVAL(TIMPLT1,U1,0.0,0.0); COVAL(TIMPLT1,RHO1,0.0,0.0)
COVAL(TIMPLT1,P1,0.0,0.0)

PATCH(TIMPLT2,PROFIL, NX,NX, 1,1, 1,1, 1,LSTEP)  
COVAL(TIMPLT2,RHO1,0.0,0.0)
COVAL(TIMPLT2,P1,0.0,0.0)

PATCH(TIMPLT3,PROFIL, 3*NX/4,3*NX/4, 1,1, 1,1, 1,LSTEP)  
COVAL(TIMPLT3,U1,0.0,0.0); COVAL(TIMPLT3,RHO1,0.0,0.0)
COVAL(TIMPLT3,P1,0.0,0.0)

 
  ** settings enabling PHOTON to plot every eighth time step
 IDISPA=NSTEPS/8;IDISPB=1;IDISPC=LSTEP  
STOP