水文預報課程設計河海

1、第一章 基本任務本次課程設計從2013年7月8日至2012年7月12日，主要任務是對廣東省東江一級支流西枝江白盆珠水庫的上游寶口流域編制預報方案與產(chǎn)匯流計算。其基本任務為：任務一：根據(jù)已給的資料、參數(shù)及做過的習題，自己編寫程序，將流域作為整體進行產(chǎn)流量計算；將計算年徑流與實測年徑流進行比較；每人計算兩年。任務二：根據(jù)已給設計暴雨資料、參數(shù)及做過的習題,自己編寫程序，將流域作為整體進行次洪產(chǎn)流量、劃分水源、直接徑流匯流、地下徑流匯流計算；繪出直接徑流過程、地下徑流過程、總的流量過程。第二章 基本資料2.1 流域概況白盆珠水庫位于廣東省東江一級支流西枝江的上游，壩址以上集雨面積856km2 。流域

2、地處粵東沿海的西部，海洋性氣候顯著，氣候溫和，雨量豐沛。暴雨成因主要是鋒面雨和臺風雨，常受熱帶風暴影響。降雨年際間變化大，年內(nèi)分配不均，多年平均降雨量為1800mm，實測年最大降雨量為3417mm，汛期49月降雨量占年降雨量的81左右；徑流系數(shù)0.5-0.7。流域內(nèi)地勢平緩，土壤主要有黃壤和砂壤，具有明顯的腐殖層，淀積層和母質(zhì)土等層次結構，透水性好。臺地、丘陵多生長松、杉、樟等高大喬木；平原則以種植農(nóng)作物和經(jīng)濟作物為主，植被良好。 流域上游有寶口水文站，流域面積為553km2，占白盆珠水庫壩址以上集雨面積的64.6。白盆珠水庫有6年逐日入庫流量資料、逐日蒸發(fā)資料和時段入庫流量資料：流域內(nèi)有7個

3、雨量站，其中寶口以上有4個。雨量站分布較均勻，有6年逐日降水資料和時段降水資料；寶口水文站具有6年以上水位、流量資料；流域?qū)偕絽^(qū)性小流域且受到地形、地貌等下墊面條件影響，洪水陡漲緩落，匯流時間一般23小時（h），有時更短；一次洪水總歷時25天（d）。2.2 基本資料2.2.1 任務一參數(shù)及相關實測資料計算流域為白盆珠水庫上游的寶口流域，其流域面積為553km2。該流域內(nèi)有四個雨量站：禾多布、馬山、高潭、寶口，其權重系數(shù)分別為：0.33、0.14、0.33、0.20。本次課設提供的資料為1987年至1992年該流域徑流量、蒸發(fā)量及各雨量站的日實測資料，其中我計算年份為1987年與1988年，所使

4、用的相關參數(shù)值如下表所示：WM(mm)WUM(mm)WLM(mm)BCIMPW(mm)WU(mm)WL(mm)WD(mm)14020600.20.160.001110104060其中蒸散發(fā)折算系數(shù)Kc需編程優(yōu)選，優(yōu)選范圍為0.90-1.30。優(yōu)選原則為：計算的2年內(nèi)每年的年徑流相對誤差盡可能不超過5%。2.2.2 任務二參數(shù)及相關實測資料任務二產(chǎn)流參數(shù)Kc為任務一中日模優(yōu)選的值，F(xiàn)C為10.0mm/3h，產(chǎn)流參數(shù)Kc用日模優(yōu)選的值，其他參數(shù)取任務1中各組的參數(shù)，初始張力水蓄量取各組的容量值。任務二匯流參數(shù)CG為0.968，QG0=55.3m3/s，單位線序號1至11對應的Q（m3/s）分別為：

5、0、40、80、130、100、80、48、20、10、5、02.3 產(chǎn)流方式論證寶口流域地處我國南方濕潤地區(qū)，氣候暖熱，雨量充沛，多年平均降雨量為1800mm>1000mm，年徑流系數(shù)在0.50.7之間，大于0.4；流域內(nèi)土質(zhì)疏松，植被良好，不易超滲；一次洪水的流量過程陡漲緩落，持續(xù)時間25d左右。從流域的氣象條件，下墊面條件和流量過程的分析知，該流域降雨徑流關系具有蓄滿產(chǎn)流的特點，可以按蓄滿產(chǎn)流建立產(chǎn)流量預報方案。第三章 計算公式3.1 降雨量計算由已知資料知，該地區(qū)雨量站分布均勻，且寶口以上有四個雨量站：禾多布、馬山、高潭、寶口，其權重系數(shù)分別為：0.33、0.14、0.33、0.

