xaj.txt

新安江模型是分散型性的水文流域模型

Tfrm_XinAnJiangModel = class(TForm)
  private
    { Private declarations }
  public
    { Public declarations }
    {--------------------------一些常用计算公式----------------------------}
    function RunoffToFlow(F,Hour ouble) ouble;
    function NetPre(K,P,EM ouble) :Double;     //NetPrecipitationEvaporation
    function MaxWSCap(WM,B,IMP :Double) :Double;//maximum water storage capacity
    function MaxFWSCap(SM,EX :Double) :Double;//maximum free water storage capacity
    function QRGroundWater(KKG,QRG0,RG,U :Double) :Double; //GroundWater地下水的计算
    function AveRainFall(StationRF,Weight :array of double) :Double;
    {------------------------蓄满产流函数---------------------------------}
    procedure xmcl(K,WM,B,IMP :Double; P,EM,W :Double; var R :Double); //蓄满产流
    procedure zsf(K,WM,WUM,WLM,C :Double; P,EM,WU0,WL0,WD0,R :Double;var WU1,WL1,WD1,E :Double);//蒸散发模型
    procedure syhf(K,IMP,WM,B,KSS,KG,SM,EX :Double; P,EM,W,R,S0 :Double; SDCount :Integer;
      var RS,RSS,RG,S1 :Double);//水源划分
    procedure UnitLine(n,K,Hour :Double; var U :array of double);//时段单位线
    procedure Confluence(I,U :array of Double; var Q :array of double);  //Confluence----汇流计算
    procedure SoilHumidity(rq :array of TDateTime; K,K4,K5,K6,K7,K8,WM,B,IMP,WUM,WLM,C,SM,EX,KSS,KG :Double;
      P,EM :array of double;SDCount :Integer; U0 :Double; var WU,WL,WD,S :array of double;
      var E,R,RS,RSS,RG :array of double); //计算土壤蓄水量及净雨量
    procedure XinAnJiang(rq :array of TDateTime; K,K4,K5,K6,K7,K8,WM,B,IMP,WUM,WLM,C,SM,EX,KSS,KG,KKG,NN,NK :Double; N:Integer;
        P,EM :array of double; SDCount :Integer; U0 :Double; var QRG0,FR,S,WU,WL,WD :Double; var Q :array of double);//新安江模型
    procedure XinAnJiang1(rq :array of TDateTime; K,K4,K5,K6,K7,K8,WM,B,IMP,WUM,WLM,C,SM,EX,KSS,KG,KKG,NN,NK :Double; N:Integer;
        P,EM :array of double; SDCount :Integer; U0 :Double; var QRG0,FR :Double; var S,WU,WL,WD,E,R,RS,RSS,RG,Q :array of double);//新安江模型
    procedure XinAnJiang2(rq :array of TDateTime; K,K4,K5,K6,K7,K8,WM,B,IMP,WUM,WLM,C,SM,EX,KSS,KG,KKG :Double; UH :array of double;
        N:Integer,EM :array of double; SDCount :Integer; U0,QJL :Double; var QRG0,FR,S,WU,WL,WD :Double; var Q :array of double);//新安江模型
    {-------------------------计算Γ(X)函数------------------------------}
    function Gamma(XX :Double) :Double; //计算Γ(X)
    function GammaLn(XX :Double) :Double; //计算lnΓ(X)
    procedure GammaSeries(var GammaSer :Double; A,X :Double; var GamLN :Double);  //Series级数   Γ分布函数的级数展式解
    procedure GammaConFrac(var GammaCF :Double; A,X :Double; var GamLN :Double);  //ConFrac-ContinuedFraction连分式 //Γ分布函数的连分式展式解
    function GammaP(A,X :Double) :Double;  //Γ分布函数P(a,x) a--n  x--t/K
    function GammaQ(A,X :Double) :Double;  //Γ分布函数余函数Q(a,x)    
  end;

var
  frm_XinAnJiangModel: Tfrm_XinAnJiangModel;

implementation
uses unt_csld;
{$R *.dfm}

{ TForm4 }

function Tfrm_XinAnJiangModel.AveRainFall(StationRF, Weight: array of double): Double;
var
  I :Integer;
  Total :Double;
