function [XMEAS,YP,YY] = tefunc(t,idv,u,x,flag)

global XMEAS XMV YP0 NN YY IDV
     global  UCLR UCVR UTLR UTVR...
     XLR XVR ETR ESR...
     TCR TKR DLR ...
     VLR VVR VTR ...
     PTR PPR...
     CRXR RR RH ...
     FWR TWR QUR HWR UAR...
     UCLS UCVS UTLS UTVS...
     XLS XVS ETS ESS...
     TCS TKS DLS...
     VLS VVS VTS...
     PTS PPS...
     FWS TWS QUS HWS...
     UCLC UTLC XLC...
     ETC ESC TCC DLC...
     VLC VTC QUC...
     UCVV UTVV XVV...
     ETV ESV TCV TKV...
     VTV PTV...
     VCV VRNG VTAU...
     FTM...
     FCM XST XMWS...
     HST TST SFR...
     CPFLMX CPPRMX CPDH...
     TCWR TCWS...
     HTR AGSP...
     XDEL XNS...
     TGAS TPROD VST IVST;
	 global...
     ADIST...
     BDIST...
     CDIST...
     DDIST...
     TLAST...
     TNEXT...
     HSPAN...
     HZERO...
     SSPAN...
     SZERO...
     SPSPAN...
     IDVWLK;
     global...
     AVP BVP CVP...
     AH BH CH...
     AG BG CG...
     AV...
     AD BD CD...
     XMW;
global G
  NN = 50;
  TIME = t;
  YY=x;
  XMV=u;

% check for initialization 
if flag==1
	teinit;
else
end

for I=1:20
      if idv(I,1)>0
		  IDV(I,1)=1;
	  else
		  IDV(I,1)=0;
	  end
end
  


	  
%     TEFUNC(NN,TIME,YY,YP);
%     TEFUNCVAL = TEFUNC(NN,TIME,YY,YP)
%
% Substitution of TEFUNC for call statement. MWB
%
%                Tennessee Eastman Process Control Test Problem
% 
%                     James J. Downs and Ernest F. Vogel
% 
%                   Process and Control Systems Engineering
%                        Tennessee Eastman Company
%                               P.O. Box 511
%                           Kingsport,TN  37662
% 
% ************************************************************************
%
% 						Re-Written in Matlab 5.0 
%									by	
%								Martin Braun
%
%			Department of Chemical, Bio and Materials Engineering
%						Arizona State University
%
% ************************************************************************
%   Reference:
%     "A Plant-Wide Industrial Process Control Problem"
%     Presented at the AIChE 1990 Annual Meeting
%     Industrial Challenge Problems in Process Control,Paper #24a
%     Chicago,Illinois,November 14,1990
% 
%  Revised 4-4-91 to correct error in documentation of manipulated variables
% 
%  Subroutines:
% 
%     TEFUNC - Function evaluator to be called by integrator
%     TEINIT - Initialization
%     TESUBi - Utility subroutines, i=1,2,..,8
% 
% 
%   The process simulation has 50 states (NN=50).
%   Differences between the code and its description in the paper:
% 
%   1.  Subroutine TEINIT has TIME in the argument list.  TEINIT sets TIME
%       to zero.
% 
%   2.  There are 8 utility subroutines (TESUBi) rather than 5.
% 
%   3.  Process disturbances 14 through 20 do NOT need to be used in
%       conjunction with another disturbance as stated in the paper.  All
%       disturbances can be used alone or in any combination.
%  
%   Manipulated Variables
% 
%     XMV(1,1)     D Feed Flow (stream 2)            (Corrected Order)
%     XMV(1,2)     E Feed Flow (stream 3)            (Corrected Order)
%     XMV(1,3)     A Feed Flow (stream 1)            (Corrected Order)
%     XMV(1,4)     A and C Feed Flow (stream 4)
%     XMV(1,5)     Compressor Recycle Valve
%     XMV(1,6)     Purge Valve (stream 9)
%     XMV(1,7)     Separator Pot Liquid Flow (stream 10)
%     XMV(1,8)     Stripper Liquid Product Flow (stream 11)
%     XMV(1,9)     Stripper Steam Valve
%     XMV(1,10)    Reactor Cooling Water Flow
%     XMV(1,11)    Condenser Cooling Water Flow
%     XMV(1,12)    Agitator Speed
% 
%   Continuous Process Measurements
% 
%     XMEAS(1,1)   A Feed  (stream 1)                    kscmh
%     XMEAS(1,2)   D Feed  (stream 2)                    kg/hr
%     XMEAS(1,3)   E Feed  (stream 3)                    kg/hr
%     XMEAS(1,4)   A and C Feed  (stream 4)              kscmh
%     XMEAS(1,5)   Recycle Flow  (stream 8)              kscmh
%     XMEAS(1,6)   Reactor Feed Rate  (stream 6)         kscmh
%     XMEAS(1,7)   Reactor Pressure                      kPa gauge
%     XMEAS(1,8)   Reactor Level                         %
%     XMEAS(1,9)   Reactor Temperature                   Deg C
%     XMEAS(1,10)  Purge Rate (stream 9)                 kscmh
%     XMEAS(1,11)  Product Sep Temp                      Deg C
%     XMEAS(1,12)  Product Sep Level                     %
%     XMEAS(1,13)  Prod Sep Pressure                     kPa gauge
%     XMEAS(1,14)  Prod Sep Underflow (stream 10)        m3/hr
%     XMEAS(1,15)  Stripper Level                        %
%     XMEAS(1,16)  Stripper Pressure                     kPa gauge
%     XMEAS(1,17)  Stripper Underflow (stream 11)        m3/hr
%     XMEAS(1,18)  Stripper Temperature                  Deg C
%     XMEAS(1,19)  Stripper Steam Flow                   kg/hr
%     XMEAS(1,20)  Compressor Work                       kW
%     XMEAS(1,21)  Reactor Cooling Water Outlet Temp     Deg C
%     XMEAS(1,22)  Separator Cooling Water Outlet Temp   Deg C
% 
%   Sampled Process Measurements
% 
%     Reactor Feed Analysis (Stream 6)
%         Sampling Frequency = 0.1 hr
%         Dead Time = 0.1 hr
%         Mole %
%     XMEAS(1,23)   Component A
%     XMEAS(1,24)   Component B
%     XMEAS(1,25)   Component C
%     XMEAS(1,26)   Component D
%     XMEAS(1,27)   Component E
%     XMEAS(1,28)   Component F
% 
%     Purge Gas Analysis (Stream 9)
%         Sampling Frequency = 0.1 hr
%         Dead Time = 0.1 hr
%         Mole %
%     XMEAS(1,29)   Component A
%     XMEAS(1,30)   Component B
%     XMEAS(1,31)   Component C
%     XMEAS(1,32)   Component D
%     XMEAS(1,33)   Component E
%     XMEAS(1,34)   Component F
%     XMEAS(1,35)   Component G
%     XMEAS(1,36)   Component H
% 
%     Product Analysis (Stream 11)
%         Sampling Frequency = 0.25 hr
%         Dead Time = 0.25 hr
%         Mole %
%     XMEAS(1,37)   Component D
%     XMEAS(1,38)   Component E
%     XMEAS(1,39)   Component F
%     XMEAS(1,40)   Component G
%     XMEAS(1,41)   Component H
% 
%   Process Disturbances
% 
%     IDV(1,1)   A/C Feed Ratio, B Composition Constant (Stream 4)          Step
%     IDV(1,2)   B Composition, A/C Ratio Constant (Stream 4)               Step
%     IDV(1,3)   D Feed Temperature (Stream 2)                              Step
%     IDV(1,4)   Reactor Cooling Water Inlet Temperature                    Step
%     IDV(1,5)   Condenser Cooling Water Inlet Temperature                  Step
%     IDV(1,6)   A Feed Loss (Stream 1)                                     Step
%     IDV(1,7)   C Header Pressure Loss - Reduced Availability (Stream 4)   Step
%     IDV(1,8)   A, B, C Feed Composition (Stream 4)            Random Variation
%     IDV(1,9)   D Feed Temperature (Stream 2)                  Random Variation
%     IDV(1,10)  C Feed Temperature (Stream 4)                  Random Variation
%     IDV(1,11)  Reactor Cooling Water Inlet Temperature        Random Variation
%     IDV(1,12)  Condenser Cooling Water Inlet Temperature      Random Variation
%     IDV(1,13)  Reaction Kinetics                                    Slow Drift
%     IDV(1,14)  Reactor Cooling Water Valve                            Sticking
%     IDV(1,15)  Condenser Cooling Water Valve                          Sticking
%     IDV(1,16)  Unknown
%     IDV(1,17)  Unknown
%     IDV(1,18)  Unknown
%     IDV(1,19)  Unknown
%     IDV(1,20)  Unknown
% 
%
% =============================================================================
%        Function Evaluator
% 
%          Inputs:
% 
%            NN   = Number of differential equations
%            Time = Current time(hrs)
%            YY   = Current state values
% 
%          Outputs:
% 
%            YP   = Current derivative values
% 

	
XMNS = 0; % Initialization by MWB
% 	NOTE: I have included isd in the /PV/ common.  This is set
% 		non-zero when the process is shutting down.

% check for disturbance	 
for I = 1:20
	if IDV(I,1) > 0
    	IDV(I,1) = 1;
	else
		IDV(I,1) = 0;
	end
end

IDVWLK(1,1)=IDV(8,1);
IDVWLK(2,1)=IDV(8,1);
IDVWLK(3,1)=IDV(9,1);
IDVWLK(4,1)=IDV(10,1);
IDVWLK(5,1)=IDV(11,1);
IDVWLK(6,1)=IDV(12,1);
IDVWLK(7,1)=IDV(13,1);
IDVWLK(8,1)=IDV(13,1);
IDVWLK(9,1)=IDV(16,1);
IDVWLK(10,1)=IDV(17,1);
IDVWLK(11,1)=IDV(18,1);
IDVWLK(12,1)=IDV(20,1);
for I = 1:9
	if TIME >= TNEXT(I,1);
	HWLK=TNEXT(I,1)-TLAST(I,1);
    SWLK=ADIST(I,1)+HWLK*(BDIST(I,1)+HWLK*(CDIST(I,1)+HWLK*DDIST(I,1)));
    SPWLK=BDIST(I,1)+HWLK*(2*CDIST(I,1)+3*HWLK*DDIST(I,1));
    TLAST(I,1)=TNEXT(I,1);
	% TESUB5(SWLK,SPWLK,ADIST(I,1),BDIST(I,1),CDIST(I,1),DDIST(I,1),...
	% TLAST(I,1),TNEXT(I,1),HSPAN(I,1),HZERO(I,1),SSPAN(I,1),SZERO(I,1),...
	% SPSPAN(I,1),IDVWLK(I,1));
	% TESUB5(S,SP,ADIST,BDIST,CDIST,DDIST,
	% TLAST,TNEXT,HSPAN,HZERO,SSPAN,SZERO,...
	% SPSPAN,IDVFLAG) "VAL" added to subroutine variables.  "I" changed to "It".
	% Substitution of TESUB5 for call statement. MWB
	S5 = SWLK;
	SP5 = SPWLK;
	ADISTVAL = ADIST(I,1);
	BDISTVAL = BDIST(I,1);
	CDISTVAL = CDIST(I,1);
	DDISTVAL = DDIST(I,1);
	TLASTVAL = TLAST(I,1);
	TNEXTVAL = TNEXT(I,1);
	HSPANVAL = HSPAN(I,1);
	HZEROVAL = HZERO(I,1);
	SSPANVAL = SSPAN(I,1);
	SZEROVAL = SZERO(I,1);
	SPSPANVAL = SPSPAN(I,1);
	IDVFLAGVAL = IDVWLK(I,1);
	It = -1;
    H5 = HSPANVAL * TESUB7(It) + HZEROVAL;
    S1 = SSPANVAL * TESUB7(It) * IDVFLAGVAL + SZEROVAL;
    S1P = SPSPANVAL * TESUB7(It) * IDVFLAGVAL;
    ADISTVAL = S5;
    BDISTVAL = SP5;
    CDISTVAL = (3 * (S1 - S5) - H5 * (S1P + 2 * SP5)) / H5^2;
    DDISTVAL = (2 * (S5 - S1) + H5 * (S1P + SP5)) / H5^3;
    TNEXTVAL = TLASTVAL + H5;
	ADIST(I,1) = ADISTVAL; 
	BDIST(I,1) = BDISTVAL;
	CDIST(I,1) = CDISTVAL;
	DDIST(I,1) = DDISTVAL;
	TNEXT(I,1) = TNEXTVAL;  
	%
	% End of TESUB5
	%
	end
end

for I = 10:12
    if TIME >= TNEXT(I,1)
    HWLK = TNEXT(I,1) - TLAST(I,1);
    SWLK = ADIST(I,1) + HWLK * (BDIST(I,1) + HWLK * (CDIST(I,1) + HWLK * DDIST(I,1)));
    SPWLK = BDIST(I,1) + HWLK * (2 * CDIST(I,1) + 3 * HWLK * DDIST(I,1));
    TLAST(I,1) = TNEXT(I,1);
    	if SWLK > 0.1
      	ADIST(I,1) = SWLK;
      	BDIST(I,1) = SPWLK;
      	CDIST(I,1) = -(3 * SWLK + 0.2 * SPWLK) / 0.01;
      	DDIST(I,1) = (2 * SWLK + 0.1 * SPWLK) / 0.001;
      	TNEXT(I,1) = TLAST(I,1) + 0.1;
      	else
      	ISD = -1;
      	HWLK = HSPAN(I,1) * TESUB7(ISD) + HZERO(I,1);
     	ADIST(I,1) = 0;
      	BDIST(I,1) = 0;
      	CDIST(I,1) = (IDVWLK(I,1)) / HWLK^2;
      	DDIST(I,1) = 0;
      	TNEXT(I,1) = TLAST(I,1) + HWLK;
  		end
	end
