function [f,state0] = gas_water_flow(r)
%% 流体性质参数
% ===================== Langmuir等温吸附模型====================
% VE = VL*P/(PL+P)
gas_prop.VL = 0.021; % Langmuir volume, m^3/kg
gas_prop.PL = 3.5; % Langmuir pressure, MPa
% Psc = 1;
% Zsc = 1;
% density_rock = 2000; % 
% ===================== Knudsen 扩散系数 ====================
gas_prop.Kn = 1; % Knudsen数
c1 =1; c2 =8/9;
gas_prop.Kn_modified_factor = 1+8*c1*gas_prop.Kn+16*c2*gas_prop.Kn^2;
% ===================== 高速非达西流 Forchheimer 方程 ====================
gas_prop.beta_non_Darcy_flow = 0; % 1e-8
% ===================== 启动压力梯度 ====================
p_grad_threshold = 0.00; % MPa/m
% =======================气相体积系数、粘度===============================
density_g_sc = 0.716;
% 第一种模式，直接输入参考压力、体积系数、压缩系数、粘度
gas_model =1;
if gas_model == 1
prg = 30; % reference pressure, MPa
Bgi = 0.01; % gas phase volume factor
cg = 1e-2 ;% gas phase compressibility, 1/MPa 
vgi = 0.01; % gas viscosity, cp
cvg = 0;
[Ppr,BG,MUG] = cal_gas_prop(prg,Bgi,cg,vgi,cvg);
end
if gas_model == 2
Ppr = [  0.1013
    2.0946
    4.0878
    6.0811
    8.0743
    10.0676
    12.0608
    14.0540
    16.0473
    18.0405
    20.0338
    22.0270
    24.0203
    26.0135
    28.0068
    30.0000
    ];
% BG = 0.01*ones(size(Ppr,1),1);
BG = [   1.18297
    0.0557041
    0.0278168
    0.0182594
    0.0134665
    0.0106154
    0.00874829
    0.00744907
    0.00650639
    0.0058006
    0.0052587
    0.0048336
    0.00449378
    0.00421751
    0.0039895
    0.00379871
    ];
% MUG = 0.01*ones(size(Ppr,1),1);
MUG = [0.0127683
    0.0130191
    0.0133973
    0.013872
    0.014437
    0.015088
    0.0158182
    0.0166173
    0.0174719
    0.0183671
    0.0192882
    0.0202219
    0.0211575
    0.0220865
    0.0230024
    0.0239008
    ];
end
% ========================水相体积系数、粘度===============================
density_w_sc = 1000;
prw = 30; % reference pressure
Bwi = 1.000; % water phase volume factor
cw = 4.0e-4 ;% water phase compressibility 
vwi = 1; % viscosity
cvw = 0;
% =============================基质相渗===============================
ifpcgl = 0; % 是否考虑毛管力
RPGW=[
    0.0000 	0.0000 	1.0000 	3.8916 
0.2000 	0.0000 	1.0000 	3.8916 
0.3160 	0.0002 	0.6784 	0.5796 
0.4350 	0.0004 	0.6215 	0.3724 
0.5620 	0.0010 	0.5456 	0.2425 
0.6140 	0.0020 	0.3939 	0.0608 
0.7020 	0.0280 	0.1399 	0.0372 
0.8120 	0.1721 	0.0515 	0.0137 
0.8750 	0.3395 	0.0297 	0.0104 
0.9060 	0.4395 	0.0226 	0.0090 
0.9370 	0.5500 	0.0173 	0.0075 
0.9690 	0.6702 	0.0131 	0.0059 
1.0 	1.0000 	0.0000 	0.0000
];
SW=RPGW(:,1);KRW=RPGW(:,2);KRG=RPGW(:,3);PCGL=RPGW(:,4);

% =============================裂缝相渗===============================
PRF=[
0.0000 	0.0000 	1.0000 	3.8916 
0.2000 	0.0000 	1.0000 	3.8916 
0.3160 	0.0002 	0.6784 	0.5796 
0.4350 	0.0004 	0.6215 	0.3724 
0.5620 	0.0010 	0.5456 	0.2425 
0.6140 	0.0020 	0.3939 	0.0608 
0.7020 	0.0280 	0.1399 	0.0372 
0.8120 	0.1721 	0.0515 	0.0137 
0.8750 	0.3395 	0.0297 	0.0104 
0.9060 	0.4395 	0.0226 	0.0090 
0.9370 	0.5500 	0.0173 	0.0075 
1.0 	1.0000 	0.0000 	0.0000
];
SWF = PRF(:,1);KRWF = PRF(:,2);KRGF = PRF(:,3);PCGLF=PRF(:,4);

%% 结合流体性质参数生成相应的数据体
f = fluidPVT_gas_water_flow(Ppr, BG, MUG, Bwi, prw, cw, vwi, cvw, SW, KRG, KRW, PCGL, SWF, KRGF, KRWF, PCGLF, density_g_sc, ifpcgl, r.rpt);
f.Dwsi = density_w_sc;
f.Dgsi = density_g_sc;
f.gas_prop = gas_prop;
f.p_grad_threshold = p_grad_threshold;
%% 初值条件 Initial Condition Section
% =================压力、饱和度初值======================================= %
%压力初值要考虑重力，给出油藏下表面压力值
P = [20 * ones(r.nmc, 1); 20 * ones(r.nfc, 1)];
% plow=25; % P =plow-1e-6*800*9.8*r.z;
Sw = [0.2 * ones(r.nmc, 1); 0.2 * ones(r.nfc, 1)];
state0 = initialRS_gas_water_flow(P, Sw);
end