begin
  Total :=0;
  for I :=1 to High(StationRF) do
    Total :=Total+StationRF*Weight;
  Result :=Total;
end;

//-------------流域最大自由蓄量
function Tfrm_XinAnJiangModel.MaxFWSCap(SM, EX: Double): Double;
begin
  Result :=SM*(1+EX);
end;

//-------------流域最大蓄量
function Tfrm_XinAnJiangModel.MaxWSCap(WM, B, IMP: Double): Double;
begin
  Result :=WM*(1+B)/(1-IMP);
end;

//-----------有效降雨
function Tfrm_XinAnJiangModel.NetPre(K, P, EM: Double): Double;
begin
  Result :=P-K*EM
end;

//--------------地下水汇流公式
function Tfrm_XinAnJiangModel.QRGroundWater(KKG, QRG0, RG, U: Double): Double;
begin
  Result :=QRG0*KKG+RG*(1-KKG)*U;
end;

//-------------单位转换系数，将径流深转化成流量
function Tfrm_XinAnJiangModel.RunoffToFlow(F, Hour: Double): Double; //F流域面积，单位为平方公里
begin
  Result :=F/(3.6*Hour);
end;

procedure Tfrm_XinAnJiangModel.syhf(K, IMP, WM, B, KSS, KG, SM, EX, P, EM,
  W, R, S0: Double; SDCount: Integer; var RS, RSS, RG, S1: Double);
var
  i,D,G :integer;
  KSSD,KGD,SMM,FR,AU :Double;
  X,PE,PE0 :Double;
  SS :array of Double;
  //S,Z1 :Double;
begin
  RS:=0; RSS:=0; RG:=0; S1:=0;
  PE :=NetPre(K,P,EM);
  if PE>0 then begin
    FR :=R/PE-IMP; //当有效降雨大于0时，产流面积的计算应是净雨深/有效降雨
    if Frac(PE/5.0)<>0 then D :=Trunc(PE/5.0)+1
    else D :=Trunc(PE/5.0);
    PE0 :=PE/D;  G :=D;
    KSSD :=(1-Power(1-(KG+KSS),1.0/G))/(1+KG/KSS);     //差分计算时，KSSD的计算公式？
    KGD :=KSSD*KG/KSS;
    SMM :=MaxFWSCap(SM,EX);
    SetLength(SS,D+1);
    SS[0] :=S0;
    //S :=S0;
    for i :=1 to D do begin
      //AU :=SMM*(1-Power(1-S0/SM,1.0/(1+EX)));
      AU :=SMM*(1-Power(1-SS[i-1]/SM,1.0/(1+EX)));
      if PE0+AU>=SMM then begin
        //RS :=(PE0+S0-SM)*FR+RS;
        RS :=(PE0+SS[i-1]-SM)*FR+RS;
        RSS :=SM*KSSD*FR+RSS;
        RG :=SM*KGD*FR+RG;
        //S0 :=(1-KSSD-KGD)*SM;
        SS :=(1-KSSD-KGD)*SM;
      end else begin
        X :=SM*Power(1-(PE0+AU)/SMM,1+EX);
        //RS :=(PE0-SM+S0+X)*FR+RS;
        RS :=(PE0-SM+SS[i-1]+X)*FR+RS;
        RSS :=(SM-X)*KSSD*FR+RSS;
        RG :=(SM-X)*KGD*FR+RG;
        //S0 :=(1-KSSD-KGD)*(SM-X);
        SS :=(1-KSSD-KGD)*(SM-X);
      end;
    end;//end for
    RS :=RS+IMP*PE; //地面径流包括透水面积上直接径流和不透水面积上直接径流
    //S1 :=S0;
    S1 :=SS[D];
    //Z1 :=RS+RSS+RG+S1*FR-S*FR;
    //Z1 :=RS+RSS+RG+S1*FR-S0*FR;
  end
  else begin   //PE<=0
    if (W>WM) and (W<WM+0.00001) then W :=WM;
    FR :=(1-IMP)*(1-Power(1-W/WM,B/(1+B)));
    G :=1;
    KSSD :=(1-Power(1-(KG+KSS),1.0/G))/(1+KG/KSS);     //差分计算时，KSSD的计算公式？
    KGD :=KSSD*KG/KSS;
    RS :=0;
    RSS :=S0*KSSD*FR;
    RG :=S0*KGD*FR;
    S1 :=(1-KSSD-KGD)*S0;
    //Z1 :=RS+RSS+RG+S1*FR-S0*FR;
  end;
  SS :=nil;