end
if TIME == 0
	for I = 1:12
    ADIST(I,1)=SZERO(I,1);
    BDIST(I,1)=0;
    CDIST(I,1)=0;
    DDIST(I,1)=0;
    TLAST(I,1)=0.0;
    TNEXT(I,1)=0.1;
  	end
end

XST(1,4) = TESUB8(1,TIME) - IDV(1,1) * 0.03 - IDV(2,1) * 2.43719 * 10 ^ (-3);
XST(2,4) = TESUB8(2,TIME) + IDV(2,1) * 0.005;
XST(3,4) = 1 - XST(1,4) - XST(2,4);
TST(1,1) = TESUB8(3,TIME) + IDV(3,1) * 5;
TST(4,1) = TESUB8(4,TIME);
TCWR = TESUB8(5,TIME) + IDV(4,1) * 5;
TCWS = TESUB8(6,TIME) + IDV(5,1) * 5;
R1F = TESUB8(7,TIME);
R2F = TESUB8(8,TIME);
for I = 1:3
      UCVR(I,1)=YY(I,1);
      UCVS(I,1)=YY(I+9,1);
      UCLR(I,1)=0.0;
      UCLS(I,1)=0.0;
end
for I = 4:8
      UCLR(I,1) = YY(I,1);
      UCLS(I,1) = YY(I+9,1);
end
  
for I = 1:8
      UCLC(I,1) = YY(I+18,1);
      UCVV(I,1) = YY(I+27,1);
end
  
      ETR = YY(9,1);
      ETS = YY(18,1);
      ETC = YY(27,1);
      ETV = YY(36,1);
      TWR = YY(37,1);
      TWS = YY(38,1);

	  
%	  TWS shows an error here.  For the initialization,
%     it's value is  7.729698353000001e+01, when it should be
%     77.296983530000
	  
for I = 1:12
      VPOS(I,1) = YY(I+38,1);
end
      UTLR = 0.0;
      UTLS = 0.0;
      UTLC = 0.0;
      UTVV = 0.0;
	  
for I = 1:8
      UTLR = UTLR + UCLR(I,1);
      UTLS = UTLS + UCLS(I,1);
      UTLC = UTLC + UCLC(I,1);
      UTVV = UTVV + UCVV(I,1);
end
for I = 1:8
      XLR(I,1) = UCLR(I,1)/UTLR;
      XLS(I,1) = UCLS(I,1)/UTLS;
      XLC(I,1) = UCLC(I,1)/UTLC;
      XVV(I,1) = UCVV(I,1)/UTVV;
end
  
      ESR = ETR/UTLR;
      ESS = ETS/UTLS;
      ESC = ETC/UTLC;
      ESV = ETV/UTVV;
% Additional precision errors show up in XLR XLS XLC XVV ESR ESS ESC ESV 

%     TESUB2(XLR,TCR,ESR,0);
%	  TESUB2(Z,T,H,ITY)
% Substitution of TESUB2 for call statement. MWB
% TESUB2 appears to update T and only T.
Z = XLR;
T = TCR;
H = ESR;
ITY = 0;
	  TIN=T;
for J = 1:100
	%      TESUB1(Z,T,HTEST,ITY);
	%	   TESUB1VAL = TESUB1(Z,T,H,ITY)	
	% Substitution of TESUB1 for call statement MWB
	%
	if ITY == 0
    	HTEST = 0;
		for I = 1:8
      		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
      		HI = 1.8 * HI;
      		HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	else
      	HTEST = 0;
		for I = 1:8 
      		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
      		HI = 1.8 * HI;
      		HI = HI + AV(I,1);
      		HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	end
	if ITY == 2
      	R1 = 3.57696 / 1000000;
     	HTEST = HTEST - R1 * (T + 273.15);
  	end
	%
	% End of TESUB1 sub
	%
	
    ERR = HTEST - H;
	
	% 		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	%		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	% Substitution of TESUB3 for call statement MWB
	%
	if ITY == 0
      	DH = 0.0;
    	for I = 1:8
      		DHI = AH(I,1) + BH(I,1) * T + CH(I,1) * T^2;
      		DHI = 1.8 * DHI;
      		DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	else
      	DH=0.0;
    	for I = 1:8
      		DHI = AG(I,1) + BG(I,1) * T + CG(I,1) * T^2;
      		DHI = 1.8 * DHI;
      		DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	end
	if ITY == 2
      	R3 = 3.57696/1000000;
     	DH = DH - R3;
 	end
	%
	% End of TESUB3 substitution.
	%
    DT = -ERR / DH;
    T = T + DT;
  	if abs(DT) < 1*10^(-12)
    	break
  	end
  	if J == 100
  		T = TIN;
  	end
end
TCR = T;
% End of TESUB2.

      TKR = TCR+273.15;

%     TESUB2(XLS,TCS,ESS,0);
%	  TESUB2(Z,T,H,ITY)
% Substitution of TESUB2 for call statement. MWB
% TESUB2 appears to update T and only T.
Z = XLS;
T = TCS;
H = ESS;
ITY = 0;
TIN=T;
for J = 1:100
	%      TESUB1(Z,T,HTEST,ITY);
	%	   TESUB1VAL = TESUB1(Z,T,H,ITY)	
	% Substitution of TESUB1 for call statement MWB
	%	
	if ITY == 0
		HTEST = 0;
		for I = 1:8
			HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	else
		HTEST = 0;
		for I = 1:8 
			HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HI = HI + AV(I,1);
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	end
	if ITY == 2
		R1 = 3.57696 / 1000000;
		HTEST = HTEST - R1 * (T + 273.15);
	end
	%
	% End of TESUB1 sub
	%
 	ERR = HTEST - H;
	% 		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	%		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	% Substitution of TESUB3 for call statement MWB
	%
	if ITY == 0
   		DH = 0.0;
   		for I = 1:8
   			DHI = AH(I,1) + BH(I,1) * T + CH(I,1) * T^2;
   			DHI = 1.8 * DHI;
   			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	else
   		DH=0.0;
   		for I = 1:8
   			DHI = AG(I,1) + BG(I,1) * T + CG(I,1) * T^2;
   			DHI = 1.8 * DHI;
   			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	end
	if ITY == 2
    	R3 = 3.57696/1000000;
    	DH = DH - R3;
  	end
	%
	% End of TESUB3 substitution.
	%
    DT = -ERR / DH;
    T = T + DT;
	if abs(DT) < 1*10^(-12)
    	break
  	end
	if J == 100
  		T = TIN;
  	end
end
TCS = T;
% End of TESUB2.

      TKS = TCS+273.15;
      
%     TESUB2(XLC,TCC,ESC,0);
%	  TESUB2(Z,T,H,ITY)
% Substitution of TESUB2 for call statement. MWB
% TESUB2 appears to update T and only T.
Z = XLC;
T = TCC;
H = ESC;
ITY = 0;
TIN=T;
for J = 1:100
	%      TESUB1(Z,T,HTEST,ITY);
	%	   TESUB1VAL = TESUB1(Z,T,H,ITY)	
	% Substitution of TESUB1 for call statement MWB
	%
	if ITY == 0
		HTEST = 0;
		for I = 1:8
			HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	else
		HTEST = 0;
		for I = 1:8 
			HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HI = HI + AV(I,1);
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	end
	if ITY == 2
		R1 = 3.57696 / 1000000;
		HTEST = HTEST - R1 * (T + 273.15);
	end
	%
	% End of TESUB1 sub
	%
	ERR = HTEST - H;
	% 		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	%		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	% Substitution of TESUB3 for call statement MWB
	%
	if ITY == 0
		DH = 0.0;
    	for I = 1:8
			DHI = AH(I,1) + BH(I,1) * T + CH(I,1) * T^2;
			DHI = 1.8 * DHI;
			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	else
		DH=0.0;
		for I = 1:8
			DHI = AG(I,1) + BG(I,1) * T + CG(I,1) * T^2;
			DHI = 1.8 * DHI;
			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	end
	if ITY == 2
		R3 = 3.57696/1000000;
		DH = DH - R3;
	end
	%
	% End of TESUB3 substitution.
	%
	DT = -ERR / DH;
	T = T + DT;
	if abs(DT) < 1*10^(-12)
		break
	end
	if J == 100
		T = TIN;
	end
end
TCC = T;
% End of TESUB2.
      
%     TESUB2(XVV,TCV,ESV,2);
%	  TESUB2(Z,T,H,ITY)
% Substitution of TESUB2 for call statement. MWB
% TESUB2 appears to update T and only T.
Z = XVV;
T = TCV;
H = ESV;
ITY = 2;
TIN=T;
for J = 1:100
	%      TESUB1(Z,T,HTEST,ITY);
	%	   TESUB1VAL = TESUB1(Z,T,H,ITY)	
	% Substitution of TESUB1 for call statement MWB
	%
	if ITY == 0
		HTEST = 0;
		for I = 1:8
			HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	else
		HTEST = 0;
		for I = 1:8 
			HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HI = HI + AV(I,1);
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	end
	if ITY == 2
		R1 = 3.57696 / 1000000;
		HTEST = HTEST - R1 * (T + 273.15);
	end
	%
	% End of TESUB1 sub
	%
	ERR = HTEST - H;
	% 		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	%		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	% Substitution of TESUB3 for call statement MWB
	%
	if ITY == 0
		DH = 0.0;
    	for I = 1:8
			DHI = AH(I,1) + BH(I,1) * T + CH(I,1) * T^2;
			DHI = 1.8 * DHI;
			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	else
		DH=0.0;
    	for I = 1:8
			DHI = AG(I,1) + BG(I,1) * T + CG(I,1) * T^2;
			DHI = 1.8 * DHI;
			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	end
	if ITY == 2
		R3 = 3.57696/1000000;
		DH = DH - R3;
	end
	%
	% End of TESUB3 substitution.
	%
	DT = -ERR / DH;
	T = T + DT;
	if abs(DT) < 1*10^(-12)
		break
	end
	if J == 100
		T = TIN;
	end
end
TCV = T;
% End of TESUB2.
	  
TKV = TCV+273.15;
      
%	  TESUB4(XLR,TCR,DLR);
%	  TESUB4(X,T,R)
% Substitution of TESUB4 for call statement. MWB
X = XLR;
T4 = TCR;
R4 = DLR;
V = 0.0;
for I = 1:8
	V = V + X(I,1) * XMW(I,1) / (AD(I,1) + (BD(I,1) + CD(I,1) * T4) * T4);
end
R4 = 1.0 / V;
DLR = R4;
% End of TESUB4.

      
%	  TESUB4(XLS,TCS,DLS);
%	  TESUB4(X,T,R)
% Substitution of TESUB4 for call statement. MWB
X = XLS;
T4 = TCS;
R4 = DLS;
V = 0.0;
for I = 1:8
	V = V + X(I,1) * XMW(I,1) / (AD(I,1) + (BD(I,1) + CD(I,1) * T4) * T4);
end
R4 = 1.0 / V;
DLS = R4;
% End of TESUB4.
      
%	  TESUB4(XLC,TCC,DLC);
%	  TESUB4(X,T,R)
% Substitution of TESUB4 for call statement. MWB
X = XLC;
T4 = TCC;
R4 = DLC;
V = 0.0;
for I = 1:8
      V = V + X(I,1) * XMW(I,1) / (AD(I,1) + (BD(I,1) + CD(I,1) * T4) * T4);
end
R4 = 1.0 / V;
DLC = R4;
% End of TESUB4.
VLR = UTLR/DLR;
VLS = UTLS/DLS;
VLC = UTLC/DLC;
VVR = VTR-VLR;
VVS = VTS-VLS;
RG = 998.9;
PTR = 0.0;
PTS = 0.0;
for I = 1:3
      PPR(I,1) = UCVR(I,1)*RG*TKR/VVR;
      PTR = PTR+PPR(I,1);
      PPS(I,1) = UCVS(I,1)*RG*TKS/VVS;
      PTS = PTS+PPS(I,1);
  end
% Some precision error with PTS
for I = 4:8
      VPR=exp(AVP(I,1)+BVP(I,1)/(TCR+CVP(I,1)));
      PPR(I,1)=VPR*XLR(I,1);
      PTR=PTR+PPR(I,1);
      VPR=exp(AVP(I,1)+BVP(I,1)/(TCS+CVP(I,1)));
      PPS(I,1)=VPR*XLS(I,1);
      PTS=PTS+PPS(I,1);
  end
      PTV=UTVV*RG*TKV/VTV;
for I=1:8
      XVR(I,1)=PPR(I,1)/PTR;
      XVS(I,1)=PPS(I,1)/PTS;
end
UTVR=PTR*VVR/RG/TKR;
UTVS=PTS*VVS/RG/TKS;
for I=4:8
      UCVR(I,1)=UTVR*XVR(I,1);
      UCVS(I,1)=UTVS*XVS(I,1);
  end
  
% This is where the Matlab RR's show different values than the FORTRAN 77
% RR's.  The code below has been modified with the chop command to 
% elliminate the stray numbers that are created by the matlab code.
% Matlab:				FORTRAN 77:
% R1F
% 1						1.0000000000000
% R2F
% 1						1.0000000000000
% TKR
% 3.935500000017000e+02	393.55000000169
% RR(1,1)
% 3.181229946057000e-09	3.1812215872809D-09
% RR(2,1)
% 1.569591352140000e-10	1.5695892900646D-10
% RR(3,1)
% 7.437635682791000e-11	7.4376432806282D-11
% RR(4,1)
% 5.708298619084000e-11	5.7083044503356D-11
%
% NOTE:  I DID A COMPARISON OF THE EXP FUNCTIONS OF MATLAB
%        AND FORTRAN 77 AND THERE IS A CONSISTENT DISCREPANCY AT
%        10^-16 BETWEEN THE OUTPUT OF THE TWO LANGUAGES
%

