TALK=F;RUN( 1, 6: 1,2,3,4,5)

 TEXT(  VSM-2: Calibration mode

  DISPLAY

  This Q1 file is used for searching for the tar retention 
  coefficients of a specific filter material, such as acetate 
  cellulose, paper...

  The investigated system consists of

    (1) one segment of filter
    (2) one segment of tobacco (standard B&H tobacco)
    (3) non-porous cigarette wrap and tipping paper
  
  ENDDIS
  
   GROUP 1. RUN TITLE AND OTHER PRELIMINARIES

   *  DECLARATION OF VARIABLES, REALS AND INTEGERS
REAL(FLTLEN,DIAM,CIRCUM,PI,FLTR1,FLTR2,FLTL0,DIAM0)
REAL(FLTMAS,BLKDEN)
REAL(REFW1,RHO1IN,RHOFLT,REFQ1,QMAX)
REAL(AIRMASS,CFLT,MFLT)
REAL(C1INFM,C2INFM)
REAL(TARCON, RCLOG)
REAL(TOBLEN,TOBRES,TOBR1,TOBR2, TOBP0,TOBL0)
REAL(BUTLEN,BRNLEN,TSPUFF,TSIN,RATIO) 
INTEGER(NPUFF, IVPUFF, NTSPUF, NTSIN, NTFLAT) 
INTEGER(NZFLT0, NZFLT,NZBUT, NZBURN, IZBURN, NZTOB)
INTEGER(ITCAL)

PI=3.14159

   ** Reference dimension
FLTL0=0.01; DIAM0=0.02475/PI

   ** Read data file for smoking condition
INTRPT(R,DATA/OPTCON)
MESG(Smoking condition read from DATA/OPTCON

   * Z-direction veloctiy (maximum) at inlet 
REFW1=REFQ1/(3.14*DIAM0*DIAM0/4)

   *  Read data file for FLTID1 
INTRPT(R,DATA/FLTID1)
MESG(Filter 1 property read from DATA/FLTID1

   *  flter porosity 
BLKDEN=FLTMASS*1.E-6/(PI*DIAM**2./4*FLTLEN)
FLTR1=1.0-BLKDEN/RHOFLT
FLTR2=1.-FLTR1

     * Read data file for tobacco 
INTRPT(R,DATA/TOBID0)
MESG(Tobacco property read from DATA/TOBID0

     * Read data file for numerical setting
INTRPT(R,DATA/NUMSET) 
MESG(Numerical setting read from DATA/NUMSET

RATIO=FLTLEN/FLTL0*10.
NZFLT=NZFLT0*RATIO/10

    GROUP 2. TRANSIENCE; TIME-STEP SPECIFICATION
STEADY=F
LSTEP=NTSPUF*IVPUFF

   The first puff
TFRAC(1)=-NTSIN; TFRAC(2)=0.8/NTSIN
TFRAC(3)=NTFLAT; TFRAC(4)=0.4/NTFLAT
TFRAC(5)=NTSIN;  TFRAC(6)=0.8/NTSIN

INTEGER(KK)
IF(IVPUFF.GT.1) THEN
+  DO JJ=2,IVPUFF
+     KK=JJ-1
+     TFRAC(6*KK+1)=NTSIN;  TFRAC(6*KK+2)=0.8/NTSIN
+     TFRAC(6*KK+3)=NTFLAT; TFRAC(6*KK+4)=0.4/NTFLAT
+     TFRAC(6*KK+5)=NTSIN;  TFRAC(6*KK+6)=0.8/NTSIN
+  ENDDO
ENDIF

   GROUP 3. X-DIRECTION GRID SPECIFICATION
CARTES=F;XULAST=0.01

    GROUP 4. Y-DIRECTION GRID SPECIFICATION
YVLAST=DIAM/2.0

    GROUP 5. Z-DIRECTION GRID SPECIFICATION
NREGZ=3
IREGZ=1; GRDPWR(Z,NZFLT,FLTLEN,1.0)
IREGZ=2; GRDPWR(Z,NZBUT,BUTLEN,1.0)
BRNLEN=TOBLEN-BUTLEN; NZBURN=IZBURN*NPUFF
IREGZ=3; GRDPWR(Z,NZBURN,BRNLEN,1.0)
NZTOB=NZBURN+NZBUT

    GROUP 7. VARIABLES STORED, SOLVED & NAMED
SOLVE(P1,W1,W2,R1,R2,TAR1,TAR2)
    ** FFV VLM used for calculating tar retention in ground
STORE(VLM,ACEL,DEN1,DEN2)
ONEPHS=F

   USE WHOLE-FIELD SOLVER FOR ALL VARIABLES
SOLUTN(P1,Y,Y,Y,P,P,P);SOLUTN(W1,Y,Y,Y,P,P,P)
SOLUTN(W2,Y,Y,Y,P,P,P);SOLUTN(TAR1,Y,Y,Y,P,P,P)
SOLUTN(TAR2,Y,Y,Y,P,P,P)
 
    GROUP 8. TERMS (IN DIFFERENTIAL EQUATIONS) & DEVICES
TERMS(P1,Y,Y,Y,Y,Y,Y);TERMS(TAR1,N,Y,P,P,Y,N)
TERMS(TAR2,N,N,P,P,N,N)

    GROUP 9. PROPERTIES OF THE MEDIUM (OR MEDIA)
AIRMASS=286; PRESS0=1E+5; TEMP0=273.0
RHO1=PRESS0/AIRMASS/(TEMP0+20)
RHO2=RHOFLT;ENUL=1.E-5;RHO1IN=RHO1
PRNDTL(TAR1)=-1.0E-6;  PRNDTL(TAR2)=-1.0E-12

     ** STORE C1 FOR USE IN IN-FORM
C1INFM=9.81*0.0248**5./(4*PI*MFLT*REFQ1)
C2INFM=(FLTMAS-CFLT*FLTLEN*1.E+3/90.)*FLTR1**2./FLTLEN/CIRCUM**3.

    GROUP 10. INTER-PHASE-TRANSFER PROCESSES AND PROPERTIES
    ** SET A CONSTANT INTER-PHASE FRICTION COEFFICIENT.
CFIPS=1.E-5

    ** HEAT AND TAR TRANSFER BETWEEN PHASES IS HANDLED BY >PATCHES
CINT(TAR1)=0.0; CINT(TAR2)=0.0
 
    GROUP 11. INITIALIZATION OF VARIABLE OR POROSITY FIELDS
INIADD=F
FIINIT(P1)=0.0;FIINIT(W1)=0.0; FIINIT(W2)=0.0
FIINIT(R1)=TOBR1; FIINIT(R2)=TOBR2
FIINIT(TAR1)=0.0; FIINIT(TAR2)=0.0

PATCH(FLTR,INIVAL,1,NX,1,NY,1,NZFLT,1,1)
INIT(FLTR,R1,0.0,FLTR1);INIT(FLTR,R2,0.0,FLTR2)

    GROUP 13. BOUNDARY CONDITIONS AND SPECIAL SOURCES

     ** AIR OUTLET AT MOUTH
REAL(WMAX); WMAX=QMAX/(PI*DIAM**2./4.)
REAL(QOFFST); QOFFST=0.5/NTSIN/4.

     ** In the following do-loop, each cycle represents a puff,
        smoulding period is neglected.

     ** NZTIP is the IZ value for the burning tip at the beginning of each puff
        TMFRST, TMLAST are the ISTEPs at the beginning and end of each puff
        there are 25 time steps in each puff

     ** TARTIP is the mass fraction of TAR in gas phase at the burning tip
INTEGER(NZTIP,TMFRST,TMLAST)
REAL(TARTIP); TARTIP=0.05
     
DO II=1, IVPUFF
+ NZTIP=NZ-(II-1)*IZBURN
+ TMFRST=(II-1)*NTSPUF+1
+ TMLAST=II*NTSPUF

     ** AIR INLETT 
+ PATCH(INTOB:II:,HIGH,1,NX,1,NY,NZTIP,NZTIP,TMFRST,TMLAST)
+ COVAL(INTOB:II:,P1,1.0,0.0)
+ COVAL(INTOB:II:,TAR1,FIXVAL,TARTIP)

     * ONE-FORTH OF SIN-WAVE
+ TIMA=1
+ PATCH(TIMOU:II:1,LOW,1,NX,1,NY,1,1,TMFRST,TMFRST-1+NTSIN)
+ COVAL(TIMOU:II:1,P1,FIXFLU,GRND10)
+ SPEDAT(SET,TIMOU:II:1,PERIODP1,R,TSIN*4.)
+ SPEDAT(SET,TIMOU:II:1,AMPLITUDEP1,R,-RHO1IN*WMAX)
+ SPEDAT(SET,TIMOU:II:1,MEANP1,R,0.0)
+ SPEDAT(SET,TIMOU:II:1,OFFSETP1,R,QOFFST)
+ SPEDAT(SET,TIMOU:II:1,TSTARTP1,R,(II-1)*2.0)

    * CONSTANT FLOWRATE
+ PATCH(TIMOU:II:2,LOW,1,NX,1,NY,1,1,TMFRST+NTSIN,$
                                    TMFRST+NTSIN+NTFLAT-1)
+ COVAL(TIMOU:II:2,P1,FIXFLU,(-1.)*RHO1IN*WMAX)

    * ONE FORTH OF COS-WAVE
+ PATCH(TIMOU:II:3,LOW,1,NX,1,NY,1,1,TMLAST-NTSIN+1,TMLAST)
+ COVAL(TIMOU:II:3,P1,FIXFLU,GRND6)
+ SPEDAT(SET,TIMOU:II:3,PERIODP1,R,TSIN*4.)
+ SPEDAT(SET,TIMOU:II:3,AMPLITUDEP1,R,-RHO1IN*WMAX)
+ SPEDAT(SET,TIMOU:II:3,MEANP1,R,0.0)
+ SPEDAT(SET,TIMOU:II:3,OFFSETP1,R,QOFFST)
+ SPEDAT(SET,TIMOU:II:3,TSTARTP1,R,(II-1)*2.0+1.2)

ENDDO

      ** SOLID is stationary
PATCH(SOLID,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP)
COVAL(SOLID,W2,FIXVAL,0.0)
 
     ** adding FLOW RESISTANCE OF TOBACCO
PATCH(TOBACO,VOLUME,1,NX,1,NY,NZFLT+1,NZ,1,LSTEP)
COVAL(TOBACO,W1,TOBRES*TOBR1,0.0)

    ** FLOW RESISTANCE OF FILTER
PATCH(FILT,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)
COVAL(FILT,W1,C1INFM*C2INFM,0.)
 
  ** SOURCE OF SOLID-PHASE AND GAS-PHASE TAR
PATCH(SCTAR1,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)
 (SOURCE OF TAR1 AT SCTAR1 is (999001)*R2^(999003)*$
 ((999002)*TAR2-TAR1) WITH LINE)

PATCH(SCTAR2,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)
 (SOURCE OF TAR2 AT SCTAR2 is (999001)*R2^(999003)*$
 (TAR1-(999002)*TAR2) WITH LINE)

IG(1)=NZTOB;  IG(2)=IZBURN 
IG(3)=IVPUFF; IG(4)=NPUFF;  IG(5)=NTSPUF; IG(6)=ITCAL 
RG(1)=TOBLEN; RG(2)=BUTLEN; RG(3)=TARTIP

    GROUP 15. TERMINATION OF SWEEPS
LSWEEP=75; RESFAC=0.01

    GROUP 17. UNDER-RELAXATION DEVICES
REAL(DTF);DTF=ZWLAST/NZ/REFW1*0.001
RELAX(P1,LINRLX,0.3);RELAX(W1,FALSDT,DTF)
RELAX(R1,LINRLX,0.0);RELAX(R2,LINRLX,0.0)
RELAX(TAR1,LINRLX,0.1);RELAX(TAR2,LINRLX,0.1)

   GROUP 22. SPOT-VALUE PRINT-OUT
NPRMON=20;IZMON=2;NZPRIN=1;NTPRIN=1

    GROUP 23. FIELD PRINT-OUT AND PLOT CONTROL
OUTPUT(W2,N,N,N,N,N,N);OUTPUT(C1,N,N,N,N,N,N)
OUTPUT(ACEL,N,N,N,N,N,N);OUTPUT(R2,N,N,Y,N,N,N)
NOWIPE=T;TSTSWP=-1
NSAVE=P001
IDISPA=4;IDISPB=4;IDISPC=LSTEP
CSG1=A

STOP

INTRPT(R,DATA/FLTID2)
MESG(Filter 2 property read from DATA/FLTID2
#PAUSE

RATIO=FLTLEN/FLTL0*10.
NZFLT=NZFLT0*RATIO/10

IREGZ=1; GRDPWR(Z,NZFLT,FLTLEN,1.0)

DIAM=CIRCUM/PI; YVLAST=DIAM/2.

   *  flter porosity 
BLKDEN=FLTMASS*1.E-6/(PI*DIAM**2./4*FLTLEN)
FLTR1=1.0-BLKDEN/RHOFLT
FLTR2=1.-FLTR1

C2INFM=(FLTMAS-CFLT*FLTLEN*1.E+3/90.)*FLTR1**2./FLTLEN/CIRCUM**3.

PATCH(FLTR,INIVAL,1,NX,1,NY,1,NZFLT,1,1)
INIT(FLTR,R1,0.0,FLTR1)
INIT(FLTR,R2,0.0,FLTR2)

WMAX=QMAX/(PI*DIAM**2./4.)

DO II=1, IVPUFF
+ NZTIP=NZ-(II-1)*IZBURN
+ TMFRST=(II-1)*NTSPUF+1
+ TMLAST=II*NTSPUF

     ** AIR INLETT 
+ PATCH(INTOB:II:,HIGH,1,NX,1,NY,NZTIP,NZTIP,TMFRST,TMLAST)

     * ONE-FORTH OF SIN-WAVE
+ TIMA=1
+ PATCH(TIMOU:II:1,LOW,1,NX,1,NY,1,1,TMFRST,TMFRST-1+NTSIN)
+ SPEDAT(SET,TIMOU:II:1,AMPLITUDEP1,R,-RHO1IN*WMAX)

    * CONSTANT FLOWRATE
+ PATCH(TIMOU:II:2,LOW,1,NX,1,NY,1,1,TMFRST+NTSIN,$
                                    TMFRST+NTSIN+NTFLAT-1)
+ COVAL(TIMOU:II:2,P1,FIXFLU,(-1.)*RHO1IN*WMAX)

    * ONE FORTH OF COS-WAVE
+ PATCH(TIMOU:II:3,LOW,1,NX,1,NY,1,1,TMLAST-NTSIN+1,TMLAST)
+ SPEDAT(SET,TIMOU:II:3,AMPLITUDEP1,R,-RHO1IN*WMAX)

ENDDO

      ** SOLID is stationary
PATCH(SOLID,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP)
 
     ** 050601: adding FLOW RESISTANCE OF TOBACCO
PATCH(TOBACO,VOLUME,1,NX,1,NY,NZFLT+1,NZ,1,LSTEP)

    ** FLOW RESISTANCE OF FILTER
PATCH(FILT,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)
COVAL(FILT,W1,C1INFM*C2INFM,0.)
 
PATCH(SCTAR1,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)
PATCH(SCTAR2,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)

IG(1)=NZTOB;  IG(2)=IZBURN 
IG(3)=IVPUFF; IG(4)=NPUFF;  IG(5)=NTSPUF; IG(6)=ITCAL 

RG(1)=TOBLEN; RG(2)=BUTLEN; RG(3)=TARTIP
NSAVE=P002

STOP

INTRPT(R,DATA/FLTID3)
MESG(Filter 3 property read from DATA/FLTID3
#PAUSE

RATIO=FLTLEN/FLTL0*10.
NZFLT=NZFLT0*RATIO/10

IREGZ=1; GRDPWR(Z,NZFLT,FLTLEN,1.0)

DIAM=CIRCUM/PI; YVLAST=DIAM/2.

   *  flter porosity 
BLKDEN=FLTMASS*1.E-6/(PI*DIAM**2./4*FLTLEN)
FLTR1=1.0-BLKDEN/RHOFLT
FLTR2=1.-FLTR1

C2INFM=(FLTMAS-CFLT*FLTLEN*1.E+3/90.)*FLTR1**2./FLTLEN/CIRCUM**3.

PATCH(FLTR,INIVAL,1,NX,1,NY,1,NZFLT,1,1)
INIT(FLTR,R1,0.0,FLTR1)
INIT(FLTR,R2,0.0,FLTR2)

WMAX=QMAX/(PI*DIAM**2./4.)

DO II=1, IVPUFF
+ NZTIP=NZ-(II-1)*IZBURN
+ TMFRST=(II-1)*NTSPUF+1
+ TMLAST=II*NTSPUF

     ** AIR INLETT 
+ PATCH(INTOB:II:,HIGH,1,NX,1,NY,NZTIP,NZTIP,TMFRST,TMLAST)

     * ONE-FORTH OF SIN-WAVE
+ TIMA=1
+ PATCH(TIMOU:II:1,LOW,1,NX,1,NY,1,1,TMFRST,TMFRST-1+NTSIN)
+ SPEDAT(SET,TIMOU:II:1,AMPLITUDEP1,R,-RHO1IN*WMAX)

    * CONSTANT FLOWRATE
+ PATCH(TIMOU:II:2,LOW,1,NX,1,NY,1,1,TMFRST+NTSIN,$
                                    TMFRST+NTSIN+NTFLAT-1)
+ COVAL(TIMOU:II:2,P1,FIXFLU,(-1.)*RHO1IN*WMAX)

    * ONE FORTH OF COS-WAVE
+ PATCH(TIMOU:II:3,LOW,1,NX,1,NY,1,1,TMLAST-NTSIN+1,TMLAST)
+ SPEDAT(SET,TIMOU:II:3,AMPLITUDEP1,R,-RHO1IN*WMAX)

ENDDO

      ** SOLID is stationary
PATCH(SOLID,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP)
 
     ** 050601: adding FLOW RESISTANCE OF TOBACCO
PATCH(TOBACO,VOLUME,1,NX,1,NY,NZFLT+1,NZ,1,LSTEP)

    ** FLOW RESISTANCE OF FILTER
PATCH(FILT,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)
COVAL(FILT,W1,C1INFM*C2INFM,0.)

 
  ** SOURCE OF SOLID-PHASE AND GAS-PHASE TAR
  PATCH(>TAR1,VOLUME,1,1,1,NY,1,NZ,1,LSTEP)
  COVAL(>TAR1,TAR2, TARCON*RCLOG, 1.0/RCLOG)

  PATCH(>TAR2,VOLUME,1,1,1,NY,1,NZ,1,LSTEP)
  COVAL(>TAR2,TAR1, TARCON, RCLOG)

PATCH(SCTAR1,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)
PATCH(SCTAR2,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)

IG(1)=NZTOB;  IG(2)=IZBURN 
IG(3)=IVPUFF; IG(4)=NPUFF;  IG(5)=NTSPUF; IG(6)=ITCAL 

RG(1)=TOBLEN; RG(2)=BUTLEN; RG(3)=TARTIP


NSAVE=P003

STOP


INTRPT(R,DATA/FLTID4)
MESG(Filter 4 property read from DATA/FLTID4
#PAUSE

RATIO=FLTLEN/FLTL0*10.
NZFLT=NZFLT0*RATIO/10

IREGZ=1; GRDPWR(Z,NZFLT,FLTLEN,1.0)

DIAM=CIRCUM/PI; YVLAST=DIAM/2.

   *  flter porosity 
BLKDEN=FLTMASS*1.E-6/(PI*DIAM**2./4*FLTLEN)
FLTR1=1.0-BLKDEN/RHOFLT
FLTR2=1.-FLTR1

C2INFM=(FLTMAS-CFLT*FLTLEN*1.E+3/90.)*FLTR1**2./FLTLEN/CIRCUM**3.

PATCH(FLTR,INIVAL,1,NX,1,NY,1,NZFLT,1,1)
INIT(FLTR,R1,0.0,FLTR1)
INIT(FLTR,R2,0.0,FLTR2)

WMAX=QMAX/(PI*DIAM**2./4.)

DO II=1, IVPUFF
+ NZTIP=NZ-(II-1)*IZBURN
+ TMFRST=(II-1)*NTSPUF+1
+ TMLAST=II*NTSPUF

     ** AIR INLETT 
+ PATCH(INTOB:II:,HIGH,1,NX,1,NY,NZTIP,NZTIP,TMFRST,TMLAST)

     * ONE-FORTH OF SIN-WAVE
+ TIMA=1
+ PATCH(TIMOU:II:1,LOW,1,NX,1,NY,1,1,TMFRST,TMFRST-1+NTSIN)
+ SPEDAT(SET,TIMOU:II:1,AMPLITUDEP1,R,-RHO1IN*WMAX)

    * CONSTANT FLOWRATE
+ PATCH(TIMOU:II:2,LOW,1,NX,1,NY,1,1,TMFRST+NTSIN,$
                                    TMFRST+NTSIN+NTFLAT-1)
+ COVAL(TIMOU:II:2,P1,FIXFLU,(-1.)*RHO1IN*WMAX)

    * ONE FORTH OF COS-WAVE
+ PATCH(TIMOU:II:3,LOW,1,NX,1,NY,1,1,TMLAST-NTSIN+1,TMLAST)
+ SPEDAT(SET,TIMOU:II:3,AMPLITUDEP1,R,-RHO1IN*WMAX)

ENDDO

      ** SOLID is stationary
PATCH(SOLID,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP)
 
     ** 050601: adding FLOW RESISTANCE OF TOBACCO
PATCH(TOBACO,VOLUME,1,NX,1,NY,NZFLT+1,NZ,1,LSTEP)

    ** FLOW RESISTANCE OF FILTER
PATCH(FILT,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)
COVAL(FILT,W1,C1INFM*C2INFM,0.)

 
  ** SOURCE OF SOLID-PHASE AND GAS-PHASE TAR
  PATCH(>TAR1,VOLUME,1,1,1,NY,1,NZ,1,LSTEP)
  COVAL(>TAR1,TAR2, TARCON*RCLOG, 1.0/RCLOG)

  PATCH(>TAR2,VOLUME,1,1,1,NY,1,NZ,1,LSTEP)
  COVAL(>TAR2,TAR1, TARCON, RCLOG)

PATCH(SCTAR1,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)
PATCH(SCTAR2,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)

IG(1)=NZTOB;  IG(2)=IZBURN 
IG(3)=IVPUFF; IG(4)=NPUFF;  IG(5)=NTSPUF; IG(6)=ITCAL 

RG(1)=TOBLEN; RG(2)=BUTLEN; RG(3)=TARTIP


NSAVE=P004

STOP


INTRPT(R,DATA/FLTID5)
MESG(Filter 5 property read from DATA/FLTID5
#PAUSE

RATIO=FLTLEN/FLTL0*10.
NZFLT=NZFLT0*RATIO/10

IREGZ=1; GRDPWR(Z,NZFLT,FLTLEN,1.0)

DIAM=CIRCUM/PI; YVLAST=DIAM/2.

   *  flter porosity 
BLKDEN=FLTMASS*1.E-6/(PI*DIAM**2./4*FLTLEN)
FLTR1=1.0-BLKDEN/RHOFLT
FLTR2=1.-FLTR1

C2INFM=(FLTMAS-CFLT*FLTLEN*1.E+3/90.)*FLTR1**2./FLTLEN/CIRCUM**3.

PATCH(FLTR,INIVAL,1,NX,1,NY,1,NZFLT,1,1)
INIT(FLTR,R1,0.0,FLTR1)
INIT(FLTR,R2,0.0,FLTR2)

WMAX=QMAX/(PI*DIAM**2./4.)

DO II=1, IVPUFF
+ NZTIP=NZ-(II-1)*IZBURN
+ TMFRST=(II-1)*NTSPUF+1
+ TMLAST=II*NTSPUF

     ** AIR INLETT 
+ PATCH(INTOB:II:,HIGH,1,NX,1,NY,NZTIP,NZTIP,TMFRST,TMLAST)

     * ONE-FORTH OF SIN-WAVE
+ TIMA=1
+ PATCH(TIMOU:II:1,LOW,1,NX,1,NY,1,1,TMFRST,TMFRST-1+NTSIN)
+ SPEDAT(SET,TIMOU:II:1,AMPLITUDEP1,R,-RHO1IN*WMAX)

    * CONSTANT FLOWRATE
+ PATCH(TIMOU:II:2,LOW,1,NX,1,NY,1,1,TMFRST+NTSIN,$
                                    TMFRST+NTSIN+NTFLAT-1)
+ COVAL(TIMOU:II:2,P1,FIXFLU,(-1.)*RHO1IN*WMAX)

    * ONE FORTH OF COS-WAVE
+ PATCH(TIMOU:II:3,LOW,1,NX,1,NY,1,1,TMLAST-NTSIN+1,TMLAST)
+ SPEDAT(SET,TIMOU:II:3,AMPLITUDEP1,R,-RHO1IN*WMAX)

ENDDO

      ** SOLID is stationary
PATCH(SOLID,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP)
 
     ** 050601: adding FLOW RESISTANCE OF TOBACCO
PATCH(TOBACO,VOLUME,1,NX,1,NY,NZFLT+1,NZ,1,LSTEP)

    ** FLOW RESISTANCE OF FILTER
PATCH(FILT,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)
COVAL(FILT,W1,C1INFM*C2INFM,0.)

 
  ** SOURCE OF SOLID-PHASE AND GAS-PHASE TAR
  PATCH(>TAR1,VOLUME,1,1,1,NY,1,NZ,1,LSTEP)
  COVAL(>TAR1,TAR2, TARCON*RCLOG, 1.0/RCLOG)

  PATCH(>TAR2,VOLUME,1,1,1,NY,1,NZ,1,LSTEP)
  COVAL(>TAR2,TAR1, TARCON, RCLOG)

PATCH(SCTAR1,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)
PATCH(SCTAR2,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)

IG(1)=NZTOB;  IG(2)=IZBURN 
IG(3)=IVPUFF; IG(4)=NPUFF;  IG(5)=NTSPUF; IG(6)=ITCAL 

RG(1)=TOBLEN; RG(2)=BUTLEN; RG(3)=TARTIP


NSAVE=P005

STOP



INTRPT(R,DATA/FLTID6)
MESG(Filter 6 property read from DATA/FLTID6
#PAUSE

RATIO=FLTLEN/FLTL0*10.
NZFLT=NZFLT0*RATIO/10

IREGZ=1; GRDPWR(Z,NZFLT,FLTLEN,1.0)

DIAM=CIRCUM/PI; YVLAST=DIAM/2.

   *  flter porosity 
BLKDEN=FLTMASS*1.E-6/(PI*DIAM**2./4*FLTLEN)
FLTR1=1.0-BLKDEN/RHOFLT
FLTR2=1.-FLTR1

C2INFM=(FLTMAS-CFLT*FLTLEN*1.E+3/90.)*FLTR1**2./FLTLEN/CIRCUM**3.

PATCH(FLTR,INIVAL,1,NX,1,NY,1,NZFLT,1,1)
INIT(FLTR,R1,0.0,FLTR1)
INIT(FLTR,R2,0.0,FLTR2)

WMAX=QMAX/(PI*DIAM**2./4.)

DO II=1, IVPUFF
+ NZTIP=NZ-(II-1)*IZBURN
+ TMFRST=(II-1)*NTSPUF+1
+ TMLAST=II*NTSPUF

     ** AIR INLETT 
+ PATCH(INTOB:II:,HIGH,1,NX,1,NY,NZTIP,NZTIP,TMFRST,TMLAST)

     * ONE-FORTH OF SIN-WAVE
+ TIMA=1
+ PATCH(TIMOU:II:1,LOW,1,NX,1,NY,1,1,TMFRST,TMFRST-1+NTSIN)
+ SPEDAT(SET,TIMOU:II:1,AMPLITUDEP1,R,-RHO1IN*WMAX)

    * CONSTANT FLOWRATE
+ PATCH(TIMOU:II:2,LOW,1,NX,1,NY,1,1,TMFRST+NTSIN,$
                                    TMFRST+NTSIN+NTFLAT-1)
+ COVAL(TIMOU:II:2,P1,FIXFLU,(-1.)*RHO1IN*WMAX)

    * ONE FORTH OF COS-WAVE
+ PATCH(TIMOU:II:3,LOW,1,NX,1,NY,1,1,TMLAST-NTSIN+1,TMLAST)
+ SPEDAT(SET,TIMOU:II:3,AMPLITUDEP1,R,-RHO1IN*WMAX)

ENDDO

      ** SOLID is stationary
PATCH(SOLID,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP)
 
     ** 050601: adding FLOW RESISTANCE OF TOBACCO
PATCH(TOBACO,VOLUME,1,NX,1,NY,NZFLT+1,NZ,1,LSTEP)

    ** FLOW RESISTANCE OF FILTER
PATCH(FILT,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)
COVAL(FILT,W1,C1INFM*C2INFM,0.)

 
  ** SOURCE OF SOLID-PHASE AND GAS-PHASE TAR
  PATCH(>TAR1,VOLUME,1,1,1,NY,1,NZ,1,LSTEP)
  COVAL(>TAR1,TAR2, TARCON*RCLOG, 1.0/RCLOG)

  PATCH(>TAR2,VOLUME,1,1,1,NY,1,NZ,1,LSTEP)
  COVAL(>TAR2,TAR1, TARCON, RCLOG)

PATCH(SCTAR1,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)
PATCH(SCTAR2,VOLUME,1,NX,1,NY,1,NZFLT,1,LSTEP)

IG(1)=NZTOB;  IG(2)=IZBURN 
IG(3)=IVPUFF; IG(4)=NPUFF;  IG(5)=NTSPUF; IG(6)=ITCAL 

RG(1)=TOBLEN; RG(2)=BUTLEN; RG(3)=TARTIP


NSAVE=P006

STOP