PHOTON USE
   p
   phi
   
   
   
   up z
   mirror y
   msg Oil concentration
   con oil x 1 sh;int 300
   msg Press  to continue
   pause
   msg Temperature contours
   con tem1 x 1 fil;0.01
   msg Press  to continue
   pause
   msg Velocity vectors
   vec x 1 sh
   msg Press e to END
   ENDUSE
   
Text(Oil tanker rupture; warm-water layer effect
TITLE
   DISPLAY   
   This case simulates an oil tanker rupture. The purpose of this 
   calculation is to show how a higher temperature in the layer
   near the sea surface affect the asending plume of the water-oil
   mixture. 
   
   The buoyancy force is modelled as a linear function of both 
   oil concentration and the temperature; and is set by In-Form.   
   
   This calculation focuses on the top layer near the surface and 
   the output from a separate parabolic calculation of the 
   turbulent plume is used as guidance to the input at the lower surface.
   The profiles of inflow velocity is therefore represented as a sine curve.
   
   The results show that the presence of a warm-water layer near 
   the sea surface causes the water-oil liquid to spread laterally,
   beneath the hot layer, perhaps never reaching the surface. 
   
   ENDDIS

NX=1; NY=40; NZ=40
xulast=0.01;yvlast=2000;zwlast=1000
 
grdpwr(y,ny,yvlast,1.1)     
grdpwr(z,nz,zwlast,1.1)     
 
CARTES=F
 
    GROUP 7. Variables stored, solved & named
SOLUTN(P1,Y,Y,Y,P,P,P); OUTPUT(P1,Y,Y,Y,Y,Y,Y)
SOLUTN(V1,Y,Y,Y,P,P,P); OUTPUT(V1,Y,Y,Y,Y,Y,Y)
SOLUTN(W1,Y,Y,Y,P,P,P); OUTPUT(W1,Y,Y,Y,Y,Y,Y)

store(tem1)
solutn(tem1,y,y,y,p,p,p);output(tem1,y,y,y,y,y,y)
solve(oil)
store(enut,len1)

turmod(kemodl)
rho1=1000.
          ** density difference between water and oil
real(dendif)
dendif=0.03
          ** solve an equation for residence time
solve(rest)
terms(rest,n,y,n,p,p,p)
patch(elapse,phasem,1,nx,1,ny,1,nz,1,lstep)
coval(elapse,rest,fixflu,1.0)

enul=1.e-6     
  
patch(ini1,inival,1,nx,1,ny,1,35,1,lstep)
(initial of tem1 at ini1 is 4.0)
 
patch(ini2,inival,1,nx,1,ny,35+1,nz,1,lstep)
(initial of tem1 at ini2 is 10.0)

        ** buoyancy force is a linear function of both 
           oil concentration and temperature at the warm-water layer
patch(upper,phasem,1,nx,1,ny,35+1,nz,1,lstep)
(source of w1 at upper is 9.81*(oil*:dendif:+1000.*1.e-4*(tem1-10.)$
))                                                                  
 
patch(lower,phasem,1,nx,1,ny,1,35,1,lstep)
(source of w1 at lower is 9.81*:dendif:*oil)

       ** Sine curve profile at the lower surface 
patch(leak,low,1,1,1,12,1,1,1,1) ! with density = 1000 and 
                                            ! maximum velocity 0.03
(source of p1 at leak is 1000*0.03*sin((532-yg)/338.7) with fixflu)
(source of w1 at leak is 0.03*sin((532-yg)/338.7) with onlyms)
coval(leak,oil,onlyms,1.e-3)
coval(leak,ke,onlyms,0.01**2)
coval(leak,ep,onlyms,fiinit(ep))
coval(leak,tem1,onlyms,4.0)

patch(bottom,cell,1,1,12+1,ny,1,1,1,1)
coval(bottom,p1,1.e3,0.0)
coval(bottom,oil,0.0,0.0)
coval(bottom,tem1,0.0,4.0)

patch(north1,cell,1,1,ny,ny,1,35,1,lstep)
coval(north1,p1,1.e3,0.0)
coval(north1,tem1,onlyms,4.0)
 
patch(north2,cell,1,1,ny,ny,35+1,nz,1,lstep)
coval(north2,p1,1.e3,.0)
coval(north2,tem1,onlyms,10.)

(stored var logo is log10(oil))

relax(p1,linrlx,0.2)
relax(oil,linrlx,0.5)
relax(v1,falsdt,0.1)
relax(w1,falsdt,0.1)
relax(tem1,falsdt,10.)
fiinit(ke)=0.01**2
fiinit(ep)=0.164*fiinit(ke)**1.5/10.0
relax(ke,linrlx,0.5)
relax(ep,linrlx,0.5)
 
 
lsweep=2000
 
tstswp=-1
 
idispa=50;csg1=t
 
izmon=25
  nyprin=1;nzprin=1
  iyprl=5;izprl=5
LIBREF=401  
STOP