Слияние кода завершено, страница обновится автоматически
using System;
using System.Collections.Generic;
using System.Text;
using System.IO;
using System.Windows.Forms;
using System.Drawing;
namespace GeoFly
{
public class RouteFlux
{
/// <summary>
/// 栅格上的入流量
/// </summary>
public double InFlux;
/// <summary>
/// 存放一个栅格上的汇流输出量
/// </summary>
public double OutFlux;
public bool bCal;
public RouteFlux()
{
InFlux = 0.0;
OutFlux = 0.0;
bCal = false;
}
}
//水文模拟总控制类
public class HydroSimulate
{
public HydroSimulate()
{
}
public HydroSimulate(DEMRiverNet river)
{
this.River = river;
}
public DEMRiverNet River;
public static GridLayer GridDem = new GridLayer();
public static GridLayerPara g_GridLayerPara = new GridLayerPara();
/// <summary>
/// 气象参数对象
/// </summary>
public static ClimatePara g_ClimatePara = new ClimatePara();
public static MidGridResultOut g_GridResultOut = new MidGridResultOut();
public static ModelRunPara g_ModelRunPara = new ModelRunPara();
public static ModelInputPara g_ModelInputPara = new ModelInputPara();
public static string path="C:\\ESSI_Data";
/// <summary>
/// 汇流次序表
/// </summary>
//public List<Cell> River.CalcuOrderList = new List<Cell>();
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ 径流过程模拟公用函数定义 +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ 控制： +
+ 1、数据读入； +
+ 2、内存初始化； +
+ 3、结果输出； +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/// <summary>
/// 读入栅格图层参数（DEM,LULC,SOIL），该函数已修改
/// </summary>
/// <returns></returns>
public void ReadInGridLayerData()
{
//读入DEM数据
string filename = g_GridLayerPara.DemFileName;
g_GridLayerPara.g_DemLayer = new GridLayer();
g_GridLayerPara.g_DemLayer.ReadASC(filename);
//读入土地利用数据
filename = g_GridLayerPara.LulcFileName;
g_GridLayerPara.g_VegLayer = new GridLayer();
g_GridLayerPara.g_VegLayer.ReadASC(filename);
//读入土壤类型数据
filename = g_GridLayerPara.SoilFileName;
g_GridLayerPara.g_SoilLayer = new GridLayer();
g_GridLayerPara.g_SoilLayer.ReadASC(filename);
}
/// <summary>
/// 读入气象数据参数
/// </summary>
/// <returns></returns>
public bool ReadInClimateData()
{
MapExtent extent = g_GridLayerPara.g_DemLayer.Extent;
double resolution = g_GridLayerPara.g_DemLayer.resolution;
g_ClimatePara.pcpdata.SetExtentInfo(extent, resolution);
g_ClimatePara.slrdata.SetExtentInfo(extent, resolution);
g_ClimatePara.hmddata.SetExtentInfo(extent, resolution);
g_ClimatePara.tmpmxdata.SetExtentInfo(extent, resolution);
g_ClimatePara.tmpmndata.SetExtentInfo(extent, resolution);
g_ClimatePara.tmpmeandata.SetExtentInfo(extent, resolution);
g_ClimatePara.WindData.SetExtentInfo(extent, resolution);
g_ClimatePara.petdata.SetExtentInfo(extent, resolution);
//读入站点位置数据和值
g_ClimatePara.pcpdata.ReadInStationPos(g_ClimatePara.PcpInfo.StationFileName);
g_ClimatePara.pcpdata.ReadInStationValue(g_ClimatePara.PcpInfo.DataFileName);
int count = g_ClimatePara.pcpdata.saDate.Count;
DateTime startDate = g_ClimatePara.pcpdata.saDate[0]; //第一天
DateTime endDate = g_ClimatePara.pcpdata.saDate[count - 1];//最后一天
g_ClimatePara.StartDate = startDate;
g_ClimatePara.EndDate = endDate;
//如果蒸散发有现成的实际测量值，则读入实际测量值
if (g_ModelRunPara.PETMethod == PET_METHOD.PET_REAL)
{
g_ClimatePara.petdata.ReadInStationPos(g_ClimatePara.PetInfo.StationFileName);
g_ClimatePara.petdata.ReadInStationValue(g_ClimatePara.PetInfo.DataFileName);
}
//如果为暴雨
if (g_ModelRunPara.RunoffSimuType == RunOffSimuType.STORM_RUNOFF_SIMULATION)
return true;
if (g_ClimatePara.SlrInfo.StationFileName != null && g_ClimatePara.SlrInfo.StationFileName != "")
{
g_ClimatePara.slrdata.ReadInStationPos(g_ClimatePara.SlrInfo.StationFileName);
g_ClimatePara.slrdata.ReadInStationValue(g_ClimatePara.SlrInfo.DataFileName);
}
if (g_ClimatePara.HmdInfo.StationFileName != null && g_ClimatePara.HmdInfo.StationFileName != "")
{
g_ClimatePara.hmddata.ReadInStationPos(g_ClimatePara.HmdInfo.StationFileName);
g_ClimatePara.hmddata.ReadInStationValue(g_ClimatePara.HmdInfo.DataFileName);
}
if (g_ClimatePara.WndInfo.StationFileName != null && g_ClimatePara.WndInfo.StationFileName!="")
{
g_ClimatePara.WindData.ReadInStationPos(g_ClimatePara.WndInfo.StationFileName);
g_ClimatePara.WindData.ReadInStationValue(g_ClimatePara.WndInfo.DataFileName);
}
if (g_ClimatePara.TempmxInfo.StationFileName != null && g_ClimatePara.TempmxInfo.StationFileName != "")
{
g_ClimatePara.tmpmxdata.ReadInStationPos(g_ClimatePara.TempmxInfo.StationFileName);
g_ClimatePara.tmpmxdata.ReadInStationValue(g_ClimatePara.TempmxInfo.DataFileName);
g_ClimatePara.tmpmndata.ReadInStationPos(g_ClimatePara.TempmxInfo.StationFileName);
g_ClimatePara.tmpmndata.ReadInStationValue(g_ClimatePara.TempmnInfo.DataFileName);
g_ClimatePara.tmpmeandata.ReadInStationPos(g_ClimatePara.TempmxInfo.StationFileName);
if (g_ClimatePara.TempmeanInfo.DataFileName != null && g_ClimatePara.TempmeanInfo.DataFileName !="")
{
g_ClimatePara.tmpmeandata.ReadInStationValue(g_ClimatePara.TempmeanInfo.DataFileName);
}
}
return true;
}
/// <summary>
/// 读入马斯京根汇流参数
/// </summary>
/// <returns></returns>
public bool ReadMuskingCoeff(string m_MuskCoeffFile)
{
if (File.Exists(m_MuskCoeffFile))
{
MessageBox.Show("马斯京根法（先演后合）汇流文件 [muskingum_coeff.txt] 不存在。");
return false;
}
List<string> saValue = MatrixFuncs.FileRead(m_MuskCoeffFile);
this.m_iNodeNum = Convert.ToInt32(saValue[0]);
if (m_iNodeNum > 1)
{
m_pX = new double[m_iNodeNum];
m_pK = new double[m_iNodeNum];
int id = 0;
for (int i = 2; i < saValue.Count; i++)
{
string[] saOut = saValue[i].Split(new char[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
m_pX[id] = Convert.ToDouble(saOut[1]);
m_pK[id] = Convert.ToDouble(saOut[2]);
id++;
}
}
return true;
}
public bool ReadWaterYearType(string WaterYearTypeFile) //读入流域丰、平、枯水年份信息
{
List<string> saValue = MatrixFuncs.FileRead(WaterYearTypeFile);
int num = saValue.Count;
if (num > 0)
{
waterYearType = new SortedList<int, WaterYearType>();
for (int i = 1; i < num; i++)
{
string[] saOut = saValue[i].Split(new char[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
int key = Convert.ToInt32(saOut[0]);
waterYearType[key] = (WaterYearType)Convert.ToInt32(saOut[1]);
}
return true;
}
else
return false;
}
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ 长时段降雨～径流过程模拟函数定义 +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/// <summary>
/// 长时段降雨径流模拟主函数
/// </summary>
public void LongTermRunoffSimulate()
{
if (g_ModelRunPara.SurfRouteMethod == SurfRouteMehtod.ROUTE_MUSK_CONGE)
{
//读入栅格汇流最优次序参数文件
this.ReadInRoutingPara();
}
this.ReadInRoutingLayerData();
if (g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_MUSKINGUM_COMBINE_FIRST ||
g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_MUSKINGUM_ROUTE_FIRST)
{
OpenFileDialog dlg = new OpenFileDialog();
if (dlg.ShowDialog() == DialogResult.OK)
{
ReadMuskingCoeff(dlg.FileName);
}
}
int m_row = g_GridLayerPara.g_DemLayer.rowCount;
int m_col = g_GridLayerPara.g_DemLayer.colCount;
int totrec = 0; //= GetTimeInterval(true);//true表明按整日算
this.MyOutletQ = new SortedList<DateTime, OutletQ>();
if (m_iNodeNum > 1)
{
this.m_pNodeSurfQ = new SortedList<DateTime, double[]>();// new double[totrec, g_GridLayerPara.m_subNum];// GridOutletInit(totrec, g_GridLayerPara.m_subNum, 0);
this.m_pNodeLatQ = new SortedList<DateTime, double[]>();// new double[totrec, g_GridLayerPara.m_subNum]; //GridOutletInit(totrec, g_GridLayerPara.m_subNum, 0);
this.m_pNodeBaseQ = new SortedList<DateTime, double[]>();// new double[totrec, g_GridLayerPara.m_subNum]; //GridOutletInit(totrec, g_GridLayerPara.m_subNum, 0);
this.m_pNodeOutQ = new SortedList<DateTime, double[]>();// new double[totrec, g_GridLayerPara.m_subNum]; //GridOutletInit(totrec, g_GridLayerPara.m_subNum, 0);
}
double dthet1334, dCp;
double dsnowfactor = 1.0;
double aetfactor = 0;
int TotalYears = 0;
int dn;
dCp = 0.5;
int StartYear, EndYear, startday, endday/*,mon,day*/;
StartYear = g_ClimatePara.StartDate.Year;
EndYear = g_ClimatePara.EndDate.Year;
startday = DateAndTime.GetDnInYear(g_ClimatePara.StartDate.Year, g_ClimatePara.StartDate.Month, g_ClimatePara.StartDate.Day);
endday = DateAndTime.GetDnInYear(g_ClimatePara.EndDate.Year, g_ClimatePara.EndDate.Month, g_ClimatePara.EndDate.Day);
DateTime curorder = new DateTime(StartYear, g_ClimatePara.StartDate.Month, startday);
TotalYears = EndYear - StartYear + 1;
if (!ReadWaterYearType(this.m_WaterYrTypeFile))
{
waterYearType = new SortedList<int, WaterYearType>();
for (int i = 0; i < TotalYears; i++)
{
int key = StartYear + i;
waterYearType[key] = WaterYearType.WATER_LOW_YEAR;
}
}
//水文过程年循环
for (int year = StartYear; year <= EndYear; year++)
{
//按年份取得该年的模拟时段
int ibeg = 0;
int iend = 1;
this.GetBEDateInYear(year, StartYear, EndYear, startday, endday, ibeg, iend);
//从该年的模拟时段的开始日到结束日
for (int j = ibeg; j <= iend; j++)
{
dn = j;
PcpDataInterp(curorder);
MeteoDataInterp(curorder);
int iMonth = DateAndTime.GetMonthByDn(year, j);
double dintensity = 0;
double dhr = 24;
double dhrIntensity = 0;
double dPE = 0;
//对每个网格做下面的处理
for (int row = 0; row < m_row; row++)
{
for (int col = 0; col < m_col; col++)
{ //如果该网格已不需要计算，则跳过
if (!IfGridNeedCal(row, col))
continue;
dsnowfactor = 1;
//如果已读取了平均气温则执行
//如果温度低于雪温，则将降水加入雪量中
if (g_ClimatePara.tmpmeandata.m_gridLayer[row, col] < g_ModelInputPara.SnowTemperature)
{
g_GridResultOut.m_SnowWater[row, col] += g_ClimatePara.pcpdata.m_gridLayer[row, col];
//降雨量设为0
g_ClimatePara.pcpdata.m_gridLayer[row, col] = 0;
dsnowfactor = 0.15;
}
else //温度高于雪温时
{
if (g_GridResultOut.m_SnowWater[row, col] > 0) //如果雪量大于0
{
//求融雪量
double smelt = DDFSnowMelt(g_ClimatePara.tmpmeandata.m_gridLayer[row, col],
g_ModelInputPara.SnowTemperature, g_ModelInputPara.DDF,
g_ModelInputPara.DailyMeanPcpTime);
//如果实际雪量小于融雪量，就把融雪量设为实际雪量
if (g_GridResultOut.m_SnowWater[row, col] < smelt)
{
smelt = g_GridResultOut.m_SnowWater[row, col];
g_GridResultOut.m_SnowWater[row, col] = 0;//雪量设为0
}
else
g_GridResultOut.m_SnowWater[row, col] -= smelt;
//降水量中加入融雪量
g_ClimatePara.pcpdata.m_gridLayer[row, col] += smelt;
dsnowfactor = 0.3;
}
}
//每日平均降水时间(小时数)
dhrIntensity = g_ModelInputPara.DailyMeanPcpTime;
//每小时降水量
dintensity = g_ClimatePara.pcpdata.m_gridLayer[row, col] / dhrIntensity;
//设置网格参数
HortonInfil hortonInfil = new HortonInfil();
hortonInfil.SetGridPara(g_GridLayerPara, row, col, g_GridResultOut.m_SoilProfileWater.Values[row,col]);
//计算Horton超渗产流量
//土壤水下渗量
m_drateinf[row, col] = hortonInfil.HortonExcessRunoff();
//计算本栅格上的土壤水赤字
int SoilID = (int)g_GridLayerPara.g_SoilLayer[row, col];
dthet1334 = g_GridLayerPara.SoilTypes[SoilID].SoilWaterDeficitContent(g_GridResultOut.m_SoilProfileWater.Values[row,col]);
GridWaterBalance gridwb = new GridWaterBalance();
gridwb.SetGridPara(row, col, dintensity,m_drateinf[row,col], curorder, dhrIntensity);
int VegID = (int)g_GridLayerPara.g_VegLayer[row, col];
double dalb = GetVegAlbedo(VegID, iMonth, 1);//获取当前月的反照率
gridwb.CalcPET(g_ModelRunPara, g_ClimatePara, dalb);
if (g_GridLayerPara.g_StrahlerRivNet[row, col] != 0)
{
g_GridResultOut.m_GridSurfQ[row, col] = g_ClimatePara.pcpdata.m_gridLayer[row, col] - gridwb.m_dPET;
if (g_GridResultOut.m_GridSurfQ[row, col] < 0)
g_GridResultOut.m_GridSurfQ[row, col] = 0;
if (g_GridResultOut.m_GridSurfQ[row, col] > 1e+10)
MessageBox.Show("有错误");//AfxMessageBox("hello");
g_GridResultOut.m_GridLateralQ[row, col] = 0;
g_GridResultOut.m_GridBaseQ[row, col] = 0;
g_GridResultOut.m_GridTotalQ[row, col] = g_GridResultOut.m_GridSurfQ[row, col];
g_GridResultOut.m_AET[row, col] = gridwb.m_dPET;
g_GridResultOut.m_CI[row, col] = 0;
g_GridResultOut.m_CIDefict[row, col] = 0;
g_GridResultOut.m_NetPcp[row, col] = g_GridResultOut.m_GridSurfQ[row, col];
g_GridResultOut.m_GridWaterYieldType[row, col] = 0;
g_GridResultOut.m_SoilProfileWater[row, col] = 0;
g_GridResultOut.m_SoilAvgWater[row, col] = 0;
}
else
{
gridwb.CalcCI();
g_GridResultOut.m_CI[row, col] = gridwb.m_dCrownInterc;
if (g_GridResultOut.m_SoilProfileWater.Values[row,col] / g_GridLayerPara.SoilTypes[SoilID].SP_Fc > 0.8)
{
if (g_ClimatePara.pcpdata.m_gridLayer[row, col] > 0)
aetfactor = 0.6;
else
aetfactor = 0.9;
}
else
{
if (g_ClimatePara.pcpdata.m_gridLayer[row, col] > 0)
aetfactor = 0.4;
else
aetfactor = 0.6;
}
//*************对蒸散发处理的特殊代码段 -- 计算实际蒸散发**************//
if (g_ModelRunPara.PETMethod == PET_METHOD.PET_REAL)
{
gridwb.m_dAET = g_ClimatePara.petdata.m_gridLayer[row, col] * aetfactor;
g_GridResultOut.m_AET[row, col] = gridwb.m_dAET;
}
else
{
gridwb.CalcAET(g_ModelRunPara, g_ClimatePara, dalb);
if (gridwb.m_dAET > g_ClimatePara.petdata.m_gridLayer[row, col]) // FRZZ MODI 2005-8-18
g_GridResultOut.m_AET[row, col] = g_ClimatePara.petdata.m_gridLayer[row, col] *
Math.Exp(-1 * g_ClimatePara.petdata.m_gridLayer[row, col] / gridwb.m_dAET);
if (gridwb.m_dAET < 0)
g_GridResultOut.m_AET[row, col] = g_ClimatePara.petdata.m_gridLayer[row, col] * 0.1;
}
//**********************对蒸散发处理的特殊代码段**********************//
gridwb.CalcNetRain();
g_GridResultOut.m_CIDefict[row, col] = gridwb.m_dCIDeficit;
g_GridResultOut.m_NetPcp[row, col] = gridwb.m_dNetRain;
g_GridResultOut.m_GridWaterYieldType[row, col] = gridwb.CalcRunoffElement(g_GridLayerPara, g_ModelInputPara,g_GridResultOut.m_SoilProfileWater[row, col]);
g_GridResultOut.m_SoilProfileWater[row, col] = gridwb.Cell_SP_Sw; //g_GridLayerPara.SoilTypes[SoilID].SP_Sw;
g_GridResultOut.m_SoilAvgWater[row, col] = g_GridLayerPara.SoilTypes[SoilID].SoilAvgWater(gridwb.Cell_SP_Sw);
g_GridResultOut.m_GridTotalQ[row, col] = gridwb.m_dTotalQ;
g_GridResultOut.m_GridSurfQ[row, col] = gridwb.m_dSurfQ;
g_GridResultOut.m_GridLateralQ[row, col] = gridwb.m_dLateralQ;
g_GridResultOut.m_GridBaseQ[row, col] = gridwb.m_dBaseQ;
if (g_GridResultOut.m_GridSurfQ[row, col] > 1e+10)
MessageBox.Show("hello2");
}
}
}
if (m_iNodeNum == 1 || g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_PURE_LAG)
{
PureLagGridRouting(g_GridResultOut.m_GridSurfQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_SURFQ, curorder, dsnowfactor, waterYearType[year]);
PureLagGridRouting(g_GridResultOut.m_GridLateralQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_LATERALQ, curorder, dsnowfactor, waterYearType[year]);
PureLagGridRouting(g_GridResultOut.m_GridBaseQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_BASEQ, curorder, dsnowfactor, waterYearType[year]);
this.MyOutletQ[curorder].m_pOutletDeepBaseQ = DeepBaseQSim(dn, g_ModelInputPara.DeepBaseQ);
this.MyOutletQ[curorder].CalcOutletQ();
}
else
{
PureLagGridRouting_Node(g_GridResultOut.m_GridSurfQ, m_pNodeSurfQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_SURFQ, curorder, totrec, dsnowfactor, waterYearType[year]);
PureLagGridRouting_Node(g_GridResultOut.m_GridLateralQ, m_pNodeLatQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_LATERALQ, curorder, totrec, dsnowfactor, waterYearType[year]);
PureLagGridRouting_Node(g_GridResultOut.m_GridBaseQ, m_pNodeBaseQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_BASEQ, curorder, totrec, dsnowfactor, waterYearType[year - StartYear]);
double dSurf, dLat, dBase;
dSurf = dLat = dBase = 0;
//for(int i=0;i<g_GridLayerPara.m_subNum;i++)
//{
// dSurf += m_pNodeSurfQ[curorder,i] * g_ModelInputPara1.SurfQLinearFactor;
// dLat += m_pNodeLatQ[curorder,i];
// dBase += m_pNodeBaseQ[curorder,i];
//}
this.MyOutletQ[curorder].m_pOutletSurfQ = dSurf;
this.MyOutletQ[curorder].m_pOutletLatQ = dLat;
this.MyOutletQ[curorder].m_pOutletBaseQ = dBase;
this.MyOutletQ[curorder].m_pOutletDeepBaseQ = DeepBaseQSim(dn, g_ModelInputPara.DeepBaseQ);
if (g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_MUSKINGUM_COMBINE_FIRST)
{
for (int i = 0; i < g_GridLayerPara.m_subNum; i++)
{
m_pNodeOutQ[curorder][i] = m_pNodeSurfQ[curorder][i] * g_ModelInputPara.SurfQLinearFactor
+ m_pNodeLatQ[curorder][i] + m_pNodeBaseQ[curorder][i];
}
MuskingumRiverRouting(24, m_pNodeOutQ, m_pX, m_pK, g_GridLayerPara.m_subNum, curorder);
this.MyOutletQ[curorder].m_pOutletQ = pRiverRoute[g_GridLayerPara.m_subNum - 1].pRoute.OutFlux
+ this.MyOutletQ[curorder].m_pOutletDeepBaseQ;
}
else if (g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_MUSKINGUM_ROUTE_FIRST)
{
;
}
else
{
this.MyOutletQ[curorder].CalcOutletQ();
}
}
TimeSpan DaySpan = new TimeSpan(1, 0, 0, 0);
curorder = curorder + DaySpan;
DateTime dateTime = new DateTime(year, iMonth, j);
MidGridResultOut(dateTime, true);
}
}
RiverOutletQ(true);
m_iTotalRec = totrec;
}
/// <summary>
/// 计算栅格上的snowfactor
/// </summary>
/// <param name="row"></param>
/// <param name="col"></param>
/// <returns></returns>
private double CalcSnowFactor(int row, int col)
{
double SnowFactor = 1;
//如果温度低于雪温，则将降水加入雪量中
if (g_ClimatePara.tmpmeandata.m_gridLayer[row, col] < g_ModelInputPara.SnowTemperature)
{
g_GridResultOut.m_SnowWater[row, col] += g_ClimatePara.pcpdata.m_gridLayer[row, col];
//降雨量设为0
g_ClimatePara.pcpdata.m_gridLayer[row, col] = 0;
SnowFactor = 0.15;
}
else //温度高于雪温时
{
if (g_GridResultOut.m_SnowWater[row, col] > 0) //如果雪量大于0
{
//求融雪量
double smelt = DDFSnowMelt(g_ClimatePara.tmpmeandata.m_gridLayer[row, col],
g_ModelInputPara.SnowTemperature, g_ModelInputPara.DDF,
g_ModelInputPara.DailyMeanPcpTime);
//如果实际雪量小于融雪量，就把融雪量设为实际雪量
if (g_GridResultOut.m_SnowWater[row, col] < smelt)
{
smelt = g_GridResultOut.m_SnowWater[row, col];
g_GridResultOut.m_SnowWater[row, col] = 0;//雪量设为0
}
else
g_GridResultOut.m_SnowWater[row, col] -= smelt;
//降水量中加入融雪量
g_ClimatePara.pcpdata.m_gridLayer[row, col] += smelt;
SnowFactor = 0.3;
}
}
return SnowFactor;
}
/// <summary>
/// 长时段降雨径流模拟主函数(刘永和修改)
/// </summary>
public void LongTermRunoffSimulate_Liu()
{
if (g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_MUSKINGUM_COMBINE_FIRST ||
g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_MUSKINGUM_ROUTE_FIRST)
{
OpenFileDialog dlg = new OpenFileDialog();
if (dlg.ShowDialog() == DialogResult.OK)
{
ReadMuskingCoeff(dlg.FileName);
}
}
int rowCount = g_GridLayerPara.g_DemLayer.rowCount;
int colCount = g_GridLayerPara.g_DemLayer.colCount;
int totrec = 0; //= GetTimeInterval(true);//true表明按整日算
this.MyOutletQ = new SortedList<DateTime, OutletQ>();
if (m_iNodeNum > 1)
{
this.m_pNodeSurfQ = new SortedList<DateTime, double[]>();
this.m_pNodeLatQ = new SortedList<DateTime, double[]>();
this.m_pNodeBaseQ = new SortedList<DateTime, double[]>();
this.m_pNodeOutQ = new SortedList<DateTime, double[]>();
}
int TotalYears = g_ClimatePara.StartDate.Year - g_ClimatePara.EndDate.Year + 1;
if (!this.ReadWaterYearType(this.m_WaterYrTypeFile))
{
waterYearType = new SortedList<int, WaterYearType>();
for (int year = g_ClimatePara.StartDate.Year; year <= g_ClimatePara.EndDate.Year; year++)
{
waterYearType[year] = WaterYearType.WATER_LOW_YEAR;
}
}
//水文过程年循环
for (DateTime curDate = g_ClimatePara.StartDate; curDate <= g_ClimatePara.EndDate; curDate.Add(new TimeSpan(1, 0, 0, 0)))
{
PcpDataInterp(curDate);
MeteoDataInterp(curDate);
//int iMonth = curDate.Month;
double dintensity = 0;
double dhr = 24;
double dhrIntensity = 0;
double dPE = 0;
double dsnowfactor = 1.0;
//对每个网格做下面的处理
for (int row = 0; row < rowCount; row++)
{
for (int col = 0; col < colCount; col++)
{ //如果该网格已不需要计算，则跳过
if (!IfGridNeedCal(row, col))
continue;
dsnowfactor = this.CalcSnowFactor(row, col);
//每日平均降水时间(小时数)
dhrIntensity = g_ModelInputPara.DailyMeanPcpTime;
//每小时降水量
dintensity = g_ClimatePara.pcpdata.m_gridLayer[row, col] / dhrIntensity;
HortonInfil hortonInfil = new HortonInfil();
//设置网格参数
hortonInfil.SetGridPara(g_GridLayerPara, row, col,g_GridResultOut.m_SoilProfileWater.Values[row,col]);
//土壤水下渗量
m_drateinf[row, col] = hortonInfil.HortonExcessRunoff();
//计算本栅格上的土壤水赤字
int SoilID = (int)g_GridLayerPara.g_SoilLayer[row, col];
double dthet = g_GridLayerPara.SoilTypes[SoilID].SoilWaterDeficitContent(g_GridResultOut.m_SoilProfileWater.Values[row,col]);
GridWaterBalance gridwb = new GridWaterBalance();
gridwb.SetGridPara(row, col, dintensity,m_drateinf.Values[row,col], curDate, dhrIntensity);
int VegID = (int)g_GridLayerPara.g_VegLayer[row, col];
double dalb = GetVegAlbedo(VegID, curDate.Month, 1);//获取当前月的反照率
double PET = gridwb.CalcPET(g_ModelRunPara, g_ClimatePara, dalb);
//如果当前栅格为河道
if (g_GridLayerPara.g_StrahlerRivNet[row, col] != 0)
{
g_GridResultOut.m_GridSurfQ[row, col] = g_ClimatePara.pcpdata.m_gridLayer[row, col] - gridwb.m_dPET;
if (g_GridResultOut.m_GridSurfQ[row, col] < 0)
g_GridResultOut.m_GridSurfQ[row, col] = 0;
if (g_GridResultOut.m_GridSurfQ[row, col] > 1e+10)
MessageBox.Show("有错误");
g_GridResultOut.m_GridLateralQ[row, col] = 0;
g_GridResultOut.m_GridBaseQ[row, col] = 0;
g_GridResultOut.m_GridTotalQ[row, col] = g_GridResultOut.m_GridSurfQ[row, col];
g_GridResultOut.m_AET[row, col] = PET;
g_GridResultOut.m_CI[row, col] = 0;
g_GridResultOut.m_CIDefict[row, col] = 0;
g_GridResultOut.m_NetPcp[row, col] = g_GridResultOut.m_GridSurfQ[row, col];
g_GridResultOut.m_GridWaterYieldType[row, col] = 0;
g_GridResultOut.m_SoilProfileWater[row, col] = 0;
g_GridResultOut.m_SoilAvgWater[row, col] = 0;
}
else //当前栅格不是河道
{
g_GridResultOut.m_CI[row, col] = gridwb.CalcCI();
double AetFactor = 0;
if (g_GridResultOut.m_SoilProfileWater.Values[row,col] / g_GridLayerPara.SoilTypes[SoilID].SP_Fc > 0.8)
{
if (g_ClimatePara.pcpdata.m_gridLayer[row, col] > 0)
AetFactor = 0.6;
else
AetFactor = 0.9;
}
else
{
if (g_ClimatePara.pcpdata.m_gridLayer[row, col] > 0)
AetFactor = 0.4;
else
AetFactor = 0.6;
}
//*************对蒸散发处理的特殊代码段 -- 计算实际蒸散发**************//
if (g_ModelRunPara.PETMethod == PET_METHOD.PET_REAL)
{
gridwb.m_dAET = g_ClimatePara.petdata.m_gridLayer[row, col] * AetFactor;
g_GridResultOut.m_AET[row, col] = gridwb.m_dAET;
}
else
{
gridwb.CalcAET(g_ModelRunPara, g_ClimatePara, dalb);
if (gridwb.m_dAET > g_ClimatePara.petdata.m_gridLayer[row, col]) // FRZZ MODI 2005-8-18
g_GridResultOut.m_AET[row, col] = g_ClimatePara.petdata.m_gridLayer[row, col] * Math.Exp(-1 * g_ClimatePara.petdata.m_gridLayer[row, col] / gridwb.m_dAET);
if (gridwb.m_dAET < 0)
g_GridResultOut.m_AET[row, col] = g_ClimatePara.petdata.m_gridLayer[row, col] * 0.1;
}
//**********************对蒸散发处理的特殊代码段**********************//
gridwb.CalcNetRain();
g_GridResultOut.m_CIDefict[row, col] = gridwb.m_dCIDeficit;
g_GridResultOut.m_NetPcp[row, col] = gridwb.m_dNetRain;
g_GridResultOut.m_GridWaterYieldType[row, col] = gridwb.CalcRunoffElement(g_GridLayerPara, g_ModelInputPara,g_GridResultOut.m_SoilProfileWater[row, col]);
g_GridResultOut.m_SoilProfileWater[row, col] = gridwb.Cell_SP_Sw; //g_GridLayerPara.SoilTypes[SoilID].SP_Sw;
g_GridResultOut.m_SoilAvgWater[row, col] = g_GridLayerPara.SoilTypes[SoilID].SoilAvgWater(gridwb.Cell_SP_Sw);
g_GridResultOut.m_GridTotalQ[row, col] = gridwb.m_dTotalQ;
g_GridResultOut.m_GridSurfQ[row, col] = gridwb.m_dSurfQ;
g_GridResultOut.m_GridLateralQ[row, col] = gridwb.m_dLateralQ;
g_GridResultOut.m_GridBaseQ[row, col] = gridwb.m_dBaseQ;
if (g_GridResultOut.m_GridSurfQ[row, col] > 1e+10)
MessageBox.Show("hello2");
}
}
}
if (m_iNodeNum == 1 || g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_PURE_LAG)
{
PureLagGridRouting(g_GridResultOut.m_GridSurfQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_SURFQ, curDate, dsnowfactor, waterYearType[curDate.Year]);
PureLagGridRouting(g_GridResultOut.m_GridLateralQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_LATERALQ, curDate, dsnowfactor, waterYearType[curDate.Year]);
PureLagGridRouting(g_GridResultOut.m_GridBaseQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_BASEQ, curDate, dsnowfactor, waterYearType[curDate.Year]);
int dn = DateAndTime.GetDnInYear(curDate.Year, curDate.Month, curDate.Day);
this.MyOutletQ[curDate].m_pOutletDeepBaseQ = DeepBaseQSim(dn, g_ModelInputPara.DeepBaseQ);
this.MyOutletQ[curDate].m_pOutletQ = this.MyOutletQ[curDate].m_pOutletSurfQ * g_ModelInputPara.SurfQLinearFactor
+ this.MyOutletQ[curDate].m_pOutletLatQ + this.MyOutletQ[curDate].m_pOutletBaseQ + this.MyOutletQ[curDate].m_pOutletDeepBaseQ;
}
else
{
PureLagGridRouting_Node(g_GridResultOut.m_GridSurfQ, m_pNodeSurfQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_SURFQ, curDate, totrec, dsnowfactor, waterYearType[curDate.Year]);
PureLagGridRouting_Node(g_GridResultOut.m_GridLateralQ, m_pNodeLatQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_LATERALQ, curDate, totrec, dsnowfactor, waterYearType[curDate.Year]);
PureLagGridRouting_Node(g_GridResultOut.m_GridBaseQ, m_pNodeBaseQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_BASEQ, curDate, totrec, dsnowfactor, waterYearType[curDate.Year]);
double dSurf, dLat, dBase;
dSurf = dLat = dBase = 0;
for (int i = 0; i < g_GridLayerPara.m_subNum; i++)
{
//dSurf += m_pNodeSurfQ[curorder, i] * g_ModelInputPara1.SurfQLinearFactor;
//dLat += m_pNodeLatQ[curorder, i];
//dBase += m_pNodeBaseQ[curorder, i];
}
this.MyOutletQ[curDate].m_pOutletSurfQ = dSurf;
this.MyOutletQ[curDate].m_pOutletLatQ = dLat;
this.MyOutletQ[curDate].m_pOutletBaseQ = dBase;
int dn = DateAndTime.GetDnInYear(curDate.Year, curDate.Month, curDate.Day);
this.MyOutletQ[curDate].m_pOutletDeepBaseQ = DeepBaseQSim(dn, g_ModelInputPara.DeepBaseQ);
if (g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_MUSKINGUM_COMBINE_FIRST)
{
for (int i = 0; i < g_GridLayerPara.m_subNum; i++)
{
m_pNodeOutQ[curDate][i] = m_pNodeSurfQ[curDate][i] * g_ModelInputPara.SurfQLinearFactor
+ m_pNodeLatQ[curDate][i] + m_pNodeBaseQ[curDate][i];
}
MuskingumRiverRouting(24, m_pNodeOutQ, m_pX, m_pK, g_GridLayerPara.m_subNum, curDate);
this.MyOutletQ[curDate].m_pOutletQ = pRiverRoute[g_GridLayerPara.m_subNum - 1].pRoute.OutFlux
+ this.MyOutletQ[curDate].m_pOutletDeepBaseQ;
}
else if (g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_MUSKINGUM_ROUTE_FIRST)
{
;
}
else
{
this.MyOutletQ[curDate].m_pOutletQ = this.MyOutletQ[curDate].m_pOutletSurfQ +
this.MyOutletQ[curDate].m_pOutletLatQ +
this.MyOutletQ[curDate].m_pOutletBaseQ +
this.MyOutletQ[curDate].m_pOutletDeepBaseQ;
}
}
MidGridResultOut(curDate, true);
}
RiverOutletQ(true);
m_iTotalRec = totrec;
}
/// <summary>
/// 通过标准日期类型转化为字符串(Julian历格式)
/// the result maybe like 1998001,1998002,...1998365,
/// </summary>
/// <param name="date">日期</param>
/// <returns>格式字符串</returns>
public string SetLongDateByDn(DateTime date)
{
int dn = DateAndTime.GetDnInYear(date.Year, date.Month, date.Day);
return dn.ToString();
}
/// <summary>
/// 根据年份和Julian历日期，返回日期字符串
/// </summary>
/// <param name="year">年份（4位）</param>
/// <param name="dn">dn -- Julian历日期</param>
/// <returns>日期字符串（格式：1998-1-1）</returns>
public DateTime DnToDate(int year, int dn)
{
int date = 0;
int[] DnNum = new int[12];
for (int i = 1; i <= 12; i++)
DnNum[i - 1] = DateAndTime.GetEachMonthDn(year, i);
DateTime datetime = new DateTime();
if (dn <= DnNum[0])
{
datetime = new DateTime(year, 1, dn);
}
else
{
for (int i = 1; i < 12; i++) //other month
{
if (dn <= DnNum[i] && dn > DnNum[i - 1])
{
date = dn - DnNum[i - 1];
datetime = new DateTime(year, i + 1, date);
}
}
}
return datetime;
}
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ 暴雨～径流过程模拟函数定义 +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/// <summary>
/// 单独计算一个栅格的水量平衡
/// </summary>
/// <param name="row">行号</param>
/// <param name="col">列号</param>
/// <param name="mon">月份</param>
/// <param name="dhr">当前时段的小时数</param>
/// <param name="curorder">当前执行序号</param>
/// <param name="dn">年内的日序号</param>
public void CalcCell_Horton(int row, int col, double dhr, DateTime curorder, int dn)
{
//如果该网格点不需要计算，则跳过
if (this.IfGridNeedCal(row, col) == false)
return;
int VegID = (int)g_GridLayerPara.g_VegLayer[row, col];
double dalb = this.GetVegAlbedo(VegID, curorder.Month, 1);
//如果栅格点上有河流级别，执行下面
if (River.m_pStrahlerOrd[row, col] != 0)
{
//降水量减去蒸散发，取得净降水量
if (g_ClimatePara.pcpdata.m_gridLayer[row, col] < 0 || g_ClimatePara.petdata.m_gridLayer[row, col] < 0)
throw new Exception("降雨或蒸发量数据不能为负值");
g_GridResultOut.m_GridSurfQ[row, col] = g_ClimatePara.pcpdata.m_gridLayer[row, col] - g_ClimatePara.petdata.m_gridLayer[row, col];
//如果得到的差为负值，则地表径流设为零
if (g_GridResultOut.m_GridSurfQ[row, col] < 0)
{
g_GridResultOut.m_GridSurfQ[row, col] = 0;
}
//栅格壤中流和栅格地下径流都设为0(这是因为河道内可以说无壤中流和潜流)
g_GridResultOut.m_GridLateralQ[row, col] = 0;
g_GridResultOut.m_GridBaseQ[row, col] = 0;
}
else//如果该栅格无河流级别
{
//计算降雨强度(以每小时计算(mm/h))
double dintensity = g_ClimatePara.pcpdata.m_gridLayer[row, col] / dhr;
//计算Horton超渗产流
HortonInfil hortonInfil= new HortonInfil();
hortonInfil.SetGridPara(g_GridLayerPara, row, col, g_GridResultOut.m_SoilProfileWater[row,col]);
//计算本栅格土壤水实际下渗量
m_drateinf[row, col] = hortonInfil.HortonExcessRunoff();
//下面采用水平衡类进行计算
GridWaterBalance gridwb = new GridWaterBalance();
//设置栅格上用于计算水平衡的参数
gridwb.SetGridPara(row, col, dintensity,m_drateinf.Values[row,col], curorder, dhr);
double aetfactor = 1.0;
int SoilID = (int)g_GridLayerPara.g_SoilLayer[row, col];
if (g_GridResultOut.m_SoilProfileWater[row,col] / g_GridLayerPara.SoilTypes[SoilID].SP_Fc > 0.8)
{
if (g_ClimatePara.pcpdata.m_gridLayer[row, col] > 0)
aetfactor = 0.75;
else
aetfactor = 0.95;
}
else
{
if (g_ClimatePara.pcpdata.m_gridLayer[row, col] > 0)
aetfactor = 0.45;
else
aetfactor = 0.65;
}
//*************对蒸散发处理的特殊代码段 -- 计算实际蒸散发**************//
gridwb.m_dAET = g_ClimatePara.petdata.m_gridLayer[row, col] * aetfactor;
g_GridResultOut.m_AET[row, col] = gridwb.m_dAET;
//计算净雨量
gridwb.CalcNetRain();
g_GridResultOut.m_NetPcp[row, col] = gridwb.m_dNetRain;
g_GridResultOut.m_GridWaterYieldType[row, col] = gridwb.CalcRunoffElement(g_GridLayerPara, g_ModelInputPara,g_GridResultOut.m_SoilProfileWater[row, col]);
g_GridResultOut.m_SoilProfileWater[row, col] = gridwb.Cell_SP_Sw;
g_GridResultOut.m_SoilAvgWater[row, col] = g_GridLayerPara.SoilTypes[SoilID].SoilAvgWater(gridwb.Cell_SP_Sw);
g_GridResultOut.m_GridTotalQ[row, col] = gridwb.m_dTotalQ;
g_GridResultOut.m_GridSurfQ[row, col] = gridwb.m_dSurfQ;
g_GridResultOut.m_GridLateralQ[row, col] = gridwb.m_dLateralQ;
g_GridResultOut.m_GridBaseQ[row, col] = gridwb.m_dBaseQ;
//Console.WriteLine(gridwb.m_dLateralQ);
}
}
public void CalcCell_GreenAmpt(int row, int col, double dhr, DateTime curorder, int dn)
{
int VegID = (int)g_GridLayerPara.g_VegLayer[row, col];
double dalb = this.GetVegAlbedo(VegID, curorder.Month, 1);
string key = g_ClimatePara.petdata.saStationName[0];
//**********************对蒸散发处理的特殊代码段**********************//
//如果当前栅格的strahler分级为河流
if (River.m_pStrahlerOrd[row, col] != 0)
{
g_GridResultOut.m_GridSurfQ[row, col] = g_ClimatePara.pcpdata.m_gridLayer[row, col] - g_ClimatePara.petdata.m_gridLayer[row, col];
if (g_GridResultOut.m_GridSurfQ[row, col] < 0.0)
g_GridResultOut.m_GridSurfQ[row, col] = 0.0;
g_GridResultOut.m_GridLateralQ[row, col] = 0.0;
g_GridResultOut.m_GridBaseQ[row, col] = 0.0;
}
else //strahler分级不是河流
{
double dintensity = g_ClimatePara.pcpdata.m_gridLayer[row, col] / dhr;//降雨强度
double dCp = g_GridLayerPara.VegTypes[VegID].CoverDeg[curorder.Month - 1];
GreenAmptInfil GreenAmpt = new GreenAmptInfil(g_GridLayerPara);
GreenAmpt.SetGridPara(row, col, 1.0, m_drateinf[row, col], m_dcumr[row, col], m_drintns[row, col], m_dcuminf[row, col], m_dexcum[row, col], g_ClimatePara.pcpdata.m_gridLayer[row, col], dCp);
GreenAmpt.GreenAmptExcessRunoff(g_GridResultOut.m_SoilProfileWater.Values[row, col]);
m_drateinf[row, col] = GreenAmpt.m_dPreRateinf;
m_dcumr[row, col] = GreenAmpt.m_dPrecumr;
m_drintns[row, col] = GreenAmpt.m_dPrerintns;
m_dcuminf[row, col] = GreenAmpt.m_dPrecuminf;
m_dexcum[row, col] = GreenAmpt.m_dPreexcum;
int SoilID = (int)g_GridLayerPara.g_SoilLayer[row, col];
GridWaterBalance gridwb = new GridWaterBalance();
gridwb.SetGridPara(row, col, dintensity, m_drateinf.Values[row, col], curorder, dhr);
//**********************对蒸散发处理的特殊代码段**********************//
if (g_ClimatePara.PetInfo.StationFileName == "" && g_ClimatePara.PetInfo.DataFileName != "")
{
key = g_ClimatePara.petdata.saStationName[0];
gridwb.m_dAET = g_ClimatePara.petdata.pMeteoData[key][curorder];
g_GridResultOut.m_AET[row, col] = gridwb.m_dAET;
}
//**********************对蒸散发处理的特殊代码段**********************//
gridwb.CalcNetRain();
g_GridResultOut.m_NetPcp[row, col] = gridwb.m_dNetRain;
g_GridResultOut.m_GridWaterYieldType[row, col] = gridwb.CalcRunoffElement(g_GridLayerPara, g_ModelInputPara, g_GridResultOut.m_SoilProfileWater[row, col]);
g_GridResultOut.m_SoilProfileWater[row, col] = gridwb.Cell_SP_Sw;
g_GridResultOut.m_SoilAvgWater[row, col] = g_GridLayerPara.SoilTypes[SoilID].SoilAvgWater(gridwb.Cell_SP_Sw);
g_GridResultOut.m_GridTotalQ[row, col] = gridwb.m_dTotalQ;
g_GridResultOut.m_GridSurfQ[row, col] = gridwb.m_dSurfQ;
g_GridResultOut.m_GridLateralQ[row, col] = gridwb.m_dLateralQ;
g_GridResultOut.m_GridBaseQ[row, col] = gridwb.m_dBaseQ;
}
}
/// <summary>
/// 进行汇流计算
/// </summary>
/// <param name="curorder">当前的执行时段序号</param>
/// <param name="dhr">当前时段的小时数</param>
/// <param name="totrec">总时段数</param>
/// <param name="dn">年内的日序</param>
public void CalcRiverRouting_Horton(DateTime curorder, double dhr, int dn)
{
//下面进行坡面漫流汇流
if (g_ModelRunPara.SurfRouteMethod == SurfRouteMehtod.ROUTE_MUSK_CONGE)
{
this.MuskCungeGridRouting(dhr, g_GridResultOut.m_GridSurfQ, g_ModelInputPara.SMCTimeWeight, g_ModelInputPara.SMCGridRTravelTime);
for (int k = 0; k < RoutePara.Length; k++)
{
int row = RoutePara[k].pRow;
int col = RoutePara[k].pCol;
g_GridResultOut.m_GridRoutingQ[row, col] = pSurfQ[k].pRoute.OutFlux;
}
//把路径上最后一个像素的出流量作为总出流量
if (this.MyOutletQ.ContainsKey(curorder) == false)
this.MyOutletQ[curorder] = new OutletQ();
this.MyOutletQ[curorder].m_pOutletSurfQ = pSurfQ[RoutePara.Length - 1].pRoute.OutFlux;
}
else
{
this.PureLagGridRouting(g_GridResultOut.m_GridSurfQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_SURFQ, curorder);
}
//下面进行壤中流汇流
if (g_ModelRunPara.LatRouteMethod == SurfRouteMehtod.ROUTE_MUSK_CONGE)
{
MuskCungeGridRouting(dhr, g_GridResultOut.m_GridLateralQ, g_ModelInputPara.LMCTimeWeight, g_ModelInputPara.LMCGridRTravelTime);
for (int k = 0; k < RoutePara.Length; k++)
{
g_GridResultOut.m_GridRoutingQ.Values[RoutePara[k].pRow, RoutePara[k].pCol] += pLatQ[k].pRoute.OutFlux;
}
if (this.MyOutletQ.ContainsKey(curorder) == false)
this.MyOutletQ[curorder] = new OutletQ();
this.MyOutletQ[curorder].m_pOutletLatQ = pLatQ[RoutePara.Length - 1].pRoute.OutFlux;
}
else
{
PureLagGridRouting(g_GridResultOut.m_GridLateralQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_LATERALQ, curorder);
}
//下面进行地下径流的汇流
if (g_ModelRunPara.BaseRouteMethod == SurfRouteMehtod.ROUTE_MUSK_CONGE)
{
this.MuskCungeGridRouting(dhr, g_GridResultOut.m_GridBaseQ, g_ModelInputPara.BMCTimeWeight, g_ModelInputPara.BMCGridRTravelTime);
for (int k = 0; k < RoutePara.Length; k++)
{
g_GridResultOut.m_GridRoutingQ.Values[RoutePara[k].pRow, RoutePara[k].pCol] += pBaseQ[k].pRoute.OutFlux;
}
if (this.MyOutletQ.ContainsKey(curorder) == false)
this.MyOutletQ[curorder] = new OutletQ();
this.MyOutletQ[curorder].m_pOutletBaseQ = pBaseQ[RoutePara.Length - 1].pRoute.OutFlux;
}
else
{
PureLagGridRouting(g_GridResultOut.m_GridBaseQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_BASEQ, curorder);
}
//下面计算深层地下径流的汇流
this.MyOutletQ[curorder].m_pOutletDeepBaseQ = DeepBaseQSim(dn, g_ModelInputPara.DeepBaseQ);
this.MyOutletQ[curorder].m_pOutletQ = this.MyOutletQ[curorder].m_pOutletSurfQ * g_ModelInputPara.SurfQLinearFactor
+ this.MyOutletQ[curorder].m_pOutletLatQ + this.MyOutletQ[curorder].m_pOutletBaseQ + this.MyOutletQ[curorder].m_pOutletDeepBaseQ;
}
/// <summary>
/// 格林－安普特下渗曲线计算暴雨过程主函数
/// </summary>
/// <returns></returns>
public void StormRunoffSim_GreenAmpt_Liu()
{
int rowCount = g_GridLayerPara.g_DemLayer.rowCount;
int colCount = g_GridLayerPara.g_DemLayer.colCount;
this.CellOutSurfFlux = new double[rowCount, colCount];
this.CellOutLatFlux = new double[rowCount, colCount];
this.CellOutBaseFlux = new double[rowCount, colCount];
GridLayerInit();
GridLayerInit_GreenAmpt();
this.MyOutletQ = new SortedList<DateTime, OutletQ>();
int n = River.CalcuOrderList.Count;
CurrentSurfQ = new RouteFlux[n];
PreSurfQ = new RouteFlux[n];
for (int i = 0; i < n; i++)
{
CurrentSurfQ[i] = new RouteFlux();
PreSurfQ[i] = new RouteFlux();
}
StreamWriter sw = new StreamWriter(HydroSimulate.path + "\\Output\\OutletQ_Liu.txt");
sw.WriteLine("Order\tTotalQ\tSurfQ\tLatQ\tBaseQ\tDeepBaseQ");
//对各降雨时段进行循环
int count = g_ClimatePara.pcpdata.saDate.Count;
for (int i = 0; i < count; i++)
{
//功能：给定年、月、日，计算该日的Julian历日
DateTime curorder = g_ClimatePara.pcpdata.saDate[i];
//对i时段的降雨量进行内插
MeteoDataInterp(curorder); //气象数据内插
//间隔的小时数
double dhr = 1.0;
int dn = DateAndTime.GetDnInYear(curorder.Year, curorder.Month, curorder.Day);
for (int row = 0; row < rowCount; row++)
{
for (int col = 0; col < colCount; col++)
{
if (!IfGridNeedCal(row, col))
continue;
this.CalcCell_GreenAmpt(row, col, 1.0, curorder, dn);
}
}
this.CalcRiverRouting_Liu(curorder, 1, dn);
sw.Write(curorder + "\t");
sw.Write(this.MyOutletQ[curorder].m_pOutletQ.ToString("0.0000") + "\t");
sw.Write(this.MyOutletQ[curorder].m_pOutletSurfQ.ToString("0.0000") + "\t");
sw.Write(this.MyOutletQ[curorder].m_pOutletLatQ.ToString("0.0000") + "\t");
sw.Write(this.MyOutletQ[curorder].m_pOutletBaseQ.ToString("0.0000") + "\t");
sw.Write(this.MyOutletQ[curorder].m_pOutletDeepBaseQ.ToString("0.0000") + "\t");
sw.WriteLine();
sw.Flush();
MidGridResultOut(g_ClimatePara.pcpdata.saDate[i], false);
}
//输出断面流量结果
this.RiverOutletQ(true);
// true;
}
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ 流域汇流模拟函数定义 +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/// <summary>
/// 读入栅格汇流最优次序参数文件
/// </summary>
/// <returns></returns>
public void ReadInRoutingPara()
{
m_OutRow = m_OutCol = 0;
string strFileStub = g_GridLayerPara.DemFileName.Substring(0, g_GridLayerPara.DemFileName.Length - 4);
string strFileName = strFileStub + "_gud.txt";
if (!File.Exists(strFileName))
throw new Exception("栅格汇流文件" + strFileName + "不存在");
List<string> saIn = MatrixFuncs.FileRead(strFileName);
int LineCount = saIn.Count; //行数
RoutePara = new MuskCungeRoutePara[LineCount];
pRouteQ = new double[LineCount];
pSurfQ = new MuskRouteFlux[LineCount];
pLatQ = new MuskRouteFlux[LineCount];
pBaseQ = new MuskRouteFlux[LineCount];
ProgressBar bar = new ProgressBar();
bar.Text = "正在加载全栅格汇流参数文件";
bar.Show();
for (int row = 0; row < LineCount; row++)
{
string strLine = saIn[row];
string[] saLine = strLine.Split(new char[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
RoutePara[row] = new MuskCungeRoutePara();
RoutePara[row].pGridNum = Convert.ToInt32(saLine[0]);
RoutePara[row].pGridRouteOrd = Convert.ToInt32(saLine[1]);
for (int i = 0; i < 8; i++)
{
RoutePara[row].pInGrid[i] = Convert.ToInt32(saLine[i + 2]);
}
RoutePara[row].pOutGrid = Convert.ToInt32(saLine[10]);
RoutePara[row].pGridSlope = Convert.ToDouble(saLine[11]);
RoutePara[row].pGridRLength = Convert.ToDouble(saLine[12]);
RoutePara[row].pGridRiverOrd = Convert.ToInt32(saLine[13]);
RoutePara[row].pRow = Convert.ToInt32(saLine[14]);
RoutePara[row].pCol = Convert.ToInt32(saLine[15]);
pSurfQ[row] = new MuskRouteFlux();
pSurfQ[row].pRoute = new RouteFlux();
pSurfQ[row].pPreRoute = new RouteFlux();
pSurfQ[row].pPreRoute.bCal = true;
pLatQ[row] = new MuskRouteFlux();
pLatQ[row].pRoute = new RouteFlux();
pLatQ[row].pPreRoute = new RouteFlux();
pLatQ[row].pPreRoute.bCal = true;
pBaseQ[row] = new MuskRouteFlux();
pBaseQ[row].pRoute = new RouteFlux();
pBaseQ[row].pPreRoute = new RouteFlux();
pBaseQ[row].pPreRoute.bCal = true;
bar.progressBar1.Value = (int)(row * 1.0 / LineCount * 100);
}
bar.Close();
m_OutRow = RoutePara[LineCount - 1].pRow;
m_OutCol = RoutePara[LineCount - 1].pCol;
}
/// <summary>
/// 读入栅格汇流相关参数，包括全流域边界、子流域边界、河道等级矩阵、地表径流汇流参数矩阵
/// 壤中径流汇流参数矩阵、地下径流汇流参数矩阵、地下径流参数矩阵
/// </summary>
/// <returns></returns>
public void ReadInRoutingLayerData()
{
string strFileStub = g_GridLayerPara.DemFileName.Substring(0, g_GridLayerPara.DemFileName.Length - 4);
string strFileName = strFileStub + "_wby.asc";
g_GridLayerPara.g_BasinBoundary = new GridLayer();
g_GridLayerPara.g_BasinBoundary.ReadASC(strFileName);
strFileName = strFileStub + "_sws.asc";
g_GridLayerPara.g_SubWaterShed = new GridLayer();
g_GridLayerPara.g_SubWaterShed.ReadASC(strFileName);
g_GridLayerPara.g_SubWaterShed.CalMaxMinValue();
g_GridLayerPara.m_subNum = (int)g_GridLayerPara.g_SubWaterShed.MaxValue;
strFileName = strFileStub + "_sor.asc";
g_GridLayerPara.g_StrahlerRivNet = new GridLayer();
g_GridLayerPara.g_StrahlerRivNet.ReadASC(strFileName);
g_GridLayerPara.g_StrahlerRivNet.CalMaxMinValue();
g_GridLayerPara.m_MaxStrahlerOrd = (int)g_GridLayerPara.g_StrahlerRivNet.MaxValue;
if (g_ModelRunPara.SurfRouteMethod == SurfRouteMehtod.ROUTE_PURE_LAG)
{
strFileName = strFileStub + "_gst.asc";
g_GridLayerPara.g_RouteSurfQTime.ReadASC(strFileName);
}
if (g_ModelRunPara.LatRouteMethod == SurfRouteMehtod.ROUTE_PURE_LAG)
{
strFileName = strFileStub + "_glt.asc";
g_GridLayerPara.g_RouteLatQTime = new GridLayer();
g_GridLayerPara.g_RouteLatQTime.ReadASC(strFileName);
}
if (g_ModelRunPara.BaseRouteMethod == SurfRouteMehtod.ROUTE_PURE_LAG)
{
strFileName = strFileStub + "_gbt.asc";
g_GridLayerPara.g_RouteBaseQTime = new GridLayer();
g_GridLayerPara.g_RouteBaseQTime.ReadASC(strFileName);
}
}
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ 径流过程模拟私有函数定义 +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
public void GridLayerInit()
{
GridLayer DemLayer=g_GridLayerPara.g_DemLayer;
g_GridResultOut.m_GridTotalQ = DemLayer.AttributesCopy();
g_GridResultOut.m_GridSurfQ = DemLayer.AttributesCopy();
g_GridResultOut.m_GridLateralQ = DemLayer.AttributesCopy();
g_GridResultOut.m_GridBaseQ = DemLayer.AttributesCopy();
g_GridResultOut.m_GridWaterYieldType = DemLayer.AttributesCopy();
g_GridResultOut.m_SoilAvgWater = DemLayer.AttributesCopy();
g_GridResultOut.m_GridRoutingQ = DemLayer.AttributesCopy();
g_GridResultOut.m_NetPcp = DemLayer.AttributesCopy();
g_GridResultOut.m_AET = DemLayer.AttributesCopy();
g_GridResultOut.m_AI = DemLayer.AttributesCopy();
g_GridResultOut.m_CI = DemLayer.AttributesCopy();
g_GridResultOut.m_CIDefict = DemLayer.AttributesCopy();
g_GridResultOut.m_NetPcp = DemLayer.AttributesCopy();
g_GridResultOut.m_GridWaterYieldType = DemLayer.AttributesCopy();
m_drateinf = DemLayer.AttributesCopy(0);
g_GridResultOut.m_SoilProfileWater = DemLayer.AttributesCopy();
for (int row = 0; row < DemLayer.rowCount; row++)
{
for (int col = 0; col < DemLayer.colCount; col++)
{
if (River.m_pGridRouteLength[row, col] > 0)
{
int SoilID =(int)g_GridLayerPara.g_SoilLayer[row, col];
SoilInfo soilInfo = g_GridLayerPara.SoilTypes[SoilID];
g_GridResultOut.m_SoilProfileWater[row, col] = soilInfo.Initial_SP_Sw;
}
}
}
}
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ +
+ 气象数据空间插值 +
+ +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/// <summary>
/// 降雨量空间插值（有时段性问题，所以和其它气象数据空间插值分开）
/// </summary>
/// <param name="ord"></param>
public void PcpDataInterp(DateTime timeOrd)
{
g_ClimatePara.pcpdata.SpatialGridOut(timeOrd);
}
public void MeteoDataInterp(DateTime timeOrd)
{
GridLayer pcp= g_ClimatePara.pcpdata.SpatialGridOut(timeOrd);
//pcp.WriteASC(path + "\\Test\\pcp.txt");
//g_ClimatePara.tmpmxdata.SpatialGridOut(timeOrd, 2, false);
//g_ClimatePara.tmpmndata.SpatialGridOut(timeOrd, 2, false);
//g_ClimatePara.slrdata.SpatialGridOut(timeOrd, 2, false);
//g_ClimatePara.hmddata.SpatialGridOut(timeOrd, 2, false);
//g_ClimatePara.WindData.SpatialGridOut(timeOrd, 2, false);
g_ClimatePara.petdata.SpatialGridOut(timeOrd);
//g_ClimatePara.tmpmeandata.SpatialGridOut(timeOrd, 2, false);
//else
//{
// if (g_ClimatePara.bReadTmpMX == true && g_ClimatePara.bReadTmpMN == true)
// {
// g_ClimatePara.tmpmeandata.m_gridLayer.extent = g_ClimatePara.tmpmndata.m_gridLayer.extent.Copy();
// g_ClimatePara.tmpmeandata.m_gridLayer.resolution = g_ClimatePara.tmpmndata.m_gridLayer.resolution;
// g_ClimatePara.tmpmeandata.m_gridLayer.initGridStruct();
// for (int j = 0; j < g_ClimatePara.tmpmeandata.m_gridLayer.rowCount; j++)
// {
// for (int k = 0; k < g_ClimatePara.tmpmeandata.m_gridLayer.colCount; k++)
// {
// double dvalue;
// dvalue = (g_ClimatePara.tmpmxdata.m_gridLayer.Values[j, k] + g_ClimatePara.tmpmndata.m_gridLayer.Values[j, k]) / 2.0;
// g_ClimatePara.tmpmeandata.m_gridLayer.Values[j, k] = dvalue;
// }
// }
// }
//}
}
/// <summary>
/// 网格中的中间结果输出
/// </summary>
/// <param name="yr"></param>
/// <param name="dn"></param>
/// <param name="ord"></param>
/// <param name="bDate"></param>
public void MidGridResultOut(DateTime dateTime, bool bDate)
{
DateTime strCurDate = dateTime;
string TimeString = strCurDate.Year + "_" + strCurDate.Month + "_" + strCurDate.Day + "_" + strCurDate.Hour;
string strtmp;
if (HydroSimulate.g_GridResultOut.iPcp)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\pcp_" + TimeString + ".asc";
HydroSimulate.g_ClimatePara.pcpdata.m_gridLayer.WriteASC(strtmp);
}
if (g_GridResultOut.iTempMax)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\TempMax_" + TimeString + ".asc";
HydroSimulate.g_ClimatePara.tmpmxdata.m_gridLayer.WriteASC(strtmp);
}
if (g_GridResultOut.iTempMin)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\TempMin_" + TimeString + ".asc";
HydroSimulate.g_ClimatePara.tmpmndata.m_gridLayer.WriteASC(strtmp);
}
if (g_GridResultOut.iTempMean)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\TempMean_" + TimeString + ".asc";
HydroSimulate.g_ClimatePara.tmpmeandata.m_gridLayer.WriteASC(strtmp);
}
if (g_GridResultOut.iSlr)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\Solar_" + TimeString + ".asc";
HydroSimulate.g_ClimatePara.slrdata.m_gridLayer.WriteASC(strtmp);
}
if (g_GridResultOut.iWnd)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\Wnd_" + TimeString + ".asc";
HydroSimulate.g_ClimatePara.WindData.m_gridLayer.WriteASC(strtmp);
}
if (g_GridResultOut.iHmd)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\Hmd_" + TimeString + ".asc";
HydroSimulate.g_ClimatePara.hmddata.m_gridLayer.WriteASC(strtmp);
}
if (g_GridResultOut.iPET)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\PET_" + TimeString + ".asc";
HydroSimulate.g_ClimatePara.petdata.m_gridLayer.WriteASC(strtmp);
}
if (g_GridResultOut.iAET)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\AET_" + TimeString + ".asc";
HydroSimulate.g_GridResultOut.m_AET.WriteASC(strtmp);
}
if (g_GridResultOut.iCI)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\CI_" + TimeString + ".asc";
HydroSimulate.g_GridResultOut.m_CI.WriteASC(strtmp);
}
if (g_GridResultOut.iSnowWater)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\Snow_" + TimeString + ".asc";
HydroSimulate.g_GridResultOut.m_SnowWater.WriteASC(strtmp);
}
if (g_GridResultOut.iAI)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\AI_" + TimeString + ".asc";
HydroSimulate.g_GridResultOut.m_AI.WriteASC(strtmp);
}
if (g_GridResultOut.iRouteOut)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\GridRouteQ_" + TimeString + ".asc";
HydroSimulate.g_GridResultOut.m_GridRoutingQ.WriteASC(strtmp);
}
if (g_GridResultOut.iSurfQ)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\SurfQ_" + TimeString + ".asc";
HydroSimulate.g_GridResultOut.m_GridSurfQ.WriteASC(strtmp);
}
if (g_GridResultOut.iWaterYieldType)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\WYT_" + TimeString + ".asc";
HydroSimulate.g_GridResultOut.m_GridWaterYieldType.WriteASC(strtmp);
}
if (g_GridResultOut.iLatQ)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\LatQ_" + TimeString + ".asc";
HydroSimulate.g_GridResultOut.m_GridLateralQ.WriteASC(strtmp);
}
if (g_GridResultOut.iBaseQ)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\BaseQ_" + TimeString + ".asc";
HydroSimulate.g_GridResultOut.m_GridBaseQ.WriteASC(strtmp);
}
if (g_GridResultOut.iInfilRate)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\InfilRate_" + TimeString + ".asc";
HydroSimulate.m_drateinf.WriteASC(strtmp);
}
if (g_GridResultOut.iProfileSoilWater)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\SoilProfileWater_" + TimeString + ".asc";
HydroSimulate.g_GridResultOut.m_SoilProfileWater.WriteASC(strtmp);
}
if (g_GridResultOut.iAvgSoilWater)
{
strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\SoilAvgWater_" + TimeString + ".asc";
HydroSimulate.g_GridResultOut.m_SoilAvgWater.WriteASC(strtmp);
}
}
/// <summary>
/// 输出所有流量值
/// </summary>
/// <param name="recnum">时段</param>
/// <param name="bDate"></param>
public void RiverOutletQ(bool bDate)
{
if(!Directory.Exists(HydroSimulate.path+"\\Output"))
Directory.CreateDirectory(HydroSimulate.path+"\\Output");
StreamWriter sw = new StreamWriter(HydroSimulate.path + "\\Output\\OutQ.txt");
sw.WriteLine("Order\tTotalQ\tSurfQ\tLatQ\tBaseQ\tDeepBaseQ");
if (bDate)
{
foreach (DateTime timeKey in this.MyOutletQ.Keys)
{
sw.Write(timeKey + "\t");
sw.Write(this.MyOutletQ[timeKey].m_pOutletQ.ToString("0.0000") + "\t");
sw.Write(this.MyOutletQ[timeKey].m_pOutletSurfQ.ToString("0.0000") + "\t");
sw.Write(this.MyOutletQ[timeKey].m_pOutletLatQ.ToString("0.0000") + "\t");
sw.Write(this.MyOutletQ[timeKey].m_pOutletBaseQ.ToString("0.0000") + "\t");
sw.Write(this.MyOutletQ[timeKey].m_pOutletDeepBaseQ.ToString("0.0000") + "\t");
sw.WriteLine();
}
}
sw.Close();
}
/// <summary>
/// 判断当前网格是否需要计算
/// </summary>
/// <param name="row"></param>
/// <param name="col"></param>
/// <returns>true-需要计算，false-不需要计算</returns>
public bool IfGridNeedCal(int row, int col)
{
//如果该网格点是流域边界，则不需要计算
//if (g_GridLayerPara.g_BasinBoundary[row, col] != 1)
// return false;
if (this.River.m_pGridRouteLength[row, col] == 0)
return false;
//如果网格点上的植被数据为空
if (g_GridLayerPara.g_VegLayer[row, col] == g_GridLayerPara.g_VegLayer.nodata)
return false;
//如果网格点上的土壤数据为空
if (g_GridLayerPara.g_SoilLayer[row, col] == g_GridLayerPara.g_SoilLayer.nodata)
return false;
//如果网格点上的DEM数据为空
if (g_GridLayerPara.g_DemLayer[row, col] == g_GridLayerPara.g_DemLayer.nodata)
return false;
//不符合以上条件的则需要计算
return true;
}
/// <summary>
/// 度－日因子模型计算日融雪量
/// </summary>
/// <param name="dtav"></param>
/// <param name="dtThresh"></param>
/// <param name="ddf"></param>
/// <param name="dtlen"></param>
/// <returns></returns>
public double DDFSnowMelt(double dtav, double dtThresh, double ddf, double dtlen)
{
return ddf * (dtav - dtThresh) * dtlen;
}
public double DeepBaseQSim(int dn, double dmin)
{
double dret = 0;
if (dmin <= 0)
return dret;
else
{
dret = (dmin + 3 * dmin) / 2 + (3 * dmin - dmin) / 2 * Math.Sin(2 * Math.PI * (dn - 82) / 365);
}
return dret;
}
//public void GridLayerInit_LongTerm()
//{
// m_drateinf.initGridStruct();
// m_drintns.initGridStruct();
// m_dcumr.initGridStruct();
// m_dcuminf.initGridStruct();
// m_dexcum.initGridStruct();
// g_GridResultOut.m_AI.initGridStruct();
// g_GridResultOut.m_CI.initGridStruct();
// g_GridResultOut.m_SnowWater.initGridStruct();
// g_GridResultOut.m_CIDefict.initGridStruct();
// g_GridResultOut.m_AET.initGridStruct();
// g_ClimatePara.petdata.m_gridLayer.initGridStruct();
// g_GridResultOut.m_NetPcp.initGridStruct();
//}
/// <summary>
/// 按年份取得该年的模拟时段
/// </summary>
/// <param name="yr">指定年份</param>
/// <param name="startyr">起始年</param>
/// <param name="endyr">结束年</param>
/// <param name="startday">起始天</param>
/// <param name="endday">结束天</param>
/// <param name="ibeg">存储起始日</param>
/// <param name="iend">存储结束日</param>
public void GetBEDateInYear(int yr, int startyr, int endyr, int startday, int endday, int ibeg, int iend)
{
bool bLeap = DateAndTime.CheckLeapYear(yr);
if (bLeap)
iend = 366;
else
iend = 365;
if (yr == startyr)
ibeg = startday;
else
ibeg = 1;
if (yr == endyr)
iend = endday;
}
public void GridLayerInit_GreenAmpt()
{
m_drintns = g_GridLayerPara.g_DemLayer.AttributesCopy(0);
m_dcumr = g_GridLayerPara.g_DemLayer.AttributesCopy(0);
m_dcuminf = g_GridLayerPara.g_DemLayer.AttributesCopy(0);
m_dexcum = g_GridLayerPara.g_DemLayer.AttributesCopy(0);
g_GridResultOut.m_AI = g_GridLayerPara.g_DemLayer.AttributesCopy();
g_GridResultOut.m_CI = g_GridLayerPara.g_DemLayer.AttributesCopy();
g_GridResultOut.m_CIDefict = g_GridLayerPara.g_DemLayer.AttributesCopy();
g_GridResultOut.m_AET = g_GridLayerPara.g_DemLayer.AttributesCopy();
g_ClimatePara.petdata.m_gridLayer = g_GridLayerPara.g_DemLayer.AttributesCopy();
}
/// <summary>
/// 马斯京根-康吉法汇流计算(先合后演)
/// </summary>
/// <param name="dT">时间的离散值</param>
/// <param name="GridQ"></param>
/// <param name="pRoute"></param>
/// <param name="pPreRoute"></param>
/// <param name="x"></param>
/// <param name="k"></param>
/// <returns></returns>
public void MuskCungeGridRouting(double dT, GridLayer GridQ, double x, double k)
{
ProgressBar bar = new ProgressBar();
bar.Text = "正在计算汇流";
bar.Show();
//按照汇流次序对逐个栅格进行汇流
for (int i = 0; i < RoutePara.Length; i++)
{
//先统计当前栅格(即第i个)的总入流量
double SumQ = 0;
for (int j = 0; j < 8; j++)
{
if (RoutePara[i].pInGrid[j] == 0) //第i个栅格在j方向上无入流，则绕过
continue;
if (pSurfQ[RoutePara[i].pInGrid[j] - 1].pRoute.bCal == true)
SumQ += pSurfQ[RoutePara[i].pInGrid[j] - 1].pRoute.OutFlux;
else
throw new Exception("routing grid uncalculated");
}
pSurfQ[i].pRoute.InFlux = SumQ;
//再使用Cunge法计算出流量
pSurfQ[i].pRoute.OutFlux = GridRoutingOut(i, RoutePara[i].pGridSlope, RoutePara[i].pGridRLength, dT,
pSurfQ[i].pPreRoute.InFlux, pSurfQ[i].pRoute.InFlux, pSurfQ[i].pPreRoute.OutFlux,
GridQ.Values[RoutePara[i].pRow, RoutePara[i].pCol] * Math.Pow(GridQ.resolution, 2) / (1000 * dT * 3600),//栅格面入流
x, k);
pSurfQ[i].pRoute.bCal = true;
bar.progressBar1.Value = (int)(i * 100.0 / RoutePara.Length);
}
bar.Close();
//标识新的入流和出流关系
for (int i = 0; i < RoutePara.Length; i++)
{
pSurfQ[i].pPreRoute.InFlux = pSurfQ[i].pRoute.InFlux;
pSurfQ[i].pPreRoute.OutFlux = pSurfQ[i].pRoute.OutFlux;
pSurfQ[i].pRoute.bCal = false;
}
}
/// <summary>
/// 计算序号为ord的栅格的汇流结果，采用Cunge法
/// </summary>
/// <param name="ord">栅格的序号</param>
/// <param name="gridslp">网格坡度</param>
/// <param name="deltaX">河长的离散值</param>
/// <param name="deltaT">时间的离散值</param>
/// <param name="Q11">差分值</param>
/// <param name="Q12">差分值</param>
/// <param name="Q21">差分值</param>
/// <param name="dGridQ">当前时段的侧向流量</param>
/// <param name="x">河段</param>
/// <param name="k">时段</param>
/// <returns>流出流量</returns>
public double GridRoutingOut(int ord, double gridslp, double deltaX, double deltaT, double Q11, double Q12, double Q21, double dGridQ, double x, double k)
{
double dVflow = 0;
RiverType iRiverType;
double dB = 0;//用于存放河道或山的宽度
double dUnitQ = 0;
double Q22 = 0;
//如果是河道，计算河道流速
if (RoutePara[ord].pGridRiverOrd > 0)
{
iRiverType = g_ModelInputPara.RiverProType;
dB = g_ModelInputPara.RiverProWidth;
dUnitQ = dGridQ / dB;//当前时段的侧向单宽入流
double slope = RoutePara[ord].pGridSlope;//取出坡度
int row = RoutePara[ord].pRow;
int col = RoutePara[ord].pCol;
int VegID = (int)g_GridLayerPara.g_VegLayer[row, col];
double dManning = g_GridLayerPara.VegTypes[VegID].dManning_n;
dVflow = 0.489 * Math.Pow(dUnitQ, 0.25) * Math.Pow(slope, 0.375) / Math.Pow(dManning, 0.75);
}
else //如果是坡面，计算坡面流速
{
iRiverType = g_ModelInputPara.HillProType;
dB = g_ModelInputPara.HillProWidth;
if (dB == 0)
dB = g_GridLayerPara.g_DemLayer.resolution;
//计算当前时段的侧向单宽入流
dUnitQ = dGridQ / dB;
//计算坡面流速
int row = RoutePara[ord].pRow;
int col = RoutePara[ord].pCol;
int VegID = (int)g_GridLayerPara.g_VegLayer[row, col];
double dManning = g_GridLayerPara.VegTypes[VegID].dManning_n;
dVflow = Math.Pow(dUnitQ, 0.4) * Math.Pow(RoutePara[ord].pGridSlope, 0.3) / Math.Pow(dManning, 0.6);
}
//准备Cunge法所需要的参数
MuskCunge.SetRoutingPara(dVflow, gridslp, deltaX, deltaT, iRiverType, Q11, Q12, Q21, dUnitQ, dB, x, k);
//计算Cunge法的洪水波流量
Q22 = MuskCunge.RoutingOutQ();
return Q22;
}
/// <summary>
/// 计算序号为ord的栅格的汇流结果，采用Cunge法
/// </summary>
/// <param name="ord">栅格的序号</param>
/// <param name="gridslp">网格坡度</param>
/// <param name="deltaX">河长的离散值</param>
/// <param name="deltaT">时间的离散值</param>
/// <param name="Q11">差分值</param>
/// <param name="Q12">差分值</param>
/// <param name="Q21">差分值</param>
/// <param name="dGridQ">当前时段的侧向流量</param>
/// <param name="TimeWeight">河段</param>
/// <param name="TravelTime">时段</param>
/// <returns>流出流量</returns>
public double GridRoutingOut_Liu(int ord, double deltaX, double deltaT, double Q11, double Q12, double Q21, double dGridQ, double TimeWeight, double TravelTime)
{
double dVflow = 0;
RiverType iRiverType;
double dB = 0;//用于存放河道或山的宽度
double dUnitQ = 0;
//如果是河道，计算河道流速
int row = River.CalcuOrderList[ord].row;
int col = River.CalcuOrderList[ord].col;
double slope = this.River.m_pSlope[row, col];//取出坡度
double dManning = 0.01;
int VegID = (int)g_GridLayerPara.g_VegLayer[row, col];
if (VegID!=-9999 && g_GridLayerPara.VegTypes[VegID] != null)
dManning = g_GridLayerPara.VegTypes[VegID].dManning_n;
if (River.m_pStrahlerOrd[row, col] > 0)
{
iRiverType = g_ModelInputPara.RiverProType;
dB = g_ModelInputPara.RiverProWidth;
dUnitQ = dGridQ / dB;//当前时段的侧向单宽入流
dVflow = 0.489 * Math.Pow(dUnitQ, 0.25) * Math.Pow(slope, 0.375) / Math.Pow(dManning, 0.75);
}
else //如果是坡面，计算坡面流速
{
iRiverType = g_ModelInputPara.HillProType;
dB = g_ModelInputPara.HillProWidth;
if (dB == 0)
dB = g_GridLayerPara.g_DemLayer.resolution;
//计算当前时段的侧向单宽入流
dUnitQ = dGridQ / dB;
//计算坡面流速
dVflow = Math.Pow(dUnitQ, 0.4) * Math.Pow(slope, 0.3) / Math.Pow(dManning, 0.6);
}
//准备Cunge法所需要的参数
double x = TimeWeight;
double k = TravelTime/60.0; //把分钟转化为小时
double denominator = k * (1 - x) + 0.5 * deltaT;
double C1 = (x * k + 0.5 * deltaT) / denominator;
double C2 = (0.5 * deltaT - x * k) / denominator;
double C3 = (k * (1 - x) - 0.5 * deltaT) / denominator;
double C4 = deltaT * deltaX / denominator;
//计算Cunge法的有限差分值
double Q22= C1 * Q11 + C2 * Q12 + C3 * Q21 + C4 * dUnitQ;
return Q22;
}
/* 马斯京根法(按子流域划分河段) */
public bool MuskingumRiverRouting(double dTLen, SortedList<DateTime, double[]> pNodeQ, double[] pX, double[] pK, int NodeNum, DateTime curOrder)
{
bool bret = true;
for (int i = 0; i < NodeNum; i++)
{
double QOut = 0;
if (i == 0)
QOut = pNodeQ[curOrder][i];
else
QOut = pNodeQ[curOrder][i] + pRiverRoute[i - 1].pRoute.OutFlux;
pRiverRoute[i].pRoute.InFlux = QOut;
pRiverRoute[i].pRoute.OutFlux = RiverRoutingOut(dTLen, pRiverRoute[i].pPreRoute.InFlux, pRiverRoute[i].pRoute.InFlux,
pRiverRoute[i].pPreRoute.OutFlux, pX[i], pK[i]);
pRiverRoute[i].pRoute.bCal = true;
}
for (int i = 0; i < NodeNum; i++)
{
pRiverRoute[i].pPreRoute.InFlux = pRiverRoute[i].pRoute.InFlux;
pRiverRoute[i].pPreRoute.OutFlux = pRiverRoute[i].pRoute.OutFlux;
pRiverRoute[i].pRoute.bCal = false;
}
return bret;
}
public double RiverRoutingOut(double deltaT, double dInFlux1, double dInFlux2, double dOutFlux1, double x, double k)
{
double dOutFlux2 = 0;
Muskingum.SetRoutingPara(deltaT, dInFlux1, dInFlux2, dOutFlux1, x, k);
dOutFlux2 = Muskingum.RoutingOutQ();
return dOutFlux2;
}
/// <summary>
/// 滞时演算法(单流域)
/// </summary>
/// <param name="GridQ"></param>
/// <param name="pRoute"></param>
/// <param name="dTlen"></param>
/// <param name="QType"></param>
/// <param name="curorder"></param>
/// <param name="totorder"></param>
/// <returns></returns>
public void PureLagGridRouting(GridLayer GridQ, double dTlen, RunOffSourceType QType, DateTime curorder)
{
this.PureLagGridRouting(GridQ, dTlen, QType, curorder, 1.0, WaterYearType.WATER_LOW_YEAR);
}
/// <summary>
/// 滞时演算法(单流域)
/// </summary>
/// <param name="GridQ">待处理的数据</param>
/// <param name="pRoute">路径像素集合</param>
/// <param name="dTlen">时段按小时计算的长度</param>
/// <param name="QType">径流的类型，如是壤中流，还是基流</param>
/// <param name="curorder">当前时次</param>
/// <param name="totorder"></param>
/// <param name="snowfactor"></param>
/// <param name="WaterYrType"></param>
/// <returns></returns>
public void PureLagGridRouting(GridLayer GridQ, double dhr, RunOffSourceType QType, DateTime curorder, double snowfactor, WaterYearType WaterYrType)
{
DateTime LagOrder;
double dGridOut = 0;
double dSurfQLoss = 0;
if (WaterYrType == WaterYearType.WATER_HIGH_YEAR)
{
dSurfQLoss = g_ModelInputPara.HighWaterSurfQLoss;
}
else if (WaterYrType == WaterYearType.WATER_MID_YEAR)
{
dSurfQLoss = g_ModelInputPara.MidWaterSurfQLoss;
}
else if (WaterYrType == WaterYearType.WATER_LOW_YEAR)
{
dSurfQLoss = g_ModelInputPara.LowWaterSurfQLoss;
}
int m_row = g_GridLayerPara.g_DemLayer.rowCount;
int m_col = g_GridLayerPara.g_DemLayer.colCount;
for (int row = 0; row < m_row; row++)
{
for (int col = 0; col < m_col; col++)
{
if (IfGridNeedCal(row, col) == false)
continue;
if (QType == RunOffSourceType.RUNOFF_ELEMENT_SURFQ)
{
int LagTime = (int)(g_GridLayerPara.g_RouteSurfQTime[row, col] / dhr); //滞后时间
TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
LagOrder = curorder.Add(LagSpan);//当前时段加上滞时
DateTime endTime = g_ClimatePara.pcpdata.saDate[g_ClimatePara.pcpdata.saDate.Count - 1];
if (LagOrder < endTime)
{
dGridOut = snowfactor * GridQ[row, col] / (dSurfQLoss * Math.Pow(g_GridLayerPara.g_RouteSurfQTime[row, col], 1));
if (this.MyOutletQ.ContainsKey(LagOrder))
this.MyOutletQ[LagOrder].m_pOutletSurfQ += dGridOut;
}
}
else if (QType == RunOffSourceType.RUNOFF_ELEMENT_LATERALQ)
{
int LagTime = (int)(g_GridLayerPara.g_RouteLatQTime[row, col] / dhr);
TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
LagOrder = curorder.Add(LagSpan);
DateTime endTime = g_ClimatePara.pcpdata.saDate[g_ClimatePara.pcpdata.saDate.Count - 1];
if (LagOrder < endTime)
{
dGridOut = snowfactor * GridQ[row, col] / (g_ModelInputPara.LatQLoss * Math.Pow(g_GridLayerPara.g_RouteLatQTime[row, col], 1));
if (this.MyOutletQ.ContainsKey(LagOrder))
this.MyOutletQ[LagOrder].m_pOutletLatQ += dGridOut;
}
}
else if (QType == RunOffSourceType.RUNOFF_ELEMENT_BASEQ)
{
int LagTime = (int)(g_GridLayerPara.g_RouteBaseQTime[row, col] / dhr);
TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
LagOrder = curorder.Add(LagSpan);
DateTime endTime = g_ClimatePara.pcpdata.saDate[g_ClimatePara.pcpdata.saDate.Count - 1];
if (LagOrder < endTime)
{
dGridOut = snowfactor * GridQ[row, col] / (g_ModelInputPara.BaseQLoss * Math.Pow(g_GridLayerPara.g_RouteBaseQTime[row, col], 1));
if (this.MyOutletQ.ContainsKey(LagOrder))
this.MyOutletQ[LagOrder].m_pOutletBaseQ += dGridOut;
}
}
}
}
}
/// <summary>
/// 滞时演算法(单流域)
/// </summary>
/// <param name="GridQ">待处理的数据</param>
/// <param name="pRoute">路径像素集合</param>
/// <param name="dTlen">时段按小时计算的长度</param>
/// <param name="QType">径流的类型，如是壤中流，还是基流</param>
/// <param name="curorder">当前时次</param>
/// <param name="totorder"></param>
/// <param name="snowfactor"></param>
/// <param name="WaterYrType"></param>
/// <returns></returns>
public void PureLagGridRouting_Liu(GridLayer GridQ, double dhr, RunOffSourceType QType, DateTime curorder, double snowfactor, WaterYearType WaterYrType)
{
DateTime LagOrder;
double dGridOut = 0;
double dSurfQLoss = 0;
if (WaterYrType == WaterYearType.WATER_HIGH_YEAR)
{
dSurfQLoss = g_ModelInputPara.HighWaterSurfQLoss;
}
else if (WaterYrType == WaterYearType.WATER_MID_YEAR)
{
dSurfQLoss = g_ModelInputPara.MidWaterSurfQLoss;
}
else if (WaterYrType == WaterYearType.WATER_LOW_YEAR)
{
dSurfQLoss = g_ModelInputPara.LowWaterSurfQLoss;
}
int m_row = g_GridLayerPara.g_DemLayer.rowCount;
int m_col = g_GridLayerPara.g_DemLayer.colCount;
for (int row = 0; row < m_row; row++)
{
for (int col = 0; col < m_col; col++)
{
if (IfGridNeedCal(row, col) == false)
continue;
if (QType == RunOffSourceType.RUNOFF_ELEMENT_SURFQ)
{
int LagTime = (int)(River.m_pSurfQRouteTime[row, col] / dhr); //滞后时间
TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
LagOrder = curorder.Add(LagSpan);//当前时段加上滞时
DateTime endTime = g_ClimatePara.pcpdata.saDate[g_ClimatePara.pcpdata.saDate.Count - 1];
if (LagOrder < endTime)
{
if (River.m_pSurfQRouteTime[row, col] <= 0)
continue;
dGridOut = snowfactor * GridQ[row, col] / (dSurfQLoss * Math.Pow(River.m_pSurfQRouteTime[row, col], 1));
if (this.MyOutletQ.ContainsKey(LagOrder))
this.MyOutletQ[LagOrder].m_pOutletSurfQ += dGridOut;
else
{
this.MyOutletQ[LagOrder] = new OutletQ();
this.MyOutletQ[LagOrder].m_pOutletSurfQ = dGridOut;
}
}
}
else if (QType == RunOffSourceType.RUNOFF_ELEMENT_LATERALQ)
{
int LagTime = (int)(River.m_pLatQRouteTime[row, col] / dhr);
TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
LagOrder = curorder.Add(LagSpan);
DateTime endTime = g_ClimatePara.pcpdata.saDate[g_ClimatePara.pcpdata.saDate.Count - 1];
if (LagOrder < endTime)
{
//无汇流
if (River.m_pLatQRouteTime[row, col] <= 0)
continue;
dGridOut = snowfactor * GridQ[row, col] / (g_ModelInputPara.LatQLoss * Math.Pow(River.m_pLatQRouteTime[row, col], 1));
//Console.Write(River.m_pBaseQRouteTime[row, col]+" ");
//Console.Write(g_ModelInputPara.LatQLoss+" ");
//Console.WriteLine(dGridOut);
if (this.MyOutletQ.ContainsKey(LagOrder))
this.MyOutletQ[LagOrder].m_pOutletLatQ += dGridOut;
else
{
this.MyOutletQ[LagOrder] = new OutletQ();
this.MyOutletQ[LagOrder].m_pOutletLatQ = dGridOut;
}
}
}
else if (QType == RunOffSourceType.RUNOFF_ELEMENT_BASEQ)
{
int LagTime = (int)(River.m_pBaseQRouteTime[row, col] / dhr);
TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
LagOrder = curorder.Add(LagSpan);
DateTime endTime = g_ClimatePara.pcpdata.saDate[g_ClimatePara.pcpdata.saDate.Count - 1];
if (LagOrder < endTime)
{
//无汇流
if (River.m_pBaseQRouteTime[row, col] <= 0)
continue;
dGridOut = snowfactor * GridQ[row, col] / (g_ModelInputPara.BaseQLoss * Math.Pow(River.m_pBaseQRouteTime[row, col], 1));
if (this.MyOutletQ.ContainsKey(LagOrder))
this.MyOutletQ[LagOrder].m_pOutletBaseQ += dGridOut;
else
{
this.MyOutletQ[LagOrder] = new OutletQ();
this.MyOutletQ[LagOrder].m_pOutletBaseQ = dGridOut;
}
}
}
}
}
}
/* 滞时演算法(分子流域) */
public bool PureLagGridRouting_Node(GridLayer GridQ, SortedList<DateTime, double[]> pRoute, double dTlen,
RunOffSourceType QType,
DateTime curorder,
int totorder,
double snowfactor,
WaterYearType WaterYrType)
{
DateTime LagOrder;
int LagTime = 0;
double dGridOut = 0;
int subID = 0;
double dSurfQLoss = 0;
switch (WaterYrType)
{
case WaterYearType.WATER_HIGH_YEAR:
{
dSurfQLoss = g_ModelInputPara.HighWaterSurfQLoss;
break;
}
case WaterYearType.WATER_MID_YEAR:
{
dSurfQLoss = g_ModelInputPara.MidWaterSurfQLoss;
break;
}
case WaterYearType.WATER_LOW_YEAR:
{
dSurfQLoss = g_ModelInputPara.LowWaterSurfQLoss;
break;
}
}
int m_row = g_GridLayerPara.g_DemLayer.rowCount;
int m_col = g_GridLayerPara.g_DemLayer.colCount;
for (int i = 0; i < m_row; i++)
{
for (int j = 0; j < m_col; j++)
{
if (g_GridLayerPara.g_BasinBoundary.Values[i, j] == 1.0)
{
subID = (int)g_GridLayerPara.g_SubWaterShed.Values[i, j] - 1;
switch (QType)
{
case RunOffSourceType.RUNOFF_ELEMENT_SURFQ:
{
LagTime = (int)(g_GridLayerPara.g_RouteSurfQTime.Values[i, j] / dTlen);
TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
LagOrder = curorder.Add(LagSpan);
//if (LagOrder < totorder)
//{
if (g_GridLayerPara.g_RouteSurfQTime.Values[i, j] <= 0)
g_GridLayerPara.g_RouteSurfQTime.Values[i, j] = 0.1;
dGridOut = snowfactor * GridQ.Values[i, j] / (dSurfQLoss * Math.Pow(g_GridLayerPara.g_RouteSurfQTime.Values[i, j], 1));
pRoute[LagOrder][subID] += dGridOut;
//}
}
break;
case RunOffSourceType.RUNOFF_ELEMENT_LATERALQ:
{
LagTime = (int)(g_GridLayerPara.g_RouteLatQTime.Values[i, j] / dTlen);
TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
LagOrder = curorder.Add(LagSpan);
//if (LagOrder < totorder)
//{
if (g_GridLayerPara.g_RouteLatQTime.Values[i, j] <= 0)
g_GridLayerPara.g_RouteLatQTime.Values[i, j] = 0.1;
dGridOut = snowfactor * GridQ.Values[i, j] / (g_ModelInputPara.LatQLoss * Math.Pow(g_GridLayerPara.g_RouteLatQTime.Values[i, j], 1));
pRoute[LagOrder][subID] += dGridOut;
//}
}
break;
case RunOffSourceType.RUNOFF_ELEMENT_BASEQ:
{
LagTime = (int)(g_GridLayerPara.g_RouteBaseQTime.Values[i, j] / dTlen);
TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
LagOrder = curorder.Add(LagSpan);
//if (LagOrder < totorder)
//{
if (g_GridLayerPara.g_RouteBaseQTime.Values[i, j] <= 0)
g_GridLayerPara.g_RouteBaseQTime.Values[i, j] = 0.1;
dGridOut = snowfactor * snowfactor * GridQ.Values[i, j] / (g_ModelInputPara.BaseQLoss * Math.Pow(g_GridLayerPara.g_RouteBaseQTime.Values[i, j], 1));
pRoute[LagOrder][subID] += dGridOut;
//}
}
break;
default:
return false;
}
}
}
}
return true;
}
/// <summary>
/// 计算植被逐日Albedo
/// </summary>
/// <param name="mon"></param>
/// <param name="day"></param>
/// <returns></returns>
public double GetVegAlbedo(int VegID, int mon, int day)
{
double result = 0.23; //取默认值为0.23
//取真实测量的值
if (mon >= 1 || mon <= 12)
{
result = g_GridLayerPara.VegTypes[VegID].Albedo[mon - 1];
}
return result;
}
public void LoadSystemParaFile(string path)
{
g_ClimatePara.ReadClimateFile(path + "\\SysPara\\Climate.txt");
g_GridLayerPara.ReadGridIOFile(path);
//g_GridLayerPara.ReadBasinExtent(m_BasinFile);
g_ModelRunPara.ReadPara(m_RunParaFile);
g_ModelInputPara.ReadPara(m_InputParaFile, g_ModelInputPara);
g_GridResultOut.ReadMidOutFile(m_MidOutFile);
}
/// <summary>
/// 用于计算河道流速
/// </summary>
/// <param name="dUnitQ"></param>
/// <param name="slp">坡度</param>
/// <param name="manning_n">曼宁系数</param>
/// <returns>流速</returns>
public double ChannelQVelocity(double dUnitQ, double slp, double manning_n)
{
double dVflow;
dVflow = 0.489 * Math.Pow(dUnitQ, 0.25) * Math.Pow(slp, 0.375) / Math.Pow(manning_n, 0.75);
return dVflow;
}
/// <summary>
/// 计算坡面流速
/// </summary>
/// <param name="dUnitQ"></param>
/// <param name="slp">坡度坡面流速</param>
/// <param name="manning_n"></param>
/// <returns></returns>
public double HillSideQVelocity(double dUnitQ, double slp, double manning_n)
{
double dVflow;
dVflow = Math.Pow(dUnitQ, 0.4) * Math.Pow(slp, 0.3) / Math.Pow(manning_n, 0.6);
return dVflow;
}
/// <summary>
/// 霍顿下渗曲线计算暴雨过程主函数(刘永和)
/// </summary>
public virtual void StormRunoffSim_Horton_Liu()
{
int rowCount = g_GridLayerPara.g_DemLayer.rowCount;
int colCount = g_GridLayerPara.g_DemLayer.colCount;
this.CellOutSurfFlux = new double[rowCount, colCount];
this.CellOutLatFlux = new double[rowCount, colCount];
this.CellOutBaseFlux = new double[rowCount, colCount];
GridLayerInit();
//下面这段程序由刘永和所加
int n = River.CalcuOrderList.Count;
CurrentSurfQ = new RouteFlux[n];
CurrentLatQ = new RouteFlux[n];
CurrentBaseQ = new RouteFlux[n];
PreSurfQ = new RouteFlux[n];
PreLatQ = new RouteFlux[n];
PreBaseQ = new RouteFlux[n];
for (int i = 0; i < n; i++)
{
CurrentSurfQ[i] = new RouteFlux();
CurrentLatQ[i] = new RouteFlux();
CurrentBaseQ[i] = new RouteFlux();
PreSurfQ[i] = new RouteFlux();
PreLatQ[i] = new RouteFlux();
PreBaseQ[i] = new RouteFlux();
}
//计算间隔时段长度
StreamWriter sw = new StreamWriter(HydroSimulate.path + "\\Output\\OutletQ.txt");
sw.WriteLine("Order\tTotalQ\tSurfQ\tLatQ\tBaseQ\tDeepBaseQ");
//初始化指定出口断面计算结果内存
this.MyOutletQ = new SortedList<DateTime, OutletQ>();
//对各降雨时段进行循环
for (int i = 0; i < g_ClimatePara.pcpdata.saDate.Count; i++)
{
DateTime dateTime = g_ClimatePara.pcpdata.saDate[i];
//取得该降雨时段的dn历
int dn = DateAndTime.GetDnInYear(dateTime.Year, dateTime.Month, dateTime.Day);
//对该时段的降雨和气象进行插值
this.MeteoDataInterp(dateTime);
//计算该时段的小时数
double dhr = 1;
ProgressBar bar = new ProgressBar();
bar.Text = "正在计算净雨量";
bar.Show();
for (int row = 0; row < rowCount; row++)
{
for (int col = 0; col < colCount; col++)
{
//计算本单元的水量平衡
this.CalcCell_Horton(row, col, dhr, dateTime, dn);
}
bar.progressBar1.Value = (int)(row * 100.0 / rowCount);
}
bar.Close();//关闭计算净雨时的bar
//下面开始汇流
this.CalcRiverRouting_Liu(dateTime, dhr, dn);
sw.Write(dateTime + "\t");
sw.Write(this.MyOutletQ[dateTime].m_pOutletQ.ToString("0.0000") + "\t");
sw.Write(this.MyOutletQ[dateTime].m_pOutletSurfQ.ToString("0.0000") + "\t");
sw.Write(this.MyOutletQ[dateTime].m_pOutletLatQ.ToString("0.0000") + "\t");
sw.Write(this.MyOutletQ[dateTime].m_pOutletBaseQ.ToString("0.0000") + "\t");
sw.Write(this.MyOutletQ[dateTime].m_pOutletDeepBaseQ.ToString("0.0000") + "\t");
sw.WriteLine();
sw.Flush();
MidGridResultOut(dateTime, false);
}
sw.Close();
RiverOutletQ(true);
}
/// <summary>
/// 进行汇流计算
/// </summary>
/// <param name="curorder">当前的执行时段序号</param>
/// <param name="dhr">当前时段的小时数</param>
/// <param name="totrec">总时段数</param>
/// <param name="dn">年内的日序</param>
public void CalcRiverRouting_Liu(DateTime curorder, double dhr, int dn)
{
//下面进行坡面漫流汇流
if (g_ModelRunPara.SurfRouteMethod == SurfRouteMehtod.ROUTE_MUSK_CONGE)
{
RouteFlux[] surfQ = this.MuskCungeSurfRouting_Liu(dhr, g_ModelInputPara.SMCTimeWeight, g_ModelInputPara.SMCGridRTravelTime);
for (int k = 0; k < River.CalcuOrderList.Count; k++)
{
int row =River.CalcuOrderList[k].row;
int col = River.CalcuOrderList[k].col;
g_GridResultOut.m_GridRoutingQ[row, col] = surfQ[k].OutFlux;
}
//把路径上最后一个像素的出流量作为总出流量
if (this.MyOutletQ.ContainsKey(curorder) == false)
this.MyOutletQ[curorder] = new OutletQ();
this.MyOutletQ[curorder].m_pOutletSurfQ = surfQ[River.CalcuOrderList.Count - 1].OutFlux;
}
else
{
this.PureLagGridRouting_Liu(g_GridResultOut.m_GridSurfQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_SURFQ, curorder, 1.0, WaterYearType.WATER_LOW_YEAR);
}
//下面进行壤中流汇流
if (g_ModelRunPara.LatRouteMethod == SurfRouteMehtod.ROUTE_MUSK_CONGE)
{
RouteFlux[] latQ = MuskCungeLatRouting_Liu(dhr, g_ModelInputPara.LMCTimeWeight, g_ModelInputPara.LMCGridRTravelTime);
for (int k = 0; k < River.CalcuOrderList.Count; k++)
{
int row = River.CalcuOrderList[k].row;
int col = River.CalcuOrderList[k].col;
g_GridResultOut.m_GridRoutingQ[row, col] += latQ[k].OutFlux;
}
if (this.MyOutletQ.ContainsKey(curorder) == false)
this.MyOutletQ[curorder] = new OutletQ();
this.MyOutletQ[curorder].m_pOutletLatQ = latQ[River.CalcuOrderList.Count - 1].OutFlux;
}
else
{
PureLagGridRouting_Liu(g_GridResultOut.m_GridLateralQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_LATERALQ, curorder, 1.0, WaterYearType.WATER_LOW_YEAR);
}
//下面进行地下径流的汇流
if (g_ModelRunPara.BaseRouteMethod == SurfRouteMehtod.ROUTE_MUSK_CONGE)
{
RouteFlux[] baseQ = this.MuskCungeBaseRouting_Liu(dhr, g_ModelInputPara.BMCTimeWeight, g_ModelInputPara.BMCGridRTravelTime);
for (int k = 0; k < River.CalcuOrderList.Count; k++)
{
int row = River.CalcuOrderList[k].row;
int col = River.CalcuOrderList[k].col;
g_GridResultOut.m_GridRoutingQ[row, col] += baseQ[k].OutFlux;
}
if (this.MyOutletQ.ContainsKey(curorder) == false)
this.MyOutletQ[curorder] = new OutletQ();
this.MyOutletQ[curorder].m_pOutletBaseQ = baseQ[River.CalcuOrderList.Count - 1].OutFlux;
}
else
{
PureLagGridRouting_Liu(g_GridResultOut.m_GridBaseQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_BASEQ, curorder, 1.0, WaterYearType.WATER_LOW_YEAR);
}
//下面计算深层地下径流的汇流
if (!this.MyOutletQ.ContainsKey(curorder))
this.MyOutletQ[curorder] = new OutletQ();
this.MyOutletQ[curorder].m_pOutletDeepBaseQ = DeepBaseQSim(dn, g_ModelInputPara.DeepBaseQ);
this.MyOutletQ[curorder].m_pOutletQ = this.MyOutletQ[curorder].m_pOutletSurfQ * g_ModelInputPara.SurfQLinearFactor
+ this.MyOutletQ[curorder].m_pOutletLatQ + this.MyOutletQ[curorder].m_pOutletBaseQ + this.MyOutletQ[curorder].m_pOutletDeepBaseQ;
}
/// <summary>
/// 马斯京根-康吉法汇流计算(先合后演)
/// </summary>
/// <param name="dT">时间的离散值</param>
/// <param name="GridQ"></param>
/// <param name="pRoute"></param>
/// <param name="pPreRoute"></param>
/// <param name="x"></param>
/// <param name="k"></param>
/// <returns></returns>
public RouteFlux[] MuskCungeSurfRouting_Liu(double dT, double TimeWeight, double TravelTime)
{
ProgressBar bar = new ProgressBar();
bar.Text = "正在计算汇流";
bar.Show();
int n = River.CalcuOrderList.Count;
for (int i = 0; i < n; i++)
{
int row=River.CalcuOrderList[i].row;
int col=River.CalcuOrderList[i].col;
//先统计当前栅格(即第i个)的总入流量
double SumQ = 0;
for (int index = 0; index < 9; index++)
{
if (index == 4)
continue;
int R = row + index / 3 - 1;
int C = col + index % 3 - 1;
if (R < 0 || R >= River.demgrid.rowCount)
continue;
if (C < 0 || C >= River.demgrid.colCount)
continue;
if (River.m_pDirection[R, C] == 8 - index)
{
SumQ += this.CellOutSurfFlux[R, C];
}
}
//if (SumQ < 0)
// throw new Exception("");
CurrentSurfQ[i].InFlux = SumQ;
//再使用Cunge法计算出流量
double dX=River.demgrid.resolution;//栅格空间上河流长度值
if(River.m_pDirection[ row, col]%2==0)
dX*=1.41421;
//如果当前栅格无数据
if (g_GridResultOut.m_GridSurfQ[row, col] < 0)
throw new Exception("流量值为空");
double tempGridQ = g_GridResultOut.m_GridSurfQ[row, col] * Math.Pow(g_GridResultOut.m_GridSurfQ.resolution, 2) / (1000 * dT * 3600);//栅格面内入流(单位：立方米/秒)
CurrentSurfQ[i].OutFlux = this.GridRoutingOut_Liu(i, dX, dT, PreSurfQ[i].InFlux, CurrentSurfQ[i].InFlux, PreSurfQ[i].OutFlux, tempGridQ, TimeWeight, TravelTime);
this.CellOutSurfFlux[row, col] = CurrentSurfQ[i].OutFlux;
bar.progressBar1.Value = (int)(100.0 * i / n);
}
//将CurrentQ和PreQ交换
RouteFlux[] result = CurrentSurfQ;
RouteFlux[] temp = PreSurfQ;
PreSurfQ = CurrentSurfQ;
CurrentSurfQ = temp;
bar.Close();
return result;
}
public RouteFlux[] MuskCungeLatRouting_Liu(double dT, double TimeWeight, double TravelTime)
{
ProgressBar bar = new ProgressBar();
bar.Text = "正在计算汇流";
bar.Show();
int n = River.CalcuOrderList.Count;
for (int i = 0; i < n; i++)
{
int row = River.CalcuOrderList[i].row;
int col = River.CalcuOrderList[i].col;
//先统计当前栅格(即第i个)的总入流量
double SumQ = 0;
for (int index = 0; index < 9; index++)
{
if (index == 4)
continue;
int R = row + index / 3 - 1;
int C = col + index % 3 - 1;
if (R < 0 || R >= River.demgrid.rowCount)
continue;
if (C < 0 || C >= River.demgrid.colCount)
continue;
if (River.m_pDirection[R, C] == 8 - index)
{
if (this.CellOutLatFlux[R, C] < 0)
throw new Exception("");
SumQ += this.CellOutLatFlux[R, C];
}
}
if (SumQ < 0)
throw new Exception("");
CurrentLatQ[i].InFlux = SumQ;
//再使用Cunge法计算出流量
double dX = River.demgrid.resolution;//栅格空间上河流长度值
if (River.m_pDirection[row, col] % 2 == 0)
dX *= 1.41421;
//如果当前栅格无数据
if (g_GridResultOut.m_GridLateralQ[row, col] < 0)
throw new Exception("流量值为空");
double tempGridQ = g_GridResultOut.m_GridLateralQ[row, col] * Math.Pow(River.demgrid.resolution, 2) / (1000 * dT * 3600);//栅格面内入流(单位：立方米/秒)
CurrentLatQ[i].OutFlux = this.GridRoutingOut_Liu(i, dX, dT, PreLatQ[i].InFlux, CurrentLatQ[i].InFlux, PreLatQ[i].OutFlux, tempGridQ, TimeWeight, TravelTime);
if (CurrentLatQ[i].OutFlux < 0)
{
if (CurrentLatQ[i].OutFlux < -0.1)
throw new Exception("");
CurrentLatQ[i].OutFlux = 0;
}
this.CellOutLatFlux[row, col] = CurrentLatQ[i].OutFlux;
bar.progressBar1.Value = (int)(100.0 * i / n);
}
//将CurrentQ和PreQ交换
RouteFlux[] result = CurrentLatQ;
RouteFlux[] temp = PreLatQ;
PreLatQ = CurrentLatQ;
CurrentLatQ = temp;
bar.Close();
return result;
}
public RouteFlux[] MuskCungeBaseRouting_Liu(double dT, double TimeWeight, double TravelTime)
{
ProgressBar bar = new ProgressBar();
bar.Text = "正在计算汇流";
bar.Show();
int n = River.CalcuOrderList.Count;
for (int i = 0; i < n; i++)
{
int row = River.CalcuOrderList[i].row;
int col = River.CalcuOrderList[i].col;
//先统计当前栅格(即第i个)的总入流量
double SumQ = 0;
for (int index = 0; index < 9; index++)
{
if (index == 4)
continue;
int R = row + index / 3 - 1;
int C = col + index % 3 - 1;
if (R < 0 || R >= River.demgrid.rowCount)
continue;
if (C < 0 || C >= River.demgrid.colCount)
continue;
if (River.m_pDirection[R, C] == 8 - index)
{
if (this.CellOutBaseFlux[R, C] < 0)
throw new Exception("");
SumQ += this.CellOutBaseFlux[R, C];
}
}
if (SumQ < 0)
throw new Exception("");
CurrentBaseQ[i].InFlux = SumQ;
//再使用Cunge法计算出流量
double dX = River.demgrid.resolution;//栅格空间上河流长度值
if (River.m_pDirection[row, col] % 2 == 0)
dX *= 1.41421;
//如果当前栅格无数据
if (g_GridResultOut.m_GridBaseQ[row, col] < 0)
throw new Exception("流量值为空");
double tempGridQ = g_GridResultOut.m_GridBaseQ[row, col] * Math.Pow(River.demgrid.resolution, 2) / (1000 * dT * 3600);//栅格面内入流(单位：立方米/秒)
CurrentBaseQ[i].OutFlux = this.GridRoutingOut_Liu(i, dX, dT, PreBaseQ[i].InFlux, CurrentBaseQ[i].InFlux, PreBaseQ[i].OutFlux, tempGridQ, TimeWeight, TravelTime);
if (CurrentBaseQ[i].OutFlux < 0)
{
if (CurrentBaseQ[i].OutFlux < -0.1)
throw new Exception("");
CurrentBaseQ[i].OutFlux = 0;
}
this.CellOutBaseFlux[row, col] = CurrentBaseQ[i].OutFlux;
bar.progressBar1.Value = (int)(100.0 * i / n);
}
//将CurrentQ和PreQ交换
RouteFlux[] result = CurrentBaseQ;
RouteFlux[] temp = PreBaseQ;
PreBaseQ = CurrentBaseQ;
CurrentBaseQ = temp;
bar.Close();
return result;
}
/// <summary>
/// 汇流计算变量
/// </summary>
public MuskCungeRoutePara[] RoutePara;
public MuskingumCunge MuskCunge = new MuskingumCunge();
public MuskRouteFlux[] pSurfQ;
public MuskRouteFlux[] pLatQ;
public MuskRouteFlux[] pBaseQ;
public double[] pRouteQ;
public Muskingum Muskingum = new Muskingum();
public MuskRouteFlux[] pRiverRoute;
public int m_iNodeNum;
public double[] m_pX;
public double[] m_pK;
/// <summary>
/// 水量平衡变量
/// </summary>
//public GridWaterBalance gridwb;
//public GreenAmptInfil GreenAmptInfil;
//public HortonInfil HortonInfil;
/// <summary>
/// 土壤水下渗
/// </summary>
public static GridLayer m_drateinf;
public static GridLayer m_drintns;
public static GridLayer m_dcumr;
public static GridLayer m_dcuminf;
public static GridLayer m_dexcum;
//计算结果输出
public int m_OutRow;
public int m_OutCol;
public SortedList<DateTime, OutletQ> MyOutletQ;
/// <summary>
/// 与汇流次序表对应的当前流量表(刘永和加)
/// </summary>
public RouteFlux[] CurrentSurfQ;
public RouteFlux[] CurrentLatQ;
public RouteFlux[] CurrentBaseQ;
/// <summary>
/// 与汇流次序表对应的前一时刻流量表(刘永和加)
/// </summary>
public RouteFlux[] PreSurfQ;
public RouteFlux[] PreLatQ;
public RouteFlux[] PreBaseQ;
public double[,] CellOutSurfFlux;
public double[,] CellOutLatFlux;
public double[,] CellOutBaseFlux;
//按子流域计算存储
SortedList<DateTime, double[]> m_pNodeSurfQ;
SortedList<DateTime, double[]> m_pNodeLatQ;
SortedList<DateTime, double[]> m_pNodeBaseQ;
SortedList<DateTime, double[]> m_pNodeOutQ;
/// <summary>
/// 这一时段的时长，如24小时或1小时，单位：分钟
/// </summary>
public SortedList<DateTime, double> m_pRainInterval;
int m_iTotalRec;
//int m_row;
//int m_col;
int m_iVegOrd;
int m_iSoilOrd;
public bool m_bDate;
string m_strBeg;
string m_strEnd;
/// <summary>
/// 丰水年或平水年或少水年(其中键为年份值)
/// </summary>
SortedList<int, WaterYearType> waterYearType;
public string m_GridFile; // GridIO.txt
public string m_BasinFile; // Basin.txt
public string m_ClimateFile; // Climate.txt
public string m_RunParaFile; // RunPara.txt
public string m_InputParaFile; // InputPara.txt
public string m_MuskCoeffFile; // muskingum_coeff.txt
public string m_MidOutFile; // MidGridOut.txt
public string m_WaterYrTypeFile; // WaterYearType.txt
}
}
Вы можете оставить комментарий после Вход в систему
Неприемлемый контент может быть отображен здесь и не будет показан на странице. Вы можете проверить и изменить его с помощью соответствующей функции редактирования.
Если вы подтверждаете, что содержание не содержит непристойной лексики/перенаправления на рекламу/насилия/вульгарной порнографии/нарушений/пиратства/ложного/незначительного или незаконного контента, связанного с национальными законами и предписаниями, вы можете нажать «Отправить» для подачи апелляции, и мы обработаем ее как можно скорее.
Опубликовать ( 0 )