home *** CD-ROM | disk | FTP | other *** search
/ Hall of Fame / HallofFameCDROM.cdr / oilfield / spe-26.lzh / FLOOD.SRC < prev    next >
Text File  |  1987-09-15  |  33KB  |  561 lines
  1. SOURCE
  2. PRECISION= 14
  3. AUTODEF=ON
  4. OPTION BASE=0
  5. ERL=OFF
  6. ERRORMODE=LOCAL
  7. RESUME=LINE
  8. FORMODE=BB
  9. PRINTMODE=GW
  10. SCOPE=ON
  11. PROCS=14
  12. INTEGER: NBuc,NLay,NStr,NMax,ITyp,IDevice
  13. REAL#: VolP,V,SatOi
  14. REAL#: SatOr,SatWi,SatIW,Viso,Visw,FvfO,FvfW,PrmrOi
  15. REAL#: PrmrWr,PrmXpO,PrmXpW,Tim
  16. REAL# ARRAY(11): InjW,DelTim,TimMx
  17. REAL# ARRAY(101): WtrIQ,OilPQ,WtrPQ
  18. REAL# ARRAY(101): PrsQ
  19. REAL# ARRAY(101): RteOQ,RteWQ
  20. REAL# ARRAY(101): PrmThk,Thk,Ara,StrmFn
  21. REAL# ARRAY(1001): Vol,RteCon,WtrI
  22. STRING: FileIn$[13],FileOut$[13]
  23. STRING: Out$[?],Ans$[?]
  24.  
  25. PROCEDURE: DataInput
  26. END PROCEDURE
  27.  
  28. PROCEDURE: Buckley
  29. END PROCEDURE
  30.  
  31. PROCEDURE: Layer
  32. END PROCEDURE
  33.  
  34. PROCEDURE: Stream
  35. END PROCEDURE
  36.  
  37. PROCEDURE: Calculate
  38. END PROCEDURE
  39.  
  40. REAL FUNCTION: Phi
  41.    REAL ARG: Z
  42. END FUNCTION
  43.  
  44. REAL FUNCTION: PhiInv
  45.    REAL ARG: Z
  46. END FUNCTION
  47.  
  48. PROCEDURE: DataOutput
  49. END PROCEDURE
  50.  
  51. PROCEDURE: OutPut()
  52.    STRING ARG: Out$
  53.    INTEGER ARG: IDevice/OPT=0
  54.    INTEGER ARG: IControl/OPT=0
  55. END PROCEDURE
  56.  
  57.  
  58. PROCEDURE: DataInput
  59.    EXTERNAL: NBuc,NLay,NStr,ITyp,VolP,V,SatOi,SatOr
  60.    EXTERNAL: SatWi,SatIW,Viso,Visw,FvfO,FvfW,PrmrOi,PrmrWr
  61.    EXTERNAL: PrmXpO,PrmXpW,InjW(),DelTim(),TimMx(),FileIn$,FileOut$,IDevice
  62.    INTEGER: I
  63.    REAL#: Arg
  64.    STRING: Data$[?],Out$[?]
  65.  
  66.    PROCEDURE: GetData()
  67.       REAL ARG: X/VAR
  68.    END PROCEDURE
  69.  
  70.    PROCEDURE: GetData
  71.       EXTERNAL: Data$
  72.       INTEGER: I,J
  73.          10 '----------------------ROUTINE GETDATA----------------------------------------
  74.          20 '.....this is the routine that actually reads the data from the input file
  75.          30 IF LEN(Data$)>0 GOTO 210
  76.          40 '...read a line from the input data file
  77.          50 IF EOF(1) THEN X=-999: GOTO 260
  78.          60 LINE INPUT #1 Data$
  79.          70 '...now sort thru Data$ and delete all characters after the first non-numeric or delimiter character
  80.          80 I=0
  81.          90 I=I+1
  82.         100 IF I>LEN(Data$) GOTO 180
  83.         110 J=ASC(MID$(Data$,I,1)): IF J>47 AND J<58 GOTO 90:'...numeric data are ok
  84.         120 IF J=45 OR J=46 OR J=69 OR J=101 GOTO 90:'decimal point, E for exponential, minus sign are ok
  85.         130 IF J=44 OR J=32 OR J=58 OR J=59 OR J=9 GOTO 150:'delimiters are comma, space, semicolon, colon, or tab
  86.         140 IF I=1 GOTO 40: ELSE Data$=LEFT$(Data$,I-1): GOTO 180
  87.         150 MID$(Data$,I,1)=",": IF I=1 THEN Data$=RIGHT$(Data$,LEN(Data$)-1): GOTO 100:'set all delimiters to commas
  88.         160 IF MID$(Data$,I-1,1)="," THEN Data$=LEFT$(Data$,I-1)+RIGHT$(Data$,LEN(Data$)-I): GOTO 100
  89.         170 GOTO 90
  90.         180 '...now remove trailing delimiters, if any
  91.         190 IF LEN(Data$)=0 GOTO 40
  92.         200 J=ASC(RIGHT$(Data$,1)): IF J=44 THEN Data$=LEFT$(Data$,LEN(Data$)-1): GOTO 190
  93.         210 '...now get data
  94.         220 I=0
  95.         230 I=I+1: IF I>LEN(Data$) THEN X=VAL(Data$): Data$="": GOTO 260
  96.         240 IF MID$(Data$,I,1)="," THEN X=VAL(LEFT$(Data$,I-1)): Data$=RIGHT$(Data$,LEN(Data$)-I): GOTO 260
  97.         250 GOTO 230
  98.         260 EXIT
  99.    END PROCEDURE
  100.       10 '-----------------ROUTINE DATAINPUT-------------------------------------------
  101.       20 '.....this routine does the data entry 
  102.       30 LINE INPUT "ENTER  INPUT FILENAME: "; FileIn$: IF FileIn$="" THEN FileIn$="FLOOD.INP"
  103.       40 LINE INPUT "ENTER OUTPUT FILENAME: "; FileOut$:IF FileOut$="" THEN FileOut$="FLOOD.OUT"
  104.       50 OPEN "I",1,FileIn$
  105.       60 GetData (VolP):      '...pore volume
  106.       70 GetData (NBuc):      '...number of buckley-leverett increments (3 for pistonlike displacement)
  107.       80 GetData (NLay):      '...number of layers (0 for pseudo-relative permeabilities)
  108.       90 GetData (NStr):      '...number of streamlines
  109.      100 GetData (V):         '...coefficient of permeability variation
  110.      110 GetData (ITyp):      '...pattern type (5,7,or 9 spot, -5 or -7 for isolated patterns, 0 for no areal effects)
  111.      120 GetData (SatOi):     '...initial oil saturation
  112.      130 GetData (SatWi):     '...initial water saturation
  113.      140 GetData (SatOr):     '...residual oil saturation
  114.      150 GetData (SatIW):     '...irreducible water saturation
  115.      160 GetData (Viso):      '...oil viscosity
  116.      170 GetData (Visw):      '...water viscosity
  117.      180 GetData (PrmrOi):    '...oil relative perm at irreducible water saturation
  118.      190 GetData (PrmrWr):    '...water relative perm at residual oil saturation
  119.      200 GetData (PrmXpO):    '...oil relative perm exponent,   Prmro=Prmroi*((SatO-SatOr)/(1-SatIW-SatOr))^PrmXpo
  120.      210 GetData (PrmXpW):    '...water relative perm exponent, Prmrw=PrmrWr*((SatW-SatIW)/(1-SatOr-SatIW))^PrmXpW
  121.      220 GetData (FvfO):      '...oil formation volume factor
  122.      230 GetData (FvfW):      '...water formation volume factor
  123.      240 FOR I=1 TO 10: InjW(I)=0: DelTim(I)=0: TimMx(I)=0: NEXT I: I=0
  124.      250 I=I+1
  125.      260 GetData (InjW(I)):   '...water injection rate.  units should be consistent with units of pore volume and time
  126.      270 GetData (DelTim(I)): '...time step
  127.      280 GetData (TimMx(I)):  '...maximum time at time step DelTim(i)
  128.      290 IF TimMx(I)>0 AND I<10 GOTO 250
  129.      300 CLOSE #1
  130.      310 '.....check data for consistency.....
  131.      320 IF TimMx(1)<0 THEN PRINT "INSUFFICIENT DATA IN INPUT FILE FOR ANALYSIS.  EXECUTION DELETED. <<<<<": STOP
  132.      330 IF VolP<0 THEN VolP=1:Out$="pore volume":Arg=VolP: GOSUB 600
  133.      340 IF NLay=1 GOTO 370: '...if only one layer, V doesn't matter
  134.      350 IF V<=0 THEN V=1.E-5:Out$="coefficient of permeability variation":Arg=V:GOSUB 600
  135.      360 IF V>=1 THEN V=.99999:Out$="coefficient of permeability variation":Arg=V:GOSUB 600
  136.      370 IF NBuc<3 THEN NBuc=3:Out$="number of buckley-leverett increments":Arg=NBuc:GOSUB 600
  137.      380 IF NBuc>100 THEN NBuc=100:Out$="number of buckley-leverett increments":Arg=NBuc:GOSUB 600
  138.      390 IF NStr<1 THEN NStr=1:Out$="number of streamlines":Arg=NStr:GOSUB 600
  139.      400 IF NStr>20 THEN NStr=20:Out$="number of streamlines":Arg=NStr:GOSUB 600
  140.      410 I=100: IF I>500/NStr THEN I=500/NStr
  141.      420 IF NLay>I THEN NLay=I:Out$="number of layers":Arg=NLay:GOSUB 600
  142.      430 IF SatOi<0 THEN SatOi=0:Out$="initial oil saturation":Arg=SatOi:GOSUB 600
  143.      440 IF SatOi>1 THEN SatOi=1:Out$="initial oil saturation":Arg=SatOi:GOSUB 600
  144.      450 IF SatWi<0 THEN SatWi=0:Out$="initial water saturation":Arg=SatOi:GOSUB 600
  145.      460 IF SatWi>1-SatOi THEN SatWi=1-SatOi:Out$="initial water saturation":Arg=SatOi:GOSUB 600
  146.      470 IF SatOr<0 THEN SatOr=0:Out$="residual oil saturation":Arg=SatOr:GOSUB 600
  147.      480 IF SatOr>1 THEN SatOr=1:Out$="residual oil saturation":Arg=SatOr:GOSUB 600
  148.      490 IF SatIW<0 THEN SatIW=0:Out$="irreducible water saturation":Arg=SatIW:GOSUB 600
  149.      500 IF SatIW>1 THEN SatIW=1:Out$="irreducible water saturation":Arg=SatIW:GOSUB 600
  150.      510 IF Viso<0 THEN Viso=1:Out$="oil viscosity":Arg=Viso:GOSUB 600
  151.      520 IF Visw<0 THEN Visw=1:Out$="water viscosity":Arg=Visw:GOSUB 600
  152.      530 IF PrmrOi<0 THEN PrmrOi=1:Out$="oil relative perm at irreducible water saturation":Arg=PrmrOi:GOSUB 600
  153.      540 IF PrmrWr<0 THEN PrmrWr=1:Out$="water relative perm at residual oil saturation":Arg=PrmrWr:GOSUB 600
  154.      550 IF FvfO<0 THEN FvfO=1:Out$="oil formation volume factor":Arg=FvfO:GOSUB 600
  155.      560 IF FvfW<0 THEN FvfW=1:Out$="water formation volume factor":Arg=FvfW:GOSUB 600
  156.      570 '.....now open file for output......
  157.      580 OPEN "O",1,FileOut$: IDevice=1
  158.      590 EXIT
  159.      600 '.....error encountered................................
  160.      610 '.....print error message and reset data.................
  161.      620 PRINT ".....ERROR encountered with ";Out$
  162.      630 PRINT SPACE$(5);Out$;" set equal to ";Arg
  163.      640 RETURN
  164. END PROCEDURE
  165.  
  166. PROCEDURE: Buckley
  167.    EXTERNAL: NBuc,NLay,V,SatOi,SatOr,SatWi,SatIW,Viso
  168.    EXTERNAL: Visw,PrmrOi,PrmrWr,PrmXpO,PrmXpW,WtrIQ(),OilPQ(),RteOQ()
  169.    EXTERNAL: RteWQ(),WtrPQ(),PrsQ(),Phi,PhiInv
  170.    INTEGER: I,N1,N2,I1
  171.    INTEGER: MaxIter,IFillup
  172.    REAL#: SatWMn,SatWMx,PrmrO,PrmrW,FrcWI,DFrcW,Mbr
  173.    REAL#: FrcW,DFrcW0,Mbr0,SatW,F,SatWf,FrcWf,DFrcWf
  174.    REAL#: Mbrf,Tolerance,D2FrcW,SatWInf
  175.  
  176.    PROCEDURE: Flow()
  177.       REAL ARG: SatW/VAR,PrmrO/VAR,PrmrW/VAR,FrcW/VAR,DFrcW/VAR,Mbr/VAR,D2FrcW/VAR
  178.    END PROCEDURE
  179.  
  180.    PROCEDURE: Flow
  181.       EXTERNAL: NLay,V,SatOr,SatIW,Viso,Visw,PrmrOi,PrmrWr
  182.       EXTERNAL: PrmXpO,PrmXpW,Phi,PhiInv
  183.       REAL#: DelSatW,PrmrWp,PrmrWm,PrmrOp,PrmrOm
  184.       REAL#: SatO
  185.  
  186.       REAL FUNCTION: FPrmrO
  187.          REAL ARG: SatO
  188.       END FUNCTION
  189.  
  190.       REAL FUNCTION: FPrmrW
  191.          REAL ARG: SatW
  192.       END FUNCTION
  193.  
  194.       REAL FUNCTION: FPrmrO
  195.          EXTERNAL: NLay,V,SatOr,SatIW,PrmrOi,PrmXpO,Phi,PhiInv
  196.          REAL#: SatQ,Sdv,Arg
  197.             10 '-------------------------FUNCTION FPRMRO------------------------------
  198.             20 '.....function to compute oil relative permeability
  199.             30 '..... PrmrO=Prmroi*((SatO-SatOr)/(1-SatIW-SatOr))^PrmXpO
  200.             40 IF SatO<=SatOr THEN RESULT=0: EXIT
  201.             50 IF NLay<=0 GOTO 90: '... Hiatt-Hearn pseudo-rel perms
  202.             60 SatQ=(SatO-SatOr)/(1-SatOr-SatIW): IF SatQ>1 THEN RESULT=PrmrOi: EXIT
  203.             70 RESULT=PrmrOi*SatQ^PrmXpO: EXIT
  204.             80 '.....the following expression calculates pseudo-rel perms using the Hiatt-Hearn model
  205.             90 SatQ=(1-SatO-SatIW)/(1-SatOr-SatIW): Sdv=-LOG(1-V):Arg=PhiInv(SatQ)+Sdv: RESULT=PrmrOi*Phi(-Arg): EXIT
  206.       END FUNCTION
  207.  
  208.       REAL FUNCTION: FPrmrW
  209.          EXTERNAL: NLay,V,SatOr,SatIW,PrmrWr,PrmXpW,Phi,PhiInv
  210.          REAL#: SatQ,Sdv,Arg
  211.             10 '-------------------------FUNCTION FPRMRW------------------------------
  212.             20 '.....function to compute water relative permeability
  213.             30 '.....PrmrW=PrmrWr*((SatW-SatIW)/(1-SatOr-SatIW))^PrmXpW
  214.             40 IF SatW<=SatIW THEN RESULT=0: EXIT
  215.             50 IF NLay<=0 AND SatW<1-SatOr GOTO 80:'...go do pseudo rel perms
  216.             60 SatQ=(SatW-SatIW)/(1-SatOr-SatIW): RESULT=PrmrWr*SatQ^PrmXpW: EXIT
  217.             70 '.....the folowing equation computes pseudo-rel perms with the Hiatt-Hearn method
  218.             80 SatQ=(SatW-SatIW)/(1-SatIW-SatOr): Sdv=-LOG(1-V): Arg=PhiInv(SatQ)+Sdv: RESULT=PrmrWr*Phi(Arg): EXIT
  219.       END FUNCTION
  220.          10 '-------------------------PROCEDURE FLOW-------------------------------
  221.          20 'routine to calculate fractional flow curves for a liquid filled system
  222.          30 DelSatW=1.E-5: SatW=SatW+DelSatW: SatO=1-SatW:PrmrOp=FPrmrO(SatO): PrmrWp=FPrmrW(SatW)
  223.          40 SatW=SatW-2*DelSatW:SatO=1-SatW: PrmrOm=FPrmrO(SatO): PrmrWm=FPrmrW(SatW): SatW=SatW+DelSatW
  224.          50 SatO=1-SatW: PrmrO=FPrmrO(SatO): PrmrW=FPrmrW(SatW)
  225.          60 IF SatW<=SatIW OR PrmrW<=0 THEN FrcW=0: DFrcW=0: GOTO 100
  226.          70 IF SatW>=1-SatOr OR PrmrO<=0 THEN FrcW=1: DFrcW=0: GOTO 100
  227.          80 FrcW=1/(1+PrmrO/PrmrW*Visw/Viso): DFrcW=FrcW*FrcW*Visw/Viso*PrmrO/PrmrW/2/DelSatW*((PrmrWp-PrmrWm)/PrmrW-(PrmrOp-PrmrOm)/PrmrO)
  228.          90 D2FrcW=0: IF FrcW>0 THEN D2FrcW=2/FrcW*DFrcW^2-FrcW^2*Visw/Viso*PrmrO/PrmrW*((PrmrOp-2*PrmrO+PrmrOm)/PrmrO-(PrmrWp-2*PrmrW+PrmrWm)/PrmrW+2*(PrmrWp-PrmrWm)/PrmrW*((PrmrWp-PrmrWm)/PrmrW-(PrmrOp-PrmrOm)/PrmrO)/4)/DelSatW^2
  229.         100 Mbr=PrmrO/Viso+PrmrW/Visw: EXIT
  230.    END PROCEDURE
  231.       10 '-------------------------PROCEDURE BUCKLEY------------------------------
  232.       20 ' Routine to compute 1-D Buckley Leverett displacement
  233.       30 ' Bisection is used to determine the frontal saturation
  234.       40 ' Arrays are set up for oil prod(DN), water prod(DW), pressure drop(DP),
  235.       50 '     as a function of water injection (DI).  The first element in the
  236.       60 '     arrays is for fillup, the second is for oil bank breakthrough,
  237.       70 '     and then use NBUC-2 equal saturation increments between the frontal
  238.       80 '     saturation and ultimate recovery.
  239.       90 IF NBuc>3 GOTO 160
  240.      100 '.....calculate displacement for pistonlike front
  241.      110 SatWMn=SatWi: SatWMx=1-SatOr: IF NBuc<3 THEN NBuc=3
  242.      120 Flow(SatWMn,PrmrO,PrmrW,FrcWI,DFrcW,Mbr,D2FrcW): Flow(SatWMx,PrmrO,PrmrW,FrcW,DFrcW0,Mbr0,D2FrcW)
  243.      130 WtrIQ(1)=1-SatOi-SatWi: OilPQ(1)=0.: WtrPQ(1)=0.: PrsQ(0)=1/Mbr: RteOQ(1)=1-FrcWI: RteWQ(1)=FrcWI
  244.      140 WtrIQ(2)=(1-SatWi-SatIW)/(1-FrcWI): OilPQ(2)=SatOi-SatOr: WtrPQ(2)=FrcWI*(WtrIQ(2)-WtrIQ(1)): PrsQ(2)=1/Mbr0: PrsQ(1)=PrsQ(0)+(PrsQ(2)-PrsQ(0))*WtrIQ(1)/WtrIQ(2): RteOQ(2)=0: RteWQ(2)=1
  245.      150 WtrIQ(NBuc)=1.E25: OilPQ(NBuc)=SatOi-SatOr: WtrPQ(NBuc)=WtrIQ(NBuc)-OilPQ(NBuc): PrsQ(NBuc)=1./Mbr0: RteOQ(NBuc)=0: RteWQ(NBuc)=1: EXIT
  246.      160 '.....calculate displacement for Buckley-Leverett front
  247.      170 MaxIter=20: Tolerance=1.E-6 : IFillup=0
  248.      180 SatWMn=SatWi: SatWMx=1-SatOr
  249.      190 Flow(SatWMn,PrmrO,PrmrW,FrcWI,DFrcW,Mbr,D2FrcW): Flow(SatWMx,PrmrO,PrmrW,FrcW,DFrcW0,Mbr0,D2FrcW)
  250.      200 WtrIQ(1)=1-SatOi-SatWi: OilPQ(1)=0: WtrPQ(1)=0: PrsQ(0)=1/Mbr: RteOQ(1)=1-FrcWI: RteWQ(1)=FrcWI
  251.      210 WtrIQ(NBuc)=1.E25: OilPQ(NBuc)=SatOi-SatOr: WtrPQ(NBuc)=WtrIQ(NBuc)-OilPQ(NBuc): PrsQ(NBuc)=1./Mbr0: RteOQ(NBuc)=0: RteWQ(NBuc)=1: I=0
  252.      220 '.....find the maximum water sat'n at which an oil bank can form.................................................
  253.      230 SatWMn=SatIW: SatWMx=1-SatOr: I=0
  254.      240 I=I+1: SatW=(SatWMn+SatWMx)/2: Flow(SatW,PrmrO,PrmrW,FrcW,DFrcW,Mbr,D2FrcW)
  255.      250 F=D2FrcW:IF ABS(F)<=Tolerance THEN SatWInf=SatW: GOTO 280
  256.      260 IF F>0 THEN SatWMn=SatW ELSE SatWMx=SatW
  257.      270 IF I<MaxIter GOTO 240 ELSE SatWInf=(SatWMx+SatWMn)/2
  258.      280 IF SatWi>=SatWInf GOTO 370
  259.      290 '... now find frontal saturation.................................................................................
  260.      300 SatWMn=SatWInf: SatWMx=1-SatOr: I=0
  261.      310 I=I+1: SatW=(SatWMn+SatWMx)/2: Flow(SatW,PrmrO,PrmrW,FrcW,DFrcW,Mbr,D2FrcW)
  262.      320 F=DFrcW*(SatW-SatWi)-(FrcW-FrcWI):IF ABS(F)<=Tolerance THEN SatWf=SatW: GOTO 350
  263.      330 IF F>0 THEN SatWMn=SatW ELSE SatWMx=SatW
  264.      340 IF I<MaxIter GOTO 310 ELSE SatWf=(SatWMx+SatWMn)/2
  265.      350 Flow(SatWf,PrmrO,PrmrW,FrcWf,DFrcWf,Mbrf,D2FrcW)
  266.      360 IF WtrIQ(1)*DFrcWf<=1 GOTO 440
  267.      370 '.....calculated breakthrough occurs prior to fillup -- model not applicable.....................................
  268.      380 PRINT "INSUFFICIENT OIL SATURATION FOR OIL BANK FORMATION.   EXECUTION DELETED.  <<<<<<<<<<":STOP
  269.      390 I=0: SatWMn=SatWi: SatWMx=1-SatOi: IFillup=1
  270.      400 I=I+1: SatW=(SatWMn+SatWMx)/2: Flow(SatW,PrmrO,PrmrW,FrcW,DFrcW,Mbr,D2FrcW)
  271.      410 F=SatW-SatWi-(1-SatWi-SatOi)*FrcW:IF ABS(F)<=Tolerance THEN SatWf=SatW: GOTO 420
  272.      420 IF F>0 THEN SatWMn=SatW: ELSE SatWMx=SatW
  273.      430 IF I<MaxIter GOTO 400 ELSE SatWf=(SatWMn+SatWMx)/2
  274.      440 '.....now set up table lookups for inj, oil prod, wtr prod, delta-press, etc.....................................
  275.      450 '.....do table at equal saturation increments between the frontal saturation and residual conditions
  276.      460 SatWMn=SatWf: Flow(SatWf,PrmrO,PrmrW,FrcWf,DFrcWf,Mbrf,D2FrcW): SatWMx=1-SatOr: N1=NBuc-1: N2=N1-1
  277.      470 FOR I1=2 TO N1: I=N1-I1+2: SatW=SatWMn+(SatWMx-SatWMn)*(I-2)/N2
  278.      480   Flow(SatW,PrmrO,PrmrW,FrcW,DFrcW,Mbr,D2FrcW): PrsQ(I)=(1/Mbr+1/Mbr0)*(DFrcW-DFrcW0)/2
  279.      490   WtrIQ(I)=1/DFrcW: OilPQ(I)=WtrIQ(I)*(1-FrcW)+SatW+SatOi-1: IF OilPQ(I)<0 THEN OilPQ(I)=0 ELSE IF OilPQ(I)>WtrIQ(I) THEN OilPQ(I)=WtrIQ(I)
  280.      500   WtrPQ(I)=WtrIQ(I)*FrcW+SatWi-SatWf: IF WtrPQ(I)<0 THEN WtrPQ(I)=0: ELSE IF WtrPQ(I)>WtrIQ(I)-OilPQ(I) THEN WtrPQ(I)=WtrIQ(I)-OilPQ(I)
  281.      510 DFrcW0=DFrcW: Mbr0=Mbr: RteOQ(I)=1-FrcW: RteWQ(I)=FrcW: NEXT I1
  282.      520 '.....now do presure drop integration............................................................................
  283.      530 FOR I1=2 TO N2: I=N2-I1+2: PrsQ(I)=PrsQ(I+1)+PrsQ(I): NEXT I1
  284.      540 FOR I=2 TO N1: PrsQ(I)=PrsQ(I)*WtrIQ(I): NEXT I
  285.      550 IF IFillup=1 THEN PrsQ(0)=1/Mbrf: RteOQ(1)=RteOQ(2): RteWQ(1)=RteWQ(2)
  286.      560 PrsQ(1)=PrsQ(0)+(PrsQ(2)-PrsQ(0))*WtrIQ(1)/WtrIQ(2)
  287.      570 EXIT
  288. END PROCEDURE
  289.  
  290. PROCEDURE: Layer
  291.    EXTERNAL: NLay,V,PrmThk(),Thk(),Phi,PhiInv
  292.    INTEGER: I
  293.    REAL#: Sdv,P0,Y
  294.    REAL#: Xp,P
  295.       10 '------------------------PROCEDURE LAYER-------------------------------
  296.       20 'Routine to compute layer distribution
  297.       30 IF NLay<=1 THEN PrmThk(1)=1: Thk(1)=1: EXIT
  298.       40 GOTO 150:'...normally use different perm and thickness for each layer
  299.       50 '.........................
  300.       60 '...this section uses equal PrmThk layers
  301.       70 FOR I=1 TO NLay: PrmThk(I)=1/NLay: NEXT I
  302.       80 Sdv=-LOG(1-V): P0=0: FOR I=1 TO NLay-1: Y=I/NLay: Xp=PhiInv(Y)-Sdv: P=Phi(Xp): Thk(I)=P-P0: P0=P: NEXT I
  303.       90 Thk(NLay)=1-P: EXIT
  304.      100 '.........................
  305.      110 '...this section uses equal thickness layers
  306.      120 FOR I=1 TO NLay: Thk(I)=1/NLay: NEXT I
  307.      130 Sdv=-LOG(1-V): P0=0: FOR I=1 TO NLay-1: Y=I/NLay: Xp=PhiInv(Y)-Sdv: P=Phi(Xp): PrmThk(I)=P-P0: P0=P: NEXT I
  308.      140 PrmThk(NLay)=1-P: EXIT
  309.      150 '.........................
  310.      160 '...this section uses thickness varying linearly with the layer number
  311.      170 '...which gives best results for small number of layers
  312.      180 FOR I=1 TO NLay: Thk(I)=(NLay+1-I)/NLay/(NLay+1)*2: NEXT I
  313.      190 Sdv=-LOG(1-V): P0=0: Y=0:FOR I=1 TO NLay-1: Y=Y+Thk(I): Xp=PhiInv(Y)-Sdv: P=Phi(Xp): PrmThk(I)=P-P0: P0=P: NEXT I
  314.      200 PrmThk(NLay)=1-P: EXIT
  315. END PROCEDURE
  316.  
  317. PROCEDURE: Stream
  318.    EXTERNAL: NStr,ITyp,Ara(),StrmFn()
  319.    INTEGER: I,J,K,NStDat
  320.    REAL# ARRAY(?): Strm
  321.    REAL#: Pi,U
  322.    STRING: S$[?]
  323.       10 '---------------------------PROCEDURE STREAM------------------------------
  324.       20 'Routine to compute streamline distribution
  325.       30 DEFINT I-N: DEFDBL A-H,O-Z: NStDat=20: DIM Strm(5*NStDat+1): Pi=4*ATN(1)
  326.       40 DATA .035907,.071925,.108167,.144751,.181799,.219443,.257826,.297108,.337467,.379109,.422276,.467255,.514403,.564170,.617150,.674168,.736455,.806073,.887306,1: '5-SPOT PATTERN
  327.       50 DATA .037201,.074505,.112014,.149838,.188092,.226888,.266366,.306667,.347955,.390417,.434271,.479778,.527259,.577119,.629897,.686317,.747499,.815288,.893521,1: '7-SPOT PATTERN
  328.       60 DATA .025839,.051862,.078262,.105246,.133048,.161943,.192265,.224440,.259037,.296864,.339169,.388145,.448794,.551506,.649352,.729621,.802082,.870117,.935637,1: '9-SPOT PATTERN
  329.       70 DATA .052422,.105216,.158763,.213493,.269853,.328375,.389686,.454549,.523927,.599070,.681657,.774034,.879623,1.003700,1.154983,1.349288,1.619447,2.050128,2.976872,4.251752:'ISOLATED 5-SPOT
  330.       80 DATA .045382,.090990,.137056,.183827,.231575,.280606,.331274,.384005,.439323,.497889,.560576,.628570,.703556,.788064,.886149,1.004921,1.158573,1.381894,1.802689,2.305500: 'ISOLATED 7-SPOT
  331.       90 FOR I=1 TO 5*NStDat: READ S$: Strm(I)=VAL(S$): NEXT I
  332.      100 IF NStr<2 THEN NStr=1: Ara(1)=1: StrmFn(1)=1: GOTO 270
  333.      110 '... first setup streamline to proper flood pattern
  334.      120 IF ITyp= 5 THEN J=       0: GOTO 180
  335.      130 IF ITyp= 7 THEN J=  NStDat: GOTO 180
  336.      140 IF ITyp= 9 THEN J=2*NStDat: GOTO 180
  337.      150 IF ITyp=-5 THEN J=3*NStDat: GOTO 180
  338.      160 IF ITyp=-7 THEN J=4*NStDat: GOTO 180
  339.      170 NStr=1: Ara(1)=1: GOTO 270: '... no pattern effect considered
  340.      180 FOR I=1 TO NStDat: K=I+J: Strm(I)=Strm(K): NEXT I: IF ITyp>0 GOTO 190
  341.      190 FOR I=1 TO NStr: U=I/NStr*NStDat: J=INT(U): IF J>=NStDat THEN J=NStDat-1
  342.      200   IF I=NStr AND ITyp<0 THEN U=(I-0.5)/NStr*NStDat: J=INT(U): IF J>=NStDat THEN J=NStDat-1: 'for unconfined patterns, adjust largest streamline
  343.      210   Ara(I)=Strm(J)+(Strm(J+1)-Strm(J))*(U-J)
  344.      220   IF J=NStDat-1 AND ITyp<0 THEN Ara(I)=(Strm(J)+Strm(1))*(SIN(Pi/NStDat)/SIN(Pi*(1-U/NStDat)))^(-2/(ITyp+1))-Strm(1)*(1-U/NStDat)
  345.      230 NEXT I
  346.      240 IF ITyp<0 THEN Ara(NStr)=Ara(NStr)^2/Ara(NStr-1)
  347.      250 FOR I=1 TO NStr-1: J=NStr-I: Ara(J+1)=Ara(J+1)-Ara(J): NEXT I
  348.      260 FOR I=1 TO NStr: StrmFn(I)=1/NStr: NEXT I
  349.      270 ERASE Strm(): EXIT
  350. END PROCEDURE
  351.  
  352. PROCEDURE: Calculate
  353.    EXTERNAL: NBuc,NLay,NStr,NMax,VolP,FvfO,FvfW,InjW()
  354.    EXTERNAL: DelTim(),TimMx(),WtrIQ(),OilPQ(),RteOQ(),RteWQ(),WtrPQ(),PrsQ()
  355.    EXTERNAL: PrmThk(),Thk(),Ara(),StrmFn(),Vol(),RteCon(),WtrI(),IDevice
  356.    EXTERNAL: OutPut
  357.    INTEGER: I,J,ICount,ITim,L
  358.    REAL#: OilPQj,WtrPQj,FrcW,Tim
  359.    REAL#: CumInj,OilP,WtrP,FacWO,Prs,WtrIEnd,RteW,RteO
  360.    REAL#: FrcO,RteNew,RteAv,WtrPQi,Rte,Rte0,Rte1,Rte2
  361.    REAL#: Rte3
  362.    REAL# ARRAY(?): Rte0j,Rte1j,Rte2j,Rte3j
  363.    STRING: Out$[?],CR$[2]
  364.  
  365.    PROCEDURE: Displace()
  366.       REAL ARG: WInj
  367.       REAL ARG: OilPQj/VAR,WtrPQj/VAR,RteQj/VAR,FrcW/VAR,FrcO/VAR
  368.    END PROCEDURE
  369.  
  370.    PROCEDURE: Displace
  371.       EXTERNAL: NBuc,WtrIQ(),OilPQ(),RteOQ(),RteWQ(),WtrPQ(),PrsQ()
  372.       INTEGER: I
  373.       REAL#: U
  374.          10 '---------------------PROCEDURE DISPLACE---------------------------------
  375.          20 'routine to lookup 1-D displacement results from tables created by BUCKLEY routine
  376.          30 IF WInj<=0 THEN OilPQj=0: WtrPQj=0: FrcW=0: RteQj=1/PrsQ(0): EXIT:'... initial injectivity is based on the permeability and streamfunction
  377.          40 IF NBuc=3 GOTO 320:' for pistonlike displacements
  378.          50 IF WInj>WtrIQ(1) GOTO 110
  379.          60 '.....production before fillup
  380.          70 OilPQj=0: WtrPQj=0: RteQj=1/(PrsQ(1)*WInj/WtrIQ(1))
  381.          80 FrcO=0: FrcW=0: IF WInj>WtrIQ(1) AND WtrPQ(2)>0 THEN FrcW=1/(1+OilPQ(2)/WtrPQ(2))
  382.          90 EXIT
  383.         100 '... this section is for WInj>0
  384.         110 IF WInj>WtrIQ(2) GOTO 190
  385.         120 '... this section is for WInj>0
  386.         130 I=2
  387.         140 U=(WInj-WtrIQ(I-1))/(WtrIQ(I)-WtrIQ(I-1))
  388.         150 OilPQj=OilPQ(I-1)+(OilPQ(I)-OilPQ(I-1))*U: WtrPQj=WtrPQ(I-1)+(WtrPQ(I)-WtrPQ(I-1))*U
  389.         160 RteQj=1/(PrsQ(I-1)+(PrsQ(I)-PrsQ(I-1))*U)
  390.         170 FrcO=RteOQ(1): FrcW=RteWQ(1)
  391.         180 GOTO 310
  392.         190 IF WInj<WtrIQ(NBuc-1) GOTO 250
  393.         200 OilPQj=OilPQ(NBuc)-(OilPQ(NBuc)-OilPQ(NBuc-1))*WtrIQ(NBuc-1)/WInj
  394.         210 WtrPQj=WtrPQ(NBuc-1)+WInj-WtrIQ(NBuc-1)-OilPQj+OilPQ(NBuc-1)
  395.         220 RteQj=1/(PrsQ(NBuc)-(PrsQ(NBuc)-PrsQ(NBuc-1))*WtrIQ(NBuc-1)/WInj)
  396.         230 FrcO=RteOQ(NBuc-1)*(WtrIQ(NBuc-1)/WInj)^2: FrcW=1-FrcO
  397.         240 GOTO 310
  398.         250 I=2
  399.         260 I=I+1: IF WInj>WtrIQ(I) AND I<NBuc GOTO 260
  400.         270 U=(WInj-WtrIQ(I-1))/(WtrIQ(I)-WtrIQ(I-1))
  401.         280 OilPQj=OilPQ(I-1)+(OilPQ(I)-OilPQ(I-1))*U: WtrPQj=WtrPQ(I-1)+(WtrPQ(I)-WtrPQ(I-1))*U
  402.         290 RteQj=1/(PrsQ(I-1)+(PrsQ(I)-PrsQ(I-1))*U)
  403.         300 FrcO=RteOQ(I-1)+(RteOQ(I)-RteOQ(I-1))*U: FrcW=RteWQ(I-1)+(RteWQ(I)-RteWQ(I-1))*U
  404.         310 EXIT
  405.         320 '.........pistonlike front.......................
  406.         330 IF WInj<WtrIQ(1) THEN OilPQj=0: WtrPQj=0: FrcO=0: FrcW=0: RteQj=1/(PrsQ(1)*WInj/WtrIQ(1)): EXIT
  407.         340 IF WInj>WtrIQ(2) THEN OilPQj=OilPQ(2): WtrPQj=WInj-WtrIQ(2)+WtrPQ(2): FrcO=0: FrcW=1: RteQj=1/PrsQ(3): EXIT
  408.         350 U=(WInj-WtrIQ(1))/(WtrIQ(2)-WtrIQ(1))
  409.         360 OilPQj=OilPQ(2)*U: WtrPQj=WtrPQ(2)*U
  410.         370 FrcO=OilPQ(2)/(WtrIQ(2)-WtrIQ(1)): FrcW=WtrPQ(2)/(WtrIQ(2)-WtrIQ(1))
  411.         380 RteQj=1/(PrsQ(1)+(PrsQ(2)-PrsQ(I))*U): EXIT
  412.    END PROCEDURE
  413.       10 '---------------------PROCEDURE CALCULATE--------------------------------
  414.       20 '.....Routine to do calculations
  415.       30 '.....first set up streamtubes
  416.       40 DIM Rte0j(501),Rte1j(501),Rte2j(501),Rte3j(501)
  417.       50 L=0: Rte0=0: FOR I=1 TO NLay: FOR J=1 TO NStr: L=L+1: Vol(L)=Thk(I)*Ara(J)*VolP: RteCon(L)=PrmThk(I)*StrmFn(J): WtrI(L)=0
  418.       60 Displace(WtrI(L),OilPQj,WtrPQj,Rte,FrcW,FrcO):Rte0j(L)=Rte*RteCon(L): Rte0=Rte0+Rte0j(L): NEXT J: NEXT I
  419.       70 '.....now type headings for output
  420.       80 NMax=L: Tim=0: ITim=0: CumInj=0: ICount=1: OilP=0: WtrP=0: FacWO=0: Prs=0: IF Rte0>0 THEN Prs=1/(PrsQ(NBuc)*Rte0)
  421.       90 CR$=CHR$(13)+CHR$(10):Out$=SPACE$(27)+"STREAMTUBE  PERFORMANCE  MODELLING": OutPut(Out$,IDevice)
  422.      100 Out$=SPACE$(27)+STRING$(10,45)+"  "+STRING$(11,45)+"  "+STRING$(9,45)+CR$: OutPut(Out$,IDevice)
  423.      110 Out$=SPACE$(39): PRINT TO Out$,DATE$;SPACE$(3);TIME$
  424.      120 Out$=Out$+CR$: OutPut(Out$,IDevice)
  425.      130 Out$=SPACE$(28)+"Relative": OutPut(Out$,IDevice)
  426.      140 Out$=SPACE$(15)+"Cumulative   Pressure    Cum Oil    Cum Water    Water-Oil": OutPut(Out$,IDevice)
  427.      150 Out$="       Time    Injection      Drop     Production  Production     Ratio": OutPut(Out$,IDevice)
  428.      160 '.......................................................................
  429.      170 '.....relative pressure drop is the pressure drop relative to that at residual oil saturation
  430.      180 '.....water-oil ratio is the instantaneous producing water-oil ratio
  431.      190 Out$=" ": FOR I=1 TO 6:Out$=Out$+"  "+STRING$(10,45): NEXT I: OutPut(Out$,IDevice)
  432.      200 Out$=CR$+" ":PRINT TO Out$ USING "#######.### ";Tim;CumInj;Prs;OilP;WtrP;FacWO: OutPut(Out$,IDevice)
  433.      210 '.....now get time step and injection rate data
  434.      220 ITim=ITim+1: IF ITim>10 OR InjW(ITim)<=0 THEN ERASE Rte0j(),Rte1j(),Rte2j(),Rte3j(): EXIT
  435.      230 '.....this section adds the time step and increments the calculation.....
  436.      240 '.....a fourth-order Runge-Kutta method is used
  437.      250 '.....     the injection in each tube at the start of the step is known
  438.      260 '.....     it is incremented and the injectivity in each tube is recalculated
  439.      270 '.....     using the injectivity, the injection is reallocated between streamtubes
  440.      280 Tim=Tim+DelTim(ITim): ICount=ICount+1: Rte1=0
  441.      290 FOR L=1 TO NMax: WtrIEnd=WtrI(L)+Rte0j(L)/2*DelTim(ITim)/Vol(L)/Rte0*InjW(ITim)*FvfW: Displace(WtrIEnd,OilPQj,WtrPQj,Rte,FrcW,FrcO): Rte1j(L)=Rte*RteCon(L): Rte1=Rte1+Rte1j(L): NEXT L
  442.      300 OilP=0: WtrP=0: FacWO=0: RteW=0: RteO=0: RteNew=0
  443.      310 Rte2=0: FOR L=1 TO NMax: WtrIEnd=WtrI(L)+Rte1j(L)/2*DelTim(ITim)/Vol(L)/Rte1*InjW(ITim)*FvfW: Displace(WtrIEnd,OilPQj,WtrPQj,Rte,FrcW,FrcO): Rte2j(L)=Rte*RteCon(L): Rte2=Rte2+Rte2j(L): NEXT L
  444.      320 Rte3=0: FOR L=1 TO NMax: WtrIEnd=WtrI(L)+Rte2j(L)*DelTim(ITim)/Vol(L)/Rte2*InjW(ITim)*FvfW: Displace(WtrIEnd,OilPQj,WtrPQj,Rte,FrcW,FrcO): Rte3j(L)=Rte*RteCon(L): Rte3=Rte3+Rte3j(L): NEXT L
  445.      330 FOR L=1 TO NMax: RteAv=(Rte0j(L)/Rte0+2*Rte1j(L)/Rte1+2*Rte2j(L)/Rte2+Rte3j(L)/Rte3)/6*InjW(ITim): WtrPQi=RteAv*DelTim(ITim)/Vol(L): WtrI(L)=WtrI(L)+WtrPQi*FvfW
  446.      340   Displace(WtrI(L),OilPQj,WtrPQj,Rte,FrcW,FrcO): Rte0j(L)=Rte*RteCon(L): RteNew=RteNew+Rte0j(L): RteW=RteW+Rte0j(L)*FrcW: RteO=RteO+Rte0j(L)*FrcO: OilP=OilP+OilPQj*Vol(L)/FvfO: WtrP=WtrP+WtrPQj*Vol(L)/FvfW
  447.      350 NEXT L
  448.      360 Rte0=RteNew: Prs=0: IF RteNew>0 THEN Prs=1/(PrsQ(NBuc)*RteNew)
  449.      370 '.....Runge-Kutta solution completed.  Now set up output.
  450.      380 FacWO=0: IF RteW>0 AND RteO>0 THEN FacWO=RteW/FvfW/RteO*FvfO: IF FacWO>999 THEN FacWO=999.999
  451.      390 IF WtrP<0 THEN WtrP=0
  452.      400 IF OilP<0 THEN OilP=0
  453.      410 CumInj=CumInj+InjW(ITim)*DelTim(ITim)
  454.      420 Out$=CR$+" ":PRINT TO Out$ USING "#######.### ";Tim;CumInj;Prs;OilP;WtrP;FacWO: OutPut(Out$,IDevice)
  455.      430 IF ICount MOD 6 = 0 THEN Out$="": OutPut(Out$,IDevice)
  456.      440 IF Tim>=(TimMx(ITim)+1.E-5-DelTim(ITim)/2) GOTO 220 ELSE GOTO 280
  457. END PROCEDURE
  458.  
  459. REAL FUNCTION: Phi
  460.    REAL#: PH,P,B1,B2,B3,B4,B5,Y
  461.    REAL#: T,F,ZZ
  462.       10 '---------------------FUNCTION PHI---------------------------------------
  463.       20 'Function to compute the normal probability integral using the approximation in Abramowitz and Stegun
  464.       30 P=0.2316419: B1=0.319381530: B2=-0.356563782: B3=1.781477937: B4=-1.821255978: B5=1.330274429
  465.       40 Y=Z: IF Z<0 THEN Y=-Z
  466.       50 IF Y>20 THEN Y=20: 'AVOID GETTING OUT OF RANGE
  467.       60 T=1/(1+P*Y)
  468.       70 F=T*((((B5*T+B4)*T+B3)*T+B2)*T+B1): ZZ=EXP(-0.5*Y*Y)/2.506628272
  469.       80 PH=ZZ*F: IF Z>0 THEN PH=1-PH
  470.       90 RESULT=PH
  471. END FUNCTION
  472.  
  473. REAL FUNCTION: PhiInv
  474.    EXTERNAL: Phi
  475.    REAL#: C0,C1,C2,D1,D2,D3,PF,Y
  476.    REAL#: F,T,YP,YP1
  477.       10 '---------------------FUNCTION PHIINV------------------------------------
  478.       20 'Function to compute the inverse normal probability integral 
  479.       30 '   The approximation presented in Abramowitz and Stegun is used for an initial estimate,
  480.       40 '   which is improved by applying 3 iterations of Newton's method to the normal probability integral.
  481.       50 C0=2.515517: C1=0.802853: C2=0.010328
  482.       60 D1=1.432788: D2=0.189269: D3=0.001308: PF=0
  483.       70 IF Z>=1 THEN PF=8.5: GOTO 150
  484.       80 Y=Z: IF Z>0.5 THEN Y=1-Z
  485.       90 IF Y<=0 THEN PF=20.: GOTO 140
  486.      100 F=-2*LOG(Y): T=SQR(F): YP=T-(C0+C1*T+C2*T*T)/(1+T*(D1+T*(D2+T*D3)))
  487.      110 YP1=-YP:YP=YP+2.506628*EXP(0.5*YP*YP)*(Phi(YP1)-Y)
  488.      120 YP1=-YP:YP=YP+2.506628*EXP(0.5*YP*YP)*(Phi(YP1)-Y)
  489.      130 YP1=-YP:PF=YP+2.506628*EXP(0.5*YP*YP)*(Phi(YP1)-Y)
  490.      140 IF Z<0.5 THEN PF=-PF
  491.      150 RESULT = PF
  492. END FUNCTION
  493.  
  494. PROCEDURE: DataOutput
  495.    EXTERNAL: NBuc,NLay,NStr,ITyp,VolP,V,SatOi,SatOr
  496.    EXTERNAL: SatWi,SatIW,Viso,Visw,FvfO,FvfW,PrmrOi,PrmrWr
  497.    EXTERNAL: PrmXpO,PrmXpW,IDevice,OutPut
  498.    STRING: CR$[2],OUT$[?]
  499.       10 '------------------------PROCEDURE DATAOUTPUT---------------------------
  500.       20 '.....this routine outputs the basic data used in the analysis
  501.       30 CR$=CHR$(13)+CHR$(10): OUT$=CR$+CR$+STRING$(75,95): OutPut(OUT$,IDevice)
  502.       40 OUT$=CR$+"  ..........................  INPUT DATA  ............................."+CR$: OutPut (OUT$,IDevice)
  503.       50 OUT$=" ": PRINT TO OUT$ USING "######.###                   Pore Volume";VolP: OutPut(OUT$,IDevice)
  504.       60 OUT$=CR$+" ": PRINT TO OUT$ USING "######                       Number of Saturation Increments";NBuc: OutPut(OUT$,IDevice)
  505.       70 OUT$=" ": PRINT TO OUT$ USING "######                       Number of Layers (0=pseudo rel perms)";NLay: OutPut(OUT$,IDevice)
  506.       80 OUT$=" ": PRINT TO OUT$ USING "######                       Number of Streamlines";NStr: OutPut(OUT$,IDevice)
  507.       90 OUT$=" ": PRINT TO OUT$ USING "######.###                   Coefficient of Permeability Variation";V: OutPut(OUT$,IDevice)
  508.      100 OUT$=CR$+" ": IF NStr<2 GOTO 120 ELSE IF ITyp>0 THEN PRINT TO OUT$ USING "###### spot                  Pattern Type (5, 7, or 9 spot or other)";ITyp: OutPut(OUT$,IDevice): GOTO 130
  509.      110 IF ITyp<0 THEN PRINT TO OUT$ USING "###### spot isolated         Pattern Type (5, 7, or 9 spot or other)";-ITyp: OutPut(OUT$,IDevice): GOTO 130
  510.      120 PRINT TO OUT$  "     ---                     No Pattern Effects Included": OutPut(OUT$,IDevice)
  511.      130 OUT$="": OutPut(OUT$,IDevice)
  512.      140 OUT$=" ": PRINT TO OUT$ USING "######.###   ######.###      ";SatOi;SatWi;:OUT$=OUT$+"Initial Oil & Water Saturations": OutPut(OUT$,IDevice)
  513.      150 OUT$=CR$+" ": PRINT TO OUT$ USING "######.###   ######.###      ";SatOr;SatIW;:OUT$=OUT$+"Residual Oil & Irreducible Water Saturations": OutPut(OUT$,IDevice)
  514.      160 OUT$=CR$+" ": PRINT TO OUT$ USING "######.###   ######.###      ";Viso;Visw;:OUT$=OUT$+"Oil & Water Viscosities": OutPut(OUT$,IDevice)
  515.      170 OUT$=CR$+" ": PRINT TO OUT$ USING "######.###   ######.###      ";PrmrOi;PrmrWr;:OUT$=OUT$+"Oil & Water Endpoint Relative Permeabilities": OutPut(OUT$,IDevice)
  516.      180 OUT$=CR$+" ": PRINT TO OUT$ USING "######.###   ######.###      ";PrmXpO;PrmXpW;:OUT$=OUT$+"Oil & Water Relative Permeability Exponents": OutPut(OUT$,IDevice)
  517.      190 OUT$=CR$+" ": PRINT TO OUT$ USING "######.###   ######.###      ";FvfO;FvfW;:OUT$=OUT$+"Oil & Water Formation Volume Factors": OutPut(OUT$,IDevice)
  518.      200 OUT$=CR$+CR$: OutPut(OUT$,IDevice)
  519.      210 EXIT
  520. END PROCEDURE
  521.  
  522. PROCEDURE: OutPut
  523.    STRING: Out1$[?]
  524.       10 '------------------------PROCEDURE OUTPUT-------------------------------
  525.       20 '.....This is the routine that actually does the output
  526.       30 Out1$=Out$:'.....first need to remove extraneous carriage returns and line feeds
  527.       40 IF LEFT$(Out1$,2)=CHR$(13)+CHR$(10) THEN Out1$=RIGHT$(Out1$,LEN(Out1$)-2)
  528.       50 IF RIGHT$(Out1$,1)=CHR$(13) THEN Out1$=LEFT$(Out1$,LEN(Out1$)-1)
  529.       60 PRINT Out1$;: IF Out$="" THEN PRINT: '..... output results to terminal
  530.       70 IF IDevice<>7 THEN PRINT #IDevice, Out$: EXIT:'.....output to file
  531.       80 IF IDevice=7 THEN LPRINT Out$: EXIT:'..............output to printer
  532. END PROCEDURE
  533.  
  534. 'MAIN Program:
  535.  
  536.  
  537.    10 '---------------------------MAIN PROGRAM--------------------------------
  538.    20 'Program FLOOD.BAS to model waterflood performance using streamtube methods.
  539.    30 'This program computes a one-dimensional Buckley-Leverett displacement for
  540.    40 'individual streamtube performance, with streamtubes being derived for
  541.    50 'various flood networks, and permeability heterogeneity being included in
  542.    60 'the determination of the streamtubes.  Pseudo-relative permeabilities
  543.    70 'using Hearn's method can be incorporated if desired.
  544.    80 '
  545.    90 'Program prepared by Dave O. Cox   4-1987
  546.    91 'Program version 1.1 revised 9-15-87
  547.   100 Tim=TIMER
  548.   105 PRINT "Program FLOOD, Version 1.1"
  549.   110 DataInput
  550.   120 Layer
  551.   130 Stream
  552.   140 Buckley
  553.   150 Calculate
  554.   160 DataOutput
  555.   170 PRINT "Execution Complete.";: Out$=CHR$(13)+CHR$(10)+" ": PRINT TO Out$ USING " Used #####.# seconds.";TIMER-Tim: OutPut(Out$,IDevice)
  556.   180 CLOSE #1: PRINT
  557.   190 PRINT: INPUT "DO YOU WANT TO MAKE ANOTHER RUN [Y/N]" Ans$
  558.   200 Ans$=LEFT$(Ans$,1): IF Ans$="Y" OR Ans$="y" GOTO 100 ELSE END
  559.  
  560. ENDFILE
  561.