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
- 相关评论
- 我要评论
-