%R1F = chop(R1F,13);
%R2F = chop(R2F,13);
%TKR = chop(TKR,13);
RR(1,1)=(exp(31.5859536-40000.0/1.987/TKR)*R1F);
RR(2,1)=(exp(3.00094014-20000.0/1.987/TKR)*R2F);
RR(3,1)=(exp(53.4060443-60000.0/1.987/TKR));
RR(4,1)=(RR(3,1)*0.767488334);

% This next section has the following output:
% Matlab:				FORTRAN 77:
% PPR(1,1)
% 5.717769731191162e+03	5717.7697311912     * good
% PPR(3,1)
% 4.159943934675613e+03 4159.9439346756		* good
% R1F
% 2.174415112221978e+04 21744.151122220		* good
% R2F
% 2.247652526817436e+01 22.476525268174		* good
% RR(1,1)
% 3.517133031511259e-01 0.35171237901398	
% RR(2,1)
% 2.861611337089543e-01 0.28616075776025

if PPR(1,1) > 0.0 & PPR(3,1) > 0.0
      R1F=PPR(1,1)^1.1544;
      R2F=PPR(3,1)^0.3735;
      RR(1,1)=RR(1,1)*R1F*R2F*PPR(4,1);
      RR(2,1)=RR(2,1)*R1F*R2F*PPR(5,1);
  else
	  
      RR(1,1)=0.0;
      RR(2,1)=0.0;
  end
%
% This next section has the following output:
% Matlab:				FORTRAN 77:
% RR(1,1)
% 2.516069042865577e+02 251.60624318195
% RR(2,1)
% 2.047125210635105e+02 204.71225211921
% RR(3,1)
% 1.134874248863408e+00 1.1348754081821
% RR(4,1)
% 5.281888235868223e-02 5.2818936315249D-02
%
RR(3,1)=RR(3,1)*PPR(1,1)*PPR(5,1);
RR(4,1)=RR(4,1)*PPR(1,1)*PPR(4,1);
for I=1:4
    RR(I,1)=RR(I,1)*VVR;
end
%
% This next section has the following output:
% Matlab:					FORTRAN 77:
% CRXR(1,1)
% -4.574542995989316e+02	-457.45337070934
% CRXR(3,1)
% -4.563194253500682e+02	-456.31849530116
% CRXR(4,1)
% -2.516861326100957e+02	-251.68547158642
% CRXR(5,1)
% -2.058473953123739e+02	-205.84712752739
% CRXR(6,1)
% 1.187693131222090e+00		1.1876943444973
% CRXR(7,1)
% 2.516069042865577e+02		251.60624318195
% CRXR(8,1)
% 2.047125210635105e+02		204.71225211921
% RH
% 2.759494513807820e+01		27.594886078737
CRXR(1,1)=-RR(1,1)-RR(2,1)-RR(3,1);
CRXR(3,1)=-RR(1,1)-RR(2,1);
CRXR(4,1)=-RR(1,1)-1.5*RR(4,1);
CRXR(5,1)=-RR(2,1)-RR(3,1);
CRXR(6,1)=RR(3,1)+RR(4,1);
CRXR(7,1)=RR(1,1);
CRXR(8,1)=RR(2,1);
RH=RR(1,1)*HTR(1,1)+RR(2,1)*HTR(2,1);
XMWS(1,1)=0.0;
XMWS(2,1)=0.0;
XMWS(6,1)=0.0;
XMWS(8,1)=0.0;
XMWS(9,1)=0.0;
XMWS(10,1)=0.0;
%
% XST & XMWS seem to be calculated correcly in the next 
% section with little error.
%
for I=1:8
	  XST(I,6)=XVV(I,1);
	  XST(I,8)=XVR(I,1);
	  XST(I,9)=XVS(I,1);
      XST(I,10)=XVS(I,1);
      XST(I,11)=XLS(I,1);
      XST(I,13)=XLC(I,1);
      XMWS(1,1)=XMWS(1,1)+XST(I,1)*XMW(I,1);
      XMWS(2,1)=XMWS(2,1)+XST(I,2)*XMW(I,1);
      XMWS(6,1)=XMWS(6,1)+XST(I,6)*XMW(I,1);
      XMWS(8,1)=XMWS(8,1)+XST(I,8)*XMW(I,1);
      XMWS(9,1)=XMWS(9,1)+XST(I,9)*XMW(I,1);
      XMWS(10,1)=XMWS(10,1)+XST(I,10)*XMW(I,1);
  end
TST(6,1)=TCV;
TST(8,1)=TCR;
TST(9,1)=TCS;
TST(10,1)=TCS;
TST(11,1)=TCS;
TST(13,1)=TCC;

% The HST's from the following substitutions for TESUB1 are 
% computed correctly.

%
%      TESUB1(XST(1,1),TST(1,1),HST(1,1),1);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,1);
T=TST(1,1);
H=HST(1,1);
ITY=1;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
	R1 = 3.57696 / 1000000;
	H = H - R1 * (T + 273.15);
end
HST(1,1) = H;
% end of TESUB1


%      TESUB1(XST(1,2),TST(2,1),HST(2,1),1);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,2);
T=TST(2,1);
H=HST(2,1);
ITY=1;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(2,1) = H;
% end of TESUB1

%      TESUB1(XST(1,3),TST(3,1),HST(3,1),1);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,3);
T=TST(3,1);
H=HST(3,1);
ITY=1;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(3,1) = H;
% end of TESUB1

%      TESUB1(XST(1,4),TST(4,1),HST(4,1),1);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,4);
T=TST(4,1);
H=HST(4,1);
ITY=1;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(4,1) = H;
% end of TESUB1

%      TESUB1(XST(1,6),TST(6,1),HST(6,1),1);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,6);
T=TST(6,1);
H=HST(6,1);
ITY=1;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(6,1) = H;
% end of TESUB1

%      TESUB1(XST(1,8),TST(8,1),HST(8,1),1);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,8);
T=TST(8,1);
H=HST(8,1);
ITY=1;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(8,1) = H;
% end of TESUB1

%      TESUB1(XST(1,9),TST(9,1),HST(9,1),1);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,9);
T=TST(9,1);
H=HST(9,1);
ITY=1;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(9,1) = H;
% end of TESUB1

      HST(10,1)=HST(9,1);
	  
%      TESUB1(XST(1,11),TST(11,1),HST(11,1),0);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,11);
T=TST(11,1);
H=HST(11,1);
ITY=0;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
 		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(11,1) = H;
% end of TESUB1

%      TESUB1(XST(1,13),TST(13,1),HST(13,1),0);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,13);
T=TST(13,1);
H=HST(13,1);
ITY=0;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(13,1) = H;
% end of TESUB1

% This next section is good.
FTM(1,1)=VPOS(1,1)*VRNG(1,1)/100.0;
FTM(2,1)=VPOS(2,1)*VRNG(2,1)/100.0;
FTM(3,1)=VPOS(3,1)*(1-IDV(6,1))*VRNG(3,1)/100.0;
FTM(4,1)=VPOS(4,1)*(1-IDV(7,1)*0.2)*VRNG(4,1)/100.0+1.E-10;
FTM(11,1)=VPOS(7,1)*VRNG(7,1)/100.0;
FTM(13,1)=VPOS(8,1)*VRNG(8,1)/100.0;
UAC=VPOS(9,1)*VRNG(9,1)*(1+TESUB8(9,TIME))/100.0;
FWR=VPOS(10,1)*VRNG(10,1)/100.0;
FWS=VPOS(11,1)*VRNG(11,1)/100.0;
AGSP=(VPOS(12,1)+150.0)/100.0;
DLP=PTV-PTR;

if DLP<0 
	DLP=0.0;
end
FLMS=1937.6*sqrt(DLP);
FTM(6,1)=FLMS/XMWS(6,1);
DLP=PTR-PTS;
if DLP<0 
	DLP=0.0;
end
FLMS=4574.21*sqrt(DLP)*(1-0.25*TESUB8(12,TIME));
FTM(8,1)=FLMS/XMWS(8,1);
DLP=PTS-760.0;
if DLP < 0 
	DLP=0.0;
end
FLMS=VPOS(6,1)*0.151169*sqrt(DLP);
FTM(10,1)=FLMS/XMWS(10,1);
PR=PTV/PTS;
if PR < 1
	PR=1.0;
end
if PR > CPPRMX
	PR=CPPRMX;
end
FLCOEF=CPFLMX/1.197;
FLMS=CPFLMX+FLCOEF*(1.0-PR^3);
CPDH=FLMS*(TCS+273.15)*1.8E-6*1.9872*(PTV-PTS)/(XMWS(9,1)*PTS);
DLP=PTV-PTS;
if DLP < 0.0
	DLP=0.0;
end
FLMS=FLMS-VPOS(5,1)*53.349*sqrt(DLP);
if FLMS < 1.E-3
	FLMS=1.E-3;
end
FTM(9,1)=FLMS/XMWS(9,1);
HST(9,1)=HST(9,1)+CPDH/FTM(9,1);
for I=1:8
      FCM(I,1)=XST(I,1)*FTM(1,1);
      FCM(I,2)=XST(I,2)*FTM(2,1);
      FCM(I,3)=XST(I,3)*FTM(3,1);
      FCM(I,4)=XST(I,4)*FTM(4,1);
      FCM(I,6)=XST(I,6)*FTM(6,1);
      FCM(I,8)=XST(I,8)*FTM(8,1);
      FCM(I,9)=XST(I,9)*FTM(9,1);
      FCM(I,10)=XST(I,10)*FTM(10,1);
      FCM(I,11)=XST(I,11)*FTM(11,1);
      FCM(I,13)=XST(I,13)*FTM(13,1);
  end

if FTM(11,1) > 0.1
	if TCC > 170
		TMPFAC=TCC-120.262;
	elseif TCC < 5.292
		TMPFAC=0.1;
	else
		TMPFAC=363.744/(177.-TCC)-2.22579488;
	end
	VOVRL=FTM(4,1)/FTM(11,1)*TMPFAC;
	SFR(4,1)=8.5010*VOVRL/(1.0+8.5010*VOVRL);
	SFR(5,1)=11.402*VOVRL/(1.0+11.402*VOVRL);
	SFR(6,1)=11.795*VOVRL/(1.0+11.795*VOVRL);
	SFR(7,1)=0.0480*VOVRL/(1.0+0.0480*VOVRL);
	SFR(8,1)=0.0242*VOVRL/(1.0+0.0242*VOVRL);
else
	SFR(4,1)=0.9999;
	SFR(5,1)=0.999;
	SFR(6,1)=0.999;
	SFR(7,1)=0.99;
	SFR(8,1)=0.98;
end
for I=1:8
      FIN(I,1)=0.0;
      FIN(I,1)=FIN(I,1)+FCM(I,4);
      FIN(I,1)=FIN(I,1)+FCM(I,11);
  end

FTM(5,1)=0.0;
FTM(12,1)=0.0;
for I=1:8
      FCM(I,5)=SFR(I,1)*FIN(I,1);
      FCM(I,12)=FIN(I,1)-FCM(I,5);
      FTM(5,1)=FTM(5,1)+FCM(I,5);
      FTM(12,1)=FTM(12,1)+FCM(I,12);
end
for I=1:8
      XST(I,5)=FCM(I,5)/FTM(5,1);
      XST(I,12)=FCM(I,12)/FTM(12,1);
end
TST(5,1)=TCC;
TST(12,1)=TCC;
%      TESUB1(XST(1,5),TST(5,1),HST(5,1),1);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,5);
T=TST(5,1);
H=HST(5,1);
ITY=1;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(5,1) = H;
% end of TESUB1

%      TESUB1(XST(1,12),TST(12,1),HST(12,1),0);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,12);
T=TST(12,1);
H=HST(12,1);
ITY=0;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
	R1 = 3.57696 / 1000000;
	H = H - R1 * (T + 273.15);
end
HST(12,1) = H;
% end of TESUB1

FTM(7,1)=FTM(6,1);
HST(7,1)=HST(6,1);
TST(7,1)=TST(6,1);

for I=1:8
	
      XST(I,7)=XST(I,6);
      FCM(I,7)=FCM(I,6);
end

if VLR/7.8 > 50.0
	UARLEV=1.0;
elseif VLR/7.8 < 10.0
    UARLEV=0.0;
else
      UARLEV=0.025*VLR/7.8-0.25;
  end
  
UAR=UARLEV*(-0.5*AGSP^2+2.75*AGSP-2.5)*855490.E-6;
QUR=UAR*(TWR-TCR)*(1-0.35*TESUB8(10,TIME));
UAS=0.404655*(1.0-1.0/(1.0+(FTM(8,1)/3528.73)^4));
QUS=UAS*(TWS-TST(8,1))*(1-0.25*TESUB8(11,TIME));
QUC=0;

if TCC < 100
	QUC=UAC*(100.0-TCC);