end;

procedure Tfrm_XinAnJiangModel.xmcl(K, WM, B, IMP, P, EM, W: Double;
  var R: Double);
var
  WMM,A,PE :Double;
begin
  //FR :=(1-IMP)*(1-Power(1-W/WM,B/(1+B))); //某时段透水面积的产流面积
  PE :=NetPre(K,P,EM);
  if PE>0 then
  begin
    WMM :=MaxWSCap(WM,B,IMP);
    if (W>WM) and (W<WM+0.00001) then W :=WM;
    A :=WMM*(1-Power(1-W/WM,1.0/(1+B)));    //防止由于计算机的精度问题使W可能微微大于120
    if (A+PE)<WMM then R :=PE-WM+W+WM*Power(1-(PE+A)/WMM,1+B)
      else R :=PE-(WM-W);
   //FR1 :=R/PE;   //表示平均净雨深面积
   // Memo1.Lines.Add('FR1='+floattostr(FR1));
  end else
  begin
    R :=0;
  end;
end;

procedure Tfrm_XinAnJiangModel.zsf(K, WM, WUM, WLM, C, P, EM, WU0, WL0,
  WD0, R: Double; var WU1, WL1, WD1, E: Double);
var
  //W0 :Double;
  EU,EL,ED,PE :Double;
begin
  //W0 :=WU0+WL0+WD0;
  if (WU0>WUM+0.00001) or (WL0>WLM+0.00001) or (WD0>WM-WUM-WLM+0.00001) then begin
    showmessage('WU0='+Floattostr(WU0)+' WL0='+floattostr(WL0)+' WD0='+floattostr(WD0));
    showmessage('WUM='+Floattostr(WUM)+' WLM='+floattostr(WLM)+' WDM='+floattostr(WM-WUM-WLM));
    showmessage('各层土壤蓄水量超过其容量，请核实');
    exit;
  end;
  if R> then begin
    showmessage('净雨深计算错误，其值应小于P');
    exit;
  end;
  PE :=NetPre(K,P,EM);
  if PE>0 then begin
  {PE>0}
    EU :=K*EM;
    EL :=0;
    ED :=0;//各层蒸散发值
    //E :=EU+EL+ED;
    //计算各层土壤蓄水量
    if WU0+PE-R>WUM then
      if WU0+WL0+PE-R-WUM>WLM then
      begin
        WU1 :=WUM;
        WL1 :=WLM;
        WD1 :=WU0+WL0+WD0+PE-R-WUM-WLM;
        if WD1<0 then //------------------------------------土壤蓄水量小于0，不可能发生------------------------------------
      end
      else begin
        WU1 :=WUM;
        WL1 :=WU0+WL0+PE-R-WUM;
        WD1 :=WD0;
      end
    else begin
      WU1 :=WU0+PE-R;
      WL1 :=WL0;
      WD1 :=WD0;
    end;  //PE>0蒸散发计算完毕
  {PE>0}
  end else begin
  {PE<0}
    if R>0 then begin showmessage('PE<0,这时R应该等于0,退出');exit;end;
    if WU0>ABS(PE) then begin
      EU :=K*EM;
      EL :=0;
      ED :=0;
      WU1 :=WU0+PE;
      WL1 :=WL0;
      WD1 :=WD0;
    end else begin
      EU :=WU0+P;
      WU1 :=0;
      if WL0>C*WLM then begin //WL0>C*WLM并不代表EL将小于WL0，当EM很大时，EL将大于WL0，所得WL1将为负值
        EL :=(K*EM-EU)*WL0/WLM;   //只有当(K*EM-EU)>WLM才会出现WL0<EL，此种情况几乎不可能发生
        ED :=0;
        WL1 :=WL0-EL;
        if WL1<0 then showmessage('WL1<0');
        WD1 :=WD0;
      end else begin
        if WL0 >C*(K*EM-EU) then begin//由于EM通常在20以内，故C*(K*EM-EU)值很小，通常不大于4
          EL :=C*(K*EM-EU);
          ED :=0;
          WL1 :=WL0-EL;
          WD1 :=WD0;
        end else begin //此种情况很难发生，除非EM(10-20)很大，WU0，WL0很小(5mm以下)
          EL :=WL0;
          ED :=C*(K*EM-EU)-EL;
          WL1 :=0;
          WD1 :=WD0-ED;
          if WD1<0 then WD1:=0;  //------------------------------------土壤蓄水量小于0，不可能发生--------------------------------------------------------------
        end;//end WL0>C*(K*EM-EU)
      end;//end WL0>C*WLM
    end;//end WU0>ABS(PE)
  {PE<0}
  end; //蒸散发模型计算完毕 返回EU,EL,ED,WU1,WL1,WD1
  E :=EU+EL+ED;
