计算程序

c fengfor五层准三维流（潜/弱/承/弱/承，对弱透水层只考虑越流不考虑弹性释放），调参用

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个；filee8：定义8字符的文件名

character file1（igs）8，file2（igs）8，filee8

write（，23）

23 format（1x，’zz1＝？zz2＝？迭代因子：zz1潜水0—1，zz2承压水1—2’）

c zz1：潜水为亚松弛系数，一般取085为好；zz2：承压水为超松弛系数，当ids＜300时取12-13，当ids＞300时取13-15为好

zz1＝085

zz2＝125

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：25万图面的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：25万图面的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）＝00

yo（i）＝00

25 d（i，j）＝00

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）eqiu）goto 94

if（idyh（i1，i2）ne0）goto 96

if（idyh（i1，i2）eq0）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））05

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）＝-tsss/30+e（i）

e（j）＝-tsss/30+e（j）

e（k）＝-tsss/30+e（k）

c越流矩阵元素，第ip单元1m水头差时的越流量（m3/d/m）均分给三节点i，j，k，随着ip＝ms1+1，ms的单元循环而对有关单元求其和

yo（i）＝yo（i）+rmkss/30

yo（j）＝yo（j）+rmkss/30

yo（k）＝yo（k）+rmkss/30

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/40

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（jeqidyh（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）＝00

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）＝00

yo（i）＝00

2050 d（i，j）＝00

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）eqiu）goto 9415

if（idyh（i1，i2）ne0）goto 9615

if（idyh（i1，i2）eq0）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））05

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（ngheqjbiao（jb））goto 778

777 continue

w（ngh，ikv）＝w（ngh，ikv）-qxia（ikv）ss/3

778 w（ngh，ikv）＝w（ngh，ikv）+x（ikv）/dtdykss/3+gq（ikv）dyk1ss/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）＝-tsss/30+e（i）

e（j）＝-tsss/30+e（j）

e（k）＝-tsss/30+e（k）

yo（i）＝yo（i）+rmkss/30

yo（j）＝yo（j）+rmkss/30

yo（k）＝yo（k）+rmkss/30

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/40

do 1061 k＝1，9

if（jeqidyh（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＝00

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（ieqifdh（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（keq0）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）＝reszz1+hed（i）

c 把1-n1个节点中最大的res挑出，并把其点号记到imax1中：

if（abs（res）leamax）goto 4002

amax＝abs（res）

imax1＝i

c 循环返回，继续下一个节点：

4002 continue

988 amax＝00

232 continue

write（，）’amax1＝’，amax，imax1

if（amaxgt9999999e-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＝-heddyo（i）+hedddyo（i+ids1）-e（i）（hed（i）-h0（i））/dt

50 if（ieqifdh（iqq））then

res＝res+q（iqq，ikv）

iqq＝iqq+1

endif

do 300 j＝1，9

k＝idyh（i，j）

if（keq0）goto 300

res＝res+d（i，j）hed（k）

300 continue

res＝res/（-d（i，1）+e（i）/dt）

hed（i）＝reszz2+hed（i）

if（abs（res）leamax）goto 400

amax＝abs（res）

imax2＝i

400 continue

write（，）’amax2＝’，amax，imax2

if（amaxgt9999999e-03）goto 988

c 至此，中层水全部节点末刻水头迭代完一轮，节点误差最大者如果 ＞001m，则988 句进行下一轮迭代，如果≤001m，则迭代结束，中层水本时段末刻流场计算完，以下开始深层水本时段末刻流场的计算：

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＝-heddyo（i）-e（i）（hed（i）-h0（i））/dt

5011 if（ieqifdh（iqqq））then

res＝res+q（iqqq，ikv）

iqqq＝iqqq+1

endif

do 301 j＝1，9

k＝idyh（i，j）

if（keq0）goto 301

res＝res+d（i，j）hed（k）

301 continue

res＝res/（-d（i，1）+e（i）/dt）

hed（i）＝reszz2+hed（i）

if（abs（res）leamax）goto 401

amax＝abs（res）

imax3＝i

401 continue

write（，）’amax3＝’，amax，imax3

if（amaxgt9999999e-03）goto 989

c至此，深层水全部节点本时段末刻水头迭代完一轮，节点误差最大者如果＞001m，则返回989句进行下一轮迭代，如果≤001m，则迭代结束，深层水本时段末刻流场计算完

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，2f70，2f72）

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＝’gandat’）

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，10f72）

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，f103，3f153，2x）

2030 format（1x，f103，3f153，2x，’＂’，a8，’＂’）

stop

end

1．一个C程序的执行是从（ A ）。

A) 本程序的main函数开始，到main函数结束

B) 本程序文件的第一个函数开始，到本程序文件的最后一个函数结束

C) 本程序的main函数开始，到本程序文件的最后一个函数结束

