PHOTON USE
    p
    parphi



    up z
    msg Oil concentration (LOG10(oil)
    con logo x 1 fil;0.01
    mirror y
    msg Press  to continue
    pause
    cl
    msg Velocity vectors
    vec x 1 sh
    msg Press  to continue
    pause
    cl
    msg Effective viscosity
    con enut x 1 fil;0.01
    pause
    msg non-dimensional velocity
    cl
    con wdim x 1 fil;0.01
    msg Press e to END
    ENDUSE

    GROUP 1. Run title and other preliminaries
TEXT(Oil tanker rupture - Turbulent plume
TITLE

   DISPLAY
   This case simulates an oil tanker rupture. The calculation is stead 
   with a uniform oil release rate which can be varied by changing 
   the parameter, "rate". The parabolic mode with expanding 
   grid is employed for both economy and accuracy. The calculation
   includes a residence-time equation. 
   In-Form is used to set the buoyance force and to calculate the 
   non-dimensional velocity and oil concentration.   
   
   The result shows that the flow is a kind of turbulent plume. 
   ENDDIS
xulast=0.1 
    GROUP 4. Y-direction grid specification

CARTES=F
       ** The oil release rate can be varied by changing the 
          parameter, rate below
real(rate,rate1,radiu)
rate=20
rate1=rate/(3.6*24*20)
rate
radiu=(1.0**2*rate/2000)**0.5
radiu
NY=40; YVLAST=radiu
YFRAC(1)=-40.;YFRAC(2)=1.0/40.
   *** Linear grid expansion with slope DYGDZ
AZYV=1.0                 ! this dictates that yvlast is linear in z
REAL(DYGDZ); DYGDZ=0.33  ! dygdz is  d(yvlast)/dz
ZWADD=YVLAST/DYGDZ       ! so YVLAST = YVLAST at the inlet
                            + DYGDZ * (ZWLAST + ZWADD)
NZ=200

ZWLAST=3500
    GROUP 5. Z-direction grid specification
PARAB=T
    
   AZDZ=PROPY       ! propy = grnd2; means dz is proportional to y
   DZW1 = 0.077     ! the proportionality factor; so dz=dzw1 * yvlast
GRDPWR(Z,NZ,ZWLAST,1.)
    GROUP 7. Variables stored, solved & named
STORE(ENUT,LEN1)
SOLVE(P1,V1,W1,OIL)
           ** solve an equation for oil residence time
SOLVE(REST)
TERMS(REST,N,Y,N,P,P,P)
PATCH(ELAPSE,PHASEM,1,NX,1,NY,1,NZ,1,1)
COVAL(ELAPSE,REST,FIXFLU,1.0)
    
    GROUP 8. Terms (in differential equations) & devices
DIFCUT=0.0
 
    GROUP 9. Properties of the medium (or media)
      ** water density =1000 and oil density = 999.97
RHO1=1000.
TURMOD(KEMODL);KELIN=3
ENUL=1.e-6
FIINIT(KE)=0.1**2;FIINIT(EP)=0.1643*FIINIT(KE)**1.5/0.1
      ** dendif is the density difference between water and oil
REAL(DENDIF)
DENDIF=0.03

    GROUP 13. Boundary conditions and special sources
          ** buoyancy force 
PATCH(WHOLE,PHASEM,1,NX,1,NY,1,NZ,1,1)
(SOURCE OF W1 AT WHOLE IS 9.81*OIL*DENDIF WITH FIXFLU) 

       ** free stream
PATCH(NORTH,NORTH,1,1,NY,NY,1,NZ,1,1)
COVAL(NORTH,P1,FIXVAL,0.0)
coval(NORTH,REST,1.0,SAME)
       ** oil release point        
PATCH(LEAK,CELL,1,1,1,NY/2,1,1,1,1)
COVAL(LEAK,KE,FIXFLU,FIINIT(KE))
COVAL(LEAK,EP,FIXFLU,FIINIT(EP))
COVAL(LEAK,OIL,FIXFLU,:RATE1:)
       ** bottom of the sea
PATCH(OUTSIDE,cell,1,1,NY/2+1,NY,1,1,1,1)
COVAL(OUTSIDE,P1,1.E3,0.0)
COVAL(OUTSIDE,W1,ONLYMS,0.0)
    GROUP 14. Downstream pressure for PARAB=T
IPARAB=1
    GROUP 16. Termination of iterations
LITHYD=50
RELAX(V1,FALSDT,1.)
RELAX(W1,FALSDT,1.)
RELAX(REST,FALSDT,2.0)
    GROUP 19. Data communicated by SATELLITE to GROUND
  ** Select strain-rate for use in Mixing-Length model
DWDY=T
    GROUP 21. Print-out of variables
OUTPUT(P1,Y,Y,Y,Y,Y,Y);OUTPUT(V1,Y,Y,Y,Y,Y,Y)
OUTPUT(W1,Y,Y,Y,Y,Y,Y)
    GROUP 22. Monitor print-out
IZMON=10;IYMON=1;ITABL=1;NPLT=1;IPLTL=LITHYD;TSTSWP=-3
     ** parabolic file dumping

TSTSWP=-5;IDISPA=10;IDISPB=1;IDISPC=NZ
    GROUP 23. Field print-out and plot control
    
libref=402       
(stored var logo is log10(oil))

(store var wdim is w1*(zzw/:zwlast:)^0.3333)
(store var odim is oil*(zzw/:zwlast:)^1.6667)

patch(axis,profil,1,1,1,1,1,nz,1,1) 
coval(axis,w1,0.0,0.0)
coval(axis,oil,0.0,0.0)
coval(axis,logo,0.0,0.0)
coval(axis,wdim,0.0,0.0)
coval(axis,odim,0.0,0.0)
coval(axis,rest,0.0,0.0)