end;

//----------------计算土壤蓄水量及净雨量
//P、EM-N个时段降雨和蒸发能力值，U0-将净雨深转化为流量，S-初期自由蓄水量  SDCount-将一天分为几个时段
//WU、WL、WD-初始值为初期土壤各层蓄水量、返回末期土壤各层蓄水量，QRG0-起始时刻流量，返回末时段流量
procedure Tfrm_XinAnJiangModel.SoilHumidity(rq :array of TDateTime; K,K4,K5,K6,K7,K8, WM, B, IMP, WUM, WLM, C, SM,
  EX, KSS, KG: Double; P, EM: array of double; SDCount: Integer; U0: Double;
  var WU, WL, WD, S, E, R, RS, RSS, RG: array of double);
var
  Year,Month,Day :Word;
  i :Integer;
  //RS,RSS,RG :array of double;
  S1,K1 :Double;
begin
  //SetLength(RS,High(P)+1);    SetLength(RSS,High(P)+1);   SetLength(RG,High(P)+1);
  for i :=Low(P)+1 to High(P) do begin
    DecodeDate(rq,Year,Month,Day);
    case Month of
      4 :K1:=K4;
      5 :K1 :=K5;
      6 :K1:=K6;
      7 :K1:=K7;
      8 :K1:=K8;
      else K1:=K;
    end;
    xmcl(K1,WM,B,IMP,P,EM,WU[i-1]+WL[i-1]+WD[i-1],R);
    zsf(K1,WM,WUM,WLM,C,P,EM,WU[i-1],WL[i-1],WD[i-1],R,WU,WL,WD,E);
    syhf(K1,IMP,WM,B,KSS,KG,SM,EX,P,EM,WU[i-1]+WL[i-1]+WD[i-1],R,S[i-1],SDCount,RS,RSS,RG,S1);
    S :=S1;
  end; //end if
end;

procedure Tfrm_XinAnJiangModel.XinAnJiang(rq :array of TDateTime; K,K4,K5,K6,K7,K8, WM, B, IMP, WUM, WLM, C, SM,
  EX, KSS, KG, KKG, NN, NK: Double; N: Integer; P, EM: array of double;
  SDCount: Integer; U0: Double; var QRG0, FR, S, WU, WL, WD: Double;
  var Q: array of double);
var
  Year,Month,Day :Word;
  i,m :Integer;
  R,RS,RSS,RG :array of double;
  WU0,WL0,WD0,E :array of double;
  IU,QRS,QRG :array of double; //IU=RS+RSS--净雨
  Hour,S0,S1,QJL,K1 :Double; //QJL-基流
  UH :array of double;
begin
  //----计算各时段的RS、RSS、RG
  SetLength(R,N+1);     SetLength(RS,N+1);    SetLength(RSS,N+1);   SetLength(RG,N+1);
  SetLength(WU0,N+2);   SetLength(WL0,N+2);   SetLength(WD0,N+2);   SetLength(E,N+1);
  SetLength(IU,N+1);    SetLength(QRS,N+1);   SetLength(QRG,N+1);
  S0:=S; WU0[1]:=WU; WL0[1]:=WL; WD0[1]:=WD;
  for i :=1 to N do begin
    DecodeDate(rq,Year,Month,Day);
    case Month of
      4 :K1:=K4;
      5 :K1 :=K5;