D) 本程序文件的第一个函数开始，到本程序的main函数结束

2．以下叙述正确的是（ C ）。

A) 在C程序中，main函数必须位于程序的最前面

B) C程序的每行中只能写一条语句

C) C语言本身没有输入输出语句

D) 在对一个C程序进行编译的过程中，可发现注释中的拼写错误

3．以下叙述不正确的是（ D ）。

A) 一个C源程序可由一个或多个函数组成

B) 一个C源程序必须包含一个main函数

C) C程序的基本组成单位是函数

D) 在C程序中，注释说明只能位于一条语句的后面

4．C语言规定，在一个源程序中，main函数的位置（ C ）。

A) 必须在最开始 B) 必须在系统调用的库函数的后面

C) 可以任意 D) 必须在最后

5．一个C语言源程序是由（ B ）。

A) 一个主程序和若干子程序组成 B) 函数组成

C) 若干过程组成 D) 若干子程序组成

二 填空题

1．C源程序的基本单位是 函数 。

2．一个C源程序中至少应包括一个 main函数 。

3．在一个C源程序中，注释部分两侧的分界符分别为 / 和 / 。

4．一个函数是由两部分组成的，它们分别是 函数名 和 函数体 。

5．在C语言中，一个函数的函数体一般包括 声明部分 和 执行部分 。

第三章 数据类型、运算符与表达式

一 单项选择题

1．下面四个选项中，均不是c语言关键字的选项是（ A ）。

A) define IF Type B) getc char printf

C) include scanf case D) while go pow

2．C语言中的标识符只能由字母、数字和下划线三种字符组成，且第一个字符（ C ）。

A) 必须为字母 B) 必须为下划线

C) 必须为字母或下划线 D) 可以是字母、数字和下划线中任一种字符

3．下列四个选项中，均是C语言关键字的选项是（ B ）。

A) auto enum include B) switch typedef continue

C) signed union scanf D) if struct type

4．下面四个选项中，均是不合法的用户标识符的选项是（ C ）。

A) A P_0 do B) float la0 _A

C) b-a goto int D) _123 temp int

5．以下有4组用户标识符，其中合法的一组是（ C ）。

A) For -sub Case B) 4d DO Size

C) f2_G3 IF abc D) WORD void define

二 填空题

1．若x和a均是int型变量，则计算表达式(1)后的x值为 12 ，计算表达式(2)后的x值为 4 。

(1) x=(a=4, 62)

(2) x=a=4, 62

2．若a是int型变量，则表达式(a=45, a2), a+6的值为 26 。

3．若b是int型变量，则表达式b=25/3%3的值为 2 。

4．若s是int型变量，且s=6 则下面表达式s%2+(s+1)%2的值为 1 。

5．若a,b和c均是int型变量，则计算表达式a=(b=4)+(c=2)后，a值为 6 ，b值为 4 ，c值为 2 。

6．若x和n均是int型变量，且x和n的初值均为5，则计算表达式x+=n++后x的值为 10 ，n的值为 6 。

7．若有定义：int b=7;float a=25,c=47;则表达式a+(int)(b/3(int)(a+c)/2)%4的值为 45 。

8．若有定义：int a=2,b=3;float x=35,y=25;则表达式(float)(a+b)/2+(int)x%(int)y的值为 35 。

9．若有定义：int x=3,y=2;float a=25,b=35;则表达式(x+y)%2+(int)a/(int)b的值为 1 。

10．若x和n均是整型变量，且x的初值为12，n的初值为5，则计算表达式x%=(n%=2)后x的值为 0 。

第四章 最简单的C程序设计——顺序程序设计

一 单项选择题 1～5 BDDAB

1．printf函数中用到格式符%5s，其中数字5表示输出的字符串占用5列，如果字符串长度大于5，则输出按方式（ B ）。

A) 从左起输出该字符串，右补空格

B) 按原字符长从左向右全部输出

C) 右对齐输出该字串，左补空格

D) 输出错误信息

2．已有定义int a= -2; 和输出语句 printf(“%8x”,a); 以下正确的叙述是（ D ）。

A) 整型变量的输出形式只有%d一种

B) %x是格式符的一种，它可以适用于任何一种类型的数据

C) %x是格式符的一种，其变量的值按十六进制输出，但%8x是错误的

D) %8x不是错误的格式符，其中数字8规定了输出字段的宽度

3．若x, y均定义成int型，z定义为double型，以下不合法的scanf函数调用语句是

（ D ）。

A) scanf(“%d %x, %le”, &x, &y, &z);

B) scanf(“%2d %d, %lf”, &x, &y, &z);

C) scanf(“%x %d %o”, &x, &y);

D) scanf(“%x %o%62f”, &x, &y, &z);

4．以下程序的输出结果是（ D ）。

main( )

