c feng.for五層準三維流（潛/弱/承/弱/承，對弱透水層只考慮越流不考慮彈性釋放），調參用
c 使潛、承、承上下節點、上下單元完全重合，且總單元數ms＝3×ms1，總節點數ids＝3×ids1
c 潛單元編號由1-ms1，中承單元編號由ms1+1-ms2，深承單元編號由ms2+1-ms
c 10-13行：可調數組的各個參數預先賦值，不同計算區只要改變這些參數即可用此程序，其中潛單元個數＝ms1個，中承單元個數＝ms2-ms1個，深承單元個數＝ms-ms2個；潛節點編號由1-ids1，中承節點編號由ids1+1ids2，深承節點編號由ids2+1-ids，n：未知水頭節點數；mi：開采節點數；igs：擬合點總數；igs1：潛擬合點數，igs2igs1：中承擬合點數，icqs：參數區個數：ihv：計算時段數，jbn：純外引地表水灌溉節點個數；nn1：潛已知水頭節點個數；nn：潛+中+深已知水頭節點總數
parameter（ms1＝108，ids1＝67，n1＝59，mi1＝17，igs1＝11，icqs1＝8）
parameter（ms2＝216，ids2＝134，n2＝126，mi2＝47，igs2＝25，icqs2＝16）
parameter（ms＝324，ids＝201，n＝193，mi＝69，igs＝35，icqs＝24）
parameter（ihv＝24，jbn＝9，nn1＝8，nn＝24）
c h0：時段初水頭（m）；hed：時段末水頭（m）；in：各單元三節點編號（必須按小中大順序排）
c fh：擬合點各時段計算水位（m）；sh：擬合點各時段實測水位（m）
dimension h0（ids），hed（ids），in（ms，3），fh（igs，ihv），sh（igs，ihv）
c idyh：導水矩陣工作單元；zb：各節點x，y坐標（輸入坐標為圖面mm數）；igdh：擬合點節點號；q：開采井各時段水量（抽為正，註為負，流出邊界為正，流入為負）（m3/d）
dimension idyh（ids，9），zb（ids，2），igdh（igs），q（mi，ihv）
c d：各節點導水矩陣；ifdh：開采井節點號；s：擬合點擬合誤差（m）；mqh：各單元的參數區號；e：各節點釋（儲）水矩陣dimension ifdh（mi），d（ids，9），s（igs，ihv），mqh（ms），e（ids）
c fqc：各參數區最多4個參數之值；sss：各節點水頭總變幅，由t0時刻起算，（m）；h：初始流場，gh潛水層獨有的灌溉入滲系數
dimension fqc（icqs，4），sss（ids），h（ids），gh（icqs1）
c yo：各節點越流矩陣工作單元；jbiao：純引地表水灌溉的節點編號，***有jbn個節點；qxia：灌溉抽取地下水均鋪在灌溉區的水層厚度；gq：灌溉抽出的地下水水層厚度；gq＝qxia；w：各時段潛水各節點的降水灌溉回歸水量
dimension yo（ids），jbiao（jbn），qxia（ihv），gq（ihv），w（ids1，ihv）
c zzz：潛水各節點底板標高（m）；ep0：潛水各節點含水層厚（m）；ep1：潛水各單元含水層厚均值（m）；x：各時段降水量（m/時段）
dimension zzz（ids1），ep0（ids1），ep1（ms1），x（ihv）
c tl：各時段步長值（天）；bc：三角單元幾何量（bi，ci，bj，cj，bk，ck）
dimension tl（ihv），bc（3，2）
c 分別是潛（h0nn1），中（h0nn2），深（h0nn3）層已知變水頭節點各時段的水頭
dimension h0nn1（nn1，ihv），h0nn2（nn1，ihv），h0nn3（nn1，ihv）
c file1（igs）*8：定義8字符的字符串igs個；filee*8：定義8字符的文件名
character file1（igs）*8，file2（igs）*8，filee*8
write（*，23）
23 format（1x，’zz1＝？zz2＝？叠代因子：zz1潛水0—1，zz2承壓水1—2’）
c zz1：潛水為亞松弛系數，壹般取0.85為好；zz2：承壓水為超松弛系數，當ids＜300時取1.2-1.3，當ids＞300時取1.3-1.5為好
zz1＝0.85
zz2＝1.25
write（*，24）ms1，ids1，n1，mi1，igs1，icqs1
write（*，24）ms2，ids2，n2，mi2，igs2，icqs2，ihv
write（*，24）ms，ids，n，mi，igs，icqs
24 format（5x，12i5）
open（1，file＝’fqc’，status＝’old’）
c 首先從’fqc’中讀取潛水的區號m10（空讀），參數1（K），參數2（μ），參數3潛-中的越流系數（k’/m’），參數4降水入滲系數（α），gh灌溉入滲系數
read（1，*）（m10，fqc（i，1），fqc（i，2），fqc（i，3），fqc（i，4），gh（i），
* i＝1，icqs1）
c 接著從’fqc’中讀取中層水的區號m10（空讀），參數1（T），參數2（μ*）
read（1，*）（m10，fqc（i，1），fqc（i，2），i＝icqs1+1，icqs2）
c 接著從’fqc’中讀取深層水的區號m10（空讀），參數1（T），參數2（μ*），參數3：中-深的越流系數（k’/m’）
read（1，*）（m10，fqc（i，1），fqc（i，2），fqc（i，3），i＝icqs2+1，icqs）
close（1）
c 參數3：中-深的越流系數（k’/m’）同樣用於中層水
do 54 i＝1，icqs1
fqc（icqs1+i，3）＝fqc（i，3）
54 continue
c qxia：灌溉抽取地下水平均鋪在灌溉區的水層厚度
open（1，file＝’gq’，status＝’old’）
read（1，*）（qxia（i），i＝1，ihv）
close（1）
c gq：灌溉抽出的地下水水層厚度
open（1，file＝’gq’，status＝’old’）
read（1，*）（gq（i），i＝1，ihv）
close（1）
c jbiao：純引地表水灌溉的節點編號，***有jbn個節點
open（1，file＝’jbiao’，status＝’old’）
read（1，*）（jbiao（i），i＝1，jbn）
close（1）
c’in’：各單元三節點編號文件（必須按小中大順序排），ia為單元號空讀
open（1，file＝’in’，status＝’old’）
read（1，*）（ia，（in（i，j），j＝1，3），i＝1，ms）
close（1）
c igdh：擬合點節點編號
open（1，file＝’igdh’，status＝’old’）
read（1，*）（igdh（i），i＝1，igs）
close（1）
c mqh：各單元的參數區號，ia為單元號空讀
open（1，file＝’mqh’，status＝’old’）
read（1，*）（ia，mqh（i），i＝1，ms）
close（1）
c ifdh：開采井節點號
open（1，file＝’ifdh’，status＝’old’）
read（1，*）（ifdh（i），i＝1，mi）
close（1）
c tl（i）：各時段值（天/時段，每月為1個時段）
open（1，file＝’tl1’，status＝’old’）
read（1，*）（tl（i），i＝1，ihv）
close（1）
c 每個時段的大氣降水量（m），x（i）/dt則換算為（m/d）
open（1，file＝’x’，status＝’old’）
read（1，*）（x（i），i＝1，ihv）
close（1）
c’zb’：各節點x，y坐標文件（此處為1：2.5萬圖面的mm數），ia為節點號空讀
open（1，file＝’zb’，status＝’old’）
read（1，*）（ia，（zb（i，j），j＝1，2），i＝1，ids）
close（1）
c 把各節點x，y坐標1：2.5萬圖面的mm數換算為實地的m數
do 888 i＝1，ids
do 888 j＝1，2
zb（i，j）＝zb（i，j）*25
888 continue
c 潛水開采井各時段的開采量（m3/d），空讀ia開采井個數，空讀ib開采井所在節點編號
open（1，file＝’qq’，status＝’old’）
read（1，*）ia，（（ib，q（i，j），i＝1，mi1），j＝1，ihv）
close（1）
c 中層水開采井第壹時段的開采量（m3/d），空讀ia開采井個數，空讀ib開采井所在節點編號
open（1，file＝’qz1’，status＝’old’）
read（1，*）ia，（ib，q（i，1），i＝mi1+1，mi2）
close（1）
c 中層水開采井第13時段的開采量（m3/d），空讀ia開采井個數，空讀ib開采井所在節點編號
open（1，file＝’qz2’，status＝’old’）
read（1，*）ia，（ib，q（i，13），i＝mi1+1，mi2）
close（1）
c 深層水開采井第壹時段的開采量（m3/d），空讀ia開采井個數，空讀ib開采井所在節點編號
open（1，file＝’qs1’，status＝’old’）
read（1，*）ia，（ib，q（i，1），i＝mi2+1，mi）
close（1）
c 深層水開采井第13時段的開采量（m3/d），空讀ia開采井個數，空讀ib開采井所在節點編號
open（1，file＝’qs2’，status＝’old’）
read（1，*）ia，（ib，q（i，13），i＝mi2+1，mi）
close（1）
c 中、深層水開采井第2-12時段的開采量等於第壹時段的開采量
do 887 j＝2，12
do 886 i＝mi1+1，mi
q（i，j）＝q（i，1）
886 continue
887 continue
c 中、深層水開采井第14-ihv時段的開采量等於第13時段的開采量
do 884 j＝14，ihv
do 883 i＝mi1+1，mi
q（i，j）＝q（i，13）
883 continue
884 continue
c sh：擬合點各時段實測水位（m），空讀ia擬合點所在節點編號
open（3，file＝’sh’，status＝’old’）
read（3，*）（ia，（sh（i，j），j＝1，ihv），i＝1，igs）
close（3）
c qc：潛水初始流場，空讀is節點編號
open（1，file＝’qc’，status＝’old’）
read（1，*）（is，h（i），i＝1，ids1）
close（1）
c zc：中層水初始流場，空讀is節點編號
open（1，file＝’zc’，status＝’old’）
read（1，*）（is，h（i），i＝ids1+1，ids2）
close（1）
c sc：深層水初始流場，空讀is節點編號
open（1，file＝’sc’，status＝’old’）
read（1，*）（is，h（i），i＝ids2+1，ids）
close（1）
c h0nn：分別是潛（h0nn1），中（h0nn2），深（h0nn3）層已知變水頭節點各時段的水頭，空讀ii節點編號
open（1，file＝’h0nn’）
read（1，*）（ii，（h0nn1（i，ikv），ikv＝1，ihv），i＝1，nn1）
read（1，*）（ii，（h0nn2（i，ikv），ikv＝1，ihv），i＝1，nn1）
read（1，*）（ii，（h0nn3（i，ikv），ikv＝1，ihv），i＝1，nn1）
close（1）
c zzz：潛水各節點底板標高（m），空讀ii節點編號
open（1，file＝’zzz’）
read（1，*）（ii，zzz（i），i＝1，ids1）
close（1）
c’file1’存放igs個字符串（擬合節點名稱占4字符，及觀測井的原編號再占4字符）
open（1，file＝’file1’，status＝’old’）
read（1，110）（file1（i），i＝1，igs）
close（1）
c’file2’存放igs個字符串（擬合節點名稱占4字符，加.dat後綴再占4字符）
open（1，file＝’file2’，status＝’old’）
read（1，110）（file2（i），i＝1，igs）
close（1）
110 format（10a8）
c 計算開始前先把全部節點的時段末刻水頭hed（i），時段初刻水頭h0（i）賦初值h（i），以使開始叠代時hed（i），h0（i）不為零
do 1993 i＝1，ids
hed（i）＝h（i）
1993 h0（i）＝h（i）
c 對中、深層承壓水導水矩陣工作單元idyh，釋水矩陣e，越流矩陣yo，導水矩陣d，賦0
do 25 i＝1+ids1，ids
do 25 j＝1，9
idyh（i，j）＝0
e（i）＝0.0
yo（i）＝0.0
25 d（i，j）＝0.0
c sum2用來累計計算區總面積（用承壓水總面積代替），先賦零
sum2＝0
c 對承壓水逐個單元計算幾何量及導水、釋水、越流矩陣
c 對承壓水第ip單元的三個節點號依次賦給i，j，k及i1，j1，k1
do 80 ip＝ms1+1，ms
i＝in（ip，1）
j＝in（ip，2）
k＝in（ip，3）
i1＝i
j1＝j
k1＝k
c idyh（i1，9）存放與i1點同單元的所有節點號，最多9個，可以用不完，即i1點的idyh可以有幾個idyh（i1，i2）＝0；當計算第ip單元時，i1點的idyh由i1占第1個（i2＝1）位置，j1，k1只能占i2＝2，3，…，9 中的位置；且先占者排前，193行：使j1，k1分兩輪到idyh中找位置；當j1找時，195行發現i1的idyh中已有j1，則跳到j1對i1的96循環體頭，換為k1對i1循環；當196行沒發現已有j1，且196行判斷此位不空時返回96循環體頭找下壹個位置，當碰到第1個空位時，由j1占據，然後跳到j1對i1循環體頭，換為k1對i1循環
do 90 iv＝1，3
idyh（i1，1）＝i1
do 94 iu＝j1，k1，k1-j1
do 96 i2＝2，9
if（idyh（i1，i2）.eq.iu）goto 94
if（idyh（i1，i2）.ne.0）goto 96
if（idyh（i1，i2）.eq.0）idyh（i1，i2）＝iu
goto 94
96 continue
94 continue
c 把ip單元的下壹個節點j提為i1，k提為j1，i降為k1，然後返回90循環頭處理ip單元的下壹個節點，循環3次則ip單元中i，j，k 3個點的idyh全部找完
iuu＝i1
i1＝j1
j1＝k1
k1＝iuu
90 continue
c 第ip單元中i，j，k三節點的x坐標賦給xi，xj，xk，y坐標賦給yi，yj，yk
xi＝zb（i，1）
xj＝zb（j，1）
xk＝zb（k，1）
yi＝zb（i，2）
yj＝zb（j，2）
yk＝zb（k，2）
c 第ip單元三角形面積ss
ss＝abs（（xj-xi）*（yk-yi）-（xk-xi）*（yj-yi））*0.5
c 累加各單元面積（這裏的sum2累加之後是中、深層兩層面積之和）
sum2＝sum2+ss
c 第ip單元的bi～bc（1，1），ci～bc（1，2），bj～b（2，1），cj～b（2，2），bk～b（3，1），ck～b（3，2）
bc（1，1）＝yj-yk
bc（1，2）＝xj-xk
bc（2，1）＝yk-yi
bc（2，2）＝xk-xi
bc（3，1）＝yi-yj
bc（3，2）＝xi-xj
c 第ip單元所在參數區號賦給jv，參數1（T）賦給txy，參數2（μ）賦給ts，參數3（k’/m’）賦給rmk
jv＝mqh（ip）
txy＝fqc（jv，1）
ts＝fqc（jv，2）
rmk＝fqc（jv，3）
c 第ip單元三節點i，j，k的釋水矩陣元素c（ii，ii）及c（ii，jj）……式中沒包括1/（2Δt）時，隨著ip＝ms1+1，ms的單元循環而對有關單元求其和
e（i）＝-ts*ss/3.0+e（i）
e（j）＝-ts*ss/3.0+e（j）
e（k）＝-ts*ss/3.0+e（k）
c越流矩陣元素，第ip單元1m水頭差時的越流量（m3/d/m）均分給三節點i，j，k，隨著ip＝ms1+1，ms的單元循環而對有關單元求其和
yo（i）＝yo（i）+rmk*ss/3.0
yo（j）＝yo（j）+rmk*ss/3.0
yo（k）＝yo（k）+rmk*ss/3.0
c 計算ip單元各節點的導水矩陣元素d（i，j），（i＝1，2，3）（j＝1，2，3）
do 100 iu＝1，3
do 104 iuu＝1，3
c 當iu＝1時，即ip單元的i節點，分別計算iuu＝1，2，3，即ip單元的i，j，k節點對i節點的ai，ai即公式**的沒求和部分
i＝in（ip，iu）
j＝in（ip，iuu）
ai＝txy*（bc（iu，1）*bc（iuu，1）+bc（iu，2）*bc（iuu，2））/ss/4.0
c 在ip單元的三個節點中，排除第3點，只讓第1，2點（i，j點）的ai加入i點的導水矩陣元素d（i，j）中，第243行的j可分別輪到i，j，k三點，但第245句的1～9個中，僅有i，j點在i點的idyh中，此句排除了第3點加入i點的導水矩陣元素d（i，j）中隨著ip＝ms1+1，ms的循環，到251句時承壓水導水矩陣已完全形成
do 106 k＝1，9
if（j.eq.idyh（i，k））d（i，k）＝ai+d（i，k）
106 continue
104 continue
100 continue
80 continue
c 至此，中、深層幾何量計算完，以下開始時段循環；時段循環中把潛水也加進來
c 在時段循環中，潛水的導水矩陣，釋水矩陣，越流矩陣與含水層厚度有關，所以在時段循環中完成，另外，承壓水與時段有關的部分也在時段循環中完成
do 280 ikv＝1，ihv
write（*，*）’ikv＝’，ikv
c 該時段步長賦給dt
dt＝tl（ikv）
c 以下開始計算潛水末刻流場
c 潛水已知水頭節點的已知水頭賦給本時段初、末刻
do 32 i＝1，nn1
hed（n1+i）＝h0nn1（i，ikv）
h0（n1+i）＝h0nn1（i，ikv）
32 continue
c 先計算潛水幾何量，因為潛水的含水層厚度M 隨時段而變，所以放在時段循環內
do 33 i＝1，ids1
h0（i）＝hed（i）
c 本時段潛水各節點的含水層厚度ep0（i）用時段初刻水位＝hed（i）減去底板標高zzz（i）求得
ep0（i）＝hed（i）-zzz（i）
c本時段潛水各節點的灌溉回歸水量加降水回歸水量w（i，ikv）先賦零，以備求和之用
w（i，ikv）＝0.0
33 continue
c 該時段潛水各單元的含水層厚度用三節點之均值
do 299 i＝1，ms1
299 ep1（i）＝（ep0（in（i，1））+ep0（in（i，2））+ep0（in（i，3）））/3.
do 2050 i＝1，ids1
c 同前，對潛水導水矩陣工作單元idyh，釋水矩陣e，越流矩陣yo，導水矩陣d，先賦0
do 2050 j＝1，9
idyh（i，j）＝0
e（i）＝0.0
yo（i）＝0.0
2050 d（i，j）＝0.0
do 8015 ip＝1，ms1
i＝in（ip，1）
j＝in（ip，2）
k＝in（ip，3）
i1＝i
j1＝j
k1＝k
do 9015 iv＝1，3
idyh（i1，1）＝i1
do 9415 iu＝j1，k1，k1-j1
do 9615 i2＝2，9
if（idyh（i1，i2）.eq.iu）goto 9415
if（idyh（i1，i2）.ne.0）goto 9615
if（idyh（i1，i2）.eq.0）idyh（i1，i2）＝iu
goto 9415
9615 continue
9415 continue
iuu＝i1
i1＝j1
j1＝k1
k1＝iuu
9015 continue
xi＝zb（i，1）
xj＝zb（j，1）
xk＝zb（k，1）
yi＝zb（i，2）
yj＝zb（j，2）
yk＝zb（k，2）
ss＝abs（（xj-xi）*（yk-yi）-（xk-xi）*（yj-yi））*0.5
jvv＝mqh（ip）
dyk＝fqc（jvv，4）
dyk1＝gh（jvv）
do 776 ig＝1，3
ngh＝in（ip，ig）
do 777 jb＝1，jbn
if（ngh.eq.jbiao（jb））goto 778
777 continue
w（ngh，ikv）＝w（ngh，ikv）-qxia（ikv）*ss/3
778 w（ngh，ikv）＝w（ngh，ikv）+x（ikv）/dt*dyk*ss/3+gq（ikv）*dyk1*ss/3
776 continue
bc（1，1）＝yj-yk
bc（1，2）＝xj-xk
bc（2，1）＝yk-yi
bc（2，2）＝xk-xi
bc（3，1）＝yi-yj
bc（3，2）＝xi-xj
jv＝mqh（ip）
txy＝fqc（jv，1）*ep1（ip）
ts＝fqc（jv，2）
rmk＝fqc（jv，3）
e（i）＝-ts*ss/3.0+e（i）
e（j）＝-ts*ss/3.0+e（j）
e（k）＝-ts*ss/3.0+e（k）
yo（i）＝yo（i）+rmk*ss/3.0
yo（j）＝yo（j）+rmk*ss/3.0
yo（k）＝yo（k）+rmk*ss/3.0
do 1001 iu＝1，3
do 1041 iuu＝1，3
i＝in（ip，iu）
j＝in（ip，iuu）
ai＝txy*（bc（iu，1）*bc（iuu，1）+bc（iu，2）*bc（iuu，2））/ss/4.0
do 1061 k＝1，9
if（j.eq.idyh（i，k））d（i，k）＝ai+d（i，k）
1061 continue
1041 continue
1001 continue
8015 continue
c 至此，該時段潛水幾何量計算完
c 時段末水頭叠代計數器ikv2，最大誤差記錄amax
ikv2＝0
9881 amax＝0.0
ikv2＝ikv2+1
write（*，34）ikv，ikv2
34 format（3x，i4，20x，’ikv2＝’，i4）
c 開采井點錄用計數器iq
iq＝1
c 對潛水n1個未知水頭節點逐點計算該時段末刻水頭
do 4002 i＝1，n1
c i節點常數項res：減該時段i點降水灌溉回歸量，減i點該時段來自承壓水的越流量（越流量計算采用時段初承壓水水頭，時段末潛水水位），減i點該時段儲存量的增加量
res＝-w（i，ikv）-（hed（i+ids1）-hed（i））*yo（i）-e（i）*（hed（i）-h0（i））/dt
if（i.eq.ifdh（iq））then
c 如果i節點恰是開采節點，則i節點常數項res中再加上該時段i點的開采量，然後開采井點錄用計數器iq加1，開采井點只能被1個節點錄用
res＝res+q（iq，ikv）
iq＝iq+1
endif
c 把與i節點***單元的k節點的導水矩陣元素（k＝1，2，…，最多可9，當k＝1時即i點自己）依次加進i節點常數項res中，k節點的導水矩陣元素采用上壹輪叠代出的hed（k）計算
do 3002 j＝1，9
k＝idyh（i，j）
if（k.eq.0）goto 3002
res＝res+d（i，j）*hed（k）
3002 continue
c 把常數項除以：（i節點對i節點自身導水矩陣元素d（i，1）的負數，加，i節點儲水陣元素e（i）/dt）
res＝res/（-d（i，1）+e（i）/dt）
c 把i節點的res乘以亞松弛系數zz1，加給上壹輪叠代出的hed（i）中，作為這壹輪叠代出的hed（i）：
hed（i）＝res*zz1+hed（i）
c 把1-n1個節點中最大的res挑出，並把其點號記到imax1中：
if（abs（res）.le.amax）goto 4002
amax＝abs（res）
imax1＝i
c 循環返回，繼續下壹個節點：
4002 continue
988 amax＝0.0
232 continue
write（*，*）’amax1＝’，amax，imax1
if（amax.gt.9.999999e-03）goto 9881
c 至此，潛水本時段末刻流場計算完，以下開始中層水本時段末刻流場的計算：
do 232 i＝1，nn1
hed（n2+i）＝h0nn2（i，ikv）
h0（n2+i）＝h0nn2（i，ikv）
iqq＝mi1+1
do 400 i＝ids1+1，n2
c 其中的越流流入量采用（潛水本時段末水位hed（i-ids1）減承壓水本時段末水位hed（i））之差hedd計算
hedd＝hed（i-ids1）-hed（i）
heddd＝hed（i）-hed（i+ids1）
res＝-hedd*yo（i）+heddd*yo（i+ids1）-e（i）*（hed（i）-h0（i））/dt
50 if（i.eq.ifdh（iqq））then
res＝res+q（iqq，ikv）
iqq＝iqq+1
endif
do 300 j＝1，9
k＝idyh（i，j）
if（k.eq.0）goto 300
res＝res+d（i，j）*hed（k）
300 continue
res＝res/（-d（i，1）+e（i）/dt）
hed（i）＝res*zz2+hed（i）
if（abs（res）.le.amax）goto 400
amax＝abs（res）
imax2＝i
400 continue
write（*，*）’amax2＝’，amax，imax2
if（amax.gt.9.999999e-03）goto 988
c 至此，中層水全部節點末刻水頭叠代完壹輪，節點誤差最大者如果 ＞0.01m，則988 句進行下壹輪叠代，如果≤0.01m，則叠代結束，中層水本時段末刻流場計算完，以下開始深層水本時段末刻流場的計算：
989 amax＝0
iqqq＝mi2+1
do 332 i＝1，nn1
hed（n+i）＝h0nn3（i，ikv）
h0（n+i）＝h0nn3（i，ikv）
332 continue
do 401 i＝ids2+1，n
hedd＝hed（i-ids1）-hed（i）
res＝-hedd*yo（i）-e（i）*（hed（i）-h0（i））/dt
5011 if（i.eq.ifdh（iqqq））then
res＝res+q（iqqq，ikv）
iqqq＝iqqq+1
endif
do 301 j＝1，9
k＝idyh（i，j）
if（k.eq.0）goto 301
res＝res+d（i，j）*hed（k）
301 continue
res＝res/（-d（i，1）+e（i）/dt）
hed（i）＝res*zz2+hed（i）
if（abs（res）.le.amax）goto 401
amax＝abs（res）
imax3＝i
401 continue
write（*，*）’amax3＝’，amax，imax3
if（amax.gt.9.999999e-03）goto 989
c至此，深層水全部節點本時段末刻水頭叠代完壹輪，節點誤差最大者如果＞0.01m，則返回989句進行下壹輪叠代，如果≤0.01m，則叠代結束，深層水本時段末刻流場計算完
c以下開始把本時段末刻水位潛，中，深層水全部節點本時段末刻水頭hed賦給下時段初刻水位h0，三層擬合點本時段末刻水位hed存入fh
189 do 281 i＝1，ids
h0（i）＝hed（i）
281 continue
do 190 j＝1，igs
mp＝igdh（j）
fh（j，ikv）＝hed（mp）
190 continue
280 continue
c 至此，三層全包括的時段循環完
c 屏幕輸出計算區總面積sum2/2，（m2）；sum2是中、深層兩層的累加和：
write（*，*）’sum＝’，sum2/2，’m2’
do 1926 i＝1，ids
1926 sss（i）＝hed（i）-h（i）
c’ok1’輸出潛水的最終流場：節點號i、坐標zb、最終流場hed、該節點0-ihv時間（最終流場比初始流場）的水位總變幅：
open（1，file＝’ok1’）
write（1，1927）（i，zb（i，1），zb（i，2），hed（i），sss（i），i＝1，ids1）
close（1）
c’ok2’輸出中層水的最終流場：節點號i、坐標zb、最終流場hed、該節點0-ihv時間（最終流場比初始流場）的水位總變幅：
open（1，file＝’ok2’）
write（1，1927）（i，zb（i，1），zb（i，2），hed（i），sss（i），i＝ids1+1，ids2）
close（1）
c’ok3’輸出深層水的最終流場：節點號i、坐標zb、最終流場hed、該節點0-ihv時間（最終流場比初始流場）的水位總變幅：
open（1，file＝’ok3’）
write（1，1927）（i，zb（i，1），zb（i，2），hed（i），sss（i），i＝ids2+1，ids）
close（1）
1927 format（i4，2f7.0，2f7.2）
c 計算出擬合點各時段誤差（計算水位-實測水位）s（i，iuv）
do 8 iuv＝1，ihv
do 2003 i＝1，igs
s（i，iuv）＝fh（i，iuv）-sh（i，iuv）
2003 continue
8 continue
c 輸出擬合點各時段誤差
open（1，file＝’gan.dat’）
do 3511 i＝1，igs，5
write（1，2993）file1（i），file1（i+1），file1（i+2），file1（i+3），file1（i+4）
2993 format（5（3x，a8，3x））
do 3778 iv＝1，ihv
f0＝fh（i，iv）
f1＝fh（i+1，iv）
f2＝fh（i+2，iv）
f3＝fh（i+3，iv）
f4＝fh（i+4，iv）
s0＝s（i，iv）
s1＝s（i+1，iv）
s2＝s（i+2，iv）
s3＝s（i+3，iv）
s4＝s（i+4，iv）
write（1，3556）iv，f0，s0，f1，s1，f2，s2，f3，s3，f4，s4
3778 continue
3556 format（i3，10f7.2）
3511 continue
close（1）
c 輸出擬合節點計算的歷時水位fh（i，iuv）、實測的歷時水位sh（i，iuv）、擬合誤差s（i，iuv），擬合節點名稱file1（i）do 2013 i＝1，igs
tt＝0
filee＝file2（i）
open（1，file＝filee）
do 2018 iuv＝1，ihv-1
tt＝tt+tl（iuv）
write（1，2020）tt，fh（i，iuv），sh（i，iuv），s（i，iuv）
2018 continue
write（1，2030）tt，fh（i，iuv），sh（i，iuv），s（i，iuv），file1（i）
close（1）
2013 continue
2020 format（1x，f10.3，3f15.3，2x）
2030 format（1x，f10.3，3f15.3，2x，’＂’，a8，’＂’）
stop
end