TALK=T;RUN( 1, 1)
 ** LOAD(x213) from the x Input Library
     ** PHOENICS VALIDATION CASE
TEXT(CHEN-KIM_2D ABRUPT PIPE EXPANS  :T213
TITLE
  DISPLAY
  The problem considered is the calculation of incompressible
  turbulent flow and heat transfer in a sudden pipe expansion 
  with a diameter ratio Do/Di=2.5. The Reynolds number of the
  larger pipe is 4.075E4, and the laminar Prandtl number is 0.7. 
  Downstream of the expansion, the fluid is heated by a uniform
  heat flux through the outer wall. This case was studied 
  experimentally by Baugh et al [1984] to determine the local 
  Nusselt number distributions for various Reynolds numbers. 
  
  Following the expansion, the flow enters the larger pipe in the 
  form of a circular jet, and then separates from the expansion
  corner to generate to a primary recirculation zone, whose reverse 
  flow is driven by the adverse pressure gradient associated with
  the expansion geometry. At the point of reattachment, the heat
  transfer coefficients (htc) are known to be several times greater 
  than the corresponding fully-developed values. For this case, the 
  measurements show that the peak htc downstream of the expansion 
  is roughly 4 times the fully-developed htc for the same Reynolds 
  number. These high htcs occur in the region of the reattachment 
  of the shear layer to the tube wall.
 
  The calculations are started at z= -4h and terminated at z=41h,
  where a fixed-pressure boundary condition is applied. Here, h is
  the step height given by h=0.5*(Do-Di).The calculation may be 
  performed with one of the following low-Re turbulence models: the 
  Lam-Bremhorst, Chen-Kim k-e & 2-Layer k-e models; the Wilcox 
  1988 & 2008 k-w models; and the Menter & SST k-w models. The 
  turbulent Prandtl number is taken as 0.86. The numerical 
  integration is taken right down to the wall, and a non-uniform 
  mesh of NY=110 by NX=200 is used which concentrates grid cells 
  close to the wall. Typically, about 5000 sweeps are required 
  for a converged solution with this mesh density. 
  ENDDIS
 
  For simplicity, uniform profiles of w, T, k and e (or w) are specified 
  at the inlet to the calculation domain. The main results of the 
  simulations are compared with the measured values below:
                         
                   LB    CK    2L   k-w  k-w-R  k-w-M  k-w-SST Data
  Xr/H         =   8.4  11.0 10.6  11.0   17.0   10.7  12.5  10.0
  Nu/NuDB,max  =  16.8  15.7  3.1   3.5    2.8    3.7   3.4   4.2
 
  Computations made using the LB and CK low-Re k-e models predict rates 
  of heat transfer that are too high by up to a factor of 4 as a result 
  of too large predicted levels of near-wall length scale. The other 
  models fare much better, but these solutions need to be checked for
  grid dependency. It is evident that the computed levels of the local 
  htc are especially sensitive to the near-wall model.
  
   Baughn, J. W., Hoffman, M. A., Takahashi, R. and Launder, B. E.,
   "Local Heat Transfer Downstream of an Abrupt Expansion in a
    Circular Channel With Constant Wall Heat Flux", ASME Journal of Heat 
	Transfer, Vol. 106, 1984, pp.789-786.
  
   Zemanick, P. P., and Dougall, R. S., "Local Heat Transfer Downstream
   of Abrupt Circular Channel Expansion," ASME Journal of Heat Transfer,
   Vol. 92, 1970, p. 53.
 
    This AUTOPLOT sequence provides a plot of the axial velocity W1 along 
	the outer surface of the solution domain versus normalised axial 
	distance Z. The axial coordinate 0.0 corresponds to the step location. 
	The reattachment point corresponds the axial location where W1 changes 
	from negative to positive.
 
   AUTOPLOT USE
   file
   phida 3
 
 
   da 1 w1 y m
   shift x -1.2 1
   divide x .3 1
   col9 1
   level y 0;level x 0
   scale x 0 40
   msg Velocity (W1) profile
   pause
   cl
   da 1 ypls y m
   shift x -1.2 1
   divide x .3 1
   col3 1
   level y 0;level x 0
   scale x 0 40
   pause
   cl
   da 1 ypls y 50
   divide x .3 1
   col3 1
   level y 0;level x 0
   scale x 0 5
   pause
   cl
   da 1 nuss y m
   shift x -1.2 1
   divide x .3 1
   plot 1
   level y 0;level x 0
   scale x 0 40  
   pause
   cl
   da 1 nusc y m
   shift x -1.2 1
   divide x .3 1
   plot 1
   level y 0;level x 0
   scale x 0 40   
   msg Press e to END
   ENDUSE
 
BOOLEAN(KWMOD);KWMOD=F
CHAR(CTURB,TLSC);INTEGER(TMODEL,NYS,NYS2,NZS,NZS2,NZS3,NZS4)
REAL(HGHT,DIAMI,DIAMO,PL_UP,REYDO,GRADI,RADO,TIN,PLDOWN)
REAL(PLD_S1,PLD_S2,PLD_S3,REYDI,NUZD,NUZDS,CP_BC,PIND)
REAL(PLEN,WIN,WOUT,TKEIN,EPSIN,VFLOW,AIN,QIN,AREAO,RHOIN)
REAL(PRLAM,EXRAT,MFLOW,SPHEAT,QWALL,DTRISE,NUDB,DRAT)
REAL(FRIC,WSTAR,KEMAX,EPMAX,ENUTIN,OMEGIN)
   ** Reynolds number, Prandtl number, Expansion Ratio  
REYDO=4.075E4
PRLAM=0.7;EXRAT=2.5 
REYDO;PRLAM  
RHOIN=1.0 ! uniform density

     ** Calculation of domain specifications
DIAMO=1.0;DIAMI=DIAMO/EXRAT;WIN=10.0
DRAT=DIAMI/DIAMO

WOUT=WIN*(DRAT)**2
ENUL=(WOUT*DIAMO)/REYDO

RADO=0.5*DIAMO;GRADI=0.5*DIAMI;HGHT=RADO-GRADI

DRAT;ENUL;WIN;WOUT;DIAMO;HGHT

  ** PLEN is the total pipe length
PL_UP =4.*HGHT
PLD_S1=HGHT
PLD_S2=15*HGHT
PLD_S3=25*HGHT
PLDOWN=PLD_S1+PLD_S2+PLD_S3
PLEN=PL_UP+PLDOWN

PL_UP;PLEN;PLDOWN

   ** Inlet Reynolds number
REYDI=WIN*DIAMI/ENUL
REYDI
     ** Estimate friction velocity
FRIC=1./(1.82*LOG10(REYDI)-1.64)**2;WSTAR=WIN*(FRIC/8.)**0.5

FRIC;WSTAR

   ** Inlet turbulence values
TKEIN=0.25*WIN*WIN*FRIC;EPSIN=0.1643*TKEIN**1.5/(0.09*GRADI)
ENUTIN=0.09*TKEIN*TKEIN/EPSIN;OMEGIN=EPSIN/(0.09*TKEIN)
TKEIN;EPSIN;ENUTIN;OMEGIN
  ** Inlet temperature
TIN=20.
   ** Nusselt number for fully-developed pipe flow (Dittus Boelter)
NUDB=0.023*(REYDO**0.8)*(PRLAM**0.4)
NUDB 
  ** Zemanick & Dougall (1970) Nu,max correlation
NUZD=0.2*REYDI**(2./3.);NUZDS=NUZD/NUDB
NUZDS
  ** Inlet temperature
TIN=20.
  ** Borda-Carnot Pressure coefficient - maximum pressure rise 
     in the absence of any wall friction
CP_BC=2.*(DRAT**2)*(1.-DRAT**2)
CP_BC  
  ** dimensionless total-pressure loss due to pipe expansion is 
    DP,LOSS= (1.-DRAT**2)**2
  ** Inlet Dynamic pressure
PIND=0.5*RHOIN*WIN*WIN

  ** pressure coefficient
(STORED of CP is (P1-P1[1,1,1])/PIND with IMAT<100) 
  ** temperature difference
(STORED of TDIF is (TEM1-:TIN:) with IMAT<100)
    GROUP 3. X-direction grid specification
CARTES=F;XULAST=0.1
    ** inlet flow area and volume flow rate
AIN=DIAMI*DIAMI*XULAST/8.;VFLOW=WIN*AIN
    GROUP 4. Y-direction grid specification
    GROUP 5. Z-direction grid specification

NYS=50;NYS2=60

NREGY=2
    ** -ve no.of cells means symmetric
IREGY=1;GRDPWR(Y,NYS,-GRADI,-1.08)
IREGY=2;GRDPWR(Y,-NYS2,-HGHT,1.10)

   ** -ve distance means Geometric Progression
NZS=45;NZS2=35;NZS3=80;NZS4=40
NREGZ=4
IREGZ=1;GRDPWR(Z,NZS,-PL_UP,-1.05)
IREGZ=2;GRDPWR(Z,NZS2,-PLD_S1,1.06)
IREGZ=3;GRDPWR(Z,NZS3,-PLD_S2,1.015)
IREGZ=4;GRDPWR(Z,NZS4,-PLD_S3,1.02)

    GROUP 6. Body-fitted coordinates or grid distortion
    GROUP 7. Variables stored, solved & named
SOLVE(P1,W1,V1,TEM1);SOLUTN(P1,Y,Y,Y,N,N,N)
SOLUTN(W1,P,P,P,P,P,N);SOLUTN(V1,P,P,P,P,P,N)
SOLUTN(TEM1,Y,Y,Y,N,N,N);STORE(ENUT)
STORE(TDIF,CP)
MESG( Enter the required turbulence model:
MESG(  LB  - Low-Re Lam-Brem. k-e model
MESG(  CK  - Low-Re Lam-Brem. k-e model 
MESG(  2L  - Low-Re 2-layer   k-e model
MESG(  KW  - Wilcox 1988 Low-Re k-w model
MESG(  KWR - Wilcox 2008 Low-Re k-w model
MESG(  KWM - Menter Low-Re k-w model
MESG(  KWS - Low-Re k-w SST model (default)
MESG(
READVDU(CTURB,CHAR,KWS)
CTURB
CASE :CTURB: OF
WHEN LB,2
+ TEXT(LAM-BRE K-E_2D ABRUPT PIPE EXPANS  :T213
+ MESG(Low-Re Lam-Bem. k-e model
+ TMODEL=1;KELIN=1;TLSC=EP
+ TURMOD(KEMODL-LOWRE);STORE(REYN)
   + SOLUTN(V1,Y,Y,Y,P,P,P);SOLUTN(W1,Y,Y,Y,P,P,P)
WHEN CK,2
+ TEXT(CHEN-KIM_2D ABRUPT PIPE EXPANS  :T213
+ MESG(Low-Re Chen-Kim k-e model
+ TMODEL=2;KELIN=1;TLSC=EP
+ TURMOD(KECHEN-LOWRE);STORE(REYN)
   + SOLUTN(V1,Y,Y,Y,P,P,P);SOLUTN(W1,Y,Y,Y,P,P,P)
WHEN 2L,2
+ TEXT(2-LAYER K-E_2D ABRUPT PIPE EXPANS  :T213
+ MESG(Low-Re 2-layer k-e model
+ TMODEL=3;KELIN=1;TLSC=EP
+ TURMOD(KEMODL-2L);STORE(REYN)
WHEN KW,2
+ TEXT(Wilcox 1988 k-w_2D ABRUPT PIPE EXPANS  :T213
+ MESG(Wilcox 1988 low-Re k-w model
+ TMODEL=4;TLSC=OMEG;KWMOD=T
+ TURMOD(KWMODL-LOWRE);STORE(REYT)
WHEN KWR,3
+ TEXT(Wilcox 2008 k-w_2D ABRUPT PIPE EXPANS  :T213
+ MESG(Wilcox 2008 low-Re k-w model
+ TMODEL=5;TLSC=OMEG;KWMOD=T
+ TURMOD(KWMODLR-LOWRE);FIINIT(FBP)=1.0
WHEN KWM,3
TEXT(Menter k-w_2D ABRUPT PIPE EXPANS  :T213
+ MESG(Menter 1992 k-w model
+ TURMOD(KWMENTER-LOWRE)
+ TMODEL=6;TLSC=OMEG;KWMOD=T
+ FIINIT(BF1)=1.0
WHEN KWS,3
TEXT(SST k-w__2D ABRUPT PIPE EXPANS  :T213
+ MESG(Menter k-w SST model
+ TURMOD(KWSST-LOWRE)
+ SOLUTN(V1,Y,Y,Y,P,P,P);SOLUTN(W1,Y,Y,Y,P,P,P)
+ TMODEL=7;TLSC=OMEG;KWMOD=T
+ FIINIT(BF1)=1.0;FIINIT(BF2)=1.0
+ RELAX(BF1,LINRLX,0.05);RELAX(BF2,LINRLX,1.0) 
ENDCASE

STORE(YPLS,STRS,KOND,ENUL,SPH1,SKIN,LEN1)
   ** Properties
RHO1=RHOIN;ENUL=(WOUT*DIAMO)/REYDO
SPHEAT=1.E3; CP1=SPHEAT  

   ** Mass inflow rate   
MFLOW = RHO1*VFLOW   
   ** Heat input; set Qin to raise temperature from 0 to 1.
DTRISE=1.0   
QIN=MFLOW*SPHEAT*DTRISE
   ** surface heat transfer area downstream of pipe expansion
AREAO=XULAST*RADO*PLDOWN
QWALL=QIN/AREAO
    Nusselt Number
   ** CNH1 = mass flow rate through high face of a cell
   ** Set CNH1 at NZ to CNH1 at NZ-1
(STORED of CNH1 is CNH1[,,-1] with if(IZ.EQ.NZ))
  ** Declaration of auxiliary In-Form variables: TSUM and ASUM
(MAKE TSUM is 0.)
(MAKE ASUM is 0.)
 
PATCH(PATCH1,CELL,1,NX,1,NY,1,NZ,1,LSTEP) ! One PATCH per domain
    ** Sum mass flow rate*specific heat for each IZ slab
(STORE1 ASUM at PATCH1 is SSUM(CNH1*CP1)!IMAT<100)
    ** Sum energy flux for each IZ slab
(STORE1 TSUM at PATCH1 is SSUM(CNH1*CP1*TEM1)!IMAT<100)
    ** Bulk temperature for each IZ slab
(STORED TAVE at PATCH1 is TSUM/ASUM)
    ** Compute wall temperature 
(STORED TWAL at HEATIN is TEM1+:QWALL:*0.5*DYV/KOND)
    ** Heat transfer coefficient in (W/m^2.degC)
(STORED HTCB at HEATIN is :QWALL:/(TWAL-TAVE+TINY))
    ** Nusselt number
(STORED NUSS at HEATIN is (HTCB*:DIAMO:/KOND))	
    ** Scaled Nusselt number
(STORED NUSC at HEATIN is NUSS/:NUDB:)	
    GROUP 8. Terms (in differential equations) & devices
TERMS(TEM1,N,P,P,P,P,P)
    GROUP 9. Properties of the medium (or media)
PRT(TEM1)=0.86;PRNDTL(TEM1)=PRLAM
    GROUP 11. Initialization of variable or porosity fields
STORE(PRPS)
PATCH(STEP,INIVAL,1,NX,(NYS+1),NY,1,NZS,1,1)
INIT(STEP,PRPS,0.,198)
EGWF=T
     ** Initial values
FIINIT(W1)=WIN;FIINIT(V1)=0.0;FIINIT(P1)=1.3E-4
FIINIT(KE)=TKEIN;FIINIT(TEM1)=TIN
IF(KWMOD) THEN
+ FIINIT(OMEG)=OMEGIN
ELSE
+ FIINIT(EP)=EPSIN
ENDIF
    GROUP 13. Boundary conditions and special sources
     ** Inlet
INLET(INLET,LOW,1,NX,1,NYS,1,1,1,1);VALUE(INLET,P1,RHO1*WIN)
VALUE(INLET,W1,WIN);VALUE(INLET,KE,TKEIN);VALUE(INLET,TEM1,TIN)
IF(KWMOD) THEN
+VALUE(INLET,OMEG,OMEGIN)
ELSE
+VALUE(INLET,EP,EPSIN)
ENDIF
     ** Exit
PATCH(OUTLET,HIGH,1,NX,1,NY,NZ,NZ,1,1);COVAL(OUTLET,P1,1.0E5,0.0)
COVAL(OUTLET,TEM1,ONLYMS,SAME)
     ** N-wall
WALL(WFUNNORT,NORTH,1,NX,NY,NY,NZS+1,NZ,1,1)

PATCH(HEATIN,NORTH,1,NX,NY,NY,NZS+1,NZ,1,1)
COVAL(HEATIN,TEM1,FIXFLU,QWALL)
    GROUP 15. Termination of sweeps
LSWEEP=5000;TSTSWP=-1
    GROUP 17. Under-relaxation devices
REAL(DTF,DTFKE);DTF=0.5*ZWLAST/NZ/WIN
RELAX(P1,LINRLX,1.0); RELAX(W1,FALSDT,DTF); RELAX(V1,FALSDT,DTF)
IF(KWMOD) THEN
+ RELAX(W1,FALSDT,DTF); RELAX(V1,FALSDT,DTF)
+ RELAX(KE,FALSDT,DTF/4.); RELAX(:TLSC:,FALSDT,DTF/4.)
  ** optional output of OMEG residuals and corrections
(STORED OF OMCR is CORR(OMEG) with IMAT<100)
(STORED OF OMRS is RESI(OMEG) with IMAT<100)
ELSE
+ KELIN=1
+ RELAX(KE,FALSDT,DTF/10.); RELAX(:TLSC:,FALSDT,DTF/10.)
  ** limit turbulence levels produced during course of 
     convergence to prevent divergence of LB & CK k-e models
+ KEMAX=(2.*WOUT)**2;EPMAX=WSTAR**4/ENUL
+ VARMAX(KE)=KEMAX;VARMAX(EP)=3.*EPMAX
  ** optional output of EP residuals and corrections
(STORED OF EPCR is CORR(EP) with IMAT<100)
(STORED OF EPRS is RESI(EP) with IMAT<100)
ENDIF
RELAX(TEM1,FALSDT,1.0)
    GROUP 18. Limits on variables or increments to them
LITER(P1)=50;LITER(TEM1)=50	  
VARMAX(TEM1)=1.E9;VARMIN(TEM1)=TIN-1.0;VARMIN(ENUT)=1.E-10
VARMAX(V1)=5.*WIN;VARMIN(V1)=-5.*WIN
VARMAX(W1)=5.*WIN;VARMIN(W1)=-5.*WIN
    GROUP 21. Print-out of variables
    GROUP 22. Monitor print-out
IZMON=NZS+2;IPLTL=2000;IYMON=NYS+2
ITABL=3;NPLT=50;NPRMON=10000;WALPRN=F
  GROUP 23. Field print-out and plot control
OUTPUT(ENUT,P,P,P,P,Y,Y)  
   ** Output Hot & Cold wall Nusselt numbers to .csv file	
PATCH(NUWALL,PROFIL,1,NX,NY,NY,NZS+1,NZ,1,LSTEP)
COVAL(NUWALL,NUSS,0.0,0.0);COVAL(NUWALL,NUSC,0.0,0.0)
COVAL(NUWALL,YPLS,0.0,0.0)
   ** Activate all 3 convergence history plots
SPEDAT(SET,GXMONI,PLOTALL,L,T)

DISTIL=T
CASE :CTURB: OF
WHEN LB,2
+EX(P1  )=6.138E+00;EX(V1  )=6.578E-02
+EX(W1  )=3.452E+00;EX(KE  )=1.237E+00
+EX(PRPS)=8.773E-01;EX(NUSC)=4.514E-02
+EX(NUSS)=4.390E+00;EX(HTCB)=2.462E-01
+EX(TWAL)=1.564E-01;EX(TAVE)=1.775E+01
+EX(CNH1)=2.008E-04;EX(SPH1)=8.773E+02
+EX(ENUL)=3.445E-05;EX(KOND)=4.921E-02
+EX(STRS)=9.007E-04
+EX(REYN)=3.802E+03;EX(LTLS)=2.731E-02
+EX(WDIS)=1.313E-01;EX(ENUT)=2.121E-02
+EX(TEM1)=1.792E+01;EX(EPRS)=7.718E-07
+EX(EPCR)=5.023E-05;EX(EP  )=2.182E+01
+EX(EPRS)=7.105E-07;EX(EPCR)=4.402E-05
+EX(LEN1)=3.899E-02
+EX(SKIN)=1.806E-01;EX(YPLS)=2.282E-02
+EX(TDIF)=3.882E-01;EX(CP  )=1.003E-01
WHEN CK,2
+EX(P1  )=6.948E+00;EX(V1  )=5.521E-02
+EX(W1  )=3.651E+00;EX(KE  )=1.023E+00
+EX(EP  )=2.000E+01;EX(PRPS)=8.773E-01
+EX(TWAL)=1.591E-01;EX(TAVE)=1.775E+01
+EX(CNH1)=2.044E-04
+EX(SPH1)=8.773E+02;EX(ENUL)=3.445E-05
+EX(KOND)=4.921E-02
+EX(YPLS)=2.048E-02;EX(REYN)=3.637E+03
+EX(LTLS)=2.731E-02;EX(WDIS)=1.313E-01
+EX(TEM1)=1.802E+01;EX(TDIF)=4.793E-01
+EX(NUSC)=4.028E-02;EX(NUSS)=3.917E+00
+EX(HTCB)=2.197E-01;EX(STRS)=7.984E-04
+EX(EPRS)=1.621E-06;EX(EPCR)=9.204E-05 
+EX(LEN1)=3.308E-02
+EX(SKIN)=6.523E-01;EX(ENUT)=1.659E-02 
+EX(CP  )=9.382E-02 
WHEN 2L,2
+EX(P1  )=6.415E+00;EX(V1  )=7.046E-02
+EX(W1  )=3.480E+00;EX(KE  )=1.205E+00
+EX(EP  )=2.473E+01;EX(PRPS)=8.773E-01
+EX(NUSC)=1.305E-02;EX(NUSS)=1.269E+00
+EX(HTCB)=7.117E-02;EX(TWAL)=1.724E-01
+EX(TAVE)=1.775E+01;EX(CNH1)=2.049E-04
+EX(SPH1)=8.773E+02;EX(ENUL)=3.445E-05
+EX(KOND)=4.921E-02;EX(CP  )=1.012E-01
+EX(STRS)=8.089E-04;EX(YPLS)=1.928E-02
+EX(REYN)=3.812E+03;EX(LTLS)=2.731E-02
+EX(WDIS)=1.313E-01;EX(ENUT)=1.789E-02
+EX(TEM1)=1.804E+01;EX(TDIF)=4.899E-01
+EX(EPRS)=1.046E+00;EX(EPCR)=2.691E-03
+EX(LEN1)=3.231E-02
+EX(SKIN)=6.179E-01
WHEN KW,2
+EX(P1  )=6.574E+00;EX(V1  )=6.217E-02
+EX(W1  )=3.500E+00;EX(KE  )=1.141E+00
+EX(EP  )=2.584E+01;EX(PRPS)=8.773E-01
+EX(NUSC)=1.451E-02;EX(NUSS)=1.411E+00
+EX(HTCB)=7.912E-02;EX(TWAL)=1.733E-01
+EX(TAVE)=1.775E+01;EX(SPH1)=8.773E+02
+EX(ENUL)=3.445E-05;EX(KOND)=4.921E-02
+EX(STRS)=6.844E-04
+EX(YPLS)=1.587E-02;EX(OMEG)=5.470E+02
+EX(TEM1)=1.808E+01;EX(TDIF)=5.381E-01
+EX(CNH1)=2.034E-04
+EX(OMRS)=2.121E-06;EX(OMCR)=2.725E-04 
+EX(LEN1)=3.578E-02
+EX(SKIN)=1.676E+00;EX(REYT)=4.551E+02 
+EX(ENUT)=1.769E-02;EX(CP  )=9.825E-02 
WHEN KWR,3
+EX(P1  )=7.792E+00;EX(V1  )=4.744E-02
+EX(W1  )=3.797E+00;EX(KE  )=8.805E-01
+EX(EP  )=2.037E+01;EX(PRPS)=8.773E-01
+EX(NUSC)=1.219E-02;EX(NUSS)=1.186E+00
+EX(HTCB)=6.650E-02;EX(TWAL)=1.779E-01
+EX(TAVE)=1.775E+01;EX(SPH1)=8.773E+02
+EX(ENUL)=3.445E-05;EX(KOND)=4.921E-02
+EX(STRS)=6.314E-04
+EX(YPLS)=1.459E-02;EX(DWDZ)=1.031E+00
+EX(DWDY)=8.636E+01;EX(DVDZ)=1.720E-01
+EX(DVDY)=9.311E-01;EX(DUDX)=2.371E-01
+EX(GEN1)=1.627E+05;EX(FBP )=8.064E-01
+EX(OMEG)=6.081E+02;EX(CP  )=8.539E-02
+EX(TEM1)=1.825E+01;EX(CNH1)=2.101E-04
+EX(OMRS)=1.832E-06;EX(OMCR)=1.146E-03
+EX(LEN1)=2.301E-02
+EX(SKIN)=6.089E-01;EX(XWP )=1.190E-01
+EX(ENUT)=1.049E-02;EX(TDIF)=7.243E-01
WHEN KWM,3
+EX(P1  )=6.541E+00;EX(V1  )=6.218E-02
+EX(W1  )=3.488E+00;EX(KE  )=1.127E+00
+EX(EP  )=2.083E+01;EX(PRPS)=8.773E-01
+EX(NUSC)=1.468E-02;EX(NUSS)=1.428E+00
+EX(HTCB)=8.009E-02;EX(TWAL)=1.769E-01
+EX(TAVE)=1.775E+01;EX(SPH1)=8.773E+02
+EX(ENUL)=3.445E-05;EX(KOND)=4.921E-02
+EX(STRS)=6.675E-04
+EX(YPLS)=1.530E-02;EX(LTLS)=2.731E-02
+EX(WDIS)=1.313E-01;EX(BF1 )=6.856E-01
+EX(OMEG)=5.648E+02;EX(ENUT)=1.824E-02
+EX(TEM1)=1.814E+01;EX(TDIF)=5.959E-01
+EX(CNH1)=2.025E-04
+EX(OMRS)=1.815E-06;EX(OMCR)=3.229E-04
+EX(LEN1)=3.607E-02;EX(SKIN)=1.529E+00
+EX(CP  )=9.979E-02
WHEN KWS,3
+EX(P1  )=6.900E+00;EX(V1  )=5.786E-02 
+EX(W1  )=3.583E+00;EX(KE  )=1.087E+00 
+EX(EP  )=2.015E+01;EX(PRPS)=8.773E-01 
+EX(NUSC)=1.408E-02;EX(NUSS)=1.369E+00 
+EX(HTCB)=7.681E-02;EX(TWAL)=1.778E-01 
+EX(TAVE)=1.775E+01
+EX(SKIN)=1.659E+00;EX(SPH1)=8.773E+02 
+EX(ENUL)=3.445E-05;EX(KOND)=4.921E-02 
+EX(STRS)=6.137E-04 
+EX(YPLS)=1.447E-02;EX(LTLS)=2.731E-02 
+EX(WDIS)=1.313E-01;EX(GEN1)=5.891E+07 
+EX(BF2 )=7.881E-01;EX(BF1 )=6.633E-01 
+EX(OMEG)=5.686E+02;EX(ENUT)=1.600E-02 
+EX(TEM1)=1.818E+01;EX(TDIF)=6.343E-01 
+EX(CNH1)=2.053E-04;EX(LEN1)=3.439E-02
+EX(OMRS)=1.067E-06;EX(OMCR)=1.415E-04
+EX(LEN1)=3.346E-02;EX(CP  )=9.523E-02
ENDCASE

STOP