end
      XMEAS(1,1)=FTM(3,1)*0.359/35.3145;
      XMEAS(2,1)=FTM(1,1)*XMWS(1,1)*0.454;
      XMEAS(3,1)=FTM(2,1)*XMWS(2,1)*0.454;
      XMEAS(4,1)=FTM(4,1)*0.359/35.3145;
      XMEAS(5,1)=FTM(9,1)*0.359/35.3145;
      XMEAS(6,1)=FTM(6,1)*0.359/35.3145;
      XMEAS(7,1)=(PTR-760.0)/760.0*101.325;
      XMEAS(8,1)=(VLR-84.6)/666.7*100.0;
      XMEAS(9,1)=TCR;
      XMEAS(10,1)=FTM(10,1)*0.359/35.3145;
      XMEAS(11,1)=TCS;
      XMEAS(12,1)=(VLS-27.5)/290.0*100.0;
      XMEAS(13,1)=(PTS-760.0)/760.0*101.325;
      XMEAS(14,1)=FTM(11,1)/DLS/35.3145;
      XMEAS(15,1)=(VLC-78.25)/VTC*100.0;
      XMEAS(16,1)=(PTV-760.0)/760.0*101.325;
      XMEAS(17,1)=FTM(13,1)/DLC/35.3145;
      XMEAS(18,1)=TCC;
      XMEAS(19,1)=QUC*1.04E3*0.454;
      XMEAS(20,1)=CPDH*0.0003927E6;
      XMEAS(20,1)=CPDH*0.29307E3;
      XMEAS(21,1)=TWR;
      XMEAS(22,1)=TWS;
	 ISD=0;
if XMEAS(7,1) > 3000.0
	ISD=1;
	end
if VLR/35.3145 > 24.0
	ISD=1;
	end
if VLR/35.3145 < 2.0
	ISD=1;
	end
if XMEAS(9,1) > 175.0
	ISD=1;
	end
if VLS/35.3145 > 12.0
	ISD=1;
	end
if VLS/35.3145 < 1.0
	ISD=1;
	end
if VLC/35.3145 > 8.0
	ISD=1;
	end
if VLC/35.3145 < 1.0
	ISD=1;
	end
if TIME > 0.0 & ISD == 0
	for I=1:22
		%      TESUB6(XNS(I,1),XMNS)
		%      TESUB6(STD,X)
		% Substitution of TESUB6 for call statement. MWB
		STD = XNS(I,1);
		X = XMNS;
		X = 0;
		for Iter =1:12
			X = X + TESUB7(Iter);
		end
		X = (X - 6) * STD;
		XMNS = X;
		XMNS=0;%turn off noise
		% End of TESUB6.

		XMEAS(I,1)=XMEAS(I,1)+XMNS;
	end
end
      XCMP(23,1)=XST(1,7)*100.0;
      XCMP(24,1)=XST(2,7)*100.0;
      XCMP(25,1)=XST(3,7)*100.0;
      XCMP(26,1)=XST(4,7)*100.0;
      XCMP(27,1)=XST(5,7)*100.0;
      XCMP(28,1)=XST(6,7)*100.0;
      XCMP(29,1)=XST(1,10)*100.0;
      XCMP(30,1)=XST(2,10)*100.0;
      XCMP(31,1)=XST(3,10)*100.0;
      XCMP(32,1)=XST(4,10)*100.0;
      XCMP(33,1)=XST(5,10)*100.0;
      XCMP(34,1)=XST(6,10)*100.0;
      XCMP(35,1)=XST(7,10)*100.0;
      XCMP(36,1)=XST(8,10)*100.0;
      XCMP(37,1)=XST(4,13)*100.0;
      XCMP(38,1)=XST(5,13)*100.0;
      XCMP(39,1)=XST(6,13)*100.0;
      XCMP(40,1)=XST(7,13)*100.0;
      XCMP(41,1)=XST(8,13)*100.0;
	  
if TIME == 0
	for I=23:41
		XDEL(I,1)=XCMP(I,1);
		XMEAS(I,1)=XCMP(I,1);
	end

	TGAS=0.1;
	TPROD=0.25;
end
if TIME >= TGAS
	for I=23:36
		XMEAS(I,1)=XDEL(I,1);
		%      TESUB6(XNS(I,1),XMNS);
		%      TESUB6(STD,X)
		% Substitution of TESUB6 for call statement. MWB
		STD = XNS(I,1);
		X = XMNS;
		X = 0;
		for Iter=1:12
			X = X + TESUB7(Iter);
		end
		X = (X - 6) * STD;
		XMNS = X;
		% End of TESUB6
		XMNS=0;%turn off noise
		XMEAS(I,1)=XMEAS(I,1)+XMNS;
		XDEL(I,1)=XCMP(I,1);
	end
	TGAS=TGAS+0.1;
end
 
if TIME >= TPROD
	for I=37:41
		XMEAS(I,1)=XDEL(I,1);
      
		%      TESUB6(XNS(I,1),XMNS)
		%      TESUB6(STD,X)
		% Substitution of TESUB6 for call statement. MWB
		STD = XNS(I,1);
		X = XMNS;
		X = 0;
		for Iter=1:12
			X = X + TESUB7(Iter);
		end
		X = (X - 6) * STD;
		XMNS = X;
		XMNS=0;%turn off noise
		% End of TESUB6

		XMEAS(I,1)=XMEAS(I,1)+XMNS;
		XDEL(I,1)=XCMP(I,1);
	end
	TPROD=TPROD+0.25;
end
% 
for iter = 1:8
      YP(iter,1) = FCM(iter,7) - FCM(iter,8) + CRXR(iter,1);
      YP(iter+9,1) = FCM(iter,8) - FCM(iter,9) - FCM(iter,10) - FCM(iter,11);
      YP(iter+18,1) = FCM(iter,12) - FCM(iter,13);
      YP(iter+27,1) = FCM(iter,1) + FCM(iter,2) + FCM(iter,3) + FCM(iter,5) + FCM(iter,9) - FCM(iter,6);
  end
  
% 		Here is the "correct" version of the separator energy balance:
% 
% 	YP(18)=HST(8)*FTM(8)-(HST(9)*FTM(9)-cpdh)-HST(10)*FTM(10)-HST(11)*FTM(11)+QUS
% 
%   	Here is the original version
	  YP(9,1)=HST(7,1) * FTM(7,1) - HST(8) * FTM(8) + RH + QUR;
      YP(18,1)=HST(8,1)*FTM(8,1)-HST(9,1)*FTM(9,1)-HST(10,1)*FTM(10,1)-HST(11,1)*FTM(11,1)+QUS;
      YP(27,1)=HST(4,1)*FTM(4,1)+HST(11,1)*FTM(11,1)-HST(5,1)*FTM(5,1)-HST(13,1)*FTM(13,1)+QUC;
      YP(36,1)=HST(1,1)*FTM(1,1)+HST(2,1)*FTM(2,1)+HST(3,1)*FTM(3,1)+HST(5,1)*FTM(5,1)+HST(9,1)*FTM(9,1)-HST(6,1)*FTM(6,1);
      YP(37,1)=(FWR*500.53*(TCWR-TWR)-QUR*1000000/1.8)/HWR;
      YP(38,1)=(FWS*500.53*(TCWS-TWS)-QUS*1000000/1.8)/HWS;
      IVST(10,1)=IDV(14,1);
      IVST(11,1)=IDV(15,1);
      IVST(5,1)=IDV(19,1);
      IVST(7,1)=IDV(19,1);
      IVST(8,1)=IDV(19,1);
      IVST(9,1)=IDV(19,1);
for I=1:12
	if TIME == 0 | abs(VCV(I,1)-XMV(I,1)) > VST(I,1)*IVST(I,1)
		VCV(I,1)=XMV(I,1);
	end
	if VCV(I,1) < 0.0
		VCV(I,1)=0.0;
	end
	if VCV(I,1) > 100.0
		VCV(I,1)=100.0;
	end
	YP(I+38,1)=(VCV(I,1)-VPOS(I,1))/VTAU(I,1);
end
if TIME > 0.0 & ISD ~= 0
	for I=1:NN
		YP(I,1)=0.0;
	end
end
% This is the final output of the subroutine:
% Everything up to this point (except what has been
% documented in the code) runs very similar to the FORTRAN 77
% version. MWB 
%  7/12/98
%  Matlab:						FORTRAN77:
%  YY
% 	   1.040491389000000e+01     1.040491389000000e+01
%      4.363996017000000e+00     4.363996017000000e+00
%      7.570059737000000e+00     7.570059737000000e+00
%      4.230042431000000e-01     4.230042431000000e-01
%      2.415513437000000e+01     2.415513437000000e+01
%      2.942597645000000e+00     2.942597645000000e+00
%      1.543770655000000e+02     1.543770655000000e+02
%      1.591865960000000e+02     1.591865960000000e+02
%      2.808522723000000e+00     2.808522723000000e+00
%      6.375581199000000e+01     6.375581199000000e+01
%      2.674026066000000e+01     2.674026066000000e+01
%      4.638532432000000e+01     4.638532432000000e+01
%      2.464521543000000e-01     2.464521543000000e-01
%      1.520484404000000e+01     1.520484404000000e+01
%      1.852266172000000e+00     1.852266172000000e+00
%      5.244639459000000e+01     5.244639459000000e+01
%      4.120394008000000e+01     4.120394008000000e+01
%      5.699317760000000e-01     5.699317760000000e-01
%      4.306056376000000e-01     4.306056376000000e-01
%      7.990620078300001e-03     7.990620078300001e-03
%      9.056036089000000e-01     9.056036089000000e-01
%      1.605425821600000e-02     1.605425821600000e-02
%      7.509759687000001e-01     7.509759687000001e-01
%      8.858285595500000e-02     8.858285595500000e-02
%      4.827726193000000e+01     4.827726193000000e+01
%      3.938459028000000e+01     3.938459028000000e+01
%      3.755297257000000e-01     3.755297257000000e-01
%      1.077562698000000e+02     1.077562698000000e+02
%      2.977250546000000e+01     2.977250546000000e+01
%      8.832481135000000e+01     8.832481135000000e+01
%      2.303929507000000e+01     2.303929507000000e+01
%      6.285848794000000e+01     6.285848794000000e+01
%      5.546318688000000e+00     5.546318688000000e+00
%      1.192244772000000e+01     1.192244772000000e+01
%      5.555448243000000e+00     5.555448243000000e+00
%      9.218489761999999e-01     9.218489761999999e-01
%      9.459927549000000e+01     9.459927549000000e+01
%      7.729698353000001e+01     7.729698353000001e+01
%      6.305263039000000e+01     6.305263039000000e+01
%      5.397970677000000e+01     5.397970677000000e+01
%      2.464355755000000e+01     2.464355755000000e+01
%      6.130192144000000e+01     6.130192144000000e+01
%      2.221000000000000e+01     2.221000000000000e+01
%      4.006374673000000e+01     4.006374673000000e+01
%      3.810034370000000e+01     3.810034370000000e+01
%      4.653415582000000e+01     4.653415582000000e+01
%      4.744573456000000e+01     4.744573456000000e+01
%      4.110581288000000e+01     4.110581288000000e+01
%      1.811349055000000e+01     1.811349055000000e+01
%      5.000000000000000e+01     5.000000000000000e+01
%   YP:
%     -9.311351132055279e-04    -2.245521557142600e-06 <- Sig Dif
%     -3.680818849716161e-07    -3.680818849716200e-07
%     -9.313903218526320e-04    -1.341411575595000e-06 <- Sig Dif
%     -6.614342966599907e-04    -4.106221069832800e-07 <- Sig Dif
%     -2.689011108429895e-04    -1.116128657940900e-06 <- Sig Dif
%     -1.312917689233473e-06    -9.964246872051801e-08 <- Sig Dif
%      6.609202557683602e-04    -1.843536381329600e-07 <- Sig Dif
%      2.688889922808357e-04    -5.530856128643800e-08 <- Sig Dif
%      5.904527482059052e-05    -1.406638006074000e-08 <- Sig Dif
%     -3.216924339355387e-06    -3.216923774473900e-06
%     -1.320602298626739e-06    -1.320602072141200e-06
%     -2.292647105051060e-06    -2.292646652080100e-06
%     -1.172087313872083e-07    -1.172087100709300e-07
%     -1.564706110457337e-06    -1.564705769396800e-06
%     -2.046840421598972e-07    -2.046840066327600e-07
%     -3.064993734369637e-07    -3.064992597501300e-07
%     -1.538568028536247e-07    -1.538566607450800e-07
%     -3.970321582613678e-08    -3.970320694435300e-08
%     -2.569517931760856e-10    -2.569517931760900e-10
%     -2.359501483084614e-13    -2.359501483084600e-13
%     -1.614086642121038e-10    -1.614086642121000e-10
%     -1.656258463711424e-11    -1.656258463711400e-11
%     -1.189860210359939e-09    -1.189860210359900e-09
%     -1.897565438113702e-11    -1.897565438113700e-11
%     -4.671457531912893e-08    -4.671457531912900e-08
%     -1.934833449013240e-08    -1.934833449013200e-08
%     -2.352361638813250e-10    -2.352361638813200e-10
%      5.300687917042524e-06     5.300687462295200e-06
%      1.688940415078832e-06     1.688940187705200e-06
%      3.474454388197046e-06     3.474453933449700e-06
%      4.445469699021487e-07     4.445469699021500e-07
%      2.613836613818421e-06     2.613836272757900e-06
%      3.049074734917667e-07     3.049074308592000e-07
%      5.932505473538185e-07     5.932504336669800e-07
%      3.048364050073360e-07     3.048362486879300e-07
%      5.976254513484491e-08     5.976253802941801e-08
%     -4.380748794538104e-07    -4.380756709460800e-07
%      3.318385615567584e-07     3.318384779400800e-07
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0


% End of TEFUNC
%disp ('initialized')

x0  = YY;
YP0 = YP;
%
% str is always an empty matrix
%
str = [];

%
% initialize the array of sample times
%
ts  = [0 0];



% end mdlInitializeSizes

%
%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%
function sys=mdlDerivatives(t,x,u)
global YP0 YY YP
TIME = t;
if t == 0
	sys = YP0;