6、20，所以寶口流域平均降雨量P=0.33*P1+0.14*P2+0.33*P3+0.20*P4。3.2 蒸散發(fā)量計算蒸散發(fā)計算采用三層蒸發(fā)計算模式，即：上層蒸發(fā)量：EU=Ep下層蒸發(fā)量：EL=Ep*WL/WLM深層蒸發(fā)量：ED=C*Ep總蒸發(fā)量：E=EU+EL+ED式中：Ep為流域蒸發(fā)能力（mm）；WL為下層土壤含水量（mm）；WLM為下層土壤含水量（mm）；C為蒸發(fā)擴散系數(shù)。三層蒸發(fā)模式按照先上層后下層的次序，具體計算為：1)當WU+P>=EP時，EU=Ep,EL=0,ED=02)當WU+P<EP，WL>=C*WLM時，EU=WU+P,EL=(EP-EU)*WL/WLM,E

7、D=03)當WU+P<EP, C*(EP-EU)<=WL<C*WLM時，EU=WU+P,EL=C*(EP-EU),ED=04)當WU+P<EP, WL<C*(EP-EU)時，EU=WU+P,EL=WL,ED=C*(EP-EU)-EL3.3 降雨產(chǎn)流量計算1） aPEWMM時，R=PEWWMWM（1)b+1 2） aPEWMM時，R=PE（WMW）3.4 水源劃分計算通過穩(wěn)滲率fc可劃分產(chǎn)流中的直接徑流和地下徑流。次洪的各水源分量為： 3.5 匯流計算根據(jù)流域凈雨和流域徑流單位線，采用卷積的差分形式算出流域出口的流量過程。其計算公式：直接徑流過程：QS(i)=RS(

8、i)*U地下徑流過程：QG(i)=CG*QG(i-1)+(1-CG)*RG(i)*U總的流量過程：Q(i)=QS(i)+QG(i)第四章 計算結果4.1 任務一結果經(jīng)過程序優(yōu)選，Kc=0.98年份實測R（mm）計算R（mm）絕對誤差（mm）相對誤差（%）19871080110019.31.791988824813-10.2-1.234.2 任務二結果t(i)P(i)R(i)Rd(i)Rg(i)QS(i)QG(i)Q(i)2004/9/23 12:001412.72.710069.969.92004/9/23 15:0014.112.92.91010.884.194.92004/9/23 18:

9、00119.709.733.297.31302004/9/23 21:0018.717.47.41058.31111692004/9/24 0:0029.728.518.51094.31232182004/9/24 3:003635.125.1101841363202004/9/24 6:0038.337.527.5103811485292004/9/24 9:007.86.906.96451547992004/9/24 12:0030.529.619.6107991669652004/9/24 15:0042.641.731.71087417710502004/9/24 18:0093.89

10、2.982.91086418810502004/9/24 21:0084.183.373.310122019814202004/9/25 0:0047.646.836.810177020819802004/9/25 3:0056.455.345.310237021825902004/9/25 6:0050.449.339.310265022728702004/9/25 9:0016.515.45.410260023628402004/9/25 12:005.84.804.8236023726002004/9/25 15:008.37.207.2186024121002004/9/25 18:0

11、02.61.501.5125023614802004/9/25 21:000.200077422810002004/9/26 0:0000003962216172004/9/26 3:0000001682143822004/9/26 6:00000072.82072802004/9/26 9:00000025.02002252004/9/26 12:0000002.701941972004/9/26 15:0000000.001881882004/9/26 18:0000000.001821822004/9/26 21:0000000.00176176第五章 誤差統(tǒng)計與分析5.1 精度評定從計