{ int k=17；

printf("%d，%o，%x\n"，k，k，k)；

}

A)17，021，0x11 B)17，17，17

B)17，0x11，021 D)17，21，11

5．下列程序的运行结果是（ B ）。

#include <stdioh>

main()

{ int a=2,c=5;

printf("a=%d,b=%d\n",a,c);

}

A) a=%2,b=%5 B) a=2,b=5

C) a=d,b=d D) a=2,c=5

二 填空题

1．下面程序的运行结果是 i:dec=4,oct=4,hex=4,unsigned=65531 。

main()

{

short i;

i= -4;

printf(“\ni: dec=%d, oct=%o, hex=%x, unsigned=%u\n”, i, i, i, i);

}

2．若想通过以下输入语句使a=50，b=4，c=3，则输入数据的形式应该是 50，4，c=3 。

int b，c； float a；

scanf("%f，%d，c=%d"，&a，&b，&c)；

3．下列程序的输出结果是1600，请填空。

main()

{ int a=9, b=2;

float x= 66 , y=11,z;

z=a/2+bx/y+1/2;

printf("%52f\n", z );

}

4．在printf格式字符中，只能输出一个字符的格式字符是 c ；用于输出字符串的格式字符是 s ；以小数形式输出实数的格式字符是 f ；以标准指数形式输出实数的格式字符是 e 。

三 编程题

1 若a=3, b=4, c=5, x=12, y=24, z= -36, u=51274, n=128765, c1=’a’, c2=’b’。想得到以下的输出格式和结果，请写出完整的程序（包括定义变量类型和设计输出）。

要求输出的结果如下：

a= <空格>3 <空格><空格> b= <空格> 4<空格><空格> c= <空格>5

x=1200000, y=2400000, z= -3600000

x+y=<空格>360<空格><空格> y+z = -120<空格><空格>z+x= -240

u=<空格>51274<空格><空格>n=<空格><空格><空格>128756

c1= a <空格>or<空格>97(ASCII)

c2= B <空格>or<空格>98(ASCII)

/

void main()

{

int a, b, c;

unsigned u;

long int n;

float x, y, z;

char c1,c2;

a=3, b=4, c=5, x=12, y=24, z= -36, u=51274, n=128765, c1='a', c2='b';

printf("a=%2d b=%2d c=%2d\n", a, b, c);

printf("x=%f, y=%f, z=%+f\n");

printf("x+y=%52f y+z =%+52f z+x=%+52f\n", x+y, y+z, z+x);

printf("u=%6d n=%9d\n");

printf("c1=%c or %d(ASCII)\n", c1, c1);

printf("c2=%c or %d(ASCII)\n", c2-32, c2);

}

/

2 输入一个华氏温度，要求输出摄氏温度。公式为

输出要有文字说明，取2位小数。

/

void main()

{

float c, f;

printf("输入一个华氏温度：\n");

scanf("%f", &f);

c = 50/90(f-32);

printf("\n输入的华氏温度 %2f，摄氏温度为： %2f\n", f, c);

}

/

3 编程序，用getchar函数读入两个字符给c1、c2，然后分别用putchar函数和printf函数输出这两个字符。

/

void main()

{

char c1, c2;

c1 = getchar();

c2 = getchar();

putchar(c1);

putchar(c2);

printf("\nc1 = %c; c2 = %c\n", c1, c2);

}

/

第五章 选择结构程序设计

一 单项选择题

1．逻辑运算符两侧运算对象的数据类型是（ D ）。

A) 只能是0或1

B) 只能是0或非0正数

C) 只能是整型或字符型数据

D) 可以是任何类型的数据

2．已知x=43, ch='A', y=0；则表达式(x>=y&&ch<'B'&&!y)的值是（ C ）。

A) 0 B) 语法错 C) 1 D) “假”

3．已知 int x=10, y=20, z=30; 以下语句执行后x, y, z的值是（ B ）。

if(x>y)

z=x;x=y;y=z;

A) x=10,y=20,z=30 B) x=20,y=30,z=30

C) x=20,y=30,z=10 D) x=20,y=30,z=20

4．执行下列语句后a的值为（ 2 ），b的值为（ 1 ）。（ C ）

int a, b, c;

a=b=c=1;

++a|| ++b && ++c;

A) 错误 1 B) 2 2 C) 2 1 D) 1 1

5．若希望当A的值为奇数时，表达式的值为“真”，A的值为偶数时，表达式的值为“假”，则以下不能满足要求的表达式是（ C ）。

A) A%2==1 B) !(A%2==0) C) !(A%2) D) A%2

二 填空题

1．若从键盘输入58,则以下程序输出的结果是 585858 。

main()

{ int a;

scanf(“%d”,&a);

if(a>50) printf(“%d”,a);

if(a>40) printf(“%d”,a);

if(a>30) printf(“%d”,a);

}