else
	global XMEAS XMV  
	global NN
	IDV = u(1:20,1);
    XMV = u(21:32,1);
  
  
		  
%     TEFUNC(NN,TIME,YY,YP);
%     TEFUNCVAL = TEFUNC(NN,TIME,YY,YP)
%
% Substitution of TEFUNC for call statement. MWB
%
%                Tennessee Eastman Process Control Test Problem
% 
%                     James J. Downs and Ernest F. Vogel
% 
%                   Process and Control Systems Engineering
%                        Tennessee Eastman Company
%                               P.O. Box 511
%                           Kingsport,TN  37662
% 
% ************************************************************************
%
% 						Re-Written in Matlab 5.0 
%									by	
%								Martin Braun
%
%			Department of Chemical, Bio and Materials Engineering
%						Arizona State University
%
% ************************************************************************
%   Reference:
%     "A Plant-Wide Industrial Process Control Problem"
%     Presented at the AIChE 1990 Annual Meeting
%     Industrial Challenge Problems in Process Control,Paper #24a
%     Chicago,Illinois,November 14,1990
% 
%  Revised 4-4-91 to correct error in documentation of manipulated variables
% 
%  Subroutines:
% 
%     TEFUNC - Function evaluator to be called by integrator
%     TEINIT - Initialization
%     TESUBi - Utility subroutines, i=1,2,..,8
% 
% 
%   The process simulation has 50 states (NN=50).
%   Differences between the code and its description in the paper:
% 
%   1.  Subroutine TEINIT has TIME in the argument list.  TEINIT sets TIME
%       to zero.
% 
%   2.  There are 8 utility subroutines (TESUBi) rather than 5.
% 
%   3.  Process disturbances 14 through 20 do NOT need to be used in
%       conjunction with another disturbance as stated in the paper.  All
%       disturbances can be used alone or in any combination.
%  
%   Manipulated Variables
% 
%     XMV(1,1)     D Feed Flow (stream 2)            (Corrected Order)
%     XMV(1,2)     E Feed Flow (stream 3)            (Corrected Order)
%     XMV(1,3)     A Feed Flow (stream 1)            (Corrected Order)
%     XMV(1,4)     A and C Feed Flow (stream 4)
%     XMV(1,5)     Compressor Recycle Valve
%     XMV(1,6)     Purge Valve (stream 9)
%     XMV(1,7)     Separator Pot Liquid Flow (stream 10)
%     XMV(1,8)     Stripper Liquid Product Flow (stream 11)
%     XMV(1,9)     Stripper Steam Valve
%     XMV(1,10)    Reactor Cooling Water Flow
%     XMV(1,11)    Condenser Cooling Water Flow
%     XMV(1,12)    Agitator Speed
% 
%   Continuous Process Measurements
% 
%     XMEAS(1,1)   A Feed  (stream 1)                    kscmh
%     XMEAS(1,2)   D Feed  (stream 2)                    kg/hr
%     XMEAS(1,3)   E Feed  (stream 3)                    kg/hr
%     XMEAS(1,4)   A and C Feed  (stream 4)              kscmh
%     XMEAS(1,5)   Recycle Flow  (stream 8)              kscmh
%     XMEAS(1,6)   Reactor Feed Rate  (stream 6)         kscmh
%     XMEAS(1,7)   Reactor Pressure                      kPa gauge
%     XMEAS(1,8)   Reactor Level                         %
%     XMEAS(1,9)   Reactor Temperature                   Deg C
%     XMEAS(1,10)  Purge Rate (stream 9)                 kscmh
%     XMEAS(1,11)  Product Sep Temp                      Deg C
%     XMEAS(1,12)  Product Sep Level                     %
%     XMEAS(1,13)  Prod Sep Pressure                     kPa gauge
%     XMEAS(1,14)  Prod Sep Underflow (stream 10)        m3/hr
%     XMEAS(1,15)  Stripper Level                        %
%     XMEAS(1,16)  Stripper Pressure                     kPa gauge
%     XMEAS(1,17)  Stripper Underflow (stream 11)        m3/hr
%     XMEAS(1,18)  Stripper Temperature                  Deg C
%     XMEAS(1,19)  Stripper Steam Flow                   kg/hr
%     XMEAS(1,20)  Compressor Work                       kW
%     XMEAS(1,21)  Reactor Cooling Water Outlet Temp     Deg C
%     XMEAS(1,22)  Separator Cooling Water Outlet Temp   Deg C
% 
%   Sampled Process Measurements
% 
%     Reactor Feed Analysis (Stream 6)
%         Sampling Frequency = 0.1 hr
%         Dead Time = 0.1 hr
%         Mole %
%     XMEAS(1,23)   Component A
%     XMEAS(1,24)   Component B
%     XMEAS(1,25)   Component C
%     XMEAS(1,26)   Component D
%     XMEAS(1,27)   Component E
%     XMEAS(1,28)   Component F
% 
%     Purge Gas Analysis (Stream 9)
%         Sampling Frequency = 0.1 hr
%         Dead Time = 0.1 hr
%         Mole %
%     XMEAS(1,29)   Component A
%     XMEAS(1,30)   Component B
%     XMEAS(1,31)   Component C
%     XMEAS(1,32)   Component D
%     XMEAS(1,33)   Component E
%     XMEAS(1,34)   Component F
%     XMEAS(1,35)   Component G
%     XMEAS(1,36)   Component H
% 
%     Product Analysis (Stream 11)
%         Sampling Frequency = 0.25 hr
%         Dead Time = 0.25 hr
%         Mole %
%     XMEAS(1,37)   Component D
%     XMEAS(1,38)   Component E
%     XMEAS(1,39)   Component F
%     XMEAS(1,40)   Component G
%     XMEAS(1,41)   Component H
% 
%   Process Disturbances
% 
%     IDV(1,1)   A/C Feed Ratio, B Composition Constant (Stream 4)          Step
%     IDV(1,2)   B Composition, A/C Ratio Constant (Stream 4)               Step
%     IDV(1,3)   D Feed Temperature (Stream 2)                              Step
%     IDV(1,4)   Reactor Cooling Water Inlet Temperature                    Step
%     IDV(1,5)   Condenser Cooling Water Inlet Temperature                  Step
%     IDV(1,6)   A Feed Loss (Stream 1)                                     Step
%     IDV(1,7)   C Header Pressure Loss - Reduced Availability (Stream 4)   Step
%     IDV(1,8)   A, B, C Feed Composition (Stream 4)            Random Variation
%     IDV(1,9)   D Feed Temperature (Stream 2)                  Random Variation
%     IDV(1,10)  C Feed Temperature (Stream 4)                  Random Variation
%     IDV(1,11)  Reactor Cooling Water Inlet Temperature        Random Variation
%     IDV(1,12)  Condenser Cooling Water Inlet Temperature      Random Variation
%     IDV(1,13)  Reaction Kinetics                                    Slow Drift
%     IDV(1,14)  Reactor Cooling Water Valve                            Sticking
%     IDV(1,15)  Condenser Cooling Water Valve                          Sticking
%     IDV(1,16)  Unknown
%     IDV(1,17)  Unknown
%     IDV(1,18)  Unknown
%     IDV(1,19)  Unknown
%     IDV(1,20)  Unknown
% 
%
% =============================================================================
%        Function Evaluator
% 
%          Inputs:
% 
%            NN   = Number of differential equations
%            Time = Current time(hrs)
%            YY   = Current state values
% 
%          Outputs:
% 
%            YP   = Current derivative values
% 
%   MEASUREMENT AND VALVE COMMON BLOCK

global XMEAS XMV  

%    DISTURBANCE VECTOR COMMON BLOCK

global NN
	
XMNS = 0; % Initialization by MWB
% 	NOTE: I have included isd in the /PV/ common.  This is set
% 		non-zero when the process is shutting down.
     global  UCLR UCVR UTLR UTVR...
     XLR XVR ETR ESR...
     TCR TKR DLR ...
     VLR VVR VTR ...
     PTR PPR...
     CRXR RR RH ...
     FWR TWR QUR HWR UAR...
     UCLS UCVS UTLS UTVS...
     XLS XVS ETS ESS...
     TCS TKS DLS...
     VLS VVS VTS...
     PTS PPS...
     FWS TWS QUS HWS...
     UCLC UTLC XLC...
     ETC ESC TCC DLC...
     VLC VTC QUC...
     UCVV UTVV XVV...
     ETV ESV TCV TKV...
     VTV PTV...
     VCV VRNG VTAU...
     FTM...
     FCM XST XMWS...
     HST TST SFR...
     CPFLMX CPPRMX CPDH...
     TCWR TCWS...
     HTR AGSP...
     XDEL XNS...
     TGAS TPROD VST IVST;
	 global...
     ADIST...
     BDIST...
     CDIST...
     DDIST...
     TLAST...
     TNEXT...
     HSPAN...
     HZERO...
     SSPAN...
     SZERO...
     SPSPAN...
     IDVWLK;
     global...
     AVP BVP CVP...
     AH BH CH...
     AG BG CG...
     AV...
     AD BD CD...
     XMW;
	 
for I = 1:20
	if IDV(I,1) > 0
    	IDV(I,1) = 1;
	else
		IDV(I,1) = 0;
	end
end

IDVWLK(1,1)=IDV(8,1);
IDVWLK(2,1)=IDV(8,1);
IDVWLK(3,1)=IDV(9,1);
IDVWLK(4,1)=IDV(10,1);
IDVWLK(5,1)=IDV(11,1);
IDVWLK(6,1)=IDV(12,1);
IDVWLK(7,1)=IDV(13,1);
IDVWLK(8,1)=IDV(13,1);
IDVWLK(9,1)=IDV(16,1);
IDVWLK(10,1)=IDV(17,1);
IDVWLK(11,1)=IDV(18,1);
IDVWLK(12,1)=IDV(20,1);
for I = 1:9
	if TIME >= TNEXT(I,1);
	HWLK=TNEXT(I,1)-TLAST(I,1);
    SWLK=ADIST(I,1)+HWLK*(BDIST(I,1)+HWLK*(CDIST(I,1)+HWLK*DDIST(I,1)));
    SPWLK=BDIST(I,1)+HWLK*(2*CDIST(I,1)+3*HWLK*DDIST(I,1));
    TLAST(I,1)=TNEXT(I,1);
	% TESUB5(SWLK,SPWLK,ADIST(I,1),BDIST(I,1),CDIST(I,1),DDIST(I,1),...
	% TLAST(I,1),TNEXT(I,1),HSPAN(I,1),HZERO(I,1),SSPAN(I,1),SZERO(I,1),...
	% SPSPAN(I,1),IDVWLK(I,1));
	% TESUB5(S,SP,ADIST,BDIST,CDIST,DDIST,
	% TLAST,TNEXT,HSPAN,HZERO,SSPAN,SZERO,...
	% SPSPAN,IDVFLAG) "VAL" added to subroutine variables.  "I" changed to "It".
	% Substitution of TESUB5 for call statement. MWB
	S5 = SWLK;
	SP5 = SPWLK;
	ADISTVAL = ADIST(I,1);
	BDISTVAL = BDIST(I,1);
	CDISTVAL = CDIST(I,1);
	DDISTVAL = DDIST(I,1);
	TLASTVAL = TLAST(I,1);
	TNEXTVAL = TNEXT(I,1);
	HSPANVAL = HSPAN(I,1);
	HZEROVAL = HZERO(I,1);
	SSPANVAL = SSPAN(I,1);
	SZEROVAL = SZERO(I,1);
	SPSPANVAL = SPSPAN(I,1);
	IDVFLAGVAL = IDVWLK(I,1);
	It = -1;
    H5 = HSPANVAL * TESUB7(It) + HZEROVAL;
    S1 = SSPANVAL * TESUB7(It) * IDVFLAGVAL + SZEROVAL;
    S1P = SPSPANVAL * TESUB7(It) * IDVFLAGVAL;
    ADISTVAL = S5;
    BDISTVAL = SP5;
    CDISTVAL = (3 * (S1 - S5) - H5 * (S1P + 2 * SP5)) / H5^2;
    DDISTVAL = (2 * (S5 - S1) + H5 * (S1P + SP5)) / H5^3;
    TNEXTVAL = TLASTVAL + H5;
	ADIST(I,1) = ADISTVAL; 
	BDIST(I,1) = BDISTVAL;
	CDIST(I,1) = CDISTVAL;
	DDIST(I,1) = DDISTVAL;
	TNEXT(I,1) = TNEXTVAL;  
	%
	% End of TESUB5
	%
	end
end
for I = 10:12
    if TIME >= TNEXT(I,1)
    HWLK = TNEXT(I,1) - TLAST(I,1);
    SWLK = ADIST(I,1) + HWLK * (BDIST(I,1) + HWLK * (CDIST(I,1) + HWLK * DDIST(I,1)));
    SPWLK = BDIST(I,1) + HWLK * (2 * CDIST(I,1) + 3 * HWLK * DDIST(I,1));
    TLAST(I,1) = TNEXT(I,1);
    	if SWLK > 0.1
      	ADIST(I,1) = SWLK;
      	BDIST(I,1) = SPWLK;
      	CDIST(I,1) = -(3 * SWLK + 0.2 * SPWLK) / 0.01;
      	DDIST(I,1) = (2 * SWLK + 0.1 * SPWLK) / 0.001;
      	TNEXT(I,1) = TLAST(I,1) + 0.1;
      	else
      	ISD = -1;
      	HWLK = HSPAN(I,1) * TESUB7(ISD) + HZERO(I,1);
     	ADIST(I,1) = 0;
      	BDIST(I,1) = 0;
      	CDIST(I,1) = (IDVWLK(I,1)) / HWLK^2;
      	DDIST(I,1) = 0;
      	TNEXT(I,1) = TLAST(I,1) + HWLK;
  		end
	end