12、算結果可見，年產(chǎn)流量絕對誤差均小于100mm，所產(chǎn)流量的相對誤差均小于5%。精度統(tǒng)計表明，率定的模型參數(shù)是基本合理的。但由于課設時間限制與任務要求，每位同學只用兩年的資料來率定Kc，不滿足水文情報預報規(guī)范中規(guī)定：洪水預報方案要求使用樣本數(shù)量不少于10年的水文氣象資料，其中應包括大、中、小水各種代表性年份，并保證有足夠代表性的場次洪水資料。顯然Kc=0.98的結果還是存在一定問題的。5.2 誤差分析影響流域降雨徑流過程的因素很多，利用蓄滿產(chǎn)流新安江模型的結構與參數(shù)能夠很好反映濕潤地區(qū)降雨徑流過程的主要規(guī)律和特點，因而能獲得較好的精度。但是模型本身以及模型計算中有很多概化，會造成誤差。造成寶口流域

13、產(chǎn)匯流計算方案誤差來源主要有以下幾個方面：1.資料代表性的影響。我們每位同學只用兩年實測資料來率定Kc，資料必不能滿足時期要求及代表性要求。實際操作中，應要求有12年以上的連續(xù)的對未來有代表性的實測資料，其中10年為率定期，2年為檢驗期，且這種資料應具備豐、平、枯水年代表性，資料系列前后一致，受人類影響較小。所以僅用兩年的資料率定勢必造成Kc一定的誤差，而且誤差可能還很大。2.量測誤差。實測的降雨、蒸發(fā)、徑流量等水文氣象信息及河流、湖泊、地形等下墊面信息是研制預報模型或編制洪水預報方案或進行作業(yè)預報的主要依據(jù)，在現(xiàn)有站網(wǎng)、儀器設備、觀測技術條件下，各種信息的時空變化是難以準確反映的，加上受自然

14、因素等客觀條件影響，勢必造成各種信息的量測誤差。3.模型結構誤差。在該蓄滿產(chǎn)流新安江模型中，有很多將非線性現(xiàn)象概化為線性現(xiàn)象或者將某些隨機因子近似作為確定因子描述等都會帶來誤差。4.模型參數(shù)誤差。模型參數(shù)是根據(jù)輸入，通過模型計算輸出，再將輸出過程與實測過程進行比較，用系統(tǒng)識別的方法作優(yōu)化調(diào)試的，上述所率定出的模型參數(shù)Kc可能不是最優(yōu)。5.人類活動的影響。隨著社會經(jīng)濟的快速發(fā)展，人類活動的影響加劇，流域內(nèi)可能新建了一些大中型水庫或其他工程措施。影響了該地區(qū)原有的產(chǎn)流特點，如汛期提前。5.3 實時校正模型計算值與實測值直接總是存在一定的誤差。造成兩者間誤差的因素很多，若針對每一個單一因素是難于描述

15、或預見的，一般采用實時校正模型來解決。實時校正模型常用的有卡爾曼濾波、自回歸模型等。第六章 編寫程序FORM 1Private Sub Command1_Click(Index As Integer) Form2.ShowEnd SubPrivate Sub Command2_Click() Form3.ShowEnd SubFORM 2 Private Sub Command1_Click() Form1.Show Form2.Hide End SubPrivate Sub Command2_Click()For m = 0 To 40 Kc = 0.9 + 0.01 * mReDim Ro

16、bs(2), Rcal(2), jderror(2), xderror(3) Dim W0 As Single, WU0 As Single, WL0 As Single, WD0 As Single WU0 = Text10.Text: WL0 = Text11.Text: WD0 = Text12.Text W0 = WL0 + WD0 + WU0 For i = 1 To 2 intYear = 1986 + i '用于定位 If intYear Mod 4 <> 0 And xuYear Mod 100 <> 0 Then intDays = 365 E

17、lse intDays = 366 End If '/讀入數(shù)據(jù) Open App.Path & ".Data" & intYear & ".csv" For Input As #1 ReDim t(intDays) ReDim Qobs(intDays) ReDim Eo(intDays) ReDim P1(intDays) ReDim P2(intDays) ReDim P3(intDays) ReDim P4(intDays) For j = 1 To 365 Input #1, t(j), Qobs(j), Eo(j

18、), P1(j), P2(j), P3(j), P4(j) Next j Close #1 '/參數(shù)賦值 Wum = Text2.Text: Wlm = Text3.Text: C = Text4.Text Wm = Text5.Text: B = Text6.Text: Im = Text7.Text ReDim P(intDays), EP(intDays), PE(intDays) ReDim E(intDays), EU(intDays), EL(intDays), ED(intDays) ReDim W(intDays), WU(intDays), WL(intDays),

