** MB-FGE Test: Gas leak calculations.
  **************************************************************
  DISPLAY
   In this example we have a gas leak from a small breach in
   a pipe, within a large equipment room.  To "capture"  the
   scale of the  gas leak (of  the order of  1mm) within the
   large room (characteristic length 1m) is difficult to  do
   with a single  grid.  Here  a sequence of  4 successively
   embedded grids are used to resolve both the leak aperture
   itself and the resulting jet.
 
   The  gas  is  heavier  than  air,  and the model includes
   bouyancy effects.  A fixed turbulent viscosity based upon
   experimental results for the round turbulent jet has been
   used.
 
   NOTE !!! The size of F-array should be at least 2000000.
 
   This Q1 file includes Photon Use information.
   ---------------------------------------------------------
  ENDDIS
L(PAUSE
  **************************************************************
  PHOTON USE
   p
 
 
 
 
   use patgeo
   vi x
   set con scale range on
   set vec av off
   set vec ref 10
   mvec 1 i 1 sh
   mvec 2 i 1 sh
   pause
   cl
   use patgeo
   con cgas i 1 fi
   0 0.25
   0.001
   pause
   vi 1 1 1
   mvec 1 j 1 sh
  ENDUSE
  **************************************************************
TEXT(FGEM, Gas Leak Example
TITLE
INTEGER(III,JJ,NXS,NYS,NZS,ICGAS)
INTEGER(NXA,NXB,NXC,NYA,NYB,NYC,NZA,NZB,NZC)
REAL(PI,FSCA,MGAS,VGAS,RHOG,MASGL,MASGH,MOMFLX,DFL,DPO,DPI,DUM1)
REAL(WIND,RHOA,VAIR,LFL,LGAP,LCL,ZCOMP,ZDOM,GRAVTY)
REAL(ZRECP,ZTURB,XPIPE,YPIPE,ZPIPE)
REAL(XMAX,YMAX,ZMAX,XA1,XA4,XB1,XB4,YA1,YA4,YB1,YB4,ZA1)
REAL(ZA4,ZB1,ZB4,XX1,XX4,YY1,YY4,ZZ1,ZZ4,DELX,DELY,DELZ,GAREA)
PI  = 3.14159;  GRAVTY= 9.81;  RHOG=2.5;  RHOA= 1.17
  ** 1kg/min release rate
MGAS= 1.0/60;   VAIR  = 0.2
  ** Grid Parameters
XA4 = 0.;   XB4=1.0;  XMAX=XB4
YA4 = 0.;   YB4=2.0; YMAX=YB4
ZA4 = 0.;   ZB4=3.0;   ZMAX=ZB4
BFC=T
MESG( Generating Grid 1 (Coarse)
XX1= XA4; XX4= XB4; YY1= YA4; YY4= YB4; ZZ1= ZA4; ZZ4= ZB4
NXS= 12;  NYS= 12;  NZS= 18
GSET(D,NXS,NYS,NZS,XMAX,YMAX,ZMAX)
CALL GNGRD2
DUMPC(MBGR1)
MESG( Generating Grid 2 (Fine)
XX1= XC(1,1,1);   XX4= XC(4,1,1)
YY1= YC(1,1,1);   YY4= YC(1,9,1)
ZZ1= ZC(1,1,7);   ZZ4= ZC(1,1,15)
NXS= 12; NYS= 16; NZS= 24
GSET(D,NXS,NYS,NZS,XMAX,YMAX,ZMAX)
CALL GNGRD2
DUMPC(MBGR2)
MESG( Generating Grid 3 (Finer)
XX1= XC(1,1,1);   XX4= XC(4,1,1)
YY1= YC(1,8,1);   YY4= YC(1,16,1)
ZZ1= ZC(1,1,2);   ZZ4= ZC(1,1,11)
NXS= 12; NYS= 16; NZS= 18
GSET(D,NXS,NYS,NZS,XMAX,YMAX,ZMAX)
CALL GNGRD2
DUMPC(MBGR3)
MESG( Generating Grid 4 (Finest)
XX1= XC(1,1,1);   XX4= XC(4,1,1)
YY1= YC(1,8,1);   YY4= YC(1,16,1)
ZZ1= ZC(1,1,2);   ZZ4= ZC(1,1,11)
NXS= 12; NYS= 16; NZS= 18
GSET(D,NXS,NYS,NZS,XMAX,YMAX,ZMAX)
CALL GNGRD2
DELX= XC(10,1,1)-XC(9,1,1)
DELY= YC(1,10,1)-YC(1,9,1)
DELZ= ZC(1,1,10)-ZC(1,1,9)
MESG( Finest cell size is :delx: x :dely: x :delz:
DUMPC(MBGR4)
NUMBLK= 4; READCO(MBGR+)
XULAST=XMAX; YVLAST=YMAX; ZWLAST=ZMAX
MBLINK(2,IN,1)
MBLINK(3,IN,2)
MBLINK(4,IN,3)
STORE(DEN1,VPOR); NAME(C4)=CGAS; ICGAS=19
SOLVE(P1,U1,V1,W1,CGAS)
L($F150)
  * properties
RHO1=RECSCAL; ENUL=1.589e-05
RHO1A=1./RHOA; RHO1B=1./RHOG - 1./RHOA; RHO1C=ICGAS
DPI=2*DELZ; GAREA=3*DELX*DPI; VGAS=MGAS/GAREA/RHOG
MESG( vgas=:vgas: vair=:vair:
  ** Turb viscosity based on turb round jet
ENUT=0.008*VGAS*DPI
  *boundary conditions
MPATCH(1,WINDIN,LOW,1,12,1,12,1,1,1,1)
COVAL(WINDIN,P1,FIXFLU,RHOA*VAIR)
COVAL(WINDIN,UC1,ONLYMS,0.0)
COVAL(WINDIN,VC1,ONLYMS,0.0)
COVAL(WINDIN,WC1,ONLYMS,VAIR)
COVAL(WINDIN,CGAS,ONLYMS,0.0)
CONPOR(PBLOK,0.0,CELL,1,3,  15,15, 70,71)
MPATCH(4,GASINN,NORTH,1,3,  14,14,  7,8, 1,LSTEP)
COVAL(GASINN,P1,FIXFLU,RHOG*VGAS)
COVAL(GASINN,UC1,ONLYMS,0.0)
COVAL(GASINN,VC1,ONLYMS,-VGAS)
COVAL(GASINN,WC1,ONLYMS,0.0)
COVAL(GASINN,CGAS,ONLYMS,1.0)
MPATCH(1,EXIT,HIGH,1,12,1,12,18,18,1,1)
COVAL(EXIT,P1,FIXP,0.0)
COVAL(EXIT,UC1,ONLYMS,0.0)
COVAL(EXIT,VC1,ONLYMS,0.0)
COVAL(EXIT,WC1,ONLYMS,0.0)
COVAL(EXIT,CGAS,ONLYMS,0.0)
MPATCH(1,NBWALL,NWALL,1,12,12,12,1,18,1,1)
MPATCH(1,SBWALL,SWALL,1,12,1,1,1,18,1,1)
MPATCH(2,S2BWALL,SWALL,1,15,1,1,1,24,1,1)
MPATCH(1,EBWALL,EWALL,12,12,1,12,1,18,1,1)
COVAL(NBWALL,UC1,1.0,0.0)
COVAL(NBWALL,WC1,1.0,0.0)
COVAL(SBWALL,UC1,1.0,0.0)
COVAL(SBWALL,WC1,1.0,0.0)
COVAL(S2BWALL,UC1,1.0,0.0)
COVAL(S2BWALL,WC1,1.0,0.0)
COVAL(EBWALL,VC1,1.0,0.0)
COVAL(EBWALL,WC1,1.0,0.0)
BUOYB= -9.81;  BUOYE= RHOA
FIINIT(WC1)=VAIR; FIINIT(DEN1)=RHOA; FIINIT(CGAS)=0.0
LITER(P1)=100; OUTPUT(CGAS,Y,Y,Y,N,Y,N)
TSTSWP=-1; LSWEEP=500
SELREF=T;  RESFAC=1.0E-3
VARMAX(CGAS)=1.0; VARMIN(CGAS)=0.0
REAL(RLXTIM); RLXTIM=0.3
MESG( False time step relaxation = :rlxtim:
RELAX(P1,LINRLX,0.3);     RELAX(UC1,FALSDT,RLXTIM)
RELAX(VC1,FALSDT,RLXTIM); RELAX(WC1,FALSDT,RLXTIM)
  * monitor cell is at the jet footprint on floor
IXMON=1;IYMON=1;IZMON=27
ENDMAIN ! 
SUBROUTINE GNGRD2
  ** Subroutine to simplify grid embedding
 GSET(P,P1 ,XX1,0.,ZZ1); GSET(P,P4 ,XX1,0.,ZZ4)
 GSET(P,P7 ,XX4,0.,ZZ4); GSET(P,P10,XX4,0.,ZZ1)
 GSET(L,L1 ,P1 ,P4 ,NZS,1.); GSET(L,L2 ,P4 ,P7 ,NXS,1.)
 GSET(L,L3 ,P7 ,P10,NZS,1.); GSET(L,L4 ,P10,P1 ,NXS,1.)
 GSET(F,XFRAM,P1,-,P4,-,P7,-,P10,-)
 GSET(M,XFRAM,+K+I,1,1,1,TRANS)
 GSET(C,J1,F,J1,+,0.,YY1,0.)
 III=NYS+1
 GSET(C,J:III:,F,J1,+,0.,YY4-YY1,0.,INC,1.0)
 GSET(I)
ENDSUB