end
if TIME == 0
	for I = 1:12
    ADIST(I,1)=SZERO(I,1);
    BDIST(I,1)=0;
    CDIST(I,1)=0;
    DDIST(I,1)=0;
    TLAST(I,1)=0.0;
    TNEXT(I,1)=0.1;
  	end
end

XST(1,4) = TESUB8(1,TIME) - IDV(1,1) * 0.03 - IDV(2,1) * 2.43719 * 10 ^ (-3);
XST(2,4) = TESUB8(2,TIME) + IDV(2,1) * 0.005;
XST(3,4) = 1 - XST(1,4) - XST(2,4);
TST(1,1) = TESUB8(3,TIME) + IDV(3,1) * 5;
TST(4,1) = TESUB8(4,TIME);
TCWR = TESUB8(5,TIME) + IDV(4,1) * 5;
TCWS = TESUB8(6,TIME) + IDV(5,1) * 5;
R1F = TESUB8(7,TIME);
R2F = TESUB8(8,TIME);
for I = 1:3
      UCVR(I,1)=YY(I,1);
      UCVS(I,1)=YY(I+9,1);
      UCLR(I,1)=0.0;
      UCLS(I,1)=0.0;
end
for I = 4:8
      UCLR(I,1) = YY(I,1);
      UCLS(I,1) = YY(I+9,1);
end
  
for I = 1:8
      UCLC(I,1) = YY(I+18,1);
      UCVV(I,1) = YY(I+27,1);
end
  
      ETR = YY(9,1);
      ETS = YY(18,1);
      ETC = YY(27,1);
      ETV = YY(36,1);
      TWR = YY(37,1);
      TWS = YY(38,1);
	  
%	  TWS shows an error here.  For the initialization,
%     it's value is  7.729698353000001e+01, when it should be
%     77.296983530000
	  
for I = 1:12
      VPOS(I,1) = YY(I+38,1);
end
      UTLR = 0.0;
      UTLS = 0.0;
      UTLC = 0.0;
      UTVV = 0.0;
	  
for I = 1:8
      UTLR = UTLR + UCLR(I,1);
      UTLS = UTLS + UCLS(I,1);
      UTLC = UTLC + UCLC(I,1);
      UTVV = UTVV + UCVV(I,1);
end
for I = 1:8
      XLR(I,1) = UCLR(I,1)/UTLR;
      XLS(I,1) = UCLS(I,1)/UTLS;
      XLC(I,1) = UCLC(I,1)/UTLC;
      XVV(I,1) = UCVV(I,1)/UTVV;
end
  
      ESR = ETR/UTLR;
      ESS = ETS/UTLS;
      ESC = ETC/UTLC;
      ESV = ETV/UTVV;
% Additional precision errors show up in XLR XLS XLC XVV ESR ESS ESC ESV 

%     TESUB2(XLR,TCR,ESR,0);
%	  TESUB2(Z,T,H,ITY)
% Substitution of TESUB2 for call statement. MWB
% TESUB2 appears to update T and only T.
Z = XLR;
T = TCR;
H = ESR;
ITY = 0;
	  TIN=T;
for J = 1:100
	%      TESUB1(Z,T,HTEST,ITY);
	%	   TESUB1VAL = TESUB1(Z,T,H,ITY)	
	% Substitution of TESUB1 for call statement MWB
	%
	if ITY == 0
    	HTEST = 0;
		for I = 1:8
      		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
      		HI = 1.8 * HI;
      		HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	else
      	HTEST = 0;
		for I = 1:8 
      		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
      		HI = 1.8 * HI;
      		HI = HI + AV(I,1);
      		HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	end
	if ITY == 2
      	R1 = 3.57696 / 1000000;
     	HTEST = HTEST - R1 * (T + 273.15);
  	end
	%
	% End of TESUB1 sub
	%
	
    ERR = HTEST - H;
	
	% 		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	%		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	% Substitution of TESUB3 for call statement MWB
	%
	if ITY == 0
      	DH = 0.0;
    	for I = 1:8
      		DHI = AH(I,1) + BH(I,1) * T + CH(I,1) * T^2;
      		DHI = 1.8 * DHI;
      		DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	else
      	DH=0.0;
    	for I = 1:8
      		DHI = AG(I,1) + BG(I,1) * T + CG(I,1) * T^2;
      		DHI = 1.8 * DHI;
      		DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	end
	if ITY == 2
      	R3 = 3.57696/1000000;
     	DH = DH - R3;
 	end
	%
	% End of TESUB3 substitution.
	%
    DT = -ERR / DH;
    T = T + DT;
  	if abs(DT) < 1*10^(-12)
    	break
  	end
  	if J == 100
  		T = TIN;
  	end
end
TCR = T;
% End of TESUB2.

      TKR = TCR+273.15;

%     TESUB2(XLS,TCS,ESS,0);
%	  TESUB2(Z,T,H,ITY)
% Substitution of TESUB2 for call statement. MWB
% TESUB2 appears to update T and only T.
Z = XLS;
T = TCS;
H = ESS;
ITY = 0;
TIN=T;
for J = 1:100
	%      TESUB1(Z,T,HTEST,ITY);
	%	   TESUB1VAL = TESUB1(Z,T,H,ITY)	
	% Substitution of TESUB1 for call statement MWB
	%	
	if ITY == 0
		HTEST = 0;
		for I = 1:8
			HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	else
		HTEST = 0;
		for I = 1:8 
			HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HI = HI + AV(I,1);
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	end
	if ITY == 2
		R1 = 3.57696 / 1000000;
		HTEST = HTEST - R1 * (T + 273.15);
	end
	%
	% End of TESUB1 sub
	%
 	ERR = HTEST - H;
	% 		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	%		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	% Substitution of TESUB3 for call statement MWB
	%
	if ITY == 0
   		DH = 0.0;
   		for I = 1:8
   			DHI = AH(I,1) + BH(I,1) * T + CH(I,1) * T^2;
   			DHI = 1.8 * DHI;
   			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	else
   		DH=0.0;
   		for I = 1:8
   			DHI = AG(I,1) + BG(I,1) * T + CG(I,1) * T^2;
   			DHI = 1.8 * DHI;
   			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	end
	if ITY == 2
    	R3 = 3.57696/1000000;
    	DH = DH - R3;
  	end
	%
	% End of TESUB3 substitution.
	%
    DT = -ERR / DH;
    T = T + DT;
	if abs(DT) < 1*10^(-12)
    	break
  	end
	if J == 100
  		T = TIN;
  	end
end
TCS = T;
% End of TESUB2.

      TKS = TCS+273.15;
      
%     TESUB2(XLC,TCC,ESC,0);
%	  TESUB2(Z,T,H,ITY)
% Substitution of TESUB2 for call statement. MWB
% TESUB2 appears to update T and only T.
Z = XLC;
T = TCC;
H = ESC;
ITY = 0;
TIN=T;
for J = 1:100
	%      TESUB1(Z,T,HTEST,ITY);
	%	   TESUB1VAL = TESUB1(Z,T,H,ITY)	
	% Substitution of TESUB1 for call statement MWB
	%
	if ITY == 0
		HTEST = 0;
		for I = 1:8
			HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	else
		HTEST = 0;
		for I = 1:8 
			HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HI = HI + AV(I,1);
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	end
	if ITY == 2
		R1 = 3.57696 / 1000000;
		HTEST = HTEST - R1 * (T + 273.15);
	end
	%
	% End of TESUB1 sub
	%
	ERR = HTEST - H;
	% 		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	%		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	% Substitution of TESUB3 for call statement MWB
	%
	if ITY == 0
		DH = 0.0;
    	for I = 1:8
			DHI = AH(I,1) + BH(I,1) * T + CH(I,1) * T^2;
			DHI = 1.8 * DHI;
			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	else
		DH=0.0;
		for I = 1:8
			DHI = AG(I,1) + BG(I,1) * T + CG(I,1) * T^2;
			DHI = 1.8 * DHI;
			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	end
	if ITY == 2
		R3 = 3.57696/1000000;
		DH = DH - R3;
	end
	%
	% End of TESUB3 substitution.
	%
	DT = -ERR / DH;
	T = T + DT;
	if abs(DT) < 1*10^(-12)
		break
	end
	if J == 100
		T = TIN;
	end
end
TCC = T;
% End of TESUB2.
      
%     TESUB2(XVV,TCV,ESV,2);
%	  TESUB2(Z,T,H,ITY)
% Substitution of TESUB2 for call statement. MWB
% TESUB2 appears to update T and only T.
Z = XVV;
T = TCV;
H = ESV;
ITY = 2;
TIN=T;
for J = 1:100
	%      TESUB1(Z,T,HTEST,ITY);
	%	   TESUB1VAL = TESUB1(Z,T,H,ITY)	
	% Substitution of TESUB1 for call statement MWB
	%
	if ITY == 0
		HTEST = 0;
		for I = 1:8
			HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	else
		HTEST = 0;
		for I = 1:8 
			HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HI = HI + AV(I,1);
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	end
	if ITY == 2
		R1 = 3.57696 / 1000000;
		HTEST = HTEST - R1 * (T + 273.15);
	end
	%
	% End of TESUB1 sub
	%
	ERR = HTEST - H;
	% 		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	%		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	% Substitution of TESUB3 for call statement MWB
	%
	if ITY == 0
		DH = 0.0;
    	for I = 1:8
			DHI = AH(I,1) + BH(I,1) * T + CH(I,1) * T^2;
			DHI = 1.8 * DHI;
			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	else
		DH=0.0;
    	for I = 1:8
			DHI = AG(I,1) + BG(I,1) * T + CG(I,1) * T^2;
			DHI = 1.8 * DHI;
			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	end
	if ITY == 2
		R3 = 3.57696/1000000;
		DH = DH - R3;
	end
	%
	% End of TESUB3 substitution.
	%
	DT = -ERR / DH;
	T = T + DT;
	if abs(DT) < 1*10^(-12)
		break
	end
	if J == 100
		T = TIN;
	end
end
TCV = T;
% End of TESUB2.
	  
TKV = TCV+273.15;
      
%	  TESUB4(XLR,TCR,DLR);
%	  TESUB4(X,T,R)
% Substitution of TESUB4 for call statement. MWB
X = XLR;
T4 = TCR;
R4 = DLR;
V = 0.0;
for I = 1:8
	V = V + X(I,1) * XMW(I,1) / (AD(I,1) + (BD(I,1) + CD(I,1) * T4) * T4);
end
R4 = 1.0 / V;
DLR = R4;
% End of TESUB4.

      
%	  TESUB4(XLS,TCS,DLS);
%	  TESUB4(X,T,R)
% Substitution of TESUB4 for call statement. MWB
X = XLS;
T4 = TCS;
R4 = DLS;
V = 0.0;
for I = 1:8
	V = V + X(I,1) * XMW(I,1) / (AD(I,1) + (BD(I,1) + CD(I,1) * T4) * T4);
end
R4 = 1.0 / V;
DLS = R4;
% End of TESUB4.
      
%	  TESUB4(XLC,TCC,DLC);
%	  TESUB4(X,T,R)
% Substitution of TESUB4 for call statement. MWB
X = XLC;
T4 = TCC;
R4 = DLC;
V = 0.0;
for I = 1:8
      V = V + X(I,1) * XMW(I,1) / (AD(I,1) + (BD(I,1) + CD(I,1) * T4) * T4);
end
R4 = 1.0 / V;
DLC = R4;
% End of TESUB4.
VLR = UTLR/DLR;
VLS = UTLS/DLS;
VLC = UTLC/DLC;
VVR = VTR-VLR;
VVS = VTS-VLS;
RG = 998.9;
PTR = 0.0;
PTS = 0.0;
for I = 1:3
      PPR(I,1) = UCVR(I,1)*RG*TKR/VVR;
      PTR = PTR+PPR(I,1);
      PPS(I,1) = UCVS(I,1)*RG*TKS/VVS;
      PTS = PTS+PPS(I,1);
  end
% Some precision error with PTS
for I = 4:8
      VPR=exp(AVP(I,1)+BVP(I,1)/(TCR+CVP(I,1)));
      PPR(I,1)=VPR*XLR(I,1);
      PTR=PTR+PPR(I,1);
      VPR=exp(AVP(I,1)+BVP(I,1)/(TCS+CVP(I,1)));
      PPS(I,1)=VPR*XLS(I,1);
      PTS=PTS+PPS(I,1);
  end
      PTV=UTVV*RG*TKV/VTV;
for I=1:8
      XVR(I,1)=PPR(I,1)/PTR;
      XVS(I,1)=PPS(I,1)/PTS;
end
UTVR=PTR*VVR/RG/TKR;
UTVS=PTS*VVS/RG/TKS;
for I=4:8
      UCVR(I,1)=UTVR*XVR(I,1);
      UCVS(I,1)=UTVS*XVS(I,1);
  end
  
% This is where the Matlab RR's show different values than the FORTRAN 77
% RR's.  The code below has been modified with the chop command to 
% elliminate the stray numbers that are created by the matlab code.
% Matlab:				FORTRAN 77:
% R1F
% 1						1.0000000000000
% R2F
% 1						1.0000000000000
% TKR
% 3.935500000017000e+02	393.55000000169
% RR(1,1)
% 3.181229946057000e-09	3.1812215872809D-09
% RR(2,1)
% 1.569591352140000e-10	1.5695892900646D-10
% RR(3,1)
% 7.437635682791000e-11	7.4376432806282D-11
% RR(4,1)
% 5.708298619084000e-11	5.7083044503356D-11
%
% NOTE:  I DID A COMPARISON OF THE EXP FUNCTIONS OF MATLAB
%        AND FORTRAN 77 AND THERE IS A CONSISTENT DISCREPANCY AT
%        10^-16 BETWEEN THE OUTPUT OF THE TWO LANGUAGES
%