19、WD(intDays) ReDim R(intDays) For j = 1 To intDays W(0) = W0: WU(0) = WU0: WL(0) = WL0: WD(0) = WD0 P(j) = 0.33 * P1(j) + 0.14 * P2(j) + 0.33 * P3(j) + 0.2 * P4(j) EP(j) = Kc * Eo(j) PE(j) = P(j) - EP(j) If PE(j) > 0 Then EU(j) = EP(j) '/蒸發(fā) EL(j) = 0 ED(j) = 0 Dim A As Single, Wmm As Single Wm

20、m = Wm * (1 + B) / (1 - Im) If Abs(W(j - 1) - Wm) < 0.01 Then A = Wmm Else A = Wmm * (1 - (1 - W(j - 1) / Wm) (1 / (1 + B) End If A = A + PE(j) If A <= Wmm Then R(j) = PE(j) + W(j - 1) - Wm + Wm * (1 - A / Wmm) (1 + B) Else R(j) = PE(j) + W(j - 1) - Wm End If If (WU(j - 1) + PE(j) - R(j) >

21、Wum Then '/第一層 WU(j) = Wum If (WU(j - 1) + PE(j) - R(j) - Wum) + WL(j - 1) > Wlm Then '/第二層 WL(j) = Wlm WD(j) = W(j - 1) + PE(j) - R(j) - WU(j) - WL(j) Else WL(j) = (WU(j - 1) + PE(j) - R(j) - Wum) + WL(j - 1) WD(j) = WD(j - 1) End If Else WU(j) = WU(j - 1) + PE(j) - R(j) WL(j) = WL(j - 1

22、) WD(j) = WD(j - 1) End If Else '/不產(chǎn)流，按三層蒸發(fā)算 If WU(j - 1) + P(j) >= EP(j) Then EU(j) = EP(j) EL(j) = 0 ED(j) = 0 Else EU(j) = WU(j - 1) + P(j) If WL(j - 1) > Wlm * C Then EL(j) = (EP(j) - EU(j) * (WL(j - 1) / Wlm) ED(j) = 0 Else If WL(j - 1) >= C * (EP(j) - EU(j) Then EL(j) = C * (EP(j)

23、 - EU(j) ED(j) = 0 Else EL(j) = WL(j - 1) ED(j) = C * (EP(j) - EU(j) - EL(j) End If End If End If WU(j) = WU(j - 1) + P(j) - EU(j) WL(j) = WL(j - 1) - EL(j) WD(j) = WD(j - 1) - ED(j) End If E(j) = EU(j) + EL(j) + ED(j) W(j) = WU(j) + WL(j) + WD(j) Next j W0 = W(intDays): WU0 = WU(intDays): WL0 = WL(

24、intDays): WD0 = WD(intDays) For j = 1 To intDays - 1 Robs(i) = Robs(i) + Qobs(j) * 3.6 * 24 / 553 '/年實測徑流深 Next j For j = 1 To intDays Rcal(i) = Rcal(i) + R(j) '/年計算徑流深 Next j Robs(i) = Int(Robs(i) * 10 + 0.5) / 10 Rcal(i) = Int(Rcal(i) * 10 + 0.5) / 10 jderror(i) = Int(Abs(Rcal(i) - Robs(i)

25、 * 10 + 0.5) / 10 xderror(i) = Int(Abs(Rcal(i) - Robs(i) / Robs(i) * 100 + 0.5) / 100 If xderror(i) <= 0.05 Then If intYear = 1987 Then Picture1.Print Kc; xderror(i) Else Picture2.Print Kc; xderror(i) End If End If Next i If xderror(i - 2) <= 0.05 And xderror(i - 1) <= 0.05 Then xderror(i)

26、= (xderror(i - 2) + xderror(i - 1) / 2 'xderror(i)為兩年的平均值 Picture3.Print Kc; xderror(i) End IfNext mText1.Text = InputBox("請輸入您選擇的Kc") '使用者選用自己認為合適的KcKc = Text1.TextReDim Robs(2), Rcal(2), jderror(2), xderror(3)WU0 = Text10.Text: WL0 = Text11.Text: WD0 = Text12.TextW0 = WL0 + WD0 +

27、 WU0 For i = 1 To 2 intYear = 1986 + i '/用于定位 If intYear Mod 4 <> 0 And xuYear Mod 100 <> 0 Then intDays = 365 Else intDays = 366 End If '/讀入數(shù)據(jù) Open App.Path & ".Data" & intYear & ".csv" For Input As #1 ReDim t(intDays) ReDim Qobs(intDays) ReDim Eo

28、(intDays) ReDim P1(intDays) ReDim P2(intDays) ReDim P3(intDays) ReDim P4(intDays) For j = 1 To 365 Input #1, t(j), Qobs(j), Eo(j), P1(j), P2(j), P3(j), P4(j) Next j Close #1 Wum = Text2.Text: Wlm = Text3.Text: C = Text4.Text Wm = Text5.Text: B = Text6.Text: Im = Text7.Text ReDim P(intDays), EP(intDa

29、ys), PE(intDays) ReDim E(intDays), EU(intDays), EL(intDays), ED(intDays) ReDim W(intDays), WU(intDays), WL(intDays), WD(intDays) ReDim R(intDays) For j = 1 To intDays W(0) = W0: WU(0) = WU0: WL(0) = WL0: WD(0) = WD0 P(j) = 0.33 * P1(j) + 0.14 * P2(j) + 0.33 * P3(j) + 0.2 * P4(j) EP(j) = Kc * Eo(j) P

30、E(j) = P(j) - EP(j) If PE(j) > 0 Then EU(j) = EP(j) '/蒸發(fā) EL(j) = 0 ED(j) = 0 Wmm = Wm * (1 + B) / (1 - Im) If Abs(W(j - 1) - Wm) < 0.01 Then A = Wmm Else A = Wmm * (1 - (1 - W(j - 1) / Wm) (1 / (1 + B) End If A = A + PE(j) If A <= Wmm Then R(j) = PE(j) + W(j - 1) - Wm + Wm * (1 - A / Wm

31、m) (1 + B) Else R(j) = PE(j) + W(j - 1) - Wm End If If (WU(j - 1) + PE(j) - R(j) > Wum Then '/第一層 WU(j) = Wum If (WU(j - 1) + PE(j) - R(j) - Wum) + WL(j - 1) > Wlm Then '/第二層 WL(j) = Wlm WD(j) = W(j - 1) + PE(j) - R(j) - WU(j) - WL(j) Else WL(j) = (WU(j - 1) + PE(j) - R(j) - Wum) + WL(

32、j - 1) WD(j) = WD(j - 1) End If Else WU(j) = WU(j - 1) + PE(j) - R(j) WL(j) = WL(j - 1) WD(j) = WD(j - 1) End If Else '/不產(chǎn)流，按三層蒸發(fā)算 If WU(j - 1) + P(j) >= EP(j) Then EU(j) = EP(j) EL(j) = 0 ED(j) = 0 Else EU(j) = WU(j - 1) + P(j) If WL(j - 1) > Wlm * C Then EL(j) = (EP(j) - EU(j) * (WL(j -

33、1) / Wlm) ED(j) = 0 Else If WL(j - 1) >= C * (EP(j) - EU(j) Then EL(j) = C * (EP(j) - EU(j) ED(j) = 0 Else EL(j) = WL(j - 1) ED(j) = C * (EP(j) - EU(j) - EL(j) End If End If End If WU(j) = WU(j - 1) + P(j) - EU(j) WL(j) = WL(j - 1) - EL(j) WD(j) = WD(j - 1) - ED(j) End If E(j) = EU(j) + EL(j) + E

34、D(j) W(j) = WU(j) + WL(j) + WD(j) Next j W0 = W(intDays): WU0 = WU(intDays): WL0 = WL(intDays): WD0 = WD(intDays) '/輸出文件 Open App.Path & "" & intYear & ".out" For Output As #2 For j = 1 To intDays Write #2, t(j), P(j), Eo(j), EP(j), EU(j), EL(j), ED(j), E(j), WU(j

35、), WL(j), WD(j), W(j), R(j) Next j Close #2 '/年統(tǒng)計 For j = 1 To intDays - 1 Robs(i) = Robs(i) + Qobs(j) * 3.6 * 24 / 553 '/年實測徑流深 Next j For j = 1 To intDays Rcal(i) = Rcal(i) + R(j) '/年計算徑流深 Next j Robs(i) = Format(Robs(i), "0.00") Rcal(i) = Format(Rcal(i), "0.00") jd

36、error(i) = Rcal(i) - Robs(i) xderror(i) = (Rcal(i) - Robs(i) / Robs(i) * 100 jderror(i) = Format(jderror(i), "0.00") xderror(i) = Format(xderror(i), "0.00") Print "實測徑流(" & i & ")=" & Robs(i) Print "計算徑流(" & i & ")=" &am

37、p; Rcal(i) Print "絕對誤差(" & i & ")=" & jderror(i) Print "相對誤差(" & i & ")=" & xderror(i) Next iEnd SubFORM 3Private Sub Command1_Click()'/任務二 '/讀入數(shù)據(jù) Open App.Path & ".Datapdata.csv" For Input As #1 Do While Not EOF(1

38、) intHours = intHours + 1 ReDim Preserve t(intHours) ReDim Preserve EM(intHours) ReDim Preserve P1(intHours) ReDim Preserve P2(intHours) ReDim Preserve P3(intHours) ReDim Preserve P4(intHours) Input #1, t(intHours), EM(intHours), P1(intHours), P2(intHours), P3(intHours), P4(intHours) Loop Close #1 &

39、#39;/參數(shù)賦值 Kc = Text1.Text: Wum = 20: Wlm = 60: C = 0.16 Wm = 140: B = 0.2: Im = 0.001 ReDim P(intHours), EP(intHours), PE(intHours) ReDim E(intHours), EU(intHours), EL(intHours), ED(intHours) ReDim W(intHours), WU(intHours), WL(intHours), WD(intHours) ReDim R(intHours) '/初始狀態(tài) W(0) = Val(Text2.Te

40、xt) + Val(Text3.Text) + Val(Text4.Text): WU(0) = Val(Text2.Text): WL(0) = Val(Text3.Text): WD(0) = (Text4.Text) For j = 1 To intHours P(j) = Int(0.33 * P1(j) + 0.14 * P2(j) + 0.33 * P3(j) + 0.2 * P4(j) * 10 + 0.5) / 10 EP(j) = Kc * EM(j) PE(j) = P(j) - EP(j) If PE(j) > 0 Then '/產(chǎn)流計算 EU(j) = E

41、P(j) EL(j) = 0 ED(j) = 0 Dim A As Single, Wmm As Single Wmm = (1 + B) * Wm / (1 - Im) If Abs(W(j - 1) - Wm) < 0.1 Then A = Wmm Else A = Wmm * (1 - (1 - W(j - 1) / Wm) (1 / (1 + B) End If A = A + PE(j) If A <= Wmm Then R(j) = PE(j) + W(j - 1) - Wm + Wm * (1 - A / Wmm) (1 + B) Else R(j) = PE(j)

42、+ W(j - 1) - Wm End If R(j) = Int(R(j) * 10 + 0.5) / 10 If (WU(j - 1) + PE(j) - R(j) > Wum Then '/由上至下，逐層填蓄 WU(j) = Wum If (WU(j - 1) + PE(j) - R(j) - Wum) + WL(j - 1) > Wlm Then WL(j) = Wlm WD(j) = W(j - 1) + PE(j) - R(j) - WU(j) - WL(j) Else WL(j) = (WU(j - 1) + PE(j) - R(j) - Wum) + WL(

43、j - 1) WD(j) = WD(j - 1) End If Else WU(j) = WU(j - 1) + PE(j) - R(j) WL(j) = WL(j - 1) WD(j) = WD(j - 1) End If Else If WU(j - 1) + P(j) >= EP(j) Then EU(j) = EP(j) EL(j) = 0 ED(j) = 0 Else EU(j) = WU(j - 1) + P(j) If WL(j - 1) / Wlm > C Then EL(j) = (EP(j) - EU(j) * WL(j - 1) / Wlm ED(j) = 0 Else If WL(j - 1) >= C * (EP(j) - EU(j) Then EL(j) = C * (EP(j) - EU(j) ED(j) = 0 Else EL(j) = WL(j - 1) ED(j) = C * (EP(j) - EU(j) - EL(j) End If End If End If WU(j) = WU(j - 1) + P(j) - EU(j) '/平衡 WL(j) = WL(j - 1) - EL(j) WD(j) = WD(j - 1) - ED(j) End If E(j) = EU(j) + EL(j) + ED(j) W(