计算程序

c name：qun-w.for，无界承压含水层干扰井群法预测开采降深场，无潜水越流进入

c 计算公式：s（j）＝∑Qi/（4πT）*W（ui），ui＝Sir（i，j）**2/（4Tt），均质，各向同性s：降深（m），S：储水系数m：预测点个数n：抽水井个数ihv：预测时段个数ak：全区K值（m/d），W：全区S值

parameter（m＝4，n＝2，ihv＝60，ihv1＝36）

double precision u，ww

dimension xn（n），yn（n），r（n，m）

dimension xm（m），ym（m）

dimension q（n），s（m，ihv），it（ihv）

character file3（m）*8，fileb*8，file2（m）*4

c 抽水井点坐标（n）

open（1，file＝’xy-jing.txt’）

read（1，*）（xn（i），yn（i），i＝1，n）

close（1）

c 预测点坐标（m）

open（1，file＝’xy-yuce.txt’）

read（1，*）（xm（i），ym（i），i＝1，m）

close（1）

open（1，file＝’kmw.txt’）

read（1，*）ak，am，w

close（1）

akm＝ak*am

open（1，file＝’q.txt’）

read（1，*）（q（i），i＝1，n）

close（1）

c’yucehao.txt’存放m个字符串，每串4字符，为预测点输出预测的历时水位曲线标注字符用

open（1，file＝’yucehao.txt’）

read（1，42）（file2（i），i＝1，m）

close（1）

42 format（5a4）

c’yucewen.txt’存放m个字符串，为’预测号.dat’，每串10字符，为预测点输出预测的历时水位曲线准备文件名

open（1，file＝’yucewen.txt’）

read（1，44）（file3（i），i＝1，m）

close（1）

44 format（5a8）

open（1，file＝’tl1.txt’）

read（1，*）（it（i），i＝1，ihv）

close（1）

it0＝0

c 各预测点j至抽水井点i之距r（m）

do 130 j＝1，m

s（j，0）＝0

do 120 i＝1，n

r（i，j）＝sqrt（（xn（i）-xm（j））**2+（yn（i）-ym（j））**2）

if（r（i，j）.eq.0）r（i，j）＝0.15

120 continue

130 continue

do 160 k＝1，ihv1

write（*，*）’k＝’，k，’t＝’，it（k）

c 对m个预测点逐点计算1～n眼抽水井在k时刻对每个预测点水位降之和s（j，k）

do 150 j＝1，m

s（j，k）＝0

c 计算k时刻i～j泰斯井函数自变量u

do 140 i＝1，n

u＝r（i，j）**2*w/（4*akm*it（k））

c 调用子程序计算k时刻i～j泰斯井函数值ww

call wu（u，ww）

c 用泰斯公式计算k时刻第i眼抽水井对第j个预测点的水位降ss

ss＝q（i）/（4*3.1415926*akm）*ww

if（ss.lt.0）ss＝0

c i＝1～n眼抽水井在k时刻对第j预测点水位降求和s（j，k）

s（j，k）＝s（j，k）+ss

140 continue

150 continue

160 continue

open（1，file＝’sids.txt’）

do 190 k＝1，ihv1

write（1，170）k

170 format（/2x，’t＝’，i4，’月’）

do 200 j＝1，m

200 write（1，180）xm（j），ym（j），-s（j，k），j

180 format（2f10.1，f8.2，i6）

190 continue

close（1）

do 300 j＝1，m

fileb＝file3（j）

open（1，file＝fileb）

write（1，310）it0，s（j，0），file2（j）

do 320 k＝1，ihv1

write（1，380）it（k），-s（j，k）

320 continue

300 close（1）

380 format（i10，f8.2）

310 format（i10，f8.2，1x，a4）

cq＝0

do 90 i＝1，n

90 cq＝cq+q（i）

open（1，file＝’qq.txt’）

write（1，142）（i，q（i），i＝1，n）

write（1，*）

write（1，*）’总采量＝’，cq，’（m3/d）’

close（1）

142 format（5（i4，f8.0））

open（1，file＝’okk.dat’）

do 290 j＝1，m

290 write（1，80）xm（j），ym（j），-s（j，ihv1），j

close（1）

80 format（2f10.1，f9.3，i4）

stop

end

c 泰斯井函数检验程序，与泰斯井函数表对照，计算完全正确

subroutine wu（u，w）

double precision u，w

dimension ni（100）

c 泰斯井函值W计算

w＝0

ni（1）＝1

do 30m＝2，10

30 ni（m）＝ni（m-1）*m

do 60 n＝1，10

60 w＝w+（-1）**（n+1）*u**n/（n*ni（n））

w＝w-0.577216-dlog（u）

10 continue

if（w.lt.0）w＝0

return

end

用time()函数，不要用clock()！

#include <time.h>

time_t begin = time(0)

time_t end = time(0)

time_t cost = end - begin

printf("%d seconds\n", cost)

就可以了

---