2．下面程序的运行结果是 -1 。

main()

{

int a=2, b=3, c;

c=a;

if(a>b) c=1;

else if(a==b) c=0;

else c=-1;

printf(“%d\n”,c);

}

3．以下程序实现：输入三个整数，按从大到小的顺序进行输出。请填空。

main()

{

int x, y, z, c;

scanf(“%d %d %d”, &x, &y, &z);

if( x<y )

if( x<z )

if( y<z )

printf(“ %d %d %d”, x, y, z);

}

4．如果运行时输入字符Q，则下面程序的运行结果是 Q 。

main()

{

char ch;

scanf(“%c”, &ch);

ch=(ch>=‘A’&& ch<=‘Z’) (ch+32):ch;

ch=(ch>=‘a’&& ch<=‘z’) (ch-32):ch;

printf(“%c”, ch);

}

5．若x为int类型,请以最简单的形式写出与逻辑表达式!x等价的C语言关系表达式 x==0 。

第六章 循环控制

一 单项选择题

1．下面程序段的运行结果是（ C ）。

int n=0;

while(n++<=2);

printf(“%d”,n);

A) 2 B) 3 C) 4 D) 有语法错

2．设有程序段：

t=0;

while(printf(“”))

{t++;

if(t<3) break;

}

下面描述正确的是（ D ）。

A) 其中循环控制表达式与0等价

B) 其中循环控制表达式与'0'等价

C) 其中循环控制表达式是不合法的

D) 以上说法都不对

3．下面程序的功能是将从键盘输入的一对数，由小到大排序输出，当输入一对相等数时结束循环，请选择填空（ B ）。

#include <stdioh>

main()

{int a,b,t;

scanf(“%d%d”,&a,&b);

while( ___________ )

{if(a>b)

printf(“%d,%d\n”,a,b);

scanf("%d%d",&a,&b);

}

}

A) !a=b B) a!=b C) a==b D) a=b

4．C语言中，while和do-while循环的主要区别是（ A ）。

A) do-while的循环体至少无条件执行一次

B) while的循环控制条件比do-while的循环控制条件更严格

C) do-while允许从外部转到循环体内

D) do-while的循环体不能是复合语句

5．对以下程序段描述正确的是（ C ）。

x=-1;

do

while(!x);

A) 是死循环 B) 循环执行二次

C) 循环执行一次 D) 有语法错误

二 填空题

1．以下程序运行后的输出结果是 52 。

main()

{ int i=10, j=0;

do

while(i>2);

printf("%d\n",j);

}

2．设有以下程序:

main()

{ int n1,n2;

scanf("%d",&n2);

while(n2!=0)

{ n1=n2%10;

n2=n2/10;

printf("%d",n1);

}

}

程序运行后，如果从键盘上输入1298，则输出结果为 8921 。

3．若输入字母b，程序输出结果为 b, B ；若输入字符，程序将怎样 等待重新输入 。

#include<stdioh>

main()

{char c1, c2;

c1=getchar();

while(c1<97||c1>122)

c1=getchar();

c2=c1-32;

printf("%c, %c\n",c1, c2);

}

4．用以下程序计算1到100的整数的累加和。

main()

{ int i=1, sum= 0 ;

for( ; i<=100 ; )

{sum+=i;

i++;

}

printf("sum=%d\n", sum );

}

5．以下程序的功能是：从键盘上输入若干个学生的成绩，统计并输出最高成绩和最低成绩，当输入负数时结束。请填空。

main()

{ float x, amax, amin;

scanf(“%f”,&x);

amax=x;

amin=x;

while( x>=0 )

{if(x>amax) amax=x;

if( x<amin ) amin=x;

scanf(“%f”,&x);

}

printf("amax=%f\namin=%f\n",amax, amin);

}

方式一：

反归一化需要知道原矩阵的两个参数。设原矩阵为a，归一化的过程后得到b：

amin=min(a(:));

amax=max(a(:));

b=(a-amin)/(amax-amin);

则反归一化需要知道amax和amin是多少，用一个语句可求：

a=b(amax-amin)+amin

方式二：

使用matlab进行运算的时候，为了避免不同数量级数据之间的影响，需要把数据进行归一化，具体方法为：

1、首先打开电脑上的“matlab”软件，主界面如下图所示，在命令行输入代码即可运行。

2、matlab的归一化函数为mapminmax，此处以A=[100 200 300 400 500]为例，使用mapminmax函数进行归一化，调用格式为[A1,PS]=mapminmax(A)。A1为归一化后的数值。

3、归一化函数的对应关系为y=(ymax-ymin)(x-xmin)/(xmax-xmin)+ymin。此处以200验证某一个元素的值，那么y=2(200-100)/(500-100)+(-1)=-1/2=-05。

4、可以利用这个对应关系PS对其他数值进行归一化。