%R1F = chop(R1F,13);
%R2F = chop(R2F,13);
%TKR = chop(TKR,13);
%RR(1,1)=chop(exp(31.5859536-40000.0/1.987/TKR)*R1F,13);
%RR(2,1)=chop(exp(3.00094014-20000.0/1.987/TKR)*R2F,13);
%RR(3,1)=chop(exp(53.4060443-60000.0/1.987/TKR),13);
%RR(4,1)=chop(RR(3,1)*0.767488334,13);

% This next section has the following output:
% Matlab:				FORTRAN 77:
% PPR(1,1)
% 5.717769731191162e+03	5717.7697311912     * good
% PPR(3,1)
% 4.159943934675613e+03 4159.9439346756		* good
% R1F
% 2.174415112221978e+04 21744.151122220		* good
% R2F
% 2.247652526817436e+01 22.476525268174		* good
% RR(1,1)
% 3.517133031511259e-01 0.35171237901398	
% RR(2,1)
% 2.861611337089543e-01 0.28616075776025

if PPR(1,1) > 0.0 & PPR(3,1) > 0.0
      R1F=PPR(1,1)^1.1544;
      R2F=PPR(3,1)^0.3735;
      RR(1,1)=RR(1,1)*R1F*R2F*PPR(4,1);
      RR(2,1)=RR(2,1)*R1F*R2F*PPR(5,1);
  else
	  
      RR(1,1)=0.0;
      RR(2,1)=0.0;
  end
%
% This next section has the following output:
% Matlab:				FORTRAN 77:
% RR(1,1)
% 2.516069042865577e+02 251.60624318195
% RR(2,1)
% 2.047125210635105e+02 204.71225211921
% RR(3,1)
% 1.134874248863408e+00 1.1348754081821
% RR(4,1)
% 5.281888235868223e-02 5.2818936315249D-02
%
RR(3,1)=RR(3,1)*PPR(1,1)*PPR(5,1);
RR(4,1)=RR(4,1)*PPR(1,1)*PPR(4,1);
for I=1:4
    RR(I,1)=RR(I,1)*VVR;
end
%
% This next section has the following output:
% Matlab:					FORTRAN 77:
% CRXR(1,1)
% -4.574542995989316e+02	-457.45337070934
% CRXR(3,1)
% -4.563194253500682e+02	-456.31849530116
% CRXR(4,1)
% -2.516861326100957e+02	-251.68547158642
% CRXR(5,1)
% -2.058473953123739e+02	-205.84712752739
% CRXR(6,1)
% 1.187693131222090e+00		1.1876943444973
% CRXR(7,1)
% 2.516069042865577e+02		251.60624318195
% CRXR(8,1)
% 2.047125210635105e+02		204.71225211921
% RH
% 2.759494513807820e+01		27.594886078737
CRXR(1,1)=-RR(1,1)-RR(2,1)-RR(3,1);
CRXR(3,1)=-RR(1,1)-RR(2,1);
CRXR(4,1)=-RR(1,1)-1.5*RR(4,1);
CRXR(5,1)=-RR(2,1)-RR(3,1);
CRXR(6,1)=RR(3,1)+RR(4,1);
CRXR(7,1)=RR(1,1);
CRXR(8,1)=RR(2,1);
RH=RR(1,1)*HTR(1,1)+RR(2,1)*HTR(2,1);
XMWS(1,1)=0.0;
XMWS(2,1)=0.0;
XMWS(6,1)=0.0;
XMWS(8,1)=0.0;
XMWS(9,1)=0.0;
XMWS(10,1)=0.0;
%
% XST & XMWS seem to be calculated correcly in the next 
% section with little error.
%
for I=1:8
	  XST(I,6)=XVV(I,1);
	  XST(I,8)=XVR(I,1);
	  XST(I,9)=XVS(I,1);
      XST(I,10)=XVS(I,1);
      XST(I,11)=XLS(I,1);
      XST(I,13)=XLC(I,1);
      XMWS(1,1)=XMWS(1,1)+XST(I,1)*XMW(I,1);
      XMWS(2,1)=XMWS(2,1)+XST(I,2)*XMW(I,1);
      XMWS(6,1)=XMWS(6,1)+XST(I,6)*XMW(I,1);
      XMWS(8,1)=XMWS(8,1)+XST(I,8)*XMW(I,1);
      XMWS(9,1)=XMWS(9,1)+XST(I,9)*XMW(I,1);
      XMWS(10,1)=XMWS(10,1)+XST(I,10)*XMW(I,1);
  end
TST(6,1)=TCV;
TST(8,1)=TCR;
TST(9,1)=TCS;
TST(10,1)=TCS;
TST(11,1)=TCS;
TST(13,1)=TCC;

% The HST's from the following substitutions for TESUB1 are 
% computed correctly.

%
%      TESUB1(XST(1,1),TST(1,1),HST(1,1),1);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,1);
T=TST(1,1);
H=HST(1,1);
ITY=1;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
	R1 = 3.57696 / 1000000;
	H = H - R1 * (T + 273.15);
end
HST(1,1) = H;
% end of TESUB1


%      TESUB1(XST(1,2),TST(2,1),HST(2,1),1);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,2);
T=TST(2,1);
H=HST(2,1);
ITY=1;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(2,1) = H;
% end of TESUB1

%      TESUB1(XST(1,3),TST(3,1),HST(3,1),1);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,3);
T=TST(3,1);
H=HST(3,1);
ITY=1;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(3,1) = H;
% end of TESUB1

%      TESUB1(XST(1,4),TST(4,1),HST(4,1),1);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,4);
T=TST(4,1);
H=HST(4,1);
ITY=1;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(4,1) = H;
% end of TESUB1

%      TESUB1(XST(1,6),TST(6,1),HST(6,1),1);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,6);
T=TST(6,1);
H=HST(6,1);
ITY=1;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(6,1) = H;
% end of TESUB1

%      TESUB1(XST(1,8),TST(8,1),HST(8,1),1);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,8);
T=TST(8,1);
H=HST(8,1);
ITY=1;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(8,1) = H;
% end of TESUB1

%      TESUB1(XST(1,9),TST(9,1),HST(9,1),1);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,9);
T=TST(9,1);
H=HST(9,1);
ITY=1;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(9,1) = H;
% end of TESUB1

      HST(10,1)=HST(9,1);
	  
%      TESUB1(XST(1,11),TST(11,1),HST(11,1),0);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,11);
T=TST(11,1);
H=HST(11,1);
ITY=0;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
 		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(11,1) = H;
% end of TESUB1

%      TESUB1(XST(1,13),TST(13,1),HST(13,1),0);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,13);
T=TST(13,1);
H=HST(13,1);
ITY=0;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(13,1) = H;
% end of TESUB1

% This next section is good.
FTM(1,1)=VPOS(1,1)*VRNG(1,1)/100.0;
FTM(2,1)=VPOS(2,1)*VRNG(2,1)/100.0;
FTM(3,1)=VPOS(3,1)*(1-IDV(6,1))*VRNG(3,1)/100.0;
FTM(4,1)=VPOS(4,1)*(1-IDV(7,1)*0.2)*VRNG(4,1)/100.0+1.E-10;
FTM(11,1)=VPOS(7,1)*VRNG(7,1)/100.0;
FTM(13,1)=VPOS(8,1)*VRNG(8,1)/100.0;
UAC=VPOS(9,1)*VRNG(9,1)*(1+TESUB8(9,TIME))/100.0;
FWR=VPOS(10,1)*VRNG(10,1)/100.0;
FWS=VPOS(11,1)*VRNG(11,1)/100.0;
AGSP=(VPOS(12,1)+150.0)/100.0;
DLP=PTV-PTR;

if DLP<0 
	DLP=0.0;
end
FLMS=1937.6*sqrt(DLP);
FTM(6,1)=FLMS/XMWS(6,1);
DLP=PTR-PTS;
if DLP<0 
	DLP=0.0;
end
FLMS=4574.21*sqrt(DLP)*(1-0.25*TESUB8(12,TIME));
FTM(8,1)=FLMS/XMWS(8,1);
DLP=PTS-760.0;
if DLP < 0 
	DLP=0.0;
end
FLMS=VPOS(6,1)*0.151169*sqrt(DLP);
FTM(10,1)=FLMS/XMWS(10,1);
PR=PTV/PTS;
if PR < 1
	PR=1.0;
end
if PR > CPPRMX
	PR=CPPRMX;
end
FLCOEF=CPFLMX/1.197;
FLMS=CPFLMX+FLCOEF*(1.0-PR^3);
CPDH=FLMS*(TCS+273.15)*1.8E-6*1.9872*(PTV-PTS)/(XMWS(9,1)*PTS);
DLP=PTV-PTS;
if DLP < 0.0
	DLP=0.0;
end
FLMS=FLMS-VPOS(5,1)*53.349*sqrt(DLP);
if FLMS < 1.E-3
	FLMS=1.E-3;
end
FTM(9,1)=FLMS/XMWS(9,1);
HST(9,1)=HST(9,1)+CPDH/FTM(9,1);
for I=1:8
      FCM(I,1)=XST(I,1)*FTM(1,1);
      FCM(I,2)=XST(I,2)*FTM(2,1);
      FCM(I,3)=XST(I,3)*FTM(3,1);
      FCM(I,4)=XST(I,4)*FTM(4,1);
      FCM(I,6)=XST(I,6)*FTM(6,1);
      FCM(I,8)=XST(I,8)*FTM(8,1);
      FCM(I,9)=XST(I,9)*FTM(9,1);
      FCM(I,10)=XST(I,10)*FTM(10,1);
      FCM(I,11)=XST(I,11)*FTM(11,1);
      FCM(I,13)=XST(I,13)*FTM(13,1);
  end

if FTM(11,1) > 0.1
	if TCC > 170
		TMPFAC=TCC-120.262;
	elseif TCC < 5.292
		TMPFAC=0.1;
	else
		TMPFAC=363.744/(177.-TCC)-2.22579488;
	end
	VOVRL=FTM(4,1)/FTM(11,1)*TMPFAC;
	SFR(4,1)=8.5010*VOVRL/(1.0+8.5010*VOVRL);
	SFR(5,1)=11.402*VOVRL/(1.0+11.402*VOVRL);
	SFR(6,1)=11.795*VOVRL/(1.0+11.795*VOVRL);
	SFR(7,1)=0.0480*VOVRL/(1.0+0.0480*VOVRL);
	SFR(8,1)=0.0242*VOVRL/(1.0+0.0242*VOVRL);
else
	SFR(4,1)=0.9999;
	SFR(5,1)=0.999;
	SFR(6,1)=0.999;
	SFR(7,1)=0.99;
	SFR(8,1)=0.98;
end
for I=1:8
      FIN(I,1)=0.0;
      FIN(I,1)=FIN(I,1)+FCM(I,4);
      FIN(I,1)=FIN(I,1)+FCM(I,11);
  end

FTM(5,1)=0.0;
FTM(12,1)=0.0;
for I=1:8
      FCM(I,5)=SFR(I,1)*FIN(I,1);
      FCM(I,12)=FIN(I,1)-FCM(I,5);
      FTM(5,1)=FTM(5,1)+FCM(I,5);
      FTM(12,1)=FTM(12,1)+FCM(I,12);
end
for I=1:8
      XST(I,5)=FCM(I,5)/FTM(5,1);
      XST(I,12)=FCM(I,12)/FTM(12,1);
end
TST(5,1)=TCC;
TST(12,1)=TCC;
%      TESUB1(XST(1,5),TST(5,1),HST(5,1),1);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,5);
T=TST(5,1);
H=HST(5,1);
ITY=1;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
      R1 = 3.57696 / 1000000;
      H = H - R1 * (T + 273.15);
end
HST(5,1) = H;
% end of TESUB1

%      TESUB1(XST(1,12),TST(12,1),HST(12,1),0);
%	   TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,12);
T=TST(12,1);
H=HST(12,1);
ITY=0;
if ITY == 0
	H = 0;
	for I = 1:8
		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
else
	H = 0;
	for I = 1:8 
		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
		HI = 1.8 * HI;
		HI = HI + AV(I,1);
		H = H + Z(I,1) * XMW(I,1) * HI;
	end
end
if ITY == 2
	R1 = 3.57696 / 1000000;
	H = H - R1 * (T + 273.15);
end
HST(12,1) = H;
% end of TESUB1

FTM(7,1)=FTM(6,1);
HST(7,1)=HST(6,1);
TST(7,1)=TST(6,1);

for I=1:8
	
      XST(I,7)=XST(I,6);
      FCM(I,7)=FCM(I,6);
end

if VLR/7.8 > 50.0
	UARLEV=1.0;
elseif VLR/7.8 < 10.0
    UARLEV=0.0;
else
      UARLEV=0.025*VLR/7.8-0.25;
  end
  
UAR=UARLEV*(-0.5*AGSP^2+2.75*AGSP-2.5)*855490.E-6;
QUR=UAR*(TWR-TCR)*(1-0.35*TESUB8(10,TIME));
UAS=0.404655*(1.0-1.0/(1.0+(FTM(8,1)/3528.73)^4));
QUS=UAS*(TWS-TST(8,1))*(1-0.25*TESUB8(11,TIME));
QUC=0;

if TCC < 100
	QUC=UAC*(100.0-TCC);
end
      XMEAS(1,1)=FTM(3,1)*0.359/35.3145;
      XMEAS(2,1)=FTM(1,1)*XMWS(1,1)*0.454;
      XMEAS(3,1)=FTM(2,1)*XMWS(2,1)*0.454;
      XMEAS(4,1)=FTM(4,1)*0.359/35.3145;
      XMEAS(5,1)=FTM(9,1)*0.359/35.3145;
      XMEAS(6,1)=FTM(6,1)*0.359/35.3145;
      XMEAS(7,1)=(PTR-760.0)/760.0*101.325;
      XMEAS(8,1)=(VLR-84.6)/666.7*100.0;
      XMEAS(9,1)=TCR;
      XMEAS(10,1)=FTM(10,1)*0.359/35.3145;
      XMEAS(11,1)=TCS;
      XMEAS(12,1)=(VLS-27.5)/290.0*100.0;
      XMEAS(13,1)=(PTS-760.0)/760.0*101.325;
      XMEAS(14,1)=FTM(11,1)/DLS/35.3145;
      XMEAS(15,1)=(VLC-78.25)/VTC*100.0;
      XMEAS(16,1)=(PTV-760.0)/760.0*101.325;
      XMEAS(17,1)=FTM(13,1)/DLC/35.3145;
      XMEAS(18,1)=TCC;
      XMEAS(19,1)=QUC*1.04E3*0.454;
      XMEAS(20,1)=CPDH*0.0003927E6;
      XMEAS(20,1)=CPDH*0.29307E3;
      XMEAS(21,1)=TWR;
      XMEAS(22,1)=TWS;
	 ISD=0;
if XMEAS(7,1) > 3000.0
	ISD=1;
	end
if VLR/35.3145 > 24.0
	ISD=1;
	end
if VLR/35.3145 < 2.0
	ISD=1;
	end
if XMEAS(9,1) > 175.0
	ISD=1;
	end
if VLS/35.3145 > 12.0
	ISD=1;
	end
if VLS/35.3145 < 1.0
	ISD=1;
	end
if VLC/35.3145 > 8.0
	ISD=1;
	end
if VLC/35.3145 < 1.0
	ISD=1;
	end
if TIME > 0.0 & ISD == 0
	for I=1:22
		%      TESUB6(XNS(I,1),XMNS)
		%      TESUB6(STD,X)
		% Substitution of TESUB6 for call statement. MWB
		STD = XNS(I,1);
		X = XMNS;
		X = 0;
		for Iter =1:12
			X = X + TESUB7(Iter);
		end
		X = (X - 6) * STD;
		XMNS = X;
		% End of TESUB6.

		XMEAS(I,1)=XMEAS(I,1)+XMNS;
	end
end
      XCMP(23,1)=XST(1,7)*100.0;
      XCMP(24,1)=XST(2,7)*100.0;
      XCMP(25,1)=XST(3,7)*100.0;
      XCMP(26,1)=XST(4,7)*100.0;
      XCMP(27,1)=XST(5,7)*100.0;
      XCMP(28,1)=XST(6,7)*100.0;
      XCMP(29,1)=XST(1,10)*100.0;
      XCMP(30,1)=XST(2,10)*100.0;
      XCMP(31,1)=XST(3,10)*100.0;
      XCMP(32,1)=XST(4,10)*100.0;
      XCMP(33,1)=XST(5,10)*100.0;
      XCMP(34,1)=XST(6,10)*100.0;
      XCMP(35,1)=XST(7,10)*100.0;
      XCMP(36,1)=XST(8,10)*100.0;
      XCMP(37,1)=XST(4,13)*100.0;
      XCMP(38,1)=XST(5,13)*100.0;
      XCMP(39,1)=XST(6,13)*100.0;
      XCMP(40,1)=XST(7,13)*100.0;
      XCMP(41,1)=XST(8,13)*100.0;
	  
if TIME == 0
	for I=23:41
		XDEL(I,1)=XCMP(I,1);
		XMEAS(I,1)=XCMP(I,1);
	end

	TGAS=0.1;
	TPROD=0.25;
end
if TIME >= TGAS
	for I=23:36
		XMEAS(I,1)=XDEL(I,1);
		%      TESUB6(XNS(I,1),XMNS);
		%      TESUB6(STD,X)
		% Substitution of TESUB6 for call statement. MWB
		STD = XNS(I,1);
		X = XMNS;
		X = 0;
		for Iter=1:12
			X = X + TESUB7(Iter);
		end
		X = (X - 6) * STD;
		XMNS = X;
		% End of TESUB6

		XMEAS(I,1)=XMEAS(I,1)+XMNS;
		XDEL(I,1)=XCMP(I,1);
	end
	TGAS=TGAS+0.1;
end
 
if TIME >= TPROD
	for I=37:41
		XMEAS(I,1)=XDEL(I,1);
      
		%      TESUB6(XNS(I,1),XMNS)
		%      TESUB6(STD,X)
		% Substitution of TESUB6 for call statement. MWB
		STD = XNS(I,1);
		X = XMNS;
		X = 0;
		for Iter=1:12
			X = X + TESUB7(Iter);
		end
		X = (X - 6) * STD;
		XMNS = X;
		% End of TESUB6

		XMEAS(I,1)=XMEAS(I,1)+XMNS;
		XDEL(I,1)=XCMP(I,1);
	end
	TPROD=TPROD+0.25;
end
% 
for iter = 1:8
      YP(iter,1) = FCM(iter,7) - FCM(iter,8) + CRXR(iter,1);
      YP(iter+9,1) = FCM(iter,8) - FCM(iter,9) - FCM(iter,10) - FCM(iter,11);
      YP(iter+18,1) = FCM(iter,12) - FCM(iter,13);
      YP(iter+27,1) = FCM(iter,1) + FCM(iter,2) + FCM(iter,3) + FCM(iter,5) + FCM(iter,9) - FCM(iter,6);
  end
  
% 		Here is the "correct" version of the separator energy balance:
% 
% 	YP(18)=HST(8)*FTM(8)-(HST(9)*FTM(9)-cpdh)-HST(10)*FTM(10)-HST(11)*FTM(11)+QUS
% 
%   	Here is the original version
	  YP(9,1)=HST(7,1) * FTM(7,1) - HST(8) * FTM(8) + RH + QUR;
      YP(18,1)=HST(8,1)*FTM(8,1)-HST(9,1)*FTM(9,1)-HST(10,1)*FTM(10,1)-HST(11,1)*FTM(11,1)+QUS;
      YP(27,1)=HST(4,1)*FTM(4,1)+HST(11,1)*FTM(11,1)-HST(5,1)*FTM(5,1)-HST(13,1)*FTM(13,1)+QUC;
      YP(36,1)=HST(1,1)*FTM(1,1)+HST(2,1)*FTM(2,1)+HST(3,1)*FTM(3,1)+HST(5,1)*FTM(5,1)+HST(9,1)*FTM(9,1)-HST(6,1)*FTM(6,1);
      YP(37,1)=(FWR*500.53*(TCWR-TWR)-QUR*1000000/1.8)/HWR;
      YP(38,1)=(FWS*500.53*(TCWS-TWS)-QUS*1000000/1.8)/HWS;
      IVST(10,1)=IDV(14,1);
      IVST(11,1)=IDV(15,1);
      IVST(5,1)=IDV(19,1);
      IVST(7,1)=IDV(19,1);
      IVST(8,1)=IDV(19,1);
      IVST(9,1)=IDV(19,1);
for I=1:12
	if TIME == 0 | abs(VCV(I,1)-XMV(I,1)) > VST(I,1)*IVST(I,1)
		VCV(I,1)=XMV(I,1);
	end
	if VCV(I,1) < 0.0
		VCV(I,1)=0.0;
	end
	if VCV(I,1) > 100.0
		VCV(I,1)=100.0;
	end
	YP(I+38,1)=(VCV(I,1)-VPOS(I,1))/VTAU(I,1);
end
if TIME > 0.0 & ISD ~= 0
	for I=1:NN
		YP(I,1)=0.0;
	end
end
% This is the final output of the subroutine:
% Everything up to this point (except what has been
% documented in the code) runs very similar to the FORTRAN 77
% version. MWB 
%  7/12/98
%  Matlab:						FORTRAN77:
%  YY
% 	   1.040491389000000e+01     1.040491389000000e+01
%      4.363996017000000e+00     4.363996017000000e+00
%      7.570059737000000e+00     7.570059737000000e+00
%      4.230042431000000e-01     4.230042431000000e-01
%      2.415513437000000e+01     2.415513437000000e+01
%      2.942597645000000e+00     2.942597645000000e+00
%      1.543770655000000e+02     1.543770655000000e+02
%      1.591865960000000e+02     1.591865960000000e+02
%      2.808522723000000e+00     2.808522723000000e+00
%      6.375581199000000e+01     6.375581199000000e+01
%      2.674026066000000e+01     2.674026066000000e+01
%      4.638532432000000e+01     4.638532432000000e+01
%      2.464521543000000e-01     2.464521543000000e-01
%      1.520484404000000e+01     1.520484404000000e+01
%      1.852266172000000e+00     1.852266172000000e+00
%      5.244639459000000e+01     5.244639459000000e+01
%      4.120394008000000e+01     4.120394008000000e+01
%      5.699317760000000e-01     5.699317760000000e-01
%      4.306056376000000e-01     4.306056376000000e-01
%      7.990620078300001e-03     7.990620078300001e-03
%      9.056036089000000e-01     9.056036089000000e-01
%      1.605425821600000e-02     1.605425821600000e-02
%      7.509759687000001e-01     7.509759687000001e-01
%      8.858285595500000e-02     8.858285595500000e-02
%      4.827726193000000e+01     4.827726193000000e+01
%      3.938459028000000e+01     3.938459028000000e+01
%      3.755297257000000e-01     3.755297257000000e-01
%      1.077562698000000e+02     1.077562698000000e+02
%      2.977250546000000e+01     2.977250546000000e+01
%      8.832481135000000e+01     8.832481135000000e+01
%      2.303929507000000e+01     2.303929507000000e+01
%      6.285848794000000e+01     6.285848794000000e+01
%      5.546318688000000e+00     5.546318688000000e+00
%      1.192244772000000e+01     1.192244772000000e+01
%      5.555448243000000e+00     5.555448243000000e+00
%      9.218489761999999e-01     9.218489761999999e-01
%      9.459927549000000e+01     9.459927549000000e+01
%      7.729698353000001e+01     7.729698353000001e+01
%      6.305263039000000e+01     6.305263039000000e+01
%      5.397970677000000e+01     5.397970677000000e+01
%      2.464355755000000e+01     2.464355755000000e+01
%      6.130192144000000e+01     6.130192144000000e+01
%      2.221000000000000e+01     2.221000000000000e+01
%      4.006374673000000e+01     4.006374673000000e+01
%      3.810034370000000e+01     3.810034370000000e+01
%      4.653415582000000e+01     4.653415582000000e+01
%      4.744573456000000e+01     4.744573456000000e+01
%      4.110581288000000e+01     4.110581288000000e+01
%      1.811349055000000e+01     1.811349055000000e+01
%      5.000000000000000e+01     5.000000000000000e+01
%   YP:
%     -9.311351132055279e-04    -2.245521557142600e-06 <- Sig Dif
%     -3.680818849716161e-07    -3.680818849716200e-07
%     -9.313903218526320e-04    -1.341411575595000e-06 <- Sig Dif
%     -6.614342966599907e-04    -4.106221069832800e-07 <- Sig Dif
%     -2.689011108429895e-04    -1.116128657940900e-06 <- Sig Dif
%     -1.312917689233473e-06    -9.964246872051801e-08 <- Sig Dif
%      6.609202557683602e-04    -1.843536381329600e-07 <- Sig Dif
%      2.688889922808357e-04    -5.530856128643800e-08 <- Sig Dif
%      5.904527482059052e-05    -1.406638006074000e-08 <- Sig Dif
%     -3.216924339355387e-06    -3.216923774473900e-06
%     -1.320602298626739e-06    -1.320602072141200e-06
%     -2.292647105051060e-06    -2.292646652080100e-06
%     -1.172087313872083e-07    -1.172087100709300e-07
%     -1.564706110457337e-06    -1.564705769396800e-06
%     -2.046840421598972e-07    -2.046840066327600e-07
%     -3.064993734369637e-07    -3.064992597501300e-07
%     -1.538568028536247e-07    -1.538566607450800e-07
%     -3.970321582613678e-08    -3.970320694435300e-08
%     -2.569517931760856e-10    -2.569517931760900e-10
%     -2.359501483084614e-13    -2.359501483084600e-13
%     -1.614086642121038e-10    -1.614086642121000e-10
%     -1.656258463711424e-11    -1.656258463711400e-11
%     -1.189860210359939e-09    -1.189860210359900e-09
%     -1.897565438113702e-11    -1.897565438113700e-11
%     -4.671457531912893e-08    -4.671457531912900e-08
%     -1.934833449013240e-08    -1.934833449013200e-08
%     -2.352361638813250e-10    -2.352361638813200e-10
%      5.300687917042524e-06     5.300687462295200e-06
%      1.688940415078832e-06     1.688940187705200e-06
%      3.474454388197046e-06     3.474453933449700e-06
%      4.445469699021487e-07     4.445469699021500e-07
%      2.613836613818421e-06     2.613836272757900e-06
%      3.049074734917667e-07     3.049074308592000e-07
%      5.932505473538185e-07     5.932504336669800e-07
%      3.048364050073360e-07     3.048362486879300e-07
%      5.976254513484491e-08     5.976253802941801e-08
%     -4.380748794538104e-07    -4.380756709460800e-07
%      3.318385615567584e-07     3.318384779400800e-07
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0
%                          0                         0

% End of TEFUNC

sys = YP;
end

% end mdlDerivatives


%
%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================
%
function sys = mdlOutputs(t,x,u)
global XMEAS YY
YY = x;
sys = XMEAS;

% end mdlOutputs

